基于Kalman濾波的大規模GNSS網參數估計方法

王 崢，周 劍

(武漢大學測繪學院，湖北 武漢 430079)

基于Kalman濾波的大規模GNSS網參數估計方法

王 崢，周 劍

(武漢大學測繪學院，湖北 武漢 430079)

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

Kalman濾波；大規模GNSS網；系統誤差；基準定義

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

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

本文以Kalman濾波理論為基礎，研究濾波算法在大規模GNSS網參數(測站坐標、衛星軌道、極移參數等)估計中的作用，并根據實測數據對方法進行驗證。結果表明，利用Kalman濾波理論可以高精度、高效率地進行GNSS網各類參數的估計。

1 參數估計方法

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

假設δlk=lk-l0k為tk時刻線性化后的準觀測值向量。其中l0k為外部提供的參數先驗信息；lk為基線解算階段所得參數估值。則觀測方程為

δlk=Akδxk+εk

(1)

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

δxk+1=Skδxk+qk

(2)

式中，Sk為狀態轉移矩陣；qk為方差為Qk的隨機過程噪聲向量。為估計tk+1時刻的狀態，首先預測tk+1時刻的參數估值及其方差-協方差陣

(3)

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

然后利用虛擬觀測值對預測的參數狀態進行修正

(4)

其中

(5)

上述即為利用Kalman濾波進行GNSS網參數估計的過程，其流程如圖1所示。

圖1 基于Kalman濾波的參數估計流程

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

為消除各時段/子網解之間的系統誤差，需在濾波函數模型中通過擴充狀態空間引入系統誤差參數。系統誤差可分為位置基準的系統誤差、時間演變基準的系統誤差、尺度基準的系統誤差和空間方位基準的系統誤差4類。對于粗差可用殘差加權平方和進行探測。殘差的加權平方和定義為

(6)

其中

(7)

濾波過程中每步均計算殘差加權平方和，若其超限則說明基線解算階段所得參數估值中可能含有粗差，或對參數的先驗約束過緊。經驗表明，殘差加權平方和的限差設為100較為適宜。

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

2 數據分析

為驗證上述理論結果，本文利用實測數據進行了試驗。文中基線解算采用GAMIT軟件完成。

2.1 方案設計

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

圖2 測站分布

基線解算階段的參數設置見表1。

基線解算時為防止先驗信息不準確引入誤差，對各類參數的約束均較松，導致結果整體存在平移、旋轉和縮放，在濾波時需根據參數的先驗信息確定基準。濾波以基線解算所得100個單日松弛解為對象，將其作為準觀測值輸入濾波器。

2.2 結果分析

本文主要通過Kalman濾波實現測站位置參數、衛星軌道參數、極移參數等的估計。

估計測站位置參數時以兩個子網的100個單日松弛解為輸入。由于對測站位置的估計可利用Helmert轉換的方式引入基準，故Kalman濾波時測站位置、衛星軌道、ERP參數的先驗約束均松弛，以濾波不發散為目的。Helmert轉換中用于基準定義的測站為圖2中11個公共站點，基準測站坐標采用IGS公布的結果。本文以IGS結果為準檢驗其余測站位置估值精度。

圖3為各測站坐標分量估值與IGS結果之差(不包括11個公共站點)。由圖3可看出，本文對測站位置參數的估值精度較高；各測站坐標分量與IGS結果的較差基本在2 cm以內。統計結果顯示，各分量較差均值分別為0.24、0.35、-0.14 cm，RMS值分別為0.85、1.1、1.21 cm。

圖3 測站坐標估值與IGS結果之差

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

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

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

圖5 衛星軌道估值與IGS精密星歷之差

ERP參數是不斷變化的，因此本文每天估計一套對應于測段中間歷元的ERP參數。在估計ERP參數時需對測站位置及衛星軌道緊約束。同時，由于日長變化(length of day，LOD)與衛星軌道的升交點赤經Ω強相關，因此利用GNSS技術不能實現LOD的直接估計。本文僅針對極移參數進行分析。

圖6為極移估值與IERS所公布結果的較差。由圖6可知，本文所得極移參數估值與IERS結果之差基本在0.1 mas以內，約對應地球表面赤道處3.1 mm的位移，精度較高。但由圖6可知估值存在系統誤差，原因尚需進一步分析。

圖6 極移參數估值與IERS結果之差

3 結束語

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

[1] 李征航, 張小紅. 衛星導航定位新技術及高精度數據處理方法[M]. 武漢: 武漢大學出版社, 2009.

[2] 楊凱. 大規模GNSS基準站網數據處理關鍵技術研究與實現[D]. 武漢: 武漢大學, 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] 施闖. 大規模高精度GPS網平差與分析理論及其應用[M]. 北京: 測繪出版社, 2002.

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

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

[8] 付夢印, 鄧志紅, 閆莉萍. Kalman濾波理論及其在導航系統中的應用[M]. 2版. 北京:科學出版社, 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網參數估計中基準定義方法研究[J]. 大地測量與地球動力學, 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濾波的大規模GNSS網參數估計方法[J].測繪通報，2017(1)：18-21.

10.13474/j.cnki.11-2246.2017.0004.

2016-04-18

國家重點基礎研究發展計劃(973計劃)(2013CB733301)；國家自然科學基金(41210006；41374022)

王 崢(1986—)，女，博士生，研究方向為衛星大地測量。E-mail: zhengwang@whu.edu.cn