火星極移時變分析與預報

魏二虎 王恒毅 白小宇 李巖林 劉經南

1 武漢大學測繪學院,武漢市珞喻路129號,430079 2 中國地質大學(武漢)地理與信息工程學院,武漢市魯磨路388號,430079 3 武漢大學衛星導航定位技術研究中心,武漢市珞喻路129號,430079

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

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

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

1 星歷數據及處理

1.1 Horizons System簡介

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

1.2 星歷數據預處理

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

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

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

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

2 火星極移時變分析

2.1 離散傅里葉變換(DFT)

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

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

(1)

(2)

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

(3)

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

2.2 周期項結果與分析

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

圖2 NP.ang和NP.ds頻譜

通過圖2可以發現,存在5個較為明顯的周期項,其對應的頻率、周期及幅值如表1所示。

表1 NP.ang和 Np.ds周期項統計

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

θ=

(4)

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

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

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

(5)

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

3 火星極移預報

3.1 短期預報模型

3.1.1 ARMA模型

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

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

(6)

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

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

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

(7)

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

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

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

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

(8)

即當期值Yt不僅與滯后項有關,還與隨機擾動有關。

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

3.1.2 ARIMA模型

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

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

3.1.3 短期預報

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

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

圖3 自相關與偏自相關

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

3.2 中長期預報模型

3.2.1 最小二乘外推模型

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

(9)

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

3.2.2 中長期預報

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

3.3 預報結果與分析

3.3.1 短期預報結果與分析

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

圖4 ARIMA模型與ARMA模型預報誤差

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

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

圖5 未來20 d的ARIMA模型預報誤差

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

3.3.2 中長期預報結果與分析

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

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

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

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

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

(10)

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

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

使用該組周期項進行擬合與預報,得到未來5 a的NP.ds中長期預報值與實際值對比圖(圖(7))。可以看出,2條曲線幾乎完全重合,只有在峰值附近(圖8)才能發現預報結果仍有一定的誤差。

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

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

4 結 語

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