商 海,倪受东,苏智勇
(南京工业大学 机械与动力工程学院,江苏 南京 211816)
机器视觉目前已经应用于很多方面,比如工业、农业、医学、科研等。然而,因为摄像机本身的精度问题,致使其无法应用于许多精度要求较高的场所中。
摄像机拍摄的原理是将空间中的物体映射到像素二维坐标系上,而清楚世界坐标系映射到像平面坐标系上的过程与关系就是摄像机标定的目的[1]。目前常用的摄像机标定技术大致有三种:直接线性变换法(也称DLT标定法)、两步法(也称Tsai标定法)以及张正友标定法。DLT标定法虽然标定方法简单、容易实现标定,但是却忽略了摄像机成像过程中产生的切向畸变和径向畸变[2]。正因为这两种对畸变系数的影响,导致使用DLT标定法精度比较低。Tsai标定法相比于DLT标定法虽然考虑了摄像机成像过程中产生的畸变系数,并且精度相比于DLT标定法高,但是此标定方法需要利用到参照物,如果在不同环境下,其参照物需要重新标定,故此方法比较难以应用[3]。张正友标定法优点是操作简单、设备要求低、精度高,缺点是稳健性较差[4]。
正因为目前常用的摄像机标定方法都有着其优点和缺点,使得近年来越来越多的科研人员致力于寻找精度高、成本低并且容易实现的标定方法[5]。目前对于摄像机标定的提高也采取了很多措施,其中利用各种算法对摄像机标定进行重新标定的方法使得精度越来越高,比如郭彤颖等人利用粒子群优化算法对摄像机进行标定[6];王子豪等人利用改进LM算法对摄像机进行标定[7];张宏峰等人利用麻雀搜索算法(SSA)对摄像机进行标定[8]。虽然这些标定方法简单方便,并且成本也相对较低,但是总的来说标定精度还是不高。为此,这里对麻雀算法进行了改进,然后进行摄像机标定,进而得到比较高的标定精度[9]。
所谓摄像机标定的原理就是将三维空间里面的某一个物体或者某一特征点的坐标转化为二维图像上的坐标从而获取摄相机参数。在这一映射条件下,主要涉及世界坐标系、摄像机坐标系、像平面坐标系和图像坐标系[10]。
空间中的点在映射过程中四种坐标系之间的关系如图1所示。
图1 四种坐标系之间的关系
图1可以简洁地反映出摄像机成像过程,实际上是通过世界坐标系(O-XYZ)至摄像机坐标系(o-xyz)、摄像机坐标系至像平面坐标系(o'-x'y')、像平面坐标系至图像系(O0-uv)三种转换完成的,可以表示成:
(1)
式中,u0、v0表示摄像机的主点在像平面系下的横、纵坐标;1/dx、1/dy分别表示图像坐标系中横、纵坐标的单位长度上的像素个数,单位为mm/pixel;f表示摄像机的焦距,可以手动设定;R3×3、T3×3分别表示采用摄像机的旋转矩阵和平移向量;(X,Y,Z)表示世界坐标系中的点;(u,v)表示像平面坐标系中的点;z表示摄像机坐标系中z轴坐标。
上述的计算公式是在没有任何干扰情况下的成像结果。然而,在实际工作当中,某一物体或者某一特征点真实成像的位置与理想成像位置之间存在一定程度上的偏差,这种偏差是摄像机的非线性光学畸变造成的[11]。图2就是由于畸变产生的位置偏移。
图2 由于畸变产生的位置偏移现象
从图2可以看出,世界坐标系上某一点(x,y,z)在像素坐标系上的理想成像点为(Xc,Yc),但是因为摄像机畸变,造成点成像位置移动到了(Xd,Yd)。为了让摄像机畸变产生偏移的点尽可能地接近理想状态下的映射点,需要对摄像机产生的畸变系数进行修正[10]。对应关系为:
(2)
式中,δx、δy、θx、θy为摄像机在x、y轴上的径向畸变和切向畸变。而畸变模型公式如下:
(3)
(4)
由式(2)~式(4)可得特征点理想位置与真实位置之间的关系为:
(5)
故根据公式(1)与(5)可得世界坐标系上的点(x,y,z),经过坐标变换与畸变之后编程实际像素坐标(Xd,Yd)。
麻雀在自然界中作为一种群居鸟类,无论是在发现食物、找寻食物还是发现天敌中都具有十分明显的分工。在整个麻雀种群中主要分为发现者、加入者与侦察者三种麻雀。在麻雀种群中,发现者主要起到发现食物并带领其他麻雀前往觅食地的作用,剩余麻雀则通过前者来获取食物。与此同时,侦察者在麻雀种群中起到观察觅食环境是否存在危险并提醒麻雀种群的作用。
麻雀搜索算法(SSA)正是因为麻雀种群中存在着这样一系列的分工与合作而产生的。在SSA中,发现者相比于种群中的其他麻雀而言具有较好的适应度值,并且其需要为整个麻雀种群提供觅食方向。发现者位置更新表达式为:
(6)
式中,n代表目前发现者位置更新所迭代的次数;Nmax代表发现者位置更新所设定的最大迭代次数;Xi1,j1表示此时第i1个麻雀在第j1维中的位置;α为其中的一个随机数,范围为α∈(0,1];R2代表预警值,其范围为R2∈[0,1];Ma代表安全值,其范围为Ma∈[0.5,1];Q是一个随机数,其服从正态分布;L表示一个元素全是1的1×d的矩阵。当R2 加入者的位置更新如下: (7) 式中,X1代表在加入者加入种群后发现者所处的最优位置;X2代表加入者加入种群后全局最差位置;A代表一个1×d的矩阵,且A中所有元素都随机在1或-1中选择,A矩阵有着A+=AT(AAT)-1的对应关系;n表示麻雀种群的大小。当i>n/2时,表示加入者没有找到此刻整个麻雀种群的最优位置,所以需要去别的地方继续寻找最优位置。 位于整个麻雀种群外围的麻雀起到侦察的作用,也就是三种麻雀中的侦察者,主要用于侦察周围环境是否存在危险。侦察者位置在每次迭代过程中所更新的表达式为: (8) 式中,X3表示目前麻雀种群中的最佳位置;β与K都是步长控制参数,但二者取值不同,β为服从方差为1、均值为0的正态分布的随机值,而K∈[-1,1];fi3表示迭代次数为i3时的适应度值;f1代表当前全局最优的适应度值;f2代表当前全局最差的适应度值;ε表示最小常数。 麻雀搜索算法简易流程如图3所示。 图3 麻雀搜索算法流程 根据上述发现者、加入者更新公式可以看出来,麻雀种群中麻雀个体收敛于全局最优解是直接从当前位置跳到最优解的附近,虽然最后得到的结果正确,但是直接跳跃到最优解附近就会导致除此之外的解出现偏差,也就是局部最优,从而导致精度降低。为了提高其精度,这里将鸟群算法和SSA进行结合和改进,利用改进后的算法可以避免出现上述的问题。 鸟群算法源于鸟群的三种行为:觅食行为、警戒行为和飞行行为。鸟群算法的目的为解决优化问题,其一些理想化规则如下: (1)鸟群中的每只鸟都可以在觅食与警戒行为中随意切换,且其觅食或者警戒行为都为随机策略。 (2)鸟群中的每只鸟都会在觅食过程中及时记录当前最佳的觅食位置,并且将其信息传递给鸟群中的其他个体,同时记录下此时鸟群的最佳觅食位置。 (3)当鸟群中每只鸟处于警戒行为下的时候,其都会做出往鸟群中心位移的趋势,但是这种行为会引起种群中竞争,因此储备量高的鸟相比于储备量低的鸟有着更容易接近鸟群中心的机会。 (4)鸟群会在一定时间内飞往另外一个地点。当飞往另外一个地点时,食物储备量高的个体会变成发现者,与此同时食物储备量低的个体会变成加入者,食物储备量位于这两者之间的个体将会随机变成发现者和加入者。 (5)发现者会在这片区域内努力寻找食物,而加入者则会随机跟随发现者进行觅食。 鸟群中每个鸟的各种行为可以描述为: (1)觅食行为:种群中的每个鸟随机进行觅食,随机生成一个数,若其均匀分布在(0,1)之间,则这只鸟就会进行觅食,否则其将保持警戒,其更新位置表达式为: (9) (2)警戒行为:处于警戒行为下的鸟类会向鸟群中心位移,其位置更新公式为: (10) (12) 式中,k(k≠i2)代表从(1,N)中随机选择的正整数;α1与α2是介于(0,2)之间的两个常数;pFiti2表示第i2只鸟的最佳适应度值;sumFit表示鸟群个体适应度值的总和;ε是最小常数;meanj2表示鸟群在j2维的平均位置。 (3)飞行行为:种群飞往另外一个地域时发现者与加入者则分别进行觅食,更新位置公式为: (13) (14) 式中,randn(0,1)表示一个服从标准差为1、均值为0的高斯分布的值;k∈(1,2,…,N),k≠i4;FL(FL∈(0,2))代表加入者将跟从发现者去觅食。 由于鸟群算法中鸟类个体位移的方式是逐步位移,相对于麻雀搜索算法中的跳跃位移可以更加精准地判断每次位移后的位置是否为最佳位置,故将此类更新位置的算法与麻雀搜索算法中发现者与加入者的位置更新公式相结合,从而达到避免出现局部最优的情况。 故根据改进后的麻雀搜索算法将发现者位置更新公式改为: (15) 式中,Xi,j表示此时第i个麻雀在第j维中的位置;Q是一个随机数,服从正态分布;R2代表预警值,范围为R2∈[0,1];Ma代表安全值,范围为Ma∈[0.5,1];公式具体意思参照公式(6)与公式(13)。 将加入者位置更新公式改为: (16) 式中,X2代表加入者加入种群后全局最差位置;k代表加入者,范围为:k∈[1,N];公式具体意思参照公式(7)与公式(14)。 利用ISSA对摄像机进行标定。具体操作过程如下: (1)角点坐标的提取。使用MATLAB中的Camera Calibrator工具箱对已经拍摄好的棋盘格图片进行角点提取[12],如图4所示。 图4 角点的提取 (2)得到初始标定值并确定优化的参数范围[13]。使用MATLAB软件中自带的摄像机标定工具箱,对已经拍摄好的棋盘格图片进行预处理,从而得出四个内部参数fx、fy、u0、v0;外部参数R3×3、T3×1和畸变参数k1、k2、p1、p2。接着将所有参数定为初始标定值,并以此为基准,确定初始值的优化范围。 (3)建立适应度函数。采用ISSA进行摄像机标定,根据公式(1)计算得到像素坐标。而目标函数则是实际的像素坐标与计算得到的像素坐标之间的距离建立的,故其表达式为: (17) 式中,N为图中棋盘格的角点数。 (4)麻雀种群初始化。随机生成n只麻雀,改进麻雀搜索算法中需要使用的参数为m(生产者数量)、p(侦察者数量)、M(安全值)以及Gmax(最大迭代次数)。迭代次数初始值为0。 (5)更新发现者的位置。 (6)更新加入者的位置。 (7)更新侦察者的位置。 (8)获取当前最新的位置。如果当前最新的位置适应度比之前较优,则将之前的位置更新为目前的位置;反之则保留之前的位置。 (9)重复步骤(4)~(8),直到迭代次数达到Gmax。 为了测试改进麻雀搜索算法(ISSA)的精度,用棋盘光刻板作为标定板,搭建了一个简易的实验平台。该实验平台包括一个用作拍摄照片的彩色工业摄像机、8列的光学标定板(方格边长为4 mm)、一个供彩色摄像机与一个计算机。 通过该实验平台的彩色摄像机拍摄20张标定板上棋盘格的图片,其20张图片需要不同位置和倾斜角度。利用MATLAB自带的标定工具箱对所有图片进行预标定,预标定结果如表1与表2所示,其中fx、fy、u0、v0、k1、k2、p1、p2为第一张标定图像所对应的摄像机内参;α、β、θ、tx、ty、tz为第一张标定图像所对应的摄像机外参;n、m、P、M、Gmax为改进麻雀搜索算法的各种初始值。而将这个预标定的值作为ISSA的初始值。 表1 ISSA初始参数(1) 表2 ISSA初始参数(2) 根据ISSA,通过MATLAB的编程,对目标函数f(x)进行迭代寻优,得到摄像机所有参数在此算法下的目标值,如表3所示。 表3 基于ISSA得到的摄像机标定参数 为了探究文中算法的精度与稳定性是否相对于改进前的算法有所提升,现利用基于ISSA的摄像机标定方法、基于SSA算法的摄像机标定方法、基于BAS算法的摄像机标定方法分别对摄像机进行标定。 首先,从采集到的20张棋盘格图片中随机选取10张图片,分别使用SSA和ISSA进行50次重复标定[14]。通过计算这两种标定方法得出的结果的平均值以及平均重投影误差,得出表4与图5。 表4 摄像机参数标定结果 由图5可知,基于ISSA的摄像机标定方法在迭代40次之后逐渐趋于平稳,基于SSA的标定方法则是大约在迭代50次之后趋于平稳,而基于BAS的标定方法大约在迭代60此之后趋于平稳,并且达到最优值。此外,通过图中的信息还可以得出,基于ISSA算法的标定方法所得出的平均重投影误差最小。由表5可知,基于ISSA的算法得出的平均重投影误差相对于基于BAS的算法得出的平均重投影误差小了0.001 2 pixel,相对于基于SSA算法小了0.002 pixel,即所提的ISSA标定方法相较于基础SSA与BAS的标定方法精度更优[15]。 图5 平均重投影误差曲线 表5 ISSA与SSA平均重投影误差 除了利用上面的平均重投影误差曲线可以看出ISSA算法相比于SSA与BAS来说精度更高,还可以通过反投影法来证明此结论的正确性。反投影法是通过计算所得的像素坐标反推世界坐标。这里以反投影法所得出来的世界坐标系中的横坐标与棋盘格实际横坐标差值为x轴,以所得出来的世界坐标系中的纵坐标与棋盘格实际纵坐标为y轴,建立坐标系。从拍摄到的20张棋盘格图片中随机抽取5张,每一张随意取10个点来计算反投影误差,如图6所示。从中可以看出,基于ISSA的标定方法得到的方形点主要分布在原点附近且比较密集,基于BAS的标定方法得到的三角形点分布相对于ISSA比较分散,基于SSA的标定方法得到的圆形点分布最远,这说明基于ISSA的标定方法相比于基于BAS和SSA的标定方法精度更高。 图6 反投影误差散点图 为了验证提出的ISSA的鲁棒性相比于SSA和BAS更好,在拍摄出的20张棋盘格图片上分别加入均值为0,方差分别为0.02、0.04、0.06、0.08和0.1的高斯噪声,分别使用提出的ISSA算法、SSA算法和BAS算法进行50次摄像机标定的实验,以噪声等级为x轴,相对重投影误差为y轴建立如图7所示的坐标系。从图中可以比较清楚发现,随着高斯噪声等级的增加,相对重投影误差逐渐增加,也就是说ISSA、SSA和BAS的标定精度都有所下降,但是改进后的麻雀搜索算法相比于另外两种算法下降幅度更小。即在一定的噪声等级变换范围内,ISSA具有比较好的鲁棒性。 图7 不同噪声等级下标定结果 针对目前一些常用的摄像机标定方法存在精度较低、稳定性差、鲁棒性差等问题,对麻雀搜索算法进行了改进,利用ISSA对摄像机进行标定。主要阐述了摄像机标定原理、SSA原理以及改进麻雀搜索算法(ISSA)的原理,并运用ISSA对摄像机内外参数进行寻优,最后将提出的ISSA算法与SSA算法和BAS算法进行一系列的实验对比,表明基于ISSA的摄像机标定方法的精度以及鲁棒性相对于基于SSA算法、BAS算法的摄像机标定方法都有着一定的提升。3 麻雀搜索算法的改进原理
3.1 鸟群算法(BAS)
3.2 麻雀搜索算法的改进
4 基于ISSA的摄像机参数优化
5 实验与分析结果
5.1 实验环境与参数设置
5.2 精度与稳定性分析
6 结束语