CN116071433A - 相机标定方法和系统、计算机可读存储介质 - Google Patents

Abstract

一种相机标定方法和系统、计算机可读存储介质，该方法包括：获取相机拍摄的标定板图像及标定板图像中所有像素点对应的空间点的世界坐标；计算单应性矩阵；根据单应性矩阵与相机内参之间的约束关系计算相机的内参；根据所得内参计算像素点对应的空间点的相机坐标；判断所有空间点是否共面，若共面则移动世界坐标系的原点至空间点的重心，根据移动后世界坐标系相对于相机坐标系的变换关系计算相机的外参；若不共面，则根据空间点的相机坐标和世界坐标的对应关系求解相机的外参；对相机的内参、外参和畸变系数进行全局优化。该方法可适应像素点对应的空间点共面和不共面的情况，能够提高相机标定的精度，获得更准确的内参、外参和畸变系数。

Description

相机标定方法和系统、计算机可读存储介质

技术领域

本发明涉及相机标定技术领域，具体涉及相机标定方法和系统、计算机可读存储介质。

背景技术

在图像测量过程以及机器视觉应用中，为了确定空间物体表面某点的三维几何位置，必须建立相机成像的几何模型，也就是确定空间物体表面某点的三维几何位置与其在图像上的对应点的对应关系，这样在获得相机所拍摄图像的图像坐标后，根据相机成像的几何模型就可以推知相应的三维空间坐标。几何模型中的参数就是相机的参数，而确定相机的参数的过程称之为相机标定。相机参数的标定是非常关键的环节，标定结果的精度及标定算法的稳定性直接影响相机工作产生结果的准确性。因此，做好相机标定是做好后续工作的前提，有必要对相机标定方法进行改进，以提高相机标定的精度，获得更准确的相机参数。相机标定常利用标定板来进行，标定板在机器视觉、图像测量、摄影测量、三维重建等方面有广泛应用，通过相机拍摄带有固定间距图案阵列的标定板的图像，经过标定算法的计算，可以得出相机成像的几何模型，从而得到高精度的测量和重建结果。目前，通常采用带有棋盘格图案或实心圆阵列图案的标定板进行相机标定，其中棋盘格标定板通过定位棋盘格角点来获取特征点，圆阵列标定板通过定位圆点中心来获取特征点，在确定特征点的坐标以及与世界坐标的对应关系后即可以进行后续的标定工作。

发明内容

本发明提供一种相机标定方法和系统、计算机可读存储介质，旨在提高相机标定的精度，获得更准确的相机内参、外参和畸变系数。

根据第一方面，一种实施例中提供一种相机标定方法，包括：

获取相机拍摄的标定板图像；

获取标定板图像中所有像素点对应的空间点的世界坐标，并获取标定板图像中的特征点，以及特征点的图像坐标；

根据特征点的图像坐标和对应的空间点的世界坐标计算单应性矩阵H_a；

根据单应性矩阵与相机内参之间的约束关系，利用单应性矩阵H_a计算相机的内参；

根据计算得到的内参计算所有像素点对应的空间点的相机坐标；

根据标定板图像中所有像素点对应的空间点的世界坐标的协方差，判断所述空间点是否共面，若共面，则移动世界坐标系的原点至所述空间点的重心，根据像素点对应的空间点的相机坐标和移动后的世界坐标，计算移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c，根据移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c计算相机的外参；若不共面，则根据像素点对应的空间点的相机坐标和世界坐标的对应关系求解相机的外参；

根据预设的目标函数，对相机的畸变系数和计算得到的内参、外参进行全局优化，得到相机最终的内参、外参和畸变系数。

一种实施例中，所述相机的内参包括等效焦距f_x和f_y；所述根据单应性矩阵与相机内参之间的约束关系，利用单应性矩阵计算相机的内参，包括：

根据以下公式由单应性矩阵H_a计算单应性矩阵H_b：

H₀＝h₀-h₆c_x，H₁＝h₁-h₇c_x，H₂＝h₂-h₈c_x，

H₃＝h₃-h₆c_y，H₄＝h₄-h₇c_y，H₅＝h₅-h₈c_y，

H₆＝h₆，H₇＝h₇，H₈＝h₈，

其中(c_x,c_y)为相机主光轴点的坐标，h₀、h₁、h₂、h₃、h₄、h₅、h₆、h₇和h₈为单应性矩阵H_a中的元素，

H₀、H₁、H₂、H₃、H₄、H₅、H₆、H₇和H₈为单应性矩阵H_b中的元素，

根据以下约束关系由单应性矩阵H_b计算相机的等效焦距f_x和f_y：

一种实施例中，所述根据标定板图像中所有像素点对应的空间点的世界坐标的协方差，判断所述空间点是否共面，包括：

对所有像素点对应的空间点的世界坐标的协方差矩阵进行奇异值分解，对分解得到的特征值按从大到小排列为λ₁、λ₂、λ₃；

若λ₃/λ₁小于设定阈值，则判断所述空间点共面，否则判断所述空间点不共面。

一种实施例中，所述相机的外参包括世界坐标系相对于相机坐标系的旋转矩阵R和平移向量t；所述移动世界坐标系的原点至所述空间点的重心，根据像素点对应的空间点的相机坐标和移动后的世界坐标，计算移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c，根据移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c计算相机的外参，包括：

根据对所有像素点对应的空间点的世界坐标的协方差矩阵进行奇异值分解得到的特征向量计算得到旋转矩阵R_v＝[v₁,v₂,v₁×v₂]，其中v₁和v₂分别为特征值λ₁和λ₂对应的特征向量；

根据以下公式计算得到平移向量t_c：

其中(x,y,z)为像素点对应的空间点的世界坐标，

和

为平移向量t_c中的元素，

移动世界坐标系的原点至所述空间点的重心，根据以下公式计算得到旋转矩阵R_c：

其中(x_c,y_c,z_c)为像素点对应的空间点的相机坐标，(x′,y′,z′)为移动世界坐标系的原点至所述空间点的重心后，像素点对应的空间点的世界坐标；

根据以下公式计算得到旋转矩阵R和平移向量t：

t＝H_bt_c，R＝H_bR_c。

一种实施例中，所述相机的外参包括世界坐标系相对于相机坐标系的旋转矩阵R和平移向量t；所述根据像素点对应的空间点的相机坐标和世界坐标的对应关系求解相机的外参，包括：

根据以下方程计算参数a₀、a₁、a₂、a₃、a₄、a₅、a₆、a₇、a₈、a₉、a₁₀和a₁₁：

其中λ为比例参数，(x_c,y_c,z_c)为像素点对应的空间点的相机坐标，(x,y,z)为像素点对应的空间点的世界坐标，旋转矩阵

平移向量

一种实施例中，所述根据预设的目标函数，对相机的畸变系数和计算得到的内参、外参进行全局优化，得到相机最终的内参、外参和畸变系数，包括：

预先设置畸变系数的初始值，并将计算得到的内参、外参作为内参、外参的初始值，根据以下目标函数迭代求解最优解，得到相机最终的内参、外参和畸变系数：

其中n_m为标定板图像中特征点的数量，n₀为相机拍摄的标定板图像的数量，q_jk为第k幅标定板图像中第j个特征点的图像坐标，p_j为特征点对应的空间点的世界坐标，e_k表示根据第k幅标定板图像得到的外参，f_k表示根据第k幅标定板图像得到的内参，K表示畸变系数，φ(p_j,e_k,f_k,K)表示将世界坐标p_j根据已得的内参、外参和畸变系数进行变换得到的图像坐标。

一种实施例中，根据公式Q_l+1＝Q_l+δ进行迭代求解最优解，其中Q_l表示第l次迭代时相机的内参、外参和畸变系数组成的向量，δ由公式Jδ＝ε确定，其中ε为当前次迭代时所有标定板图像的所有特征点对应的q_jk-φ(p_j,e_k,r_k,K)的值组成的向量，矩阵J具体为

其中J_e表示特征点的图像坐标对于外参的偏导数，具体为

其中(u,v)为特征点的图像坐标，r＝[r_x,r_y,r_z]为基于旋转矩阵R的等效旋转轴，t＝[t_x,t_y,t_z]为平移向量；

J_e,k表示第k幅标定板图像中特征点的图像坐标对于外参的偏导数；

J_k表示第k幅标定板图像中特征点的图像坐标对于内参和畸变系数的偏导数，具体为

其中f_x和f_y为等效焦距，(c_x,c_y)为相机主光轴点的坐标，k₀、k₁、k₂、k₃、k₄、k₅、k₆和k₇为畸变系数。

根据第二方面，一种实施例中提供一种相机标定系统，包括：

标定板图像获取模块，用于获取相机拍摄的标定板图像；

图像和世界坐标获取模块，用于获取标定板图像中所有像素点对应的空间点的世界坐标，并获取标定板图像中的特征点，以及特征点的图像坐标；

单应性矩阵计算模块，用于根据特征点的图像坐标和对应的空间点的世界坐标计算单应性矩阵H_a；

内参计算模块，用于根据单应性矩阵与相机内参之间的约束关系，利用单应性矩阵H_a计算相机的内参；

相机坐标获取模块，用于根据计算得到的内参计算所有像素点对应的空间点的相机坐标；

外参计算模块，用于根据标定板图像中所有像素点对应的空间点的世界坐标的协方差，判断所述空间点是否共面，若共面，则移动世界坐标系的原点至所述空间点的重心，根据像素点对应的空间点的相机坐标和移动后的世界坐标，计算移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c，根据移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c计算相机的外参；若不共面，则根据像素点对应的空间点的相机坐标和世界坐标的对应关系求解相机的外参；

全局优化模块，用于根据预设的目标函数，对相机的畸变系数和计算得到的内参、外参进行全局优化，得到相机最终的内参、外参和畸变系数。

一种实施例中，所述相机的内参包括等效焦距f_x和f_y；所述内参计算模块具体用于：

根据以下公式由单应性矩阵H_a计算单应性矩阵H_b：

H₀＝h₀-h₆c_x，H₁＝h₁-h₇c_x，H₂＝h₂-h₈c_x，

H₃＝h₃-h₆c_y，H₄＝h₄-h₇c_y，H₅＝h₅-h₈c_y，

H₆＝h₆，H₇＝h₇，H₈＝h₈，

其中(c_x,c_y)为相机主光轴点的坐标，h₀、h₁、h₂、h₃、h₄、h₅、h₆、h₇和h₈为单应性矩阵H_a中的元素，

H₀、H₁、H₂、H₃、H₄、H₅、H₆、H₇和H₈为单应性矩阵H_b中的元素，

根据以下约束关系由单应性矩阵H_b计算相机的等效焦距f_x和f_y：

根据第三方面，一种实施例中提供一种计算机可读存储介质，所述介质上存储有程序，所述程序能够被处理器执行以实现上述第一方面所述的相机标定方法。

依据上述实施例的相机标定方法和系统、计算机可读存储介质，其中的方法首先获取相机拍摄的标定板图像，之后获取标定板图像中所有像素点对应的空间点的世界坐标，并获取标定板图像中的特征点，以及特征点的图像坐标，根据特征点的图像坐标和对应的空间点的世界坐标计算单应性矩阵，根据单应性矩阵与相机内参之间的约束关系，利用单应性矩阵计算相机的内参，再根据计算得到的内参计算所有像素点对应的空间点的相机坐标，根据标定板图像中所有像素点对应的空间点的世界坐标的协方差，判断空间点是否共面，若共面，则移动世界坐标系的原点至空间点的重心，根据移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c计算相机的外参；若不共面，则根据像素点对应的空间点的相机坐标和世界坐标的对应关系求解相机的外参，最后根据预设的目标函数，对相机的畸变系数和计算得到的内参、外参进行全局优化，得到相机最终的内参、外参和畸变系数。该方法可适应像素点对应的空间点共面和不共面的情况，对于共面和不共面的情况针对性地采用不同的方式求解外参，并且对计算得到的内参、外参和初始化的畸变系数根据预设的目标函数进行全局优化，构建全局最优的分析方法，从而能够提高相机标定的精度，获得更为准确的相机内参、外参和畸变系数。

附图说明

图1为针孔相机模型中各坐标系的变换示意图；

图2为一种实施例的相机标定方法的流程图；

图3为一种实施例的棋盘格标定板的示意图；

图4为一种实施例的圆阵列标定板的示意图；

图5为一种实施例的相机标定系统的结构示意图。

具体实施方式

下面通过具体实施方式结合附图对本发明作进一步详细说明。其中不同实施方式中类似元件采用了相关联的类似的元件标号。在以下的实施方式中，很多细节描述是为了使得本申请能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他元件、材料、方法所替代。在某些情况下，本申请相关的一些操作并没有在说明书中显示或者描述，这是为了避免本申请的核心部分被过多的描述所淹没，而对于本领域技术人员而言，详细描述这些相关操作并不是必要的，他们根据说明书中的描述以及本领域的一般技术知识即可完整了解相关操作。

另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书和附图中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。

本文中为部件所编序号本身，例如“第一”、“第二”等，仅用于区分所描述的对象，不具有任何顺序或技术含义。而本申请所说“连接”、“联接”，如无特别说明，均包括直接和间接连接(联接)。

在高精度测量场景中，由于镜头的设计或者安装的原因，会导致产生桶形畸变或者枕形畸变，为了提升测量的精度，这个时候需要完成畸变的校正。畸变是在由镜头到相机的像平面的变化过程中产生的，因此需要获取相机的内参信息，来对畸变进行一个校正运算和处理。相机镜头的畸变通常是一种非线性畸变，相机的线性模型并不能准确的描述相机的成像几何关系，需要对这种非线性畸变进行修正，目前常采用的方法有固定畸变模型(如切向畸变、径向畸变等)进行优化计算、无固定畸变模型智能算法等，这些算法在进行全局优化时，无法兼具计算精度与效率。

相机标定的目的就是为了获得相机的内参、外参和畸变系数。当前相机标定主要根据张正友标定法进行计算，张正友标定法主要包括以下几个计算步骤：

(1)根据标定板中特征点的图像坐标和对应的空间点的世界坐标，获得单应性矩阵；对于空间中某一点P_w，若其经相机成像后得到像素点q_i，则称空间点P_w为像素点q_i对应的空间点，像素点q_i为空间点P_w对应的像素点；

(2)对单应性矩阵进行分解，计算得到内参和/或外参的初始参数；

(3)采用LM(Levenberg-Marquardt，列文伯格-马夸尔特)算法对初始参数进行非线性优化，迭代计算内参、外参和畸变系数，获得最终的标定结果。

为了能更清楚地理解本发明的技术方案，下面对相机成像过程中的坐标变换关系进行介绍。如图1所示，相机的成像过程中各坐标系的投影变换关系可以用针孔相机模型来表示。世界坐标系(WCS，World Coordinate System)中的点(即空间点)P_w通过镜头投影中心投影到成像平面上的点P，为了得到点P_w投影到成像平面上的图像坐标q_i，需要首先将其转换到相机坐标系(CCS，Camera Coordinate System)中。相机坐标系的x轴和y轴分别平行于图像的u轴和v轴，z轴垂直于图像所在的成像平面并且z轴的方向设置将使相机前的所有点的z坐标为正数，其中图像的u轴方向为图像的水平方向，v轴方向为图像的垂直方向。图1中x_c轴、y_c轴和z_c轴分别表示相机坐标系的x轴、y轴和z轴。由世界坐标系到相机坐标系的变换可以用公式p_c＝^cH_wp_w来表示，其中p_c＝(x_c,y_c,z_c)^T为相机坐标系中的坐标，p_w＝(x_w,y_w,z_w)^T为世界坐标系中的坐标，^cH_w可以用旋转矩阵R和平移向量t来表示。

用矩阵形式表示为：

记

因此可得：

x_c＝R₀x_w+R₁y_w+R₂z_w+t₀

y_c＝R₃x_w+R₄y_w+R₅z_w+t₁。             (1)

z_c＝R₆x_w+R₇y_w+R₈z_w+t₂

在将世界坐标系转换到相机坐标系后，需要将其转换到像平面坐标系上，这是一个由3D坐标转换到2D坐标的过程。这一变换可以表示为：

用矩阵形式表示为：

其中f表示相机镜头的焦距，(x_I,y_I)^T表示像平面坐标系中的坐标。

结合公式(1)可得：

在投影到成像平面后，镜头的畸变将导致坐标q_c＝(x_I,y_I)^T发生变化，使得在成像平面上形成的是畸变的坐标p_d＝(x_d,y_d)^T，这个变化可以单独在成像平面上建立模型，也就是说这里不需要三维信息。畸变包括径向畸变和切向畸变，可以用以下模型描述畸变：

其中

k₀、k₁、k₂、k₃、k₄、k₅、k₆、k₇为畸变系数，包括了径向畸变的畸变系数和切向畸变的畸变系数。用矩阵形式表示为：

其中矩阵[K]表示按以上模型根据畸变系数k₀、k₁、k₂、k₃、k₄、k₅、k₆、k₇进行的畸变变换。

根据以上模型，在已知畸变系数和畸变后坐标p_d＝(x_d,y_d)^T的情况下，利用牛顿法可以求解得到未畸变的坐标q_c＝(x_I,y_I)^T。

最后将像平面坐标系转换到图像坐标系(ICS，Image Coordinate System)中，用公式表示为：

用矩阵形式表示为：

其中s_x和s_y分别为相机水平方向上和垂直方向上的像素尺寸，(c_x,c_y)为相机主光轴点，一般为图像的中心。

因此上述的整个变换如果不考虑畸变可以表示为：

这便是相机标定所依据的数学模型。其中

为相机的内参部分，旋转矩阵R和平移向量t为外参部分。

在上述模型的基础上，本发明提供一种相机标定方法，请参考图2，一种实施例中该方法包括步骤110～190，下面具体说明。

步骤110：获取相机拍摄的标定板图像。

标定板可以是棋盘格标定板(如图3所示)、圆阵列标定板(如图4所示)等。在进行标定时，可以根据经验将相机置于多个位姿(位姿也即相机相对于标定板的位置和角度)，并在处于每个位姿时对标定板进行拍摄，从而得到多个不同的标定板图像用于进行标定。

步骤120：获取标定板图像中所有像素点对应的空间点的世界坐标，并获取标定板图像中的特征点，以及特征点的图像坐标。

对于标定板图像来说，其是由标定板成像而得到，因此标定板图像中像素点对应的空间点为标定板上的点。世界坐标系为一空间坐标系，可以根据标定板的参数信息来构建世界坐标系从而得到空间点的世界坐标，标定板的参数信息可以包括标定板的尺寸、棋盘格的尺寸、圆形特征点的半径、特征点之间的间距等。对于棋盘格标定板来说特征点是棋盘格的角点，对于圆阵列标定板来说特征点是圆阵列中圆形特征点的重心，圆形特征点即圆阵列标定板上的圆形图案。关于像素点对应的空间点的世界坐标、标定板图像中的特征点以及特征点的图像坐标的获取，可以通过现有技术实现，在此不进行赘述。

步骤130：根据特征点的图像坐标和对应的空间点的世界坐标计算单应性矩阵H_a。可以理解，利用多个特征点的图像坐标和对应的空间点的世界坐标可以利用最小二乘法计算出单应性矩阵，记为H_a。

步骤140：根据单应性矩阵与相机内参之间的约束关系，利用单应性矩阵H_a计算相机的内参。

令等效焦距f_x＝f/s_x，f_y＝f/s_y，则内参部分可表示为

建立世界坐标系时，通常认为标定板上的点位于平面z＝0上，因此可以忽略z方向上的旋转和平移，因此外参部分可以表示为

其中r₁、r₂、r₃、r₄、r₅和r₆为旋转矩阵的元素，t_x和t_y分别为平移向量的x分量和y分量。

将单应性矩阵H_a中的元素用h₀、h₁、h₂、h₃、h₄、h₅、h₆、h₇和h₈表示，则

根据单应性矩阵和内外参的关系，可得：

因此可得h₀＝f_xr₁+c_xr₃，h₁＝f_xr₄+c_xr₆，h₂＝f_xt_x+c_x，h₃＝f_yr₂+c_yr₃，h₄＝f_yr₅+c_yr₆，h₅＝f_yt_y+c_y，h₆＝r₃，h₇＝r₆，h₈＝1。

若将图像坐标系原点设为图像中心，则所得单应性矩阵H_b为：

因此可得H₀＝f_xr₁，h₁＝f_xr₄，H₂＝f_xt_x，H₃＝f_yr₂，H₄＝f_yr₅，H₅＝f_yt_y，H₆＝r₃，H₇＝r₆，H₈＝1。因此可根据以下公式由单应性矩阵H_a计算单应性矩阵H_b：

H₀＝h₀-h₆c_x，H₁＝h₁-h₇c_x，H₂＝h₂-h₈c_x，

H₃＝h₃-h₆c_y，H₄＝h₄-h₇c_y，H₅＝h₅-h₈c_y，

H₆＝h₆，H₇＝h₇，H₈＝h₈，

其中相机主光轴点坐标(c_x,c_y)可以通过相机手册得知。此时

由旋转矩阵中各向量的单位正交性约束可得单应性矩阵H_b与相机内参之间的约束关系，具体地，根据单位正交性可建立方程：

将式(2)代入可得约束关系：

则根据以上约束关系可由单应性矩阵H_b计算得到等效焦距f_x和f_y。

步骤150：根据计算得到的内参计算所有像素点对应的空间点的相机坐标。

由上文对相机成像过程中坐标的变换关系可知，根据相机内参可由像素点的图像坐标得到对应的相机坐标系中的坐标，也就是像素点对应的空间点的相机坐标。

步骤160：根据标定板图像中所有像素点对应的空间点的世界坐标的协方差，判断空间点是否共面，若共面，则执行步骤170，若不共面则执行步骤180。

由于在建立世界坐标系时可能存在误差，导致所获得的空间点的世界坐标并不共面，因此本发明在计算外参前首先判断像素点对应的空间点是否共面，可适应像素点对应的空间点共面和不共面的情况，并且对于共面和不共面的情况针对性地采用不同的方式求解外参，有利于提高外参估计的准确性。

一种实施例中，可以根据所有像素点对应的空间点的世界坐标的协方差矩阵进行共面与否判断，协方差矩阵可以表示为：

其中(x_i,y_i,z_i)为第i个像素点对应的空间点的世界坐标，c_x、c_y、c_z分别表示所有像素点对应的空间点的世界坐标的x坐标、y坐标、z坐标的平均值，n表示像素点的数量。

判断空间点是否共面，具体地，可以对协方差矩阵进行奇异值分解(Single ValueDecomposition，SVD)，则可以得到三个特征值，按从大到小排列为λ₁、λ₂、λ₃，对应的特征向量记为v1、v2、v3，若λ₃/λ₁小于设定阈值，则判断空间点共面，否则判断空间点不共面。

步骤170：移动世界坐标系的原点至空间点的重心，根据像素点对应的空间点的相机坐标和移动后的世界坐标，计算移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c，根据移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c计算相机的外参。

一种实施例中，可以首先根据对上述协方差矩阵进行奇异值分解得到的特征向量计算得到一旋转矩阵R_v＝[v₁,v₂,v₁×v₂]，之后根据以下公式计算得到平移向量t_c：

其中(x,y,z)为像素点对应的空间点的世界坐标，

，即

和

为平移向量t_c中的元素，分别为平移向量t_c的x分量、y分量和z分量。

然后移动世界坐标系的原点至空间点的重心，根据以下公式计算得到旋转矩阵R_c：

其中(x_c,y_c,z_c)为像素点对应的空间点的相机坐标，(x′,y′,z′)为移动世界坐标系的原点至空间点的重心后，像素点对应的空间点的世界坐标。空间点的重心可以是上述坐标平均值c_x、c_y、c_z构成的坐标点(c_x,c_y,c_z)。

最后根据以下公式计算得到旋转矩阵R和平移向量t：

t＝H_bt_c，R＝H_bR_c。

步骤180：根据像素点对应的空间点的相机坐标和世界坐标的对应关系求解相机的外参。

对于一个像素点，已经从上述步骤获得其对应的空间点的相机坐标和世界坐标，可获得一个相机坐标和世界坐标的坐标对，利用多个这样的坐标对则可计算得到旋转矩阵R和平移向量t。

一种实施例中，可以根据以下方程计算参数a₀、a₁、a₂、a₃、a₄、a₅、a₆、a₇、a₈、a₉、a₁₀和a₁₁：

则旋转矩阵

平移向量

其中λ为比例参数，(x_c,y_c,z_c)和(x,y,z)分别为同一像素点对应的空间点的相机坐标和世界坐标。需要说明的是，比例参数λ并不一定需要预先设置或者求解出具体数值，只要能求解出参数a₀、a₁、a₂、a₃、a₄、a₅、a₆、a₇、a₈、a₉、a₁₀和a₁₁即可。

本发明提供一种求解上述方程的方法，具体地，对上式展开后，可得：

λx_c＝a₀x+a₁y+a₂z+a₃

λy_c＝a₄x+a₅y+a₆z+a₇，

λz_c＝a₈x+a₉y+a₁₀z+a₁₁

将第三行代入前两行可得：

a₀xz_c+a₁yz_c+a₂zz_c+a₃z_c-a₈xx_c-a₉yx_c-a₁₀zx_c-a₁₁x_c＝0

a₄xz_c+a₅yz_c+a₆zz_c+a₇z_c-a₈xy_c-a₉yy_c-a₁₀zy_c-a₁₁y_c＝0，

写成矩阵形式为：

其中X＝[a₀,a₁,a₂,a₃,a₄,a₅,a₆,a₇,a₈,a₉,a₁₀,a₁₁]^T，X左边的矩阵为系数矩阵，记为S，若根据N个像素点进行计算，则系数矩阵S为一2N×12的矩阵。将矩阵S^TS做SVD分解，则所得的最小特征值对应的特征向量为X的解。

步骤190：根据预设的目标函数，对相机的畸变系数和计算得到的内参、外参进行全局优化，得到相机最终的内参、外参和畸变系数。

本步骤主要对内参、外参和畸变系数进行全局优化，同时也是对畸变系数进行非线性优化，实现畸变校正。全局优化依据预设的目标函数，对内参、外参和畸变系数的初始值进行迭代求解最优解，当达到预设的停止条件时停止迭代，得到相机最终的内参、外参和畸变系数。内参、外参的初始值可以为上述步骤计算得到的内参、外参，而畸变系数的初始值可以根据经验预先设置。

一种实施例中，可以根据以下目标函数迭代求解最优解，得到相机最终的内参、外参和畸变系数：

其中n_m为标定板图像中特征点的数量，n₀为相机拍摄的标定板图像的数量，其中不同的标定板图像为相机和/或标定板处于不同位姿时所拍摄得到，q_jk为第k幅标定板图像中第j个特征点的图像坐标，p_j为特征点对应的空间点的世界坐标，e_k表示根据第k幅标定板图像得到的外参，f_k表示根据第k幅标定板图像得到的内参，K表示畸变系数，φ(p_j,e_k,f_k,K)表示将世界坐标p_j根据已得的内参、外参和畸变系数进行变换得到的图像坐标，由上文可知这其中包括利用外参将世界坐标变换到相机坐标，再利用内参和畸变系数将相机坐标变换到图像坐标的过程。

一种实施例中，可以利用LM算法进行迭代计算，迭代过程中参数的更新可以表示为Q_l+1＝Q_l+δ，其中Q_l表示第l次迭代时相机的内参、外参和畸变系数组成的向量，δ由公式Jδ＝ε确定，其中ε为当前次迭代时所有标定板图像的所有特征点对应的q_jk-φ(p_j,e_k,r_k,K)的值组成的向量，矩阵J具体可以为

其中J_e表示特征点的图像坐标对于外参的偏导数，具体为

其中(u,v)为特征点的图像坐标，r＝[r_x,r_y,r_z]为基于旋转矩阵R的等效旋转轴，t＝[t_x,t_y,t_z]为平移向量。旋转矩阵R可以用等效旋转轴r和等效轴角θ来表示，因此等效旋转轴r和等效轴角θ可视为外参，这里对于旋转矩阵R部分的偏导数使用对于等效旋转轴r的偏导数表示。

J_e,k则表示第k幅标定板图像中特征点的图像坐标对于外参的偏导数。

J_k表示第k幅标定板图像中特征点的图像坐标对于内参和畸变系数的偏导数，具体为

由式(3)可以看到，本发明更新相机参数的方式，所使用的矩阵J为一稀疏矩阵，从而提高了计算速度。

下面说明各部分偏导数的计算方法。

对于内参部分，特征点的图像坐标相对于等效焦距的偏导数为：

特征点的图像坐标相对于相机主光轴点的偏导数为：

特征点的图像坐标相对于畸变系数的偏导数可分为几部分计算，相对于k₀和k₁的偏导数为：

相对于k₂和k₃的偏导数为：

相对于k₄的偏导数为：

相对于k₅的偏导数为：

相对于k₆的偏导数为：

相对于k₇的偏导数为：

对于外参部分，特征点的图像坐标相对于平移向量的偏导数为：

其中

d₁＝1+k₀(x_c ²+y_c ²)+k₁(x_c ²+y_c ²)²+k₄(x_c ²+y_c ²)³，

令

可得：

则d₁对平移向量的偏导数为：

d₂对平移向量的偏导数为：

特征点的图像坐标相对于等效旋转轴的偏导数为：

其中

上述部分计算过程用链式法则可以表示为：

其中m＝(u,v)，c＝(c_x,c_y)，f＝(f_x,f_y)，K表示畸变系数。

至此则完成了对矩阵J的求解。

在上述相机标定方法的基础上，本发明还提供一种相机标定系统，请参考图5，一种实施例中该系统包括标定板图像获取模块1、图像和世界坐标获取模块2、单应性矩阵计算模块3、内参计算模块4、相机坐标获取模块5、外参计算模块6和全局优化模块7，下面分别说明。

标定板图像获取模块1用于获取相机拍摄的标定板图像。

图像和世界坐标获取模块2用于获取标定板图像中所有像素点对应的空间点的世界坐标，并获取标定板图像中的特征点，以及特征点的图像坐标。关于像素点对应的空间点的世界坐标、标定板图像中的特征点以及特征点的图像坐标的获取，可以通过现有技术实现，在此不进行赘述。

单应性矩阵计算模块3用于根据特征点的图像坐标和对应的空间点的世界坐标计算单应性矩阵H_a。

内参计算模块4用于根据单应性矩阵与相机内参之间的约束关系，利用单应性矩阵H_a计算相机的内参。一种实施例中，内参计算模块4具体用于：

根据以下公式由单应性矩阵H_a计算单应性矩阵H_b：

H₀＝h₀-h₆c_x，H₁＝h₁-h₇c_x，H₂＝h₂-h₈c_x，

H₃＝h₃-h₆c_y，H₄＝h₄-h₇c_y，H₅＝h₅-h₈c_y，

H₆＝h₆，H₇＝h₇，H₈＝h₈，

再根据以下约束关系由单应性矩阵H_b计算相机的等效焦距f_x和f_y：

相机坐标获取模块5用于根据计算得到的内参计算所有像素点对应的空间点的相机坐标。

外参计算模块6用于根据标定板图像中所有像素点对应的空间点的世界坐标的协方差，判断空间点是否共面，若共面，则移动世界坐标系的原点至空间点的重心，根据像素点对应的空间点的相机坐标和移动后的世界坐标，计算移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c，根据移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c计算相机的外参；若不共面，则根据像素点对应的空间点的相机坐标和世界坐标的对应关系求解相机的外参。

一种实施例中，外参计算模块6具体用于对所有像素点对应的空间点的世界坐标的协方差矩阵进行奇异值分解，对分解得到的特征值按从大到小排列为λ₁、λ₂、λ₃，若λ₃/λ₁小于设定阈值，则判断所述空间点共面，否则判断所述空间点不共面。

一种实施例中，若标定板图像中所有像素点对应的空间点共面，则外参计算模块6具体用于：根据对所有像素点对应的空间点的世界坐标的协方差矩阵进行奇异值分解得到的特征向量计算得到旋转矩阵R_v＝[v₁,v₂,v₁×v₂]，根据以下公式计算得到平移向量t_c：

移动世界坐标系的原点至空间点的重心，根据以下公式计算得到旋转矩阵R_c：

根据以下公式计算得到旋转矩阵R和平移向量t：

t＝H_bt_c，R＝H_bR_c。

一种实施例中，若标定板图像中所有像素点对应的空间点不共面，则外参计算模块6具体用于：根据以下方程计算参数a₀、a₁、a₂、a₃、a₄、a₅、a₆、a₇、a₈、a₉、a₁₀和a₁₁：

则旋转矩阵

平移向量

全局优化模块7用于根据预设的目标函数，对相机的畸变系数和计算得到的内参、外参进行全局优化，得到相机最终的内参、外参和畸变系数。

一种实施例中，全局优化模块7具体用于：

预先设置畸变系数的初始值，并将计算得到的内参、外参作为内参、外参的初始值，根据以下目标函数迭代求解最优解，得到相机最终的内参、外参和畸变系数：

其中n_m为标定板图像中特征点的数量，n₀为相机拍摄的标定板图像的数量，其中不同的标定板图像为相机和/或标定板处于不同位姿时所拍摄得到，q_jk为第k幅标定板图像中第j个特征点的图像坐标，p_j为特征点对应的空间点的世界坐标，e_k表示根据第k幅标定板图像得到的外参，f_k表示根据第k幅标定板图像得到的内参，K表示畸变系数，φ(p_j,e_k,f_k,K)表示将世界坐标pj根据已得的内参、外参和畸变系数进行变换得到的图像坐标，由上文可知这其中包括利用外参将世界坐标变换到相机坐标，再利用内参和畸变系数将相机坐标变换到图像坐标的过程。

一种实施例中，全局优化模块7还用于：根据公式Q_l+1＝Q_l+δ进行迭代求解最优解，其中Q_l表示第l次迭代时相机的内参、外参和畸变系数组成的向量，δ由公式Jδ＝ε确定，其中ε为当前次迭代时所有标定板图像的所有特征点对应的q_jk-φ(p_j,e_k,r_k,K)的值组成的向量，矩阵J具体为

其中J_e表示特征点的图像坐标对于外参的偏导数，具体为

其中(u,v)为特征点的图像坐标，r＝[r_x,r_y,r_z]为基于旋转矩阵R的等效旋转轴，t＝[t_x,t_y,t_z]为平移向量。

则J_e,k表示第k幅标定板图像中特征点的图像坐标对于外参的偏导数。

J_k表示第k幅标定板图像中特征点的图像坐标对于内参和畸变系数的偏导数，具体为

关于各偏导数的计算方法可参考上文关于相机标定方法的相应内容，在此不再赘述。

依据上述实施例的相机标定方法和系统，在进行内参估计时只需要一张标定板图像，直接根据该标定板图像即可计算出内参估计值，而现有技术需要多于一张的标定板图像，因此本发明的相机标定方法和系统低了对图像数量的要求，提高了效率。此外，本发明的相机标定方法和系统可适应像素点对应的空间点共面和不共面的情况，对于共面和不共面的情况针对性地采用不同的方式求解外参，并且对计算得到的内参、外参和初始化的畸变系数根据预设的目标函数进行全局优化，构建全局最优的分析方法，在进行优化时可以将径向畸变的畸变系数和切向畸变的畸变系数一起进行优化，不仅可以消除镜像畸变，也可以消除切向畸变。综上所述，本发明的技术方案能够提高相机标定的精度，获得更为准确的相机内参、外参和畸变系数。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。

以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。

Claims (10)

1.一种相机标定方法，其特征在于，包括：

获取相机拍摄的标定板图像；

获取标定板图像中所有像素点对应的空间点的世界坐标，并获取标定板图像中的特征点，以及特征点的图像坐标；

根据特征点的图像坐标和对应的空间点的世界坐标计算单应性矩阵H_a；

根据单应性矩阵与相机内参之间的约束关系，利用单应性矩阵H_a计算相机的内参；

根据计算得到的内参计算所有像素点对应的空间点的相机坐标；

根据标定板图像中所有像素点对应的空间点的世界坐标的协方差，判断所述空间点是否共面，若共面，则移动世界坐标系的原点至所述空间点的重心，根据像素点对应的空间点的相机坐标和移动后的世界坐标，计算移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c，根据移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c计算相机的外参；若不共面，则根据像素点对应的空间点的相机坐标和世界坐标的对应关系求解相机的外参；

根据预设的目标函数，对相机的畸变系数和计算得到的内参、外参进行全局优化，得到相机最终的内参、外参和畸变系数。

2.如权利要求1所述的相机标定方法，其特征在于，所述相机的内参包括等效焦距f_x和f_y；所述根据单应性矩阵与相机内参之间的约束关系，利用单应性矩阵计算相机的内参，包括：

根据以下公式由单应性矩阵H_a计算单应性矩阵H_b：

H₀＝h₀-h₆c_x，H₁＝h₁-h₇c_x，H₂＝h₂-h₈c_x，

H₃＝h₃-h₆c_y，H₄＝h₄-h₇c_y，H₅＝h₅-h₈c_y，

H₆＝h₆，H₇＝h₇，H₈＝h₈，

其中(c_x,c_y)为相机主光轴点的坐标，h₀、h₁、h₂、h₃、h₄、h₅、h₆、h₇和h₈为单应性矩阵H_a中的元素，

H₀、H₁、H₂、H₃、H₄、H₅、H₆、H₇和H₈为单应性矩阵H_b中的元素，

根据以下约束关系由单应性矩阵H_b计算相机的等效焦距f_x和f_y：

3.如权利要求1所述的相机标定方法，其特征在于，所述根据标定板图像中所有像素点对应的空间点的世界坐标的协方差，判断所述空间点是否共面，包括：

对所有像素点对应的空间点的世界坐标的协方差矩阵进行奇异值分解，对分解得到的特征值按从大到小排列为λ₁、λ₂、λ₃；

若λ₃/λ₁小于设定阈值，则判断所述空间点共面，否则判断所述空间点不共面。

4.如权利要求3所述的相机标定方法，其特征在于，所述相机的外参包括世界坐标系相对于相机坐标系的旋转矩阵R和平移向量t；所述移动世界坐标系的原点至所述空间点的重心，根据像素点对应的空间点的相机坐标和移动后的世界坐标，计算移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c，根据移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c计算相机的外参，包括：

根据对所有像素点对应的空间点的世界坐标的协方差矩阵进行奇异值分解得到的特征向量计算得到旋转矩阵R_v＝[v₁,v₂,v₁×v₂]，其中v₁和v₂分别为特征值λ₁和λ₂对应的特征向量；

根据以下公式计算得到平移向量t_c：

其中(x,y,z)为像素点对应的空间点的世界坐标，

和

为平移向量t_c中的元素，

移动世界坐标系的原点至所述空间点的重心，根据以下公式计算得到旋转矩阵R_c：

其中(x_c,y_c,z_c)为像素点对应的空间点的相机坐标，(x′,y′,z′)为移动世界坐标系的原点至所述空间点的重心后，像素点对应的空间点的世界坐标；

根据以下公式计算得到旋转矩阵R和平移向量t：

t＝H_bt_c，R＝H_bR_c。

5.如权利要求1所述的相机标定方法，其特征在于，所述相机的外参包括世界坐标系相对于相机坐标系的旋转矩阵R和平移向量t；所述根据像素点对应的空间点的相机坐标和世界坐标的对应关系求解相机的外参，包括：

根据以下方程计算参数a₀、a₁、a₂、a₃、a₄、a₅、a₆、a₇、a₈、a₉、a₁₀和a₁₁：

其中λ为比例参数，(x_c,y_c,z_c)为像素点对应的空间点的相机坐标，(x,y,z)为像素点对应的空间点的世界坐标，旋转矩阵

平移向量

6.如权利要求1至5中任一项所述的相机标定方法，其特征在于，所述根据预设的目标函数，对相机的畸变系数和计算得到的内参、外参进行全局优化，得到相机最终的内参、外参和畸变系数，包括：

预先设置畸变系数的初始值，并将计算得到的内参、外参作为内参、外参的初始值，根据以下目标函数迭代求解最优解，得到相机最终的内参、外参和畸变系数：

其中n_m为标定板图像中特征点的数量，n₀为相机拍摄的标定板图像的数量，q_jk为第k幅标定板图像中第j个特征点的图像坐标，p_j为特征点对应的空间点的世界坐标，e_k表示根据第k幅标定板图像得到的外参，f_k表示根据第k幅标定板图像得到的内参，K表示畸变系数，φ(p_j,e_k,f_k,K)表示将世界坐标p_j根据已得的内参、外参和畸变系数进行变换得到的图像坐标。

7.如权利要求6所述的相机标定方法，其特征在于，根据公式Q_l+1＝Q_l+δ进行迭代求解最优解，其中Q_l表示第l次迭代时相机的内参、外参和畸变系数组成的向量，δ由公式Jδ＝ε确定，其中ε为当前次迭代时所有标定板图像的所有特征点对应的q_jk-φ(p_j,e_k,r_k,K)的值组成的向量，矩阵J具体为

其中J_e表示特征点的图像坐标对于外参的偏导数，具体为

其中(u,v)为特征点的图像坐标，r＝[r_x,r_y,r_z]为基于旋转矩阵R的等效旋转轴，t＝[t_x,t_y,t_z]为平移向量；

J_e,k表示第k幅标定板图像中特征点的图像坐标对于外参的偏导数；

J_k表示第k幅标定板图像中特征点的图像坐标对于内参和畸变系数的偏导数，具体为

其中f_x和f_y为等效焦距，(c_x,c_y)为相机主光轴点的坐标，k₀、k₁、k₂、k₃、k₄、k₅、k₆和k₇为畸变系数。

8.一种相机标定系统，其特征在于，包括：

标定板图像获取模块，用于获取相机拍摄的标定板图像；

图像和世界坐标获取模块，用于获取标定板图像中所有像素点对应的空间点的世界坐标，并获取标定板图像中的特征点，以及特征点的图像坐标；

单应性矩阵计算模块，用于根据特征点的图像坐标和对应的空间点的世界坐标计算单应性矩阵H_a；

内参计算模块，用于根据单应性矩阵与相机内参之间的约束关系，利用单应性矩阵H_a计算相机的内参；

相机坐标获取模块，用于根据计算得到的内参计算所有像素点对应的空间点的相机坐标；

外参计算模块，用于根据标定板图像中所有像素点对应的空间点的世界坐标的协方差，判断所述空间点是否共面，若共面，则移动世界坐标系的原点至所述空间点的重心，根据像素点对应的空间点的相机坐标和移动后的世界坐标，计算移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c，根据移动后的世界坐标系相对于相机坐标系的旋转矩阵R_c和平移向量t_c计算相机的外参；若不共面，则根据像素点对应的空间点的相机坐标和世界坐标的对应关系求解相机的外参；

全局优化模块，用于根据预设的目标函数，对相机的畸变系数和计算得到的内参、外参进行全局优化，得到相机最终的内参、外参和畸变系数。

9.如权利要求8所述的相机标定系统，其特征在于，所述相机的内参包括等效焦距f_x和f_y；所述内参计算模块具体用于：

根据以下公式由单应性矩阵H_a计算单应性矩阵H_b：

H₀＝h₀-h₆c_x，H₁＝h₁-h₇c_x，H₂＝h₂-h₈c_x，

H₃＝h₃-h₆c_y，H₄＝h₄-h₇c_y，H₅＝h₅-h₈c_y，

H₆＝h₆，H₇＝h₇，H₈＝h₈，

其中(c_x,c_y)为相机主光轴点的坐标，h₀、h₁、h₂、h₃、h₄、h₅、h₆、h₇和h₈为单应性矩阵H_a中的元素，

H₀、H₁、H₂、H₃、H₄、H₅、H₆、H₇和H₈为单应性矩阵H_b中的元素，

根据以下约束关系由单应性矩阵H_b计算相机的等效焦距f_x和f_y：

10.一种计算机可读存储介质，其特征在于，所述介质上存储有程序，所述程序能够被处理器执行以实现如权利要求1-7中任一项所述的相机标定方法。

Images