燃料最优是目前世界各国发射火星探测器的一项重要指标。自Howard于1962年首次提出利用气动力辅助轨道转移以来^[1]，人们开始研究利用大气阻力对火星探测器进行减速，使其成为环火星探测器，这种方式能够大大减少火星探测器在捕获过程中的燃料消耗。火星探测中的气动辅助技术主要有气动捕获和气动刹车2种，其中气动捕获技术是指探测器仅在一次穿越火星大气后即可在大气阻力的作用下到达目标轨道，而气动刹车则需要多次穿越大气。相比于气动刹车来说，气动捕获可以大大缩短轨道转移的时间。气动捕获轨道的优化问题是关系到气动捕获能否成功以及最终效果的关键问题，其实际上是探测器在穿越火星大气时，以燃料最优或能量最优为性能指标，满足一系列状态约束和边界约束的非线性优化问题。将火星探测器由捕获轨道以接近最小的燃料消耗和相应的进入条件转移到指定的目标轨道，根据轨道力学和可控集理论可以求得入口条件和出口条件^[2]。气动捕获轨道的优化问题一般来说没有解析解，现有的轨迹优化方法主要包括间接法、直接法等^[3]。

在利用间接法解决轨迹优化问题时，首先需要依据Pontryagin极大值原理将其转化为一个Hamiltonian两点边值问题，然后进行求解^[4-5]。采用间接法求解轨迹优化问题能够保证最优解满足一阶最优性必要条件，但是间接法存在求解过程复杂、收敛域小、对初值依赖性高且难以估计等缺点^[6]。

相对于间接法来说，直接法在轨迹优化中的应用更为广泛。其利用离散变换将问题转化为非线性规划问题，然后根据优化目标对其进行寻优。其中配点法是广泛应用的间接法，它同时离散控制变量和状态变量，并通过多项式插值对控制变量和状态变量进行拟合，通过优化配点处的控制变量搜寻满足约束的最优轨迹^[7-8]。正交配点法是最常用的一种配点法，又称作伪谱法。配点一般选择正交多项式的根，利用全局插值多项式得到状态的表达式，通过对插值多项式求导将微分方程转化为代数约束，通过积分得到终端状态^[9-11]。根据所选用的插值多项式，可以将伪谱法分为Legendre伪谱法、Gauss伪谱法、Chebyschev伪谱法、Radau伪谱法^[6]。

近年来，启发式智能优化算法凭借其鲁棒性、适应性以及全局优化性能受到越来越多的关注，已经被应用到再入轨迹优化和高超声速飞行器上升轨迹的优化等问题中，其无需使用哈密顿函数、也不用计算其导数即可获得接近最佳的结果^[12-13]。另外，也有学者将启发式智能优化算法和配点法相结合，使用较高阶多项式来近似控制量，并将此控制变量在配点处进行离散，使用启发式智能算法优化正交配点的值^[14-15]。

基于鸽子独特的导航能力，文献[16]提出了一种新型的群体智能优化方法：鸽群算法。文献[16]利用所提出的鸽群算法给出了空间机器人路径规划的优化方法。文献[17]引入捕食逃逸机制对鸽群算法进行改进，并将其用于无人机紧密编队协同控制的内环控制器的优化设计。文献[18]提出了基于自适应种群变异的鸽群算法，并将其应用于航天器集群的轨道优化，提高了种群演化深度和收敛速度。文献[19]则将鸽群算法与正交配点法相结合，用于再入航天器的轨道设计。文献[20]针对基本鸽群算法的寻优缺陷利用量子行为规则对其进行了改进，并利用改进后的鸽群算法用于无人机紧密编队控制。

本文考虑火星探测器的气动捕获轨道的优化问题。以火星探测器穿越火星大气时的动力学方程为基础，分析并确定了能量最优气动捕获轨道，给出了优化的性能指标。以倾侧角为控制序列，将气动捕获轨道优化问题转化为多参数优化问题。为了解决该多参数优化问题，本文通过引入参数对基本鸽群算法进行了改进，并将其用于火星探测器的气动捕获轨道优化。

1 气动捕获轨道优化问题建模 1.1 大气内飞行动力学方程的建立

假设火星大气相对于火星静止，并且认为飞行过程中探测器的质量变化可忽略，则探测器在火星大气内飞行的动力学方程为^[21]

$ \left\{ {\begin{array}{*{20}{l}} {\dot r = v{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }\\ {\dot v = - D - {g_{\rm{M}}}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }\\ {\dot \gamma = \frac{1}{v}L{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \sigma + \left( {\frac{v}{r} - \frac{{{g_{\rm{M}}}}}{v}} \right){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }\\ {\dot \lambda = \frac{{v{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi }}{{r{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta }}}\\ {\dot \theta = \frac{v}{r}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi }\\ {\dot \psi = \frac{{L{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \sigma }}{{v{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }} + \frac{v}{r}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma {\rm{tan}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta } \end{array}} \right. $

（1）

式中：探测器的状态由6个变量r、λ、θ、v、γ、ψ(分别为火心距、经度、纬度、飞行速度、飞行路径角及偏航角)确定；g_M为火星重力加速度；σ为滚转角；D为阻力加速度；L为升力加速度。D、L、g_M的表达式分别为

$ \left\{ {\begin{array}{*{20}{l}} {D = \rho {v^2}{S_{\rm{r}}}{C_D}/(2m)}\\ {L = \rho {v^2}{S_{\rm{r}}}{C_L}/(2m)}\\ {{g_{\rm{M}}} = \frac{\mu }{{{r^2}}}} \end{array}} \right. $

式中：S_r为探测器的有效截面积；m为探测器的质量；C_D为阻力系数；C_L为升力系数；μ为火星引力常数；ρ为火星大气密度，不同高度的火星大气密度不同，火星大气密度与高度之间的关系为

$ \rho = {\rho _0}{\rm{exp}}( - h/{h_{\rm{s}}}) $	（2）

式中：h_s=8 805.7 m为火星大气标高；ρ₀=0.0147 4 kg/m³为火星大气参考密度。又由式(1)可知，火心距、飞行速度和飞行路径角不改变轨道平面，并且其变化与其他3个状态无关，所以当不考虑轨道平面的变化时，上述大气内飞行动力学方程可以简化为

$ \left\{ {\begin{array}{*{20}{l}} {\dot r = v{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }\\ {\dot v = - D - {g_{\rm{M}}}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }\\ {\dot \gamma = \frac{1}{v}L{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \sigma + \left( {\frac{v}{r} - \frac{{{g_{\rm{M}}}}}{v}} \right){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma } \end{array}} \right. $

（3）

1.2 能量最优气动捕获轨道的确定

利用火星大气对探测器进行减速，通过调整姿态进而改变所受升力的大小来控制飞行轨迹，可以认为此过程是不消耗能量的无动力飞行过程，这里的能量最优指的是由捕获轨道转移到目标轨道时所消耗的能量最少。由捕获轨道转移到目标轨道所需要的速度增量为^[2]

$ \Delta {V_1} = \sqrt {2\mu } \left( {\sqrt {\frac{1}{{{r_{{\rm{ca}}}}}} - \frac{1}{{{r_{{\rm{ca}}}} + {r_{{\rm{tp}}}}}}} - \sqrt {\frac{1}{{{r_{{\rm{ca}}}}}} - \frac{1}{{{r_{{\rm{ca}}}} + {r_{{\rm{cp}}}}}}} } \right) $	（4）

$ \Delta {V_2} = \sqrt {2\mu } \left( {\sqrt {\frac{1}{{{r_{{\rm{tp}}}}}} - \frac{1}{{{r_{{\rm{ta}}}} + {r_{{\rm{tp}}}}}}} - \sqrt {\frac{1}{{{r_{{\rm{tp}}}}}} - \frac{1}{{{r_{{\rm{ca}}}} + {r_{{\rm{tp}}}}}}} } \right) $	（5）

总速度增量为

$ \Delta V = \Delta {V_1} + \Delta {V_2} $	（6）

则

$ \frac{{\partial (\Delta V)}}{{\partial ({r_{{\rm{ca}}}})}} = \frac{{\partial (\Delta {V_1})}}{{\partial ({r_{{\rm{ca}}}})}} + \frac{{\partial (\Delta {V_2})}}{{\partial ({r_{{\rm{ca}}}})}} $	（7）

又

$ \begin{array}{*{20}{c}} {\frac{{\partial (\Delta {V_1})}}{{\partial ({r_{{\rm{ca}}}})}} = \frac{{\sqrt {2\mu } }}{2}\sqrt {{r_{{\rm{cp}}}}} \frac{{2{r_{{\rm{ca}}}} + {r_{{\rm{cp}}}}}}{{{{(r_{{\rm{ca}}}^2 + {r_{{\rm{cp}}}}{r_{{\rm{ca}}}})}^{\frac{3}{2}}}}} - }\\ {\frac{{\sqrt 2 \mu }}{2}\sqrt {{r_{{\rm{tp}}}}} \frac{{2{r_{{\rm{ca}}}} + {r_{{\rm{tp}}}}}}{{{{(r_{{\rm{ca}}}^2 + {r_{{\rm{tp}}}}{r_{{\rm{ca}}}})}^{\frac{3}{2}}}}}} \end{array} $

（8）

$ \frac{{\partial (\Delta {V_2})}}{{\partial ({r_{{\rm{ca}}}})}} = \left\{ {\begin{array}{*{20}{l}} { - \frac{{\sqrt {2\mu } }}{2} \cdot \frac{{\sqrt {{r_{{\rm{tp}}}}} }}{{\sqrt {{r_{{\rm{ca}}}}} {{({r_{{\rm{ca}}}} + {r_{{\rm{tp}}}})}^{\frac{3}{2}}}}}}&{{r_{{\rm{ca}}}} < {r_{{\rm{ta}}}}}\\ {\frac{{\sqrt {2\mu } }}{2} \cdot \frac{{\sqrt {{r_{{\rm{tp}}}}} }}{{\sqrt {{r_{{\rm{ca}}}}} {{({r_{{\rm{ca}}}} + {r_{{\rm{tp}}}})}^{\frac{3}{2}}}}}}&{{r_{{\rm{ca}}}} > {r_{{\rm{ta}}}}} \end{array}} \right. $

（9）

式中：r_ca为捕获轨道远火点半径；r_cp为捕获轨道近火点半径；r_ta为目标轨道远火点半径；r_tp为目标轨道近火点半径。由式(9)可知，对于给定的r_cp，ΔV在区间r_ca∈[0, r_ta]内单调递减，在区间r_ca∈[r_ta, +∞]内单调递增。因此，当r_ca=r_ta时，ΔV取得最小值，即当探测器穿越大气后的后捕获轨道的远地点高度与目标轨道高度相等，也就是二者相切时，轨道转移所需的能量消耗最少。

由动量矩守恒原理和能量守恒原理可知

$ {r_{{\rm{ta}}}}{v_{{\rm{ca}}}} = {r_{\rm{a}}}{v_{\rm{f}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}} $	（10）

式中：v_f为探测器飞出大气时的速度；γ_f为探测器飞出大气时的航迹角；r_a为火星大气的高度；v_ca为捕获轨道远火点处速度。根据能量守恒有

$ \frac{{v_{{\rm{ca}}}^2}}{2} - \frac{\mu }{{{r_{{\rm{ta}}}}}} = \frac{{v_{\rm{f}}^2}}{2} - \frac{\mu }{{{r_{\rm{a}}}}} $	（11）

联立式(10)和式(11)可得

$ \frac{1}{2}{\left( {\frac{{{r_{\rm{a}}}}}{{{r_{{\rm{ta}}}}}}{v_{\rm{f}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}}} \right)^2} - \frac{\mu }{{{r_{{\rm{ta}}}}}} = \frac{{v_{\rm{f}}^2}}{2} - \frac{\mu }{{{r_{\rm{a}}}}} $	（12）

由式(12)可得

$ {v_{\rm{f}}} = \sqrt {\frac{{2\mu {r_{{\rm{ta}}}}({r_{{\rm{ta}}}} - {r_{\rm{a}}})}}{{{r_{\rm{a}}}[r_{{\rm{ta}}}^2 - {{({r_{\rm{a}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}})}^2}]}}} $	（13）

因此，当捕获轨道远火点与目标轨道远火点相切时有

$ \Delta V = {v_{{\rm{ta}}}} - \sqrt {\frac{{2\mu {r_{{\rm{ta}}}}({r_{{\rm{ta}}}} - {r_{\rm{a}}})}}{{{r_{{\rm{ta}}}}[r_{{\rm{ta}}}^2 - {{({r_{\rm{a}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}})}^2}]}}} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}} $	（14）

式中:v_ta为目标轨道远火点处的速度，当目标轨道确定时其为一常值。由式(14)可知，探测器到达目标轨道高度时所需脉冲变轨的大小与飞出角度γ_f密切相关。可见，飞出角度γ_f越小，在到达目标轨道高度时所需的脉冲ΔV也就越小。

气动捕获轨道优化问题实质上就是对探测器穿越火星大气时的轨迹进行控制，使得探测器在飞出大气时恰好到达或接近能量最优气动捕获轨道。由前述分析可知，要使气动捕获轨道为最优气动捕获轨道，则气动捕获轨道的远地点应与目标轨道相切，飞出时的航迹角尽量小，且飞出速度与飞出时的航迹角满足式(13)。取最优控制指标为

$ J = {\left\{ {{v_{{\rm{ta}}}} - \sqrt {\frac{{2\mu {r_{\rm{a}}}({r_{{\rm{ta}}}} - {r_{\rm{a}}})}}{{{r_{{\rm{ta}}}}(r_{{\rm{ta}}}^2 - {{({r_{\rm{a}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}})}^2})}}} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}}} \right\}^2} + {({v_{\rm{f}}} - v_{\rm{f}}^*)^2} $	（15）

式中：

$ v_{\rm{f}}^* = \sqrt {\frac{{2\mu {r_{{\rm{ta}}}}({r_{{\rm{ta}}}} - {r_{\rm{a}}})}}{{{r_{\rm{a}}}(r_{{\rm{ta}}}^2 - {{({r_{\rm{a}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}})}^2})}}} $

为对应飞出航迹角满足式(13)的飞出速度。进入条件为

$ \left\{ {\begin{array}{*{20}{l}} {{r_0} = {r_{\rm{a}}}}\\ {{v_0} = v_0^*}\\ {{\gamma _0} = \gamma _0^*} \end{array}} \right. $	（16）

式中：v₀^*为进入大气的速度；γ₀^*为进入航迹角。出口条件为

$ \left\{ {\begin{array}{*{20}{l}} {{r_{\rm{f}}} = {r_{\rm{a}}}}\\ {0 < {\gamma _{\rm{f}}} < {\gamma _{{\rm{fmax}}}}} \end{array}} \right. $	（17）

式中：γ_fmax为保证捕获轨道与目标轨道相切的最大飞出航迹角。探测器大气内飞行时的总飞行约束为

$ \left\{ {\begin{array}{*{20}{l}} {n - {n_{{\rm{max}}}} \le 0}\\ {\dot Q - {{\dot Q}_{{\rm{max}}}} \le 0}\\ {{h_{{\rm{max}}}} - h \le 0}\\ {{\gamma _{\rm{f}}} - {\gamma _{{\rm{ fmax }}}} \le 0} \end{array}} \right. $	（18）

为将气动捕获轨道优化问题转化为多参数优化问题，选取倾侧角σ作为控制变量，将控制变量在配点处进行离散，并以这些离散点作为节点构造Lagrange插值多项式来近似拟合控制变量，选取Chebyschev多项式的根作为配点，有

$ {\tau _1} = {\rm{cos}}(\pi l/n)\quad l = 0,1, \cdots ,n $	（19）

式中：n为插值多项式的阶数，并且插值点落在[－1, 1]之间，对t∈(t₀, t_f)进行时域转换

$ t = \left( {\frac{{{t_0} - {t_{\rm{f}}}}}{2}} \right)\tau + \left( {\frac{{{t_0} + {t_{\rm{f}}}}}{2}} \right) $	（20）

在这里由于大气内飞行的时间未知，所以将飞出时间t_f也作为待优化的变量。由于控制量σ是关于时间t的连续函数，很难得到其精确表达式。为此，按照式(19)和式(20)的时间转换方式对控制量σ进行离散化，即

$ {\sigma _l} = \sigma (t({\tau _l}))\quad l = 0,1, \cdots ,n $

因此优化参数为

$ \mathit{\boldsymbol{\xi }} = {\left[ {\begin{array}{*{20}{l}} {{t_{\rm{f}}}}&{{\sigma _0}}&{{\sigma _1}}& \cdots &{{\sigma _n}} \end{array}} \right]^{\rm{T}}} $	（21）

于是，对于气动捕获轨道的优化问题就转化为对大气内飞行时间t_f以及控制序列σ在配点处的离散值的参数优化问题。针对这个多参数优化问题，第2节将提出一种改进鸽群算法对上述待优化的参数向量进行优化。

2 鸽群算法的分析与改进 2.1 原始鸽群算法及其存在的问题

鸽群算法是受鸽子觅食的启发提出的一种仿生学算法^[16]。当鸽子距离目标较远时，通过感知磁场在脑海中形成地图，来不断接近目标；随着鸽子不断接近目的地，其导航工具由磁场变为目的地附近的地标。熟悉地标的鸽子直接飞往目的地，其他鸽子则跟随那些熟悉地标的鸽子向目的地飞行。

设某优化问题的目标函数为J(ξ)，其中优化变量为ξ。鸽群算法在处理优化问题时模拟鸽子导航的过程。假设在优化时需要N只鸽子，每只鸽子的位置X_i, i=1, 2, …, N表示优化变量ξ的值，在优化时化还需要对应鸽子的速度V_i, i=1, 2, …, N假设优化变量是m维的，则每只鸽子的位置和速度分别记作

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{X}}_i} = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{x}}_{i1}}}&{{\mathit{\boldsymbol{x}}_{i2}}}& \cdots &{{\mathit{\boldsymbol{x}}_{im}}} \end{array}} \right]}\\ {{\mathit{\boldsymbol{V}}_i} = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{v}}_{i1}}}&{{\mathit{\boldsymbol{v}}_{i2}}}& \cdots &{{\mathit{\boldsymbol{v}}_{im}}} \end{array}} \right]} \end{array}} \right. $

在采用鸽群算法优化时，设鸽群中有N只鸽子，它们的位置和速度分别为X_i和V_i。该算法求解优化问题时包括2部分迭代过程，分别称作地图迭代和地标迭代。在地图迭代中，鸽群中的每一个体通过种群中的最优信息来更新自身的位置和速度，更新公式为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{V}}_i}(t) = {\mathit{\boldsymbol{V}}_i}(t - 1){{\rm{e}}^{ - Rt}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{rand }} \cdot [\mathit{\boldsymbol{G}}(t - 1) - {\mathit{\boldsymbol{X}}_i}(t - 1)]}\\ {{\mathit{\boldsymbol{X}}_i}(t) = {\mathit{\boldsymbol{X}}_i}(t - 1) + {\mathit{\boldsymbol{V}}_i}(t)} \end{array}} \right. $

（22）

式中：R为迭代因子；rand为[0, 1]区间上的随机值；t=1, 2, …, t_max为当前迭代次数，t_max为迭代总次数；i=1, 2, …, N代表第i只鸽子；V_i(t)为第i个鸽子在第t次迭代的速度；X_i(t)为第i个鸽子在第t次迭代的位置；G(t－1)为X_i(t－1) (i=1, 2, …, N)中的某一个，其满足：

$ J(\mathit{\boldsymbol{G}}(t - 1)) = \mathop {{\rm{min}}}\limits_{i \in \{ 1,2, \ldots ,N\} } J({\mathit{\boldsymbol{X}}_i}(t - 1)) $

当迭代次数到达最大迭代次数的75%以后，进入地标迭代的过程，选取X_c(t)为第t次迭代鸽群的中心点，作为鸽子飞行的参考方向对每一只鸽子进行更新，更新公式为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{X}}_{\rm{c}}}(t) = \frac{{\sum\nolimits_{i = 1}^{N(t)} {{\mathit{\boldsymbol{X}}_i}} (t - 1){\rm{fitness}}[{\mathit{\boldsymbol{X}}_i}(t - 1)]}}{{N(t)\sum\nolimits_{i = 1}^{N(t)} {{\rm{fitness}}} [{\mathit{\boldsymbol{X}}_i}(t - 1)]}}}\\ {{\mathit{\boldsymbol{X}}_i}(t) = {\mathit{\boldsymbol{X}}_i}(t - 1) + {\rm{rand}} [{\mathit{\boldsymbol{X}}_{\rm{c}}}(t) - {\mathit{\boldsymbol{X}}_i}(t - 1)]} \end{array}} \right. $

（23）

式中：

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{ fitness }}[{\mathit{\boldsymbol{X}}_i}(t)] = \frac{1}{{J({\mathit{\boldsymbol{X}}_i}(t)) + \varepsilon }}}\\ {N(t) = \frac{{N(t - 1)}}{2}} \end{array}} \right. $

其中：fitness[X_i(t)]为当前迭代下鸽子的适应度；N(t)为第t次迭代下种群的数量；ε为防止分母为零而附加的充分小的数。

可见，在地标迭代过程中，每次迭代舍去部分个体，既保证了算法的较优信息又提高了迭代收敛的速率。原始鸽群算法的全局搜索能力主要体现在地图迭代的过程中，动态权值e^－RN体现了算法的探索能力，在鸽群向当前最优值收敛的过程中保持鸽群的多样性。然而这种方式具有局限性，在鸽群数目一定的情况下，R值的选取在很大程度上决定了全局搜索的结果。

图 1为R在不同取值下动态因子的变化情况。可以看出，当R取较大值0.3时，迭代到20次e^-Rt就已衰减接近为零，之后的迭代过程主要体现在向最优值靠近，失去了探索能力，算法易早熟；当R取较小值0.002 6时，在迭代完成时e^-Rt仍然具有较大的取值，算法不易收敛；为了平衡二者的关系，可以通过指定迭代完成时期望的动态因子的取值自适应地选取R值。令迭代完成时期望的动态因子的值为α, 有e^-RT=α, 则R=－(1/T)lnα, 其中T为迭代次数。

通过自适应的方式选取R值，虽然能够在一定程度上缓解易早熟和不易收敛的问题，但是却不能从根本上很好地平衡全局搜索与局部搜索之间的关系。在实际应用中，在迭代初期，往往更注重于全局搜索能力，而在迭代后期则希望局部搜索能力更强一些。由图 1可以看出，当取α=0.01时，在迭代后期动态因子较小，算法主要在最优值附近进行局部搜索，然而在迭代初期，动态因子下降很快，并不能很好地进行全局搜索；当取α=0.4时，虽然使得动态因子下降速度减缓，全局搜索能力增强，但是却削弱的算法的局部搜索能力。

2.2 改进的鸽群算法

本节在原有鸽群算法的负指数动态权值的基础上，为了平衡全局搜索与局部搜索提出一种广义的动态权值项。即将式(23)修改为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{V}}_i}(t) = {\mathit{\boldsymbol{V}}_i}(t - 1)f(t) + {\rm{rand}} [G(t - 1) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{X}}_i}(t - 1)]}\\ {{\mathit{\boldsymbol{X}}_i}(t) = {\mathit{\boldsymbol{X}}_i}(t - 1) + {\mathit{\boldsymbol{V}}_i}(t)} \end{array}} \right. $

（24）

式中：

$ f(t) = \frac{1}{{k + {{\rm{e}}^{\omega (t - b)}}}} $	（25）

且有

$ 0 \le f(t) \le 1 $	（26）

为了更加方便直观地确定参数ω、b的取值，引入另外2组参数(t₁, α₁)、(t₂, α₂)，令

$ \left\{ {\begin{array}{*{20}{l}} {f({t_1}) = {\alpha _1}}\\ {f({t_2}) = {\alpha _2}} \end{array}} \right. $	（27）

表示当迭代次数为t₁时，动态权值取值为α₁；当迭代次数为t₂时，动态权值取值为α₂。即

$ \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{k + {{\rm{e}}^{\omega ({t_1} - b)}}}} = {\alpha _1}}\\ {\frac{1}{{k + {{\rm{e}}^{\omega ({t_2} - b)}}}} = {\alpha _2}} \end{array}} \right. $	（28）

求解可得

$ \left\{ {\begin{array}{*{20}{l}} {\omega = \frac{1}{{{t_2} - {t_1}}}\left[ {{\rm{ln}}\left( {\frac{1}{{{\alpha _2}}} - k} \right) - {\rm{ln}}\left( {\frac{1}{{{\alpha _1}}} - k} \right)} \right]}\\ {b = \frac{1}{\omega } \cdot \frac{{{t_1}}}{{{t_2} - {t_1}}}\left[ {{\rm{ln}}\left( {\frac{1}{{{\alpha _2}}} - k} \right) - {\rm{ln}}\left( {\frac{1}{{{\alpha _1}}} - k} \right)} \right] - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{\omega }{\rm{ln}}\left( {\frac{1}{{{\alpha _1}}} - k} \right)} \end{array}} \right. $

（29）

因此，可以通过指定k、(t₁, α₁)、(t₂, α₂)来唯一确定动态权值。

下面对各参数的取值情况进行讨论。对式(25)求导可得

$ \dot f = \frac{{ - {{\rm{e}}^{\omega (t - b)}}\omega }}{{{{[k + {{\rm{e}}^{\omega (t - b)}}]}^2}}} $	（30）

由式(29)可知，当取t₁≤t₂且α₂≤α₁时，ḟ < 0，式(25)单调递减。

当k≥1时，对式(25)取极限可得

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{li}}{{\rm{m}}_{t \to \infty }}\frac{1}{{k + {{\rm{e}}^{\omega (t - b)}}}} = 0}\\ {{\rm{li}}{{\rm{m}}_{t \to - \infty }}\frac{1}{{k + {{\rm{e}}^{\omega (t - b)}}}} = \frac{1}{k} \le 1} \end{array}} \right. $	（31）

此时k、(t₁, α₁)、(t₂, α₂)可任意取值都可以满足式(26)的约束。

当k < 1时，此时k的取值与(t₁, α₁)、(t₂, α₂)相互制约，由式(26)可得

$ 0 < \frac{1}{{k + {{\rm{e}}^{ - \omega b}}}} \le 1 $	（32）

即

$ 1 - {{\rm{e}}^{ - \omega b}} \le k $	（33）

又由式(29)可得

$ {{\rm{e}}^{ - \omega b}} = \frac{{{{\left( {\frac{1}{{{\alpha _1}}} - k} \right)}^{1 + c}}}}{{{{\left( {\frac{1}{{{\alpha _2}}} - k} \right)}^c}}} $	（34）

式中: $ c = \frac{{{t_1}}}{{{t_2} - {t_1}}} $，将式(34)代入式(33)得

$ k \ge 1 - \frac{{{{\left( {\frac{1}{{{\alpha _1}}} - k} \right)}^{1 + c}}}}{{{{\left( {\frac{1}{{{\alpha _2}}} - k} \right)}^c}}} $	（35）

不妨取t₁=T/2、t₂=T，则有

$ k \ge 1 - \frac{{{{\left( {\frac{1}{{{\alpha _1}}} - k} \right)}^2}}}{{\frac{1}{{{\alpha _2}}} - k}} $	（36）

因为k < 1且0≤α₂≤1，所以 $ \frac{1}{{{a_2}}} - k > 0 $，有

$ \frac{1}{{{\alpha _2}}}k - {k^2} \ge \frac{1}{{{\alpha _2}}} - k - \left[ {{{\left( {\frac{1}{{{\alpha _1}}}} \right)}^2} - \frac{{2k}}{{{\alpha _1}}} + {k^2}} \right] $	（37）

即

$ \left( {1 - \frac{2}{{{\alpha _1}}} + \frac{1}{{{\alpha _2}}}} \right)k \ge \frac{1}{{{\alpha _2}}} - {\left( {\frac{1}{{{\alpha _1}}}} \right)^2} $	（38）

又

$ \left\{ {\begin{array}{*{20}{l}} {k < 1}\\ {1 - \frac{2}{{{\alpha _1}}} + \frac{1}{{{\alpha _2}}} > \frac{1}{{{\alpha _2}}} - {{\left( {\frac{1}{{{\alpha _1}}}} \right)}^2}} \end{array}} \right. $	（39）

当 $ 1 - \frac{2}{{{a_1}}} + \frac{1}{{{a_2}}} > 0 $时，有

$ \left\{ {\begin{array}{*{20}{l}} {{\alpha _1} > \frac{{2{\alpha _2}}}{{1 + {\alpha _2}}}}\\ {k \ge \frac{{\frac{1}{{{\alpha _2}}} - {{\left( {\frac{1}{{{\alpha _1}}}} \right)}^2}}}{{1 - \frac{2}{{{\alpha _1}}} + \frac{1}{{{\alpha _2}}}}}} \end{array}} \right. $	（40）

当 $ 1 - \frac{2}{{{a_1}}} + \frac{1}{{{a_2}}} < 0 $时，有

$ \left\{ {\begin{array}{*{20}{l}} {{\alpha _1} < \frac{{2{\alpha _2}}}{{1 + {\alpha _2}}}}\\ {k < 1} \end{array}} \right. $	（41）

取α₂=α₁²时，有

$ {\alpha _1} - \frac{{2{\alpha _2}}}{{1 + {\alpha _2}}} = \frac{{{\alpha _1}}}{{1 + \alpha _1^2}}{(1 - {\alpha _1})^2} > 0 $	（42）

由式(40)可知此时k≥0。当取k=0时，由式(29)有b=0，又由式(25)可得

$ f(t) = \frac{1}{{{{\rm{e}}^{\omega t}}}} $	（43）

此时改进算法退化为原始的鸽群算法，也就是说原始鸽群算法是改进后算法的一种特殊情况。

改进后的鸽群算法能够通过调整动态权值很好地平衡全局搜索与局部搜索之间的关系，既保证了算法的收敛性，又大大加强了算法的探索能力。图 2为改进后的鸽群算法在k、(t₁=N/2, α₁)、(t₂=N, α₂=0.1)的不同取值下动态权值的变化情况，其中N=200。通过调节α₁和k的值可得到不同的动态权值曲线。

可以看出，当k=0且b=0时，算法即为原始的鸽群算法，即原算法是改进后的一种特殊情况。当k≠0时，可以通过调节(t₁, α₁)、(t₂, α₂)选择合适的动态权值曲线，在搜索前期使动态权值缓慢下降，保持搜索活力，提高全局搜索的能力；在搜索后期，使得动态权值取得较小的值，保证算法的收敛。相比于原算法，改进后的算法不仅兼容了原算法中指数形式下降的动态权值，而且对其进行了扩展，在形式选择上提供了更大的自由度，提升了全局搜索与局部搜索之间的平衡。

3 气动捕获轨道优化

采用改进的鸽群算法对气动捕获轨道进行优化。鸽群中的每一只鸽子所在的位置代表所求问题的一个可能解，其初值在搜索空间内随机生成。选择鸽群中鸽子的数目为m，则鸽群的位置X可表示为

$ \mathit{\boldsymbol{X}} = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{X}}_1}}&{{\mathit{\boldsymbol{X}}_2}}& \cdots &{{\mathit{\boldsymbol{X}}_m}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{t_{f1}}}&{{t_{f2}}}& \cdots &{{t_{fm}}}\\ {{\sigma _{01}}}&{{\sigma _{02}}}& \cdots &{{\sigma _{0m}}}\\ {{\sigma _{11}}}&{{\sigma _{12}}}& \cdots &{{\sigma _{1m}}}\\ \vdots & \vdots & \ddots & \vdots \\ {{\sigma _{n1}}}&{{\sigma _{n2}}}& \cdots &{{\sigma _{nm}}} \end{array}} \right] $

（44）

式中：X_k(k=1, 2, …, m)为每一只鸽子当前所得到的待优化变量的值，包括探测器飞出大气的时间t_f，以及以倾侧角为控制序列在配点处的离散值。鸽群的更新速度V与X同维，首先对每一只鸽子所优化得到的控制量进行Lagrange插值得到控制序列，然后根据所得到的控制序列和优化得到的飞出时间t_f对动力学方程式(3)进行积分得到各个状态的值，计算此时的目标函数值，然后进行迭代寻优。

鸽群算法本身并不具备处理约束的能力，在进行迭代寻优的过程中为了保证满足过程约束和终端约束，将这些约束条件以罚函数的方式加到目标函数当中，经过修正的目标函数为

$ \begin{array}{*{20}{l}} {J = {p_1}{{\left\{ {{v_{{\rm{ta}}}} - \sqrt {\frac{{2\mu ({r_{{\rm{ta}}}} - {r_{\rm{a}}})}}{{{r_{{\rm{ta}}}}[r_{{\rm{ta}}}^2 - {{({r_{\rm{a}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}})}^2}]}}} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}}} \right\}}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {p_2}{{({v_{\rm{f}}} - v_{\rm{f}}^*)}^2} + {p_3}\sum\limits_{i = 1}^4 {{\rm{max}}} (0,{g_i}(x))} \end{array} $

式中：p₁、p₂和p₃为惩罚系数；g_i(x)为式(18)所示的飞行约束。如果取目标轨道为圆轨道，则有 $ {\upsilon _{{\rm{ta}}}} = \sqrt {\mu /{r_{{\rm{ta}}}}} $，此时目标函数式可简化为

$ \begin{array}{*{20}{l}} {J = {p_1}{{\left[ {1 - \sqrt {\frac{{2{r_{\rm{a}}}({r_{{\rm{ta}}}} - {r_{\rm{a}}})}}{{r_{{\rm{ta}}}^2 - {{({r_{\rm{a}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}})}^2}}}} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _{\rm{f}}}} \right]}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {p_2}{{({v_{\rm{f}}} - v_{\rm{f}}^*)}^2} + {p_3}\sum\limits_{i = 1}^4 {{\rm{max}}} (0,{g_i}(x))} \end{array} $

（45）

4 仿真分析

通过本章提出的改进鸽群算法，对火星探测器的气动捕获轨道进行优化。根据文献选择火星探测器的参数如表 1所示，其中m为探测器的质量，r_n为探测器曲率半径，S为探测器的有效截面积，C_L为升力系数，C_D为阻力系数；火星的参数如表 2所示。选取鸽群数目为N=150，迭代次数k=120，最大热流率 $ {{\dot Q}_{\max }} = 7{\rm{ \times }}{10^6}{\rm{W}}/{{\rm{m}}^2} $，最大过载n_max=4.5 g₀，最小飞行高度h_min=35 000 m，选取目标轨道为圆轨道，目标轨道高度r_t=3 895 000 m，初始进入条件r₀=3 520 000 m，v₀=5 000 m/s，γ₀=-0.17 rad，飞出条件r_f=3 520 000 m。

采用改进鸽群算法，取倾侧角为五阶插值多项式时，改进鸽群算法的迭代收敛曲线以及优化结果如图 3所示，图 3(a)为探测器飞行高度随时间变化的情况，可见进入高度和飞出高度相等，都等于火星大气高度，并且最小飞行高度大于所要求的35 000 m。图 3(b)为飞行速度随时间变化的曲线，飞出速度与飞出航迹角基本满足式(13)。图 3(c)为航迹角随时间变化的曲线，图 3(d)和图 3(e)为路径约束，分别表示过载和热流率随时间的变化情况，可见大气内飞行过程可以满足过程约束，图 3(f)为控制变量随时间变化的曲线。综上，通过改进鸽群算法，以倾侧角为控制变量可以得到一条接近能量最优的气动捕获轨道。

5 结论

1) 在对火星探测器的气动捕获过程分析的基础上，从脉冲变轨的角度给出了气动捕获轨道的优化指标。

2) 对原始鸽群算法的局限性进行分析，提出一种能够平衡收敛速度与全局搜索关系的改进算法。

3) 将气动捕获轨道优化问题转化为多参数优化问题，并利用提出的改进个鸽群算法对其进行优化。