潜在威胁的小行星与地球轨道在接近或交叉的情况下，存在着撞击地球的风险。历史上地球曾多次遭到小行星撞击，部分撞击对地球造成了严重破坏，6 500万年前的希克苏鲁伯撞击事件曾导致全球60% ~ 80%的物种灭绝^[1]，因此围绕潜在威胁的小行星开展行星防御对人类长远生存发展具有重要的意义。评估核爆、动能撞击等防御手段需要精确的小行星星历与小行星质量、密度、自转等信息^[2]。受时空与天体运动的限制，现有的地基、天基观测手段所获取的小行星信息具有一定的局限性，对于庞大的小行星群体，通过实施小行星探测任务，人类仅对少数小行星的物理化学性质有过详细的探查^[3]。

Vetrisano等^[4]的研究表明，近距离抵近飞越是获取小行星特性信息的一种有效手段。飞越任务通常利用光学、光谱与红外远程成像系统对小行星进行探测，飞越探访可获取小行星表面的图像，测量小行星光谱，反演其物质组成、大小、体积、质量、密度和结构等来了解小行星的物理化学性质。针对有潜在撞地威胁的小行星，获取其精确的星历与质量、形状数据可准确计算撞击、偏转效应，对于制定有效的小行星防御策略具有重要的意义^[5]。联合国空间任务规划咨询组（Space Mission Planning Advisory Group，SMPAG）也在2021年呼吁加强对小行星飞越探测任务的支持，以获取小行星详细的目标特性信息^[6]。

“隼鸟号”（Hayabusa）、“黎明号”（Dawn）^[7-8]探测器对特定小行星进行了探测。小行星飞越探测可在单次任务中飞越数量可观的小行星，开展多目标特性普查，通过合理规划访问顺序与优化轨道可显著地降低其任务成本。小行星探测任务也由单一目标向多目标发展，美国国家航空航天局（National Aeronautics Space and Administration，NASA）于2021年10月发射了“露西号”（Lucy）探测器，计划在12 a的任务周期内飞越探访木星特洛伊族的8颗小行星，将创下单个航天器访问小行星数目最多的记录^[9]。日本宇宙航空研究开发机构（Japan Aerospace eXploration Agency，JAXA）开发的“命运﹢”（DESTINY+）探测器计划飞越多颗近地小天体并最终于2028年飞越双子座流星雨母体Phaethon 3200（TB 1983）^[10]。欧洲航天局（European Space Agency，ESA）计划于2029年发射彗星拦截器（Comet Interceptor），探测器将停泊在日地L2点周期轨道以期择机飞越新发现的长周期彗星^[11]。

既往小行星探测轨道设计多聚焦于交会探访任务，飞越探测仅作为整个任务的衍生。Gao等^[12]研究了以地球特洛伊族小行星2010TK7为主要探测目标的交会轨道设计，在抵达2010TK7前探测器飞越其它2颗小行星；夏炎等^[13]研究了主带小行星的可达性，提出了探测器最终伴飞小行星2004VZ60并在途中飞越4颗主带小行星的轨道设计；Qiao等^[14]研究了对近地小行星进行采样返回并在途中对不同光谱类型小行星进行飞越探测的轨道设计。McNutt等^[15]提出利用连续推力探测器对近地小行星进行近距离飞越探测的侦查兵（Near Earth Asteroid Scout，NEA Scout）计划。Greco等^[16]提出利用立方星双星系统对选定的近地天体飞越进行探访的任务，认为在近地天体过黄道面时对其进行飞越探访可有效节约改变轨道倾角、偏心率所导致的燃料消耗。Giuseppe等^[17]提出了针对Apollo型小行星的黄道面内多目标飞越序列规划方法，采取迭代扩大搜索区域的飞越序列拓展策略，一旦搜索区域中出现可行候选飞越目标就终止搜索，可导致解集的过度削减并遗漏最优解，需研究全局性更强、范围更广的多目标小行星飞越探测任务轨道优化设计方法。

多目标探测任务轨道优化的主要难点在于探测序列搜索。目前常用的路径搜索^[18-19]有：分支–限界算法、智能优化算法与混合优化方法。分支–限界算法通常以树结构表示解序列，通过不同策略或约束条件对树结构剪枝，可在较低计算成本下得出较好的优化结果。贪婪算法可视为分支–限界算法的一种特殊情况，在树结构的每层均选择局部最优值的策略易使其陷入局部最优。近年来，波束搜索算法作为分支–限界算法的一种，其易于部署、计算复杂度可由波束宽度控制的特点，在国际轨道优化竞赛当中展示了良好的性能^[20-21]。以蚁群算法为代表的智能优化算法需要设计者对信息素矩阵、启发式信息矩阵及算法参数设计有较深的认识，该方法在多目标交会任务当中应用广泛。混合优化方法，通常将多种路径搜索算法、人工筛选方法与结合专家知识的设计方法混合使用，部分还涉及到求解最优序列后进一步的局部优化^[22]，要求具有较强的专家经验。

不同于多目标交会任务，飞越探测任务在飞越时刻探测器与目标的速度并不相等，这使得两条不同飞越序列的相同子序列如$ C\to A\to B\to D $与$ E\to A\to B\to D $当中，由$ A\to B $的转移速度增量不一致。因此，在交会类任务当中被广泛使用的多种轨道转移速度增量近似计算方法，Edelbaum等^[23]估算方法无法在飞越类任务当中使用。转移速度增量的不恒定也会导致以启发式信息矩阵、信息素矩阵为群体智能算法进化依据的蚁群算法在多目标飞越问题上表现不佳。在交会任务中常用的分割轨迹^[24]方法，将完整的飞越任务分为多个子飞越任务，会遇到无法确定子轨迹速度初始值的问题。

本文本文的潜在威胁小行星数据来自喷气推进实验室（Jet Propulsion Laboratory，JPL）小行星数据库，候选目标为全体2 073个潜在威胁小行星。在Greco等^[16]工作的基础上，将黄道面内飞越基本策略进一步应用于全体潜在威胁小行星，在分析潜在威胁小行星的轨道特点的情况下，设计了在黄道面内多目标飞越潜在威胁小行星的轨道优化模型；针对任务时间长、候选目标多导致的序贯飞越解空间巨大的问题，将黄道面内飞越策略与波束选择树搜索算法相结合，对多目标飞越任务的轨道进行了快速求解。研究了不同策略的性能表现，分析了该类任务发射时间窗口的普遍存在性；对比了波束选择树搜索算法与迭代扩大搜索区域算法^[17]在相同任务场景下的性能，并进一步对比分析了蚁群算法在黄道面内飞越探测模型下的性能表现及原因；结合运载火箭发射能力、卫星燃料质量比等工程约束，对2024—2028年不同发射时间进行优化计算，得到探测器飞越数颗小行星的序列集合，就有代表性的序列进行结果分析与展示，可为未来飞越探测任务提供设计参考。

近地天体是轨道近日点在1.3 AU以内的太阳系小行星。包含小行星和彗星两类，近地小行星占绝大多数。与地球轨道的最近距离在0.05 AU以内且绝对星等小于22等（等效直径大于140 m）的近地小行星定义为潜在威胁小行星，其轨道可近距离接近地球，并且大到足以在发生撞击时造成严重的区域破坏。潜在威胁小行星倾角主要分布在0° ～ 30°，轨道半长轴决定其轨道周期分布在0.5～4 a，如图1所示。

图  1  潜在威胁小行星轨道半长轴、倾角分布

Figure 1.  Distribution of semi-major axis and inclination angle of potentially hazardous asteroids' orbit

由于探测器改变轨道倾角的机动需要消耗的燃料远高于轨道面内机动，且无论何种轨道构型的潜在威胁小行星其运动轨迹总会经过黄道面，可得到一种确定性的多目标小行星飞越基本策略为探测器仅在固定轨道平面内进行轨道面内机动，调整自身轨道与穿过该面的小行星进行序贯飞越的基本策略。考虑到探测器需从地月系逃逸与潜在威胁小行星过黄道面的位置分布情况，黄道面是该固定轨道平面的理想选择，确定潜在威胁小行星在特定时间段穿越黄道面时刻与位置是该算法的基础。

计算潜在威胁小行星飞越黄道面的时刻和位置时，只考虑太阳中心引力体的引力作用，小行星轨道状态采用开普勒轨道根数表示$ \left\{a,e,i,\varOmega ,\omega ,\nu \right\} $，分别为半长轴、偏心率、轨道倾角、升交点赤经、近日点幅角与真近点角。潜在威胁小行星过黄道面位置、时刻求解利用Giuseppe等^[17]在飞越任务中使用的计算方法，借助近焦点坐标系向日心–黄道坐标系转换可以快速求取。

潜在威胁小行星过黄道面位置计算见式（1），其中$ r\left(\nu \right)=a(1-{e}^{2})/(1+e\mathrm{cos}\nu ) $，令$ z=0 $，有$ {\mathrm{tan}\nu }_{\mathrm{tr}}=-\mathrm{tan}\omega $，用$ {\nu }_{\mathrm{t}\mathrm{r}}^{1} $，$ {\nu }_{\mathrm{t}\mathrm{r}}^{2} $代表小行星两次穿过黄道面对应真近点角（$ {\nu }_{\mathrm{t}\mathrm{r}}^{1} < {\nu }_{\mathrm{t}\mathrm{r}}^{2} $）。

将求得的$ {\nu }_{\mathrm{t}\mathrm{r}}^{1} $，$ {\nu }_{\mathrm{t}\mathrm{r}}^{2} $带回式（1）可得到小行星从2个方向穿越黄道面位置$r_\mathrm{tr}^1$，$r_\mathrm{tr}^2$、速度$ {\nu }_{\mathrm{t}\mathrm{r}}^{1} $，$ {\nu }_{\mathrm{t}\mathrm{r}}^{2} $，由真近点角可求取潜在威胁小行星飞越黄道面的时刻${t}_{\mathrm{t}\mathrm{r}}^{\mathrm{I}}$，${t}_{\mathrm{t}\mathrm{r}}^{\mathrm{II}}$。全体潜在威胁小行星穿越黄道面位置如图2所示，其主要分布在距日心3 AU的距离内。

图  2  潜在威胁小行星穿越黄道面位置分布

Figure 2.  Location distribution of potentially hazardous asteroids crossing the ecliptic plane

为简化卫星转移弧段的计算，本文设定卫星在飞越上一颗小行星后立即变轨，至此在潜在威胁小速度增量可以利用兰伯特定理^[25]解出探测器在出发时刻、位置$Sat \langle t,r,v\rangle$到目标小行星过面时刻，位置$As{t}_{\mathrm next}\langle t,r,v\rangle$的转移速度增量$ \Delta v $来表示。但需要注意的是，由于重复计算兰伯算法消耗大量程序时间，在计算时进行合理的矩阵化与并行化可以显著提高程序运算效率。飞越黄道面位置时刻算法最终输出满足探测器最大机动速度增量约束的候选小行星即

其中：$ {v}_{\mathrm{i}\mathrm{n}} $为机动前速度；$ {v}_{\text{af}} $为机动后速度。

探测器在单次飞越完成后即开始对全体潜在威胁小行星进行单圈兰伯特转移计算，搜索下一次黄道面内飞越的候选小行星，除此之外不再对单次转移时间做其他约束。因为小行星过黄道面位置、时刻存在周期关系避免了重复计算，所以确定的黄道面内多目标小行星飞越的基本策略比传统依靠时间窗口确定兰伯特转移始末时刻的大规模计算大幅降低了复杂度。

确定黄道面内飞越小行星的基本原则与兰伯特转移弧段表示的单次飞越轨迹后，转化为如何在给定任务的约束条件下，飞越尽可能多的潜在威胁小行星。小行星飞越问题的可行解是一个典型的序贯结构，可利用决策树生成。树节点表示探测器从出发小行星经过兰伯特转移到飞越目标小行星，每个节点包括的信息：出发小行星编号${A}_{\text {Ns}}$、目标小行星编号${A}_{\text {Nt}}$，出发时刻$ {t}_{\text{Ns}} $、飞越时刻$ {t}_{\text{Nf}} $，本弧段速度增量$ {v}_{\text{Narc}} $，剩余任务时间$ {t}_{\mathrm{Re}} $，剩余机动速度增量$ {v}_{\text{Re}} $。

对于树结构的解探索，有基于宽度或深度优先搜索的全局遍历方法与基于波束生成的选择性方法。全局遍历方法在构造序列解受到多重约束的特定条件下，存在着在可接受的时间、空间复杂度下完成遍历的情况。本文将转移搜索区域调整为整个黄道面，搜索区域内的候选目标较多，且候选解随着序列长度的增加呈指数增长，难以在可接受的时间内完成对可行解的完全遍历，因此采用基于宽度优先的波束树搜索算法^[21]作为解优化的方法。

波束选择树搜索算法的子树拓展策略对于算法整体表现有着决定性的作用。树节点值函数$ r $作为子树拓展策略的直接参考，一般可由探测器飞越下一个子节点所需的速度增量$ {v}_{\text{Narc}} $或子节点飞越时刻$ {t}_{\text{Nf}} $定义。考虑到直接定义下的$ {v}_{\text{Narc}} $与$ {t}_{\text{Nf}} $在各弧段互相独立缺乏全局性，对以上两种树节点值函数做无量纲化^[26]的定义：剩余速度增量比$ {r}_{v}={v}_{\mathrm{R}\mathrm{e}}/{v}_{\mathrm{a}\mathrm{l}\mathrm{l}} $与剩余任务时间比$ {r}_{t}={t}_{\mathrm{Re}}/{t}_{\mathrm{a}\mathrm{l}\mathrm{l}} $，其中${v}_{\mathit{\text {Re}}}$为剩余速度增量，$ {v}_{\mathrm{a}\mathrm{l}\mathrm{l}} $为总速度增量，$ {t}_{\mathrm{Re}} $为任务剩余时间，$ {t}_{\mathrm{a}\mathrm{l}\mathrm{l}} $为任务总时间。

在只考虑速度增量或任务时间的拓展策略基础上，融合两种评价指标的拓展策略及其变体在波束树搜索处理空间序贯决策任务当中也有出色的效果^[21]。波束选择树搜索算法的子树拓展策略值函数可有表1的不同定义方式，其中$ {r}_{1} $、$ {r}_{2} $分别为仅考虑燃料的燃料最省策略与仅考虑剩余任务时间的时间最省策略，$ {r}_{3} $、$ {r}_{4} $、$ {r}_{5} $为将剩余速度增量与剩余任务时间综合考虑，分别求取平均值、加权值与最小值的拓展策略。

不同的子树拓展规则本质是对剩余资源的选择分配，策略的优劣需要结合任务场景分析，优先选择在当前任务节点中更为稀缺的资源，避免出现一种资源耗尽而另一种资源剩余较多，从而影响最终的结果。

波束选择树搜索通过限制生成树拓展子树时的宽度来达到剪枝的效果，本文采用的波束宽度${bw}$为10～10 000，在构造每一层的解序列后，根据节点的$ r $值排列所有可行解，若满足约束的单次转移兰伯特弧段数量小于波束宽度上限则拓展全部子弧段，若单次转移兰伯特弧段数量大于${bw}$则依据子树拓展规则对特定的子弧段进行宽度优先拓展，其余节点被剪枝删去。波束选择树搜索算法如表2所示。

1）生成根节点，根节点为地球E，设定发射时间与任务约束等任务指标；

2）拓展生成子节点，从根节点出发计算向目标小行星飞越弧段，如果飞越目标小行星满足速度增量约束则生成子节点，所有可行的子节点构成了搜索树的新一层解集；

3）排列与剪枝，对新一层子节点解集按照$ r $值进行排列，保留${bw}$个子节点，对该层其余子节点进行剪枝删除；

4）重复步骤2与3，直到违反剩余机动速度增量约束、任务总时间约束或无可行目标小行星。

波束选择树搜索算法流程如图3所示，最终根据总速度增量约束与总任务时间约束来终止子树的生成，待完整搜索树生成后，根据目标函数筛选出最优子树分支，本例目标函数设定为

图  3  波束选择树搜索算法流程图

Figure 3.  Flowchart of beam selection tree search algorithm

即选择最长的子树分支作为最优解，$ N $为探测器飞行单条轨迹飞越的小行星总数，但在单次空间任务生成的决策树中可能得到数条飞越小行星数目相同的轨迹，因此可用总速度增量或任务总时长占比作为除飞越目标个数外的评价指标供多个解飞越数量一致时筛选使用。

假设卫星总速度增量3.5 km/s，单次机动最大脉冲速度增量不大于1 km/s，任务周期10 a，两次飞越间隔时间不做约束。任务目标为在一次任务中飞越探测小行星数目最多，多个解同数目情况下按照式（4）选择剩余燃料多者，飞越最后一颗小行星时间作为任务终止时间。

假定发射时间选择在2025年1月1日，波束宽度设置为1 000。采用5种不同策略最终形成的解飞越个数、目标函数值、任务终止时间与剩余燃料的对比结果如表3所示。

结果显示，在该问题场景下以飞越目标个数为评价指标，相比任务时间，燃料无疑是更加“稀缺”的资源，策略2时间优先原则下的结果明显差于考虑燃料的其他原则，策略1、5的最优飞越结果一致也证明了这种观点；同时在波束宽度取值较大时，出现了不同策略最优序列相同的情况，说明该波束宽度对整个决策树的探索较为完全，决策树波束宽度${bw}$取值合理。

同时以波束宽度为10、50、200、500、5 000、10 000等为例在策略3下对比说明波束宽度的选择对结果的影响，发射时间设置为2025年1月1日，结果见表4。当波束宽度选择较小时，对于决策树的解探索不完全，得到的最优飞越个数相较于大波束宽度有一定差距；波束宽度增大，在飞越数量一致情况下最优解序列剩余速度增量变大，最终解序列则有一定的收敛趋势。

由于不同时刻潜在威胁小行星位置分布的差异，发射时间对于黄道面内多目标飞越小行星任务决策树的解生成与优化也有着重要影响，本文自2024年1月1日—2029年1月1日，以每6个月为一个发射窗口，选择策略$ {r}_{4} $，波束宽度设置为1 000，对比验证发射时间对于解影响，结果如表5所示。在不同发射窗口下的最优序列飞越小行星个数在18 ～ 21颗，该任务对发射窗口要求不敏感，不同发射时间的该类任务均可以飞越数量可观的小行星。

在表5的基础上，对发射窗口按月份进一步密集化搜索，发现于2027年9月1日发射的任务序列$ L $最终也可飞越21颗小行星，但剩余燃料多，为0.194 5 km/s。该条轨道飞越各个小行星时刻与各次机动速度增量见表6，其中小行星直径、小行星轨道距地球轨道最近距离由ESA小行星数据库获得^[27]，探测器轨迹与飞越小行星位置如图4所示。

图  4  探测器轨迹与飞越探访小行星

Figure 4.  Trajectories of spacecraft and visited asteroids

本文以序列$ L $为例考虑“长征三号乙”（CZ-3B）运载火箭发射探测器，探究单次任务中飞越小行星总数与剩余质量的关系，结果如图5所示。发射能量为零时“长征三号乙”运载火箭可提供约3 284 kg发射质量^[28]，结合工程需要，本方法可快速估计干重与飞越总数的关系，若在此发射条件下，干重为1 000 kg的探测器可飞越21颗小行星，干重为1 500 kg则可飞越16颗小行星。

图  5  探测器剩余质量与飞越总数关系图

Figure 5.  Relationship between spacecraft residual mass and total number of flybys

为比较本文方法与Giuseppe等^[17]迭代扩大搜索区域算法性能，在原方法搜索得到最优序列的同等任务条件下利用本文方法仿真，结果对比如表7所示。

——任务时间：2024年1月1日—2029年1月1日；

——单次机动最大脉冲速度增量：$ \Delta {v}_{\mathrm{max}}=0.5\; \mathrm{k}\mathrm{m}/\mathrm{s} $；

——单次任务总速度增量：${v}_{\mathrm{t}\mathrm{o}\mathrm{t}}\text{ = 5}.5 \; \mathrm{k}\mathrm{m}/\mathrm{s}$；

利用本文方法搜索，设置波束宽度$ bw $为1 000，选择策略$ {r}_{4} $得到的针对Apollo型小行星的最优多目标飞越轨道共可飞越25颗小行星，多于其方法的21颗小行星，所需总速度增量为5.12 km/s也少于其所需的5.30 km/s，本文方法在个人电脑单次仿真用时为1 277.81 s^[3]。仿真结果表明本文方法在保证快速性的同时，相比原方法搜索性能更强，更具全局性。

其次本文参考西安飞行控制中心在GTOC9（Global Trajectory optimization Competition 9th）中的方法^[29]，针对飞越场景改造的蚁群算法对黄道面内多目标小行星飞越模型进行优化。本文所用方法与燃料最优的贪婪算法、蚁群算法优化结果的对比如表8所示。

以2027年9月1日发射的序列$ L $为对照，多数情况下蚁群算法生成的序列飞越总数在17颗左右，需要重复调参并多次运行才得到飞越21颗小行星的解序列。而波束树搜索作为一种确定性算法，在确定节点拓展策略与波束宽度后结果稳定收敛。因此，对于多目标小行星飞越问题，本算法收敛稳定性优于蚁群算法。此外，蚁群算法收敛至飞越21颗小行星的序列所需总速度增量也在3.5～3.9 km/s，远大于本文方法。

蚁群算法在多目标小行星飞越问题上表现不佳，原因可能为飞越探访转移速度增量不恒定导致启发式信息矩阵不断变化、目标群体数量较大导致信息素矩阵过于稀疏更新迭代慢及其swap、2-opt等局部搜索算子作用有限。

本文研究了在单次任务中对多颗潜在威胁小行星进行多目标飞越探测的轨道优化设计。针对该任务背景解空间大、任务轨迹计算耗时长的问题，利用黄道面内飞越潜在威胁小行星的基本策略简化模型，通过波束选择树搜索算法建立了一种快速有效求解多目标小行星飞越任务轨道的优化方法。仿真结果显示，波束宽度与树节点选择策略共同影响算法优化结果，在时间、空间复杂度允许的前提下，加大波束宽度减少剪枝范围有助于算法寻优；在本文任务场景下，燃料相比任务时间为更稀缺资源，仅考虑任务时间的树节点选择策略结果差于综合考虑燃料、飞越时间的树节点选择策略。不同发射窗口最优序列飞越潜在威胁小行星的个数有差异，仿真结果显示2024—2028年发射的任务至少可飞越探测18颗潜在威胁小行星，进一步细分时间窗口优化可得到数条最优飞越21颗小行星的序列。与迭代扩大区域搜索算法、蚁群算法的对比证明了本方法对多目标小行星飞越问题良好的求解性能。黄道面内多目标小行星飞越轨道设计方法，可为未来实施对潜在威胁小行星物理化学性质飞越普查任务设计提供参考。

未来将从黄道面内飞越小行星模型入手，将探测器飞越小行星后即刻变轨改为其可在原轨道运行一段时间后变轨；不再限定探测器仅在黄道面内进行飞越，允许探测器变倾角机动对在单次任务当中飞越探测更多潜在威胁小行星进行下一步的研究。