CN101936761A - 大型料场料堆视觉测量方法 - Google Patents

Abstract

本发明涉及的是一种非接触测量方法，具体是一种大型料场料堆视觉测量方法。一种大型料场料堆视觉测量方法，包括以下步骤：1)构建料场视觉测量***，2)安装双目视觉子***，3)建立CCD摄像机几何模型，4)CCD摄像机参数标定，5)图像采集，6)图像处理，7)角点的寻找与匹配，8)物点三维坐标值的确定，9)在世界坐标系中进行曲面插值并建立高度图，10)体积测量。本发明的视觉测量方法自动消除了广角镜头引起的图像畸变，能对左、右视图自动进行角点的精确匹配，利用图像数据变换关系自动生成物点的三维空间坐标，优化立体角点分布、构建世界坐标系高度图，使得几何不规则料堆也能简捷、直观地实现其体积的精确计算。

Description

大型料场料堆视觉测量方法

技术领域

本发明涉及的是一种非接触测量方法，具体是一种大型料场料堆视觉测量方法。

背景技术

至今为止，钢铁、化工、电厂等企业所拥有的大型原料场的堆取料操作与进出料规划，基本上还是依靠以人工操作为主的管理与施工模式，因此极大程度地影响了这些企业的整体自动化和信息化技术水平的提升。

由于此类料场中所堆放的原料均为散装物质，如铁矿沙、焦炭、石灰等，在生产过程中，其体积量始终处于一种动态变化状态，一般情况下，企业对这些原料的库存量丈量及其动态管理均依靠人工测量的方法进行，效率低下，每隔一定时间还要停机进行人工测量影响企业生产的连续性。因此，要实现料场的自动化管理势必首先实现料场的数字化技术，而料场数字化技术的关键点又在于如何实现料场体积的自动丈量和动态观测。在现有技术中已经开始采用激光测距技术实现对料场体积的自动丈量和动态观测，由于激光技术在该领域的应用存在装置复杂、对工作环境要求苛刻、价格高昂等条件限制，使得该技术的推广有着很大的局限性。基于机器视觉的测量技术是当前的一种新兴技术，但是由于其技术尚不完善，加之在料场数字化***中的应用尚无先例，因此没有现成的成熟技术可以借鉴。

经对现有技术文献的检索发现，中国发明专利号200510026200.5公开了一种料场测量的计算机视觉***自动测量的方法，该发明对选定的原料料堆设计CCD摄像机布局，CCD摄像机所采集到的图像信号经由视频电缆与图像采集卡输至图像处理前端；图像处理前端中图像处理软件对图像信号实时进行背景淡化与边缘提取处理，完成图像的数字化信息表达；图像处理前端数字化表达的信息通过IEEE1394协议接口与上位计算机连接，实现图像处理前端与上位计算机的信息交互；上位计算机对图像的数字化表达信息进行再加工，求取公共角点重构对象的三维图像，与此同时运算输出料堆的空间几何尺度测量结果。可以肯定，该发明技术方法在实现料场的计算机自动测量过程是有效的，但是，该技术过于依靠“CCD摄像机布局，满足作为高位视觉采集点，采集料堆全景图像”进行对料场图像信息的采集；同时，该发明尚不具备对图像畸变等影响因素自动消除能力。

发明内容

本发明的目的在于提供一种大型料场料堆视觉测量方法，该测量方法基于机器视觉的料场体积测量方法，利用料场现有的堆取料机行走机构，在堆取料机作业过程中自动完成料场图像采集与处理，最终实现料堆体积的精确计算。

一种大型料场料堆视觉测量方法，包括以下步骤：

步骤一、构建料场视觉测量***，由第一双目视觉子***、第二双目视觉子***和服务器三部分组成，服务器的图像信息处理模块接收第一双目视觉子***、第二双目视觉子***的信号后进行处理和运算后输出最终对图像信息的识别、理解和测量结果；

步骤二、安装双目视觉子***，将第一双目视觉子***安装于堆料机靠近料堆一侧机架上，将第二双目视觉子***安装于取料机靠近料堆一侧机架上，第一和第二双目视觉子***中的两个CCD摄像机的安装应该遵循：两个摄像机的光心处于同一水平面并且两个光轴与两光心连线形成一个等腰三角形；

步骤三、建立CCD摄像机几何模型，包括：世界坐标系与摄像机坐标系的变换、理想透视投影变换、归一化虚平面图像坐标系中图像畸变矫正数学模型、归一化虚平面图像坐标系畸变点坐标到计算机图像坐标系的转换，最后获得从世界坐标系到计算机图像坐标系的映射关系方程；

步骤四、CCD摄像机参数标定，包括内部参数和外部参数，采用分步标定方法，先标定图像的纵横比和主点坐标，再利用径向平行约束，求解大部分的外部参数，最后再引入畸变模型，线性求解内部参数和平移量；

步骤五、图像采集，两个双目视觉子***分别随着堆取料机的行走在线采集料堆及料场图像；每一个双目视觉子***在每一个采样时刻能够通过其中的两个CCD摄像机同时采集两幅图像，在前端信号处理器中保存图像序列，并将该两幅图像经过预处理后，同时接受由远程控制程序控制指令要求堆取料机上的前端信号处理器向上位服务器传送数字图像；服务器在每一个采样时刻能够同时接收到四幅图像，包括来自第一、第二双目视觉子***的各两幅图像，每两幅图像又分左、右两幅图像，在整个采样期间，所有采集到的图像分别形成第一、第二双目视觉子***的两组左右视图对序列；

步骤六、图像处理，消除图像畸变，由计算机图像坐标系像素点坐标(u，v)通过归一化虚平面图像坐标系畸变点坐标到计算机图像坐标系的转换方程计算求取对应点(u，v)在归一化虚平面图像坐标系畸变点坐标(x_d，y_d)，再将(x_d，y_d)的值输至代入归一化虚平面图像坐标系中图像畸变矫正数学模型进行求逆运算，因此求得矫正畸变后的归一化虚平面图像坐标系理想点坐标(x_u，y_u)，再由理想点坐标(x_u，y_u)替代(x_d，y_d)的值代入归一化虚平面图像坐标系中图像畸变矫正数学模型求出对应点(u，v)新的坐标值，即新的像素点位置；计算机图像坐标系中的所有像素点坐标(u，v)逐一经过上述计算过程后，就能够获得像素点在画面上的理想排列，即真实景物图像的恢复，或者说，获得一幅反映真实景物的理想图像；

步骤七、角点的寻找与匹配，在采样图像中寻找若干个角点作为代表最具特征的边缘点，采用基于图像灰度的角点检测算法，主要通过计算曲率及梯度来达到检测角点的目的；针对已经寻找出来的角点在外极线约束的指导下在左、右视图中寻找与之最为匹配的点；实际操作过程需要从左、右视图的数字图像中取得足够多的立体角点；

步骤八、物点三维坐标值的确定，在左、右视图中相互匹配的角点坐标(u_l，v_l)、(u_r，v_r)通过从世界坐标系到计算机图像坐标系的映射关系方程能够转换为对应的世界坐标(X_Wl，Y_Wl，Z_Wl)和(X_Wr，Y_Wr，Z_Wr)，采用求取两条直线公垂线中值的方法最后获得物点空间坐标的最优近似；

步骤九、在世界坐标系中进行曲面插值并建立高度图，对图像在全黑背景之上显示的有限个离散像素点，必须在相邻离散点之间采用贝赛尔曲面插值算法进行插值以形成一个接近实际的连续高度图；这个高度图作为图像处理***的主要结果，包含了关于料场的全部信息以便用于后续的深层分析和数据库的数据存储；

步骤十、体积测量，首先，定义“空地”，即未堆料的地面高度，从头至尾寻找空地，根据空地连续出现的位置判断应该将料场划分为几个料堆，并排除很小的起伏造成的伪料堆和杂讯，得到每一个料堆的起始和终止的位置；其次，从料场高度图中截取出料堆高度个图；令高度图中某一个像素的灰度值所代表的实际高度为h_i，该处料堆的体积可以看成由有限个高度为h_i、底面积为δ_i的微长方体累加而成(i为参与计算的像素点，i＝1，2，K，n，n为参与计算的像素点总数)，因此，所有几何不规则的料堆体积V均可按照如下公式计算：

所述的构建料场视觉测量***中每个双目视觉子***，包括：两个CCD摄像机、图像采集卡、前端信号处理器、两个光电转换器和光纤；左、右CCD摄像机在采样时刻将所采集到的左、右视图并行输入至图像采集卡；图像采集卡将原始图像转换成数字图像信号并输至前端信号处理器；前端信号处理器对数字图像进行去除噪声和强化预处理；经过预处理后的数字图像信号经第一个光电转换器将电信号转换成光信号，光信号通过光纤传输到第二个光电转换器；第二个光电转换器再将图像光信号转换成电信号，最后，电信号通过服务器的输入接口输入至服务器的图像信息处理模块；第一、第二双目视觉子***分别将各自采集到的左、右视图模拟信号并行输至本子***中的图像采集卡，图像采集卡将图像模拟信号转换成数字图像信号并输至前端信号处理器，前端信号处理器对数字图像进行去除噪声和强化预处理，经过预处理后的数字图像信号经第一个光电转换器将电信号转换成光信号，光信号通过光纤传输到第二个光电转换器，第二个光电转换器再将图像光信号转换成电信号，最后，电信号通过服务器的第一输入接口输入至服务器的图像信息处理模块，服务器将数字图像信号进行处理和运算后输出最终对图像信息的识别、理解和测量结果。

所述的CCD摄像机参数标定，具体如下：

CCD摄像机需要标定的外部参数：包括R和t，即

$\left[\begin{array}{cc}R& t\end{array}\right]=\left[\begin{array}{cccc}{r}_1& r_2& r_3& t_x\\ r_4& r_5& r_6& t_y\\ r_7& r_8& r_9& t_z\end{array}\right]$

共12个，但因R为单位正交矩阵，必须满足6个正交约束，故实际上只有6个外部参数需要标定；

内部参数：包括f、k₁、s_x和c_x、c_y，共5个；

首先确定图像的纵横比，垂直拍摄一个圆环，然后计算其在水平方向和垂直方向上的像素直径比即为纵横比s_x；

其次用相机的不同焦距拍摄同一景物，然后计算其缩放中心得到主点坐标(c_x，c_y)；

最后通过设定过渡参数求解R和t_x、t_y，然后根据R的正交性计算得到有效焦距f、径向畸变系数k₁和平移量t_z。

所述在世界坐标系中进行曲面插值并建立高度图，具体如下：

对图像在全黑背景之上显示足够多的离散像素点，如1400个，在相邻离散点之间采用贝赛尔曲面插值算法进行插值以形成一个接近实际的连续高度图，该高度图作为图像处理***的主要结果，包含了关于料场的全部信息以便用于后续的深层分析和数据库的数据存储，高度图是一个平面灰度图像，水平方向和世界坐标系的X轴方向一致，竖直方向和Z轴方向一致，某个点的像素灰度(0～255)和该点的高度(Y分量)成正比，将最低点0米对应于0(全黑)，最高点16米作为255(全白)，中间的高度对应于相应的灰色；

每一个像素代表一个面积，对应整幅图像就代表着一个特定长×宽的原料场，如：像素代表0.5×0.5m²，一个1440×96pixel的图像就代表720m×48m的一个完整料场，其中的任何一点对应着相应的位置高度；

数字图像以BMP格式进行存放以供下一步的分析，以及用于数据库管理***的信息交互。

所述体积测量，具体如下：

①定义：高度低于0.5米的被定义为“空地”，即未堆料；首先从头至尾寻找空地，根据空地连续出现的位置判断应该将料场划分为几个料堆，同时排除很小的起伏造成的伪料堆和杂讯，得到每一个料堆的起始和终止的位置；

②从料场高度图中截取出包含要求计算料堆的矩形个图，令高度图中某一个像素的灰度值所代表的实际高度为h_i、底面积为δ_i，i为参与计算的像素点，i＝1，2，K，n，n为参与计算的像素点总数，因此，所有几何不规则的料堆体积V均可按照公式计算，当δ_i＝const＝0.25m²时，还可以直接表达为

当定义0.5m以下为“空地”时，以包含料堆的矩形底面区域进行体积运算，公式还可以表达为：

式中，h_kl为第k列、第l行的像素点对应高度，k＝1，2，K，n，n为像素总列数，l＝1，2，K，m，m为像素总行数，dx、dy分别为像素点行向与列向的间隔所对应的实际距离，单位为m。

本发明的视觉测量方法采用料场网格高度图的方法建立料堆体积计算公式并在定义“空地”的条件下，将运算进一步简洁表达为

实现料堆和料场体积的精确计算；自动消除了广角镜头引起的图像畸变，能对左、右视图自动进行角点的精确匹配，利用图像数据变换关系自动生成物点的三维空间坐标，优化立体角点分布、构建世界坐标系高度图，使得几何不规则料堆也能简捷、直观地实现其体积的精确计算。

附图说明

图1为摄像机几何模型图；

图2为外极几何图；

图3为双目视觉***中左、右CCD摄像头水平安装图；

图4为左、右CCD摄像机在堆、取料机上的工作状态图。

图中：11第一左CCD摄像机、12第一右CCD摄像机、21第二左CCD摄像机、22第二右CCD摄像机。

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

如图1所示，实施例中采用的摄像机几何模型：

①计算机图像坐标系o_uvuv：原点o_uv位于CCD图像平面的左上角，u和v分别表示像素位于数组的列数和行数，单位为pixel(像素)。

②物理图像坐标系oxy：原点(即主点o)定义在摄像机光轴与物理图像平面的交点，x、y轴分别与u、v轴平行。p_u(x_u，y_u)表示针孔模型下P点的理想图像坐标，p_d(x_d，y_d)是由透镜径向畸变引起的偏离p_u(x_u，y_u)的实际图像坐标。

③摄像机坐标系o_cx_cy_cz_c：原点o_c定义在摄像机的光心，x_c、y_c轴分别平行于x、y轴，z_c轴与光轴重合。o_co为摄像机的有效焦距f，(x_c，y_c，z_c)表示物点P在摄像机坐标系下的三维坐标。

④世界坐标系O_WX_WY_WZ_W：根据具体环境来定，(X_W，Y_W，Z_W)表示物点P在世界坐标系下的三维坐标。

如图2所示，所述的外极几何约束，即针对双目视觉***从两个角度观测同一个点，寻求两个成像点之间的外极线约束的相关几何学问题，被称为外极几何。图中，左、右摄像机交叉摆放，P为左摄像机和右摄像机观测的同一个三维点，O_l、O_r分别为左、右摄像机的光心；P点在左、右摄像机的归一化虚成像平面上的成像点分别为p_l[x_l y_l 1]^T和p_r[x_r y_r 1]^T；e_l、e_r分别为左、右像面上的外极中心，对于共同观测点P的左、右两条外极线分别通过点P_l、e_l和P_r、e_r。

本实施例包括如下具体步骤：

步骤一、构建料场视觉测量***

由第一双目视觉子***、第二双目视觉子***和服务器三部分组成料场视觉测量***。每个双目视觉子***，包括：两个CCD摄像机、图像采集卡、前端信号处理器、两个光电转换器和光纤。左、右CCD摄像机在采样时刻将所采集到的左、右视图并行输入至图像采集卡；图像采集卡将原始图像转换成数字图像信号并输至前端信号处理器；前端信号处理器对数字图像进行去除噪声和强化等预处理；经过预处理后的数字图像信号经第一个光电转换器将电信号转换成光信号，光信号通过光纤传输到第二个光电转换器；第二个光电转换器再将图像光信号转换成电信号，最后，电信号通过服务器的输入接口输入至服务器的图像信息处理模块。

第一、第二双目视觉子***分别将各自采集到的左、右视图模拟信号并行输至本子***中的图像采集卡，图像采集卡将图像模拟信号转换成数字图像信号并输至前端信号处理器，前端信号处理器对数字图像进行去除噪声和强化等预处理，经过预处理后的数字图像信号经第一个光电转换器将电信号转换成光信号，光信号通过光纤传输到第二个光电转换器，第二个光电转换器再将图像光信号转换成电信号，最后，电信号通过服务器的第一输入接口输入至服务器的图像信息处理模块，服务器将数字图像信号进行处理和运算后输出最终对图像信息的识别、理解和测量结果。

步骤二、安装双目视觉子***

如图3所示，第一双目视觉子***安装于堆料机靠近料堆一侧机架上，第二双目视觉子***安装于取料机靠近料堆一侧机架上，第一和第二双目视觉子***中的两个CCD摄像机的安装应该遵循：两个摄像机的光心处于同一水平面原则和两个光轴与两光心连线形成一个等腰三角形原则。

如图4所示，第一双目视觉子***中的左、右CCD摄像机11、12安装在堆料机上，两者之间的几何关系确定方法如下：

①第一左CCD摄像机11与第一右CCD摄像机12之间的水平间距，即两者的摄像机坐标系原点之间距离：根据堆料机平台最大可延伸宽度大约在20m左右，因此可以确定两个镜头水平间距为20m。

②镜头垂直方向角：根据料堆为16m高、像头的垂直位置离料脚的纵向距离约为4m、料堆的安息角为26°、堆料机的固定平台最高可利用高度约为6m，镜头只能安装于6m高度的平台上，垂直视场的上边缘线与料场斜坡面平行，垂直视场的下边缘线与料场底线的交点距离镜头垂直点大约在2.926m左右，因此可以预算出镜头的下俯视角为19°。

③镜头水平方向角度：按照立体视觉的一般规律，双目视觉器中的两个摄像机和所测量的物点应该组成一个等腰三角形的原则，根据上述水平间距和镜头垂直方向角，堆料机上镜头水平方向角度为32°。

第二双目视觉子***中的左、右CCD摄像机21、22安装在取料机上，两者之间的几何关系确定方法如下：

①第二左CCD摄像机21与第二右CCD摄像机22之间的水平间距，即两者的摄像机坐标系原点之间距离：根据取料机平台最大可延伸宽度大约在10m左右，受其所限，只能确定两个镜头水平间距为10m。

②镜头垂直方向角与堆料机上的一样。

③镜头水平方向角度：按照立体视觉的一般规律，双目视觉器中的两个摄像机和所测量的物点应该组成一个等腰三角形的原则，根据上述水平间距和镜头垂直方向角，取取料机上镜头水平方向角度为17.35°。

步骤三、建立CCD摄像机几何模型

全部几何模型包括：

1)世界坐标系与摄像机坐标系的变换，用旋转矩阵R和平移向量t来描述：

$\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]=R\left[\begin{array}{c}{X}_W\\ Y_W\\ Z_W\end{array}\right]+t;$ $R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right],$ $t=\left[\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right]$ (公式1)

[X_W Y_W Z_W]^T为世界坐标系点坐标；[x_c y_c z_c]^T为摄像机坐标系点坐标。

2)理想透视投影变换，即点坐标从摄像机坐标系到归一化虚平面图像坐标系的变换：

x_u＝fx_c/z_c，y_u＝fy_c/z_c            (公式2)

式中，f为摄像机的成像焦距，(x_u，y_u)为归一化虚平面图像坐标系理想点坐标。

3)归一化虚平面图像坐标系中图像畸变矫正数学模型，仅考虑一阶径向畸变时的畸变方程为：

x_d＝(1+k₁r²)x_u，y_d＝(1+k₁r²)y_u    (公式3)

式中，

k1为一阶径向畸变系数。

4)归一化虚平面图像坐标系畸变点坐标到计算机图像坐标系的转换：

u＝x_d/dx+c_x，v＝y_d/dy+c_y，s_x＝dy/dx    (公式4)

式中，(u，v)为计算机图像坐标系像素点坐标，(c_x，c_y)为主点o的像素坐标即主点坐标，(dx，dy)分别为图像平面上x、y方向上单位像素间的距离，sx为图像纵横比。

由公式1～4得

$\left[\begin{array}{c}p\\ 1\end{array}\right]=\frac{1}{s}K\left[\begin{array}{cc}R& t\end{array}\right]\left[\begin{array}{c}P\\ 1\end{array}\right]$ (公式5)

公式5称为世界坐标系到计算机图像坐标系的映射关系方程，其中，p＝[u v]^T，P＝[X_W Y_W Z_W]^T，s为世界坐标系中空间点映射到CCD摄像机坐标系(x_c，y_c，z_c)中z_c轴上的分量，其数值等于上述映射关系方程式右边计算结果所得到的第三列向量中的第三元素值；

为三维点坐标从归一化虚成像平面到物理成像平面的等比例缩放矩阵，其缩放的比例(逻辑焦距)和实际焦距f有关，单位为“像素/m”，即

$f_x=\frac{f(1+k_1{r}_d^2)}{\mathrm{dx}};$ $f_y=\frac{f(1+k_1{r}_d^2)}{\mathrm{dy}}$ (公式6)

这表示，虽然只有一个物理焦距，但是表现在图像上可能有2个不同的逻辑焦距。

所述归一化虚平面图像坐标系，是指：沿光轴在与摄像机光学成像平面形成对称的平面称为虚平面，在其上所建立的坐标系称为虚平面图像坐标系，将图像像素特征量进行归一化后在虚平面坐标系上的表达称为图像的归一化虚平面坐标系表达，使用归一化表达的虚平面坐标系所处平面，又称为归一化虚平面。

步骤四、CCD摄像机参数标定

CCD摄像机需要标定的参数包括内部参数和外部参数。

外部参数：包括R和t，即

$\left[\begin{array}{cc}R& t\end{array}\right]=\left[\begin{array}{cccc}{r}_1& r_2& r_3& t_x\\ r_4& r_5& r_6& t_y\\ r_7& r_8& r_9& t_z\end{array}\right]$

共12个，但因R为单位正交矩阵，必须满足6个正交约束，故实际上只有6个外部参数需要标定。

内部参数：包括f、k₁、s_x和c_x、c_y，共5个。

标定过程采用分步标定方法，先标定图像的纵横比s_x和主点坐标(c_x，c_y)，再利用径向平行约束，求解大部分的外部参数，最后再引入畸变模型，线性求解内部参数f、k₁和平移量t_z，具体步骤如下：

首先、图像纵横比的确定

由于摄像机焦距同时在x、y方向上缩放图像，因此，垂直拍摄一个圆环，然后计算其在水平方向和垂直方向上的像素直径比即为纵横比s_x。

其次、主点坐标的标定

当摄像机的有效焦距变化时，相应的图像发生缩放，而光轴不变，故光轴与图像平面的交点(主点)不变。因此，可以用相机的不同焦距拍摄同一景物，然后计算其缩放中心即得主点坐标。

假定摄像机的焦距由f₁变为f₂，由公式4得：

$\frac{x_2-c_x}{x_1-c_x}=\frac{y_2-c_y}{y_1-c_y}\⇒\left[\begin{array}{cc}{y}_1-y_2& x_2-x_1\end{array}\right]\×{\left[\begin{array}{cc}{c}_x& c_y\end{array}\right]}^T=x_2{y}_1-x_1{y}_2$ (公式7)

式中，(x₁，y₁)为在有效焦距f₁下某特征点的像素坐标，(x₂，y₂)为同一特征点在有效焦距f₂下的像素坐标。选择多个特征点，利用最小二乘法即可线性解出(c_x，c_y)。

再是、旋转矩阵R和平移分量t_x、t_y的求解，具体过程如下：

①先求取过渡参数

由于我们使用的是共面点，所以Z_W＝0，由径向平行约束有

$\frac{x_d}{y_d}=\frac{x_c}{y_c}=\frac{r_1{X}_W+r_2{Y}_W+t_x}{r_4{X}_W+r_5{Y}_W+t_y};$ $x_d=\frac{\mathrm{dy}(u-c_x)}{s_x},$ y_d＝dy(v-c_y)  (公式8)

考虑到透镜变形对相机有效焦距和到标定平面距离的影响，故对上式右边分子分母同除以t_y，并整理可得到以下矩阵形式的方程：

$\left[\begin{array}{ccccc}{y}_d{X}_W& y_d{Y}_W& y_d& {-x}_d{X}_W& {-x}_d{Y}_W\end{array}\right]\×{\left[\begin{array}{ccccc}\frac{r_1}{t_y}& \frac{r_2}{t_y}& \frac{t_x}{t_y}& \frac{r_4}{t_y}& \frac{r_5}{t_y}\end{array}\right]}^T=x_d$ (公式9)

式中，列向量中的5个元素为未知的，对于每一个标定点，如果已知其世界坐标和相应的图像坐标，就可以列出如上的一个方程，取多个标定点(大于5)，用最小二乘法求该超定方程，即得过渡参数：

l₁＝r₁/t_y、l₂＝r₂/t_y、l₃＝t_x/t_y、l₄＝r₄/t_y、l₅＝r₅/t_y         (公式10)

②根据过渡参数求解R和t_x、t_y

求出上述过渡参数后，就可以根据R的正交性，求出t_y的大小，当参数l1、l2、l₃、l₄、l₅不两两同时为0时，

$t_y^2=[A-{(A^2-4B)}^{1/2}]/2B;$ $A=l_1^2+l_2^2+l_4^2+l_5^2,B={(l_1{l}_5-l_2{l}_4)}^2$ (公式11)

否则为其余两个参数t_x与t_z的平方和的倒数。

判断t_y的符号，在运用径向平行约束时，包括两平行向量同向和反向的两种情况，故t_y的符号有两种可能。可采用下列方法判断：

先设t_y为正，由过渡参数求出r₁、r₂、r₄、r₅、t_x，由此将标定点再投影到图像平面上，计算出对应点的图像坐标x_c、y_c，如果x_d与x_c同号，且y_d与y_c同号，那么t_y为正，否则t_y为负，则上面求出的其余外部参数都作相应的改变。

根据R的正交性，计算R的其余参数

$r_3=\sqrt{1-{r_1}^2-r_2^2},$ $r_6=\sqrt{1-r_4^2-r_5^2},$ r₇＝r₂r₆-r₃r₅、

r₈＝r₃r₄-r₁r₆、r₉＝r₁r₅-r₂r₄                        (公式12)

其中，如果r₁r₄+r₂r₅的符号为正，则在r₆前加负号；如果下面计算的f为负号，则r₃、r₆、r₇、r₈的符号与上面相反。

有效焦距f，径向畸变系数k₁和平移量t_z的求解，由公式3得：

$x_u=f\frac{x_c}{z_c}=\frac{x_d}{1+k_1{r}_d^2}=f\frac{r_1{X}_W+r_2{Y}_W+t_x}{r_7{X}_W+r_8{Y}_W+t_z},$

$y_u=f\frac{y_c}{z_c}=\frac{y_d}{1+k_1{r}_d^2}=f\frac{r_4{X}_W+r_5{Y}_W+t_y}{r_7{X}_W+r_8{Y}_W+t_z}$ (公式13)

整理后可得如下含未知数f、t_z、k的线性方程组：

$\left[\begin{array}{ccc}E& {\mathrm{Er}}_d^2& {-x}_d\end{array}\right]\×{\left[\begin{array}{ccc}f& k& t_z\end{array}\right]}^T=x_{d}G,$

$\left[\begin{array}{ccc}F& {\mathrm{Fr}}_d^2& {-y}_d\end{array}\right]\×{\left[\begin{array}{ccc}f& k& t_z\end{array}\right]}^T=y_{d}G$ (公式14)

式中，E＝r₁X_W+r₂Y_W+t_x、F＝r₄X_W+r₅Y_W+t_y、G＝r₇X_W+r₈Y_W、k＝fk₁，用最小二乘法联合解这个线性方程组，即可得f、t_z、k的解。进一步判断f的符号，如果f为负，则将f反号。进一步解得k₁＝k/f，其余参数符号按上面所述作相应的改变。

至此，通过分步标定，采用线性方法解出了摄像机的全部内外参数，因此避免了非线性优化的繁琐和不稳定。

步骤五、图像采集

双目视觉子***采用120°广角光学镜头，两个双目视觉子***分别随着堆取料机的行走在线采集料堆及料场图像。每一个双目视觉子***在每一个采样时刻能够通过其中的两个CCD摄像机同时采集两幅图像，按照1帧/米的速度在堆取料机上的信号处理器中保存图像序列，并将该两幅图像经过预处理后，同时接受由远程控制程序控制指令要求堆取料机上的前端信号处理器向上位服务器传送数字图像。服务器在每一个采样时刻能够同时接收到四幅图像，包括来自第一、第二双目视觉子***的各两幅图像，每两幅图像又分左、右两幅图像，在整个采样期间，所有采集到的图像分别形成第一、第二双目视觉子***的两组左右视图对序列。

目前通信光缆对无压缩全幅(即768*576)彩色数字图像的传输速度可以达到5帧/秒，堆取料机的行进速度为1米/秒，因此，***设定每隔1米距离采集1幅图像完全能够满足通信光缆对图像网络传输速率的要求。

现场图像以BMP的格式保存在文件中，供后续处理和数据库子***使用。按照左镜头、右镜头分别存储在两个文件夹下，如df、dr等，分别表示堆料机左、右图像文件夹；文件名即为地址数据，单位为米，如201.bmp、202.bmp等，分别表示201、202米处采集到的图像。图像分辨率和摄像机分辨率一致，对PAL制式的摄像机为768*576。

步骤六、图像处理

消除图像畸变，由计算机图像坐标系像素点坐标(u，v)通过公式4计算求取对应点(u，v)在归一化虚平面图像坐标系畸变点坐标(x_d，y_d)，再将(x_d，y_d)的值输至代入公式3进行求逆运算，因此求得矫正畸变后的归一化虚平面图像坐标系理想点坐标(x_u，y_u)，再由理想点坐标(x_u，y_u)替代(x_d，y_d)的值代入公式4求出对应点(u，v)新的坐标值，即新的像素点位置。

计算机图像坐标系中的所有像素点坐标(u，v)逐一经过上述计算过程后，就能够获得像素点在画面上的理想排列，即真实景物图像的恢复，或者说，获得一幅反映真实景物的理想图像。

步骤七、角点的寻找与匹配

在采样图像中寻找若干个角点作为代表最具特征的边缘点，采用基于图像灰度的角点检测算法，主要通过计算曲率及梯度来达到检测角点的目的。具体实现为：确定图像中的一个局部窗口，通过窗口在图像中的移动，确定图像的灰度变化值。在移动的过程中有三种可能：如果窗口处在灰度均匀的地方，则计算结果是一个较小的值；如果窗口跨越了边缘，那么在沿边缘线运动的结果是一个较小的值，而在垂直于边缘的地方将是一个较大的灰度变化；如果窗口处在角点或者是独立的点上，那么在该点的每个方向上都是灰度变化较大的地方，角点可以迅速被检测出来。

针对已经寻找出来的角点在外极线约束的指导下在左、右视图中寻找与之最为匹配的点。采用“外极线”约束对所获空间点在两幅图像上的投影点是否匹配进行判断，即空间点P在图像上的投影点P_l、P_r必须在空间任意一点P与左右透视中心O_l、O_r构成一个平面与两个成像平面的交线O_lP、O_rP上。

对于交叉摆放的CCD，O_l、O_r连线分别与左、右视图像面相交于e_l、e_r，即左。右像面上的外极中心；左、右两条外极线分别通过点P_l、e_l和P_r、e_r。由立体视觉理论可知，左、右视图像面上任意一对对应点必然在它们相应的外极线上。

令左、右摄像机坐标系分别为(x_l，y_l，z_l)和(x_r，y_r，z_r)，左、右摄像机坐标系的变换关系为：

$\left[\begin{array}{c}{x}_l\\ y_l\\ z_l\end{array}\right]=\left[\begin{array}{ccc}{r}_{\mathrm{rl}1}& r_{\mathrm{rl}2}& r_{\mathrm{rl}3}\\ r_{\mathrm{rl}4}& r_{\mathrm{rl}5}& r_{\mathrm{rl}6}\\ r_{\mathrm{rl}7}& r_{\mathrm{rl}8}& r_{\mathrm{rl}9}\end{array}\right]\left[\begin{array}{c}{x}_r\\ y_r\\ z_r\end{array}\right]+\left[\begin{array}{c}{t}_{\mathrm{rlx}}\\ t_{\mathrm{rly}}\\ t_{\mathrm{rlz}}\end{array}\right];$ $\left[\begin{array}{c}{x}_r\\ y_r\\ z_r\end{array}\right]=\left[\begin{array}{ccc}{r}_{\mathrm{lr}1}& r_{\mathrm{lr}2}& r_{\mathrm{lr}3}\\ r_{\mathrm{lr}4}& r_{\mathrm{lr}5}& r_{\mathrm{lr}6}\\ r_{\mathrm{lr}7}& r_{\mathrm{lr}8}& r_{\mathrm{lr}9}\end{array}\right]\left[\begin{array}{c}{x}_l\\ y_l\\ z_l\end{array}\right]+\left[\begin{array}{c}{t}_{\mathrm{lrx}}\\ t_{\mathrm{lry}}\\ t_{\mathrm{lrz}}\end{array}\right]$ (公式15)

将左、右摄像机坐标系各自转换为世界坐标系表达式后，联立方程，即可求取公式15中的所有参数值。

左摄像机透视中心O_l在右像面的投影点为右外极中心同理右摄像机透视中心O_r在左像面的投影点为左外极中心

左像面任意一点

在右像面的外极线

上，同理右像面任意一点

在左像面的外极线

上，左、右像面上所有的外极线均通过左、右外极中心，因此可根据左右摄像机坐标对应关系求解外极线方程。设左像面上外极线斜率为k_l，并过e_l点；在右像面上的外极线斜率为k_r，并过e_r点，推导出k_l到k_r的变换关系和k_r到k_l的变换关系：

$k_r=-\frac{t_{\mathrm{lrz}}(r_{\mathrm{lr}4}+r_{\mathrm{lr}5}{k}_l)-t_{\mathrm{lry}}(r_{\mathrm{lr}7}+r_{\mathrm{lr}8}{k}_l)}{t_{\mathrm{lrx}}(r_{\mathrm{lr}7}+r_{\mathrm{lr}8}{k}_l)-t_{\mathrm{lrz}}(r_{\mathrm{lr}1}+r_{\mathrm{lr}2}{k}_l)};$

$k_l=-\frac{t_{\mathrm{rlz}}(r_{\mathrm{rl}4}+r_{\mathrm{rl}5}{k}_r)-t_{\mathrm{rly}}(r_{\mathrm{rl}7}+r_{\mathrm{rl}8}{k}_r)}{t_{\mathrm{rlx}}(r_{\mathrm{rl}7}+r_{\mathrm{rl}8}{k}_r)-t_{\mathrm{rlz}}(r_{\mathrm{rl}1}+r_{\mathrm{rl}2}{k}_r)}$ (公式16)

对于求解后的曲线，可利用外极特性沿曲线寻找对应点。假定左视图曲线图像为参考图像，则右视图为待匹配图像，首先求出沿左视图某标号曲线的任一点

与左外极中心e_l连线的斜率k_l，再求

在右视图外极线的斜率k_r，由斜率k_r和右外极线中心e_r求出外极线。再利用模糊相似度方法实现左、右两视图中角点的精确匹配，图8示出在左、右两视图中寻找角点并实现匹配的效果，图中红点为正确线上匹配点，黄点为错误的全局匹配点。

实际操作过程需要从左、右视图的数字图像中取得足够多的立体角点，如：在700米长的料场中，共建立了1400个立体角点(每米堆料机和取料机各一个点)。这些立体角点的坐标是基于统一的料场世界坐标系，通过前述的方法得到立体坐标点，并且根据行程定位编码器提供的料机平移坐标进行校正，规范到统一的料场世界坐标系。料场世界坐标系的X_W轴垂直于地面，Y_W轴平行于地面，Z_W轴垂直于轨道，坐标系原点在料场平面图的左下角。

第八步、物点三维坐标值的确定

将左、右视图中相互匹配的角点坐标(u_l，v_l)、(u_r，v_r)转换为对应的世界坐标(X_Wl，Y_Wl，Z_Wl)和(X_Wr，Y_Wr，Z_Wr)。

由于摄像机成像模型并非理想的小孔成像以及CCD成像时噪音的影响和标定参数的微小误差，都会使得左、右视图中的物点换算坐标与透视中心的连接直线未必相交于空间一点，而这两条直线又不可能平行，所以这两条直线势必是异面直线。换句话说，一般情况下，(X_Wl，Y_Wl，Z_Wl)与(X_Wr，Y_Wr，Z_Wr)不会完全重合，因此需要采用在上述两条直线的公垂线上求取中值的方法最后获得物点空间坐标的最优近似。

在几何上，使用异面直线公垂线的中点逼近空间点的误差较小。具体计算中，首先计算两条异面直线的最短距离，即公垂线的长度。如果公垂线的长度比较短，我们就假设公垂线的中点为待求点(当公垂线的长度为0时，两异面直线为相交直线，交点即为待求点)；如果公垂线太长，那么可以断定所获取的左、右视图中的匹配点在匹配时误差太大，该点应该予以舍弃。

第九步、在世界坐标系中进行曲面插值并建立高度图

对图像在全黑背景之上显示的1400个离散像素点，必须在相邻离散点之间采用贝赛尔曲面插值算法进行插值以形成一个接近实际的连续高度图。高度图作为图像处理***的主要结果，包含了关于料场的全部信息以便用于后续的深层分析和数据库的数据存储。高度图是一个平面灰度图像，水平方向和世界坐标系的X轴方向一致，竖直方向和Z轴方向一致，某个点的像素灰度(0～255)和该点的高度(Y分量)成正比，将最低点0米对应于0(全黑)，最高点16米作为255(全白)，中间的高度对应于相应的灰色。

每一个像素代表0.5×0.5m²，这是一个1440×96pixel的图像，能够表达720m×48m的一个完整料场，其中的任何一点对应着相应的位置高度。

数字图像以BMP格式进行存放以供下一步的分析，以及用于数据库管理***的信息交互。

第十步、体积测量

①定义：高度低于0.5米的被定义为“空地”，即未堆料。首先从头至尾寻找空地，根据空地连续出现的位置判断应该将料场划分为几个料堆，同时排除很小的起伏造成的伪料堆和杂讯，得到每一个料堆的起始和终止的位置。

②从料场高度图中截取出料堆高度个图。

从理论上讲，这是一个求体积积分的运算，颇为复杂，但是利用高度图中的数据来计算料堆的体积，就显得十分简捷、直观。

令高度图中某一个像素的灰度值所代表的实际高度为h_i、底面积为δ_i，i为参与计算的像素点，i＝1，2，K，n，n为参与计算的像素点总数，因此，所有几何不规则的料堆体积V均可按照如下公式计算：

$V=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i{h}_i$ (公式17)

当δ_i＝const＝0.25m²时，公式17可以直接表达为

$V=0.25\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{h}_i$ (公式18)

定义0.5m以下为“空地”时，以包含料堆的矩形底面区域进行体积运算，公式还可以表达为：

$V=d_x{d}_y[\underset{l}{\overset{m}{\mathrm{\Σ}}}\underset{k}{\overset{n}{\mathrm{\Σ}}}{h}_{\mathrm{kl}}-0.5m\×n]$ (公式19)

式中，h_kl为第k列、第l行的像素点对应高度，k＝1，2，K，n，n为像素总列数，l＝1，2，K，m，m为像素总行数，dx、dy分别为像素点行向与列向的间隔所对应的实际距离，单位为m。采用公式19进行体积计算，在高度图中截取的矩形个图只要包含了要求计算的料堆，可以不必去顾及该料堆底面边界线的分布规律，因此能够简化诸多计算过程、节省CPU运算时间，对计算机编程极为方便、简捷。

Claims (5)

1.一种大型料场料堆视觉测量方法，其特征是，包括以下步骤：

步骤一、构建料场视觉测量***，由第一双目视觉子***、第二双目视觉子***和服务器三部分组成，服务器的图像信息处理模块接收第一双目视觉子***、第二双目视觉子***的信号后进行处理和运算后输出最终对图像信息的识别、理解和测量结果；

步骤二、安装双目视觉子***，将第一双目视觉子***安装于堆料机靠近料堆一侧机架上，将第二双目视觉子***安装于取料机靠近料堆一侧机架上，第一和第二双目视觉子***中的两个CCD摄像机的安装应该遵循：两个摄像机的光心处于同一水平面并且两个光轴与两光心连线形成一个等腰三角形；

步骤三、建立CCD摄像机几何模型，包括：世界坐标系与摄像机坐标系的变换、理想透视投影变换、归一化虚平面图像坐标系中图像畸变矫正数学模型、归一化虚平面图像坐标系畸变点坐标到计算机图像坐标系的转换，最后获得从世界坐标系到计算机图像坐标系的映射关系方程；

步骤四、CCD摄像机参数标定，包括内部参数和外部参数，采用分步标定方法，先标定图像的纵横比和主点坐标，再利用径向平行约束，求解大部分的外部参数，最后再引入畸变模型，线性求解内部参数和平移量；

步骤五、图像采集，两个双目视觉子***分别随着堆取料机的行走在线采集料堆及料场图像；每一个双目视觉子***在每一个采样时刻能够通过其中的两个CCD摄像机同时采集两幅图像，在前端信号处理器中保存图像序列，并将该两幅图像经过预处理后，同时接受由远程控制程序控制指令要求堆取料机上的前端信号处理器向上位服务器传送数字图像；服务器在每一个采样时刻能够同时接收到四幅图像，包括来自第一、第二双目视觉子***的各两幅图像，每两幅图像又分左、右两幅图像，在整个采样期间，所有采集到的图像分别形成第一、第二双目视觉子***的两组左右视图对序列；

步骤六、图像处理，消除图像畸变，由计算机图像坐标系像素点坐标(u，v)通过归一化虚平面图像坐标系畸变点坐标到计算机图像坐标系的转换方程计算求取对应点(u，v)在归一化虚平面图像坐标系畸变点坐标(x_d，y_d)，再将(x_d，y_d)的值输至代入归一化虚平面图像坐标系中图像畸变矫正数学模型进行求逆运算，因此求得矫正畸变后的归一化虚平面图像坐标系理想点坐标(x_u，y_u)，再由理想点坐标(x_u，y_u)替代(x_d，y_d)的值代入归一化虚平面图像坐标系中图像畸变矫正数学模型求出对应点(u，v)新的坐标值，即新的像素点位置；计算机图像坐标系中的所有像素点坐标(u，v)逐一经过上述计算过程后，就能够获得像素点在画面上的理想排列，即真实景物图像的恢复，或者说，获得一幅反映真实景物的理想图像；

步骤七、角点的寻找与匹配，在采样图像中寻找若干个角点作为代表最具特征的边缘点，采用基于图像灰度的角点检测算法，主要通过计算曲率及梯度来达到检测角点的目的；针对已经寻找出来的角点在外极线约束的指导下在左、右视图中寻找与之最为匹配的点；实际操作过程需要从左、右视图的数字图像中取得足够多的立体角点；

步骤八、物点三维坐标值的确定，在左、右视图中相互匹配的角点坐标(u_l，v_l)、(u_r，v_r)通过从世界坐标系到计算机图像坐标系的映射关系方程能够转换为对应的世界坐标(X_Wl，Y_Wl，Z_Wl)和(X_Wr，Y_Wr，Z_Wr)，采用求取两条直线公垂线中值的方法最后获得物点空间坐标的最优近似；

步骤九、在世界坐标系中进行曲面插值并建立高度图，对图像在全黑背景之上显示的有限个离散像素点，必须在相邻离散点之间采用贝赛尔曲面插值算法进行插值以形成一个接近实际的连续高度图；这个高度图作为图像处理***的主要结果，包含了关于料场的全部信息以便用于后续的深层分析和数据库的数据存储；

步骤十、体积测量，首先，定义“空地”，即未堆料的地面高度，从头至尾寻找空地，根据空地连续出现的位置判断应该将料场划分为几个料堆，并排除很小的起伏造成的伪料堆和杂讯，得到每一个料堆的起始和终止的位置；其次，从料场高度图中截取出料堆高度个图；令高度图中某一个像素的灰度值所代表的实际高度为h_i，该处料堆的体积可以看成由有限个高度为h_i、底面积为δ_i的微长方体累加而成(i为参与计算的像素点，i＝1，2，K，n，n为参与计算的像素点总数)，因此，所有几何不规则的料堆体积V均可按照如下公式计算：

2.根据权利要求1所述的大型料场料堆视觉测量方法，其特征是，所述的构建料场视觉测量***中每个双目视觉子***，包括：两个CCD摄像机、图像采集卡、前端信号处理器、两个光电转换器和光纤；左、右CCD摄像机在采样时刻将所采集到的左、右视图并行输入至图像采集卡；图像采集卡将原始图像转换成数字图像信号并输至前端信号处理器；前端信号处理器对数字图像进行去除噪声和强化预处理；经过预处理后的数字图像信号经第一个光电转换器将电信号转换成光信号，光信号通过光纤传输到第二个光电转换器；第二个光电转换器再将图像光信号转换成电信号，最后，电信号通过服务器的输入接口输入至服务器的图像信息处理模块；第一、第二双目视觉子***分别将各自采集到的左、右视图模拟信号并行输至本子***中的图像采集卡，图像采集卡将图像模拟信号转换成数字图像信号并输至前端信号处理器，前端信号处理器对数字图像进行去除噪声和强化预处理，经过预处理后的数字图像信号经第一个光电转换器将电信号转换成光信号，光信号通过光纤传输到第二个光电转换器，第二个光电转换器再将图像光信号转换成电信号，最后，电信号通过服务器的第一输入接口输入至服务器的图像信息处理模块，服务器将数字图像信号进行处理和运算后输出最终对图像信息的识别、理解和测量结果。

3.根据权利要求1所述的大型料场料堆视觉测量方法，其特征是，所述的CCD摄像机参数标定，具体如下：

CCD摄像机需要标定的外部参数：包括R和t，即

$\left[\begin{array}{cc}R& t\end{array}\right]=\left[\begin{array}{cccc}{r}_1& r_2& r_3& t_x\\ r_4& r_5& r_6& t_y\\ r_7& r_8& r_9& t_z\end{array}\right]$

共12个，但因R为单位正交矩阵，必须满足6个正交约束，故实际上只有6个外部参数需要标定；

内部参数：包括f、k₁、s_x和c_x、c_y，共5个；

首先确定图像的纵横比，垂直拍摄一个圆环，然后计算其在水平方向和垂直方向上的像素直径比即为纵横比s_x；

其次用相机的不同焦距拍摄同一景物，然后计算其缩放中心得到主点坐标(c_x，c_y)；

最后通过设定过渡参数求解R和t_x、t_y，然后根据R的正交性计算得到有效焦距f、径向畸变系数k₁和平移量t_z。

4.根据权利要求1～3中任意一权利要求所述的大型料场料堆视觉测量方法，其特征是，所述在世界坐标系中进行曲面插值并建立高度图，具体如下：

对图像在全黑背景之上显示足够多的离散像素点，如1400个，在相邻离散点之间采用贝赛尔曲面插值算法进行插值以形成一个接近实际的连续高度图，该高度图作为图像处理***的主要结果，包含了关于料场的全部信息以便用于后续的深层分析和数据库的数据存储，高度图是一个平面灰度图像，水平方向和世界坐标系的X轴方向一致，竖直方向和Z轴方向一致，某个点的像素灰度(0～255)和该点的高度(Y分量)成正比，将最低点0米对应于0(全黑)，最高点16米作为255(全白)，中间的高度对应于相应的灰色；

每一个像素代表一个面积，对应整幅图像就代表着一个特定长×宽的原料场，如：像素代表0.5×0.5m²，一个1440×96pixel的图像就代表720m×48m的一个完整料场，其中的任何一点对应着相应的位置高度；

数字图像以BMP格式进行存放以供下一步的分析，以及用于数据库管理***的信息交互。

5.根据权利要求1～3中任意一权利要求所述的大型料场料堆视觉测量方法，其特征是，所述体积测量，具体如下：

①定义：高度低于0.5米的被定义为“空地”，即未堆料；首先从头至尾寻找空地，根据空地连续出现的位置判断应该将料场划分为几个料堆，同时排除很小的起伏造成的伪料堆和杂讯，得到每一个料堆的起始和终止的位置；

②从料场高度图中截取出包含要求计算料堆的矩形个图，令高度图中某一个像素的灰度值所代表的实际高度为h_i、底面积为δ_i，i为参与计算的像素点，i＝1，2，K，n，n为参与计算的像素点总数，因此，所有几何不规则的料堆体积V均可按照公式

计算，当δ_i＝const＝0.25m²时，还可以直接表达为

当定义0.5m以下为“空地”时，以包含料堆的矩形底面区域进行体积运算，公式还可以表达为：

式中，h_kl为第k列、第l行的像素点对应高度，k＝1，2，K，n，n为像素总列数，l＝1，2，K，m，m为像素总行数，dx、dy分别为像素点行向与列向的间隔所对应的实际距离，单位为m。

Images