CN106780628A - 基于混合畸变模型的高精度摄像机标定方法 - Google Patents

Abstract

基于混合畸变模型的高精度摄像机标定方法，包括如下步骤：S1.建立摄像机畸变模型；S2.标定图案选取、标定特征提取及特征点自动编号；S3.摄像机标定。本发明全面考虑了各种畸变，并基于提出的准则对伪特征点进行滤除，而后基于提出的圆心自动编号方法实现标志点的自动定位。

Description

基于混合畸变模型的高精度摄像机标定方法

技术领域

本发明属于信号处理领域，涉及标定技术领域的高精度锡膏测量***摄像机标定方法。

背景技术

结构光三维测量***是通过采集物体的二维图像信息，在二维信息的基础之上建立相应物体的三维形貌，三维形貌的准确性是以规范的和精确的二维图像信息为基础的。因此，规范计算机视觉***参数对于结构光投影三维测量前期的二维图像采集及后期的三维形貌的重构都十分重要。在大多数情况下，这些参数必须通过实验和计算才能得到，这个过程称为摄像机标定。摄像机标定作为机器视觉的主要内容之一，是机器视觉实现的前提和关键步骤。随着机器视觉技术的不断普及，对摄像机标定技术进行深入研究有着重要的意义，并受到了国内外学者的高度重视和广泛关注。在过去二十多年的发展中，摄像机的标定技术得到了很大的发展，许多学者根据实际应用具体需要提出了一些不同的摄像机模型和标定方法，基于不同的出发点和思路取得了一系列的成果。目前摄像机标定领域的学术研究思想非常的活跃，新技术新方法也不断涌现。

对于摄像机的标定的研究有很多，主要分为传统摄像机标定法、自标定法以及基于主动视觉的摄像机标定法三类。利用最为广泛的摄像机标定方法为Tsai法和张正友法。Tsai的方法依赖于对确定好参考坐标系的外部3-D校准物的精确3-D测量。所有位置和方位的测量都是相对于这个参考坐标系的。该方法己被广泛地用于在多摄像头***。张的方法采用的是平面的校准物。该方法的优点是，与3-D的校准物相比平面校准物更容易构造，而且相机的外在(相机与世界或相机与相机构成的几何形状)和固有参数(内部参数)可以在同一***中进行标定。Sturm和Maybank也提出类似的方法，对标定结果的奇异性进行了讨论。Brown提出了用以补偿镜头畸变的8参数模型，后来研究学者有发展了10参数模型和多参数模型。张正友以及Tsai使用的模型都是较为简单的参数模型。2001、2002年Luhmann和Hastedt等提出了一种基于有限元的混合畸变差模型。2004年冯文颜教授等提出基于二维直接线性变换的数字畸变差模型。

摄像机标定技术不仅受到国外学者的研究，国内的研究者们也极度关注和广泛重视摄像机标定技术。近年来，中国科学院自动化所的马颂德等人对摄像机标定技术有着标志性的见解，代表了我们国内的摄像机标定技术的研究水平。虽然在国际应用领域还没有一定的影响力，但表示国内已经具备了开展摄像机标定领域研究的人才和技术条件，并且为推动摄像机标定技术的研究和应用起到了模范的作用。此外，国内很多重点实验室也对摄像机标定做了广泛的研究，比如清华大学的智能技术与***国家重点实验室、中国科学院的模式识别国家重点实验室和天津大学的精密测试技术与仪器国家重点实验室等。

如何提高***标定的精度和效率一直都是众多学者所研究的热点和难点。尽管不少研究学者提出了不同的摄像机标定方法。但是大多数现有的摄像机标定方法所采用的畸变模型都是比较简单的畸变模型，还有些学者直接采用的线性模型根本都不考虑畸变，很少有人使用比较复杂的畸变模型对摄像机进行标定。

发明内容

本发明目的在于解决贴近物理现实条件的高精度锡膏测量***摄像机标定问题，提出了一种基于混合畸变模型的高精度摄像机标定方法，技术方案如下：

一种基于混合畸变模型的高精度摄像机标定方法，包括如下步骤：

S1.建立摄像机畸变模型；

S2.标定图案选取、标定特征提取及特征点自动编号；

S3.摄像机标定。

进一步的，所述建立摄像机畸变模型，包括

1)小孔成像模型：在该模型中建立四个基本坐标系，包括：世界坐标系O_wX_wY_wZ_w，相机坐标系O_cX_cY_cZ_c，图像坐标系xy，像素坐标系uv，摄像机坐标系(X_c,Y_c,Z_c)以摄像机镜头光心O_c为坐标原点，X_c,Y_c轴平行于图像平面，Z_c轴垂直于图像平面，于图像平面的交点在图像坐标系上的坐标为(u₀,v₀)，其摄像机的主点，摄像机镜头光心与主点之间的距离就是焦距f；

基于上述坐标系，可得，空间一点P(X_w,Y_w,Z_w)的位置在摄像机坐标系下可表示为：

其中为3×3正交旋转矩阵，为3×1平移矩阵；

基于小孔成像模型，得摄像机坐标(X_c,Y_c,Z_c,1)^T在像素坐标系下可表示为：

其中，(u,v,1)^T为图像像素坐标，ρ_x,ρ_y分别为像平面中x,y方向上单位距离的像素数，(u₀,v₀)为摄像机主点坐标，其是像平面中心的像素坐标，α为描述像素坐标轴垂直误差的偏斜因子；

由以上两式，得图像坐标与世界坐标的关系式：

其中，为摄像机内参数矩阵，E＝[R T]为外参数矩阵；

2)构建摄像机畸变模型：采用基于矩形有限元方法的混合畸变模型，此模型描述如下：

其中，dx,dy分别为像点x坐标和y坐标的偏差；

其中，x_i,y_i分别为符合小孔成像的理想无畸变像平面坐标，K_c1,K_c2,K_c3为三阶径向畸变系数；

其中，P_c1,P_c2为二阶切向畸变系数；

其中，S_c1,S_c2为二阶薄棱镜畸变系数；

其中，δx,δy分别为像点在矩形网格x,y方向上的长度比。

进一步的，所述步骤S2中标定图案选取，其采用的标定板是铝制标定板，使用其对摄像机进行参数标定，标定板上有9×11个圆形标定点，其中小圆直径为4mm，大圆直径为8mm，五个大圆用于在标定过程中标识标定板方向。

进一步的，标定特征提取包括以下步骤：图像采集、图像灰度化、图像滤波、图像二值化、轮廓提取、伪点滤除、椭圆拟合、圆心提取；

利用OpenCV函数库里的函数直接将彩色图像读入变为灰色图像，紧接着对获取的灰色图像进行滤波，使用函数cvSmooth实现图像滤波，图像二值化中，阈值经试验可设为100，轮廓提取采用基于迭代的改进Harris角点检测法，此方法描述如下：考虑角点q邻域里的任一点p，▽I(p)表示p点处的图像梯度矢量，则由梯度性质可得：

▽I(p)(q-p)＝0

考虑角点q邻域里所有p点，则构成一个超定方程组，从而利用最小均方误差(MMSE)准则进行求解；

提取轮廓后，利用如下规则消除伪特征点：

面积准则：提取特征对象的面积根据***配置进行范围确定：

S_min＜S＜S_max

圆度准则：用下式衡量形状与圆的接近程度：

其中T为圆周长，

误差准则：将边缘轮廊点到拟合椭圆的代数距离作为误差，同时采用以下准则进一步滤除误差较大的圆形特征点：

其中表示轮廓点的平均误差，ε_max表示最大允许误差，ε_mean表示允许平均误差；

进一步的，伪特征点消除后，采用加权最小二乘法进行轮廓椭圆拟合提取亚像素级精度的圆心坐标。

所述步骤S2中特征点自动编号：利用标定板上的5个实心圆，每个实心圆都对应一个确定不变的编号，进行如下步骤：

(1)先用鼠标选取空间某矩形区域在图像上对应的四个角点，获取四个角点的图像坐标并经亚像素处理；

(2)对选定区域的圆进行自动记数，并根据计数结果生成网格图。

所述摄像机标定的方法如下：

1)单应性矩阵求解

假设模板位于世界坐标系Z＝0的位置，由针孔模型，得

sm＝HM

其中，s为尺度系数，m,M分别为图像平面以及对应的模板上的点，H＝K_c[R T]；

基于平面模板各角点空间坐标，通过求解如下的MMSE估计，得到H的估计值，

其中， 为H矩阵的第i个行向量；

采用改进的Leverberg-Marquardt法求解，用非线性优化方法求解需要有一个用来迭代的合适初始值，通过求解下面方程获得：

给定n＞6对点，就形成超定方程组，则通过奇异值分解(SVD)进行求解；

2)内外参数求解

基于内参数矩阵R为正交矩阵的性质，通过求解如下超定方程组得到摄像机的内外参数：

Vb＝0

此方程组为n≥3次观察得到的如下方程的叠加：

式中v_ij＝[h_1ih_1j,h_1ih_2j+h_2ih_1j,h_2ih_2j,h_3ih_1j+h_1ih_3j,h_3ih_2j+h_2ih_3j,h_3ih_3j]^T，b＝[B₁₁,B₁₂,B₂₂,B₁₃,B₂₃,B₃₃]^T，

而后，通过最大似然估计对此解进行精细化，以得到更好的优化解，似然估计问题表示如下：

式中，m_ij为第i幅图像上第j个像点；为空间点M_j在第i幅图像的投影点；

3)畸变参数的求解；

4)有限元模型畸变参数求解；

5)改进LM算法。

有益效果：

本发明针对现有标定算法所采用畸变模型较为简单，从而不能实现高精度标定的问题，提出了基于混合畸变模型的摄像机标定方法。此方法首先基于小孔成像建立成像模型，而后针对三种畸变以及平面不均等物理误差建立了混合畸变模型。为标定摄像机内外参数以及畸变参数，本发明确定了可实现自动定位的标定模板，而后基于此模板实现标志点的亚像素定位。在获得标志点的精确定位条件下，本发明首先通过求解超定方程组得到单应矩阵的初始解，而后利用所提出的改进LM方法进行高效优化；基于此单应矩阵，可求得摄像机的内外参数，并可利用最大似然估计实现进一步优化；在求得内外参数后，通过求解超定方程而后利用改进的LM方法进一步优化，可分别获得摄像机的三种畸变参数以及有限元畸变参数。与传统的标定方法相比，该本发明全面考虑了各种畸变，并基于提出的准则对伪特征点进行滤除，而后基于提出的圆心自动编号方法实现标志点的自动定位，此外，在求解标定参数的过程中，基于改进的LM算法进行优化，从而显著加快了收敛速度并提高了标定精度。基于以上讨论可知，本发明所提方法可为机器视觉、三维测量等工程应用中高精度相机标定提供坚实的理论与实现依据。

附图说明

图1为本发明实现流程图；

图2为小孔成像原理图；

图3为本发明采用的标定图案；

图4为特征点提取流程图；

图5为本专利拍摄标定图像；

图6为边缘检测流程图；

图7为标志点编号图；

图8为特征点自动编号结果；

图9摄像机标定误差分布图。

具体实施方式

实施例1：

基于混合畸变模型的高精度摄像机标定方法，首先通过求解超定方程组得到单应矩阵的初始解，而后利用所提出的改进LM方法进行高效优化；基于此单应矩阵，可求得摄像机的内外参数，并可利用最大似然估计实现进一步优化；在求得内外参数后，通过求解超定方程而后利用改进的LM方法进一步优化，可分别获得摄像机的三种畸变参数以及有限元畸变参数。从而获得相机混合畸变参数的精确估计，进而实现摄像机的高精度标定。

实现本发明的基本思路是，首先基于小孔成像建立成像模型，而后针对三种畸变以及平面不均等物理误差建立了混合畸变模型。为标定摄像机内外参数以及畸变参数，本发明确定了可实现自动定位的标定模板，而后基于此模板实现标志点的亚像素定位。在获得标志点的精确定位条件下，本发明首先通过求解超定方程组得到单应矩阵的初始解，而后利用所提出的改进LM方法进行高效优化；基于此单应矩阵，可求得摄像机的内外参数，并可利用最大似然估计实现进一步优化；在求得内外参数后，通过求解超定方程而后利用改进的LM方法进一步优化，可分别获得摄像机的三种畸变参数以及有限元畸变参数。与传统的标定方法相比，该本发明全面考虑了各种畸变，并基于提出的准则对伪特征点进行滤除，而后基于提出的圆心自动编号方法实现标志点的自动定位，此外，在求解标定参数的过程中，基于改进的LM算法进行优化，从而显著加快了收敛速度并提高了标定精度。具体步骤包含如下：

1.建立摄像机畸变模型

1)构建小孔成像模型

相机模型是光学成像几何关系的简化，许多相机成像几何模型都是根据小孔成像原理得到的。由于小孔成像的原理较简单而使针孔成像模型成为最简单的相机模型，它是相机标定算法的基本模型。若再适当地考虑其镜头本身的畸变因素，就能满足很多应用场合所需要的精度。本专利中采用的是如图2所示的理想情况下的针孔模型，在此基础上再考虑镜头畸变的影响，就能较精确地反映成像过程。世界坐标系(X_w,Y_w,Z_w)是在环境中选择的一个基准坐标系，用来描述摄像机的位置，可以根据描述和计算的方便等原则来自由选取。在此模型中建立了四个基本坐标系，包括：世界坐标系O_wX_wY_wZ_w，相机坐标系O_cX_cY_cZ_c，图像坐标系xy，像素坐标系uv。由空间物点的三维坐标到图像的成像过程就是这几个坐标系逐步转换的过程。摄像机坐标系(X_c,Y_c,Z_c)以摄像机镜头光心O_c为坐标原点，X_c,Y_c轴平行于图像平面，Z_c轴垂直于图像平面，于图像平面的交点在图像坐标系上的坐标为(u₀,v₀)，即摄像机的主点。需要注意的是该点一般位于图像平面的中心处，但由于摄像机的制作原因，有时也会偏离，因此摄像机主点坐标一般也是需要标定的两个参数。摄像机镜头光心与主点之间的距离就是焦距f。

基于上述坐标系，可得，空间一点P(X_w,Y_w,Z_w)的位置在摄像机坐标系下可表示为：

其中为3×3正交旋转矩阵，为3×1平移矩阵。

基于摄像机针孔成像模型，可得摄像机坐标(X_c,Y_c,Z_c,1)^T在像素坐标系下可表示为：

其中，(u,v,1)^T为图像像素坐标，ρ_x,ρ_y分别为像平面中x,y方向上单位距离的像素数，(u₀,v₀)为摄像机主点坐标，即像平面中心的像素坐标，α为描述像素坐标轴垂直误差的偏斜因子。

由以上两式，可得图像坐标与世界坐标的关系式，

其中，为摄像机内参数矩阵，E＝[R T]为外参数矩阵。

2)摄像机畸变模型

上一部分摄像机模型为不考虑镜头畸变的线性模型，但是由于镜头设计、制作、装配等因素的影响，实际成像***不可能严格的满足小孔成像原理，会存在畸变使光线发生细微偏移。且有大量实验表明，上一部分给出的摄像机针孔模型不能准确地描述摄像机的成像几何关系，尤其是在使用广角镜头时，图像中远离图像中心的像点会有很大的畸变。摄像机镜头的非线性光学畸变主要有三类：径向畸变(Radial Distortion)、偏心畸变(Decentering Distortion)以及薄棱镜畸变(Prism Distortion)。其中径向畸变使像点沿着像主点径向向内或向外移动，切向畸变是使像点沿切向移动，偏心畸变及薄棱镜畸变既会产生径向位置的畸变又会产生切向位置的畸变。

参数模型对于径向畸变、切向畸变以及薄棱镜畸变的描述非常准确，且未知参数较少，求解相对简单，但无法描述像平面不平等因素引起的畸变。有限元模型可以描述各种畸变，但是考虑的未知参数却很多，求解相对困难。而有限元模型综合了两种模型的优势，对于给畸变造成影响较大的三大畸变可以准确的解算出，对于像平面等畸变同样可以进行补偿。因此，本发明采用基于矩形有限元方法的混合畸变模型，此模型可描述如下：

其中，dx,dy分别为像点x坐标和y坐标的偏差；

其中，x_i,y_i分别为符合小孔成像的理想无畸变像平面坐标，K_c1,K_c2,K_c3为三阶径向畸变系数；

其中，P_c1,P_c2为二阶切向畸变系数；

其中，S_c1,S_c2为二阶薄棱镜畸变系数；

其中，δx,δy分别为像点在矩形网格x,y方向上的长度比。

2.标定图案选取、标定特征提取及特征点自动编号

1)标定图案选取

通常情况下二维平面标定物的设计应满足如下基本要求：第一，标定物的图像特征应易于识别，即图像上的标志点在色彩上与背景之间应存在鲜明的差别，并要尽可能的清晰和均匀；第二，标定物上的标志点在进行图像处理时应要易于提取，即能够提供质量较髙的轮廓；第三，标定物应具有较高的稳定性，即其标志点不应随着摄像机位置的改变而产生较大范围的畸变，且要易于进行测量；第四，要实现自动化标定，设计的标志点应易于实现自动编码。

基于上述要求，本专利采用图3所示的铝制高精度标定板对摄像机进行参数标定，标定板上有9×11个圆形标定点，其中小圆直径为4mm，大圆直径为8mm，五个大圆用于在标定过程中标识标定板方向，从而可实现标定自动化。

2)标定特征提取

要准确提取出标定板图像中圆心坐标,图像处理技术是必不可少的。圆心提取过程的结构框图如图4所示。主要包括以下步骤：图像采集，图像灰度化，图像滤波，图像二值化，轮廓提取，伪点滤除，椭圆拟合，圆心提取。

本专利利用摄像机对标定板在不同位置及角度釆集了12幅图像，如图5所示。

本专利利用OpenCV函数库里的函数直接将彩色图像读入变为灰色图像，紧接着对获取的灰色图像进行滤波，函数cvSmooth可以有效的实现图像滤波。图像二值化中，阈值经试验可设为100。

轮廓提取可采用本专利提出的基于迭代的改进Harris角点检测法，此方法可描述如下：考虑角点q邻域里的任一点p，▽I(p)表示p点处的图像梯度矢量，则由梯度性质可得

▽I(p)(q-p)＝0

考虑角点q邻域里所有p点，则可构成一个超定方程组，从而可以利用最小均方误差(MMSE)准则进行求解。

令q^k＝[x_k,y_k]^T为角点初始估计点，则

解此优化问题，可得

q^k+1＝F-¹W

其中， X,Y为一阶灰度梯度，w_u,v为高斯滤波系数。

重复迭代，则可更新q^k+1，直至||q^k+1-q^k||不再减小。上述边缘检测方法如流程图6所示。

提取轮廓后，可利用如下规则消除伪特征点：

1、面积准则

由于成像大小仅与摄像机的分辨率和测量距离有关，因此提取特征对象的面积也可根据***配置进行范围确定，即

S_min＜S＜S_max

2、圆度准则

可用下式衡量形状与圆的接近程度：

其中T为圆周长，C越大说明形状与圆差别越大。本专利的圆度准则为

3、误差准则

本专利将边缘轮廊点到拟合椭圆的代数距离作为误差，同时采用以下准则进一步滤除误差较大的圆形特征点。

其中表示轮廓点的平均误差，ε_max表示最大允许误差，ε_mean表示允许平均误差。本专利设误差相关值如下：ε_max＝0.6pixel，ε_mean＝0.3pixel。

伪特征点消除后，本发明采用加权最小二乘法进行轮廓椭圆拟合，从而提取亚像素级精度的圆心坐标。

3)特征点自动编号

前面部分的操作已经提取感兴趣的椭圆中心坐标，但是无法将椭圆中心的坐标与标定板上的圆形标志点圆心坐标对应起来。要实现摄像机标定，将相应的椭圆中心坐标与标志点圆心对应起来是必要的，并且圆形标志点的自动编号是实现摄像机自动化标定的关键所在。

为了实现圆形标志点的自动编号，可利用标定板上的5个大实心圆，标定板中各实心圆的编号如图7所示，每个实心圆都对应一个确定不变的编号。

要实现对所有圆形标志点进行正确自动编号，需要经过以下步骤：

(3)先用鼠标选取空间某矩形区域在图像上对应的四个角点，获取四个角点的图像坐标并经亚像素处理。为了更多地利用图像信息，该区域应尽可能大，以包含更多的控制点。设此四点相应的空间平面模板上的点的齐次坐标为(0,0,1),(1,0,1),(1,1,1),(0,1,1)，这实际上是对空间平面模板上的点做了一次坐标变换。这一步同时建立起了世界坐标系。建立世界坐标系的方法是，将选取的第一个角点作为世界坐标系原点，Z轴垂直于方块格所在平面指向平面模板外，确定Z轴后，以右手系为标准确定X和Y轴。

(4)对选定区域的圆进行自动记数，并根据计数结果生成网格图。自动计数方法如下：a)大圆定位：1、提取5个大圆：首先筛选轮廓信息，只保留椭圆信息，然后根据椭圆拟合得到的椭圆参数计算椭圆的面积及周长，选择面积最大的5个。2、计算5个大圆相互之间的距离。3、大圆分组：距离最近的为A组26，27号，最远为B组46，52号，剩下的即为71号。4、大圆编号：分别计算组B中与组A中两圆的中心距离之和，大者为52号，另一则为46号，同理可得组A中大者为27号，另一位26号。b)小圆定位：在大圆定位基础上，可进行小圆定位。1、划定区域：由于过27号、71号以及46、52号大圆圆心的实际坐标以得到，可计算此二直线轨迹，然后沿该轨迹上下扩展形成局部图像(上下扩展比大圆大3个像素)。2、计算灰度曲线：计算此区域每个像素点灰度之和。3、滤除噪声：为防止曲线过零点有抖动，可进行滤波。3、根据小圆半径、***配置等先验知识，结合灰度曲线，即可实现小圆的自动编号。编号结果如图8所示，其中十字标标示的为椭圆圆心坐标。

3.摄像机标定

由于本专利采用的是混合畸变模型，因此摄像机的畸变主要包括两部分：利用参数描述的各类畸变和利用有限元模型描述的像平面不平等畸变。由于参数模型可以很准确的描述径向畸变、切向畸变和薄棱镜畸变，有限元模型描述的为每个像点由于像平面不平等因素引起的在x和y方向上的畸变,为此可将标定过程分成一般的参数模型的求解以及有限元模型的求解。

1)单应性矩阵求解

假设模板位于世界坐标系Z＝0的位置，由针孔模型，可得

sm＝HM

其中，s为尺度系数，m,M分别为图像平面以及对应的模板上的点，H＝K_c[R T]。

基于平面模板各角点空间坐标，可通过求解如下的MMSE估计，得到H的估计值，即

其中， 为H矩阵的第i个行向量。

此问题是一个非线性最小二乘问题，可以用梯度法、高斯-牛顿法或Leverberg-Marquardt法等多种非线性优化来解决。此处采用本专利提出的改进的Leverberg-Marquardt法求解(详见5))用非线性优化方法求解需要有一个用来迭代的合适初始值，可通过求解下面方程获得：

给定n＞6对点，就可以形成超定方程组，则可以通过奇异值分解(SVD)进行求解。

2)内外参数求解

基于内参数矩阵R为正交矩阵的性质，可通过求解如下超定方程组得到摄像机的内外参数：

Vb＝0

此方程组为n≥3次观察得到的如下方程的叠加：

式中v_ij＝[h_1ih_1j,h_1ih_2j+h_2ih_1j,h_2ih_2j,h_3ih_1j+h_1ih_3j,h_3ih_2j+h_2ih_3j,h_3ih_3j]^T，b＝[B₁₁,B₁₂,B₂₂,B₁₃,B₂₃,B₃₃]^T，

而后，可通过最大似然估计对此解进行精细化，以得到更好的优化解，似然估计问题可表示如下：

式中，m_ij为第i幅图像上第j个像点；为空间点M_j在第i幅图像的投影点。

3)畸变参数的求解

由上述畸变模型可得：

式中，(x_r,y_r)为实际的像平面坐标，(x_i,y_i)为理想的无畸变像平面坐标。

又知，

式中，(u_r,v_r)为实际的图像坐标，(u_i,v_i)为理想的无畸变图像坐标。

若每幅图像具有m个标志点，且采集了n幅图像，则基于以上两式，可得以下超定方程组：

GC＝D

其中，

利用最小二乘方法(LS)可求得向量C：

C＝(G^TG)^-1G^TD

求得此初始解后，可利用非线性方法进一步优化，从而可得畸变参数的精确解。

4)有限元模型畸变参数求解

根据上述基于有限元的畸变模型，利用获得的数据，可得超定方程组，从而可利用LS方法求得有限元模型畸变参数的初始解，而后利用改进的Levenberg-Marquartdt算法进行非线性优化。

5)改进LM算法

目前常用的非线性最小二乘法有三种：最速下降法、高斯-牛顿法和LM法。其中LM方法中，迭代增量公式如下：

Δθ＝(J^TJ+μI)^-1J^Tε_i

式中，μ为阻尼因子，J为偏导数矩阵，ε_i为估计误差。其中，μ的取值很大程度决定了LM方法的收敛速度以及求解的精确度。然而，实际应用中，μ的最优取值无法确定。为改善LM算法的收敛速度以及求解的精确度，本发明提出一种改进因子，此因子可随迭代次数变化而变化，使得迭代初期可放宽对初始值，而后期则可加大收敛速度，从而使得LM算法在迭代初期具有最速下降法，而在最优解附近又具有高斯-牛顿法的优势，此下降因子为：

1-(1-μ)ⁱ⁺¹

从而增量迭代公式可重新表示为：

Δθ＝(J^TJ+(1-(1-μ)ⁱ⁺¹)I)^-1J^Tε_i

实施例2：

一种基于混合畸变模型的高精度摄像机标定方法：包括

1.建立摄像机畸变模型

1)构建小孔成像模型

相机模型是光学成像几何关系的简化，许多相机成像几何模型都是根据小孔成像原理得到的。由于小孔成像的原理较简单而使针孔成像模型成为最简单的相机模型，它是相机标定算法的基本模型。若再适当地考虑其镜头本身的畸变因素，就能满足很多应用场合所需要的精度。本专利中采用的是如图2所示的理想情况下的针孔模型，在此基础上再考虑镜头畸变的影响，就能较精确地反映成像过程。世界坐标系(X_w,Y_w,Z_w)是在环境中选择的一个基准坐标系，用来描述摄像机的位置，可以根据描述和计算的方便等原则来自由选取。在此模型中建立了四个基本坐标系，包括：世界坐标系O_wX_wY_wZ_w，相机坐标系O_cX_cY_cZ_c，图像坐标系xy，像素坐标系uv。由空间物点的三维坐标到图像的成像过程就是这几个坐标系逐步转换的过程。摄像机坐标系(X_c,Y_c,Z_c)以摄像机镜头光心O_c为坐标原点，X_c,Y_c轴平行于图像平面，Z_c轴垂直于图像平面，于图像平面的交点在图像坐标系上的坐标为(u₀,v₀)，即摄像机的主点。需要注意的是该点一般位于图像平面的中心处，但由于摄像机的制作原因，有时也会偏离，因此摄像机主点坐标一般也是需要标定的两个参数。摄像机镜头光心与主点之间的距离就是焦距f。

基于上述坐标系，可得，空间一点P(X_w,Y_w,Z_w)的位置在摄像机坐标系下可表示为：

其中为3×3正交旋转矩阵，为3×1平移矩阵。

基于摄像机针孔成像模型，可得摄像机坐标(X_c,Y_c,Z_c,1)^T在像素坐标系下可表示为：

其中，(u,v,1)^T为图像像素坐标，ρ_x,ρ_y分别为像平面中x,y方向上单位距离的像素数，(u₀,v₀)为摄像机主点坐标，即像平面中心的像素坐标，α为描述像素坐标轴垂直误差的偏斜因子。

由以上两式，可得图像坐标与世界坐标的关系式，

其中，为摄像机内参数矩阵，E＝[R T]为外参数矩阵。

2)摄像机畸变模型

上一部分摄像机模型为不考虑镜头畸变的线性模型，但是由于镜头设计、制作、装配等因素的影响，实际成像***不可能严格的满足小孔成像原理，会存在畸变使光线发生细微偏移。且有大量实验表明，上一部分给出的摄像机针孔模型不能准确地描述摄像机的成像几何关系，尤其是在使用广角镜头时，图像中远离图像中心的像点会有很大的畸变。摄像机镜头的非线性光学畸变主要有三类：径向畸变(Radial Distortion)、偏心畸变(Decentering Distortion)以及薄棱镜畸变(Prism Distortion)。其中径向畸变使像点沿着像主点径向向内或向外移动，切向畸变是使像点沿切向移动，偏心畸变及薄棱镜畸变既会产生径向位置的畸变又会产生切向位置的畸变。

参数模型对于径向畸变、切向畸变以及薄棱镜畸变的描述非常准确，且未知参数较少，求解相对简单，但无法描述像平面不平等因素引起的畸变。有限元模型可以描述各种畸变，但是考虑的未知参数却很多，求解相对困难。而有限元模型综合了两种模型的优势，对于给畸变造成影响较大的三大畸变可以准确的解算出，对于像平面等畸变同样可以进行补偿。因此，本发明采用基于矩形有限元方法的混合畸变模型，此模型可描述如下：

其中，dx,dy分别为像点x坐标和y坐标的偏差；

其中，x_i,y_i分别为符合小孔成像的理想无畸变像平面坐标，K_c1,K_c2,K_c3为三阶径向畸变系数；

其中，P_c1,P_c2为二阶切向畸变系数；

其中，S_c1,S_c2为二阶薄棱镜畸变系数；

其中，δx,δy分别为像点在矩形网格x,y方向上的长度比。

2.标定图案选取、标定特征提取及特征点自动编号

1)标定图案选取

通常情况下二维平面标定物的设计应满足如下基本要求：第一，标定物的图像特征应易于识别，即图像上的标志点在色彩上与背景之间应存在鲜明的差别，并要尽可能的清晰和均匀；第二，标定物上的标志点在进行图像处理时应要易于提取，即能够提供质量较髙的轮廓；第三，标定物应具有较高的稳定性，即其标志点不应随着摄像机位置的改变而产生较大范围的畸变，且要易于进行测量；第四，要实现自动化标定，设计的标志点应易于实现自动编码。

基于上述要求，本专利采用图3所示的铝制高精度标定板对摄像机进行参数标定，标定板上有9×11个圆形标定点，其中小圆直径为4mm，大圆直径为8mm，五个大圆用于在标定过程中标识标定板方向，从而可实现标定自动化。

2)标定特征提取

要准确提取出标定板图像中圆心坐标,图像处理技术是必不可少的。圆心提取过程的结构框图如图4所示。主要包括以下步骤：图像采集，图像灰度化，图像滤波，图像二值化，轮廓提取，伪点滤除，椭圆拟合，圆心提取。

本专利利用摄像机对标定板在不同位置及角度釆集了12幅图像，如图5所示。

本专利利用OpenCV函数库里的函数直接将彩色图像读入变为灰色图像，紧接着对获取的灰色图像进行滤波，函数cvSmooth可以有效的实现图像滤波。图像二值化中，阈值经试验可设为100。

轮廓提取可采用本专利提出的基于迭代的改进Harris角点检测法，此方法可描述如下：考虑角点q邻域里的任一点p，▽I(p)表示p点处的图像梯度矢量，则由梯度性质可得

▽I(p)(q-p)＝0

考虑角点q邻域里所有p点，则可构成一个超定方程组，从而可以利用最小均方误差(MMSE)准则进行求解。

令q^k＝[x_k,y_k]^T为角点初始估计点，则

解此优化问题，可得

q^k+1＝F-¹W

其中， X,Y为一阶灰度梯度，w_u,v为高斯滤波系数。

重复迭代，则可更新q^k+1，直至||q^k+1-q^k||不再减小。上述边缘检测方法如流程图6所示。

提取轮廓后，可利用如下规则消除伪特征点：

4、面积准则

由于成像大小仅与摄像机的分辨率和测量距离有关，因此提取特征对象的面积也可根据***配置进行范围确定，即

S_min＜S＜S_max

5、圆度准则

可用下式衡量形状与圆的接近程度：

其中T为圆周长，C越大说明形状与圆差别越大。本专利的圆度准则为

6、误差准则

本专利将边缘轮廊点到拟合椭圆的代数距离作为误差，同时采用以下准则进一步滤除误差较大的圆形特征点。

其中表示轮廓点的平均误差，ε_max表示最大允许误差，ε_mean表示允许平均误差。本专利设误差相关值如下：ε_max＝0.6pixel，ε_mean＝0.3pixel。

伪特征点消除后，本发明采用加权最小二乘法进行轮廓椭圆拟合，从而提取亚像素级精度的圆心坐标。

3)特征点自动编号

前面部分的操作已经提取感兴趣的椭圆中心坐标，但是无法将椭圆中心的坐标与标定板上的圆形标志点圆心坐标对应起来。要实现摄像机标定，将相应的椭圆中心坐标与标志点圆心对应起来是必要的，并且圆形标志点的自动编号是实现摄像机自动化标定的关键所在。

为了实现圆形标志点的自动编号，可利用标定板上的5个大实心圆，标定板中各实心圆的编号如图7所示，每个实心圆都对应一个确定不变的编号。

要实现对所有圆形标志点进行正确自动编号，需要经过以下步骤：

(5)先用鼠标选取空间某矩形区域在图像上对应的四个角点，获取四个角点的图像坐标并经亚像素处理。为了更多地利用图像信息，该区域应尽可能大，以包含更多的控制点。设此四点相应的空间平面模板上的点的齐次坐标为(0,0,1),(1,0,1),(1,1,1),(0,1,1)，这实际上是对空间平面模板上的点做了一次坐标变换。这一步同时建立起了世界坐标系。建立世界坐标系的方法是，将选取的第一个角点作为世界坐标系原点，Z轴垂直于方块格所在平面指向平面模板外，确定Z轴后，以右手系为标准确定X和Y轴。

(6)对选定区域的圆进行自动记数，并根据计数结果生成网格图。自动计数方法如下：a)大圆定位：1、提取5个大圆：首先筛选轮廓信息，只保留椭圆信息，然后根据椭圆拟合得到的椭圆参数计算椭圆的面积及周长，选择面积最大的5个。2、计算5个大圆相互之间的距离。3、大圆分组：距离最近的为A组26，27号，最远为B组46，52号，剩下的即为71号。4、大圆编号：分别计算组B中与组A中两圆的中心距离之和，大者为52号，另一则为46号，同理可得组A中大者为27号，另一位26号。b)小圆定位：在大圆定位基础上，可进行小圆定位。1、划定区域：由于过27号、71号以及46、52号大圆圆心的实际坐标以得到，可计算此二直线轨迹，然后沿该轨迹上下扩展形成局部图像(上下扩展比大圆大3个像素)。2、计算灰度曲线：计算此区域每个像素点灰度之和。3、滤除噪声：为防止曲线过零点有抖动，可进行滤波。3、根据小圆半径、***配置等先验知识，结合灰度曲线，即可实现小圆的自动编号。编号结果如图8所示，其中十字标标示的为椭圆圆心坐标。

3.摄像机标定

由于本专利采用的是混合畸变模型，因此摄像机的畸变主要包括两部分：利用参数描述的各类畸变和利用有限元模型描述的像平面不平等畸变。由于参数模型可以很准确的描述径向畸变、切向畸变和薄棱镜畸变，有限元模型描述的为每个像点由于像平面不平等因素引起的在x和y方向上的畸变,为此可将标定过程分成一般的参数模型的求解以及有限元模型的求解。

1)单应性矩阵求解

假设模板位于世界坐标系Z＝0的位置，由针孔模型，可得

sm＝HM

其中，s为尺度系数，m,M分别为图像平面以及对应的模板上的点，H＝K_c[R T]。

基于平面模板各角点空间坐标，可通过求解如下的MMSE估计，得到H的估计值，即

其中， 为H矩阵的第i个行向量。

此问题是一个非线性最小二乘问题，可以用梯度法、高斯-牛顿法或Leverberg-Marquardt法等多种非线性优化来解决。此处采用本专利提出的改进的Leverberg-Marquardt法求解(详见5))用非线性优化方法求解需要有一个用来迭代的合适初始值，可通过求解下面方程获得：

给定n＞6对点，就可以形成超定方程组，则可以通过奇异值分解(SVD)进行求解。

2)内外参数求解

基于内参数矩阵R为正交矩阵的性质，可通过求解如下超定方程组得到摄像机的内外参数：

Vb＝0

此方程组为n≥3次观察得到的如下方程的叠加：

式中v_ij＝[h_1ih_1j,h_1ih_2j+h_2ih_1j,h_2ih_2j,h_3ih_1j+h_1ih_3j,h_3ih_2j+h_2ih_3j,h_3ih_3j]^T，b＝[B₁₁,B₁₂,B₂₂,B₁₃,B₂₃,B₃₃]^T，

而后，可通过最大似然估计对此解进行精细化，以得到更好的优化解，似然估计问题可表示如下：

式中，m_ij为第i幅图像上第j个像点；为空间点M_j在第i幅图像的投影点。

3)畸变参数的求解

由上述畸变模型可得：

式中，(x_r,y_r)为实际的像平面坐标，(x_i,y_i)为理想的无畸变像平面坐标。

又知，

式中，(u_r,v_r)为实际的图像坐标，(u_i,v_i)为理想的无畸变图像坐标。

若每幅图像具有m个标志点，且采集了n幅图像，则基于以上两式，可得以下超定方程组：

GC＝D

其中，

C＝[K_c1,K_c2,K_c3,P_c1,P_c2,S_c1,S_c2]^T，D＝[d₁d₂…d_mn]^T。

利用最小二乘方法(LS)可求得向量C：

C＝(G^TG)^-1G^TD

求得此初始解后，可利用非线性方法进一步优化，从而可得畸变参数的精确解。

4)有限元模型畸变参数求解

根据上述基于有限元的畸变模型，利用获得的数据，可得超定方程组，从而可利用LS方法求得有限元模型畸变参数的初始解，而后利用改进的Levenberg-Marquartdt算法进行非线性优化。

5)改进LM算法

目前常用的非线性最小二乘法有三种：最速下降法、高斯-牛顿法和LM法。其中LM方法中，迭代增量公式如下：

Δθ＝(J^TJ+μI)^-1J^Tε_i

式中，μ为阻尼因子，J为偏导数矩阵，ε_i为估计误差。其中，μ的取值很大程度决定了LM方法的收敛速度以及求解的精确度。然而，实际应用中，μ的最优取值无法确定。为改善LM算法的收敛速度以及求解的精确度，本发明提出一种改进因子，此因子可随迭代次数变化而变化，使得迭代初期可放宽对初始值，而后期则可加大收敛速度，从而使得LM算法在迭代初期具有最速下降法，而在最优解附近又具有高斯-牛顿法的优势，此下降因子为：

1-(1-μ)ⁱ⁺¹

从而增量迭代公式可重新表示为：

Δθ＝(J^TJ+(1-(1-μ)ⁱ⁺¹)I)^-1J^Tε_i

本发明的效果可通过以下仿真进一步说明：

仿真条件：

基于以上标定算法，利用采集的图像及标定图案的提取和自动编号的结果对摄像机进行标定。最后得到的标定结果的标准偏差为[0.04356,0.06792]，单位为像素。标定结果的误差分布如图9所示，由图可知标定的结果稳定可靠，误差分布均匀。

综上所述，本发明针对现有标定算法所采用畸变模型较为简单，从而不能实现高精度标定的问题，提出了基于混合畸变模型的摄像机标定方法。此方法首先基于小孔成像建立成像模型，而后针对三种畸变以及平面不均等物理误差建立了混合畸变模型。为标定摄像机内外参数以及畸变参数，本发明确定了可实现自动定位的标定模板，而后基于此模板实现标志点的亚像素定位。在获得标志点的精确定位条件下，本发明首先通过求解超定方程组得到单应矩阵的初始解，而后利用所提出的改进LM方法进行高效优化；基于此单应矩阵，可求得摄像机的内外参数，并可利用最大似然估计实现进一步优化；在求得内外参数后，通过求解超定方程而后利用改进的LM方法进一步优化，可分别获得摄像机的三种畸变参数以及有限元畸变参数。与传统的标定方法相比，该本发明全面考虑了各种畸变，并基于提出的准则对伪特征点进行滤除，而后基于提出的圆心自动编号方法实现标志点的自动定位，此外，在求解标定参数的过程中，基于改进的LM算法进行优化，从而显著加快了收敛速度并提高了标定精度。基于以上讨论可知，本发明所提方法可为机器视觉、三维测量等工程应用中高精度相机标定提供坚实的理论与实现依据。

以上所述，仅为本发明创造较佳的具体实施方式，但本发明创造的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明创造披露的技术范围内，根据本发明创造的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明创造的保护范围之内。

Claims (6)

1.一种基于混合畸变模型的高精度摄像机标定方法，包括如下步骤：

S1.建立摄像机畸变模型；

S2.标定图案选取、标定特征提取及特征点自动编号；

S3.摄像机标定。

2.如权利要求1所述的基于混合畸变模型的高精度摄像机标定方法，所述建立摄像机畸变模型，包括

1)小孔成像模型：在该模型中建立四个基本坐标系，包括：世界坐标系O_wX_wY_wZ_w，相机坐标系O_cX_cY_cZ_c，图像坐标系xy，像素坐标系uv，摄像机坐标系(X_c,Y_c,Z_c)以摄像机镜头光心O_c为坐标原点，X_c,Y_c轴平行于图像平面，Z_c轴垂直于图像平面，于图像平面的交点在图像坐标系上的坐标为(u₀,v₀)，其摄像机的主点，摄像机镜头光心与主点之间的距离就是焦距f；

基于上述坐标系，可得，空间一点P(X_w,Y_w,Z_w)的位置在摄像机坐标系下可表示为：

$\left[\begin{array}{c}{X}_c\\ Y_c\\ Z_c\end{array}\right]=\left[\begin{array}{cc}R& T\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\end{array}\right]$

其中为3×3正交旋转矩阵，为3×1平移矩阵；

基于小孔成像模型，得摄像机坐标(X_c,Y_c,Z_c,1)^T在像素坐标系下可表示为：

$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}{\mathrm{\ρ}}_x& {\mathrm{\α\ρ}}_x& u_0\\ 0& {\mathrm{\ρ}}_y& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{X}_c\\ Y_c\\ Z_c\\ 1\end{array}\right]$

其中，(u,v,1)^T为图像像素坐标，ρ_x,ρ_y分别为像平面中x,y方向上单位距离的像素数，(u₀,v₀)为摄像机主点坐标，其是像平面中心的像素坐标，α为描述像素坐标轴垂直误差的偏斜因子；

由以上两式，得图像坐标与世界坐标的关系式：

$Z_c\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}{f}_{\mathrm{\α}}& {\mathrm{\αf}}_{cx}& u_0\\ 0& f_{cy}& v_0\\ 0& 0& 1\end{array}\right]\[\begin{array}{cc}R& T\end{array}\]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]=K_{c}E\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]$

其中，为摄像机内参数矩阵，E＝[R T]为外参数矩阵；

2)构建摄像机畸变模型：采用基于矩形有限元方法的混合畸变模型，此模型描述如下：

$\left\{\begin{array}{c}dx={\mathrm{dx}}_r+{\mathrm{dx}}_d+{\mathrm{dx}}_p+{\mathrm{dx}}_u\\ dy={\mathrm{dy}}_r+{\mathrm{dy}}_d+{\mathrm{dy}}_p+{\mathrm{dy}}_u\end{array}\right.$

其中，dx,dy分别为像点x坐标和y坐标的偏差；

$\left\{\begin{array}{c}d{x}_r=x_i\[K_{c1}{R}^2+K_{c2}{R}^4+K_{c3}{R}^6\]\\ d{y}_r=y_i\[K_{c1}{R}^2+K_{c2}{R}^4+K_{c3}{R}^6\]\end{array}\right.$

其中，x_i,y_i分别为符合小孔成像的理想无畸变像平面坐标，K_c1,K_c2,K_c3为三阶径向畸变系数；

$\left\{\begin{array}{c}{\mathrm{dx}}_d=P_{c1}(R^2+2x_i^2)+2P_{c2}{x}_i{y}_i\\ {\mathrm{dy}}_d=P_{c2}(R^2+2y_i^2)+2P_{c1}{x}_i{y}_i\end{array}\right.$

其中，P_c1,P_c2为二阶切向畸变系数；

$\left\{\begin{array}{c}{\mathrm{dx}}_p=S_{c1}(x_i^2+y_i^2)=S_{c1}{R}^2\\ {\mathrm{dy}}_p=S_{c2}(x_i^2+y_i^2)=S_{c2}{R}^2\end{array}\right.$

其中，S_c1,S_c2为二阶薄棱镜畸变系数；

$\left\{\begin{array}{c}{\mathrm{dx}}_u=(1-\mathrm{\δ}x)(1-\mathrm{\δ}y){\mathrm{dx}}_{i,j}+\mathrm{\δ}x(1-\mathrm{\δ}y){\mathrm{dx}}_{i,j+1}+(1-\mathrm{\δ}x){\mathrm{\δydx}}_{i+1,j}+{\mathrm{\δx\δydx}}_{i+1,j+1}\\ {\mathrm{dy}}_u=(1-\mathrm{\δ}x)(1-\mathrm{\δ}y){\mathrm{dx}}_{i,j}+\mathrm{\δ}x(1-\mathrm{\δ}y){\mathrm{dx}}_{i,j+1}+(1-\mathrm{\δ}x){\mathrm{\δydx}}_{i+1,j}+{\mathrm{\δx\δydx}}_{i+1,j+1}\end{array}\right.$

其中，δx,δy分别为像点在矩形网格x,y方向上的长度比。

3.如权利要求1所述的基于混合畸变模型的高精度摄像机标定方法，其特征在于，所述步骤S2中标定图案选取，其采用的标定板是铝制标定板，使用其对摄像机进行参数标定，标定板上有9×11个圆形标定点，其中小圆直径为4mm，大圆直径为8mm，五个大圆用于在标定过程中标识标定板方向。

4.如权利要求1所述的基于混合畸变模型的高精度摄像机标定方法，其特征在于，标定特征提取包括以下步骤：图像采集、图像灰度化、图像滤波、图像二值化、轮廓提取、伪点滤除、椭圆拟合、圆心提取；

利用OpenCV函数库里的函数直接将彩色图像读入变为灰色图像，紧接着对获取的灰色图像进行滤波，使用函数cvSmooth实现图像滤波，图像二值化中，阈值经试验可设为100，轮廓提取采用基于迭代的改进Harris角点检测法，此方法描述如下：考虑角点q邻域里的任一点p，▽I(p)表示p点处的图像梯度矢量，则由梯度性质可得：

▽I(p)(q-p)＝0

考虑角点q邻域里所有p点，则构成一个超定方程组，从而利用最小均方误差(MMSE)准则进行求解；

提取轮廓后，利用如下规则消除伪特征点：

面积准则：提取特征对象的面积根据***配置进行范围确定：

S_min＜S＜S_max

圆度准则：用下式衡量形状与圆的接近程度：

$C=\frac{T^2}{S}$

其中T为圆周长，

误差准则：将边缘轮廊点到拟合椭圆的代数距离作为误差，同时采用以下准则进一步滤除误差较大的圆形特征点：

$\left\{\begin{array}{c}\mathrm{\&epsiv;}<{\mathrm{\&epsiv;}}_{\max }\\ \stackrel{\&OverBar;}{\mathrm{\&epsiv;}}<{\mathrm{\&epsiv;}}_{mean}\end{array}\right.$

其中表示轮廓点的平均误差，ε_max表示最大允许误差，ε_mean表示允许平均误差；

伪特征点消除后，采用加权最小二乘法进行轮廓椭圆拟合提取亚像素级精度的圆心坐标。

5.如权利要求1所述的基于混合畸变模型的高精度摄像机标定方法，其特征在于，所述步骤S2中特征点自动编号：利用标定板上的5个实心圆，每个实心圆都对应一个确定不变的编号，进行如下步骤：

(1)先用鼠标选取空间某矩形区域在图像上对应的四个角点，获取四个角点的图像坐标并经亚像素处理；

(2)对选定区域的圆进行自动记数，并根据计数结果生成网格图。

6.如权利要求2所述的基于混合畸变模型的高精度摄像机标定方法，其特征在于，所述摄像机标定的方法如下：

1)单应性矩阵求解

假设模板位于世界坐标系Z＝0的位置，由针孔模型，得

sm＝HM

其中，s为尺度系数，m,M分别为图像平面以及对应的模板上的点，H＝K_c[R T]；

基于平面模板各角点空间坐标，通过求解如下的MMSE估计，得到H的估计值，

$\hat{H}=\arg cmin\underset{i}{\&Sigma;}\left|\right|m_i-{\hat{m}}_i|{|}^2$

其中， 为H矩阵的第i个行向量；

采用改进的Leverberg-Marquardt法求解，用非线性优化方法求解需要有一个用来迭代的合适初始值，通过求解下面方程获得：

$\left[\begin{array}{ccc}{\hat{M}}^T& 0^T& -u{\hat{M}}^T\\ 0^T& {\hat{M}}^T& -u{\hat{M}}^T\end{array}\right]\left[\begin{array}{c}{h}_1^T\\ h_2^T\\ h_2^T\end{array}\right]=0$

给定n＞6对点，就形成超定方程组，则通过奇异值分解(SVD)进行求解；

2)内外参数求解

基于内参数矩阵R为正交矩阵的性质，通过求解如下超定方程组得到摄像机的内外参数：

Vb＝0

此方程组为n≥3次观察得到的如下方程的叠加：

$\left[\begin{array}{c}{v}_{12}^T\\ v_{11}^T-v_{22}^T\end{array}\right]b=0$

式中v_ij＝[h_1ih_1j,h_1ih_2j+h_2ih_1j,h_2ih_2j,h_3ih_1j+h_1ih_3j,h_3ih_2j+h_2ih_3j,h_3ih_3j]^T，b＝[B₁₁,B₁₂,B₂₂,B₁₃,B₂₃,B₃₃]^T，

而后，通过最大似然估计对此解进行精细化，以得到更好的优化解，似然估计问题表示如下：

$\underset{i=1}{\overset{n}{\&Sigma;}}\underset{j=1}{\overset{m_i}{\&Sigma;}}\left|\right|m_{ij}-\hat{m}(K_c,R_i,T_i,M_j)|{|}^2$

式中，m_ij为第i幅图像上第j个像点；为空间点M_j在第i幅图像的投影点；

3)畸变参数的求解；

4)有限元模型畸变参数求解；

5)改进LM算法。