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

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

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

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

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

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

1 CEEMDAN方法

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

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

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

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

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

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

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

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

分解得到k階殘差為:

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

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

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

2 CEEMDAN閾值濾波去噪模型

2.1 磁場誤差模型的建立

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

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

感應干擾磁場是指軟磁性材料在地磁場的作用下產生的磁場，其磁性隨外界磁場的變化而變化，呈現一定的比例關系，可表達為:

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

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

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

根據地磁場理論，在某一空間內地磁場的磁場強度矢量M的標量大小為恒定值，設為M0，所以當不考慮噪聲影響時，彈體在這一空間內做任意姿態運動時，雙軸磁強計所測得的磁場分量滿足:

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

2.2 利用自相關系數和傅里葉變換確定分界點

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

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

圖1 隨機信號與其歸一化自相關函數

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

圖2 正弦信號與其歸一化自相關函數

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

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

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

2.3 小波閾值去噪算法

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

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

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

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

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

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

2.4 去噪效果評價指標

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

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

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

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

3 實例分析

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

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

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

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

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

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

圖5 IMF5~IMF11與余項的自相關系數

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

圖6 IMF5~IMF8階的頻譜分析結果

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

得到重構信號和真實信號如圖7所示。

圖7 原始地磁信號與重構地磁信號

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

表1 去噪精度指標

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

圖8 不同信號構成斜橢圓

磁場誤差補償的效果可以用重構后的信號構成的橢圓距真實磁場數據構成的橢圓的幾何距離L來衡量

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

表2 幾何距離指標

4 結論

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