CN104036512A - 一种基于正交消隐点的新型Tsai’s摄像机标定改进方法 - Google Patents

Abstract

本发明提出一种基于正交消隐点的新型Tsai’s摄像机标定改进方法，首先拍摄若干幅棋盘格标定板平面图像，并对所有图像进行预处理，对每一幅图像的角点进行亚像素坐标提取，然后将角点坐标划分成四组，应用最小二乘指标函数拟合每个方向上的平行直线束，计算每个方向上的最佳消隐点坐标，接着建立摄像机内部参数的约束方程组，求解内部参数f_x,f_y,u₀,v₀，将这些计算值作为初始值，代入外部参数的方程组中，求解所有的外部参数；最后引入一阶、二阶径向畸变系数，在实际图像平面坐标系上，建立图像残差当量最小化寻优函数，对所有的参数进行优化。本发明大大减少了计算量，速度快、精度高，可广泛应用于机器视觉研究、工业三维测量等多个领域的摄像机标定。

Description

一种基于正交消隐点的新型Tsai’s摄像机标定改进方法

技术领域

本发明属于光学测量技术领域，涉及一种基于正交消隐点的新型Tsai’s摄像机标定改进方法。

背景技术

在计算机视觉测量领域中，摄像机标定是从2D图像信息获取3D空间信息必不可少的步骤，具体是根据获取的2D图像信息来确定摄像机内部参数，畸变参数和外部参数的过程。

目前，学者们在传统标定方法的理论基础上，已经研究出了众多简单灵活的摄像机新型标定方法。这些新型标定方法：一是在标定模板上进行不断创新，从最初工艺要求极高的3D精密标定立方块到现今的2D平面模板和1D线性模板，其中2D的主要有圆点靶标、十字叉靶标、棋盘格靶标等，1D的要属单条结构光最为广泛应用；二是在标定的理论方法上改进与创新，自标定方法中的消隐点、圆环点、Kruppa方程、绝对二次曲线、绝对二次曲面等。其中，B.CAPRILE和P.Beardsley等人首次提出了利用消隐点的几何特性来标定摄像机。

近几年，研究者从消隐点的产生模型、方法、拍摄图像的数量等方面不断改进与创新，出现了一些便捷灵活的新标定方法。比如在消隐点的产生模型有两对正交的平行直线、非平行的矩形对、圆柱、长方块、正方形等，这些模型虽然选取简易灵活，但不适合于精度要求高的场合，且有一个共同的问题：仅仅只凭借一对平行直线产生消隐点，这样会出现消隐点退化的情形，即在图像上不能够形成消隐点。

而在方法上，将消隐点作为计算圆环点的过度量，应用于建立标定的约束方程；或者通过消隐点与摄像机光心轨迹的连线所形成的球，建立约束方程进行标定等，这些方法并没有直接利用正交消隐点的性质来建立与摄像机内部参数的约束方程，导致计算量较大，过程复杂。

另外，对Tsai标定法改进的现有技术中，普遍采用Tsai法的计算思路，一是人为的约定内部参数的主点坐标就是拍摄的图像中心，这在现今的CCD摄像机中不再适合，由于镜头制造工艺存在偏差，导致主点会偏离图像中心；二是在计算方程组的过程中，依然采用先计算t_y这个未知量，这样会引入判定未知量正负号的麻烦，增大了计算量，降低了运算效率。

发明内容

本发明目的是解决上述背景技术中提出的三个技术问题，为此公开了一种基于正交消隐点的新型Tsai’s摄像机标定改进方法。

本发明的技术方案为：

所述一种基于正交消隐点的新型Tsai’s摄像机标定改进方法，其特征在于：采用以下步骤：

Step1：搭建光学测量系统，建立CCD摄像机模型；

Step2：使摄像机光轴与棋盘格标定板平面所成夹角在0°～80°之间，拍摄若干幅不同位姿的棋盘格标定板平面图像：并对拍摄的图像进行中值滤波、锐化预处理；棋盘格标定板的图案规格为m×n，其中m＝n或m＝n-1；

Step3：应用Harris角点提取算法对预处理后的图像提取角点亚像素坐标，角点亚像素坐标为s_i,j＝[u_i,j,v_i,j]^T,i＝1,2,...m,j＝1,2,...n；

Step4：将Step3中提取的角点划分成水平、垂直、主对角线、副对角线四组S＝{^gs|g＝hor,ver,mad,vid}；hor对应水平，ver对应垂直，mad对应主对角线，vid对应副对角线：

4.1若m＝n，则执行步骤4.2→4.6；若m＝n-1，则执行步骤4.7→4.11；

4.2定义棋盘格标定板的左下角第一个角点为s_1,1，水平向右依次为s_1,1,s_1,2,...,s_1,n，垂直向上依次为s_1,1,s_2,1,...,s_n,1；

4.3第一组：水平方向的角点坐标集合为

^hors＝{^hors¹,^hors²,...,^horsⁱ,...,^horsⁿ},i＝1,2,...,n，

其中子集合^horsⁱ＝{s_i,1,s_i,2,...s_i,j,...s_i,n},j＝1,2,...n；

4.4第二组：垂直方向的角点坐标集合为

^vers＝{^vers¹,^vers²,...,^vers^j,...,^versⁿ},j＝1,2,...,n，

其中子集合^vers^j＝{s_1,j,s_2,j,...s_i,j,...s_n,j},i＝1,2,...n；

4.5第三组：主对角线方向的角点坐标集合为

^mads＝{^mads¹,^mads²,...,^mads^k,...,^mads^2n-3},k＝1,2,...,2n-3；

其中：

1)当k＝1时，子集合^mads¹＝{s_2,1,s_1,2}，表示第2行第1列的角点和第1行第2列的角点组成的集合；

2)当k＝2时，子集合^mads²＝{s_3,1,s_2,2,s_1,3}，表示第3行第1列、第2行第2列和第1行第3列的角点组成的集合；

3)当k＝n-1时，子集合^mads^n-1＝{s_n,1,s_n-1,2,...,s_n-i,j,...,s_1,n},i＝0,1,...,n-1；j＝1,2,...,n；

4)当k＝n时，子集合^madsⁿ＝{s_n,2,s_n-1,3,...,s_n-i,j,...,s_2,n},i＝0,1,...,n-2；j＝2,3,...,n；

5)当k＝2n-4时，子集合^mads^2n-4＝{s_n,n-2,s_n-1,n-1,s_n-2,n}；

6)当k＝2n-3时，子集合^mads^2n-3＝{s_n,n-1,s_n-1,n}；

4.6第四组：副对角方向的角点坐标集合为

^vids＝{^vids¹,^vids²,...,^vids^q,...,^vids^2n-3},q＝1,2,...,2n-3；

其中：

1)当q＝1时，子集合^vids¹＝{s_n-1,1,s_n,2}，表示第n-1行第1列的角点和第n行第2列的角点组成的集合；

2)当q＝2时，子集合^vids²＝{s_n-2,1,s_n-1,2,s_n,3}，表示第n-2行第1列、第n-1行第2列和第n行第3列的角点组成的集合；

3)当q＝n-1时，子集合^vids^n-1＝{s_1,1,s_2,2,...,s_i,j},i,j＝1,2,...,n；

4)当q＝n时，子集合^vidsⁿ＝{s_1,2,s_2,3,...,s_i,j},i＝1,2,...,n-1；j＝2,3,...,n；

5)当q＝2n-4时，子集合^vids^2n-4＝{s_1,n-2,s_2,n-1,s_3,n}；

6)当q＝2n-3时，子集合^vids^2n-3＝{s_1,n-1,s_2,n}；

4.7定义棋盘格标定板的左下角第一个角点为s_1,1，水平向右依次为s_1,1,s_1,2,...,s_1,n，垂直向上依次为s_1,1,s_2,1,...,s_m,1；

4.8第一组：水平方向的角点坐标集合为

^hors＝{^hors¹,^hors²,...,^horsⁱ,...,^hors^m},i＝1,2,...,m，

其中子集合^horsⁱ＝{s_i,1,s_i,2,...s_i,j,...s_i,n},j＝1,2,...n；

4.9第二组：垂直方向的角点坐标集合为

^vers＝{^vers¹,^vers²,...,^vers^j,...,^versⁿ},j＝1,2,...,n，

其中子集合^vers^j＝{s_1,j,s_2,j,...s_i,j,...s_m,j},i＝1,2,...m；

4.10第三组：主对角线方向的角点坐标集合为

^mads＝{^mads¹,^mads²,...,^mads^k,...,^mads^m+n-3},k＝1,2,...,m+n-3；

其中：

1)当k＝1时，子集合^mads¹＝{s_2,1,s_1,2}，表示第2行第1列的角点和第1行第2列的角点组成的集合；

2)当k＝2时，子集合^mads²＝{s_3,1,s_2,2,s_1,3}，表示第3行第1列、第2行第2列和第1行第3列的角点组成的集合；

3)当k＝m-1时，子集合^mads^m-1＝{s_m,1,s_m-1,2,...,s_m-i,j},i＝0,1,...,m-1；j＝1,2,...,n-1；

4)当k＝m时，子集合^mads^m＝{s_m,2,s_m-1,3,...,s_m-i,j},i＝0,1,...,m-1；j＝2,3,...,n；

5)当k＝m+n-4时，子集合^mads^m+n-4＝{s_m,n-2,s_m-1,n-1,s_m-2,n}；

6)当k＝m+n-3时，子集合^mads^m+n-3＝{s_m,n-1,s_m-1,n}；

4.11第四组：副对角方向的角点坐标集合为

^vids＝{^vids¹,^vids²,...,^vids^q,...,^vids^m+n-3},q＝1,2,...,m+n-3；

其中：

1)当q＝1时，子集合^vids¹＝{s_m-1,1,s_m,2}，表示第m-1行第1列的角点和第m行第2列的角点组成的集合；

2)当q＝2时，子集合^vids²＝{s_m-2,1,s_m-1,2,s_m,3}，表示第m-2行第1列、第m-1行第2列和第m行第3列的角点组成的集合；

3)当q＝n-1时，子集合^vids^n-1＝{s_1,1,s_2,2,...,s_i,j},i,j＝1,2,...,n；

4)当q＝n时，子集合^vidsⁿ＝{s_1,2,s_2,3,...,s_i,j},i＝1,2,...,m；j＝2,3,...,n；

5)当q＝m+n-4时，子集合^vids^m+n-4＝{s_1,n-2,s_2,n-1,s_3,n}；

6)当q＝m+n-3时，子集合^vids^m+n-3＝{s_1,n-1,s_2,n}；

Step5：应用Step4中角点的划分结果，采用最小二乘法计算hor、ver、mad、vid方向上的平行直线束参数：

若m＝n，则平行直线束参数为：

^horlⁱ[^horaⁱ,-1,^horbⁱ],i＝1,2,...,n

^verl^j[^vera^j,-1,^verb^j],j＝1,2,...,n

^madl^k[^mada^k,-1,^madb^k],k＝1,2,...,2n-3

^vidl^q[^vida^q,-1,^vidb^q],q＝1,2,...,2n-3

若m＝n-1，则平行直线束参数为：

^horlⁱ(^horaⁱ,-1,^horbⁱ),i＝1,2,...,m

^verl^j(^vera^j,-1,^verb^j),j＝1,2,...,n

^madl^k(^mada^k,-1,^madb^k),k＝1,2,...,m+n-3

^vidl^q(^vida^q,-1,^vidb^q),q＝1,2,...,m+n-3

Step6：建立hor、ver、mad、vid方向的最佳消隐点目标函数：

$F_{\mathrm{opt}}^{\mathrm{hor}}=\underset{l_i^{\mathrm{hor}}}{\arg \min }\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{hor}},v_{\mathrm{va}}^{\mathrm{hor}}),l[a_i^{\mathrm{hor}},-1,b_i^{\mathrm{hor}}\left]\right]=\underset{l_i^{\mathrm{hor}}}{\arg \min }\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\frac{|a_i^{\mathrm{hor}}u_{\mathrm{va}}^{\mathrm{hor}}-v_{\mathrm{va}}^{\mathrm{hor}}+b_i^{\mathrm{hor}}|}{\sqrt{{\left(a_i^{\mathrm{hor}}\right)}^2+{(-1)}^2}};$

$F_{\mathrm{opt}}^{\mathrm{ver}}=\underset{l_j^{\mathrm{ver}}}{\arg \min }\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{ver}},v_{\mathrm{va}}^{\mathrm{ver}}),l[a_j^{\mathrm{ver}},-1,b_j^{\mathrm{ver}}\left]\right]=\underset{l_j^{\mathrm{ver}}}{\arg \min }\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{|a_j^{\mathrm{ver}}u_{\mathrm{va}}^{\mathrm{ver}}-v_{\mathrm{va}}^{\mathrm{ver}}+b_j^{\mathrm{ver}}|}{\sqrt{{\left(a_j^{\mathrm{ver}}\right)}^2+{(-1)}^2}};$

$F_{\mathrm{opt}}^{\mathrm{mad}}=\underset{l_k^{\mathrm{mad}}}{\arg \min }\underset{k=1}{\overset{m+n-3}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{mad}},v_{\mathrm{va}}^{\mathrm{mad}}),l[a_k^{\mathrm{mad}},-1,b_k^{\mathrm{mad}}\left]\right]=\underset{l_k^{\mathrm{hor}}}{\arg \min }\underset{k=1}{\overset{m+n-3}{\mathrm{\Σ}}}\frac{|a_k^{\mathrm{mad}}u_{\mathrm{va}}^{\mathrm{mad}}-v_{\mathrm{va}}^{\mathrm{mad}}+b_k^{\mathrm{mad}}|}{\sqrt{{\left(a_k^{\mathrm{mad}}\right)}^2+{(-1)}^2}};$

$F_{\mathrm{opt}}^{\mathrm{vid}}=\underset{l_q^{\mathrm{vid}}}{\arg \min }\underset{q=1}{\overset{m+n-3}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{vid}},v_{\mathrm{va}}^{\mathrm{vid}}),l[a_q^{\mathrm{vid}},-1,b_q^{\mathrm{vid}}\left]\right]=\underset{l_q^{\mathrm{vid}}}{\arg \min }\underset{q=1}{\overset{m+n-3}{\mathrm{\Σ}}}\frac{|a_q^{\mathrm{vid}}u_{\mathrm{va}}^{\mathrm{vid}}-v_{\mathrm{va}}^{\mathrm{vid}}+b_q^{\mathrm{vid}}|}{\sqrt{{\left(a_q^{\mathrm{vid}}\right)}^2+{(-1)}^2}};$

并应用Levenberg-Marquard算法对最佳消隐点目标函数寻优计算出最佳的消隐点坐标^gV＝[^gu_va,^gv_va]^T；

Step7：利用消隐线，验证消隐点计算的是否合理：

7.1对Step6中4个消隐点进行随机组合，生成种组合，并拟合出条直线，计算这些直线两两间的夹角，选择出最大的夹角θ_max；

7.2判断若θ≤θ_max，则这4个消隐点共线，形成了消隐线；反之，转到Step3，重新提取角点亚像素坐标，其中θ为阈值；

Step8：建立摄像机内部参数几何约束方程组：

$\left\{\begin{array}{c}\hat{V}_{\mathrm{hor}}\mathrm{\ω}\hat{V}_{\mathrm{ver}}=0\\ \hat{V}_{\mathrm{mad}}\mathrm{\ω}\hat{V}_{\mathrm{vid}}=0\end{array}\right.,\hat{V}_g={[u_{\mathrm{va}}^g,v_{\mathrm{va}}^g,1]}^T$

其中：

$\mathrm{\ω}=K^{-T}{K}^{-1}=\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& -\frac{u_0}{f_x^2}\\ 0& \frac{1}{f_y^2}& -\frac{v_0}{f_y^2}\\ -\frac{u_0}{f_x^2}& -\frac{v_0}{f_y^2}& \frac{u_0^2}{f_x^2}+\frac{v_0^2}{f_y^2}+1\end{array}\right]=\left[\begin{array}{ccc}{W}_{11}& 0& W_{13}\\ 0& W_{22}& W_{23}\\ W_{13}& W_{23}& W_{33}\end{array}\right]$

f_x,f_y为CCD摄像机模型u轴和v轴上的归一化像素焦距，(u₀,v₀)是CCD摄像机模型成像平面的中心坐标；根据摄像机内部参数几何约束方程组，采用奇异值(SVD)分解和Zhang’s的方法计算该内部参数f_x,f_y,u₀,v₀；

Step9：计算CCD摄像机模型中摄像机外部参数平移向量T的Z_C轴方向的分量t_z为：

9.1建立摄像机外部参数的方程组：

$\left\{\begin{array}{c}{u}_{i,j}=f_x\frac{r_{1}X_{i,j}^W+r_{2}Y_{i,j}^W+t_x}{r_{7}X_{i,j}^W+r_{8}Y_{i,j}^W+t_z}+u_0\\ v_{i,j}=f_y\frac{r_{4}X_{i,j}^W+r_{5}Y_{i,j}^W+t_y}{r_{7}X_{i,j}^W+r_{8}Y_{i,j}^W+t_z}+v_0\end{array}\right.,i=\mathrm{1,2},...,m;j=\mathrm{1,2},...,n$

其中：摄像机外部参数的平移向量

$T=\left[\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right]$

摄像机外部参数的旋转矩阵

$R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right]$

^WX_i,j和^WY_i,j为棋盘格标定板中角点的空间坐标；

9.2将摄像机外部参数的方程组变形为矩阵形式为：

$\left[\begin{array}{cccccccc}{f}_{x}X_{i,j}^W& f_{x}Y_{i,j}^W& f_x& 0& 0& 0& {-U}_{x}X_{i,j}^W& {-U}_{x}Y_{i,j}^W\\ 0& 0& 0& f_{y}X_{i,j}^W& f_{y}Y_{i,j}^W& f_y& -V_{x}X_{i,j}^W& {-V}_{x}Y_{i,j}^W\end{array}\right]b=\left[\begin{array}{c}{U}_x\\ V_x\end{array}\right]$

式中：向量b＝[b₁ b₂ b₃ b₄ b₅ b₆ b₇ b₈]^T；U_x＝u_i,j-u₀；V_x＝v_i,j-v₀；

9.3应用最小二乘法计算9.2中矩阵形式，得到向量b＝[b₁ b₂ b₃ b₄ b₅ b₆ b₇ b₈]^T；并根据公式

t_z＝(t_z1+t_z2+t_z3)/3

$t_{z1}=\sqrt{({\mathrm{Sr}}_1-\sqrt{{{\mathrm{Sr}}_1}^2-{4\mathrm{Se}}_1})/{2\mathrm{Se}}_1},{\mathrm{Sr}}_1=b_1^2+b_2^2+b_4^2+b_5^2,{\mathrm{Se}}_1={(b_1{b}_5-b_2{b}_4)}^2$

$t_{z2}=\sqrt{({\mathrm{Sr}}_2-\sqrt{{{\mathrm{Sr}}_2}^2-{4\mathrm{Se}}_2})/{2\mathrm{Se}}_2},{\mathrm{Sr}}_2=b_4^2+b_5^2+b_7^2+b_8^2,{\mathrm{Se}}_2={(b_5{b}_7-b_4{b}_8)}^2$

$t_{z3}=\sqrt{({\mathrm{Sr}}_3-\sqrt{{{\mathrm{Sr}}_3}^2-{4\mathrm{Se}}_3})/{2\mathrm{Se}}_3},{\mathrm{Sr}}_3=b_1^2+b_2^2+b_7^2+b_8^2,{\mathrm{Se}}_3={(b_1{b}_8-b_2{b}_7)}^2$

得到t_z；

Step10：计算旋转矩阵R中的分量r₁,r₂,r₃,r₄,r₅,r₆,r₇,r₈,r₉和平移向量中的分量t_x,t_y；

Step11：引入一、二阶径向畸变系数k₁,k₂，对内部参数f_x,f_y,u₀,v₀和外部参数R,T进行非线性全局优化：

在摄像机的实际图像平面坐标系UV内建立图像残差当量最小化寻优函数：

式中：：3D空间点^WS从世界坐标系X_WY_WZ_W投影到物理图像平面坐标系UV的畸变2D点实际图像平面坐标；：3D空间点^WS从像素图像平面坐标系uv逆向投影到物理图像平面坐标系UV的畸变2D点实际图像坐标。

有益效果

本发明的实现过程为：首先搭建光学测量系统，建立摄像机模型，使摄像机光轴与标定板平面所成夹角在0°～80°之间，拍摄若干幅不同位姿的棋盘格标定板平面图像，并对所有图像进行中值滤波、锐化预处理，且应用Harris算法对预处理的每一幅图像的角点进行亚像素坐标提取；然后按照棋盘格的水平、垂直、主对角线和副对角线四个方向的纹理特征，将这些角点坐标划分成四组，应用最小二乘指标函数拟合每个方向上的平行直线束，进而建立最佳消隐点目标函数，计算每个方向上的最佳消隐点坐标，并通过消隐线来判定消隐点的正确性；接着应用正交消隐点的性质，直接建立摄像机内部参数的约束方程组，求解内部参数f_x,f_y,u₀,v₀，将这些计算值作为初始值，代入外部参数的方程组中，通过先计算t_z来避免学者们在计算过程中引入的判定未知量正负号的步骤，进而求解所有的外部参数；最后引入一阶、二阶径向畸变系数，在实际图像平面坐标系上，建立图像残差当量最小化寻优函数，对所有的参数进行优化。相比于现有技术，本发明大大减少了计算量，速度快、精度高，可广泛应用于机器视觉研究、工业三维测量等多个领域的摄像机标定。

附图说明

图1一种基于新型正交消隐点的Tsai’s摄像机标定改进方法的流程图；

图2摄像机模型示意图；

图3常见的消隐点产生模型；

其中，(a)正交的平行线；(b)非平行的矩形；(c)正方形；(d)圆柱；(e)长方体(正方体)；(f)退化；

图4新型的消隐点产生模型；

图5实际拍摄的3幅棋盘格标定板图像；

图6划分后的3幅角点坐标图像；

图7实验计算出的消隐点坐标及形成的消隐线图像；

图8非线性全局优化模型。

具体实施方式

下面结合具体实施例描述本发明：

Step1：搭建光学测量系统，建立CCD摄像机模型；

上述Step1的具体实现方法是：

1.1不考虑摄像机镜头畸变，建立线性CCD摄像机针孔模型：

式中：

K＝M₂M₁＝[f_x 0 u₀；0 f_y v₀；0 0 1]：摄像机的内部参数矩阵；

(f_x,f_y)＝(f/d_x,f/d_y)：u轴和v轴上的归一化像素焦距；

f：有效的物理焦距；

(d_x,d_y)：CCD传感器单元的尺寸大小

(u₀,v₀)：成像平面的中心坐标；

$R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right]$ ：摄像机外部参数的旋转矩阵；

$T=\left[\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right]$ ：摄像机外部参数的平移向量；

^WS＝[^WX,^WY,^WZ]^T：3D空间中的点； ^WS的齐次坐标；

s＝[u,v]^T：像素图像平面坐标系uv的2D点；s的齐次坐标；

λ：不确定因子；

1.2引入一、二阶径向畸变，建立畸变的数学模型：

式中：

k₁,k₂：一、二阶径向畸变系数；

$\mathrm{\Δr}=\sqrt{{(U-U^{*})}^2+{(V-V^{*})}^2};$

S_I＝[U,V]^T：物理图像平面坐标系UV的无畸变2D点坐标；

S^*＝[U^*,V^*]^T：畸变后的图像平面坐标；

Step2：使摄像机光轴与棋盘格标定板平面所成夹角在0°～80°之间，拍摄若干幅不同位姿的棋盘格标定板平面图像：并对拍摄的图像进行中值滤波、锐化预处理；棋盘格标定板的图案规格为m×n，其中m＝n或m＝n-1；

本实施例中：

(1)设置摄像机的拍摄参数，见表1：

表1 摄像机拍摄参数

(2)使用的棋盘格标定板图案规格为12×13，格子间距为6mm，为了避免与摄像机光轴垂直，使摄像机光轴与标定板平面所成夹角在0°～80°之间，拍摄3幅不同位姿的棋盘格标定板平面图像，见图5。

(3)对这3幅图像进行中值滤波、锐化预处理。

Step3：应用Harris角点提取算法对预处理后的图像提取角点亚像素坐标，角点亚像素坐标为s_i,j＝[u_i,j,v_i,j]^T,i＝1,2,…m,j＝1,2,…n；

本实施例中，在Matlab平台上，采用Harris原理提取格子的角点，Image1的部分角点坐标见表2；

表2 Image1的部分角点坐标

Step4：将Step3中提取的角点划分成水平、垂直、主对角线、副对角线四组S＝{^gs|g＝hor,ver,mad,vid}；hor对应水平，ver对应垂直，mad对应主对角线，vid对应副对角线：

4.1若m＝n，则执行步骤4.2→4.6；若m＝n-1，则执行步骤4.7→4.11；

4.2定义棋盘格标定板的左下角第一个角点为s_1,1，水平向右依次为s_1,1,s_1,2,…,s_1,n，垂直向上依次为s_1,1,s_2,1,…,s_n,1；

4.3第一组：水平方向的角点坐标集合为

^hors＝{^hors¹,^hors²,…,^horsⁱ,…,^horsⁿ},i＝1,2,…,n，

其中子集合^horsⁱ＝{s_i,1,s_i,2,…s_i,j,…s_i,n},j＝1,2,…n；

4.4第二组：垂直方向的角点坐标集合为

^vers＝{^vers¹,^vers²,…,^vers^j,…,^versⁿ},j＝1,2,…,n，

其中子集合^vers^j＝{s_1,j,s_2,j,…s_i,j,…s_n,j},i＝1,2,…n；

4.5第三组：主对角线方向的角点坐标集合为

^mads＝{^mads¹,^mads²,…,^mads^k,…,^mads^2n-3},k＝1,2,…,2n-3；

其中：

1)当k＝1时，子集合^mads¹＝{s_2,1,s_1,2}，表示第2行第1列的角点和第1行第2列的角点组成的集合；

2)当k＝2时，子集合^mads²＝{s_3,1,s_2,2,s_1,3}，表示第3行第1列、第2行第2列和第1行第3列的角点组成的集合；

3)当k＝n-1时，子集合^mads^n-1＝{s_n,1,s_n-1,2,…,s_n-i,j,…,s_1,n},i＝0,1,…,n-1；j＝1,2,…,n；

4)当k＝n时，子集合^madsⁿ＝{s_n,2,s_n-1,3,…,s_n-i,j,…,s_2,n},i＝0,1,…,n-2；j＝2,3,…,n；

5)当k＝2n-4时，子集合^mads^2n-4＝{s_n,n-2,s_n-1,n-1,s_n-2,n}；

6)当k＝2n-3时，子集合^mads^2n-3＝{s_n,n-1,s_n-1,n}；

4.6第四组：副对角方向的角点坐标集合为

^vids＝{^vids¹,^vids²,…,^vids^q,…,^vids^2n-3},q＝1,2,…,2n-3；

其中：

1)当q＝1时，子集合^vids¹＝{s_n-1,1,s_n,2}，表示第n-1行第1列的角点和第n行第2列的角点组成的集合；

2)当q＝2时，子集合^vids²＝{s_n-2,1,s_n-1,2,s_n,3}，表示第n-2行第1列、第n-1行第2列和第n行第3列的角点组成的集合；

3)当q＝n-1时，子集合^vids^n-1＝{s_1,1,s_2,2,…,s_i,j},i,j＝1,2,…,n；

4)当q＝n时，子集合^vidsⁿ＝{s_1,2,s_2,3,…,s_i,j},i＝1,2,…,n-1；j＝2,3,…,n；

5)当q＝2n-4时，子集合^vids^2n-4＝{s_1,n-2,s_2,n-1,s_3,n}；

6)当q＝2n-3时，子集合^vids^2n-3＝{s_1,n-1,s_2,n}；

4.7定义棋盘格标定板的左下角第一个角点为s_1,1，水平向右依次为s_1,1,s_1,2,…,s_1,n，垂直向上依次为s_1,1,s_2,1,…,s_m,1；

4.8第一组：水平方向的角点坐标集合为

^hors＝{^hors¹,^hors²,…,^horsⁱ,…,^hors^m},i＝1,2,…,m，

其中子集合^horsⁱ＝{s_i,1,s_i,2,…s_i,j,…s_i,n},j＝1,2,…n；

4.9第二组：垂直方向的角点坐标集合为

^vers＝{^vers¹,^vers²,…,^vers^j,…,^versⁿ},j＝1,2,…,n，

其中子集合^vers^j＝{s_1,j,s_2,j,…s_i,j,…s_m,j},i＝1,2,…m；

4.10第三组：主对角线方向的角点坐标集合为

^mads＝{^mads¹,^mads²,…,^mads^k,…,^mads^m+n-3},k＝1,2,…,m+n-3；

其中：

1)当k＝1时，子集合^mads¹＝{s_2,1,s_1,2}，表示第2行第1列的角点和第1行第2列的角点组成的集合；

2)当k＝2时，子集合^mads²＝{s_3,1,s_2,2,s_1,3}，表示第3行第1列、第2行第2列和第1行第3列的角点组成的集合；

3)当k＝m-1时，子集合^mads^m-1＝{s_m,1,s_m-1,2,…,s_m-i,j},i＝0,1,…,m-1；j＝1,2,…,n-1；

4)当k＝m时，子集合^mads^m＝{s_m,2,s_m-1,3,…，s_m-i,j},i＝0,1,…,m-1；j＝2,3,…,n；

5)当k＝m+n-4时，子集合^mads^m+n-4＝{s_m,n-2,s_m-1,n-1,s_m-2,n}；

6)当k＝m+n-3时，子集合^mads^m+n-3＝{s_m,n-1,s_m-1,n}；

4.11第四组：副对角方向的角点坐标集合为

^vids＝{^vids¹,^vids²,…,^vids^q,…,^vids^m+n-3},q＝1,2,…,m+n-3；

其中：

1)当q＝1时，子集合^vids¹＝{s_m-1,1,s_m,2}，表示第m-1行第1列的角点和第m行第2列的角点组成的集合；

2)当q＝2时，子集合^vids²＝{s_m-2,1,s_m-1,2,s_m,3}，表示第m-2行第1列、第m-1行第2列和第m行第3列的角点组成的集合；

3)当q＝n-1时，子集合^vids^n-1＝{s_1,1,s_2,2,…,s_i,j},i,j＝1,2,…,n；

4)当q＝n时，子集合^vidsⁿ＝{s_1,2,s_2,3,…,s_i,j},i＝1,2,…,m；j＝2,3,…,n；

5)当q＝m+n-4时，子集合^vids^m+n-4＝{s_1,n-2,s_2,n-1,s_3,n}；

6)当q＝m+n-3时，子集合^vids^m+n-3＝{s_1,n-1,s_2,n}；

本实施例中，利用发明内容中的步骤4.7→4.11，将这些角点划分水平、垂直、主对角、副对角4组，附图6是部分角点坐标的划分，坐标见表3；

表3 Image1的部分角点坐标的划分

Step5：应用Step4中角点的划分结果，采用最小二乘法计算hor、ver、mad、vid方向上的平行直线束参数：

若m＝n，则平行直线束参数为：

^horlⁱ[^horaⁱ,-1,^horbⁱ],i＝1,2,…,n

^verl^j[^vera^j,-1,^verb^j],j＝1,2,…,n (3)

^madl^k[^mada^k,-1,^madb^k],k＝1,2,…,2n-3

^vidl^q[^vida^q,-1,^vidb^q],q＝1,2,…,2n-3

若m＝n-1，则平行直线束参数为：

^horlⁱ(^horaⁱ,-1,^horbⁱ),i＝1,2,…,m

^verl^j(^vera^j,-1,^verb^j),j＝1,2,…,n

^madl^k(^mada^k,-1,^madb^k),k＝1,2,…,m+n-3 (4)

^vidl^q(^vida^q,-1,^vidb^q),q＝1,2,…,m+n-3

Step5的具体实现方法是：

若m＝n时：

1)对水平方向上的n条平行直线组成的平行直线束建立第i行的直线最小二乘误差指标函数：

$E_{\mathrm{LS}}^{\mathrm{hor}}=\mathrm{\ϵ}(a_i^{\mathrm{hor}},b_i^{\mathrm{hor}})=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(\right[a_i^{\mathrm{hor}},-1,b_i^{\mathrm{hor}}\left]{[u_{i,j},v_{i,j},1]}^T\right)}^2=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(a_i^{\mathrm{hor}}{u}_{i,j}+b_i^{\mathrm{hor}}-v_{i,j})}^2---\left(5\right)$

2)对垂直方向上的n条平行直线组成的平行直线束建立第j列的直线最小二乘误差指标函数：

$E_{\mathrm{LS}}^{\mathrm{ver}}=\mathrm{\ϵ}(a_j^{\mathrm{ver}},b_j^{\mathrm{ver}})=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\right[a_j^{\mathrm{ver}},-1,b_j^{\mathrm{hor}}\left]{[u_{i,j},v_{i,j},1]}^T\right)}^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(a_j^{\mathrm{ver}}{u}_{i,j}+b_j^{\mathrm{ver}}-v_{i,j})}^2---\left(6\right)$

3)对主对角方向上的2n-3条平行直线组成的平行直线束建立第k条的直线最小二乘误差指标函数：

$E_{\mathrm{LS}}^{\mathrm{mad}}=\mathrm{\ϵ}(a_k^{\mathrm{mad}},b_k^{\mathrm{mad}})=\underset{k=1}{\overset{2n-3}{\mathrm{\Σ}}}{\left(\right[a_k^{\mathrm{mad}},-1,b_k^{\mathrm{mad}}\left]{[u_{i,j},v_{i,j},1]}^T\right)}^2=\underset{j=1}{\overset{2n-3}{\mathrm{\Σ}}}{(a_k^{\mathrm{mad}}{u}_{i,j}+b_k^{\mathrm{mad}}-v_{i,j})}^2---\left(7\right)$

4)对副对角方向上的2n-3条平行直线组成的平行直线束建立第q条的直线最小二乘误差指标函数：

$E_{\mathrm{LS}}^{\mathrm{vid}}=\mathrm{\ϵ}(a_q^{\mathrm{vid}},b_q^{\mathrm{vid}})=\underset{q=1}{\overset{2n-3}{\mathrm{\Σ}}}{\left(\right[a_q^{\mathrm{vid}},-1,b_j^{\mathrm{vid}}\left]{[u_{i,j},v_{i,j},1]}^T\right)}^2=\underset{q=1}{\overset{2n-3}{\mathrm{\Σ}}}{(a_q^{\mathrm{vid}}{u}_{i,j}+b_q^{\mathrm{vid}}-v_{i,j})}^2---\left(8\right)$

若m＝n-1时：

1)对水平方向上的m条平行直线组成的平行直线束建立第i行的直线最小二乘误差指标函数：

$E_{\mathrm{LS}}^{\mathrm{hor}}=\mathrm{\ϵ}(a_i^{\mathrm{hor}},b_i^{\mathrm{hor}})=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(\right[a_i^{\mathrm{hor}},-1,b_i^{\mathrm{hor}}\left]{[u_{i,j},v_{i,j},1]}^T\right)}^2=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(a_i^{\mathrm{hor}}{u}_{i,j}+b_i^{\mathrm{hor}}-v_{i,j})}^2---\left(9\right)$

2)对垂直方向上的n条平行直线组成的平行直线束建立第j列的直线最小二乘误差指标函数：

$E_{\mathrm{LS}}^{\mathrm{ver}}=\mathrm{\ϵ}(a_j^{\mathrm{ver}},b_j^{\mathrm{ver}})=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left(\right[a_j^{\mathrm{ver}},-1,b_j^{\mathrm{hor}}\left]{[u_{i,j},v_{i,j},1]}^T\right)}^2=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(a_j^{\mathrm{ver}}{u}_{i,j}+b_j^{\mathrm{ver}}-v_{i,j})}^2---\left(10\right)$

3)对主对角方向上的m+n-3条平行直线组成的平行直线束建立第k条的直线最小二乘误差指标函数：

$E_{\mathrm{LS}}^{\mathrm{mad}}=\mathrm{\ϵ}(a_k^{\mathrm{mad}},b_k^{\mathrm{mad}})=\underset{k=1}{\overset{m+n-3}{\mathrm{\Σ}}}{\left(\right[a_k^{\mathrm{mad}},-1,b_k^{\mathrm{mad}}\left]{[u_{i,j},v_{i,j},1]}^T\right)}^2=\underset{j=1}{\overset{m+n-3}{\mathrm{\Σ}}}{(a_k^{\mathrm{mad}}{u}_{i,j}+b_k^{\mathrm{mad}}-v_{i,j})}^2---\left(11\right)$

4)对副对角方向上的m+n-3条平行直线组成的平行直线束建立第q条的直线最小二乘误差指标函数：

$E_{\mathrm{LS}}^{\mathrm{vid}}=\mathrm{\ϵ}(a_q^{\mathrm{vid}},b_q^{\mathrm{vid}})=\underset{q=1}{\overset{2n-3}{\mathrm{\Σ}}}{\left(\right[a_q^{\mathrm{vid}},-1,b_j^{\mathrm{vid}}\left]{[u_{i,j},v_{i,j},1]}^T\right)}^2=\underset{q=1}{\overset{2n-3}{\mathrm{\Σ}}}{(a_q^{\mathrm{vid}}{u}_{i,j}+b_q^{\mathrm{vid}}-v_{i,j})}^2---\left(12\right)$

本实施例中，利用式(9)→(12)，计算水平、垂直、主对角线和副对角线可得到最佳的直线参数^horlⁱ[^horaⁱ,-1,^horbⁱ]，^verl^j[^vera^j,-1,^verb^j]，^madl^k[^mada^k,-1,^madb^k]，^vidl^q[^vida^q,-1,^vidb^q]，选取Image1图像其各个方向上的直线束参数见表4。

表4 Image1各个方向上的直线束参数

Step6：建立hor、ver、mad、vid方向的最佳消隐点目标函数：

$F_{\mathrm{opt}}^{\mathrm{hor}}=\underset{l_i^{\mathrm{hor}}}{\arg \min }\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{hor}},v_{\mathrm{va}}^{\mathrm{hor}}),l[a_i^{\mathrm{hor}},-1,b_i^{\mathrm{hor}}\left]\right]=\underset{l_i^{\mathrm{hor}}}{\arg \min }\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\frac{|a_i^{\mathrm{hor}}u_{\mathrm{va}}^{\mathrm{hor}}-v_{\mathrm{va}}^{\mathrm{hor}}+b_i^{\mathrm{hor}}|}{\sqrt{{\left(a_i^{\mathrm{hor}}\right)}^2+{(-1)}^2}}---\left(13\right)$

$F_{\mathrm{opt}}^{\mathrm{ver}}=\underset{l_j^{\mathrm{ver}}}{\arg \min }\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{ver}},v_{\mathrm{va}}^{\mathrm{ver}}),l[a_j^{\mathrm{ver}},-1,b_j^{\mathrm{ver}}\left]\right]=\underset{l_j^{\mathrm{ver}}}{\arg \min }\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{|a_j^{\mathrm{ver}}u_{\mathrm{va}}^{\mathrm{ver}}-v_{\mathrm{va}}^{\mathrm{ver}}+b_j^{\mathrm{ver}}|}{\sqrt{{\left(a_j^{\mathrm{ver}}\right)}^2+{(-1)}^2}}---\left(14\right)$

$F_{\mathrm{opt}}^{\mathrm{mad}}=\underset{l_k^{\mathrm{mad}}}{\arg \min }\underset{k=1}{\overset{m+n-3}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{mad}},v_{\mathrm{va}}^{\mathrm{mad}}),l[a_k^{\mathrm{mad}},-1,b_k^{\mathrm{mad}}\left]\right]=\underset{l_k^{\mathrm{hor}}}{\arg \min }\underset{k=1}{\overset{m+n-3}{\mathrm{\Σ}}}\frac{|a_k^{\mathrm{mad}}u_{\mathrm{va}}^{\mathrm{mad}}-v_{\mathrm{va}}^{\mathrm{mad}}+b_k^{\mathrm{mad}}|}{\sqrt{{\left(a_k^{\mathrm{mad}}\right)}^2+{(-1)}^2}}---\left(15\right)$

$F_{\mathrm{opt}}^{\mathrm{vid}}=\underset{l_q^{\mathrm{vid}}}{\arg \min }\underset{q=1}{\overset{m+n-3}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{vid}},v_{\mathrm{va}}^{\mathrm{vid}}),l[a_q^{\mathrm{vid}},-1,b_q^{\mathrm{vid}}\left]\right]=\underset{l_q^{\mathrm{vid}}}{\arg \min }\underset{q=1}{\overset{m+n-3}{\mathrm{\Σ}}}\frac{|a_q^{\mathrm{vid}}u_{\mathrm{va}}^{\mathrm{vid}}-v_{\mathrm{va}}^{\mathrm{vid}}+b_q^{\mathrm{vid}}|}{\sqrt{{\left(a_q^{\mathrm{vid}}\right)}^2+{(-1)}^2}}---\left(16\right)$

并应用Levenberg-Marquard算法对最佳消隐点目标函数寻优计算出最佳的消隐点坐标^gV＝[^gu_va,^gv_va]^T；

本实施例中，利用式(13)→(16)，计算每一幅图像上的4个消隐点坐标，附图7是Image1、Image2、Image3的消隐点坐标，坐标见表5。

表5 消隐点坐标(单位：pixel)

Step7：利用消隐线，验证消隐点计算的是否合理：

7.1对Step6中4个消隐点进行随机组合，生成种组合，并拟合出条直线，计算这些直线两两间的夹角，选择出最大的夹角θ_max；

7.2判断若θ≤θ_max，则这4个消隐点共线，形成了消隐线；反之，转到Step3，重新提取角点亚像素坐标，其中θ为阈值；

本实施例中，设定阈值θ＝0.01。这3幅图像上对应的4条直线的斜率见表6。从表中，可以观察出：最大误差在给定的阈值之内，计算出的消隐点共线。

表6 直线的斜率一览表

Step8：计算摄像机的内部参数f_x,f_y,u₀,v₀：

$v_0=-\frac{W_{23}}{W_{22}},\mathrm{\δ}=W_{33}-\frac{W_{13}^2-v_0{W}_{11}{W}_{23}}{W_{11}},f_x=\sqrt{\frac{\mathrm{\δ}}{W_{11}}},f_y=\sqrt{\frac{\mathrm{\δ}}{W_{22}}},u_0=-\frac{W_{13}{f}_x^2}{\mathrm{\δ}}---\left(17\right)$

Step8的具体实现方法是：

8.1建立摄像机内部参数几何约束方程组

$\left\{\begin{array}{c}\hat{V}_{\mathrm{hor}}\mathrm{\ω}\hat{V}_{\mathrm{ver}}=0\\ \hat{V}_{\mathrm{mad}}\mathrm{\ω}\hat{V}_{\mathrm{vid}}=0\end{array}\right.,\hat{V}_g={[u_{\mathrm{va}}^g,v_{\mathrm{va}}^g,1]}^T,$ g同Step4 (18)

其中：

$\mathrm{\ω}=K^{-T}{K}^{-1}=\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& -\frac{u_0}{f_x^2}\\ 0& \frac{1}{f_y^2}& -\frac{v_0}{f_y^2}\\ -\frac{u_0}{f_x^2}& -\frac{v_0}{f_y^2}& \frac{u_0^2}{f_x^2}+\frac{v_0^2}{f_y^2}+1\end{array}\right]=\left[\begin{array}{ccc}{W}_{11}& 0& W_{13}\\ 0& W_{22}& W_{23}\\ W_{13}& W_{23}& W_{33}\end{array}\right]---\left(19\right)$

f_x,f_y为CCD摄像机模型u轴和v轴上的归一化像素焦距，(u₀,v₀)是CCD摄像机模型成像平面的中心坐标；

8.2根据摄像机内部参数几何约束方程组，采用奇异值(SVD)分解和Zhang’s的方法计算该内部参数f_x,f_y,u₀,v₀；

本实施例中，利用式(17)、(18)、(19)，计算摄像机内部参数f_x,f_y,u₀,v₀，标定结果见表7。

表7 摄像机的内部参数标定结果

Step9：计算CCD摄像机模型中摄像机外部参数平移向量T的Z_C轴方向的分量t_z为：

t_z＝(t_z1+t_z2+t_z3)/3 (20)

式中：

$t_{z1}=\sqrt{({\mathrm{Sr}}_1-\sqrt{{{\mathrm{Sr}}_1}^2-{4\mathrm{Se}}_1})/{2\mathrm{Se}}_1},{\mathrm{Sr}}_1=b_1^2+b_2^2+b_4^2+b_5^2,{\mathrm{Se}}_1={(b_1{b}_5-b_2{b}_4)}^2---\left(21\right)$

$t_{z2}=\sqrt{({\mathrm{Sr}}_2-\sqrt{{{\mathrm{Sr}}_2}^2-{4\mathrm{Se}}_2})/{2\mathrm{Se}}_2},{\mathrm{Sr}}_2=b_4^2+b_5^2+b_7^2+b_8^2,{\mathrm{Se}}_2={(b_5{b}_7-b_4{b}_8)}^2---\left(22\right)$

$t_{z3}=\sqrt{({\mathrm{Sr}}_3-\sqrt{{{\mathrm{Sr}}_3}^2-{4\mathrm{Se}}_3})/{2\mathrm{Se}}_3},{\mathrm{Sr}}_3=b_1^2+b_2^2+b_7^2+b_8^2,{\mathrm{Se}}_3={(b_1{b}_8-b_2{b}_7)}^2---\left(23\right)$

上述Step9的具体实现方法是：

9.1建立摄像机外部参数的方程组：

$\left\{\begin{array}{c}{u}_{i,j}=f_x\frac{r_{1}X_{i,j}^W+r_{2}Y_{i,j}^W+t_x}{r_{7}X_{i,j}^W+r_{8}Y_{i,j}^W+t_z}+u_0\\ v_{i,j}=f_y\frac{r_{4}X_{i,j}^W+r_{5}Y_{i,j}^W+t_y}{r_{7}X_{i,j}^W+r_{8}Y_{i,j}^W+t_z}+v_0\end{array}\right.,i=\mathrm{1,2},...,m;j=\mathrm{1,2},...,n---\left(24\right)$

其中：摄像机外部参数的平移向量

$T=\left[\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right]$

摄像机外部参数的旋转矩阵

$R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right]$

^WX_i,j和^WY_i,j为棋盘格标定板中角点的空间坐标；

$\left[\begin{array}{cccccccc}{f}_{x}X_{i,j}^W& f_{x}Y_{i,j}^W& f_x& 0& 0& 0& {-U}_{x}X_{i,j}^W& {-U}_{x}Y_{i,j}^W\\ 0& 0& 0& f_{y}X_{i,j}^W& f_{y}Y_{i,j}^W& f_y& -V_{x}X_{i,j}^W& {-V}_{x}Y_{i,j}^W\end{array}\right]b=\left[\begin{array}{c}{U}_x\\ V_x\end{array}\right]---\left(25\right)$

式中：向量b＝[b₁ b₂ b₃ b₄ b₅ b₆ b₇ b₈]^T；U_x＝u_i,j-u₀；V_x＝v_i,j-v₀；

9.3应用最小二乘法计算9.2中矩阵形式，得到向量b＝[b₁ b₂ b₃ b₄ b₅ b₆ b₇ b₈]^T；进而得到t_z。

本实施例利用式(20)→(25)，计算得到t_z＝447.4675mm；

Step10：计算旋转矩阵R中的分量r₁,r₂,r₃,r₄,r₅,r₆,r₇,r₈,r₉和平移向量中的分量t_x,t_y；

上述Step10的具体实现方法是：

10.1向量b中的各分量为

b₁＝r₁/t_z,b₂＝r₂/t_z,b₃＝t_x/t_z,b₄＝r₄/t_z,b₅＝r₅/t_z,b₆＝t_y/t_z,b₇＝r₇/t_z,b₈＝r₈/t_z将t_z带入各分量便可求得r₁,r₂,r₄,r₅,r₇,r₈,t_x,t_y。

10.2对旋转矩阵R的前两列叉乘，边可得到r₃,r₆,r₉：

[r₃ r₆ r₉]^T＝[r₁ r₄ r₇]^T×[r₂ r₅ r₈]^T (26)

本实施例计算结果见表8：

表8 外部参数[R T]

Step11：引入一、二阶径向畸变系数k₁,k₂，对内部参数f_x,f_y,u₀,v₀和外部参数R,T进行非线性全局优化：

在摄像机的实际图像平面坐标系UV内建立图像残差当量最小化寻优函数：

式中：：3D空间点^WS从世界坐标系X_WY_WZ_W投影到物理图像平面坐标系UV的畸变2D点实际图像平面坐标；：3D空间点^WS从像素图像平面坐标系uv逆向投影到物理图像平面坐标系UV的畸变2D点实际图像坐标。

本实施例中引入畸变系数之后，利用式(27)，优化的内部参数、外部参数标定结果见表9-a和表9-b：

表9-a 内部参数优化结果

Claims (1)

1.一种基于正交消隐点的新型Tsai’s摄像机标定改进方法，其特征在于：采用以下步骤：

Step1：搭建光学测量系统，建立CCD摄像机模型；

Step2：使摄像机光轴与棋盘格标定板平面所成夹角在0°～80°之间，拍摄若干幅不同位姿的棋盘格标定板平面图像：并对拍摄的图像进行中值滤波、锐化预处理；棋盘格标定板的图案规格为m×n，其中m＝n或m＝n-1；

Step3：应用Harris角点提取算法对预处理后的图像提取角点亚像素坐标，角点亚像素坐标为s_i,j＝[u_i,j,v_i,j]^T,i＝1,2,...m,j＝1,2,...n；

Step4：将Step3中提取的角点划分成水平、垂直、主对角线、副对角线四组S＝{^gs|g＝hor,ver,mad,vid}；hor对应水平，ver对应垂直，mad对应主对角线，vid对应副对角线：

4.1若m＝n，则执行步骤4.2→4.6；若m＝n-1，则执行步骤4.7→4.11；

4.2定义棋盘格标定板的左下角第一个角点为s_1,1，水平向右依次为s_1,1,s_1,2,...,s_1,n，垂直向上依次为s_1,1,s_2,1,...,s_n,1；

4.3第一组：水平方向的角点坐标集合为

^hors＝{^hors¹,^hors²,...,^horsⁱ,...,^horsⁿ},i＝1,2,...,n，

其中子集合^horsⁱ＝{s_i,1,s_i,2,...s_i,j,...s_i,n},j＝1,2,...n；

4.4第二组：垂直方向的角点坐标集合为

^vers＝{^vers¹,^vers²,...,^vers^j,...,^versⁿ},j＝1,2,...,n，

其中子集合^vers^j＝{s_1,j,s_2,j,...s_i,j,...s_n,j},i＝1,2,...n；

4.5第三组：主对角线方向的角点坐标集合为

^mads＝{^mads¹,^mads²,...,^mads^k,...,^mads²ⁿ-³},k＝1,2,...,2n-3；

其中：

1)当k＝1时，子集合^mads¹＝{s_2,1,s_1,2}，表示第2行第1列的角点和第1行第2列的角点组成的集合；

2)当k＝2时，子集合^mads²＝{s_3,1,s_2,2,s_1,3}，表示第3行第1列、第2行第2列和第1行第3列的角点组成的集合；

3)当k＝n-1时，子集合^mads^n-1＝{s_n,1,s_n-1,2,...,s_n-i,j,...,s_1,n},i＝0,1,...,n-1；j＝1,2,...,n；

4)当k＝n时，子集合^madsⁿ＝{s_n,2,s_n-1,3,...,s_n-i,j,...,s_2,n},i＝0,1,...,n-2；j＝2,3,...,n；

5)当k＝2n-4时，子集合^mads^2n-4＝{s_n,n-2,s_n-1,n-1,s_n-2,n}；

6)当k＝2n-3时，子集合^mads^2n-3＝{s_n,n-1,s_n-1,n}；

4.6第四组：副对角方向的角点坐标集合为

^vids＝{^vids¹,^vids²,...,^vids^q,...,^vids^2n-3},q＝1,2,...,2n-3；

其中：

1)当q＝1时，子集合^vids¹＝{s_n-1,1,s_n,2}，表示第n-1行第1列的角点和第n行第2列的角点组成的集合；

2)当q＝2时，子集合^vids²＝{s_n-2,1,s_n-1,2,s_n,3}，表示第n-2行第1列、第n-1行第2列和第n行第3列的角点组成的集合；

3)当q＝n-1时，子集合^vids^n-1＝{s_1,1,s_2,2,...,s_i,j},i,j＝1,2,...,n；

4)当q＝n时，子集合^vidsⁿ＝{s_1,2,s_2,3,...,s_i,j},i＝1,2,...,n-1；j＝2,3,...,n；

5)当q＝2n-4时，子集合^vids^2n-4＝{s_1,n-2,s_2,n-1,s_3,n}；

6)当q＝2n-3时，子集合^vids^2n-3＝{s_1,n-1,s_2,n}；

4.7定义棋盘格标定板的左下角第一个角点为s_1,1，水平向右依次为s_1,1,s_1,2,...,s_1,n，垂直向上依次为s_1,1,s_2,1,...,s_m,1；

4.8第一组：水平方向的角点坐标集合为

^hors＝{^hors¹,^hors²,...,^horsⁱ,...,^hors^m},i＝1,2,...,m，

其中子集合^horsⁱ＝{s_i,1,s_i,2,...s_i,j,...s_i,n},j＝1,2,...n；

4.9第二组：垂直方向的角点坐标集合为

^vers＝{^vers¹,^vers²,...,^vers^j,...,^versⁿ},j＝1,2,...,n，

其中子集合^vers^j＝{s_1,j,s_2,j,...s_i,j,...s_m,j},i＝1,2,...m；

4.10第三组：主对角线方向的角点坐标集合为

^mads＝{^mads¹,^mads²,...,^mads^k,...,^mads^m+n-3},k＝1,2,...,m+n-3；

其中：

1)当k＝1时，子集合^mads¹＝{s_2,1,s_1,2}，表示第2行第1列的角点和第1行第2列的角点组成的集合；

2)当k＝2时，子集合^mads²＝{s_3,1,s_2,2,s_1,3}，表示第3行第1列、第2行第2列和第1行第3列的角点组成的集合；

3)当k＝m-1时，子集合^mads^m-1＝{s_m,1,s_m-1,2,...,s_m-i,j},i＝0,1,...,m-1；j＝1,2,...,n-1；

4)当k＝m时，子集合^mads^m＝{s_m,2,s_m-1,3,...,s_m-i,j},i＝0,1,...,m-1；j＝2,3,...,n；

5)当k＝m+n-4时，子集合^mads^m+n-4＝{s_m,n-2,s_m-1,n-1,s_m-2,n}；

6)当k＝m+n-3时，子集合^mads^m+n-3＝{s_m,n-1,s_m-1,n}；

4.11第四组：副对角方向的角点坐标集合为

^vids＝{^vids¹,^vids²,...,^vids^q,...,^vids^m+n-3},q＝1,2,...,m+n-3；

其中：

1)当q＝1时，子集合^vids¹＝{s_m-1,1,s_m,2}，表示第m-1行第1列的角点和第m行第2列的角点组成的集合；

2)当q＝2时，子集合^vids²＝{s_m-2,1,s_m-1,2,s_m,3}，表示第m-2行第1列、第m-1行第2列和第m行第3列的角点组成的集合；

3)当q＝n-1时，子集合^vids^n-1＝{s_1,1,s_2,2,...,s_i,j},i,j＝1,2,...,n；

4)当q＝n时，子集合^vidsⁿ＝{s_1,2,s_2,3,...,s_i,j},i＝1,2,...,m；j＝2,3,...,n；

5)当q＝m+n-4时，子集合^vids^m+n-4＝{s_1,n-2,s_2,n-1,s_3,n}；

6)当q＝m+n-3时，子集合^vids^m+n-3＝{s_1,n-1,s_2,n}；

Step5：应用Step4中角点的划分结果，采用最小二乘法计算hor、ver、mad、vid方向上的平行直线束参数：

若m＝n，则平行直线束参数为：

^horlⁱ[^horaⁱ,-1,^horbⁱ],i＝1,2,...,n

^verl^j[^vera^j,-1,^verb^j],j＝1,2,...,n

^madl^k[^mada^k,-1,^madb^k],k＝1,2,...,2n-3

^vidl^q[^vida^q,-1,^vidb^q],q＝1,2,...,2n-3

若m＝n-1，则平行直线束参数为：

^horlⁱ(^horaⁱ,-1,^horbⁱ),i＝1,2,...,m

^verl^j(^vera^j,-1,^verb^j),j＝1,2,...,n

^madl^k(^mada^k,-1,^madb^k),k＝1,2,...,m+n-3

^vidl^q(^vida^q,-1,^vidb^q),q＝1,2,...,m+n-3

Step6：建立hor、ver、mad、vid方向的最佳消隐点目标函数：

$F_{\mathrm{opt}}^{\mathrm{hor}}=\underset{l_i^{\mathrm{hor}}}{\arg \min }\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{hor}},v_{\mathrm{va}}^{\mathrm{hor}}),l[a_i^{\mathrm{hor}},-1,b_i^{\mathrm{hor}}\left]\right]=\underset{l_i^{\mathrm{hor}}}{\arg \min }\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\frac{|a_i^{\mathrm{hor}}u_{\mathrm{va}}^{\mathrm{hor}}-v_{\mathrm{va}}^{\mathrm{hor}}+b_i^{\mathrm{hor}}|}{\sqrt{{\left(a_i^{\mathrm{hor}}\right)}^2+{(-1)}^2}};$

$F_{\mathrm{opt}}^{\mathrm{ver}}=\underset{l_j^{\mathrm{ver}}}{\arg \min }\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{ver}},v_{\mathrm{va}}^{\mathrm{ver}}),l[a_j^{\mathrm{ver}},-1,b_j^{\mathrm{ver}}\left]\right]=\underset{l_j^{\mathrm{ver}}}{\arg \min }\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{|a_j^{\mathrm{ver}}u_{\mathrm{va}}^{\mathrm{ver}}-v_{\mathrm{va}}^{\mathrm{ver}}+b_j^{\mathrm{ver}}|}{\sqrt{{\left(a_j^{\mathrm{ver}}\right)}^2+{(-1)}^2}};$

$F_{\mathrm{opt}}^{\mathrm{mad}}=\underset{l_k^{\mathrm{mad}}}{\arg \min }\underset{k=1}{\overset{m+n-3}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{mad}},v_{\mathrm{va}}^{\mathrm{mad}}),l[a_k^{\mathrm{mad}},-1,b_k^{\mathrm{mad}}\left]\right]=\underset{l_k^{\mathrm{hor}}}{\arg \min }\underset{k=1}{\overset{m+n-3}{\mathrm{\Σ}}}\frac{|a_k^{\mathrm{mad}}u_{\mathrm{va}}^{\mathrm{mad}}-v_{\mathrm{va}}^{\mathrm{mad}}+b_k^{\mathrm{mad}}|}{\sqrt{{\left(a_k^{\mathrm{mad}}\right)}^2+{(-1)}^2}};$

$F_{\mathrm{opt}}^{\mathrm{vid}}=\underset{l_q^{\mathrm{vid}}}{\arg \min }\underset{q=1}{\overset{m+n-3}{\mathrm{\Σ}}}\mathrm{dis}[V(u_{\mathrm{va}}^{\mathrm{vid}},v_{\mathrm{va}}^{\mathrm{vid}}),l[a_q^{\mathrm{vid}},-1,b_q^{\mathrm{vid}}\left]\right]=\underset{l_q^{\mathrm{vid}}}{\arg \min }\underset{q=1}{\overset{m+n-3}{\mathrm{\Σ}}}\frac{|a_q^{\mathrm{vid}}u_{\mathrm{va}}^{\mathrm{vid}}-v_{\mathrm{va}}^{\mathrm{vid}}+b_q^{\mathrm{vid}}|}{\sqrt{{\left(a_q^{\mathrm{vid}}\right)}^2+{(-1)}^2}};$

并应用Levenberg-Marquard算法对最佳消隐点目标函数寻优计算出最佳的消隐点坐标^gV＝[^gu_va,^gv_va]^T；

Step7：利用消隐线，验证消隐点计算的是否合理：

7.1对Step6中4个消隐点进行随机组合，生成种组合，并拟合出条直线，计算这些直线两两间的夹角，选择出最大的夹角θ_max；

7.2判断若θ≤θ_max，则这4个消隐点共线，形成了消隐线；反之，转到Step3，重新提取角点亚像素坐标，其中θ为阈值；

Step8：建立摄像机内部参数几何约束方程组：

$\left\{\begin{array}{c}\hat{V}_{\mathrm{hor}}\mathrm{\ω}\hat{V}_{\mathrm{ver}}=0\\ \hat{V}_{\mathrm{mad}}\mathrm{\ω}\hat{V}_{\mathrm{vid}}=0\end{array}\right.,\hat{V}_g={[u_{\mathrm{va}}^g,v_{\mathrm{va}}^g,1]}^T$

其中：

$\mathrm{\ω}=K^{-T}{K}^{-1}=\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& -\frac{u_0}{f_x^2}\\ 0& \frac{1}{f_y^2}& -\frac{v_0}{f_y^2}\\ -\frac{u_0}{f_x^2}& -\frac{v_0}{f_y^2}& \frac{u_0^2}{f_x^2}+\frac{v_0^2}{f_y^2}+1\end{array}\right]=\left[\begin{array}{ccc}{W}_{11}& 0& W_{13}\\ 0& W_{22}& W_{23}\\ W_{13}& W_{23}& W_{33}\end{array}\right]$

f_x,f_y为CCD摄像机模型u轴和v轴上的归一化像素焦距，(u₀,v₀)是CCD摄像机模型成像平面的中心坐标；根据摄像机内部参数几何约束方程组，采用奇异值(SVD)分解和Zhang’s的方法计算该内部参数f_x,f_y,u₀,v₀；

Step9：计算CCD摄像机模型中摄像机外部参数平移向量T的Z_C轴方向的分量t_z为：

9.1建立摄像机外部参数的方程组：

$\left\{\begin{array}{c}{u}_{i,j}=f_x\frac{r_{1}X_{i,j}^W+r_{2}Y_{i,j}^W+t_x}{r_{7}X_{i,j}^W+r_{8}Y_{i,j}^W+t_z}+u_0\\ v_{i,j}=f_y\frac{r_{4}X_{i,j}^W+r_{5}Y_{i,j}^W+t_y}{r_{7}X_{i,j}^W+r_{8}Y_{i,j}^W+t_z}+v_0\end{array}\right.,i=\mathrm{1,2},...,m;j=\mathrm{1,2},...,n$

其中：摄像机外部参数的平移向量

$T=\left[\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right]$

摄像机外部参数的旋转矩阵

$R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right]$

^WX_i,j和^WY_i,j为棋盘格标定板中角点的空间坐标；

9.2将摄像机外部参数的方程组变形为矩阵形式为：

$\left[\begin{array}{cccccccc}{f}_{x}X_{i,j}^W& f_{x}Y_{i,j}^W& f_x& 0& 0& 0& {-U}_{x}X_{i,j}^W& {-U}_{x}Y_{i,j}^W\\ 0& 0& 0& f_{y}X_{i,j}^W& f_{y}Y_{i,j}^W& f_y& -V_{x}X_{i,j}^W& {-V}_{x}Y_{i,j}^W\end{array}\right]b=\left[\begin{array}{c}{U}_x\\ V_x\end{array}\right]$

式中：向量b＝[b₁ b₂ b₃ b₄ b₅ b₆ b₇ b₈]^T；U_x＝u_i,j-u₀；V_x＝v_i,j-v₀；

9.3应用最小二乘法计算9.2中矩阵形式，得到向量b＝[b₁ b₂ b₃ b₄ b₅ b₆ b₇ b₈]^T；并根据公式

t_z＝(t_z1+t_z2+t_z3)/3

$t_{z1}=\sqrt{({\mathrm{Sr}}_1-\sqrt{{{\mathrm{Sr}}_1}^2-{4\mathrm{Se}}_1})/{2\mathrm{Se}}_1},{\mathrm{Sr}}_1=b_1^2+b_2^2+b_4^2+b_5^2,{\mathrm{Se}}_1={(b_1{b}_5-b_2{b}_4)}^2$

$t_{z2}=\sqrt{({\mathrm{Sr}}_2-\sqrt{{{\mathrm{Sr}}_2}^2-{4\mathrm{Se}}_2})/{2\mathrm{Se}}_2},{\mathrm{Sr}}_2=b_4^2+b_5^2+b_7^2+b_8^2,{\mathrm{Se}}_2={(b_5{b}_7-b_4{b}_8)}^2$

$t_{z3}=\sqrt{({\mathrm{Sr}}_3-\sqrt{{{\mathrm{Sr}}_3}^2-{4\mathrm{Se}}_3})/{2\mathrm{Se}}_3},{\mathrm{Sr}}_3=b_1^2+b_2^2+b_7^2+b_8^2,{\mathrm{Se}}_3={(b_1{b}_8-b_2{b}_7)}^2$

得到t_z；

Step10：计算旋转矩阵R中的分量r₁,r₂,r₃,r₄,r₅,r₆,r₇,r₈,r₉和平移向量中的分量t_x,t_y；

Step11：引入一、二阶径向畸变系数k₁,k₂，对内部参数f_x,f_y,u₀,v₀和外部参数R,T进行非线性全局优化：

在摄像机的实际图像平面坐标系UV内建立图像残差当量最小化寻优函数：

式中：：3D空间点^WS从世界坐标系X_WY_WZ_W投影到物理图像平面坐标系UV的畸变2D点实际图像平面坐标；：3D空间点^WS从像素图像平面坐标系uv逆向投影到物理图像平面坐标系UV的畸变2D点实际图像坐标。