张氏标定法原理及其改进1

张正友算法原理及其改进由于世界坐标系的位置可以任意选取，我们可以假定世界坐标系和摄像机坐标系重合，故定义模板平面落在世界坐标系的Z=0平面上。用r表示R的每一列向量，那么对平面上的每一点，有：=A[r1XWYW01=A[=A[r1XWYW01=A[r1XWYW1(3.1)这样，在模板平面上的点和它的像点之间建立了一个单应性映射H，又称单应性矩阵或投影矩阵。如果已知模板点的空间坐标和图像坐标和M，可以求解单应性矩阵H。(X,Y,1)(u,v,1)那么就已知muXWV=HYW11，其中H=因为swh11h21h31wh12h22h32h13h231，可推出:su=hX+hY+h(3.2)11W12W13<sv=hX+hY+h21W22W23s=hX+hY+131W32W(3.2)故，(3.3)u=h11Xw+VW+h13h31XW02。+1hX+hY+hV=_W2^-W23(3.3)h31XWfl+1将分母乘到等式左边，即有〔VXWh+vY〔VXWh+vYh+v31W32=h21X+hWY+h22W23又令h'=h11h12h13h21h22h23h31h32L，则WW00100XW0YW01-uXW-vXW-uY-W-vYWh'=uv12YW+h13多个对应点的方程叠加起来可以看成Sh=d。利用最小二乘法求解该方程，即h'=(StS)-1STd，进而得到H。|uXh+uYh+u=hX+hW31W3211W(3.4)(3.5)摄像机内部参数求解在求取单应性矩阵后，我们进一步要求得摄像机的内参数。首先令气表示H的每一列向量，需要注意到上述方法求得的H和真正的单应性矩阵之间可能相差一个比例因子，则H可写成：[h1h2h3]=XARrt]

又因为"和七是单位正交向量，所以有hrA-rA-ih2=01hrA-tA-ih=hrA-tA-ih这样就为内参数的求解提供了两个约束方程。下面，令TOC\o"1-5"\h\zBB12"isB=A-rA-i=BBB212223堕i"32"331以2以2Pyv-Pu以2Py(yv0-Pu0)v以2P2P2(yv-Pu以2Py(yv0-Pu0)v以2P2P2(yv-Pu)2v2——0~-~0—+i+i以2P2P2—y2i—————+-以2P以2P2P2yv-Puy(yv-Pu)v

°c0-以2P以2P2P设H的第i列向量为h=[h,h,h]r，因此有iiii2i3hrBh=Vrb(3.7)(3.8)(3.9)(3.10)中V..(3.7)(3.8)(3.9)(3.10)中V..=[hjhh+hh,hh,

iij2i2jii2j2hh+hh,hh+hh,

i3jiiij3i3j2i2j3h.3h」r，那么，就可以将内参数的两个约束写成关于b的两个方程为:Vr(3.11)(3.12)i2b=0Vr-V(3.11)(3.12)如果有〃幅图像的话，把它们的方程式叠加起来，得到Vb=0其中，V是一个2nx6的矩阵。当n>3时，一般情况下，b可以在相差一个尺度因子的意义下唯一确定；当n=2时，此时的方程的个数少于未知数的个数，我们可以加上一个附加约束y=0，即Bi2=0，因此可用[0i0000〕b=0作为式(3.12)的一个附加方程。方程(3.12)的最小二乘解即是VV的最小特征值对应的特征向量，将该向量归一化即得到要求的b，进而得到B；当n=i时，两个方程只能解两个未知数，我们可以假定光心投影在图像的中心，从而求出摄像机在水平和垂直方向上的最大倍数。一旦b被求出，就能根据下面两种方法计算出摄像机的内参数矩阵A:(1)由b构造出B，再利用Cholesky矩阵分解算法求解出A-i，再求逆得到A。(2)由b构造出B,在相差一个尺度因子的意义下(B=人A-rA-i,其中为人尺度因子)，由绝对二次曲线的性质，很容易求出摄像机的内部参数：TOC\o"1-5"\h\zv=(BB—BB)/(BB—B2)012131123112212X=B-[B2+v(BB-BB)]/B221301213112311f=^X7b~x11f=&/B(BB-B2))y'11112212a=-B12f2f/Xu=av/f-Bf2/X00x13x摄像机外部参数求解由每幅图像的单应性矩阵H和上一节的计算结果就可以求得每幅图像的外部参数。一旦A求得后，根据式(3.6)，每幅图像的外部参数很容易求出：r=XA-1hr=XA-1hr=rxrt=XA-1h11223123这里的尺度因子X=1/||A-1hJ|=1/||A-1hJ|o当然，由于图像必然有噪声，因此这样解得的R=[r,r,r]并不能完全满足旋转矩阵的性质，所以要从一个给定123的矩阵求解一个最佳的旋转矩阵。非线性优化(优化内参)以上我们所得到的摄像机的内参数矩阵和每幅图像对应的外参数矩阵都只是一个粗糙解，没有具体的物理意义，可以通过最大似然估计对所有参数进行非线性优化，进一步求精。在这里可以假定有〃幅关于模板平面的图像，模板平面上有m个标定点，那么可建立评价函数：(3.13)C=廿i=1j=1其中m是第i幅图像中的第j个像点，R是第i幅图坐标系的旋转矩阵，t是''!(3.13)第，幅图坐标系的平移向量，M是第j个点的空间坐标，m(A,R,t,M)是通过这些已知量求得的像点坐标。由于旋转矩阵有9个参量但是只有三个自由度，因此可用三个参量的矢量来表示，即一个旋转可由一个三维向量即旋转向量来表示，他的方向就是旋转轴的方向，他的模等于旋转角。由三个欧拉角参数确定，r(r,r,r)r是旋转矩阵的罗德里克(Rodrigues)表123示，R与r之间的关系由公式给出：R=I+sin0/0[r]+(1-cos0)/02[r]2其中，旋转向量r=(x,y,z)r，定义由它构成的反对称矩阵为其中9=32+y2+z2=11rII是旋转角。VJ使评价函数最小的A,R,ti，心j就是这个问题的最优解。这是一个经典的非线性最小二乘问题，对式(3.13)求极小值仍采用Levenberg-Marquardt算法来求解，其初始估计可利用上面线性求解的结果。很明显可以看出，计算顺序依次是投影矩阵、内部参数、外部参数，最后进行优化。对径向畸变处理(u,v)理想像素坐标，(U,而为实际的像素坐标，同样(x,y)和(云y)为理想和实际的图像坐标。x=x+x[k(x2+y2)+k(x2+y2)2]

y=y+y[k](x2+y2)+k2(x2+y2)2]其中，k1,k2为径向畸变，对于中心点畸变同样适用:u=u+(u一u)[k(x2+y2)+k(x2+y2)2]v=v+(v一v)[k(x2+y2)+k(x2+y2)2]可以通过下面的方法求解畸变系数:kk2」(u-u)(x2+y2)(u-u)(x2+y2)2(v-v)(x2+y2)(v-v)(x2+y2kk2」可以利用极大似然估计来得到畸变系数:空m||i=1J=12m一m(A,k,k,R,t,M空m||i=1J=1非线性优化的过程:采用高斯一牛顿方法，而没有用Levenberg一Marquardt(LM)法是因为在实际应用中发现，如果没有先验知识，LM法阻尼因子u不知道去什么值合适，而随意的取值使得方程在求解时收烟速度慢。而高斯一牛顿法在求解的附近一般收敛速度快，虽然对初始值要求严格，当初始值与实际值偏离的太远，其系数阵JTJ容易出现病态解，但是，所以参数的初始值必须要做优化。项(KA,*,八虬)t=tj=]

这里m（K,K,R,T,M）是点M.按照投影公式在第i幅图片上的投影，Kciijj是摄像机标定内参数矩阵，K是畸变系数矩阵，Ri、T.是第i副图像的外参数矩阵，其中R通过Rodrigues公式用三个参数的旋转向量r来表达优化的步骤：（1）、首先确定目标函数。这里的优化的目的是使优化后的外参数，按照本文建立的摄像机模型将空间点投影到图像平面上时，投影点与实际获取的图像控制点之间的残差平方和最小。不妨设投影函数为"侦，3：匚虬）其中k是摄像机内参数矩阵，"Tuid是畸变系数矩阵，丁是旋转向量，T是平移矩阵，"二1」一⑴是模板平面上控制点Mi的坐标。"-是投影点图像像素坐标，该点实际的图像坐标是\共n个控制点。要优化外参数，则最后需要的偏导数系数Jacobian矩阵dm./dr和dm./dT得组合，而函数、））显然是个复合函数，因此这实际上是个复合函数求偏导的问题。下面针对第一次迭代即r=r0、T=T0为初值时的一个点..广。1来分步求解第一步，是r与R的转换。由Rodrigues公式可以求出I-cos^?r「PT&2&=Rodrigues^）和。&/粉",因为］＜有9个元而r有3个元，显然两伽"是个9号的矩阵，如（4一32、(4.32)I-cos^?r「PT&2(4.32)(4.33)(434)

艾小奇..芸.卜：55这里偏导数矩阵中的下标i表示求出偏导数矩阵后，再将当前点（即第i点）各相应坐标代入以求得偏导数系数矩阵。设机=区,％力）七m=（y，tz）「，则很容易求得就仞…电瓦=流挣'(434)艾小奇..芸.卜：55这里偏导数矩阵中的下标i表示求出偏导数矩阵后，再将当前点（即第i点）各相应坐标代入以求得偏导数系数矩阵。设机=区,％力）七m=（y，tz）「，则很容易求得就仞…电瓦=流挣'…电％按洌-记依「％0。匕0oZ前。=°n00匕00孔〔0°Y00么0D(4拓)因为虬有3个元而R有9个元,00ZHJJ11是…次函数，所以匕L弋的结果是很明祢J同时可以看出，若直接以R为参量，因为复合函数求偏导是偏导矩阵相乘的关系，则由于是3*9矩阵而'村.『西是3、；3中；'：。第三步，将厉Mci转化为归一化图像坐标。设M,的归一化图像坐标为M可得Mni的计算公式及其偏导数系数阵：ci第三步，将厉Mci转化为归一化图像坐标。设M,的归一化图像坐标为M可得Mni的计算公式及其偏导数系数阵：cici•，nix_XCi(4.36)商dr(1=20注意：(4、37)bi\''(8)/=i,,b.■-IL勺】表示qxp系数矩阵B将与矩阵A中的所有分块矩阵七（i=1,2,,n,j=1,2,,m）均相乘，其中分块矩阵A^（i=1,2,,n,j=1,2,,m）的行数等于系数矩阵B的列数p。

第四步，计算考虑有图像畸变的归一化图像坐标。设「•《'」..'.『,’，则含有畸变的归一化图像坐标和偏导数系数阵可由归一化公式（4.38）式来计算（4.39）（4.40）：「七1M祁-.438）「孩］|~又［11■"卜顷乂尸］义2前』-乩（3丸\疔）I_=（439）顷」［用I+a'+馈如馈'］++3丸’）+%戒"「&丽』二仲"「如源』-泠ST」|gE［腕儿*37liI（4.40）=«叮亟亟］［c印3T」Ln并且可以算得其中八［（1+标尸+左,+灼汽+；（2*+4顼+庶尸）+（2内下+6成）|j卜氏丸（2耳+"+6垢，）+（2p£+6a丸）Lc=［xiiyn（2k1+4"J6上尸）+（2/低十6夕疣）d=［（1+佑尸+顷+"与+衬（2肉+4侦'+6"』）+（%丸十2p成）.对畸变系数的优化:dmnd=dKcxdmnd=dKcxr2n—一2yr2nxnr4yr4jl,⑷心)8mdmndfajl,⑷心)8mdmndfau0fV8m8m8maT=(~m]cndc8m=8m=xnddK~[00ynd0100ynd1'到此我们得到了第一次送代求得扇像素点坐标和迭代增届公式中所需的Jacobiinitial即优化外部参数：山』里^11L济泣」卜=仆5优化全局参数：j=[8m8m8m8mi8K'8K'8r'8T1，r0，T0第六步，求第一次迭代产生的偏差，并由式:网三{尸J)\尸虬{444)A6=(JtJ)-1Jt&i计算参数迭代增量，为下一次迭代做准备。q=m,.-mj/_IBP=八（4.46）或

LJ然.希从第•步升始，用参数初始值/*+彼、咋、本+如替代gL进仃第：次迭代,人断循沼，百到得到合适的参量值后结束需要注明的是对巳44)称为高斯-牛顿法.若0-(J『J顼)*0则为Levenberg-Marquardvi改进的方法具体介绍(1)、改进的平面模板两步法标定摄像机一一毛剑飞，邹细勇，诸静该方法张正友提出用平面模板两步法来标定摄像机。该方法能高精度地标定摄像机，且简便易行，可以说是两步法标定的代表，但其算法模型仅考虑了镜头径向畸变，其对切向畸变较大的场合，如鱼眼镜头成像，则不适用。为此，通过改进其摄像机模型，提出了一种改进的两步法，该方法先用图像中心附近点求取初值，由于图像中心附近点畸变很小，故求取的初值能很好地逼近准确值，然后采用一种基于内部映射牛顿法的子空间置信域法通过求精来得到所有参数。由于两步法要求一个比较好的初值，若初值选择不当，则算法难以收敛或只能收敛到局部最小，从而大大降低标定精度。zhang的方法缺点在求取初值中，虽先不考虑各种畸变，而是先将所有点代入求解，但由于远离图像中心的像点畸变很大，如将这些点也看作没有畸变的像点代入，显然会加大求解初值的误差，但由于考虑到图像中心附近点的畸变很小，因此可以先利用图像中心附近点求取初值。本文算法正是基于这点考虑的，由于本文算法准确地计算了初值，因此接下来的全面考虑各参数的非线性最小化计算就能很快收敛。仿真摄像机的特性：分辨率为768X576像素，真实实验使用三星SCC-421P型彩色摄像机，分辨率为576X768,模板是用激光打印机打印7X10的国际象棋图案，并贴在硬塑料板上制成的标定模板，每格边长2.5000cm。因为整个包括角点检测的标定程序已经编好，所以只要让模板对着摄像机随意转动几次，即可完成标定，时