甘 雨 隋立芬 張清華 桑 渤

1)信息工程大學(xué)地理空間信息學(xué)院，鄭州 450052

2)72946 部隊，淄博255000

1 引言

GNSS/INS 組合導(dǎo)航應(yīng)用Kalman 濾波進行數(shù)據(jù)融合，經(jīng)典Kalman 濾波要求有可靠的函數(shù)模型及隨機模型。先驗濾波模型的確定不僅困難，而且有時并不能完全合理地描述整個導(dǎo)航過程，因此許多學(xué)者研究了多種自適應(yīng)濾波理論［5－6］，并對經(jīng)典Kalman 濾波進行了擴展，解決了許多實際應(yīng)用問題，特別是抗差自適應(yīng)濾波，可以同時控制狀態(tài)擾動和觀測異常的影響。

然而，應(yīng)用Kalman 濾波時，先驗信息的確定往往具有一定的隨意性和較強的經(jīng)驗性，缺乏理論基礎(chǔ)。如果先驗?zāi)Ｐ团c真實模型相差較大，不僅影響經(jīng)典Kalman 濾波，也會降低自適應(yīng)濾波的可靠性。在組合導(dǎo)航系統(tǒng)中，狀態(tài)隨機模型的確定是難點?，F(xiàn)有文獻極少涉及先驗?zāi)Ｐ痛_定問題，影響了Kalman 濾波在組合導(dǎo)航中的實際應(yīng)用。

組合導(dǎo)航系統(tǒng)本質(zhì)上是連續(xù)線性系統(tǒng)，具有功率譜密度形式的頻域隨機模型。提出利用功率譜密度分析、自相關(guān)函數(shù)分析、Allan 方差分析等時頻特性分析方法，由慣性元件數(shù)據(jù)計算各狀態(tài)噪聲成分的功率譜密度，為組合導(dǎo)航提供較為準確的先驗信息，避免實際應(yīng)用Kalman 濾波時繁瑣的調(diào)試。算例驗證了時頻分析精化先驗濾波模型的有效性。

2 GNSS/INS 組合導(dǎo)航的隨機模型

假設(shè)線性系統(tǒng)用微分方程可表示為

式中，X(t)為系統(tǒng)狀態(tài)，u(t)為系統(tǒng)驅(qū)動白噪聲過程，F(xiàn)(t)和G(t)為系數(shù)陣。u(t)滿足:

q 為u(t)的功率譜密度陣，又稱方差強度陣。問題的關(guān)鍵是求取狀態(tài)噪聲譜密度陣q。

陀螺誤差可表示為［7，8］:

其中，εc為非白噪聲部分，wε為白噪聲部分。εc可以描述為隨機常值、高斯馬爾可夫過程、角速率隨機游走或者它們的組合形式。

同理，加速度誤差▽可表示為:

其中，▽c為非白噪聲部分，w▽為白噪聲部分。

設(shè)εc和▽c的驅(qū)動白噪聲分別為ηε和η▽，則系統(tǒng)的狀態(tài)噪聲u 為:

由此可知，精化組合導(dǎo)航的狀態(tài)隨機模型，就是由陀螺和加速度計數(shù)據(jù)計算wε、w▽、ηε、η▽的功率譜密度，然后形成譜密度陣q。

3 基于時頻特性分析的先驗濾波模型精化

以X 軸陀螺為例說明利用時頻特性分析提取噪聲成分功率譜密度信息的方法。

3.1 角速率白噪聲的模型

設(shè)角速率白噪聲對應(yīng)于陀螺誤差中的wεx(wε的X 軸分量)。對于該誤差項，只需研究其隨機模型，即確定wεx的功率譜密度qωx共有三種方法。

1)功率譜密度法

白噪聲的譜密度為常數(shù)，對應(yīng)于雙對數(shù)功率譜密度(logSω(f)－logf)斜率為0 的部分。先計算靜態(tài)陀螺數(shù)據(jù)的功率譜密度Sω(f)，再對雙對數(shù)譜密度(logSω(f)－logf)擬合出斜率為0 的直線，與縱軸交點為logqωx，qωx即為角速率白噪聲的譜密度。若計算的是單邊譜密度［9］而非雙邊譜密度Sω(f)，雙對數(shù)單邊譜擬合直線與縱軸交點為log2qωx。

圖1 角速率白噪聲譜密度擬合Fig.1 Fitting the PSD of angular rate white noise

2)Allan 方差法

角度隨機游走系數(shù)與功率譜密度的轉(zhuǎn)換關(guān)系為［10，11］:

3)方差法

功率譜密度對頻率的積分即為方差，對于某一純粹的白噪聲過程x(t)，Sx(f)為常數(shù)，則有:

其中，Δf 為頻率帶寬。功率譜密度為

如果Δf 未知，可取其上限Fs，F(xiàn)s 為采樣頻率。則有qwx的計算公式為

式中，ΔT 為采樣間隔。

式(9)是純粹白噪聲過程的功率譜密度求取方法。對于許多陀螺數(shù)據(jù)，我們無法求取其中白噪聲的方差，只能求取陀螺數(shù)據(jù)整體的方差，而這一方差還包括其他有色噪聲成分的影響，因此，該方法在白噪聲起主導(dǎo)作用的陀螺上比較有效。

3.2 高斯馬爾可夫過程的模型

常用一階高斯馬爾可夫(G-M)過程對陀螺的相關(guān)噪聲成分進行建模，設(shè)一階G-M 模型表示的X 軸陀螺誤差成分為εcgx，則有:

其中，Tc為相關(guān)時間，ηεx即為該G-M 誤差成分的驅(qū)動白噪聲，作為組合系統(tǒng)狀態(tài)噪聲的一部分，其譜密度為

其中，σ2為εcgx的方差。

將εcgx引入到組合導(dǎo)航的狀態(tài)方程中，需要求得其方差σ2及相關(guān)時間Tc。雖然通過Allan 方差可以得到一些G-M 過程的信息，但是Allan 方差對G-M 過程的表達不夠精確［12］，更理想的分析方法是自相關(guān)函數(shù)法。

利用自相關(guān)函數(shù)經(jīng)過式(12)可以求取σ2及Tc，即可求得驅(qū)動白噪聲ηεx的譜密度qηx。相關(guān)時間Tc=300 s，方差σ2=1.0 的一階G-M 過程的自相關(guān)函數(shù)如圖2 所示。

圖2 一階G-M 過程的自相關(guān)函數(shù)Fig.2 Autocorrelation function of 1st order G-M process

對于高階G-M 過程，其模型比一階G-M 過程略微復(fù)雜，建模后狀態(tài)參數(shù)的個數(shù)更多，基本分析方法與一階G-M 過程類似。

許多文獻使用AR 模型描述慣性元件誤差，實際上，任意階G-M 過程都可以用一定階數(shù)的AR 過程表示。這里側(cè)重于用連續(xù)的隨機過程表達系統(tǒng)的物理特征。

3.3 角速率隨機游走的模型

隨機游走可以看成是相關(guān)時間很長的一階G-M過程。當Tc→∞時，一階G-M 過程變成，即隨機游走。設(shè)角速率隨機游走表示的X 軸陀螺誤差成分為εcrx，則:

式中，ξηx為該過程的驅(qū)動白噪聲，譜密度為qξx。易知，隨機游走的方差隨時間而增長，而其本質(zhì)上為非平穩(wěn)過程，因此無法對其進行自相關(guān)分析或功率譜分析，需要使用其它方法求取ξηx的功率譜密度qξx。

3.3.1 Allan 方差法

角速率隨機游走的功率譜密度為［10］:

K 為速率隨機游走系數(shù)。根據(jù)線性系統(tǒng)傳遞函數(shù)的理論可由驅(qū)動白噪聲ξηx的譜密度qξx求得隨機游走εcrx功率譜密度為［13］:

顯然有:

對Allan 方差進行最小二乘擬合或分段最小二乘擬合，可以求取K 的值，如果譜密度的單位取deg2/h4/Hz，可將式(16)用下式代替:

3.3.2 方差法

Maybeck 給出了隨機游走過程的驅(qū)動白噪聲譜密度所滿足的關(guān)系式［14］，對ξηx而言即:

變形可得:

εcrx(t+ΔT)－εcrx(t)/ΔT 實際上近似為隨機游走過程的微分——白噪聲ξηx。顯然，式(19)本質(zhì)上與式(9)相同。同理，該方法在角速率隨機游走為主要噪聲源的陀螺上比較有效。

3.4 應(yīng)用中的模型與方法選擇

若用εcbx表示X 軸陀螺隨機常值分量，則有:

在實際應(yīng)用中，如果對各軸陀螺和加速度計建模，并不一定將上述四種模型都引入到濾波的狀態(tài)方程中。仍以X 軸陀螺為例，如果同時引入wεx、εcgx、εcrx、εcbx，除白噪聲wεx外，其他三種屬于非白噪聲εc的模型都將增加狀態(tài)向量的維數(shù)，加重濾波解算的負擔。一般情況下，白噪聲wεx要引入到狀態(tài)模型，它不會增加濾波負擔。εcbx反映一些非隨機信息，一般要引入，但在需要保證濾波效率的情況下也可不引入。εcgx和εcrx的作用都是描述噪聲的時相關(guān)性，實用中只選其一。有文獻認為隨機游走過程并不是描述角速率誤差的現(xiàn)實模型，因為隨機游走的方差無限增大，而角速率輸出中一般會含有指數(shù)相關(guān)的一階馬爾可夫過程和寬帶白噪聲［15］。

在模型精化的各種方法中，Allan 方差法和自相關(guān)函數(shù)法都需要對長時間的靜態(tài)慣性數(shù)據(jù)進行分析才能達到滿意的精度，方差法在噪聲成分比較單一的條件下才滿足其理論條件，否則計算的方差不能反映該類噪聲的特點?？梢?，方法的選擇與數(shù)據(jù)的長度與特性有關(guān)，應(yīng)根據(jù)實際的分析結(jié)果進行判斷取舍。

4 計算與分析

采用的GPS/INS 實驗數(shù)據(jù)共約30 分鐘，其中初始靜止狀態(tài)約5 分鐘。IMU 采樣頻率100 Hz，GPS 采樣頻率1 Hz。IMU 的XYZ 三軸分別指向下、左、前方向，標稱陀螺漂移1 deg/h，標稱加速度偏置10－4g。對GPS/INS 數(shù)據(jù)進行松組合導(dǎo)航。采用雙差C/A 碼解算的位置和Doppller 解算的速度作為GPS 輸出值與INS 輸出值構(gòu)成觀測量，以載波相位差分獲得的位置作為參考解。

組合導(dǎo)航濾波狀態(tài)變量中，加速度計誤差中可忽略時相關(guān)誤差，這是由于該分量相對較小，同時也為了使濾波器的維數(shù)盡量低些［8］，因而這里加速度計誤差用隨機常值和白噪聲進行建模。陀螺誤差也包括了隨機常值和白噪聲，但還使用高斯馬爾可夫模型描述時相關(guān)誤差。整個系統(tǒng)的狀態(tài)噪聲為u=［wεηεw▽］T，其譜密度矩陣為q，有

濾波初始參數(shù)及觀測值精度設(shè)置如下:初始姿態(tài)失準角100.0 s、100.0 s、500.0 s，陀螺隨機常值和高斯馬爾可夫過程初始標準差1.0 deg/h，加速度計隨機常值的初始標準差10－4g;GPS 水平方向位置標準差1.0 m，速度標準差0.05 m/s，垂直方向位置標準差2.0 m，速度標準差0.10 m/s。

對于狀態(tài)隨機模型相關(guān)的濾波參數(shù)q，用時頻特性分析方法進行計算。

由于條件限制，只有5 分鐘的靜止數(shù)據(jù)。因為數(shù)據(jù)太短，Allan 方差和自相關(guān)函數(shù)的估計精度不高，因而這里主要用功率譜密度法提取白噪聲譜密度qωx和q▽x，對于高斯馬爾可夫噪聲，這里的相關(guān)時間的估計精度也受數(shù)據(jù)長度限制，只能根據(jù)元件規(guī)格取經(jīng)驗值Tc=600 s，再由噪聲方差通過式(12)計算高斯馬爾可夫驅(qū)動白噪聲的譜密度qηx。分析計算的結(jié)果如表1 所示，將其代入式(24)即可得到頻域譜密度矩陣q。

表1 功率譜密度計算結(jié)果Tab.1 Calculating results of PSD

對GPS/INS 組合數(shù)據(jù)分別進行INS 解算、雙差C/A 碼GPS 解算和GPS/INS 組合解算。圖3 所示為INS 緯度及參考緯度，虛線為參考值，點線為INS值;圖4 所示為GPS/INS 的緯度誤差和GPS 的緯度誤差，“+”為GPS 結(jié)果，粗實線為GPS/INS 結(jié)果，RMS 及最大誤差(MAX)如表2 所示;圖5 所示為283560.0－283566.0 s 內(nèi)的GPS 及GPS/INS 緯度值，“+”為GPS 結(jié)果，實線為GPS/INS 結(jié)果。

表2 GPS 和GPS/INS 的RMS 和MAXTab.2 RMS and MAX of GPS and GPS/INS

圖3 INS 緯度與參考緯度Fig.3 INS latitude and the reference latitude

圖4 GPS 與GPS/INS 的緯度誤差Fig.4 Latitude error of GPS and GPS/INS

由上述結(jié)果可以看出:

1)受到慣性元件誤差、初始誤差等因素的影響，單獨INS 的誤差隨時間積累，需要借助其他導(dǎo)航手段如GNSS 對其誤差進行抑制。但是當GNSS 失鎖時，只能依賴于INS，此時可通過對INS 誤差進行消噪處理來削弱誤差。

2)使用時頻分析方法對實測數(shù)據(jù)分析計算得到的濾波參數(shù)能夠合理地平衡INS 和GPS 這兩大系統(tǒng)的關(guān)系，可以在不需要過多經(jīng)驗信息和復(fù)雜調(diào)試的條件下精化濾波先驗?zāi)Ｐ?。當先驗?zāi)Ｐ洼^為準確時，GPS/INS 對GPS 具有類似于擬合的效果，可以平滑掉部分GPS 誤差。

3)組合導(dǎo)航系統(tǒng)在保證精度的前提下具有更高的結(jié)果輸出頻率，在本例中，輸出頻率為100 Hz，結(jié)果的精度優(yōu)于單獨的GPS 系統(tǒng)或INS 系統(tǒng)。

4)在少數(shù)歷元內(nèi)，組合導(dǎo)航系統(tǒng)的誤差大于GPS 的誤差，這是因為運動中存在一些狀態(tài)突變，導(dǎo)致了與先驗濾波模型之間的偏差，需要用自適應(yīng)濾波理論加以處理。

圖5 GPS 與GPS/INS 緯度結(jié)果對比Fig.5 Comparison between the latitude results of GPS and GPS/INS

5 結(jié)語

GNSS/INS 組合導(dǎo)航的狀態(tài)噪聲由陀螺儀和加速度計的誤差所構(gòu)成，建立先驗濾波模型的前提是對陀螺和加速度計中的誤差建模，尤其是以功率譜密度為代表的頻域隨機模型的計算。

由于連續(xù)系統(tǒng)的頻域隨機模型和離散系統(tǒng)的時域隨機模型可以互相轉(zhuǎn)換，也可以由數(shù)據(jù)分析計算出幾種常見誤差源的時域隨機模型，但是其過程更為復(fù)雜，而且頻域的描述更貼近于系統(tǒng)的本質(zhì)，因此這里主要使用頻域模型。

利用時頻分析精化先驗濾波模型的方法，既減少了確定濾波模型時的主觀性和經(jīng)驗性，也避免了實際濾波應(yīng)用前的繁瑣調(diào)試。它基于實際的數(shù)據(jù)進行分析，該方法不僅可以用于組合導(dǎo)航，對其他的Kalman 濾波相關(guān)應(yīng)用領(lǐng)域都具有一定參考價值。

1 Jazwinski A H.Stochastic processes and filtering theory［M］.New York:Mathematics in Science and Engineering，Academic Press，1970.

2 Mohamed A H and Schwarz K P.Adaptive Kalman Filtering for INS/GPS［J］.Journal of Geodesy，1999，73(4):193－203.

3 夏啟軍，孫優(yōu)賢，周春暉.漸消卡爾曼濾波器的最佳自適應(yīng)算法及其應(yīng)用［J］.自動化學(xué)報，1990，16(3)，210－216.(Xia Qijun，Sun Youxian and Zhou Chunhui.An optimal adaptive algorithm for fading Kalman filter and its application［J］.Acta Automatica Sinaca，1990，16(3)，210－216)

4 歐吉坤，柴艷菊，袁運斌.自適應(yīng)選權(quán)濾波［A］.大地測量與地球動力學(xué)進展［C］.武漢:湖北科學(xué)技術(shù)出版社，2004，816－823.(Ou Jikun，Chai Yanju and Yuan Yunbin.Adaptive flitering by selecting the parameter weight factor［A］.Progress in geodesy and geodynamics［C］.Wuhan:Hubei Science and Technology Press，2004，816－823 )

5 Yang Yuanxi，He Haibo and Xu Guochang.Adaptively robust filtering for kinematic geodetic positioning［J］.Journal of Geodesy，2001，75(2/3):109－116.

6 楊元喜，何海波，徐天河.論動態(tài)自適應(yīng)濾波［J］.測繪學(xué)報，2001，30(4):293－298.(Yang Yuanxi，He Haibo and Xu Tianhe.Adaptive robust filtering for kinematic GPS positioning［J］.Acta Geodaetica et Cartographica Sinica，30(4):293－298)

7 Gelb A.Applied optimal estimation［M］.USA:The Analytic Science Corporation，1974.

8 秦永元，張洪鉞，汪叔華.卡爾曼濾波與組合導(dǎo)航原理［M］.西 安:西 北 工 業(yè) 大 學(xué) 出 版 社，1998.(Qin Yongyuan，Zhang Hongyue and Wang Shuhua.Kalman filter and the principle of integrated navigation［M］.Xian:Northwest Industrial Press，1998)

9 何正嘉，等.現(xiàn)代信號處理及工程應(yīng)用［M］.西安:西安交通大學(xué)出版社，2007.(He Zhengjia，et al，Modern signal processing and its application［M］.Xi’an:Xi’an Jiaotong University Press，2007)

10 毛奔，林玉榮.慣性器件測試與建模［M］.哈爾濱:哈爾濱工程大學(xué)出版社，2008.(Mao Ben and Lin Yurong.Inertial sensor testing and modeling［M］.Harbin:Harbin Engineering University Press，2008)

11 IEEE Std 952-1997.IEEE standard specification format guide and test procedure for single-axis interferometric fiber optic gyros［S］.New York:IEEE Standard Board，1997.

12 Flenniken W S.Modeling inertial measurement units and analyzing the effect of their errors in navigation applications［D］.Graduate Faculty of Auburn University，2005.

13 熊凱，雷擁軍，曾海波.基于Allan 方差法的光纖陀螺建模與仿真［J］.空間控制技術(shù)與應(yīng)用，2010，36(3)，7－13.(Xiong Kai，Lei Yongjun and Zeng Haibo.Modeling and simulation of fiber optic gyros based on Allan variance method［J］.Aerospace Control and Application，2010，36(3):7－13)

14 Maybeck P S.Stichastic models，estimation and control，volume 1［M］.New York:Academic Press，1979.

15 Gebre-Egziabher D.Design and performance analysis of a low-cost aided dead reckoning navigator［D］.USA:Stanford University，2004.