密集匹配点云下的矿石块度提取

崔欣， 王井利*， 吴冬， 赵鑫

(1. 沈阳建筑大学交通与测绘工程学院， 沈阳 110168； 2. 中铁十九局集团矿业投资有限公司， 北京 100161)

目前，对矿石块度测定方法的研究大多分直接测量法与间接测量法。直接测量方法包括筛分法、二次爆破统计法、爆堆直接测量法等。间接测量方法包括相关数据测量法、视觉分析法、摄影测量法[1]。文献[2-3]使用三维高分辨率激光扫描技术实现了爆破产生的矿石的全尺度量化。程雄等[4]采用数字图像处理技术实现对矿石块度参数的在线测量，并确定矿石块度参数的分布。张建立等[5]基于工业摄像机获取到的传送带上的矿石颗粒图像，利用矿石颗粒像素标定以及像素面积检测的方法研究了矿石的粒度分布。Ma等[6]通过对矿石破碎过程的分析，利用计算机视觉和图像处理技术，实现了矿石粒度的在线检测和粒度分布的计算。同样的，文献[7-9]也采用了图像处理、图像识别等图像分析技术对矿石的块度进行测定与统计。此类方法对获取的矿石影像要求较高，矿石块度提取精度受图像分辨率的限制。Maron[10]提出了一种通过安装新的粒度跟踪(particle size tracking，PST)测量技术来估算矿石粒径的方法，实现了在线粒度测量，由于数据量较大，算法处理实时性较差。吴开兴等[11]使用水洗筛分法研究了单颗粒粒度分布特征，该算法不适用于大量矿石参数测定。Krauze等[12]基于震动信号处理方法开发一种基于振动的松散固体粒度和流动在线测量方法，实现了颗粒粒径和流速的估算。在利用图像识别块度领域，Noy等[13]引入矿石堆的立体相机成像为基于图像的测量算法提供了改进的自动粒子分割；孙厚广等[14]对爆破后的矿山形态进行表面全景垂直拍照，通过软件将照片划分规则格网，通过对规则格网的分析，来评价爆破的块度问题。此类基于图像处理技术的块度测定方法对相机分辨率要求较高。喻战江等[15]利用三维激光扫描技术获取爆破碎块的几何尺寸，实现了块度分级，但并未真正意义上达到块度识别的目的，方法精度有待提高。谢博等[16]利用3D激光扫描仪获得岩块表面模型，在内业处理中，通过不同的超体素聚类的阈值设置和权重因子分配，得到最佳的爆堆点云超体素聚类效果，之后利用区域生长法的思想进行分割达到块度识别的目的。上述方法多采用图像识别或者三维激光扫描技术进行矿石块度测定，其主要问题在于矿石块度精度低、提取效率差、仅适用于少量矿石块度测定。为此，现采用无人机摄影测量密集匹配点云进行矿石块度估算，旨在提高矿石块度估算的精度与效率。

1 数据预处理

采用大疆Phantom 4 RTK无人机获取了新疆磁海铁矿某矿石堆的影像数据，利用ContextCapture Master对影像进行空中三角测量以及密集匹配点云模型的构建，构建出的三维模型如图1所示。在完成三维点云模型构建之后针对矿石点云与地面点云相互粘黏问题，采用布料模拟滤波算法对点云进行滤波处理，实现地面点与非地面点的分离。经实验发现在该算法中当布料分辨率为0.5，分类阈值为0.2，最大迭代次数取500时滤波效果最为理想。为了便于后期能够快速地提取出每个矿石堆的矿石块度，在分离出的矿石堆点云中采用基于欧氏距离的点云分割算法实现单个矿石堆点云的有效分割，最终经过预处理后的矿石堆分割结果如图2所示。

图1 研究区三维点云模型Fig.1 Three-dimensional point cloud model of the study area

图2 矿石堆分割结果Fig.2 Results of ore pile segmentation

2 矿石堆坡度估算

坡度(slope)是地表单元陡缓的程度，通常把坡面的垂直高度h和水平方向的距离d的比叫作坡度，θ为该坡度下对应的坡度角，如图3所示。其数学表达式为

(1)

图3 地形坡度Fig.3 Topographic slope

依据坡度的定义式可知，寻找出矿石堆某一坡面处矿石点的最高点以及最低点即可完成对该坡面的坡度估计。对于同一矿石堆而言，最高点是唯一的，即矿石堆中高程最大值所对应的点。由于三维点云中点的三维坐标是通过(x，y，z)进行表示的，在点云处理过程中默认将三维坐标中的z分量作为高程值。所以本文研究在进行最高点确定时，首先遍历所有坐标点的z值；其次比较了所有坐标点对应的z值，将z值最大的坐标点(xi，yi，zmax)作为矿石堆的最高点P，其中xi、yi表示点云中z值最大点对应的x与y轴分量。

矿石堆不同坡面的最低点并不是唯一确定的，最低点的选取与坡面有关。在本节中最低点是指在计算矿石堆某一坡向处矿石堆坡度时，该坡面的坡脚所对应的高程最低点。

为确定矿石堆不同坡面所对应的最低点，笔者以矿石堆点云的最高点为起始点，基于k近邻搜索以及坐标方位角进行最低点搜寻，其整个实现步骤如下。

步骤1利用k近邻搜索法搜寻距离最高点P的距离在[r，r+Δr]的所有离散点Ω。

步骤2将获取到的离散点Ω与最高点P均投影至XOY平面上得到点P′和Ω′，在平面中计算环形点云Ω′中每一点与P′点连线的坐标方位角。

步骤3选取步骤2中计算的方位角与0°0′0″最为接近的点，以该点的投影前的坐标点作为下一搜索的中心点P1，重复步骤1～步骤3继续进行最低点搜索，当同一方位角搜索出来的前后两坐标点Pn、Pn+1的高程值(即z值)之差小于设定的阈值ε，即|Pn+1z-Pnz|<ε时，停止搜索，点Pn即为矿石堆在0°0′0″方向处的最低点。

步骤4同理选取步骤2中计算的方位角与kΔα(k=1，2，…且kΔα≤360°，Δα取值越小其后期分割结果越理想，本文中Δα取6°)最为接近的点，类似于寻找在0°0′0″方向处的最低点寻找出各个方位处矿石堆的最低点。

整个最低点的搜索过程剖面示意图与俯视图分别如图4和图5所示。

图4 最低点搜寻剖面示意图Fig.4 Schematic diagram of the lowest point search profile

图5 最低点搜寻俯视图Fig.5 The top view of the lowest point search

在完成矿石堆点云最高点以及最低点提取后，根据两最值点三维坐标计算出矿石堆坡度i。计算公式为

(2)

式(2)中：xi、yi、zmax表示点云中最高点对应的x、y、z轴分量；xj、yj、zmin表示搜寻到的最低点对应的x、y、z轴分量。研究中将坡度计算结果与人工测量的坡度结果进行对比，分析坡度计算的误差，其误差分布直方图如图6所示。从图6中可以看出计算得到大多数坡度的绝对误差在0.1以内，且测定误差呈现正态分布。

图6 矿石堆坡度计算误差直方图Fig.6 Histogram of calculation error of ore pile slope

3 基于坡度信息的矿石点云分割

通过对矿石堆坡面进行逐层切割寻找出相邻矿石间的交界点，以此获得矿石点云边界线，最后依据得到的边界线完成对矿石点云的分割。

3.1 基于坡度信息的矿石堆点云切割

矿石堆点云的切割平面主要由矿石堆的坡度i以及最高点与最低点的三维坐标所决定。其基本原理如图7所示。

利用求得矿石堆点云最高点P、最低点Q以及坡度i确定一个过最高点P与最低点Q同时与地面夹角为坡度角α的平面β0。将平面β0作为初始矿石堆点云裁剪的平面。其中平面β0的表达式为

(3)

式(3)中：xP、yP、zP分别为最高点P对应的三维坐标x、y、z分量；xQ、yQ、zQ分别为最高点Q对应的三维坐标x、y、z分量。

图7 裁剪平面确定示意图Fig.7 Schematic diagram of determining the clipping plane

但是由于矿石堆的坡面上会存在着多个凸起或凹陷的矿石，利用单个裁剪平面β0进行裁剪无法保证坡面上所有的矿石点云都被裁剪下来。另一方面，单一平面无法保证矿石点云与该平面相交所截出的平面为该矿石最大的截面，从而影响矿石块度计算的精度。针对上述问题，提出了多层切片的裁剪方式。通过沿着平面β0的法向量方向上下平移平面β0得到更多的裁剪平面。

平面β0上下移动的距离以及最终平移的最大距离并非无限制的。而平面β0最终平移的最大距离由坡面表面矿石隆起的最大高度以及矿石凹陷的最大深度所决定的。本文以该坡面上所有离散点到平面β0的距离作为平面β0最终平移的最大距离判定依据。其实现步骤如下。

步骤1计算坡面内所有离散点(xk，yk，zk)，k=1，2，…到平面β0的距离dis：

(4)

式(4)中：xP、yP、zP为最高点P对应的三维坐标x、y、z分量；xQ、yQ、zQ为最高点Q对应的三维坐标x、y、z分量；xk、yk、zk(k=1，2，…)为第k个离散点对应的三维坐标x、y、z分量。

步骤2判断离散点与平面的位置关系，并对距离dis进行方向标定，如式(5)所示。其中“+”表示离散点位于平面之上，“-”表示离散点位于平面之下。

(5)

式(5)中：zk(k=1，2，…)为第k个离散点对应的三维坐标z分量；z′k(k=1，2，…)为第k个离散点的x、y分量对应于平面该点处的高程值。

步骤3分别求取平面上下采样点到平面的最大距离max_dis+、max_dis-，将两个最大距离作为平面β0沿着法向量方向上下平移的极限值max_l+、max_l-。

在已知平面β0沿着法向量方向上下平移的极限值max_l+、max_l-后，可计算得到最高点及最低点上下移动的极值max_h+、max_h-分别为

(6)

根据最高点及最低点上下移动的极值max_h+、max_h-可知，最高点高程值的变化范围为(zP- max_h-，zP+max_h+)。

由此可以得出第k个裁剪平面经过点(xP，yP，zP- max_h-+kΔh)。故第k个裁剪平面的平面方程为

[z-(zP-max_h-+kΔh)]=0

(7)

式(7)中：xP、yP、zP为最高点P对应的三维坐标x、y、z分量；xQ、yQ、zQ为最高点Q对应的三维坐标x、y、z分量；k为第k个裁剪平面；Δh为初始平面每次平移量；max_h+为平面上移的最大距离。

点云裁剪平面确定完成后即可根据点云中离散点的空间坐标以及裁剪平面的空间位置进行点云裁剪。

3.2 矿石点云交界线粗提取

以估算得到的坡度i为基础对矿石堆点云进行多层切割处理。对前后两次不同裁剪平面裁剪得到的点云进行分割处理，同时计算分割得到的每个矿石点云的中心点，其结果如图8所示。依据前后两次计算得到的中心点的数量以及位置关系进行矿石边界线粗提取。

在图8中，若两次裁剪后的矿石点云中中心点数量以及位置关系相等或相近，则判断前后两次分割点云均来自同一个矿石，如图8(a)～图8(c)所示3个不同裁剪平面对同一坡面点云进行裁剪，其裁剪后的矿石点云中心点均为4个且4个中心点的位置一一对应[图8(e)～图8(h)]，则认为这3次裁剪出来的点云中不含有边界点云；若两个裁剪出的矿石点云中心点数量不一致，同时中心点位置也无法一一对应，则认为两次切割点云中含有交界点，如图8(c)～图8(d)所示，两个裁剪平面裁剪后的矿石点云分别形成了4个和3个点云，且点云中心点存在不对应现象，如图8(g)和图8(h)，表明图8(d)裁剪出的点云中包含两个矿石点云交界点。

图8 不同裁剪平面切割出的点云Fig.8 Point clouds cut by different cutting planes

为了确定出矿石的交界点，将8(g)～图8(h)结果中中心点不对应的点云进行作差处理，其结果如图9所示，图9中红色点云即为作差得到的结果Ω。

两点云作差结果中不仅包括了两矿石的交界点，还包括矿石非交界点。为了剔除非交界点，计算图8(h)结果中的两中心点P1、P2到点云Ω中每一点的距离，并按从小到大进行排序，分别选取距离P1、P2点距离在某一阈值(实验中选取为前60%)的点Ω1、Ω2，如图10所示，图10中红色点云表示距离P1点距离在前60%的点Ω1，蓝色点云表示距离P2点距离在前60%的点Ω2。最后对Ω1、Ω2进行求交运算，求交结果即为粗提取得到的两矿石交界点，图10中橘黄色点即为交界点。

图9 点云作差结果Fig.9 Point cloud difference result

图10 粗提取得到的交界点Fig.10 Junction points obtained by rough extraction

3.3 矿石点云交界线精提取

在3.2节中，受阈值选取的影响，提取的交界点中仍存在着非交界点，为提高矿石的分割精度，对交界点进行了进一步的精提取。

3.3.1 灰度信息提取

受光照条件的影响，矿石点云中矿石交界线处点云照比其他点云较黑，故在研究中将点云的红(R)、绿(G)、蓝(B)三个颜色通道信息转换成灰度信息，其灰度值Gray的计算公式为

Gray=0.299R+0.587G+0.114B

(8)

式(8)中：Gray为点云灰度值；R、G、B为点云的红、绿、蓝3个颜色通道分量。

将计算得到的点云灰度值进行可视化，图11(a)为矿石原始点云，图11(b)为点云灰度可视化结果。从图11(b)可以看出两矿石交界处的点云呈现蓝色或暗蓝色，故点云的灰度信息可以用于矿石分割的依据。

3.3.2 邻域特征分析

通过对两矿石交界处点云进行分析发现，交界点处的点云高程值小于其周边的大多数点云高程。

图11 点云灰度可视化前后对比Fig.11 Comparison before and after point cloud gray scale visualization

依据这一特征可以对提取出的矿石交界处点云进行验证，具体方法如下。

首先分别以每个粗提取出的矿石点云边界点为中心进行K邻近搜索寻找邻近点云；其次比较近邻点云高程与中心点高程的大小，记录高程大于中心的邻近点数量；若大于中心的邻近点数量占邻近点总数量的比值超过某一阈值，则判定该点为矿石交界点，反之则为非交界点。

在完成矿石交界点提取后，根据矿石交界点的连线完成对单个矿石的分割。图12(a)为最终提取出的矿石点云的交界点，图12(b)为依据提取出的矿石交界点分割得到的单个矿石结果。从图12(b)中可以看出分割结果良好，但也存在着少量欠分割和过分割现象。

4 块度测定及分布规律研究

矿石块度(lumpiness)是指采用爆破等手段所开采的矿石碎块的几何尺寸的大小，常用矿石碎块两端最大距离表示。在图像处理领域，学者们多采用等效匹配法对矿石参数进行估算，该方法中最为常用的算法包括最小外接矩形以及最小拟合椭圆等。

将分割得到的矿石点云进行投影变换得到近似于矿石最大截面的点云，采用了基于形状特征的矿石块度参数测定方法。即根据矿石截面的形状特征选用不同的图形对矿石截面边界点云进行等效匹配，使得计算出的截面周长与截面面积误差最小。本文采用了圆度的概念对矿石截面的形状进行描述。圆度(circ_rity)是指矿石最大截面接近理论圆的程度。当圆度circ_rity越接近1时，表示截面形状越近似于圆，此时采用最小拟合椭圆对矿石的最大截面进行等效匹配；否则采用最小外接矩形进行等效匹配。

图12 基于矿石交界线的点云分割Fig.12 Point cloud segmentation based on ore boundary

圆度的定义式为

(9)

式(9)中：S为矿石截面面积；C为矿石截面周长。

实验中选取了一个矿石堆中的32个矿石进行分析，计算出了32个矿石利用最小拟合椭圆和最小外接矩形进行等效匹配时面积、周长误差与矿石最大截面圆度的关系。其结果如图13(a)和图13(b)所示。

从图13(a)可知，当圆度circ_rity≥0.55时，利用椭圆拟合计算得到的周长误差普遍小于利用最小外接矩形计算得到的结果，且大部分误差在10%以内。从图13(b)中可以看出，当circ_rity≥0.55时，利用椭圆拟合计算得到的面积误差普遍小于利用最小外接矩形计算得到的结果。综上所述，在计算矿石参数前，首先进行圆度计算，当圆度circ_rity≥0.55时，利用椭圆拟合的方式进行等效匹配；否则

图13 周长、面积误差随圆度变化Fig.13 Circumference and area error change with roundness

利用最小外接矩形的方式进行等效匹配。并将拟合椭圆的长半轴或者最小外接矩形的长度作为矿石块度l：

(10)

式(10)中：C为矿石最大截面的周长；S为矿石最大截面的面积；circ_rity为矿石最大截面圆度。

将提取出的矿石最大截面的长度作为最终提取出的矿石块度，并利用统计学方式对矿石块度进行统计分析，其结果如表1所示。同时绘制了矿石长度分布直方图，如图14所示。

从表1以及图14中可以看出，矿石块度介于0.5～1.0 m的比例最大，约为41.63%。块度在4 m以上的矿石数量较少，最大矿石的块度为7.674 m。

表1 矿石块度统计表Table 1 Statistical table of ore lumpiness

图14 矿石块度分布直方图Fig.14 Histogram of ore lumpiness distribution

5 结论

针对矿石块度测定困难、精度低等问题，提出了基于密集匹配点云的矿石块度计算方法。通过对本文方法进行试验验证，得到以下结论。

(1)利用搜寻矿石堆最高点、最低点的方式计算出的矿石堆坡度具有一定的可靠性，其计算误差在0.1以内，且满足正态分布。

(2)利用矿石堆的坡度信息、点云灰度信息以及点云邻域特征可以较好地实现单个矿石点云的有效分割，且分割结果满足后期矿石块度计算的要求。

(3)研究区内的矿石块度较为均匀，矿石块度大多介于0.5～1.0 m。