0 引言

可控源电磁法(Controlled Source Electromagnetic Method，CSEM)具有抗干扰能力强、工作效率高、测量精度高等诸多优点，是研究中浅部电性结构的主要地球物理方法之一，广泛应用于固体矿产、水文、石油天然气普查、地热田勘探、环境地质调查及环境与工程地球物理勘查等领域^[1-3]。目前，实际勘探任务大多以复杂三维地质构造为主，传统的一维、二维假设条件下的资料解释技术不能满足目前的地质解释要求，迫切需要开发三维可控源电磁资料解释技术。正演是反演的基础，因此国内外学者都致力于任意复杂地形的地电结构CSEM三维正演模拟研究，寻求快速、高精度的三维正演方法。

针对CSEM三维正演问题，学者们做了大量研究工作。首先针对场源奇异性导致正演精度损失的问题，将CSEM正演分解为总场算法和二次场算法。二次场算法能够去除场源奇异性，但无法实现任意复杂地形条件下地电模型的正演模拟^[4-9]。基于总场算法的求解方案早期主要是以伪delta函数刻画场源问题，一定程度上降低了场源奇异性，但无法从根本上解决场源的奇异性问题^[10-12]，甚至会加剧场源加载的难度。为此，得益于近年来非结构化网格离散技术的发展，广泛应用局部加密技术降低场源奇异性，结合delta函数场源积分能够避免场源加载的困难^[13-17]。随后，基于自适应的有限元正演技术被广泛应用于CSEM三维正演模拟，进一步弥补了场源奇异性带来的精度损失^{[14, 18]}。在实现大规模的电磁快速正演求解方面，早期基于节点有限元结合散度矫正的技术被应用，虽然弥补了电场方程迭代求解不收敛现象，但同时降低了正演求解的精度^[7]。之后，基于电场双旋度方程的矢量有限元直接求解技术被提出，虽大幅度提高了正演方程的求解精度^{[13, 16, 19]}，但无法避免双旋度结构带来的空解问题，使常规Krylov子空间迭代求解困难，导致直接求解的三维反演算法对内存的需求更高，限制了该算法的应用。为了使电场双旋度方程可用于迭代计算，Schwarzbach等^[14]开展了E-Φ方程的CSEM正演模拟(E和Φ分别表示电场和电标量位)，取得了一定效果，但亦无法从根本上解决迭代收敛慢等问题。之后，有学者提出基于势位(如A-Φ势^{[13, 20-25]}，A表示磁矢量位)的三维电磁正演求解系统，在一定程度上扩展了三维电磁求解系统的种类，计算精度和效率得到一定的改善。

综上所述，目前三维电磁感应模拟求解系统常用公式包括电场公式、磁场公式及势位公式(如矢量位和标量位)。在低频段，磁场公式对于空气空间没有意义，导致磁场公式很难得到正确的解^[26]。关于电场公式，双旋度算子离散系统矩阵会导致空解，因此最后系统矩阵会出现较高的条件数，需利用稳健的直接求解器PARDISO进行求解^[27]。目前应用较广的两类位场公式是T-Ω^[28]公式(T和Ω分别表示电矢量位和磁标量位)和A-Φ公式^{[13, 20-25]}。T-Ω位场公式需假设总电流密度在空气中为零，因此很难处理高频电磁感应问题。与之相反，A-Φ公式对于宽频电磁感应问题的研究非常有效，在低频电磁感应研究中得到广泛应用。基于A-Φ公式的有限元求解技术可大致分解成两类：一类是采用矢量基函数和节点基函数的离散求解系统，即基于高斯定理，用矢量位A表示磁场，用矢量位A和标量位Φ表示电场。通过Faraday’s定理推导基于双旋度结构的A-Φ耦合势方程，并采用矢量有限元法进行求解；另一类是通过强加Coulomb规范条件，使双旋度结构算子转化为拉普拉斯结构算子，从而得到拉普拉斯结构的A-Φ耦合势方程。Coulomb规范条件被强加▽·A=0，以保证A在整个求解区域是连续的，采用节点有限元求解整个系统。目前未见文献对基于电场方程及前述两种A-Φ方程的CSEM三维正演进行较系统的对比分析。因此，本文对基于电场方程的矢量有限元及两种结构的A-Φ耦合势CSEM求解系统进行对比研究，为其后续应用提供参考。

1 基本理论 1.1 CSEM偏微分方程

在准静态条件下，可控源电磁法满足Maxwell方程组，取时间因子e^-iωt，其频率域表达式为

$ \nabla \times \boldsymbol{E}=-\zeta \boldsymbol{H} $

(1)

$ \nabla \times \boldsymbol{H}=\chi \boldsymbol{E}+\boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}} $

(2)

$ \nabla \cdot \boldsymbol{B}=0 $

(3)

$ \nabla \cdot \chi \boldsymbol{E}=-\nabla \cdot \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}} $

(4)

式中：H为磁场；B=μH为磁感应强度，μ为磁导率；J_e^s表示电性源的电流分布；ζ=-iωμ表示阻抗，ω为角频率，${\rm{i}} = \sqrt { - 1} $为虚数单位；χ=σ-iωε为导纳率，σ表示电导率，ε为介电常数。

为消除磁场参数，将式(1)代入式(2)，可得

$ \nabla \times \frac{1}{\zeta} \nabla \times \boldsymbol{E}+\chi \boldsymbol{E}=-\boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}} $

(5)

上式即是三维CSEM需满足的电场方程。

有限元技术的基函数经历了节点基函数到矢量基函数的发展。为了消除伪解并遵循电场法向不连续性特征，基于矢量基函数的有限元技术求解电场方程成为目前主流方法。但是，可控源电磁法中的梯度电场是由地表和地下电导率不连续界面上的积累电荷产生的，导致基于电场方程的矢量有限元法很难实现梯度电场的模拟，而该梯度电场对于可控源电磁法非常重要，这一矛盾被称为矢量有限元的空解难题。因此，有学者将目光聚焦到A-Φ耦合势方程，进行了较系统的研究。

电场E和磁感应强度B可用矢量位A和标量位Φ分别表示为

$ \boldsymbol{B}=\nabla \times \boldsymbol{A} $

(6)

$ \boldsymbol{E}=\mathrm{i} \omega(\boldsymbol{A}-\nabla \varPhi) $

(7)

将式(6)和式(7)代入式(5)可得A-Φ耦合势方程

$ \nabla \times \nabla \times \boldsymbol{A}+\zeta \chi(\boldsymbol{A}-\nabla \varPhi)=\mu \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}} $

(8)

根据电荷守恒以及欧姆定律，可得辅助方程

$ \nabla \cdot[\zeta \chi(\boldsymbol{A}-\nabla \varPhi)]=\nabla \cdot\left(\mu \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}}\right) $

(9)

联立式(8)和式(9)，可得基于双旋度结构的A-Φ耦合势可控源电磁方程

$ \left\{\begin{array}{l} \nabla \times \nabla \times \boldsymbol{A}+\zeta \chi(\boldsymbol{A}-\nabla \varPhi)=\mu \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}} \\ \nabla \cdot[\zeta \chi(\boldsymbol{A}-\nabla \varPhi)]=\nabla \cdot\left(\mu \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}}\right) \end{array}\right. $

(10)

上式采用矢量基函数和节点基函数分别离散矢量位A和标量位Φ，通过矢量基函数实现对矢量位A进行隐式散度规范(▽·A=0)。除此之外，亦可通过显式强加Coulomb规范条件▽·A=0，对式(10)的矢量位A的双旋度结构进行化简，即

$ \nabla \times \nabla \times \boldsymbol{A}=\nabla(\nabla \cdot \boldsymbol{A})-\nabla^{2} \boldsymbol{A} $

(11)

最后形成系统方程

$ \left\{\begin{array}{l} -\nabla^{2} \boldsymbol{A}+\zeta \chi(\boldsymbol{A}-\nabla \varPhi)=\mu \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}} \\ \nabla \cdot[\zeta \chi(\boldsymbol{A}-\nabla \varPhi)]=\nabla \cdot\left(\mu \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}}\right) \end{array}\right. $

(12)

式(12)的求解系统采用节点基函数对矢量位A和标量位Φ进行离散，可获得拉普拉斯结构的A-Φ可控源电磁法求解方程。另外，本文通过扩边加载截断边界条件，实现电磁场在无穷远处消失为零。

1.2 有限元分析

以式(5)为控制方程，应用Galerkin法推导有限元方程。设余量为

$ \boldsymbol{r}_{\mathrm{e}}=\nabla \times \frac{1}{\zeta} \nabla \times \boldsymbol{E}+\chi \boldsymbol{E}-\boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}} $

(13)

令r_e在计算区域T满足

$ \iiint_{T} \boldsymbol{N} \cdot \boldsymbol{r}_{\mathrm{e}} \mathrm{d} T=0 $

(14)

式中N为矢量基函数。将式(13)代入式(14)，可得

$ \iiint_{T} \boldsymbol{N} \cdot\left(\nabla \times \frac{1}{\zeta} \nabla \times \boldsymbol{E}+\chi \boldsymbol{E}-\boldsymbol{J}_{\mathrm{s}}^{\mathrm{e}}\right) \mathrm{d} T=0 $

(15)

对式(15)应用第一矢量格林定理，可将其简化为

$ \begin{aligned} &\iiint_{T} \frac{1}{\zeta}(\nabla \times \boldsymbol{E} \cdot \nabla \times \boldsymbol{N}+\zeta \chi \boldsymbol{E} \cdot \boldsymbol{N}) \mathrm{d} T+ \\ &\iint_{\varGamma} \boldsymbol{N} \cdot\left(-\boldsymbol{n} \times \frac{1}{\zeta} \nabla \times \boldsymbol{E}\right) \mathrm{d} \varGamma=\iiint_{T} \boldsymbol{N} \cdot \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}} \mathrm{d} T \end{aligned} $

(16)

式中：Γ=∂T表示求解区域外边界；n为外边界单位法向量。另外，扩边处理使得求解区域足够大，因此式(16)中的面积分项可忽略不计。

同理，对式(10)应用Galerkin法进行推导并化简，最终的有限元积分方程为

$ \begin{aligned} &\iiint_{T}\left(\nabla \times \boldsymbol{N} \cdot \nabla \times \boldsymbol{A}+\zeta \chi{\boldsymbol{N}} \cdot \boldsymbol{A}\right) \mathrm{d} T- \\ &\quad\iiint_{T} \boldsymbol{N} \cdot \zeta \chi \nabla \varPhi \mathrm{d} T+\iint_{\varGamma} \boldsymbol{N} \cdot(\boldsymbol{n} \times \boldsymbol{B}) \mathrm{d} \varGamma \\ &\quad=\iiint_{T} \boldsymbol{N} \cdot\left(\mu \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}}\right) \mathrm{d} T \end{aligned} $

(17)

$ \begin{aligned} &\zeta \chi \iiint_{T}(\nabla L \cdot \boldsymbol{A}+\nabla L \cdot \nabla \varPhi) \mathrm{d} T+ \\ & \quad \mu\chi \iint_{\varGamma} L(\boldsymbol{E} \cdot \boldsymbol{n}) \mathrm{d} \varGamma=-\iiint_{T} L\left(\nabla \cdot \mu \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}}\right) \mathrm{d} T \end{aligned} $

(18)

式中L表示有限元中的节点基函数。

同理，对式(12)应用Galerkin法进行公式推导及化简，得到最终的有限元积分表达式为

$ \begin{array}{r} \iiint_{T}\left(\nabla L \cdot \nabla A_{x}+\zeta \chi L A_{x}-\zeta \chi L \frac{\partial \varPhi}{\partial x}\right) \mathrm{d} T \\ =\iint_{\varGamma} L\left(\nabla A_{x} \cdot \boldsymbol{n}\right) \mathrm{d} \varGamma+\iiint_{T} L J_{\mathrm{s}}^{x} \mathrm{~d} T \end{array} $

(19)

$ \begin{array}{r} \iiint_{T}\left(\nabla L \cdot \nabla A_{y}+\zeta \chi L A_{y}-\zeta \chi L \frac{\partial \varPhi}{\partial y}\right) \mathrm{d} T \\ =\iint_{\varGamma} L\left(\nabla A_{y} \cdot \boldsymbol{n}\right) \mathrm{d} \varGamma+\iiint_{T} L J_{\mathrm{s}}^{y} \mathrm{~d} T \end{array} $

(20)

$ \begin{array}{r} \iiint_{T}\left(\nabla L \cdot \nabla A_{z}+\zeta \chi L A_{z}-\zeta \chi L \frac{\partial \varPhi}{\partial z}\right) \mathrm{d} T \\ =\iint_{\varGamma} L\left(\nabla A_{z} \cdot \boldsymbol{n}\right) \mathrm{d} \varGamma+\iiint_{T} L J_{\mathrm{s}}^{z} \mathrm{~d} T \end{array} $

(21)

式中：A_x、A_y、A_z分别表示矢量位A的三个分量；J_s^x、J_s^y、J_s^z分别表示电性源J_e^s的三个分量。

下面分别叙述如何选择上述三类方程有限元基函数。

(1) 对电场方程采用矢量有限元进行离散。对于每一个四面体单元e，任意点x, y, z处的电场可由各条棱边的电场E_i和矢量基函数N_i表示为

$ \boldsymbol{E}=\sum\limits_{i=1}^{6} E_{i} \boldsymbol{N}_{i} $

(22)

式中矢量基函数N_i可用节点基函数^[29]表示为

$ \boldsymbol{N}_{i}=\left(L_{j_{1}} \nabla L_{j_{2}}-L_{j_{2}} \nabla L_{j_{1}}\right) l_{i} $

(23)

式中：l_i为棱边i的长度；L_j1与L_j2表示该棱边两个节点的基函数，具体标号如图 1所示。

(2) 对于双旋度结构的A-Φ耦合势CSEM方程，矢量位A和标量位Φ可分别采用矢量基函数和节点基函数进行描述

$ \left\{\begin{array}{l} \boldsymbol{A}=\sum\limits_{i=1}^{6} A_{i} \boldsymbol{N}_{i} \\ \varPhi=\sum\limits_{j=1}^{4} \varPhi_{j} L_{j} \end{array}\right. $

(24)

(3) 拉普拉斯结构的A-Φ耦合势CSEM方程，矢量位A和标量位Φ可采用节点基函数表示为

$ \left\{\begin{array}{l} A_{x}=\sum\limits_{j=1}^{4} L_{j} A_{x, j} \\ A_{y}=\sum\limits_{j=1}^{4} L_{j} A_{y, j} \\ A_{z}=\sum\limits_{j=1}^{4} L_{j} A_{z, j} \\ \varPhi=\sum\limits_{j=1}^{4} L_{j} \varPhi_{j} \end{array}\right. $

(25)

因此，通过以上不同的离散形式可得到三种不同的CSEM求解系统，即电场方程、双旋度结构的A-Φ耦合势和拉普拉斯结构的A-Φ耦合势，本文结合稳健的直接求解器PARDISO和Krylov子空间算法对这三种求解系统进行对比分析^{[27, 30]}。

1.3 场源加载

对于CSEM三维正演来说，源项加载不准确和场源的奇异性会严重影响正演的求解精度。针对这两种问题，常用的场源处理技术主要包含两类：①二次场算法，其目的是解决源的奇异性问题，但需要计算具有解析表达式的一次场，因而对复杂地电模型模拟的局限性较大；②基于总场算法并结合偶极源积分实现场源加载，具有较好的适用性，但仍无法避免场源奇异性问题。为此，本文采用偶极源积分结合非结构化网格局部细化技术，不仅能有效降低场源奇异性，还能实现场源的统一表示形式。

对于电场方程的三维CSEM有限元系统，场源加载可表示为

$ \left\{\begin{array}{l} S_{1}^{\mathrm{e}}=\iiint_{T}\left(\boldsymbol{N}_{\mathrm{i}} \cdot \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}}\right) \mathrm{d} T \\ S_{1}^{\mathrm{m}}=\iiint_{T}\left[\boldsymbol{N}_{\mathrm{i}} \cdot\left(\nabla \times \boldsymbol{J}_{\mathrm{m}}^{\mathrm{s}}\right)\right] \mathrm{d} T \end{array}\right. $

(26)

式中：S₁^e、S₁^m分别表示电性源和磁性源的加载积分形式；J_m^s表示磁性源的电流分布。

对于x方向电性源来说，源线段置于四面体棱边上，因此有

$ \begin{aligned} &S_{1}^{\mathrm{e}}=\iiint_{T} \boldsymbol{N}_{i} \cdot \boldsymbol{J}_{\mathrm{e}}^{\mathrm{s}} \mathrm{d} T=\frac{I}{\left(6 V^{\mathrm{e}}\right)^{2}} \times \\ &\ \ \iiint_{T}\left(N_{x} \boldsymbol{i}_{x}+N_{y} \boldsymbol{i}_{y}+N_{z} \boldsymbol{i}_{z}\right) \cdot\left[\varTheta\left(x_{i+1}\right)-\right. \\ &\ \ \left.\varTheta\left(x_{i}\right)\right] \delta\left(y-y_{0}\right) \delta\left(z-z_{0}\right) \boldsymbol{i_{x}} \mathrm{~d} T \end{aligned} $

(27)

式中：Θ表示Heaviside函数；δ是delta函数；N_x、N_y、N_z分别表示矢量基函数N的三个分量；i_x、i_y、i_z表示笛卡尔坐标系下沿坐标轴方向的单位向量；x₀, y₀, z₀表示偶极源的中心坐标。

垂直磁偶源M的表达式为

$ \boldsymbol{M}=m \delta\left(x-x_{0}\right) \delta\left(y-y_{0}\right) \delta\left(z-z_{0}\right) \boldsymbol{i}_{z} $

(28)

式中m表示磁矩。磁性源的电流分布J_m^s的表达式为

$ \boldsymbol{J}_{\mathrm{m}}^{\mathrm{s}}=\mathrm{i} \omega \mu \boldsymbol{M} $

(29)

因此，S₁^m的场源积分表达式为

$ S_{1}^{\mathrm{m}}=\mathrm{i} \omega \mu \iiint_{T} \boldsymbol{N}_{i} \cdot(\nabla \times \boldsymbol{M}) \mathrm{d} T $

(30)

由上式可得

$ S_{1}^{\mathrm{m}}=\mathrm{i} \omega \mu \iiint_{T} \boldsymbol{M} \cdot\left(\nabla \times \boldsymbol{N}_{i}\right) \mathrm{d} T-\iint_{\varGamma} \boldsymbol{N}_{i} \times \boldsymbol{M} \cdot \boldsymbol{n} \mathrm{d} \varGamma $

(31)

将场源置于四面体单元内，上式中的面积分项(公式右边第二项)等于零，因此积分式只剩下右边第一项。若采用电性源的形式加载磁性源，场源的积分公式与电性源一样，在此不再叙述。

2 算例分析

为了分析上述三种求解系统对三维CSEM的模拟效果，设计了两个模型进行正演计算。一个模型是存在解析解的均匀半空间，另外一个是无解析解的低阻地电模型。测试本文开发算法的计算机配置为：内存32G，intel(R) Core^TM i7-9700 CPU@ 3.00GHz。

2.1 模型1：均匀半空间 2.1.1 算法验证

对均匀半空间模型加载一个正方形的垂直磁偶源，磁偶源是大小为1.0m×1.0m的线框，位于求解区域的中心。地下背景电阻率为100Ω·m，空气电阻率为10⁸Ω·m。求解区域大小为[-30km, 30km]³。为了降低场源的奇异性，对场源进行局部加密，被离散为27个小的线段，如图 2所示。

计算3Hz的垂直磁偶源电磁响应。整个求解区域被剖分成403789个四面体单元，471329条边和66521个节点。沿x方向等间距设置观测点，点距为120m，共50个观测点。

采用Krylov子空间的GMRES算法结合ILU(0)预条件因子分别对前述两种结构的A-Φ耦合势有限元方程进行求解，采用PARDISO直接求解器对电场方程进行求解，得到磁场垂直分量H_z的实部和虚部，并与解析解进行对比(图 3)。两种结构的A-Φ耦合势方程组求解的收敛曲线见图 4。

从图 3可知，三种求解系统的结果与解析解都高度吻合，误差基本都小于1.5%。从图 4的收敛曲线可知，双旋度结构的A-Φ耦合势方程组求解需要迭代3000次，相对残差为2.31×10^-15，而拉普拉斯结构的A-Φ耦合势方程组求解仅需200次，相对残差为1.9×10^-15。

本文还计算了100Hz时垂直磁偶源的电磁响应。求解区域被离散为814255个四面体单元、675519条边单元和138736个节点，并沿x方向等间距设置观测点，间距为20m，共50个观测点。图 5展示了100Hz时的电磁响应。两种结构的A-Φ耦合势收敛曲线如图 6所示。从图 6可知，双旋度结构的A-Φ耦合势迭代3000次未达到预设的求解精度，而拉普拉斯结构的A-Φ耦合势经迭代2400次即达到预设精度。虽然前者得到的垂直磁场H_z与解析解吻合较好，但求解的迭代次数较多。对比分析表明，本文开发的三种CSEM求解方程算法准确，结果可靠。

2.1.2 求解系统特点分析

为了分析三种CSEM求解方程的特点，仍然采用上节的均匀半空间地电模型。沿x方向、在地面放置一个电偶源，源长度为1m，发射电流为1A，发射频率为10Hz。求解区域的大小为[-30km, 30km]³，分别采用前述三种求解方程对该模型进行求解，地面上各场分量等值线图见图 7~图 9。

由图 7~图 9可以看出，基于双旋度结构的A-Φ耦合势计算的矢量和标量位无明显规律，基于拉普拉斯结构的A-Φ耦合势计算的矢量和标量位特征符合预期。通过两种A-Φ耦合势计算的电场E_x等值线与基于电场方程计算的电场E_x等值线具有高度一致性，表明矢量基函数离散的矢量位A违背了连续性条件，导致矢量位A切向不连续，使得计算得到的矢量和标量位无明显规律，影响其收敛性。

2.2 模型2：低阻模型

设计图 10所示的长方体低阻异常体模型^[23]。块状异常体(5Ω·m)被置于电阻率为50Ω·m的均匀半空间，异常体尺寸为120m×200m×400m，中心点坐标为(1000m, 0m, 300m)。沿x方向布设有限长导线源，源长度为100m，中心点坐标为(50m, 0m, 0m)，发射电流为1A。沿着x方向布设测线，起点和终点位置分别为(400m, 0m, 0m)、(1400m, 0m, 0m)，点距为20m。利用Krylov子空间算法的GMRES和BICGSTAB迭代求解器分别对三种求解系统进行求解，分析其收敛性。3Hz和10Hz对应的矢量位分量A_x、标量位▽_xΦ及电场分量E_x计算结果见图 11和图 12。

从图 11和图 12可见，两种结构的A-Φ耦合势方程计算的电场E_x实部和虚部分量与电场方程的结果高度吻合。然而，两者矢量位分量的实部和虚部存在明显的差异，基于双旋度结构的A-Φ耦合势方程计算的A_x分量的实部和虚部曲线不光滑，标量位▽_xΦ的虚部也存在类似现象。同样的规律也发生10Hz的电磁响应上。从上述的分析可知，基于有限元矢量基函数自然可以保证每个单元的散度为零，但由于数值解误差的存在，不能完全保证散度自由，从而导致类似现象出现。基于拉普拉斯结构的A-Φ耦合势强加库伦规范条件保证了方程的唯一性，使得方程的解较准确且曲线光滑。

3Hz所对应的基于三种求解系统的不同求解器以及预条件因子的收敛性能统计结果见表 1。由表可知，基于电场方程的CSEM矢量有限元方程在常规迭代求解器和预条件因子下均不能得到满意的计算结果。在同等精度、迭代求解器条件下，基于拉普拉斯结构的A-Φ耦合势CSEM求解方程收敛性总体上优于基于电场方程和双旋度结构的A-Φ耦合势方程。对于同一种求解方程，GMRES和BICGSTAB迭代求解器的求解效果相差不大。除此之外，GMRES和BICGSTAB迭代求解器结合SOR(超松弛)预条件因子对两种结构的A-Φ耦合势CSEM求解方程进行求解，收敛性差异不明显；而GMRES和BICGSTAB迭代求解器结合其他预条件因子进行求解的收敛性明显变差，所需迭代次数明显增加。这说明选择合适的预条件因子可提高方程组的求解效率。

表 2是不同求解方程/求解器组合下CSEM求解内存消耗、求解时间、迭代次数及相对残差的统计。对于方程组未知数在百万级以内的情况来说，直接求解器求解电场方程的优势明显，小型计算机即可满足内存需求；若需求解的未知数超过100万，计算机内存需求则超过15GB；当需求解的未知数继续增加，利用直接求解器求解线性方程组的计算机内存需求则难以满足。对于A-Φ耦合势求解方程，当未知数超过100万，内存需求不超过6GB，即A-Φ耦合势方程的迭代解法在内存需求方面的优势非常明显。但是，在时间消耗方面，对于未知数不超过100万的方程组求解，直接求解器相比迭代求解器优势更明显。

3 结论

本文采用非结构化网格对求解区域进行离散，实现了三种求解系统的CSEM正演计算。建立均匀半空间模型和低阻异常体地电模型，从计算精度、效率以及内存消耗三个方面分析基于两种A-Φ耦合势方程和电场方程正演问题的求解特性。得出以下几点认识。

(1) 本文开发的三种正演求解算法程序准确且计算精度较高，采用的统一场源处理技术能够较好地扩展场源的适用性，同时结合局部网格加密技术可有效降低场源奇异性问题。

(2) 在同等条件下，基于非结构网格的矢量有限元电场方程不适合于采用常规的Krylov子空间算法进行方程组的求解；相比于双旋度结构，基于拉普拉斯结构的A-Φ耦合势更适合于电磁场的迭代求解，其适用性更强。

(3) 相比于A-Φ耦合势方程，基于电场方程的直接解法在方程组未知数不高于100万时，在内存需求和求解效率方面具有明显的优势；但是，若方程组未知数超过100万级，相比直接解法，A-Φ耦合势方程迭代算法在内存需求方面优势更明显。另外，若网格密度增加，常规的Krylov子空间迭代算法求解大型方程组所需要的时间成本会逐渐增加。

总之，本文系统对比了三种频率域CSEM正演求解系统，并对正演模拟的内存需求及效率进行了探讨。对于小尺度三维CSEM正演问题，电场方程的直接解法能满足需求；对于超大尺度三维CSEM正演问题，需要借助迭代算法开展线性方程组求解。