汤文武 柳建新 叶益信 张 华(东华理工大学地球物理与测控技术学院，江西南昌 330013； 中南大学地球科学与信息物理学院，湖南长沙 410083)

1 引言

电磁法是根据电磁感应原理研究天然电磁场或人工(可控)场源在大地中激励的交变电磁场分布， 并由观测到的电磁场分布研究地下目标体的电性及地质特征的一种地球物理方法[1]，已被广泛地应用于金属矿勘探[2，3]、地下水资源调查[4]、油气田勘探[5]、油气动态监测[6]等领域。根据场源的类型可分为天然源及可控源两类，现今可控源电磁法发展迅速，但场源的引入给其资料处理和解释带来了诸多困难。作为资料解释的重要手段，正反演技术的发展对进一步推动可控源电磁法的广泛应用具有重要作用，尤其是三维正、反演方法，已成为研究热点。

针对三维可控源电磁法的正演问题，近年来国内外学者开展了大量研究。Newman等[7]采用基于二次电场的Helmholtz 方程，将全空间或层状半空间电磁场响应作为一次场引入到方程右端，避免了直接模拟电偶源或磁偶源的困难。利用交错网格有限差分法对Helmholtz方程进行离散，得到的复系数线性方程组通过雅可比预处理的QMR(quasi-minimal residual method)迭代法进行求解，提高了正演效率。Weiss等[8]利用有限差分法对三维各向异性介质的电磁感应问题进行了正演模拟，并采用全张量表示电导率参数，更有利于解决实际问题。Streich[9]利用MUMPS直接求解器对有限差分法离散频率域可控源电磁方程得到的线性方程组进行求解，对于中等规模的求解问题精度高，且适用于多激发源的正演计算。Badea等[10]利用有限单元法计算了三维可控源的电磁响应，通过二次耦合势的引入消除了源的奇异性，并通过局部网格加密提高了求解精度。Stalnaker[11]利用基于非结构化四面体网格的有限元对三维可控源电磁进行正演模拟，采用基于二次耦合势的方式避免了直接模拟场源。基于非结构化的四面体网格能够适应更复杂的地形及地质体结构，且模拟精度很高。Farquharson等[12]采用矢量有限元进行大地电磁三维正演，并在各个单元接触面上施加散度校正以提高正演计算速度。Puzyrev等[13]利用节点型有限单元法求解基于二次耦合势的三维可控源电磁方程，并考虑了电导率的各向异性。Tong等[14]利用有限单元法对三维大地电磁法进行了正演模拟，通过将散度校正引入到有限元方程中，避免了伪解的出现。徐志峰等[15]基于磁矢量势及电标量势将Maxwell方程转化为位势的类Helmholtz方程，通过引入罚项及稳定化方法避免了伪解及数值不稳定，采用有限元方法计算了三维可控源的电磁响应。另外，付长民等[16]、殷长春等[17]、王超[18]、韩波等[19]、汤文武等[20]、严波等[21]、李建慧等[22]分别对三维可控源电磁正演问题开展了研究。

有限单元法由于其对地形及不规则异常体具有很好的适应性受到广泛关注[23]。对于有限单元法在三维可控源电磁中的应用，早期主要采用节点型有限元，其概念及实现相对简单，但可能存在伪解问题。近年来，矢量有限元法逐步获得更多关注，究其原因在于矢量有限元在电导率分界面上自动满足电磁场的边界条件，可消除伪解。但目前并未见两种正演方法实际应用的详细对比。因此，基于前人研究成果，本文开展关于节点有限元及矢量有限元在三维频率域可控源电磁正演计算的应用对比，旨在为两种方法的合理选用提供参考。

2 理论基础

频率域可控源电磁法的三维正演可基于电场开展，亦可依据耦合势进行，文中的三维正演是基于二次电场的有限元正演，可避免源的奇异性，与总场正演方式相比，同等条件下具有可采用更小网格的优势。基于节点有限元与矢量有限元的三维正演，其基本方程是一致的，但具体处理边值问题时会有所差异。文中矢量有限元三维正演实现是基于笔者已发表的成果[20]，为了叙述完整，以下简要讨论其基本原理。

2.1 基于二次电场的偏微分方程

电磁场理论是可控源电磁法的基本理论。频率域中，在准静态条件下(忽略位移电流)，电磁场满足以下麦克斯韦方程组

(1)

(2)

式中：E为电场矢量(单位为V/m)；H为磁场矢量(单位为A/m)；Jc为外电流密度矢量(单位为A/m2)；μ0为真空中的磁导率；σ为地下电性介质的电导率。在式(1)和式(2)中，时谐因子设为e-iω t。联立方程式(1)和式(2)，可得到以下关于电场矢量的双旋度方程

(3)

为了消除场源引起的奇异性，常将总场分解为一次场和二次场。通过预先选定某个简单参考模型，利用计算得到的一次场Ep求解二次场Es，并最终得到总场，即

E=Ep+Es

(4)

对于一次场Ep，同样可得以下双旋度方程

(5)

式中σp为预先选定的参考模型的电导率分布，一般可选为全空间、均匀半空间或水平层状介质模型以简化一次场的计算。用式(3)减去式(5)，并化简，得到

(6)

式中σa为异常电导率分布，且有σa=σ-σp。

2.2 边界条件

采用一次场、二次场分离的方式计算总场时，边界条件可以取得比较简单。当边界距离目标研究区域足够远时，二次场可近似为零。因此，采用如下狄利克雷边界条件

n×Es=0

(7)

式(6)、式(7)便构成了可控源电磁三维正演的边值问题。这里，选定均匀半空间模型计算一次场，其水平电偶极的谐波场存在解析积分公式[24]，具体表达式见附录A。积分表达式中都是含有第一类贝塞尔函数J0或J1项的汉克尔变换，可利用数字线性滤波器快速求解，本文采用的是Guptasarma等[25]提出的数值滤波算法。

3 有限元分析

对于上节建立的边值问题，可采用节点有限元或矢量有限元方法求解。为简单起见，这里仅以长方体剖分单元为例展开讨论。当采用节点有限元时，电磁场分量置于节点上采样，而矢量有限元则将电磁场分量置于各个单元的棱边上采样(图1)。图1a为节点有限元中各个电场分量的分布，每个节点上都布置有电场三分量，因此每个单元有24个自由度；而图2b为矢量有限元各个电场分量的布置，每条棱边上都有与之平行的电场分量，每个单元的自由度为12。

图1 电场分量在节点有限元(a)和矢量有限元(b)的剖分单元中的布置示意图

对于以上边值问题，基于变分原理可推导相应的变分问题如下

iωμ0σaEs·Ep]dV

(8)

(9)

式中λ是调节散度项在以上泛函表达式的比重因子，当取值为1时可获得较满意结果。

当采用矢量有限元时，由于矢量元在棱边上取样，在各个单元内自动满足散度条件，因此不需要在各个剖分单元内施加散度条件。但进一步研究发现，基于矢量有限元求解时获得的有限元方程的刚度矩阵条件数很大，采用迭代法收敛很慢，原因是在不同单元接触面上可能不满足散度条件，因此有必要进一步校正以加速迭代法求解过程。

4 模型实例对比

为了对比节点有限元与矢量有限元在三维可控源电磁中的模拟效果，设计了两个模型开展正演计算，包括一个存在解析解的水平层状介质模型及一个无解析解的低阻体模型。所有模拟结果都是设定两种方法在同一剖分网格及线性方程组求解器的情况下获得的。采用的计算平台(laptop)为：Windows 7；CPU，Intel i7-3610qm；8GB RAM。

4.1 水平层状介质模型

水平层状介质模型的电磁响应具有解析表达式，这里选用如图2所示的三层水平介质模型对正演计算结果进行对比验证。场源取为单位电偶源并置于原点处，分别计算1～8192Hz按对数等间隔分布的14个频率的电磁响应。采用的网格为91×76×50，其中z方向包含20个空气层单元和30个地下介质层单元。图3展示了坐标为(0，-2300m，0)的测深点处卡尼亚视电阻率及相位的解析解、节点有限元及矢量有限元卡尼亚视电阻率及相位的数值解，其中图3a为卡尼亚视电阻率曲线的对比，图3b为相位曲线的对比。

图2 水平三层介质模型

由图3a定性分析可知：矢量有限元获得的卡尼亚视电阻率曲线与解析解吻合很好，而节点有限元获得的卡尼亚视电阻率数值大体上与解析解吻合，但在个别低频及高频端数据误差较大。进一步由计算结果统计可知，对于卡尼亚视电阻率而言，矢量有限元数值解所有频点数据误差在1%以下，平均误差为0.17%；节点有限元数值解高频、低频段数据误差大，中间频段数据误差小，基本在3%以下。对图3b做类似分析可知，矢量有限元数值解整体上与解析解吻合较好，平均误差为4.4%；而节点有限元在低频、高频段误差大，中间频段数据误差小，基本在5%以下。

进一步通过模拟实验可知，缩小中间网格剖分间距可以改善节点有限元在高频段误差大的现象，因此高频段的误差大主要跟网格剖分有关，相比矢量有限元而言需要更密的网格才能达到合理的精度。但低频段的误差通过网格剖分控制基本无法进一步改善，可能的原因应该是与节点有限元方法本身有关。

由于本文算法是基于二次场开展的，因此，为了更准确地对比矢量有限元与节点有限元的计算精度，有必要直接针对二次场进行对比。图4所示是坐标为(0，-2300m，0)的测深点处的二次电场Exs的对比图，其中二次场的解析解由图2中三层介质模型及另一个电阻率为100Ω·m的均匀半空间模型的解析解之差获得。图4a代表二次电场分量Exs实部绝对值(以|Re(Exs)|表示)的对比，图4b代表二次电场分量Exs虚部绝对值(以|Im(Exs)|表示)的对比。由图可知，矢量有限元及节点有限元求解获得的二次场与解析解基本吻合，但相比而言，矢量有限元精度更高。

图3 卡尼亚视电阻率(a)及相位曲线(b)对比图

图4 二次电场实部分量(a)和虚部分量(b)曲线对比图

4.2 单一低阻体模型

以一个单一低阻体模型为实例，分别进行基于节点有限元及矢量有限元法的三维可控源电磁响应模拟，并作对比。图5是其断面(图5a)及平面(图5b)示意图。该模型描述了一个置于电阻率为100Ω·m的均匀半空间的长方体低阻体，激发源为单位电偶源，距离低阻体的最大水平距离为8.8km。该低阻体长、宽、高分别为0.6、0.5、0.2km，电阻率为10Ω·m。

统一采用91×116×50的网格对该模型进行正演计算，图6分别给出了采用节点有限元及矢量有限元法的模拟结果。图6a为节点有限元视电阻率不同频率下的三维等值线切片图； 图6b为矢量有限元法正演结果对应的切片。5个不同切片分别对应64、32、16、8、4Hz。由频率测深原理可知，5个不同切片分别对应由浅至深的地下介质的电阻率分布及变化情况。由于采用相同色标，因此颜色的变化即体现了视电阻率的变化。从图中等值线分布情况来看，两种正演方法得到的不同频率下的视电阻率分布形态基本一致，因此相互验证了各自正演算法的正确性。

具体分析各张切片图：64Hz切片，尽管从等值线分布看两种方法差别较大，但实际上视电阻率基本等于均匀大地的电阻率(100Ω·m)，此时趋肤深度最小，低阻体的影响基本可忽略； 从32Hz逐步过渡到8Hz，随着趋肤深度增大，低阻体影响显现出来，切片图上可见蓝色低值区，视电阻率最低约为86Ω·m； 4Hz对应的切片未见明显低值区，视电阻率普遍高于100Ω·m，此时已基本进入近区。因此，分析不同频率下切片图的变化也在一定程度上验证了正演算法的正确性。

对于高维电磁正演计算，有限元分析得到的往往是一个大型、稀疏、对称复系数线性方程组，其求解是正演计算重点考虑的一个问题。采用预处理的QMR迭代法求解线性方程组时，统计了两种不同算法在不同频率的迭代残差。如图7所示，图7a为节点有限元方程在不同频率下的正演收敛曲线，图7b为矢量有限元方程在不同频率下的正演收敛曲线。

图5 低阻体模型断面(a)和平面(b)示意图

图6 节点有限元法(a)和矢量有限元法(b)模拟的不同频率下视电阻率三维等值线切片图

图7 不同频率下的迭代收敛曲线图(a)节点有限元法； (b)矢量有限元法

分析图7a可知，不同频率的迭代收敛曲线基本一致，大约80次迭代可使得相对误差降到10-6；分析图7b可知，不同频率下线性方程组求解的收敛曲线区别较大，总体规律是：频率越低，迭代收敛越慢，达到同等精度需要更多迭代次数。因此，当采用迭代法求解有限元方程时，矢量有限元正演计算速度受频率影响较大，且迭代收敛较节点有限元明显要慢。表1给出了两种正演计算在各个频点的计算时间，总体而言，节点有限元在不同频率下计算时间基本一致，而矢量有限元低频时正演计算时间长。值得注意的是，尽管从自由度的角度来分析，矢量有限元的自由度小于节点有限元的自由度，同等条件下节点有限元离散得到的刚度矩阵更大。但由于节点有限元线性方程组迭代求解收敛较矢量有限元快，矢量有限元总体正演时间较节点有限元长，平均到各个频率基本上矢量有限元计算速度较节点有限元慢一半。

表1 各频率下正演计算时间对比

5 结论

本文对节点有限元和矢量有限元可控源电磁三维正演问题进行对比、分析，在控制并确保网格剖分一致并采用相同线性方程组求解器的情况下，得出以下认识和结论：

(1)同等情况下，矢量有限元求解精度较节点有限元精度高，在全频率都能取得满意的计算结果；而节点有限元在低频段正演计算精度不高；

(2)采用迭代法求解有限元方程时，矢量有限元正演计算速度跟频率有较大关系，通常来说，频率越低，正演速度越慢；而强加散度条件的节点有限元正演计算速度与频率关系不大；

(3)同等情况下，节点有限元正演计算速度较矢量有限元法快1倍以上。

总之，在计算效率优先的情况下可考虑选择节点有限元开展正演计算；当计算精度是第一指标时应当选用矢量有限元开展正演计算。

附录A 水平电偶极源的谐波场

在水平电偶极源激发下，均匀半空间模型在地面上的电磁场响应可利用已有的积分公式计算，但由于基于二次场的三维正演模拟要利用地下空间的一次场，因此还需另外推导相应的积分公式。

考夫曼等[24]通过频率域和时间域电磁测深法推导了均匀半空间模型上水平电偶极源在地面上的电磁场响应。基于此，本文进一步推导了地下空间的电磁场响应的积分表达式。

图A-1为一个均匀半空间模型，x轴沿正北方向，y轴沿正东方向，z轴竖直向下。时谐电流为I，长为dl的水平电偶极子位于坐标原点。均匀半空间的电导率为σb，设电磁场待求点为P，P与原点的距离为r，r与x轴的夹角为φ。假设时谐因子为e-iω t。均匀半空间电磁场矢量可由一个矢势及一个标势来表示(不考虑位移电流作用)

(A-1)

(A-2)

图A-1 水平电偶极激发下的均匀半空间模型

(A-3)

(A-4)

且电场矢量可用矢势单独表示

(A-5)

依据文献[24]的推导，假设Ay=0并求解式(A-4)，可得地下空间中矢势的其他二分量积分表达式

(A-6)

(A-7)

(A-8)

采用圆柱坐标系，将式(A-6)～式(A-8)代入式(A-5)并整理，可得电场三分量表达式

(A-9)

(A-10)

(A-11)

同理，结合式(A-1)、式(A-6)和式(A-7)亦可得到磁场三分量表达式。

参考文献

[1] 汤井田，何继善．可控源音频大地电磁法及其应用．湖南长沙：中南大学出版社，2005.

Tang Jingtian，He Jishan．The Controlled-Source Audio-frequency Magnetotellurics and Its Application．Central South University Press，Changsha，Hunan，2005.

[2] Nabighian M，Corbett J．Electromagnetic Methods in Applied Geophysics：Theory．SEG，1988.

[3] 杨怀杰，潘和平，孟庆鑫等．导电围岩对井中三维瞬变电磁响应的影响规律研究．石油物探，2016，55(2)：288-293.

Yang Huaijie，Pan Heping，Meng Qingxin et al．Influence laws of conductive host on borehole 3D transient electromagnetic responses．GPP，2016，55(2)：288-293.

[4] 柳建新，麻昌英，孙丽影等．可控源音频大地电磁测深法在地热勘探中的应用．工程地球物理学报，2014，11(3)：319-325.

Liu Jianxin，Ma Changying，Sun Liying et al．The application of CSAMT to the geothermal exploration．Chinese Journal of Engineering Geophysics，2014，11(3)：319-325.

[5] He Z，Liu X，Qiu W et al．Mapping reservoir boundary by borehole-surface TFEM：Two case studies．The Leading Edge，2005，24(9)：896-900.

[6] Wirianto M，Mulder W A，Slob E C．A feasibility study of land CSEM reservoir monitoring in a complex 3-D model．Geophysical Journal of the Royal Astronomical Society，2010，181(2)：741-755.

[7] Newman G A，Alumbaugh D L．Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences．Geophysical Prospecting，1995，43(8)：1021-1042.

[8] Weiss C，Newman G．Electromagnetic induction in a generalized 3D anisotropic earth，Part 2：The LIN preconditioner．Geophysics，2003，68(3)：922-930.

[9] Streich R．3D finite-difference frequency-domain mode-ling of controlled-source electromagnetic data：Direct solution and optimization for high accuracy．Geophy-sics，2009，74(5)：F95-F105.

[10] Badea E A，Everett M E，Newman G A et al．Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials．Geophy-sics，2001，66(3)：786-799.

[11] Stalnaker J L．A Finite Element Approach to the 3D CSEM Modeling Problem and Applications to the Study of the Effect of Target Interaction and Topo-graphy [D]．Texas A & M University，2004.

[12] Farquharson C G，Miensopust M P．Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction．Journal of Applied Geophysics，2011，75(4)：699-710.

[13] Puzyrev V，Koldan J，Puente Jdl et al．A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling．Geophysical Journal International，2013，193(2)：678-693.

[14] Tong X Z，Liu J X，Xie W et al．Three-dimensional forward modeling for magnetotelluric sounding by finite element method．Journal of Central South University of Technology，2009，16(1)：136-142.

[15] 徐志锋，吴小平．可控源电磁三维频率域有限元模拟．地球物理学报，2010，53(8)：1931-1939.

Xu Zhifeng，Wu Xiaoping．Controlled source electromagnetic 3-D modeling in frequency domain for finite element method．Chinese Journal of Geophysics，2010，53(8)：1931-1939.

[16] 付长民，底青云，王妙月．海洋可控源电磁法三维数值模拟．石油地球物理勘探，2009，44(3)：358-363.

Fu Changmin，Di Qingyun，Wang Miaoyue．The 3D numerical modeling for MCSEM．OGP，2009，44(3)：358-363.

[17] 殷长春，贲放，刘云鹤等．三维任意各向异性介质中海洋可控源电磁法正演研究．地球物理学报，2014，57(12)：4110-4122.

Yin Changchun，Ben Fang，Liu Yunhe et al．MCSEM 3D modeling for arbitrarily anisotropic media．Chinese Journal of Geophysics，2014，57(12)：4110-4122.

[18] 王超．海洋可控源电磁法三维正演研究[学位论文]．北京：中国地质大学(北京)，2014.

Wang Chao．Marine CSEM 3D Forward Modeling [D]．China University of Geosciences (Beijing)，Beijing，2014.

[19] 韩波，胡祥云，Schultz A等．复杂场源形态的海洋可控源电磁三维正演．地球物理学报，2015，58(3)：1059-1071.

Han bo，Hu Xiangyun，Schultz A et al．Three dimensional forward modeling of marine controlled source electromagnetic field of complex source geometries．Chinese Journal of Geophysics，2015，58(3)：1059-1071.

[20] 汤文武，李耀国，柳建新等．基于二次电场的可控源电磁法三维矢量有限元正演模拟．石油物探，2015，54(6)：665-673.

Tang Wenwu，Li Yaoguo，Liu Jianxin et al．Three-dimensional controlled-source electromagnetic forward modeling by edge-based finite element using secondary electrical field．GPP，2015，54(6)：665-673.

[21] 严波，李予国，韩波等．任意方位电偶源的MCSEM电磁场三维正演．石油地球物理勘探，2017，52(4)：859-868.

Yanbo，Li Yuguo，Hanbo et al．The 3D forward mode-ling of the EM field for MCSEM using the electrical dipole source with arbitrary azimuth．OGP，2017，52(4)：859-868.

[22] 李建慧，胡祥云，陈斌等．复杂形态回线源激发电磁场的矢量有限元解．石油地球物理勘探，2017，52(6)：1324-1332.

Li Jianhui，Hu Xiangyun，Chen Bin et al．The vector finite element solution for the EM field induced by the loop source with complex form．OGP，2017，52(6)：1324-1332.

[23] 赵慧，刘颖，李予国．自适应有限元海洋大地电磁场二维正演模拟．石油地球物理勘探，2014，49(3)：578-585.

Zhao Hui，Liu Ying，Li Yuguo．The adaptive finite ele-ment forward modeling for 2D marine MT．OGP，2014，49(3)：578-585.

[24] 考夫曼 A A，凯勒 G V．频率域和时间域电磁测深．北京：地质出版社，1987.

Kaufman A A，Keller G V．Frequency and Transient Sounding．Geological Publishing House，Beijing，1987.

[25] Guptasarma D，Singh B．New digital linear filters for Hankel J0and J1transforms．Geophysical Prospecting，1997，45(5)：745-762.