基于并行化直接解法的频率域可控源电磁三维正演

韩波， 胡祥云， 黄一凡， 彭荣华， 李建慧， 蔡建超

中国地质大学(武汉)地球物理与空间信息学院，地球内部多尺度成像湖北省重点实验室， 武汉 430074



基于并行化直接解法的频率域可控源电磁三维正演

韩波， 胡祥云*， 黄一凡， 彭荣华， 李建慧， 蔡建超

中国地质大学(武汉)地球物理与空间信息学院，地球内部多尺度成像湖北省重点实验室， 武汉 430074

电磁法的三维数值模拟是一个对数值算法和计算机硬件要求都非常高的问题.对常用的微分类方法如有限单元法和有限差分法而言，求解最后所得的大型线性方程组是至关重要的一步，直接影响到正演算法的实用性.如何高效、稳定且准确地解线性方程长期以来一直是被探讨的问题.本文实现了基于线性系统直接求解技术的频率域可控源电磁(CSEM)三维正演.使用交错网格有限体积法(FV)来离散化关于二次电场的Helmholtz方程；使用直接解法取代传统的迭代解法来求解离散线性系统，即对系统矩阵进行完全LU分解，具体通过调用大规模并行矩阵直接求解器(MUMPS)来实现.基于理论模型做了一系列数值实验，首先证明了直接解法的高精度和稳定性，并考察了其内存需求、计算时间和并行可伸缩性等主要计算性能， 最后检验了所开发的算法快速模拟多场源CSEM问题的能力以及对常规海洋和陆地CSEM模拟的有效性.

可控源电磁法； 三维模拟； LU分解； 直接解法

1 引言

频率域可控源电磁法(controlled-source electromagnetic methods, CSEM)通过观测人工场源在地层、空气和海水等介质中激发的电磁场来获取介质的电性分布，长期以来一直是金属矿勘探中最重要的地球物理方法(底青云和王若，2008；Hu et al.，2013)，并在最近十多年逐渐成为海洋油气勘探中一种不可或缺的地球物理技术(Constable，2010).为了通过CSEM观测数据尽可能真实地还原出介质的电性结构以提升勘探效果，实施三维勘探并对数据进行三维模拟是必需的.正演模拟是反演模拟的基础.CSEM中最常用的三维正演方法包括积分方程法(IE)(Avdeev and Knizhnik，2009；王若等，2009；陈桂波等，2009；Zaslavsky et al.，2011)，有限单元法(FE)(徐志锋和吴小平，2010；Schwarzbach et al.，2011；Li et al.，2011；da Silva et al.，2012；Puzyrev et al.，2013)，有限差分法(FD)(Newman and Alumbaugh，1995；沈金松，2003；Streich，2009；邓居智，2011)以及与FD密切相关的有限体积法(FV)(Weiss and Constable，2006；杨波等，2012；Jahandari and Farquharson，2014；韩波等，2015).对这几种方法的评述可以在很多文献中找到(Avdeev，2005；吴桂桔等，2010； Börner，2010).

使用FE、FD和FV法对Maxwell方程离散化后一般会得到大型、稀疏线性系统方程，在三维情况下待求解的未知数很容易达到上百万个，要高效地求解这样的方程并不是一个容易的问题.普遍采取的求解手段是利用Krylov子空间迭代技术(Saad，2003)来逐步逼近方程的真实解.这一类方法不必显式存储系数矩阵的元素，而只需存储矩阵与向量的乘积(Weiss，2001)，因此对计算机内存的需求很小，完全可以在现代的PC机上执行，并且计算速度很快(在收敛性好的情况下).然而，迭代解法的收敛性能很大程度上取决于系数矩阵谱的性质.如果谱的性质不佳，就需要大量的迭代来收敛到满足精度要求的解，但也可能会无法收敛甚至发散，求解失败.此时需要对原始方程进行"预优"处理来改善系数矩阵的谱的性质，再应用迭代解法才有可能快速收敛.在电磁法数值模拟中，空气、地层和海水之间巨大的电导率差异，以及不同区域的离散网格尺寸的巨大差异等往往导致最后的线性系统方程谱性质极差(比如具有很大的谱条件数)，必须要预优.预优因子的选取也是一个问题，好的预优因子往往计算起来较困难，并且没有一种预优因子能保证迭代解法对于任意复杂的模型都能稳定收敛.

基于矩阵分解的直接解法是求解线性方程的另一个途径.与迭代解法相比，直接解法分解矩阵要占用大量的内存，且计算时间一般也更长.尤其是对于系统矩阵极其庞大的3D正演，直接解法的这两个缺点变得非常突出，令普通PC机难以承受，因此在以往基本不被考虑.然而直接解法相对于迭代解法也有显著的优点，首先是求解精度明显更高，且计算时间主要取决于系数矩阵的大小，与其条件数关系不大；此外，求解多个系数矩阵相同而右端项不同的方程组，只需要进行一次矩阵分解，因此计算量与解一个方程组相当(Grayver，2013)，这一点不仅对于多场源的CSEM模拟非常有利(Oldenburg et al.，2008)，而且有利于电磁反演中灵敏度的快速构建(韩波等，2012).近年来，直接解法在2D/2.5D电磁模拟中已有应用(Abubakar et al.，2008；Streich et al.，2011；Luo et al.，2014).随着数值算法的发展，在科学计算领域出现了许多性能优良且免费的矩阵直接求解器软件包，其中包括针对大规模稀疏矩阵作了高度优化的并行求解器，如MUMPS(Amestoy et al.，2006)和PARDISO(Schenk and Gärtner，2004)，加上计算机性能的快速提升，使得在小型服务器、工作站和PC机群上利用直接解法求解3D 电磁问题成为可能.Streich(2009)、da Silva等(2012)以及Jahandari和Farquharson(2014)分别使用FD、FE和FV法进行CSEM三维正演计算，他们均利用MUMPS作为线性求解器；Schwarzbach等(2011)则在其高阶FE算法中调用PARDISO解线性方程.

本文使用交错网格FV法来进行频率域CSEM的三维正演计算，尝试利用大规模并行矩阵求解器MUMPS解线性系统方程.通过大量的理论模型正演测试来考察MUMPS的求解精度、稳定性、内存需求、计算时间以及并行可伸缩性等计算性能，并验证所开发的算法和程序的有效性和广泛适用性.

2 方法理论

2.1 控制方程

对于一般CSEM勘探所使用的频率范围和大部分已知的传播介质的电阻率变化范围，位移电流的影响可以忽略，则在电流型场源激发下关于电场的二阶频率域微分方程为(取时谐因子为eiω t)

(1)

其中，E为电场强度(V/m)，σ为电导率(S/m)，ω为角频率，μ0为真空中的磁导率(H/m)，J为外加场源的电流密度(A·m-2).J为冲激函数，在场源点具有奇异性(无穷大)(WeissandConstable，2006)，使得数值求解方程(1)有一定的困难.为此，可假设有一个参考模型σp，在场源J的激发下同样满足上述形式的方程：

(2)

式(1)与式(2)相减，可得

(3)

若将实际模型的场看做总场，参考模型的场看做一次场，令二次场Es=E-Ep，式(3)可重写为

(4)

其中Js=(σ-σp)Ep，可看做激发二次场的电流密度.方程(4)便是关于二次电场的二阶控制方程，也是本文使用FV法直接求解的方程.与总场方程(1)相比，由于Js不是奇异函数，方程(4)求解起来更加容易.对于一次场，一般将参考模型选为所需计算模型的一维“背景模型”，可以通过解析法求得；一次场与二次场相加便得到总场.

2.2 有限体积法

图1 (a)电磁场的交错采样和(b)采样点处的控制体积Fig.1 (a) Sampling positions of EM fields on a staggered grid and (b) the control volume for a sampling site

求解形如式(4)的Helmholtz方程，一种简单而有效的方法是交错网格有限差分法(Staggered Finite-Difference, SFD)，其将研究区域离散化为一系列规则长方体单元，然后在长方体单元上交错地对电场和磁场采样，如图1a所示.该方法在电磁法3D正演中被广泛使用，其基本理论比较简单，可以在很多文献中找到(e.g.，Newman and Alumbaugh，1995)，这里不再赘述.使用SFD法对(4)进行离散化直接得到的线性方程系数矩阵是不对称的，为了使其对称以便于方程的求解，每个采样点的线性方程两端需要同时乘以一个对应于该采样点的“尺度因子”V，V具有某种体积的意义，即SFD线性方程实际上对应的微分方程为

(5)

对SFD法稍加改动，便可以实现交错网格有限体积法(SFV)，即在每个采样点处定义一个包围该采样点的控制体积，由共享该采样点所在棱边的相邻4个网格单元体积的1/4组成.以x分量为例，如图1b所示.对式(4)在某个体积单元Ω内求积分：

(6)

然后对式(6)求解.事实上，式(5)中的尺度因子刚好是SFV中的控制体积值，因此式(5)与式(6)在物理意义上的区别在于前者用某个离散点处的变量值代表包围该点的一片区域内每一处的变量值，后者则考虑到变量在区域内的变化，对区域进行体积分.式(6)左端第一项可以完全展开(杨波等，2012)，并且与(5)左端第一项具有完全相同的离散形式.对于(6)的左端第二项，通常认为二次电场在控制体积范围内的变化不大，因此Es可提取到积分号以外；采样点处的电导率可与SFD中一样，使用采样点所在棱边的相邻4个单元电导率的体积加权平均值，也不必求积分.经过这样的简化，SFV的左端项就与SFD完全相同，因此二者的线性方程系数矩阵相同，唯一不同的是右端项.由于CSEM的电磁场可能在很小的空间范围内发生较大变化(如场源附近)，因此SFV比SFD更能准确地代表电磁场的分布.

CSEM的场源是局部的，因此可让剖分域的边界离场源足够远，从而运用齐次的Dirichlet边界条件，即令边界上的切向电场为0；SFV法右端项的一次场体积积分可使用数值积分实现，例如Gauss-Legendre积分(DavisandRabinowitz，2007)；一次场的计算为一维正演问题，使用基于Hankel变换的准解析法(Key，2009).

2.3 大型线性系统方程的求解

Helmholtz方程经SFD或SFV离散化后得到线性系统方程：

Ae=b,

(7)

其中A为大型、对称、稀疏复矩阵(仅对角元为复数)，其行(列)数也即待求的未知数个数约为3Nx·Ny·Nz(Nx、Ny和Nz分别为三个方向的网格数)，每行最多只有13个非零元.作为一个简单示例，图2中显示了一个网格数为5×6×7的模型的SFD/SFV离散线性方程系数矩阵的稀疏结构.

本文使用基于矩阵三角分解的直接解法求解方程(7)，具体通过一个开源的软件包——多波前大规模并行求解器(MultifrontalMassivelyParallelSolver，MUMPS)(Amestoyetal.，2006)来实现.MUMPS专门针对大型稀疏矩阵，并且还能利用矩阵的对称性，就内存需求和运行时间作了高度优化，因此非常适合解方程(7)；此外MUMPS还能通过消息传递接口(MPI)实现并行化求解.MUMPS解方程包括3个阶段：①矩阵的分析与预处理；②矩阵的LU(LDLT)分解；③三角方程的求解.其中矩阵的分解占用绝大部分计算时间和内存.

图2 网格数为5×6×7的模型的SFD/SFV离散线性方程系数矩阵A的稀疏结构Fig.2 Sparsity pattern of A arising from the SFD/SFV discretization on a staggered-grid with size 5×6×7

对照控制方程(4)或(6)来看方程(7)，在模型剖分不变的情况下，改变频率，则系数矩阵A与右端项b都会改变，需要重新分解A并求解；若模型剖分和频率均固定，只改变发射源的方位或位置，则一次场会改变，导致b发生改变，但A不会受到影响，因此无需重新进行矩阵分解，而可以利用原有的分解因子，对新的右端项重新进行解三角方程就够了.这样就大大节省了计算量.在实际CSEM勘探中(Oldenburgetal.，2008；Constable，2010；Grayver，2013)，为了更全面地获取地下信息，常需要频繁改变发射源的极化方向(方位)和空间位置，此时利用直接解法可以快速计算多方位、多空间位置发射源的电磁响应.

2.4 计算流程

根据前文的描述，本文的SFV算法程序的计算流程如图3所示，从中可以清晰地看出MUMPS的各个阶段需要执行的次数与频率数、场源数之间的关系.

图3 本文SFV三维正演算法的计算流程Fig.3 Workflow of the presented SFV modeling scheme

3 直接解法的性能测试

为了测试MUMPS求解大型线性系统的性能，基于两个较典型的层状地电模型(图4)进行了一系列正演测试.图4a为海洋CSEM中常用的模拟高阻油气薄层的模型(Constable and Weiss，2006)，图4b为一个陆地三层H型模型.本节的测试主要与模型参数和频率有关，而与一些具体的收发设置如测点设置无关，因此将不提及这些设置.本文的所有计算均在中科曙光“天阔W580I-G10”工作站上完成，该工作站包含两颗Intel XEON E5-2603型CPU，每颗CPU为4核心4线程，主频1.8 GHz.为了应对矩阵分解的巨大内存需求，配置了128 GB的内存.操作系统为Ubuntu 12.04.

3.1 解方程的精度与稳定性

从基本原理上看，直接解法未作任何近似，因此其求解结果在理论上是精确的，但实际执行时由于计算机的字长有限，对浮点数运算会有舍入误差，导致直接解法也不可避免地存在误差，但这种误差与迭代解法中由于对数学问题本身的近似带来的误差有本质区别.

图4 用来测试直接解法性能的(a)海洋与(b)陆地一维层状模型Fig.4 The (a) marine layered 1D model and (b) land layered 1D model used to test the performance of the direct solver

方程(7)中，系数矩阵A的数值特性与频率、空气电导率等密切相关，因此方程的求解也可能会受这些因素的影响，尤其是迭代解法往往对这些因素非常敏感.为此，对以上两个1D模型正演时，大范围改变频率和空气电导率来测试直接解法的稳定性.令空气电导率从10-5S/m变化到10-11S/m，频率在实际勘探常用的范围内变化，表1和表2给出了直接解法解方程的残差和矩阵分解耗时.这里残差的定义式为‖Ae-b‖，相对残差为‖Ae-b‖/‖b‖；两个模型的网格剖分数分别为68×40×43和64×42×46；使用了8个线程并行计算.可以看出，空气电导率的变化对相对残差的影响非常小；频率由高到低变化，相对残差会有明显的增加，但最大的相对残差也仅仅在10-11数量级，而迭代解法收敛的相对残差门槛值一般设在10-5～10-8之间(Grayver，2013).在矩阵分解耗时方面，频率和空气电导率变化对海洋1D模型的影响很小，对陆地模型的影响稍大，但时间的波动相对于整体耗时并不大.由此可见，在实际勘探常用的频率范围内，对于包含很大电性差异的模型，直接解法都具有相当高的精度，计算速度也很稳定.

表1 对海洋1D模型，频率与空气层电导率在较大范围内变化时直接解法的精度和矩阵分解时间Table 1 The accuracy and factorization time of the direct solver for the marine 1D model for a wide range of frequencies & air conductivities

3.2 内存需求与计算时间

再来考察直接解法的内存需求与计算时间对模型剖分网格数的敏感程度.对海洋1D模型，发射频率为1 Hz，空气电导率取10-11S/m(本文从此处开始，空气电导率值均取该值).设计了从小到大21种网格，其中最小的网格数为7×8×8，最大的网格数为98×80×57.图5显示了直接解法的内存需求与计算时间(MUMPS的三个阶段之和)随未知数个数(约为3Nx·Ny·Nz)的变化.在调用MUMPS时执行了“Out-Of-Core”选项，即利用硬盘来存储矩阵的分解因子，这样能节约很大一部分内存.对于较大的网格数，使用了多个进程并行计算(进程的增加会导致内存需求的增加，见下一小节).可以看到，内存需求与计算时间均随模型网格数的增加而快速增加.当网格数为78×50×57时(对应的未知数个数约为650000)，内存需求在16GB左右，而当网格数达到98×80×57时，内存需求则达到了接近120GB.虽然无法给出网格数与内存需求之间的具体关系式，但不难看出整体上内存需求增加的速度要远快于网格数增加的速度.

表2 对陆地1D模型，频率与空气层电导率在较大范围内变化时直接解法的精度和矩阵分解时间Table 2 The accuracy and factorization time of the direct solver for the land 1D model for a wide range of frequencies & air conductivities

图5 对海洋1D模型，对于不同的模型网格剖分数直接解法的(a)内存需求与(b)解方程时间Fig.5 (a) Memory requirement and (b) solving time of the direct solver for the marine 1D model for various grid sizes

3.3 并行可伸缩性

在并行计算中，可伸缩性(Scalability)是指系统通过改变可用计算资源和调度方式来动态调整自身计算性能的能力.为了测试MUMPS的可伸缩性，依次使用1～8个进程对海洋1D模型进行正演计算，频率为1 Hz，使用3种网格剖分，统计了解方程的时间、相对于单进程运算的加速因子、总内存消耗和平均每个进程的内存消耗等，如图6所示.可以看到，计算时间随进程数的增加稳定下降，当进程数达到8时，计算时间是单进程的1/4左右.但从加速因子来看，还远不能达到理论上最理想的并行(加速因子等于进程数)，这是由矩阵分解的可并行程度有限所决定的.对于同一种网格，总体内存消耗随进程数的增加呈上升趋势(个别情况除外)，8进程运行的总内存消耗是单进程的4倍左右，这是由于进程之间需要共享很多变量.因此在试图增加进程数来提升计算速度时，还要考虑到进程增多带来的内存需求增大.尽管如此，内存消耗平均到每个进程上，还是会随进程数增加而下降的，这一点对机群系统比较有利.

4 三维正演程序的有效性检验

接下来通过具体正演结果的比较与分析来检验所开发的SFV三维正演程序的有效性.

4.1 与解析解的对比

首先，还是基于第3节中的两个层状模型，比较SFV数值解与解析解.所有的解析解均使用开源的Dipole1D程序(Key，2009)获得.

对于海洋层状模型，观测系统设置为：发射电偶极沿x方向，位于海底上方50 m，取偶极中心在x-y平面内的投影为坐标原点，发射频率为1 Hz；接收器位于海底表面，沿x轴方向布设(即inline观测，参考图4a).图7给出了网格数为86×50×57时Ex分量和Hy分量的振幅-收发距(MVO)曲线与相位-收发距(PVO)曲线的SFV解与解析解对比，图8则为在3种不同网格剖分下SFV解相对于解析解的误差.可以看出，SFV解的误差整体上随收发距的增加而增大；网格为54×32×27时误差非常明显，Hy分量振幅的最大相对误差可达20%以上；当网格加密为68×40×43时，数值解的精度得到了极大的提高，振幅最大相对误差在8%左右，相位最大相差6°左右；进一步细化网格至86×50×57时，精度继续得以提升，但并不显著，这也许跟网格加密的方式有一定关系(不在本文讨论范围以内).

图6 对海洋层状模型，直接解法解方程的时间与内存需求随所使用的进程数的改变Fig.6 Solving time & memory requirement of the direct solver vary as the number of computing processors changes for the marine layered model

图7 对海洋层状模型，网格为86×50×57时水平(a, b)电场和(c, d)磁场的SFV数值解(圆圈)与解析解(线条)对比Fig.7 Comparison of horizontal (a, b) electric fields and (c, d) magnetic fields obtained from the SFV method (circles) and the analytic solution (line) for the marine 1D model with grid size 86×50×57

图8 对于海洋层状模型，在3种不同网格剖分下SFV解相对于解析解的误差Fig.8 Relative errors of SFV solutions compared to analytic solutions for the marine 1D model for three different grid sizes

图9 对陆地层状模型，坐标为(8000, 0)的测点的水平(a, b)电场、(c, d)磁场和(e, f)卡尼亚视电阻率的SFV解(圆圈)与解析解(线条)对比，以及(g, h)视电阻率的SFV解相对于解析解的误差Fig.9 Comparison of (a, b) electric fields, (c, d) magnetic fields and (e, f) Cagniard apparent resistivities of the site with coordinate (8000, 0) obtained from the SFV method (circles) and the analytic solution (line) for the land 1D model. (g, h) are the apparent resistivity errors of SFV solutions compared to analytic solutions

对于陆地层状模型，采用地面可控源音频大地电磁法(CSAMT)(底青云和王若，2008)的观测方式，观测系统设置为：发射源为沿y方向的1000 m的接地导线，取其中心为坐标原点，发射0.125～8192 Hz之间按对数等间隔分布的17个频率；测线平行于发射源布设，中点到场源中心的最大距离为8 km(参考图4b).网格数为66×80×46.图9给出了坐标为(8000, 0)的测点上所有频率的Ey分量、Hx分量和卡尼亚视电阻率ρyx的SFV解与解析解对比，以及ρyx的误差.其中ρyx的计算式为

(8)

可以看出，数值解与解析解吻合较好.频率位于64～256 Hz之间时误差稍大，其余频率的视电阻率相对误差均在8%以内，阻抗相位误差均在±2°以内.

4.2 多场源问题的模拟

考虑海洋CSEM勘探中需要发射源在多个位置激发的情况，对上述海洋1D模型，令发射源沿x方向在-2000 m到8000 m之间移动，其余观测参数与4.1节中相同.如2.3节中描述的，频率不变时，发射源每改变一次激发的位置(或方位)，就会给三维正演模拟中线性系统方程增加一个右端项，但不会增加系数矩阵.图10为右端项数目由1增加到101时，SFV正演计算的总时间及各主要部分消耗的时间变化.可以看到，当右端项比较少时，矩阵分解占用了绝大部分计算时间，而解三角方程(Solution)这一步消耗的时间几乎可以忽略不计；随着右端项的增多，由于SFV法每构建一次右端项都要在控制体积内进行数值积分，构建右端项所占用的时间比例越来越大；当右端项增加到101个时，构建右端项占用的时间比例已经很显著了，但解三角方程所占比重仍然很小，总计算时间也仅仅只比1个右端项时增加了约1倍.

4.3 海洋三维模型

若将前面的海洋层状模型中高阻层在水平方向的边界截断，则得到更接近真实情况的三维油气储层模型，如图11所示.已有不少其他作者研究过类似的模型(e.g.，Constable and Weiss，2006；Streich，2009)，本文的模型与这些作者略有差别，即高阻体在x-y平面内的截面为正方形而非圆形.

分别考虑海水无穷深、不含空气层和海水1000 m深、包含空气层这两种情况.参数设置为：场源沿x方向，位于高阻体x负向边界中点的正上方，海底上方50 m，坐标为(0, 0, -50 m)；发射频率为1 Hz；接收器沿x方向布设、分布于海底y=0的测线上(inline观测).对高阻体截面边长(设为a)为1 km、2 km和5 km的三种情况进行了正演.所得到的水平电场Ex分量如图12所示，其中还给出了不含高阻体的半空间(a=0)和高阻体在水平方向无限延伸(a=∞)这两种极限情况的响应(使用解析法获得).

看图12a和12b，Ex响应的分段特征非常明显：在小收发距内，经海水传播的直达波占主导地位，振幅大约按1/r3衰减；随着收发距的增大，地层波开始起主导作用，含高阻层(a>0)与不含高阻层(a=0)的响应逐渐区分开来，振幅曲线在对数坐标系中的斜率趋于常数，说明振幅趋于按指数衰减.a=1 km与半空间模型的响应几乎重合，说明该尺寸的异常体难以分辨；a=2 km与半空间的响应虽然在大收发距时的衰减速度保持一致，但二者的绝对差值较为明显，说明该尺寸的异常体容易分辨；a=5 km与半空间的响应差别非常明显，且在收发距大于高阻体边长时与高阻体无限延伸(a=∞)的响应区分开来，衰减速度逐渐趋于半空间模型的衰减速度.根据以上分析，可以再次确认Constable和Weiss

图10 对于海洋层状模型，SFV正演总体计算时间及各主要阶段消耗的时间随右端项数目的变化Fig.10 The total running time and running time of the significant tasks in SFV modeling vs the numbers of right hand sides for the marine 1D model

(2006)的结论：对于此类水平薄板状高阻体，只有当其水平尺寸至少为埋深的2倍时，才可以利用CSEM将其分辨出来.

再来看包含空气层的情况.对比图12c、d和a、b，可以看到，空气波的影响非常明显，在收发距超过5 km时，a=1 km、a=2 km和半空间模型的响应中空气波开始起主导作用，电场振幅衰减极其缓慢；而a=5 km的模型直到收发距达到8 km左右时才开受空气波影响，a=∞的模型在所展示的收发距内还未被空气波影响.

图11 海洋三维油气储层模型Fig.11 The marine 3D reservoir model

图12 三维储层模型的储层截面边长变化时的Ex响应(a, b)为不考虑空气层、假设海水无限深的情况，(c, d) 为海水1000 m深、包含空气层的情况.Fig.12 Ex responses for the 3D reservoir model with various horizontal edge lengths of the reservoir(a, b) is for the case air layer is not presented and the seawater is arbitrarily thick, while (c, d) is for the case the seawater is 1000 m thick and air layer is included.

4.4 陆地三维模型

最后对一个陆地三维模型进行正演.如图13所示，在100 Ωm的均匀半空间中植入一个400 m×500 m×400 m的低阻块体，埋深为200 m，电阻率5 Ωm.使用两个分别沿y方向和x方向的正交场源交替发射来模拟张量CSAMT观测，发射导线长度均为1000 m，发射0.125～8192 Hz之间按对数等间隔分布的17个频率；发射源的中心为坐标原点，到低阻体中心的水平距离为6000 m；网格数为62×44×53.

由于有两个正交场源，因此可以使用类似于大地电磁法(MT)中的方式来计算张量阻抗，进而计算卡尼亚视电阻率响应，具体计算式见谭捍东等(2004).图14和图15分别为沿x方向和沿y方向穿过低阻体正上方的两条测线上的响应，图16则为频率为4 Hz的响应切片，均显示了XY和YX两种模式的响应.整体上来看，这些对CSEM的电磁场用MT的方式处理得到的响应与MT响应有一定的相似度，异常体的空间轮廓在正演响应中有良好的反映.XY和YX模式的响应特征有所不同，XY响应对x方向的电性边界反映较好，YX响应则对y方向的电性边界反映较好；在与低阻体对应的位置、对应的频率处，XY视电阻率比YX视电阻率在数值上更接近低阻体的电阻率值.在频率低于2 Hz时，视电阻率值均超过了300 Ωm，甚至可达1000 Ωm，这是因为频率较低时，电磁波在背景地层中的趋肤深度较大，测点所处的位置远远不能满足远区条件，此时卡尼亚视电阻率公式完全不适用.

图13 陆地三维低阻体模型以及发射源的布置示意图Fig.13 The land 3D conductive brick model and the source geometry

5 结论与讨论

本文利用多波前大规模并行求解器软件包MUMPS，实现了基于矩阵三角分解的频率域可控源电磁三维有限体积正演.

针对MUMPS直接解法的计算性能做了一系列测试，结果表明：(1)直接解法非常稳定，不易被系数矩阵的病态程度所影响，在实际CSEM勘探常用的频率范围内，对于包含很大电性差异的模型，都能获得高精度的解，计算速度也很稳定；(2)直接解法的内存需求与计算时间取决于模型网格数，且均随模型网格数的增加而急剧增加，内存需求很容易超出单台普通PC机的物理内存；(3)使用多进程并行执行MUMPS时，其计算速度能得到大幅度提升，也能有效地减少单节点的内存消耗，但由于矩阵分解本身的可并行程度有限，MUMPS的并行可伸缩性也比较有限.

对理论模型正演结果的比较与分析表明：本文的SFV三维正演算法是有效的，适用于常规的陆地和海洋CSEM勘探，并且由于利用了直接解法，非常适合于CSEM中典型的多场源问题，场源的增加只会带来极小的额外计算量.

就目前最普通、最可用的计算设备如个人电脑、工作站和小型服务器的平均性能而言，电磁法的三维数值模拟仍然是一个对计算资源需求很高的问题，使用直接解法时更是如此，因此模型剖分的网格数要严格控制(本文3.2节中对此给出了参考)，这无疑会影响计算结果的准确性.本文的研究便受限于计算设备的性能，在并行计算的进程数、地面CSAMT模拟中的收发距等方面受到了较大的限制.另外，在计算设备的性能足够强劲的情况下，对多频率的正演计算实施“多层次”并行化无疑能进一步提高效率，即将不同频率的正演任务分配到不同的节点组，然后每个节点组对单个频率进行并行正演计算.这些都是有待改进和完善的地方.

随着数值算法的发展以及计算机性能的快速提升，可以预见，直接解法由于其特殊的优点，必将越来越多地被应用于3D电磁模拟中.

致谢 感谢MUMPS软件包和Dipole1D的开发者们，没有他们的开源精神，本文的研究将无法完成；特别感谢美国俄勒冈州立大学(OSU)Adam Schultz教授的指导.

图14 陆地三维模型在x=6000 m的测线上的卡尼亚视电阻率响应拟断面图Fig.14 The Cagniard apparent resistivity response pseudo-sections of the survey line at x=6000 m for the land 3D model

图15 陆地三维模型在y=0 m的测线上的卡尼亚视电阻率响应拟断面图Fig.15 The Cagniard apparent resistivity response pseudo-sections of the survey line at y=0 m for the land 3D model

图16 陆地三维模型在频率f=4 Hz时的卡尼亚视电阻率响应频率切片图Fig.16 The Cagniard apparent resistivity response frequency slices at f=4 Hz for the land 3D model

Abubakar A, Habashy T M, Druskin V L, et al. 2008. 2.5D forward and inverse modeling for interpreting low-frequency electromagnetic measurements.Geophysics, 73(4): F165-F177.

Amestoy P R, Guermouche A, L′Excellent J, et al. 2006. Hybrid scheduling for the parallel solution of linear systems.ParallelComputing, 32(2): 136-156.

Avdeev D B. 2005. Three-dimensional electromagnetic modelling and inversion from theory to application.SurveysinGeophysics, 26(6): 767-799. Avdeev D, Knizhnik S. 2009. 3D integral equation modeling with a linear dependence on dimensions.Geophysics, 74(5): F89-F94. Börner R U. 2010. Numerical modelling in Geo-Electromagnetics: advances and challenges.SurveysinGeophysics, 31(2): 225-245. Chen G B, Wang H N, Yao J J, et al. 2009. Modeling of electromagnetic responses of 3-D electrical anomalous body in a layered anisotropic earth using integral equations.ChineseJ.Geophys. (in Chinese), 52(8): 2174-2181, doi: 10.3969/j.issn.0001-5733.2009.08.028.

Constable S C, Weiss C J. 2006. Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling.Geophysics, 71(2): G43-G51.

Constable S C. 2010. Ten years of marine CSEM for hydrocarbon exploration.Geophysics, 75(5): 75A67-75A81.

da Silva N V, Morgan J V, MacGregor L, et al. 2012. A finite element multifrontal method for 3D CSEM modeling in the frequency domain.Geophysics, 77(2): E101-E115.

Davis P J, Rabinowitz P. 2007. Methods of Numerical Integration. 2nd ed. New York: Dover Publications.

Deng J Z. 2011. Studies of three dimensional stagged-grid finite difference for controlled source audio-frequency magnetotelluric numerical simulation [Ph. D. thesis] (in Chinese). Beijing: China University of Geosciences (Beijing).

Di Q Y, Wang R. 2008. CSAMT Forward Modeling and Inversion and Its Application (in Chinese). Beijing: Science Press. Grayver A V. 2013. Three-dimensional controlled-source electromagnetic inversion using modern computational concepts [Ph. D. thesis]. Berlin: Free University of Berlin.

Han B, Hu X Y, He Z X, et al. 2012. Mathematical classification of magnetotelluric inversion methods.OilGeophysicalProspecting(in Chinese), 47(1): 177-187.

Han B, Hu X Y, Schultz A, et al. 2015. Three-dimensional forward modeling of the marine controlled-source electromagnetic field with complex source geometries.ChineseJ.Geophys. (in Chinese), 58(3): 1059-1071, doi: 10.6038/cjg20150330.

Hu X Y, Peng R H, Wu G J, et al. 2013. Mineral Exploration using CSAMT data: Application to Longmen region metallogenic belt, Guangdong Province, China.Geophysics, 78(3): B111-B119. Jahandari H, Farquharson C G. 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids.Geophysics, 79(6): E287-E302. Key K. 2009. 1D inversion of multicomponent, multifrequency marine CSEM data: methodology and synthetic studies for resolving thin resistive layers.Geophysics, 74(2): F9-F20.

Li J H, Zhu Z Q, Liu S C, et al. 2011. 3D numerical simulation for the transient electromagnetic field excited by the central loop based on the vector finite-element method.JournalofGeophysicsandEngineering, 8(4): 560-567. Luo M, Li Y, Liu Y. 2014. 2.5D controlled-source electromagnetic modeling by an adaptive finite-element method with an arbitrarily oriented finite length dipole source. 22nd EM Induction Workshop. Abstracts. Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences.GeophysicalProspecting, 43(8): 1021-1042.

Oldenburg D W, Haber E, Shekhtman R. 2008. Forward modelling and inversion of multi-source TEM data. ∥ 78th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 559-563. Puzyrev V, Koldan J, de la Puente J, et al. 2013. A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling.Geophys.J.Int., 193(2): 678-693. Saad Y. 2003. Iterative Methods for Sparse Linear Systems. 2nd ed. Philadelphia: SIAM.

Schenk O, Gärtner K. 2004. Solving unsymmetric sparse systems of linear equations with PARDISO.FutureGenerationComputerSystems, 20(3): 475-487.

Schwarzbach C, Börner R U, Spitzer K. 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—a marine CSEM example.Geophys.J.Int., 187(1): 63-74. Shen J S. 2003. Modeling of 3-D electromagnetic responses in frequency domain by using staggered-grid finite difference method.ChineseJ.Geophys. (in Chinese), 46(2): 281-289.

Streich R. 2009. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy.Geophysics, 74(5): F95-F105.

Streich R, Becken M, Ritter O. 2011. 2.5D controlled-source EM modeling with general 3D source geometries.Geophysics, 76(6): F387-F393.

Tan H D, Wei W B, Deng M, et al. 2004. General use formula in MT tensor impedance.OilGeophysicalProspecting(in Chinese), 39(1): 113-116. Wang R, Di Q Y, Wang M Y, et al. 2009. Research on the effect of 3D body between transmitter and receivers on CSAMT response using Integral Equation method.ChineseJ.Geophys. (in Chinese), 52(6): 1573-15821, doi: 10.3969/j.issn.0001-5733.2009.06.019.

Weiss C J. 2001. A matrix-free approach to solving the fully 3D electromagnetic induction problem. ∥71th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1451-1454. Weiss C J, Constable S C. 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part II — Modeling and analysis in 3D.Geophysics, 71(6): G321-G332.

Wu G J, Hu X Y, Liu H. 2010. Progress in CSAMT three-dimensional forward numerical simulation.ProgressinGeophys. (in Chinese), 25(5): 1795-1801, doi: 10.3969/j.issn.1004-2903.2010.05.037.

Xu Z F, Wu X P. 2010. Controlled source electromagnetic 3-D modeling in frequency domain by finite element method.ChineseJ.Geophys. (in Chinese), 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019.

Yang B, Xu Y X, He Z X, et al. 2012. 3D frequency-domain modeling of marine controlled source electromagnetic responses with topography using finite volume method.ChineseJ.Geophys. (in Chinese), 55(4): 1390-1399, doi: 10.6038/j.issn.0001-5733.2012.04.035.

Zaslavsky M, Druskin V, Davydycheva S, et al. 2011. Hybrid finite-difference integral equation solver for 3D frequency domain anisotropic electromagnetic problems.Geophysics, 76(2): F123-F137.

附中文参考文献

陈桂波, 汪宏年, 姚敬金等. 2009. 用积分方程法模拟各向异性地层中三维电性异常体的电磁响应. 地球物理学报, 52(8): 2174-2181, doi: 10.3969/j.issn.0001-5733.2009.08.028.

邓居智. 2011. 可控源音频大地电磁法三维交错采样有限差分数值模拟研究[博士论文]. 北京: 中国地质大学(北京).

底青云, 王若. 2008. 可控源音频大地电磁数据正反演及方法应用. 北京: 科学出版社.

韩波, 胡祥云, 何展翔等. 2012. 大地电磁反演方法的数学分类. 石油地球物理勘探, 47(1): 177-187.

韩波, 胡祥云, Schultz A等. 2015. 复杂场源形态的海洋可控源电磁三维正演. 地球物理学报, 58(3): 1059-1071, doi: 10.6038/cjg20150330.

沈金松. 2003. 用交错网格有限差分法计算三维频率域电磁响应. 地球物理学报, 46(2): 281-289.

谭捍东, 魏文博, 邓明等. 2004. 大地电磁法张量阻抗通用计算公式. 石油地球物理勘探, 39(1): 113-116.

王若, 底青云, 王妙月等. 2009. 用积分方程法研究源与勘探区之间的三维体对CSAMT观测曲线的影响. 地球物理学报, 52(6): 1573-1582, doi: 10.3969/j.issn.0001-5733.2009.06.019.

吴桂桔, 胡祥云, 刘慧. 2010. CSAMT三维正演数值模拟研究进展. 地球物理进展, 25(5): 1795-1801, doi: 10.3969/j.issn.1004-2903.2010.05.037.

徐志锋, 吴小平. 2010. 可控源电磁三维频率域有限元模拟. 地球物理学报, 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019.

杨波, 徐义贤, 何展翔等. 2012. 考虑海底地形的三维频率域可控源电磁响应有限体积法模拟. 地球物理学报, 55(4): 1390-1399, doi: 10.6038/j.issn.0001-5733.2012.04.035.

(本文编辑 胡素芳)

3-D frequency-domain CSEM modeling using a parallel direct solver

HAN Bo， HU Xiang-Yun*， HUANG Yi-Fan， PENG Rong-Hua， LI Jian-Hui， CAI Jian-Chao

HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China

Three-dimensional modeling of electromagnetic data is a computationally demanding problem. For frequently-used numerical techniques such as finite-element and finite-difference methods, solving the large linear systems arising from the discretization of Maxwell's equations is a key step which has a major impact on the applicability of the solution, and it has always been a research topic to solve the linear equations efficiently, robustly and accurately. A 3D modeling scheme based on direct solutions of the linear system is presented for frequency domain controlled-source electromagnetic (CSEM) surveys. The Helmholtz equation in terms of secondary electric fields is discretized using a finite-volume (FV) method over a staggered grid. Taking advantage of recent developments in numerical algorithms and the availability of computational resources, the resulting linear system of FV equations is solved directly using the massively parallel solver, namely MUMPS, instead of the most commonly used linear solvers, i.e. Krylov subspace iterative techniques. The direct solver carries out an LU (and possibly LDLT) decomposition of the system matrix and then computes solutions efficiently by applying forward and backward substitutions.To evaluate the computational performance of the direct solver, a series of numerical tests based on synthetic 1D models were conducted, and the results indicate that (1) Normalized residuals of solutions are almost independent of the conductivity value assigned to air layers but increase rapidly as the frequency value decreases. Nevertheless, the order of magnitude of the largest normalized residual is as small as 10-11. At the same time, although the matrix factorization time varies as either the air conductivity or the frequency changes, the variation is only a fraction of the total run time. (2) Both the execution time and required memory increase rapidly (more than linearly) with increasing grid sizes. (3) By executing MUMPS in parallel over multiple processors, not only the total run time but also the average memory used per processor can be reduced a lot. However, the total memory requirement increases with the number of processes. The scalability of MUMPS is limited. Additional numerical experiments considering specific survey settings were done to demonstrate the reliability and effectiveness of the code, and the results are as follows: (1) The FV numerical solutions show excellent agreement with semi-analytic solutions for the 1D models. (2) The computation time of a multitransmitter problem is comparable to that of a single-transmitter problem. (3) Reasonable modeling results of 3D models can be obtained for both typical land and marine survey scenarios.In summary, compared with iterative linear solvers, the direct solvers generally benefit CSEM modeling in three aspects. The first is they often provide more accurate solutions. The second is that direct solvers are much more stable for ill-conditioned linear systems, which are almost inevitable because of large electrical conductivity contrasts and/or non-uniform grid. The last is in multitransmitter problems, only a single matrix factorization is necessary, and multiple solutions can be achieved very easily by reusing the factors. The presented 3D CSEM modeling scheme, which employs the MUMPS direct solver, possesses all these advantages. In addition, solving linear systems can be executed in parallel to speed up the computation and to reduce the average memory used per node, although the parallel scalability of MUMPS is limited. In spite of the fact that matrix factorizations for large models can entail tremendous computational cost, it can be anticipated that direct solvers will be used more and more widely as the development of both numerical algorithms and computers.

Controlled-source electromagnetics; 3-D modeling; LU factorization; Direct solver

国家自然科学基金项目(41274077、41474055)和中国地质调查局项目(12120113101800)联合资助.

韩波，男，1987年生，博士研究生，研究方向为电磁法的正演与反演模拟.E-mail：hanbo1735@163.com

*通讯作者胡祥云，男，1966年生，教授，主要从事电磁法的理论与应用研究.E-mail：xyhu@cug.edu.cn

10.6038/cjg20150816.

10.6038/cjg20150816

P631

2015-02-12，2015-06-16收修定稿

韩波，胡祥云，黄一凡等. 2015. 基于并行化直接解法的频率域可控源电磁三维正演.地球物理学报，58(8):2812-2826,

Han B, Hu X Y, Huang Y F, et al. 2015. 3-D frequency-domain CSEM modeling using a parallel direct solver.ChineseJ.Geophys. (in Chinese),58(8):2812-2826,doi:10.6038/cjg20150816.