光纖陀螺基于本征模態(tài)函數(shù)篩選的閾值濾波算法

肖文健，陳志斌，馬東璽，張 勇，劉先紅

（軍械工程學(xué)院 軍械技術(shù)研究所，石家莊 050000）

光纖陀螺基于本征模態(tài)函數(shù)篩選的閾值濾波算法

肖文健，陳志斌，馬東璽，張 勇，劉先紅

（軍械工程學(xué)院 軍械技術(shù)研究所，石家莊 050000）

為降低光纖陀螺隨機噪聲，提高其測量精度，利用周期圖法辨識光纖陀螺的隨機噪聲特征參數(shù)，針對其噪聲特征，提出了基于本征模態(tài)函數(shù)篩選的微分經(jīng)驗?zāi)B(tài)分解閾值濾波算法。以本征模態(tài)函數(shù)和原始信號二者的概率密度函數(shù)的空間距離為判別依據(jù)，對所有本征模態(tài)函數(shù)進(jìn)行篩選，根據(jù)已估計的噪聲參數(shù)計算閾值大小，采用時間序列閾值的方法對篩選出的本征模態(tài)函數(shù)進(jìn)行處理。仿真和實驗結(jié)果表明，該濾波算法能夠在跟蹤光纖陀螺信號變化的同時，使其零偏不穩(wěn)定性下降90.35%，角隨機游走下降93.75%，對隨機噪聲有較好的抑制能力。

光纖陀螺；隨機噪聲；微分經(jīng)驗?zāi)B(tài)分解；本征模態(tài)函數(shù)篩選

當(dāng)前光纖陀螺（FOG）在航空航天等各個領(lǐng)域有著十分廣泛和重要的應(yīng)用，然而受內(nèi)部結(jié)構(gòu)和外界環(huán)境等因素的影響，在實際使用過程中其輸出信號中往往伴有大量隨機噪聲。

光纖陀螺隨機噪聲的濾波方法主要分為建模濾波和非建模濾波。經(jīng)典的建模濾波方法就是根據(jù)隨機噪聲時間序列建立ARMA模型，然后利用卡爾曼濾波實時更新該模型的輸出[1-2]。這種建模濾波方法在隨機噪聲模型建立準(zhǔn)確的情況下，具有較好的濾波效果，然而在實際應(yīng)用中很難得到準(zhǔn)確的隨機噪聲模型。非建模濾波則不需要建立誤差模型。傳統(tǒng)非建模濾波主要利用FIR、IIR等數(shù)字低通濾波器來濾除高頻噪聲，當(dāng)信號和噪聲的頻帶相近時效果較差。后來小波變換被應(yīng)用到了光纖陀螺隨機噪聲的濾波中，小波變換具有良好的時頻局部化特性[3]。然而小波分解與傅里葉分解類似，分解所用的基函數(shù)固定，不能自適應(yīng)的跟蹤信號變化。1998年Huang提出經(jīng)驗?zāi)B(tài)分解（EMD），該方法與傳統(tǒng)信號分解方法的不同之處在于它不用選擇基函數(shù)，而是據(jù)信號本身特征自適應(yīng)地將信號分解成若干本征模態(tài)函數(shù)(IMF)[4]。傳統(tǒng)EMD去噪方法需要分辨哪些IMF是信號哪些IMF是噪聲，再進(jìn)行信號重構(gòu)。而實際中多數(shù) IMF既含信號又含噪聲，Kopsinis等人受小波閾值濾波的啟發(fā)提出了EMD閾值濾波[5]，然而該方法只考慮了白噪聲，沒有分析有色噪聲的濾波效果。Yu Gan等人指出慣性傳感器的隨機噪聲是分形噪聲，并通過改進(jìn)閾值計算方法，使EMD閾值濾波在分形噪聲濾波方面的得到應(yīng)用[6]。在利用EMD濾波時，多數(shù)文獻(xiàn)都是直接把EMD分解的第一個IMF分量當(dāng)作噪聲，然后對后面每個IMF分量分別進(jìn)行閾值處理，而對IMF分量的篩選問題很少研究，這樣不僅計算量大而且很難排除噪聲干擾。

考慮到光纖陀螺信號與某些噪聲頻率相近，EMD分解出的IMF頻率區(qū)段可能不會嚴(yán)格按照從高到低排列，產(chǎn)生波形混疊的現(xiàn)象。為了有效抑制光纖陀螺的隨機噪聲，本文在辨識光纖陀螺隨機噪聲特征基礎(chǔ)上，提出一種基于IMF篩選的微分經(jīng)驗?zāi)B(tài)分解（DEMD）閾值濾波，在抑制噪聲的同時，提高了DEMD閾值濾波算法的運行效率。

1 光纖陀螺隨機噪聲特征辨識

在 IEEE關(guān)于光纖陀螺測試標(biāo)準(zhǔn)中指出，光纖陀螺的隨機噪聲包括量化噪聲、角隨機游走、速率隨機游走、零偏不穩(wěn)定性和速率斜坡噪聲[7]。這 5類噪聲指標(biāo)通常用來評價陀螺性能，而在實際對某一光纖陀螺的隨機噪聲進(jìn)行濾波時，只需已知陀螺隨機噪聲表現(xiàn)出的總體噪聲特征即可。根據(jù)當(dāng)前國內(nèi)外相關(guān)的研究，光纖陀螺隨機噪聲表現(xiàn)為分形噪聲的特征。

分形噪聲是白噪聲的一般化，它的統(tǒng)計特征取決于只有一個參數(shù) H（0＜H＜1）的二階矩。分形噪聲被定義為一個零均值穩(wěn)定高斯過程，其自相關(guān)序列為[9]

式中，τ為時延，σ2為RH(0)的方差。

辨識光纖陀螺隨機噪聲特性的主要任務(wù)就是計算參數(shù)H。目前參數(shù)H的估計方法有聚合序列方差法、R/S法、周期圖法、Whittle法等[8]，根據(jù)光纖陀螺信號特征，本文選用周期圖法。對于一個離散時間序列，將其自相關(guān)函數(shù)進(jìn)行離散傅里葉變換，可以得到其功率譜密度函數(shù)為

式中：f表示頻率，C表示常數(shù)。其功率譜密度函數(shù)可以近似為

對公式(3)兩邊取對數(shù)可以得到：

本文選用的實驗對象為俄羅斯Fizoptika公司生產(chǎn)的VG910型開環(huán)光纖陀螺，將其靜置于隔振平臺上采集其隨機噪聲并繪制周期圖log- log曲線以及擬合直線如圖1所示，其分形噪聲參數(shù)H=0.879。

圖1 光纖陀螺周期圖log–log曲線及擬合直線Fig.1 Log–log periodogram and fitting line of FOG

2 光纖陀螺隨機噪聲濾波

2.1 DEMD濾波概述

對于傳統(tǒng)EMD方法，若分解的信號頻率比小于1.5或振幅比例小時，EMD 將不能區(qū)分兩種成分，從而造成混波現(xiàn)象。DEMD能提高混合信號的振幅比，因而在一定程度上改善EMD中的混波現(xiàn)象[9]。信號通過DEMD分解流程如圖2所示。

圖2 DEMD分解流程圖Fig.2 Flowchart of differential empirical mode decomposition

DEMD算法將原信號經(jīng)過對時間一次求導(dǎo)，使新信號的能量盡可能地按頻率從高到低遞減，提高了對信號頻率的分辨能力，這樣就為DEMD在濾波方面的應(yīng)用提供了可能。重構(gòu)信號時，選取不同位置的IMF分量會產(chǎn)生不同的濾波效果，因此在DEMD濾波中，IMF篩選是非常關(guān)鍵的一個步驟。

2.2 IMF篩選準(zhǔn)則

因為DEMD分解出的IMF分量從前到后頻率逐漸減小，所以IMF分量按前后順序可以分為只含噪聲IMF、信號與噪聲混合IMF以及只含信號IMF。目前在大多數(shù)文獻(xiàn)中均把第一個或者前兩個IMF分量直接當(dāng)成純噪聲，然后再對后續(xù)的IMF分量進(jìn)行相應(yīng)處理。很顯然，這種處理方法比較粗糙，不對所有IMF分量進(jìn)行有效篩選，全部交由后續(xù)處理無疑會增加后續(xù)處理的工作量。所以采用DEMD濾波首先要對IMF分量進(jìn)行篩選，找出信號與噪聲混合 IMF同只含信號IMF分量的分界點m。這樣對于前兩個IMF分量當(dāng)作噪聲直接舍去，對于分界點m以后的IMF分量不做處理，只需要對中間處于3～m-1(m≥4)的IMF進(jìn)行閾值處理，最后可按照公式(5)對信號進(jìn)行重構(gòu)。

Ali Komaty等人從概率密度函數(shù)（PDF）的角度詳細(xì)分析了每個IMF分量與信號之間的關(guān)系，并指出每個IMF分量PDF與原始信號PDF之間的相似度能夠反映該IMF分量所含信號的比重[10]。受此啟發(fā)，本文將IMF分量PDF與原始信號PDF之間的相似度用來作為 IMF篩選的標(biāo)準(zhǔn)。假設(shè)原始信號的 PDF為P( k)，某一IMF分量的PDF為Q( k)，那么它們的相似度定義為

式中，N為濾波窗口內(nèi)數(shù)據(jù)長度。分別計算每個IMF分量與原始信號的相似度即可得到每個 IMF分量與原始信號x( t)的相似度序列D( i), i=1～L。 D( i)序列中第一個局部極大值后開始減小的點即為要求取的分界點m，即：

2.3 DEMD閾值濾波

在DEMD分解得到的IMF分量中，有一部分既包含噪聲又包含信號。對于這一部分IMF分量，如果當(dāng)作噪聲直接舍去會丟失部分有用信息，如果當(dāng)作信號直接采用則會引入部分噪聲。對于這部分IMF分量可以對其進(jìn)行閾值處理。

參照小波閾值濾波的思想，常用的閾值函數(shù)主要包括硬閾值函數(shù)和軟閾值函數(shù)。本文通過實驗發(fā)現(xiàn)，在對光纖陀螺隨機噪聲濾波過程中，使用軟閾值濾波會丟失信號許多細(xì)節(jié)信息，與原始信號偏差較大，因此采用如下硬閾值函數(shù)：

式中，Ti表示第i個IMF分量的閾值。

另外，對于每個IMF分量都有許多過零點，而且排位越靠前的過零點越多。考慮到每個過零點附近的IMF值比較小，低于設(shè)定的閾值，但這些數(shù)據(jù)也會包含有用信息，如果直接把低于閾值的數(shù)據(jù)置為零則會造成濾波后的信息部分缺失。因此本文采用區(qū)間閾值處理的方法來保護(hù)部分低于閾值的有用信息。對某一IMF分量，以相鄰零點之間的時間序列為一個區(qū)間，若此區(qū)間內(nèi)的極值大于閾值則保留，小于閾值則置零。如圖3所示，左側(cè)為原始信號，中間為采用直接硬閾值處理的結(jié)果，右側(cè)為區(qū)間硬閾值處理的結(jié)果。從圖 3中可以清楚地看出，采用時間序列區(qū)間硬閾值處理既能夠較好地去除噪聲干擾，又能較完整地保留有用信息。

圖3 直接閾值與時間序列區(qū)間閾值Fig.3 Direct thresholding and interval thresholding

對于每個IMF分量閾值的選擇，比較經(jīng)典的計算方法為

如前文所述，前兩個IMF分量可認(rèn)為只含噪聲，所以其標(biāo)準(zhǔn)差為

由于后續(xù)IMF分量既包含噪聲又包含信號，所以后續(xù)IMF分量的閾值不能直接用上式計算。參考文獻(xiàn)[6]中給出了分形噪聲經(jīng)過EMD分解后每個IMF分量PSD之間的關(guān)系為

式中，k′＞k≥2，ρH近似滿足

根據(jù)信號PSD與方差的關(guān)系，結(jié)合公式(12)可以得出：

由式(13)即可計算后續(xù)IMF中包含噪聲的方差為

這樣，利用公式(10)和公式(14)就可以得到整個IMF分量的閾值為

綜上所述，本文提出的濾波算法具體過程為：

① 分析光纖陀螺隨機噪聲類型，利用周期圖法確定分形噪聲參數(shù)H；

② 對光纖陀螺原始信號進(jìn)行DEMD分解，得到IMF分量；

③ 對所有IMF分量進(jìn)行篩選，確定信號與噪聲混合IMF同只含只含信號IMF的分界點m；

④ 計算每個IMF分量的噪聲閾值；

⑤ 對位于3～m-1的 IMF進(jìn)行時間序列區(qū)間的硬閾值處理得到

⑥ 利用公式(8)重構(gòu)光纖陀螺輸出信號。

3 實驗與結(jié)果分析

為了驗證本文所提出算法，分別進(jìn)行了數(shù)據(jù)仿真實驗和實際光纖陀螺濾波實驗。

3.1 數(shù)據(jù)仿真實驗

在Matlab軟件中的標(biāo)準(zhǔn)測試信號中，Doppler信號是一種連續(xù)變化信號，頻率覆蓋范圍大，相對其它測試信號更接近光纖陀螺實際信號，因此本文選用Doppler信號進(jìn)行仿真。在Doppler信號基礎(chǔ)上，分別加入3種不同參數(shù)的分形噪聲以驗證本文所提算法對不同參數(shù)噪聲的適應(yīng)能力。仿真所用分形噪聲的具體參數(shù)為 H=0.2, 0.5, 0.8，所合成含噪信號的信噪比（SNR）均設(shè)定為5 dB。

以H=0.8的仿真過程為例對整個仿真過程進(jìn)行說明。圖4中粗實線表示Doppler標(biāo)準(zhǔn)信號，細(xì)虛線表示疊加了分形噪聲的含噪信號。含噪信號經(jīng)DEMD分解為9個IMF分量；然后利用本文提出的篩選準(zhǔn)則對IMF分量進(jìn)行篩選，確定臨界點m=5；最后對含噪聲IMF分量進(jìn)行閾值濾波。為了說明該算法的濾波效果，這里將本文所述算法與小波閾值去噪以及傳統(tǒng)基于EMD部分重構(gòu)濾波進(jìn)行了對比，如圖 5所示。對于H=0.5和H=0.2時，使用同樣步驟對信號進(jìn)行濾波。分別計算在不同噪聲參數(shù)下三種方法濾波后信號的均方根（RMS）誤差，見表1。

分析仿真結(jié)果可以發(fā)現(xiàn)：

① 當(dāng) H=0.5時，小波閾值去噪方法與本文提出的方法濾波效果接近；當(dāng) H≠0.5時，小波閾值去噪方法的濾波效果降低，而本文提出的方法濾波效果基本不變。這是因為傳統(tǒng)小波閾值對白噪聲去噪效果較好，當(dāng)噪聲參數(shù)H=0.5時，信號的噪聲為白噪聲，所以此時濾波效果較好；而當(dāng)H≠0.5時，信號的噪聲為分形噪聲，此時小波閾值去噪的效果則會下降。

圖4 標(biāo)準(zhǔn)Doppler信號與含噪Doppler信號Fig.4 Standard doppler signal and noised doppler signal

圖5 三種方法濾波結(jié)果對比Fig.5 Filtering results of three denoising methods

表1 三種方法濾波后信號的RMS誤差Tab.1 RMS Errors of the signal after filtering

② 傳統(tǒng)基于EMD的部分重構(gòu)濾波效果受噪聲參數(shù)變化的影響不大，但由于該方法在對信號進(jìn)行重構(gòu)時IMF仍帶有部分隨機噪聲，所以該方法的濾波效果總體較差。

總體而言，仿真結(jié)果可以說明本文提出經(jīng)過IMF篩選的DEMD閾值濾波效果較好，而且受噪聲參數(shù)變化的影響很小。

3.2 光纖陀螺濾波實驗

為了驗證該濾波算法對光纖陀螺動態(tài)信號的跟蹤與降噪能力，本文通過轉(zhuǎn)臺的搖擺運動向光纖陀螺輸入動態(tài)信號。根據(jù)上文結(jié)果，實驗所用光纖陀螺的分形噪聲參數(shù)。對于光纖陀螺實際輸出選用滑動窗口的方式濾波，濾波窗口N=500，濾波窗口的長度需綜合考慮對濾波性能的要求，從而預(yù)先確定。設(shè)定搖擺運動的角位置滿足，那么對其求導(dǎo)即可得到光纖陀螺的輸入角速度實驗中以100 Hz的采樣頻率采集10個搖擺周期的光纖陀螺輸出數(shù)據(jù)。光纖陀螺的輸出信號經(jīng)DEMD分解為12個IMF分量；然后利用本文所述IMF篩選準(zhǔn)則對其進(jìn)行篩選，確定分界點m=8；再對篩選后的IMF進(jìn)行閾值濾波處理，最終得到的濾波結(jié)果。光纖陀螺的信號原始信號與濾波后的信號如圖6所示。

圖6 濾波前后信號光纖陀螺信號Fig.6 FOG signal before and after filtering

為了定量的說明本文所述方法的濾波效果，分別計算濾波前后光纖陀螺信號 RMS誤差以及用 Allan方差分析法獲得濾波前后光纖陀螺的角隨機游走和零偏不穩(wěn)定性，如表2所示。

表2 光纖陀螺濾波前后噪聲參數(shù)Tab.2 Noise parameters of FOG before and after filtering

從對光纖陀螺實際濾波的結(jié)果來看，本文所述算法在有效跟蹤其信號變化的前提下，對其隨機噪聲具有較好的濾波效果。經(jīng)過濾波后，該光纖陀螺的零偏不穩(wěn)定性和角隨機游走都有顯著的下降，其中 RMS誤差下降了 93.84%，零偏不穩(wěn)定性下降了 90.35%，角隨機游走下降了93.75%。

4 結(jié) 論

本文首先對光纖陀螺的噪聲特征參數(shù)進(jìn)行了辨識。針對其噪聲特征提出了基于IMF篩選的DEMD時間序列閾值濾波算法。該算法首先以IMF和原始信號PDF的相似度為依據(jù)對IMF進(jìn)行篩選，減少了IMF的處理個數(shù)；然后采用時間序列閾值處理的方法，對信號的噪聲進(jìn)行了有效濾除。采用滑動窗口的方法使得該算法能夠在線運行。實驗結(jié)果表明，本文所提出的算法相比傳統(tǒng)的小波閾值去噪以及EMD部分重構(gòu)去噪具有一定優(yōu)勢，在跟蹤光纖陀螺信號變化的同時對于光纖陀螺的隨機噪聲具有良好的抑制能力。

（References）：

[1] Li Xue, Wang Qin. A novel Kalman filter for combining outputs of MEMS gyroscope array[J]. Measurement, 2012, 45(4): 745-754.

[2] 曾慶化, 黃磊, 劉建業(yè), 等. 基于ARMA模型的光纖陀螺隨機噪聲濾波方法[J]. 中國慣性技術(shù)學(xué)報, 2015, 23 (1): 120-124. Zeng Qing-hua, Huang Lei, Liu Jian-ye, et al. Real-time filtering methods of FOG random noise based on ARMA model[J]. Journal of Chinese Inertial Technology, 2015, 23(1): 120-124.

[3] 黨淑雯, 田蔚風(fēng), 錢峰. 基于提升小波的光纖陀螺分形噪聲濾除方法[J]. 中國激光, 2009, 36(3): 625-629. Dang Shu-wen, Tian Wei-feng, Qian Feng. De-noising fractional noise in fiber optic gyroscopes based on lifting wavelet[J]. Chinese Journal of Lasers, 2009, 36(3): 625-629.

[4] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[C]//Proc. R. Soc. Lond. A, 1998, 454(3): 56-78.

[5] Kopsinis Y, McLaughlin S. Development of EMD-based denoising methods inspired by wavelet thresholding[J]. IEEE Transactions on Signal Process, 2009, 57(4): 1351-1362.

[6] Gan Yu, Sui Li-fen, Wu Jiang-fei, et al. An EMD threshold de-noising method for inertial sensors[J]. Measurement, 2014, 49: 34-41.

[7] IEEE Std 952-1997. IEEE standard speciation format guide and test procedure for single-axis interferometric fiber optic gyros[S]. New York: IEEE Standards Board, 1997.

[8] 韓忠明, 趙慶展, 李偉. Hurst參數(shù)估計方法的性能評價與分析[J]. 計算機應(yīng)用與軟件, 2010, 27(9): 56-58. Han Zhong-ming, Zhao Qing-zhan, Li Wei. Performance evaluation and analysis of Hurst parameter estimate methods[J]. Computer Applications and Software, 2010, 27(9): 56-58.

[9] Li Ming, Li Fu-cai, Jing Bei-bei, et al. Multi-fault diagnosis of rotor system based on differential-based empirical mode decomposition[J]. Journal of Vibration and Control, 2013, 20(3): 1-17.

[10] Komaty A, Boudraa A O, Augier B, et al. EMD-based filtering using similarity measure between probability density functions of IMFs[J]. IEEE Transactions on Instrumentation and Measurement, 2014, 63(1): 27-34.

De-noising method of FOG with DEMD threshold based on IMF criterion

XIAO Wen-jian, CHEN Zhi-bin, MA Dong-xi, ZHANG Yong, LIU Xian-hong
(Ordnance Technology Institute, Ordnance Engineering College, Shijiazhuang 050000, China)

To reduce the random noise of fiber optic gyro (FOG) and improve its measurement precision, the characteristic parameters of FOG random noise are identified by the method of periodogram. According to the noise characteristics, a de-noising method with differential empirical mode decomposition (DEMD) interval thresholding is proposed based on intrinsic mode function (IMF) criterion. The relevant IMFs are selected based on the space distance of the probability density functions (PDF) between the original signal and each IMF. The filtering thresholdings are calculated by the estimated noise parameters, and the relevant IMFs are filtered by the DEMD interval thresholding. Simulation and test results show that the FOG’s bias instability is decreased by 90.35% and the FOG’s angle random walk is decreased by 93.75% after filtering by the proposed method.

FOG; random noise; differential empirical mode decomposition; intrinsic mode functions criterion

V249.3

：A

2015-06-05；

：2015-09-21

國家自然科學(xué)基金（51308075）

肖文健（1989—），男，博士研究生，從事光電測量與信號處理研究。E-mail：xiao_wen_jian@163.com

聯(lián) 系 人：陳志斌（1965—），男，研究員，博士生導(dǎo)師。E-mail：shangxinboy@163.com

1005-6734(2015)05-0685-05

10.13695/j.cnki.12-1222/o3.2015.05.022