光纖陀螺隨機漂移的補償方法研究

汪徐勝,張瑞民,信 東,李 冰

摘 要:在研究高精度光纖陀螺時,尤其是在捷聯(lián)慣性導航系統(tǒng)中,隨機誤差是光纖陀螺誤差中不可忽視的部分,對光纖陀螺隨機誤差的補償就顯得非常必要。這里基于對光纖陀螺隨機漂移建模的方法,首先采用ARIMA方法對光纖陀螺儀隨機漂移進行建模;然后采用強跟蹤卡爾曼濾波器進行濾波補償,并利用實測數(shù)據(jù)進行了實驗驗證。實驗結果證明,這種方法能夠較好地補償光纖陀螺的隨機漂移。

關鍵詞:光纖陀螺;捷聯(lián)慣性導航系統(tǒng);隨機誤差;卡爾曼濾波器;自適應濾波器

中圖分類號:TP274文獻標識碼:A

文章編號:1004-373X(2009)12-095-04

Research of Compensation Method for Stochastic Shifting of Fiber Optic Gyros

WANG Xusheng1,ZHANG Ruimin1,XIN Dong1,3,LI Bing2

(1.The Second Artillery Officer School,Qingzhou,262500,China;

2.The Second Artillery Officer Stationed in the Forth Space Flight Academe,Xi′an,710025,China;

3.The Second Artillery Engineering College,Xi′an,710025,China)

Abstract:When doing research in optical fiber gyros,especially in strapdown inertial navigation system,stochastic error is an assignable part of the errors of optical fiber gyros.It is very necessary to compensate stochastic error of optical fiber gyros.Basing on the first method,the article builds a model for the stochastic shifting,then adopts strong tracing Kalman filter to filter and compensate.At last,the methods are testified by experimental data.The results prove that the method is effective.

Keywords:optical fiber gyros;strapdown inertial navigation system;stochastic error;Kalman filter;self-adapting filter

0 引 言

在慣性技術中,通常將慣性敏感元件的誤差模型分為靜態(tài)誤差模型、動態(tài)誤差模型和隨機誤差模型三類[1,2]。靜態(tài)誤差模型和動態(tài)誤差模型可以比較容易地用代數(shù)方程來表示,因此比較容易補償。隨機誤差是指誤差中隨機變化的部分,要用統(tǒng)計規(guī)律描述,誤差要通過濾波方法來減小。研究低精度光纖陀螺時,由于陀螺的隨機漂移較大,嚴重影響了陀螺的測試精度,必須進行隨機誤差補償。光纖陀螺應用于慣性導航中,由于工作時間比較長,積分結果將會對導彈造成不可容忍的誤差。因此,對光纖陀螺的漂移進行補償具有重要的意義。

人們對這一部分的研究也相對比較成熟[3-5],但在建立隨機誤差模型時,大都局限于采用簡單的自回歸模型,影響了誤差模型的精度;同時,濾波大都采用普通kalman濾波方法,由于這種濾波方法的局限性(噪聲模型必須準確已知),容易導致濾波的不穩(wěn)定甚至發(fā)散,本文采用能夠跟蹤信號變化的強跟蹤kalman濾波方法[6],有效地克服了這一問題。

1 時間序列建模方法研究

1.1 建模前的數(shù)據(jù)檢驗

一般,光纖陀螺的漂移數(shù)據(jù)可以近似看作時間序列來處理,而常用的時間序列建模方法[7]有自回歸模型法(AR)、滑動模型法(MA)、自回歸滑動模型法(ARMA)。通常ARMA模型更能夠準確地反映數(shù)據(jù)的隨機特性。無論在用那種方法進行建模,首先時間序列必須是平穩(wěn)的。但實際情況是,由于受各種干擾的影響,所測得的時間序列一般不是平穩(wěn)的,這就需要進行數(shù)據(jù)的檢驗和平穩(wěn)化處理。主要包括平穩(wěn)性檢驗、正態(tài)性檢驗以及零均值檢驗[8,9],在建模前,要對數(shù)據(jù)進行這3種檢驗。由于零均值檢驗和正態(tài)性檢驗在參考文獻中有詳細的介紹,文中僅對平穩(wěn)性檢驗進行介紹。

時間序列的平穩(wěn)性是建模的重要前提,因此需要對測試數(shù)據(jù)進行平穩(wěn)性檢驗。這里采用一種工程中常用的逆序檢驗法,其具體過程為:

設樣本序列x1,x2,…,xN足夠長,把樣本序列分成k個子序列,即N=lM,M是一個較大的正整數(shù);l也為正整數(shù)。

逆序檢驗法就是檢驗各子列均值的差異性,把各子列均值μj構成一個序列μ1,μ2,…,μl。當i>j時,每出現(xiàn)一次μi>μj,定義Aj為μj的一個逆序數(shù)。對于某一μj,定義μj的逆序數(shù)Aj為μi>μj(i>j)出現(xiàn)的次數(shù),所有逆序數(shù)Aj的總和為序列的逆序總數(shù),其表達式為:

A=∑lj=1Aj

由于{xk}是隨機過程的一個樣本,所以均值μj也應該是隨機的。當{xk}是平穩(wěn)序列時,μ1后面的l-1個隨機數(shù)μj大于和小于μ1的概率是相等的,故μ1的逆序數(shù)A1的理論平均值E[A1]=(l-1)/2;同理,E[A2]=(l-2)/2,…,E[Ak]=(l-k)/2(k=1,2,…,l-1),從而逆序數(shù)的理論平均值:

E[A]=∑l-1j=1EAj=l(l-1)/4

逆序總數(shù)的理論方差為:

σ2A=l(2l2+3l-5)/72

構造統(tǒng)計量:

u=(A+0.5-E[A])/σA

這樣構造的統(tǒng)計量滿足標準正態(tài)分布,因此當顯著性水平為0.05時,如果|u|≤1.96,則可確定μj間無顯著差異,即序列是平穩(wěn)的。

1.2 ARIMA時間序列建模方法

在許多實際問題中,所考慮的時間序列{xk}并不是平穩(wěn)的,它不能用ARMA模型來表示。但是,如果將{xk} 進行有限次差分后,得到的時間序列{x′k} 是平穩(wěn)的,{x′k}可以用ARMA模型來描述。原時間序列{xk}可以用ARIMA(n,d,m)(基于序列平穩(wěn)化的ARMA建模)表示。平穩(wěn)序列的ARMA(n,m)模型表達式為:

Φ(B)xk=Θ(B)ak(1)

式中:Φ(B)=1-∑ni=1ΦiBi;Θ(B)=1-∑mj=1θjBj;n是AR部分的階次;m是MA部分的階次;B是后移算子,Bxk=xk-1;ak是殘差,滿足ak～(0,σ2a)。

如果序列不是平穩(wěn)的,則要對序列進行d階差分使序列平穩(wěn)化。差分是將時間序列逐項遞減而得到一個新的序列,其目的是為了消除原始序列的趨勢項,將非平穩(wěn)序列化為平穩(wěn)序列。然后建立差分后序列的ARMA(n,m)模型,原序列則可以表示為ARIMA(n,d,m)模型:

Φ(B)(1-B)dxk=Θ(B)ak(2)

通過以上分析,建立ARIMA(n,d,m)是建立在ARMA(n,m)模型基礎之上的。在建模前要進行數(shù)據(jù)平穩(wěn)化的檢驗,如果隨機序列是平穩(wěn)的,取d=0。ARIMA(n,d,m)模型簡化為ARMA(n,m)。

1.3 建立ARMA模型

1.3.1 擬合模型

按照ARMA模型參數(shù)估計中采用的方法和理論,模型參數(shù)估計方法分為3種[2]。由于時序理論估計法算比較簡單,而且可以保證較高的精度,所以得到比較廣泛的應用。在此采用的也是這種參數(shù)估計方法。

ARMA(n,m)模型[7]可以表示為:

xk=∑ni=1Φixk-1-∑mj=1θjak-j+ak(3)

對自回歸模型參數(shù)Φi和滑動回歸模型參數(shù)θj可以采用長自回歸模型法。這種方法的思路是:基于觀測序列建立起來的模型AR(n),MA(m),ARMA(n,m)都是等價系統(tǒng)的數(shù)學模型,具有相同的傳遞函數(shù)?；谶@個原則,可以先估計出AR模型,然后利用傳遞函數(shù)的等價性,估計出ARMA的參數(shù)自回歸模型參數(shù)Φi和滑動回歸模型參數(shù)θj。下面就長自回歸模型法[3,5]做簡單介紹。

隨機序列{xk}可以用模型AR(p)表示,也可以用ARMA(n,m)表示,AR(p)描述的等價系統(tǒng)傳遞函數(shù)為:

1/Φp(B)=1/(1-∑pi=1IiBi)(4)

式中:B為后移算子,Ii是逆函數(shù),其取值等于AR(p)模型的參數(shù)Φi(i=1,2,…,p)。用ARMA(n,m)描述的等價系統(tǒng)傳遞函數(shù)為:

θm(B)/Φn(B)=(1-∑mj=1θjBj)/(1-∑ni=1ΦiBi)(5)

由于不同形式的傳遞函數(shù)描述的系統(tǒng)是等價的,所以上述兩式應該相等。把兩式展開,得到等式:

(1-I1B-I2B2-…-IpBp)(1-θ1B-θ2B2-

…-θmBm)=1-Φ1B-Φ2B2-…-ΦnBn

所以,對應B的相同階次的系數(shù)應該相等。由此得到:

Φ1=θ1+I1

Φ2=θ2-θ1I1+I2

Φ3=θ3-θ2I1-θ1I2+I3

…

Φn=-θmIn-m-θm-1In-m-1-…-θ1In-1+In

0=-θmIk-m-θm-1Ik-m-1-…-θ1Ik-1+Ik,

k>n(6)

用最后一個方程,當k取n+1,n+2,…,n+m的時候,可以用以下矩陣求出θj(j=1,2,…,m):

In+1In+2驣n+m=InIn-1In-2…In-m+1In+1InIn-1…In-m-2

In+m-1In+m-2In+m-3…Inθ1θ2螃萴(7)

在求出θj(j=1,2,…,m)以后,利用式(1)可以求出Φi(i=1,2,…,n)。關于Φi(i=1,2,…,n)的線性方程為:

Φ1Φ2螃祅=θ1θ2螃萵+100…0-θ110…0

-θn-θn-1-θn-2…0I1I2驣n(8)

求解上述線性方程,便得到Φi(i=1,2,…,n)的值,這就是長自回歸模型法的計算原理。長自回歸模型法在估計Φi和θj的過程中都是解線性方程組。它將ARMA(n,m)模型參數(shù)估計的非線性問題轉化為線性問題來解決,大大簡化了問題。

2 干涉式光纖陀螺儀隨機漂移的建模

這里以一組VG951型光纖陀螺儀實測數(shù)據(jù)為例,建立光纖陀螺儀的隨機漂移模型,這組數(shù)據(jù)采樣間隔為1 s,總體樣本長度為3 000,數(shù)據(jù)中已減去了地球常值的影響。圖1為實測數(shù)據(jù)曲線。

圖1 光纖陀螺儀的實測數(shù)據(jù)

由圖1可以看出序列均值不為零,且有下降的趨勢。用逆序法檢驗平穩(wěn)性,取l=6,得到u=-2.630 1,因為|u|>1.96,所以序列不是平穩(wěn)的。其平均值為3.846 2,漂移值為2.081 6,因此建模前要對序列進行預處理。

首先對原序列進行一階差分,得到一新的序列,檢驗新序列的統(tǒng)計特性,得到如下數(shù)據(jù):

平穩(wěn)性檢驗:u=0.751 5;

正態(tài)性檢驗:ξ=0.014 3;γ=3.004 1;

經(jīng)過檢驗,處理后的數(shù)據(jù)基本上符合平穩(wěn)性和正態(tài)性的檢驗,可以對其建立ARMA(n,m)模型。

分析處理后的序列,其可以用ARMA(n,m)模型表示。建模采用BIC準則(見參考文獻[5])確定模型的階次,用長自回歸模型法辨識模型的系數(shù)。經(jīng)計算ARMA(3,1)模型為適用模型,建模過程為:

首先確定AR(3)模型的系數(shù),模型可以表示為(其中v(k)為白噪聲):

y(k)=a1y(k-1)+a2y(k-2)+a3y(k-3)+

(1-a1-a2-a3)+v(k)

求解關于模型參數(shù)的矩陣方程,得系數(shù)值:

123=-0.760 4-0.521 8-0.257 5

從而得到ARMA(2,1)的系數(shù):

θ1=0.493 5;Φ1=-0.267 0;Φ2=-0.146 6

代入差分方程,得到ARIMA模型的系數(shù):

Φ=0.7330.120 40.146 6;θ=[0.499 2]

從而,原隨機序列模型的表達式為:

xk=0.733xk-1+0.120 4xk-2+0.146 6xk-3-

ak+0.493 5ak-1(9)

3 濾波算法

強跟蹤卡爾曼濾波是一種適合計算機計算的遞推濾波方法,它成功地采用了狀態(tài)空間概念,把信號過程視為白噪聲作用下一個線性系統(tǒng)的輸出。這種方法不需要存儲過去的觀測數(shù)據(jù),當新的數(shù)據(jù)被觀測后,只要根據(jù)新的數(shù)據(jù)與前一個時刻的估計量,再借助于信號過程本身的狀態(tài)轉移方程,按照遞推公式,即可計算出該時刻的估計量。

光纖陀螺輸入/輸出數(shù)據(jù)的卡爾曼濾波器的設計流程,如圖2所示。

圖2 卡爾曼濾波的總體流程

離散系統(tǒng)的卡爾曼最優(yōu)濾波的遞推方程為:

(k|k)=(k|k-1)+

K(k)[Z(k)-H(k)(k|k-1)]

(k|k-1)=Φ(k,k-1)(k-1|k-1)

K(k)=P(k|k-1)HT(k)?

[H(k)P(k|k-1)HT(k)+Rk]-1

P(k|k-1)=Φ(k,k-1)P(k-1|k-1)ΦT?

(k,k-1)+Γ(k,k-1)Qk-1ΓT(k,k-1)

P(k|k)=[I-K(k)H(k)]P(k|k-1)?

[I-K(k)H(k)]T+K(k)RkKT(k)

對第2部分隨機序列建模的結果,可以直接運用卡爾曼濾波方法減小隨機漂移,但是這種方法不能跟蹤變化的信號,在仿真中給出的是一種能夠跟蹤變化信號的濾波方法。

4 仿真驗證

此前建立的ARIMA模型可作為卡爾曼濾波器的狀態(tài)方程,建立觀測方程得到隨機漂移的觀測數(shù)據(jù)。目前,得到的數(shù)據(jù)只有光纖陀螺輸出值?？紤]到光纖陀螺的輸出由常值漂移、角速度敏感量、隨機漂移和噪聲組成,可以假設相鄰兩時刻的角加速度變化率相等。光纖陀螺輸入/輸出的關系式可以表示為:

w=(1/ks)(y-y0-x-ε)

對上式進行二階微分可以得到:

d2wdt2=1ks(-d2ydt2+d2xdt2+d2εdt2)=0

差分后得到結果:

yk-2yk-1+yk-2=xk-2xk-1+xk-2+

εk-2εk-1+εk-2

設:

νk=εk-2εk-1+εk-2

令Zk=yk-2yk-1+yk-2作為隨機噪聲的觀測量,可以寫出觀測方程:

Zk=xk-2xk-1+xk-2+νk

利用上式作為觀測方程,與前面所得到的模型方程一起構成典型卡爾曼濾波器。

令Xk=[xk,xk-1,xk-2],將所建隨機模型(1)和觀測方程寫成矩陣形式為:

Xk=ΦXk-1+GWkZk=HXk+Vk(10)

式中,Φ系統(tǒng)狀態(tài)轉移矩陣,其值為:

Φ=0.7330.120 40.146 6100010

G是系統(tǒng)噪聲陣,可表示為:

G=1-0.493 50000000

H是系統(tǒng)觀測陣,表示為:H=[1-21]。

式中:Wk是系統(tǒng)噪聲;Vk是觀測噪聲,兩者是互不相關的零均值高斯白噪聲,在濾波時取為標量。圖3為補償后的光纖陀螺隨機漂移數(shù)據(jù)。

圖3 補償后的隨機漂移

經(jīng)補償,序列的隨機漂移為0.745 5 °/h。可以看出,依據(jù)漂移的模型進行濾波能夠有效地抑制隨機漂移。這里的關鍵問題是得到光纖陀螺的隨機漂移模型。但在一些實際問題中,光纖陀螺的漂移模型是時變的,事先建立的漂移模型不一定能夠準確反映當前的陀螺漂移。因此,在實時補償時,會存在一定的危險。為了達到較好的補償效果,需要在線建立光纖陀螺的隨機漂移模型,這樣才能達到較好的動態(tài)補償效果。分析影響光纖陀螺隨機漂移的各種因素,采用自適應建模的方法,建立光纖陀螺的自適應補償模型,根據(jù)檢測的環(huán)境變量和采樣數(shù)據(jù)的統(tǒng)計特性自動修改模型參數(shù),能夠有效地克服光纖陀螺漂移模型的時變性,可以達到較理想的效果,其中采用橫向自適應濾波器是一種較好的解決方案。

5 結 語

靜止時光纖陀螺輸出的是一隨機序列,在分析了光纖陀螺隨機漂移特性以及時間序列的隨機建模的基礎上,利用實測數(shù)據(jù)建立了光纖陀螺的ARIMA漂移模型。然后利用強跟蹤卡爾曼濾波進行隨機漂移的最優(yōu)估計,用以對輸出進行補償。對檢測數(shù)據(jù)的處理結果證明,依據(jù)ARIMA模型對光纖陀螺的隨機漂移進行補償具有較好的效果,大大降低了光纖陀螺的隨機漂移。采用的ARIMA建模方法以及強跟蹤Kalman濾波方法,不管是在精度上,還是在算法的簡單性和可行性上,都比以前的隨機建模及濾波方法有一定的提高,這也是本文的創(chuàng)新點所在。

參考文獻

[1]王珍熙.捷聯(lián)式慣性導航系統(tǒng)慣性元件的設置與可靠性[J].中國慣性技術學報,1996,4(1):61-65.

[2]郭秀中.慣導系統(tǒng)陀螺儀理論[M].北京:國防工業(yè)出版社,1996.

[3]王新龍.光纖陀螺隨機誤差模型分析[J].北京航空航天大學學報,2006,32(7):769-772.

[4]吉訓生,王壽榮.MEMS陀螺儀隨機漂移誤差研究[J].宇航學報,2006,27(4):74-76.

[5]李愛華.閉環(huán)光纖陀螺的建模方法[J].電光與控制,2006,13(4):43-45.

[6]劉銘,周東華.殘差歸一化的強跟蹤濾波器及其應用[J].中國電機工程學報,2005,25(2):74-78.

[7]王志賢.最優(yōu)估計與系統(tǒng)辨識[M].西安:西北工業(yè)大學出版社,2003.

[8]秦永元.卡爾曼濾波與組合導航原理[M].西安:西北工業(yè)大學出版社,2004.

[9]陸光華.隨機信號處理[M].西安:西安電子科技大學出版社,2003.