１月地转移轨道设计约束条件分析
月心段返回轨道相对月心为双曲线轨道，但为了节省燃料，地心段轨道一般采用地心椭圆轨道。如图１所示。本文所使用的边界条件有：１）目标矢量ｒＴ。ｒＴ是从地球质心到月球影响球（ＳｐｈｅｒｅｏｆＩｎｆｌｕｅｎｃｅ）上一点的矢量，代表航天器出现在月球影响球上的位置。因此，该矢量与航天器穿过月球影响球的历元ＴＳｏＩ相关。２）再入点矢量ｒｉ。ｒｉ是地心到再入点的矢量。忽略地球半径在两极和赤道之间大约２０ｋｍ的差别及其他差异，可近似认为再入高度为一常数，本文将再入高度取为１２０ｋｍ。３）返回轨道相对地球赤道面的倾角ｉＲ。从ｒＴ到ｒｉ的轨道有２条，θ＜π或θ＞π，对应图１中的轨道１和轨道２，倾角分别为ｉＲ和１８０°－ｉＲ。由于逆行返回轨道相对大气的速度会增大，而顺行返回轨道则相反，因此逆行返回轨道不适合工程应用。４）再入面的飞行速度倾角γｉ。γｉ（即再入角）总是为负值，一般取绝对值。设计中要求其在一个范围（即再入走廊）内，γｉｍｉｎ＜γｉ＜γｉｍａｘ，其中γｉｍｉｎ为允许的最小再入角，γｉｍａｘ为允许的最大再入角。５）出月球影响球方向。必须选择是降轨出月球影响球还是升轨出月球影响球，否则将有２个不同的轨道面包含目标矢量。６）飞行圈数。从目标矢量到再入矢量可以单圈到达，也可以飞行Ｎ圈（Ｎ＞０）后到达。飞行圈数大于０的情况在航天器交会对接中经常使用。在月地转移轨道中，为了尽快返回，只考虑单圈返回的情况，即取Ｎ＝０。７）从月球影响球到再入的飞行时间ＴＦＲ。为了尽快返回地球，脉冲转移方式下，月地转移时间一般为３～５天。如果忽略月球重力场，并将月地转移轨道近似为半长轴等于地月平均距离（３８５０００ｋｍ）的椭圆轨道，那么该轨道半个周期约为５天（～１２０ｈ），该值可以看作是月地转移飞行时间上限的近似值［１１］。研究结果表明，当月球离地球最近时，从月球影响球开始的最大飞行时间大约为３．５天；而当月球离地球最远时，从月球影响球开始的最大飞行时间大约为４天。抛物线飞行时间可以作为ＴＦＲ的下限。可以用Ｂａｒｋｅｒｓ方程得到２个矢径（ｒＴ＜ｒｉ）间抛物线轨道飞行时间的表达式，见式（１）。该方程假设ｒＴ和ｒｉ定义的点之间没有近拱点。根据经验，对月球转移轨道而言，当月球离地球最近时ＴＦＲ大约为１．７天，最远时大约为２．１天。
２摄动项分析和轨道动力学模型
２．１摄动力及其量级对低轨（高度为１００～３００ｋｍ）月球卫星来说，所受的摄动源包含下列１０类，对应的摄动量级见表１［１２］。对高轨地球卫星，比如ＧＥＯ（ＧｅｏｓｔａｔｉｏｎａｒｙＥａｒｔｈＯｒｂｉｔ，地球静止轨道）卫星而言，所受的摄动力及其量级见表２。从上表可以看出，对于一般的轨道分析，只要考虑非球形引力摄动和第三体引力摄动，这２类摄动因素的影响可体现航天器轨道变化的主要特征。２．２轨道动力学方程根据摄动力及其量级分析结果和月地转移轨道的精度要求可以看出，在实际工作中，对于一个具体问题，并不需要考虑所有的摄动力因素，而是可以对具体的摄动力进行取舍，采用合适的有限精度模型。下面以月球影响球为界，采用２个受摄二体问题建立月地转移轨道的轨道力学运动方程。当月球探测器运动于月球影响球中时，月球引力作为主要力源，地球引力和太阳引力作为摄动力，根据摄动力的大小，月心段轨道运动考虑的摄动因素主要为地球引力摄动和太阳引力摄动，运动方程为¨ｒ＝－μｍｒ３ｒ－μｅ（ｒｅｄｒ３ｅｄ＋ｒｅｒ３ｅ）－μｓ（ｒｓｄｒ３ｓｄ＋ｒｓｒ３ｓ）其中μｅ，μｍ和μｓ分别是地球、月球和太阳的引力常数；ｒｅ，ｒｓ和ｒ分别为地球、太阳和月球探测器在月心天球坐标系中的位置矢量；ｒｅｄ表示探测器到地球的矢径；ｒｓｄ表示探测器到太阳的矢径；地球和太阳的星历从ＤＥ４０５中获取。当月球探测器出月球影响球后，地球引力成为主要力源，月球和太阳的引力为摄动力。根据摄动力的大小，地心段轨道考虑的摄动因素有地球非球形引力摄动、月球引力摄动和太阳引力摄动，建立如下运动方程¨ｒ＝－μｅｒ３ｒ－μｍ（ｒｍｄｒ３ｍｄ＋ｒｍｒ３ｍ）－μｓ（ｒｓｄｒ３ｓｄ＋ｒｓｒ３ｓ）＋ａｅ其中ｒｍ，ｒｓ和ｒ分别为月球、太阳和卫星相对于地心的矢径；ｒｍｄ表示探测器到月球的矢径；ａｅ为地球非球形摄动加速度；月球和太阳的星历从ＤＥ４０５中获取。
３月地转移轨道精确设计方法
３．１微分改正法以基于双二体模型得到的出月球影响球时刻和相应的位置、速度为初值，通过迭代修正出月球影响球的速度，使从月球影响球出发积分到再入时刻的位置与指定再入点的位置偏差满足精度要求。该迭代过程中，初始参数为出月球影响球的速度矢量，终端参数为再入点的位置矢量。终端参数可以表示成初始参数的函数，记为ｑ＝ｆ（ｐ）其中ｑ为终端参数，ｑ＝［ｑｘｑｙｑｚ］Ｔ为再入点的位置；ｐ为初始参数，ｐ＝Ｖ＝［ｖｘｖｙｖｚ］Ｔ为出月球影响球的速度。按线性摄动法，初始参数的小偏差Δｐ传播至再入时刻带来的终端参数的变化Δｑ可近似表示为Δｑ＝ｆｐΔｐ（２）记误差传递矩阵为Ｋ，Ｋ＝ｆｐ为３×３的矩阵，也可以表示为３个偏导数矩阵的乘积Ｋ＝ｆｐ＝ｑｘ１ｘ１ｘ０ｘ０ｐ其中ｘ０和ｘ１分别为探测器的初始和终端状态，即位置和速度。利用以上关系，一个初始参数可以得到对应的唯一一个终端参数。但对转移轨道来说，根本无法找出Ｋ的解析表达式，而我们的目的主要是得到它的具体数值。因此本文采用有限差分法进行求解，其基本原理就是在初始参数的每个分量上增加一个小量，计算终端参数ｑｎ。计算偏差Δｐｎ＝ｑｎ－ｑｔａｒｇｅｔ其中ｑｔａｒｇｅｔ为要求的目标终端参数；ｎ为迭代次数。解式（２），得到出月球影响球速度的修正量ΔＶΔｐ＝ΔＶ＝Ｋ－１Δｑ然后更新初始参数ｐｎ＋１＝ｐｎ＋Δｐ重复该过程，直到ΔＶ为零，得到最终解。３．２精确轨道设计的计算流程本文采用了一种双向嵌套循环搜索算法，求解同时满足两端约束条件的精确月地转移轨道。以出月球影响球的时刻和位置、速度为中间变量，一方面采用前向数值积分和微分改正法搜索满足地球再入端的轨道，另一方面采用后向数值积分并进行倾角和近月距修正得到满足月球端的轨道。由于月地转移轨道与月球影响球交点的位置初值是在无摄动的情况下得到的，这样，一旦用受摄模型反向积分到近月点时，得到的近月距和倾角将不是原先设定的月球停泊轨道近月距和倾角，所以必须对中间变量进行倾角和近月距修正。一般修正２～３次就可满足容差要求。通过这种双向嵌套循环搜索算法，可以求解同时满足两端约束条件的精确轨道，并使得两段轨道在月球影响球边界处的位置和速度连续，从而获得一条完整的满足两端约束条件的月地转移精确轨道。月地转移轨道精确设计的计算流程如图２所示。
４设计算例
４．１约束条件和参数设置选取２０１７年１月２６日出月球影响球作为返回窗口，下面给出该窗口下月地转移轨道精确设计的一个算例。表３为该算例对应的约束条件设置。根据表３给出的约束条件，以快速轨道设计结果作为计算初值，采用精确轨道动力学模型和数值积分法，通过３．２节的设计流程，最后可以得到４条满足两端约束条件的月地转移轨道，分别为：降轨出月球影响球顺行轨道、降轨出月球影响球逆行轨道、升轨出月球影响球顺行轨道和升轨出月球影响球逆行轨道。由于逆行轨道返回所需的速度增量远大于顺行轨道，速度增量差可以达到１１３ｍ／ｓ，因此在工程实际中不考虑逆行返回。４．２设计结果顺行的月地转移轨道入轨所需的速度增量小，适合工程使用，这里仅给出其中一条升轨出月球影响球的顺行精确轨道的最后计算结果，见表４。在ＳＴＫ中采用同样的模型对该轨道进行仿真，结果见图３～图６。图３给出了从１００ｍ高的环月停泊轨道经一次制动进入月地转移轨道的惯性轨迹，图４为该轨道在３Ｄ下的全景，图５和图６分别给出了地面和月面上的星下点轨迹。可以看出，ＳＴＫ的仿真结果与本文程序计算结果完全一致。５结束语本文以我国后续月球探测任务的实际工程需求为背景，以基于Ｌａｍｂｅｒｔ算法的快速轨道设计结果为初值，开展了精确轨道设计研究。通过对整个返回飞行阶段的摄动项和量级分析，建立了适合于月地转移轨道精确设计的轨道动力学方程。精确轨道设计以出月球影响球的时刻和位置、速度为中间变量，将轨道分为地心段和月心段分别进行计算。采用微分改正法和前向、后向数值积分的双向嵌套循环搜索算法，求解同时满足两端约束条件的精确轨道。通过这种双向嵌套循环搜索算法，获得分别满足月球端约束条件的月心段轨道和满足再入端约束条件的地心段轨道，并使得两段轨道在月球影响球边界处的位置和速度连续，从而获得一条完整的满足两端约束条件的月地转移精确轨道。最后以２０１７年１月２６日出月球影响球作为返回窗口，给出了具体的设计算例，并将程序计算结果在ＳＴＫ中进行了仿真验证，进一步证明了其正确性。
作者：郑爱武 周建平 单位：航天飞行动力学技术重点实验室 北京航天飞行控制中心 中国载人航天工程办公室
