木星探测器轨道计算中的动力学

曹建峰，黄 勇，段建锋，秦松鹤，张 宇，李 勰

(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。