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

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

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

受多路徑效應(yīng)、地質(zhì)構(gòu)造活動、測站相關(guān)誤差以及外界環(huán)境影響，GNSS坐標(biāo)時間序列中不可避免地存在粗差[1-3]，從而對建模精度和結(jié)果的可靠性產(chǎn)生重大影響。因此，探測并剔除GNSS坐標(biāo)時間序列中的粗差就顯得尤為重要。

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

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

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

1 原理及方法

1.1 MAD粗差探測方法

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

|xi-K|>n·MAD

(1)

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

1.2 WT-MAD粗差探測方法

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

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

(2)

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

(3)

(4)

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

|ri-K|>n·MAD

(5)

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

2 實驗與分析

2.1 模擬數(shù)據(jù)

為驗證WT-MAD方法的有效性，采用由趨勢項、周期項和噪聲項所組成的函數(shù)模型[15]對GNSS坐標(biāo)時間序列進(jìn)行建模，構(gòu)造模擬時間序列，函數(shù)表達(dá)式如下：

y(t)=b+v0ti+

(6)

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

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

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

圖1 模擬的高程時間序列

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

表1 相鄰分解層數(shù)均方根誤差之比

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

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

表2 各方法粗差探測結(jié)果對比

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

2.2 實測數(shù)據(jù)

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

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

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

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

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

表3中粗差探測結(jié)果表明，3個IGS站均含有粗差，其中LHAZ站探測出的粗差最少，為100個；BJFS站最多，為104個。為驗證WT-MAD方法的有效性，以LHAZ站為例，得到該方法剔除粗差后的高程時間序列，并與SOPAC提供的“Raw”數(shù)據(jù)進(jìn)行對比，結(jié)果如圖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站原始高程時間序列的粗差后，通過小波變換可得到去除趨勢項和周期項后的殘差，結(jié)果如圖7(a)所示。圖7(b)為SOPAC提供的LHAZ站去除粗差、趨勢項和周期項后的殘差。對比圖7(a)和圖7(b)可知，圖7(b)中歷元2014.0 a和歷元2021.0 a附近仍含有粗差，而圖7(a)中這兩個歷元附近均無異常值，且殘差序列近似為白噪聲序列，表明WT-MAD法可更有效地剔除原始時間序列中的粗差。

圖7 LHAZ站殘差序列

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

表4 4種方法粗差探測結(jié)果及精度統(tǒng)計

3 結(jié) 語

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

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