0 概述

随着计算机技术的发展, 计算机视觉被广泛应用于各领域, 其得以实施的前提是通过相机标定(即相机内外参数的求解过程^[1])建立相机坐标系到真实世界坐标系的映射。传统标定方法利用靶标上的控制点来辅助标定。文献[2]提出一种由3个正交平面组成的高精度立体靶标进行相机标定的方法, 但由于成像模型过于简单, 效果不理想。文献[3]提出一种自标定方法, 根据相机参数本身的约束关系进行标定。相比传统标定方法, 自标定方法不依赖于靶标, 但要求相机进行刚体运动, 而且精度相对传统标定方法差。文献[4]提出张正友标定法, 使用单平面棋盘格作为靶标, 制作简便, 精度明显优于自标定方法, 且对普通相机可以保证较高的精度, 是当时的主流方法。

但张正友标定法通过LM算法^[5]求解相机参数, 因此存在传统非线性优化方法共有的局限性, 即初值设置不当容易使其陷入局部最优。针对该问题, 文献[6]提出自适应阻尼因子的LM算法, 文献[7]使用改进的遗传算法代替LM算法对内参数求解进行优化, 文献[8]在数据归一化前对成像中心附近的点使用非线性优化方法预先优化参数。另外, 张正友标定法需要采集标定板不同方位的多张图片, 无法应用于运动受限的场景, 对此, 一些学者进行了相关的改进, 如文献[9]使用立体靶标的LM迭代优化, 实现单张图像的标定。

随着全景视觉的发展, 鱼眼镜头得到了广泛应用^[10-11]。鱼眼镜头视场角一般超过140°甚至180°^[12], 具有比普通镜头大得多的径向畸变。传统多项式畸变模型对较大畸变的描述效果不好, 因此文献[13]引入除法模型^[14]对鱼眼镜头进行自标定, 相比其他模型, 除法模型具备描述较大畸变的能力^[14], 但由于自标定方法受场景、运动等因素影响, 精度较低, 难以应用在高精度测量中。

针对鱼眼镜头高精度标定的需求, 本文提出一种基于除法模型的标定方法, 利用立体标靶, 并结合畸变中心的特性, 避免局部极小值的干扰, 实现快速鲁棒的标定。

1 相机成像模型

通常真实镜头引入的畸变无法用透视投影来建模^[15], 因此为准确描述相机成像模型, 本文采用如图 1所示的针孔透视模型与畸变除法模型对相机成像进行描述。

在图 1中, O_cX_cY_cZ_c为相机(camera)空间坐标系, Oxy为像平面坐标系, O_wX_wY_wZ_w为真实世界(world)三维坐标系, $\widehat O\widehat x\widehat y$为原点移至畸变中心的像平面坐标系。文献[16-17]实验证明畸变中心可能会明显偏离相机的成像中心或者主点, 因此在该模型中畸变中心$\widehat O$和像主点O_p不一定重合。

令P_w为观测点在真实世界的三维坐标, 则P_u为P_w在像平面上的无畸变的(undistorted)透视投影, P_d为P_w在像平面上的带畸变的(distorted)真实投影。

设P_w的坐标为V_w=(X_w  Y_w  Z_w  1)^T, P_u的坐标为v_u=(x_u  y_u  1)^T, 则有:

$ \rho {\mathit{\boldsymbol{v}}_{\rm{u}}} = \mathit{\boldsymbol{A}}\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{R}}&\mathit{\boldsymbol{t}} \end{array}} \right]{\mathit{\boldsymbol{V}}_{\rm{w}}} $

(1)

其中, ρ为放缩因子, A为3×3的内参数矩阵, R为3×3的旋转矩阵, t为3×1的平移向量。

令投影矩阵为:

$ \mathit{\boldsymbol{P}} = \mathit{\boldsymbol{A}}\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{R}}&\mathit{\boldsymbol{t}} \end{array}} \right] $

(2)

式(1)可简化为:

$ \rho {\mathit{\boldsymbol{v}}_{\rm{u}}} = \mathit{\boldsymbol{P}}{\mathit{\boldsymbol{V}}_{\rm{w}}} $

(3)

设P_d的坐标为v_d=(x_d  y_d  1)^T, 畸变中心$\widehat O$的坐标为e=(du₀  dv₀  1)^T, 则真实像点P_d到$\widehat O$的径向(radial)距离为:

$ {r_{\rm{d}}} = \sqrt {{{\left( {{x_{\rm{d}}} - d{u_0}} \right)}^2} + {{\left( {{y_{\rm{d}}} - d{v_0}} \right)}^2}} $

(4)

引入除法模型描述带畸变点P_u和真实点P_d的关系:

$ {\mathit{\boldsymbol{v}}_{\rm{u}}} - \mathit{\boldsymbol{e}} = \frac{{{\mathit{\boldsymbol{v}}_{\rm{d}}} - \mathit{\boldsymbol{e}}}}{{1 + {\lambda _1}r_{\rm{d}}^2 + {\lambda _2}r_{\rm{d}}^4 + \cdots }} $

(5)

2 标定方法 2.1 畸变中心估算

对单个像点, 由式(5)可知:

$ {\mathit{\boldsymbol{v}}_{\rm{d}}} = \mathit{\boldsymbol{e}} + k\left( {{\mathit{\boldsymbol{v}}_{\rm{u}}} - \mathit{\boldsymbol{e}}} \right) $

(6)

其中, k=1+λ₁r_d²+λ₂r_d⁴+…称为畸变因子。

在式(6)两边左乘向量e的叉乘反对称矩阵e_×, 由e_×·e=e×e=0可知:

$ {\mathit{\boldsymbol{e}}_ \times }{\mathit{\boldsymbol{v}}_{\rm{d}}} = k{\mathit{\boldsymbol{e}}_ \times }{\mathit{\boldsymbol{v}}_{\rm{u}}} $

(7)

其中:

$ {\mathit{\boldsymbol{e}}_ \times } = \left( {\begin{array}{*{20}{c}} 0&{ - 1}&{d{v_0}}\\ 1&0&{ - d{u_0}}\\ { - d{v_0}}&{d{u_0}}&0 \end{array}} \right) $

(8)

将式(7)两边左乘v_d^T, 由反对称矩阵的性质可知, 等式左边v_d^Te_×v_d为0, 结合式(3)可得:

$ \mathit{\boldsymbol{v}}_{\rm{d}}^{\rm{T}}{\mathit{\boldsymbol{e}}_ \times }\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{V}}_{\rm{w}}} = 0 $

(9)

令:

$ {\mathit{\boldsymbol{G}}_{\rm{r}}} = {\mathit{\boldsymbol{e}}_ \times }\mathit{\boldsymbol{P}} $

(10)

式(10)可转换为:

$ \mathit{\boldsymbol{v}}_{\rm{d}}^{\rm{T}}{\mathit{\boldsymbol{G}}_{\rm{r}}}{\mathit{\boldsymbol{V}}_{\rm{w}}} = 0 $

(11)

其中, G_r为3×4矩阵, 因此仅需最少12组已知点的坐标v_d和V_w即可求解G_r。

对式(10)两边左乘e^T得:

$ {\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{G}}_{\rm{r}}} = {\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{e}}_ \times }\mathit{\boldsymbol{P}} = {\bf{0}} $

(12)

将求取的G_r代入式(12)即可求得畸变中心的坐标e。

2.2 投影矩阵求解

由于e_×是奇异矩阵, 因此通过式(10)分解得到的投影矩阵P不唯一。为方便求解, 将像平面坐标系原点移至畸变中心$\widehat O$。此时, 畸变中心坐标为$\mathit{\boldsymbol{\widehat e}}$=(0 0 1)^T, 相应的像点坐标变为:

$ {{\mathit{\boldsymbol{\hat v}}}_{\rm{d}}} = \left( {\begin{array}{*{20}{c}} {{{\hat x}_{\rm{d}}}}\\ {{{\hat y}_{\rm{d}}}}\\ 1 \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{x_{\rm{d}}} - d{u_0}}\\ {{y_{\rm{d}}} - d{v_0}}\\ 1 \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 1&0&{ - d{u_0}}\\ 0&1&{ - d{v_0}}\\ 0&0&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{x_{\rm{d}}}}\\ {{y_{\rm{d}}}}\\ 1 \end{array}} \right) $

(13)

令转换矩阵为:

$ \mathit{\boldsymbol{T}} = \left( {\begin{array}{*{20}{c}} 1&0&{ - d{u_0}}\\ 0&1&{ - d{v_0}}\\ 0&0&1 \end{array}} \right) $

(14)

式(13)简化为:

$ {{\mathit{\boldsymbol{\hat v}}}_{\rm{d}}} = \mathit{\boldsymbol{T}}{\mathit{\boldsymbol{v}}_{\rm{d}}} $

(15)

结合式(11), 得到矩阵$\widehat {{\mathit{\boldsymbol{G}}_{\rm{r}}}}$满足:

$ \mathit{\boldsymbol{v}}_{\rm{d}}^{\rm{T}}{\mathit{\boldsymbol{T}}^{\rm{T}}}{{\mathit{\boldsymbol{\hat G}}}_{\rm{r}}}{\mathit{\boldsymbol{V}}_{\rm{w}}} = 0 $

(16)

推导出:

$ {{\mathit{\boldsymbol{\hat G}}}_{\rm{r}}} = {\mathit{\boldsymbol{T}}^{\rm{T}}}{\mathit{\boldsymbol{G}}_{\rm{r}}} = \left( {\begin{array}{*{20}{c}} 1&0&0\\ 0&1&0\\ { - d{u_0}}&{ - d{v_0}}&1 \end{array}} \right){\mathit{\boldsymbol{G}}_{\rm{r}}} $

(17)

又由于$\mathit{\boldsymbol{\widehat e}} = {(0\;\;0\;\;1)^{\rm{T}}}$, 因此:

$ {{\mathit{\boldsymbol{\hat G}}}_{\rm{r}}} = {{\mathit{\boldsymbol{\hat e}}}_ \times }\mathit{\boldsymbol{\hat P}} = \left( {\begin{array}{*{20}{c}} 0&{ - 1}&0\\ 1&0&0\\ 0&0&0 \end{array}} \right)\mathit{\boldsymbol{\hat P}} $

(18)

令$\mathit{\boldsymbol{\widehat g}}_i^{\rm{T}}$为$\widehat {{\mathit{\boldsymbol{G}}_{\rm{r}}}}$的第i行, $\mathit{\boldsymbol{\widehat p}}_i^{\rm{T}}$为$\mathit{\boldsymbol{\widehat P}}$的第j行, 则有:

$ \mathit{\boldsymbol{g}}_1^{\rm{T}} = - \mathit{\boldsymbol{p}}_2^{\rm{T}},\mathit{\boldsymbol{g}}_2^{\rm{T}} = \mathit{\boldsymbol{\hat p}}_1^{\rm{T}} $

(19)

由于e_×第3行为0, 因此$\mathit{\boldsymbol{\widehat P}}_3^{\rm{T}}$无法由$\mathit{\boldsymbol{\widehat g}}_i^{\rm{T}}$表示。在坐标系$\widehat O\widehat x\widehat y$下, 由式(3)和式(5)得到:

$ \frac{\rho }{{\hat k}}\left( {\begin{array}{*{20}{c}} {{{\hat x}_{\rm{d}}}}\\ {{{\hat y}_{\rm{d}}}}\\ {\hat k} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat p}}_1^{\rm{T}}}\\ {\mathit{\boldsymbol{\hat p}}_2^{\rm{T}}}\\ {\mathit{\boldsymbol{\hat p}}_3^{\rm{T}}} \end{array}} \right){\mathit{\boldsymbol{V}}_{\rm{w}}} $

(20)

该坐标系下新的畸变因子为:

$ \hat k = 1 + {\lambda _1}\hat r_{\rm{d}}^2 + {\lambda _2}\hat r_{\rm{d}}^4 + \cdots $

(21)

将式(21)代入式(20)中, 简化后可得:

$ \begin{array}{l} \left( {\begin{array}{*{20}{l}} {{{\hat x}_{\rm{d}}}\left( {\mathit{\boldsymbol{\hat p}}_3^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}} \right) - \left( {\mathit{\boldsymbol{\hat p}}_1^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}} \right) \cdot \left( {{\lambda _1}\hat r_{\rm{d}}^2 + {\lambda _2}\hat r_{\rm{d}}^4 + \cdots } \right)}\\ {{{\hat y}_{\rm{d}}}\left( {\mathit{\boldsymbol{\hat p}}_3^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}} \right) - \left( {\mathit{\boldsymbol{\hat p}}_2^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}} \right) \cdot \left( {{\lambda _1}\hat r_{\rm{d}}^2 + {\lambda _2}\hat r_{\rm{d}}^4 + \cdots } \right)} \end{array}} \right) = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat p}}_1^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}}\\ {\mathit{\boldsymbol{\hat p}}_2^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}} \end{array}} \right) \end{array} $

(22)

将式(22)代入式(19)转换成矩阵乘法的形式:

$ \left( {\begin{array}{*{20}{c}} {{{\hat x}_{\rm{d}}}{\mathit{\boldsymbol{V}}_{\rm{w}}}}&{\left( { - \mathit{\boldsymbol{\hat g}}_2^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}} \right)\left( {\begin{array}{*{20}{c}} {\hat r_{\rm{d}}^2}&{\hat r_{\rm{d}}^4}&1 \end{array}} \right)}\\ {{{\hat y}_{\rm{d}}}{\mathit{\boldsymbol{V}}_{\rm{w}}}}&{\left( {\mathit{\boldsymbol{\hat g}}_1^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}} \right)\left( {\begin{array}{*{20}{c}} {\hat r_{\rm{d}}^2}&{\hat r_{\rm{d}}^4}&1 \end{array}} \right)} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat p}}}_3}}\\ {{\lambda _1}}\\ {{\lambda _2}}\\ \vdots \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat g}}_2^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}}\\ { - \mathit{\boldsymbol{\hat g}}_1^{\rm{T}}{\mathit{\boldsymbol{V}}_{\rm{w}}}} \end{array}} \right) $

(23)

将$\mathit{\boldsymbol{\widehat P}}_3^{\rm{T}}$作为未知参数进行分离, 通过若干点的坐标v_d和V_w, 即可一次性线性求解$\mathit{\boldsymbol{\widehat P}}_3^{\rm{T}}$和λ_i, 点数N满足2N≥n+3, 其中n为除法模型中系数λ_i的数量。结合式(3)和式(15)得到:

$ \mathit{\boldsymbol{P}} = {\mathit{\boldsymbol{T}}^{ - 1}}\mathit{\boldsymbol{\hat P}} = \left( {\begin{array}{*{20}{c}} 1&0&{d{u_0}}\\ 0&1&{d{v_0}}\\ 0&0&1 \end{array}} \right)\mathit{\boldsymbol{\hat P}} $

(24)

将P代入$\mathit{\boldsymbol{\widehat P}}$可求得像平面坐标系Oxy下的投影矩阵P。

2.3 内外参数分离

设旋转矩阵R=(r₁ r₂ r₃), 由式(2)可得:

$ {\mathit{\boldsymbol{r}}_i} = {\mathit{\boldsymbol{A}}^{ - 1}}{\mathit{\boldsymbol{p}}_i} $

(25)

由于R为正交矩阵, 因此有:

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{p}}_i^{\rm{T}}{\mathit{\boldsymbol{A}}^{ - {\rm{T}}}}{\mathit{\boldsymbol{A}}^{ - 1}}{\mathit{\boldsymbol{p}}_j} = 0}\\ {\mathit{\boldsymbol{p}}_i^{\rm{T}}{\mathit{\boldsymbol{A}}^{ - {\rm{T}}}}{\mathit{\boldsymbol{A}}^{ - 1}}{\mathit{\boldsymbol{p}}_i} = \mathit{\boldsymbol{p}}_j^{\rm{T}}{\mathit{\boldsymbol{A}}^{ - {\rm{T}}}}{\mathit{\boldsymbol{A}}^{ - 1}}{\mathit{\boldsymbol{p}}_j}} \end{array}} \right. $

(26)

其中, 1≤i, j≤3且i≠j, 将i和j代入式(26)可以列出5个约束方程构成方程组。

内参数矩阵A可表示为:

$ \mathit{\boldsymbol{A}} = \left( {\begin{array}{*{20}{c}} {{f_x}}&s&{{x_0}}\\ 0&{{f_y}}&{{y_0}}\\ 0&0&1 \end{array}} \right) $

(27)

式(27)仅包含5个内参数:f_x和f_y为x轴和y轴方向上的焦距, s为描述x轴和y轴倾斜角的参数, (x₀ y₀)为像主点的坐标。由于未知内参数数量等于式(26)的方程数量, 因此可以直接求解B=A^TA^－1, 进而求解5个内参数。将A代入式(2)即可求解外参数矩阵R和平移向量t。

3 实验结果与分析

本文程序运行环境为:操作系统Windows10, 内存6 GB, CPU Intel i5-4200M, Matlab2016a。使用7.5 mm鱼眼镜头APS-C画幅摄像机采集图像, 靶标为两正交的棋盘格平面组成。获取的靶标图像如图 2所示。

通过角点提取方法可以得到每个角点(棋盘格顶点)在像平面的亚像素坐标v_d。建立三维坐标系如图 3所示, 根据棋盘格实际大小, 容易得到各角点(比如点P)在本体坐标系的三维坐标V_w。

将v_d和V_w作为式(11)和式(23)的输入, 即可完成标定。表 1展示了本文方法和文献[9]方法的标定结果对比。2种方法都采用立体靶标, 只需1张图像即可完成标定, 区别在于本文方法采用除法畸变模型结合线性求解方法, 文献[9]方法采用多项式畸变模型, 通过LM算法进行非线性优化求解。由表 1可见, 本文方法和文献[9]方法得到的焦距相差2.0像素~3.0像素, 主点坐标相差约9.0像素(由于文献[9]仅将张正友方法推广到三维靶标, 仍未分离畸变中心, 因此对主点精度影响较大), 从而验证本文方法的有效性。

本文通过比较5张相片的重投影误差来评估标定精度, 结果如表 2所示。另外, 将每张相片所有角点的x轴和y轴方向的重投影误差画在坐标系上, 如图 4、图 5所示。两种方法的重投影误差处于一个数量级。但是本文方法的重投影误差分布更集中, 略优于文献[9]方法。

由于本文提出的线性求解方法无需迭代, 因此运行效率比文献[9]的非线性优化方法显著提升。图 6显示了两种方法的标定时间比较, 其中文献[9]方法的平均标定时间为281.252 4 ms, 本文方法的平均标定时间为25.355 2 ms, 可见, 本文方法的标定速度是文献[9]方法的10倍左右。

4 结束语

本文针对鱼眼镜头的特性, 引入除法模型, 利用畸变中心的解耦投影矩阵和畸变系数, 提出鲁棒的鱼眼镜头线性标定方法。实验结果表明, 该方法避免了传统非线性标定方法对初值敏感的问题, 且求解过程为线性运算, 极大提高了运行效率。在实际标定过程中, 仅需已知靶标特征的空间分布, 通过1次成像即可完成高精度的标定, 适用于需要快速精确标定的场景。但由于本文线性标定方法是针对较大畸变的镜头, 因此下一步将对畸变较小的镜头进行标定研究, 使得该方法更具普适性。