CN106679634A - 一种基于立体视觉的空间非合作目标位姿测量方法 - Google Patents

Abstract

本发明公开了一种基于立体视觉的空间非合作目标位姿测量方法，无需任何辅助测量工具，也无需知道目标的几何尺寸，直接识别非合作目标自身的固有特征实现对目标位置姿态参数的测量；采用“自适应平滑滤波+Wallis滤波”的图像预处理框架，解决空间复杂环境造成的图像退化问题，提高非合作目标识别与特征提取精度；采用一种基于空间几何约束的特征匹配算法，将像平面内的极线约束推广到空间，在完成特征匹配的同时直接获取其三维信息，降低图像处理算法的复杂度。

Description

一种基于立体视觉的空间非合作目标位姿测量方法

技术领域

本发明属于空间在轨服务技术领域，具体涉及一种基于立体视觉的空间非合作目标位姿测量方法。

背景技术

近年来，人类航天活动不断增加，每年发射的航天器也越来越多。随着时间的消逝，众多航天器中就会产生一些废弃航天器、故障航天器以及空间碎片等，从而影响其他航天器的安全；另外，空间中一些特殊轨道资源是有限的，如地球静止轨道GEO，一旦航天器出现故障或报废，不仅会造成经济损失，同时还会占用相应的轨道资源。因此，以航天器在轨维护、废弃航天器清理、空间攻防等为目的的空间在轨服务技术已成为航天领域的重点研究方向。

根据在轨服务的目标能否提供有效的合作信息，空间在轨服务技术分为合作目标在轨服务和非合作目标在轨服务两类。目前，国内外对非合作目标还没有统一的定义。总结目前在研的空间非合作目标，其非合作特性的表征和程度有所不同，但大都具有以下特点：1)无法利用星间链路来直接传输其位置和姿态信息；2)没有安装用于辅助测量的合作光学标志器和特征块等；3)运动状态、空间结构等信息完全未知或部分未知。从目前的航天技术发展情况来看，空间合作目标的在轨服务技术已相对成熟，并已成功应用于一些航天器在轨维护项目，如日本的ETS-VII、美国的“轨道快车”等；而空间非合作目标的在轨服务技术还不是很成熟，现在还没有国家能够实现空间非合作目标在轨服务的演示验证。

非合作目标在轨服务是一项十分艰巨的任务，摆在面前的难点是如何在目标运动情况和空间结构未知的情况下精确获取其位置和姿态信息。国内外学者针对非合作目标位姿测量问题展开了大量研究，当前主要的非合作目标位姿测量手段有微波雷达测量、激光雷达测量以及光学成像测量等。与其他几种方法比较，基于视觉的光学成像测量以其体积小、质量低、寿命长和稳定性好等优点，已成为逼近阶段非合作目标位姿获取的重要途径。然而，基于视觉的非合作目标位姿测量仍然面临以下问题：1)空间环境比较复杂，获取的非合作目标图像存在一定程度的退化，影响目标识别及特征提取精度。2)图像处理算法复杂度较高，影响算法实时性。

发明内容

有鉴于此，本发明的目的是提供一种基于立体视觉的空间非合作目标位姿测量方法，以有效解决空间复杂环境造成的图像退化问题以及降低图像处理算法的复杂度，实现空间非合作目标位姿的高精度高效测量。

本发明的一种基于立体视觉的空间非合作目标位姿测量方法，采用两个摄像机对空间非合作目标进行拍摄，将两个摄像机分别定义为左摄像机和右摄像机，获得的目标图像分别定义为左图像和右图像，包括如下步骤：

步骤1、对左图像和右图像分别进行预处理，以减少噪声干扰，增强图像特征信息；

步骤2、对其中左图像进行目标特征识别，得到左图像中的特征点；

步骤3、立体匹配与三维重构，具体为：

S301、针对左图像中的任意一个特征点，定义为p_l，连接p_l与左摄像机的光心O，得到光线在该光线上找到特征点p_l在空间中对应的物点，定义为空间点P₀；并得到该空间点在世界坐标系X_wY_wZ_w中的三维坐标空间三维坐标P₀(x₀,y₀,z₀)，再假设空间点P₀的Z_w轴坐标值容差为Δz，在光线上取Z_w轴坐标分别为(z₀-Δz)和(z₀+Δz)的两个空间点P_min和P_max；

S302、在左图像中，选取以特征点p_l为中心的窗口W；在线段上任意选择一个空间点，作为参照空间点，将窗口W投影到右图像上，获得p_l在右图像中的相应的匹配窗口；获得窗口W和对应的匹配窗口之间的相关系数NCC，如下式所示：

其中，m×n表示窗口W的大小，s为位于左图像的窗口W中的像点，s(z_w)为右图像中相应匹配窗口中的像点；I_l和I_r分别表示左图像和右图像中像点的灰度；表示窗口W中的像点的平均灰度值，表示窗口W的匹配窗口内像点的平均灰度值；

S303、按照S302的方法，遍历线段上所有空间点，将每一个空间点分别作为参照空间点，在右图像中分别找到窗口W的匹配窗口，并分别计算相关系数NCC，则相关系数NCC最大时对应的空间点即为特征点p_l的匹配点，同时获得特征点p_l对应的空间点的空间三维坐标；

步骤4、位姿参数解算，具体为：

S401、在左图像中提取的所有目标特征点中任意选取三个不共线的目标特征点，分别定义为P₁,P₂,P₃；

S402、建立目标坐标系，取P₁为目标坐标系的原点，为目标坐标系的y轴方向，特征点P₁,P₂,P₃所在平面的法向量为目标坐标系的z轴方向，根据右手定则确定目标坐标系的x轴方向；

S403、根据步骤S301的方法得到目标特征点P₁,P₂,P₃在世界坐标系中的三维坐标后，则获得目标坐标系与世界坐标系之间的旋转矩阵R与平移向量T；

S404、姿态欧拉角通过旋转矩阵R求出，相对位置由平移向量T给出；由此，得到非合作目标的位置和姿态。

所述图像预处理包括进行自适应平滑滤波与Wallis滤波，其中，自适应平滑滤波算法为迭代算法，设迭代次数为T，则基本迭代步骤为：

(a)以图像I中各像点I(x,y)为中心分别确定窗口Q，计算像点I(x,y)的梯度：

(b)计算窗口权系数w(x,y)：

其中，k为平滑参数，取值为：

k²＝E(I(x,y)-E(I(x,y)))²；

(c)根据以下公式对图像中的像点灰度值进行更新；

(d)根据步骤(c)中更新后的图像I的像点灰度值，采用步骤(a)和(b)的方法，得到梯度和窗口权系数，再利用步骤(c)的方法再次对图像中的像点灰度值进行更新，反复执行步骤(a)、(b)和(c)，直到迭代次数满足要求。

目标特征识别的方法为：利用Canny算法对目标边缘信息进行提取，然后利用Hough变换与最小二乘拟合算法提取目标直线特征和椭圆边界特征，最终结合提取的直线特征与椭圆边界特征实现对目标特征的识别。

本发明具有如下有益效果：

(1)本发明无需任何辅助测量工具，也无需知道目标的几何尺寸，直接识别非合作目标自身的固有特征实现对目标位置姿态参数的测量；(2)采用“自适应平滑滤波+Wallis滤波”的图像预处理框架，解决空间复杂环境造成的图像退化问题，提高非合作目标识别与特征提取精度；(3)采用一种基于空间几何约束的特征匹配算法，将像平面内的极线约束推广到空间，在完成特征匹配的同时直接获取其三维信息，降低图像处理算法的复杂度。

附图说明

图1为本发明基于立体视觉的空间非合作目标位姿测量流程图；

图2为空间几何约束示意图；

图3为双目立体视觉特征三维重构原理图；

图4为本发明基于空间几何约束的图像特征匹配图。

图5为本发明的目标坐标系与世界坐标系相互关系示意图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

基于立体视觉的空间非合作目标位姿测量方法，对左、右摄像机同步采集到的图像进行处理，并计算非合作目标的相对位置姿态参数，主要步骤如下：

步骤1、图像预处理：分别对左、右摄像机图像进行自适应平滑滤波与Wallis滤波，以减少噪声干扰，增强图像特征信息。

自适应平滑滤波算法为迭代算法，设迭代次数为T，则基本迭代步骤为：

(a)以影像I中各像点I(x,y)为中心分别确定窗口Q，计算像点I(x,y)的梯度，如下式所示：

(b)计算窗口权系数w(x,y)

其中，k为平滑参数，k的大小决定了平滑过程中图像细节的保留程度。如果k的取值过大，则图像中的各个边缘也将被平滑掉，这将和传统的加权平均平滑算法没有什么区别；如果k的取值过小，图像中所有的边缘都将被保留，这其中也包括噪声，相当于没有对图像进行平滑。k的取值可用图像的方差来计算：

k²＝E(I(x,y)-E(I(x,y)))²

(c)卷积计算

根据上式对图像中的像点灰度值进行更新；

(d)根据步骤(c)中更新后的图像I的像点灰度值，采用步骤a、b的方法，得到梯度和窗口权系数，利用c的方法再次对图像中的像点灰度值进行更新，反复执行步骤(a)、(b)和(c)，直到迭代次数满足要求。

Wallis滤波的实现步骤如下：

(a)把数字图像分为若干矩形区域，每个矩形区域互不重叠，并且每个矩形区域都有自己要增强的纹理模式尺度；

(b)在各矩形区域中计算其相应的灰度均值与方差；

(c)为各矩形区域灰度均值与方差设定期望的数值，其中方差的值应随着区域尺度的减小而减小，以防止大量像素的灰度值被饱和(即落于[0,255]之外)，然后计算出各矩形区域的Wallis滤波器的乘性系数r₁和加性系数r₀；

(d)由于各矩形区域互不重叠，因此数字图像的任一像素的乘性系数r₁和加性系数r₀均可采用双线性内插得到，并根据下式计算出所有像素新的灰度值：

f(x,y)＝g(x,y)r₁+r₀

r₁＝(cs_f)/(cs_g+s_f/c),r₀＝bm_f+(1-b-r₁)m_g

其中，g(x,y),f(x,y)分别表示滤波前后的图像，参数r₁,r₀分别为乘性系数和加性系数。

步骤2、对左图像进行目标特征识别：利用Canny算法对目标边缘信息进行提取，然后利用Hough变换与最小二乘拟合算法提取目标直线特征和椭圆边界特征，最终结合提取的直线特征与椭圆边界特征实现对目标被测特征如太阳能帆板支架、远地点发动机以及星箭对接环等的识别。

Canny边缘提取算法的步骤如下：(a)用高斯滤波器平滑图像；(b)用一阶偏导的有限差分来计算梯度的幅值和方向；(c)对梯度幅值应用非极大值抑制；(d)用双阈值算法检测和连接边缘。

在利用Canny算法提取出目标边缘信息后，利用Hough变换提取目标直线特征。Hough变换算法的步骤如下：(a)适当地量化参数空间；(b)假定参数空间的每一个单元都是一个累加器，把累加器初始化为零；(c)对图像空间的每一点，在其满足的参数方程对应的累加器上加1；(d)累加器阵列的最大值对应模型的参数。

同时，利用最小二乘拟合算法提取目标椭圆边界特征。最小二乘拟合算法的步骤如下：(a)将椭圆参数的确定转化为多变量约束最优化问题的求解；(b)利用Newton迭代法求解上述最优化问题，确定椭圆最优参数；(c)利用椭圆判别式判断是否为椭圆。

步骤3、立体匹配与三维重构：为了降低图像处理算法的复杂度，在进行目标特征匹配与三维重构时，采用一种基于空间几何约束的特征匹配算法。该算法将像平面内的极线约束推广到空间，放弃传统的基于像平面的匹配策略，采用基于像-物空间关系的匹配策略。在匹配过程中运用空间几何约束条件引导匹配，减少匹配时间，提高匹配的可靠性和精度，并且可在完成匹配的同时获取特征的三维信息，提高效率。

如图2所示，在双目立体视觉中，左图像中的像点pl对应的空间点位于该像点与左摄像机光心连线上，设空间点用P表示，可将空间点P投影到右图像上，得到像点pr，当沿着左图像像点的光线移动空间点，即升降空间点的z轴坐标，同时把空间点投影到右图像上，从而得到一系列像点的轨迹，将其称之为空间投影曲线。空间投影曲线与像平面内的极线有类似的性质，将其称为空间几何约束。

在立体视觉中，特征三维重构原理如图3所示，重构公式如下式所示：

其中，(x_w，y_w，z_w)表示特征点对应的空间点在世界坐标系中的坐标，(X_l,Y_l)表示左图像中特征点的图像坐标，(X_r,Y_r)表示特征点在右图像中的图像坐标，表示摄像机标定外参数旋转矩阵；表示摄像机标定外参数平移向量；f_l和f_r分别表示左、右摄像机焦距。

在左图像中，选取以特征点p_l为中心的窗口W。对于给定的空间点，根据上述空间几何约束，将窗口W投影到右图像上，获得p_l在右图像中的相应的匹配窗口。当移动空间点，则可以在右图像中得到不同的匹配窗口。关于点p_l在左图像I_l和右图像I_r中相对应的匹配窗口之间的相似性测度定义为归一化相关系数NCC，如下式所示：

其中，m×n表示窗口W的大小，窗口W根据匹配精度和匹配时间共同确定。精度要求高，则窗口W取值大，反之，窗口W取值小；当要求匹配时间短时，要求窗口W取值小，反之亦然。s为位于左图像的窗口W中的像点，s(z_w)为右图像中相应匹配窗口中的像点；I_l和I_r分别表示左图像和右图像中像点的灰度；表示窗口W中的像点的平均灰度值，表示窗口W的匹配窗口内像点的平均灰度值。

由归一化相关系数NCC公式可以看出，这里NCC不同于传统的归一化相关系数，它是像点p_l和其空间z_w(z_w∈[z₀-Δz,z₀+Δz])坐标的函数。由此，给定左图像上的特征点p_l和其相应的空间z_w坐标的概略值z₀和Δz，使NCC达到最大值时的z_w正是给定点p_l的正确空间z_w坐标，同时也得到了p_l在右图像上的最佳匹配位置。点p_l在右图像上的搜索范围由初始容差Δz确定。初始容差Δz根据匹配时间共同确定。要求匹配时间短时，要求Δz取值小，反之Δz取值大。为保证匹配的可靠性，在匹配过程中为归一化相关系数NCC设定一定的阈值，只有匹配像对的NCC大于给定的阈值才认为该匹配像对是一对正确的匹配，否则，认为该匹配像对是误匹配，抛弃该匹配像对。

左右摄像机图像构成立体像对I_l-I_r，p_l为左图像I_l上的一个给定特征点，根据透视投影模型可知，光线一定通过p_l在空间中的对应点(O表示左摄像机光心)。对于像点p_l，假设其对应空间点三维坐标中z_w的估计值为z₀，则根据重构公式可得p_l点的空间三维坐标P₀(x₀,y₀,z₀)，再假设p_l点给定z_w值的容差为Δz，可得光线上的两个空间点P_min和P_max，则像点p_l所对应的空间点一定位于线段上。接下来通过相关摄像机标定参数将线段投影到右图像I_r上，得到左图像中点p_l在右图像中的近似极线线段。根据极线的性质，点p_l在右图像I_r中的匹配点一定位于该线段上，如图4所示。然后，利用归一化相关系数NCC公式，使NCC取值最大的点即为点p_l的匹配点，同时获得点p_l的空间z坐标，进而根据重构公式即可得到点p_l相应的空间三维坐标。

步骤4、位姿参数解算：在完成目标特征匹配与三维重构后，在所有提取的目标特征点中任意选取三个不共线的目标特征点P₁,P₂,P₃建立目标坐标系，取P₁为目标坐标系的原点，为目标坐标系的y轴方向，特征点P₁,P₂,P₃所在平面的法向量为目标坐标系的z轴方向，根据右手定则确定目标坐标系的x轴方向。建立的目标坐标系与世界坐标系的关系如图5所示。根据步骤3的方法得到目标特征点P₁,P₂,P₃在世界坐标系中的三维坐标后，则可获得目标坐标系与世界坐标系之间的旋转矩阵R与平移向量T。非合作目标的位置姿态参数可通过旋转矩阵R与平移向量T获得。

设(x_wi,y_wi,z_wi)(i＝1,2,3)表示目标特征点在世界坐标系中的坐标，r_ij(i,j＝1,2,3)表示旋转矩阵R的元素，则旋转矩阵R与平移向量T的计算公式如下：

其中， 与分别为在目标坐标系x轴与y轴上的分量。

姿态欧拉角可通过旋转矩阵R求出，相对位置由平移向量T给出。至此，非合作目标的位置姿态均已测出。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于立体视觉的空间非合作目标位姿测量方法，采用两个摄像机对空间非合作目标进行拍摄，将两个摄像机分别定义为左摄像机和右摄像机，获得的目标图像分别定义为左图像和右图像，其特征在于，包括如下步骤：

步骤1、对左图像和右图像分别进行预处理，以减少噪声干扰，增强图像特征信息；

步骤2、对其中左图像进行目标特征识别，得到左图像中的特征点；

步骤3、立体匹配与三维重构，具体为：

S301、针对左图像中的任意一个特征点，定义为p_l，连接p_l与左摄像机的光心O，得到光线在该光线上找到特征点p_l在空间中对应的物点，定义为空间点P₀；并得到该空间点在世界坐标系X_wY_wZ_w中的三维坐标空间三维坐标P₀(x₀,y₀,z₀)，再假设空间点P₀的Z_w轴坐标值容差为Δz，在光线上取Z_w轴坐标分别为(z₀-Δz)和(z₀+Δz)的两个空间点P_min和P_max；

S302、在左图像中，选取以特征点p_l为中心的窗口W；在线段上任意选择一个空间点，作为参照空间点，将窗口W投影到右图像上，获得p_l在右图像中的相应的匹配窗口；获得窗口W和对应的匹配窗口之间的相关系数NCC，如下式所示：

$NCC(p_l,z_w)=\frac{\underset{s\∈W}{\Σ}\left(I_l\right(s)-{\stackrel{\‾}{I}}_l)\×\left(I_r\right(s\left(z_w\right))-{\stackrel{\‾}{I}}_r)}{\sqrt{\underset{s\∈W}{\Σ}{\left(I_l\right(s)-{\stackrel{\‾}{I}}_l)}^2}\sqrt{\underset{s\∈W}{\Σ}\left(I_r\right(s\left(z_w\right))-{\stackrel{\‾}{I}}_r)}}$

其中，m×n表示窗口W的大小，s为位于左图像的窗口W中的像点，s(z_w)为右图像中相应匹配窗口中的像点；I_l和I_r分别表示左图像和右图像中像点的灰度；表示窗口W中的像点的平均灰度值，表示窗口W的匹配窗口内像点的平均灰度值；

S303、按照S302的方法，遍历线段上所有空间点，将每一个空间点分别作为参照空间点，在右图像中分别找到窗口W的匹配窗口，并分别计算相关系数NCC，则相关系数NCC最大时对应的空间点即为特征点p_l的匹配点，同时获得特征点p_l对应的空间点的空间三维坐标；

步骤4、位姿参数解算，具体为：

S401、在左图像中提取的所有目标特征点中任意选取三个不共线的目标特征点，分别定义为P₁,P₂,P₃；

S402、建立目标坐标系，取P₁为目标坐标系的原点，为目标坐标系的y轴方向，特征点P₁,P₂,P₃所在平面的法向量为目标坐标系的z轴方向，根据右手定则确定目标坐标系的x轴方向；

S403、根据步骤S301的方法得到目标特征点P₁,P₂,P₃在世界坐标系中的三维坐标后，则获得目标坐标系与世界坐标系之间的旋转矩阵R与平移向量T；

S404、姿态欧拉角通过旋转矩阵R求出，相对位置由平移向量T给出；由此，得到非合作目标的位置和姿态。

2.如权利要求1所述的一种基于立体视觉的空间非合作目标位姿测量方法，其特征在于，所述图像预处理包括进行自适应平滑滤波与Wallis滤波，其中，自适应平滑滤波算法为迭代算法，设迭代次数为T，则基本迭代步骤为：

(a)以图像I中各像点I(x,y)为中心分别确定窗口Q，计算像点I(x,y)的梯度：

$\left\{\begin{array}{c}{G}_x(x,y)=I(x+1,y)-I(x-1,y)\\ G_y(x,y)=I(x,y+1)-I(x,y-1)\end{array}\right.$

(b)计算窗口权系数w(x,y)：

$w(x,y)=\exp (-\frac{G_x^2(x,y)+G_y^2(x,y)}{2k^2})$

其中，k为平滑参数，取值为：

k²＝E(I(x,y)-E(I(x,y)))²；

(c)根据以下公式对图像中的像点灰度值进行更新；

$I(x,y)=\frac{\underset{i\∈Q}{\Σ}\underset{j\∈Q}{\Σ}I(x+i,y+j)w(x+i,y+i)}{\underset{i\∈Q}{\Σ}\underset{j\∈Q}{\Σ}w(x+i,y+i)}---\left(6\right)$

(d)根据步骤(c)中更新后的图像I的像点灰度值，采用步骤(a)和(b)的方法，得到梯度和窗口权系数，再利用步骤(c)的方法再次对图像中的像点灰度值进行更新，反复执行步骤(a)、(b)和(c)，直到迭代次数满足要求。

3.如权利要求1所述的一种基于立体视觉的空间非合作目标位姿测量方法，其特征在于，目标特征识别的方法为：利用Canny算法对目标边缘信息进行提取，然后利用Hough变换与最小二乘拟合算法提取目标直线特征和椭圆边界特征，最终结合提取的直线特征与椭圆边界特征实现对目标特征的识别。