基于粒子群优化的月球捕获制动方法研究

郭继峰，白成超，贺国平

（哈尔滨工业大学航天工程系，哈尔滨150001）

1 引言

月球探测是航天技术发展的主要领域之一。目前，我国已经成功实施了月球探测一期工程（绕）和二期工程（落），三期工程（回）也已立项实施，而轨道设计是探月任务的关键技术。月球轨道分为停泊、地月转移、捕获制动、环绕和着陆五类，其中地月转移轨道相对于月球为双曲线轨道，在探月卫星飞至近月点时必须进行减速机动以被月球引力捕获，否则卫星将飞离地月系统且不会返回［1］。因此，捕获轨道是轨道设计中的核心之一，直接关系到轨道器任务实施的成败。我国嫦娥一号和嫦娥二号探月卫星采用近月点三次减速机动的方案实现月球捕获，最终进入平均高度分别为200 km和100 km的极月圆形目标轨道［2⁃3］。 嫦娥三号探测器变推力发动机经过361 s制动后，顺利进入距月面平均高度约100 km的环月轨道［3］。

探月卫星轨道的捕获制动策略可分为固定推力方向制动和变推力方向制动，它对捕获时燃料消耗、捕获后目标轨道确定等任务规划和探月总体设计有直接的影响［4］。另外，由于探测器可携带燃料有限，对轨道捕获策略进行优化设计，使完成捕获任务的燃料消耗最小也十分重要，常用的优化算法主要包括粒子群算法、蚁群算法、模拟退火算法、遗传算法、差分进化算法等。并且，Myatt证明了相比于遗传算法和粒子群算法，差分进化算法在深空引力辅助机动轨道设计中的优势［5］；而Bessette比较了粒子群算法和差分进化算法在轨道设计问题中的适用性［6］；Vasile进一步提出了一种结合差分进化算法和单纯形算法中启发式的混合算法，并仿真验证了该混合算法在深空轨道设计中的优越性［7］。

目前，国内很多学者也致力于捕获制动段探月轨道设计与优化策略，并且取得了一定得成果。蒋小勇提出用动力学分解方法来快速获得小推力轨道的近似解析解［8］。黄文德以载人登月轨道和中止轨道一体化设计为指导，提出了以快速设计结果为初值的载人登月摄动轨道设计方法，同时针对设计参数多、对初值比较敏感等特点，采取参数群体优化的方法进行轨道设计［9］。李军锋研究了捕获制动结束后卫星绕月球轨道半长轴和偏心率给定情况下的燃料最优问题，通过求解最优控制对应的两点边值问题，得到了燃料最优情况下的初始条件、捕获策略和燃料消耗，验证了最优捕获策略为捕获发动机开机后以最大幅值工作至捕获结束，推力的方向近似沿着速度相反方向［4］；他还与其合作者基于有限推力模型，以绕月轨道的半长轴和偏心率为目标约束，分析了燃料最优捕获、姿态匀速转动和姿态惯性定向捕获策略，并分别利用最优控制和粒子群优化算法针对不同捕获策略进行了求解［10］。王云财等对双曲线轨道的近心点高度进行了优化处理，并在此基础上对整个捕获过程中的推力方向和发动机的开关机时间进行优化，最终得到整个捕获过程中的控制参数，进一步降低了捕获过程中的燃料消耗［11］。

对于月球捕获制动段轨道设计优化问题，以上提到的方法存在求解分析过程复杂的不足，本文在考虑地球摄动的前提下，基于固定推力方向制动策略进行月球捕获制动研究，以燃料消耗以及最终飞行器近月点高度偏差为指标，采用简单有效的粒子群优化算法，进行捕获制动轨迹优化。

2 月球捕获制动问题描述

2.1 动力学模型

在月球J2000惯性坐标系下，考虑月球引力和地球的摄动力，探测器轨道动力学方程如式（1）［3］：

其中，m为探测器的质量；Isp为推力器比冲；g为地球重力加速度；r和v分别为探测器在月球J2000惯性坐标系下的位置矢量和速度矢量，且r＝‖r‖；rE为地球在月球J2000惯性坐标系下的位置矢量，rE＝‖rE‖，考虑地球对探测器的摄动时，不考虑地球扁率的影响；rpE为探测器相对地球的位置矢量，且rpE＝ r－ rE，μL、μE分别为月球和地球的引力常数；F为捕获制动时推力器产生的推力。忽略太阳光压微小摄动力。

2.2 固定推力方向制动策略

固定推力方向制动策略工程上实现简单，是一种很可取的捕获制动方式。定义探测器近月球点的位置和速度单位矢量分别为＾rp和＾vp，常推力方向制动中，令推力矢量保持在＾rp和＾vp确定的平面内，推力方向角定义为α，从＾vp反方向向上测量为正，从＾rp正方向向下测量为负，如图1所示，则存在推力方向平衡关系如式（2）：

开始点火的位置矢量与预测近月点的位置矢量的夹角为点火角，如图2所示。

图1 推力方向及推力方向角示意图Fig.1 Schematic diagram of thrust direction and thrust angle

图2 捕获制动角度关系与点火角定义Fig.2 Angle relations of braking process and igni⁃tion angle

假设开始进行连续推力捕获制动时，飞行器的轨道倾角满足指标要求，仅需进行平面内制动。由此可以确定需要计算的参数包括：初始点火角θ、推力方向角α以及点火时长Tp。

3 粒子群优化

3.1 粒子群优化算法

在PSO算法中，优化问题在所有搜索范围内的解向量都对应着多维度的有限空间里的一个点，我们称之为“位置”，而有一定数量的占据随机位置的初始解向量，被称为“粒子”，每一个粒子都占据一个位置，每一个的粒子所在的位置都对应一个适应值，称为“适应度”，每一个粒子还都有一个决定下一次迭代时位置移动量的值，称为“速度”，在迭代过程中，粒子们就根据上一次迭代时的速度、自身经历过的适应度最高的位置和所有粒子经历过的适应度最高的位置决定的速度对位置进行更新［12］。

在一个D维的搜索空间中，由n个粒子组成的种群X可以表示为式（3）：

其中第i个粒子为一个如式（4）所示的D维向量：

该空间里，第i个粒子的速度可以表示为式（5）：

在迭代搜索过程中，定义如下两个目标位置作为迭代的参考位置：

1）单个粒子自身在迭代过程经历过的适应度最高的位置pi（局部最优），如式（6）所示：

2）全部粒子在迭代过程中经历过的适应度最高的位置pg（全局最优），如式（7）所示，pg在每次迭代中都是唯一的：

所示，速度和位置更新公式就如式（8）～（9）所示［12］：

其中 i ＝ 1，2，…，n，d ＝ 1，2，…，D；k 为当前迭代次数；ω为惯性权重；c1，c2为加速常数；ξ，η为 0，1[ ]区间的随机数为了防止粒子的盲目搜索，一般建议将位置或者速度限制在一定的搜索区间。PSO算法流程如图3所示。

图3 PSO算法流程Fig.3 Flow chart of PSO algorithm

如前所述，对于捕获制动轨道优化问题，需要计算的参数包括初始点火角θ、推力方向角α以及点火时长Tp三个。由于制动过程中轨道偏心率单调下降，所以优化过程中，将轨道偏心率的值作为制动过程结束的判断条件，故点火时长不作为粒子群优化的变量，优化变量只有点火角和推力方向角，而约束只剩下近月点高度一个。因子此处D ＝ 2，Xi＝ [ Xi1Xi2]，Xi1表示推力方向角，Xi2表示点火角；vi＝ [ vi1vi2]，对应为两个优化变量的变化速度。

3.2 惯性权重的选择

惯性权重ω体现的是粒子继承先前的速度的能力，一个较大的惯性权值有利于全局搜索，而一个较小的惯性权值更有利于局部搜索。为了更好的平衡算法的全局搜索与局部搜索能力，采用式（10）所示惯性权重选择方法［12］：

其中ωstart为初始惯性权重，ωend为迭代至最大次数时的惯性权重，k为当前迭代次数，Tmax为最大迭代次数。一般来说，惯性权值ωstart＝0.9，ωend＝0.4左右算法性能最好。这样，随着迭代的进行，惯性权重由0.9线性递减至0.4，迭代初期较大的惯性权重使算法保持了较强的全局搜索能力，而迭代后期较小的惯性权重有利于算法进行更精确的局部搜索。

4 基于粒子群优化的固定推力制动方法验证

4.1 参数设置

在Visual Studio 2015下，编写C＋＋仿真主程序，进行仿真。根据月球半径及其影响球半径，假设给定以下仿真参数［13］：

末端近月点距离约束：rpf＝1938 km；

最终绕月椭圆轨道离心率约束：e＝0.6845 ；

优化指标：J＝min( Tp)。

固定推力制动策略中，主要参数有点火角，推力方向角与点火时间（制动时间）三个，因此仿真设置了3个优化变量，2个等式约束，1个优化指标。由于制动过程中轨道偏心率单调下降，所以优化过程中，将轨道偏心率的值作为制动过程结束的判断条件，故点火时长不作为粒子群优化的变量，所以优化变量只有点火角和推力方向角，而约束只剩下近月点高度一个。使用PSO寻优算法优化，通过罚函数将等式约束去掉。粒子群优化基本参数设定如表1所示。

表1 参数搜索区间Table 1 Search interval of parameters

从而得到新的粒子群优化指标如式（11）：

假定初始双曲线轨道参数如表2所示，在理想姿态无误差情况下，探测器飞入月球主引力区并制动，制动前探测器质量为120 t，主发动机推力大小为200 kN，比冲为900 s，优化轨道制动参数，考察制动结果和过程。

表2 初始双曲线轨道近月点参数Table 2 Initial hyperbolic orbit near point parameters

计算地球引力时，使用JPL DE421星历表查询地月间的相对运动参数，并假定起始时刻为2017年4月21日0时0分0秒。

4.2 仿真结果

优化得到的制动变轨三要素如表3所示：

表3 固定推力方向制动策略参数Table 3 Parameters of braking strategy in fixed thrust direction

制动后轨道六根数如表4所示：

表4 制动后轨道六根数Table 4 Six orbital elements after braking

制动后指标变量如表5所示。

表5 制动结束后指标变量Table 5 Indicator variables after braking

近月点处的距离、速度、离心率、轨道倾角的设定值、真实值及设定值和真实值之间的偏差如表6所示。

表6 多种参数的偏差Table 6 Deviation of parameters

由表6可看出，制动结束后，近月点高度误差为0.000 27 km，小于100 m；近月点速度偏差为 0.000 513 km／s，小于 1 m／s；离心率的偏差为0.001 058，误差百分率约为1%；轨道倾角偏差为0.000 11°，小于 0.05°。 以上结果说明，固定推力方向制动策略可以通过优化，在消耗推进剂尽可能少的同时满足指定的终端状态约束。

从飞行器进入月球影响球开始到捕获制动结束后飞行器绕月球飞行近一圈为止飞行器的运动轨迹仿真图如图4所示。其中红色曲线为制动前的轨迹，蓝色曲线为制动过程中的轨迹，绿色曲线为捕获制动结束后的绕月球飞行的椭圆轨迹。

图4 从进入月球影响球开始的整体捕获制动轨迹Fig.4 Capture trajectory after entering the influence ball of the moon

图5 为捕获制动过程轨迹的放大图，同样，红色曲线为制动前的轨迹，蓝色曲线为制动过程中的轨迹，绿色曲线为捕获制动结束后的绕月球飞行的椭圆轨迹的一部分。

制导过程中以离心率减小到一定值为终止制动条件，离心率随时间的变化如图6所示。

图5 捕获制动过程轨迹放大图Fig.5 Enlarged view of capture process trajectory

图6 捕获制动过程中离心率随时间的变化Fig.6 Change of eccentricity over time during cap⁃ture process

优化指标是使燃耗最少，因此，捕获制动过程中，对飞行器质量变化的追踪也至关重要。飞行器质量随时间变化的过程如图7所示。

图7 飞行器质量随时间变化Fig.7 Mass change of the spacecraft over time

捕获制动过程中，轨道能量和轨道偏心率是单调下降的，由于是平面内制动，轨道倾角只由摄动力影响，所以轨道倾角变化很小，由轨道近月点高度在固定推力方向制动完成时正好达到200 km的目标近月点高度。

5 结论

本文在考虑地球摄动的情况下，基于固定推力方向制动策略，以燃料消耗以及最终飞行器近月点高度偏差为指标，以最终绕月轨道偏心率以及最终绕月轨道近心点高度为约束，采用粒子群算法进行捕获制动段的轨道优化。仿真结果表明，本文方法可有效完成月球对探测器的捕获，探月卫星初始以双曲线轨道进入月球影响球，通过捕获制动，最终到达绕月大椭圆轨道上。最终状态约束误差表明，本方法能在满足所有约束的情况下最小化燃耗，这验证了基于粒子群优化的固定推力方向制动方法的有效性。此外，本文没有复杂的方程分析，这说明了本文方法的简单与易实现性。

目前在考虑地球对探测器的摄动影响时，没有考虑地球的扁率，认为地球是一个规则的球体，此外没有考虑太阳及光压等摄动项的影响，在今后的工作中，将进一步考虑这些摄动影响，建立更为精确的动力学模型。同时，将对多种捕获制动策略展开研究，并进行对比分析。

［1］ 李革非，韩潮.月球捕获控制分析［C］／／全国第十二届空间及运动体控制技术学术会议，2006：240⁃244.

Li Gefei， Han Chao.Analysis on Lunar capture control［C］ ／／Twelfth National Conference on Space and Movement Control Technology， 2006： 240⁃244.（in Chinese）

［2］ 王宏，董光亮，胡小工，等.月球探测器捕获分析及评估［J］.测绘科学技术学报， 2007，24（6）： 406⁃409.

Wang Hong， Dong Guangliang， Hu Xiaogong， et al.Analysis and assessment of the Lunar probe captured by Moon［J］.Journal of Zhengzhou Institute of Surveying and Mapping，2007， 24（6）： 406⁃409.（in Chinese）

［3］ 盛伟强.火星近火点捕获制动姿轨一体化控制［D］.哈尔滨：哈尔滨工业大学，2015.

Sheng Weiqiang.Study on Coupling Control of Attitude and Orbit for Mars Orbit Insertion at Periareon［D］.Harbin： Har⁃bin Institute of Technology， 2015.（in Chinese）

［4］ 李军锋.基于精确动力学模型的火星探测轨道设计研究［D］.北京：清华大学，2012.

Li Junfeng.Research on Orbit Design for Mars Explorer by the Accurate Dynamic［D］.Beijing： Tsinghua University，2012.（in Chinese）

［5］ Myatt D R，Becerra V M，Nasuto SJ，et al.Advanced global optimization tools for mission analysis and design［R］.Final Report of ESA Ariadna ITT AO4532 ／18138 ／04 ／NL／MV，2004.

［6］ Bessette C R，Spencer D B.Performance comparison of sto⁃chastic search algorithms on the interplanetary gravity⁃assist trajectory problem［J］.Journal of Spacecraft and Rockets，2015， 44（3）： 722⁃724.

［7］ Vasile M，Minisci E.Analysis of some global optimization al⁃gorithms for space trajectory design［J］.Journal of Spacecraft and Rockets， 2010， 47（2）： 334⁃344.

［8］ 蒋小勇.深空探测器小推力轨道优化设计方法及其应用研究［D］.长沙：国防科学技术大学，2014.

Jiang Xiaoyong.Study on Low⁃Thrust Trajectory Design and Optimization Method for Deep Space Explorer and Application［D］.Changsha： National University of Defense Technology，2014.（in Chinese）

［9］ 黄文德.载人登月中止轨道的特性分析与优化设计［D］.长沙：国防科学技术大学，2011.

Huang Wende.Characteristic Analysis and Optimal Design of Abort Trajectories for Manned Lunar Landing Mission［D］.Changsha： National University of Defense Technology， 2011.（in Chinese）

［10］ 李军锋，龚胜平.有限推力模型火星探测捕获策略分析［J］.中国科学： 物理学 力学 天文学，2013（06）：781⁃786.

Li Junfeng， Gong Shengping.Analysis of capture strategies for Mars explorer with finite⁃thrust［J］.Sci Sin⁃Phys Mech As⁃tron， 2013（06）： 781⁃786.（in Chinese）

［11］ 王云财，刘辉，胡晓赛.一种改进的火星探测捕获段轨道设计方法［C］／／全国航天动力学与控制学术会议.2015.

Wang Yuncai， Liu Hui， Hu Xiaosai.An improved Martian exploration capture trajectory design method［C］／／National Conference on Space Dynamics and Control， 2015.（in Chi⁃nese）

［12］ 史峰，王辉，郁磊，等.MATLAB智能算法30个案例分析［M］.北京：北京航空航天大学出版社，2011：102⁃107.

Shi Feng， Wang Hui， Yu Lei， et al.30 Cases Analysis of MATLAB Intelligent Algorithm［M］.Beijing： Beihang Uni⁃versity Press， 2011： 102⁃107.（in Chinese）

［13］ 郗晓宁，曾国强，任萱，等.月球探测器轨道设计［M］.北京：国防工业出版社，2001：240⁃248.

Xi Xiaoning， Zeng Guoqiang， Ren Xuan， et al.Lunar Explo⁃ration Orbit Design［M］.Beijing： National Defence Industry Press， 2001： 240⁃248.（in Chinese）