基于改進的經驗模態分解和SSA聯合算法的GPS坐標時間序列重構

張祎磊 邊家文 丁開華 冉佳諾 劉文平

1 中國地質大學(武漢)數學與物理學院，武漢市魯磨路388號，430074 2 中國地質大學(武漢)地理與信息工程學院，武漢市魯磨路388號，430074 3 湖北經濟學院信息管理與統計學院，武漢市楊橋湖大道8號，430205

GPS坐標時間序列主要包括長期趨勢項、周期項以及噪聲項，對GPS坐標時間序列進行精確重構，不僅可以提取該站點的周期性信號，還可以準確估計噪聲振幅、速度不確定度等參數。近年來，國內外學者提出多種方法估計時變振幅的周期信號，如奇異譜分析法(singular spectrum analysis, SSA)[1]、小波分解法(wavelet decomposition, WD)[2]、滑動最小二乘法(moving ordinary least squares, MOLS)[3]以及經驗模態分解法(empirical mode decomposition, EMD)[4]等。其中，EMD法本身存在模態混疊，且在提取閃爍噪聲環境下的周期信號時存在吸收噪聲的現象，會導致重構精度降低。

針對EMD算法本身的局限性，有學者提出自適應噪聲完備集成經驗模態分解(CEEMD with adaptive noise, CEEMDAN)[5]和改進的自適應噪聲完備集成經驗模態分解(improved CEEMDAN, ICEEMDAN)[6]等算法。ICEEMDAN方法由CEEMDAN方法改進而來，采用局部均值的估計值替換模態估計值以減少冗余模態的影響。k階模態由信號的局部均值來提取而不直接使用白噪聲，可降低分解后模態的噪聲殘余，基本可解決模態混疊問題，可以穩定提取低頻信息。然而實際的GPS坐標時間序列中往往存在振幅較小的弱周期信號[7]，這些弱周期信號對應的奇異值與噪聲的奇異值往往較為接近，在提取時容易被噪聲掩蓋而難以順利提取。基于上述分析，本文首先利用ICEEMDAN方法提取GPS坐標時間序列中不同時間尺度下的趨勢以及周期信號，然后將ICEEMDAN方法分解得到的各個本征模態分量排成Hankel矩陣并利用SSA對分解得到的每個本征模態分量進行重構，最后獲得整個GPS坐標時間序列的重構結果。模擬結果表明，本文提出的ICEEMDAN-SSA聯合算法在重構精度上優于多個現有方法，對于閃爍噪聲具有更好的壓制效果，可以更好地分離出趨勢和周期信號。

1 具有時變振幅的GPS坐標時間序列模型

GPS坐標時間序列一般由長期趨勢、周期項以及噪聲項組成，一般建模如下：

y(t)=h0+vt+s(t)+ε(t)，t=1,2,…,N

(1)

式中，h0為初始位置，v為速度，s(t)為周期項，ε(t)為噪聲項。其中，周期項s(t)通常認為由年周期信號和半年周期信號構成。研究表明[7-8]，對流層和地磁活動中存在準兩年振蕩，可能會對日長和極移產生潛在影響，從而影響GPS觀測結果。準兩年信號會使速度估計產生偏差，因此有必要在周期信號模型中考慮準兩年周期信號。由于準兩年周期與兩年周期的頻率信息較為接近，本文模型由半年、年以及較小振幅的兩年周期信號構成，即

bi(t)sin(2πfit)

(2)

式中，ai(t)、bi(t)為時變振幅項，數學形式上表現為振幅項隨著時間的變化而波動，時變振幅選取以下形式[1]：

ai(t)=ai+epsin(t)，bi(t)=bi+eqsin(t)

(3)

式中，ai、bi(i=1,2,3)為周期項振幅的常數項，p、q為周期波動參數。

2 基于ICEEMDAN-SSA的GPS坐標時間序列分解和重構

基于ICEEMDAN-SSA算法的信號分解和重構具體步驟如下：

1)輸入GPS坐標時間序列x，通過EMD法進行第i次迭代(i=1,2,…,I)，迭代公式如下：

x(i)=x+β0E1(w(i))

(4)

式中，E1為對輸入信號進行EMD分解并獲得的第一個模態分量；i為迭代次數；w(i)為標準正態分布的白噪聲；β0為噪聲系數，即加入白噪聲的標準差和輸入信號的標準差的比值。通過計算其局部均值獲得第一個殘差序列r1：

(5)

式中，M為產生信號局部均值的算子。第一個模態計算公式如下：

IMF1=x-r1

(6)

2)第二個殘差序列r2以及第二個模態IMF2計算公式如下：

(7)

IMF2=r1-r2

(8)

3)循環計算IMFk和rk，公式如下：

(9)

IMFk=rk-1-rk

(10)

式中，k=3,4,5,…,Q,Q為滿足分解結束條件時共產生的模態分量個數。循環計算式(9)、式(10)直至殘差序列極值點小于2，分解步驟結束。

4)得到Q個不同時間尺度下的模態{IMF1,IMF2,…,IMFQ}，通過計算各模態分量與原GPS信號的相關系數和方差貢獻率篩選出噪聲模態和信號模態，此時有：

(11)

式中，q為噪聲模態的數量。

5)將各信號模態分別排成Hankel矩陣Ds(L,K)形式，其中s為第s個信號模態分量(s=q+1,q+2,…,Q)，L=N-K+1≤N/2：

(12)

計算其協方差矩陣Cj：

(13)

(14)

6)計算重建成分：

(15)

(16)

3 模擬實驗分析

3.1 模擬實驗設置

在模擬實驗與實際站點實驗中，ICEEMDAN-SSA方法的噪聲系數βj=5，迭代次數i=1 000，滑動窗口長度L=5。對于SSA方法，窗口長度選為3 a[1]。WD方法采用Meyer小波濾波器，將信號分解后得到的第7道和第8道信號作為提取的周期信號[2]。對于MOLS方法，在模擬實驗中對長度為6 000 d的GPS坐標時間序列進行分段處理，每3 a為一小段，相鄰兩段重疊時間段為1 a，取均值作為重疊部分進行周期信號重構。

本文選取4個評價指標來衡量各個算法對GPS坐標時間序列周期信號的重構精度，并對比各算法來檢驗ICEEMDAN-SSA方法的有效性。評價指標包括：1)Misfit：估計信號與不帶噪的真實周期信號的標準差；2)譜指數κ：噪聲在功率譜密度對數形式下的斜率[9]；3)殘差振幅Apl：殘差估計振幅；4)速度不確定度：由譜指數κ和殘差振幅Apl計算得出，具體公式見文獻[10]。

3.2 實驗結果分析

圖1 ICEEMDAN方法提取的周期信號Fig.1 Periodic signals extracted by theICEEMDAN method

由圖1可知，在ICEEMDAN方法的處理階段，已初步分解出含有半年、一年、兩年周期信號的模態分量。但由于GPS信號的有色噪聲規模較大，初步分解后得到的模態分量仍有部分噪聲殘余。

由圖2可知，在SSA方法的處理階段，通過對每一個周期信號選取貢獻率更高的特征值來重構模態分量，可以對各模態分量進一步去噪，提取出精度更高的各周期分量以及趨勢信號。將各個重構的趨勢、周期模態分量累加即可得到重構后的時間序列。將ICEEMDAN-SSA與SSA、WD、MOLS、ICEEMDAN方法進行比較，對比結果如圖3、圖4、表1所示。

圖2 基于ICEEMDAN和ICEEMDAN-SSA提取周期信息與趨勢信息Fig.2 Periodic information and trend information extracted by ICEEMDAN and ICEEMDAN-SSA

圖4 閃爍噪聲振幅為時各個方法得到的PSDFig.4 PSD obtained by various estimation methodswith noise amplitude of 10

表1 閃爍噪聲振幅為時各方法的指標估計值

由圖3可知，SSA、MOLS、ICEEMDAN-SSA方法對周期信號均有較好的擬合精度，其中ICEEMDAN-SSA方法的Misfit指標最小，說明該方法的擬合效果最好。結合圖4各方法殘差序列的PSD可知，在噪聲振幅較大時，WD吸收噪聲情況嚴重。ICEEMDAN方法在2 cpy處與噪聲PSD較為相近，在0.5 cpy與1 cpy處均吸收大量噪聲，說明ICEEMDAN方法對閃爍噪聲去噪效果有限。SSA與MOLS方法效果相近，這兩種方法均在0.5 cpy處與噪聲PSD較為吻合，在1 cpy以及2 cpy處對噪聲PSD曲線擬合較差，說明兩者在估計年周期和半年周期信號時存在吸收噪聲的情況。ICEEMDAN-SSA方法在2 cpy處與噪聲曲線幾乎完全吻合，但在0.5 cpy以及1 cpy 處與噪聲曲線存在少量偏差，說明該方法對半年周期信號的估計精度較高，在估計年、兩年周期信號時會吸收少量噪聲。

4 實際數據分析

選取昆明(KUNM)GPS觀測站的帶趨勢數據進行分析，同時也將ICEEMDAN-SSA聯合算法與其他幾種經典方法進行對比來檢驗其提取時變振幅周期信號的有效性。以下真實站點數據可通過IGS全球數據中心SOPAC網站下載。本文選取的KUNM站點信息如表2所示。

表2 GPS站點信息

為盡可能地減少系統誤差，只選取各站點最近3 000 d的數據進行實驗。由于GPS坐標時間序列存在數據缺失以及野值干擾的問題，實驗開始前需要對數據進行偏移量檢測和缺失數據插值等預處理，本文根據3σ準則去除野值[11]，然后通過KKF插值方法進行插值[12]。

對于插值后的時間序列，使用ICEEMDAN-SSA方法和其他傳統方法分別對KUNM站點N、E、U方向的時間序列進行重構，結果如圖5所示。

與模擬實驗不同的是，在處理實際數據時無法得到真實無噪的干凈信號，因此Misfit指標不予計算，在評價各方法時僅考慮速度不確定度、譜指數κ以及殘差振幅Apl(表3)。由表可知，相比于MOLS方法，ICEEMDAN-SSA方法的重構曲線更加平滑，與其他傳統方法在擬合效果上更為接近；ICEEMDAN-SSA方法的殘差振幅與SSA方法較為接近；與ICEEMDAN方法相比，ICEEMDAN-SSA方法在重構穩定性方面有所提高，不會出現噪聲分離失敗的現象。

本文通過計算各種方法得到的重構信號之間的相關系數來了解各個估計方法的有效性(表4)。由表可知，SSA方法和ICEEMDAN方法與ICEEMDAN-SSA聯合重構算法的重構結果更為接近。除此之外，本文涉及的各算法重構結果之間的相關系數在N方向和E方向上均為0.99左右，U方向均不小于0.91，說明本文涉及的所有方法均對實際GPS觀測信號具備有效性。

圖5 KUNM站各方向時間序列重構結果Fig.5 Reconstruction results of time series in each direction of KUNM station

表3 通過不同方法得到KUNM站點各指標值

表4 ICEEMDAN-SSA在KUNM站N、E、U方向的重建結果與SSA、WD、MOLS、ICEEMDAN重建結果之間的相關系數

5 結 論

本文利用ICEEMDAN算法多尺度分解的優勢，結合SSA提出基于ICEEMDAN-SSA的GPS坐標時間序列聯合估計方法。該方法利用ICEEMDAN在提取周期、趨勢時可以將信號和噪聲分離的優勢，同時可保留不同振幅大小的周期信號；對于單個分離信號，利用SSA精確提取時間序列中的主要特征，可以更好地提升GPS坐標時間序列的重構精度。通過對比模擬數據和實際觀測數據，將ICEEMDAN-SSA方法與SSA、WD、MOLS三種傳統方法進行對比。

模擬實驗表明，本文提出的ICEEMDAN-SSA方法相比上述方法對GPS坐標時間序列的重構效果相對更優。實驗結果表明，ICEEMDAN-SSA方法的重構信號與各方法的重構信號均具備較高的相關度，這也表明本文方法具備有效性。綜上所述，ICEEMDAN-SSA方法對于GPS坐標時間序列的周期、趨勢信號具有更好的估計精度與重構效果。