奇異譜迭代區(qū)間四分法在GPS坐標(biāo)時(shí)間序列插補(bǔ)中的應(yīng)用

邱榮海 成英燕 王 虎 王曉明 曹炳強(qiáng)

1 山東科技大學(xué)測(cè)繪科學(xué)與工程學(xué)院，青島市前灣港路579號(hào)，266590

2 中國(guó)測(cè)繪科學(xué)研究院，北京市蓮花池西路28號(hào)，100830

3 皇家墨爾本理工大學(xué)數(shù)學(xué)與空間科學(xué)學(xué)院，澳大利亞墨爾本市

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

1 基于奇異譜分析的插補(bǔ)理論

奇異譜分析是一種從時(shí)間序列的動(dòng)力重構(gòu)出發(fā)并與經(jīng)驗(yàn)正交函數(shù)相聯(lián)系的統(tǒng)計(jì)技術(shù)。奇異譜分析是建立在相空間重構(gòu)基礎(chǔ)上的EOF分解，它首先計(jì)算了原始數(shù)據(jù)的滯后自協(xié)方差矩陣的特征向量空間，然后將原始數(shù)據(jù)滯后排列的矩陣向該正交空間進(jìn)行投影，最后對(duì)重構(gòu)相空間以主成分的方法完成數(shù)據(jù)的內(nèi)在結(jié)構(gòu)分析［6，8－9］。由于篇幅所限，對(duì)于奇異譜分析原理與方法在這里不作詳細(xì)介紹，具體可參考文獻(xiàn)［6，10］。

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

1）將時(shí)間序列分成3個(gè)部分：訓(xùn)練數(shù)據(jù)、交叉驗(yàn)證數(shù)據(jù)和待插補(bǔ)數(shù)據(jù)。

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

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

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

平均細(xì)分區(qū)間［1，M］為1、a、b、c、M（對(duì)應(yīng)的均方根誤差用f（1）、f（a）、f（b）、f（c）、f（M）表示），可根據(jù)以下5個(gè)條件判斷最小值屬于哪個(gè)區(qū)間或哪兩個(gè)相鄰區(qū)間：

①當(dāng)f（a）≤min（f（1），f（b））時(shí)，最小值位于［1，b］之間；

②當(dāng)f（b）≤min（f（a），f（c））時(shí)，最小值位于［a，c］之間；

③當(dāng)f（c）≤min（f（b），f（M））時(shí)，最小值位于［b，M］之間；

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

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

2 實(shí)例分析

選取ARTU 測(cè)站共1 461個(gè)單日解坐標(biāo)，圖1中原始時(shí)間序列分成3類不同的數(shù)據(jù)類型，訓(xùn)練數(shù)據(jù)長(zhǎng)度為1 201個(gè)，交叉驗(yàn)證數(shù)據(jù)共200個(gè)，待插補(bǔ)數(shù)據(jù)共60個(gè)，模擬的數(shù)據(jù)缺失率為17.8%。

圖1 ARTU 測(cè)站時(shí)間序列Fig.1 The time seires of ARTU station

2.1 插補(bǔ)效率

以N方向插補(bǔ)來(lái)具體說(shuō)明。取嵌入窗口長(zhǎng)度M為400，計(jì)算插值階數(shù)K為1、100、200、300、400，對(duì)應(yīng)的交叉驗(yàn)證均方根誤差分別為6.10、9.34、10.13、11.28、11.37mm，均方根誤差最小值落在K∈［1，100］之間；再細(xì)分區(qū)間［1，100］，計(jì)算K為25、50、75 的交叉驗(yàn)證均方根誤差分別為9.18、9.25、9.32mm，可見最小值位于K∈［1，25］之間；繼續(xù)計(jì)算K為6、12、18時(shí)對(duì)應(yīng)的均方根誤差分別為7.33、8.93、8.96mm，再次縮小最小值所在的范圍位于［1，6］之間；最后計(jì)算K為2、3、4、5的均方根誤差分別為4.49、3.35、2.15、2.80mm。由此可見，當(dāng)嵌入窗口為400時(shí)，插值階數(shù)4為最佳。

在計(jì)算量上，采用奇異譜迭代的區(qū)間四分法插值只需計(jì)算15次，均方根誤差就能夠找到最佳的K值，而用普通奇異譜迭代插值需要計(jì)算400次，前者在計(jì)算效率上得到很大的提高。對(duì)于插值結(jié)果的正確性，本文采用普通的奇異譜插值計(jì)算400次均方根誤差來(lái)尋找最佳的插值階數(shù)，結(jié)果如圖2所示。可以看出，當(dāng)嵌入窗口為400，均方根誤差最小值對(duì)應(yīng)的插值階數(shù)為4時(shí)，其結(jié)果與采用奇異譜迭代的區(qū)間四分法一致。

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

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

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

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

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

表2 IGS站最佳參數(shù)值對(duì)應(yīng)的交叉驗(yàn)證均方誤差Tab.2 The optimal parameter values corresponding to the mean square error of cross validation in IGS station

2.2 精度分析

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

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

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

圖3給出了兩種插值法的殘差圖。奇異譜插補(bǔ)的殘差值水平方向在4 mm 以內(nèi)，高程方向在10mm 以內(nèi)；拉格朗日插補(bǔ)殘差值水平方向上達(dá)10mm 左右，高程方向部分點(diǎn)達(dá)到20mm 以上，波動(dòng)較大，插值不穩(wěn)定。

圖4為最終插補(bǔ)的時(shí)間序列圖。由圖可見，插補(bǔ)效果較好。

其余4個(gè)IGS 站兩種插補(bǔ)法的精度比較結(jié)果如表4所示（單位：mm）。

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

圖4 奇異譜插補(bǔ)時(shí)序圖Fig.4 The time series diagram about singular spectrum interpolation

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

從表4可以看出，無(wú)論是國(guó)內(nèi)還是國(guó)際IGS站，奇異譜迭代的區(qū)間四分法插值精度在水平方向均優(yōu)于3mm，高程方向均優(yōu)于5mm；拉格朗日插值精度在10mm 左右，精度次于奇異譜迭代的區(qū)間四分法插值結(jié)果。

3 結(jié) 語(yǔ)

1）奇異譜迭代的區(qū)間四分法是一種新穎且有很大應(yīng)用前景的缺損數(shù)據(jù)插補(bǔ)方法，對(duì)于數(shù)據(jù)量大的時(shí)間序列插補(bǔ)效果較好。本文插補(bǔ)5個(gè)IGS測(cè)站，每個(gè)站的總數(shù)據(jù)量為1 461個(gè)，數(shù)據(jù)缺失率達(dá)到17.8%，插值精度水平方向優(yōu)于3mm，高程方向優(yōu)于5mm，插值精度高。

2）采用奇異譜迭代的區(qū)間四分法能夠有效地提高插值效率。以ARTU 測(cè)站N方向?yàn)槔?，采用一般奇異譜迭代插值需要計(jì)算400次均方根誤差，而采用奇異譜迭代的區(qū)間四分法只需計(jì)算15次，大幅度降低了計(jì)算量，并且能夠準(zhǔn)確地找出最佳的插值階數(shù)。

3）在與拉格朗日插值比較時(shí)發(fā)現(xiàn)，奇異譜迭代的區(qū)間四分法能夠從整體上完成缺失數(shù)據(jù)的插值，提高插值的真實(shí)可靠性。而拉格朗日插值因受階數(shù)的限制，只能從局部數(shù)據(jù)上完成插值，過(guò)度依賴局部數(shù)據(jù)，整體性不強(qiáng)。

［1］成英燕，丁繼新.我國(guó)高精度GPS陸海垂直運(yùn)動(dòng)監(jiān)測(cè)網(wǎng)的數(shù)據(jù)處理［J］，測(cè)繪科學(xué)，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］張鵬，蔣志浩，秘金鐘，等.我國(guó)GPS跟蹤站數(shù)據(jù)處理與時(shí)間序列特征分析［J］，武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，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高程時(shí)間序列粗差的抗差探測(cè)與插補(bǔ)研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，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］武艷強(qiáng)，黃立人.時(shí)間序列處理的新插值方法［J］.大地測(cè)量與地球動(dòng)力學(xué)，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］肖蒙，李軍.切比雪夫多項(xiàng)式及其插值法在檢測(cè)中的應(yīng)用研究［J］.自動(dòng)化與儀器儀表，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坐標(biāo)框架動(dòng)態(tài)特性及非線性建模方法研究［D］.北京：中國(guó)測(cè)繪科學(xué)研究院，2013（Wang Xiaoming.The Study on Movement Characteristics and Non－Linear Model of CGCS2000Framework［D］.Beijing：Chinese Academy of Surveying and Mapping，2013）

［7］王輝贊，張韌，劉巍，等.奇異譜迭代插補(bǔ)的改進(jìn)算法及其在缺失數(shù)據(jù)恢復(fù)中的應(yīng)用［J］.應(yīng)用數(shù)學(xué)和力學(xué)，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］吳洪寶，吳蕾.氣候變率診斷和預(yù)測(cè)方法［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］陳兆林，張書畢，佟瑞菊.用拉格朗日多項(xiàng)式內(nèi)插計(jì)算GPS衛(wèi)星位置［J］.全球定位系統(tǒng)，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）