基于Hausdorff距離的光纖陀螺閾值濾波算法

曹 健 楊志強 劉晨晨

（長安大學地質工程與測繪學院，陜西西安710054）

采礦機械的導航定位與姿態控制是實現采礦自動化的關鍵技術之一，慣性導航以導航信息的實時性與自主性使其成為一種在復雜采礦工作面可行的采礦機械定位與姿態控制技術。光纖陀螺是一種基于Sagnac效應的全固態角速度傳感器，是慣性導航系統（Inertial navigation system，INS）的核心器件[1]，但其隨機漂移是制約INS精度的重要因素。為了抑制光纖陀螺的隨機漂移，傳統的濾波方法主要有兩種，一種是進行時間序列分析，建立自回歸滑動平均模型（Auto-Regressive Moving Average Model，ARMA），并設計Kalman濾波器進行補償[2]。而Kalman濾波需要完整的信號和噪聲的統計特性，其初始參數確認較復雜，同時在實際情況中，光纖陀螺會受到測試系統噪聲及外界環境等各種因素干擾，使輸出信號表現出非平穩、弱非線性的特征，因此難以建立準確的數學模型。第二種是使用小波變換的降噪方法[3-5]。小波具有良好的時頻局部性和多分辨分析能力，但小波消噪性能非常依賴于基函數的選擇和分解層數的預先設定，這些先驗信息的確定使去噪方案失去了自適應性，無法保證時變信號小波去噪的效果。

Huang等[6]針對非平穩、非線性信號提出了一種新型的信號處理方法——經驗模態分解。EMD是一種不需要任何信號先驗信息的自適應數據處理方法，能夠根據信號的本身特征，自適應地將信號分解為若干個由高頻到低頻的本征模態函數及一個余項。

在利用EMD濾波時，很多文獻都是通過經驗判斷或者多次測試，直接把分解后的前若干個IMF分量當作噪聲進行剔除，這類方法過于主觀，缺乏嚴謹的理論框架支撐，且很難恰當地排除噪聲干擾。將EMD與時間序列分析相結合，可以從全頻率角度對隨機誤差進行分析建模[7]，但該方法需要對所有IMF分別進行建模與去噪處理，過程繁瑣，計算量較大。當采用局部重構方法時，需要對提取的IMF進行物理解釋，以確定IMF是純噪聲、純信號，還是兩者都包含[8]。基于相關系數的EMD部分重構方法（EMD-Correlation Coefficient，EMD-Cor），使用相關系數來區分相關和不相關的IMF[9]。對于不同信噪比的噪聲信號，在某些情況下，由于各IMF與原始信號的相關性太強或太弱，使得該方法非常不穩定。基于連續均方誤差準則（EMD-Consecutive Mean Square Error，EMD-CMSE）的濾波方法[10]對IMF在能量第一次發生重大變化處開始進行部分重構，取得了不錯的效果，但是在某些情況下EMD-CMSE準則可能會被限制在局部極小值中。Ali Komaty等[11]提出表示數據分布形狀的概率密度函數，可以反映兩種信號的差異。根據輸入信號與各IMF的PDF之間的相似性來選擇相關IMF。這種相似性值越大，說明所比較的信號越具有相似的特征，因此重構信號中應包含相關的IMF。而基于信息理論度量（Kullback-Leibler Divergence，KLD）的概率相似性測量方法對于分布的形式（高斯、均勻、雙峰等）非常敏感，在某些情況下可能會造成誤判。

基于以上問題，本研究提出了一種基于幾何度量的概率相似性測量方法——Hausdorff距離[12]與硬閾值相結合的EMD信號濾波算法。

1 濾波方法原理

EMD將任意隨時間t變化的信號x()t分解為一組IMF，每個IMF具有不同的頻率尺度[6]。信號x(t)可分解為L個從高頻到低頻排列的IMF，以及一個余項rL(t)，即：

式中，IMFi(t)為第i階IMF分量。

根據經驗模態分解原理，對于混有隨機噪聲的信號，低階IMF分量通常對應于信號的高頻噪聲[13]。若除去若干個低階IMF，把剩余的相關IMF與余項重構為一個新信號，可以削弱隨機噪聲的影響。相關IMF的選擇需要根據一個給定的準則，該準則能夠識別攜帶了輸入信號信息的IMF分量。

假設y(t)為無噪的純信號，n(t)為高斯白噪聲，則光纖陀螺輸出信號x(t)模型可表示為

去噪的目標是找到x(t)的估計y?(t)，來提高信號的信噪比，即：

式中，kth為局部重構起始分量的階數，使用該式的關鍵是如何確定參數kth。

由于PDF包含其對應IMF的完整信息，因此可以使用PDF相似性量度來識別帶有y(t)特征的IMF。對原始信號進行EMD分解之后，使用核密度函數估計IMF的PDF。所以就要尋找這樣一種度量方法，利用PDF來區分相關和不相關的IMF，也就確定了參數kth。

1.1 Hausdorff距離

Hausdorff距離是描述兩個點集之間距離的方法，該方法可作為兩個點集之間的相似性量度。在本研究中，HD不是用作兩個幾何形狀間的相似性量度，而是用在兩個一維PDF之間。HD相較于其他距離度量方法的優點是對異常值非常敏感，同時顧及目標的整體形狀（或邊界），而不是將兩個空間目標的距離簡單地表達為兩個點之間的距離。

令d(a,b)=‖a-b‖為元素a、b之間的歐氏距離。元素a到有限集合B={b1,…,bN}之間的距離定義為

兩個有限集A={a1,…,aM}和B={b1,…,bN} 之間的單向HD可以定義為

類似地，兩個有限集B={b1,…,bN}和A={a1,…,aM}之間的單向HD可定義為

單向HD具有非對稱性，這是極大極小函數的特性。絕對HD為2個單向HD的最大值，即：

如果A、B兩組點是相似的（距離較近），但B中存在一點b，有遠離A的任意一點，那么該距離就是由這個點決定的。這一距離可以表明IMF的分布曲線有多少尖銳和突出的部分。

1.2 IMF的選擇及閾值濾波

為了選擇相關的IMF，首先將信號x(t)分解為若干個IMF，之后估計x(t)及各IMF的PDF，分別記為PDF[x(t)]與PDF[IMFi(t)]。將估計的分布看作度量空間（歐式空間）的子集，用HD測量相似程度。采用（5）式、（6）式和（7）式計算x(t)與各IMF之間PDF的HD值，記為HD(i):

重構所選擇的起始IMF分量是HD在第一個局部最大值之后開始減小的第一個IMF。記這個指數為kHD，即：

得到kHD后，則可確定信號的局部直接重構信號y1(t)：

若直接重構，則重構信號中仍然含有部分噪聲。所以在重構之前需要對相關IMF進一步去噪，方法可以借鑒小波閾值去噪。閾值去噪方法有硬閾值和軟閾值[14]，兩種閾值后的各階IMF分別見式（11）、式（12）。

式中，kHD≤i≤L；Ti表示第i個IMF的閾值，Ti的選擇可采用小波閾值估計方法[15-16]。

軟閾值的優點是具有良好的連續性和平滑性，但它是對所有系數進行壓縮，如果閾值估計產生誤差，該誤差會擴散到所有系數上，導致有偏的結果。而硬閾值只對低于閾值的部分進行處理，對高于閾值的部分予以保留，可以免受閾值估計誤差的影響，故而本研究采用硬閾值濾波方法。

由于使用了硬閾值的消噪策略，隨著i的增加，IMF系數的大部分幅值都高于閾值，閾值的作用會逐漸降低。所以無需判別噪聲信號混合的IMF與有用信號為主的高階IMF的邊界，不加區分地對所有需要被重構的模態進行閾值濾波。所以最終的重構信號為

綜上所述，本研究提出的基于Hausdorff距離的EMD閾值濾波（EMD-Hausdorff Distance&Threshold Filtering，EMD-HD&TF）算法的具體實現步驟如下：

（1）對原始信號進行EMD分解，得到若干個IMF分量。

（2）使用核密度函數，估計原始信號及各IMF的PDF。

（3）使用HD方法比較IMF分量與原始信號PDF的相似性，確定重構的起始分量——第kHD個IMF。

（4）計算各需要重構的IMF的噪聲閾值Ti，按照式（13）重構輸出信號。

2 試驗與結果分析

2.1 數據仿真試驗

為驗證硬閾值的優越性，對測試信號（以Doppler信號為例）加入不同大小的白噪聲，輸入信噪比為-6～6 dB，分別使用硬閾值與軟閾值對所有IMF進行濾波并重構，結果如圖1。硬閾值的去噪效果在各信噪比下都優于軟閾值，故本研究選擇HD篩選準則與硬閾值濾波相結合。

為測試本研究濾波方法的有效性，對4個代表性的測試信號（Doppler、Heavysine、Bumps和Blocks）分別加入以上噪聲，使用EMD-Cor、EMD-CMSE、EMDHD（EMD-Hausdorff Distance）及 EMD-HD&TF 4種方法進行去噪，對去噪結果進行比較研究。

以輸入信噪比為-4 dB的Doppler含噪信號為例，其經過EMD分解后得到的10個IMF見圖2。各階IMF的HD取值見圖3，由此可得kHD=6。在EMDHD&TF方法中，對IMF6～10分別進行閾值濾波，結果如圖4所示。

最終按式（13）將濾波后的IMF與余項重構。從圖5中可以看出，相比于其他幾種方法，由EMDHD&TF方法得到的信號更能保留Doppler標準信號的特征，初步判斷為最佳的去噪結果。表1給了的4種濾波結果的信噪比（Signal to Noise Ratio，SNR）及均方根誤差（Root Mean Square Error，RMSE）。

從表1可知，4種方法都對白噪聲具有一定的去噪能力，但EMD-Cor方法效果相對較差。兩種基于HD準則的EMD局部重構方法的去噪性能明顯優于另外兩種方法，在此基礎上使用硬閾值方法進行濾波可以進一步提升去噪效果。

統計4種信號在不同輸入信噪比下濾波后的信噪比及均方根誤差，如圖6與圖7所示。由此可知，基于EMD的濾波方法對于不同類型的信號都表現出不同程度的去噪效果，具有一定的適用性。其中EMD-Cor方法在輸入信噪比較低時性能較差，但會隨著輸入信噪比的提高而改善，逐漸達到與其他方法相近的水平。但在多次試驗中，經常會有重構后的信號與原信號相似度過高或者過低的現象發生，導致無法進行篩選，所以該方法穩定性較差。EMD-HD&TF在EMD-HD基礎上加入硬閾值，雖然濾波效果在本研究中只是略有提升，但在絕大多數情況下都能得到最佳的濾波效果。上述去噪效果不僅支持了EMD-HD和閾值濾波混合方法的有效性，而且證明了IMF中所攜帶的信號信息可以通過它們的PDF來表現，同時HD比EMD-Cor和EMD-CMSE法更有效地量化了信號與相關IMF之間的相互作用，進而使用HD進行篩選更為有效。

2.2 陀螺濾波試驗

在常溫（25℃）下，對某型號陀螺漂移信號進行了試驗測試。采樣周期為100 ms，為了避免啟動過程中的溫度變化，選取中間3 h的數據，樣本總數為108 000個。對該組數據分別采用文中的4種方法進行處理。其中EMD-Cor、EMD-CMSE、與兩種HD方法重構的起始IMF分量分別為第5階、第7階和第9階。為了驗證降噪效果，采用了慣性器件誤差分析中最常用的艾倫（Allan）方差。Allan方差的雙對數曲線圖可以識別陀螺的各種隨機誤差成分，包括量化噪聲（Q）、角度隨機游走（N）、零偏不穩定性（B）、角速率隨機游走（K）和速率斜坡（R），由試驗所得的Allan方差圖如圖8所示。不同的誤差具有不同的斜率特性，可由Allan方差擬合得到[17-18]，擬合得到的各種誤差成分如表2所示。

由圖8和表2可知，4種方法對角速率隨機游走（K）和速率斜坡（R）均有良好的抑制效果。但經EMD-Cor和EMD-CMSE方法處理后，量化噪聲（Q），角度隨機游走（N）和零偏不穩定性（B）不但沒有減少，甚至還有小幅度增加。基于HD方法得到的重構IMF較為恰當，其抑制5種噪聲的效果優于前兩者。在此基礎上加入閾值濾波的EMD-HD&TF方法又能進一步濾除噪聲，使各項噪聲得到了進一步降低，取得了最佳濾波效果。

?

3 結 論

（1）本研究提出了一種新的EMD局部重構與閾值相結合的EMD-HD&TF方法，該方法首先以各階IMF和原始信號PDF的HD為依據對IMF進行篩選，減少了IMF的重構個數；然后采用硬閾值處理方法，對信號的噪聲進行了進一步濾除。該方法的優點是對異常值敏感、復雜度低，在濾波應用中較為實用。

（2）為了驗證該方法的有效性，在不同的信噪比情況下進行了數值模擬。與EMD-Cor、EMD-CMSE和單純的EMD-HD方法相比，該方法具有更好的去噪性能。最后，將該方法應用于對實測陀螺漂移信號進行試驗。從Allan方差曲線可以看出，EMD-HD&TF方法大大降低了各項隨機誤差。因此，本研究提出的光纖陀螺消噪算法在采礦機械導航定位與姿態控制的自動化方面具有廣泛的應用前景。

?