基于改进风驱动算法的摄像机标定优化方法

任久斌，曹中清

(西南交通大学 机械工程学院，四川 成都 610031)

0 引 言

现有的摄像机标定方法主要有直接线性变换法、Tsai两步标定法和张正友标定法[1-4]。由于张正友标定法使用的是平面标定板，操作简单，成本低廉，所以它在实际摄像机标定中被广泛应用。但因其采用的非线性优化算法具有局限性，常常陷入局部最优解，所以得到的摄像机标定结果精度较差。近年来，许多研究者提出了基于遗传算法(GA)、模拟退火算法(SA)、粒子群算法(PSO)等的改进方法，尽管这些算法在一定程度上提升了标定精度，但是其改进效果并不明显[5-8]。

风驱动优化算法是近年来提出的一种全局优化算法，它具有开发能力强、收敛速度快、寻优效率高的特点，非常适合于求解多维多模态的问题[9-11]。然而基本风驱动算法中固有参数的人为设置非常困难，不同的设置值通常会带来差异性的结果。基于此，本文提出基于协方差矩阵自适应进化策略的改进算法，通过利用协方差矩阵中不同维度的约束关系，实现参数采样空间的选择和重组，确保参数的采样分布朝着适应度值减小的方向移动。在获得自适应选择的固有参数后，再利用基本风驱动算法优化摄像机内参。实验结果表明，改进算法的收敛速度、收敛精度以及稳定性都优于其它算法。

1 摄像机成像模型

为了表达场景中的目标点与其投影点之间的关系，需要首先建立摄像机的成像模型[3]。摄像机成像模型中最常用的是线性模型(即针孔模型)，如图1所示。图中Ow-XwYwZw表示世界坐标系，OC-XCYCZC表示摄像机坐标系，O0-uv表示图像像素坐标系，O1-xy表示图像物理坐标系。

图1 摄像机线性模型

假设空间点P在世界坐标系中的坐标为(Xw,Yw,Zw)，其在图像像素坐标系中的投影点p的坐标为(u,v)，它们之间的关系可以表示为

(1)

式中：αx=f/dx，αy=f/dy，它们分别表示摄像机在x方向和y方向上以像素为单位的焦距，γ为坐标轴的不垂直度系数；M1为3×3的矩阵，由摄像机的内部参数αx、αy、u0、v0、γ构成，称为摄像机的内参矩阵；M2为3×4的矩阵，由旋转矩阵R和平移向量t构成，描述了摄像机坐标系相对于世界坐标系的关系，称为摄像机的外参矩阵；M为3×4的矩阵，即内参矩阵和外参矩阵的乘积，称为摄像机的投影矩阵。

(2)

其中，δx和δy分别表示图像中某一点在x方向和y方向的畸变值。它们的形式如式所示

(3)

式(3)中3项分别是径向畸变、切向畸变和薄棱镜畸变，p1、p2、k1、k2、s1、s2为非线性畸变系数。本文只考虑径向畸变和切向畸变的前两项，即k1、k2、p1、p2，因为这两项足以描述高精度测量中的透镜畸变。

2 摄像机模型参数初值的求解

根据张正友的摄像机标定方法，为了得到高精度的标定结果，需要首先为优化算法计算提供内参初值。通过将世界坐标系固定在平面标定板上，令标定板平面为Zw=0的平面，则上述摄像机的成像模型可以简化为

(4)

其中，r1和r2为旋转矩阵R的第一列和第二列。再利用图像平面和标定板平面间的单应性关系可得

(5)

(6)

其中，λ为尺度因子。由于旋转矩阵的两个列向量r1和r2正交，于是有

(7)

然后再将式(6)和式(7)结合，整理变换后可得如下两个约束

(8)

由于摄像机内参矩阵有5个参数，而每个单应性矩阵只能够提供两个约束，因此至少需要3个不同的单应性矩阵才能求解，这里可以通过摄像机采集至少3幅不同方位的标定板图像实现[8]。

为了便于计算，令

(9)

由式(9)可以看到B是一个对称矩阵，它可以用一个6维的向量表示为

(10)

hiTBhj=vijTb

(11)

其中，vij为

(12)

由式(8)的约束关系可得

(13)

在采集到标定板的n(n≥3)幅不同图像后，利用最小二乘法可以求解出b，进而求解得到摄像机内参的初始值。在得到单应性矩阵和内参矩阵的基础上，还可以求得摄像机的外部参数

(14)

3 基于改进风驱动算法的摄像机参数优化

3.1 基本风驱动优化算法

风驱动优化算法(wind driven optimization，WDO)是由Bayraktar提出的一种基于群体的启发式迭代搜索算法，在求解全局优化问题时具有简单易行和高效的特点。Bayraktar首先通过对地球表面对流层空气运动状况的分析，抽象出简化的大气模型，然后结合牛顿第二定律和理想气体状态方程，推导出WDO算法中空气微团的速度更新公式和位置更新公式[11]。

对某一单位质量的空气微团受力分析，由牛顿第二定律可得

ρa=∑Fi

(15)

式中：ρ为空气的密度；a为空气微团运动的加速度，其大小为a=Δu/Δt，Δu表示Δt时间内的速度变化；Fi为作用在空气微团上的力。

由空气动力学知识，式(15)中主要包含了4种主要的作用力，使得空气微团以一定速度向某个方向移动或者使它偏离现有路径，这4种力分别是：气压梯度力(FPG)、摩擦力(FF)、重力(FG)和科里奥利力(FC)。利用它们的简化表示，式(15)可以变换为

(16)

气压梯度的变化是因为空气温度差异引起的，为了将空气微团的密度和温度联系起来，可以利用如下所示的理想气体状态方程

P=ρRT

(17)

其中，P为大气压强，R为理想气体常数，T为气体的温度。

为了简化式(16)，令δV=1， Δt=1，并将式(17)代入式(16)后得

(18)

式中：Pcur表示当前位置空气微团的压力值；Δu=unew-ucur，ucur表示空气微团当前迭代的速度，unew表示空气微团下一次迭代的速度。

气压梯度是使空气微团从当前位置xcur移向最优位置xopt的原因，其大小为|Popt-Pcur|，因此气压梯度可表示为

(19)

其中，Popt表示最优位置的压力值。将重力加速度矢量g改写为g=|g|(0-xcur)，代入到式(18)中得

(20)

(21)

其中，c=-2|Ω|RT。式(21)中使用了实际的压力值，当压力值非常大的时候，空气微团更新的速度也会变得不切实际的大，从而削弱了WDO的搜索性能。取而代之的是使用基于压力值排序的方法，即通过对所有空气微团的压力值按升序排序，用排序后当前位置空气微团的序号代替上述公式中的压力值，这样式(21)可以改写为

(22)

式中：i代表当前位置xcur空气微团压力值Pcur排序后的序号，最优位置xopt处的最优压力值Popt排序序号始终为1。

在得到空气微团的速度更新方程以后，根据式(23)可以得出空气微团的位置更新方程

xnew=xcur+(unewΔt)

(23)

式中：通常假设Δt=1。

式(22)和式(23)即为最终得到的基本WDO算法的速度更新方程和位置更新方程。

3.2 协方差矩阵自适应进化策略

基本WDO算法中包含了4个参数α、g、RT和c，不同的优化问题通常需要将它们设置为不同的值，而且在整个优化迭代过程中它们固定不变[11]。为了减轻手动设置的负担，避免人为选择参数带来的不确定性，本文提出基于协方差矩阵自适应进化策略(CMAES)自动地选择WDO算法的4个参数。

CMAES将α、g、RT和c作为待优化参数，群体的采样分布用协方差矩阵描述，通过自适应地改变采样分布的位置和形状，引导搜索群体朝着适应度值减小的方向进化，直到达到最大迭代次数为止[13-15]。CMAES算法的具体操作步骤如下：

(1)参数初始化。设置群体的大小λ，重组操作的子群大小μ=λ/2，重组权值wi(i=1,2,…,μ)以及自适应调整所需的常量cσ、μeff、dσ、cc、c1见表1，设置算法的最大迭代次数G和问题的维度N(N=4，即待优化的4个参数)。在搜索空间中随机选择采样分布的初始均值m(0)，采样初始步长σ(0)设置为最大参数范围的1/3，设置初始步长进化路径Pσ(0)和初始协方差矩阵进化路径Pc(0)为0，将初始协方差矩阵C(0)设置为单位矩阵I，进化代数g的初值设置为0。

表1 CMAES内部常量参数的计算公式

(2)采样操作。按照下式从采样分布空间中抽取λ个个体形成搜索种体

(24)

(25)

(3)选择和重组操作。计算群体中所有个体的适应度值(即空气微团的压力值)，并对其进行升序排序，选择适应度较小的前μ个个体形成子群，然后对子群个体加权求和更新采样分布的均值

(26)

(4)采样步长调整

(27)

(28)

其中，协方差矩阵的逆(C(g))-1/2=B(g)(D(g))-1(B(g))T。

(5)协方差矩阵调整，引导群体搜索全局最优解。先对协方差矩阵进化路径Pc进行更新

(29)

Pc的进化方向保证了采样分布沿均值更新方向移动。协方差矩阵自适应更新公式为

(30)

其中，第三项利用了适应度较好的μ个个体，因此采样分布朝着适应度减小的方向进化[15]。

(6)若g达到最大迭代次数G，则退出循环；否则，g=g+1，并跳转到步骤(2)继续执行。

3.3 摄像机内参优化过程

为了避免过多的参数优化导致算法寻优能力下降，本文采用外参固定的优化模型，在第2节计算的初值基础之上，优化摄像机的内部参数。假设采集了n幅标定板图像，每幅图像有m个特征点，定义优化的目标函数(即重投影误差)如下所示

(31)

基于改进风驱动算法的摄像机内参优化过程如下：

(1)算法初始化。设置空气微团种群的大小M，优化问题的维度D(即待优化的9个摄像机内参：αx、αy、u0、v0、γ、k1、k2、p1、p2)，最大进化代数Gen，将当前进化代数gen初始化为0。在摄像机内参初值的基础上产生空气微团的位置种群和速度种群，并为每个空气微团随机初始化一个参数(α、g、RT和c)，所有参数组成参数种群；

(2)利用目标函数(31)评估初始空气微团种群的适应度值，并对它们进行升序排序，记录下适应度值最小的个体gbest；

(3)根据参数种群提供的WDO算法参数，由式(22)和式(23)更新空气微团种群中每个个体的速度和位置，再次利用式(31)计算种群更新后的适应度值并排序，更新全局最优个体gbest；

(4)利用CMAES自适应更新参数种群，具体实现如3.2节所示，其中参数种群的大小λ和空气微团种群大小M相同，维度N=4；

(5)令gen=gen+1，若gen

4 实验与结果分析

为了验证本文改进算法的性能，下面进行了摄像机标定实验。首先，用摄像机从不同方位拍摄棋盘格标定板(9×12，每个方格大小为10 mm)的20幅图像，由于世界坐标系固定在标定板上(Zw=0)，因此可以很方便地计算出标定板上每个角点的世界坐标，而对于标定板图像中角点的像素坐标，本文利用Harris角点检测算法提取。然后再根据第2节所描述的方法计算出摄像机内参的初值，所得结果见表2。

表2 摄像机内参初值

在得到各个参数的初值后，利用下面的方式对它们进行简单地扩展，以此创建每个参数的搜索边界：[αxαyu0v0]±100，γ∈[-1,1]， [k1k2p1p2]∈[-2,2]。实验中，本文设置的种群大小为50，最大迭代次数为100。

接下来分别利用张正友标定法(Zhang)、粒子群算法(PSO)、基本风驱动算法(WDO)和本文改进的风驱动算法(CMAES-WDO)进行实验，各算法分别运行50次，将各次所得的参数优化结果和重投影误差值取平均，结果见表 3。

从表3中可以看出，张正友标定法的重投影误差均值为0.217 69，经过PSO算法优化的标定重投影误差均值为0.179 32，经过基本WDO算法优化的标定重投影误差均值为0.158 12，而经过本文改进的风驱动算法的标定重投影误差均值仅为0.079 61。多次实验的结果表明，尽管经过基本WDO算法优化标定的摄像机内参精度有所改善，但因为是人为选择的算法参数，所以改进的程度并不高。

表3 各算法的优化结果

图2～图5中展示了4种方法的标定重投影误差分布。从图2和图3的对比可以看出，PSO优化标定的重投影误差分布相较于Zhang标定法的重投影误差分布更加密集，说明PSO算法在一定程度上能够提高摄像机的标定精度。从图3和图4的对比可以看出，基本WDO算法相较于PSO算法的性能更优，但标定精度的改善程度并不大。从图5中可以看出，经过CMAES-WDO算法优化标定后的重投影误差在x方向和y方向均集中在[-0.15,0.15]之间，和其它3种方法相比其分布密度最小，所以它最终得到的摄像机标定参数精度最高。图6中展示了PSO算法、基本WDO算法和CMAES-WDO算法用于摄像机标定优化时目标值的变化过程，可以发现CMAES-WDO算法迭代约10次就可以完全收敛，而且收敛精度约为0.08，相较于其它两种算法，CMAES-WDO算法的收敛速度和收敛精度都更优。

图2 Zhang标定法的重投影误差分布

图3 PSO算法的重投影误差分布

图4 WDO算法的重投影误差分布

图5 CMAES-WDO算法的重投影误差分布

图6 目标值优化过程

5 结束语

本文利用协方差矩阵自适应进化策略自动地确定了基本风驱动算法中的固有参数，避免了因为人为参数选择带来的差异性结果，并将改进后的算法用于摄像机参数的标定优化中。实验结果表明，改进风驱动算法相较于基本风驱动算法的性能更好，摄像机参数标定优化后的精度更高，得到的摄像机内参结果更准确。通过与粒子群算法和基本风驱动算法对比，改进算法的收敛速度更快，收敛精度更高。此外，多次重复实验的结果表明，改进算法的稳定性也更好。