小行星动能撞击自主导航与制导方法研究

黄翔宇，徐 超，胡荣海，郭敏文

（1.北京控制工程研究所，北京 100094；2.空间智能控制技术重点实验室，北京 100094）

引言

太阳系中存在着大量的小天体，距离地球较近的小天体可能撞击地球，从而引发巨大的灾难。自20世纪90年代以来，近地天体撞击地球的潜在威胁受到了关注，各国开始实施小行星防御活动。在近地小行星撞击防御方面，提出了核爆、动能撞击、引力拖车等多种技术方案[1]。其中动能撞击是最为成熟的技术手段。由美国实施的“双小行星重定向测试”（Double Asteroid Redirection Test，DART）任务，以撞击双小行星系统中的一颗小行星，目前正在开展小行星动能撞击在轨技术验证[2]。在“十四五”期间，基于深空探测工程基础，中国也将论证实施近地小行星动能撞击在轨处置技术验证任务。

与以往小行星探测任务不同的是，小行星撞击任务是通过相对速度很大的撞击器直接撞击小行星以改变小行星的动量，其一般只适用于小尺寸的小行星防御[3]，DART任务目标小行星直径160 m，其对撞击精度要求极高，且在撞击段可用的小行星测量手段有限，一般仅有光学相机。基于光学图像的小行星探测自主导航方法一般可分为基于小行星中心视线测量和基于小行星表面陆标图像测量的方法。Chavez等[4-5]给出了基于光学陆标图像的小行星探测自主导航方法，但该方法只适用于距离小行星较近时的导航任务；采用小行星中心视线测量的自主导航方法由于缺少距离测量信息，通常需要辅以其他测量才能实现探测器状态的完全可观。Kubotaa等[6]给出了小行星中心视线测量结合测距信息的接近段自主导航方法；Takao等[7]给出了小行星尺寸信息辅助中心视线测量的导航方法；Jia等[8-9]给出了探测器间距离测量辅助的小行星接近段自主导航方法。这些方法均针对小行星抵近探测任务，其特点是利用光学相机得到的视线测量结合其他敏感器测量以视线探测器相对小行星位置和速度的完全估计。小行星撞击任务，由于相对速度可达10 km/s，采用测距敏感器获取距离测量信息的方法并不适用，主要是受限于测距敏感器的测量范围有限，可获取测距信息的次数极少，而且由于撞击任务中小行星尺寸小，在较远距离时要想获得距离测量需要保证极高的测距指向精度，工程实现难度大。同样小行星尺寸小，在撞击段的大部分时间里小行星在图像中均为点目标，而且通常缺少小行星精确的尺寸信息，利用小行星尺寸信息辅助的方法进行相对导航对小行星撞击任务也意义不大。对于采用探测器间距离测量辅助的方法需要同时发射多个探测器才能实施，并不是一个通用的方法，本文中并不考虑该方法。对于小行星撞击任务，本文提出仅利用小行星中心视线测量进行相对导航，将撞击器相对小行星的位置和速度中不可观方向进行分离，利用其可精确估计的位置和速度分量即可保证撞击精度，基于该方案撞击器以10 km/s的相对速度实现了50 m级近地小行星10 m级的撞击精度。

1 小行星撞击段GNC任务分析与方案设计

小行星撞击段定义为距小行星3万km～ 0 m的任务飞行段。针对50 m级近地小行星10 m级精度撞击任务，小行星动能撞击器在撞击段的自主导航、制导与控制（Guidance Navigation &Control，GNC）系统需要具备的功能要求：高精度的小行星观测、成像和中心提取、高精度的相对导航以及制导与控制。

考虑远距离50 m级小行星观测的需求，选择窄视场相机、星敏和惯性测量单元（Inertial Measurement Unit，IMU）作为撞击段的主要导航敏感器，其中窄视场相机需具备在3万km探测视15等星的能力，星敏提供撞击器的姿态基准，IMU在轨道机动时提供角速度和加速度测量。

对于小行星撞击段，窄视场相机测量的小行星中心视线方向是唯一可用的相对于小行星的测量。针对高精度小行星中心视线测量的需求，设计了远距离暗弱目标图像处理和提取方法以及图像面目标高精度中心提取方法，实现了小行星中心视线优于3角秒的提取精度。针对撞击段高精度自主相对导航的需求，设计了一种基于小行星中心视线测量的不完全可观相对自主导航方法，可观性分析表明该方法可对垂直于视线方向的相对位置和速度进行精确估计，从而可满足精确撞击制导任务需求。

在基于小行星中心视线的自主相对导航基础上给出一种迭代预测制导方法。考虑到在撞击段开始后需要尽快消除初始的速度偏差，以及为减小最后一次轨道修正所需的速度增量，在最后一次轨道修正前还应增加一次轨道修正，因此方案选择进行3次轨道修正。综合考虑轨道修正执行所需时间约束以及总的速度增量大小约束，方案设计分别在距小行星24 000、6 000、300 km处进行轨道修正，撞击目标位置设为视线方向与小行星表面交点处，最终实现10 m级的撞击精度。

2 小行星撞击段自主导航

2.1 小行星中心提取算法

图像处理在小天体动能撞击任务中具有非常重要的作用，是视觉导航的前提，其处理结果直接影响撞击器的状态估计及后续的轨道机动。然而，从首次捕获小天体的视觉图像到高速接近并对其进行撞击的整个过程当中，目标在导航相机中的像素面积、亮度大小和细节信息的丰富程度存在巨大的差异，故需要分别使用不同的图像处理技术来获取对应的导航信息，撞击器与小行星相距2 000～300 km的情况下，直径50 m的小行星在导航相机（视场角：1°，成像像元阵的大小：2 048 × 2 048）中的图像，如图1所示。

图1 不同距离下，导航相机拍摄的小行星图像Fig.1 Asteroid images captured by navigation camera at different distances

当小天体与撞击器距离较远（大于1 000 km）时，其在图像中表现为点目标，约占几个像素。此时，通过提取目标在图像中的一阶矩（即亮度中心）来粗略估计其中心的视线方向，然后通过解析函数拟合法估计其像素坐标。当撞击器不断接近小天体，小天体在图像中逐渐由点目标拓展为面目标（相对距离小于1 000 km，撞击前100 s以内）。值得注意的是，受太阳相位角（太阳-目标小行星-导航相机之间的角度）的影响，小行星在图像中的亮度分布将随着其自旋而发生周期性的变化，且相对距离越近，这种特性越明显。此时，通过计算目标在图像中的一阶矩来估计其中心将会造成较大误差，需要采用相位角修正。

令导航相机在k时刻观测到的图像为Ik。首先，利用中值滤波器去除图像Ik中散点噪声。然后，设置亮度阈值T以分割Ik中的目标小行星，得到二值图像Bk为

则图像Bk中非零元素所构成的集合Tk即为目标的像素区域，计算Tk在图像Bk中的一阶矩M以粗略地估计目标的中心[10]

其中：∑Tk表示计算Tk中元素的个数；分别表示Tk中第i个元素在目标二值图像Bk中的x坐标和y坐标。

当目标在图像中为几个像素大小的点目标时，利用点扩散函数（Point Spread Function，PSF）做进一步处理，对其亮度进行建模，并通过迭代的方法计算精度更高的目标中心。PSF取二维的高斯函数[11-12]

其中：s为亮度比例因子；[x0y0]为G（x，y）的峰值坐标；σ为标准差；b为背景亮度（pixels）。定义x=[s x0y0σb]T为待求参数。

令图像Ik中坐标为[u v]的点所对应的像素亮度为I(u,v)，即观测量。而通过在一个像素宽度内对式（3）积分可以估计出图像中坐标为[u v]的像素点所对应的亮度值

则观测残差（代价函数）为

推导e（u，v）相对于待求参数x的偏导数，可以得到其雅可比矩阵，然后通过最优化方法进行迭代求解。

为确定优化过程的初值，令x0=mx，y0=my，b为Ik中非目标区域的平均亮度值，s和σ由Tk的等效半径及M处的亮度值共同确定。

当小行星在图像中扩展为面目标且仍无法提取出清晰的轮廓时，解析函数拟合法将不再适用。同时，由于太阳相位角的影响，目标在图像中的一阶矩与其真实的质心位置存在很大的偏差，需要对一阶矩进行修正。

根据文献[13]，对于太阳相位角α以及Tk的等效半径R，偏移因子γ可通过下式计算得到

须指出，式（6）适用于球形天体。然而，当不规则的小行星在距离较远的时候，仍然可用于修正其亮度中心。修正后的坐标如下

其中：θ为太阳光在图像平面内的投影与图像x轴之间的夹角。

2.2 基于小行星中心视线的自主相对导航方法

在日心J2000惯性系下建立小行星的动力学方程为

其中：rA和vA分别表示小行星在日心J2000惯性系下的位置和速度；µs和µE分别表示太阳引力常数和地球引力常数；rE为地球在日心J2000惯性系下的位置；rA=‖rA‖，rE=‖rE‖；nA表示小行星动力学模型误差项，假设为高斯白噪声。本文小行星选为近地小行星，动力学模型只考虑太阳引力和地球引力摄动项影响，不考虑其他天体摄动及太阳光压等的影响。对于小行星动能撞击段来说，撞击器相对小行星速度快，撞击段持续时间短，因此上述假设是合理的。

在日心J2000惯性系下建立撞击器的动力学方程为

其中：rsc和vsc分别表示撞击器在日心J2000惯性系下的位置和速度，rsc=‖rsc‖；aSRP表示太阳光压引起的加速度；假设撞击器飞行过程中采用速度脉冲的形式进行轨道修正，对应加速度建模为狄拉克函数形式Δviδ(t-τi)[14]，Δvi为第i次轨道修正的速度增量，τi为第i次轨道修正的时间；nsc表示撞击器动力学模型误差项，假设为高斯白噪声。

对于小行星动能撞击任务来说，由于撞击段时间持续短（小于1 h），撞击器近似直线飞向小行星，因此为简化星上计算过程，导航系统中撞击器相对于小行星的动力学方程可近似描述为

其中：rscA和vscA分别表示撞击器在小行星J2000惯性系下的位置和速度。

在撞击器动能撞击飞行过程中利用2.1节中方法可提取小行星中心在像平面的位置（uc，vc），建立的观测模型为

其中：f为相机焦距。

其中：CCI为J2000惯性系到相机系的转换阵，其可通过星敏感器及相机安装等确定。

将式（10）给出的相对动力学方程作为导航系统的状态方程，将式（11）给出的小行星中心视线测量作为导航系统的观测方程，采用扩展卡尔曼滤波（Extended Kalman Filter，EKF）导航滤波方法即可对撞击器状态进行估计。

2.3 可观性分析

根据文献[15]可得到基于小行星视线矢量导航系统的可观测矩阵

式中：Hk和Φk,1分别为观测阵和误差状态转移阵。

其中：rscA1和vscA1分别为时刻撞击器在小行星J2000惯性系下的位置和速度。由式(16)可知MkN=02k×1，即为系统的不可观方向。

当k=2时

由式（1 4）可知Hc1和Hc2的秩均为2，即rank（Hc1）=rank（Hc2）=2，可得rank（M2）=4，则可观测矩阵M2零空间的维数为2，即对于两次观测情况下导航系统存在两个不可观方向[16]。由式(16)可知为 零空间的一组基底。令易知M2N2=04×1，显然N1，N2线性无关，故M2零空间可由N1，N2扩张而成，即对于两次观测情况下导航系统的两个不可观方向分别为N1，N2。

综上可知基于小行星视线矢量导航系统，当速度方向与相对小行星位置方向不一致时有且仅有一个不可观方向，即初始位置和速度不可观；当速度方向与位置方向一致时存在两个不可观方向，即小行星视线方向的位置和速度不可观。对于小行星动能撞击任务的撞击段来说，基于小行星视线矢量导航系统撞击器视线方向的位置和速度不可观。需要注意的是，对于撞击任务首先需要保证的是成功撞击到小行星预定位置，因此只要垂直于视线方向的位置和速度可观即可满足撞击任务控制需求，在撞击段自主导航过程中也只需利用小行星视线测量对垂直于视线方向的位置和速度进行导航修正。

3 小行星撞击段自主制导

设t时刻由光学相机测量得到小行星中心在相平面的位置为（ucAt,vcAt），则在小行星J2000惯性系下撞击器到小行星中心的视线矢量lAt可表示为

其中：R为小行星半径，本文中取25 m。设t时刻导航估计的撞击器在日心J2000惯性系下的位置和速度分别为估计的小行星位置和速度分别为则预测飞行时间可计算为

利用迭代制导计算t时刻轨道修正的速度增量Δvt为

利用式（20）通过多次迭代即可得到满足撞击精度要求的轨道修正速度增量。

一般来说，制导控制的目标为精确的撞击位置及撞击时间，而要实现这一目标需要导航系统能够精确给出撞击器的位置估计，考虑到小行星撞击段撞击器近似沿视线方向飞向小行星，在本节的制导方法中垂直于视线方向上的位置和速度估计决定了撞击的位置控制精度，沿视线方向上的位置和速度估计则决定了撞击时间的控制精度。由可观性分析结果可知，基于视线矢量的小行星撞击段自主相对导航系统，垂直于视线方向上的位置和速度是可被精确估计的，但沿视线方向的位置和速度不可观，导航系统未对其进行修正，制导控制过程中采用的沿视线方向上的位置和速度主要依赖先验信息。对于小行星撞击任务，重点关注撞击精度，而且沿视线方向也即撞击器飞行方向可通过无线电测控等手段获取较为精确的先验信息，因此利用本文提出的基于视线矢量的自主相对导航方法及在视线方向上位置和速度存在误差情况下的迭代制导方法能够满足小行星撞击任务的需求。

4 仿真校验

本文仿真初始数据如下：撞击器在小行星轨道系下的初始位置[30 000 × sin20°，0，30 000 × cos20°] km，初始速度[-10 × sin20°，0.001，-10 × cos20°] km/s；小行星在日心J2000惯性系下的初始位置[-80 768 079.149，-137 382 451.608，2 507 154.394] km，初始速度[24.885 9，-12.070 2，-3.779 1] km/s；鉴于沿撞击器飞行方向上地面定轨精度较高，设置初始导航误差：垂直于视线方向上位置随机误差100 km（1σ）；速度随机误差各方向10 m/s（1σ）；沿视线方向的位置随机误差10 km（1σ）；速度随机误差各方向1m/s（1σ）；导航相机视场1°；图像分辨率2 048 × 2 048；轨道修正速度增量控制方向误差0.55°（3σ），大小误差3%。

在撞击器接近小行星的过程当中，不同距离下目标中心的部分提取结果如图2所示，目标中心提取的测量误差如图3所示。

图2 不同距离下，目标小行星的中心提取结果Fig.2 Extracted centers of target asteroid under different distances

图3 不同距离下，目标小行星的中心提取误差Fig.3 Measurement errors of target center under different distances

当撞击器与小行星距离较远时（30 000～1 500 km），中心测量误差较小（0.2 pixels）并周期性地变化；因为小行星的形状不规则，同时存在自旋现象，导致其在图像当中的形状产生周期性变化，结果如图2（b）与图2（c）所示。基于该方案撞击器以10 km/s的相对速度实现了50 m级近地小行星10 m级的撞击精度。也会周期性地影响中心测量结果。随着相对距离的降低（1 500～300 km），u方向的像素测量噪声有增大的趋势；因为此时小行星在图像当中逐渐扩展为面目标，中心提取结果受太阳相位角的影响越来越大。末端（300 km）的中心提取误差为[0.446，0.073] pixels。

撞击段撞击器到小行星中心的距离变化曲线如图4所示。基于小行星视线测量的自主导航撞击器相对与小行星位置和速度估计误差在视线坐标系下的表示如图5和图6所示。其中视线坐标系定义：以式(17)中的lAt为视线坐标系zLOS轴指向，任取不平行于lAt的方向m，本文中取m=则视线坐标系yLOS=lAt×m，根据右手准则可得视线坐标系xLOS轴。

图4 撞击段撞击器到小行星距离Fig.4 Distance between impactor and asteroid during impact phase the impact phase

由图5和图6可知基于小行星视线测量的自主导航方法能够准确估计垂直与视线方向的位置和速度，对于沿视线方向上的位置和速度导航过程中未作修正。撞击器相对于小行星的位置和速度估计误差在J2000惯性系下表示的对比结果如图7和图8所示，由于受到视线方向位置误差的影响，在J2000惯性系下撞击器相对于小行星的三轴位置和速度估计均不准确。

图5 视线坐标系下撞击器相对于小行星的位置估计误差Fig.5 Estimation error of relative position between impactor and asteroid in LOS frame

图6 视线坐标系下撞击器相对于小行星的速度估计误差Fig.6 Estimation error of relative velocity between impactor and asteroid in LOS frame

图7 J2000惯性系下撞击器相对于小行星位置估计误差Fig.7 Estimation error of relative position between impactor and asteroid in J2000 inertial frame

图8 J2000惯性系下撞击器相对于小行星速度估计误差Fig.8 Estimation error of relative velocity between the impactor and the asteroid in J2000 inertial frame

根据本文给出的方案，在小行星撞击段分别在距小行星24 000、6 000和300 km处进行3次轨道修正，每次轨道修正的速度增量如图9所示，3次轨道修正的速度增量分别为1.21、0.23和1.48 m/s。在撞击段飞行过程中撞击器相对于小行星的速度方向与相对位置方向夹角如图10所示。

从图9和图10可知：轨道修正后相对速度与相对位置方向夹角均被修正到小于5×10-4°，且距小行星越近时撞击器相对于小行星的位置和速度方向夹角增大速率越快。因此，最终的撞击精度由最后一次轨道修正决定，但在最后一次轨道修正前，尤其是在距离小行星较近时增加轨道修正频次，可以极大减少轨道修正所需的速度增量。

图9 在J2000惯性系下撞击器轨道修正速度增量Fig.9 Velocity increment shown in J2000 inertial frame for impactor’s orbit corrections

图10 撞击器相对于小行星的速度方向与位置方向夹角Fig.10 Angle between relative position and relative velocity between impactor and asteroid

垂直于视线方向的撞击器相对位置估计误差500次打靶仿真结果如图11～13所示。其中在视线坐标系下垂直于视线方向上撞击器相对位置和速度估计误差随与小行星距离变化的曲线如图11和图12所示。500次打靶中实际撞击点与目标撞击点位置误差在垂直于视线平面上的散布如图13所示，可见撞击精度均优于4 m。

图11 500次打靶垂直于视线方向的撞击器相对位置估计误差Fig.11 Relative position errors in the plane perpendicular to the LOS in 500 Monte Carlo simulation runs

图12 500次打靶垂直于视线方向的撞击器相对速度估计误差Fig.12 The relative velocity errors in the plane perpendicular to the LOS in 500 Monte Carlo simulation runs

图13 500次打靶撞击点位置误差Fig.13 Impact site position errors in 500 Monte Carlo simulation runs

5 结论

针对50 m级近地小行星10 m级撞击任务，提出了一种小行星撞击段高精度自主GNC方案。针对高精度小行星中心视线测量的需求，提出了点目标和面目标的中心提取方法，首先利用矩算法获取初值，然后通过函数拟合以及相位修正来进一步提高精度，实现了中心视线的测量精度优于3″。针对撞击段高精度自主相对导航的需求，设计了一种基于小行星中心视线测量的相对自主导航方法，分析了仅利用视线测量的自主导航系统可观性，实现了对垂直于视线方向的相对位置和速度的精确估计，满足了精确撞击制导的要求。针对精确撞击制导的需求，基于相对导航的状态估计给出了一种迭代预测制导方法，分析了制导修正时机对轨道修正速度增量的影响，实现了10 m级撞击精度任务要求。