大地电磁测深法(Magnetotelluric sounding, MT)是20世纪50年代初由A.N.Tikhonov和L. Cagnird分别提出的天然电磁场探测方法(Tikhonov,1950; Cagniard,1953)，它以天然的平面电磁波为场源，通过在地表观测相互正交的电磁场分量来获取地下电性构造信息.由于天然场中含有从高频到低频丰富的频率成分，而不同频率成分的电磁波具有不同的穿透(趋肤)深度，因而大地电磁法能达到测深的目的(陈乐寿，王光鄂, 1990).由于它具有探测深度大、不受高阻屏蔽的影响、对低阻层反应灵敏等优点，该方法始终为人们所关注和研究.早在20世纪70年代,国内、外就开始了大地电磁测深正、反演的研究工作.近年来,随着计算机技术的迅速发展，大地电磁测深二维正、反演算法的研究已经趋于成熟；二维反演方法已成为大地电磁测深生产与科学研究中主要的反演计算方法.而水平地形条件下大地电磁测深三维正演算法研究的进展，也使三维反演算法日渐趋于实际应用(Avdeev, 2005; Börner, 2010)；但是，要真正实现大地电磁测深三维反演的实用化，则必须解决复杂地形条件下大地电磁场的三维正演模拟和反演算法问题.

20世纪80年代以来，尽管国内、外许多学者在起伏地形条件下大地电磁场数值模拟问题上做了大量研究(Wannamaker et al., 1986; 徐世浙等, 1997; 陈伯舫等, 1998; Han et al., 2008)，但对于带地形三维数值模拟，尤其是带地形三维反演方法的研究仍然很少，仅有的几种方法也未见公开代码(Ledo, 2005).由此可见，目前大地电磁场数值模拟领域一个亟待解决的问题是如何在大地电磁测深三维正演模拟和反演计算中有效加入复杂地形要素(Siripunvaraporn, 2011).我国地域辽阔，东、西部高程变化大，地形极为复杂，无论是大地电磁测深长剖面观测或是区域性阵列观测都难免跨越多个地形单元；因此，为了提高中国大陆大地电磁测深数据的解释精度，在进行大地电磁测深数据反演时，也同样需要考虑复杂地形的影响；为此，我们也迫切需要研究和发展一种高效的，能够进行带地形要素三维正、反演的大地电磁测深数值模拟软件.

在本项研究中，我们实现了基于交错网格剖分的有限差分算法的大地电磁测深带地形三维正演，以及基于非线性共轭梯度算法的带地形三维反演方法，本文将具体阐述所取得的研究成果.

众所周知，地球物理数值模拟的正演问题是反演方法的基础和前提，反演系统是否准确和有效很大程度上取决于正演方法的准确度和计算效率.对于大地电磁测深三维正演数值模拟算法，前人已经进行了大量的工作，常见的方法有积分方程法(IE)(Ting and Hohmann, 1981; Hohmann, 1983),有限差分法(FD) (Mackie et al., 1994; Smith 1996; Newman and Alumbaugh,1997;谭捍东等,2003), 有限单元法(FE)(Pridmore et al., 1981; Mitsuhata and Uchida, 2004; Nam et al., 2007),边界元法(BE)(徐世浙等, 1997; 阮百尧, 王友学, 2005)等.有限差分法因其具有计算精度和效率较高，不需要专门的网格离散化工具等优点，是目前大地电磁测深正、反演方法采用较多的数值模拟方法.

本研究中使用的正演方法为基于交错网格剖分的有限差分方法(谭捍东等, 2003).在大地电磁测深法研究的频率范围内，通常位移电流的作用可以被忽略.取电磁场的时谐因子为e^－iωt，有麦克斯韦方程组的微分形式为

其中μ₀是真空磁导率，σ是电导率，而ω为角频率，介质内没有电流源或磁场源.将方程组(2)式代入(1)式中，可以导出仅基于电场与仅基于磁场的麦克斯韦方程分别为

要计算出感兴趣的异常体中的电场 E 或磁场 H ，就需要将上述方程离散化；在有限差分法中，通常把求解模型域划分为若干长方体差分网格，用有 限个网格节点代替连续的模型域.在交错网格方法中，剖分单元中的电场和磁场并不取在网格单元的同一位置，而是取各网格单元边缘中点的电(或磁)场，与各网格面元中心的磁(或电)场相匹配，这样将可以最大限度的保证电磁场分布能够满足能量守恒定律(Mitsuhata and Uchida, 2004).

其中 C 为 E 或 H .假设模型域介质可以离散为 M =N_x×N_y×N_z个长方体网格单元，假设单元体中磁场为沿网格边元平均分布，而电场为沿网格面元均匀分布，如图1a所示，由平面波假设，大地电磁场的电磁波可以分解为相互垂直的线性偏振波的 组合的形式，设每一个网格有x,y,z三个方向的 电场/磁场物理量，单元体的长宽高分别为d_x(i)(i=1,2,3,…,N_x),d_y(j) (j=1,2,3,…,N_y)和d_z(k)(k=1,2,3,…,N_z)，单元体的电阻率为ρ=ρ(i,j,k). 对于感应磁场，(1)式的积分形式∮ E ·dl= ∫∫iωμ₀ H ·d S 的x,y,z

三方向分量可以按照图1b离散为

图1

Fig. 1

图1 交错网格有限差分单元网格离散化示意图 (a) 磁场沿网格边元平均分布，而电场沿网格面元均匀分布; (b) 感生磁场离散化示意图; (c) 感生电场离散化示意图. Fig.1 Discretion of stagger-grid finite element of a unit cell (a) Electric and magnetic field components are assumed to be edge-averaged and face-averaged, respectively; (b) Elementary loop of electrically induced magnetic field; (c) Elementary loop of magnetically induced electric field.

对于给定的频率，消去电场或磁场分量，这一问题就可以转换为求解线性方程组

其中， A 为3M×3M阶大型稀疏矩阵， x 为电/磁场的3M阶未知矢量，而 b 则为边界条件与场源的3M阶矢量；在这里，我们选择消去电场分量，求取磁场分量的线性方程组的解.与水平地面三维正演模拟相同，对于上边界的边界条件，其场值可以直接取采样点所在的源场值；而四个侧边界的边界条件则可以通过二维正演得到，在边界处取与三维介质边界层网格相同的电阻率和网格剖分建立二维模型，并针对不同的极化模式进行二维正演计算，即可得到四个侧边界上边界场值的分布；对于下边界的边界条件，一般可以设下边界以下半空间介质为均匀半空间或者层状介质，通过一维正演来进行计算(Mackie et al., 1994).由于正演问题中大型线性方程组的规模和复杂性，直接求解这一方程组几乎不可能，一般采用迭代计算的方法进行求解(Avdeev, 2005; Siripunvaraporn, 2011)；为此，我们使用双共轭梯度(Bi-Conjugate Gradient) 最优化方法对其进行迭代求解(Smith, 1996).在已知磁场的分布后，就可以根据法拉第定律通过式(2)求得电场的分布.

在进行地形的数值模拟时，最重要的问题是如何处理空气与地下介质的关系.而在交错网格有限差分法中，由于要求满足介质渐变假设，空气电导率一般不能设置为零，而是设为一个很小的非零电导率(Nam, et al., 2007).而这样可能会使得系统方程组产生严重的病态(ill-conditioned)(Mackie, et al., 1994; Mitsuhata and Uchida, 2004).当ω→0时，由于方程组(5)或(6)的右侧趋于零，而左侧保持不变，因而这一点在计算低频响应时尤为严重.本研究中，我们采用了散度校正(divergence correction)(Mackie, et al., 1994; Smith, 1996)来对这一问题进行处理.根据库仑规范(Coulumb gauge)，在地下介质中，显式地要求散度为零(即不存在场源)，对于电势e或磁势h分别有：

而在空气中，由于要求不存在电荷的积累，可以将设电势设为

利用求解出的势场，对电场或磁场进行散度校正：

在求解(15)的迭代计算过程中，每隔若干次迭代 进行一次场的散度计算，当散度不能视为零时(＞10^-10)，将使用(19)对磁场进行散度校正.由梯 度的旋度 Δ × Δφ=0，对于磁场的校正，有 Δ × H c= Δ ×( H － Δφ)= Δ × H =σ E ，这表明磁场的散度校正对于电场并无影响.类似的，对于电场我们有 Δ × E c= Δ ×( E － Δφ)= Δ × E =iωμ₀ H ，表明电场的散度校正对于磁场同样无影响.

此外，由于大地电磁场源的平面波假设要求离散网格的面元要垂直于波阵面，当使用有限差分法模拟地形时，很难直接模拟倾斜的地形表面；因而，多采用“台阶状”网格逼近地形的起伏变化(如图4c所示).当在地形起伏界面(空气-地面界面)上加大剖分网度，使其剖分足够精细时，大地电磁场有限差分模拟响应将与理论场值相近，甚至趋于一致(陈伯舫等, 1998).

具体计算时，对于构筑起伏地形介面的“台阶状”网格单元顶面的大地电磁场值，假设有自顶部向下第k个网格为空气-地形介面，可以通过直接计算第k个网格的波场来确定；若网格剖分足够精细，可以认为“台阶状”网格单元顶面的大地电磁场值与空间位置的关系近似为线性关系.当测点在网格顶面正中心时，测点磁场值可以简单的取相邻网格磁场的平均值；而当测点不在网格顶面正中心时，可以通过线性插值的方法计算出测点位置的大地电磁场值.由于我们之前假设单元体中磁场为沿网格边元平均分布，而电场为沿网格面元均匀分布，设测点位置为x,y,z，测点所在单元为i,j,k，则对于H_y只需在x方向上进行线性插值：

H_x则只需在y方向上进行线性插值：

H_z则需要进行x和y方向的双线性插值(注意，这样H_z的变化将不是严格线性的)，对于介面所在的第k个网格，有

而对于与介面上方临接的网格(k-1)，可以类似的得到

因此在第(i,j)个网格，空气-地下介面处(k和k-1之间)H_z可表示为

即可由磁场值计算电流密度水平x和y方向的电流为

分别对于J_x和J_y做Y和X方向的线性插值，可以得到测点不在网格中心时的水平电流为

再由(12—13)式可得到空气-地下介面处水平方向电场为

对于水平地形的三维电磁数值模拟，由于空气与地下介质的电阻率间断面的存在，再考虑到大地电磁测深的频率范围，一般强制设置E_z=0.而若考虑垂直方向电场E_z对介面处电流密度的影响，不仅垂直方向的电场将不可被忽略，水平磁场也会受到当前网格(i,j)及周边四个网格(i-1,j),(i+1,j),(i,j-1),(i,j+1)的垂直电场的影响发生改变，从而进一步改变水平方向的电场，如图2所示.

图2

Fig. 2

图2 交错网格有限差分法空气-地形介面电磁场离散化示意图 Fig.2 Discretionof electromagnetic fields at air-earth interface of stagger-grid finite difference method

对于第(i,j)网格所在单元，垂直方向的电流密度和电场为

对于第(i-1,j)网格所在单元，有：

类似的可以求出第(i+1,j)网格，第(i,j-1)网格以及第(i,j+1)网格所在单元的电流密度和电场强度值，垂直电场对x方向磁场的贡献为

而垂直电场对y方向磁场的贡献为

取之前得到的E_sx和H_sy的值，有x方向水平电场可以写成

类似由E_sy和H_sx的值的可以得到y方向水平电场为

计算出电/磁场值之后，即可继而计算出大地电磁阻抗张量响应(Mackie et al., 1994; 谭捍东等, 2003).设大地电磁场的平面偏振波沿x方向的偏振分量为H_x1,E_x1,H_y1,E_y1，y方向的偏振分量分别为H_x2,E_x2,H_y2,E_y2，则阻抗张量可以由下式进行计算(Mackie et al., 1993)：

总体而言，现有的多种反演方法都只是最优化方法的不同表述.水平地面三维大地电磁测深反演方法，包括高斯-牛顿法(GN)(Sasaki, 2001; Haber et al., 2004)，共轭梯度高斯-牛顿法(GN-CG) (Mackie et al., 1993; 吴小平，徐果明, 2000; 林昌洪等, 2008)， 拟牛顿法(QN)(Avdeev and Avdeeva, 2009)，非线性共轭梯度法(NLCG)(Mackie et al., 2001; Newman and Alumbaugh, 2002; 林昌洪等, 2008)以及OCCAM法(Siripunvaraporn et al., 2005)等，其最终目标都是寻找一个“最优”模型，它既可以拟合观测数据，又在地质上最为“现实”.两者都通过寻求各自目标函数的最小化这一过程来实现.其中，NLCG方法因其不必对最优化问题进行线性化，在计算过程中也不需显式的计算与存储Jacobian矩阵，在反演时间和内存的耗费上相对较低，成为近年来大地电磁测深二维和三维反演的主流方法之一.

对于大地电磁测深反演而言，非线性最优化问题的表述可以简单地写为

其中 d 表示数据向量，在这里即为实测阻抗张量(Z)；而 m 表示模型向量，在这里即为电导率(log₁₀σ)； e 为观测误差向量； F 则表示上一节给出的带地形大地电磁场三维正演算子.对于非约束最优化数值模拟问题有目标函数为

目标函数的第一部分ψ_d表征数据的拟合情况；第二部分ψ_m则为模型的一阶或二阶光滑度，其中 C _d为控制数据误差的权重的协方差矩阵，λ则是平衡数据拟合和模型光滑度的正则化因子.而对于某个给定的模型 m ，目标函数的导数形式 Δ Ψ可以写为(Newman and Alumbaugh, 2002)

其中，第二部分ψ_m，即关于模型第k分量的偏导数则可以简单的写为

而第一部分ψ_d，即关于模型第k分量的偏导数为

注意这里计算目标函数的导数 g 不需要显式的计算和存储雅可比矩阵，而只需要其转秩和(加权的)误差分量的乘积 J ( m ) ^CT ^－1_ed 即可(林昌洪等, 2008).其中 J 为雅可比矩阵，e_n=Z_n－ F ( m )为第n个数据误差， ε n为协方差矩阵 C _d对于第n个数据的分量，*为复共轭.

对于 Z /m_k，具体到 Z 的某个分量(Z_xx,Z_xy,Z_yx,Z_yy)，根据式(42—45)，我们可以得到其相对m_k的导数，以Z_xy为例，有

其中a=H_x1H_y2－H_x2H_y1(Newman and Alumbaugh, 2002).注意到式(51)中包含电磁场分量关于模型分量的偏导数，因此要计算 Z /m_k还需得到这些偏导数的值.对于给定位置j的测站，有第i个分量的磁场/电场相对于第k个模型分量的偏导数为

其中A即为式(15)中的线性关系，^h g ^T_ji为在给定测点j处的磁场的第i个分量(H_x,H_y)基于该给定磁场 H 的插值向量(interpolation vector)，而^e g ^T_ji为在给定测点j处的电场的第i个分量(E_x,E_y)基于该给定磁场 H 的插值向量.将式(52—53)代入到式(51)中，可以得到(林昌洪等, 2008)

其中：

类似的可以求出 Z _xx/m_k, Z _yx/m_k, Z _yy/m_k的偏导数(Newman and Alumbaugh, 2002)，由于篇幅所限，这里不再赘述.于是(50)式即可表示为

设

进而如果设

又由 A = A T，我们只需要求解如下线性方程，相当于两次正演计算：

就可以得到¹ v 和² v ，进而得到目标函数的导数 Δ Ψ(Newman and Alumbaugh, 2002).

若由任意给定的模型初始值 m ₀开始，第k次迭代的模型的改变量可以由

进行迭代计算，其中 α k为第k次迭代的搜索步长，可以通过线搜索(line search)得到；而 d k为第k次迭代的搜索方向，可以递归的由

进行计算，其中 g _k= Δ F( m _k)T，即为函数 f 对于 x _k的导数向量.观察 d ₀可以发现，对于共轭梯度法，第一次进行最优化搜索的方向就是- g ₀，即为函数值下降最快的方向，这一点和所谓最速下降法或梯度法是相同的.式中β_k为共轭梯度方向更新因子；而由前人各自发展的不同的因子就确定了一系列不同的共轭梯度方法.在大地电磁数值模拟中，使用较多的是Hestenes and Stiefel(1952)的经典CG方法 β^HS_k= g ^T_k+1( g _k+1－ g _k) d T_k( g _k+1－ g _k) (Mackie et al., 1993)，或者Polak and Ribiere(1969)的方法β^PR_k= g ^T_k M _k( g _k－ g _k－1) ( g ^T_k－1 M _k－1－ g _k－1) (Mackie et al., 2001).在这里，经过对比试验，我们使用了Dai and Yuan的方法 β^DY_k= g ^T_k M _k( g _k) g ^T_k M _k－1( g _k－ g _k－1) 以提高全局收敛性(1999).

对于给定的搜索方向 d _k， 我们可以通过线搜索来得到模型搜索步长 α _k， 有最优化问题

即寻找目标函数在 α _k 所在的超平面上的最小值.实际上，这就变成了一个单分量 α _k 函数的最优化问 题，对于函数 f ( x + αd )，将其进行泰勒级数分解，有

因此有

其中 f′ ( x )为梯度矩阵，而 f″ ( x )为海森矩阵(Hessian)(Ypma,1995).为了避免计算海森矩阵，我们将二次导数近似的表示为一次导数在给定的邻域 α =(0,σ)的微分，有

这样，利用泰勒级数展开，有

同样的当上式趋于零时， f ( x + αd )达到最小值(Press et al., 1992)，这时有：

需要注意的是，以上的方法只能分辨出二阶导数为零的情况，而并不能分辨出最大值和最小值(Presset al., 1992).因此，在NLCG方法(或者任意一个依赖于线性搜索的方法)中，初值的选择有时可能会影响收敛速度，甚至出现不收敛(有可能收敛到局部最大而不是最小值)的情况；当搜索不到最小值时，我们将把线性搜索的方向设为- g _k， 重置反演线性搜索方向.

一旦计算出了 α _k和 d _k，下个模型 m _k+1 就可以不断的通过(63)式以迭代的方式计算出来，直到模型拟合差达到期望值或者迭代次数超出最大次数为止.对于带地形反演，我们不希望空气层的电阻率随 着反演的进行而发生变化，因此需要锁定地形中空气层所在的电阻率， 为此我们引入模型控制矩阵 C _m为

从而(40)式可表示为 m _k+1= m _k+ C _m· α _k d _k，这样，模型改变量将不会应用到空气层当中.

由于地球物理问题中的线性系统常常是大型稀疏而病态的(Smith, 1996)，为改善矩阵条件数(condition number)，提高反演的收敛速度，我们引入预条件矩阵 M .对于线性问题 Ax = b ，我们可以通过解

间接的进行计算，如果有COND( M ^－1)≤COND( A )，我们就可以比原系统更快的收敛到所求解的邻域.对于第k次迭代，对于目标函数(47)式而言，效率最高的 M 自然应该是 A ^－1=[λ^－1( d － F ［ m _k])］T C _d( d － F [ m _k])+ L T L ］^－1，然而，这样将导致求解 M ^－1 Ax = M ^－1 b 和 Ax = b 完全一样而失去了 意义.在 M 的选择上，一种常见的做法是取目标函数的海森矩阵，即二次导数矩阵的近似.然而，由于海森矩阵的计算量过大，Newman等人(2002)提出了仅仅使用海森矩阵的对角元素作为预调节矩阵的方法，而其对角元素可以依靠目标函数的导数和上一次共轭梯度迭代的结果迭代地给出

反演的主要流程如图3所示.

本项研究实现了基于有限差分法的带地形三维大地电磁测深反演算法，并使用一系列三维理论模型及由正演理论响应分别构筑的模拟实测数据，验证了所研发的带地形三维正、反演算法的准确性和有效性.

文中给出的正演模型中设计的地形条件虽然较为极端，但仍然可以有效的证明在地形起伏较大的情况下，普通水平地形的三维反演方法在处理带地形数据时，虽然也可以得到异常体的主要特征，但其形态与产状都受到了较强的畸变，可能是由于无法模拟地形产生的复杂电性畸变，反演程序对数据进行了过度拟合，从而产生了一系列明显的虚假异常；而带地形三维反演方法可以得到比较接近地下真实导电性结构的反演模型.因此，在地形起伏较大的情况下，应优先使用带地形方法进行大地电磁反演解释；当忽略地形进行反演时，需要依实际情况进行谨慎处理和分析，以免得出不当解释和推断.此外，由 于模型需要模拟地形起伏，而有限差分方法又要求在空气界面处剖分足够精细，这样不可避免的会造成模型网格数出现较大的增加，因此内存耗费以及反演CPU时间增长等问题难以避免，进一步的研 究将考虑采用并行算法以解决内存及计算效率问题.