CN103150724B - 一种基于分段模型的摄像机标定方法 - Google Patents

Abstract

一种基于分段模型的摄像机标定方法属于摄像机标定方法领域，该标定方法计算标定模型参数的图像特征点分别在像素平面中划分的三个区域内，标定模型是分段函数；其用畸变较少的圆形区域内的特征点计算线性模型，求解并在后续标定模型中延用这些准确性高的内部参数矩阵A和外部参数R、t；仅需对圆环区域模型中的畸变系数求解，并对畸变系数的优化初始值进行求解，降低了待求解的维数，避免了计算结果的不收敛情况并大幅缩短了计算时间；其将第三区域所建立的畸变模型转化为连续的样条平滑函数以修正第三区域内的全部像素点畸变，使整幅图像中距图像中线较远位置的特殊畸变都能够被精确修正，大幅度地提高了标定精度。

Description

一种基于分段模型的摄像机标定方法

技术领域

本发明属于摄像机标定方法领域，具体涉及一种基于分段模型的摄像机标定方法。

背景技术

摄像机标定是由已知特征点的像素坐标和世界坐标去求解标定模型参数的过程，从而建立图像像素位置与三维场景点坐标之间的对应投影关系。

高精度摄像机标定是数字图像测量和三维重构技术的基础，摄像机两步标定法是目前传统标定方法中应用最为广泛的方法，具有较高的标定精度。当应用场合所要求的精度很高且摄像机的参数不经常变化时，传统标定方法中的两步标定法通常会作为首选方法。

两步标定法的第一步是先利用直接线性变换方法求解摄像机参数，第二步考虑畸变因素，将第一步中求得的参数作为初始值，利用优化算法重新对标定模型中的参数进行非线性求解，从而获得更准确的标定模型。

目前，两步标定法中有代表性的是Tsai，Heikkila和张正友的标定方法，这三种方法都采用了统一形式的畸变函数来描述整个镜头的畸变，其中Tsai和张正友的模型中考虑了径向畸变的影响，而Heikkila在此基础上又增加了切向畸变。

但是，由于镜头制造误差存在随机性，各种畸变误差在图像中的分布规律并不完全一致，因此，在实际应用中很难找到能够精准体现图像扭曲程度的畸变模型。如果想更准确和更详尽的描述镜头畸变，模型中的参数就会增加，复杂的模型反而会引起数值计算的不稳定性。

基于上述因素，现有的标定方法只能选择一些具体的畸变形式。也就是说，图像中会存在一些畸变类型，并未被选定的畸变模型所体现，最终会影响摄像机标定的精度。而且，现有的标定方法中待求解的维数往往较大，导致在优化过程中容易出现结果不收敛或计算时间过长的情况。

此外，现有方法中的畸变模型都试图利用一个统一的函数来描述整幅图像中的畸变，最终，标定模型中的参数是利用图像中所有特征点计算获得的。然而，位于距离图像中心较远的区域中的点通常带有不确定型畸变，这些点通常位于距离图像中心较远的区域中，在优化求解过程中会将大量误差带入模型中，最终降低和损害整幅图像的标定精度。

Tsai R Y在《A versatile camera calibration technique for high-accuracy 3Dmachine vision metrology using off-the-shelf TV cameras and lenses》[J].IEEEJournal of Robotics and Automation,1987,3(4):323-344中提出了一种摄像机标定法使用的转换坐标系，即：建立世界坐标系、摄像机坐标系、图像物理坐标系和图像像素平面坐标系。

张正友在Zhang Z.(2000).A flexible new technique for camera calibration.IEEE Trans.on Pattern Analysis and Machine Intelligence,22(11):1330-1334.提出了一种相机标定方法，其中包括一种具体的线性模型求解和标定方法。

发明内容

为了解决现有摄像机两步标定法标定模型中的参数都是利用图像中所有特征点计算获得并都是采用一个统一的函数来描述整幅图像中的畸变，图像中会存在一些畸变类型并未被选定的畸变模型所体现，以及现有的标定方法虽然大都能获得较好的优化初值，但是由于待求解的维数较大，在优化过程中的大量非线性计算容易出现计算结果不收敛或计算时间过长的情况的技术问题，本发明提供一种基于分段模型的摄像机标定方法。

本发明解决技术问题所采取的技术方案如下：

一种基于分段模型的摄像机标定方法，包括如下步骤：

步骤一：固定摄像机并将棋盘格标定板的中心置于摄像机光轴附近，使棋盘格标定板在一个相对固定的成像位置围绕摄像机镜头变换不同的倾角和姿态，用摄像机对标定板的每一种姿态拍摄一幅照片图像，最终取得多幅照片图像；

步骤二：按照Tsai R Y的坐标转换方法确立转换坐标系，即：建立世界坐标系、摄像机坐标系、图像物理坐标系和图像像素平面坐标系，世界坐标系内的三维坐标点坐标(X_w,Y_w,Z_w)被依次转换为摄像机点坐标(X_c,Y_c,Z_c)、理想图像点坐标(x_u,y_u)、真实图像点坐标(x_d,y_d)和像素点坐标(x_p,y_p)；另外，参照张正友的摄像机标定方法将世界坐标系的Z_w轴坐标设为0，即Z_w=0；此时，世界坐标系内任意点的坐标表示为(X_w,Y_w,0)，或仅表示为世界坐标系XY轴平面的二维点坐标(X_w,Y_w)；

步骤三：利用亚像素角点检测方法分别对步骤一所取得的各幅照片图像提取对应的角点像素坐标；

步骤四：将步骤二所述的图像像素坐标平面划分成三个相邻的区域：圆形区域、圆环形区域和第三区域；

步骤五：将步骤三所取得的各幅照片图像的角点像素坐标按照其在步骤四所属的三个不同区域分类，分别得到属于圆形区域内的图像角点像素检测坐标、属于圆环形区域内的图像角点像素检测坐标和属于第三区域内的图像角点像素检测坐标；

步骤六：按照步骤二所建立的转换坐标系在步骤四所取得的圆形区域建立线性标定模型，将棋盘格标定板上方格角点在世界坐标系对应的三维点坐标(X_w,Y_w,0)转换为属于圆形区域内的图像像素坐标(x_p,y_p)，忽略畸变的影响，所建立的圆形区域内对应平面的标定线性模型如下：

$s\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=A\left[\begin{array}{cccc}{r}_1& r_2& r_3& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 0\\ 1\end{array}\right]=A[R,t]\left[\begin{array}{c}{X}_w\\ Y_w\\ 0\\ 1\end{array}\right]...\left(1\right)$

当世界坐标系的Z_w轴坐标设为0时，即世界坐标系内任意点的坐标表示为(X_w,Y_w)时，式（1）可变换为下式

$s\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=A\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]=H\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]...\left(1.1\right)$

式（1）和式（1.1）中，s表示比例因子，R=[r1,r2,r3]是世界坐标的旋转矩阵，r_i（i=1,2,3）表示旋转矩阵R的第i列，t表示平移向量，R和t统称为标定模型的外部参数；

$A=\left[\begin{array}{ccc}\mathrm{\α}& \mathrm{\γ}& u_0\\ 0& \mathrm{\β}& v_0\\ 0& 0& 1\end{array}\right]$ 是内部参数矩阵，H表示3×3的单应性矩阵；α和β表示像素平面中U轴和V轴的比例因子，γ表示像素平面两坐标轴的不垂直因子，(u₀,v₀)表示相机光轴与图像平面的交点在像素坐标系上的坐标；

步骤七：利用步骤五所述圆形区域内的角点像素检测坐标(x_p,y_p)以及这些角点在步骤二所述世界坐标系中对应的二维点坐标值(X_w,Y_w)，并使用张正友的相机标定方法中线性模型求解方法，求解其对应在步骤六中所述的平面标定线性模型，完成对步骤四所述圆形区域的线性标定；在上述求解过程中，同时求得步骤六所述的内部参数A，并且同时求得各组步骤一所述每幅照片图像对应在步骤六所述的外部参数R和t；

步骤八：利用步骤七求出的内部参数A以及外部参数R和t，对步骤四所述属于圆环形区域建立非线性畸变模型如下：

$s\left[\begin{array}{c}{x}_u\\ y_u\\ 1\end{array}\right]=\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]...\left(2\right)$

$\left[\begin{array}{c}{x}_u\\ y_u\end{array}\right]=(1+k_1{r}^2+k_2{r}^4)\left[\begin{array}{c}{x}_d\\ y_d\end{array}\right]+\left[\begin{array}{c}2p_1{x}_d{y}_d+p_2(r^2+2{x_d}^2)\\ p_1(r^2+2{y_d}^2)+2p_2{x}_d{y}_d\end{array}\right]...\left(3\right)$

$\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=\left[\begin{array}{ccc}\mathrm{\α}& \mathrm{\γ}& u_0\\ 0& \mathrm{\β}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_d\\ y_d\\ 1\end{array}\right]...\left(4\right)$

式（3）中的 $r=\sqrt{{x_u}^2+{y_u}^2};$

k₁,k₂,p₁,p₂表示图像物理坐标系中径向与切向畸变函数的系数；

步骤九：求解步骤八中的畸变系数k₁,k₂,p₁,p₂，并根据已求得的畸变系数k₁,k₂,p₁,p₂代入步骤八所述的由式（2）、式（3）和式（4）构成的整体非线性畸变模型，以完成对步骤四所述圆环形区域的标定；

步骤十：获取步骤四所述属于第三区域内的图像的角点的理想像素坐标(u_d,v_d)；将步骤一所取得的每一幅照片图像中属于第三区域内的图像的角点在世界坐标系中对应的二维点坐标值(X_w,Y_w)作为转换目标，利用在步骤七获得的圆形区域内对应平面的标定线性模型的内部参数矩阵A和每幅照片图像对应的外部参数R和t对前述的二维点坐标值(X_w,Y_w)投影到像素坐标系平面，以获得它们对应的理想参考像素坐标(u_d,v_d)；

步骤十一：将步骤十所述理想参考像素坐标(u_d,v_d)与步骤五属于第三区域内的图像角点像素检测坐标(u_p,v_p)作差比较，获得属于第三区域内的图像角点像素检测坐标在像素平面中的离散偏差分布，如下式：

δ_u(u_p,v_p)=u_d-u_p

δ_v(u_p,v_p)=v_d-v_p……(9)

式（9）中δ_u和δ_v分别表示在角点像素坐标的像素平面中U轴方向和V轴方向的离散偏差分布；

步骤十二：拟合获得连续的样条平滑函数：利用matlab中的样条平滑算法对式（9）所述的两个离散偏差分布进行运算处理，获得步骤四所述属于像素平面中第三区域内图像像素点分别在U轴方向和V轴方向上的样条平滑修正函数：

δ_u=f_u(u,v)

δ_v=f_v(u,v)……(10)

式（10）中(u,v)表示第三区域内任意像素点的坐标，f_u表示像素点在像素平面U轴方向上的偏差修正函数，f_v表示像素点在像素平面V轴方向上的偏差修正函数；

步骤十三：在像素坐标平面内，利用步骤十二中取得的连续的样条平滑函数，对第三区域内的图像像素点的畸变进行修正；修正完成后，可获得步骤十所述理想参考像素坐标(u_d,v_d)，最终可利用步骤七获得的内部参数A以及每幅图像对应外部参数R和t将上述理想坐标点反投影成对应的世界坐标点，完成步骤四所述属于像素平面中第三区域的标定，进而完成摄像机的整个标定过程。

上述步骤四还包括如下步骤：

4-1）、确定圆形区域：设摄像机拍摄分辨率为M列×N行，将图像像素平面分辨率的中心点设定为圆心，并在像素平面上定义一个R₁为半径的圆形区域；

4-2）、确定圆环形区域：以R₁为内径，R₂为外径，在图像像素平面上定义一个圆环形区域，其圆心与圆形区域同心；

4-3）、确定第三区域：在图像像素平面上，将圆环形区域以外的图像区域定义为第三区域。

上述步骤九还包括如下步骤：

9-1）、通过使步骤四中的圆环形区域内角点像素坐标反投影的世界坐标值与其对应在真实世界坐标值的差值的平方和最小，建立四个畸变系数k₁,k₂,p₁,p₂的目标函数：

式（5）中M表示角点的世界坐标，表示角点像素坐标的反投影世界坐标，m表示标定图像的数目，n表示每幅标定图像中所用角点的数目；

9-2）、通过以下步骤求解四个畸变系数k₁,k₂,p₁,p₂的优化初值；

9-2-1）、利用线性最小二乘法将步骤八中的式（3）展开，将步骤八所述的畸变模型写成如下矩阵形式：

$\left[\begin{array}{cccc}{x}_d(x_d^2+y_d^2)& x_d{(x_d^2+y_d^2)}^2& 2x_d{y}_d& {(x_d^2+y_d^2)}^2+2x_d^2\\ y_d(x_d^2+y_d^2)& y_d{(x_d^2+y_d^2)}^2& {(x_d^2+y_d^2)}^2+2y_d^2& 2x_d{y}_d\end{array}\right]\left[\begin{array}{c}{k}_1\\ k_2\\ p_1\\ p_2\end{array}\right]=\left[\begin{array}{c}{x}_u-x_d\\ y_u-y_d\end{array}\right]...\left(6\right)$

上式中，四个畸变系数k₁,k₂,p₁,p₂是待求解变量；

9-2-2）、利用步骤七中获得的外部参数R和t以及步骤五所述圆环形区域内的角点像素检测坐标在其世界坐标系中对应的二维点坐标值(X_w,Y_w)，代入步骤八所述式（2）中，即可获得步骤五所述圆环形区域内的角点的理想图像坐标(x_u,y_u)；

9-2-3）、通过将步骤五所述属于圆环形区域内的图像角点像素检测坐标(x_p,y_p)和步骤七求得的内部参数A代入步骤八的式（4）中，即可获得步骤五所述属于圆环形区域内的图像角点的真实图像坐标(x_d,y_d)；

9-2-4）、设步骤一所述过程共取得m幅标定照片图像，每幅图像中包含n_i个步骤四所述属于圆环形区域内的图像角点，则根据式（6）可获得2m×n₀个方程，则这些方程的矩阵可写成如下表达形式：

DK=d……(7)

式（7）中D表示2mn₀×4的系数矩阵，d表示2mn₀维的向量，K=[k₁ k₂ p₁ p₂]^T,T为向量转置符号；

9-2-5）、利用线性最小二乘法将式（7）变形为

K=(D^TD)^-1D^Td……(8)

通过对式（8）求解可以解出K的初始值：

9-3）根据步骤9-2-5）所取得的K的最优初始值，利用matlab中的优化算法对步骤9-1）所述的目标函数F(k₁,k₂,p₁,p₂)进行最小寻优计算，得到使目标函数F(k₁,k₂,p₁,p₂)最小的k₁,k₂,p₁,p₂解。

本发明的有益效果是：该摄像机标定方法计算标定模型参数的图像特征点分别在像素平面中划分的三个区域内，标定模型是分段函数；该方法采用畸变较少的圆形区域内的特征点计算线性模型，求解并在后续标定模型中延用这些准确性高的内部参数矩阵A和外部参数R、t；仅需对圆环区域模型中的畸变系数求解，并对畸变系数的优化初始值进行求解，降低了待求解的维数，避免了计算结果的不收敛情况并大幅缩短了计算时间；该方法将第三区域所建立的畸变模型转化为连续的样条平滑函数以修正第三区域内的全部像素点畸变，使整幅图像中距图像中线较远位置的特殊畸变都能够被精确修正，大幅度地提高了标定精度。

附图说明

图1是本发明一种基于分段模型的摄像机标定方法的流程图；

图2是本发明中的转换坐标系示意图；

图3是本发明中对像素平面划分为圆形区域、圆环形区域和第三区域的示意图；

图4是本发明中第三区域内像素点在U轴方向上的偏差曲面示意图；

图5是本发明中第三区域内像素点在V轴方向上的偏差曲面示意图。

具体实施方式

下面结合附图对本发明做进一步详细说明。

如图1所示，本发明基于分段模型的摄像机标定方法包括如下步骤：

步骤一：固定摄像机并将棋盘格标定板的中心置于摄像机光轴附近，使棋盘格标定板在一个相对固定的成像位置围绕摄像机镜头变换不同的倾角和姿态，用摄像机对标定板的每一种姿态拍摄一幅照片图像，最终取得多幅照片图像；

步骤二：如图2所示，按照Tsai R Y的坐标转换方法确立转换坐标系，即：建立世界坐标系、摄像机坐标系、图像物理坐标系和图像像素平面坐标系，世界坐标系内的三维坐标点坐标(X_w,Y_w,Z_w)被依次转换为摄像机点坐标(X_c,Y_c,Z_c)、理想图像点坐标(x_u,y_u)、真实图像点坐标(x_d,y_d)和像素点坐标(x_p,y_p)。另外，参照张正友的摄像机标定方法将世界坐标系的Z_w轴坐标设为0，即Z_w=0。此时，世界坐标系内任意点的坐标表示为(X_w,Y_w,0)，或仅表示为世界坐标系XY轴平面的二维点坐标(X_w,Y_w)；

步骤三：利用亚像素角点检测方法分别对步骤一所取得的各幅照片图像提取对应的角点像素坐标；

步骤四：如图3所示，将步骤二所述的图像像素坐标平面划分成三个相邻的区域：圆形区域1、圆环形区域2和第三区域3，具体步骤如下：

4-1）、确定圆形区域：设摄像机拍摄分辨率为M列×N行，将图像像素平面分辨率的中心点设定为圆心，并在像素平面上定义一个R₁为半径的圆形区域1；

4-2）、确定圆环形区域：以R₁为内径，R₂为外径，在图像像素平面上定义一个圆环形区域2，其圆心与圆形区域1同心；

4-3）、确定第三区域：在图像像素平面上，将圆环形区域2以外的图像区域定义为第三区域3；

步骤五：将步骤三所取得的各幅照片图像的角点像素坐标按照其在步骤四所属的三个不同区域分类，分别得到属于圆形区域1内的图像角点像素检测坐标、属于圆环形区域2内的图像角点像素检测坐标和属于第三区域3内的图像角点像素检测坐标；

步骤六：按照步骤二所建立的转换坐标系在步骤四所取得的圆形区域1）建立线性标定模型，将棋盘格标定板上方格角点在世界坐标系对应的三维点坐标(X_w,Y_w,0)转换为属于圆形区域1内的图像像素坐标(x_p,y_p)，忽略畸变的影响，所建立的圆形区域1内对应平面的标定线性模型如下：

$s\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=A\left[\begin{array}{cccc}{r}_1& r_2& r_3& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 0\\ 1\end{array}\right]=A[R,t]\left[\begin{array}{c}{X}_w\\ Y_w\\ 0\\ 1\end{array}\right]...\left(1\right)$

当世界坐标系的Z_w轴坐标设为0时，即世界坐标系内任意点的坐标表示为(X_w,Y_w)时，式（1）可变换为下式：

$s\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=A\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]=H\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]...\left(1.1\right)$

式（1）和式（1.1）中，s表示比例因子，R=[r1,r2,r3]是世界坐标的旋转矩阵，r_i（i=1,2,3）表示旋转矩阵R的第i列，t表示平移向量，R和t统称为标定模型的外部参数；

$A=\left[\begin{array}{ccc}\mathrm{\α}& \mathrm{\γ}& u_0\\ 0& \mathrm{\β}& v_0\\ 0& 0& 1\end{array}\right]$ 是内部参数矩阵，H表示3×3的单应性矩阵，α和β表示像素平面中U轴和V轴的比例因子，γ表示像素平面两坐标轴的不垂直因子，(u₀,v₀)表示相机光轴与图像平面的交点在像素坐标系上的坐标；

步骤七：利用步骤五所述圆形区域1内的角点像素检测坐标(x_p,y_p)以及这些角点在步骤二所述世界坐标系中对应的二维点坐标值(X_w,Y_w)，并使用张正友的相机标定方法中线性模型求解方法，求解其对应在步骤六中所述的平面标定线性模型，完成对步骤四所述圆形区域1的线性标定。在上述求解过程中，同时求得步骤六所述的内部参数A，并且同时求得各组步骤一所述每幅照片图像对应在步骤六所述的外部参数R和t；

步骤八：为保证像素平面中的点被反投影到相同的世界坐标平面内，仍然利用步骤七求出的内部参数A以及外部参数R和t，对步骤四所述属于圆环形区域2建立非线性畸变模型如下：

$s\left[\begin{array}{c}{x}_u\\ y_u\\ 1\end{array}\right]=\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]...\left(2\right)$

$\left[\begin{array}{c}{x}_u\\ y_u\end{array}\right]=(1+k_1{r}^2+k_2{r}^4)\left[\begin{array}{c}{x}_d\\ y_d\end{array}\right]+\left[\begin{array}{c}2p_1{x}_d{y}_d+p_2(r^2+2{x_d}^2)\\ p_1(r^2+2{y_d}^2)+2p_2{x}_d{y}_d\end{array}\right]...\left(3\right)$

$\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=\left[\begin{array}{ccc}\mathrm{\α}& \mathrm{\γ}& u_0\\ 0& \mathrm{\β}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_d\\ y_d\\ 1\end{array}\right]...\left(4\right)$

式（3）中的k₁,k₂,p₁,p₂表示图像物理坐标系中径向与切向畸变函数的系数；

步骤九：求解步骤八中的畸变系数k₁,k₂,p₁,p₂，包括如下步骤：

9-1）、通过使步骤四中的圆环形区域2内角点像素坐标反投影的世界坐标值与其对应在真实世界坐标值的差值的平方和最小，建立四个畸变系数k₁,k₂,p₁,p₂的目标函数：

式（5）中M表示角点的世界坐标，表示角点像素坐标的反投影世界坐标，m表示标定图像的数目，n表示每幅标定图像中所用角点的数目；

9-2）、通过以下步骤求解四个畸变系数k₁,k₂,p₁,p₂的优化初值；

9-2-1）、利用线性最小二乘法将步骤八中的式（3）展开，将步骤八所述的畸变模型写成如下矩阵形式：

$\left[\begin{array}{cccc}{x}_d(x_d^2+y_d^2)& x_d{(x_d^2+y_d^2)}^2& 2x_d{y}_d& {(x_d^2+y_d^2)}^2+2x_d^2\\ y_d(x_d^2+y_d^2)& y_d{(x_d^2+y_d^2)}^2& {(x_d^2+y_d^2)}^2+2y_d^2& 2x_d{y}_d\end{array}\right]\left[\begin{array}{c}{k}_1\\ k_2\\ p_1\\ p_2\end{array}\right]=\left[\begin{array}{c}{x}_u-x_d\\ y_u-y_d\end{array}\right]...\left(6\right)$

上式中，四个畸变系数k₁,k₂,p₁,p₂是待求解变量；

9-2-2）、利用步骤七中获得的外部参数R和t以及步骤五所述圆环形区域2内的角点像素检测坐标在其世界坐标系中对应的二维点坐标值(X_w,Y_w)，代入步骤八所述式（2）中，即可获得步骤五所述圆环形区域2内的角点的理想图像坐标(x_u,y_u)；

9-2-3）、通过将步骤五所述属于圆环形区域2内的图像角点像素检测坐标(x_p,y_p)和步骤七求得的内部参数A代入步骤八的式（4）中，即可获得步骤五所述属于圆环形区域2内的图像角点的真实图像坐标(x_d,y_d)；

9-2-4）、设步骤一所述过程共取得m幅标定照片图像，每幅图像中包含n_i个步骤四所述属于圆环形区域2内的图像角点，则根据式（6）可获得2m×n₀个方程，则这些方程的矩阵可写成如下表达形式：

DK=d……(7)

式（7）中D表示2mn₀×4的系数矩阵，d表示2mn₀维的向量，K=[k₁ k₂ p₁ p₂]^T,T为向量转置符号；

9-2-5）、利用线性最小二乘法将式（7）变形为

K=(D^TD)^-1D^Td……(8)

通过对式（8）求解可以解出K的初始值；

9-3）将步骤9-2-5）所取得的K的最优初始值输入matlab软件，并利用matlab软件中Levenberg–Marquardt程序的优化算法对步骤9-1）所述的目标函数F(k₁,k₂,p₁,p₂)进行最小寻优计算，得到使目标函数F(k₁,k₂,p₁,p₂)最小的k₁,k₂,p₁,p₂解；

根据已求得的畸变系数k₁,k₂,p₁,p₂代入步骤八所述的由式（2）、式（3）和式（4）构成的整体非线性畸变模型，以完成对步骤四所述圆环形区域2的标定；

步骤十：获取步骤四所述属于第三区域3内的图像的角点的理想像素坐标(u_d,v_d)。将步骤一所取得的每一幅照片图像中属于第三区域3内的图像的角点在世界坐标系中对应的二维点坐标值(X_w,Y_w)作为转换目标，利用在步骤七获得的圆形区域1内对应平面的标定线性模型的内部参数矩阵A和每幅照片图像对应的外部参数R和t对前述的二维点坐标值(X_w,Y_w)投影到像素坐标系平面，以获得它们对应的理想参考像素坐标(u_d,v_d)；

步骤十一：将步骤十所述理想参考像素坐标(u_d,v_d)与步骤五属于第三区域3内的图像角点像素检测坐标(u_p,v_p)作差比较，获得属于第三区域3内的图像角点像素检测坐标在像素平面中的离散偏差分布，如下式：

δ_u(u_p,v_p)=u_d-u_p

δ_v(u_p,v_p)=v_d-v_p……(9)

式（9）中δ_u和δ_v分别表示在角点像素坐标的像素平面中U轴方向和V轴方向的离散偏差分布。

步骤十二：拟合获得连续的样条平滑函数：利用matlab中的样条平滑算法对式（9）所述的两个离散偏差分布进行运算处理，获得步骤四所述属于像素平面中第三区域3内图像像素点分别在U轴方向和V轴方向上的样条平滑修正函数：

δ_u=f_u(u,v)

δ_v=f_v(u,v)……(10)

式（10）中(u,v)表示第三区域内任意像素点的坐标，f_u表示像素点在像素平面U轴方向上的偏差修正函数，f_v表示像素点在像素平面V轴方向上的偏差修正函数；

步骤十三：在像素坐标平面内，利用步骤十二中取得的连续的样条平滑函数，对第三区域3内的图像像素点的畸变进行修正。修正完成后，可获得步骤十所述理想参考像素坐标(u_d,v_d)，最终仍然利用步骤七获得的内部参数A以及每幅图像属于圆形区域内的部分所对应外部参数R和t将上述理想坐标点反投影成对应的世界坐标点，完成步骤四所述属于像素平面中第三区域3的标定，进而完成摄像机的整个标定过程。

实施例：

以JAI CV-M4+CL型摄像机的标定过程为例，标定过程中使用的CCD摄像机的型号是：JAI CV-M4+CL，其感光单元面积为6.45μm×6.45μm；所使用的光学镜头型号为：ComputarM2514-MP；焦距是：25mm；所使用的棋盘格型标定板型号是：NANO CBC 75mm-2.0，其精度为1μm，方格大小为2×2mm，其被拍摄形成的一幅照片图像如图3中标号4所示。

在对该相机进行标定的过程中：进行步骤4-1）确定圆形区域时，设摄像机拍摄分辨率为1376×1024像素，将图像像素平面分辨率的中心点（688，512）设定为圆心，并在像素平面上定义一个R1=200像素为半径的圆形区域1。

进行步骤4-2）确定圆环形区域时，以R₁=200像素为内径，R₂=400像素为外径，在图像像素平面上定义一个圆环形区域2，其圆心与圆形区域1同心。

进行步骤4-3）确定第三区域时，在图像像素平面上，将圆环形区域2以外的图像区域定义为第三区域3。

进行步骤七时，求得的内部参数矩阵为：

$A=\left[\begin{array}{ccc}4189.05& 0.06& 680.99\\ 0& 4188.84& 501.89\\ 0& 0& 1\end{array}\right]$

进行步骤9-3）时，求得的K=[-0.3169,49.6091,-0.0005,-0.0001]；

进行步骤十二拟合获得连续的样条平滑函数时，利用在matlab软件中的csaps样条平滑算法程序对式（9）所述的两个离散偏差分布进行运算处理，以获得步骤四所述属于像素平面中第三区域3内图像像素点分别在U轴方向和V轴方向上的样条平滑修正函数。

进行步骤十二时，如图4所示，所获得的f_u表示像素点在像素平面U轴方向上的偏差修正函数，其在matlab软件中可显示为第三区域内像素点在U轴方向上的偏差曲面。如图5所示，所获得的f_v表示像素点在像素平面V轴方向上的偏差修正函数，其在matlab软件中可显示为第三区域内像素点在V轴方向上的偏差曲面。

本发明的摄像机标定方法计算标定模型参数的图像特征点分别在像素平面中划分的三个区域内，标定模型是分段函数；本发明采用畸变较少的圆形区域内的特征点计算线性模型，求解并在后续标定模型中延用这些准确性高的内部参数矩阵A和外部参数R、t；仅需对圆环区域模型中的畸变系数求解，并对畸变系数的优化初始值进行求解，降低了待求解的维数，避免了计算结果的不收敛情况并大幅缩短了计算时间；本发明将第三区域所建立的畸变模型转化为连续的样条平滑函数以修正第三区域内的全部像素点畸变，使整幅图像中距图像中线较远位置的特殊畸变都能够被精确修正，大幅度地提高了标定精度。

Claims (3)

1.一种基于分段模型的摄像机标定方法，其特征在于，该方法包括如下步骤：

步骤一：固定摄像机并将棋盘格标定板的中心置于摄像机光轴附近，使棋盘格标定板在一个相对固定的成像位置围绕摄像机镜头变换不同的倾角和姿态，用摄像机对标定板的每一种姿态拍摄一幅照片图像，最终取得多幅照片图像；

步骤二：按照TsaiRY的坐标转换方法确立转换坐标系，即：建立世界坐标系、摄像机坐标系、图像物理坐标系和图像像素平面坐标系，世界坐标系内的三维坐标点坐标(X_w,Y_w,Z_w)被依次转换为摄像机点坐标(X_c,Y_c,Z_c)、理想图像点坐标(x_u,y_u)、真实图像点坐标(x_d,y_d)和像素点坐标(x_p,y_p)；另外，参照张正友的摄像机标定方法将世界坐标系的Z_w轴坐标设为0，即Z_w＝0；此时，世界坐标系内任意点的坐标表示为(X_w,Y_w,0)，或仅表示为世界坐标系XY轴平面的二维点坐标(X_w,Y_w)；

步骤三：利用亚像素角点检测方法分别对步骤一所取得的各幅照片图像提取对应的角点像素坐标；

步骤四：将步骤二所述的图像像素坐标平面划分成三个相邻的区域：圆形区域(1)、圆环形区域(2)和第三区域(3)；

步骤五：将步骤三所取得的各幅照片图像的角点像素坐标按照其在步骤四所属的三个不同区域分类，分别得到属于圆形区域(1)内的图像角点像素检测坐标、属于圆环形区域(2)内的图像角点像素检测坐标和属于第三区域(3)内的图像角点像素检测坐标；

步骤六：按照步骤二所建立的转换坐标系在步骤四所取得的圆形区域(1)建立线性标定模型，将棋盘格标定板上方格角点在世界坐标系对应的三维点坐标(X_w,Y_w,0)转换为属于圆形区域(1)内的图像像素坐标(x_p,y_p)，忽略畸变的影响，所建立的圆形区域(1)内对应平面的标定线性模型如下：

$s\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=A\left[\begin{array}{cccc}{r}_1& r_2& r_3& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 0\\ 1\end{array}\right]=A[R,t]\left[\begin{array}{c}{X}_w\\ y_W\\ 0\\ 1\end{array}\right]......\left(1\right)$

当世界坐标系的Z_w轴坐标设为0时，即世界坐标系内任意点的坐标表示为(X_w,Y_w)时，式(1)可变换为下式

$s\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=A\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]=H\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]......\left(1.1\right)$

式(1)和式(1.1)中，s表示比例因子，R＝[r1,r2,r3]是世界坐标的旋转矩阵，r_i(i＝1,2,3)表示旋转矩阵R的第i列，t表示平移向量，R和t统称为标定模型的外部参数；

$A=\left[\begin{array}{ccc}\mathrm{\α}& \mathrm{\γ}& u_0\\ 0& \mathrm{\β}& v_0\\ 0& 0& 1\end{array}\right]$ 是内部参数矩阵，H表示3×3的单应性矩阵；α和β表示像素平面中U轴和V轴的比例因子，γ表示像素平面两坐标轴的不垂直因子，(u₀,v₀)表示相机光轴与图像平面的交点在像素坐标系上的坐标；

步骤七：利用步骤五所述圆形区域(1)内的角点像素检测坐标(x_p,y_p)以及这些角点在步骤二所述世界坐标系中对应的二维点坐标值(X_w,Y_w)，并使用张正友的相机标定方法中线性模型求解方法，求解其对应在步骤六中所述的平面标定线性模型，完成对步骤四所述圆形区域(1)的线性标定；在上述求解过程中，同时求得步骤六所述的内部参数A，并且同时求得各组步骤一中每幅照片图像对应在步骤六所述的外部参数R和t；

步骤八：利用步骤七求出的内部参数A以及外部参数R和t，对属于步骤四所述圆环形区域(2)建立非线性畸变模型如下：

$s\left[\begin{array}{c}{x}_u\\ y_u\\ 1\end{array}\right]=\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]......\left(2\right)$

$\left[\begin{array}{c}{x}_u\\ y_u\end{array}\right]=(1+k_1{r}^2+k_2{r}^4)\left[\begin{array}{c}{x}_d\\ y_d\end{array}\right]\left[\begin{array}{c}2p_1{x}_d{y}_d+p_2(r^2+2{x_d}^2)\\ p_1(r^2+2{y_d}^2)+2p_2{x}_d{y}_d\end{array}\right]......\left(3\right)$

$\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=\left[\begin{array}{ccc}\mathrm{\α}& \mathrm{\γ}& u_0\\ 0& \mathrm{\β}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_d\\ y_d\\ 1\end{array}\right]......\left(4\right)$

式(3)中的 $r=\sqrt{{x_u}^2+{y_u}^2};$

k₁,k₂,p₁,p₂表示图像物理坐标系中径向与切向畸变函数的系数；

步骤九：求解步骤八中的畸变系数k₁,k₂,p₁,p₂，并根据已求得的畸变系数k₁,k₂,p₁,p₂代入步骤八的由式(2)、式(3)和式(4)构成的整体非线性畸变模型，以完成对步骤四所述圆环形区域(2)的标定；

步骤十：获取步骤四所述属于第三区域(3)内的图像的角点的理想像素坐标(u_d,v_d)；将步骤一所取得的每一幅照片图像中属于第三区域(3)内的图像的角点在世界坐标系中对应的二维点坐标值(X_w,Y_w)作为转换目标，利用在步骤七获得的圆形区域(1)内对应平面的标定线性模型的内部参数矩阵A和每幅照片图像对应的外部参数R和t对前述的二维点坐标值(X_w,Y_w)投影到像素坐标系平面，以获得它们对应的理想参考像素坐标(u_d,v_d)；

步骤十一：将步骤十所述理想参考像素坐标(u_d,v_d)与步骤五属于第三区域(3)内的图像角点像素检测坐标(u_p,v_p)作差比较，获得属于第三区域(3)内的图像角点像素检测坐标在像素平面中的离散偏差分布，如下式：

δ_u(u_p,v_p)＝u_d-u_p

                    ……(9)

δ_v(u_p,v_p)＝v_d-v_p

式(9)中δ_u和δ_v分别表示在角点像素坐标的像素平面中U轴方向和V轴方向的离散偏差分布；

步骤十二：拟合获得连续的样条平滑函数：利用matlab中的样条平滑算法对式(9)所述的两个离散偏差分布进行运算处理，获得属于所述步骤四像素平面中第三区域(3)内图像像素点分别在U轴方向和V轴方向上的样条平滑修正函数：

δ_u＝f_u(u,v)

                  ……(10)

δ_v＝f_v(u,v)

式(10)中(u,v)表示第三区域内任意像素点的坐标，f_u表示像素点在像素平面U轴方向上的偏差修正函数，f_v表示像素点在像素平面V轴方向上的偏差修正函数；

步骤十三：在像素坐标平面内，利用步骤十二中取得的连续的样条平滑函数，对第三区域(3)内的图像像素点的畸变进行修正；修正完成后，可获得步骤十所述理想参考像素坐标(u_d,v_d)，最终可利用步骤七获得的内部参数A以及每幅图像对应外部参数R和t将上述理想坐标点反投影成对应的世界坐标点，完成属于所述步骤四像素平面中第三区域(3)的标定，进而完成摄像机的整个标定过程。

2.根据权利要求1所述的一种基于分段模型的摄像机标定方法，其特征在于，所述步骤四还包括如下步骤：

4-1)、确定圆形区域：设摄像机拍摄分辨率为M列×N行，将图像像素平面分辨率的中心点设定为圆心，并在像素平面上定义一个R₁为半径的圆形区域(1)；

4-2)、确定圆环形区域：以R₁为内径，R₂为外径，在图像像素平面上定义一个圆环形区域(2)，其圆心与圆形区域(1)同心；

4-3)、确定第三区域：在图像像素平面上，将圆环形区域(2)以外的图像区域定义为第三区域(3)。

3.根据权利要求1所述的一种基于分段模型的摄像机标定方法，其特征在于，所述步骤九还包括如下步骤：

9-1)、通过使步骤四中的圆环形区域(2)内角点像素坐标反投影的世界坐标值与其对应在真实世界坐标值的差值的平方和最小，建立四个畸变系数k₁,k₂,p₁,p₂的目标函数：

式(5)中M表示角点的世界坐标，表示角点像素坐标的反投影世界坐标，m表示标定图像的数目，n表示每幅标定图像中所用角点的数目；

9-2)、通过以下步骤求解四个畸变系数k₁,k₂,p₁,p₂的优化初值；

9-2-1)、利用线性最小二乘法将步骤八中的式(3)展开，将步骤八所述的畸变模型写成如下矩阵形式：

$\left[\begin{array}{cccc}{x}_d(x_d^2+y_d^2)& x_d{(x_d^2+y_d^2)}^2& 2x_d{y}_d& {(x_d^2+y_d^2)}^2+2x_d^2\\ y_d(x_d^2+y_d^2)& y_d{(x_d^2+y_d^2)}^2& {(x_d^2+y_d^2)}^2+2y_d^2& 2x_d{y}_d\end{array}\right]\left[\begin{array}{c}{k}_1\\ k_2\\ p_1\\ p_2\end{array}\right]=\left[\begin{array}{c}{x}_u-x_d\\ y_u-y_d\end{array}\right]......\left(6\right)$

上式中，四个畸变系数k₁,k₂,p₁,p₂是待求解变量；

9-2-2)、利用步骤七中获得的外部参数R和t以及步骤五所述圆环形区域(2)内的角点像素检测坐标在其世界坐标系中对应的二维点坐标值(X_w,Y_w)，代入步骤八所述式(2)中，即可获得步骤五所述圆环形区域(2)内的角点的理想图像坐标(x_u,y_u)；

9-2-3)、通过将步骤五所述属于圆环形区域(2)内的图像角点像素检测坐标(x_p,y_p)和步骤七求得的内部参数A代入步骤八的式(4)中，即可获得步骤五所述属于圆环形区域(2)内的图像角点的真实图像坐标(x_d,y_d)；

9-2-4)、设步骤一所述过程共取得m幅标定照片图像，每幅图像中包含n_i个属于所述步骤四中圆环形区域(2)内的图像角点，则根据式(6)可获得2m×n₀个方程，则这些方程的矩阵可写成如下表达形式：

DK＝d……(7)

式(7)中D表示2mn₀×4的系数矩阵，d表示2mn0维的向量，K＝[k₁ k₂ p₁ p₂]^T,T为向量转置符号；

9-2-5)、利用线性最小二乘法将式(7)变形为

K＝(D^TD)^-1D^Td……(8)

通过对式(8)求解可以解出K的初始值：

9-3)根据步骤9-2-5)所取得的K的最优初始值，利用matlab中的优化算法对步骤9-1)所述的目标函数F(k₁,k₂,p₁,p₂)进行最小寻优计算，得到使目标函数F(k₁,k₂,p₁,p₂)最小的k₁,k₂,p₁,p₂解。