基于SSA-Sn法的GNSS坐標時間序列粗差探測

徐華卿 魯鐵定 賀小星 周世健 陶 蕊

1 東華理工大學測繪工程學院，南昌市廣蘭大道418號，330013 2 江西理工大學土木與測繪工程學院，江西省贛州市紅旗大道86號，341000 3 南昌航空大學校長辦公室，南昌市豐和南大道696號，330063

以GNSS為主的衛星導航定位基準站網近年來發展迅速，為大地測量的應用提供了高精度的空間基準基礎設施[1]。GNSS坐標時間序列蘊含豐富的物理信息，包括測站長期受區域內構造應力影響產生的線性變化以及地球物理效應與其他外部環境等導致的非線性變化[2]。由于測站會受到電離層延遲、多路徑效應、環境負載、軌道異常等多種誤差影響[3]，GNSS坐標時間序列在呈現非線性變化的同時也存在異常值[4]。這些異常值在整體或局部數據中表現出較大的波動，無法準確反映出時間序列的變化特征，對數據分析與應用產生較大的負面影響。

經典的粗差探測法有3σ法[5]和MAD(median absolute deviation)法[6-7]。3σ法的原理是利用最小二乘法求得觀測序列的殘差，再對殘差序列進行粗差探測。由于最小二乘法本身含有粗差，因此3σ法的殘差中誤差估計值不夠準確，從而影響粗差探測的可靠性[8-11]。MAD方法的運算崩潰污染率可達50%，其異常值探測相比于3σ法更加穩健，但MAD法需要滿足樣本對稱分布的條件[12-13]。基于上述討論，諸多學者在此基礎上引入新的準則和方法以提高粗差探測的效率和準確率。明鋒等[8]提出L1范數與IQR組合探測方法，使模型化后的殘差探測準確率更高；吳浩等[5]提出一種利用小波改進的3σ粗差探測方法，既能準確提取變形趨勢，也能提高3σ在坐標數據中的適用性；Wang等[14]發現有效提取趨勢項和周期項至關重要， SSA在這方面具有一定優勢，能準確提取出殘差項；蔡曉軍等[15]研究表明，在提取趨勢項和周期項方面SSA優于小波分析法。

綜合以上異常值探測方法的優勢與不足并結合GNSS坐標時間序列的數據特征，建立一種SSA與Sn估計量[14]相結合的異常值探測方法。選用模擬的時間序列數據以及中國部分IGS站的實測數據進行實驗，結果表明，在粗差探測方面SSA-Sn法的探測效果優于3σ法和MAD法。

1 原理及方法

1.1 SSA數學模型

(1)

C=VΛVT

(2)

將C的特征值進行降序排列，得到λk，Λ為λk構成的對角矩陣，正交矩陣V的行向量為vk。主成分矩陣為：

(3)

A的第k行ak為第k個主成分，其中ak,i為：

(4)

式中，vk,j是vk的第j個元素。時間序列中第k個重構分量為：

(5)

(6)

1.2 交叉驗證法確定參數

采用交叉驗證法確定SSA需要的窗口長度L與主成分個數P兩個參數。首先給出兩個參數的初始值，并在原始數據中隨機挑選部分數據用于驗證，驗證數據處用0填充，即可得到訓練數據。隨后將訓練數據進行奇異譜分析，計算驗證部分模型導出值與參考值的均方根誤差RMSE，RMSE最小時的窗口長度L與主成分個數P即為最優參數。Wang 等[14]表明，2～3 a的窗口長度適合提取年周期及半年周期分量，因此本文采用2 a作為窗口長度的初值。

1.3 Sn尺度估計量

Sn是一個由Rousseeuw和Croux提出的穩健統計量，Sn估計量具有與MAD法相同的穩健性，并且不依賴于對稱分布。假設一組數據為{x1,x2,…,xn}，則Sn估計量被定義為：

Sn=cmedi{medi|xi-xj|}

(7)

式中，常數c是為了滿足Fisher一致性的調節因子，其值為1.192 6；外中位數是一個秩為[(n+1)/2]階的順序統計量，而內中位數是一個秩為[n/2]+1的順序統計量。

1.4 SSA-Sn粗差探測原理

在對GNSS坐標時間序列進行粗差探測時，可將時間序列數據視為一組被噪聲和粗差污染的趨勢信號。SSA相較于小波分析能夠更加準確地提取趨勢項、周期項等主要的時間序列特征，從而分離出殘差項；Sn與MAD法對異常值探測的穩健性相同且都優于3σ法，估計量崩潰點均為50%，不易受到粗差的污染，但Sn無需依賴數據的對稱分布，因此本文采用SSA-Sn法進行粗差探測。具體步驟如下：

1)輸入坐標時間序列原始數據集X(t)，使用交叉驗證法確定最優的窗口長度L與主成分個數P；

3)原始數據與重構數據作差，得到含有噪聲以及粗差的殘差序列：

是“煙富3號”的早熟株變。2015年通過貴州省品種審定委員會審定。其主要的優良變異體現在成熟期上，在貴州長順8月中下旬成熟。果實近圓形，果形指數0.87，平均單果重210～250克。果面平滑，稀有果粉，果皮底色淡黃，條紅，著色面積85%以上。果肉淡黃，質地硬脆，肉質細，汁液多，風味酸甜，香氣濃，品質上等。可溶性固形物14.7%～15.9%。

(8)

式中，v為殘差項；

4)采用Sn估計量對殘差序列v進行粗差探測，判別式為：

|vi-median(v)|>3Sn

(9)

式中，median(v)為殘差序列的中值。

2 實驗結果與分析

2.1 模擬GNSS站坐標時間序列

為檢驗粗差探測方法的有效性，根據GNSS單站、單分量坐標時間序列函數模型和隨機模型構造包含趨勢項、周期項和噪聲項的模擬時間序列，其函數表達式為：

(10)

式中，ti為觀測時間，a為測站時間序列的起始位置，b為測站運動的線性速度，c、d、e、f分別為測站周年運動和半周年運動的振幅，vi為GPS坐標時間序列中的噪聲。為更加真實地模擬GPS時間序列噪聲，采用式(11)生成白噪聲和有色噪聲：

(11)

式中，w(ti)是均值為0、方差為1的白噪聲，f1和f2均為0.2。

表1 模擬數據統計

圖1 模擬時間序列Fig.1 Simulated time series

2.2 模擬數據窗口長度L和主成分個數P的確定

圖2 模擬數據的RMSE隨窗口長度L與主成分個數P的變化情況Fig.2 The RMSE of the simulated data varies with the length of the window and the number of principal components

由圖2可見，模擬數據的RMSE對主成分個數P的變化比對窗口長度L的變化更加敏感，在主成分個數0～10范圍內每條曲線的RMSE都達到了最小值。經交叉驗證后，選取最優窗口長度L為830、主成分個數P為3。

2.3 模擬數據結果分析

圖3 原始序列與重構序列Fig.3 Original sequence and reconstructed sequence

圖4 提取的殘差序列Fig.4 The extracted residual sequence

利用3σ法、MAD法和SSA-Sn法對模擬數據進行粗差探測(圖5)。由圖可見，由于數據標準差σ會受到粗差的污染，不具有穩健性，因此3σ法只能探測出偏離較大的粗差，對偏離程度較小粗差的探測效果有限。相較于3σ法，MAD法具有50%的抗差穩健性，能夠探測出絕大部分粗差。為進一步體現粗差探測的效果，進行3種方法的探測效果統計(表2)。由表可知，3σ法的探測數為43、準確率為38.7%；MAD法的探測數為104、準確率為93.6%；SSA-Sn法的探測數為112，雖然該方法存在1個誤判，但能夠探測出所有的粗差點，其準確率相較于3σ法和MAD法分別提升61.3百分點和6.4百分點，因此SSA-Sn法對于偏離程度較大或較小的粗差都具有很好的探測效果。

圖5 3種方法的粗差探測結果Fig.5 Gross error detection results of three methods

表2 3種方法粗差探測對比

圖6 原始序列與重構序列 and reconstructed sequence

圖7 提取的殘差序列Fig.7 The extracted residual sequence

圖8 SSA-Sn法粗差探測結果Fig.8 SSA-Sn method gross error detection results

2.4 實測數據實驗結果分析

為進一步驗證SSA-Sn法探測粗差的有效性，使數據更加貼合實際，選用SOPAC GPS提供的原始時間序列數據，選取BJFS、SHAO、URUM三個測站高程U方向上2009～2015年共計6 a的實測數據進行處理，圖9為3個測站U方向的原始時間序列數據。由圖9(a)和圖9(b)可見，原始數據中粗差在時域和大小上分布范圍較廣。

圖9 3個測站原始數據Fig.9 Raw data of three stations

分別采用3σ法、MAD法和SSA-Sn法對3個測站的數據進行粗差探測(圖10)。由BJFS站結果可知，2009～2010年的5個數據與整體數據存在明顯偏離，SSA-Sn法成功探測出2009～2010年處的粗差，但3σ法和MAD法未能全部探測出；由URUM站結果可知，3σ法和MAD法將2011年處的數據判定為粗差，但SSA-Sn法并未進行判定。結合圖像分析后得知，URUM站2011年處的數據符合整體數據的周期趨勢，不存在偏離，因此不屬于粗差。綜上所述，3σ法和MAD法均能準確探測出偏離程度較大的粗差，而SSA-Sn法對于偏離程度較小的粗差數據的探測效果更加優越。由于利用剔除粗差個數來評判粗差探測方法的優劣具有一定的局限性，本文采用SSA重構后的時間序列作為參照，以3種方法剔除粗差后的RMSE來評判粗差探測的效果(表3)。由表可見，經SSA-Sn法剔除粗差后3個測站數據的RMSE分別為4.30 mm、4.03 mm和4.90 mm，均小于MAD法和3σ法處理后3個測站數據的RMSE。因此3種方法都能探測出一定數量的粗差，但SSA-Sn法的粗差探測效果最優。

圖10 3個測站3種方法的粗差探測結果Fig.10 Gross error detection results of three methods at three stations

表3 3種方法的粗差探測數量與剔除后的RMSE

3 結 語

1)3σ法和MAD法的粗差探測率分別為38.7%和93.6%，而SSA-Sn法雖然存在1個誤判，卻能將模擬的粗差全部探測出來，并且對于噪聲模型為白噪聲+閃爍噪聲+隨機游走噪聲這種貼合實際情況的模擬數據同樣具有適用性。

2)采用SOPAC GPS提供的BJFS、SHAO、URUM三個測站高程U方向的原始數據進行實驗，從探測粗差個數以及剔除粗差后數據的RMSE兩方面證明，SSA-Sn法的粗差探測效果優于3σ法和MAD法。

今后需要結合GNSS時間序列和奇異譜分析構建系統的奇異譜分析方法，從而達到更高效的粗差探測效果。