CN104376572A - 移动终端相机标定方法 - Google Patents

Abstract

移动终端相机标定方法，包括以下步骤：采用黑白棋盘格图像作为标定物，在标定物上选择多个特征点；利用待标定的移动终端相机拍摄至少2张标定物影像；识别特征点在拍摄的标定物影像上的像方坐标，并确定其对应的物方坐标；根据特征点的像方坐标(u,v)、物方坐标(X,Y,Z)的投影变化关系式解算标定参数：主距分量(f_x,f_y)、偏心畸变系数(c_x,c_y)、径向畸变系数(k₁,k₂,k₃)、切向畸变系数(p₁,p₂)。本发明是专门针对移动终端相机物理特性设计的标定方法，实现移动终端相机标定参数解算，解决传统标定法工作时间长、操作不便、价格昂贵、不适用于移动终端等问题，该方法解算结果精度高、稳定性好，为有效实现移动终端相机标定，获取高质量影像提供了有力保证。

Description

移动终端相机标定方法

技术领域

本发明涉及相机标定技术领域，具体地，涉及一种移动终端相机标定方法。

背景技术

数字摄影测量作为摄影测量技术发展第三个阶段，以其原始数据与产品成果数字化、处理过程自动化等特点，已广泛用于三维测图、数字城市、文物保护以及三维导航等多个方面。

在数字摄影测量中，影像的质量是决定后期数据处理与最终成果精度与准确性的重要因素之一。然而，由于镜头制作工艺、原材料质量等因素的限制，影像不可避免地包含诸如直线弯曲、矩形错切等不可忽略的几何畸变。这些畸变将导致后期处理中产生诸如精度损失、影像变形等问题。因此，数字摄影测量在获取影像前，往往需确定相机畸变系数等标定参数，即进行相机标定，并利用标定结果对影像进行重采样，消除影像的畸变。

随着电子技术的发展、手机、平板等移动终端逐渐普及并广泛应用于摄影摄像中，其拍摄影像也同样存在畸变的问题。

目前，国内外许多学者针对不同类型相机、不同标定方法均进行了大量的研究，但目前还没有一种适用于所有相机的标定算法和流程。传统的相机标定是在获取标定物影像的基础上，利用如CDW(Close-Range Digital Workstation)等专业软件，进行相机内方位元素与畸变系数等标定参数的解算。此类方法一般需要专门设备、指定程序以及特殊标定环境的支持，不但标定设备价格昂贵，也不能及时获取标定参数，缺乏一定的灵活性，不适用于移动终端。

发明内容

本发明的目的就在于克服上述现有技术的缺点和不足，提供一种移动终端相机标定方法，以准确、高效、灵活、低成本地消除移动终端拍摄影像的几何畸变，提高影像质量，解决传统标定法标定设备昂贵、操作不便等问题。

本发明解决上述问题所采用的技术方案是：

移动终端相机标定方法，包括以下步骤：

S1选用标定物：采用黑白棋盘格图像作为标定物，在标定物上选择多个特征点；

S2拍摄待解算影像：利用待标定的移动终端相机拍摄至少2张标定物影像；

S3识别特征点坐标：识别特征点在拍摄的标定物影像上的像方坐标，并确定其对应的物方坐标；

S4解算标定参数：根据特征点的像方坐标(u,v)和物方坐标(X,Y,Z)解算标定参数，所述标定参数包括主距分量(f_x,f_y)、偏心畸变系数(c_x,c_y)、径向畸变系数(k₁,k₂,k₃)、切向畸变系数(p₁,p₂)；特征点的像方坐标(u,v)、物方坐标(X,Y,Z)的投影变化关系满足以下两式：

$s\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]={\mathrm{\λ}}^{\′}\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{R}_1& R_2& R_3& T\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]---\left(1\right);$

$\left[\begin{array}{c}u\\ v\end{array}\right]=(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)\left[\begin{array}{c}{u}^{\′}\\ v^{\′}\end{array}\right]+\left[\begin{array}{c}{2p}_1{u}^{\′}{v}^{\′}+p_2(r^2+{2u}^{\′2})\\ {2p}_2{u}^{\′}{v}^{\′}+p_1(r^2+{2v}^{\′2})\end{array}\right],r^2=u^{\′2}+v^{\′2}---\left(2\right).$

公式(1)中，R₁、R₂、R₃为三个旋转向量，用于表示空间点物方坐标与像方坐标的旋转变化；T为平移向量，用于表示空间点物方坐标与像方坐标的平移变化；s为齐次坐标系数；λ'为图像变化时的缩放系数，用于表示经旋转、平移以及内方位元素投影变换后的空间点物方坐标与像方坐标的缩放变化；

公式(2)中,(u',v')是将物方坐标(X,Y,Z)代入公式(1)计算得到的像方坐标理论值。

作为本发明的进一步改进，步骤S4中根据特征点的像方坐标(u,v)和物方坐标(X,Y,Z)解算标定参数的过程包括以下步骤：

S41、将特征点的像方坐标(u,v)及其对应的物方坐标(X,Y,Z)代入公式(1)求解主距分量(f_x,f_y)、偏心畸变系数(c_x,c_y)，本步骤中至少代入4个特征点的像方坐标(u,v)和物方坐标(X,Y,Z)，这些特征点至少来自两张影像，每张影像至少选取两个特征点；

S42、将解出的f_x、f_y、c_x、c_y的值代入公式(1)，并将特征点的物方坐标(X,Y,Z)代入公式(1)中计算出的特征点的像方坐标理论值，该像方坐标理论值用(u′,v′)表示；

S43、将特征点的像方坐标理论值(u′,v′)以及像方坐标(u,v)代入公式(2)，采用最小二乘法求解径向畸变系数(k₁,k₂,k₃)、切向畸变系数(p₁,p₂)。

进一步，步骤S41包括以下步骤：

S411、首先利用标定物黑白棋盘格图像为平面，即所有特征点物方坐标Z＝0的特点，将公式(1)简化为：

$s\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]={\mathrm{\λ}}^{\′}A\left[\begin{array}{ccc}{R}_1& R_2& T\end{array}\right]\left[\begin{array}{c}X\\ Y\\ 1\end{array}\right],A=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]---\left(3\right)$

S412、建立单应性矩阵H，令单应性矩阵H为：

H＝[h₁h₂h₃]＝λ′A[R₁R₂T]   (4)

S413、将特征点的像方坐标(u,v)及其对应的物方坐标(X,Y)带入式(3)，在最小二乘原理下，采用平面坐标转换参数确定方法求得λ′A[R₁R₂T]即H的具体值，从而可求得h₁、h₂和h₃的具体值；

S414、由于R₁,R₂为旋转向量，利用其相互正交且模相等两个特点，构建约束方程：

$\left\{\begin{array}{c}{h}_1^T{\mathrm{\λ}}^2{A}^{-T}{A}^{-1}{h}_2=0\\ h_1^T{\mathrm{\λ}}^2{A}^{-T}{A}^{-1}{h}_1=h_2^T{\mathrm{\λ}}^2{A}^{-T}{A}^{-1}{h}_2\end{array}\right.---\left(5\right)$

其中， $\mathrm{\λ}=\frac{1}{{\mathrm{\λ}}^{\′}};$

S415、将步骤S413中求得的h₁、h₂的具体值代入公式(5)中求出λ²A^-TA^-1的值；

S416、求解f_x、f_y、c_x、c_y：

由于 $A=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right],$ 则 ${\mathrm{\λ}}^2{A}^{-T}{A}^{-1}={\mathrm{\λ}}^2\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& \frac{{-c}_x}{f_x^2}\\ 0& \frac{1}{f_y^2}& \frac{{-c}_y}{f_y^2}\\ \frac{-c_x}{f_x^2}& \frac{{-c}_y}{f_y^2}& \frac{c_x}{f_x^2}+\frac{c_y}{f_y^2}+1\end{array}\right]$

假定步骤S415中求得的λ²A^-TA^-1的值为 $\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{12}& B_{22}& B_{23}\\ B_{13}& B_{23}& B_{33}\end{array}\right],$ 等量代换得

${\mathrm{\λ}}^2{A}^{-T}{A}^{-1}=\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& \frac{-c_x}{f_x^2}\\ 0& \frac{1}{f_y^2}& \frac{{-c}_y}{f_y^2}\\ \frac{{-c}_x}{f_x^2}& \frac{{-c}_y}{f_y^2}& \frac{c_x}{f_x^2}+\frac{c_y}{f_y^2}+1\end{array}\right]=\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{12}& B_{22}& B_{23}\\ B_{13}& B_{23}& B_{33}\end{array}\right]---\left(6\right)$

根据公式(6)可求得：

$\left\{\begin{array}{c}{f}_x=\sqrt{\mathrm{\λ}/B_{11}}\\ f_y=\sqrt{\mathrm{\λ}/(B_{11}{B}_{22}-B_{12}^2)}\\ c_x=-B_{12}{f}_x^2/\mathrm{\λ}\\ c_y=(B_{12}{B}_{13}-B_{11}{B}_{23})/(B_{11}{B}_{22}-B_{12}^2)\\ \mathrm{\λ}=B_{33}-[B_{13}^2+c_y(B_{12}{B}_{13}-B_{11}{B}_{23})]/B_{11}\end{array}\right.---\left(7\right).$

进一步，步骤S41中还利用剩余的特征点的像方坐标(u,v)和物方坐标(X,Y,Z)采用最小二乘法对主距分量(f_x,f_y)、偏心畸变系数(c_x,c_y)进行优化，优化过程中特征点满足的条件为：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{m}_{\mathrm{ij}}-\hat{m}{(A,{\stackrel{\‾}{R}}_i,{\stackrel{\‾}{T}}_i,M_{\mathrm{ij}})}^2=\min ---\left(8\right),$

式中，n为待解算影像张数，m为每张待解算影像包含的特征点个数，m_ij为第i张待解算影像中第j个角点的像方坐标(u,v)，为第i张待解算影像的旋转矩阵；M_ij为第i张待解算影像中第j个角点的物方坐标(X,Y,Z)，为第i张待解算影像的平移向量；为通过公式(1)计算出的第i张待解算影像中第j个角点的像方坐标理论值(u′,v′)。

进一步，步骤S43中采用最小二乘法求解径向畸变系数(k₁,k₂,k₃)、切向畸变系数(p₁,p₂)时满足的条件为：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{m}_{\mathrm{ij}}-\hat{m}{(A,k_1,k_2,k_3,p_1,p_2,{\stackrel{\‾}{R}}_i,{\stackrel{\‾}{T}}_i,M_{\mathrm{ij}})}^2=\min ---\left(9\right)$

式中，n为待解算影像张数，m为每张待解算影像包含的特征点个数，m_ij为第i张待解算影像中第j个角点的像方坐标(u,v)，为第i张待解算影像的旋转矩阵；M_ij为第i张待解算影像中第j个角点的物方坐标(X,Y,Z)，为第i张待解算影像的平移向量；为通过公式(1)和公式(2)计算出的第i张待解算影像中第j个角点的最终理论像方坐标。

进一步，步骤S1中，将所述黑白棋盘格图像的每四个黑白交错的方格的中心点作为特征点，选择特征点分布为6×9或7×10的黑白棋盘格图像作为标定物。

进一步，步骤S2中，利用待标定的移动终端相机拍摄5张标定物影像。

进一步，步骤S2中，利用待标定的移动终端相机从不同的角度拍摄5张标定物影像，拍摄5张标定物影像时的移动终端相机与标定物的垂直距离相等，移动终端到标定物中心的直线与标定物成40-50度夹角。

进一步，步骤S2中，利用待标定的移动终端相机从不同的角度拍摄5张标定物影像，1张正对拍摄，另外四张分别从上下左右四个方向以移动终端到标定物中心的直线与标定物成40-50度夹角进行拍摄。

进一步，步骤S3包括以下步骤：

S31特征点探测：获取每张待解算影像的各特征点的像素级的影像坐标(u⁰,v⁰)，判定特征点的个数和排布情况是否与标定物一致，若一致则视拍摄影像正确并予以保留，若不一致，则将特征点存在错误的影像剔除；

S32精化坐标：采用亚像素检测方法获取特征点的亚像素级影像坐标(u″,v″)，提供后期标定参数解算用的像方坐标数据；

S33建立特征点对：确定影像中特征点对应的物方坐标，并关联每个特征点的物方坐标和像方坐标。

综上，本发明的有益效果是：

1、本发明是专门针对移动终端相机物理特性设计的移动终端相机标定方法，实现移动终端相机标定参数解算。解决传统标定法工作时间长、操作不便、价格昂贵、不适用于移动终端等问题，为有效实现移动终端相机标定，获取高质量影像提供了有力保证。

2、本发明使用6×9、7×10的黑白棋盘格图像制作标定板，能够满足成像清晰、完整的条件，而且特征点数合理，能有效提高解算结果精度。

3、本发明拍摄待解算图像时，从左、右、上、下45°拍摄，待解算图像相互之间的区别程度大，单应性矩阵不相近，包含了更为丰富的解算信息，使得其在校正结果精度方面，相比于随机拍摄有更高的精度与更好的稳定性，重采样影像效果也更好。

附图说明

图1是移动终端相机标定方法的一种实施流程图；

图2是本发明的实施例中使用的标定物图像。

具体实施方式

下面结合实施例及附图，对本发明作进一步地的详细说明，但本发明的实施方式不限于此。

移动终端相机标定方法，其标定流程如图1所示，包括以下步骤：

(1)选用标定物、设计标定板

利用Photoshop绘图软件，绘制方格大小相同黑白相间的“国际象棋棋盘”图像，作为待拍摄的标准图像即标定物，如图2所示。每四个黑白交错的方格中间有个重合的点，是这4个方格的中心点，称为角点，也作为后期标定参数解算时使用的特征点。在标定物影像中角点分布稀疏则标定参数解算数据不足，过密则影响后期角点识别精度。考虑到移动终端相机镜头尺寸小的特点，为保证拍摄影待解算影像的质量，标定板图像中角点的分布情况定位为6×9或7×10，6×9表示6排9列，共有6排角点，每排9个角点共54个角点，即54个特征点；同样7×10表示7排10列共70个特征点。实际应用中，也可以采用现成的满足前述要求的黑白棋盘格图像作为标定物，不必自行绘制。

本实施例中，利用电脑显示屏或平板电脑等成像设备无畸变地展示黑白棋盘格图像，将电脑显示屏或平板电脑屏幕作为标定板。实际应用中，也不仅限于使用前述电脑显示屏或平板电脑等成像设展示黑白棋盘格图像，也可以将黑白棋盘格图像制成纸板、木板等板状结构进行展示。

(2)拍摄待解算影像

使用待标定的移动终端在固定对焦模式下，从展示黑白棋盘格图像的成像设备屏幕的正面、左侧面、右侧面、上侧面以及下侧面对标定板影像进行拍摄，获取当前视场下共计5张待解算影像，作为后期标定参数解算的起算数据。拍摄过程中，移动终端与展示黑白棋盘格图像的成像设备屏幕的距离大致相等，且为减少重复角点坐标数据对解算精度的影响，4张侧面影像在拍摄时，移动终端与成像设备屏幕近似成45度。

(3)解算标定参数

需解算的标定参数包括主距分量(f_x,f_y)、偏心畸变系数(c_x,c_y)两个内方位元素以及径向畸变系数(k₁,k₂,k₃)、切向畸变系数(p₁,p₂)两个畸变系数。考虑上述标定参数，被摄物体的特征点在真实空间与像平面空间变化时，其物方坐标(X,Y,Z)与像方坐标(u,v)满足的投影变化关系如式(1)、式(2)所示：

$s\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]={\mathrm{\λ}}^{\′}\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{R}_1& R_2& R_3& T\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]---\left(1\right)$

$\left[\begin{array}{c}u\\ v\end{array}\right]=(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)\left[\begin{array}{c}{u}^{\′}\\ v^{\′}\end{array}\right]+\left[\begin{array}{c}{2p}_1{u}^{\′}{v}^{\′}+p_2(r^2+{2u}^{\′2})\\ {2p}_2{u}^{\′}{v}^{\′}+p_1(r^2+{2v}^{\′2})\end{array}\right],r^2=u^{\′2}+v^{\′2}---\left(2\right)$

式(1)式(2)即为包含标定参数的共线方程，其中(u,v)为未消除几何畸变的特征点像方坐标，R₁、R₂、R₃分别为三个旋转向量，用于表示空间点物方坐标与像方坐标的旋转变化，前述空间点即物方空间的点；T为平移向量，用于表示空间点物方坐标与像方坐标的平移变化，s为齐次坐标系数，λ'为图像变化时的缩放系数，用于表示经旋转、平移以及内方位元素投影变换后空间点物方坐标与像方坐标的缩放变化。(u',v')是将物方坐标(X,Y,Z)代入公式(1)计算得到的像方坐标理论值。标定参数解算的核心就是将特征点的像方坐标(u,v)与物方坐标(X,Y,Z)作为已知量代入包含相机畸变系数的共线方程，完成对未知标定参数的求解。本发明采用张正友提出的解算方法，实现参数数值的解算。该步骤中具体的执行过程包括以下步骤：

1)角点探测

选取解算影像(至少包括从正面、左侧面、右侧面、上侧面以及下侧面拍摄的解算影像)，并通过识别角点，判断待解算影像的正确性。首先获取每张待解算影像的各角点的像素级(单位：像素)像方坐标(u⁰,v⁰)，再判定角点的个数与排布情况是否与标定物图像一致，若一致则视拍摄影像正确并予以保留，若不一致，则将于高光、抖动等客观拍摄因素所导致角点存在错误的影像，从起算数据中剔除，从新选取待解算图像，若拍摄影像数量不够，重新拍摄影像用于解算标定参数。正确识别角点的个数与分布情况是后期正确建立约束方程的保证。当影像出现高亮等异常时，会引起角点识别错误，进而影响校正解算正常执行。故对角点的个数与分布情况进行探测，剔除无效影像，可以避免后期错误解算，保证标定参数的准确性。

2)精化坐标

即优化角点像素级图像坐标。该步骤中主要利用计算机视觉领域中的亚像素检测方法，对待解算图像中的棋盘角点坐标进行迭代搜索、循环优化以达到获取角点的亚像素级(单位：1/10像素)像方坐标(u″,v″)的目的，为后期标定参数的解算提供高精度的图像坐标数据。后期带入公式(1)、(2)中解算的特征点的坐标(u,v)全部采用该步骤中精化后的亚像素级像方坐标(u″,v″)的值。

3)建立特征点对

即确定待解算影像中角点对应的物方坐标，并与之对应的图像坐标相关联。因标定物黑白棋盘格图像为平面，角点的物方坐标中(X,Y)参照标定物图像的黑白格行列号定义，Z值则统一为0。对于6×9的黑白棋盘格图像标定物，其左下角的坐标为(0,0,0)，相邻特征点的坐标相差1，正右方紧邻特征点(0,0,0)的角点的坐标为(0,1,0)，正上方紧邻特征点(0,0,0)的角点的坐标为(1,0,0)，以此类推，右上角的坐标为(6,9,0)。

4)参数解算

即标定参数的数值解算。本实施例中，以6×9的黑白棋盘格图像标定物为例说明标定参数的数值解算过程。

首先利用标定物黑白棋盘格图像为平面，所有角点物方坐标的Z＝0的条件将共线方程简化为：

$s\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]={\mathrm{\λ}}^{\′}A\left[\begin{array}{ccc}{R}_1& R_2& T\end{array}\right]\left[\begin{array}{c}X\\ Y\\ 1\end{array}\right],A=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]---\left(3\right)$

建立单应性矩阵H，并令单应性矩阵H为：

H＝[h₁h₂h₃]＝λ′A[R₁R₂T]    (4)

利用角点的精化后像方坐标(u″,v″)与物方坐标(X,Y)在最小二乘原理下，采用平面坐标转换参数确定(四参数方程参数确定)方法求得λ′A[R₁R₂T]即H的具体值。假定H的具体值为 $\left[\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{12}& h_{22}& h_{23}\\ h_{13}& h_{23}& h_{33}\end{array}\right],$ 则h₁、h₂、h₃也求得具体值，那么 $h_1=\left[\begin{array}{c}{h}_{11}\\ h_{12}\\ h_{13}\end{array}\right]={\mathrm{\λ}}^{\′}{\mathrm{AR}}_1,$ $h_2=\left[\begin{array}{c}{h}_{12}\\ h_{22}\\ h_{32}\end{array}\right]={\mathrm{\λ}}^{\′}{\mathrm{AR}}_2,$ 令 $\mathrm{\λ}=\frac{1}{{\mathrm{\λ}}^{\′}},$ 则R₁＝λA^-1h₁、R₂＝λA^-1h₂。

由于R₁,R₂为旋转向量，利用其相互正交(即)且模相等(即)两个特点，构建约束方程：

$\left\{\begin{array}{c}{h}_1^T{\mathrm{\λ}}^2{A}^{-T}{A}^{-1}{h}_2=0\\ h_1^T{\mathrm{\λ}}^2{A}^{-T}{A}^{-1}{h}_1=h_2^T{\mathrm{\λ}}^2{A}^{-T}{A}^{-1}{h}_2\end{array}\right.---\left(5\right)$

将上述求得的h₁、h₂的具体值代入公式(5)中可求出λ²A^-TA^-1的值；假定求得的λ²A^-TA^-1的值为 $\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{12}& B_{22}& B_{23}\\ B_{13}& B_{23}& B_{33}\end{array}\right];$

又由于 $A=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right],$ 则 $A^{-T}{A}^{-1}=\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& \frac{{-c}_x}{f_x^2}\\ 0& \frac{1}{f_y^2}& \frac{{-c}_y}{f_y^2}\\ \frac{-c_x}{f_x^2}& \frac{{-c}_y}{f_y^2}& \frac{c_x}{f_x^2}+\frac{c_y}{f_y^2}+1\end{array}\right],$

从而

${\mathrm{\λ}}^2{A}^{-T}{A}^{-1}={\mathrm{\λ}}^2\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& \frac{{-c}_x}{f_x^2}\\ 0& \frac{1}{f_y^2}& \frac{{-c}_y}{f_y^2}\\ \frac{-c_x}{f_x^2}& \frac{{-c}_y}{f_y^2}& \frac{c_x}{f_x^2}+\frac{c_y}{f_y^2}+1\end{array}\right];$

等量代换得：

${\mathrm{\λ}}^2{A}^{-T}{A}^{-1}=\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& \frac{-c_x}{f_x^2}\\ 0& \frac{1}{f_y^2}& \frac{{-c}_y}{f_y^2}\\ \frac{{-c}_x}{f_x^2}& \frac{{-c}_y}{f_y^2}& \frac{c_x}{f_x^2}+\frac{c_y}{f_y^2}+1\end{array}\right]=\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{12}& B_{22}& B_{23}\\ B_{13}& B_{23}& B_{33}\end{array}\right]---\left(6\right)$

方程(6)的解即为所需的内方位元素，由式(6)可知A所包含未知数包括(f_x,f_y,c_x,c_y,λ)在内一共5个，但比例系数λ与其他未知数线性相关，因而实际需要求解的未知参数个数为4，而每张影像的单应性矩阵唯一且只能提供2个方程。为求解该方程，至少需要2张待解算影像。当待解算影像张数为2，即方程适定时，该方程有唯一线性解。

求解方程(6)得f_x、f_y、c_x、c_y的解为

$\left\{\begin{array}{c}{f}_x=\sqrt{\mathrm{\λ}/B_{11}}\\ f_y=\sqrt{\mathrm{\λ}/(B_{11}{B}_{22}-B_{12}^2)}\\ c_x=-B_{12}{f}_x^2/\mathrm{\λ}\\ c_y=(B_{12}{B}_{13}-B_{11}{B}_{23})/(B_{11}{B}_{22}-B_{12}^2)\\ \mathrm{\λ}=B_{33}-[B_{13}^2+c_y(B_{12}{B}_{13}-B_{11}{B}_{23})]/B_{11}\end{array}\right.---\left(7\right)$

由于拍摄的待解算影像张数为5，存在多余特征点坐标数据，故将这些多余坐标数据依次带入方程(1)，在最小二乘准则下，对内方位元素的解算值进行进一步优化。针对6×9的标定板图像，需要满足的条件为：

$\underset{i=1}{\overset{5}{\mathrm{\Σ}}}\underset{j=1}{\overset{54}{\mathrm{\Σ}}}{m}_{\mathrm{ij}}-\hat{m}{(A,{\stackrel{\‾}{R}}_i,{\stackrel{\‾}{T}}_i,M_{\mathrm{ij}})}^2=\min$

上式中，

m_ij为第i张待解算影像中第j个角点的像方坐标(u,v)；计算时其值是由计算机识

别并完成精化的角点像方坐标(u″,v″)；

为第i张待解算影像的旋转矩阵，即R_i1、R_i2、R_i3为第i张待解算影像的旋转向量；

M_ij为第i张待解算影像中第j个角点的物方坐标(X,Y,Z)；

为第i张待解算影像的平移向量；

为通过公式(1)计算出的第i张待解算影像中第j个角点的像方坐标理论值(u′,v′)。

最后，将求解的内方位元素f_x、f_y、c_x、c_y作为已知量代入公式(1)和公式(2)，同样在最小二乘准则下求解其余的畸变系数，此时，需要满足的条件为：

$\underset{i=1}{\overset{5}{\mathrm{\Σ}}}\underset{j=1}{\overset{54}{\mathrm{\Σ}}}{m}_{\mathrm{ij}}-\hat{m}{(A,k_1,k_2,k_3,p_1,p_2,{\stackrel{\‾}{R}}_i,{\stackrel{\‾}{T}}_i,M_{\mathrm{ij}})}^2=\min$

式中，m_ij为第i张待解算影像中第j个角点的像方坐标(u,v)，计算时其值是由计

算机识别并完成精化的角点像方坐标(u″,v″)；

为第i张待解算影像的旋转矩阵；

M_ij为第i张待解算影像中第j个角点的物方坐标(X,Y,Z)；

为第i张待解算影像的平移向量；

为通过公式(1)和公式(2)计算出的第i张待解算影像中第j个角点的最终理论像方坐标，计算过程中物方坐标首先带入公式(1)计算出像方坐标理论值(u'，v')，再将(u'，v')带入公式(2)计算最终理论像方坐标。

经在最小二乘准则下的反复优化，完成所有特征点的带入求解后，至此完成所有标定参数的解算。

(3)成果输出

完成上述标定处理后，输出当前视场的内方位元素(f_x,f_y)、(c_x,c_y)与畸变系数(k₁,k₂,k₃)、(p₁,p₂)的解算值。

本实施例中，还对将特征点为5×8、6×9、7×10、8×11、9×11的黑白棋盘格图像制作的标定板进行了对比测试，分别测试计算了前述标定板在相同成像分辨率以及成像尺寸的条件下的解算结果精度和稳定性，计算结果表明按前述秩序，解算结果精度先下降，后趋于平稳。对于5×8标定板而言，由于特征点数较少，虽满足解算条件，但用于迭代解算的数据少，因而在解算时精度不高。6×9与7×10的标定板相比，两者在解算时的精度大致相近，但在解算的稳定性方面，7×10标定板高于6×9标定板，8×11与9×11标定板则在解算结果精度上，有降低并渐于稳定的趋势。这与特征点数过多有关：过多的特征点分布于同面积的屏幕时，棋盘格的尺寸缩小，点间距离也相应变短，计算机特征点识别难度增加，降低了角点坐标的识别精度。另一方面是因为标定参数解算是一个迭代过程，当特征点数达到一定时结算过程收敛，解算结果趋于稳定。因此，在成像屏幕大小一定的情况下，若设计标定板中包含过少的特征点，则不能够提供足够多的解算信息，难以保证解算结果精度；反之，过多特征点则会增加特征点的识别难度，降低特征点的坐标精度，同样将会引起解算结果精度的降低。因此，本实施例中使用的6×9、7×10的黑白棋盘格图像制作的标定板在满足成像清晰、完整的条件下，特征点数合理，能有效提高解算结果精度。

此外，本实施例中拍摄待解算图像时，从展示黑白棋盘格图像的成像设备屏幕的正面、左侧面、右侧面、上侧面以及下侧面对标定板影像进行拍摄，拍摄过程中，移动终端与展示黑白棋盘格图像的成像设备屏幕的距离大致相等，移动终端与成像设备屏幕近似成45度。从左右上下45°拍摄的影像，相互之间区别程度大，单应性矩阵不相近，包含了更为丰富的解算信息，使得其在校正结果精度方面，相比于随机拍摄有更高的精度与更好的稳定性，重采样影像效果也更好。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (10)

1.移动终端相机标定方法，其特征在于，包括以下步骤：

S1选用标定物：采用黑白棋盘格图像作为标定物，在标定物上选择多个特征点；

S2拍摄待解算影像：利用待标定的移动终端相机拍摄至少2张标定物影像；

S3识别特征点坐标：识别特征点在拍摄的标定物影像上的像方坐标，并确定其对应的物方坐标；

S4解算标定参数：根据特征点的像方坐标(u,v)和物方坐标(X,Y,Z)解算标定参数，所述标定参数包括主距分量(f_x,f_y)、偏心畸变系数(c_x,c_y)、径向畸变系数(k₁,k₂,k₃)、切向畸变系数(p₁,p₂)；特征点的像方坐标(u,v)、物方坐标(X,Y,Z)的投影变化关系满足以下两式：

$s\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]={\mathrm{\λ}}^{\′}\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{R}_1& R_2& R_3& T\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]---\left(1\right);$

$\left[\begin{array}{c}u\\ v\end{array}\right]=(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)\left[\begin{array}{c}{u}^{\′}\\ v^{\′}\end{array}\right]+\left[\begin{array}{c}2p_1{u}^{\′}{v}^{\′}+p_2(r^2+{2u}^{\′2})\\ {2p}_2{u}^{\′}{v}^{\′}+p_1(r^2+{2v}^{\′2})\end{array}\right],r^2=u^{\′2}+v^{\′2}---\left(2\right).$

公式(1)中，R₁、R₂、R₃为三个旋转向量，用于表示空间点物方坐标与像方坐标的旋转变化；T为平移向量，用于表示空间点物方坐标与像方坐标的平移变化；s为齐次坐标系数；λ'为图像变化时的缩放系数，用于表示经旋转、平移以及内方位元素投影变换后的空间点物方坐与像方坐标的缩放变化；

公式(2)中，(u',v')是将物方坐标(X,Y,Z)代入公式(1)计算得到的像方坐标理论值。

2.根据权利要求1所述的移动终端相机标定方法，其特征在于，步骤S4中根据特征点的像方坐标(u,v)和物方坐标(X,Y,Z)解算标定参数的过程包括以下步骤：

S41、将特征点的像方坐标(u,v)及其对应的物方坐标(X,Y,Z)代入公式(1)求解主距分量(f_x,f_y)、偏心畸变系数(c_x,c_y)，本步骤中至少代入4个特征点的像方坐标(u,v)和物方坐标(X,Y,Z)，这些特征点至少来自两张影像，每张影像至少选取两个特征点；

S42、将解出的f_x、f_y、c_x、c_y的值代入公式(1)，并将特征点的物方坐标(X,Y,Z)代入公式(1)中计算出的特征点的像方坐标理论值，该像方坐标理论值用(u′,v′)表示；

S43、将特征点的像方坐标理论值(u′,v′)以及像方坐标(u,v)代入公式(2)，采用最小二乘法求解径向畸变系数(k₁,k₂,k₃)、切向畸变系数(p₁,p₂)。

3.根据权利要求2所述的移动终端相机标定方法，其特征在于，步骤S41包括以下步骤：

S411、首先利用标定物黑白棋盘格图像为平面，即所有特征点物方坐标Z＝0的特点，将公式(1)简化为：

$s\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]={\mathrm{\λ}}^{\′}A\left[\begin{array}{ccc}{R}_1& R_2& T\end{array}\right]\left[\begin{array}{c}X\\ Y\\ 1\end{array}\right],A=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]$

S412、建立单应性矩阵H，令单应性矩阵H为：

H＝[h ₁h₂ h₃]＝λ′A[R₁ R₂ T]     (4)

S413、将特征点的像方坐标(u,v)及其对应的物方坐标(X,Y)带入式(3)，在最小二乘原理下，采用平面坐标转换参数确定方法求得λ′A[R₁ R₂ T]即H的具体值，从而可求得h₁、h₂和h₃的具体值；

S414、由于R₁,R₂为旋转向量，利用其相互正交且模长相等两个特点，构建约束方程：

$\left\{\begin{array}{c}{h}_1^T{\mathrm{\λ}}^2{A}^{-T}{A}^{-1}{h}_2=0\\ h_1^T{\mathrm{\λ}}^2{A}^{-T}{A}^{-1}{h}_1=h_2^T{\mathrm{\λ}}^2{A}^{-T}{A}^{-1}{h}_2\end{array}\right.---\left(5\right)$

其中， $\mathrm{\λ}=\frac{1}{{\mathrm{\λ}}^{\′}};$

S415、将步骤S413中求得的h₁、h₂的具体值代入公式(5)中求出λ²A^-TA^-1的值；

S416、求解f_x、f_y、c_x、c_y：

由于 $A=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right],$ 则 ${\mathrm{\λ}}^2{A}^{-T}{A}^{-1}={\mathrm{\λ}}^2\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& \frac{-c_x}{f_x^2}\\ 0& \frac{1}{f_y^2}& \frac{-c_y}{f_y^2}\\ \frac{{-c}_x}{f_x^2}& \frac{{-c}_y}{f_y^2}& \frac{c_x}{f_x^2}+\frac{c_y}{f_y^2}+1\end{array}\right]$

假定步骤S415中求得的λ²A^-TA^-1的值为 $\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{12}& B_{22}& B_{23}\\ B_{13}& B_{23}& B_{33}\end{array}\right],$ 等量代换得

${\mathrm{\λ}}^2{A}^{-T}{A}^{-1}={\mathrm{\λ}}^2\left[\begin{array}{ccc}\frac{1}{f_x^2}& 0& \frac{-c_x}{f_x^2}\\ 0& \frac{1}{f_y^2}& \frac{-c_y}{f_y^2}\\ \frac{{-c}_x}{f_x^2}& \frac{{-c}_y}{f_y^2}& \frac{c_x}{f_x^2}+\frac{c_y}{f_y^2}+1\end{array}\right]=\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{12}& B_{22}& B_{23}\\ B_{13}& B_{23}& B_{33}\end{array}\right]---\left(6\right)$

根据公式(6)可求得：

$\left\{\begin{array}{c}{f}_x=\sqrt{\mathrm{\λ}/B_{11}}\\ f_y=\sqrt{\mathrm{\λ}/(B_{11}{B}_{22}-B_{12}^2)}\\ c_x=-B_{12}{f}_x^2/\mathrm{\λ}\\ c_y=(B_{12}{B}_{13}-B_{11}{B}_{23})/(B_{11}{B}_{22}-B_{12}^2)\\ \mathrm{\λ}=B_{33}-[B_{13}^2+c_y(B_{12}{B}_{13}-B_{11}{B}_{23})]/B_{11}\end{array}\right.---\left(7\right).$

4.根据权利要求2所述的移动终端相机标定方法，其特征在于，步骤S41中还利用剩余的特征点的像方坐标(u,v)和物方坐标(X,Y,Z)采用最小二乘法对主距分量(f_x,f_y)、偏心畸变系数(c_x,c_y)进行优化，优化过程中特征点满足的条件为：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{m}_{\mathrm{ij}}-\hat{m}{(A,{\stackrel{\‾}{R}}_i,{\stackrel{\‾}{T}}_i,M_{\mathrm{ij}})}^2=\min ---\left(8\right),$

式中，n为待解算影像张数，m为每张待解算影像包含的特征点个数，m_ij为第i张待解算影像中第j个角点的像方坐标(u,v)，为第i张待解算影像的旋转矩阵；M_ij为第i张待解算影像中第j个角点的物方坐标(X，Y，Z)，为第i张待解算影像的平移向量；为通过公式(1)计算出的第i张待解算影像中第j个角点的像方坐标理论值(u′,v′)。

5.根据权利要求2所述的移动终端相机标定方法，其特征在于，步骤S43中采用最小二乘法求解径向畸变系数(k₁,k₂,k₃)、切向畸变系数(p₁,p₂)时满足的条件为：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{m}_{\mathrm{ij}}-\hat{m}{(A,k_1,k_2,k_3,p_1,p_2,{\stackrel{\‾}{R}}_i,{\stackrel{\‾}{T}}_i,M_{\mathrm{ij}})}^2=\min ---\left(9\right)$

式中，n为待解算影像张数，m为每张待解算影像包含的特征点个数，m_ij为第i张待解算影像中第j个角点的像方坐标(u,v)，为第i张待解算影像的旋转矩阵；M_ij为第i张待解算影像中第j个角点的物方坐标(X,Y,Z)，为第i张待解算影像的平移向量；为通过公式(1)和公式(2)计算出的第i张待解算影像中第j个角点的最终理论像方坐标。

6.根据权利要求1至5任一所述的移动终端相机标定方法，其特征在于，步骤S1中，将所述黑白棋盘格图像的每四个黑白交错的方格的中心点作为特征点，选择特征点分布为6×9或7×10的黑白棋盘格图像作为标定物。

7.根据权利要求1至5任一所述的移动终端相机标定方法，其特征在于，步骤S2中，利用待标定的移动终端相机拍摄5张标定物影像。

8.根据权利要求1至5任一所述的移动终端相机标定方法，其特征在于，步骤S2中，利用待标定的移动终端相机从不同的角度拍摄5张标定物影像，拍摄5张标定物影像时的移动终端相机与标定物的垂直距离相等，移动终端到标定物中心的直线与标定物成40-50度夹角。

9.根据权利要求8所述的移动终端相机标定方法，其特征在于，步骤S2中，利用待标定的移动终端相机从不同的角度拍摄5张标定物影像，1张正对拍摄，另外四张分别从上下左右四个方向以移动终端到标定物中心的直线与标定物成40-50度夹角进行拍摄。

10.根据权利要求1至5任一所述的移动终端相机标定方法，其特征在于，步骤S3包括以下步骤：

S31特征点探测：获取每张待解算影像的各特征点的像素级的影像坐标(u⁰,v⁰)，判定特征点的个数和排布情况是否与标定物一致，若一致则视拍摄影像正确并予以保留，若不一致，则将特征点存在错误的影像剔除；

S32精化坐标：采用亚像素检测方法获取特征点的亚像素级影像坐标(u″,v″)，提供后期标定参数解算用的像方坐标数据；

S33建立特征点对：确定影像中特征点对应的物方坐标，并关联每个特征点的物方坐标和像方坐标。