木星探測器軌道計算中的動力學

曹建峰，黃 勇，段建鋒，秦松鶴，張 宇，李 勰

(1.航天飛行動力學技術重點實驗室，北京100094；2. 北京航天飛行控制中心，北京100094；3. 中國科學院上海天文臺，上海200030；4. 中國科學院大學，北京 100049)

0 引 言

木星是太陽系中體積最大、自轉最快的氣態巨行星，質量大約是地球的318倍。對木星的系統性探測有助于更好的了解太陽系起源，以及大行星的形成演化過程。至今成功的木星探測活動共計9次，其中7次飛越探測，2次專門探測。飛越探測始于1972年發射的先驅者10號，它于1973年12月3日發回了第一組木星的近距離拍攝圖像[1]。對木星的2次專門探測為1989年發射的伽利略號飛船[2]， 2011年發射的朱諾號(JUNO)[3]。通過對木星的多次探測，基本認識了木星的系統起源、內部結構、大氣層以及磁層等[4-6]。但是相對于太陽系其他大行星的認知，對木星的研究仍顯任重道遠。

美國、歐洲將木星探測作為深空探測活動的中長期戰略目標。2008年，NASA和ESA組成了木星系探測聯合研究組，并提出在2020年后實施“木衛二木星系統任務”和“土衛六土星系統任務”，美歐將集中資源聯合開展木星系、土星系探索任務。2012年5月，ESA確定了將在2022年發射“木星冰月探測器”，探索木星衛星存在生命的可能性。

《2016-2030空間科學規劃研究報告》明確提出了木星探測計劃，以帶動我國空間技術、科學及科學應用的發展。軌道動力學是木星探測任務測控系統必須掌握的一項關鍵技術，直接決定飛行任務能否順利實施；是開展軌道計算、任務規劃的基礎，為科學應用提供探測器精確空間位置信息。2018年，陳略等[7]組織中國深空網開展了木星探測器開環測量試驗，提取了噪聲水平約為10 mHz的多普勒測量，為木星探測任務的測軌積累了有益經驗，深空網的系統建設也推動著著自主木星探測的進程[8-9]。

現有木星探測任務主要由美國國家航空航天局(NASA)實施，實測數據在美國噴氣推進實驗室(JPL)處理。當前的公開文獻在木星探測的動力學問題上都是簡單提及，不具備實操性，距離定軌軟件的研制、工程任務的實施仍有很大差距。本文針對自主木星探測任務軌道計算問題，系統性的梳理了環木星探測器的受力問題，詳盡的整理了相關的計算公式及具體參數，給出了繞飛階段軌道計算需要考慮的時空參考系、動力學模型，并利用木星探測器JUNO的星歷數據進行驗證。

1 木星概況

木星是太陽系最大的行星，與地球存在很大的差異。木星距離太陽更為遙遠，自轉周期更小；木星的直徑是地球的11倍，但是密度較小；木星有著太陽系內最厚的行星大氣層，跨越的高度超過5000 km(由于木星不存在固態表面，大氣層基準通常定義是大氣壓力等于1 bar之處)；木星的質量是太陽系其他行星質量總和的2.5倍，該特征使得太陽系的質心落在太陽表面之外[10-11]。表1～2給出了木星的軌道參數及基本物理參數。

木星擁有眾多的衛星，至2018年底發現的天然衛星達79顆，其中最大的衛星有4顆：木衛一(Io)，木衛二(Europa)，木衛三(Ganymede)，木衛四(Callisto)。這4顆衛星是伽利略于1610年首次發現，因此也稱為伽利略衛星，表3給出了幾顆木星衛星的軌道和物理參數。

表1 木星的軌道特性Table 1 Jupiter’s Orbital Characteristics

表2 木星的物理參數Table 2 Jupiter’s Physical Parameters

表3 木星主要衛星的軌道和物理參數Table 3 Jupiter Satellite Orbital and Physical Characteristics

木星具有強大的磁場，在太陽風作用下形成輻射帶。因此，木星探測需要考慮輻射帶分布的特點，設計盡量避開輻射帶的軌道開展探測活動。如JUNO的設計軌道近木點位于木星赤道附近，減輕輻射影響。不同于地球，木星內部釋放的熱能與木星接收的太陽輻射能量相當[12-13],對于環繞探測器，木星反照壓及紅外輻射壓也更為顯著。

2 木星探測器軌道動力學

2.1 時空參考系

木星探測器定軌計算中，軌道動力學面臨的首要問題是相對論時空參考系的使用。木星探測器定軌中涉及到的時間坐標系包括協調世界時(Universal time coordinated,UTC)，質心動力學時(Bargcentric dynamic time,TDB)，地球時(Terrestrial time,TT)等，坐標系統包括地心天球坐標系，地球參考系，質心天球坐標系，木星質心天球坐標系和木星固聯坐標系。和木星相關特有的時空坐標系有木星動力學時(Jupiter dynamical time, TDJ)、木星質心天球坐標系和木星固聯坐標系，下面進行具體描述。

嚴格意義上，類似于地球有必要建立一套局部的木星天球參考系，時間系統為TDJ。Hellings給出的Lorentz變換式可進行太陽系質心時空參考系中的空間坐標到局部參考系中的空間坐標的轉換[14]。截斷至1階參數化后牛頓項(PPN)，TDJ與TDB的關系式為，

(1)

式中:c表示光速，v表示木星公轉速度，mk表示木星之外其他天體(k)的質量，rk表示天體(k)與木星的距離。

圖1給出了以2000- 01- 01為起點，50年內TDJ與TDB的差異(TDJ與TDB的起點均為2000- 01- 01T00:00:00.000)，TDJ除了木星公轉相關的周期項外，還存在較大的線性項，50年的長期漂移達到近5 s。另外，TDJ-TDB的周期項的幅值也較TT-TDB高出近1個量級。

圖1 TDJ與TDB的差異Fig.1 Difference between TDJ and TDB

木星探測器的軌道坐標(Jovicentric coordinate)在木星局部參考系中描述，其與質心天球參考系(Barycentric celestial reference system, BCRS)的轉換關系為，

(2)

位置轉換關系對時間進行求導，可以計算速度轉換關系式為，

(3)

2.2 定向參數模型

木星引力場在木星固聯坐標系中描述，而軌道積分需采用的是局部慣性坐標系(木星天球參考系)，因而力模型的計算需要解決木星天球參考系與木星固聯參考系的相互轉換，即木星的定向參數模型。

IAU地圖坐標及旋轉參數工作組(Working group on cartographic coordinates and rotational elements) 大概每3年發布1次太陽系大天體的指向參數模型報告。報告中采用國際天球參考系(ICRS)，通過簡單關系式將各天體的天球坐標系與固聯坐標系聯系起來[15-17]。

圖2描述了木星天球參考系與木星固聯坐標系之間的轉換關系。這里涉及到3個坐標系。

1)木星固聯坐標系(Jupiter-centered jupiter mean equator and prime meridian of date)

木星固聯坐標系為星體坐標系，基于IAU木星定向參數模型定義。木星平赤道面為固聯坐標系的參考平面，基本方向為木星本初子午線與參考平面的交點。IAU指定在木星平赤道平面上自西向東，從Q點開始至本初子午線角度為W。Q點定義為木星平赤道相對于地球平赤道的升交點。

W=284.95+870.5360000d

式中:d為從J2000.0歷元起算的地球天數。

2)木星平赤道坐標系(Jupiter-centered jupiter mean equator and IAU-vector of epoch)

木星歷元平赤道為坐標系的參考平面，IAU矢量為參考方向，IAU矢量由木星質心指向木星歷元平赤道與地球J2000.0歷元平赤道面的交點Q。

3)木星天球參考系(Jupiter-centered earth mean equator and equinox of epoch)

木星天球坐標系坐標原點選取為木星質心，參考平面為地球歷元平赤道，參考方向為歷元平春分點方向。該坐標系與地心天球參考系完全對應，不同之處為坐標原點由地球質心平移至木星質心。

圖2 木星坐標系轉換關系Fig.2 Relationship between prime meridian and IAU-vector

IAU定義的木星天極在J2000.0地球平赤道平春分點坐標系中的方向，其形式為，

α=268.056595-0.006499T+0.000117sinJa+

0.000938 sinJb+0.001432 sinJc+

0.000030 sinJd+0.002150 sinJe

δ=64.495303+0.002413T+0.000050cosJa+

0.000404cosJb+0.000617cosJc-

0.000013cosJd+0.000926cosJe

(4)

(5)

式中：T表示從J2000.0起算的儒略世紀數(36 525 d)，α,δ的單位為(°)。

木星天球參考系至木星平赤道坐標系轉換關系為，

rIAU=Rx(90-δ)Rz(90+α)rcrs

(6)

木星平赤道坐標系至木星固聯坐標系的旋轉矩陣，

Rbf=Rz(W)rIAU

(7)

2.3 動力學模型

木星環繞型探測器的軌道計算本質上仍然是一個受攝運動二體問題。環繞木星飛行的探測器受到的主要作用力是木星產生的質點引力，其余各類作用力相對于木星質點引力都是一個小量。環繞探測的運動方程可以表示為[18]，

(8)

式中:r表示探測器在木星天球參考系中的位置矢量，p表示其他攝動源產生的攝動加速度，μJ為木星引力常數。

木星探測器的主要攝動力包括：

1)非球形引力

木星非球形引力勢可以展開為諧系數表達形式，

[Clmcos(mλ)+Slmsin(mλ)]Plm(sinθ)

(9)

式中:G為萬有引力常數，M為木星質量，RJ為木星的赤道半徑，(λ,θ)為航天器在木星固聯坐標系下的經緯度，(Clm,Slm)為木星重力場系數，Plm為締合勒讓德多項式。

當前發布的最新木星重力場模型基于JUNO測軌數據解算，其系數如下[19]:

表4中給出的是非歸一化諧系數，歸一化諧的系數可以根據如下公式計算，

(10)

表4 木星重力場模型系數Table 4 Jovian gravity field parameters

(11)

2)第三體引力攝動

太陽、大行星及木星的自然衛星所產生的加速度統稱為第三體引力攝動，在木星天球參考系中，第三體引力加速度為，

(12)

式中:m′為第三體質量，r′為第三體相對木星的位置矢量，Δ為航天器相對第三體的位置矢量，即Δ=r-r′。

木星擁有眾多的衛星，幾個主要的自然衛星會對木星環繞探測器產生較大的攝動力，對于探測器的精密軌道計算，有必要考慮其攝動影響。網站(https://ssd.jpl.nasa.gov/horizons.cgi)提供了各小天體的軌道參數。

3)太陽輻射壓

太陽電磁輻射于航天器表面，產生的作用力稱為太陽光壓，與航天器距離太陽的距離及航天器受照部件的屬性相關。航天器各部件對太陽輻射具有吸收與反射能力，其作用力可以分為3部分：吸收太陽輻射產生的作用力為dFa=-KrsuncosθdS，鏡面反射產生的作用力為，

dFs=-2Kcos2θ(1-)(1-D)ndS

(13)

漫反射產生的作用力為，

dFd=-Kcosθ(1-)

(14)

式中:dS為表面積，K為1AU處太陽的輻射流量，θ為太陽矢量rsun與面元法向n在機械坐標系下的夾角。為表面吸收系數，D為表面反射率。1AU處太陽的輻射流量的通常取值為K=4.56×10-6N·m-2。

對于工程定軌計算，太陽輻射壓產生的加速度可以采用簡化的固定面質比模型，

(15)

式中:Δs為由太陽指向航天器的方向矢量，A為等效面積，Cr為太陽輻射壓系數，aU為天文單位。

4)木星反照壓和紅外輻射壓

木星接收太陽輻射的能量后，為保持自身的熱平衡狀態，會以兩種不同的方式將輻射能量釋放出去：光學輻射和紅外輻射。紅外輻射是木星表面以長波形式的二次輻射。木星反照壓與紅外輻射壓所產生的加速度取決于探測器可視的木星表面區域。

木星表面輻射區域p所產生的反照加速度計算公式為，

(16)

A一般以球諧系數展開，以表示不同的表面區域。

(17)

式中:(λ,θ)為輻射區域P的經緯度，Plm為締合勒讓德多項式，

文獻[11]給出的各系數參考值為，

紅外輻射壓產生的加速度，與反照壓具有相同的形式，可以表示為，

(18)

式中:各符號的意義同太陽反照壓。E為木星輻射系數，以諧系數表示。

(19)

5)潮汐動力學攝動

在太陽、大行星、以及木星自然衛星的引力位作用下，木星內部質量分布隨著時間發生變化，改變木星的引力場。根據“流體靜力平衡”理論，當天體外部存在引潮力時，天體將發生形變，使得天體內部在重力、彈性力、粘滯力和引潮力的作用下處于平衡狀態，稱為平衡潮。并假定，這種形變使天體內部等勢面發生變化。這種流體靜力平衡理論是一種理想狀態，與真實的形變還有一定差別。勒夫數就是平衡潮與真實的潮汐形變之間的比例系數，潮汐所產生的引力勢可以通過勒夫數表示，文獻[3]給出的k2=0.7，文獻[18]給出的k22=0.625。

(20)

潮汐攝動產生的加速度一般與非球形引力攝動加速度一起計算，可以展開表示為諧系數的形式，該諧系數是對木星重力場系數的修正。

(21)

6)相對論效應攝動

廣義相對論效應是對牛頓力學的修正，木星探測所采用的木星天球參考系基于廣義相對論建立，因而與牛頓理論所定義的木星非旋轉坐標系仍存在差異。在衛星動力學方面，相對論效應使得探測器的運動方程增加了相對論效應加速度，包括：Schwarzschild項，測地歲差，以及Lense-Thirring歲差。相對論效應攝動引起的加速度為，

(22)

3 動力學模型檢驗

JUNO由是美國“新疆界計劃”的第2個項目，于2011年發射，歷經5年進入木星軌道，其運行管理由JPL負責。JUNO位于近極軌的大偏心率橢圓軌道上，軌道周期大約為53天，近木點距離大約為木星半徑的1.06倍。

為了驗證力模型的正確性，選用JUNO探測器的軌道數據進行動力學擬合。圖3給出了JUNO在軌飛行約一個軌道周期內的高度與速度變化，圖4為近木點附近的局部放大圖。近木點高度3400 km，速度大約為58 km/s。

圖3 JUNO一個軌道周期近木點距離及速度的變化曲線Fig.3 Variation of height and velocity for one orbital period

圖4 JUNO近木點距離及速度的變化曲線Fig.4 Variation of height and velocity during perijove pass

JPL在利用JUNO的測軌數據反演重力場時，使用近木點前后大約6～8小時的數據進行定軌計算。本文選取2017年5月9日2時至10時共計8小時的弧段的星歷進行擬合，位置偏差小于10m，速度偏差在近木點小于6 mm/s(如圖5所示)。選取2017年4月23日之后連續10天的星歷進行擬合，位置偏差小于50m，速度偏差小于1 mm/s(如圖6所示)。JPL提供的星歷文件，每隔2～3天速度數據會出現0.1～0.2 mm/s跳變，這可能是星歷文件生成中考慮了隨機脈沖進行誤差補償，在本檢驗程序中未考慮該速度變化項，因而軌道擬合效果會降低。

圖5 近木點星歷擬合后的位置、速度偏差Fig.5 Post-fit deviation of position and velocity during perjove pass

圖6 星歷擬合后的位置、速度偏差Fig.6 Post-fit deviation of position and velocity

4 結 論

回顧木星的基本狀況，分析木星探測器軌道計算所涉及的時空參考系，動力學模型，并給出了木星坐標系的相互轉換關系，木星探測器軌道計算所需動力學攝動的具體數學模型。利用JUNO探測器的星歷數據進行力模型檢驗，近木點擬合位置偏差小于10 m，速度偏差小于6 mm/s；利用遠木點10天數據進行擬合，位置偏差小于50 m，速度偏差小于1 mm/s。