CN104748750B

CN104748750B - 一种模型约束下的在轨三维空间目标姿态估计方法及系统

Info

Publication number: CN104748750B
Application number: CN201310740553.6A
Authority: CN
Inventors: 张天序; 王亮亮; 周钢; 李明; 杨卫东; 刘宽; 郑亚云
Original assignee: Huazhong University of Science and Technology
Current assignee: Huazhong University of Science and Technology
Priority date: 2013-12-28
Filing date: 2013-12-28
Publication date: 2015-12-02
Anticipated expiration: 2033-12-28
Also published as: CN104748750A; US20170008650A1; WO2015096508A1

Abstract

本发明公开一种在轨三维空间目标姿态估计方法及系统，包括离线特征库构建和在线姿态估计步骤，所述离线特征库构建步骤为：根据空间目标三维模型获取目标多视点特性视图，从中提取几何特征构成几何特征库，包含目标主体高宽比、目标纵向对称度、目标横向对称度和目标主轴倾斜角；所述在线姿态估计步骤为：对待测在轨目标图像预处理和提取特征，将提取的特征在几何特征库中匹配，匹配结果对应的特性视图所表征的目标姿态即为姿态估计结果。本发明匹配用的几何特征具有尺度不变性，只要在三维建模阶段准确获取目标各部件间的尺寸比例和位置关系，就可保证后续较高的匹配精度。整个方法实现简单、鲁棒性好、姿态估计精度高，受成像条件影响小，适用性好。

Description

一种模型约束下的在轨三维空间目标姿态估计方法及系统

技术领域

本发明属于航天技术与模式识别交叉领域，具体涉及一种适用于卫星、空间飞行器等在轨三维空间目标姿态估计方法及系统。

背景技术

国内外发射的大量通信卫星、资源卫星等空间目标可用于网络通信、航空摄影、大地测量等应用场所。对这些空间目标进行地基光电观察，分析、调整其姿态是该类系统中不可或缺的部分。由于地基望远镜系统的空间分辨率限制以及大气环境对长距离光学成像的随机干扰，使得地基传感器获取的图像容易出现目标边界模糊不清的现象。当成像目标边界模糊不清时，传统的基于特征点匹配的姿态估计及三维重建算法准确性往往会随着目标的模糊程度增加而迅速下降。姿态估计意在从二维相机坐标系获取的目标投影图像中计算出目标在三维目标坐标系下的俯仰角α和偏航角度β，一对(α,β)的角度值对应一个姿态。姿态估计的准确性对分析空间目标的部件尺寸、部件相对位置关系以及空间目标的功能属性等具有重要意义。因此，有必要研究地基长距离光学成像条件下稳健的姿态估计算法。

国内外学者对这类成像下的空间目标姿态估计算法进行了详细的研究，并且取得了相关成果。如，赵汝进、张启衡、徐智勇的“一种基于特征点间线段倾角的姿态测量方法”，见《光子学报》，2010年2月，第39卷，第2期。研究了一种基于目标特征点间倾角角度信息的目标3维姿态迭代解算方法，该方法适用于远距离弱透视成像目标和相机内参量未知条件下目标姿态求解。但算法精度严重依赖于提取到的边缘、直线、角点精度，且在迭代初值偏离真实姿态误差较大时，算法需要较多的迭代次数，运算量大，且可能出现迭代不收敛情况。地基长距离光学成像目标边界容易出现模糊不清的现象，影响特征点定位精度，因此，算法精度差。王鲲鹏、张小虎、于起峰的“基于目标特征点比例的单站图像定姿方法”，见《应用光学》2009年11月，第30卷，第6期，提出了一种针对实况记录图像的单站定姿方法，利用目标特征点坐标的比例信息及目标成像模型及坐标系间位姿参数关系，采用非线性方程组的迭代求解获得目标姿态参数。该算法求解精度高、鲁棒性好，但需要事先知道目标上的标记点，不适合非合作目标和非标记目标的姿态求解，算法适应能力差。FISHLERMA，FISHLERMA，BOLLESRC.”Randomsampleconsensus:aparadigmformodelfittingwithapplicationstoimageanalysisandautomatedcartography”[J].CommunicationsoftheACM,1981,24(6):381-395.通过提取目标及其投影图像上的大量点对，采用一致交叉验证的方式选择最少的特征点进行姿态的三维重建。该算法需要提取大量的特征点对，运算量大，特征点对存在误匹配时，算法误差大。上述研究成果都对该类问题的特殊情况提出了自己的解决方案，各个方案具有自己的算法特点。但是算法都存在运算量大、精度差、或适应性差等问题。

发明内容

为解决现有方法运算量大、精度差、或适应性差的问题，本发明提供一种在轨三维空间目标姿态估计方法及系统，能有效地从空间目标二维图像中估计出目标的三维空间姿态信息，精度高，计算量小，适应性好。

一种在轨三维空间目标姿态估计方法，包括离线特征库构建步骤和在线姿态估计步骤：

所述离线特征库构建步骤，具体为：

（A1）根据空间目标三维模型获取表征空间目标各种姿态的目标多视点特性视图；

（A2）从各空间目标多视点特性视图中提取几何特征构成几何特征库；所述几何特征包含目标主体高宽比T_i,1、目标纵向对称度T_i,2、目标横向对称度T_i,3和目标主轴倾斜角T_i,4；所述目标主体高宽比T_i,1是指目标最小外接矩形的高宽比；所述目标纵向对称度T_i,2是指在目标最小外接矩形所围成的矩形区域内，目标上半部分面积与下半部分面积之比；所述目标横向对称度T_i,3是指在目标最小外接矩形所围成的矩形区域内，目标左半部分面积与右半部分面积之比；所述目标主轴倾斜角T_i,4是指特性视图的目标柱体主轴线与视图水平方向的夹角；

所述在线姿态估计步骤，具体为：

（B1）对待测在轨空间目标图像预处理；

（B2）对预处理后的待测图像提取特征，该特征与步骤（2）提取的特征相同；

（B3）将对待测图像提取的特征在几何特征库中进行匹配，匹配结果对应的特性视图所表征的空间目标姿态即为待测图像中的目标姿态。

进一步地，所述特征目标主体高宽比T_i,1的提取方式为：

（A2.1.1）对特性视图F_i使用最大类间方差阈值准则得到阈值T_i，将特性视图F_i中大于阈值T_i的像素灰度值f_i(x,y)设置为255，小于或等于阈值T_i的像素灰度值f_i(x,y)设置为零，由此得到二值图像G_i，G_i为宽度n、高度m的像素矩阵，g_i(x,y)为G_i中点（x,y）处像素灰度值；

（A2.1.2）对二值图像G_i按照从上往下、从左往右的顺序进行扫描，若当前点像素值g_i(x,y)等于255，则记录当前像素横坐标x=Topj、纵坐标y=Topi，停止扫描；

（A2.1.3）对二值图像G_i按照从下往上、从左往右的顺序进行扫描，若当前点像素值g_i(x,y)等于255，则记录当前像素横坐标x=Bntj、纵坐标y=Bnti，停止扫描；

（A2.1.4）对二值图像G_i按照从左往右，从上往下的顺序进行扫描，若当前点像素值g_i(x,y)等于255，则记录当前像素横坐标x=Leftj、纵坐标y=Lefti，停止扫描；

（A2.1.5）对二值图像G_i按照从右往左，从上往下的顺序进行扫描，若当前点像素值g_i(x,y)等于255，则记录当前像素横坐标x=Rightj、纵坐标y=Righti，停止扫描；

（A2.1.6）特性视图F_i的目标主体高宽比H_i＝|Topi-Bnti|，W_i＝|Leftj-Rightj|，符号|V|表示变量V的绝对值。

进一步地，所述特征目标纵向对称度T_i,2的提取方式为：

（A2.2.1）计算特性视图F_i的中心点横坐标纵坐标符号表示对变量V取整数部分；

（A2.2.2）统计二值图像G_i中的，1≤横坐标x≤n，1≤纵坐标y≤C_iy的区域内，灰度值为255的像素点个数，即为特性视图F_i的目标上半部分面积ST_i；

（A2.2.3）统计二值图像G_i中，1≤横坐标x≤n，C_iy+1≤纵坐标y≤m的区域内，灰度值为255的像素点个数，即为特性视图F_i的目标下半部分面积SD_i；

（A2.2.4）计算特性视图F_i的目标纵向对称度

进一步地，所述特征目标横向对称度T_i,3的提取方式为：

（A2.3.1）统计二值图像G_i中，1≤横坐标x≤C_ix，1≤纵坐标y≤m的区域内，灰度值为255的像素点个数，即为特性视图F_i的目标左半部分面积SL_i；

（A2.3.2）统计二值图像G_i中，C_ix+1≤横坐标x≤n，1≤纵坐标y≤m的区域内，灰度值为255的像素点个数，即为特性视图F_i的目标右半部分面积SR_i；

（A2.3.3）计算特性视图F_i的目标横向对称度

进一步地，所述特征目标主轴倾斜角T_i,4的提取方式为：

（A2.4.1）计算特性视图F_i对应的二值图像G_i重心横坐标x_i0、纵坐标y_i0：

\{\begin{matrix} x_{i 0} = M_{i} (1,0) / M_{i} (0,0) \\ y_{i 0} = M_{i} (0,1) / M_{i} (0,0) \end{matrix},

式中，

M_{i} (k, j) = Σ_{x = 1}^{n} Σ_{y = 1}^{m} x^{k} y^{j} f_{i} (x, y), k = 0,1, j = 0,1;

（A2.4.2）计算特性视图F_i对应的二值图像G_i对应的p+q阶中心矩μ_i(p,q)：

μ_{i} (p, q) = Σ_{x = 1}^{n} Σ_{y = 1}^{m} {(x - x_{i 0})}^{p} {(y - y_{i 0})}^{q} g_{i} (x, y), p = 0,1,2, q = 0,1,2;

（A2.4.3）构建实对称矩阵

Mat = [\begin{matrix} μ_{i} (2,0), μ_{i} (1,1) \\ μ_{i} (1,1), μ_{i} (0,2) \end{matrix}],

计算矩阵Mat的特征值V₁，V₂，及其对应的特征向量

S_{1} = [\begin{matrix} S_{1 y} \\ S_{1 x} \end{matrix}],

S_{2} = [\begin{matrix} S_{2 y} \\ S_{2 x} \end{matrix}];

（A2.4.4）计算特性视图F_i目标主轴倾斜角T_i4：

T_{i, 4} = \{\begin{matrix} a \tan 2 (| S_{1 x} |, S_{1 y}) * 180 / π, & V_{1} &GreaterEqual; V_{2}, S_{1 x} \leq 0 \\ 180 - a \tan 2 (| S_{1 x} |, S_{1 y}) * 180 / π, & V_{1} &GreaterEqual; V_{2}, S_{1 x} > 0 \end{matrix};

T_{i, 4} = \{\begin{matrix} a \tan 2 (| S_{2 x} |, S_{2 y}) * 180 / π, & V_{1} < V_{2}, S_{1 x} \leq 0 \\ 180 - a \tan 2 (| S_{2 x} |, S_{2 y}) * 180 / π, & V_{1} < V_{2}, S_{1 x} > 0 \end{matrix};

式中，符号π表示圆周率，符号atan2表示反正切函数。

进一步地，还对所述步骤（A2）构建的几何特征库进行归一化处理，对所述步骤（B2）提取的待测图像特征进行归一化处理。

进一步地，所述步骤（A1）根据空间目标三维模型获取表征目标各种姿态的目标多视点特性视图的具体实现方式为；

按俯仰角α每隔角度γ、偏航角β每隔γ将高斯观测球划分为K个二维平面，α＝-180°～0°，β＝-180°～180°，K＝360*180/γ²；

将空间目标三维模型O_T置于高斯观测球球心，从球心将三维模型O_T分别向K个二维平面进行正投影，共得到K个三维模板目标的多视点特性视图F_i；各特性视图F_i为宽度n、高度m的像素矩阵，f_i(x,y)为F_i在点（x,y）处像素灰度值，1≤横坐标x≤n，1≤纵坐标y≤m，i＝1、2、…、K。

进一步地，所述步骤（B1）首先采用非局部均值滤波对待测图像先进行噪声抑制，再选用最大似然估计算法进行去模糊。

进一步地，所述（B3）的具体实现方式为：

（B3.1）遍历整个几何特征库SMF，计算待测图像的四种几何特征{SG₁,SG₂,SG₃,SG₄}与几何特征库SMF中各行向量之间的欧式距离记为D₁、…、D_K，K为目标多视点特性视图的数量；

（B3.2）从各欧式距离D₁、…、D_K中选取四个最小值D_S、D_t、D_u、D_v，对其对应的四个目标姿态求取算术平均，即为待测图像中的目标姿态。

一种在轨三维空间目标姿态估计系统，包括离线特征库构建模块和在线姿态估计模块：

所述离线特征库构建模块，具体包括：

第一子模块，用于根据空间目标三维模型获取表征空间目标各种姿态的目标多视点特性视图；

第二子模块，用于从各空间目标多视点特性视图中提取几何特征构成几何特征库；所述几何特征包含目标主体高宽比T_i,1、目标纵向对称度T_i,2、目标横向对称度T_i,3和目标主轴倾斜角T_i,4；所述目标主体高宽比T_i,1是指目标最小外接矩形的高宽比；所述目标纵向对称度T_i,2是指在目标最小外接矩形所围成的矩形区域内，目标上半部分面积与下半部分面积之比；所述目标横向对称度T_i,3是指在目标最小外接矩形所围成的矩形区域内，目标左半部分面积与右半部分面积之比；所述目标主轴倾斜角T_i,4是指特性视图的目标柱体主轴线与视图水平方向的夹角；

所述在线姿态估计模块，具体包括：

第三子模块，用于对待测在轨空间目标图像预处理；

第四子模块，用于对预处理后的待测图像提取特征，该特征与第二子模块提取的特征相同；

第五子模块，用于将对待测图像提取的特征在几何特征库中进行匹配，匹配结果对应的特性视图所表征的空间目标姿态即为待测图像中的目标姿态。

本发明的技术效果体现在：

本发明步骤（A1）～（A2）为离线训练阶段，利用三维模板目标模型获取目标多视点特性视图，提取特性视图几何特征，进而建立模板目标几何特征库；步骤（B1）～（B2）为待测图像姿态在线估计阶段，将待测图像几何特征与模板目标几何特征库对比，进而估计出待测图像的姿态。本发明匹配所特定使用的几何特征具有尺度不变性，只要在三维建模阶段准确获取目标各部件的相对尺寸比例和位置关系，就可保证后续较高的匹配精度。整个方法实现简单、鲁棒性好、姿态估计精度高，受成像条件影响小，适用性好。

作为优化，对提取的几何特征进行归一化处理，可有效均衡各特征量对姿态估计的影响；对待测图像进行预处理操作，优选非局部均值滤波和最大似然估计算法对待测图像进行去噪、去模糊处理，提高了算法在湍流模糊成像条件下姿态估计精度；对姿态估计结果进行算术加权平均，提高了姿态估计算法的稳定性。

附图说明

图1为姿态估计示意图；

图2为本发明流程示意图；

图3为高斯观测球示意图；

图4为哈勃望远镜三维模型示意图；

图5(a)为哈勃望远镜在俯仰角α=0°、偏航角β=0°情况下投影特性视图；

图5(b)为哈勃望远镜在俯仰角α=0°、偏航角β=90°情况下投影特性视图；

图5(c)为哈勃望远镜在俯仰角α=-90°、偏航角β=90°情况下投影特性视图；

图5(d)为哈勃望远镜在俯仰角α=-180°、偏航角β=90°情况下投影特性视图；

图6(a)为哈勃望远镜某一帧特征视图F_i；

图6(b)为对图6(a)使用最大类间方差准则阈值分割结果；

图6(c)为哈勃望远镜目标高宽比示意图，矩形框ABCD为特性视图F_i的最小外接矩形，|AC|为特性视图F_i的目标主体高度H_i，|CD|为特性视图F_i的目标主体宽度W_i；

图6(d)为哈勃望远镜目标纵向对称度示意图，矩形框abcd所围成的区域为特性视图F_i的目标上半部分，矩形cdef所围成的区域为特性视图F_i的目标下半部分；

图6(e)为哈勃望远镜目标横向对称度示意图，矩形框hukv所围成的区域为特性视图F_i的目标左半部分，矩形ujvl所围成的区域为特性视图F_i的目标右半部分；

图6(f)为哈勃望远镜目标主轴倾斜角示意图，向量为特性视图F_i的目标柱体主轴线，其与水平方向的夹角∠QOR即为目标主轴倾斜角即哈勃望远镜卫星平台的主轴倾斜角；

图7(a)为仿真的哈勃望远镜图像，其对应的俯仰角α、偏航角度β为（α，β）=（-40°,-125°）；

图7(b)为对图7(a)非局部均值滤波结果；

图7(c)为对7(b)最大似然估计算法(MAP)算法校正结果；

图7(d)为对图7(c)姿态估计结果1（α，β）=（-40°,-130°）；

图7(e)为对图7(c)姿态估计结果2（α，β）=（-40°,-140°）；

图7(f)为对图7(c)姿态估计结果3（α，β）=（-40°,-120°）；

图7(g)为对图7(c)为姿态估计结果4（α，β）=（-40°,-150°）；

图7(h)为图7(d)～图7(g)算术平均值结果，作为图7(c)最终的姿态估计结果（α，β）=（-40°,-135°）。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明中在轨三维空间目标为在轨哈勃望远镜，其卫星平台结构为圆柱体；卫星平台上主要搭载了两块矩形太阳能帆板，所需估计的目标姿态是指卫星平台在三维目标坐标系下的姿态。如图1所示为姿态估计示意图，地心坐标系中，X轴指向本初子午线，Z轴指向正北，按照右手螺旋法则确定Y轴方向。目标坐标系中，目标卫星质心始终指向地心，X_s轴平行于地心坐标系中的Y轴，Y_s轴平行于地心坐标系中的Z轴，姿态估计意在从相机坐标系中的目标卫星投影图像中估计出三维目标卫星在目标坐标系下的俯仰角α，即∠N′O_SN、偏航角β即∠N′O_SX_S；其中，O_SN为圆柱体卫星平台轴线；O_SN′为卫星平台轴线O_SN在平面X_SO_SY_S上的投影；相机平面X_mO_mY_m平行于目标坐标系中的X_SO_SY_S平面，同时也平行于地心坐标系中的YOZ平面。

下面以图4所示目标结构为例对本发明做进一步详细说明。以下结合附图和实施例对本发明进一步说明。

本发明流程如图2所示，具体实施方法包括以下步骤。包括：获取模板目标多视点特性视图步骤、建立模板目标几何特征库步骤、计算待测图像几何特征步骤、目标姿态估计步骤；

（A1）获取模板目标多视点特性视图步骤，包括下述子步骤：

（A1.1）建立模板目标三维模型步骤：

对于合作空间目标，如卫星目标，可精确获得卫星平台、卫星搭载的负荷、以及卫星各部件之间的相对位置关系等详细的三维结构及相对位置关系；对于非合作的空间目标，以从目标的多视点投影图像中推断出目标各部件大致的几何结构及相对位置关系。利用目标卫星在轨道上运行时卫星平台质心与地心连线垂直于卫星平台，其太阳能帆板始终指向太阳光入射方向等先验知识，进一步确定卫星各部件之间的空间位置关系。使用三维建模工具MultigenCreator建立目标卫星三维模型；如图4为使用MultigenCreator建立的哈勃望远镜三维模型示意图；

（A1.2）获取模板目标多视点特性视图步骤：

如图3所示，按俯仰角α每隔γ、偏航角β每隔γ将高斯观测球划分为2592个二维平面，α＝-180°～0°，β＝-180°～180°，3°＜γ＜10°，本实例γ=5°；

本发明以哈勃望远镜仿真卫星作为模板目标，如图4所示，将三维模板目标哈勃望远镜O_T置于高斯观测球球心，从球心将三维模板目标O_T分别向2592个二维平面进行正投影，共得到2592个三维模板目标的多视点特性视图F_i。如图5(a)为仿真的哈勃望远镜俯仰角、偏航角（α，β）=（0°，0°）对应的特性视图；图5(b)为（α，β）=（0°，90°）对应的特性视图；图5(c)为（α，β）=（-90°，90°）对应的特性视图；图5(d)为（α，β）=（-180°，90°）对应的特性视图；各特性视图F_i为宽度n=500、高度m=411的像素矩阵，f_i(x,y)为F_i在点（x,y）处像素灰度值，1≤横坐标x≤500，1≤纵坐标y≤411，i＝1、2、…、K，K=2592；

（A2）建立模板目标几何特征库步骤，包括下述子步骤：

以2592帧特性视图中的i=1886帧为例说明本实例：

（A2.1）计算各特性视图F_i的目标主体高宽比T_i,1：

（A2.1.1）对输入的特性视图F_i如图6(a)所示，其对应的俯仰角、偏航角（α，β）=（-50°，-115°）使用最大类间方差阈值准则，得到阈值T_i=95。将像素矩阵F_i中大于95的像素灰度值f_i(x,y)设置为255，小于或等于95的像素灰度值f_i(x,y)设置为零，得到二值图像G_i如图6(b)所示,g_i(x,y)为像素矩阵G_i在点（x,y）处像素灰度值。

（A2.1.2）对二值图像G_i按照从上往下、从左往右的顺序进行扫描，若当前点像素值g_i(x,y)等于255，则记录当前像素横坐标x=Topj、纵坐标y=Topi，停止扫描，本实例中Topj=272，Topi=87。

（A2.1.3）对二值图像G_i按照从下往上、从左往右的顺序进行扫描，若当前点像素值g_i(x,y)等于255，则记录当前像素横坐标x=Bntj、纵坐标y=Bnti，停止扫描，本实例中Bntj=330，Bnti=315。

（A2.1.4）对二值图像G_i按照从左往右，从上往下的顺序进行扫描，若当前点像素值g_i(x,y)等于255，则记录当前像素横坐标x=Leftj、纵坐标y=Lefti，停止扫描，本实例中Leftj=152，Lefti=139。

（A2.1.5）对二值图像G_i按照从右往左，从上往下的顺序进行扫描，若当前点像素值g_i(x,y)等于255，则记录当前像素横坐标x=Rightj、纵坐标y=Righti，停止扫描，本实例中Rightj=361，Righti=282。

（A2.1.6）特性视图F_i的目标主体高宽比定义为目标的高度H_i与目标的宽度W_i的比值H_i＝|Topi-Bnti|，W_i＝|Leftj-Rightj|，符号|V|表示变量V的绝对值。如图6(c)所示，目标主体高宽比T_i,1为目标主体高度AC与目标主体宽度CD的比值。本实例中T_i,1=1.0909,H_i=228，W_i=209。

（A2.2）计算各特性视图F_i的目标纵向对称度T_i,2：

（A2.2.1）计算特性视图F_i的中心点横坐标纵坐标符号表示对变量V取整数部分，本实例中，C_ix=256，C_iy=201。

（A2.2.2）统计二值图像G_i中，1≤横坐标x≤500，1≤纵坐标y≤201的区域内，g_i(x,y)的灰度值为255的像素点个数，即为特性视图F_i的目标上半部分面积ST_i。本实例中，即为图6(d)中矩形框abcd所围成的区域面积，ST_i=10531。

（A2.2.3）统计二值图像G_i中，1≤横坐标x≤500，202≤纵坐标y≤411的区域内，g_i(x,y)的灰度值为255的像素点个数，即为特性视图F_i的目标下半部分面积SD_i。本实例中，即为图6(d)中矩形框cdef所围成的区域面积，SD_i=9685。

（A2.2.4）特性视图F_i的目标纵向对称度

特性视图F_i的目标纵向对称度定义为目标最小外接矩形所围成的矩形区域内，目标上半部分面积ST_i与下半部分面积SD_i之比，本实例中，T_i,2=1.0873。

（A2.3）计算各特性视图F_i的目标横向对称度T_i,3：

（A2.3.1）统计二值图像G_i中，1≤横坐标x≤C_ix，1≤纵坐标y≤m的区域内，g_i(x,y)的灰度值为255的像素点个数，即为特性视图F_i的目标左半部分面积SL_i。本实例中，即为图6(e)中矩形框hukv所围成的区域面积，SL_i=10062。

（A2.3.2）统计二值图像G_i中，C_ix+1≤横坐标x≤n，1≤纵坐标y≤m的区域内，g_i(x,y)的灰度值为255的像素点个数，即为特性视图F_i的目标右半部分面积SR_i。本实例中，即为图6(e)中矩形框ujvl所围成的区域面积，SR_i=10154。

（A2.3.3）计算特性视图F_i的目标横向对称度

特性视图F_i的目标横向对称度定义为目标最小外接矩形所围成的矩形区域内，目标左半部分面积SL_i与右半部分面积SR_i之比。本实例中T_i,3=0.9909。

（A2.4）计算特性视图F_i的目标主轴倾斜角T_i,4：

目标主轴倾斜角定义为特性视图F_i的目标柱体轴线与图像水平方向的夹角θ。该特征最明显地反映了目标的姿态特征，其取值范围为：0°～180°，采用一维浮点数表示。

如图6(f)所示为哈勃望远镜目标主轴倾斜角示意图，向量为特性视图F_i的目标柱体主轴线（本示例中为哈勃望远镜卫星平台主轴线），其与水平方向的夹角∠QOR即为目标主轴倾斜角；

（A2.4.1）计算各特性视图F_i对应的二值图像G_i重心横坐标x_i0、纵坐标y_i0：本实例中，x_i0=252，y_i0=212。

（A2.4.2）计算特性视图F_i对应的二值图像G_i的p+q阶中心矩μ_i(p,q)：

（A2.4.3）构建实对称矩阵

Mat = [\begin{matrix} μ_{i} (2,0), μ_{i} (1,1) \\ μ_{i} (1,1), μ_{i} (0,2) \end{matrix}],

计算矩阵Mat的特征值V₁，V₂，及其对应的特征向量

S_{1} = [\begin{matrix} S_{1 y} \\ S_{1 x} \end{matrix}],

S_{2} = [\begin{matrix} S_{2 y} \\ S_{2 x} \end{matrix}],

本实例中，

Mat = [\begin{matrix} 1.3385 \times 10^{10}, - 8.4494 \times 10^{9} \\ - 8.4494 \times 10^{9}, 1.6366 \times 10^{10} \end{matrix}],

特征值V₁＝6.2955×10⁹，V₂＝2.3455×10¹⁰，特征向量

S_{1} = [\begin{matrix} - 0.7661 \\ - 0.6427 \end{matrix}],

S_{2} = [\begin{matrix} - 0.6427 \\ 0.7661 \end{matrix}] .

（A2.4.4）使用如下公式计算特性视图F_i，如图6(a)目标主轴倾斜角T_i4：

T_{i, 4} = \{\begin{matrix} a \tan 2 (| S_{1 x} |, S_{1 y}) * 180 / π, & V_{1} &GreaterEqual; V_{2}, S_{1 x} \leq 0 \\ 180 - a \tan 2 (| S_{1 x} |, S_{1 y}) * 180 / π, & V_{1} &GreaterEqual; V_{2}, S_{1 x} > 0 \end{matrix};

T_{i, 4} = \{\begin{matrix} a \tan 2 (| S_{2 x} |, S_{2 y}) * 180 / π, & V_{1} < V_{2}, S_{1 x} \leq 0 \\ 180 - a \tan 2 (| S_{2 x} |, S_{2 y}) * 180 / π, & V_{1} < V_{2}, S_{1 x} > 0 \end{matrix};

式中，符号π表示圆周率，符号atan2表示反正切函数。

本实例中目标主轴倾斜角T_i4=50.005°。

（A2.5）构建模板目标多视点特性视图F_i几何特征库MF：

MF = \{\begin{matrix} T_{1,1}, T_{1,2}, T_{1,3}, T_{1,4} \\ T_{2,1}, T_{2,2}, T_{2,3}, T_{2,4} \\ \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \\ T_{i, 1}, T_{i, 2}, T_{i, 3}, T_{i, 4} \\ \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \\ T_{K, 1}, T_{K, 2}, T_{K, 3}, T_{K, 4} \end{matrix}\},

式中，第i行{T_i，1,T_i，2,T_i，3,T_i，4}，表示第i帧特征视图F_i的几何特征；本实例中如图6(a)，{T_i，1,T_i，2,T_i，3,T_i，4}={1.0909,1.0873,0.9909,50.005}。

（A2.6）归一化处理步骤：

对于模板目标多视点特性视图F_i几何特征库MF进行归一化处理，得到模板目标归一化几何特征库SMF：

SMF = \{\begin{matrix} {ST}_{1,1}, {ST}_{1,2}, {ST}_{1,3}, {ST}_{1,4} \\ {ST}_{2,1}, {ST}_{2,2}, {ST}_{2,3}, {ST}_{2,4} \\ \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \\ {ST}_{i, 1}, {ST}_{i, 2}, {ST}_{i, 3}, {ST}_{i, 4} \\ \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \\ {ST}_{K, 1}, {ST}_{K, 2}, {ST}_{K, 3}, {ST}_{K, 4} \end{matrix}\},

式中，Vec_j＝max{T_1，j,T_2，j,…T_i，j…,T_K，j}，i＝1、2、…、K，j=1、2、3、4；符号max{V}表示取集合V的最大值。

在线姿态估计步骤，具体为：

（B1）计算待测图像几何特征步骤，包括下述子步骤：

（B1.1）待测图像预处理步骤

空间目标成像数据噪声很大、信噪比低，而且模糊严重。因此在对成像数据进行后续处理之前，必须首先对成像数据进行预处理，即首先对成像数据进行去噪，随后针对成像数据的特点，利用有效的校正算法对空间目标图像进行图像恢复处理。本实例中，选用非局部均值滤波，选用如下参数：相似性窗口大小为5×5，搜索窗口大小为15×15，衰减参数为15）对待测图像先进行噪声抑制，如图7(a)所示为仿真的哈勃望远镜地基长距离光学成像数据，其对应的俯仰角α、偏航角度β为（α，β）=（-40°,-125°）。如图7(b)为对图7(a)使用非局部均值滤波噪声抑制结果；再选用最大似然估计算法，本实例选用如下参数：外循环次数为8次，估计点扩展函数和目标图像的内循环次数都设置为3次）进行去模糊，得到预处理后的图像g(x,y)，如图7(c)所示为对图7(b)使用最大似然估计算法去模糊结果，即预处理后的图像g(x,y)。

（B2）待测图像几何特征提取步骤

将预处理后的图像g(x,y)代替f_i(x,y)，进行子步骤（2.1）～（2.4），得到待测图像几何特征{G₁,G₂,G₃,G₄}，对几何特征{G₁,G₂,G₃,G₄}进行归一化处理，得到待测图像归一化几何特征{SG₁,SG₂,SG₃,SG₄}：

SG_j＝G_j/Vec_j，j＝1、2、3、4；

（B3）目标姿态估计步骤，包括下述子步骤：

（B3.1）遍历整个模板目标几何特征库SMF，计算待测图像几何特征{SG₁,SG₂,SG₃,SG₄}与SMF中各行向量之间的欧式距离D₁、…、D_K；

（B3.2）从各欧式距离D₁、…、D_K中选取4个最小值D_S、D_t、D_u、D_v，待测图像的姿态设为D_S、D_t、D_u、D_v所代表的模式姿态的算术平均值。如图7(d)～图7(g)所示为D_S、D_t、D_u、D_v所代表的模式姿态，其对应的俯仰角α、偏航角度β分别为（α，β）=（-40°,-130°），（α，β）=（-40°,-140°），（α，β）=（-40°,-120°），（α，β）=（-40°,-150°）。图7(h)为对图7(d)～图7(g)进行算术平均得到的姿态估计结果，（α，β）=（-40°,-135°），即为对图7(a)进行姿态估计的结果。

结果表明俯仰角估计结果的精度误差为零度，偏航角度β估计结果的精度误差为10度。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims

1.一种在轨三维空间目标姿态估计方法，包括离线特征库构建步骤和在线姿态估计步骤：

所述离线特征库构建步骤，具体为：

(A1)根据空间目标三维模型获取表征空间目标各种姿态的目标多视点特性视图；

(A2)从各空间目标多视点特性视图中提取几何特征构成几何特征库；所述几何特征包含目标主体高宽比T_i,1、目标纵向对称度T_i,2、目标横向对称度T_i,3和目标主轴倾斜角T_i,4；所述目标主体高宽比T_i,1是指目标最小外接矩形的高宽比；所述目标纵向对称度T_i,2是指在目标最小外接矩形所围成的矩形区域内，目标上半部分面积与下半部分面积之比；所述目标横向对称度T_i,3是指在目标最小外接矩形所围成的矩形区域内，目标左半部分面积与右半部分面积之比；所述目标主轴倾斜角T_i,4是指特性视图的目标柱体主轴线与视图水平方向的夹角；

所述在线姿态估计步骤，具体为：

(B1)对待测在轨空间目标图像预处理；

(B2)对预处理后的待测图像提取特征，该特征与步骤(A2)提取的特征相同；

(B3)将对待测图像提取的特征在几何特征库中进行匹配，匹配结果对应的特性视图所表征的空间目标姿态即为待测图像中的目标姿态。

2.根据权利要求1所述的在轨三维空间目标姿态估计方法，其特征在于，所述特征目标主体高宽比T_i,1的提取方式为：

(A2.1.1)对特性视图F_i使用最大类间方差阈值准则得到阈值T_i，将特性视图F_i中大于阈值T_i的像素灰度值f_i(x,y)设置为255，小于或等于阈值T_i的像素灰度值f_i(x,y)设置为零，由此得到二值图像G_i，G_i为宽度n、高度m的像素矩阵，g_i(x,y)为G_i中点(x,y)处像素灰度值；

(A2.1.2)对二值图像G_i按照从上往下、从左往右的顺序进行扫描，若当前点像素灰度值g_i(x,y)等于255，则记录当前像素横坐标x＝Topj、纵坐标y＝Topi，停止扫描；

(A2.1.3)对二值图像G_i按照从下往上、从左往右的顺序进行扫描，若当前点像素灰度值g_i(x,y)等于255，则记录当前像素横坐标x＝Bntj、纵坐标y＝Bnti，停止扫描；

(A2.1.4)对二值图像G_i按照从左往右，从上往下的顺序进行扫描，若当前点像素灰度值g_i(x,y)等于255，则记录当前像素横坐标x＝Leftj、纵坐标y＝Lefti，停止扫描；

(A2.1.5)对二值图像G_i按照从右往左，从上往下的顺序进行扫描，若当前点像素灰度值g_i(x,y)等于255，则记录当前像素横坐标x＝Rightj、纵坐标y＝Righti，停止扫描；

(A2.1.6)特性视图F_i的目标主体高宽比H_i＝|Topi-Bnti|，W_i＝|Leftj-Rightj|，符号|V|表示变量V的绝对值。

3.根据权利要求2所述的在轨三维空间目标姿态估计方法，其特征在于，所述特征目标纵向对称度T_i,2的提取方式为：

(A2.2.1)计算特性视图F_i的中心点横坐标纵坐标符号表示对变量V取整数部分；

(A2.2.2)统计二值图像G_i中的，1≤横坐标x≤n，1≤纵坐标y≤C_iy的区域内，灰度值为255的像素点个数，即为特性视图F_i的目标上半部分面积ST_i；

(A2.2.3)统计二值图像G_i中，1≤横坐标x≤n，C_iy+1≤纵坐标y≤m的区域内，灰度值为255的像素点个数，即为特性视图F_i的目标下半部分面积SD_i；

(A2.2.4)计算特性视图F_i的目标纵向对称度

4.根据权利要求3所述的在轨三维空间目标姿态估计方法，其特征在于，所述特征目标横向对称度T_i,3的提取方式为：

(A2.3.1)统计二值图像G_i中，1≤横坐标x≤C_ix，1≤纵坐标y≤m的区域内，灰度值为255的像素点个数，即为特性视图F_i的目标左半部分面积SL_i；

(A2.3.2)统计二值图像G_i中，C_ix+1≤横坐标x≤n，1≤纵坐标y≤m的区域内，灰度值为255的像素点个数，即为特性视图F_i的目标右半部分面积SR_i；

(A2.3.3)计算特性视图F_i的目标横向对称度

5.根据权利要求4所述的在轨三维空间目标姿态估计方法，其特征在于，所述特征目标主轴倾斜角T_i,4的提取方式为：

(A2.4.1)计算特性视图F_i对应的二值图像G_i重心横坐标x_i0、纵坐标y_i0：

\{\begin{matrix} x_{i 0} = M_{i} (1, 0) / M_{i} (0, 0) \\ y_{i 0} = M_{i} (0, 1) / M_{i} (0, 0) \end{matrix},

式中，

M_{i} (k, j) = Σ_{x = 1}^{n} Σ_{y = 1}^{m} x^{k} y^{j} f_{i} (x, y),

k＝0、1，j＝0、1；

(A2.4.2)计算特性视图F_i对应的二值图像G_i对应的p+q阶中心矩μ_i(p,q)：

μ_{i} (p, q) = Σ_{x = 1}^{n} Σ_{y = 1}^{m} {(x - x_{i 0})}^{p} {(y - y_{i 0})}^{q} g_{i} (x, y),

p＝0、1、2，q＝0、1、2；

(A2.4.3)构建实对称矩阵

M a t = [\begin{matrix} μ_{i} (2, 0), μ_{i} (1, 1) \\ μ_{i} (1, 1), μ_{i} (0, 2) \end{matrix}],

计算矩阵Mat的特征值V₁，V₂，及其对应的特征向量

S_{1} = [\begin{matrix} S_{1 y} \\ S_{1 x} \end{matrix}], S_{2} = [\begin{matrix} S_{2 y} \\ S_{2 x} \end{matrix}];

(A2.4.4)计算特性视图F_i目标主轴倾斜角T_i，4：

T_{i, 4} = {\begin{matrix} a t a n 2 (| S_{1 x} |, S_{1 y}) * 180 / π, & V_{1} &GreaterEqual; V_{2}, S_{1 x} \leq 0 \\ 180 - a t a n 2 (| S_{1 x} |, S_{1 y}) * 180 / π, & V_{1} &GreaterEqual; V_{2}, S_{1 x} > 0 \end{matrix};

T_{i, 4} = {\begin{matrix} a t a n 2 (| S_{2 x} |, S_{2 y}) * 180 / π, & V_{1} < V_{2}, S_{2 x} \leq 0 \\ 180 - a t a n 2 (| S_{2 x} |, S_{2 y}) * 180 / π, & V_{1} < V_{2}, S_{2 x} > 0 \end{matrix};

式中，符号π表示圆周率，符号atan2表示反正切函数。

6.根据权利要求1～4任意一项所述的在轨三维空间目标姿态估计方法，其特征在于，还对所述步骤(A2)构建的几何特征库进行归一化处理，对所述步骤(B2)提取的待测图像特征进行归一化处理。

7.根据权利要求1～4任意一项所述的在轨三维空间目标姿态估计方法，所述步骤(A1)根据空间目标三维模型获取表征目标各种姿态的目标多视点特性视图的具体实现方式为；

将空间目标三维模型O_T置于高斯观测球球心，从球心将三维模型O_T分别向K个二维平面进行正投影，共得到K个三维模板目标的多视点特性视图F_i；各特性视图F_i为宽度n、高度m的像素矩阵，f_i(x,y)为F_i在点(x,y)处像素灰度值，1≤横坐标x≤n，1≤纵坐标y≤m，i＝1、2、…、K。

8.根据权利要求1～4任意一项所述的在轨三维空间目标姿态估计方法，所述步骤(B1)首先采用非局部均值滤波对待测图像先进行噪声抑制，再选用最大似然估计算法进行去模糊。

9.根据权利要求1～4任意一项所述的在轨三维空间目标姿态估计方法，所述(B3)的具体实现方式为：

(B3.1)遍历整个几何特征库SMF，计算待测图像的四种几何特征{SG₁,SG₂,SG₃,SG₄}与几何特征库SMF中各行向量之间的欧式距离记为D₁、…、D_K，K为目标多视点特性视图的数量；

(B3.2)从各欧式距离D₁、…、D_K中选取四个最小值D_S、D_t、D_u、D_v，对其对应的四个目标姿态求取算术平均，即为待测图像中的目标姿态。

10.一种在轨三维空间目标姿态估计系统，包括离线特征库构建模块和在线姿态估计模块：

所述离线特征库构建模块，具体包括：

所述在线姿态估计模块，具体包括：

第三子模块，用于对待测在轨空间目标图像预处理；

第五子模块，用于将对待测图像提取的特征在几何特征库中进行匹配，匹配结果对应的特性视图所表征的空间目标姿态即为待测图像中的目标姿。