面向復雜條件的北斗/GNSS數(shù)據(jù)粗差處理方法

李雪珍 章浙濤 劉歡 袁海軍

0 引言

全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)包括全球的、區(qū)域的衛(wèi)星導航系統(tǒng)以及增強系統(tǒng)，如美國的全球定位系統(tǒng)(GPS)、中國的北斗衛(wèi)星導航系統(tǒng)(BDS)、俄羅斯的全球衛(wèi)星導航系統(tǒng)(GLONASS)、歐洲的伽利略衛(wèi)星定位系統(tǒng)(Galileo)、日本的準天頂系統(tǒng)(QZSS)等.GNSS在全球或區(qū)域范圍內(nèi)提供定位、導航和授時(PNT)服務(wù)，被廣泛應(yīng)用于自然災(zāi)害監(jiān)測、變形監(jiān)測、航空運輸、防震減災(zāi)、車輛駕駛等領(lǐng)域.BDS是中國自主研發(fā)、獨立運行的衛(wèi)星導航系統(tǒng)，其發(fā)展分為3個階段：北斗衛(wèi)星導航試驗系統(tǒng)(BDS-1)、區(qū)域北斗系統(tǒng)(BDS-2)和全球北斗系統(tǒng)(BDS-3)[1-2].如今BDS-3已經(jīng)建成，由3顆地球靜止軌道(GEO)衛(wèi)星、3顆傾斜地球同步軌道(IGSO)衛(wèi)星和24顆中圓地球軌道(MEO)衛(wèi)星組成，于2020年向全球用戶提供導航定位服務(wù)[3].BDS-3在繼承BDS-2原有B1I和B3I信號的同時，新增了B1C、B2a和B2b頻點[4].目前BDS能接收更多的可見衛(wèi)星信號，得到更優(yōu)的有效觀測數(shù)據(jù)，實現(xiàn)良好的衛(wèi)星幾何構(gòu)型，進而提升定位精度.BDS全球服務(wù)水平和高程定位精度達到10 m，在亞太地區(qū)可達到5 m，系統(tǒng)服務(wù)可用性優(yōu)于95%[5].

GNSS實際應(yīng)用環(huán)境是復雜且惡劣的.在城市峽谷環(huán)境中，高大建筑物的表面及街道上的樹木和車輛等物體使衛(wèi)星信號產(chǎn)生一定的反射和折射現(xiàn)象[6].衛(wèi)星信號的衰減、中斷導致可見衛(wèi)星數(shù)量的減少和衛(wèi)星幾何分布不理想，出現(xiàn)較多較為嚴重的粗差.少量的粗差，能對參數(shù)估計的結(jié)果造成干擾，進而嚴重影響GNSS定位的精度和可靠性.目前用戶對于定位的實時性、精度、可靠性的需求急速增加.在多頻多系統(tǒng)和城市峽谷環(huán)境等復雜條件下，如何有效抵御粗差的影響是GNSS數(shù)據(jù)處理中的棘手問題，更是數(shù)據(jù)質(zhì)量控制的關(guān)鍵點.為剔除粗差，提高定位精度，北斗/GNSS數(shù)據(jù)的粗差處理主要可概括為兩大類方法：基于均值漂移的粗差處理方法[7-8]和基于方差膨脹的粗差處理方法[9].

本文以最常用的數(shù)據(jù)探測及穩(wěn)健估計為例闡述兩類粗差處理方法的流程，并研究給出了適用于復雜條件的閾值方案，最后利用實測數(shù)據(jù)的GPS+BDS組合定位進行驗證.

1 基于均值漂移的粗差處理

GNSS函數(shù)模型中的粗差通常會導致最小二乘估計存在偏差.基于均值漂移的粗差處理即粗差納入函數(shù)模型，認為含有粗差的觀測值，其期望發(fā)生平移[10].最早提出且目前使用最為廣泛的基于均值漂移的粗差處理方法是數(shù)據(jù)探測(data snooping)[7].

具體而言，數(shù)據(jù)探測依賴于在原假設(shè)模型和一組備選假設(shè)模型之間進行假設(shè)檢驗，要求檢驗統(tǒng)計量必須遵循正態(tài)分布、τ分布、x2分布和F分布[11-12].此外，由于觀測模型的幾何形狀[11]、假設(shè)之間的可分離性[13]、選擇的統(tǒng)計量[14]、預(yù)先確定的臨界值[15]等原因，不可避免地會產(chǎn)生漏檢、誤報和錯誤識別，因此并不能完全消除偏差并確保最終參數(shù)估計的無偏性.數(shù)據(jù)探測最初只是針對一個粗差研究的.在復雜條件下,一組觀測值包含多個粗差時，數(shù)據(jù)探測只能一個接一個地探測，難以滿足實際的應(yīng)用需求.而多維粗差的同時定位與定值(ELGE法)不僅能確定多個粗差的位置，還能求出各個粗差的數(shù)值大小[16].部分最小二乘平差把觀測值是否含有粗差分成兩組，在不含粗差的一組里實施最小二乘平差，同樣能求出多個粗差的位置與數(shù)值[17].與上述以殘差為對象的方法不同的是擬準檢定法，它依據(jù)觀測值的真誤差判斷粗差的位置[18].在觀測值相互獨立且等權(quán)時，以上幾種方法可以認為是等價的[19-20].

Baarda[7]最早提出測量系統(tǒng)的可靠性理論.可靠性分析為粗差探測提供了相應(yīng)的理論基礎(chǔ).可靠性分析的3項重要指標分別是最小可檢測偏差(MDB)、內(nèi)部可靠性和外部可靠性[21].其中MDB指在一定正確檢驗概率條件下可被檢測到的最小偏差的絕對值.內(nèi)部可靠性指的是發(fā)現(xiàn)粗差的能力，外部可靠性是指不可發(fā)現(xiàn)粗差對平差結(jié)果的影響程度[22].對于多個備選假設(shè)的情況，最小可分離偏差(MSB)和最小可識別偏差(MIB)被研究和定義，它們可衡量與另一個備選假設(shè)成功分離或在多個備選假設(shè)中以一定的概率正確識別的最小偏差的大小[13,23].將MDB代入?yún)?shù)解得到外部可靠性，以衡量未檢測到的粗差對估計結(jié)果的影響顯著性.其他相關(guān)的可靠性測量也被研究，如可靠性數(shù)據(jù)指標[24]和可控性(區(qū)域三角)[25].目前，基于均值漂移的粗差處理方法已被廣泛應(yīng)用于各種領(lǐng)域，如大地測量網(wǎng)的質(zhì)量控制和完好性監(jiān)測等.

2 基于方差膨脹的粗差處理

另一類基于方差膨脹的粗差處理方法認為含有粗差的觀測值的方差發(fā)生變化，但期望不變[5].基于方差膨脹的粗差處理方法的研究主要集中在穩(wěn)健估計(抗差估計)[22,26].穩(wěn)健估計將粗差納入隨機模型，通過逐次迭代平差結(jié)果，不斷改變觀測值的權(quán)，最終使含有粗差的觀測值權(quán)為零.穩(wěn)健估計主要依賴等價權(quán)函數(shù)來檢測和排除GNSS觀測數(shù)據(jù)中的粗差[27].穩(wěn)健估計主要包括M估計(廣義極大似然估計)[28]、R估計(秩檢驗估計)[29]、L1范數(shù)估計[30]、最小截平方和估計[31]、符號約束的穩(wěn)健估計[23].M估計的抗差性和有效性主要取決于所采用的參數(shù)的可靠性、等價權(quán)函數(shù)及其臨界值的合理性[32]，因此不同權(quán)函數(shù)消除或減弱粗差影響的能力不盡相同.已有許多不同的權(quán)函數(shù)能高效和準確處理獨立觀測值和相關(guān)觀測值，如雙重法[33]、Huber法[27]、Hampel法[34]、Danish法[35]、IGGI法[36]和IGGIII法[9,37].其中IGGI作為國內(nèi)早期經(jīng)典的質(zhì)量控制算法，它依據(jù)殘差的假設(shè)檢驗結(jié)果調(diào)整權(quán)陣來減弱含粗差的觀測值對定位的影響.在IGGI的基礎(chǔ)上，利用等價權(quán)原理，進一步發(fā)展到相關(guān)觀測值的IGGIII方法.同樣基于等價權(quán)，Yang等[38]利用方差分量的穩(wěn)健估計來減弱粗差的影響.早期，在隨機模型中的觀測值的方差被認為是相等的[39].后來采用異方差假設(shè)來代替GNSS觀測的不切實際的同方差.通常采用衛(wèi)星高度角和載噪比作為觀測值異方差假設(shè)的兩個主要指標.其中高度角和觀測值精度之間的關(guān)系可以用余弦函數(shù)或指數(shù)函數(shù)來充分描述[40].由于載噪比與GNSS觀測值記錄在同一個跟蹤環(huán)路中，因此載噪比與GNSS觀測值精度高度一致[41].

此外，雙因子方差膨脹模型和雙因子等價權(quán)模型同樣有效地控制粗差對參數(shù)估值的影響[42].可從污染誤差模型入手，研究最小均方差準則的參數(shù)穩(wěn)健估計[43].楊元喜[44]根據(jù)多種抗差濾波及其性質(zhì)構(gòu)建抗差自適應(yīng)濾波理論體系.同樣地，基于方差膨脹的粗差處理方法經(jīng)常被用于大地測量數(shù)據(jù)處理，例如變形分析[29]、基準變換[45]、GNSS完好性監(jiān)測[46]、卡爾曼濾波[47].但是它們在模型結(jié)構(gòu)不強、粗差較多的情況下，并不能完全消除粗差的影響.因此在多系統(tǒng)多頻和城市峽谷環(huán)境等復雜條件下，對粗差處理方法的研究顯得尤為重要.

3 兩種粗差處理方法

3.1 數(shù)據(jù)探測

在復雜條件下，GNSS觀測值通常不可避免地包含粗差.采用數(shù)據(jù)探測的方法對粗差進行后驗處理，具體分為以下3個步驟[7-8]：

1)探測

采用整體檢驗，構(gòu)造Tq統(tǒng)計量：

(1)

2)識別

采用標準化殘差構(gòu)造wi檢驗統(tǒng)計量:

(2)

3)調(diào)節(jié)

最后利用剩余的觀測值重新進行最小二乘估計得到可靠的結(jié)果，并將其輸出.

3.2 穩(wěn)健估計

選權(quán)迭代法作為M估計的重要方法，同時是穩(wěn)健估計中最為常用且計算簡便的方法，其基本步驟如下：

1)建立數(shù)學模型

(3)

2)按最小二乘法求解參數(shù)估值及其殘差

(4)

(5)

3)求解觀測值的等價權(quán)矩陣，應(yīng)用抗差最小二乘迭代計算，即

(6)

選定某一值ω作為迭代閾值(后續(xù)實驗所討論的閾值)，當?shù)趉次與第k-1次迭代所得的估值之差的絕對值小于等于該值，則停止迭代.

(7)

4)最后求解

(8)

(9)

選權(quán)迭代法的關(guān)鍵在于選擇權(quán)函數(shù)，最經(jīng)典的是IGGⅢ方法.IGGⅢ方法的權(quán)函數(shù)pi確定如下:

(10)

(11)

(12)

式中，q為多余觀測數(shù).

本文中討論了面向復雜條件時ω、k0、k1的選取，由于過多地剔除觀測值，會導致定位精度降低，因此可以通過迭代閾值(ω)來控制迭代次數(shù)，使穩(wěn)健估計既具抗差性，又高效.

4 實際數(shù)據(jù)分析

為驗證本文提出的粗差處理方案的有效性并評價定位精度，本文選取2021年年積日第288天24小時的5組西南地區(qū)滑坡監(jiān)測數(shù)據(jù)進行實驗.5組數(shù)據(jù)分別命名為Test1、Test2、Test3、Test4和Test5，其基線長度分別約92 m、143 m、53 m、133 m和75 m.實驗主要利用GPS的L1+L2觀測值和BDS的B1+B2 觀測值進行處理，電離層延遲和對流層延遲分別利用Klobuchar模型和Saastamoinen模型進行改正，模糊度固定采用最小二乘模糊度降相關(guān)平差(LAMBDA)方法.具體的解算策略如表1所示.

表1 解算策略

表2 粗差處理方案

表3 觀測數(shù)據(jù)的模糊度固定

限于篇幅，本文僅給出Test1的詳細數(shù)據(jù)處理結(jié)果分析.在衛(wèi)星定位中，可視衛(wèi)星數(shù)越多，衛(wèi)星幾何構(gòu)型越穩(wěn)定.圖1為Test1觀測數(shù)據(jù)的GPS+BDS的可視衛(wèi)星數(shù)和定位精度因子(PDOP).可以看出，觀測數(shù)據(jù)的可視衛(wèi)星數(shù)平均值為21，但衛(wèi)星數(shù)存在明顯的波動，接收衛(wèi)星信號不太穩(wěn)定，表明觀測質(zhì)量不佳.所有歷元的PDOP值都大于1，其PDOP均值約1.6，有些歷元的PDOP值超過3，因此整體看衛(wèi)星的空間幾何分布良好，但有時存在較差的空間結(jié)構(gòu).從整體看，隨著衛(wèi)星數(shù)的減少，PDOP的值隨之增加.圖2為Test1觀測數(shù)據(jù)的衛(wèi)星天空視圖.可以發(fā)現(xiàn)GPS和BDS存在部分低高度角衛(wèi)星，導致衛(wèi)星信號存在反射、衍射和被阻擋的現(xiàn)象.因此在數(shù)據(jù)處理過程中設(shè)置了截止高度角為15°，剔除了不符合要求的衛(wèi)星.

圖1 Test1數(shù)據(jù)的衛(wèi)星數(shù)和PDOP值Fig.1 Satellite number and PDOP value of Test1

圖2 Test 1數(shù)據(jù)的衛(wèi)星天空視圖Fig.2 Skyplot of Test 1

圖3和圖4分別給出了A方案的雙差偽距殘差和雙差載波殘差.由圖3可知，GPS和BDS的偽距殘差絕大多數(shù)集中在10 m之內(nèi)，其中GPS L2的偽距殘差較大，可以看到在某些歷元存有偽距粗差.由圖4可看出，GPS和BDS的載波殘差絕大多數(shù)集中在0.1 m之內(nèi).其中載波粗差沒有完全剔除，而且較大的載波粗差會影響歷元其他觀測值，對定位結(jié)果產(chǎn)生嚴重影響.總之，需要加入合適的粗差處理方法來提高定位精度.

圖3 A方案的雙差偽距殘差Fig.3 Double-differenced pseudorange residuals of scheme A

圖5、6、7分別描繪了各個方案的定位誤差結(jié)果.需要說明的是，定位結(jié)果中包含了固定解和浮點解的結(jié)果，以進一步整體評估其定位性能.由圖5可以看出，無任何粗差處理方法的A方案定位結(jié)果是不太理想的，尤其是U方向的波動較為明顯.對比圖5，圖6的3個數(shù)據(jù)探測方案都能探測到觀測值的粗差，并進行粗差剔除，其中B方案在U方向上絕大多數(shù)歷元的定位誤差都比其他兩個方案小.圖7中大迭代閾值(ω)的穩(wěn)健估計方案比小迭代閾值(ω)的穩(wěn)健估計方案在相同的閾值(k0、k1)下有更為明顯的收斂.其中在小閾值(k0、k1)時，E1方案在U方向上明顯比E方案的誤差小.同樣E1方案對比其他穩(wěn)健估計方案有更為明顯的粗差處理效果，在衛(wèi)星定位中表現(xiàn)了更好的定位性能.因此小閾值(k0、k1)結(jié)合大迭代閾值(ω)的穩(wěn)健估計方案更適用于復雜條件的粗差處理.

圖4 A方案的雙差載波殘差Fig.4 Double-differenced carrier residuals of scheme A

圖5 A方案的定位誤差Fig.5 Positioning error of scheme A

圖6 數(shù)據(jù)探測方案的定位誤差Fig.6 Positioning error of data snooping

表5 穩(wěn)健估計方案的定位RMSE統(tǒng)計

為驗證各方案粗差處理的有效性，表6和表7分別列出數(shù)據(jù)探測與穩(wěn)健估計的定位可用性，可以看到數(shù)據(jù)探測方案的可用性在任何的水平分量都優(yōu)于A方案.同樣穩(wěn)健估計的E、F方案的定位可用性略差于A方案，但其他方案的定位可用性比A方案有所提升，其中E1方案在水平分量為2.0 m時，定位可用性達到96.50%.上述結(jié)果同樣證明大迭代閾值(ω)配合小閾值(k0、k1)的穩(wěn)健估計方案可以更好地識別粗差，對觀測值權(quán)重的調(diào)整更為合理.

圖7 穩(wěn)健估計方案的定位誤差Fig.7 Positioning error of robust estimation

表6 數(shù)據(jù)探測方案的定位可用性

表7 穩(wěn)健估計方案的定位可用性

5 結(jié)語

觀測數(shù)據(jù)中含有異常粗差時，會對定位的精度和模糊度收斂產(chǎn)生嚴重的影響.本文系統(tǒng)研究了面向復雜條件時粗差處理的方法，針對復雜條件下的北斗/GNSS數(shù)據(jù)提出多種粗差處理方案，并以五組實驗數(shù)據(jù)分析各方案的模糊度收斂速率，著重以一組數(shù)據(jù)對比分析各方案的粗差處理效果和定位性能.

綜上，數(shù)據(jù)探測和穩(wěn)健估計的處理粗差效果相當，但由于數(shù)據(jù)探測的粗差處理相對更精細化，因此在本文實驗中數(shù)據(jù)探測可以表現(xiàn)出更好的適用性.