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

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

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

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

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

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

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

參數(shù)合并后的單頻PPP模型為

(1)

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

(2)

(3)

式中，k表示第k歷元。

可將觀測方程寫為

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

(4)

式中，Qk為對角陣；Hg,k和Hy,k是對應(yīng)未知向量gk和yk的偏導(dǎo)數(shù)矩陣。

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

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

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

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

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

2.1 觀測方程

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

(5)

滑動窗內(nèi)的觀測方程可簡化為

Y=Ax+ε

(6)

(7)

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

若靜、動態(tài)情況下解唯一，則可用衛(wèi)星數(shù)與聯(lián)合的歷元數(shù)需要滿足的條件為

(8)

式中，nstatic為靜態(tài)情況下的可用衛(wèi)星數(shù)；nkinematic為動態(tài)情況下的可用衛(wèi)星數(shù)；INT(·)為取整函數(shù)。

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

2.2 隨機模型

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

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

(9)

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

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

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

(10)

因此，自相關(guān)函數(shù)中元素的自相關(guān)系數(shù)可以導(dǎo)出為

(11)

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

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

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

表1 歷元間鐘差的相關(guān)系數(shù)

表2 歷元間電離層延遲的相關(guān)系數(shù)

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

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

2.3 參數(shù)估計方法

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

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

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

V=Ba+Cd-l

(12)

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

約束條件方程寫為

Dd-E=0

(13)

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

表3 各參數(shù)約束情況

表3共聯(lián)合了5個歷元的觀測值，靜態(tài)情況下只對接收機鐘差與電離層延遲進行約束，接收機鐘差與電離層延遲的約束歷元數(shù)都為4個，分別表示對前4個歷元的接收機鐘差與電離層延遲進行約束；動態(tài)情況下只對坐標(biāo)與接收機鐘差進行約束，其約束歷元數(shù)都為3個。

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

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

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

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

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

3 試驗分析

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

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

圖3 測站分布Fig.3 Stations distribution

3.1 靜態(tài)定位

本文通過對5個歷元的觀測數(shù)據(jù)進行聯(lián)合解算得到如圖4所示的結(jié)果。

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

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

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

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

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

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

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

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

表4 靜態(tài)定位平均RMS、STD和收斂時間

圖7統(tǒng)計了SHAO站多天及所有站DOY 91的靜態(tài)定位精度。由圖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優(yōu)于0.02 m，只有BAKO站超過了0.02 m，這是因為BAKO站在赤道附近電離層運動活躍，影響了定位精度；N分量整體優(yōu)于0.015 m，且大多數(shù)測站不超過0.01 m；U分量的RMS結(jié)果相比E與N分量略差，僅有CONZ站的高程精度達到0.07 m，是由于該測站離海邊較近，海潮改正的精度較差從而影響了其高程精度，但其余測站的RMS都優(yōu)于0.08 m。整體來看，N分量定位精度均優(yōu)于E分量，這是由于單頻PPP只能使用浮點解，不能將模糊度固定為整數(shù)。

3.2 動態(tài)定位

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

表5 兩種方案下動態(tài)定位平均RMS

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

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

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

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

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

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

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

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

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

表6 動態(tài)定位的最大偏差和RMS

4 結(jié) 論

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

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

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

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