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

Abstract

本发明提供了一种基于双目光流的空间目标旋转轴及质心估计方法，包括：依次分别提取双目序列图像中的特征点，计算形成的双目光流；从特征点的双目光流中剔除非有效重建三维运动轨迹的区域，得到有效区域约束的双目光流，重建出空间目标的三维运动轨迹；多次迭代剔除重建三维运动矢量误差较大的光流区域，根据得到每个特征点的三维运动矢量序列估计旋转轴，将多个特征点估计结果进行加权平均，得到估计的旋转轴的空间方程,根据估计的两条旋转轴联立求解出质心目标的空间坐标。本方法流程可用于观测自旋运动或调姿状态中的空间目标，准确有效重建这些空间目标的三维旋转轨迹，估计自旋运动的空间目标的旋转轴和质心。

Description

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

技术领域

本发明属于数字图像处理与空间科学技术交叉领域，具体涉及一种鲁棒的基于双目光流的空间目标旋转轴及质心估计方法。

背景技术

从相机获取的二维图像信息出发计算三维空间中目标物体的结构信息，位置信息和运动信息，并由此重建，识别和跟踪目标是计算机视觉的基本任务之一，空间目标表面特征点的运动轨迹与其在序列图像中对应特征点之间存在着对应关系。基于空间任务的要求，比如航天器空间交会对接，双星编队绕飞、伴飞，多星编队伴飞导航，需要对空间目标进行探测、识别和跟踪等，并对空间目标进行运动分析和姿态估计，迫切需要具有精度高、不依赖外部信息、较好的感知复杂环境的视觉系统和算法，重建空间目标的姿态，运动轨迹。为此，需要对目标的结构、运动轨迹、空间位置和三维姿态等进行高精度的定量测量。基于序列图像的分析是获取这些信息的重要手段，但是现有的重建方法存在很多局限性，由于离散化的误差，双目光流的重建过程中，存在着非有效重建三维运动矢量的光流区域。目前如何提高对空间目标观测，跟踪，重建的准确性和稳定性，我们在这里提出了一种鲁棒的基于双目光流的空间目标旋转轴及质心估计方法，可以有效地提高对于空间目标的轨迹进行跟踪和重建的精度，准确的重建空间目标的运动信息。

运动结构重建是研究从场景的图像序列中提取场景中目标物体的结构、位置和运动信息的理论和方法。这项技术被广泛的应用机器人导航、汽车智能交通、非接触物体测量等领域。在未知环境中，自主导航系统面临的技术困难已不仅仅是规避静止物体，而且需要适应变化的环境，跟踪、观测、规避运动的物体。比如机器人导航、太空卫星跟踪编队等。很多场合需要视觉系统自动地跟踪和观测未知物体。对未知目标进行近距离的特性观测和相对导航，需要处理的核心问题就是估计视觉系统与目标物体之间的相对位置和相对姿态，建立相对导航的运动方程。为此，首先就要估计目标物体的质心或旋转轴。

针对专利《一种基于双目光流的空间目标旋转轴及质心估计方法》存在的问题，提出了一种有效区域约束的双目光流重建三维运动矢量，并根据得到的三维运动矢量估计空间目标旋转轴和质心的方法，专利《一种基于双目光流的空间目标旋转轴及质心估计方法》在重建的过程中可能存在误差，而没有分析误差产生的原因和没有相应的解决的方法，本发明则分析了这些存在的问题，提出了解决方法。

发明内容

本发明的目的是在未知目标物体其它先验信息的情况下，提供一种鲁棒的基于双目光流估计非合作物体旋转轴和质心的方法。该方法通过双目相机捕获得到基于时间的场景图像序列，从场景的图像序列中提取场景中目标物体的特征点，得到目标跟踪点的双目光流，剔除非有效重建三维运动矢量的光流轨迹，得到有效区域约束的光流轨迹段，多次迭代重建目标的实际三维运动矢量，目的在于估计目标物体的旋转轴和质心。空间物体的旋转可以发生在两种情况，一种是自旋稳定的卫星，另一种是卫星在调姿过程中。

假定在两个不同时间段T1和T2，物体绕两个不同的轴旋转。按动力学的规律，这两个旋转轴是通过物体质心的。我们估计两个时间段的旋转轴位置，求得这两个轴的交点，就可得到物体的质心。步骤如下：

（1）分别提取立体图像对中的特征点；

（2）依次计算序列图像中各特征点形成的双目光流；

（3）从各特征点的双目光流中剔除非有效重建三维运动轨迹的区域，得到有效区域约束的双目光流；

（4）根据（3）得到的有效区域约束的双目光流，重建出各特征点的三维运动轨迹；

（5）计算出各特征点的三维运动轨迹平面的法矢量，依次计算各特征点法矢量与各特征点每段运动轨迹方向矢量的余弦距离，剔除各特征点的余弦距离模最大的那段轨迹；

（6）将各特征点的剔除后的轨迹重新计算轨迹平面的法矢量，求解各特征点的两次法矢量的余弦距离，判断各特征点两次法矢量的余弦距离是否大于设置的余弦距离阈值，如果大于则返回步骤（5）根据各特征点剔除后的轨迹重新计算；

（7）估计各特征点的运动轨迹所在平面，以及各特征点的运动轨迹的圆心和半径；

（8）根据各特征点的运动轨迹的半径，将各轨迹平面的法矢量以及轨迹圆心进行加权平均，得到最终的旋转轴的空间方程；

（9）当空间目标的旋转轴改变后，重复步骤（1）～（8），得到第二条旋转轴的空间方程，与前一次求出的旋转轴的方程合并，计算出空间目标的质心坐标。

本发明提出并设计了一种估计空间目标自旋运动旋转轴和质心的方法。针对空间目标的运动轨迹进行准确的重建，不需要在空间目标物体上设置标记，用边缘检测，线检测结合角点检测算法来确定目标物体上的特征点作为光流的跟踪点。目标物体在旋转过程中在左右相机中形成不同的光流，本发明提取出重建三维运动矢量的有效双目光流，根据有效区域约束的光流重建出特征点的三维运动矢量，根据一系列三维运动矢量拟合出特征点的三维运动轨迹，再多次迭代剔除误差较大的光流轨迹，重新重建出特征点的三维运动矢量拟合的三维运动轨迹。物体旋转时特征点在两个传感器成像面上形成的光流是有差异的，本发明根据这些光流差异重建出特征点的三维运动矢量，根据一系列三维运动矢量拟合出特征点的三维运动轨迹，再根据圆弧运动轨迹半径的不同，将多个特征点的三维运动轨迹所在平面的法矢量以及圆心作加权平均，确定本次自旋运动的旋转轴的空间方程，再利用同样方法估计另一时间段的不同的旋转轴。根据动力学的规律，物体旋转时，旋转轴是通过质心的，计算两个时间段的旋转轴交点就可得到物体的质心的空间位置坐标。

附图说明

图1为本发明空间目标旋转轴及质心估计方法流程图；

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

图3为特征点圆弧轨迹示意图；

图4为双目光流重建空间目标三维运动矢量轨迹示意图；

图5为光流轨迹无法重建的区域形成的拐点示意图；

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

图7为T1时间段，检测得到的4个特征点示意图；

图8(a)为目标绕第一个轴旋转时，由隔1帧计算的光流估计的A特征点的运动轨迹（实线）与实际圆弧轨迹（虚线）的比较示意图；

图8(b)为目标绕第一个轴旋转时，由隔1帧计算的光流估计的B特征点的运动轨迹（实线）与实际圆弧轨迹（虚线）的比较示意图；

图8(c)为目标绕第一个轴旋转时，由隔1帧计算的光流估计的C特征点的运动轨迹（实线）与实际圆弧轨迹（虚线）的比较示意图；

图8(d)为目标绕第一个轴旋转时，由隔1帧计算的光流估计的D特征点的运动轨迹（实线）与实际圆弧轨迹（虚线）的比较示意图；

图9为目标绕第一个轴旋转时，最终估计的旋转轴L’（虚线）与实际旋转轴L（实线）比较示意图；

图10为空间目标在旋转过程中形成的理想的三维运动轨迹示意图；

图11(a)为双目相机拍摄到的场景序列图像中的一帧左图像；

图11(b)为双目相机拍摄到的场景序列图像中的一帧右图像；

图12为立体相机在目标物体旋转时得到的序列图像的光流轨迹以及双目光流无法重建的区域示意图；

图13为双目光流存在非有效重建区域情况下，根据双目光流轨迹重建出的目标物体三维运动轨迹的运动矢量示意图；

图14为在未剔除双目光流非有效重建的区域的情况下，双目光流的重建轨迹的X,Y,Z分量的分析示意图；

图15为在剔除双目光流的非有效重建区域后，根据双目光流轨迹重建出的目标物体三维运动轨迹的运动矢量示意图；

图16为剔除无法重建的双目光流段和误差较大的光流段时，迭代计算重建出的光流轨迹示意图；

图17为立体相机在目标物体旋转时得到的序列图像的光流轨迹以及双目光流段中存在多个无法重建的区域示意图；

图18为在未分割出双目光流的多个无法重建区域情况下，根据双目光流轨迹重建出的目标物体三维运动轨迹的运动矢量示意图；

图19为双目光流段中存在多个无法重建区域的情况下，重建轨迹的X,Y,Z分量的分析示意图；

图20为剔除多个非有效重建的双目光流段和误差较大的光流段时，迭代计算重建出的光流轨迹示意图；；

图21为目标旋转轴与质心估计结果示意图；。

具体实施方式

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

为清楚地理解本发明，首先对本发明中有关概念和术语加以说明：机器视觉中的三个坐标系：图像坐标系(U,V)、摄像机坐标系和世界坐标系。机器视觉采集得到的图像数据是直接用图像坐标系表示的，即图像上的像素行和列位置。而被测物体对象的客观位置、大小、相互关系必须通过世界坐标系才能正确描述。这两个坐标系间的转换又要通过摄像机坐标系。

1.图像坐标系

摄像机采集的数字图像在计算机内可以存储为数组，数组中的每一个元素的值即是图像点的亮度。在图像上定义直角坐标系U-V，每一象素的坐标(U,V)分别是该象素在数组中的列数和行数。故(U,V)是以象素为单位的图像坐标系坐标.在UV坐标系中。原点O定义在摄像机光轴和图像平面的相交处，成为图像的主点。

2.摄像机坐标系

摄像机成像几何关系如图2表示，其中O1，O2分别为左右摄像机的光心，X1与Y1轴与成像平面坐标系的U1轴和V1轴平行，Z1轴为摄像机的光轴与图像平面垂直，光轴与图像平面的交点为图像主点O1'，X2与Y2轴与成像平面坐标系的U2轴和V2轴平行，Z2轴为摄像机的光轴与图像平面垂直，光轴与图像平面的交点为图像主点O2'，点O1与X1、Y1、Z1轴组成的直角坐标系称为左摄像机坐标系，O1O1'为左摄像头的焦距，点O2与X2、Y2、Z2轴组成的直角坐标系称为右摄像机坐标系，O2O2'为右摄像头的焦距。

3.世界坐标系

在环境中选择一个参考坐标系来描述摄像机和物体的位置，该坐标系称为世界坐标系.摄像机坐标系和世界坐标系之间的关系可用旋转矩阵与平移向量来描述.本文中选择左摄像机的坐标系作为世界坐标系。

本发明实施例中，利用可见光双目相机组成的双目立体视觉系统，拍摄自转的目标，总体流程图如图1所示：

本发明中目标的特征点（比如角点，边缘点，拐点，暗区域的亮点以及亮区域的暗点）会在双目相机的成像平面上形成空间目标运动轨迹的二维投影，这些三维运动轨迹在像面上的二维投影称为“光流”。同一个点在两个相机上形成的光流称作光流立体像对，本发明的主要方法就是根据自转目标上多个特征点在双目相机上的光流立体像对，精确的估计空间目标的运动轨迹。

（1）分别提取立体图像对中的特征点

空间目标表面结构丰富，包含大量点基元，线基元，如星主体有棱柱、立方体、正八面体等几何形状；太阳能帆板一般是长方体结构。许多目标由直线或曲线构成，尤其是直线，是组成目标最基础的几何单元之一；许多人造目标是由平行线和垂直线构成的方体、柱体等基本几何形状组合而成的。可以用边缘算子提取场景图像中目标的边缘信息，再根据Hough变换检测直线段，并计算相邻直线段的交点，检测出目标物体的特征点，也可利用Harris角点检测或者SUSAN角点检测的算法辅助检测目标角点，以提高特征点检测的精度。

（2）依次计算序列图像中特征点形成的双目光流

不失一般性，将目标物体做匀速转动，得到基于时间的场景序列图像，将步骤（1）中得到的特征点作为光流跟踪点，依次计算出场景序列图像对中的后续每一帧图像的特征点，得到序列图像中特征点形成的双目光流，控制目标物体匀速转动的速度，可以得到较大的光流或较小的光流。在实际操作时，可以根据具体的情况调整目标物体的自转速度。

（3）根据特征点的双目光流剔除不能重建出三维运动轨迹的区域，得到有效区域约束的双目光流

某些时刻，空间目标运动形成的三维运动矢量在两相机成像面投影得到的光流，无法准确有效重建出三维运动矢量，或者信息的损失较大，导致重建误差也较大，局部误差严重影响后续对于目标运动的分析与估计，这段区域为非有效重建三维运动矢量的双目光流区域，因此在重建三维运动矢量的空间轨迹时，就需要剔除非有效重建三维运动矢量的光流段，得到的就是有效区域约束的双目光流。

如图2，P点在左相机成像平面的投影为P₁，P′点在左相机成像平面的投影为P₁′，那么三维运动矢量

在左相机成像平面的投影形成的光流即为

同理，

在右相机成像平面的投影形成的光流即为

在左相机成像平面坐标系O₁′U₁V₁中，记P₁的坐标为(U₁,V₁)，P₁′的坐标为(U′₁,V₁′)；在右相机成像平面坐标系O′₂U₂V₂中，记P₂的坐标为(U₂,V₂)，P₂′的坐标为(U′₂,V₂′)。

光流所在直线的斜率和截距也可以计算得到，设l₁的斜率为k₁，截距为m₁，l₂的斜率为k₂，截距为m₂（这里定义斜率为

），所以l₁，l₂也可以表示为：l₁:U₁-k₁V₁＝m₁，l₂:U₂-k₂V₂＝m₂。

当光流在U方向上没有视差时，Δu＝0，得k₁＝k₂＝0，当光流在V方向上没有视差时，Δv＝0，得k₁＝k₂＝∞，在这两种情况下，无法通过光流段重建出空间目标特征点的运动矢量,如图5，在没有剔除非有效重建三维运动轨迹的双目光流区域时重建出来的轨迹，存在明显的拐点，图中A标记的区域即为非有效重建三维运动轨迹的双目光流区域。

对于目标物体特征点P的n段光流，在左图像中的光流为P₁P₁'，P'₁P₁''，P₁''P₁'''……P₁ ^(n-1)P₁ ⁽ⁿ⁾，在右图像中的光流为P₂P₂'，P'₂P₂''，P₂''P₂'''……P₂ ^(n-1)P₂ ⁽ⁿ⁾,分别计算每一段光流P₁ ^(i-1)P₁ ⁽ⁱ⁾，P₂ ^(i-1)P₂ ⁽ⁱ⁾，(i＝1,2...n)其在U方向和V方向的差值，当

或P₂ ^(i-1)P₂ ⁽ⁱ⁾在U或V方向上的差值为0时，剔除该段光流。

（4）根据（3）得到的有效双目光流，重建出空间目标的三维运动轨迹

三维运动矢量到双目成像系统的投影模型如图2所示，在该模型中，目标物体绕旋转轴做旋转运动，实线所示的目标物体为旋转之前的状态，虚线所示的目标物体是旋转之后的状态，旋转之前的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₁′，那么三维运动矢量

在左相机成像平面的投影形成的光流即为P₁P₁'，同理，

在右相机成像平面的投影形成的光流即为P₂P₂'。

在左相机成像平面坐标系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′)。

由相机针孔模型可得：

$\begin{array}{cc}\left\{\begin{array}{c}{U}_1=X\frac{f_{x1}}{Z}\\ V_1=Y\frac{f_{y1}}{Z}\end{array}\right.& \left\{\begin{array}{c}{U}_1^{\′}=X^{\′}\frac{f_{x1}}{Z^{\′}}\\ V_1^{\′}=Y^{\′}\frac{f_{y1}}{Z^{\′}}\end{array}\right.\end{array}---\left(1\right)$

$\begin{array}{cc}\left\{\begin{array}{c}{U}_2=(X-b_0)\frac{f_{x2}}{Z}\\ V_2=Y\frac{f_{y2}}{Z}\end{array}\right.& \end{array}\left\{\begin{array}{c}{U}_2^{\′}=(X^{\′}-b_0)\frac{f_{x2}}{Z^{\′}}\\ V_2^{\′}=Y^{\′}\frac{f_{y2}}{Z^{\′}}\end{array}\right.---\left(2\right)$

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

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

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

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

空间三维运动矢量所在直线在空间坐标系中的方程可以写为：

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

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

所在的直线为l₁，右相机成像平面上的光流P₂P'₂所在直线为l₂，根据式(3)、(4)及(5)，l₁，l₂分别在各自的成像平面坐标系的方程可以表示为：

l₁:qf_y1U₁-pf_x1V₁＝(aq-bp)f_x1f_y1    (6)

l₂:qf_y2U₂-(p-b₀)f_x2V₂＝[aq-b(p-b₀)]f_x2f_y2    (7)

其中，b₀,f_x1,f_y1,f_x2,f_y2为两个相机之间的基线长度和左右相机X,Y方向的等效焦距。

由于光流可以直接在图像序列中检测得到，因此光流所在直线的斜率和截距也可以计算得到，设l₁的斜率为k₁，截距为m₁，l₂的斜率为k₂，截距为m₂（这里定义斜率为

），所以l₁，l₂也可以表示为：

l₁:U₁-k₁V₁＝m₁    (8)

l₂:U₂-k₂V₂＝m₂    (9)

根据式(6)～式(9)可得：

$\left\{\begin{array}{c}{k}_1=\frac{p{f}_{x1}}{q{f}_{y1}}\\ m_1=\frac{\mathrm{aq}-\mathrm{bp}}{q}{f}_{x1}\\ k_2=\frac{p-b_0}{q{f}_{y2}}{f}_{x2}\\ m_2=\frac{\mathrm{aq}-\mathrm{bp}+b{b}_0}{q}{f}_{x2}\end{array}\right.---\left(10\right)$

根据式(10)，解出a,b,p,q：

$\left\{\begin{array}{c}a=(m_2{f}_{x1}-m_1{f}_{x2})/(k_1{f}_{x2}{f}_{y1}-k_2{f}_{x1}{f}_{y2})\\ b=(m_2{f}_{x1}-m_1{f}_{x2})/(k_1{f}_{x2}{f}_{y1}-k_2{f}_{x1}{f}_{y2})\\ p=k_1{b}_0{f}_{x2}{f}_{y1}/(k_1{f}_{x2}{f}_{y1}-k_2{f}_{x1}{f}_{y2})\\ q=b_0{f}_{x1}{f}_{x2}/(k_1{f}_{x2}{f}_{y1}-k_2{f}_{x1}{f}_{y2})\end{array}\right.---\left(11\right)$

由式(11)，再根据式(5)就可以得到三维运动矢量在空间坐标系中所在直线的方程。

（5）计算出轨迹平面的法矢量，依次计算法矢量与每段轨迹方向矢量的内积，剔除内积最大的那段轨迹；将剔除后的轨迹重新计算轨迹平面的法矢量，求解两次法矢量的余弦距离

在左右图像中得到某特征点的n段光流，可以重建出空间目标上该特征点绕旋转轴旋转的n段运动矢量，记该特征点为P，其在左图像中的投影为P₁，在右图像中的投影为P₂，左图像中形成的n段光流从P₁到P₁ ⁽ⁿ⁾，右图像中形成的n段光流从P₂到P₂ ⁽ⁿ⁾，重建出的运动矢量为

如图4，计算出当前光流轨迹平面的法矢量，记为

第（4）步计算出的每一段光流的运动矢量，空间三维运动矢量

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

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

则其每一段的运动矢量的方向为

余弦距离可用来衡量两个向量方向的差异，在这里我们选择余弦距离来估算运动矢量的相似性度量，计算每一段运动矢量与光流轨迹平面的法矢量的余弦距离cos(θ_i),(i＝1,2...n)，其方程为：

$\cos \left({\mathrm{\θ}}_i\right)=\frac{{\stackrel{\→}{w_t}}^T\·\stackrel{\→}{w_i}}{\left|\stackrel{\→}{w_t}\right|\left|\stackrel{\→}{w_i}\right|},(i=\mathrm{1,2}...n)---\left(13\right)$

余弦距离取值范围为[-1,1]。夹角余弦越大表示两个向量的夹角越小，夹角余弦越小表示两向量的夹角越大。当两个向量的方向重合时夹角余弦取最大值1，当两个向量的方向完全相反，夹角余弦取最小值-1。计算出最大的余弦距离的模E(j)

E(j)＝max(|cos(θ_i)|),(i＝1,2...n)    (14)

（6）重复迭代第（4）和（5）步中的计算，直到两次的轨迹平面的法矢量的余弦距离不变或者变化很小

在得到的k段双目光流中，剔除余弦距离的模最大的一段双目光流，得到剩余的k-1段双目光流，重新计算k-1段双目光流估计的法矢量计算两次法矢量的余弦距离

$\cos \left(\mathrm{\α}\right)=\frac{{\stackrel{\→}{w_t}}^T\·\stackrel{\→}{w_t^{\′}}}{\left|\stackrel{\→}{w_t}\right|\left|\stackrel{\→}{w_t^{\′}}\right|}---\left(15\right)$

设置余弦距离的阈值为L,当两次法矢量的余弦距离大于阈值时，重复第（4）步和第（5）步，迭代计算，直到余弦距离小于或等于阈值。

某些时刻，当三维运动矢量在两相机成像面投影得到的光流大小远小于三维运动矢量本身的大小时（即三维运动矢量与成像平面的角度较大时），由于光流的双目视差较小，所以在重建三维运动矢量时，深度信息的损失较大，所以重建误差也较大，因此在重建三维运动矢量的空间轨迹时，就会表现出在某些时刻重建的轨迹中会出现明显的拐点（如图6中A、B两处所示）。基于以上分析，在对旋转轴进行估计之前，应首先估计重建的三维运动矢量序列的双目光流的有效区域，若不对双目光流进行处理，会严重影响旋转轴的估计精度。图6中，L为实际旋转轴，L’为估计的旋转轴，P点为估计错误的圆弧轨迹的圆心。

（7）估计各特征点的运动轨迹所在平面，以及轨迹的圆心。

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

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

与的矢量夹角为

那么

与

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

因目标旋转的角速度为

且

那么由式(16)可得：

$\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(17\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(18\right)$

其中，为

的单位向量，那么 $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(19\right)$

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

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

$M=\left[\begin{array}{ccc}-(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(21\right)$

令A点的坐标为(x_A,y_A,z_A)^T，O点的坐标为(x_O,y_O,z_O)^T，向量

可由双目立体光流场重建得到，设

那么：

$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(22\right)$

将参数作如下变量代换：

$Q_1=(1-\cos \mathrm{\Ωt})(1-{\mathrm{\ω}}_x^{\′2}),$

$Q_2=(1-\cos \mathrm{\Ωt})(1-{\mathrm{\ω}}_y^{\′2}),$

$Q_3=(1-\cos \mathrm{\Ωt})(1-{\mathrm{\ω}}_z^{\′2}),$

Q₄＝(1-cosΩt)ω'_xω'_y,

Q₅＝(1-cosΩt)ω'_xω'_z,

Q₆＝(1-cosΩt)ω'_yω'_z,

Q₇＝ω'_xsinΩt,

Q₈＝ω'_ysinΩt,

Q₉＝ω'_zsinΩ,

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₅-y₀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]\·\left[\begin{array}{c}{Q}_1\\ \·\\ \·\\ \·\\ Q_{12}\end{array}\right]=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]---\left(23\right)$

如果只观测一个特征点在一个时间段的空间位置变化，可以得到上述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]=\left[\begin{array}{c}{Q}_{10}\\ Q_{11}\\ Q_{12}\end{array}\right]---\left(24\right)$

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

而ω'_x:ω'_y:ω'_z＝Q₇:Q₈:Q₉    (25)

又由于

是个单位矢量，所以有结合式（25）可以解出(ω'_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(26\right)$

半径为 $R^{\′}=\frac{|(x_O-x_A,y_O-y_A,z_O-z_A)\×(x,y,z)|}{\sqrt{x^2+y^2+z^2}}$

另外，还可以根据 $\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^{\′}}$ 求出物体旋转的角速度的大小。

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

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

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

令每个圆弧轨迹的圆心坐标为(x_i,y_i,z_i)，半径为r_i，轨迹拟合平面的法线方向为(a_i,b_i,c_i)，那么最终估计的圆心坐标

和法线方向

的计算公式如下所示：

$\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(27\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(28\right)$

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

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

（9）在T₂、T₃时间段，当目标的旋转轴改变后，重复步骤（1）～（8），得到另一条旋转轴的空间方程，与前一次求出的旋转轴的方程联立，求解出目标的质心空间坐标。

本发明通过序列图像处理和双目光流计算，得到有效区域约束的双目光流区域，多次迭代重建出空间目标的三维运动轨迹，来估计空间目标的旋转轴及质心，用于高精度地定量分析空间目标的运动轨迹和姿态估计。应用于观测卫星对在太空中处于自旋运动或姿态调节过程中的空间物体的特性测量及观测卫星与空间物体在逐渐靠近过程中的自主图像导航。

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

在T₁时段，仿真卫星与相机相距20m，相机的视场角为10度，卫星的旋转轴为图7中的L，检测得到的特征点为图7中的A、B、C、D。分别由特征点A,B,C,D隔1帧计算的光流估计的旋转轴的结果和真实运动场的对比如图8(a)-(d)（图中直线代表理想转轴，星号标注的是根据步骤(3)得到的圆弧轨迹的圆心估计结果），实线代表特征点的运动轨迹，虚线代表实际圆弧轨迹，L为实际旋转轴，P为估计运动轨迹的圆心。

图9为以相邻帧计算光流后，由B点三维运动矢量估计旋转轴的情况，其中实线L为实际旋转轴，虚线L’为估计旋转轴，P点为估计的圆弧轨迹的圆心。

在实际的实验过程中，将神州飞船的模型作为我们的观测目标，将神州飞船的太阳能帆板上的特征点作为跟踪点，图11给出了立体相机获取的场景序列图像中的一帧图像，通过提取场景序列图像中的特征点在左右图像对中的坐标，形成双目光流，图12表示目标物体在运动过程中在双目相机上形成的投影光流，并标记出了存在无法重建区域的光流段，在未分割出双目光流的无法重建区域情况下，根据双目光流轨迹重建出的目标物体三维运动轨迹的运动矢量如图13所示，在无法重建区域段重建出的轨迹存在明显的拐点，在图14中，分析了双目光流的重建轨迹的每一段运动矢量相交的交点的X,Y,Z分量的变化，在图像光流U方向或者V方向上差异较小时，会形成较大的误差，如图12中椭圆标识的区域，图像光流在V方向没有差异，X,Y,Z分量均为0，在剔除双目光流的非有效重建三维运动矢量的区域后，根据有效区域约束的双目光流重建出的目标物体三维运动矢量，如图15所示，将重建出来的轨迹进行多次迭代计算，剔除误差较大的光流段，如图16所示。

考虑到存在多个非有效重建三维运动矢量的双目光流区域的情况，分析了在U方向和V方向上（如图17）均存在死区的情况，在两个方向上都存在没有差值的光流段，在没有剔除非有效重建三维运动矢量的双目光流区域的情况下，可以看出重建的三维运动矢量在非有效重建区域内存在明显的偏差（如图18），在图19中，给出了双目光流段中存在多个非有效重建三维运动矢量的双目光流区域的情况下，重建轨迹的X,Y,Z分量的重建结果，图20给出了剔除多个无法重建的光流段和误差较大的光流段时，迭代计算重建出的光流轨迹。图21分别是目标旋转轴在两个方向上旋转重建出的光流轨迹，实际的光流轨迹以及估计出的旋转轴和实际的旋转轴，以及通过估计出的两方向的旋转轴求解出的空间目标的质心以及目标的实际质心。

Claims (7)

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

（1）分别提取立体图像对中的特征点；

（2）依次计算序列图像中各特征点形成的双目光流；

（3）从各特征点的双目光流中剔除非有效重建三维运动轨迹的区域，得到有效区域约束的双目光流；

（4）根据（3）得到的各特征点有效区域约束的双目光流，重建出各特征点的三维运动轨迹；

（5）计算出各特征点的三维运动轨迹平面的法矢量，依次计算各特征点法矢量与各特征点每段运动轨迹方向矢量的余弦距离，剔除各特征点的余弦距离模最大的那段轨迹；

（6）将各特征点的剔除后的轨迹重新计算轨迹平面的法矢量，计算各特征点的两次法矢量的余弦距离，判断各特征点两次法矢量的余弦距离是否大于设置的余弦距离阈值，如果大于则返回步骤（5）根据各特征点剔除后的轨迹重新计算；

（7）估计各特征点的运动轨迹所在平面，以及各特征点的运动轨迹的圆心和半径；

（8）根据各特征点的运动轨迹的半径，将各轨迹平面的法矢量以及轨迹圆心进行加权平均，得到最终的旋转轴的空间方程；

（9）当空间目标的旋转轴改变后，重复步骤（1）～（8），得到第二条旋转轴的空间方程，与前一次求出的旋转轴的方程合并，计算出空间目标的质心坐标。

2.根据权利要求1所述的鲁棒的基于双目光流的空间目标旋转轴及质心估计方法，其特征在于，所述步骤（3）中从各特征点的双目光流中剔除非有效重建三维运动轨迹的区域具体为：

对于目标物体特征点P的n段光流，在左图像中的光流为P₁P₁'，P'₁P₁''，P₁''P₁'''……P₁ ^(n-1)P₁ ⁽ⁿ⁾，在右图像中的光流为P₂P₂'，P'₂P₂''，P₂''P₂'''……P₂ ^(n-1)P₂ ⁽ⁿ⁾，分别计算每一段光流P₁ ^(i-1)P₁ ⁽ⁱ⁾，P₂ ^(i-1)P₂ ⁽ⁱ⁾，(i＝1,2...n)在U方向和V方向的差值，当P₁ ^(i-1)P₁ ⁽ⁱ⁾或P₂ ^(i-1)P₂ ⁽ⁱ⁾在U或V方向上的差值为0时，剔除该段光流。

3.根据权利要求1或2所述的鲁棒的基于双目光流的空间目标旋转轴及质心估计方法，其特征在于，所述步骤（4）中根据步骤（3）得到的各特征点有效区域约束的双目光流，重建出各特征点的三维运动轨迹具体为：

在左右图像中得到特征点P的n段光流，重建出空间目标上该特征点绕旋转轴旋转的n段运动矢量，其在左图像中的投影为P₁，在右图像中的投影为P₂，左图像中形成的n段光流从P₁到P₁ ⁽ⁿ⁾，右图像中形成的n段光流从P₂到P₂ ⁽ⁿ⁾，重建出的运动矢量为

从而得到该特征点的三维运动轨迹，

其中运动矢量i＝1,2...n，在空间坐标系中的方程为：

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

$\left\{\begin{array}{c}a=(m_2{f}_{x1}-m_1{f}_{x2})/(k_1{f}_{x2}{f}_{y1}-k_2{f}_{x1}{f}_{y2})\\ b=(m_2{f}_{x1}-m_1{f}_{x2})/(k_1{f}_{x2}{f}_{y1}-k_2{f}_{x1}{f}_{y2})\\ p=k_1{b}_0{f}_{x2}{f}_{y1}/(k_1{f}_{x2}{f}_{y1}-k_2{f}_{x1}{f}_{y2})\\ q=b_0{f}_{x1}{f}_{x2}/(k_1{f}_{x2}{f}_{y1}-k_2{f}_{x1}{f}_{y2})\end{array}\right.,$

其中，b₀,f_x1,f_y1,f_x2,f_y2为两个相机之间的基线长度和左右相机分别在X,Y方向的等效焦距，k₁和m₁分别为该特征点P在左相机成像平面上的光流

所在的直线l1的斜率和截距，k₂和m₂分别为该特征点P在右相机成像平面上的光流

所在的直线l₂的斜率和截距。

4.根据权利要求1至3任一项所述的鲁棒的基于双目光流的空间目标旋转轴及质心估计方法，其特征在于，步骤（5）具体包括：

根据步骤（4）中得到的各特征点的每一段空间目标的三维运动矢量

计算出各特征点当前光流轨迹平面的法矢量，记为

记每一段的运动矢量为 $\stackrel{\→}{w_i}=(a,b,1),(i=\mathrm{1,2}...n);$

计算各特征点每一段运动矢量与光流轨迹平面的法矢量的余弦距离cos(θ_i),(i＝1,2...n)，其方程为：

剔除各特征点余弦距离的模最大的一段双目光流，各特征得到剩余的k-1段双目光流。

5.根据权利要求4所述的鲁棒的基于双目光流的空间目标旋转轴及质心估计方法，其特征在于，所述步骤（6）具体包括：

重新计算各特征点k-1段双目光流估计的法矢量

计算各特征点两次法矢量的余弦距离

当两次法矢量的余弦距离大于余弦距离的阈值L时，则返回步骤（5）根据各特征点剔除后的轨迹重新计算。

6.根据权利要求1至5任一项所述的鲁棒的基于双目光流的空间目标旋转轴及质心估计方法，其特征在于，所述步骤（7）中：

各特征点的运动轨迹的圆心坐标为(x_O,y_O,z_O)^T，其利用下式计算：

$\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],$

各特征点的运动轨迹的法矢量方向为其计算公式如下：

ω'_x:ω'_y:ω'_z＝Q₇:Q₈:Q₉；

各特征点的运动轨迹的半径为 $R^{\′}=\frac{|(x_O-x_A,y_O-y_A,z_O-z_A)\×(x,y,z)|}{\sqrt{x^2+y^2+z^2}}$

其中：

$Q_1=(1-\cos \mathrm{\Ωt})(1-{\mathrm{\ω}}_x^{\′2}),$

$Q_2=(1-\cos \mathrm{\Ωt})(1-{\mathrm{\ω}}_y^{\′2}),$

$Q_3=(1-\cos \mathrm{\Ωt})(1-{\mathrm{\ω}}_z^{\′2}),$

Q₄＝(1-cosΩt)ω'_xω'_y,

Q₅＝(1-cosΩt)ω'_xω'_z,

Q₆＝(1-cosΩt)ω'_yω'_z,

Q₇＝ω'_xsinΩt,，

Q₈＝ω'_ysinΩt,

Q₉＝ω'_zsinΩ,

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₅-y₀Q₆-y₀Q₇+x₀Q₈.

(x_A,y_A,z_A)^T为该特征点运动轨迹的起点坐标，(x,y,z)^T为该特征点运动轨迹所对应的光流，Ω为目标旋转的角速度矢量

的模。

7.根据权利要求1至6任一项所述的鲁棒的基于双目光流的空间目标旋转轴及质心估计方法，其特征在于，所述步骤（8）具体包括：

最终的圆心坐标

和法线方向

的计算公式如下所示：

$\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),$

$\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),$

最终的旋转轴的方程为：

$\frac{x-\stackrel{\‾}{x}}{\stackrel{\‾}{a}}=\frac{y-\stackrel{\‾}{y}}{\stackrel{\‾}{b}}=\frac{z-\stackrel{\‾}{z}}{\stackrel{\‾}{c}},$

其中第i个特征点运动轨迹的圆心坐标为(x_i,y_i,z_i)，半径为r_i，轨迹平面的法线方向为(a_i,b_i,c_i)，i＝1,2...n。

Images