基于滑動L2優化準則的GPS坐標時間序列重構

周 侗，邊家文*，丁開華，陳保周，伍浩琛

（1.中國地質大學 數學與物理學院，湖北 武漢 430074；2.中國地質大學 地理與信息工程學院，湖北 武漢 430074）

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

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

1 模型與重構算法

1.1 GPS坐標時間序列的模型

GPS坐標時間序列的一般模型為：

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

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

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

記：則式（2）可寫為：

式（2）可轉化為式（1），相應的矩陣表示形式為：

式中，A1X1為趨勢項；A2X2為周期項；ε為噪聲項。

1.2 基于L2優化準則的GPS坐標時間序列重構

1.2.1 優化模型

為了估計周期項和趨勢項，建立的優化模型為：

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

1.2.2 ADMM算法

ADMM算法既有乘子法的強收斂性質又有對偶上升法的分解性，因此該方法比較適合求解大規模的優化問題。目前，該方法已廣泛應用于壓縮感知、圖像處理、機器學習等諸多領域。ADMM算法的優點是僅有一個參數ρ，且在一定條件下，其對任意參數均收斂。對于優化問題式（8），可使用梯度下降法求解，但該方法對步長的選擇非常敏感，選擇不合適的步長甚至會導致算法不收斂；而ADMM算法對步長具有較強的魯棒性，在一定條件下，能保證對所有的步長都收斂[12]。因此，本文采用ADMM算法求解優化問題式（8）。

構造優化問題式（8）的增廣拉格朗日函數，則有：

式中，Z為對偶變量并具有適當的縮放尺寸；ρ為懲罰參數[13]。

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

對式（10）的優化目標函數進行簡化，可得：

對rk+1求導，令結果為0，則有：

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

對式（11）的優化目標函數進行簡化，則有：

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

ADMM算法的收斂性由其殘差來刻畫，其停止準則為：

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

1.2.3 算法評價指標和計算公式

本文采用4個評價指標全面衡量重構信號的效果，通過Misfit衡量重構信號和真實信號的誤差，噪聲譜指數和噪聲振幅衡量殘差和真實噪聲的接近程度，速度不確定度衡量殘差對速度估計的影響。

1）Misfit為真實信號和重構信號的標準差。2）噪聲譜指數（κ）。功率譜的計算公式為[14]：

兩邊取對數，則有：

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

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

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

ε的協方差矩陣C的計算公式為[16]：

式中，I為單位矩陣；Jκ為譜指數為κ的噪聲的單位協方差矩陣。

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

式中，ΔT為時間間隔；Γ為Gamma函數。

1.2.4 算法流程

通過上述分析，基于L2優化準則的GPS坐標時間序列重構方法的具體步驟為：①根據式（2）產生時變周期的GPS坐標時間序列；②根據式（3）建立矩陣A1和A2；③根據信號的長度N，選取合適的分段長度W，并逐年滑動得到M段信號（M=([N/365]-W)+2）；④針對每段信號建立式（8）的優化模型，再按照迭代步驟估計參數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；⑩根據式（22）和式（23）分別計算rpri和sdual；?若rpri和sdual滿足式（22）和式（23）則終止循環，否則繼續重復步驟⑥～⑨；?重復步驟④～?，直到得到每段時間序列的估計；?依據X1和X2逐年恢復信號（重復的年份取平均），并得到殘差；?根據式（25），利用最小二乘法得到殘差的譜指數估計；?根據式（26）～（30），利用最大似然法得到殘差的振幅估計；?計算重構信號的Misfit，并根據式（31）計算速度不確定度。

2 模擬實驗

本文采用兩組模擬數據來檢驗本文提出的算法性能，并與MOLS、SSA和WD方法進行對比。實驗中MOLS每段時間長度為1 100 d，滑動間隔為1 a，SSA的窗口大小為4 a（窗口大小為1 460 d），WD方法采用的是Meyer小波。本文利用Misfit、噪聲譜指數、噪聲振幅和速度不確定度來衡量重構信號的精度。

模擬數據的參數為：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。為了驗證算法在有色噪聲干擾情形下的有效性，本文模擬實驗的噪聲取為閃爍噪聲。閃爍噪聲較小時，閃爍噪聲振幅取為1 mm/a-1/4；閃爍噪聲較大時，閃爍噪聲振幅取為10 mm/a-1/4。實驗結果如表1、2所示。

表1 噪聲振幅為1 mm/a-1/4的實驗結果

表2 噪聲振幅為10 mm/a-1/4的實驗結果

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

圖1 4種方法的GPS坐標時間序列重構比較

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

3 真實數據分析

本文選取的是北京房山站1999-2017年的數據，實際GPS坐標時間序列缺失部分數據，且含有野值（數據觀測錯誤），因此需對數據進行預處理。首先利用3σ準則（拉依達準則）去除野值，然后利用線性插值對缺失的部分進行補全，最后進行GPS坐標時間序列的重構。為了說明ML2方法的有效性，本文通過計算相關系數來比較4種方法之間的相似性。

實驗結果如表3所示，房山站U方向4種方法的相關性如表4所示，可以看出，ML2方法介于幾種方法之間，有著與其他方法相同的表現；ML2方法與MOLS方法和SSA方法的相關性較高與WD方法相關性較低；與MOLS方法相比，ML2方法減少了過擬合，從而具有更好的效果。4種方法U方向的功率譜如圖3所示，可以看出，ML2方法的重構效果優于MOLS方法，表明ML2方法能通過減少過擬合來提高GPS坐標時間序列的重構精度。

表3 房山站3個方向的實驗結果

表4 房山站U方向4種方法的相關性

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

4 結 語

本文研究了坐標時間序列的重構，將時變的GPS坐標時間序列模型轉化為一般的GPS坐標時間序列模型。為了重構GPS坐標時間序列，本文采用ML2方法建立了優化模型，以防止過擬合，在噪聲較大時能獲得更好的重構效果；利用ADMM算法求解優化模型，并采用滑動的方法重構GPS坐標時間序列。通過模擬數據和真實數據的實驗得到以下結論：

1）在模擬數據實驗中，ML2方法能重構GPS坐標時間序列，與MOLS、SSA和WD方法相比具有最小的Misfit。從模擬實驗的結果來看，ML2方法比MOLS方法的效果更好，避免了過擬合問題，使得重構的精度更高；從真實數據實驗的結果來看，ML2方法的效果介于幾種方法之間，證明了ML2的有效性，且與MOLS方法相比具有更好的重構效果，這也表明了ML2方法能避免過擬合問題，從而使重構精度更高。2）SSA和WD方法在信號太弱的情形下難以準確分離信號和噪聲，而ML2方法則在該情形下能夠得到較好的重構效果，這也體現了ML2方法相較于其他重構方法的優越性。