空间-波数域三维大地电磁场积分方程法数值模拟

2022-06-02戴世坤陈轻蕊凌嘉宣李昆

地球物理学报 2022年6期

戴世坤, 陈轻蕊*, 凌嘉宣, 李昆

1 中南大学地球科学与信息物理学院, 长沙 410083 2 中南大学有色金属成矿预测与地质环境监测教育部重点实验室, 长沙 410083 3 西南石油大学地球科学与技术学院, 成都 610500

0 引言

电磁法勘探广泛应用于探测地壳物质结构、普查石油天然气、煤田、地热和寻找地下水和金属矿产等资源中.频率域电磁法可分为天然场源和人工源法.电磁法正演方法主要有有限差分法、有限单元法、有限体积法和积分方程法四种.有限差分法原理简单，是以差分和差商代替求导的一种数值方法，目前应用得比较多的是Yee(1966)创立的交错网格，能够很好的处理场在介质内部的不连续性，在求解大地电磁各向异性介质中的场(殷长春等，2014；薛帅等，2017)、带地形(Mackie et al.，1993，1994；陈伯舫等，1998；董浩等，2014)和算法并行(Tan et al.，2006；李焱等，2012；秦策等，2017)等方面都有突破， Varilsüha和Candansayar(2018)对比了基于不同规范下大地电磁控制方程的有限差分正演精度和速度，总结了direct EM, the Coulomb-gauged, the ungauged, the Lorenz, 和the axial-gauged formulations这四种方法在高频和低频迭代计算时的一些规律.罗威等(2019)基于交错网格有限差分法研究了球坐标系下的大地电磁三维正演；有限单元法是依据变分原理或加权余量法推导出的与定解问题相等价的积分弱解形式，是求解边值问题中数学理论最为完备的数值方法，在地球物理中应用的较多的可分为节点有限元和矢量有限元法.节点有限元(Zyserman and Santos，2000；Mitsuhata and Uchida，2004；Farquharson and Miensopust，2011；冯建新等，2012；Grayver and Bürg，2014)不满足电场法向分量不连续性，需要进行散度校正.矢量有限元可以确保不同物性边界处切线方向场的连续性，且自动满足无源区散度为零的条件，不需要散度校正，可以克服传统节点有限元出现伪解的不足，由于这个特性，矢量有限元(Liu et al.，2008；顾观文等，2014；Ren et al.，2014；Kordy et al.，2016；Prihantoro et al.，2016)在大地电磁数值模拟中得到了很好的应用，已逐渐代替节点有限元，成为地球物理电磁场数值模拟的标准方式(汤井田等，2015).有限体积法又称控制体积法，先将整个计算区域离散成一系列不重叠的控制网格，然后将微分方程在每一个控制体积进行体积分，单元合成得到线性方程组.有限体积法(Haber and Ascher, 2000; Haber and Ruthotto, 2014；Weiss，2013；Jahandari and Farquharson，2014；陈辉等，2018)是有限差分和有限单元法的一种中间产物，近十几年才发展比较迅速，相对于有限单元法，计算精度略低但效率更高.上述三种算法需要对整个计算区域进行剖分，对于大规模使得形成的矩阵方程大、要求解的未知量个数多，对计算机的性能要求高.

相比之下，积分方程法有如下优点：(1)只需对异常体进行剖分，计算量小，占用内存少；(2)积分方程的解析解具有半解析解的精度，常常被用于测试新开发算法的精度；(3)基于积分方程的反演算法精度和效率高(任政勇等，2017)，所以开发高精度高效率的积分方程正演算法仍具有一定的研究价值.Raiche(1974)、Hohmann(1971，1975)及Wannamaker和Hohmann(1984)为了积分方程技术在三维电磁模拟中的应用做了很多基础性的工作；Berdichevsky和Zhdanov(1984)介绍了大地电磁场谱域满足的方程和表达式；Wannamaker(1991)对大地电磁积分方程法的精度和效率进行了初探；陈久平等(1990)，殷长春和朴化荣(1994)分别对半空间、两层及多层大地介质中的3D异常体的电磁场积分方程法正演进行了研究.鲍光淑等(1999)进行了均匀半空间频率域三维电磁响应的数值模拟,深入研究了均匀导电半空间频域三维电磁散射问题；Zhdanov和Fang(1997)，Hursán和Zhdanov(2002)使用积分方程法模拟了3D地电结构的电磁响应，改进了格林函数算子；Zhdanov等(2006)提出了一种新的积分方程(IE)方法，用于非均匀背景电导率(IBC)复杂结构的三维电磁(EM)建模.美国犹他大学(2006)推出了一个基于积分方程法的三维正演软件INTEM3D，该软件可以对水平层状介质中的三维地电结构用积分方程法进行频率域电磁模拟；后来学者研究了(Farquharson et al.，2006；Zhdanov et al.，2007)用于大电导率对比模型三维电磁建模的新方法；徐凯军等(2006)提出利用结合连分式展开的高斯求积代替常规的快速汉克尔变换方法，进行了半空间大地电磁正演模拟；张辉等(2006)用直接求解体积分方程的方法模拟了电偶源激发时均匀导电半空间频率域三维电磁响应，张量格林函数采用差分近似的方法求解，结合连分式展开的高斯求积求解含有贝塞尔积分，取得了较好的计算效果；阮百尧等(2007)提出了一种用边界元法计算大地电磁场三维地形影响的数值模拟方法；张罗磊等(2010)基于MNS技术进行了大地电磁三维正演模拟，提出在空间-波数域计算张量格林函数，大大提高了计算效率；陈桂波等(2009)采用积分方程法进行了各向异性地层电磁场三维数值模拟研究，分析了各向异性对三维电磁响应的影响特征和规律；Zaslavsky等(2011)采用有限差分和积分方程混合的方法(称为有限差分积分方程法)，降低了每次迭代的计算成本；任政勇等(2017)采用四面体单元、解析的并矢Green函数奇异积分表达式，借助于PARDISO高性能并行直接求解器，进行了三维大地电磁积分方程正演.

但是上述方法最终都归结于大型线性方程组的求解，三维电磁场空间域数值模拟计算量大，存储要求高，限制了现有方法的大规模应用.针对这一问题，本文提出一种空间-波数域三维电磁场数值模拟方法，该方法利用电磁场积分方程是一个三维卷积的特性，沿水平方向进行二维傅里叶变换，将三维空间域卷积问题转换为多个不同波数的空间垂向一维积分问题，一维积分相对独立，并行度好，由此大大减少了计算量和存储需求；保留垂向为空间域，将一维积分垂向可离散为多个单元积分之和，每个单元采用二次形函数表征散射电流变化，得出单元积分的解析表达式，计算精度高、效率高，垂向单元网格灵活，可以准确模拟任意复杂模型，兼顾计算精度与计算效率；利用压缩算子，将电磁场采用迭代法求解，占用内存小，计算速度快.该方法充分利用快速傅里叶变换的高效性和一维形函数积分的高精度特性，实现了三维电磁场高效高精度数值模拟.与常规积分方程法不同，本文方法适用于地下任意复杂介质，计算效率高；且相较于其他算法，随着计算规模增大，算法优势越明显.

1 三维电磁积分方程

1.1 基本原理

E(r)=Eb(r)+∭vG(r,r′)Δσ(r′)E(r′)dv,

(1)

式中r(x,y,z)表示观测点位置，r′(xs,ys,zs)表示异常体位置，v为异常体体积，E(r)表示r处的总电场，Eb(r)为r处的背景电场，式中G(r,r′)为r′处的源在r处的电场张量格林函数，式中Δσ是异常导电率，Δσ(r′)=σ(r′)-σb(r′)，σ(r′)是异常体的导电率，σb(r′)是背景导电率.

设散射电场的表达式为

(2)

式中Es(r)为散射电场，J(r′)为散射电流.

电场与磁场之间的关系为

(3)

将式(2)进行二维傅里叶变换，可得

(4)

式(4)采用迭代法求解.定义一个线性算子：

G(·)=G(Δσ·E).

(5)

式(1)写为

E=Eb+G(Δσ·E).

(6)

对于式(6)，理论上可以采用迭代法逐次逼近进行求解，即可以写为

E(n+1)=Eb+G(Δσ·E(n)).

(7)

E(n+1)和E(n)分别表示第(n+1)和n次计算得的总场，n≥0，由泛函分析中的Banach定理可知，要使得式(7)收敛的条件是

‖G(Δσ·(E(n+1)-E(n)))‖<κ‖Δσ·(E(n+1)-E(n))‖,

(8)

式中κ<1，‖·‖表示2-范数，使得上式成立的条件是

‖G‖<1.

(9)

基于Singer(1995)、Pankratov等(1997)，Zhdanov和Fang(1997)和Avdeev等(2002)提出的一种波恩级数法， Hursán和Zhdanov(2002)从能量不等式出发，对算子G进行了修正，构造了满足式(9)的压缩算子，该算子能使得迭代稳定收敛(Gao and Torres-Verdin,2006).其迭代格式如下:

E(n)=Eb+G(Δσ·E(n-1)),

(10)

E(n)=αE(n)+βE(n-1),

(11)

式(11)中左端项E(n)是通过压缩算子更新的第n迭代的总场，它将作为第(n+1)次迭代计算的初值，式中α,β是与背景电导率σb、异常体电导率与背景电导率的差Δσ有关的张量:

(12)

(13)

利用水平方向二维傅里叶变换，将散射电流满足的三重卷积转化为多个空间域垂向的一维积分，不同波数之间的积分相互独立，并行性好；垂向一维积分离散为多个单元积分之和，每个单元采用二次形函数表征散射电流，求得单元积分的解析表达式，垂向网格剖分灵活，能兼顾计算精度和计算效率；用迭代求解电磁场，占用内存小，计算速度快.

1.2 波数域格林函数

张量格林函数表示点源(单元偶极子源)的响应，是一个3×3的矩阵，在直角坐标系中表示为

(14)

式中G的第一列表示x方向单位电偶极子源产生电场，第二列表示y方向单位电偶极子源产生的电场，第三列表示z方向单位电偶极子源产生的电场.求解张量格林函数可以转化为求解x、y、z三个方向电偶极子源产生的电场问题.

本文采用的电导率背景模型为均匀层状介质模型，推导均匀层状介质波数域张量格林函数的基本思路是：从Maxwell方程组出发，引入矢量位和标量位，将方程转换成矢量位和标量位的方程，再带入洛伦兹规范条件，得到矢量位所满足的Helmholtz方程.对亥姆霍兹方程进行x、y方向的傅里叶变换，将三维偏微分方程降为一维常微分方程，带入边界条件，得到位的波数域表达式，利用波数域位与场之间的关系，求解波数域场的表达式，分别求解x、y、z三个方向电偶极子源产生的电磁场得波数域表达式，即可合成求得波数域的张量格林函数.推导过程和表达式详见附录A.

根据傅里叶变换的微分性质，将推导空间域张量格林函数的思路应用到波数域张量格林函数的推导，波数域张量格林函数表达式中无复杂的汉克尔积分，将空间域奇异点的计算转化为零波数的计算，方法简单，计算量大大减少，效率高；能大大提升积分方程的正演计算速度.

1.3 一维积分计算

式(4)中，垂向积分为空间域，将zs方向的一维积分离散为多个单元之和，波数域张量格林函数的指数项表达式中存在变量zs，将波数域格林函数中与z，zs无关的系数提取到积分外，剩余积分的表达式可归纳为通式

(15)

2 算法

利用迭代法求电场，计算电场的流程如图1所示，利用如图所示流程求得电场数值后，利用式(3)可求得磁场.

图1 空间波数域算法流程图Fig.1 Flow chart of the space-wavenumber domain algorithm to calculate electric field

图2 模型平面示意图Fig.2 Model schematic plane

3 算例

3.1 算法验证

将平面波作为一次场，进行大地电磁场数值模拟，用美国犹他大学开发的基于积分方程法的三维正演软件INTEM3D的计算结果为参照，验证方法的正确性.测试的计算机为Intel(R) Core(TM) i7-6700HQ CPU主频为2.60 GHz,内存为16 GB、64位win10系统，算法在Microsoft Visual Studio 2015开发平台上运行.

模型XOY平面投影如图2所示，背景为均匀半空间介质，上半空间为空气，导电率σ1=10-12S·m-1，下半空间导电率σb=0.01 S·m-1，将频率分别设为0.01 Hz、1 Hz、100 Hz和10000 Hz，进行大地电磁场三维数值模拟，计算范围x方向-1000～500 m，y方向-1000～500 m，剖分网格节点个数101×101×101，三个方向均匀剖分，Δx、Δy均为10 m，异常体范围x方向-100～100 m，y方向-200～200 m，异常体导电率σ=0.1 S·m-1.基于趋肤深度的考虑，频率为0.01 Hz、1 Hz和100 Hz的z方向计算范围设为0～1000 m,异常体范围z方向200～400 m，Δz为10 m；当频率为10000 Hz时，z方向计算范围设为0～400 m,异常体范围z方向40～120 m，Δz为4 m.傅里叶变换采用4个高斯点的Gauss-FFT法(Wu and Tian, 2014)实现.

设迭代终止的条件为：相邻两次迭代所有节点电场模的总和的相对误差小于10-4.表达式为

(16)

式中|En|表示第n次迭代的总场的模，|En+1|表示第n+1次迭代的总场的模.

图3和图4分别是频率为1 Hz本文算法(space wavenumber domain integral electromagnetic method, SWIEM)的数值解与软件INTEM3D计算的地面视电阻率和相位等值线图，从图中可看出，数值解和传统积分方程解的等值线吻合程度高，ρxx、ρxy、ρyx和ρyy分量的相对均方根误差(Wu, 2016)分别为1.12%、0.052%、0.062%、1.13%，φxx、φxy、φyx和φyy分量的相对均方根误差分别为3.12%、0.0024%、0.0011%、5.13%，ρxx、ρyy、φxx和φyy数量级小，舍入误差稍大，ρxy、φxy、φyx和ρyx误差小于1‰，表明算法正确.

图3 地面视电阻率等值线图Fig.3 The profile contour map of the apparent resistivity on the ground

图4 地面相位等值线图Fig.4 The profile contour map of the phase on the ground

图5为y=0 km不同频率视电阻率和相位对应的相对误差曲线.图中，不同频率视电阻率和相位相对误差均小于0.5%. 通过对比不同频率的电阻率和相位相对误差，进一步表明本文空间波数域电磁三维积分方程数值模拟方法对天然大地电磁场的适应性、稳定性和可靠性.

图5 地面y=0 km不同频率视电阻率和相位的相对误差曲线Fig.5 The relative errors of apparent resistivity and phase for different frequencies in line y=0 km on the ground

3.2 迭代算法收敛性

利用3.1节低频验证模型，分别改变异常体大小(异常体顶面埋深不变，计算频率为10 Hz)、异常体顶面深度(异常体大小不变，计算频率为1 Hz)和计算频率，记录电场三分量达到迭代终止条件所需的迭代次数，结果如表1、表2和表3所示.

表1 异常体大小与迭代次数Table 1 Sizes of anomalies and iteration times

表2 异常体顶面深度与迭代次数Table 2 Depths of anomalies and iteration times

表3 计算频率与迭代次数Table 3 Frequencies and iteration times

综合表1—3可知，分别改变异常体大小、顶面埋深和频率大小，电场三分量达到计算精度的迭代次数相差不大，说明异常体大小、顶面埋深和频率对算法的收敛速度影响较小，几乎不影响算法的迭代收敛性.

而背景导电率与异常体导电率的差异对算法收敛的速度影响大.研究导电率差异对算法收敛速度的影响，保持异常体大小和埋深不变，设计与背景导电率不同差异异常体，记录其迭代次数.利用3.1节模型，设场源为x方向极化，频率为1 Hz，分别设高阻和低阻异常体的导电率为背景导电率5倍，10倍，50倍，100倍和1000倍，分别记录电场三分量达到迭代终止条件所需的迭代次数.

表4和表5分别为不同导电率差异的低阻异常体和高阻异常体的电场三分量达到计算精度的迭代次数.表中，高阻和低阻异常的迭代次数都随着异常导电率与背景导电率差异的增大而增多；高阻达到计算精度的迭代次数比低阻少，收敛更快；因一次场为x方向极化，所以电场Ex分量的收敛慢，迭代次数比其他两个分量多，Ey分量收敛最快；若将一次场改为y方向极化，则Ey分量收敛慢，Ez次之，Ex分量收敛最快.

表4 低阻异常与迭代次数Table 4 Low-resistivity anomalies and iteration times

表5 高阻异常与迭代次数Table 5 High-resistivity anomalies and iteration times

3.3 计算效率

利用3.1节的低频模型，设频率为10 Hz，进行数值模拟正演，记录不同网格剖分个数算法迭代一次的耗时和占用内存，探究算法的计算效率.

表6是不同网格采用标准FFT迭代一次的耗时与内存占用，图6为其相应的变化曲线图，结合图、表发现，随着网格个数的增多，迭代一次所用时间和所占内存近似呈线性增长，当网格数为200×200×100(4080501个节点)时，本文算法串行迭代计算一次的时间约为4 s，正演计算总耗时约50 s，内存占用约1.3 G，计算速度快，占用内存少，在普通的笔记本上也能高效计算.

表6 不同规模串行迭代一次耗时与内存占用Table 6 Time consumption and memory consumption of one iteration for the different grids

图6 不同网格迭代一次所需时间和内存Fig.6 Time and memory required to iterate through different grids

3.4 复杂模型适应性

利用陈辉等(2018)设计的复杂模型，如图7所示，计算范围22.6 km×30.3 km×20.6 km，背景采用均匀半空间介质，上半空间为空气，下半空间背景电阻率为100 Ωm，计算频率为0.01 Hz，三个异常体电阻率分别为ρ1=10 Ωm，ρ2=1 Ωm，ρ3=1000 Ωm.图7为该模型在YOZ剖面和XOY平面示意图.水平x、y方向均匀剖分，z方向非均匀剖分.以x极化方向的电磁场为例，记录不同网格算法迭代一次耗时与计算总耗时，并对比文献(陈辉等，2018)中AGMG-GCR算法的计算效率.表7为本文算法不同网格的计算效率.表8为陈辉等(2018)中AGMG-GCR算法的效率，因本文与文献(陈辉等，2018)中的算法平台和程序运行环境一致，能进行对比.

图7 复杂模型示意图Fig.7 Schematic diagram of the complex model

表7 本文算法计算效率Table 7 Computational efficiency of the SWIEM

表8 AGMG-GCR算法不同网格计算效率Table 8 Computational efficiency of the AGMG-GCR

图8为不同算法计算不同未知量个数时x极化方向电磁场时的耗时，图(a)表示迭代一次的耗时，图(b)表示计算总耗时.由图8可以看出，两种方法的计算耗时随着未知量个数的增加均呈现出线性增长的规律.在相同的算法平台和程序运行环境下，对比AGMG-GCR方法，本文算法效率高，且随着网格个数的增多，本文算法耗时的优势越来约明显，待求解未知量约为4000000个时，AGMG-GCR算法迭代一次耗时约3 s，本文串行算法约为0.5 s，速度快5倍；待求解未知量约为6000000个时，AGMG-GCR方法计算x极化方向电场总耗时约1000 s，本文串行算法耗时约100 s，速度快9倍.相比陈辉等(2018)中计算效率最高的AGMG-GCR方法，本文算法耗时更短，为大规模问题的大地电磁三维数值模拟提供重要的技术保障.

图8 x极化方向电磁场计算耗时(a) 迭代一次耗时; (b) 计算总耗时.Fig.8 Computational efficiency of electromagnetic field in x polarization direction(a) Computation time to iterate once; (b) Computation time to calculate total.

图9为该模型用不同数值模拟方法计算得到的阻抗Zxy和Zyx的实部与虚部的解，测点位于地面y=0 m，x方向-4000～4000 m，共21个测点，每个测点间距400 m.其中FDM表示有限差分法计算的结果，是采用大地电磁三维正反演代码ModEM(Kelbert et al., 2014; Dong and Egbert,2019)计算得到；FEM为大地电磁三维有限元直接解法的计算结果(Xiong et al.,2018)，SWIEM为本文算法的计算结果.从图中可以看出，三种算法计算得到的阻抗实部和虚部吻合得很好.本文算法与有限元结果的最大相对误差为1.1%，与有限差分法结果的最大相对误差为1.5%，误差均在可接受的范围内，表明本文算法能在满足精度要求的前提下达到高效率.

图9 不同数值模拟方法的计算结果对比Fig.9 Comparison of the results of different numerical simulation methods

图10为该复杂模型在地面的视电阻率和相位图，组合异常体电阻率差异大，在异常体位置异常明显，本文算法对复杂模型的适应性强，适用于大规模三维地下任意复杂模型大地电磁快速正演模拟.

图10 地面视电阻率和相位图Fig.10 Apparent resistivity and phase diagram of the ground

4 结论

利用空间域电磁场积分方程为三重卷积的特点，本文提出一种空间-波数域数值模拟方法，该方法借助水平方向二维傅里叶变换，将三维空间域卷积问题转换为多个波数之间相对独立的空间垂向一维积分问题，由此大大减少了计算量和存储需求并且算法高度并行.一维积分垂向可离散为多个单元积分之和，每个单元采用二次形函数表征电流变化，可得出单元积分的解析表达式.保留垂向为空间域，优势之一在于可根据实际情况灵活调整单元疏密程度，兼顾计算精度与计算效率；优势之二在于用形函数拟合求得积分的解析解，计算精度和效率高.该方法充分利用一维形函数积分的高效和高精度、不同波数一维积分之间相对独立高度并行和快速傅里叶变换的高效性, 实现电磁场高效高精度的三维数值模拟.将积分方程软件INTEM3D的数值解与本文算法的数值解对比验证了算法正确；异常体与背景场的导电率差异越大，算法所需迭代次数越多，高阻异常体比低阻异常体收敛速度快；随着网格个数的增多，算法耗时和所占内存近似呈线性增长，与目前主流方法相比，速度快一个数量级以上，且随着计算规模的增大，算法优势越明显. 本文的正演算法为大规模高效、高精度的反演研究奠定了基础.

本文算法在不同波数一维积分计算和不同深度节点电磁场正、反傅里叶变换均具有高度并行性，采用并行计算，计算效率将进一步提升.此外，本文仅研究大地电磁场三维数值模拟，若将背景场改为均匀层状介质可控源电磁场，即可进行可控源电磁场三维数值模拟.

值得注意的是，本文采用的压缩算子在低阻异常体导电率与背景导电率对比度大(>100)时，迭代次数需上百次，影响算法效率；在复杂地质条件下，本文算法需要精细剖分来拟合地下异常体，增加算法的计算量.但由于这两个问题也普遍存在于常规空间域数值模拟方法中，且在满足计算精度的前提下，相比常规算法，本文算法效率高，此时新算法依然有一定的优势.此外，本文算例的傅里叶变换采用标准FFT，网格均匀剖分，下一步将研究非均匀网格条件下的空间-波数域正演，进一步提高算法的适用性和效率.

致谢感谢中国地质大学(武汉)罗天涯博士后提供复杂模型有限单元法的数值解数据，感谢中南大学王旭龙博士生提供复杂模型ModEM软件有限差分法的数值解数据.感谢审稿专家和编辑提出的宝贵意见.

附录A 波数域张量格林函数公式推导

将空间-波数域电场张量格林函数的求解转化为x,y,z三个方向的电偶源产生的场，频率域Maxwell为