CN103426168B - 基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法 - Google Patents

Abstract

一种基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法，该方法有四大步骤：步骤一：摄像机内参的初始化；步骤二：摄像机外参的初始化；步骤三：优化目标函数的建立；步骤四：摄像机内外参的非线性优化。由于现有的标定算法不能够满足鱼眼镜头摄像机立体视觉系统的准确性与方便性，因此本发明即在此通用相机模型的基础上设计了对于普通、广角、鱼眼镜头摄像机通用的立体视觉系统标定方法。通过Matlab仿真和真实的实验可以验证本发明是可行的并且有效的。它在机器视觉技术领域里有较好的实用价值和应用前景。

Description

基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法

技术领域

本发明涉及一种基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法，该发明属于机器视觉领域。

背景技术

在机器视觉领域，立体视觉是重要的组成部分，被广泛应用于导航、三维测量及场景建模等领域。当前大部分的研究是集中于标定普通和广角立体摄像机系统，这些方法并不能直接用于标定鱼眼摄像机。另外，普通和广角摄像机系统只能提供比较有限的公共视场，而鱼眼立体摄像机系统的公共视场却能覆盖摄像机系统整个区域。因此，为了更好的发挥鱼眼摄像机的功能，提出一种精确、高效、方便的标定方法是非常重要和有意义的。

发明内容

本发明提供了一种基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法。它解决了鱼眼立体摄像机系统乃至多摄像机混合系统的标定复杂以及灵活性不高的问题。

本发明中采用的通用摄像机模型如下：

r(θ)=k₁θ+k₂θ³(1)

这一模型仅有两个参数k₁和k₂，由于r(θ)在[0，π/2]区间单调递增，因此模型可以反向求解。

如图1所示，一个三维空间点P通过一个鱼眼摄像机成像点为p，若通过针孔摄像机成像点则是p′，p′′为P在以O_c为球心的半球体的投影点。O_c-X_cY_cZ_c表示摄像机坐标系，o-xy表示图像坐标系。通过下式计算p的图像像素坐标为：

其中r(θ)可由式(1)获得，是径向方向与o-xy的x正轴方向的夹角，(x,y)为图像坐标（单位：mm），(u₀,v₀)是主点像素坐标，dx和dy分别是水平和垂直方向的像素大小。p′′点的坐标可表示为对应的准齐次坐标p′′′为：

对于任意空间点及其成像像素坐标可通过下列计算式求得θ和

$\left(\begin{array}{c}x\\ y\end{array}\right)={\left[\begin{array}{cc}1/\mathrm{dx}& 0\\ 0& 1/\mathrm{dy}\end{array}\right]}^{-1}\left(\begin{array}{c}u-u_0\\ v-v_0\end{array}\right)---\left(4\right)$

得到径向距离为：

$r\left(\mathrm{\θ}\right)=\sqrt{x^2+y^2}---\left(5\right)$

由式(1)通过卡丹公式求解一元三次方程组选择符合θ∈[0，π/2]的解。

径向夹角可通过反向求解获得

对于立体视觉而言，需要标定的参数包括各摄像机内参和摄像机间外参。摄像机内参即为摄像机模型中{k₁,k₂,u₀,v₀,dx,dy}六个参数，摄像机间外参即为摄像机坐标系之间的位置和姿态的关系，一般用3×3旋转矩阵R和3×1平移向量T表示：

$R=\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ R_{21}& R_{22}& R_{23}\\ R_{31}& R_{32}& R_{33}\end{array}\right],T={(t_x,t_y,t_z)}^T---\left(7\right)$

本发明提出了一种基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法，其基本思路为：首先利用摄像机产品说明书和摄像机成像特征点对信息计算出立体摄像机内外参数，然后利用计算得到的摄像机内外参数作为优化初始值，对建立的目标函数进行两步优化，从而得到摄像机内外参数的最优值。

本发明一种基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法，该方法具体步骤如下：

步骤一：摄像机内参的初始化

如上所述，每个摄像机的内参包括{k₁,k₂,u₀,v₀,dx,dy}六个参数。

1、主点(u₀,v₀)是光轴通过图像坐标系交点的坐标，一般可直接用图像中心点坐标初始化。

2、像素大小dx和dy可以摄像机生产商提供值作为初始值。

3、对于普通和广角摄像机，k₁和k₂可以通过摄像机生产商提供的标称焦距f和最大视角θ_max利用透视投影模型r=ftanθ来得到。具体实现方法为：将0～θ_max等分取m(m≥2)个θ_i(i=1,…,m)，分别代入透视投影模型r=ftanθ和(1)有

$f\tan \left({\mathrm{\θ}}_i\right)=k_1{\mathrm{\θ}}_i+k_2{\mathrm{\θ}}_i^3(i=1,\·\·\·m)---\left(8\right)$

此时，即得到关于k₁和k₂的m个线性方程，可通过线性最小二乘方法求得k₁和k₂。

$\left[\begin{array}{c}{k}_1\\ k_2\end{array}\right]={\left(A_t^T{A}_t\right)}^{-1}{A}_t^T{b}_r---\left(9\right)$

其中

$A_t=\left[\begin{array}{cc}{\mathrm{\θ}}_1& {\mathrm{\θ}}_1^3\\ \·& \·\\ \·& \·\\ \·& \·\\ {\mathrm{\θ}}_m& {\mathrm{\θ}}_m^3\end{array}\right],b_r=\left[\begin{array}{c}f\tan \left({\mathrm{\θ}}_1\right)\\ \·\\ \·\\ \·\\ f\tan \left({\mathrm{\θ}}_m\right)\end{array}\right]$

而对于鱼眼相机而言，则可令k₁=f,k₂=0。

步骤二：摄像机外参的初始化

摄像机外参的初始化可通过对摄像机系统的本质矩阵E进行奇异值分解得到。

1、对初始特征点进行归一化处理

给定n对两幅图像中的特征点集，根据式(4)～(6)以及各摄像机初始内参可得各特征点成像时对应的θ_i和通过式(3)计算出特征点在各摄像机模型的半球面投影点的准齐次坐标为其中X_i=(x_i,y_i,1)^T，X_i′=(x_i′,y_i′,1)^T。

令

$Y_i=\left[\begin{array}{c}{u}_i\\ v_i\\ 1\end{array}\right]=T_{\mathrm{norm}}{X}_i,{Y_i}^{\′}\left[\begin{array}{c}{u_i}^{\′}\\ {v_i}^{\′}\\ 1\end{array}\right]=T_{\mathrm{norm}}^{\′}{X}_i^{\′}---\left(10\right)$

其中，

$T_{\mathrm{norm}}=\left[\begin{array}{ccc}\sqrt{2}/s& 0& -\sqrt{2}\stackrel{\‾}{x}/s\\ 0& \sqrt{2}/s& -\sqrt{2}\stackrel{\‾}{y}/s\\ 0& 0& 1\end{array}\right]---\left(11\right)$

$\stackrel{\‾}{x}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}{x}_i,\stackrel{\‾}{y}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}{y}_i,$ $s=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[{(x_i-\stackrel{\‾}{x})}^2+{(y_i-\stackrel{\‾}{y})}^2]}.$ 同理，T′_norm也可类似求出。

2、计算本质矩阵

将归一化后的数据Y_i和Y_i′(i=1,…,n)代入

$A_e=\left[\begin{array}{ccccccccc}{u}_1^{\′}{u}_1& u_1^{\′}{v}_1& u_1^{\′}& v_1^{\′}{u}_1& v_1^{\′}{v}_1& v_1^{\′}& u_1& v_1& 1\\ \·& \·& \·& \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·& \·& \·& \·\\ u_n^{\′}{u}_n& u_n^{\′}{v}_n& u_n^{\′}& v_n^{\′}{u}_n& v_n^{\′}{v}_n& v_n^{\′}& u_n& v_n& 1\end{array}\right]---\left(12\right)$

对n×9的矩阵A_e进行奇异值分解得

$A_e=U_e\left(\begin{array}{c}{S}_e\\ 0_{(n-9)\×9}\end{array}\right)V_e^T---\left(13\right)$

其中S_e=diag(σ₁,…,σ₉)，且σ₁≥…≥σ₉，diag表示对角阵，0_(n-9)×9为(n-9)×9的全0矩阵。V_e为9×9矩阵，可由V_e的第9列元素得到：

$E_{\mathrm{norm}}=\left[\begin{array}{ccc}{V}_e\left(\mathrm{1,9}\right)& V_e\left(\mathrm{2,9}\right)& V_e\left(\mathrm{3,9}\right)\\ V_e\left(\mathrm{4,9}\right)& V_e\left(\mathrm{5,9}\right)& V_e\left(\mathrm{6,9}\right)\\ V_e\left(\mathrm{7,9}\right)& V_e\left(\mathrm{8,9}\right)& V_e\left(\mathrm{9,9}\right)\end{array}\right]---\left(14\right)$

采用奇异值分解将E_norm分解为

$E_{\mathrm{norm}}=U_n\left[\begin{array}{ccc}{D}_1& 0& 0\\ 0& D_2& 0\\ 0& 0& D_3\end{array}\right]V_n^T---\left(15\right)$

其中，D₁≥D₂≥D₃。令D₃=0，则有：

$E_{\mathrm{norm}}^{\′}=U_n\left[\begin{array}{ccc}{D}_1& 0& 0\\ 0& D_2& 0\\ 0& 0& 0\end{array}\right]V_n^T---\left(16\right)$

可得本质矩阵为：

$E={\left(T_{\mathrm{norm}}^{\′}\right)}^T{E}_{\mathrm{norm}}^{\′}{T}_{\mathrm{norm}}---\left(17\right)$

3、利用本质矩阵求取摄像机外参的可能解R′和T′

当通过式(17)获得本质矩阵E后，摄像机外参的可能解R′和T′即可通过本质矩阵E的奇异值分解来获得。令

E=USV^T(18)

此时R′和T′有四个可能的组合：

$(R^{\′},T^{\′})=\left\{\begin{array}{c}({\mathrm{UWV}}^T,U{\left(\mathrm{0,0,1}\right)}^T)\\ ({\mathrm{UWV}}^T,-U{\left(\mathrm{0,0,1}\right)}^T)\\ ({\mathrm{UW}}^T{V}^T,U{\left(\mathrm{0,0,1}\right)}^T)\\ ({\mathrm{UW}}^T{V}^T,U{\left(\mathrm{0,0,1}\right)}^T)\end{array}\right.---\left(19\right)$

这里

$W=\left[\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 1\end{array}\right]---\left(20\right)$

为了最终得到正确的组合，需遵循以下准则：1）R′的行列式值为1；2）每一个重建的三维点X必然位于两个摄像机的前方，也就是说，X的Z坐标必须全为正值。如因噪声影响不全为正，应取使得重建的X的Z坐标为正值最多的组合作为最终正确组合。

4、重建三维点得到真实解

假设三维点M在左右两个摄像机O₀-X₀Y₀Z₀坐标系和O₁-X₁Y₁Z₁坐标系的坐标分别为M₀=(X₀,Y₀,Z₀)^T和M₁=RM₀+T，且点在两相机的成像像素坐标点为(u_m,v_m)和(u′_m,v′_m)，由式(4)～(6)及初始内参可求得对应的θ₀,θ₁,

此时即可利用线性最小二乘方法可直接求解出X₀,Y₀,Z₀：

$\left[\begin{array}{c}{X}_0\\ Y_0\\ Z_0\end{array}\right]={\left(A_r^T{A}_r\right)}^{-1}{A}_r^T{b}_x---\left(21\right)$

其中

因此可将(19)中四组可能解分别替换(R,T)代入式(21)中进行三维点重建计算，利用判断准则对可能解进行筛选，直到得到最终的唯一解，此时有R=R′。由于||T′||=1，因此计算得到的平移向量T'与真实的平移向量T相差一个常值因子。令L′表示两个重建三维特征点的距离，L表示真实距离。为了最小化计算误差，常值因子λ需要使用两个相机所有的图像点对计算得到，此时有

$\mathrm{\λ}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{L}{L^{\′}}---\left(22\right)$

其中N表示图像点对的数量。最终求得的转移向量为：

步骤三：目标函数的建立

本发明中一维标定杆上特征点的个数并没有限制，只要能保证每一对图像上同一三维点的不同成像点能够准确匹配就能完成标定。在本发明中在一维标定杆上装载了三个特征点，如图2所示，若因应用需要添加特征点，依次类推计算即可。由于噪声，一维标定杆上的特征点重建以后会产生误差。令d(A,B)=||A-B||（||·||表示欧几里得向量范数），得到距离和共线误差可表示为：

$g_1^k(x,y)=|L_1-d(A_k^r,B_k^r)|---\left(24\right)$

$g_2^k(x,y)=|L_2-d(B_k^r,C_k^r)|---\left(25\right)$

$g_3^k(x,y)=|L-d(A_k^r,C_k^r)|---\left(26\right)$

$g_4^k(x,y)=|d(A_k^r,C_k^r)-d(A_k^r,B_k^r)-d(B_k^r,C_k^r)|---\left(27\right)$

其中 k₁,k₂,u₀,v₀,dx,dy和k₁′,k₂′,u₀′,v₀′,dx′,dy′分别为左右摄像机的内部参数；可以通过Rodrigues公式与旋转矩阵R互相转化，即 $R=e^{{\left[r\right]}_{\×}}=I+\sin \mathrm{\ϵ}{\left[\hat{r}\right]}_{\×}+(1-\cos \mathrm{\ϵ}){\left[\hat{r}\right]}_{\×}^2,$ $r=\mathrm{\ϵ}\hat{r},\mathrm{\ϵ}=\left|\right|r\left|\right|,$ 为单位向量；为第k幅图像对中通过(21)计算的A,B,C的重建点；L,L₁,L₂分别为真实特征点AC、AB及BC之间的实际距离。

假设a_0,k,b_0,k,c_0,k和a_1,k,b_1,k,c_1,k是采集第k幅图像点对时通过式(3)求得的三维特征点A,B,C在左右摄像机半球面模型上投影点对应的准齐次坐标，有：

$g_5^k(x,y)=\left|a_{1,k}^T{\left[T\right]}_{\×}{\mathrm{Ra}}_{0,k}\right|---\left(28\right)$

$g_6^k(x,y)=\left|b_{1,k}^T{\left[T\right]}_{\×}{\mathrm{Rb}}_{0,k}\right|---\left(29\right)$

$g_7^k(x,y)=\left|c_{1,k}^T{\left[T\right]}_{\×}{\mathrm{Rc}}_{0,k}\right|---\left(30\right)$

其中，[T]_×为由T定义的反对称矩阵：

${\left[T\right]}_{\×}=\left[\begin{array}{ccc}0& -t_z& t_y\\ t_z& 0& -t_x\\ -t_y& t_x& 0\end{array}\right]$

结合式(24)～(30)，得到的优化目标函数为

$g(x,y)=\min ({\mathrm{\ρ}}_1\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(x,y)}^2+g_2^k{(x,y)}^2+g_3^k{(x,y)}^2)+$ (31)

${\mathrm{\ρ}}_2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}(g_4^k{(x,y)}^2+g_5^k{(x,y)}^2+g_6^k{(x,y)}^2+g_7^k{(x,y)}^2))$

其中ρ₁∈(800,1200),ρ₂∈(0.5,1.5)。

步骤四：内外参数的非线性优化

由步骤三知，共有18个参数需要优化，可以将优化过程分为两步：

1、固定y为y₀，即y=y₀，y₀可由步骤一的初始值给定。设定x的估计为初始值，利用LM（Levenberg-Marquardt）方法来求解如下优化问题

$\underset{x}{\min g}(x,y_0)---\left(32\right)$

得到最优解x^*。

2、设定y₀为初始值，利用LM求解如下优化问题

$\underset{y}{\min g}(x^{*},y)---\left(33\right)$

得到最优解y^*。

此时，即得到了摄像机内外参数的最优解x^*和y^*。

其中，在步骤四中所述的“LM方法”是一种使用广泛的非线性最小二乘算法，中文为列文伯格-马夸尔特法，即Levenberg-Marquardt法。

优点与功效：

本发明一种基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法，其优点在于，解决了当前标定算法操作的复杂性以及针对鱼眼相机标定的不精确性等问题。本发明通过采用一维标定杆作为标定特征物以及通用的相机模型，使得在标定中实现了对不同相机进行同时标定，大大减小了实际应用中标定的工作量，同时具有精度高，操作方便，简单易行等特点。

附图说明

图1：鱼眼相机模型

图2：球形立体视觉示意图

图3：不同视觉系统标定仿真结果图

图4：实际试验使用的一维标定杆图

图5：为本发明所述方法流程图

图中符号说明如下：

图1中的符号说明：O表示图像坐标系原点，x,y为图像坐标系的坐标轴。O_c表示摄像机坐标系的原点，X_c,Y_c,Z_c为摄像机坐标系的坐标轴。P为空间一个三维点，p为P在鱼眼摄像机下的成像点，p′为P在针孔摄像机下的成像点，p′′为P在以O_c为球心的半球体上的投影点。f为摄像机的焦距，r为鱼眼摄像机成像点p距离主点的径向距离。θ为入射线与摄像机光轴之间的夹角，是径向距离r与图像坐标系x正轴方向的夹角。

图2中的符号说明：A,B,C（A_k,B_k,C_k分别表示它们在第k幅图像对中的位置）分别表示一维标定杆上的三个共线特征点；O₀-X₀Y₀Z₀和O₁-X₁Y₁Z₁分别表示左右两个摄像机所处的摄像机坐标系，(R,T)分别为两摄像机之间的旋转矩阵和平移向量；a_0,k,b_0,k,c_0,k和a_1,k,b_1,k,c_1,k分别表示三维空间点A_k,B_k,C_k分别投影到以O₀为球心的单位半球体上的和以O₁为球心的单位半球体上的投影点。

具体实施方式

本发明一种基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法，见图1～图5所示。

仿真过程是在主频3.07GHz、内存4.00GB的计算机上，WindowsXP环境下的Matlab2009b上进行的。实验对象采用的是BaslerscA640-120gc相机，视场角为185°的FujinonFE185C057HA-1鱼眼镜头，视场角为86.77°的PentaxC60402KP广角镜头。

首先介绍一下仿真中立体视觉摄像机的各个参数，这些参数根据真实摄像机参数假设。采集图像对数量N=300，图像主点u₀=u₀′=320,v₀=v₀′=240，像素大小dx=dx′=0.0056mm，dy=dy′=0.0056mm。普通镜头的视场角为45°，标称焦距为5.5mm；广角镜头的视场角为100°，标称焦距为2.2mm；鱼眼镜头的视场角为185°，标称焦距为1.8mm。仿真实验分为四组不同的立体视觉仿真：1）普通镜头立体摄像机组；2）广角镜头立体摄像机组；3）鱼眼镜头立体摄像机组；4）广角镜头与鱼眼镜头混合组成的立体摄像机组；

（1）本发明一种基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法，见图5，该方法具体步骤如下：

步骤一：摄像机内参的初始化

根据仿真参数，所有系统中u₀=u₀′=320,v₀=v₀′=240，dx=dx′=0.0056，dy=dy′=0.0056；

对于普通相机，根据式(9)得k₁=5.3727,k₂=2.6344；对于广角相机，根据式(9)得k₁=1.4586,k₂=2.2714；

对于鱼眼相机，根据仿真设置有k₁=f=1.8,k₂=0；

步骤二：摄像机外参的初始化

利用随机序列产生N对图像对仿真实际过程中摄像机采集到的一维标定杆特征点，根据该图像对数据和初始化内参利用式(4)～(6)计算得到对应的θ和，其次利用(10)～(17)求得本质矩阵E，通过式(18)～(21)对本质矩阵E进行奇异值分解并利用立体视觉特征筛选出正确的摄像机内外参数(R′,T′)，此时有R=R′，通过式(22)～(23)最小化计算误差，得到平移向量T；

步骤三：目标函数的建立

利用式(24)～(30)计算各类误差项，其中L₁=200mm，L₂=400mm，L=600mm，建立的目标函数为：

$g(x,y)=\min \left({\mathrm{\ρ}}_1\underset{k=1}{\overset{N}{\mathrm{\Σ}}}(g_1^k{(x,y)}^2+g_2^k{(x,y)}^2+g_3^k{(x,y)}^2\right)+$ (34)

${\mathrm{\ρ}}_2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}(g_4^k{(x,y)}^2+g_5^k{(x,y)}^2+g_6^k{(x,y)}^2+g_7^k{(x,y)}^2))$

其中ρ₁=1000，ρ₂=1；

步骤四：内外参数的非线性优化

给定上述图像对数据以及内外参数的初始值，对步骤三中的目标函数进行两步优化，直接使用matlab的lsqnonlin函数，得到最终的最优内外参数x^*和y^*。

（2）仿真结果分析

仿真过程中，一维标定杆在每个立体视觉系统中随机挥舞300次，得到了N=300组图像对数据。为了描述偏差，在得到的图像点中添加均值μ=0、标准差σ=0～3个像素的高斯白噪声，每次试验执行10次去均值作为结果。图3显示了不同的立体视觉系统的标定结果，这里RMS表示线段AC的真实值与所有图像对重建值的均方根误差。

由图可知，即使当噪声标准差σ=3个像素时，RMS误差依然小于4%，表明本发明提出的算法有高的准确性和稳定性。

（3）实物试验验证

为了进一步验证本发明的可行性，我们在真实的立体视觉系统中进行了测试。真实试验主要完成了三个不同的立体视觉系统：1）鱼眼镜头立体视觉系统；2）广角镜头立体视觉系统；3）鱼眼和广角镜头混合立体视觉系统。试验采用的一维标定杆如图4所示，它是一根镶有三个共线LED灯的空心棒，这三个LED灯之间的距离分别为600mm，400mm和200mm。同样以实际标志点距离AC与标定后重建距离作为衡量标准，不同立体视觉系统的实际实验结果如下表1‐3：

表1.鱼眼镜头立体视觉系统中测量值与实际值的比较

表2.广角镜头立体视觉系统中测量值与实际值的比较

表3.混合镜头立体视觉系统中测量值与实际值的比较

Claims (1)

1.一种基于一维标定杆的普通、广角、鱼眼立体摄像机的通用标定方法，其特征在于：该方法具体作法如下：

本通用标定方法所采用的通用摄像机模型如下：

r(θ)＝k₁θ+k₂θ³(1)

这一模型仅有两个参数k₁和k₂，由于r(θ)在[0,π/2]区间单调递增，因此模型是反向求解；

一个三维空间点P通过一个鱼眼摄像机成像点为p，若通过针孔摄像机成像点则是p′，p″为P在以O_c为球心的半球体的投影点，O_c-X_cY_cZ_c表示摄像机坐标系，o-xy表示图像坐标系，通过下式计算p的图像像素坐标为：

其中r(θ)由式(1)获得，是径向方向与o-xy的x正轴方向的夹角，(x,y)为图像坐标，(u₀,v₀)是主点像素坐标，dx和dy分别是水平和垂直方向的像素大小；p″点的坐标表示为对应的准齐次坐标p″′为：

对于任意空间点及其成像像素坐标通过下列计算式求得θ和

$\left(\begin{array}{c}x\\ y\end{array}\right)={\left[\begin{array}{cc}1/dx& 0\\ 0& 1/dy\end{array}\right]}^{-1}\left(\begin{array}{c}u-u_0\\ v-v_0\end{array}\right)---\left(4\right)$

得到径向距离为：

$r\left(\mathrm{\θ}\right)=\sqrt{x^2+y^2}---\left(5\right)$

由式(1)通过卡丹公式求解一元三次方程组选择符合θ∈[0,π/2]的解；

径向夹角通过反向求解获得

对于立体视觉而言，需要标定的参数包括各摄像机内参和摄像机间外参，摄像机内参即为摄像机模型中{k₁,k₂,u₀,v₀,dx,dy}六个参数，摄像机间外参即为摄像机坐标系之间的位置和姿态的关系，一般用3×3旋转矩阵R和3×1平移向量T表示：

$R=\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ R_{21}& R_{22}& R_{23}\\ R_{31}& R_{32}& R_{33}\end{array}\right],T={(t_x,t_y,t_z)}^T---\left(7\right)$

该方法具体步骤如下：

步骤一：摄像机内参的初始化

每个摄像机的内参包括{k₁,k₂,u₀,v₀,dx,dy}六个参数；

1、主点(u₀,v₀)是光轴通过图像坐标系交点的坐标，一般直接用图像中心点坐标初始化；

2、像素大小dx和dy以摄像机生产商提供值作为初始值；

3、对于普通和广角摄像机，k₁和k₂通过摄像机生产商提供的标称焦距f和最大视角θ_max利用透视投影模型r＝ftanθ来得到；具体实现方法为：将0～θ_max等分取m(m≥2)个θ_i(i＝1,…,m)，分别代入透视投影模型r＝ftanθ和式(1)，得到

$ftan\left({\mathrm{\θ}}_i\right)=k_1{\mathrm{\θ}}_i+k_2{\mathrm{\θ}}_i^3(i=1,...,m)---\left(8\right)$

此时，即得到关于k₁和k₂的m个线性方程，通过线性最小二乘方法求得k₁和k₂；

$\left[\begin{array}{c}{k}_1\\ k_2\end{array}\right]={\left(A_t^T{A}_t\right)}^{-1}{A}_t^T{b}_r---\left(9\right)$

其中

$A_t=\left[\begin{array}{cc}{\mathrm{\θ}}_1& {\mathrm{\θ}}_1^3\\ \·& \·\\ \·& \·\\ \·& \·\\ {\mathrm{\θ}}_m& {\mathrm{\θ}}_m^3\end{array}\right],b_r=\left[\begin{array}{c}ftan\left({\mathrm{\θ}}_1\right)\\ \·\\ \·\\ \·\\ ftan\left({\mathrm{\θ}}_m\right)\end{array}\right]$

而对于鱼眼相机而言，则令k₁＝f,k₂＝0；

步骤二：摄像机外参的初始化

摄像机外参的初始化通过对摄像机系统的本质矩阵E进行奇异值分解得到；

1、对初始特征点进行归一化处理

给定n对两幅图像中的特征点集，根据式(4)～(6)以及各摄像机初始内参得各特征点成像时对应的θ_i和通过式(3)计算出特征点在各摄像机模型的半球面投影点的准齐次坐标为其中X_i＝(x_i,y_i,1)^T，X′_i＝(x′_i,y′_i,1)^T；

令

$Y_i=\left[\begin{array}{c}{u}_i\\ v_i\\ 1\end{array}\right]=T_{norm}{X}_i,Y_i^{\′}=\left[\begin{array}{c}{u}_i^{\′}\\ v_i^{\′}\\ 1\end{array}\right]=T_{norm}^{\′}{X}_i^{\′}---\left(10\right)$

其中，

$T_{norm}=\left[\begin{array}{ccc}\sqrt{2}/s& 0& -\sqrt{2}\stackrel{\‾}{x}/s\\ 0& \sqrt{2}/s& -\sqrt{2}\stackrel{\‾}{y}/s\\ 0& 0& 1\end{array}\right]---\left(11\right)$

$\stackrel{\‾}{x}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}{x}_i,\stackrel{\‾}{y}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}{y}_i,s=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\[{(x_i-\stackrel{\‾}{x})}^2\mathrm{+(}{y}_i-\stackrel{\‾}{y}{)}^2\]};$ 同理，T′_norm也类似求出；

2、计算本质矩阵

将归一化后的数据Y_i和Y′_i(i＝1,…,n)代入

$A_e=\left[\begin{array}{ccccccccc}{u}_1^{\′}{u}_1& u_1^{\′}{v}_1& u_1^{\′}& v_1^{\′}{u}_1& v_1^{\′}{v}_1& v_1^{\′}& u_1& v_1& 1\\ \·& \·& \·& \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·& \·& \·& \·\\ u_n^{\′}{u}_n& u_n^{\′}{v}_n& u_n^{\′}& v_n^{\′}{u}_n& v_n^{\′}{v}_n& v_n^{\′}& u_n& v_n& 1\end{array}\right]---\left(12\right)$

对n×9的矩阵A_e进行奇异值分解得

$A_e=U_e\left(\begin{array}{c}{S}_e\\ 0_{(n-9)\×9}\end{array}\right)V_e^T---\left(13\right)$

其中S_e＝diag(σ₁,…,σ₉)，且σ₁≥…≥σ₉，diag表示对角阵，0_(n-9)×9为(n-9)×9的全0矩阵；V_e为9×9矩阵，由V_e的第9列元素得到：

$E_{norm}=\left[\begin{array}{ccc}{V}_e(1,9)& V_e(2,9)& V_e(3,9)\\ V_e(4,9)& V_e(5,9)& V_e(6,9)\\ V_e(7,9)& V_e(8,9)& V_e(9,9)\end{array}\right]---\left(14\right)$

采用奇异值分解将E_norm分解为

$E_{norm}=U_n\left[\begin{array}{ccc}{D}_1& 0& 0\\ 0& D_2& 0\\ 0& 0& D_3\end{array}\right]V_n^T---\left(15\right)$

其中，D₁≥D₂≥D₃，令D₃＝0，则有：

$E_{norm}^{\′}=U_n\left[\begin{array}{ccc}{D}_1& 0& 0\\ 0& D_2& 0\\ 0& 0& 0\end{array}\right]V_n^T---\left(16\right)$

得到本质矩阵为：

E＝(T′_norm)^TE′_normT_norm(17)

3、利用本质矩阵求取摄像机外参的可能解R′和T′

当通过式(17)获得本质矩阵E后，摄像机外参的可能解R′和T′即通过本质矩阵E的奇异值分解来获得；令

E＝USV^T(18)

此时R′和T′有四个组合：

$(R^{\′},T^{\′})=\left\{\begin{array}{c}({\mathrm{UWV}}^T,U{(0,0,1)}^T)\\ ({\mathrm{UWV}}^T,-U{(0,0,1)}^T)\\ ({\mathrm{UW}}^T{V}^T,U{(0,0,1)}^T)\\ ({\mathrm{UW}}^T{V}^T,U{(0,0,1)}^T)\end{array}\right.---\left(19\right)$

这里

$W=\left[\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 1\end{array}\right]---\left(20\right)$

为了最终得到正确的组合，需遵循以下准则：1)R′的行列式值为1；2)每一个重建的三维点X必然位于两个摄像机的前方，也就是说，X的Z坐标必须全为正值，如因噪声影响不全为正，应取使得重建的X的Z坐标为正值最多的组合作为最终正确组合；

4、重建三维点得到真实解

假设三维点M在左右两个摄像机O₀-X₀Y₀Z₀坐标系和O₁-X₁Y₁Z₁坐标系的坐标分别为M₀＝(X₀,Y₀,Z₀)^T和M₁＝RM₀+T，且点在两相机的成像像素坐标点为(u_m,v_m)和(u′_m,v′_m)，由式(4)～(6)及初始内参求得对应的θ₀,θ₁,

此时利用线性最小二乘方法直接求解出X₀,Y₀,Z₀：

$\left[\begin{array}{c}{X}_0\\ Y_0\\ Z_0\end{array}\right]={\left(A_r^T{A}_r\right)}^{-1}{A}_r^T{b}_x---\left(21\right)$

其中

因此将式(19)中四组可能解分别替换(R,T)代入式(21)中进行三维点重建计算，利用判断准则对可能解进行筛选，直到得到最终的唯一解，此时有R＝R′；由于||T′||＝1，因此计算得到的平移向量T'与真实的平移向量T相差一个常值因子；令L′表示两个重建三维特征点的距离，L表示真实距离；为了最小化计算误差，常值因子λ需要使用两个相机所有的图像点对计算得到，此时有

$\mathrm{\λ}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}\frac{L}{L^{\′}}---\left(22\right)$

其中N表示图像点对的数量，最终求得的转移向量为：

步骤三：目标函数的建立

本方法中一维标定杆上特征点的个数并没有限制，只要能保证每一对图像上同一三维点的不同成像点能够准确匹配就能完成标定；本方法中在一维标定杆上装载了三个特征点，若因应用需要添加特征点，依次类推计算即行；由于噪声，一维标定杆上的特征点重建以后会产生误差，令d(A,B)＝||A-B||，||·||表示欧几里得向量范数，得到距离和共线误差表示为：

$g_1^k(x,y)=|L_1-d(A_k^r,B_k^r)|---\left(24\right)$

$g_2^k(x,y)=|L_2-d(B_k^r,C_k^r)|---\left(25\right)$

$g_3^k(x,y)=|L-d(A_k^r,C_k^r)|---\left(26\right)$

$g_4^k(x,y)=|d(A_k^r,C_k^r)-d(A_k^r,B_k^r)-d(B_k^r,C_k^r)|---\left(27\right)$

其中k₁,k₂,u₀,v₀,dx,dy和k′₁,k′₂,u′₀，v′₀,dx′,dy′分别为左右摄像机的内部参数；通过Rodrigues公式与旋转矩阵R互相转化，即 $R=e^{{\[r\]}_{\×}}=I+sin\mathrm{\ϵ}{\[\hat{r}\]}_{\×}+(1-cos\mathrm{\ϵ}){\[\hat{r}\]}_{\×}^2,$ $r=\mathrm{\ϵ}\hat{r},\mathrm{\ϵ}=\left|\right|r\left|\right|,$ 为单位向量；为第k幅图像对中通过式(21)计算的A,B,C的重建点；L,L₁,L₂分别为真实特征点AC、AB及BC之间的实际距离；

假设a_0,k,b_0,k,c₀,_k和a_1,k,b_1,k,c_1,k是采集第k幅图像点对时通过式(3)求得的三维特征点A,B,C在左右摄像机半球面模型上投影点对应的准齐次坐标，有：

$g_5^k(x,y)=\left|a_{1,k}^T{\[T\]}_{\×}{\mathrm{Ra}}_{0,k}\right|---\left(28\right)$

$g_6^k(x,y)=\left|b_{1,k}^T{\[T\]}_{\×}{\mathrm{Rb}}_{0,k}\right|---\left(29\right)$

$g_7^k(x,y)=\left|c_{1,k}^T{\[T\]}_{\×}{\mathrm{Rc}}_{0,k}\right|---\left(30\right)$

其中，[T]_×为由T定义的反对称矩阵：

${\[T\]}_{\×}=\left[\begin{array}{ccc}0& -t_z& t_y\\ t_z& 0& -t_x\\ -t_y& t_x& 0\end{array}\right]$

结合式(24)～(30)，得到的优化目标函数为

$\begin{array}{c}g(x,y)=\min \left({\mathrm{\ρ}}_1\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\right(g_1^k{(x,y)}^2+g_2^k{(x,y)}^2+g_3^k{(x,y)}^2)+\\ {\mathrm{\ρ}}_2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}(g_4^k{(x,y)}^2+g_5^k{(x,y)}^2+g_6^k{(x,y)}^2+g_7^k{(x,y)}^2))\end{array}---\left(31\right)$

其中ρ₁∈(800,1200),ρ₂∈(0.5,1.5)；

步骤四：内外参数的非线性优化

由步骤三知，共有18个参数需要优化，将优化过程分为两步：

1、固定y为y₀，即y＝y₀，y₀由步骤一的初始值给定，设定x的估计为初始值，利用LM方法来求解如下优化问题

$\underset{x}{\min }g(x,y_0)---\left(32\right)$

得到最优解x^*；

2、设定y₀为初始值，利用LM求解如下优化问题

$\underset{y}{\min }g(x^{*},y)---\left(33\right)$

得到最优解y^*；

此时，即得到了摄像机内外参数的最优解x^*和y^*。