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

(1.中國(guó)科學(xué)院國(guó)家授時(shí)中心,西安 710600；2.中國(guó)科學(xué)院大學(xué),北京 100049；3.中國(guó)科學(xué)院精密導(dǎo)航定位與定時(shí)技術(shù)重點(diǎn)實(shí)驗(yàn)室,西安 710600)

0 引 言

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

當(dāng)前，對(duì)粗差的質(zhì)量控制，主要有兩種方法，第一種是將粗差歸入函數(shù)模型，將含粗差的觀測(cè)值看作與其他同類觀測(cè)值具有相同的方差，但數(shù)學(xué)期望發(fā)生改變，通過構(gòu)造檢驗(yàn)統(tǒng)計(jì)量來對(duì)粗差進(jìn)行判斷和剔除[7]；第二種是將粗差歸入隨機(jī)模型的方差膨脹模型，將含粗差的觀測(cè)值看作與其他同類觀測(cè)值具有相同的數(shù)學(xué)期望，但其方差發(fā)生變化，即抗差估計(jì)[8].通過等價(jià)權(quán)來對(duì)粗差觀測(cè)值進(jìn)行控制，以消除或減弱異常誤差的影響，如IGG 方案、雙因子等價(jià)權(quán)以及等價(jià)方差-協(xié)方差函數(shù)等[9-11]，但抗差估計(jì)處理對(duì)于海量觀測(cè)數(shù)據(jù)會(huì)增加計(jì)算負(fù)擔(dān)[12].

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

1 觀測(cè)值域增強(qiáng)PPP 和綜合改正數(shù)質(zhì)量控制方法

1.1 觀測(cè)值域增強(qiáng)PPP 算法

單頻相位觀測(cè)值可表達(dá)為

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

已知參考站坐標(biāo)，引入衛(wèi)星星歷與衛(wèi)星鐘差產(chǎn)品，消除與測(cè)站相關(guān)的可模型化誤差和站星幾何距離后，以載波觀測(cè)值為例，綜合改正數(shù)可表達(dá)為

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

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

式中，A、B、C為內(nèi)插系數(shù)，計(jì)算方法在此不在贅述.

在站網(wǎng)范圍不大、大氣變化不太劇烈的情況下可以認(rèn)為內(nèi)插后的改正數(shù)與流動(dòng)站真實(shí)值差距很小，因此內(nèi)插后的綜合改正數(shù)應(yīng)用在流動(dòng)站且忽略觀測(cè)噪聲和多徑效應(yīng)后，單頻載波觀測(cè)值可表達(dá)為

其中：

1.2 觀測(cè)值域綜合改正數(shù)質(zhì)量控制方法

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

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

其中，因子0.674 5 使得MAD 等同于正態(tài)分布數(shù)據(jù)的標(biāo)準(zhǔn)偏差.當(dāng)數(shù)據(jù)在 [m-n·MAD,m+n·MAD] 范圍之外時(shí)，被識(shí)別為異常值，整數(shù)n為常量.若用過于嚴(yán)格的判別標(biāo)準(zhǔn)導(dǎo)致的錯(cuò)判會(huì)影響解的精度和效率，因此本文取n為10.

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

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

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

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

2 試驗(yàn)分析

為分析綜合改正數(shù)異常對(duì)定位結(jié)果影響以及本文所提方法對(duì)異常值的探測(cè)情況，作者選用香港CORS 測(cè)站2021年第141 天HKST、HKSL、HKOH和HKSC 測(cè)站數(shù)據(jù)，HKST、HKSL 和HKOH 作為參考站構(gòu)成參考網(wǎng)，平均邊長(zhǎng)約為26 km，HKSC 作為流動(dòng)站，測(cè)站分布如圖2所示.用戶站與參考站平均距離為15 km，利用多個(gè)參考站內(nèi)插得到的用戶端大氣延遲與真實(shí)大氣延遲空間相關(guān)性強(qiáng)，實(shí)驗(yàn)中忽略大氣延遲殘余誤差的影響.

圖2 測(cè)站分布圖

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

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

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

圖4 仿動(dòng)態(tài)PPP 增強(qiáng)解算中PDOP 值和衛(wèi)星數(shù)量

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

圖5 加入0.5 周較小異常時(shí)PPP 增強(qiáng)定位結(jié)果序列

圖6 加入2 周較大異常時(shí)PPP 增強(qiáng)定位結(jié)果序列

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

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

表1 異常值探測(cè)情況

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

圖8 未進(jìn)行異常值識(shí)別與控制時(shí)PPP 增強(qiáng)定位結(jié)果序列

圖9 異常值控制前后PPP 增強(qiáng)定位結(jié)果比對(duì)

3 結(jié)論及展望

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

致謝：感謝科廷大學(xué)GNSS-SPAN 和香港特別行政區(qū)政府地政總署測(cè)繪處開源的CORS 數(shù)據(jù)，本文的研究得到了國(guó)家自然科學(xué)基金項(xiàng)目(41974032)的支持.