三维广域电磁法自适应有限元正演

周峰，徐锦通，张志勇，蓝泽鸾，陈辉，汤洪志

（1.东华理工大学核资源与环境国家重点实验室，江西南昌 330013； 2.江西省防震减灾与工程地质灾害探测工程研究中心，江西南昌 330013； 3.中南大学有色金属成矿预测与地质环境监测教育部重点实验室，湖南长沙 410083）

0 引言

广域电磁法是由何继善院士及其团队潜心研究多年提出的一种频率域有源电磁方法［1-4］。相比可控源音频大地电磁法（CSAMT），广域电磁法仅需单分量电磁场即可定义视参数，而无需正交电场和磁场，野外工作效率较高［5-8］。另外，由于CSAMT 沿用了卡尼亚视电阻率计算公式，需要利用正交电场和磁场计算阻抗，导致过渡区和近区的数据不满足平面波的假设，因而这部分数据不能参与资料解译工作，这在一定程度上影响了CSAMT 勘探效果。相反，广域电磁法摒弃了远区的限制，扩展了频率域可控场源观测区范围，具有勘探深度大、观测范围广、工作效率高等优点，对地下地质结构的探测效果明显。近年来，广域电磁法已经在油气、页岩气等矿产资源、地下水资源以及环境工程等领域得到了广泛应用，取得了较好的应用效果［1，9-10］。

现阶段广域电磁法资料解译工作主要集中在一维、二维［6，11］，迫切需要开展三维广域电磁法资料解译技术研究。三维正演模拟是资料解译的重要环节，同时又是认识复杂地质结构响应特征的必要途径。因此，研究高效率、高精度的三维广域电磁法正演方法成为目前研究的热点。

三维数值模拟方法主要包含积分方程法［12-16］、有限差分法［17-19］、有限体积法［20-22］、谱元法［23-24］以及有限单元法［25-29］。每一种数值模拟方法都具有不同的应用特点，如：积分方程（IE）法是一种纯二次场算法，只需要对异常体进行剖分，具有较高的计算精度，但该方法形成的系统矩阵是密实矩阵，求解效率低，需要开展相应的加速算法（快速FFT）［14，30］；有限差分法数值模拟方法理论简单，网格拓扑关系简洁［17］，但刻画复杂地形及模型较困难；有限单元法［31-32］是目前应用最广的数值模拟方法之一，随着非结构网格离散技术的发展，对起伏地形和复杂地质结构模型的剖分更加精细［33］。随着自适应技术的发展，基于后验误差估计策略的自适应有限元技术得到飞速发展［34-35］。自适应加密技术主要采用两类策略：一类基于全局误差估计的自适应加密策略；另一类是面向目标的自适应加密策略。前者主要关注整体数值误差；后者重点考虑的则是观测点的数值误差，对非重点区域不做过多的考虑，避免了不必要的网格离散。当前，自适应有限元正演技术已在频率域电磁法数值模拟中得到较广泛的应用［36-38］。但是，关于广域电磁法的自适应正演算法鲜有报道。此外，基于有限体积法的电磁法数值模拟近年来得到了更多关注，尤其是与非结构网格相结合的应用发展较快。谱元法是有限元和谱方法相结合而发展起来的一种数值模拟方法，兼顾了谱方法的高精度和有限元的灵活性［23-24］，近年来在电磁法数值模拟方面得到了广泛关注。目前，该方法网格离散大多采用结构化网格，而关于非结构化的谱方法数值模拟技术报道较少［39-40］。

综上所述，为了更好地模拟复杂地电模型及提高广域电磁法正演精度，本文采用自适应有限元技术并结合非结构四面体单元实现三维广域电磁法的高精度、高效率正演计算。首先，从广域电磁法满足的微分方程出发，推导出三维广域电磁法满足的电场方程，并采用矢量有限元对电场方程进行离散，形成线性方程组；其次，开展以电流密度法向连续后验误差估计策略为依据的自适应有限元算法进行高精度电磁场分量计算，并采用一种稳健的Brent-Dekker算法计算广域视电阻率；最后，设计地电模型验证自适应算法的正确性，并对广域视电阻率响应特征进行分析。

1 方法理论

1.1 广域电磁法边值问题

假定负谐时间因子为e-iωt，电场E和磁场H满足Maxwell方程

式中：σ表示电导率；μ表示磁导率；σ-iωε为复电导率，其中ω为角频率，ε为介电常数；Js为外部电流密度。对式（1）两边取旋度，可得双旋度结构的电场方程

以式（5）为控制方程，利用Galerkin 有限元推导有限元方程。设余量

令r在每个计算单元Ωe内满足

式中Ne表示矢量基函数。将式（6）代入式（7），经矢量恒等变换，式（7）可改写为

采用非结构四面体单元离散整个求解区域，引入矢量基函数N，则四面体单元内任意一点的电场可由各条棱边的电场和矢量基函数表示

式中：K=C+M为刚度矩阵，其中；E为未知电场；为外部场源积分项。本文将场源置于四面体单元的各边［41］，假定源的中点坐标为(x0，y0，z0)，外部电流密度Js可表示为

式中：I为电流密度；δ为Delta函数。为了对长导线源进行积分，可将有限长源看作多个偶极源的叠加。此外，本文采用扩边技术加载截断边界条件，即E|Γ=0，这里Γ为无穷远边界。

1.2 基于后验误差估计的自适应正演策略

为保证有限元数值解的精度，需人为对网格单元进行优化，这一过程耗时且费力。近年来，基于后验误差估计的自适应有限元技术逐步应用于地球物理数值模拟，提高了正演计算的精度。为了实现三维多场源广域电磁法自适应正演模拟，需要构建三维广域电磁法问题对应的对偶问题

式中：T表示测点所涉及四面体单元个数；Vk为第k个单元的体积；I0=(1，1，1)为三分量的单位源；W表示电场对应的对偶场；κ表示权重系数，取值为观测点到场源中心之间的距离的立方，其作用是平衡场源对每个观测点的影响［33］。式（13）与式（5）具有相同的结构形式，采用Galerkin矢量有限元对式（13）进行离散，可得到相同的刚度矩阵K，即对偶问题形成的线性方程组为

后验误差分析是自适应网格加密的核心引擎［33-35］。为此，本文采用以电流密度法向连续为后验估计的网格自适应细化方案［33，36-37］，单元的后验误差e可根据下式求得

式中：∂Δi表示四面体的第i个面；J=σ·E表示电流密度，J1、J2分别表示共面∂Δi的面电流密度；n表示面法向。因此，本文利用式（15）评价单元的数值误差、标记需要细化的四面体单元。同时，为了提高测点计算精度，对对偶问题线性方程（式（14）） 进行求解，对偶问题所对应的单元误差w可根据下式求得

式中：W1、W2分别表示共面∂Δi的对偶场；σ1、σ2分别表示相邻四面体单元的电导率。然后，采用加权思想，最终获得加权误差估计为

通过上述过程可以标记每个单元的加权误差，然后按照一定比例对四面体进行细化处理，从而得到新的网格单元；然后，重复上述过程，直到有限元解的精度满足要求或达到最大网格细化迭代次数为止；最后，通过得到的电场或磁场分量计算广域视电阻率。需要指出，可控源电磁法在场源位置不满足电流密度法向连续的条件，因而包含场源的单位误差会很大。为此，本文利用局部加密技术细化场源，在计算单元误差时将源所在单元剔除，以降低源的影响，保证细化过程朝预设目标进行。

1.3 广域视电阻率的计算

广域电磁法的特点是采用电磁场单分量进行计算，利用均匀半空间条件下电磁场分量关于电阻率的非线性方程定义广域电磁视电阻率，克服了卡尼亚电阻率定义中对远区的要求，拓展了人工场源的观测范围，简化了野外观测装备，可显著提高工作效率。为了阐述广域视电阻率的优势，以x方向布设的电偶源为例，均匀半空间的视电阻率为

式中：r表示收发距；k表示波数；I表示电流强度；Ex是沿观测方向x两点M、N 间的电场；为M、N 间的电位差，其中表示观测电极M、N 之间的距离；，其中dL表示偶极源的长度；FE-Ex(ikr)=1-3 sin2ϕ+e-ikr(1+ikr)，其中ϕ表示测点和偶极源中心连线与偶极源之间的夹角，此项涉及电阻率、频率、场源与观测点等相关参数。为提高广域电磁视电阻率计算过程的稳定性，本文将Brent-Dekker求根算法［42-43］引入广域视电阻率的计算，并与二分法和直接迭代法［44］计算结果进行对比。

假设一个电阻率为100 Ω·m 的均匀半空间，在原点设置一个沿x方向的场源，观测点位于(0，5000 m，0)。文献［43］对基于Brent-Dekker 算法的求解过程做了简单描述，具体过程如下。

首先，假设在区间(a0，b0)内求解函数f(x)，若f(a0)与f(b0)符号相反，则说明在区间(a0，b0)存在根，否则无根。利用该算法进行迭代求解会涉及三个参数：ak-1，bk-1，bk-2，其中ak-1是第k-1 次迭代所得区域的下边界，bk-1表示第k-1次迭代所得到的迭代点，bk-2是上一次迭代的迭代点或预估值。反复迭代，直到f(ak-1)与f(bk-1)符号相反，以保障(ak-1，bk-1)内存在解，并满足|f(ak-1)|＜|f(bk-1)|。

然后，对每一步迭代进行逆二次插值

如果f(bk-1) 与f(bk-2)相等，则使用割线法

式中bk-2取值为ak-1。若要保证sk被接纳，需要满足两个条件：①sk在或内；②若上次迭代采用二分法计算，则必须满足，若上次迭代采用逆二次插值法或割线法，则必须满足，否则，不能保证求解过程的稳定性。表1所示为均匀半空间模型不同算法广域视电阻率计算结果。

表1 均匀半空间模型不同算法广域视电阻率计算结果

从表1 可知，Brent-Dekker 算法相较于直接迭代法，在全区迭代次数相对稳定、相同精度条件下，在过渡带（8 Hz）和近区（1 Hz）的迭代次数明显少于二分法。因此，为稳定求解三维模型视电阻率，本文采用Brent-Dekker算法。

2 算例分析

2.1 层状模型

为了验证算法的正确性，设计图1 所示的三层层状介质模型。场源沿x方向布设，长度为10 m，中心坐标为（0，0，0），发射电流为1 A，设置14个发射频率，范围为0.1～4096 Hz。沿x轴-6000～6000 m 范围布设了一条观测剖面，测点间距为200 m。求解区域为 ［-40 km，40 km］3。该模型的卡尼亚视电阻率和广域视电阻率剖面见图2，0.1 Hz和8 Hz所对应的视电阻率和相对误差error曲线见图3和图4。

图1 层状模型示意图

图2 层状模型卡尼亚视电阻率（ 上）和广域视电阻率（下）剖面

图3 层状模型0.1 Hz 时1-th、3-th 以及13-th 网格下的卡尼亚（上）和广域电磁（下）视电阻率（a）和相对误差（b）曲线

图4 层状模型8 Hz 时1-th、3-th 以及13-th 网格下的卡尼亚（上）和广域电磁（下）视电阻率（a）和相对误差（b）曲线

图3 和图4 分别展示了不同层次（1-th，3-th，13-th）网格下得到的卡尼亚和广域电磁视电阻率曲线及误差分布。由图可见，随着网格密度的增加，计算精度明显提升，相对误差在4%以内。13-th网格条件下得到视电阻率与解析解吻合最好，误差相对较小。还可以看出，卡尼亚视电阻率相对广域电磁视电阻率受非平面波的影响更大。图5为以电流密度法向连续的自适应有限技术细化三种层次网格（1-th、3-th、13-th）后的加权误差β，可见随着网格密度增加，误差整体出现降低趋势，符合预期。

图5 层状模型不同网格剖分下的加权误差β

2.2 梯形山模型

为了分析广域视电阻率受地形影响的特征，设计图6 所示的梯形山模型。地下半空间的电阻率为100 Ω·m，空气部分的电阻率为1×108Ω·m；求解区域为 ［-10 km，10 km］3； 场源沿y方向布设，场源中心坐标为（5 km，0，0），长度为100 m，发射电流为1 A。在地面沿y方向布设了一条观测剖面，观测范围为［-1500 m，1500 m］。分别计算0.1、2、4 Hz的电场、磁场、卡尼亚视电阻率以及广域视电阻率（图7和图8）。

图6 梯形山地形模型示意图

图7 梯形山地形模型电场分量Ey（上）和磁场分量Hx（下）的实部（a）和虚部（b）曲线

图8 梯形山地形模型卡尼亚和广域视电阻率曲线

图7 为不同频率下电场分量Ey和磁场分量Hx的实部和虚部，可见地形会引起电、磁场的变化，在地形起伏的拐点位置畸变最明显。该模型的正地形相当于一个低阻异常体，对电流有吸引作用，因而地形正上方的曲线向下弯曲。

图8 是不同频率卡尼亚视电阻率和广域视电阻率曲线，可见随着频率增高，卡尼亚视电阻率曲线出现抬升特征，数据进入了近区，无法真实客观地反映地下介质电阻率信息。然而，整个频率段的广域视电阻率值受场源影响较小，一定程度上能够避免卡尼亚视电阻率因场源影响而被抬升的情况。频率为0.1 Hz时此现象最明显，几乎在每一个观测点上这两种视电阻率的差值均达到一个数量级。与卡尼亚视电阻率曲线相比，广域视电阻率值更接近真实情况，说明了后者具有更大勘探深度，具备近区勘探能力。

2.3 球状异常体模型

设计图9 所示的球状异常体模型。计算区域为[-10 km，10 km]3。沿x方向布设线性源，源中心坐标为(0，-5 km，0)，长度为1 m，发射电流为1 A。在地面设置一个正方形观测区域：x=[-1000 m，1000 m]，y= [-1000 m，1000 m]。分别计算1、2、8、16 Hz 的电场分量Ex、磁场分量Hy、卡尼亚视电阻率以及Ex的广域视电阻率，结果见图10～图13。

图9 球状异常体模型示意图

图10 球状异常体模型1 Hz（a）和2 Hz（b）电场Ex（上）和磁场Hy（下）分布

图11 球状异常体模型1 Hz（a）和2 Hz（b）卡尼亚视电阻率（上）和广域视电阻率（下）分布

图12 球状异常体模型8 Hz（a）和16 Hz（b）电场分量Ex（上）和磁场分量Hy（下）分布

图13 球状异常体模型8 Hz（a）和16 Hz（b）卡尼亚视电阻率（上）和广域视电阻率（下）分布

由图10和图11可见，越靠近场源，电磁场越强，符合物理规律。磁场对地下异常体的敏感度弱于电场，因而电场Ex受低阻异常体的影响出现等值线弯曲现象，而磁场Hy则基本不受影响，即：低阻异常体能够吸引电流，而对磁场则没有影响。图中两种视电阻率等值线形态区别较大，主要是受场源的影响，卡尼亚视电阻率在靠近场源区域已经进入近区或过渡区，因而电阻率明显高于背景电阻率，其电阻率信息不能反映真实的地下电阻率情况；广域视电阻率则完全不同，在靠近场源的区域不会过早进入近区或过渡区，因而其视电阻率信息能够较真实地反映地下介质的地电信息。因此，广域电磁法能在一定程度上避免远区的限制，具有更大的勘探深度。

图12展示了球状异常体模型在8、16 Hz时的电场分量Ex、磁场分量Hy分布，图13是对应的卡尼亚视电阻率和广域视电阻率平面等值线图。从图12可以看出，电场和磁场平面分布特征与图7类似，不同之处在于随着频率增高，电场幅值增加，而磁场幅值减小。图13中的视电阻率等值线分布特征与图8类似，不同的是卡尼亚视电阻率受非平面波的影响变弱，仅较少的测点进入了近区，大部分测点的视电阻率接近真实值。因此，广域视电阻率能够很好地反映低阻异常体的存在，在异常体正上方呈现低阻特征。

总之，在同等收发距的前提下，广域电磁法摒弃了“远区”的限制，不受非平面波的影响，与CSAMT相比具有更大的勘探范围和深度，得到视电阻率能够较真实地反映地下介质的电性分布，适用性更强，准确度更高。

3 结论

本文采用基于电流密度法向连续为后验策略的自适应有限元算法，实现了三维广域电磁法正演计算。通过模型测试得出以下几点认识：

（1）针对广域视电阻率的计算，采用了一种稳健的Brent-Dekker 求根算法，该算法结合二分法、割线法等算法的优势，使求根过程更稳健，迭代次数相对适中，确保了计算结果的正确性。

（2）将电流密度法向连续的后验误差估计策略的自适应有限元算法引入三维广域电磁法正演模拟，提高了广域视电阻率的计算精度。通过与解析解对比，验证了算法的正确性。该算法测试发现，场源附近单元受场源影响较大，会造成场源附近单元的过度细化，影响了场源积分的求解效率。

（3）对球状异常体和梯形山地形等典型模型进行正演计算，分析了广域电磁法的探测深度与范围。在同等情况下，CSAMT 因受到非平面波的影响，视电阻率曲线发生抬升现象，不能很好地反映真实的地下介质信息，而广域视电阻率摒弃了近区的限制，得到的视电阻率更真实可靠，拓展了有源电磁法的观测范围，同时增加了勘探深度。模型计算结果表明，地形对广域电磁法和CSAMT 都有影响，尤其在地形拐点，畸变尤为明显。因此，在进行资料解释时需要考虑地形的影响，尽量降低地形影响带来的解释误差。