一种未知形状小天体的撞击中心计算方法及系统

1.本发明属于导航信息处理领域，具体涉及一种未知形状小天体的撞击中心计算方法及系统。


背景技术：

2.小天体撞击中心计算是小天体防御撞击偏转光学导航的重要内容。现有的技术包括矩估计方法、轮廓扫描法、椭圆拟合法与互相关方法等，它们在计算复杂度、导航精度、适用天体类型、所需先验信息等方面各有特点。
3.矩估计方法简称矩方法，是最基本、也是最有效的中心提取算法。矩方法是利用矩算法对图像中目标天体的亮度中心进行估计。目标天体所占的像素范围通常由背景分割阈值决定，背景分割的目标是目标天体亮度中心的估计不被背景亮度以及其它的亮点所干扰而引起偏差。通过矩估计找到亮心以后，需要修正太阳相角所引起的亮心与天体中心的偏差效应。
4.轮廓扫描法是光学导航中常用的方法。首先要在图像中提取轮廓点，然后基于相机与小天体的相对位姿的先验信息与小天体模型表面的三维点进行对应。基于二维点与三维点的对应关系，可以构建pnp问题使用最小二乘法对图像中目标的中心位置进行估计更新。该方法需要小天体精确的三维模型先验信息。
5.椭圆拟合法适用于三维模型能用三轴椭球体良好描述的目标。它的基本思想是：当一个椭球体投影到二维图像平面上时，形成一个椭圆；这个椭圆可以用来估计相机相对于小天体的位姿。显然该类方法需要目标的先验信息，且前提是小天体需要能用三轴椭球体良好描述，但是对于未知的小尺寸小天体很难成立。
6.互相关方法通过预测的像元亮度值模板与真实图像的相关系数来确定目标中心，并可以通过拟合二维抛物面来进一步细化相关性峰值以确定表面的峰值亚像素位置。这种方法的应用很广泛，既可以用于静止光斑的中心计算，也可用于目标拖尾图像中心提取。但是这类方法需要基于小天体模型和成像场景的先验知识渲染小天体的预测视图来生成用于执行互相关的模板。
7.上述方法中只有矩估计方法是不需要先验模型的，但是矩估计方法在目标三维模型非常不规则，且非凸时，估计的中心位于三维模型的边缘，中心错失的风险较大。而对于轮廓扫描法、椭圆拟合法与互相关方法等其它方法，要么假设目标天体形状模型是三轴椭球体，要么假设目标天体形状的凸的半规则椭球体，要么假设目标是不规则的，但是这些模型都是精确已知的，由解析模型或三角片模型或球谐模型来表示。然而，在小天体防御的背景下以及随着越来越多的小尺寸小天体被发现的背景下，目标小天体的形状是非常不规则的，且是未知的，这让这些方法统统失效。因此，需要提出一种新的小天体撞击中心计算方法，利用该方法能够从安全性出发有效地提取目标小天体的图像中心，进而为小天体撞击提供瞄准中心，降低撞击中心错失风险。


技术实现要素：

8.本发明的主要目的是提供一种未知形状小天体的撞击中心计算方法，旨在解决现有同类方法依赖于小天体形状先验信息和撞击中心错失风险较高的技术问题。
9.为实现上述目的，本发明提出一种未知形状小天体的撞击中心计算方法，包括以下步骤：步骤s1，小天体的可见部分轮廓提取：基于捕获的目标未知形状小天体的图像，获得图像中小天体的可见部分轮廓；步骤s2，将所述小天体的可见部分轮廓进行照明真轮廓与照明伪轮廓分割：计算太阳光照方向到图像坐标系的y轴方向的旋转矩阵，通过计算的所述旋转矩阵将包含小天体轮廓的图像进行旋转，然后对旋转后图像的每一列像素从下至上进行扫描，并按照从左至右的方向对所有列像素进行逐列扫描；对于每一列像素，取第一个遇到的轮廓点作为照明真轮廓点，所有的照明真轮廓点构成所述照明真轮廓；如果有第二个及之后遇到的轮廓点，则将其作为照明伪轮廓点，所有的照明伪轮廓点构成所述照明伪轮廓；步骤s3，识别所述照明真轮廓的凹凸特性：对所述照明真轮廓进行旋转，然后对旋转后的照明真轮廓进行二次曲线拟合，根据所述二次曲线拟合的结果来判断所述照明真轮廓的凹凸特性；步骤s4，选择估计方法获取小天体的撞击中心：如果判断到所述照明真轮廓是凹的，则采用基于骨架/中轴线的中心估计方法获取小天体的撞击中心；如果判断到所述照明真轮廓是凸的，则采用矩估计方法获取小天体的撞击中心；步骤s5，根据步骤s4获取的撞击中心计算所述撞击中心的安全边界距离。
10.进一步的，所述步骤s1中，对捕获的图像依次进行图像平滑、二值化、闭运算处理，并进行外轮廓提取，获得图像中小天体的可见部分轮廓。
11.可选的，所述图像平滑采用高斯滤波器，所述二值化处理采用otsu阈值分割方法，所述外轮廓提取采用行扫描法。
12.进一步的，所述步骤s2中，所述图像坐标系为建立在相机像平面内的直角坐标系o
‑
xy，原点位于像平面的中心o，x轴指向像平面的水平方向，y轴指向像平面的竖直方向；所述太阳光照方向到图像坐标系y轴方向的旋转矩阵为：式中，sin表示正弦函数，cos表示余弦函数，表示太阳光照方向在图像中的投影单位矢量与图像坐标系x轴的夹角且表达式为：其中，是太阳光照方向在图像中的投影单位矢量在图像坐标系x轴方向的坐标分量，是太阳光照方向单位向量在图像坐标系y轴方向的坐标分量，是符号函数，arctan表示反正切函数。
13.进一步的，基于所述照明真轮廓的轮廓点拟合得到的二次曲线方程为：进一步的，基于所述照明真轮廓的轮廓点拟合得到的二次曲线方程为：，
则定义照明真轮廓拟合误差平方和为式中，分别是基于所述照明真轮廓的轮廓点拟合得到的二次曲线方程的常数项、一次项、二次项的系数，分别是照明真轮廓的轮廓点的x坐标与y坐标分量，下标，其中表示照明真轮廓的轮廓点数量；根据所述照明真轮廓拟合误差平方和，构建照明真轮廓对应的最优拟合二次曲线方程的多项式系数的计算公式如下其中，分别为照明真轮廓对应的最优拟合二次曲线方程的常数项、一次项、二次项的系数，arg min表示最优化问题取极小值点的数学符号；通过判断所述多项式系数中的符号确定照明真轮廓的凹凸性：如果，则照明真轮廓是凸的，否则照明真轮廓是凹的。
14.进一步的，所述步骤s5中，所述安全边界距离为小天体的撞击中心距离真实小天体轮廓的符号距离，安全边界距离根据下式进行计算：，式中，表示估计轮廓中心；与分别是在图像坐标系o
‑
xy下的横坐标与纵坐标；表示小天体的真实轮廓且，其中，分别为小天体的外轮廓点在图像坐标系o
‑
xy中的横坐标与纵坐标，是小天体真实轮廓的轮廓点数量。
15.此外，本发明还提供一种未知形状小天体的撞击中心计算系统，包括计算机设备，该计算机被编程或配置以执行所述未知形状小天体的撞击中心计算方法的步骤，或该计算机设备的存储器上存储有被编程或配置以执行所述未知形状小天体的撞击中心计算方法的计算机程序。
16.进一步的，本发明还提供一种卫星，该卫星带有所述未知形状小天体的撞击中心计算系统。
17.此外，本发明还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行所述未知形状小天体的撞击中心计算方法的计算机程序。
18.与现有技术相比，本发明技术方案带有的有益技术效果主要有：1.相比于轮廓扫描法、椭圆拟合法与互相关方法等，本发明方法的小天体撞击中心的计算不依赖于小天体形状的先验信息，增加了本发明方法在实际任务中的可行性；2.相比于矩方法，本发明方法能够减少复杂太阳光照条件对于撞击中心确定的影
响。本发明方法基于光照方向分割出图像中小天体的照明真轮廓，并判断其凹凸特性，选择骨架中心估计和矩方法中合适的方法，其受光照等因素影响更小，对于不规则的小天体更加稳定。
19.3.对于照明真轮廓呈现凹特性时，本发明方法能够降低小天体撞击中心的错失风险，提高小天体防御撞击偏转方案的可靠性。本发明方法在面对照明真轮廓呈现凹特性时将采用基于骨架的中心估计方法，该方法能够保证提取的撞击中心不会错失目标天体。
附图说明
20.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其它的附图。
21.图1为本发明所述未知形状小天体的撞击中心计算方法的流程图；图2为小孔成像模型与小天体可见部分轮廓示意图；图3是照明真轮廓、照明伪轮廓的特写示意以及将太阳光照方向旋转到y轴方向的示意图；图4为照明真轮廓凹凸特性示意图，其中（a）为凸特性示意图，（b）为凹特性示意图；图5为基于骨架/中轴线的中心估计方法、矩估计方法两种方法的示意图；图6为撞击中心安全边界距离示意图；图7为本发明所提方法与矩估计方法两种方法对小天体撞击中心估计的安全边界距离对比结果图。
22.附图标记：1——小天体；2——像平面；3——照明真轮廓；4——照明伪轮廓；5——列像素；6——提取的轮廓；7——骨架；8——骨架中心；9——矩中心；10——质心；11——小天体的真实轮廓；d1——骨架中心的安全边界距离；d2——矩中心的安全边界距离。
具体实施方式
23.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例不是本发明的全部实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。需要说明的是，在本发明中，除非另有明确的规定和限定，“平移”、“旋转”等用于描述空间相对位置关系的术语应做广义理解。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
24.本发明所述的未知形状小天体的撞击中心计算方法与传统的方法不同，其核心思路包括图1所示的5个步骤。下面结合图1至图7，详细说明本发明提供的未知形状小天体的撞击中心计算方法的步骤流程以及利用所述方法而取得的有益效果。
25.步骤s1，小天体的可见部分轮廓提取：基于捕获的目标未知形状小天体的图像，依次进行图像平滑、二值化、闭运算处理，并进行外轮廓提取，获得图像中小天体的可见部分
轮廓；假设执行小天体撞击任务的航天器携带有导航相机，导航相机能够捕获目标未知形状小天体1的图像，该航天器再利用相机所捕获的图像进行小天体撞击的导航、制导与控制。图2为小孔成像模型与小天体可见部分轮廓示意图，在图2中建立了3种坐标系，分别为：相机坐标系o
‑
xyz、图像坐标系o
‑
xy和像素坐标系o
’‑
uv。相机坐标系原点取相机光心o，z轴与相机光轴重合，且取成像方向为正向，即z轴方向为相机视轴方向，x轴与y轴、z轴构成右手直角坐标系。正成像平面，以下简称像平面2，与实际成像平面关于光学o中心对称，分别位于z=f和z=
‑
f的平面内，其中f为相机焦距。图像坐标系o
‑
xy建立在像平面2内，原点位于像平面的中心o，x轴指向像平面2的水平方向，y轴指向像平面2的竖直方向，图像坐标系中的x轴和y轴与相机坐标系的x轴和y轴平行。像素坐标系通过像元建立与图像坐标系的联系，像素坐标是给定像元尺寸下图像坐标系的数字离散值，且像素坐标系的原点在图像的左上角。
26.像平面上的组成轮廓的外轮廓点在图像坐标系中的坐标为，对应的像素坐标下的坐标为，其中是轮廓点数量。基于相机投影关系，图像坐标系中的坐标为的外轮廓点对应的是小天体1表面上的点，该点在相机坐标系中的坐标为，显然从光心出发，经过图像外轮廓点的光线与小天体表面相切，相切点正是。
27.基于捕获的图像依次进行如下操作：1、通过高斯滤波器对导航相机捕获的图像进行平滑；2、采用otsu阈值分割方法对平滑后的图像进行二值化，得到二值图像；3、对二值图像进行闭运算，填平图像中的小洞，弥合小裂缝，但目标小天体1在图像中总的位置和形状不变；4、采用行扫描法对图像进行轮廓提取，获取外轮廓，即图像中小天体的可见部分轮廓，提取的外轮廓参见图2。所述行扫描法可参见参考文献“使用行扫描法提取图像轮廓”（作者：郭鑫, 刘浩；期刊名称：机械制造与自动化；年份：2012；卷期：41(005)；页码：104
‑
109）。
28.步骤s2，将所述小天体的可见部分轮廓进行照明真轮廓与照明伪轮廓分割：计算太阳光照方向到图像坐标系的y轴方向的旋转矩阵，通过计算的所述旋转矩阵将包含小天体轮廓的图像进行旋转，然后对旋转后图像的每一列像素从下至上进行扫描，并按照从左至右的方向对所有列像素进行逐列扫描；对于每一列像素，取第一个遇到的轮廓点作为照明真轮廓点，所有的照明真轮廓点构成所述照明真轮廓；如果有第二个及之后遇到的轮廓点，则将其作为照明伪轮廓点，所有的照明伪轮廓点构成所述照明伪轮廓；参考图2所示的小孔成像模型，小天体1在阳光照射下，在相机的像平面2上成像。受太阳相角的影响，小天体1在像平面上可能只是部分可见。这意味小天体1的真实外轮廓在像平面2上只是部分可见，真实外轮廓的可见部分是受到光照的部分，称为照明真轮廓，只是像平面2上成像轮廓的一部分。像平面上成像轮廓的另一部分则称为照明伪轮廓，其实是小天体1表面中间的点在像平面2上的投影，它们是小天体1的可见与不可见部分的分界线，不可见部分由于处于阴影区，与昏暗的背景融为了一体。
29.通过太阳敏感器可以获得相机成像时太阳光照方向。这里给定太阳光照方向在图像中的投影单位矢量，其与图像坐标系x轴的夹角为：
(1)其中，是太阳光照方向在图像中的投影单位矢量在图像坐标系x轴方向的坐标分量，是太阳光照方向单位向量在图像坐标系y轴方向的坐标分量，是符号函数，arctan表示反正切函数。
30.参见图3，假设旋转矩阵将光照方向旋转到图像坐标系y轴方向，则可以计算旋转矩阵为(2)式中，sin表示正弦函数，cos表示余弦函数。
31.以太阳光照方向在图像中的投影单位方向向量旋转到y轴的旋转矩阵作用于图像中小天体可见部分轮廓，得到旋转后的轮廓。即(3)参见图3的右半部分，随后沿着旋转到y轴方向的太阳光照方向对每一列像素从下至上进行扫描，取第一个遇到的轮廓的点作为旋转后的照明真轮廓点，并按照图3中标注的列扫描方向对所有列像素进行逐列扫描，将所有的旋转后的照明真轮廓点构成旋转后的照明真轮廓，其中是照明真轮廓的轮廓点数量。如果有第二个及之后遇到的轮廓点，则作为旋转后的照明伪轮廓点，同样地，将对所有列进行扫描得到的这些旋转后的照明伪轮廓点一起构成旋转后的照明伪轮廓，其中是照明伪轮廓的轮廓点数量，且：。
32.之后可以得到旋转前的照明真轮廓：与照明伪轮廓：，计算如下：(4)(5)照明真轮廓如图3中的实线轮廓所示，照明伪轮廓如图3中的虚线轮廓所示。
33.步骤s3，识别所述照明真轮廓的凹凸特性：对所述照明真轮廓进行旋转，然后对旋
转后的照明真轮廓进行二次曲线拟合，根据所述二次曲线拟合的结果来判断所述照明真轮廓的凹凸特性；如图4所示，本发明将轮廓分为凸、凹两种类型，其中，图4（a）为凸特性示意图，图4（b）为凹特性示意图。凹凸分类是通过旋转后轮廓的拟合二次曲线的凹凸性来进行确定的。可以将轮廓点进行二次曲线拟合，拟合的二次曲线方程可以表示如下(6)其中，a0,a1,a2分别是二次曲线方程常数项、一次项、二次项的系数。二次曲线的凹凸特性由二次项系数符号决定，如果大于0，则二次曲线是凸的，如果如果小于0，则二次曲线是凹的。
34.假设基于照明真轮廓的轮廓点拟合得到的二次曲线方程为假设基于照明真轮廓的轮廓点拟合得到的二次曲线方程为，则定义照明真轮廓拟合误差平方和为(7)其中，分别是照明真轮廓的轮廓点拟合的二次曲线方程常数项、一次项、二次项的系数，是照明真轮廓的轮廓点的坐标，其中是照明真轮廓的轮廓点数量。
35.假设基于照明伪轮廓的轮廓点拟合得到的二次曲线为：则定义照明伪轮廓拟合误差平方和为：(8)其中，分别是照明伪轮廓的轮廓点拟合的二次曲线方程常数项、一次项、二次项的系数，照明伪轮廓的轮廓点的坐标，其中是照明伪轮廓的轮廓点数量。
36.则照明真轮廓和照明伪轮廓对应的最优拟合二次曲线方程的多项式系数的计算公式分别如下(9)(10)其中，是照明真轮廓对应的最优拟合二次曲线方程的多项式系数（分别为常数项、一次项、二次项系数），是照明伪轮廓对应的最优拟合二次曲线方程的多项式系数（分别为常数项、一次项、二次项系数），arg min表示最优化问题取极小值点的数学符号。上式可以通过最小二乘解析法或者牛顿下降法等数值方法进行求解。
37.通过判断的符号即可确定照明真轮廓的凹凸性，如果，则照明真轮廓是
凸的，否则照明真轮廓是凹的；同理判断的符号即可确定照明伪轮廓的凹凸性，如果，则照明伪轮廓是凸的，否则照明伪轮廓是凹的。但照明伪轮廓的凹凸特性并不采用。
38.步骤s4，选择估计方法获取小天体的撞击中心：如果判断到所述照明真轮廓是凹的，则采用基于骨架/中轴线的中心估计方法获取小天体的撞击中心；如果判断到所述照明真轮廓是凸的，则采用矩估计方法获取小天体的撞击中心；图5为基于骨架/中轴线的中心估计方法、矩估计方法两种方法的示意图，在图5中，6表示提取的轮廓；7为骨架；8为骨架中心；9为矩中心；10为质心；11表示小天体的真实轮廓。由图5可以看出：提取的轮廓6呈具有一定弧度的长条状，无法如实体现小天体的真实轮廓。在这种情况下，基于提取的轮廓6得到的矩中心9非常靠近小天体的边缘位置，将矩中心作为撞击中心具有非常大的风险，这本质上由于提取的轮廓6的特殊形状决定的；而基于提取的轮廓6得到的骨架中心8则相对地远离了小天体的边缘位置，作为撞击中心更加安全。
39.矩估计方法常用于估计质心，但是如图5所看到的，在面对凹形状时，错失目标的风险很大。相反，基于骨架/中轴线的中心估计方法的优势是能够保证中心位于目标上，但是其无法保证质心估计的准确性；将这两种方法组成撞击中心估计的方法池，通过判断形状凹凸特性，对撞击中心估计方法进行智能选择，能够弥补矩估计方法的不足。
40.本发明通过对照明真轮廓的凹凸性的判断来确定采用基于骨架/中轴线的中心估计方法还是采用矩估计方法。如步骤s3中所述，本发明对于首先对照明真轮廓的凹凸性进行判断，如果照明真轮廓是凹的，则采用基于骨架/中轴线的中心估计方法，基于骨架/中轴线的中心估计方法详见参考文献——徐超, 肖潇, 骆燕,等. 基于距离变换的新型骨架提取方法[j]. 仪器仪表学报, 2012, 33(012):2851
‑
2856.。如果照明真轮廓是凸的，就采用矩估计方法，矩估计方法详见参考文献——chang s ,grover c p . centroid detection based on optical correlation[j]. optical engineering, 2002, 41(10):2479
‑
2486。切换逻辑如图1所示。
[0041]
步骤s5，根据步骤s4获取的撞击中心计算所述撞击中心的安全边界距离。
[0042]
如图6所示，首先给出安全边界距离的定义。
[0043]
给定小天体真实轮廓，与小天体提取轮廓，其中，是小天体真实轮廓的轮廓点数量，是小天体提取轮廓的轮廓点数量。
[0044]
假设真实轮廓中心为，估计轮廓中心为，则撞击中心的安全边界距离是撞击中心距离真实小天体轮廓的符号距离，根据下式进行计算：(11)式中，表示估计轮廓中心；与分别是在图像坐标系下的横纵坐标，
表示小天体的真实轮廓，且，其中，是小天体真实轮廓的轮廓点数量。
[0045]
图6中，d1表示骨架中心的安全边界距离；d2表示矩中心的安全边界距离。其中，。
[0046]
小天体防御任务中，当面对尺寸较小、形状不规则且不确定的目标小天体，最重要的考虑就是任务的安全性。从安全保守性出发，在小天体形状非常不规则，且形状未知的情况下，除了继续评估小天体估计中心的偏差距离以外，还需要评估小天体错失撞击目标的可能性。本发明定义了安全边界距离，即小天体撞击中心到小天体真实边界的距离，安全边界距离越小，则小天体错失撞击目标的可能性就越大，这种情况当然是并不想看到的。因此，本发明以安全边界距离来评估本发明所提方法相比其它方法的性能优越性。
[0047]
为了验证所提方法的有效性，通过仿真对本发明所述方法与传统的矩估计方法的性能进行了测试，对比结果如图7所示，由图7可以看出，对于小天体的中心提取，所提方法的平均安全边界距离和最小安全边界距离在任何太阳相角的时候都大于矩估计方法的平均中心偏差距离和最大中心偏差距离，这表明所提方法的中心提取的安全性在非常不规则小天体上要远高于矩估计方法。可见本发明所提方法能够从安全性出发提取目标小天体的图像中心，进而为小天体撞击提供瞄准中心，降低撞击中心错失风险。
[0048]
综合以上分析，本发明所述未知形状小天体的撞击中心计算方法的流程图如图7所示。本发明以最小化中心错失风险为优化目标，设计了以小天体可见部分推测小天体撞击中心位置的算法。第一步，对小天体图像进行轮廓提取，并基于太阳光照方向的先验信息进行照明真轮廓与照明伪轮廓的分割；第二步，对小天体照明真轮廓与照明伪轮廓的凹凸特性进行识别；第三步，基于小天体照明真轮廓与照明伪轮廓的凹凸特性，智能切换为矩估计、骨架/中轴线中心估计两种撞击中心估计模式之一，并对中心错失风险进行估计。本发明采用不同规则度小天体模型，在多种航天器与小天体相对位姿条件下进行了小天体撞击中心提取仿真实验，并研究了光照条件对小天体撞击中心提取方法的影响。仿真结果表明即使在非常恶劣的光照与成像姿态条件下，所提方法也能够提供可靠且安全的小行撞击中心估计。由于小天体防御撞击偏转任务的面临着小天体形状模型不确定、小天体存在未知自转、光照条件复杂等问题，所提的高可靠高鲁棒小天体撞击中心估计算法将发挥巨大作用。
[0049]
以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效中心计算方法，或直接/间接运用在其它相关的技术领域均包括在本发明的专利保护范围内。