LOESS用于GNSS垂向坐標(biāo)時間序列噪聲識別與提取

陳 祥 楊志強 楊 兵 楊偉華

1 長安大學(xué)地質(zhì)工程與測繪學(xué)院，西安市雁塔路126號，710054

全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system, GNSS)作為空間大地測量的主要技術(shù)手段之一，已被廣泛應(yīng)用于地球內(nèi)部動力學(xué)機(jī)制、板塊運動和地殼形變監(jiān)測等研究[1]。目前，中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)(陸態(tài)網(wǎng)絡(luò))已建成由近300個GNSS基準(zhǔn)站和2 000多個GNSS流動站構(gòu)成的觀測網(wǎng)絡(luò)[2]。由于受數(shù)據(jù)處理策略、鐘差、電離層延遲、對流層延遲和多路徑效應(yīng)等因素的影響，解算得到的GNSS垂向坐標(biāo)時間序列存在誤差，從而影響GNSS定位的精度和可靠性，導(dǎo)致某些地球物理現(xiàn)象出現(xiàn)錯誤解釋[3]。因此，對GNSS垂向坐標(biāo)時間序列進(jìn)行降噪研究具有非常重要的科學(xué)意義和現(xiàn)實意義[4]。

目前用于GNSS坐標(biāo)時間序列分解的方法主要有小波[5]、經(jīng)驗?zāi)B(tài)分解[6]及STL(seasonal-trend decomposition procedure based on LOESS)[7]等，其中小波和經(jīng)驗?zāi)B(tài)分解方法側(cè)重于去噪分析，而STL方法主要用于季節(jié)性變化的提取，較少涉及去噪效果分析[8]。本文采用LOESS方法對GNSS垂向坐標(biāo)時間序列進(jìn)行降噪分析，并對其有效性及降噪效果進(jìn)行驗證和評價。

1 LOESS原理及降噪方法

1.1 LOESS方法基本原理

LOESS是一種基于局部回歸分析的非參數(shù)降噪方法[7]，與線性回歸、多項式回歸等方法相比，該方法未對預(yù)測點進(jìn)行限制，而是直接基于數(shù)據(jù)進(jìn)行分析。該方法需要預(yù)先設(shè)置一個窗口，根據(jù)窗口內(nèi)的點計算擬合點的擬合值，具體過程如下：

2)設(shè)(xk,yk)為窗口在最左端時的中心點，重復(fù)步驟1)，對(x1,y1)到(xk,yk)進(jìn)行擬合，此時窗口位置固定，但從點(xk+1,yk+1)開始需要以擬合點為中心移動窗口位置，并逐個計算各點的擬合值，參與擬合值計算的仍為窗口內(nèi)所有點。

3)當(dāng)窗口移動到坐標(biāo)時間序列的最右端時，與最左端的情況類似，此時窗口中心點(xm-k+1,ym-k+1)到右端點(xm,ym)(m為坐標(biāo)時間序列的總長度)均采用同一窗口內(nèi)的所有點進(jìn)行擬合。

4)通過上述步驟對坐標(biāo)時間序列中所有點進(jìn)行擬合，將擬合后的點連接成完整的LOESS回歸曲線，該曲線即為降噪后的坐標(biāo)時間序列。

該方法的關(guān)鍵在于設(shè)置窗口寬度，經(jīng)過多次實驗，本文將窗寬設(shè)置為75。關(guān)于窗口內(nèi)所有點權(quán)重的計算，可通過窗口內(nèi)各點與擬合點的距離進(jìn)行確定，該距離指X軸方向上的距離。隨著距離的增加，權(quán)重減小，因此需要采用三次權(quán)重函數(shù)進(jìn)行轉(zhuǎn)化：

(1)

此時還需要確定距離擬合點最遠(yuǎn)的點，并計算最大距離，對其他距離進(jìn)行歸一化處理，從而確定各點的權(quán)重：

(2)

式中，x為擬合點橫坐標(biāo)，xi為窗口內(nèi)某點的橫坐標(biāo)，Δ(x)為最遠(yuǎn)點與擬合點之間的距離，υi(x)為該點相對于x的權(quán)重。

當(dāng)目標(biāo)模型為非線性模型時，局部加權(quán)回歸算法相對于線性回歸算法更合適，該算法選擇窗口內(nèi)的點而不是全部點進(jìn)行回歸。加權(quán)最小二乘法的目標(biāo)函數(shù)為：

(3)

式中，J(θ)為損失函數(shù)，(xi,yi)為窗口內(nèi)某點，θ為最佳回歸系數(shù)。

1.2 技術(shù)路線

用LOESS方法對陸態(tài)網(wǎng)絡(luò)289個測站的垂向坐標(biāo)時間序列進(jìn)行降噪分析，具體流程如下：1)降噪前先剔除坐標(biāo)時間序列中的粗差；2)用LOESS方法對剔除粗差后的序列進(jìn)行降噪處理，得到降噪后的坐標(biāo)時間序列及噪聲序列；3)確定噪聲模型，并利用極大似然估計(MLE)方法解算GNSS坐標(biāo)時間序列諧波模型對應(yīng)的各項運動參數(shù)；4)檢驗降噪后噪聲序列的自相關(guān)性和降噪前后各指標(biāo)的顯著性；5)綜合數(shù)學(xué)指標(biāo)和GNSS測站運動模型參數(shù)對降噪結(jié)果進(jìn)行評價。技術(shù)路線如圖1所示。

圖1 技術(shù)路線

2 數(shù) 據(jù)

2.1 數(shù)據(jù)來源

考慮到坐標(biāo)時間序列數(shù)據(jù)缺失的問題，本文采用陸態(tài)網(wǎng)絡(luò)中289個GNSS基準(zhǔn)站原始的垂向坐標(biāo)時間序列進(jìn)行降噪分析。GNSS坐標(biāo)時間序列數(shù)據(jù)來源于中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(http:∥www.cgps.ac.cn)，數(shù)據(jù)處理方法與精度指標(biāo)可參考該平臺提供的說明文檔(ftp:∥ftp.cgps.ac.cn/doc/processing_manual.pdf)。

2.2 數(shù)據(jù)預(yù)處理

由于受各種外界因素的影響，原始數(shù)據(jù)中含有一定的粗差，在分析前需要對序列中的粗差進(jìn)行探測與剔除。本文采用IQR方法[9]消除GNSS原始坐標(biāo)時間序列中的異常值。

3 算例及結(jié)果分析

3.1 算例分析

采用§1.1的方法對選取的289個GNSS站進(jìn)行降噪處理。為驗證LOESS方法對GNSS坐標(biāo)時間序列降噪的可靠性，以XJHT站為例進(jìn)行說明，圖2為XJHT站原始坐標(biāo)時間序列、降噪后時間序列及噪聲序列。由表1可知，XJHT站垂向坐標(biāo)時間序列降噪后的信噪比達(dá)48.47 dB，降噪效果較好，標(biāo)準(zhǔn)差減小了1.11 mm，白噪聲和閃爍噪聲分別從1.81 mm、11.27 mm·a-0.25降至0.01 mm、2.40 mm·a-0.25，表明序列中包含的噪聲明顯降低，速度不確定度從0.48 mm·a-1降至0.10 mm·a-1。

圖2 XJHT站垂向GPS原始坐標(biāo)時間序列、降噪后時間序列和噪聲項

表1 XJHT站降噪前后評價指標(biāo)

3.2 自相關(guān)性檢驗

為進(jìn)一步分析降噪后噪聲序列是否存在自相關(guān)性，本文采用Durbin-Watson(DW)檢驗對降噪后的噪聲序列進(jìn)行自相關(guān)性檢驗。DW檢驗是一種檢驗序列自相關(guān)性的方法[10]，采用該方法對各測站降噪后的噪聲序列進(jìn)行檢驗可得到檢驗統(tǒng)計量。由判斷準(zhǔn)則可知，若DW值約為2，即可推斷序列完全不相關(guān)。

采用§1.1的方法對289個GNSS站的垂向坐標(biāo)時間序列進(jìn)行降噪處理后，對噪聲序列進(jìn)行DW自相關(guān)性檢驗，結(jié)果如圖3所示。結(jié)果表明，DW平均值為1.64±0.19(n=289)，其中78%介于1.50～2.10之間，表明采用LOESS方法進(jìn)行降噪處理后，各測站的噪聲序列不存在自相關(guān)性，可進(jìn)行后續(xù)的評價工作。

圖3 DW檢驗結(jié)果

3.3 評價指標(biāo)

本文采用以下4種指標(biāo)來定量評價LOESS方法對GNSS坐標(biāo)時間序列降噪的可靠性和效果：1)信噪比(signal-noise ratio, SNR)；2)標(biāo)準(zhǔn)差(standard deviation, STD)；3)白噪聲(white noise, WN)和閃爍噪聲(flicker noise, FN)；4)速度不確定度。

3.3.1 信噪比

SNR計算公式可表示為：

(4)

3.3.2 GNSS測站運動模型

研究表明[11]，白噪聲(WN)和閃爍噪聲(FN)組成的噪聲模型(WN+FN)是GNSS垂向坐標(biāo)時間序列的最優(yōu)噪聲模型，因此本文采用WN+FN噪聲模型。

作為常用的GNSS坐標(biāo)時間序列分析軟件，Hector軟件采用極大似然估計方法求解GNSS測站的運動模型參數(shù)：

(5)

3.3.3 指標(biāo)改正率

參照文獻(xiàn)[12-14]使用的指標(biāo)，本文計算評價指標(biāo)的改正率(correction rate, CR)，以進(jìn)一步分析評價LOESS方法對GNSS垂向坐標(biāo)時間序列的降噪效果：

(6)

式中，Ibefore和Iafter分別為各評價指標(biāo)在時間序列降噪前后的指標(biāo)值，CR為正表示LOESS方法能夠有效削弱GNSS垂向坐標(biāo)時間序列中包含的噪聲，CR值越大則表示降噪效果越明顯。

3.4 降噪效果評價

采用LOESS方法對各測站進(jìn)行降噪，得到降噪后的時間序列及噪聲部分，并計算降噪前后SNR、STD、噪聲和速度不確定度等4個評價指標(biāo)值，分析其變化情況。為驗證降噪前后序列的各項指標(biāo)是否存在顯著性差異，本文采用非參數(shù)Wilcoxon秩和檢驗[15]對各項指標(biāo)在降噪前后進(jìn)行顯著性檢驗，若檢驗結(jié)果中P值均小于0.05，則表明在95%的置信區(qū)間內(nèi)，指標(biāo)在降噪前后存在顯著性差異。

3.4.1 SNR與STD

由式(4)可知，原始序列的SNR均為0。圖4為各測站坐標(biāo)時間序列降噪后的信噪比，測站分布直方圖如圖5所示。由圖4和5可知，垂直分量坐標(biāo)時間序列降噪后的信噪比均較高，最高達(dá)48.47 dB，均值為8.14 dB，80%以上測站介于0～15 dB之間；部分測站的信噪比為負(fù)，其原因可能為使用LOESS方法降噪時設(shè)置的窗寬較小。結(jié)果表明，采用LOESS方法進(jìn)行降噪后信噪比總體有顯著提高。

圖4 各測站坐標(biāo)時間序列降噪后的信噪比

圖5 降噪后信噪比的測站分布

圖6為各測站垂向坐標(biāo)時間序列降噪前后的標(biāo)準(zhǔn)差變化，標(biāo)準(zhǔn)差改正率的測站分布情況如圖7所示。由圖6可知，采用LOESS方法降噪后的STD均有所減小，平均減小21.35%，最高減小48.77%，且Wilcoxon檢驗表明，降噪前后標(biāo)準(zhǔn)差具有顯著性差異(p<0.05,n=289)。由圖7可知，22%的測站改正率介于30%～50%之間，34%的測站改正率介于20%～30%之間，因此采用LOESS方法能夠顯著降低時間序列的標(biāo)準(zhǔn)差。

圖6 降噪前后各測站序列的標(biāo)準(zhǔn)差

圖7 標(biāo)準(zhǔn)差改正率的測站分布

3.4.2 運動模型參數(shù)

圖8(a)和8(b)分別為各測站垂向坐標(biāo)時間序列降噪前后白噪聲和閃爍噪聲變化情況，2種噪聲改正率的測站分布情況分別見圖9(a)和9(b)。由圖8(a)可知，降噪后序列的噪聲均明顯降低，特別是白噪聲均降低至0.03 mm以下，且Wilcoxon檢驗表明，降噪前后序列的白噪聲具有顯著性差異(p<0.05,n=289)。由圖9(a)可知，94.81%的測站白噪聲改正率在95%以上。從圖8(b)可以看出，閃爍噪聲降低至0～5 mm·a-0.25范圍內(nèi)，改正率最低為68.35%，最高為85.28%，且Wilcoxon檢驗表明，降噪前后序列的閃爍噪聲具有顯著性差異(p<0.05,n=289)。由圖9(b)可知，約34%的測站閃爍噪聲改正率在80%～86%之間，約64%的測站在70%～80%之間。結(jié)果表明，LOESS方法可有效降低時間序列中的噪聲，其中HNCS站FN由0變?yōu)?1.20 mm·a-0.25，其原因可能是該測站的噪聲并不符合WN+FN噪聲模型，而降噪會將非整數(shù)譜指數(shù)冪律噪聲PL變成FN。

圖8 降噪前后各測站坐標(biāo)時間序列的白噪聲和閃爍噪聲

圖9 白噪聲和閃爍噪聲改正率的測站分布

圖10為各測站垂向坐標(biāo)時間序列降噪前后的速度不確定度，其改正率的測站分布情況如圖11所示。由圖10可知，降噪后速度不確定度也明顯降低，平均改正率為78.68%，且Wilcoxon檢驗表明，降噪前后序列的速度不確定度具有顯著性差異(p<0.05,n=289)。由圖11可知，35%的測站改正率在80%以上，63%的測站改正率在70%～80%之間。分析表明，LOESS方法可顯著降低序列的速度不確定度。

圖10 降噪前后各測站坐標(biāo)時間序列的速度不確定度

圖11 速度不確定度改正率的測站分布

4 結(jié) 語

本文將LOESS方法引入GNSS垂向坐標(biāo)時間序列的降噪研究中，并使用數(shù)學(xué)指標(biāo)、統(tǒng)計檢驗和GNSS測站運動模型參數(shù)進(jìn)行效果分析，結(jié)果表明，LOESS方法可有效提高GNSS原始坐標(biāo)時間序列的精度。

分析表明，降噪后信噪比均值提高至8.14 dB，最高達(dá)48.47 dB；標(biāo)準(zhǔn)差平均減小1.38 mm，平均改正率為21.35%，最高達(dá)48.77%。對于降噪后的序列，白噪聲和閃爍噪聲均值分別由2.685 mm、12.415 mm·a-0.25降至0.006 mm、2.723 mm·a-0.25，平均改正率為78.63%；同時，速度不確定度得到明顯改善，平均改正率達(dá)78.68%。

本文在以下方面仍需進(jìn)一步研究：1)未將本文方法與其他降噪方法進(jìn)行對比，以分析本文方法的優(yōu)勢和不足；2)本文直接采用了普遍認(rèn)為的最佳噪聲模型WN+FN，未考慮GNSS基準(zhǔn)站復(fù)雜的噪聲特性。