靳凱迪，柴洪洲，宿楚涵，向民志

信息工程大學地理空間信息學院，河南 鄭州 450001

由于元件構造和環(huán)境影響，光纖陀螺(fiber optic gyroscope，F(xiàn)OG)會產(chǎn)生復雜的噪聲源，直接影響慣性導航系統(tǒng)初始對準和導航解算整個過程的精度。光纖陀螺噪聲主要由系統(tǒng)性噪聲與隨機噪聲組成，系統(tǒng)性噪聲可通過試驗標校有效補償，隨機噪聲由于性質(zhì)復雜、隨機性強、補償難度大，成為衡量光纖陀螺精度的重要指標[1-3]。因此，對光纖陀螺隨機噪聲進行數(shù)學建模并通過濾波等手段加以補償，對準確分析光纖慣性導航系統(tǒng)誤差和提高導航解算精度等方面具有重要意義。目前，用于光纖陀螺隨機誤差建模的方法主要包括Allan方差分析、功率譜密度分析和時間序列分析等[4-8]。文獻[9]使用自回歸滑動平均模型(auto-regressive moving average，ARMA)建模光纖陀螺隨機噪聲，并給出原始觀測序列的Kalman濾波方程，但并未考慮MA模型存在時狀態(tài)噪聲為有色噪聲的情況。文獻[10]使用高階AR模型計算ARMA模型的噪聲估值，達到白化有色噪聲的目的，但其作為一種近似方法并不嚴密。后續(xù)研究中或是將隨機噪聲建模為AR模型[11]，或是沿用文獻[9—10]的思想[12-13]。針對隨機噪聲濾波模型先驗信息不準確的問題，文獻[11]使用Sage-Husa自適應濾波(Sage-Husa adaptive Kalman filter，SHAKF)在線估計系統(tǒng)噪聲和量測噪聲。文獻[12]使用動態(tài)Allan方差分析法(dynamic Allan variance，DAVAR)估計量測噪聲。文獻[7]針對自適應濾波參數(shù)耦合問題，使用DAVAR和SHAKF分別估計狀態(tài)噪聲和量測噪聲。然而，SHAKF易喪失矩陣正定性或非正定性[14-16]，尤其在ARIMA模型收斂初期，此現(xiàn)象更加顯著。

為此，本文使用并擴展了Harvey法,將光纖陀螺差分自回歸滑動平均模型(auto-regressive integrated moving average，ARIMA)構造為原始觀測序列的Kalman濾波方程，實現(xiàn)MA參數(shù)存在時的有色噪聲白化。針對濾波模型及先驗噪聲不準確的問題，在使用序貫平差更新ARIMA模型參數(shù)的同時，計算狀態(tài)噪聲遞推解。針對模型更新參數(shù)與狀態(tài)噪聲相互耦合，分析了動態(tài)Allan方差估計量測噪聲的不準確性，使用變分貝葉斯自適應濾波(variational Bayesian adaptive Kalman filter,VBAKF)修正狀態(tài)噪聲和量測噪聲，實現(xiàn)Kalman濾波函數(shù)模型與隨機模型的實時更新。實測數(shù)據(jù)表明，該方法可有效降低光纖陀螺隨機誤差，提高陀螺輸出精度。

1 光纖陀螺隨機噪聲建模

ARIMA模型實質(zhì)上是差分運算與ARMA模型的融合。ARIMA建模光纖陀螺隨機噪聲序列的主要內(nèi)容包括噪聲隨機序列平穩(wěn)性檢驗、ARIMA模型階數(shù)確定、模型參數(shù)估計及適用性檢驗等[17-20]。ARMA模型要求平穩(wěn)的時間序列，因此，首先要對采樣光纖陀螺原始隨機噪聲進行平穩(wěn)性和隨機性檢驗，如檢驗不通過，則需進行差分直至滿足平穩(wěn)隨機過程。使用中海達iPos光纖慣導以100 Hz采樣頻率采集靜態(tài)數(shù)據(jù)如圖1所示。采用單位根檢驗(包括ADF檢驗與KPSS檢驗)對信號檢測通過，因此，ARIMA模型差分階數(shù)為0。設ARIMA(p,0,q)模型如下[21]

圖1 FOG隨機噪聲原始輸出信號Fig.1 Original random noise signals of FOG

(1)

式中，xk為歷元k時隨機噪聲值；ai為自回歸項系數(shù)；bj為滑動平均項系數(shù)；εk為未知方差的白噪聲；c為常值。

使用前2000個歷元數(shù)據(jù)建立ARIMA模型，計算過程分為模型階數(shù)(p,q)確定和模型系數(shù)(ai，bj)計算。(p,q)可以通過時間序列的偏自相關函數(shù)圖確定，但由偏自相關函數(shù)圖得到的ARIMA模型階數(shù)往往不準確。因此采用赤池信息準則(akaike information criterion，AIC)與貝葉斯信息準則(Bayesian information criterionm，BIC)確定隨機噪聲序列ARIMA階數(shù)[6]，計算結果為p=3，q=3。模型階數(shù)確定后，通過最小二乘法計算模型系數(shù)值，所得ARIMA模型見表1。

表1 光纖陀螺隨機噪聲ARIMA(p,0,q)模型參數(shù)Tab.1 Parameters for ARIMA(p,0,q) of FOG random nosie

圖2為對所建模型采用標準化殘差檢驗結果。由圖2可知，所得殘差序列符合正態(tài)分布。Durbin-Watson統(tǒng)計量是常用的自相關統(tǒng)計量，設殘差序列為vk(k=1,2,…,n)，構造DW統(tǒng)計量為[22]

圖2 標準化殘差檢驗Fig.2 Standardized residuals test

(2)

式中，ρ為自相關因子。ρ越接近0，等價于DW統(tǒng)計量越接近2，序列越不存在一階相關性，本次殘差序列相關性檢驗所得結果為1.953 9，認為殘差序列符合高斯正態(tài)分布。綜上，認為所建模型較為準確，通過適用性檢驗。

2 ARIMA模型有色噪聲Kalman濾波

在光纖陀螺隨機噪聲建模為ARIMA模型后，使用Kalman濾波是進行隨機噪聲降噪的有效方法。以往處理存在MA系數(shù)的ARIMA模型時，往往忽視噪聲的相關性或使用噪聲估值進行白化。為此，建立基于Harvey法的Kalman濾波模型白化ARIMA有色噪聲。

2.1 ARIMA模型濾波方程

獲得如式(1)所示的隨機噪聲ARIMA(p,0,q)模型后，令狀態(tài)向量Xk=[xkxk-1…xk-p+1c]T，建立Kalman濾波狀態(tài)方程和量測方程為[9]

(3)

2.2 Harvey法有色噪聲濾波方程

當存在MA系數(shù)時，式(3)中的狀態(tài)噪聲為有色噪聲，此時，Kalman濾波將會導致結果發(fā)散或失真[14，21]。傳統(tǒng)上，通常將有色噪聲建模為一階AR模型[23-24]，但ARIMA的有色噪聲數(shù)學模型是明確的，對噪聲再進行時間序列建模并不完全準確。在標準Harvey法濾波方程中[25]，針對如下所示ARIMA模型

(4)

令y1,k=xk，有

a1y1,k-1+y2,k-1+εk

(5)

考慮y2,k

a2y1,k-1+y3,k-1+b1εk

(6)

遞推可得如下關系

(7)

式中，bi=0(i=q+1,q+2,…,p-1)；m=1,2，…,p-1。構建狀態(tài)方程為

(8)

式中，F(xiàn)為狀態(tài)轉(zhuǎn)移矩陣；G為噪聲驅(qū)動矩陣?？梢娛?8)僅適用于p>q時，且建模時，沒有考慮常數(shù)項c。對于本文所述的q=p，且含常數(shù)項的光纖陀螺隨機噪聲ARIMA模型，需要對Harvey方程進行擴展。考慮常數(shù)項c，對于式(7)遞推至m=p-1時，有

(9)

繼續(xù)進行遞推

(10)

(11)

式中

(12)

式中，εk和Vk均為白噪聲且兩者不相關，滿足隨機模型

(13)

擴展后的Harvey方程多出了一維由“常值+白噪聲”構成的隨機常值狀態(tài)參數(shù)。實現(xiàn)了附帶常數(shù)項，且q=p的ARIMA模型狀態(tài)方程的構建。通過Harvey法構造的狀態(tài)噪聲為高斯白噪聲，滿足標準Kalman濾波條件。

3 自適應Kalman濾波

3.1 ARIMA模型及狀態(tài)噪聲在線更新

光纖陀螺隨機噪聲復雜且時變，實際應用中使用有限數(shù)量或不同時期的數(shù)據(jù)所得到的ARIMA模型隨時間的推移，誤差會逐漸增大。因此，有必要依據(jù)當前觀測數(shù)據(jù)進行ARIMA參數(shù)在線更新。將式(1)改寫為

(14)

式中

(15)

使用序貫平差實現(xiàn)模型參數(shù)在線更新的主要流程如下

(16)

(17)

3.2 變分貝葉斯自適應濾波

基于Allan方差的濾波器是一種帶通濾波器，可直接濾除部分低頻系統(tǒng)噪聲的特性，傳統(tǒng)上，將Allan方差值近似為寬頻量測噪聲的方差[1]，選擇取樣間隔為最短采樣時間的Allan方差為

(18)

對于上文所述的量測模型，量測值Zk：N(xk,R)，其中N(μ,∑)表示均值為μ和協(xié)方差陣為∑的高斯概率密度函數(shù)，由式(18)可得

(19)

(20)

可見，當且僅當序列真值X=[X1X2…Xk]為均值時，即量測序列為均值+白噪聲時，偏差項等于0，最小時間間隔的Allan方差才與真實噪聲方差相等，而對于時變序列如本文所建模的ARIMA模型序列，將會存在上述偏差項。

VBAKF是一種近似的貝葉斯法，它利用已知的模型信息、觀測信息和先驗信息來獲得狀態(tài)向量和未知參數(shù)聯(lián)合后驗密度的近似解。可同時修正Qk與Rk的VBAKF的主要流程[15]如下。

(1) 狀態(tài)更新過程

(21)

(22)

式中，IW分布的兩參數(shù)分別為自由度參數(shù)和逆尺度矩陣。

(2) 變分量測迭代更新過程。

返回②繼續(xù)迭代。

4 試驗與分析

對圖1所示的光纖陀螺隨機噪聲進行試驗分析，原始陀螺隨機噪聲Allan方差曲線如圖3所示，可見該陀螺隨機漂移主要由角度隨機游走、零偏不穩(wěn)定性組成。圖4為通過序貫平差建立的ARIMA系數(shù)在線更新，在試驗條件下，光纖陀螺隨機噪聲序列的ARIMA模型較為穩(wěn)定，估值很快收斂，因此，在實際測量中，可僅于初始階段進行模型參數(shù)更新，以提升計算效率。

圖3 光纖陀螺隨機噪聲Allan方差Fig.3 Allan-variance of FOG random noise

圖4 光纖陀螺隨機噪聲ARIMA系數(shù)更新Fig.4 ARIMA coefficients update for FOG random noise

將動態(tài)Allan方差應用于量測噪聲自適應估計是一種近似方法。為對比DAVAR與VBAKF對量測噪聲的估計能力，分別在原始訓練數(shù)據(jù)中，增加均方差為1×10-3、1.5×10-3、2×10-3(°)/s的高斯白噪聲進行解算，量測噪聲估計結果如圖5所示，應認為原始訓練數(shù)據(jù)與所建ARIMA模型相符性較好，因此，量測噪聲與添加噪聲大小近似相同，在圖5中第1000歷元，VBAKF量測噪聲估值分別收斂至1.05×10-3、1.44×10-3、1.95×10-3(°)/s，DAVAR量測噪聲估值收斂至1.46×10-3、1.84×10-3、2.36×10-3(°)/s。因此，可認為VBAKF較DAVAR對量測噪聲有更好的估計性能。使用Harvey法構造Kalman濾波方程，結合VBAKF(以Harvey VBAKF代替)對光纖陀螺隨機噪聲進行實時濾波處理。同時，分別使用不考慮有色噪聲的濾波方程(以KF代替)及不使用VBAKF的Harvey法濾波方程(以Harvey代替)對同一數(shù)據(jù)進行處理，原始數(shù)據(jù)與濾波結果如圖6所示，數(shù)據(jù)統(tǒng)計特性見表2。

圖5 DAVAR與VBAKF噪聲自適應對比Fig.5 Adaptive comparison of DAVAR and VBAKF measurement noise

表2 濾波結果統(tǒng)計特性對比Tab.2 Comparison of statistical results for filtering data (°)/s

圖6 濾波輸出數(shù)據(jù)比較Fig.6 Comparison of filtering output data

比較圖1、圖6可知，總體上，光纖陀螺隨機噪聲經(jīng)Kalman濾波后明顯減小。3種濾波策略結果的均值與原始數(shù)據(jù)相同(表2)，文獻[5]使用滑動平均消除常值結果與原始數(shù)據(jù)均值相差在10-5量級，因此認為本文濾波結果的均值沒有變化。均方差上，相較原始數(shù)據(jù)，傳統(tǒng)Kalman濾波、Harvey法濾波、Harvey法結合VBAKF濾波分別降低了29%、41%、54%，隨機噪聲顯著減小，可見在保證無偏估計的同時，Harvey法建立的Kalman濾波方程有效消除了ARIMA模型濾波中有色噪聲影響，VBAKF實現(xiàn)了對ARIMA模型與狀態(tài)噪聲的解耦，可以有效估計狀態(tài)噪聲與量測噪聲，它對光纖陀螺隨機噪聲抑制效果更加明顯。

5 結 論

使用Harvey法構造的光纖陀螺ARIMA濾波模型可以有效減弱有色噪聲影響，VBAKF具有良好的收斂性和準確的噪聲估計精度。本文將Harvey法構造ARIMA濾波方差引入光纖陀螺隨機誤差降噪，并針對其中存在的狀態(tài)噪聲耦合及先驗信息不準確問題，使用變分貝葉斯自適應卡爾曼濾波實現(xiàn)對狀態(tài)噪聲和量測噪聲的同時修正，并且驗證了VBAKF對量測噪聲有比動態(tài)Allan方差更準確的估計精度。該方法在光纖陀螺隨機誤差建模與補償中可提供一定的參考。