1 引言

航空电磁法是航空物探常用的测量方法之一，具有速度快、成本低、探测范围广、可用于海域探测等优势.20世纪90年代以来，时间域直升机电磁探测系统以其操作灵活机动，水平分辨率高得到了迅猛发展，并广泛应用于矿产资源探查及地质填图^[1~3]等，如加拿大的HeliGEOTEM，THEM，VTEM，AeroTEM，南非的ExplorHEM，美国的NewTEM，澳大利亚的HoisTEM，以及丹麦的SkyTEM^[4~7]等.

我国山区较多，地形复杂，西部荒漠人烟稀少，工作条件恶劣，给地面布线、勘探带来很多困难，因而航空电磁法作为一种快速经济的勘查工具在我国有着广泛的应用前景，如我国西部地区以铜为主的金属矿床，我国西北、华北等干旱地区的区域水文物探调查等.我国航空电磁法研究始于20世纪50年代末，发展并不顺利，与国际发展水平还有相当大的差距.近年来，国家大力资助具有自主知识产权的航空物探仪器系统的研制以及航空物探理论的研究.

电导率深度成像是一种近似反演方法，利用在时间域，无论介质的形状和电性如何，其电磁响应都可以用衰减指数总和表示这一特性，计算出随深度变化的视电导率.而反演是求目标函数最小值的过程，通过初始模型和迭代计算来实现.尽管一维反演在反演精度上优于电导率深度成像，但由于航空电磁数据数据量大，一维反演速度太慢，而电导率深度成像因其计算速度的优势以及可接受的成像精度，被世界许多先进的航空地球物理勘探公司，如Geotech，Aeroquest，Fugro等广泛使用.

目前，从时间域获得电导率深度成像（CDI，Conductivity Depth Imaging）的方法主要有; Maxwell镜像源计算法^[8~10]、基于最大电流层在半空间模型中深度的近似成像法^[11~14]、采用卷积或解卷积的方法、从任意发射电流波形off-time数据计算得到阶跃发射电流响应方法^[15~19]，采用薄板差分方法^{[15, 20, 21]}、以及基于假层半空间模型的查表法（Huang，2008）^[22]等.

本文采用Huang^[22]的假层半空间模型，基于直接数值积分与G-S变换的正演计算方法，通过神经网络的构建、训练及测试，研究时间域航空电磁数据的电导率深度成像算法，并给出反演结果.

2 一维正演算法

时间域直升机航空电磁探测采用如图 1所示的中心回线方式，吊舱（发射与接收线圈）距直升机35m，吊舱（发射线圈Tx、接收线圈Rx）距地面h（单位m），线圈中心为坐标原点，第i层大地的厚度、电导率以及磁导率分别为h_i，σ_i，μ_i.

2.1 中心回线方式垂直磁场分量的计算公式

通有e^iωt电流半径为a的发射线圈置于层状大地上方h处（图 1），线圈中心处垂直磁场分量^[23]为

(1)

其中，第i层的输入阻抗，固有阻抗，底层无限半空间的固有阻抗，其中，σ_i，μ_i分别为第i层的电导率和磁导率，ω为角频率.

用拉氏变换变量s代替（1）式中的iω，得到频率域响应的拉氏表示H_z（s）=H_z（0，h，s），除以s可得到上阶跃瞬态响应H_z（s）=H_z（s）/s.

2.2 接收线圈感应电动势的计算

野外实测的感应电动势V与磁场分量H_z对时间t的一阶导数成正比

(2)

式中B为磁感应强度（单位T），S为接收线圈的有效面积（单位m2），空气的磁导率μ₀=4π×10^-7（单位H/m），式中（单位）又称为归一化感应电动势.

根据拉氏变换的时域微分性，=-Sμ₀H_z（s），通过拉氏逆变换，可得到感应电动势的上阶跃响应，取其负值，得到感应电动势的下阶跃响应.因此发射电流关断后归一化的感应电动势

(3)

为求得瞬态电磁响应，采用G-S变换^[24]计算（3）式的拉氏逆变换，即

(4)

其中，

J是取决于计算机字长的偶整数（这里取16），M=J/2，m是（j+1）/2的整数部分.

将（1）式带入（4）式，交换求和、积分顺序，对被积函数求和计算G-S变换后，再采用高斯数值积分法^[25]，计算汉克尔变换.

由此，周期为T的双极性方波产生的电磁响应V（t）为

其中t_d为脉冲宽度，N为一个足够大的整数，这里取2.

2.3 正演计算结果

根据目前国内正在研制的吊舱式直升机时间域航空电磁探测系统的技术指标，计算了接收线圈内的归一化感应电动势，使用参数如下; （1）测量装置为中心回线方式，接收线圈位于发射线圈中心，发射线圈直径15m，发射电流1A，接收线圈等效面积1m2；（2）发射电流为25Hz正负方波，脉宽10ms；（3）方波的断电延迟时间为零，采样时间按等比序列分布于发射电流关断后的0.005~8ms之间，15个采样点.

图 2为吊舱高度30m时，半空间模型（0.025~20S/m）的1~7道电磁响应，可以看到，每一道的响应形态相近，在低电导率段，电磁响应随电导率增加而增大，在高电导率段，电磁响应随电导率增加而减小.正是由于这种二值性，单一通道数据反演具有很大的不确定性.图 3给出了5μs时刻，电导率为0.02S/m均匀半空间模型在不同吊舱高度情况下的电磁响应，当吊舱高度从30m变为34m（与H=30m曲线太近，略去，可从图中看出趋势）时，电磁响应从24640nT/s衰减到19365nT/s，也就是，吊舱高度增大了13%，归一化感应电动势减小了21%.由此可见，吊舱高度对航空电磁响应的影响较大，这也是航空电磁探测与地面探测方法重要区别之一，因此在飞行测量过程中，应准确记录吊舱高度，并且在反演解释中引入吊舱视高度.

3 基于神经网络的电导率深度成像算法 3.1 人工神经网络模型

人工神经网络方法因其学习记忆和非线性逼近能力在地球物理反演中得到广泛应用^[26~29].BP反向传播网络（Back_PropagationNetwork）^[30]是对非线性可微分函数进行权值训练的多层前向网络，分为输入层、隐含层和输出层，各层之间为权系数连接，且每一层的输出都是下一层的输入，可含有一个或多个隐含层.不同层有不同的节点（神经元）数目，每个神经元都带有一个输入为1的偏差值（亦称阈值）θ.以三层网络为例（图 4），输入向量为X=（x₁，x₂，x₃，x₄）^T，节点数n=4，隐含层输出向量为Y=（y₁，y₂，y₃）^T，节点数m=3，输出层输出向量为O=（o₁）^T，节点数k=1，输入层到隐层之间的权值矩阵用V表示，隐层到输出层之间的权值矩阵用W表示，对于输出层有; o_k=f（net_k），net_k=，对于隐层有; y_j=f（net_j），net_j=，转移函数f（x）采用双极性Sigmoid函数; f（x）=.网络误差是各层权值wjk，vij的函数，因此调整权值可改变误差，误差梯度下降为权值调整原则.

3.2 基于BP神经网络的电导率深度成像算法

Huang（2008）^[22]基于假层半空间模型，采用查表方式，实现了航空电磁数据的电导率深度成像.吊舱高度对航空电磁响应的影响较大，假层的引入解决了高度计以及地形变化引起的误差，比半空间反演模型具有很大的优势（见文献[22]中图 5a）.

本文电导率深度成像算法均采用Huang的假层半空间反演模型，如图 5所示，电导率为零、厚度为Δ的假想层代表树木、楼房等人文圈层.假想层的厚度Δ为反演解释吊舱高度h与飞机高度计值a之差.

电磁响应的二值性（图 2）导致单一通道数据反演成像具有很大的不确定性，这里将任意两道（第i道，第i+k道）的航空电磁数据按式（1）、（2）变换为时间常数τ_i和响应幅度α_i，

(1)

(2)

其中，t_i是第i道的中心时间，A_i为第i道的归一化感应电动势，k为采样道差.

利用前述正演算法，形成分别在23个吊舱高度（h=20~60m）下，64种电导率半空间模型（σ=0.0025~20S/m）在关断后0.005~8ms的15道时间域电磁响应集.若道差k=1，根据式（1）、（2），可形成14个“时间常数τ_i和响应幅度α_i-半空间模型”的正演样本集，每个样本集包含23×64对τ_i和α_i以及与之对应的吊舱高度h和大地的电导率σ.

为实现自动CDI成像，设计了14个具有双隐层的BP神经网络，网络输入变量分别为时间常数τ_i和响应幅度α_i，网络输出即为吊舱视高度h_i与大地视电导率σ_i.14个正演样本集分别用于训练这14个神经网络，根据正演样本集中网络输入、输出变量的变化规律，对网络输入响应幅度α_i以及网络输出视电导率σ_i分别采用以10为底和以2为底的对数归一化方法，网络输入时间常数τ_i和网络输出视吊舱高度h_i采用线性归一化方法.并对14个含有23×64个样本的正演样本集进行欧氏距离筛选，选出相关性小的训练样本，以提高神经网络训练效率.在比较各训练算法后，选用了RPROP训练算法，收敛快且精度高.当测试样本反演精度在5%以内，认为网络训练成功.

这样，任意15道航电数据计算出的14对时间常数τ_i和响应幅度α_i，分别输入至已训练好的14个神经网络，将得到14对网络输出吊舱视高度h_i和大地视电导率σ_i.视深度d_i由扩散深度δ_i（式（3））与假层视厚度Δi（Δi=h_i-a），根据Huang的实验规律（文献9中的图 4）确定.

(3)

3.3 电导率深度成像结果

为评价神经网络电导率深度成像算法的反演效果，设计了与文献[9]中相同的两个三层模型（表 1）.计算两个模型的15道电磁数据，并通过式（1）、（2）得到14对时间常数τ_i和响应幅度α_i，依次送入已训练好的神经网络，并根据扩散深度δ_i和吊舱视高度h_i确定视深度d_i.图 6、7给出了模型1、2采用神经网络法与查表法的CDI反演结果，其中，直线为真实模型，实心点线为基于假层半空间反演模型获取的CDI结果，空心点线为基于半空间反演模型的CDI结果.从图 6可以看到，查表法和神经网络法得到的CDI对中低阻覆盖层以及地下高导层的成像能力接近，对200mS/m的高导层，查表法计算视电阻率结果约为80mS/m，神经网络法的计算结果约为100mS/m，而对高阻基岩的分辨神经网络法优于查表法，对2mS/m的基岩，查表法计算视电阻率的结果约为30mS/m，而神经网络法的计算结果为10mS/m.模型2的CDI结果进一步地说明神经网络法计算视电阻率的优势，图 7a中查表法的CDI反演结果几乎是一条直线，很难看出20m厚的高阻夹层，图 7b中采用神经网络方法的CDI反演结果能看出地下高阻层，但仍与实际高阻层的电导率存在很大差异.这也与时间域电磁法探测原理有关，该方法对高导体分辨能力强.

从图 6、7可以看到，神经网络方法得到的CDI结果存在两个问题; （1）CDI曲线不光滑；（2）成像反演深度偏大.CDI曲线不光滑是由于神经网络训练程度不同所致，不同深度的视电导率来自于不同的神经网络输出.反演深度偏大，因为本文视深度的计算采用Huang（2008）^[22]针对发射基频90Hz，发射线圈直径5m时得到的经验公式，本文发射基频为25Hz，发射线圈直径15m，因而造成反演深度有一定的偏差，需要通过大量合成数据的成像实验对Huang的视深度经验规律进行系数调整.本文经过30多组模型数据成像实验，初步得到调整系数为0.92.若采用该系数对图 6b，7b中的成像深度进行调整后，成像结果将与真实模型更为接近，但因该系数仍需通过大量合成数据的成像实验进行研究验证，因此这里成像结果中并没有采用此修正参数.

神经网络方法成像是通过大量正演样本训练神经网络，使其“输入-输出”关系逼近“航空电磁响应-地球物理模型”这一非线性映射关系，也即逼近麦克斯韦偏微分方程组的映射关系，反演精度很大程度上取决于训练样本的选择以及网络训练成熟程度.如果不考虑海量数据反演的计算时间，一维反演是最好的解释方法，在某种程度上神经网络方法成像“类似”于一维反演，但计算速度上却有很大不同，一维反演在反演过程中需要不断调整模型参数进行多次拟合迭代，计算速度较慢，而神经网络成像方法在成像前需要大量时间进行神经网络训练，当测试样本达到成像精度后，网络训练结束，电磁数据输入神经网络并自动进行反演成像，计算速度快；查表法没有逼近“响应-模型”关系的过程，单纯去查找相近“响应”值并通过插值确定“模型”，其精度取决于数据库搭建精度，以及插值算法，并且每次反演成像都需要进行查询和插值.

CDI适用于具有层状地电结构大地的反演成像，但当电磁数据为二维或三维时，CDI具有鲁棒性，仍能得出正确成像结果（文献[22]的图 8，图 9）.神经网络方法的CDI对非一维电磁数据也具有鲁棒性，因为BP神经网络具有鲁棒性，训练神经网络的“响应-模型”数据即为查表法表格中的数据，具有鲁棒性，用鲁棒性的数据训练鲁棒性的神经网络，得到的神经网络也具有鲁棒性.BP神经网络的鲁棒性，是在构建神经网络时通过选用较为平缓的节点激发函数，充分训练网络以减小输出误差，以及合理的隐层数目设置等条件来保证的.本文没有构造文献9中的图 8、图 9描述的准二维模型，因其即为一系列一维模型的组合，而一维电磁数据成像在本文的图 6、图 7中已有论述.

4 结论

本项研究基于Huang的假层半空间反演模型，采用BP神经网络方法，解决了航空电磁数据电导率深度成像问题，主要结论;

（1）直升机时间域电磁响应随时间呈指数状衰减，对于低端电导率的均匀半空间模型，电导率越大，电磁响应相应越大，而电导率较高时，电导率越大，电磁响应反而越小，电磁响应具有二值性，并且电导率越大，衰减越缓慢；

（2）引入电导率为零、厚度为Δ的假想层代表树木、楼房等人文圈层，解决了高度计以及地形变化引起的误差；

（3）神经网络CDI成像具有逼近“响应-模型”非线性关系的特性，也即逼近电磁响应满足的麦克斯韦偏微分方程组的映射关系，因此能更准确地进行视电导率计算.

本文以假层半空间模型为基础，采用神经网络实现了航空电磁数据的电导率深度成像，能够提供较为准确可靠的地下电导率分布信息，简便易行.从合成数据的成像效果来看，神经网络方法得到的视电导率比查表法的更为接近实际模型，但在网络训练程度控制以及视深度的计算等方面仍需要进一步改进.