瞬变电磁测深法(TEM)近年来在国内外得到广泛应用^[1~5]，在国民经济建设中发挥了重要的作用，成为国内外引人瞩目的地球物理勘探方法之一.瞬变电磁法正演计算方法相对比较成熟^[6~9].该方法的反演计算虽然经过大量的研究^[10~15]，并解决了一些地质问题，但是，随着精细勘探任务需求的增加，难以满足实际工程技术的需要.必需研究新的解释方法.近年来，瞬变电磁拟地震成像技术成了一个前沿的研究方向.

可以通过一个经验公式^[16]，把TEM 数据等效转换成平面波场数据^[17]，借用MT 数据的拟地震思路，求取反射系数序列进行成像，但是成像精度相对较低.Lee等^[18~20]根据无耗媒质中电磁场所满足的无阻尼波动方程与高导电媒质中电磁场所满足的扩散方程之间的相似性，推导出了以上两种媒质中的电磁场量所满足的变换式，称为波场变换式.随后，瞬变磁法地震解释方法研究不断得到推进^[21~23].李貅^[24]于2005 年提出波场转换的优化算法，薛国强^[25]于2007年运用这种优化算法成功地转换出虚拟地震子波.

虽然变换出的虚拟波场在形式上与地震波场一样都满足波动方程，但是由于两种波场的物理背景不同，它们之间存在一些重要区别，前者是与感应的瞬变电磁衰减曲线相对应的“反射"子波，是虚拟的，而后者是在弹性介质中传播的客观存在的地震子波.另外，虚拟波场在每一介质中的传播速度不仅与本介质的电导率有关，而且还受相邻介质电导率的影响，这与地震波场是不同的.并且存在严重的波形展宽效应，使得计算得到的虚拟波场数值分辨率降低，成为瞬变电磁法拟地震成像解释的主要障碍之一.显然波形展宽，这一现象严重影响着TEM 成像的空间分辨能力.所以，有必要讨论这种虚拟的波动场波形展宽的原因及对策.

本文首先给出了波形展宽的现象，然后讨论了虚拟波形展宽的原因，认为不是由于波在介质中传播时能量损耗所致，而是波场变换式中高斯分布的核函数随虚拟时间的增加分布范围增大的分布特点所决定的.最后对虚拟波场的离散数据求取反褶积，消除了波场变换的波形展宽效应.通过对理论模型计算处理，证明了该方法可以压缩子波的宽度，并增强瞬变电磁法识别地下电性分界面的能力.

Lee^[18]于1989 年建立了满足时域扩散方程的电磁场与虚拟波场的关系式，即波场到时域场的波场变换式：

(1)

式中，自变量τ 是时间平方根的量纲，函数U(r，τ)是以波速传播的虚拟波动场.H_m(r，τ)为时域瞬变响应扩散场.

这一变换过程称为正问题.如果反过来，已知时域场求波场，则称为反问题.为此可以利用文献[24, 25]中的方法，采用两次最优化技术和正则化算法，解决(1)式的反变换问题，实现了从时域场到波场的变换.

以三层地电断面(K 型、A 型、Q 型)为例进行波场变换，K 型断面地电参数为ρ₁=5Ωm，h1=60m，ρ₂=50Ωm，h2=60m，ρ₃=5Ωm.A 型断面地电参数为ρ₁=5Ωm，h1=60 m，ρ₂=50Ωm，h2=60m，ρ₃=500Ωm.Q 型断面地电参数为ρ₁=25Ωm，h₁=60m，ρ₂=5Ωm，h₂=60 m，ρ₃=1Ωm.图 1为以上三个模型的波场变换结果.对三个模型的转换结果进行分析，其共同的特点就是波形较宽，正负脉冲之间的时间间隔较大，导致该方法在拟地震成像处理中对电性界面的分辨能力下降.

图 1

Fig. 1

图 1 地电模型的波场转换结果 (a)K-型曲线波场转换结果；（b) Q-型曲线波场转换结果；（c）A-型曲线波场转换结果. Fig. 1 Wave-field transformed result of geo-electric models (a) Result of K-type wave-field transformation； (b)Result of Q-type wave-field transformation； (c) Result of A-l:ype wave-field transformation.

在瞬变电磁场转化为虚拟波场的过程中，所得虚拟波场波形会随着时间的增大而不可避免的展宽，严重影响成像的分辨率.原因在于波场变换式(1)中核函数的分布特点，如下的核函数

(2)

为高斯分布，以真实的物理时间为方差，因此随着时间的增加，时域场结束的时间愈加滞后，波场波形也就愈加宽阔.为了说明核函数的分布宽度随时间增大而增大，按照(2)式，分别计算t=80μs、325μs、800μs、2.1ms、8.7 ms、27 ms、81 ms时核函数值随虚拟时间τ 变化时的分布形态.

由图 2可见，t=80μs时，核函数对应的τ值的最大取值范围为0~0.05$\sqrt{s}$；t=325μs时，核函数对应的τ值的最大取值范围为0~0.1$\sqrt{s}$；t=800μs时，核函数对应的τ值的最大取值范围为0~0.15$\sqrt{s}$；t=2.1 ms时，核函数对应的τ 值的最大取值范围为0~0.25$\sqrt{s}$；t=8.7ms时，核函数对应的τ 值的最大取值范围为0~0.5$\sqrt{s}$；t=27ms时，核函数对应的τ 值的最大取值范围为0~0.85$\sqrt{s}$；t=81 ms时，核函数对应的τ值的最大取值范围为0~1.5$\sqrt{s}$.随着时间t的推移，核函数a(t，τ)随虚拟时间τ 的增加分布范围明显增大，这就导致经反变换获得的虚拟波场出现了随着时间的推移而展宽的现象.

另外，在自由空间中，电磁波的波速为c₀ =，而导电媒质中的虚拟波速为两者相差约6个数量级.由于导电媒质中的虚拟波速非常慢，这样，传播同样的距离差，就需要更多的时间.所以，其波形的展宽效应就变得很显著了.

尽管虚拟波形展宽这一现象是由数学变换产生的，但是与地震勘探中的“大地滤波作用"非常相似.在地震勘探中，由于大地滤波的作用，加上噪声和多次波的叠加，使本来可以清晰地呈现出的地层层序的尖脉冲序列“模糊化"，降低了地层的纵向分辨率.可以通过反褶积运算来抵消大地滤波的作用.反褶积作用的结果就是压缩子波，可以通过求取一反滤波因子，把先前得到展宽的波场能量集中化，达到精细分层的效果，这样分辨率可望大幅提高.

Peacock和Treitel^[26]在1969年提出预测反褶积，此项技术在地震资料处理中得到了非常广泛的应用.这项技术同样可以用来处理转换出的虚拟地震子波.

假设，y(t)为实际转换出的拟地震子波，y(t)=(y(0)，y(1)，…，y(n))，我们设计一个滤波因子h(t)，h(t)=(h(-m0)，h(-m0+1)，h(-m0+2)，…，h(-m0+m)，h(t)的起始时刻为-m0，其延续长度为(m+1)，对其进行如下计算：

(3)

图 2

Fig. 2

图 2 核函数值分布性态 (a)t=80μs;(b)t=325μs;(c)t=800μs;(d)t=2.1ms;(e)t=8.7ms;(f)t=27ms;(g)t=81ms. Fig. 2 Distribution property of kernel function

然后，得到延迟时刻a后的子波$\hat{y}$(t+a)，并且使$\hat{y}$(t+a)与(t+a)时刻的虚拟地震子波值误差Q最小，即：

(4)

把(3)式代入(4)式，得：

(5)

要使Q为最小，数学上就是求Q的极值问题.将Q对滤波因子求偏导、化简，并令其为0，即

(6)

(7)

是虚拟子波的自相关，而是子波与期望输出的互相关，所以方程(7)可写成：

(8)

将(8)写成矩阵形式为

(9)

通过解方程(9)，可以求出反褶积因子h(t)，求得h(t)后，就可以通过(3)式，求出转换虚拟波所对应的反褶积处理计算结果.

反滤波算子的长度原则上可以任意选择.地震中反滤波算子的长度一般选择80、120、160 ms、200、240ms等，但我们在计算时要将它们换算为采样点数.要根据插值后数据的总个数来进行选择，当达到一定长度时，反褶积效果趋于稳定.在本文的计算中，反滤波算子的长度基本选取在400~600之间.

图 3为运用了反褶积技术的结果，图 3a是反褶积处理前波形，在反射界面处虚拟反射波波形明显较宽，分辨力下降；图 3b反褶积处理后波形，出现两组虚拟反射波，波形相对较窄，由图可见，通过反褶积处理，可以将电性分界面处产生的被拉长的虚拟反射波锐化，使用瞬变电磁法的分辨能力有较大的提高.

图 3

Fig. 3

图 3 反褶积计算前后的H-型模型子波对比 (a)压缩前子波波形;(b)压缩后子波波形(ρ1=ρ3=100Ωm,ρ2=5Ωm,h1=h2=100m). Fig. 3 Comparison of H-type model wave before and after de-convolution processing (a) Waveform before de-convolutionprocessing; (b) Waveform afterde-convolution processing (H-typemodel: ρ1 =ρ3 =100Ωm,ρ2=5Ωm,h1=h2=100m).

为了验证反褶积的处理效果，首先考虑一维地电模型，将对图 3 中的地电模型(H 型模型：ρ₁ =ρ₃=100Ωm，ρ₂=5Ωm，h1=h2=100 m)的反褶积处理后的虚拟波场进行拟地震成像，图 4给出了H型模型的成像的结果，理论上，地电模型的两个电性界面分别位于地下100m 和200m 处.由图 4可知，浅层电性界面所产生的波场能量集中于100m 左右，而较深界面能量高峰产生于200m 左右.说明在界面深度定位方面与模型预先给定情况大体一致，效果令人满意.

图 4

Fig. 4

图 4 H型模型拟地震成像图 (a) H型模型成像图；（b)H型模型成像切片图. Fig. 4 Result of H-type model pseudo-seismic (a) II-type model imaging map；(b) Section of II-type model.

为了验证反褶积的处理效果，我们还设计了三维地电模型(图 5a，取均匀半空间的电阻率为ρ₁=10Ωm.异常体为30m×30m×50m 的块体，电阻率为ρ₂=300Ωm，埋深为70m.在平面上，选取与模型水平延伸方向为测线方向，测线长度为100m，(图 5b)，采用均匀半空间下三维积分方程法进行模拟计算，按文中提出的方法对正演计算结果进行虚拟子波提取并进行反褶积处理.图 5c是反褶积前的结果，图 5d为反褶积处理后的结果.从图中可以看到，经反滤波处理后，波形得到锐化，图 5e为拟地震成像剖面.从该图可以看出，上部波形能量集中在70m 处，下部波形能量集中在120 m 处，与模型设置吻合很好.特别是下部界面，是从高阻到低阻，本身是弱反射面，波形展宽很严重，经反褶积后，波形得到了锐化.

图 5

Fig. 5

图 5 反褶积计算前后三维模型子波对比 (a)三维地电模型；（b)计算场点平面分布；(c)压缩前子波波形；（d）压缩后子波波形；（e）拟地震成像剖面. Fig. 5 Comparison of 3-D body wave before and after de-convolution processing (a) 3-D geo-electric model； (b) Plane distribution of calculated field point；(c) Waveform before de-convolution processing； (d) Waveform after de-convolution processing；(e)Imaging section by the de-conversion wave data.

由于导电媒质中电磁场所满足的扩散方程与无耗媒质电磁场所满足的无阻尼波动之间存在数学转换关系式，从而可以实现瞬变电磁法拟地震资料解释.这一变换的实质是提取瞬变电磁响应中与传播有关的特征，压制或去除电磁波传播过程中与频散和衰减有关的特征.但是转换出的虚拟波往往波形较宽，使得分辨能力下降.

波形展宽的主要原因是由于波场变换式中高斯分布的核函数以真实的物理时间为方差，随着时间的推移，分布范围增大，虚拟时间愈加滞后，波场波形也就愈加宽阔.

通过对一维地电模型和三维模型的虚拟波场信号进行反褶积计算，所得到的结论为：可以利用反褶积技术把宽幅的虚拟子波汇聚到相对较窄的区域内，消除波场变换的波形展宽效应.证明瞬变电磁拟地震方法可以增强瞬变电磁法识别地下电性分界面的能力.研究成果对进一步推广瞬变电磁法的理论与应用具有重要意义.