基于LS+AR模型的三種方式預(yù)報地球自轉(zhuǎn)參數(shù)的比較

雷雨，高玉平，蔡宏兵

（1.中國科學(xué)院國家授時中心，西安 710600；2.中國科學(xué)院時間頻率基準(zhǔn)重點實驗室，西安 710600；3.中國科學(xué)院大學(xué)，北京100049）

基于LS+AR模型的三種方式預(yù)報地球自轉(zhuǎn)參數(shù)的比較

雷雨1，2，3，高玉平1，2，蔡宏兵1，2

（1.中國科學(xué)院國家授時中心，西安 710600；2.中國科學(xué)院時間頻率基準(zhǔn)重點實驗室，西安 710600；3.中國科學(xué)院大學(xué)，北京100049）

為研究最小二乘（least squares，LS）模型和自回歸（autoregression，AR）模型的組合（LS+AR）方法用于地球自轉(zhuǎn)參數(shù)（Earth rotation parameters，ERP）的預(yù)報時，不同的預(yù)報方式對預(yù)報結(jié)果的影響，我們采用遞推、迭代和間隔這3種預(yù)報方式對ERP進行預(yù)報。結(jié)果表明，這3種方式對日長變化（length of day，LOD）所有跨度預(yù)報的精度相當(dāng)，而遞推方式在極移所有跨度的預(yù)報上表現(xiàn)出精度優(yōu)勢，間隔方式次之，迭代方式最差。在數(shù)據(jù)利用率和計算速度方面，遞推和迭代方式的數(shù)據(jù)利用率高，但前者的計算量明顯小于后者，而間隔方式的數(shù)據(jù)利用率低，但計算速度最快。

地球自轉(zhuǎn)參數(shù)；預(yù)報；LS+AR模型；遞推方式；迭代方式；間隔方式

0 引言

地球自轉(zhuǎn)參數(shù)（Earth rotation parameters，ERP）表征了地球自轉(zhuǎn)的實時轉(zhuǎn)動運動狀態(tài)[1]，包括日長變化（length of day，LOD）和極移（polar motion，PM）X，Y分量PMX，PMY。ERP是實現(xiàn)地球參考坐標(biāo)系與天球參考坐標(biāo)系相互轉(zhuǎn)換的必需參數(shù)，在深空探測、衛(wèi)星精密定軌和天文地球動力學(xué)研究等領(lǐng)域都有重要應(yīng)用[2-3]。衛(wèi)星激光測距（satellite laser ranging，SLR）、甚長基線干涉（very long baseline interferometry，VLBI）和全球定位系統(tǒng)（global positioning system，GPS）等現(xiàn)代測地技術(shù)被廣泛用于地球自轉(zhuǎn)變化的常規(guī)測量中，提供了高時空分辨率和高精度的ERP觀測資料[3-4]。但是由于復(fù)雜的數(shù)據(jù)處理過程，由現(xiàn)代測地技術(shù)獲取的ERP存在幾個小時至幾天的延遲，而衛(wèi)星精密定軌等領(lǐng)域?qū)崟rERP有重要需求[5]，因此，研究高精度的ERP預(yù)報方法具有十分重要的現(xiàn)實意義。

目前預(yù)報ERP的主要方法有：1）小波分解（wavelet decomposition，WD）與自協(xié)方差（auto-covariance，AC）模型的組合（WD+AC）[6]；2）大氣角動量（atmospheric angular momentum，AAM）預(yù)報方法[7]；3）卡爾曼濾波（Kalman filter）方法[8]；4）最小二乘（least squares，LS）模型與自回歸（autoregressive，AR）模型的組合（LS+AR）[9-10]；5）LS與人工神經(jīng)網(wǎng)絡(luò)（artificial neural networks，ANN）模型的組合（LS+ANN）[10-12]；6）頻譜分析與LS模型的組合[13]。

為了檢驗各種方法的預(yù)報效果，維也納大學(xué)大地測量與地球物理研究所從2005年10月1日至2008年2月28日組織了全球性的地球定向參數(shù)預(yù)報比較研究項目（Earth orientation parameters prediction comparison campaign，EOP PCC）[14]，通過2年多的預(yù)報比較，EOP PCC獲得一個重要結(jié)論：沒有一種預(yù)測方法既適合于ERP所有分量的預(yù)報，又適合于所有跨度的預(yù)報，上述的前3種方法更適合于LOD預(yù)報，而后3種方法則對PMX，PMY預(yù)報稍有優(yōu)勢。

EOP PCC比較了不同預(yù)報方法的預(yù)報精度情況，但沒有分析同一預(yù)報方法下不同預(yù)報方式的適用特點。預(yù)報方式一般分為迭代預(yù)報和間隔預(yù)報兩種，間隔預(yù)報方式實現(xiàn)簡單，計算速度快，但該方式要求較長的基礎(chǔ)序列，如果序列較短，隨著預(yù)報跨度的增大，數(shù)據(jù)利用率顯著降低，對ERP中長期預(yù)報精度不理想，且對規(guī)律項不明顯的序列預(yù)報效果不好；而迭代預(yù)報方式中，數(shù)據(jù)利用不受預(yù)報跨度的影響，數(shù)據(jù)利用率高，且對規(guī)律項不明顯的序列預(yù)報效果好。這兩種方式被廣泛用于ERP預(yù)報中，如王琪潔等采用LS+ANN模型間隔方式預(yù)報LOD[11]；嚴(yán)鳳等采用LS+AR模型迭代方式預(yù)報極移和LOD[15]；許雪晴等對比了AR模型迭代方式和間隔方式的ERP預(yù)報精度，該研究表明，迭代方式適合LOD所有跨度預(yù)報和極移中長期預(yù)報，而間隔方式則較適合極移短期預(yù)報[16]。

本文基于經(jīng)典的LS+AR預(yù)報模型，對實測ERP數(shù)據(jù)分別進行遞推方式、迭代方式和間隔方式的預(yù)報，并與實測值進行比較，對比這3種方式的預(yù)報精度，分析各自的適用性。

1 LS+AR模型及其算法流程

1.1 LS模型

地球自轉(zhuǎn)變化包含趨勢項和周期項等復(fù)雜過程，其中趨勢項和周期項部分采用LS模型進行擬合，然后外推預(yù)報[9]。LS模型的基本形式如下：

式（1）中，xt（t=1，2，…，N）為觀測序列，a0和a1為趨勢項擬合系數(shù)，m為周期項個數(shù)，Ai，Pi和φi分別為周期項振幅、周期和相位，t為UTC時。對于PMX，PMY的擬合，該模型包含3個周期項，周期分別為1.183，1和0.5周年，而關(guān)于LOD的擬合，該模型包含4個周期項，周期分別為18.6，9.3，1和0.5周年。式（1）中的模型參數(shù)可根據(jù)最小二乘原理求得。

1.2 AR模型

AR模型是指根據(jù)時間序列變量自身過去的規(guī)律來建立預(yù)報模型[17]，其數(shù)學(xué)模型為

式（2）中，zt（t=1，2，…，N）為平穩(wěn)時間序列，p為模型階次，φi（i=1，2，…，p）為模型系數(shù)，εt為白噪聲。

AR模型使用的關(guān)鍵是模型階次p的確定，常用的定階準(zhǔn)則有最終預(yù)報誤差準(zhǔn)則（final prediction error criterion，F(xiàn)PE）、傳遞函數(shù)準(zhǔn)則和信息量準(zhǔn)則等3種[18]，這3種準(zhǔn)則在實際應(yīng)用中是等效的。本文采用應(yīng)用相對廣泛的赤池信息量準(zhǔn)則（Akaike information criterion，AIC）來確定AR模型的階次p[19]，其準(zhǔn)則函數(shù)為

1.3 算法流程

LS+AR模型建模預(yù)報的算法流程如圖1所示。其預(yù)報過程為：首先利用LS模型擬合趨勢項和周期項，然后采用AR模型對LS擬合殘差序列進行建模預(yù)報，最后對LS模型的外推值和AR模型的預(yù)報值疊加獲得最終的ERP預(yù)報值。

圖1　LS+AR模型算法流程圖

2 預(yù)報方式

2.1 迭代預(yù)報方式

然后把zN+1加入原有序列zt（t=1，2，…，N）中組成新序列zt（t=1，2，…，N +1），再重新建立AR模型，進而外推得到第2天的預(yù)報值zN+2，依次類推，即可獲得第M天的預(yù)報值zN+M：

迭代預(yù)報方式充分利用了基礎(chǔ)序列信息，并且在預(yù)報過程中反復(fù)建立AR模型，因此最終確定的AR模型會受到預(yù)報誤差的影響。對于跨度為M天的預(yù)報，這種預(yù)報方式需要對基礎(chǔ)序列作M次AR模型運算，所以計算量比較大，耗費的時間會長一些。迭代預(yù)報方式的流程圖可用圖2表示。對于M天的預(yù)報，迭代預(yù)報方式的計算復(fù)雜度為O（ N）+O（ N+1）+…+O（ N+M -1），其中O（ N）表示序列長度為N時AR模型的計算復(fù)雜度。

圖2　AR模型迭代預(yù)報方式流程圖

2.2 遞推預(yù)報方式

遞推預(yù)報方式的流程如圖3所示。迭代預(yù)報方式與遞推預(yù)報方式的區(qū)別在于，對于M天的預(yù)報，前者需要建立M個AR模型，故需要作M次AR模型計算，而后者只建立1個AR模型，故只需作1次AR模型計算，計算量大為減少，其計算復(fù)雜度為（O）N。

圖3　AR模型遞推預(yù)報方式流程圖

2.3間隔預(yù)報方式

式（7）中，zN+M為第M天的預(yù)報值。間隔預(yù)報方式實現(xiàn)起來簡單，對所有跨度預(yù)報均只需作1次AR模型運算，其計算復(fù)雜度為O（int（（N-1）/M）），計算量比迭代方式和遞推方式要少，但隨著預(yù)報時長的增大，基礎(chǔ)預(yù)報序列長度迅速縮短，AR模型建模誤差大。間隔預(yù)報方式的流程如圖4所示。

圖4　AR模型間隔預(yù)報方式流程圖

3 實驗結(jié)果與分析

根據(jù)上述預(yù)報方式，利用LS+AR模型分別對PMX，PMY和LOD作1～90 d不同跨度的預(yù)報，對比分析3種預(yù)報方式的預(yù)報效果以及適用情況。本文所用ERP數(shù)據(jù)來自國際地球自轉(zhuǎn)與參考系服務(wù)（International Earth Rotation and Reference Service，IERS）的EOP 05C04序列，數(shù)據(jù)采樣間隔為1 d，選用1990年1月1日0時至1999年12月31日0時（UTC時）共10 a的ERP序列作為初始基礎(chǔ)序列，預(yù)報從2000年1月1日0時到2005年2月6日0時為止的ERP值，每隔5 d預(yù)報1次，預(yù)報起始時刻每5 d更新1次，基礎(chǔ)序列長度始終保持不變，這樣共產(chǎn)生356組預(yù)報序列，每組預(yù)報序列包含1～90 d跨度的預(yù)報值。

需要說明的是，與極移預(yù)報過程不同，對于LOD預(yù)報，首先需要對其進行數(shù)據(jù)預(yù)處理，即首先根據(jù)IERS規(guī)范扣除LOD序列中周期為5 d～18.6 a的62個固體地球帶諧潮汐項[20]，再對改正潮汐項的LOD（稱為LODR）序列進行建模預(yù)報，最后將固體地球帶諧潮汐項加入LODR的預(yù)報序列，得到LOD預(yù)報值。

為了評估預(yù)報效果，選用平均絕對誤差（mean absolute error，MAE）作為預(yù)報結(jié)果的精度評定指標(biāo)，預(yù)報跨度為M的預(yù)報結(jié)果的平均絕對誤差MAEM的計算公式為

表1分別給出了3種預(yù)報方式的ERP預(yù)報精度。為了更加直觀地對比3種方式的預(yù)報結(jié)果，圖5～ 7給出了3種方式在不同跨度情況下的預(yù)報精度比較圖。

結(jié)合表1和圖5～7可以看出：

1）對于LOD預(yù)報而言，3種方式的預(yù)報精度大體相當(dāng)，其中遞推方式較另外2種方式稍有精度優(yōu)勢，而間隔方式的預(yù)報精度則隨著預(yù)報時長的增大表現(xiàn)出不穩(wěn)定的趨勢；

2）遞推方式對PMX和PMY所有跨度預(yù)報的精度明顯高于另外2種方式，間隔方式的預(yù)報精度居中，在超短期（1～10 d）的預(yù)報中精度略遜于遞推方式，而迭代方式的預(yù)報誤差則隨著預(yù)報時長的增大表現(xiàn)出明顯的增大趨勢，這種現(xiàn)象在LOD預(yù)報中也有所體現(xiàn)，只是不如極移預(yù)報那么明顯，這是由于用預(yù)報序列反復(fù)建立AR模型而導(dǎo)致的誤差累積效應(yīng)。

以下根據(jù)第1節(jié)介紹的預(yù)報方法和第2節(jié)闡述的預(yù)報方式算法流程，分析3種預(yù)報方式的數(shù)據(jù)利用率和計算耗時。對長度為N的基礎(chǔ)序列作跨度為M天的預(yù)報，遞推預(yù)報方式利用的建模序列長度始終為N，所以其數(shù)據(jù)利用率始終保持不變，迭代預(yù)報方式可利用的建模序列長度隨著預(yù)報跨度的增大而增加，故數(shù)據(jù)利用率非常高，而間隔預(yù)報方式可利用的建模序列長度為int（（N-1）/M），int（）代表取整，因此其數(shù)據(jù)利用率隨著預(yù)報跨度的增加而迅速降低，故當(dāng)預(yù)報跨度較大時，對基礎(chǔ)序列長度有一定要求。從第2節(jié)所述的計算復(fù)雜度不難得出，間隔方式的計算速度最快，最省時，遞推方式的計算速度相對于間隔方式而言稍慢些，但由于該方式只需要進行1次AR模型的運算，故仍然有較高的計算效率，而迭代模型由于需要反復(fù)建立AR模型，所以計算會相當(dāng)耗時。

表1　基于LS+AR模型用遞推、迭代和間隔預(yù)報方式預(yù)報ERP的精度（MAEM）比較

圖5　基于LS+AR模型用3種方式預(yù)報PMX的精度比較

圖6　基于LS+AR模型用3種方式預(yù)報PMY的精度比較

圖7　基于LS+AR模型用3種方式預(yù)報LOD的精度比較

4 結(jié)語

結(jié)合經(jīng)典的LS+AR預(yù)報模型，闡述了遞推預(yù)報方式、迭代預(yù)報方式和間隔預(yù)報方式的算法流程。選取IERS 05C04 ERP數(shù)據(jù)作為基礎(chǔ)預(yù)報序列，分別利用上述3種預(yù)報方式對ERP進行預(yù)報，并與實際觀測值進行比較，統(tǒng)計這3種方式的預(yù)報精度。預(yù)報結(jié)果表明，這3種方式對LOD所有跨度預(yù)報的精度基本相當(dāng)，但遞推方式在PMX和PMY所有跨度的預(yù)報上較另外2種方式有明顯的精度優(yōu)勢，迭代方式的預(yù)報誤差隨著預(yù)報時長的增大表現(xiàn)出明顯的增大趨勢，精度最差。在數(shù)據(jù)利用率和計算量方面，迭代預(yù)報方式的數(shù)據(jù)利用率最高，但計算量相對較大，而間隔預(yù)報方式的數(shù)據(jù)利用率低，對基礎(chǔ)序列長度有一定要求，計算速度最快，遞推預(yù)報方式在保證數(shù)據(jù)利用率的同時又具有較快的計算速度。因此，間隔預(yù)報方式適合于對實時性要求較高的場合，例如ERP的實時快速預(yù)報。

實際應(yīng)用中可以將間隔方式和遞推方式兩者結(jié)合起來用于ERP預(yù)報，即在超短期的實時快速預(yù)報中采用間隔預(yù)報方式，而在短期和中長期預(yù)報中采用遞推預(yù)報方式。

[1] EUBANKS T M.Variations in the Orientation of the Earth[M].Washington DC: AGU Press, 1993: 1-54.

[2] 魏二虎, 史青, 嚴(yán)韋, 等.ERP精度對“嫦娥一號”差分VLBI定位精度的影響[J].武漢大學(xué)學(xué)報(信息科學(xué)版), 2011, 36(11): 1324-1327.

[3] GAMBIS D, LUZUM B, LEWANDOWSKI W, et al.Earth rotation monitoring, UT1 determination and prediction[J].Metrologia, 2011, 48(4): 165-170.

[4] 何戰(zhàn)科, 楊旭海, 李志剛, 等.利用GPS觀測資料解算地球自轉(zhuǎn)參數(shù)[J].時間頻率學(xué)報, 2010, 33(1): 69-76.

[5] 張衛(wèi)星, 劉萬科, 龔曉穎.EOP預(yù)報誤差對自主定軌結(jié)果影響分析[J].大地測量與地球動力學(xué), 2011, 31(5): 106-110.

[6] KOSEK W, KALARUS M, JOHNSON T J, et al.A comparison of LOD and UT1-UTC forecasts by different combined prediction techniques[J].Artificial Satellites, 2005, 40(2): 119-125.

[7] JOHNSON T J, LUZUM B J, RAY J R.Improved near-term UT1R predictions using forecasts of atmospheric angular momentum[J].Journal of Geodynamics, 2005, 39(3): 209-221.

[8] GROSS R S, EUBANKS T M, STEPPE J A, et al.A Kalman filter-based approach to combining independent Earth-orientation series[J].Journal of Geodesy, 1998, 72(4): 215-235.

[9] TOMASZ N, KOSEK W.Prediction of UT1-UTC, LOD and AAM3χ by combination of least-squares and multivariate stochastic methods[J].Journal of Geodesy, 2008, 82(2): 83-92.

[10] 許雪晴, 周永宏.地球定向參數(shù)高精度預(yù)報方法研究[J].飛行器測控學(xué)報, 2010, 29(2): 70-76.

[11] 王琪潔, 廖得春, 周永宏.地球自轉(zhuǎn)速率變化的實時快速預(yù)報[J].科學(xué)通報, 2007, 52(15): 1728-1731.

[12] SCHUH H, ULRICH M, EGGER D, et al.Prediction of Earth orientation parameters by artificial neural networks[J].Journal of Geodesy, 2002, 76(5): 247-258.

[13] AKULENKO L D, KUMAKSHEV S A, MARKOV Y G, et al.Forecasting the polar motions of the deformable Earth[J].Astronomy Reports, 2002, 46(10): 858-865.

[14] KALARUS M, SCHUH H, KOSEK W, et al.Achievements of the Earth orientation parameters prediction comparison campaign[J].Journal of Geodesy, 2010, 84(10): 587-596.

[15] 嚴(yán)鳳, 姚宜斌.地球自轉(zhuǎn)參數(shù)短期預(yù)報方法及其實現(xiàn)[J].大地測量與地球動力學(xué), 2012, 32(4): 71-75.

[16] 許雪晴, 周永宏.AR模型間隔方式和迭代方式預(yù)報地球定向參數(shù)對比[J].中國科學(xué)院上海天文臺年刊, 2012, 33: 20-27.

[17] 劉大杰, 陶本藻, 朱建軍, 等.實用測量數(shù)據(jù)處理方法[M].北京: 測繪出版社, 2000: 83-89.

[18] 張善文, 雷英杰, 馮有前, 等.MATLAB在時間序列分析中的應(yīng)用[M].西安: 西安電子科技大學(xué)出版社, 2007: 130-133.

[19] AKAIKE H, MINATO KU.A new look at the statistical model identification[J].IEEE Transactions on Automatic Control, 1974, 19(6): 716-723.

[20] GERARD P, BRAIN L.IERS Conventions (2010)[S].Germany: Verlag des Bundesamts fü r Kartographie und Geodasie, 2011: 123-131.

Comparison of three modes of LS+AR model for predicting Earth rotation parameters

LEI Yu1，2，3，GAO Yu-ping1，2，CAI Hong-bing1，2
（1.National Time Service Center，Chinese Academy of Sciences，Xi′an 710600，China；2.Key Laboratory of Time and Frequency Primary Standards，National Time Service Center，Chinese Academy of Sciences，Xi′an 710600，China；3.University of Chinese Academy of Sciences，Beijing 100049，China）

In order to study the influences of the different modes of combined LS（least square）+AR（autoregression） model on the prediction of ERP（Earth rotation parameters），we adapted three modes，i.e.the recursion mode，the iteration mode and the interval mode，to predict ERP.The results show that the three modes have equivalent prediction accuracies in forecasting the length of day（LOD） for all prediction spans，while in forecasting the polar motion the recursion mode has the highest prediction accuracy，and the interval mode outperforms the iterative mode.In addition，the analyses indicate that the data-utilization rates for the recursionand iterative modes are higher than that for the interval mode，however，the interval mode shows the fastest computation speed compared with the other two modes，and the amount of calculation for recursion mode is less than that for iteration mode.

Earth rotation parameters（ERP）； prediction； LS+AR model； recursion mode； iteration mode；interval mode

P227

A

1674-0637（2015）01-0052-09

10.13875/j.issn.1674-0637.2015-01-0052-09

2014-05-27

中國科學(xué)院“西部之光”人才培養(yǎng)計劃“聯(lián)合學(xué)者”資助項目

雷雨，男，博士研究生，助理研究員，主要從事地球自轉(zhuǎn)變化監(jiān)測與預(yù)報研究。