洪菊,涂銳,3,張世旋,張鵬飛,王思遙,李芳馨,劉明玥
(1.中國科學(xué)院國家授時中心,西安 710600;2.中國科學(xué)院大學(xué),北京 100049;3.中國科學(xué)院精密導(dǎo)航定位與定時技術(shù)重點實驗室,西安 710600)
實際觀測中,全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)觀測值由于受到觀測信號、傳播路徑誤差、接收機等影響不可避免的存在粗差.若不對異常值進(jìn)行剔除或控制,勢必會造成參數(shù)估值的扭曲[1].因此,異常值的探測與消除是GNSS 質(zhì)量控制的重要工作.同時,對于未準(zhǔn)確探測出來的小周跳也會影響GNSS精密定位的結(jié)果.社會生產(chǎn)和生活對高精度導(dǎo)航定位的需求促進(jìn)了多系統(tǒng)GNSS 以及多種高精度、高可靠性定位技術(shù)的發(fā)展[2-6].觀測值域增強精密單點定位(PPP) 技術(shù)利用參考站提取的綜合改正數(shù)實現(xiàn)快速高精度的PPP 定位,具有收斂快、精度高等特點而受到廣泛關(guān)注.該項技術(shù)中,參考站的綜合改正數(shù)作為虛擬觀測值直接參與用戶端PPP 解算,綜合改正數(shù)中的異常對用戶端定位性能會產(chǎn)生直接影響.
當(dāng)前,對粗差的質(zhì)量控制,主要有兩種方法,第一種是將粗差歸入函數(shù)模型,將含粗差的觀測值看作與其他同類觀測值具有相同的方差,但數(shù)學(xué)期望發(fā)生改變,通過構(gòu)造檢驗統(tǒng)計量來對粗差進(jìn)行判斷和剔除[7];第二種是將粗差歸入隨機模型的方差膨脹模型,將含粗差的觀測值看作與其他同類觀測值具有相同的數(shù)學(xué)期望,但其方差發(fā)生變化,即抗差估計[8].通過等價權(quán)來對粗差觀測值進(jìn)行控制,以消除或減弱異常誤差的影響,如IGG 方案、雙因子等價權(quán)以及等價方差-協(xié)方差函數(shù)等[9-11],但抗差估計處理對于海量觀測數(shù)據(jù)會增加計算負(fù)擔(dān)[12].
基于這些前期的理論基礎(chǔ)和觀測值域增強PPP 技術(shù)中綜合改正數(shù)特性,本文在數(shù)據(jù)預(yù)處理階段,對綜合改正數(shù)進(jìn)行頻率間和二階歷元間差分,對得到的差分組合值采用中位數(shù)法探測可能存在的粗差或者周跳,并對使用異常綜合改正數(shù)的衛(wèi)星進(jìn)行模糊度重初始化、降權(quán)或剔除控制處理,從而有效控制綜合改正數(shù)異常對定位結(jié)果的影響.
單頻相位觀測值可表達(dá)為
式中:下標(biāo)b和i分別為測站和頻率;上標(biāo) s 為衛(wèi)星;Rsb為幾何距離;Lsb,i為相位觀測值;dts為衛(wèi)星鐘差;c為光速;λi為頻率i的波長;dtb為接收機鐘差;為i頻點的模糊度;為i頻點電離層延遲誤差;ZWD為對流層天頂濕延遲;db,i和dis分別為接收機和衛(wèi)星i頻點相位硬件延遲偏差;為觀測值噪聲以及多路徑效應(yīng)和未模型化誤差.對流層干延遲、相位纏繞誤差、天線相位誤差以及潮汐誤差等均使用相應(yīng)模型改正,未在公式中體現(xiàn).精密衛(wèi)星鐘差產(chǎn)品包含消電離層組合形式的衛(wèi)星偽距硬件延遲偏差,即
已知參考站坐標(biāo),引入衛(wèi)星星歷與衛(wèi)星鐘差產(chǎn)品,消除與測站相關(guān)的可模型化誤差和站星幾何距離后,以載波觀測值為例,綜合改正數(shù)可表達(dá)為
式中:orbsb和sclksb分別為衛(wèi)星軌道和衛(wèi)星鐘差殘余量;omcsb,i為觀測值減去計算值(observed minus computed,omc),稱之為觀測值域綜合改正數(shù).
流動站上的改正數(shù)可以根據(jù)其周圍參考站估計的非差omc線性插值得到,當(dāng)選擇三個參考站(b1、b2、b3)時非差omc內(nèi)插值omcsL,i,b可表達(dá)為[13]
式中,A、B、C為內(nèi)插系數(shù),計算方法在此不在贅述.
在站網(wǎng)范圍不大、大氣變化不太劇烈的情況下可以認(rèn)為內(nèi)插后的改正數(shù)與流動站真實值差距很小,因此內(nèi)插后的綜合改正數(shù)應(yīng)用在流動站且忽略觀測噪聲和多徑效應(yīng)后,單頻載波觀測值可表達(dá)為
其中:
相位電離層殘差組合二階歷元間差分受電離層殘差等各項誤差影響較小,因此常被用來探測周跳[14].同樣,不同頻率相位觀測值域改正值之差可以表示為(以第1、2 頻率為例)
由式(6)可知,該組合值消去了星歷誤差、對流層延遲等與頻率無關(guān)的誤差,只包含頻率間電離層延遲、模糊度參數(shù)、多路徑效應(yīng)和觀測噪聲、衛(wèi)星端未校準(zhǔn)相位延遲(UPD)以及接收機端UPD.硬件延遲和多徑隨時間變化比較緩慢[15],因此,對頻間差分改正數(shù)進(jìn)行歷元間二次差后得到的二階歷元間差分組合值 Δ?omcsL,GF,b在連續(xù)弧段下主要受觀測值精度和電離層殘差的影響,在零值附近波動.因此,通過對上述差分組合值序列檢驗識別綜合改正數(shù)中存在的異常.考慮到觀測噪聲和電離層殘余量會隨著歷元間隔變化而不同,采用中位數(shù)(MAD)探測方法,定義為[16]
其中,因子0.674 5 使得MAD 等同于正態(tài)分布數(shù)據(jù)的標(biāo)準(zhǔn)偏差.當(dāng)數(shù)據(jù)在 [m-n·MAD,m+n·MAD] 范圍之外時,被識別為異常值,整數(shù)n為常量.若用過于嚴(yán)格的判別標(biāo)準(zhǔn)導(dǎo)致的錯判會影響解的精度和效率,因此本文取n為10.
根據(jù)式(4)和誤差傳播定律得
根據(jù)經(jīng)驗值[17],載波測量中誤差均為mφ=±0.01周,以4倍檢測量中誤差為限差作為檢驗異常的最小閾值判斷條件,綜合中位數(shù)判別條件在數(shù)據(jù)預(yù)處理階段對存在的異常值進(jìn)行判別.
結(jié)合上述理論方法,本文采用的異常值識別、定位與控制流程如圖1所示.首先計算參考站第j歷元綜合改正數(shù)omcsL,GF,b(j),并得到二階歷元間差分組合值Δ?omcsL,GF,b(j).利用公式(7)~(9)確定異常識別閾值,若未超出閾值且上一歷元無異常,則正常解算,否則根據(jù)上一歷元組合值特性分為兩種情況控制異常值:(a)當(dāng)上一歷元 Δ?omcsL,GF,b(j-1) 為異常值時,若Δ?omcsL,GF,b(j) 與 Δ?omcsL,GF,b(j-1) 大小相近、符號相反,判斷第j-1 個歷元存在未探測出的周跳,參數(shù)估計時對該衛(wèi)星模糊度進(jìn)行重新初始化,同時從第j個歷元開始重新計算差分組合值序列,反之,采用將Δ?omcsL,GF,b(j-1)設(shè)為未存在異常的相鄰歷元差分組合值等方法反算 Δ?omcsL,GF,b(j),利用公式(7)~(9)確定異常識別閾值,若超出閾值則標(biāo)記為異常,對衛(wèi)星進(jìn)行降權(quán)或者剔除處理;(b)當(dāng)上一歷元Δ?omcsL,GF,b(j-1)無異常時,對衛(wèi)星進(jìn)行降權(quán)或剔除處理.
圖1 綜合改正數(shù)異常值識別、定位與控制流程
為分析綜合改正數(shù)異常對定位結(jié)果影響以及本文所提方法對異常值的探測情況,作者選用香港CORS 測站2021年第141 天HKST、HKSL、HKOH和HKSC 測站數(shù)據(jù),HKST、HKSL 和HKOH 作為參考站構(gòu)成參考網(wǎng),平均邊長約為26 km,HKSC 作為流動站,測站分布如圖2所示.用戶站與參考站平均距離為15 km,利用多個參考站內(nèi)插得到的用戶端大氣延遲與真實大氣延遲空間相關(guān)性強,實驗中忽略大氣延遲殘余誤差的影響.
圖2 測站分布圖
對不同采樣率下的相位omc二階歷元間差分組合值特性分析以探究異常值探測水平.圖3展示了該組數(shù)據(jù)采樣間隔為1 s、15 s、30 s 時相位omc二次歷元間差分組合值序列.1 s、15 s 和30 s 采樣間隔下各衛(wèi)星組合值始終在零值附近波動,其平均標(biāo)準(zhǔn)差分別為0.009 周、0.021 周和0.025 周.受觀測噪聲和電離層殘差等影響,1 s、15 s、30 s 采樣率下組合值的標(biāo)準(zhǔn)差有所不同,隨著采樣間隔的變長歷元間差分組合值變大.
圖3 1 s、15 s、30 s 采樣間隔下相位omc 二階歷元間差分組合值時間序列
以采樣間隔30 s 數(shù)據(jù)為例,分析不同粗差大小對仿動態(tài)增強PPP 定位結(jié)果的影響.選用單天數(shù)據(jù),每8 h 重新解算.隨機選取每個時段第960 個歷元某顆衛(wèi)星按以下方案進(jìn)行測試.方案一,人為在第一頻點L1 改正數(shù)上加入0.5 周粗差,分析較小粗差對GPS動態(tài)增強PPP 的影響情況;方案二,人為在第一頻點L1 改正數(shù)上加入2 周粗差,分析較大粗差對GPS 動態(tài)增強PPP 的影響情況;方案三,利用無異常值的綜合改正數(shù)進(jìn)行增強PPP 解算.圖4展示了參與解算的衛(wèi)星以及位置精度因子(PDOP)值,衛(wèi)星數(shù)量平均為7.4 顆,PDOP 值平均為2.04,衛(wèi)星幾何分布較好.
圖4 仿動態(tài)PPP 增強解算中PDOP 值和衛(wèi)星數(shù)量
由圖5可知,相比于平面方向,加入0.5 周小粗差對天頂(U)方向影響較大,但是同樣大小粗差對定位結(jié)果的影響不同,第1 和第3 時段對U 方向影響在1~2 cm,在第2 個時段對U 方向影響約為1 dm,這與觀測值的權(quán)有關(guān).一般認(rèn)為高度角較大時,觀測值多路徑效應(yīng)、對流層延遲及噪聲等誤差較小,GNSS 數(shù)據(jù)處理中通常采用衛(wèi)星高度角相關(guān)的隨機模型[18-20].第1~3 時段所選衛(wèi)星分別為G28、G13 和G29,加入誤差時對應(yīng)高度角分別為17°、75°和20°,所以加入相同誤差時第2 個時段對定位結(jié)果影響較大.當(dāng)加入2 周粗差后,結(jié)果如圖6所示,對東(E)和U 方向定位精度均有影響,對U 方向定位精度影響最大,第2 個時段影響最大可達(dá)0.3 m.因此選擇合適的質(zhì)量控制方法是十分必要的.
圖5 加入0.5 周較小異常時PPP 增強定位結(jié)果序列
圖6 加入2 周較大異常時PPP 增強定位結(jié)果序列
中位數(shù)法中常量取值不同對異常探測能力不同,本文放寬檢驗標(biāo)準(zhǔn),n取10.圖7展示了不同衛(wèi)星在不同時段的探測閾值,圖7(a)、(b)分別代表閾值上限和下限,若omc二階歷元間差分組合值超過閾值被判定為異常,即內(nèi)插后的綜合改正數(shù)存在的異常引起差分組合值超過閾值時可被識別.圖中不同顏色代表不同衛(wèi)星,由圖可知,不同衛(wèi)星在不同時段異常探測能力不同,當(dāng)n取10 時可探測出差分組合值中大部分1 周以上的較大粗差和部分1 周以內(nèi)的小粗差.表1結(jié)果表明利用該方法能夠有效探測出加入的0.5 周和2 周粗差.但是根據(jù)圖7所示,某些時段閾值較大,不能完全探測出0.5 周小粗差,對于未探測出的異??煽紤]調(diào)整n大小或者采用抗差估計消除或減弱異常誤差的影響.
圖7 異常探測方法閾值(不同顏色代表不同衛(wèi)星)
表1 異常值探測情況
為了進(jìn)一步檢驗上述方法的有效性,筆者選用科廷大學(xué)CORS 網(wǎng)2018年第182 天CUT0-CUT2 零基線數(shù)據(jù)進(jìn)行仿動態(tài)實驗.以CUT2 為參考站,圖8展示了CUT0 測站未進(jìn)行異常值探測時PPP 增強解算的結(jié)果.可以看出,后半段定位結(jié)果存在分米量級的異常抖動.經(jīng)分析,G21 衛(wèi)星綜合改正數(shù)存在異常,較大異常值處的歷元二次歷元間差分組合值為-0.284 周,下一歷元其值為0.289 周,因此判斷異常處發(fā)生周跳.其對當(dāng)前歷元以及后面所有歷元定位結(jié)果造成高達(dá)3 dm 的誤差.因此,本文對該異常值處衛(wèi)星降權(quán)處理并對下一歷元模糊度重新初始化.圖9展示的是未進(jìn)行異常值控制與進(jìn)行異常值控制后PPP 增強定位結(jié)果的對比,由圖9可知,對綜合改正數(shù)進(jìn)行質(zhì)量控制后,PPP 增強定位結(jié)果精度有了明顯提升.
圖8 未進(jìn)行異常值識別與控制時PPP 增強定位結(jié)果序列
圖9 異常值控制前后PPP 增強定位結(jié)果比對
綜合改正數(shù)質(zhì)量直接影響PPP 增強定位精度,因此,對綜合改正數(shù)中異常值的探測與控制在觀測值域增強PPP 定位中具有重要作用.本文對觀測值域增強改正數(shù)質(zhì)量控制方法展開研究,針對綜合改正數(shù)特性,提出利用經(jīng)頻間和二階歷元間差分后的差分組合值,作為檢驗量進(jìn)行異常值識別,并對使用該異常值的衛(wèi)星采用模糊度重新初始化、降權(quán)或剔除方法進(jìn)行處理.差分組合值受觀測值精度的影響,因此可以通過計算組合觀測值精度探測綜合改正數(shù)中可能存在的異常,但隨著歷元間隔變化觀測噪聲和電離層殘差會變化同時避免過大的異常觀測值對檢驗所產(chǎn)生的破壞性影響,選用中位數(shù)法設(shè)置探測閾值.中位數(shù)法常量設(shè)置決定了探測水平,為適當(dāng)放寬檢驗標(biāo)準(zhǔn),本文常量取值為10.經(jīng)實驗驗證該方法能夠有效探測出差分組合值中大部分1 周以上的較大異常和部分1 周以內(nèi)的異常,并針對實驗中出現(xiàn)的異常通過判斷異常類型采用模糊度重新初始化方法抑制了異常對定位的影響.在GNSS 定位解算中,通常采用高度角定權(quán)模型,同時經(jīng)過實驗分析可知,不同高度角衛(wèi)星的綜合改正數(shù)存在的異常對定位結(jié)果影響不同.因此,針對不同高度角衛(wèi)星存在的異常值處理方法仍需進(jìn)一步探究.上述方法存在一定的局限性,比如無法識別不同頻點改正數(shù)存在的相同異常,下一步繼續(xù)研究結(jié)合MW 組合用于異常識別的方法.
致謝:感謝科廷大學(xué)GNSS-SPAN 和香港特別行政區(qū)政府地政總署測繪處開源的CORS 數(shù)據(jù),本文的研究得到了國家自然科學(xué)基金項目(41974032)的支持.