2017, Vol. 38
2. 海信集团有限公司，山东 青岛 266000;
3. 国家电网公司东北分部，辽宁 沈阳 110180
2. Hisense Group Co., Ltd., Qingdao 266000, China;
3. Northeast Branch of State Grid Corporation of China, Shenyang 110180, China
摄像机获取的二维图像与对应的三维物体的位置关系是由摄像机成像的几何模型决定的, 几何模型中的参数就是摄像机参数, 由试验和计算求解这些参数的过程就是摄像机标定.摄像机标定是三维重建技术的基本问题, 同时也是立体视觉研究的重要组成部分.摄像机标定结果的好坏直接影响三维空间测量的精度和重建的结果[1].由于摄像机镜头制造工艺的原因, 使得成像系统不能在整个视场上完全满足针孔成像模型, 从而导致投影射线变为曲线, 造成成像平面上像素点实际坐标与理论坐标不一致.因此在拍摄图像时往往会产生非线性畸变, 必须进行畸变校正来提高图像几何位置精度.
目前, 对摄像机镜头进行畸变校正的方法主要分为两大类.一类是将畸变系数与摄像机内外参数融合在一起, 畸变系数的求解过程就是非线性的优化过程.其中两步法采用线性方法求解大部分参数, 采用迭代方法求解少数参数, 因而迭代少, 计算速度较快, 但畸变模型简单, 对畸变问题的解决不够理想[2];基于二维平面靶标的摄像机标定方法不需要特殊的标定物, 因此, 标定过程比较灵活, 但是计算过程相对复杂[3].另一类畸变校正方法是将畸变系数从摄像机模型中分离出来, 从而避免将其重复标定.其中基于共线点的畸变分离模型[4]仅考虑一阶径向畸变, 且要求图像中心点已知; 用拟合圆来求解图像中心点和一阶径向畸变系数的方法[5]虽然可以求出图像中心点, 但是由于采用单参数除式模型, 不能求解二阶径向畸变系数.
针对以上方法存在的不足, 并考虑具有二阶径向畸变的摄像机模型, 本文首先采用模拟退火[6]与粒子群[7]相结合的优化算法 (SA-PSO) 标定畸变系数和图像中心点, 将畸变形成的围线面积作为畸变测度, 然后利用求得的畸变系数和图像中心点坐标标定摄像机的其他参数.
1 摄像机畸变模型与参数摄像机的针孔成像模型示意图如图 1所示, 其中OwXwYwZw是用来表示摄像机位置和实际物体位置的世界坐标系, O1XY是摄像机获取图像的成像平面坐标系, O0uv是以像素为单位表示图像的像素坐标系, OcXcYcZc是与摄像机镜头的光心、光轴相关的摄像机坐标系.
|
图 1 摄像机畸变模型 Fig.1 Camera model with distortion |
为了标定方便, 将平面模板设立在Zw=0的二维世界坐标系平面上.由摄像机的针孔成像模型可得, 二维平面上的点与对应的图像中的点之间的关系可以写成[2]:
|
(1) |
式中: K为摄像机内部参数矩阵, 其中au, av分别为u轴、v轴上的尺度因子; (u0, v0) 是摄像机光轴与图像平面交点 (称为图像中心点) 的坐标.而r1=[r11 r21 r31]T, r2=[r12 r22 r32]T分别指Xw轴、Yw轴在摄像机坐标系OcXcYcZc中的方向; t=[t1 t2 t3]T是OwXwYwZw的坐标原点在摄像机坐标系OcXcYcZc中的位置, 因此[R t]称为摄像机外部参数矩阵.
对于一般的定焦摄像机来说, 考虑二阶径向畸变就能满足标定的精度要求[8].本文考虑摄像机二阶径向畸变误差, 如图 1中P为三维空间中任意一点, 则Pu(xu, yu) 和Pd(xd, yd) 分别为P点在成像平面上的理想坐标和畸变图形下的实际坐标, 二者之间的关系可表示为
|
(2) |
其中δx和δy是非线性畸变值, 可表示为
|
(3) |
其中:r2=xd2+yd2; k1, k2即为待标定的二阶径向畸变系数.
2 粒子群算法的改进与应用粒子群算法是一种迭代的启发式优化算法, 由于算法简单、容易理解，并且参数少, 现已广泛应用于多目标问题、约束问题、动态问题等实际应用问题中.粒子群优化算法的原理:一个由m个粒子组成的群体在D维搜索空间中以确定的速度飞行, 每个粒子在搜索时记录自己搜索到的历史最优点和群体内以及邻域内其他粒子的历史最优点, 然后在此基础上进行位置变化.
在第t代, 粒子群的每个粒子都是由位置、历史最优位置和速度3个D维向量组成, 第i个粒子的3个向量分别表示为
|
(4) |
|
(5) |
|
(6) |
其中：i=1, 2, …, m, d=1, 2, …, D.
在算法每一次迭代过程中, 对当前位置xidt与历史最优位置进行比较.如果当前位置优于历史最优位置, 那么就用当前位置替换历史最优位置.另外将整个粒子群中迄今为止搜索到的最优位置记为gbest, 即全局最优解.
每个粒子的位置和速度更新公式为
|
(7) |
|
(8) |
其中:w是惯性权重, 其值对算法的全局搜索能力有很大影响; c1, c2是学习因子, 分别影响全局搜索的最大步长和最优粒子的飞行方向; λ1, λ2是 (0~1) 之间的随机数.由于粒子位置优劣的判断需要评测函数, 因此将改进的粒子群用于摄像机畸变参数和图像中心点标定时, 需要定义评测函数.
2.1 定义评测函数由透视投影理论可知, 对棋盘格平面模板提取特征点时, 任意一条特征点形成的直线在像平面所成的像仍然是一条直线.但实际上, 摄像机镜头存在的畸变会导致直线所成的像变成曲线[9-10].当定焦摄像机仅考虑二阶径向畸变时, 特征点所成特征直线通常会变成凸线或凹线, 畸变成S型线的情况非常少见, 可以忽略.因此可以连接凸线或者凹线的两个端点, 通过计算围线的面积来判断畸变校正的程度.
本文提出一种简便的计算围线面积的方法.以凸线为例, 如图 2所示, L1是畸变特征点所形成的特征曲线, L2是连接曲线端点的直线.A (x1, y1), B (x2, y2) 是L1上任意相连的两个特征点.过A, B两点分别做x轴的垂线, 与L2相交于C, D两点.c1, c2分别为C, D两点纵坐标.E, F是位于曲线两个端点处的特征点.
|
图 2 围线和畸变曲线 Fig.2 Contour and distorted curve |
图 2中阴影部分的面积可表示为
|
(9) |
当特征曲线是凹线时, 相应的阴影面积为
|
(10) |
由式 (9) 和式 (10) 可以看出, 特征曲线分别为凸线和凹线时, 各点的相对位置发生变化, 阴影面积计算结果的表达式只差一个负号.由于特征直线为凹线时, (y1+y2) 小于 (c1+c2), 所以不管特征曲线是凹线还是凸线, 面积都是正值.
假设棋盘格模板有M条竖直线, 每条直线上有N个特征点, xij和yij分别表示第i条特征曲线的第j个特征点的横坐标和纵坐标.过第i条特征曲线的第j个特征点做Y轴的平行线与连接特征曲线两端点的直线交于一点, 这点的纵坐标就是cij.
|
(11) |
|
(12) |
则畸变评测函数可定义为
|
(13) |
将畸变评测函数作为粒子群优化算法的目标函数, 则畸变评测函数越小代表畸变校正的精度越高.为了求得畸变评测函数的最小值, 本文提出一种模拟退火与粒子群相结合的优化算法 (SA-PSO).
2.2 基于SA-PSO的畸变系数标定本文使用的摄像机畸变模型中, 畸变系数以及图像中心点的求解属于非线性优化问题.由于传统粒子群优化算法的惯性权重w和学习因子c1, c2是确定的值, 使得算法收敛速度慢, 且极易陷入局部最优.而模拟退火算法以概率的方式进行搜索, 增加了搜索过程的灵活性, 具有很强的全局搜索能力.因此为了克服粒子群算法的缺点, 本文用模拟退火原理优化粒子群算法中的惯性权重w和学习因子c1, c2这3个参数, 然后用优化后的粒子群算法求解摄像机的畸变系数和图像中心点.
结合模拟退火的改进粒子群算法步骤如下:
1) 粒子群优化算法初始化.
① 设定种群规模m=50, 最大迭代次数Tmax=1000, 畸变参数k1, k2与图像中心点u0, v0共有4个待求参数, 所以粒子维数D=4, 每个粒子的位置对应着待求的4个参数.
② 设置粒子群的初始速度vid0和初始位置xid0.
③ 将式 (13) 作为目标函数, 计算初始粒子群的目标函数值f (xi0), 然后求出相应的初始个体最优pid和初始群体最优位置gbest.
2) 模拟退火初始化.
① 将第t代目标函数最小值记为f (xit*), 设置初始温度
② 令评价函数c (s)=f (xi0*).
3) 由摄动装置生成新的s′(w′, c′1, c′2).
4) 由式 (7)、式 (8) 更新粒子的位置, 其中, w, c1, c2按s′取值, 更新位置后, 计算目标函数.
5) 令c (s′)=f (xit*), Δc=c (s′)－c (s).若Δc < 0或e－Δc/T > rand (0, 1), 则c (s)=c (s′), s=s′, T=f (xit*)/ln5, 并接受由s′(w′, c′1, c′2) 所更新的速度和位置.否则, 不接受s′的值, 用s (w, c1, c2) 更新速度和位置, 并计算目标函数值.
6) 根据目标函数值, 更新pid, gbest.
7) 若满足结束条件, 输出最优位置, 若不满足, 转到步骤3).
算法最后输出的最优位置gbest=(g1, g2, g3, g4) 中的数值分别对应摄像机的畸变系数k1, k2和两个内参数u0, v0.
2.3 摄像机标定用改进的粒子群算法 (SA-PSO) 已经求出摄像机的畸变系数和图像中心点坐标.其他的内部参数和外部参数可用如下的方法标定.
设靶标平面与摄像机之间的单映性矩阵为H, 则式 (1) 可以进一步写成:
|
(14) |
由文献[2]可得, 单映性矩阵对内参数约束为
|
(15) |
对式 (15) 展开, 可以求出另外两个摄像机内参数au和av的值[2].
由内参数矩阵K和单映性矩阵H的求解结果, 可计算出每幅图像的外部参数为[2]
|
(16) |
将本文提出的标定算法与经典的张正友标定方法[2]进行对比分析, 实验用到的CCD摄像机, 镜头的光圈系数 (F) 为1.4, 焦距 (f) 为25 mm, 图像分辨率为1024×768像素.不失一般性, 对标定模板从不同角度拍摄3幅图像, 如图 3所示.
|
图 3 畸变图像 Fig.3 Distorted images (a)—角度1；(b)—角度2; (c)—角度3. |
用Harris角点检测[11]提取特征点, 然后用本文提出的改进方法求解畸变系数和图像中心点以及由式 (15) 求出的内部参数，结果如表 1所示.
|
表 1 本文方法标定的摄像机内部参数 Table 1 Camera intrinsic parameters for the proposed |
而用张正友标定方法求得的摄像机内部参数如表 2所示.
|
|
表 2 张正友方法标定的摄像机内部参数 Table 2 Camera intrinsic parameters for Zhang's method |
采用本文方法, 由式 (16) 求出的对应于图 3中3幅图像的摄像机外部参数和张正友方法标定的外部参数结果如表 3所示.
|
|
表 3 摄像机外部参数标定结果 Table 3 Calibration results of camera extrinsic parameters |
为了验证两种标定算法的标定精度, 随机从第一幅图中抽取10组检测点, 计算图像坐标的反投影坐标值.经式 (1) 计算, 本文方法的图像坐标的反投影误差在u轴和v轴上都小于1个像素, 均值为0.514个像素.而张正友方法计算图像坐标的反投影误差在u轴和v轴上均值为1.189个像素.这说明本文提出的方法标定精度很高, 具有一定的可行性.
为了证明本文方法的鲁棒性, 在实验中, 对图像提取的各个特征点分别加入均匀分布的噪声, 噪声级从0(像素) 以0.25的步长依次增加至2.0(像素), 分别用本文提出的方法和张正友提出的方法计算每幅图像在当前噪声级下相应世界坐标的平均绝对误差, 结果分别如图 4、图 5所示.
|
图 4 SA-PSO方法的误差 Fig.4 Error of SA-PSO method |
|
图 5 张正友方法的误差 Fig.5 Error of Zhang's method |
从图 4、图 5可以看出, 随着噪声的增大, 误差结果也随着变大, 但仍能维持在一个很小的误差范围内, 在相同的噪声级下, 本文提出的方法计算的误差比张正友方法计算的误差要小, 说明本文提出的方法稳定性较好.
4 结论1) 为了避免将畸变系数和图像中心点纳入摄像机模型中进行多次重复标定, 本文提出一种畸变分离的摄像机模型, 研究出新的评判函数.
2) 利用经过模拟退火优化的粒子群算法找到最优的畸变系数和图像中心点, 然后标定摄像机剩余的参数.
3) 与经典的张正友方法相比, 本文方法更加灵活, 适用范围更广, 稳定性和精度都有所提高.仿真实验表明该方法需要的已知参数少, 容易实现, 而且图像坐标的平均反投影误差明显减小.
| [1] | Liu J Y, Li Y F, Chen S Y. Robust camera calibration by optimal localization of spatial control points[J]. Transactions on Instrumentation and Measurement, 2014, 63(12): 3076–3087. DOI:10.1109/TIM.2014.2324792 |
| [2] | Tsai R Y. A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf TV cameras and lenses[J]. Robotics and Automation, 1987, 3(4): 323–344. |
| [3] | Zhang Z. A flexible new technique for camera calibration[J]. Pattern Analysis and Machine Intelligence, 2000, 22(11): 1330–1334. DOI:10.1109/34.888718 |
| [4] |
周富强, 胡坤, 张广军.
基于共线特征点的摄像机镜头畸变校正[J]. 机械工程学报, 2006, 42(9): 174–177.
( Zhou Fu-qiang, Hu Kun, Zhang Guang-jun. Correction distortion of camera lens with collinear points[J]. The Chinese Journal of Mechanical Engineering, 2006, 42(9): 174–177. ) |
| [5] | Fitzgibbon A W.Simultaneous linear estimation of multiple view geometry and lens distortion[C]//Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition.Hawaii, 2001:125-132. |
| [6] | Gomez D, Prieto F, Guzman M. Nearest neighbors by adaptive simulated annealing[J]. IEEE Latin America Transactions, 2015, 13(7): 2398–2404. DOI:10.1109/TLA.2015.7273804 |
| [7] | Wang H, Yen G G. Adaptive multiobjective particle swarm optimization based on parallel cell coordinate system[J]. Evolutionary Computation, 2015, 19(1): 1–18. DOI:10.1109/TEVC.2013.2296151 |
| [8] | Ma L, Chen Y Q, Moore K L. Analytical piecewise radial distortion model for precision camera calibration[J]. Vision, Image and Signal Processing, 2006, 153(4): 468–474. DOI:10.1049/ip-vis:20045035 |
| [9] | Rahman T, Krouglicof N. An efficient camera calibration technique offering robustness and accuracy over a wide range of lens distortion[J]. Transactions on Image Processing, 2013, 21(2): 626–637. |
| [10] | Lourenco M, Barreto J P, Vasconcelos F. sRD-SIFT:key point detection and matching in images with radial distortion[J]. Transactions on Robotics, 2012, 28(3): 752–760. DOI:10.1109/TRO.2012.2184952 |
| [11] | Hsiao P, Lu C L, Fu L C. Multilayered image processing for multiscale harris corner detection in digital realization[J]. Transactions on Industrial Electronics, 2010, 57(5): 1799–1850. DOI:10.1109/TIE.2010.2040556 |