高梧桐，謝 攀，鄢建國

（1.武漢大學測繪學院，武漢，430070；2.上海衛(wèi)星工程研究所，上海，201100；3.武漢大學測繪遙感信息工程國家重點實驗室，武漢 430070）

引 言

“火衛(wèi)1”是火星兩顆自然衛(wèi)星中質量較大的一顆，自1877年被發(fā)現(xiàn)以來，人們對“火衛(wèi)1”進行了許多研究建模與觀測?，F(xiàn)有資料表明[1]，“火衛(wèi)1”質量為1.0659× 1016kg，平均半徑為11 km。軌道周期約為0.318 9 個地球日，平均軌道速度達2.138 km/s。近年來火星及其衛(wèi)星的探測成為各國航天任務的熱點。2003年，歐洲航天局（European Space Agency，ESA）發(fā)射的“火星快車”（Mars EXpress，MEX）進行了飛越“火衛(wèi)1”探測并收集了大量有價值的數(shù)據(jù)[2]。我國首次獨立自主的火星探測任務將于2020年進行，任務期間有可能開展“火衛(wèi)1”飛掠測量。日本宇宙航空研究開發(fā)機構（Japan Aerospace Exploration Agency，JAXA）也已經確定火衛(wèi)探測任務及“火衛(wèi)1”采樣返回計劃MMX（Martian Moons eXploration）[3]。更長遠來說，美國國家航空航天局（National Aeronautics and Space Agency，NASA）計劃在21世紀30年代將人類送往火星[4]。探究“火衛(wèi)1”的軌道，對于未來的火星探測有著重要意義。

由于“火衛(wèi)1”軌道的重要性，國際上多個機構開展了這方面的研究。目前“火衛(wèi)1”歷表可以從由噴氣動力實驗室（Jet Propulsion Laboratory，JPL）的火星衛(wèi)星星歷MAR097、俄羅斯的EPM2017 和法國的INPOP17a 獲取。MAR097 星歷是JPL 發(fā)表的火星衛(wèi)星星歷[5]，通過對火星衛(wèi)星的數(shù)值積分，使用地球觀測數(shù)據(jù)和探測器影像觀測數(shù)據(jù)進行擬合。該星歷的動力學模型中主要為火星重力場模型MGS95J[6]（包括8 階帶協(xié)項與5 階次田協(xié)項），“火衛(wèi)1”和“火衛(wèi)2”的互相影響，太陽、地月系統(tǒng)、木星和土星的攝動，“火衛(wèi)1”的天平動，“火衛(wèi)1”引起的火星潮汐。但由于計算軟件限制，并未使用重力場模型MGS95J對應的火星自轉模型，而是使用了對應的ICRF（International Terrestrial Reference Frame）參考框架，計算過程中通過星歷DE421 獲取太陽、行星及其衛(wèi)星的位置和常數(shù)。

EPM（Ephemeris of Planets and the Moon）是俄羅斯科學院應用天文學研究所（Institute of Applied Astronomy of the Russian Academy of Sciences，IAARAS）發(fā)表的一套高精度星歷數(shù)據(jù)集[7]，目前EPM最新版本為EPM2017及其加長版本EPM2017H，包括太陽及其行星（包括冥王星）、各行星衛(wèi)星、小行星等星歷數(shù)據(jù)，涵蓋時段分別為1787-2214 和10107BC-AD3036。 EPM2017 （H）參考系選取ICRF2，定向精度高于0.2 mas，時間系統(tǒng)選取質心動力學時（Barycentric Dynamical Time，TDB），在動力學模型方面，引入了太陽扁率、301 顆大型小行星、30 顆大型外海王星天體，小行星帶、柯伊伯帶等攝動影響。行星星歷部分，共解算了約270 個參數(shù)。通過IAARAS 官網的星歷計算服務[8]，選擇EPM2017星歷，可以計算“火衛(wèi)1”的星歷。

INPOP17a 是由法國天體力學與星歷計算研究所（Institute of Mechanics Celestial and Calculation of Ephemeris，IMCCE）發(fā)表的行星及衛(wèi)星星歷[9]，它提供各大天體在ICRF（International Terrestrial Reference Frame）坐標框架下軌道數(shù)據(jù)和TT-TDB 時間轉換參數(shù)。INPOP17a 在火星星歷部分擬合了當時最新的MEX、Mars Odyssey 探測器的觀測值信息，并在動力模型上引入了168顆小行星的攝動影響。與EPM一樣，通過IMCCE 官網計算服務[10]，以INPOP17a為行星星歷，可以獲取“火衛(wèi)1”的星歷數(shù)據(jù)。

本文結合最新的“火衛(wèi)1”攝動力模型，通過高精度數(shù)值積分方法計算“火衛(wèi)1”的軌道，并將計算結果與多個歷表進行驗證，確保模型的精度和可靠性。通過這一工作，可以分析不同的攝動力源對“火衛(wèi)1”軌道的影響，將對后續(xù)“火衛(wèi)1”探測任務的開展提供參考。

1 積分器及動力學模型

本文主要通過對“火衛(wèi)1”運動方程進行數(shù)值積分來實現(xiàn)對“火衛(wèi)1”的軌道預報。積分過程中初始狀態(tài)向量從現(xiàn)有星歷中獲取，方法采用精度較高的Rouge-Kutta-7(8)積分器[11]進行起步，然后使用計算性能較好的Adams-Bashforth-Moulton 積分器[12]進行后續(xù)積分。

經對“火衛(wèi)1”的受力分析可知，中心天體火星的引力是“火衛(wèi)1”受力的主要來源，計算過程中首先考慮了最基本的質點模型，但是火星內部質量分布并不均勻，形狀也不規(guī)則，考慮非球形攝動，將擾動位函數(shù)進行球諧展開，得到[13]

其中：Re為火星的半徑；為完全正規(guī)化的二階勒讓德函數(shù)；r、ψ、λ為“火衛(wèi)1”在火星固定坐標框架下的3個球坐標分量；為完全正規(guī)化后的重力場模型球諧系數(shù)。

對于中心天體以外的天體，其“火衛(wèi)1”的引力攝動用質點模型表示。引入了潮汐模型，考慮火星的二階固體潮Love 數(shù)，對于某一質量為Mj、半徑為ae的攝動天體，在火星上產生固體潮，對“火衛(wèi)1”的引力位函數(shù)為[14]

其中：r為“火衛(wèi)1”相對火星的距離；rj為“火衛(wèi)1”相對攝動天體的距離。

另外，若視“火衛(wèi)1”運行軌道為圓軌道，某一質量為Mj的天體，受廣義相對論效應的影響，對衛(wèi)星的引力加速度產生相應的修正[15]為

相對論效應修正量等于該天體對“火衛(wèi)1”的引力加速度與3v2/c2的乘積，對于一般的衛(wèi)星，3v2/c2約為3× 10-10，考慮足夠長時間時，相對論效應引起的誤差不容忽視。

該算法參考時間系統(tǒng)為以太陽系質心為中心的TDB 時間系統(tǒng)，起始時刻為J2000.0（2000－01－01 12:00:00）。參考框架選取火星慣性系J2000，計算火星重力場過程中涉及到的火星固定坐標系，由對應的重力場模型中給出的火星自轉模型求得。攝動天體的位置通過行星星歷DE430[16]獲得?；鹦呛汀盎鹦l(wèi)1”的質量和半徑從對應重力場模型獲取。此外，其它天體的質量、半徑等常數(shù)，UTC 到TDB 時間轉換過程中的閏秒數(shù)據(jù)，以及用于對比的“火衛(wèi)1”的星歷也從對應的SPICE Kernel 文件讀取，這些數(shù)據(jù)均可在NASA的NAIF機構官網[17]獲得。

2 動力學模型影響分析

對“火衛(wèi)1”軌道計算時，所引入的動力模型越完善，計算效果越好，但同時計算時間成本和空間成本隨之增大。數(shù)值積分時，既要盡可能保證力模型的全面性，又要考慮計算性能和效率。本節(jié)在考慮某項動力學模型和不考慮該模型的情況下，在一段時間內計算各自的軌道位置，由3 個位置分量的差異dx、dy、dz，計算位置差異dr，分析其影響量級，進而考慮是否引入該模型。

2.1 重力場階數(shù)

隨著采用重力場模型階數(shù)的增大，軌道預報結果精度將會提高。但是對于每一步積分求值，重力場模型計算的時間復雜度是O(n2)，其中n為重力場階數(shù)。應合理選擇重力場模型的階數(shù)，在合理的精度下，控制計算的時間成本。計算過程中選取重力場模型MRO120D 及其對應的火星自轉模型[18]進行重力場的計算，最高階數(shù)計算到20 階，后續(xù)階次的增大對計算結果沒有影響。火星引力常數(shù)為42 828.375 815 756 1 km3/s2，參考半徑為3 396 km。為確保中心引力模型的一致性，算法同時適配了一階重力場模型，即將中心天體視為質點的重力模型，同時比較了引入重力場模型與質點模型對計算結果影響的量級。試驗過程中隨機選取1982－1992 和2000－2010 兩個時間段，在不同階數(shù)重力場模型下計算軌道結果，并與MAR097星歷比較，計算不同階數(shù)10年最大距離誤差dr，結果如圖1 所示。

圖1 不同階數(shù)MRO120D重力場模型引起的最大軌道誤差Fig.1 Max distance difference caused by MRO120D gravity field model in different degrees

兩個時間段的比較結果具有一致性，均說明1階重力場模型（質點引力模型）能夠對“火衛(wèi)1”產生數(shù)千千米的軌道偏差，因此實際計算過程中不能簡單將火星視為質點。在重力場模型6階以后距離誤差的量級便穩(wěn)定下來，10 階以后重力場計算結果的差異僅在幾米到十幾米?？紤]計算成本，重力場模型階數(shù)選擇6～10階即可。

2.2 N體攝動

本節(jié)比較了太陽及除火星以外各天體系統(tǒng)的攝動對“火衛(wèi)1”軌道的影響，比較時段為10年。首先在兩次積分過程中，分別引入和不引入太陽的引力攝動，比較太陽攝動對“火衛(wèi)1”位置的影響，結果如圖2所示。

可以看到，在10年這一長周期的積分時段內，太陽攝動引起的“火衛(wèi)1”軌道差異最高可達30 km，并隨時間呈線性增長。曲線放大后，可以發(fā)現(xiàn)存在振幅為數(shù)千米的高頻振蕩。因此可以得到結論，太陽攝動無論在10年的長周期，還是在較短的時間段內，都會對“火衛(wèi)1”的軌道產生不可忽略的影響。

圖2 太陽攝動引起的軌道誤差Fig.2 Distance difference with and without considering perturbation of Sun

按照相同的思路，對太陽系各行星系統(tǒng)（包括其衛(wèi)星）攝動對“火衛(wèi)1”軌道的影響進行分析，需要注意的是，各行星系統(tǒng)的位置矢量指向的是對應的行星系統(tǒng)質心，而非單個行星的質心，質量也應考慮行星與衛(wèi)星的質量總和。對各行星系統(tǒng)的引力攝動影響分析結果如圖3所示。結果表明，在10年的時間尺度計算“火衛(wèi)1”軌道，行星系統(tǒng)引力攝動引起的軌道差異量級僅在幾米到十幾米不等，大多數(shù)攝動天體的作用可以忽略。而太陽攝動引起的距離差異能夠占所有三體攝動的99%以上，因此，計算過程中必須要考慮太陽的攝動。對于其它天體攝動，在長時間段內引起的距離誤差量級都是可以忽略的，是否引入可以依照差異量級的大小和考慮攝動模型的完善性進行選擇。

圖3 太陽系及其他行星系統(tǒng)攝動引起的軌道誤差Fig.3 Distance difference of other planets in solar system

2.3 潮汐攝動

對于每一個天體引起的潮汐模型的分析，與N體攝動分析方法一致。由于太陽與“火衛(wèi)1”是火星固體潮汐攝動的主要來源[19]，因此在潮汐模型中不再引入其余天體。在考慮火星二階“Love 數(shù)”的情況下，太陽和“火衛(wèi)1”使得火星固體部分產生彈性形變，進而改變了火星的質量和重力場分布，最終對“火衛(wèi)1”軌道產生的影響如圖4 所示。由圖4 可知，太陽和“火衛(wèi)1”固體潮汐在10年的時間段對“火衛(wèi)1”位置的影響最大分別為300 m和75 m左右，要大于行星系統(tǒng)的攝動引力影響。從誤差大小來看，可以在保證計算效率的前提下引入太陽和“火衛(wèi)1”的潮汐影響。

2.4 相對論效應

一般來說，天體相對論效應產生的對衛(wèi)星的引力修正項的大小取決于天體的施瓦西半徑（2GM/c2）和衛(wèi)星的速度。所以，與N體攝動、潮汐攝動一樣，所有對“火衛(wèi)1”產生相對論效應的天體中，太陽占主導作用。太陽引起的相對論效應對“火衛(wèi)1”軌道的影響結果見圖5。在短時間內相對論效應引起的誤差量級可以忽略，但是對于10年的積分間隔，相對論效應將會對“火衛(wèi)1”軌道產生最大5 km 的影響，誤差隨時間大致線性增長。綜上所述，最終的軌道預報模型，采用的動力模型及參數(shù)如表1所示。

圖4 太陽和“火衛(wèi)1”潮汐引起的軌道誤差Fig.4 Distance difference caused by tidal effect of Sun and Phobos

圖5 太陽相對論效應引起的軌道誤差Fig.5 Distance difference caused by relativity effect of Sun

表1 “火衛(wèi)1”軌道預報采用的動力模型及參數(shù)Table 1 Force models and parameters in orbit prediction of Phobos

3 星歷對比

本節(jié)將計算出的“火衛(wèi)1”軌道與3 種星歷進行差異比較，分別選取JPL 的MAR097 火星衛(wèi)星星歷，以及由俄羅斯EMP2017 和法國INPOP17a 求得的“火衛(wèi)1”星歷。比較過程中，首先由從各自星歷中獲取初始狀態(tài)向量，分別計算軌道數(shù)據(jù)，并與對應星歷相比較，計算距離差。其次，將狀態(tài)向量轉換為軌道根數(shù)（半長軸SA、離心率E、軌道傾角I、升交點赤經RAAN、近地點角距AOP 和平近點角MA），再分別做差比較。利用上述的動力模型及積分器，設定初始時間為1982－03－26 UTC，從3 種星歷中獲取的初始狀態(tài)向量如表2所示。

表2 “火衛(wèi)1”1982-03-26 UTC時刻下在3種星歷中的狀態(tài)向量Table 2 State Vectors of Phobos in three different ephemeris at 1982-03-26 UTC

與MAR097 星歷比較時，首先取一年的積分時長，比較短期內的距離誤差與軌道根數(shù)誤差，評價該算法模型的精度（結果見圖6 和圖7）；再取10年的積分時長，驗證模型的穩(wěn)定性（結果見圖8 和圖9）。按照同樣的方法，分別比較其軌道距離差異和軌道根數(shù)差異，對“火衛(wèi)1”軌道預報模型進行進一步驗證（結果分別見圖10～11和圖12～13）。

圖6 與圖7 的結果表明，該模型在1年內的軌道計算結果與MAR097 星歷相比，最大距離差異為900 m 左右，誤差圖像高頻振幅也隨時間緩慢增長，而軌道長半軸的差異不超過±0.25 m，其余軌道根數(shù)的差異量級也足夠微小。將散點放大后可以發(fā)現(xiàn)，軌道傾角、升交點赤經兩個軌道根數(shù)中，高頻部分相對的振幅較小且隨時間變化不大，而長半軸、離心率、近地點角距和平近點角的誤差圖像的振幅隨時間線性增長。

圖6 與MAR097相比1年軌道距離誤差Fig.6 Distance difference in 1 year compare with MAR097

圖7 與MAR097相比1年軌道根數(shù)誤差Fig.7 Orbital difference in 1 year compare with MAR09

通過圖8 和圖9 可知，將比較時段拓展到10年，最終距離誤差與一年的結果相比，誤差增長趨勢基本一致，圖像振幅逐漸減小，最大距離差異為8 km 左右。10年的軌道根數(shù)差異中，軌道傾角和升交點赤經差異開始呈現(xiàn)出周期性變化趨勢，差異都大致收斂在±10-5量級內，而長半軸、離心率、近地點角距和平近點角在誤差圖像的增長趨勢上與1年的比較結果相一致，僅在數(shù)值上相比有一個量級的增加。

圖8 與MAR097相比10年軌道距離誤差Fig.8 Distance difference in 10 years compare with MAR097

圖9 與MAR097相比10年軌道距離誤差Fig.9 Orbital difference in 1 year compare with MAR097

圖10 與圖11 為模型計算結果與星歷EPM2017在10年內的比較結果，與MAR09710年的比較結果有所不同的是，計算結果與EPM2017 的軌道距離誤差幅度變化較大，一年內差異達到了5 km，但是基本不隨時間增長，逐漸收斂。軌道根數(shù)中，近地點角距的誤差圖像呈現(xiàn)出周期性變化的趨勢，圖像中高頻振幅逐漸增大，與MAR097 對應的圖像存在差異。從數(shù)值上來看，長半軸、離心率、近地點角距和平近點角的差異與MAR09710年軌道根數(shù)比較結果相似，都在同一量級，但是升交點赤經和軌道傾角的差異增大了兩個量級。

圖10 與EPM2017相比10年軌道距離誤差Fig.10 Distance difference in 10 years compare with EPM2017

圖11 與EPM2017相比10年軌道距離誤差Fig.11 Orbital difference in 1 year compare with EPM2017

圖12與圖13為模型與星歷 INPOP17a 在 10年內的比較結果，其無論是軌道距離差異還是軌道根數(shù)差異都與MAR097的比較結果相似，差異也都隨時間大致呈線性增長。其根本原因是MAR097 中“火衛(wèi)1”星歷信息與INPOP17a 求得的“火衛(wèi)1”星歷互相差異本身就十分小，在10年內不到2 km，而且誤差也不隨時間增長，呈收斂趨勢。

可以看到，該數(shù)值模型得到的“火衛(wèi)1”軌道預報結果與MAR097、EPM2017 與INPOP17a 三種星歷對比，10年內距離誤差在10 km 以內。對比分析可知，計算結果對JPL 星歷MAR097 的擬合效果最好，誤差最小。主要原因是在動力學模型計算過程中涉及到的各天體的位置信息都從JPL 星歷DE430 中獲取。若要進一步對星歷EPM2017 和INPOP17a 進行擬合，可以改為從其各自對應歷表獲取其他天體位置信息。

圖12 與INPOP17a相比10年軌道距離誤差Fig.12 Distance difference in 10 years compare with INPOP17a

圖13 與INPOP17a相比10年軌道距離誤差Fig.13 Orbital difference in 1 year compare with INPOP17a

為了比較模型計算結果的精度，圖14 比較了以上3 種星歷各自之間在10年相同計算時段下“火衛(wèi)1”軌道的位置差異，可以看到，不同星歷之間1年內便出現(xiàn)了數(shù)千米的距離差異。但由于這些星歷都使用了實測數(shù)據(jù)對“火衛(wèi)1”進行精密定軌，所以其相互間的差異不會隨時間增長，并且在觀測數(shù)據(jù)密集的時間段誤差相對較小。從數(shù)值上來看，差異最大的是EPM2017 與INPOP17a，最大值6 km，最小值的是MAR097 與INPOP17a，僅為1.6 km。這是因為INPOP 系列星歷本身就與JPL 的DE 系列星歷在觀測資料方面體現(xiàn)出較好的一致性，所以計算出來的火衛(wèi)一星歷與MAR097 相差較小。而EPM2017 與其他兩者差異較大，這可能是觀測值數(shù)據(jù)處理和加權方案的不同引起的。

在未引入觀測值的條件下，該軌道預報模型的計算結果，在不同時段與不同來源的星歷比較后體現(xiàn)出的差異，和星歷自身之間的差異性相比，量級一致，僅在數(shù)值上稍大。可以充分說明該模型的預報結果在短期具有足夠的精度，長期也有一定的可靠性。

圖14 星歷MAR097，EPM2017和INPOP17a之間的“火衛(wèi)1”軌道距離差異Fig.14 Distance difference between MAR097，EPM2017 and INPOP17a

4 結 論

本文通過對“火衛(wèi)1”的動力學影響分析，確定了數(shù)值積分計算“火衛(wèi)1”軌道的動力模型，最終與不同星歷進行了比較分析。可以得到結論，該算法模型計算和預報“火衛(wèi)1”的軌道是切實有效、可靠的。實驗結果也存在需要改進的地方。從長時段的星歷比較結果來看，只利用動力模型計算出的數(shù)值結果和真實星歷存在10 km左右的差異，其中重要的原因是數(shù)值積分本身的誤差累計。真實星歷在數(shù)值積分的同時使用了許多對“火衛(wèi)1”的觀測值信息，對衛(wèi)星軌道進行改正，有更良好、穩(wěn)定的精度[20]。在動力模型部分，后續(xù)考慮“火衛(wèi)1”的天平動有望取得更好的計算結果[21]。

“火衛(wèi)1”的軌道預報，可用于研究其軌道的性質及其隨時間的變化。接下來的工作將以此為基礎，結合對“火衛(wèi)1”的觀測信息，特別是火星快車飛掠期間觀測的大量針對“火衛(wèi)1”的影像數(shù)據(jù)，建立高精度的星歷?，F(xiàn)有的火衛(wèi)歷表尚未處理這一數(shù)據(jù)。聯(lián)合“火星快車”地面深空站數(shù)據(jù)、“火星快車”對“火衛(wèi)1”的影像數(shù)據(jù)、地基“火衛(wèi)1”影像數(shù)據(jù)，進行“火星快車”和“火衛(wèi)1”聯(lián)合定軌，有望進一步改進“火衛(wèi)1”的軌道精度，進而在未來的“火衛(wèi)1”探測任務中發(fā)揮重要的作用。