李肅義，林 君，陽貴紅，田培培，王 遠，于生寶，嵇艷鞠

地球信息探測儀器教育部重點實驗室，吉林大學儀器科學與電氣工程學院，長春 130061

1 引 言

目前，我國礦資源供需矛盾日益突出，巨大的資源短缺制約經(jīng)濟社會發(fā)展，“攻深探盲”戰(zhàn)略是解決此矛盾的主要途徑也是資源勘查所面臨的首要難題.時間域地空電磁法（Time－domain Ground－Airborne Electromagnetic Method）是融合地面 TEM 與航空TEM優(yōu)勢的一種新型勘探方法，通常將發(fā)射系統(tǒng)放置于地面，并鋪設幾公里的長接地導線或大定源回線源，然后將接收系統(tǒng)、傳感器安裝在直升機、無人機或飛艇上進行飛行測量.此方法不僅具有航空TEM的空間分辨率高、野外布線方便快捷、探測高效等優(yōu)勢，還具有地面TEM的大發(fā)射磁矩、信噪比高、勘探深度大的優(yōu)勢，是深部礦產(chǎn)資源快速勘查的重要技術手段［1－5］.

國外從上世紀90年代初開始研究地空電磁勘探技術［2－5］，Tohru［6］，將基于直升機的地空電磁探測系統(tǒng)成功應用于 Mount Bandai火山結構勘查；Hisatoshi等［7］成功應用直升機地空電磁系統(tǒng)探測到了800m深處的地下水資源.我國地空電磁探測研究起步較晚，吉林大學的研究人員于2009年開始研究地空電磁探測方法，由于直升機勘探費用高、飛行員及飛機在飛行過程中存在風險，所以，研究小組研發(fā)了基于無人飛艇的地空電磁勘探系統(tǒng)［8－9］.基于飛艇的地空電磁勘探方法不僅可以解決直升機航空飛行勘探困難的問題，而且適用于地形復雜的山區(qū)資源探測，在我國深部礦產(chǎn)資源探測和地質(zhì)普查應用中具有廣闊的發(fā)展前景.

但是，在基于飛艇的地空電磁勘探中，飛艇會受風向、大氣氣流、地形、地面局部溫度場變化等影響，導致飛行高度、航跡、姿態(tài)等發(fā)生變化.這些變化使得固定在艇囊前端的接收線圈切割大地磁場，引起磁通量改變而產(chǎn)生感應電動勢.這部分感應電動勢疊加在大地的電磁感應電動勢曲線上，會令測得的電磁信號發(fā)生嚴重的基線漂移，降低電磁數(shù)據(jù)的信噪比，影響電磁數(shù)據(jù)的電阻率成像質(zhì)量.線圈運動引起的噪聲具有幅值較大、分布范圍廣、頻率低的特點，是地空電磁信號中的主要噪聲之一.此外，各種天然磁場及人文設施也會產(chǎn)生較大的噪聲與干擾，如天電噪聲、50Hz的工頻干擾等，這些噪聲會致使數(shù)據(jù)反演解釋時出現(xiàn)假異常，干擾地下電性結構的正確識別，進而影響數(shù)據(jù)解釋的準確度與精度.

目前時間域地空電磁探測方法仍處于起步階段，直接針對地空電磁信號去噪的相關研究較少.不過，近年來，隨著計算機以及信號處理技術的飛速發(fā)展，信號去噪新理論、新方法不斷涌現(xiàn)，并且國內(nèi)外學者在多種航空電磁數(shù)據(jù)噪聲去除方法的研究方面取得了很大進展［10－14］.Fugro航空地球物理探測公司則在硬件上增加線圈運動補償，改善懸吊穩(wěn)定性，并在數(shù)據(jù)處理中采用數(shù)據(jù)疊加以及高通濾波器等方法對運動噪聲進行去除.但是，當采用以上方法時，需要增加分量測量或姿態(tài)記錄硬件電路；Lemire等［12］采用樣條插值和拉格朗日優(yōu)化方法，對時間域航空電磁數(shù)據(jù)中的基線漂移進行了校正，但計算過于繁瑣；對于天電噪聲與白噪聲，Abderrezak Bouchedda等［13］采用平穩(wěn)小波變換對航空電磁數(shù)據(jù)中的天電噪聲進行了抑制；Reninger等［14］利用奇異值分解方法，對時間域航空電磁數(shù)據(jù)中混入的天電噪聲及白噪聲進行了處理.

本文針對基于飛艇的時間域地空電磁數(shù)據(jù)中的主要噪聲，借鑒時間域航空電磁數(shù)據(jù)的消噪方法，在分析了地空電磁信號特點與主要噪聲時頻特性的基礎上，提出一種綜合小波消噪法：首先使用小波高尺度近似分量估計基線漂移，極大限度地消除基線漂移的干擾，再使用小波閾值收縮法有效地消除其余噪聲的影響.通過異常環(huán)模型的理論響應和實測數(shù)據(jù)實驗，結果表明該方法對多種噪聲均有很好的抑制作用，是一種實用有效的時間域地空電磁數(shù)據(jù)消噪方法.

2 綜合小波消噪法

無論飛行載體是直升機還是無人機或飛艇，在空中進行飛行測量時，測量的電磁信號中均會混入多種噪聲，主要包括：飛行時搖擺導致接收線圈產(chǎn)生的低頻類基線漂移、分布范圍廣的隨機噪聲、天電噪聲及工頻干擾等.白噪聲、天電噪聲及工頻干擾可以通過小波閾值收縮各尺度上的小波系數(shù)進行抑制，但基線漂移的能量主要集中在高尺度的近似分量上，所以使用閾值收縮法往往不能很好地去除基線.為了彌補這一不足，根據(jù)基線漂移類噪聲特點，有兩種改善方案，一是提高閾值法中的分解級數(shù)，二是使用兩種算法相結合的思路.

增加閾值收縮法中的小波分解尺度，使用最高尺度的近似分量估計基線漂移，然后收縮各級的小波系數(shù)，重構時使用處理過的小波系數(shù)，并將最高尺度的近似分量置0，這種算法可以達到既消除基線漂移又消除其他噪聲的目的.但是，隨分解尺度增多，閾值處理尺度也會增多，有用信息損失的風險就會增加，重構時誤差也會增大，甚至消噪后的信噪比不升反降.并且，尺度增多計算量會增大，處理速度也會降低.所以，通過理論與實踐證明，高尺度分解的閾值收縮法并不適合地空信號的消噪.

因此，本文提出使用兩種方法相結合的綜合消噪法：首先使用一種方法極大限度地消除基線漂移的干擾，再使用另一種方法有效地消除其余噪聲的影響.

2.1 小波基的選取

不同的小波基，具有不同的時頻特性，所以使用不同的小波基進行地空電磁信號的消噪，結果可能相差甚遠.

根據(jù)Mallat的最優(yōu)小波基選取準則：只需用較少數(shù)目的非零小波系數(shù)即可精確地重構信號［15］.我們使用常用小波基對地空電磁信號做了多次分解重構實驗，并結合它們的數(shù)學特性以及要處理的地空電磁信號特點，最終選取sym8作為綜合小波消噪算法中的小波基.

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

基線漂移的能量主要集中在高尺度的近似分量上，幾乎沒有奇異點存在，而閾值收縮法的分解尺度不能過大，并且主要是對小波系數(shù)進行處理，不對近似系數(shù)進行處理，所以不適合校正信號中的基線漂移.

利用小波的多分辨率分析特性可以將信號逐級分解成低頻信息與高頻信息，將交織在一起的各種不同頻率組成的混合信號分解成不同頻帶的子信號.低頻信息對應著信號的近似部分，即尺度系數(shù)；高頻信息對應著細節(jié)部分，即小波系數(shù)［11－13］.本文地空數(shù)據(jù)的采樣率為51.2kHz，根據(jù)采樣定理，信號中包含的頻率信息為0～25.6KHz.由于基線漂移頻率較低，根據(jù)本文中每道采集的數(shù)據(jù)點數(shù)，確定分解級數(shù)為12級，即利用sym8小波12級分解所得的近似分量作為基線漂移的估計，從原始信號中減去這部分低頻信息，所得信號即為校正信號.

2.3 基于小波收縮法的地空信號去噪算法

小波收縮法對于頻率范圍分布廣的白噪聲，以及與有效信號發(fā)生頻率重疊的噪聲具有很好的去除效果.將實際信號X寫成向量形式有

S為有效信號，e為噪聲，目的是從X中去除噪聲e，估計出S.若e為滿足獨立同分布（i.i.d.）或 N（0，σ2）的噪聲，則可以使用正交離散小波變換矩陣W來估計信號S：

其中w，θ，與z分別是X，S與e的正交小波變換.由于e已被假設為白噪聲，所以變換z也為白噪聲，此時，小波系數(shù)可表示為

我們的目標是估計θ具有最小的L2風險，最小均方誤差MSE，即，找到最佳的估計具有最小L2風險：

（1）分解層數(shù)的選取

閾值收縮法中分解層數(shù)的選取直接影響到消噪效果.隨著小波分解層數(shù)的增多，得到的細節(jié)信息就越豐富.但是，并不是層數(shù)越多越好，因為小波變換中分解合成實際上是對頻帶的劃分，分解層數(shù)越多，產(chǎn)生的子帶越多，頻帶劃分得越細，由于上一級頻帶分解的信號輸出又作為下一級頻帶分解的輸入，層數(shù)增加意味著級間的濾波器增多，造成信號移位也越大；另一方面小波分解合成都要進行邊界延拓，層數(shù)越多引起邊界失真越大.而這些損失的信息量均是小波逆變換不能恢復的損失，因此基于小波分解的層數(shù)不宜過高，根據(jù)地空有效信號的頻帶，本文選用5層小波分解.

（2）閾值與閾值收縮策略

根據(jù)電磁信號特點，對比了多種閾值策略和閾值收縮函數(shù)后，本文應用極小極大（Minimaxi）原理對每一尺度上的小波系數(shù)計算相應的Minimaxi閾值.Minimaxi原理是一種找出失敗的最大可能性中最小值的算法.Minimaxi閾值相比其他閾值略小，采用這個閾值策略主要考慮讓信號盡可能小的減少有效信息的損失.閾值收縮函數(shù)也各有利弊，由于電磁信號具有衰減特性，為此我們采用硬閾值收縮函數(shù)（Hard），Hard收縮函數(shù)最適于保持信號的奇異點.

3 理論模型算法測試

由于閉合異常環(huán)模型不僅可以模擬有限導體的晚期電磁響應，而且在野外便于鋪設.所以本文首先基于閉合異常環(huán)模型計算的理論響應，對算法壓制噪聲和校正基線的有效性進行驗證.

采用閉合異常環(huán)模型，進行電性源地空電磁探測實驗的示意圖如圖1所示.

閉合異常環(huán)用一個等效電感L和等效電阻R進行近似，當發(fā)射線圈中通有交變電流時，在閉合異常環(huán)中產(chǎn)生的感應電流可以寫為

式中iAL為閉合異常環(huán)的感應電流，LAL為線圈的電感，MTL為發(fā)射導線與閉合異常環(huán)的互感.τL為閉合異常環(huán)的時間常數(shù)，與線圈的材料、電氣性能有關，其表達式為

式中μ0為磁導率，n為閉合異常環(huán)匝數(shù)，rAL為線圈半徑，ra為導線半徑，ρL為導線電阻率.當發(fā)射電流為梯形波時，在閉合異常環(huán)的取樣單元中產(chǎn)生的感應電壓信號為

圖1 基于異常環(huán)模型的電性源地空電磁探測示意圖Fig.1 Schematic diagram of electrical－source groundairborne electromagnetic detection system based on anomaly loop model

式中MRL為接收線圈與閉合異常環(huán)的互感，I為發(fā)射電流幅值，T1為電流上升時間，T2為電流下降時間，T3為發(fā)射電流寬度.

對于任意位置、邊長為2a×2b的矩形異常環(huán)與半徑為r2的圓形接收線圈之間的互感系數(shù)，可以寫為

式中（Tx，Ty，Tz）為接收線圈的圓心坐標.

對于長導線與異常環(huán)的互感可以通過簡化公式（9）得到.基于吉林大學自主研制的電磁探測系統(tǒng)，工作參數(shù)如表1所示，采用高斯數(shù)值積分方法，利用公式（6）和（8），就可以計算出閉合異常環(huán)的電磁響應［17－18］.

表1 異常環(huán)地空實驗工作參數(shù)表Table 1 Experimental work parameters of ground－airborne electromagnetic detection based on anomaly loop model

將計算的異常環(huán)電磁響應作為有效信號s（n），如圖2a所示.采用信噪比為3dB的Gaussian白噪聲與低頻信號來模擬電磁噪聲和運動噪聲e（n），在s（n）上疊加e（n）來模擬含噪信號x（n），由于線圈運動引起的噪聲呈現(xiàn)低頻特性，為了清晰看到基線漂移，我們同時顯示了多個周期的信號，如圖2b所示.

圖2 電性源地空測量時異常環(huán)的電磁響應（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）或誤差范數(shù)（Error norm）來量化評估消噪方法的有效性.MSE定義為

MSE的值越小，則兩信號的相似度就越高；Error norm的值越小，表明兩信號的誤差越小.表2給出了模擬含噪信號與綜合法消噪信號的均方誤差值與誤差范數(shù).可以看到，綜合法消噪后有效降低了MSE值與Error norm值，消噪后的信號與理論信號的相似度高、誤差小.

表2 信號消噪前后MSE與Error norm比較Table 2 Comparison of MSE and Error norm of noisy signal and de－noised signal

4 野外測數(shù)據(jù)消噪驗證

2010年11月，在長春市農(nóng)安縣燒鍋鎮(zhèn)化東屯，采用吉林大學自主研制的電性源發(fā)射系統(tǒng)進行地面發(fā)射，利用吊車將地空接收系統(tǒng)和接收線圈吊在空中5m高，開展了低阻異常模型的飛行模擬實驗.具體工作參數(shù)同表1，接收系統(tǒng)采用24位AD采集、100KHz采樣頻率、全波形連續(xù)采集、GPS同步控制以及無線遙控方式.實驗中，我們首先測量了背景場，然后鋪設了異常環(huán)（匝數(shù)為10m×10m×5匝、時間常數(shù)為4.62ms），位置如圖4所示，最后沿發(fā)射電性源方向測量了4次，分別如圖4中R1－R4測線所示.

圖5a為異常環(huán)模型的野外實測數(shù)據(jù)，圖5b為使用綜合小波法去噪后的信號，通過對比兩信號，可以看到，該方法在失真度很小的前提下，很好地校正了基線.

圖6為R3測線未減去背景場的原始數(shù)據(jù)剖面曲線，圖7為減去背景場并進行綜合小波消噪后的剖面曲線，在測線的0～55m、78～120m均為背景場區(qū)域，由于沒有異常存在，所以各道數(shù)據(jù)基本重合；但在測線的56～77m處，對應異常環(huán)位置，剖面曲線變化明顯.圖8為理論計算的異常環(huán)電磁響應剖面曲線，可以看到，圖7所示的剖面曲線與圖8中理論計算的剖面曲線變化規(guī)律一致.但是，從圖7和圖8中可以觀察到異常寬度存在一定的差異，主要原因是無法記錄實驗采用吊車在測量時的轉(zhuǎn)動速度、吊車轉(zhuǎn)動通過異常環(huán)的時間，只能估算測點的距離，因此導致了實測數(shù)據(jù)中的異常寬度與理論計算的異常寬度存在誤差，此誤差可以在計算視電阻率斷面圖時進行修正.

圖4 加異常線圈測試場地工區(qū)圖Fig.4 Schematic diagram of wild field laid an anomaly loop

在對實測數(shù)據(jù)進行視電阻率計算時，根據(jù)場源等效思想，將長導線源剖分為若干個電偶極源，基于均勻半空間模型求解了地面LOTEM的全區(qū)視電阻率，然后進行了視電阻率－深度計算［9］.圖9為采用原始數(shù)據(jù)計算的視電阻率斷面圖，圖10為采用綜合小波消噪法后計算的視電阻率斷面圖.燒鍋鎮(zhèn)化東屯地質(zhì)概況屬于長春地區(qū)第四系結構，地表為黃土覆蓋層，基巖為白堊系地層，以砂礫巖和泥巖為主，沙層、泥巖和泥沙巖等交替出現(xiàn)，膠結微弱孔隙較發(fā)育.距地表13m為黃土狀亞粘土低阻蓋層，地下50m至70m之間為含水砂巖，電阻率值在20～50Ωm之間變化，電阻率變化與鉆井資料一致.從圖9中可以看到由于基線漂移與噪聲的存在，導致解釋結果中異常傾斜，而且即使減掉大地的電磁響應，也出現(xiàn)了虛假的低阻異常；而圖10中由異常環(huán)引起的低阻異常清晰，且在淺層和深層異常一致.

5 結 論

5.1 本文根據(jù)地空電磁信號和主要噪聲的特點，通過分析常用小波系的數(shù)學特性與比較信號重構誤差，提出了一種基于sym8小波基的綜合小波去噪法對地空電磁信號中混入的多種噪聲進行消除.針對基線漂移類低頻噪聲，先基于小波多分辨率分析，將原始信號進行多層小波分解，通過計算最高分解層的近似分量來模擬電磁數(shù)據(jù)中的基線漂移，再從原始數(shù)據(jù)中去除該近似分量以達到消除基線漂移的目的.針對頻率分布廣的白噪聲，應用MiniMaxi閾值配合硬收縮函數(shù)進行5級小波收縮進行消噪，既保證了數(shù)據(jù)波形的失真度小，又較好地抑制了具有白噪聲特征的多種電磁噪聲.

5.2 基于異常環(huán)模型進行理論電磁響應仿真，比較了綜合小波去噪法消噪前后信號的MSE與波形失真度.并將算法用于了長春市燒鍋鎮(zhèn)野外異常環(huán)模型實測數(shù)據(jù)的處理，比較了消噪前后疊加所成的視電阻率斷面圖.結果表明這種綜合消噪法對多種噪聲均有很好的抑制作用，是一種實用有效的時間域地空電磁數(shù)據(jù)消噪方法.

（References）

［1］Fountain D.60years of airborne em－focus on the last decade.AEM2008－5th International Conference on Airborne Electromagnetics Haikko Manor，F(xiàn)inland，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信號提取.煤田地質(zhì)與勘探，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系統(tǒng)測試和標定方法研究.地球物理學報，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.