多系統GNSS RTK的改進卡爾曼濾波算法

范紅平，周志峰，王永泉

(1. 上海工程技術大學機械工程學院，上海 201620； 2. 上海司南衛星導航技術股份有限公司，上海 201801)

實時動態差分法(real-time kinematic,RTK)是基于載波相位觀測值來提高GNSS位置數據精度的技術，而多系統GNSS RTK融合了GPS、GLONASS、BDS和Galileo的4大衛星導航系統的定位數據，提供更多的可見衛星，實現更高的RTK定位精度[1]。Kalman濾波應用于RTK是通過動態推導，不存儲過去的觀測數據，只根據新的數據和前一時刻的估計量，借助狀態轉移方程不斷迭代，逐步實現預測位置收斂于實際位置坐標即穩態值[2]。多系統GNSS增加了可見衛星數，使狀態向量的維數急劇增加，這導致Kalman濾波方程乘法次數增加[3]。研究學者對其進行了多方面改進，趙占祥等提出一種矩陣外積法減少乘法運算次數，同時避免了矩陣求逆，占原算法CPU 1/8的時間[4]；胡朱林引用平方根協方差濾波，實現Q-R矩陣分解的自適應濾波，證明其可行性[5]；郭樹人等利用稀疏矩陣、矩陣對稱性及其求逆降維等方法，提出了快速Kalman算法，降低了1/6的乘法運算復雜度，并且CPU耗時降至原算法的1/3左右[6]。以上方法均證實有效，但對矩陣計算方法不足。

本文利用狀態變量解算每顆衛星的整周模糊度和電離層延時，利用衛星數目與矩陣維數成正比[7]，提出一個改進Kalman濾波計算方法，推導出稀疏狀態轉移矩陣F，并且對協方差矩陣分塊求解，最后利用對稱性減少n(n+1)/2次乘法(n為矩陣維數)，大幅縮短CPU運行時間。

1 GNSS RTK的Kalman濾波模型

1.1 線性高斯系統的Kalman濾波方程

建立以下Kalman濾波的系統動態方程，式(1)、式(2)為一步預測過程，得出下一時刻最優預測值[8]。式(3)、式(4)和式(5)為一步更新過程，先計算出卡爾曼增益K，然后計算狀態估計值。

(1)

(2)

(3)

(4)

Qk,k=[1-KkHk]Qk,k-1

(5)

1.2 GNSS RTK的雙差觀測方程

假設用戶接收機u和基準站接收機r同時跟蹤衛星i和衛星j，則GNSS第k歷元在L1和L2兩個載波頻率上雙差載波相位和雙差偽距的觀測方程分別為

(6)

(7)

1.3 優化Kalman濾波運算的必要性

當前理想狀態下，上海地區可見衛星數(高度角10°)為36顆(GPS、GLONASS、BDS和Galileo分別為12、8、12和4顆)，Kalman濾波方程的狀態轉移矩陣和協方差矩陣維數將達到117維[12]。僅一次迭代，式(2)的協方差矩陣的求解達到27 378次乘法，占用82 ms(32位系統一次乘法浮點運算占3 μs)。此外隨可見衛星數的變化，Xk維數也隨之變化，尤其當衛星數增加時，為保證濾波的連續性，應動態更新協方差陣Pk,k[13]。本文考慮以上因素，從兩方面改進Kalman濾波方法。

2 改進型Kalman濾波方法

2.1 改進1：狀態轉移矩陣F模型

本文取三維直角坐標構建定常加速度動態模型，假設載體處于穩定加速狀態，則狀態參數需要增加3個加速度分量，則9+3i維狀態變量為

(8)

(9)

式中，τ表示k-1到k歷元的時差；In為n維單位矩陣。

2.2 改進2：協方差矩陣Q計算

根據式(2)與式(8)，Qk,k為9+3i維矩陣[14]，可將其分為Q11、Q12、Q21和Q224塊(階數分別為9×9、9×3i、3i×9和3i×3i)，故Qk,k-1一步迭代計算

(10)

3 算例及分析

3.1 理論分析

(11)

(12)

依據普通矩陣乘法[15]，Qk,k乘法次數為H2=2×(9+3i)3次， 以H1/H2百分比為縱坐標，衛星數為橫坐標繪制柱狀圖如圖1所示。

圖1 改進算法優化比率

如圖1所示，可見衛星多于4顆，其改進算法乘法次數H1均低于普通矩陣H2的1.2%，且隨可見衛星數i增加，乘法次數提高到0.05%，達到了百倍的提升。

3.2 試驗平臺分析

為直觀顯示改進算法的特性，本次采用車載GNSS動態測量試驗數據。2017年6月用3臺SinoGNSS M600 GNSS接收機進行動態測量試驗，如圖2所示，其中一臺接收機放置在基準站Location1，其他2號和3號放在速度為80 km/h的汽車車內(固定位置不變)。

圖2 試驗平臺

本文采用CPU為2.9 GHz的聯想計算機進行Matlab編程，對比兩算法解算耗時分析。表1列出了不同的動態GNSS應用中，兩算法的Kalman濾波循環一萬次取平均值。

表1 不同參數下CPU運行時間

T1、T2分別代表普通算法和改進算法的Kalman濾波一次迭代所用時間，改進算法所用耗時均低于普通算法10%，最低可以達到4.31%。實際CPU處理器采用四級流水線技術，將幾條指令并行處理，其對Kalman濾波算法總體上加快了指令流速度，縮短了程序執行時間，故實際Kalman濾波運算耗時高于一次迭代測量值T2/T1，故表1僅概略反應算法效率。

4 結 論

Kalman濾波迭代收斂過程中，經處理器的流水線技術處理，使實際計算耗時高于測量值。因此表1為估計數值，但可以確定CPU計算耗時低于普通算法的10%。該改進算法乘法次數降低到1%，故實現了算法的高效性。在動態GNSS數據處理中，通常Kalman濾波算法會采取適當優化，本文基于非優化算法求解雙差整周模糊度，尤其在多系統產生的可見衛星數目增加情況下，有一定的參考價值。

致謝：上海司南衛星導航技術股份有限公司對本文的工作提供了試驗平臺和技術支持，在此表示衷心感謝。

參考文獻：

[1] 任曉東,張柯柯,李星星,等.BeiDou、Galileo、GLONASS、GPS多系統融合精密單點[J].測繪學報,2015，44(12):1307-1313.

[2] MOURI A，KUBO Y，SUGIMOTO S.Detection and Correction of Doppler Biases in Kalman Filter-based Positioning[J].Proceedings of the ISCIE International Symposium on Stochastic Systems Theory and its Applications,2015(10):156-164.

[3] WANG D，LV H，WU J.Augmented cubature Kalman Filter for Nonlinear RTK/MIMU Integrated Navigation with Non-additive Noise[J].Measurement，2017,97(2):111-125.

[4] 趙占祥,李興國,婁國偉,等.矩陣外積法Kalman濾波器在動態GPS定位中的應用[J].火力與指揮控制,2006，31(4):27-29.

[5] 胡朱林.基于矩陣分解的卡爾曼濾波技術分析及應用[J].數字技術與應用,2016，(11):222.

[6] 郭樹人,郭海榮,何海波,等.GPS動態數據處理中的快速Kalman濾波算法[J].測繪科學技術學報,2006,23(3):171-173.

[7] SAHMOUDI M，LANDRY R.A Nonlinear Filtering Approach for Robust Multi-GNSS RTK Positioning in Presence of Multipath and Ionospheric Delays[J].IEEE Journal of Selected Topics in Signal Processing,2009,3(5):764-776.

[8] LI B，SHEN Y，FENG Y.GNSS Ambiguity Resolution with Controllable Failure Rate for Long Baseline Network RTK[J].Journal of Geodesy,2014,88(2):99-112.

[9] 孫紅星,李德仁.使用雙頻相關法單歷元解算GPS整周模糊值[J].測繪學報，2003,32(3):208-212.

[10] 董緒榮,陶大欣.一個快速Kalman濾波方法及其在GPS動態數據處理中的應用[J].測繪學報,1997(3):35-41.

[11] PROCHNIEWICZ D，SZPUNAR R，WALO J.A New Study of Describing the Reliability of GNSS Network RTK Positioning with the Use of Quality Indicators[J].Measurement Science and Technology,2017,28(1):12-15.

[12] 郭海榮,楊元喜,何海波,等.導航衛星原子鐘Kalman濾波中噪聲方差-協方差的確定[J].測繪學報,2010,39(2):146-150.

[13] Bartonek D，BURES J，SVABENSKY O.Optimization of Process Field Measurement GNSS-RTK for Railway Infrastructure[J].Solid State Phenomena，2017,25(8):481-484.

[14] PAZIEWSKI J，SIERADZKI R.Integrated GPS+BDS Instantaneous Medium Baseline RTK Positioning：Signal Analysis，Methodology and Performance Assessment[J].Advances in Space Research，2017，60(12):2561-2573.

[15] PAN V Y.Nearly Optimal Computations with Structured Matrices[J].Theoretical Computer Science，2017：953-962.