采用自回歸滑動平均模型的地球自轉參數短期預報

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

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

1 引言

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

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

2 長自回歸白噪化方法

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

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

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

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

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

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

如果AR（pN）擬合的好，就應當是白噪聲序列，如果擬合的不好，就不是白噪聲序列。是否為白噪聲序列，可以直接用 “周期圖法”進行檢驗。如果不是白噪聲序列，就要改變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。

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

3 計算策略

3.1 日長變化的預報

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

圖1 日長變化預報的流程圖

3.2 極移的預報

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

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

圖2 極移序列預報的流程圖

4 計算與分析

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

圖3 PX殘差序列預報7d的結果

圖4 PX二次差分后殘差序列預報7d的結果

圖5 PY殘差序列預報7d的結果

圖6 PY二次差分后殘差序列預報7d的結果

圖7 UT1－UTC殘差序列預報7d的結果

圖8 UT1－UTC二次差分后殘差序列預報7d的結果

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

表1 基于ARMA算法的ERP 短期預報結果

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

5 結論與分析

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

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

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

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

［7］ 楊叔子.時間序列分析的工程應用［M］.2版.武漢：華中科技大學出版社，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.