GNSS觀測值域增強PPP定位中綜合改正數質量控制算法

洪菊，涂銳,3，張世旋，張鵬飛，王思遙，李芳馨，劉明玥

(1.中國科學院國家授時中心,西安 710600；2.中國科學院大學,北京 100049；3.中國科學院精密導航定位與定時技術重點實驗室,西安 710600)

0 引 言

實際觀測中，全球衛星導航系統(GNSS)觀測值由于受到觀測信號、傳播路徑誤差、接收機等影響不可避免的存在粗差.若不對異常值進行剔除或控制，勢必會造成參數估值的扭曲[1].因此，異常值的探測與消除是GNSS 質量控制的重要工作.同時，對于未準確探測出來的小周跳也會影響GNSS精密定位的結果.社會生產和生活對高精度導航定位的需求促進了多系統GNSS 以及多種高精度、高可靠性定位技術的發展[2-6].觀測值域增強精密單點定位(PPP) 技術利用參考站提取的綜合改正數實現快速高精度的PPP 定位，具有收斂快、精度高等特點而受到廣泛關注.該項技術中，參考站的綜合改正數作為虛擬觀測值直接參與用戶端PPP 解算，綜合改正數中的異常對用戶端定位性能會產生直接影響.

當前，對粗差的質量控制，主要有兩種方法，第一種是將粗差歸入函數模型，將含粗差的觀測值看作與其他同類觀測值具有相同的方差，但數學期望發生改變，通過構造檢驗統計量來對粗差進行判斷和剔除[7]；第二種是將粗差歸入隨機模型的方差膨脹模型，將含粗差的觀測值看作與其他同類觀測值具有相同的數學期望，但其方差發生變化，即抗差估計[8].通過等價權來對粗差觀測值進行控制，以消除或減弱異常誤差的影響，如IGG 方案、雙因子等價權以及等價方差-協方差函數等[9-11]，但抗差估計處理對于海量觀測數據會增加計算負擔[12].

基于這些前期的理論基礎和觀測值域增強PPP 技術中綜合改正數特性，本文在數據預處理階段，對綜合改正數進行頻率間和二階歷元間差分，對得到的差分組合值采用中位數法探測可能存在的粗差或者周跳，并對使用異常綜合改正數的衛星進行模糊度重初始化、降權或剔除控制處理，從而有效控制綜合改正數異常對定位結果的影響.

1 觀測值域增強PPP 和綜合改正數質量控制方法

1.1 觀測值域增強PPP 算法

單頻相位觀測值可表達為

式中：下標b和i分別為測站和頻率；上標 s 為衛星；Rsb為幾何距離；Lsb,i為相位觀測值；dts為衛星鐘差；c為光速；λi為頻率i的波長；dtb為接收機鐘差；為i頻點的模糊度；為i頻點電離層延遲誤差；ZWD為對流層天頂濕延遲；db,i和dis分別為接收機和衛星i頻點相位硬件延遲偏差；為觀測值噪聲以及多路徑效應和未模型化誤差.對流層干延遲、相位纏繞誤差、天線相位誤差以及潮汐誤差等均使用相應模型改正，未在公式中體現.精密衛星鐘差產品包含消電離層組合形式的衛星偽距硬件延遲偏差，即

已知參考站坐標，引入衛星星歷與衛星鐘差產品，消除與測站相關的可模型化誤差和站星幾何距離后，以載波觀測值為例，綜合改正數可表達為

式中：orbsb和sclksb分別為衛星軌道和衛星鐘差殘余量；omcsb,i為觀測值減去計算值(observed minus computed，omc)，稱之為觀測值域綜合改正數.

流動站上的改正數可以根據其周圍參考站估計的非差omc線性插值得到，當選擇三個參考站(b1、b2、b3)時非差omc內插值omcsL,i,b可表達為[13]

式中，A、B、C為內插系數，計算方法在此不在贅述.

在站網范圍不大、大氣變化不太劇烈的情況下可以認為內插后的改正數與流動站真實值差距很小，因此內插后的綜合改正數應用在流動站且忽略觀測噪聲和多徑效應后，單頻載波觀測值可表達為

其中：

1.2 觀測值域綜合改正數質量控制方法

相位電離層殘差組合二階歷元間差分受電離層殘差等各項誤差影響較小，因此常被用來探測周跳[14].同樣，不同頻率相位觀測值域改正值之差可以表示為(以第1、2 頻率為例)

由式(6)可知，該組合值消去了星歷誤差、對流層延遲等與頻率無關的誤差，只包含頻率間電離層延遲、模糊度參數、多路徑效應和觀測噪聲、衛星端未校準相位延遲(UPD)以及接收機端UPD.硬件延遲和多徑隨時間變化比較緩慢[15]，因此，對頻間差分改正數進行歷元間二次差后得到的二階歷元間差分組合值 Δ?omcsL,GF,b在連續弧段下主要受觀測值精度和電離層殘差的影響，在零值附近波動.因此，通過對上述差分組合值序列檢驗識別綜合改正數中存在的異常.考慮到觀測噪聲和電離層殘余量會隨著歷元間隔變化而不同，采用中位數(MAD)探測方法，定義為[16]

其中，因子0.674 5 使得MAD 等同于正態分布數據的標準偏差.當數據在 [m-n·MAD,m+n·MAD] 范圍之外時，被識別為異常值，整數n為常量.若用過于嚴格的判別標準導致的錯判會影響解的精度和效率，因此本文取n為10.

根據式(4)和誤差傳播定律得

根據經驗值[17]，載波測量中誤差均為mφ=±0.01周，以4倍檢測量中誤差為限差作為檢驗異常的最小閾值判斷條件，綜合中位數判別條件在數據預處理階段對存在的異常值進行判別.

結合上述理論方法，本文采用的異常值識別、定位與控制流程如圖1所示.首先計算參考站第j歷元綜合改正數omcsL,GF,b(j)，并得到二階歷元間差分組合值Δ?omcsL,GF,b(j).利用公式(7)～(9)確定異常識別閾值，若未超出閾值且上一歷元無異常，則正常解算，否則根據上一歷元組合值特性分為兩種情況控制異常值：(a)當上一歷元 Δ?omcsL,GF,b(j-1) 為異常值時，若Δ?omcsL,GF,b(j) 與 Δ?omcsL,GF,b(j-1) 大小相近、符號相反，判斷第j-1 個歷元存在未探測出的周跳，參數估計時對該衛星模糊度進行重新初始化，同時從第j個歷元開始重新計算差分組合值序列，反之，采用將Δ?omcsL,GF,b(j-1)設為未存在異常的相鄰歷元差分組合值等方法反算 Δ?omcsL,GF,b(j)，利用公式(7)～(9)確定異常識別閾值，若超出閾值則標記為異常，對衛星進行降權或者剔除處理；(b)當上一歷元Δ?omcsL,GF,b(j-1)無異常時，對衛星進行降權或剔除處理.

圖1 綜合改正數異常值識別、定位與控制流程

2 試驗分析

為分析綜合改正數異常對定位結果影響以及本文所提方法對異常值的探測情況，作者選用香港CORS 測站2021年第141 天HKST、HKSL、HKOH和HKSC 測站數據，HKST、HKSL 和HKOH 作為參考站構成參考網，平均邊長約為26 km，HKSC 作為流動站，測站分布如圖2所示.用戶站與參考站平均距離為15 km，利用多個參考站內插得到的用戶端大氣延遲與真實大氣延遲空間相關性強，實驗中忽略大氣延遲殘余誤差的影響.

圖2 測站分布圖

對不同采樣率下的相位omc二階歷元間差分組合值特性分析以探究異常值探測水平.圖3展示了該組數據采樣間隔為1 s、15 s、30 s 時相位omc二次歷元間差分組合值序列.1 s、15 s 和30 s 采樣間隔下各衛星組合值始終在零值附近波動，其平均標準差分別為0.009 周、0.021 周和0.025 周.受觀測噪聲和電離層殘差等影響，1 s、15 s、30 s 采樣率下組合值的標準差有所不同，隨著采樣間隔的變長歷元間差分組合值變大.

圖3 1 s、15 s、30 s 采樣間隔下相位omc 二階歷元間差分組合值時間序列

以采樣間隔30 s 數據為例，分析不同粗差大小對仿動態增強PPP 定位結果的影響.選用單天數據，每8 h 重新解算.隨機選取每個時段第960 個歷元某顆衛星按以下方案進行測試.方案一，人為在第一頻點L1 改正數上加入0.5 周粗差，分析較小粗差對GPS動態增強PPP 的影響情況；方案二，人為在第一頻點L1 改正數上加入2 周粗差，分析較大粗差對GPS 動態增強PPP 的影響情況；方案三，利用無異常值的綜合改正數進行增強PPP 解算.圖4展示了參與解算的衛星以及位置精度因子(PDOP)值，衛星數量平均為7.4 顆，PDOP 值平均為2.04，衛星幾何分布較好.

圖4 仿動態PPP 增強解算中PDOP 值和衛星數量

由圖5可知，相比于平面方向，加入0.5 周小粗差對天頂(U)方向影響較大，但是同樣大小粗差對定位結果的影響不同，第1 和第3 時段對U 方向影響在1～2 cm，在第2 個時段對U 方向影響約為1 dm，這與觀測值的權有關.一般認為高度角較大時，觀測值多路徑效應、對流層延遲及噪聲等誤差較小，GNSS 數據處理中通常采用衛星高度角相關的隨機模型[18-20].第1～3 時段所選衛星分別為G28、G13 和G29，加入誤差時對應高度角分別為17°、75°和20°，所以加入相同誤差時第2 個時段對定位結果影響較大.當加入2 周粗差后，結果如圖6所示，對東(E)和U 方向定位精度均有影響，對U 方向定位精度影響最大，第2 個時段影響最大可達0.3 m.因此選擇合適的質量控制方法是十分必要的.

圖5 加入0.5 周較小異常時PPP 增強定位結果序列

圖6 加入2 周較大異常時PPP 增強定位結果序列

中位數法中常量取值不同對異常探測能力不同，本文放寬檢驗標準，n取10.圖7展示了不同衛星在不同時段的探測閾值，圖7(a)、(b)分別代表閾值上限和下限，若omc二階歷元間差分組合值超過閾值被判定為異常，即內插后的綜合改正數存在的異常引起差分組合值超過閾值時可被識別.圖中不同顏色代表不同衛星，由圖可知，不同衛星在不同時段異常探測能力不同，當n取10 時可探測出差分組合值中大部分1 周以上的較大粗差和部分1 周以內的小粗差.表1結果表明利用該方法能夠有效探測出加入的0.5 周和2 周粗差.但是根據圖7所示，某些時段閾值較大，不能完全探測出0.5 周小粗差，對于未探測出的異常可考慮調整n大小或者采用抗差估計消除或減弱異常誤差的影響.

圖7 異常探測方法閾值(不同顏色代表不同衛星)

表1 異常值探測情況

為了進一步檢驗上述方法的有效性，筆者選用科廷大學CORS 網2018年第182 天CUT0-CUT2 零基線數據進行仿動態實驗.以CUT2 為參考站，圖8展示了CUT0 測站未進行異常值探測時PPP 增強解算的結果.可以看出，后半段定位結果存在分米量級的異常抖動.經分析，G21 衛星綜合改正數存在異常，較大異常值處的歷元二次歷元間差分組合值為-0.284 周，下一歷元其值為0.289 周，因此判斷異常處發生周跳.其對當前歷元以及后面所有歷元定位結果造成高達3 dm 的誤差.因此，本文對該異常值處衛星降權處理并對下一歷元模糊度重新初始化.圖9展示的是未進行異常值控制與進行異常值控制后PPP 增強定位結果的對比，由圖9可知，對綜合改正數進行質量控制后，PPP 增強定位結果精度有了明顯提升.

圖8 未進行異常值識別與控制時PPP 增強定位結果序列

圖9 異常值控制前后PPP 增強定位結果比對

3 結論及展望

綜合改正數質量直接影響PPP 增強定位精度，因此，對綜合改正數中異常值的探測與控制在觀測值域增強PPP 定位中具有重要作用.本文對觀測值域增強改正數質量控制方法展開研究，針對綜合改正數特性，提出利用經頻間和二階歷元間差分后的差分組合值，作為檢驗量進行異常值識別，并對使用該異常值的衛星采用模糊度重新初始化、降權或剔除方法進行處理.差分組合值受觀測值精度的影響，因此可以通過計算組合觀測值精度探測綜合改正數中可能存在的異常，但隨著歷元間隔變化觀測噪聲和電離層殘差會變化同時避免過大的異常觀測值對檢驗所產生的破壞性影響，選用中位數法設置探測閾值.中位數法常量設置決定了探測水平，為適當放寬檢驗標準，本文常量取值為10.經實驗驗證該方法能夠有效探測出差分組合值中大部分1 周以上的較大異常和部分1 周以內的異常，并針對實驗中出現的異常通過判斷異常類型采用模糊度重新初始化方法抑制了異常對定位的影響.在GNSS 定位解算中，通常采用高度角定權模型，同時經過實驗分析可知，不同高度角衛星的綜合改正數存在的異常對定位結果影響不同.因此，針對不同高度角衛星存在的異常值處理方法仍需進一步探究.上述方法存在一定的局限性，比如無法識別不同頻點改正數存在的相同異常，下一步繼續研究結合MW 組合用于異常識別的方法.

致謝：感謝科廷大學GNSS-SPAN 和香港特別行政區政府地政總署測繪處開源的CORS 數據，本文的研究得到了國家自然科學基金項目(41974032)的支持.