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

(1.吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院， 長(zhǎng)春 130012； 2.地球信息探測(cè)儀器教育部重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)春 130012)

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

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

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

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

1 理論基礎(chǔ)

1.1 地空頻率域電磁探測(cè)原理

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

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

(1)

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

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

(2)

(3)

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

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

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

(4)

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

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

圖2 測(cè)線信號(hào)幅值圖Fig.2 Amplitude diagram of line signal

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

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

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

2 基于Robust-M估計(jì)的干擾抑制方法

2.1 Robust-M估計(jì)算法

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

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

(5)

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

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

(6)

為方便求解，定義

(7)

則通過計(jì)算，可求得yi的估計(jì)真實(shí)值xi。

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

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

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

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

(8)

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

(9)

求解目標(biāo)函數(shù)，當(dāng)目標(biāo)函數(shù)達(dá)到最小時(shí)，求解得到的BZ(i)即為估計(jì)結(jié)果。結(jié)合理論分析結(jié)果以及大地電磁信號(hào)具有連續(xù)性的特點(diǎn)，測(cè)點(diǎn)數(shù)據(jù)與測(cè)點(diǎn)前后數(shù)據(jù)具有連貫性，計(jì)算包含B(i)在內(nèi)的B(i-3)～B(i+3)七個(gè)測(cè)量值的中位數(shù)Bimd、標(biāo)準(zhǔn)差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估計(jì)中，φ(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，通過對(duì)測(cè)點(diǎn)與測(cè)點(diǎn)前后數(shù)據(jù)應(yīng)用改進(jìn)后的方法，估計(jì)測(cè)點(diǎn)響應(yīng)真實(shí)值。

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

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

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

(13)

(14)

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

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

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

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

3.2 實(shí)測(cè)數(shù)據(jù)測(cè)試

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

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

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

測(cè)區(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。共四條測(cè)線，測(cè)線間隔100 m，長(zhǎng)2 800 m。四條測(cè)線相互平行，中軸線重復(fù)性測(cè)量。飛機(jī)飛行速度為6 m/s。

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

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

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

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

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

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

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

4 結(jié)論與展望

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

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

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