魏二虎 王恒毅 白小宇 李巖林 劉經(jīng)南

1 武漢大學(xué)測繪學(xué)院,武漢市珞喻路129號,430079 2 中國地質(zhì)大學(xué)(武漢)地理與信息工程學(xué)院,武漢市魯磨路388號,430079 3 武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心,武漢市珞喻路129號,430079

火星是類地行星,具有與地球相近的“黃赤交角”和自轉(zhuǎn)周期,也存在歲差、章動與極移等現(xiàn)象[1]。和地球極移一樣,火星極移是指火星自轉(zhuǎn)軸相對于火星參考系極軸的運動,與火星的內(nèi)核運動及表面季節(jié)性變化有關(guān)。通過對火星極移時間序列各周期項進行研究,可以更加深入地了解火星液核、地幔等結(jié)構(gòu)及火星大氣、潮汐運動的特征與規(guī)律。而對火星極移進行預(yù)報,則可以為火星探測器的精密定軌及著陸器的定位提供所需的定向參數(shù),對太陽系起源與演化研究及人類航天事業(yè)發(fā)展都有著十分重要的意義[2]。

目前對火星定向參數(shù)的研究仍是將地球自轉(zhuǎn)運動的各項結(jié)果拓展到火星上[3],而對于火星極的空間運動,既可以用理論進行模擬,也需要用由火星探測獲得的數(shù)據(jù)資料進行約束[4]。將火星極運動的理論值與由火星探測得到的觀測值進行比較,是檢驗當(dāng)下火星模型的重要手段,也是為改進與火星定向參數(shù)有關(guān)的各項理論提供依據(jù)的有效途徑。因此,開展火星極移參數(shù)預(yù)報模型的相關(guān)研究至關(guān)重要。

本文首先利用NASA提供的NP.ang和NP.ds火星極移序列進行時變分析與預(yù)報,并分別利用自回歸滑動平均(autoregressive moving average model,ARMA)模型和ARIMA模型對NP.ds進行短期預(yù)報實驗;再利用最小二乘外推模型對NP.ds進行中長期預(yù)報實驗,在對離散傅里葉變換結(jié)果進行分析后,改進用于擬合的周期項數(shù)值,以實現(xiàn)更高精度的中長期預(yù)報。

1 星歷數(shù)據(jù)及處理

1.1 Horizons System簡介

本文采用的星歷數(shù)據(jù)來自于NASA噴氣推進實驗室太陽系動力學(xué)小組提供的Horizons System。Horizons System是一個在線太陽系數(shù)據(jù)和星歷計算服務(wù)系統(tǒng),可以向用戶提供對太陽系關(guān)鍵數(shù)據(jù)的訪問途徑,還可以針對太陽系內(nèi)的對象靈活生成高精度星歷表,這些對象包括1 204 271顆小行星、3 794顆彗星、211顆行星的衛(wèi)星(包括地球的衛(wèi)星及矮行星冥王星)、太陽、太陽系內(nèi)全部8顆行星、經(jīng)過挑選的航天器、拉格朗日點L1、L2及太陽系質(zhì)心等[5]。

1.2 星歷數(shù)據(jù)預(yù)處理

由于本文需要對得到的火星極移數(shù)據(jù)時間序列進行離散傅里葉變換,以提取其中的周期項,因此需要盡可能多地獲取火星極移數(shù)據(jù),以提高離散傅里葉變換結(jié)果的精度。本文獲取公元1600-01-02～2500-01-01長達900 a的星歷數(shù)據(jù),數(shù)據(jù)步長設(shè)置為1 d。

由于本文的研究對象為火星極移的時間序列,歷表設(shè)置選定的星歷數(shù)據(jù)類型均為與火星北極位置變化有關(guān)的數(shù)據(jù)類型,包括NP.ang和NP.ds。其中,NP.ang指在數(shù)據(jù)輸出時刻的火星北極位置角(即相對于真北天極方向的角度,單位(°)),NP.ds指在數(shù)據(jù)輸出時刻火星北極與火星視中心的角距(單位(″))。

通過觀察原始星歷數(shù)據(jù)時間序列發(fā)現(xiàn),NP.ang和NP.ds存在跳變現(xiàn)象。在查閱Horizons System關(guān)于輸出數(shù)據(jù)類型的介紹后了解到,NP.ang的跳變是因為該時刻北極位置角小于0°,由于系統(tǒng)沒有設(shè)定負角度,位置角會被自動加上360°。而NP.ds的跳變是因為當(dāng)火星北極位于遠離觀測者的那部分半球時,其與視中心的角距將變?yōu)樨撝礫6]。為方便后續(xù)時變分析與預(yù)報計算,對這2組數(shù)據(jù)中發(fā)生跳變的部分進行相應(yīng)處理,結(jié)果如圖1所示。

圖1 1600～1610年NP.ang 和 Np.ds時間序列

2 火星極移時變分析

2.1 離散傅里葉變換(DFT)

離散傅里葉變換利用測量到的原始信號,以累加的方式計算該信號中不同正弦波信號的頻率、振幅和相位,最終得到的頻譜圖中振幅大的頻率代表周期性強。

離散傅里葉變換的定義如式(1)所示,其中k為分析的特定采樣次數(shù)。雖然理論上離散周期信號的頻譜是無窮多的,但由于傅里葉級數(shù)具有周期性,一般只取主值區(qū)間0≤k

(1)

(2)

N個采樣點經(jīng)離散傅里葉變換后的結(jié)果為N個復(fù)數(shù),每個復(fù)數(shù)對應(yīng)1個頻率(第n≤N/2個點對應(yīng)的頻率為(n-1)/Nfs),該復(fù)數(shù)的模值表示該頻率的振幅特征。為提取周期項,需要求出相關(guān)參數(shù)頻率和幅值,頻率計算公式見式(3),式中,fs為采樣頻率:

(3)

快速傅里葉變換(FFT)是指計算機高效、快速計算離散傅里葉變換(DFT)的方法。離散傅里葉變換實現(xiàn)時間的復(fù)雜度為O(n2),而快速傅里葉變換時間復(fù)雜度縮小到了O(nlogn)。

2.2 周期項結(jié)果與分析

通過使用C#語言編寫FFT程序,分別對NP.ang與NP.ds的時間序列進行處理,得到2組數(shù)據(jù)的頻譜(圖2)。

圖2 NP.ang和NP.ds頻譜

通過圖2可以發(fā)現(xiàn),存在5個較為明顯的周期項,其對應(yīng)的頻率、周期及幅值如表1所示。

表1 NP.ang和 Np.ds周期項統(tǒng)計

由表1可見,NP.ang與NP.ds的周期項長度基本一致,且與火星年具有很強的相關(guān)性:第1組周期項的周期接近1個火星年,第2組周期項接近0.5個火星年,第3組接近1/3個火星年,第4組接近1/4個火星年,第5組則接近1/5個火星年。假設(shè)在沒有火星極移的情況下,即在火星北極嚴格指向平北天極的情況下,火星北極點在ICRS/J2000框架下的太陽系質(zhì)心慣性坐標為(X,Y,Z),火星圓盤視中心坐標為(XC,YC,ZC),火星質(zhì)心到太陽系質(zhì)心的距離減去火星赤道半徑為R,則NP.ds(記為θ)滿足公式:

θ=

(4)

由于火星公轉(zhuǎn)軌道是橢圓軌道,太陽系質(zhì)心位于橢圓軌道的一個焦點上,這使得R存在因公轉(zhuǎn)造成的周期性變化(周期為1個火星年,約687 d)。在不考慮極移的情況下,影響NP.ds周期性變化的因素只有R的周期性變化。也就是說,如果不存在極移,NP.ds的時間序列應(yīng)當(dāng)有且僅有長度為1個火星年的周期項。同樣,在不考慮火星極移的情況下,NP.ang也只與R的周期性變化相關(guān),其時間序列也應(yīng)當(dāng)有且僅有長度為1個火星年的周期項。然而,排除掉NP.ang與NP.ds時間序列中接近1個火星年的周期項后,仍存在4個較為明顯的周期項(表2),由此可以認定這些周期項應(yīng)當(dāng)與火星極移有關(guān)。

表2 可能包含極移信息的周期項

通過查閱與火星定向參數(shù),尤其是火星極移有關(guān)的文獻資料了解到,火星極移主要由兩部分組成:其一是季節(jié)性變化項,主要和火星大氣與其表面質(zhì)量交換等原因?qū)е碌淖赞D(zhuǎn)角動量變化有關(guān);其二是與地球極移相似的錢德勒(Chandler)擺動,這是非剛體火星的自由擺動,與火星的自由核章動(周期約為230 d,接近1/3個火星年)有關(guān)[7]。因此,火星極移參數(shù)可表示為:

(5)

式中,前4項的fm=m/687代表由火星季節(jié)性變化導(dǎo)致的極移變化頻率,而最后1項f5表示火星錢德勒擺動的頻率[8]。據(jù)此可以確定,表2中前3項均為對應(yīng)季節(jié)性變化項。

3 火星極移預(yù)報

3.1 短期預(yù)報模型

3.1.1 ARMA模型

自回歸滑動平均(ARMA)模型是以自回歸(AR)模型與滑動平均(MA)模型為基礎(chǔ)混合得到的模型,其中AR模型是變量對變量自身的滯后期項進行的回歸。具體來說,p階自回歸模型(AR(p))的當(dāng)期值Yt滿足方程:

Yt=φ1Yt-1+φ2Yt-2+…+φpYt-p+et

(6)

即Yt為距離t期最近的p階滯后項與隨機擾動et的線性組合,其中et包括了被擬合的平穩(wěn)序列在t期無法用過去值來解釋的所有新信息。

而MA模型是將時間序列寫成一系列不相關(guān)的隨機變量線性組合。具體地講,q階自回歸模型(MA(q))的當(dāng)期值Yt滿足方程:

Yt=μ+et-θ1et-1-θ2et-2-…-θqet-q

(7)

即Yt與以前各期的序列值無關(guān),而是基于前q期隨機擾動et-1、et-2、…、et-q的線性回歸模型。式中,μ為被擬合的平穩(wěn)序列均值。

ARMA模型綜合了二者的特點。階數(shù)分別為p和q的自回歸滑動平均模型(ARMA(p,q))的當(dāng)期值Yt滿足方程:

Yt=φ1Yt-1+φ2Yt-2+…+φpYt-p+

et-θ1et-1-θ2et-2-…-θqet-q

(8)

即當(dāng)期值Yt不僅與滯后項有關(guān),還與隨機擾動有關(guān)。

ARMA模型階數(shù)的確定主要通過計算平穩(wěn)序列的自相關(guān)系數(shù)(ACF)與偏自相關(guān)系數(shù)(PACF),并生成對應(yīng)的序列自相關(guān)圖與偏自相關(guān)圖,根據(jù)兩幅圖中系數(shù)開始衰減的滯后期數(shù)來確定ARMA模型的階數(shù)。

3.1.2 ARIMA模型

差分整合移動平均自回歸(ARIMA)模型的基本思想與ARMA模型相似:首先假定預(yù)測對象隨時間的推移形成的數(shù)據(jù)與其過去的數(shù)據(jù)有關(guān),這樣就可以建立自回歸模型及滑動平均模型,并將模型外推來進行預(yù)報。對于非平穩(wěn)時間序列,ARIMA模型會采用差分方法使其顯示出平穩(wěn)的性質(zhì),再用于擬合與預(yù)報。

ARIMA模型具有3個參數(shù),即在ARMA模型的基礎(chǔ)上增加了代表時間序列從非平穩(wěn)序列變?yōu)槠椒€(wěn)序列所要進行的差分階數(shù)d[9],對應(yīng)的ARIMA模型可記為ARIMA(p,d,q),其中AR項p與MA項q的確定方法與ARMA模型類似,而差分階數(shù)d則需要通過ADF檢驗與KPSS檢驗來確定。

3.1.3 短期預(yù)報

考慮到本文使用的星歷數(shù)據(jù)是JPL通過對火星衛(wèi)星的數(shù)值進行積分,并使用地球觀測數(shù)據(jù)和火星探測器影像觀測數(shù)據(jù)進行擬合得來的,過晚或過早時期的星歷數(shù)據(jù)不宜作為預(yù)報的擬合值,應(yīng)采用與開展火星定向參數(shù)觀測任務(wù)相近時期的星歷數(shù)據(jù),故本文短期預(yù)報部分采用2000-01-02～2010-01-01星歷數(shù)據(jù),共計3 653 d。

因為在時變分析時已經(jīng)發(fā)現(xiàn),NP.ang具有的周期項NP.ds同樣具備,而且各小周期項幅值相對于火星周年項更大,加之NP.ds具有更高的精度,所以僅以NP.ds為例進行ARIMA模型的擬合。經(jīng)過3次差分后發(fā)現(xiàn),所得序列的ADF檢驗與KPSS檢驗結(jié)果已滿足要求,因此,差分階數(shù)d定為3。對二階差分序列繪制自相關(guān)圖與偏自相關(guān)圖,根據(jù)圖3可以判斷,p為21、q為4,因此選用的擬合模型應(yīng)為ARIMA(21,3,4)。

圖3 自相關(guān)與偏自相關(guān)

此外,為了對比ARIMA模型與ARMA模型在非平穩(wěn)時間序列擬合與預(yù)報上的差距,在進行ARIMA模型擬合與預(yù)報的同時也進行ARMA模型的擬合與預(yù)報。預(yù)報結(jié)束后,將二者的預(yù)報值分別與實際值作差,并繪制預(yù)報誤差曲線,進而直觀地對比二者在非平穩(wěn)時間序列擬合與預(yù)報上的差距。

3.2 中長期預(yù)報模型

3.2.1 最小二乘外推模型

最小二乘外推模型在研究地球極移長期預(yù)報問題的線性模型預(yù)報法中占有重要地位[10],本文選用的NP.ds時間序列的周期性要好于由IERS提供的地球極移參數(shù)結(jié)果,因此使用最小二乘外推模型對NP.ds時間序列進行擬合并進行中長期預(yù)報是可行的。在§2.2中已經(jīng)確定了5個較為明顯的周期項,因此本文最小二乘外推模型為:

(9)

式中,a0為常數(shù)項,a1為線性趨勢項系數(shù),a2～a11為周期項系數(shù)。

3.2.2 中長期預(yù)報

中長期預(yù)報采用1991-01-02～2010-01-01星歷數(shù)據(jù),共計6 940 d。在確定用于擬合最小二乘外推模型的星歷數(shù)據(jù)后,按照式(9)進行中長期預(yù)報。

3.3 預(yù)報結(jié)果與分析

3.3.1 短期預(yù)報結(jié)果與分析

對所選時間序列分別進行ARMA模型及ARIMA模型擬合,并對未來50 d的極移進行預(yù)報,ARIMA模型的預(yù)報值與實際值的偏差及ARMA模型的預(yù)報值與實際值的偏差如圖4所示。

圖4 ARIMA模型與ARMA模型預(yù)報誤差

由圖4可明顯看出,前10 d兩種模型預(yù)報誤差的差異不太顯著,但在第10 d后,ARIMA模型的預(yù)報效果顯著優(yōu)于ARMA模型,驗證了ARIMA模型在處理非平穩(wěn)時間序列時的優(yōu)越性。

當(dāng)把預(yù)報時間縮短至未來20 d時,ARIMA模型的預(yù)報值與實際值的偏差如圖5所示。由圖可見,使用ARIMA模型對未來20 d的NP.ds數(shù)據(jù)進行預(yù)報,其誤差基本在1.0 mas以內(nèi),小于星歷數(shù)據(jù)本身的誤差,這說明使用ARIMA模型對NP.ds進行20 d內(nèi)的短期預(yù)報是可行的。

圖5 未來20 d的ARIMA模型預(yù)報誤差

此外,通過觀察圖4后半段的預(yù)報誤差趨勢發(fā)現(xiàn),隨著時間的推移,ARIMA模型的預(yù)報誤差逐漸增大,尤其是第23 d后預(yù)報誤差迅速增大,說明ARIMA模型不適合對NP.ds進行中長期預(yù)報。

3.3.2 中長期預(yù)報結(jié)果與分析

對所選時間序列進行最小二乘外推模型(選用的周期項為§2.2中得到的5個周期項)擬合,并對未來5 a的極移進行預(yù)報,NP.ds中長期預(yù)報值與實際值的對比如圖6所示。由圖6(a)可見,預(yù)報值整體趨勢與實際值基本一致,除了2個峰值差距較大外,2條曲線基本重合。事實上,即使是未來5 a中峰值的最大預(yù)報誤差也在20 mas以內(nèi),這說明使用最小二乘外推模型進行NP.ds數(shù)據(jù)的中長期預(yù)報是可行的。

圖6 未來5 a中長期預(yù)報值與實際值對比及預(yù)報誤差

圖7 未來5 a中長期預(yù)報值與實際值對比

此外,通過進一步觀察預(yù)報值與實際值差值的時間序列(圖6(b))可以發(fā)現(xiàn),預(yù)報誤差具有明顯的周期性特征。

初步推測,這是由于最小二乘外推模型中采用的通過離散傅里葉變換獲取的周期項結(jié)果與時間序列中蘊含的真正周期有一定差距,導(dǎo)致預(yù)報結(jié)果出現(xiàn)周期性誤差。事實上,由離散傅里葉變換提取的周期項均應(yīng)滿足公式:

(10)

式中,N為樣本總數(shù);k為頻率,是0～(N-1)之間的整數(shù)?？梢钥闯?在頻率較小時,周期項之間周期長度的變化跨度較大,由離散傅里葉變換提取的長周期項結(jié)果與真實值的差距理論上要大于短周期項。換言之,就是短周期項的結(jié)果更接近時間序列中蘊含的真實周期。考慮到前文得到的各周期項與火星年具有很強的相關(guān)性,以其中最短的周期項(137.39 d)作為參考,得到一組新的用于擬合最小二乘外推模型的周期項(表3)。

表3 用于擬合最小二乘外推模型的周期項

使用該組周期項進行擬合與預(yù)報,得到未來5 a的NP.ds中長期預(yù)報值與實際值對比圖(圖(7))?？梢钥闯?2條曲線幾乎完全重合,只有在峰值附近(圖8)才能發(fā)現(xiàn)預(yù)報結(jié)果仍有一定的誤差。

圖8 中長期預(yù)報值與實際值在峰值附近的對比

不過,這一結(jié)果仍驗證了通過使用準確度更高的周期項進行最小二乘外推模型擬合,可顯著提高對NP.ds中長期預(yù)報的精度,也進一步證明使用最小二乘外推模型對NP.ds進行中長期預(yù)報具有可行性。

4 結(jié) 語

本文利用NASA提供的NP.ang和NP.ds等包含火星極移信息的星歷數(shù)據(jù)時間序列進行時變分析與預(yù)報,結(jié)果表明,NP.ang與NP.ds兩組數(shù)據(jù)的周期項長度基本一致,且與火星年具有很強的相關(guān)性;短期預(yù)報的實驗結(jié)果表明,使用ARIMA模型對NP.ds進行20 d內(nèi)的短期預(yù)報是可行的;中長期預(yù)報實驗結(jié)果表明,使用最小二乘外推模型進行NP.ds的中長期預(yù)報是可行的,且使用準確度更高的周期項進行最小二乘外推模型擬合,可顯著提高對NP.ds中長期預(yù)報的精度。本文結(jié)果可為火星極移的預(yù)報研究提供參考。