基于ARIMA模型的地球自轉參數(shù)預報研究

楊　杰，葉修松，曾　光，朱　俊

(1.宇航動力學國家重點實驗室，西安　710043；2.西安衛(wèi)星測控中心，西安　710043)

?

基于ARIMA模型的地球自轉參數(shù)預報研究

楊杰1，2，葉修松1，2，曾光1，2，朱俊1，2

(1.宇航動力學國家重點實驗室，西安710043；2.西安衛(wèi)星測控中心，西安710043)

摘要：為了提高基于ARIMA模型的地球自轉參數(shù)預報精度，著重分析了ERP數(shù)據(jù)序列預處理中周期項成分的精確獲取技術、ARMA模型階次的優(yōu)化選擇問題，以及ERP數(shù)據(jù)序列不同差分次數(shù)對模型精度的影響。利用國際地球自轉和參考系統(tǒng)服務組織公布的ERP最終產(chǎn)品，365組預報結果表明，UT1-UTC預報30 d精度可達3.398 6 ms，10 d以內(nèi)預報應選取2次差分模型，10 d以上選取1次差分模型；日長預報30 d精度可達0.337 5 ms，應選取零次差分預報模型；Xp預報30 d精度為11.189 1 mas，Yp預報30 d精度為7.932 1 mas，兩者均選擇1次差分模型。

關鍵詞：差分自回歸滑動平均模型；地球自轉參數(shù)預報；差分；模型階次；優(yōu)化

0引言

地球定向參數(shù)是連接天球和地球參考框架的關鍵參數(shù)，是衛(wèi)星導航應用的關鍵因素，廣泛應用于衛(wèi)星、飛船及空間飛行器的地基測定軌、地面點或船只的定位、彈道導彈的軌道設計、目標的瞄準發(fā)射以及來犯導彈的攔截、空基定位的飛行器星下點坐標的測定和對地姿態(tài)的測控等方面。同時，地球定向參數(shù)也是天文地球動力學、天體測量學、地球物理學和大地測量學研究的最主要觀測量和基礎科學數(shù)據(jù)，為研究地球自轉變化的機制、板塊運動、固體潮及地球內(nèi)部結構等地球動力學問題提供了極大的便利，也是研究地球整體旋轉、表面物質運動、大氣環(huán)流變化甚至地質災害預報的基礎信息，極大推動了天文地球動力學的發(fā)展[1]。

20世紀60年代以來，甚長基線干涉、激光測衛(wèi)、全球定位系統(tǒng)等現(xiàn)代大地測量技術的快速發(fā)展，極大地提高了地球定向參數(shù)的觀測精度。以國際地球自轉和參考系統(tǒng)服務組織(International Earth Rotation and Reference Systems Service，IERRSS)為代表的國際組織機構定期發(fā)布精確的地球自轉參數(shù)(Earth rotation parameters，ERP)，但無法適應精密大地測量、航天測控工程等領域的實時應用，為此需要利用IERRSS的預報產(chǎn)品或基于IERRSS的解算數(shù)據(jù)開發(fā)自主的ERP預報產(chǎn)品。常見的ERP預報算法有最小二乘法、協(xié)方差分析法、神經(jīng)網(wǎng)絡法、自回歸滑動平均法、小波分析法、卡爾曼濾波法，以及其他組合預報法等[2-9]。

采用差分自回歸滑動平均模型(auto-regressive integrated moving average，ARIMA)對ERP參數(shù)進行外推預報。為使預處理后的ERP更符合自回歸滑動平均(auto-regressive moving average，ARMA)模型對隨機過程的要求，可從2個方面進行分析：①對數(shù)據(jù)序列進行精確的頻譜分析，獲取準確的周期項成分；②差分ERP數(shù)據(jù)序列以使其滿足隨機過程的要求。此外，為了提高ARMA模型的預報精度，需要選取合適的模型階次，可通過數(shù)據(jù)序列的自吻合誤差確定最優(yōu)的模型階次。在開發(fā)ERP預報產(chǎn)品時，可通過比較不同差分次數(shù)對測試樣本數(shù)據(jù)序列預報精度的影響以選擇合適的差分預報模型。

1ERP預報流程

圖1　UT1-UTC的預報算法

圖2　日長和極移的預報算法

ERP預報的主要流程見圖1和圖2[10]，主要包含預測模型的建立和外推估計，針對大氣角動量和海洋角動量建模復雜的特點，在預報過程中先不予考慮。在預測模型建立過程中，針對周日變化，需要減去跳秒、潮汐、UT2-UT1、趨勢項和周期項成分等，而日長(length of day，LOD)和極移參數(shù)預測模型的建立則相對簡單，僅需進行趨勢項和周期項成分的分離。針對預處理后的ERP，需要進一步進行差分運算，以獲取滿足ARMA模型特性的數(shù)據(jù)時間序列，再進行ARMA模型參數(shù)估計和外推預測。

為了使預處理后的ERP數(shù)據(jù)序列滿足ARMA模型的特性要求，需要充分剔除數(shù)據(jù)序列中的周期項。同時，可通過差分數(shù)據(jù)序列進一步提高ERP數(shù)據(jù)序列預處理的精度。圖1和圖2中的最小二乘模型表示為

(1)

式(1)中，ωi表示ERP數(shù)據(jù)序列中的周期項頻率，比如錢德勒周期項和年周期項等。

2ARIMA模型的預報算法

ARIMA模型預報算法的主要流程如圖3所示[11]

圖3　ARIMA模型預報算法

針對圖3中的預報流程，主要解決3個技術問題，即ERP數(shù)據(jù)序列的差分次數(shù)、ARMA模型的階次選擇、ERP數(shù)據(jù)序列不同差分次數(shù)的預報精度等。圖3中，各個環(huán)節(jié)的計算過程表示如下。

2.1樣本自相關和偏相關函數(shù)的計算

給定樣本序列{X1,X2,…,XN}， 自相關和互相關函數(shù)計算結果分別表示為

(2)

(3)

2.2ARIMA模型參數(shù)的矩估計

Step 1：自回歸參數(shù)向量φ的矩估計結果為

(4)

Step 2：向量Xt=Xt-φ1Xt-1-…-φpXt-p自協(xié)方差函數(shù)的矩估計結果表示為

(k=0,1,2,…,q)

(5)

Step3：滑動平均系數(shù)向量θ的矩估計結果為

(6)

2.3ARIMA模型不同差分階次預報算法

2.3.1直接預報

(7)

(8)

2.3.2差分預報

差分后的序列預報結果見式(6)和式(7)，那么，一次差分預報結果為

(9)

二次差分預報結果為

(10)

2.3.3ARMA(p，q)序列格林函數(shù)的計算

格林函數(shù)計算過程表示為

(11)

結合式(11)，針對不同的模型階次組合，格林函數(shù)計算過程分別表示為

1)ARMA(1，1)模型及其格林函數(shù)表示為

Xt-φ1Xt-1=at-θ1at-1

(12)

(13)

2)ARMA(1，2)模型及其格林函數(shù)表示為

Xt-φ1Xt-1=at-θ1at-1-θ2at-2

(14)

(15)

3)ARMA(1，3)模型及其格林函數(shù)表示為

Xt-φ1Xt-1=at-θ1at-1-θ2at-2-θ3at-3

(16)

(17)

4)ARMA(2，1)模型及其格林函數(shù)表示為

Xt-φ1Xt-1-φ2Xt-2=at-θ1at-1

(18)

(19)

5)ARMA(2，2)模型及其格林函數(shù)表示為

Xt-φ1Xt-1-φ2Xt-2=at-θ1at-1-θ2at-2

(20)

(21)

6)ARMA(2，3)模型及其格林函數(shù)表示為

Xt-φ1Xt-1-φ2Xt-2=at-θ1at-1-θ2at-2-θ3at-3

(22)

(23)

7)ARMA(3，1)模型及其格林函數(shù)表示為

Xt-φ1Xt-1-φ2Xt-2-φ3Xt-3=at-θ1at-1

(24)

(25)

8)ARMA(3，2)模型及其格林函數(shù)表示為

Xt-φ1Xt-1-φ2Xt-2-φ3Xt-3=at-θ1at-1-θ2at-2

(26)

(27)

9)ARMA(3，3)模型及其格林函數(shù)表示為

Xt-φ1Xt-1-φ2Xt-2-φ3Xt-3=

at-θ1at-1-θ2at-2-θ3at-3

(28)

(29)

3ERP數(shù)據(jù)序列頻譜分析結果

利用IERRSS最終公布的2000年～2012年的數(shù)據(jù)序列，通過快速傅立葉變換(fast Fourier transformation，F(xiàn)FT)變換對去除趨勢項后的各參數(shù)序列進行頻譜分析，求取ERP序列中的周期項頻譜分量，結果表示如下。

3.1UT1-UTC頻譜分析結果

UT1-UTC頻譜分析結果見圖4。

圖4　UT1-UTC頻譜分析結果

根據(jù)圖4獲取的UT1-UTC時間序列周期信號頻譜分量見表1。

表1　UT1-UTC序列周期項

根據(jù)表1可知，UT1-UTC序列中的周期分量不僅包含錢德勒周期、0.5 a、1 a周期項，還包括其它高頻低頻周期項。

3.2LOD頻譜分析結果

LOD頻譜分析結果見圖5。

圖5　LOD頻譜分析結果

根據(jù)圖5獲取的LOD時間序列周期信號頻譜分量見表2。

表2　LOD序列周期項

根據(jù)表2可知，LOD序列中的周期分量不僅包含錢德勒周期、0.5 a、1 a周期項，還包括其它高頻低頻周期項。

3.3Xp及Yp頻譜分析結果

Xp及Yp頻譜分析結果見圖6和圖7。

圖6　Xp頻譜分析結果

圖7　Yp頻譜分析結果

根據(jù)圖6和圖7獲取的極移序列周期信號頻譜分量見表3。

表3　Xp及Yp序列周期項

根據(jù)表3可知，極移序列中的周期分量主要為錢德勒周期和1 a周期項。

4ARMA模型階次優(yōu)化選擇及ERP預報精度分析

采用2000年～2012年的數(shù)據(jù)序列，訓練樣本共4 565個數(shù)據(jù)，采用滑動窗口的方式，每組訓練樣本數(shù)據(jù)序列分別預報30 d的ERP參數(shù)，共統(tǒng)計365組預報結果。每組預報誤差為30 d誤差絕對值的最大值。

每組訓練樣本在ARMA模型參數(shù)矩估計中，模型階次(p，q)依次取為(1，1)，(1，2)，(1，3)，(2，1)，(2，2)，(2，3)，(3，1)，(3，2)，(3，3)，選取每組訓練樣本自擬合誤差的最小偏差為評價準則進行模型階次優(yōu)化選擇。這里，僅給出前30組UT1-UTC數(shù)據(jù)序列的最優(yōu)階次，見表4。

表4　前30組UT1-UTC數(shù)據(jù)序列ARMA模型最優(yōu)階次

4.1UT1-UTC預報結果

零次差分預報結果較差。這里僅給出1次和2次差分預報30 d UT1-UTC的365組預報誤差分布和平均預報誤差，見圖8和圖9。

圖8　UT1-UTC差分預報誤差分布

圖9　UT1-UTC差分預報平均誤差

根據(jù)圖8和圖9可知，UT1-UTC的5 d最好預報結果為0.427 ms，10 d為1.117 ms，30 d為3.399 ms。當預報天數(shù)小于10 d時，二次差分精度略高；當預報天數(shù)大于10 d時，一次差分預報精度明顯高于二次差分。

4.2LOD預報結果

這里給出不同差分次數(shù)預報30 d LOD的365組預報誤差分布和平均預報誤差，見圖10和圖11。

圖10　LOD差分預報誤差分布

圖11　LOD差分預報平均誤差

根據(jù)圖10和圖11可知，LOD的5 d最好預報結果為0.134 ms，10 d為0.205 ms，30 d為0.338 ms。零次差分預報精度略優(yōu)于一次差分，明顯優(yōu)于二次差分。

4.3Xp預報結果

零次差分預報結果較差。這里僅給出1次和2次差分預報30 dXp的365組預報誤差分布和平均預報誤差，見圖12和圖13。

圖12　Xp差分預報誤差分布

圖13　Xp差分預報平均誤差

根據(jù)圖12和圖13可知，5 dXp最好的預報結果為2.287 mas，10 d為5.941 mas，30 d為11.189 mas。預報天數(shù)小于5 d時，2種差分模型預報精度相當；預報天數(shù)大于5 d時，零次差分預報精度較差，一次差分預報精度明顯優(yōu)于二次差分。

4.4Yp預報結果

零次差分預報結果較差。這里僅給出1次和2次差分預報30 dYp的365組預報誤差分布和平均預報誤差，見圖14和圖15。

圖14　Yp差分預報誤差分布

圖15　Yp差分預報平均誤差

根據(jù)圖14和圖15可知，5 dYp最好的預報結果為1.506 mas，10 d為2.765 mas，30 d為7.852 mas。預報天數(shù)小于5 d時，2種差分模型預報精度相當；預報天數(shù)大于5 d時，零次差分預報精度較差，一次差分預報精度明顯優(yōu)于二次差分。

5結束語

本文分析了ARIMA模型預報的3項關鍵技術。分別研究了ERP數(shù)據(jù)序列中周期項成分的精確獲取問題、ERP模型階次的優(yōu)化問題，以及不同差分次數(shù)對ERP預報精度的分析。為此，可得到如下結論：

1)訓練樣本數(shù)據(jù)序列的頻譜分析可準確獲取周期項成分，從而使得預處理后的ERP數(shù)據(jù)序列具有更符合ARMA模型要求的隨機過程特性。

2)訓練樣本數(shù)據(jù)序列的最小自擬合誤差均方差可建立模型階次優(yōu)化選擇的評價準則。

3)測試樣本數(shù)據(jù)序列的最小預報誤差有利于模型差分次數(shù)的優(yōu)化選擇，可用于不同天數(shù)ERP預報產(chǎn)品的開發(fā)。

4)通過上述分析方法可知，30 d UT1-UTC的預報精度為3.399 ms，30 d LOD預報精度為0.338 ms，30 dXp預報精度為11.189 mas，30 dYp預報精度為7.852 mas。

致謝：作者感謝國際GNSS監(jiān)測評估系統(tǒng)(iGMAS)提供的觀測數(shù)據(jù)(www.igmas.org)。

參考文獻

[1]鄭大偉，虞南華.地球自轉及其和地球物理現(xiàn)象的聯(lián)系：I日長變化[J].地球物理學進展，1996，11(2)：34-37.

[2]AKAIKE H.Utoregressive model fitting for control[J].Annals of the Institute of Statistical Mathematics，1971，23(1)：162-180.

[3]BARNES R T H，HIDE R，WHITE A A，et al.Atmospheric angular momentum fluctuations，length-of-day changes and polar motion[C]//Royal Society.Proceedings of Mathematical and Physical Sciences.London，UK：Royal Society，1983：31-73.

[4]HAMDAN K，SUNG L Y.Stochastic modeling of length of day and universal time[J].Journal of Geodesy，1996，70(1)：307-320.

[5]KOSEK W，MCCARTHY D，LUZUM B J.Possible improvement of Earth orientation forecast using autocovariance prediction procedures[J].Journal of Geodesy，1998，72(12)：189-199.

[6]GROSS R，EUBANKS T M，STEPPE J A.A Kalman filter-based approach to combining independent Earth-orientation series[J].Journal of Geodesy，1998，72(9)：215-235.

[7]SCHUH H，ULRICH M，EGGER D，et al.Prediction of Earth orientation parameters by artificial neural networks[J].Journal of Geodesy，2002，76(3)：247-258.

[8]JOHNSON T，LUZUM B J，RAY J R.Improved near-term Earth rotation predictions using atomospheric angular momentum analysis and forecasts[J].Journal of Geodesy，2005，79(6)：209-221.

[9]HUBER P J.Modeling the length of day and extrapolating the rotation of the Earth[J].Journal of Geodesy，2006，80(12)：283-303.

[10]NIEDZIELSKI T，KOSEK W.Prediction of UT1-UTC，LOD and AAM χ3 by combination of least-squares and multivariate stochastic methods[J].Journal of Geodesy，2008，82(2)：83-92.

[11]BROCKWELL P J，DAVIS R A.Introduction to time series and forecasting[M].New York，USA：Springer，1979：212-220.

Research of the ARIMA-based Prediction of Earth Rotation Parameters

YANGJie1，2，YEXiusong1，2，ZENGGuang1，2，ZHUJun1，2

(1.State Key Laboratory of Astronautic Dynamics，Xi’an 710043，China；2.Xi’an Satellite Control Center，Xi’an 710043，China)

Abstract：Some critical problems are successfully resolved to improve the Earth rotation parameters (ERP) prediction accuracy based on the auto-regressive integrated moving average (ARIMA) model，which are the precise extraction of periodic components in the ERP，the order optimization of ARIMA model and the effects of model order on the ERP prediction accuracy.The prediction experiments are implemented on the final ERP products released by the International Earth Rotation and Reference Systems Service (IERRSS).Some key conclusions can be drawn from the total 365 groups of prediction results.Firstly，the 30-day prediction error of UT1-UTC is 3.398 6 ms.The two-order difference model can be selected for the prediction in less than 10 d，but the one-order difference model is selected for the prediction in more than 10 d.Secondly，the 30-day prediction error of length of day (LOD) is 0.337 5 ms and the difference model order is zero.Finally，the 30-day prediction error is 11.189 1 mas and 7.932 1 mas for the x-axis and y-axis component of polar motion respectively.Additionally，the difference model order is one for the above two cases.

Key words：auto-regressive integrated moving average model；Earth rotation parameters prediction；difference；model order；optimization

中圖分類號：P228

文獻標識碼：A

文章編號：2095-4999(2016)-01-0068-08

作者簡介：第一楊杰(1985—)，男，山西運城人，工程師，博士，現(xiàn)主要從事航天器精密定軌和定位方面工作。

收稿日期：2015-06-15

引文格式：楊杰，葉修松，曾光，等.基于ARIMA模型的地球自轉參數(shù)預報研究[J].導航定位學報，2016，4(1)：68-74,87.(YANG Jie，YE Xiusong，ZENG Guang，et al.Research of the ARIMA-based Prediction of Earth Rotation Parameters[J].Journal of Navigation and Positioning，2016，4(1)：68-74,87.)DOI：10.16547/j.cnki.10-1096.20160114.