CN103247053B - 基于双目显微立体视觉的零件精确定位方法 - Google Patents

Abstract

本发明基于双目显微立体视觉的零件精确定位方法属于计算机视觉测量技术领域，涉及一种基于双目显微立体视觉的精密零件的精确定位方法。采用双目显微立体视觉系统，利用两个CCD摄像机采集被测零件的图像，立体显微镜将被测零件上的待测区域的图像信息放大；采用棋盘格标定板对两个CCD摄像机进行标定；采用Harris角点检测算法及亚像素提取算法进行特征点的提取。将提取后的特征点进行初匹配和匹配点对的校正，将特征点图像坐标输入到标定好的系统中得到特征点的空间实际坐标。本发明解决了由于待测目标区域小、定位精度要求高、非接触等问题所产生的测量难题。采用基于双目显微立体视觉的非接触式测量方法，很好的完成了精密零件的精确定位。

Description

基于双目显微立体视觉的零件精确定位方法

技术领域

本发明属于计算机视觉测量技术领域，特别涉及一种基于双目显微立体视觉的精密零件的精确定位方法。

背景技术

双目立体视觉技术通过模拟人类双眼的处理方式，具有获取待测物体深度信息的能力，进而能够获得被测物体的空间位置信息，同时又具有无损测量及实时测量的优点，在各行各业得到了广泛的应用。双目显微立体视觉技术是建立在双目立体视觉技术和体式显微镜技术之上的，通过左右两个CCD摄像机分别捕捉立体显微镜两条光路的图像，通过捕捉的图像存在同双目视觉测量时类似的角度差，能够实现对目标的三维测量。双目显微立体视觉技术在保证非接触性、实时性等优点的同时能够精确地获得目标点的空间位置信息。其可应用于三维高精度伺服控制、高精度空间定位以及高精度三维测量等方面。零件的加工定位精度直接影响着零件的加工精度，而对于精密零件定位精度有可能直接影响加工零件是否可用，同时针对精密超精密零件的加工定位往往要求非接触的特点，基于双目显微立体视觉的精密零件的精确定位方法应运而生。

重庆大学的汤宝平等人申请的发明专利CN102567989A“基于双目立体视觉的空间定位方法”针对空间定位存在的摄像机标定过程繁琐、定位精度不高问题提出了一种新的基于双目立体视觉的空间定位方法，该方法通过对摄像机进行标定求得摄像机的单应性矩阵和畸变系数，再通过标定结果获得目标点的理论坐标，最后结合目标点的实际坐标计算出目标点的空间坐标。然而，由于双目显微立体视觉系统具有更加复杂的光路，焦深小，畸变因素较多，视场窄，而目前广泛采用的摄像机标定方法都是建立在大视场下实现的，因此上述方法并不能解决基于双目显微立体视觉的精密零件的精确定位。

发明内容

本发明要解决的技术难题是克服现有技术的缺陷，发明一种基于双目显微立体视觉的精密零件的精确定位方法，解决由于待测目标区域小、定位精度要求高、非接触等问题所产生的测量难题。采用基于双目显微立体视觉来实现精密零件的精确定位，拓宽了传统的基于双目立体视觉的空间定位方法的应用范围，并提高了测量的精度，解决了待测零件区域小，定位精度要求高，难以测量的问题。

本发明所采用的技术方案是基于双目显微立体视觉的精密零件的精确定位方法，其特征是，采用双目显微立体视觉系统，利用左右两个CCD摄像机2、2’采集被测零件4的图像位置信息，立体显微镜3将被测零件4上的待测区域6的图像位置信息进行放大；放大后的图像经过图像采集卡传到工业控制计算机9内，采用高精度较的棋盘格标定板对左右两个CCD摄像机进行标定；采用Harris角点检测算法及亚像素提取算法进行待测特征点的提取；将提取后的待测特征点进行初匹配和匹配点对的校正；将已经匹配好的特征点图像坐标输入到标定好的系统中得到特征点的空间实际坐标；测量方法的具体步骤如下：

(1)左右两个CCD摄像机的标定

左右两个CCD摄像机的标定包括摄像机内参数和摄像机外参数；摄像机内参数包括尺度因子α、β，主点坐标u₀，v₀，以及垂直因子γ；在标定过程中，需要先得到摄像机的五个内参数，在求得内参数矩阵的基础上，求解摄像机外参数矩阵；尺度因子是空间点在经过平移旋转变换之后，它在摄像机坐标系中的坐标与它在图像坐标系中的坐标的尺度关系，本发明对尺度因子的标定采用张氏摄像机标定法；根据摄像机的针孔模型可建立图像坐标系与世界坐标系的关系公式为：

$Z_c\left[\begin{array}{c}u\\ v\\ 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}{cc}R& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]---\left(1\right)$

其中：α和β就是需要标定的尺度因子，X_w,Y_w,Z_w是空间中一点P的三维空间坐标，u,v是P点在图像上的图像坐标，R,t代表了摄像机坐标系相对于世界坐标系的旋转和平移矩阵；若令X′_w＝(X_w,Y_w,Z_w)^T，Z_c＝s,则

$s{\stackrel{\‾}{x}}_w={\mathrm{HX}}_w^{\′}---\left(2\right)$

其中，H也称为单应性矩阵，且有

H＝K[r₁,r₂,t]    (3)

每一个位置的标定板图像都独立对应一个单应性矩阵，令H＝[h₁,h₂,h₃]则可由上式推得：

[h₁ h₂ h₃]＝λK[r₁ r₂ t]    (4)

其中，λ为任意的比例因子，K摄像机内参数矩阵，因r₁和r₂是单位正交向量，即r₁ ^Tr₁＝r₂ ^Tr₂＝1且r₁ ^Tr₂＝0则：

$h_1^T{K}^{-T}{K}^{-1}{h}_1=h_2^T{K}^{-T}{K}^{-1}{h}_2=1---\left(5\right)$

h₁K^-TK^-1h₂＝0    (6)

根据公式5公式6，计算出立体显微镜左右光路的摄像机尺度因子；针对主点坐标，采用变倍率法对摄像机的主点进行标定；假设空间中任意一点P在摄像机坐标系中的坐标为x₁,y₁,z₁，在不考虑非线性畸变以及图像坐标垂直度的情况下，该点的图像平面投影的坐标为：

u₁＝rx₁+u₀

                                                              (7)

v₁＝ry₁+v₀

其中u₁,v₁是该点的图像坐标，r是任一放大倍率，u₀,v₀是主点坐标；将式中的r消去，可得直线方程：

$\frac{u_1-u_0}{x_1}=\frac{v_1-v_0}{y_1}---\left(8\right)$

即在任何放大倍率下，点x₁,y₁,z₁的图像坐标都在同一直线上，且该直线一定经过摄像机主点u₀,v₀；为了求取摄像机主点u₀,v₀，选定12个在图像投影平面上不重合的点，令这12个点分别在不同放大倍率下0.7、1、2、2.4、3、5、5.8、8进行投影，利用最小二乘法对获得的12条直线进行拟合，并用最小二乘法求取12条直线的交点坐标即摄像机主点u₀,v₀；由于垂直因子γ对成像精度影响不大，因此在初始参数标定时，将γ设为0，通过优化手段得到γ具体值；针对摄像机外参数旋转矩阵R和平移矩阵t；由公式2知，在标定完内参数矩阵K和单应性矩阵H后可以求得摄像机旋转矩阵R＝[r₁,r₂,r₃]和平移向量t如下：

r₁＝λK^-1h₁

r₂＝λK^-1h₂    (9)

r₃＝r₁×r₂

t＝λK^-1h₃

在初步标定完左右两个CCD摄像机的主参数后，下一步是对主参数和畸变参数的优化；为了简化优化参数，本发明采用四元数法将旋转矩阵中的九个未知数缩减为四个，公式如下：

$\begin{array}{c}{q}_1=\frac{\sqrt{1+r_{11}+r_{22}+r_{33}}}{2},q_2=\frac{r_{23}-r_{32}}{{4q}_1},q_3=\frac{r_{31}-r_{13}}{{4q}_1},q_4=\frac{r_{12}-r_{21}}{{4q}_1}\\ q_2=\frac{\sqrt{1+r_{11}-r_{22}-r_{33}}}{2},q_1=\frac{r_{23}-r_{32}}{{4q}_2},q_3=\frac{r_{12}+r_{21}}{{4q}_2},q_4=\frac{r_{31}-r_{13}}{{4q}_2}\\ q_3=\frac{\sqrt{1-r_{11}+r_{22}-r_{33}}}{2},q_1=\frac{r_{31}-r_{13}}{{4q}_3},q_2=\frac{r_{12}-r_{21}}{{4q}_3},q_4=\frac{r_{23}-r_{32}}{{4q}_3}\\ q_4=\frac{\sqrt{1-r_{11}-r_{22}+r_{33}}}{2},q_1=\frac{r_{12}-r_{21}}{{4q}_4},q_2=\frac{r_{31}+r_{13}}{{4q}_4},q_3=\frac{r_{23}-r_{32}}{{4q}_4}\end{array}---\left(10\right)$

从公式(10)计算得到新的旋转矩阵：

$R=\left[\begin{array}{ccc}-1+{2q}_2^2+{2q}_1^2& {2q}_2{q}_3-{2q}_1{q}_4& {2q}_1{q}_3+{2q}_2{q}_4\\ {2q}_2{q}_3+{2q}_4{q}_1& -1+{2q}_1^2+{2q}_3^2& -{2q}_1{q}_2+{2q}_3{q}_4\\ {-2q}_1{q}_3+{2q}_4{q}_2& {2q}_1{q}_2+{2q}_3{q}_4& -1+{2q}_1^2+{2q}_4^2\end{array}\right]---\left(11\right)$

由于四元数q₁,q₂,q₃,q₄的平方和等于一，因此旋转矩阵的待优化参数简化为三个；将已求得的左右两个CCD摄像机标定参数值作为优化初始值，利用优化算法对摄像机主参数和畸变参数进行整体优化；对左右两个CCD摄像机主参数和畸变参数的优化是基于最大似然估计准则，对于给定的n个标定板图片提供的m×n个标定点坐标，左右两个CCD摄像机参数的最优化问题通过下式的最小化问题来表达:

$S\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_i^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{[y_i-m_i(C_j,p_j,X_j)]}^2---\left(12\right)$

其中j是指参与计算的第j个摄像机，i是指第j个摄像机获取的第i个点，X_i是输入的空间点坐标，y_i是空间第i点的图像坐标，C_j是固定不变的摄像机参数向量，其长度为n₀，p_j是需要调整的摄像机参数向量，其长度为n₁，n₀+n₁即为摄像机所有参数向量长度，m_i(C_j,p_j,X_i)是摄像机的成像方程；本发明采用基于LM算法的光束平差法来解决上式的最小化优化问题；利用光束平差法的编写优化程序时除了要提供左右两个CCD摄像机的标定参数初值外还要提供一定数量的和形式的三维空间坐标点以及这些三维坐标点的相应图像坐标，因此为了精确得到整个视场空间内高精密标定板的角点坐标，采用数控机床z轴带动显微镜纵向移动的方法；并采用激光干涉仪对z轴进行测量得到z轴移动的实际距离，保证所采集角点坐标点的z向实际精度；在光束平差法中，优化算法的迭代因子是由迭代参数的一阶导数组成的雅克比矩阵构成的，因此需要在计算过程中引入待成像方程中调整参数的雅克比矩阵；

$J=[\frac{\mathrm{df}}{\mathrm{dp}\left(1\right)},\frac{\mathrm{df}}{\mathrm{dp}\left(2\right)},....,\frac{\mathrm{df}}{\mathrm{dp}\left(n_1\right)}]---\left(13\right)$

光束平差法以其待调整参数的雅克比矩阵组成迭代因子，对模型参数进行调整，得到使成像误差S(θ)达到最小的模型参数最优解。

(2)待测区域特征点的提取

本发明针对待测区域特征点的提取采用的是Harris角点检测算法；Harris提出了用高斯函数代替方形区域计算梯度的算法，若目标点的像素坐标为x,y，其在x方向和y方向的位移分别为u和v,则该目标点的灰度变化由如下公式表示：

$\begin{array}{c}E(x,y)=\underset{u,v}{\mathrm{\Σ}}\mathrm{\ω}(u,v){[I(x+u,y+v)-I(x,y)]}^2\\ \≈\underset{u,v}{\mathrm{\Σ}}\mathrm{\ω}(u,v)(u,v)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]\left(\begin{array}{c}u\\ v\end{array}\right)\end{array}---\left(14\right)$

令 $M(x,y)=\underset{u,v}{\mathrm{\Σ}}\mathrm{\ω}(u,v)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$ 则：

$E(x,y)\≈(u,v)M(x,y)\left(\begin{array}{c}u\\ v\end{array}\right)---\left(15\right)$

计算矩阵M的特征值λ₁,λ₂，如果两个特征值都比较大，说明目标点的图像灰度自相关函数的两个正交方向上的值均较大，则该点即为特征点；通过Harris算法对特征点进行初提取后；采用亚像素角点提取算法提取出高精度的角点坐标；对于理想角点，它附近的像素点的灰度梯度方向均垂直于该点与理想角点的连线；这一特点可以用公式表达为：

$\▿\stackrel{\→}{H}(\stackrel{\→}{\mathrm{\α}}-\stackrel{\→}{\mathrm{\β}})=0---\left(16\right)$

其中是理想角点的灰度梯度方向，向量为图像原点指向理想角点的坐标，向量为图像原点指向理想角点附近任一边缘点的坐标；实际上，由于受到图像噪声的影响，公式(16)通常都不为零，可以将其视为误差e，即

$e=\▿\stackrel{\→}{H}(\stackrel{\→}{\mathrm{\α}}-\stackrel{\→}{\mathrm{\β}})---\left(17\right)$

在以角点为中心的邻域内，对所有点按上式计算，误差和为E,则

$E=\underset{i}{\mathrm{\Σ}}e=\underset{i}{\mathrm{\Σ}}\▿\stackrel{\→}{H}(\stackrel{\→}{\mathrm{\α}}-\stackrel{\→}{\mathrm{\β}})---\left(18\right)$

其中，i是邻域内第i个点；这样求使误差和E最小的点即是亚像素级别的角点坐标。

(3)特征点的立体匹配

采用归一化交叉算法(NCC)实现特征点的初匹配，以CCD摄像机2拍得的图片作为基准图片，并以特征点p_l为中心,构造一个N*N的局部图像块作为模版图T，并令模版图T在待匹配图像S的包含有特征点的一定范围内进行遍历，搜索模版所覆盖的子图记做S^i,j，i,j为子图中心点在匹配图像S中的像素坐标，通过使用归一化交叉算法(NCC)来计算模版图T和搜索子图S^i,j之间的归一化互相关系数，算法如下:

$R(i,j)=\frac{\underset{m,n}{\mathrm{\Σ}}|S^{i,j}(m,n)-E\left(S^{i,j}\right)||T(m,n)-E\left(T\right)|}{\sqrt{\underset{m,n}{\mathrm{\Σ}}{(S^{i,j}(m,n)-E\left(S^{i,j}\right))}^2\underset{m,n}{\mathrm{\Σ}}{(T(m,n)-E\left(T\right))}^2}}---\left(19\right)$

其中m,n代表各个像素的坐标，E(S^i,j)和E(T)分别是搜索子图S^i,j与模版图T的灰度平均值；互相关系数R(i,j)的值越大，则搜索子图与模板图的匹配程度越高；设定互相关系数R(i,j)的阈值后，选出特征点p_l的候选匹配点；使用NCC算法进行初匹配后还需要利用外极限约束和距离约束来对匹配的特征点对进行校正；通过使用Longguet-Higgins提出的归一化8点算法在摄像机标定时计算左右两个CCD摄像机的基本矩阵F，根据外极线理论，对基准图像特征点p_l，在待匹配图像S中对应的外极线l_r可以表示为：

l_r＝F·p_l    (20)

初匹配的候选特征点若在外极线l_r附近则可进一步确定该候选特征点为匹配点；最后统计每一候选特征点到其他特征点的距离和基准图像中特征点到其他特征点点距离相等的个数num；num最大值时对应的候选特征点则为基准图像中特征点的正确匹配点。

(4)三维坐标的求取

首先计算出双目图像中特征点经亚像素算法提取出的图像坐标值，并大致估算出特征点在世界坐标系中的估计值；在优化程序中将左右两个CCD摄像机参数对应雅克比矩阵设为零，同时在优化程序计算迭代因子部分中嵌入关于空间坐标x、y、z的雅克比矩阵方程；将已经优化好的左右两个CCD摄像机参数以及特征点图像坐标值和世界坐标估计值输入光束平差法优化程序，就可以对空间点的三维坐标进行优化，得到特征点的三维空间测量值，从而完成三维坐标点的求取并最终实现特征点的定位。

本发明的有益效果是解决了双目显微视觉标定时透镜多、视场狭窄、景深小、具有一定放大倍率并包含有较多的畸变因素的技术困难，实现了双目显微立体视觉的高精度标定，从而完成了精密零件的精确定位。

附图说明

图1所示为基于双目显微立体视觉的精密零件定位的装置模型图。其中，1—为双目显微立体视觉系统提供照明的LED光源，2—左CCD摄像机，2’—右CCD摄像机，3—立体显微镜，4—待测量的工件；5—可倾数控转台；6—待测量区域；7—精密数控位移平台；8—床身；9—工业控制计算机。

图2所示为基于双目显微立体视觉的精密零件定位测量方法的流程图。

图3所示为基于光束平差算法的双目显微立体视觉的标定和三维坐标点的求取的流程图。

具体实施方式

以下结合技术方案和附图详细叙述本发明的具体实施方式。附图1为基于双目显微立体视觉的精密零件定位的装置模型图。本装置通过左右两个CCD摄像机2、2’采集待测工件中待测量区域的位置信息，通过已建立好的图像坐标系和世界坐标系之间的关系，找到待测量点的三维空间坐标即实现精密定位。其装置的安装方式如下：机床床身8置于地面；精密数控位移平台7通过导轨与机床床身安装在一起，其可以满足标定实验中对棋盘格标定板的高精度平移和角度旋转需求；可倾数控转台5通过螺栓固定在铸铁平台的T形槽上；待测工件4通过专用夹具安装在可倾数控转台上，其待测区域6位于双目显微立体视觉系统的视场内；立体显微镜3通过专用的显微镜夹具安装在机床床身上；LED光源1固定在夹具上，夹具通过螺栓固连在床身上；左右两个CCD摄像机2和2’通过螺纹连接安装在立体显微镜上；左右两个CCD摄像机采集到的图像分别通过1394图像采集卡传递到工业控制计算机9内进行图像数据处理。

附图1为本发明的一个实施例，采用左右两个CCD摄像机2、2’拍摄待测物体的图像位置信息，左右两个CCD摄像机采用的是奥林巴斯DP26摄像机，图像分辨率：2448*1920，采集速度：7fps，芯片尺寸：2/3英寸。立体显微镜3采用的是奥林巴斯生产的SZX-16研究级立体显微镜。变焦倍率：0.7-11.5，工作距离：60mm，采集图像最大视场：12.5mm×16.6mm，采集图像最小视场：0.76mm×1.02mm。为实现双目显微立体视觉系统的标定，标定所需的标定板采用的是深圳科创时代生产的CG-050-T-0.5型玻璃基底棋盘格标定板。该型标定板方格大小为0.5mm×0.5mm，棋盘格图案整体宽度为51mm×51mm，制造精度为1μm，能够满足立体显微镜标定的需要。

附图2表示了基于双目显微立体视觉的精密零件定位测量方法的流程图，精密定位的主要流程包括左右两个CCD摄像机2、2’的标定，双目图像特征点的提取，双目图像特征点的匹配与校正，双目图像特征点的三维坐标求取。其中，左右两个CCD摄像机的标定是利用高精密棋盘格标定板来实现对左右两CCD摄像机2、2’的主参数和畸变参数求解的，主要包括采用张氏标定法对尺度因子的标定、采用变倍率法实现主点坐标的标定和采用光束平差法对主参数、畸变参数的优化；双目图像的特征点利用Harris角点检测算法和亚像素提取算法进行提取；双目图像的标记点匹配是通过对双目图像的标记点进行初匹配及特征点对的校正来实现的，最终通过光束平差法实现三维坐标点的求取。

附图3所示为基于光束平差算法的双目显微立体视觉的标定和三维坐标点求取的流程图，通过拍摄得到的不同倍率下的图像使用变倍率法获得主点坐标，通过30副各种角度的图像使用张氏标定法获得尺度因子，并通过已获得的内参数求取左右两个CCD摄像机的外参数，将内参数和外参数作为光束平差法优化的初始值，并使用激光干涉仪获取高精度的立体点阵作为优化的约束条件，以待调整参数的雅克比矩阵组成迭代因子，通过光束平差法获得优化后的参数，将优化后的参数带入到成像模型中得到最终成像模型，最后把成像模型、待测点的图像坐标及三维坐标估计值输入到光束平差法的程序中求取待测点的空间三维坐标。

测量方法的具体步骤如下：

(1)左右两个CCD摄像机的标定

本发明采用摄像机相对固定的方式，以精密加工的棋盘格标定板对摄像机进行标定的方法。张氏摄像机标定法要求不同标定板图片有较大的角度关系，通常最佳角度为45°；而我们限于立体显微镜的景深因素，当标定板与水平面角度超过20°，视场内会出现许多模糊的角点，为了能够稳定的得到摄像机的尺度因子，并对张氏标定的实验过程进行总结，发现标定板图片数量越多，相互平行的标定板图片数量越少，就越能稳定的求得摄像机的尺度因子，因此对左右两CCD摄像机2、2’的尺度因子的标定的具体步骤如下：将标定板旋转分为沿水平轴旋转和竖直轴旋转两种方式，竖直轴旋转每次旋转40°，这样标定板在就存在九个水平角位置。在每个水平角位置进行水平轴旋转，将标定板水平轴旋转位置分为+10°，0°，10°三个位置，在每个位置上拍摄一张标定板图像。整个对尺度因子的标定过程共拍摄54张标定板图片。针对主点坐标，考虑到主点实际上是光轴与CCD成像面的交点，而在显微镜放大倍率发生改变时，光轴的位置实际上是不变的也就是摄像机的主点坐标始终保持不变，在任何放大倍率下，点(x₁,y₁,z₁)的图像坐标都在同一直线上，且该直线一定经过摄像机主点(u₀,v₀)；这样，n个在图像平面上的投影不重合的空间点，通过在不同放大倍率下进行投影，就能够形成n条相交于同一点的直线；这个交点的坐标就是摄像机主点坐标。采用变倍率法实现主点坐标的标定具体操作步骤：确定12个在图像投影平面上不重合的点，令这12个点分别在不同放大倍率下0.7、1、2、2.4、3、5、5.8、8进行投影，利用最小二乘法对获得的12条相交于一点的直线，求取交点坐标值。外参数矩阵在标定完内参数及单应性矩阵利用公式7进行求解。采用光束平差法对内外标定参数的优化的过程：在优化程序中输入成像方程、空间点图像坐标观察值、待优化的参数初值，程序的输出即为能使S(θ)输出最小的优化后的标定参数值。

为了验证标定结果的精度，利用标定后的模型参数，对标定板提供的三维空间点进行三维重建，得到三维空间点的测量值。通过将测量值与三维空间点坐标的实际值进行比较分析，来评价参数的标定精度。标定结果显示，在1倍率下，x方向、y方向、z方向的平均标定精度分别为：2.3um、1.4um、7.9um；在2倍率下，x方向、y方向、z方向的平均标定精度分别为：1.5um、1.2um、5.4um；在5倍率下，x方向、y方向、z方向的平均标定精度分别为：0.7um、0.5um、3.2um。标定结果表明标定后的模型参数具有较高的横向重建精度和较好的纵向重建精度。

(2)特征点的提取

对左右两个CCD摄像机2和2’采集到的图片通过使用采用Harris角点检测算法进行特征点X₁,X₂……X_n初提取，Harris角点检测算法是通过计算特征值λ₁,λ₂，如果两个特征值都比较大，说明目标点的图像灰度自相关函数的两个正交方向上的值均较大，则该点即为特征点。以提取后的特征点X₁,X₂……X_n为中心，分别计算每一个特征点的误差和E，求使误差和E最小的点即是亚像素级别的角点位置。通过这种方式，就能够在Harris角点提取的基础上得到亚像素角点坐标，提高特征点的提取精度。

(3)特征点的匹配

采用归一化交叉算法(NCC)实现特征点的初匹配，根据算法公式(19),以CCD摄像机2拍得的图片作为基准图片，并以特征点p_l为中心,构造一个N*N的局部图像块作为模版图T，并令模版图T在待匹配图像S的包含有特征点的一定范围内进行遍历，搜索模版所覆盖的子图记做S^i,j，i,j为子图中心点在匹配图像S中的像素坐标，通过使用归一化交叉算法(NCC)来计算模版图T和搜索子图S^i,j之间的归一化互相关系数。互相关系数R(i,j)的值越大，则搜索子图与模板图的匹配程度越高。设定互相关系数R(i,j)的阈值后，选出特征点p_l的候选匹配点。使用NCC算法进行初匹配后还需要利用外极线约束和距离约束来对匹配的特征点对进行校正；通过使用Longguet-Higgins提出的归一化8点算法在摄像机标定时计算左右两个CCD摄像机的基本矩阵F，根据外极线理论，对基准图像特征点p_l，使用公式(20)计算出基准图像中的特征点p_l在待匹配图像中的外极线位置，初匹配的候选特征点若在外极线附近则可进一步确定该候选特征点为匹配点。最后是通过距离约束进一步校正匹配点对。采用归一化交叉算法(NCC)进行特征点的初匹配，外极线约束和距离约束来校正匹配的特征点对，左右两个图像匹配的特征点对分别为：X_l,X_r'，X_2l,X_2r'……X_3l,X_3r'。

(4)三维坐标的求取

光束平差法除了能够对摄像机参数进行优化外，还能同时进行三维坐标点的优化；此处三维坐标点的优化即是三维坐标的求取；求取过程如下：将优化程序中摄像机参数对应雅克比矩阵设为零，同时在优化程序计算迭代因子部分中嵌入关于空间坐标x、y、z的雅克比矩阵方程；将已经优化好的左右两个CCD摄像机参数以及空间点的图像坐标和世界坐标估计值输入光束平差法优化程序，就可以对空间点的三维坐标进行优化，得到三维空间点的测量值，从而完成三维坐标点的求取。

本发明较好的解决了由于待测目标区域小、定位精度要求高、非接触等问题所产生的测量难题。采用基于双目显微立体视觉的非接触式测量方法，很好的完成了精密零件的精确定位。

Claims (1)

1.一种基于双目显微立体视觉的零件精确定位方法，其特征是，采用双目显微立体视觉系统，利用左右两个CCD摄像机(2、2’)采集被测零件(4)的图像位置信息，立体显微镜(3)将被测零件(4)上的待测区域(6)的图像位置信息进行放大；放大后的图像经过图像采集卡传到工业控制计算机(9)内，采用精度较高的棋盘格标定板对左右两个CCD摄像机进行标定；采用Harris角点检测算法及亚像素提取算法进行待测特征点的提取；将提取后的待测特征点进行初匹配和匹配点对的校正；将已经匹配好的特征点图像坐标输入到标定好的系统中得到特征点的空间实际坐标；测量方法的具体步骤如下：

(1)左右两个CCD摄像机的标定

左右两个CCD摄像机的标定包括摄像机内参数和摄像机外参数；摄像机内参数包括尺度因子α、β，主点坐标u₀，v₀，以及垂直因子γ；在标定过程中，需要先得到摄像机的五个内参数，在求得内参数矩阵的基础上，求解摄像机外参数矩阵；尺度因子是空间点在经过平移旋转变换之后，它在摄像机坐标系中的坐标与它在图像坐标系中的坐标的尺度关系，本发明对尺度因子的标定采用张氏摄像机标定法；根据摄像机的针孔模型可建立图像坐标系与世界坐标系的关系公式为：

$Z_c\left[\begin{array}{c}u\\ v\\ 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}{cc}R& t\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]---\left(1\right)$

其中：α和β就是需要标定的尺度因子，X_w,Y_w,Z_w是空间中一点P的三维空间坐标，u,v是P点在图像上的图像坐标，R,t代表了摄像机坐标系相对于世界坐标系的旋转和平移矩阵；若令X′_w＝(X_w,Y_w,Z_w)^T，Z_c＝s,则

$s{\stackrel{\‾}{x}}_w=H{X}_w^{\′}---\left(2\right)$

其中，H也称为单应性矩阵，且有

H＝K[r₁,r₂,t]                  (3)

每一个位置的标定板图像都独立对应一个单应性矩阵，令H＝[h₁,h₂,h₃]则可由上式推得：

[h₁ h₂ h₃]＝λK[r₁ r₂ t]        (4)

其中，λ为任意的比例因子，K摄像机内参数矩阵，因r₁和r₂是单位正交向量，即 ${r_1}^T{r}_1={r_2}^T{r}_2=1$ 且 ${r_1}^T{r}_2=0$ 则：

$h_1^T{K}^{-T}{K}^{-1}{h}_1=h_2^T{K}^{-T}{K}^{-1}{h}_2=1---\left(5\right)$

h₁K^-TK^-1h₂＝0               (6)

根据公式(5)公式(6)，计算出立体显微镜左右光路的摄像机尺度因子；针对主点坐标，采用变倍率法对摄像机的主点进行标定；假设空间中任意一点P在摄像机坐标系中的坐标为x₁,y₁,z₁，在不考虑非线性畸变以及图像坐标垂直度的情况下，该点的图像平面投影的坐标为：

u₁＝rx₁+u₀

                       (7)

v₁＝ry₁+v₀

其中u₁,v₁是该点的图像坐标，r是任一放大倍率，u₀,v₀是主点坐标；将式中的r消去，可得直线方程：

$\frac{u_1-u_0}{x_1}=\frac{v_1-v_0}{y_1}---\left(8\right)$

即在任何放大倍率下，点x₁,y₁,z₁的图像坐标都在同一直线上，且该直线一定经过摄像机主点u₀,v₀；为了求取摄像机主点u₀,v₀，选定12个在图像投影平面上不重合的点，令这12个点分别在不同放大倍率下0.7、1、2、2.4、3、5、5.8、8进行投影，利用最小二乘法对获得的12条直线进行拟合，并用最小二乘法求取12条直线的交点坐标即摄像机主点u₀,v₀；由于垂直因子γ对成像精度影响不大，因此在初始参数标定时，将γ设为0，通过优化手段得到γ具体值；针对摄像机外参数旋转矩阵R和平移矩阵t；由公式(2)知，在标定完内参数矩阵K和单应性矩阵H后可以求得摄像机旋转矩阵R＝[r₁,r₂,r₃]和平移向量t如下：

r₁＝λK^-1h₁

r₂＝λK^-1h₂                           (9)

r₃＝r₁×r₂

t＝λK^-1h₃

在初步标定完左右两个CCD摄像机的主参数后，下一步是对主参数和畸变参数的优化；为了简化优化参数，本发明采用四元数法将旋转矩阵 $R=[r_1,r_2,r_3]=\left[\begin{array}{ccc}{r}_{11}& r_{21}& r_{31}\\ r_{12}& r_{22}& r_{32}\\ r_{13}& r_{23}& r_{33}\end{array}\right]$ 中的九个未知数缩减为四个，公式如下：

$q_1=\frac{\sqrt{1+r_{11}+r_{22}+r_{33}}}{2},q_2=\frac{r_{23}-r_{32}}{4q_1},q_3=\frac{r_{31}-r_{13}}{4q_1},q_4=\frac{r_{12}-r_{21}}{4q_1}$

$q_2=\frac{\sqrt{1+r_{11}-r_{22}-r_{33}}}{2},q_1=\frac{r_{23}-r_{32}}{4q_2},q_3=\frac{r_{12}+r_{21}}{4q_2},q_4=\frac{r_{31}-r_{13}}{4q_2}$

                         (10)

$q_3=\frac{\sqrt{1-r_{11}+r_{22}-r_{33}}}{2},q_1=\frac{r_{31}-r_{13}}{4q_3},q_2=\frac{r_{12}-r_{21}}{4q_3},q_4=\frac{r_{23}-r_{32}}{4q_3}$

$q_4=\frac{\sqrt{1-r_{11}-r_{22}+r_{33}}}{2},q_1=\frac{r_{12}-r_{21}}{4q_4},q_2=\frac{r_{31}+r_{13}}{4q_4},q_3=\frac{r_{23}-r_{32}}{4q_4}$

从公式(10)计算得到新的旋转矩阵：

$R=\left[\begin{array}{ccc}-1+2q_2^2+2q_1^2& 2q_2{q}_3-2q_1{q}_4& 2q_1{q}_3+2q_2{q}_4\\ 2q_2{q}_3+2q_4{q}_1& -1+2q_1^2+2q_3^2& -2q_1{q}_2+2q_3{q}_4\\ -2q_1{q}_3+2q_4{q}_2& 2q_1{q}_2+2q_3{q}_4& -1+2q_1^2+2q_4^2\end{array}\right]---\left(11\right)$

由于四元数q₁,q₂,q₃,q₄的平方和等于一，因此旋转矩阵的待优化参数简化为三个；将已求得的左右两个CCD摄像机标定参数值作为优化初始值，利用优化算法对摄像机主参数和畸变参数进行整体优化；对左右两个CCD摄像机主参数和畸变参数的优化是基于最大似然估计准则，对于给定的n个标定板图片提供的m×n个标定点坐标，左右两个CCD摄像机参数的最优化问题通过下式的最小化问题来表达:

$S\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_i^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{[y_i-m_i(C_j,p_j,X_i)]}^2---\left(12\right)$

其中j是指参与计算的第j个摄像机，i是指第j个摄像机获取的第i个点，X_i是输入的空间点坐标，y_i是空间第i点的图像坐标，C_j是固定不变的摄像机参数向量，其长度为n₀，p_j是需要调整的摄像机参数向量，其长度为n₁，n₀+n₁即为摄像机所有参数向量长度，m_i(C_j,p_j,X_i)是摄像机的成像方程；本发明采用基于LM算法的光束平差法来解决上式的最小化优化问题；利用光束平差法的编写优化程序时除了要提供左右两个CCD摄像机的标定参数初值外还要提供一定数量的和形式的三维空间坐标点以及这些三维坐标点的相应图像坐标，因此为了精确得到整个视场空间内高精密标定板的角点坐标，采用数控机床z轴进给带动显微镜纵向移动的方法；并采用激光干涉仪(测量精度可达到0.1um)对z轴进行测量得到z轴移动的实际距离，保证所采集角点坐标点的z向实际精度；在光束平差法中，优化算法的迭代因子是由迭代参数的一阶导数组成的雅克比矩阵构成的，因此需要在计算过程中引入待成像方程中调整参数的雅克比矩阵；

$J=[\frac{\mathrm{df}}{\mathrm{dp}\left(1\right)},\frac{\mathrm{df}}{\mathrm{dp}\left(2\right)},\·\·\·,\frac{\mathrm{df}}{\mathrm{dp}\left(n_1\right)}]---\left(13\right)$

其中，f为成像方程，p(1)…p(n₁)是需要调整的摄像机参数向量元素；

光束平差法以其待调整参数的雅克比矩阵组成迭代因子，对模型参数进行调整，得到使成像误差S(θ)达到最小的模型参数最优解；

(2)待测区域特征点的提取

本发明针对待测区域特征点的提取采用的是Harris角点检测算法；Harris提出了用高斯函数代替方形区域计算梯度的算法，若目标点的像素坐标为x,y，其在x方向和y方向的位移分别为u和v，ω(u,v)为窗口函数，则该目标点的灰度变化由如下公式表示：

$\begin{array}{c}E(x,y)=\underset{u,v}{\mathrm{\Σ}}\mathrm{\ω}(u,v){[I(x+u,y+v)-I(x,y)]}^2\\ \≈\underset{u,v}{\mathrm{\Σ}}\mathrm{\ω}(u,v)(u,v)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]\left(\begin{array}{c}u\\ v\end{array}\right)\end{array}---\left(14\right)$

令 $M(x,y)=\underset{u,v}{\mathrm{\Σ}}\mathrm{\ω}(u,v)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$ 则：

$E(x,y)\≈(u,v)M(x,y)\left(\begin{array}{c}u\\ v\end{array}\right)---\left(15\right)$

计算矩阵M(x,y)的特征值λ₁,λ₂，如果两个特征值都比较大，说明目标点的图像灰度自相关函数的两个正交方向上的值均较大，则该点即为特征点；通过Harris算法对特征点进行初提取后；采用亚像素角点提取算法提取出精度更高的角点坐标；对于理想角点，它附近的像素点的灰度梯度方向均垂直于该点与理想角点的连线；这一特点可以用公式表达为：

$\▿\stackrel{\→}{H}(\stackrel{\→}{\mathrm{\α}}-\stackrel{\→}{\mathrm{\β}})=0---\left(16\right)$

其中是理想角点的灰度梯度方向，向量为图像原点指向理想角点的坐标，向量为图像原点指向理想角点附近任一边缘点的坐标；实际上，由于受到图像噪声的影响，公式(10)通常都不为零，可以将其视为误差e，即

$e=\▿\stackrel{\→}{H}(\stackrel{\→}{\mathrm{\α}}-\stackrel{\→}{\mathrm{\β}})---\left(17\right)$

在以角点为中心的邻域内，对所有点按上式计算，误差和为E,则

$E=\underset{i}{\mathrm{\Σ}}e=\underset{i}{\mathrm{\Σ}}\▿\stackrel{\→}{H}(\stackrel{\→}{\mathrm{\α}}-\stackrel{\→}{\mathrm{\β}})---\left(18\right)$

其中，i是邻域内第i个点；这样求使误差和E最小的点即是亚像素级别的角点坐标；

(3)特征点的立体匹配

采用归一化交叉算法(NCC)实现特征点的初匹配，以CCD摄像机2拍得的图片作为基准图片，并以特征点p_l为中心,构造一个图像像素N*N,的局部图像块作为模版图T，并令模版图T在待匹配图像S的包含有特征点的一定范围内进行遍历，搜索模版所覆盖的子图记做S^i,j，i,j为子图中心点在匹配图像S中的像素坐标，通过使用归一化交叉算法(NCC)来计算模版图T和搜索子图S^i,j之间的归一化互相关系数，算法如下:

$R(i,j)=\frac{\underset{m,n}{\mathrm{\Σ}}|S^{i,j}(m,n)-E\left(S^{i,j}\right)||T(m,n)-E\left(T\right)|}{\sqrt{\underset{m,n}{\mathrm{\Σ}}{(S^{i,j}(m,n)-E\left(S^{i,j}\right))}^2\underset{m,n}{\mathrm{\Σ}}{(T(m,n)-E\left(T\right))}^2}}---\left(19\right)$

其中：m,n代表各个像素的坐标，T(m,n)为点m,n处的灰度值，

E(S^i,j)和E(T)分别是搜索子图S^i,j与模版图T的灰度平均值；互相关系数R(i,j)的值越大，则搜索子图与模板图的匹配程度越高；设定互相关系数R(i,j)的阈值后，选出特征点p_l的候选匹配点；使用NCC算法进行初匹配后还需要利用外极限约束和距离约束来对匹配的特征点对进行校正；通过使用Longguet-Higgins提出的归一化8点算法在摄像机标定时计算左右两个CCD摄像机的基本矩阵F，根据外极线理论，对基准图像特征点p_l，在待匹配图像S中对应的外极线l_r可以表示为：

l_r＝F·p_l              (20)

初匹配的候选特征点若在外极线l_r附近则可进一步确定该候选特征点为匹配点；最后统计每一候选特征点到其他特征点的距离和基准图像中特征点到其他特征点点距离相等的个数num；num最大值时对应的候选特征点则为基准图像中特征点的正确匹配点；

(4)三维坐标的求取

首先计算出双目图像中特征点经亚像素算法提取出的图像坐标值，并大致估算出特征点在世界坐标系中的估计值；在优化程序中将左右两个CCD摄像机参数对应雅克比矩阵设为零，同时在优化程序计算迭代因子部分中嵌入关于空间坐标x、y、z的雅克比矩阵方程；将已经优化好的左右两个CCD摄像机参数以及特征点点图像坐标值和世界坐标估计值输入光束平差法优化程序，就可以对空间点的三维坐标进行优化，得到特征点的三维空间测量值，从而完成三维坐标点的求取并最终实现特征点的定位。