CN110728718A - 改进摄像机标定各参数的方法 - Google Patents

Abstract

本发明涉及一种改进摄像机标定各参数的方法，利用参数待定的摄像机采集数张不同位置高精度标定板的标定图像；将采集的标定图像导入软件进行自适应伽马变换得到自适应对比度增强的标定板图像，对自适应对比度增强的标定板图像进行亚像素角点提取；将提取的标定板图像亚像素角点信息送入MATLAB标定工具箱，得到摄像机参数的线性解，此为初步参数；采用自适应权重和变异粒子群优化算法对摄像机参数的线性解进行非线性优化，得到精确标定参数。此方法获得较为精确的摄像机参数，标定误差小于0.1像素，满足工程中对摄像机高精度的要求，可以用于三维重建以及双目视觉立体标定中摄像机内外参数的确定。

Description

改进摄像机标定各参数的方法

技术领域

本发明涉及一种摄像机技术，特别涉及一种改进摄像机标定各参数的方法。

背景技术

近年来，计算机视觉在许多领域得到了广泛的应用。计算机视觉研究的一个重点是如何从摄像机获取的图像信息中得到对应物体的三维空间几何信息，由此实现三维空间物体的重建。摄像机标定是实现三维重建的关键。因此摄像机标定精度决定着空间点三维坐标的准确性。标定方法主要有传统标定法、主动视觉标定法和自标定法。传统标定方法标定精度高，但标定过程复杂，需要高精度的已知结构信息。主动视觉标定法通常可以线性求解，鲁棒性比较高，但是不适用于摄像机运动未知的场合。自标定方法仅需要建立图像之间的对应关系，灵活性强，但是鲁棒性不高。Zhang的平面标定方法是介于传统标定方法和自标定方法之间的一种方法。它既避免了传统方法设备要求高，操作繁琐等缺点，又较自标定方法精度高,因此广泛应用于摄像机标定中。基于Zhang的标定法主要步骤：首先，采集不同方位的标定板的图像；提取图像中的角点坐标；然后，利用线性技术产生内外参数的近似解，最后，采用非线性优化的方法改进线性解，得到优化后的摄像机参数。

许多优秀的研究者们为了进一步提高摄像机标定的精度做了很多的努力，并取得了良好的效果。但是大多数改进摄像机标定的方法仅局限于改进摄像机标定过程中的某一步骤，并没有考虑对标定过程进行整体改进，因此限制了标定精度的提高。大多数改进算法的标定误差仍在0.1像素以上，对一些对摄像机精度要求较高的场合，仍需要进一步改进。

发明内容

本发明是针对提高摄像机标定精度的问题，提出了一种改进摄像机标定各参数的方法，从图像采集，角点提取和摄像机参数非线性优化三个方面对传统的摄像机标定工具箱存在的缺陷进行改进。

本发明的技术方案为：一种改进摄像机标定各参数的方法，利用参数待定的摄像机采集数张不同位置高精度标定板的标定图像；将采集的标定图像导入软件进行自适应伽马变换得到自适应对比度增强的标定板图像，对自适应对比度增强的标定板图像进行亚像素角点提取；将提取的标定板图像亚像素角点信息送入MATLAB标定工具箱，得到摄像机参数的线性解；采用自适应权重和变异粒子群优化算法对摄像机参数的线性解进行非线性优化，得到精确标定参数。

所述得到自适应对比度增强的标定板图像具体步骤如下：

1.1)采用自适应伽马变换，调整图像亮度以及合理地增强黑白棋盘格的对比度，其基本的伽玛矫正的公式如下：

s＝c·r^γ,r∈[0,1]

其中r表示输入图像的灰度值，s表示输出图像的灰度值，c和γ为常数，设置c＝1；

1.2)为了自适应调整输入图像灰度值，利用图像的概率密度函数pdf调整γ值，设输入图像的灰度直方图为H(k_l,k_u)，其中k_l和k_u分别表示输入图像灰度值的下限和上限；

H(k)＝n_k,k＝k_l,...,k_u；

n_k表示图像中灰度级为k的个数，图像的概率密度函数pdf为：

其中N表示输入图像中的像素总数；

然后,自适应γ值表示为:

因此，自适应伽马变换可以表示为：

所述进行亚像素角点提取具体实现方法：采用基于单应性矩阵映射的亚像素角点提取，首先利用MATLAB标定工具箱选取棋盘格图像中的四个角点并得到角点的初始坐标；

采用亚像素坐标检测的方法得到精度较高的四个角点的坐标：假设角点的精确位置在q点,那么任意一个从q点到指向其邻域中的另一点p_i的向量,都与p_i处的图像灰度梯度垂直，由于噪声的存在,这两个正交矢量的点积不等于零,而存在误差，误差ε_i具体公式为：

式中，为p_i处的图像灰度梯度；p_i为p_i点的坐标；q点的值应该是使ε_i最小的那点的坐标；因此，建立q点邻域内所有p_i的方程，得到更精确的q坐标；以新的q为区域中心，继续使用该方法迭代将获得更高的精度；利用以上方法，由摄像机工具箱获取的标定板四个角点的初始坐标(x_t,y_t),t＝1,2,3,4，更新为(x_c,y_c),c＝1,2,3,4；利用以下式子可以求出单应性矩阵，

αm＝HM

其中m表示在图像中提取的角点坐标；M＝[X_W,Y_W,1]^T表示在标定板平面上对应角点的世界坐标，H∈R_3*3即为单应性矩阵，α为常数因子，因此当选取的每个图像中角点的数目大于等于4时，可以求出对应的单应性矩阵H；

通过单应性矩阵H，将在世界坐标系下指定的提取范围内的所有角点坐标映射到像素坐标系下，得到所有角点的像素坐标。

所述采用自适应权重和变异粒子群优化算法对摄像机参数的线性解进行非线性优化的具体步骤如下：

4.1)利用摄像机参数线性解作为N个粒子初始位置的均值，使粒子均匀分布在线性解周围，设置速度上下限，并在其范围内随机产生N个粒子的速度，设置最大迭代次数T以及加速因子c₁和c₂；

4.2)构建适应度函数：设定有a张标定图像，每张图像中有b个角点，每个角点的实际像素点坐标为(U_ij,V_ij),重投影的像素点坐标为(u_ij,v_ij)，其中下标i、j代表着第i张图片的第j个角点,适应度函数表达为：

4.3)为了衡量每个粒子的全局和局部搜索能力，定义动态个体搜索来动态调整惯性权重w,其具体公式如下：

其中x_id是第i个粒子在第d维的位置，pbest是自身最优解，gbest是群体最优解，ε是一个接近于0的正数，D代表求解的维数,每一维粒子的权重被重新定义为:

其中β是[0,1]之间的常数,当β＝0.3,该算法可以得到较佳的性能；

4.4)更新种群位置和速度，其具体公式如下：

其中t为迭代次数；w_id为惯性权重；c1和c2为学习因子；r1和r2为均匀分布于0和1之间的随机数；v_id(t)代表粒子第t次迭代的速度；v_id(t+1)代表粒子第t+1次迭代的速度；pbest是当前某个粒子发现的的最佳位置；gbest是当前整个种群中发现的所有粒子的最佳位置；

4.5)为了避免种群陷入局部最优，采用基于群体适应度方差确定最佳粒子gbest的变异概率，因此，变异操作增强算法跳出局部最优解的能力，粒子群的群体适应度方差定义为：

其中σ²为群体适应度方差，f是归一化因子，其作用是限制σ²的大小，f_i为第i个粒子的适应度，f_avg为q全体种群目前的平均适应度；群体适应度方差σ²反映了粒子群中所有粒子的“收敛”程度，σ²越小，则粒子群趋于收敛；反之，粒子群处于随机搜索阶段；当σ²＝0，粒子群算法陷入早熟收敛或则达到全局收敛；因此当σ²＝0时,并且此时得到的最优解不是理论最优解或者期望最优解，则可以断定粒子群陷入局部收敛，算法将出现早熟收敛；对满足变异条件的gbest按一定的概率p_m变异，p_m的计算公式如下：

其中k取[0.1,0.3]之间的任意数值，在此算法中p_m＝0.1；对于gbest变异操作采用增加随机扰动的方法；其具体公式为：

gbest_d＝(1+0.01*η)×gbest_d

其中gbest_d为gbest的第d维取值，η是服从Gauss(0，l)分布的随机变量；4.6)确定是否满足结束条件:如果达到最大迭代次数，则循环结束并输出结果；否则，跳到步骤4.3)继续迭代。

本发明的有益效果在于：本发明改进摄像机标定各参数的方法，此方法获得较为精确的摄像机参数，标定误差小于0.1像素，满足工程中对摄像机高精度的要求，可以用于三维重建以及双目视觉立体标定中摄像机内外参数的确定。

附图说明

图1为本发明改进摄像机标定各参数的方法步骤与验证示意图；

图2为本发明方法中在图像采集时，对标定图像进行自适应伽马变换前后的对比图；

图3为本发明方法中在标定图像角点提取时，与传统Harris角点提取局部放大对比图；

图4为本发明方法中在相机参数非线性优化时，与基本粒子群算法(PSO)、差分进化算法(DE)、遗传算法(GA)收敛性能的对比图；

图5为本发明方法中在利用重投影误差验证标定有效性时，实际角点像素坐标与重投影像素坐标的相对误差图。

具体实施方式

以下结合附图说明对本发明的实施例作进一步详细描述，但本实施例并不用于限制本发明，凡是采用本发明的相似结构及其相似变化，均应列入本发明的保护范围，本发明中的顿号均表示和的关系。

如图1所示改进摄像机标定各参数的方法步骤与验证示意图，软硬件所需要的配置有参数待定的摄像机、高精度标定板、MATLAB标定工具箱，具体步骤如下：利用摄像机采集12张不同位置高精度标定板的标定图像；将采集的标定图像导入到MATLAB软件中进行自适应对比度增强以及标定图像亚像素角点的提取；然后采用MATLAB标定工具箱实现摄像机参数(光心、焦距、畸变系数)的初步确定，最后采用非线性优化的方法确定更为精确的摄像机参数。

1、自适应对比度增强

1.1、考虑到在采集标定图片时，由于自然光照的不确定性，可能出现所采集的图像过度曝光或者曝光不足的现象。采用自适应伽马变换，调整图像亮度以及合理地增强黑白棋盘格的对比度。其基本的伽玛矫正的公式如下：

s＝c·r^γ,r∈[0,1]

其中r表示输入图像的灰度值，s表示输出图像的灰度值，c和γ为常数，设置c＝1。

1.2、为了自适应调整输入图像灰度值，利用图像的概率密度函数(pdf)调整γ值。设输入图像的灰度直方图为H(k_l,k_u)，其中k_l和k_u分别表示输入图像灰度值的下限和上限。

H(k)＝n_k,k＝k_l,...,k_u

n_k表示图像中灰度级为k的个数，图像的概率密度函数(pdf)为：

其中N表示输入图像中的像素总数。然后,自适应γ值可以表示为:

因此，自适应伽马变换可以表示为：

不同的γ值，对图像处理的效果不同：

A、γ＝1时，输出图像灰度值与输入图像灰度值相同；

B、γ<1时，过暗的像素得到增强；

C、γ>1时，过亮的像素起到抑制。

如图2所示在图像采集时，对标定图像进行自适应伽马变换前后的对比图，左图为原标定图像，右图为自适应伽马变换后的图像，经过自适应伽马变换，合理地调节了图像亮度，增强了标定图像中黑白网格的对比度，为接下来角点的提取奠定了基础。

2、对标定板图像进行亚像素角点提取：

为了避免Harris角点提取的漏检和错检的情况，采用基于单应性矩阵映射的亚像素角点提取。首先利用标定工具箱选取棋盘格图像中的四个角点并得到角点的初始坐标。为了进一步得到精度较高的四个角点的坐标，采用亚像素坐标检测的方法。假设角点的精确位置在q点,那么任意一个从q点到指向其邻域中的另一点p_i的向量,都与p_i处的图像灰度梯度垂直。由于噪声的存在,这两个正交矢量的点积不等于零,而存在误差。误差ε_i具体公式为：

式中，

为p_i处的图像灰度梯度；p_i为p_i点的坐标；q点的值应该是使ε_i最小的那点的坐标。因此，建立q点邻域内所有p_i(假设n个点)的方程，得到更精确的q坐标。以新的q为区域中心，继续使用该方法迭代将获得更高的精度。利用以上方法，由摄像机工具箱获取的标定板四个角点的初始坐标(x_t,y_t),(t＝1,2,3,4)更新为(x_c,y_c),(c＝1,2,3,4)。利用以下式子可以求出单应性矩阵。

αm＝HM

其中m表示在图像中提取的角点坐标；M＝[X_W,Y_W,1]^T表示在标定板平面上对应角点的世界坐标。H∈R_3*3即为单应性矩阵。α为常数因子。因此当选取的每个图像中角点的数目大于等于4时，可以求出对应的单应性矩阵H。

通过单应性矩阵H，将在世界坐标系下指定的提取范围内的所有角点坐标映射到像素坐标系下，得到所有角点的像素坐标。

3、将提取的标定板图像亚像素角点信息送入MATLAB标定工具箱，得到摄像机参数的线性解。

4、摄像机参数的非线性优化

采用自适应权重和变异粒子群优化算法(AWMPSO)对摄像机参数的线性解进行非线性优化。

4.1)利用摄像机参数线性解作为N个粒子初始位置的均值，使粒子均匀分布在线性解周围。设置速度上下限，并在其范围内随机产生N个粒子的速度.设置最大迭代次数T以及加速因子c₁和c₂。

4.2)构建适应度函数：设定有a张标定图像，每张图像中有b个角点，每个角点的实际像素点坐标为(U_ij,V_ij),重投影的像素点坐标为(u_ij,v_ij)，其中下标i、j代表着第i张图片的第j个角点。适应度函数表达为：

4.3)为了衡量每个粒子的全局和局部搜索能力，定义动态个体搜索(ISA)来动态调整惯性权重w。其具体公式如下：

其中x_id是第i个粒子在第d维的位置，pbest是自身最优解，gbest是群体最优解，ε是一个接近于0的正数，D代表求解的维数。每一维粒子的权重被重新定义为:

其中β是[0,1]之间的常数。当β＝0.3,该算法可以得到较佳的性能。

4.4)更新种群位置和速度，其具体公式如下：

其中t为迭代次数；w_id为惯性权重；c1和c2为学习因子；r1和r2为均匀分布于0和1之间的随机数；v_id(t)代表粒子第t次迭代的速度；v_id(t+1)代表粒子第t+1次迭代的速度。pbest是当前某个粒子发现的的最佳位置。gbest是当前整个种群中发现的所有粒子的最佳位置。

4.5)为了避免种群陷入局部最优，采用基于群体适应度方差确定最佳粒子gbest的变异概率。因此，变异操作增强算法跳出局部最优解的能力。粒子群的群体适应度方差定义为：

其中σ²为群体适应度方差，f是归一化因子，其作用是限制σ²的大小，f_i为第i个粒子的适应度，f_avg为q全体种群目前的平均适应度。群体适应度方差σ²反映了粒子群中所有粒子的“收敛”程度，σ²越小，则粒子群趋于收敛；反之，粒子群处于随机搜索阶段。当σ²＝0，粒子群算法陷入早熟收敛或则达到全局收敛。因此当σ²＝0时,并且此时得到的最优解不是理论最优解或者期望最优解，则可以断定粒子群陷入局部收敛，算法将出现早熟收敛。对满足变异条件的gbest按一定的概率p_m变异，p_m的计算公式如下：

其中k可以取[0.1,0.3]之间的任意数值，在此算法中p_m＝0.1。对于gbest变异操作采用增加随机扰动的方法。其具体公式为：

gbest_d＝(1+0.01*η)×gbest_d

其中gbest_d为gbest的第d维取值，η是服从Gauss(0，l)分布的随机变量。

4.6)确定是否满足结束条件:如果达到最大迭代次数，则循环结束并输出结果；否则，跳到步骤4.3)继续迭代。

通过以上三种方法得到了摄像机的内外参数，为了验证标定的效果，利用标定图像中角点重投影坐标与实际角点坐标的相对误差进行衡量。其中采用本发明方法的适应度值最后稳定在0.037，PSO的适应度值稳定在0.319。

如图5是发明实施例的改进摄像机标定技术中，在利用重投影误差验证标定有效性时，实际角点像素坐标与重投影像素坐标的相对误差图。其中不同灰度的十字架(+)，代表不同标定图片的重投影相对误差，可以看出大部分的误差范围在[-0.25,0.2]之间，更多的误差范围控制在[-0.1,0.1]之间，从而证明了采用以上方法可以提高摄像机标定参数的参数。

Claims (4)

1.一种改进摄像机标定各参数的方法，其特征在于，利用参数待定的摄像机采集数张不同位置高精度标定板的标定图像；将采集的标定图像导入软件进行自适应伽马变换得到自适应对比度增强的标定板图像，对自适应对比度增强的标定板图像进行亚像素角点提取；将提取的标定板图像亚像素角点信息送入MATLAB标定工具箱，得到摄像机参数的线性解；采用自适应权重和变异粒子群优化算法对摄像机参数的线性解进行非线性优化，得到精确标定参数。

2.根据权利要求1所述改进摄像机标定各参数的方法，其特征在于，所述得到自适应对比度增强的标定板图像具体步骤如下：

1.1)采用自适应伽马变换，调整图像亮度以及合理地增强黑白棋盘格的对比度，其基本的伽玛矫正的公式如下：

s＝c·r^γ,r∈[0,1]

其中r表示输入图像的灰度值，s表示输出图像的灰度值，c和γ为常数，设置c＝1；

1.2)为了自适应调整输入图像灰度值，利用图像的概率密度函数pdf调整γ值，设输入图像的灰度直方图为H(k_l,k_u)，其中k_l和k_u分别表示输入图像灰度值的下限和上限；

H(k)＝n_k,k＝k_l,...,k_u；

n_k表示图像中灰度级为k的个数，图像的概率密度函数pdf为：

其中N表示输入图像中的像素总数；

然后,自适应γ值表示为:

因此，自适应伽马变换可表示为：

3.根据权利要求1所述改进摄像机标定各参数的方法，其特征在于，所述进行亚像素角点提取具体实现方法：采用基于单应性矩阵映射的亚像素角点提取，首先利用MATLAB标定工具箱选取棋盘格图像中的四个角点并得到角点的初始坐标；

采用亚像素坐标检测的方法得到精度较高的四个角点的坐标：假设角点的精确位置在q点,那么任意一个从q点到指向其邻域中的另一点p_i的向量,都与p_i处的图像灰度梯度垂直，由于噪声的存在,这两个正交矢量的点积不等于零,而存在误差，误差ε_i具体公式为：

式中，

为p_i处的图像灰度梯度；p_i为p_i点的坐标；q点的值应该是使ε_i最小的那点的坐标；因此，建立q点邻域内所有p_i的方程，得到更精确的q坐标；以新的q为区域中心，继续使用该方法迭代将获得更高的精度；利用以上方法，由摄像机工具箱获取的标定板四个角点的初始坐标(x_t,y_t),t＝1,2,3,4，更新为(x_c,y_c),c＝1,2,3,4；利用以下式子可以求出单应性矩阵，

αm＝HM

其中m表示在图像中提取的角点坐标；M＝[X_W,Y_W,1]^T表示在标定板平面上对应角点的世界坐标，H∈R_3*3即为单应性矩阵，α为常数因子，因此当选取的每个图像中角点的数目大于等于4时，可以求出对应的单应性矩阵H；

通过单应性矩阵H，将在世界坐标系下指定的提取范围内的所有角点坐标映射到像素坐标系下，得到所有角点的像素坐标。

4.根据权利要求1所述改进摄像机标定各参数的方法，其特征在于，所述采用自适应权重和变异粒子群优化算法对摄像机参数的线性解进行非线性优化的具体步骤如下：

4.1)利用摄像机参数线性解作为N个粒子初始位置的均值，使粒子均匀分布在线性解周围，设置速度上下限，并在其范围内随机产生N个粒子的速度，设置最大迭代次数T以及加速因子c₁和c₂；

4.2)构建适应度函数：设定有a张标定图像，每张图像中有b个角点，每个角点的实际像素点坐标为(U_ij,V_ij),重投影的像素点坐标为(u_ij,v_ij)，其中下标i、j代表着第i张图片的第j个角点,适应度函数表达为：

4.3)为了衡量每个粒子的全局和局部搜索能力，定义动态个体搜索来动态调整惯性权重w,其具体公式如下：

其中x_id是第i个粒子在第d维的位置，pbest是自身最优解，gbest是群体最优解，ε是一个接近于0的正数，D代表求解的维数,每一维粒子的权重被重新定义为:

其中β是[0,1]之间的常数,当β＝0.3,该算法可以得到较佳的性能；

4.4)更新种群位置和速度，其具体公式如下：

其中t为迭代次数；w_id为惯性权重；c1和c2为学习因子；r1和r2为均匀分布于0和1之间的随机数；v_id(t)代表粒子第t次迭代的速度；v_id(t+1)代表粒子第t+1次迭代的速度；pbest是当前某个粒子发现的的最佳位置；gbest是当前整个种群中发现的所有粒子的最佳位置；

4.5)为了避免种群陷入局部最优，采用基于群体适应度方差确定最佳粒子gbest的变异概率，因此，变异操作增强算法跳出局部最优解的能力，粒子群的群体适应度方差定义为：

其中σ²为群体适应度方差，f是归一化因子，其作用是限制σ²的大小，f_i为第i个粒子的适应度，f_avg为q全体种群目前的平均适应度；群体适应度方差σ²反映了粒子群中所有粒子的“收敛”程度，σ²越小，则粒子群趋于收敛；反之，粒子群处于随机搜索阶段；当σ²＝0，粒子群算法陷入早熟收敛或则达到全局收敛；因此当σ²＝0时,并且此时得到的最优解不是理论最优解或者期望最优解，则可以断定粒子群陷入局部收敛，算法将出现早熟收敛；对满足变异条件的gbest按一定的概率p_m变异，p_m的计算公式如下：

其中k取[0.1,0.3]之间的任意数值，在此算法中p_m＝0.1；对于gbest变异操作采用增加随机扰动的方法；其具体公式为：

gbest_d＝(1+0.01*η)×gbest_d

其中gbest_d为gbest的第d维取值，η是服从Gauss(0，l)分布的随机变量；

4.6)确定是否满足结束条件:如果达到最大迭代次数，则循环结束并输出结果；否则，跳到步骤4.3)继续迭代。

Images