邢雯淋, 曹 辉, 曹礼刚, 李小翠, 魏统彪, 张朝晖

(成都理工大学 地球勘探与信息技术教育部重点实验室，成都 610059)

0 引言

大地电磁测深(Magnetotelluric Sounding，MT)是由20世纪50年代苏联学者Tikhonov[1]和法国学者L. Cagnird[2]提出的以天然的平面电磁波为场源，根据趋肤效应的原理研究地球电性结构的一种地球物理勘探方法。目前已被应用于各种领域，在油气勘探、构造研究、矿产探测、地震预报、地热勘探等方向均有相应的研究[3-5]。

随着目前计算设备和运算软件的高速发展以及数值计算方法的进步，三维正反演发展取得了巨大的进展[6]。但目前大部分关于大地电磁三维正反演的文献主要集中在数值方法的阐述和讨论上，而很少关注如何有效地对数值方法进行编程，仅有几种基于不同数值方法并充分采用了现代计算技术的三维编码，被学术界广泛应用于现场实测数据，如Siripunvaraporn等[7]采用数据空间的方法使矩阵的维数取决于数据集的大小，而不是模型参数的数量；Kelbert等[8]设计了基于电磁反演计算机代码的模块化系统(称为ModEM)，这些代码在效率、精度和处理多尺度结构和复杂地形的能力方面表现出了出色的性能;Geraskin[9]提出，通常三维电磁代码中只有一小部分专门用于密集计算(一般采用Fortran和C，即静态语言)，大部分的琐碎工作使用Python和MATLAB(即动态语言)更加便捷。这说明一般的三维电磁代码的编写主要还是使用 “双重编程语言”来完成，而一个崭新语言的Julia正好可以简化“双重编程语言”带来的复杂性，它是在一定的性能需求下，牺牲一些不那么重要的动态性寻找一个更优的平衡点，可以仅使用一种语言来进行编程开发，因此我们可以大大减少编码期间所需的工作，并使代码更容易理解和维护[10]。

因此，笔者尝试使用一个全新的语言，Julia来实现由有限体积法(MFV)建立的基于各向同性介质的大地电磁测深法三维正演模型，并对其进行数值模拟。

1 Julia语言简介

由于程序员在编写高级而复杂的算法逻辑时习惯使用简单的动态语言(如Python、MATLAB)，但为了追求计算速度，编写完后经常会将需要密集计算的部分转换为静态语言(如C、Fortran)。而编写两种语言的代码需要去考虑两种语言之间的类型转换和内存管理，其复杂性可能更甚于仅使用其中一种语言，为此一种新的高级编程语言Julia被开发出来了[11]。

Julia是一种高级的、高性能的、动态的语言，在科学计算领域出类拔萃，其语法像MATLAB或Python一样简单易懂，并且与传统的动态语言相比Julia的核心语言很小，有着丰富的基础类型，具有可选类型标注和多重派发这两个特性，还有接近C语言的性能。Bezanson等[11]将Julia语言与其他六种语言(C++、Python、MATLAB、Octave、R 和 JavaScript)的速度进行了定量分析比较(表1)，结果表明，Julia语言在各项测试中的速度均较快，因为Julia语言在设计之初就考虑如何利用现有的技术去高效加速动态语言，它有一个与其系统结合在一起的即时编译器(JIT)，提供了比Python，MATLAB等更高效地交互式编程和动态性的同时，也可以达到静态编译语言(如C和Fortran)一般地性能，相当于在一定程度上对“双重编程语言”问题进行了优化[10]。

表1 微基准测试结果(相对于C++的时间)[11]Tab.1 Microbenchmark results (times relative to C++)[11]

正是Julia语言的这些特性和它在数值计算上的专业性，才使其能有效降低编程工作的复杂性，并且拥有宽泛的应用范围，能被笔者运用于实现大地电磁三维正演模拟。

2 大地电磁三维正演算法

2.1 边值问题

大地电磁法遵循麦克斯韦方程，取时间因子为e-iωt，可满足如下方程[3]：

▽×E=iωμH

(1)

▽×H=σ-iωεE

(2)

其中：E为电场；H为磁场；ω=2πf为角频率；f为频率；ε为介电常数；μ为真空中的磁导率；σ是介质的电导率。

将式 (1) 代入式 (2) 可导出电场E所应满足的微分方程为式(3)。

▽×(▽×E)-k2E=0

(3)

式中，k2=iωμσ+ω2εμ≈iωμσ，由于在大地电磁勘探低频时位移电流可忽略不计。

为了确保控制方程有确定解，必须设定所谓的“自然”边界条件，由于大地电磁场的源都来自地球的外部，可以认为初始大地电磁场是平面波场，初始电场E的偏振方向沿X轴(图1)。由于空气中的电磁场受到地形和不均匀体的影响而不均匀的分布，因此在数值模拟中必须包括空气层和地下介质[6]。如图1所示，模拟区域设为Ω，由空气层和地下介质两部分组成，Ω=Ω1∪Ω2。在模拟区域中，设空气层的电导率为10-5S/m，这时空气和地下介质的接触面(地表)变成内部边界。

图1 正演三维模型示意图Fig.1 Schematic diagram of forward 3D model

取狄利克雷边界条件：

E(x,y,z)|∂Ω=g(x,y,z)|∂Ω

(4)

式中：g为边界上的矢量电场;Ω为模拟区域，可以使用一维或二维的大地电磁计算值[12]。联立式(3)和式(4)建立了大地电磁三维正演算法的边值问题。

2.2 矢量有限体积分析

这里使用有限体积法求解式(1)和式(2)的解，为了处理任意地形和复杂电导率结构模型，我们将电场和磁场在Yee[13-14]提出的交错网格上进行离散化 (图2)。

图2 Yee交错网格用于确定电场和磁场位置Fig.2 Yee staggered grid is used to determine the position of electric and magnetic field

参照Haber[15]使用有限体积法对式(1)和式(2)离散化完整详细的处理，可得到如下离散的麦克斯韦方程组的类比：

CURLe+iωb=0

(5)

CURLTMf(μ)b-Me(σ)e=s

(6)

A(σ)e=(CURLTMf(μ)CURL+

iωMe(σ)e=-iωs=q

(7)

式中：Me(σ)是唯一包含电导率的项。Haber[16]给出了各向同性电导率单元的Me(σ)的显式表达式为式(8)。

(8)

式中：矩阵Aej为在j方向上单元中心到边缘的平均变量；v是单元体积的矢量；⊙是点方向的哈达玛积。

2.3 大型线性方程组的求解

在大地电磁的数值模拟中，最终都会得出一个大型稀疏线性方程组，式(7)最终可以表示为式(9)[16]。

A(σ)e=K+iωMe(σ)e=q

(9)

式中，K是由▽×μ-1▽×算子离散化得到的对称半正定矩阵。笔者使用直接求解法通过直接分解系数矩阵来求解该稀疏线性方程组，其得出的结果更具高精度和高稳定性[17]。

近年来，在计算成本方面已经对大型稀疏线性方程组的直接求解技术进行了高度优化，已有开源的MPI分布式并行求解器MUMPS程序包和OpenMP共享内存式并行求解器Pardiso程序包等[18-19]。随着计算机性能的不断升级，目前已经可以使用直接法在个人计算机上建立适度大小的模型进行大地电磁的三维正演模拟[20]。综上所述，笔者选择采用MUMPS求解器，对大地电磁三维正演中大型线性方程组进行求解。

3 Julia语言与大地电磁的结合

近年来，在三维地球物理电磁模拟和实用软件的开发中Julia得到了一定的运用，如Peng等[21]完成了基于非精确高斯牛顿法的海洋可控源电磁三维反演其在Julia中的数值实现；Börner等[22]使用Julia语言完成了瞬变电磁数据的三维反演；Ruthotto等[23]为了解决从噪声和间接测量中估计偏微分方程(PDEs)的参数而求解不适定逆问题所需的大量计算，他们提出了一个用Julia语言编写的开源软件jInv，用于解决具有多种测量值的参数估计问题的并行算法；Han 等[20]用Julia语言实现采用插值方法将电场的不同分量平均到同一位置的一般各向异性三维大地电磁正演的有限体积算法，编写过程非常高效，并产生了具有良好可读性、可维护性和可扩展性的代码。由此可知,Julia语言可高效便捷地应用于三维地球物理电磁模拟。

Julia提供了用于所有基本运算的简单函数，并对模块和复合类型有着良好的支持，这里参照Han 等[20]所使用的代码模块化用于三维大地电磁各向同性建模，其结构如图3所示。其中，线性求解模块使用了第三方的Julia包 MUMPS.jl[24]，提供了MUMPS求解器。

图3 三维大地电磁各向同性正演建模代码的主要结构[20]Fig.3 The main structure of the 3D magnetotelluric isotropic forward modeling code [20]

4 Julia语言正演准确性验证

为检验本文采用Julia语言编写的大地电磁三维正演算法程序模拟结果的准确性，选取由MTNet国际论坛上发布的全球公共模型——Dublin测试模型1 (DTM1)进行正确性验证[25]。

DTM1是一个三维正演模型，图4是在均匀的半空间中由3个各向同性的异常体组合而成，均匀半空间的电阻率为100 Ω·m，各个异常体的电阻率如表2所示。

表2 DTM1模型异常体相关数据Tab.2 Data related to abnormal bodies in DTM1 model

图4 DTM1模型三维示意图Fig.4 3D schematic diagram of DTM1 model(a)xz轴垂直横截面；(b)yz轴垂直横截面；(c)xy轴水平横截面

将模型域离散为Nx×Ny×Nz=50×55×82网格单元(包括7个空气层)，最小单元尺寸为2×1×0.5 km，选用0. 1和10 000 s之间21个对数间隔周期进行计算，并对上述频率在(0，0，0)观测点的正演结果与其他研究者模拟的结果进行分析比对(图5)。

图5 本文Julia代码正演结果(My)与其他研究者正演结果在(0,0,0) 观测点视电阻率和相位的对比Fig.5 Comparison of apparent resistivity and phase at (0,0,0) observation point between forward modeling results of Julia code (My) in this paper and those of other researchers(a)视电阻率_XY模式；(b)相位_XY模式；(c)视电阻率_YX模式;(d)相位_YX模式(e)视电阻率_XX模式；(f)相位_XX模式；(g)视电阻率_YY模式;(h)相位_YY模式

由图5可知，本文方法计算出视电阻率和相位的曲线与其他研究者得出的结果基本相吻合，只有在XX和YY模式时的相位高频段上拟合较差(由于视电阻率数值太小，相位容易受到精度影响)，其他模式的视电阻率和相位无论在低频还是在高频上拟合的效果都很好。

对基于 Julia语言编写的大地电磁正演程序的结果进行精度分析，利用视电阻率和相位误差计算公式：

(10)

(11)

式中：ερ为视电阻率相对误差值；ρJulia为Julia程序模拟视电阻率值；ρOther为其他研究者的模拟视电阻率值；n为频点的个数；εφ为相位相对误差值；φJulia为Julia程序模拟相位值；φOther为其他研究者的模拟视相位值；n为频点的个数。表3为本文Julia程序与其他研究者XY和YX视电阻率和相位的模拟结果的相对误差值。

表3 与其他研究者正演结果在(0,0,0) 观测点视电阻率和相位的相对误差Tab.3 Relative error of apparent resistivity and phase at (0,0,0) observation point between forward modeling results and those of other researchers

由表3可以看出，除了与wsinv3dmt的模拟结果相对误差较大，与其他结果的视电阻率相对误差值均小于7%，相位相对误差值均小于0.02°，符合电法勘探的精度要求。

5 模型试算

为了验证基于 Julia语言编写的大地电磁正演程序的可行性，图6设计了一个高低阻组合模型，其中模型1为2 km×2 km×1.7 km，电阻率为10 Ω·m，埋深为0.3 km；模型2为2 km×1 km×1.7 km，电阻率为1 Ω·m，埋深为0.3 km；模型3为4 km×2 km×1 km，电阻率为1 000 Ω·m，埋深为2 km；背景场电阻率为100 Ω·m，空气层电阻率设为10-8Ω·m。并根据图6模型设计了图7的两种模型与之对比。

图6 高低阻组合模型三维示意图Fig.6 3D schematic diagram of combination model of high and low resistance(a)正视图；(b)俯视图

图7 对比模型Fig.7 Comparison model(a)低阻模型;(b)低阻-高阻模型

先将模型离散为Nx×Ny×Nz= 64×40×47网格单元(包括7个空气层)，然后在地表布置测线范围为X=[-4 000 m，4 000 m] 和Y=[-4 000 m，4 000 m]，测点和测线之间的距离均为 200 m，最后对 0.1 Hz、1 Hz、10 Hz 和 100 Hz 的不同频率的正演结果进行计算，选取结果见图8和图9。对比模型0.1 Hz的结果如图10所示。

图8 XY和YX模式在Y=0测线上各频率的对比图Fig.8 Comparison of each frequencies of XY and YX modes on the survey line Y=0(a)视电阻率_XY模式；(b)相位_XY模式;(c)视电阻率_YX模式;(d)相位_YX模式

图9 XY和YX模式的视电阻率和相位的等值面图Fig.9 Isosurface plots of apparent resistivity and phase for XY and YX modes(a)0.1 Hz_XY_Resistivity/Ω·m;(b)0.1 Hz_XY_Phase/°;(c)0.1 Hz_XY_Resistivity/Ω·m;(d)0.1 Hz_XY_Phase/°;(e)1 Hz_XY_Resistivity/Ω·m;(f)1 Hz_XY_Phase/°;(g)1 Hz_XY_Resistivity/Ω·m;(h)1 Hz_XY_Phase/°;(i)10 Hz_XY_Resistivity/Ω·m;(j)10 Hz_XY_Phase/°;(k)10 Hz_XY_Resistivity/Ω·m;(l)10 Hz_XY_Phase/°(m)100 Hz_XY_Resistivity/Ω·m;(n)100 Hz_XY_Phase/°;(o)100 Hz_XY_Resistivity/Ω·m;(p)100 Hz_XY_Phase/°

图10 对比模型0.1 Hz时XY 和YX模式的视电阻率和相位的等值面图Fig.10 Isosurface of apparent resistivity and phase of XY and YX modes at 0.1 Hz for the comparison model(a)低阻模型;(b)低阻-高阻模型

从图8与图9可以看出，0.1 Hz时XY模式和YX模式均对模型中浅部低阻异常体在水平面上的数值大小和位置显示较准确(其中YX模式的位置显示更为准确)，并且两个低阻体之间较小的电阻差异也能体现出来，且视电阻率值也与异常体实际电阻率较为接近；而对深部高阻体时，XY模式可以较好地显示出水平面上浅部低阻体和深部高阻体重叠的综合反映，而YX模式的显示容易对高阻体的位置产生误判。1 Hz时XY模式和YX模式的显示与0.1 Hz较相近，均对模型浅部低阻体数值大小和位置显示较好，XY模式视电阻率对深部高阻体的数值大小显示较0.1 Hz略差，但也能明显看出有高阻体的存在。

10 Hz时XY和YX模式对模型浅部低阻体数值大小和位置均有适度的显示，但除了XY模式相位图对深部高阻体略有显现，其余视电阻率和相位图均已完全无法显示出深部高阻体的存在；100 Hz时XY和YX模式不仅对模型浅部低阻体数值大小显示较差，并且完全也未能显示出深部高阻体的存在，但对浅部低阻体的位置略有显示。

对比模型在0.1 Hz时XY与YX结果可知，低阻模型与低阻-高阻模型，均对模型中浅部低阻异常体在水平面上的数值大小和位置显示较准确，其中低阻-高阻模型的结果，在深部高阻体处只有XY模式与低阻模型在数值上有所差异，与高低阻组合模型的结果情况基本一致。

综上所述，对本文设计的高低阻组合模型进行正演模拟，不同频率下的正演结果对于浅部低阻异常体的位置反映均较好，但低频时对浅部低阻异常体的模拟结果均较高频时准确，表明频率越低异常体的分辨率越高，且由于低频的勘探深度大，只有在低频时的正演模拟结果才能显示出对深部高阻体的响应。

5 总结

这里通过Julia语言编写的正演代码实现由有限体积法的大地电磁各向同性的三维正演，并进行了正确性验证，得到如下结论：

1)笔者参考全球公共模型DTM1进行正确性验证，并对高低阻组合模型进行正演试算，结果表明，在正确性验证中与多个研究者的正演结果进行对比分析结果拟合较好；在模型试算中各频率对浅层低阻体的正演效果较好，而对深部高阻体只有在低频时才能得到较好的正演结果，这与大地电磁的相关理论相符，并且正演结果基本能将异常体的形态和位置特征体现出来。综上所述两个模型的正演结果,均证明了使用基于Julia语言的大地电磁正演程序是准确、可靠、可行的。

2)对于通常需要进行正反演模拟的研究者和学生来说，Julia作为一种高性能性和高动态性的新型编程语言，提供了跟Python，MATLAB等相比更加高效地交互式编程和动态性，并也达到了一般的静态编译语言(C、Fortran)的性能，它语言结构简单明了容易掌握，并在在数值计算上有着极强的专业性，可以更快更高效更简易的进行我们正反演代码的编写工作。