一种双目立体视觉相机标定参数的优化方法

1.本发明属于双目立体视觉标技术领域，更具体的讲，涉及一种双目立体视觉相机标定参数的优化方法。


背景技术：

2.随着经济和科技的迅速发展，智能制造技术越来越受到重视，其中，体现人工智能的双目立体视觉技术显得尤为重要。基于双目视觉的三维重建可以恢复场景的三维信息，在智能生产、交通监控、机器导航、航天航空、医学建模、文物复原以及非接触高精度测量等方面，能够有效地提高生产效率和产品质量，降低运营成本和资源浪费。
3.相机标定在双目立体视觉领域有着非常重要的作用，是实现三维场景重建的首要步骤，并且相机标定的精度直接影响重建的结果。相机标定就是由已知空间位置的点与其图像对应点求解相机内外参数的过程。其中，待优化的相机标定参数如下：
4.左相机内参：
5.左相机畸变系数：d
l
＝[k
1l k
2l p
1l p
2l k
3l
]；
[0006]
右相机内参：
[0007]
右相机畸变系数：dr＝[k
1r k
2r p
1r p
2r k
3r
]：
[0008]
左右相机之间的旋转参数：以及平移向量：t＝[t
1 t
2 t3]
t
。
[0009]
其中，f
xl
、f
yl
分别为x、y方向上的1个像素在左相机感光板上的物理长度，u
0l
、v
0l
分别为左图像的中心像素坐标和图像原点像素坐标之间相差的横向和纵向像素数，k
1l
、k
2l
、k
3l
为左相机的径向畸变系数，p
1l
、p
2l
为左相机的切向畸变系数；同样地，f
xr
、f
yr
分别为x、y方向上的1个像素在右相机感光板上的物理长度，u
0r
、v
0r
分别为右图像的中心像素坐标和图像原点像素坐标之间相差的横向和纵向像素数，k
1r
、k
2r
、k
3r
为右相机的径向畸变系数，p
1r
、p
2r
为右相机的切向畸变系数；r1,...,r9均为x、y、z轴旋转角的函数，t1、t2、t3分别为x、y、z方向的平移尺寸。
[0010]
双目相机标定过程坐标转换如图1所示，标定系统中包括了世界坐标系、相机坐标系、图像坐标系和像素坐标系，在标定的过程中，标定坐标通过由相机内外参数构成的变换矩阵转换到不同的坐标系下。其中，张正友相机标定法是最常用的一种相机标定方法，鲁棒性和精度都很高，但此方法只考虑了2d像素约束且进行参数估计时引入了一定的假设，使得解析参数存在一定的误差，其参数精度仍有提高的空间。


技术实现要素：

[0011]
本发明的目的在于克服现有技术的不足，提供一种双目立体视觉相机标定参数的优化方法，在已知标定后的相机内外参数的条件下，综合多方面约束，对标定后的相机标定参数进行优化，以提高相机标定参数的标定精度，强化相机标定结果的准确性和实用性。
[0012]
为实现上述发明目的，本发明双目立体视觉相机标定参数的优化方法，其特征在于，包括：
[0013]
(1)、获取b幅标定板图像，并根据b幅标定板图像依据相机标定方法获得相机标定参数初始值并作为待优化的相机标定参数，其中，标定板角点规格中，每一行的角点个数为n，每一列的角点个数为m；
[0014]
(2)、构建2d像素约束
[0015]
2.1)、对于第q,q＝1,2,...,b幅标定板图像的第i,i＝1,2,...m行第j,j＝1,2,...n列个角点qp
ij
，其世界坐标系空间坐标为qp
ij
＝[qx
ij
,
qyij
,
qzij
]
t
，qx
ij
、
qyij
、
qzij
分别为角点qp
ij
的x、y、z轴坐标值，直接对第q幅标定板图像进行角点提取后，得到的角点qp
ij
在左右相机图像上角点的坐标位置分别为的坐标位置分别为的坐标位置分别为分别为角点的x、y轴的坐标值，分别为角点的x、y轴的坐标值；
[0016]
根据针孔相机模型，分别通过左、右相机标定参数将角点qp
ij
投影到左、右相机图像平面上，得到左相机2d投影点和右相机投影点和右相机投影点分别为投影点的x、y轴的坐标值，分别为投影点的x、y轴的坐标值；
[0017]
将左相机坐标系作为世界坐标系，则角点qp
ij
在左相机上的参数转换关系式为：
[0018][0019]
在右相机上的参数转换关系式为：
[0020][0021]
其中，w
l
和wr为未知的尺度因子；
[0022]
2.2)、将所有的投影点的坐标与直接对图像进行角点提取的坐标进行比较，计算其重投影误差，构成2d像素约束，即：
[0023][0024]
(3)、构建3d重建约束
[0025]
3.1)、通过左相机外参t
lwc
＝[0 0 0]
t
将世界坐标系下的标定板角点的坐标点集ωw＝{qp
ij
|[qx
ij
,
qyij
,
qzij
]
t
}，q＝1,2,...,b，i＝1,2,...m，j＝1,2,...n变换到左相机坐标系下，得到一系列在左相机坐标系下的理想标定板角点三维坐标点集
分别为点的x、y、z轴坐标值，坐标转换关系式如下：
[0026][0027]
3.2)、通过相机畸变系数对左、右图像中提取的标定板角点的像素坐标进行去畸变，记理想无畸变的左图像像素坐标为和右图像像素坐标为和右图像像素坐标为分别为像素点的x、y轴坐标值，分别为像素点的x、y轴坐标值，则有左图像的径向畸变矫正公式如下：
[0028][0029]
分别为左图像上理想角点径向畸变后的x、y轴坐标值；
[0030]
右图像的径向畸变矫正公式如下：
[0031][0032]
分别为右图像上理想角点径向畸变后的x、y轴坐标值；
[0033]
其中，均为图像像素点到图像中心点的距离；
[0034]
左图像的切向畸变矫正公式如下：
[0035][0036]
分别为左图像上角点切向畸变后的x、y轴坐标值；
[0037]
右图像的切向畸变矫正公式如下：
[0038][0039]
分别为右图像上角点切向畸变后的x、y轴坐标值；
[0040]
结合左右相机的径向畸变和切向畸变矫正公式，有如下关系式：
[0041][0042]
反解上式得到最终理想无畸变的左图像像素坐标和右图像像素坐标
[0043]
3.3)、通过直接线性变换法计算得到重建的三维点坐标，记左、右相机的投影矩阵分别为p
l
、pr，其表示公式如下：
[0044][0045]
其中，p
l1
,p
l2
,p
l3
为p
l
的行向量，p
r1
,p
r2
,p
r3
为pr的行向量；
[0046]
根据直接线性变换法有以下等式：
[0047][0048]
使用奇异值分解方法解上式得到重建的三维点，记所有重建点的坐标集为ω＝{qp
ij
|[qx
ij
，
qyij
，
qzij
]
t
}，qx
ij
、
qyij
、
qzij
分别为重建点qp
ij
的x、y、z轴坐标值；
[0049]
3.4)、计算当前残差，得到构造3d重建约束为：
[0050][0051]
(4)、构建3d点的共面约束
[0052]
4.1)、第q幅标定板图像中理想标定板角点所在平面的法向量用相机外参的旋转矩阵r的第三列表示，即向量
[0053]
4.2)、第q幅标定板图像重建出的三维角点点集ωq＝{qpi′
|[qxi′
,qyi′
,qzi′
]
t
}中，记从角点qpi′
，i
′
＝1,2,...,m
×
n出发到角点qpj′
，j
′
＝1,2,...,m
×
n,j
′
≠i
′
的向量的单位方向向量为即：
[0054][0055]
4.3)、将向量和单位方向向量内积的绝对值求和，作为衡量面上向量与面的法向量垂直程度的误差值，构成共面的约束项：
[0056]
[0057]
(5)、构建3d点的共线约束
[0058]
5.1)、对于第q幅标定板图像，找到其重建的三维角点点集ωq＝{qpi′
|[qxi′
,qyi′
,qzi′
]
t
}中位于标定板在空间中第r行的首端点qp
rs
和尾端点qp
re
，将首尾端点连接得共线向量即
[0059]
5.2)、对于第r行的第k个中间点qp
rk
，同理，将该中间点qp
rk
与首端点qp
rs
连接得到第二个向量
[0060]
5.3)、接着，计算向量在第r行的共线向量上的投影点，计算公式为：其中，是该行的共线向量的单位方向向量；
[0061]
5.4)、然后，计算中间点qp
rk
与其投影点间的距离，即中间点qp
rk
与第r行共线向量的距离，作为重建的共线误差中的一个残差项，对所有行与列，和所有点进行相同的处理，计算得到最终的共线误差项e
line
：
[0062][0063]
其中，qpk′c、分别为第q幅标定板图像在空间中第c列的第k
′
个中间点及其在该列共线向量上的投影点；
[0064]
(6)、构建3d点的径向约束
[0065]
6.1)、第q幅标定板图像的三维角点点集ωq＝{qpi′
|[qxi′
,qyi′
,qzi′
]
t
}，其编号为从左到右、从上到下，选取四个角点qp1、qpn、qpm×n、qp
(m-1)
×
n+1
；
[0066]
6.2)、则径向约束的两项误差为：
[0067][0068]
其中，radius为理想标定板径向长度；
[0069]
(7)、构建3d点的角度约束
[0070]
7.1)、将步骤6.1)选出的四个角点qp1、qpn、qpm×n、qp
(m-1)
×
n-1
分别记为点a、b、c、d，并构成一个四边形，将线段ab、ad、bc、dc的方向向量分别写为有：
[0071][0072]
7.2)、对该四边形的四个角度进行计算，并将每一个角度值与直角进行残差计算，对于所有的标定图像，该项角度误差可以写为：
[0073]
[0074]
(8)、考虑上述所有6个约束项，得到一个新约束即优化目标函数为：
[0075]
min e
total
＝e
pix
+e
pnt
+e
pl
+e
line
+e
diag
+e
ang
[0076]
将待优化的相机标定参数写为如下形式：
[0077][0078]
将优化目标函数简写为：
[0079][0080]
式中，θ表示上述待优化的相机标定参数的集合{a
l
,ar,d
l
,dr,r,t}，e是所有约束项函数平方根的集合
[0081]
(9)、初始化相机标定参数变量的初始值θ
←
θ0，采用列文伯格-马夸尔特方法，依据优化目标函数对待优化的相机标定参数进行最优化，得到相机标定参数的最优解θ
*
，其中，θ0为标定后相机标定参数。
[0082]
本发明的发明目的是这样实现的：
[0083]
要提高相机标定参数精度，还需要使得标定板三维重建后的结果与实际三维场景保持一致，即重建点的坐标应与实际空间坐标一致且点与点之间的相对位置不变，同时标定板在空间中的几何特性应经过标定重建后保持不变。综合考虑上述情况，本发明双目立体视觉相机标定参数的优化方法对不同情况构造了相应的几何约束，并将所有约束项进行叠加构成最终的优化目标函数，最后采用列文伯格-马夸尔特方法对目标函数进行最优化，得到优化后的相机标定参数。本发明提出的相机标定参数优化算法相较于张正友相机标定法有更高的准确度，使用优化后的参数进行标定后，重投影误差和三维重建测量结果误差大大降低。
附图说明
[0084]
图1是双目相机标定过程坐标转换的示意图；
[0085]
图2是本发明双目立体视觉相机标定参数的优化方法一种具体实施方式的流程图；
[0086]
图3是共面约束示意图；
[0087]
图4是共线表征示意图；
[0088]
图5是径向约束示意图；
[0089]
图6是双目捕获图像：(a)左相机；(b)右相机；
[0090]
图7是迭代过程中，左、右相机的部分内参数变化曲线；
[0091]
图8是为迭代过程中，外参数的变化曲线，由于旋转矩阵的参数较多，为了便于显
示本发明将其变换为四元数的表示；
[0092]
图9是优化算法中的参数变化曲线，其中，(a)为增益比ρ变化曲线，(b)为修正因子ε变化曲线，(c)为步长调整因子v变化曲线，(d)为梯度值||j(θu)||变化曲线；
[0093]
图10是迭代过程中目标函数的变化曲线；
[0094]
图11是各个单幅标定位姿的重投影误差对比图，横坐标表示标定位姿，纵坐标表示重投影误差，下方曲线为本发明的单幅标定位姿的重投影误差，上方曲线为张正友标定法的单幅标定位姿的重投影误差；
[0095]
图12为重投影点的坐标误差；图(a)-(b)为传统法的结果；图(c)-(d)为本发明算法的结果；图(a)、(c)横坐标轴表示的是重投影点与理想角点像素坐标在横坐标上的差值的绝对值，纵坐标表示的是重投影点与理想角点像素坐标在纵坐标轴上的差值的绝对值；图(b)、(d)横坐标表示重投影点坐标值跟理想坐标的差值，纵坐标为重投影点数；
[0096]
图13是重建点云一实例图。
[0097]
具体实施方式
[0098]
下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。
[0099]
图1是本发明双目立体视觉相机标定参数的优化方法一种具体实施方式的流程图。
[0100]
在本实施例中，如图2所示，本发明双目立体视觉相机标定参数的优化方法包括以下步骤：
[0101]
步骤s1：获取标定板图像以及标定后相机标定参数初始值
[0102]
获取b幅标定板图像，并根据b幅标定板图像依据相机标定方法获得相机标定参数初始值并作为待优化的相机标定参数，其中，标定板角点规格中，每一行的角点个数为n，每一列的角点个数为m。在具体实施过程中，相机标定方法可以采用张正友相机标定法以获得相机标定参数初始值。
[0103]
步骤s2：构建2d像素约束
[0104]
步骤s2.1：对于第q,q＝1,2,...,b幅标定板图像的第i,i＝1,2,...m行第j,j＝1,2,...n列个角点qp
ij
，其世界坐标系空间坐标为qp
ij
＝[qx
ij
,
qyij
,
qzij
]
t
，qx
ij
、
qyij
、
qzij
分别为角点qp
ij
的x、y、z轴坐标值，直接对第q幅标定板图像进行角点提取后，得到的角点qp
ij
在左右相机图像上角点的坐标位置分别为的坐标位置分别为分别为角点的x、y轴的坐标值，的x、y轴的坐标值，分别为角点的x、y轴的坐标值。
[0105]
根据针孔相机模型，分别通过左、右相机标定参数将角点qp
ij
投影到左、右相机图像平面上，得到左相机2d投影点和右相机投影点和右相机投影点分别为投影点的x、y轴的坐标值，分别为投影点的x、y轴的坐标值。
[0106]
将左相机坐标系作为世界坐标系，则角点qp
ij
在左相机上的参数转换关系式为：
[0107][0108]
在右相机上的参数转换关系式为：
[0109][0110]
其中，w
l
和wr为未知的尺度因子。
[0111]
步骤s2.2：将所有的投影点的坐标与直接对图像进行角点提取的坐标进行比较，计算其重投影误差，构成2d像素约束，即：
[0112][0113]
步骤s3：构建3d重建约束
[0114]
步骤s3.1：通过左相机外参t
lwc
＝[0 0 0]
t
将世界坐标系下的标定板角点的坐标点集ωw＝{qp
ij
|[qx
ij
,
qyij
,
qzij
]
t
}，q＝1,2,...,b，i＝1,2,...m，j＝1,2,...n变换到左相机坐标系下，得到一系列在左相机坐标系下的理想标定板角点三维坐标点集变换到左相机坐标系下，得到一系列在左相机坐标系下的理想标定板角点三维坐标点集分别为点的x、y、z轴坐标值，坐标转换关系式如下：
[0115][0116]
步骤s3.2：通过相机畸变系数对左、右图像中提取的标定板角点的像素坐标进行去畸变，记理想无畸变的左图像像素坐标为和右图像像素坐标为和右图像像素坐标为分别为像素点的x、y轴坐标值，分别为像素点的x、y轴坐标值，则有左图像的径向畸变矫正公式如下：
[0117][0118]
分别为左图像上理想角点径向畸变后的x、y轴坐标值；
[0119]
右图像的径向畸变矫正公式如下：
[0120][0121]
分别为右图像上理想角点径向畸变后的x、y轴坐标值；
[0122]
其中，均为图像像素点到图像中心点的距离；
[0123]
左图像的切向畸变矫正公式如下：
[0124][0125]
分别为左图像上角点切向畸变后的x、y轴坐标值；
[0126]
右图像的切向畸变矫正公式如下：
[0127][0128]
分别为右图像上角点切向畸变后的x、y轴坐标值；
[0129]
结合左右相机的径向畸变和切向畸变矫正公式，有如下关系式：
[0130][0131]
反解上式得到最终理想无畸变的左图像像素坐标和右图像像素坐标
[0132]
步骤s3.3：通过直接线性变换法计算得到重建的三维点坐标，记左、右相机的投影矩阵分别为p
l
、pr，其表示公式如下：
[0133][0134]
其中，p
l1
,p
l2
,p
l3
为p
l
的行向量，p
r1
,p
r2
,p
r3
为pr的行向量；
[0135]
根据直接线性变换法有以下等式：
[0136][0137]
使用奇异值分解方法解上式得到重建的三维点，记所有重建点的坐标集为ω＝
{qp
ij
|[qx
ij
,
qyij
,
qzij
]
t
}，qx
ij
、
qyij
、
qzij
分别为重建点qp
ij
的x、y、z轴坐标值。
[0138]
步骤s3.4：计算当前残差，得到构造3d重建约束为：
[0139][0140]
步骤s4：构建3d点的共面约束
[0141]
步骤s4.1：第q幅标定板图像中理想标定板角点所在平面的法向量用相机外参的旋转矩阵r的第三列表示，即向量
[0142]
步骤s4.2：第q幅标定板图像重建出的三维角点点集ωq＝{qpi′
|[qxi′
,qyi′
,qzi′
]
t
}中，记从角点qpi′
，i
′
＝1,2,...,m
×
n出发到角点qpj′
，j
′
＝1,2,...,m
×
n,j
′
≠i
′
的向量的单位方向向量为即：
[0143][0144]
步骤s4.3：如图3所示，理想情况下，面的法向量与面上任意两点间的向量因为垂直，而内积应该等于0，在误差的情况下该值将非0，因此，将向量和单位方向向量内积的绝对值求和，作为衡量面上向量与面的法向量垂直程度的误差值，构成共面的约束项：
[0145][0146]
步骤s5：构建3d点的共线约束
[0147]
如图4所示，理想的重建标定板三维点在空间中应保持同一行的点在同一条直线上，而同一列的点也在同一条直线上，该特性在经过刚体变换和投影变换与逆变换之后不应发生变化。但由于标定参数的不准确，该项将被引入误差，导致原本在同一列或同一行上的点不在严格满足共线约束，如图4左图所示。反过来共线的程度可以表征该标定参数的准确性，而作为标定参数优化的一项约束。为了表征本应在同一条线上的重建标定板角点的共线程度，本发明构建3d点的共线约束模型步骤如下：
[0148]
步骤s5.1：对于第q幅标定板图像，找到其重建的三维角点点集ωq＝{qpi′
|[qxi′
,qyi′
,qzi′
]
t
}中位于标定板在空间中第r行的首端点qp
rs
和尾端点qp
re
，将首尾端点连接得共线向量即
[0149]
步骤s5.2：对于第r行的第k个中间点qp
rk
，同理，将该中间点qp
rk
与首端点qp
rs
连接得到第二个向量
[0150]
步骤s5.3：接着，计算向量在第r行的共线向量上的投影点，计算公式为：其中，是该行的共线向量的单位方向向量。
[0151]
步骤s5.4：然后，计算中间点qp
rk
与其投影点间的距离，即中间点qp
rk
与第r行共线向量的距离，作为重建的共线误差中的一个残差项，对所有行与列，和所有点进行相同的处理，计算得到最终的共线误差项e
line
：
[0152][0153]
其中，qpk′c、分别为第q幅标定板图像在空间中第c列的第k
′
个中间点及其在该列共线向量上的投影点。
[0154]
步骤s6：构建3d点的径向约束
[0155]
重建的标定板角点还应满足如图5所示的径向约束。重建点阵集的点与点之间的最远距离应保持一致。为方便表示，将标定板的左上角标记为a、右上角b、右下角c以及左下角d。ad与bc的距离应与理想值保持一致。
[0156]
步骤s6.1：第q幅标定板图像的三维角点点集ωq＝{qpi′
|[qxi′
,qyi′
,qzi′
]
t
}，其编号为从左到右、从上到下，选取四个角点qp1、qpn、qpm×n、qp
(m-1)
×
n+1
，即图4中，左上角a、右上角b、右下角c以及左下角d。
[0157]
步骤s6.2：则径向约束的两项误差为：
[0158][0159]
其中，radius为理想标定板径向长度。
[0160]
步骤s7：构建3d点的角度约束
[0161]
步骤s7.1：将步骤s6.1选出的四个角点qp1、qpn、qpm×n、qp
(m-1)
×
n-1
分别记为点a、b、c、d，并构成一个四边形，将线段ab、ad、bc、dc的方向向量分别写为有：
[0162][0163]
步骤s7.2：对该四边形的四个角度进行计算，并将每一个角度值与直角进行残差计算，对于所有的标定图像，该项角度误差可以写为：
[0164][0165]
步骤s8：构建优化目标函数
[0166]
考虑上述所有6个约束项，得到一个新约束即优化目标函数为：
[0167]
min e
total
＝e
pix
+e
pnt
+e
pl
+e
line
+e
diag
+e
ang
[0168]
将待优化的相机标定参数写为如下形式：
[0169][0170]
将优化目标函数简写为：
[0171][0172]
式中，θ表示上述待优化的相机标定参数的集合{a
l
,ar,d
l
,dr,r,t}，e是所有约束项函数平方根的集合
[0173]
步骤s9：以标定后相机标定参数作为初值，采用列文伯格-马夸尔特方法，依据优化目标函数进行优化。
[0174]
初始化相机标定参数变量的初始值θ
←
θ0，采用列文伯格-马夸尔特方法，依据优化目标函数对待优化的相机标定参数进行最优化，得到相机标定参数的最优解θ
*
，其中，θ0为标定后相机标定参数。
[0175]
由于相机的畸变模型是非线性的，因此，在偏微分后，可采用列文伯格-马夸尔特(l-m，levenberg-marquardt)方法对目标函数进行最优化。根据l-m法的原理，本发明相机标定参数在优化过程的更新规则具有以下的形式：
[0176][0177]
式中，u表示的是迭代次数，εu为第u次迭代的步长自适应参数，εui是正定修正项，每次迭代时进行调整，以确保减少残差。当εu值较大时，收敛方式与梯度下降法类似，迭代次数减少，θu快速下降至极小值附近；εu值较小乃至接近于0时，与高斯牛顿法类似，步长缩小，稳定收敛。εu的调整需要引入阻尼调整因子ρu：
[0178][0179]
ρu表示近似模型与实际模型函数下降量的差异，分子为实际模型的下降量，分母为近似模型下降量。ρu值越大，越接近于1，证明近似模型可靠，收敛至极小值附近，可减小εu值，进行类似于高斯牛顿收敛；如果ρu值较小时，说明还处于类似梯度下降阶段，距离极值较远，可增大εu值；而ρu值为负时，意味着收敛方向相反，需要增大εu值使ρu变为正值，纠正收敛方向。j(θu)为本系统的雅可比矩阵，可写为：
[0180][0181]
记标定参数的初始估计值θ0∈r
30
，优化目标函数项e＝{e1,...,e6}∈r
30
→
r是从30维实数空间到1维实数空间的映射，迭代步长自适应变化的参数u
max
,v
in
,τ，和终止条件的判断阈值σ1,σ2,σ3。接着，求优化算法系统的解，即θ
*
∈r
30
标定参数，使得目标函数最小。
[0182]
优化算法的详细步骤如下：
[0183]
步骤s9.1：初始化标定参数变量的初始值θ
←
θ0，迭代步长自适应变化的初始参数v
←vin
和迭代次数u
←
0；
[0184]
步骤s9.2：为避免死循环，对于第u次迭代，当相当于梯度的幅值j(θu)
t
e(θu)收敛到一个接近0的阈值时，说明最优化算法已到达极小值，则可根据下式计算终止条件的标志η：
[0185][0186][0187]
当η＝true时，算法的迭代过程将终止；
[0188]
步骤s9.3：判断迭代是否终止，终止条件为：η＝true或u≥u
max
，若成立，则跳转到步骤s9.10；否则，进入下一步骤s9.4；
[0189]
步骤s9.4：由最优化公式可知，lm法是高斯牛顿法的修正，为保证每次迭代产生的方向是下降方向，即保证系统的hessian矩阵(j(θu)
t
j(θ
new
)+εui)是正定的，引入了正定修正项因子εu。此时，由雅可比矩阵的转置和本身的乘积可以得到半正定矩阵，选取此半正定矩阵对角向量的最大值可计算得修正因子εu，计算过程如下：
[0190][0191]
则修正因子εu可以根据输入的修正参数τ进行赋值：
[0192][0193]
步骤s9.5：计算标定参数在当前第u次迭代的更新变化值δθu。根据正定修正因子，将步骤s9.4中计算获得的半正定矩阵修正成为正定矩阵，以保证总体解严格最小。并将其逆矩阵乘上步骤s9.2中的判断矩阵，获得更新变化值，其过程用公式表示为：
[0194][0195]
其中，δθu的维度为30，第一个元素δf
xl-u
：
[0196][0197]
同理可求其他元素的更新变化值；
[0198]
步骤s9.6：判断更新变化值是否满足迭代终止条件，当更新变化值过小时，希望迭代终止，因此，如果max(δf
xl-u
,δf
yl-u
,...,δt
3-u
)≤σ2·
max(f
xl-u
,f
yl-u
,...,t
3-u
)，那么s＝true，进入步骤s9.10；否则，进入步骤s9.7；
[0199]
步骤s9.7：根据变化值更新标定参数θ
new
←
θu+δθu，并且计算当前的增益比ρu，该参数用以判断当前步长的质量，并为自适应步长更新提供指标。ρu的计算方式如下：
[0200]
步骤s9.7.1：首先，根据步骤s9.2中的判断矩阵，计算：
[0201][0202]
步骤s9.7.2：然后，计算增益比：
[0203][0204]
步骤s9.8：由步骤s9.4可知j(θu)
t
j(θu)是半正定的，在计算(j(θu)
t
j(θu))-1
的过程中，j(θu)
t
j(θu)可能是奇异矩阵或病态矩，会导致参数增量不稳定甚至算法不收敛。因此，l-m方法在高斯-牛顿法的基础上引入了信赖域，可用信赖域限制根据增益比的值ρu更新信赖域区间。当ρu＞0时，表明模型的质量较差，需要减小信赖域，且步长的自适应调整因子设置为初始值v＝v
in
；否则，可以通过自适应步长调整因子v适当增加搜索半径，依次执行：ε
u+1
＝εu·
v；v＝v
in
·
v；
[0205]
步骤s9.9：判断终止条件是否满足，当||j(θu)
t
e(θu)||
∞
≤σ1或e2(θ
u+1
)≤σ3时，将终止标志η的值设置为true，进入步骤s9.10；否则，令迭代次数u＝u+1，并返回步骤s9.2进行新的循环；
[0206]
步骤s9.10：输出此时的标定参数值θ
*
＝θu作为最优解。
[0207]
实例
[0208]
在本实例中，对一个汇聚式双目立体视觉系统进行相机标定参数的优化。采用的标定板规格为11列x8行(即m＝11，n＝8)，板上共88个内角，每个内角点之间的距离为15.00mm。本实例使用的汇聚式双目摄像机的分辨率均为4096x3000pixel，实验过程中一共切换了10个标定板的位姿(b＝10)，因此左相机捕获的图像i
lk
和右相机捕获的图像i
rk
共有20张图像，如图6所示。
[0209]
初始的相机标定参数估计值即传统标定算法的相机标定参数结果表1所示。
[0210][0211]
表1
[0212]
使用上述标定结果对源数据进行重投影，得到平均重投影误差的值为1.61367pixel。
[0213]
根据本发明介绍的各项约束项的表达式建立整体的优化目标函数：
[0214][0215]
将表1中的标定结果作为初始参数值θ0输入到标定优化算法中，设置优化算法的最大迭代次数为u
max
＝50，初始阻尼系数值为v
in
＝2，τ＝0.8，终止阈值参数梯度容差σ1＝1.0e-10
，参数容差σ2＝1.0e-8
，函数容差σ3＝1.0e-6
，并将初始值代入到步骤(9)中进行循环迭代求解。
[0216]
本实例一共迭代了27次，结合图7、图8和图9可知，在开始的几次迭代中，目标函数值(如图7和图8所示)和梯度值||j(θu)||变化较大，因此信赖域ε也应该逐渐增加，由图9(a)可知，此时ρ值的最大值出现，步长自适应调整因子v仍为初始值。随着算法的逐渐迭代循环，目标函数的变化值逐渐减小，梯度值也逐渐趋向于0，在此过程中，信赖域逐渐减小而逼近最优解。
[0217]
目标函数的输出变化曲线图10所示，可以看到目标函数的值最终稳定在20附近，这是由于目标函数为0仅是理想情况下的假设，实际情况中存在一定的误差，导致目标函数无法精确为0，其值最终收敛稳定说明了本发明算法的可行性。本发明获得的优化后的相机标定参数如表2所示。
[0218][0219]
表2
[0220]
为了验证本发明双目立体视觉相机标定参数的优化方法的有效性，计算平均重投影误差为0.649198pixel，如图11的重投影误差对比图所示，该结果与传统标定算法的平均重投影误差相比，减小了59.77％。
[0221]
为更直观的比较本发明与传统张氏标定算法的优化结果，将每个标定板角点的重投影结果与理想值进行比较，得到实验结果如图12所示。由图12(a)可见，传统张氏标定算法的重投影结果将95％的点用一半径为r的圆包裹在内，坐标差值集中的区域面积用s＝πr2表示，在图12(c)中，本发明的重投影结果的95％的点可以仅用一个半径为r的扇形区域覆盖，该区域的面积为s＝0.25πr2＜s，本发明的输出结果与传统张氏相机标定参数方法相比，重投影点的坐标误差分布更为集中，且在原点(0,0）附近，统计直方图12(b)、12(d)也较为直观的印证了上述结论。
[0222]
相机标定参数的准确度需要进一步在三维重建中验证，本实例采用张氏标点算法的参数对三维角点进行重建，并将角点之间的距离取平均，作为最终的重建结果，值为：14.9984mm。再采用本发明进行优化输出相机标定参数进行角点重建，求得平均角点距离为：15.0006mm。已知源标定板的内焦点距离为：15.0000mm。因此，两种相机标定算法的重建误差分别为：-0.0016mm和0.0006mm，与张式相机标定算法相比，本发明双目立体视觉相机标定参数的优化方法的重建误差减少了62.50％。
[0223]
综上所述，与传统的张氏标定法相比，本发明提出的双目立体视觉标定参数优化算法具有更高的准确度，体现在实验数据下的重投影误差降低了59.77％，三维重建的测量结果误差降低了62.50％。
[0224]
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术邻域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术邻域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

技术特征：
1.一种双目立体视觉相机标定参数的优化方法，其特征在于，包括：(1)、获取b幅标定板图像，并根据b幅标定板图像依据相机标定方法获得相机标定参数初始值并作为待优化的相机标定参数，其中，标定板角点规格中，每一行的角点个数为n，每一列的角点个数为m；(2)、构建2d像素约束2.1)、对于第q,q＝1,2,
…
,b幅标定板图像的第i,i＝1,2,
…
m行第j,j＝1,2,
…
n列个角点
q
p
ij
，其世界坐标系空间坐标为
q
p
ij
＝[
q
x
ij,q
y
ij,q
z
ij
]
t
，
q
x
ij
、
q
y
ij
、
q
z
ij
分别为角点
q
p
ij
的x、y、z轴坐标值，直接对第q幅标定板图像进行角点提取后，得到的角点
q
p
ij
在左右相机图像上角点的坐标位置分别为的坐标位置分别为的坐标位置分别为分别为角点的x、y轴的坐标值，分别为角点的x、y轴的坐标值；根据针孔相机模型，分别通过左、右相机标定参数将角点
q
p
ij
投影到左、右相机图像平面上，得到左相机2d投影点和右相机投影点和右相机投影点分别为投影点的x、y轴的坐标值，分别为投影点的x、y轴的坐标值；将左相机坐标系作为世界坐标系，则角点
q
p
ij
在左相机上的参数转换关系式为：在右相机上的参数转换关系式为：其中，w
l
和w
r
为未知的尺度因子；2.2)、将所有的投影点的坐标与直接对图像进行角点提取的坐标进行比较，计算其重投影误差，构成2d像素约束，即：(3)、构建3d重建约束3.1)、通过左相机外参t
lwc
＝[0 0 0]
t
将世界坐标系下的标定板角点的坐标点集ω
w
＝{
q
p
ij
|[
q
x
ij,q
y
ij,q
z
ij
]
t
}，q＝1,2,
…
,b，i＝1,2,
…
m，j＝1,2,
…
n变换到左相机坐标系下，得到一系列在左相机坐标系下的理想标定板角点三维坐标点集机坐标系下，得到一系列在左相机坐标系下的理想标定板角点三维坐标点集分别为点的x、y、z轴坐标值，坐标转换关系式如下：
3.2)、通过相机畸变系数对左、右图像中提取的标定板角点的像素坐标进行去畸变，记理想无畸变的左图像像素坐标为和右图像像素坐标为和右图像像素坐标为分别为像素点的x、y轴坐标值，分别为像素点的x、y轴坐标值，则有左图像的径向畸变矫正公式如下：如下：分别为左图像上理想角点径向畸变后的x、y轴坐标值；右图像的径向畸变矫正公式如下：右图像的径向畸变矫正公式如下：分别为右图像上理想角点径向畸变后的x、y轴坐标值；其中，均为图像像素点到图像中心点的距离；左图像的切向畸变矫正公式如下：左图像的切向畸变矫正公式如下：分别为左图像上角点切向畸变后的x、y轴坐标值；右图像的切向畸变矫正公式如下：右图像的切向畸变矫正公式如下：分别为右图像上角点切向畸变后的x、y轴坐标值；结合左右相机的径向畸变和切向畸变矫正公式，有如下关系式：
反解上式得到最终理想无畸变的左图像像素坐标和右图像像素坐标3.3)、通过直接线性变换法计算得到重建的三维点坐标，记左、右相机的投影矩阵分别为p
l
、p
r
，其表示公式如下：其中，p
l1
,p
l2
,p
l3
为p
l
的行向量，p
r1
,p
r2
,p
r3
为p
r
的行向量；根据直接线性变换法有以下等式：使用奇异值分解方法解上式得到重建的三维点，记所有重建点的坐标集为ω＝{
q
p
ij
|[
q
x
ij,q
y
ij,q
z
ij
]
t
}，
q
x
ij
、
q
y
ij
、
q
z
ij
分别为重建点
q
p
ij
的x、y、z轴坐标值；3.4)、计算当前残差，得到构造3d重建约束为：(4)、构建3d点的共面约束4.1)、第q幅标定板图像中理想标定板角点所在平面的法向量用相机外参的旋转矩阵r的第三列表示，即向量4.2)、第q幅标定板图像重建出的三维角点点集ω
q
＝{
q
p
i
′
|[
q
x
i
′
,q
y
i
′
,q
z
i
′
]
t
}中，记从角点
q
p
i
′
，i
′
＝1,2,
…
,m
×
n出发到角点
q
p
j
′
，j
′
＝1,2,
…
,m
×
n,j
′
≠i
′
的向量的单位方向向量为即：4.3)、将向量和单位方向向量内积的绝对值求和，作为衡量面上向量与面的法向量垂直程度的误差值，构成共面的约束项：(5)、构建3d点的共线约束5.1)、对于第q幅标定板图像，找到其重建的三维角点点集ω
q
＝{
q
p
i
′
|[
q
x
i
′
,q
y
i
′
,q
z
i
′
]
t
}中位于标定板在空间中第r行的首端点
q
p
rs
和尾端点
q
p
re
，将首尾端点连接得共线向量即5.2)、对于第r行的第k个中间点
q
p
rk
，同理，将该中间点
q
p
rk
与首端点
q
p
rs
连接得到第二
个向量5.3)、接着，计算向量在第r行的共线向量上的投影点，计算公式为：其中，是该行的共线向量的单位方向向量；5.4)、然后，计算中间点
q
p
rk
与其投影点间的距离，即中间点
q
p
rk
与第r行共线向量的距离，作为重建的共线误差中的一个残差项，对所有行与列，和所有点进行相同的处理，计算得到最终的共线误差项e
line
：其中，
q
p
k
′
c
、分别为第q幅标定板图像在空间中第c列的第k
′
个中间点及其在该列共线向量上的投影点；(6)、构建3d点的径向约束6.1)、第q幅标定板图像的三维角点点集ω
q
＝{
q
p
i
′
|[
q
x
i
′
,q
y
i
′
,q
z
i
′
]
t
}，其编号为从左到右、从上到下，选取四个角点
q
p1、
q
p
n
、
q
p
m
×
n
、
q
p
(m-1)
×
n+1
；6.2)、则径向约束的两项误差为：其中，radius为理想标定板径向长度；(7)、构建3d点的角度约束7.1)、将步骤6.1)选出的四个角点
q
p1、
q
p
n
、
q
p
m
×
n
、
q
p
(m-1)
×
n-1
分别记为点a、b、c、d，并构成一个四边形，将线段ab、ad、bc、dc的方向向量分别写为有：7.2)、对该四边形的四个角度进行计算，并将每一个角度值与直角进行残差计算，对于所有的标定图像，该项角度误差可以写为：(8)、考虑上述所有6个约束项，得到一个新约束即优化目标函数为：min e
total
＝e
pix
+e
pnt
+e
pl
+e
line
+e
diag
+e
ang
将待优化的相机标定参数写为如下形式：
将优化目标函数简写为：式中，θ表示上述待优化的相机标定参数的集合{a
l
,a
r
,d
l
,d
r
,r,t}，e是所有约束项函数平方根的集合(9)、初始化相机标定参数变量的初始值θ
←
θ0，采用列文伯格-马夸尔特方法，依据优化目标函数对待优化的相机标定参数进行最优化，得到相机标定参数的最优解θ
*
，其中，θ0为标定后相机标定参数。

技术总结
本发明公开了一种双目立体视觉相机标定参数的优化方法，对不同情况构造了相应的几何约束，使得标定板三维重建后的结果与实际三维场景保持一致，即重建点的坐标应与实际空间坐标一致且点与点之间的相对位置不变，同时标定板在空间中的几何特性应经过标定重建后保持不变，并将所有约束项进行叠加构成最终的优化目标函数，最后采用列文伯格-马夸尔特方法对目标函数进行最优化，得到优化后的相机标定参数。本发明提出的相机标定参数优化算法相较于张正友相机标定法有更高的准确度，使用优化后的参数进行标定后，重投影误差和三维重建测量结果误差大大降低，提高了相机标定参数精度。提高了相机标定参数精度。


技术研发人员：殷春 王泽琪 杨澳琳 陈凯 王文 朱丹丹 程玉华
受保护的技术使用者：电子科技大学
技术研发日：2023.05.09
技术公布日：2023/8/13