基于Lorenz规范的空间-波数域大地电磁三维正演

戴世坤，陈轻蕊*，凌嘉宣，李昆，赵东东，张莹，张钱江

1 中南大学有色金属成矿预测与地质环境监测教育部重点实验室，长沙 410083 2 中南大学地球科学与信息物理学院，长沙 410083 3 西南石油大学，成都 610500 4 桂林电子科技大学电子工程与自动化学院，桂林 541004 5 桂林理工大学，桂林 541004

0 引言

大地电磁法利用天然电磁场进行地下电性结构的勘探，由于其野外工作简单、成本低、受野外地质条件约束小和探测深度广等优点被广泛应用于矿产资源勘查、深部地质构造研究和工程和环境勘查等领域（王家映,1997;Pedersen et al.,2006;闫永利等,2007;Tezkan and Saraev,2008;Becken et al.,2011;徐光晶等,2015）.正演是反演的基础，高效高精度的正演是提高反演速度和效果的关键，也是准确进行实测数据资料处理和地质地球物理解释的重要工具.

空间-波数域算法，是一种结合傅里叶变换算法和空间域方程的方法，有基于积分方程和微分方程两条思路，在重磁位场（李昆等，2019；Dai et al.，2019）、直流电正演（Dai et al.,2021）计算中已成功应用，文献（戴世坤等，2022）基于Coulomb规范的微分方程实现了空间-波数域三维电磁场数值模拟，体现了这种方法的低内存占用、高效高精度的特性.本文基于Lorenz规范实现了空间-波数域三维大地电磁场数值模拟.用二次场方法，利用沿水平方向的二维傅里叶变换，将Lorenz矢量位满足的三维控制方程转换为多个波数之间相互独立的常微分方程组，用有限单元法求解一维常微分方程.引入压缩算子，构成了稳定收敛的迭代格式.设计模型验证了算法的正确性、收敛性和适应性，试验结果表明，新方法具有计算精度高、占用内存少和计算效率高的特性，相比基于Coulomb规范的空间域数值模拟方法（戴世坤等，2022），计算效率进一步提高.

1 理论方法

设时谐因子为eiω t，无源频率域Maxwell方程组可写为

（1）

（2）

（3）

（4）

（5）

将式（5）代入式（1）可得

（6）

由于标量的梯度的旋度为0，引入标量位Φ，可以将E重新写为

（7）

（8）

式（8）即为矢量位和标量位满足的微分方程（Weiss，2013；陈辉等，2016）.

基于二次场方法进行三维大地电磁数值模拟，将场分为背景场和二次场，背景介质设为均匀层状介质，此时，背景场和总场满足的方程可写为

（9）

（10）

（11）

（12）

式（12）即为二次场矢量位和标量位满足的微分方程.

1.1 基于Lorenz规范的空间域矢量位控制方程

引入二次场矢量位与二次场标量位之间的Lorenz规范，表达式为

（13）

式（12）可改写为

（14）

（15）

（16）

（17）

利用电磁场的边界条件可得三个矢量位的边界条件，将式（16）满足的偏微分方程加载相应的边界条件，即可对三维大地电磁场进行求解.

采用传统空间域数值模拟算法求解式（16）的偏微分方程时，将构成一个大型稀疏线性方程组，计算量大，存储要求高，且规模越大，计算效率越低.空间-波数域方法对偏微分方程进行水平方向二维傅里叶变换，将一个三维问题分解成不同波数下相对独立的一维常微分方程问题，一维方程简单且计算量小，能大大减少计算量和内存需求，提高三维数值模拟的效率.

1.2 空间-波数域矢量位控制方程

（18）

上边界

（19）

下边界

（20）

Lorenz规范下，利用电磁场水平分量在边界连续，垂直方向电流密度连续的特点，当背景为突变介质时，突变位置的边界条件如下：

（21）

综合式（18）—（21），即为基于Lorenz规范下的空间-波数域三维电磁场矢量位的边值问题.

1.3 有限单元法

采用基于二次插值的一维有限单元法对式（18）—（21）组成的边值问题进行求解，利用伽辽金方法可转化为有限元方程：

（22）

（23）

1.4 空间-波数域电磁场计算

利用有限单元法求得空间-波数域矢量位后，利用矢量位与电磁场之间的关系式（17）在空间-波数域中的表达式可求得空间-波数域二次电磁场：

（24）

（25）

式（24）和（25）中z方向的求导均采用差分求导，将空间-波数域电磁场进行水平方向二维反傅里叶变换，即可求得到空间域电磁场.

1.5 迭代求解格式

将式（22）右端项中的总场采用背景场替代，得到的解为born近似解，在异常体电导率与背景场电导率差异很小（2倍以内）时，能近似作为真解.当电导率对比度差异大时，born近似解误差大，本文采用迭代法逐次逼近真解，其核心思想是构造一个稳定收敛的迭代算子.在基于积分方程法的三维电磁场数值模拟中，Singer（1995）、Pankratov等（1997），Zhdanov和Fang（1997）和Avdeev等（2002）构建一种收敛的波恩级数，把电磁场积分方程中格林算子G，利用坡印廷定理和能量不等式（Hursán and Zhdanov，2002）构造了格林算子的线性变换，获得了一个范数小于1的修正格林算子‖Gm‖<1.构造了满足任意有损耗介质的电磁场迭代计算稳定收敛的格式（Gao and Torres-Verdin,2006）：

E（n）=αE（n）+βE（n-1），

（26）

式（26）中，左端项E（n）是通过压缩算子更新的第n次正演总场值，右端E（n）和E（n-1）分别表示第n和n-1次正演计算得到的总场，n=1时，E0=Eb，式中α，β的表达式与背景电导率σb、异常体电导率与背景电导率的差Δσ有关，写为

（27）

（28）

若相邻两次迭代的计算结果满足给定的精度要求，则停止迭代，输出最新一次的正演总场值；否则将利用式（26）更新的第n次正演总场值作为第n+1迭代的初始总场，代入式（22）中右端项进行新一轮的正演求解，直到满足精度要求为止.

基于积分法与微分法的统一性，本文沿用积分方程法中构造的迭代格式（27）—（28），测试结果证明，对于任意有损耗的介质，算法的迭代过程稳定收敛.

综上，本文提出的基于Lorenz规范的空间-波数域三维电磁场数值模拟方法的要点主要有以下四点：

（1） 从Maxwell方程组出发，引入Lorenz规范，将电磁场满足的偏微分方程转化为关于Lorenz矢量位的亥姆霍兹方程组，基于二次场的数值模拟方法，得到关于二次场Lorenz矢量位的亥姆霍兹方程组；

（2） 利用水平方向二维傅里叶变换，将二次场矢量位满足的三维空间域偏微分方程转化为多个一维空间-波数域常微分方程组，一维方程组计算量小，求解简单，大大减少了计算量和存储需求；

（3） 常微分方程组采用有限单元法求解，单元采用二次插值形函数，构成三个五对角方程，采用追赶法求解，且不同波数之间的方程组求解相互独立，并行性好，垂向网格剖分灵活，进一步提升了算法的计算精度和计算效率；

（4） 采用迭代法逐次逼近进行求解，借用积分方程法中构造的压缩算子，算法稳定收敛，占用内存小，结合空间-波数域三维电磁场数值模拟算法，进一步提高了算法的计算效率.

2 算法与算例

图1所示为本文算法的流程图.图中ε〈〉为迭代终止的误差计算公式，εmin为设定的最大迭代误差，本文设置的迭代终止条件为：相邻两次迭代所有节点电场模的总和的相对误差小于εmin=10-4，表达式为

（29）

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

傅里叶变换的精度受波数选取、网格剖分的影响，本文模型水平方向网格均匀剖分，波数选取规律满足采样定理.测试表明，综合模型复杂程度和频率趋肤深度的影响选择合适的网格剖分间距，并利用采样定理选取的离散波数进行数值模拟，数值解基本能满足精度要求.具体的波数选取规则详见文献（陈龙伟等，2016），本文不再赘述.

本节设计模型进行大地电磁场数值模拟，验证了算法的正确性，研究了算法的收敛性，并通过对比其他算法分析了算法效率.本文模型采用扩边FFT进行计算，介质磁导率μ均设为真空磁导率μ0，介电常数ε均设为真空介电常数ε0.本文测试的计算机为Intel（R） Core（TM） i9-7980XE CPU 主频为2.60 GHz，内存为64 GB.

2.1 正确性验证

用犹他大学开发的基于积分方程法（Integral Equation algorithm,IE）的三维正演软件INTEM3D的计算结果为参照，验证本文算法（Space-wavenumber domain finite element method,SWFEM）在不同频率和不同异常体电导率情况下视电阻率和相位的计算精度.模型为半空间内的长方体，如图2所示，空气、下半空间介质的电导率分别是10-12S·m-1和0.01 S·m-1，地面z=0，模型计算范围：x方向-500～500 m，y方向-500～500 m，z方向0～500 m.异常体范围：x方向-100～100 m，y方向-100～100 m，z方向50～250 m.剖分网格节点个数101×101×101，三个方向均匀剖分，Δx、Δy均为10 m，Δz为5 m.观测面z=0，测点数据101×101个.

表1 不同频率积分方程法与本文算法的视电阻率和相位均方根误差Table 1 The RRMS of apparent resistivity and phase between SWFEM algorithm and IE algorithm at different frequencies

图1 基于Lorenz规范空间-波数域大地电磁三维正演流程图Fig.1 Flow chart of the 3D MT modeling based on Lorenz gauge in the space-wavenumber domain

图2 长方体模型Fig.2 The sketch of prism model

图3 不同频率积分方程法与本文算法视电阻率和相位对比曲线Fig.3 Apparent resistivity and phase of SWFEM algorithm and IE algorithm at different frequencies

图4 不同频率积分方程法与本文算法视电阻率和相位的相对误差曲线Fig.4 Relative errors of apparent resistivity and phase between SWFEM algorithm and IE algorithm at different frequencies

测试分为两组，第一组测试给定异常体电导率为0.1 S·m-1，频率分别设为0.01 Hz、1 Hz、100 Hz和10000 Hz，地面y=0测线上视电阻率和相位的计算结果及其误差如表1、图3和图4所示.由图表可知，不同频率本文算法与积分方程算法视电阻率和相位数据曲线拟合程度好，相对误差均小于1%，均方根误差均小于0.1%，计算精度高，验证了算法的正确性，且表明算法在不同频率下均能计算正确.

第二组测试给定频率为10 Hz，异常体电导率分别设为0.0001 S·m-1、0.001 S·m-1、0.1 S·m-1和1 S·m-1，地面y=0测线上视电阻率和相位的计算结果及其误差如表2、图5和图6所示.综合图表可知，不同电导率异常体模型本文算法与积分方程算法结果吻合，当异常体电导率为1 S·m-1时，与背景电导率对比度为100，测线上视电阻率最大相对误差小于2%，相位小于0.3%；其他异常体模型两种数值解的相对误差均小于1%.两种数值解地面视电阻率和相位的均方根误差均小于0.15%，表明算法对于不同电导率对比度异常体模型的适应性好.

2.2 收敛性分析

本节研究算法的收敛性，测试了不同电导率异常体迭代计算的收敛速度，记录了电场三个分量的迭代误差随迭代次数的变化.模型为半空间内有一正方体异常体，如图7所示，空气、下半空间电导率分别为10-12S·m-1和0.01 S·m-1，计算频率为10 Hz.模型计算范围：x方向-1000～1000 m，y方向-1000～1000 m，z方向0～1000 m，地面z=0.剖分网格节点个数101×101×101，三个方向均匀剖分，Δx、Δy均为20 m，Δz为10 m.异常体大小400×400×400 m3，范围：x方向-200～200 m，y方向-200～200 m，z方向200～600 m.将异常体电导率分别设为高阻和低阻，与背景电导率对比度均分别为5、10、20、50和100倍，高阻异常体电导率分别为0.0001 S·m-1、0.0002 S·m-1、0.0005 S·m-1、0.001 S·m-1和0.005 S·m-1；低阻异常体电导率分别为0.05 S·m-1、0.1 S·m-1、0.2 S·m-1、0.5 S·m-1和1 S·m-1.进行大地电磁正演，以x极化模式下电场迭代收敛情况为例，研究电导率不同对比度情况下电场三分量的收敛速度.

图6 不同异常体电导率积分方程法与本文算法计算的视电阻率和相位的相对误差曲线Fig.6 Relative errors of apparent resistivity and phase between SWFEM algorithm and IE algorithm at different conductivity of anomalies.

表2 不同电导率异常体积分方程法与本文算法视电阻率和相位均方根误差Table 2 The RRMS of apparent resistivity and phase between SWFEM algorithm and IE algorithm at different conductivity of anomalies

图8和图9分别为不同高阻异常体和低阻异常体电场三分量的迭代误差变化，图中黑色横线为本文设置的迭代终止误差，εmin=10-4.由图中可知，随着迭代的增加，电场三分量的误差均逐渐减小，且与背景电导率差异越大，达到计算精度所需迭代次数越多；相比高阻，低阻电场Ex和Ey分量达到计算精度所需的迭代次数更多，Ez分量在电导率对比度更大时迭代收敛曲线的震荡更剧烈，总体来说，低阻异常体比高阻异常体的收敛慢，但算法均能稳定收敛，表明压缩算子适用于本文算法.

2.3 Dublin模型测试

采用Dublin（DTM1）模型测试方法的效率.均匀半空间中有三个长方体异常，模型计算范围x方向-37.5～37.5 km；y方向-25～25 km；z方向0～55 km，三个方向均采用均匀剖分.半空间地下电导率为σ0=0.01 S·m-1，空气电导率为10-12S·m-1.模型如图10所示，异常体的范围和电导率参数如表3所示.计算频率为0.01 Hz.

文献（戴世坤等，2022）基于Coulomb规范实现了空间-波数域三维电磁场数值模拟，本文采用Lorenz规范空间-波数域方法进行三维电磁场数值模拟.二者采用规范不同导致最后得到的空间-波数域控制方程不同，Coulomb规范有4个未知量（三个矢量位、一个标量位），各未知量满足的方程相互耦合；Lorenz规范有3个未知量（三个矢量位），在各向同性介质中，未知量方程能相对独立.将DTM1 模型采用不同网格剖分，对比Coulomb规范和Lorenz规范空间-波数域算法三维大地电磁正演效率，表4为统计结果，记录了两种算法不同网格剖分情况下单次迭代时间、正演总时间和内存占用，图11和12为相应的变化曲线.

图7 正方体模型Fig.7 The sketch of cube model

图8 不同高阻异常体的电场三分量迭代误差曲线Fig.8 Iterative error curves of three electric fields for high resistance anomaly

图9 不同低阻异常体的电场三分量迭代误差曲线Fig.9 Iterativeerror curves of three electric fields for low resistance anomaly

图10 DTM1 模型示意图Fig.10 Schematic diagram of DTM1 model

表3 DTM1 模型三个长方体异常的分布范围和电导率Table 3 Range and conductivity of three cuboid anomalies in DTM1 model

综合图表可知，两种空间-波数域算法耗时和内存占用随计算节点数增多均近似线性增长，且相同网格剖分情况下，本文基于Lorenz规范下的空间-波数域正演方法比Coulomb规范的耗时更短、内存需求更小.对于872289个节点数模型，本文算法单次迭代时间约0.47 s，基于Coulomb规范算法约1.91 s；本文算法总耗时约240.09 s，基于Coulomb规范算法约962.43 s；本文算法内存占用约458.5 MB，基于Coulomb规范算法588.3 MB，两种算法内存需求差别小，但本文算法耗时减少近3倍.

DTM1模型是一个复杂模型，其计算规模大、异常体范围大，异常体电导率与背景电导率对比度差异大.而电导率差异大会导致本文算法迭代终止时迭代次数较多，影响效率.文献（Long and Farquharson,2019）采用基于径向基函数的无网格单元法进行三维大地电磁数值模拟，计算DTM1模型范围与本文相同，在16核Intel Xeon E5-2670 CPU、主频为2.6 GHz的电脑上，计算120598个离散节点，正演一次耗时约797.7 s.在节点数为222589时，Coulomb规范的空间-波数域算法总耗时312.10 s，本文算法总耗时约69.82 s.节点数比文献（Long and Farquharson,2019）多，在配置更差的电脑上，空间-波数域算法效率更高，且本文算法的优势更明显.

图13为节点数为61×41×89的模型本文算法正演结果与Siripunvaraporn等（2002）的对比图，数据来源于文献（Miensopust et al.,2013），测线y=0 m上视电阻率的相对误差最大为3%，平均相对误差为1.2%，相位的相对误差均在1%以下，满足精度要求，表明本文算法对复杂模型适应性好，且算法效率高.

表4 Dublin模型不同网格剖分计算效率Table 4 Calculation efficiency of Dublin model at different grids

图11 不同网格剖分Coulomb规范和Lorenz规范算法的计算时间Fig.11 Calculation time of Coulomb gauge algorithm and Lorenz gauge algorithm of space-wavenumber domain at different grids

图12 不同网格剖分Coulomb规范和Lorenz规范算法的内存占用Fig.12 Memory of Coulomb gauge algorithm and Lorenz gauge algorithm of space-wavenumber domain at different grids

3 结论

本文提出一种基于Lorenz规范的空间-波数域三维大地电磁数数值模拟方法.引入Lorenz规范，将Maxwell方程转换成矢量位的偏微分方程组.采用二次场方法，将总场拆分为背景场和二次场，处理二次场时，利用水平方向二维傅里叶变换，将二次场矢量位满足的偏微分方程转换为多个波数的常微分方程，常微分方程采用追赶法求解，计算量小、存储需求小、计算速度快.最后利用压缩算子，采用迭代法求解空间域电磁场.取得以下结论：

（1）设计不同频率和不同电导率的异常体验证了算法的正确性和适应性；

（2）分析异常体电导率与背景电导率不同对比度时算法的迭代收敛性，对比度越大，收敛越慢，相同对比度的高阻异常体比低阻异常体的迭代收敛速度更快；

（3）采用DTM1 模型对比了基于Coulomb规范和Lorenz规范的空间-波数域三维大地电磁正演效率，随着计算规模的增大，两种算法耗时与存储均呈近似线性增长，但Lorenz规范的算法耗时更短、内存需求更少.相比传统算法，空间-波数域算法效率高，且本文算法优势更明显，非常适合大规模三维大地电磁正演.

图13 地面y=0 m测线上视电阻率和相位Fig.13 The apparent resistivity and phase in line y=0 m on surface

本文仅研究大地电磁场三维数值模拟，若加上外加源项，即可进行带源的电磁场三维数值模拟，当频率为0时，可以退化为直流电模拟.为大规模三维电磁法数值模拟的快速计算提供了一种高效的新方法，为精细、高效反演成像提供新工具.下一步将研究带地形的处理和各向异性介质的数值模拟.

致谢感谢审稿专家和编辑提出的宝贵意见．

附录A 单元积分形式

式（22）种存在三种形函数积分，查表（徐世浙,1994）可得三种积分的单元积分表达式.

（1）第一类积分

（A1）

（2）第二类积分

（A2）

（3）第三类积分

（A3）

其中，

以上三种类型的积分涵盖了式（22）全部积分类型，其他积分项均可以在这三种积分类型中找到相应的单元积分.