频率域可控源电磁法(controlled-source electromagnetic method，CSEM)由于具有勘探深度大，分辨力高且野外抗干扰能力强等优点，经过30多年快速发展，目前已被广泛应用于包括陆地、海洋及航空等领域内的油气、矿产资源勘探和近地表地质勘查.随着电磁仪器、数据采集技术以及解释方法的快速发展，可控源电磁法在复杂地质环境下的三维勘探成为当前电磁勘探的重要趋势及研究热点.在电磁数据的处理解释中，反演是必不可少的步骤，而正演计算是反演的核心.因此，电磁数据的定量解释首先需要发展有效的正演算法.从数值模拟的角度来说，目前常用于三维电磁正演的数值方法主要包括：积分方程法(Avdeev et al., 2002；Hursán and Zhdanov, 2002；Zhdanov et al., 2006；Avdeev and Knizhnik, 2009)、有限差分法(Mackie et al., 1994；Newman and Alumbaugh, 1995；Smith, 1996；Weiss and Newman, 2002；沈金松, 2003；Streich, 2009)、有限体积法(Haber and Ascher, 2001；Weiss and Constable, 2006；杨波等, 2012；Jahandari and Farquharson, 2014；韩波等, 2015a)和有限单元法(Badea et al., 2001；徐志锋和吴小平, 2010；Schwarzbach et al., 2011；Puzyrev et al., 2013；Ansari and Farquharson, 2014；李勇等, 2015；杨军等, 2015).另外，其他方法如伪谱法(Huang et al., 2010; Liu et al., 2013)、Lanczos谱分解法(Knizhnerman et al., 2009)也被广泛用于三维电磁正演问题的求解.有关电磁法三维正演方面进展的综述可以详见(Avdeev, 2005；Börner, 2010).总体来说，目前三维正演算法研究的重点在于提高精度和计算效率.

CSEM正演中常采用散射场方法来避免场源的奇异性，然而构建散射场的场源项需要计算几乎每个离散网格采样点处的一次电场或磁场.对于三维问题，离散网格中的电场或磁场采样点个数大约是网格数的3倍，很容易达到几百万甚至上千万.因此构建散射场场源往往需要大量的时间，随着发射场源的增加，甚至有可能超过解线性方程需要的时间(韩波等, 2015b).相较之下，总场方法直接对物理场源离散化，不存在这个问题.但总场方法中一般需要在物理场源附近加密网格来削弱场源奇异性的影响，大量网格数同样会导致大量的计算时间.除了加密网格以外，还可采用场源校正技术等手段来降低场源奇异性的影响.

正演计算中的一个关键步骤是求解离散化后产生的大型线性方程组.在早期受限于计算机内存，三维正演中大型线性方程组几乎都采用对内存需求极低的Krylov子空间迭代法.随着大型稀疏矩阵分解算法的不断优化(Amestoy et al., 2001, 2006；Li, 2005；Schenk and Görtner, 2006)以及计算机硬件技术的快速发展，利用直接解法求解三维正演的大型线性方程组成为可能.与迭代解法相比，对于同一线性方程组，直接解法通常需要消耗更多的计算内存和更长的计算时间.尽管如此，直接解法相对于迭代解法有两大显著优点：一是求解时间和精度基本不受系数矩阵的条件数影响，因此网格的剖分方式或介质的电性差异对正演算法影响很小，使得求解精度更高，正演算法更稳定；另一个是直接解法的主要开销集中在系数矩阵的分解阶段，对于单一频率多个场源位置的问题，只需进行一次矩阵分解加上多次计算量极小的前代-回代过程即可，而不像迭代解法对每个频率每个场源都要单独求解，因此直接解法特别适合具有多个场源的CSEM三维模拟(Oldenburg et al., 2013).正因为这些优点，基于矩阵分解的直接解法最近几年开始应用于CSEM三维正演计算(Streich, 2009；Chung et al., 2014；杨军等, 2015).

本文采用拟态有限体积法(Mimetic Finite Volume, MFV)，通过求解离散化的三维电场矢量Helmholtz方程，实现了频率域CSEM三维正演算法.为获得稳定且高精度的数值结果，采用直接矩阵分解法来求解离散所得到的大型稀疏线性方程组；为消弱场源奇异性及边界条件的影响，采用场源校正技术来提高求解精度.对于多频率正演问题，采用分频并行策略来降低正演计算时间；最后通过多个理论合成模型正演计算测试，验证了本文所开发的正演算法对频率域CSEM模拟计算的准确性及有效性.

在CSEM所采用的频率范围内，位移电流的影响可以忽略不计.假定时谐因子取为e^iwt，则频率域Maxwell方程组可表示为

(1)

(2)

其中，E为电场强度(V·m^-1)，B为磁感应强度(T)，ω为角频率(rad/s).σ为介质电导率(S·m^-1)，对于各向同性介质，σ为标量.μ₀为真空中磁导率(H·m^-1). J_s为外部激发场源的电流密度(A·m^-2).对(1)式取旋度并将其带入(2)式，消去磁感应强度B，得到关于电场E的二阶矢量Helmholtz方程为

(3)

为求解方程(3)，本文采用拟态有限体积法(MFV)来对连续的Maxwell方程组进行离散化.MFV方法的最大特点是能够对连续的微分算符进行精确的离散模拟，并确保离散化后的矢量场仍然满足其对应的连续形式的矢量性质及物理特性(Hyman and Shashkov, 1999a, b).如离散化的电磁场的能量守恒能够自动满足，有效避免了伪解的产生.对于正交规则网格及各向同性介质来说，该离散化方法与常用的交错网格有限差分法(Staggered Finite-Difference, SFD)相一致(Smith, 1996).尽管最近也有将非正交交错网格有限差分法应用于电磁数值模拟的研究(邱稚鹏等，2013)，但MFV适用范围更广，其不仅对于非正交网格同样适用(Hyman and Shashkov, 1997)，而且对于具有高度电性差异及各向异性介质，MFV能够自然地得到对称的离散化形式(Hyman et al., 1997；Haber and Ruthotto, 2014)，这对于复杂三维介质的模拟十分有效.

为获得Maxwell方程组的离散形式，首先考虑方程(1)和(2)的弱解形式为(Hyman and Shashkov, 1999a)：

(4)

(5)

其中W和F分别为与电场E和磁感应强度B属于相同Hilbert空间的任意测试函数(Hyman and Shashkov, 1999a).对(5)式左端项进行分部积分，可得公式为

(6)

其中n为计算区域边界的外向单位法向向量.在可控源电磁法中，由于人工场源的局部性，一般采用理想导体边界条件.即假设场源所激发的电磁场在距离目标靶区足够远的区域内可以忽略不计，此时计算区域边界上的切向电场分量均取为零.本文采用的理想导体边界条件为(Perfect Electric Conductor)

(7)

将(6-7)式代入到(5)式中，得到公式为

(8)

由于MFV是在离散单元的控制体积内进行积分，因此正交结构网格(Weiss and Constable, 2006；韩波等, 2015a)、半结构网格(Haber and Heldmann, 2007)和非结构网格(Jahandari and Farquharson, 2014)都能用来对(4)式和(8)式的积分进行离散化.本文采用基于Yee网格(Yee, 1966)的交错网格.图 1展示了交错网格电磁场采样方式及电场分量和磁场分量各自的控制体积.其中电场定义在单元棱边的中心，磁感应强度定义在单元面的中心，介质电导率和磁导率定义在单元中心.

图 1

Fig. 1

图 1 电磁场分量交错采样方式及其控制体积 (a) 电场x分量采样位置及其控制体积； (b) 磁场x分量采样位置及其控制体积. Fig. 1 Staggered discretization of MFV (a) Integration volume for x-component of electric field； (b) Integration volume for x-component of magnetic field.

通过在控制体积内进行近似积分，基于交错网格离散化的(4)式和(8)式可以简化成矩阵形式为

(9)

(10)

其中向量为所有网格中心的电导率值，向量为所有网格单元的体积，为所有网格棱边中心采样点的电场分量，为对应的棱边网格函数.为所有网格单元面中心采样点的磁场分量，为对应的单元面中心网格函数.Curl为旋度算子的离散式.Θ为Hadamard矩阵乘积算子，表示两个矩阵相同位置的元素相乘.diag[·]为对角矩阵算符.A_ec和A_fc为平均矩阵，分别将变量从网格单元棱边中心及网格单元面中心映射到网格单元中心.

由于F和W为任意的网格测试函数，将(9)式代入(10)式并进行简化可得(3)式的离散形式为

(11)

其中

(11)式可简化成大型线性方程组为

(12)

其中系数矩阵A为稀疏、正定、对称复矩阵.图 2给出了网格数为5×5×5的模型的MFV离散所得到的系数矩阵A的稀疏结构.

图 2

Fig. 2

图 2 网格数为5×5×5的模型的MFV离散线性方程系数矩阵的稀疏结构 Fig. 2 Sparse structure of coefficient matrix A resulted from MFV discretization on a 5×5×5 staggered grid

在可控源电磁法中，由于外加场源的存在，待求解的电磁场在场源附近会发生急剧变化，场源的奇异性会极大降低数值计算结果的精度.目前主要有两种处理方法来降低场源奇异性的影响.第一种方法是采用散射场公式(Newman and Alumbaugh, 1995；Streich, 2009；韩波等, 2015a)，将待求解的电磁场分解为由参考模型(一般为均匀空间或一维层状模型)在外加场源激发下所产生的一次场和三维模型所产生的二次场的叠加.由于外加场源已包含在一次场的计算中，二次场的场源不再具有奇异性.另一种是基于电磁场总场公式，通过对场源点附近进行网格加密的方法来减少场源奇异性的影响(Plessix et al., 2007).

在散射场公式中，为避免计算二次场时出现场源的奇异性，必须确保选择的参考模型的电阻率与场源点附近的电阻率一致.对于海洋CSEM测量来说，由于发射场源一般位于海水中，一般容易满足该条件.但对于陆地CSEM测量来说，通常CSEM勘探会在较大测区内布设多个发射位置(Streich et al., 2011)，当场源处的电阻率发生变化时，需要求解不同参考模型的一次电磁场响应，从而增加计算量.另外，在散射场公式中，每个场源处都需要计算三维空间内所有电磁场采样点处的一次电磁场来构建场源项或进行总场合成.随着发射场源位置的增加，场源项的构建所消耗的时间会快速增加，甚至会超过系数矩阵分解所消耗的时间(韩波等, 2015b).与之对比，在总场公式中，网格加密的策略虽然能够降低场源奇异性的影响，但与此同时会显著地增加待求解的未知量个数，从而增大计算开销，降低求解效率.另外，为满足边界条件，无论是采用散射场公式还是总场公式，通常都需要将计算区域的外边界设置得足够远，这使得待求解的方程组会非常巨大.

为降低场源奇异性及边界条件对数值精度的影响，本文采用虚拟场源校正技术(Pidlisecky et al., 2007)来对总场控制方程(11)式中的源项进行校正，避免了在场源点附近进行网格局部加密所造成的计算量的增大.Pidlisecky等(2007)、Pidlisecky和Knight(2008)将该校正技术运用到直流电阻率的正演计算中，取得了很好的效果.为得到校正后的场源项，考虑存在解析解的均匀半空间模型σ₀，利用上节发展的有限体积法对其离散化可得公式为

(13)

其中系数矩阵A(σ₀)和右端项S₀均由(11)式确定.假设E₀是模型σ₀在外加场源J_s激励下的网格采样点处的电场响应，可通过解析求解很容易得到(Ward and Hohmann, 1988)，由于离散化误差及边界条件的影响，待求解的电场分量E与真实的电场响应E₀会有一定差异.此时，我们可以得到场源校正量为

(14)

将(14)式所得到的场源校正量添加到(12)式右端，得公式为

(15)

通过求解(15)式便可以得到待求解的电场值.

由于离散后所得到的系数矩阵为稀疏、正定、对称复矩阵，本文通过调用基于多波前分解算法的MUMPS (Amestoy et al., 2001, 2006)线性运算库对其进行LDL^H分解，然后计算所有发射场源的CSEM电磁响应.有关MUMPS的运算效率及与迭代法的对比可参考Streich(2009)和Oldenburg等(2013)，另外关于MUMPS的求解精度、稳定性及内存需求等性能的详细分析，可参考韩波等(2015b).

为增大勘探深度及提高数据空间分辨率，CSEM测量中通常会采集多个频率的数据.对于三维正演来说，为降低计算时间，通常需考虑对其进行并行计算.一种思路是利用并行矩阵分解算法对多个频率依次进行求解.由于矩阵并行分解算法涉及到不同进程之间大量的通信需求，并行效率不甚理想；并且，随着并行进程数的增加，总的内存消耗也会迅速增加(Pardo et al., 2012；韩波等, 2015b).对于频率域CSEM来说，考虑到不同频率之间的正演计算是相互独立的特征，在计算资源允许的条件下，可以对多个频率同时进行并行矩阵分解计算.Grayver等(2013)的测试结果表明，分频并行策略比采用矩阵并行分解算法依次对多个频率进行分解计算的效率要高很多.因此本文采用分频并行策略，将所有频点尽可能均匀地分配到所有进程中进行并行计算.

并行程序采用主从并行模式.其中主进程作为控制节点，负责模型网格、收发参数的发送以及计算任务的分配，并对各子进程的计算结果进行回收以及结果的输出，不参与计算；各子进程负责从主进程处接受计算任务并将计算结果发回给主进程.在整个并行计算过程中，各子进程只与主进程进行通信.由于本文所使用的并行计算平台计算资源(可调用的计算节点数和计算内存)的限制，在子进程中，对于每个频点的计算，并未采用矩阵并行分解算法.当分配有多个频点时，需依次对多个频率进行分解计算.整个并行算法的计算流程如图 3所示.

图 3

Fig. 3

图 3 CSEM 正演并行计算流程图 Fig. 3 Flowchart of 3D parallel CSEM forward modeling

本节通过对多个理论模型的正演数值结果的对比分析来检验本文所开发的正演算法对频率域CSEM模拟计算的准确性及有效性.

为了验证本文所开发的CSEM三维正演算法的准确性，首先考虑具有准解析解的一维层状模型.图 4a是海洋CSEM中用于模拟高阻油气层的标准模型(Constable and Weiss, 2006).为模拟强电性差异，油气层电阻率设为1000 Ωm，厚度为100 m，顶部距海底1000 m；海底沉积层电阻率为1 Ωm，海水层厚度为1000 m，电阻率为0.33 Ωm.图 4b为典型陆地一维层状模型.电阻率为100 Ωm沉积层中包含一层电阻率为10 Ωm的低阻目标层，该目标层顶部埋深200 m，厚度为400 m.所有模型的空气层电阻率均取10⁹ Ωm.本文所有计算均在一个小型并行机上完成.该小型并行集群系统配置有8个计算节点，每个计算节点含有2个4核Intel^®Xeon^®E5410型处理器，主频为2.33 GHz，每个计算节点配置有32G内存.操作系统为CentOS 5.5.层状模型的拟解析解均通过Dipole1D程序(Key，2009)求得.

图 4

Fig. 4

图 4 海洋(a)及陆地(b)一维层状模型 Fig. 4 Sketches of (a) marine-based and (b) land-based 1D layered models

对于海底层状模型(图 4a)，其观测系统设置为：电偶极子位于海底上方100 m，发射频率取0.25 Hz和1 Hz，电流强度为1A.接收器位于海底表面，布设方向与偶极方向一致(inline观测方式).模型剖分网格数为132×52×52.图 5给出了不同频率下E_x分量和B_y分量的振幅-收发距(Magnitude Versus Offset, MVO)曲线与相位-收发距(Phase Versus Offset, PVO)曲线的MFV三维数值解与解析解的对比.从图中可以看出，不同频率的电磁场分量的MFV数值解与解析解均高度一致.图 6则给出了MFV解相对于解析解的误差情况.当未进行场源校正时，可以看到在发射偶极附近(收发距小于1000 m)，电磁场的各个分量均具有较大的误差，振幅最大误差可达10%左右，相位最大误差可达2°，这主要是由于场源的奇异性造成的.需要说明的是，为降低待求解的线性方程组的尺寸，减少直接求解的计算量，本文并未在场源点附近进行精细的网格加密来降低场源奇异性的影响.当采用场源校正技术后，在场源附近，场源奇异性对于数值精度的影响会迅速降低，振幅最大误差降为5%以内，相位最大误差降为1.5°以内.另外，随着收发距的增加，场源奇异性的影响迅速降低，误差也随之迅速减小，电磁场分量的振幅误差降为2%以内，相位误差在2°以内.对于海洋CSEM来说，近场区内(收发距小于1000 m)一次场占据主导地位，一般并不包含海底目标层信息(Zhdanov et al., 2014).另外，从数值计算结果对比来看，场源奇异性对于远场区的计算精度影响有限，因此采用总场公式对于海洋CSEM的正演计算是完全可行的.

图 5

Fig. 5

图 5 海洋层状模型CSEM数值解与拟解析解对比 (a)(c)为电磁场分量的振幅，(b)(d)为对应电磁场的相位. Fig. 5 Accuracy comparison between MFV numerical solutions and quasi-analytic solutions for marine 1D layered model (a)(c)are amplitudes of field components,(b)(d)are corresponding phases.

图 6

Fig. 6

图 6 海洋层状模型CSEM数值解相对于拟解析解的误差 (a)(c)为电磁场分量振幅的相对误差，(b)(d)为对应电磁场的相位差异. Fig. 6 Errors between MFV numerical solutions and quasi-analytic solutions for marine 1D layered model (a)(c)indicate relative errors in amplitude,(b)(d)are corresponding phase difference.

对于陆地层状模型(图 4b)，观测系统采用可控源音频大地电磁法(CSAMT)(底青云和王若，2008)的观测方式.其中发射源为1000 m的接地导线，发射频点数为14个，频率在0.25~4096 Hz之间呈对数等间隔分布；测线平行于发射源布设.模型剖分网格数为99×56×40.图 7给出了与场源中心距离为6 km的测点上所有频率的E_y分量和B_x分量的MFV数值解与解析解对比.从图中可以看出，所有频点的数值解与解析解均非常吻合.图 8则给出了数值解相对于解析解的误差情况.从图中可以看出，当未进行场源校正时，大多数频点的电磁场分量的振幅误差为6%~8%；与之对比，采用场源校正技术后，电磁场分量的振幅误差降为5%以内.另外，在进行场源校正后，电磁场分量的相位误差略有增加，但仍保持在2°以内.在野外CSAMT测量中，由于各种人文噪声的影响，视电阻率误差和相位误差一般会比较大(Hu et al., 2013)，经过场源校正后，采用总场公式来进行CSAMT的正演模拟完全能够满足实测精度要求.

图 7

Fig. 7

图 7 层状模型CSEM数值解相对于拟解析解的误差 (a)(c)为电磁场分量的振幅，(b)(d)为对应电磁场的相位. Fig. 7 Accuracy comparison between MFV numerical solutions and quasi-analytic solutions for land 1D layered model (a)(c)are amplitudes of field components,(b)(d)are corresponding phases.

图 8

Fig. 8

图 8 层状模型CSEM数值解相对拟解析解的误差 (a)(c)为电磁场分量振幅的相对误差，(b)(d)为对应电磁场的相位差异. Fig. 8 Errors between MFV numerical solutions and quasi-analytic solutions for marine 1D layered model (a)(c)indicate relative errors in amplitude,(b)(d)are corresponding phase difference.

在1D层状模型测试的基础上，本节对海洋3D模型进行了正演测试，并与已经发表的采用散射场公式的正演算法(韩波等, 2015b)所得到的结果进行了对比验证.

图 9为典型海洋三维油气藏模型.油气藏模型尺寸为4 km×4 km×0.1 km，电阻率为100 Ωm，其顶部埋深距海底1000 m.均匀海底沉积层电阻率为1 Ωm，海水层厚度为1 km，电阻率为0.33 Ωm.观测系统设置为：发射偶极子沿x方向，位于海底上方100 m处，坐标为(0, 0, -100 m).接收站沿y=0的测线布设于海底，测线方向与偶极方向一致.发射频率为1 Hz.

图 9

Fig. 9

图 9 海洋三维油气藏模型平面(圆点为接收站位置，箭头表示发射源位置) (a) 水平切片图z=1050 m； (b) 垂直剖面图y=0 m. Fig. 9 Plan views of 3D marine hydrocarbon reservoir (a) Horizontal slice at z=1050 m； (b) Vertical section at y=0 m.

图 10给出了对于海洋三维模型，采用本文算法所得到的数值结果与散射场算法的正演结果对比图.从图中可以看出，无论是电场分量(图 10a, b)还是磁场分量(图 10c, d)，利用本文算法所得到的正演结果均与散射场算法的结果吻合的非常好，尽管在大偏移距处(>5 km)，两者的结果有略微的偏差.

图 10

Fig. 10

图 10 海洋三维模型数值结果对比 (a)(c)电磁场分量振幅； (b)(d)电磁场分量相位.其中方框为本文正演所得到的结果，圆圈为韩波等(2015a)计算的结果. Fig. 10 Comparison of data obtained from our modeling scheme to data resulted from Han et al. (2015a) for 3D marine model (a) Ex amplitudes； (b) Ex phases； (c) By amplitudes； (d) By phases.

为了检验本文所开发CSEM三维正演算法的并行效率，对陆地1D层状模型的多频率正演响应进行了并行计算测试，按对数等间隔在0.25~4096 Hz之间取12个频率.表 1给出了使用不同数量的计算节点时的运行时间和并行加速比.可以看出，随着计算节点数的增加，总的计算时间迅速降低.由于各个频率之间的计算是完全独立的，主进程与各个子进程只涉及到非常少的通信(计算参数的分配及结果回收)，从并行加速比来看，并行效率非常接近于理想情况(加速比等于计算进程数).

表 1(Table 1)

表 1 陆地1D模型，CSEM正演并行计算时间统计表 Table 1 Statistics of runtime of CSEM forward modeling for land-based 1D model 程序运行 模式 计算 节点数 每个节点 分配的 频点个数 模型剖分 网格大小 程序运行 时间(s) 并行 加速比 串行模式 1 12 99×56×40 5315.4 无 并行模式 1 12 99×56×40 5357.2 0.99 2 6 99×56×40 2688.3 1.98 4 3 99×56×40 1372.6 3.87 6 2 99×56×40 935.9 5.68

表 1 陆地1D模型，CSEM正演并行计算时间统计表 Table 1 Statistics of runtime of CSEM forward modeling for land-based 1D model

(1) 本文采用直接解法来求解离散化的电场矢量Helmholtz方程，实现了基于有限体积法的频率域CSEM三维正演计算.对多个理论模型的正演测试结果表明本文的算法对于典型海洋和陆地CSEM测量是有效的.

(2) 由于采用直接解法来求解线性方程组，一次矩阵分解的结果可以用于不同发射源位置的计算，特别适合具有多发射源的CSEM测量；此外，在矩阵分解中，系数矩阵的病态程度对于矩阵分解结果影响很小，使得基于矩阵分解的直接解法对于包含强电性差异的地电模型都可以获得稳定且高精度的解.

(3) 采用了总场公式，场源的模拟采用直接离散化的方式.与基于散射场公式的正演算法相比，场源项的构建并不需要多次求解均匀模型或1D模型的一次电磁场，发射源位置的增加只引起计算量的少量增加，同样有利于多发射源的CSEM测量.与此同时，为避免对场源点附近进行网格局部加密所引起的计算量的增加，采用虚拟场源校正技术来降低场源奇异性及边界条件对数值精度的影响，取得了较好的效果.

(4) 三维电磁模拟对于计算资源要求很高，并且往往花费较长的计算时间.因此提高计算效率对于快速的电磁三维正反演至关重要.随着个人工作站和小型服务器的普及，一个有效的策略是考虑并行计算.对于多频率CSEM测量来说(如CSAMT)，基于分频并行的策略能够有效降低总的计算时间.另外，在计算资源允许的条件下，还可考虑对矩阵分解算法进行并行，从而进一步提高并行效率.

感谢加拿大英属哥伦比亚大学(UBC)Eldad Haber教授的指导及与GIF组内同学进行的有益讨论，感谢MUMPS线性运算库的开发者们，感谢两名匿名审稿人对本文提出的宝贵修改意见.