一种基于正则化修正方程的三维大地电磁数值模拟方法

1.本发明涉及地球物理技术领域，具体涉及一种三维大地电磁正演数值模拟方法。


背景技术：

2.大地电磁测深法是利用天然大地电磁场作为场源探测地下介质的电性结构进行深部地质构造研究的一种频率域电磁测深方法。其应用离不开大地电磁数据解释，特别是高效、精确的正演算法。但是由于电磁双旋度扩散方程具有丰富的零空间，且当频率较低时无法模拟电荷积累的情况，因此，采用传统的迭代方法求解大地电磁正演问题时收敛较慢。为了提高收敛效率，通常采用预处理器对方程系数矩阵的条件数进行改善，并采用散度校正来保证电流的连续性，但是效果并不明显。而将散度校正方程整合到天然电磁场的双旋度方程中时，方程会退化为拉普拉斯方程，从而容易求解。
3.综上所述，急需一种三维大地电磁数值模拟方法，以解决现有技术中存在的问题。


技术实现要素：

4.本发明目的在于提供一种基于正则化修正方程的三维大地电磁数值模拟方法，能够有效解决在异常体电导率不同的情况下单元体之间电流不连续、散度不为零的问题，具体技术方案如下：
5.一种基于正则化修正方程的三维大地电磁数值模拟方法，具体包括如下步骤：
6.根据勘探研究目标体的几何特征和电阻率特征圈定研究区域；在坐标系条件下将所述研究区域剖分成一系列单元体，构建研究所用的模型；基于麦克斯韦方程组以及目标体的几何特征和电阻率特征，构建天然电磁场的控制方程；
7.用模型的电导率σ乘以电场强度e，再连续求散度和梯度，并带入天然电磁场的控制方程，则获得新的电磁场控制方程：其中：是旋度，是对电场强度e求双旋度，λ是修正因子，ω是角频率，μ是磁导率；
8.对模型的内部棱边ii的电场值进行修正，则得到修正后的控制方程为：a
neweii
＝b；其中：a
new
是新的电磁场控制方程的系数矩阵，b为待求方程的右端项，e
ii
是模型的内部棱边ii上的待求电场；上式中：a
new
＝a
ii-gd，b＝-a
ibeib
；其中：a
ii
是模型的内部棱边ii对应的系数矩阵，a
ib
是模型的外部棱边ib对应的系数矩阵，e
ib
是模型的外部棱边ib上的待求电场，矩阵gd是修正项的系数矩阵；
9.对矩阵gd采用加权余量法得到加权余量方程，然后用节点有限元法进行离散；
10.对修正后的控制方程进行迭代求解，分别获得xy-极化模式下和yx-极化模式下各个测点的电场分量，并求解各个测点相应的磁场分量；
11.根据各个测点的电场分量和磁场分量计算各个测点的视电阻率和相位。
12.优选的，所述修正因子λ为模型的电导率σ的倒数，即
13.优选的，用矩阵形式表示天然电磁场的控制方程为：ae＝0；
14.其中：矩阵a是天然电磁场的控制方程的系数矩阵；e为模型所有棱边的电场值，所述棱边包括内部棱边ii和外部棱边ib。
15.优选的，所述天然电磁场的控制方程的矩阵形式表示为：
16.a
iieii
+a
ibeib
＝0；
17.其中：a
iieii
＝-a
ibeib
，右端项b＝-a
ibeib
。
18.优选的，电流密度的散度为零的条件为：
[0019][0020]
其中：j是电流密度且j＝σe，φ是静电势。
[0021]
优选的，采用节点有限元法离散的方程为：
[0022][0023]
其中：ne是单元总数，是未知量，μh是节点有限元的插值基函数，dv是体积分。
[0024]
优选的，基于各个测点的电场分量根据如下公式求得各个测点的磁场分量h：
[0025][0026]
优选的，各个测点的视电阻率和相位的计算包括如下步骤：
[0027]
s1、根据如下公式计算在xy-极化模式下和yx-极化模式下的阻抗：
[0028]
以及
[0029]
其中：z
xy
是在xy-极化模式下的阻抗；z
yx
是在yx-极化模式下的阻抗；
[0030]
s2、根据如下公式计算在xy-极化模式下的视电阻率和相位：
[0031]
以及
[0032]
其中：ρ
xy
是在xy-极化模式下的视电阻率；是在xy-极化模式下的相位；
[0033]
根据如下公式计算在yx-极化模式下的视电阻率和相位：
[0034]
以及
[0035]
其中：ρ
yx
是在yx-极化模式下的视电阻率；是在yx-极化模式下的相位。
[0036]
优选的，所述单元体为规则的六面体单元。
[0037]
优选的，所述单元体的电阻率恒定。
[0038]
应用本发明的技术方案，具有的有益效果为：采用本发明的三维大地电磁数值模拟方法，能够有效解决在异常体电导率不同的情况下单元体接触表面电流不连续的问题，对单元体的静态势进行计算，并且从近似电场中减去静态势的梯度的散度，就能符合电流密度的散度为零的条件；由于矢量有限元法中单元内部符合电流散度为零的条件，因此正则化修正只需要考虑单元的面积分，简化了积分的难度，能够满足电磁数据精细、反演成像快速的需求。
[0039]
除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。
附图说明
[0040]
构成本技术的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
[0041]
图1为本发明优选实施例提供的基于正则化修正方程的三维大地电磁数值模拟方法的流程图；
[0042]
图2为本发明优选实施例提供的含有异常体的模型示意图；
[0043]
图3为采用本发明优选实施例提供的三维大地电磁数值模拟方法和传统有限元法在xy-极化模式下计算的视电阻率拟合图；
[0044]
图4为采用本发明优选实施例提供的三维大地电磁数值模拟方法和传统有限元法在yx-极化模式下计算的视电阻率拟合图；
[0045]
图5为采用本发明优选实施例提供的三维大地电磁数值模拟方法和传统有限元法在xy-极化模式下计算的相位拟合图；
[0046]
图6为为采用本发明优选实施例提供的三维大地电磁数值模拟方法和传统有限元法在 yx-极化模式下计算的相位拟合图；
[0047]
图7为采用现有技术中未进行正则化修正的方法(即传统有限元方法)以及本实施例中方法进行迭代次数的对比图；
[0048]
图8为采用现有技术中未进行正则化修正的方法(即传统有限元方法)以及本实施例中方法进行迭代计算时间对比图。
具体实施方式
[0049]
以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。
[0050]
实施例：
[0051]
参见图1，本实施例公开一种优选的基于正则化修正方程的三维大地电磁数值模拟方法的流程图，具体包括：
[0052]
根据勘探目标体的几何特征和电阻率特征，圈定研究范围，构建相应模型：
[0053]
以三维大地异常体正上方的地表作为坐标原点，建立直角坐标系。根据勘探研究目标体的几何特征和电阻率特征圈定研究区域，并在坐标系条件下将所述研究区域剖分成一系列单元体，构建研究用的模型。将模型沿着x方向、y方向以及z方向离散为m个六面体的小单元，离散的具体步骤如下：
[0054]
对于x方向和y方向，距离目标体较近的区域剖分较密，随着与目标体之间的距离增加，网格越来越疏，直至无穷远处；
[0055]
对于z方向，地表采用较密的网格，随着深度的增加，网格越来越疏，以此在保证精度的前提下有效的减少单元数目，从而进行快速正演；
[0056]
对于空气层，采用指数方式进行剖分，空气层的厚度一般设置为106m；
[0057]
对于单元体，每个单元体的几何特征用dx、dy以及dz表示，模型在x方向、y方向以
及z方向上的单元数量分比为n
x
、ny以及nz，则m＝n
x
×
ny×
nz。
[0058]
对上述m个六面体小单元以及其棱边、节点以及表面进行统一编号，根据单元体棱边的位置关系，将棱边分为内部棱边ii和外部棱边ib。根据拟定的电阻率分布，对所有的小单元的电阻率进行赋值，异常体由若干个小单元的电阻率值来体现。
[0059]
根据麦克斯韦方程组构建基于电场的双旋度控制方程：
[0060][0061]
其中：是旋度，是对电场强度e求双旋度，待测量e是电场强度，λ是修正因子，ω是角频率，μ是磁导率数值。
[0062]
对控制方程进行矢量有限元离散并给定边界条件：
[0063]
采用galerkin加权余量法得到加权余量方程后，用矢量有限元法求解：
[0064][0065]
其中：ω是研究区域；δe为点乘上式，是一个整体，物理上是电场强度e的变分。
[0066]
相应的，根据矢量恒等式及散度定理有：
[0067][0068]
其中：k是传播系数，efgh表示研究区域的下底面，且研究区域的下底面是均匀介质；假设电导率为σ0，则e
x
是电场在x方向上的分量。
[0069]
相应的，采用矢量有限元法有：
[0070][0071]
其中：k
1e
是12
×
12的矩阵，其由九个小的矩阵组成(k
11-k
33
)，每个小的矩阵都是维度为4
×
4的方阵，e是一个单元，插值基函数不同会导致k
1e
的值具有差异；是δe 在一个单元上的分量，是e在一个单元上的分量；上标t是转置。
[0072]
为了更好地说明本实施例，下面将举例说明插值基函数下的k
1e
：
[0073][0074]
其中：dx是单元体在x方向上的长度，dy是单元体在y方向上的长度，dz是单元体在z方向上的长度；n1、n2以及n3没有实际意义，仅是三个维度为的4
×
4方阵，值与插值基函数有关。
[0075]
相应的，采用矢量有限元法有：
[0076]
[0077]
其中：k
2e
是12
×
12的矩阵，其由三个小的矩阵(k
11
,k
22
,k
33
)和大量零元素组成，每个小的矩阵都是维度为4
×
4的方阵。
[0078]
为了更好地说明本实施例，下面将举例说明插值基函数下的k
2e
：
[0079][0080]
其中：ne是特定矩阵，其值不会随着插值基函数的变化而变化。
[0081]
用矩阵形式表示天然电磁场的控制方程为：ae＝0；
[0082]
其中：矩阵a是天然电磁场的控制方程的系数矩阵；e为模型所有棱边的电场值，所述棱边包括内部棱边ii和外部棱边ib。
[0083]
相应的，将天然电磁场的控制方程的矩阵形式表示为：a
iieii
+a
ibeib
＝0；
[0084]
其中：a
ii
是模型的内部棱边ii对应的系数矩阵；a
ib
是模型的外部棱边ib对应的系数矩阵；e
ii
是模型的内部棱边ii的待求电场；e
ib
是模型的外部棱边ib的待求电场；
[0085]
相应的，a
iieii
＝-a
ibeib
，右端项b＝-a
ibeib
。
[0086]
对于e
ib
，将模型分为上表面、下表面以及测表面：对于上表面，使用单位“1”来表示；对于侧表面，使用大地电磁的一维有限元解来表示；对于下表面，使用零来表示。
[0087]
采用节点有限元离散电流散度方程并对原始控制方程进行正则化修正：
[0088]
原始控制方程即为天然电磁场的控制方程。对于天然电磁场的控制方程的修正，理论上是使方程满足电流的散度为零的条件，即j是电流密度。本实施例充分考虑矢量有限元的特点，使得单元体内部的电流散度为零，此时，只需要考虑单元表面的积分，以此简化积分的难度。
[0089]
用模型的电导率σ乘以电场强度e，再连续求散度和梯度，并带入天然电磁场的控制方程，则获得新的电磁场控制方程：
[0090][0091]
其中：修正因子λ与模型的电导率σ相关，优选的，修正因子λ为模型的电导率σ的倒数，即
[0092]
对于天然电磁场的控制方程，当频率较低时，-iωμσe的影响较小，此时电场的解不再满足电流散度为零的条件，即导致收敛速度变慢，其中：j＝σe。用矩阵gd表示修正项的系数矩阵，由于外部棱边ib的电场值是人为给定的，故只需要对内部棱边ii的电场值进行修正，修正后的控制方程为：a
neweii
＝b；其中：a
new
是新的电磁场控制方程的系数矩阵，且a
new
＝a
ii-gd。
[0093]
因为矢量有限元的基函数无散，且每个单元体内部的电导率是常数，因此在每个单元体内部的电流的散度都为零。但是，近似电流密度在单元体面之间并不连续，因此在单元体面上的电流散度不为零。此时，对单元体面上的电荷的静态势进行计算，并从近似电场中减去静态势的梯度φ，则可以达到电流密度的散度为零的条件，即：
[0094][0095]
相应的，采用galerkin加权余量法得到加权余量方程后，用节点有限元法离散：
[0096][0097]
其中：ne是单元总数，是未知量，μh是节点有限元的插值基函数，dv是体积分。
[0098]
相应的，根据第一标量格林定理，在一个单元p内有：
[0099][0100]
其中：v
p
是单元体的表面，s
p
是单元的表面，是表面的法向向量，ds是面积分。
[0101]
一般情况下，面积分的贡献会被来自共享该面的单元的相同贡献而抵消。然而，在近似电流密度不连续的情况下，同一面的两个单元的贡献之间存在着非零的差异。
[0102]
相应的，有：
[0103][0104]
其中：是迭代求解的电场值，nk是矢量有限元的插值基函数。
[0105]
令面积分每个单元在x方向的电场分量e
x
、y方向的电场分量,ey,以及z方向的电场分量ez都对应有两个面，一个面的法向与场的方向一样，另一个相反，下面将举例说明：
[0106]
对于x方向的电场分量e
x
，单元两侧的y表面和z表面都需要计算积分。在一个单元中，我们应该计算六个面的面积分s，s为8
×
12的矩阵。
[0107]
相应的，有：
[0108]
s＝(s
hk
)8×
12
＝[s1,s2,s3]＝σ[t1,t2,t3]；
[0109]
其中：s
hk
为s的分量，且有s1表示x方向的面积分，s2表示 y方向的面积分，s3表示z方向的面积分，且每一项面积分都由同一方向的两个面积分累加而成；t1、t2以及t3均为8
×
4的矩阵，其值均不唯一，由矢量插值基函数和节点插值基函数共同决定。
[0110]
求解两种极化模式下的电磁场方程组：
[0111]
对修正后的控制方程进行迭代求解，分别获得在xy-极化模式下和yx-极化模式下各个测点的电场分量，并基于各个测点的电场分量根据如下公式求得各个测点的磁场分量h：
[0112][0113]
根据阻抗以及视电阻率、相位公式求解相应的视电阻和相位，具体包括如下步骤：
[0114]
s1、根据如下公式计算在xy-极化模式下和yx-极化模式下的阻抗：
[0115]
以及
[0116]
其中：z
xy
是在xy-极化模式下的阻抗；z
yx
是在yx-极化模式下的阻抗；
[0117]
s2、根据如下公式计算在xy-极化模式下的视电阻率和相位：
[0118]
以及
[0119]
其中：ρ
xy
是在xy-极化模式下的视电阻率；是在xy-极化模式下的相位；
[0120]
根据如下公式计算在yx-极化模式下的视电阻率和相位：
[0121]
以及
[0122]
其中：ρ
yx
是在yx-极化模式下的视电阻率；是在yx-极化模式下的相位。
[0123]
为了更好地说明本实施例的优势和目的，下面举例说明：
[0124]
参见图2，本实施例公开一种优选的含有异常体的模型示意图。
[0125]
所述模型示意图包含低阻异常体的复杂模型。此时，以异常体的正上方的地面为坐标原点，建立三维空间坐标系，异常体在x方向和y方向均从-8km到8km，z方向从-8km到
ꢀ‑
24km。因此，异常体的大小为16km
×
16km
×
16km。此时，采用等距剖分所述异常体，每个网格的大小均为2km
×
2km
×
2km。在异常体的外部，对网格进行适当的拉升，使得模拟区域尽量大，研究区域一共被剖分为70
×
70
×
64个单元，总的模型区域为160km
×
160km
×
140km。空气层的层数为8层，采用指数方式进行剖分，厚度为103km。大地的电阻率100ω
·
m，异常体的电阻率为10ω
·
m，空气层的电阻率设置为10
10
ω
·
m。
[0126]
将本实施例基于正则化修正方程的三维大地电磁数值模拟方法与传统有限元方法进行对比，对比项包括视电阻率和阻抗相位，并同时考虑xy-极化模式和yx-极化模式，计算周期设置为100s。比对结果参见图3至图6。
[0127]
图3至图6表明，本实施例采用的三维大地电磁数值模拟方法与传统有限元方法进行视电阻率和相位拟合时，均可以达到高度的拟合，因此证明了本实施例采用的三维大地电磁数值模拟方法的正确性。
[0128]
将本实施例基于正则化修正方程的三维大地电磁数值模拟方法与传统有限元方法的迭代次数和计算时间进行对比，采用xy-极化模式，计算容差为10-10
，计算周期分别为100s、200s、500s和1000s。比对结果参见图7与图8。图7与图8表明，本实施例采用的三维大地数值模拟方法的周期较传统有限元方法方法相比更稳定，且迭代次数和计算时间均优于传统有限元方法。随着计算周期的增大，本实施例采用的三维大地电磁数值模拟方法的优势更加明显。
[0129]
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。