1 引言

瞬变电磁拟地震成像是实现瞬变电磁三维反演的有效手段，其中关键问题之一是实现由瞬变电磁扩散场到虚拟波场的转换.只有进行了波场逆变换，才有可能实现基于波动方程的电磁法拟地震解释.

所谓的电磁法拟地震解释主要有三种，第一种是根据电磁波在导电介质中和地震波在弹性介质中传播的相似性以及反射函数表达式之间的一致性把地震勘探中的一些基本方法应用于电磁测深资料解释.主要成果有：Levy（1988）根据弹性波与大地电磁场之间的类比性，获得了与反射地震类似的拟反射函数的脉冲响应时间断面图，将大地电磁资料用线性规划技术对一维和简单的二维地电断面进行了成像计算，并推导了直流问题的拟脉冲响应函数^[1]；王家映根据电磁波在导电介质和地震波在弹性介质中传播的相似性及其反射函数表达式之间的一致性，在国内率先提出了对大地电磁资料进行拟地震解释的方法和思路^[2-5]；郭文波、薛国强等基于中心回线装置观测成果与MT成果的一致性，借助大地电磁拟地震解释，进行了TEM拟地震成像解释^[6-7].第二种是借助偏移成像的广义概念在扩散场条件下实现电磁偏移.主要成果有：前苏联学者Kjahos吸取了“偏移成像”的广义概念，在电磁法中确立了正则偏移和解析延拓偏移两种方法^[8-10]；Zhdanov等^[11-13]借鉴了地震勘探中的逆时偏移概念，对时间域的瞬变电磁场和频率域的大地电磁场进行了逆时偏移成像的深入研究，提出了偏移电磁场的概念，并且建立在逆时偏移电磁场基础上对二、三维反演问题也开展系列研究.第三种方法是基于Lavrent'ev（1980）的研究即瞬变电磁扩散场与虚拟波场之间存在着准确的数学变换关系，实现扩散场到波场的变换，在虚拟波场的条件下引入地震勘探处理和解释的一些基本理论进行电磁成像^[14].主要成果有陈本池、李金铭等利用奇异值分解实现瞬变电磁波场变换，利用有限差分进行二维瞬变电磁拟地震成像解释^[15-16]；李貅等对瞬变电磁时间信号进行了分段处理，结合正则化算法分段实现了瞬变电磁波场变换，利用Kirchhoff积分实现了瞬变电磁三维曲面延拓成像^[17-21]，并基于此开展了提高虚拟波场分辨率的系列研究^[22-24].但是，目前这方面的研究在算法和方法适用性方面仍存在不少问题.在对复杂介质进行大地电磁拟地震解释时遇到一些困难，如电磁波速度求取问题；基于扩散方程的逆时偏移成果仅限于单个或者多个孤立异常体的成像；利用扩散场到波场的变换关系可以实现复杂介质的成像，但是由于波场反变换方程是第一类Fredholm积分，属于典型的不适定问题^[25-27]，李貅等早期研究为了降低问题的病态程度，采取了时间分段手段，有效降低了矩阵的病态程度，但这也引入了段与段的衔接问题^{[18, 20, 28]}.

本文在李貅教授原有的研究基础上进行瞬变电磁的波场变换算法，采用预条件正则化共轭梯度法实现全时域的波场变换算法，既有效避免了分段，又实现了由扩散场到波动场的变换.利用虚拟波场可以在波动方程下进行虚拟波场偏移成像.

2 波场变换的基本原理

早在1972年，Kunetz^[29]在对大地电磁解释和反演问题的研究中发现，电磁场满足的扩散方程与波动方程之间具有相似关系.随后，Lavrent'ev于1980年，在《Ill-posed problems of mathematical physics and analysis》一书中对扩散方程与波动方程之间的对应关系进行了深入研究，并给出了两个方程之间准确的数学变换关系式^[14].Lee在1987年与1989年的研究成果中又进一步证明了这种数学变换适用于任意的场分量^[30-31].

在导电介质中，忽略位移电流，瞬变电磁场满足扩散方程.为了不失一般性，取f（x，y，z，t）为瞬变电磁场的电场或磁场分量函数，其满足的扩散方程如式（1）所示：

(1)

其中μ为介质磁导率，σ为电导率.引入函数U（x，y，z，τ），其满足波动方程如式（2）所示：

(2)

根据文献[14, 16, 21, 29]可知，由扩散方程到波动方程转换的对应关系表达式为

(3)

（3）式为第一类Fredholm型积分方程，由扩散场求波场是典型的不适定问题.将其离散后得到的线性代数方程组是病态的，且随着阶数的增加，矩阵条件数急剧增大，病态性更加严重，因此必须选用可靠的离散方式和稳定的数值方法.本文采用预条件正则化共轭梯度法实现波场反变换的计算^[32-33].

3 波场变换的数值离散

将式（3）进行离散，写成数值积分形式

(4)

其中a（t_i，τ_j）=τ_je^{-τ _j2 /4t_i}，w_j为积分系数.核函数a（t_i，τ_j）是随着τ的增加达到某一极值后快速衰减到零的函数（如图 1所示），对于不同的时刻t_i，核函数曲线的形态基本一致，主要差别在峰值大小和展布范围.了解积分核函数的这些特征，有利于开展后面的波场正变换和波场反变换的研究.

计算式（4）的关键是如何选取一组τ_j（j=1，2，…n）以及对应的一组w_j（j=1，2，…n）使其最好的满足式（4）.对于虚拟时间与积分系数的选择，文献[18, 20]通过给定的特殊积分，根据两步最优化选定一组最优的虚拟时间和一组最优的积分系数.但是文中矩阵A的条件数较大，矩阵方程严重病态.虽然如此，我们仍可以借助特殊积分选取合适的离散方式，提取最优的虚拟时间τ和最优的积分系数w.

我们对比了一些常规的离散方式，如等间距离离散，对数等间距离离散，等面积离散和高度等间距离离散四种方式，得到不同离散方式的虚拟时间分布示意图（图 2）.根据不同离散方式选取虚拟时间，采用数值积分方法进行积分，比较各数值积分精度和离散后各矩阵条件数，结果如表 1所示.通过比较选取均方误差最小，形成系数矩阵条件数最好的等间距离散.

4 预条件正则化共轭梯度法实现波场反变换

考虑到文中方程组系数矩阵较大，单纯的正则化共轭梯度法（RCG）或者预条件共轭梯度法（PCG），都不能很好地解决问题.因此，将两种方法相结合形成预条件正则化共轭梯度法（PRCG）.由于在实际应用当中主要采集磁场的垂直分量，因此文中以磁场垂直分量为例进行分析.

将式（4）写成矩阵形式为

(5)

其中A=[a_ijw_j]_m×n，系数矩阵A包含积分系数w_j，U=[u_j]_n×1为虚拟子波，F=[h_i]_m×1为接收的瞬变场时间信号.

为了利用共轭梯度迭代，将式（5）转化为

(6)

只要A是列满秩矩阵，A^TA就是对称正定矩阵，因此可以利用共轭梯度法.但是A^TA的条件数较A的条件数更大，使方程的病态更加严重.

为了进一步降低矩阵的条件数，改善方程的病态程度，在进行正则化共轭梯度之前，对系数矩阵进行预条件.对于预条件矩阵的构造采用超松弛预条件法，这是因为对称超松弛预条件是一种较为有效的预条件方法，不仅预条件子容易求得，而且有效降低矩阵的条件数^[34-37].数学上已经证明经过超松弛预条件后，矩阵条件数降为原来的平方根^[36].

设矩阵

(1)

可分解为

(2)

其中

(3)

(4)

K，C_l和C_u分别为S的对角元、下三角元和上三角元，ω为（0，2）内的参数.于是我们可以选择预条件矩阵P如式（7）所示：

(7)

假设根据矩阵A（v）构造的预条件子为M（v），构造新的方程如式（8），矩阵M（v）^-1A（v）接近单位阵，因此迭代很快收敛.具体计算过程如图 3所示，其中k_max为外层循环最大迭代次数，ε为正则化共轭梯度法迭代终止条件，l_max为内层循环最大的迭代次数，ξ为内层共轭梯度迭代终止条件，v为选定的正则化参数.

(8)

其中，x^k为第k次迭代的x值；x的初值x^（0）选为单位向量；A（v）=vI+A^TA.

正则化参数v的选择非常重要，正则化参数v（δ）使U在近似性与稳定性之间进行优化选择. Zhdanov等^[38-39]提出正则化因子不断递减的自适应算法.并与L曲线法做了比较，认为采用自适应算法所得反演结果不逊于L曲线法所得结果，由于L曲线法需要通过多次反演计算来确定最佳的正则化因子，计算量将成倍增长，而自适应算法只需在每次迭代前确定参与当次反演的正则化因子，因此可以大大减少计算时间.在文献[40]中，王彦飞依据偏差原理提出了重开始共轭梯度法（RSCG）计算最佳正则化因子，并与CGNR和CGER方法进行比较，认为RSCG方法更加稳定^[40].正则化方法在大的电磁反演中已有广泛的应用，如文献[41-44]等均利用正则化方法实现了大地电磁反演，并给出了相应的正则化因子的确定方法.但是，瞬变电磁场较大地电磁场时间范围与动态范围均更广，不适定性更加严重，不能简单地引用大地电磁反演中使用的正则化方法.

经分析，文献[40]中所述的方法更适合本文情况.因此，文中依照文献[40]重开始共轭梯度（RSCG）计算方法，根据部分先验信息来选择最优的正则化参数值.正则化因子v为一渐变的量，其初始值v₀为数据拟合泛函与稳定泛函的比值，在前期大量模拟计算的基础上基本可以确定正则化参数的初值，根据经验可以取v₀等于0.00005为分段最大的正则化因子.在此后的迭代过程中，如果数据拟合残差随迭代次数逐渐变小，正则化因子可保持不变，否则按照（9）式进行选择：

(9)

上式中ξ为经验系数，有ξ＞1，k表示预条件正则化共轭梯度（PRCG）的迭代过程中第k次迭代.

在完成积分离散以及正则化参数选择后，对方程（4）式运用预条件正则化共轭梯度法进行求解计算.取波场理论值为U（x，y，z，τ）=1，这时方程为一简单积分，求得函数值表 2所示为适用时间区间[0.000078，0.006280]的正则化参数.利用（PRCG）法求得反变换结果如图 4a所示，最大误差小于2%，可知文中所述方法满足要求.

为了便于比较, 给出采用正则化共轭梯度法(RCG)获得的反变换结果, 如图 4b所示.比较图 4a与图 4b可知，采用预条件处理后，正则化参数可以选取较小（v=0.0011285）就能达到很好的收敛效果；不采用预条件处理，正则化参数选择较大时（v=0.008），其收敛效果尚不如采用预条件技术的结果.由此可知预条件处理是必要的.

众所周知，正则化参数需要在稳定性与分辨率之间平衡，正则化参数越大，分辨率越低.为了进一步说明问题，现给出同一模型两种方法的波场变换结果，且两方法选择的正则化参数一致，为v=0.003，如图 5所示，图 5a为采用PRCG方法获得的波场变换结果，图 5b为采用RCG方法获得波场变换结果.由图可知模型具有两个电性分界面，且两个界面之间为一薄层，采用预条件正则化共轭梯度方法的结果可以分辨出薄层以及两个电性界面，而正则化共轭梯度方法获得的波场结果不能区别出两个界面.由此可以说明，文中所述的PRCG反变换方法比RCG方法的分辨率高.

文献[45]中，张军、李貅（2009）等对不同时刻进行了分段，为了便于比较，截取相同时段计算结果与其进行对比，如图 6所示.对比可知，在相同时段内，本文方法计算精度高于分段正则化算法精度.

5 瞬变电磁虚拟波场分析 5.1 不同地电模型波场变换响应

上文根据特殊子波对方法进行了验证，说明方法的正确性.为了验证方法对复杂介质的适用性，利用文中波场变换方法分别对两层模型和三层模型响应进行波场变换分析.瞬变电磁响应反映地下介质的电阻率分布情况，当介质为高电阻率时信号衰减较快，当介质为低电阻率时信号衰减较慢.根据瞬变电磁衰减信号的特征，我们可以判断层状介质的电阻率相对变化情况.以高斯脉冲为子波给出虚拟波场值，通过式（3）可以获得与之对应的瞬变电磁衰减信号，并判断地下介质电阻率分布情况.根据瞬变电磁衰减信号可知，含有一个子波的为两层模型，含有两个子波的为三层模型，且子波为正表示下层电阻率较上一层电阻率低，子波为负表示下层电阻率较上一层电阻率高.

根据特殊子波U（x，y，z，τ）=1可以确定正则化参数的合理范围.选择正则化参数后，对给定的地电模型进行模拟.对两层模型的反变换结果如图 7、8所示，由图可知，反变换所得波场与已知虚拟波场基本吻合.

三层模型的波场反变换结果如图 9、10、11、12所示，由图可知，反变换结果与已知波场基本吻合，说明方法适合于复杂介质.图中第二波峰较第一个波峰振幅有很大衰减，子波宽度略有增加.这与波在有耗介质中的传播规律相符，从而证明了算法的正确性.另一方面，模型的反演结果充分说明了瞬变电磁方法对浅部异常分辨率较高.随着深度的增加，电磁波高频部分损失严重，响应主要为低频响应，因此对深部异常的分辨率降低.

5.2 含噪信号的波场变换分析

野外测量数据受各种因素影响难免会有干扰，因此有必要研究含噪信号的波场变换.文中以单脉冲响应的波场变换为例进行分析，加性噪声分别为1%、2%、3%、5%.如图 13、14、15、16所示为不同信噪比的D型模型时域响应信号和与之对应的波场变换结果；图 17、18、19、20为不同信噪比的G型地电模型时间域响应和与之对应的波场变换结果.

由图 16b与20b可以明显看出噪声使得时间域响应曲线变得很不平滑.由波场反变换结果来看，当信噪比为小于5%时，反演值与期望值吻合得较好，当信噪比超过5%时，反演值波动剧烈，与预期结果差距较大.因此对于野外信号采集必须控制噪音为5%以下，并且对于局部不光滑数据不能直接进行处理，需要对数据进行平滑去噪后才能进行波场变换.因此本文方法完全适用于实际数据的处理.当然为提高波场变换的适用性也可以采取其他手段，如合成孔径方法来提高波场的抗噪能力.

5.3 实测数据分析

在前文中已经利用模型对方法进行了验证，现对实测数据进行处理，进一步对算法进行验证.以某煤田采空区实测数据为例进行分析.仪器选用GDP-32电法工作站，用中心回线方式进行测量，测线为1320、1360、1400，测点300-900.

测区视电阻率剖面如图 21所示.由图可知，区域视电阻率变化较大，变化范围从几十到几百欧姆米，总体表现为地表黄土覆盖视电阻率高，基底视电阻率较高.在点号300-500深度120~350 m左右范围内，有一低阻异常，结合工区资料，推断为富水采空区.根据视电阻率断面图可大概推断采空区分布在点号（300，500）区间，由于深层电阻率变化缓慢，不容易估计采空区底界面，从视电阻率断面图可初步估计底界面在300~400m范围内.

为了验证方法效果，给出测线1320、1360、1400，点号为300到500间11个测点数据的全域波场变换结果.如图 22所示，其横轴为点数，纵坐标为虚拟时间，单位为，在虚拟时间为100处有一明显的电性分界面，在虚拟时间为400左右也有一电性分界面分布.根据图 21可知地下介质类似于三层模型，与图 22波场变换结果相符.但是波场变换结果显示的界面与视电阻率界面并不一致，这是由于图 22纵轴为虚拟时间，要想得到真实的地下介质起伏形态还需要进行速度分析与延拓成像，虚拟波场速度是一个与地下介质电阻率有关的量，经过波场延拓后，得到深度剖面，其显示结果应与视电阻率起伏形态相似.综上所述，本文方法正确有效，能够用于实际数据处理，并且对推断多层采空具有明显的优势.

6 结论

文中主要研究由扩散方程到波动方程转换的数值计算方法.首先，对前人推导的波场表达式进行数值离散，分析了多种离散方式的优劣，选用效果较好的等间距剖分进行离散.然后，在前人工作基础上，采用预条件正则化共轭梯度法，实现了全时域波场反变换.正则化共轭梯度法适合病态方程的求解，采用SSOR预条件子进行预条件，将系数矩阵条件数降为原条件数的平方根，有效改善方程的病态性，这样选择较小的正则化参数就可以得到较高的精度，有利于运用正则化共轭梯度实现方程的迭代计算.

实现波场反变换后，对算法进行了验证，并对瞬变电磁虚拟波场的性质进行了详细讨论.利用均匀模型（即子波为常数时），对反变换计算方法进行验证，得到的子波与理论子波误差小于2%，说明PRCG算法是有效的.对两层（子波为单脉冲）和三层（子波为双脉冲）模型进行反变换，得到的子波与理论子波相似，由此可以证明该方法适合于复杂介质.对两层模型和三层模型进行了详细的分析，可知虚拟波场对浅层目标分辨较好，随着深度的增加，虚拟波场幅值迅速衰减，子波展宽严重，分辨率逐渐降低.

利用文中所述方法和选择的正则化参数，在选定时段对于含有噪声的瞬变电磁信号进行波场变换，当干扰在5%以下时波场变换结果与理论值吻合较好.但是当干扰大于5%时，波场变换结果与理论值相差较大.因此瞬变电磁信号的前处理，即信号平滑与去噪是必要的.当然我们也可以增大正则化参数，但这样做势必影响波场的分辨率.

经过波场反变换后，虚拟波场不仅满足波动方程，而且与地震子波类似，具有波的传播特性.这为后续利用波动方程偏移成像方法提供了良好的基础.