利用時頻分析法精化GNSS/INS組合導航先驗濾波模型*

甘 雨 隋立芬 張清華 桑 渤

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

2)72946 部隊，淄博255000

1 引言

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

然而，應用Kalman 濾波時，先驗信息的確定往往具有一定的隨意性和較強的經驗性，缺乏理論基礎。如果先驗模型與真實模型相差較大，不僅影響經典Kalman 濾波，也會降低自適應濾波的可靠性。在組合導航系統中，狀態隨機模型的確定是難點。現有文獻極少涉及先驗模型確定問題，影響了Kalman 濾波在組合導航中的實際應用。

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

2 GNSS/INS 組合導航的隨機模型

假設線性系統用微分方程可表示為

式中，X(t)為系統狀態，u(t)為系統驅動白噪聲過程，F(t)和G(t)為系數陣。u(t)滿足:

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

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

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

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

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

設εc和▽c的驅動白噪聲分別為ηε和η▽，則系統的狀態噪聲u 為:

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

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

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

3.1 角速率白噪聲的模型

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

1)功率譜密度法

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

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

2)Allan 方差法

角度隨機游走系數與功率譜密度的轉換關系為［10，11］:

3)方差法

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

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

如果Δf 未知，可取其上限Fs，Fs 為采樣頻率。則有qwx的計算公式為

式中，ΔT 為采樣間隔。

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

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

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

其中，Tc為相關時間，ηεx即為該G-M 誤差成分的驅動白噪聲，作為組合系統狀態噪聲的一部分，其譜密度為

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

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

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

圖2 一階G-M 過程的自相關函數Fig.2 Autocorrelation function of 1st order G-M process

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

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

3.3 角速率隨機游走的模型

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

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

3.3.1 Allan 方差法

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

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

顯然有:

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

3.3.2 方差法

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

變形可得:

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

3.4 應用中的模型與方法選擇

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

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

在模型精化的各種方法中，Allan 方差法和自相關函數法都需要對長時間的靜態慣性數據進行分析才能達到滿意的精度，方差法在噪聲成分比較單一的條件下才滿足其理論條件，否則計算的方差不能反映該類噪聲的特點。可見，方法的選擇與數據的長度與特性有關，應根據實際的分析結果進行判斷取舍。

4 計算與分析

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

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

濾波初始參數及觀測值精度設置如下:初始姿態失準角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。

對于狀態隨機模型相關的濾波參數q，用時頻特性分析方法進行計算。

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

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

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

表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

由上述結果可以看出:

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

2)使用時頻分析方法對實測數據分析計算得到的濾波參數能夠合理地平衡INS 和GPS 這兩大系統的關系，可以在不需要過多經驗信息和復雜調試的條件下精化濾波先驗模型。當先驗模型較為準確時，GPS/INS 對GPS 具有類似于擬合的效果，可以平滑掉部分GPS 誤差。

3)組合導航系統在保證精度的前提下具有更高的結果輸出頻率，在本例中，輸出頻率為100 Hz，結果的精度優于單獨的GPS 系統或INS 系統。

4)在少數歷元內，組合導航系統的誤差大于GPS 的誤差，這是因為運動中存在一些狀態突變，導致了與先驗濾波模型之間的偏差，需要用自適應濾波理論加以處理。

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

5 結語

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

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

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

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 夏啟軍，孫優賢，周春暉.漸消卡爾曼濾波器的最佳自適應算法及其應用［J］.自動化學報，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 歐吉坤，柴艷菊，袁運斌.自適應選權濾波［A］.大地測量與地球動力學進展［C］.武漢:湖北科學技術出版社，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 楊元喜，何海波，徐天河.論動態自適應濾波［J］.測繪學報，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 秦永元，張洪鉞，汪叔華.卡爾曼濾波與組合導航原理［M］.西 安:西 北 工 業 大 學 出 版 社，1998.(Qin Yongyuan，Zhang Hongyue and Wang Shuhua.Kalman filter and the principle of integrated navigation［M］.Xian:Northwest Industrial Press，1998)

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

10 毛奔，林玉榮.慣性器件測試與建模［M］.哈爾濱:哈爾濱工程大學出版社，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］.空間控制技術與應用，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.