0 引言

大地电磁三维正演模拟是理解大地电磁物理现象及三维地质体电磁响应规律的有效手段。目前，大地电磁二维正、反演已经比较成熟，而大地电磁三维正、反演还不完善，对于复杂地质体的响应规律还需要进一步认识。因此，高效高精度的三维正演成为研究热点和难点。近年来，三维正演数值方法发展迅速^[1-5]。主流的数值方法有积分方程法、有限差分法、有限单元法和有限体积法。积分方程法由于只需对异常体进行剖分、不需考虑截断边界、可以快速模拟简单的地质异常体、计算精度高，在模拟简单地电模型时具有其他方法难以比拟的优势^[2]；但是，对复杂的大尺寸异常体模型以及起伏地形进行模拟时，张量格林函数的计算变得相当复杂。传统的有限差分法虽然理论简单，易于实现，但在处理界面突变处的场间断问题时较为困难。相对于传统的有限差分法，交错有限差分法和有限体积方法既允许由于电性差异引起的电场法向分量不连续性现象存在，也保证了电场和磁场的通量守恒，从而得到广泛应用。Xiong等^[3]将交错差分技术引入到地球物理电磁正演中。Smith^[4]研究了频率域的交错网格差分实现方法，并采用预条件技术和散度校正技术使得方程收敛速度明显提升。谭捍东等^[5]、陈辉等^[6]将交错差分算法应用到三维大地电磁中，论述了实现的细节，并检验了算法的正确性和计算精度。但是交错有限差分法要求网格剖分十分规则，遇到复杂形态异常体或者起伏地形时，难以做到高精度模拟。

有限元法具有内边界条件自动满足、网格剖分灵活等优点，适合模拟复杂形态异常体。但是三维大地电磁中电磁场的各个分量互相耦合，使得有限元方程中总体系数矩阵的阶数远远大于二维，方程组求解困难。直到近十年，三维大地电磁有限元算法才得到快速发展。为了解决内部边界处电磁场的突变问题，Badea等^[7]采用涡流场中常用的磁矢位和电标位即A－Φ系统，引入库伦规范，避免了场的不连续性，提高了算法的稳定性。Um等^[8]直接运用洛伦兹规范对矢量位的三分量进行解耦，使得双旋度方程化为赫姆霍兹方程，大大提高了计算效率。Mitsuhata等^[9]采用电矢位和磁标位即T－Ω系统，并利用六面体单元及矢量基函数和标量基函数分别对电矢位和磁标位进行离散，并验证了算法的精确性。徐志峰等^[10]同样利用磁矢位和电标位进行了可控源三维有限元正演，并研究了矢量位和电标位的意义。利用势函数求解大地电磁虽然可以避免场的不连续性问题，可是转化为电场或者磁场需要求场的导数，将导致精度损失。王若等^[11]尝试了单分量三维可控源有限元模拟研究，提高了三维有限元的计算速度。张继锋^[12]使用散度校正的标量有限元模拟了可控源的电磁响应。徐凌华等^[13]、童孝忠等^[14]从广义变分原理出发进行了大地电磁三维标量有限元正演。但是，传统的标量有限元法在直接求解电场或磁场时，存在散度校正和难以处理不连续边界条件的缺点^[15]。

矢量有限元由于其基函数的无散性和允许场法向分量的突变，在直接计算电场时具有较高的精度。Nam等^[16]采用不规则六面体矢量有限元直接计算电场，研究了起伏地形下大地电磁电阻率和相位的变化规律。王烨^[17]进行了大地电磁矢量有限元理论的总结。刘长生^[18]基于非结构化剖分，采用网格局部加密技术对大地电磁进行了矢量有限元正演模拟，并取得了较好的结果。顾观文等^[19]采用矢量有限元法研究了地形起伏条件下三维阻抗张量的变化规律。

在前人基础上，本文推导了矢量有限元算法的离散形式，并对强加边界条件的方法进行了改进，采用SSOR(symmetric successive over relaxation)预处理技术的双共轭稳定梯度法求解方程组，以提高计算效率。为了验证算法的有效性，计算了国际模型COMMEMI 3D-1的大地电磁响应，研究了典型三维低阻模型的视电阻率及相位分布规律及响应特征。

1 矢量有限元算法 1.1 三维大地电磁泛函表达式

设时谐因子为e^－iωt，频率域电磁场的麦克斯韦方程组为(在大地电磁频率范围内，忽略位移电流)

式中：E为电场强度；ω为角频率；μ为磁导率；H为磁场强度；σ为电导率；J为电流密度；ρ为电荷密度；D为电位移矢量。介质的本构方程为

式中：B为磁感应强度；ε为介电常数。

对式(1)两边取旋度，将式(2)代入，得到频率域的双旋度方程：

式(3)是电磁场的微分控制方程，在二维情况下可以做简化，转化为单分量赫姆霍兹方程。但是，在计算三维模型时，由于在异常体界面产生场的突变，导致电场散度在边界处无法计算；因此，式(3)不能表示为单分量赫姆霍兹方程，这使得三维运算量远大于二维。对式(3)运用广义变分原理^[16]可建立大地电磁方程的泛函：

式中，V表示三维数值计算的区域。由于采用线性有限元插值，式(4)中双旋度项为0，需利用第一矢量格林定理降低微分阶数，式(4)化为

式中：S表示外边界面；(∇×E)·(∇×E)对应磁场能量；－iωμσE²对应电场能量；∫_S(E×∇×E)·dS代表能量的流动，值得说明的是，电磁场中场分量的内边界条件与外边界条件均体现在这一项中。式(5)的整个泛函表达式是电磁场的能量总和，而广义变分法的原理就是求取使得能量泛函最小的电磁场分布。

1.2 矢量有限元离散表达式

将计算区域(图 1a)离散为六面体单元；在每个单元内，将自由度放在棱边上，采用Whtiney矢量基函数N_j^[15]近似电场，棱边编码如图 1b所示。单元内的电场表示为

式中，E_j^e为单元e中第j条棱边的电场值。

将式(6)代入式(5)，由于式(5)中∫_S(E×∇×E)·dS代表电磁能量沿垂直交界面方向流动，而且电磁场内部区域无场源，因此在单元界面两侧无能量奇点，对内部区域的任一面流入和流出的能量相等，其面积分互相抵消，只剩下区域外边界；在外边界处设置合适的边界条件，即可求得电场分布。式(5)泛函化为

式中，n_e为单元数。

为了求式(7)泛函的极小值，令式(7)中单元泛函在第i条棱边E_i^e(i=1,…,12)处的偏导数为0，即：

式中，单元系数矩阵K^e=∫_Ω[∇×N_i·∇×N_j－iωμσN_i·N_j]dV。令式(8)=0，将单元棱边的局部编码替换成全局编码，同时组装所有的n_e个单元的K^e，并强加边界条件，得到有限元的线性方程组：

式中：K为大型稀疏复对称矩阵，$K=\sum\limits_{e=1}^{{{n}_{e}}}{K{{~}^{e}}}$；E=(E₁,E₂，…，E_np)，n_p为所有棱边总数；b由强加边界条件决定。

1.3 边界条件与加载方法

在未加入边界条件之前，式(9)中总体系数矩阵K为奇异矩阵，方程(9)无唯一解。为了确定唯一的电磁场，需在计算区域V的外边界指定电场或者磁场的切向分量，即加载合适的边界条件。

利用混合积定理将式(5)中的面积分表示为∫_S(E^e·n_S×∇×E^e)dS/2，其中，n_S为单元对应面上的法向量。本文在四个侧面及顶面上采用狄利克雷边界条件，则上述面积分为常数，其对未知棱边电场的导数为0，在总体系数矩阵中涉及边界棱边的元素并没有影响。大地电磁场中，入射场只有水平分量，垂直分量是由异常体电磁感应产生的二次场。假设异常体距计算外边界足够远，异常体所产生的二次场未影响到边界，即计算区域的侧面边界处E_z为0。由此，侧面及顶面边界条件如下。

顶面边界，电场设为

z、y侧边界，边界处的场值可直接设为背景场：

同理，在 x、z侧面，边界条件为：

式中，E_x0为层状介质的背景电场，由一维有限元正演计算得出。

为了方便计算层状介质，在底面边界对E_x采用一阶边界条件，则式(5)中的面积分需要计算；底面边界处的边界条件化为

式中，e_x为x方向的单元向量。

式(10)—式(12)为狄利克雷边界条件，式(13)为一阶吸收边界条件。将式(13)代入式(5)中的面积分得

同理，通过将式(5)代入式(14)，对面积分项进行离散，可得到边界项的单元矩阵。

假设大地电磁场的场源一般位于高空^[20]，不在计算区域中，导致外边界条件的加载十分重要；外边界条件同时也代表入射场在截断边界处的分布，其加载对结果的精度影响较大。数值实验表明，在侧面边界的垂直分量强加为0非常必要。

有限元法中一般采用乘大数法^[21]强加狄利克雷边界条件，这种加载方法简便但是会有精度损失。同时，乘以大数，总体系数矩阵内部数量级差别变大，条件数大大增加，增加了迭代法求解难度。乘大数法强加之后，当罚系数过大时，总体系数矩阵数量级差别较大，导致结果精度不高甚至出现错误结果；罚系数太小时，会影响边界条件的强加效果。因此，本文采用直接法^[22]强加边界条件，可以减少矩阵阶数，处理过程直观，且适合用迭代法求解以保证精度，不过需要对未知棱边重新编码。其原理如下。

将总体系数矩阵分块表示：

式中：a为未知棱边(外边界棱边)的编号；b为已知棱边的编号；E_b为外部边界上的已知电场值；E_a为待求的未知电场值。在计算泛函导数时，由于外边界处电场由狄利克雷边界条件确定，故不需求泛函在外边界上棱边电场的导数，即K_ba、K_bb无用。则式(15)可化为

因此，只需计算一次矩阵和向量的乘积就可以得到式(9)右端项b=－K_abE_b。

1.4 大型稀疏病态线性方程组求解

由于有限元离散后得到的总体系数矩阵性态较差而且阶数很大，在求解方程组时需采用预处理技术进行迭代求解，以提高收敛速度。二维情况下通常采用ILU(0)预条件技术，可取得很好的效果；但在三维情况下，随着网格剖分的不均匀和频率的降低，ILU(0)预条件的收敛效果会明显变差。故本文采用较为稳定的SSOR预处理技术^[23]。

为了提高收敛速度，采用一维线性有限元的数值解作为迭代解法的初始解，当异常体和围岩的电阻率差别不大时，可以有效提高收敛速度。当计算层状介质时，迭代次数为0，这验证了算法计算一维层状介质的正确性。

常用的解有限元方程组迭代方法有：复共轭梯度平方法(CGS)、双共轭稳定梯度法(BICGSTAB)、广义最小残量方法(GMRES)、双共轭梯度法(BICG)和准最小残差法(QMR)。在SSOR预处理下，通过采用Matlab求解函数，对比了针对同一模型三维大地电磁矢量有限元形成的系数线性方程组，不同求解方法的收敛时间见表 1(x,y,z三个方向网格数目分别为27，27，31)。由表 1可知，高频收敛速度快，低频收敛速度慢。这种现象被称为低感应数问题^[24]。Smith^[4]认为，频率变低时，式(3)中的电导率所起的作用变小，造成矩阵性态变差，进而导致低频时，方程收敛速度变慢。

通过对比不同求解方法可以看到：双共轭稳定梯度法(BICGSTAB)在高频情况下计算速度较快，而复共轭梯度平方法(CGS)在计算低频时消耗时间较少，其他求解方法在求解时间上均较长。但是，CGS方法收敛残差跳跃大，不稳定。本文采用收敛平滑的BICGSTAB方法。矢量有限元形成的系数矩阵的性态受网格剖分影响较大：当网格尺寸较为接近时，收敛较快；当网格尺寸差异较大时，收敛速度明显变慢。

2 三维大地电磁矢量有限元正演算法正确性验证 2.1 国际模型对比

为了检验算法及程序的正确性，对国际模型^[25]COMMEMI3D-1A进行数值模拟，并与文献^[25]的参考结果进行了对比。x、y、z方向网格数分别为33,31,29(5层空气层)，计算区域为16 600 m×20 200 m×24 300 m，空气层厚度为3 550 m，空气层电阻率设为10⁸ Ω·m，容许误差设为10^-6。模型示意图如图 2所示。

计算频率为10.0 Hz和0.1 Hz，测线位置为y=0，计算结果如图 3所示。由图 3可知，10.0 Hz时，本文算法得到的ρ_xy和ρ_yx都较COMMEMI参考结果小。这种现象在文献^[10]中曾出现。这是因为COMMEMI参考结果是许多算法的平均值，计算结果差异较大。本文算法的结果也在国际模型参考结果的误差之内。0.1 Hz时，不同算法的参考结果较为一致，本文算法和COMMEMI的结果吻合较好。

为了研究矢量有限元法在模拟分界面处电场法向分量的突变效果，在地表上(-300 m,-600 m)处，分析了矢量有限元和标量有限元E_z分量沿着z轴的结果(图 4)。由模型示意图(图 2)可知，在地表及其下部存在三个电性突变界面，位置分别为z=0 m(地-空界面)，z=250 m(低阻模型的上界面)，z=2 250 m(低阻模型的下界面)。从图 4可以看出：标量有限元和矢量有限元的E_z分量结果总体相似，尤其在低阻体中，两者的E_z分量都接近于0；但是矢量有限元中E_z在分界面处明显存在突变，而标量有限元中E_z在分界面处是平缓过渡。依据电流密度的法向分量连续可以推出，理论上E_z在界面两端是突变的。由此可以看出，矢量有限元可以比较精确地模拟法向分量突变这一物理现象，而标量有限元由于节点上只赋予一个自由度，不可避免法向分量的连续。

2.2 交错差分算法对比

为了进一步验证算法的可靠性，我们与交错差分算法^[6]进行对比，其模型示意图及结果见图 5。由图 5b、c看出，无论是视电阻率还是相位，矢量有限元的结果在高、中、低频率上均与交错差分算法的结果吻合较好，这进一步验证了本文算法的正确性。

3 三维异常体视电阻率和相位响应分析

为研究三维情况下异常体模型的视电阻率和相位断面图的响应特征，本文将图 6所示低阻体模型的视电阻率与相位和二维模式做了细致的对比。具体模型参数为：背景介质电阻率为100 Ω·m，异常体电阻率为10 Ω·m，尺寸为700 m×700 m×700 m，埋深400 m。正演采用网格数：31×31×33 (x，y，z方向)，计算频率：1~1 000 Hz，计算区域为: 21 100 m×21 100 m×15 250 m。在地表上沿着x方向放置三条测线，分别位于y= 0 m，y=200 m，y=400 m。

根据大地电磁探测深度计算公式^[26]δ=1 186 $\sqrt{\rho /f}$，计算出1 Hz下背景场的理论探测深度约为11 860 m。在低频时上述计算区域中，二次场为0的边界条件依旧可以近似满足。10 Hz时，地表张量阻抗的分布如图 7所示。

由图 7看出，由于模型对称，Re(Z_xy)与Re(Z_yx)也是关于x=y对称的，Re(Z_xx)与Re(Z_yy)的形态也基本一致。由此看出，当长方体异常体长宽比为1∶1或者形状对称时，Re(Z_xy)与Re(Z_yx)也是互相对称的。说明在三维情况下，已经不存在单一的TE和TM模式。由图 7看出，Re(Z_xx)与Re(Z_yy)分为四瓣，且较Re(Z_xy)与Re(Z_xy)小；可以认为Re(Z_xx)和Re(Z_yy)反映了边界累积电荷的影响，Re(Z_xx)与Re(Z_yy)四瓣的中心反映了异常体的边界。根据张量阻抗理论可知，当构造为二维构造时，Z_xx和Z_yy为0。即当异常体走向方向越长时，Z_xy与Z_yx越小，Z_xy与Z_yx的差异也越大。由此，从张量阻抗对异常体的形态进行初步估计，为精确反演提供依据。

图 8为不同测线的视电阻率和相位断面图。由于测线沿着x方向，所以三维情况下中心测线上的XY模式(E_x入射，磁场平行于走向)对应二维的TM模式(图 8中TM)，YX模式对应二维的TE模式(图 8中TE)。由中心测线上的视电阻率断面图可以看出，YX模式与XY模式视电阻率和相位形态均与二维TM模式的视电阻率形态相似，ρ_yx在低频部分更宽缓，其变化幅度较ρ_xy小，ρ_xy在低频部分更尖锐。同时随着测线远离中心测线，视电阻率和相位异常均逐渐变弱。由图 8中TM模式与三维XY模式的响应看出，在中心测线上XY模式与TM模式视电阻率响应基本一致；而YX模式在高频时结果和二维TE模式比较一致，低频时，与TE模式差异较大。这是由于在二维情况下，TE模式的电流方向和异常体走向平行，在边界上没有累积电荷产生；而在三维情况下，电流必然会垂直穿过异常体分界面，在分界面上产生累积电荷，进而产生异常电场。高频时，异常体有一定埋深，电流集中于浅层，导致产生的累积电荷和异常电场也很小，这时以电流之间的感应为主，因此中心测线上的XY模式和YX模式的电阻率响应均与二维较一致；低频时，电场衰减慢，会在异常体分界面处产生较强的累积电荷，而TM模式本身已经较大，因此YX模式和TE模式的视电阻率和相位差别较大。考虑电荷的影响，所以中心测线处的的XY模式和TM模式电阻率响应基本一致。

4 结论

1) 本文采用直接法强加边界条件，避免了传统乘大数法增加系数矩阵条件数的缺点，既保证了总体系数矩阵的对称性又保证了求解精度，降低了矩阵阶数，同时指出四个侧面的电场z分量需置为零。

2) 针对同一模型同一网格，对比了不同解方程迭代方法的求解时间，通过采用SSOR预处理的BICGSTAB求解方程组，同时以背景场值作为迭代的初始解，提高了求解速度。

3) 对比国际标准模型，并与国际参考结果进行对比，同时与标量有限元的结果进行了对比，论证了矢量有限元模拟分界面处电场法向分量突变的优势。

4) 分析了典型模型的波阻抗响应形态特征，Z_xx与Z_yy的大小可以在一定程度上反应异常体的简单特征。分析了异常体的视电阻率和相位的响应规律，由于边界累积电荷的影响，XY模式和YX模式的视电阻率响应均与二维的TM模式接近，低频情况下YX模式相位及视电阻率与TE模式差异较大。