专利名称：一种多摄像机阵列的联合标定方法
技术领域：
本发明属于图像处理技术领域，具体涉及一种属于摄像机标定领域的不同分辨率摄像机阵列的标定方法。
背景技术：
使用多摄像机阵列系统进行测量，必须首先进行摄像机阵列的标定，只有获得了鲁棒和高精度的摄像机参数值，才能保证整个系统的测量精度。标定所需得到的参数一般为各摄像机的内部参数以及摄像机之间的转化关系。
传统的多摄像机标定方法一般需要电子经纬仪等设备，以用于建立世界总体坐标系，其缺点是增加了系统的成本和复杂性。在计算各摄像机间转换关系时，现有方法多是将单摄像头标定的数据直接计算得到，这种计算方法会存在参数一致性的问题。而当系统中含有不同分辨率的摄像机时，这种方法的系统标定精度受低分辨率摄像机影响，标定精度较低，且不能保证鲁棒性。
发明内容
为克服现有技术的上述缺陷，本发明的目的在于提供一种多摄像机阵列系统的联合标定方法，这种方法的标定精度高，保证了鲁棒性，且适用于具有不同分辨率的摄像机阵列系统。
本发明的设计思路如下由于多摄像机阵列系统中的摄像机拥有不同的分辨率，故本发明人把阵列系统视作一个整体进行考虑，以系统优化的思路对阵列中的各摄像机进行联合优化，解决参数的一致性问题，提高系统的鲁棒性，达到阵列总体误差最小意义下的最优标定。另外，由于低分辨率摄像机会对系统整体精度的产生影响，本发明还采用高分辨摄像机与低分辨率摄像机分开处理的策略。即先标定好高分辨率的摄像机，然后对高分辨率摄像机和低分辨率摄像机进行联合优化，利用高分辨率摄像机的结果提高低分辨率摄像机的参数精度。通过使用本方法，阵列系统标定的整体误差大大减小，各参数的鲁棒性显著提高。
为此，本发明拟采用的技术方案如下 一种多摄像机阵列的联合标定方法，它由不同分辨率的摄像机组成多摄像机阵列系统，假设共有w个高分辨率摄像机，各摄像机对一个含有m个特征点的平面模板同步地拍摄n幅图像，进行单摄像机标定，获得了各自的内、外参数，外参数共有n组，分别对应拍摄的n幅图像，并且已获得特征点在各摄像机成像平面上像素坐标的测量值。该联合标定方法包括如下步骤 步骤1、对高分辨率摄像机两两间进行标定，以得到各高分辨率摄像机的内参数初值和各高分辨率摄像机间转换参数的初值； 步骤2、对各高分辨率摄像机进行联合标定，以得到各高分辨率摄像机的内参数和各高分辨率摄像机间转换参数的最终值； 步骤3、对高分辨率摄像机与低分辨率摄像机进行联合标定，以得到低分辨率摄像机内参数和低分辨率摄像机与各高分辨率摄像机间转换参数的最终值。
所述步骤1的对高分辨率摄像机两两间进行标定包括如下步骤 步骤11、对各高分辨率摄像机进行编号，记为摄像机1、摄像机2、…、摄像机w，w为摄像机数； 步骤12、任取一对高分辨率摄像机，将编号小的摄像机作为左摄像机，将编号大的摄像机作为右摄像机； 步骤13、将预先通过单摄像机标定获得的内、外参数作为左摄像机、右摄像机内、外参数的初值； 步骤14、对左、右摄像机的参数进行非线性最小二乘优化； 步骤15、对每一对摄像机重复执行步骤12至步骤14； 步骤16、对得到的参数进行鲁棒性处理； 经过上述步骤的优化，得到每个高分辨率摄像机的内参数的初值，摄像机1的外参数的初值，摄像机1到其他高分辨率摄像机的转换矢量的初值。
所述步骤14的对左、右摄像机的参数进行非线性最小二乘优化包括如下步骤 步骤141、由左、右摄像机外参数的初值，计算其转换关系，作为左摄像机到右摄像机的转换参数的初值； 步骤142、由左、右摄像机的内、外参数的初值，计算模板特征点在左、右摄像机成像平面上的像素坐标，并将此计算值与特征点在左右摄像机成像平面上像素坐标的测量值相减求差； 步骤143、对左、右摄像机内参数和左摄像机外参数，以及步骤141得到的左摄像机到右摄像机的转换参数，进行非线性最小二乘优化。
所述步骤141包括如下步骤 步骤1411、将左、右摄像机的外参数分别用Rl、tl和Rr、tr表示，其中Rl、Rr表示旋转矩阵，tl、tr表示平移矩阵，则左摄像机到右摄像机的转换参数即旋转矩阵和平移矩阵Rl2r、tl2r为 tl2r＝tr-Rl2rtl (1) 步骤1412、对拍摄的n幅图像，即n组外参数，重复步骤1411； 步骤1413、对步骤1412得到的n组Rl2r、tl2r参数进行中值滤波，将结果作为左摄像机到右摄像机的转换参数的初值。
所述步骤142包括如下步骤 步骤1421、由左、右摄像机的外参数，计算特征点在左、右摄像机坐标系下的坐标，方法如下 设XX＝[X；Y；Z]代表点P在世界坐标系中的坐标，XXC＝[XC；YC；ZC]代表点P在摄像机坐标系下的坐标，那么XX与XXC的关系可以用如下刚体运动转化来表示 XXc＝rodrigues(om)*XX+T(2) 其中，rodrigues(.)表示对括号内容输入Rodrigues公式； om、T分别表示摄像机外参数中的旋转矢量和平移矢量； 步骤1422、计算特征点归一化针孔投影的坐标Xn 步骤1423、计算考虑镜头畸变影响后的归一化点坐标xd 其中，畸变包括径向畸变与偏心畸变的切向分量，用5×1的矢量kc表示，其中k(1)、k(2)、k(5)分别是径向畸变的二次、四次和六次项，k(3)和k(4)分别是两个切向畸变系数； 步骤1424、计算特征点P在图像平面上最终的像素坐标pixel＝[xp，yp] 其中KK表示摄像机的内参数矩阵 其中f、c、alpha分别表示焦距、主点坐标、倾斜因子； 步骤1425、将特征点在左右摄像机成像平面上像素坐标的计算值与测量值相减求差。
所述步骤143包括 对左、右摄像机内参数和左摄像机外参数，以及左摄像机到右摄像机的转换参数的初值进行非线性最小二乘优化，使如下代价函数 最小，即 其中，param为需要优化的摄像机参数 param＝(fl；cl；alphal；kl；fr；cr；alphar；kr；oml2r；Tl2r；oml1；Tl1；...；omln；Tln)；
为第j个特征点Mj在在第k个摄像机的第i幅图像上的投影的测量值，m(param，Mj，k，i)为相应的计算值，即为步骤142计算得到的差值；(fl，cl，alphal，kl)(fr，cr，alphar，kr)分别表示左、右摄像机的内参数；oml2r，Tl2r表示左摄像机到右摄像机的转换参数；omli，Tli表示左摄像机第i幅图像对应的外参数； 优化的方法采用Gauss-Newton迭代，迭代的步长应满足 JTWJδ＝-αJTWε 其中阻尼系数α取为1；加权矩阵W取单位阵I；J是ε对param的一阶导数； 其中， 所述的步骤16包括如下步骤 步骤161、记录摄像机1到其他摄像机的转换参数； 步骤162、对摄像机1的w-1组外参数进行中值滤波； 步骤163、对每一个摄像机的w-1组内参数进行中值滤波。最终得到的参数为各摄像机的内参数、摄像机1的外参数、摄像机1到其他摄像机的转换参数。
所述步骤2包括如下步骤 步骤21由步骤11-16所得到的各参数初值，计算特征点在各高分辨率摄像机成像平面上的像素坐标，并将此计算值与特征点在各高分辨率摄像机成像平面上像素坐标的测量值相减求差； 步骤22对各个参数进行非线性最小二乘优化，得到各高分辨率摄像机的内参数、摄像机1的外参数和摄像机1到其他摄像机的转换参数的最终值。
至此，高分辨率摄像机的标定工作完成，我们得到了高精度、具有较好鲁棒性的摄像机参数，包括各高分辨率摄像机的内参数，摄像机1的外参数，摄像机1到其他摄像机的转换参数。这些参数在后面的步骤中不再做调整。
所述的步骤21，包括 步骤211由摄像机1的外参数与摄像机1到其他摄像机的转换参数，计算其他各摄像机的外参数； 设两个摄像机C1和C2的外参数分别用R1、t1和R2、t2表示，摄像机C1到C2的旋转矢量和平移矢量用R、t表示，则 R2＝RR1 t2＝t+Rt1 (8) 步骤212由各摄像机的内参数初值和步骤211得到的各摄像机外参数，计算特征点在各摄像机成像平面上的像素坐标，并将此计算值与特征点在左右摄像机成像平面上像素坐标的测量值相减求差。计算的方法与步骤142所述方法相同。
所述步骤22包括 对各摄像机内参数和摄像机1的外参数，以及摄像机1到其他各摄像机的转换参数进行非线性最小二乘优化，使如下代价函数 最小，即 其中param为需要优化的摄像机参数 param＝(f1；c1；alpha1；k1；...；fw-1；cw-1；alphaw；kw；om1；T1；...； omw-1；Tw-1；om11，T11；...om1n；T1n) 其中
为第j个特征点Mj在在第k个摄像机的第i幅图像上的投影的测量值，m(param，Mj，k，i)为相应的计算值，即为步骤21计算得到的差值。fk，ck，alphak，kk表示第k个摄像机的内参数；omj，Tj表示摄像机1到第j+1个像机的转换参数；omli，Tli表示摄像机1的第i幅图像对应的外参数； 优化的方法采用Gauss-Newton迭代，迭代的步长应满足 JTWJδ＝-αJTWε J是e对param的一阶导数， 所述的步骤3，包括 步骤31、计算每个高分辨率摄像机到其他高分辨率摄像机的转换参数； 步骤32、将预先通过单摄像机标定获得的数据作为低分辨率摄像机内、外参数的初值； 步骤33、由低分辨率摄像机的内、外参数初值，计算特征点在低分辨率摄像机成像平面上的像素坐标，并将此计算值与特征点在低分辨率摄像机成像平面上像素坐标的测量值相减，求得差值，记为e_l，计算的方法与步骤142所述方法相同； 步骤34、进行低分辨率摄像机与高分辨率摄像机的联合标定。
所述步骤34包括如下步骤 步骤341、任取一个高分辨率摄像机，记为Cam_h； 步骤342、利用低分辨率摄像机的参数和Cam_h的参数重新计算各高分辨率摄像机的外参数； 步骤343、用各高分辨率摄像机的内参数与步骤342得到的各高分辨率摄像机的外参数，计算特征点在各高分辨率摄像机成像平面上的像素坐标，并将此计算值与特征点在各高分辨率分辨率摄像机成像平面上像素坐标的测量值相减，求得差值，记为e_h，计算的方法与步骤142所述方法相同； 步骤344、综合考虑e_l和e_h，对低分辨率摄像机的内、外参数，以及低分辨率摄像机到各高分辨率摄像机的转换参数进行非线性最小二乘优化； 步骤345、将剩下的w-1个高分辨率摄像机依次取为Cam_h，重复执行步骤342至步骤344； 步骤346、对得到的w组参数进行鲁棒性处理，得到最终参数，包括低分辨率摄像机的内、外参数，以及低分辨率摄像机到各高分辨率摄像机的转换参数。
所述的步骤342包括如下步骤 步骤3421、用低分辨率摄像机的外参数初值和Cam_h的外参数，计算低分辨率摄像机到Cam_h的转换参数，计算方法与步骤1411所述方法相同； 步骤3422、用低分辨率摄像机到Cam_h的转换参数和Cam_h到其他高分辨率摄像机的转换参数，计算低分辨率摄像机到其他高分辨率摄像机的转换参数； 假设摄像机C1到摄像机C2的转换参数为(R12，t12)，摄像机C2到摄像机C3的转换参数为(R23，t23)，则摄像机C1到摄像机C3的转换参数(R13，t13)为 R13＝R23R12 t13＝R23t12+t23 (10) 步骤3423、由低分辨率摄像机外参数的初值和步骤3422得到的低分辨率摄像机到其他高分辨率摄像机的转换参数，重新计算除Cam_h以外的各高分辨率摄像机的外参数。计算方法与步骤211所述方法相同。
所述步骤344包括 对低分辨率摄像机的内、外参数，以及低分辨率摄像机到各高分辨率摄像机的转换参数进行非线性最小二乘优化，使如下代价函数 最小，即 其中param为需要优化的摄像机参数 param＝(fl；cl；alphal；kl；oml；Tl；...；omw；Tw；oml1；Tl1；...；omln；Tln) 其中
为第j个特征点Mj在在第k个摄像机的第i幅图像上的投影的测量值，m(param，Mj，k，i)为相应的计算值，即为e＝[e_l，e_h]，(fl，cl，alphal，kl)表示低分辨率摄像机的内参数；omj，Tj表示低分辨率摄像机到第j个高分辨率摄像机的转换参数；omli，Tli表示低分辨率摄像机的第i幅图像对应的外参数； 优化的方法采用Gauss-Newton迭代，迭代的步长应满足JTWJδ＝-αJTWε J是e对param的一阶导数， 所述的步骤346包括对得到的w组低分辨率摄像机的内、外参数，以及低分辨率摄像机到各高分辨率摄像机的转换参数，分别进行中值滤波。
本发明的有益效果是1、使用本方法进行多摄像机阵列系统的标定，除平面模板外无需其他设备，节省了成本，且方便易行。
2、本方法适用于含有不同分辨率的摄像机阵列系统。
3、使用本方法可以显著提高标定的精度。使用各标定参数对模板特征点进行重投影，可以看出使用本方法前，各高分辨率摄像机和低分辨率摄像机均存在较大误差；而使用本方法后，高分辨率摄像机的误差达到了亚像素级别，低分辨率摄像机的精度也有了很大的提高。例如，对由D1分辨率的高分辨率摄像机和QVGA分辨率的低分辨率摄像机组成的阵列系统，运用本联合标定方法的结果如下所示 1、用高分辨率摄像机的参数进行重投影的误差 表1各摄像机x方向的平均误差 表1是用高分辨率摄像机的参数进行重投影得到的各摄像机x方向的平均误差，单位为像素。
表2各摄像机y方向的平均误差 表2是用高分辨率摄像机的参数进行重投影得到的各摄像机x方向的平均误差，单位为像素。
2、用低分辨率摄像机的参数进行重投影的误差 表3各摄像机x方向的平均误差 表3是用低分辨率摄像机的参数进行重投影得到的各摄像机x方向的平均误差，单位为像素。
表4各摄像机y方向的平均误差 表4是用低分辨率摄像机的参数进行重投影得到的各摄像机y方向的平均误差，单位为像素。
从表中可以看出，使用本联合标定方法后，用高分辨率摄像机和低分辨率摄像机参数进行重投影的误差都大为减小，标定精度有了很大提高。
图1为本发明一种多CCD摄像机与TOF摄像机组成的阵列系统实施例结构图； 图2为本发明一种多摄像机阵列联合标定方法实施例的总体框图； 图3为本发明一种高分辨率摄像机两两间标定方法实施例的流程图； 图4为本发明一种各高分辨率摄像机联合标定方法实施例的流程图； 图5为本发明一种低分辨率摄像机与高分辨率摄像机联合标定方法实施例的流程图。
具体实施例方式 下面，结合附图进一步说明本发明的具体实施方式
。
以由6个高分辨率的CCD摄像机和1个低分辨率的TOF(Time of Flight)摄像机组成的阵列系统为例，进行具体说明。该阵列系统的结构如图1所示。该阵列的所有摄像机位于同一平面上。假设平面模板上有100个特征点，各CCD摄像机和TOF摄像机同步地对模板拍摄了10张图像，进行单摄像机标定，获得了各自的内、外参数，并且已获得特征点在各摄像机成像平面上像素坐标的测量值。
如图2所示，本方法先使用单摄像机标定得到的数据对CCD摄像机和TOF摄像机的内、外参数进行初始化。然后先对CCD摄像机进行处理，包括CCD摄像机两两间的标定和各CCD摄像机的联合标定。最后对TOF摄像机和各CCD摄像机进行联合标定。
下面，分阶段对各步骤进行具体介绍。
第一阶段，CCD摄像机两两间的标定。其总体流程如图3所示。
设CCD摄像机已编号为Cam1，Cam2，…，Cam6。
步骤1任取6个CCD摄像机中的一对，将编号小的摄像机作为左摄像机，将编号大的摄像机作为右摄像机。由单标定得到的数据对左、右摄像机参数进行初始化，包括fl，cl，alphal，kl，omli，Tli和fr，cr，alphar，kr，omri，Tri，i＝1，...，10。
步骤2根据公式(1)，计算左摄像机到右摄像机的转换参数(oml2r，Tl2r)。
步骤3对10张图像，重复步骤1至步骤2。写成伪代码，即是 for i＝1n，％此处n＝10 Ri＝rodrigues(omri)*inv(rodrigues(omli)) omi＝rodrigues(Ri) Ti＝Tri-RiTli end； 步骤4对步骤3得到的10组(omi，Ti)参数进行中值滤波，得到左摄像机到右摄像机的转换参数初值。写成伪代码，即是 oml2r＝median((om1，om2，...，omn)T)T (n＝10) Tl2r＝median((T1，T2，...，Tn)T)T 步骤5由左、右摄像机的内、外参数的初始值，计算模板特征点在左、右摄像机成像平面上的像素坐标。
设XX＝[X；Y；Z]代表点P在世界坐标系中的坐标，XXc＝[Xc；Yc；Zc]代表点P在摄像机坐标系下的坐标，按公式(2)(3)(4)(5)进行计算。
步骤6对参数 param＝(fl；cl；alphal；kl；fr；cr；alphar；kr；oml2r；Tl2r；oml1；Tl1；...；om1n；T1n)，按公式(6)进行迭代优化。其中fl，cl，alphal，kl为左摄像机内参数(分别为摄像机焦距、主点坐标、倾斜因子、畸变系数)；fr，cr，alphar，kr为右摄像机内参数；oml2r，Tl2r为左摄像机到右摄像机的旋转矢量和平移矢量；omli，Tli，i＝1，...，10为左摄像机的外参数。
步骤7对每一对摄像机重复步骤1至步骤5。此处共有6*(6-1)/2＝15对。
步骤8对每一个摄像机，通过步骤6得到了5组内参数，对其进行中值滤波。
步骤9纪录Cam1到其他摄像机的转换参数，并对Cam1的5组外参数进行中值滤波。
经过这一阶段的优化，便得到了param的较好初始值。
第二阶段，各CCD摄像机间的联合标定。其总体流程如图4所示。
步骤1用第一阶段获得的较好的初值对各参数进行初始化。
步骤2用Cam1的外参数与Cam1到其他5个CCD摄像机的转换矢量(omltoi，Tltoi)，i＝2，3，...，6，依据公式(9)计算Cam2至Cam6的外参数。
步骤3用步骤2获得的外参数与第一阶段获得的内参数，根据公式(2)(3)(4)(5)，计算世界坐标系下各特征点在各高分辨率成像平面上的像素坐标。
步骤4对参数按公式(10)(7)进行迭代优化，其中fk，ck，alphak，kk表示摄像机Camk的内参数；omj，Tj表示Cam1到Camj的转换参数；omli，Tli表示Cam1的第i幅图像对应的外参数。
即为 ε＝((1)ε，...，(10)ε)T 其中，(i)e ＝((i)e1l，...，(i)e1m；(i)e2l，...，(i)e2m；...；(i)ekj；...；(i)e6l，...，(i)e6m)T，i＝1，...，10；m＝100。
其中，(i)ekj＝((i)x’kj-(i)xkj，(i)y’kj-(i)ykj)T，i＝1，...，10；k＝1，...，6；j＝1，...，100。
其中，((i)x’kj，(i)y’kj)T表示第k个摄像机第i帧中的第j个特征点的计算值在图像平面的坐标；((i)xkj，(i)ykj)T表示第k个摄像机第i帧中的第j个特征点的测量值在图像平面的坐标。
第三阶段，CCD摄像机与TOF摄像机的联合标定。其总体流程如图5所示。
步骤1计算出6个CCD摄像机间的转换参数。选取Cam1为Cam_CCD，并用单标定数据对TOF摄像机的内、外参数进行初始化，包括(ftof，ctof，alphatof，ktof，omtofi，Ttofi)，i＝1，...，10。
步骤2根据公式(1)计算TOF到Cam_CCD的转换参数，进而根据公式(11)计算出TOF摄像机到其他CCD摄像机的转换参数。
步骤3根据公式(9)，由TOF摄像机的外参数与步骤2得到的转换参数，重新计算各CCD摄像机(除了Cam_CCD)的外参数。
步骤4由步骤3得到的各CCD摄像机外参数与第二阶段得到的各CCD摄像机内参数，根据公式(2)(3)(4)(5)重新计算特征点在各CCD摄像机成像平面上的坐标，并将计算值与测量值相减得到e_ccd。
步骤5由TOF摄像机的内外参数，根据公式(2)(3)(4)(5)重新计算特征点在各TOF摄像机成像平面上的坐标，并将计算值与测量值相减得到e_tof。
步骤6对参数按公式(12)(7)进行迭代优化，其中ftof，ctof，alphatof，ktof)表示TOF摄像机的内参数；omj，Tj表示TOF摄像机到第j个CCD摄像机的转换参数；omtofi，Ttofi表示TOF摄像机的第i幅图像对应的外参数。
为TOF与CCD两方面误差之和，即为e＝[e_tof，e_ccd]，即ε＝((1)ε，...，(n)ε)T。
其中，(i)ε＝((i)etof1，...，(i)etofm，(i)eccd11，...，(i)eccd1m，...，(i)eccd6l...，(i)eccd6m)T， i＝1，...，10；m＝100。
其中
和((i)xccdkj，(i)yccdkj)表示第k个CCD摄像机的第i帧的第j个特征点的图像坐标计算值和测量值。
和
表示TOF摄像机的第i帧的第j个特征点的图像坐标计算值和测量值。
步骤7将剩余5个CCD摄像机依次取为Cam_CCD，重复步骤2至步骤6。
步骤8将步骤7得到的6组param进行中值滤波，得到最终的鲁棒性的参数。
权利要求
1、一种多摄像机阵列的联合标定方法，它由不同分辨率的摄像机组成多摄像机阵列系统，其技术特征在于，包括如下步骤
步骤1、对高分辨率摄像机两两间进行标定，以得到各高分辨率摄像机的内参数初值和各高分辨率摄像机间转换参数的初值；
步骤2、对各高分辨率摄像机进行联合标定，以得到各高分辨率摄像机的内参数和各高分辨率摄像机间转换参数的最终值；
步骤3、对高分辨率摄像机与低分辨率摄像机进行联合标定，以得到低分辨率摄像机与各高分辨率摄像机间转换参数的最终值。
2、如权利要求1所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤1包括如下步骤
步骤11、对各高分辨率摄像机进行编号，记为摄像机1、摄像机2、…、摄像机w，w为摄像机数；
步骤12、任取一对高分辨率摄像机，将编号小的摄像机作为左摄像机，将编号大的摄像机作为右摄像机；
步骤13、将预先通过单摄像机标定获得的内、外参数作为左摄像机、右摄像机内、外参数的初值；
步骤14、对左、右摄像机的参数进行非线性最小二乘优化；
步骤15、对每一对摄像机重复执行步骤12至步骤14；
步骤16、对得到的参数进行鲁棒性处理；
经过上述步骤的优化，得到每个高分辨率摄像机的内参数的初值，摄像机1的外参数的初值，摄像机1到其他高分辨率摄像机的转换矢量的初值。
3、如权利要求2所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤14包括如下步骤
步骤141、由左、右摄像机外参数的初值，计算其转换关系，作为左摄像机到右摄像机的转换参数的初值；
步骤142、由左、右摄像机的内、外参数的初值，计算模板特征点在左、右摄像机成像平面上的像素坐标，并将此计算值与特征点在左右摄像机成像平面上像素坐标的测量值相减求差；
步骤143、对左、右摄像机内参数和左摄像机外参数，以及步骤141得到的左摄像机到右摄像机的转换参数，进行非线性最小二乘优化。
4、如权利要求3所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤141包括如下步骤
步骤1411、将左、右摄像机的外参数分别用Rl、tl和Rr、tr表示，其中Rl、Rr表示旋转矩阵，tl、tr表示平移矩阵，则左摄像机到右摄像机的转换参数即旋转矩阵和平移矩阵Rl2r、tl2r为
tl2r＝tr-Rl2rtl (1)
步骤1412、对拍摄的n幅图像，即n组外参数，重复步骤1411；
步骤1413、对步骤1412得到的n组Rl2r、tl2r参数进行中值滤波，将结果作为左摄像机到右摄像机的转换参数的初值。
5、如权利要求3所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤142包括如下步骤
步骤1421、由左、右摄像机的外参数，计算特征点在左、右摄像机坐标系下的坐标，方法如下
设XX＝[X；Y；Z]代表点P在世界坐标系中的坐标，XXC＝[XC；YC；ZC]代表点P在摄像机坐标系下的坐标，那么XX与XXC的关系可以用如下刚体运动转化来表示
XXc＝rodrigues(om)*XX+T (2)
其中，rodrigues(.)表示对括号内容输入Rodrigues公式；
om、T分别表示摄像机外参数中的旋转矢量和平移矢量；
步骤1422、计算特征点归一化针孔投影的坐标Xn
r2＝x2+y2 (3)
步骤1423、计算考虑镜头畸变影响后的归一化点坐标xd
其中，畸变包括径向畸变与偏心畸变的切向分量，用5×1的矢量kc表示，其中k(1)、k(2)、k(5)分别是径向畸变的二次、四次和六次项，k(3)和k(4)分别是两个切向畸变系数；
步骤1424、计算特征点P在图像平面上最终的像素坐标pixel＝[xp，yp]
其中KK表示摄像机的内参数矩阵
其中f、c、alpha分别表示焦距、主点坐标、倾斜因子；
步骤1425、将特征点在左右摄像机成像平面上像素坐标的计算值与测量值相减求差。
6、如权利要求3所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤143包括
对左、右摄像机内参数和左摄像机外参数，以及左摄像机到右摄像机的转换参数的初值进行非线性最小二乘优化，使如下代价函数
最小，即
其中，param为需要优化的摄像机参数
param＝(fl；cl；alphal；kl；fr；cr；alphar；kr；oml2r；Tl2r；oml1；Tl1；...；omln；Tln)
为第j个特征点Mj在在第k个摄像机的第i幅图像上的投影的测量值，m(param，Mj，k，i)为相应的计算值；(fl，cl，alphal，kl)(fr，cr.alphar.kr)分别表示左、右摄像机的内参数；oml2r，Tl2r表示左摄像机到右摄像机的转换参数；omli，Tli表示左摄像机第i幅图像对应的外参数；
优化的方法采用Gauss-Newton迭代，迭代的步长应满足
JTWJδ＝-αJTWε
其中阻尼系数α取为1；加权矩阵W取单位阵I；J是ε对param的一阶导数；
其中，
7、如权利要求2所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤16包括如下步骤
步骤161、记录摄像机1到其他摄像机的转换参数；
步骤162、对摄像机1的w-1组外参数进行中值滤波；
步骤163、对每一个摄像机的w-1组内参数进行中值滤波。
8、如权利要求2所述的多摄像机阵列的联合标定方法，其特征在于，所述的对各高分辨率摄像机进行联合标定步骤如下
步骤21由权利要求2所得到的各参数初值，计算特征点在各高分辨率摄像机成像平面上的像素坐标，并将此计算值与特征点在各高分辨率摄像机成像平面上像素坐标的测量值相减求差；
步骤22对各个参数进行非线性最小二乘优化，得到各高分辨率摄像机的内参数、摄像机1的外参数华和摄像机1到其他摄像机的转换参数的最终值。
9、如权利要求8所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤21，包括
步骤211由摄像机1的外参数与摄像机1到其他摄像机的转换参数，计算其他各摄像机的外参数；
设两个摄像机C1和C2的外参数分别用R1、t1和R2、t2表示，摄像机C1到C2的旋转矢量和平移矢量用R、t表示，则
R2＝RR1
t2＝t+Rt1 (8)
步骤212由各摄像机的内参数初值和步骤211得到的各摄像机外参数，计算特征点在各摄像机成像平面上的像素坐标，并将此计算值与特征点在左右摄像机成像平面上像素坐标的测量值相减求差。
10、如权利要求8所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤22包括
对各摄像机内参数和摄像机1的外参数，以及摄像机1到其他各摄像机的转换参数进行非线性最小二乘优化，使如下代价函数
最小，即
其中param为需要优化的摄像机参数
param＝(f1；c1；alpha1；k1；...；fw-1；cw-1；alphaw；kw；om1；T1；...；
omw-1；Tw-1；oml1，Tl1...omlx；Tlx)
其中
为第j个特征点Mj在在第k个摄像机的第i幅图像上的投影的测量值，m(param，Mj，k，i)为相应的计算值；fk，ck，alphak，kk表示第k个摄像机的内参数；omj，Tj表示摄像机1到第j+1个像机的转换参数；omli，Tli表示摄像机1的第i幅图像对应的外参数；
优化的方法采用Gauss-Newton迭代，迭代的步长应满足
JTWJδ＝-αJTWε
J是e对param的一阶导数，
11、如权利要求1所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤3，包括
步骤31、计算每个高分辨率摄像机到其他高分辨率摄像机的转换参数；
步骤32、将预先通过单摄像机标定获得的数据作为低分辨率摄像机内、外参数的初值；
步骤33、由低分辨率摄像机的内、外参数初值，计算特征点在低分辨率摄像机成像平面上的像素坐标，并将此计算值与特征点在低分辨率摄像机成像平面上像素坐标的测量值相减，求得差值，记为e1；
步骤34、进行低分辨率摄像机与高分辨率摄像机的联合标定。
12、如权利要求11所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤34包括如下步骤
步骤341、任取一个高分辨率摄像机，记为Cam_h；
步骤342、利用低分辨率摄像机的参数和Cam_h的参数重新计算各高分辨率摄像机的外参数；
步骤343、用各高分辨率摄像机的内参数与步骤342得到的各高分辨率摄像机的外参数，计算特征点在各高分辨率摄像机成像平面上的像素坐标，并将此计算值与特征点在各高分辨率分辨率摄像机成像平面上像素坐标的测量值相减，求得差值，记为e_h；
步骤344、综合考虑e_l和e_h，对低分辨率摄像机的内、外参数，以及低分辨率摄像机到各高分辨率摄像机的转换参数进行非线性最小二乘优化；
步骤345、将剩下的w-1个高分辨率摄像机依次取为Cam_h，重复执行步骤342至步骤344；
步骤346、对得到的w组参数进行鲁棒性处理，得到最终参数，包括低分辨率摄像机的内、外参数，以及低分辨率摄像机到各高分辨率摄像机的转换参数。
13、如权利要求12所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤342包括如下步骤
步骤3421、用低分辨率摄像机的外参数初值和Cam_h的外参数，计算低分辨率摄像机到Cam_h的转换参数；
步骤3422、用低分辨率摄像机到Cam_h的转换参数和Cam_h到其他高分辨率摄像机的转换参数，计算低分辨率摄像机到其他高分辨率摄像机的转换参数；
假设摄像机C1到摄像机C2的转换参数为(R12，t12)，摄像机C2到摄像机C3的转换参数为(R23，t23)，则摄像机C1到摄像机C3的转换参数(R13，t13)为
R13＝R23R12
t13＝R23t12+t23(10)
步骤3423、由低分辨率摄像机外参数的初值和步骤3422得到的低分辨率摄像机到其他高分辨率摄像机的转换参数，重新计算除Cam_h以外的各高分辨率摄像机的外参数。
14、如权利要求12所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤344包括
对低分辨率摄像机的内、外参数，以及低分辨率摄像机到各高分辨率摄像机的转换参数进行非线性最小二乘优化，使如下代价函数
最小，即
其中param为需要优化的摄像机参数
param＝(fl；cl；alphal；kl；oml；Tl；...；omw；Tw；omll；Tll；...；omln；Tln)
其中
为第j个特征点Mj在在J投影的测量值，m(param，Mj，k，i)为相应的计算值， 为e＝[e_l，e_h]，(fl，cl，alphal，kl)表示低分辨率摄像机的内参数；omj，Tj表示低分辨率摄像机到第j个高分辨率摄像机的转换参数；omli，Tli表示低分辨率摄像机的第i幅图像对应的外参数；
优化的方法采用Gauss-Newton迭代，迭代的步长应满足
JTWJδ＝-αJTWε
J是e对param的一阶导数，
15、如权利要求12所述的多摄像机阵列的联合标定方法，其特征在于，所述的步骤346包括对得到的w组低分辨率摄像机的内、外参数，以及低分辨率摄像机到各高分辨率摄像机的转换参数，分别进行中值滤波。
全文摘要
本发明公开了一种多摄像机阵列的联合标定方法，它由不同分辨率的摄像机组成多摄像机阵列系统，包括如下步骤1.对高分辨率摄像机两两间进行标定，以得到各高分辨率摄像机的内参数初值和各高分辨率摄像机间转换参数的初值；2.对各高分辨率摄像机进行联合标定，以得到各高分辨率摄像机的内参数和各高分辨率摄像机间转换参数的最终值；3.对高分辨率摄像机与低分辨率摄像机进行联合标定，以得到低分辨率摄像机与各高分辨率摄像机间转换参数的最终值。采用本发明标定方法的优点是标定精度高，保证了鲁棒性，且适用于具有不同分辨率的摄像机阵列系统。
文档编号H04N5/222GK101420520SQ20081016252
公开日2009年4月29日 申请日期2008年12月1日 优先权日2008年12月1日
发明者晴 雷, 刘云海, 于慧敏 申请人:浙江大学