王 崢，周 劍

(武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079)

基于Kalman濾波的大規(guī)模GNSS網(wǎng)參數(shù)估計(jì)方法

王 崢，周 劍

(武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079)

針對(duì)傳統(tǒng)方法存在的缺陷，研究了利用Kalman濾波技術(shù)進(jìn)行大規(guī)模GNSS網(wǎng)參數(shù)(主要包括測(cè)站位置參數(shù)、衛(wèi)星軌道參數(shù)及極移參數(shù))估計(jì)的理論方法與關(guān)鍵技術(shù)，并利用40個(gè)全球均勻分布的IGS站多天的觀測(cè)數(shù)據(jù)對(duì)理論成果進(jìn)行了驗(yàn)證。結(jié)果表明，本文估計(jì)得到的測(cè)站位置參數(shù)與IGS結(jié)果各分量較差的RMS值分別為0.85、1.1、1.21 cm，得到的衛(wèi)星軌道參數(shù)外推1 h后與IGS最終星歷各分量較差的RMS值分別為9.8、8.6、7.2 cm，得到的極移參數(shù)與IERS結(jié)果的較差基本在1 mas之內(nèi)；該方法具有較高的估值精度，可有效地用于GNSS網(wǎng)各類參數(shù)的估計(jì)。

Kalman濾波；大規(guī)模GNSS網(wǎng)；系統(tǒng)誤差；基準(zhǔn)定義

隨著GNSS技術(shù)的發(fā)展與完善，多個(gè)國家、地區(qū)、行業(yè)紛紛建立了滿足自身需求的連續(xù)運(yùn)行GNSS基準(zhǔn)站網(wǎng)[1]。隨著基準(zhǔn)站網(wǎng)數(shù)量的不斷增加，其數(shù)據(jù)處理方法也越來越引起人們的重視。如何形成一套有效的數(shù)據(jù)處理策略，在保證解算精度與可靠性的前提下節(jié)省解算時(shí)間和存儲(chǔ)空間，一直是人們所研究的重點(diǎn)[2]。

將全部測(cè)站、所有時(shí)段的觀測(cè)數(shù)據(jù)一并處理雖可保證數(shù)學(xué)模型的嚴(yán)密性，但仍存在諸多難點(diǎn)。如目前全球公認(rèn)的高精度GNSS數(shù)據(jù)解算軟件GAMIT、GIPSY、Bernese、EPOS等中，GAMIT只能解算不超過100個(gè)測(cè)站的數(shù)據(jù)[3]，Bernese、EPOS等雖可同時(shí)解算100個(gè)以上測(cè)站的數(shù)據(jù)，但解算時(shí)間過長，占用內(nèi)存較多。為解決上述問題，目前常用的參數(shù)估計(jì)方法一般分3步完成：①將整個(gè)觀測(cè)網(wǎng)分為若干子網(wǎng)，將各個(gè)子網(wǎng)、各個(gè)時(shí)段的觀測(cè)數(shù)據(jù)單獨(dú)解算得到松約束下的參數(shù)估值及其方差-協(xié)方差矩陣；②結(jié)合所有松約束結(jié)果得到站坐標(biāo)、速度場(chǎng)、地球自轉(zhuǎn)參數(shù)及衛(wèi)星軌道等的統(tǒng)一估值；③根據(jù)外部信息引入基準(zhǔn)，得到參數(shù)最終估值[4]。此算法大幅節(jié)省了計(jì)算時(shí)間，也可最大限度地保證理論嚴(yán)密。上述步驟②一般通過法方程疊加完成。法方程疊加的本質(zhì)為序貫最小二乘，理論嚴(yán)密，且省略諸多矩陣求逆運(yùn)算，效率較高。以基線向量為虛擬觀測(cè)值，施闖研究了基于法方程疊加的大規(guī)模GNSS網(wǎng)平差方法，并編制了相關(guān)軟件[5]。姚宜斌研究了坐標(biāo)模式的法方程疊加方法，推導(dǎo)了相應(yīng)公式，對(duì)各類參數(shù)的估計(jì)方法進(jìn)行了總結(jié)[6-7]。但法方程疊加法對(duì)狀態(tài)參數(shù)的處理較復(fù)雜。起源于20世紀(jì)60年代的Kalman濾波理論將狀態(tài)空間的概念引入到隨機(jī)估計(jì)問題中，把信號(hào)過程視為白噪聲作用下線性系統(tǒng)的輸出，并用狀態(tài)方程描述這種關(guān)系，可方便處理狀態(tài)參數(shù)，且不需存儲(chǔ)中間過程，大幅節(jié)省了存儲(chǔ)空間，可有效地應(yīng)用于GNSS基準(zhǔn)站網(wǎng)數(shù)據(jù)處理[8]。

本文以Kalman濾波理論為基礎(chǔ)，研究濾波算法在大規(guī)模GNSS網(wǎng)參數(shù)(測(cè)站坐標(biāo)、衛(wèi)星軌道、極移參數(shù)等)估計(jì)中的作用，并根據(jù)實(shí)測(cè)數(shù)據(jù)對(duì)方法進(jìn)行驗(yàn)證。結(jié)果表明，利用Kalman濾波理論可以高精度、高效率地進(jìn)行GNSS網(wǎng)各類參數(shù)的估計(jì)。

1 參數(shù)估計(jì)方法

采用Kalman濾波算法進(jìn)行參數(shù)估計(jì)，首先應(yīng)根據(jù)未知參數(shù)的先驗(yàn)信息及狀態(tài)方程預(yù)測(cè)下一時(shí)刻的參數(shù)估值及其方差-協(xié)方差矩陣；再以下一時(shí)刻所有子網(wǎng)在基線解算階段所得參數(shù)估值及其方差-協(xié)方差陣為虛擬觀測(cè)值對(duì)預(yù)測(cè)值進(jìn)行修正，并根據(jù)修正后的參數(shù)信息預(yù)測(cè)下時(shí)刻的狀態(tài)。該過程循環(huán)漸進(jìn)，依次進(jìn)行，直至所有信息均考慮在內(nèi)。利用Kalman濾波算法結(jié)合多個(gè)子網(wǎng)的松約束結(jié)果與結(jié)合多個(gè)時(shí)段的解算結(jié)果方法類似，實(shí)際上前者可看作是后者在時(shí)段間隔為0時(shí)的特殊情況。

假設(shè)δlk=lk-l0k為tk時(shí)刻線性化后的準(zhǔn)觀測(cè)值向量。其中l(wèi)0k為外部提供的參數(shù)先驗(yàn)信息；lk為基線解算階段所得參數(shù)估值。則觀測(cè)方程為

δlk=Akδxk+εk

(1)

式中，δxk為待估參數(shù)向量，包含所有測(cè)站位置參數(shù)、所有參與解算衛(wèi)星軌道參數(shù)，以及各測(cè)段中間歷元的ERP參數(shù)及其線性速度；Ak為設(shè)計(jì)矩陣；εk是方差為Pk的零均值噪聲向量。狀態(tài)方程為

δxk+1=Skδxk+qk

(2)

式中，Sk為狀態(tài)轉(zhuǎn)移矩陣；qk為方差為Qk的隨機(jī)過程噪聲向量。為估計(jì)tk+1時(shí)刻的狀態(tài)，首先預(yù)測(cè)tk+1時(shí)刻的參數(shù)估值及其方差-協(xié)方差陣

(3)

式中，Ck為基線解算階段所得參數(shù)估值的方差-協(xié)方差矩陣。

然后利用虛擬觀測(cè)值對(duì)預(yù)測(cè)的參數(shù)狀態(tài)進(jìn)行修正

(4)

其中

(5)

上述即為利用Kalman濾波進(jìn)行GNSS網(wǎng)參數(shù)估計(jì)的過程，其流程如圖1所示。

圖1 基于Kalman濾波的參數(shù)估計(jì)流程

在Kalman濾波過程中，濾波起始階段的先驗(yàn)信息對(duì)于濾波是否收斂及收斂速度具有重要影響。在GNSS網(wǎng)參數(shù)估計(jì)中，一般選取基線解算階段所得參數(shù)估值為先驗(yàn)值啟動(dòng)濾波器。濾波時(shí)需在狀態(tài)方程中加入隨機(jī)過程噪聲以削弱諸如測(cè)站震后蠕變、天線高量測(cè)誤差、ERP參數(shù)變化等的影響。隨機(jī)游走過程(random walk process，RWP)較多地被應(yīng)用其中。測(cè)站坐標(biāo)的隨機(jī)游走噪聲參數(shù)需根據(jù)各測(cè)站數(shù)據(jù)采集過程中的具體情況確定。根據(jù)IERS的結(jié)果，極移的隨機(jī)游走噪聲方差一般可設(shè)為22.8 mas2/a，表示允許極移參數(shù)每天變動(dòng)0.25 mas。結(jié)果表明，允許時(shí)段之間參數(shù)的隨機(jī)變化可有效吸收模型誤差及測(cè)站位置變動(dòng)等對(duì)估計(jì)結(jié)果的影響，提高狀態(tài)參數(shù)的估值精度[9]。

為消除各時(shí)段/子網(wǎng)解之間的系統(tǒng)誤差，需在濾波函數(shù)模型中通過擴(kuò)充狀態(tài)空間引入系統(tǒng)誤差參數(shù)。系統(tǒng)誤差可分為位置基準(zhǔn)的系統(tǒng)誤差、時(shí)間演變基準(zhǔn)的系統(tǒng)誤差、尺度基準(zhǔn)的系統(tǒng)誤差和空間方位基準(zhǔn)的系統(tǒng)誤差4類。對(duì)于粗差可用殘差加權(quán)平方和進(jìn)行探測(cè)。殘差的加權(quán)平方和定義為

(6)

其中

(7)

濾波過程中每步均計(jì)算殘差加權(quán)平方和，若其超限則說明基線解算階段所得參數(shù)估值中可能含有粗差，或?qū)?shù)的先驗(yàn)約束過緊。經(jīng)驗(yàn)表明，殘差加權(quán)平方和的限差設(shè)為100較為適宜。

若在濾波起始階段對(duì)參數(shù)的約束較松，則濾波完成后得到的結(jié)果是一無約束或松約束網(wǎng)形，需引入外部信息，以確定整網(wǎng)基準(zhǔn)[10]。對(duì)測(cè)站位置基準(zhǔn)的定義可通過HERMERT參數(shù)轉(zhuǎn)換的方法實(shí)現(xiàn)。該方法最大限度地保留了GNSS觀測(cè)的原始精度，防止引入的基準(zhǔn)信息誤差破壞參數(shù)估計(jì)結(jié)果。已有研究表明，利用Helmert參數(shù)轉(zhuǎn)換的方法定義基準(zhǔn)在只估計(jì)3個(gè)旋轉(zhuǎn)參數(shù)的情況下GNSS網(wǎng)形不發(fā)生任何變化；若同時(shí)顧及3個(gè)平移參數(shù)會(huì)導(dǎo)致網(wǎng)形畸變，但對(duì)上萬千米的基線，畸變量在3 mm以內(nèi)[11]。

2 數(shù)據(jù)分析

為驗(yàn)證上述理論結(jié)果，本文利用實(shí)測(cè)數(shù)據(jù)進(jìn)行了試驗(yàn)。文中基線解算采用GAMIT軟件完成。

2.1 方案設(shè)計(jì)

本文采用的數(shù)據(jù)為2013年1月1日(年積日1)至2月19日(年積日50)全球均勻分布的40個(gè)IGS跟蹤站的觀測(cè)數(shù)據(jù)，站點(diǎn)分布如圖2所示。文中將40個(gè)跟蹤站分為兩個(gè)子網(wǎng)分別進(jìn)行基線解算。子網(wǎng)1包含26個(gè)測(cè)站(圖2中灰色點(diǎn)與黑色點(diǎn))；子網(wǎng)2包含25個(gè)測(cè)站(圖2中灰色點(diǎn)與淺灰色點(diǎn))。兩個(gè)子網(wǎng)擁有11個(gè)公共站點(diǎn)(圖2中灰色點(diǎn))。公共站點(diǎn)基本均勻分布于全球，用于濾波時(shí)聯(lián)結(jié)兩個(gè)子網(wǎng)的解算結(jié)果。

圖2 測(cè)站分布

基線解算階段的參數(shù)設(shè)置見表1。

基線解算時(shí)為防止先驗(yàn)信息不準(zhǔn)確引入誤差，對(duì)各類參數(shù)的約束均較松，導(dǎo)致結(jié)果整體存在平移、旋轉(zhuǎn)和縮放，在濾波時(shí)需根據(jù)參數(shù)的先驗(yàn)信息確定基準(zhǔn)。濾波以基線解算所得100個(gè)單日松弛解為對(duì)象，將其作為準(zhǔn)觀測(cè)值輸入濾波器。

2.2 結(jié)果分析

本文主要通過Kalman濾波實(shí)現(xiàn)測(cè)站位置參數(shù)、衛(wèi)星軌道參數(shù)、極移參數(shù)等的估計(jì)。

估計(jì)測(cè)站位置參數(shù)時(shí)以兩個(gè)子網(wǎng)的100個(gè)單日松弛解為輸入。由于對(duì)測(cè)站位置的估計(jì)可利用Helmert轉(zhuǎn)換的方式引入基準(zhǔn)，故Kalman濾波時(shí)測(cè)站位置、衛(wèi)星軌道、ERP參數(shù)的先驗(yàn)約束均松弛，以濾波不發(fā)散為目的。Helmert轉(zhuǎn)換中用于基準(zhǔn)定義的測(cè)站為圖2中11個(gè)公共站點(diǎn)，基準(zhǔn)測(cè)站坐標(biāo)采用IGS公布的結(jié)果。本文以IGS結(jié)果為準(zhǔn)檢驗(yàn)其余測(cè)站位置估值精度。

圖3為各測(cè)站坐標(biāo)分量估值與IGS結(jié)果之差(不包括11個(gè)公共站點(diǎn))。由圖3可看出，本文對(duì)測(cè)站位置參數(shù)的估值精度較高；各測(cè)站坐標(biāo)分量與IGS結(jié)果的較差基本在2 cm以內(nèi)。統(tǒng)計(jì)結(jié)果顯示，各分量較差均值分別為0.24、0.35、-0.14 cm，RMS值分別為0.85、1.1、1.21 cm。

圖3 測(cè)站坐標(biāo)估值與IGS結(jié)果之差

估計(jì)衛(wèi)星軌道參數(shù)時(shí)可以廣播星歷為先驗(yàn)值。為防止軌道估計(jì)的“末端效應(yīng)”，本文利用年積日1日至3日連續(xù)3天共6個(gè)單日松弛解文件估計(jì)衛(wèi)星軌道，并取2日的軌道與IGS軌道對(duì)比以檢驗(yàn)參數(shù)估計(jì)精度。在基線解算階段所有3天的衛(wèi)星軌道參數(shù)均為中間一天中間時(shí)刻的密切元素以及太陽光壓、Y軸偏差等。在濾波過程中緊約束測(cè)站位置及ERP參數(shù)，以提供基準(zhǔn)。測(cè)站坐標(biāo)的先驗(yàn)值采用IGS的結(jié)果，ERP參數(shù)的先驗(yàn)值采用IERS的結(jié)果。

圖4為作為初始值的廣播星歷與IGS軌道在參考時(shí)刻之差。圖5為估計(jì)的衛(wèi)星軌道參數(shù)外推1小時(shí)得到的衛(wèi)星位置與IGS軌道之差。與IGS精密星歷相比，本文所得軌道在3個(gè)方向的較差均值分別為-0.6、-0.9、-1.1 cm；各分量較差的RMS值分別為9.8、8.6、7.2 cm。作為初始值的廣播星歷3個(gè)方向的較差均值分別為5.2、-5.1、11.5 cm；各分量較差的RMS值分別為66.6、77.8、60.7 cm。由此可知，利用Kalman濾波，融合多天觀測(cè)數(shù)據(jù)以估計(jì)衛(wèi)星軌道的方法大幅提高了軌道精度。本文的主要目的為方法驗(yàn)證，因此采用的測(cè)站數(shù)目與觀測(cè)時(shí)間有限。可以預(yù)見，如果利用更多測(cè)站、更長的觀測(cè)時(shí)間估計(jì)衛(wèi)星軌道，應(yīng)該可以得到更好的結(jié)果。

圖4 廣播星歷與IGS精密星歷之差

圖5 衛(wèi)星軌道估值與IGS精密星歷之差

ERP參數(shù)是不斷變化的，因此本文每天估計(jì)一套對(duì)應(yīng)于測(cè)段中間歷元的ERP參數(shù)。在估計(jì)ERP參數(shù)時(shí)需對(duì)測(cè)站位置及衛(wèi)星軌道緊約束。同時(shí)，由于日長變化(length of day，LOD)與衛(wèi)星軌道的升交點(diǎn)赤經(jīng)Ω強(qiáng)相關(guān)，因此利用GNSS技術(shù)不能實(shí)現(xiàn)LOD的直接估計(jì)。本文僅針對(duì)極移參數(shù)進(jìn)行分析。

圖6為極移估值與IERS所公布結(jié)果的較差。由圖6可知，本文所得極移參數(shù)估值與IERS結(jié)果之差基本在0.1 mas以內(nèi)，約對(duì)應(yīng)地球表面赤道處3.1 mm的位移，精度較高。但由圖6可知估值存在系統(tǒng)誤差，原因尚需進(jìn)一步分析。

圖6 極移參數(shù)估值與IERS結(jié)果之差

3 結(jié)束語

隨著GNSS技術(shù)的發(fā)展，其逐漸成為獲取各類參數(shù)的重要手段。本文主要研究了利用Kalman濾波結(jié)合多個(gè)子網(wǎng)、多個(gè)時(shí)段的松弛解文件進(jìn)行大規(guī)模GNSS網(wǎng)參數(shù)估計(jì)的理論，并通過實(shí)測(cè)數(shù)據(jù)對(duì)方法進(jìn)行了驗(yàn)證。結(jié)果表明，利用Kalman濾波技術(shù)可進(jìn)行地面測(cè)站位置參數(shù)、導(dǎo)航衛(wèi)星軌道參數(shù)及地球自轉(zhuǎn)參數(shù)(主要為極移參數(shù))的高精度估計(jì)，從而為利用GNSS技術(shù)估計(jì)各類參數(shù)提供了一種有效的技術(shù)手段，可以在科研生產(chǎn)中推廣應(yīng)用。

[1] 李征航, 張小紅. 衛(wèi)星導(dǎo)航定位新技術(shù)及高精度數(shù)據(jù)處理方法[M]. 武漢: 武漢大學(xué)出版社, 2009.

[2] 楊凱. 大規(guī)模GNSS基準(zhǔn)站網(wǎng)數(shù)據(jù)處理關(guān)鍵技術(shù)研究與實(shí)現(xiàn)[D]. 武漢: 武漢大學(xué), 2011.

[3] MIT. GAMIT Reference Manual Release 10.4[EB/OL]. [2014-11-15].http:∥www-gpsg.mit.edu/～simon/gtgk/GAMIT_Ref.pdf.

[4] DONG D, HERRING T A, KING R W. Estimating Regional Deformation from a Combination of Space and Terrestrial Geodetic Data[J]. Journal of Geodesy, 1998, 72(4): 200-214.

[5] 施闖. 大規(guī)模高精度GPS網(wǎng)平差與分析理論及其應(yīng)用[M]. 北京: 測(cè)繪出版社, 2002.

[6] 姚宜斌, 劉經(jīng)南, 陶本藻, 等. 基于坐標(biāo)模式的廣義網(wǎng)平差模型研究[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2006, 31(1):16-18.

[7] 姚宜斌. GPS精密定位定軌后處理[M]. 北京: 測(cè)繪出版社, 2008.

[8] 付夢(mèng)印, 鄧志紅, 閆莉萍. Kalman濾波理論及其在導(dǎo)航系統(tǒng)中的應(yīng)用[M]. 2版. 北京:科學(xué)出版社, 2010.

[9] HERRING T A, DAVIS J L, SHAPRIO I I. Geodesy by Radio Interferometry: The Application of Kalman Filtering to the Analysis of Very Long Baseline Interferometry Data[J]. Journal of Geophysical Research, 1990, 95(B8): 12561-12581.

[10] HEFLIN M, BERTIGER W, BLEWITT G, et al. Global Geodesy Using GPS without Fiducial Sites[J]. Geophysical Research Letters, 1992, 19(2): 131-134.

[11] 肖玉鋼, 劉鴻飛, 王崢. GNSS網(wǎng)參數(shù)估計(jì)中基準(zhǔn)定義方法研究[J]. 大地測(cè)量與地球動(dòng)力學(xué), 2012, 32(4):79-82.

[12] DONG D, BOCK Y. Global Positioning System Network Analysis with Phase Ambiguity Resolution Applied to Crustal Deformation Studies in California[J]. Journal of Geophysical Research, 1989, 94(B4): 3949-3966.

[13] SCHMID R, ROTHACHER M, THALLER D, et al. Absolute Phase Center Corrections of Satellite and Receiver Antennas[J]. GPS Solutions, 2005, 9(4): 283-293.

ResearchandCaseAnalysisofEstimationMethodBasedonKalmanFilteringforGeodeticParametersfromGNSSNetworks

WANG Zheng,ZHOU Jian

(School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China)

Focusing on the shortcomings of traditional methods, theory and key technologies with using Kalman filtering were researched for estimating geodetic parameters from GNSS network, such as station coordinates, satellite orbits and polar motion parameters. Then the achievements were verified with days of observations from 40 globally distributed IGS stations. The results indicated that the RMS values of the difference between IGS station coordinates estimated here and that advised by IGS were 0.85, 1.1 and 1.21 cm inX,Y,Zdirection respectively. The RMS values achieved in the similar way for the IGS final orbits and that acquired by extrapolating for an hour with the orbit parameters estimated here were 9.8, 8.6 and 7.2 cm. The difference between the estimated polar motion parameters and those advised by IERS were mostly within 1 mas. The results showed that the algorithm used here was precise and could be utilized to estimate the geodetic parameters from GNSS networks.

Kalman filtering; large scale GNSS network; systematic error; datum definition

P228

A

0494-0911(2017)01-0018-04

王崢，周劍.基于Kalman濾波的大規(guī)模GNSS網(wǎng)參數(shù)估計(jì)方法[J].測(cè)繪通報(bào)，2017(1)：18-21.

10.13474/j.cnki.11-2246.2017.0004.

2016-04-18

國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)(2013CB733301)；國家自然科學(xué)基金(41210006；41374022)

王 崢(1986—)，女，博士生，研究方向?yàn)樾l(wèi)星大地測(cè)量。E-mail: zhengwang@whu.edu.cn