1 引言

大地电磁法是一种利用天然交变电磁场源，来探测地球内部电性结构的地球物理方法(Nabighian, 1987)，被广泛应用于矿产资源勘查、深部地质构造探测和地球动力学研究(Xu et al., 2016; Newman et al., 2008; Unsworth, 2003; Tang et al., 2013).

作为数据解释的基础，大地电磁正演一直是国内外学者的研究热点.大地电磁正演方法主要有：有限差分法、有限单元法、有限体积法和积分方程法.有限差分法具有实现较易的优点(谭捍东等, 2003; 李焱等, 2012)，但是不适合模拟复杂异常体，且需要对整个区域进行网格离散.有限单元法能够模拟复杂模型(Mitsuhata and Uchida, 2004; Liu et al., 2008; Ansari and Farquharson, 2014; Ren et al., 2014; 殷长春等, 2017)，但是需要对全空间进行网格剖分，未知数个数较多.有限体积法适合处理高对比度问题(Du et al., 2016)，但是，与有限差分法和有限单元相同，需要对全空间进行网格离散，未知数个数较多.相比之下，积分方程法有如下的优点：1)只需要对异常体进行单元剖分，未知数少；2)积分方程结果具有半解析解的精度，常常被用于测试新开发的算法的精度；3)由于研究目标体为异常体，基于积分方程的反演算法具有非常高的反演效率和精度(Kruglyakov et al., 2016).因此，开发高精度的三维大地电磁积分方程正演算法仍然具有一定的研究价值.

目前，积分方程法主要用于可控源电磁法的正演计算(Zhdanov et al., 2006; Gao and Torres-Verdín, 2006; 陈桂波等, 2009, 2011, 2013)，只有为数不多的研究人员从事了大地电磁问题的积分方程求解，主要集中在基于规则六面体单元的并矢Green函数计算、密实矩阵的快速迭代解法和积分方程近似解法.Hohmann在1975年基于六面体单元使用积分方程法解决了三维电磁正演问题，打开了三维解释的大门，之后积分方程法被广泛应用到三维电磁数值模拟中(Hohmann, 1975).继承着Hohmann教授的工作，Ting采用结构化网格计算了三维大地电磁法的积分方程数值解(Ting and Hohmann, 1981).Wannamaker给出了背景为层状水平地层的积分方法解(Wannamaker et al., 1984).1991年，Wannamaker对积分方程法在地球物理电磁感应的可能应用进行了综合评述(Wannamaker, 1991).张辉等2005年使用高斯求积和连分式展开计算电磁并矢Green函数积分(张辉等, 2005).鲁来玉等研究了电阻率随位置线性变化的三维大地电磁模拟(鲁来玉等, 2003).Xiong，Mackie等人着重讨论了积分密实矩阵的快速迭代求法(Xiong, 1992; Mackie et al., 1993).Zhadanov，Abubakar，Kruglyakov等人深入研究了采用扩展Born近似、广义扩展Born近似、准线性近似等方法来求解积分方程(Kruglyakov et al., 2016; Zhdanov et al., 2000; Abubakar and Habashy, 2005)，从而大大加快了积分方程的计算速度，但是降低了求解的精度，难以处理电导率高对比度问题.上述研究都是基于规则的六面体单元进行网格剖分，随之计算机技术的不断发展和网格剖分技术的推进，结构化网格在模拟地下复杂异常体带来的离散误差的不足日益突出.更为重要的是，上述大部分研究通常采用等效体积法进行并矢Green函数中存在奇异性的积分，进一步降低了积分方程的计算精度.

本文采用一套新的三维大地电磁积分方程求解策略.采用四面体单元离散地下异常体，减少网格离散带来的误差，提高积分方程法处理地下复杂异常体大地电磁正演问题的能力.采用一种无奇异性的解析表达式计算并矢Green函数奇异值积分，解决并矢Green函数强奇异性问题，提高已有积分方程法的计算精度.最后，基于国际标准测试模型3D-1，将本文计算结果与面向目标的自适应有限单元法、T-Ω有限单元法、E-φ有限单元法和已有积分方程法(COMMEMI)五种结果进行对比，验证本文方法的精度.对不同电导率对比度模型进行模拟，将结果与面向目标自适应有限单元法进行比较，进一步验证本文方法的精度，同时说明本文方法能够处理高对比度问题.以棱柱体模型为例，验证本文方法的收敛性，并说明本文方法适合处理地下复杂模型.

2 积分方程理论 2.1 正演理论

对于大地电磁正演问题，忽略位移电流，电磁场可以由如下形式的Maxwell方程组表示(时间因子为e^iωt)(Nabighian, 1987)：

(1)

(2)

其中，E为电场强度，单位为V·m^-1；H为磁场强度，单位为A·m^-1；μ为磁导率，单位为H·m^-1；σ为电导率，单位为S·m^-1；ω=2πf为角频率，f为频率，单位为为虚数单位.

将电磁场(E、H)分为背景场(E^P、H^P)和二次场(E^S、H^S)，背景场可以表示为

(3)

(4)

其中，μ₀为均匀半空间的介电常数，σ₁为背景模型的电导率.通过式(1)—(3)，式(3)—(4)，并整理得

(5)

(6)

其中，J^S=(σ₂－σ₁)E，为散射电流，只存在于异常体中，σ₂为任意形状的异常体的电导率.

二次场可以表示为(Ting and Hohmann, 1981)

(7)

其中，A和φ分别为矢量位和标量位，可以表示为

(8)

(9)

其中，G为标量格林函数，可以表示为

(10)

其中，R=|r－r′|，表示点r=(x, y, z)到点r′=(x′, y′, z′)的距离.对于均匀半空间中存在异常体，需要加上额外项来考虑空气中的电流项.由式(7)-式(9)可知，二次场由电流和电荷产生，电荷出现在散射电流J^S不连续界面.加入入射场和二次场，可以得到方程：

(11)

式(11)可以写为

(12)

其中，G为并矢Green函数.需要指出的是，背景模型为均匀半空间时，G为均匀半空间并矢Green函数(Ting and Hohmann, 1981)，背景模型为层状介质时，G为层状介质并矢Green函数(Wannamaker et al., 1984).

为了进行数值求解，采用非结构化网格(Tetgen(Si, 2015))将异常体离散成N个四面体单元，假设每个单元内部的电导率和电场为一常数，那么式(12)可以写成

(13)

其中，Γ为并矢Green函数的体积分，需要指出的是，Γ在r=r′时存在奇异性，需要特殊处理.那么，每个四面体单元m中点处的电场可以表示为

(14)

重新排序，得

(15)

其中，

(16)

其中，I为单位张量，O为零张量.将式(15)写成矩阵的形式，有

(17)

通过式(17)，可以计算得到异常体内的电场值，单元矩阵可以表示为

(18)

2.2 并矢Green函数计算

计算L_mn时，需要计算并矢Green函数的体积分.当r=r′，G(r, r′)和∇∇′G(r, r′)存在奇异性，因此，积分∫_vG(r, r′)dv′和∫_v∇∇′G(r, r′)dv′存在奇异性，传统积分方程法采用等效体积法进行近似计算，降低了积分的计算精度，本文基于前期在重力无奇异解析计算的工作基础，利用散度定理转换强奇异值体积分为一系列弱奇异性的面积分公式，借助弱奇异性的面积分公式的解析计算表达式(Ren et al., 2012, 2017a, 2017b)，从而高精度计算并矢Green函数中存在奇异性的积分.首先，令

(19)

那么积分∫_vG(r, r′)dv′可以写为

(20)

对于式(20)中的第一项，由于积分核函数一阶连续可导，可以通过高斯数值积分高精度计算，第二项由于核函数在r=r′时存在奇异性，本文使用散度定理将体积分转换为面积分.忽略系数，将其写成统一的形式∫_vR^qdv′，其中q=±1.结合恒等式(Ren et al., 2017c)：

(21)

那么，积分可以写为

(22)

其中，n_s为四面体各个面的法向方向，h_j=n_sj·(r-r′)表示点r到第j个面的距离.需要注意的是，h_j可以为负值.那么只需要考虑如何解析计算K^q=∫_SjR^qds′.

同样地，积分∫_v∇∇′G(r, r′)dv′可以写为

(23)

式(23)中第一项的积分核函数一阶连续可导，可以通过高斯数值积分高精度计算.第二项的积分核函数在r=r′时存在奇异性，本文借助其解析计算表达式来高精度计算.忽略系数，将其写成统一的形式，其中q=±1.积分和的解析计算表达式见附录A.借助这些弱奇异性的面积分公式的解析计算表达式，从而精确计算了L系统矩阵.采用Pardiso求解器(Kalinkin et al., 2015)求解式便得到异常体内电场分布.

当获得异常体内的电场值后，任意一点的电场值只需要通过式求解计算.任意一点处的磁场通过式得到

(24)

其中，H^P为电导率为σ₁的均匀半空间的磁场.利用观察站的电场值和磁场值，获得测点处的视电阻率和相位.

3 测试结果 3.1 正确性及高对比测试

测试平台为：Intel(R) Core(TM) i5 CPU 3.30 GHz，8G内存.为验证程序的正确性，采用COMMEMI 3D-1标准测试模型(Zhdanov et al., 1997)进行验证，模型如图 1所示.均匀半空间的电阻率为ρ₁=100 Ωm，低阻异常体的电阻率为ρ₂=0.5 Ωm.异常体埋藏在地下x=[－0.5 km, 0.5 km]、y=[－1 km, 1 km]和z=[0.25 km, 2.25 km]处.边界处的电场变化较大，异常体边界附近采用密网格进行剖分，内部采用相对较粗的网格进行剖分，剖分成2401个四面体单元.计算频率为0.1 Hz.61个测点均匀分布在地表(x=[－3000 m, 3000 m]).

将本文方法计算的视电阻率和相位结果与其他4种方法的计算结果进行对比，如面向目标的自适应有限单元法(Ren et al., 2013)、T-Ω有限单元法(Mitsuhata and Uchida, 2004)、E-φ有限单元法(Farquharson and Miensopust, 2011)和已有积分方程法(COMMEMI)(Zhdanov et al., 1997).由于计算结果沿着剖面是对称的，所以只需要对比半条剖面的计算结果.面向目标的自适应有限单元的结果是基于网格依次迭代产生，因此具有非常高的精度，可被用于做参考值.对比图 2可知，本文结果和其他几种方法的计算十分吻合，与面向目标的自适应有限单元法计算结果非常一致，视电阻率相对误差小于3.1%，相位绝对误差小于0.1°.而传统COMMEMI积分结果则和面向目标的自适应有限单元法结果具有较大偏差.因此，本文的积分方程法计算结果对于传统积分方程法计算结果存在较大的改善.改善的主要原因是本文采用了精确的并矢量Green函数奇异值积分解析解.本文只使用了7203个未知数，便得到了高精度的计算结果，突出了积分方程法只需要对异常体剖分、未知数少的特点.但是，需要指出的是，由于积分方程法合成的是密实方程组，求解方程组比较耗时.

电导率高对比度模型是积分方程法的固有缺陷.采用如图 1所示模型，测试积分方程对不同对比度模型的模拟能力，均匀半空间电阻率为ρ₁=100 Ωm，低阻异常体电阻率分别为ρ₂=20 Ωm、0.5 Ωm、0.1 Ωm，对应的对比度分别为50:1、200:1、1000:1.

图 3中，离散点为本文积分方程正演计算结果，实线为采用面向目标的自适应有限单元法结果.两种方法的计算结果十分吻合，视电阻率相对误差和相位绝对误差如表 1所示，对比度为1000:1的正演模型，ρ_a^xy的相对误差小于3.6%，ρ_a^yx的相对误差小于4.8%，φ^xy绝对误差小于0.11°，φ^yx的绝对误差小于0.12°，表明本文方法能够处理电导率高达1000:1模型问题.图 3还可以看出，随着异常体电阻率的减小，视电阻率的最小值变小.对比度达到1000:1时的ρ_a^xy与对比度为200:1差异不是特别大.相对于ρ_a^xy，ρ_a^yx更能够接近异常体的真实电阻率，但是，由于体积效应，二者的视电阻率值都大于异常体的真实电阻率.

3.2 收敛性测试

为了体现模拟复杂模型的能力和测试收敛性，建立了具有复杂几何结构的六棱柱低阻模型，如图 4所示.均匀半空间电阻率为ρ₁=100 Ωm，六棱柱电阻率为ρ₂=1 Ωm.计算频率为1 Hz，测点个数为41，均匀分布在地表(x=[－2000 m, 2000 m]).

对六棱柱区域进行局部网格加密，获得单元个数为834、1541、2667和3840个的四种从疏到密的网格，如表 2及图 5所示.

图 6展示了积分方程法采用上述四种网格剖分计算得到的视电阻率和相位结果.以最密的网格(3840个四面体)为参照，稀疏网格剖分正演结果视电阻率相对误差和相位绝对误差如表 2所示，最稀的网格(834个四面体)的ρ_a^xy的相对误差小于0.05%，ρ_a^yx的相对误差小于0.035%，φ^xy绝对误差小于0.075°，φ^yx的绝对误差小于0.005°.总体来说，四种网格结果十分接近，可以说明本文提出的积分方程法具有很好的收敛性，同时说明，积分方程法具有非常快的收敛性，当结果达到一定精度时，加密网格不会有较大幅度精度的提高.

4 结论

(1) 本文采用四面体单元进行网格剖分，实现了地下复杂异常体的三维大地电磁正演模拟，借助奇异值积分的解析计算方法，给出了基于四面体单元的并矢Green函数无奇异性高精度计算方法，提高了传统积分方程的计算精度，并通过算例验证了本文方法的正确性和适合处理复杂模型的能力.

(2) 本文的研究成果为地下复杂模型的大地电磁积分方程快速、高精度正演研究提供了基础.同时，本文采用的奇异值积分解析计算方法同样适用于直流电法、可控源电磁法等的积分计算.

附录A

下面给出积分和K^q=∫_Sj∇R^qds′的解析计算表达式(Yla-Oijala and Taskinen, 2003).首先，有如下已知积分表达式：

(A1)

(A2)

其中，∂_iS, i=1, 2, 3表示面S的三条边，R_i⁺和R_i^－分别表示点r到∂_iS两个端点的距离，s_i^－和s_i⁺分别表示点r在∂_iS的两个端点的垂向距离，K^－3为式K^q=∫_SjR^qds′在q=－3时的一个特殊情况.β_i的计算表达式如下：

(A3)

其中，|ω₀|为点r到三角面S的距离，需要注意的是，ω₀可以为负值.

符号说明如图A1所示.以四面体的一个三角面S为例，三角形的面积为T，面外方向方向为n.建立局部坐标系，以面S的一个顶点为坐标原点，以改点所在的一条边为坐标轴，建立局部坐标系uov，点d为点r在面S上的投影，三角面S的三个顶点q₁、q₂、q₃可以表示为(0, 0)、(l₃, 0)和(u₃, v₃)，每边的长为l_i, i=1, 2, 3.每条边的外方向方向为m_i, i=1, 2, 3，每边的切向方向为t_i, i=1, 2, 3，有n=m_i×t_i, i=1, 2, 3.点r到三条边的距离表示为t_i⁰, i=1, 2, 3.具体计算表达式如下：

借助如下迭代计算公式：

(A4)

(A5)

那么，积分K^q=∫_SjR^qds′和K^q=∫_Sj∇R^qds′可以表示为

(A6)

(A7)