基于EEMD 和DWT 的弱磁異常信號(hào)提取方法

寧文茜 ,王艷華 ,樊黎明,3* ,張曉峻 ,謝志臻

(1.西北工業(yè)大學(xué) 航海學(xué)院,陜西 西安,710077;2.自然資源部 海洋環(huán)境探測(cè)技術(shù)與應(yīng)用重點(diǎn)實(shí)驗(yàn)室,廣東 廣州,510300;3.西北工業(yè)大學(xué) 青島研究院,山東 青島,266200;4.哈爾濱工程大學(xué) 物理與光電學(xué)院,黑龍江 哈爾濱,150001)

0 引言

磁性目標(biāo)產(chǎn)生的磁場(chǎng)疊加在背景磁場(chǎng)上,使得其周?chē)臻g的磁場(chǎng)分布發(fā)生變化,從而形成磁異常。磁異常探測(cè)通過(guò)揭示環(huán)境磁場(chǎng)中的異常來(lái)檢測(cè)鐵磁物體,是一種廣泛應(yīng)用的被動(dòng)目標(biāo)探測(cè)方法[1]。該技術(shù)已在許多領(lǐng)域得到應(yīng)用,如未爆炸彈藥(unexploded ordnance,UXO)檢測(cè)[2-4]、海底電纜檢測(cè)[5]、交通監(jiān)視[6]及人體醫(yī)學(xué)研究[7]等。然而,磁異常會(huì)隨著距離的增加而迅速衰減。

近年來(lái),學(xué)者們提出多種磁異常檢測(cè)的方法。根據(jù)特點(diǎn),這些方法主要可以分為兩類(lèi)。一類(lèi)是基于信號(hào)分析的磁異常檢測(cè)方法[8-16]。Sheinker 等[8]提出標(biāo)準(zhǔn)正交基函數(shù)(orthonormal basis function,OBF)的匹配濾波,探測(cè)深埋在白高斯噪聲中的磁異常信號(hào)。為了應(yīng)對(duì)功率譜密度為1/fα(0 ＜α ＜2)的背景磁噪聲,其設(shè)計(jì)了白化濾波器,能夠提升磁異常的探測(cè)性能。該方法可以被認(rèn)為是自回歸(autoregressive,AR)過(guò)程,并且可以在高階過(guò)濾器下有效工作。Liu 等[11]分析了空間分布系數(shù)的約束,并提出一種改進(jìn)的基于空間投影系數(shù)變換的OBFs 方法用于磁異常檢測(cè)。Fan 等[13]提出磁異常梯度信號(hào)的OBF 探測(cè)方法,該方法能夠有效降低探測(cè)的虛警概率。假設(shè)磁噪聲的模式變化是由磁異常引起,最小磁熵法[9]可以揭示出該模式的變化,從而達(dá)到檢測(cè)磁異常的目的。該方法的優(yōu)點(diǎn)是不需要對(duì)目標(biāo)進(jìn)行先驗(yàn)假設(shè),計(jì)算復(fù)雜度較低。然而,該方法中噪聲概率密度函數(shù)通過(guò)統(tǒng)計(jì)靜止?fàn)顟B(tài)下的背景噪聲獲得。針對(duì)運(yùn)動(dòng)平臺(tái)的探測(cè)需求,Fan 等[14]提出基于經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)和最小熵的磁異常探測(cè)方法。該方法通過(guò)經(jīng)驗(yàn)?zāi)Ｊ椒纸饪煽焖儆?jì)算運(yùn)動(dòng)狀態(tài)下的近似噪聲概率密度函數(shù),實(shí)現(xiàn)對(duì)目標(biāo)的探測(cè)。為了提高低信噪比條件下的檢測(cè)性能,Wan 等[10]提出了一種基于隨機(jī)共振的磁異常檢測(cè)方法。然而,隨機(jī)共振參數(shù)的選取影響了該方法的有效性。此外,學(xué)者們還提出了混合算法來(lái)檢測(cè)磁異常[15]。

另一類(lèi)磁異常檢測(cè)方法是基于機(jī)器學(xué)習(xí)的磁異常檢測(cè)[17-20]。該類(lèi)方法的主要思想是將磁異常檢測(cè)視為分類(lèi)問(wèn)題,即確定磁信號(hào)中是否存在異常信號(hào)。Liu 等[19]提出了一種使用全連接神經(jīng)網(wǎng)絡(luò)的磁異常檢測(cè)方法,其中選擇了2 種類(lèi)型的特征作為輸入。Hu 等[17]提出了一種基于卷積神經(jīng)網(wǎng)絡(luò)的磁異常檢測(cè)方法,將特征信息的提取和分類(lèi)任務(wù)結(jié)合在一個(gè)單一體系中。Fan 等[18]提出了一種使用支持向量機(jī)的自適應(yīng)磁異常檢測(cè)方法,其中OBF 能量和磁熵作為磁異常特征輸入到支持向量機(jī)。

通過(guò)磁噪聲抑制、磁信號(hào)增強(qiáng)或磁異常多特征信息融合等方式,上述方法能實(shí)現(xiàn)在有效探測(cè)范圍內(nèi)識(shí)別磁異常信號(hào),判斷是否存在磁性目標(biāo)。當(dāng)檢測(cè)到磁性目標(biāo)后,在進(jìn)一步定位與識(shí)別磁性目標(biāo)時(shí),需要獲取目標(biāo)的弱磁異常信號(hào)。針對(duì)低信噪比下弱磁異常信號(hào)的獲取問(wèn)題,文中提出基于集合經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition,EEMD)和離散小波變換(discrete wavelet transformation,DWT)的弱磁異常提取方法。

1 弱磁異常信號(hào)分析

當(dāng)磁性目標(biāo)與磁傳感器的距離大于目標(biāo)最大尺度的3 倍時(shí),可以把磁性目標(biāo)視為一個(gè)磁偶極子。磁性目標(biāo)產(chǎn)生的磁場(chǎng)可以表示為[21]

在實(shí)際測(cè)量中,測(cè)量磁場(chǎng)Bm由地磁場(chǎng)Be和目標(biāo)產(chǎn)生的磁場(chǎng)Ba組成。由于磁場(chǎng)隨距離的快速衰減,當(dāng)目標(biāo)遠(yuǎn)離磁傳感器時(shí),Ba遠(yuǎn)小于Be。當(dāng)利用標(biāo)量磁力儀探測(cè)遠(yuǎn)距離的目標(biāo)時(shí),磁場(chǎng)的標(biāo)量形式可以表示為

由式(2),目標(biāo)產(chǎn)生磁異常的標(biāo)量形式可以表示為

由式(1)和式(3)可知,目標(biāo)產(chǎn)生的磁異常信號(hào)隨距離的3 次方衰減。隨著探測(cè)距離的增大,目標(biāo)磁異常信號(hào)將會(huì)變得微弱。因此,目標(biāo)弱磁異常提取是開(kāi)展遠(yuǎn)距離目標(biāo)定位與識(shí)別的基礎(chǔ)。

2 弱磁異常提取方法

2.1 方法概述

當(dāng)探測(cè)遠(yuǎn)距離磁性目標(biāo)時(shí),由于磁性目標(biāo)產(chǎn)生的磁場(chǎng)隨著距離的增加快速衰減,導(dǎo)致目標(biāo)的弱磁異常淹沒(méi)在磁噪聲中。通過(guò)磁異常探測(cè)方法確定目標(biāo)存在后,需要對(duì)弱磁異常信號(hào)進(jìn)行有效的提取,從而為后續(xù)的目標(biāo)定位與識(shí)別提供高精度的數(shù)據(jù)基礎(chǔ)。通常認(rèn)為,磁異常信號(hào)的頻率位于極低頻范圍[22]。

為有效抑制背景磁噪聲并提取目標(biāo)產(chǎn)生的弱磁異常信號(hào),提出基于EEMD 和DWT 的弱磁異常提取方法。該方法框架圖如圖1 所示。首先,利用EEMD 將磁信號(hào)分解成一系列固有模態(tài)函數(shù)(intrinsic mode function,IMF)項(xiàng)和殘差項(xiàng)。分解出來(lái)的各IMF 分量包含了原信號(hào)不同時(shí)間尺度的局部特征信號(hào),包括從高到低的不同頻段。根據(jù)噪聲評(píng)估結(jié)果,IMF 被分為噪聲域和信號(hào)域。噪聲域的IMF 需要進(jìn)行降噪處理,獲得部分磁異常信號(hào)。其次,利用DWT 將噪聲域IMF 分解為低頻帶的近似系數(shù)(AC)和相對(duì)較高頻帶上的一系列細(xì)節(jié)系數(shù)(DC)[23]。保留低頻帶的近似系數(shù),利用逆離散小波變換(inverse discrete wavelet transform,IDWT)重新構(gòu)建噪聲域的磁異常。最后,將信號(hào)域磁異常信號(hào)項(xiàng)和重構(gòu)的噪聲域磁異常信號(hào)合成去噪的磁異常信號(hào),從而實(shí)現(xiàn)微弱磁異常信號(hào)的提取,提高磁異常的信噪比。

圖1 基于EEMD-DWT 的弱磁異常提取方法框圖Fig.1 Block diagram of weak magnetic anomaly extraction method based on EEDM-DWT

2.2 EEMD

EMD 是由Huang 等[24]提出的一種處理非線(xiàn)性非平穩(wěn)信號(hào)方法。該方法可將復(fù)雜的時(shí)間序列信號(hào)分解為少數(shù)幾個(gè)具有不同特征尺度的時(shí)間序列,稱(chēng)為 IMF 和一個(gè)趨勢(shì)項(xiàng)。因此,信號(hào)可以表示為

式中:IIMF,i(t)為第i個(gè)本征模函數(shù);r(t)為趨勢(shì)項(xiàng)。

EMD 方法從數(shù)據(jù)本身特征出發(fā),不依賴(lài)于任何外部函數(shù)和參數(shù),具有自適應(yīng)性。然而,由于某些信號(hào)具有間歇性,分解時(shí)常常出現(xiàn)模態(tài)混合現(xiàn)象,即單獨(dú)一個(gè)IMF 分量包含著不同尺度的信號(hào)或一種尺度的信號(hào)分布在不同的IMF 分量中,導(dǎo)致分解出的IMF 分量缺乏物理意義。

為克服這個(gè)問(wèn)題,Wu 等[25]發(fā)展了EEMD 方法。該方法是通過(guò)在原信號(hào)中添加白噪聲序列,再利用EMD 方法將添加了白噪聲序列的信號(hào)分解成若干個(gè)IMF,重復(fù)多次且每次添加不同的白噪聲序列,最后取多次分解得到的相應(yīng)IMF 的算術(shù)平均值作為最終結(jié)果。該方法依據(jù)數(shù)據(jù)自身的時(shí)間尺度特征來(lái)進(jìn)行信號(hào)分解,可以應(yīng)用于任何類(lèi)型信號(hào)的分解,在處理非平穩(wěn)及非線(xiàn)性數(shù)據(jù)上,具有非常明顯的優(yōu)勢(shì),適合于分析非線(xiàn)性、非平穩(wěn)信號(hào)序列,具有很高的信噪比。

通過(guò)EEMD 所分解出來(lái)的各IMF 分量包含了原信號(hào)的不同時(shí)間尺度的局部特征信號(hào)。根據(jù)噪聲評(píng)估結(jié)果,IMF 被分為噪聲域和信號(hào)域。信號(hào)域的IMF 需保留,噪聲域中的IMF 需要進(jìn)行降噪處理。

基于閾值信息,開(kāi)展IMF 噪聲域和信號(hào)域的區(qū)分。Kopsinis 等[26]提出了一種魯棒的閾值選擇策略,并推導(dǎo)了分解中分形高斯噪聲的方差傳遞模型。根據(jù)噪聲方差,推導(dǎo)出噪聲能量在各階IMF中的傳遞模型,然后根據(jù)各階IMF 中的噪聲能量確定相應(yīng)的閾值。最后,根據(jù)閾值信息,確定對(duì)應(yīng)的IMF 為信號(hào)域或噪聲域。

通常,第1 層IMF 主要成分是噪聲,噪聲的能量可以通過(guò)噪聲的標(biāo)準(zhǔn)差來(lái)估計(jì),即

其他IMF 的噪聲能量通過(guò)傳遞模型獲得,即

式中,CH=(E1/βH),文獻(xiàn)[27]給出了 βH和 ρH的取值說(shuō)明。根據(jù)噪聲能量,可通過(guò)構(gòu)造閾值函數(shù)獲得閾值。常用的閾值計(jì)算模型為

式中:Tk為第k層IMF 對(duì)應(yīng)的閾值;Ek為第k層IMF 的噪聲能量;N表示信號(hào)長(zhǎng)度。

根據(jù)閾值信息,區(qū)分信號(hào)域和噪聲域IMF 的步驟如下:

2.3 DWT 方法

DWT 具備處理信號(hào)低頻部分的能力,同時(shí)能夠較好地分解高頻部分,收集信號(hào)各個(gè)頻段的有用信息,進(jìn)而提高去噪精度。因此,DWT 可將信號(hào)x(t)的時(shí)域表示分解為母小波的多個(gè)時(shí)移和縮放。數(shù)學(xué)上,小波ψa,b(t)由它的母小波ψ(·)通過(guò)縮放和平移得到,即

式中:a∈R+為尺度參數(shù);b∈R為平移參數(shù)。將參數(shù)a和b離散化得到DWT。在DWT 中,根據(jù)2 的冪次選擇參數(shù):a=2j,b=k2j,其中k,j∈Z。因此,信號(hào)的DWT 可寫(xiě)成

通常,小波細(xì)節(jié)系數(shù)表征高頻信號(hào)部分,小波近似系數(shù)表征低頻信號(hào)部分。由于磁異常信號(hào)的頻率位于極低頻范圍。因此,對(duì)噪聲域部分xn(t)需要進(jìn)行DWT,獲得近似系數(shù)和細(xì)節(jié)系數(shù)。通過(guò)利用近似系數(shù)可重構(gòu)出低頻的弱磁異常信號(hào)實(shí)現(xiàn)多噪聲域IMF 中的高頻磁噪聲抑制。

2.4 弱磁異常信號(hào)提取

低信噪比的弱磁異常信號(hào)經(jīng)過(guò)EEMD 和DWT的去噪處理,獲得去噪后的磁異常信號(hào)為

為了量化該方法的性能,使用信噪比 (signalnoise-ratio,SNR)和均方根誤差(root mean squared error,RMSE)進(jìn)行評(píng)價(jià)。評(píng)價(jià)指標(biāo)的數(shù)學(xué)表達(dá)式為

式中:s(t)為目標(biāo)的磁異常信號(hào);x(t)為含有噪聲的磁異常信號(hào);x～(t)為噪聲抑制后的磁異常信號(hào)。

3 實(shí)驗(yàn)驗(yàn)證

3.1 仿真實(shí)驗(yàn)

為了有效評(píng)價(jià)該方法的性能,設(shè)計(jì)1 個(gè)一維傳感器陣列(如圖2 所示)來(lái)收集磁背景噪聲。磁傳感器為高靈敏度的光泵磁力儀,其本征噪聲約為選擇哈爾濱市金沙灘濕地公園采集了真實(shí)的磁噪聲樣本,該區(qū)域環(huán)境磁活動(dòng)非常低。地磁場(chǎng)局地傾角和偏角分別為63.3°和-10.1°。地磁場(chǎng)的振幅約為55 200 nT。磁力儀的采樣頻率設(shè)置為10 Hz。

圖2 一維磁傳感器陣列Fig.2 One-dimensional magnetic sensor array

典型的磁異常信號(hào)是由仿真目標(biāo)產(chǎn)生的。目標(biāo)以6.4 m/s 的恒定速度平行于X軸運(yùn)動(dòng),從 (191,40,5) m 開(kāi)始,在 (-192,40,5) m 結(jié)束。目標(biāo)的磁矩在運(yùn)動(dòng)過(guò)程中保持不變,磁矩的大小為 200 A·m2。通過(guò)仿真可以得到典型的磁異常,然后將傳感器陣列測(cè)量的真實(shí)磁噪聲添加到模擬信號(hào)中,可以獲得具有真實(shí)噪聲和模擬磁異常的合成磁信號(hào),如圖3 所示。

圖3 合成的磁異常信號(hào)Fig.3 Synthetic magnetic anomaly signal

基于EEMD-DWT 方法和基于EEMD 方法用于合成磁異常信號(hào)的信號(hào)提取。在基于EEMDDWT 方法中,噪聲能量傳遞模型中相關(guān)參數(shù)設(shè)置如下:H=0.5,βH=0.719,ρH=2.010。嘗試多組小波函數(shù),確定最優(yōu)的小波函數(shù)為db4。根據(jù)磁力儀的采樣頻率和磁異常的頻率范圍,確定最優(yōu)的小波分解層數(shù)為5 層?；贓EMD 方法中的參數(shù)設(shè)置與基于EEMD-DWT 方法中的參數(shù)一致,分別利用這2 種方法對(duì)合成磁異常信號(hào)進(jìn)行處理,提取磁異常的結(jié)果如圖4 所示。通過(guò)圖4 可知,基于EEMD-DWT 方法提取的磁異常信號(hào)更接近于目標(biāo)信號(hào)。在基于EEMD 方法中,信噪比為1.72 dB。在基于EEMD-DWT 方法中,信噪比為5.27 dB,比基于EEMD 方法高出3.56 dB。

圖4 不同方法的弱磁異常提取結(jié)果Fig.4 Results of weak magnetic anomaly signal extraction using different methods

2 種方法獲得的RMSE 結(jié)果如圖5 所示。基于EEMD-DWT 方法的RMSE 值為0.000 9,磁異常最大值處偏差的誤差為4.36%?；贓EMD 方法的RMSE 值為0.002 4,磁異常最大值處偏差的誤差為29.20%。因此,基于EEMD-DWT 方法能夠有效提取弱磁異常信號(hào)。

圖5 不同方法的RMSE 結(jié)果Fig.5 Results of RMSE using different methods

3.2 外場(chǎng)試驗(yàn)

為了驗(yàn)證所提方法在實(shí)際磁測(cè)中的有效性,在試驗(yàn)中將其應(yīng)用于實(shí)際測(cè)量的弱磁異常提取,如圖6 所示。在磁噪聲采集試驗(yàn)中,磁測(cè)量系統(tǒng)由2 個(gè)銫光泵磁力儀(CS-L 和Scintrex)、1 個(gè)數(shù)據(jù)處理單元(MMS-4,PICO)和1 臺(tái)計(jì)算機(jī)組成。該測(cè)量系統(tǒng)以10 Hz 的采樣率用于獲取磁場(chǎng)信號(hào)。磁測(cè)量系統(tǒng)放置在河道岸邊,河道中行駛的拉沙船、摩托艇等作為磁性目標(biāo)。

圖6 河道岸邊磁場(chǎng)測(cè)量現(xiàn)場(chǎng)Fig.6 Scene of magnetic field measurement at river bank

拉沙船的體積大,磁性材料多,其產(chǎn)生的磁異常信號(hào)較大,背景磁噪聲對(duì)其磁異常信號(hào)影響有限。試驗(yàn)不對(duì)拉沙船的磁異常信息進(jìn)行處理。由于摩托艇的體積較小,磁性材料少,產(chǎn)生的磁異常較弱,磁噪聲對(duì)其磁異常信號(hào)影響大。因此,選擇摩托艇產(chǎn)生的磁異常信號(hào)進(jìn)行分析和處理。

利用基于EEMD-DWT 方法對(duì)摩托艇產(chǎn)生的磁異常信號(hào)進(jìn)行處理,獲得目標(biāo)的磁異常信號(hào),如圖7 所示。試驗(yàn)結(jié)果顯示,該方法能夠有效抑制背景磁噪聲的影響,獲得目標(biāo)的磁異常信號(hào)。相比于基于EEMD 方法,所提出的基于EEMD-DWT方法在弱磁異常信號(hào)提取中具有明顯的優(yōu)勢(shì)。

圖7 實(shí)測(cè)磁異常信號(hào)提取Fig.7 Measured magnetic anomaly signal extraction

4 結(jié)束語(yǔ)

為獲取遠(yuǎn)距離目標(biāo)產(chǎn)生的微弱磁異常,針對(duì)弱磁異常信號(hào)的極低頻特性,文中提出基于EEMDDWT 的弱磁異常提取方法。采用EEMD 將微弱磁異常信號(hào)分解為信號(hào)域信號(hào)和噪聲域信號(hào)。噪聲域信號(hào)通過(guò)DWT 進(jìn)行低頻信號(hào)的提取。利用信號(hào)域信號(hào)和經(jīng)DWT 處理噪聲域信號(hào)的疊加,獲得目標(biāo)的弱磁異常信號(hào)。在仿真實(shí)驗(yàn)中,基于EEMDDWT 方法提取弱磁異常與真實(shí)磁異常的RMSE值為0.000 9,磁異常最大值處的誤差為4.36%。外場(chǎng)試驗(yàn)中的結(jié)果顯示基于EEMD-DWT 方法提取弱磁異常信號(hào)與實(shí)測(cè)數(shù)據(jù)的變化趨勢(shì)一致。該方法提取的磁異常信號(hào)可為后續(xù)目標(biāo)的定位與識(shí)別提供必要的數(shù)據(jù)基礎(chǔ)。