可变参数的相机标定

1、可变参数的初标定

在进行摄像机拍摄时,可能会对焦距进行改变从而影响相机内参的确定,我们把影响相机内参的标定称为可变参数标定,相机的标定是为了获取相机内外参数,消除畸变.

可变参数的标定主要采用平面模板标定法,但是平面模板标定需要一个比较好的初值,若初值寻则不当,则算法难以收敛或者只能收敛到局部最小,从而大大降低标定的精度.在采用平面模板标定法时,虽先不考虑各种畸变,而是先对所有点代入求解,但由于远离图像中心的像点畸变很大,如将这些点也看做没有畸变的像点代入，显然会加大求解初值的误差。新算法考虑到图像中心附近点的畸变很小，因此可以先利用图像中心附近点求取初值，接下来全面考虑各可变参数的非线性优化计算时就能很快收敛。

1.1单应性矩阵求解

单应性矩阵是同射影几何中的直射变换的概念是一致的，表示从空间Pⁿ到Pⁿ的投影变换，习惯将二维欧式空间的直射变换称为单应性。

由摄像机的针孔模型知，图像坐标与世界坐标的转换可以表示为

${ m{s}}left[ {egin{array}{*{20}{c}}
u\
v\
1
end{array}} ight] = { m{A}}left[ {{ m{R}},{ m{T}}} ight]left[ {egin{array}{*{20}{c}}
{{X_w}}\
{{Y_w}}\
{{Z_w}}\
1
end{array}} ight]$

式中，s是任意比例因子，A是摄像机内参矩阵，[R,T]是外参矩阵。

对于平面模型而言，我们假定世界坐标系的Z_w分量为0，并令r_i表示旋转矩阵R的列向量，则理想的针孔透视模型可以写为

${ m{s}}left[ {egin{array}{*{20}{c}}
u\
v\
1
end{array}} ight] = { m{A}}left[ {{{ m{r}}_1},{{ m{r}}_2},{{ m{r}}_3},{ m{T}}} ight]left[ {egin{array}{*{20}{c}}
{{X_w}}\
{{Y_w}}\
0\
1
end{array}} ight] = { m{A}}left[ {{{ m{r}}_1},{{ m{r}}_2},{ m{T}}} ight]left[ {egin{array}{*{20}{c}}
{{X_w}}\
{{Y_w}}\
1
end{array}} ight]$

（1.1）

在式1.1中，令m=[u,v]^T,${ m{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}} over m} }}$=[u,v,1]^T,分别表示空间内点的坐标及其齐次坐标（n维变成n+1维，用于将笛卡尔坐标系转换到透视坐标系，其次坐标常用于坐标系的转换），令M=[X_w,Y_w]^T,${ m{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}} over M} }}$=[X_w,Y_w,1]^T表示空间点的世界坐标系及其齐次坐标，则图像平面上的点${ m{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}} over m} }}$和世界坐标（平面模板上的点）${ m{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}} over M} }}$之间就可以用单应性矩阵H来联系

$smathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over m} = Hmathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} $

（1.2）

其中单应性矩阵H=A[r₁,r₂,T]是一个3x3的矩阵。因为平移向量T是从世界坐标系原点到摄像机光心的矢量, r₁,r₂是图像平面两个坐标轴在世界坐标系中的方向矢量,,故T与r₁,r₂不在一个平面上，且由于旋转矩阵R的正交性，所以r₁与r₂正交，由此可得det(H)≠0。

如何求解H？

通过图像处理获得的空间点的世界坐标系坐标与图像坐标理论上能满足式1.2，但是由于噪声点的存在，式1.2并不能完全得到满足，那么H的极大似然估计应该满足最小二乘估计

${ m{min}}mathop sum limits_i {left| {{m_i} - {{hat m}_i}} ight|^2}$   （1.3）

式1.3中${hat m_i} = left[ {frac{{{{hat h}_1}{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }_i}}}{{{{hat h}_3}{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }_i}}},frac{{{{hat h}_2}{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }_i}}}{{{{hat h}_3}{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }_i}}}} ight]$,是根据式1.2消去s后得到的图像坐标${ m{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}} over m} }}$的期望值（参考齐次矩阵归一化），${hat h_i}$表示矩阵H的第i个行向量。

这是一个非线性最小二乘问题，可以用最小二乘法进行非线性优化计算，用非线性优化方法求解需要有一个合适的初始值来迭代，可通过下面的方法解决。

${hat h_3}mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} left[ {egin{array}{*{20}{c}}
u\
v
end{array}} ight] - left[ {egin{array}{*{20}{c}}
{{{hat h}_1}{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }_i}}\
{{{hat h}_2}{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }_i}}
end{array}} ight] = 0$

（1.4）

或进一步改写成

$left{ {egin{array}{*{20}{c}}
{{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }^T}h_1^T - u{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }^T}h_3^T = 0}\
{{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }^T}h_2^T - v{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }^T}h_3^T = 0}
end{array}} ight.$

（1.5）

令X=[h₁^T，h₂^T ，h₃^T]^T，从而式1.5改写为

[left[ {egin{array}{*{20}{c}}
{{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }^T}}&{{0^T}}&{u{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }^T}}\
{{0^T}}&{{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }^T}}&{v{{mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} }^T}}
end{array}} ight]{ m{X}} = 0]

（1.6）

如果有n对给定的点我们就能得到n个上述方程,联立这n个方程组,我们会得到形如LX=0的矩形方程，其中L是一个2n×9的矩阵。显然这是一个齐次线性方程组，且X具有8个未知数，那么就需要4对点，若n＞4就形成了一个超定方程组，这个方程的求解就变成了一个最优问题。S这类问题通常采用奇异值分解的方法来进行（SVD）,SVD分解经常用来解决最小平方误差问题，其解是L矩阵的最小奇异值对应的右奇异向量（或者矩阵L^TL的最小特征值对应的特征向量）。

这里首先采用SVD分解法,若n>4，则将SVD分解得到的解作为初值，再用最小二乘法优化。

在数学上，奇异值分解实际上是“对称矩阵正交相似于对角矩阵”的推广，它的基本结论是，设任意一个实矩阵A，A^T表示A的转置矩阵，r(A)表示A的秩，则

（1）   A的特征值为非负数

（2）   存在M阶正交阵U、m×n广义对角阵∑、n阶正交阵V，使A=U∑V^T。

其中（1）可以直接验证得到，（2）的证明思路为：设r(A)=r，A^TA为n阶对称阵，由（1）知，A^TA的特征值为l₁≥l₂≥…≥l_r≥0,实际上l₁≥l₂≥…≥l_r>0，l_r+1=…=l_n=0。令，${delta _i} = sqrt {{lambda _i}} $（i=1,2,…,r）${delta _i}$（被称为A的奇异值），由于A^TA为n阶对称阵，故有n阶正交阵V使得${delta _i} = sqrt {{lambda _i}} $（i=1,2,…,r）${delta _i}$（被称为A的奇异值），由于A^TA为n阶对称阵，故有n阶正交阵V使得

${V^T}{A^T}AV = left[ {egin{array}{*{20}{c}}
{{lambda _1}}&{}&{}\
{}& ddots &{}\
{}&{}&{{lambda _n}}
end{array}} ight]$

(1.7)

设V₁,V₂,…,V_n分别为V的列向量，取

${U_i} = frac{1}{{{delta _i}}}A{V_i}$   ,i=1,2,…,r   (1.8)

设U₁,U₂,…,U_r是一单位正交向量组，扩充成R^m的一组标准正交基U₁,U₂,…,U_r,…,U_m并设U=[U₁,U₂,…,U_r,…,U_m],则有A=UDV^T，其中U_i,V_i分别称为矩阵A的左奇异向量和右奇异向量，D为非标准对角矩阵，D的非零对角元是${delta _i}$(i=1,2,…,r)，证明完毕。

针对摄像机标定过程中的实际情况，利用SVD分解可以得到矩阵方程LX=0的最小配置解与超定解。若n=4，L是一个8×9的矩阵，此时L 的奇异值只有8个，而右奇异向量有9个，所以给定这个最小数目的对应点时，能够得到精确解，即空间的点准确地投影到它们被测量的图像上。该解是对应奇异值为零的右奇异向量，这个解称为最小配置解。

对于矩阵L的维数较大的情况(L是一个2n×9的矩阵，n>4)，LX=0将不存在精确解，X的解可以通过最小化一个代数或几何误差来获得。为了减少运算时间并防止计算机内存消耗过多，则针对9×9方阵L^TL作SVD分解得

[U,D,V]=SVD(L^TL)

此时D是标准对角矩阵

${ m{D}} = left[ {egin{array}{*{20}{c}}
{{delta _1}}&{}&{}\
{}& ddots &{}\
{}&{}&{{delta _9}}
end{array}} ight]$

式中，${delta _1}$(i=1,2,…,9)是矩阵L的特征值。满足方程LX=0的解是最小特征值${delta _9}$对应的特征向量V₉，即取出矩阵最后一列即为方程的解，这个解称为超定解。

不过在实际求解过程中，观察L矩阵的构成可知，其元素有的是像素值有的是世界坐标系，有的是这二者的乘积，即构成矩阵L的各元素数量值相差很大，因此L是一个非常病态的方程，需要对其进行改造。研究D的构成发现，${delta _1}/{delta _9}$的结果非常大。根据SVD分解理论可知，${delta _1}/{delta _9}$被称为方程的条件数，而条件数的大小又是线性问题稳定性分析中最关键的一个因素，也即L中的数据即使很小的扰动也可能造成结果很大的变动。可以说正是由于L的条件数很大，从而使得其解得误差可能更大。举例说明，一个典型的图像坐标可能是（100,100,1）左右，即元素大小的数量级为10²~1之间，则矩阵各元素大小的数量级将在10⁴~1之间，而D矩阵对角元素大小的数量级就在10⁸~1之间。解决这个问题方法是对L矩阵对角元素进行SVD分解前对L进行改造，这里采用Hartley的八点算法对矩阵L进行改造。改造方法如下。

坐标原点平移：将空间坐标及其对应点图像坐标的原点从原来位置平移到质心点，并以相对量来表示坐标，这样坐标将以质心点为中心呈均匀分布。所谓质心：就是同维数（对世界坐标来说就是X_w轴与Y_w轴，对图像坐标来说就是u，v轴）所有坐标点的平均值。

坐标归一化处理：有各向同性和各向异性两种方法，对二维系统来说，各向同性是两坐标轴上数据缩放的比例相同，经过变换后各坐标对原点的距离平均值为，各向异性是两坐标轴上数据缩放比例不等，但处理后使得各点到原点的平均距离为1。采用各向同性的好处是在很大程度上可以避免后期数据处理时的迭代不收敛的情况。

数据归一化不仅提高了结果的精度，它还提供了第二个好处，即对初始数据归一化的算法将对任何尺度缩放和坐标原点的选择不变，原因是归一化步骤通过为测量数据选择有效的标准坐标系，预先消除了坐标变换的影响。因此由于代数最小化在一个固定的标准坐标系中进行，使算法实际上关于相似变换不变。

这里对图像坐标和世界坐标进行归一化处理以改善矩阵L的病态性能，设各图像的点坐标为(u_i,v_i)(i=1,2,…,n),各空间点的世界坐标系坐标为(X_wi,Y_wi)(i=1,2,…,n)

a. 分别求出u，v两个轴各自坐标的平均值

${m_u} = frac{{mathop sum olimits_{i = 1}^n {u_i}}}{n}$ ,${m_v} = frac{{mathop sum olimits_{i = 1}^n {v_i}}}{n}$   

(1.9)

则各图像点相对于平均值点的相对值为

$Delta {u_i} = {u_i} - {m_u}$

$Delta {v_i} = {v_i} - {m_v}$

(1.10)

同样对世界坐标系坐标进行处理可得

${m_x} = frac{{mathop sum olimits_{i = 1}^n {X_{wi}}}}{n}$

${m_y} = frac{{mathop sum olimits_{i = 1}^n {Y_{wi}}}}{n}$

$Delta {X_{wi}} = {X_{wi}} - {m_x}$

$Delta {Y_{wi}} = {Y_{wi}} - {m_y}$

(1.11)

b. 分别求出u,v两轴上以及X_wi,Y_wi两轴上所有相对坐标到质心的平均值

${{ m{d}}_{uv}} = frac{{mathop sum olimits_{i = 1}^n sqrt {{{left( {Delta {u_i}} ight)}^2} + {{left( {Delta {v_i}} ight)}^2}} }}{n}$

${{ m{d}}_{xy}} = frac{{mathop sum olimits_{i = 1}^n sqrt {{{left( {Delta {X_{wi}}} ight)}^2} + {{left( {Delta {Y_{wi}}} ight)}^2}} }}{n}$

(1.12)

c. 分别求出u，v两轴上以及X_wi,Y_wi两轴上比例缩放因子

${{ m{s}}_{uv}} = frac{{sqrt 2 }}{{{{ m{d}}_{uv}}}}$

${{ m{s}}_{xy}} = frac{{sqrt 2 }}{{{{ m{d}}_{xy}}}}$

(1.13)

d. u，v两轴上以及X_wi,Y_wi两轴上的变换关系可以用矩阵表示

${T_{uv}} = left[ {egin{array}{*{20}{c}}
{{{ m{s}}_{uv}}}&0&{ - {{ m{s}}_{uv}}{m_u}}\
0&{{{ m{s}}_{uv}}}&{ - {{ m{s}}_{uv}}{m_v}}\
0&0&1
end{array}} ight]$

${T_{xy}} = left[ {egin{array}{*{20}{c}}
{{{ m{s}}_{xy}}}&0&{ - {{ m{s}}_{xy}}{m_x}}\
0&{{{ m{s}}_{xy}}}&{ - {{ m{s}}_{xy}}{m_y}}\
0&0&1
end{array}} ight]$

(1.14)

e. 设经过归一化处理的图像坐标与世界坐标系分别为$mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}} over m} ' = {[egin{array}{*{20}{c}}{Delta u}&{Delta v}&1end{array}]^T}$与$mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}} over M} ' = {[egin{array}{*{20}{c}}{Delta {X_{ m{w}}}}&{Delta {Y_{ m{w}}}}&1end{array}]^T}$,则有

$mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over m} ' = {T_{uv}}mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over m} $

[mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} ' = {T_{xy}}mathord{uildrel{lower3pthbox{$scriptscriptstylesmile$}}
over M} ]

(1.15)

由式1.10和1.11求出经过归一化处理的图像坐标和世界坐标系坐标的单应性矩阵为$ ilde H$，则原坐标系中图像平面和空间坐标系之间的单应性矩阵为

${ m{H}} = T_{uv}^{ - 1} ilde H{T_{xy}}$