基于SSA-Sn法的GNSS坐標(biāo)時(shí)間序列粗差探測(cè)

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

1 東華理工大學(xué)測(cè)繪工程學(xué)院，南昌市廣蘭大道418號(hào)，330013 2 江西理工大學(xué)土木與測(cè)繪工程學(xué)院，江西省贛州市紅旗大道86號(hào)，341000 3 南昌航空大學(xué)校長(zhǎng)辦公室，南昌市豐和南大道696號(hào)，330063

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

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

綜合以上異常值探測(cè)方法的優(yōu)勢(shì)與不足并結(jié)合GNSS坐標(biāo)時(shí)間序列的數(shù)據(jù)特征，建立一種SSA與Sn估計(jì)量[14]相結(jié)合的異常值探測(cè)方法。選用模擬的時(shí)間序列數(shù)據(jù)以及中國(guó)部分IGS站的實(shí)測(cè)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，結(jié)果表明，在粗差探測(cè)方面SSA-Sn法的探測(cè)效果優(yōu)于3σ法和MAD法。

1 原理及方法

1.1 SSA數(shù)學(xué)模型

(1)

C=VΛVT

(2)

將C的特征值進(jìn)行降序排列，得到λk，Λ為λk構(gòu)成的對(duì)角矩陣，正交矩陣V的行向量為vk。主成分矩陣為：

(3)

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

(4)

式中，vk,j是vk的第j個(gè)元素。時(shí)間序列中第k個(gè)重構(gòu)分量為：

(5)

(6)

1.2 交叉驗(yàn)證法確定參數(shù)

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

1.3 Sn尺度估計(jì)量

Sn是一個(gè)由Rousseeuw和Croux提出的穩(wěn)健統(tǒng)計(jì)量，Sn估計(jì)量具有與MAD法相同的穩(wěn)健性，并且不依賴于對(duì)稱分布。假設(shè)一組數(shù)據(jù)為{x1,x2,…,xn}，則Sn估計(jì)量被定義為：

Sn=cmedi{medi|xi-xj|}

(7)

式中，常數(shù)c是為了滿足Fisher一致性的調(diào)節(jié)因子，其值為1.192 6；外中位數(shù)是一個(gè)秩為[(n+1)/2]階的順序統(tǒng)計(jì)量，而內(nèi)中位數(shù)是一個(gè)秩為[n/2]+1的順序統(tǒng)計(jì)量。

1.4 SSA-Sn粗差探測(cè)原理

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

1)輸入坐標(biāo)時(shí)間序列原始數(shù)據(jù)集X(t)，使用交叉驗(yàn)證法確定最優(yōu)的窗口長(zhǎng)度L與主成分個(gè)數(shù)P；

3)原始數(shù)據(jù)與重構(gòu)數(shù)據(jù)作差，得到含有噪聲以及粗差的殘差序列：

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

(8)

式中，v為殘差項(xiàng)；

4)采用Sn估計(jì)量對(duì)殘差序列v進(jìn)行粗差探測(cè)，判別式為：

|vi-median(v)|>3Sn

(9)

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

2 實(shí)驗(yàn)結(jié)果與分析

2.1 模擬GNSS站坐標(biāo)時(shí)間序列

為檢驗(yàn)粗差探測(cè)方法的有效性，根據(jù)GNSS單站、單分量坐標(biāo)時(shí)間序列函數(shù)模型和隨機(jī)模型構(gòu)造包含趨勢(shì)項(xiàng)、周期項(xiàng)和噪聲項(xiàng)的模擬時(shí)間序列，其函數(shù)表達(dá)式為：

(10)

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

(11)

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

表1 模擬數(shù)據(jù)統(tǒng)計(jì)

圖1 模擬時(shí)間序列Fig.1 Simulated time series

2.2 模擬數(shù)據(jù)窗口長(zhǎng)度L和主成分個(gè)數(shù)P的確定

圖2 模擬數(shù)據(jù)的RMSE隨窗口長(zhǎng)度L與主成分個(gè)數(shù)P的變化情況Fig.2 The RMSE of the simulated data varies with the length of the window and the number of principal components

由圖2可見，模擬數(shù)據(jù)的RMSE對(duì)主成分個(gè)數(shù)P的變化比對(duì)窗口長(zhǎng)度L的變化更加敏感，在主成分個(gè)數(shù)0～10范圍內(nèi)每條曲線的RMSE都達(dá)到了最小值。經(jīng)交叉驗(yàn)證后，選取最優(yōu)窗口長(zhǎng)度L為830、主成分個(gè)數(shù)P為3。

2.3 模擬數(shù)據(jù)結(jié)果分析

圖3 原始序列與重構(gòu)序列Fig.3 Original sequence and reconstructed sequence

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

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

圖5 3種方法的粗差探測(cè)結(jié)果Fig.5 Gross error detection results of three methods

表2 3種方法粗差探測(cè)對(duì)比

圖6 原始序列與重構(gòu)序列 and reconstructed sequence

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

圖8 SSA-Sn法粗差探測(cè)結(jié)果Fig.8 SSA-Sn method gross error detection results

2.4 實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)結(jié)果分析

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

圖9 3個(gè)測(cè)站原始數(shù)據(jù)Fig.9 Raw data of three stations

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

圖10 3個(gè)測(cè)站3種方法的粗差探測(cè)結(jié)果Fig.10 Gross error detection results of three methods at three stations

表3 3種方法的粗差探測(cè)數(shù)量與剔除后的RMSE

3 結(jié) 語(yǔ)

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

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

今后需要結(jié)合GNSS時(shí)間序列和奇異譜分析構(gòu)建系統(tǒng)的奇異譜分析方法，從而達(dá)到更高效的粗差探測(cè)效果。