賁放，黃威，路寧，韓飛，鄭紅閃，丁志強，李軍峰

(1.自然資源部 地球物理電磁法探測技術重點實驗室，河北 廊坊 065000； 2.中國地質科學院 地球物理地球化學勘查研究所，河北 廊坊 065000； 3.國家現代地質勘查工程技術研究中心，河北 廊坊 065000； 4.中水東北勘測設計研究有限責任公司，吉林 長春 130026)

0 引言

隨著航空電磁正反演技術及儀器采集技術的進步，航空電磁法已廣泛應用于諸多領域[1]，實現了數據密集多通道采集，采集的原始數據需要進行數據預處理為數據解釋做準備。原始數據中包含了多種噪聲，其中對數據影響較大的為天電噪聲和運動噪聲。降低天電噪聲影響可以提高數據信噪比，并有利于增大勘探深度[2-4]和提高小異常的分辨能力[5]。

天電主要由大氣中的帶電電荷經過磁暴或者雷電作用發(fā)生電磁輻射而產生，并通過空氣進行傳播。磁暴引起的天電噪聲幅值大且不可預測，而雷電引起的天電噪聲有一定的規(guī)律。全球雷電發(fā)生的頻率約在100次/秒，在電磁信號中表現為瞬時噪聲(尖脈沖或振蕩)，測量數據中表現為突跳現象。天電干擾的數量、頻率、強度、持續(xù)時間及幅度隨時間、位置和源距的不同均有所變化，對各頻段數據的影響依賴于電離層情況。天電噪聲的頻帶范圍為5～100 kHz，在低頻范圍內幅值的變化和頻率近似成反比，高頻部分幅值迅速增長[6]，在500 Hz～2.5 kHz的范圍內對數據影響較小，而在低于500 Hz和2.5 kHz～10 kHz的范圍對數據影響較大[5]。赤道或熱帶地區(qū)的天電干擾比中緯度地區(qū)高200多倍，而中緯度地區(qū)的天電干擾夏天比冬天高100多倍；中午的天電干擾比早晚高10多倍[7]。在較近的勘探區(qū)，由于天電直接傳播到測點導致數據受天電干擾較大，使得航空電磁測量能力下降甚至無法分辨；在較遠的勘探區(qū)，天電通過電磁波在地面—電離層波導中多次反射到測點，相比較而言對數據影響小。經過反射的天電由多個閃電疊加而成[8]。天電噪聲的水平極化導致其對數據垂直分量影響非常小，但由于作業(yè)過程中線圈會發(fā)生姿態(tài)變化，因此數據垂直磁場分量同樣也受到天電噪聲的影響[5]。在背景噪聲較小地區(qū)，水平分量和垂直分量數據中分別有0.05%和 0.005%的噪聲[9]，在天電噪聲頻繁情況下觀測的數據，水平分量數據有2%含有天電噪聲，垂直分量則有0.5%[4]含有天電噪聲。

通過對天電特性的分析，研究去除天電噪聲的方法。早在1984年，Macnae等利用裁剪法(pruning)或錐形疊加的方法對非平穩(wěn)或相干噪聲進行去除，達到提高信噪比的目的，其認為高頻率采樣的航空電磁數據量充足，且經過研究含有天電的數據量很小，裁剪掉含有天電影響的連續(xù)半周期數據不會對整體數據質量造成影響[9]。同年，Bednar和Watt提出利用α-trimmed均值/中值濾波的方法對噪聲進行濾波，該方法可以對大規(guī)模數據進行快速濾波，對天電等突跳特點的噪聲具有較好的濾除效果[10]。1989年，Sutarno和Vozoff[11]提出了回歸M值—估計法，同時對最小平方法的最佳給出了證明。隨后，Buselli和Cameron[4]利用該方法對天電噪聲進行去除，同時對中值濾波、魯棒統(tǒng)計法、均值濾波、M值—估計法進行了比較，經過研究當天電噪聲符合高斯分布特性時，M值—估計法去除天電噪聲較其他方法更好。1998年，Buselli等通過實測數據中天電和地面基站觀測信號的相關性對天電噪聲進行識別，闡述了局部噪聲預測濾波和遠噪聲濾波，利用預測濾波對高頻部分的天電噪聲進行去除，去除過程中應用了神經網格。Cull[12]提出利用線性疊加對天電噪聲進行去除，Lane等[13]通過對天電噪聲的來源及特征進行分析后，選擇使用非線性濾波對天電噪聲進行去除。2000年，Leggatt等通過增大發(fā)射功率和利用更大計算能力的計算機來阻擋天電噪聲對數據的干擾。2010年，Bouchedda等[5]利用小波提取和小波變換對天電噪聲進行去除，小波提取技術是基于平穩(wěn)小波變換，通過將時間域航空電磁數據分解為不同頻率成分的細節(jié)和近似量[14]，在一階小波細節(jié)系數中利用能量檢測方法識別天電噪聲，將數據中與天電噪聲相關的小波細節(jié)系數進行去除，通過剩余小波細節(jié)系數利用反變換重構方法對信號進行重構，從而實現去除天電噪聲的目的。雖然該方法可以有效地對天電噪聲進行去除，但仍存在一定的不足，當天電噪聲位于發(fā)射電流開、關期間內時，由于天電噪聲中低頻成分占主導降低了該方法對天電去除的有效性。為解決該問題，Bouchedda等提出了小波變換的方法，利用靜態(tài)小波變換具有線性和位移不變性的特征，在小波域中進行小波細節(jié)系數疊加，其中利用了算術平均、中值平均及α-trimmed均值疊加技術，通過對合成數據和實測數據分析說明了該方法非常適合去除晚期道低頻天電干擾。Reninger等[15]利用奇異值分解法對天電噪聲進行有效去除。

對于天電噪聲，去除的方法眾多，各有優(yōu)劣。裁剪法固然有效，但需要對天電噪聲進行識別后方可剪裁，當天電噪聲較多時，直接剪掉數據會影響疊加次數，從而對后期的數據預處理產生影響。直接對含有天電噪聲數據進行中值濾波或均值濾波會使數據曲線中心產生偏移，例如當窗口滑動到含有天電成分的數據時，因數據窗內部含有尖峰信號，濾波后數值會遠大于準確值，從而產生了因數據預處理而出現的假異常，且均值濾波是對信號進行平滑處理，對噪聲去除效果不明顯[16]。M值—估計法的疊加窗較短不適合天電噪聲較多的情況，而實測數據中天電噪聲含量是無法預測和估計的。利用神經網絡對天電噪聲進行去除時，需要利用不同時間、不同位置的大量數據對神經網絡進行訓練，嚴重降低了其天電噪聲去除的速度[16]。小波變換和奇異值分解法去除天電噪聲耗時較多。在實測數據中去除天電噪聲，過程中不僅需要考慮去除效果，同時需要考慮處理速度等因素，而Bednar和Watt提出的α-trimmed均值/中值濾波方法不僅可以快速去除天電噪聲，且去除效果較好。

綜上，為保證原始數據量的不變性，參考α-trimmed均值/中值濾波方法進行改進，可稱之為剪切均值濾波。原方法中描述不同的α(0～0.5)取值，對應了均值濾波或中值濾波，其中對[αN]進行判斷和取整，α、N的取值對數據處理效果有較大影響，而剪切均值濾波直接給定濾波的窗寬(整數)、單邊剪切長度(整數)，無需判斷[αN]的情況進行取整，并通過對單窗口濾波參數(窗口長度，單邊剪切窗寬)進行研究，提出對實測數據進行混合窗口濾波。

1 天電噪聲去除原理

電磁工作的主要頻率范圍在5～25 kHz[8]，在實際數據采集過程中經常受到各種不同噪聲的干擾，由于噪聲的來源廣泛，增加了實測數據預處理的難度，且不同于地面瞬變電磁，航空電磁[17]數據采集過程中多了天電噪聲和運動噪聲，是對數據影響較大的兩種噪聲?；讦?trimmed均值/中值濾波方法[10]改進的剪切均值濾波，取消了對[αN]的判斷求整，直接給出濾波的窗寬、單邊剪切窗寬:

其中：Y是經過剪切均值濾波后一個窗寬的數據；N是濾波的窗寬(奇數)；Lw是單邊剪切窗寬，其中需要滿足2Lw

1) 從原始數據中選取對應濾波窗寬(N)長度的數據；

2) 對窗口內數據進行升降排序；

3) 對排序數列的兩端數據進行剪切，單邊剪切長度為Lw；

4) 將剩余數據取均值，作為最后剪切均值濾波的結果；

5) 滑動窗口對所有待處理數據進行剪切均值濾波。

以上給出了剪切均值濾波的實現過程。當對測線數據進行濾波處理時，需要考慮邊界數據的處理方法。本文考慮幾種常用的處理方法，如周期擴邊(圖2a)、鏡像擴邊(圖2b)，兩種方法均可以較好地對邊界數據進行處理，保持數據的邊界變化趨勢。

圖1 剪切均值濾波過程示意Fig.1 The diagram of Trim-mean filter processing

圖2 兩種擴邊方法Fig.2 Two extension methods

2 天電去除研究

剪切均值濾波涉及到兩個參數：濾波的窗寬、單邊剪切窗寬，兩個參數的取值直接影響濾波效果。下面對兩個參數的選取進行研究，文中實測數據是由基頻25 Hz的時間域航空電磁系統(tǒng)采集，采樣率為100 kHz，半正弦發(fā)射波形的占空比為1∶4。

2.1 濾波窗寬

當濾波窗口長度N取值過小時，會出現某一窗口內均為天電噪聲數據，導致處理結果中天電噪聲去除不徹底；當濾波窗口長度N取值過大時，不僅會增加計算時間，且處理結果中供電和斷電區(qū)間內的數據過于平均，導致供電窗口數據提前和斷電窗口數據延后的特征，甚至將斷電處有用信號濾平?，F截取部分含有天電噪聲的數據，參數中單邊剪切長度設置為8，改變?yōu)V波窗口大小來說明去噪效果。

圖3中分別給出了不同濾波窗口長度(N=21、41、101)的去噪效果。當N值較小(N=21、41)時，去噪結果雖然保持待處理數據的變化趨勢(圖3中b、c處)，但天電噪聲濾除不徹底(圖3中a處藍色實線)；當濾波窗口長度較大(N=101)時，可以較好地濾除天電噪聲(圖3中a處紫色實線)，但會使得供電處、斷電處及峰值數據偏離原數據方向，出現數據過于平均而丟失有用信號(圖3中c處紫色實線)或者供電數據提前(圖3中b處紫色實線)、斷電數據延后(圖3中c處紫色實線)、峰值減小的情況。從本組處理結果可以看出小窗口濾波可以很好地保持供電處、峰值、斷電處的數據變化情況，大窗口濾波可以較好地去除數據中的天電噪聲。

圖3 不同濾波窗口長度對噪聲去除效果比較Fig.3 Comparison of noise removal effects of different filter window lengths

2.2 剪切長度

剪切長度的大小即為窗口兩端去除數據量的大小，它決定了窗口內噪聲數據去除的比例。剪切長度過小或過大，均可能導致窗口被平均化，或供電處、斷電處濾波后數據偏移正確的方向。

圖4中分別給出了不同單邊剪切長度的去噪效果，其中濾波窗口N=101，單邊剪切長度Lw=8、30、45。經過濾波后的供電處、峰值、斷電處，Lw=8、30均出現了數據被嚴重平均化，導致了供電窗口提前、峰值變小、斷電窗口延遲等問題(圖4中b、c處深藍和天藍實線)，Lw=8時以上現象較嚴重。當Lw=45時，供電和斷電處濾波效果相對較好，但峰值出現了過于平均化現象(圖4中b處紫色實線)。圖4中斷電數據含有天電噪聲，由于天電噪聲位于斷電數據中部且為大窗口剪切均值濾波，因此剪切長度為8和30時均可以較好地去除噪聲，當剪切長度取值過大(Lw=45)時，濾波后數據波動明顯(圖4中a處紫色實線)。可以看出濾波窗口長度和剪切長度要相互參考適當取值。

2.3 混合窗口天電噪聲濾波

天電噪聲是不可預見的，可考慮利用閾值法先識別天電噪聲再進行去除，具體可參考文獻[8]、[9]和[16]。根據前面討論的情況可知，剪切均值濾波法的小窗口濾波適合數據趨勢變化明顯部分，如供電處、峰值和斷電處，大窗口濾波適合斷電后數據的處理，因此本文提出采用混合窗口濾波來對天電噪聲進行去除，即供電處、峰值和斷電處數據采用小窗口剪切均值濾波，斷電后衰減數據使用大窗口剪切均值濾波。

圖4 不同單邊剪切長度對噪聲去除效果比較Fig.4 Comparison of noise removal effects of different trim lengths

圖5給出了利用混合窗口對數據進行天電去除的效果，其中小窗口濾波參數為：N=11、Lw=3，大窗口濾波參數為N=101、Lw=30。從圖中可以明顯看出，小窗口濾波后數據保持了原數據供電處、峰值和斷電處的變化趨勢(圖5中a處)；大窗口濾波可以較好的濾除天電噪聲(圖5中b處供電階段、c處斷電階段)。采用混合窗口剪切均值濾波即可以保持數據變化趨勢，又可以達到有效去除天電噪聲的目的。

圖5 混合窗口剪切均值濾波去噪效果Fig.5 Effects of hybrid window trim-mean filter

3 結論

本文利用剪切均值濾波對時間域航空電磁數據進行天電噪聲去除。航空電磁測線數據中包含了快速變化的供電數據及斷電后變化相對緩慢的數據，在研究過程中發(fā)現，如對測線數據進行同時間道剪切均值濾波，則需要先去除運動噪聲，否則會出現數據整體趨勢的均值化。在對測線數據進行剪切均值濾波時，濾波效果受濾波參數影響。經過研究，小窗口參數濾波會造成數據中天電噪聲去除不徹底，大窗口參數濾波會造成數據過度平均化，損失原始信號嚴重。為此，本文利用混合窗口剪切均值濾波對天電噪聲進行有效去除。混合窗口方法的使用不僅保留了供電處、峰值和斷電處有用信號的變化趨勢和穩(wěn)定性，且有效去除了天電噪聲，是一種快速、高效去除天電噪聲的方法。天電噪聲對水平分量影響大，因此希望通過本文研究，不僅可以有效地去除天電噪聲提高數據信噪比和增大勘探深度，且可以實現數據水平和垂直分量的聯合應用。

致謝：向中科院上海微系統(tǒng)所助理研究員裴易峰和吉林大學電磁“千人計劃”研究團隊成員提供的幫助表示感謝。