光纖陀螺基于本征模態函數篩選的閾值濾波算法

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

（軍械工程學院 軍械技術研究所，石家莊 050000）

光纖陀螺基于本征模態函數篩選的閾值濾波算法

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

（軍械工程學院 軍械技術研究所，石家莊 050000）

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

光纖陀螺；隨機噪聲；微分經驗模態分解；本征模態函數篩選

當前光纖陀螺（FOG）在航空航天等各個領域有著十分廣泛和重要的應用，然而受內部結構和外界環境等因素的影響，在實際使用過程中其輸出信號中往往伴有大量隨機噪聲。

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

考慮到光纖陀螺信號與某些噪聲頻率相近，EMD分解出的IMF頻率區段可能不會嚴格按照從高到低排列，產生波形混疊的現象。為了有效抑制光纖陀螺的隨機噪聲，本文在辨識光纖陀螺隨機噪聲特征基礎上，提出一種基于IMF篩選的微分經驗模態分解（DEMD）閾值濾波，在抑制噪聲的同時，提高了DEMD閾值濾波算法的運行效率。

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

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

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

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

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

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

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

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

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

2 光纖陀螺隨機噪聲濾波

2.1 DEMD濾波概述

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

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

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

2.2 IMF篩選準則

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

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

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

2.3 DEMD閾值濾波

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

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

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

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

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

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

如前文所述，前兩個IMF分量可認為只含噪聲，所以其標準差為

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

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

根據信號PSD與方差的關系，結合公式(12)可以得出：

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

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

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

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

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

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

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

⑤ 對位于3～m-1的 IMF進行時間序列區間的硬閾值處理得到

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

3 實驗與結果分析

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

3.1 數據仿真實驗

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

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

分析仿真結果可以發現：

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

圖4 標準Doppler信號與含噪Doppler信號Fig.4 Standard doppler signal and noised doppler signal

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

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

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

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

3.2 光纖陀螺濾波實驗

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

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

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

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

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

4 結 論

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

（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] 曾慶化, 黃磊, 劉建業, 等. 基于ARMA模型的光纖陀螺隨機噪聲濾波方法[J]. 中國慣性技術學報, 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] 黨淑雯, 田蔚風, 錢峰. 基于提升小波的光纖陀螺分形噪聲濾除方法[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參數估計方法的性能評價與分析[J]. 計算機應用與軟件, 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

國家自然科學基金（51308075）

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

聯 系 人：陳志斌（1965—），男，研究員，博士生導師。E-mail：shangxinboy@163.com

1005-6734(2015)05-0685-05

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