LS+AR模型預報極移的一種修正算法

郭忠臣，孫 朋

（宿州學院環(huán)境與測繪工程學院，安徽宿州234000）

地球在自轉過程中，自轉軸的位置相對于地球本身會發(fā)生變化，進而引發(fā)地極運動，這種現(xiàn)象稱為極移。極移是天文研究、大地測量、衛(wèi)星導航等領域不可缺少的參數(shù)[1-3]。當前主要有全球導航衛(wèi)星系統(tǒng)（GNSS）、衛(wèi)星激光測距（SLR）、甚長基線干涉測量（VLBI）等技術來測定極移，但由于高精度極移參數(shù)解算復雜，不能實時獲取，所以只能通過預報的方式獲得高精度極移數(shù)據[4-7]。

極移預報模型較多，主要有最小二乘外推（LS）模型、最小二乘外推和自回歸組合（LS+AR）模型、神經網絡模型和聯(lián)合角動量模型等[8-12]。但由地球自轉參數(shù)預測方案比較大會戰(zhàn)（Earth Orientation Parameters Prediction Comparison Campaign,EOP PCC）可知，沒有一種模型適合所有間隔的預報，但LS+AR模型對短期預報效果較優(yōu)[13]。

利用LS+AR模型對極移進行短期預報時，將極移分解為趨勢項、周期項和殘差項，分別進行預報，各項預報之和即為極移預報值。由于地球上物質分布的不均勻性，導致了極移部分周期項的振幅和相位變化具有一定的無規(guī)律性[14]，同時考慮到Kalman濾波的諸多優(yōu)點，本文將其引入到LS+AR組合模型中對預報模型進行修正，結果表明，修正后的模型預報精度優(yōu)于LS+AR模型。

1 極移預報模型及修正策略

LS+AR模型是EOP PCC活動所認可的極移短期預報較優(yōu)的方法，首先使用LS模型對極移序列的趨勢項及周期項進行擬合，與原始序列對比得到殘差項，然后再使用AR模型對殘差項進行建模預報，兩者預報之和即為最終預報結果[15]。

1.1 LS模型

極移受到較多激發(fā)源的影響，對極移序列進行分析可知，其主要由趨勢項和周期項組成，趨勢項滿足線性變化，周期項主要受Chandler項、周年項和半周年項影響，故對極移序列可構建LS模型：

其中，ai、bi為趨勢項擬合系數(shù)，ci、di、ei為周期項的振幅，P1、P2、P3分別為Chandler項、周年項和半周年項的周期，t為時間。

式（1）可簡化為

其中，L表示極移序列值，B表示系數(shù)矩陣，X表示待求參數(shù)。

B、X分別表示為

利用最小二乘計算公式即可求得參數(shù)X：

1.2 AR模型

AR模型建模簡單，主要利用序列內部的相關性進行建模預報，其基本形式可表示為

其中，z為原始序列，φi為AR模型系數(shù)，εt為t時刻的白噪聲，P為模型階數(shù)。

AR模型預報精度受模型階數(shù)p的影響較大，一般情況下，模型階數(shù)越高，模型包含的信息越多，預報精度越好，但當階數(shù)增加到一定程度時，會出現(xiàn)較多無用信息導致預報精度下降，因此需要一種合理的方法來確定模型階數(shù)。

常用定階方法有AIC準則和BIC準則，AIC準則和BIC準則均在模型中引入了與模型參數(shù)個數(shù)相關的懲罰項，以達到降低模型參數(shù)、避免過擬合的目的。但是AIC準則在樣本數(shù)據量較大時，通常會失效，因為似然函數(shù)值較大，超過了模型參數(shù)的影響，而在BIC準則中加大了與模型參數(shù)個數(shù)相關的懲罰力度，在樣本數(shù)量較大時，可有效控制模型的復雜度，即模型參數(shù)的個數(shù)。因此，本文采用BIC準則確定模型階數(shù)。

BIC準則定義為[16]：

其中，p為模型參數(shù)個數(shù)，n為樣本量，L為似然函數(shù)。

1.3 Kalman濾波修正策略

Kalman濾波計算簡單，可對包含隨機誤差的線性動態(tài)系統(tǒng)進行最優(yōu)估計[17]，本文將Kalman濾波帶入到AR模型中，對AR模型預報結果進行修正。

由Kalman濾波狀態(tài)預測方程和公式（4）可得：

其中，Xt表示t時刻的狀態(tài)向量，A為狀態(tài)轉移矩陣，?表示系統(tǒng)噪聲向量，z表示觀測向量，，ε表示觀測噪聲向量。

由LS外推之后得到的殘差項為平穩(wěn)的隨機序列，且觀測噪聲向量和系統(tǒng)噪聲向量均為白噪聲，故狀態(tài)轉移矩陣A和噪聲向量的方差矩陣D?、Dε均可定義成單位陣。

狀態(tài)向量最優(yōu)估值的計算過程可表示如下。

（1）計算誤差矩陣預測方差：

（2）計算濾波增益：

（3）狀態(tài)向量最優(yōu)估值計算：

（4）誤差矩陣更新：

公式（6）～（10）表示了Kalman濾波對AR模型參數(shù)求解過程的修正方式，利用修正之后的參數(shù)即可對LS模型的殘差項實時構建AR模型進行預報，該模型記為LS+AR+KF模型。

1.4 精度評價指標

本文采用平均絕對誤差（Mean Absolute Error,MAE）對極移短期預報精度進行評價，其計算公式為

其中，i為預報跨度，j=1,2，…，n為預報期數(shù)，P為極移預報值，O為極移觀測值。

2 實驗分析

本文實驗所用數(shù)據來源于國際地球自轉服務（International Earth Rotation Service,IERS）發(fā)布的EOP 08C04序列（http://hpiers.obspm.fr/eoppc/eop/eopc04/eopc04.62-now），序列包含從1962年至今、數(shù)據采樣間隔為1天的全部極移數(shù)據。考慮到極移測定精度的不斷提高，選取1998年01月01日至2018年01月01日期間的極移序列作為實驗數(shù)據，以10年長度的極移序列作為原始序列建模預報之后30天的極移數(shù)據，每期預報間隔30天，共預報121期，采用MAE對預報精度進行評價。為驗證加入Kalman濾波修正之后的預報精度，將其與LS+AR模型及EOP PCC活動預報精度對比，圖1給出了EOP PCC活動極移的短期預報結果，圖2和圖3給出了LS+AR模型和LS+AR+KF模型的預報結果。

通過分析可知：

（1）LS+AR+KF模型預報精度優(yōu)于LS+AR模型，且隨著預報時間的增加，精度提高明顯；

（2）EOP PCC活動對極移X方向（PMX）預報1天的精度多在0.5～1.0mas內，對極移Y方向（PMY）預報1天的精度多在0.4～1.0mas內，LS+AR模型與LS+AR+KF模型對PMX和PMY預報1天的精度分別為0.283mas和0.281mas，精度明顯優(yōu)于EOP PCC活動結果，原因在于所用建模數(shù)據的精度較高，且建模序列長度選用10年的，建模精度較高；

（3）EOP PCC活動預報的PMX在10天跨度內，精度多在3.3～5.0mas內，預報的PMY在10天跨度內，精度多在1.6～3.6mas內，LS+AR模型與LS+AR+KF模型對PMX和PMY預報10天的精度分別為 3.353mas、3.287mas和 2.176mas、2.035mas，LS+AR+KF模型預報精度優(yōu)于EOP PCC活動中的大多數(shù)方法；

（4）EOP PCC活動預報的PMX在30天跨度內，精度多在8.0～14.0mas內，預報的PMY在30天跨度內，精度多在4.0～12.0mas內，LS+AR模型與LS+AR+KF模型對PMX和PMY預報30天的精度分別為11.514mas、8.527mas 和7.818mas、5.094mas，LS+AR+KF模型預報精度優(yōu)于EOP PCC活動中的大多數(shù)方法，且LS+AR+KF模型與LS+AR模型相比較，預報30天的PMX和PMY精度分別提高25.94%和34.85%。

圖1 EOP PCC活動預報結果

圖2 兩種方案預報10天精度對比

圖3 兩種方案預報30天精度對比

3 總結

針對極移序列具有時變性的特點，在傳統(tǒng)LS+AR模型中引入Kalman濾波對預報結果進行修正，實驗結果表明，添加Kalman濾波之后的修正模型較LS+AR模型預報精度有所提高，且隨著預報時間的增加，精度提高明顯，并且優(yōu)于EOP PCC活動中大部分模型的預報精度。