王民頓 尚俊娜

1 杭州電子科技大學通信工程學院，杭州市白楊街道2號大街1158號，310018

常用的GNSS時間序列粗差剔除方法主要有3σ法、中位數(shù)(MAD)法、四分位距(inter quartile range，IQR)法等[1-6]，但這3種方法都有一個共同的缺陷——數(shù)據(jù)剔除的效果在很大程度上受限于數(shù)據(jù)長度，以至于無法把握真實的數(shù)據(jù)趨勢[7]。為解決這個問題，可以采用滑動窗口的方式一段一段地剔除數(shù)據(jù)，但這會增加更多的限制，而且還會受到窗口選取的影響。

本文針對窗口難以選取的問題提出一種基于小波分析的一階導(dǎo)數(shù)粗差剔除法，經(jīng)過仿真模型及實際驗證，解決了傳統(tǒng)粗差探測算法在數(shù)據(jù)處理中的過剔除現(xiàn)象，在小波分析過程中能準確提取到信號的真實形變趨勢，盡可能多地區(qū)分出粗差點；同時，在粗差點采用廣義延拓插值補點，兼顧了監(jiān)測數(shù)據(jù)的連續(xù)性。

1 GNSS形變監(jiān)測時間序列的粗差數(shù)據(jù)處理方法

1.1 基于小波分析的一階導(dǎo)數(shù)粗差剔除法

基于小波分析的一階導(dǎo)數(shù)粗差剔除法步驟如下：

1)首先引入一階導(dǎo)數(shù)分析信號趨勢項。將信號求一階導(dǎo)數(shù)后，會得到類似高尺度下的小波系數(shù)，借用小波閾值的思想，原本小波閾值函數(shù)是將信號中低幅值的小波系數(shù)進行過濾，現(xiàn)在則認為大于閾值的一階導(dǎo)數(shù)點是粗差點。

2)利用小波原有的minimaxi法則求出閾值，以分離異常導(dǎo)數(shù)值，這里采用經(jīng)驗方程。如果要達到更好的效果，可以在閾值的設(shè)定方面進行拓展，這與小波閾值的設(shè)計相同，將其標志成為異常點，大部分異常點會在此被剔除。閾值計算公式為：

(1)

3)將剔除后的數(shù)據(jù)只進行1次3σ剔除，此處也可以選取大一點的判決門限(如5σ等)防止剔除了有用信號。

4)最后借用小波分析的手段對數(shù)據(jù)進行粗差剔除，其主要方式是通過小波分解得到低頻趨勢項，將原信號和低頻趨勢項作差得到殘差?？紤]到重復(fù)計算趨勢項會大幅增加計算量，這里只求1次趨勢項，對獲取到的殘差進行拉依達準則(3σ準則)計算。將異常點置0(這里的異常點指的是偏離趨勢項的異常值)，殘差就會慢慢向0收縮，最后將殘差為0對應(yīng)的原始信號點進行剔除。

1.2 形變監(jiān)測缺失數(shù)據(jù)的廣義延拓插值法

粗差剔除后，原始信號中會存在信號間斷及數(shù)據(jù)缺失，為使監(jiān)測數(shù)據(jù)連續(xù)，需要進行插值運算，本文采用廣義延拓插值法進行數(shù)據(jù)填補。廣義延拓吸取現(xiàn)有插值法和擬合法的長處與特點，采用分片光滑的做法，利用延拓域構(gòu)建單元域擬合函數(shù)，并鎖定單元域邊界節(jié)點，以逼近最好效果。

圖1 延拓域及相應(yīng)函數(shù)值Fig.1 Continuation domain and correspondingfunction value

在延拓域構(gòu)建逼近函數(shù)ye(x)：

(2)

建立廣義延拓逼近內(nèi)插模型：

(3)

式中，I(a1,a2,…,aj)為逼近函數(shù)與先驗值的誤差。

綜上，本文GNSS形變時間序列粗差數(shù)據(jù)處理流程見圖2。

圖2 數(shù)據(jù)處理流程Fig.2 Data processing flow

2 實驗及分析

2.1 模擬驗證

為驗證新算法的有效性，模擬GNSS形變的坐標時間序列，其原始表達式為：

(4)

采樣頻率為1 Hz，取4 096個歷元，考慮到實際形變監(jiān)測過程中粗差數(shù)量多、分布范圍廣，在原始信號中添加10倍標準差的粗差點，粗差點數(shù)為500個，占比為12.22%(圖3)。采用3σ法、MAD法、四分位距法和一階導(dǎo)數(shù)小波剔除法進行對比，小波基選擇為db1～db5，分別作1～8層分解并計算剔除率，仿真次數(shù)為10 000次，結(jié)果見表1和2。

圖3 原始數(shù)據(jù)及全部粗差Fig.3 Raw data and all gross errors

表1 常規(guī)粗差剔除法的剔除率

表2 一階導(dǎo)數(shù)小波剔除法在不同小波分解情況的剔除率

從表2可見，除db1外其他幾種小波使用效果基本相似，說明對粗差影響最大的是分解層數(shù)。由于增加小波分解層數(shù)會大幅增加運算量，即可以隨意選取一組小波基作一層分解，本文后續(xù)都采用基函數(shù)db4進行一層分解。

4種粗差剔除法效果對比見圖4，可以看出，一階導(dǎo)數(shù)的小波剔除率遠高于其他3種傳統(tǒng)粗差剔除算法，幾乎能探測到所有的粗差點。最后將圖4(d)剔除后的信號分別進行廣義延拓插值、分段三次樣條插值(spline)、相鄰非缺失值的線性插值(linear)、保形分段三次樣條插值(pchip)及修正Akima三次Hermite插值(makima)，并對比原始信號計算RMSE來評價插值效果，結(jié)果如表3所示。可以看出，廣義延拓插值效果優(yōu)于其余插值方法，故將其應(yīng)用于后續(xù)計算。

圖4 4種粗差剔除法剔除效果Fig.4 Four gross error elimination methods

表3 插值精度對比

2.2 工程實例測試

本文使用富陽市金鑫廣場桁架張拉測試時的數(shù)據(jù)，時間跨度為2020-10-26～27，共47 561個歷元。監(jiān)測拆除施工輔助支架桁架下沉時高度數(shù)據(jù)的變化情況，現(xiàn)場遍布高壓電線和信號基站，使監(jiān)測設(shè)備終端接收信號產(chǎn)生干擾，造成原始數(shù)據(jù)存在許多粗差點，真實波形難以觀察。由于真實信號未知，對原始含噪數(shù)據(jù)采用移動中位數(shù)進行平滑處理，得到近似的真實信號，窗口長度為200。

從圖5(a)～5(c)可以看出，傳統(tǒng)粗差算法難以準確還原信號的原有趨勢，而一階導(dǎo)數(shù)的小波剔除法保留了信號的有用成分。為達到較好的剔除效果，3種傳統(tǒng)算法每次需要剔除1個點后再重新將剩余數(shù)據(jù)進行相同的操作，4種粗差剔除算法運行時間結(jié)果見表4。

表4 算法執(zhí)行速度

由于MAD法和四分位距法需要計算中位數(shù)和百分位數(shù)，所以花費時間較多，而本文提出的新算法所用時間僅是3σ法的0.01倍。粗差剔除的正確性會影響到插值算法的效果，故將4組算法得到的剔除數(shù)據(jù)與原始數(shù)據(jù)進行精度對比，并挑選精度最高的進行廣義延拓插值，以補充剔除位置的數(shù)據(jù)。

由于真實信號采用的是原始信號平滑濾波后的數(shù)據(jù)，因此除了RMSE外，同時使用均值偏離量u作為精度評價指標，其計算公式為：

(5)

式中，s為平滑濾波后的真實信號近似，s′為經(jīng)過粗差剔除后剩余的信號，q和p分別為2個信號的總長度。

圖5 金鑫廣場桁架形變監(jiān)測數(shù)據(jù)及粗差剔除效果Fig.5 Jinxin square truss deformation monitoring data and gross error removal effect

從表5可見，前4種算法的RMSE較為接近，說明這4種算法都可以剔除影響較大的誤差。3σ法、MAD法和四分位距法存在過剔除現(xiàn)象，所以剔除后信號的均值與真實值存在較大的偏離，而一階導(dǎo)數(shù)小波剔除法能夠避免窗口的選取，經(jīng)過廣義延拓插值后與真實值僅存在約0.03 mm的偏離。

表5 精度對比

3 結(jié) 語

針對GNSS形變監(jiān)測的傳統(tǒng)粗差探測算法性能受限于數(shù)據(jù)長度的問題，提出基于小波分析的一階導(dǎo)數(shù)粗差剔除法。采用原始信號的一階導(dǎo)數(shù)來區(qū)分正常波動和粗差異常，可以有效避免窗口數(shù)據(jù)的選擇，而且只進行1次一階導(dǎo)數(shù)的剔除，大幅減少了計算量；同時，在粗差點采用廣義延拓插值補點，保證了監(jiān)測數(shù)據(jù)的連續(xù)性。經(jīng)實例驗證，本文算法的計算效率和準確度都較傳統(tǒng)算法有較大提升。