航空瞬变电磁法是航空物探中近几年发展起来的方法，具有速度快、成本低、通行性好、覆盖面积大、可用于山区和海域勘探等优势^[1-3].目前航空瞬变电磁法解释仍在一维水平^[4-8]，主要用于地质填图、确定地质构造等，将其用于详细地质勘查还有一定的距离，不能满足矿产资源精细探测的需要.这就要求提出新的解释技术提高航空电磁法的解释水平，使其达到甚至超过地面勘探精度，从而代替传统的地面瞬变电磁法勘察，以满足高原或海洋区域先进的探测解释技术.

随着瞬变电磁波场变换算法的提出^[9-10]，实现了由具有扩散特征的瞬变场向虚拟波场的转变，为实现航空瞬变电磁法的合成孔径成像创造了条件.瞬变电磁合成孔径成像技术是借用合成孔径雷达成像的思路^[11-16]，即利用机载真实孔径发射线圈与目标的相对运动，把尺寸较小的真实天线孔径用数据处理的方法合成一较大的等效孔径的发射线圈，使其分辨能力更高、穿透能力更强.对于航空瞬变电磁法而言其观测方式与机载合成孔径雷达十分相似，完全可以借助合成孔径雷达的成像思想，实现虚拟波场条件下的航空瞬变电磁合成孔径成像.

虽然合成孔径技术在雷达成像中早已广泛应用，但与合成孔径雷达信号不同，瞬变电磁场满足的是扩散方程，场的变化特征是随着时间衰减的，显然不利于相关叠加处理，也不利于合成成像.为此，必须利用瞬变电磁场与虚拟波场间存在的数学关系式^[10]进行波场变换，通过优化算法提取出虚拟波场信号^[17-19].数字模拟和模型试验的结果都证明相邻位置上同一地质体的反射回波具有较好的相关性^[20]，因此根据不同位置信号的相关系数生成不同的权值函数，相邻各列信号在做相关叠加时以权函数进行加权，加强重建的地质异常体信号，从而提高信噪比，达到突出弱异常的目的，进而提高分辨率，加大勘探深度.在分析合成孔径雷达算法的基础上结合瞬变电磁信号的特点，对采样信号进行相关加权叠加形成瞬变电磁合成孔径数据体，对该数据体进行克希霍夫偏移成像，得到地质体高清晰度的数字图像.

通过模型数据的计算，并进行合成孔径成像处理，说明了该方法的有效性，对以往地面数据的再处理，并与原来的处理结果相比较，说明了该方法在提高瞬变电磁分辨力方面的优越性.该方法的成功应用，将对用高分辨的航空瞬变电磁法代替海面和山区地面勘探具有重要意义.

瞬变电磁场满足如下扩散方程：

(1)

接收线圈中接收到的电磁响应

(2)

的变化特征是随着时间单调衰减的，显然不利于相邻点信号的相关叠加计算，更无法进行偏移成像.较好的解决办法是将瞬变场变换成虚拟波场，利用来自同一反射点的信号，对有效波的相关特性，实现合成孔径的相关叠加计算，提高分辨率，加大勘探深度，最终利用克希霍夫积分公式进行偏移成像计算.

瞬变电磁场与虚拟波场之间存在如下对应关系^[10]：

(3)

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

(3)式是典型的第一类Fredholm 型算子方程，已知瞬变场求波场显然具有不适定性^[21].采用最优化技术和正则化算法求解其反变换问题，可以实现从瞬变场到波场的变换^[17].

经波场变换获得的虚拟子波存在严重的波形展宽现象，使得计算得到的虚拟波场数值分辨率降低，严重影响着TEM 成像的空间分辨能力.为此，对虚拟波场的离散数据求取反褶积，消除波场变换的波形展宽效应^[22].

设U(r，τ)为实际转换出的虚拟子波，希望通过一个反褶积滤波因子h(t)，进行如下计算，求取一个宽度得到压缩后的新的子波：

(4)

其中，反褶积滤波因子h(t)可用最小平方反褶积求得^[22].

瞬变电磁场经过波场变换，已经把原来的感应场转换成了波场，每一点的数据相当于自激自收的波动场.而前人已经通过实验分析确定瞬变电磁场在多激励源情况下存在场的相关叠加性^[20]，基于瞬变电磁场的上述特点，采用相关叠加的方法来进行合成孔径.其合成示意图如图 1所示.

图 1

Fig. 1

图 1 合成孔径示意图 Fig. 1 Schematic diagram of synthetic aperture

首先选取一个中心点，取为第i点，此点的波场值可表示为U(r_i，τ)，其中r_i为i点到-N，…，N内某点的距离，τ 为相对时移量.然后选定2N+1 个测点的长度为合成孔径的长度，即选取i点左右两侧从-N到N的测点分别与中心点做相关，其归一化的互相关系数

(5)

其中m为每一测点的时间道数.

互相关系数ρ(r_i，τ)表示了两列波场的相关程度，但它与相对时移量τ 值有关，通过改变时移量τ值，找出最大相关系数ρ^max(r_i，τ^m)对应的时移量值，称为最佳延时，记为τ^m.由此可以得到2N+1个最大相关系数ρ^max(rk，τ^m_k)和最佳延时τ_k^m(k=-N，…，-1，0，1，…，N).之后，用相关计算得到的最大相关系数作为权系数分别乘上各点的波场值，叠加到中心点，最终可得到中心点的合成值为

(6)

在一条测线上依次移动，则可得到第i+1，i+2，i+3，…个点为中心点的合成值.

为了说明原理，以上只给出了剖面上的一维合成孔径算法，也可以推广到平面上的二维合成孔径计算，进一步还可以发展自聚焦算法等.

在获得了观测平面上的合成孔径波场值后，借助地震勘探中的克希霍夫偏移成像方法，依据瞬变电磁虚拟波场满足的波动方程，可给出Kirchhoff积分解^[22]：

(7)

由于偏移是一个求取记录的逆过程，所以可设自激自收的上行波为G（x，y，z₀，t），则上式可写为

(8)

这即是波场的向下延拓.

在实际求取时，因为采集的数据都是离散数据，计算时需要对克希霍夫积分进行离散化，采用三维边界元技术对上述积分方程进行求解，通过给出向下的偏移距离，可以得到地下的三维偏移成像结果^[23].

为了验证本文提出方法的有效性，设计了半空间中赋存高阻块状异常体和低阻块状异常体的三维模型，采用中心回线激发，对所得的电磁响应进行波场变换，经合成孔径计算得到合成波场值，最终用三维边界元法实现克希霍夫偏移成像计算.

其模型三维示意图如图 2a所示，异常体为30m×30m×50m 的块体，发射边长为100m.地表测点分布图如图 2b所示，图中所示主剖面上共布置11 个测点，点距为10m.

图 2

Fig. 2

图 2 模型三维示意图 (a)三维模型立体图；（b)主剖面平面图. Fig. 2 Schematic of three-dimensional model(a) 3D model Stereogram； (b) Main section plan.

取均匀半空间的电阻率为ρ₁=10Ωm，高阻块状异常体电阻率为ρ₂=300Ωm，顶板埋深为70m.首先，用正演方法获得了主剖面上的视电阻率断面图，如图 3所示，可以看出高阻异常的中心与模型的位置并不对应，出现了向上偏移的现象，从电阻率断面图中已经很难确定异常体的位置和大小.

图 3

Fig. 3

图 3 高阻异常体视电阻率剖面 Fig. 3 Apparent resistivity protiles

利用瞬变电磁合成孔径算法，对该模型计算了11个测点的数据.对比图 4(a，b)可以看出，在相关性较好的中心部分合成后的波形得到了一定的增强.在相关性不好的边缘部分，波形得到弱化，合成孔径的有效性得到验证.最后，经过速度分析、波场克希霍夫偏移成像，得到主剖面上的成像结果，如图 5 所示.从图中可以看到，上下两组波峰分别对应在70m 和120m 处，这和模型设置吻合很好，纵向分辨率得到了提高.第二界面波形幅值相对第一界面减小了很多，应该是从高阻到低阻本身是弱反射面所致.和以上电阻率断面图相比，它具有更高的分辨率.

图 4

Fig. 4

图 4 高阻异常体模型计算结果()合成孔径成像前电磁波场变换结果；（b)合成孔径后电磁波场变换结果. Fig. 4 Model of processed results(a) Results before synthetic aperture wave-field transformation；b) Results after synthetic aperture wave-field transformation.

图 5

Fig. 5

图 5 高阻异常体主剖面偏移成像断面图 Fig. 5 Main profile migration imaging section diagram

均匀半空间的电阻率取为25Ωm.低阻异常体电阻率为5Ωm，模型大小尺寸不变，将其顶板埋深加大到120m.主剖面上的视电阻率断面图如图 6所示，可以看出低阻异常的中心与模型的位置也不对应，出现了向下偏移的现象，从电阻率断面图确定异常体的位置和大小也会出现偏差.经波场变换、相关合成等一系列处理，得到合成前后的主剖面对比图，如图 7所示，可以看出，在相关性较好的中心部分合成后的波形得到了增强.在相关性不好的边缘部分，波形得到弱化，同样也证明了合成孔径的有效性.

图 6

Fig. 6

图 6 低阻异常体视电阻率剖面 Fig. 6 Apparent resistivity profiles

图 7

Fig. 7

图 7 低阻异常体模型计算结果 (a)合成孔径成像前波场变换结果；b)合成孔径后波场变换结果. Fig. 7 Model of processed Results (a) Results before synthetic aperture wave-field transformation； (b) Results after synthetic aperture wave-field transformation.

图 8给出了偏移成像结果，可以看出两个界面的存在，但是由于深度加大又是低阻模型体波形出现展宽现象，位置向下略有偏移.

图 8

Fig. 8

图 8 低阻异常体主剖面偏移成像断面图 Fig. 8 Main profile migration imaging section diagram

山西很多地方煤矿在开采过程中，长期受到煤矿积水、陷落柱、采空区等矿井地质灾害的困扰，严重影响着矿井的生产建设.采用瞬变电磁法对煤矿的采空区做系统的研究，探明其范围和埋深.

工作中使用的仪器为GDP-32II电法工作站.考虑到工区的实际地质情况及所需探测的地质体，采用瞬变电磁法进行探测.工作中根据实际试验情况，结合瞬变电磁特点，选择采用大定源回线装置.其发射边框采用边长较大的矩形回线(边长300m)，而接收采用专用的瞬变电磁探头接收(接收面积2000m²).发送频率25 Hz，时间范围：0.087~8ms其中测线间距20m，测点间距20m.

图 9a是视电阻率平面图，由图可知，在测区的西北部、东南部出现两个比较明显的低阻异常区，推测是由于采空区充水引起.参考全区视电阻率断面图成果，以65Ωm 等值线圈定两个低阻区，即A 和B区.图中划红色的正方形框为三维成像数据体采集区域.图 9b为该区域的主剖面视电阻率等值线断面图，红线左侧为三维成像数据体采集区域，从断面图中可以看出，在浅层，电阻率较高，达到100Ωm 以上.从120 m 左右开始，电阻率快速下降至150 m处达到极小，约55Ωm.之后随深度增加，电阻率又开始缓慢提高，至500m 又恢复到100Ωm 左右.从图中可以隐约看出两个电性分界面，第一层位于130~160m 深度附近，第二层界面由于变化缓慢不易确定，初步估计在250~350m 深附近.

图 9

Fig. 9

图 9 工区成果图 (a)视电阻率平面图；（b)典型剖面视电阻率等值线断面图. Fig. 9 Work area results map (a) Apparent resistivity plan；b) Typical profile of apparent resistivity contour map section.

将工区1300~1500号线的300~500点的数据组成一个三维数据体，经波场变换和适当均衡处理后，进行合成孔径计算，得到整个区域的合成孔径波场三维像图，其镂空效果图如图 10 所示，可以看出在经过相关叠加后，突出了有用信号，增强了信噪比，成像效果明显.将合成后的数据进行偏移成像，得到以深度为纵轴的三维成像图，如图11 所示.图中蓝色代表波场负值区，蓝色越深则表示负值越大；红色代表波场正值区；而黄色､绿色为过渡色，表示波动幅值较小的地方.从图中可以看到，和电阻率的断面图相比，此时的分层已经非常明显，通过对成像结果进行解释，认为有两层采空区分布，上层采空区顶板位于地下140m 处；下层采空区顶板位于地下280 m 处，由于充水的影响会使异常产生放大作用，采空区的分布范围也会有一定的放大.经验证解释结果与实际地质情况吻合较好，充分证明了波场变换后合成孔径方法的有效性.

图 10

Fig. 10

图 10 典型地段成像三维镂空图 Fig. 10 Lot imaging 3D diagram of typical hollow

图 11

Fig. 11

图 11 三维波场成像效果图 (a) X方向切片；（b) y方向切片. Fig. 11 3-D Wavefield imaging resutt (a) X direction slice； (b) Y direction slice.

瞬变电磁法在我国已应用与发展了多年，但是仍有很多亟待解决的问题.近几年，瞬变电磁拟地震偏移成像成为一个研究热点.瞬变电磁波场变换算法的提出，实现了由扩散的瞬变场向虚拟的波动场的转变，使得瞬变电磁拟地震成像技术成为可能.本文是在以往波场转换研究的基础上，做了进一步的尝试与研究.

为进一步提高成像分辨率，从现有数据资料中提取更多的有用信息，借鉴合成孔径雷达的思想，提出了瞬变电磁合成孔径算法.为验证所提方法的有效性，特别设计了高、低阻块状异常体模型.从不同深度模型的处理中发现，相关叠加合成确实具有增强有用信号、提高信噪比、提高分辨率的诸多优点.对野外实测资料的处理也表明了相关合成方法的有效性.

将传统的以单点为主的处理方式发展成为以测点为中心的在一定范围内多点相关叠加合成的处理方法，以期达到提高信噪比，突出弱异常，进而提高分辨率、加大勘探深度的目的.合成孔径思想引入到瞬变电磁法中，是一种大胆的尝试，由于该项研究是具有探索性的，很多工作有待进一步深入的研究.

该方法如果进一步完善并加以应用，将会对高分辨率航空瞬变电磁勘探具有重要意义.