時(shí)間域航空電磁的天電噪聲去除研究

賁放，黃威，路寧，韓飛，鄭紅閃，丁志強(qiáng)，李軍峰

(1.自然資源部 地球物理電磁法探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，河北 廊坊 065000； 2.中國(guó)地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所，河北 廊坊 065000； 3.國(guó)家現(xiàn)代地質(zhì)勘查工程技術(shù)研究中心，河北 廊坊 065000； 4.中水東北勘測(cè)設(shè)計(jì)研究有限責(zé)任公司，吉林 長(zhǎng)春 130026)

0 引言

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

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

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

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

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

1 天電噪聲去除原理

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

其中：Y是經(jīng)過(guò)剪切均值濾波后一個(gè)窗寬的數(shù)據(jù)；N是濾波的窗寬(奇數(shù))；Lw是單邊剪切窗寬，其中需要滿足2Lw

1) 從原始數(shù)據(jù)中選取對(duì)應(yīng)濾波窗寬(N)長(zhǎng)度的數(shù)據(jù)；

2) 對(duì)窗口內(nèi)數(shù)據(jù)進(jìn)行升降排序；

3) 對(duì)排序數(shù)列的兩端數(shù)據(jù)進(jìn)行剪切，單邊剪切長(zhǎng)度為L(zhǎng)w；

4) 將剩余數(shù)據(jù)取均值，作為最后剪切均值濾波的結(jié)果；

5) 滑動(dòng)窗口對(duì)所有待處理數(shù)據(jù)進(jìn)行剪切均值濾波。

以上給出了剪切均值濾波的實(shí)現(xiàn)過(guò)程。當(dāng)對(duì)測(cè)線數(shù)據(jù)進(jìn)行濾波處理時(shí)，需要考慮邊界數(shù)據(jù)的處理方法。本文考慮幾種常用的處理方法，如周期擴(kuò)邊(圖2a)、鏡像擴(kuò)邊(圖2b)，兩種方法均可以較好地對(duì)邊界數(shù)據(jù)進(jìn)行處理，保持?jǐn)?shù)據(jù)的邊界變化趨勢(shì)。

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

圖2 兩種擴(kuò)邊方法Fig.2 Two extension methods

2 天電去除研究

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

2.1 濾波窗寬

當(dāng)濾波窗口長(zhǎng)度N取值過(guò)小時(shí)，會(huì)出現(xiàn)某一窗口內(nèi)均為天電噪聲數(shù)據(jù)，導(dǎo)致處理結(jié)果中天電噪聲去除不徹底；當(dāng)濾波窗口長(zhǎng)度N取值過(guò)大時(shí)，不僅會(huì)增加計(jì)算時(shí)間，且處理結(jié)果中供電和斷電區(qū)間內(nèi)的數(shù)據(jù)過(guò)于平均，導(dǎo)致供電窗口數(shù)據(jù)提前和斷電窗口數(shù)據(jù)延后的特征，甚至將斷電處有用信號(hào)濾平?，F(xiàn)截取部分含有天電噪聲的數(shù)據(jù)，參數(shù)中單邊剪切長(zhǎng)度設(shè)置為8，改變?yōu)V波窗口大小來(lái)說(shuō)明去噪效果。

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

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

2.2 剪切長(zhǎng)度

剪切長(zhǎng)度的大小即為窗口兩端去除數(shù)據(jù)量的大小，它決定了窗口內(nèi)噪聲數(shù)據(jù)去除的比例。剪切長(zhǎng)度過(guò)小或過(guò)大，均可能導(dǎo)致窗口被平均化，或供電處、斷電處濾波后數(shù)據(jù)偏移正確的方向。

圖4中分別給出了不同單邊剪切長(zhǎng)度的去噪效果，其中濾波窗口N=101，單邊剪切長(zhǎng)度Lw=8、30、45。經(jīng)過(guò)濾波后的供電處、峰值、斷電處，Lw=8、30均出現(xiàn)了數(shù)據(jù)被嚴(yán)重平均化，導(dǎo)致了供電窗口提前、峰值變小、斷電窗口延遲等問(wèn)題(圖4中b、c處深藍(lán)和天藍(lán)實(shí)線)，Lw=8時(shí)以上現(xiàn)象較嚴(yán)重。當(dāng)Lw=45時(shí)，供電和斷電處濾波效果相對(duì)較好，但峰值出現(xiàn)了過(guò)于平均化現(xiàn)象(圖4中b處紫色實(shí)線)。圖4中斷電數(shù)據(jù)含有天電噪聲，由于天電噪聲位于斷電數(shù)據(jù)中部且為大窗口剪切均值濾波，因此剪切長(zhǎng)度為8和30時(shí)均可以較好地去除噪聲，當(dāng)剪切長(zhǎng)度取值過(guò)大(Lw=45)時(shí)，濾波后數(shù)據(jù)波動(dòng)明顯(圖4中a處紫色實(shí)線)?？梢钥闯鰹V波窗口長(zhǎng)度和剪切長(zhǎng)度要相互參考適當(dāng)取值。

2.3 混合窗口天電噪聲濾波

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

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

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

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

3 結(jié)論

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

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