顧及三維方向的GNSS坐標時間序列粗差探測方法

許 煒 呂志平 董 行 鄺英才 楊凱淳 黃 輝

1 信息工程大學地理空間信息學院，鄭州市科學大道62號，450001 2 中國人民解放軍96712部隊，江西省景德鎮市，333300 3 哈爾濱工業大學(深圳)空間科學與應用技術研究院，深圳市平山一路6號，518055 4 中國人民解放軍61206部隊，北京市，100042 5 中國人民解放軍96711部隊，安徽省池州市，247100

針對GNSS坐標時間序列粗差探測與剔除的研究很多，Pernetti[1]利用DIA(detection identification adaptation)方法對GPS坐標時間序列的粗差和中斷進行檢測，結論顯示，早期的GPS坐標時間序列中大約有6%的數據為粗差；Mao等[2]及李曉光等[3]采用3倍中誤差(3σ)準則進行粗差探測，但3σ準則易受粗差的影響；為增強粗差探測的穩健性，四分位距(interquartile range, IQR)準則得到了更加廣泛的應用，其對粗差有一定的抵御能力。IQR準則也衍生出不少新的算法，如經典的最小二乘結合IQR準則的方法(LS-IQR)[4]、抗差1范數結合IQR準則的方法(L1-ModIQR)[5]、數據驅動的奇異譜分析技術結合IQR準則的方法(SSA-IQR)[6]及在頻域上與小波分析結合的方法(WA-IQR)[7]等；崔太岷等[8]還提出基于假設檢驗的方法對觀測數據進行粗差探測。近年來，隨著人工智能的廣泛應用，甚至發展出了利用機器學習和深度學習對各種時間序列進行異常探測的算法[9]。

本文通過檢驗常用的3σ準則和IQR準則的探測效果，對三維GNSS大地坐標時間序列進行粗差探測，并分析三維方向的探測效果。

1 GNSS坐標時間序列粗差的特點

三維GNSS坐標時間序列的粗差是空間位置的粗差，如圖1(單位mm)所示，假設對某點進行多次GNSS定位，并模擬其中4次定位結果的粗差。

圖1 三維GNSS定位

為方便說明，以二維平面坐標為例。假設圖2(單位mm)為GNSS多次測量某點的定位結果，A、B、C、D為4個模擬粗差。

圖2 二維GNSS定位

以3σ準則為例，如圖2所示，假設通過平面定位某點100次，N方向和E方向的中誤差均為1 mm，并服從正態分布，且分布于3σ之外的數據不大于1%，其中1σ區域為中間淡藍色區域，3σ為外圈淺黃色區域，淺紅色正方形區域為通過N、E方向判定不是粗差而分布于3σ之外的區域。A、B、C、D為4個人為加入的偏離超過3σ的粗差點，如僅考慮N方向，A與B可被識別為粗差，但C與D將不能被正確識別；同樣，如果僅考慮E方向，B和D無法被正確探測出來。雖然2個方向的探測均不能識別出在3σ之外的D點，但這樣的區域占比很小，僅同時探測N、E方向可大大簡化計算量，且不損失過多的探測效率。

類似于平面位置的3σ準則時間序列粗差探測，三維GNSS坐標的粗差探測應探測各個方向的偏差是否大于某個統計量指標，需要求出模型化后三維各方向的統計量信息。但這種處理方法計算復雜，若簡化為僅考慮N、E、U三個相互正交方向的偏差，計算量可被大大簡化，且能覆蓋大多數情形，所以本文采用各個分量分別探測出來的時間并集作為粗差發生歷元。

為更有效地檢測GNSS坐標時間序列中的粗差，可在時域或頻域上對時間序列進行建模，然后針對得到的殘差序列進行探測。由于在頻域上對時間序列進行分析要求采樣均勻[10]，但GNSS時間序列中存在不少缺失數據，破壞了其均勻性，進行插值則可能會引入新的粗差。所以，本文分別依據3σ準則和IQR準則，利用迭代最小二乘法求取參數和剔除粗差。

2 綜合N、E、U三方向的迭代最小二乘粗差探測方法

在實際應用中，時間序列的自相關性隨觀測歷元間隔的增大而降低[11]，坐標時間序列中存在與時間相關的有色噪聲成分，而且對于跨度較長的坐標時間序列，其前后的年周期和半年周期振幅往往不一樣，所以一般采用移動開窗法對粗差進行探測。考慮到GNSS坐標時間序列中的年周期信號最明顯[12]，所以窗口及步長均設為1 a，不足1 a的數據以離目標數據最近的1 a數據作為探測窗口。

僅靠1次探測在很多情況下不能成功探測出所有粗差，所以在單方向上采取迭代求取參數和剔除粗差的方式。迭代的優勢在于可以逐步剔除粗差，并消除粗差對最小二乘法估計參數的影響，最終得到正確的殘差序列。而三維粗差探測為N、E、U方向探測出的粗差并集，可大幅降低單方向漏判的概率，具體探測步驟如下：

1)輸入觀測方程和權。

2)利用最小二乘法估計出GNSS時間序列模型，求取參數和殘差。

3)按3σ準則或IQR準則探測和剔除N、E、U方向的粗差。

3σ準則是指當大量GNSS坐標時間序列數據模型化后的殘差服從正態分布時，誤差分布于±3σ之內的數據占99.73%，所以將分布于3σ之外的數據認為是粗差造成的，其誤判比例很小。具體判別式為：

(1)

(2)

4)重復步驟1)～3)，直至探測出的粗差個數為0。

5)求取N、E、U三方向探測所得的粗差并集。

3 算例分析

3.1 模擬數據分析

對于單站單方向GNSS坐標時間序列，一般采取式(3)進行建模[10]：

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

(3)

式中，ti為序列的歷元，單位a；a對應于截距，其值與預設的初始坐標位置有關；b為站點坐標線速度；c、d和e、f分別為年周期項和半年周期項系數，能夠完整描述周期函數的振幅和初始相位；gi為由各種原因引起的階躍式坐標突變；Thj為階躍發生的開始時刻；H為階躍函數，發生突變前H值為0，發生突變后H值為1；vti為觀測噪聲，可表示成不同噪聲模型的組合。

為說明粗差探測方法的有效性并驗證三維GNSS坐標時間序列的粗差剔除效果，對模擬數據進行處理。為使問題不過于復雜，又具有一定的代表性，實驗根據北京房山(BJFS)站的線速度、周年和半周年振幅模擬了2009～2019年共10 a的三維GNSS坐標時間序列，截距(a)在初始位置N、E、U方向坐標均設為0，加入中誤差均為3 mm的白噪聲，且不添加階躍信號，依據式(3)的參數模擬結果如表1(單位mm)所示。

表1 模擬時間序列參數值

粗差模擬過程如下：1)將三維方向的噪聲按照勾股定理求取空間位置誤差的模長，將得到的模長放大6倍，并從中隨機選取200個大于3倍空間位置誤差模長的時間點作為粗差，粗差占比為5.5%。粗差的天頂距φ在0～π中均勻分布，水平方向θ在0～2π中均勻分布。2)將粗差的模長分別乘以cosφ、sinφ·cosθ及sinφ·sinθ，并加入U、N、E方向選取的200個粗差時間點，得到各方向被粗差污染的時間序列。原序列和被粗差污染的序列如圖3所示。

圖3 模擬數據及粗差

為檢驗2種準則的檢測性能和GNSS三維位置粗差的探測效果，設計模擬實驗對2種方案進行對比：1)利用最小二乘結合3σ準則方法(LS-3σ)進行單方向及綜合三方向探測；2)利用最小二乘結合IQR準則方法(LS-IQR)進行單方向及綜合三方向探測。對比各方法探測的粗差個數和誤判個數。

各方案粗差探測結果如表2所示，3σ準則和IQR準則在各方向的探測結果如圖4所示，綜合三方向的探測結果如圖5所示。

表2 模擬GNSS坐標時間序列粗差探測效果

圖4 3σ準則與IQR準則各方向粗差探測效果

圖5 3σ準則與IQR準則綜合三方向粗差探測誤判情況

結合表2和圖4可以看出，若僅考慮單方向，2種準則均能較好地識別出偏離較大的粗差，但無法識別偏離較小的粗差，所以單方向的探測率均不高。其中，3σ準則水平方向的探測率在55%左右，存在5%左右的誤判，而垂直方向的探測率在80%左右，誤判率也在5%左右；IQR準則探測率較低，水平方向在45%左右，垂直方向約為75%，不過IQR準則在本次實驗中沒有出現誤判的情況，可認為其誤判率顯著低于3σ準則。低誤判率的主要原因為IQR統計量并未進行標準化，相當于對粗差判定的置信區間進行了壓縮。2種方法在垂直方向的探測效率均好于水平方向，主要原因是在粗差的模擬過程中，垂直分量粗差的占比要大于水平分量。將各方向探測得到的粗差求取并集后，3σ準則可探測所有粗差的發生時刻，但誤判率也達到了13.5%；IQR準則的探測率超過98%，沒有出現誤判。圖5中的粗差誤判和漏判歷元均是在3個方向上的，粗差與模擬正常值沒有明顯差別。綜合來看，IQR準則要好于3σ準則，其以稍低的探測率換取了低誤判率。

3.2 實測數據分析

根據模擬實驗結果，實測數據采用IQR準則進行粗差探測。實測數據為SOPAC(Scrips Orbit and Permanent Array Center)利用GAMIT軟件處理的全球參考站坐標時間序列，選取其中9個近年來數據相對完整且穩定的站點作為分析對象，時間跨度為2009.001 4～2018.998 6(2009-01-01～2018-12-31)，具體歷元數因各站點的數據缺失程度不一而不同。

SOPAC提供4類坐標時間序列：1)僅去除較大粗差的原始序列(raw)，即平差前對水平方向偏離大于1 000 mm、垂直方向大于3 000 mm的觀測值進行剔除；2)移除了均值、粗差、階躍的干凈時間序列(cleaned trended)；3)在干凈時間序列的基礎上移除趨勢的時間序列(cleaned detrended)；4)對干凈時間序列進行主成分分析濾波得到的過濾時間序列(filtered)。其中，cleaned trended數據的粗差剔除策略是觀測值減去計算值(O-C)，以檢測模型化后的殘差是否大于所設置的閾值，其閾值在N、E、U方向分別為25.0 mm、25.0 mm和35.0 mm。

本文對SOPAC提供的raw數據進行粗差探測，并與提供的cleaned trended數據進行對比來進行驗證。為保證模型的一致性，在進行粗差探測前，通過查詢cleaned trended數據中的階躍信息進行模型改正，再利用IQR準則對時間序列的三方向進行粗差探測，并求得并集，探測結果見表3。由表3可知，每個測站都存在數據缺失的情況，缺失比例在2%～32%，不滿足均勻采樣的條件，表明實測的IGS站原始序列數據不適合直接在頻域上進行分析。利用本文方法對9個IGS站數據進行粗差探測得出，原始坐標時間序列中含有的粗差占比在0.1%～4.7%之間。

表3 IGS參考站粗差探測效果

為詳細說明探測效果，以BUCU站為例，對其三維位置粗差探測效果進行分析，檢驗探測結果的有效性，如圖6所示。可以看出：

圖6 BUCU站粗差探測效果

1)由于水平方向存在非常明顯的線性趨勢，局部特征容易被線性趨勢掩蓋，僅靠目測不易判斷其中存在的粗差。在僅考慮垂直方向的情況下，2個水平方向的明顯粗差未能被識別出來。綜合探測3個方向后，識別的粗差個數明顯增多，包括很多從圖上看起來不像粗差的值。

2)根據模擬實驗的結果，在多個方向上同時探測出某個歷元存在粗差的比例大概在70%左右，而BUCU站的實測數據中僅有10%左右，說明存在不小比例的誤判。其原因可能是存在一些未探明的階躍信號，如2010.053 4～2010.072 6連續8 d被標記為粗差，通過檢查N方向數據發現，其很有可能是某種類型的階躍，表明在cleaned trend數據中還存在沒有被探明的階躍。

3)坐標時間序列中的噪聲并不是單純的白噪聲，其先驗隨機模型不準確，導致其與模擬實驗結果存在差異。為更明顯地對粗差剔除效果進行對比，實驗通過最小二乘法估計出利用IQR準則剔除粗差后序列的趨勢及截距，將得到的去趨勢數據與SOPAC提供的cleaned detrend數據進行對比，結果如圖7所示。從圖7可以看出，相比于cleaned detrend數據，IQR準則綜合三方向剔除粗差后的時間序列不僅剔除了孤立的異常點，細節特征也更加明顯。

圖7 IQR準則綜合三方向剔除粗差后時間序列與cleaned detrend時間序列對比

4 結 語

本文對比了最小二乘分別結合3σ準則和IQR準則對三維GNSS坐標時間序列進行粗差探測的方法。由模擬實驗結果可見，綜合三方向的探測結果明顯好于僅考慮單方向的結果；3σ準則的探測效率高，但存在較多的誤判；IQR準則探測效率略低，但沒有出現誤判的情況。總體來看，IQR準則要好于3σ準則。

然而實測數據情況要更為復雜，其中含有大量的階躍信號，如果未準確建模，將會對參數估計產生很大影響，進而影響探測效果。而在SOPAC提供的干凈時間序列中依然存在很多沒有被探明的階躍，造成了探測中粗差與階躍的混淆，因此將階躍和粗差進行區分十分重要。

另外，GNSS坐標時間序列中的噪聲成分復雜，其隨機模型不準確會使參數估值和殘差的計算受到很大影響，進而影響粗差探測的效果。獲取更加準確的函數模型和隨機模型也有助于提高粗差探測的準確性。

GNSS坐標時間序列中粗差產生的原因常被籠統地歸結為儀器故障或測量條件的限制，對粗差產生的原因進行定性和定量分析，有助于在利用GNSS觀測數據獲取坐標時間序列階段避免或減弱粗差的影響。