一種基于MAD改進的GNSS高程時間序列粗差探測方法

魯鐵定 何錦亮 徐華卿 賀小星

1 東華理工大學測繪工程學院，南昌市廣蘭大道418號，330013 2 江西理工大學土木與測繪工程學院，江西省贛州市紅旗大道86號，341000

受多路徑效應、地質構造活動、測站相關誤差以及外界環境影響，GNSS坐標時間序列中不可避免地存在粗差[1-3]，從而對建模精度和結果的可靠性產生重大影響。因此，探測并剔除GNSS坐標時間序列中的粗差就顯得尤為重要。

經典的統計量粗差探測方法有3σ準則、四分位距法(inter-quartile range，IQR)以及中位數絕對偏差法(median absolute deviation，MAD)。其中，3σ法是利用最小二乘法求得觀測序列的殘差，然后對殘差進行粗差探測。但最小二乘本身難以抵抗粗差，使得3σ法對殘差的中誤差估計有偏，從而影響粗差探測的可靠性[2]。IQR法相比3σ準則受粗差影響較小，具有更好的穩健性[4]；但對于離散度較大的數據會使估計值偏大，從而影響粗差探測的效果[5]。MAD法具有抗差性強、對粗差敏感的優點，相比3σ準則和IQR法具有更好的穩健性[6-7]，已廣泛應用于粗差探測領域[7-9]。

在GNSS坐標時間序列粗差探測中，首先需構建合適的時間序列模型，獲取殘差序列[2]。近年來，國內部分學者利用小波分析[10](wavelet analysis，WA)、奇異譜分析[1,5](singular spectrum analysis, SSA)、L1范數[3](L1-norm)獲取GNSS坐標時間序列的殘差向量，并構造粗差判別統計量進行粗差探測。這些方法的關鍵在于如何準確地提取原始序列的趨勢項和周期項，以分離出含粗差的殘差項。其中，奇異譜分析在選取滯后窗口時具有一定主觀性，不同的滯后窗口對信號的提取影響較大[11]；而文獻[3]中L1范數與IQR組合算法可能會對粗差產生“誤判”或“漏判”。小波分析具有局部化時頻分析能力和多分辨率分析特性，能準確提取時間序列的趨勢項和周期項，進而反映出時間序列的內在特征[12-13]。

因此，本文在MAD法基礎上結合小波分析，建立一種WT-MAD粗差探測方法，以提高傳統MAD法對GNSS高程時間序列中偏離度較小的粗差的探測能力。

1 原理及方法

1.1 MAD粗差探測方法

黃博華等[7]在傳統MAD方法基礎上對其數學表達式進行改進，得到式(1)形式的表達式。對于服從正態分布的觀測序列Xn={x1,x2,…,xi,…xn}，將觀測數據xi與其中位數K及中位數絕對偏差(MAD)數倍之和進行比較，若xi滿足

|xi-K|>n·MAD

(1)

則認為xi為粗差點。式中，K=Med{xi}，MAD=Med{|xi-K|/0.674 5}(Med表示求中位數)，常數n取值按實際需求確定，n值較大不易剔除偏離度較小的粗差，而n值較小可能產生誤判，一般取3～5[7,9]。在探測出粗差后，將粗差值設為0或其他值予以剔除。

1.2 WT-MAD粗差探測方法

GNSS高程時間序列可看作由不同頻率部分組成的信號，趨勢項、周期項主要集中在低頻信號中，噪聲和粗差干擾項主要集中在高頻信號中[14]。本文在傳統MAD方法基礎上，利用小波分析對其進行改進，主要思路如下：利用小波多尺度分解提取原始時間序列的低頻系數，將其重構為趨勢項和周期項；原始時間序列減去重構序列得到殘差序列，再利用MAD法對殘差序列進行粗差探測，進而確定粗差位置。具體步驟如下：

1)利用3σ準則剔除原始時間序列X(t)中偏離度較大的粗差，以避免小波分解和重構受到影響。對于剔除粗差后的缺失值，采用線性插值法插補得到X′(t)。

(2)

式中，RMSEj表示分解層數為j時原始信號與降噪信號間的均方根誤差，可由式(3)求得；r為相鄰層數均方根誤差之比，當r接近于1時，則取最佳層數為j或j+1。

(3)

(4)

4)利用MAD估計值探測殘差序列中的粗差，當殘差ri滿足

|ri-K|>n·MAD

(5)

時，則認為ri為粗差點，即原始數據中xi為粗差點。式中，K為殘差序列的中位數。步驟1)和4)中探測的粗差即為WT-MAD法探測出的所有粗差。

2 實驗與分析

2.1 模擬數據

為驗證WT-MAD方法的有效性，采用由趨勢項、周期項和噪聲項所組成的函數模型[15]對GNSS坐標時間序列進行建模，構造模擬時間序列，函數表達式如下：

y(t)=b+v0ti+

(6)

式中，i為坐標歷元時刻標識，y(ti)為GNSS測站某一分量ti時刻的坐標，b為橫軸截距，v0為線性速度，m0為周期性信號個數，fm為周期項頻率，am和bm表示頻率為fm的周期項對應的振幅，rti為隨機噪聲。

將利用式(6)模擬的高程方向坐標時間序列作為原始數據，各參數設置為：b=5.0，v0=2，m0=2，a1=b1=5，a2=b2=3，σwn=4。

向原始數據中加入粗差：首先模擬服從正態分布且標準差為6σwn(σwn為噪聲標準差)的隨機誤差序列；然后提取大于3σwn的數據，并將其隨機插入到原始數據中，得到一組被粗差污染的GNSS高程時間序列。模擬時間序列的時間跨度為2016～2020年，歷元間隔為1 d，總歷元數為1 827，粗差總數為115，粗差占總歷元比例為6.29%(圖1)。

圖1 模擬的高程時間序列

模擬實驗采用4種方法進行對比：LS-3σ法、LS-IQR法、LS-MAD法、WT-MAD法。前3種方法均基于式(6)，利用最小二乘法去除趨勢項和周期項以獲取殘差序列，再構造粗差判別式進行粗差探測。在LS-MAD法中，經篩選取n=3；在WT-MAD法中，經篩選取n=3，選用db4小波基函數，并利用式(2)求解最佳小波分解層數，相鄰層數均方根誤差之比見表1。由表可知，第4層和第5層間的均方根誤差之比最小，因此取小波分解層數為5。粗差探測結果如圖2和表2所示。

表1 相鄰分解層數均方根誤差之比

從圖2可以看出，LS-3σ法可探測出大部分粗差，但對偏離度較小的粗差探測效果有限，而LS-IQR法、LS-MAD法和WT-MAD法均可探測出絕大部分粗差。由表2可知，WT-MAD法粗差剔除率為98.3%，高于其他3種方法；雖然LS-IQR法和LS-MAD法均存在1個誤判，但這兩種方法的粗差探測率與WT-MAD法相近，這主要是因為模擬時間序列未考慮季節項、階躍等非線性變化，且加入的噪聲僅為白噪聲。在先驗模型準確的情況下，利用最小二乘法同樣能準確獲取模擬數據的殘差，從而得到較好的粗差探測效果。在利用WT-MAD法剔除粗差后，通過小波分析可得到去除趨勢項和周期項后的殘差，結果如圖3所示。從圖3可以看出，殘差序列為白噪聲序列，無明顯異常值，表明WT-MAD法可有效探測出模擬時間序列中的粗差。

圖2 LS-3σ法、LS-IQR法、LS-MAD法和WT-MAD法粗差探測結果

表2 各方法粗差探測結果對比

圖3 WT-MAD法剔除粗差后的殘差序列

2.2 實測數據

本文選用SOPAC(scrips orbit and permanent array center)提供的“Raw”和“Clean”類型的單天高程(U)時間序列作為實測數據，其中“Raw”為原始時間序列，“Clean”為刪除異常值的時間序列。

本次實驗選取LHAZ、BJFS、TWTF三個IGS站2006～2020年共15 a的“Raw”數據進行分析。在IGS站中，GNSS坐標時間序列數據缺失的情況較為常見[16]，而小波分析要求原始數據均勻采樣，因此在粗差探測前先利用線性插值法插補缺失值。得到完整時間序列后，以LHAZ站為例，分別利用小波變換與LS法獲取去除趨勢項和周期項后的殘差序列(圖4)。由圖可知，小波變換所得殘差趨近為白噪聲，而LS法獲取的殘差序列含有殘余的周期性信號，表明LS法難以抵御粗差和有色噪聲的影響，而小波變換在無需數據先驗信息的情況下，能夠更為準確地提取原始數據的趨勢項和周期項。

圖4 LHAZ站高程方向殘差序列

利用WT-MAD法對各IGS站的時間序列進行粗差探測，與模擬實驗相同，采用db4小波基函數，由式(2)求得各站最佳小波分解層數均為5，取n=3，探測結果見表3。

表3 LHAZ、BJFS、TWTF三個IGS站粗差探測結果

表3中粗差探測結果表明，3個IGS站均含有粗差，其中LHAZ站探測出的粗差最少，為100個；BJFS站最多，為104個。為驗證WT-MAD方法的有效性，以LHAZ站為例，得到該方法剔除粗差后的高程時間序列，并與SOPAC提供的“Raw”數據進行對比，結果如圖5～6所示。

圖5 LHAZ站“Raw”高程時間序列

圖6 WT-MAD法剔除粗差后的LHAZ站高程時間序列

從圖5～6可以看出，原始時間序列在歷元2006.5 a、歷元2009.0 a、歷元2014.0 a、歷元2016.0 a以及歷元2019.0 a附近均含有偏離度較大的粗差，而WT-MAD法可剔除這些粗差。

在利用WT-MAD法剔除LHAZ站原始高程時間序列的粗差后，通過小波變換可得到去除趨勢項和周期項后的殘差，結果如圖7(a)所示。圖7(b)為SOPAC提供的LHAZ站去除粗差、趨勢項和周期項后的殘差。對比圖7(a)和圖7(b)可知，圖7(b)中歷元2014.0 a和歷元2021.0 a附近仍含有粗差，而圖7(a)中這兩個歷元附近均無異常值，且殘差序列近似為白噪聲序列，表明WT-MAD法可更有效地剔除原始時間序列中的粗差。

圖7 LHAZ站殘差序列

為進一步驗證WT-MAD方法的有效性，分別利用LS-3σ法、LS-IQR法、LS-MAD法和WT-MAD法對表3中3個IGS站的數據進行粗差探測，并將利用小波分析提取的趨勢項和周期項作為參照，計算原始時間序列剔除粗差前后的RMSE，以此作為評判標準，結果如表4所示。從表4可以看出，3個測站中WT-MAD法探測出的粗差數量均多于LS-3σ法、LS-IQR法和LS-MAD法，且RMSE均小于其他3種方法，表明WT-MAD法可得到較好的探測效果。

表4 4種方法粗差探測結果及精度統計

3 結 語

本文針對GNSS高程時間序列非線性、不平穩性導致粗差探測困難的問題，構建一種基于MAD改進的WT-MAD粗差探測方法。該方法在進行粗差探測時不易受趨勢項和周期項等影響，同時可降低數據非對稱分布對MAD法的影響，從而可有效探測粗差。模擬實驗和實測結果均表明，相比LS-3σ法、LS-IQR法和LS-MAD法，WT-MAD法可得到較好的探測效果，說明本文方法具有適用性和有效性。

在進行小波分析時，小波基函數及分解層數根據經驗來選取，且需要進行重復實驗確定，這會影響粗差探測效率，且在一定程度上具有主觀性；若分解層數過大可能會造成粗差誤判，過小則可能造成漏判。因此，如何快速、準確地選取小波基函數及分解層數還需進一步研究。