刁云云，高金耀*，吳國(guó)超，蔡曉仙，岳 梅

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

0 引言

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

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

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

1.1 數(shù)據(jù)來(lái)源

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

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

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

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

(1)

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

1.2 形態(tài)學(xué)濾波

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

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

(2)

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

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

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

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

EMD分解的地磁信號(hào)，通過(guò)Hilbert變換[21]獲得時(shí)頻分析，公式為

(3)

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

(4)

(5)

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

對(duì)各IMF含噪信號(hào)進(jìn)行自相關(guān)分析，公式為

(6)

式中：E為數(shù)學(xué)期望，σ為標(biāo)準(zhǔn)差。

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

1.4 濾波效果檢驗(yàn)

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

(7)

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

2 結(jié)果和分析

2.1 形態(tài)學(xué)濾波

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

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

2.2 EMD分解結(jié)果

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

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

圖3 EMD分解結(jié)果Fig.3 Decomposition of total field signal by EMD

圖4展示了IMF各分量自相關(guān)結(jié)果，在零時(shí)刻，IMF1～I(xiàn)MF5各分量幅值達(dá)到最大，自相關(guān)程度最大，但在偏離零時(shí)刻的其他時(shí)間段，幅值急劇衰減到0值附近，表明這些分量的自相關(guān)程度小，說(shuō)明這些分量是由噪聲主導(dǎo)的；偏離零時(shí)刻的IMF6～I(xiàn)MF11分量幅值變化緩慢(越低頻的分量，幅值變化越平緩)，且大部分不接近0，自相關(guān)程度高，表明IMF6～I(xiàn)MF11分量由地磁信號(hào)主導(dǎo)。

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

2.3 希爾伯特譜分析

IMF1～I(xiàn)MF6分量的希爾伯特譜如圖5所示，IMF1～I(xiàn)MF2、IMF3～I(xiàn)MF4以及IMF5分量對(duì)應(yīng)的希爾伯特譜(圖5 a～5c)，其頻率范圍逐步變低，依次分別是：0.030～0.460 Hz、0.010～0.100 Hz、0.005～0.030 Hz。IMF6分量信號(hào)頻率整體接近0 Hz，但仍含有和0 Hz接近的噪聲成分，從IMF7開(kāi)始，信號(hào)都為0 Hz (圖略)，說(shuō)明從IMF6分量開(kāi)始，信號(hào)由相對(duì)穩(wěn)定的地磁場(chǎng)構(gòu)成。圖5e為地磁總場(chǎng)的希爾伯特譜分布，在整個(gè)時(shí)間段，噪聲和地磁場(chǎng)信號(hào)混疊在一起。

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

2.4 形態(tài)學(xué)-HHT濾波后的地磁總場(chǎng)

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

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

圖5 各個(gè)IMF分量及總場(chǎng)信號(hào)頻譜圖Fig.5 Spectrum of each IMF and total field signal注：(a)代表IMF1和IMF2分量疊加的希爾伯特譜；(b)代表IMF3和IMF4分量疊加的希爾伯特譜；(c)和(d)分別為IMF5、IMF6對(duì)應(yīng)的希爾伯特譜；(e)為總場(chǎng)信號(hào)的希爾伯特譜；(f)為總場(chǎng)信號(hào)的傅里葉譜。圖a和b中紅色細(xì)線(xiàn)指示分量信號(hào)最低頻率界限。各分量圖下方顏色條帶代表振幅能量，單位：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 兩種濾波處理后的總場(chǎng)數(shù)據(jù)Fig.6 The total field signal data after filtering by two different methods

2.5 濾波效果的真實(shí)性檢驗(yàn)

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

圖7 兩種濾波處理后的船磁補(bǔ)償結(jié)果對(duì)比Fig.7 The comparison of ship magnetic compensation by two different filtering methods

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

3 結(jié)論

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