附加歷元間約束的滑動窗單頻實時精密單點定位算法

楊凱淳，呂志平,2，李林陽,3，鄺英才，許 煒，鄭 茜

1. 信息工程大學地理空間信息學院，河南 鄭州 450001； 2. 哈爾濱工業大學(深圳)空間科學與應用技術研究院，廣東 深圳 518055； 3. 中國科學院精密測量科學與技術創新研究院大地測量與地球動力學國家重點實驗室，湖北 武漢 430077； 4. 西安科技大學測繪科學與技術學院，陜西 西安 710054

單頻接收機由于價格低廉、操作方便，已廣泛應用于變形監測、低軌衛星定軌及地面網解算等領域[1-3]，然而單頻精密單點定位(precise point positioning,PPP)的定位精度與多頻PPP相比，還存在一定差距，最大的難點在于電離層延遲的改正[4-6]。不少學者對其開展了很多研究，得到了有效的處理方法[7-9]。

目前，主要采取以下4種方法改正單頻PPP的電離層延遲：①采用電離層模型，但利用GPS導航電文中的Klobuchar模型，只能消除50%～60%的電離層延遲，而IGS格網數據精度為2 TECU，在電離層活躍及IGS站很少的區域其精度更低[9-10]；②構成半合觀測值，該方法只消除了電離層延遲的一階項，且由于模糊度與接收機鐘差高度相關，無法獲得真實的鐘差估值[11-12]；③將電離層斜延遲作為未知參數估計[13-16]，目前普遍采用的是文獻[13]中的三參數法，同時估計對天頂延遲、水平梯度、映射函數，但模型過于簡單，無法較好地描述電離層的活動，特別是在電離層延遲劇烈變化的區域；④通過多頻接收機建立區域電離層模型，對單頻接收機進行改正[17-18]，該方法在距離較遠且條件不好時，不方便使用。若不使用電離層模型也不擬合電離層斜延遲，而是直接將其作為參數求解，不僅可以獲取高精度的電離層延遲信息與測站坐標，還能在一定程度上提高收斂速度[14,19]。

單歷元單頻PPP算法在同時解算電離層延遲時會出現秩虧問題[11,15,19]，而使用多歷元聯合求解足以克服這一難題，本文提出了一種附加歷元間約束的滑動窗算法(簡稱多歷元遞推單頻PPP算法)。該方法無須額外的電離層延遲先驗信息與模型，直接將各衛星斜路徑上的電離層延遲作為參數估計，降低了電離層延遲對定位精度的影響。另外，本文采用附加約束的平方根信息濾波(square root information filter,SRIF)進行解算，可以充分利用多歷元遞推過程中各個參數的時間相關性建立內部約束，提高定位精度和收斂速度。

1 單歷元單頻精密單點定位算法

參數合并后的單頻PPP模型為

(1)

除了三參數法[13]和附加約束的單頻PPP算法[14]，還可以添加模糊度約束[15]。文獻[15]提出了同時估計電離層延遲的方法，將式(1)的未知參數分為兩類

(2)

(3)

式中，k表示第k歷元。

可將觀測方程寫為

Lk=Hy,kyk+Hg,kgk+vk,vk～N(0,Qk)

(4)

式中，Qk為對角陣；Hg,k和Hy,k是對應未知向量gk和yk的偏導數矩陣。

該方法通過參數約化減少解算時間，但定位精度仍受限于約束方程的準確性。

2 多歷元單頻定位模型及解算方法

相對于用模型改正電離層延遲，直接對電離層延遲進行估計定位精度更高。約束條件中，先驗信息的準確性及約束建立的客觀性都會導致精度下降及收斂時間的增加。多頻非差非組合PPP對電離層延遲的處理方式分為3種：①僅通過電離層模型對其進行修復[20]；②直接將電離層斜延遲作為參數進行估計[21]；③添加電離層延遲約束，并對電離層斜延遲進行參數估計[22]。

方法①與方法③已在單頻PPP中得到應用，但秩虧導致無法直接采用方法②。在多頻非差非組合PPP中，方法③通過添加電離層延遲約束，分離硬件延遲與電離層延遲參數，定位精度比方法①與方法②更高，而單頻PPP中電離層延遲參數只含有衛星端硬件延遲，可通過DCB(differential code biases)文件進行改正，不添加電離層延遲約束就達到較高的定位精度，避免外部先驗約束的影響。

對于方法②造成的秩虧問題，假設可用衛星數為n，則單歷元觀測的方程個數為2×n個，待定參數包括三維坐標、接收機鐘差、對流層天頂濕延遲、n個站星斜徑分量電離層延遲以及n個模糊度參數，總共5+2×n個參數，秩虧數為5。此時，可通過聯合多個歷元的觀測值，在解決秩虧問題的同時不引入額外的誤差。

2.1 觀測方程

假定有m個歷元，n顆衛星，線性化后的多歷元單頻PPP模型可表示為

(5)

滑動窗內的觀測方程可簡化為

Y=Ax+ε

(6)

(7)

未發生周跳時，整周模糊度數值不變，因此，在多歷元聯合解算時，每顆衛星只需解算一個模糊度參數，此時通過聯合多個歷元，即可得到唯一解。假設有n顆衛星，聯合m個歷元，則共有2×n×m個觀測值，靜態情況下，需要估計3個坐標參數、1個對流層濕延遲參數、m個鐘差參數、n×m個電離層參數及n個模糊度參數，待求參數有4+m+n×m+n個，其自由度為n×m-1)-m-4；動態情況下有4×m+1+n×m+n個參數，自由度為n×(m-1)-4×m-1。此時，除坐標參數需要估計3×m個以外，其余參數設置均與靜態情況相同。

若靜、動態情況下解唯一，則可用衛星數與聯合的歷元數需要滿足的條件為

(8)

式中，nstatic為靜態情況下的可用衛星數；nkinematic為動態情況下的可用衛星數；INT(·)為取整函數。

由式(8)可知，當聯合的歷元數增加，所需衛星數在逐漸變少?？紤]到計算效率及計算精度，本文選擇聯合5個歷元觀測值進行計算，此時靜態情況下至少需要3顆衛星，動態情況下至少需要6顆衛星。

2.2 隨機模型

利用衛星高度角和觀測值噪聲來確定偽距和載波觀測值的權，其中σP=0.1 m，σφ=0.001 m。

另外，本文采用多歷元的解算模式，可以考慮歷元間觀測值的相關性，并將其體現在方差-協方差矩陣中，此時該矩陣不再是對角陣，表示如下

(9)

式中，QYmYm為第m個歷元觀測值之間的方差-協方差陣；QY1Ym指第1個歷元的觀測值與第m個歷元觀測值之間的方差-協方差陣。

自相關系數法可以描述不同歷元觀測值之間的時間相關性[23]，本文通過自相關系數法來確定不同歷元間觀測值的方差-協方差陣。該方法定義如下

ρτ=ρij=ρjiτ=|i-j| (i，j=1,2，…,m)

(10)

因此，自相關函數中元素的自相關系數可以導出為

(11)

對于參數的隨機模型，靜態情況下，坐標固定且先驗約束為100 m，動態情況下，坐標采用白噪聲模型且先驗約束也為100 m；接收機鐘差采用白噪聲模型；對流層干延遲采用Saastamoinen模型改正，其濕延遲采用隨機游走模型估計；電離層延遲采用隨機游走模型；模糊度當作常數處理，且為了克服接收機鐘差、電離層延遲和相位模糊度之間高度相關導致的列秩虧問題(秩虧數為1)，選擇第1個滑動窗的第1組模糊度作為S基準，并將其約束到近似值，該近似值通過單頻的相位觀測值減去偽距觀測值得到[24-25]。

對于不同歷元間參數的先驗方差-協方差陣，不同測站確定的先驗方差-協方差陣不同，這里以SHAO站DOY 91的觀測數據為例，首先通過同時估計電離層單頻PPP算法100個歷元的計算結果，確定相鄰歷元間鐘差及電離層參數的相關系數，結果見表1、表2；然后根據其確定多歷元遞推單頻PPP算法的先驗方差-協方差陣。

表1中，i=1,2,…,5；cDt 1表示第1個歷元的鐘差參數，Iono 1表示第1個歷元的電離層延遲參數。由于本文方法在靜態定位解算時，只有電離層延遲與鐘差參數需要每個歷元都設置1個，因此，表1和表2中只給出了這兩個參數的相關關系。經反復試驗后發現，在動態定位中也同樣需要考慮鐘差與電離層延遲參數在歷元間的相關性。

表1 歷元間鐘差的相關系數

表2 歷元間電離層延遲的相關系數

由表1和表2可知，不同歷元間鐘差之間的相關系數基本保持在0.6以上。而反觀電離層延遲，僅僅在相鄰歷元間相關系數大于0.6，其余歷元間的相關系數都較低，Iono 1與Iono 3之間的相關系數為0.31，Iono 1分別與Iono 4和Iono 5之間的相關系數都低于0.2，可認為相關性較弱。相鄰歷元間的電離層延遲之間以及接收機鐘差之間的相關性都較大，而不相鄰歷元間的電離層延遲相關性較小，接收機鐘差相關性仍較大。根據以上分析，本文方法必須考慮不同歷元間的電離層以及鐘差的相關性，并體現在參數的方差-協方差陣中。

在確定先驗方差-協方差陣后，進行多歷元遞推單頻PPP的解算，圖1給出了不同歷元間各個鐘差及電離層延遲參數之間的相關系數。將圖1與表1、表2進行對比，可以看出多歷元單頻PPP比單歷元單頻PPP參數之間的時間相關性要強，在進行多歷元聯合解算的過程中必須要考慮。

2.3 參數估計方法

多歷元聯合解算通常使用最小二乘法，而批處理十分占內存，不滿足動態實時系統狀態估計的需求。為了充分利用之前的解算結果，本文提出一種多歷元遞推與SRIF算法相結合的方法，建立滑動窗口，滑動窗口的大小等于聯合解算的歷元數，每次解算時，將其往下滑動一個歷元，直到最后一個歷元完成解算。由于在相鄰的兩次解算結果中，部分參數是等價的，此時可以通過這個等價關系添加歷元間的參數約束，提高解算精度、增加解的穩定性。

圖1 歷元間電離層延遲與接收機鐘差的相關系數Fig.1 Correlation coefficient of ionospheric delay and receiver clock correction between epochs

將觀測方程式(6)中的參數按被約束與不被約束分組，寫成誤差方程形式

V=Ba+Cd-l

(12)

式中，a為不被約束的m個歷元的參數；d為被約束的參數；B與C分別為a與d的系數行滿秩矩陣；l為誤差向量。

約束條件方程寫為

Dd-E=0

(13)

考慮到計算效率以及計算結果的精度，預先對不同約束組合進行試驗，并根據試驗結果分別在動靜態情況下對各參數進行不同的約束，見表3。

表3 各參數約束情況

表3共聯合了5個歷元的觀測值，靜態情況下只對接收機鐘差與電離層延遲進行約束，接收機鐘差與電離層延遲的約束歷元數都為4個，分別表示對前4個歷元的接收機鐘差與電離層延遲進行約束；動態情況下只對坐標與接收機鐘差進行約束，其約束歷元數都為3個。

圖2給出了本文算法的流程，分為如下4步。

(1) 預處理，通過多普勒積分法與相位偽距組合法探測并修復周跳[26]，利用數據探測法探測并修復粗差，探測與修復鐘跳，并剔除衛星數據不好的時段。

(2) 聯合預處理后的多歷元觀測信息，建立觀測方程并將其線性化。

(3) 根據自相關系數法確定觀測值之間的方差-協方差陣，通過同時估計電離層單頻PPP算法100個歷元的計算結果確定參數的先驗方差-協方差陣。

(4) 進行參數解算，選擇第1個滑動窗的第一組模糊度作為S基準，并使用帶約束的SRIF算法進行解算。

3 試驗分析

試驗選取全球范圍內的15個IGS站，采用2019年3月19日(DOY 78)至2019年4月1日(DOY 91)共14 d，采樣間隔為30 s的觀測數據，分別進行動靜態PPP測試分析，圖3為所選的測站分布圖。以下分別將各方法的PPP解算結果與參考真值做差，獲得E、N、U 3個分量上的坐標偏差，以分析同時估計電離層延遲的單頻PPP、雙頻無電離層PPP及多歷元遞推單頻PPP算法的收斂時間和定位精度。本文中的濾波收斂定義為E、N、U各向定位偏差均優于10 cm。為確保結果的可靠性，同時檢查首次收斂時刻后續20個歷元的位置偏差，只有當連續20個歷元的偏差都在限值以內時，才認為濾波在當前歷元收斂[27]。

圖2 多歷元遞推單頻PPP算法流程Fig.2 The flow chart of multi-epoch recursive single-frequency PPP algorithm

圖3 測站分布Fig.3 Stations distribution

3.1 靜態定位

本文通過對5個歷元的觀測數據進行聯合解算得到如圖4所示的結果。

圖4為多歷元算法與單歷元算法在單天內的PDOP值及可用衛星數的變化。在多數情況下，多歷元比單歷元的可用衛星數多且衛星幾何分布更好，但由于多歷元算法對數據質量要求較高，需要剔除數據質量不好的衛星，在部分歷元，多歷元算法反而比單歷元算法的PDOP值低且可見衛星數少。另外，對未經預處理的原始數據進行試驗，發現多歷元遞推算法更易受到周跳及粗差的影響。

圖4 多歷元與單歷元算法的PDOP值及可見衛星數對比Fig.4 Comparison of PDOP value and visible satellites between multi-epoch and single-epoch algorithm

經過以上分析可知，多歷元算法比單歷元算法更好。使用傳統的間接平差對多歷元聯合算法進行數據處理，可以得到圖5中的定位結果。由圖5可知，定位結果較差，E分量最大達到了1 m，N分量在0.5 m左右波動，U分量大部分在1 m左右波動。根據2.3節的分析可知，間接平差只處理當前存儲的數據，并沒有充分利用上一次解算的結果，若當前解算的m個歷元的數據質量較差就會導致偏差很大，在U分量上更為明顯。

圖5 SHAO站多歷元聯合的單頻PPP(間接平差)Fig.5 SHAO station single-frequency PPP based on multi-epoch (indirect adjustment)

采用本文算法，設置方案(a)、方案(b)、方案(c) 3種處理方案，分別為同時估計電離層延遲的單頻PPP，多歷元遞推單頻PPP(平方根信息濾波)以及雙頻無電離層PPP，處理得到3種方案下SHAO站的靜態定位精度，如圖6所示。相比同時估計電離層的單頻PPP算法，多歷元遞推單頻PPP算法在3個分量的精度更高，特別是在高程分量，且該方法與雙頻無電離層PPP的定位精度基本相同。

對比圖6(b)與圖5可知，相比使用間接平差的多歷元聯合單頻PPP算法，采用多歷元遞推單頻PPP算法擁有更高的定位精度，水平分量的精度能夠達到厘米級，且在高程分量上有較大改善。

剔除部分存在跳變的歷元，綜合15個測站14 d的定位結果，統計了靜態定位的RMS、STD及平均收斂歷元數，見表4。方案(b)中多歷元遞推算法的E、N和U分量上的RMS分別為0.017、0.008和0.028 m，優于方案(a)。方案(b)的收斂時間大大縮短，與方案(c)雙頻無電離層PPP的收斂速度基本保持一致。

圖6 3種方案下SHAO站的靜態定位精度Fig.6 Static positioning accuracy of SHAO station under three different schemes

表4 靜態定位平均RMS、STD和收斂時間

圖7統計了SHAO站多天及所有站DOY 91的靜態定位精度。由圖7(a)可知，除U分量波動較大外，其余分量的RMS基本沒有變化。E分量的RMS在0.015 m左右浮動；N分量的RMS略小于E分量，基本在0.005 m左右；U分量的RMS比E和N分量大，大約為0.025 m，但最大也僅為0.032 m。圖7(b)中，水平分量的變化不大，E分量的RMS優于0.02 m，只有BAKO站超過了0.02 m，這是因為BAKO站在赤道附近電離層運動活躍，影響了定位精度；N分量整體優于0.015 m，且大多數測站不超過0.01 m；U分量的RMS結果相比E與N分量略差，僅有CONZ站的高程精度達到0.07 m，是由于該測站離海邊較近，海潮改正的精度較差從而影響了其高程精度，但其余測站的RMS都優于0.08 m。整體來看，N分量定位精度均優于E分量，這是由于單頻PPP只能使用浮點解，不能將模糊度固定為整數。

3.2 動態定位

采用方案(a)、方案(b)分別為同時估計電離層延遲的單頻PPP，多歷元遞推單頻PPP(平方根信息濾波)，對多站多天的數據進行模擬動態定位試驗，表5統計了平均RMS。

表5 兩種方案下動態定位平均RMS

對比表5與表4中動靜態RMS結果可知，動態情況下兩種方法精度的差別更為明顯，其原因是靜態結果統計是對各時段收斂后最后1個歷元的定位偏差計算RMS，而動態結果統計是從各時段的收斂時刻開始對偏差序列計算RMS。

圖7 多天及多站的靜態定位RMSFig.7 Static positioning RMS for multi-days and multi-stations

圖8給出了SHAO站多天定位及多站的動態定位精度。圖8(a)中各天的變化較大，3個分量的變化量都達到了分米級，特別是在高程分量，其最大的RMS為0.16 m，最小的RMS約為0.12 m，該方法比雙頻無電離層PPP算法受到動態噪聲的影響大，但比同時估計電離層的單頻PPP算法受到的影響小。

圖8(b)中，對于不同的IGS站，由于所在的緯度不同，電離層的活躍程度不同，使用多歷元遞推算法進行解算得到的平均RMS差別較大。其中CONZ站的高程精度達到了0.4 m，由于其位于海邊，在動態情況下受到的海潮改正誤差及觀測噪聲的影響更嚴重，導致其高程精度不理想。其余測站E、N和U分量上大部分分別優于0.15、0.12和0.23 m。

采用2011年10月26日7∶30—9∶30在德國莫里茨湖采集的船載觀測數據進行試驗，采樣率為1 s，并與GNSSer軟件的RTK固定解比較，圖9給出了差值。

圖8 多天及多站的動態定位RMSFig.8 Kinematic positioning RMS for multi-days and multi-stations

圖9 船載觀測數據的定位誤差Fig.9 Positioning errors of Shipborne data

由圖9可知，PPP算法的定位結果比RTK算法差，相比模擬動態試驗，海上實測動態試驗的定位偏差更大，但動態情況下誤差在E和N分量上都不超過0.5 m，在U分量上不超過 0.7 m。

表6統計了船載數據動態定位的RMS和最大偏差，可以看出水平分量上的最大值為0.416 m，高程分量上的最大值超過了0.6 m，3個分量上的RMS不超過0.2 m，驗證了本文方法在實測動態情況下的可靠性。

表6 動態定位的最大偏差和RMS

4 結 論

本文研究了一種聯合多歷元觀測值數據進行解算的單頻PPP算法，該算法不需要外部電離層延遲信息，也不需要對觀測值進行組合或求差，通過多歷元聯合求解的方式解決了添加電離層延遲參數之后的秩虧問題，主要結論如下。

(1) 使用多歷元聯合解算，不需要對電離層進行建模，且能夠在保留原始觀測信息的同時，解算電離層延遲。另外，該算法還考慮了參數之間及觀測值之間的相關關系，提高了收斂速度及定位精度和穩定性。

(2) 根據多歷元遞推觀測方程的特殊性，通過一種帶約束的SRIF算法，充分利用了相鄰歷元間參數的相關關系，提高解算精度，減少收斂時間。

(3) 采用多歷元觀測值聯合求解，大大提高了可用衛星數，靜態定位精度優于3 cm，仿動態解約為1.5 dm，相比同時估計電離層延遲的單頻PPP方法，收斂速度提高了24%，甚至與雙頻無電離層PPP算法的收斂速度基本一致，定位精度提高了30%，尤其是高程分量定位精度提高更為明顯。