CN101365140A - 基于矩阵分解的摄像机阵列标定方法 - Google Patents

Abstract

本发明公开了一种基于矩阵分解的摄像机阵列标定方法。包括以下步骤：1)对摄像机编号并对标定板拍照，同时将所得的图像也按照拍照顺序及其对应的摄像机编号；2)利用标定板上特征点的实际物理坐标数据和拍照获得的图像数据求取摄像机的单应性矩阵；3)利用求得的单应性矩阵及其对应的摄像机和标定板的编号构造测量矩阵，并调整单应性矩阵的比例因子；4)根据摄像机阵列的排列情况——密集或稀疏，选择对应的方法对存在数据缺失的测量矩阵

Description

基于矩阵分解的摄像机阵列标定方法

技术领域

本发明属于计算机多媒体技术领域，特别涉及立体视觉，3DTV(三维电视)，IBR(基于图像的绘制技术)中的摄像机阵列标定技术。

背景技术

摄像机标定，就是求取摄像机的成像几何模型的过程，该几何模型的参数就是摄像机参数。其中包括摄像机的内部几何光学参数(内部参数)A_i以及摄像机相对于世界坐标系的位置关系(外部参数)，包括旋转矩阵R_i和平移向量t_i。其中A_i形式如下：

$A_i=\left(\begin{array}{ccc}{\mathrm{\α}}_i& c_i& u_{0i}\\ 0& {\mathrm{\β}}_i& v_{0i}\\ 0& 0& 1\end{array}\right)$

对于单个摄像机的标定而言，目前已有很多很成熟的算法，包括：线性法，非线性优化法，两步法，双平面法等。

但随着多媒体技术的发展，单个摄像机已逐渐地不能满足人们的需求。无论是目前正在兴起的3DTV(三维电视)还是通过多幅图像的绘制来获得更高质量图像的IBR技术，都需要多个摄像机，甚至摄像机阵列同时工作。因此，如何精确，简单的获得摄像机阵列中各个摄像机的内外参数，并进而由此获取物体的空间信息，深度信息，成为摄像机阵列应用的核心内容。

目前，国外对摄像机阵列有着比较广泛的研究。MIT早期曾建立过一个3-Droom实验环境，通过从不同的角度对房间内部场景同时拍摄，来实现内部场景的三维重建。Stanford大学组建的摄像机阵列拥有约100多个摄像机，并将阵列按照排列的密集程度的不同来实现不同的应用需求和图像处理。此外，CMU大学的摄像机阵列是一个可以自动配置和移动的摄像机阵列，从而使整个系统在使用IBR技术进行场景渲染时能取得更好的效果。对类似以上几个这样大规模的摄像机阵列的应用，如果我们还使用单个摄像机逐个标定的方法，将造成：一方面，算法的复杂度很高；另一方面，系统的精度也难以达到要求。

摄像机阵列的标定，现有的方法主要有：(1)张正友算法的扩展。首先对摄像机逐一进行标定，求取摄像机的内参。然后利用前一步得出的摄像机和标定板之间的相互位置关系，通过非线性优化求出摄像机与摄像机之间的相互关系；(2)Peter Sturm，Andrei Zaharescu，Henrik Aanaes等人也先后提出了在摄像机已单独标定好的情况下，利用摄像机和平面板标定物之间的旋转角度矩阵的秩3约束，通过矩阵分解求取摄像机与摄像机之间的旋转角度和位移；(3)ToshioUeshiba提出的利用摄像机和平面标定板之间的单应性矩阵构造测量矩阵(Measurement Matrix)，而该测量矩阵又同时满足秩4约束，从而将测量矩阵分解为摄像机参数矩阵和标定板参数矩阵，来对摄像机阵列进行标定。同第一种方法相比，方法(3)使标定所求得的摄像机参数之间具有更好的一致性。这种一致性将使得后面在进行深度提取和3DTV中的虚拟视点重建具有更高的精确度。同第二种方法相比，方法(3)不再要求在求取摄像机的空间位置姿态时摄像机已经事先标定准确的限制，使得标定过程更为灵活。

但是在Toshio Ueshiba的方法中，假定摄像机之间的视野是基本重合的，或者说标定过程中对不同位置的标定板进行拍照时，这些标定板必须被所有的摄像机同时看到。然而，当我们的摄像机阵列比较庞大，或者摄像机阵列排列的比较稀疏，比较特殊(如：视频监控系统，MIT的3-D ROOM)时，我们就不可能保证所有位置的标定板都能被所有的摄像机拍摄到。即便是摄像机阵列比较密集的情况，我们也很难做到每个位置的标定板都能被所有的摄像机拍摄到，并且图像质量能达到进行标定的要求。

当摄像机拍摄不到某个标定板，或该标定板的图像质量不能达到要求时，我们就无法求出该摄像机和标定板之间的单应性矩阵，此时Toshio方法采用SVD对测量矩阵进行分解就不再可行了。此外为了对测量矩阵进行分解，需要对单应性矩阵的比例因子做相应的调整，以保证测量矩阵分解时的秩4约束。但当某一位置标定板不在某摄像机的视野中时，这个缺失的单应性矩阵有可能会影响到其它单应性矩阵比例因子的调整，这时测量矩阵就不再满足秩4约束，矩阵的分解也不再可行。

发明内容

本发明的目的是为了克服现有标定技术的不足，充分解决摄像机阵列比较庞大，或摄像机的分布比较分散而造成的摄像机视野相差很大的情况，而提供一种基于矩阵分解的摄像机阵列标定方法

本发明的基于矩阵分解的摄像机阵列标定方法，包括以下步骤：

1)将摄像机编号并对标定板拍照，同时将所得的图像也按照拍照顺序及其对应的摄像机编号；

2)利用标定板上特征点的实际物理坐标数据和拍照获得的图像数据求取摄像机的单应性矩阵

3)将求得的单应性矩阵按照其对应的摄像机和标定板的编号构造测量矩阵

4)根据摄像机阵列的密集或稀疏排列情况，调整单应性矩阵的比例因子，并对丢失数据的测量矩阵

进行矩阵分解；

5)添加度量约束，求解摄像机的内外参数；

6)对所得内外参数进行非线性优化。

本发明中，所说的利用标定板上特征点的实际物理坐标数据和拍照获得的图像数据求取摄像机的单应性矩阵的步骤如下：

a)摄像机的成像模型如式(1)

其中，X表示空间某点在世界坐标系中的坐标，(uv)为该空间点投影到摄像机图像坐标系中的坐标，A_i表示第i个摄像机的内参矩阵，R_i，t_i为摄像机坐标系在世界坐标系中的旋转和位移，

指等式两边相差一个比例因子，

标定板上特征点坐标与世界坐标系之间的转换关系如式(2)

以上两式中，P_i为摄像机参数矩阵，Q^j为标定板的位置参数矩阵；

b)对图像数据和已知的标定板上特征点的实际坐标数据做规范化处理，并按式(3)求取摄像机的单应性矩阵 ${\tilde{H}}_i^j={\mathrm{\λ}}_i^j{H}_i^j,$

其中的

表示第i个摄像机和第j个标定板图像之间的真实单应性矩阵。

本发明中，所说的将求得的单应性矩阵按照其对应的摄像机和标定板的编号构造测量矩阵的方法是，将步骤2)中得到的摄像机和标定板之间的单应性矩阵

依照其上下标按式(4)所示方式排列，得到测量矩阵

即矩阵的“行标”表示摄像机的标号，“列标”表示标定板图像的标号，每一个

都是一个(3×3)的矩阵，

$\tilde{W}=\left[\begin{array}{ccc}\widetilde{H_1^1}& \·\·\·& \widetilde{H_1^J}\\ \·& & \·\\ \·& & \·\\ \·& & \·\\ \widetilde{H_I^1}& \·\·\·& \widetilde{H_I^J}\end{array}\right]---\left(4\right)$

其中得

表示实际计算出来的单应性矩阵与真实的单应性矩阵相差一个比例因子即： ${\tilde{H}}_i^j={\mathrm{\λ}}_i^j{H}_i^j.$

本发明中，所说的根据摄像机阵列的密集或稀疏排列情况，调整单应性矩阵的比例因子，并对丢失数据的测量矩阵进行矩阵分解，其步骤按密集或稀疏排列情况分别如下：

a)当摄像机阵列排列密集时

(1)调整单应性矩阵的比例因子

调整 ${\tilde{H}}_i^j={\mathrm{\λ}}_i^j{H}_i^j$ 中的比例因子

方法如下：

对第l和第f个摄像机来说，它们分别关于第r和第p个标定板图像的相关单应性矩阵

形式如下：

$N_{\mathrm{l\; f}}^{\mathrm{r\; p}}=H_f^r{\left(H_l^r\right)}^{-1}{H}_l^p{\left(H_f^p\right)}^{-1}$

则计算得出的带有比例因子的相关单应性矩阵为：

${\tilde{N}}_{\mathrm{l\; f}}^{\mathrm{r\; p}}={\tilde{H}}_f^r{\left({\tilde{H}}_l^r\right)}^{-1}{\tilde{H}}_l^p{\left({\tilde{H}}_f^p\right)}^{-1}=\left({\mathrm{\λ}}_f^r{\mathrm{\λ}}_l^p\right)/\left({\mathrm{\λ}}_l^r{\mathrm{\λ}}_f^p\right)N_{\mathrm{l\; f}}^{\mathrm{r\; p}}---\left(5\right)$

具有一个值为

的二重特征根如式(6)

$u_{\mathrm{l\; l}}^{\mathrm{r\; p}}=\left({\mathrm{\λ}}_f^r{\mathrm{\λ}}_l^p\right)/\left({\mathrm{\λ}}_l^r{\mathrm{\λ}}_f^p\right)---\left(6\right)$

将步骤3)构造的测量矩阵

中的元素按行和列排序，使其前三行对应的摄像机能拍摄到所有的标定板，前三列对应的标定板能被所有的摄像机同时拍摄，并固定此时的测量矩阵中的第一个摄像机和第一个标定板作为计算相关单应性矩阵的参考摄像机和参考标定板，则第一个摄像机和第一个标定板与第i个摄像机的第j个标定板图像的相关单应性矩阵为： ${\tilde{N}}_{l1}^{r1}=u_{l1}^{r1}{N}_{l1}^{r1},$ 且

通过

的二重特征根计算得出，记

为

对

中的每一个元素做如下调整，

令 ${\tilde{H}}_i^j=u_i^j{\tilde{H}}_i^j=u_i^j{\mathrm{\λ}}_i^j{H}_i^j=({\mathrm{\λ}}_i^1{\mathrm{\λ}}_1^j/{\mathrm{\λ}}_1^1)P_i{Q}^j---\left(7\right)$

则调整后的

为：

$\tilde{W}=\left[\begin{array}{ccc}\widetilde{H_1^1}& \·\·\·& \widetilde{H_1^J}\\ \·& & \·\\ \·& & \·\\ \·& & \·\\ \widetilde{H_I^1}& \·\·\·& \widetilde{H_I^J}\end{array}\right]$

(2)丢失数据的测量矩阵

的分解

按以下三个步骤进行：

1)对一个秩为4，大小为3I×3J的测量矩阵给定任意初值的3I×4矩阵P₀；

2)在P_k-1已知的情况下，按式 $B{(j,:)}^T=\left(\tilde{W}{(:,j)}^{T}A\right){\left(A^{T}A\right)}^{-1}(j=1\·\·\·3J)$ 求取3J×4矩阵Q_k，并对Q_k的列向量做正交化处理：Q_k＝Q_k×N_k，使得最小，ij代表对

存在的因子求和；

3)在Q_k已知的情况下，按式 $A{(i,:)}^T=\left(\tilde{W}(i,:)B\right){\left(B^{T}B\right)}^{-1}(i=1\·\·\·3I)$ 求取3I×4矩阵P_k，使得 $\underset{\mathrm{ij}}{\mathrm{\Σ}}{|{\tilde{W}}_{\mathrm{ij}}-{\left(P_k{Q}_k^T\right)}_{\mathrm{ij}}|}^2$ 最小，ij代表对

存在的因子求和；

4)重复以上2)，3)两个步骤直到

收敛，并对P，Q逆排序后，得到测量矩阵

的分解形式如式(8)：

$\tilde{W}=(1/{\mathrm{\λ}}_1^1)\left[\begin{array}{c}{\mathrm{\λ}}_1^1{P}_1\\ \·\\ \·\\ \·\\ {\mathrm{\λ}}_I^1{P}_I\end{array}\right]\left[\begin{array}{ccc}{\mathrm{\λ}}_1^1{Q}^{1^T}& \·\·\·& {\mathrm{\λ}}_1^J{Q}^{J^T}\end{array}\right]---\left(8\right)$

当迭代不收敛或者误差

超过阈值时，更换迭代初值P₀重新迭代即可。

b)当摄像机阵列的排列稀疏时

(1)对步骤3)中构造的测量矩阵

重排序，使得

的前三行对应的摄像机能拍摄到最多数目的标定板图像，前三列对应的标定板能被最多数目的摄像机同时拍摄到，并固定此时的测量矩阵中的第一个摄像机和第一个标定板作为计算相关单应性矩阵的参考摄像机和参考标定板，则第一个摄像机和第一个标定板与第i个摄像机的第j个标定板图像的相关单应性矩阵为： ${\tilde{N}}_{l1}^{r1}=u_{l1}^{r1}{N}_{l1}^{r1},$ 且

通过

的二重特征根计算得出，记

为

对

中前三行和前三列不为空的

和

对应的那些行和列的

做如下调整，

令 ${\tilde{H}}_i^j=u_i^j{\tilde{H}}_i^j=u_i^j{\mathrm{\λ}}_i^j{H}_i^j=({\mathrm{\λ}}_i^1{\mathrm{\λ}}_1^j/{\mathrm{\λ}}_1^1)P_i{Q}^j---\left(7\right)$

则调整后的

为：

$\tilde{W}=\left[\begin{array}{ccc}\widetilde{H_1^1}& \·\·\·& \widetilde{H_1^J}\\ \·& & \·\\ \·& & \·\\ \·& & \·\\ \widetilde{H_I^1}& \·\·\·& \widetilde{H_I^J}\end{array}\right]$

(2)丢失数据的测量矩阵的分解

按以下几个步骤进行：

1)对一个秩为4，大小为3I×3J的矩阵

给定任意初值的3I×4矩阵P₀；

2)设每一个

的比例因子 ${\mathrm{\λ}}_i^j=1;$

3)利用当前的

构建测量矩阵

4)在P_k-1已知的情况下，按式 $Q_k{(j,:)}^T=\left(\tilde{W}{(:,j)}^T{P}_{k-1}\right){\left({P_{k-1}}^T{P}_{k-1}\right)}^{-1}(j=1\·\·\·3J)$ 求取3J×4矩阵Q_k，并对Q_k的列向量做正交化处理：Q_k＝Q_k×N_k，使得

最小，ij代表对

存在的因子求和；

5)在Q_k已知的情况下，按式 $P_k{(i,:)}^T=\left(\tilde{W}(i,:)Q_k\right){\left({Q_k}^T{Q}_k\right)}^{-1}(i=1\·\·\·3I)$ 求取3I×4矩阵P_k，使得

最小，ij代表对存在的因子求和；

6)调整

中前三行和前三列为空的

和

对应的那些行和列的

的比例因子如下： ${\mathrm{\λ}}_i^j=\mathrm{\Σ}{P}_i(3,:)\×Q_j{(3,:)}^T,$ 其中P_i(3，：)代表P_i的第三行，Q_j(3，：)代表Q^j的第三列；

7)重复步骤3)到6)直到

收敛，并对P，Q逆排序后，得到的分解形式如式(8)：

$\tilde{W}=(1/{\mathrm{\λ}}_1^1)\left[\begin{array}{c}{\mathrm{\λ}}_1^1{P}_1\\ \·\\ \·\\ \·\\ {\mathrm{\λ}}_I^1{P}_I\end{array}\right]\left[\begin{array}{ccc}{\mathrm{\λ}}_1^1{Q}^{1^T}& \·\·\·& {\mathrm{\λ}}_1^J{Q}^{J^T}\end{array}\right]---\left(8\right)$

当迭代不收敛或者误差

超过阈值时，更换迭代初值P₀重新迭代即可。

本发明中，所说的添加度量约束，求解摄像机的内外参数的步骤如下：

a)对P，Q去规范化，得到测量矩阵分解形式如下：

b)求取摄像机内参

真实的摄像机参数矩阵P和实际分解得到的P＇关系如(9)式：

其中T为一非奇异矩阵，可按如下方式确定：

对第一个摄像机来说，R₁为单位阵，t₁的各分量均为0，则P₁＝A₁R₁ ^T[I-t₁]＝[K₁　0]，同时令分解后的矩阵

具有形式 $P_1^{\′}=\left[\begin{array}{cc}I& 0\end{array}\right],$ 于是可得 $T=\left[\begin{array}{cc}{K}_1^{-1}& 0\\ h^T& h\end{array}\right],$ 其中h^T是一个3×1的矢量，h是一个比例因子，

设式(9)中的比例因子为β^j，则式(9)可写为：TQ^j′＝β^jQ^j，

对每个标定板对应的Q^j，由式(2)中的p，q的正交性得出式(10)：

p^j′Tω₁p^j′＝q^j′Tω₁q^j′＝(β^j)²，p^j′Tω₁q^j′＝⁰　　　　(10)

其中 ${\mathrm{\ω}}_1=K_1^{-T}{K}_1^{-1},$ 按式(10)求出第一个摄像机能拍摄到的所有标定板的正交方程，并叠加所有方程，求出K₁，从而得：

${\mathrm{\β}}^j=\sqrt{(p^{j\′T}{\mathrm{\ω}}_1{p}^{j\′}+q^{j\′T}{\mathrm{\ω}}_1{q}^{j\′})/2}---\left(11\right)$

又，每一个标定板均满足方程(12)：

[h^T　h]Q^j′＝[0　0　β^j]　　(12)

叠加上述方程(12)，得[h^T　h]，进而获得变换矩阵T，于是有第一个摄像机参数矩阵 $P_1=P_1^{\′}{T}^{-1},$ 其它摄像机的参数也可按 $P_i=P_i^{\′}{T}^{-1}$ 获得。

c)求取摄像机的外参

摄像机外参数按如下方式求得：

P_i＝A_iR_i ^T[I-t_i]＝[M|b]      (13)

于是： $R_i^{\′}=M^T{A}_i^{-T},$ 并对

做SVD分解后求得旋转矩阵R_i如式(14)，

$R_i^{\′}={\mathrm{USV}}^T\⇒R_i={\mathrm{UV}}^T---\left(14\right)$

摄像机间的位移：t_i＝-M^-1 _ib　　　　　　(15)。

本发明中，所说的对所得内外参数进行非线性优化的步骤：

由式(16)进行非线性优化：

$\min \underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{j=1,i}{\overset{J}{\mathrm{\Σ}}}\underset{n}{\overset{\mathrm{Nj}}{\mathrm{\Σ}}}{\left|\right|m_{\mathrm{in}}^j-\hat{m}(K_i,R_i,t_i,Q^j,M_n^j)\left|\right|}^2---\left(16\right)$

式中，是将标定板上某一特征点按求得的摄像机模型重新投影到图像坐标系中的坐标

是

在实际拍摄的图像中的位置坐标。

本发明适用于以任意方式排列的摄像机阵列。很好的解决了当摄像机视野差别很大或某些图像的质量不能达到摄像机标定的要求时，将测量矩阵正确分解为摄像机参数矩阵和标定板参数矩阵的实现方法，进而精确的求出摄像机的内外参数。此外，通过对矩阵分解形式的观察，我们还提出了通过迭代有效的调整摄像机单应性矩阵的比例因子，来保证测量矩阵的秩4约束，并优化标定结果。同目前已有的其它摄像机阵列标定方法相比，本发明具有求得的摄像机参数之间的一致性较好，精度高，抗噪声性能强，实现简单，容易，应用广泛等诸多优点。

附图说明

图1是以密集方式排列的摄像机阵列的示意图；

图2是以稀疏方式排列的摄像机阵列的示意图；

图3是摄像机参数矩阵和标定板参数矩阵示意图；

图4是摄像机焦距计算误差与噪声系数关系图；

图5是摄像机中心点u0计算误差与噪声系数关系图；

图6是摄像机中心点v0计算误差与噪声系数关系图；

图7是摄像机中心点u0计算误差与噪声系数关系图；

图8是摄像机中心点v0计算误差与噪声系数关系图；

具体实施方式

以下结合实施例进一步说明本发明。

基于矩阵分解的摄像机阵列标定方法，包括以下步骤：

1.假设有六台摄像机对标定板在24个不同位置拍照(如图1所示)，同时将所得的图像也按照拍照顺序及其对应的摄像机编号；

在后面的表示中我们将用下标i表示摄像机编号，上标j表示标定板图像的编号。同时每个摄像机至少拍摄到4幅标定板的图像。

2.利用标定板上特征点的实际物理坐标数据和拍照获得的图像数据求取摄像机的单应性矩阵

a)摄像机的成像模型如式(1)

其中，X表示空间某点在世界坐标系中的坐标，(u　v)为该空间点投影到摄像机图像坐标系中的坐标，A_i表示第i个摄像机的内参矩阵，R_i，t_i为摄像机坐标系在世界坐标系中的旋转和位移，指等式两边相差一个比例因子。

标定板上特征点与世界坐标系之间的转换关系如式(2)

式(2)中p^j、q^j分别为标定板坐标系的单位矢量，d^j为标定板与世界坐标系之间的距离，以上两式中，P_i为摄像机参数矩阵，Q^j为标定板的位置参数矩阵；其具体物理意义可参照图2，图中，x，y为标定板坐标系，x_i、y_i、z_i为摄像机i的摄像机坐标系，x₀、y₀、z₀为第一个摄像机的摄像机坐标系并与世界坐标系重合。

b)我们使用的标定板具有12行12列，共144个特征点，对这144特征点在标定板上的实际坐标(x　y)及其摄像机上的成像坐标数据(u　v)做规范化处理后，(规范化的步骤请参考文献：“In defense of the 8-point algorithm”.R.Hartley.In Proc.5th International Conference on Computer Vision，pages 1064-1070，Boston，MA，June 995；)按式(3)求取摄像机的单应性矩阵 ${\tilde{H}}_i^j={\mathrm{\λ}}_i^j{H}_i^j,$

其中

表示第i个摄像机和第j个标定板图像之间的真实单应性矩阵。

3.将求得的单应性矩阵按照其对应的摄像机和标定板的编号构造测量矩阵

将步骤2)中得到的摄像机和标定板之间的单应性矩阵

依照其上下标按下式所示方式排列，得到测量矩阵

即矩阵的“行标”表示摄像机的标号，“列标”表示标定板的标号。注意：每一个

都是一个(3×3)的矩阵。

$\tilde{W}=\left[\begin{array}{ccc}\widetilde{H_1^1}& \·\·\·& \widetilde{H_1^J}\\ \·& & \·\\ \·& & \·\\ \·& & \·\\ \widetilde{H_I^1}& \·\·\·& \widetilde{H_I^J}\end{array}\right]---\left(4\right)$

其中得

表示实际计算出来的单应性矩阵与真实的单应性矩阵

相差一个比例因子即： ${\tilde{H}}_i^j={\mathrm{\λ}}_i^j{H}_i^j.$

将因为标定板不在摄像机视野范围而无法计算的单应性矩阵按空矩阵对待，如图1中，摄像机1看不到位置在标定板2处的标定板，则为空矩阵，不填任何数据。

4.根据摄像机阵列的密集或稀疏排列情况，调整单应性矩阵的比例因子，并对丢失数据的测量矩阵进行矩阵分解；

按照式(3)中H和P，Q的关系，真实的测量矩阵W可分解如下：

即对一个3I×3J的测量矩阵矩阵W(I表示摄像机的个数，J表示拍摄的标定板图像的个数)，可以分解为大小为3I×4摄像机参数矩阵P和大小为3J×4标定板参数矩阵Q。

a)当摄像机阵列排列密集时

这里的密集是指至少有一个摄像机能拍摄到所有位置的标定板，并且至少在一个位置上该标定板能同时被所有的摄像机同时拍摄到。此种情况比较适用于平面平行排列的摄像机阵列，如图一所示。Stanford大学的摄像机阵列和CMU大学的摄像机阵列都属于这一类型。

(1)调整单应性矩阵的比例因子

本例中摄像机个数I为6，拍摄的标定板的图像个数J为24，如图1所示。

调整 ${\tilde{H}}_i^j={\mathrm{\λ}}_i^j{H}_i^j$ 中的比例因子

方法如下：

定义相关单应性矩阵：对第l和第f个摄像机来说，它们关于第r和第p个标定板图像的相关单应性矩阵

形式如下：

$N_{\mathrm{l\; f}}^{\mathrm{r\; p}}=H_f^r{\left(H_l^r\right)}^{-1}{H}_l^p{\left(H_f^p\right)}^{-1}$

则计算得出的带有比例因子的相关单应性矩阵

为：

${\tilde{N}}_{\mathrm{l\; f}}^{\mathrm{r\; p}}={\tilde{H}}_f^r{\left({\tilde{H}}_l^r\right)}^{-1}{\tilde{H}}_l^p{\left({\tilde{H}}_f^p\right)}^{-1}=\left({\mathrm{\λ}}_f^r{\mathrm{\λ}}_l^p\right)/\left({\mathrm{\λ}}_l^r{\mathrm{\λ}}_f^p\right)N_{\mathrm{l\; f}}^{\mathrm{r\; p}}---\left(5\right)$

由于具有一个单位二重根，因此

具有一个值为

的二重特征根如式(6)

$u_{\mathrm{l\; l}}^{\mathrm{r\; p}}=\left({\mathrm{\λ}}_f^r{\mathrm{\λ}}_l^p\right)/\left({\mathrm{\λ}}_l^r{\mathrm{\λ}}_f^p\right)---\left(6\right)$

将步骤3)构造的测量矩阵

中的元素按行和列排序，使其前三行对应的摄像机能拍摄到所有的标定板，前三列对应的标定板能被所有的摄像机同时拍摄，并固定此时的测量矩阵中的第一个摄像机和第一个标定板作为计算相关单应性矩阵的参考摄像机和参考标定板，则第一个摄像机和第一个标定板与第i个摄像机的第j个标定板图像的相关单应性矩阵为： ${\tilde{N}}_{i1}^{j1}=u_{i1}^{j1}{N}_{i1}^{j1},$ 且

可通过

的二重特征根计算得出，记为

对中的每一个元素做如下调整，

令 ${\tilde{H}}_i^j=u_i^j{\tilde{H}}_i^j=u_i^j{\mathrm{\λ}}_i^j{H}_i^j=({\mathrm{\λ}}_i^1{\mathrm{\λ}}_1^j/{\mathrm{\λ}}_1^1)P_i{Q}^j---\left(7\right)$

则调整后的

为：

$\tilde{W}=\left[\begin{array}{ccc}\widetilde{H_1^1}& \·\·\·& \widetilde{H_1^J}\\ \·& & \·\\ \·& & \·\\ \·& & \·\\ \widetilde{H_I^1}& \·\·\·& \widetilde{H_I^J}\end{array}\right]$

在密集排列方式下，对于

中的缺失元素(如：

)，不做调整。

(2)对丢失数据的测量矩阵

进行矩阵分解。

按以下三个步骤进行：

1)对一个秩为4，大小为3I×3J的测量矩阵

给定任意初值的3I×4矩阵P₀；

2)在P_k-1已知的情况下，按式 $B{(j,:)}^T=\left(\tilde{W}{(:,j)}^{T}A\right){\left(A^{T}A\right)}^{-1}(j=1\·\·\·3J)$ 求取3J×4矩阵Q_k，并对Q_k的列向量做正交化处理：Q_k＝Q_k×N_k，使得

最小，ij代表对

存在的因子求和；

3)在Q_k已知的情况下，按 $A{(i,:)}^T=\left(\tilde{W}(i,:)B\right){\left(B^{T}B\right)}^{-1}(i=1\·\·\·3I)$ 求取3I×4矩阵P_k，使得 $\underset{\mathrm{ij}}{\mathrm{\Σ}}{|{\tilde{W}}_{\mathrm{ij}}-{\left(P_k{Q}_k^T\right)}_{\mathrm{ij}}|}^2$ 最小，ij代表对

存在的因子求和；

4)重复以上2)，3)两个步骤直到

收敛，并对P，Q逆排序后，得到测量矩阵

的分解形式如式(8)：

$\tilde{W}=(1/{\mathrm{\λ}}_1^1)\left[\begin{array}{c}{\mathrm{\λ}}_1^1{P}_1\\ \·\\ \·\\ \·\\ {\mathrm{\λ}}_I^1{P}_I\end{array}\right]\left[\begin{array}{ccc}{\mathrm{\λ}}_1^1{Q}^{1^T}& \·\·\·& {\mathrm{\λ}}_1^J{Q}^{J^T}\end{array}\right]---\left(8\right)$

当迭代不收敛或者误差

超过阈值时，更换迭代初值P₀重新迭代即可。

b)当摄像机阵列的排列稀疏时

这里的稀疏指的是当摄像机的排列比较零散，或者摄像机之间的旋转角度比较悬殊，造成摄像机的视野差别很大的情况。主要应用如视频监控系统，MIT的3D ROOM等。如图3所示，8个摄像机围成一圈，此时，标定板不论摆放在什么位置，都将存在某几个摄像机不能拍摄到该标定板的正面，则这些摄像机和标定板之间对应的单应性矩阵

就不能求出。

(1)调整单应性矩阵的比例因子

本例中摄像机个数I为8，拍摄的标定板的图像个数J为60，如图3所示，拍照过程中，应尽量使得每个摄像机能拍摄到的标定板正面板的数目比较均匀。

对步骤3)中构造的测量矩阵

重排序，使得

的前三行对应的摄像机能拍摄到最多数目的标定板图像，前三列对应的标定板能被最多数目的摄像机同时拍摄到，并固定此时的测量矩阵中的第一个摄像机和第一个标定板作为计算相关单应性矩阵的参考摄像机和参考标定板，则第一个摄像机和第一个标定板与第i个摄像机的第j个标定板图像的相关单应性矩阵为： ${\tilde{N}}_{l1}^{r1}=u_{l1}^{r1}{N}_{l1}^{r1},$ 且通过

的二重特征根计算得出，记

为

对

中前三行和前三列不为空的

和

对应的那些行和列的

做如下调整，

令 ${\tilde{H}}_i^j=u_i^j{\tilde{H}}_i^j=u_i^j{\mathrm{\λ}}_i^j{H}_i^j=({\mathrm{\λ}}_i^1{\mathrm{\λ}}_1^j/{\mathrm{\λ}}_1^1)P_i{Q}^j---\left(7\right)$

因稀疏排列的限制，对无法按照式(7)调整的

的比例因子，我们利用测量矩阵

形式上的特点在矩阵分解迭代过程中获得。

(2)丢失数据的测量矩阵

的分解

按以下几个步骤进行：

1)对一个秩为4，大小为3I×3J的矩阵

给定任意初值的3I×4矩阵P₀；

2)设每一个

的比例因子 ${\mathrm{\λ}}_i^j=1;$

3)利用当前的

构建测量矩阵

4)在P_k-1已知的情况下，按式 $Q_k{(j,:)}^T=\left(\tilde{W}{(:,j)}^T{P}_{k-1}\right){\left({P_{k-1}}^T{P}_{k-1}\right)}^{-1}(j=1\·\·\·3J)$ 求取3J×4矩阵Q_k，并对Q_k的列向量做正交化处理：Q_k＝Q_k×N_k，使得

最小，ij代表对存在的因子求和；

5)在Q_k已知的情况下，按式 $P_k{(i,:)}^T=\left(\tilde{W}(i,:)Q_k\right){\left({Q_k}^T{Q}_k\right)}^{-1}(i=1\·\·\·3I)$ 求取3I×4矩阵P_k，使得

最小，ij代表对

存在的因子求和；

6)调整

中前三行和前三列为空的

和

对应的那些行和列的

的比例因子如下： ${\mathrm{\λ}}_i^j=\mathrm{\Σ}{P}_i(3,:)\×Q_j{(3,:)}^T,$ 其中P_i(3，：)代表P_i的第三行，Q_j(3，：)代表Q^j的第三列；

7)重复步骤3)到6)直到收敛，并对P，Q逆排序后，得到的分解形式如式(8)：

$\tilde{W}=(1/{\mathrm{\λ}}_1^1)\left[\begin{array}{c}{\mathrm{\λ}}_1^1{P}_1\\ \·\\ \·\\ \·\\ {\mathrm{\λ}}_I^1{P}_I\end{array}\right]\left[\begin{array}{ccc}{\mathrm{\λ}}_1^1{Q}^{1^T}& \·\·\·& {\mathrm{\λ}}_1^J{Q}^{J^T}\end{array}\right]---\left(8\right)$

当迭代不收敛或者误差

超过阈值时，更换迭代初值P₀重新迭代即可。

5.添加度量约束，求解摄像机的内外参数

a)对P，Q去规范化，得到测量矩阵

分解形式如下：

b)求取摄像机内参

真实的摄像机参数矩阵P和实际分解得到的P′关系如(9)式：

  

其中T为一非奇异矩阵，可按如下方式确定：

对第一个摄像机来说，R₁为单位阵，t₁的各分量均为0，则P₁＝A₁R₁ ^T[I-t₁]＝[K₁ 0]，同时令分解后的矩阵

具有形式 $P_1^{\′}=\left[\begin{array}{cc}I& 0\end{array}\right],$ 于是可得 $T=\left[\begin{array}{cc}{K}_1^{-1}& 0\\ h^T& h\end{array}\right],$ 其中h^T是一个3×1的矢量，h是一个比例因子。

设式(9)中的比例因子为β^j，则式(9)可写为：TQ^j′＝β^jQ^j。

对每个标定板对应的Q^j，由式(2)中的p，q的正交性得出式(10)：

p^j′Tω₁p^j′＝q^j′Tω₁q^j′＝(β^j)²，p^j′Tω₁q^j′＝0　　　　　(10)

其中 ${\mathrm{\ω}}_1=K_1^{-T}{K}_1^{-1},$ 按式(10)求出第一个摄像机能拍摄到的所有标定板的正交方程，并叠加所有方程，求出K₁，从而得：

${\mathrm{\β}}^j=\sqrt{(p^{j\′T}{\mathrm{\ω}}_1{p}^{j\′}+q^{j\′T}{\mathrm{\ω}}_1{q}^{j\′})/2}---\left(11\right)$

又，每一个标定板均满足方程(12)：

[h^T　h]Q^j′＝[0 0 β^j]    (12)

叠加上述方程(12)，得[h^T h]，进而获到变换矩阵T，于是有第一个摄像机参数矩阵 $P_1=P_1^{\′}{T}^{-1},$ 其它摄像机的参数也可按 $P_i=P_i^{\′}{T}^{-1}$ 获得。

c)求取摄像机的外参

摄像机外参按如下方式求取：

P_i＝A_iR_i ^T[I-t_i]＝[M|b]           (13)

于是有： $R_i^{\′}=M^T{A}_i^{-T},$ 并对

做SVD分解后求得旋转矩阵R_i如式(14)，

$R_i^{\′}={\mathrm{USV}}^T\⇒R_i={\mathrm{UV}}^T---\left(14\right)$

摄像机间的位移：t_i＝-M^-1 _ib。            (15)

6.非线性优化

为进一步提高计算得出的摄像机内外参数的准确性，对这些参数按式(16)进行非线性优化：

$\min \underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{j=1,i}{\overset{J}{\mathrm{\Σ}}}\underset{n}{\overset{\mathrm{Nj}}{\mathrm{\Σ}}}{\left|\right|m_{\mathrm{in}}^j-\hat{m}(K_i,R_i,t_i,Q^j,M_n^j)\left|\right|}^2---\left(16\right)$

式中，

是将标定板上某一特征点

按求得的摄像机模型重新投影到图像坐标系中的坐标，

是

在实际拍摄的图像中的位置坐标。

图4到图8是对本发明的方法进行仿真后的误差性能图。仿真中所使用的摄像机阵列包含6个摄像机，其中，每次实验中摄像机的内外参数设置如下：

1)摄像机外参数设置：摄像机与摄像机在x方向的间隔在(300mm，600mm)之间均匀分布，y，z方向的间隔在(-5mm，5mm)之间均匀分布，各摄像机关于世界坐标系(x，y，z)坐标轴的旋转角度也在(0°，18°)之间均匀分布。

2)摄像机内参数设置：各摄像机的焦距符合均值为790，方差为20的正态分布，摄像机的中心点坐标为(300，300)。

仿真中，标定板约在摄像机阵列500mm前自由移动，其关于(x，y，z)各坐标轴的旋转角度在(0°，60°)间均匀分布。同时，为了测试本发明方法的信噪比，实验中，给每个图像点叠加均值为0，方差从0到2之间以0.2为步长递进的高斯噪声，并对每一个噪声系数做100次实验取平均结果。

图4为摄像机的焦距误差和噪声系数的关系曲线，图中分别只列出了第一，三，五个摄像机各自焦距的平均误差以及所有摄像机的累计平均焦距误差与噪声系数之间的关系曲线，其中x轴为叠加的噪声系数值，y轴为摄像机焦距误差值。

图5为摄像机中心点横坐标u0的误差和噪声系数的关系曲线，图中只分别列出了第一，三，五个摄像机各自横坐标u0的平均误差以及所有摄像机的累计平均u0误差同仿真中的噪声系数的关系曲线，其中x轴为叠加的噪声系数值，y轴为摄像机横坐标u0误差值。

图6为摄像机中心点纵坐标v0的误差和噪声系数的关系曲线，图中只分别列出了第一，三，五个摄像机各自纵坐标v0的平均误差以及所有摄像机的累计平均v0误差同仿真中的噪声系数的关系曲线，其中x轴为叠加的噪声系数值，y轴为摄像机纵坐标v0误差值。

图7为摄像机相互距离位置关系误差和噪声系数的关系曲线，图中只分别列出了第二，三，五个摄像机各自同第一个摄像机的相互距离位置误差以及所有摄像机同第一个摄像机的相互距离位置误差的累计平均误差与仿真中的噪声系数的关系曲线，其中x轴为叠加的噪声系数值，y轴为摄像机相互距离误差值。

图8为摄像机相互间旋转角度误差和噪声系数的关系曲线，图中只分别列出了第二，三，五个摄像机各自同第一个摄像机的旋转角度误差以及所有摄像机同第一个摄像机的旋转角度误差的累计平均误差与仿真中的噪声系数的关系曲线，其中x轴为叠加的噪声系数值，y轴为摄像机相互间旋转角度误差值。

由图4到图8所示的误差性能曲线可见，利用本发明方法所得的各摄像机内外参数准确性较高，抗噪声能力强，能够满足对摄像机阵列进行标定的要求。

Claims (6)

1.基于矩阵分解的摄像机阵列标定方法，其特征在于包括以下步骤：

1)将摄像机编号并对标定板拍照，同时将所得的图像也按照拍照顺序及其对应的摄像机编号；

2)利用标定板上特征点的实际物理坐标数据和拍照获得的图像数据求取摄像机的单应性矩阵

；

3)将求得的单应性矩阵按照其对应的摄像机和标定板的编号构造测量矩阵

；

4)根据摄像机阵列的密集或稀疏排列情况，调整单应性矩阵的比例因子，并对丢失数据的测量矩阵

进行矩阵分解；

5)添加度量约束，求解摄像机的内外参数；

6)对所得内外参数进行非线性优化。

2.根据权利要求1所述的基于矩阵分解的摄像机阵列标定方法，其特征在于所说的利用标定板上特征点的实际物理坐标数据和拍照获得的图像数据求取摄像机的单应性矩阵的步骤如下：

a)摄像机的成像模型如式(1)

其中，X表示空间某点在世界坐标系中的坐标，(u v)为该空间点投影到摄像机图像坐标系中的坐标，A_i表示第i个摄像机的内参矩阵，R_i，t_i为摄像机坐标系在世界坐标系中的旋转和位移，

指等式两边相差一个比例因子，

标定板上特征点坐标与世界坐标系之间的转换关系如式(2)

以上两式中，P_i为摄像机参数矩阵，Q^j为标定板的位置参数矩阵；

b)对图像数据和已知的标定板上特征点的实际坐标数据做规范化处理，并按式(3)求取摄像机的单应性矩阵 ${\tilde{H}}_i^j={\mathrm{\λ}}_i^j{H}_i^j,$

其中的

表示第i个摄像机和第j个标定板图像之间的真实单应性矩阵。

3.根据权利要求1所述的基于矩阵分解的摄像机阵列标定方法，其特征在于所说的将求得的单应性矩阵按照其对应的摄像机和标定板的编号构造测量矩阵的方法是，将步骤2)中得到的摄像机和标定板之间的单应性矩阵

依照其上下标按式(4)所示方式排列，得到测量矩阵

即矩阵的“行标”表示摄像机的标号，“列标”表示标定板图像的标号，每一个

都是一个(3×3)的矩阵，

$\tilde{W}=\left[\begin{array}{ccc}\widetilde{H_1^1}& \·\·\·& \widetilde{H_1^J}\\ \·& & \·\\ \·& & \·\\ \·& & \·\\ \widetilde{H_I^1}& \·\·\·& \widetilde{H_I^J}\end{array}\right]---\left(4\right)$

其中得

表示实际计算出来的单应性矩阵与真实的单应性矩阵

相差一个比例因子

即： ${\tilde{H}}_i^j={\mathrm{\λ}}_i^j{H}_i^j.$

4.根据权利要求1所述的基于矩阵分解的摄像机阵列标定方法，其特征在于所说的根据摄像机阵列的密集或稀疏排列情况，调整单应性矩阵的比例因子，并对丢失数据的测量矩阵进行矩阵分解，其步骤按密集或稀疏排列情况分别如下：

a)当摄像机阵列排列密集时

(1)调整单应性矩阵的比例因子

调整 ${\tilde{H}}_i^j={\mathrm{\λ}}_i^j{H}_i^j$ 中的比例因子方法如下：

对第l和第f个摄像机来说，它们分别关于第r和第p个标定板图像的相关单应性矩阵

形式如下：

${N_l^r}_f^p=H_f^r{\left(H_l^r\right)}^{-1}{H}_l^p{\left(H_f^p\right)}^{-1}$

则计算得出的带有比例因子的相关单应性矩阵

为：

${{\tilde{N}}_l^r}_f^p={\tilde{H}}_f^r{\left({\tilde{H}}_l^r\right)}^{-1}{\tilde{H}}_l^p{\left({\tilde{H}}_f^p\right)}^{-1}=\left({\mathrm{\λ}}_f^r{\mathrm{\λ}}_l^p\right)/\left({\mathrm{\λ}}_l^r{\mathrm{\λ}}_f^p\right){N_l^r}_f^p---\left(5\right)$

具有一个值为

的二重特征根如式(6)

${u_l^r}_l^p=\left({\mathrm{\λ}}_f^r{\mathrm{\λ}}_l^p\right)/\left({\mathrm{\λ}}_l^r{\mathrm{\λ}}_f^p\right)---\left(6\right)$

将步骤3)构造的测量矩阵

中的元素按行和列排序，使其前三行对应的摄像机能拍摄到所有的标定板，前三列对应的标定板能被所有的摄像机同时拍摄，并固定此时的测量矩阵中的第一个摄像机和第一个标定板作为计算相关单应性矩阵的参考摄像机和参考标定板，则第一个摄像机和第一个标定板与第i个摄像机的第j个标定板图像的相关单应性矩阵为： ${\tilde{N}}_{l1}^{r1}=u_{l1}^{r1}{N}_{l1}^{r1},$ 且

通过

的二重特征根计算得出，记

对

中的每一个元素做如下调整，

令 ${\tilde{H}}_i^j=u_i^j{\tilde{H}}_i^j=u_i^j{\mathrm{\λ}}_i^j{H}_i^j=({\mathrm{\λ}}_i^1{\mathrm{\λ}}_1^j/{\mathrm{\λ}}_1^1)P_i{Q}^j---\left(7\right)$

则调整后的为：

$\tilde{W}=\left[\begin{array}{ccc}\widetilde{H_1^1}& \·\·\·& \widetilde{H_1^J}\\ \·& & \·\\ \·& & \·\\ \·& & \·\\ \widetilde{H_I^1}& \·\·\·& \widetilde{H_I^J}\end{array}\right]$

(2)丢失数据的测量矩阵

的分解

按以下三个步骤进行：

1)对一个秩为4，大小为3I×3J的测量矩阵

给定任意初值的3I×4矩阵P₀；

2)在P_k-1已知的情况下，按式 $B{(j,:)}^T=\left(\tilde{W}{(:,j)}^{T}A\right){\left(A^{T}A\right)}^{-1}(j=1\·\·\·3J)$ 求取3J×4矩阵Q_k，并对Q_k的列向量做正交化处理：Q_k＝Q_k×N_k，使得

最小，ij代表对存在的因子求和；

3)在Q_k已知的情况下，按式 $A{(i,:)}^T=\left(\tilde{W}(i,:)B\right){\left(B^{T}B\right)}^{-1}(i=1\·\·\·3I)$ 求取3I×4矩阵P_k，使得

最小，ij代表对

存在的因子求和；

4)重复以上2)，3)两个步骤直到

收敛，并对P，Q逆排序后，得到测量矩阵

的分解形式如式(8)：

$\tilde{W}=(1/{{\mathrm{\λ}}_1}^1)\left[\begin{array}{c}{{\mathrm{\λ}}_1}^1{P}_1\\ \·\\ \·\\ \·\\ {\mathrm{\λ}}_I^1{P}_I\end{array}\right]\left[\begin{array}{ccc}{{\mathrm{\λ}}_1}^1{Q}^{1^T}& \·\·\·& {{\mathrm{\λ}}_1}^J{Q}^{J^T}\end{array}\right]---\left(8\right)$

当迭代不收敛或者误差

超过阈值时，更换迭代初值P₀重新迭代即可。

b)当摄像机阵列的排列稀疏时

(1)对步骤3)中构造的测量矩阵

重排序，使得

的前三行对应的摄像机能拍摄到最多数目的标定板图像，前三列对应的标定板能被最多数目的摄像机同时拍摄到，并固定此时的测量矩阵中的第一个摄像机和第一个标定板作为计算相关单应性矩阵的参考摄像机和参考标定板，则第一个摄像机和第一个标定板与第i个摄像机的第j个标定板图像的相关单应性矩阵为： ${\tilde{N}}_{l1}^{r1}=u_{l1}^{r1}{N}_{l1}^{r1},$ 且

通过

的二重特征根计算得出，记为

对

中前三行和前三列不为空的

和

对应的那些行和列的

做如下调整，

令 ${\tilde{H}}_i^j=u_i^j{\tilde{H}}_i^j=u_i^j{\mathrm{\λ}}_i^j{H}_i^j=({\mathrm{\λ}}_i^1{\mathrm{\λ}}_1^j/{\mathrm{\λ}}_1^1)P_i{Q}^j---\left(7\right)$

则调整后的为：

$\tilde{W}=\left[\begin{array}{ccc}\widetilde{H_1^1}& \·\·\·& \widetilde{H_1^J}\\ \·& & \·\\ \·& & \·\\ \·& & \·\\ \widetilde{H_I^1}& \·\·\·& \widetilde{H_I^J}\end{array}\right]$

(2)丢失数据的测量矩阵

的分解

按以下几个步骤进行：

1)对一个秩为4，大小为3I×3J的矩阵

给定任意初值的3I×4矩阵P₀；

2)设每一个

的比例因子 ${\mathrm{\λ}}_i^j=1;$

3)利用当前的

构建测量矩阵

4)在P_k-1已知的情况下，按式 $Q_k{(j,:)}^T=\left(\tilde{W}{(:,j)}^T{P}_{k-1}\right){\left({P_{k-1}}^T{P}_{k-1}\right)}^{-1}(j=1\·\·\·3J)$ 求取3J×4矩阵Q_k，并对Q_k的列向量做正交化处理：Q_k＝Q_k×N_k，使得

最小，ij代表对

存在的因子求和；

5)在Q_k已知的情况下，按式 $P_k{(i,:)}^T=\left(\tilde{W}(i,:)Q_k\right){\left({Q_k}^T{Q}_k\right)}^{-1}(i=1\·\·\·3I)$ 求取3I×4矩阵P_k，使得

最小，ij代表对

存在的因子求和；

6)调整

中前三行和前三列为空的

和对应的那些行和列的

的比例因子如下： ${\mathrm{\λ}}_i^j=\mathrm{\Σ}{P}_i(3,:)\×Q_j{(3,:)}^T,$ 其中P_i(3，：)代表P_i的第三行，Q_j(3，：)代表Q^j的第三列；

7)重复步骤3)到6)直到

收敛，并对P，Q逆排序后，得到的分解形式如式(8)：

$\tilde{W}=(1/{{\mathrm{\λ}}_1}^1)\left[\begin{array}{c}{{\mathrm{\λ}}_1}^1{P}_1\\ \·\\ \·\\ \·\\ {\mathrm{\λ}}_I^1{P}_I\end{array}\right]\left[\begin{array}{ccc}{{\mathrm{\λ}}_1}^1{Q}^{1^T}& \·\·\·& {{\mathrm{\λ}}_1}^J{Q}^{J^T}\end{array}\right]---\left(8\right)$

当迭代不收敛或者误差

超过阈值时，更换迭代初值P₀重新迭代即可。

5.根据权利要求1所述的基于矩阵分解的摄像机阵列标定方法，其特征在于所说的添加度量约束，求解摄像机的内外参数的步骤如下：

a)对P，Q去规范化，得到测量矩阵

分解形式如下：

b)求取摄像机内参

真实的摄像机参数矩阵P和实际分解得到的P′关系如(9)式：

    

其中T为一非奇异矩阵，可按如下方式确定：

对第一个摄像机来说，R₁为单位阵，t₁的各分量均为0，则 $P_1=A_1{R_1}^T[I-t_1]=\left[\begin{array}{cc}{K}_1& 0\end{array}\right],$ 同时令分解后的矩阵

具有形式 $P_1^{\′}=\left[\begin{array}{cc}I& 0\end{array}\right],$ 于是可得 $T=\left[\begin{array}{cc}{K}_1^{-1}& 0\\ h^T& h\end{array}\right],$ 其中h^T是一个3×1的矢量，h是一个比例因子，

设式(9)中的比例因子为β^j，则式(9)可写为：TQ^j′＝β^jQ^j，

对每个标定板对应的Q^j，由式(2)中的p，q的正交性得出式(10)：

p^j′Tω₁p^j′＝q^j′Tω₁q^j′＝(β^j)²，p^j′Tω₁q^j′＝0            (10)

其中 ${\mathrm{\ω}}_1=K_1^{-T}{K}_1^{-1},$ 按式(10)求出第一个摄像机能拍摄到的所有标定板的正交方程，并叠加所有方程，求出K₁，从而得：

${\mathrm{\β}}^j=\sqrt{(p^{j\′T}{\mathrm{\ω}}_1{p}^{j\′}+q^{j\′T}{\mathrm{\ω}}_1{q}^{j\′})/2}---\left(11\right)$

又，每一个标定板均满足方程(12)：

[h^T h]Q^j′＝[0 0 β^j]           (12)

叠加上述方程(12)，得[h^T h]，进而获得变换矩阵T，于是有第一个摄像机参数矩阵 $P_1=P_1^{\′}{T}^{-1},$ 其它摄像机的参数也可按 $P_i=P_i^{\′}{T}^{-1}$ 获得。

c)求取摄像机的外参

摄像机外参数按如下方式求得：

$P_i=A_i{R}_i^T[I-t_i]=\left[M|b\right]---\left(13\right)$

于是： $R_i^{\′}=M^T{A}_i^{-T},$ 并对

做SVD分解后求得旋转矩阵R_i如式(14)，

       $R_i^{\′}={\mathrm{USV}}^T\⇒R_i={\mathrm{UV}}^T---\left(14\right)$

摄像机间的位移：t_i＝-M^-1 _ib     (15)。

6.根据权利要求1所述的一种基于矩阵分解的摄像机阵列标定方法，其特征在于所说的对所得内外参数进行非线性优化的步骤：

由式(16)进行非线性优化：

$\min \underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{j=1,i}{\overset{J}{\mathrm{\Σ}}}\underset{n}{\overset{\mathrm{Nj}}{\mathrm{\Σ}}}{\left|\right|m_{\mathrm{in}}^j-\hat{m}(K_i,R_i,t_i,Q^j,M_n^j)\left|\right|}^2---\left(16\right)$

式中，

是将标定板上某一特征点

按求得的摄像机模型重新投影到图像坐标系中的坐标，

是

在实际拍摄的图像中的位置坐标。

Images