本發(fā)明涉及一種基于自適應(yīng)Kalman濾波的陀螺隨機(jī)噪聲AR建模方法，屬于捷聯(lián)慣性導(dǎo)航
技術(shù)領(lǐng)域：
。
背景技術(shù)：
：陀螺儀輸出誤差的建模和補(bǔ)償是提高捷聯(lián)慣性導(dǎo)航系統(tǒng)精度的重要方法。陀螺儀的輸出誤差可以分為確定性誤差和隨機(jī)誤差。確定性誤差可以通過實(shí)驗(yàn)標(biāo)定的方法予以大部分消除。對陀螺儀隨機(jī)誤差的處理一般是從信號處理角度，通過數(shù)學(xué)工具進(jìn)行建模和補(bǔ)償，其中時(shí)間序列建模方法得到了廣泛的應(yīng)用。常用的時(shí)間序列模型有ARMA(AutoRegressiveMovingAverage，自回歸滑動平均模型)、AR(AutoRegressive，自回歸模型)和MA(MovingAverage，滑動平均模型)三種。由Wold定理和Kolmogorov-Szego定理可知，這三種模型可以相互轉(zhuǎn)化，都具有普遍適用性。但如果選擇的模型合適，就可以用低階的模型獲得相對高的建模精度，從而使系統(tǒng)更加簡單。實(shí)際工程中不同型號和批次的陀螺儀其隨機(jī)噪聲特性并不相同。有些型號批次陀螺儀的噪聲數(shù)據(jù)其自相關(guān)和偏相關(guān)特性采用AR模型建模更加適合。傳統(tǒng)陀螺隨機(jī)噪聲AR建模方法存在著算法收斂速度慢，所需樣本多等缺點(diǎn)。技術(shù)實(shí)現(xiàn)要素：為解決傳統(tǒng)陀螺隨機(jī)噪聲建模方法的不足，本發(fā)明提出一種基于自適應(yīng)Kalman濾波技術(shù)的陀螺隨機(jī)噪聲AR模型建模方法。相比傳統(tǒng)方法，本發(fā)明方法建模速度快，所需樣本少，建模精度高。本發(fā)明采用的技術(shù)方案為：一種基于自適應(yīng)Kalman濾波技術(shù)的陀螺隨機(jī)噪聲建模方法，包括以下步驟：(1)首先對陀螺儀隨機(jī)噪聲進(jìn)行平穩(wěn)隨機(jī)檢驗(yàn)。采用數(shù)字微分方法保證隨 機(jī)噪聲數(shù)據(jù)滿足平穩(wěn)隨機(jī)過程的要求。(2)將AR模型的參數(shù)做為系統(tǒng)狀態(tài)量。以AR(2)模型為例：z(k)＝a1z(k-1)+a2z(k-2)+ε(k)；式中a1、a2是AR模型的自回歸系數(shù)；z(k-1)是一階自回歸項(xiàng)；z(k-2)是二階自回歸項(xiàng)；ε(k)是均值為0零，方差未知(常數(shù))的白噪聲；a1，a2是待定的模型參數(shù)。將a1，a2做為系統(tǒng)的狀態(tài)量X：X＝[a1(k)，a2(k)]T；(3)建立系統(tǒng)觀測方程：z(k)＝H(k)X(k)+v(k)；式中H(k)是系統(tǒng)觀測陣：H(k)＝[z(k-1)，z(k-2)]；v(k)是系統(tǒng)觀測噪聲，它主要由AR(2)模型的白噪聲ε(k)構(gòu)成。(4)建立系統(tǒng)狀態(tài)方程：X(k+1)＝X(k)；(5)采用未知定常的觀測噪聲統(tǒng)計(jì)估值器估計(jì)觀測噪聲的均值和方差，在此基礎(chǔ)上運(yùn)用自適應(yīng)Kalman濾波器對系統(tǒng)狀態(tài)量X，也即AR模型的參數(shù)進(jìn)行估計(jì)。進(jìn)一步，上述步驟(1)所述數(shù)據(jù)預(yù)處理包括以下步驟：(a)周期性數(shù)據(jù)檢驗(yàn)和濾除；(b)平穩(wěn)隨機(jī)檢驗(yàn)；(c)分析數(shù)據(jù)的自相關(guān)和偏相關(guān)特性；當(dāng)自相關(guān)系數(shù)函數(shù)呈現(xiàn)拖尾性而偏相關(guān)系數(shù)函數(shù)呈現(xiàn)截尾性時(shí)，則選擇AR模型建模。且偏相關(guān)系數(shù)函數(shù)截?cái)嗵幘褪茿R模型的階數(shù)。且步驟(b)所述的平穩(wěn)隨機(jī)檢驗(yàn)還包括差分處理過程。進(jìn)一步，用下標(biāo)“k”表示Kalman濾波的迭代次數(shù)，則步驟(5)所述的自適應(yīng)Kalman濾波器的迭代過程如下所示：狀態(tài)一步預(yù)測：x^k|k-1=Φk|k-1x^k-1|k-1=x^k-1|k-1;]]>式中是Kalman濾波器根據(jù)：k-1迭代時(shí)得到的系統(tǒng)狀態(tài)量X的估計(jì)值對第k次迭代時(shí)系統(tǒng)狀態(tài)量X估計(jì)值的預(yù)測，稱為狀態(tài)一步預(yù)測值；Φk|k-1是系統(tǒng)狀態(tài)轉(zhuǎn)移陣，當(dāng)對AR(2)進(jìn)行參數(shù)估計(jì)時(shí)，它為2×2維的單位矩陣；是k-1次迭代時(shí)得到的狀態(tài)量X的估計(jì)值；一步預(yù)測誤差方差陣：Pk|k-1=Φk|k-1Pk-1|k-1Φk|k-1T+Γk|k-1Qk-1Γk|k-1T=Pk-1|k-1;]]>式中Pk|k-1是一步預(yù)測誤差方差陣；Q是系統(tǒng)激勵噪聲的方差，Q＝0；Pk-1|k-1是k-1次迭代時(shí)得到的估計(jì)誤差方差陣；濾波增益矩陣：Kk=Pk|k-1HkT[hkPk|k-1HkT+R^k-1]-1]]>式中Kk是Kalman濾波器第k次迭代時(shí)的濾波增益矩陣；Pk|k-1是一步預(yù)測誤差方差陣；Hk是第k次迭代時(shí)的系統(tǒng)觀測矩陣；是第k-1次迭代時(shí)得到的系統(tǒng)觀測噪聲方差的估計(jì)；新息：ϵk=zk-Hkx^k|k-1-r^k-1;]]>式中是第k次迭代時(shí)的新息；zk是第k次迭代時(shí)的系統(tǒng)觀測量，也即陀螺隨機(jī)噪聲的采樣；Hk是第k次迭代時(shí)的系統(tǒng)觀測矩陣；是狀態(tài)一步預(yù)測；是第k-1次迭代時(shí)得到的系統(tǒng)觀測噪聲均值的估計(jì)；狀態(tài)估計(jì)：x^k|k=x^k|k-1+Kkϵk;]]>式中是Kalman濾波在迭代k次時(shí)得到的狀態(tài)量X的估計(jì)值；是狀態(tài)一步預(yù)測；Kk是第k次迭代時(shí)的濾波增益矩陣；是第k次迭代時(shí)的新息；估計(jì)誤差方差陣：Pk|k＝[In-KkHk]Pk|k-1；式中Pk|k是第k次迭代時(shí)的估計(jì)誤差方差陣；In是n×n維的單位矩陣，當(dāng)對AR(2)進(jìn)行參數(shù)估計(jì)時(shí)，n＝2；Kk是第k次迭代時(shí)的濾波增益矩陣；Hk是第k次迭代時(shí)的系統(tǒng)觀測矩陣；Pk|k-1是一步預(yù)測誤差方差陣；未知定常的觀測噪聲均值估計(jì)：r^k=(1-1k)r^k-1+1k(zk-Hkx^k|k-1);]]>式中為濾波器第k次迭代時(shí)觀測噪聲均值的估計(jì)值；zk是第k次迭代時(shí)的系統(tǒng)觀測量；是第k次迭代時(shí)狀態(tài)量的一步預(yù)測，Hk是第k次迭代時(shí)的系統(tǒng)觀測陣；未知定常的觀測噪聲方差估計(jì)：R^k=(1-1k)R^k-1+1k(ϵkϵkT-HkPk|k-1HkT);]]>式中為濾波器第k次迭代時(shí)觀測噪聲方差的估計(jì)值；是第k次迭代時(shí)的新息；Pk|k-1是第k次迭代時(shí)的一步預(yù)測誤差方差陣，Hk是第k次迭代時(shí)的系統(tǒng)觀測陣。只要AR模型階數(shù)選取合適，通過自適應(yīng)Kalman濾波器迭代，則可迅速得到AR模型參數(shù)的最優(yōu)估計(jì)。本發(fā)明的有益之處在于：本發(fā)明可以有效減少采樣樣本數(shù)和采樣時(shí)間，建模速度快。當(dāng)新的噪聲樣本到來時(shí)，自適應(yīng)Kalman濾波器能夠保證建立的AR模型參數(shù)可以及時(shí)進(jìn)行更新，具有快速實(shí)時(shí)建模的優(yōu)點(diǎn)。附圖說明圖1為某國產(chǎn)光纖陀螺的原始噪聲數(shù)據(jù)；圖2為該陀螺原始噪聲經(jīng)一階差分后滿足平穩(wěn)隨機(jī)過程要求的數(shù)據(jù)；圖3為該陀螺隨機(jī)噪聲數(shù)據(jù)的自相關(guān)特性分析；圖4為該陀螺隨機(jī)噪聲數(shù)據(jù)的偏相關(guān)特性分析；圖5為發(fā)明方法對該陀螺隨機(jī)噪聲進(jìn)行AR(2)建模結(jié)果a1；圖6為發(fā)明方法對該陀螺隨機(jī)噪聲進(jìn)行AR(2)建模結(jié)果a2；圖7是為本發(fā)明方法步驟過程示意圖。具體實(shí)施方式以下通過具體實(shí)施例對本發(fā)明作具體的介紹。下面以對某國產(chǎn)光纖陀螺儀輸出隨機(jī)噪聲的AR(2)建模為例，首先對陀螺隨機(jī)噪聲的采樣數(shù)據(jù)進(jìn)行預(yù)處理，保證其滿足平穩(wěn)隨機(jī)過程的要求。數(shù)據(jù)預(yù)處理包括3步：1)周期性數(shù)據(jù)檢驗(yàn)和濾除；2)平穩(wěn)隨機(jī)檢驗(yàn)；3)根據(jù)自相關(guān)和偏相關(guān)特性選擇合適的時(shí)間序列模型建模。(a)周期性數(shù)據(jù)檢驗(yàn)和濾除如果信號的采樣含有周期性或準(zhǔn)周期性的數(shù)據(jù)，則在功率譜上會出現(xiàn)尖峰，尖峰處對應(yīng)的角頻率也即周期項(xiàng)的角頻率，而不含周期項(xiàng)的隨機(jī)數(shù)據(jù)的功率譜沒有明顯尖峰。信號功率譜和自相關(guān)函數(shù)的定義如下：Sx(ω)=limN→∞12N+1|X(ejω,N)|2,Rx(m)=1NΣn=0N-m-1x(n)x(n+m)---(1)]]>如果陀螺的噪聲數(shù)據(jù)含有周期為T的周期項(xiàng)，則使用T步差分是去除周期項(xiàng)的一種有效辦法。(b)平穩(wěn)隨機(jī)檢驗(yàn)經(jīng)過趨勢項(xiàng)和周期項(xiàng)檢驗(yàn)和濾除后的陀螺噪聲數(shù)據(jù)可以認(rèn)為是隨機(jī)序列。但還需檢驗(yàn)其平穩(wěn)性。如不滿足平穩(wěn)隨機(jī)序列的要求，也需通過差分處理的方法使之平穩(wěn)化。實(shí)際工程應(yīng)用中一般用非參數(shù)法，例如游程檢驗(yàn)法來檢驗(yàn)時(shí)間序列的平穩(wěn)性。游程是指在保持隨機(jī)序列原有順序的情況下，隨機(jī)序列中具有相同符號的序列。這種檢驗(yàn)法里把觀測值分成兩類：例如和是序列均值。分別用“+”和“-”表示，稱為正負(fù)符號。游程檢驗(yàn)法是以正負(fù)符號的數(shù)目N1和N2做為統(tǒng)計(jì)檢驗(yàn)量，根據(jù)正態(tài)分布表來確定檢驗(yàn)的接受域。例如，以置信度95％為可接受，即假設(shè)檢驗(yàn)的顯著性水平α＝0.05。按照雙邊假設(shè)檢驗(yàn) (2σ)查找標(biāo)準(zhǔn)正態(tài)分布表，則檢驗(yàn)統(tǒng)計(jì)量U在[-1.96，1.96]范圍時(shí)為假設(shè)接受域?？梢哉J(rèn)為該序列是平穩(wěn)的。檢驗(yàn)統(tǒng)計(jì)量的計(jì)算公式為：U=γ-μγσγ---(2)]]>式中U是檢驗(yàn)統(tǒng)計(jì)量；γ是游程數(shù)；μγ是游程的期望數(shù)；σγ是游程的標(biāo)準(zhǔn)差；計(jì)算方法如下：μγ=2N1N2N+1,σγ=[2N1N2(2N1N2-N)N2(N-1)]12,N=N1+N2---(3)]]>式中N1、N2分別是游程中正、負(fù)符號的數(shù)目；μγ是游程的期望數(shù)；σγ是游程的標(biāo)準(zhǔn)差。不管多復(fù)雜的非平穩(wěn)隨機(jī)過程，經(jīng)過一定次數(shù)的差分后，總可以變成平穩(wěn)隨機(jī)過程。而對于陀螺的隨機(jī)噪聲數(shù)據(jù)來說，一般經(jīng)過一階或二階差分處理，即可基本滿足平穩(wěn)性要求。(c)根據(jù)現(xiàn)代時(shí)間序列分析理論，時(shí)間序列模型選擇的依據(jù)在于平穩(wěn)隨機(jī)數(shù)據(jù)的自相關(guān)系數(shù)函數(shù)和偏相關(guān)系數(shù)函數(shù)的截尾和拖尾特性的不同，具體區(qū)分如下表所示：表1三種時(shí)間序列模型的自相關(guān)和偏相關(guān)特性其中，自相關(guān)系數(shù)函數(shù)p(h)和偏自相關(guān)系數(shù)函數(shù)的數(shù)學(xué)定義如下：式中x(k)、x(k-h)分別代表第k和第k-h個(gè)樣本數(shù)據(jù)；“Cov”代表自協(xié)方差函數(shù)；“Var”代表方差函數(shù)；和分別是利用x(n+1)，x(n+2)，…，x(n+k-1)對x(n)和x(n+k)做的最佳線性估計(jì)。附圖1是地面靜態(tài)實(shí)驗(yàn)測得的某國產(chǎn)光纖陀螺的原始噪聲數(shù)據(jù)(x軸，共3600s)。采樣頻率100Hz，單位：°/s，陀螺坐標(biāo)軸指向：x-東，y-北，z-天。首先按對其進(jìn)行數(shù)據(jù)預(yù)處理，通過一階差分，最終得到滿足平穩(wěn)隨機(jī)過程要求的數(shù)據(jù)(見附圖2)。分析其自相關(guān)系數(shù)函數(shù)p(h)和偏相關(guān)系數(shù)函數(shù)定義如式(4)所示。從附圖3、附圖4中可以看出，該陀螺儀隨機(jī)噪聲數(shù)據(jù)的自相關(guān)系數(shù)函數(shù)明顯呈現(xiàn)拖尾性，而偏自相關(guān)系數(shù)函數(shù)呈現(xiàn)截尾特性，且在“2”后截?cái)唷谋?可知應(yīng)優(yōu)先選用AR模型建模。且偏相關(guān)系數(shù)函數(shù)的截?cái)嗵帯?”就是AR模型的階數(shù)。因此選用二階AR模型進(jìn)行建模。傳統(tǒng)AR模型建模方法有最小二乘法、Yule-Walker法等。但這些傳統(tǒng)方法需要的樣本數(shù)較多，算法收斂速度較慢。為了解決傳統(tǒng)陀螺隨機(jī)噪聲AR建模方法中存在的問題，本發(fā)明提出了一種新的基于自適應(yīng)Kalman濾波的陀螺隨機(jī)噪聲AR模型建模方法，如圖7所示。首先對陀螺儀隨機(jī)噪聲進(jìn)行平穩(wěn)隨機(jī)檢驗(yàn)，通過數(shù)字微分保證隨機(jī)噪聲數(shù)據(jù)滿足平穩(wěn)隨機(jī)過程的要求后，將AR模型的參數(shù)做為系統(tǒng)狀態(tài)量；假設(shè)陀螺儀隨機(jī)噪聲的AR模型階數(shù)為AR(2)：z(k)＝a1z(k-1)+a2z(k-2)+ε(k)(5)式中的ε(k)是均值為0方差未知(常數(shù))的白噪聲。取狀態(tài)量：X(k)＝[a1(k)，a2(k)T，則可以把上式轉(zhuǎn)化為系統(tǒng)觀測方程：z(k)＝H(k)X(k)+v(k)(6)式中的H是系統(tǒng)觀測陣：H(k)＝[z(k-1)，z(k-2)]，v(k)＝ε(k)?？紤]到當(dāng)陀螺儀的隨機(jī)噪聲特性相對穩(wěn)定時(shí)，當(dāng)滿足一定的樣本數(shù)后，所建立的AR模型的參數(shù)估計(jì)值應(yīng)當(dāng)穩(wěn)定收斂于真值，不再隨樣本數(shù)量的增加而改變。即式中上標(biāo)“^”表示參數(shù)的估計(jì)值；可以建立系統(tǒng)的狀態(tài)方程：X(k+1)＝X(k)(7)從上式可以看出，系統(tǒng)噪聲w(k)為0，狀態(tài)轉(zhuǎn)移陣Φ＝I。觀測噪聲v(k)主要由AR模型中的白噪聲ε(k)引起。有：Ewk=0,E[wk,wj]=0Evk=rk,E[vk,vj]=Rk,δkjE[wkvjT]=0---(8)]]>式中“E[]”代表均值函數(shù)；rk、Rk是系統(tǒng)觀測噪聲v(k)的均值陣和方差陣；δkj是克羅內(nèi)克函數(shù)。實(shí)際系統(tǒng)由于不可避免存在建模誤差，所以rk不一定為0。又由于AR模型中白噪聲ε(k)的方差特性是未知定常的，因此Rk也是未知定常型的。對于未知定常型的觀測噪聲，通常利用Sage-Husa的次優(yōu)無偏MAP(MaximumAPosteriori極大后驗(yàn))噪聲統(tǒng)計(jì)估值器估計(jì)未知定常觀測噪聲的統(tǒng)計(jì)特性。對應(yīng)的Kalman濾波器稱為自適應(yīng)Kalman濾波。未知定常觀測噪聲統(tǒng)計(jì)估值器如下式：r^k=(1-1k)r^k-1+1k(zk-Hkx^k|k-1)R^k=(1-1k)R^k-1+1k(ϵz,k|k-1ϵz,k|k-1T-HkPz,k|k-1|HkT)---(9)]]>對該光纖陀螺隨機(jī)噪聲數(shù)據(jù)采用發(fā)明提出的自適應(yīng)Kalman濾波方法進(jìn)行建模，濾波初值選取X^(0|0)=[0,0]T,]]>P(0)=10000100,]]>r^k(0)=0,R^k(0)=0.01.]]>建模過程如附圖5-6中的實(shí)曲線所示。圖中虛直線是采用傳統(tǒng)“Y-W”法(1小時(shí)采樣數(shù)據(jù)360000樣本點(diǎn))對該陀螺隨機(jī)噪聲數(shù)據(jù)進(jìn)行AR(2)建模的結(jié)果：z(k)＝-0.671z(k-1)-0.334z(k-2)+ε(k)(10)在濾波迭代時(shí)，可以根據(jù)系統(tǒng)的建模精度要求設(shè)定閾值θ。如有連續(xù)10次運(yùn)行結(jié)果中各參數(shù)估計(jì)值的最大值和最小值之間的差值都小于閾 值θ，則做為迭代的退出條件。這里選取θa1＝θb1＝θb2＝0.001。程序在運(yùn)行980次迭代后退出，狀態(tài)量估計(jì)收斂于：X^(980)=a^1(980)a^2(980)=-0.669-0.338---(11)]]>對應(yīng)的AR(2)模型為：z(k)＝-0.669z(k-1)-0.338z(k-2)+ε(k)(12)如果仍采用傳統(tǒng)Y-W法，但只對前980個(gè)陀螺隨機(jī)噪聲數(shù)據(jù)進(jìn)行建模，其建模結(jié)果為：z(k)＝-0.682z(k-1)-0.349z(k-2)+ε(k)(13)分別將發(fā)明的基于自適應(yīng)Kalman濾波的陀螺隨機(jī)噪聲建模結(jié)果(980點(diǎn)數(shù)據(jù))，與傳統(tǒng)Y-W法建模的結(jié)果(360000點(diǎn)數(shù)據(jù)和980點(diǎn)數(shù)據(jù))列于下表：表2發(fā)明方法與傳統(tǒng)Y-W法的建模結(jié)果對比從表2可以看出，只用980點(diǎn)陀螺隨機(jī)噪聲數(shù)據(jù)進(jìn)行建模時(shí)，發(fā)明的基于自適應(yīng)Kalman濾波的陀螺隨機(jī)噪聲AR建模方法的建模結(jié)果(式(12))與傳統(tǒng)的Y-W法建模結(jié)果(式(13))相比，更接近1小時(shí)完整數(shù)據(jù)的建模結(jié)果(式(10))。說明發(fā)明方法在快速建模方面具有優(yōu)勢。此外，當(dāng)有新的陀螺隨機(jī)噪聲樣本數(shù)據(jù)到來時(shí)，傳統(tǒng)的AR建模方法必須對全部歷史數(shù)據(jù)重新計(jì)算一次才能得到新的模型參數(shù)。而自適應(yīng)Kalman濾波器則可以利用新的噪聲觀測數(shù)據(jù)對AR模型的參數(shù)估值進(jìn)行修正。既解決了重復(fù)計(jì)算問題，也可以使建立的模型能夠?qū)崟r(shí)跟蹤陀螺隨機(jī)噪聲特性的變化，建模的精度高。全部歷史數(shù)據(jù)重新計(jì)算一次才能得到新的模型參數(shù)。而自適應(yīng)Kalman濾波器則可以利用新的噪聲觀測數(shù)據(jù)對AR模型的參數(shù)估值進(jìn)行修正。既解決了重復(fù)計(jì)算問題，也可以使建立的模型能夠?qū)崟r(shí)跟蹤陀螺隨機(jī)噪聲特性的變化，建模的精度高。發(fā)明方法在實(shí)際應(yīng)用時(shí)，選用的AR模型的階數(shù)并不唯一，本發(fā)明方法同樣適用于對陀螺隨機(jī)噪聲其它階AR模型的建模，例如AR(1)，AR(3)等，前提是該陀螺隨機(jī)噪聲數(shù)據(jù)的自相關(guān)系數(shù)函數(shù)和偏相關(guān)系數(shù)函數(shù)的特性滿足該模型的要求；但不同階數(shù)的AR模型在進(jìn)行自適應(yīng)Kalman濾波建模時(shí)，濾波參數(shù)的維數(shù)需要隨之改變；例如如果陀螺隨機(jī)噪聲數(shù)據(jù)適合AR(3)模型建模，則系統(tǒng)狀態(tài)量將變?yōu)椋篨＝[a1(k)，a2(k)，a3(k)]T；系統(tǒng)觀測陣將變?yōu)椋篐(k)＝[z(k-1)，z(k-2)，z(k-3)]；一步預(yù)測誤差方差陣Pk|k-1、估計(jì)誤差方差陣Pk|k和濾波增益矩陣Kk也將由原來的2階變?yōu)?階。上述實(shí)施例不以任何形式限制本發(fā)明，凡采用等同替換或等效變換的方式所獲得的技術(shù)方案，均落在本發(fā)明的保護(hù)范圍內(nèi)。當(dāng)前第1頁1 2 3