刘世兴，史东华，陈菊，易中贵，王本亮，龚自正，郭永新*

(1.辽宁大学物理学院，沈阳 110036；2.辽宁大学空间科学与技术研究院，沈阳 110036；3.北京理工大学数学与统计学院，北京 100081；4.北京理工大学宇航学院，北京 100081；5.北京卫星环境工程研究所，北京 100094)

1 引言

“杞人忧天”出自《列子·天瑞》：杞国有人忧天地崩坠，身亡所寄，废寝食者。虽然这里的“杞人”是被人嘲讽的对象，但实际上，小行星撞击地球的事件在历史上时有发生。夜空划过的美丽流星，每年11月份出现的狮子座流星雨等都是流星颗粒进入地球大气层和空气摩擦而形成的。我们将围绕太阳运行、尺寸在1m～800km内，且不容易释放出气体和尘埃的天体，称作小行星。太阳系中存在大量的小行星，主要分布在火星与木星轨道之间的小行星带和海王星外的柯伊伯带，如图1所示。在天文学上，将距离地球最小距离在0.3AU(1AU=1.496×108km)范围内的小行星称为近地小行星。地球曾发生过22次不同程度的生物灭绝事件，其中至少10次以上是由小行星撞击地球所引起[1]，如6500万年前的恐龙灭绝事件，1908年的通古斯空爆事件，2013年的俄罗斯车里雅宾斯克地区空爆事件等[2-5]。仅2021年内全球已发生近地小行星飞掠地球事件1605次，观测到29颗小行星进入地球大气层发生火流星事件。我国境内也经常发生火流星事件，比如，2014年11月5日在内蒙古锡林郭勒，2017年10月4日在云南香格里拉，2018年在云南西双版纳，2019年10月11日在吉林松原地区，2020年在青海玉树地区，2021年在河南驻马店等发生了火流星事件。在2022年3月11日，编号为2022 EB5，直径大约3m的小行星在冰岛北部上空进入了地球大气层。虽然上述事件均未造成灾难，但已威胁到了人们的生产和生活安全。由此可见，近地小行星撞击风险概率虽小，但危害极大，但近地小行星撞击地球的危害可测可防。防范近地小行星撞击风险，具有非常迫切的现实需求和深远的战略意义[6-10]。

图1 太阳系中的八大行星和近地小行星带分布图

应对近地小行星撞击风险，给国际天文界、航天界等多个领域带来重大科学和技术挑战。21世纪以来，国际社会高度重视小行星撞击风险防御工作，联合国成立多个机构专门负责协调应对近地小行星撞击风险，各主要航天国家也积极致力于相关工作，并有力牵引、推动了空间科学研究、空间技术探索和空间资源利用。我国近地小行星撞击风险应对工作起步较晚、基础薄弱、国际贡献率低、国际话语权小，这与我国航天大国地位不符[8]。近几年我国已经开始开展近地小行星防御问题研究，受到国家高度重视，并已上升为国家战略[7]，且在小行星的轨道预测、领航监测器的轨道设计与控制、小行星的在轨处置和监测预警等方面取得了一定的研究成果[11-24]。应对小行星撞击风险主要涉及两个方面：(1)监测预警。搜索发现、跟踪定轨和数据更新、物性测量、计算撞击概率、预估撞击风险走廊和落点，预报撞击风险。即发现小行星和确定撞击风险。(2)在轨处置。依据监测预警信息，发射防御航天器，改变小行星轨道避免其撞击地球(即偏转，Deflection)，或将小行星分裂成碎片避免或降低其对地球的危害(即摧毁，Disruption)。这里无论是小行星轨道精确预报、在大天体摄动影响下的轨道不确定性，还是在人为外力作用下(如动能撞击等瞬时作用，或拖曳等长期作用下)的小行星轨道偏移、防御航天器自身的深空自主导航与控制技术等，都离不开对小行星运行轨道的精确计算测定，离不开精确的轨道设计与动力学研究。小行星的轨道动力学问题是防御的基础、前提和核心问题。精确的轨道设计与计算需要采用基于约束系统的几何动力学方法进行建模，同时采用高精度的几何数值积分方法进行数值模拟。而轨道动力学一直是分析力学发展进程中的核心问题之一，但它针对航天工程提出的更复杂的轨道动力学问题研究的不够深入，特别是航天工程中需要考虑到小行星的非规则形状和自转、引力源非球对称性对牛顿引力场的修正、太阳引力场的相对论修正、小行星运行环境的随机性等问题，给小行星轨道动力学建模、轨道高精度计算以及在轨处置航天器(群)的控制提出了更高要求，也是分析力学面临的新课题。

2 几何力学与几何数值积分的发展现状

现代分析力学最具标志性的成果是几何力学的发展，以及在此基础上构造的保结构算法——几何数值积分方法，将二者与控制结合起来所兴起的计算几何力学与控制使得若干动力学系统实现了长时间、高精度的数值建模分析，也将为小行星运行轨道的高精度数值模拟与测量、航天器(群)的精确轨道设计、动力学建模与控制提供一种理论基础。同时，小行星轨道动力学建模和高精度计算也必将进一步完善分析力学的高精度算法，特别是Hamel场变分积分子算法。

几何力学是基于流形、纤维丛、Lie群与Lie代数等现代微分几何工具，将经典分析力学升级成无坐标语言的整体表述形式，形式更加简洁优美，为多体系统、流体、弹性力学、场论、量子系统及几何控制理论等提供了统一的框架，它在空间探索的轨道设计、机器人运动控制、等离子体、超流、洋流动力学和气候建模、计算解剖学等诸多领域已经得到成功的应用[25,26]。利用现代微分几何语言表述Lagrange力学和Hamilton力学有以下几点优势：首先，适用于位形流形上的一般力学系统，以及由纤维丛上联络所刻画的约束力学系统[27,28]。联络可以自然地描述由非完整约束所决定的不可积分布，通过它将系统的基底变量和纤维变量联系起来，并具有鲜明的物理意义，进而采用规范理论的观点统一且有力地处理非完整力学和非线性控制理论中的诸多问题，包括力学系统的几何表示和几何结构、稳定性和镇定、最优控制等，在车辆、机器人等工程领域有重要的应用；其次，方便用于描述系统的大范围运动，也便于进行几何类比而将有限维连续系统的结构推广至无穷维和离散力学系统[29,30]。特别地，对于位形空间是Lie群情况，采用对称性约化理论[31]可实现代数空间上的线性化，还可以得到判别稳定性的能量—动量方法。最后，基于离散Hamilton原理的变分积分子构造保持力学系统几何结构的数值格式——几何数值积分算法，为获得长时间高精度的动力系统模拟提供有效的方法[32-36]。将计算几何力学同几何控制理论结合，得到能长时间保持高精度的数值控制算法，这一研究将为精确、高效、实时地控制航天器、潜水器、各类机器人等现代工程系统提供合理的设计原理。

目前，在无穷维系统中的几何力学和几何数值积分算法是国内外相关领域的研究热点，这既是分析力学研究的逻辑延伸，也是若干工程科学的实际需要，这其中Hamel形式及Hamel场变分积分子在动力学系统建模和高精度计算中发挥核心作用。采用无坐标的几何语言，可以通过类比，将具有对称性的有限维力学系统的成熟方法推广到无穷维系统，特别是借助于几何类比，可将分析力学中不依赖于坐标的活动标架法推广到无穷维系统和场论情况，通过直接对非物质速度进行变分的方式构造Hamel场变分积分子[36,37]，这个算法的能量保持显著优于Shake和Rattle算法的数值特性[38]；对基于SE(3)的几何精确梁进行动力学分析，验证了该算法能够长时间保持几何精确梁的能量与动量[39,40]；从协变场论观点与时空离散角度研究Hamel场变分积分子算法的保动量性质，在证明该算法保几何结构性质的同时，间接说明该方法具有长时间数值稳定的性质[41]。此外，Hamel场变分积分子在计算无穷维非完整力学系统时具有明显优势[42]，该方法表述Chaplygin-Timoshenko雪橇模型时表述简单，便于分析该动力学方程解的适定性和解释其特有的非线性行为。Hamel场变分积分子已用在工业机械臂、球面上移动机器人等的最优控制，展现出了较好的数值特性[39,43,44]。

尽管Hamel场变分积分子算法在无穷维系统的动力学建模和高精度计算中得到了验证，并且有理由相信不依赖于坐标选取的Hamel形式适于描述不规则形状和自转的小行星轨道动力学和姿-轨耦合问题，但针对小行星的精确轨道计算和预测，还需要考虑各种实际非理想因素，如小行星不规则形状、后牛顿效应的摄动等对构造Hamel变分积分子的影响。此外，小行星的运行环境是动态变化的，甚至是随机的，尽管在有关随机微分方程的处理上，几何力学和几何数值积分方法具有明显的优势[45]，但是Hamel场变分积分子对随机问题和基于学习的离散场论保结构算法的拓展还缺乏研究。

3 近地小行星轨道动力学高精度建模研究进展

小行星在太阳系中的运动非常复杂，存在各种摄动因素的影响，如各大行星和小行星带上质量较大小行星的质点引力摄动、太阳和大行星扁率产生的引力摄动、太阳引力场的相对论效应、小行星自转和形状不规则以及表面热物理环境的差异导致的亚尔科夫斯基(Yarkovsky)效应所产生的非引力摄动项等。因此，小行星的运动不是简单的Kepler问题，要精确计算和预测小行星运行轨道，一方面需要考虑诸多影响小行星运行的摄动因素，包括小行星在中心天体(太阳)和第三天体(大行星)的不规则形状造成的引力摄动[46-49]。为精确计算小行星轨道，并考虑大天体(八大行星和月球)与主带67颗质量较大的小行星的质点引力摄动、地球扁率项的引力摄动、太阳质心的后牛顿效应、太阳辐射压对小行星轨道的摄动，通过引入测量随机误差、利用运动方程的积分计算、观测量的理论归算等求解小天体运行轨道[12, 14]。文献[14]考虑各大行星的引力摄动、后牛顿效应等因素，采用改进的显式辛算法和嵌套的RKF7(8)积分器对43颗已命名的近地小行星轨道演化进行数值研究，在一定程度上弥补了对于大偏心率问题采用定步长的不足。除此之外，由于小行星自转、形状不规则以及表面热物理环境的差异导致了Yarkovsky效应也要产生非引力摄动项。另一方面在计算方法上，构造包含非保守力项的高阶Galerkin变分积分子能够处理小天体姿-轨耦合、太阳辐射力等非保守力项的情况[50]，而使用Chebyshev控制点与Lagrange插值方法构造的高阶谱变分积分子比常规的Galerkin变分积分子具有更好的稳定性[51]。采用变步长的数值算法，构造精细的力学模型并减少计算机截断误差，以达到小行星轨道计算精度和效率要求[11, 13]。虽然以上成果对影响小行星运行轨道的各种摄动因素有所研究，但是针对不同的小行星，对这些引力摄动和非引力效应如何对比观测数据进行系统考虑，以提高计算精度和效率，还缺乏深入研究。现有的小行星运行轨道的计算方法难以在动态环境下实现长时间高精度的效果，缺乏对复杂摄动因素的系统分析。此外，小行星运行环境瞬息万变，各种摄动要素的影响交替变化，具有一定的随机性。为此需要采用分析力学的几何数值积分方法，特别是Hamel场变分积分子算法，确保在引力和非引力摄动效应下的长时间高精度计算需要加强几何数值积分方法与随机问题的结合研究。

4 近地小行星探测/在轨处置航天器(群)轨道动力学与控制研究进展

近地小行星具有弱引力场、形状不规则的特点，从而使其外部引力场环境极其复杂，从而导致其运动形态复杂多样。这些复杂性导致轨道优化算法对初始值极其敏感，为探测/在轨处置航天器(群)的轨道设计和控制带来了极大困难和挑战[52-54]。

4.1 小行星处置航天器轨道设计与最优控研究进展

目前，近地小行星处置航天器的轨道设计通常利用Lagrange点L1和L2 周围的周期轨道的稳定和不稳定流形，揭示相空间中的管道网络，这解释了圆型限制性三体问题(CR3BP)中的迁移现象[55-59]。从不变流形理论出发，利用行星际高速公路来设计返回的转移轨道[15]，基于三体问题的不变流形设计低成本登月轨道[16,17]。最近发现了太阳系中的一个从小行星带延伸到天王星及更远地方的流形管道结构[60]；利用Lagrange拟序结构(LCS)理论[61]计算椭圆型限制性三体问题(ER3BP)中的迁移势垒，证明了ER3BP相空间中存在周期性脉动形式的Lagrange拟序结构，同时限制性四体问题中也存在依赖时间的不变流形[19-20]。将不变流形理论和离散力学与最优控制方法(DMOC)结合，成功计算出从地球到月球的具有脉冲推力的优化轨迹[62]。文献[18]综述了利用不变流形理论和弱稳定边界理论研究三体系统中低能量转移与捕获轨道设计的研究进展。在航天器群重构制导方面，基于笛卡尔坐标的重构技术，采用随机树搜索、遗传算法、非线性优化技术、序列二次规划、相关状态Riccati方程(SDRE)技术和脑风暴优化算法等，计算了跟随航天器从给定初态到给定末态的最省能量控制轨迹[63-68]。

4.2 伴飞航天器轨道动力学的综合建模研究进展

考虑伴飞航天器受太阳光压、其它星体引力摄动、随机扰动等综合影响的高精度轨道动力学模型具有重要意义。已有工作主要考虑类地J2摄动而忽略其它的引力项对小行星弱引力场中的摄动轨道建立引力摄动模型[69]，也有仅考虑太阳光压效应而不考虑其它因素，建立绕小行星的轨道摄动模型，用于小行星航天器群编队控制开发[70]。为实现从密切状态到平均状态的转换，现有的转换公式往往也是需要类地的关于引力场J2项的假设, 利用一个基于J2幂展开的Hamilton量来转换密切根数和平均根数，需假设J2比余纬向势大几个数量级，有的工作甚至忽略了除J2之外的任何摄动项[71]；对于中心天体的旋转速率等于或大于轨道平均运动，使用频率分解从相关信号中剔除密切(短周期)频率的做法计算代价大，且需频繁测量来检测短周期频率[72]。另一种思路是设计消除密切效应的滤波器，利用Kalman滤波器，考虑二阶引力势从平均绝对根数中重构密切绝对根数[73]。

4.3 探测/在轨处置航天器群编队制导与控制研究进展

近地小行星防御需要航天器群在深空轨道高精度编队飞行，与绕地轨道的编队场景(星链(图2)、GPS(图3))相比，具有更大的挑战性，描述其相对运动的方程更为复杂，必须考虑摄动影响，编队高精度控制的连续小推力往往控制在纳牛到微牛的量级[21]。要完成航天器群小行星防御任务需要解决两个基本问题：编队制导和编队控制，前者指集群轨道设计问题，后者为集群在轨保持或重构队形。在编队制导方面，常见的方法通常使用线性规划、迭代学习控制、线性二次调节器(LQR)、相关状态Riccati方程(SDRE)、Lagrange乘子、模型预测控制和感知框架等技术，通过跟踪预定义的参考点实现[74-81]。另外，基于相对轨道根数的动力学模型是将编队制导描述为跟踪特定相对轨道根数的形式[82-84]。在编队控制方面，可利用人工势场来实现重构[85,86]，也有使用多种生物优化技术来确定最优的重构方案[87,88]，还有提出一种基于相对轨道根数的解析式燃料最优脉冲编队重构策略[89]，以及基于相对轨道根数的ΔV最优方案和相关轨道动力学的精确线性化方法[90,91]。绕地球运行的集群编队控制方法，需要考虑大气阻力和J2对相对运动的影响[82-92]。为避免集群中个体的相对漂移、内部的碰撞，并减少整个任务中使用的推进剂，可以采用几种新的开环控制方法[91]。另外，编队保持方法也可针对特定轨道根数来设置[93,94]。

图2 美国星链计划[95]

图3 GPS星座[96]

5 问题与结论

小行星防御涉及数学、物理、化学、力学和天体力学等多个学科领域，是一个复杂的系统工程，但其防御离不开四步：一来、二去、三绕、四控。小行星防御过程如图4所示。

图4 小行星防御过程示意图

首先，所谓的“一来”，即如何对近地小行星轨道进行精确计算与预测问题。尽管国内外对影响小行星运行轨道的各种摄动因素有所研究，包括各大行星和小行星带上质量较大小行星的质点引力摄动、太阳和大行星扁率产生的引力摄动、太阳引力场的相对论效应、小行星自转和形状不规则以及表面热物理环境的差异导致的Yarkovsky效应所产生的非引力摄动项等，但是针对不同的小行星，对这些引力摄动和非引力效应如何对比观测数据进行系统考虑，以提高计算精度和效率，还缺乏深入研究。同时，现有的小行星运行轨道的计算方法难以在动态环境下实现长时间高精度的效果，影响了复杂摄动因素的系统分析，为此需要采用分析力学的几何数值积分方法，特别是Hamel场变分积分子算法，确保在引力和非引力摄动效应下的长时间高精度计算。此外，小行星运行环境瞬息万变，各种摄动要素的影响交替变化，具有一定的随机性，需要加强几何数值积分方法与随机问题的结合研究。

其次，所谓的“二去”，即如何设计发射近地小行星处置航天器的轨道与实现最优控制问题。根据时变系统的Lagrange拟序结构理论，通过有限时间Lyapunov指数(FTLE)可以找到时变系统中的分界线，以此可以区分流场中有不同动力学特征的区域。然而，关于小行星防御航天器的轨道设计方面，多采用不变流形理论来设计轨道，缺乏采用Lagrange拟序结构为最优控制算法提供初始值的研究。另一方面，航天器群制导未利用领航器的优化轨道，且基于笛卡尔坐标的重构技术计算量很大，此外，采用DMOC方法计算编队飞行卫星群重构的最优开环控制律，利用Halo轨道，限制了对小行星轨道的应用。

再次，所谓的“三绕”，即如何对近地小行星伴飞航天器进行综合建模问题。在构建绕小行星的受扰Kepler轨道模型时，主要考虑类地球小行星的J2摄动建立引力摄动模型，无法达到任意小行星所需精度水平，也缺乏考虑相对运动中太阳光压、其它第三星体引力摄动、随机扰动等综合影响，因而这样的轨道动力学模型精度不高。此外，天体力学N体问题中存在一种规范对称性，忽略其规范自由度可能造成轨道演化动力学在数值上的不稳定性[72]。由于现有模型主要考虑类地球系统的引力场J2项，密切状态到平均状态转换方法计算效率和计算精度较低，不适用于任意非地球小行星和其它行星体。有的剔除密切(短周期)频率的频率分解法具有计算代价大且需频繁检测的缺点，也不适于旋转速度慢的小行星。即使采用Kalman滤波器，仅考虑低阶引力场并不能准确得到小行星中心轨道的密切或平均轨道根数。

最后，所谓的“四控”，即如何对探测/处置航天器群实现有效的编队制导与控制，从而对具有潜在威胁的小行星实行干预或清除问题。目前，关于小行星防御航天器群队形重构和保持研究的多数技术是简单地推广近地编队飞行的方法和成果，不能满足集群任务关键要求。而且，航天器群制导未利用领航器的优化轨道，且基于笛卡尔坐标的重构技术计算量很大，缺乏计算量少、能有效利用ΔV、基于相对轨道根数的大规模编队重构控制方案。当前的一些方法只能考虑一个轨道周期的短时间跨度，动力学的线性化只在极短的时间内是准确的，而频繁的线性化又会增加计算成本。基于相对轨道根数的编队重构控制方案通常局限于两个航天器，并不能解决大群体的避碰问题。另外，缺乏解决在小行星特有的动态环境下航天器群全姿态最优控制方法，忽略航天器姿态可能会影响任务的执行。

综上所述，目前针对小行星防御问题的轨道预测、偏移与设计、航天器群的编队最优制导与控制的研究还不完善，尚缺乏综合考虑各种摄动效应的精确的数值计算方法和合理的小行星动力学模型，以及带随机扰动的在轨航天器(群)全位姿动力学模型；同时，围绕地球以外小行星特有动态环境下的航天器群编队最优制导与控制方法仍需进一步发展。上述问题的解决可为小行星的轨道偏转、天基监测系统的建立、引力拖车等编队任务提供理论指导和技术支持。鉴于此，采用基于几何力学和几何数值积分建立动态环境下的高精度Hamel变分积分子，通过高精度积分子研究长期作用力下小行星轨道的演化规律，对小行星的轨道偏转问题进行动力学建模，实现移步易景，实时计算小行星的运动状态，做到对小行星偏转轨道的精准预测，给小行星的灾害防御提供精准的动态数据。在此基础上，构建处置航天器基于Lagrange拟序结构的最优控制算法；建立考虑各种摄动的小行星在轨航天器全位姿动力学模型；针对航天器群编队的构建、保持和重构等任务，研究其场论框架下柔性集群编队控制理论和算法，为小行星防御的在轨处置提供理论支撑。