基于Robust-M估計的地空頻率域電磁探測數(shù)據(jù)干擾抑制方法

劉長勝， 白江蘭， 康寧， 李茁維， 梁潔， 周海根*， 闞紹佑

(1.吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院， 長春 130012； 2.地球信息探測儀器教育部重點實驗室， 長春 130012)

20世紀70年代初，國外出現(xiàn)了基于地面發(fā)射、空中接收的電磁探測系統(tǒng)TURAIR，有較好的分辨率和電阻率鑒別率，應(yīng)用于探測深度大或地面樹林高、地形起伏大等地區(qū)[1-2]。地空電磁探測將地面電磁探測與航空電磁探測相結(jié)合，早期的地空電磁探測使用直升機搭載探測，成本較高[3]；中國較早進行地空電磁探測時，采用無人飛艇搭載電磁接收系統(tǒng)，主要應(yīng)用于沙漠、海陸交互帶以及山區(qū)等復(fù)雜地形區(qū)域，具有探測深度大、分辨率高等優(yōu)點[4-5]。目前中國無人機行業(yè)快速發(fā)展，同時也帶動了地空電磁探測的發(fā)展，使用無人機搭載接收系統(tǒng)進行探測，降低了飛行成本，在空中進行飛行測量，相比較在地面測量，節(jié)省了大量的人力物力[6]。

地空電磁探測分為地空時間域電磁探測[7]和地空頻率域電磁探測，各有所長，為了綜合時域和頻域電磁探測方法，簡化探測過程，通過時頻融合傳輸方法，同時激發(fā)大地產(chǎn)生瞬態(tài)和穩(wěn)態(tài)響應(yīng)，從而提高電磁探測效率[8]。目前中外地空頻率域電磁探測主要是對垂直磁場的測量，為了提高探測的準確性，也在不斷向著三分量磁場采集的方向發(fā)展[9]，同時進行姿態(tài)數(shù)據(jù)的采集，對磁場進行姿態(tài)校正[10]。地空頻率域電磁探測，通過接地導(dǎo)線向大地發(fā)射2n序列偽隨機波形，在空中采用線圈進行磁場數(shù)據(jù)采集，在不降低效率的情況下，還能保持一定的分辨率[11]。相比于地面電磁探測，其探測效率高，可快速獲得大范圍內(nèi)的地電信息，在復(fù)雜地形環(huán)境中適用性強[12]。地空頻率域電磁探測系統(tǒng)的探測范圍和深度與發(fā)射頻率、發(fā)射功率、實際環(huán)境噪聲、系統(tǒng)參數(shù)有關(guān)[13]。在復(fù)雜地形區(qū)域探測時，單源發(fā)射導(dǎo)致能量分散，接收端信號微弱，信噪比低，為了進行快速準確探測，Zhou等[14-15]提出了一種雙源發(fā)射的地空頻率域電磁探測方法，提高了探測精度，進而對多源發(fā)射的實現(xiàn)提供了理論基礎(chǔ)。通過多源發(fā)射的實際應(yīng)用，探測結(jié)果顯示接收信號明顯增強，從而使地空電磁探測更加適用于噪聲更強的探測區(qū)域。

空中接收電磁系統(tǒng)因受到風(fēng)向、飛行顛簸等因素的影響，導(dǎo)致數(shù)據(jù)存在低頻運動噪聲、姿態(tài)噪聲、隨機噪聲、尖峰噪聲等，通過小波變換可以有效去除低頻運動噪聲。朱凱光等[16]研究了基于BP(back propagation)神經(jīng)網(wǎng)絡(luò)的線圈姿態(tài)校正算法，提高了數(shù)據(jù)質(zhì)量。穩(wěn)健M估計常用于時間域大地電磁探測[17]，壓制高斯隨機噪聲和尖峰脈沖噪聲[18]。對離群值進行剔除或降權(quán)[19]，實現(xiàn)對噪聲的壓制。

地空頻率域電磁探測，針對不同噪聲采取不同的噪聲抑制方法，具有一定的效果[20-21]。但在數(shù)據(jù)處理中，噪聲壓制仍然存在很多問題，因受到飛行因素與環(huán)境因素的影響，通過無人機搭載線圈進行探測，電磁數(shù)據(jù)并不平穩(wěn)，信噪比出現(xiàn)忽高忽低現(xiàn)象，影響數(shù)據(jù)質(zhì)量。通過引入Robust-M估計，使信號趨于穩(wěn)定，提高數(shù)據(jù)質(zhì)量。

1 理論基礎(chǔ)

1.1 地空頻率域電磁探測原理

地空頻率域電磁探測系統(tǒng)，如圖1所示。

由發(fā)射系統(tǒng)通過接地導(dǎo)線向大地發(fā)射2n序列偽隨機波形，在測區(qū)采用無人機搭載線圈，進行數(shù)據(jù)采集。接收系統(tǒng)采集數(shù)據(jù)為電信號V，通過數(shù)據(jù)預(yù)處理得到垂直磁場B。即

(1)

圖1 地空頻率域電磁探測示意圖Fig.1 GAFED Schematic diagram

U=-jnBSω=-jnBS·2πf

(2)

(3)

式中：φ為磁通量；n為線圈匝數(shù)；S為線圈有效面積；f為發(fā)射頻率。

1.2 傳統(tǒng)時窗選取

傳統(tǒng)的數(shù)據(jù)處理是對數(shù)據(jù)加窗分段后，進行傅里葉變換，進而得到各個頻率下的磁感應(yīng)強度。數(shù)據(jù)處理的時窗寬度，通常選取為1 s、2 s、3 s、5 s、10 s、…，通過數(shù)據(jù)信噪比占比統(tǒng)計結(jié)果，對分段時窗寬度進行調(diào)整。對野外探測的數(shù)據(jù)，使用不同時窗寬度進行處理，當(dāng)時窗為2 s時，其中一條測線的磁感應(yīng)強度幅值曲線，如圖2所示，通過計算測線數(shù)據(jù)信噪比，統(tǒng)計測區(qū)信噪比占比，結(jié)果如圖3所示。測區(qū)具體情況詳見3.2節(jié)。信噪比為

(4)

式(4)中：SNR為信噪比；S為信號幅值；N為噪聲幅值[22]。

因野外測量，采用無人機搭載探測，而受到環(huán)境以及飛行因素的影響，數(shù)據(jù)存在隨機噪聲與尖峰噪聲，導(dǎo)致信噪比忽高忽低，數(shù)據(jù)質(zhì)量差。而時窗寬度與數(shù)據(jù)信噪比有著直接的關(guān)系，由圖3可知，隨著時窗的不斷增大，測區(qū)SNR≥3數(shù)據(jù)占比不斷增大，隨后趨于穩(wěn)定；測區(qū)1.5

圖2 測線信號幅值圖Fig.2 Amplitude diagram of line signal

圖3 信噪比占比圖Fig.3 SNR ratio diagram

通過增大時窗，在一定程度上可以提高數(shù)據(jù)質(zhì)量，但同時分辨率有所降低。為在保持一定分辨率的基礎(chǔ)上，進一步提高數(shù)據(jù)質(zhì)量，本文中引入Robust-M估計方法。

Robust-M估計方法通過對大量數(shù)據(jù)進行迭代、剔除異常值，從而估計出一個準確值，來代表真實值。大地電磁探測通過在固定地點長時間進行數(shù)據(jù)采集，使得同一個測點處有大量數(shù)據(jù)，因此，Robust-M估計方法常用于大地電磁探測。而地空電磁探測采用無人機搭載接收系統(tǒng)進行數(shù)據(jù)采集，針對單個測點，僅有少量測量數(shù)據(jù)，并不能直接使用Robust-M估計方法估計真實值，所以需要對Robust-M估計方法進行改進，使其適用于地空頻率域電磁探測的數(shù)據(jù)處理。

2 基于Robust-M估計的干擾抑制方法

2.1 Robust-M估計算法

Robust-M估計是一種穩(wěn)健估計法，類似于最大似然估計，用于解決統(tǒng)計問題時[23]，通過觀測得到原始數(shù)據(jù)，表達式為

yi=xi+εi,i=1,2,…，m

(5)

式(5)中：m為原始數(shù)據(jù)個數(shù)；xi為原始數(shù)據(jù)的估計真實值；εi為誤差序列。

根據(jù)最大似然估計準則，選用分段連續(xù)的可微凸函數(shù)ρ(z)，構(gòu)造目標函數(shù)

(6)

為方便求解，定義

(7)

則通過計算，可求得yi的估計真實值xi。

2.2 權(quán)重的優(yōu)化設(shè)計

為確定固定大小的異常體對探測數(shù)據(jù)的影響范圍，進行仿真計算。異常體大小為100 m×200 m×1 250 m，埋深200 m，模型的yoz剖面如圖4所示，測線沿y軸方向，地空頻率域電磁探測仿真計算結(jié)果，如圖5所示，相對異常大于10%時視為有效異常。通過均勻大地與帶異常體仿真計算結(jié)果對比分析，可以發(fā)現(xiàn)測線方向100 m長的異常體，使異常體前后近500 m范圍內(nèi)的數(shù)據(jù)受到有效影響。所以在進行數(shù)據(jù)處理時，將測點前后數(shù)據(jù)聯(lián)合，進行測點響應(yīng)估計。

地空頻率域電磁探測通過無人機搭載線圈，進行測線數(shù)據(jù)采集，數(shù)據(jù)預(yù)處理后，得到磁感應(yīng)強度為

B(i)=BZ(i)+ε(i)

(8)

式(8)中：B(i)為磁感應(yīng)強度測量值；BZ(i)為磁感應(yīng)強度預(yù)估真實值；ε(i)為誤差序列；i為測點數(shù)。根據(jù)最大似然估計準則，對探測數(shù)據(jù)進行加窗分段后，構(gòu)造目標函數(shù)

(9)

求解目標函數(shù)，當(dāng)目標函數(shù)達到最小時，求解得到的BZ(i)即為估計結(jié)果。結(jié)合理論分析結(jié)果以及大地電磁信號具有連續(xù)性的特點，測點數(shù)據(jù)與測點前后數(shù)據(jù)具有連貫性，計算包含B(i)在內(nèi)的B(i-3)～B(i+3)七個測量值的中位數(shù)Bimd、標準差Bisd，令

(10)

BZ(i)=Bimd+φ(i)

(11)

圖4 異常體大小及位置Fig.4 Size and location of abnormal body

圖5 理論模擬結(jié)果Fig.5 Theoretical simulation results

穩(wěn)健M估計中，φ(x)函數(shù)有：中值函數(shù)、Huber函數(shù)、Hampel函數(shù)等多種類型。選用Hampel函數(shù)，即

(12)

式(12)中：a、b、c為調(diào)節(jié)常數(shù)[24]，定義a=1.2，b=3.5，c=8，通過對測點與測點前后數(shù)據(jù)應(yīng)用改進后的方法，估計測點響應(yīng)真實值。

3 應(yīng)用效果測試

3.1 模擬數(shù)據(jù)測試

如圖1所示，以發(fā)射源為坐標原點，建立空間直角坐標系，x軸沿發(fā)射導(dǎo)線方向，z軸垂直向下為正。求解帶有邊界條件的麥克斯韋方程組，計算空間位置(x,y,z)處的垂直電磁響應(yīng)。計算公式為

(13)

(14)

式中：μ0為真空磁導(dǎo)率；I為發(fā)射電流；L為導(dǎo)線長度；rTE為層狀大地反射系數(shù)；λ為漢克爾積分變量；J1為1階貝塞爾函數(shù)；ε0為真空中的介電常數(shù)；f為發(fā)射頻率。

通過均勻大地磁場表達式，運用MATLAB進行長導(dǎo)線源一維層狀介質(zhì)正演，獲得均勻大地垂直磁感應(yīng)強度。通過Comsol進行均勻大地仿真、帶異常體仿真。發(fā)射極矩40 000 A·m，測線沿y軸方向，大地電阻率100 Ω·m，異常體電阻率1 Ω·m，異常體大小為100 m×200 m×2 500 m，如圖6所示。

圖6 異常體位置及大小Fig.6 Size and location of abnormal body

電偶極源仿真計算32 Hz垂直磁場響應(yīng)幅度和相對異常曲線，如圖7所示。通過均勻大地電磁響應(yīng)與帶異常體電磁響應(yīng)對比，可以發(fā)現(xiàn)，測線方向100 m長的異常體將影響異常體前后約500 m范圍內(nèi)的磁場響應(yīng)，且磁場響應(yīng)在500 m范圍內(nèi)進行一個緩慢連續(xù)的變化；利用MATLAB產(chǎn)生服從正態(tài)分布的高斯隨機噪聲，如圖8所示，通過帶異常體32 Hz正演數(shù)據(jù)加噪前后對比圖，可以發(fā)現(xiàn)，隨機噪聲、尖峰噪聲在短時間內(nèi)會發(fā)生突變，并不符合異常體存在特點。對加噪后的帶異常體32 Hz正演數(shù)據(jù)分別進行平滑濾波與基于Robust-M估計的時窗疊加濾波處理，可以發(fā)現(xiàn)，如若在短時間內(nèi)發(fā)生多次突變，平滑濾波處理將會產(chǎn)生一段距離內(nèi)的數(shù)據(jù)起伏，與存在異常體的響應(yīng)特點相近，從而形成假異常，而基于Robust-M估計的時窗疊加濾波方法，結(jié)合測點前后數(shù)據(jù)情況，可以很好地將突變及短時間內(nèi)的多次突變進行降權(quán)處理，以降低突變點的影響，優(yōu)化數(shù)據(jù)。

3.2 實測數(shù)據(jù)測試

為了驗證基于Robust-M估計的時窗疊加濾波方法的有效性，在重慶市城口縣進行地空頻率域電磁探測，使用無人機搭載單分量線圈進行數(shù)據(jù)采集。測區(qū)概況如圖9所示。

圖7 理論模擬結(jié)果Fig.7 Theoretical simulation results

圖8 不同方法對加噪后信號處理結(jié)果Fig.8 Different methods to add noise signal processing results

測區(qū)位于重慶東北方向329 km，海拔1 500 m。發(fā)射電流20 A，接地導(dǎo)線1 000 m。發(fā)射頻率為32、64、128、256、512、1 024、2 048、4 096、8 192 Hz。共四條測線，測線間隔100 m，長2 800 m。四條測線相互平行，中軸線重復(fù)性測量。飛機飛行速度為6 m/s。

圖9 測區(qū)概況圖Fig.9 Survey area overview map

通過圖2中32～512 Hz頻率下的磁場響應(yīng)可以看出，因受到飛行因素及外部環(huán)境影響，實測信號內(nèi)存在隨機噪聲與尖峰噪聲，上下波動較大，并不連續(xù)，信號質(zhì)量較差。對測量數(shù)據(jù)分別進行時窗2 s、時窗10 s、時窗2 s-濾波處理、時窗2 s-Robust-M估計處理，對比應(yīng)用效果。32 Hz實測數(shù)據(jù)處理時域結(jié)果，如圖10所示。因大地的連續(xù)性，測得的磁場響應(yīng)也是連續(xù)的，將短時間內(nèi)出現(xiàn)的突變視為干擾。

通過圖10(a)、圖10(b)對比，可以發(fā)現(xiàn)，隨著時窗的增大，分辨率降低，可以有效減少隨機噪聲、尖峰噪聲；通過圖10(a)、圖10(c)、圖10(d)對比，可以發(fā)現(xiàn)，在時窗一定的情況下，平滑濾波與Robust-M估計方法對數(shù)據(jù)均具有一定的優(yōu)化作用，而平滑濾波根據(jù)理論數(shù)據(jù)處理結(jié)果，會因短時間內(nèi)的多次干擾而造成假異常，Robust-M估計可對異常值進行剔除或降權(quán)，時域圖與時窗10 s時的時域圖更為接近。

分別計算不同處理方法下的測區(qū)數(shù)據(jù)信噪比，統(tǒng)計結(jié)果，如表1所示。

通過時窗2 s處理結(jié)果與時窗10 s處理結(jié)果對比分析，可知，降低分辨率，可以有效減少不合格數(shù)據(jù)量，提高數(shù)據(jù)質(zhì)量。通過時窗2 s處理結(jié)果與時窗2 s-平滑濾波處理結(jié)果對比分析，可知，平滑濾波處理效果甚微。通過時窗2 s處理結(jié)果與時窗2 s-Robust-M估計處理結(jié)果對比，可知，基于Robust-M估計的時窗疊加處理方法，可以有效去除隨機噪聲、尖峰噪聲，提高數(shù)據(jù)質(zhì)量。

表1 測區(qū)信噪比統(tǒng)計結(jié)果Table 1 Statistical results of SNR in detection area

圖10 32 Hz實測數(shù)據(jù)處理時域圖Fig.10 32 Hz Time domain diagram of measured data processing

4 結(jié)論與展望

地空頻率域電磁探測采集數(shù)據(jù)中，存在眾多噪聲，需針對不同噪聲特點，進行抑制，從而提高數(shù)據(jù)質(zhì)量，提高探測的準確性。提出基于Robust-M估計的時窗疊加濾波方法，針對隨機噪聲、尖峰噪聲，具有很好的抑制效果。

(1)通過采用基于Robust-M估計的時窗疊加濾波處理，可以有效剔除數(shù)據(jù)中過大或過小的異常點。時窗10 s處理結(jié)果與時窗2 s-Robust-M估計處理結(jié)果對比，可知，SNR<1.5占比幾乎一致。通過對比不同處理方法的結(jié)果，基于Robust-M估計的時窗疊加濾波方法在保持一定分辨率的基礎(chǔ)上，有效提高了數(shù)據(jù)質(zhì)量。

(2)雖然Robust-M估計在抑制隨機噪聲與尖峰噪聲方面取得了較好的效果，但是當(dāng)測點數(shù)據(jù)存在較大的姿態(tài)噪聲時，需對電磁數(shù)據(jù)進行姿態(tài)校正后，再應(yīng)用Robust-M估計濾波方法，從而得到更準確的電磁數(shù)據(jù)。