改進的LS+AR模型在極移短期預測中的應用研究

張軍杰

（上海市巖土地質研究院有限公司，上海 200072）

0.引言

極移是指地球的極點在地球表面發生有規律的緩慢變化的現象，極移參數作為地固坐標系和天球坐標系的轉換參數之一，在深空探測、衛星定軌和激光測月等領域中均有著廣泛地應用[1]。但是由于復雜的數據處理過程，得到的最終極移參數存在著時間上的滯后性，無法滿足某些要求實時或準實時極移參數的領域[2]。尤其是隨著北斗二代衛星導航系統的逐步拓展，對極移參數的時效性和精度的需求日益迫切[3]，因此，如何提高極移參數的短期預報精度成為亟需解決的問題。

近年來，國內外學者對于極移參數預報做了大量有益的研究工作，提出了很多非常實用的模型。神經網絡模型（Artificial Neural Network,ANN）是常見的模型之一。人工神經網絡是由大量處理單元互聯組成的非線性、自適應信息處理系統，王琪潔等人研究和探索了應用非線性的神經網絡技術，將大氣角動量時間序列引入到地球自轉變化預報中，改善地球自轉參數（ERP）的預報精度，以及應用神經網絡技術預報El Nino南方濤動（ENSO）事件，結果證實了神經網絡具有良好的非線性預報能力。張志、廖瑛等人利用傅里葉分析法，獲得了插值基礎序列的周期，驗證了極移參數基礎序列重采樣的可行性，然后提取插值后的基礎序列數據的趨勢項，利用多輸入—單輸出反向傳播（Back Propagation，BP）神經網絡進行建模，預報不同跨度的殘差序列項,最后合并趨勢項和殘差序列得到最終的極移參數預報。

最小二乘和自回歸預報模型（Least squares and autoregressive,LS+AR）[4,5]也較常用。雷雨等人分別利用遞推、迭代和間隔3種預報方法對地球自轉參數進行預報。預報結果顯示，這3種方式對日長變化（length of day，LOD）參數的所有跨度預報的精度相當,遞推方式在極移所有跨度的預報結果上精度最高，迭代方式預報結果最差。姚宜斌等人通過對LS+AR模型短期預報殘差的時間序列統計分析，發現臨期的模型預報殘差具有極高的相關性，由此提出利用上一期模型預報殘差和經驗調節矩陣對下一期預報結果進行修正，從而提高了極移參數的超短期和短期預報精度。張昊、王琪潔等人利用最小二乘外推與ARIMA（P，1，0）組合模型對極移參數進行預報，ARIMA的全稱為求和自回歸移動模型，是d階差分運算（對原始時間序列進行差分的階數）和ARIMA（p,q）模型的結合，p為AR模型的階數，q為WA模型的階數。ARIMA（p,d,q）進行預報的思路是首先對原始極移參數時間序列進行d階差分，然后運用ARMA（p,q）對差分項進行預報，最后逆運算得到預報值結果。最小二乘外推與ARIMA(P,1,0)組合模型對于提高超短期極移參數預報精度有一定幫助。王小輝等人利用LS-SVM模型，將赤道向大氣角動量（AAM）和赤道向海洋角動量（OAM）時間序列引入到極移參數序列的預報中，改善了極移參數的預報精度。張昊根據最小二乘模型中錢德勒周期項的時變特性，提出了基于錢德勒參數時變修正的CLS模型,并利用CLS+AR模型對極移參數序列進行預測，CLS+AR模型在極移參數預報精度上較LS+AR模型有所改善。趙丹寧等人利用經驗模態分解對LS+AR模型進行改進，首先利用經驗模態分解方法對極移序列進行分解,獲得極移的高頻分量和低頻分量；然后采用最小二乘外推模型對極移低頻分量進行擬合，獲得最小二乘擬合殘差；之后利用自回歸模型對極移高頻分量和最小二乘擬合殘差之和進行建模預報；最后將最小二乘模型和自回歸模型外推值相加獲得極移參數的預報值。結果表明可以有效改善極移參數的預報精度。為檢驗各種模型的極移預報精度，維也納理工大學舉行了全球性的EOP預報對比運動（Earth Orientation Parameters Prediction Comparison Campaign，EOP PCC），確定了目前預報精度最高的模型是諧波最小二乘外推+自回歸預報模型，即LS+AR模型。

本文將基于國際上公認為最優的LS+AR模型，對其進行改進，以期進一步提高極移參數的預報精度，達到國際先進水平。

1.LS+AR模型

1.1 LS模型

自1900年國際緯度局利用目視天頂儀觀測極移以來，天文工作者對極移已進行了近百年的研究。根據長期的觀測和研究發現，極移參數存在多個基本的周期項，如，錢德勒項、周年項和半周年項等周期項。在利用最小二乘方法對極移參數序列進行建模時，首先要確定其固定周期項和趨勢項，在本文中包含趨勢項、錢德勒項、周年項和半周年項等四項。LS模型的具體數學表達式如式（1）、式（2）所示：

式（1）中，a0為極移長趨勢項；t為極移參數對應的時間，（單位為年）；a1為趨勢項系數；a2、a3為錢德勒項系數，對應的T1為錢德勒周期項的周期時間1.183a；a4、a5為周年項系數；對應的T2為周年項周期1a；a6、a7為半周年項的系數，對應的T3為半周年項周期0.5a。極移Y分量對應參數含義與式（1）中的一致。

LS模型的參數具體解算方法如式（3）～式（6）：

式（3）～式（6）中，X為最小二乘的模型參數矩陣；B為模型參數的系數矩陣；L為X分量的觀測值序列組成的矩陣。極移Y分量的計算方法與X分量的算法一致。

1.2 AR模型

AR模型是依據自身過去的時間序列的變化規律來預測下一時刻的值，要求建模的時間序列是零均值的平穩隨機序列。所以建模之前需要對LS模型擬合后的殘差序列進行預處理，通常是先進行一階差分后再進行建模。最終的極移預報值為LS模型的外推值與AR模型預報值之和。AR模型的數學表達式可表示為隨機序列zt（t＝1，2，…，n）與t時刻以前的規律性變化和t時刻的白噪聲的關系：

式（7）中，η1，η2，…，ηp為模型參數；at為白噪聲序列；p為模型階數；式（7）中zt稱為p階自回歸模型，簡記為AR（P）。

建立AR模型就必須首先確定模型的階，然后再確定AR模型參數。定階方法主要有3種，分別是最終預測誤差準則、信息論準則、傳遞函數準則。理論上這三種方法是等效的，本文采用最終預測誤差準則來確定AR模型的階。

式（9）中，PM為用AR模型擬合zt（t=1，2，3，…，n）序列殘差的均方差，當FPE（M）為最小值時取M作為AR模型的階。定階也可以利用Matlab高階統計量工具箱中的函數來解決，計算十分方便。

1.3 精度評定

為了便于評估預測精度，本文采用國際上通用的平均絕對誤差（Mean Absolute Error,MAE）作為精度評定標準，數學公式如式（10）所示：

式（10）中，n為預報次數；k為預報長度；εi，j為預報值與真值的差值。

2.改進的LS+AR模型

極移是用來描述地極在地球表面位置隨時間而變化的現象，其自身激發機制十分復雜，包括太陽輻射、月球引力、海洋潮汐和地下水分布，還有季節的變化引起的降雪變化和大氣流動變化等因素。極移現象激發因素的不確定性，導致極移分量的周期項和趨勢項具有時變性的特征。為了提高極移參數的預報精度，本文在獲得基本的極移參數后，利用基本的趨勢項和德勒項、周年項和半周年項等周期項進行初步擬合，在獲得一次擬合后的殘差項后，利用頻譜分析探測殘差項中隱含的殘差周期項，此時殘差項依然包含比較明顯的未知趨勢項和周期項，利用LS模型建模對殘差項再次進行擬合，獲得第一次殘差項的趨勢項和周期項，然后對二次殘差項進行一階差分，對差分結果進行AR建模，然后利用AR模型對差分項進行預報，最后根據LS模型得到的兩次趨勢項和周期項以及AR模型預報的殘差差分項獲得最終的預報極移參數。具體流程（如圖1所示）：

圖1 極移參數預報流程

3.基礎序列的選取

極移運動成因復雜，不同長度的基礎序列對短期預報結果的精度也是不同的，本文首先以2a、4a、6a、8a和10a作為基礎序列長度對極移參數進行預報，尋找更加合適的基礎序列長度。基礎序列長度并非越長越好，主要是因為隨著時間間隔的增加，地球季節性變化、海洋潮汐、月球引力、太陽輻射等各種偶然物理影響因素也會隨之增加，極移序列之間的相關性會逐步下降。因此在利用LS+AR模型對極移參數進行建模時，基礎序列長度不宜過長，這樣不僅可以提高預報效率，還可以保障預報精度的穩定性。對比基礎序列長度與預報精度的關系（如圖2、圖3所示），結果顯示當基礎序列為4a時預報精度最佳，所以本次預報采用4a長度作為預報極移參數的基礎序列。

圖2 Xpole參數2、4、6、8年的預報精度

圖3 Ypole參數2、4、6、8年的預報精度

4.預報結果與精度分析

本次預報所用到的數據來自國際地球自轉和參考系服務組織（IERS）所發布的EOP 08C04序列（https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html）。該序列包含從1962年1月1日至今的極移參數的Xpole、Ypole分量，時間間隔為1天。本文自2005年1月1日（MJD：55317）起開始預報，每7天預報一次，一共預報了350期，預報長度為50天。為進一步削弱一次擬合后的殘差項里隱含周期項對預報精度的影響，首先用LS模型對極移參數進行擬合，采用傅里葉分析法對一次擬合后的殘差項時間序列進行頻譜分析，目的是將時間域中的隱藏波形轉變為頻率域的頻譜，顯示其隱藏周期。每140天對極移X和Y分量的殘差項的時間序列去趨勢項和功率譜分析，對極移參數第一次擬合后的殘差結果的頻譜分析值（如表1所示），在18次的頻譜分析值中，每次頻譜分析的周期值變化差異較大，說明了十分有必要對其進行二次擬合，削弱明顯周期項的影響。

表1 擬合殘差的頻譜分析周期表

利用頻譜分析得到的殘差周期再次對極移殘差序列進行擬合，圖4和圖5中藍色為一次擬合后的殘差序列，紅色為兩次擬合后的殘差序列。可以明顯看出極移分量兩次擬合后的殘差時間序列更加趨向于平穩、零均值，不包含明顯的周期項。然后對得到的兩次擬合后的殘差項進行一階差分，利用AR模型對差分后的極移殘差序列進行建模，得到一階差分項的預報值，最后將兩次擬合值和一階差分值預報值逆運算結果相加，獲得最終極移參數預報值。

圖4 Xpole參數第二次擬合前后殘差對比

圖5 Ypole參數第二次擬合前后殘差對比

為便于比較最小二乘和自回歸預報（LS+AR）模型改進后的預報精度，本文將采用三種預報方式進行短期預報。方案一：最小二乘和自回歸預報模型（LS+AR）；方案二：首先利用最小二乘進行擬合，然后對殘差項進行一階差分，最后對差分項進行自回歸預報（LS+一階差分+AR）；方案三：首先利用最小二乘進行擬合，然后對殘差項進行頻譜分析，獲得一次擬合后的殘差周期項，再次進行擬合，最后對二次擬合后的殘差項進行一階差分，利用自回歸模型對差分項進行預報（LS+LS（頻譜分析）+一階差分+AR）。最終試驗的預報結果（如圖6、圖7所示）：

圖6 Xpole參數預報精度比較

圖7 Ypole參數預報精度比較

由圖6和圖7可以看出：前10天三種預報方案預報精度差別不大，10～50天之間方案2和方案3的極移參數預報精度比方案1的極移參數預報精度有著明顯提高，再次驗證了AR模型要求建模的時間序列必須是平穩的隨機序列。方案2和方案3的極移參數預報精度整體比較接近，圖8是方案3較之方案2的精度提高百分比，可以明顯看出極移X分量的預報精度提高較為明顯，極移Y分量除了在3～6日有少許下降外，其他時間段預報精度均有明顯提升。由此得出改進后的LS+AR模型是三種方案里預報精度最高的，驗證了本文所提出改進方法的可行性。

圖8 精度提高百分比

5.結束語

本文考慮到極運動周期時變性的特征，對LS+AR預報模型進行改進，對首次擬合后的極移參數殘差序列進行頻譜分析，利用頻譜分析得到的周期再次進行擬合，然后再對殘差項進行一階差分，最后對一階差分項進行自回歸預報，有效地提高了極移參數的預報精度。該方法僅從數學角度進行極移參數短期預報，并未分析頻譜分析得到周期的實質的物理意義。研究極移序列潛在周期的具體激發源，并進一步改進將是下一步的工作重點。此外，本文推薦以4年作為極移基礎序列進行短期預報較為合適。