1 引言

频率域电磁法(Frequency-Domain Electromagnetic，FDEM)通过观测和分析不同频率的人工场源或天然场源激发地下介质产生的二次场或总电磁场分布规律以探明地下电性分布特征.该方法种类繁多，按照场源性质(Nabighian，1988)可分为平面波场的大地电磁法(MT)、音频大地电磁法(AMT)、甚低频电磁法(VLF)；接地导线激发的可控源电磁法(包括地面可控源音频大地电磁CSAMT和海洋可控源电磁法MCSEM)；不接地回线激法的电磁剖面法(包括地面电磁和航空电磁法AEM)等.这些方法在深部地球结构探测(Selway，2014)、地热勘探(Muñoz，2014)、矿产勘查(Smith，2014)、油气勘探(Strack，2014)以及坏境与工程勘探(Sheard et al.，2005)等各领域起到关键作用.

频率域电磁法三维正、反演已成为目前电磁勘探研究的热点和难点问题.其中，MT三维正反演技术已基本成熟，并开始走向实用化(Siripunvaraporn，2012；Miensopust et al.，2013).Avdeev(2005)，Börner(2010)和Newman(2014)对电磁法三维正反演技术发展和挑战进行综述，指出在正演算法中积分方程法(Wannamaker，1991; Song et al.，2002; Farquharson et al.，2006; Zhdanov et al.，2006; Avdeev and Knizhnik，2009; Bakr and Mannseth，2009)通常利用Green函数求解关于二次电场的Fredholm积分式，方法已基本成熟，但对复杂模型的正演精度不高；有限单元法(Yoshimura and Oshiman，2002; Mukherjee and Everett，2011; Pardo et al.，2011; Yahya et al.，2012; Schwarzbach and Haber，2013；Cai et al.，2014)通常利用变分原理或加权余量法求解关于电场、磁场或耦合势的微分方程，能够模拟复杂地形和地下结构，但计算量巨大、求解速度慢；有限差分或有限体积法既能模拟复杂模型，又具有较快计算速度，已成为三维电磁正反演中主流算法.Mackie等(1993)使用有限差分法实现了三维大地电磁正演模拟并将计算结果与积分方程法进行了对比，验证了方法的正确性；Smith(1996a，1996b)对交错采样有限差分法的原理、误差分析以及加快计算速度的双共轭梯度法进行了详细阐述；谭捍东等(2003)，Tan等(2006)系统论述了消去电场分量得到关于磁场的大地电磁三维交错采样有限差分数值模拟算法，并对实现过程中交错网格剖分、积分公式离散化、边界条件、方程组求解、三维张量阻抗的计算等内容做出研究，并实现了并行计算；陈辉等(2011)在此基础上提出利用磁场散度实校正方法(RRCM)提高三维正演精度和加快正演速度；沈金松(2003)利用交错网格有限差分法消去磁场分量得到关于电场的线性方程组，实现了三维频率域电磁响应正演模拟，并且在文中提出在迭代过程中施加了散度校正的方法；Weiss和Constable(2006)利用交错网格有限差分求解二次电场的双旋度方程实现海洋电磁法三维数值模拟；邓居智等(2011)利用交错网格有限差分求解二次磁场的双旋度方程实现可控源音频大地电磁模拟；殷长春等(2014)利用交错网格有限差分求解二次电场的双旋度方程实现海洋电磁法各向异性三维正演模拟；Haber等(2000)针对基于电或磁场的双旋度方程在空气电阻率无穷大引起非椭圆方程和地面强烈耦合作用造成求解不准问题，提出利用有限体积法求解Coulomb规范条件下矢势和标势耦合方程，以实现频率域电磁法三维正演，其离散线性方程组是正定的，但非对称.Hou等(2006)和张烨等(2012)利用有限体积法求解Coulomb规范条件下二次矢势和标势耦合方程，实现感应测井三维各向异性模拟；周建美等(2014)利用该方法实现海洋电磁法三维各向异性正演模拟.综上所述，为了解决人工场源源项处理和场源畸变问题，通常采用数值算法求解二次电场和磁场双旋度方程；考虑到空气层和地下介质耦合问题，人们近年提出求解Coulomb规范条件下二次场标势和矢势耦合方程，但该算法具有网格节点数巨大，求解一次场非常耗时及Coulomb规范条件产生的耦合方程非对称等缺点.

本文从麦克斯韦方程组出发，在Lorenz规范条件下建立磁矢势和标势耦合天然对称方程，采用交错采样网格和有限体积法离散磁矢势和标势方程，并将不同场源类型分解成一系列短导线(电性)源组合，对离散线性方程组采用Jacobi迭代的预处理和拟最小残差法(QMR)求解，从而实现任意场源激发下频率域三维电磁正演.最后，我们通过与各种不同类型一维理论解进行对比验证了该算法的有效性，并进一步对比分析不同场源激发下典型地电模型的电磁响应特征.

2 磁矢势和标势耦合方程及边界条件

在电磁勘探中，对大多数地下岩石我们可以忽略位移电流.假定时间变化因子为e^iωt，则准静态条件下麦克斯韦方程式可写为

(1a)

(1b)

(1c)

(1d)

式中，E为电场，H为磁场，J_e为电流密度，μ为磁导率，ε为介电常数，σ为电导率.

式(1d)表明，磁场是无源场.因此我们引入一个磁矢势B=∇×A.根据式(1a)，电场可以用磁矢势A和一个标势φ表示为

(2)

假设介质磁导率μ和介电常数ε均为常数，将式(2)代入式(1b)可得

(3)

为了保证磁矢势A和标势φ唯一及方程的对称性，我们引入洛伦兹(Lorenz)规范条件：

(4)

其中c=iωμ.将式(4)代入式(3)可得

(5)

对式(3)两边取散度，并代入Lorenz规范条件(4)可得

(6)

联立式(5)和式(6)，我们建立了基于Lorenz规范条件的磁矢势和标势耦合方程组：

(7)

可以看出，该耦合方程具有天然对称性，对离散后形成的线性方程组求解精度和稳定性有一定提升.由于地下介质为有损介质，电磁场随传播距离增加呈指数衰减，在正演模拟中可选择足够大的求解区域Ω，在求解区域Ω外的电磁场值非常小.因此可以采用简单的截断边界条件.与其对应的Lorenz规范条件下磁矢势和标势的边界条件可表示为

(8)

3 频率域三维有限体积电磁正演模拟

为了求解磁矢势和标势耦合方程(7)，我们将系统阐述有限体积法离散过程、场源处理及线性方程组合成求解等.

3.1 区域离散

首先，将求解区域Ω采用一系列平行于x、y、z坐标轴的三组平面剖分成若干个小的六面体网格单元.假设沿x轴方向被剖分成N_x段，含有l个节点(l=N_x+1)；节点编号沿x轴方向递增i=1，2，…，l，网格间距为dx_i(i=1，…，N_x)；沿y轴方向被剖分成N_y段，含有m个节点(m=N_y+1)，节点编号沿y轴方向递增，j=1，2，…，m，网格间距为dy_j(j=1，…，N_y)；沿z轴方向被剖分成N_z段，含有n个节点(n=N_z+1)，节点编号沿z轴方向递增，k=1，2，…，n，网格间距为dz_k(k=1，…，N_z);共计有N_x×N_y×N_z个长方体剖分单元.

根据电场与标势和矢势的关系式(2)，我们发现电场为磁矢势A与标势梯度∇φ之和.因此，对任意一个长方体剖分单元(如图 1)，我们定义磁矢势分量在长方体的边上，记录点为边的中点，而标势定义在长方体的节点上，电导率定义在整个长方体单元内.这种交错采样能够保证磁矢势A和标势梯度∇φ在同一个记录点上自然满足电场与势的转换关系，同时保证相邻单元间电场切向分量连续和磁场法向分量的连续性.

3.2 方程离散

为实现有限体积法离散，需要定义磁矢势三个分量和标势体积单元大小.我们以A_x分量为例.参见图 2a，体积单元沿x方向边长为dx_i，沿y和z方向的边长分别为(dy_j-1+dy_j)/2和(dz_k-1+dz_k)/2.因此，该单元体积为

(9)

按照等效电阻率原理，该体积单元电导率可以表示为

(10)

同理可得A_y，A_z分量体积单元表达式.对于标势φ体积单元(参见图 2b)，沿x，y和z方向的边长分别为(dx_i-1+dx_i)/2，(dy_j-1+dy_j)/2，(dz_k-1+dz_k)/2.由此，单元体积为

(11)

而其体积单元等效电导率为

(12)

对磁矢势和标势耦合方程组式(7)第一式A_x分量，在某个A_x体积单元上积分可得

(13)

对于(13)式左端四项采用有限体积离散，并结合(10)和(12)式的单元等效电阻率可得到：

(14)

(15)

(16)

(17)

对于磁矢势和标势耦合方程组(7)第一式磁矢势其他两个分量A_y，A_z推导过程基本类似，只是右端项场源不同.本文不再详述.

对于磁矢势和标势耦合方程组(2)第二式，在某个标势φ体积单元上积分可得

(18)

仿照前面离散过程，(18)式中左端四项的离散公式为

(19)

(20)

(21)

(22)

3.3 场源设置

对于方程组(7)右端源项J_e和∇·J_e的处理，根据前面有限体积思想可以分别看成在某个剖分单元内电流密度的体积分和面积分，可进一步等效成各剖分单元内电偶极子和剩余电流大小.频率域电磁法场源分为有限长导线的电性源、不接地回线磁性源以及平面电磁波天然场源.为了建立统一的求解系统，本文将不同场源都看成由若干短导线电流源组成，并且在求解区域离散中把其设置在每个剖分单元的棱边上.对于有限长导线的电性源(图 3a)，我们可以将其看成同一方向许多短导线电流源合成；而不接地回线磁性源(图 3b)可以看成许多短导线电流源合成四个方向有限长导线组成的回路.对于平面电磁波天然场源，可等效为沿着某一方向电流源极化的平面电磁波(图 3c)，因而可将其分解成整个平面内无数个同向短导线电流源组成.

3.4 离散线性方程组合成及求解

将(14)—(17)式代入方程组(7)可以建立磁矢势A_x分量和标势φ的离散表达式(见图 4a)，同理将(19)—(22)式代入方程组(7)可建立标势φ和磁矢势A_x、A_y、A_z分量的离散表达式(见图 4b).将方程组(7)离散后A_x、A_y、A_z、φ四个等式，并且结合边界条件(8)式可以合成总体线性方程组：

(23)

式中C_x，C_y，C_z，C_φ及D_x，D_y，D_z为系数项，而r_x，r_y，r_z，r_φ为源项，待求未知数个数为

(24)

对于线性方程组求解目前主要采用直接求解和迭代求解.迭代求解通常在Krylov子空间内进行.首先采用不完全LU分解、不完全Cholesky分解等预处理方法降低系数矩阵的条件数，然后采用共轭梯度法(CG)，双共轭梯度法(BiCG)，广义最小残量法(GMRES)，拟最小残差法(QMR)等迭代方法进行求解.该算法需要的内存小，计算速度快，适合求解场源较少的正反演问题.直接求解方法一般直接对系数矩阵进行LU分解，然后利用各种成熟的求解库MUMPS、PARDISO等进行求解.该方法所需内存大、计算速度慢，但特别适合多源问题.本文离散线性方程组(23)的系数矩阵为大型对称稀疏矩阵，由于仅考虑单一发射源问题，我们选用Jacobi迭代的QMR算法进行求解.

4 数值模拟与结果分析

为了验证本文的频率域三维电磁正演算法的可行性和有效性，我们利用本文算法结果和一维正演结果进行对比，研究不同场源激发条件下频率域视电阻率和相位响应特征.

4.1 算法验证

选取水平层状三层H模型(如图 5a).第一层厚度500 m，电阻率100 Ωm；第二层厚度1000 m，电阻率10 Ωm；第三层电阻率为1000 Ωm.求解区域大小为10000 m×10000 m×10000 m，网格剖分为100×100×100，分别向外延拓16层，延拓指数为1.5.

首先进行平面波场X极化和Y极化模式下的电磁场三维模拟，然后计算大地电磁阻抗视电阻率和相位.计算频率为1～1000 Hz.图 5b为水平层状介质三维和一维正演结果响应对比.由图可见三维模拟的视电阻率和相位与一维结果完全吻合，视电阻率相对误差最大0.3%，相位最大相对误差2.9%，说明本文算法对大地电磁法三维模拟具有较高精度.

为验证本文算法的广谱有效性，我们进一步对有限长接地线源的CSAMT进行三维正演模拟.网格参数设置和上面模型一致.假设接地长导线发射源沿x方向置于地面，长度为2000 m，中心点相对投影(0，0，0)，发射频率为10 Hz.接收点设在位于发射源中垂线的正向y轴上，点距间隔为100 m.图 6为层状介质CSAMT三维和一维模拟结果对比.从图可以看出，E_x分量的实部和虚部三维模拟结果几乎与一维结果重合，最大相对误差不超过3%；H_z分量实部和虚部也基本与一维正演结果吻合，除了H_z分量实部受场源的畸变影响在发射源附近两个测点相对误差达到10%，其他测点相对误差均小于5%.考虑到CSAMT通常在远区或过渡区进行观测，在场源附近发生的畸变并不会影响数据解释精度.

最后对大回线源激发条件下的三维电磁场进行正演模拟.网格参数设置和上面模型一致.设置发射源为3000 m×3000 m方形回线，置于地面，中心点相对投影(0，0，0)，发射频率为10 Hz.接收点位于发射源中垂线的正向y轴上，点距100 m.图 6为层状介质大回线激发下的电磁场三维和一维模拟结果对比.从图中可以看出，E_x分量的实部和虚部三维模拟和一维模拟结果重合，最大相对误差不超过3%；H_z分量实部和虚部也基本与一维正演结果吻合，除了H_z分量实部受场源的畸变影响在靠近发射源1.5 km附近相对误差可达10%，但其余各测点相对误差均小于5%.因此，本文算法对大回线激发下的电磁场三维模拟具有较高精度.

4.2 典型地电模型

通过水平三层模型不同类型场源的数值模拟结果与解析解的对比，证明了本文算法适合各类场源的频域电磁法三维正演，并且具有较高精度.下面，我们研究不同类型场源复杂模型的频率域电磁响应特征，特别是进行视电阻率和相位响应结果的对比.如图 8，我们设计一个低阻异常体模型大小为800 m×800 m×500 m，顶部埋深为500 m，异常体电阻率为10 Ωm，围岩电阻率为100 Ωm.将10000 m×10000 m×10000 m计算区域剖分为100×100×100个单元，沿x、y、z三个方向的扩展层数为16，空气扩展层为12.将异常体剖分成8×8×5个单元，计算频率为1～1000 Hz.人工源频率域电磁法采用赤道观测装置，有线长接地导线长度为2000 m，磁偶源为600 m×600 m回线，收发距都为6 km.

图 9为不同场源激发下低阻模型中心点视电阻率和相位结果对比.从图可以看出，在100 Hz以上四种不同类型场源(XY模式、YX模式、接地长导线、大回线磁偶源)的高频段视电阻率和相位曲线基本重合；随着频率降低XY模式、YX模式大地电磁视电阻率响应受三维异常体的影响出现分离，而相位曲线基本重合；其中XY模式受三维异常体影响要大于YX模式.对于接地长导线、大回线源，随频率降低观测场由远区逐渐进入近区，电阻率呈现快速上升，相位快速下降趋势.其中大回线磁偶源进入近区频率要高于接地长导线源，而且幅度大于长导线源.因此，实际电磁勘探工作中，为便于在远区观测通常采用接地长导线源作为发射源.

图 10为低阻异常体模型不同类型场源激发条件下视电阻率和相位拟断面图.由图可以看出，当频率高于100 Hz时，即处于远区条件下，CSAMT、回线源赤道装置的视电阻率响应基本相同，在异常体上方近地表位置有一个不太明显的小范围高阻区，下方为低阻异常区；低阻异常区的分布范围与模型中低阻异常体的实际宽度基本相同；相位特征也具有类似特征，呈现高值异常.随着频率降低，XY模式视电阻率呈现低阻视电阻率异常和高值相位异常；YX模式呈现低阻加两个对称高阻视电阻率异常和高值相位外加两个对称低值异常；然而接地长导线和大回线源的赤道装置受近场影响，呈现高阻假异常和低值相位假异常.其中大回线源进入近场频率要高于接地长导线，而且异常响应范围和幅值都要小.因此开发不同场源的频率域电磁三维反演算法可克服三维异常体和场源引起的虚假异常，有助于提高数据解释质量.

5 结论

本文从麦克斯韦方程组出发，在Lorenz规范条件下建立磁矢势和标势耦合对称方程.采用交错网格和有限体积法离散磁矢势和标势耦合方程，并将不同场源类型分解成一系列短导线电性源的组合；对离散线性方程组采用Jacobi迭代的预处理拟最小残差法(QMR)进行求解，我们成功实现任意场源下频率域三维电磁正演模拟.通过对水平层状模型不同类型场源的三维电磁模拟可以发现，本文开发的算法能够实现不同类型场源各种不同应用范围频率域三维电磁正演模拟，避免耗时的一次场计算，计算精度较高.由于采用特殊的源处理技术，使得频率域三维电磁正反演统一数据处理和解释平台成为可能.对不同场源激发下典型地电模型模拟发现，不同类型场源的视电阻率和相位特征和规律不同，在实际勘探中需要考虑场源的选择和设置，同时在电磁数据解释中应考虑场源的畸变效应.