多辐射场源地空瞬变电磁法是基于单辐射场源地空瞬变电磁法的一种电磁勘探新方法，该法一方面延续了单辐射场源地空瞬变电磁法工作效率高、观测信号信噪比高、勘探深度大的优点，另一方面又在单辐射场源地空瞬变电磁法的基础上做了进一步的发展，由于采用多个发射源同时发射大功率瞬变电磁场，通过调整发射源的位置及电流方向，该法可有效加强不同分量采集信号强度、削弱随机噪声、减少电性源体积效应的影响、更全面反映地下异常体信息、进一步提高勘探深度(张莹莹等，2015).

单辐射场源地空瞬变电磁法目前已在地热调查、火山结构调查、地下水盐渍化调查、地下水监测等领域得到成功应用(Mogi et al.，1998，2009；Allah et al.，2011；Wang et al.，2013；Ito et al.，2014).根据Mogi等(2009)，目前已知的单辐射场源地空瞬变电磁法的最大探测深度已达800 m，若采用多辐射场源，还可进一步提高该法的探测深度.到目前为止，单辐射场源地空瞬变电磁法的解释方法中最常用到的是视电阻率成像，但是视电阻率的分层效果不够好，导致解释精度仍不够高.李貅等结合TEM三维拟地震解释技术，建立了地空逆合成孔径成像的三维解释理论体系，完善了地空电磁探测体系，但是该算法的计算效率较低，无法进行快速实时成像(李貅等，2015).为了推广该法，丰富地空瞬变电磁法解释体系，研究多辐射源地空瞬变电磁法的快速成像方法是有必要的.

视纵向电导参数对时间的灵敏度高，不仅可以在较早的时间范围内以较高的灵敏度和信噪比将电性层分辨清楚，对深层位电性特征在晚期部分也能反映出来(Sidorov and Tickshaev，1969；Liu and Asten，1993；Tartaras and Zhdanov，1996).近年来，大量学者对视纵向电导解释方法做了大量的研究工作，Tartaras等(2000)基于S反演方法，建立了瞬变电磁数据快速成像算法；Zhdanov等(2002)基于S反演方法提出局部S反演方法，将各向异性考虑在反演中，并将该法应用 在地面和航空瞬变电磁数据中，取得了良好的效果；严良俊等(2003)从“薄板理论”出发，实现了全区视纵向电导与视深度的反演计算，所求取的全区视纵向电导曲线不受过渡期影响，且具有良好的分层性，算法简单易行且有效；李貅等(2003)采用B样条函数实现了电导参数的微分成像，该法能够以“电性同相轴”的形式直观、形象地刻画电性界面的分布 形态；李貅等(2006)将瞬变电磁法引入隧道超前地质预报中，采用视纵向电导参数进行解释，该法对水体不良地质体反映灵敏，取得了良好的应用效果；薛国强和李貅(2008)以二次电导微分参数为特征量建立了隧道超前预报成像系统，以二次电导微分参数绘制成像剖面，该法对隧道掌子面前方水体病害预报效果明显；王文忠和刘继东(2008)将视纵向电导微分成像方法应用于矿井地质灾害调查，结果表明该法比视电阻率曲线的解释结果更加可靠；苏茂鑫等(2010)利用视纵向电导对电性界面识别能力强的特点，将视电阻率等值线图与视纵向电导二次微分成像联合起来，提高了隧道超前地质预报的精度及准确度.但是，上述方法大都是针对中心回线方法，无法解决带偏移距数据的快速成像解释问题，且都是以单点数据为依据进行解释处理，没有较好的利用不同测点数据间的相关性.

本文结合等效导电平面方法和相关叠加合成算法，实现了多辐射场源地空瞬变电磁法的快速成像解释算法.基于等效导电平面方法推导的电偶极子磁场响应公式不仅可用于解决单辐射场源、多辐射场源地空瞬变电磁法的快速解释问题，利用电偶极子叠加原理还可实现带偏移距的地面瞬变电磁法快速解释，在拟地震解释中还可用于速度分析实现精细成像；相关叠加算法利用了地空瞬变电磁法数据量大的特点，将传统的单点解释拓展到多点联合解释，巧妙的解决了二次微分带来的假异常现象，能够有效压制噪声、突出真异常；且由于该法计算速度快，可进行地空瞬变电磁法快速实时成像.

等效导电平面法是建立在全空间理论基础上的一种瞬变电磁解释方法，是根据视纵向电导曲线的特征值直观划分地层的一种近似解释方法(朴化荣，1990；李貅，2002).根据电磁理论，可用一导电平面来代替地下均匀介质，然后根据镜像法求出空间任一点的感应电磁场.当电源断开时，由于电磁感应，导电平面上产生涡流，为了求出涡流产生的电磁场，在导电平面的下方对称位置放置虚源代替导电平面中的涡流，在空间任一点某一时刻观测到的二次场响应信号就等同于对应某一深度的虚源产生的信号，在另一时刻观测到的感应信号又可以等同于对应另外一个深度的虚源产生的信号，这样，随着时间的推移，导电平面会上下“浮动”，从而实现空间任一点的感应电磁场计算.

设在地表铺设有m个电性源，AiBi表示第i个电性源，电流强度为I_i，源的长度为d_si，源A₁B₁沿X轴放置，中点位于坐标原点O，Z轴向下，整个XYZ坐标系满足右手坐标系，飞行器与地表的距离为z，在地表的投影为P，导电平面电导率为σ，距地表的距离为h，虚源Ai′Bi′关于导电平面与地表的电性源AiBi对称，见图 1.

图 1

Fig. 1

图 1 多辐射场源地空瞬变电磁法等效导电平面示意图 Fig. 1 Sketch of floating plate theory of multi-source ground-airborne transient electromagnetic method

对于水平层状介质，大地表面m个电性源在空中任一点处时间域磁场垂直分量响应的近似计算公式为

由(1)式可知，对于给定的g值，给出不同的研究深度H，就可确定一系列H_z(t)和t，进而得到H_z(t)近似计算曲线.可见，等效导电平面算法近似计算层状模型空中瞬变响应的计算精度主要由参数g决定.中心回线法的等效导电平面近似算法中，g值是根据解析形式的晚期瞬变响应公式确定的，通过这种方法得到的g值，计算晚期瞬变响应的精度尚可，但早期瞬变响应的误差很大(与解析解的相对误差高达20%)，最终影响了等效导电平面近似算法的计算精度(李貅，2002).

由于多辐射场源可视作单辐射场源情况的叠加，而单辐射场源又可视作电偶极子源情况的叠加，接下来以电偶极子源为例讨论g值的选取.在求解g的过程中我们发现，不同的g值，正演得到的瞬变场在双对数坐标系下与计算精度较高的滤波解呈现出不同的相对关系.图 2所示为均匀半空间模型下，改变g值后得到的瞬变响应与滤波解对比图，模型参数如下：源AB的长度10 m，电流大小100 A，均匀半空间模型电阻率为100 Ω·m，测点位置(300 m，400 m，-50 m)，图(a)为测点及源所在坐标系示意图，AB表示电性源，M表示空中的测点，P表示M在地表的投影，z是空中测点M距地表的高度；从图(b)中可以看出，改变g值对早期响应的影响较小，对晚期响应的影响较大，若想在全时域提高等效导电平面近似计算的精度，仅靠调节g值是不够的.但通过调整g值，可以找到一个在双对数坐标系下与滤波解线性程度最高的瞬变响应曲线，如果采用这个g值，再通过校正方法将对应的瞬变响应曲线校正到滤波解处，从而提高等效导电平面算法的正演计算精度.为此，将(1)式改写为

图 2

Fig. 2

图 2 均匀大地模型改变值对等效导电平面正演算法的影响 (a)测点及源所在坐标系示意图; (b)改变值对正演瞬变响应的影响. Fig. 2 Forward modeling with floating plate method of homogeneous earth with different g (a) Sketch of survey point and sources in the coordinate system; (b) Forward modeling results with different g.

这样，等效导电平面方法近似计算水平层状介质地表多辐射场源在空中任一点H_z(t)的问题就转化为(2)式中参数g和c的求解问题.阻尼最小二乘方法是一种得到广泛应用的最优化方法，这里将采用该法获取参数g和c.基于数据拟合和压制随机误差的考虑，建立如下所示目标函数

从图 2中可以大致判断出g和c的取值范围，但每次求解得到的解未必满足这一范围，为使优化参数满足约束条件，尽量减少多解性的影响，引入约束条件：给出g和c的上界和下界(此范围称之为可行域)，当解不在可行域范围内时，取与解更接近的g和c的上界或下界作为此次迭代得到的解.这样，每次得到的解都在可行域范围内，如此反复迭代运算，直至求出最优化的g和c，带约束的阻尼最小二乘法流程见图 3.

图 3

Fig. 3

图 3 阻尼最小二乘法求取最优化g和c流程图 Fig. 3 Flow chart for optimized g and c with damped least square method

最终得到的最优化参数g和c分别为：g=4.3，c=1.458.对应的多辐射场源地空瞬变电磁法等效导电平面算法近似计算磁场强度垂直分量的公式为

设计了四种典型层状模型对(4)式所示的等效导电平面正演算法进行验证，图 4为A、H、K、Q型模型的等效导电平面解与滤波解对比以及误差分布图，计算采用的源长度均为100 m，电流大小100 A，为不失一般性，两个电性源之间的夹角为30°，接收机高度-50 m，测点位置(300 m，400 m，-50 m)，A模型参数：ρ₁=100 Ω·m，ρ₂=200 Ω·m，ρ₃=300 Ω·m，h₁=30 m，h₂=20 m；H模型参数：ρ₁=100 Ω·m，ρ₂=50 Ω·m，ρ₃=100 Ω·m，h₁=30 m，h₂=20 m；K模型参数：ρ₁=100 Ω·m，ρ₂=500 Ω·m，ρ₃=100 Ω·m，h₁=30 m，h₂=20 m； Q模型参数：ρ₁=500 Ω·m，ρ₂=200 Ω·m，ρ₃=100 Ω·m，h₁=30 m，h₂=20 m；测点坐标示意图及模型俯视图见(a)和(b).由图 3可见，对于这四种典型地电模型，由上述等效导电平面法正演得到的近似解与滤波解吻合很好，相对误差均不足3%.相比较目前中心回线法等效导电平面近似算法高达20%的 相对误差，本文提出的算法有效提高了等效导电平面算法正演的精度.

图 4

Fig. 4

图 4 四种典型层状模型等效导电平面正演算法验证 (a)测点坐标示意图;(b)模型坐标俯视图;(c)A型模型等效导电平面解与滤波解对比图;(d)H型模型等效导电平面解与滤波解对比图;(e)K型模型等效导电平面解与滤波解对比图;(f)Q型模型等效导电平面解与滤波解对比图;(g)A型模型等效导电平面解与滤波解相对误差图;(h)H型模型等效导电平面解与滤波解相对误差图;(i)K型模型等效导电平面解与滤波解相对误差图;(j)Q型模型等效导电平面解与滤波解相对误差图. Fig. 4 Test of floating plate theory with four typical three-layered models (a) Sketch of survey point in the coordinate system; (b) Top view of sources and survey point; (c) Comparison between results of floating plate method and those of digital filtering method of A model; (d) Comparison between results of floating plate method and those of digital filtering method of H model; (e) Comparison between results of floating plate method and those of digital filtering method of K model; (f) Comparison between results of floating plate method and those of digital filtering method of Q model; (g) Relative error of forward modeling results of A model; (h) Relative error of forward modeling results of H model; (i) Relative error of forward modeling results of K model; (j) Relative error of forward modeling results of Q model.

根据(4)式，可得随时间变化的视纵向电导函数为

当该函数取极小时，值取得最优解.获得值后，结合(5)和(6)式即可求得S_τ(t)，再根据下式得到等效导电平面与地表之间的距离为

求得S_τ(t)和h_τ(t)后，对视纵向电导数据在深度进行一次和二次微分，得到视纵向电导微分成像数据.

视纵向电导的二次微分对地电模型的界面有很好的识别能力，曲线起跳的方向反映了相邻界面的电阻率变化趋势(高阻进入低阻时起跳为正，低阻进入高阻时起跳为负)，振幅则反映了相邻界面电阻率差异的大小，差异越大振幅越大(李貅等，2013；戚志鹏等，2015).经过二次求导的微分电导成像方法计算速度快，对电性界面反映灵敏，但是从数学的角度分析，微分运算是放大变化的过程，在这一过程中误差也会被放大.由于实际工程数据往往存在干扰和误差，视纵向电导微分成像会产生许多假异常，因此该法受干扰和误差的影响较大.

为了压制视纵向电导微分成像带来的假异常，突出真异常界面的振幅，可以采用相关叠加合成方法(李貅等，2012).该法对空中观测点周边一定范围内的观测信号进行相关叠加处理，对于相邻的测点i-1，i，i+1，得到的瞬变电磁信号具有相关性，测点间距越近，相关性越大.相关叠加合成算法正是利用相邻测点的相关性对视纵向电导微分成像曲线进行相关叠加，加强反映地电结构的真实异常，而干扰和误差的相关性较差，在计算中会得到压制，最终达到有效压制噪声、提高信噪比、加强有用信号、提高分辨率的目的.

两列数据的归一化互相关系数为

图 5是相关叠加合成示意图，以i点为中心，阐述n点合成的过程，首先分别计算i点与测点的归一化互相关系数，通过改变τ值得到最大相关系数，此时的即为最佳偏移深度，将与测点的数据相乘得到新的视纵向电导微分数据，i点的视纵向电导数据记作W₀，i点的自相关系数为1，偏移为0.n点合成即为按照最佳偏移深度叠加到中心点W₀上得到新的i点数据，即

图 5

Fig. 5

图 5 相关叠加合成示意图 Fig. 5 Sketch of correlation stack

为了验证多辐射场源地空瞬变电磁法微分电导成像算法，设计了如下三维模型：在电阻率100 Ω·m的均匀大地中，赋存两个电阻率为10 Ω·m的块状体，两个异常体的顶板埋深均为80 m，异常体尺寸大小均为60 m×60 m×100 m，块状体之间的距离为80 m，两个互相平行的发射源长度均为100 m，距两个块状体中轴线的距离均为100 m，电流大小30 A，电流方向相反，接收机位于空中100 m处，三维正演采用矢量有限元方法完成.模型示意图及俯视图见图 6，结合(a)和(b)可确定块状异常体的具体位置.

图 6

Fig. 6

图 6 三维模型示意图 (a)异常体模型示意图; (b)模型俯视图. Fig. 6 Model sketch (a) Anomalous body sketch (b) Top view of the model.

Line33(y=0 m)位于两个块状异常体的中轴线上，对这条剖面的正演数据进行全域视电阻率计算(张莹莹等，2015)，由图 7可见，由于三维模型中设计的两个块状异常体埋深较大且间距较小，在全域视电阻率断面图上只能看出存在低阻异常体，但无法分辨出这两个低阻异常体，也无法确定低阻异常体的位置.

图 7

Fig. 7

图 7 Line33剖面全域视电阻率成像结果 Fig. 7 Full field apparent resistivity imaging results of Line 33

Line33(y=0 m)测线的微分电导成像图见图 8，从图 8a中可以看出微分电导成像结果显示存在三个较明显的界面，第一个界面指示地表，由于空气和大地的电性差异很大，因此该界面的反应比较强烈，且遍布整个区域；第二个界面的范围较窄，指示两个低阻异常体的顶板位置；第三个界面应指示两个低阻体的底板位置，由于异常埋深较大，深度出现了偏差且微分电导波形出现了展宽现象.图 8b显示的是对三维正演数据加10%白噪后的微分电导成像结果，对比图 8a可以看出，微分电导结果受噪声的影响较大，只有指示地表的第一个界面的信号较强且完整，指示低阻异常体的信号受到较大的干扰，只能依稀从微分电导结果中分辨出异常体顶板和底板位置，此外，噪声还会带来诸多假异常干扰.

图 8

Fig. 8

图 8 Line33剖面微分电导成像结果 (a)无噪声;(b)加10%白噪. Fig. 8 Differential conductance imaging results of Line 33 (a) Without noise; (b) With 10% white noise.

采用五点合成对图 8所示的微分电导成像结果进行相关叠加合成处理，结果见图 9.由图可见，由于计算误差等带来的干扰、随机噪声等不具有相关性，因此进行相关叠加合成后均得到压制，而与电性差异变化较大的界面相对应的信号均得到了加强.

图 9

Fig. 9

图 9 Line33剖面相关叠加合成成像结果 (a)无噪声;(b)加10%白噪. Fig. 9 Correlation stack imaging results of Line 33 (a) Without noise; (b) With 10% white noise.

5.1    本文基于等效导电平面算法和相关叠加合成算法，介绍了一种适用于多辐射场源地空瞬变电磁法的快速成像方法.该法利用了视纵向电导二次微分对电性突变界面敏感以及计算速度快的特点，从电偶极子等效导电平面法原理出发，建立了多辐射场源地空瞬变电磁微分电导成像方法.针对地空瞬变电磁法数据量大、利于进行剖面解释的特点，引入相关叠加合成算法，将传统的单点解释拓展到多点联合解释；由于实际工程数据往往存在干扰和误差，视纵向电导微分成像会产生许多假异常，但来自同一地质体的信号具有很好的相关性，而噪声本身相关性很差，相关叠加合成算法利用信号的这一特点巧妙的解决了假异常现象，能够有效压制噪声、突出真异常.

5.2      文中基于等效导电平面方法推导的电偶极子磁场响应公式不仅可用于单辐射场源、多辐射场源地空瞬变电磁法快速解释；利用叠加原理，还可用于带偏移距的地面瞬变电磁法快速解释；在拟地震解释中还可用于速度分析解决精细成像问题；将均匀大地的电偶极子响应公式推广到层状模型中，就可以实现空中、地下、地表任意一点瞬变场的快速近似计算.

5.3     本文介绍的多辐射场源地空瞬变电磁法快速成像方法计算速度快、抗噪能力强、对电性突变界面反映灵敏，可实现实时成像，但文中仅对同一测线上的数据进行了相关叠加合成处理，在今后的工作中可以考虑进行不同测线、不同剖面之间的相关叠加合成，进行三维相关叠加合成计算，进一步提高分辨率.