频率域可控源电磁法(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的求解精度、稳定性、内存需求、计算时间以及并行可伸缩性等计算性能，并验证所开发的算法和程序的有效性和广泛适用性.

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

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

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

图 1

Fig. 1

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

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

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

图 2

Fig. 2

图 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

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

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

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

图 3

Fig. 3

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

为了测试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.

图 4

Fig. 4

图 4 用来测试直接解法性能的(a)海洋与(b)陆地一维层状模型 Fig. 4 The(a)marine layered 1D model and (b)l and 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 Table 1

表 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 频率/Hz 空气层电导率(S/m) 残差 相对残差 矩阵分解耗时/s 0.01 10-5 9.10×10-21 6.12×10-11 1122 10-8 9.09×10-21 6.12×10-11 1138 10-11 8.99×10-21 6.06×10-11 1133 0.1 10-5 9.36×10-21 5.57×10-12 1123 10-8 9.58×10-21 5.70×10-12 1146 10-11 9.37×10-21 5.57×10-12 1123 1 10-5 7.20×10-21 4.58×10-13 1119 10-8 7.30×10-21 4.63×10-13 1138 10-11 7.32×10-21 4.65×10-13 1125 10 10-5 2.25×10-22 3.52×10-14 1111 10-8 2.21×10-22 3.46×10-14 1118 10-11 2.23×10-22 3.49×10-14 1134

表 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

表 2 Table 2

表 2 对陆地1D模型，频率与空气层电导率在较大范围内变化时直接解法的精度和矩阵分解时间 Table 2 The accuracy and factorization time of the direct solver for the l and 1D model for a wide range of frequencies & air conductivities 频率 /Hz 空气层电导率(S/m) 残差 相对残差 矩阵分解耗时/s 0.0625 10-5 1.76×10-17 1.81×10-11 1874 10-8 1.66×10-17 1.70×10-11 2074 10-11 1.69×10-17 1.74×10-11 2115 1 10-5 2.13×10-17 1.37×10-12 1689 10-8 2.07×10-17 1.33×10-12 1949 10-11 2.07×10-17 1.33×10-12 2024 16 10-5 2.04×10-17 8.27×10-14 1474 10-8 1.89×10-17 7.65×10-14 1857 10-11 1.90×10-17 7.69×10-14 2011 256 10-5 2.48×10-17 6.55×10-15 1428 10-8 2.33×10-17 6.16×10-15 1744 10-11 2.22×10-17 5.85×10-15 1945 4096 10-5 3.88×10-17 2.29×10-15 1397 10-8 4.19×10-17 2.47×10-15 1547 10-11 5.09×10-17 3.00×10-15 1876

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

再来考察直接解法的内存需求与计算时间对模 型剖分网格数的敏感程度.对海洋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)，内存需求在16 GB左右，而当网格数达到98×80×57时，内存需求则达到了接近120GB.虽然无法给出网格数与内存需求之间的具体关系式，但不难看出整体上内存需求增加的速度要远快于网格数增加的速度.

图 5

Fig. 5

图 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

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

图 6

Fig. 6

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

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

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

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

图 7

Fig. 7

图 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

Fig. 8

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

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

图 9

Fig. 9

图 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 l and 1D model.(g，h)are the apparent resistivity errors of SFV solutions compared to analytic solutions

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

图 10

Fig. 10

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

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

图 11

Fig. 11

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

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

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

图 12

Fig. 12

图 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.

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

最后对一个陆地三维模型进行正演.如图 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.

图 13

Fig. 13

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

由于有两个正交场源，因此可以使用类似于大地电磁法(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，这是因为频率较低时，电磁波在背景地层中的趋肤深度较大，测点所处的位置远远不能满足远区条件，此时卡尼亚视电阻率公式完全不适用.

图 14

Fig. 14

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

图 15

Fig. 15

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

图 16

Fig. 16

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

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

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

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

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

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

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