一种磁性源瞬变电磁法时域矢量有限元正演方法及装置

1.本发明涉及地球物理技术领域，更具体涉及一种磁性源瞬变电磁法时域矢量有限元正演方法及装置。


背景技术：

2.电磁学的发展体现在包括地球物理领域在内的各种应用学科当中，由于历史原因，电磁场理论的重点更多的是放在频域电磁学中，特别是求解时谐稳态电磁场的理论，目前已经相当成熟。随着电子计算机的诞生和发展，时域电磁场分析的数值方法逐渐占据主要地位，成为地球物理电磁法领域的研究热点。
3.在地球物理领域，电磁法是不可或缺的一类方法，其中包括瞬变电磁法，也称时间域电磁法。由于理论的复杂性和仪器的制约，瞬变电磁法目前三维正演过程以有限差分方法和频率域有限元方法为主。求解瞬变电磁场主要有两种思路：一种是间接求解，先将时间域电磁场的计算式变换到变换域中，在变换域中求得电磁响应，然后再转换到时间域；另一种是直接在时域求解电磁场，如常用的是时间域的有限差分方法。
4.目前对于磁性源瞬变电磁的主要的研究工作在第一种思路。然而第一种思路在时频转换的过程是存在误差的，根据傅里叶变化的定义，只有足够多的频点才能进行反变换到时间域以保证精度，如此一来，效率问题和精度问题会受到所用的数值积分方法的影响，其次，这种方法对源的要求过于严格，只能解决几种理想的激励源，然而瞬变电磁在实际应用中，激励源的激发电流波形往往不是理想的梯形或斜阶跃波形，间接法显然无法准确的模拟任意发射波形的数值模拟。
5.在第二种思路中，时域数值方法有很多频域方法所不具备的优点，首先时域方法是直接模拟物理过程的方法，不仅具有更形象的模拟过程，并且具有实际的物理意义。其次，在一些高频、瞬态场的模拟中，时域相较于频域更具优势。尽管如此，瞬变电磁法的时域有限元三维模拟仍存在阻碍，其中最重要的一点是效率问题，瞬变电磁法需要模拟快速关断电流的扩散场，这要求时间步长是较小的，因此传统的有限元方法需要进行上万步迭代，并且在时间步长改变时，需要求解上百、上千万阶的大型稀疏方程组，这极大的影响了导致瞬变电磁法的三维正演效率较低。


技术实现要素：

6.本发明所要解决的技术问题在于提供了一种磁性源瞬变电磁法时域矢量有限元正演方法及装置，以提高瞬变电磁法的三维正演效率。
7.本发明是通过以下技术方案解决上述技术问题的：
8.本发明提供了一种磁性源瞬变电磁法时域矢量有限元正演方法，所述方法包括：
9.建立模拟地质模型，对模拟地质模型进行网格剖分，得到网格剖分属性信息，其中，模拟地质模型为待求解的地质区域；且网格剖分属性信息包括：节点编号、棱边编号、单元编号以及节点坐标信息；
10.对网格剖分后的模拟地质模型施加激励源，获取模拟地质模型的麦克斯韦方程组，从时间域麦克斯韦方程组中得到电场控制方程；
11.采用加权余量法离散控制方程并选择第一类whitney基函数建立有限元方程；
12.选择时间差分格式对有限元方程进行时间离散；
13.组装总体刚度矩阵，并调用padiso求解器求解大型线性稀疏方程组；
14.根据激励源的电流关断特征确定迭代的时间步长，再以可变的时间步长进行迭代得到电场各个分量，并通过基函数得到各个节点电场值；根据法拉第电磁感应定律，求得瞬变电磁观测值和视电阻率参数。
15.可选的，所述采用加权余量法离散控制方程，包括：
16.利用公式，采用加权余量法离散控制方程，其中，
17.为微分算符；v为积分域；e为时间域电场强度；nj为第一类whitney基函数；t为时间；μ为磁导率：σ为计算域中的介质的电导率；j为发射电流的电流密度；下角标j为权系数的序号。
18.可选的，所述选择时间差分格式对有限元方程进行时间离散，包括：
19.利用公式，采用无条件稳定的隐式欧拉法离散有限元方程，其中，
20.n为第n次时间步的迭代。
21.可选的，在迭代过程中，将电场值附在剖分后得到的网格单元棱边上。
22.可选的，在迭代过程中，在电流关断之后，保持迭代的时间步长不变。
23.可选的，在迭代过程中，在激励源的激发电流上升期，时间步长保持10-7s；
24.在激发电流持续期，先按照指数增大，随后保持不变；在激发电流下降期，步长保持10-7s不变，；在电流关断之后，按照对数域等间隔增大，并每次增大一次保持一段时间步长不变。
25.本发明还提供了一种磁性源瞬变电磁法时域矢量有限元正演装置，所述装置包括：
26.剖分模块，用于建立模拟地质模型，对模拟地质模型进行网格剖分，得到网格剖分属性信息，其中，网格剖分属性信息包括：节点编号、棱边编号、单元编号以及节点坐标信息；
27.激励模块，用于对网格剖分后的模拟地质模型施加激励源，获取模拟地质模型的麦克斯韦方程组，从时间域麦克斯韦方程组中得到电场控制方程；
28.离散模块，用于采用加权余量法离散控制方程并选择第一类whitney基函数建立有限元方程；
29.选择时间差分格式对有限元方程进行时间离散；
30.求解模块，用于组装总体刚度矩阵，并调用padiso求解器求解大型线性稀疏方程组；
31.迭代模块，用于根据激励源的电流关断特征确定迭代的时间步长，再以可变的时间步长进行迭代得到电场各个分量，并通过基函数得到各个节点电场值；根据法拉第电磁感应定律，求得瞬变电磁观测值和视电阻率参数。
32.本发明相比现有技术具有以下优点：
33.本发明为一种时域矢量有限元瞬变电磁正演方法。本发明首先从时间域麦克斯韦方程组中得到电场控制方程，并采用加权余量法离散控制方程得到有限元方程；其次，采用结构化六面体单元剖分计算域，选择第一类whitney基函数进行单元分析；随后组装总体刚度矩阵，并调用padiso求解器求解大型线性稀疏方程组，最后按照以可变的时间步长进行迭代，得到电场各个分量，并通过基函数得到各个节点电场值，相对于现有技术中的固定迭代步长，极大的减少了迭代次数，即能保持较高的计算精度，同时节省了计算时间，而且计算机模拟结果也证明了这一点。
34.另外，本发明具有良好的通用性，可用于航空、浅海、地面以及巷道中任何形状的磁性源瞬变电磁法正演数值模拟
附图说明
35.图1为本发明实施例提供的一种磁性源瞬变电磁法时域矢量有限元正演方法的流程示意图；
36.图2为本发明实施例提供的一种磁性源瞬变电磁法时域矢量有限元正演方法的原理示意图；
37.图3为本发明实施例1采用的网格剖分策略示意图；
38.图4为本发明实施例1采用的第一类whitney基函数示意图；
39.图5为实施例1中a型地层模型的数值模拟结果图；
40.图6为实施例1中均匀半空间中的低阻异常体模型x-z剖面和x-y平面示意图；
41.图7为实施例1中均匀半空间中的低阻体模型，在关断之后不同时刻的波场快照。
具体实施方式
42.下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
43.本发明实施例优选适用于利用瞬变电磁法进行地面、大气、浅海等场景下的勘测。本发明所要解决的问题在于提供一种自适应步长的磁性源瞬变电磁正演的方法。特点在于能够在较少的时间步之内，实现对任意发射源、任何规则异常体的地面、航空、浅海瞬变电磁勘探方法的数值模拟。
44.实施例1
45.图1为本发明实施例提供的一种磁性源瞬变电磁法时域矢量有限元正演方法的流程示意图；图2为本发明实施例提供的一种磁性源瞬变电磁法时域矢量有限元正演方法的
原理示意图，如图1和图2所示，方法包括：
46.s101：建立模拟地质模型，对模拟地质模型进行网格剖分，得到网格剖分属性信息，其中，网格剖分属性信息包括：节点编号、棱边编号、单元编号以及节点坐标信息。
47.在实际应用中，模拟地质模型是将目标地域的实际地质环境映射到数字空间后得到的数字模型。对模拟地质模型进行正交剖分后得到若干个网格，每一个网格均为六面体单元，这种剖分方式网格剖分容易，单元、节点编号及棱边编号规律性强，易于编程实现。为了对网格剖分方式进行详细说明，图3为本发明实施例1采用的网格剖分策略示意图，如图3所示，图3中整个立方体即为地质模型的尺寸，地质模型为边长10000厘米的立方体。相应的，激励源沿立方体6个面的中线设置，每一个激励源为边长为10000厘米的方形口字形线路，合计设置3个激励源。在靠近激励源的位置采用最小网格步长为5m，随着远离激励源，每增加一个网格的距离，下一个网格宽度变为上一个网格宽度的1.5倍，例如，与激励源临近的第一个网格的步长为5米，则在远离激励源方向上，第一个网格的下一个网格的步长为7.5米，第三个网格的步长为7.5*1.5＝11.25米，依次类推。异常体以及观测点位置附近采用同样的剖分策略，在异常体所在的区域采用最小为异常体几何尺寸1/20的网格步长，并按照1.5倍增大，在远离异常体、激发源以及测量点时网格逐渐增大，适当提高网格增大倍数，最大满足相邻网格步长不超过两倍关系。
48.s102：对网格剖分后的模拟地质模型施加激励源，获取模拟地质模型的麦克斯韦方程组，从时间域麦克斯韦方程组中得到电场控制方程。
49.具体的，可以忽略高频波动项，如位移电流；并对模拟地质模型施加激励源与边界条件。由于宏观世界中的电磁问题几乎都服从于麦克斯韦方程组，因此在在各向同性、且均匀的介质中，对时间域麦克斯韦方程组的旋度方程(似稳态)：对公式两边取旋度；
50.再联合公式，就可以从麦克斯韦方程组得到电磁场扩散方程。可以理解的是，基于时间域麦克斯韦方程组得到电磁场扩散方程的数学过程为现有技术，本发明实施例在此并不对其进行赘述。
51.然后，再将电磁场扩散方程作为电场的控制方程，得到的控制方程如下：
[0052][0053]
其中，为微分算符；r为距离；t为时间；μ为磁导率：e(r,t)为时间域的电场磁场描述，j(r,t)为发射电流的电流密度，σ为计算域中的介质的电导率；ε为空间中的介电常数；μ为空间中的磁导率。
[0054]
时间域控制方程是包含位置与时间的四维方程，因此为了求解该方程需要在空间和时间维度上均采取离散来满足电子计算机的计算条件，这一点与频率域只考虑空间剖分的方式有所不同。
[0055]
进一步的，为了使本步骤中的控制方程具有唯一解，需对时间域的电场磁场描述施加dirichlet边界条件：
[0056]n×
e(r,t)＝0
[0057]
其中，n为边界的外法向单位向量。
[0058]
s103：图4为本发明实施例1采用的第一类whitney基函数示意图，如图4所示，采用加权余量法离散控制方程并选择第一类whitney基函数建立有限元方程。
[0059]
例如，可以定义表述待解变值问题的余量为：r＝a(φ')-p，其中，r为近似解与真实解的余量，a为微分算符，p为解析解的响应。
[0060]
则所定义的r满足以下积分：
[0061][0062]
其中，ωj为一组权系数；ω为求解域；a为微分算符；αi为基函数的加权系数；下角标j为权系数的序号；下角标i为基函数的序号；fi为一组基函数。
[0063]
令fi＝ωj，则s102所述控制方程的加权余量应当满足：
[0064]
其中，
[0065]
v为积分域。
[0066]
因此，考虑到边界条件，可以将s102步骤中的控制方程的离散为：
[0067][0068]
其中，nj为第一类whitney基函数(相当于上式fi)，v为积分域；e为时间域电场强度。
[0069]
s104：在时间维度上，选择时间差分格式对有限元方程进行时间离散，这是时间域有限元法的关键步骤，通过不同时刻的电流幅值可以模拟任意的发射波形以及模拟物理场随时间的变化。
[0070]
为方便起见，先将s103步骤中得到的离散后的控制方程改写为：
[0071]
其中，
[0072]
再采用无条件稳定的隐式欧拉法离散s103步骤中得到的有限元方程：
[0073][0074]
其中，n为第n次时间步的迭代。
[0075]
s105：组装总体刚度矩阵，并调用padiso求解器求解大型线性稀疏方程组。具体的，在某个时间步内，对每个单元网格依次进行单元分析，读取节点编号、棱边编号、单元编
号以及节点坐标计算单元刚度矩阵。
[0076]
首先，将步骤s104中离散后的限元方程可改写为：ke·
{ee}＝be，其中，
[0077][0078]
且ke为系数矩阵，{ee}为待求电磁场量，即单元棱边组成的列向量，be为右端已知项，其中包含着场源信息，且改写后的公式中，
[0079][0080][0081]
其中，其中，为对角阵的x分量的子矩阵；[m1]、[m2]、[m3]是由六面体单元决定的系数矩阵；为对角阵的y分量的子矩阵；为对角阵的z分量的子矩阵；和为分块矩阵的子矩阵；l
x
,ly,lz分别为剖分后得到的网格单元的棱边长度。
[0082]
则单元内任意一点的电场值可由与其相关的4条平行棱边上的电场进行表示：
[0083]
其中，是电磁场量分别在x，y和z方向上的电场分量。
[0084]
使用表示单元的几何中心坐标，使用l
x
、ly和lz分别表示网格单元单元在不同方向上的棱边。则上式中的基函数ne可以改写为：
[0085][0086]
其中，
[0087]
分别为单元内的4条x方向棱边上的基函数；分别为单元内的4条x方向棱边上的基函数；为单元的几何中心坐标；x、y、z为计算的场点坐标；为为单元内的4条y方向棱边上的基函数；为分别为单元内的4条z方向棱边上的基函数。
[0088]
在当前次迭代过程中，组装所有单元的刚度矩阵组成总体刚度矩阵，得到大型线性稀疏方程组，采用求解器padiso直接求解大型线性方程组，得到当前次迭代在在当前时刻的空间所有计算点的电场、磁场的三个分量。在求解过程中，将电场值附在单元棱边上，而非顶点上，对于s101所选结构化六面体单元，节点有限元一共具有24个未知数，而棱边基函数只有12个待求未知量，计算成本远小于节点有限元方法。其次，节点有限元方法中顶点是相邻几个单元的公用顶点，但在顶点处，两个相邻单元取值却不一定相同，需要做散度矫正。本发明实施例中的矢量有限元将场值定义在棱边上，自然满足电磁场的连续性条件，因而，相对于现有技术不需要散度矫正，提高了计算效率。
[0089]
s105：根据激励源的电流关断特征确定迭代的时间步长，再以可变的时间步长进行迭代得到电场各个分量，并通过基函数得到各个节点电场值；根据法拉第电磁感应定律，求得瞬变电磁观测值和视电阻率参数。
[0090]
当前次迭代对应的时间初始值可以从0时刻开始，此时计算域中的激励源还没有开始工作，电磁场值均为0，迭代的时间步长为：
[0091][0092]
其中，λ为计算频段的最高频率的波长，c为光速。
[0093]
然后，将0
△
t时刻作为当前次迭代的下一次迭代的起始时刻，然后将当前次迭代的下一次迭代作为当前次迭代，并返回执行s105的步骤，直至收敛，或者迭代次数达到设定次数。
[0094]
在本步骤中，为了进一步减少迭代步数，提高正演效率，本方法采用两种方式处理时间步长：
[0095]
其一，由于控制方程右端项仅与源和时间步长有关，因此在电流关断之后，激励源消失，右端项变为仅与时间步长相关，若保持时间步长不变，右端项系数便是相同的矩阵，这样就可以在每次时间步长改变之后，分解一次矩阵便可以重复回代，极大的提高了瞬变电磁法的三维正演迭代过程的计算效率。
[0096]
其二，需要适时的增大时间步长。具体操作如下：在激励源的激发电流上升期，时间步长保持10-7
s，在激发电流持续期，先按照指数形式增大，比如前一个时间步长是10-7
s，下一个时间步长就是10-6
。随后在激发电流持续期内迭代步长保持不变，在激发电流下降期，步长保持10-7
s不变。在电流关断之后，首先在关断早期保持较小的时间步长
△
t，并在迭代n次之后记录场值，然后再计算一次(
△
t
×
n)步长的场值，按照物理规律此时相同的时刻物理场应相同，考虑到数值计算的误差，给出一个阈值σ，如果两次计算的场值满足阈值，则认为增大之后的时间步稳定，此时选择增大一次时间步长，然后再次进行这样的迭代过程，直至下一次时间步长的后验误差满足精度要求即可再次增大步长。
[0097]
再根据基函数计算得到单元节点上的电场值，并根据法拉第电磁感应定律，求得瞬变电磁观测值和视电阻率参数。
[0098]
为了验证本发明中方法的正确性，设计一个a型地层模型，模型参数如下所示：ρ1＝50ω
·
m、ρ2＝100ω
·
m、ρ3＝200ω
·
m、h1＝100ω
·
m、h2＝200ω
·
m。观测点选择原点位置。图5为实施例1中a型地层模型的数值模拟结果图，从图5双对数坐标系可以看出，图5中的(a)图反应的是计算的感应电动势的衰减曲线，二者衰减曲线一致。图5中的(b)图反应的是本发明实施例的计算结果与理论解进行对比的误差曲线，本发明所述方法的正演模拟结果与理论解误差在绝大多数时间内，误差均在1％之下，而仅仅在极早期误差较大，验证了本发明所述方法的有效性。
[0099]
为进一步验证本文所述方法的效果，设计一个低阻异常体模型，图6为实施例1中均匀半空间中的低阻异常体模型x-z剖面和x-y平面示意图，如图6中的(a)图所示，均匀半空间背景电阻率为100ω
·
m，建立一个50m
×
50m
×
50m的低阻立方体模型，模型电阻率为500ω
·
m，其在三维空间中的位置如图6所示。异常体的顶面埋深是50米。计算域为80000m
×
80000m
×
60000m，采用非等间距的网格，在异常体处剖分间距略小，最小间距为5m。如图6中的(b)图所示，发射源为边长为100m的矩形线框，线框中心位于计算域原点位置；发射波形仍采用梯形全波形，上升沿与下降沿为10-6
s，激励源的激发的激发电流持续期为4ms。
[0100]
图7为实施例1中均匀半空间中的低阻体模型，在关断之后不同时刻的波场快照，如图7所示，从图7中可以看出，磁性激发源激发的电流环在未扩散至高阻异常体前是对称的，而扩散至异常体时，电场传播的更快，之后，有高阻异常体的一侧，电磁波传播速度更
快，左侧涡流环中心相较于右侧向下传播速度更快，而到了晚期，异常体所在位置显示出一个场强较强的区域，这是因为，在高阻体内部传播速度较快，而外部速度较慢，因此，使得异常体的位置形成了一个“突出”的场强较大的区域。此结果与物理规律相符，证明了本发明所述方法的正确性。
[0101]
本发明所述方法提供了一种自适应步长的磁性源瞬变电磁三维正演的方法，该方法能够提高瞬变电磁时间域三维正演效率，并适用于任何发射电流波形，进而尽可能的保留了时域信息。本发明为快速、高效的三维时间域反演提供了正演算子，进而可以更好的用于航空、浅海、地面和巷道瞬变电磁问题中。本发明的另一个有益效果是，本发明所述方法，采用对数等间隔增大的时间步长，并在每次增大之后保持时间步长不变，极大的减少了迭代次数，即能保持较高的计算精度，同时节省了计算时间，提高了瞬变电磁法的三维正演过程的效率。
[0102]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。