范紅平，周志峰，王永泉

(1. 上海工程技術(shù)大學(xué)機(jī)械工程學(xué)院，上海 201620； 2. 上海司南衛(wèi)星導(dǎo)航技術(shù)股份有限公司，上海 201801)

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

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

1 GNSS RTK的Kalman濾波模型

1.1 線性高斯系統(tǒng)的Kalman濾波方程

建立以下Kalman濾波的系統(tǒng)動(dòng)態(tài)方程，式(1)、式(2)為一步預(yù)測(cè)過(guò)程，得出下一時(shí)刻最優(yōu)預(yù)測(cè)值[8]。式(3)、式(4)和式(5)為一步更新過(guò)程，先計(jì)算出卡爾曼增益K，然后計(jì)算狀態(tài)估計(jì)值。

(1)

(2)

(3)

(4)

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

(5)

1.2 GNSS RTK的雙差觀測(cè)方程

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

(6)

(7)

1.3 優(yōu)化Kalman濾波運(yùn)算的必要性

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

2 改進(jìn)型Kalman濾波方法

2.1 改進(jìn)1：狀態(tài)轉(zhuǎn)移矩陣F模型

本文取三維直角坐標(biāo)構(gòu)建定常加速度動(dòng)態(tài)模型，假設(shè)載體處于穩(wěn)定加速狀態(tài)，則狀態(tài)參數(shù)需要增加3個(gè)加速度分量，則9+3i維狀態(tài)變量為

(8)

(9)

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

2.2 改進(jìn)2：協(xié)方差矩陣Q計(jì)算

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

(10)

3 算例及分析

3.1 理論分析

(11)

(12)

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

圖1 改進(jìn)算法優(yōu)化比率

如圖1所示，可見(jiàn)衛(wèi)星多于4顆，其改進(jìn)算法乘法次數(shù)H1均低于普通矩陣H2的1.2%，且隨可見(jiàn)衛(wèi)星數(shù)i增加，乘法次數(shù)提高到0.05%，達(dá)到了百倍的提升。

3.2 試驗(yàn)平臺(tái)分析

為直觀顯示改進(jìn)算法的特性，本次采用車載GNSS動(dòng)態(tài)測(cè)量試驗(yàn)數(shù)據(jù)。2017年6月用3臺(tái)SinoGNSS M600 GNSS接收機(jī)進(jìn)行動(dòng)態(tài)測(cè)量試驗(yàn)，如圖2所示，其中一臺(tái)接收機(jī)放置在基準(zhǔn)站Location1，其他2號(hào)和3號(hào)放在速度為80 km/h的汽車車內(nèi)(固定位置不變)。

圖2 試驗(yàn)平臺(tái)

本文采用CPU為2.9 GHz的聯(lián)想計(jì)算機(jī)進(jìn)行Matlab編程，對(duì)比兩算法解算耗時(shí)分析。表1列出了不同的動(dòng)態(tài)GNSS應(yīng)用中，兩算法的Kalman濾波循環(huán)一萬(wàn)次取平均值。

表1 不同參數(shù)下CPU運(yùn)行時(shí)間

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

4 結(jié) 論

Kalman濾波迭代收斂過(guò)程中，經(jīng)處理器的流水線技術(shù)處理，使實(shí)際計(jì)算耗時(shí)高于測(cè)量值。因此表1為估計(jì)數(shù)值，但可以確定CPU計(jì)算耗時(shí)低于普通算法的10%。該改進(jìn)算法乘法次數(shù)降低到1%，故實(shí)現(xiàn)了算法的高效性。在動(dòng)態(tài)GNSS數(shù)據(jù)處理中，通常Kalman濾波算法會(huì)采取適當(dāng)優(yōu)化，本文基于非優(yōu)化算法求解雙差整周模糊度，尤其在多系統(tǒng)產(chǎn)生的可見(jiàn)衛(wèi)星數(shù)目增加情況下，有一定的參考價(jià)值。

致謝：上海司南衛(wèi)星導(dǎo)航技術(shù)股份有限公司對(duì)本文的工作提供了試驗(yàn)平臺(tái)和技術(shù)支持，在此表示衷心感謝。

參考文獻(xiàn)：

[1] 任曉東,張柯柯,李星星,等.BeiDou、Galileo、GLONASS、GPS多系統(tǒng)融合精密單點(diǎn)[J].測(cè)繪學(xué)報(bào),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] 趙占祥,李興國(guó),婁國(guó)偉,等.矩陣外積法Kalman濾波器在動(dòng)態(tài)GPS定位中的應(yīng)用[J].火力與指揮控制,2006，31(4):27-29.

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

[6] 郭樹(shù)人,郭海榮,何海波,等.GPS動(dòng)態(tài)數(shù)據(jù)處理中的快速Kalman濾波算法[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào),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，F(xiàn)ENG Y.GNSS Ambiguity Resolution with Controllable Failure Rate for Long Baseline Network RTK[J].Journal of Geodesy,2014,88(2):99-112.

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

[10] 董緒榮,陶大欣.一個(gè)快速Kalman濾波方法及其在GPS動(dòng)態(tài)數(shù)據(jù)處理中的應(yīng)用[J].測(cè)繪學(xué)報(bào),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] 郭海榮,楊元喜,何海波,等.導(dǎo)航衛(wèi)星原子鐘Kalman濾波中噪聲方差-協(xié)方差的確定[J].測(cè)繪學(xué)報(bào),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.