基于改進(jìn)的經(jīng)驗(yàn)?zāi)B(tài)分解和SSA聯(lián)合算法的GPS坐標(biāo)時(shí)間序列重構(gòu)

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

1 中國(guó)地質(zhì)大學(xué)(武漢)數(shù)學(xué)與物理學(xué)院，武漢市魯磨路388號(hào)，430074 2 中國(guó)地質(zhì)大學(xué)(武漢)地理與信息工程學(xué)院，武漢市魯磨路388號(hào)，430074 3 湖北經(jīng)濟(jì)學(xué)院信息管理與統(tǒng)計(jì)學(xué)院，武漢市楊橋湖大道8號(hào)，430205

GPS坐標(biāo)時(shí)間序列主要包括長(zhǎng)期趨勢(shì)項(xiàng)、周期項(xiàng)以及噪聲項(xiàng)，對(duì)GPS坐標(biāo)時(shí)間序列進(jìn)行精確重構(gòu)，不僅可以提取該站點(diǎn)的周期性信號(hào)，還可以準(zhǔn)確估計(jì)噪聲振幅、速度不確定度等參數(shù)。近年來(lái)，國(guó)內(nèi)外學(xué)者提出多種方法估計(jì)時(shí)變振幅的周期信號(hào)，如奇異譜分析法(singular spectrum analysis, SSA)[1]、小波分解法(wavelet decomposition, WD)[2]、滑動(dòng)最小二乘法(moving ordinary least squares, MOLS)[3]以及經(jīng)驗(yàn)?zāi)B(tài)分解法(empirical mode decomposition, EMD)[4]等。其中，EMD法本身存在模態(tài)混疊，且在提取閃爍噪聲環(huán)境下的周期信號(hào)時(shí)存在吸收噪聲的現(xiàn)象，會(huì)導(dǎo)致重構(gòu)精度降低。

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

1 具有時(shí)變振幅的GPS坐標(biāo)時(shí)間序列模型

GPS坐標(biāo)時(shí)間序列一般由長(zhǎng)期趨勢(shì)、周期項(xiàng)以及噪聲項(xiàng)組成，一般建模如下：

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

(1)

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

bi(t)sin(2πfit)

(2)

式中，ai(t)、bi(t)為時(shí)變振幅項(xiàng)，數(shù)學(xué)形式上表現(xiàn)為振幅項(xiàng)隨著時(shí)間的變化而波動(dòng)，時(shí)變振幅選取以下形式[1]：

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

(3)

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

2 基于ICEEMDAN-SSA的GPS坐標(biāo)時(shí)間序列分解和重構(gòu)

基于ICEEMDAN-SSA算法的信號(hào)分解和重構(gòu)具體步驟如下：

1)輸入GPS坐標(biāo)時(shí)間序列x，通過(guò)EMD法進(jìn)行第i次迭代(i=1,2,…,I)，迭代公式如下：

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

(4)

式中，E1為對(duì)輸入信號(hào)進(jìn)行EMD分解并獲得的第一個(gè)模態(tài)分量；i為迭代次數(shù)；w(i)為標(biāo)準(zhǔn)正態(tài)分布的白噪聲；β0為噪聲系數(shù)，即加入白噪聲的標(biāo)準(zhǔn)差和輸入信號(hào)的標(biāo)準(zhǔn)差的比值。通過(guò)計(jì)算其局部均值獲得第一個(gè)殘差序列r1：

(5)

式中，M為產(chǎn)生信號(hào)局部均值的算子。第一個(gè)模態(tài)計(jì)算公式如下：

IMF1=x-r1

(6)

2)第二個(gè)殘差序列r2以及第二個(gè)模態(tài)IMF2計(jì)算公式如下：

(7)

IMF2=r1-r2

(8)

3)循環(huán)計(jì)算IMFk和rk，公式如下：

(9)

IMFk=rk-1-rk

(10)

式中，k=3,4,5,…,Q,Q為滿足分解結(jié)束條件時(shí)共產(chǎn)生的模態(tài)分量個(gè)數(shù)。循環(huán)計(jì)算式(9)、式(10)直至殘差序列極值點(diǎn)小于2，分解步驟結(jié)束。

4)得到Q個(gè)不同時(shí)間尺度下的模態(tài){IMF1,IMF2,…,IMFQ}，通過(guò)計(jì)算各模態(tài)分量與原GPS信號(hào)的相關(guān)系數(shù)和方差貢獻(xiàn)率篩選出噪聲模態(tài)和信號(hào)模態(tài)，此時(shí)有：

(11)

式中，q為噪聲模態(tài)的數(shù)量。

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

(12)

計(jì)算其協(xié)方差矩陣Cj：

(13)

(14)

6)計(jì)算重建成分：

(15)

(16)

3 模擬實(shí)驗(yàn)分析

3.1 模擬實(shí)驗(yàn)設(shè)置

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

本文選取4個(gè)評(píng)價(jià)指標(biāo)來(lái)衡量各個(gè)算法對(duì)GPS坐標(biāo)時(shí)間序列周期信號(hào)的重構(gòu)精度，并對(duì)比各算法來(lái)檢驗(yàn)ICEEMDAN-SSA方法的有效性。評(píng)價(jià)指標(biāo)包括：1)Misfit：估計(jì)信號(hào)與不帶噪的真實(shí)周期信號(hào)的標(biāo)準(zhǔn)差；2)譜指數(shù)κ：噪聲在功率譜密度對(duì)數(shù)形式下的斜率[9]；3)殘差振幅Apl：殘差估計(jì)振幅；4)速度不確定度：由譜指數(shù)κ和殘差振幅Apl計(jì)算得出，具體公式見文獻(xiàn)[10]。

3.2 實(shí)驗(yàn)結(jié)果分析

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

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

由圖2可知，在SSA方法的處理階段，通過(guò)對(duì)每一個(gè)周期信號(hào)選取貢獻(xiàn)率更高的特征值來(lái)重構(gòu)模態(tài)分量，可以對(duì)各模態(tài)分量進(jìn)一步去噪，提取出精度更高的各周期分量以及趨勢(shì)信號(hào)。將各個(gè)重構(gòu)的趨勢(shì)、周期模態(tài)分量累加即可得到重構(gòu)后的時(shí)間序列。將ICEEMDAN-SSA與SSA、WD、MOLS、ICEEMDAN方法進(jìn)行比較，對(duì)比結(jié)果如圖3、圖4、表1所示。

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

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

表1 閃爍噪聲振幅為時(shí)各方法的指標(biāo)估計(jì)值

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

4 實(shí)際數(shù)據(jù)分析

選取昆明(KUNM)GPS觀測(cè)站的帶趨勢(shì)數(shù)據(jù)進(jìn)行分析，同時(shí)也將ICEEMDAN-SSA聯(lián)合算法與其他幾種經(jīng)典方法進(jìn)行對(duì)比來(lái)檢驗(yàn)其提取時(shí)變振幅周期信號(hào)的有效性。以下真實(shí)站點(diǎn)數(shù)據(jù)可通過(guò)IGS全球數(shù)據(jù)中心SOPAC網(wǎng)站下載。本文選取的KUNM站點(diǎn)信息如表2所示。

表2 GPS站點(diǎn)信息

為盡可能地減少系統(tǒng)誤差，只選取各站點(diǎn)最近3 000 d的數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。由于GPS坐標(biāo)時(shí)間序列存在數(shù)據(jù)缺失以及野值干擾的問(wèn)題，實(shí)驗(yàn)開始前需要對(duì)數(shù)據(jù)進(jìn)行偏移量檢測(cè)和缺失數(shù)據(jù)插值等預(yù)處理，本文根據(jù)3σ準(zhǔn)則去除野值[11]，然后通過(guò)KKF插值方法進(jìn)行插值[12]。

對(duì)于插值后的時(shí)間序列，使用ICEEMDAN-SSA方法和其他傳統(tǒng)方法分別對(duì)KUNM站點(diǎn)N、E、U方向的時(shí)間序列進(jìn)行重構(gòu)，結(jié)果如圖5所示。

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

本文通過(guò)計(jì)算各種方法得到的重構(gòu)信號(hào)之間的相關(guān)系數(shù)來(lái)了解各個(gè)估計(jì)方法的有效性(表4)。由表可知，SSA方法和ICEEMDAN方法與ICEEMDAN-SSA聯(lián)合重構(gòu)算法的重構(gòu)結(jié)果更為接近。除此之外，本文涉及的各算法重構(gòu)結(jié)果之間的相關(guān)系數(shù)在N方向和E方向上均為0.99左右，U方向均不小于0.91，說(shuō)明本文涉及的所有方法均對(duì)實(shí)際GPS觀測(cè)信號(hào)具備有效性。

圖5 KUNM站各方向時(shí)間序列重構(gòu)結(jié)果Fig.5 Reconstruction results of time series in each direction of KUNM station

表3 通過(guò)不同方法得到KUNM站點(diǎn)各指標(biāo)值

表4 ICEEMDAN-SSA在KUNM站N、E、U方向的重建結(jié)果與SSA、WD、MOLS、ICEEMDAN重建結(jié)果之間的相關(guān)系數(shù)

5 結(jié) 論

本文利用ICEEMDAN算法多尺度分解的優(yōu)勢(shì)，結(jié)合SSA提出基于ICEEMDAN-SSA的GPS坐標(biāo)時(shí)間序列聯(lián)合估計(jì)方法。該方法利用ICEEMDAN在提取周期、趨勢(shì)時(shí)可以將信號(hào)和噪聲分離的優(yōu)勢(shì)，同時(shí)可保留不同振幅大小的周期信號(hào)；對(duì)于單個(gè)分離信號(hào)，利用SSA精確提取時(shí)間序列中的主要特征，可以更好地提升GPS坐標(biāo)時(shí)間序列的重構(gòu)精度。通過(guò)對(duì)比模擬數(shù)據(jù)和實(shí)際觀測(cè)數(shù)據(jù)，將ICEEMDAN-SSA方法與SSA、WD、MOLS三種傳統(tǒng)方法進(jìn)行對(duì)比。

模擬實(shí)驗(yàn)表明，本文提出的ICEEMDAN-SSA方法相比上述方法對(duì)GPS坐標(biāo)時(shí)間序列的重構(gòu)效果相對(duì)更優(yōu)。實(shí)驗(yàn)結(jié)果表明，ICEEMDAN-SSA方法的重構(gòu)信號(hào)與各方法的重構(gòu)信號(hào)均具備較高的相關(guān)度，這也表明本文方法具備有效性。綜上所述，ICEEMDAN-SSA方法對(duì)于GPS坐標(biāo)時(shí)間序列的周期、趨勢(shì)信號(hào)具有更好的估計(jì)精度與重構(gòu)效果。