0 引 言

1954年频率域航空电磁法最早在加拿大首次用于生产，并在铜矿床的勘探中发挥了重要作用.经历了近70年的发展，目前国际上常用的频域航空电磁系统分为固定翼和吊舱式两种，其中吊舱式系统是应用较广的一种.国土资源部物化探研究所1975年研制成功Y5飞机HDY202单频航电系统；吉林大学(原长春地质学院)在七八十年代研制了时间域系统；70年代末国土资源部航遥中心从加拿大引进TRIDEM三频航电系统，2004年又引进IMPULSE吊舱式直升机频域电磁系统.IMPUSLE系统已经在小范围内进行了试验和生产工作.然而，制约该系统进一步应用的瓶颈就是缺乏对实测数据高精度的解释，特别是能够针对复杂地质条件下的三维航空电磁响应的数值模拟方法还存在很多问题：

(1)频率范围跨度较大，一般最低频率几赫兹，最高频率可达几十万赫兹.低频计算中位移电流的忽略将使高频计算变得不稳定性.

(2)在低频电磁响应的计算中，存在低频迭代发散的问题，需要进行散度校正.

(3)三维计算模型较大，受计算机容量的限制.同时在计算磁偶极子场源附近的响应时，极大的场值变化梯度就要求将网格加密；另一方面，由于计算机有限的存储，这会使得加密后的计算模型缩小，从而达不到可靠精度下的应用要求.

(4)航空电磁法数值计算不仅需要考虑地下介质，地面和空中的电磁影响也需要考虑，譬如飞机本底响应.这些都给计算带来一定问题.

李小康(2011)对基于MPI的频率域航空电磁法有限元二维正演并行计算进行了研究，谭林(2010)利用有限元方法研究开发了频率域航空电磁法2.5维模拟软件.前者提高了正演计算的速度，后者则丰富了解释的维度.但是对于当今航空电磁面积性的勘探，二维或者2.5维的解释仍旧很难达到高精度勘探的目的.同时对复杂地质体的数值模拟方法来说，有限差分法较有限元方法更加简单快捷.Newman和Alumbaugh(1995)，Sasaki 等(2003)利用三维有限差分法模拟了二维简单地形的频域航空电磁响应，表明有限差分方法在解决频域航空电磁问题方面比有限元方法更加方便简洁易行.然而Newman等研究所选的频率最低只有900 Hz，对于实际应用中的更低频率并未给出计算结果，而且在计算过程中都提及低频情形下的方程数值求解的发散的情况，并未给出讨论.Mackie 等(1994)讨论了低频发散的原因，并提出相应的校正方法——散度校正(Mackie et al., 1994；沈金松，2003; 陈辉等，2011).该方法在迭代的过程中还需要再进行大型稀疏线性代数方程组的求解，从而增加了迭代计算的时间.为了避免这种附加的迭代，本文在Newman三维有限差分方法的基础上，将大型稀疏复线性代数方程组的求解方法作了修改：将方程组两边实虚部分离，得到关于实数的大型稀疏线性代数方程组，利用稳定的双共轭梯度法求解方程组.计算结果表明该方法不仅可以提高方程求解速度，而且保证方程在较低频率下解的收敛性，满足航空电磁法宽频带的要求.最后根据航空电磁法本身特点，提出了将初值进行优化的双重网格迭代方法，可以提高50%左右的计算速度，从而为三维反演打好基础.

1 三维频率域航空电磁有限差分方法

假设电磁波传播过程中的时谐因子为e^iωt，其中i=，并且频率f=ω/2π.在低频电磁法勘探中位移电流一般认为可以被忽略，麦克斯韦方程组为

为了避免场源附近场值变化梯度过大的问题，选择麦克斯韦方程组的散射场计算，公式为

利用交错网格有限差分方法可以求(8)式，继而得到地下介质的响应.在将其离散的过程中，介质被离散为一个个长方体单元，每个长方体单元代表一个电性体.本文将电场分量定义在单元体边的中点，而磁场分量定义在单元面的中心，如图 1.散射场可以根据不等步长的网格计算出来.在使用第一边界条件时，为了减少网格边界的反射，可以在靠近网格边界的地方采用更大的步长，使模型边界扩大到距离目标区域无穷远的地方.另外，由于有限差分计算的是一个平均化了的一个点，同样我们需要对某一点所代表的电阻率或电导率取平均.在计算电阻率时，同一点的电阻率在各个方向按照相邻单元在该方向截面大小进行加权平均(Alumbaugh et al., 1996).

图 1

Fig. 1

图 1 交错网格中电磁场各分量分布 Fig. 1 Components of the electromagnetic field distribution in staggered grid

利用(8)式，按照图 1所示的有限差分方法，可得Maxwell方程组的离散结果(见附录).如若求得各个单元电场值，可以用插值的方法求得空间任意点电场，同时利用(3)式中的电场旋度方程可以求得任意点的磁场值.

2 大型稀疏线性代数方程组的求解

将附录中离散格式整理得到一个形如

A x= B

的方程组，其中A是大型稀疏复对称矩阵

(9)式矩阵中各个子元 M _ij(i，j=x，y)同样也是稀疏对称矩阵，其中 M _xj为附录中求 x分量时的系数，M _xy对应于 C_6-9，M _xz对应于 C_10-13，而M _xx对应于 C _1-5，其余两个分量按照此规则依次对应. x为欲求解场值，B 为和背景场相关的向量.由离散结果可知，对于这样的大型方程组的问题往往利用迭代法来求解，主要计算方法有基于Lanczos分解原理(Lanczos，1950; Paig，1976; Druskin and Knizhnerman， 1994)的双共轭梯度法和基于剩余模最小的Amoldi原理(Fujino and Zhang， 1996)最小残差方法.本文比较了两种方法各自具有代表性的迭代算法：带有稳定因子的双共轭梯度法BiCGStab2(St2)(张继锋等，1996；王志刚等，2006；张永杰和孙秦，2007)和拟最小残差法(QMR)(李欣和朱景福，2009).

为了提高对大型稀疏复线性代数方程组的求解速度(谭捍东等，2003；付长民等，2009; 邓居智，2011；邓居智等，2011)，本文将复数矩阵化为实数矩阵，把方程组中每个矩阵分成实部和虚部两部分为

在利用迭代方法求解方程组时，往往需要给出一个迭代初始值，选择一个好的初始值可以大大缩短迭代次数和迭代时间.由于频域航空电磁法对每个频率进行多点剖面测量，飞机在飞行过程中每隔2~3 m就进行采样测量，在实际模拟过程中若采样较密集，则可以将前一次计算场值的空间分布作为正在进行迭代计算的初始值，进而实现初始值的优化处理，同时将前一场源点在求解方程组中的搜索方向保存至当前计算点.

另外，计算精度的提高还需要网格数的加密，也就是计算步长的减小.在计算过程中，对于固定大小的区域，划分网格数的多少将直接影响计算速度和计算结果.网格数越少，计算越快，但是计算精度越低；反之则反之.这又构成计算中的一个矛盾，计算步长和计算时间之间的矛盾.这一矛盾的解决将直接影响最终的计算效果和效率.本文在此提出解决该矛盾的初步想法：如果在初始迭代时，用较粗的网格，经过少量的迭代后即可达到一定精度；然后将此精度下的结果插值后作为初始值代入精细网格下进行迭代求解，这样即可处理好以上矛盾.

3 数值计算准确性检验及算例

为了证明方法以及按照该方法所编程序的正确性，本文将交错网格有限差分的结果和解析解进行了对比.选择纳比吉安(1992)所给出的经典模型，将磁偶极子置于电阻率为100 Ωm的均匀大地表面，观测点距磁偶极子100 m，磁偶极子磁矩为1 A·m².

由图 2和3可知用有限差分方法得到的数值解和解析解两者很好地拟合，而且由St2方法求得的Hx分量比QMR求得的结果在低频时拟合效果更好，Hz分量在较高频率(曲线变化尖端处)拟合效果也是前者优于后者.

图 2

Fig. 2

图 2 稳定双共轭梯度法与解析解求解结果比较图 Fig. 2 Comparison of results between St2 and analytical solution

图 3

Fig. 3

图 3 QMR方法和解析解求解结果比较图 Fig. 3 Comparison of results between QMR and analytical solution

为了进一步验证程序的准确性和稳定性，下面用三种不同步长下的模型去计算同一地质情况下的航空电磁响应.图 4为均匀大地模型示意图，三种网格划分依次为： 41×41×41，50×50×50，60×60×60，分别记为情况1、2、3，计算选用频率4000 Hz，磁偶极矩为1 A·m².计算时将发射点置于400 m处，在发射点两边对称地每隔10 m取一记录点，从340 m到460 m形成12个记录点.发射装置分垂直磁偶极子和水平磁偶极子两种，图 5和图 6即为垂直磁偶极子和水平磁偶极子共发射点记录在三种情况下的结果.从最后结果可以看到，不管是水平磁偶极子还是垂直磁偶极子在三种网格下的结果都取得很好的一致性.而且垂直磁偶极子Z分量关于发射点对应坐标镜像对称，而X分量关于该坐标中心对称，水平磁偶极子恰好与之相反.

图 4

Fig. 4

图 4 均匀大地模型，其中探头距离地面保持20 m， 大地电阻率为100 Ωm，接受装置沿测线对称地 设在400 m处的发射装置两边 Fig. 4 Half space model，the bird is kept 20 m above the earth surface，the resistivity of the earth is 100 Ωm， the transmitter is hold at 400 m with receivers set symmetrically at both sides of the transmitter along X direction

图 5

Fig. 5

图 5 垂直磁偶极子共发射点记录，其中菱形为情况1，圆圈为情况2，方形为情况3， 上面两幅为磁场x分量实虚部，下面两幅为磁场z分量实虚部 Fig. 5 Responses of the Vertical Magnetic Dipole(VMD)with the same transmitter above earth，diamond st and s for case 1，circle st and s for case 2，square st and s for case3，up two figures are the records of Hx and the down two figures are Hz’s responses

图 6

Fig. 6

图 6 水平磁偶极子共发射点记录，其中菱形为情况1，圆圈为情况2，方形为情况3， 上面两幅为磁场x分量实虚部，下面两幅为磁场z分量实虚部 Fig. 6 Responses of the Horizontal Magnetic Dipole(HMD)with the same transmitter above earth diamond st and s for case 1，circle st and s for case 2，square st and s for case 3，up two figures are the records of Hx and the down two figures are Hz’s responses

在图 4 所示的均匀半空间模型中加一个低阻体，其宽40 m，顶面埋深20 m，底面埋深125 m，走向延伸400 m，电阻率为1 Ωm.下面比较在迭代过程中初值优化与否对计算时间的影响.表 1 所示是在4000 Hz下的计算比较结果，相对误差均控制在1%以内.由表可知初值优化的迭代方法均可以将迭代时间缩减50%以上.图 7展示的是在1 kHz发射频率时相对误差在10^-4下的迭代结果，其发射点置于45 m处，对于一般迭代过程，经过717次迭代后达到设定精度，耗时89 s；而先在40 m处设置一发射点，经过一般迭代后得到全区的场值，再将上述场值作为45 m处的初始值进行迭代，则经过498次即收敛，耗时54 s.以上计算均在普通笔记本电脑上完成.由此可见，对于45 m处的场值计算，时间缩短39%.

图 7

Fig. 7

图 7 初值优化前后迭代次数比较图 Fig. 7 Comparison of the iterative times of two methods，blue line represents general iteration and the red st and s for iteration with initial optimization

表 1(Table 1)

表 1 不同网格规模下的迭代次数比较 Table 1 Comparison of iterative times among different kinds of grids 网格数 41×41×41 60×60×60 72×72×72 普通迭代时间/迭代次数(s/次) 22/165 172/440 424/596 优化迭代时间/迭代次数(s/次) 9/62 61/141 205/291 相对缩短时间(%) 59% 64% 52%

表 1 不同网格规模下的迭代次数比较 Table 1 Comparison of iterative times among different kinds of grids

为了说明程序的适定性，下面给出了一个起伏地表下的垂直磁偶极子响应.并对三种情况下的大地电阻率：50 Ωm、 100 Ωm和500 Ωm，进行了讨论.计算区域为X方 向1000 m，Y方向和Z方向均为800 m，其中地形为三角形山峰，坡角23度，斜坡的一面在水平方向的投影为69 m，且两坡关于400 m处对称.为了和IMPULSE系统的参数尽量接近，发射频率4650 Hz，收发距保持6.5 m，探头离地面高度适中设为30 m.最终图 8呈现的是对一次场归一化了的结果.从最终结果可以看出，垂直磁偶极子的地形响应和地形起伏存在镜像对称的特点，而且大地的电导率越高，响应越大.这说明航空电磁法同其他的电磁方法一样对低阻体比较敏感，利用该法方法可以探测地下隐伏低阻体，同时提醒航空电磁工作者地形影响是不可忽略的.

图 8

Fig. 8

图 8 不同大地电阻率下三角山模型的 垂直磁偶极子响应比较图 Fig. 8 Comparison of the responses of VMP with the triangle-hill model of different resistivity

(1)利用本文中提出的将大型复线性代数方程组转换为求解实方程组的稳定双共轭梯度方法较直接求解复方程组的QMR方法，可以更好满足三维航空电磁宽频带数值模拟的要求，保证在低频情形下的收敛以及高频下的计算精度.

(2)依据航空电磁法本身测量特点，提出的初值优化的思想，以及双重网格迭代的方法，在合适条件下可以明显提高迭代速度，缩短整体计算处理时间.

(3)需要注意的是该方法还受地下介质复杂程度的影响.地下介质越复杂，则地下各处场值差异越大，即使在两个较近的相邻发射点上，整个空间场值对应点的分布前后也会有差异，从而使初值优化的效率降低.

致 谢 感谢吉林大学地球探测科学与技术学院曾昭发教授、李桐林教授以及殷长春教授在课题完成过程中给予的指导和帮助，以及范翠松博士在程序编写中提供的建议.

附录

对于X方向，可形成如下离散形式：

C₁ E_x(i+1/2,j,k-1)+C₂E_x(i+1/2,j-1,k)+C₃E_x(i+1/2,j,k)+C₄E_x(i+1/2,j+1,k)+C₅E_x(i+1/2,j,k+1)+C₆E_y(i,j-1/2,k)+C₇E_y(i+1,j-1/2,k)+C₈E_y(i,j+1/2,k)+C₉E_y(i+1,j+1/2,k)+C₁₀E_z(i,j,k-1/2)+C₁₁E_z(i+1,j,k-1/2)+C₁₂E_z(i,j,k+1/2)+C₁₃E_z(i+1,j,k+1/2)=C_oE_xp(i+1/2,j,k),

对于Y方向，可形成如下离散形式：

C₁₀=C₁₃=Δx′ _{i ，}

C₁₁=C₁₂=-Δx′ _{i ，}

C₀=-iωμ₀(σ(i,j+1/2,k)-σ_p(i,j+1/2,k))Δx′ _{iΔyjΔz′ k ，}

C₁ E_x(i-1/2,j,k)+C₂E_x(i+1/2,j,k)+C₃E_x(i-1/2,j,k+1)+C₄E_x(i+1/2,j,k+1)+C₅E_y(i,j-1/2,k)+C₆E_y(i,j+1/2,k)+C₇E_y(i,j-1/2,k+1)+C₈E_y(i,j+1/2,k+1)+C₉E_z(i,j-1,k+1/2)+C₁₀E_z(i-1,j,k+1/2)+C₁₁E_z(i,j,k+1/2)+C₁₂E_z(i+1,j,k+1/2)+C₁₃E_z(i,j+1,k+1/2)=C_oE_zp(i,j,k+1/2) ，