High-precision camera calibration method considering projected circular edge blur and eccentricity error

通过立体视觉测量系统可以获得物体在世界坐标系下的位姿信息，其中一个重要的前提是已知该系统中各个相机的内外参数，相机标定就是获取该参数的过程^[1]，相机标定得到的参数精度决定着三维重构后的效果^[2-3]。传统相机标定方法依赖于高精度的空间参照物，基于针孔相机模型以及世界坐标系和像素坐标系的对应关系建立相机参数约束并求解。传统标定方法中的张正友标定法因其有较好的鲁棒性且不需要昂贵的合作标靶进而备受关注^[4]，但是基于棋方格的标定中使用的角点特征易受到光照以及噪声的影响，不利于在实际工业检测中的应用。而圆心特征作为一种重要的特征点，比角点有更高的稳定性，对噪声以及图像模糊的情况有更好的抗干扰能力^[5-6]，因此圆形标靶得到了广泛的应用。

利用圆形标靶进行相机标定时，由于标靶平面与相机光轴不垂直以及成畸变等因素，空间圆形特征在成像平面的投影一般是椭圆，且圆心的投影点与投影的中心并不重合^[7]。传统的做法是将投影椭圆的拟合圆心直接作为中心的投影点^[8]，由于存在投影误差所以标定得到的参数并不精确。魏振忠^[7]等人对投影误差进行了建模仿真，分析了影响该偏差大小的因素，并得到了误差补偿的数学模型，该模型指出误差大小与物像平面距离以及圆形特征平面和成像平面的夹角有关，但是实际应用中无法获取这些精确数据，故而限制了此法在工程领域的应用。解则晓^[9]等人提出了一种基于对偶二次曲线提取平面标靶圆形投影的方法，在圆形特征平面与成像平面存在大倾角时也可得到高精度的投影中心坐标，但此法依赖于特定的圆形特征几何关系，每一步迭代均需要三个非共线圆形标记点，标定过程中增加了计算量。卢晓东^[10]等人则通过反投影虚拟标靶进行反复迭代以降低投影偏心误差对标定精度的影响，基于此法标定的相机参数精度较之传统方法有了提升但是标定的计算量却是传统方法的数倍。此外，基于同心圆，椭圆锥等方法^[11-12]利用几何关系约束条件，可以得到高精度的圆心投影坐标，但是对于不存在此类几何关系的圆形标靶并不适用。

此外，相机标定的过程中由于环境干扰以及图像采集硬件噪声的影响，得到的圆形特征图像存在边缘模糊现象^[13]。传统的边缘定位方法主要有插值法，拟合法，矩方法^[14]。黄霄霄^[15]采用Zernike矩对存在边缘模糊问题的红外图像进行亚像素定位，比传统算法的定位精度提升了0.19个像素等级。陈阔^[16]等人对Canny边缘检测算法进行改进，对边缘的水平与垂直方向梯度进行修正最后使用高斯拟合获取高精度的边缘坐标，但是抗噪能力较差。王宁^[17]等人通过构建Bertrand曲面来对边缘像素进行拟合定位，算法的定位精度高，抗干扰能力强，但是由于该模型的计算量较大，不适用于快速边缘检测。

基于上述分析，影响基于圆形特征的相机标定精度的主要问题为圆形特征边缘模糊和投影偏差。文中充分考虑了传统圆形特征点中心坐标提取中存在的问题，引入高斯误差函数，构建了非线性的边缘灰度分布函数来精确描述圆形特征的模糊边缘，并基于此边缘灰度模型计算图像的改进Zernike矩，由于Zernike矩有抗噪能力强、定位精度高的特点，可以对圆形特征进行高精度的定位。然后利用求解的边缘坐标进行迭代拟合获得圆形投影的中心坐标，与此同时，文中分析了拟合的投影圆心与真实中心投影偏差大小的影响因素，对获得的拟合中心坐标进行了偏差补偿以逼近真实的中心投影，最终使用高精度圆心坐标对双目视觉系统进行标定进而获取精确的相机参数。实验结果表明，文中方法提取的投影中心更接近真值，可以提高双目立体视觉系统的标定精度。

参考文献[18]建立了一种理想的边缘灰度模型，利用Zernike矩的正交性以及旋转不变性对理想阶跃边缘模型进行描述。图1为亚像素边缘检测理想模型，其中直线$l$被单位圆包围的区域代表理想边缘，其两侧的灰度值分别为$h$和$h + k$,$k$代表背景与前景图像之间的灰度差，$l$为图像坐标系原点到理想边缘的距离,$\varphi $为$l$和$x$轴的夹角，图1(b)为图1(a)绕原点顺时针旋转$\varphi $角度的结果。

Figure 1.  Ideal step edge model

针对一幅图像，其Zernike矩的形式为：

式中：${V_{nm}}\left( {\rho ,\theta } \right)$表示Zernike多项式；$n$与$m$表示Zernike矩的阶数；“*”表示共轭变换。

Zernike多项式在极坐标下可以表示为：

Zernike矩具有旋转不变性：

根据图1(b)可得，经旋转角度$\varphi $的模型关于$x$轴对称，即质心落在$x$轴上，满足公式(4)的关系：

根据上述等式关系以及Zernike矩的性质，可以列出如下公式：

理想边缘模型的各个参数可以通过公式(5)求得：

因此图像边缘的子像素边缘位置可以表示为：

式中：$\left( {{x_s},{y_s}} \right)$与$\left( {x,y} \right)$ 分别表示在图像像素坐标系下子像素边缘坐标与Sobel算子提取的边缘坐标。

对于一幅图像，Zernike矩是由模板与图像进行卷积计算得到，对于大小为 $N \times N$ 的模板，考虑其放大的作用，边缘像素点可以进一步表示为：

上述基于Zernike矩的亚像素边缘定位方法将实际的边缘简化为三个参数，分别为背景灰度值，边缘灰度阶跃以及边界位置，作为理想化的边缘模型，其并未考虑图像采样过程中的相机透射畸变，传感器噪声以及环境不均匀光照分布等因素，因此不能精确描述图像中边缘的灰度变化。实际图像中的边缘可用连续变化的灰度值表示，且沿着不同的方向灰度值遵循不同变化规律，如图2所示。

Figure 2.  Grayscale distribution of projected ellipse

考虑到边缘处的灰度连续变换，引入分段斜坡灰度模型以更精确地描述灰度值变化规律：

式中：$h$为背景灰度值；$k$为背景与目标灰度差值；${l_1}$为边缘的初始位置；${l_2}$为边缘的终止位置，图中的边缘位置$l$为${l_1}$，${l_2}$的均值。对于${l_1} < x < {l_2}$的边界区域，考虑到灰度变化应为非线性，故对斜坡模型的线性部分使用公式(10)所示的高斯误差函数(erf，error function or Gauss error function)进行优化，得到公式(11)所示的高斯边缘灰度模型：

最终可得到描述边缘灰度连续不均匀变换的高斯灰度模型：

公式(12)考虑了由于光照以及相机采样角度等因素引起的图像边缘灰度不均匀分布，$\sigma \left( {{x_i},{y_i}} \right)$为坐标$\left( {{x_i},{y_i}} \right)$为中心的邻域内像素灰度标准差，基于参数估计的方法构建了精确的边缘灰度模型，图3为边缘模型在二维坐标下的表示。

Figure 3.  Two-dimensional representation of edge gray distribution model

将上述边缘模型代入Zernike矩中即可得到考虑边缘模糊的改进Zernike矩，如下式：

则有$Z_{00}'$、$Z_{11}'$、$Z_{20}'$、$Z_{31}'$的值为：

从而进一步得到新的边缘像素定位表达式：

式中的$\Delta k$为：

文中提出的改进Zernike矩考虑了图像中边界像素灰度分布情况，将边界的区域表示为灰度值的连续变化而非灰度值阶跃的形式，更符合实际的情况，同时结合相机成像系统硬件噪声以及环境对相机成像的干扰等因素，利用高斯误差函数进行边界过渡区域的灰度值估计，精确地描述了边界处像素灰度分布情况，利用此边缘模型计算的Zernike矩可以求解得到准确的边缘像素坐标，为后续圆形特征圆心的高精度定位提供了良好的数据基础。

对得到的边缘像素坐标进行椭圆拟合，椭圆方程具有如下形式：

进一步表示为矩阵形式：

其中：

拟合椭圆的最优化问题可表示为：

其中：

对于上述的约束目标，引入拉格朗日乘子，利用拉格朗日乘子法进行求解即可求得各项参数得到拟合椭圆。但是由于图像采集时环境的干扰以及设备自身硬件的噪声，经边缘定位的边界点包含目标边缘以及部分噪声，为了剔除这部分数据，文中采取迭代拟合椭圆的方法，对每次迭代得到二次曲线进行残差分析，滤除残差较大的边界点进行二次拟合，重复这些步骤直到所有的保留边界点的残差值全部小于设定的阈值则认定拟合椭圆为投影椭圆。

根据针孔成像模型，当空间的圆形特征物理平面不严格垂直于相机的光轴时，空间圆的透视投影成像应为一个椭圆，此变换过程的模型如图4所示。

Figure 4.  Spatial circular projection model

图4中，p点为空间圆的中心，其在像平面的投影点为${p_t}$，而投影椭圆的几何中心为${p_d}$，在一般情况下${p_t}$与${p_d}$不重合，且两点的像素坐标分别为：

公式(19)中各符号的意义为：

式中：$({X_i},{Y_i})$为空间圆上各点的世界坐标；${\eta _i}$为半径；$t$为相机外参中的平移分量；${r_i}$为旋转矩阵中的元素。又已知这两点并不重合，以两点之间的偏移量作为圆形特征投影的畸变误差：

分析上式所述畸变误差，根据参考文献[7]魏振忠等人的仿真分析结果，此畸变项受空间圆心和像平面距离以及空间圆平面和像平面的夹角影响，误差会随着物象距离的增大而减小，随着物象夹角的减小而减小，误差最小可达$5.2$ μm，而在文中的相机标定过程中，物象距离较远且圆特征与相机成像平面的夹角较小，因此获得高精度逼近投影椭圆几何中心的拟合圆心坐标对于获得更好的标定效果很有必要。

对得到和拟合圆心进一步处理逼近真实的投影椭圆几何中心。处理边缘坐标，获取距离最远的两个点，其连线认定为椭圆长轴，中点为椭圆圆心，比较中点与拟合圆心的距离，此距离小于设定阈值将两点的中点作为投影圆心坐标，否则计算两点连线与投影椭圆的交点$({u_1},{v_1})$和$\left( {{u_2},{v_2}} \right)$，根据射影变换的简比不变与直线不变原则，得到高精度的圆心投影点坐标$\left( {{u_o},{v_o}} \right)$为：

式中：R的为空间圆形特征的半径。

通过以上方法可以对得到的拟合椭圆圆心进行补偿误差，得到高精度逼近真实投影点的坐标值，算法流程图如图5所示。

Figure 5.  Projection ellipse fitting algorithm

基于亚像素的灰度图像边缘定位方法有：灰度矩法，拟合法，Zernike矩法等，为了对文中所提出的改进Zernike矩算法进行定位精度的评估，将其分别与上述三种方法在仿真实验中进行比较，并对改进算法的鲁棒性和定位精度进行分析。

为了对上述几种亚像素边缘定位算法的定位精度和鲁棒性进行比较，文中基于绘图工具绘制了大小为1 680$ \times $1 680 pixel的黑底白圆模拟图像，如图6所示，圆心坐标为(256,256)，半径$r = {\rm{5}}0$ pixel，为了模拟实际情况，将方差为0，0.002，0.004的高斯白噪声加入到模拟图像中，通过分析上述三种传统方法以及文中算法对模拟图像边缘提取的误差分布和圆心定位偏差大小来验证所提方法的有效性和鲁棒性。

Figure 6.  Simulative image

图7为模拟图像加入了方差为0.002的高斯白噪声时，用四种方法对圆形特征边缘定位与边缘真值的误差曲线对比图。从图中可以看出，对于相同噪声水平的模拟图像，文中所提算法的定位误差明显小于其他三种传统算法，误差分布更接近0，由此可见文中算法对图像噪声的敏感度较低，有较好的抗干扰能力。

Figure 7.  Error distribution of ${\sigma ^{\rm{2}}}{\rm{ = 0}}{\rm{.002}}$

对每次实验得到的边缘点坐标进行随机抽样，获取50个样本数据并计算边缘点与圆心的距离，与半径真值作差并计算误差绝对值，然后计算误差绝对值的均值与方差，得到表1的数据。可以看出，在没有噪声干扰的情况下，四种拟合方法都有较高的定位精度，但是随着加入的噪声增强，四种算法的定位误差均逐渐增大，但是文中算法计算的边缘坐标误差均值最小，且相对于其他算法的抗干扰性更强，在不同的噪声等级下误差均值与方差都是最小的，因此可以看出文中算法的抗噪性能最好，精度最高。

得到了模拟图像的边缘像素坐标之后，利用这些边界点进行迭代拟合获取其圆心坐标。基于文中算法对模拟图像处理，得到的边缘点以及拟合椭圆如图8所示。

Figure 8.  Edge and center localization based on improved Zernike moments

计算每次提取得到的拟合圆心与模拟图像圆心真值的距离偏差值，得到表2中的数据。从表2中数据可以看出，在相同噪声水平下，文中提出的算法对于圆形特征的圆心定位精度是最高的，相较Zernike矩，文中算法由于使用了更精确的高斯边缘模型计算得到的改进Zernike矩，精度较之提升了约25%。与此同时，由于文中在拟合计算投影椭圆曲线时采用了迭代拟合，滤除野点的方法，可以降低噪声对定位的干扰，进一步增强算法的抗干扰能力。

相机标定的精度会受到圆形标志圆心提取误差，标定板精度，镜头畸变等因素的影响，通过仿真实验可以看出文中提出的算法对于圆形特征中心有较高的定位精度，且具有一定的噪声抑制能力，可以在标定过程中得到精确的圆形标志点图像坐标为双目视觉系统的标定提供良好的数据基础。

通过仿真实验验证了文中算法能够对圆心特征投影中心进行高精度的提取，为了进一步验证算法对相机标定效果的提升，设计实际标定实验计算相机参数并通过标准杆来对相机标定精度进行检验。构建双目立体视觉系统架构如图9所示，采用两台由Mikrotron生产的MC3010相机，相机图像分辨率为1680×1720 pixel，像元尺寸为${\rm d}u \times {\rm d}v =$7.4 μm×7.4 μm，相机镜头型号均为AF Zoom-Nikkor 24-85 mm/1:2.8-4D，用其拍摄99圆形标定板进行标定完成文中标定方法的验证。

Figure 9.  Acquisition equipment of stereo camera system image

设置圆形标定板位于双目中线方向约1200 mm的位置，通过图9中所示的双目立体成像系统对其进行图像采集，并对标定板上的99个圆形标志点进行圆心提取与编号，如图10所示。

Figure 10.  Extraction and coding of circular mark points

调整平面标定板的姿态，空间位置，使其占据视觉系统的视场大部分区域，通过文中算法对标定板的圆形特征圆心进行高精度提取，并将得到的圆心依次进行编码，图11为部分处于不同位姿状态的标定板图像及编码情况。

Figure 11.  Target image acquisition and center extraction

此次标定实验中双目立体视觉系统左右相机分别采集了25张标定板图像，对所有图像中的圆形标志点利用文中方法进行圆心的高精度定位求解图像坐标，并依据图像和标定板的对应关系以及标定板各个圆形特征的几何关系计算双目相机系统的相机参数。表3为通过文中算法解算出的相机参数以及双目之间的位姿关系。图12为基于标定结果重构的标定板特征点空间三维分布效果图以及双目相机和标定板的空间位姿关系。

Figure 12.  Stereo vision system calibration results

然后利用如图13所示的标准杆对双目视觉系统的标定结果进行精度验证，将标准杆放置于立体视觉系统的视场范围内并采集图像，分别利用灰度矩，拟合法，Zernike矩以及Hough变换等传统做法与文中算法得到的相机参数对标准杆采样图像进行重构，获得标准杆的两端标志点的坐标，考虑到物体在世界坐标系下的实际坐标不易获得，因此设计实验使用文中算法标定获取的相机参数对标准杆件进行长度测量，然后与标准杆长度真值比较并计算测量误差，各个方案多次测量得到的平均结果如表4所示。

Figure 13.  Calibration rod

表4中的数据显示，基于文中算法优化的标定方法得到的双目视觉系统参数将测量误差值控制在0.2 mm以内，明显小于其他四种传统做法的测量结果，测量精度提升约30%，可以看出文中算法在提升标定精度上具有一定的优势。这是因为传统的标定做法中，忽略了圆形特征投影时存在的偏心误差以及圆形特征投影的边缘模糊问题，在计算标定板圆形特征的圆心时引入了较大的误差；而文中算法通过建立了高斯边缘模型精确描述了圆形特征的边缘并且补偿了偏心误差，因此能够实现更高精度的相机标定。

文中提出的算法考虑了由于干扰引起的圆形特征投影边缘模糊以及投影偏心问题，提出了更适合描述此类边缘灰度分布的边缘模型并基于此模型计算图像的Zernike矩实现亚像素定位，最终利用高精度的定位坐标数据进行拟合并补偿圆心误差，以此圆心坐标进行相机标定得到的相机参数相对传统方法有更高精度，但是此算法也因为加入了更多的计算，在标定过程中相比传统做法存在运行时间长的问题，这也将是未来工作需要着力去改进的地方。

文中考虑了在相机标定过程中由于相机本身的传感器噪声以及环境光照分布不均等因素引起的圆形特征投影边缘模糊现象，引入高斯误差函数建立了高斯边缘灰度分布模型，并基于此模型计算图像的Zernike矩来实现边缘的亚像素定位，然后利用高精度的亚像素坐标迭代拟合得到投影椭圆的拟合圆心坐标，并通过分析投影几何圆心和真实投影中心的偏离关系，对拟合圆心进行补偿误差，最终基于此圆心坐标对立体视觉系统进行标定得到相机参数。通过仿真实验和实测实验可以验证该算法的有效性和准确性，文中提出的算法在边缘像素坐标提取和圆心拟合上均优于四种传统算法，对于标准杆的三维重建测量精度也比传统算法提升了30%，能够满足高精度的三维测量任务，具有实际的工程应用价值。

下一步工作将针对文中算法计算量大、标定时间久的问题进行改进，可以通过GPU进行并行运算，基于异构加速提高算法的运行效率。