孫張振，徐天河，范朋飛

(1．長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪工程學(xué)院，陜西西安710054;2．地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西西安710054;3．西安測(cè)繪研究所，陜西西安710054)

一、引 言

地球的自轉(zhuǎn)運(yùn)動(dòng)表征著地球整體的運(yùn)動(dòng)狀況，由地球各圈層地核、地幔地殼和大氣之間的相互作用引起的，可用地球定向參數(shù)(earth orientation para meters，EOP)描述[1]。EOP 包括 3 部分：① 章動(dòng)和歲差;② 極移;③ 日長(zhǎng)變化。其中，極移是表征地球運(yùn)動(dòng)的重要參數(shù)，它用來(lái)描述地球瞬時(shí)自轉(zhuǎn)軸在地球本體的運(yùn)動(dòng)而導(dǎo)致地球極點(diǎn)在地球表面上的位置發(fā)生緩慢變化的現(xiàn)象。高精度的極移觀測(cè)和預(yù)測(cè)數(shù)據(jù)對(duì)于天球參考框架和地球參考框架具有重要意義，是提高空間導(dǎo)航定位技術(shù)精度的必要部分。針對(duì)極移的預(yù)報(bào)，國(guó)內(nèi)外許多學(xué)者建立了各種預(yù)報(bào)模型對(duì)其進(jìn)行研究、預(yù)報(bào)和監(jiān)測(cè)[2-8]。總的來(lái)說(shuō)，這些模型大部分是將極移分為有規(guī)律的確定性部分和無(wú)規(guī)律的隨機(jī)性部分，即首先利用最小二乘方法(least squares，LS)來(lái)擬合極移時(shí)間序列中的趨勢(shì)項(xiàng)，這些有規(guī)律的確定性部分包括極移的長(zhǎng)期趨勢(shì)(LT)、錢德勒擺動(dòng)(CW)、周年擺動(dòng)(AW)和半周年擺動(dòng)(SAW)，并對(duì)其進(jìn)行外推預(yù)測(cè);然后利用AR 模型[6]、

ARMA模型或者人工神經(jīng)網(wǎng)絡(luò)[9]等確定性或者隨機(jī)性模型對(duì)其殘差部分進(jìn)行預(yù)測(cè)，最后得到預(yù)報(bào)的極移結(jié)果。在上述預(yù)測(cè)方法中，利用最小二乘方法擬合其趨勢(shì)項(xiàng)時(shí)，國(guó)內(nèi)外學(xué)者在周期項(xiàng)的具體數(shù)值上雖大體一致但普遍不同，國(guó)外學(xué)者在錢德勒周期中選擇在433～435 d之間較多。國(guó)內(nèi)學(xué)者如許雪晴在極移預(yù)報(bào)時(shí)選取錢德勒周期約為435 d[10]，張昊等錢德勒周期選取為 1．183 a(432．08 d)[11]，在周年項(xiàng)和半周年項(xiàng)上國(guó)內(nèi)外選取比較一致，都為1 a(365．24 d)和0．5 a(182．62 d)。由于在極移觀測(cè)數(shù)據(jù)中其周期項(xiàng)和趨勢(shì)項(xiàng)都有時(shí)變性的特點(diǎn)[12]，為了更加精確地進(jìn)行極移預(yù)報(bào)，本次在LS+AR預(yù)報(bào)模型的基礎(chǔ)上，在利用最小二乘擬合數(shù)據(jù)趨勢(shì)項(xiàng)之前，首先對(duì)62 a來(lái)的極移數(shù)據(jù)做頻譜分析[13]，由計(jì)算結(jié)果顯示其周期項(xiàng)中錢德勒周期為435．17 d，周年項(xiàng)為365．54 d，半周年項(xiàng)峰值對(duì)于X、Y分量有兩個(gè)峰值，卻有微小不同，鑒于其波峰值較小，故仍統(tǒng)一采用182．62 d。因此，此次在進(jìn)行LS擬合時(shí)，應(yīng)用頻譜分析得出的周期項(xiàng)參數(shù)，聯(lián)合AR模型預(yù)報(bào)極移并與其他學(xué)者所用周期項(xiàng)預(yù)報(bào)結(jié)果相比較。計(jì)算結(jié)果表明，應(yīng)用本文所選取的周期項(xiàng)進(jìn)行預(yù)報(bào)時(shí)，其預(yù)報(bào)精度最高。

二、LS+AR模型

LS+AR模型為一種組合的模型，在對(duì)極移進(jìn)行預(yù)報(bào)時(shí)，首先應(yīng)用最小二乘擬合趨勢(shì)項(xiàng)，扣除后得到殘差部分，再利用AR模型對(duì)殘差部分進(jìn)行預(yù)測(cè)，結(jié)合最小二乘外推部分即可得到預(yù)測(cè)值。

1．LS模型

在本文中，利用最小二乘外推模型中包含有線性項(xiàng)、半周年項(xiàng)、周年項(xiàng)和錢德勒項(xiàng)，模型公式為

式中，a0為常數(shù)項(xiàng);a1為趨勢(shì)項(xiàng);B1、B2、C1、C2、D1、D2…為各周期項(xiàng)擬合系數(shù);R1、R2、R3…分別為周期項(xiàng)中的周期;t為UTC時(shí)間。

2．AR模型

式(2)中，φ1，φ2，…，φp為模型參數(shù);at為白噪聲序列;p為模型階數(shù);式(2)稱為p階自回歸模型，簡(jiǎn)記為AR(p)。

AR模型要求建模的時(shí)間序列是平穩(wěn)隨機(jī)序列，即滿足平穩(wěn)、正態(tài)、零均值的條件。運(yùn)用AR模型的關(guān)鍵是確定AR模型的階數(shù)p，一般采用最終誤差準(zhǔn)則、信息論準(zhǔn)則、和傳遞函數(shù)準(zhǔn)則，實(shí)際中這3種方法是等效的。本文采取最終誤差準(zhǔn)則來(lái)確定AR模型的階數(shù)p：

最終誤差準(zhǔn)則

其中

式(4)中使FPE(p)達(dá)到極小值的p就取作AR模型的階數(shù)。φ1，φ2，…，φp等模型參數(shù)可以通過(guò)求解Yule-Walker方程來(lái)確定。

3．精度評(píng)定

為了便于評(píng)估預(yù)測(cè)精度，采用平均絕對(duì)誤差(mean absolute error，MAE)作為精度評(píng)定標(biāo)準(zhǔn)，其計(jì)算公式為

式中，Pi為第i期的預(yù)報(bào)值;Oi為相應(yīng)的觀測(cè)值;Ei為跨度第i期預(yù)報(bào)的真誤差(假定觀測(cè)值為真值);n為總預(yù)報(bào)期數(shù);MAEj為預(yù)報(bào)跨度為j的平均絕對(duì)誤差。

三、數(shù)據(jù)計(jì)算和結(jié)果分析

1．?dāng)?shù)據(jù)說(shuō)明

本次計(jì)算所用到的數(shù)據(jù)來(lái)自國(guó)際地球自轉(zhuǎn)和參考系服務(wù)IERS發(fā)布的EOP 05C04序列(http：∥hpiers．obspm．fr/eoppc/eop/)，該序列包含1962 年1 月1日至今的極移分量Xp和Yp、日長(zhǎng)變化LOD、地球自轉(zhuǎn)速率UT1-UTC等，時(shí)間間隔為1 d。首先對(duì)極移分量Xp和Yp做頻譜分析，圖1為極移頻譜周期分布圖，從圖上可明顯看出極移中錢德勒周期約在435 d和周年周期項(xiàng)約為365 d。結(jié)合實(shí)際計(jì)算數(shù)據(jù)，可知其中錢德勒周期為435．17 d，周年項(xiàng)為365．54 d。

圖1 極移頻譜分析圖

2．?dāng)?shù)據(jù)處理

為了確定極移中具體周期對(duì)極移預(yù)報(bào)精度的影響，本文此次在進(jìn)行最小二乘擬合時(shí)根據(jù)以往所采用的周期經(jīng)驗(yàn)值和本文中所確定的周期值分別對(duì)極移進(jìn)行預(yù)報(bào)，設(shè)計(jì)方案如下。

方案一：周期項(xiàng)選取錢德勒周期項(xiàng)432．08 d，周年項(xiàng)365．24 d，半周年項(xiàng)182．62 d。

方案二：周期項(xiàng)選取錢德勒周期項(xiàng)433 d，周年項(xiàng)365．24 d，半周年項(xiàng)182．62 d。

方案三：周期項(xiàng)選取錢德勒周期項(xiàng)435 d，周年項(xiàng)365．24 d，半周年項(xiàng)182．62 d。

方案四：周期項(xiàng)選取錢德勒周期項(xiàng)435．17 d，周年項(xiàng)365．24 d，半周年項(xiàng)182．62 d。

方案五：周期項(xiàng)選取錢德勒周期項(xiàng)435．17 d，周年項(xiàng)365．54 d，半周年項(xiàng)182．62 d。

在對(duì)極移(PM)進(jìn)行預(yù)報(bào)時(shí)，其觀測(cè)基礎(chǔ)序列長(zhǎng)度對(duì)預(yù)報(bào)精度有很大的影響，為了提高預(yù)報(bào)精度，本次預(yù)報(bào)所應(yīng)用的都是最優(yōu)基礎(chǔ)序列長(zhǎng)度，極移分量X和Y基礎(chǔ)序列長(zhǎng)度都為10 a觀測(cè)序列。預(yù)報(bào)精度結(jié)果統(tǒng)計(jì)是從2004年1月1日至2012年4月8日，分別對(duì)其進(jìn)行跨度為1～360 d的預(yù)報(bào)。圖2～4給出了以上5種方案的預(yù)報(bào)精度對(duì)比圖。

從圖2可以看出，錢德勒周期的具體數(shù)據(jù)選取對(duì)極移的預(yù)報(bào)精度有很大關(guān)系，雖然其變化只在3 d內(nèi)，但越是接近錢德勒周期真值，其極移預(yù)報(bào)就會(huì)越精確。從圖3中可知，方案三和方案四預(yù)報(bào)精度非常接近，方案四更精確一些，也說(shuō)明本文所用的錢德勒周期數(shù)據(jù)更加接近于真值。而從圖4可以看出，方案五的預(yù)報(bào)精度要高于方案四，這也可以說(shuō)明在極移數(shù)據(jù)中，周年項(xiàng)的具體真值并不是只局限于365．24 d，本文中所應(yīng)用的周年項(xiàng)數(shù)值365．54 d更接近與真值。

圖2 方案一、二、三預(yù)報(bào)精度對(duì)比

圖3 方案三、四預(yù)報(bào)精度對(duì)比

圖4 方案四、五預(yù)報(bào)精度對(duì)比

四、結(jié)束語(yǔ)

本文基于LS+AR預(yù)測(cè)模型對(duì)極移進(jìn)行高精度的預(yù)報(bào)，由于在極移數(shù)據(jù)中其周年項(xiàng)、錢德勒項(xiàng)等周期項(xiàng)具有時(shí)變性的特點(diǎn)，并不是只局限于某些具體的數(shù)值，故在預(yù)報(bào)之前首先對(duì)極移數(shù)據(jù)做頻譜分析，確定出所應(yīng)用觀測(cè)數(shù)據(jù)中的最優(yōu)周年項(xiàng)和錢德勒項(xiàng)數(shù)值，然后應(yīng)用所確定的具體數(shù)值來(lái)進(jìn)行最小二乘的擬合進(jìn)而預(yù)報(bào)未來(lái)極移數(shù)據(jù)，通過(guò)此種方法預(yù)報(bào)的極移精度要比以往國(guó)內(nèi)外學(xué)者應(yīng)用同種方法的具體周期數(shù)所得的預(yù)報(bào)精度都要高，因此，建議在進(jìn)行地球定向參數(shù)預(yù)報(bào)之前，應(yīng)首先確定所用數(shù)據(jù)的具體周期頻譜值，在極移預(yù)報(bào)中加以考慮，可進(jìn)一步提高預(yù)報(bào)的精度和可靠性。

[1]鄭大偉，虞南華．地球自轉(zhuǎn)及其地球物理現(xiàn)象的聯(lián)系：I日長(zhǎng)變化[J]．地球物理學(xué)進(jìn)展，1996，11(2)：81-101．

[2]AKYILMAZ O，KUTTERER H．Prediction of Earth Rotation Parameters by Fuzzy Inference Systems[J]．Journal of Geodesy，2004(78)：82-93．

[3]SCHUH H，ULRICH M，EGGER D，et al．Prediction of Earth Orientation Parameters by Artificial Neural Networks[J]．Journal of Geodesy，2002(76)：247-258．

[4]KOSEKW，KALARUSM，NIEDZIELSKIT．Forecasting of the Earth Orientation Parameters—Comparison of Different Algorithms[C]∥ Proceedings of the Journèes 2007，SystèMes De Référence Spatio-temporels“The Celestial Reference Frame for the Future”．Paris：[s．n．]，2007：155-158．

[5]KOSEKW．Future Improvements in EOP Prediction[C]∥Proceedings of the IAG 2009，“Geodesy for Planet Earth”．Argentina：[s．n．]，2009．

[6]KOSEK W．Causes of Prediction Errors of Pole Coordinates Data[C]∥Proceedings of the 6th Orlov’s Conference，“The Study of the Earth As a Planet by Methods of Geophysics，Geodesy and Astronomy”Ukraine：[s．n．]，2009：96-103．

[7]KALARUSM，SCHCH H，KOZEK W，et al．Achievements of Earth Orientation Parameters Prediction Comparison Campaign[J]．Journal of Geodesy，2010(84)：587-596．

[8]SCHUH H，NAGEL S，SEITZ T．Linear Drift and Periodic Variations Observed in Long Time Series of Polar Motion[J]．Journal of Geodesy，2001(74)：701-710．

[9]王琪潔．基于神經(jīng)網(wǎng)絡(luò)技術(shù)的地球自轉(zhuǎn)變化預(yù)報(bào)[D]．上海：中國(guó)科學(xué)院上海天文臺(tái)，2007．

[10]XU X Q，ZHOU Y H，LIAO X H．Short-term Earth O-rientation Parameters Predictions by Combination of the Least-squares，AR Model and Kalman Filter[J]．Journal of Geodynamics，2011(1107)：4．

[11]張昊，王琪潔，朱建軍，等．加權(quán)最小二乘法與AR組合模型在極移預(yù)測(cè)中的應(yīng)用研究[J]．天文學(xué)進(jìn)展，2011，29(3)：343-352．

[12]KOSEK W，MCCARTHY D D，LUZUM B J．El Nino Impact on Polar Motion Prediction Errors[J]．Stadia Geophysica et Geodaetica，2001(45)：347-361．

[13]DENNISDM，BRIAN JL．Prediction of Earth Orientation[J]．Journal of Geodesy，1991(65)：18-21．