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

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

1 引言

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

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

2 長自回歸白噪化方法

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

長自回歸白噪化方法的步驟為：

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

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

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

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

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

（3）對于 ARMA（p，q）序列：

有

記

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

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

3 計算策略

3.1 日長變化的預(yù)報

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

圖1 日長變化預(yù)報的流程圖

3.2 極移的預(yù)報

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

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

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

4 計算與分析

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

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

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

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

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

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

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

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

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

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

5 結(jié)論與分析

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

［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）模型在極移實時預(yù)測中的應(yīng)用［J］，大地測量與地球動力學(xué)，2010，30（1）：88－91.

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

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

［7］ 楊叔子.時間序列分析的工程應(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.