电性源时域地空电磁数据小波去噪方法研究

李肃义，林 君，阳贵红，田培培，王 远，于生宝，嵇艳鞠

地球信息探测仪器教育部重点实验室，吉林大学仪器科学与电气工程学院，长春 130061

1 引 言

目前，我国矿资源供需矛盾日益突出，巨大的资源短缺制约经济社会发展，“攻深探盲”战略是解决此矛盾的主要途径也是资源勘查所面临的首要难题.时间域地空电磁法（Time－domain Ground－Airborne Electromagnetic Method）是融合地面 TEM 与航空TEM优势的一种新型勘探方法，通常将发射系统放置于地面，并铺设几公里的长接地导线或大定源回线源，然后将接收系统、传感器安装在直升机、无人机或飞艇上进行飞行测量.此方法不仅具有航空TEM的空间分辨率高、野外布线方便快捷、探测高效等优势，还具有地面TEM的大发射磁矩、信噪比高、勘探深度大的优势，是深部矿产资源快速勘查的重要技术手段［1－5］.

国外从上世纪90年代初开始研究地空电磁勘探技术［2－5］，Tohru［6］，将基于直升机的地空电磁探测系统成功应用于 Mount Bandai火山结构勘查；Hisatoshi等［7］成功应用直升机地空电磁系统探测到了800m深处的地下水资源.我国地空电磁探测研究起步较晚，吉林大学的研究人员于2009年开始研究地空电磁探测方法，由于直升机勘探费用高、飞行员及飞机在飞行过程中存在风险，所以，研究小组研发了基于无人飞艇的地空电磁勘探系统［8－9］.基于飞艇的地空电磁勘探方法不仅可以解决直升机航空飞行勘探困难的问题，而且适用于地形复杂的山区资源探测，在我国深部矿产资源探测和地质普查应用中具有广阔的发展前景.

但是，在基于飞艇的地空电磁勘探中，飞艇会受风向、大气气流、地形、地面局部温度场变化等影响，导致飞行高度、航迹、姿态等发生变化.这些变化使得固定在艇囊前端的接收线圈切割大地磁场，引起磁通量改变而产生感应电动势.这部分感应电动势叠加在大地的电磁感应电动势曲线上，会令测得的电磁信号发生严重的基线漂移，降低电磁数据的信噪比，影响电磁数据的电阻率成像质量.线圈运动引起的噪声具有幅值较大、分布范围广、频率低的特点，是地空电磁信号中的主要噪声之一.此外，各种天然磁场及人文设施也会产生较大的噪声与干扰，如天电噪声、50Hz的工频干扰等，这些噪声会致使数据反演解释时出现假异常，干扰地下电性结构的正确识别，进而影响数据解释的准确度与精度.

目前时间域地空电磁探测方法仍处于起步阶段，直接针对地空电磁信号去噪的相关研究较少.不过，近年来，随着计算机以及信号处理技术的飞速发展，信号去噪新理论、新方法不断涌现，并且国内外学者在多种航空电磁数据噪声去除方法的研究方面取得了很大进展［10－14］.Fugro航空地球物理探测公司则在硬件上增加线圈运动补偿，改善悬吊稳定性，并在数据处理中采用数据叠加以及高通滤波器等方法对运动噪声进行去除.但是，当采用以上方法时，需要增加分量测量或姿态记录硬件电路；Lemire等［12］采用样条插值和拉格朗日优化方法，对时间域航空电磁数据中的基线漂移进行了校正，但计算过于繁琐；对于天电噪声与白噪声，Abderrezak Bouchedda等［13］采用平稳小波变换对航空电磁数据中的天电噪声进行了抑制；Reninger等［14］利用奇异值分解方法，对时间域航空电磁数据中混入的天电噪声及白噪声进行了处理.

本文针对基于飞艇的时间域地空电磁数据中的主要噪声，借鉴时间域航空电磁数据的消噪方法，在分析了地空电磁信号特点与主要噪声时频特性的基础上，提出一种综合小波消噪法：首先使用小波高尺度近似分量估计基线漂移，极大限度地消除基线漂移的干扰，再使用小波阈值收缩法有效地消除其余噪声的影响.通过异常环模型的理论响应和实测数据实验，结果表明该方法对多种噪声均有很好的抑制作用，是一种实用有效的时间域地空电磁数据消噪方法.

2 综合小波消噪法

无论飞行载体是直升机还是无人机或飞艇，在空中进行飞行测量时，测量的电磁信号中均会混入多种噪声，主要包括：飞行时摇摆导致接收线圈产生的低频类基线漂移、分布范围广的随机噪声、天电噪声及工频干扰等.白噪声、天电噪声及工频干扰可以通过小波阈值收缩各尺度上的小波系数进行抑制，但基线漂移的能量主要集中在高尺度的近似分量上，所以使用阈值收缩法往往不能很好地去除基线.为了弥补这一不足，根据基线漂移类噪声特点，有两种改善方案，一是提高阈值法中的分解级数，二是使用两种算法相结合的思路.

增加阈值收缩法中的小波分解尺度，使用最高尺度的近似分量估计基线漂移，然后收缩各级的小波系数，重构时使用处理过的小波系数，并将最高尺度的近似分量置0，这种算法可以达到既消除基线漂移又消除其他噪声的目的.但是，随分解尺度增多，阈值处理尺度也会增多，有用信息损失的风险就会增加，重构时误差也会增大，甚至消噪后的信噪比不升反降.并且，尺度增多计算量会增大，处理速度也会降低.所以，通过理论与实践证明，高尺度分解的阈值收缩法并不适合地空信号的消噪.

因此，本文提出使用两种方法相结合的综合消噪法：首先使用一种方法极大限度地消除基线漂移的干扰，再使用另一种方法有效地消除其余噪声的影响.

2.1 小波基的选取

不同的小波基，具有不同的时频特性，所以使用不同的小波基进行地空电磁信号的消噪，结果可能相差甚远.

根据Mallat的最优小波基选取准则：只需用较少数目的非零小波系数即可精确地重构信号［15］.我们使用常用小波基对地空电磁信号做了多次分解重构实验，并结合它们的数学特性以及要处理的地空电磁信号特点，最终选取sym8作为综合小波消噪算法中的小波基.

2.2 基于多分辨率分析的基线校正算法

基线漂移的能量主要集中在高尺度的近似分量上，几乎没有奇异点存在，而阈值收缩法的分解尺度不能过大，并且主要是对小波系数进行处理，不对近似系数进行处理，所以不适合校正信号中的基线漂移.

利用小波的多分辨率分析特性可以将信号逐级分解成低频信息与高频信息，将交织在一起的各种不同频率组成的混合信号分解成不同频带的子信号.低频信息对应着信号的近似部分，即尺度系数；高频信息对应着细节部分，即小波系数［11－13］.本文地空数据的采样率为51.2kHz，根据采样定理，信号中包含的频率信息为0～25.6KHz.由于基线漂移频率较低，根据本文中每道采集的数据点数，确定分解级数为12级，即利用sym8小波12级分解所得的近似分量作为基线漂移的估计，从原始信号中减去这部分低频信息，所得信号即为校正信号.

2.3 基于小波收缩法的地空信号去噪算法

小波收缩法对于频率范围分布广的白噪声，以及与有效信号发生频率重叠的噪声具有很好的去除效果.将实际信号X写成向量形式有

S为有效信号，e为噪声，目的是从X中去除噪声e，估计出S.若e为满足独立同分布（i.i.d.）或 N（0，σ2）的噪声，则可以使用正交离散小波变换矩阵W来估计信号S：

其中w，θ，与z分别是X，S与e的正交小波变换.由于e已被假设为白噪声，所以变换z也为白噪声，此时，小波系数可表示为

我们的目标是估计θ具有最小的L2风险，最小均方误差MSE，即，找到最佳的估计具有最小L2风险：

（1）分解层数的选取

阈值收缩法中分解层数的选取直接影响到消噪效果.随着小波分解层数的增多，得到的细节信息就越丰富.但是，并不是层数越多越好，因为小波变换中分解合成实际上是对频带的划分，分解层数越多，产生的子带越多，频带划分得越细，由于上一级频带分解的信号输出又作为下一级频带分解的输入，层数增加意味着级间的滤波器增多，造成信号移位也越大；另一方面小波分解合成都要进行边界延拓，层数越多引起边界失真越大.而这些损失的信息量均是小波逆变换不能恢复的损失，因此基于小波分解的层数不宜过高，根据地空有效信号的频带，本文选用5层小波分解.

（2）阈值与阈值收缩策略

根据电磁信号特点，对比了多种阈值策略和阈值收缩函数后，本文应用极小极大（Minimaxi）原理对每一尺度上的小波系数计算相应的Minimaxi阈值.Minimaxi原理是一种找出失败的最大可能性中最小值的算法.Minimaxi阈值相比其他阈值略小，采用这个阈值策略主要考虑让信号尽可能小的减少有效信息的损失.阈值收缩函数也各有利弊，由于电磁信号具有衰减特性，为此我们采用硬阈值收缩函数（Hard），Hard收缩函数最适于保持信号的奇异点.

3 理论模型算法测试

由于闭合异常环模型不仅可以模拟有限导体的晚期电磁响应，而且在野外便于铺设.所以本文首先基于闭合异常环模型计算的理论响应，对算法压制噪声和校正基线的有效性进行验证.

采用闭合异常环模型，进行电性源地空电磁探测实验的示意图如图1所示.

闭合异常环用一个等效电感L和等效电阻R进行近似，当发射线圈中通有交变电流时，在闭合异常环中产生的感应电流可以写为

式中iAL为闭合异常环的感应电流，LAL为线圈的电感，MTL为发射导线与闭合异常环的互感.τL为闭合异常环的时间常数，与线圈的材料、电气性能有关，其表达式为

式中μ0为磁导率，n为闭合异常环匝数，rAL为线圈半径，ra为导线半径，ρL为导线电阻率.当发射电流为梯形波时，在闭合异常环的取样单元中产生的感应电压信号为

图1 基于异常环模型的电性源地空电磁探测示意图Fig.1 Schematic diagram of electrical－source groundairborne electromagnetic detection system based on anomaly loop model

式中MRL为接收线圈与闭合异常环的互感，I为发射电流幅值，T1为电流上升时间，T2为电流下降时间，T3为发射电流宽度.

对于任意位置、边长为2a×2b的矩形异常环与半径为r2的圆形接收线圈之间的互感系数，可以写为

式中（Tx，Ty，Tz）为接收线圈的圆心坐标.

对于长导线与异常环的互感可以通过简化公式（9）得到.基于吉林大学自主研制的电磁探测系统，工作参数如表1所示，采用高斯数值积分方法，利用公式（6）和（8），就可以计算出闭合异常环的电磁响应［17－18］.

表1 异常环地空实验工作参数表Table 1 Experimental work parameters of ground－airborne electromagnetic detection based on anomaly loop model

将计算的异常环电磁响应作为有效信号s（n），如图2a所示.采用信噪比为3dB的Gaussian白噪声与低频信号来模拟电磁噪声和运动噪声e（n），在s（n）上叠加e（n）来模拟含噪信号x（n），由于线圈运动引起的噪声呈现低频特性，为了清晰看到基线漂移，我们同时显示了多个周期的信号，如图2b所示.

图2 电性源地空测量时异常环的电磁响应（a）理论信号；（b）模拟含噪信号.Fig.2 The electromagnetic response of electrical source ground-airborne based on anomaly loop model（a）Theoretical signal；（b）Simulated noisy signal.

图3 消噪方法比较（a）小波多分辨率分析法校正基线；（b）小波阈值收缩法消噪；（c）小波综合法消噪.Fig.3 Comparison of different wavelet-based de-noising methods applied to the electromagnetic response（a）Baseline drift correction by using estimation of high-level approximation；（b）De-noising by wavelet shrinkage；（c）Denoising by the combined wavelet method.

对于图2b中的噪声信号，如果单独使用小波多分辨率分析法校正基线，则校正后的信号如图3（a）所示；单独使用小波阈值收缩法消噪后的信号如图3b所示；使用小波综合法校正基线、消除噪声后的结果如图3c所示.通过图3我们可以看到，综合法不仅很好地校正了基线，而且较好地压制了大部分的噪声.

此外，我们还通过均方误差（Mean Squared Error，MSE）或误差范数（Error norm）来量化评估消噪方法的有效性.MSE定义为

MSE的值越小，则两信号的相似度就越高；Error norm的值越小，表明两信号的误差越小.表2给出了模拟含噪信号与综合法消噪信号的均方误差值与误差范数.可以看到，综合法消噪后有效降低了MSE值与Error norm值，消噪后的信号与理论信号的相似度高、误差小.

表2 信号消噪前后MSE与Error norm比较Table 2 Comparison of MSE and Error norm of noisy signal and de－noised signal

4 野外测数据消噪验证

2010年11月，在长春市农安县烧锅镇化东屯，采用吉林大学自主研制的电性源发射系统进行地面发射，利用吊车将地空接收系统和接收线圈吊在空中5m高，开展了低阻异常模型的飞行模拟实验.具体工作参数同表1，接收系统采用24位AD采集、100KHz采样频率、全波形连续采集、GPS同步控制以及无线遥控方式.实验中，我们首先测量了背景场，然后铺设了异常环（匝数为10m×10m×5匝、时间常数为4.62ms），位置如图4所示，最后沿发射电性源方向测量了4次，分别如图4中R1－R4测线所示.

图5a为异常环模型的野外实测数据，图5b为使用综合小波法去噪后的信号，通过对比两信号，可以看到，该方法在失真度很小的前提下，很好地校正了基线.

图6为R3测线未减去背景场的原始数据剖面曲线，图7为减去背景场并进行综合小波消噪后的剖面曲线，在测线的0～55m、78～120m均为背景场区域，由于没有异常存在，所以各道数据基本重合；但在测线的56～77m处，对应异常环位置，剖面曲线变化明显.图8为理论计算的异常环电磁响应剖面曲线，可以看到，图7所示的剖面曲线与图8中理论计算的剖面曲线变化规律一致.但是，从图7和图8中可以观察到异常宽度存在一定的差异，主要原因是无法记录实验采用吊车在测量时的转动速度、吊车转动通过异常环的时间，只能估算测点的距离，因此导致了实测数据中的异常宽度与理论计算的异常宽度存在误差，此误差可以在计算视电阻率断面图时进行修正.

图4 加异常线圈测试场地工区图Fig.4 Schematic diagram of wild field laid an anomaly loop

在对实测数据进行视电阻率计算时，根据场源等效思想，将长导线源剖分为若干个电偶极源，基于均匀半空间模型求解了地面LOTEM的全区视电阻率，然后进行了视电阻率－深度计算［9］.图9为采用原始数据计算的视电阻率断面图，图10为采用综合小波消噪法后计算的视电阻率断面图.烧锅镇化东屯地质概况属于长春地区第四系结构，地表为黄土覆盖层，基岩为白垩系地层，以砂砾岩和泥岩为主，沙层、泥岩和泥沙岩等交替出现，胶结微弱孔隙较发育.距地表13m为黄土状亚粘土低阻盖层，地下50m至70m之间为含水砂岩，电阻率值在20～50Ωm之间变化，电阻率变化与钻井资料一致.从图9中可以看到由于基线漂移与噪声的存在，导致解释结果中异常倾斜，而且即使减掉大地的电磁响应，也出现了虚假的低阻异常；而图10中由异常环引起的低阻异常清晰，且在浅层和深层异常一致.

5 结 论

5.1 本文根据地空电磁信号和主要噪声的特点，通过分析常用小波系的数学特性与比较信号重构误差，提出了一种基于sym8小波基的综合小波去噪法对地空电磁信号中混入的多种噪声进行消除.针对基线漂移类低频噪声，先基于小波多分辨率分析，将原始信号进行多层小波分解，通过计算最高分解层的近似分量来模拟电磁数据中的基线漂移，再从原始数据中去除该近似分量以达到消除基线漂移的目的.针对频率分布广的白噪声，应用MiniMaxi阈值配合硬收缩函数进行5级小波收缩进行消噪，既保证了数据波形的失真度小，又较好地抑制了具有白噪声特征的多种电磁噪声.

5.2 基于异常环模型进行理论电磁响应仿真，比较了综合小波去噪法消噪前后信号的MSE与波形失真度.并将算法用于了长春市烧锅镇野外异常环模型实测数据的处理，比较了消噪前后叠加所成的视电阻率断面图.结果表明这种综合消噪法对多种噪声均有很好的抑制作用，是一种实用有效的时间域地空电磁数据消噪方法.

（References）

［1］Fountain D.60years of airborne em－focus on the last decade.AEM2008－5th International Conference on Airborne Electromagnetics Haikko Manor，Finland，2008：1－1.

［2］Elliott P J.New airborne electromagnetic method provides fast deep－target data turnaround.The Leading Edge，1996，15（4）：309－303.

［3］Elliott P J.The principles and practice of FLAIRTEM （Fixed Loop Airborne Electromagnetics）.Exploration Geophysics，1998，29（2）：58－60.

［4］Smith R S，Annan A P，McGowanz P D.A comparison of data from airborne，semi－airborne，and ground electromagnetic systems.Geophysics，2001，66（5）：1379－1385.

［5］Mogi T，Tanaka Y，Kusunoki K，et al.Development of grounded electrical source airborne transient EM （GREATEM）.Exploration Geophysics，1998，29（2）：61－64.

［6］Mogi T，Kusunoki K，Kaieda H，et al.Grounded electricalsource airborne transient electromagnetic （GREATEM）survey of Mount Bandai，north－eastern Japan.Exploration Geophysics，2009，40（1）：1－7.

［7］Ito H，Mogi T，Jomori A，et al.Further investigations of underground resistivity structures in coastal areas using grounded－source airborne electromagnetics.Earth，Planets and Space，2011，63（8）：e9－e12.

［8］嵇艳鞠，林君，许洋铖.大定源时间域吊舱式半航空电磁勘探理论研究.第九届中国国际会地球电磁学术讨论会论文集，2009：121－126.Ji Y J，Lin J，Xu Y C.Semi－airborne time domain electromagnetic exploration theoretical study using a largefixed－loop （in Chinese）.Collection of the Ninth China International Geo－Electromagnetic Workshop，2009：121－126.

［9］Ji Y J，Yang G H，Guan S S，et al.Interpretation research on electrical source of time domain ground－airborne electromagnetic data.World Automation Congress （WAC），2012，4（6）：1537－1540.

［10］李貅，宋建平，马宇等.基于小波分析的TEM信号提取.煤田地质与勘探，2005，33（2）：72－75.Li X，Song J P，Ma Y，et al.The abstract of TEM signal based on the wavelet analysis.Coal Geology ＆ Exploration（in Chinese），2005，33（2）：72－75.

［11］武粤，孟小红，李淑玲.小波分析及其在我国地球物理学研究中的应用进展.地球物理学进展，2012，27（2）：750－760.Wu Y，Meng X H，Li S L.Wavelet analysis and its application in geophysics of China.Progress in Geophysics（in Chinese），2012，27（2）：750－760.

［12］Lemire D.Baseline asymmetry，Tau projection，B－field estimation and automatic half－cycle rejections： THEM Geophysics Inc.Technical Report，2001：1－13.

［13］Bouchedda A，Chouteau M，Keating P，et al.Sferics noise reduction in time－domain electromagnetic systems：application to MegaTEMIIsignal enhancement.Exploration Geophysics，2010，41（4）：225－239.

［14］Reninger P A，Martelet G，Deparis J，et al.Singular value decomposition as a denoising tool for airborne time domain electromagnetic data.Journal of Applied Geophysics，2011，75（2）：264－276.

［15］Mallat S.A Wavelet Tour of Signal Processing （Second Edition）.Beijing：China Machine Press，2003.

［16］Donoho D L. De－noising by soft－thresholding.IEEE Transactions on Information Theory，1995，41（3）：613－627.

［17］嵇艳鞠，李肃义，于生宝等.基于异常线圈的时间域AEM系统测试和标定方法研究.地球物理学报，2011，54（10）：2690－2697.Ji Y J，Li S Y，Yu S B，et al.A study on time－domain AEM testing and calibration method based on anomaly loop.Chinese J.Geophys.（in Chinese），2011，54（10）：2690－2697.

［18］嵇艳鞠，栾卉，李肃义等.全波形时间域航空电磁探测分辨率.吉林大学学报 （地球科学版），2011，41（3）：885－891.Ji Y J，Luan H，Li S Y，et al.Resolution of full－waveform airborne TEM.Journal of Jilin University （Earth Science Edition）（in Chinese），2011，41（3）：885－891.