CN104616292B - 基于全局单应矩阵的单目视觉测量方法 - Google Patents

Abstract

一种基于全局单应矩阵的单目视觉测量方法。本发明提出了一种图像平面与所观测地面之间的全局单应矩阵标定方法，从而获得图像平面与整个场景平面之间的一一映射关系。首先，将标定板放置于地面不同位置处，获得多个标定板与对应图像平面直接的局部单应矩阵，然后将多个局部单应矩阵进行信息融合，从而得到全局意义上的映射关系，即全局单应矩阵。同时，本发明对关联高度信息的单应矩阵进行了标定，从而能够对任意已知高度的待测平面进行视觉测量。本发明无需摄像机内参数，且标定精度较高。标定结果成功地应用于室内移动机器人位姿测量。对比实验结果表明，在整个摄像机视野范围内，相比局部单应矩阵的标定方法，本发明具有更高的视觉测量精度。

Description

基于全局单应矩阵的单目视觉测量方法

技术领域

本发明属于单目视觉测量与单应矩阵标定的技术领域，特别是涉及一种基于全局单应矩阵的单目视觉测量方法。

背景技术

作为典型的非接触式传感器，视觉传感器具有信息丰富、可靠性高等优点，从而得到了广泛的应用。为了测量物体的位姿，单目视觉和立体视觉通常被认为是比较有效的测量工具[1]。通常来说，单目视觉测量方法是利用目标点间固有的几何约束关系进行求解[2]，而双目立体视觉测量方法，是通过对同一目标点在两台摄像机中的成像进行交会三角计算，从而得到目标点在测量坐标系下的三维坐标[3,4]。双目视觉测量设备受到视场范围的限制[5]，只能测量两摄像机视野重叠的部分，测量范围小并且短。而单目视觉测量系统不受视场范围的限制，可以通过变换焦距，在所需测量范围以及测量距离内测量靶标的空间位置姿态[6]。此外，单目视觉测量具有结构简单、标定步骤少等优点。

单目视觉测量方法可以利用多种特征点，例如点特征、线特征以及其他的高级几何特征[7]。相对而言，目前对基于点特征的单目视觉测量方法研究较多。点特征定位又称为PnP(Perspective n-Points Problem)问题，它是在1981年首先由Fischler和Bolles[8]提出的。给定n(n≥3)个三维参考点和它们相应的二维图像坐标，PnP问题的目标就是获得摄像机的位置和姿态[9][10]。在过去的十几年里，针对不同数目的特征点，很多学者对PnP问题进行了大量的研究工作(从至少3个点，到一般化的n个点)。文献[11]、文献[12]等对P3P问题进行了系统的研究，文献[13]对P4P和P5P问题也进行了研究，文献[14]、[15]、[16]等对于一般化的PnP(n≥4)问题提出了一些线性方法。在文献[17]、[18]中，消影点和平行直线被用于解决视觉测量问题，但是这些方法会受到场景的限制，因为很多场景中不包含足够多的平行直线。利用场景约束，文献[19,20]利用一幅未标定的图像进行三维场景的重建。

在一些应用中，静止或者运动的被测目标点均位于特定高度的同一水平面上。在这种情况下，可以利用该平面约束，进行单应矩阵标定来解决视觉测量问题。文献[9]指出，单应矩阵表示两平面间的可逆齐次变换[21]。一旦该单应矩阵被标定，我们可以直接将一个平面上图像坐标映射到目标平面上的坐标。

为了对图像平面与视野平面之间的单应矩阵进行标定，我们通常将一个标定板放置在摄像机视野范围内。文献[22]利用一个特定的模板来标定摄像机外参数，但是标定精度有待于提高。文献[23]虽然获得了较高的标定精度，但是它需要复杂的实验设备来保证标定性能。根据标定板与其在图像平面上的点的对应关系，通过线性估计和非线性优化[24]，可以标定得到一个单应矩阵。为了保证标定精度，通常需要标定板足够大，最好能够覆盖整个摄像机视野范围。然而，在很多应用中，摄像机的视野范围远远大于标定板的大小。在这种情况下，由于缺乏整个视野范围内的标定数据，得到的单应矩阵仅仅在标定板放置的局部区域内标定精度高。而且，由于摄像机镜头的畸变，不同区域内，图像畸变情况不一样。通过以上分析，我们知道，如果我们只将标定板放在一个位置处标定得到一个单应矩阵，其在整个视野范围内的平均测量误差会比较大。

发明内容

本发明的目的是解决现有局部单应矩阵存在的上述不足，提供一种基于全局单应矩阵的单目视觉测量方法。

为了减小整个视野范围内的测量误差，本发明提出了将多个局部单应矩阵进行信息融合与优化得到全局单应矩阵的方法。首先，我们将标定板放置在视野范围内的不同位置，并标定得到不同的局部单应矩阵。然后，通过分析各局部单应矩阵之间的关系，并将其进行融合、优化得到全局单应矩阵。在整个视野范围内，本发明提出的全局单应矩阵比单一的局部单应矩阵具有更高的测量精度。另外，为了得到具有一定高度平面上点的测量结果，我们定义并标定了一个关联高度信息的单应矩阵，它可以将任意已知高度平面上的点的图像坐标转换到相应零平面上的点的图像坐标。在此基础上,再根据全局单应矩阵得到被测点在待测平面上的二维世界坐标(第三维世界坐标为平面的高度)。大量实验结果验证了本发明的有效性。

本发明提供的基于全局单应矩阵的单目视觉测量方法包括：

第1，全局单应矩阵的标定

定义“局部单应矩阵”为标定板在待测地面单个位置处，标定得到的与图像平面之间的单应矩阵。这里，采用“局部”，是因为标定板的尺寸比较小，只能覆盖整个视野范围内待测地面的小部分区域。相应地，我们定义“全局单应矩阵”，用来表示将标定板放置在不同位置得到多幅图像，并将这些位置的数据进行融合得到的单应矩阵。“局部单应矩阵”反映了图像平面和局部区域之间的映射关系，而“全局单应矩阵”，利用不同位置的数据，更准确地反映了图像平面和整个场景平面之间的映射关系。

图1是一个单目视觉测量系统，将一个未标定的摄像机固定在某高度的金属杆顶端，将标定板放置在摄像机视野范围内待测地面(零平面)的不同位置处。以位置i处的标定板左下角为原点建立坐标系不失一般性，选择标定板在第一个位置处所建立的坐标系为参考世界坐标系。

第1.1，建立各局部单应矩阵之间的关系

用局部单应矩阵H_i建立特征点图像坐标与世界坐标之间的关系如下：

p_ik＝λ_ikH_iP_ik (3)

其中，λ_ik表示归一化比例因子，P_ik＝[x_ik y_ik 1]^T表示标定板上特征点在X_wY_w平面上的二维齐次世界坐标，p_ik＝[u_ik v_ik 1]^T表示相应的齐次图像坐标，i＝1,2,…,N；k＝1,2,…,N_p，N表示标定板的放置次数，N_p表示标定板上特征点的个数。

根据空间几何关系，推导出在参考坐标系处获得的局部单应矩阵H₁与在其他位置处获得的局部单应矩阵H_i之间的关系为：

其中，为坐标系在坐标系下的变换矩阵，¹θ_i表示两坐标系之间的旋转角度，[¹t_xi ¹t_yi]^T表示两坐标系之间的平移向量，表示归一化比例因子。

这样，根据任意一个局部单应矩阵和参考坐标系处的单应矩阵，可以求出它们之间的变换矩阵¹M_i，进而求出旋转与平移参数¹θ_i、¹t_xi和¹t_yi。

然后我们通过变换矩阵¹M_i，将坐标系下的坐标转换到坐标系下，从而可以建立多组约束方程：

由于i＝1,2,…,N；k＝1,2,…,N_p，根据式(8)可以得到N·N_p组约束方程，构成约束方程组。

第1.2，利用非线性最小二乘算法求解全局单应矩阵

采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少1个局部单应矩阵最小化如下目标函数J，得到全局单应矩阵的解：

其中

其中，为归一化比例因子，为参数¹θ_i,¹t_xi,¹t_yi的估计量。非线性优化的初始值选择为：优化后，得到融合N个局部单应矩阵的全局单应矩阵

第2，关联高度信息的单应矩阵的估计

由第1步，得到了全局单应矩阵利用该全局单应矩阵，能够根据图像计算得到零平面(高度为0的待测地面)上特征点在参考世界坐标系中的世界坐标。然而，在很多实际应用中，被测量的特征点并不在零平面上，而是在某一特定高度的平面上。为进一步解决这个问题，我们估计出了关联高度信息的单应矩阵，该单应矩阵描述了位于零平面(平面0)上的点的图像坐标与位于高度为h的平面(平面h)上的点的图像坐标之间的关系，如图2所示。一旦该关联高度信息的单应矩阵被标定，给定高度h，我们就可以通过该单应矩阵将位于平面h上的点的图像坐标x_hi转换到平面0上的点的图像坐标x_0i，然后利用第1步得到的全局单应矩阵即可求得位于零平面上点的二维世界坐标。

第2.1，建立关联高度信息的单应矩阵与高度的关系

考虑静态场景中的N_h个特征点，X_0i和X_hi分别表示平面0和平面h上对应点的齐次世界坐标：

X_0i＝[x_i y_i 0 1]^T,X_hi＝[x_i y_i h 1]^T (13)

相应的图像齐次坐标x_0i和x_hi为：

x_0i＝[u_0i v_0i 1]^T,x_hi＝[u_hi v_hi 1]^T (14)

推导出平面0和平面h上的点的图像坐标之间的关系为：

其中

c_ij(i＝1,2,3；j＝1,2,4)为矩阵的元素，它们通过第1步已经优化得到，而c₁₃、c₂₃和c₃₃为待求参数。若高度h已知，给定至少两对点，可以得到至少4个约束方程，即可对参数c₁₃、c₂₃和c₃₃进行优化求解。

第2.2，利用非线性最小二乘算法对关联高度信息的单应矩阵进行估计

采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少2组空间点对最小化如下目标函数J_h，得到c₁₃,c₂₃,c₃₃的解：

其中

λ_i为归一化比例因子，为对变量c₁₃,c₂₃,c₃₃进行非线性优化后得到的最终结果。

最后利用c₁₃,c₂₃,c₃₃得到然后再结合全局单应矩阵将h设为变量，即得到关联高度信息的单应矩阵因此，对任意已知高度h，便可代入得到G(h)，将位于平面h上的点的图像坐标x_hi转换到平面0上的点的图像坐标x_0i，然后利用第1步得到的全局单应矩阵即可求得位于零平面上点的二维世界坐标。

本发明方法的理论依据及推导过程

第1，全局单应矩阵的标定

图1是一个单目视觉测量系统，一个未标定的摄像机固定在某高度的金属杆顶端，将标定板放置在摄像机视野范围内待测地面(零平面)的不同位置处。以位置i处的标定板左下角为原点建立坐标系不失一般性，选择标定板在第一个位置处所建立的坐标系为参考世界坐标系。

第1.1，建立各局部单应矩阵的之间的关系

假设一块标定板上有N_p个特征点，将该标定板放置在N个不同的位置。定义P_ik(i＝1,2,…,N；k＝1,2,…,N_p)为特征点在坐标系下的三维齐次世界坐标：

P_ik＝[x_ik y_ik 1]^T (1)

相应的齐次图像坐标为：

p_ik＝[u_ik v_ik 1]^T (2)

在位置i处，图像坐标和世界坐标之间可以通过一个局部单应矩阵建立联系：

p_ik＝λ_ikH_iP_ik (3)

其中，H_i表示位置i处图像坐标和世界坐标之间的3×3维局部单应矩阵，λ_ik表示归一化比例因子。通过线性估计和非线性优化等方法，利用至少4个特征点可以求解局部单应矩阵H_i。

通过以上分析，将标定板放置在参考坐标系上，可以得到局部单应矩阵H₁。然后，我们把其他局部单应矩阵H_i(i≠1)与H₁进行融合，得到整个场景平面(坐标系为)和图像平面之间的全局单应矩阵

为了将这些局部单应矩阵进行融合，首先我们需要知道各单应矩阵之间的关系。我们知道，P_ik表示标定板上的特征点在坐标系下的三维齐次世界坐标。令¹P_ik表示特征点在坐标系下的齐次世界坐标。通过几何分析，得到：

其中

¹θ_i表示坐标系和之间的旋转角度，[¹t_xi ¹t_yi]^T表示坐标系和之间的平移向量，表示归一化比例因子。

根据对局部单应矩阵的定义，得到：

p_ik＝λ_ikH_iP_ik (6)

p_ik＝λ_ikH₁·¹P_ik (7)

将式(4)代入式(7)，得：

结合式(6)和(8)，得：

此时，通过局部单应矩阵H₁和H_i可以求解矩阵¹M_i。因此，只要估计出局部单应矩阵，就可以求解矩阵¹M_i中的参数¹θ_i、¹t_xi和¹t_yi，即：

¹t_xi＝m₁₃,¹t_yi＝m₂₃ (11)

其中m_ij(i＝1,2,3；j＝1,2,3.)为矩阵¹M_i的第i行第j列元素。

第1.2，利用非线性最小二乘算法求解全局单应矩阵

利用1.1中的局部单应矩阵之间的关系以及相应的位姿参数，通过非线性优化方法即可将多个局部单应矩阵进行融合得到全局的单应矩阵，采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少1个局部单应矩阵最小化如下的目标函数J(·)，得到的解：

其中

为归一化比例因子。

优化后，得到融合N个局部单应矩阵的全局单应矩阵

第1.3全局单应矩阵标定过程

●局部单应矩阵的估计：将标定板放置在N个不同的位置，在每一个位置，利用传统的线性估计与非线性优化方法标定得到一个局部单应矩阵H_i。

●全局单应矩阵以及相关位姿参数的初值：将在第一位置处标定得到的局部单应矩阵H₁作为全局单应矩阵的初始值，的初始值通过式(9)、(10)、(11)计算得到。

●非线性优化：利用Levenberg-Marquardt(LM)算法[25]对目标函数(12)进行非线性优化，的最终值即为优化得到的全局单应矩阵。

第2，关联高度信息的单应矩阵的估计

第2.1，建立关联高度信息的单应矩阵与高度的关系

考虑静态场景中的N_h个特征点，X_0i和X_hi分别表示零平面和高度为h平面上对应点的齐次世界坐标：

X_0i＝[x_i y_i 0 1]^T,X_hi＝[x_i y_i h 1]^T (13)

相应的齐次图像坐标x_0i和x_hi为：

x_0i＝[u_0i v_0i 1]^T,x_hi＝[u_hi v_hi 1]^T (14)

根据摄像机小孔成像模型得：

x_0i＝λ_0iCX_0i (15)

x_hi＝λ_hiCX_hi (16)

其中λ_0i，λ_hi为归一化比例因子，C为3×4维矩阵，表示摄像机内外参数的乘积：

其中A_3×3表示摄像机内参数矩阵，R和t表示摄像机坐标系和世界坐标系之间的旋转平移关系。由于X_0i的第三列坐标为零，得到：

其中X_0i′＝[x_i y_i 1]^T，为图像平面和平面0之间的射影单应矩阵：

其中即为第1步中得到的全局单应矩阵，即c₁₁,c₁₂,c₁₄,c₂₁,c₂₂,c₂₄,c₃₁,c₃₂,c₃₄为已知参数。

根据式(13)中X_0i和X_hi的表达式，容易得到：

X_hi＝TX_0i (20)

其中，

根据式(16)和(20)，容易得到：

x_hi＝λ_hiCX_hi＝λ_hiCTX_0i＝λ_hiC′X_0i (22)

其中

因为X_0i的第三列坐标为零，重写式(22)得到：

x_hi＝λ_hiH′_hX_0i′ (24)

其中

根据式(18)和(24)得到：

其中G(h)为平面0和平面h上图像坐标间的关联高度信息的单应矩阵，其表达形式为：

其中全局单应矩阵可由第1步得到，因此矩阵C中除了参数c₁₃,c₂₃,c₃₃，其他均为已知。将式(19)和(25)代入式(27)，可知，关联高度信息的单应矩阵G(h)中只有三个未知参数，即c₁₃,c₂₃,c₃₃。

通过式(26)可知，一组点对(x_0i,x_hi)可以产生两组关于c₁₃,c₂₃,c₃₃的约束方程，所以给定高度h，至少两组点对就可以求解参数c₁₃,c₂₃,c₃₃。

第2.2，利用非线性最小二乘算法对关联高度信息的单应矩阵进行估计

采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少2组空间点对最小化如下目标函数J_h(·)，得到c₁₃,c₂₃,c₃₃的解：

其中

λ_i表示归一化比例因子，为c₁₃,c₂₃,c₃₃的最终优化结果。

最后利用c₁₃,c₂₃,c₃₃得到然后再结合全局单应矩阵将h设为变量，即得到关联高度信息的单应矩阵因此，对任意已知高度h，便可代入得到G(h)，将位于平面h上的点的图像坐标x_hi转换到平面0上的点的图像坐标x_0i，然后利用第1步得到的全局单应矩阵即可求得位于零平面上点的二维世界坐标。

第2.3关联高度信息的单应矩阵标定过程

为了得到关联高度信息的单应矩阵G(h)，我们用非线性优化的方法得到G(h)中的三个未知参数。给定高度h后，待优化的参数为c₁₃,c₂₃,c₃₃。标定过程如下：

●三个未知参数的初值：将高度为h₀的标定杆竖直放置在平面0上的两个不同的位置处，从而得到两组对应点的图像坐标x_0i和x_hi。将得到的图像坐标带入式(26)中，得到四个方程，求解得到未知参数的初始值。

●非线性优化：利用Levenberg-Marquardt(LM)算法对目标函数(28)进行非线性优化，当优化得到三个参数后，给定任意一个高度h，关联高度信息的单应矩阵G(h)就得到了。

得到平面h上特征点的图像坐标后，我们利用G(h)将其投影到相应平面0上点的图像坐标，然后再根据第1步得到的全局单应矩阵获得特征点的二维世界坐标x、y。

本发明的优点和有益效果

本发明提出了一种基于全局单应矩阵的单目视觉测量方法。本发明的主要特点如下：(1)在图像平面和场景平面之间，通过将多个局部单应矩阵进行融合获得精度更高的全局单应矩阵；(2)通过一个竖直杆对关联高度信息的单应矩阵进行标定，从而将某一高度平面上点的图像坐标转换到相应零平面上点的图像坐标。本发明的优点有：(1)本发明提出了一种精确的视觉测量方法，其精度接近图像中每像素所对应的实际物理距离；(2)在实际应用中，本发明简单可行；(3)本发明不需要摄像机内参数。实验结果证明本发明测测量精度较高并且容易实现。

附图说明：

图1为全局单应矩阵的标定过程示意图；

图2为关联高度信息的单应矩阵与各平面点的图像坐标之间的关系示意图；

图3为视觉测量系统的实物图；

图4为贴有反光膜的棋盘标定板；

图5为某一位置处的标定板图像；

图6为利用单应矩阵计算得到的水平杆的长度(实际长度为300mm)，其中■表示由全局单应矩阵计算得到的长度，*表示由局部单应矩阵计算得到的长度；

图7为利用单应矩阵计算得到的水平杆的长度(实际长度为1000mm)，其中■表示由全局单应矩阵计算得到的长度，*表示由局部单应矩阵计算得到的长度；

图8为用不同全局单应矩阵(由不同数目的局部单应矩阵融合得到)计算得到的水平杆的长度误差的均值和方差，其中■表示误差的均值，*表示误差的方差；

图9为改变水平杆放置位置的多组实验(标定板放置在不同的三个位置)中，水平杆测量误差的均值和方差，其中■表示误差的均值，*表示误差的方差，实线表示误差均值的均值，虚线表示误差方差的均值；

图10为利用关联高度信息的单应矩阵计算得到的投影点误差的均值和方差，其中■表示竖直杆高度为500mm时的情况，*表示竖直杆高度为600mm时的情况，◇表示竖直杆高度为700mm时的情况；

图11为改变标定杆放置位置的多组实验(标定杆放置在不同的四个位置)中，利用关联高度信息的单应矩阵计算得到的投影点平均误差，其中实线、虚线和点划线分别表示不同高度的投影点平均误差；

图12为对移动机器人的位姿进行三角测量和视觉测量的示意图；

图13为由三个标记点计算得到机器人位姿的过程示意图；

图14为由不同方法得到的移动机器人位置和姿态，其中○或□表移动机器人的位置，箭头表示移动机器人的方向，○表示用三角形的测量方法得到的结果，□表示本发明所提供的视觉测量方法得到的结果。

具体实施方式：

实施例1：

图1和图3展示了本发明用到的视觉测量系统，该系统由一个千兆以太网黑白摄像机组成，其视觉传感器为1/4-in(英寸)CMOS传感器，具有良好的动态性能，每秒可以采集120帧大小为640×480像素的图像。为了使目标更容易被检测到，我们选择了能够发射850nm红外波的红外发射器，并在目标位置处贴上反光膜，它可以将照射过来的红外光线反射回去，红外发射器紧临摄像机固定。同时，在摄像机镜头前端安装一个滤光片，保证只有波长为850nm的光线可以通过该摄像头。整个装置固定在三米高的金属杆顶端，如图3所示。

第1，全局单应矩阵的标定

第1.1各局部单应矩阵的标定

图4为标定过程中所使用的6×10的棋盘格，每一个格子的四角处被贴上圆形反光膜作为特征点。把标定板放置在摄像机视野范围内的25个不同的位置，并获得相应特征点的图像坐标。由于使用了红外装备，我们可以很方便地得到标定板上特征点的图像坐标，图5即为某位置处提取出的标定板上的特征点图像。在坐标系下，通过线性估计和非线性优化的方法，求解出25个不同位置处的局部单应矩阵H₁、H₂、…、H₂₅。

第1.2，建立各局部单应矩阵的之间的关系

根据空间几何关系，推导出在参考坐标系处获得的局部单应矩阵H₁与在其他位置处获得的局部单应矩阵H_i之间的关系为：

然后我们通过变换矩阵¹M_i，将坐标系下的坐标转换到坐标系下，从而可以建立多组约束方程：

由于N＝25、N_p＝77，根据式(8)可以得到N·N_p组约束方程，构成约束方程组。

第1.3，利用非线性最小二乘算法求解全局单应矩阵

采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少1个局部单应矩阵最小化如下目标函数J，得到全局单应矩阵的解：

其中

其中，为归一化比例因子，为参数¹θ_i,¹t_xi,¹t_yi的估计量。非线性优化的初始值选择为：优化后，得到融合25个局部单应矩阵的全局单应矩阵，如下：

第2，关联高度信息的单应矩阵的估计

取一高度为700mm的杆，其顶端和底端固定了两个反光点。首先，在摄像机视野范围内，将木杆竖直放置在20个不同的位置，得到20组顶端和底端的图像点对。

第2.1，建立关联高度信息的单应矩阵与高度的关系

考虑静态场景中的20个特征点，X_0i和X_hi分别表示零平面和高度为h＝700mm平面上对应点的齐次世界坐标：

X_0i＝[x_i y_i 0 1]^T,X_hi＝[x_i y_i h 1]^T (13)

相应的图像齐次坐标x_0i和x_hi为：

x_0i＝[u_0i v_0i 1]^T,x_hi＝[u_hi v_hi 1]^T (14)

推导出平面0和平面h上的点的图像坐标之间的关系为：

其中

c_ij(i＝1,2,3；j＝1,2,4)为矩阵的元素，它们通过第1步已经优化得到，而c₁₃、c₂₃和c₃₃为待求参数。

第2.2，利用非线性最小二乘算法对关联高度信息的单应矩阵进行估计

采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少2组空间点对最小化如下目标函数J_h，得到c₁₃,c₂₃,c₃₃的解：

其中

λ_i为归一化比例因子，为对变量c₁₃,c₂₃,c₃₃进行非线性优化后得到的最终结果，如下：

c₁₃＝-0.087,c₂₃＝-0.197,c₃₃＝0

代入式(27)，得到关联高度信息的单应矩阵如下：

当h＝700mm时，关联高度信息的单应矩阵为

第3，实验效果描述

第3.1，全局单应矩阵实验效果描述

第3.1.1全局单应矩阵的精度

为了验证全局单应矩阵比局部单应矩阵具有更高的准确度，我们分别采用单应矩阵和H₁对一定长度的木杆进行测量，该木杆的两端固定了两个反光点。我们将该杆水平放置在摄像机视野范围内的N_l个不同的位置，并得到杆上两端点对应的图像坐标，我们分别利用全局单应矩阵和局部单应矩阵得到杆上端点在坐标系下的世界坐标，并计算出该杆的长度。

首先我们取一个长度为300mm的杆，并将其放置在41个不同的位置，通过全局单应矩阵和局部单应矩阵，分别得到41个长度。通过计算，我们得到这些数据的绝对误差的均值和方差，如表1所示。图6给出了通过不同单应矩阵计算得到的不同长度，x轴表示杆放置的不同位置，y轴表示计算长度。在前10次放置过程中，我们将杆放置在坐标系原点附近，即在坐标系所在的局部单应矩阵区域，在后面的放置过程中，我们将杆放置在整个视野范围的其他位置。根据局部单应矩阵计算得到的长度，在坐标系原点附近，测量精度较好，但是在整个视野范围内，平均测量误差比根据全局单应矩阵计算得到的要大。因此，利用全局单应矩阵，平均误差更小，全局性能更好。

表1视觉测量绝对误差的均值和方差

为了进一步验证全局单应矩阵的性能，我们取另一个长度为1000mm的杆，然后分别利用局部单应矩阵和全局单应矩阵计算其长度，计算结果如图7所示。可以很清楚看到，在整个视野范围内，全局单应矩阵的计算精度比局部单应矩阵的高。

第3.1.2优化全局单应矩阵所需局部单应矩阵个数

在优化全局单应矩阵过程中，只需要融合三个局部单应矩阵即可得到高精度的全局单应矩阵。为了表明算法性能，进行了如下实验，首先将标定板放置在15个不同的位置，并得到15个不同的局部单应矩阵，通过将N(N＝1,2,…15)个局部单应矩阵进行融合，得到相应的全局单应矩阵为了验证这些全局单应矩阵的测量精度，同样利用长度为300mm的杆进行验证，我们将其水平放置在41个不同的位置，利用优化得到的全局单应矩阵进行测量，并计算其误差，其误差的均值和方差如图8所示。从图中可以看出，当局部单应矩阵的个数大于或者等于3时，误差的均值和方差均趋于稳定。为了说明该问题的一般性，我们又做了15组实验，每组实验中均将标定板放置在相距较远的三个位置，并用标定得到的3个局部单应矩阵进行优化，进而得到全局单应矩阵。同样利用上述实验所用的41组数据进行验证，结果如图9所示，从图中可以看出，测量误差的均值和方差均趋于稳定。因此，为了得到全局单应矩阵，只需将标定板放置在摄像机视野范围内的相距较远的三个位置即可，该方法简单可行。

第3.2，关联高度信息的单应矩阵实验效果描述

为了验证关联高度信息的单应矩阵的精度，取一高度为700mm的杆，其顶端和底端固定了两个反光点。将杆竖直放置在与标定位置不同的30个位置，得到每个位置顶端和底端点的图像坐标，用这些数据验证关联高度信息的单应矩阵的精度。为了验证用高度为700mm的杆标定出来的单应矩阵也同样适用于其它高度，又采用了高度为600mm和500mm的两个杆来进行验证。用式(26)计算出来的测量误差为：

其中表示顶端和底端点在x、y方向的误差。根据式(29)，我们很容易能够得到坐标误差的欧氏距离，利用该距离来评估关联高度信息的单应矩阵的精度：

进一步地，我们计算出该误差的平均误差和最大误差，如表2所示。从表中可以看出，平均测量误差大概在1像素附近波动，可以满足大部分视觉测量任务。

表2平均误差和最大误差

注：N_h表示标定杆放置位置的次数，E_A表示平均误差，E_M表示最大误差，单位：pixel

为了说明本发明提出的方法简单可行，我们将标定杆分别放置2次、3次、……、20次，并计算平均测量误差，结果如图10所示。从图中可以看出，当标定杆放置次数大于等于4时，平均测量误差趋于稳定。类似地，为了验证其一般性，我们又做了15组实验，每次实验中，将标定杆放置在不同的四个位置，平均测量误差如图11所示，可以看出，将标定杆放置在4个位置标定得到的关联高度信息的单应矩阵，其平均测量误差均近似趋于稳定。因此，在实际标定过程中，该方法易于实现，只需要将标定杆放置在摄像机视野范围内的4个不同位置，便能够标定得到较为准确的关联高度信息的单应矩阵。

第3.3，移动机器人的视觉位姿测量实验

前两部分实验分别对全局单应矩阵和关联高度信息的单应矩阵的精度进行了验证，这部分将前两部分结合在一起，对移动机器人进行视觉位姿测量。实验对象为先锋3-DX移动机器人，并在其高度为H_R的顶层平面放置三个反光点(三个点组成等腰三角形)。通过操作手柄将移动机器人移动到不同的位置，并用本发明的方法对其位置和姿态进行测量。

为了得到移动机器人的真实位姿，我们在参考坐标系的原点和其x轴上的某一点放置高度为H_R的基准点，使移动机器人上的任意一个反光点与两个基准点构成一个三角形，通过测量三角形三边的长度，来计算该反光点在坐标系下的世界坐标，图12描述了移动机器人上某一反光点的测量过程。具体来说，在点[0 0 H_R]^T和点[1000 0 H_R]^T(单位：mm)位置处放置两个基准点。三个反光点与两个基准点之间的距离a_i,b,c_i(i＝1,2,3)可以用直尺测量，然后反光点在坐标系下的x、y坐标为：

x_i＝a_i cosα_i (32)

y_i＝a_i sinα_i (33)

其中，a_i和c_i分别为反光点到两个基准点的距离，b为两基准点之间的距离，即b＝1000mm，α_i为三角形的边a_i和边b之间的夹角。

同时，我们通过本发明提供的视觉测量系统得到三个反光点的图像坐标，然后利用关联高度信息的单应矩阵将高度为H_R平面上的这三个点的图像坐标转换到零平面上，最后利用全局单应矩阵将零平面上的图像坐标转换为参考坐标系下的世界坐标。

由于移动机器人上的三个反光点组成了一个等腰三角形，所以得到三个点的世界坐标之后，可以得到移动机器人的位置和姿态角，如图13所示，其结果表达式为：

其中s＝[x_t y_t θ_t]^T为移动机器人的位姿。

实验得到了移动机器人在10个不同位置处的位姿结果，并计算出其误差的均值和方差如表3所示，并将通过本发明得到的移动机器人位姿和通过三角形计算得到的移动机器人位姿用图14描述出来。从图中可以看出，由本发明得到的位姿结果与真实测量的结果非常接近。由于图像分辨率的限制，摄像机每像素的测量精度A(单位：mm/pixel)可以由下式计算得到：

表3移动机器人测量位姿的绝对误差的均值和方差

其中γ_H×γ_V为测量范围(单位：mm×mm)，R_H×R_V为图像分辨率(单位：pixel×pixel)。采用运算符max{·,·}是为摄像机测量精度提供一个更可靠的评价标准。在本实验中，通过尺子的粗略测量，γ_H×γ_V＝4700mm×3500mm，摄像机的分辨率为R_H×R_V＝640pixel×480pixel，则摄像机每像素的测量精度为：

通过表3可以看出，本发明所提出方法的测量精度与摄像机每像素的测量精度非常接近。对于一般移动机器人，通过本发明测量得到的位置和姿态精度已经足够高了。如果需要更高的精度需求，可以选择一个更高分辨率或者更小视野范围的摄像机。

参考文献

1.C.Chen,D.Schonfeld,Pose estimation from multiple cameras based onsylvesters equation,Computer Vision and Image Understanding 114(6)(2010)652–666.

2.W.Wolfe,D.Mathis,C.Sklair,M.Magee,The perspective view of threepoints,IEEE Trans.Pattern Analysis and Machine Intelligence 13(1)(1991)66–73.

3.J.Kelsey,J.Byrne,M.Cosgrove,S.Seereeram,R.Mehra,Vision-basedrelative pose estimation for autonomous rendezvous and docking,in:Proceedingsof 2006 IEEE Conference on Aerospace Conference,2006.

4.M.Knorr,W.Niehsen,C.Stiller,Online extrinsic multi-cameracalibration using ground plane induced homographies,in:IEEE IntelligentVehicles Symposium(IV),2013,pp.236–241.

5.F.Zhu,L.He,Y.Hao,et al.,Contrast of calculated accuracy betweenmonocular vision and binocular vision without modeling,Chinese Journal ofScientific Instrument 28(4)(2007)165–169.

6.X.Xiao,Study on the monocular vision measurement methodPh.D.thesis,Tianjin University(2009).

7.J.Li,Q.Ding,P.Wang,Monocular vision measurement method of objectposition and orientation based on quaternions,Aeronautical ComputingTechnique 6(2011)006.

8.M.A.Fischler,R.C.Bolles,Random sample consensus:a paradigm formodel fitting with applications to image analysis and automated cartography,Communications of the ACM 24(6)(1981)381–395.

9.R.Hartley,A.Zisserman,Multiple view geometry in computer vision,Cambridge university press,2003.

10.R.Lakemond,C.Fookes,S.Sridharan,Evaluation of two-view geometrymethods with automatic ground-truth generation,Image and Vision Computing 31(12)(2013)921–934.

11.X.Gao,X.Hou,J.Tang,H.Cheng,Complete solution classification forthe perspective-three-point problem,IEEE Trans.Pattern Analysis and MachineIntelligence 25(8)(2003)930–943.

12.L.Kneip,D.Scaramuzza,R.Siegwart,A novel parametrization of theperspective-three-point problem for a direct computation of absolute cameraposition and orientation,in:Proceedings of 2011 IEEE Conference on ComputerVision and Pattern Recognition(CVPR),2011,pp.2969–2976.

13.Z.Hu,F.Wu,A note on the number of solutions of the noncoplanar p4pproblem,IEEE Trans.Pattern Analysis and Machine Intelligence 24(4)(2002)550–555.

14.L.Quan,Z.Lan,Linear n-point camera pose determination,IEEETrans.Pattern Analysis and Machine Intelligence 21(8)(1999)774–780.

15.A.Ansar,K.Daniilidis,Linear pose estimation from points or lines,IEEE Trans.Pattern Analysis and Machine Intelligence 25(5)(2003)578–589.

16.D.Xu,Y.F.Li,M.Tan,A general recursive linear method and uniquesolution pattern design for the perspective-n-point problem,Image and VisionComputing 26(6)(2008)740–750.

17.T.Scoleri,Video metrology without the image-to-ground homography,in:Preceedings of 2010 International Conference on Digital Image Computing:Techniques and Applications (DICTA),2010,pp.335–342.

18.G.Wang,Y.Wu,Z.Hu,A novel approach for single view based planemetrology,in:Proceedings of 16th International Conference on PatternRecognition,Vol.2,2002,pp.556–559.

19.G.Wang,Z.Hu,F.Wu,H.-T.Tsui,Single view metrology from sceneconstraints,Image and Vision Computing 23(9)(2005)831–840.

20.G.Wang,H.-T.Tsui,Z.Hu,F.Wu,Camera calibration and 3dreconstruction from a single view based on scene constraints,Image and VisionComputing 23(3)(2005)311–323.

21.C.Meng,N.Zhou,X.Xue,Y.Jia,Homography-based depth recovery withdescent images,Machine Vision and Applications 24(5)(2013)1093–1106.

22.J.Li,X.Xie,Study of monocular measuring technique based onhomography matrix,in:Computing and Intelligent Systems,2011,pp.316–324.

23.C.Sun,Z.Zhang,P.Wang,Measuring method for the object pose based onmonocular vision technology,in:Photonics Asia 2010,2010,pp.78550Z–78550Z.

24.T.Elamsy,A.Habed,B.Boufama,Self-calibration of stationarynonrotating zooming cameras,Image and Vision Computing 32(3)(2014)212–226.

25.K.Madsen,H.B.Nielsen,O.Tingleff,Methods for non-linear leastsquares problems,1999.

Claims (1)

1.一种基于全局单应矩阵的单目视觉测量方法，其特征在于该方法包括：

第1，全局单应矩阵的标定

定义“局部单应矩阵”为标定板在待测地面单个位置处，标定得到的与图像平面之间的单应矩阵；这里，采用“局部”，是因为标定板的尺寸比较小，只能覆盖整个视野范围内待测地面的小部分区域；相应地，我们定义“全局单应矩阵”，用来表示将标定板放置在不同位置得到多幅图像，并将这些位置的数据进行融合得到的单应矩阵；“局部单应矩阵”反映了图像平面和局部区域之间的映射关系，而“全局单应矩阵”，利用不同位置的数据，更准确地反映了图像平面和整个场景平面之间的映射关系；

将单目视觉测量系统中的一个未标定的摄像机固定在某高度的金属杆顶端，将标定板放置在摄像机视野范围内待测地面即零平面的不同位置处；以位置i处的标定板左下角为原点建立坐标系选择标定板在第一个位置处所建立的坐标系为参考世界坐标系；

第1.1，建立各局部单应矩阵之间的关系

用局部单应矩阵H_i建立特征点图像坐标与世界坐标之间的关系如下：

p_ik＝λ_ikH_iP_ik

其中，λ_ik表示归一化比例因子，P_ik＝[x_ik y_ik 1]^T表示标定板上特征点在参考世界坐标系X_wY_w平面上的二维齐次世界坐标，p_ik＝[u_ik v_ik 1]^T表示相应的齐次图像坐标，i＝1,2,…,N；k＝1,2,…,N_p，N表示标定板的放置次数，N_p表示标定板上特征点的个数；

根据空间几何关系，推导出在参考坐标系处获得的局部单应矩阵H₁与在其他位置处获得的局部单应矩阵H_i之间的关系为：

${M_1}_i=\frac{1}{{\mathrm{\λ}}_{M_i}}{H_1}^{-1}{H}_i$

其中，为坐标系在坐标系下的变换矩阵，¹θ_i表示两坐标系之间的旋转角度，[¹t_xi ¹t_yi]^T表示两坐标系之间的平移向量，表示归一化比例因子；

这样，根据任意一个局部单应矩阵和参考坐标系处的单应矩阵，可以求出它们之间的变换矩阵¹M_i，进而求出旋转与平移参数¹θ_i、¹t_xi和¹t_yi；

然后通过变换矩阵¹M_i，将坐标系下的坐标转换到坐标系下，从而建立多组约束方程：

$p_{ik}={\mathrm{\λ}}_{ik}{\mathrm{\λ}}_{M_i}{H}_1\·{M_1}_i\·P_{ik}$

由于i＝1,2,…,N；k＝1,2,…,N_p，根据式(8)可以得到N·N_p组约束方程，构成约束方程组；

第1.2，利用非线性最小二乘算法求解全局单应矩阵

采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少1个局部单应矩阵最小化如下目标函数J，得到全局单应矩阵的解：

$J({\hat{H}}_g,{\widehat{\mathrm{\θ}}_1}_i,{\hat{t}_1}_{xi},{\hat{t}_1}_{yi})=\underset{i=1}{\overset{N}{\Σ}}\underset{k=1}{\overset{N_p}{\Σ}}{(p_{ik}-s_{ik}{\hat{H}}_g\·{M_1}_i({\widehat{\mathrm{\θ}}_1}_i,{\hat{t}_1}_{xi},{\hat{t}_1}_{yi})\·P_{ik})}^2$

其中

${M_1}_i({\widehat{\mathrm{\θ}}_1}_i,{\hat{t}_1}_{xi},{\hat{t}_1}_{yi})=\left[\begin{array}{ccc}\cos \left({\widehat{\mathrm{\θ}}_1}_i\right)& -\sin \left({\widehat{\mathrm{\θ}}_1}_i\right)& {\hat{t}_1}_{xi}\\ \sin \left({\widehat{\mathrm{\θ}}_1}_i\right)& \cos \left({\widehat{\mathrm{\θ}}_1}_i\right)& {\hat{t}_1}_{yi}\\ 0& 0& 1\end{array}\right]$

其中，为归一化比例因子，为参数¹θ_i,¹t_xi,¹t_yi的估计量；非线性优化的初始值选择为：优化后，得到融合N个局部单应矩阵的全局单应矩阵

第2，关联高度信息的单应矩阵的估计

由第1步，得到了全局单应矩阵利用该全局单应矩阵，能够根据图像计算得到零平面上特征点在参考世界坐标系中的世界坐标；然而，在很多实际应用中，被测量的特征点并不在零平面上，而是在某一特定高度的平面上；为进一步解决这个问题，应估计出关联高度信息的单应矩阵，该单应矩阵描述了位于零平面即平面0上的点的图像坐标与位于高度为h的平面即平面h上的点的图像坐标之间的关系；一旦该关联高度信息的单应矩阵被标定，给定高度h，就可以通过该单应矩阵将位于平面h上的点的图像坐标x_hi转换到平面0上的点的图像坐标x_0i，然后利用第1步得到的全局单应矩阵即可求得位于零平面上点的二维世界坐标；

第2.1，建立关联高度信息的单应矩阵与高度的关系

考虑静态场景中的N_h个特征点，X_0i和X_hi分别表示平面0和平面h上对应点的齐次世界坐标：

X_0i＝[x_i y_i 0 1]^T,X_hi＝[x_i y_i h 1]^T

相应的图像齐次坐标x_0i和x_hi为：

x_0i＝[u_0i v_0i 1]^T,x_hi＝[u_hi v_hi 1]^T

推导出平面0和平面h上的点的图像坐标之间的关系为：

$x_{0i}=\frac{{\mathrm{\λ}}_{0i}}{{\mathrm{\λ}}_{hi}}G\left(h\right)x_{hi}$

其中，λ_0i，λ_hi为归一化比例因子，

$\begin{array}{c}G\left(h\right)={\hat{H}}_g{H}_h^{\′-1}\\ =\left[\begin{array}{ccc}{c}_{11}& c_{12}& c_{14}\\ c_{21}& c_{22}& c_{24}\\ c_{31}& c_{32}& c_{34}\end{array}\right]{\left[\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}h+c_{14}\\ c_{21}& c_{22}& c_{23}h+c_{24}\\ c_{31}& c_{32}& c_{33}h+c_{34}\end{array}\right]}^{-1}\end{array}$

c_ij为矩阵的元素，i＝1,2,3；j＝1,2,4，它们通过第1步已经优化得到，而c₁₃、c₂₃和c₃₃为待求参数；若高度h已知，给定至少两对点，可以得到至少4个约束方程，即可对参数c₁₃、c₂₃和c₃₃进行优化求解；

第2.2，利用非线性最小二乘算法对关联高度信息的单应矩阵进行估计

采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少两对点最小化如下目标函数J_h，得到c₁₃,c₂₃,c₃₃的解：

$J_h({\hat{c}}_{13},{\hat{c}}_{23},{\hat{c}}_{33})=\underset{i=1}{\overset{n}{\Σ}}{(x_{0i}-{\mathrm{\λ}}_i{\hat{H}}_g{\hat{H}}_h^{\′-1}({\hat{c}}_{13},{\hat{c}}_{23},{\hat{c}}_{33})\·x_{hi})}^2$

其中

${\hat{H}}_h^{\′}({\hat{c}}_{13},{\hat{c}}_{23},{\hat{c}}_{33})=\left[\begin{array}{ccc}{c}_{11}& c_{12}& {\hat{c}}_{13}h+c_{14}\\ c_{21}& c_{22}& {\hat{c}}_{23}h+c_{24}\\ c_{31}& c_{32}& {\hat{c}}_{33}h+c_{34}\end{array}\right]$

λ_i为归一化比例因子，为对变量c₁₃,c₂₃,c₃₃进行非线性优化后得到的最终结果；

最后利用c₁₃,c₂₃,c₃₃得到然后再结合全局单应矩阵将h设为变量，即得到关联高度信息的单应矩阵因此，对任意已知高度h，便可代入得到G(h)，将位于平面h上的点的图像坐标x_hi转换到平面0上的点的图像坐标x_0i，然后利用第1步得到的全局单应矩阵即可求得位于零平面上点的二维世界坐标。