一种用于深空自主导航的不规则小天体图像质心确定方法

摘要

一种用于深空自主导航的不规则小天体图像质心确定方法，根据地面数据所生成的小天体三维形状模型确定其质心坐标系，根据不同方位的成像，建立小天体不同方位下的Hu氏7类归一化不变中心矩，在该方位下给出对应当前形心的质心修正系数，形成模型库；根据当前时刻拍摄图像，对小天体的轮廓进行提取；对轮廓二值图像计算一阶矩，即形心作为中心，计算Hu氏的归一化一类和二类二阶不变中心矩；计算二阶中心矩，进行后续匹配操作，如果匹配不成功，则继续计算三类至七类三阶中心矩，然后进行补充匹配；将计算出的中心矩与模型库中存储的中心矩进行匹配，当相似度满足设定阈值条件时，则认为当前位置匹配完成，将当前的形心利用质心系数进行修正。本发明保证在不同方位下质心输出的一致性。

说明书

技术领域

本发明涉及一种用于深空自主导航的不规则小天体图像质心确定方法，属 于航天技术中的图像处理领域。

背景技术

小天体是太阳系形成后的剩余物质，是太阳系的重要组成部分。它们保持 了太阳系起源时的原始状态，蕴含太阳系形成和演化的重要信息。研究小天体 对揭示太阳系的起源及演化、探索生命起源都具有非常重要的科学意义，对探 索和开发太阳系和应对可能的撞击地球事件具有重要的现实意义。

但是小天体不一定类似于球体，很可能是不规则的几何体。在飞越或交会 的过程中，飞行器利用光学敏感器对小天体进行成像，需要提供小天体的一个 固定中心点作为观测量，进行指向信息的实时更新。但是因为小天体的自旋和 相对运动，不同时刻的成像对应着小天体不同的方位地形。不仅外形轮廓不一 样，地面的纹理特征也不一致，特征点、灰度质心和形状中心都会发生变化。 所以怎样在图像中提取不变的中心点成为小天体探测的一大难题。

目前国内研究利用小天体指向信息进行导航的方法很多，但是大多假设小 天体的指向信息已获得，具体的获得方法不详。国外文献调研中，小天体的中 心就定义为图像中最大亮斑的中心，再进行补偿，但未详细介绍补偿系数的选 取和计算方法。

发明内容

本发明要解决的技术问题是：针对小天体外形不规则的特点，提取成像时 刻当前方位下小天体的轮廓信息，利用其形状中心求取多阶中心矩，采用不变 中心矩的平移、比例和旋转不变性，与模型库中存储的小天体不同方位的矩特 征进行匹配，得到小天体的当前方位，从而根据系数修正其质心坐标，保证在 不同方位下质心的一致性。

本发明的技术解决方案包括以下步骤：一种用于深空自主导航的不规则小 天体图像质心确定方法，实现步骤如下：

第一步，根据地面数据所生成的小天体三维形状模型确定其质心坐标系， 根据不同方位的成像，建立小天体不同方位下的Hu氏7类归一化不变中心矩 和比例因子。在该方位下给出对应当前形心的质心修正系数；

第二步，根据当前时刻拍摄图像，对小天体的轮廓进行提取。将图像变换 为仅含轮廓信息的二值图像。对轮廓图像计算一阶矩，即形心作为中心，计算 Hu氏的归一化一类和二类二阶不变中心矩。由于计算的需要，初始仅计算二阶 中心矩，然后进行后续匹配操作，如果匹配阈值不满足要求，则认为匹配失败， 再继续计算三类至七类三阶中心矩，然后进行补充匹配；

第三步，将计算出的中心矩与模型库中存储的中心矩进行匹配，当相似度 接近设定阈值时，则认为当前位置匹配完成，将当前的形心利用质心系数进行 修正，给出质心的指向信息。

所述第一步模型库的计算过程为：

(11)将小天体模型在等效轨道高度上进行投影成像，生成样本数据；

小天体三维形状模型已知，在标称的环绕轨道高度上对小天体进行投影成 像，直接形成轮廓图像作为样本图像，同时将质心(x_o1，y_o1)投影到像平面上。根 据指向精度的要求，确定样本图像的数量。例如精度要求为1度，则样本图像 的间隔角度要小于0.5度。

(12)计算Hu氏7类归一化不变中心矩和比例因子；

假设输入的轮廓二值图像为f(x，y)，x＝1，2，...，M，y＝1，2，...，N，图像的大小为 M×N，首先计算其一阶矩，即形心(x_c1，y_c1)：

$x_{c1}=\bar{x}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\mathrm{mf}(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)},$$y_{c1}=\bar{y}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\mathrm{nf}(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)}.$

然后计算二维中心矩：

${\eta }_{20}=\frac{u_{20}}{u_{00}^1}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^{2}f(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)},$${\eta }_{02}=\frac{u_{02}}{u_{00}^1}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{y})}^{2}f(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)},$

${\eta }_{11}=\frac{u_{11}}{u_{00}^1}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}(m-\bar{x})(n-\bar{y})f(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)}.$

其中$u_{\mathrm{pq}}=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^p{(n-\bar{y})}^{q}f(m,n).$

接着计算独立的一类二阶不变矩σ₁和二类二阶不变矩σ₂：

$\left(\begin{array}{c}{\sigma }_1={\eta }_{20}+{\eta }_{02}\\ {\sigma }_2={({\eta }_{20}-{\eta }_{02})}^2+4{\eta }_{11}^2\end{array}\right),$然后计算比例不变因子：

最后计算三阶矩：

${\eta }_{30}=\frac{u_{30}}{u_{00}^{3/2}}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^{3}f(m,n)}{{\left(\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)\right)}^{3/2}},$${\eta }_{03}=\frac{u_{03}}{u_{00}^{3/2}}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{y})}^{3}f(m,n)}{{\left(\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)\right)}^{3/2}},$

${\eta }_{12}=\frac{u_{12}}{u_{00}^{3/2}}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}(m-\bar{x}){(n-\bar{y})}^{2}f(m,n)}{{\left(\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)\right)}^{3/2}},$${\eta }_{21}=\frac{u_{21}}{u_{00}^{3/2}}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^2(n-\bar{y})f(m,n)}{{\left(\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)\right)}^{3/2}}$

由上述量计算独立的归一化三类、四类、五类、六类和七类三阶中心不变 矩存入数据库：

σ₃＝(η₃₀-3η₁₂)²+(η₀₃-3η₂₁)²

σ₄＝(η₃₀+η₁₂)²+(η₀₃+η₂₁)²

σ₅＝(η₃₀-3η₁₂)(η₃₀+η₁₂)[(η₃₀+η₁₂)²-3(η₀₃+η₂₁)²]

+(η₀₃-3η₂₁)(η₀₃+η₂₁)[(η₀₃+η₂₁)²-3(η₃₀+η₁₂)²]

σ₆＝(η₂₀-η₀₂)[(η₃₀+η₁₂)²-(η₀₃+η₂₁)²]+4η₁₁(η₃₀+η₁₂)(η₀₃+η₂₁)

σ₇＝(η₃₀-3η₂₁)(η₀₃+η₂₁)[(η₀₃+η₂₁)²-3(η₃₀+η₁₂)²]

-(η₀₃-3η₂₁)(η₃₀+η₁₂)[(η₃₀+η₁₂)²-3(η₀₃+η₂₁)²]

将计算出来的κ_i、σ₁到σ₇进行对应位置保存。

(13)给出当前方位形心到质心的修正系数；

计算图像中形心和质心的修正常量，$\left(\begin{array}{c}\mathrm{\Delta x}=x_{c1}-x_{o1}\\ \mathrm{\Delta y}=y_{c1}-y_{o1}\end{array}\right),$Δx和Δy为两个方向的修 正常量，对应每个样本进行保存。

所述第二步具体实现为：

(21)采用典型的边缘提取算法对拍摄的图像进行轮廓边缘提取；

(22)计算轮廓图像的一类和二类二阶不变矩；

假设输入的轮廓二值图像为f(x，y)，x＝1，2，…，M，y＝1，2，…，N，图像的大小为 M×N，首先计算其一阶矩，即形心：

$\bar{x}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\mathrm{mf}(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)},$$\bar{y}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\mathrm{nf}(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)}.$

然后计算二维中心矩：。

${\eta }_{20}=\frac{u_{20}}{u_{00}^1}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^{2}f(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)},$${\eta }_{02}=\frac{u_{02}}{u_{00}^1}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{y})}^{2}f(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)},$

${\eta }_{11}=\frac{u_{11}}{u_{00}^1}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}(m-\bar{x})(n-\bar{y})f(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)}.$

其中$u_{\mathrm{pq}}=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^p{(n-\bar{y})}^{q}f(m,n).$

最后计算独立的一类二阶不变矩和二类二阶不变矩

进行第三步的匹配过程，如果匹配不上，则进行下一步。

(23)计算其它高阶不变矩；

如果仅凭和无法确认匹配，则针对小天体轮廓图像，继续计算其三阶 矩：

${\eta }_{30}=\frac{u_{30}}{u_{00}^{3/2}}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^{3}f(m,n)}{{\left(\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)\right)}^{3/2}},$${\eta }_{03}=\frac{u_{03}}{u_{00}^{3/2}}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{y})}^{3}f(m,n)}{{\left(\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)\right)}^{3/2}},$

${\eta }_{12}=\frac{u_{12}}{u_{00}^{3/2}}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}(m-\bar{x}){(n-\bar{y})}^{2}f(m,n)}{{\left(\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)\right)}^{3/2}},$${\eta }_{21}=\frac{u_{21}}{u_{00}^{3/2}}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^2(n-\bar{y})f(m,n)}{{\left(\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)\right)}^{3/2}}$

然后计算独立的归一化三类、四类、五类、六类和七类三阶中心不变矩：

所述的第三步具体实现：

(31)比较一类和二类二阶不变矩的步骤；

计算比例不变因子：

将在模型库寻求相近的比例不变因子κ，|κ-κ_i|＜ε_κ，i从1到num，num为 模型库中存储的位置个数，ε_κ为设定的相似阈值。对得到的满足要求的n_κ个相 似位置再进行进一步精确匹配。计算模型库与真实图像的缩放因子σ_1i为 模型库n_κ个相似位置中第i个位置的一类二阶不变矩，i从1到n_κ。计算 其中σ_2i表示模型库n_κ个相似位置中第i个位置的二类 二阶不变矩。如果最小值和次小值之间差值的比值大于阈值ε，即表示次小值)，则认为其中的最小值对应的形状即为匹配形状，匹配正确。

(32)比较其它不变矩的步骤；

如果不能通过和完成匹配，则利用计算的和采用 下列公式计算或或重 新进行匹配确认，取其最小值，并判断最小和次小值的倍数是否满足ε。如果仍 无法确认，则输出和并给出此为形心的标识。

(33)对形心进行修正的过程。

根据匹配的结果，得到模型库中的修正系数，进行如下计算：

$\left(\begin{array}{c}{x}_o=x_c+\mathrm{\lambda \Delta x}\\ y_o=y_c+\mathrm{\lambda \Delta y}\end{array}\right),$x_c和y_c为小天体的形心，Δx和Δy为两个方向的修正常量， λ＝1/λ_i为缩放因子。

则质心指向计算如下：(x_odx，y_ody，f)，dx、dy为像元横、纵尺寸，f为相机 焦距。进行归一化后如下：

$(\frac{x_o\mathrm{dx}}{\sqrt{{\left(x_o\mathrm{dx}\right)}^2+{\left(y_o\mathrm{dy}\right)}^2+f^2}},\frac{y_o\mathrm{dy}}{\sqrt{{\left(x_o\mathrm{dx}\right)}^2+{\left(y_o\mathrm{dy}\right)}^2+f^2}},\frac{f}{\sqrt{{\left(x_o\mathrm{dx}\right)}^2+{\left(y_o\mathrm{dy}\right)}^2+f^2}}).$

本发明的优点是：

(1)本发明利用Hu氏归一化中心矩来表示小天体的轮廓形状，与模型库 的数据进行匹配，解决了小天体外形不规则，形心与质心不重合，提取不准确 的问题。

(2)采用7类不变矩作为模型库的存贮内容，极大优化了模型库的数据量， 提高了匹配的效率。

附图说明

图1为本发明质心提取方法的流程图；

图2为本发明中数学仿真生成的小天体图像；

图3为本发明中小天体的轮廓提取效果图。

具体实施方式

如图1所示，本发明具体实现如下：

第一步，小天体的轮廓提取。

采用Sobel算子模板：

$h1=\left(\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right)$$h2=\left(\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right)$

原图像矩阵F使用Sobel边缘提取有：

$E_1=F\otimes h1,$$E_2=F\otimes h2$

$E=\sqrt{E_1·E_1+E_2·E_2}$

其中为卷积运算符，E阵取阈值就可以得到二值化的边缘图，如图2为 原始图像，图3为图2对应的轮廓图像，白色的表示提取出来的像素轮廓。

第二步，小天体的归一化中心矩求取。

假设输入的轮廓二值图像为f(x，y)，x＝1，2，…，M，y＝1，2，…，N，图像的大小为 M×N，首先计算其一阶矩，即形心：

$\bar{x}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\mathrm{mf}(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)},$$\bar{y}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\mathrm{nf}(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)}.$

然后计算二维中心矩：。

${\eta }_{20}=\frac{u_{20}}{u_{00}^1}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^{2}f(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)},$${\eta }_{02}=\frac{u_{02}}{u_{00}^1}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{y})}^{2}f(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)},$

${\eta }_{11}=\frac{u_{11}}{u_{00}^1}=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}(m-\bar{x})(n-\bar{y})f(m,n)}{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}f(m,n)}.$

其中$u_{\mathrm{pq}}=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{(m-\bar{x})}^p{(n-\bar{y})}^{q}f(m,n).$

最后计算独立的一类二阶不变矩和二类二阶不变矩

图3中的轮廓计算得到的

第三步，与模型库归一化中心矩的匹配确认。

计算比例不变因子：

图3中的轮廓计算得到κ＝0.2661。

将在模型库寻求相近的比例不变因子κ，|κ-κ_i|＜ε_κ，i从1到num，num为 模型库中存储的位置个数，ε_κ为设定的相似阈值。对得到的满足要求的n_κ个相 似位置再进行进一步精确匹配。计算模型库与真实图像的缩放因子σ_1i为 模型库n_κ个相似位置中第i个位置的一类二阶不变矩，i从1到n_κ。计算 其中σ_2i表示模型库n_κ个相似位置中第i个位置的二类 二阶不变矩。如果最小值和次小值之间差值的比值大于阈值ε，即表示次小值)，则认为其中的最小值对应的形状即为匹配形状，得到其位置存贮 的质心修正系数，取其缩放因子λ，对中心进行修正。

对图3的轮廓图像而言，根据κ的取值，与模型库中的κ＝0.27和κ＝0.30的 图像最为近似，其中，κ＝0.27的图像σ₁₁＝2653.9329，σ₂₁＝519088.4199；κ＝0.30的 图像σ₁₂＝8415.2522，σ₂₂＝6697671.0698。计算λ₁＝0.1073，λ₂＝0.3401，且则认为κ＝0.27的位置与图3 匹配。

第四步，利用质心修正系数，修正中心，给出指向矢量。

根据匹配的结果，得到模型库中的修正系数，进行如下计算：

$\left(\begin{array}{c}{x}_o=x_c+\mathrm{\lambda \Delta x}\\ y_o=y_c+\mathrm{\lambda \Delta y}\end{array}\right),$x_c和y_c为小天体的形心，Δx和Δy为两个方向的修正常量， λ＝1/λ_i为缩放因子。

对应图3的结果为$\left(\begin{array}{c}{x}_o=251.79+\mathrm{\Delta x}/0.1073=245.09\\ y_o=247.5+\mathrm{\Delta y}/0.1073=253.45\end{array}\right)$

则质心指向计算如下：(x_odx，y_ody，f)，dx、dy为像元横纵尺寸，f为相机焦 距。进行归一化后如下：

$(\frac{x_o\mathrm{dx}}{\sqrt{{\left(x_o\mathrm{dx}\right)}^2+{\left(y_o\mathrm{dy}\right)}^2+f^2}},\frac{y_o\mathrm{dy}}{\sqrt{{\left(x_o\mathrm{dx}\right)}^2+{\left(y_o\mathrm{dy}\right)}^2+f^2}},\frac{f}{\sqrt{{\left(x_o\mathrm{dx}\right)}^2+{\left(y_o\mathrm{dy}\right)}^2+f^2}}).$

得到图3小天体的指向为(-0.0007439，-0.0001739，0.99999971)。

经过上述算法，得到小天体在不同方位和距离下成像，较为一致的质心指 向。