CN102750704B - 一种摄像机分步自标定方法 - Google Patents

Abstract

本发明提供一种摄像机分步自标定方法，包括：第1阶段，对同一场景拍摄两张在不同焦距下的图像，对图像进行特征提取和匹配并剔除误匹配点，利用匹配点对求解摄像机的主点坐标；第2阶段，从不同角度对同一场景进行拍摄获得三幅可供匹配的图像，对图像进行特征提取和匹配并剔除误匹配点，然后基于Kruppa方程，将已获得的主点坐标代入方程实现倾斜因子以及摄像机在成像平面u、v轴方向上的尺度因子三个参数的求解。通过本发明的方法，摄像机的主点坐标的标定精度较高，并减小了Kruppa方程系数矩阵规模，降低了求解运算量，具有实时性。该方法适用于视觉系统的摄像机标定，可在三维测量、三维重建、机器导航以及增强现实等领域使用。

Description

一种摄像机分步自标定方法

技术领域

本发明涉及计算机视觉领域，更具体的说，涉及摄像机自标定方法。

背景技术

随着信息和通信技术的不断发展，计算机视觉中三维重建以及增强现实技术得到了越来越广泛的应用。由于计算机获取的都是二维的图像图形信息，缺乏必要的真实感，三维重建技术则可以通过利用场景的两幅、多幅图像或者一段视频序列恢复场景的三维立体结构，使人们能够更加全面地认识和把握真实的事物；增强现实则是通过计算机技术将虚拟物体准确的放置在真实的环境中，并且借助各种显示设备实现虚拟世界与现实世界的无缝融合，从而呈现给用户一个视觉效果真实的新环境。三维重建以及增强现实技术在实际中应用范围十分广大，如在现代工业中，轿车车身、凸轮轴、旋转机械叶片等的表层曲面难于测量设计，使用三维重建技术则能够方便快速的重构出这些结构和曲面；增强现实则可以在电视转播的过程中对视频数据进行补充，添加虚拟的广告信息，然后将融合后的视频数据呈现在电视屏幕上，从而有效降低广告的成本。一个三维重建、增强现实系统的实现，摄像机标定是其中的重要环节，是三维重建、增强现实实现的基础。在三维重建中，摄像机自标定将确定如何从多幅图像中恢复摄像机的运动参数(包括位姿状态)、自身成像有关参数(如焦距等)和场景结构的方法，而在增强现实中，摄像机自标定可为虚实配准提供准确的数据来源，是一个增强现实系统必不可少的关键技术。

摄像机标定的基本任务之一是从摄像机获取的图像信息出发计算三维空间中物体的几何信息，并由此重建和识别物体，而空间物体表面某点的三维几何位置与其在图像中对应点之间的相互关系是由摄像机成像几何模型决定的，这些几何模型参数就是摄像机参数。在大多数条件下，这些参数必须通过实验与计算才能得到，这个过程称为摄像机标定，标定过程就是确定摄像机的几何和光学参数以及摄像机相对于世界坐标系的方位。摄像机标定技术十分重要，摄像机标定的参数精度直接影响后期三维重建的精度，因此摄像机标定技术一直是计算机视觉领域里研究的基本任务之一。

摄像机标定通常采用经典的针孔摄像机成像模型，该模型的数学表达如下公式所示：

其中，(x_w y_w z_w 1)^T是空间中任一点P的齐次坐标，(u v 1)^T是图像上对应空间点(x_w y_w z_w 1)^T像点的齐次坐标，α是一比例因子，矩阵K为内参数矩阵，(u₀ v₀)^T是摄像机主点的坐标，即摄像机光轴与图像平面交点的图像坐标，f_u，f_v是摄像机图像u方向和v方向上物理坐标到图像像素坐标的比例系数，R是一个3×3的单位正交矩阵，t是一个平移向量，而矩阵M描述了摄像机相对于世界坐标系的位姿状态。

一般来说，摄像机标定方法根据有无标定物可分为传统的摄像机标定方法和自标定方法。传统的标定方法主要是利用景物上点的世界坐标和其图像上对应点坐标之间的对应关系，建立约束方程，然后通过求解该约束方程和优化方法得到摄像机内外参数。传统的标定方法中有“Geometric cameracalibration using circular control points”(Heikkila J. IEEETransactions on Pattern Analysis and Machine Intelligence，2000，22(10)：1066-1077)”，“Accurate Camera Calibration from Multi-ViewStereo and Bundle Adjustment”(Furukawa Y，Ponce J.InternationalJournal of Computer Vision，2009，84(3)：257-268)，“Single View PointOmnidirectional Camera Calibration from Planar Grids”(Mei C，RivesP.IEEE International Conference on Robotics and Automation，Italy，2007：3945-3950)。由于传统的标定方法需要使用标定物，因而对环境条件有所限制。如果有办法可以直接对三维场景进行标定，那么应用的范围将会进一步扩大。于是Faugeras、Hartley等在20世纪90年代初首次提出了摄像机自标定的概念，使得摄像机参数的标定可以不受场景位置和摄像机运动的限制。摄像机自标定方法是仅依靠多幅图像之间的对应关系，这种关系可用Kruppa方程进行数学建模表示，然后通过求解Kruppa方程得到摄像机内参数。自标定方法中有“Camera Self-Calibration for SequentialBayesian Structure from Motion”(Javier C，Diana R B，Andrew J D etal.IEEE International Conference on Robotics and Automation，Japan，2009：3411-3416)，“Continuous Stereo Self-Calibration by CameraParameter Tracking”(Thao D，Christian H，Christoph S.IEEETransactions on Image Processing，2009，18(7)：1536-1550)，“CameraSelf-Calibration from Bivariate Polynomials Derived from Kruppa′sEquations”(Habed A.Pattern Recognition，2008，41(8)：2484-2492)等。

国内的计算机视觉研究起步相对较晚，但近年来在基于主动视觉的摄像机自标定技术上发表了大量研究成果。其中，马颂德(A Self-calibrationTechnique for Active Vision System，IEEE Transactions on Robotics andAutomation，1996，12(1)：114-120)提出了一种利用摄像机做两组三正交平移运动的自标定方法，这种方法对摄像机运动要求苛刻，在普通的Pan-Tilt-Translation主动视觉平台上很难实现；杨长江等(一种基于主动视觉的摄像机内参数自标定方法，计算机学报，1998，21(5)：428-435)对此进行改进，提出了一种利用四组平移运动(其中每组运动包括两次相互正交的平移运动)线性标定摄像机内参数的方法；雷成(一种新的基于主动视觉系统的摄像机自标定方法，计算机学报，2000，23(11)：1130-1136)则指出利用摄像机做四次平移运动(其中任意三次不在同一平面)即可线性求解摄像机的全部五个内参数；此外曲学军等提出了一种“基于平行直线束图像序列的摄像机标定”方法(曲学军，张璐.计算机测量与控制，2010，18(6)：1421-1427)；朱嘉等利用摄像机做一组二维平移运动标定参数的方法(摄像机的一种主动视觉标定方法，光学学报，2010，30(5)：1297-1303)。由于这些方法需要对摄像机运动做严格的限制，故实际应用范围并不广泛。

综上，传统的摄像机标定方法虽然标定精度高，但标定物的需要影响了应用范围。基于主动视觉的自标定技术需对摄像机的运动做严格的限制，因而无法满足实时标定的要求。而所有的自标定方法本质上都是基于摄像机内参数自身存在的约束关系而对摄像机进行标定的，这种约束关系只与摄像机内参数有关，而与摄像机的运动姿态以及所在场景的结构无关，因此自标定方法要比传统的标定方法更为灵活。

基于Kruppa方程的摄像机自标定技术主要是通过建立图像序列的Kruppa方程组，然后根据Kruppa方程的约束条件将摄像机的五个内参数统一求解出来。现有自标定方法大多是通过优化方法求解Kruppa方程，Kruppa方程的求解仍是一个非线性问题，且标定的结果易受图像噪声的干扰，标定过程耗时较长，且标定结果的鲁棒性易受噪声的影响。

发明内容

针对上述问题，本发明的研究重点在于基于Kruppa方程的摄像机自标定方法。

本发明提出了一种摄像机分步自标定方法，所述摄像机采取的成像模型矩阵为： $K=\left(\begin{array}{ccc}{f}_u& s& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right),$ 其中，(u₀，v₀)是摄像机主点坐标，s是倾斜因子，f_u、f_v分别为摄像机在成像平面u轴、v轴方向上的尺度因子；其特征在于，所述方法分为2个阶段进行：

第1阶段，标定摄像机的主点坐标(u₀，v₀)；

第2阶段，结合第1阶段获得的主点坐标，标定摄像机其他三个参数，即倾斜因子s以及摄像机在成像平面u轴、v轴方向上的尺度因子f_u、f_v；

其中在第1阶段中，通过调节摄像机的焦距，同时保持摄像机的位姿状态不变，对同一场景拍摄两张在不同焦距下的图片，对图像进行特征提取和匹配，并剔除误匹配点，最后利用匹配点对求解摄像机的主点坐标；在第2阶段中，控制摄像机的内参数不变，从不同角度对同一场景进行拍摄获得三幅含有足够多特征点的可供匹配的图像，对图像进行特征提取和匹配，并剔除误匹配点，然后基于Kruppa方程，确定Kruppa方程的比例因子初始值并利用遗传算法对其进行优化，通过优化后的比例因子以及第1阶段获得的主点坐标代入Kruppa方程实现摄像机其他三个参数的求解。

进一步，所述第1阶段的具体流程包括：

步骤1.1：保持摄像机的位姿状态不变，调整摄像机的焦距，获得摄像机对同一场景分别在焦距f₁和f₂条件下的两张图像；

步骤1.2：利用SURF算法对获取的两张图像进行特征提取，然后对两图像所含的特征点进行初匹配操作，得到两图像特征点的匹配对集合；

步骤1.3：采用对极几何约束的RANSAC算法对步骤1.2中的匹配对集合进行误匹配剔除操作，将误匹配点对从匹配点对集合中去除，然后将该匹配点对集合保存起来；

步骤1.4：判断步骤1.3中所获得的匹配点对数P是否大于第一阈值T1若大于T1，则进行步骤1.5操作，否则转入步骤1.1，其中T1≥8；

步骤1.5：将步骤1.4所获得的P个匹配特征点对写成如下的方程组形式：

$\left\{\begin{array}{c}({y_1}^{\′}-y_1)u_0+(x_1-{x_1}^{\′})v_0=x_1{y_1}^{\′}-{x_1}^{\′}{y}_1\\ ({y_2}^{\′}-y_2)u_0+(x_2-{x_2}^{\′})v_0=x_2{y_2}^{\′}-{x_2}^{\′}{y}_2\\ .......\\ ({y_P}^{\′}-y_P)u_0+(x_P-{x_P}^{\′})v_0=x_P{y_P}^{\′}-{x_P}^{\′}{y}_P\end{array}\right.,$

利用最小二乘法求解该方程组，得到的解(u₀，v₀)便是摄像机的主点坐标，其中，(x_i，y_i)和(x′_i，y′_i)是两图像的匹配特征点对集合中第i对特征点坐标，i＝1，2，...，P。

进一步骤，所述第2阶段的具体流程包括：

步骤2.1：保持摄像机内参数不变，通过调整摄像机的位姿状态，从不同角度对同一场景拍摄图像，获得至少三幅图像；

步骤2.2：利用SURF算法对步骤2.1获得的图像进行特征提取，并对两两图像中的特征点进行初始匹配操作；

步骤2.3：采用对极几何约束的RANSAC算法对步骤2.2中的特征点对集合进行误匹配剔除操作，并将正确的匹配点对保存起来；

步骤2.4：判断步骤2.3中所获得的两两图像匹配点对数p是否大于第二阈值T2，若大于T2，则进行步骤2.5操作，否则舍掉该图像对，并返回步骤2.1，重复执行2.1至步骤2.4得到满足条件的图像对，T2≥8；

步骤2.5：利用步骤2.3中所获得的两两图像匹配点对计算该图像对间的基本矩阵和外极点的反对称矩阵，然后构造该图像对的Kruppa方程，进而构造所有图像对的Kruppa方程组；

步骤2.6：利用公式求出所有图像对Kruppa方程的比例因子初始值，其中 是图像对间基本矩阵F的两个特征值，是外极点e′的反对称矩阵的特征值；

步骤2.7：利用遗传算法对所有图像对的Kruppa方程比例因子进行优化；

步骤2.8：将优化后的比例因子以及第1阶段获得的摄像机主点坐标代入步骤2.5构造的Kruppa方程组，然后利用奇异值分解法以及cholesky分解，最终求出摄像机内其他三个内参数。

进一步，所述T1＝10，T2＝10。

进一步，所述对极几何约束的RANSAC算法具体包括：

步骤1：给定两幅图像I₁、I₂，将通过所述匹配方法获得的图像点对集合进行坐标归一化处理，并初始化循环次数i＝0；

步骤2：从点对集合中随机取8对点，然后利用8点算法计算基本矩阵F_i；

步骤3：使用sampson距离计算公式 $d=\frac{m^{\′T}{F}_{i}m}{{\left(F_{i}m\right)}_1^2+{\left(F_{i}m\right)}_2^2+{\left(F_i^T{m}^{\′}\right)}_1^2+{\left(F_i^T{m}^{\′}\right)}_2^2},$ 在匹配点对集合中寻找所有符合d＜t条件的点对，将它们做为内点，记录每一次的内点数量Q_i，其中t取(0.001，0.01)之间的数；

步骤4：i＝i+1，若i＜N，则重复执行步骤2和步骤3，并记录每一次的内点数量，否则转入步骤5；其中P代表所期望达到的概率，e是正确的数据比例，每次计算时不断更新；

步骤5：从步骤3求得的Q₁，Q₂，……，Q_N中选取最大的内点数Q_i，将其对应的基本矩阵F_i作为最终要求解的基本矩阵F，再以sampson距离为依据，寻找匹配点对集合中所有满足d＜t的点对，将它们做为最终的正确的匹配点对，其它的点对予以剔除。

通过本发明提出的方案，该方法无需任何标定物，只需通过摄像机获取的图像序列，就可以完成摄像机内参数的标定。本发明的方法对摄像机的主点进行单独标定，只需两幅图像，过程极为简单且主点的标定精度较高，同时在第2阶段中，由于代入了已获得的主点坐标，Kruppa方程系数矩阵规模减小，降低了求解运算量，因此整个标定过程较直接求解Kruppa方程的自标定技术耗时更短，具有实时性，且通过仿真实验结果表明本发明的标定结果精度较高，鲁棒性好，具有很好的实用性。该方法适用于视觉系统的摄像机标定，可在三维测量、三维重建、机器导航以及增强现实等领域使用。

附图说明

图1为本发明的摄像机分步自标定方法的流程图。

图2为线性摄像机成像几何模型的示意图。

图3为摄像机主点坐标的标定流程图。

图4为摄像机其它三个内参数的标定流程图。

图5为两视图极几何成像模型的示意图。

图6为对极几何约束的RANSAC算法流程图。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

本文采用摄像机成像模型矩阵如下：

$K=\left(\begin{array}{ccc}{f}_u& s& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right)$

其中，(u₀，v₀)是摄像机主点坐标，s是倾斜因子，f_u、f_v分别为摄像机在成像平面u轴、v轴方向上的尺度因子，它们统称为摄像机的内参数。

本发明提供一种摄像机分步自标定方法，包括：第1阶段，对同一场景拍摄两张在不同焦距下的图像，对图像进行特征提取和匹配并剔除误匹配点，利用匹配点对求解摄像机的主点坐标；第2阶段，从不同角度对同一场景进行拍摄获得三幅可供匹配的图像，对图像进行特征提取和匹配并剔除误匹配点，然后基于Kruppa方程，将已获得的主点坐标代入方程实现倾斜因子以及摄像机在成像平面u、v轴方向上的尺度因子三个参数的求解。

一种较为具体的实施例为，本发明从基于Kruppa方程的摄像机自标定技术入手，提出了一种分步自标定方法，其将摄像机内参数标定过程分为两步。首先标定摄像机的主点坐标，通过调节摄像机的焦距，同时保持摄像机的位姿状态不变，对同一场景拍摄两张在不同焦距下的图像，然后利用SURF算法提取两图像的特征点及匹配，再通过对极几何约束的RANSAC算法剔除误匹配点，最后利用匹配点对求解摄像机的主点坐标，此过程只需要两幅图片即可完成主点的标定过程，操作简单，方便实用；其次标定摄像机其他三个内参数，控制摄像机的内参数不变，从不同角度对同一场景进行拍摄，然后选取三幅含有足够多特征点的可供匹配的图像，同样采用SURF算法对图像进行特征提取和匹配，并通过对极几何约束的RANSAC算法剔除误匹配点，然后从基于Kruppa方程的摄像机自标定技术入手，确定Kruppa方程的比例因子初始值并利用遗传算法对其进行优化，通过优化后的比例因子实现摄像机其他三个内参数的求解，该过程的特点在于首先将第一步标定的主点代入Kruppa方程中，使得Kruppa方程的未知数由5个减少至3个，Kruppa方程的系数矩阵也由3×6变为3×4，降低了运算量。

以下结合附图来予以说明。图1示出了本发明摄像机分步自标定的方法。具体包括：

第1阶段，标定摄像机的主点坐标(u₀，v₀)；

第2阶段，标定摄像机其他三个参数，即倾斜因子s以及摄像机在成像平面u轴、v轴方向上的尺度因子f_u、f_v。

首先，说明本发明的第1阶段，摄像机主点坐标(u₀，v₀)的标定。

图2示出了线性摄像机成像几何模型，其中P(u₀，v₀)是摄像机的主点，f是摄像机的焦距。

由摄像机成像模型可知，当焦距f变化时，图像会发生缩放，但摄像机光轴是不变的，因此图像的主点也不会发生变化，所以可通过调节摄像机焦距得到同一场景在不同焦距下拍摄的两幅图像，然后根据下式计算摄像机主点坐标：

$\frac{x^{\′}-u_0}{x-u_0}=\frac{y^{\′}-v_0}{y-v_0}$

$\⇒(y^{\′}-y)u_0+(x-x^{\′})v_0=x{y}^{\′}-x^{\′}y$

其中(x，y)和(x′，y′)分别是两幅图像上匹配的特征点对坐标。上式为线性方程，可选择多个特征点对，然后利用最小二乘法求解摄像机主点坐标(u₀，v₀)。

由此，通过调节摄像机的焦距，同时保持摄像机的位姿状态不变，对同一场景拍摄两张在不同焦距下的图片，然后利用SURF算法提取两图像的特征点及匹配，再通过对极几何约束的RANSAC算法剔除误匹配点，最后利用匹配点对求解摄像机的主点坐标。此过程只需要两幅图片即可完成主点的标定过程，操作简单，方便实用。具体流程如图3所示，包括如下步骤：

步骤1.1：保持摄像机的位姿状态不变，调整摄像机的焦距，获得摄像机对同一场景分别在焦距f₁和f₂条件下的两张图像；

步骤1.2：利用SURF算法对获取的两张图像进行特征提取，然后对两图像所含的特征点进行初匹配操作，得到两图像特征点的匹配对集合；

步骤1.3：采用对极几何约束的RANSAC算法对步骤1.2中的匹配对集合进行误匹配剔除操作，将误匹配点对从匹配点对集合中去除，然后将该匹配点对集合保存起来；

步骤1.4：判断步骤1.3中所获得的匹配点对数P是否大于第一阈值T1若大于T1，则进行步骤1.5操作，否则转入步骤1.1，其中T1≥8，优选地，可以取T1＝10；

步骤1.5：将步骤1.4所获得的P个匹配特征点对写成如下的方程组形式：

$\left\{\begin{array}{c}({y_1}^{\′}-y_1)u_0+(x_1-{x_1}^{\′})v_0=x_1{y_1}^{\′}-{x_1}^{\′}{y}_1\\ ({y_2}^{\′}-y_2)u_0+(x_2-{x_2}^{\′})v_0=x_2{y_2}^{\′}-{x_2}^{\′}{y}_2\\ .......\\ ({y_P}^{\′}-y_P)u_0+(x_P-{x_P}^{\′})v_0=x_P{y_P}^{\′}-{x_P}^{\′}{y}_P\end{array}\right.,$

其中，(x_i，y_i)和(x′_i，y′_i)是两图像的匹配特征点对集合中第i(i＝1，2，...，P)对特征点坐标。此时可利用最小二乘法求解该方程组，最终得到的解(u₀，v₀)便是摄像机的主点坐标。

完成摄像机主点的标定后，本发明进入第2阶段，对摄像机其他三个内参数的标定工作。这一过程主要是根据基于Kruppa方程的摄像机自标定方法进行的。首先控制摄像机的内参数不变，从不同角度对同一场景进行拍摄，然后选取三幅含有足够多特征点的可供匹配的图像，同样采用SURF算法对图像进行特征提取和匹配，并通过对极几何约束的RANSAC算法剔除误匹配点，然后从基于Kruppa方程的摄像机自标定技术入手，确定Kruppa方程的比例因子初始值并利用遗传算法对其进行优化，通过优化后的比例因子实现摄像机其他三个内参数的求解，该过程的特点在于首先将第一步标定的主点代入Kruppa方程中，使得Kruppa方程的未知数由5个较少至3个，Kruppa方程的系数矩阵也由3×6变为3×4，降低了运算量。这一过程较为重要的技术环节还有：Kruppa方程的比例因子初始值的确定和利用遗传算法优化比例因子，将在稍后具体介绍。

对摄像机其他三个内参数的标定的具体流程如图4所示，包括如下步骤：

步骤2.1：保持摄像机内参数不变，通过调整摄像机的位姿状态，从不同角度对同一场景拍摄图像，获得至少三幅图像；

步骤2.2：利用SURF算法对步骤2.1获得的图像进行特征提取，并对两两图像中的特征点进行初始匹配操作；

步骤2.3：采用两视图极几何理论中对极几何约束的RANSAC算法对步骤2.2中的特征点对集合进行误匹配剔除操作，并将正确的匹配点对保存起来；

步骤2.4：判断步骤2.3中所获得的两两图像匹配点对数p是否大于第二阈值T2，若大于T2，则进行步骤2.5操作，否则舍掉该图像对，并返回步骤2.1，重复执行2.1至步骤2.4得到满足条件的图像对，T2≥8，优选地，可以取T2＝10；

步骤2.5：利用步骤2.3中所获得的两两图像匹配点对计算该图像对间的基本矩阵和外极点的反对称矩阵，然后构造该图像对的Kruppa方程，进而构造所有图像对的Kruppa方程组；

步骤2.6：利用公式求出所有图像对Kruppa方程的比例因子初始值，其中和是图像对间基本矩阵F的两个特征值，是外极点e′的反对称矩阵的特征值；

步骤2.7：利用遗传算法对所有图像对的Kruppa方程比例因子进行优化；

步骤2.8：将优化后的比例因子以及第1阶段获得的摄像机主点坐标代入步骤2.5构造的Kruppa方程组，然后利用奇异值分解法以及cholesky分解，最终求出摄像机内其他三个内参数。

在获取图像的过程中，要保持摄像机的位姿状态不同，且获得的图像要有足够的特征可供提取和匹配。优选的，在使用对极几何约束RANSAC算法剔除误匹配时，设定t＝0.005，遗传算法优化Kruppa方程的过程中，设定种群数为100，生成代数为200，共运行10次，设定σ_init＝0.3，θ＝0.2。

下面，具体对上述步骤中涉及的其它技术予以详细的说明。

1、利用SURF算法提取图像序列特征点

由于基于Kruppa方程的自标定方法需要计算两图像间的基本矩阵，故需要对两图像进行特征点提取和匹配，故本发明选取SURF算法(Bay H，Tuvtellars T，Gool L Van.SURF：Speeded Up Robust Features.LectureNotes in Computer Science，2006，3951：404-417)对图像进行特征提取。SURF算法是Herbert Bay首次提出的一种特征提取算法，全称叫“SpeededUp Robust Features”。该算法是对SIFT算法(Scale Invariant FeatureTransform)(David G L.Object Recognition from Local Scale-invariantFeatures.International Conference on Computer Vision，1999，1150-1157)的改进，不仅满足尺度不变性和旋转不变性，而且与SIFT算法相比，大大提高了计算速度。

2、基本矩阵和外极点的计算

基本矩阵描述了同一摄像机在不同角度下拍摄的两幅图像之间的极几何约束关系。图5示出了两视图极几何成像模型，其中O，O′分别是摄像机的光心在两个视点下的位置，M是空间中某一点，m，m′分别是M在两视图的像点，是一对对应点，记为连线OO′分别交两视图像平面于e，e′两点，这两点称为两视图下的极点，点m和e之间的连线为l，点m′和点e′之间的连线为l′，交线l、l′为一对极线对应，记为

基本矩阵从代数的角度上描述了极几何约束关系，它是不可逆的，秩为2，自由度为7的3×3矩阵。基本矩阵可通过图像间对应的点求出，具体关系如下：

m′^TFm＝0

此外基本矩阵与两图像的极点还满足：Fe＝0，F^Te′＝0。计算基本矩阵的方法有很多，如线性方法中的7点算法、8点算法，基于Newton-Raphson或Levenberg-Marquardt的迭代方法，基于M-estimators、LMedS(least-median-of-square)、RANSAC等算法的鲁棒估计方法。本发明采用较为简单的8点算法，具体过程如下：

给定n对图像点对应和它们应满足m_i′^TFm_i＝0，令基本矩阵 $F=\left(\begin{array}{ccc}{f}_{11}& f_{12}& f_{13}\\ f_{21}& f_{22}& f_{23}\\ f_{31}& f_{32}& f_{33}\end{array}\right)$ f＝(f₁₁，f₁₂，f₁₃，f₂₁，f₂₂，f₂₃，f₃₁，f₃₂，f₃₃)则展开m_i′^TFm_i＝0有：

$\mathrm{Mf}=\left(\begin{array}{ccccccccc}{u}_1^{\′}{u}_1& u_1^{\′}{v}_1& u_1^{\′}& v_1^{\′}{u}_1& v_1^{\′}{v}_1& v_1^{\′}& u_1& v_1& 1\\ \·& \·& \·& \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·& \·& \·& \·\\ u_n^{\′}{u}_n& u_n^{\′}{v}_n& u_n^{\′}& v_n^{\′}{u}_n& v_n^{\′}{v}_n& v_n^{\′}& u_n& v_n& f\end{array}\right)f=0$

当n≥8时就可以线性求解上述方程。对M进行奇异值分解，M＝UDV^T，则f为矩阵M的最小奇异值所对应的特征向量V₉，即矩阵V的最后一列。求出f后，构建矩阵F，然后对F做奇异值分解，F＝Udiag(t₁，t₂，t₃)V^T，则 就是最终基本矩阵的估计。

3、利用对极几何约束的RANSAC算法剔除误匹配点

利用前面的方法对特征点进行匹配往往存在明显的误匹配点对，如果误匹配点对参与摄像机的标定过程，则会严重影响参数计算精度，因此要消除误匹配点对。本发明使用对极几何约束的RANSAC算法剔除误匹配点对，具体流程图如图6所示，包括：

步骤1：给定两幅图像I₁、I₂，将通过前面的匹配方法获得的图像点对集合进行坐标归一化处理，并初始化循环次数i＝0；

步骤2：从点对集合中随机取8对点，然后利用8点算法计算基本矩阵F_i；

步骤3：使用sampson距离计算公式 $d=\frac{m^{\′T}{F}_{i}m}{{\left(F_{i}m\right)}_1^2+{\left(F_{i}m\right)}_2^2+{\left(F_i^T{m}^{\′}\right)}_1^2+{\left(F_i^T{m}^{\′}\right)}_2^2},$ 在匹配点对集合中寻找所有符合d＜t条件的点对，将它们做为内点，记录每一次的内点数量Q_i，其中t一般取(0.001，0.01)之间的数；需要说明的是，在理论上，匹配点对之间的sampson距离d应该为0，此处为避免浮点运算的不精确，选取一接近于0的小数t，与距离d进行比较，如果d＜t，则认为d＝0，选取该匹配点对作为内点。

步骤4：i＝i+1，若i＜N，则重复执行步骤2和步骤3，并记录每一次的内点数量，否则转入步骤5；其中N可通过下式计算：其中P代表所期望达到的概率，如设定P＝0.99，e是正确的数据比例，每次计算时不断更新；

步骤5：从步骤3求得的Q₁，Q₂，……，Q_N中选取最大的内点数Q_i，将其对应的基本矩阵F_i作为最终要求解的基本矩阵F，再以sampson距离为依据，寻找匹配点对集合中所有满足d＜t的点对，将它们做为最终的正确的匹配点对，其它的点对予以剔除。

4、Kruppa方程比例因子初始值的确定

给定一对图像，通过SURF特征提取、匹配及对极几何约束的RANSAC算法剔除误匹配得到两图像间的基本矩阵F，可建立关于摄像机内参数的约束方程，如下式：

${\mathrm{FCF}}^T=\mathrm{\μ}{\left(e^{\′}\right)}_{\×}C{\left(e^{\′}\right)}_{\×}^T$

其中，C＝KK^T，F为两图像间的基本矩阵，(e′)_×是第二幅图像外极点的反对称矩阵，则若获得n对图像，可得到如下的n个Kruppa方程：

$\left\{\begin{array}{c}{F}_1{\mathrm{CF}}_1^T={\mathrm{\μ}}_1{\left[e_1^{\′}\right]}_{\×}C{\left[e_1^{\′}\right]}_{\×}^T\\ F_2{\mathrm{CF}}_2^T={\mathrm{\μ}}_2{\left[e_2^{\′}\right]}_{\×}C{\left[e_2^{\′}\right]}_{\×}^T\\ \·\·\·\·\·\·\\ F_n{\mathrm{CF}}_n^T={\mathrm{\μ}}_n{\left[e_n^{\′}\right]}_{\×}C{\left[e_n^{\′}\right]}_{\×}^T\end{array}\right.$

根据上述方程组，选取每个Kruppa方程提供的三个约束，可重写上述方程组为如下形式：

A_3n×6(μ)c＝0

其中，c＝(c₁，c₂，c₃，c₄，c₅，1)^T，μ＝(μ₁，μ₂，......，μ_n)^T。

当摄像机图像主点坐标已知时， $C=K{K}^T=\left(\begin{array}{ccc}{f}_u^2+{\mathrm{\γ}}^2+u_0^2& {\mathrm{sf}}_v+u_0{v}_0& u_0\\ {\mathrm{sf}}_v+u_0{v}_0& f_v^2+v_0^2& v_0\\ u_0& v_0& 1\end{array}\right),$ 令c₂＝γf_v，C可写为 $C=\left(\begin{array}{ccc}{c}_1+u_0^2& c_2+u_0{v}_0& u_0\\ c_2+u_0{v}_0& c_3+v_0^2& v_0\\ u_0& v_0& 1\end{array}\right),$ 假设c＝(c₁，c₂，c₃，1)^T，则Kruppa方程组由A_3n×6(μ)c＝0变为A_3n×4(μ)c＝0，Kruppa方程方程组的系数矩阵由A_3n×6(μ)变为A_3n×4(μ)，系数矩阵的规模有所减小。

其中，A_3n×4(μ)是关于μ的函数，如下所示：

$A_{3n\×4}\left(\mathrm{\μ}\right)=\left[\begin{array}{c}{A}_{3\×3}^1\left({\mathrm{\μ}}_1\right),A_{3\×1}^1\left({\mathrm{\μ}}_1\right)\\ A_{3\×3}^2\left({\mathrm{\μ}}_2\right),A_{3\×1}^2\left({\mathrm{\μ}}_2\right)\\ \·\·\·\·\·\·\\ A_{3\×3}^n\left({\mathrm{\μ}}_n\right),A_{3\×1}^n\left({\mathrm{\μ}}_n\right)\end{array}\right]$

其中

$A_{3\×3}^k\left({\mathrm{\μ}}_k\right)=\left(\begin{array}{ccc}{\left(f_{11}^k\right)}^2& 2f_{11}^k{f}_{12}^k& {\left(f_{12}^k\right)}^2-{\mathrm{\μ}}_k\\ f_{11}^k{f}_{21}^k& f_{11}^k{f}_{22}^k+f_{12}^k{f}_{21}^k+{\mathrm{\μ}}_k& f_{22}^k{f}_{12}^k\\ f_{11}^k{f}_{31}^k& f_{12}^k{f}_{31}^k+f_{11}^k{f}_{32}^k-{\mathrm{\μ}}_k{e}_2^k& f_{12}^k{f}_{32}^k+{\mathrm{\μ}}_k{e}_1^k\end{array}\right)$

$A_{3\×1}^k\left({\mathrm{\μ}}_k\right)=\left(\begin{array}{c}{\left(f_{11}^k\right)}^2{u}_0^2+2f_{11}^k{f}_{12}^k{u}_0{v}_0+(2{\mathrm{\μ}}_k-{\mathrm{\μ}}_k)v_0^2+2f_{11}^k{f}_{13}^k{\mathrm{\μ}}_0+(2f_{12}^k{f}_{13}^k+{\left(f_{12}^k\right)}^2+2{\mathrm{\μ}}_k{e}_2^k)v_0+{\left(f_{13}^k\right)}^2-{\mathrm{\μ}}_k{e}_2^2\\ f_{11}^k{f}_{21}^k{u}_0^2+(f_{11}^k{f}_{22}^k+f_{21}^k{f}_{12}^k+{\mathrm{\μ}}_k)u_0{v}_0+f_{22}^k{f}_{12}^k{v}_0^2+(f_{23}^k{f}_{11}^k+f_{21}^k{f}_{13}^k-{\mathrm{\μ}}_k{e}_2^k)u_0+(f_{23}^k{f}_{12}^k+f_{22}^k{f}_{13}^k-{\mathrm{\μ}}_k{e}_1^k)v_0+f_{23}^k{f}_{13}^k+{\mathrm{\μ}}_k{e}_1^k{e}_2^k\\ f_{11}^k{f}_{31}^k{u}_0^2+(f_{12}^k{f}_{31}^k+f_{11}^k{f}_{32}^k-{\mathrm{\μ}}_k{e}_2^k)u_0{v}_0+(f_{12}^k{f}_{32}^k+{\mathrm{\μ}}_k{e}_1^k)v_0^2+(f_{13}^k{f}_{31}^k+f_{11}^k{f}_{33}^k+{\mathrm{\μ}}_k{\left(e_2^k\right)}^2)u_0+(f_{13}^k{f}_{32}^k+f_{12}^k{f}_{33}^k-{\mathrm{\μ}}_k{e}_1^k{e}_2^k)v_0+f_{13}^k{f}_{33}^k\end{array}\right).$

此时若系数矩阵A_3n×4(μ)已知时，就可以线性求解方程组A_3n×4(μ)c＝0。

由于系数矩阵A_3n×4(μ)含有未知的比例因子μ，故本发明采用遗传算法优化求解最佳的比例因子μ，此时需要给定比例因子μ的初始值，本发明采用下式确定所给n个Kruppa方程的比例因子初始值：

${\mathrm{\μ}}_k\≈\frac{{\mathrm{\λ}}_1^{F_k}{\mathrm{\λ}}_2^{F_k}}{{\left({\mathrm{\λ}}_{{\left[e_k^{\′}\right]}_{\×}}\right)}^2}$

其中，和是第K对图像基本矩阵F_k的两个特征值，是第K对图像外极点e′_k对应的反对称矩阵的特征值。

5、利用遗传算法优化比例因子

通常，遗传算法比传统的优化算法具有更好的全局搜索能力，故本发明采用遗传算法对Kruppa方程的比例因子μ_k进行优化。实现过程如下：

(1)种群初始化：种群中每个遗传个体为n个比例因子的组合，令为种群第1代第i个个体，则其中为第j个方程的比例因子，且其值为(1-σ_init，1+σ_init)μ_k的一个随机值，σ_init的值可根据结果进行调整。

(2)个体评价：由于Kruppa方程的系数矩阵秩为3，故可选取使系数矩阵A_3n×4(μ)第4个奇异值σ₄最小的个体为最佳个体，个体的适应度函数采用如下公式计算：其中为第j代个体中第i个比例因子。

(3)选择运算：根据矩阵C的正定性，可将那些根据比例因子求解出的矩阵C不满足正定性的个体适应度赋给一个很小的权值，使其选择的概率降低。然后通过“赌轮”策略，使那些适应度较高的个体有机会将基因传到下一代去。

(4)交叉运算：这里采用“单点”交叉策略，交叉点位置随机产生。

(5)变异运算：在变异运算中，将个体中需要变异的比例因子乘以一个在(1-θ，1+θ)之间的随机数。变异的概率需根据遗传个体适应度的大小进行自适应调节。

(6)终止条件：在预定的停止准则满足时终止算法优化过程，并选取所有代中适应度值最大的个体，得到优化后的比例因子组合(μ₁，μ₂，......μ_n)^T。

7、线性求解Kruppa方程

在使用遗传算法优化后，得到Kruppa方程比例因子的最终结果，将比例因子的最终结果代入Kruppa方程中，此时Kruppa方程的系数矩阵A_3n×4(μ)变为已知，则可对其进行奇异值分解，即A＝UDV^T，对矩阵V的最后一列元素做齐次化处理即为c＝(c₁，c₂，c₃，1)^T，求解出c后即可得到矩阵 $C=\left(\begin{array}{ccc}{c}_1+u_0^2& c_2+u_0{v}_0& u_0\\ c_2+u_0{v}_0& c_3+v_0^2& v_0\\ u_0& v_0& 1\end{array}\right),$ 再利用cholesky分解C＝KK^T，其中K为上三角阵，即为摄像机的内参数矩阵。

在具体分析完各个步骤之后，以仿真实验来证明本发明所提出的方法的有效性。在仿真实验中，通过在虚拟的三维空间中建立世界坐标系，然后在一定区域内随进生成一定数量的点，接着设定一个虚拟摄像机，根据随机生成的旋转向量和平移向量，生成空间点相应的图像，最后根据本发明求解该虚拟摄像机矩阵。

设定该虚拟摄像机的内参数矩阵为 $\left(\begin{array}{ccc}1000& 0& 500\\ 0& 1000& 500\\ 0& 0& 1\end{array}\right),$ 在中心坐标为(0，0，0)，长宽高均为200的正方体区域内随机产生100个点，并随机产生100个旋转向量和平移向量，且能使产生的100个点都能投影到摄像机图像平面1000×1000内，表1示出了本发明的实验结果。

表1本文方法的实验结果

其中，在标定主点的过程中需要固定摄像机位姿，然后调节摄像机焦距如下：

$\left(\begin{array}{ccc}1500& 0& 500\\ 0& 1500& 500\\ 0& 0& 1\end{array}\right),$

然后利用最小二乘法求得的。从表1的结果可以看出，本发明的标定结果与理论值接近，表明了本发明在摄像机内参数标定方面的可靠性。

表2是本文方法的鲁棒性检验，对所得到的100个点坐标增加均值为0，方差为1的高斯噪声扰动，试验次数为8。

表2鲁棒性检验

从表2可以看出，在对图像点增加一定的高斯噪声时，本发明的标定结果并没有出现较大的浮动，因此具有一定的鲁棒性。

本发明提出了一种分步的基于Kruppa方程的自标定方法。通过与传统的经典标定方法实验对比，该方法的标定结果精度高，且本发明在图像中存在一定噪声的情况下，鲁棒性较好，误差在允许的范围内。本发明中采用的方法和技术均可通过软件实现且运行条件容易满足，易于在各种增强现实及三维重建系统中使用。

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

Claims (5)

1.一种摄像机分步自标定方法，所述摄像机采取的成像矩阵为： $K=\left(\begin{array}{ccc}{f}_u& s& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right),$ 其中，(u₀，v₀)是摄像机主点坐标，s是倾斜因子，f_u、f_v分别为摄像机在成像平面u轴、v轴方向上的尺度因子；其特征在于，所述方法分为2个阶段进行：

第1阶段，标定摄像机的主点坐标(u₀，v₀)；

第2阶段，结合第1阶段获得的主点坐标，标定摄像机其他三个参数，即倾斜因子s以及摄像机在成像平面u轴、v轴方向上的尺度因子f_u、f_v；

其中在第1阶段中，通过调节摄像机的焦距，同时保持摄像机的位姿状态不变，对同一场景拍摄两张在不同焦距下的图像，对图像进行特征提取和匹配，并剔除误匹配点，最后利用匹配点对求解摄像机的主点坐标；在第2阶段中，控制摄像机的内参数不变，从不同角度对同一场景进行拍摄获得三幅的可供匹配的图像，对图像进行特征提取和匹配，并剔除误匹配点，然后基于Kruppa方程，确定Kruppa方程的比例因子初始值并利用遗传算法对其进行优化，通过优化后的比例因子以及第1阶段获得的主点坐标代入Kruppa方程实现摄像机其他三个参数的求解；其中所述第1阶段的具体流程包括：

步骤1.1：保持摄像机的位姿状态不变，调整摄像机的焦距，获得摄像机对同一场景分别在焦距f₁和f₂条件下的两张图像；

步骤1.2：利用SURF算法对获取的两张图像进行特征提取，然后对两图像所含的特征点进行初匹配操作，得到两图像特征点的匹配对集合；

步骤1.3：采用对极几何约束的RANSAC算法对步骤1.2中得到的匹配对集合进行误匹配剔除操作，将误匹配点对从匹配点对集合中去除，然后将剩余的匹配点对集合保存起来；

步骤1.4：判断步骤1.3中所获得的匹配点对数P是否大于第一阈值T1，若大于T1，则进行步骤1.5操作，否则转入步骤1.1，其中T1≥8；

步骤1.5：将步骤1.4所获得的P个匹配特征点对写成如下的方程组形式：

$\left\{\begin{array}{c}({y_1}^{\′}-y_1)u_0+(x_1-{x_1}^{\′})v_0=x_1{y_1}^{\′}-{x_1}^{\′}{y}_1\\ ({y_2}^{\′}-y_0)u_0+(x_2-{x_2}^{\′})v_0=x_2{y_2}^{\′}-{x_2}^{\′}{y}_2\\ .......\\ ({y_P}^{\′}-y_P)u_0+(x_P-{x_P}^{\′})v_0=x_P{y_P}^{\′}-{x_P}^{\′}{y}_P\end{array}\right.,$

利用最小二乘法求解该方程组，得到的解(u₀，v₀)便是摄像机的主点坐标，其中，(x_i，y_i)和(x′_i，y′_i)是两图像的匹配特征点对集合中第i对特征点坐标，i＝1，2，…，P。

2.根据权利要求1所述的方法，所述第2阶段的具体流程包括：

步骤2.1：保持摄像机内参数不变，通过调整摄像机的位姿状态，从不同角度对同一场景拍摄图像，获得至少三幅图像；

步骤2.2：利用SURF算法对步骤2.1获得的图像进行特征提取，并对两两图像中的特征点进行初始匹配操作；

步骤2.3：采用对极几何约束的RANSAC算法对步骤2.2中的得到的匹配对集合进行误匹配剔除操作，将误匹配点对从匹配点对集合中去除，然后将剩余的匹配点对集合保存起来；

步骤2.4：判断步骤2.3中所获得的两两图像匹配点对数p是否大于第二阈值T2，若大于T2，则进行步骤2.5操作，否则舍掉该图像对，并返回步骤2.1，重复执行步骤2.1到2.4，直到得到满足条件的图像对，T2≥8；

步骤2.5：利用步骤2.3中所获得的两两图像匹配点对计算该图像对间的基本矩阵和外极点的反对称矩阵，然后构造该图像对的Kruppa方程，进而构造所有图像对的Kruppa方程组；

步骤2.6：利用公式求出所有图像对Kruppa方程的比例因子初始值，其中是图像对间基本矩阵F的两个特征值，是外极点e′的反对称矩阵的特征值；

步骤2.7：利用遗传算法对所有图像对的Kruppa方程比例因子进行优化；

步骤2.8：将优化后的比例因子以及第1阶段获得的摄像机主点坐标代入步骤2.5构造的Kruppa方程组，然后利用奇异值分解法以及cholesky分解，最终求出摄像机内其他三个内参数。

3.根据权利要求1所述的方法，所述T1＝10。

4.根据权利要求2所述的方法，所述T2＝10。

5.根据前述任一权利要求所述的方法，所述对极几何约束的RANSAC算法具体包括：

步骤1：给定两幅图像I₁、I₂，将通过所述匹配方法获得的图像点对集合进行坐标归一化处理，并初始化循环次数a＝0；

步骤2：从点对集合中随机取8对点，然后利用8点算法计算基本矩阵F_i；

步骤3：使用sampson距离计算公式在匹配点对集合中寻找所有符合d＜t条件的点对，将它们做为内点，记录每一次的内点数量Q_i，其中t取(0.001，0.01)之间的数；

步骤4：a＝a+1，若a＜N，则重复执行步骤2和步骤3，并记录每一次的内点数量，否则转入步骤5；其中G代表所期望达到的概率，e是正确的数据比例，每次计算时不断更新；

步骤5：从步骤3求得的Q₁，Q₂，……，Q_N中选取最大的内点数Q_i，将其对应的基本矩阵F_i作为最终要求解的基本矩阵F，再以sampson距离为依据，寻找匹配点对集合中所有满足d＜t的点对，将它们做为最终的正确的匹配点对，其它的点对予以剔除。