基于ARIMA模型的地球自轉(zhuǎn)變化預(yù)報

李英冰

(武漢大學(xué)測繪學(xué)院,湖北 武漢 420079)

0 引 言

地球自轉(zhuǎn)運動表征著地球的整體運動狀況,以及地核、地幔、地殼、海洋和大氣等各圈層之間的相互作用過程[7]。地球自轉(zhuǎn)運動主要包括兩部分:一是地球自轉(zhuǎn)速率的變化,用周日變化(UT1-UTC)表征[1-7]。圖1是1992年1月1日至2008年7月30日UT1-UTC的時間序列圖,數(shù)據(jù)采樣率為1天,該數(shù)據(jù)從地球自轉(zhuǎn)服務(wù)(IERS)處下載。二是地球瞬時自轉(zhuǎn)軸相對于地球表面的運動,簡稱為極移(PM)[8-10],PM的主要激發(fā)因素是地球上物質(zhì)分布的變化[1,8]。圖2是PM 的時間序列圖,圖中虛線是x分量,實線為y分量。

圖1 UT1-UTC的時間序列

圖2 極移分量的時間序列

地球自轉(zhuǎn)變化的預(yù)報具有重要的科學(xué)意義和實用價值:

1)地球自轉(zhuǎn)參數(shù)是精密大地測量必不可少的基本數(shù)據(jù),如GMAIT、BERNESE、GIPSY 等精密解算軟件必須實時更新準(zhǔn)確的地球自轉(zhuǎn)變化數(shù)據(jù)。

2)用于月球、火星等星球探測的空間飛行器的正常運行需要準(zhǔn)確的地球自轉(zhuǎn)變化數(shù)據(jù)。在航天器導(dǎo)航和跟蹤中,地球自轉(zhuǎn)變化的不確定性是主要的誤差源。如果地球自轉(zhuǎn)參數(shù)不準(zhǔn)確,有可能導(dǎo)致衛(wèi)星飛行器不能準(zhǔn)確入軌。

近幾十年來,隨著甚長基線干涉、人衛(wèi)激光測距和全球定位系統(tǒng)等現(xiàn)代測地技術(shù)的飛速發(fā)展,極大地推動了地球自轉(zhuǎn)運動的研究。美國海軍天文臺的McCarthy和美國噴氣推進實驗室的Gross等專家長期專注于地球自轉(zhuǎn)變化的預(yù)報研究[3,5]。常用方法有最小二乘法、協(xié)方差法、神經(jīng)網(wǎng)絡(luò)法、自回歸移動平均法、小波分析、卡爾曼濾波等方法[2-5,9],預(yù)報準(zhǔn)確度已經(jīng)達到較高水平。

1 大氣與海洋變化對地球自轉(zhuǎn)變化的影響

全球流體變化會對地球自轉(zhuǎn)變化產(chǎn)生影響[1]。自轉(zhuǎn)運動與大氣和海洋活動的研究受到廣泛關(guān)注,1995年,IUGG/IAG專門設(shè)立了大氣和海洋與地球自轉(zhuǎn)相互作用的專題研究組。1998年,IERS組建了全球地球物理流體中心(GGFC),其下轄有大氣、海洋等7個分中心,各分中心定期公布最新的觀測數(shù)據(jù)和模型研究成果[1,10]。

對大氣角動量(AAM)和海洋角動量(OAM)進行積分運算,能夠分別計算出大氣和海洋變化對地球自轉(zhuǎn)變化的貢獻。其中AAM數(shù)據(jù)可以從美國國家環(huán)境預(yù)報中心(NCEP)獲取,OAM數(shù)據(jù)可從GGFC下載,數(shù)據(jù)采樣率都為6 h。

通過對AAM 和OAM的chi3分量積分,可得大氣和海洋變化對UT1-UTC的貢獻,積分結(jié)果表示為UTAAM+UTOAM。通過對AAM和OAM的軸分量積分,可得大氣和海洋變化對PM的貢獻,表示為PMAAM+PMOAM。積分公式為

圖3是UTAAM+UTOAM的時間序列,時間區(qū)間為2001年1月1日至2008年7月30日。圖4是PMAAM+PMOAM的時間序列,實線是x分量,虛線為y分量。

圖3 UTAAM+UTOAM/(ms)

圖4 PMAAM+PMOAM/(mas)

2 周日變化的預(yù)報

周日變化的預(yù)報算法如圖5所示,主要包括建立預(yù)測模型和外推預(yù)測兩個部分。預(yù)測模型建立的步驟包括減去周跳、潮汐、UT2-UT1、趨勢項和周期項等。外推預(yù)測主要包括基于周日殘差序列應(yīng)用ARIMA進行外推預(yù)報,加上周期項、趨勢項、UT2-UT1、潮汐、周跳等預(yù)測值。

建立預(yù)測模型主要步驟為

1)數(shù)據(jù)準(zhǔn)備,從www.iers.org下載Finals.all(IAU2000)文件,抽取 UT1-UTC數(shù)據(jù)列,如圖1所示。

2)UT1-UTC是不連續(xù)的,存在跳秒問題。因而需要移除跳秒,產(chǎn)生連續(xù)的UT1-TAI數(shù)據(jù)列。

式中 :詳見文獻[6],進行潮汐改正后的結(jié)果記為UT1R-TAI。

(4)減去UT2-UT1,其中UT2-UT1的計算公式為

計算結(jié)果記為UT2R-TAI。以上數(shù)據(jù)的時間序列如圖6所示。

(5)用最小二乘法計算趨勢項和周期項,公式為

式中 :A0、A1、A2是趨勢項的擬合系數(shù),B0、B1、B2為周期項的擬合系數(shù),φ1、φ2、φ3 為初始相位,ω1、ω2、ω3為周期項的頻率,分別對應(yīng)于半年、一年和Chandler周期項。移除趨勢項和周期項后的數(shù)據(jù)序列記為 ΔΔ(UT2R-TAI)。

(6)利用AAM和OAM數(shù)據(jù)計算大氣和海洋 變化對周日變化的貢獻。

圖5 UT1-UTC的預(yù)報算法

周日變化的外推預(yù)報過程包括

1)基于差分自回歸移動平均模型(ARIMA)對序列(UT2R-TAI)進行外推預(yù)報。公式為

式中:X(t)=a1Xt-1+…+apXt-p+et+b1et-1+…+bqet-q,a1,……,ap是P階自回歸系統(tǒng),b1,……,bq為q階移動平均系數(shù),d為差分階數(shù)。當(dāng)UTAAM和UTOAM參與計算時,需要采用多元回歸模型(MAR)。

2)根據(jù)周期項的擬合系數(shù),外推周期項的預(yù)報序列,并與 ΔΔ(UT2R-TAI)的預(yù)報值相加,得到殘差和周期項Δ(UT2R-TAI)的預(yù)報值。

3)根據(jù)多項式的擬合系數(shù),外推趨勢項,將它與前面的預(yù)報結(jié)果相加,得到UT2R-TAI的預(yù)測值。

4)由公式(3)計算UT2-UT1的外推估值,并與前面的預(yù)測值求和,得到UT1R-TAI的預(yù)測值。

5)加上由公式(2)得到的潮汐預(yù)測估值,得到UT1-TAI的預(yù)測值。

6)加上跳秒,得到UT1-UTC的預(yù)測值。

利用上述算法進行實際計算,預(yù)測未來365天的UT1-UTC的數(shù)值。共進行實際預(yù)報720組,將預(yù)測值與實測值比較,統(tǒng)計預(yù)測誤差絕對值的平均值(MAE),計算公式為

式中:εi,n是預(yù)報誤差。圖7中的虛線是周日變化的MAE統(tǒng)計曲線,實線是對IERS官方預(yù)報序列的MAE統(tǒng)計。我們的預(yù)報效果總體優(yōu)于IERS成果,特別是近期預(yù)報。

3 極移變化的預(yù)報

極移變化的預(yù)報算法如圖8所示,包括建立預(yù)測模型和外推預(yù)測兩個部分。

圖8 極移的預(yù)報算法

在建立預(yù)測建模時,用最小二乘法計算出趨勢項和周期項的擬合系數(shù),計算公式為

上式中的符號含義與公式(4)一樣。對殘差序列的預(yù)報采用季節(jié)性自回歸移動平均(SARMA)模型進行預(yù)報,公式為

式中S為季節(jié)項,其他參數(shù)含義與公式(5)相似。

圖9 極移變化X分量的MAE統(tǒng)計圖形

利用上述算法進行實際計算,預(yù)報未來365天的極移變化。進行實際預(yù)報720組,將預(yù)測值與實測值比較,并統(tǒng)計出MAE序列。圖9和圖10分別給出了極移X和Y分量的MAE統(tǒng)計曲線及其與IERS成果的對比。統(tǒng)計結(jié)果表明:一周內(nèi)的超短期預(yù)報優(yōu)于IERS成果,但整體結(jié)果要比IERS的成果略差。

圖10 極移變化X分量的MAE統(tǒng)計圖形

4 結(jié) 論

采用最小二乘法和ARIMA模型對周日變化進行預(yù)測,總體的預(yù)報效果要優(yōu)于IERS發(fā)布的產(chǎn)品。采用最小二乘法和SARMA對極移變化進行預(yù)報,總體預(yù)報效果沒有IERS的預(yù)報效果好。當(dāng)采用AAM和OAM數(shù)據(jù)時,對于周日預(yù)報有輕微改善,但對極移的預(yù)報沒有得到改善。

致謝:感謝美國俄亥俄州立大學(xué)的郭俊義和C.K.Shum教授提供的訪美經(jīng)費支持。

[1]Chao B F,Dehant V,Gross R S,and et al.Space geodesy monitors mass transports in global geophysical fluids[J].Eos Trans.,2000,81(22):247,249-250.

[2]Dickman S R.Determination of oceanic dynamic barometer corrections to atmospheric excitation of Earth rotation[J].J.Geophys.Res.,1998,103(B7):15127-15144.

[3]Gross R S,Eubanks T M,Steppe J A,and et al.AKalman-filter-based approach to combining independent Earth-orientation series[J].J.Geodesy,1998,72(4):215-235.

[4]Li Y B,Guo J Y,Shum C K,Johnson T.Empirical predictions of UT1 and polar motion using least squares fit and ARIMA[C]//AGU full meeting,San Francisco,2008.

[5]McCarthy,Luzum.Prediction of Earth Oriention[J].Bull Geod.,1991,65(1):18-22.

[6]McCarthy,D D,Petit G,IERS Conventions[R].Paris:Observatoire de Paris,2003.

[7]鄭大偉,虞南華.地球自轉(zhuǎn)及其和地球物理現(xiàn)象的聯(lián)系:I日長變化[J].地球物理學(xué)進展,1996,11(2):81-101.

[8]虞南華,鄭大偉.地球自轉(zhuǎn)及其和地球物理現(xiàn)象的聯(lián)系:II地極運動[J].地球物理學(xué)進展,1996,11(3):71-81.

[9]王琪潔.基于神經(jīng)網(wǎng)絡(luò)技術(shù)的地球自轉(zhuǎn)變化預(yù)報[D].中國科學(xué)院研究生院博士學(xué)位論文,2007.

[10]周永宏,鄭大偉,虞南華,廖新浩.地球自轉(zhuǎn)運動與大氣、海洋活動[J].科學(xué)通報,2000,45(24):2258-2597.