奇異譜迭代區間四分法在GPS坐標時間序列插補中的應用

邱榮海 成英燕 王 虎 王曉明 曹炳強

1 山東科技大學測繪科學與工程學院，青島市前灣港路579號，266590

2 中國測繪科學研究院，北京市蓮花池西路28號，100830

3 皇家墨爾本理工大學數學與空間科學學院，澳大利亞墨爾本市

GPS 坐標時間序列分析在變形監測領域發揮著重要作用［1－2］。GPS時間序列分析之前，需對所研究的時間序列進行預處理，數據處理中難免會遇到數據缺失和數據質量不合格的情況。因此如何對缺損數據進行插補，提取有效信息，是挖掘和拓展數據信息資源的重要途徑。很多學者用不同方法進行插值，如三次樣條插值法、切比雪夫多項式插值法等。但是三次樣條插值無法應用于缺失數據較多的時間序列，在缺失數據嚴重的坐標時間序列插值中應用受限［3－4］。切比雪夫插值受插值階數的限制，階數高時會出現病態矩陣［5］。本文基于奇異譜迭代的思想，對GPS坐標時間序列進行插值。雖然有些學者利用該思想完成了一些數據實驗，也得到了較好的結果［6］，但不同時間序列內部數據結構不同，插值效果可能存在差異。同時，奇異譜迭代插值效率低也會影響到奇異譜插值法的推廣。本文采用奇異譜迭代的區間四分法［7］完成5個IGS站數據的插值，并與拉格朗日插值法進行比較。

1 基于奇異譜分析的插補理論

奇異譜分析是一種從時間序列的動力重構出發并與經驗正交函數相聯系的統計技術。奇異譜分析是建立在相空間重構基礎上的EOF分解，它首先計算了原始數據的滯后自協方差矩陣的特征向量空間，然后將原始數據滯后排列的矩陣向該正交空間進行投影，最后對重構相空間以主成分的方法完成數據的內在結構分析［6，8－9］。由于篇幅所限，對于奇異譜分析原理與方法在這里不作詳細介紹，具體可參考文獻［6，10］。

奇異譜迭代的區間四分法具體流程如下［7］：

1）將時間序列分成3個部分：訓練數據、交叉驗證數據和待插補數據。

2）將訓練數據去中心化并記錄其平均值，交叉驗證數據和待插補數據用0填補，得到新的時間序列x（n）。

3）對時間序列x（n）作奇異譜分解，根據具體情況選擇不同的嵌入窗口M，將區間［1，M］平均分成4個小區間，取K為1、M／4、M／2、3M／4、M，分別取前K個主要成分得到重構時間序列x1（n），將x（n）序列中缺失的值用x1（n）對應值替代，得到新的時間序列x2（n）。

4）如果max｜x2（n）－x1（n）｜≤0.000 1，則退出循環，計算不同嵌入窗口M和插值階數K對應的均方根誤差，判斷最小均方誤差屬于哪個區間或者哪兩個相鄰的區間，把縮小范圍后最小值所在的區間再細分為4個小區間，重復第3）步，再判斷最小值屬于哪個更小的區間。依次下去，最后確定最小的均方根誤差所對應的插值階數K值，M和K的插值結果作為最終插值。

平均細分區間［1，M］為1、a、b、c、M（對應的均方根誤差用f（1）、f（a）、f（b）、f（c）、f（M）表示），可根據以下5個條件判斷最小值屬于哪個區間或哪兩個相鄰區間：

①當f（a）≤min（f（1），f（b））時，最小值位于［1，b］之間；

②當f（b）≤min（f（a），f（c））時，最小值位于［a，c］之間；

③當f（c）≤min（f（b），f（M））時，最小值位于［b，M］之間；

④若前3條都不滿足，當f（a）＜f（b）時，最小值位于［1，a］之間；

⑤若前4條都不滿足，當f（b）＞f（c）時，最小值位于［c，M］之間。

2 實例分析

選取ARTU 測站共1 461個單日解坐標，圖1中原始時間序列分成3類不同的數據類型，訓練數據長度為1 201個，交叉驗證數據共200個，待插補數據共60個，模擬的數據缺失率為17.8%。

圖1 ARTU 測站時間序列Fig.1 The time seires of ARTU station

2.1 插補效率

以N方向插補來具體說明。取嵌入窗口長度M為400，計算插值階數K為1、100、200、300、400，對應的交叉驗證均方根誤差分別為6.10、9.34、10.13、11.28、11.37mm，均方根誤差最小值落在K∈［1，100］之間；再細分區間［1，100］，計算K為25、50、75 的交叉驗證均方根誤差分別為9.18、9.25、9.32mm，可見最小值位于K∈［1，25］之間；繼續計算K為6、12、18時對應的均方根誤差分別為7.33、8.93、8.96mm，再次縮小最小值所在的范圍位于［1，6］之間；最后計算K為2、3、4、5的均方根誤差分別為4.49、3.35、2.15、2.80mm。由此可見，當嵌入窗口為400時，插值階數4為最佳。

在計算量上，采用奇異譜迭代的區間四分法插值只需計算15次，均方根誤差就能夠找到最佳的K值，而用普通奇異譜迭代插值需要計算400次，前者在計算效率上得到很大的提高。對于插值結果的正確性，本文采用普通的奇異譜插值計算400次均方根誤差來尋找最佳的插值階數，結果如圖2所示?？梢钥闯?，當嵌入窗口為400，均方根誤差最小值對應的插值階數為4時，其結果與采用奇異譜迭代的區間四分法一致。

圖2 ARTU 站N 方向不同插值階數對應的均方根誤差Fig.2 The different interpolation order corresponding to the root mean square error in N direction of ARTU station

吳洪寶等［11］建議，當窗口長度為M時，能夠較好地提取周期為M／5～M的震蕩。又因為坐標序列中周期成分主要以年周期項與半年周期項為主，在插值中為了能夠有效地插入信息的主要成分，M值選擇上不宜過小。但是窗口長度M又要考慮客觀條件，為了防止最大滯后M的自協方差估計的統計誤差不超過估計量本身，M不宜超過總長度的1／3。綜合上述理由，結合插值方法，窗口長度M和所對應的最佳插值階數K取值如表1所示，同時給出奇異譜迭代的區間四分法與普通奇異譜迭代法計算次數的比較。表1中改進前表示采用普通的奇異譜迭代插值法的計算次數，改進后表示采用奇異譜迭代的區間四分法計算次數。

由表1，奇異譜迭代的區間四分法在計算次數上要比普通奇異譜迭代插值法少很多。普通奇異譜迭代插值需要幾百次的計算，而奇異譜迭代的區間四分法只需十幾次運算就能準確找到最佳參數值，計算效率得到很大提高。同時，N方向上的嵌入窗口M為400，插值階數為4，E方向的嵌入窗口為450，插值階數為5，U方向的嵌入窗口為450，插值階數是7，此為最佳參數值。

采用該方法尋找其余4 個IGS 站的最佳參數值，見表2。

表1 ARTU 站不同參數值的計算次數Tab.1 The number of calculation about different parameter values in ARTU station

表2 IGS站最佳參數值對應的交叉驗證均方誤差Tab.2 The optimal parameter values corresponding to the mean square error of cross validation in IGS station

2.2 精度分析

插值精度反映了插值效果的優越性。本文采用拉格朗日插值法［12］與奇異譜迭代的區間四分法作比較。拉格朗日插值數據段為交叉驗證數據和待插補數據兩部分，選擇的最佳插值階數為10階，插值點位于插值數據的中間。兩種方法的均方誤差和平均誤差如表3所示（單位：mm）。

表3 ARTU 站兩種插補方法精度信息Tab.3 The precision of two interpolation methods in ARTU station

從表3可以看出，奇異譜迭代的區間四分法水平方向上均方誤差小于2mm，高程方向不到5 mm，平均誤差水平方向在2mm 以下，高程方向在3mm 左右。而拉格朗日插值法插值誤差結果均比奇異譜插值大。

圖3給出了兩種插值法的殘差圖。奇異譜插補的殘差值水平方向在4 mm 以內，高程方向在10mm 以內；拉格朗日插補殘差值水平方向上達10mm 左右，高程方向部分點達到20mm 以上，波動較大，插值不穩定。

圖4為最終插補的時間序列圖。由圖可見，插補效果較好。

其余4個IGS 站兩種插補法的精度比較結果如表4所示（單位：mm）。

圖3 ARTU 站兩種插值法殘差對比圖Fig.3 The two interpolation error diagram about ARTU station

圖4 奇異譜插補時序圖Fig.4 The time series diagram about singular spectrum interpolation

表4 IGS站兩種插補方法插值精度Tab.4 The precision of two interpolation methods in IGS stations

從表4可以看出，無論是國內還是國際IGS站，奇異譜迭代的區間四分法插值精度在水平方向均優于3mm，高程方向均優于5mm；拉格朗日插值精度在10mm 左右，精度次于奇異譜迭代的區間四分法插值結果。

3 結 語

1）奇異譜迭代的區間四分法是一種新穎且有很大應用前景的缺損數據插補方法，對于數據量大的時間序列插補效果較好。本文插補5個IGS測站，每個站的總數據量為1 461個，數據缺失率達到17.8%，插值精度水平方向優于3mm，高程方向優于5mm，插值精度高。

2）采用奇異譜迭代的區間四分法能夠有效地提高插值效率。以ARTU 測站N方向為例，采用一般奇異譜迭代插值需要計算400次均方根誤差，而采用奇異譜迭代的區間四分法只需計算15次，大幅度降低了計算量，并且能夠準確地找出最佳的插值階數。

3）在與拉格朗日插值比較時發現，奇異譜迭代的區間四分法能夠從整體上完成缺失數據的插值，提高插值的真實可靠性。而拉格朗日插值因受階數的限制，只能從局部數據上完成插值，過度依賴局部數據，整體性不強。

［1］成英燕，丁繼新.我國高精度GPS陸海垂直運動監測網的數據處理［J］，測繪科學，2000，25（3）：40－41（Cheng Yingyan，Ding Jixin.Data Processing of High Precision GPS Vertical Movement Monitoring Network in Land and Sea［J］.Science of Surveying and Mapping，2000，25（3）：40－41）

［2］張鵬，蔣志浩，秘金鐘，等.我國GPS跟蹤站數據處理與時間序列特征分析［J］，武漢大學學報：信息科學版，2007，32（3）：251－254（Zhang Peng，Jiang Zhihao，Bei Jinzhong，et al.Data Processing and Time Series Analysis for GPS Fiducial Station in China［J］.Geomatics and Information Science of Wuhan University，2007，32（3）：251－254）

［3］張恒璟，程鵬飛.基于GPS高程時間序列粗差的抗差探測與插補研究［J］.大地測量與地球動力學，2011，31（4）：71－75（Zhang Hengjing，Cheng Pengfei.Study on Robust Detection and Interpolation from Gross Errors of GPS Height Series［J］.Journal of Geodesy and Geodynamics，2011，31（4）：71－75）

［4］武艷強，黃立人.時間序列處理的新插值方法［J］.大地測量與地球動力學，2004，24（4）：43－47（Wu Yanqiang，Huang Liren.A New Interpolation Method in Time Series Analyzing［J］.Journal of Geodesy and Geodynamics，2004，24（4）：43－47）

［5］肖蒙，李軍.切比雪夫多項式及其插值法在檢測中的應用研究［J］.自動化與儀器儀表，2006（3）：13－16（Xiao Meng，Li Jun.Study on Measurement Using Tchebyshev Polynomial and Its Interpolating Algorithm［J］.Journal of Automation and Instrumentaton，2006（3）：13－16）

［6］王曉明.CGCS2000坐標框架動態特性及非線性建模方法研究［D］.北京：中國測繪科學研究院，2013（Wang Xiaoming.The Study on Movement Characteristics and Non－Linear Model of CGCS2000Framework［D］.Beijing：Chinese Academy of Surveying and Mapping，2013）

［7］王輝贊，張韌，劉巍，等.奇異譜迭代插補的改進算法及其在缺失數據恢復中的應用［J］.應用數學和力學，2008，29（10）：1 227－1 236（Wang Huizan，Zhang Ren，Liu Wei，et al.Improved Interpolation Method Based on Singular Spectrum Analysis Iteration and Its Application in Missing Data Recovery［J］.Applied Mathematics and Mechanics，2008，29（10）：1 227－1 236）

［8］Vautard R，Yiou P，Ghil M.Singular－Spectrum Analysis：A Toolkit for Short，Noisy Chaotic Signals［J］.Physica D Nonlinear Phenomena，1992，58（1－4）：95－126

［9］Schoellhamer，David H.Singular Spectrum Analysis for Time Series with Missing Data［J］.Geophysical Research Letters，2001，28（16）：3 187－3 190

［10］吳洪寶，吳蕾.氣候變率診斷和預測方法［M］.北京：氣象出版社（Wu Hongbao，Wu Lei.Methods for Diagnosing and Forecasting Climate Variability［M］.China Meteorological Press，2005）

［11］吳洪寶.奇異譜和多通道奇異譜分析［J］.氣象教育與科技，1997（4）：1－9（Wu Hongbao.Singular and Multichannel Spectrum Analysis［J］.Meteorological Education and Science and Technology，1997（4）：1－9）

［12］陳兆林，張書畢，佟瑞菊.用拉格朗日多項式內插計算GPS衛星位置［J］.全球定位系統，2007，32（2）：33－35（Chen Zhaolin，Zhang Shubi，Tong Ruiju.Computing the Location of GPS Satellites by Lagrange Polynomial［J］.GNSS World of China，2007，32（2）：33－35）