大地电磁法由于利用天然交变电磁场源，具有数据采集方便、勘探深度大等优点.目前该方法已经被广泛应用于地质普查、矿产勘查、油气勘探、地壳和上地幔电性结构研究等方面.大地电磁正演模拟作为反演和数据解释的基础，已经得到了国内外众多学者的重视.传统基于结构化网格的大地电磁正演模拟技术已十分成熟(Ting and Hohmann，1981；Mackie et al.，1993；Siripunvaraporn et al.，2002；Mitsuhata and Uchida，2004；Nam et al.，2007；Tong et al.，2009).随着计算机的发展以及人们对正演模拟精度和模型复杂程度要求的提高，近年来基于非结构化网格的大地电磁正演模拟算法得到了良好的发展(Frank et al.，2007a；Liu et al.，2008；Usui，2015).三维电磁正演通常模型规模较大，需要巨大的计算资源.合理的网格可以利用较少的单元达到较高的正演模拟精度，达到节省计算资源的目的.为此，自适应网格剖分技术近年来被逐渐应用于求解电磁正演模拟问题(Key and Weiss，2006；Frank et al.，2007b；Ren et al.，2013).

电磁正演模拟方法通常包括有限体积法(杨波等，2012; Jahandari and Farquhason，2014)、有限差分法(Commer and Newman，2004; Streich，2009)、积分方程法(Zhdanov et al.，2006)和有限单元法(Ren et al.，2013; 殷长春等，2015)等.这些方法中，有限差分法必须使用阶梯网格来模拟复杂介质边界，为了得到较高精度的正演模拟响应，必须将网格划分得很密，从而耗费大量计算时间.积分方程法因其矩阵方程是密实的，无法对复杂模型响应进行模拟.有限体积法和有限单元法可以采用非结构化网格对模型进行任意剖分，正在被广泛应用于复杂模型正演模拟.本文选择使用基于非结构化网格的有限单元法对大地电磁进行正演模拟.有限单元法按照插值函数可分为节点(标量)有限元法和棱边(矢量)有限元法.标量有限元法因插值函数在单元内不满足散度为零条件，容易出现非物理解.此外，由于该方法在物性分界面处无法保证场值的连续性，强加边界条件十分不便(金建铭，1998).矢量有限单元法在单元内能够保证形函数自然满足散度条件，且在物性分界面满足场值连续性，无需强加边界条件.本文选择使用矢量有限单元法对大地电磁正演问题进行求解.

自适应算法可分为面向目标的自适应算法和全局自适应算法.面向目标的自适应算法与全局自适应算法相比，由于使用加权的后验误差判断网格是否需要被优化，因而能够对计算区域进行更合理的划分.自适应算法目前较为常用的算法主要有两种.一种是基于场的梯度的超收敛性，通过比较直接求解出的梯度与修正后的梯度之间的差距估算后验误差(Key and Weiss，2006；Schwarzbach et al.，2011)；另一种是基于场和电流的连续性条件，通过计算求解出的场或电流在界面处与连续性条件之间的差距估算后验误差(Ren et al.，2013).本文利用垂向电流密度在物性分界面的连续性对后验误差进行估计；而加权系数则通过求解正演问题的对偶问题计算得出.由于对偶问题的解反映了电场分布，而上述方法后验误差估计使用的是电流密度，因此本文对加权系数进行改进，使用电导率与对偶问题解的乘积作为最终的加权系数以确保两者之间的统一.

为了验证本文算法的精度，我们模拟了均匀半空间模型的单点大地电磁响应，分析了误差下降规律，并比较了本文面向目标的自适应算法与全局自适应算法的正演模拟效果.为了验证本文算法对起伏地表问题的有效性，我们模拟并比较了Nam等(2007) 给出的经典模型的大地电磁响应.随后，通过计算起伏地表下埋有异常体模型的大地电磁响应，展示本文算法不仅能够对接收点附近的起伏地表进行合理加密，同时还能对具有一定埋深的地下介质分界面进行有效加密.最后，通过对比未改进的加权系数算法和本文算法对起伏地表下埋有异常体模型的网格剖分效果证明本文改进加权系数算法的优越性.

本文讨论的面向目标大地电磁自适应正演模拟算法可分为以下三个部分：1) 基于非结构化网格的有限元大地电磁正演模拟算法； 2) 加权后验误差估计算法；3) 网格优化算法.

对于大地电磁正演问题，电磁场可由如下形式的Maxwell方程组表示：

(1)

(2)

式中σ是电导率，E是电场，H是磁场，μ是磁导率，ω是角速度.对公式(1) 和(2) 消去磁场，可以得到三维大地电磁对应的电场Helmholtz方程：

(3)

为保证(3) 式具有唯一解，我们在计算区域外边界添加了边界条件.假设在计算区域上边界，存在沿x方向的极化电场，则上边界面的边界条件为：

(4)

式中Ω表示计算区域，E_source表示上边界面的场值.对于其他边界面，有限元方程中无需做任何处理(徐世浙，1994).(3) 和(4) 式组成了大地电磁正演的边值问题.定义希尔伯特矢量空间：H(curl，Ω)={Φ∈L₂(Ω)，Δ×Φ∈L₂(Ω)}，其中L₂是希尔伯特空间下的平方可积函数，并在希尔伯特空间下定义内积：

(5)

(6)

将公式(3) 与矢量插值函数做内积，并应用格林公式，可以得到：

(7)

其中Φ∈H(curl，Ω)是矢量插值函数.将计算区域离散成小单元，并在每个单元内使用有限维函数

(8)

代替E，最终可以得到如下形式的大型线性方程组：

(9)

求解上述线性方程组，可以得到计算区域内各个棱边的电场值，区域内任意位置的电场值可以通过公式(8) 得出，而磁场则可以通过Faraday′s定律得出

(10)

仅由假设上边界存在沿x方向极化源计算出的电场和磁场无法确定阻抗张量：

(11)

为得到阻抗张量，还需要对假设上边界存在沿y方向极化源的情况进行正演模拟，此时阻抗张量可由(12) —(15) 式确定：

(12)

(13)

(14)

(15)

式中下角标1和2表示假设上界面存在沿x方向和y方向电场源的响应.阻抗张量公式的推导见李焱等(2012) .

定义τ为单元Γ_i和Γ_j之间的面，使用τ两侧单元分别计算出的面上电场值为E₁和E₂.在界面处，法向电流满足连续性(纳比吉安，1992)：

(16)

(17)

其中n_τ表示τ面的法方向，σ₁和σ₂表示面两侧单元的电导率，J₁和J₂表示从两侧流过面的电流密度.在无限维希尔伯特空间中，电流密度严格遵守公式(16) ，然而在有限维希尔伯特空间中，电流密度无法保证严格遵守(16) .本文利用有限维希尔伯特空间中电流密度在单元界面处不连续估计后验误差.后验误差的定义形式如下(Ren et al.，2013)：

(18)

其中Γ_i表示单元的第i个面.

使用上述后验误差对网格进行自适应优化称为全局自适应算法，该方法对所有后验误差较大的单元进行加密.然而在大地电磁正演模拟中，观测点往往仅包含于少量单元内，这些观测点处的数值模拟精度是人们关心的焦点.全局自适应算法除了加密对观测点模拟精度影响较大的单元外，对其他单元也进行了加密，浪费了大量的计算资源.面向目标的自适应算法重点对能使观测点数值模拟精度提高的网格进行加密，对不影响观测点处数值模拟精度的网格加密力度小.与全局自适应算法相比，面向目标的自适应算法可以使用更少的网格达到更高的计算精度.

定义评价观测点所在单元内电场大小的线性函数L(E)，则自适应算法的目标为减小误差函数L(e)，其中e=E－E_h表示正演模拟解与精确解之间的差距.上述误差函数与加权系数w存在如下关系(Oden and Vemaganti，2000)：

(19)

其中B^*(w，Φ)表示B(w，Φ)的对偶形式，B(w，Φ)具体表达式如下：

(20)

从公式(19) 可以发现，w近似表达的是电场分布，而估计的后验误差与电流密度对应.为了使两者统一，本文重新定义加权系数为：w′=σw.w在空间的变化规律与在接收点处放置源的广义格林函数相似；而本文新定义的w′由于考虑了电导率的影响，与w相比在电场衰减更快的低阻介质内数值更大，从而使该区域的网格更容易被加密.为了表达方便，在后文将上角标撇号省略.根据公式(19) 以及Cauchy-Schwartz不等式，可以对L(e)的上限进行估计(Oden and Prudhomme，2001)：

(21)

式中表示正的常数值，N_all表示与L(e)有关的单元个数.式(21) 中w可由式(19) 得出，然而由于电场的精确解无法得到，误差e难以被直接计算.由于(18) 式反映了因使用有限维插值函数代替无限维函数而带来的计算误差，我们使用η_e代替‖e‖_L2.此时加权后验误差变为：

(22)

其中η_w，n=‖w‖_L2.由于L()是衡量物理量在单元内大小的线性函数，定义公式(19) 中的L(Φ)为(Oden and Vemaganti，2000)：

(23)

式中V_j表示第j个单元的体积，I为单位向量，t表示目标测点所在的单元个数，目标测点的定义将在后文给出.

在上面的内容中，我们描述了大地电磁正演模拟算法以及加权后验误差的计算方法，下面将根据上述内容整体介绍本文的面向目标自适应算法流程.

本文提出的自适应算法迭代终止准则是基于各个测点模拟结果的收敛性，即比较本次正演结果与之前正演结果之间的差距，若两者相差很小，则认为该测点模拟结果已经具有较高精度(收敛准则)，该测点对应的网格无需进行加密；否则认为该测点模拟精度不能达到要求，对应的网格仍需继续加密.不满足上述收敛准则的测点定义为目标测点，否则定义为收敛测点.本文面向目标自适应算法过程如下：

(1) 设定正演问题的最大迭代次数和最大棱边数；

(2) 如果迭代次数小于最大迭代次数，且棱边个数小于最大棱边数，则执行步骤(3) ，否则执行步骤(6) ；

(3) 对每个单元定义相对加权后验误差β_n=，并定义β_n的阀值为β_t.为保证程序的稳定性，对所有单元添加了最小体积约束V_min.若一个单元满足β_n＞β_t和V＞V_min，则对该单元进行优化；

(4) 使用重新划分的网格进行正演模拟；

(5) 比较本次正演结果与前两次正演结果的差距，并记录目标测点和收敛测点.若所有测点均为收敛测点，则执行步骤(6) ，否则执行步骤(2) ；

(6) 结束循环.

本文数值模拟计算环境为Intel(R)Xeon(R)CPU E5-2650 V2 @ 2.60GHz 2.60 GHz(2处理器)，内存为128 GB.为了验证本文算法的精度，并比较面向目标的自适应算法与全局自适应算法的效果，本文使用这两种算法对均匀半空间大地电磁响应进行数值模拟，并比较两种算法的精度随自适应迭代次数增加的变化规律.均匀半空间模型电阻率为100 Ωm，大地电磁信号频率为2 Hz.本文所有模型均假设介质磁导率为μ=μ₀.图 1a和1b给出了两种算法响应的相对误差曲线，图 1c和1e给出了两种算法的网格剖分，而图 1d和1f分别给出了图 1c和1e由白色方框标记区域网格剖分细节展示.从图 1a和1b给出的相对误差曲线可以发现，面向目标的自适应算法与全局自适应算法相比能够使用更少的网格达到更高的计算精度；面向目标自适应算法的视电阻率相对误差可下降到1%以下，而相位误差可下降到0.1°.从图 1c—1f 给出的网格划分结果可以看出，面向目标的自适应算法在观测点附近将网格划分得很密，而在远离观测点处，网格剖分比较稀疏.相比之下，全局自适应算法在整个地表均对网格进行了整体加密，而在接收点附近网格相对粗糙.两种方法的加密效果对比表明面向目标的自适应算法对网格的加密方式更为合理.

图 1

Fig. 1

图 1 均匀半空间模型网格剖分及响应误差 (a)面向目标及全局自适应算法视电阻率相对误差；(b)面向目标及全局自适应算法相位误差；(c)面向目标自适应算法网格剖分结果；(d)展示图(c)中白色方框内网格剖分细节；(e)全局自适应算法网格剖分结果；(f)展示图(e)中白色方框内网格剖分细节. Fig. 1 Meshes and results for a homogeneous half-space model (a)Relative errors of apparent resistivity calculated by goal-oriented and global adaptive methods;(b)Errors of phase calculated by goal-oriented and global adaptive methods;(c)Final mesh for goal-oriented adaptive method;(d)Local view of the mesh marked by a white rectangular in(c);(e)Final mesh for global adaptive method;(f)Local view of the mesh marked by a white rectangular in(e).

从图 1a和1b中可以发现误差曲线虽然整体呈下降趋势，但存在波动.这是由于上述算法仅考虑了数值模拟的精度，却忽略了插值带来的误差.随着迭代次数的增加，虽然数值模拟的计算精度有所提高，但由于包含观测点的单元一直在改变，因此插值误差也一直在变化.这种误差对于简单的半空间模型随迭代次数的增加会被压制，然而对于带有起伏地表的复杂模型，这种误差会被明显表现出来.为解决这一问题，我们人为地在初始网格剖分文件中对每个观测点附近添加了一个点，以减小包含观测点网格的体积，进而减小插值误差.

为了验证本文算法对起伏地表模型的有效性，我们对Nam(2007) 使用有限元方法计算的大地电磁模型进行了模拟并对比响应曲线.图 2a和2b分别给出了模型的俯视图和剖面图；图 2c—2f分别给出了视电阻率和相位响应对比曲线；而图 2g和2h给出了自适应网格剖分结果.从图 2c，2d，2e，2f可以看出，本文模拟结果与Nam模拟结果吻合较好，从而证明了本文算法对起伏地表模型具有较高的计算精度.从图 2g和2h可以进一步看出，本文算法对接收点附近的起伏地表进行了很好的网格加密，随着地下介质距离观测点和地表物性分界面的增加，网格逐渐变得稀疏.

图 2

Fig. 2

图 2 本文模拟结果与Nam 等 (2007)模拟结果对比 (a) 模型俯视图；(b) 模型切面图；(c) 和(d) 视电阻率响应结果对比；(e) 和(f) 相位响应结果对比；(g) 自适应算法网格剖分结果；(h) 展示图(g) 所示白色方框内网格剖分细节. Fig. 2 Comparison of results from this paper with those from Nam et al.(2007) (a) Plane view of the model; (b) Sectional view of the model; (c) and (d) Comparison of the apparent resistivities ρxya and ρyxa ; (e) and (f) Comparison of phases φxy and φyx; (g) Final mesh for the topography model; (h) Local view of the mesh marked by a white rectangular in (g).

从图 1和图 2的模型结果可以发现，本文提出的算法具有较高计算精度和较强处理起伏地表的能力，对网格的自适应优化也较为合理.然而，对于大地电磁正演问题，由于观测点位于具有较大物性参数对比的空气-大地分界面上，算法容易实现对地表附近网格的加密；但当起伏地表下埋有异常体时，合理的自适应网格不仅要对接收点附近的地表分界面进行加密，同时对与接收点有一定距离的异常体分界面也应该予以加密.本文通过引入电导率加权改进的后验误差权系数有效地实现地表分界面和地下介质分界面的合理自适应加密.

图 3给出起伏地表下埋有异常体模型的网格剖分及大地电磁响应.起伏地表模型参数与图 2相同，异常体埋深200 m，规模为800 m×800 m×800 m，电阻率为1 Ωm(图中红色方框).图 3a和3b给出直接使用(19) 式中w作为加权因子的网格剖分结果，而图 3c和3d给出利用本文改进的加权因子的网格剖分结果.表 1给出了本文改进算法每次迭代使用的网格信息和求解方程耗费的时间.从图 3a和3b可以发现，直接利用w作为加权因子的自适应算法没能将低阻异常体进行有效加密，这与电磁波在低阻介质中衰减梯度大，需要加密网格的规律不符.从图 3c和3d给出的利用本文加权因子的网格剖分结果可以发现，本文算法不仅在观测点附近的起伏地表处对网格进行了很好的加密，而且对具有一定埋深的异常体表面也进行了网格加密.此外，在异常体的顶部，网格加密程度要高于底部，这是由于随着埋深的增加，介质对观测点影响的程度逐渐变小，因而降低网格剖分的精细程度.比较图 3a—3d的剖分结果，可以明显看出本文改进算法的优势.从图 3e—3h给出的模拟结果可以看出，本文算法的响应曲线光滑且对称性较好.

图 3

Fig. 3

图 3 起伏地表下异常体模型网格剖分及大地电磁响应 (a) 未改进自适应算法网格剖分结果；(b) 展示图(a) 中白色方框内网格剖分细节；(c) 本文自适应算法网格剖分结果；(d) 展示图(c) 中白色方框内网格剖分细节；(e) 和(f) 视电阻率响应剖面曲线；(g) 和(h) 相位响应剖面曲线.图中红框表示异常体. Fig. 3 Adaptive meshes and MT responses on a topographic model with an embedded abnormal body The red rectangular denotes the abnormal body, (a) Final unimproved adaptive meshes; (b) Local view of the mesh marked by a white rectangular in (a); (c) Final adaptive meshes; (d) Local view of the mesh marked by a white rectangular in (c); (e) and (f) Apparent resistivities ρxya and ρyxa; (g) and (h) Phases φxy and φyx.

表 1(Table 1)

表 1 起伏地表下异常体模型网格参数及PARDISO方程求解器求解方程耗费时间 Table 1 Mesh parameters for the topographic model with an embedded abnormal body and consumed time of PARDISO solver 迭代次数 单元个数 棱边条数 方程求解时间(s) 1 21893 25677 0.691 3 101017 117423 2.978 5 244697 283859 8.765 7 484239 560991 21.86

表 1 起伏地表下异常体模型网格参数及PARDISO方程求解器求解方程耗费时间 Table 1 Mesh parameters for the topographic model with an embedded abnormal body and consumed time of PARDISO solver

基于非结构有限元方法可以模拟任意起伏地表，而自适应算法可以为非结构有限元方法提供合理网格的特点，本文成功地将面向目标的自适应算法应用于大地电磁非结构有限元正演模拟中.针对均匀半空间模型模拟结果表明，本文算法具有较高精度；面向目标的自适应算法比全局自适应算法能够更合理地对网格进行加密.对经典起伏地表模型的模拟结果表明，本文算法可以有效地对起伏地表等复杂物性界面进行加密，从而为取得较高精度的正演模拟结果提供保证.通过引入电导率加权改进的后验误差权系数，成功实现了对含有接收点的起伏地表分界面和地下具有一定埋深且物性对比较小的异常体表面进行有效加密.本文算法对复杂模型大地电磁模拟和反演具有指导意义.