周 侗，邊家文*，丁開(kāi)華，陳保周，伍浩琛

（1.中國(guó)地質(zhì)大學(xué) 數(shù)學(xué)與物理學(xué)院，湖北 武漢 430074；2.中國(guó)地質(zhì)大學(xué) 地理與信息工程學(xué)院，湖北 武漢 430074）

GPS監(jiān)測(cè)技術(shù)已廣泛應(yīng)用于地殼形變、板塊運(yùn)動(dòng)監(jiān)測(cè)，地震、海嘯預(yù)警，參考框架建立與維持等領(lǐng)域，因此精確地分析和估計(jì)GPS坐標(biāo)時(shí)間序列具有重要意義。在精確分析和重構(gòu)GPS坐標(biāo)時(shí)間序列之前，需對(duì)GPS坐標(biāo)時(shí)間序列的噪聲類型進(jìn)行分析，國(guó)內(nèi)外學(xué)者對(duì)此進(jìn)行了較為全面的研究，如MAO A L[1]等分析了全球23個(gè)GPS站點(diǎn)長(zhǎng)達(dá)3 a的時(shí)間序列，結(jié)果表明噪聲的最佳模型是“閃爍噪聲+白噪聲”；李昭[2]等分析了中國(guó)區(qū)域GPS坐標(biāo)時(shí)間序列的噪聲類型，結(jié)果表明測(cè)站最佳噪聲模型組合為“閃爍噪聲+白噪聲”和“帶通冪律噪聲+白噪聲”；在GPS坐標(biāo)時(shí)間序列重構(gòu)方面，WANG X M[3]和CHEN Q[4]等提出了利用奇異譜分析法（SSA）提取序列周期的方法，SSA方法可以較好地從含噪聲的時(shí)間序列中提取信息，但在分離噪聲時(shí)會(huì)吸收部分有色噪聲；張雙成[5]等利用經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）方法對(duì)GPS坐標(biāo)時(shí)間序列進(jìn)行了降噪和重構(gòu)，EMD方法能合理地分離序列中的信號(hào)和噪聲，有效削弱噪聲對(duì)GPS坐標(biāo)時(shí)間序列的影響，進(jìn)一步提高GPS坐標(biāo)時(shí)間序列的精度，但在噪聲較大時(shí)，會(huì)出現(xiàn)模態(tài)混疊現(xiàn)象，導(dǎo)致重構(gòu)效果較差；董偉[6]采用小波分解（WD）方法對(duì)GPS坐標(biāo)時(shí)間序列進(jìn)行了多分解、去噪和重構(gòu)，WD方法可為任何調(diào)制信號(hào)建模時(shí)變多樣的季節(jié)曲線，可在限制的頻率之間模擬信號(hào)和噪聲，從而確定分辨率，但不進(jìn)行信號(hào)和噪聲之間的分離；LIU N[7]和Didova O[8]等利用卡爾曼濾波方法（KF）對(duì)缺失的GPS坐標(biāo)時(shí)間序列進(jìn)行重構(gòu)，KF方法的最大缺點(diǎn)是必須用到無(wú)限過(guò)去的數(shù)據(jù)。

近年來(lái)，Anna等提出利用滑動(dòng)普通最小二乘法（MOLS）來(lái)重構(gòu)信號(hào)。MOLS方法通過(guò)選取合適的窗口，滑動(dòng)利用最小二乘得到每段的估計(jì)，最終得到整段GPS坐標(biāo)時(shí)間序列的估計(jì)，速度較快，但由于忽略了噪聲特性，導(dǎo)致精度較低，且會(huì)產(chǎn)生過(guò)擬合現(xiàn)象[9]。本文通過(guò)L2優(yōu)化方法來(lái)防止過(guò)擬合，由于在整段時(shí)變周期的GPS坐標(biāo)時(shí)間序列上采用L2優(yōu)化方法會(huì)導(dǎo)致較大的誤差，因此將信號(hào)分段，逐年滑動(dòng)，在每段上使用L2優(yōu)化方法，即滑動(dòng)L2優(yōu)化估計(jì)方法（ML2），并在每段采用交替迭代乘子法（ADMM）求解。ADMM算法具有分布式運(yùn)算、并行處理的特點(diǎn)，且在迭代次數(shù)較少時(shí)能得到合理的精度[10]，使GPS坐標(biāo)時(shí)間序列的重構(gòu)精度更高。

1 模型與重構(gòu)算法

1.1 GPS坐標(biāo)時(shí)間序列的模型

GPS坐標(biāo)時(shí)間序列的一般模型為：

式中，n=1,2,…,N，N為信號(hào)長(zhǎng)度；y(tn)為信號(hào)；x0和vx分別為初始位置和速度；ai和bi分別是角速度為ωi的周期信號(hào)的正弦和余弦項(xiàng)的常數(shù)；εt為噪聲，本文采用“閃爍噪聲+白噪聲”的噪聲模型；t0為參考時(shí)間；tn為時(shí)間。

在實(shí)際信號(hào)處理中周期信號(hào)的振幅ai和bi往往會(huì)隨著時(shí)間發(fā)生變化[11]，因此本文利用的時(shí)變周期的GPS坐標(biāo)時(shí)間序列模型為：

式中，ai(t)和bi(t)為時(shí)變振幅，假設(shè)它們由其均值和隨機(jī)擾動(dòng)組成，即ai(t)=ai+μi(t)，bi(t)=bi+vi(t)，μi(t)和vi(t)為隨機(jī)擾動(dòng)，因此需要估計(jì)時(shí)變振幅的均值。

記：則式（2）可寫為：

式（2）可轉(zhuǎn)化為式（1），相應(yīng)的矩陣表示形式為：

式中，A1X1為趨勢(shì)項(xiàng)；A2X2為周期項(xiàng)；ε為噪聲項(xiàng)。

1.2 基于L2優(yōu)化準(zhǔn)則的GPS坐標(biāo)時(shí)間序列重構(gòu)

1.2.1 優(yōu)化模型

為了估計(jì)周期項(xiàng)和趨勢(shì)項(xiàng)，建立的優(yōu)化模型為：

式中，r=Y-A1X1-A2X2；α和β為正則項(xiàng)的參數(shù)。本文采用L2范數(shù)約束參數(shù)X1和X2，防止參數(shù)X1和X2過(guò)擬合。

1.2.2 ADMM算法

ADMM算法既有乘子法的強(qiáng)收斂性質(zhì)又有對(duì)偶上升法的分解性，因此該方法比較適合求解大規(guī)模的優(yōu)化問(wèn)題。目前，該方法已廣泛應(yīng)用于壓縮感知、圖像處理、機(jī)器學(xué)習(xí)等諸多領(lǐng)域。ADMM算法的優(yōu)點(diǎn)是僅有一個(gè)參數(shù)ρ，且在一定條件下，其對(duì)任意參數(shù)均收斂。對(duì)于優(yōu)化問(wèn)題式（8），可使用梯度下降法求解，但該方法對(duì)步長(zhǎng)的選擇非常敏感，選擇不合適的步長(zhǎng)甚至?xí)?dǎo)致算法不收斂；而ADMM算法對(duì)步長(zhǎng)具有較強(qiáng)的魯棒性，在一定條件下，能保證對(duì)所有的步長(zhǎng)都收斂[12]。因此，本文采用ADMM算法求解優(yōu)化問(wèn)題式（8）。

構(gòu)造優(yōu)化問(wèn)題式（8）的增廣拉格朗日函數(shù)，則有：

式中，Z為對(duì)偶變量并具有適當(dāng)?shù)目s放尺寸；ρ為懲罰參數(shù)[13]。

按照順序更新變量r、X1和X2，再更新對(duì)偶變量Z，即

對(duì)式（10）的優(yōu)化目標(biāo)函數(shù)進(jìn)行簡(jiǎn)化，可得：

對(duì)rk+1求導(dǎo)，令結(jié)果為0，則有：

可得rk+1的迭代公式為：

對(duì)式（11）的優(yōu)化目標(biāo)函數(shù)進(jìn)行簡(jiǎn)化，則有：

式中，I2和I4分別為2階單位矩陣和4階單位矩陣。

ADMM算法的收斂性由其殘差來(lái)刻畫，其停止準(zhǔn)則為：

式中，rpri為原始誤差；sdual為對(duì)偶誤差；εpri和εdual為給定的閾值。

1.2.3 算法評(píng)價(jià)指標(biāo)和計(jì)算公式

本文采用4個(gè)評(píng)價(jià)指標(biāo)全面衡量重構(gòu)信號(hào)的效果，通過(guò)Misfit衡量重構(gòu)信號(hào)和真實(shí)信號(hào)的誤差，噪聲譜指數(shù)和噪聲振幅衡量殘差和真實(shí)噪聲的接近程度，速度不確定度衡量殘差對(duì)速度估計(jì)的影響。

1）Misfit為真實(shí)信號(hào)和重構(gòu)信號(hào)的標(biāo)準(zhǔn)差。2）噪聲譜指數(shù)（κ）。功率譜的計(jì)算公式為[14]：

兩邊取對(duì)數(shù)，則有：

式中，P(f)為功率譜密度（本文采用經(jīng)典周期圖法計(jì)算）；f為頻率；P0和f0為常量。本文采用最小二乘法估計(jì)噪聲譜指數(shù)。

3）噪聲振幅（σ）。本文采用的噪聲模型為：

式中，α為白噪聲；e為白噪聲的振幅；βκ為譜指數(shù)為κ的有色噪聲（κ=-1為閃爍噪聲）；σ為譜指數(shù)為κ的有色噪聲的振幅[15]。

ε的協(xié)方差矩陣C的計(jì)算公式為[16]：

式中，I為單位矩陣；Jκ為譜指數(shù)為κ的噪聲的單位協(xié)方差矩陣。

4）速度不確定度（Uv）[17-19]。

式中，ΔT為時(shí)間間隔；Γ為Gamma函數(shù)。

1.2.4 算法流程

通過(guò)上述分析，基于L2優(yōu)化準(zhǔn)則的GPS坐標(biāo)時(shí)間序列重構(gòu)方法的具體步驟為：①根據(jù)式（2）產(chǎn)生時(shí)變周期的GPS坐標(biāo)時(shí)間序列；②根據(jù)式（3）建立矩陣A1和A2；③根據(jù)信號(hào)的長(zhǎng)度N，選取合適的分段長(zhǎng)度W，并逐年滑動(dòng)得到M段信號(hào)（M=([N/365]-W)+2）；④針對(duì)每段信號(hào)建立式（8）的優(yōu)化模型，再按照迭代步驟估計(jì)參數(shù)X1，X2；⑤初始化α、β、ρ、r、X1、X2、Z、εpri、εdual、rpri、sdual；⑥利用式（16）得到r第k次的更新rik；⑦利用式（19）得到X1的第k次更新Xk1,i；⑧利用式（20）得到X2的第k次更新Xk2,i；⑨利用式（21）得到Z的第k次更新Zik；⑩根據(jù)式（22）和式（23）分別計(jì)算rpri和sdual；?若rpri和sdual滿足式（22）和式（23）則終止循環(huán)，否則繼續(xù)重復(fù)步驟⑥～⑨；?重復(fù)步驟④～?，直到得到每段時(shí)間序列的估計(jì)；?依據(jù)X1和X2逐年恢復(fù)信號(hào)（重復(fù)的年份取平均），并得到殘差；?根據(jù)式（25），利用最小二乘法得到殘差的譜指數(shù)估計(jì)；?根據(jù)式（26）～（30），利用最大似然法得到殘差的振幅估計(jì)；?計(jì)算重構(gòu)信號(hào)的Misfit，并根據(jù)式（31）計(jì)算速度不確定度。

2 模擬實(shí)驗(yàn)

本文采用兩組模擬數(shù)據(jù)來(lái)檢驗(yàn)本文提出的算法性能，并與MOLS、SSA和WD方法進(jìn)行對(duì)比。實(shí)驗(yàn)中MOLS每段時(shí)間長(zhǎng)度為1 100 d，滑動(dòng)間隔為1 a，SSA的窗口大小為4 a（窗口大小為1 460 d），WD方法采用的是Meyer小波。本文利用Misfit、噪聲譜指數(shù)、噪聲振幅和速度不確定度來(lái)衡量重構(gòu)信號(hào)的精度。

模擬數(shù)據(jù)的參數(shù)為：t0=0，x0=10 mm，vx=5 mm/a，a1,b1～N(3,1)，a2,b2～N(2,0.5)，ω1=2π/365，ω2=4π/365。為了驗(yàn)證算法在有色噪聲干擾情形下的有效性，本文模擬實(shí)驗(yàn)的噪聲取為閃爍噪聲。閃爍噪聲較小時(shí)，閃爍噪聲振幅取為1 mm/a-1/4；閃爍噪聲較大時(shí)，閃爍噪聲振幅取為10 mm/a-1/4。實(shí)驗(yàn)結(jié)果如表1、2所示。

表1 噪聲振幅為1 mm/a-1/4的實(shí)驗(yàn)結(jié)果

表2 噪聲振幅為10 mm/a-1/4的實(shí)驗(yàn)結(jié)果

4種方法重構(gòu)的GPS坐標(biāo)時(shí)間序列對(duì)比如圖1所示，可以看出，ML2比其他3種方法具有更好的重構(gòu)效果，且MOLS方法過(guò)擬合導(dǎo)致重構(gòu)GPS坐標(biāo)時(shí)間序列吸收了部分噪聲，而ML2方法則避免了這一點(diǎn)，因此效果更好。由表1、2可知，ML2方法擁有最小的Misfit，且噪聲譜指數(shù)和噪聲振幅與真實(shí)的最接近；SSA和MOLS方法的Misfit較大，且在噪聲較小時(shí)SSA和MOLS方法會(huì)吸收閃爍噪聲，導(dǎo)致噪聲譜指數(shù)的估計(jì)較差；WD方法會(huì)吸收大量的噪聲導(dǎo)致誤差較大，尤其是在噪聲較小時(shí)，WD方法與其他方法相比誤差更大，重構(gòu)效果最差；MOLS方法在重構(gòu)GPS坐標(biāo)時(shí)間序列時(shí)會(huì)吸收部分噪聲，從而導(dǎo)致過(guò)擬合，該現(xiàn)象在噪聲較小時(shí)尤為明顯，ML2方法則有效避免了這一缺點(diǎn)。4種方法的頻譜圖如圖2所示，period為真實(shí)周期，可以看出，這4種方法均能提取周期，MOLS、SSA和ML2方法均有較好表現(xiàn)，而WD方法則吸收了較多噪聲，效果較差，尤其在噪聲振幅較小時(shí)，WD方法失真。圖2也證實(shí)了ML2方法能減少過(guò)擬合，從而有更好的重構(gòu)效果。

圖1 4種方法的GPS坐標(biāo)時(shí)間序列重構(gòu)比較

圖2 4種方法的殘差功率譜比較

3 真實(shí)數(shù)據(jù)分析

本文選取的是北京房山站1999-2017年的數(shù)據(jù)，實(shí)際GPS坐標(biāo)時(shí)間序列缺失部分?jǐn)?shù)據(jù)，且含有野值（數(shù)據(jù)觀測(cè)錯(cuò)誤），因此需對(duì)數(shù)據(jù)進(jìn)行預(yù)處理。首先利用3σ準(zhǔn)則（拉依達(dá)準(zhǔn)則）去除野值，然后利用線性插值對(duì)缺失的部分進(jìn)行補(bǔ)全，最后進(jìn)行GPS坐標(biāo)時(shí)間序列的重構(gòu)。為了說(shuō)明ML2方法的有效性，本文通過(guò)計(jì)算相關(guān)系數(shù)來(lái)比較4種方法之間的相似性。

實(shí)驗(yàn)結(jié)果如表3所示，房山站U方向4種方法的相關(guān)性如表4所示，可以看出，ML2方法介于幾種方法之間，有著與其他方法相同的表現(xiàn)；ML2方法與MOLS方法和SSA方法的相關(guān)性較高與WD方法相關(guān)性較低；與MOLS方法相比，ML2方法減少了過(guò)擬合，從而具有更好的效果。4種方法U方向的功率譜如圖3所示，可以看出，ML2方法的重構(gòu)效果優(yōu)于MOLS方法，表明ML2方法能通過(guò)減少過(guò)擬合來(lái)提高GPS坐標(biāo)時(shí)間序列的重構(gòu)精度。

表3 房山站3個(gè)方向的實(shí)驗(yàn)結(jié)果

表4 房山站U方向4種方法的相關(guān)性

圖3 房山站U方向4種方法的功率譜

4 結(jié) 語(yǔ)

本文研究了坐標(biāo)時(shí)間序列的重構(gòu)，將時(shí)變的GPS坐標(biāo)時(shí)間序列模型轉(zhuǎn)化為一般的GPS坐標(biāo)時(shí)間序列模型。為了重構(gòu)GPS坐標(biāo)時(shí)間序列，本文采用ML2方法建立了優(yōu)化模型，以防止過(guò)擬合，在噪聲較大時(shí)能獲得更好的重構(gòu)效果；利用ADMM算法求解優(yōu)化模型，并采用滑動(dòng)的方法重構(gòu)GPS坐標(biāo)時(shí)間序列。通過(guò)模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)的實(shí)驗(yàn)得到以下結(jié)論：

1）在模擬數(shù)據(jù)實(shí)驗(yàn)中，ML2方法能重構(gòu)GPS坐標(biāo)時(shí)間序列，與MOLS、SSA和WD方法相比具有最小的Misfit。從模擬實(shí)驗(yàn)的結(jié)果來(lái)看，ML2方法比MOLS方法的效果更好，避免了過(guò)擬合問(wèn)題，使得重構(gòu)的精度更高；從真實(shí)數(shù)據(jù)實(shí)驗(yàn)的結(jié)果來(lái)看，ML2方法的效果介于幾種方法之間，證明了ML2的有效性，且與MOLS方法相比具有更好的重構(gòu)效果，這也表明了ML2方法能避免過(guò)擬合問(wèn)題，從而使重構(gòu)精度更高。2）SSA和WD方法在信號(hào)太弱的情形下難以準(zhǔn)確分離信號(hào)和噪聲，而ML2方法則在該情形下能夠得到較好的重構(gòu)效果，這也體現(xiàn)了ML2方法相較于其他重構(gòu)方法的優(yōu)越性。