电磁法是一种广泛应用于矿产资源勘查、水文地质和环境地质调查等的地球物理勘探方法，依据响应性质可分为频率域电磁法和时间域电磁法.瞬变电磁法(Transient Electromagnetic Method)属于时间域电磁法，经过60多年的发展，其在理论研究和应用水平方面均取得了较大进步(薛国强等， 2007，2008；闫述等，2011).但由于多维正反演耗时久、对硬件要求高，因此前人研究不多，野外数据资料的解释仍停留在一维水平(熊彬和罗延钟，2006).2.5维问题，与三维相比较，大大减少了计算量，与纯二维比较又能更接近野外实际地质情况.

对于瞬变电磁法高维正演，常采用的方法为积分方程法( San Filipo and Hohmann， 1985；Xiong，1992)、有限单元法(Sugeng et al., 1993； Sugeng，1998；王华军和罗延钟，2003；熊彬，2005；周俊杰，2011)和时域有限差分法(Oristaglio and Hohmann， 1984；Adhidjaja and Hohmann， 1989；Leppin，1992；Wang and Hohmann， 1993；闫述等， 2001，2002；岳建华等，2007;闫述等，2009).积分方程法主要适用于三维正演且只适合于简单模型，这种算法首先频率域内利用基于张量格林函数的体积分方程法计算电磁场各分量的频率域响应，之后利用快速数字滤波技术将计算结果转换到时间域；有限单元法主要是在频率域中进行异常场的求解，之后通过傅氏变换转换到时间域，这种方法适用于二维、三维、2.5维问题，但理论较有限差分法复杂；时域有限差分法直接把含时间变量的Maxwell旋度方程或波动方程(扩散方程)，在Yee氏网格或其他离散网格中转换成差分方程.每个网格点上的场值均与其相邻的场值有关.在每个时间步计算网格空间各点上的场值，随着时间的推进，便能直接模拟电磁波的传播及其与物体的相互作用.由于时域有限差分法的直接出发点是概括电磁场普遍规律的麦克斯韦方程，因此这一方法具有广泛的适用性.在时间域直接求解，避免了更多数学工具的使用，使得该算法简单易于实现(王长清和祝西里，1994)，且能够直观的反应瞬变场随时间的扩散情况，因此在瞬变电磁的二维、三维正演中得到了广泛的应用.本文尝试利用有限差分方法进行瞬变电磁法2.5维正演模拟，通过与解析解的对比，证明了算法的正确性与可靠性.

瞬变电磁场在各项同性的三维无源介质中近似为扩散方程

在直角坐标系中磁场强度各分量同样满足上述扩散方程，有

式中μ，σ分别为地下介质的磁导率与电导率，H_z为磁场强度垂直分量.

瞬变电磁法观测的数据主要为垂直感应电动势，可通过磁场对时间的导数来进行求解

式中 E _MF为感应电动势，μ为磁导率，n为接收线圈匝数，R_x为接收线圈的面积.因此瞬变电磁法2.5维正演将从(2)式出发.假设地下介质的电导率沿y方向均匀分布，则地下介质为二维结构，而场源为三维场源，则利用傅里叶变换可以将三维问题转换为二维问题.现引入如下傅氏变换

将式中磁场强度沿y方向做傅氏变换可以得到

为了简单起见，上式和下文将用H来表示磁场的垂直分量H_z，上标波浪线表示波数域场值.

2.5维有限差分问题中，空间网格的离散与二维相同，只是图中磁场值表示的是波数域的磁场值.如图 1，下标(i，j)表示网格点的空间坐标，σ_i，j表示各单元内的电导率值.Δx和Δz为网格空间步长，即Δx_i=x_i+1－x_i，Δz_j=z_j+1－z_j，Δt为时间步长，n为时间步长的个数，则Δt_n=t_n+1－t_n.

图 1

Fig. 1

图 1 2.5维模型网格离散及波数域Hz分布 Fig. 1 Grid discretization of 2.5-D model and the distribution of Hz in wavenumber domain

空间项的离散，对矩形ABCD积分并应用格林定理得

其中σ _i，j为(x_i，z_j)节点周围四个单元电导率的面积加权平均值，即

对时间导数的最简单近似为t=nΔt和t= n+1 Δt两个时刻的向前差分，即

这样利用(7)(8)和(9)式便可以推导出均匀网格情况下，按时间向前推算扩散方程的显示欧拉方法，即

其中r_i，j= Δt μσ _i，jΔ²，称为局部网度.

不难证明当网度或小于1/(4+Δ²k²_y)时，欧拉法是稳定的(米萨克N纳比吉安等，1992).故对于欧拉法最大的时间步长为

其中min(σ _i，j)表示加权平均电导率的最小值.

针对不同的电导率，电磁场的扩散情况也不同，这对欧拉法是一个难题.在一个典型的地电断面中，min(σ _i，j)可能低到10^-3 S/m.若网格大小为10 m，式给出的最大时间步长约为0.025 μs.

现考虑较利用具有更高精度的中心差分来离散时间导数

此外利用如下近似

最终得磁场的差分格式如下

应用Dufort-Frankel法得到差分格式(14)是对任何时间步长无条件稳定的显示差分格式.Oristaglio and Hohmann(1984)给出的最大实用时间步长为

在电磁法模拟中取t=0.1 ms较为合理，因为它是大多数瞬变电磁系统开始记录的时间.当σ_min为10^-3 S/m，网格大小为10 m时，(15)式给出的最大时间步长为1.77 μs.

根据三维有限差分计算的经验将初始时间选取为(Oristaglio and Hohmann(1984))

其中μ和σ分别为地表的磁导率和电导率，Δ₁为地表第一层的网格间隔.

本文采用矩形回线源作为激发源，将矩形回线源在均匀半空间上形成的波数域瞬变电磁场作为TEM2.5维有限差分正演模拟的初始条件.均匀半空间上空矩形回线源形成的拉氏域波数域垂直磁场表达式为(王华军和罗延钟，2003；周俊杰，2011)

空气层(z≤z₀)中

地下介质层(z>z₀)中

在勘查地球物理中，观测的是负阶跃函数的响应，及观测稳定电流切断后场的衰减，负阶跃响应可以表示为(米萨克N纳比吉安，1992)

因此将所求的的频率域的响应通过逆拉氏变换到时间域，得到上阶跃电流的时间域响应，然后再利用(19)式来求取所需要的下阶跃响应.

(17)及(18)计算式中含有正余弦项，王华军(2004)运用成熟的汉克尔变换理论导出了正余弦变换的数值滤波法，并给出了250点的滤波系数，并详细分析了该方法的计算精度.其表达式为

式中Δ为采样间隔，x为波数域变量，ccos(nΔ)、csin(nΔ)分别为为第n个余弦变换和正弦变换滤波系数.

逆拉氏变换的数值实现有多种办法，本文采用Gaver Stehfest(G-S)变换(罗延钟和昌彦君，2000).假设拉氏空间中的象函数为F(s)，实空间的象原函数为f(t)，则拉氏变换的表达式为

式中s为拉氏域变量，t为时间域变量，N为系数个数，V_m为GS变换系数

式中INT为取整函数，min为极小值函数.

本文根据前人的经验在计算中取n=12或n=16(周俊杰，2011).

运用基于三次样条差值的反傅氏变换，对于反傅氏变换公式，三次样条插值法求取的公式及步骤如下(徐士良，1995)：

设已知的函数为y=f(x)，在给定的节点处x₁2…n－1n上的函数值y₁，y₂，…，y_n－1，y_n.以及两个端点处的一阶导数值为y′ 1(x₁)与y′n(x_n)通过下式可以计算其他节点处的导数值为 a₁=0，b₁=y′(x₁)
h_j=x_j+1－x_j，j=1，2，…，n－1
a_j=h_j－1/(h_j－1+h_j)，j=2，3，…，n－1
β_j=3(1－a_j)(y_j－y_j－1)/h_j－1+a_j(y_j+1－y_j)/h_j，j=2，3，…，n－1
a_j=－a_j/ 2+(1－a_j)a_j－1，j=2，3，…，n－1
b_j= β_j－(1－a_j)b_j－1 / 2+(1－a_j)a_j－1，j=2，3，…，n－1
y′(x_j)=a_jy′(x_j+1)+b_j，j=n…1，n－2，…，2
利用数值导数y′(x_i)(i=1，2，…，n)，求解插值节点函数值公式为

在地下边界中，简单的假设位于边界上的波数域场量为零.这样做的前提条件是网格剖分要足够大，边界离源足够远(Oristaglio and Hohmann， 1984；Adhidjaja and Hohmann， 1989；Leppin，1992；Wang and Hohmann， 1993)

已知在似稳条件下，自由空间的磁感应强度满足拉普拉斯方程

对于所研究的2.5维问题，经过傅氏变换，(25)式在波数域中变为Helmholtz方程

各分量可表示为

在求上式各分量的解的过程中规定x有界，而z趋向于无限大.且地空边界处的磁场值 B ~和磁场的导数 B ~ /z已知.上述亥姆赫兹方程的解的一部分可表示为可用观测点 x₀，z₀ 对数奇异(Leppin，1992)

其中K₀为第二类0阶虚宗量贝塞尔函数.根据格林定理，有

其中F为区域F的边界，但不包括围绕观测点的一个小圆环，法向量 n 垂直F指向外，积分路径如图 2.

图 2

Fig. 2

图 2 积分路径图 Fig. 2 Integral path

对于Neuman问题，即地表处(z=0)磁场的导数 B ~ /z已知，此时定义格林函数形式为

则 G z =0.将图中积分路径的上半部分置于无穷远，则格林公式为

如果F₀的半径足够小，则上式中此边界上的积分项中的G为－ln k_yρ +K₀ 2k_y z₀，而G/ n 为1/ρ，dl为ρdφ.因ρ→0，则由上式可以得到计算地面以上部分磁场值的公式

而对于Dirichlet问题，即在地表处已知，此时格林函数定义为

则G=0，且

其中K₀为第二类0阶虚宗量贝塞尔函数，根据格林定理，最终得到

在有限差分计算过程中，已知地表的垂直磁场，利用(34)式便可以求出地表上空一个网格处磁场的垂直分量，进而进行差分迭代.

在计算中采用均匀磁导率μ₀，因此以上延拓公式均适用于磁场强度各分量H_x，H_y，H_z.上述积分式中包含的零阶和一阶修正的贝塞尔函数，其求解方法数学上已有大量的研究，本文直接采用其近似表达式(石庆冬，2001).

对于式中无穷积分的求取，我们采用下述处理方法(李世健，2012)

上式中第一部分的有限积分选择高斯-勒让德积分进行求解.又根据高-勒让德积分区间不能过大的性质，将大区间积分改成多个小区间积分之和，即

而第二部分采用高斯—拉盖尔积分进行求解.计算表明Δ和a这两个参数的选取直接影响着最终的计算精度和速度，综合考虑二者相互制约的特性，并通过大量的试算比较，计算中给定Δ=5，a=500-3000.

以上介绍了瞬变电磁法2.5维有限差分正演模拟的理论基础.利用Fortran语言编写了瞬变电磁法2.5维有限差分正演程序，现对数值模拟结果进行分析.

均匀半空间电阻率为100 Ωm，方形发射回线边长60 m，接收线圈等效面积50 m²，供电电流1 A，采用中心回线测量.有限差分最小时间步长1×10^－8，最大延迟时间为10 ms，波数选取17个.图 3a、b给出了两种非均匀剖分网格条件下有限差分数值解与解析解的对比曲线，图a中网格规模为181×80，最小网格2 m；图b中网格规模为109×55，最小网格5 m.图c为误差曲线可以看出数值解与解析解基本吻合，当最小网格为2 m时，在10^－6～10^－3s时间范围内最大相对误差不超过10%，10^－3～10^－2 s时间范围内最大相对误差为16%；当最小网格为5 m时，整体误差增加.当最小网格大小变化时对计算精度有所影响，且在时域有限差分计算中空间网格与时间步长的选取相互制约，影响着最终的计算精度.

图 3

Fig. 3

图 3 均匀半空间数值解与解析解对比(a)垂直感应电动势(Δmin=2 m);(b)垂直感应电动势(Δmin=5 m);(c)垂直感应电动势相对误差. Fig. 3 Comparison of analytical and numerical solutions for homogeneous half-space

首先是电阻率分层均匀的H型三层地电模型，模型参数为：ρ₁=100 Ωm、ρ₂=10 Ωm、ρ₃=100 Ωm；层厚度h₁=100 m、h₂=50 m，方形发射回线边长60 m，接收线圈等小面积50 m²，供电电流1 A，采用中心回线测量.有限差分网格规模181×80，最小网格2 m，最小时间步长1×10^－8，最大延迟时间为10 ms，波数选取17个.图 4为感应电动势响应曲线，可以看出有限差分数值解与解析解基本吻合，在如图所示的10^－6～10^－2 s时间范围内，计算的相对最大相对误差为10.91%，最小相对误差为0.36%，平均相对误差为3.91%.

图 4

Fig. 4

图 4 H型层状模型数值解与解析解对比(a)垂直感应电动势;(b)垂直感应电动势相对误差. Fig. 4 Comparison of analytical and numerical solutions for H-layered model

接下来是K型三层地电模型，模型参数为：ρ₁=100 Ωm、 ρ₂=1000 Ωm、ρ₃=100 Ωm；层厚度h₁=100 m、h₂=50 m，其他参数同上.图 5给出了感应电动曲线，可以看出有限差分数值解与1D数值解基本吻合，在如图所示的10^－6～10^－2 s时间范围内，计算的相对最大相对误差为15.80%，最小相对误差为0.0009%，平均相对误差为4.88%，可见K型模型计算精度较H型模型差.且可以看出瞬变电磁法对低阻体的反映较之高阻要好.

图 5

Fig. 5

图 5 K型层状模型数值解与解析解对比(a)垂直感应电动势;(b)垂直感应电动势相对误差. Fig. 5 Comparison of analytical and numerical solutions for K-layered model

误差来源分析：瞬变电磁法测量的感应电动势从早期到晚期跨越了八、九甚至更多个数量级，而且在时域有限差分法数值模拟中涉及到多个环节，各个环节误差积累.主要表现在：

(1)地空边界条件的处理中，插值与无穷积分的计算，引入误差，地下边界处场值设为0，计算区域有限，不能够真实体现场的传播特性，引入误差.

(2)均匀半空间的数值解本身存在误差，再作为有限差分计算的初始源，将误差带入，并随着时间的推进使误差叠加增大.

(3)反傅氏变换为无穷积分，我们只选取部分波数进行计算，带入误差.

(4)空间网格、时间步长的选取带来误差，当网格增大或时间步长增大时将严重影响最终的计算精度.

首先给出本文计算结果与周俊杰(2011)有限元结果的对比，计算中方形发射回线边长50 m，其他参数上文中H型层状模型相同.图 6a为有限元(FEM)数值解与1D解析解的对比曲线，b为有限差分(FDTD)数值解与1D解析解的对比曲线，c为二者与解析解相对误差曲线.图c中可以看出，两种方法在早期和中期均与解析解相吻合，而在10^－3～10^－2 s范围内本文有限差分计算明显优于文献中有限元方法结果.

图 6

Fig. 6

图 6 H型层状介质FEM与FDTD数值解与解析解的对比(a)FEM数值解与解析解对比;(b)FDTD数值解与解析解对比;(c)FDTD与FEM数值解与解析解相对误差. Fig. 6 Comparison of analytical and FEM and FDTD numerical solutions for H-layered model

此外，熊彬和罗延钟(2006)有限元计算结果指出，对于H型模型，计算平均相对误差为-5%；K型层状模型平均相对误差为5%，本文计算误差与之相当，进一步证明本文计算的可靠性.

为了进一步验证程序的正确性，下面给出两个局部异常体模型中心回线装置剖面响应特征.图 7为模型一直立低阻板状体示意图.如图所示，围岩电阻率ρ₀=100 Ωm，低阻直立板状体电阻率ρ₁=1 Ωm，顶端埋深20 m，异常体垂向和横向延伸分别为100 m和20 m.测量装置为中心回线装置，方形发射回线边长50 m，接收线圈等效面积50 m²，供电电流1 A.

图 7

Fig. 7

图 7 模型一示意图(中心回线装置) Fig. 7 Schematic diagram of model one(central loop unit)

图 8给出了上述模型21个测深点的中心回线瞬变电磁响应.图中横轴为测深点位置，坐标原点位于异常体中心正上方，测深点位置关于坐标原点对称分布；纵轴为感应电动势，从上到下对应多个采样时间道上的响应.从图中可以看到，早期道异常表现为单峰异常，极大值出现在了异常体中心正上方，随着时间的延迟，异常表现为双峰异常，两个极大值出现在直立板两侧，而直立板上方出现极小值，在晚期到异常又表现为单峰异常，极大值出现在板状体中心正上方.异常响应特征物理模拟(牛之琏，1992)结果以及其他数值模拟结果(熊彬和罗延钟，2006)规律一致.

图 8

Fig. 8

图 8 模型一瞬变电磁测深多测道剖面图 Fig. 8 Multi-chanel profile for TEM sounding of model one

图 9为模型二双低阻体模型示意图，围岩电阻率ρ₀=100 Ωm，等轴状低阻体电阻率分别为电阻率ρ₁=1 Ωm，ρ₂=10 Ωm，规模均为20 m×20 m，顶端埋深均为20 m，两异常体相距200 m.测量装置为中心回线装置，方形发射回线边长50 m，接收线圈等效面积50 m²，供电电流1 A.

图 9

Fig. 9

图 9 模型二示意图(中心回线装置) Fig. 9 Schematic diagram of model 2(central loop unit)

图 10给出了测线方向28个测深点各采样时间道上的瞬变电磁场剖面响应曲线.平剖面测量中心点位于坐标原点，亦即量异常体中心位置.可以看到，在每个异常体上方表现出了异常特征，与板状体模型类似，早期到为单峰异常，之后变为双峰异常，最后又变为单峰异常吗，且在中期时间段，右侧异常体异常受到左侧较低组异常体的影响，单峰异常并未表现出对称形态.此外可以看出，电阻率相对较小的异常体的异常较大.这说明电磁场在低阻介质中扩散速度较慢；目标体(异常体)电阻率越低，瞬变电磁法的探测效果越好.

图 10

Fig. 10

图 10 模型二瞬变电磁测深多测道剖面图 Fig. 10 Multi-chanel profile for TEM sounding of model two

在瞬变电磁法二维有限差分正演模拟的基础上，展开了瞬变电磁法2.5维有限差分正演模拟研究.其中涉及到的关键问题一是波数域中的初始条件，二是波数域中边界条件的处理.

2.5D正演模拟中采用的是矩形回线源，只能推导出其在均匀半空间上形成的波数域、频率域瞬变电磁场表达式，之后利用正余弦变换、GS变换等数值计算方法得到波数域的瞬变场作为初始条件.

2.5D正演模拟中，地空边界条件扔采用向上严拓的方法.文中推导出了波数域的向上延拓公式.在计算中涉及到虚宗量贝塞尔函数的求解、无穷积分的处理、高斯勒让德积分、高斯拉盖尔积分、三次样条插值等，计算复杂、速度慢.

利用Fortran语言编写了瞬变电磁法2.5维有限差分正演程序，将FDTD数值模拟结果与均匀半空间模型、H型层状模型、K型层状模型的解析解进行了对比，曲线形态与解析解基本吻合，平均相对误差控制在5%以内，证明了算法的正确性与可靠性；将FDTD数值模拟结果与FEM数值模拟结果进行对比，在早期，FDTD比FEM计算结果稍差，而在晚期FDTD计算结果优于FEM计算结果；通过对两个简单模型的模拟，证明了算法的有效性.

致 谢 诚挚感谢强建科副教授、任政勇副教授在本文撰写中给予的亲切指导及宝贵意见.