0 引言

大地电磁测深法是一种频率域电磁勘探方法，以天然交变电磁波为场源研究地下介质的地电特征，是深部地球物理探测的一种重要方法和必不可少的手段。高精度、高效率的正演模拟是反演的基础，因此对复杂地质模型进行高效正演是大地电磁测深数据反演的重点和难点。

经典的大地电磁正演求解方法主要包括积分方程法(IE)^[1-2]、有限差分法(FD)^[3]、有限元法(FE)^[4]、有限体积法(FV)^[5-6]等。这些方法大多采用结构网格进行区域离散，难以精确处理大地电磁正演中的起伏地形及倾斜界面等复杂结构。非结构化网格具有剖分灵活且剖分区域边界光滑的优点，可以很好地解决前述问题。Wannamaker等^[7]通过连接非均匀正交网格的对角线将网格分为四个等腰三角形，基于规则结构的三角网格实现了二维有限元正演。Franke等^[8]对三角网格法进行扩展，提出了一种基于自适应非结构化网格的电磁场二维有限元方法。非结构化网格可以实现更灵活的结构边界映射，因此地球物理学家开展了基于非结构化网格的二维和三维大地电磁正、反演模拟^[9-12]。最常用的非结构化网格由Delaunay三角形生成。Shewchuk^[13]通过对Delaunay三角形进行最小角度和最大面积约束，构造了适合区域划分的、稳健的、对图形契合度较高的三角形分割。本文利用Shewchuk开发的开源三角网格生成器Triangle生成非结构化三角形网格。

大地电磁有限元正演计算的精度在很大程度上取决于网格剖分的合理性及精细程度。因此，要获得较高的正演精度，就需要对网格尽可能地加密。但是，过于精细的网格会带来巨大的计算量，对于非常复杂的模型甚至无法实现。自适应网格加密技术可在精度与计算量之间达到最优平衡。目前该技术已广泛应用于大地电磁正演^[14-17]及其他地球物理正演问题^[18-23]。然而，该技术所需的加密比例及最大加密次数等参数往往需要调优，以获得精度与效率间的最佳平衡。另外，有限元模拟的精度在很大程度上依赖于自适应策略，而自适应模拟的程序实现相对比较繁琐。

除了自适应技术，加速大地电磁正演的另一途径是多重网格(MG)法。数学上已严格证明MG法至少对线性椭圆型偏微分方程是最优的^[24-25]，其计算量仅与节点数成正比，但是需要在粗网格与细网格之间来回迭代，计算过程相对复杂。Bornemann等^[26]提出了一种从粗网格到细网格的单向瀑布式多重网格(CMG) 算法，只需采用插值和迭代两种运算，无需对粗网格进行修正，程序实现简便。Chen等^[27]和胡宏伶等^[28]基于有限元误差展开式提出了外推瀑布式多重网格(EXCMG)算法，相较于经典CMG算法，该算法收敛更快，构造的有限元近似解更准确。EXCMG算法成功应用于2.5/3维直流电阻率法^[29-30]和三维大地电磁法^[31]的有限元正演模拟，加速效果明显。然而，这些算法均采用结构化的正交网格(二维长方形/三维长方体)，大大限制了该方法的应用。目前尚未有基于非结构网格的EXCMG方法的应用。究其原因主要有以下几点：其一，自适应技术相对成熟，学者们更倾向于选用该技术进行网格加密及正演计算，但自适应方法忽视了嵌套网格之间有限元解的联系，导致每次生成密网都只能从零初值进行迭代，这个过程中实际存在大量的重复计算；其二，要想针对非结构网格设计EXCMG算法，首先要解决的问题就是延拓算子的构造。已有的EXCMG算法都是基于结构化正交网格，很难将其直接推广到非结构化网格。

针对以上问题，本文提出利用EXCMG算法求解二维大地电磁正演问题。首先，利用节点有限元方法在非结构网格上离散Maxwell方程；然后，构造基于三角单元的外推及插值算子，利用粗网格上的有限元解构造密网上的有限元解的近似，并结合基于不完全Cholesky分解预条件BiCGStab作为磨光算子，对离散系统加速求解。在数值实验中，详细比较了不同网格尺寸下EXCMG算法、BiCGStab、聚合型代数多网格(AGMG)法的计算速度和迭代次数，并利用国际标准模型COMMEMI-2D1、COMMEMI-2D4验证了方法的正确性和有效性，并将该方法推广到起伏地形下各向异性介质的大地电磁正演问题。

1 控制方程

假设入射到地面的天然电磁场是平面电磁波，可忽略位移电流，取时谐因子为e^－iωt，则大地电磁场方程^[32]为

$ \nabla \times \boldsymbol{E}=i \omega \mu_0 \boldsymbol{H} $

(1)

$ \nabla \times \boldsymbol{H}=\boldsymbol{\sigma} \boldsymbol{E} $

(2)

式中：σ为电导率张量；μ₀为真空磁导率；ω为角频率；E为电场强度；H为磁场强度。

对于二维各向异性异常体，假设其构造走向沿x方向，即电导率只沿y轴和z轴变化，沿x轴是不变的。假设该异常体的电性主轴i、k、m分别与x、y、z轴所成夹角为α_i、α_k、α_m，电导率张量σ可表示为^{[9, 11]}

$ \boldsymbol{\sigma}=\boldsymbol{A}_{\boldsymbol{\sigma}_0} \boldsymbol{A}^{\mathrm{T}} $

(3)

式中：A表示坐标旋转矩阵

$ \boldsymbol{A}=\boldsymbol{R}_i \cdot \boldsymbol{R}_k \cdot \boldsymbol{R}_m $

(4)

其中R_i、R_k、R_m分别为电性主轴i、k、m的坐标旋转矩阵，其表达式分别为

$ \boldsymbol{R}_i=\left[\begin{array}{ccc} \cos \alpha_i & -\sin \alpha_i & 0 \\ \sin \alpha_i & \cos \alpha_i & 0 \\ 0 & 0 & 1 \end{array}\right] $

(5)

$ \boldsymbol{R}_k=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos \alpha_k & \sin \alpha_k \\ 0 & -\sin \alpha_k & \cos \alpha_k \end{array}\right] $

(6)

$ \boldsymbol{R}_m=\left[\begin{array}{ccc} \cos \alpha_m & -\sin \alpha_m & 0 \\ \sin \alpha_m & \cos \alpha_m & 0 \\ 0 & 0 & 1 \end{array}\right] $

(7)

σ₀为对角矩阵

$ \boldsymbol{\sigma}_0=\left[\begin{array}{lll} \sigma_i & & \\ & \sigma_k & \\ & & \sigma_m \end{array}\right] $

(8)

其中σ_i、σ_k、σ_m分别为电性主轴i、k、m方向上的电导率。

对于二维介质，式(1)、式(2)可进一步解耦简化为^[13]

$ \nabla \cdot(\underset{=}{\boldsymbol{\tau}} \nabla u)+\lambda u=0 $

(9)

式中：在TE极化模式下

$ \left\{\begin{array}{c} u=E_x \\ \lambda=\sigma_i \\ \underset{=}{\boldsymbol{\tau}}=\tau=\frac{1}{\text{i} \omega \mu_0} \end{array}\right. $

(10)

在TM极化模式下

$ \left\{\begin{array}{c} u=H_x \\ \lambda=\text{i} \omega \mu_0 \\ \underset{=}{\boldsymbol{\tau}} =\left[\begin{array}{ll} \tau_{11} & \tau_{12} \\ \tau_{21} & \tau_{22} \end{array}\right] \end{array}\right. $

(11)

其中：E_x、H_x分别表示电场和磁场的x分量；$\underset{=}{\boldsymbol{\tau}}$的各元素表达式分别为

$ \begin{aligned} \tau_{11} & =\frac{1}{\sigma_k} \sin ^2 \alpha+\frac{1}{\sigma_m} \cos ^2 \alpha \\ \tau_{22} & =\frac{1}{\sigma_k} \cos ^2 \alpha+\frac{1}{\sigma_m} \sin ^2 \alpha \\ \tau_{21} & =\tau_{12}=\frac{1}{2}\left(\frac{1}{\sigma_k}-\frac{1}{\sigma_m}\right) \sin 2 \alpha \end{aligned} $

当α=0°时，式(9)退化为文献[32]中的对角各向异性的情况；进一步地，若σ_i=σ_k=σ_m，式(9)进一步退化为各向同性的情况。

2 有限元分析

对求解区域Ω离散后，式(9)可转化为等价的弱形式，即寻找u_h∈U_h，使得

$ \int\limits_{\varOmega}\left(\underset{=}{\boldsymbol{\tau}} \nabla u_{\mathrm{h}} \cdot \nabla v_{\mathrm{h}}+\lambda u_{\mathrm{h}} v_{\mathrm{h}}\right) \mathrm{d} x=0 \quad \forall v_{\mathrm{h}} \in \boldsymbol{V}_{\mathrm{h}} $

(12)

式中：V_h为检验函数空间；U_h为试探函数空间。

假设空间V_h的基函数为ϕ_j(j=1, ⋯, M)，这里M表示空间V_h的维度，则

$ u_{\mathrm{h}}(x)=\sum\limits_{j=1}^M U_j \phi_j(x) \quad \forall u_{\mathrm{h}} \in \boldsymbol{V}_{\mathrm{h}} $

(13)

将式(13)代入式(12)，取v_h=ϕ_i(1≤i≤M)，得

$ \sum\limits_{j=1}^M \int\limits_{\Omega}\left(\underset{=}{\boldsymbol{\tau}} \nabla \phi_j \cdot \nabla \phi_i+\lambda \phi_j \phi_i\right) \mathrm{d} x U_j=0 $

(14)

将求解区域Ω分解为n_w个三角单元，单元w编号记为1, 2, ⋯, n_w。将式(14)按照三角单元进行分解，得到

$ \sum\limits_{w=1}^{n_w} \sum\limits_{j=1}^M \int\limits_w\left(\underset{=}{\boldsymbol{\tau}}\nabla \phi_j \cdot \nabla \phi_i+\lambda \phi_j \phi_i\right) \mathrm{d} x U_j=0 $

(15)

将式(15)转为矩阵形式，可得有限元方程

$ \boldsymbol{K} \boldsymbol{U}=0 $

(16)

其中

$ \left\{\begin{array}{l} \boldsymbol{K}=\left(\begin{array}{ll} \boldsymbol{K}_{11} & \boldsymbol{K}_{12} \\ \boldsymbol{K}_{21} & \boldsymbol{K}_{22} \end{array}\right) \\ \boldsymbol{U}=\left(\begin{array}{l} \boldsymbol{E}_x \\ \boldsymbol{H}_x \end{array}\right) \end{array}\right. $

(17)

式中K₁₁和K₂₂是复对称矩阵，且K₁₂=K₂₁^T，因此K为稀疏复对称矩阵。U包含了所有节点上电场和磁场。求解区域Ω的边界Γ存在非零Dirichlet边界条件，因此

$ \boldsymbol{K} \boldsymbol{U}_{\varOmega / \varGamma}=-\boldsymbol{K} \boldsymbol{U}_{\varGamma} $

(18)

式中：Ω/Γ为求解区域Ω边界Γ的内部区域；U_Ω/Γ为区域Ω/Γ内所有节点的电场和磁场；U_Γ为边界Γ上所有节点的电场和磁场。

最终所需求解的线性方程组为

$ \boldsymbol{K} \boldsymbol{U}=\boldsymbol{b} $

(19)

$ \boldsymbol{b}=-\boldsymbol{K} \boldsymbol{U}_{\varGamma} $

(20)

根据式(19)求得E_x、E_y、H_x、H_y后，可得TE模式和TM模式下的视电阻率为

$ \left\{\begin{array}{l} \rho_{\mathrm{TE}}=\frac{1}{\omega \mu}\left|\frac{E_x}{H_y}\right|^2 \\ \rho_{\mathrm{TM}}=\frac{1}{\omega \mu}\left|\frac{E_y}{H_x}\right|^2 \end{array}\right. $

(21)

相位为

$ \left\{\begin{array}{l} \phi_{\mathrm{TE}}=\arctan \frac{E_x}{H_y} \\ \phi_{\mathrm{TM}}=\arctan \frac{E_y}{H_x} \end{array}\right. $

(22)

式中μ表示磁导率。

3 半结构外推瀑布式多重网格法

有限元经离散可得到大型稀疏、复线性方程组(式(19))，若网格剖分太细，直接求解会存在耗时较长且存储需求过高的问题。迭代法在内存方面较有优势，EXCMG算法利用粗网上的近似解，通过外推和高次插值技术构造密网有限元解的高阶逼近，作为迭代初值，可加速收敛。均匀正交网格上的EXCMG算法的详细介绍参见文献[29-31]。本文将EXCMG算法推广到非结构化三角网格。多网格延拓算子的具体构造过程如下。

对图 1a所示三角形进行逐层剖分，网格由粗到细分别记为Z₀(图 1a)、Z₁(图 1b)、Z₂(图 1c), 并对节点进行编号。图 1a中最粗网格Z₀的节点1~3处的有限元解记为U⁽⁰⁾={U₁⁽⁰⁾, U₂⁽⁰⁾, U₃⁽⁰⁾}，图 1b中网格Z₁的节点1~6处的有限元解记为U⁽¹⁾={U₁⁽¹⁾, U₂⁽¹⁾, U₃⁽¹⁾, U₄⁽¹⁾, U₅⁽¹⁾, U₆⁽¹⁾}，以此类推。

下面介绍构造网格Z₂上有限元解U⁽²⁾的高阶逼近$ \tilde{\boldsymbol{U}}^{(2)}$，以此作为迭代初值，加速收敛。

假设椭圆边值问题的线性有限元解U_j⁽ⁱ⁾有如下误差渐近展开式^[29]

$ \begin{aligned} e^{(i)}\left(x_j, z_j\right) & =U_j^{(i)}-u\left(x_j, z_j\right) \\ & =A\left(x_j, z_j\right) h_i^2+O\left(h_i^4\right) \end{aligned} $

(23)

式中: $\left(x_j, z_j\right)$为编号$j=1, \cdots, 15$的点对应的坐标; $h_i$为第$i$层网格的最大直径, 这里$i=0, 1, 2 ; e^{(i)}$为第$i$层网格的误差; $O\left(h_i^4\right)$表示$h_i^4$的同阶无穷小量; $A\left(x_j, z_j\right)$是一个与$h$无关的平滑函数。

利用U⁽⁰⁾和U⁽¹⁾的线性组合逼近下一层网格Z₂上点1、2、3的有限元解

$ C U_j^{(1)}+(1-C) U_j^{(0)}=U_j^{(2)}+O\left(h_0^4\right) \quad j=1, 2, 3 $

(24)

式中C是待求常数。将式(24)代入式(23)，得到C=1.25。因此，对于三角形顶点1、2和3，可利用如下外推公式得到网格Z₂上的有限元解U_j⁽²⁾的高阶逼近

$ \tilde{U}_j^{(2)}=U_j^{(1)}+\frac{1}{4}\left[U_j^{(1)}-U_j^{(0)}\right]=U_j^{(2)}+O\left(h_0^4\right) $

(25)

利用式(23)可得网格Z₂中点4、5、6的有限元近似解

$ \begin{aligned} \widetilde{U}_j^{(2)} & =U_j^{(1)}+\frac{1}{8}\left\{\left[U_i^{(1)}-U_i^{(0)}\right]+\left[U_k^{(1)}-U_k^{(0)}\right]\right\} \\ & =U_j^{(2)}+O\left(h_0^4\right) \\ j & =4, 5, 6 ; i, k=1, 2, 3 \end{aligned} $

(26)

由于已经求得$\tilde{U}_i^{(2)} $(i=1, 2, ⋯，6)这6个有限元近似解，可以对其进行二元二次完全多项式插值，构造其余9个点的有限元近似解

$ \widetilde{U}_j^{(2)}=\sum\limits_{i=1}^6 N_{i, j} \tilde{U}_i^{(2)} \quad j=7, 8, \cdots, 15 $

(27)

其中形函数N_{i, j}(i=1, 2, ⋯，6, j=7, 8, ⋯，15)为编号7~15的点所在三角形单元上的Lagrange型形状函数，求解过程如下。

对于任意一个如图 1a所示的三角形单元，三个顶点的坐标为(x_i, y_i)，其面积记为Δ₁₂₃，有

$ \begin{aligned} \varDelta_{123}= & 0.5 \times\left|\begin{array}{ccc} 1 & 1 & 1 \\ x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \end{array}\right| \\ = & 0.5 \times\left[\left(x_2 y_3-x_3 y_2\right)+\left(x_1 y_3-x_3 y_1\right)+\right. \\ & \left.\left(x_1 y_2-x_2 y_1\right)\right] \\ = & 0.5 \times\left[\left(y_2-y_3\right)\left(x_1-x_3\right)-\right. \\ & \left.\left(y_3-y_1\right)\left(x_3-x_2\right)\right] \end{aligned} $

(28)

记a_i=y_j－y_k，b_i=x_k－x_j，c_i=x_jy_k－x_ky_j，其中i、j、k按照1、2、3轮换。

取三角形中任意一点P(x, y)，将点P与三角形三个顶点相连，得到三个小三角形单元P₂₃、P₃₁、P₁₂，面积分别用Δ₁、Δ₂、Δ₃表示。令

$ L_i=\frac{\mathit{\Delta}i}{\mathit{\Delta}_{123}} \quad i=1, 2, 3 $

(29)

式中L_i为点P的面积坐标。面积坐标与直角坐标之间的转换关系为

$ L_i=\frac{1}{2 \mathit{\Delta}_{123}}\left(a_i x+b_i y+c_i\right) \quad i=1, 2, 3 $

(30)

利用式(30)可得任意三角单元上的二元二次Lagrange型形函数为

$ \left\{\begin{array}{l} N_1=\left(2 L_1-1\right) L_1 \\ N_2=\left(2 L_2-1\right) L_2 \\ N_3=\left(2 L_3-1\right) L_3 \\ N_4=2 L_2 L_3 \\ N_5=2 L_3 L_1 \\ N_6=2 L_1 L_2 \end{array}\right. $

(31)

式(27)中的形函数如式(31)的形式。将求得的$\tilde{\boldsymbol{U}}^{(2)} $作为BiCGStab的迭代初值，结合不完全Cho-lesky分解预条件，可计算得到最细网格上的近似解。下文将三角形网格节点个数称为自由度，即每层网格上未知数的个数。EXCMG算法具体描述如下。

(1) 求解最粗两层网格Z₀、Z₁上的有限元解U⁽⁰⁾、U⁽¹⁾。

(2) 对i=1, ⋯，Y(Y表示网格剖分次数)执行：

a. 利用外推公式(式(25)、式(26))和插值公式(式(27))构造密网Z_i上的初值$\tilde{U}^{(i)} $；

b. 以U⁽ⁱ⁾为初值，利用BiCGStab迭代求得近似解U⁽ⁱ⁾。

(3) 输出近似解U^(Y)，退出程序。

4 模型计算与分析

首先，对国际标准各向同性模型COMMEMI-2D1和COMMEMI-2D4进行模拟，比较EXCMG与BiCGStab、AGMG等算法的计算结果，验证EXCMG算法的精度。然后，对比EXCMG法的计算结果与文献[33]中的数值解，分析EXCMG算法对有限元正演过程的加速效果。最后，通过起伏地形和各向异性模型，分析EXCMG算法的精度及适应性。本文各个算例初始非结构化网格均利用国际常用的三角网格剖分器Triangle生成，剖分参数均为－DpqA，最小角度约束为20°。

4.1 国际标准模型COMMEMI-2D1

COMMEMI-2D1模型如图 2所示。求解区域为200 km(y)×200 km(z)的正方形，z方向空气和地下区域各占一半。对模型利用二分加密生成第2~6层网格。初始单元数为614，自由度为1259；第2层网格加密生成9824个三角形单元，自由度为4973；⋯⋯；第6层加密网格包含628736个三角形单元，自由度为1258433。从图 3可以看出，异常体附近区域经5次加密后，网格已经十分细密。

从初始网格剖分到最密网格剖分，网格数和自由度都是成倍增长，如果直接求解式(19)，那么存储量与计算量都是巨大的。但是，若采用EXCMG算法，则只需求解初始网格和第1层网格上的有限元解，并基于这两组解进行外推和线性插值，便可求得下一层网格的有限元近似解，计算量和存储需求都大大降低。

对于模型COMMEMI-2D1，对频率f=0.1 Hz利用EXCMG方法计算TE和TM极化模式下的视电阻率和相位，并计算与文献[33]中标准解的绝对误差error，结果见图 4和图 5。可以看出，TE模式下，利用本文计算方法计算的视电阻率和相位与文献[33]中标准解的最大error分别为0.0019 Ω·m、0.0025°；TM极化模式下，视电阻率和相位的最大error分别为0.0263 Ω·m、0.0054°。由此可见，EXCMG算法求解正演模型有较高的精度。

为了验证EXCMG算法对有限元方法正演求解的加速效果，对模型COMMEMI-2D1，针对频率f=0.1 Hz分别采用EXCMG、AGMG、SSOR-BiCGStab和ILU-BiCGStab等算法对有限元构造的复线性方程组(19)进行求解，统计TE和TM模式下的计算时间和迭代次数，见表 1。这四种算法均设定最大迭代次数为5×10⁵，相对误差限为1×10^－8。网格规模超过39296时，SSOR-BiCGStab算法的计算时间过长，加速效果不明显，因此这里不对其进行统计。

从表 1可以看出，在TE和TM模式下，随着网格规模的增大，EXCMG算法的迭代次数逐渐减小，而SSOR-BiCGStab和ILU-BiCGStab算法的迭代次数逐步增加，AGMG算法的迭代次数先增后减。因此，EXCMG算法所需计算量明显低于其他算法。

从表 1还可以看出，AGMG算法的计算时间明显小于SSOR-BiCGStab、ILU-BiCGStab和AGMG算法。TE和TM模式下AGMG法与EXCMG法的加速比见表 2，最大加速比可达42.31。

单元数为9824时，四种算法的归一化相对残差收敛曲线见图 6。可见不论TE模式还是TM模式下，EXCMG算法均比AGMG、SSOR-BiCGStab和ILU-BiCGStab算法相对残差收敛更快，说明EXCMG算法对大规模网格模型的计算优势明显。

针对三角形单元数为39296的情况，对比频率为0.1、1、10Hz时EXCMG与AGMG的计算时间和迭代次数，结果见表 3。由表可见，EXCMG算法的计算时间和迭代次数均小于AGMG，说明EXCMG算法对于大规模网格模型的计算优势明显。

4.2 国际标准模型COMMEMI-2D4

图 7所示国际标准模型COMMEMI-2D4为一复杂2D地电模型，其异常体个数、形状和分布均比模型COMMEMI-2D1更复杂。选择该模型验证EXCMG算法对复杂地电模型求解的适应性。求解区域为：y∈[－100 km, 100 km]，z∈[－100 km, 100 km]。测点位于地面、沿y方向[－100 km，100 km]范围内均匀布设，点距为100 m，共2001个测点。计算频率范围为1×10^－4~1×10^0.5 Hz，对数等间隔地分布，共50个频点。正演模拟的初始三角形网格单元数为4011，经两次加密后单元数为64176。异常体附近区域网格剖分情况见图 8。

基于EXCMG算法的TE模式和TM模式下的视电阻率及相位见图 9。由图可见：横向上，TE模式和TM模式的视电阻率剖面上出现两个明显的高阻异常区域；纵向上，在视电阻率剖面上存在高阻异常区域，相位也呈现明显的分层特征，说明计算频率及异常体的形状大小、分布情况、电阻率均会对视电阻率和相位产生影响。

为了检验EXCMG算法对不同频率大地电磁的正演精度，在固定网格单元数为64176的情况下，分别对0.01、1、9 Hz进行电磁场正演模拟，EXCMG计算结果与解析解^[33]最大绝对误差见表 4。视电阻率、相位最大绝对误差分别小于0.001 Ω·m和0.001°，可见EXCMG算法计算精度较高。

另外，在相同网格规模(64176)、相同频率(0.01 Hz)下，TE模式下EXCMG算法所需时间为16.41 s，直接求解和SSOR-BiCGStab则分别耗时57.23 s、60.72 s；TM模式下EXCMG算法所需时间为19.76 s，直接求解和SSOR-BICGStab耗时分别为59.16 s、62.31 s。可见使用EXCMG加速有限元方法正演模拟复杂地电模型，在保证精度的同时可明显减低计算时间。

4.3 各向同性山脊模型

上文两个标准模型中无地形起伏，为了验证EXCMG算法对带地形模型的适应性，建立图 10所示各向同性山脊模型^[34-35]。求解区域为：y∈[－40 km, 40 km]，z∈[－20 km, 20 km]。模型的空气区域中有一个电阻率为100 Ω·m的电阻率各向同性山脊，高度为2000 m，宽度为37944.8 m，山脊的形状函数为

$ z(y)=-\frac{h}{2}\left(\cos \frac{2 \pi y}{B_x}+1\right) $

(32)

式中: h=2000 m，表示山脊高度；B_x=37944.8 m，表示三角函数的波长。

在f=0.1 Hz的情况下，针对图 10模型的EXCMG算法的初始三角形网格单元数为893，第5层的网格单元数为228608。EXCMG算法和文献[34]中使用有限元直接求解所得的视电阻率及二者间的error见图 11。由图可见，TE模式下，视电阻率和相位的最大error分别为0.0048 Ω·m和0.009°；TM模式下，视电阻率和相位的最大error分别为0.0011 Ω·m和0.00064°，可见EXCMG算法的求解精度较高。TE模式下，EXCMG算法耗时6.82 s，直接求解耗时9.92 s^[34]；TM模式下，EXCMG算法耗时5.71 s，直接求解耗时9.83 s^[34]。因此，EXCMG算法求解速度快，计算精度高。

4.4 简单各向异性模型

为了进一步验证基于EXCMG方法的有限元求解算法对大地电磁各向异性模型正演计算的适应性，建立图 12所示二维简单各向异性模型。求解区域为：y∈[－10 km, 10 km]，z∈[－10 km, 10 km]。

该简单各向异性模型的初始三角形网格单元数为80，经5次均匀加密后的网格单元数为81920。计算频率为0.1 Hz。基于EXCMG算法，TE模式的计算时间为2.67 s，迭代次数为291；TM模式的计算时间为1.38 s，迭代次数为217。同样的参数情况下，AGMG算法TE模式的计算时间为8.46 s，迭代次数为448；TM模式的计算时间为2.17 s，迭代次数为240。因此，与AGMG算法相比，EXCMG算法更省时，收敛更快。

文献[11]利用有限元法模拟了该模型的电磁场，所得到的线性有限元方程组用预条件共轭梯度法求解。图 13为频率0.1 Hz时EXCMG算法和文献[11]中该模型时的视电阻率及error曲线。由图可见，TE、TM模式下的视电阻率最大error分别为1.11×10^－2、3.18×10^－2 Ω·m。该模型测试结果进一步验证了EXCMG算法的精度，同时也说明EXCMG算法对电性各向异性模型的模拟具有较强的适应性。

4.5 复杂各向异性模型

建立图 14所示包含两个各向异性异常体的复杂模型，分析本文方法对复杂模型的计算精度和适应性。测点沿y=－50~50 km均匀布设，点距为50 m，总测点数为2001；模拟频率范围为0.001~10 Hz，以对数等间隔分布，共50个频点。

该模型初始三角形单元数为7709，经过3次均匀加密后的网格单元数为123344。采用EXCMG算法计算耗时4.56 s，用时较短。图 15为TE、TM模式下采用EXCMG算法得到的视电阻率和相位的拟断面。由图可见，TE模式和TM模式的视电阻率和相位剖面上都出现不同程度的分层现象，且在视电阻率剖面上出现低阻区域。说明TE模式的视电阻率和相位剖面及TM模式的视电阻率剖面对异常体均有所反映。

4.6 带地形的各向异性模型

对图 10所示山脊模型和图 12所示各向异性模型进行组合，构建图 16所示的带地形的各向异性模型。求解区域为：y∈[－40 km, 40 km]，z∈[－20 km, 20 km]。测点沿y轴－39~39 km均匀布设，点距为250 m，共313个测点。计算频率范围为0.001~10 Hz，以对数等间隔分布，共50个频点。

模拟空间的初始三角形网格单元数为943，经过3次均匀加密，得到第4层的网格单元数为60352。异常体附近的网格剖分见图 17。基于EXCMG方法对该模型的正演模拟时间为3.53 s，耗时较短。

图 18是基于EXCMG算法的正演视电阻率及相位拟断面。由图可见，在异常体附近(y=－4~4 km)，TE、TM模式下的视电阻率剖面上均出现了低阻异常区，TE模式的相位剖面上出现极值，TM模式的相位剖面上存在明显的分层现象。

从图 18中提取f=0.1 Hz的视电阻率和相位数据，见图 19。与不带地形的各向异性模型的视电阻率(图 13)相比，由于山脊的存在，山脊范围内的视电阻率略高；与无异常体的山脊模型的视电阻率(图 11)相比，由于异常体的存在，异常体范围内视电阻率异常更大。还可以看出，受地形和各向异性异常体的双重影响，TE模式下的相位曲线在山脊范围内较TM模式变化更明显。

4.7 实际数据应用

基于贵州某山区数据^[36]建立图 20所示的带地形复杂二维各向异性地电模型，电性主轴k、m与坐标轴y、z的夹角α均为0°。求解区域为：y∈[0, 1400 m]，z∈[－900 m, 200 m]，测线及测点分布见图 21，测线长度为1400 m，点距为40 m。令测线方向为y轴，东北方向为正。

基于EXCMG算法的初始网格单元数为119，经过5次均匀加密得到第6层网格的单元数为121856。图 22为经5次加密得到的网格，可见求解区域剖分网格已较密集，故采用该加密网络进行计算。

对图 20所示模型基于EXCMG算法得到的视电阻率和相位拟断面见图 23。由图可见：TE模式下，视电阻率断面上频率为10 Hz、y=1200 m附近出现一个高阻异常区域；TM模式下，视电阻率断面上1~10 Hz、y=1200 m附近出现一个高阻异常区域。这些电性异常是地形起伏和电性异常体共同作用的结果。对比图 23和文献[36]，可见二者关于视电阻率和相位的分布特征基本一致，证明了EXCMG算法在计算含有起伏地形和复杂异常体的大地电磁正演问题中适用性较强，计算精度较高。

5 结束语

本文提出了利用基于非结构化三角形网格的外推瀑布式多重网格算法(EXCMG)快速求解电阻率各向异性介质中二维大地电磁有限元正演的大规模复线性系统。对本文的计算结果与SSOR-BiCGStab、ILU-BiCGStab、AGMG算法计算结果进行对比，可见前者可以使用更少的有限元解构造大规模网格节点上的有限元近似解，且随着网格规模的增加，迭代次数逐步降低，大幅度加快计算速度。

通过对国际标准模型、各向同性模型、各向异性模型及实际地电模型进行计算，验证了EXCMG算法计算结果的可靠性，证明了该算法计算精度高，可满足大规模复线性方程组的求解需求，是一种求解二维大地电磁复杂模型正演问题的高效方法。