一种常值推力作用的轨道机动优化设计策略

刘将辉，彭祺擘，张海联，周建平

（1.中国航天员科研训练中心，北京 100094；2.中国载人航天工程办公室，北京 100071；3.北京航空航天大学宇航学院，北京 100083）

0 引言

近年来，随着航天应用技术和深空探测的迅速发展，作为实现空间复杂任务的基础，连续推力作用下的航天器轨道机动优化设计问题在空间探测中具有越来越重要的意义［1-3］。

针对连续推力的航天器轨道机动问题，相关研究主要集中在常值推力和非常值推力这两种不同类型推力作用下的轨道机动方式。常值推力作用的轨道机动指航天器在大小恒定、方向可控的推力作用下实现轨道转移的机动方式［4-6］。文献［7］研究了连续常值推力在轨道机动中的应用，该研究的轨道机动方式仅限于共面轨道转移和共面轨道逃逸。文献［8］研究了二体条件下连续常值径向推力作用的轨道机动问题，得到了航天器在常值径向推力作用时沿圆轨道飞行的条件。文献［9］研究了多段常值推力的水滴悬停轨道控制问题。文献［10］研究了常值推力作用的航天器最优交会问题，通过分段多项式描述航天器的状态量和控制量。文献［11］采用间接法研究了连续常值推力作用的共面圆轨道最小时间轨道转移问题。文献［12］采用改进间接法研究了常值推力作用下航天器共面轨道转移的燃料优化问题，给出了最优轨迹的方向角约束条件。文献［13］分析了航天器在常值径向推力作用下的轨道运动特性。文献［14］研究了考虑通信窗口约束的常值推力椭圆轨道最优交会问题。目前有关常值推力作用的轨道研究虽然取得许多成果，但局限于面内轨道机动，使得其存在应用范围过窄的缺陷。

非常值推力作用的轨道机动研究方面，文献［15］研究了连续推力作用的人工太阳同步轨道设计问题。文献［16］分析了连续小推力作用下航天器轨道机动动力学特性。文献［17-18］将轨道转移问题转化成凸优化问题，提出了一种分布式空间系统自主重构的新方法，该方法可确保航天器编队的安全性，同时可优化燃料消耗。针对近地圆轨道航天器编队重构的轨道机动问题，文献［19］提出了考虑速度增量约束的多脉冲机动规划算法，文献［20］提出一种基于凸优化方法的最优轨迹规划方法。文献［21-22］研究了连续小推力轨道机动方式在空间碎片清除任务中的应用。文献［23］分析了将大型空间碎片移轨/转轨到处置轨道的两种策略效率。文献［24］研究了面向小行星探测的连续小推力变轨问题，提出了采用自适应交互式多模型无迹卡尔曼滤波算法。文献［25］研究了微推力对微纳卫星轨道机动的影响，文献［26］分析了小推力和大气阻力对近地轨道机动的影响。文献［27］研究了圆轨道小推力卫星的轨道机动与燃料优化问题，文献［28］研究了采用电推进的小推力地心轨道转移的燃料优化问题。针对静止轨道卫星的轨道机动问题，文献［29］提出一种4 脉冲远程快速交会规划方法。文献［30］研究了航天器在近圆轨道之间机动的混合推力问题，提出了一种改进的交互多模型算法。文献［31］首次提出基于逆向设计的形状法，假设航天器沿某一类型的曲线飞行，通过曲线拟合法求解方程中的待定系数，从而得到航天器飞行轨迹的解析表达式，实现对连续小推力航天器轨道机动问题的快速计算。文献［32］采用正弦指数模型设计了小推力星际转移轨道。针对带推力约束的航天器空间轨道机动问题，文献［33］提出了有限傅里叶级数模型，文献［34］提出了一种改进的形状设计方法。文献［35］将轨道形状方程与贝塞尔曲线方程应用到平面轨道机动问题中，提出了一种基于贝塞尔曲线的形状设计方法。

相比于常值推力作用的轨道研究，非常值推力研究拓展到三维空间的轨道机动，但存在3点不足：1）工程实际中航天器的推力大小一般是恒定的。2）空间轨道机动的研究大多采用脉冲变轨这种理想化的变轨方式，该方式仅适用于机动时长占总飞行时长比例较小的情况。3）目前非常值推力的研究大多采用基于逆向设计的形状法，该方法适合于面内圆轨道之间的轨道转移，当轨道偏心率和始末状态轨道倾角均相差较大时，相对固定的拟合曲线导致形状法很难收敛到可行解。

为了解决上述不足，本文在柱坐标系上建立了常值推力作用的航天器三维空间轨道机动动力学模型，在推力大小恒定的前提下，通过调整推力在空间中的俯仰角和方位角实现航天器在三维空间的轨道机动。本文采用正向设计法，提出了一种先轨道等待再轨道机动的优化设计策略，该策略能在轨道偏心率和始末状态轨道倾角均相差较大的情况下实现航天器在三维空间的轨道机动任务。本文方法为三维空间轨道机动优化设计问题提供了新的思路。

1 轨道机动动力学模型

首先定义空间运动相关坐标系，如图1 所示。定义EJ2000 坐标系OXIYIZI的原点O为地心，基准平面为历元平赤道面，OXI轴由地心指向平春分点ϒ，OZI轴与基准平面的法向平行，OYI轴的方向通过右手法则确定，OXIYIZI为静系。定义参考坐标系ORXRYRZR的原点OR为航天器质心，ORZR轴与OZI轴平行，ORXR轴在矢量p和ORZR轴构成的平面内且与ORZR轴垂直，ORYR轴的方向通过右手法则确定，ORXRYRZR为动系。定义柱坐标系Orθz的原点为地心，极轴χ与OXI轴重合，以历元平赤道面为基准面，r在基准面内且与极轴χ的夹角为θ，z与OZI轴平行。

图1 空间运动坐标系之间关系Fig.1 Relationship between spatial motion coordinate systems

假设航天器的位置矢量p和速度矢量v在EJ2000 中的投影分别为pI=[pxI，pyI，pzI]T速度矢量为vI=[vxI，vyI，vzI]T=由图1可知：

对式（1）求导可得：

根据位置矢量p的绝对导数与相对导数关系可得：

式中：S(ω)=为ω=[ω1，ω2，ω3]T的反对称矩阵。dpdt表示在静系中矢量p对时间t的绝对导数表示在动系中矢量p对时间t的相对导数，ω表示动系相对于静系的角速度矢量。v和p在参考坐标系ORXRYRZR中的投影如图2所示。

图2 相关矢量在参考坐标系的投影Fig.2 Projections of the relevant vectors in the reference coordinate system

根据式（3）可得：

根据速度矢量v的绝对导数与相对导数关系可得：

式中：m为航天器的质量，p=‖p‖=为航天器质心到地心的距离，μ为地心引力常数，F为航天器所受控制力。F在参考坐标系ORXRYRZR中的投影如图2 所示，其中φ为俯仰角，η为方位角。根据式（5）可得：

式中：F=‖F‖。航天器质量变化率为：

式中：ve表示航天器推进系统喷气速度。根据式（4）、式（6）和式（7）可得到航天器轨道机动动力学模型：

式中：X(t)=(r，θ，z，vr，vt，vz)描述了航天器的状态。只需要知道航天器在惯性系中描述的初始状态和终端状态，通过式（1）和式（2）将其转换为柱坐标系中描述的初始状态X(t0)和终端状态X(tf)，利用优化算法对控制力F、俯仰角φ和为方位角η进行优化求解，便可实现航天器轨道机动优化设计。在工程实际中航天器的推力一般是恒定的，F和ve是常值，质量变化率也为常值。本文假设推力方向与航天器体轴平行，通过优化航天器飞行过程中的俯仰角φ和方位角η实现轨道机动优化设计。

2 参数无量纲化

式（8）描述的轨道机动动力学模型中各参数存在数量级差别，如果直接求解容易产生奇异，且影响计算效率。因此，需要对式（8）中各参数进行无量纲化处理。本文取航天器的初始质量m0为质量基本量纲，取地球半径rE为长度基本量纲。力、速度和时间的基本量纲分别为和角度的基本量纲为1。式（8）中各参数经无量纲化处理之后的表达式为：

对式（9）～式（14）分别进行微分可得：

根据式（8）～式（21）可得到无量纲化的航天器轨道机动动力学模型：

3 优化设计策略

航天器轨道机动优化设计问题的实质是寻找一条从初始状态到目标状态的飞行轨迹，在满足相关约束条件的同时使得航天器燃料消耗最低。传统的轨道机动优化设计策略研究的是直接轨道机动问题，航天器从初始状态到目标状态的飞行过程中全程施加推力。本文提出一种先轨道等待再轨道机动的优化设计策略，航天器先在初始轨道上无动力飞行一段时间进行轨道等待，然后寻找一个最优时刻施加推力再进行轨道机动。与传统的轨道机动优化策略相比，本文不仅需要寻找飞行过程中的最优控制量，还需要寻找最优等待时刻。本文优化设计的数学描述如下式所示：

图3 优化设计策略Fig.3 Optimal design strategy

本文将多重直接打靶法（Multiple direct shooting method）和内点法相结合对优化问题进行具体求解。与直接打靶法相比，多重直接打靶法具有计算精度高且收敛性强的优点。其思路是将积分时间段划分成若干个子时间区间，将各个子时间区间的待优化量进行离散化处理。已知初始状态可以依次求出各子时间区间的状态，中间每个子时间节点的状态既是上一个子时间区间的末端状态又是下一个子时间区间的起始状态。这样，将原先复杂的轨迹优化设计问题转化为每个子时间区间的非线性规划（Non-linear programming，NLP）问题，通过内点法求解得到最优飞行轨迹、最优俯仰角、最优方位角和

4 仿真算例

本文以航天器从近地小椭圆轨道异面机动转移到远地大椭圆轨道为例进行仿真分析，航天器进行轨道机动的始末端轨道根数如表1所示。航天器质量为3 000 kg，引力常数μ=3.986 × 1014m3s2，长度的基本量纲rU 为地球半径。时间和速度的基本量纲分别为tU=806.785 6 s和vU=7 905.4 m s，质量的基本量纲为mU=3 000 kg，力的基本量纲为fU=29 396 N。无量纲化后的航天器推力和喷气速度分别为0.1 fU 和6.202 2 vU。将表1 中机动轨道始末端轨道根数转化为EJ2000 坐标系描述的始末端状态，通过坐标转换并进行无量纲化处理，得到表2所示的柱坐标系描述的无量纲化始末端状态。

表1 机动轨道始末端轨道根数Table 1 Orbital elements at the beginning and end of the maneuvering orbit

表2 无量纲化始末端状态Table 2 Dimensionless initial states and end states

本文设置了两组仿真场景以验证所提轨道机动优化设计策略，即：

1）先轨道等待再轨道机动

2）直接轨道机动

场景一：将轨道等待段和轨道机动段各划分为10个子时间区间，每个子区间设置30个离散点。俯仰角和方位角的上界和下界分别设为π 和-π。在[ -π，π]范围内，俯仰角和方位角的初值选取不影响多重直接打靶法求解的收敛性。本文中，各离散点的俯仰角和方位角初值均设为1。时刻和时刻的初值选取对多重直接打靶法求解的收敛性有影响，时刻的初值选取根据初始轨道周期来估计，时刻的初值选取根据终端轨道周期来估计。本文中，时刻和的初值分别设为5 和10，其上界分别设为7.884 2 和20.379 2，下界均设为0。仿真结果如图4～图7所示。

图4 三维空间中的位置变化Fig.4 Position changes in three-dimensional space

图4为航天器在三维空间中的位置变化，由图4可知，航天器在初始轨道先飞行一段时间进行轨道等待，在t1时刻开始进行轨道机动，t2时刻到达目标位置。图4 表明，在轨道偏心率和始末状态轨道倾角均相差较大的情况下，航天器依然能实现三维空间的轨道机动任务。

图5为航天器速度随时间的变化关系，在[t0，t1)时间段不施加作用力，航天器进行无动力飞行。在[t1，t2]时间段施加常值推力，航天器进行轨道机动，总机动时间为5.788 6 tU，整个任务的总飞行时间为8.108 tU。

图5 速度随时间的变化Fig.5 Variation of velocity over time

图6为航天器质量随时间的变化关系，在[t0，t1)时间段航天器质量不变，在[t1，t2]时间段航天器质量匀速减小，在t2时刻航天器质量为0.906 6 mU。

图6 质量随时间的变化Fig.6 Variation of mass over time

图7 为俯仰角和方位角随时间的变化关系，在常值推力条件下，航天器通过控制俯仰角和方位角实现三维空间轨道机动。

图7 俯仰角和方位角随时间的变化Fig.7 Variations of pitch angle and azimuth angle over time

场景二：将轨道机动段划分为10 个子时间区间，每个子时间区间设置30个离散点。俯仰角和方位角的设置与场景一相同。时刻的初值设为10，其上界和下界分别设为20 和0。仿真结果如图8～图11所示。将场景二与场景一进行对比可知，航天器在场景二中没有轨道等待段，在初始时刻即进行轨道机动，tf时刻到达目标位置，实现了预定轨道机动任务。场景二的总机动时间为6.933 tU，机动时间比场景一多1.144 4 tU，任务总飞行时间比场景一少1.175 tU。到达目标位置时场景二中航天器质量为0.888 1 mU，场景二的燃料消耗比场景一多0.018 5 mU。场景二仿真结果表明，在完成相同的轨道机动任务情况下，本文所提轨道机动策略的燃料消耗更低。

图8 三维空间的位置变化Fig.8 Variations of position in three-dimensional space

图9 速度随时间的变化Fig.9 Variation of velocity over time

图10 质量随时间的变化Fig.10 Variation of mass over time

图11 俯仰角和方位角随时间的变化Fig.11 Variations of pitch angle and azimuth angle over time

为了进一步分析航天器轨道机动过程中的燃料消耗，将本文方法与非连续推力的脉冲变轨方式进行对比。以脉冲施加时刻和相应的脉冲速度增量为待优化量，以总脉冲速度增量作为优化目标。各方向脉冲速度增量的上界和下界分别为0.05 vU 和-0.05 vU。根据场景二总燃料消耗折算成的总速度增量可知，最多施加10次脉冲便可完成预定的轨道机动任务，因此本文采用10次脉冲方式进行轨道机动。相邻脉冲施加时刻之间的间隔上界为7.884 2，下界为0。本文中，各方向脉冲速度增量的初值均设为0.01 vU，相邻脉冲施加时刻之间的间隔初值均设为1。轨道机动的始末端轨道参数和各物理量与连续推力情况相同，非连续推力的脉冲变轨仿真结果如图12 所示。航天器的脉冲施加时刻和各方向速度增量如表3所示。

表3 脉冲时刻和速度增量Table 3 Pulse moments and velocity increments

图12 三维空间的位置变化Fig.12 Variations of position in three-dimensional space

由图12 可知，经过8 次脉冲变轨，航天器实现了预定的轨道机动任务，航天器总飞行时间为10.551 2 tU。将图12 与图4、图8 进行对比可知，与连续推力方式的轨道机动不同，脉冲变轨方式先抬升轨道高度，在远地点附近施加脉冲以改变轨道面。在完成轨道机动任务的前提下实现燃料优化目标。整个轨道机动任务过程中，航天器的总脉冲速度增量为0.491 3 vU，经换算，到达目标位置时航天器质量为0.923 8 mU。非连续推力的脉冲变轨方式燃料消耗比场景一低0.017 2 mU，这是因为连续推力变轨过程中存在无效推力而造成额外质量损失，脉冲变轨方式是理想化的变轨方式。

为了验证本文方法的在轨道机动任务中的普适性，有必要对航天器在平面内的轨道机动进行仿真分析。平面内轨道机动的始末端轨道根数如表4所示，表5 为无量纲化处理后的柱坐标系描述的始末端状态。各仿真参数设置与场景一相同，仿真结果如图13～图16所示。

表4 机动轨道始末端轨道根数Table 4 Orbital elements at the beginning and end of the maneuvering orbit

表5 无量纲化始末端状态Table 5 The dimensionless initial states and end states

图13 三维空间中的位置变化Fig.13 Variations of position in three-dimensional space

图14 速度随时间的变化Fig.14 Variation of velocity over time

图15 质量随时间的变化Fig.15 Variation of mass over time

图16 俯仰角和方位角随时间的变化Fig.16 Variations of pitch angle and azimuth angle over time

由图13可知，航天器在初始轨道先飞行一段时间进行轨道等待，在t1时刻开始进行轨道机动，t2时刻到达目标位置，实现了平面内的轨道机动任务。航天器无动力飞行时间为3.16 tU，总机动时间为5.917 tU，在t2时刻航天器质量为0.906 6 mU。

5 结论

本文研究了常值推力作用的航天器轨道机动优化设计问题。建立了无量纲化的轨道机动动力学模型并提出了先轨道等待再轨道机动的优化设计策略。通过多重直接打靶法将复杂的轨迹优化问题转化为非线性规划问题并采用内点法进行求解。仿真结果表明，本文所提方法在轨道机动任务中的具有普适性，不仅适用于平面轨道机动任务，同样适用于空间轨道机动任务。与传统方法相比，本文所提的轨道机动优化设计策略燃料消耗更低。