基于奇異譜分解的微機械加速度計振動噪聲抑制方法

伍宗偉， 姚敏立， 馬紅光， 馬幫立， 田方浩

(1. 第二炮兵工程大學(xué) 空間工程系，西安 710025;2. 鄖陽師專物理與電子工程系，十堰 442700)

低成本姿態(tài)估計系統(tǒng)廣泛應(yīng)用于導(dǎo)航控制、光電和慣性穩(wěn)定平臺、機器人，以及無人機等系統(tǒng)中[1-5]，而微機械(Micro-Electro Mechanical Systems, MEMS)慣性器件的應(yīng)用有效地降低了由于高精度慣性器件所帶來的昂貴成本。微機械慣性器件的成本較低，但是其性能較差，不僅測量噪聲與零偏變化較大，且振動也對其造成較大的影響。因此，在實際實用中需要采用數(shù)據(jù)融合方法來克服使用單一器件帶來的較大誤差。文獻(xiàn)[3-5]代表了工程上常用的低成本姿態(tài)估計方法，即利用陀螺的短時精度和加速度計的長時精度通過濾波算法融合得到較好的姿態(tài)估計，去除陀螺的零偏以及加速度計較大的高頻測量噪聲，算法簡單，但是需要克服兩個問題：①MEMS慣性器件受到振動噪聲影響很大，雖然卡爾曼濾波算法一定程度上可以認(rèn)為是低通濾波器，但是仍然無法去除疊加在有用姿態(tài)信號上的低頻振動噪聲；②由于加速度計依靠重力矢量測量得到姿態(tài)角，因此易受載體的機動影響。因此，在融合濾波器之前需要采用預(yù)處理方法來克服這些誤差和干擾。第二個問題已經(jīng)有較多應(yīng)對方法，如判斷開關(guān)以及自適應(yīng)濾波方法[4-5]，因此如何處理低頻振動噪聲是提高姿態(tài)估計精度的關(guān)鍵。

針對MEMS慣性器件的噪聲預(yù)處理，主要方法有小波分析以及經(jīng)驗?zāi)B(tài)分析(Empirical Mode Decomposition，EMD)兩種。文獻(xiàn)[6]代表了基于小波閾值方法去除高頻噪聲的方法，并且利用機動性判斷來確定小波分解層數(shù)，一定程度上克服了小波閾值去噪分解層數(shù)難以確定的問題。但是，該方法并沒有能夠有效地去除與有用姿態(tài)信號頻域疊加的低頻振動噪聲，且在一定意義上與后續(xù)的融合算法作用重復(fù)，達(dá)不到提高姿態(tài)融合精度的目的，反而增加了運算量。文獻(xiàn)[7]指出了小波分解的不足，即合適的小波基難以選取以及小波分解層數(shù)難以確定，提出了采用EMD的方法抑制陀螺的隨機漂移，EMD是完全數(shù)據(jù)驅(qū)動的一種自適應(yīng)數(shù)據(jù)分解方法，但是對于MEMS慣性器件來說，噪聲較大且噪聲的不確定性很強，分解出來的虛假模式很多，并且各個固有模態(tài)之間信息混疊，且姿態(tài)變化與振動噪聲并沒有顯著特性，很難提取出真實的振動模式和有用姿態(tài)信號，同時，EMD自身存在端點效應(yīng)等問題，限制了該方法在MEMS慣性傳感器噪聲處理中的應(yīng)用。另外，文獻(xiàn)[8]提出結(jié)合小波和EMD來進行去噪，不僅算法復(fù)雜，而且沒有從原理上克服二者的缺陷。

本文以車載平臺的水平姿態(tài)估計為研究對象，提出采用奇異譜分析(Singular Spectrum Analysis, SSA)的方法抑制對姿態(tài)融合影響較大的振動噪聲，該方法能夠有效地提取姿態(tài)變化的趨勢項，分離與振動噪聲有關(guān)的諧振項。跑車實驗的結(jié)果表明，SSA算法能夠有效地的分離有用信號和振動噪聲，提高低成本姿態(tài)估計的精度。

1 振動噪聲分析

1.1 姿態(tài)估計原理

分析噪聲對低成本姿態(tài)估計影響的前提是對姿態(tài)估計原理與融合濾波器的分析。首先簡要介紹姿態(tài)估計原理，陀螺積分可以得到姿態(tài)角，而加速度計的重力測量分量也能夠得到姿態(tài)角，定義坐標(biāo)系為：地理坐標(biāo)系取為東北天坐標(biāo)系，車體坐標(biāo)系的x,y,z軸分別指向車右側(cè)方向、車頭方向以及車體正上方。在無機動加速度干擾的情況下，加速度計測量姿態(tài)角的公式為[3]：

(1)

式中g(shù)為重力加速度值，fx、fy分別為兩個水平放置的加速度計測量值，n為測量噪聲。對于微機械加速度來說，受到振動噪聲影響較大。而陀螺的積分能夠提供較為精確的姿態(tài)角，但是由于零偏的存在使得姿態(tài)角存在長期誤差積累，因此將二者融合即可得到不受振動干擾和時間累積誤差的精確姿態(tài)角。通常選擇Kalman濾波器作為融合濾波器，將陀螺零偏擴增為狀態(tài)變量進行估計，選擇合適的噪聲陣，能夠有效抑制加速度計的測量噪聲n。由Kalman濾波器的性質(zhì)可知，選擇較大的測量噪聲陣R，可以抑制加速度計的測量噪聲n，相當(dāng)于一個低通濾波器，但是如果選擇的R陣過大，則測量值對系統(tǒng)狀態(tài)方程的修正量過小，會造成濾波器的擾動，圖1給出了一組在實驗中的估計結(jié)果，R陣中的測量噪聲方差分別選擇為5×10-3與5×10-1的比較結(jié)果(實際的測量噪聲方差為5×10-4，但是在Kalman濾波中一般適當(dāng)放大該值以取得較好的濾除高頻誤差效果)。

圖1 測量噪聲方程陣對Kalman濾波器的影響

從圖中可以看出，較大的噪聲陣能夠有助于抑制振動噪聲，但是會對Kalman濾波器本身帶來擾動，實際的應(yīng)用中選擇的R陣要兼顧精度和抗擾動，因此，低頻的振動噪聲一般難以通過Kalman濾波器得到較好的抑制。

1.2 振動噪聲

車輛的振動噪聲主要來自于三個方面的激勵：①路面不平度；②車輛發(fā)動機；③車輪/輪胎(內(nèi)部激勵，主要由于車輪不平衡以及車輪沖擊和輪胎質(zhì)量不均勻?qū)е?[9]。因此，為了研究載體振動噪聲，分別采樣載體不同運動狀態(tài)下的傳感器測量數(shù)據(jù)進行頻譜分析。MEMS陀螺與加速度計均來源于星網(wǎng)宇達(dá)公司的微機械慣性測量單元XW-ADU5220，分析結(jié)果如下圖2所示。

從頻譜分析中我們可看出：載體的振動噪聲對MEMS陀螺影響較小，變化較為明顯是第五種情況下，即土路的顛簸不平對姿態(tài)角的影響較大，從而對陀螺的測量幅值產(chǎn)生較大的影響；而加速度計則明顯受到載體振動噪聲影響，發(fā)動機的振動噪聲以及路面不平度引起的噪聲在全頻段均體現(xiàn)出來，靜止且發(fā)動機振動情況下，有明顯的振動峰值譜以及諧振譜，而車輛在低速平坦路面行駛時同樣有頻率大致在16 Hz和40 Hz附近的振動噪聲，且隨著車速以及路面顛簸度增強，該16 Hz附近的振動有加強趨勢，低頻段的振動噪聲也同樣加強。

圖2 載體不同運動狀態(tài)下的測量數(shù)據(jù)頻譜分析

圖3 低通之后的測量數(shù)據(jù)頻譜分析

首先想到的是通過低通濾波器濾除高頻噪聲。通過實際的行車實驗中參考系統(tǒng)提供的高精度姿態(tài)測量，可知車體的姿態(tài)變化一般低于5 Hz，因此首先設(shè)計低通濾波器濾除大部分的高頻噪聲。圖3顯示了通過一個5 Hz的低通濾波器之后的加速度計的測量值與實際通過高精度慣導(dǎo)計算出的加速度計值的比較，選擇的是慢速平坦的測量環(huán)境。

通過上圖我們看到，經(jīng)過低通濾波之后加速度計仍然含有較大的低頻噪聲，這個噪聲分量主要是由于車體振動引起的。通過大量的跑車實驗，發(fā)現(xiàn)該低頻振動噪聲的最大強度在1 Hz左右，實驗中采用的是全順面包車，測量設(shè)備均固連于車頂，對于不同的載體與安裝位置，均可進行簡單的頻譜分析得到該主要低頻振動噪聲的集中頻點。

2 奇異譜分析

奇異譜分析的目的是將原始的時間序列分解為一些序列之和，使得這些序列能夠分離成為趨勢項、周期項以及其他的非周期項，或者噪聲，然后選擇合適的序列重構(gòu)，從而達(dá)到處理信號的目的。實現(xiàn)的手段類似于經(jīng)驗?zāi)B(tài)分解或者小波分解，方法類似于主成分分析。由于該方法本質(zhì)上是一種利用協(xié)同延遲矩陣進行主成分分析的方法，因此可以看作是一種從能量域分析信號的方法。該方法不同于小波分析需要所分析信號的先驗知識，含有經(jīng)驗?zāi)B(tài)分解數(shù)據(jù)驅(qū)動的優(yōu)點，同時又很適合短數(shù)據(jù)的趨勢提取與諧波分量和噪聲分量的分離，已經(jīng)成功的應(yīng)用于氣象與地理學(xué)、電信以及金融等領(lǐng)域[10]。將加速度計的測量分為三個方面：①有用的姿態(tài)變化信號，即趨勢項；②振動引起的諧波分量；③測量噪聲。通過奇異譜分解可以分別對上述的三種分量進行基于數(shù)據(jù)的自適應(yīng)提取，其基本算法分為兩個互補的步驟：分解和重構(gòu)。具體的實現(xiàn)算法為[11]：

(1) 插值構(gòu)造

插值構(gòu)造實際上是一種將單維數(shù)據(jù)映射到多維的方法。假設(shè)我們有一個一維時間序列構(gòu)成的信號矢量X=[x1,…,xN]T，我們將其擴張為多維序列信號矩陣H=[X1,…,Xk]，其中的每一個Xi=[xi,…,xi+L-1]T，K=N-L-1，都是一個所謂延遲信號矢量，這里的時間窗口長度L是一個整數(shù)，其取值范圍為2

(2)

從上面的構(gòu)成過程看，實際上構(gòu)成的H陣為一個Hankel陣，且構(gòu)成過程就是一個擴維過程，擴張的維數(shù)由窗口長度L確定，SSA根據(jù)這個軌線矩陣H來確定主要分量，后續(xù)的分析實際上就是找到H的主分量以及每一個主分量對應(yīng)的子空間，然后選擇合適的信號子空間，對軌線矩陣進行投影。

(2) 奇異值分解

計算HHT的奇異值，實際上和PCA類似，求取其特征值λ1,λ2,…λL。由于HHT為對稱矩陣，則有HHT=VΛVT，這里的Λ為以λi為對角元素的對角陣，且λi滿足λ1≥λ2≥…≥λL>0，令每一個λi對應(yīng)的分量為Hi，我們得到

(3)

(3) 分類

將上述分解得到的一些初等矩陣H進行分組，令I(lǐng)1={i1,…,ip}是一個I1分組的下標(biāo)i1,…,ip集合，則組I就定義為HI+H…Hip，總共分為s組，則原來的H陣就變?yōu)椋篐I=HI1+…+HIs。對于給定的組，其占用的權(quán)重為:

一般情況下，都是直接分為兩類，一類為噪聲，一類為趨勢項，即s=2。

(4) 對角平均

所謂對角平均即是將上述給定的分組I中的每一個初等矩陣Hi對應(yīng)的一維信號恢復(fù)出來，也就是所謂的Hankel化。假定給定的任意矩陣為Z，其元素為zij，首先將其變換為Hankel矩陣形式Z，然后從Hankel矩陣求序列，也就是說第k個數(shù)值是將所有滿足i+j=k+2的HZ的元素求平均即可。這里的H算子其定義如下：對于任意的L×K矩陣Z=(zij)，L≤K，對于i+j=c，N=L+K-1，Z的元素可以通過下面的算法進行計算:

(4)

通過SSA算法可以有效地分離出趨勢項和諧振項，諧振項可以認(rèn)為是由于振動噪聲引起的，而趨勢項則反映了姿態(tài)的實際變化，由于Kalman濾波器融合陀螺積分和加速度計測量結(jié)果得到姿態(tài)估計，因此，這里的加速度計利用的是其長期精度，故取該趨勢項足以滿足Kalman濾波的要求，同時還能夠去除疊加在有用姿態(tài)信號上的低頻振動噪聲。

注意這里的趨勢項并不一定是SSA分析得到的分離信號中的第一個譜分量，由于SSA得到的信號序列是按照能量大小進行排列的，對于MEMS加速度計來說，很多情況下都是振動噪聲的能量要遠(yuǎn)大于有用的姿態(tài)變化能量，因此，這里需要采用一個合適的方法區(qū)分出振動項和趨勢項，一個簡單的準(zhǔn)則就是過零點檢測法，由于姿態(tài)變化是慢變的趨勢項，因此過零點是最少的。較為復(fù)雜和精確的方法是采用FFT進行頻域分離，在零頻處頻譜最大值即為趨勢項。

另外，在SSA實現(xiàn)過程中，僅有一個參數(shù)需要設(shè)定，即SSA算法中的時間窗口L，根據(jù)軌線矩陣H的構(gòu)成來看，每一列都是L個采樣點，因此，該參數(shù)應(yīng)該與需要SSA分析的信號中的主要周期項的周期大致可比擬才能得到較好的諧振項與趨勢項的分離。根據(jù)實驗中測量的主要振動項大致為1 Hz，而采樣頻率為100 Hz，故L設(shè)定為100。

3 跑車實驗數(shù)據(jù)分析

實驗設(shè)備包括星網(wǎng)宇達(dá)公司的微機械慣性測量單元XW-IMU5220，對比的參考系統(tǒng)為XW-ADU7612。其中XW-IMU5220包含有3個MEMS陀螺和3個MEMS加速度計，陀螺的零偏<0.08°/s，零偏穩(wěn)定性為<0.05°/s，而加速度計的零偏<0.005g，零偏穩(wěn)定性為<0.001g。XW-ADU7612能夠以100Hz的數(shù)據(jù)更新速率輸出被測對象精確的姿態(tài)角(航向、俯仰和橫滾角)，精度可達(dá)0.1°。如圖4所示，實驗設(shè)備安裝于全順車頂。

圖4 實驗載體及設(shè)備安裝照片

實驗的行車環(huán)境分別為圖1中的后兩種運動環(huán)境，即城市公路以及鄉(xiāng)村土路，行車速度大致為30km/h。

(1)城市公路

首先分析城市公路的數(shù)據(jù)處理情況。圖5給出了車輛近似直行情況下加速度計的測量值與高精度參考得到的數(shù)據(jù)對比結(jié)果。從圖5中可以看出，加速度計的測量值存在較大的低頻振動噪聲項，但是其趨勢項與參考值大致能夠保持一致，因此，如果能夠較好的分離出趨勢項與低頻振動噪聲項，則可以在最后的融合估計中對振動噪聲項進行抑制。下面我們通過SSA對橫滾角測量數(shù)據(jù)進行詳細(xì)分析，分離出前6個能量依次遞減的譜分量如圖6所示。

從圖6中可以明顯看出趨勢項與分離的振動項，對比圖6和圖5可以看出，SSA方法在提取趨勢項方面是很有效的，另外各個諧振項也能夠分離出來，由于低通之后的測量數(shù)據(jù)低頻振動較強，因此分離結(jié)果中的2、3，以及4、5分量大體一致。通過簡單的過零點檢測方法即可提取出趨勢項，趨勢項與參考值對比如圖7所示，同時給出主振動項，對比圖7中的提取結(jié)果與圖5原始的測量數(shù)據(jù)，更進一步說明了SSA方法在趨勢項提取以及振動項分離方面的優(yōu)勢。

圖5 加速度計的測量數(shù)據(jù)

圖6 橫滾角加速度計SSA分離的信號分量

圖7 橫滾角加速度計通過SSA提取的趨勢項與振動項

圖8 俯仰角加速度計通過SSA提取的趨勢項與振動項

同理，我們給出俯仰角加速度計的SSA處理結(jié)果，如圖8所示。

利用上述提取出的趨勢項進行Kalman濾波融合得到姿態(tài)估計，具體的計算過程見參考文獻(xiàn)[3,5]，測量方差選擇為5×10-2，得到的姿態(tài)估計結(jié)果對比如圖9所示，圖中給出了采用SSA抑制振動噪聲預(yù)處理之后的姿態(tài)估計誤差與直接進行濾波融合的估計誤差結(jié)果，表明了該方法能夠提高低成本姿態(tài)估計的精度。

圖9 姿態(tài)估計結(jié)果對比

(2) 鄉(xiāng)村土路

下面給出在鄉(xiāng)村土路的一組測量數(shù)據(jù)處理結(jié)果。由于鄉(xiāng)村土路上顛簸環(huán)境，因此加速度計測量結(jié)果中振動噪聲更大，因此，采用SSA預(yù)處理的效果將體現(xiàn)的更為明顯。圖10給出了橫滾加速度計原始測量值、SSA預(yù)處理之后的測量值以及參考值的對比。

圖10 鄉(xiāng)村土路SSA-ICA預(yù)處理結(jié)果

利用SSA方法預(yù)處理之后得到的姿態(tài)估計誤差與未處理的估計誤差對比如下圖11所示，這里由于是在顛簸土路，因此測量方差選擇為10×10-2。

圖11 鄉(xiāng)村土路姿態(tài)估計結(jié)果對比

從兩組不同路況下的實驗結(jié)果來看，提出的振動噪聲抑制方法均能夠有效地提高姿態(tài)估計精度。

4 結(jié) 論

為了提高低成本MEMS慣性傳感器估計姿態(tài)的精度，抑制載體振動噪聲對MEMS加速度計帶來的影響，本文提出了采用SSA方法來去除慣性傳感器中的低頻振動噪聲。利用SSA對趨勢項以及周期項能夠有效分離的優(yōu)勢，抑制低頻振動噪聲，實驗結(jié)果表明，提出的傳感器噪聲預(yù)處理方法能夠有效地抑制低頻振動噪聲，提高通過姿態(tài)融合算法得到的姿態(tài)角的估計精度。同時，該方法有利于提高微機械慣性器件的工程實用性。

參 考 文 獻(xiàn)

[1]Maeder U, Morari M. Attitude estimation for vehicles with partial inertial measurement [J]. IEEE Transactions on Vehicular Technology, 2011, 60(4): 1496-1504.

[2]Lai Y C, Jan S S. Attitude estimation based on fusion of gyroscopes and single antenna GPS for small UAVs under the influence of vibration [J]. GPS Solutions, 2011, 15(1): 67-77.

[3]Wu Z W, Yao M L, Ma H G, et al. Low-cost attitude estimation with MIMU and two-antenna GPS for Satcom-on-the-move[J]. GPS Solutions, 2013, 17(1): 75-87.

[4]Suh Y S. Orientation estimation using a quaternion-based indirect Kalman filter with adaptive estimation of external acceleration [J]. IEEE Transactions on Instrument and Measurement, 2010, 59(12): 3296-3305.

[5]伍宗偉, 姚敏立, 馬紅光, 等. 移動衛(wèi)星通信低成本多傳感器融合姿態(tài)估計方法 [J]. 西安交通大學(xué)學(xué)報, 2012, 46(12): 55-61.

WU Zong-wei, YAO Min-li, MA Hong-guang, et al. Low-cost attitude estimation method using multi-sensor fusion for mobile satellite communication [J]. Journal of Xi’an Jiaotong University, 2012, 46(12): 55-61.

[6]沈曉衛(wèi), 姚敏立, 常瑞花. 模糊邏輯提升小波在慣性傳感器去噪中的應(yīng)用 [J]. 振動與沖擊, 2011, 30(12): 238-242.

SHEN Xiao-wei, YAO Min-li, CHANG Rui-hua. Application of fuzzy logic and lifted wavelet transformation in de-noising of inertial sensors [J]. Journal of vibration and shock, 2011, 30(12): 238-242.

[7]甘 雨，隋立芬. 基于EMD-小波隨機消噪模型的GPS /INS 組合導(dǎo)航[J]. 測繪學(xué)報, 2011, 40(6): 745-750.

GAN Yu, SUI Li-fen. De-noising method for gyro signal based on EMD [J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 745-750.

[8]王 堅, 李增科, 高井祥, 等.基于EMD-小波隨機消噪模型的GPS/INS組合導(dǎo)航 [J]. 東南大學(xué)學(xué)報: 自然科學(xué)版, 2012, 42(3): 406-412.

WANG Jian, LI Zeng-ke, GAO Jing-xiang, et al. EMD-wavelet based stochastic error reducing model for GPS /INS integrated navigation [J]. Journal of Southeast University: Natural Science Edition, 2012, 42(3): 406-412.

[9]Mitschke M, Wallentowitz H. 汽車動力學(xué) [M]. 第四版. 陳蔭三，余 強，譯. 北京：清華大學(xué)出版社，2009.

[10]Harris T J, Yuan H. Filtering and frequency interpretations of singular spectrum analysis [J]. Physica D, 2010, 239: 1958-1967.

[11]Ma H G, Jiang Q B, Liu Z Q, et al. A novel blind source separation method for single-channel signal [J]. Signal Processing, 2010, 90: 3232-3241.