基于CEEMDAN閾值濾波的磁場信號去噪模型*

董曉芬，陳國光，田曉麗，閆小龍，趙歸平

(中北大學(xué)機電工程學(xué)院，山西 太原 030051)

地磁場具有全天候、全天時、全區(qū)域的特點，磁強計具有抗干擾能力強、無積累誤差、成本低、抗沖擊機性強的優(yōu)點，使得地磁導(dǎo)航姿態(tài)測量技術(shù)得到了廣泛的應(yīng)用。但是地磁信號非常微弱，安裝在載體內(nèi)部的磁強計容易受到固定磁場、感應(yīng)磁場和隨機磁場等干擾，從而導(dǎo)致載體姿態(tài)測量精度大幅度下降，因此需要對各種干擾磁場進行補償[1-2]。

目前，國內(nèi)外學(xué)者研究磁場信號去噪的方法主要有傅里葉變換、小波分析、獨立分量法以及基于信號時域分析的形態(tài)濾波法等，單志超根據(jù)鐵磁性目標(biāo)信號的噪聲特性和時頻特性，提出了一種基于小波變換的磁異常信號背景噪聲消除方法[3]。Mandrikova采用小波閾值法和平穩(wěn)小波變換抑制噪聲，提高了地磁異常信號的分辨率[4]。李季采用形態(tài)濾波器對磁場信號進行處理，有效地消除了信號中含有的強脈沖干擾[5]。除了這些方法，經(jīng)驗?zāi)B(tài)分解EMD(Empirical Mode Decomposition)和Hilbert-Huang變換在信號處理領(lǐng)域應(yīng)用較為廣泛，它是一種對非線性和非平穩(wěn)信號自適應(yīng)分解的處理方法和時頻幅分析方法，但是使用這種方法會出現(xiàn)模態(tài)混疊、端點效應(yīng)等問題，為了改進這些問題，很多學(xué)者提出了許多改進算法，如集合經(jīng)驗?zāi)B(tài)分解(EEMD)、完備總體經(jīng)驗?zāi)B(tài)分解(CEEMD)和完整集成經(jīng)驗?zāi)B(tài)分解(CEEDMAN)[6-7]。

本文在分析磁場誤差模型的基礎(chǔ)上，借鑒EMD閾值設(shè)置方式和小波軟閾值方法，將CEEMDAN閾值濾波方法應(yīng)用于磁場信號的去噪中，對含噪聲的地磁信號進行完備分解，并將經(jīng)過閾值濾波處理后的信號與有效信號重構(gòu)，得到去噪后的信號，通過對實際磁場信號進行分析驗證了該方法的有效性。

1 CEEMDAN方法

CEEMDAN根據(jù)信號特性自適應(yīng)地向其添加不同的高斯白噪聲，通過計算的殘余信號來獲取下一級的IMF分量，將含噪信號從高頻到低頻依次分解出一組IMF分量，可有效解決使用EMD出現(xiàn)的模態(tài)混疊問題并提升分解過程的完備性及分解效率[8]。

首先定義一個運算符Ej(·)，表示為采用EMD產(chǎn)生的第j個IMF，ωi(n)為N(0，1)的白噪聲，εi為白噪聲的幅值系數(shù)，x(n)為原始數(shù)據(jù)，CEEMDAN算法的步驟如下:

①首先向原始信號x(n)中添加白噪聲ε0ωi(n)，獲得信號Xi(n)為:

式中:i為加入高斯白噪聲的次數(shù)i＝1，2，…，N。先求取第一階IMF分量IMF1，對于第i次加入白噪聲的Xi(n)，進行EMD分解，得到IMFi1，重復(fù)計算N次，對其集成平均可得到第一階IMF為:

由原始信號x(n)和第一階IMF獲取一階殘差為:

②在一階殘差r1(n)的基礎(chǔ)上添加白噪聲ε1ωi(n)，獲得信號r1(n)＋ε1ωi(n)，同樣進行EMD分解，得到IMFi2，重復(fù)計算N次，對其集成平均可得到第二階IMF為:

由一階殘差r1(n)和第二階IMF獲取二階殘差為:

③在殘差rk-1(n)的基礎(chǔ)上添加白噪聲εk-1ωi(n)，獲得信號rk-1(n)＋εk-1ωi(n)，進行EMD分解，得到IMFik，重復(fù)計算N次，則k階IMF為:

分解得到k階殘差為:

④返回步驟③繼續(xù)進行下一步分解，直到殘差不能分解為止，最終殘差為R(n)；

⑤獲得最終分解結(jié)果，原始信號x(n)可以表示為:

最終信號通過CEEMDAN分解后得到一系列從高頻到低頻的IMF分量，CEEMDAN利用噪聲自適應(yīng)的特點，可對原始信號進行精確完整的重構(gòu)。

2 CEEMDAN閾值濾波去噪模型

2.1 磁場誤差模型的建立

使用雙軸磁強計測量彈體滾轉(zhuǎn)角過程中，有很多磁場干擾源，根據(jù)其表現(xiàn)形式的不同可以分為固定干擾磁場、感應(yīng)干擾磁場、隨機干擾磁場等[9]。

固定干擾磁場是指彈體上硬磁性材料產(chǎn)生的磁場，其磁性在一定時間內(nèi)不隨外界磁場的變化而變化，因此固定磁場的磁場強度不變，可表達為:

感應(yīng)干擾磁場是指軟磁性材料在地磁場的作用下產(chǎn)生的磁場，其磁性隨外界磁場的變化而變化，呈現(xiàn)一定的比例關(guān)系，可表達為:

磁強計的制造誤差包括零位誤差、靈敏度誤差和非正交誤差，其系數(shù)分別為零位誤差系數(shù)Bo、靈敏度誤差系數(shù)Cl和非正交誤差系數(shù)Cn，大小不隨磁場的變化而變化。

隨機磁場的來源復(fù)雜，頻譜范圍寬，主要包括電路系統(tǒng)內(nèi)電子元器件內(nèi)部產(chǎn)生的熱噪聲、散粒噪聲和量子噪聲，可認(rèn)為其為白噪聲，還有舵機產(chǎn)生的影響，舵機由電機驅(qū)動，在工作的瞬間會產(chǎn)生瞬時強磁場脈沖，工作過程中，舵翼轉(zhuǎn)動的頻率大小也會產(chǎn)生不同頻率的噪聲，隨機干擾磁場用ωt表示。

設(shè)真實磁場信號為M，則磁強計輸出的誤差模型為:

根據(jù)地磁場理論，在某一空間內(nèi)地磁場的磁場強度矢量M的標(biāo)量大小為恒定值，設(shè)為M0，所以當(dāng)不考慮噪聲影響時，彈體在這一空間內(nèi)做任意姿態(tài)運動時，雙軸磁強計所測得的磁場分量滿足:

從上式可以看出，不考慮噪聲影響時，隨著彈體的姿態(tài)角運動，矢量M頂點的軌跡是一個圓心在原點，半徑為M的地磁圓[10]。在固定磁場干擾、感應(yīng)磁場干擾、自身誤差和隨機誤差的影響下，由磁強計誤差模型Ms＝CrM＋Mr＋ωt可知，雙軸磁強計所測數(shù)據(jù)的軌跡在幾何上可解釋為經(jīng)過伸縮、旋轉(zhuǎn)、平移畸變?yōu)樾睓E圓，由于隨機誤差不能用固定的參數(shù)進行表示，導(dǎo)致斜橢圓模型的固定干擾系數(shù)、感應(yīng)干擾系數(shù)、零位誤差系數(shù)、靈敏度系數(shù)、非正交系數(shù)等不能準(zhǔn)確地估計，所以需要提前將磁場信號進行濾波處理，濾除隨機干擾，然后進行斜橢圓擬合參數(shù)估計，可得到斜橢圓模型如式(13)所示:

2.2 利用自相關(guān)系數(shù)和傅里葉變換確定分界點

采用CEEMDAN方法對含隨機噪聲的磁場信號進行完備分解后，得到主成分為噪聲的IMF、有效地磁信號與噪聲混疊的IMF和有效地磁信號IMF，為了將包含噪聲的IMF分量中的噪聲信號濾除，和剩下的有效信號IMF分量進行重構(gòu)來實現(xiàn)信號的降噪。因此需要準(zhǔn)確地找到噪聲IMF分量與有用信號的IMF分量之間的分界，它直接決定CEEMDAN濾波方法的降噪效果。

文中利用噪聲的數(shù)學(xué)特性來尋找分割點，自相關(guān)函數(shù)是描述信號變化的劇烈程度的函數(shù)，也反映了同一個信號在不同時刻間的相關(guān)程度。由圖1可看出，對于理想高斯白噪聲信號，由于數(shù)據(jù)在各個時間點的隨機性，在不同時刻兩個信號之間的關(guān)聯(lián)性很差，所以隨機噪聲的自相關(guān)系數(shù)在零處取得最大值，在其他點處趨近于零[11]。

圖1 隨機信號與其歸一化自相關(guān)函數(shù)

分析圖2的正弦函數(shù)歸一化自相關(guān)函數(shù)，可觀察到自相關(guān)系數(shù)也是在零處取得最大值，但是由于正弦信號在不同時刻的數(shù)據(jù)存在著相關(guān)性，因此在其他點處自相關(guān)系數(shù)是隨著時間差的變化而變化的，沒有立即減小到零。

圖2 正弦信號與其歸一化自相關(guān)函數(shù)

已知原始信號為x(t)，其自相關(guān)函數(shù)表達式為(14):

采用歸一化自相關(guān)函數(shù)計算噪聲與信號模態(tài)分量的自相關(guān)函數(shù)值，歸一化函數(shù)的表達式為式(15):

對含有隨機噪聲的地磁信號進行CEEMDAN分解得到IMF分量，可以根據(jù)上述特點來根系各個IMF的自相關(guān)函數(shù)，進而判斷出含噪聲的IMF分量和有效信號IMF之間的界限[11]。但是由于噪聲的復(fù)雜性，自相關(guān)函數(shù)曲線的不確定性，可以對各個IMF分量進行傅里葉分析，從頻域?qū)π盘栠M行分析，根據(jù)地磁信號的低頻特性來進一步確定分界點。

2.3 小波閾值去噪算法

由于一般信號在時間域上的連續(xù)性，而隨機噪聲在時間域上是不連續(xù)性的，因此在小波域，有效信號所產(chǎn)生的小波系數(shù)其模值較大；而高斯白噪聲經(jīng)過小波變換，在小波閾仍然表現(xiàn)為很強的隨機性，其對應(yīng)的系數(shù)很小。根據(jù)該特點，需要通過設(shè)置閾值函數(shù)進行去噪處理，大于閾值的系數(shù)保留下來，而小于閾值的分解系數(shù)則通過置零進行消除[12]。

CEEMDAN通過添加有限次的自適應(yīng)白噪聲來減少重構(gòu)誤差，但是經(jīng)過分解得到的高頻IMF分量可能中含有有效分量信號，因此本文根據(jù)EMD閾值的設(shè)置方法，采用小波軟閾值的處理方式，將高頻IMF進行閾值處理后，對其濾波后的分量與有用分量進行重構(gòu)以達到信號去噪的目的[13]。

對前(k-1)個IMF噪聲分量使用閾值濾波，閾值處理形式采用軟閾值的方法，表達式為:

式中:imf′i為去噪后的第i階IMF分量，Ti為第i階IMF分量的閾值。閾值確定的方法采用EMD方法中的閾值計算模型:

式中:σ唯一常數(shù)，Ei為第i階IMF對應(yīng)的能量，其值可以使用式(18)進行估計:

對前k階IMF進行小波軟閾值處理之后，得到前k階濾波后的IMF分量為IMF′，和k階后的IMF分量可得到重構(gòu)后的信號為:

2.4 去噪效果評價指標(biāo)

為了分析CEEMDAN閾值濾波算法的降噪效果，本文采用均方根誤差和信噪比指標(biāo)對地磁信號的去噪效果進行評價。

均方根誤差表示的是真實信號與重構(gòu)信號偏差大小的一個指標(biāo)，均方根誤差越小，表示去噪后的信號與真實信號之間的偏差越小[14]，其表達式為:

信噪比為有用信號功率與噪聲功率的比值，表示的是信號中的有效成分與噪聲成分的比例關(guān)系，信噪比越大，表示信號中噪聲所占的比例越小，其表達式為:

式中:PS為有用信號的功率，PN為噪聲的功率。

3 實例分析

實驗中采用雙軸磁強計，彈體上的電子設(shè)備及供電系統(tǒng)和執(zhí)行機構(gòu)工作的時候會產(chǎn)生預(yù)期范圍內(nèi)的高頻干擾噪聲。為了驗證CEEMDAN濾波算法在地磁測量信號中的濾波效果，首先，只給地磁組件供電，將裝有磁強計的彈體放在無磁轉(zhuǎn)臺上，以1.5 Hz的頻率勻速轉(zhuǎn)動，磁強計輸出的模擬信號經(jīng)過16AD轉(zhuǎn)換后經(jīng)過DSP以50 Hz的采樣頻率進行采集，得到無外界干擾的磁場信號，然后地磁組件和其他組件一起供電，在旋轉(zhuǎn)的過程中電機以10 Hz的頻率開始工作，可得到包含噪聲的磁場信號，以采集到的Y軸的信號Ny為例，如圖3所示。

圖3 原始磁場信號與含隨機噪聲的磁場信號

文中將采集到的無外界干擾的磁場信號作為真值，使用CEEMDAN閾值濾波去噪方法對包含噪聲的磁場信號進行處理后與真值進行比較。

首先使用CEEMDAN方法對采集到的磁場信號Ny進行分解，其中參數(shù)設(shè)置如下:噪聲標(biāo)準(zhǔn)差為0.2，集成平均次數(shù)為500，允許的最大迭代次數(shù)為500。Ny經(jīng)過分解后得到11個IMF分量和一個殘余分量，其中分量IMF5~IMF11和余項12如圖4所示，通過圖4可以觀察到經(jīng)分解后的IMF分量從第八階IMF分量開始幅值的數(shù)量級發(fā)生了變化。

圖4 含隨機噪聲的磁場信號經(jīng)CEEMDAN分解得到IMF5~IMF11與余項

然后根據(jù)Ny分解得到的各階IMF的自相關(guān)系數(shù)來判斷含噪聲的IMF分量和有效信號IMF分量之間的分界點，其中第1階IMF到12階IMF的自相關(guān)系數(shù)如圖5所示:

圖5 IMF5~IMF11與余項的自相關(guān)系數(shù)

從圖5中可以看出Ny的IMF分量的前七階IMF的自相關(guān)系數(shù)和高斯白噪聲信號的自相關(guān)系數(shù)的特點類似，從第八階開始自相關(guān)函數(shù)的曲線的特點為零點處最大，在其他點隨著時間差的變化而變化，和一般信號的自相關(guān)函數(shù)相似，由此可以估計分界點為第八階。再觀察第5階到第八階各階IMF的傅里葉變換，如圖6所示，從圖5中可以觀察到第八階IMF分量的頻率集中在1.5 Hz，和轉(zhuǎn)臺轉(zhuǎn)動的頻率接近，根據(jù)地磁信號屬于低頻信號，其頻率值一般在1 Hz以下，由此可以確定含有噪聲的IMF分量和有效信號IMF分量之間的分界點為第八階IMF。

圖6 IMF5~IMF8階的頻譜分析結(jié)果

然后分別對磁場信號Ny的前八階IMF分量進行小波軟閾值處理得到IMF′，根據(jù)式(22)

得到重構(gòu)信號和真實信號如圖7所示。

圖7 原始地磁信號與重構(gòu)地磁信號

最后為了驗證該算法的可行性，本文分別采用了使用CEEMDAN強制濾波和本文提出的CEEMDAN閾值濾波算法處理磁場信號進行比較，通過表1可以發(fā)現(xiàn)分別信噪比提高了1.156 dB，均方根誤差提升了10.52%，如表1所示。

表1 去噪精度指標(biāo)

將經(jīng)過CEEMDAN小波閾值濾波得到的重構(gòu)信號Ny、Nz根據(jù)磁場誤差模型，建立斜橢圓模型，得到由真實數(shù)據(jù)、含噪數(shù)據(jù)和重構(gòu)數(shù)據(jù)得到的斜橢圓分別如圖8所示:

圖8 不同信號構(gòu)成斜橢圓

磁場誤差補償?shù)男Ч梢杂弥貥?gòu)后的信號構(gòu)成的橢圓距真實磁場數(shù)據(jù)構(gòu)成的橢圓的幾何距離L來衡量

幾何距離L越小，說明補償后的斜橢圓越接近真實斜橢圓，可得到含噪聲的磁場數(shù)據(jù)、經(jīng)過CEEMDAN強制濾波后的磁場數(shù)據(jù)和經(jīng)過CEEMDAN閾值濾波處理后的磁場數(shù)據(jù)的幾何距離，如表2所示，對比CEEMDAN強制濾波，經(jīng)過CEEMDAN閾值濾波處理后的磁場信號幾何距離下降了13.87%，表明了本文所提方法的有效性，可以為橢圓參數(shù)估計進行預(yù)處理。

表2 幾何距離指標(biāo)

4 結(jié)論

本文針對地磁信號的特點，根據(jù)CEEMDAN方法能夠?qū)π盘栠M行完備分解的基礎(chǔ)上，提出使用自相關(guān)系數(shù)和傅里葉變換確定IMF分量中噪聲IMF分量與有用信號IMF分量之間的界限，借鑒EMD閾值設(shè)置方式，結(jié)合小波軟閾值對含噪聲的IMF分量進行濾波處理，最后重構(gòu)信號以得到去噪后的信號。通過對實際磁強計輸出磁信號進行CEEMDAN閾值濾波處理之后分析均方根誤差、信噪比和橢圓間的幾何距離，結(jié)果表明:本文提出的方法要優(yōu)于CEEMDAN強制去噪方法，驗證了該方法在處理含噪聲的磁場信號上的可行性和優(yōu)越性。