-
飞跃航天器可以追溯到美国国家航空航天局(National Aeronautics and Space Administration,NASA)发射的“勘测者6号”(Surveyor VI),1967年11月17日,该探测器着陆月球后,又重新点火微调发动机2.5 s,向西飞跃了2.4 m远[1]。2007年,谷歌发起了月球X大奖赛,要求用私人资金发送一个探测器着陆到月球,在月球表面移动至少500 m,向地球传递录像、图像、数据。参赛的Next Giant Leap、Moon Express、SpacelL团队都提出了各自的飞跃探测器。其中Next Giant Leap团队的飞跃探测器由Draper实验室开发,其导航、制导与控制(Guidance,Navigation and Control,GNC)技术通过了一系列地面验证。目前飞跃探测器相关GNC技术尚未成熟,除Surveyor VI之外尚未发现应用案例,常规着陆器的GNC技术[2-3]对飞跃探测器具有一定的借鉴意义。
任务环境存在诸多不确定性、复杂性、多变性,飞跃过程受到地形约束、燃料约束、姿态约束、角速度约束,轨迹末端需考虑避障要求。特殊的任务需求与严苛的约束条件对月面飞跃探测器轨迹的适应范围以及自主规划能力提出了更高的要求。
针对上述问题,本文根据飞跃期间任务约束种类的不同,把飞跃轨迹划分为垂直上升段、程序转弯段、无动力滑行段、接近段、缓速下降段5个阶段,采用粒子群优化算法,得到了一种月面弹道式飞跃探测的轨迹生成方法。
-
在射向平面内建立二维平动方程
$$ \dot { X} = f\left( { {X},{ U}} \right) $$ (1) 其中
$$ { X} = \left[ {\begin{split} x\; \\ y \;\\ {{v_x}} \\ {{v_y}} \\ m \end{split}} \right],\;\;{ U} = \left[ {\begin{split} F \\ \theta \end{split}} \right] $$ $$ f\left( {{ X},{ U}} \right) = \left[ {\begin{split} & \quad \quad {{v_x}} \\ & \quad \quad {{v_y}} \\ & {F\cos \theta /m - {g_{\rm{m}}}} \\ & \quad {F\sin \theta /m} \\ & \quad { - F/{I_{{\rm{sp}}}}} \end{split}} \right] $$ 其中:
$x$ 、$y$ 分别为高程与航程;${v_x}$ 、${v_y}$ 分别为天向速度与航向速度大小;俯仰角$\theta $ 定义为探测器纵轴与天向夹角(头朝前为正);$m$ 为质量;F 为主发动机推力(沿纵轴指向头部);${I_{{\rm{sp}}}}$ 为主发动机比冲;${g_{\rm{m}}}$ 为月球引力加速度常数;${v_x}$ 、${v_y}$ 、$\theta $ 、F 初值均为0。考虑起飞安全、飞跃过程安全、俯仰角速度限幅、燃料消耗与主发动机推力范围,飞跃过程中需要满足的约束为
$$ \left\{ \;\; \begin{aligned} & {t_1} > {t_{\min }} \\ & x > H\left( y \right) \\ & \left| {{\omega _z}} \right| < {\omega _{\max }} \\ & m > {m_{\rm{d}}} \\ & 0 < F < {F_{\max }} \end{aligned} \right.$$ (2) 其中:
${t_1}$ 为设计参数,表示垂直上升段的结束时刻;${t_{\min }}$ 为垂直上升的最短时间;$H\left( y \right)$ 为航程$y$ 处的地形高度;${\omega _z}$ 为俯仰角速度;${\omega _{\max }}$ 为俯仰角速度的最大值;${m_{\rm{d}}}$ 为探测器干重;${F_{\max }}$ 为主发动机最大推力。 -
参考现有的着陆上升技术,本文把首次飞跃轨迹划分为垂直上升段、程序转弯段、无动力滑行段、接近段、缓速下降段5个阶段,其中垂直上升段主发动机最大推力工作,保证最短时间达到安全高度;程序转弯段主发动机仍然最大推力工作,同时按照恒定的俯仰角速度进行转弯;当探测器达到一定的飞行高度与飞行速度后,主发动机关机,进入无动力滑行段以节约燃料;无动力滑行段探测器一直在预测减速制动所需要的推力,当预测推力达到设定值后,进入接近段;接近段制导引入闭环,主发动机根据制导输出提供推力,使探测器到达安全着陆点上方;之后进入缓速下降段,探测器开始匀速垂直下降。
-
根据2.1节描述,转段条件为状态量、时间、当前飞行阶段以及一些设计参数的函数
$$ {l_{k + 1}} = {f_l}\left( {{ X},t,{l_k},{Z}} \right) $$ (3) 其中:
$l = 1, \cdots ,6$ 为飞行阶段序号,1~5分别表示5个飞行阶段,6表示结束;${l_k}$ 为当前控制周期的阶段序号;${l_{k + 1}}$ 为下一控制周期的阶段序号;$t$ 为从起飞到当前时刻的飞行时间;$Z$ 表示待设计参数${t_1}$ 等。式(3)在各段的具体形式为1)垂直上升段至程序转弯段
$$ {l_{k + 1}} = \left\{ \begin{split} 2,\;{l_k} = 1{\text{且}}t \geqslant {t_1}\\ 1,\;{l_k} = 1{\text{且}}t < {t_1} \end{split} \right. $$ (4) 2)程序转弯段至无动力滑行段
$$ {l_{k + 1}} = \left\{ \begin{split} 3,\;{l_k} = 2{\text{且}}\theta \geqslant {\theta _2}\\ 2,\;{l_k} = 2{\text{且}}\theta < {\theta _2} \end{split} \right. $$ (5) 3)无动力滑行段至接近段
$$ {l_{k + 1}} = \left\{ \begin{split} 4,\;{l_k} = 3{\text{且}}{F_{\rm{p}}} \geqslant {F_3}\\ 3,\;{l_k} = 3{\text{且}}{F_{\rm{p}}} < {F_3} \end{split} \right. $$ (6) 4)接近段至缓速下降段
$$ {l_{k + 1}} = \left\{ \begin{split} 5,\;{l_k} = 4{\text{且}}x - H\left( y \right) \leqslant {h_4}\\ 4,\;{l_k} = 4{\text{且}}x - H\left( y \right) > {h_4} \end{split} \right. $$ (7) 5)任意段至结束
$$ {l_{k + 1}} = \left\{ \begin{split} 6,\;x \leqslant H\left( y \right){\text{或}}m \leqslant {m_{\rm{d}}}\\ {l_k},\;x > H\left( y \right){\text{且}}m > {m_{\rm{d}}} \end{split} \right. $$ (8) 其中:
${\theta _2}$ 、${F_3}$ 、${h_4}$ 为设计参数,分别表示程序转弯段结束时的俯仰角、无动力滑行段结束时的预测推力大小、接近段结束时的相对高度;${F_{\rm{p}}}$ 为预测推力大小。 -
根据2.1节描述,控制量为状态量、时间、当前飞行阶段以及一些设计参数的函数
$$ U = {f_u}\left( {{ X},t,l,{ Z}} \right) $$ (9) 式(9)在各段的具体形式为
1)垂直上升段
$$ F = {F_{\max }},\;\;\;\theta = 0 $$ (10) 2)程序转弯段
$$ F = {F_{\max }},\;\;\;\theta = {\omega _2}\left( {t - {t_1}} \right) $$ (11) 其中:
${\omega _2}$ 为设计参数,表示程序转弯段俯仰角速度。3)无动力滑行段
$$ F = 0,\;\;\;\theta = {\arctan}\left( {{a_y}/{a_x}} \right) $$ (12) 令
${{a}} = {\left[ {\begin{array}{*{20}{c}} {{a_x}}&{{a_y}} \end{array}} \right]^{\rm{T}}} = {f_a}\left( X \right)$ 为接近段制导律输出的加速度。4)接近段
$$ F = m\left\| {{a}} \right\|,\;\;\;\theta = {\arctan}\left( {{a_y}/{a_x}} \right) $$ (13) 5)缓速下降段
$$ F = m{g_{\rm{m}}},\;\;\;\theta = 0 $$ (14) -
基于式(1)平动方程、式(4)~(8)转段条件、式(10)~(14)控制量,可以采用Runge-Kutta方法数值解算得到状态量的终端值。为了减少运算时间,便于在线应用,本文推导得到了一种快速运算方法。
下一控制周期的状态量为当前状态量、时间、当前飞行阶段以及一些设计参数的函数
$$ {{ X}_{k + 1}} = \varPhi \left( {{{ X}_k},{{ U}_k},{l_k},{ Z}} \right) $$ (15) 式(15)在各段的具体形式为
1)垂直上升段
该段可以解析得到转段时各状态量的表达式
$$\left\{ \begin{aligned} {x_1} & = {x_0} - {m_0}I_{_{{\rm{sp}}}}^2\ln \left( {{m_0}{I_{{\rm{sp}}}}} \right)/{F_{\max }}+ \\ & {I_{{\rm{sp}}}}{t_1}\left( {1 + \ln \left( {{m_0}{I_{{\rm{sp}}}}} \right)} \right) - 0.5{g_{\rm{m}}}t_{_1}^2 - \\ & {I_{{\rm{sp}}}}{t_1}\ln \left( {{m_0}{I_{{\rm{sp}}}} - {F_{\max }}{t_1}} \right)+ \\ & {m_0}I_{_{{\rm{sp}}}}^2\ln \left( {{m_0}{I_{{\rm{sp}}}} - {F_{\max }}{t_1}} \right)/{F_{\max }} \\ & {y_1} = {y_0} \\ & {v_{x1}} = {I_{{\rm{sp}}}}\ln \left( {{m_0}{I_{{\rm{sp}}}}} \right) - {g_{\rm{m}}}{t_1} - \\ & {I_{{\rm{sp}}}}\ln \left( {{m_0}{I_{{\rm{sp}}}} - {F_{\max }}{t_1}} \right) \\ & {v_{y1}} = 0 \\ & {m_1} = {m_0} - {F_{\max }}{t_1}/{I_{{\rm{sp}}}} \\ & {\theta _1} = 0 \\ \end{aligned} \right. $$ (16) 其中:
${x_0}$ 、${y_0}$ 、${m_0}$ 分别为初始高程、航程、质量;${x_1}$ 、${y_1}$ 、${v_{x1}}$ 、${v_{y1}}$ 、${m_1}$ 、${\theta _1}$ 分别为垂直上升段结束时的高程、航程、天向速度、航向速度、质量、俯仰角。2)程序转弯段
可以解析得到该段的质量表达式
$$ m = {m_1} - {F_{\max }}\left( {t - {t_1}} \right)/{I_{{\rm{sp}}}} $$ (17) 与俯仰角表达式
$$ \theta = {\omega _2}\left( {t - {t_1}} \right) $$ (18) 该段的终端速度与位置不能解析得到,本文通过数值积分得到。
垂向速度表达式为
$$ \begin{aligned}[b] & {v_x} = {v_{x1}} - {g_{\rm{m}}}\left( {t - {t_1}} \right)+ \\ & \;\;\;\;\; \int_0^{t - {t_1}} {{F_{\max }}\cos \left( {\omega \tau } \right)/\left( {{m_1} - {F_{\max }}\tau /{I_{{\rm{sp}}}}} \right)} {\rm{d}}\tau \\ \end{aligned} $$ (19) 航向速度表达式为
$${v_y} = {v_{y1}} + \int_0^{t - {t_1}} {{F_{\max }}\sin \left( {\omega \tau } \right)/\left( {{m_1} - {F_{\max }}\tau /{I_{{\rm{sp}}}}} \right)} {\rm{d}}\tau $$ (20) 高程表达式
$$ x = {x_1} + \int_0^{t - {t_1}} {{v_x}\left( {{t_1} + \tau } \right)} {\rm{d}}\tau $$ (21) 航程表达式
$$ y = {y_1} + \int_0^{t - {t_1}} {{v_y}\left( {{t_1} + \tau } \right)} {\rm{d}}\tau $$ (22) 该段结束时间
${t_2} = {\theta _2}/{\omega _2}$ ,对式(19)和式(20)分别做数值积分可以得到终端垂向速度与航向速度;对式(21)和式(22)分别做二次数值积分可以得到终端高程与航程。3)无动力滑行段
该段需要一直预测当前需要的主发动机推力,因此需要解算每一控制周期的状态。
主发动机关机,质量不变,故
$$ m = {m_2} $$ (23) 垂向速度表达式为
$$ {v_x} = {v_{x2}} - {g_{\rm{m}}}\left( {t - {t_2}} \right) $$ (24) 航向速度表达式为
$$ {v_y} = {v_{y2}} $$ (25) 高程表达式
$$ x = {x_2} + {v_{x2}}\left( {t - {t_2}} \right) - 0.5{g_{\rm{m}}}{\left( {t - {t_2}} \right)^2} $$ (26) 航程表达式
$$ y = {y_2} + {v_{y2}}\left( {t - {t_2}} \right) $$ (27) 其中:
${x_2}$ 、${y_2}$ 、${v_{x2}}$ 、${v_{y2}}$ 、${m_2}$ 分别为程序转弯段结束时的高程、航程、天向速度、航向速度、质量。基于当前位置、速度、质量以及目标位置、速度、加速度,可以预测当前状态对应的接近段制导律输出加速度的大小与方向,该方向就是本体纵轴方向,由此确定了无动力滑行段的俯仰姿态。
4)接近段
该段也需要一直预测当前需要的主发动机推力,因此需要解算每一控制周期的状态。
质量表达式为
$$ m = {m_3} - \int_{{t_3}}^t F /{I_{{\rm{sp}}}}{\rm{d}}\tau $$ (28) 垂向速度表达式为
$$ {v_x} = {v_{x3}} - {g_{\rm{m}}}\left( {t - {t_3}} \right) + \int_0^{t - {t_1}} {F\cos \theta /m} {\rm{d}}\tau $$ (29) 航向速度表达式为
$$ {v_y} = {v_{y3}} + \int_{{t_3}}^t {F\sin \theta /m} {\rm{d}}\tau $$ (30) 高程表达式
$$ x = {x_3} + \int_{{t_3}}^t {{v_x}} {\rm{d}}\tau $$ (31) 航程表达式
$$ y = {y_3} + \int_{{t_3}}^t {{v_y}} {\rm{d}}\tau $$ (32) 其中:
${x_3}$ 、${y_3}$ 、${v_{x3}}$ 、${v_{y3}}$ 、${m_3}$ 分别为无动力滑行段结束时的高程、航程、天向速度、航向速度、质量。主发动机推力大小
$F$ 与俯仰角$\theta $ 由接近段制导律以及估计质量得到。接近段每个控制周期都进行制导解算,同时对式(28)~(32)进行一步数值积分。
5)缓速下降段
该段可以解析得到触月时各状态量的表达式
$$ \left\{ \; \begin{aligned} & {x_5} = {x_{\rm{f}}} \\ & {y_5} = {y_4} \\ & {v_{x5}} = {v_{x4}} \\ & {v_{y5}} = {v_{y4}} \\ & {m_5} = {m_4}{{\rm{e}}^{\frac{{{g_{_{\rm m}}}{h_4}}}{{I_{_{s{\rm p}}}{v_{{x_4}}}}}}} \\ & {\theta _5} = 0 \end{aligned} \right. $$ (33) 其中:
${y_4}$ 、${v_{x4}}$ 、${v_{y4}}$ 、${m_4}$ 分别为接近段结束时的航程、天向速度、航向速度、质量;${x_5}$ 、${y_5}$ 、${v_{x5}}$ 、${v_{y5}}$ 、${m_5}$ 、${\theta _5}$ 分别为触月时的高程、航程、天向速度、航向速度、质量、俯仰角;${x_{\rm{f}}}$ 为着陆点地形高度,由于缓速下降段水平速度很小,因此${x_{\rm{f}}} \approx H\left( {{y_4}} \right)$ 。 -
参考剖面待设计参数总共有5个:
${t_1}$ 、${\omega _2}$ 、${\theta _2}$ 、${F_3}$ 、${h_4}$ 。即当这5个参数确定后,把式(2)约束条件、式(4)~(8)转段条件、式(10)~(14)控制量代入平动方程(1),可以唯一确定一条飞行轨迹。当${F_3}$ 接近${F_{\max }}$ 时,无动力滑行时间延长,有利于节约推进剂;当${h_4}$ 接近0时,匀速下降时间减小,同样有利于节约推进剂。因此${F_3}$ 和${h_4}$ 直接给定,不再参与优化搜索。优化参数选取为$$ { Z} = {\left[ {\begin{array}{*{20}{c}} {{t_1}}&{{\omega _2}}&{{\theta _2}} \end{array}} \right]^{\rm{T}}} $$ (34) 考虑飞行安全、轨迹高程差、轨迹航程差、探测器姿控性能、探测器推重比等因素,优化参数范围设定为
$$ {{ Z}_d} \underline \prec { Z} \underline \prec {{ Z}_u} $$ (35) 其中:
${{ Z}_d}$ 、${{ Z}_u}$ 分别为优化参数的下界与上界;$\underline \prec $ 代表${ Z} $ 上的向量不等式或分量不等式。 -
优化的主要目标为寻找一条飞行轨迹,保证航程误差在允许的范围内燃料消耗最少。由于优化过程引入了制导律,优化目标中不再考虑航程误差。为了保证避障敏感器工作条件以及着陆速度、姿态满足指标要求,最终选定的飞行轨迹必须包含5个飞行阶段。另外,由于地形高度的不确定性,为了保证飞行安全,本文在优化目标中引入接近段相对月面的最小高度
${h_{\min }}$ 。优化目标为$$ J\left( { Z} \right) = \left\{ \begin{split} & \Delta m + {k_1}{{\rm e}^{ - {k_2}{h_{\min }}}},\;n = 5 \\ & N,\quad \quad \quad \quad\quad n < 5 \\ \end{split} \right. $$ (36) 其中:
$\Delta m = {m_0} - {m_5}$ 为燃料消耗量;${k_1}$ 、${k_2}$ 为设计参数;$n$ 为轨迹包含的飞行阶段数;$N$ 为设定的大数。 -
至此,优化目标、优化参数、约束条件均设计完毕,参考轨迹设计问题转化为一般的优化问题,其数学描述如下
$$ \begin{split} & \mathop {\min }\limits_Z J\left( { Z} \right) \\ & {\rm{s}}{\rm{.t}}{\rm{.}} \\ & \left\{ \begin{aligned} & {{ X}_{k + 1}} = \varPhi \left( {{{ X}_k},{{ U}_k},{l_k},{ Z}} \right) \\ & {{ U}_k} = {f_u}\left( {{{ X}_k},{t_k},{l_k},{ Z}} \right) \\ & {l_{k + 1}} = {f_l}\left( {{{ X}_k},{t_k},{l_k},{ Z}} \right) \\ & {{ Z}_d} \leqslant { Z} \leqslant {{ Z}_u} \end{aligned} \right. \end{split} $$ (37) -
粒子群算法[4]示意图如图1所示:在优化参数空间中所有粒子一直并行搜索;每个控制周期搜索过程中,记录下每个粒子的迄今最优位置
${p^i}\left( k \right)$ 以及整个粒子群迄今最优位置${p^g}\left( k \right)$ ;通过速度更新与位置更新得到下一控制周期的搜索方向与位置。图1中黑点表示各粒子的当前位置,黑色带实线的箭头表示各粒子的当前搜索方向,白点表示某粒子的所有历史位置,白色箭头表示该粒子的所有历史搜索方向。本文采取速度更新的改进形式[5],速度更新与位置更新方程为
$$ \left\{ \begin{aligned} w_j^i\left( {k + 1} \right) = & w_j^i\left( k \right) + {\lambda _j}{\rm{rand}} \left( {0,1} \right)\left[ {p_j^i\left( k \right) - q_j^i\left( k \right)} \right] +\\ & {\eta _j}{\rm{rand}} \left( {0,1} \right)\left[ {p_j^g\left( k \right) - q_j^i\left( k \right)} \right] \\ q_j^i\left( {k + 1} \right) = & q_j^i\left( k \right) + w_j^i\left( {k + 1} \right) \\ \end{aligned} \right.$$ (38) 其中:
${q^i}\left( k \right)$ 表示第$i$ 个粒子当前(第$k$ 个控制周期搜索)的位置;${w^i}\left( k \right)$ 表示该粒子当前的搜索方向;${\lambda _j}$ 与${\eta _j}$ 分别为控制个体认知分量和群体社会分量相对贡献的第$j$ 维学习率;$j$ 表示搜索参数的第$j$ 维。 -
以月球南极某陨石坑探测为例,规划轨迹及各状态量与控制量如图2~8所示,可以发现:转段过程中飞行轨迹与各状态量平稳过渡。
图2中粗实线为地形曲线,虚线为规划轨迹,该图表明规划的轨迹远离月表,保证了飞行过程的安全性。图3表明由于引入了制导律,规划轨迹的航程误差很小,远远高于一般的需求指标。图4表明缓速下降段垂向速度误差很小,水平速度接近0,2.4节末尾关于xf ≈ H(y4)的假设合理。
在普通计算机Windows XP系统中用Matlab进行仿真,规划算法平均耗时40 s,能满足在线规划要求。另外,根据实际需求,在C环境中对程序进行优化可进一步提高该方法的实时性。
本文轨迹规划过程中引入制导律,搜索结果自动满足航程误差要求,简化了优化目标,提高了其余两项目标(燃料消耗、飞行安全)的搜索效率;在无动力滑行段预测转段条件,降低了优化参数的维数,提高了优化速度。
-
针对月球陨石坑探测任务,本文提出了基于粒子群算法的在线轨迹规划方法,同时考虑了燃料消耗、地形不确定性与导航误差对飞行安全的影响、姿轨控能力等,实时性也得到了仿真验证,在提高探测器自主飞行能力、应对异常事件发生等方面具有很好的工程参考价值。
Trajectory Planning for Lunar Crater Exploration Based on Particle Swarm Optimization
-
摘要: 针对月球陨石坑探测问题,提出一种基于粒子群算法的在线轨迹规划方法。该方法首次把飞行轨迹划分为垂直上升段、程序转弯段、无动力滑行段、接近段、缓速下降段5个阶段,根据各段任务特点设计转段条件与控制量形式,在此基础之上使用解析与数值计算得到一种状态快速更新方法,通过选取3个优化参数和2项优化目标,考虑各种约束,采用粒子群算法在线优化得到参考轨迹。该方法不同于一般的轨迹规划,在无动力滑行段和接近段引入制导律,以提高规划效率。仿真结果表明,该方法运算速度快,生成的轨迹适应能力强,有利于提高探测器自主生存能力,非常适合月球陨石坑探测任务。Abstract: A novel onboard trajectory planning methodology is proposed for lunar crater exploration. For the first time the flight trajectory is divided in to 5 phases: the vertical rise, the single axis rotation, the unpowered gliding, the approach and the vertical descent. According to different characters in these phases, transfer conditions and control modes are designed. Then a new method for state updating is developed in analytical and numerical ways. After that, 3 optimization parameters and 2 optimization objects are chosen. Considering the constraints, the reference trajectory is finally generated by the particle swarm optimization. Different from other planning, the guidance law is introduced during the unpowered gliding and the approach phase to improve optimization efficiency. Numerical simulations show that the methodology proposed is good real-time and the planned trajectory can meet all the requirements, which can improve the autonomous survivability of the explorer and is suitable for the lunar crater exploration.
-
Key words:
- particle swarm optimization /
- crater exploration /
- trajectory planning
Highlights● A novel onboard trajectory planning methodology is proposed for lunar crater exploration. ● For the first time the flight trajectory is divided in to 5 phases:the vertical rise,the single axis rotation,the unpowered gliding,the approach and the vertical descent. ● Different from other planning,the guidance law is introduced to improve optimization efficiency. -
[1] JAFFE L D, STEINBACHER R H. Surveyor program results final report: NASA-SP-184 69N36451[R]. Washington, DC: [s.n.], 1969. [2] 黄翔宇,张洪华,王大轶,等. “嫦娥三号”探测器软着陆自主导航与制导技术[J]. 深空探测学报,2014,1(1):52-59.HUANG X Y,ZHANG H X,WANG D Y,et al. Autonomous navigation and guidance for Chang’e-3 soft landing[J]. Journal of Deep Space Exploration,2014,1(1):52-59. [3] KOS L D, POLSGROVE T P, SOSTARIC R R, et al. Altair descent and ascent reference trajectory design and initial dispersion analyses[C]//AIAA Guidance, Navigation, and Control Conference. Toronto, Ontario, Canada: AIAA, 2010. [4] 孙增圻, 邓志东, 张再兴. 智能控制理论与技术[M]. 北京: 清华大学出版社, 2011. [5] 陈上上,何英姿,刘贺龙. 基于粒子群优化的再入飞行器在线轨迹规划[J]. 上海航天,2015,32(6):1-7. doi: 10.3969/j.issn.1006-1630.2015.06.001CHEN S S,HE Y Z,LIU H L. Onboard trajectory planning for entry vehicle based on particle swarm optimization[J]. Aerospace Shanghai,2015,32(6):1-7. doi: 10.3969/j.issn.1006-1630.2015.06.001 -