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

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

Info

Publication number
CN113628231A
CN113628231A CN202111179121.3A CN202111179121A CN113628231A CN 113628231 A CN113628231 A CN 113628231A CN 202111179121 A CN202111179121 A CN 202111179121A CN 113628231 A CN113628231 A CN 113628231A
Authority
CN
China
Prior art keywords
contour
illumination
celestial body
small celestial
center
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202111179121.3A
Other languages
English (en)
Other versions
CN113628231B (zh
Inventor
袁浩
王杰
刘望
罗青
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202111179121.3A priority Critical patent/CN113628231B/zh
Publication of CN113628231A publication Critical patent/CN113628231A/zh
Application granted granted Critical
Publication of CN113628231B publication Critical patent/CN113628231B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/12—Edge-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/60—Rotation of whole images or parts thereof
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/20—Image enhancement or restoration using local operators
    • G06T5/30—Erosion or dilatation, e.g. thinning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G06T7/73—Determining position or orientation of objects or cameras using feature-based methods

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开一种未知形状小天体的撞击中心计算方法及系统，所述小天体撞击中心计算方法包括以下步骤：步骤S1，提取图像中小行星可见部分轮廓；步骤S2，进行照明真轮廓/照明伪轮廓分割；步骤S3，识别照明真轮廓的凹凸特性；步骤S4，选择合适的估计方法获取小天体的撞击中心；步骤S5，计算撞击中心的安全边界距离。依据本发明所述方法进行未知形状小天体撞击中心计算，能够让小天体撞击中心的计算不依赖于小天体形状的先验信息；能够降低小天体撞击中心的错失风险，提高小天体防御撞击偏转方案的可靠性。

Description

一种未知形状小天体的撞击中心计算方法及系统
技术领域
本发明属于导航信息处理领域，具体涉及一种未知形状小天体的撞击中心计算方法及系统。
背景技术
小天体撞击中心计算是小天体防御撞击偏转光学导航的重要内容。现有的技术包括矩估计方法、轮廓扫描法、椭圆拟合法与互相关方法等，它们在计算复杂度、导航精度、适用天体类型、所需先验信息等方面各有特点。
矩估计方法简称矩方法，是最基本、也是最有效的中心提取算法。矩方法是利用矩算法对图像中目标天体的亮度中心进行估计。目标天体所占的像素范围通常由背景分割阈值决定，背景分割的目标是目标天体亮度中心的估计不被背景亮度以及其它的亮点所干扰而引起偏差。通过矩估计找到亮心以后，需要修正太阳相角所引起的亮心与天体中心的偏差效应。
轮廓扫描法是光学导航中常用的方法。首先要在图像中提取轮廓点，然后基于相机与小天体的相对位姿的先验信息与小天体模型表面的三维点进行对应。基于二维点与三维点的对应关系，可以构建PnP问题使用最小二乘法对图像中目标的中心位置进行估计更新。该方法需要小天体精确的三维模型先验信息。
椭圆拟合法适用于三维模型能用三轴椭球体良好描述的目标。它的基本思想是：当一个椭球体投影到二维图像平面上时，形成一个椭圆；这个椭圆可以用来估计相机相对于小天体的位姿。显然该类方法需要目标的先验信息，且前提是小天体需要能用三轴椭球体良好描述，但是对于未知的小尺寸小天体很难成立。
互相关方法通过预测的像元亮度值模板与真实图像的相关系数来确定目标中心，并可以通过拟合二维抛物面来进一步细化相关性峰值以确定表面的峰值亚像素位置。这种方法的应用很广泛，既可以用于静止光斑的中心计算，也可用于目标拖尾图像中心提取。但是这类方法需要基于小天体模型和成像场景的先验知识渲染小天体的预测视图来生成用于执行互相关的模板。
上述方法中只有矩估计方法是不需要先验模型的，但是矩估计方法在目标三维模型非常不规则，且非凸时，估计的中心位于三维模型的边缘，中心错失的风险较大。而对于轮廓扫描法、椭圆拟合法与互相关方法等其它方法，要么假设目标天体形状模型是三轴椭球体，要么假设目标天体形状的凸的半规则椭球体，要么假设目标是不规则的，但是这些模型都是精确已知的，由解析模型或三角片模型或球谐模型来表示。然而，在小天体防御的背景下以及随着越来越多的小尺寸小天体被发现的背景下，目标小天体的形状是非常不规则的，且是未知的，这让这些方法统统失效。因此，需要提出一种新的小天体撞击中心计算方法，利用该方法能够从安全性出发有效地提取目标小天体的图像中心，进而为小天体撞击提供瞄准中心，降低撞击中心错失风险。
发明内容
本发明的主要目的是提供一种未知形状小天体的撞击中心计算方法，旨在解决现有同类方法依赖于小天体形状先验信息和撞击中心错失风险较高的技术问题。
为实现上述目的，本发明提出一种未知形状小天体的撞击中心计算方法，包括以下步骤：
步骤S1，小天体的可见部分轮廓提取：基于捕获的目标未知形状小天体的图像，获得图像中小天体的可见部分轮廓；
步骤S2，将所述小天体的可见部分轮廓进行照明真轮廓与照明伪轮廓分割：计算太阳光照方向到图像坐标系的y轴方向的旋转矩阵，通过计算的所述旋转矩阵将包含小天体轮廓的图像进行旋转，然后对旋转后图像的每一列像素从下至上进行扫描，并按照从左至右的方向对所有列像素进行逐列扫描；对于每一列像素，取第一个遇到的轮廓点作为照明真轮廓点，所有的照明真轮廓点构成所述照明真轮廓；如果有第二个及之后遇到的轮廓点，则将其作为照明伪轮廓点，所有的照明伪轮廓点构成所述照明伪轮廓；
步骤S3，识别所述照明真轮廓的凹凸特性：对所述照明真轮廓进行旋转，然后对旋转后的照明真轮廓进行二次曲线拟合，根据所述二次曲线拟合的结果来判断所述照明真轮廓的凹凸特性；
步骤S4，选择估计方法获取小天体的撞击中心：如果判断到所述照明真轮廓是凹的，则采用基于骨架/中轴线的中心估计方法获取小天体的撞击中心；如果判断到所述照明真轮廓是凸的，则采用矩估计方法获取小天体的撞击中心；
步骤S5，根据步骤S4获取的撞击中心计算所述撞击中心的安全边界距离。
进一步的，所述步骤S1中，对捕获的图像依次进行图像平滑、二值化、闭运算处理，并进行外轮廓提取，获得图像中小天体的可见部分轮廓。
可选的，所述图像平滑采用高斯滤波器，所述二值化处理采用Otsu阈值分割方法，所述外轮廓提取采用行扫描法。
进一步的，所述步骤S2中，所述图像坐标系为建立在相机像平面内的直角坐标系o-xy，原点位于像平面的中心o，x轴指向像平面的水平方向，y轴指向像平面的竖直方向；
所述太阳光照方向到图像坐标系y轴方向的旋转矩阵
Figure 796507DEST_PATH_IMAGE001
为：
Figure 20815DEST_PATH_IMAGE002
式中，sin表示正弦函数，cos表示余弦函数，
Figure 82312DEST_PATH_IMAGE003
表示太阳光照方向在图像中的投影单位矢量与图像坐标系x轴的夹角且表达式为：
Figure 678247DEST_PATH_IMAGE004
其中，
Figure 90774DEST_PATH_IMAGE005
是太阳光照方向在图像中的投影单位矢量在图像坐标系x轴方向的坐标分量，
Figure 118773DEST_PATH_IMAGE006
是太阳光照方向单位向量在图像坐标系y轴方向的坐标分量，
Figure 972459DEST_PATH_IMAGE007
是符号函数，arctan表示反正切函数。
进一步的，基于所述照明真轮廓的轮廓点拟合得到的二次曲线方程为：
Figure 162132DEST_PATH_IMAGE008
Figure 61955DEST_PATH_IMAGE009
，
则定义照明真轮廓拟合误差平方和
Figure 628066DEST_PATH_IMAGE010
为
Figure 992051DEST_PATH_IMAGE011
式中，
Figure 87046DEST_PATH_IMAGE012
分别是基于所述照明真轮廓的轮廓点拟合得到的二次曲线方程的常数项、一次项、二次项的系数，
Figure 208586DEST_PATH_IMAGE013
分别是照明真轮廓的轮廓点的x坐标与y坐标分量，下标
Figure 578387DEST_PATH_IMAGE014
，其中
Figure 203403DEST_PATH_IMAGE015
表示照明真轮廓的轮廓点数量；
根据所述照明真轮廓拟合误差平方和
Figure 843201DEST_PATH_IMAGE016
，构建照明真轮廓对应的最优拟合二次曲线方程的多项式系数的计算公式如下
Figure 186458DEST_PATH_IMAGE017
其中，
Figure 94371DEST_PATH_IMAGE018
分别为照明真轮廓对应的最优拟合二次曲线方程的常数项、一次项、二次项的系数，arg min表示最优化问题取极小值点的数学符号；
通过判断所述多项式系数中
Figure 964107DEST_PATH_IMAGE019
的符号确定照明真轮廓的凹凸性：如果
Figure 932063DEST_PATH_IMAGE020
，则照明真轮廓是凸的，否则照明真轮廓是凹的。
进一步的，所述步骤S5中，所述安全边界距离为小天体的撞击中心距离真实小天体轮廓的符号距离，安全边界距离
Figure 762616DEST_PATH_IMAGE021
根据下式进行计算：
Figure 208640DEST_PATH_IMAGE022
，
式中，
Figure 808249DEST_PATH_IMAGE023
表示估计轮廓中心；
Figure 556893DEST_PATH_IMAGE024
与
Figure 874742DEST_PATH_IMAGE025
分别是
Figure 858879DEST_PATH_IMAGE026
在图像坐标系o-xy下的横坐标与纵坐标；
Figure 578573DEST_PATH_IMAGE027
表示小天体的真实轮廓且
Figure 747386DEST_PATH_IMAGE028
，其中，
Figure 552531DEST_PATH_IMAGE029
分别为小天体的外轮廓点在图像坐标系o-xy中的横坐标与纵坐标，
Figure 340358DEST_PATH_IMAGE030
是小天体真实轮廓的轮廓点数量。
此外，本发明还提供一种未知形状小天体的撞击中心计算系统，包括计算机设备，该计算机被编程或配置以执行所述未知形状小天体的撞击中心计算方法的步骤，或该计算机设备的存储器上存储有被编程或配置以执行所述未知形状小天体的撞击中心计算方法的计算机程序。
进一步的，本发明还提供一种卫星，该卫星带有所述未知形状小天体的撞击中心计算系统。
此外，本发明还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行所述未知形状小天体的撞击中心计算方法的计算机程序。
与现有技术相比，本发明技术方案带有的有益技术效果主要有：
1.相比于轮廓扫描法、椭圆拟合法与互相关方法等，本发明方法的小天体撞击中心的计算不依赖于小天体形状的先验信息，增加了本发明方法在实际任务中的可行性；
2.相比于矩方法，本发明方法能够减少复杂太阳光照条件对于撞击中心确定的影响。本发明方法基于光照方向分割出图像中小天体的照明真轮廓，并判断其凹凸特性，选择骨架中心估计和矩方法中合适的方法，其受光照等因素影响更小，对于不规则的小天体更加稳定。
3.对于照明真轮廓呈现凹特性时，本发明方法能够降低小天体撞击中心的错失风险，提高小天体防御撞击偏转方案的可靠性。本发明方法在面对照明真轮廓呈现凹特性时将采用基于骨架的中心估计方法，该方法能够保证提取的撞击中心不会错失目标天体。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其它的附图。
图1为本发明所述未知形状小天体的撞击中心计算方法的流程图；
图2为小孔成像模型与小天体可见部分轮廓示意图；
图3是照明真轮廓、照明伪轮廓的特写示意以及将太阳光照方向旋转到y轴方向的示意图；
图4为照明真轮廓凹凸特性示意图，其中（a）为凸特性示意图，（b）为凹特性示意图；
图5为基于骨架/中轴线的中心估计方法、矩估计方法两种方法的示意图；
图6为撞击中心安全边界距离示意图；
图7为本发明所提方法与矩估计方法两种方法对小天体撞击中心估计的安全边界距离对比结果图。
附图标记：
1——小天体；2——像平面；3——照明真轮廓；4——照明伪轮廓；5——列像素；6——提取的轮廓；7——骨架；8——骨架中心；9——矩中心；10——质心；11——小天体的真实轮廓；D1——骨架中心的安全边界距离；D2——矩中心的安全边界距离。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例不是本发明的全部实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。需要说明的是，在本发明中，除非另有明确的规定和限定，“平移”、“旋转”等用于描述空间相对位置关系的术语应做广义理解。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
本发明所述的未知形状小天体的撞击中心计算方法与传统的方法不同，其核心思路包括图1所示的5个步骤。下面结合图1至图7，详细说明本发明提供的未知形状小天体的撞击中心计算方法的步骤流程以及利用所述方法而取得的有益效果。
步骤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轴平行。像素坐标系通过像元建立与图像坐标系的联系，像素坐标是给定像元尺寸下图像坐标系的数字离散值，且像素坐标系的原点在图像的左上角。
像平面上的组成轮廓的外轮廓点在图像坐标系中的坐标为
Figure 180139DEST_PATH_IMAGE031
，对应的像素坐标下的坐标为
Figure 660798DEST_PATH_IMAGE032
，其中
Figure 953240DEST_PATH_IMAGE033
是轮廓点数量。基于相机投影关系，图像坐标系中的坐标为
Figure 121922DEST_PATH_IMAGE034
的外轮廓点对应的是小天体1表面上的点，该点在相机坐标系中的坐标为
Figure 550629DEST_PATH_IMAGE035
，显然从光心出发，经过图像外轮廓点
Figure 61245DEST_PATH_IMAGE036
的光线与小天体表面相切，相切点正是
Figure 106561DEST_PATH_IMAGE037
。
基于捕获的图像依次进行如下操作：1、通过高斯滤波器对导航相机捕获的图像进行平滑；2、采用Otsu阈值分割方法对平滑后的图像进行二值化，得到二值图像；3、对二值图像进行闭运算，填平图像中的小洞，弥合小裂缝，但目标小天体1在图像中总的位置和形状不变；4、采用行扫描法对图像进行轮廓提取，获取外轮廓
Figure 970612DEST_PATH_IMAGE038
，即图像中小天体的可见部分轮廓，提取的外轮廓参见图2。所述行扫描法可参见参考文献“使用行扫描法提取图像轮廓”（作者：郭鑫, 刘浩；期刊名称：机械制造与自动化；年份：2012；卷期：41(005)；页码：104-109）。
步骤S2，将所述小天体的可见部分轮廓进行照明真轮廓与照明伪轮廓分割：计算太阳光照方向到图像坐标系的y轴方向的旋转矩阵，通过计算的所述旋转矩阵将包含小天体轮廓的图像进行旋转，然后对旋转后图像的每一列像素从下至上进行扫描，并按照从左至右的方向对所有列像素进行逐列扫描；对于每一列像素，取第一个遇到的轮廓点作为照明真轮廓点，所有的照明真轮廓点构成所述照明真轮廓；如果有第二个及之后遇到的轮廓点，则将其作为照明伪轮廓点，所有的照明伪轮廓点构成所述照明伪轮廓；
参考图2所示的小孔成像模型，小天体1在阳光照射下，在相机的像平面2上成像。受太阳相角的影响，小天体1在像平面上可能只是部分可见。这意味小天体1的真实外轮廓在像平面2上只是部分可见，真实外轮廓的可见部分是受到光照的部分，称为照明真轮廓，只是像平面2上成像轮廓的一部分。像平面上成像轮廓的另一部分则称为照明伪轮廓，其实是小天体1表面中间的点在像平面2上的投影，它们是小天体1的可见与不可见部分的分界线，不可见部分由于处于阴影区，与昏暗的背景融为了一体。
通过太阳敏感器可以获得相机成像时太阳光照方向。这里给定太阳光照方向在图像中的投影单位矢量
Figure 253826DEST_PATH_IMAGE039
，其与图像坐标系x轴的夹角为：
Figure 341867DEST_PATH_IMAGE040
(1)
其中，
Figure 218688DEST_PATH_IMAGE041
是太阳光照方向在图像中的投影单位矢量在图像坐标系x轴方向的坐标分量，
Figure 152009DEST_PATH_IMAGE042
是太阳光照方向单位向量在图像坐标系y轴方向的坐标分量，
Figure 289729DEST_PATH_IMAGE043
是符号函数，arctan表示反正切函数。
参见图3，假设旋转矩阵
Figure 283093DEST_PATH_IMAGE001
将光照方向
Figure 37422DEST_PATH_IMAGE044
旋转到图像坐标系y轴方向，则可以计算旋转矩阵
Figure 774434DEST_PATH_IMAGE001
为
Figure 297819DEST_PATH_IMAGE045
(2)
式中，sin表示正弦函数，cos表示余弦函数。
以太阳光照方向在图像中的投影单位方向向量
Figure 727663DEST_PATH_IMAGE046
旋转到y轴的旋转矩阵
Figure 766027DEST_PATH_IMAGE001
作用于图像中小天体可见部分轮廓
Figure 306729DEST_PATH_IMAGE047
，得到旋转后的轮廓
Figure 950200DEST_PATH_IMAGE048
。即
Figure 285367DEST_PATH_IMAGE049
(3)
参见图3的右半部分，随后沿着旋转到y轴方向的太阳光照方向对每一列像素从下至上进行扫描，取第一个遇到的轮廓
Figure 14288DEST_PATH_IMAGE050
的点作为旋转后的照明真轮廓点，并按照图3中标注的列扫描方向对所有列像素进行逐列扫描，将所有的旋转后的照明真轮廓点构成旋转后的照明真轮廓
Figure 257882DEST_PATH_IMAGE051
，其中
Figure 755860DEST_PATH_IMAGE052
是照明真轮廓的轮廓点数量。如果有第二个及之后遇到的轮廓点，则作为旋转后的照明伪轮廓点，同样地，将对所有列进行扫描得到的这些旋转后的照明伪轮廓点一起构成旋转后的照明伪轮廓
Figure 261927DEST_PATH_IMAGE053
，其中
Figure 478145DEST_PATH_IMAGE054
是照明伪轮廓的轮廓点数量，且：
Figure 173700DEST_PATH_IMAGE055
。
之后可以得到旋转前的照明真轮廓：
Figure 260605DEST_PATH_IMAGE056
与照明伪轮廓：
Figure 937574DEST_PATH_IMAGE057
，计算如下：
Figure 109929DEST_PATH_IMAGE058
(4)
Figure 389601DEST_PATH_IMAGE059
(5)
照明真轮廓如图3中的实线轮廓所示，照明伪轮廓如图3中的虚线轮廓所示。
步骤S3，识别所述照明真轮廓的凹凸特性：对所述照明真轮廓进行旋转，然后对旋转后的照明真轮廓进行二次曲线拟合，根据所述二次曲线拟合的结果来判断所述照明真轮廓的凹凸特性；
如图4所示，本发明将轮廓分为凸、凹两种类型，其中，图4（a）为凸特性示意图，图4（b）为凹特性示意图。凹凸分类是通过旋转后轮廓的拟合二次曲线的凹凸性来进行确定的。可以将轮廓点进行二次曲线拟合，拟合的二次曲线方程可以表示如下
Figure 596591DEST_PATH_IMAGE060
(6)
其中，a 0,a 1,a 2分别是二次曲线方程常数项、一次项、二次项的系数。二次曲线的凹凸特性由二次项系数
Figure 178882DEST_PATH_IMAGE061
符号决定，如果
Figure 104113DEST_PATH_IMAGE061
大于0，则二次曲线是凸的，如果如果
Figure 62841DEST_PATH_IMAGE061
小于0，则二次曲线是凹的。
假设基于照明真轮廓的轮廓点拟合得到的二次曲线方程为
Figure 498240DEST_PATH_IMAGE062
Figure 517011DEST_PATH_IMAGE063
，则定义照明真轮廓拟合误差平方和
Figure 663959DEST_PATH_IMAGE064
为
Figure 160799DEST_PATH_IMAGE065
(7)
其中，
Figure 201436DEST_PATH_IMAGE066
分别是照明真轮廓的轮廓点拟合的二次曲线方程常数项、一次项、二次项的系数，
Figure 391109DEST_PATH_IMAGE067
是照明真轮廓的轮廓点的坐标，其中
Figure 25353DEST_PATH_IMAGE068
是照明真轮廓的轮廓点数量。
假设基于照明伪轮廓的轮廓点拟合得到的二次曲线为：
Figure 591464DEST_PATH_IMAGE069
则定义照明伪轮廓拟合误差平方和
Figure 627553DEST_PATH_IMAGE070
为：
Figure 253706DEST_PATH_IMAGE071
(8)
其中，
Figure 453874DEST_PATH_IMAGE072
分别是照明伪轮廓的轮廓点拟合的二次曲线方程常数项、一次项、二次项的系数，
Figure 558097DEST_PATH_IMAGE073
照明伪轮廓的轮廓点的坐标，其中
Figure 714271DEST_PATH_IMAGE074
是照明伪轮廓的轮廓点数量。
则照明真轮廓和照明伪轮廓对应的最优拟合二次曲线方程的多项式系数的计算公式分别如下
Figure 511326DEST_PATH_IMAGE075
(9)
Figure 120162DEST_PATH_IMAGE076
(10)
其中，
Figure 355971DEST_PATH_IMAGE077
是照明真轮廓对应的最优拟合二次曲线方程的多项式系数（分别为常数项、一次项、二次项系数），
Figure 101073DEST_PATH_IMAGE078
是照明伪轮廓对应的最优拟合二次曲线方程的多项式系数（分别为常数项、一次项、二次项系数），arg min表示最优化问题取极小值点的数学符号。上式可以通过最小二乘解析法或者牛顿下降法等数值方法进行求解。
通过判断
Figure 177352DEST_PATH_IMAGE079
的符号即可确定照明真轮廓的凹凸性，如果
Figure 7904DEST_PATH_IMAGE080
，则照明真轮廓是凸的，否则照明真轮廓是凹的；同理判断
Figure 453929DEST_PATH_IMAGE081
的符号即可确定照明伪轮廓的凹凸性，如果
Figure 319117DEST_PATH_IMAGE082
，则照明伪轮廓是凸的，否则照明伪轮廓是凹的。但照明伪轮廓的凹凸特性并不采用。
步骤S4，选择估计方法获取小天体的撞击中心：如果判断到所述照明真轮廓是凹的，则采用基于骨架/中轴线的中心估计方法获取小天体的撞击中心；如果判断到所述照明真轮廓是凸的，则采用矩估计方法获取小天体的撞击中心；
图5为基于骨架/中轴线的中心估计方法、矩估计方法两种方法的示意图，在图5中，6表示提取的轮廓；7为骨架；8为骨架中心；9为矩中心；10为质心；11表示小天体的真实轮廓。由图5可以看出：提取的轮廓6呈具有一定弧度的长条状，无法如实体现小天体的真实轮廓。在这种情况下，基于提取的轮廓6得到的矩中心9非常靠近小天体的边缘位置，将矩中心作为撞击中心具有非常大的风险，这本质上由于提取的轮廓6的特殊形状决定的；而基于提取的轮廓6得到的骨架中心8则相对地远离了小天体的边缘位置，作为撞击中心更加安全。
矩估计方法常用于估计质心，但是如图5所看到的，在面对凹形状时，错失目标的风险很大。相反，基于骨架/中轴线的中心估计方法的优势是能够保证中心位于目标上，但是其无法保证质心估计的准确性；将这两种方法组成撞击中心估计的方法池，通过判断形状凹凸特性，对撞击中心估计方法进行智能选择，能够弥补矩估计方法的不足。
本发明通过对照明真轮廓的凹凸性的判断来确定采用基于骨架/中轴线的中心估计方法还是采用矩估计方法。如步骤S3中所述，本发明对于首先对照明真轮廓的凹凸性进行判断，如果照明真轮廓是凹的，则采用基于骨架/中轴线的中心估计方法，基于骨架/中轴线的中心估计方法详见参考文献——徐超, 肖潇, 骆燕,等. 基于距离变换的新型骨架提取方法[J]. 仪器仪表学报, 2012, 33(012):2851-2856.。如果照明真轮廓是凸的，就采用矩估计方法，矩估计方法详见参考文献——Chang S ,Grover C P . Centroid detectionbased on optical correlation[J]. Optical Engineering, 2002, 41(10):2479-2486。切换逻辑如图1所示。
步骤S5，根据步骤S4获取的撞击中心计算所述撞击中心的安全边界距离。
如图6所示，首先给出安全边界距离的定义。
给定小天体真实轮廓
Figure 192395DEST_PATH_IMAGE083
，与小天体提取轮廓
Figure 775823DEST_PATH_IMAGE084
，其中，
Figure 25539DEST_PATH_IMAGE085
是小天体真实轮廓的轮廓点数量，
Figure 10812DEST_PATH_IMAGE086
是小天体提取轮廓的轮廓点数量。
假设真实轮廓中心为
Figure 382888DEST_PATH_IMAGE087
，估计轮廓中心为
Figure 453612DEST_PATH_IMAGE088
，则撞击中心
Figure 507019DEST_PATH_IMAGE089
的安全边界距离
Figure 346799DEST_PATH_IMAGE090
是撞击中心
Figure 640508DEST_PATH_IMAGE089
距离真实小天体轮廓
Figure 932949DEST_PATH_IMAGE091
的符号距离，根据下式进行计算：
Figure 790047DEST_PATH_IMAGE092
(11)
式中，
Figure 484333DEST_PATH_IMAGE093
表示估计轮廓中心；
Figure 198211DEST_PATH_IMAGE094
与
Figure 243528DEST_PATH_IMAGE095
分别是
Figure 638737DEST_PATH_IMAGE096
在图像坐标系下的横纵坐标，
Figure 187530DEST_PATH_IMAGE097
表示小天体的真实轮廓，且
Figure 275572DEST_PATH_IMAGE098
，其中，
Figure 542605DEST_PATH_IMAGE085
是小天体真实轮廓的轮廓点数量。
图6中，D1表示骨架中心
Figure 646565DEST_PATH_IMAGE099
的安全边界距离；D2表示矩中心
Figure 49864DEST_PATH_IMAGE100
的安全边界距离。其中，
Figure 43228DEST_PATH_IMAGE101
。
小天体防御任务中，当面对尺寸较小、形状不规则且不确定的目标小天体，最重要的考虑就是任务的安全性。从安全保守性出发，在小天体形状非常不规则，且形状未知的情况下，除了继续评估小天体估计中心的偏差距离以外，还需要评估小天体错失撞击目标的可能性。本发明定义了安全边界距离
Figure 531978DEST_PATH_IMAGE102
，即小天体撞击中心到小天体真实边界的距离，安全边界距离越小，则小天体错失撞击目标的可能性就越大，这种情况当然是并不想看到的。因此，本发明以安全边界距离来评估本发明所提方法相比其它方法的性能优越性。
为了验证所提方法的有效性，通过仿真对本发明所述方法与传统的矩估计方法的性能进行了测试，对比结果如图7所示，由图7可以看出，对于小天体的中心提取，所提方法的平均安全边界距离和最小安全边界距离在任何太阳相角的时候都大于矩估计方法的平均中心偏差距离和最大中心偏差距离，这表明所提方法的中心提取的安全性在非常不规则小天体上要远高于矩估计方法。可见本发明所提方法能够从安全性出发提取目标小天体的图像中心，进而为小天体撞击提供瞄准中心，降低撞击中心错失风险。
综合以上分析，本发明所述未知形状小天体的撞击中心计算方法的流程图如图7所示。本发明以最小化中心错失风险为优化目标，设计了以小天体可见部分推测小天体撞击中心位置的算法。第一步，对小天体图像进行轮廓提取，并基于太阳光照方向的先验信息进行照明真轮廓与照明伪轮廓的分割；第二步，对小天体照明真轮廓与照明伪轮廓的凹凸特性进行识别；第三步，基于小天体照明真轮廓与照明伪轮廓的凹凸特性，智能切换为矩估计、骨架/中轴线中心估计两种撞击中心估计模式之一，并对中心错失风险进行估计。本发明采用不同规则度小天体模型，在多种航天器与小天体相对位姿条件下进行了小天体撞击中心提取仿真实验，并研究了光照条件对小天体撞击中心提取方法的影响。仿真结果表明即使在非常恶劣的光照与成像姿态条件下，所提方法也能够提供可靠且安全的小行撞击中心估计。由于小天体防御撞击偏转任务的面临着小天体形状模型不确定、小天体存在未知自转、光照条件复杂等问题，所提的高可靠高鲁棒小天体撞击中心估计算法将发挥巨大作用。
以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效中心计算方法，或直接/间接运用在其它相关的技术领域均包括在本发明的专利保护范围内。

Claims (9)

1.一种未知形状小天体的撞击中心计算方法，其特征在于包括以下步骤：
步骤S1，小天体的可见部分轮廓提取：基于捕获的目标未知形状小天体的图像，提取图像中小天体的可见部分轮廓；
步骤S2，将所述小天体的可见部分轮廓进行照明真轮廓与照明伪轮廓分割：计算太阳光照方向到图像坐标系的y轴方向的旋转矩阵，通过计算的所述旋转矩阵将包含小天体轮廓的图像进行旋转，然后对旋转后图像的每一列像素从下至上进行扫描，并按照从左至右的方向对所有列像素进行逐列扫描；对于每一列像素，取第一个遇到的轮廓点作为照明真轮廓点，所有的照明真轮廓点构成所述照明真轮廓；如果有第二个及之后遇到的轮廓点则作为照明伪轮廓点，所有的照明伪轮廓点构成所述照明伪轮廓；
步骤S3，识别所述照明真轮廓的凹凸特性：对所述照明真轮廓进行旋转，然后对旋转后的照明真轮廓进行二次曲线拟合，根据所述二次曲线拟合的结果来判断所述照明真轮廓的凹凸特性；
步骤S4，选择估计方法获取小天体的撞击中心：如果判断到所述照明真轮廓是凹的，则采用基于骨架/中轴线的中心估计方法获取小天体的撞击中心；如果判断到所述照明真轮廓是凸的，则采用矩估计方法获取小天体的撞击中心；
步骤S5，根据步骤S4获取的撞击中心计算所述撞击中心的安全边界距离。
2.根据权利要求1所述的未知形状小天体的撞击中心计算方法，其特征在于，所述步骤S1中，对捕获的图像依次进行图像平滑、二值化、闭运算处理，并进行外轮廓提取，获得图像中小天体的可见部分轮廓。
3.根据权利要求2所述的未知形状小天体的撞击中心计算方法，其特征在于，所述图像平滑采用高斯滤波器，所述二值化处理采用Otsu阈值分割方法，所述外轮廓提取采用行扫描法。
4.根据权利要求1所述的未知形状小天体的撞击中心计算方法，其特征在于，所述步骤S2中，所述图像坐标系为建立在相机像平面内的直角坐标系o-xy，原点位于像平面的中心o，x轴指向像平面的水平方向，y轴指向像平面的竖直方向；
所述太阳光照方向到图像坐标系y轴方向的旋转矩阵
Figure 527477DEST_PATH_IMAGE001
为：
Figure 444618DEST_PATH_IMAGE002
式中，sin表示正弦函数，cos表示余弦函数，
Figure 224355DEST_PATH_IMAGE003
表示太阳光照方向在图像中的投影单位矢量与图像坐标系x轴的夹角且表达式为：
Figure 681881DEST_PATH_IMAGE004
其中，
Figure 496253DEST_PATH_IMAGE005
是太阳光照方向在图像中的投影单位矢量在图像坐标系x轴方向的坐标分量，
Figure 318716DEST_PATH_IMAGE006
是太阳光照方向单位向量在图像坐标系y轴方向的坐标分量，
Figure 585749DEST_PATH_IMAGE007
是符号函数，arctan表示反正切函数。
5.根据权利要求1所述的未知形状小天体的撞击中心计算方法，其特征在于，所述步骤S3的详细过程为：
基于所述照明真轮廓的轮廓点拟合得到的二次曲线方程为
Figure 784649DEST_PATH_IMAGE008
Figure 759253DEST_PATH_IMAGE009
，则定义照明真轮廓拟合误差平方和
Figure 752617DEST_PATH_IMAGE010
为
Figure 506946DEST_PATH_IMAGE011
式中，
Figure 243958DEST_PATH_IMAGE012
分别是基于所述照明真轮廓的轮廓点拟合得到的二次曲线方程的常数项、一次项、二次项的系数，
Figure 767343DEST_PATH_IMAGE013
分别是照明真轮廓的轮廓点的x坐标与y坐标分量，下标
Figure 993925DEST_PATH_IMAGE014
，其中
Figure 235551DEST_PATH_IMAGE015
表示照明真轮廓的轮廓点数量；
根据所述照明真轮廓拟合误差平方和
Figure 510674DEST_PATH_IMAGE016
，构建照明真轮廓对应的最优拟合二次曲线方程的多项式系数的计算公式如下
Figure 154145DEST_PATH_IMAGE017
其中，
Figure 489312DEST_PATH_IMAGE018
分别为照明真轮廓对应的最优拟合二次曲线方程的常数项、一次项、二次项的系数，arg min表示最优化问题取极小值点的数学符号；
通过判断所述多项式系数中
Figure 31283DEST_PATH_IMAGE019
的符号确定照明真轮廓的凹凸性：如果
Figure 110097DEST_PATH_IMAGE020
，则照明真轮廓是凸的，否则照明真轮廓是凹的。
6.根据权利要求1~5中任意一项所述的未知形状小天体的撞击中心计算方法，其特征在于，所述步骤S5中，所述安全边界距离为小天体的撞击中心距离真实小天体轮廓的符号距离，安全边界距离
Figure 608075DEST_PATH_IMAGE021
根据下式进行计算：
Figure 379722DEST_PATH_IMAGE022
式中，
Figure 330360DEST_PATH_IMAGE023
表示估计轮廓中心；
Figure 212865DEST_PATH_IMAGE024
与
Figure 362087DEST_PATH_IMAGE025
分别是
Figure 412957DEST_PATH_IMAGE026
在图像坐标系o-xy下的横坐标与纵坐标；
Figure 585313DEST_PATH_IMAGE027
表示小天体的真实轮廓且
Figure 740350DEST_PATH_IMAGE028
，其中，
Figure 212920DEST_PATH_IMAGE029
分别为小天体的外轮廓点在图像坐标系o-xy中的横坐标与纵坐标，
Figure 60790DEST_PATH_IMAGE030
是小天体真实轮廓的轮廓点数量。
7.一种未知形状小天体的撞击中心计算系统，包括计算机设备，其特征在于，该计算机被编程或配置以执行权利要求1~6中任意一项所述未知形状小天体的撞击中心计算方法的步骤，或该计算机设备的存储器上存储有被编程或配置以执行权利要求1~6中任意一项所述未知形状小天体的撞击中心计算方法的计算机程序。
8.一种卫星，其特征在于，该卫星带有权利要求7所述未知形状小天体的撞击中心计算系统。
9.一种计算机可读存储介质，其特征在于，该计算机可读存储介质上存储有被编程或配置以执行权利要求1~6中任意一项所述未知形状小天体的撞击中心计算方法的计算机程序。
CN202111179121.3A 2021-10-11 2021-10-11 一种未知形状小天体的撞击中心计算方法及系统 Active CN113628231B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111179121.3A CN113628231B (zh) 2021-10-11 2021-10-11 一种未知形状小天体的撞击中心计算方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111179121.3A CN113628231B (zh) 2021-10-11 2021-10-11 一种未知形状小天体的撞击中心计算方法及系统

Publications (2)

Publication Number Publication Date
CN113628231A true CN113628231A (zh) 2021-11-09
CN113628231B CN113628231B (zh) 2022-01-25

Family

ID=78390825

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111179121.3A Active CN113628231B (zh) 2021-10-11 2021-10-11 一种未知形状小天体的撞击中心计算方法及系统

Country Status (1)

Country Link
CN (1) CN113628231B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102054122A (zh) * 2010-10-27 2011-05-11 东南大学 基于实际测量的纹理的力触觉表达方法
US20170361875A1 (en) * 2016-06-16 2017-12-21 Honda Motor Co., Ltd. Vehicle body side structure
CN109669481A (zh) * 2019-01-24 2019-04-23 中国人民解放军国防科技大学 基于等碰撞概率面法的航天器安全接近控制方法
CN110480678A (zh) * 2019-07-19 2019-11-22 南京埃斯顿机器人工程有限公司 一种工业机器人碰撞检测方法
CN112037211A (zh) * 2020-09-04 2020-12-04 中国空气动力研究与发展中心超高速空气动力研究所 一种用于微小空间碎片撞击事件动态监测的损伤特征识别方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102054122A (zh) * 2010-10-27 2011-05-11 东南大学 基于实际测量的纹理的力触觉表达方法
US20170361875A1 (en) * 2016-06-16 2017-12-21 Honda Motor Co., Ltd. Vehicle body side structure
CN109669481A (zh) * 2019-01-24 2019-04-23 中国人民解放军国防科技大学 基于等碰撞概率面法的航天器安全接近控制方法
CN110480678A (zh) * 2019-07-19 2019-11-22 南京埃斯顿机器人工程有限公司 一种工业机器人碰撞检测方法
CN112037211A (zh) * 2020-09-04 2020-12-04 中国空气动力研究与发展中心超高速空气动力研究所 一种用于微小空间碎片撞击事件动态监测的损伤特征识别方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
王云财 等: "深空自主导航系统可视化仿真平台", 《计算机仿真》 *

Also Published As

Publication number Publication date
CN113628231B (zh) 2022-01-25

Similar Documents

Publication Publication Date Title
Shan et al. A stereovision-based crack width detection approach for concrete surface assessment
CN111899334B (zh) 一种基于点线特征的视觉同步定位与地图构建方法及装置
CN108122256B (zh) 一种逼近状态下旋转目标位姿测量的方法
CN107228860B (zh) 一种基于图像旋转周期特性的齿轮缺陷检测方法
CN108225319B (zh) 基于目标特征的单目视觉快速相对位姿估计系统及方法
CN111652896B (zh) 惯导辅助陨石坑由粗到精检测方法
JP6483168B2 (ja) ビジョンシステムで画像内のプローブを効率的に採点するためのシステム及び方法
US10043279B1 (en) Robust detection and classification of body parts in a depth map
Li et al. A system of the shadow detection and shadow removal for high resolution city aerial photo
CN112435252B (zh) 一种战斗部破片穿孔和凹坑检测方法
CN116433701B (zh) 一种工件孔洞轮廓的提取方法、装置、设备及存储介质
CN117576219A (zh) 大广角鱼眼镜头单次拍摄图像的相机标定设备及标定方法
Lv et al. A point tracking method of TDDM for vibration measurement and large-scale rotational motion tracking
CN113705564B (zh) 一种指针式仪表识别读数方法
CN105574860B (zh) 一种舵片偏转极性识别方法
CN113628231B (zh) 一种未知形状小天体的撞击中心计算方法及系统
Wang et al. Hand posture recognition from disparity cost map
CN117635874A (zh) 基于图像多特征提取与融合算法的破片场三维快速反演系统
CN117611525A (zh) 一种受电弓滑板磨耗视觉检测方法和系统
CN110738098A (zh) 一种目标的识别定位与锁定跟踪方法
CN114049380B (zh) 目标物体定位追踪方法、装置、计算机设备和存储介质
Oumer Monocular 3D pose tracking of a specular object
Alkazzaz et al. Precise pupil boundary detection using angular integral projection and active contour model
CN112446894B (zh) 基于方向空间的图像分割方法
CN118485670A (zh) 基于视觉的射击自动报靶方法、装置、设备及介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant