彭 坤,孙国江,杨 雷,徐 明,奉振球
(1.中国空间技术研究院载人航天总体部,北京100094;2.北京航空航天大学宇航学院,北京100191)
相比于冒险进入遥远的太阳系深处,增加人类在地月空间的探索活动更加实际可行。同时,月球及其邻近空间可作为深空探测技术的天然试验场和载人深空探测的中转站[1],故月球探测活动的实施可为人类深空探测打下坚实的技术基础。地月转移轨道设计是实施月球探测任务的关键。不同于中心引力体环绕轨道,地月转移轨道涉及两个中心引力体,非线性强,设计难度较大。
目前,已有大量学者对地月转移轨道设计方法进行研究。 张汉清等[2]和何巍等[3]采用引入太阳引力摄动的弱稳定边界理论,在四体问题模型下设计了地月低能转移轨道。徐明等[4]则在四体模型基础上设计了小推力变轨方式的地月转移轨道。不管是地月低能转移轨道还是小推力地月转移轨道,其转移时间较长,通常为几十到几百天。考虑到转移时间约束,月球探测任务往往采用两脉冲变轨形式的地月转移轨道,在近地轨道进行加速,到达近月点后实施制动以进入目标环月轨道,其转移时间仅为3~5天[5-6]。对于两脉冲地月转移轨道,一般设计步骤为先利用低精度模型设计轨道初值,再在高精度模型中迭代搜索精确轨道。常用的低精度模型有双二体模型、圆型限制性三体模型和伪状态模型。双二体模型是将地月转移轨道分为地心段和月心段,在月球影响球处进行拼接,并以拼接点处轨道参数作为设计变量[7-8]。该方法中设计变量物理意义不明显,迭代搜索复杂,且地月转移过程动力学模型不连续,精度稍低。目前圆型限制性三体模型主要用于求解自由返回型的地月转移轨道,根据自由返回轨道对称性以近月点参数作为设计变量[9-10]。张科等[11]利用圆型限制性三体模型求解了双脉冲地月转移轨道,但其求得的地月转移轨道转移时间较大,不适合工程应用。伪状态模型本质上是限制性三体问题的封闭近似解法[12],可较好近似地球月球双引力场作用下的轨道推演,可用于一般地月转移轨道和自由返回轨道的初值求解[13-14]。 高玉东等[15]则采用分层搜索思想,利用B平面参数直接在高精度模型中搜索地月转移轨道,但仅针对极月目标轨道进行求解。杨维廉等[16]以地月转移轨道出发速度、升交点赤经和近地点幅角为设计变量,详细推导了地月转移轨道的转移矩阵和修正方法,分析了地月转移轨道的特性,但未给出具体的初值估计方法。贺波勇等[17]给出了地月转移轨道半长轴、升交点赤经和近地点幅角的初值估计方法,且算例中估计值非常接近精确值,其初值依赖近月点月球位置,没有反映出地月转移轨道特性。
由文献[18]知,地月转移过程中航天器所受的地球引力和月球引力影响最大,因此忽略地球和月球的非球形影响,以及太阳引力摄动影响,在三体模型下设计地月转移轨道不影响其本质特性。同时,根据文献[16]的研究,与月球轨道面共面的二维地月转移轨道和三维地月转移轨道特性类似,故可通过共面转移情况分析地月转移轨道基本特性,为三维地月转移轨道设计提供依据。本文建立三体模型下地心旋转系动力学模型,并以地月转移加速速度增量和地心旋转系对准角为控制变量,提出其二体模型下初值估计解析计算公式,针对不同运行方向地月转移轨道提出相应的搜索方法,采用序列二次规划算法对地月转移轨道进行求解,分析地月转移轨道的全局特性,以及地月转移轨道与设计变量之间的关系,其结论可应用于高保真动力学模型中的精确初值猜测,提高高保真模型下地月转移轨道的设计效率。
在三体模型下,假设地球、月球和航天器均被看作质点,月球绕地球做匀速圆周运动,航天器只受到地球和月球的万有引力影响。该三体模型假设与圆型限制性三体模型假设相同。不同于圆型限制性三体模型中将坐标原点选在地月系统中心,为揭示双脉冲地月转移轨道控制变量与转移时间的关系,本文以地心为中心建立动力学方程。航天器地月转移加速时刻的地月连线为x轴,由地心指向月心为正,z轴为月球公转运动角动量方向,y轴由右手定则确定,建立地心惯性系,则该坐标系下航天器动力学方程如式(1)所示:
式中,R为航天器相对地球的位置矢量,Rm为航天器相对月球的位置矢量,RME为月球相对地心的位置矢量,μE为地心引力,μM为月心引力。与地心惯性系类似,以地心为中心,地月连线为x轴,可建立地心旋转系。将式(1)从地心惯性系转移到地心旋转系,只考虑xy面内运动,并以极坐标系代替直角坐标系可得式(2):
式中,[r,θ,vr,vθ] 为地心旋转极坐标系统(图1)下航天器的极径、极角、径向速度和横向速度,ω为月球绕地球公转角速度,d为地月距离,Rm为航天器与月球的距离。
地月转移轨道初始状态为近地停泊圆轨道,满足轨道高度约束,施加切向脉冲Δv后进入地月转移轨道;终端状态为近月点,其月心航迹角γm为0,且需满足环月轨道高度约束。设地月转移加速时刻为0,终端时刻为tf,则其初始条件和目标约束可分别表示为式(3)和(4):
式中,RE和RM分别为地球和月球的半径,he和hm分别为近地停泊轨道高度和环月轨道高度,rm( tf)和vm( tf)分别为地心旋转系下终端时刻航天器相对月球的位置和速度矢量,其表达式如式(5)所示:
由系统模型可知,在初始条件中 θ0()和vθ0()未知。本文设置初始时刻地月转移加速切向脉冲Δv,以及极径r与地心转移系-x轴的夹角α(简称地心旋转系对准角,逆时针方向为正)为控制变量,其与初始条件关系如式(6)所示:
为提高搜索速度,将航天器相对月心的航迹角作为终止条件,即式(4)的第二分式。同时增加对航迹角变化率的判断,进一步确定终端点为近月点,如式(7)所示:
转移时间作为目标约束,则相应的目标约束变为式(8):
式中,T为目标转移时间。通过以上转换,可将地月转移轨道转化为非线性规划问题,控制变量为Δv,α[],目标约束如式(8)所示,状态方程如式(2)所示,积分终止条件为式(7)。
地月转移轨道大部分受地球引力影响,在初值估计时可将其假设为二体模型下的霍曼转移轨道,即一个近地距为RE+he,远地距为d+RM+hm的地心椭圆轨道。根据二体模型近地点速度计算公式可算得Δv的值。同时二体模型假设下,地月转移轨道拱线应对准近月点时刻的月球位置,则α应为航天器地月转移过程中月球在白道面内的扫角。由以上分析可给出控制变量的初值估计解析公式,如式(9)所示。考虑月球引力的影响,初值估计中适当增加地月转移加速速度增量值。
通过式(9)可计算出控制变量Δv,α[]的初值,代入求解模型中计算出目标约束式(8)的误差,通过修正算法不断修正控制变量,从而得到满足目标约束的地月转移轨道。本文采用成熟的序列二次规划算法SNOPT[19]对此非线性规划问题进行求解。
对于二维地月转移轨道,满足相同近地轨道高度、环月轨道高度和转移时间约束的轨道共有四条,分别对应地心顺行月心顺行、地心顺行月心逆行、地心逆行月心顺行和地心逆行月心逆行,如图2所示(与地心旋转系的旋转方向相同为顺行,反之为逆行)。图2(a)对应地心均为顺行方向,月心分别为顺行和逆行方向的两条地月转移轨道;图2(b)对应月心均为顺行方向,地心分别为顺行和逆行方向的两条地月转移轨道。
对于地心顺行和逆行的情况,可通过设置初始时刻横向速度vθ方向的方法进行区分,如式(10)所示:
对于月心顺行和逆行的情况,两者在大部分阶段轨迹近似相同,只是在靠近月球段略有不同,若直接采用SNOPT对地月转移轨道进行求解时无法进行区分。为明确地月转移轨道月心运行方向,本文设置月心顺行和逆行搜索策略来控制其搜索方向。比较月心顺行和逆行的地月转移轨道可知,月心逆行轨道的近月点在月球外侧,因此其地月转移加速Δv和地心旋转系对准角α的值相对较大。根据这一特性设计如下月心顺逆行搜索策略:先按标称程序进行搜索,若满足月心顺逆行要求,则直接输出结果;若不满足月心顺逆行要求,则对当前已搜到的控制变量进行修正,再采用SNOPT进行搜索,其修正公式如式(11)所示:
设地球和月球引力常数分别为 μE=398 600.4415km3/s和μM= 4902.8000 km3/s,地球和月球的半径分别为 RE=6378.140 km和RM=1738.000 km,地月间距离 d =384 400 km,月球公转角速度 ω = 2.649 07 ×10-6rad/s。 为验证地月转移轨道求解方法,设近地轨道高度he=200 km,环月轨道高度hm=100 km[6]。
将地月转移轨道设置为地心顺行月心逆行轨道,令转移时间为3天。采用本文提出的搜索算法进行求解,仿真结果如表1所示。由表1可知,搜索算法仅迭代3步就搜索到地月转移轨道,耗时仅为5.68 s,验证了算法的可行性和快速性。同时从表1可看出,地月转移速度增量Δv的初值估计和真实值几乎相同;对准角α的初值估计和真实值仅相差约9°,在0°~360°范围内看相差也较小,说明了初值估计方法的正确性。
表1 单个地月转移轨道搜索结果Table 1 Search results of single trans-lunar trajectory
通过坐标变换可绘制地月转移轨道分别在地心旋转系和地心惯性系的轨迹,如图3所示。由图3(a)可知,地月转移轨道大部分轨迹都在月球影响球之外;由图3(b)可知,地月转移轨道拱线近似对准近月点时刻的月球位置,验证了初值估计方法中二体模型假设的合理性,同时解释了为何初值估计与真实值误差较小。
令转移时间为3天,分别搜索地心顺逆行和月心顺逆行4种类型组合的地月转移轨道,其搜索结果如表2所示。由表2可知,对于地心段运行方向,地心顺行轨道的地月转移加速Δv比相应的地心逆行轨道略小,与初值估计近似相同;地心顺行轨道的对准角α比相应的地心逆行轨道大,且大于初值估计,而地心逆行轨道对准角α小于初值估计。该结果说明了地心顺行轨道要尽量向近月点时刻月球位置前方瞄准;而地心逆行轨道要尽量向近月点时刻月球位置后方瞄准。对于月心段运行方向,月心顺行轨道地月转移加速Δv比相应月心逆行轨道小,其对准角α也比相应月心逆行轨道略小。该结果说明了与月心顺心轨道相比,月心逆行轨道的对准角要略偏向前方,速度增量略大,验证了不同类型轨道搜索策略的正确性。
表2 4种类型地月转移轨道搜索结果Table 2 Search results of four types of trans-lunar trajectories
为了验证本文搜索算法的鲁棒性,对转移时间为1~7天的4种类型地月转移轨道进行搜索,每隔0.5天搜索一次,其搜索结果如图4所示。在搜索过程中,除个别轨道搜索时间略长,大部分轨道均在10 s以内收敛,说明该算法能适应大范围转移时间变化;同时搜索方向准确,能找到指定类型的地月转移轨道。
在4种类型地月转移轨道中,地心逆行轨道不能充分利用地球自转速度从而降低火箭运载能力,一般不予采用。对于载人月球探测任务,其地月转移轨道要求设计为自由返回轨道,而自由返回轨道的近月点为月心逆行轨道。以下主要对地心顺行月心逆行的地月转移轨道特性进行分析,为载人月球探测地月转移轨道选取和设计提供依据。
转移时间是地月转移轨道的关键设计因素[20],其决定了地月转移加速速度增量和近月制动速度增量的大小。取T=[1,7]d,近地轨道高度he=200 km,环月轨道高度hm=100 km,分析此转移时间范围内地月转移轨道两端点速度增量大小,如图5所示。随着转移时间的增大,地月转移加速和近月制动速度增量逐渐减小,到3天后两速度增量变化趋于平缓,最小速度增量出现在5天左右。故对于无人地月转移轨道,要求变轨总速度增量最小,其飞行时间可取为5天转移时间;对于载人地月转移轨道,要综合考虑速度增量和转移时间,可取为3天转移时间。
在初值估计过程中,地月转移加速速度增量和对准角均基于二体模型进行估值,未考虑月球引力影响。以下比较初值估计与真实值之间的误差,分析不同转移时间下月球引力对二体模型初值估计的修正量,如图6所示。对于地月转移加速速度增量,转移时间大于3天时,初值估计与真实值相差不大,且均可近似为常值;转移时间小于3天时,转移时间越小,地月转移加速速度增量真实值与初值估计相差越大。故当转移时间小于3天时,月球引力对地月转移加速速度增量影响大,初值估计时需要进行大幅修正。对于对准角,转移时间小于2天时,对准角随转移时间递减,真实值与初值估计相差较大;转移时间大于2天时,对准角随转移时间递增,真实值与初值估计误差逐渐减小,到转移时间为3天时,误差在10°之内。因此,对于转移时间3~7天的地月转移轨道,地月转移加速速度增量近似维持不变,对准角与转移时间近似成线性关系本文的初值估计方法能对控制变量进行高精度估计。
近地点高度和近月点高度决定了月球探测任务近地停泊轨道和环月轨道的高度,其选取除考虑火箭运载约束和环月轨道稳定性外,还应考虑变轨速度增量大小。以T=3 d为例,分析不同近地点高度和近月点高度下地月转移轨道地月转移加速和近月制动速度增量变化情况。为提高收敛速度,根据转移时间影响分析结果,对对准角α初值估计进行修正,其修正公式为α=ωT+8°()。
设近地点高度he=[200,1000]km,近月点高度hm=100 km,地月转移轨道变轨速度增量变化曲线如图7所示。随着近地点高度增加,地月转移加速速度增量逐渐减小,变化范围约为200 m/s;近月制动速度增量也逐渐减小,但变化范围仅为3 m/s,可近似认为不变。
设近月轨道高度hm=[100,5000]km,近地轨道高度he=200 km,地月转移轨道变轨速度增量的变化曲线如图8所示。随着近月点轨道高度增加,地月转移加速速度增量逐渐增大,但变化范围仅为4 m/s,可近似认为不变;而近月制动速度增量逐渐减小,变化范围约为150 m/s。
本文在三体模型下建立了地心旋转系地月转移轨道求解模型,给出了控制变量初值估计解析式,提出针对地心顺逆行和月心顺逆行4种类型轨道的搜索策略,采用序列二次规划算法进行了求解。仿真结果表明该方法搜索速度快,鲁棒性好,能适应1~7天转移时间内不同类型地月转移轨道的搜索。当转移时间范围为3~7天时,初值估计方法能够对真实值进行准确估计。同时,根据初值估计与真实值的误差,能够得到月球引力对二体模型解析解的修正量,为后续高保真模型中地月转移轨道设计提供更精确的初值。
转移时间增加会导致地月转移轨道地月转移加速和近月制动速度增量的减少,当转移时间为3天时趋于平缓,可将3天~5天作为月球探测任务中地月转移轨道转移时间设计约束。近地点高度增加会导致地月转移轨道地月转移加速和近月制动速度增量的减少,主要影响地月转移加速速度增量,近月制动速度增量几乎不变。近月点高度增加会导致地月转移轨道地月转移加速速度增量的增加和近月制动速度增量的减少,主要影响近月制动速度增量,地月转移加速速度增量几乎不变。