CN101435732B - 一种基于双目光流的空间目标旋转轴及质心估计方法 - Google Patents

Abstract

一种基于双目光流的空间目标旋转轴及质心估计方法，包括①分别对左右视图进行特征点的检测；②按时序依次隔若干帧计算光流；③根据特征点的双目光流进行单位时间内三维运动矢量的重建；④根据每个特征点的三维运动矢量序列对旋转轴进行单次估计，并计算出该特征点运动轨迹半径；⑤将各轨迹平面的法矢量以及轨迹圆心进行加权平均，得到估计的旋转轴的空间方程；⑥待目标的旋转轴改变后，重复步骤①～⑤，得到第二条旋转轴的空间方程，求解出目标的质心空间坐标。随着观测卫星与空间目标的逐步逼近，本发明的估计精度逐渐增加，错误逐渐减少，可靠性逐渐增加，实时性与准确性均得到保证；可运用于在未知环境下相对于多面体空间物体的自主图像导航。

Description

一种基于双目光流的空间目标旋转轴及质心估计方法

技术领域

本发明属于数字图像处理与空间科学技术交叉的技术领域，具体涉及一种基于双目光流的空间目标旋转轴及质心估计方法。该方法通过序列图像处理和双目光流计算，来估计空间目标的旋转轴及质心，可应用于观测卫星对在太空中处于自旋运动或姿态调节过程中的空间物体的特性测量及观测卫星与空间物体在逐渐靠近过程中的自主图像导航。

背景技术

视觉运动分析是研究从场景的图像序列中提取场景中目标物体的结构、位置和运动信息的理论和方法。这项技术被广泛的应用于机器人导航、汽车智能交通、非接触物体测量等领域。

在未知环境中，自主导航系统面临的技术困难已不仅仅是规避静止物体，而且需要根据环境的变化来跟踪、观测、规避运动的物体。比如机器人导航、太空卫星跟踪编队等。很多场合需要视觉系统自动地跟踪和观测未知物体。对未知目标进行近距离的特性观测和相对导航，需要处理的核心问题就是估计视觉系统与目标物体之间的相对位置和相对姿态，建立相对导航的运动方程。为此，首先就要估计目标物体的质心或旋转轴。目前的光流技术通常是基于单目观测的，即分析单个成像传感器上的光流来估计物体的速度、距离。如Yair Barniv在文献“passive ranging using imageexpansion”(Yair Barniv，IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONICSYSTEMS VOL.31，NO.1 JANUARY 1995)中通过光流的变化来估计立方体由远及近的速度和距离。该文介绍了几种光流检测的方法，它只简单分析了利用单目光流场来估计立方体相对相机的接近速度和相对距离，并没有分析其他复杂的运动模型，如刚体的旋转，也没有结合双目立体视觉来进行相关研究。

由于在进行卫星编队飞行或者交会时需要对目标的质心进行准确定位，因此准确估计卫星的质心是非常必要的。目前国内外对空间目标进行质心估计的研究还停留在针对合作目标(能够与目标进行实时通讯交流)和具备先验知识的目标(已知目标的三维特征或其他信息)阶段。对于非合作目标和未知目标的质心估计，国内外还没有相关文献报道。而本发明中介绍的方法可适用于非合作的未知目标。

发明内容

本发明的目的是提供一种基于双目光流的空间目标旋转轴及质心估计方法，该方法可有效估计空间自旋目标的旋转轴和质心，从而解决观察卫星与目标的相对导航与定位问题。

本发明提供的基于双目光流的空间目标旋转轴及质心估计方法，其步骤包括：

(1)对左、右视图进行特征点检测，左、右视图分别由左、右相机同时拍摄，并且包括至少五个时刻的图像对；

(2)相隔至少一帧按时序依次分别计算左、右视图的光流，得到特征点的双目光流。

(3)根据特征点的双目光流进行单位时间内三维运动矢量的重建。

(4)由单个特征点的三维运动矢量序列对旋转轴进行单次估计，并计算出该特征点运动轨迹半径。

(5)根据各特征点拟合轨迹的半径，将各轨迹平面的法矢量以及轨迹圆心进行加权平均，得到最终的旋转轴的空间方程。

(6)在另外的时间段，当目标的旋转轴改变后，重复步骤(1)～(5)，得到另一条旋转轴的空间方程，与前一次求出的旋转轴的方程联立，求解出目标的质心空间坐标。

根据运动力学原理，空间自旋目标的质心必然在其旋转轴上，针对这一物理特性，本发明通过目标特征点在两个成像传感器获取的序列图像上形成的光流，估计其运动轨迹，进而估计目标的旋转轴，通过估计不同时段的多条旋转轴，从而估计目标的质心。本发明的优点在于不需要目标物体的先验信息，也不需要在目标物体上设置星标，而是通过两个成像传感器(如双目可见光相机或红外相机)获取的序列图像，根据物体运动特征点在两个成像面上形成的光流来估计目标物体的旋转轴和质心。由该方法估计的自旋运动目标(包括调姿过程中的短暂自旋运动)质心的相对误差小于10％。

附图说明

图1是本发明方法的处理流程图；

图2是对同一特征点在运动过程中逐帧计算光流；

图3是对同一特征点在运动过程中采取隔n-1帧计算光流；

图4是三维运动矢量到双目成像系统的投影模型；

图5是特征点圆弧轨迹示意图；

图6是估计的特征点空间运动轨迹有拐点时，对旋转轴估计的影响；

图7是目标在不同时刻的三个旋转轴的估计(L₁、L₂、L₃分别为目标在T₁、T₂、T₃时刻的实际旋转轴，l₁、l₂、l₃分别为相应的估计旋转轴)；

图8是T₁时间段，检测得到的4个特征点，L为旋转轴；

图9是T₁时间段，隔1帧计算的4个特征点在观测时间内的光流轨迹；

图10是目标绕第一个轴旋转时，最终估计的旋转轴L’(虚线)与实际旋转轴L(实线)比较；

图11是T₁时间段，对目标特征点B以相邻帧计算光流后，估计的旋转轴L’(虚线)与实际旋转轴L(实线)的比较；

图12是T₁时间段，对目标特征点B按时序依次隔一帧计算光流后，估计的旋转轴L’(虚线)与实际旋转轴L(实线)的比较；

图13是T₂时间段，检测得到的4个特征点，L为旋转轴；

图14是目标绕第二个轴旋转时，最终估计的旋转轴L’(虚线)与实际旋转轴L(实线)比较；

图15是T₃时间段，检测得到的4个特征点，L为旋转轴；

图16是目标绕第三个轴旋转时，最终估计的旋转轴L’(虚线)与实际旋轴L(实线)比较；

图17是表示三条估计轴与实际轴的误差角度随计算光流所隔帧数的变化(L₁、L₂、L₃分别对应T₁、T₂、T₃三个时刻)；

图18是表示三条估计轴与实际轴的误差距离随计算光流所隔帧数的变化(L₁、L₂、L₃分别对应T₁、T₂、T₃三个时刻)。

具体实施方式

下面结合附图和实例对本发明作进一步详细的说明。

本发明利用两个可见光相机形成双目立体视觉系统拍摄自转的目标，目标的形殊点(如角点)会在两个相机的像面上形成三维运动轨迹的二维投影，这些三维运动轨迹在像面上的二维投影称为“光流”。同一个点在两个相机上形成的光流称作光流立体像对，本发明的主要方法就是根据自转目标上多个形殊点在两相机上的光流立体像对，估计出目标的旋转轴，由于物体自转的旋转轴经过物体的质心，所以当目标在不同时段绕不同旋转轴旋转时，就可以根据两次(或多次)旋转的旋转轴估计出物体的质心位置。如图1所示，本发明方法包括以下步骤：

(1)对左、右视图进行特征点检测，左、右视图分别由左、右相机同时拍摄，并且包括至少五个时刻的图像对；

通常目标表面结构丰富，特别是多面体构成的卫星物体，完全可以选择符合条件的特征点来进行光流检测。可以用Canny或其它边缘算子提取目标的边缘信息，再根据Hough变换检测直线段并计算相邻直线段的交点的方法，检测出目标的角点，另外，可利用Harris角点检测或者SUSAN角点检测法辅助检测目标角点，以提高角点检测的精度。

(2)相隔至少一帧按时序依次分别计算左、右视图的光流，得到特征点的双目光流。

计算光流的图像对可以是非相邻帧的。如图2所示，相邻帧计算光流得到的每个光流较小，这使得后面重建三维运动矢量的误差较大。因此，为了提高重建三维运动矢量的精度，本发明推出了按时序依次隔至少一帧计算特征点光流的方法。如图3所示，图中的序号表示时刻数，各序号所对应的点表示特征点在该时刻在图像中的位置。那么，为了得到较大的光流，可以延迟n-1帧计算光流，这样，以第1帧为起点，第n帧为终点可得到一条光流，依此类推，第2帧为起点，第n+1为终点，可得到第二条光流……

在具体操作时，可根据实际情况调整相隔的帧数，若目标在图像中的区域较小，那么可以适当增大相隔的帧数。

(3)根据特征点的双目光流进行单位时间内三维运动矢量的重建。

三维运动矢量到双目成像系统的投影模型如图4所示，在该模型中，一立方体绕旋转轴做旋转运动，实线所示的立方体为旋转之前的状态，虚线所示的立方体是旋转之后的状态，旋转之前的P点对应于旋转之后的P’点，

即为立方体上某一点的三维运动矢量。

O₁X₁Y₁Z₁为左相机坐标系，O₂X₂Y₂Z₂为右相机坐标系，O₁′U₁V₁为左相机的成像平面坐标系，O₁′U₁V₁平面与O₁X₁Y₁平面平行，O₂′U₂V₂为右相机的成像平面坐标系，O₂′U₂V₂平面与O₂X₂Y₂平面平行。两相机均采用针孔模型，而O₁和O₂分别为左右相机的“针孔”，O₁O₁′和O₂O₂′分别为左右相机的焦距，大小都为f。两相机的光轴O₁Z₁与O₂Z₂平行，且基线距离O₁O₂大小为b₀。

不失一般性，为了计算的简便，不妨定义空间坐标系与左相机坐标系O₁X₁Y₁Z₁重合。在该定义下，左相机坐标系相对空间坐标系不发生旋转或平移运动，而右相机坐标系相对空间坐标系只发生水平方向的偏移，偏移量即为两相机的基线距离b₀。

P点在左相机成像平面的投影为P₁，P′点在左相机成像平面的投影为P₁′，那么三维运动矢量在左相机成像平面投影形成的光流即为

，同理，

在右相机成像平面投影形成的光流为图4中的

在左相机成像平面坐标系O₁′U₁V₁中，令P₁的坐标为(U₁，V₁)，P₁′的坐标为(U₁′，V₁′)；在右相机成像平面坐标系O₂′U₂V₂中，令P₂的坐标为(U₂，V₂)，P₂′的坐标为(U₂′，V₂′)；在世界坐标系中，令P点的坐标为(X，Y，Z)，P′点的坐标为(X′，Y′，Z′)。

由相机针孔模型可得：

$\left\{\begin{array}{c}{U}_1=X\frac{f}{Z}\\ V_1=Y\frac{f}{Z}\end{array}\right.,\left\{\begin{array}{c}{U_1}^{\′}=X^{\′}\frac{f}{Z^{\′}}\\ {V_1}^{\′}=Y^{\′}\frac{f}{Z^{\′}}\end{array}\right.,---\left(1\right)$

$\left\{\begin{array}{c}{U}_2=(X-b_0)\frac{f}{Z}\\ V_2=Y\frac{f}{Z}\end{array}\right.,\begin{array}{}\end{array}\left\{\begin{array}{c}{U_2}^{\′}=(X^{\′}-b_0)\frac{f}{Z^{\′}}\\ {V_2}^{\′}=Y^{\′}\frac{f}{Z^{\′}}\end{array}\right..---\left(2\right)$

O₁′U₁V₁坐标系中，所在直线的方程为：

$V=\frac{V_1-{V_1}^{\′}}{U_1-{U_1}^{\′}}(U-U_1)+V_1.---\left(3\right)$

O₂′U₂V₂坐标系中，

所在直线的方程为：

$V=\frac{V_2-{V_2}^{\′}}{U_2-{U_2}^{\′}}(U-U_2)+V_2.---\left(4\right)$

根据图4所示的三维运动矢量到双目成像系统的投影模型，可以由三维运动矢量在两个相机成像平面上投影形成的光流，通过直线重建的方法，估计出三维运动矢量所在空间直线的方程。

空间三维运动矢量

所在直线在空间坐标系中的方程可以写为：

$\left\{\begin{array}{c}x=\mathrm{az}+p\\ y=\mathrm{bz}+q\end{array}\right.,(q\≠0)---\left(5\right)$

设左相机成像平面上的光流

所在的直线为l₁，右相机成像平面上的光流

所在的直线为l₂，根据式(5)以及式(3)、(4)l₁，l₂分别在各自成像平面坐标系的方程可以写为：

l₁：qU₁-pV₁＝f(aq-bp)，(6)

l₂：qU₂-(p+b₀)V₂＝f[aq-b(p+b₀)]，(7)

其中b₀和f分别为两个相机之间的基线长度和焦距。

由于光流可直接在图像序列中检测得到，因此光流所在直线的斜率和截距也可以计算得到，设l₁的斜率为k₁，截距为m₁，l₂的斜率为k₂，截距为m₂(这里斜率的定义为 $k=\frac{\mathrm{\Δu}}{\mathrm{\Δv}}$ )分别为：

l₁：U₁-k₁V₁＝m₁，(8)

l₂：U₂-k₂V₂＝m₂.(9)

联立式(5)～式(7)可得：

$\left\{\begin{array}{c}{k}_1=p/q,\\ k_2=(q+b_0)/q,\\ m_1=f(\mathrm{aq}-\mathrm{bp})/q,\\ m_2=f\[\mathrm{aq}-b(p+b_0)\]/q.\end{array}\right.---\left(10\right)$

由(10)解出a，b，p，q，得到直线方程(适用于q≠0情况)：

$\left\{\begin{array}{c}a=\frac{m_1}{f}+\frac{k_1(m_1-m_2)}{f(k_2-k_1)},\\ b=\frac{m_1-m_2}{f(k_2-k_1)},\\ p=\frac{k_1{b}_0}{k_2-k_1},\\ q=\frac{b_0}{k_2-k_1}.\end{array}\right.---\left(11\right)$

求出了a，b，p，q后，代入方程(5)，就可求出三维运动矢量

所在的直线在空间坐标系中的方程。

当q＝0，p≠0时，空间直线在两相机成像平面的投影相同：

l₁：V₁＝fb

l₂：V₂＝fb

由于两幅图像在U方向没有差异(即不存在视差)，无法重建直线。

当p＝0，q＝0时(即空间直线经过坐标原点)，由于针孔模型为仿射投影，直线在图像上的投影为一个点，此时将不能重建空间直线。

(4)由单个特征点的三维运动矢量序列对旋转轴进行单次估计，并计算出该特征点运动轨迹半径。

如图5所示，目标上某点绕旋转轴在单位时间t内由A点运动到B点，圆弧轨迹的半径为R，旋转的角度为θ，线速度为

，单位时间t内，由A点到B点的平均速度为

与

的矢量夹角为。那么

与的大小关系就如式(12)所示：令目标旋转的角速度为 $\stackrel{\→}{\mathrm{\ω}}={({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,{\mathrm{\ω}}_z)}^T,$ 且 $\mathrm{\Ω}=\left|\stackrel{\→}{\mathrm{\ω}}\right|,$

因目标旋转的角速度为

，且 $\mathrm{\Ω}=\left|\stackrel{\→}{\mathrm{\ω}}\right|,$ 那么由式(12)可得：

$\frac{\left|\stackrel{\→}{v}\right|}{\left|\stackrel{\→}{v_0}\right|}=\frac{\sin \left(\frac{1}{2}\mathrm{\Ωt}\right)}{\frac{1}{2}\mathrm{\Ωt}}---\left(13\right)$

以上分析了与

的模的关系，矢量

可看作是矢量旋转

后再进行模的放缩得到，令R_t为

旋转到

的旋转矩阵，又令矩阵

的表达式如下：

$\widehat{{\mathrm{\ω}}^{\′}}=\left[\begin{array}{ccc}0& -{{\mathrm{\ω}}_z}^{\′}& {{\mathrm{\ω}}_y}^{\′}\\ {{\mathrm{\ω}}_z}^{\′}& 0& -{{\mathrm{\ω}}_x}^{\′}\\ -{{\mathrm{\ω}}_y}^{\′}& {{\mathrm{\ω}}_x}^{\′}& 0\end{array}\right]---\left(14\right)$

其中， $\stackrel{\→}{{\mathrm{\ω}}^{\′}}={({{\mathrm{\ω}}_x}^{\′},{{\mathrm{\ω}}_y}^{\′},{{\mathrm{\ω}}_z}^{\′})}^T,$

为

的单位向量，那么

$R_t=I+\sin \left(\frac{1}{2}\mathrm{\Ωt}\right)\·\widehat{{\mathrm{\ω}}^{\′}}+\[1-\cos \left(\frac{1}{2}\mathrm{\Ωt}\right)\]\·{\widehat{{\mathrm{\ω}}^{\′}}}^2.$

根据以上分析，可得到

与

的关系如下：

$\stackrel{\→}{v}=\frac{\sin \left(\frac{1}{2}\mathrm{\Ωt}\right)}{\frac{1}{2}\mathrm{\Ωt}}{R}_t\stackrel{\→}{v_0}---\left(15\right)$

又由于 $\stackrel{\→}{v_0}=\stackrel{\→}{\mathrm{\ω}}\×\stackrel{\→}{\mathrm{OA}}=\mathrm{\Ω}\·\widehat{{\mathrm{\ω}}^{\′}}\·\stackrel{\→}{\mathrm{OA}},$ 因此可得：

$\stackrel{\→}{\mathrm{AB}}=\stackrel{\→}{v}\·t=\frac{\sin \left(\frac{1}{2}\mathrm{\Ωt}\right)}{\frac{1}{2}\mathrm{\Ωt}}{R}_t\·\stackrel{\→}{v_0}\·t=\{\sin \left(\mathrm{\Ωt}\right)\·\widehat{{\mathrm{\ω}}^{\′}+}\[1-\cos \left(\mathrm{\Ωt}\right)\]\·{\widehat{{\mathrm{\ω}}^{\′}}}^2\}\·\stackrel{\→}{\mathrm{OA}}---\left(16\right)$

令矩阵 $M=\sin \left(\mathrm{\Ωt}\right)\·\widehat{{\mathrm{\ω}}^{\′}}+\[1-\cos \left(\mathrm{\Ωt}\right)\]\·{\widehat{{\mathrm{\ω}}^{\′}}}^2,$ 展开后可得：

$M=\left[\begin{array}{cccc}-(1-\cos \mathrm{\Ωt})(1-{{\mathrm{\ω}}_x}^{\′2})& & (1-\cos \mathrm{\Ωt}){{\mathrm{\ω}}_x}^{\′}{{\mathrm{\ω}}_y}^{\′}-{{\mathrm{\ω}}_z}^{\′}\sin \mathrm{\Ωt}& (1-\cos \mathrm{\Ωt}){{\mathrm{\ω}}_x}^{\′}{{\mathrm{\ω}}_z}^{\′}+{{\mathrm{\ω}}_y}^{\′}\sin \mathrm{\Ωt}\\ (1-\cos \mathrm{\Ωt}){{\mathrm{\ω}}_x}^{\′}{{\mathrm{\ω}}_y}^{\′}+{{\mathrm{\ω}}_z}^{\′}\sin \mathrm{\Ωt}& & -(1-\cos \mathrm{\Ωt})(1-{{\mathrm{\ω}}_y}^{\′2})& (1-\cos \mathrm{\Ωt}){{\mathrm{\ω}}_y}^{\′}{{\mathrm{\ω}}_z}^{\′}-{{\mathrm{\ω}}_x}^{\′}\sin \mathrm{\Ωt}\\ (1-\cos \mathrm{\Ωt}){{\mathrm{\ω}}_x}^{\′}{{\mathrm{\ω}}_z}^{\′}-{{\mathrm{\ω}}_y}^{\′}\sin \mathrm{\Ωt}& & (1-\cos \mathrm{\Ωt}){{\mathrm{\ω}}_y}^{\′}{{\mathrm{\ω}}_z}^{\′}+{{\mathrm{\ω}}_x}^{\′}\sin \mathrm{\Ωt}& -(1-\cos \mathrm{\Ωt})(1-{{\mathrm{\ω}}_z}^{\′2})\end{array}\right]---\left(17\right)$

令A点的坐标为(x_A，y_A，z_A)^T，O点的坐标为(x_O，y_O，z_O)^T，设 $\stackrel{\→}{\mathrm{AB}}={(x,y,z)}^T,$ 向量可由双目立体光流场重建得到，而A点的坐标为向量

的起始点，也可由双目立体光流场重建得到，那么由式(16)可得到：

$M\·\left[\begin{array}{c}{x}_A-x_O\\ y_A-y_O\\ z_A-z_O\end{array}\right]=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]---\left(18\right)$

将参数作如下变量代换：

Q₁＝(1-cosΩt)(1-ω_x′²)，

Q₂＝(1-cosΩt)(1-ω_y′²)，

Q₃＝(1-cosΩt)(1-ω_z′²)

Q₄＝(1-cosΩt)ω_x′ω_y′，

Q₅＝(1-cosΩt)ω_x′ω_z′，

Q₆＝(1-cosΩt)ω_y′ω_z′，

Q₇＝ω_x′sinΩt，

Q₈＝ω_y′sinΩt，

Q₉＝ω_z′sinΩ，

Q₁₀＝x₀Q₁-y₀Q₄-z₀Q₅-z₀Q₈+y₀Q₉，

Q₁₁＝y₀Q₂-x₀Q₄-z₀Q₆+z₀Q₇-x₀Q₉，

Q₁₂＝z₀Q₃-x₀Q₅-z₀Q₆-y₀Q₇+x₀Q₈.

可得到：

$\left[\begin{array}{cccccccccccc}-x_A& 0& 0& y_A& z_A& 0& 0& z_A& -y_A& 1& 0& 0\\ 0& -y_A& 0& x_A& 0& z_A& -z_A& 0& x_A& 0& 1& 0\\ 0& 0& -z_A& 0& x_A& y_A& y_A& -x_A& 0& 0& 0& 1\end{array}\right]\·\begin{array}{}\end{array}\left[\begin{array}{c}{Q}_1\\ \·\\ \·\\ \·\\ Q_{12}\end{array}\right]=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]---\left(19\right)$

如果只观测特征点在一个时间段的空间位置变化(由A运动到B)，可以得到上述3个方程，若观测该特征点在多个时间段的运动变化(至少四个时间段)，使方程的系数矩阵的行数大于12，即可以解出Q₁到Q₁₂。

由于：

$\left[\begin{array}{ccc}{Q}_1& Q_9-Q_4& -(Q_5+Q_8)\\ -(Q_9+Q_4)& Q_2& Q_7-Q_6\\ Q_8-Q_5& -(Q_6+Q_7)& Q_3\end{array}\right]\·\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]=\begin{array}{}\end{array}\left[\begin{array}{c}{Q}_{10}\\ Q_{11}\\ Q_{12}\end{array}\right]---\left(20\right)$

所以可解出O点的坐标(x_O，y_O，z_O)^T。

而ω_x′∶ω_y′∶ω_z′＝Q₇∶Q₈∶Q₉.(21)

又由于 $\stackrel{\→}{{\mathrm{\ω}}^{\′}}={({{\mathrm{\ω}}_x}^{\′},{{\mathrm{\ω}}_y}^{\′},{{\mathrm{\ω}}_z}^{\′})}^T$ 是个单位矢量，所以有ω_x′²+ω_y′²+ω_z′²＝1，结合式(21)可以解出(ω_x′，ω_y′，ω_z′)^T。由于旋转轴的方向与

的方向相同，又解出了该旋转轴经过点(x_O，y_O，z_O)^T，所以空间坐标系中，旋转轴的方程为：

$\frac{X-x_0}{{{\mathrm{\ω}}_x}^{\′}}=\frac{Y-y_0}{{{\mathrm{\ω}}_y}^{\′}}=\frac{Z-z_0}{{{\mathrm{\ω}}_z}^{\′}}.---\left(22\right)$

另外，还可以根据 $\mathrm{\Ω}=\frac{1}{t}\arcsin \frac{Q_7}{{{\mathrm{\ω}}_x}^{\′}}=\frac{1}{t}\arcsin \frac{Q_8}{{{\mathrm{\ω}}_y}^{\′}}=\frac{1}{t}\arcsin \frac{Q_9}{{{\mathrm{\ω}}_z}^{\′}}$ 求出物体旋转的角速度的大小。

由于已求出Ω和(ω_x′，ω_y′，ω_z′)^T，所以根据式(17)可求出式(18)中的矩阵M，而由式(20)已求出(x_O，y_O，z_O)^T，因此可根据式(18)得到：

$\left[\begin{array}{c}{x}_A\\ y_A\\ z_A\end{array}\right]=M^{-1\·}\left[\begin{array}{c}x\\ y\\ z\end{array}\right]+\left[\begin{array}{c}{x}_O\\ y_O\\ z_O\end{array}\right]---\left(23\right)$

所以可求得圆弧轨迹的半径为：

$R=\sqrt{{(x_A-x_O)}^2+{(y_A-y_O)}^2+{(z_A-z_O)}^2}---\left(24\right)$

根据上述原理，我们基于目标绕同一个旋转轴旋转的左右视图序列，获得目标上某一点在一系列时刻的离散运动场投影到左右相机的光流，并且对这一系列光流进行重建，获得重建后的三维运动矢量，对每个特征点检测4组(或4组以上)的光流立体像对，将重建后的三维运动矢量代入式(18)，可分别得到由不同的特征点的运动轨迹估计的旋转轴的方向，以及轨迹的半径和圆心。

某些时刻，当三维运动矢量在两相机成像面投影得到的光流大小远小于三维运动矢量本身的大小时(即三维运动矢量与成像平面的角度较大时)，由于光流的双目视差较小，所以在重建三维运动矢量时，深度信息的损失较大，所以重建误差也较大，因此在重建三维运动矢量的空间轨迹时，就会表现出在某些时刻重建的轨迹中会出现明显的拐点(如图6中A、B两处所示)。

基于以上分析，在单次估计过程中应特别注意：在对旋转轴进行单次估计之前，应首先估计重建的三维运动矢量序列的平滑度，若不对拐点处进行平滑处理，会严重影响旋转轴的估计精度。图6中，L为实际旋转轴，L’为估计的旋转轴，P点为估计错误的圆弧轨迹的圆心。

(5)根据各特征点的拟合轨迹的半径，将各轨迹平面的法矢量以及轨迹圆心进行加权平均，得到最终的旋转轴的空间方程。

由一个点的运动轨迹，可以确定该点的圆弧轨迹的圆心以及垂直于该圆平面的法线方向，而在理论上，而该物体旋转轴的方向和各个点的圆弧轨迹所在平面的法线方向相同，且经过圆弧的圆心，所以将不同的点的运动轨迹所确定的旋转轴的方向和圆心的坐标加权平均，就可以最终估计出一个误差较小的旋转轴的方向，和旋转轴所经过的一点(该点就是各个圆弧轨迹圆心加权平均的结果)。由于旋转半径大的点所形成的光流的检测误差相对较小，因此在加权平均时，半径较大的运动轨迹所确定的圆心坐标和法线方向所占的权重应该较大。

令每个圆弧轨迹的圆心坐标为(x_i，y_i，z_i)，半径为R_i，轨迹拟合平面的法线方向为(a_i，b_i，c_i)，那么最终估计的圆心坐标(x，y，z)和法线方向{a，b，c}的计算公式如下所示：

$\stackrel{\‾}{x}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{x}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right),\stackrel{\‾}{y}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{y}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right),\stackrel{\‾}{z}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{z}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right)---\left(25\right)$

$\stackrel{\‾}{a}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{a}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right),\stackrel{\‾}{b}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{b}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right),\stackrel{\‾}{c}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{c}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right)---\left(26\right)$

因此最终估计的旋转轴的方程为：

$\frac{x-\stackrel{\‾}{x}}{\stackrel{\‾}{a}}=\frac{y-\stackrel{\‾}{y}}{\stackrel{\‾}{b}}=\frac{z-\stackrel{\‾}{z}}{\stackrel{\‾}{c}}---\left(27\right)$

(6)在另外的时间段，当目标的旋转轴改变后，重复步骤(1)～(5)，得到另一条旋转轴的空间方程，与前一次求出的旋转轴的方程联立，求解出目标的质心空间坐标。

实例：

下面结合仿真实例来详细说明本发明方法的正确性和有效性。

假设仿真卫星在T₁、T₂、T₃三个时段分别绕三个不同的旋转轴L₁、L₂、L₃旋转。图7中给出了三个旋转轴在空间中的相对位置，其中l₁、l₂、l₃是用本方法估计的仿真卫星在三个不同时段的旋转轴。

在T₁时段，仿真卫星与相机相距20m，相机的视场角为10度，卫星的旋转轴为图8中的L，检测得到的特征点为图8中的A、B、C、D。

在T₂时段，仿真卫星与相机相距20m，相机的视场角为20度，卫星的旋转轴为图13中的L，检测得到的特征点如图13中的A、B、C、D。

在T₃时段，仿真卫星与相机相距20m，相机的视场角为20度，卫星的旋转轴为图15中的L，检测得到的特征点如图21中的A、B、C、D。

为简化运算，设帧间的时间间隔为1秒，在以上三个时段中，卫星的自转角速度均为15度/秒，即每帧相对于前一帧物体在世界坐标系中绕轴旋转了15度。

图9(a)-(d)分别画出了A、B、C、D各特征点在T₁时间段隔1帧计算得到的光流的轨迹。

根据步骤(2)和步骤(3)所述，重建卫星主体上4个特征点在空间中的三维运动场。

图10画出了在T₁时间段4个特征点的估计运动轨迹与这些轨迹拟合圆的圆心(“*”所示)，图中M点(“Δ”标注)为各个圆心加权平均后的点，该点在估计的旋转轴L’(图中虚线所示)上，实线为实际的旋转轴L。

各个点的运动轨迹所在平面的法线方向与实际旋转轴的方向的夹角记为Δθ，运动轨迹拟合圆的圆心距离实际旋转轴的距离记为d，运动轨迹拟合圆的半径记为R，下表列出了由4个特征点的运动轨迹计算出的转轴的方向矢量、Δθ、d与R。

表1由隔1帧计算的光流估计第一个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-0.08，-0.32}与轴真实的方向矢量{1，-0.07，-0.37}的误差夹角为2.2°，估计轴与真实轴的空间距离误差为0.1m。

表2由隔2帧计算的光流估计第一个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-0.08，-0.34}与轴真实的方向矢量{1，-0.07，-0.37}的误差夹角为2.2度，估计轴与真实轴的空间距离误差为0.1m。

表3由隔3帧计算的光流估计第一个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-0.08，-0.36}与轴真实的方向矢量{1，-0.07，-0.37}的误差夹角为2.0度，估计轴与真实轴的空间距离误差为0.1m。

表4由隔4帧计算的光流估计第一个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-0.07，-0.33}与轴真实的方向矢量{1，-0.07，-0.37}的误差夹角为1.9度，估计轴与真实轴的空间距离误差为0.1m。

为了说明隔若干帧计算光流能提高旋转轴的估计精度，我们将本实例中B点的光流分别用逐帧取点计算和隔一帧取点计算进行比较。

图11为以相邻帧计算光流后，由B点三维运动矢量估计旋转轴的情况，其中实线L为实际旋转轴，虚线L’为估计旋转轴，P点为估计的圆弧轨迹的圆心，用逐帧取点估计的旋转轴与实际旋转轴的误差夹角为10.3度，估计轴与旋转轴的空间距离误差为0.53m。

图12为按时序依次隔一帧计算光流后，由B点三维运动矢量估计旋转轴的情况，其中实线L为实际旋转轴，虚线L’为估计旋转轴，P点为估计的圆弧轨迹的圆心，用隔一帧帧取点估计的旋转轴与实际旋转轴的误差夹角为2.6度，估计轴与旋转轴的空间距离误差为0.15m。

由此可见，隔帧计算光流，确实能在很大程度上提高旋转轴的估计精度。

根据步骤(2)和步骤(3)所述，重建卫星帆板上4个特征点在空间中的三维运动场。

图14画出了在T₂时间段4个特征点的估计运动轨迹与这些轨迹拟合圆的圆心(“*”所示)，图中的“Δ”为各个圆心加权平均后的点，该点在估计的旋转轴(图中虚线所示)上，实线为真实的旋转轴。

各个点的运动轨迹所在平面的法线方向与实际旋转轴的方向的夹角记为Δθ，运动轨迹拟合圆的圆心距离实际旋转轴的距离记为d，运动轨迹拟合圆的半径记为R，下表列出了由4个特征点的运动轨迹计算出的转轴的方向矢量、Δθ、d与R。

表5由隔1帧计算的光流估计第二个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-3.14，3.14}，与轴真实的方向矢量{1，-3.19，3.26}的误差夹角为2.3度，估计轴与真实轴的空间距离误差为0.4m。

表6由隔2帧计算的光流估计第二个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-3.16，3.18}，与轴真实的方向矢量{1，-3.19，3.26}的误差夹角为2.2度，估计轴与真实轴的空间距离误差为0.3m。

表7由隔3帧计算的光流估计第二个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-3.15，3.21}，与轴真实的方向矢量{1，-3.19，3.26}的误差夹角为2.2度，估计轴与真实轴的空间距离误差为0.3m。

表8由隔4帧计算的光流估计第二个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-3.16，3.22}，与轴真实的方向矢量{1，-3.19，3.26}的误差夹角为2.1度，估计轴与真实轴的空间距离误差为0.2m。

根据步骤(2)和步骤(3)所述，重建卫星帆板上4个特征点在空间中的三维运动场。

图16画出了在T₃时间段4个特征点的估计运动轨迹与这些轨迹拟合圆的圆心(“*”所示)，图中的“Δ”为各个圆心加权平均后的点，该点在估计的旋转轴(图中虚线所示)上，实线为实际的旋转轴。

各个点的运动轨迹所在平面的法线方向与实际旋转轴的方向的夹角记为Δθ，运动轨迹拟合圆的圆心距离实际旋转轴的距离记为d，运动轨迹拟合圆的半径记为R，下表列出了由4个特征点的运动轨迹计算出的转轴的方向矢量、Δθ、d与R。

表9由隔1帧计算的光流估计第三个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-0.09，-1.33}，与轴真实的方向矢量{1，-0.11，-1.19}的误差夹角为4.4度，估计轴与真实轴的空间距离误差为0.3m。

表10由隔2帧计算的光流估计第三个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-0.09，-1.26}，与轴真实的方向矢量{1，-0.11，-1.19}的误差夹角为4.0度，估计轴与真实轴的空间距离误差为0.2m。

表11由隔3帧计算的光流估计第三个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-0.10，-1.26}，与轴真实的方向矢量{1，-0.11，-1.19}的误差夹角为4.0度，估计轴与真实轴的空间距离误差为0.2m。

表12由隔4帧计算的光流估计第三个旋转轴的误差分析

加权平均后，最终的估计轴的方向矢量为{1，-0.10，-1.23}，与轴真实的方向矢量{1，-0.11，-1.19}的误差夹角为3.8度，估计轴与真实轴的空间距离误差为0.2m。

图17表示三条估计轴与实际轴的误差角度随计算光流所隔帧数的变化(L₁、L₂、L₃分别对应T₁、T₂、T₃三个时刻)。图18表示三条估计轴与实际轴的误差距离随计算光流所隔帧数的变化(L₁、L₂、L₃分别对应T₁、T₂、T₃三个时刻)。

最终估计的质心与实际质心的绝对误差为0.1m。

根据以上实验结果，总结出在运用本发明对空间目标进行旋转轴估计时应注意以下几点：

(1)可适当减小相机的视场角，使目标的感兴趣部分在图像中占的比例较大，特征点运动时在图像中形成的光流较大，提高估计精度。

(2)可适当增加计算光流时所隔的帧数，以得到较大的光流，提高估计精度。

(3)选择特征点时，应选取旋转半径较大的特征点，以得到较大的光流，若旋转半径都过小，特征点的运动在图像中形成的光流较小，将会影响对旋转轴方向的估计。

Claims (4)

1.一种基于双目光流的空间目标旋转轴及质心估计方法，包括步骤：

(1)对左、右视图进行特征点检测，左、右视图分别由左、右相机同时拍摄，并且包括至少五个时刻的图像对；

(2)相隔至少一帧按时序依次分别计算左、右视图的光流，得到特征点的双目光流；

(3)根据特征点的双目光流，通过直线重建的方法，进行特征点单位时间内三维运动矢量的重建；

(4)由单个特征点的三维运动矢量序列对旋转轴进行单次估计，并计算出该特征点运动轨迹半径；

(5)将各轨迹平面的法矢量以及轨迹圆心进行加权平均，得到估计的旋转轴的空间方程；

(6)当目标的旋转轴改变后，重复步骤(1)～(5)，得到第二条旋转轴的空间方程，与前一次求出的旋转轴的方程合并，估计出目标的质心空间坐标。

2.根据权利要求1所述的基于双目光流的空间目标旋转轴及质心估计方法，其特征在于，步骤(3)具体包括下述过程：

设旋转之前的P点对应于旋转之后的P′点，

即为立方体上某一点的三维运动矢量；O₁X₁Y₁Z₁为左相机坐标系，O₂X₂Y₂Z₂为右相机坐标系，则空间三维运动矢量所在直线在空间坐标系中的方程为式(I)：

$\left\{\begin{array}{c}x=\mathrm{az}+p\\ y=\mathrm{bz}+q\end{array}\right.,q\≠0---\left(I\right)$

其中，a、b、p、q由下式(II)计算：

$\left\{\begin{array}{c}a=\frac{m_1}{f}+\frac{k_1(m_1-m_2)}{f(k_2-k_1)},\\ b=\frac{m_1-m_2}{f(k_2-k_1)},\\ p=\frac{k_1{b}_0}{k_2-k_1},\\ q=\frac{b_0}{k_2-k_1}.\end{array}\right.---\left(\mathrm{II}\right)$

其中b₀和f分别为两个相机之间的基线长度和焦距，

k₁、m₁分别为直线l₁的斜率和截距，k₂、m₂分别为直线l₂的斜率和截距，l₁为左相机成像平面上的光流

所在的直线，l₂为右相机成像平面上的光流

所在的直线。

3.根据权利要求2所述的基于双目光流的空间目标旋转轴及质心估计方法，其特征在于，步骤(4)包括下述过程：

设任一时段t内，由双目立体光流场重建得到三维运动矢量 $\stackrel{\→}{\mathrm{AB}}={(x,y,z)}^T,$ 以及A点的坐标(x_A，y_A，z_A)^T，将Q₁、Q₂……Q₁₂作为待求解参数建立如下方程组：

$\left[\begin{array}{cccccccccccc}-x_A& 0& 0& y_A& z_A& 0& 0& z_A& -y_A& 1& 0& 0\\ 0& -y_A& 0& x_A& 0& z_A& -z_A& 0& x_A& 0& 1& 0\\ 0& 0& -z_A& 0& x_A& y_A& y_A& -x_A& 0& 0& 0& 1\end{array}\right]\·\left[\begin{array}{c}{Q}_1\\ .\\ .\\ .\\ Q_{12}\end{array}\right]=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]---\left(\mathrm{III}\right)$

观测该特征点在至少四个时间段的运动变化，使方程的系数矩阵的行数大于12，解出Q₁……Q₁₂；

再由式(IV)：

$\left[\begin{array}{ccc}{Q}_1& Q_9-Q_4& -(Q_5+Q_8)\\ -(Q_9+Q_4)& Q_2& Q_7-Q_6\\ Q_8-Q_5& -(Q_6+Q_7)& Q_3\end{array}\right]\·\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]=\left[\begin{array}{c}{Q}_{10}\\ Q_{11}\\ Q_{12}\end{array}\right]---\left(\mathrm{IV}\right)$

解出该特征点的旋转中心O点的坐标(x_O，y_O，z_O)^T；令旋转轴的单位矢量 $\stackrel{\→}{{\mathrm{\ω}}^{\′}}={({{\mathrm{\ω}}_x}^{\′},{{\mathrm{\ω}}_y}^{\′},{{\mathrm{\ω}}_z}^{\′})}^T,$ 该矢量三个分量的比例关系如式(V)所示：

ω_x′∶ω_y′∶ω_z′＝Q₇∶Q₈∶Q₉    (V)

而单位矢量的分量满足方程ω_x′²+ω_y′²+ω_z′²＝1，因此结合式(V)解出(ω_x′，ω_y′，ω_z′)^T，由于旋转轴的方向与的方向相同，而(x_O，y_O，z_O)^T已由式(IV)解出，所以得到单次估计的旋转轴的方程：

$\frac{X-x_0}{{{\mathrm{\ω}}_x}^{\′}}=\frac{Y-y_0}{{{\mathrm{\ω}}_y}^{\′}}=\frac{Z-z_0}{{{\mathrm{\ω}}_z}^{\′}}.---\left(\mathrm{VI}\right)$

另外，根据求得的O点坐标(x_O，y_O，z_O)^T和由双目立体光流场重建得到的A点坐标(x_A，y_A，z_A)^T，计算圆弧轨迹的半径R为：

$R=\sqrt{{(x_A-x_O)}^2+{(y_A-y_O)}^2+{(z_A-z_O)}^2}---\left(\mathrm{VII}\right).$

4.根据权利要求1或2所述的基于双目光流的空间目标旋转轴及质心估计方法，其特征在于，步骤(5)为：

最终估计的旋转轴的方程为：

$\frac{x-\stackrel{\‾}{x}}{\stackrel{\‾}{a}}=\frac{y-\stackrel{\‾}{y}}{\stackrel{\‾}{b}}=\frac{z-\stackrel{\‾}{z}}{\stackrel{\‾}{c}}---\left(\mathrm{VIII}\right)$

其中，最终估计的圆心坐标(x，y，z)和法线方向{a，b，c}的计算公式如式(IX)和(X)所示：

$\stackrel{\‾}{x}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{x}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right),$ $\stackrel{\‾}{y}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{y}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right),$ $\stackrel{\‾}{z}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{z}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right)---\left(\mathrm{IX}\right)$

$\stackrel{\‾}{a}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{a}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right),$ $\stackrel{\‾}{b}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{b}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right),$ $\stackrel{\‾}{c}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{R_i{c}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{R}_j}\right)---\left(X\right).$

Images