昌彦君，肖明顺，孟永良，罗延钟

昌彦君，肖明顺，孟永良，罗延钟_大人地球物理科学

昌彦君¹，肖明顺²，孟永良³，罗延钟¹

（1.中国地质大学地球物理与空间信息学院，武汉，430074；2.中南冶金地质总局中南地质勘察院，武汉，430081；3.中国地质大学信息工程学院，武汉，430074）

摘　要：介绍了带地形的瞬变电磁2.5维有限元数值模拟的意义及相关内容研究的历程，建立了其解题思路和数值计算方法。首先，通过拉氏变换和傅氏变换将三维空间中时间域的四维矢量偏微分方程边值问题，转换为拉氏变换和傅氏变换域中的二维标量偏微分方程边值问题；其次，将待求解的二维标量偏微分方程边值问题变换为二维变分问题，进而用二维有限单元法求其数值解；最后，通过逆拉氏变换和反傅氏变换求得三维空间中电场的瞬变过程。该算法中，采用在自适应有限元网格中增加两对角线形成三角网格剖分，同时用高斯消元法消除矩形网格中心结点的待求未知量，既可较准确地模拟带地形起伏的任意二维复杂地电断面，又可节省计算量。采用直接计算二次场的算法，只需计算电场和磁场沿地电构造走向两个分量的一次场，不但计算精度较高，而且并不显著增加计算量。对几例典型地电模型进行了正演计算，并计算出其视电阻率曲线或拟断面图，获得了较好的效果，验证了所建立的算法的有效性。

关键词：瞬变电磁，2.5维，有限单元法，地形，自适应网格

前　言

瞬变电磁法作为时间域电磁法的一种，在金属矿和石油资源勘探、地下水、环境和工程勘查中有着广泛的应用。但目前的解释水平依然很低，主要局限于一维解释。现有的三维正反演算法计算量都很大，在没有很好地解决三维问题的速度与精度之前，它是不具备实用价值的。纯二维的正反演已有较多的研究成果，而且计算量不是很大，但由于场源的三维性质，纯二维的正反演只是粗略的近似。基于上述原因，处在二维和三维问题中间的2.5维问题得到了人们的重视。所谓2.5维，即二维模型三维场源，二维模型即三维地质体电（磁）参数在走向方向上无变化。二维模型是三维模型的特例，在许多特定的地质构造环境下，进行模型的二维近似处理是有效的和有意义的。此外，为了适用于模拟山区地形，必须考虑地形问题。因此，研究带地形的瞬变电磁2.5维有限元数值模拟问题是务实的选择。

2.5维时间域电磁场数值模拟问题，是目前国际上未妥善解决的计算地球物理的疑难问题。相关问题的研究始于20世纪70年代中期，Coggon^［1］开创性地将有限单元算法引入到电磁场数值模拟问题中。80年代，Wannamaker^［2，3］利用有限元单元法实现了带地形的二维大地电磁数值模拟。90年代开始为人们所重视，Leepin^［4］用有限差分法（FDM）求解了2.5维矩形电流回线的时间域电磁响应问题。Torres-Verdin等^［5］及Habashy等^［6］基于非均匀介质中电磁散射问题的局部非线性近似，用积分方程法研究了频率域中垂直磁偶极源的2.5维正演；积分方程法涉及到求张量格林函数和散射电流的褶积问题，只适用于模拟为数不多的局部异常体，否则计算量会成倍增加。Unsworth和Everett等^［7，8］对洋中脊模型计算过电流偶极源的电磁场瞬变响应，他们用有限单元法求解总电场和总磁场沿构造走向分量的两个联立线性二阶偏微分方程组，为消除场源处的奇异性，采用了从求解区中“挖除”场源（TM模式）的做法，实现了频率域电流偶极源电磁场的2.5维有限单元法模拟；但当异常较小时，其计算精度欠佳。Yu^［9］在计算洋中脊模型的电偶极源瞬变电磁响应时采用了二次场算法，避免了上述问题；而且它可用于围岩横向非各向同性的情况，这对海洋电磁法十分重要；在进行频率域向时间域的转换时，通过采用Gaver-Stehfest变换作逆拉氏变换，将复数计算化为纯实数计算，提高了计算速度^［9，10］。进入21世纪后，Mitsuhata^［11］撇开常规的二次场算法，定义了一种新的拟δ函数模拟电偶极子源，并利用等参单元有限元模拟地形，实现了带地形的频率域水平电偶极子源2.5维有限单元正演，并利用余弦变换得到时间域电磁场值。Kerry等^［12］针对电磁法有限元数值模拟中矩形网格的不足，以二维大地电磁为例，提出了利用自适应有限元网格的数值模拟思路。近些年，国外研究热点转移到海洋电磁法正反演数值模拟研究，Abubakar等^［13］研究了海洋可控源电磁法（CSEM）2.5维正反演问题。在国内，中国地质大学（武汉）最早在20世纪80年代便开始了瞬变电磁2.5维电磁场的数值模拟研究^{［14，15］}，在频率域和时间域2.5维电磁法数值模拟问题上做了大量的研究工作。如：二维地电构造上谐变电偶极子和电偶源CSAMT法的2.5维有限单元算法^{［16～18］}；90年代，针对时间域电磁方法，开展了时间域和频率域电偶源2.5维有限单元正演模拟^［17］以及时间谱电阻率法的2.5维正演算法的研究^{［18～22］}。在此基础上，21世纪初又取得了一系列的进展^{［23～27］}。2001年，实现了二次场算法的中心回线瞬变电磁2.5维正演^{［25，27］}，并开发了正演软件以实现人机交互联作反演；后来，又进一步完善了二次场算法的中心回线瞬变电磁2.5维正演算法^{［28，29］}，并基于McGillivray辅助场原理^［30］，以磁偶源瞬变电磁场为例，推导了2.5维瞬变电磁反演的关键内容——偏导数矩阵^［31］。2008年开展了带地形的瞬变电磁2.5维有限元数值模拟研究^［32］，并取得了较好的计算效果。

瞬变电磁法有多种工作装置形式，可进行多种电磁场分量的观测。本文的主要内容是带地形的中心回线瞬变电磁2.5维的正演模拟，拟为建立实用近区人工源2.5维电磁法算法奠定基础。

图1　地电断面模型图k₁.三维异常体内的波数；k₂.均匀大地中的波数

1　电磁场的偏微分方程定解问题

三维各向同性均匀导电介质中，电（磁）偶极子形成的

电磁场满足Maxwell方程。在时间域中，有如下方程组：

在介质分界面上（如图1所示），满足衔接条件：

式中：ε为介电常数；σ为电导率；μ为磁导率；n为分界面的法向方向矢量；E为电场矢量；H为磁场矢量；J_e为场源电流密度矢量；J_m为场源磁流密度矢量。以下标1和2表示不同的介质。介质电性分界面为Γ，外法线方向为n，逆时针转90°为切线方向τ。

不失一般性，令电（磁）偶极子源位于笛卡尔坐标原点，则有

式中：δ（x）、δ（y）和δ（z）为狄拉克源函数；AB为电偶极距矢量，方向由负极指向正极。P_e＝AB×I为源的电流偶极距，S为供电回线面积矢量，按右手螺旋法则指向回线法线方向。P_m＝S×I为源的磁流偶极距。μ₀为真空中的磁导率。f_e（t）和f_m（t）分别为描述电性源和磁性源随时间变化的波形函数，在时间域（瞬变）电磁场研究中，通常取场源波形函数为上阶跃函数。

电磁场定解问题（1）～（4）是一个空间三维（x，y，z）加上时间维（t）的四维矢量问题，求其数值解的计算量非常巨大。

为降维（消除时间变量t），利用如下拉普拉斯变换对：

对给定的拉氏变量s值，上述定解问题是三维问题。

在介质电性参数（ε，σ和μ）与坐标y无关的二维地电条件下，利用如下傅氏变换对：

沿y轴对电磁场作傅氏变换

可进一步降维（消除y坐标）。利用矢量运算规则，由（9）和（10）式出发，可导出拉氏和傅氏变换（s，m）域中电磁场分量满足的偏微分方程：

式中：e_x、e_y、e_z、h_x、h_y、h_z和j_e，x、j_e，y、j_e，z、j_m，x、j_m，y、j_m，z分别表示（s，m）域中电场、磁场和场源电流密度、场源磁流密度沿笛卡尔坐标轴的分量；i为虚数单位。经过消元，并引用符号

可导出仅与（s，m）域中电场和磁场沿地电构造走向（y）分量有关的u和υ的偏微分方程组：并且，由式（3）和式（4）出发，借助于式（18）～（23）可以导出，在与构造走向垂直的xoz坐标平面上，无场源的介质分界线Γ两侧，u和υ满足的边值条件：

式中：肩标（1）和（2）分别标识分界线Γ两侧的参数；n和τ分别表示分界线Γ的外法线和切线方向矢量。于是，电磁场的计算问题归结为求（s，m）域电磁场y分量定解问题式（26）～（31）的数值解。在算出u和v后，可借助于由（18）～（23）式导出的关系式，计算（s，m）域中的其他电磁场分量：

应该指出，对于给定的s和m值，定解问题（26）～（31）式是二维标量问题。它的数值解法相对于求解原先的四维矢量问题，大为简化。计算思路是：①分别对若干给定的拉氏和傅氏变换变量s和m值，在（s，m）域中求电磁场分量u、v及e_x、e_z、h_x和h_z的数值解；②由（s，m）域中的电磁场分量，按（8）式作逆拉氏变换，求得傅氏变换（m）域中电磁场分量的瞬变响应（时间域解）；③由第2步算得的各给定时间（t）的m域电磁场分量，按（15）式作反傅氏变换，便可算出三维空间域中电磁场分量的瞬变响应。

2　（s，m）域电磁场分量的泛函极值

用有限单元法求定解问题式（26）～（31）的数值解。基于边值条件（28）～（31）式，可以导出在有限单元法算域（网格）外边界G满足第三类边值条件：

其中，系数η^u和η^v可由均匀半空间条件下（一次场）的电磁场分量近似算得：

可以证明，与偏微分方程边值问题式（26）～（31）和式（36）、（37）等价的泛函极值问题是：

式中：u；j_，、j_，、j_，和j_，、j_，、j_，分别表示场源电_exeyezmxmymz流密度和磁流密度沿x、y、z轴的分量；Ω表示求解区（xoz平面上的二维网格）；G表示求解区外边界；j_n和j_τ分别表示在G上场源电流密度沿外法线n和切线方向τ的分量。

3　有限单元法求解

用有限单元法求解变分问题（41）式的具体步骤包括区域剖分、线性插值、单元分析、总体合成和求取变分，将泛函极值问题离散化为大型稀疏线性方程组：

式中：［A］为P×P阶方阵（刚度矩阵），｛X｝为求解区内各矩形网格节点上电磁场分量u和v（待确定的）P阶未知列矢量，右端项｛S｝为各节点场源分布的P阶场源列矢量。我们采用了自适应有限元网格，即根据研究区域和研究目标体的形状自动划分有限元网格^［12］。这种网格剖分方案既适用于模拟物性界面，又有利于减少网格单元，且网格自动生成。四边形网格剖分为4个三角形后，网格中心节点的未知量u和v，经高斯消元处理，未进入线性方程组（41）中。这样大大降低方程阶数，减少计算时间消耗（约减少到八分之一）。图2为带地电模型的剖分网格；图3为地形网格剖分。

图2　带地电模型的剖分网格

图3　地形剖分网格

4　二次场算法

求解线性方程组式（41）算出的未知电磁场x（u，v）称为总场，它是同样场源在地下电性分布均匀时形成的一次场x₁与地下不均匀体引起的二次场x₂之和。仿照（41）式，可以写出一次场x₁的线性方程组

式中：［A₁］是在地下电性与场源处电性相同的均匀大地条件下的刚度矩阵；右端（场源）列矢量｛S｝与式（41）相同，这是因为预先假定场源相同。

由式（41）和式（42），容易得到：

式中：｛X₁｝和｛X₂｝＝｛X｝-｛X₁｝为一次场和二次场列矢量。

式（43）是直接计算二次场的线性方程组，其右端项｛S₂｝可由一次场和总场刚度矩阵之差［A₁］-［A］和一次场｛X₁｝之乘积算得。前者只在电性与均匀围岩不相同的不均匀体内及其相邻节点上才不为零，并且易于由［A］算出；后者（一次场）可由显式表达式通过数值计算获得，不过，当不均匀体范围较大时计算量相当大。

在按式（43）算出二次场x₂（u₂，v₂）后，可利用与总场公式（32）～（35）相似的公式，计算其余电磁场分量的二次场：

式中：电导率异常Δσ＝σ-σ₁；磁导率异常Δμ＝μ-μ₁；电磁场分量的下标1和2分别表示一次场和二次场。

在做该项目的二维正演计算时，我们都是采用直接计算二次场的算法。这是因为与直接计算总场相比，直接计算二次场通常可获得较高的计算精度；不过，其需要计算不均匀体内的一次场。

5　空间电磁场分量瞬变响应的计算

三维空间的电磁问题，经拉普拉斯变换、傅立叶变换之后，由原先的四维边值问题转而成为求解二维的边值问题。用前面所介绍的算法对一系列给定的拉氏变量s（也称频率）和傅氏变量m（也称波数），算出的拉氏和傅氏变换（s，m）域中的电（磁）场分量二次场后，可按式（8）和式（15）作逆拉氏变换和反傅氏变换获得空间电（磁）场分量的瞬变响应。

需要指出，波数选取和反傅氏变换算法对空间电（磁）场分量的计算精度和速度的影响很大。所以，要求选用适当的数值积分算法和波数，以保证在合理的计算精度下，尽量减少计算量。选择何种数值积分算法及波数的个数和分布，决定于傅氏变换函数随波数变化的形态。前人已做了很多关于最优选取波数的研究^{［16，26，28］}。我们采用根据波数渐近线划定波数范围，在每一十进位间以等自然对数间隔取定3～4个波数的方法^［4，28］。

关于逆拉氏变换，目前已有多种多样的算法^［19］，我们采用Gaver-Stehfest概率变换算法（简称G-S变换）^［10］。

设拉氏空间的像函数为F（s），实空间的像原函数为f（t），则

这里

式中：N是决定于计算机位数的正偶数。

或许G-S变换并不是最精确或最好的算法，但它是纯实数运算，而且只需对较少的s值做计算（我们对每一瞬变时间选用12个s值），因而是一种计算速度较快的算法。特别是我们在用G-S变换法做逆拉氏变换时，利用拉氏变换延迟定理在倍增的时间间隔中插值，从而实现了对密集采样时间瞬变过程的快速计算^［25］。

6　算　例

6.1　三层H模型

模型参数：电阻率ρ1＝100Ω·m、ρ2＝10Ω·m、ρ3＝100Ω·m，厚度H1＝100m、H2＝150m。图4绘制磁感电动势二次场和总场对比曲线。考虑到在以层状模型验证有限元算法的正确性和有效性时，视电阻率曲线可以更加直观地体现算法的精度，并且实际资料中常用视电阻率来进行分析，因此，我们计算了其全区视电阻率。图5绘制了视电阻率对比曲线，2.5维和一维计算的视电阻率曲线吻合得很好。

图4　H型断面感应电动势二次场（a）和总场（b）对比曲线

图5　H型视电阻率对比曲线

6.2　二维倾斜低阻体模型正演

下面4幅图中，图6为二维局部低阻体地电模型参数图：发射回线面积为1 600m²；接收回线等效面积1m²；发射电流1A。图7为计算获得的磁感电动势多测道剖面图。图8为由磁感电动势计算出的定义全区视电阻率拟断面图。图9为不同测点处磁场和磁感电动势定义视电阻率曲线。视电阻率拟断面图较好地说明了倾斜低阻体的存在。

图6　局部低阻体模型图

6.3　地形模型1

图10绘制了地形模型1示意图，图11为磁感电动势定义全区视电阻率拟断面图。同平坦大地局部低阻体情况类似的是，在拟断面图上出现了高阻异常区域和低阻异常区域。高阻异常区域为磁感电动势定义全区视电阻率绕曲所致；低阻异常区域反映了地下的低阻体。需要注意的是，拟断面图在150m深度左右两侧的向地表的绕曲可能为地形影响的结果。

图7　磁感电动势多测道剖面图

图8　磁感电动势定义全区视电阻率拟断面图

图9　不同测点处磁感电动势定义视电阻率曲线

（a）时间-视电阻率；（b）视深度-视电阻率

图10　地形模型1

图11　磁感电动势定义全区视电阻率拟断面图

6.4　地形模型2

图12绘制了地形模型2示意图，图13为磁感电动势定义全区视电阻率拟断面图。瞬变电磁对高阻的分辨率没有低阻好。同地形模型1类似的是：拟断面图在150m深度左右两侧向地表绕曲。

图12　地形模型2

图13　磁感电动势定义全区视电阻率拟断面图

7　结　语

关于带地形的瞬变电磁2.5维有限元数值模拟，在文中建立了其基本解题思路和数值计算方法：①通过拉氏变换和傅氏变换将三维空间中时间域的四维矢量偏微分方程边值问题，转换为拉氏变换和傅氏变换域中的二维标量偏微分方程边值问题；②将待求解的二维标量偏微分方程边值问题变换为二维变分问题，进而用二维有限单元法求其数值解；③通过逆拉氏变换和反傅氏变换求得三维空间中电场的瞬变过程。该算法有以下特点：①采用二维有限单元算法，在自适应有限元网格中增加两对角线形成三角网格剖分，同时用高斯消元法消除矩形网格中心结点的待求未知量，这样，既可较准确地模拟带地形起伏的任意二维复杂地电断面；又可节省计算量；②采用直接计算二次场的算法，只需计算电场和磁场沿地电构造走向两个分量的一次场，因而，不但计算精度较高，而且并不显著增加计算量；③采用Gaver-Stehfest变换法作逆拉氏变换，并利用拉氏变换延迟定理在倍增的时间间隔中插值的方法，可实现对密集采样时间瞬变过程的快速计算。

参考文献

［1］Coggon J H.Electromagnetic and electrical modeling by the finite-element method［J］.Geophysics，1971，36（1）：132-155.

［2］Wannamaker P E，Stodt J A，Roij L.Two-dimensional topographic responses in magnetotellurics modeled using finite elements［J］.Geophysics，1986，51（11）：2 131-2 144.

［3］Wannamaker P E，Stodt J A，Roij L.A stable finite element solution for two-dimensional magnetotel-luric modeling［J］.Geophys.J.R.Astr.Sic，1987，88：277-296.

［4］Leppin M.Electromagnetic modeling of 3-D sources over 2-D inhomogeneities in the time domain［J］.Geophysics，1992，57（8）：994-1 003.

［5］Torres-Verdin C，Habashy T M.Rapid 2.5dimensional forward modeling and inversion via a new nonlinear scattering approximation［J］.Radio Science，1994，29（4）：1 051-1 079.

［6］Habashy T M，Groom R W，Spies B R.Beyond the Born and Rytov approximations：A nonlinear approach to electromagnetic scattering［J］.J.Geophys.Res.，1993，98（B2）：1 759-1 775.

［7］Unsworth M J，Traves B J，Chave A D.Electromagnetic induction by a finite electric dipole source o-ver a 2-D Earth［J］.Geophysics，1993，58：198-214.

［8］Everett M E，Edwards R N.Transient marine electro magnetics：the 2.5-D forward problem［J］.Geo-phys.J.Int.，1993，113：545-561.

［9］Yu L M.Computation of the electrical responses of mid-ocean ridged structures［Ph.D.thesis］［M］.Toronto：University of Toronto，1994.

［10］Knight J H，Raiche A P.Transient electromagnetic calculation using the Gaver-Stehfest inverse La-place transform method［J］.Geophysics，1982，47（1）：47-50.

［11］Mitsuhata Y.2-D electromagnetic modeling by finite-element method with a dipole source and topography［J］.Geophysics，2000，65（2）：465-475.

［12］Kerry Key，Chester Weiss.Adaptive finite-element modeling using unstructured grids：the 2D magnetotelluric example［J］.Geophysics，2006，71（6）：291-299.

［13］Abubakar A，Habashy T，Druskin V，et al.Two-and-half-Dimensional forward and inverse model-ing for marine CSEM problems［C］.SEG New Orleans 2006Annual Meeting，USA，2006：750-754.

［14］张胜业.磁偶极子场中二维导体电磁响应的数值解［D］.北京：武汉地质学院北京研究生部硕士论文，1981.

［15］罗延钟，孟永良.关于用有限单元法对二维构造作电阻率法模拟的几个问题［J］.地球物理学报，1986，29：613-621.

［16］罗延钟，张桂青.电子计算机在电法勘探中的应用［M］.武汉：武汉地质学院出版社，1987.

［17］曾繁京.二维地电构造上水平电偶源电磁场的有限单元算法［D］.武汉：中国地质大学硕士论文，1988.

［18］孟永良，罗延钟.二维地电构造上谐变电偶极子电磁场的有限单元算法［J］.中国地球物理学会年刊，北京：地震出版社，1994.

［19］昌彦君，张桂青.电磁场从频率域转换到时间域的几种算法比较［J］.物探化探计算技术，1995，17：25-29.

［20］孟永良，罗延钟.电偶源CSAMT法二维正演的有限单元算法［C］.罗延钟等主编：勘查地球物理勘查地球化学文集，北京：地质出版社，1996（第20集）：103-114.

［21］罗延钟，昌彦君，张桂青.层状可极化大地上偶极—偶极装置瞬变电场的异常性态［C］.叶笃正主编：赵九章纪念文集，北京：科学出版社，1997.

［22］罗延钟，孟永良，昌彦君.时间谱电阻率法的二维正演算法［C］.首届中国与北美地区华人勘探地球物理学术讨论会（论文摘要集），美国，休斯顿：1998，20-22.

［23］昌彦君.时间谱电阻率法研究［D］.武汉：中国地质大学博士论文，2000.

［24］孟永良，罗延钟，昌彦君.时间谱电阻率法的二维正演算法［J］.地球科学：中国地质大学学报，2000，25（6）：656-662.

［25］罗延钟，昌彦君.G-S变换快速算法［J］.地球物理学报，2000，43（5）：684-690.

［26］王华军.中心回线瞬变电磁2.5维正反演方法研究［D］.武汉：中国地质大学博士论文，2001.

［27］王华军，罗延钟.中心回线瞬变电磁2.5维有限元算法［J］.地球物理学报，2003，46（6）：855-862.

［28］熊彬.大回线瞬变电磁2.5维正演算法及其资料解释方法中若干关键技术研究［D］.武汉：中国地质大学博士论文，2004.

［29］熊彬，罗延钟.电导率分块均匀的瞬变电磁2.5维有限元数值模型［J］.地球物理学进展，2004，19（3）：616-620.

［30］McGillivray P R，Oldenburg D W，Ellis R G，et al.Calculation of sensitivities for the frequency-do-main electromagnetic problem［J］.Geophys.J.Int.，1994，116：1-4.

［31］熊彬，罗延钟，强建科.瞬变电磁2.5维反演中灵敏度矩阵计算方法［J］.地球物理学报，2006，49（2）：590-597.

［32］肖明顺.带地形的瞬变电磁2.5维有限元数值模拟研究［D］.武汉：中国地质大学硕士论文，2008.