基于形態(tài)學-HHT算法的船載地磁三分量信號分析與預處理

刁云云，高金耀*，吳國超，蔡曉仙，岳 梅

(1.自然資源部 第二海洋研究所，浙江 杭州 310012;2.自然資源部 海底科學重點實驗室，浙江 杭州 310012;3.浙江大學 海洋學院，浙江 杭州 310058)

0 引言

地磁場是一種基本的地球物理場，通過X、Y、Z三個方向分量可反映其空間分布特征[1]，時間上基本穩(wěn)定，變化頻率接近0 hz。拖曳式地磁測量是海洋磁力信息采集的主要手段，但在南、北極地磁的測量中，受限于水深環(huán)境、浮冰、作業(yè)效率等因素，通常采用航空或船載地磁三分量系統(tǒng)(STCM)測量。STCM為日本在1972—1975年間開發(fā)，已在多個海區(qū)采集了矢量地磁數(shù)據(jù)，精度為50±25 nT[2-7]，測得的地磁信號包含真實地磁場、感應磁場、固有磁場、渦流場以及隨機磁場，其中隨機磁場來自磁力儀的測量噪聲、船載電氣設備的高頻電子干擾、海底沉船、水面船只等鐵磁性物質(zhì)的短時強干擾[8]以及船體轉(zhuǎn)向的脈沖干擾等。

濾波處理是消除隨機磁場信號的重要手段。傳統(tǒng)的濾波方法是對信號進行快速傅里葉變換(Fast Fourier Transform,FFT)獲得頻譜圖，設置截止頻率進行去噪處理獲得目標信號。但對于頻率復雜、隨時間非線性、非平穩(wěn)變化的極地地磁信號，該方法不適用。形態(tài)學濾波是基于數(shù)學形態(tài)學理論發(fā)展起來的一種針對非線性、非穩(wěn)定變化信號的處理方法[9]，可用于地球物理信號在時域上的信噪分離。如陳輝 等[10]用該方法對地震信號進行去噪，獲得了較高信噪比的地震剖面信息。希爾伯特-黃變換(Hilbert-Huang Transform,HHT)是近年發(fā)展起來的一種新的信號處理方法[11-16]，與FFT相比，HHT具有保留原信號頻率的優(yōu)勢[17]。BATTISTA et al[12]研究表明，針對電氣設備的雜散場對地磁的干擾，HHT去噪效果較好；ZHOU et al[18]改進了HHT算法，通過引入差分磁場，增強了HHT的噪聲識別能力和去噪效果；QIAO et al[19]優(yōu)化了HHT算法中EMD分解算法；李季 等[8]發(fā)現(xiàn)HHT與形態(tài)學濾波結合可以有效降低模擬信號中隨機噪聲的干擾?？傮w上，對于實測地磁信號分析和處理的相關研究尚處于探索階段。本文構建了一種形態(tài)學濾波和HHT變換相結合的濾波處理方法，對實測船載地磁三分量信號進行去噪和分析。

1 數(shù)據(jù)和方法

1.1 數(shù)據(jù)來源

第29次南極科考在普里茲灣海域同步采集了 2 443 km 船載地磁三分量和拖曳地磁總場數(shù)據(jù)，本研究分析的P6-2測線全長約187 km，呈南北向分布，采樣時間約2.5 h。船載地磁三分量由Grad-03-500M三軸磁力梯度儀采集，儀器相關參數(shù)如表1所示。拖曳地磁總場由G880型標量地磁傳感器采集，兩者采集時間同步。P6-2測線數(shù)據(jù)見圖1。

圖1 第29次南極科考船載地磁三分量測線數(shù)據(jù)Fig.1 Three-component geomagnetic line data of the 29th Antarctic Scientific Expedition

三軸磁力梯度儀安裝在船尾部風廓儀平臺的桅桿上，平臺距二層甲板15 m高，X分量和船艏向一致，Y分量指向右舷，Z分量垂直向下。姿態(tài)儀傳感器用于采集船的姿態(tài)，與三軸磁力梯度儀相距約 10 m，兩者同步變化。

本文以總場數(shù)據(jù)為例展示信號分析過程，總場的磁場值為50 000～51 500 nT，依據(jù)下式計算得到：

(1)

式中：X、Y、Z分別代表船載地磁測量的東向、北向、垂直地心三個方向的地磁分量，單位：nT；M代表地磁總場，單位：nT。

1.2 形態(tài)學濾波

采用形態(tài)學濾波分離船載地磁三分量中的脈沖信號。結構元素的選取是形態(tài)學濾波一個重要因素，其幅值一般應大于噪聲幅值，小于信號主要輪廓高度的1%。參考陳輝 等[10]和柏林 等[20]設計的開-閉和閉-開組合濾波器，其中開運算可以壓制峰值處的脈沖，閉運算能過濾低谷的脈沖干擾，組合濾波器公式為

y(n)={OC[f(n)+COf(n)]}/2

(2)

式中：OC表示開-閉運算，CO表示閉-開運算,f(n)代表原始信號,y(n)表示形態(tài)學濾波結果。

1.3 希爾伯特-黃變換(HHT)

HHT包括經(jīng)驗模態(tài)分解(EMD)和Hilbert變換兩部分。EMD將非線性、非平穩(wěn)的地磁信號自適應地分解為若干個線性、穩(wěn)定、固有模態(tài)函數(shù)(IMF)。EMD分解判定條件有兩個：(1)給定信號中，過零點數(shù)量和極值點數(shù)量相同或者相差1個；(2)在給定信號中，上下包絡線的均值為零，具體計算方法可參考文獻[8]和[12] 。IMF階次越大代表地磁信號中低頻、穩(wěn)定成分越多；相反，IMF階次越小代表中、高頻成分占主要部分。

利用Hilbert變換計算各IMF分量的Hilbert譜，得到各信噪頻段的分布。

EMD分解的地磁信號，通過Hilbert變換[21]獲得時頻分析，公式為

(3)

根據(jù)式(3)可進一步得到：

(4)

(5)

式中：ai為瞬時振幅；φi為瞬時相位；fi為瞬時頻率，單位：Hz。

對各IMF含噪信號進行自相關分析，公式為

(6)

式中：E為數(shù)學期望，σ為標準差。

理想的白噪聲在t2=t1時，R(τ)=1，其他位置為0。地磁信號的高頻IMF分量以噪聲成分居多，其自相關函數(shù)表現(xiàn)為以零點為中心，向兩側(cè)衰減為0；低頻IMF分量以真實地磁信號為主，零點位置附近衰減緩慢[22]。

1.4 濾波效果檢驗

對形態(tài)學-HHT濾波和FFT濾波后的地磁總場進行濾波效果真實性檢驗。分別對兩種濾波后的總場數(shù)據(jù)進行船磁校正計算，以拖曳地磁總場為基準計算均方根誤差(RMSE)，公式如下

(7)

式中：n代表數(shù)據(jù)個數(shù)，Mcal,i代表船磁校正后的地磁總場，Mty,i代表拖曳地磁總場。

2 結果和分析

2.1 形態(tài)學濾波

結合圖1總場中脈沖信號的幅值和寬度情況，選取了三角結構元素為：[0，50，100，200，100，50，0],并采用式(2)的組合濾波方式，其計算結果如圖2所示。

與原始地磁總場(圖1a)相比，形態(tài)學濾波后的地磁總場(圖2)濾除了大部分的脈沖點(在0 s和7 000～9 000 s之間達到51 000 nT以上的點)，同時保留了真實地磁信號(50 000～51 000 nT)，表明形態(tài)學濾波對船載地磁三分量數(shù)據(jù)有良好的去噪能力。

2.2 EMD分解結果

將形態(tài)學濾波處理后的總場數(shù)據(jù)(圖2)經(jīng)EMD分解成11階次IMF分量和一個殘余分量(圖3)，IMF1～IMF6在整個時間序列上幅值較均勻，且主要集中在±100 nT以內(nèi)，推斷集中了大量高頻交變磁場噪聲。IMF7～IMF11分量幅值變化平緩，代表了總場信號的低頻成分，而殘余分量幅值在50 000 nT左右，為總場信號的主要成分。

圖2 總場形態(tài)學濾波結果Fig.2 The total field result after morphology filtering

圖3 EMD分解結果Fig.3 Decomposition of total field signal by EMD

圖4展示了IMF各分量自相關結果，在零時刻，IMF1～IMF5各分量幅值達到最大，自相關程度最大，但在偏離零時刻的其他時間段，幅值急劇衰減到0值附近，表明這些分量的自相關程度小，說明這些分量是由噪聲主導的；偏離零時刻的IMF6～IMF11分量幅值變化緩慢(越低頻的分量，幅值變化越平緩)，且大部分不接近0，自相關程度高，表明IMF6～IMF11分量由地磁信號主導。

圖4 IMF各分量自相關函數(shù)Fig.4 Autocorrelation function of each IMF

2.3 希爾伯特譜分析

IMF1～IMF6分量的希爾伯特譜如圖5所示，IMF1～IMF2、IMF3～IMF4以及IMF5分量對應的希爾伯特譜(圖5 a～5c)，其頻率范圍逐步變低，依次分別是：0.030～0.460 Hz、0.010～0.100 Hz、0.005～0.030 Hz。IMF6分量信號頻率整體接近0 Hz，但仍含有和0 Hz接近的噪聲成分，從IMF7開始，信號都為0 Hz (圖略)，說明從IMF6分量開始，信號由相對穩(wěn)定的地磁場構成。圖5e為地磁總場的希爾伯特譜分布，在整個時間段，噪聲和地磁場信號混疊在一起。

希爾伯特譜表明，主要地磁信號在0.030 Hz以下(圖5a中紅色直線)，噪聲范圍為0.010～0.500 Hz(圖5a～5c)。由原始三分量的總場計算的傅里葉譜(圖5f)可知，信號頻率在0～0.500 Hz之間，頻帶范圍較寬，與希爾伯特譜相比，傅里葉譜無法得知各頻率所對應的采樣時間。

2.4 形態(tài)學-HHT濾波后的地磁總場

根據(jù)EMD分解結果和譜分析，對IMF1～IMF6分量分別采用巴特沃斯低通濾波[23]，低通濾波設計參數(shù)為：通帶最大衰減為3 db、阻帶最大衰減為60 db、通帶截止頻率為0.01～0.03 Hz、阻帶截止頻率為0.10 Hz。

圖6b中的綠色曲線代表形態(tài)學-HHT濾波預處理后重構的地磁總場圖，包含了形態(tài)學濾波后經(jīng)低通濾波的IMF1～IMF6分量和未濾波低頻分量(IMF7～IMF11及殘余分量)兩個部分。與FFT濾波預處理相比，形態(tài)學-HHT濾波預處理后，數(shù)據(jù)整體較光滑，振蕩程度較小，符合地磁場變化特性，大部分脈沖被有效壓制，如0 s附近和7 000～9 000 s之間大于51 000 nT 的點(圖1)，高頻噪聲也被大幅削減。而圖6a中的黑色曲線(FFT濾波后的地磁總場)，雖然去掉了部分大于0.1 Hz的高頻噪聲(圖5f)和脈沖信號，但仍含有較多低頻噪聲成分。

圖5 各個IMF分量及總場信號頻譜圖Fig.5 Spectrum of each IMF and total field signal注：(a)代表IMF1和IMF2分量疊加的希爾伯特譜；(b)代表IMF3和IMF4分量疊加的希爾伯特譜；(c)和(d)分別為IMF5、IMF6對應的希爾伯特譜；(e)為總場信號的希爾伯特譜；(f)為總場信號的傅里葉譜。圖a和b中紅色細線指示分量信號最低頻率界限。各分量圖下方顏色條帶代表振幅能量，單位：nT。Note:(a)represents the Hilbert spectrums overlying by IMF1 and IMF2 together;(b)represents the Hilbert spectrums overlying by IMF 3 and IMF 4 together;(c)and (d)are the Hilbert spectrums corresponding to IMF5 and IMF6 respectively;(e)are the Hilbert spectrums of all IMFs signal;(f)is the Fourier spectrum of total field signal by FFT.The red thin lines in fig.a and fig.b indicate the lowest frequency limit of the component signal.The color bars below graph represent the amplitude energy,unit:nT.

圖6 兩種濾波處理后的總場數(shù)據(jù)Fig.6 The total field signal data after filtering by two different methods

2.5 濾波效果的真實性檢驗

雖然形態(tài)學-HHT濾波壓制了隨機噪聲的干擾，但與FFT濾波方法一樣，濾波處理后的地磁信號中仍包含低頻船磁的影響。圖7為對2種濾波處理后的地磁信號進行船磁校正后的地磁總場結果。船磁校正后，形態(tài)學-HHT濾波的地磁信號整體起伏平緩，僅在局部存在“震蕩”現(xiàn)象，除初始階段(0～200 s)與拖曳地磁總場相差較大外，其他時段兩者起伏趨勢一致。而FFT濾波后的地磁總場信號整體震蕩大，與拖曳地磁總場存在較大差距。分別計算船磁校正后2種濾波方法的地磁總場與拖曳地磁總場的均方根誤差(RMSE)，形態(tài)學-HHT濾波的RMSE約為194 nT，F(xiàn)FT濾波的RMSE約為600 nT，相差近3倍。

圖7 兩種濾波處理后的船磁補償結果對比Fig.7 The comparison of ship magnetic compensation by two different filtering methods

總體上，與FFT濾波相比，形態(tài)學-HHT濾波處理后地磁信號形態(tài)和幅值變化更接近拖曳地磁總場，誤差值更小，更適合用于船載三分量地磁數(shù)據(jù)的預處理。值得注意的是，船磁校正效果除了與預處理方法有關外，還與船磁校正算法、數(shù)據(jù)質(zhì)量等息息相關。

3 結論

船載三分量地磁信號中包含了除船磁外的大量隨機噪聲干擾，形態(tài)學-HHT方法能夠在時間尺度有效識別信號的頻譜成分，自適應地進行信號分解，降低工作量，提高譜分析有效性，解決脈沖噪聲干擾，對非平穩(wěn)、非線性的船載地磁信號噪聲壓制優(yōu)于FFT濾波算法，通過與拖曳地磁總場測量結果的比較也證明了該濾波方法的有效性。