基于經驗模分解的陀螺信號消噪

甘 雨，隋立芬

信息工程大學測繪學院，河南鄭州450052

基于經驗模分解的陀螺信號消噪

甘 雨，隋立芬

信息工程大學測繪學院，河南鄭州450052

陀螺隨機漂移是影響慣性導航精度的重要因素。小波消噪方法對異常噪聲效果不明顯，且對小波基和分解尺度等因素依賴性較強。提出陀螺信號經驗模分解（EMD）消噪方法，將信號進行經驗模分解得到一個本征模態(tài)函數（IMF）組，先基于2σ準則處理異常噪聲IMF分量，再利用相關系數確定高頻噪聲IMF分量個數，將噪聲分量去除以實現(xiàn)陀螺信號消噪。詳細對比小波方法與EMD方法，利用交疊式Allan方差分析兩者的消噪效果，通過慣導算例進一步驗證EMD方法的實效性。結果表明，相比小波方法，EMD消噪法能剔除異常噪聲，可以更有效地抑制陀螺漂移。

陀螺隨機漂移；小波；經驗模分解；消噪

1 引 言

陀螺隨機漂移是影響慣性導航精度的重要因素。抑制陀螺漂移的方法主要有兩種［1］：① 建立漂移模型，使用Kalman濾波等方法進行補償；②對陀螺輸出信號進行消噪處理。由于隨機漂移往往表現(xiàn)為弱非線性、非平穩(wěn)、慢時變，且易受到外部環(huán)境等多種不確定因素的影響［2］，無法建立其準確的系統(tǒng)模型，故需要采用陀螺信號消噪的方法。

目前對陀螺信號消噪主要采用小波方法。小波具有優(yōu)良的多分辨率分析特性，小波消噪不需要系統(tǒng)的誤差模型，因此被廣泛用于陀螺信號的消噪處理中［1－4］。然而，小波分解雖然也能實現(xiàn)對非平穩(wěn)信號的濾波，但其實質是帶通濾波器，限制了濾波的精確性，且小波變換中小波基一經選定，整個信號分析過程中就只能使用這一個小波基，即小波變換是非適應性的［5］。有研究表明［6］：對于如陀螺信號一類的非平穩(wěn)信號，當異常噪聲淹沒了有用信號時，采用小波消噪效果也不甚理想，對這類信號的消噪目前還未發(fā)現(xiàn)較好的方法。

文獻［7］提出一種分析非平穩(wěn)、非線性信號的自適應分解方法：經驗模分解方法（empirical mode decomposition，EMD），它將復雜的信號分解成若干個按頻率高低排列的本征模態(tài)函數（intrinsic mode function，IMF），每個IMF是一個單分量信號。該方法與小波分析的區(qū)別在于它不需要事先選定基函數，而是根據信號本身的特性自適應地產生合適的模態(tài)函數，這些模態(tài)函數能很好地反映信號在任何時間局部的頻率特征。該方法已經用于機械振動信號分析、SAR影像濾波、氣象以及GPS信號處理等領域［8－11］。

將EMD方法引入到陀螺信號的處理中，給出陀螺EMD消噪的方法，按照2σ準則剔除振幅偏大的異常噪聲，利用相關系數確定隨機高頻噪聲IMF分量的個數。與小波消噪依賴小波基、分解尺度、閾值估計方法不同，EMD消噪過程完全依賴陀螺信號本身特性。用交疊式Allan方差對比分析本文方法與小波消噪方法，通過慣性導航算例進一步驗證方法的實效性。

2 EMD消噪方法

2.1 EMD基本原理

EMD分解方法認為任何待分解信號都由一組固有振動模式構成，并據此將信號分解為若干本征模態(tài)函數IMF的和。這些本征模態(tài)函數既可以是線性的，也可以是非線性的；既可以是平穩(wěn)的，也可以是非平穩(wěn)的［11］。分解得到的IMF分量滿足：①零點數目與極值點數目相同或至多相差1；②函數由局部極大值點構成的包絡線和由局部極小值構成的包絡線的均值為零。分解過程通過一個稱為“篩選”的步驟來完成。“篩選”過程可以表示為［7］：

（1）分別由原始信號x（t）的極小值點與極大值點用三次樣條插值得到x（t）的上下包絡線，計算上下包絡的均值m1，進而計算x（t）和m1的差值h1＝x（t）－m1，判斷h1是否滿足IMF的兩個條件。若滿足，則h1為x（t）的第一個分量imf1。

（2）若不滿足，則將h1作為新的信號繼續(xù)步驟（1），得到h11，判斷h11是否滿足條件。若不滿足則繼續(xù)（1）的步驟，直到重復k次后h1k滿足IMF的條件，則有imf1＝h1k，求出原始信號與imf1的差值r1＝x（t）－imf1。

（3）將r1作為新的“原始”信號重復上述步驟，直到提取出第2個，第3個，直至第n個IMF分量imf2，imf3，…，imfn，有r2＝r1－imf2，…，rn＝rn－1－imfn，當IMF分量imfn或余項rn小于預先設定的值，或者余項rn已經成為單調函數時，整個篩選分解過程結束，本文以余項為單調函數作為分解的終止條件。

經過上述步驟后，原始信號x（t）可分解為n個IMF分量和1個余項的和

從上述的經驗模態(tài)分解方法可以看出，越早分解出來的IMF頻率越高［7］，第一個分解出來的代表原信號的最高頻率成份，各個IMF的頻率幾乎是按2的負冪次方的形式遞減［6，12］。對于混有隨機噪聲的信號，其先分解出的IMF分量通常對應于信號的高頻噪聲［8］。若IMF組去除了先分解的幾個IMF，把其余的IMF組合起來形成一個信號，可以削弱信號的噪聲。可以看出，用EMD對陀螺信號消噪的關鍵在于從信號分解出的IMF組中辨識出噪聲IMF分量。

2.2 處理異常噪聲

由于外界環(huán)境的異常變化，陀螺信號會在局部時段內受到振幅較大的異常噪聲的干擾，導致信號輸出失真，對陀螺信號的處理首先要消除這類噪聲的影響，而小波對這種噪聲的消噪效果不佳。經驗模分解將信號按照其固有振動模式分成IMF組［13］，而異常噪聲的頻率和振動模式一般既不同于隨機高頻噪聲，也不同于真實信號，對信號進行分解之后，異常噪聲一般被分解到其中一個或幾個IMF之中，只要能將相應的IMF進行識別并加以處理，就能減弱異常噪聲的影響。

利用誤差理論中制定極限誤差的2σ準則，可以對陀螺信號的異常噪聲IMF分量進行識別。設原始信號為x（t），求得其采樣標準差為σ，計算IMF各個分量imfi的最大振幅Ai，對各分量按如下判斷準則進行處理

異常噪聲作用的時間相對比較短，σ通常反映陀螺信號的整體觀測精度，若某個imfi的最大振幅大于2σ，則可認為它是異常噪聲分量。IMF分量的頻率各不相同，且如前面所述按2的負冪次方遞減，因而異常噪聲通常會被分解到少數IMF分量且是中高頻分量，非常便于進行識別和剔除。

判別和處理異常噪聲是陀螺EMD消噪方法的第一步，除非陀螺的觀測環(huán)境非常好，否則都應該在對信號進行EMD分解后進行異常噪聲的辨識和剔除。

圖1給出某次試驗中陀螺（標稱陀螺漂移為1°／h）靜基座下X軸8 000個歷元的輸出信號，在1 500～2 000、4 100～4 300、5 600～6 900歷元處受到異常噪聲的干擾。對該信號進行EMD分解，得到10個IMF分量，圖2顯示為前6個IMF以及最終的余項r10。可直觀看出，1 500～2 000、4 100～4 300歷元的異常噪聲主要分解在imf3內，5 600～6 900歷元的分量分解到imf2和imf3內。

計算陀螺信號標準差為σ＝34.02，則限差2σ＝68.04，各IMF分量的最大振幅分別為A1＝10.66，A2＝216.38，A3＝87.51，A4＝39.85，A5＝42.96，A6＝12.28。按照式（2），imf2和 imf3的振幅超限，應該將這兩個分量置零，這與圖形所示的結果一致。

圖1 陀螺原始信號Fig.1 Original gyro signal

圖2 陀螺信號EMD分解的IMF分量Fig.2 IMFs from EMD of gyro signal

2.3 確定噪聲IMF分量個數

剔除異常噪聲后，需要消除陀螺信號中的高頻噪聲。設前m級IMF分量（已經過異常噪聲的處理）為高頻噪聲分量，則消噪過程可表示為

若選取的m偏小，濾波后的信號x′（t）仍會受到噪聲干擾，若選取的m偏大，可能會把有用的信號濾去。為避免此問題，有文獻借鑒小波閾值消噪的思想提出對每個IMF分量也進行閾值處理［14］。然而，每個IMF分量都是具有特定物理意義的信號成分，反映原始信號某一方面的內在特征［7，13］，利用閾值消噪的方法可能會破壞信號特征的完整性，且閾值消噪法存在確定閾值估計方法的難點，其實用性有待進一步研究。EMD消噪的關鍵還是在于確定噪聲分量的個數m。

長度為N的信號序列x（t）和y（t）的采樣相關系數可定義為

相關系數反映兩個信號之間的相似程度或依賴程度。設前m個IMF分量之和為若這m個分量都是噪聲IMF分量，不含原信號的有用成分，則按照式（3）濾波后的信號x′m（t）與原信號x（t）之間應該具有較強的依賴程度。因此可考慮由相關系數確定噪聲分量的個數，計算

從m＝1開始計算，當ρxx′m＞c（c為常量，可取經驗值0.75～0.8）時，m＝m＋1，計算下一個ρxx′m，直到某個ρxx′m＜c，停止計算，說明第m個分量imfm為信號的有用成分，不是噪聲IMF，取M＝m－1作為噪聲分量個數。最終的EMD消噪結果為

陀螺信號的EMD消噪方法歸納為：

（1）對陀螺信號進行EMD分解，得到IMF分量；

（2）按2σ準則處理異常噪聲IMF；

（3）根據相關系數確定噪聲IMF個數，原信號減去噪聲IMF完成消噪過程。

以上的EMD過程無需事先的參數設置，只與陀螺信號本身有關。陀螺信號的EMD分解是按照信號的固有振動模式進行的，采樣標準差σ和相關系數ρxx′m都由陀螺信號及其分解的IMF分量來計算，整個消噪過程不受外在因素影響。需要注意的是，常量c的取值及2σ準則的確定是EMD消噪中的重要環(huán)節(jié)，如果選取不合適可能會影響最終的消噪結果，因此實用中需要結合實際的數據類型進行分析確定。

確定圖2中噪聲分量個數，從m＝4開始，相關系數均小于0.7，則取噪聲分量個數為3，EMD消噪的結果如圖3所示。

圖3 EMD消噪信號Fig.3 EMD de－noised signal

3 消噪方法對比

3.1 消噪信號直接比較

為對比小波消噪法與本文EMD消噪法，對圖1中的陀螺原始信號進行小波閾值消噪，采用軟閾值函數，閾值估計方法為SUREShrink閾值［3，15］。分別取不同小波基與分解尺度進行消噪試驗，消噪的效果總體一致，圖4、圖5分別為db8小波9尺度消噪結果和Haar小波7尺度消噪結果。

圖4 db8小波消噪信號Fig.4 db8wavelet de－noised signal

圖5 Haar小波消噪信號Fig.5 Haar wavelet de－noised signal

比較圖1、圖3、圖4和圖5，可以看出：

（1）EMD利用異常噪聲與有用信號和高頻噪聲不同的內在振動模態(tài)將異常噪聲分離，因此EMD消噪法既能消去隨機高頻噪聲，又能剔除異常噪聲；

（2）小波消噪方法能有效地去除高頻噪聲，但是不能抵制異常噪聲的干擾；

（3）小波分析的本質決定了它不能分離出異常噪聲，通過選取不同的小波基和分解尺度也無法克服這一缺陷。

3.2 交疊式Allan方差對比分析

采用交疊式Allan方差進一步對比兩種方法的消噪效果。Allan方差適合于分析非平穩(wěn)的隨機信號，是IEEE推薦的陀螺隨機誤差分析方法［16］，其中的交疊式Allan方差在相同的置信水平下比普通Allan方差分析方法具有更大的置信區(qū)間，是對Allan方差的改進［17］，其具體定義見文獻［17—18］。

信號的Allan方差可表示為σ2（τ）（τ＝nτ0，τ0為采樣間隔），則Allan標準差為σ（τ）～τ雙對數曲線圖可以描述陀螺的各種隨機誤差成分，不同的成分具有不同的斜率特性，這些成分包括量化噪聲，角度隨機游走，零偏不穩(wěn)定性，速率隨機游走和速率斜坡，可由Allan方差擬合得到，具體計算方法和它們在σ（τ）～τ圖中的斜率特性可參考文獻［18］。

分別求取陀螺原始信號、小波消噪信號、EMD消噪信號的交疊式Allan標準差σ（τ），σ（τ）～τ雙對數曲線如圖6所示，其中db8和Haar小波的結果接近，只顯示db8小波結果。根據計算的Allan方差擬合得到的各種誤差成分的值如表1所示。

圖6 σ（τ）～τ曲線對比Fig.6 Comparison ofσ（τ）～τcurve

表1 陀螺誤差成分對比Tab.1 Comparison of error components of gyro

分析圖6和表1的結果可知：

（1）小波對陀螺信號具有一定的消噪作用，但是由于受到異常噪聲的干擾，小波消噪后的陀螺信號中仍然含有較大的誤差成分；

（2）EMD消噪方法有效地削弱了陀螺各種誤差成分；

（3）在τ＝0.2～0.3s處，EMD及小波消噪后信號的Allan標準差與原始信號相比都沒有明顯改善，是因為受到殘余低頻有色噪聲的影響。一般可通過建立時間序列模型的方法削弱陀螺信號的有色噪聲成分。

4 慣導解算與分析

所用數據為一組靜態(tài)慣性測量單元（inertial measurement unit，IMU）數據，標稱陀螺漂移和加速度計偏置分別為1°／h和10－4g，IMU采樣頻率100Hz，取歷元282 930.83～283 232.82s（GPST，對應的周數為1 467）的數據進行試驗。利用282 930.83～283 112.83s的數據進行Kalman濾波精對準得到初始姿態(tài)角，初始航向角、俯仰角和翻滾角分別為8.916°、－0.013°、1.094°，分別采取3種方案對283 112.83～283 232.82s的IMU數據進行慣導解算，得到各歷元的速度值。各方案的加速度計數據均為原始輸出，但陀螺數據的處理方法不同。由于靜止狀態(tài)下慣導的速度真值實際上是“零”，因此，計算的速度值其實就是速度誤差。采用如下3種方案進行解算：

方案1 使用IMU輸出的各軸陀螺原始數據解算；

方案2 使用db8小波消噪后的各軸陀螺數據進行解算；

方案3 使用本文EMD方法消噪后的各軸陀螺數據進行解算。

3種方案的速度誤差見圖7～圖9，東向和北向的誤差結果類似，這里給出的是東向的結果，RMS的比較如表2所示。

圖7 方案1東向速度誤差Fig.7 Veerror of scheme 1

圖8 方案2東向速度誤差Fig.8 Veerror of scheme 2

圖9 方案3東向速度誤差Fig.9 Veerror of scheme 3

表2 3種方案RMS比較Tab.2 Comparison of RMS for three schemes（m／s）

分析解算結果，可知：

（1）受到加速度計誤差和殘余陀螺漂移的影響，單獨慣導解算的誤差隨著歷元數增加而不斷積累，3種方案的誤差均呈現(xiàn)擴大趨勢，需要引入其他的導航信息（如衛(wèi)星導航）才能予以修正。

（2）方案2計算的速度誤差明顯小于方案1的結果，說明小波具有削弱陀螺漂移、減小慣導誤差的能力。

（3）用EMD方法不僅削弱了陀螺漂移中隨機高頻噪聲的影響，而且消除了異常噪聲的干擾，相比小波又進一步減小了慣導誤差。

5 結束語

小波消噪方法能在一定程度上抑制陀螺漂移，但是它有一系列固有缺陷：不能消除異常噪聲；分解精度受測不準原理影響；對小波基、分解尺度和閾值估計方法等依賴性太大，需要繁瑣的調試才能達到好的效果。EMD消噪方法根據陀螺信號本身特性進行分解得到IMF，無需任何事先的參數設置，消噪過程中所需的閾值及相關系數由原始陀螺信號及其IMF分量計算得到，因此EMD方法的消噪過程只與陀螺信號本身有關。利用EMD方法對陀螺信號消噪，既能夠抵制異常噪聲干擾，又能削弱高頻噪聲，減少陀螺各項誤差成分，提高慣導解算的精度。

［1］ HUO Ju，WANG Shijing，YANG Ming，et al.Noise Processing of FOG Signal Based on Wavelet Threshold－value［J］.Journal of Chinese Inertial Technology，2008，16（3）：343－347.（霍炬，王石靜，楊明，等.基于小波變換閾值法處理光纖陀螺信號噪聲［J］.中國慣性技術學報，2008，16（3）：343－347.）

［2］ WAN Yanhui，QING Yongyuan.Application of Wavelet Analysis in Gyro Signal Filtering［J］.Piezoelectrics and Acoustooptics，2005，27（4）：455－457.（萬彥輝，秦永元.小波分析在陀螺信號濾波中的研究［J］.壓電與聲光，2005，27（4）：455－457.）

［3］ WU Fumei，YANG Yuanxi.GPS／INS Integrated Navigation by Adaptive Filtering Based on Wavelet Threshold De－noising［J］.Acta Geodaetica et Cartographica Sinica，2007，36（2）：124－128.（吳富梅，楊元喜.基于小波閾值消噪自適應濾波的GPS／INS組合導航［J］.測繪學報，2007，36（2）：124－128.）

［4］ TANG Wei，LI Shixin，LIU Luyuan，et al.Select of Wavelet Basis in Gyro Signal Processing［J］.Journal of Chinese Inertial Technology，2002，10（5）：28－30.（湯巍，李士心，劉魯源，等.關于陀螺信號處理中小波基選取的研究［J］.中國慣性技術學報，2002，10（5）：28－30.）

［5］ LIU Bin，JIANG Jinshui，YU Weikai.EMD De－noising Method and Its Application in the Processing for Rolling Mill Signals［J］.Acta Metrologica Sinica，2009，30（1）：73－77.（劉彬，蔣金水，于偉凱.EMD相關度消噪及其在軋機信號處理中的應用［J］.計量學報，2009，30（1）：73－77.）

［6］ JIANG Li，LI Changyun.A Study of Wavelet Threshold Filtering Based on Empirical Mode Decomposition［J］.Signal Processing，2005，21（6）：659－662.（江力，李長云.基于經驗模分解的小波閾值濾波方法研究［J］.信號處理，2005，21（6）：659－662.）

［7］ HUANG N E，SHEN Z，LONG S R，et al.The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Nonstationary Time Series Analysis［C］∥Proceedings of Royal Society.London：［s.n.］，1998：903－993.

［8］ DAI Wujiao，DING Xiaoli，ZHU Jianjun，et al.EMD Filter Method and Its Application in GPS Multipath［J］.Acta Geodaetica et Cartographica Sinica，2006，35（11）：321－327.（戴吾蛟，丁曉利，朱建軍，等.基于經驗模式分解的濾波消噪法及其在GPS多路徑效應中的應用［J］.測繪學報，2006，35（11）：321－327.）

［9］ KATHLEEN T.Stratospheric and Tropospheric Signals Extracted Using the Empirical Mode Decomposition Method［D］.Washington：University of Washington，2003.

［10］ WANG Jian，GAO Jingxiang，WANG Jinling.GPS Baseline Solution Based on Empirical Mode Decomposition［J］.Acta Geodaetica et Cartographica Sinica，2008，37（1）：10－14.（王堅，高井祥，王金嶺.基于經驗模態(tài)分解的GPS基線解算模型［J］.測繪學報，2008，37（1）：10－14.）

［11］ CHEN Jun，XU Youli.Structural Damage Identification Based on Response of Forced Vibration［J］.Journal of Vibration，Measurement and Diagnosis，2005，25（2）：101－104.（陳雋，徐幼麟.經驗模分解在信號趨勢項提取中的應用［J］.振動、測試與診斷，2005，25（2）：101－104.）

［12］ FLANDRIN P，RILLING G，GONCALVES P.Empirical Mode Decomposition as a Filter Bank［J］.IEEE Signal Processing Letters，2004，11（2）：112－114.

［13］ HUANG N E，WU Zhaohua.A Review on Hilbert－Huang Transform：Method and Its Applications to Geophysical Studies［J］.Reviews of Geophysics，2008，46：1－23.

［14］ BOUDRAA A O，CEXUS J C，SAIDI Z.EMD－based Signal Noise Reduction［J］.International Journal of Signal Processing，2004，1：33－36.

［15］ DONOHO D L，JOHNSTONE I M.Adapting to Unknown Smoothness via Wavelet Shrinkage［J］.Journal of the American Statistical Association，1995，90（12）：1200－1224.

［16］ IEEE STD 952－1997.IEEE Standard Specification Format Guide and Test Procedure for Single－Axis Interferometric Fiber Optic Gyros［S］.［S.l.］：IEEE Standard Board，1997.

［17］ IEEE 1139.Definitions of Physical Quantities for Fundamental Frequency and Time Metrology－Random Instabilities［S］.［S.l.］：IEEE Standard Board，2008.

［18］ LI Xiaoying，HU Min，ZHANG Peng，et al.Applying Overlapping Allan Variance Theory to Better Stochastic Modeling of Microgyro［J］.Journal of Northwestern Polytechnical University，2007，25（2）：225－229.（李曉瑩，胡敏，張鵬，等.交疊式Allan方差在微機械陀螺隨機誤差辨識中的應用［J］.西北工業(yè)大學學報，2007，25（2）：225－229.）

De－noising Method for Gyro Signal Based on EMD

GAN Yu，SUI Lifen
Institute of Surveying and Mapping，Information Engineering University，Zhengzhou 450052，China

Gyro random drift is a remarkable factor that can affect the precision of inertial navigation system（INS）.Wavelet de－noising method is poor in coping with exceptional noise，and it depends greatly on the selection of wavelet base and decomposition scale.Empirical mode decomposition（EMD）de－noising method for gyro signal is presented.The signal is decomposed into an intrinsic mode function（IMF）group.Based on this group，IMFs of exceptional noise are first disposed by2σcriterion and then the number of IMFs of high frequency noise is determined by correlation coefficient.The de－noising process is finally done by removing the noisy IMFs.Detailed comparison between EMD method and wavelet method is given.Overlapping Allan variance is used to analyze the effect of the two methods，and the applicable ability of EMD method is tested through an INS calculation.It is shown that EMD method outperforms wavelet method in removing exceptional noise and is more efficient in weakening random drift.

gyro random drift；wavelet；empirical mode decomposition；de－noising

GAN Yu（1988—），male，postgraduate，majors in dynamic geodetic data processing.

1001－1595（2011）06－0745－06

P228

A

國家自然科學基金（40974010）；信息工程大學測繪學院碩士學位論文創(chuàng)新與創(chuàng)優(yōu)基金（S201101）

宋啟凡）

2010－12－06

2010－12－30

甘雨（1988—），男，碩士生，主要從事動態(tài)大地測量數據處理研究。

E－mail：ganyu099＠163.com