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

张正友算法原理及其改进由于世界坐标系的位置可以任意选取，我们可以假定世界坐标系和摄像机坐标系重合，故定义模板平面落在世界坐标系的0WZ平面上。用ir表示R的每一列向量，那么对平面上的每一点，有：123120111(3.1)这样，在模板平面上的点和它的像点之间建立了一个单应性映射H，又称单应性矩阵或投影矩阵。如果已知模板点的空间坐标和图像坐标，那么就已知m和M，可以求解单应性矩阵H。)1,,(wwYX)1,,(vu因为11WWuXsvHY,其中11121321222331321hhhHhhhhh,可推出：11121321222331321(3.2)故，1112133132212223313211(3.3)将分母乘到等式左边，即有31321112133132212223(3.4)又令Thhhhhhhhh3231232221131211，则10000001(3.5)多个对应点的方程叠加起来可以看成Shd。利用最小二乘法求解该方程，即1()TThSSSd，进而得到H。摄像机内部参数求解在求取单应性矩阵后，我们进一步要求得摄像机的内参数。首先令ih表示H的每一列向量，需要注意到上述方法求得的H和真正的单应性矩阵之间可能相差一个比例因子，则H可写成：12312hhhArrt(3.6)又因为1r和2r是单位正交向量，所以有1210TThAAh(3.7)111122TTTThAAhhAAh(3.8)这样就为内参数的求解提供了两个约束方程。下面，令1112131212223313233002222000222222222000000002222222111TBBBBAABBBBBBvuvuvvuvuvuvv注意到B是一个对称矩阵，所以它可以由一个6维向量来定义，即111222132333TbBBBBBB(3.9)设H的第i列向量为123[,,]Tiiiihhhh，因此有TTijijhBhVb(3.10)其中111221223113322333,,,,,TijijijijijijijijijijVhhhhhhhhhhhhhhhhhh，那么，就可以将内参数的两个约束写成关于b的两个方程为：1211220TTTVbVV(3.11)如果有n幅图像的话，把它们的方程式叠加起来，得到0Vb(3.12)其中，V是一个26n的矩阵。当3n时，一般情况下，b可以在相差一个尺度因子的意义下唯一确定；当2n时，此时的方程的个数少于未知数的个数，我们可以加上一个附加约束0，即120B，因此可用0100000b作为式(3.12)的一个附加方程。方程(3.12)的最小二乘解即是TVV的最小特征值对应的特征向量，将该向量归一化即得到要求的b，进而得到B；当1n时，两个方程只能解两个未知数，我们可以假定光心投影在图像的中心，从而求出摄像机在水平和垂直方向上的最大倍数。一旦b被求出，就能根据下面两种方法计算出摄像机的内参数矩阵A：(1)由b构造出B，再利用Cholesky矩阵分解算法求解出1A，再求逆得到A。(2)由b构造出B，在相差一个尺度因子的意义下(1TBAA,其中为尺度因子),由绝对二次曲线的性质，很容易求出摄像机的内部参数：2012131123112212()/()vBBBBBBB2221301213112311[()]/BBvBBBBB11/xfB211112212/()yfBBBB212/xyBff20013//xxuvfBf摄像机外部参数求解由每幅图像的单应性矩阵H和上一节的计算结果就可以求得每幅图像的外部参数。一旦A求得后，根据式(3.6)，每幅图像的外部参数很容易求出：111rAh122rAh312rrr13tAh这里的尺度因子11121/1/AhAh。当然，由于图像必然有噪声，因此这样解得的123[,,]Rrrr并不能完全满足旋转矩阵的性质，所以要从一个给定的矩阵求解一个最佳的旋转矩阵。非线性优化（优化内参）以上我们所得到的摄像机的内参数矩阵和每幅图像对应的外参数矩阵都只是一个粗糙解，没有具体的物理意义，可以通过最大似然估计对所有参数进行非线性优化，进一步求精。在这里可以假定有n幅关于模板平面的图像，模板平面上有m个标定点，那么可建立评价函数：211,,,nmijijijCmmARtM(3.13)其中ijm是第i幅图像中的第j个像点，iR是第i幅图坐标系的旋转矩阵，it是第i幅图坐标系的平移向量，jM是第j个点的空间坐标，,,,ijmARtM是通过这些已知量求得的像点坐标。由于旋转矩阵有9个参量但是只有三个自由度，因此可用三个参量的矢量来表示，即一个旋转可由一个三维向量即旋转向量来表示，他的方向就是旋转轴的方向，他的模等于旋转角。由三个欧拉角参数确定，123(,,)Trrrr是旋转矩阵的罗德里克(Rodrigues)表示，R与r之间的关系由公式给出：22sin/[](1cos)/[]xxRIrr(3.14)其中，旋转向量Tzyxr),,(，定义由它构成的反对称矩阵为000][xyxzyzrx，其中||||222rzyx是旋转角。使评价函数最小的,,,ijARtM就是这个问题的最优解。这是一个经典的非线性最小二乘问题，对式(3.13)求极小值仍采用Levenberg-Marquardt算法来求解，其初始估计可利用上面线性求解的结果。很明显可以看出，计算顺序依次是投影矩阵、内部参数、外部参数，最后进行优化。对径向畸变处理(u,v)理想像素坐标，),(vu为实际的像素坐标，同样(x,y)和),(yx为理想和实际的图像坐标。])()([])()([22222212222221yxkyxkyyyyxkyxkxxx其中，21,kk为径向畸变，对于中心点畸变同样适用：])()()[(])()()[(2222221022222210yxkyxkvvvvyxkyxkuuuu可以通过下面的方法求解畸变系数：vvuukkyxvvyxvvyxuuyxuu2122202202220220))(())(())(()()(高斯一牛顿方法优化：可以利用极大似然估计来得到畸变系数：21121||),,,,,(||nimjjiiijMtRkkAmm非线性优化的过程：采用高斯一牛顿方法，而没有用Levenberg一Marquardt(LM)法是因为在实际应用中发现，如果没有先验知识，LM法阻尼因子u不知道去什么值合适，而随意的取值使得方程在求解时收烟速度慢。而高斯一牛顿法在求解的附近一般收敛速度快，虽然对初始值要求严格，当初始值与实际值偏离的太远，其系数阵JTJ容易出现病态解，但是，所以参数的初始值必须要做优化。这里),,,,(jiicMTRKKm是点jM按照投影公式在第i幅图片上的投影，K是摄像机标定内参数矩阵，cK是畸变系数矩阵，iR、iT是第i副图像的外参数矩阵，其中R通过Rodrigues公式用三个参数的旋转向量r来表达。优化的步骤：(1)、首先确定目标函数。这里的优化的目的是使优化后的外参数，按照本文建立的摄像机模型将空间点投影到图像平面上时，投影点与实际获取的图像控制点之间的残差平方和最小。不妨设投影函数为其中K是摄像机内参数矩阵，是畸变系数矩阵，是旋转向量，T是平移矩阵，是模板平面上控制点iM的坐标。是投影点图像像素坐标，该点实际的图像坐标是，共n个控制点。要优化外参数，则最后需要的偏导数系数Jacobian矩阵rmi/和Tmi/得组合，而函数)显然是个复合函数，因此这实际上是个复合函数求偏导的问题。下面针对第一次迭代即r=r0、T=T0为初值时的一个点来分步求解第一步，是r与R的转换。由Rodrigues公式可以求出.第二步，是世界坐标系到摄像机坐标系的转换。由坐标转换公式求出摄像机坐标系下的坐标及相应偏导数系数阵：这里偏导数矩阵中的下标i表示求出偏导数矩阵后，再将当前点(即第i点)各相应坐标代入以求得偏导数系数矩阵。同时可以看出，若直接以R为参量，因为复合函数求偏导是偏导矩阵相乘的关系，则由于是3*9矩阵而。第三步，将厉ciM转化为归一化图像坐标。设ciM的归一化图像坐标为niM，可得niM的计算公式及其偏导数系数阵：注意：表示qxp系数矩阵B将与矩阵A中的所有分块矩阵ijA(i=1,2,,n,j=1,2,,m)均相乘，其中分块矩阵ijA(i=1,2,,n,j=1,2,,m)的行数等于系数矩阵B的列数p。第四步，计算考虑有图像畸变的归一化图像坐标。设，则含有畸变的归一化图像坐标和偏导数系数阵可由归一化公式(4.38)式来计算(4.39)(4.40)：对畸变系数的优化：4242ryryrxrxKmnnnncnd第五步，求迭代计算出的投影点图像像素坐标和偏导数系数阵，按照(4.18)式进行计算：vundffmm0])[(cndndcKmmmKm若是对内部参数的优化为：1000001ndndndyyxKm优化外部参数：优化全局参数：00,,|],,,[TriciTmrmKmKmJ第六步，求第一次迭代产生的偏差，并由式：iTTJJJ1)(计算参数迭代增量，为下一次迭代做准备。iiimm改进的方法具体介绍(1)、改进的平面模板两步法标定摄像机——毛剑飞，邹细勇，诸静该方法张正友提出用平面模板两步法来标定摄像机。该方法能高精度地标定摄像机，且简便易行，可以说是两步法标定的代表，但其算法模型仅考虑了镜头径向畸变，其对切向畸变较大的场合，如鱼眼镜头成像，则不适用。为此，通过改进其摄像机模型，提出了一种改进的两步法，该方法先用图像中心附近点求取初值，由于图像中心附近点畸变很小，故求取的初值能很好地逼近准确值，然后采用一种基于内部映射牛顿法的子空间置信域法通过求精来得到所有参数。由于两步法要求一个比较好的初值，若初值选择不当，则算法难以收敛或只能收敛到局部最小，从而大大降低标定精度。zhang的方法缺点在求取初值中，虽先不考虑各种畸变，而是先将所有点代入求解，但由于远离图像中心的像点畸变很大，如将这些点也看作没有畸变的像点代入，显然会加大求解初值的误差，但由于考虑到图像中心附近点的畸变很小，因此可以先利用图像中心附近点求取初值。本文算法正是基于这点考虑的，由于本文算法准确地计算了初值，因此接下来的全面考虑各参数的非线性最小化计算就能很快收敛。仿真摄像机的特性：分辨率为768×576像素，真实实验使用三星SCC-421P型彩色摄像机，分辨率为576×768，模板是用激光打印机打印7×10的国际象棋图案，并贴在硬塑料板上制成的标定模板，每格边长2.5000cm。因为整个包括角点检测的标定程序已经编好，所以只要让模板对着摄像机随意转动几次，即可完成标定，时间不到1