基于空气地球分离策略的大地电磁三维正演

李袁傲, 任政勇,2,3,4*, 汤井田,2,3,4, 吴启红

1 中南大学地球科学与信息物理学院, 长沙 410083 2 中南大学有色金属成矿预测与地质环境监测教育部重点实验室, 长沙 410083 3 有色资源与地质灾害探查湖南省重点实验室, 长沙 410083 4 自然资源部覆盖区深部资源勘查工程技术创新中心， 合肥 230001 5 成都大学建筑与土木工程学院, 成都 610106

0 引言

大地电磁法(magnetotellurics，MT)是一种利用天然交变电磁场研究地球内部电性结构的地球物理方法(Cagniard,1953)，被广泛运用于矿产资源勘察、岩石圈电性结构研究以及石油天然气勘探等领域(Chave and Jones, 2012).当前，大地电磁数据采集趋向于三维方式，采集区域地质情况复杂且常具有起伏地形，因此需要研发高效的大地电磁数据三维反演方法.反演的迭代过程中需要多次调用正演算法来计算模型更新量，正演响应的准确性保证了更新模型的正确性.因此，研发快速高精度的大地电磁三维正演算法是大地电磁数据解释的研究热点.

大地电磁正演算法主要基于积分方程法，有限差分法和有限单元法.积分方程法只离散化异常体，适合于较小异常体模型的计算.由于受限于简单模型才有格林函数，积分方程法一般只能用于简单模型电磁场响应的计算与分析(Hohmann，1975; Lajoie and West，1976; Zhdanov et al.，2006; 王若等，2009).有限差分法采用差分来近似微分，具有方法简单、编程容易等优点，在积分方程法之后迅速受到了普及，被广泛用于三维地电磁感应问题的求解 (Mackie et al.，1994; Newman and Alumbaugh，2002; Siripunvaraporn et al.，2002; 沈金松，2003; 谭捍东等，2003; Streich，2009; 李展辉等，2009; Dong and Egbert，2019) .但是，由于差分近似一般适用于直线、不适合于复杂折线和曲线，有限元差分方法在处理任意起伏地形时存在一定的困难.有限元法采用在小区域内利用简单多项式来近似复杂未知场的思想，在小区域几何形状的选取上具有任意性，可以精确逼近三维复杂地质情况以及起伏地形，从而计算出高精度的电磁场响应，因此有限单元法逐渐成为了当前地电磁感应正演问题的主流算法.

Coggon(1971)率先基于有限单元法计算了地球电磁场响应.随后，地电磁感应问题的有限元计算方法在网格离散化技术、线性方程组求解技术、并行加速技术、后处理方法与技术等方面取得了长足进步(Raiche，1994; Avdeev，2005; Börner，2010; Newman，2014).与此相对的是，在微分控制方程最优化设计方向的研究成果较少.目前，地电磁感应问题的有限元计算一般采用电场方程、磁场方程和势场方程.利用电场方程，Nam等(2007)和Farquharson和Miensopust (2011)进行了大地电磁问题的有限元计算；Ren等(2013)将面向目标自适应有限元法应用到大地电磁三维正演计算中.利用磁场方程，Siripunvaraporn等(2002)和Franke等(2007)进行了大地电磁三维场分布的计算.势场方程方面，Haber(2000)基于A-φ势实现了大地电磁三维有限元模拟；Mitsuhata和Uchida(2004)实现了基于T-Ω势的大地电磁三维有限元模拟.

空解现象产生的本质原因为在空气和地球中采用了同样的控制方程，由于空气和地下具有完全不同的电性特征(电导率对比度为无穷大)，采用同样的控制方程必然导致病态特征.为解决此问题，本文提出了一种新的空气与地球分离混合公式.解决思路为：首先，将求解区域分解为空气、地下两个独立的区域，分别采用不同的控制方程；然后，利用地空分界面上电磁场连续性条件，进行界面两侧的耦合与信息交换，保证电磁场解的唯一性.基于该解决思路，空气层中磁场H分解为电矢量势T和磁标量势φ，利用Coulomb规范条件，使得矢量势T和标量势φ满足具有Laplace结构的势场方程.地下区域中，采用双旋度结构的电场方程.在地下区域采用电场双旋度方程的原因在于：地下介质的电导率大于零，这使得双旋度系统在地下区域避免了零空间问题且有良好的稳定性.

最后，我们计算了国际标准测试模型COMMEMI 3D-1、DTM-1与梯形山模型的大地电磁响应，将计算结果与其他算法进行了对比，来验证本文算法的正确性，以及处理起伏地表问题的能力.我们还进行了系统条件数与特征值的计算和对比分析，来检验本文算法对电场双旋度系统病态性的改进.

1 节点矢量混合有限元系统

利用大地电磁的场源位于求解区域外的一般性假设，在求解区域内(见图1)，大地电磁问题的电磁场响应满足频率域Maxwell方程组(选取时谐因子为e-iω t)：

(1)

其中，E为电场矢量，H为磁场矢量，i为纯虚数单位，ω为角频率，μ为磁导率，

(2)

其中，J为电流密度，σ为电导率，ε为介电常数，

(3)

(4)

图1 MT地电模型示意图Ω0表示空气区域，Ω1表示地下区域，Γ表示地空边界，Γ0表示空气外部边界，Γ1表示地球外部边界，n0为空气区域外边界的外法向向量，n1为地球区域外边界的外法向向量，n01为地空边界上空气一侧的外法向向量，方向由空气指向地球，n10为地空边界上地球一侧的外法向向量，方向由地球指向空气，n为约定的地空边界上法向向量，方向由地球指向空气.Fig.1 Schematic diagram of MT geoelectric modelΩ0 represents the air area, Ω1 represents the underground area, Γ represents the air-earth interface，Γ0 represents the outer boundary of air，Γ1 represents the outer boundary of earth，n0 is the outer normal vector of the outer boundary of the air region, n1 is the outer normal vector of the outer boundary of the earth region, n01 is the outer normal vector of the air side on the air-earth interface, and the direction is from the air to the earth, n10 is the outer normal vector on the side of the earth on the air-earth interface, and the direction is from the earth to the air, n is the normal vector on the agreed air-earth interface, with the direction from the earth to the air.

1.1 空气层电磁场积分弱形式

(5)

将公式(5)代入公式(2)中得到：

(6)

其中H0表示空气中的磁场.引入标量函数φ的负梯度来表示公式(6)中的矢量部分，可得到：

(7)

将公式(7)代入公式(2)并取旋度得到：

即：

(8)

将公式(8)中的矢量势T的双旋度结构进行化简，强加Coulomb规范条件来唯一化矢量势T：

(9)

(10)

(11)

综上，空气区域的控制方程为

(12)

为了求解公式(12)，采用Galerkin节点有限元法(Jin，2002)，先假设余量rT为

(13)

余量rφ为

(14)

采用线性插值基函数，四面体单元中任意点处的矢量势T以及标量势φ可由节点基函数表示为

(15)

其中Txi，Tyi，Tzi分别为矢量势T的三个分量，Li为节点基函数，在空气区域中满足以下条件：

∭Ω0Li·rTdv=0,

∭Ω0Li·rφdv=0.

(16)

利用Green恒等式：

(17)

(18)

综合公式(15)—(18)，给出公式(12)的积分弱形式：

高温环境下GH3536基体软化，导致其粘着趋势上升，存在外部载荷时易发生塑性变形，在磨球滑动过程中以剥层的方式不断脱落的磨屑被胶合在一起，形成层状物(见图5a)并逐渐被撕裂。另外，基体磨痕的局部区域可见沿滑动方向的犁沟，这可能是由于部分磨屑在磨损过程中嵌入到基体中，在滑动中推挤基体使之塑性流动并犁出沟槽，从而引起磨粒磨损。由图5b可见，NiAlW涂层的磨痕比基体更为光滑和平整，没有明显的单个颗粒沉积物脱落和涂层整块粘着或撕裂现象。涂层的磨损机制以微切削为主，其表面的单个颗粒沉积物凸起部分首先被切削，在外界的反复撞击与挤压下发生破碎，最后脱落并造成涂层的磨损[15-16]。

(19)

写成分量形式为

(20)

公式(19)与公式(20)即为空气中的控制方程积分弱形式.

1.2 地球电磁场积分弱形式

地下区域中，电场满足双旋度方程：

(21)

采用Galerkin矢量有限元法(Jin，2002)，求解公式(21).先假设余量RE为

(22)

采用棱边基函数，四面体单元中的电场可表示为

(23)

其中E1i为对应棱边上的电场切向分量，Ni为棱边基函数，计算区域内满足以下条件：

∭Ω0Ni·REdv=0.

(24)

(25)

其中a,b分别为矢量函数，则

(26)

(27)

(28)

公式(28)即为地下介质区域控制方程的积分弱形式.

1.3 空气-地球耦合系统

结合公式(19)和公式(28)以及图1中所定义的边界外法向向量，空气地球耦合系统的积分弱形式为

(29a)

(29b)

(29c)

根据矢量公式：

(30)

公式(29a)可以进一步改写为

(31)

将公式(29)中与地空界面有关的项移到等式左边，且在地空界面上有n=n10=-n01，可以得到：

(32)

(33)

在地空分界面上施加如下电磁场连续性条件:

(34)

(35a)

(35b)

(35c)

公式(35a)—(35c)中的右端项均可由1D层状介质的解析解给出，将公式(35)写成矩阵形式：

(36)

其中：

公式(35)和公式(36)即为空气地球分离策略的有限元离散方程.

其中i,j=x,y(x≠y)，Zij为阻抗，ρij为视电阻率，φij为相位，E1i和H1j分别为电场E1、磁场H1在对应方向上的分量.

2 数值试验

本小节所有计算均在CPU主频为2.30 GHz，RAM超过200 GB的小型AMAX工作站上完成.

2.1 COMMIME 3D-1模型

3D-1模型是国际标准测试COMMEMI项目(Zhdanov et al., 1997)的一个模型.它是均匀半空间中嵌套一个导电棱柱，其顶埋深为250 m，背景电阻率为100 Ωm，棱柱电阻率为0.5 Ωm，尺寸为1000 m×2000 m×2000 m，空气电阻率为1016Ωm,计算区域尺寸为70 km×70 km×140 km.

计算网格包含148576个四面体单元.传统双旋度方程共173335个自由度，本文提出的地空分离公式共168978个自由度，频率为0.1 Hz.我们将测试结果(Air-Earth decomposition)与高精度自适应有限元计算结果(Curl-Curl system)(Ren et al., 2013)进行对比，图2展示了视电阻率和相位的对比，图3展示了这两种方法的相对误差.空气地球分离公式计算用时 19.0 s，电场双旋度公式计算用时 21.8 s.

图2和图3显示混合公式系统与双旋度公式系统能够获得一致的视电阻率与相位曲线.视电阻率值从测线端点到y=0处逐渐降低，过y=0之后逐渐增大，曲线具有对称性.视电阻率最大相对误差0.34%，相位最大相对误差0.18%，验证了算法的正确性和精度.为了进一步研究混合系统的数值特征，我们计算并对比了特征值分布(如图4和图5).

图2 3D-1模型视电阻率与相位，测试频率为0.1 Hz, Curl-Curl system代表了自适应有限元结果(Ren et al., 2013),Air-Earth decomposition 代表了本文计算的结果Fig.2 Comparison of calculated apparent resistivity and phase on the 3D-1 model among our solution (denoted by Air-Earth decomposition) and the reference solution (Curl-Curl system, Ren et al., 2013), at a frequency of 0.1 Hz

图3 3D-1模型视电阻率与相位相对误差(高精度自适应有限元计算结果为参考值(Ren et al., 2013))，测试频率为0.1 HzFig.3 Relative error of apparent resistivity and phase by comparing to the reference solution (Curl-Curl system,Ren et al., 2013) for the 3D-1 model at a frequency of 0.1 Hz

图4表明两个系统的最大特征值分布具有相似性.图5表明地空分离公式系统有更多的非零或模长更大的最小特征值.当一个系统的特征值全部位于复数平面的左半部分时，系统是收敛的.与此相反，当特征值全部位于复数平面的右半部分时，系统则是发散的.从图4与图5可以直观的看出，空气地球分离公式的特征值更多地集中在复数平面的左半部分且距离原点更远，因而空气地球分离公式相对于电场双旋度公式有更好的稳定性.根据条件数的定义(最大特征值与最小特征值之比)，地空分离系统拥有更小的条件数.为了定量研究地空分离系统的条件数，计算了不同频率下的条件数(从0.001 Hz，0.01 Hz，…，10 Hz共六个对数等间隔分布的频点).测试结果如图6所示，结果表明混合系统的条件数相对于双旋度系统均降低，下降可达3个数量级.

图4 3D-1模型上地空分离公式和双旋度电场公式对应的前1000个最大特征值分布的对比分析，测试频率为0.1 HzFig.4 The distribution of 1000 maximum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the 3D-1 model, at a frequency of 0.1 Hz

图5 3D-1模型上地空分离公式和双旋度电场公式对应的后1000个最小特征值分布的对比分析，测试频率为0.1 HzFig.5 The distribution of 1000 minimum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the 3D-1 model at a frequency of 0.1 Hz

图6 3D-1模型地空分离公式和双旋度电场公式的条件数对比Fig.6 Comparison of condition numbers for two methods (one using the Air-Earth decomposition, and one using the Curl-Curl system) on the 3D-1 model at different frequencies

2.2 梯形山模型

为了测试对地形的适应性，计算了梯形山模型(Nam et al., 2007).如图7所示，梯形山顶部为450 m×450 m的正方形，底部为2000 m×2000 m的正方形且高为450 m.地形电阻率为100 Ωm，空气的电阻率为1016Ωm，测线方向为x=0方向，计算区域尺寸为25 km×25 km×50 km.计算区域包含186168个四面体单元，双旋度系统共227241个自由度，地空分离系统共219957个自由度，测试频率为2 Hz.地空分离公式和双旋度电场公式的模型视电阻率与相位曲线对比如图8所示，相应的相对误差见图9.空气地球分离公式计算用时 54.26 s，电场双旋度公式计算用时 59.66 s.

图7 梯形山模型Fig.7 The trapezoidal hill model

图8 梯形山模型视电阻率与相位，测试频率为2 Hz，Curl-Curl system代表了自适应有限元结果(Ren et al., 2013),Air-Earth decomposition代表了本文计算的结果Fig.8 Comparison of calculated apparent resistivity and phase on the trapezoidal hill model among our solution (denoted by Air-Earth decomposition) and the reference solution (Curl-Curl system, Ren et al., 2013) at a frequency of 2 Hz

测试结果显示：视电阻率最大相对误差0.25%，相位最大相对误差0.50%.在视电阻率曲线中，ρxy在梯形山基底部位的电阻率稍高于围岩的电阻率，而在起伏地形以及山顶部分则低于围岩电阻率，ρyx则呈现两端高中间低的电阻率变化趋势，这是由于电流的地形效应引起的.xy模式中的电流沿y方向(垂直测线方向)传播并在山体的侧面上向山顶弯曲，使得沿测线方向山脚处的电流发散而山脊和山顶处电流聚集，导致了山脚的电阻率高于围岩而山脊和山顶的电阻率低于围岩.在yx模式中，电流沿x轴方向(沿测线方向)传播，电流的传播路径发生弯曲且都在x=0方向上汇聚，因而在yx模式中，ρyx低于围岩电阻率.

为了研究混合系统与双旋度系统的数值特征，我们计算了特征值分布，测试结果如图10和图11.与3D-1模型类似，图10表明两个系统具有相似的最大特征值分布，图11表明地空分离系统拥有模长更大的最小特征值，且地空分离系统的特征值更加集中在复数平面的左半部分且距离原点更远.

图10 梯形山模型上地空分离公式和双旋度电场公式对应的前1000个最大特征值分布的对比分析，测试频率为1 HzFig.10 The distribution of 1000 maximum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the trapezoidal hill model at a frequency of 1 Hz

图11 梯形山模型上地空分离公式和双旋度电场公式对应的后1000个最小特征值分布的对比分析，测试频率为1 HzFig.11 The distribution of 1000 minimum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the trapezoidal hill model at a frequency of 1 Hz

我们在这个模型中选取0.01 Hz，0.1 Hz，…，100 Hz共五个频点测试了系统的条件数.如图12所示，模型的条件数分布与COMMEMI 3D-1模型的条件数分布特征一致，条件数有1～3个数量级的改善，在10 Hz频率下两系统条件数差异最大.

图12 梯形山模型地空分离公式和双旋度电场公式的条件数对比Fig.12 Comparison of condition numbers of the Air-Earth decomposition and the Curl-Curl system on the trapezoidal hill model at different frequencies

2.3 DTM-1模型

该模型由嵌套在均匀半空间中的三个异常体组成，其模型示意图如图13，其中ρ1=10 Ωm,ρ2=1 Ωm,ρ3=10000 Ωm，介质电阻率最大对比度可达10000∶1，模型外边界距离地面250 km.计算区域尺寸为250 km×250 km×500 km，测线沿y=0方向分布.

图13 DTM-1模型结构Fig.13 Structure of DTM-1 model

计算区域包含289178个四面体单元，双旋度系统共343978个自由度，地空分离系统共343615个自由度.选取测试频率为0.01 Hz进行计算，空气地球分离公式计算用时 70.47 s，电场双旋度公式计算用时 71.86 s.模型的电磁响应曲线如图14所示，相对误差如图15所示，结果显示两种系统的响应曲线吻合度十分高，视电阻率最大相对误差0.20%，相位最大相对误差0.43%.与梯形山模型视电阻率曲线分析类似，该模型的视电阻率的变化趋势依旧可以根据电流的传播方向来分析.由于ρ3=10000 Ωm高于围岩电阻率，使得通过异常体的电流会向周围发散，则在x=[-25000 m,0]的范围内电流向测点聚集，反应在视电阻率曲线上则是这一段范围内的视电阻率呈现下降趋势.相反，ρ1=10 Ωm,ρ2=1 Ωm低于围岩电阻率，电流向异常体聚集，所以在x=[0,25000 m]范围内测点处的电流发散，视电阻率升高.

图14 DTM-1模型视电阻率与相位，测试频率为0.01 Hz，Curl-Curl system代表了自适应有限元结果(Ren et al., 2013),Air-Earth decomposition代表了本文计算的结果Fig.14 Comparison of calculated apparent resistivity and phase on the DTM-1 model among our solution (denoted by Air-Earth decomposition) and the reference solution (Curl-Curl system, Ren et al., 2013) at a frequency of 0.01 Hz

图15 DTM-1模型电磁响应相对误差(高精度自适应有限元计算结果为参考值(Ren et al., 2013))，测试频率为0.01 HzFig.15 Relative error of our apparent resistivity and phase by comparing to the reference solution (Curl-Curl system,Ren et al., 2013) for the DTM-1 model at a frequency of 0.01 Hz

同样，我们也在该工作频率下测试了两个系统的特征值分布，如图16和图17所示,DMT-1模型的特征值分布特征与3D-1模型以及梯形山模型的测试结果相似，空气地球分离公式系统矩阵的特征值更多地分布在复数平面左半部分距离原点更远处.我们在这个基础上选择了从0.001 Hz，0.01 Hz，…，100 Hz共六个对数等间隔分布的频点进行了条件数的测试，结果如图18所示.根据图18可以得知，混合系统对双旋度系统的病态性做出了改进，在1 Hz的工作频率下，条件数的改善最明显，可达3个数量级.

图16 DTM-1模型上地空分离公式和双旋度电场公式对应的前1000个最大特征值分布的对比分析，测试频率为0.01 HzFig.16 The distribution of 1000 maximum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the DTM-1 model at a frequency of 0.01 Hz

图17 DTM-1模型上地空分离公式和双旋度电场公式对应的后1000个最小特征值分布的对比分析，测试频率为0.01 HzFig.17 The distribution of 1000 minimum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the DTM-1 model at a frequency of 0.01 Hz

图18 DTM-1模型地空分离公式和双旋度电场公式的条件数对比Fig.18 Comparison of condition numbers of the Air-Earth decomposition and the Curl-Curl system on the DTM-1 model at different frequencies

3 结论

我们推导了一种新的基于空气地球分离策略的大地电磁混合计算公式，在空气中采用Laplace算子取代双旋度算子，将双旋度系统分解为Helmholtz-Curl耦合系统，有效改善了系统矩阵在低频条件下的病态性特征.

模型测试验证了新系统的正确性与精度.在相同网格的条件下，本文算法的自由度数低于电场双旋度方程的自由度数.在计算速度上，空气地球分离公式稍高于电场双旋度公式.特征值分布与条件数分布表明新系统具有更好的特征值分布与更低的条件数，在数值上更加稳定，对双旋度系统的病态性做出了有效的改进，因而新系统应更适用于低频条件下的大地电磁正演工作.

本研究基于空气地球分离思路推导了一种具有良好条件数的混合公式，为寻求大地电磁最优化控制方程开启了新的求解思路.基于新思路，我们期待更多、性能更佳的大地电磁混合公式.