Calibration method of quad-camera measurement system based on tensor decompose

基于机器视觉的三维测量具有测量速度快、测量柔性好、测量精度高以及全场测量的特点，被广泛应用在光电探测、运动捕捉与分析等领域^[1-3]。对于运动捕捉等测量应用，由于单个摄像机的视场有限，测量系统通常由多个摄像机组成，多摄像机的标定与单个摄像机的标定在标定实现上存在的主要差别是多摄像机标定不仅需要标定每个摄像机的参数还需要确定各摄像机之间的变换关系。

根据多摄像机的安装组合方式，多摄像机可以配置为汇聚式、发散式和组网式^[4]。针对不同的配置方式所需的标定方法也不尽相同，对于大范围视觉测量任务的摄像机配置方式往往结合了几种基本配置，如基于视觉的飞行器测试床，其配置方式为汇聚式与组网式相结合的组合形式，针对这种组合式的多摄像机标定，已提出了多种方法，主要分为两类：基于平面靶标与非线性优化的标定方法^[5-6]与基于一维靶标与非线性优化的标定方法^[7-10]。基于平面靶标与非线性优化标定方法的思路为利用张正友平面标定算法分别确定两相邻摄像机之间的变换关系，然后任选一个摄像机坐标系作为世界坐标系，将所有摄像机都统一到该坐标系之下，最后利用非线性优化算法，以摄像机之间的变换关系固定为约束条件，进行非线性优化求解，这种方法的优点是测量算法思路简单，但其标定结果具有累积误差，且存在视场遮挡问题，难以得到高精度的标定结果；基于一维靶标与非线性优化的标定方法又可分为两类，一类为非自由式，其基本思路为将一维靶标一端固定，在空间中摆动另一端，利用一维靶标可以同时被多个摄像机观测到的特点实现对多摄像机的标定，并结合非线性优化算法提高标定精度，该类算法解决了多摄像机标定的视场遮挡问题，但需要对标定靶标的运动进行约束，如绕固定点或在平面内运动等，标定精度不高，且也存在标定误差累计的问题；另一类为自由式，其基本思路为将一维标定靶标在空间中自由移动，利用各摄像机采集的图像数据与一维标定杆中标记点之间的距离信息，通过分层重建的方式实现对多摄像机的标定，这类方法不存在视场遮挡的情况，且无需对标定靶标的运动进行约束，操作方便，标定精度高，但仍存在标定误差累积的问题。对于标定结果的非线性优化目前主流的多摄像机标定的非线性优化算法主要有ICP(Iterative Closet Point)最近点迭代^[11]、光束法平差(BA)^[12]、IO(Orthogonal Iteration)正交迭代^[13]，不同方法对应着不同目标函数。

文中针对运动捕捉系统中现有多摄像机标定方法存在累积误差且标定精度不高的问题，提出了一种基于多焦点张量的四摄像机之间结构参数全局标定方法，该方法利用具有两个标记点的一维靶标作为标定物，通过摄像机采集一维标定靶标在测量空间内自由运动时的图像，然后利用改进Zernike矩算法对标定点图中坐标进行亚像素定位^[14]，并由对应的图像坐标计算多摄像机之间的多焦点张量，实现多摄像机的射影标定，根据一维标定杆的距离信息实现多摄像机的欧氏标定。最后，以标定点图像与其反向投影点之间的Sampson距离为代价，建立目标函数，利用光束法平差算法实现对多摄像机标定参数的优化，实现多摄像机之间结构参数的高精度标定。

对于多摄像机测量系统，其标定任务包括每个摄像机的内部参数和各摄像机之间的旋转平移矩阵的求解。

假设2D图像坐标点定义为${\boldsymbol{m'}} = {\left( {u,v} \right)^{\rm{T}}}$（单位：pixel），${\boldsymbol{M'}} = {\left( {x,y,{\textit{z}}} \right)^{\rm{T}}}$定义为空间三维点（单位：mm），${\boldsymbol{m}} = {\left( {u,v,1} \right)^{\rm{T}}}$与${\boldsymbol{M}} = {\left( {x,y,{\textit{z}},1} \right)^{\rm{T}}}$定义为对应2D与3D点的齐次坐标，则根据透视投影成像模型，空间点的成像过程可以描述为：

式中：$s$为尺度因子；${\boldsymbol{K}}$为摄像机内部参数；$\left[ {{\boldsymbol{R}},{\boldsymbol{t}}} \right]$表示摄像机的外部参数；${\boldsymbol{R}}$为一个3×3的旋转矩阵；${\boldsymbol{t}}$为一个三维平移向量；${f_u}$与${f_v}$分别表示单位为pixel的图像坐标系$u$轴与$v$轴方向的等效焦距；${u_0}$与${v_0}$表示摄像机的主点；$\gamma $表示图像坐标系$u$轴与$v$轴的不垂直因子，非常小，可以忽略；${\boldsymbol{P}}$表示摄像机矩阵。

如图1所示，A、C为一维标定杆上的两个标记点，其中$\left\| {AC} \right\| = d$，$d$为标定杆长度，$ \left( {{a_{ij}},{c_{ij}}} \right) $，$ i = 1,2, \cdots , m $，$j = 1,2, \cdots ,n$分别为标定杆上A、C两点在第$j$次运动后在4个摄像机上的成像点，其中$m$表示测量系统中摄像机个数，$n$表示运动的次数。

Figure 1.  Geometry diagram of multi-camera calibration

若将世界坐标系设为第一个摄像机坐标系，则系统中所有摄像机相对于世界坐标系的摄像机矩阵${{\boldsymbol{P}}_i}$可以表示为：

式中：${{M}} \geqslant 2$表示摄像机的数目；${{\boldsymbol{R}}_{01}} = {\boldsymbol{I}}$，${{\boldsymbol{R}}_{\left( {m - 1} \right)m}}$与${{\boldsymbol{t}}_{n\left( {n - 1} \right)}}$分别表示两相邻摄像机之间的旋转矩阵与平移向量；${\boldsymbol{I}}$表示单位矩阵；${{\boldsymbol{K}}_i}$表示系统中每一摄像机的内部参数矩阵。

多摄像机成像模型可以表示为：

式中：${\boldsymbol{x}} = \left( {u,v} \right)$表示标定点的图像坐标；${\boldsymbol{X}} = \left( {x,y,z} \right)$表示标定点的三维坐标。

首先，根据摄像机之间的位置与摄像机的个数将摄像机分组；然后，求解每两两相邻摄像机之间的基本矩阵；最后，利用矩阵分解算法求解摄像机的外部参数。

如图1所示的4个摄像机组成的测量系统，可将摄像机分为(1,2)、(2,3)、(3,4)、(4,1)、(1,3)与(2,4)6个组，对于每组的基本矩阵可表示为${{\boldsymbol{F}}_{C_6^2}}$，利用标定杆运动后采集的图像坐标${m_{ij}}$，通过归一化8点法，求得基本矩阵${\boldsymbol{F}}$。又由极线几何可知：

式中：${\boldsymbol{E}}$为两摄像机之间本质矩阵；${{\boldsymbol{K}}_i}$为摄像机的内部参数矩阵；$i$为摄像机编号。

因为摄像机内部参数在测量过程中始终保持不变，因此可以利用二维平面靶标标定得到^[5]。

求得摄像机内部参数${{\boldsymbol{K}}_i}$后，两摄像机之间的本质矩阵可以表示为：

式中：${\boldsymbol{S}}$为由归一化平移向量${\boldsymbol{t'}}$定义的反对称矩阵。又根据本质矩阵的性质可知：

平移向量$ {\boldsymbol{t'}} $可以通过求解公式（6）得到：

因此，${\boldsymbol{t'}} = {\boldsymbol{U}}{\left( {0,0,1} \right)^{\text{T}}}$，由公式(4)可知，本质矩阵仅刻画摄像机之间的射影变换关系，与欧氏空间的平移向量${\boldsymbol{t}}$相差一个尺度因子。

根据公式(5)所示本质矩阵的形式，由参考文献[15]可知：若假设矩阵E可被分解为一个正交矩阵R与一个反对称矩阵S的乘积形式，则矩阵E可被奇异值分解（SVD）为UDV^T的形式，其中D=diag(1,1,0)，且在相差一个尺度因子的条件下有：

式中：${\boldsymbol{S}} \;\;\approx\;\; {\boldsymbol{UZ}}{{\boldsymbol{U}}^{\rm{T}}}$；${{\boldsymbol{R}}_{\text{1}}} \;\;\approx\;\; {\boldsymbol{UQ}}{{\boldsymbol{V}}^{\rm{T} }}$或${{\boldsymbol{R}}_{\text{2}}} \;\;\approx\;\; {\boldsymbol{U}}{{\boldsymbol{Q}}^{\rm{T}}}{{\boldsymbol{V}}^{\rm{T}}}$，${\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} 0&{ - 1}&0 \\ 1&0&0 \\ 0&0&0 \end{array}} \right]$，$ {\boldsymbol{Z}} = \left[ {\begin{array}{*{20}{c}} 0&1&0 \\ { - 1}&0&0 \\ 0&0&1 \end{array}} \right] $。

由上述引理可知，存在四组可能的解：

利用上述求解得到的相差一个尺度因子下的摄像机参数${\boldsymbol{R}}$、${\boldsymbol{t'}}$，结合摄像机的内部参数${\boldsymbol{K}}$求解一维靶标上标记点的坐标值，由于靶标点位于两个摄像机的前部，因此其坐标的$z$分量大于0，据此可以判断分析出四组可能解中的正确解。至此，求得了相差一个比例因子条件下的两摄像机之间的旋转平移向量。

而比例因子$k$又可以利用一维标定杆的距离信息进行求解，具体为：

式中：$d$为一维标定杆在欧氏空间的实际长度；$d'$为在$\left\| {{\boldsymbol{t'}}} \right\| = 1$的情况下的测量长度。

若以第一个摄像机坐标系作为世界坐标系，则多摄像机的标定参数可分别表示为：

上一节多摄像机标定过程具有一维靶标同时被四个摄像机观测的特点，而标定算法确将其分组，又由多视图几何可知，如图1所示的4个图像平面之间的多重线性关系可类比于双目立体视觉的基本矩阵描述两视图之间的线性变换关系，用四焦点张量进行唯一描述^[15]，且相较于基本矩阵与三焦点张量的方法^[16]，基于四焦点张量的方法因具有更多的冗余信息而具有更高的精度与更好的算法稳定性。因此，此节考虑将四焦点张量引入四摄像机测量系统的结构参数标定中。

对于一维靶标在4个摄像机同时采集的四幅图像中的同名靶标点${\boldsymbol{x}}$，${\boldsymbol{x'}}$，${\boldsymbol{x''}}$，${\boldsymbol{x'''}}$的图像齐次坐标与多摄像机的四焦点张量${Q^{pqrs}}$之间的关系可以表示为如下形式：

式中：$\varepsilon $为三维空间中的交错符号；$i$、$j$、$k$、$l$为标记点图像坐标的坐标分量标号；$w$、$x$、$y$、$z$的不同取值可表示4个齐次坐标点之间不同的多重线性关系，所有上下标的取值均在$\left( {1,2,3} \right)$中。

同时四焦点张量$ Q^{pqrs} $被定义为：

且向量${{\boldsymbol{a}}^p}$，${{\boldsymbol{b}}^q}$，${{\boldsymbol{c}}^r}$，${{\boldsymbol{d}}^s}$表示摄像机矩阵${{\boldsymbol{P}}_1}$，${{\boldsymbol{P}}_2}$，${{\boldsymbol{P}}_3}$，${{\boldsymbol{P}}_4}$的行向量，$p$，$q$，$r$，$s$为逆变指标。

展开公式(12)可知，其仅有16个线性独立的方程。利用4个摄像机采集的标定杆图像之间的点对应$\left\{ {{a_{0 j}} \leftrightarrow {a_{1 j}} \leftrightarrow \cdots {a_{mj}}} \right\}$，$\left\{ {{c_{0 j}} \leftrightarrow {c_{1 j}} \leftrightarrow \cdots {c_{mj}}} \right\}$，其中$m$表示摄像机的个数，$j = 1,2, \cdots ,n$，$n$表示标定杆运动的次数，解线性方程组(12)即可解得四焦点张量${Q^{pqrs}}$，其表示4个摄像机矩阵的本质关系。

通常来说，一组对应点可以产生$4 \times 4 = 16$个线性独立的方程，然而经分析得出，对应于两个不关联的点对应的方程组有相关性。因此一组对应点并不能求解方程式(12)。

又因为4个摄像机系统具有$11 \times 4 - 15 + 3 n$个自由度，而$m$幅图像上的$mn$个点只能提供$2 mn$个方程，所以为了估计$ {Q^{pqrs}} $至少需要满足：

式中：$round\left( x \right)$表示对$x$向前取整。

因此至少需要6个独立的三维点，然而对于距离已知的一维标定杆上的两个靶标点坐标并不是相互独立的。

又因为一个三维点具有3个自由度，而只有两个标定点的一维标定杆仅具有5个自由度，因此，若设标定杆位姿变化次数为$M$，则有：

因此，经分析标定杆在空间中的运动次数至少为三次。

最后，同样利用一维标定杆之间的距离信息，求得射影意义下的摄像机矩阵与欧氏意义下的摄像机矩阵之间的比例因子$k$，即可获得多摄像机标定参数。

若以第一个摄像机坐标系作为世界坐标系，为了从四焦点张量恢复摄像机矩阵，根据参考文献[15]提出的张量简化方法可知，摄像机矩阵可以表示为一个对角矩阵与一个由极点形成的列向量组成，则4个摄像机的标定参数可分别表示为：

式中：${b_i}$，${c_i}$，${d_i}$与${b'_i}$，${c'_i}$，${d'_i}$为利用简化后的四焦点张量，通过公式(13)解得$i = 1,2,3$。简化后的张量${\hat Q^{abcd}}$与未简化的张量${Q^{pqrs}}$之间的关系可表示为：

式中：${\boldsymbol{T}}$, ${\boldsymbol{T'}}$, ${\boldsymbol{T''}}$, ${\boldsymbol{T'''}}$为图像坐标简化矩阵的逆矩阵。

由公式（12）、（13）、（16）可知，${P_i}$的值仅与图像坐标${\boldsymbol{x}}$, ${\boldsymbol{x'}}$, ${\boldsymbol{x''}}$, ${\boldsymbol{x'''}}$相关，其中$ i=2,3,4 $，因此采用本小节方法标定得到的摄像机矩阵相互之间不会进行误差的传递与累积，另一方面，四焦点张量描述的是空间同一点在4个摄像机成像平面上点坐标之间的多重线性关系，对于空间的定位约束相较于利用两两相邻约束标定得到的摄像机参数来定位空间点坐标更加准确。因此，该方法在理论上可以获得相较于基于基本矩阵法的多摄像机标定法具有更高的标定精度。

如果图像坐标含有测量噪声，则上一节求解出的摄像机矩阵将不满足公式(18)：

式中：${{\boldsymbol{x}}_{ij}}$为第$j$个点在第$i$个摄像机中的图像坐标；${{\boldsymbol{X}}_j}$表示第$j$个三维点。在本节中取$i = 1,2,3,4$，$j = 1,2,3$。

为了对含有噪声条件下的摄像机矩阵与空间三维点进行精确的估计，建立如下目标函数：

因为测量噪声相对于测量值很小，因此$d\left( {x,y} \right)$取齐次坐标点$x$与$y$之间的Sampson距离，减小计算量。

在上述目标函数最小化的条件下利用Levenberg-Marquardt优化方法求得摄像机矩阵与空间三维点后，计算一维标定杆标记点之间的距离并将其与标准值相比较，有如下目标函数：

式中：${d_2}\left( {x,y} \right)$表示在三维空间中两点之间的几何距离；${d_1}$表示一维标定杆标定点之间的实际距离。利用Levenberg-Marquardt优化方法求解最小化目标函数公式(19)、（20），得到精确的摄像机矩阵。

为了验证文中提出算法，设计了数值仿真分析与实际标定试验。

根据对两种多摄像机标定算法的描述，多摄像机标定误差主要来源于图像坐标提取误差、摄像机标定内部参数误差。该节将通过仿真实验研究定位噪声与原理误差对多摄机标定精度的影响。

实验中，摄像机的分辨率为2048×2048 pixel，4个摄像的内部参数$ \left( {\alpha,\;\beta,\; \gamma,\;{u_0},{v_0}} \right) $均设置为：(1100, 1100,0, 1024, 1024)，一阶畸变系数${k_{i1}} = 0.01$，二阶畸变系数${k_{i2}} = 0.001$。如图2所示，将4个摄像机设置在边长为6000的正方形$ {O_1}{O_2}{O_3}{O_4} $的4个顶点处，摄像机光轴与正方形的对角线重合；一维靶标杆AC之间的距离$AC = {\text{5}}00\;{\rm{mm}}$。

Figure 2.  Schematic of multi-camera calibration simulation structure

仿真实验中利用残差与重建误差来验证文中所提算法的优越性，其中残差公式定义为：

式中：${d_x} = {x_i} - {\hat x_i}$；${d_{x'}} = {x'_i} - {\hat x'_i}$；${d_{x''}} = {x''_i} - {\hat x''_i}$；${d_{x'''}} = {x'''_i} - {\hat x'''_i}$，${x_i}$，$ {x'_i} $，${x''_i}$，${x'''_i}$表示特征点在4个摄像机中的同名成像点；${\hat x_i}$，$ {\hat x'_i} $，${\hat x''_i}$，${\hat x'''_i}$表示根据求解得到的摄像机矩阵将空间三维特征点反向投影到4个摄像机中的成像点。

重建误差公式定义为：

式中：$\left( {{x_s},{y_s},{z_s}} \right)$表示仿真生产的空间三维点坐标；$\left( {{x_r},{y_r},{z_r}} \right)$表示三维重建生产的空间点坐标。

每次标定实验时，令线段AC在正方形中心附近做20次任意的刚体运动，为了保证算法的统计特性，在每一固定成像误差下重复进行100次标定实验，且保证标定点都在4个摄像机的成像平面内。给图像坐标添噪声水平在0.0~1.1 pixel的成像误差。统计每一误差等级下标定杆长度解算误差模的均值。其算法精度对比实验结果如图3所示。

Figure 3.  Comparison curve of re-projection error of multi-camera calibration changing with noise

根据图3可知：

（1) 五种摄像机的标定误差等级都随图像坐标定位结果的信噪比的减小而增大，但平面法无法同时被4个摄像机同时观测到，因此在某些视角下无法实现四摄像机测量系统的标定；

（2) 相比与其他摄像机标定算法，当图像坐标定位精度为1 pixel时，文中提出的多摄像机标定算法的反向投影误差保持在3 mm以内；

（3) 加权基本矩阵法、三焦点张量法与文中提出的算法之间的标定大致相当，但由于文中所提出的四摄像机标定算法直接分解四焦点张量即可实现，因此其效率更高。

6426个标定点被生成且均匀分布在一个半径为2 000 mm的球体内，重建不确定性大小为1.0 pixel，均值为0。每一个点的重建误差$er{r_{rec}}$被统计，其归一化的误差统计结果如图4所示。其统计参数如表1所示。

根据图4可知，文中提出的算法误差最为集中，主要集中在2.2~2.4 mm之间，且呈现正态分布的形式，这说明文中的方法对于每一误差项的影响是均匀的。其他算法都有一个拖尾效应，在其误差分布上可以明显地看出耦合的有指数分布，这说明不同的位置，其误差因素的影响是不一样的，因此文中提出的方法更加得合理。

Figure 4.  Re-projection error histogram of spatial point location by quad-camera

表1给出了四种一维标定算法误差的均值、标准差、最大值、以及均方根，对比结果可以得出，文中提出的算法的误差均值与均方根相比于加权基本矩阵法稍大以外，其他统计指标都更小，这说明了文中提出的算法即可达到加权基本矩阵法的精度，又具有更集中的误差分布，及更快的标定效率。

为了验证算法的可行性与精度，搭建了如图5(a)所示多摄像机测量系统。该系统由4台摄像机、4个LED红外光源，4个磁盘存储阵列、一个T型一维标定杆，一个相机同步控制器和1台计算机组成；其中摄像机的型号为IO Industries Flare 4 M140，CCD分辨率为2 048×2 048 pixel，像元大小为0.0055 mm；镜头为焦距为25 mm的CHIOPT HC1605 A；红外光源波长为850 nm；一维标定杆上标定球直径为d=12.7 ，标定杆长度500 mm，其表面的红外反光材料为3 M7610；系统标定时将4个摄像机如图5(b)所示的布局方式固定在衍架上。

Figure 5.  Multi-camera measurement system

如图6所示，在测量空间中挥舞标定杆80次，使得标定点均匀地分布在测量空间内，并采集其图像。标定点图像坐标采用参考文献[14]提出的方法进行定位，其定位误差为1 pixel。

Figure 6.  Quad-camera external parameter calibration experiment based on 1D calibration rod

根据上述摄像机的内部参数以及采集到的1D标定杆图像坐标，然后利用文中提出的四摄像机标定方法进行标定，其结果（每一次挥杆两靶标点之间距离误差统计）如图7、图8所示。其中摄像机的内部参数由二维标定算法^[5]分别进行精确标定得到。系统中每一台摄像机的内部参数如表2所示。摄像机的外部参数如表3所示。

Figure 7.  Statistical error chart of spatial point locating of proposed multi-camera calibration method

Figure 8.  Statistical error chart of spatial point locating of proposed multi-camera calibration method

从图7~8可以看出：

(1) 文中提出的四摄像机标定算法对空间标准杆长度测量误差小于4 mm（3σ），满足室内运动目标的三维定位需求；

(2) 实际测量值大于仿真分析所得值的原因是因为现场光照条件复杂，标定靶标点的图像灰度分布不满足高斯分布的假设，且对摄像机的内部参数标定也存在误差，在利用这些数据进行标定时必然会将误差传递到定位结果上。

为了实现对室内运动目标位姿的精确测量，建立了一套多摄像机测量系统，并通过分析现有基于一维靶标的多摄像机标定方法的不足，提出了一种基于多焦点张量的多摄像机标定方法，并通过仿真分析与实际测试实验验证了基于四焦点张量分解方法显著提高了标定操作的效率，且标定精度相较于基于基本矩阵的方法提高了近40%，在4000 mm×4000 mm×2000 mm的测量范围内测量精度达到4 mm，解决了四摄像机测量系统现场标定操作繁琐，精度较差的问题。另外，值得指出的是文中所提出的采用四焦点张量分解的多摄像机标定方法需要对张量进行分解操作，计算比基本矩阵对方法更为复杂，因此后续将在减小四焦点张量分解计算量的同时，考虑将文中提出的方法推广到由更多的摄像机组成的超大空间摄影测量网外部旋转与平移结构参数的现场标定中。