LOESS用于GNSS垂向坐標時間序列噪聲識別與提取

陳 祥 楊志強 楊 兵 楊偉華

1 長安大學地質工程與測繪學院，西安市雁塔路126號，710054

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

目前用于GNSS坐標時間序列分解的方法主要有小波[5]、經驗模態分解[6]及STL(seasonal-trend decomposition procedure based on LOESS)[7]等，其中小波和經驗模態分解方法側重于去噪分析，而STL方法主要用于季節性變化的提取，較少涉及去噪效果分析[8]。本文采用LOESS方法對GNSS垂向坐標時間序列進行降噪分析，并對其有效性及降噪效果進行驗證和評價。

1 LOESS原理及降噪方法

1.1 LOESS方法基本原理

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

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

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

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

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

(1)

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

(2)

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

當目標模型為非線性模型時，局部加權回歸算法相對于線性回歸算法更合適，該算法選擇窗口內的點而不是全部點進行回歸。加權最小二乘法的目標函數為：

(3)

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

1.2 技術路線

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

圖1 技術路線

2 數 據

2.1 數據來源

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

2.2 數據預處理

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

3 算例及結果分析

3.1 算例分析

采用§1.1的方法對選取的289個GNSS站進行降噪處理。為驗證LOESS方法對GNSS坐標時間序列降噪的可靠性，以XJHT站為例進行說明，圖2為XJHT站原始坐標時間序列、降噪后時間序列及噪聲序列。由表1可知，XJHT站垂向坐標時間序列降噪后的信噪比達48.47 dB，降噪效果較好，標準差減小了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原始坐標時間序列、降噪后時間序列和噪聲項

表1 XJHT站降噪前后評價指標

3.2 自相關性檢驗

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

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

圖3 DW檢驗結果

3.3 評價指標

本文采用以下4種指標來定量評價LOESS方法對GNSS坐標時間序列降噪的可靠性和效果：1)信噪比(signal-noise ratio, SNR)；2)標準差(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垂向坐標時間序列的最優噪聲模型，因此本文采用WN+FN噪聲模型。

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

(5)

3.3.3 指標改正率

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

(6)

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

3.4 降噪效果評價

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

3.4.1 SNR與STD

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

圖4 各測站坐標時間序列降噪后的信噪比

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

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

圖6 降噪前后各測站序列的標準差

圖7 標準差改正率的測站分布

3.4.2 運動模型參數

圖8(a)和8(b)分別為各測站垂向坐標時間序列降噪前后白噪聲和閃爍噪聲變化情況，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范圍內，改正率最低為68.35%，最高為85.28%，且Wilcoxon檢驗表明，降噪前后序列的閃爍噪聲具有顯著性差異(p<0.05,n=289)。由圖9(b)可知，約34%的測站閃爍噪聲改正率在80%～86%之間，約64%的測站在70%～80%之間。結果表明，LOESS方法可有效降低時間序列中的噪聲，其中HNCS站FN由0變為 1.20 mm·a-0.25，其原因可能是該測站的噪聲并不符合WN+FN噪聲模型，而降噪會將非整數譜指數冪律噪聲PL變成FN。

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

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

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

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

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

4 結 語

本文將LOESS方法引入GNSS垂向坐標時間序列的降噪研究中，并使用數學指標、統計檢驗和GNSS測站運動模型參數進行效果分析，結果表明，LOESS方法可有效提高GNSS原始坐標時間序列的精度。

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

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