威胁小天体动能撞击防御误差演化分析*

张高辇，李翔宇

（北京理工大学宇航学院，北京 100081）

1 引言

近地小天体的轨道与地球轨道接近，部分轨道甚至与地球相交，存在潜在的撞击风险，是人类面临的潜在重大威胁之一。地球上曾经发生的22 次不同程度的生物灭绝事件至少有10 次以上是由小行星撞击地球引起的［1］，如2013 年俄罗斯车里雅宾斯克地区陨石坠落事件等。因此研究小行星撞击防御问题具有重要意义。近年来，针对小行星防御的研究逐渐成为航天领域研究的热点。目前现有提出的方法主要包括核爆［2］、动能撞击［3］、激光烧蚀［4］、引力拖曳［5-7］等。其中，动能撞击被认为是针对直径1 km 以下小行星防御的高效并且可行的方法［2］。2022 年，美国国家航空航天局（NASA）成功实施针对Didymos双小行星系统第二颗小行星（65803）的撞击实验DART，该实验将是采用动能撞击进行小行星防御的第一次验证，撞击结果将通过欧洲航天局（ESA）未来发射的LICIA 立方星观测得到［8］。对于威胁小行星动能撞击防御任务而言，动能撞击效果的主要评价指标是撞击后小行星近地点距离的增加量，而这一指标可近似等效于目标小行星在接近地球的B 平面上的偏转距离［9，10］。以此为基础Izzo［11］提出了撞击几何（IG，Impact Geometry）的概念用于评估撞击偏转距离，2019年，Davide Farnocchia 等［12］针对考虑误差下的行星接近分析的B 平面理论做了进一步研究。但是以往的动能撞击研究大多针对对心碰撞，动量传递效果仅通过一个系数K 决定［13］，或者仅考虑撞击后沿撞击点处溅射物对法向方向动量传递的影响，而忽略了溅射物切向方向动量的传递。清华大学的焦艺菲等［14］提出一种利用两个动量传递系数β和γ描述动能撞击的斜撞击模型，该模型结合撞击几何可实现最优斜撞击方向的快速求解，并且可应用于任意形状小行星最优撞击方向的确定。此外，中科院的王艺睿等［15］针对Apophis 偏转任务分别研究了考虑运载能力限制下进行撞击器轨道的优化设计，后续并在考虑地球借力和多脉冲机动模式下进行了改进［16］。

由于动能撞击一般在深空进行，而深空环境下撞击器的导航制导控制误差必然会带来动能撞击任务的偏差，并且小行星动能撞击防御的周期一般在10 年以上，动能撞击时刻撞击的微小误差可能会在时间的累积下产生巨大的影响，因此有必要针对威胁小天体动能撞击防御进行撞击后的误差演化分析，进而评估撞击任务设计的合理性与有效性。2010 年，Armellin R 等［17］利用微分代数方法分析了小行星Apophis 在2029 年飞掠地球时的近心点距离和近地点时刻偏差；2017年，Feldhacker 等［18］研究了动能撞击位置偏差和动量传递系数偏差对不同形状小行星撞击效果的影响，但是该文献没有分析动能撞击后偏差对小行星轨道分布的影响。状态转移张量（STT）方法是2006 年首先由Scheeres 等［19］提出的一种计算航天器误差演化的工具，该方法通过半解析地计算状态转移张量实现了误差预报的快速高效分析。杨震等［20，21］针对带有脉冲机动过程的误差演化方法进行了研究，推导出可以连续预报的状态转移张量。

本文针对威胁小天体动能撞击防御过程的误差演化问题，考虑由撞击器制导控制误差带来的撞击位置误差，根据斜撞击模型下的动量传递关系，计算出由撞击位置误差引起的小行星速度增量偏差，以及偏差的均值和方差的传递关系，然后，基于状态转移张量STT 方法，针对动能撞击带来的速度增量偏差演化到小行星抵达近地点时刻，得到小行星考虑误差情况下在地球附近的位置速度散布情况，结合蒙特卡洛方法进行仿真验证，结果表明STT 方法可以有效预报动能撞击过程的误差。

2 基于状态转移张量的小行星误差演化方法

2.1 动力学方程

下面给出航天器撞击小行星后的动力学模型，本文选择考虑八大行星和月球作为第三体引力摄动影响下的动力学方程，该方程为：

式中：dpi表示第i个行星的日心矢量，ppi表示第i个行星到小行星的矢量，µS和µpi分别表示太阳和第i个行星的引力常数，行星和月球的位置调用SPICE星历［22］。

2.2 状态转移张量的推导

将航天器的动力学写为张量形式：

式中：x={xi|i=1， …， n}，n=6 是动力学方程的维数。对于给定的初始状态x0=x（t0），因此x（t）可以看作初始状态x0的函数，进而定义隐式解为：

注意到，该隐式解满足：

可知，若获得了航天器运动状态的非线性映射关系φ，则从航天器初始状态x0到任意时刻状态x（t）的传播关系就可以解析计算。给定航天器初始参考状态及相对参考状态的状态偏差δx0，可将任意时刻航天器运动轨迹相对标称轨迹的偏差量δx（t）表示为

对其求导可得

则标称轨道与当前状态的偏差可以表示为：

式中：

同理，对状态量的导数进行Taylor展开可得

式中：

式中：i表示向量函数的第i个分量，M为泰勒级数展开的阶次，为局部动力学张量，沿标称状态计算，是从t0到t时刻的状态转移张量，这里运用了爱因斯坦求和约定。将绝对轨道偏差的一个具体样本δx（t）看作航天器实际状态x（t）相对标称状态x(t)的一个相对运动状态，式是一组非线性的相对运动方程，对公式求导可得

将公式展开并对比其与公式相同项的系数，可以得到初始相对状态到任意时刻相对状态的状态转移张量所需要满足的微分方程。状态转移张量可以通过积分微分方程得到，这样δx（t）就可以由x0解析表示，这种方法称为半解析方法。一般地，状态转移张量展开的次数越高，预报的精度越高，但是随之带来计算量将指数增加。文献［20］表明，对于航天器轨道的误差演化问题，采用2 阶模型就可以达到较高的精度。本文针对小行星状态的误差演化选取2 阶精度。

2.3 误差统计矩的分析

根据概率论的基本理论，自由向量δx的均值和协方差矩阵定义为：

其中E［·］代表期望算符。将式代入方程可以得出相应的均值和协方差演化结果，即为：

由公式可知，当各阶状态转移张量沿标称轨迹计算后，基于状态转移张量的各阶统计矩预报仅为解析的数学运算。当状态转移张量展开到M阶，预报均值需要计算初始偏差的前M阶中心距，预报协方差矩阵需要计算初始偏差的前2M阶中心距。因此，Taylor展开阶次越高、预报的统计矩阶次越高，计算公式越复杂，且计算量越大。

本文仅考虑前2 阶矩的非线性传播，预报协方差矩阵需要计算前4阶中心距，终端均值m(t0)与协方差矩阵P(t0)可简化为：

对于高斯分布，前4阶统计矩可以写为均值m和协方差P的表达式：

为了简便起见，均值和协方差误差演化结果可以简化为：

需要注意的是，在本文研究的问题中，目标威胁小行星在飞掠地球时可能处于地球影响球内部，此时小行星受力主要为地球引力，而影响球外的大部分轨道主要受太阳引力影响。因此在使用变步长积分器计算状态转移张量STT 时，需要注意在影响球内降低积分步长。

3 基于斜撞击的动能撞击过程分析

3.1 动能撞击任务设计原理

为了计算动能撞击对小行星防御的偏转效果，通常需要积分计算小行星整个任务过程中的非线性轨道力学方程。文献［11］表明，通过B 平面计算动能撞击造成小行星偏转距离的可达范围是有效的。B 平面是指过地心并且与小行星双曲线轨道渐近线相垂直的平面。即该平面是与进入双曲超速方向U 相垂直的平面，如图1 所示。相对速度U=|U|的大小定义如下［12］：

图1 小行星接近地球及其B平面定义Fig.1 Asteroid approaching the Earth and its B-plane definition

为了利用B 平面来描述小行星飞越地球时的轨迹，可以建立行星中心坐标系(ξ,η,ζ) ，如图1所示。其中η轴沿着U 矢量方向并且ζ轴指向地球速度矢量vE在B 平面投影的反向。ξ轴满足右手定则。

对于与地球存在碰撞可能的小行星，采用动能撞击方法需要最大化偏转小行星到达地球的B平面距离。文献［9］给出了近似的定量关系。

式中：K代表动量传递系数，一般被认为是给定值。aast表示小行星的半长轴，vE代表地球的速度大小，θaE表示小行星双曲线渐近线与地球速度方向夹角，μs表示太阳的引力常数，mast表示小行星的质量，ms/c表示航天器的质量，Vsa表示撞击航天器相对小行星的速度矢量，vast表示小行星的速度矢量。

通常将偏转距离Δζ作为动能撞击的优化指标，Δζ定义为地球和渐近线之间的距离。从上式可以看出为了最大化偏转距离|Δζ|，需要最大化内积Vsa·vast的绝对值，该值常被称为撞击几何。

3.2 斜撞击模型下的动量传递关系

当撞击器相对小行星的速度与小行星-撞击器连线不重合时，称为斜撞击。撞击速度vr一般表示为vsc-vast，即撞击器相对于小行星的速度［14］，其中vsc表示航天器速度，vast表示小行星速度。

图2 定义了小行星撞击坐标系Oxyz，该坐标系可以用来描述和计算不同撞击位置对小行星动量变化的影响。其中k 为vr反方向的单位矢量，i为撞击瞬时小行星角动量方向的单位矢量，j满足右手定则，i、j、k 分别为Ox、Oy 和Oz 轴的单位方向向量。由于撞击过程极短，本文忽略了小行星自转带来的影响。图3 定义了撞击点位置的极坐标角度参数θ和φ，其中θ为从i轴逆时针转到撞击位置矢量在xy平面投影所需角度，φ为撞击点位置矢量与k轴夹角。其中nˆ表示撞击点处切平面指向内部的法向单位矢量，表示单位矢量k 与单位矢量nˆ作2次叉乘得到的单位矢量。根据几何关系可以得到单位矢量nˆ和tˆ关于θ和φ的显式表达式为

图2 撞击坐标系定义Fig.2 Impact coordinate system definition

图3 小行星撞击坐标系定义Fig.3 Asteroid impact coordinate system definition

式中：0≤θ≤2π，为防止斜撞击角度过大而偏离小行星，一般限制φ∈ ［0， π /3］。

根据文献［14］中相关的定义，引入法向动量传递系数βf和切向动量传递系数γf，用以表征斜撞击带来的撞击点法向和切向动量传递效果。此时小行星的速度改变量可以表示为：

代入公式可以得到速度改变量的表达式：

撞击参数中，需要考虑的主要有撞击角度θimp以及两个动量传递系数β和γ。一般动量传递系数通过数值仿真软件得到，以光滑粒子流体动力学（SPH）方法为算法基础的软件可以作为求解手段［23］。

3.3 动能撞击后小行星的状态误差分析

由于动能撞击一般发生在深空，远离地球测控，动能撞击末制导段不可避免产生撞击误差，进而使得小行星的速度改变量发生偏差。本节将分析动能撞击过程中小行星的状态偏差传递过程，并给出动能撞击后小行星的状态误差模型。

本文主要关心动能撞击后小行星轨道的误差演化情况，整个动能撞击后小行星的运行轨迹流程如图4 所示。根据动能撞击的实际物理过程，本文将小行星撞击后的误差δx+0分为2 部分处理，即动能撞击前小行星自身的测定轨误差δx （t0），以及动能撞击位置偏差带来的小行星速度改变量的误差δxv0。

图4 动能撞击小行星示意图Fig.4 Diagram of kinetic impact on the asteroid

首先根据动能撞击的实际过程，分析动能撞击后小行星状态的偏差。定义动能撞击时刻为t0，动能撞击前一时刻t-0小行星的状态量误差为δx0，撞击带来的小行星速度脉冲的误差为δΔv0，则动能撞击后t+0时刻小行星的状态量偏差可以表示为：

其中R= [ 03×3,I3×3]T，δΔv0是小行星速度增量与标称速度增量之间的偏差，δxv0是由动能撞击偏差引起的小行星状态量的偏差，δx（ t0）是小行星在动能撞击前的定轨误差。

其次考虑动能撞击带来的小行星速度改变量的偏差。将由动能撞击带来的小行星速度改变当作脉冲增量，则可以记标称情况下小行星的速度增量为考虑小行星的真实脉冲误差为δΔv0，那么小行星真实的速度改变量Δv0为：

对于速度改变量的偏差δΔv0，参考文献［18］，本文仅考虑由航天器在垂直于相对小行星速度方向上的撞击位置偏差带来的小行星速度增量的偏差，如图5 所示。航天器撞击小行星不同位置时通过前文的斜撞击模型求解小行星的速度增量。本文假定撞击航天器在垂直于相对速度平面（即小行星撞击坐标系的x-y 平面）上撞击位置的分布为2 维高斯分布。对于球形的小行星，通过在垂直于相对速度平面求撞击速度增量Δv0相应的数学期望即可得到小行星速度改变量的均值m（Δv）和方差P （Δv）。具体公式如下：

图5 垂直相对速度方向撞击位置误差Fig.5 Impact position error perpendicular to relative velocity direction

其中，S表示球形小行星在x-y平面上的投影区域。注意到，当假设小行星为球形时存在几何关系：

其中（x，y）表示撞击位置在垂直相对速度平面上的坐标。小行星速度改变量Δv 的方差P （Δv）表达式为：

撞击位置的二维高斯分布概率密度函数为：

式中：x = [x, y]T，表示撞击位置在小行星撞击坐标系x-y 平面投影的坐标，并且其中均值m=m（x），方差P=P （x）。此时，本文给出了动能撞击后小行星位置速度的误差模型，后文将基于该模型进行数值仿真。

4 仿真校验

4.1 算例配置

假定目标小行星为2004年发现的威胁小行星Apophis，该小行星质量估计为4.0×1010kg，忽略小行星原有形状，假定其为等质量的均质圆球，假定小行星的密度为2600 kg/m3，可以计算得到小行星的等效半径约为154.288 m。Apophis小行星将于2029年4月13日飞掠地球，近地点距离将低于同步轨道卫星距离。考虑撞击器的质量为4500 kg，撞击处的相对速度大小为6 km/s，相对速度方向沿小行星速度方向，动量传递系数β和γ分别为2.5和0.5。假定动能撞击时刻为2024年4月13日21：46：07.59（TDB），关于小行星和撞击器状态量的设置详见表1，该状态是通过积分文献［16］中2019 年1 月1 日00：00：00（TDB）小行星Apophis的位置速度到动能撞击时刻得到。

表1 动能撞击时小行星与航天器日心赤道惯性系下状态量Table 1 States of the asteroid and spacecraft in the heliocentric equatorial inertial system at kinetic impact

仿真结束条件为小行星到达近地点，本文假定地球引力影响球的半径为925000 km。

在航天器与小行星动能撞击时刻，设置标称撞击点在小行星撞击坐标系下的角度θ=0°，角度φ=0°，从而得到日心赤道惯性系下小行星速度标称改变量为：［0.10678；1.57778；0.58899］mm/s。撞击前后小行星的轨道以及对应时刻地球轨道关系如图6、图7所示。

图6 初始小行星Apophis、撞击后小行星Apophis轨道和地球轨道Fig.6 Initial Apophis orbit， post-impact Apophis orbit and Earth orbit

图7 小行星Apophis到达地球影响球边界处的位置散布Fig.7 Distribution of the positions where Apophis reaching the boundary of the Earth influence sphere

4.2 仿真结果

针对小行星在与航天器动能撞击后的误差演化问题，本文分析了小行星在仅考虑撞击位置偏差带来的速度增量误差情况下，到达地球影响球边界t1时刻以及到达近地点t2时刻小行星位置速度的误差分布情况。仿真结果如下。

4.2.1 仅考虑撞击后小行星的速度增量偏差

对于仅考虑动能撞击造成小行星速度增量偏差的情况，如前所述，本文假定动量传递系数和撞击相对速度无偏差，仅考虑航天器撞击小行星表面不同位置带来的小行星速度增量偏差。设置撞击器在x-y 平面投影位置（x， y）的均值为mxy=[0, 0]T，方差为Pxy=diag（[5m, 5m]2）。

小行星Apophis 到达影响球边界时在日心赤道惯性系下的位置速度散布如图8 和图9所示。

图8 小行星到达影响球边界处的位置3σ椭球Fig.8 The 3σ ellipsoid of the positions where the asteroid reaching the boundary of the influence sphere

从上述图中可以看出，在本算例中仅考虑撞击位置偏差，演化5 年后，当小行星到达地球影响球边界时，其位置散布在10 km 量级，速度散布在0.1 cm/s量级。

如图10 和图11 所示，可以看出动能撞击前后小行星轨道在近心点附近的改变，从图11可以清晰看出动能撞击后小行星位置误差散布。对比图11和图7可以看出小行星在到达影响球边界和到达近地点时的位置偏差3σ散布差异不大，均在10 km左右。

图11 小行星到达近地点处的3σ椭球Fig.11 The 3σ ellipsoid at which the asteroid reaching its perigee

图12 小行星到达近地点处的位置3σ椭球Fig.12 The 3σ ellipsoid of positions where the asteroid reaching its perigee

从上述图中可以看出，在本算例中仅考虑撞击位置偏差，当小行星到达近地点时，其位置散布在10 km量级，速度散布在10 cm/s量级。但是对比图13和图9可以看出，小行星在到达近地点时的速度误差散布较在地球影响球边界处高2 个数量级。同时，上述结果表明，采用2 阶STT 可以对小行星的终端位置速度误差进行有效预报，预报结果与蒙特卡洛仿真结果吻合。

图13 小行星到达近地点处的速度3σ椭球Fig.13 The 3σ ellipsoid of velocities when the asteroid reaching its perigee

5 结论

本文基于STT 方法针对小行星动能撞击防御中的斜撞击进行误差演化分析。首先，基于斜撞击模型描述航天器撞击小行星过程，通过两个动量传递参数β和γ得到撞击后小行星的速度改变量；然后，针对考虑航天器撞击小行星表面的位置偏差，在斜撞击的模型基础上推导了撞击位置误差到撞击后速度增量误差之间的均值和方差传递关系；最后，针对撞击后小行星位置速度的偏差结果，通过状态转移张量STT 方法演化得到小行星进入地球影响球边界以及近地点处的位置速度误差。结果表明，动能撞击位置误差对小行星对动能撞击效果的影响较大，对于半径150 m 左右的小行星、小行星速度增量1～2 mm/s量级的动能撞击而言，3σ边界为15 m的撞击位置偏差将在5 年后使得小行星在近地点处产生10 公里左右的位置偏差，0.1 m/s量级的速度偏差散布。如果考虑到小行星质量估计的偏差、小行星材质带来的动量传递系数偏差等因素的影响，真实的动能撞击带来的位置速度散布将更大。因此在动能撞击任务设计的过程中需要考虑撞击航天器的制导控制精度的影响，避免撞击误差对撞击效果的不利影响。