考慮LS+AR模型基礎(chǔ)數(shù)據(jù)量的極移預(yù)報研究*

徐海龍 喬書波 王穆陽

(1 信息工程大學(xué)地理空間信息學(xué)院 鄭州 450001)

(2 31121部隊 南京 210018)

1 引言

地球定向參數(shù)(Earth Orientation Parameter,EOP)是實現(xiàn)天球參考框架和地球參考框架轉(zhuǎn)換的必要參數(shù), 對于高精度導(dǎo)航定位、衛(wèi)星精密定軌和深空探測等方面有著重要意義[1]. 地球定向參數(shù)的觀測模式、數(shù)據(jù)處理和解算較為復(fù)雜, 導(dǎo)致EOP的解算結(jié)果通常具有幾小時至數(shù)天的延遲, 用戶不能得到實時EOP數(shù)據(jù), 因此, EOP的短期至中長期預(yù)報數(shù)據(jù)在衛(wèi)星自主導(dǎo)航和天文動力學(xué)研究中有重要的應(yīng)用[2]. 地球定向參數(shù)由極移(Polar Motion,PM)、日長變化、歲差和章動組成, 極移的兩個分量又可以分別用PMX和PMY表示. EOP中歲差和章動模型IAU2000/IAU2006擬合精度較高, 可以實現(xiàn)長期高精度預(yù)測[3]. 所以, 當(dāng)前對地球定向參數(shù)的預(yù)報研究主要是極移的兩個分量和日長變化,這兩者又被統(tǒng)稱為地球自轉(zhuǎn)參數(shù)(Earth Rotation Parameter, ERP).

國內(nèi)外學(xué)者相繼提出了多種預(yù)報方法, 如最小二乘(Least Square, LS)外推和自回歸(AutoRegressive, AR)組合模型[4–5]、LS外推和神經(jīng)網(wǎng)絡(luò)(Artificial Neural Network, ANN)組合模型[6–8]、灰色模型[9]、模糊推理系統(tǒng)[10]、自協(xié)方差模型等以及這些模型的組合及改正算法[11–12]. 為了比較驗證當(dāng)前各種預(yù)報方法的精度和穩(wěn)定性, 進一步提高地球定向參數(shù)的預(yù)報精度, 維也納理工大學(xué)依托國際地球自轉(zhuǎn)與參考系服務(wù)(International Earth Rotation and Reference Systems Service, IERS)在2005年到2008年組織了地球定向參數(shù)預(yù)報比較活動(Earth Orientation Parameter Prediction Comparison Campaign, EOP PCC). 此次活動在全球范圍內(nèi)有8個國家共12個相關(guān)領(lǐng)域頂尖專家或小組參與, 涉及20種預(yù)報方法, 最終結(jié)論表明對于極移序列預(yù)報精度最高的為LS+AR模型[3].

基于LS+AR模型, 相關(guān)學(xué)者在模型改進和數(shù)據(jù)選取方面展開了廣泛研究, 進一步提高了模型預(yù)測精度.模型改進方面,姚宜斌等[13]通過對LS+AR模型短期預(yù)報殘差的相關(guān)性分析, 提出了利用預(yù)報殘差和調(diào)節(jié)矩陣修正結(jié)果的新模型; 加松等[14]利用Kalman濾波對AR模型修正、利用最小均方誤差自適應(yīng)濾波對LS擬合外推項修正, 兩種修正模型的預(yù)報精度均較LS+AR模型有較大提升. 張昊等[15]、Sun等[16]和孫張振[17]針對極移近期數(shù)據(jù)與未來走勢相關(guān)性更強的特點, 并結(jié)合數(shù)據(jù)解算精度, 通過不同方法設(shè)計了LS擬合時的權(quán)矩陣,有效提高了極移預(yù)報精度. 數(shù)據(jù)選取方面, 隨著EOP的測定精度不斷提高和數(shù)據(jù)長度的增加, 相關(guān)學(xué)者也注意到不同的基礎(chǔ)序列長度影響極移預(yù)報結(jié)果的精度. 姚宜斌等[13]在研究附加誤差修正的LS+AR新模型時, 選用了6、9、12和15 yr 4個不同長度的極移序列進行預(yù)報, 結(jié)果顯示在進行10 d和30 d預(yù)報時, 采用6 yr的極移序列精度優(yōu)于其他長度序列預(yù)報結(jié)果; 嚴(yán)鳳等[18]通過實驗對比發(fā)現(xiàn),對于1–10 d的ERP超短期預(yù)報, 基礎(chǔ)序列越長反而精度越差, 6 yr是超短期預(yù)報的最優(yōu)基礎(chǔ)序列長度;Xu等[19]采用LS+AR模型進行地球自轉(zhuǎn)參數(shù)預(yù)報時, 通過計算不同長度輸入序列和對應(yīng)預(yù)報跨度關(guān)系, 證明了對于預(yù)報不同跨度時需要輸入不同的最優(yōu)基礎(chǔ)序列長度.

上述針對LS+AR模型的研究均是將該組合模型看作一個整體, 數(shù)據(jù)選取時只是變化不同長度的基礎(chǔ)序列, 將預(yù)報結(jié)果對比后確定LS+AR模型的單一最優(yōu)基礎(chǔ)長度. 然而這種數(shù)據(jù)選取方法并不理想, 如果數(shù)據(jù)量太小則LS模型不能很好地提取周期信息, 如果數(shù)據(jù)量太大又容易導(dǎo)致輸入AR模型的殘差數(shù)據(jù)冗余, 使模型過擬合造成模型偏差較大從而影響精度. Wu等[20]對輸入AR模型的數(shù)據(jù)量進行了討論, 結(jié)果表明采用不同大小的LS殘差數(shù)據(jù)量對AR模型的預(yù)報精度影響較大, 但文獻中LS基礎(chǔ)序列采用的是固定10 yr長度進行實驗, 并不能代表LS+AR模型整體輸入數(shù)據(jù)序列的最佳長度. 本文在實驗中考慮LS模型和AR模型各自擬合需求的不同, 分別計算確定用于極移預(yù)報的最優(yōu)組合數(shù)據(jù)長度, 將最終預(yù)報結(jié)果同只確定一個基礎(chǔ)序列長度的預(yù)報結(jié)果對比, 并在最后給出結(jié)果分析.

2 預(yù)報模型與原理

使用LS+AR模型預(yù)報時, 首先用LS模型對極移序列的線性趨勢項和主要周期項進行擬合并外推預(yù)測, 然后使用AR模型對LS擬合殘差序列進行建模預(yù)報, 最后將LS外推值和AR模型預(yù)報值相加即為最終預(yù)報結(jié)果.

2.1 LS模型

時間序列中經(jīng)常出現(xiàn)線性趨勢項和諧波振蕩周期項, 極移序列中同樣存在趨勢項和周期項. 該確定性部分可以利用下式的最小二乘模型進行擬合外推, 并得到LS擬合殘差項.

式中,t為時間,X(t)表示t時刻的極移數(shù)值,a0和a1表示序列線性趨勢項參數(shù),i為擬合時第i個周期項,ci和di表示序列的各周期項參數(shù),Ti表示周期時長,n表示周期項個數(shù). (1)式中前兩項為線性部分,可對極移序列的線性趨勢項擬合, 后一項由三角函數(shù)構(gòu)成, 可對極移序列的周年、半周年和錢德勒項等周期項進行擬合.

2.2 AR模型

對極移序列確定性部分進行建模和預(yù)報后, 剩余殘差項采用自回歸AR模型擬合預(yù)報. AR模型是一種根據(jù)時間序列變量自身過去的規(guī)律來建立的預(yù)報模型[11], 其數(shù)學(xué)模型為:

式中,xt為殘差序列,φj為AR模型的第j個參數(shù),εt為擬合后剩余噪聲序列,p為模型階數(shù). (2)式又可稱為p階自回歸模型, 簡稱為AR(p).

使用AR模型進行預(yù)報關(guān)鍵在于模型階數(shù)的確定, 常用的定階準(zhǔn)則有最小信息量準(zhǔn)則(Akaike Information Criterion, AIC)、最終預(yù)報誤差準(zhǔn)則(Final Prediction Error, FPE)和矩陣奇異值分解(Singular Value Decomposition, SVD)等. 本文在實驗中使用AIC準(zhǔn)則進行模型定階. AIC準(zhǔn)則的思想是通過兩個方面去衡量擬合模型的優(yōu)劣[21]: 綜合考慮擬合程度的似然函數(shù)值和擬合模型中未知參數(shù)個數(shù). 通常來說, 似然函數(shù)值越大就代表該模型擬合效果越好. 模型未知參數(shù)個數(shù)越多代表模型靈活度較高、模型擬合準(zhǔn)確度較高. 而從另一個方面看, 模型參數(shù)越多也代表參數(shù)估計難度加大, 容易導(dǎo)致模型精度變差. 所以如何發(fā)現(xiàn)未知參數(shù)和擬合精度之間的最優(yōu)配置, 對最優(yōu)模型至關(guān)重要.AIC準(zhǔn)則即是擬合精度和未知參數(shù)的加權(quán)函數(shù), 使AIC函數(shù)值達到最小的模型也稱為最優(yōu)模型.

通過AIC準(zhǔn)則確定AR模型階次p, 便可進一步利用時間序列估計模型中未知參數(shù). 最經(jīng)典的方法是利用最小二乘估計對AR(p)模型參數(shù)求解[21], 首先對于(2)式中的待求參數(shù)φ1,φ2,···,φp可以用表示, 即:

則AR(p)模型擬合值Ft()為:

此時, 殘差項如(5)式所示. 最小二乘原理指出, 使殘差平方和Qt(β)達到最小的參數(shù)值即為參數(shù)~β的估計值, 可如下計算:

該解算策略充分利用了時間序列的信息, 因此最小二乘估計的精度很高. 利用此過程得到AR(p)模型參數(shù), 最后根據(jù)下列方程進行序列多步預(yù)報:

其中,N表示模型輸入基礎(chǔ)序列長度,xN為第N天的殘差序列值,L為AR模型的預(yù)報跨度,, ···,為AR模型預(yù)報結(jié)果.

3 實驗分析

3.1 精度評價標(biāo)準(zhǔn)

本文采用EOP PCC推薦的平均絕對誤差(Mean Absolute Error, MAE)作為最終精度評價標(biāo)準(zhǔn), 并增加絕對誤差(Absolute Error, AE)來進一步驗證對比預(yù)報穩(wěn)定性, 其公式分別如下:

式中, prer和obsr分別為極移PMX和PMY序列第r天的預(yù)報值和真實值,u為預(yù)報跨度,Np為預(yù)報天數(shù), MAE表示平均絕對誤差, AE表示預(yù)報值與真實值絕對差值.

3.2 實驗及結(jié)果

極移預(yù)報的基礎(chǔ)序列大多來自IERS發(fā)布的EOP C04序列, 該序列為1962年以來間隔1 d的地球定向參數(shù). 本文選用當(dāng)前最新的EOP 14C04序列中的極移數(shù)據(jù)進行實驗(https://hpiers.obspm.fr/iers/eop/eopc04). 預(yù)報起始時間為2018年1月第1周, 每周預(yù)報1期, 每期預(yù)報360 d, 至2019年12月底結(jié)束預(yù)報共計104期. 實驗時的基礎(chǔ)數(shù)據(jù)分別采用原始方式的單一數(shù)據(jù)長度和本文的組合數(shù)據(jù)方式.單一數(shù)據(jù)分別選擇3、6、9、12、15、18、21和24 yr作為實驗基礎(chǔ)序列長度, 利用LS+AR模型進行預(yù)報, 以便得出模型最優(yōu)單一基礎(chǔ)序列長度. 本文組合數(shù)據(jù)同樣分別選擇3、6、9、12、15、18、21和24 yr作為實驗預(yù)報序列長度, 首先利用LS模型進行擬合外推并得到殘差序列, 然后重新選擇1/5、2/5、3/5、4/5長度的殘差序列輸入AR模型進行預(yù)報, 兩者預(yù)報值相加即為最終預(yù)報值, 即可對比選出適合模型的最優(yōu)組合基礎(chǔ)序列長度.

圖1為極移PMX和PMY分量使用單一基礎(chǔ)序列長度時的實驗結(jié)果. 結(jié)果顯示除1–10 d的超短期預(yù)報精度比較接近外, 隨著預(yù)報跨度的增加, 預(yù)報時輸入基礎(chǔ)數(shù)據(jù)長度不同, 對應(yīng)預(yù)報精度也不同.對于PMX和PMY分量, 最優(yōu)單一基礎(chǔ)數(shù)據(jù)長度均為6 yr, 這一結(jié)論與姚宜斌等[13]、嚴(yán)鳳等[18]的工作一致. 分析發(fā)現(xiàn), 這主要是因為6 yr正好是極移序列約1.2 yr的錢德勒項周期和周年項周期的最小公倍數(shù), 一定程度上體現(xiàn)了極移的周期特征; 另一方面, AR模型作為LS+AR模型的一部分, 其處理的是LS模型的擬合殘差序列, 如果數(shù)據(jù)量太少, 可能引起AR模型擬合時參數(shù)解算方差過大; 而過多的數(shù)據(jù)可能造成AR模型過擬合, 導(dǎo)致模型偏差較大,對最終預(yù)報結(jié)果產(chǎn)生影響. LS模型和AR模型之間這種微妙的制約關(guān)系, 使得6 yr成為適合LS+AR模型的單一基礎(chǔ)序列長度. 但這種簡單使用6 yr長度的LS+AR模型實驗, 并不能代表模型的最適合數(shù)據(jù)量, 所以仍需要進行組合數(shù)據(jù)實驗驗證才能得出分別適用于LS模型和AR模型的數(shù)據(jù)長度.

圖2和圖3分別為極移PMX和PMY分量預(yù)報結(jié)果. 當(dāng)實驗選擇不同的基礎(chǔ)數(shù)據(jù)長度時, 對應(yīng)的AR模型所輸入的數(shù)據(jù)也不同, 實驗中各分量分別產(chǎn)生32種結(jié)果, 其中當(dāng)LS擬合基礎(chǔ)數(shù)據(jù)長度為21 yr, AR擬合數(shù)據(jù)量為1/5殘差總長時, PMX和PMY分量預(yù)報結(jié)果精度均最高. 進一步分析可以發(fā)現(xiàn), 當(dāng)基礎(chǔ)數(shù)據(jù)長度為21 yr時, 預(yù)報精度隨著AR模型數(shù)據(jù)量比例增大而減小, 這可能是因為采用過長或過久的殘差序列造成AR模型存在過擬合, 導(dǎo)致模型偏差較大而影響預(yù)報精度.

圖2 PMX分量不同組合數(shù)據(jù)預(yù)報精度. 其中圖(a)、(b)、(c)和(d)分別表示AR模型輸入數(shù)據(jù)量分別為LS擬合殘差序列長度的1/5、2/5、3/5和4/5.Fig.2 The prediction accuracies of combined data of PMX. In the figure, panels (a), (b), (c) and (d) indicate that the input data amount of AR model is 1/5, 2/5, 3/5 and 4/5 length of the residual sequence of LS fitting, respectively.

圖3 PMY分量不同組合數(shù)據(jù)預(yù)報精度. 其中圖(a)、(b)、(c)和(d)分別表示AR模型輸入數(shù)據(jù)量分別為LS擬合殘差序列長度的1/5、2/5、3/5和4/5.Fig.3 The prediction accuracies of combined data of PMY. In the figure, panels (a), (b), (c) and (d) indicate that the input data amount of AR model is 1/5, 2/5, 3/5 and 4/5 length of the residual sequence of LS fitting, respectively.

將組合數(shù)據(jù)模式實驗中的預(yù)報最優(yōu)結(jié)果與圖1中最優(yōu)單一數(shù)據(jù)為6 yr時極移預(yù)報結(jié)果對比. 統(tǒng)計結(jié)果如表1所示, 表中“6”即代表單一基礎(chǔ)數(shù)據(jù)長度為6 yr, “21+1/5”代表組合長度為LS選擇21 yr、AR模型選擇1/5殘差總長的最優(yōu)組合數(shù)據(jù). 表1所示為單一數(shù)據(jù)和組合數(shù)據(jù)分別進行極移兩個分量預(yù)報時的MAE精度對比.對于PMX分量,組合數(shù)據(jù)結(jié)果在1–360 d的預(yù)報跨度上比單一數(shù)據(jù)預(yù)報結(jié)果具有更高的精度, 從50 d開始, 組合數(shù)據(jù)相對單一數(shù)據(jù)的預(yù)報精度提升百分比達到14%, 最大可以達到28%, 第100 d后穩(wěn)定在20%以上. 對于PMY分量,在預(yù)報初期組合數(shù)據(jù)預(yù)報結(jié)果精度略低于單一數(shù)據(jù), 但兩者相差不大. 隨著預(yù)報跨度的增加, 至第30 d開始精度逐漸優(yōu)于單一數(shù)據(jù)結(jié)果, 精度提高比例最高時也可以達到23%. 結(jié)果表明, 對于1–360 d跨度的極移序列預(yù)報, 本文采用的組合數(shù)據(jù)預(yù)報相比單一數(shù)據(jù)預(yù)報精度提高較大.

表1 兩種數(shù)據(jù)方案預(yù)報精度統(tǒng)計(MAE: mas)Table 1 Statistic of the two kinds of data accuracies (MAE: mas)

圖1 PMX分量(a)和PMY分量(b)單一數(shù)據(jù)預(yù)報精度Fig.1 The prediction accuracies of single data of PMX (a) and PMY (b)

從預(yù)報跨度來看, 短期1–30 d的預(yù)報跨度內(nèi)提升效果不佳, 甚至PMY分量的組合數(shù)據(jù)預(yù)報結(jié)果誤差略有增加, 但隨著預(yù)報時長的增加, 組合數(shù)據(jù)預(yù)報優(yōu)勢逐漸顯現(xiàn), 較單一數(shù)據(jù)結(jié)果精度提升效果明顯. 這是因為在短期30 d內(nèi)時, 極移整體變化幅度不是很大, 數(shù)據(jù)量的改變對預(yù)報影響還不是特別明顯. 隨著預(yù)報跨度的增加, 當(dāng)進行至30–200 d左右的預(yù)報時, 組合數(shù)據(jù)預(yù)報精度穩(wěn)步提升, 在這一階段, 組合數(shù)據(jù)長度的優(yōu)勢使得LS模型和AR模型各自的預(yù)報精度都在提高, 所以預(yù)報精度提高比例也在上升. 而當(dāng)預(yù)報跨度增加到200–360 d, 預(yù)報精度提高比例趨于穩(wěn)定. 這可能是因為隨著預(yù)報跨度的增加, AR模型殘差部分預(yù)報值所發(fā)揮的作用逐漸減小, 而21 yr的LS擬合基礎(chǔ)序列, 包含信息更加豐富, 更適合中長期預(yù)報, 此時預(yù)報結(jié)果主要由趨勢項和周期項預(yù)報值主導(dǎo)所致.

為進一步研究證明組合數(shù)據(jù)結(jié)構(gòu)有利于極移預(yù)報精度的提高, 分別計算預(yù)報時段內(nèi)每期絕對誤差A(yù)E, 統(tǒng)計如圖4和圖5所示. 通過分析104期預(yù)報絕對誤差可以看出, 組合數(shù)據(jù)預(yù)報結(jié)果的絕對誤差相比單一數(shù)據(jù)也有了明顯改善, 而且當(dāng)預(yù)報跨度較大和預(yù)報誤差較大時改善效果更加明顯. 這是由分別計算決定LS+AR模型數(shù)據(jù)量的優(yōu)勢決定的. 使用21 yr的長期數(shù)據(jù)進行LS擬合, 可以充分提取序列中的周期信息, 而AR模型使用少部分的殘差序列即避免了冗長數(shù)據(jù)對AR模型參數(shù)求解的影響,也很好地利用了臨近數(shù)據(jù)的相關(guān)性, 這也進一步證明了本文數(shù)據(jù)選取方法在理論上的合理性.

圖4 PMX分量單一數(shù)據(jù)(a)和組合數(shù)據(jù)(b)絕對誤差統(tǒng)計Fig.4 Statistic of absolute error of single data (a) and combined data (b) of PMX

圖5 PMY分量單一數(shù)據(jù)(a)和組合數(shù)據(jù)(b)絕對誤差統(tǒng)計Fig.5 Statistic of absolute error of single data (a) and combined data (b) of PMY

4 結(jié)語

本文通過對LS+AR模型分別采取單一數(shù)據(jù)量和組合數(shù)據(jù)量的預(yù)報結(jié)果統(tǒng)計分析, 發(fā)現(xiàn)數(shù)據(jù)量的不同對于極移預(yù)報結(jié)果影響較大. 通過實驗分析,顯示單一數(shù)據(jù)序列在使用6 yr長度作為極移最優(yōu)數(shù)據(jù)時預(yù)報結(jié)果最好, 組合數(shù)據(jù)預(yù)報在使用21 yr長度的LS擬合序列和1/5殘差長度的AR擬合序列時預(yù)報精度最高. 對這兩種不同數(shù)據(jù)選取方式的預(yù)報結(jié)果進行對比, 結(jié)果表明, 在進行1–30 d的極移短期預(yù)報時, 改變輸入數(shù)據(jù)的選取方式對預(yù)報精度影響較小, 但在進行30 d以上中長期預(yù)報時, 組合數(shù)據(jù)預(yù)報精度相較于單一數(shù)據(jù)預(yù)報結(jié)果具有較大的提升, 其中PMX分量精度提升最高為28%, PMY分量最高為23%. 建議在極移分量的預(yù)報中, 可以采用組合數(shù)據(jù)預(yù)報, 以提高當(dāng)前極移的預(yù)報精度. 需要注意的是, 本文實驗只是探索驗證了組合數(shù)據(jù)的可行性, 所得出結(jié)果也只是數(shù)據(jù)選取的一個大概合適范圍, 但給出的結(jié)果仍可以證明組合數(shù)據(jù)預(yù)報結(jié)果在360 d的預(yù)報跨度內(nèi)精度更高.

致謝感謝參與本文評審的各位專家對文章提出的寶貴建議, 使得文章的質(zhì)量有了顯著的提高. 感謝IERS提供地球定向參數(shù)數(shù)據(jù).