基于誤差時間特性的INS/GNSS觀測噪聲方差在線計算方法

邵夢晗,高曉穎,呂建強

(1.北京航天自動控制研究所， 北京 100854； 2.宇航智能控制技術國家級重點實驗室， 北京 100854)

精確制導與精確入軌效能主要受制于導航系統性能，導航系統的精度從根本上決定了后續制導、控制的精度。目前工程上廣泛使用基于卡爾曼濾波的INS/GNSS組合導航，其濾波效果依賴于噪聲統計特性等先驗信息。但對飛行環境復雜、高機動的載體，多路徑效應、可見衛星數目及幾何分布、外界電磁波干擾等多種因素并存，造成GNSS觀測噪聲多變且難以預測，為此涌現了多種自適應濾波算法[1-4]。但大部分自適應算法都是從卡爾曼濾波算法本身出發，獲得觀測量與濾波參數的自適應調節關系，同時又引入了其他需要事先確定的參數，導致其對不同噪聲環境的適應性、可靠性不高，難以在實際工程中應用。文獻[5]提出一種組合導航濾波參數在線估計的并行架構，并將遺傳算法引入進行參數優化，但是適值函數難以確定。文獻[6]提出了一種根據傳感器不同測量特性的自適應濾波算法，簡單易行，但其只進行了理論分析，且僅限于組合導航開環校正使用。然而目前工程上廣泛使用速度、位置反饋校正，甚至全狀態量閉環校正方式。為此，本文針對組合導航閉環校正方式，依據常用慣性導航短時高精度和GNSS測量噪聲的時間不相關特性，提出一種基于小樣本統計的在線計算GNSS觀測噪聲方差陣的方法；并利用實測數據驗證了這種誤差特性，為該思路的可行性提供了數據支撐；最后給出了具體的實現步驟，并利用一組飛行仿真數據對該方法進行了驗證。

1 閉環校正下GNSS觀測噪聲方差在線計算方法

1.1 發射慣性系下組合濾波模型

在發射慣性系下，以速度、位置、姿態、加速度計零偏、陀螺零偏的誤差量作為組合濾波狀態量X(t)，以INS與GNSS的位置差、速度差作為觀測量Z(t)，構建如下濾波模型：

(1)

Z(t)=H(t)X(t)+Wgnss(t)

(2)

X(t)=(δvx,δvy,δvz,δrx,δry,δrz,

δφx,δφy,δφz,δK0x,δK0y,δK0z,δD0x,δD0y,δD0z)T

Z(t)=(δvx,δvy,δvz,δrx,δry,δrz)T

其中，δvx,δvy,δvz為速度誤差；δrx,δry,δrz為位置誤差；δφx,δφy,δφz為姿態角誤差；δK0x,δK0y,δK0z為加速度零次項誤差；δD0x,δD0y,δD0z為陀螺零次項誤差；Wins(t)、Wgnss(t)分別是系統白噪聲陣、觀測白噪聲陣；F(t)是系統矩陣；G(t)是系統噪聲驅動陣；H(t)是觀測矩陣。

1.2 慣性導航短期高精度信息利用

衛星導航速度誤差在0.1 m/s量級、位置誤差在10 m量級，且測量誤差與時間不相關，誤差不隨著時間積累；常用激光慣組陀螺零偏在0.01°/h量級、加速度計零偏在10-5g量級，慣性導航更新頻率快，短期精度高，但誤差隨著時間累積而發散。但是慣性導航在1 s、5 s乃至更長一段時間內的測量誤差都是小于衛星導航測量誤差的，為此，考慮充分挖掘慣性導航的短期高精度信息，對GNSS觀測噪聲進行在線估計。但載體實際飛行時無真實值參考，因此，聯想到利用INS、GNSS各自相鄰時刻增量信息：慣組測量噪聲遠小于GNSS接收機測量噪聲，且相鄰時刻作差能夠抵消一部分慣組不隨時間變化的誤差，導致純慣性導航短期增量精度高于衛星導航的增量精度，以此來估計GNSS觀測噪聲方差陣。

1.3 GNSS觀測噪聲方差在線計算方法

以誤差量為卡爾曼濾波器狀態量，得到的估計有兩種利用方法[7-8]：一種是將估計作為慣導系統輸出的校正量，稱為開環法；另一種是將估計反饋到慣導系統中，作為慣性導航力學編排方程中的相應參數或校正量，即閉環法。開環校正的濾波方程就是離散化后的卡爾曼濾波方程；閉環校正每次校正后需要把狀態量清零，其余與開環校正無異。

開環校正下，INS、GNSS的k時刻速度誤差、位置誤差狀態估計量：

(3)

其中，Xv,r(k)是速度誤差、位置誤差狀態量真值；fins(k)是慣性導航誤差趨勢項。

閉環校正下，組合濾波反饋修正公式：

(4)

由慣性導航更新公式：

(5)

(6)

其中，dδXv,r(k-1)是速度、位置修正量誤差。取GNSS更新周期為組合濾波周期，分別求INS、GNSS相鄰濾波時刻的狀態增量：

ΔINS(k)=Xins(k)-Xins(k-1)=[X(k)-X(k-1)]+

[fins(k-1)-fins(k)]+[Wins(k-1)-Wins(k)]+

[dδXv,r(k-2)-dδXv,r(k-1)]

(7)

ΔGNSS(k)=Xgnss(k)-Xgnss(k-1)=

[X(k)-X(k-1)]+[Wgnss(k-1)-Wgnss(k)]

(8)

式(7)、式(8)作差：

ΔINS(k)-ΔGNSS(k)=[fins(k-1)-fins(k)]+

[Wins(k-1)-Wins(k)]-[Wgnss(k-1)-Wgnss(k)]+

[dδXv,r(k-2)-dδXv,r(k-1)]

(9)

求式(9)的自協方差，并考慮Wins(k)、Wgnss(k)是零均值白噪聲，有：

[fINS(k-1)-fINS(k)]2+Q(k-1)+Q(k)+R(k-1)+

R(k)+Pv,r(k-2)+Pv,r(k-1)

(10)

其中，Q(k)是系統噪聲方差陣，R(k)是觀測噪聲方差陣，Pv,r(k-2)、Pv,r(k-1)分別是k-2、k-1時刻組合濾波速度、位置狀態量誤差方差陣，是6×6維矩陣。

考慮慣組測量噪聲遠小于GNSS接收機測量噪聲Q(k)<

[fINS(k-1)-fINS(k)]2<

(11)

故上式(10)可寫為

R(k-1)+R(k)+Pv,r(k-2)+Pv,r(k-1)

(12)

從而，可得到GNSS觀測噪聲方差陣R(k)的統計估計：

[Pv,r(k-2)+Pv,r(k-1)]}

(13)

2 實測數據驗證與定量分析

截取某飛行器某段時間的純慣性導航、衛星導航的實測數據，激光慣組和GNSS接收機均獨立、正常工作。慣性導航更新周期是0.02 s，衛星導航更新周期是1 s。為同步求增量，慣性導航需要進行1 s的遞推解算，然后分別求INS、GNSS相鄰1 s的速度增量Δv、位置增量Δr，繪制出慣性導航與衛星導航的增量曲線如圖1、圖2所示。

由圖1、圖2可以看出：GNSS結果增量在慣性導航結果增量附近波動。兩者的增量均值誤差可忽略不計，任取某時間段的樣本數據，計算純慣性導航速度增量標準差為0.009 6 m/s，位置增量標準差為0.566 1 m；GNSS速度增量標準差為0.198 4 m/s，位置增量標準差為6.236 7 m。可見，純慣性導航增量的標準差遠小于對應的GNSS結果增量的標準差，即：

[fINS(k-1)-fINS(k)]2+Q(k)+Q(k-1)<

(14)

進而說明式(12)取近似是合理的。

圖1 速度增量曲線

圖2 位置增量曲線

進一步，根據實際飛行情況進行定量分析，假設載體運動速度v=3 000 m/s，加速度a=20 m/s2，速度方向誤差增量Δθ=1′，則Δt=1 s時的速度誤差增量、位置誤差增量分別為

Δvins=a×Δθ×Δt=20×1/60/180×π=0.005 8 m/s

Δrins=v×Δθ×Δt=3 000×1×1/60/180×π=

0.872 639 m

可以看出：即使在較高速度和大誤差角增量的假設下，1 s內純慣性導航的速度誤差增量只有0.006 m/s、位置誤差增量只有0.9 m。而GNSS測速誤差一般為0.1 m/s量級，定位誤差一般為10 m量級。因此，1 s內慣性導航誤差增量的平方值遠小于GNSS定位誤差的方差[6]，所以式(11)的假設是合理的，因此式(12)是符合客觀實際的。

3 算法實現

考慮實際算法的穩定性和可行性，具體實現步驟如下：

1) 計算INS、GNSS各自相鄰時刻速度、位置狀態增量ΔINS、ΔGNSS，并計算：

B(k)=ΔINS(k)-ΔGNSS(k)

(15)

B(k)序列含有GNSS測量噪聲信息，因此可以通過對一段時間的樣本數據進行統計，得到觀測噪聲陣的估計。樣本容量越大，統計結果的可靠性越高，但當實際噪聲發生變化時，樣本容量越大，這種動態信息被平滑，即當前時刻的動態信息淹沒在大量的歷史信息中，導致統計結果不能及時跟蹤動態的噪聲信息，實時性較差。樣本容量越小，統計結果的可靠性降低。所以，通過滑動窗口利用最新一段導航信息進行方差統計，具體公式為

[Pv,r(k-2)+Pv,r(k-1)]}

(16)

其中，k為當前濾波時刻，k>M；M為滑動窗口長度，一般為20～50 s，可根據實際飛行情況選取。

2) 粗差剔除。B(k)受傳感器觀測粗差的影響較大，因此在計算前需要剔除粗差。一般認為，高精度激光慣組出現粗差的可能性很小，因此，主要是考慮GNSS出現粗差的情況。本文采用標準化新息作為假設檢驗統計量進行粗差探測，具體方法可參考文獻[9-11]。若探測到某時刻存在觀測粗差，將該時刻數據剔除，不參與統計估計。

3) 觀測噪聲方差陣R(k)的估計。為避免出現小樣本統計下R(k)變化較為劇烈的情況，引入指數加權法[8]更新R(k)陣：

(17)

其中，b為遺忘因子，一般取0.95～0.995。

4 仿真驗證

4.1 仿真設置

采用一組模擬飛行數據進行仿真驗證，在發射慣性系下，初始速度為[406,0.1,-53]m/s；初始位置為[0.1,0.1,0.1]m；初始姿態為[90,-0.1,0.1]°；陀螺零偏0.1°/h，測量噪聲0.05°/h；加速度計零偏0.002g，測量噪聲50 μg；組合濾波時間為1 302 s。仿真設置R(k)陣為對角陣，其值是GNSS觀測白噪聲的方差。為考察閉環校正下在線計算R(k)陣的效果，在280～580 s設置觀測噪聲方差發生突變，具體R(k)陣對角線元素仿真真值如表1所示。

表1 R(k)仿真真值

為進一步驗證組合導航閉環校正下在線計算方法的有效性，將其與同條件下標準卡爾曼濾波、理想方法對比分析。其中，標準卡爾曼濾波方法與在線計算方法的觀測噪聲方差陣初值相同，為充分驗證該方法的有效性，設置觀測噪聲方差陣初值為：Rv0=0.04(m/s)2，Rr0=400(m)2，與仿真真值相差較大；理想方法下觀測噪聲方差陣始終等于仿真真值，其結果可以作為參考值對比。

4.2 仿真結果與分析

將在線計算方法、標準卡爾曼濾波方法、理想方法解算的速度、位置分別與仿真軌跡真值作差，得到速度誤差、位置誤差如圖3、圖4所示。

圖3 速度誤差

圖4 位置誤差

同時統計了觀測噪聲發生突變時間段的速度合成誤差、位置合成誤差均值和均方根，如表2所示。

分析圖3、圖4，在初始觀測噪聲方差陣誤差較大情況下，本文提出的在線計算方法能夠較好地根據實際情況進行調整，并能及時跟蹤噪聲突變，精度與理想方法較為接近；而標準卡爾曼濾波算法始終采用較差的初始觀測噪聲方差陣計算，速度、位置估計出現了較大的誤差。分析表2，噪聲突變段，在線計算方法解算的速度誤差均值0.016 m/s，位置誤差均值0.567 m，速度誤差均方根0.097 m/s，位置誤差均方根1.255 m，精度與理想方法相當；而標準卡爾曼濾波得到的速度誤差均值0.1 m/s，位置誤差均值13.864 m，速度誤差均方根0.210 m/s，位置誤差均方根9.416 m。可見，本文提出的在線計算濾波方法性能明顯優于標準卡爾曼濾波算法。

表2 噪聲突變時間段(280～580 s)導航誤差均值、均方根

注：表2中方法1是本文提出的在線計算觀測噪聲方差陣方法；方法2是標準卡爾曼濾波方法；方法3是理想方法。

繪制3種方法下的每濾波時刻R(k)陣對角線元素值如圖5所示。

圖5 速度、位置觀測噪聲方差

可以看出：標準卡爾曼濾波方法測量噪聲方差始終為初始值，誤差較大；而在線計算R(k)方法是以過去一段時間的滑動窗口數據為樣本進行統計計算，同時為避免大幅震蕩采用了指數加權，因此存在一定的振蕩和滯后，但是能夠隨著真實觀測噪聲的變化而變化，始終跟蹤實際值，效果較好。

5 結論

本文提出了一種基于INS、GNSS不同誤差時間特性的在線計算GNSS觀測噪聲方差陣的方法。該方法利用滑動時間窗口內的增量樣本數據在線計算觀測噪聲方差陣，避免了傳統自適應濾波算法引入其他參數調整的問題，對工程上廣泛使用的組合導航閉環校正適用。仿真試驗表明該方法能夠準確跟蹤觀測噪聲的變化，獲得較高精度的導航狀態量，對環境適應性較強，可應用于飛行環境復雜、高機動載體的導航測量中。