采用自回歸滑動(dòng)平均模型的地球自轉(zhuǎn)參數(shù)短期預(yù)報(bào)

葉修松，何雨帆，曾 光，郭 海

（1.宇航動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710043；2.中國(guó)西安衛(wèi)星測(cè)控中心，西安 710043）

1 引言

由于地球自轉(zhuǎn)參數(shù)（Earth rotation parameters，ERP）的解算過(guò)程比較復(fù)雜，甚長(zhǎng)基線干涉測(cè)量（very long baseline interferometry，VLBI）、衛(wèi)星激光測(cè)距（satellite laser ranging，SLR）等高精度觀測(cè)手段往往需要幾天甚至更長(zhǎng)的時(shí)間才能獲取其變化值；全球定位系統(tǒng)（global positioning system，GPS）盡管解算速度較快（約3h），但GPS解算有較大的系統(tǒng)影響，長(zhǎng)期穩(wěn)定性不好，需要用VLBI和SLR解來(lái)進(jìn)行修正［1］。因此，這些參數(shù)實(shí)際上是無(wú)法實(shí)時(shí)獲得的。由于極移和日長(zhǎng)的變化受眾多激發(fā)源的影響，其變化不僅具有一定的線性趨勢(shì)及多種周期變化趨勢(shì)，同時(shí)還包括了復(fù)雜的不規(guī)則變化特性，這給預(yù)測(cè)工作帶來(lái)了極大的困難。往往是預(yù)測(cè)幾天之后，其預(yù)測(cè)誤差便遠(yuǎn)遠(yuǎn)大于其觀測(cè)誤差。因此，對(duì)ERP的預(yù)測(cè)算法進(jìn)行研究，建立更為合理、有效的預(yù)測(cè)模型，提高地球自轉(zhuǎn)參數(shù)的預(yù)報(bào)精度，具有十分重要的意義［2－4］。

為此，本文擬采用自回歸滑動(dòng)平均模型（autoregressive moving average，ARMA）對(duì)ERP進(jìn)行短期預(yù)報(bào)。該算法在于它不受地球自轉(zhuǎn)動(dòng)力學(xué)模型等因素的制約，具有計(jì)算速度快、算法穩(wěn)定可靠等優(yōu)點(diǎn)。

2 長(zhǎng)自回歸白噪化方法

ARMA（p，q）序列參數(shù)估計(jì)的較為有效且簡(jiǎn)便的方法是 “長(zhǎng)自回歸白噪化方法”。長(zhǎng)自回歸白噪化方法的基本思想是認(rèn)為隨機(jī)序列｛Xt｝ 可以用階數(shù)較低的自回歸模型擬合，擬合的殘差作為白噪聲的估計(jì)，再用最小二乘法的原理，估計(jì)出參數(shù)＝ （φ1，…，φp）和＝ （θ1，…，θp）。這一方法的特點(diǎn)是全部求解過(guò)程都是解線性方程組，避免了非線性運(yùn)算。

長(zhǎng)自回歸白噪化方法的步驟為：

設(shè)有ARMA（p，q）序列 ｛Xt｝的一段樣本X1，X2，…，XN。

（1）先用AR（pN）模型擬合一段樣本，pN取適當(dāng)?shù)牡碗A，AR（pN）的參數(shù)是φ1，…，φpN。按照AR（p）序列參數(shù)的矩估計(jì)，由樣本值X1，X2，…，XN，可以求出φ1，…，φpN的矩估計(jì)量，…，［5－7］。

（2）由于 AR（pN）模型為

將，…代入此式左端，就得到αt的估計(jì)值：

如果AR（pN）擬合的好，就應(yīng)當(dāng)是白噪聲序列，如果擬合的不好，就不是白噪聲序列。是否為白噪聲序列，可以直接用 “周期圖法”進(jìn)行檢驗(yàn)。如果不是白噪聲序列，就要改變pN的值，直到｝ 為白噪聲序列為止。

（3）對(duì)于 ARMA（p，q）序列：

有

記

在S（）的表達(dá)式中，Xt，Xt－1，…，Xt－p是已知的樣本值，αt－1，…，αt－q是由（4）式算得的αt的估計(jì)值，求和符號(hào)Σ的下標(biāo)取為p＋q＋1是由于要計(jì)算出αt－q必需要求t≥p＋q＋1。

按最小二乘法的原理來(lái)確定未知參數(shù)φ1，…，φp，φ1，…，φp的估計(jì)值，…，，，…，。即要求

3 計(jì)算策略

3.1 日長(zhǎng)變化的預(yù)報(bào)

日長(zhǎng)（length of day，LOD）與 UT1－UTC是可以相互轉(zhuǎn)換的，即扣除UT1－UTC的跳秒后，再作差分，結(jié)果的相反數(shù)即為L(zhǎng)OD。UT1為世界時(shí)（universal time）、UTC為協(xié)調(diào)世界時(shí)（coordinated universal time）、TAI為原子時(shí)（international atomic time）、UT1R－TAI為從 UT1－TAI序列中減去潮汐分量的部分為，UT2R－TAI為從UT1RTAI序列中減去地球自轉(zhuǎn)的周期性季節(jié)變化的部分。因此，對(duì)LOD的預(yù)報(bào)可以轉(zhuǎn)換為對(duì)UT1－UTC的預(yù)報(bào)［1］。UT1－UTC的預(yù)報(bào)流程如圖1所示［8］。

圖1 日長(zhǎng)變化預(yù)報(bào)的流程圖

3.2 極移的預(yù)報(bào)

極移變化的預(yù)報(bào)算法如圖2所示，包括建立預(yù)測(cè)模型和外推預(yù)測(cè)兩個(gè)部分。在建立預(yù)測(cè)建模時(shí)，用最小二乘法計(jì)算出趨勢(shì)項(xiàng)和周期項(xiàng)的擬合系數(shù)，計(jì)算公式為［9－10］

對(duì)于極移殘差序列的預(yù)報(bào)也采用ARMA算法實(shí)現(xiàn)，極移序列ARMA預(yù)報(bào)的原理類同于日長(zhǎng)變化的預(yù)報(bào)。

圖2 極移序列預(yù)報(bào)的流程圖

4 計(jì)算與分析

計(jì)算與分析的數(shù)據(jù)來(lái)源于地球自轉(zhuǎn)服務(wù)（international Earth rotation service，IERS）官方網(wǎng)站公布的ERP時(shí)間序列。該序列包含1962－01－01至今的地球定向參數(shù)值，采樣間隔為1d。實(shí)際計(jì)算選取6a極移時(shí)間作為基礎(chǔ)序列長(zhǎng)度，因?yàn)?a接近Chandler周期與周年的最小公倍數(shù)，這樣更能體現(xiàn)地球自轉(zhuǎn)參數(shù)的周期特征。采用ARMA模型要求時(shí)間序列必須為平穩(wěn)的隨機(jī)序列，通過(guò)上述的ERP計(jì)算策略分析發(fā)現(xiàn)，僅對(duì)ERP時(shí)間序列扣除周期項(xiàng)和趨勢(shì)項(xiàng)之后所得的ERP殘差子時(shí)間序列隨機(jī)序列的波動(dòng)性較明顯、平穩(wěn)性較差，而基于殘差作二次差分后的子時(shí)間序列的波動(dòng)更趨于平穩(wěn)，更接近零均值。因此，本文基于自主編寫的ERP短期預(yù)報(bào)軟件分析了ERP殘差子時(shí)間序列與殘差作二次差分后的子時(shí)間序列對(duì)ARMA模型的適用性。

圖3 PX殘差序列預(yù)報(bào)7d的結(jié)果

圖4 PX二次差分后殘差序列預(yù)報(bào)7d的結(jié)果

圖5 PY殘差序列預(yù)報(bào)7d的結(jié)果

圖6 PY二次差分后殘差序列預(yù)報(bào)7d的結(jié)果

圖7 UT1－UTC殘差序列預(yù)報(bào)7d的結(jié)果

圖8 UT1－UTC二次差分后殘差序列預(yù)報(bào)7d的結(jié)果

從圖3～圖8可以得出，基于殘差二次差分后的預(yù)報(bào)結(jié)果明顯要優(yōu)于直接使用殘差預(yù)報(bào)的結(jié)果，其主要原因是剔除周期項(xiàng)和趨勢(shì)項(xiàng)的殘差（包括一次殘差差分）序列沒(méi)有明顯的規(guī)律特征，其平穩(wěn)性還不能滿足平穩(wěn)隨機(jī)序列的條件。而基于殘差二次差分后的殘差子時(shí)間序列的波動(dòng)更趨于平穩(wěn)，更接近零均值。

表1 基于ARMA算法的ERP 短期預(yù)報(bào)結(jié)果

從表1可知，直接采用殘差時(shí)間序列預(yù)報(bào)1d的ERP精度約為：PX＝10mas，PY＝0.9mas，UT1－UTC＝0.5ms；采用殘差殘二次差分后的時(shí)間序列預(yù)報(bào)1d的ERP精度約為：PX＝0.1ms，PY＝0.4mas，UT1－UTC＝0.003ms。綜上所述，本文采用的基于非線性自回歸的ERP短期預(yù)報(bào)精度與IERS公布的結(jié)果相當(dāng)。

5 結(jié)論與分析

本文采用線性自回歸方法對(duì)地球自轉(zhuǎn)參數(shù)進(jìn)行短期預(yù)報(bào)，并與EOPPCC的結(jié)果進(jìn)行比較，以驗(yàn)證預(yù)報(bào)模型的有效性。本文提出的基于ARMA的ERP短期預(yù)報(bào)技術(shù)已達(dá)到國(guó)際水平。本文還通過(guò)試驗(yàn)發(fā)現(xiàn)，基于ARMA方法進(jìn)行ERP時(shí)間序列預(yù)報(bào)時(shí)，要求ERP時(shí)間序列不能過(guò)長(zhǎng)或過(guò)短，其主要原因是數(shù)據(jù)精度不一致，從而影響預(yù)報(bào)精度。本文同時(shí)還驗(yàn)證了文獻(xiàn) ［2］提出的采用6a極移時(shí)間作為基礎(chǔ)序列長(zhǎng)度的有效性。另外，本文僅采用ERP時(shí)間序列本身來(lái)進(jìn)行預(yù)報(bào)，沒(méi)有考慮與之有極強(qiáng)相關(guān)性的大氣角動(dòng)量，這也是本文的不足之處。

［1］ ARTZ T，BERNHARD L，NOTHNAGEL A，et al.Methodology for the Combination of Sub－daily Earth Rotation from GPS and VLBI Observations［J］.Journal of Geodesy，2012，86（3）：221－239.

［2］ KOSEK W，MCCARTHY D D，LUZUM B J.Possible Improvement of Earth Orientation Forecast Using Auto－covariance Prediction Procedures［J］.Journal of Geodesy，1998，72（4）：189－199.

［3］ 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.

［4］ 徐君毅.基于截距修正的 GM（1，1）模型在極移實(shí)時(shí)預(yù)測(cè)中的應(yīng)用［J］，大地測(cè)量與地球動(dòng)力學(xué)，2010，30（1）：88－91.

［5］ 魏武雄，WILLIAM W S.時(shí)間序列分析—單變量和多變量方法［M］.易丹輝，劉超，賀學(xué)強(qiáng)，等，譯.2版.北京：中國(guó)人民大學(xué)出版社，2009.

［6］ 潘國(guó)榮.基于時(shí)間序列分析的動(dòng)態(tài)變形預(yù)測(cè)模型研究［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2005，30（6）：483－487.

［7］ 楊叔子.時(shí)間序列分析的工程應(yīng)用［M］.2版.武漢：華中科技大學(xué)出版社，2007.

［8］ NIEDZIELSKI T，KOSEK W.Prediction of UT1－UTC，LOD and AAMχ3by Combination of Least－squares and Multivariate Stochastic Methods［J］.Journal of Geodesy.2008，82（2）：83－92.

［9］ THALLER D，KRüGEL M，ROTHACHER M，et al.Combined Earth Orientation Parameters Based on Homogeneous and Continuous VLBI and GPS Data［J］.Journal of Geodesy.2007，81（6－8）：529－541.

［10］ CERVEIRA P J M.Earth Rotation Observed by Very Long Baseline Interferometry and Ring Laser［J］.Pure and Applied Geophysics，2009，166（8－9）：1499－1517.