王植, 安世缘, 邹俊, 张紫瑞
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2020-10-12
基金项目：中央高校基本科研业务费专项资金资助项目(N170113027)。
作者简介：王植(1979-)，男，辽宁沈阳人，东北大学副教授。

摘要：台阶线信息对于露天开采具有重要价值, 现有获取台阶线的方法工作量大、效率低、精度差, 降低了矿山的生产效率和验收精度.因此, 本文基于序列无人机影像生成的露天矿密集点云数据, 研究并提出了一种自动提取露天矿台阶线的方法.该方法利用渐进形态学滤波算法对点云进行预处理, 提出一种顾及邻域几何属性的三维边缘检测与曲率指数加权方法提取出台阶线特征点, 并使用移动最小二乘法精确拟合出台阶线.实验结果表明该算法可以自动、高效、精确地提取出露天矿台阶线, 生成露天开采现状图, 对于露天矿生产和安全具有重要的应用价值.
关键词：露天矿台阶线三维点云曲率自动提取
Step Line Extraction from Point Cloud Data of Open-Pit Mine
WANG Zhi, AN Shi-yuan, ZOU Jun, ZHANG Zi-rui
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: WANG Zhi, E-mail: wangzhi@mail.neu.edu.cn.

Abstract: The information of step line is of great importance to open-pit mining. The existing method of obtaining step line has large workload, low efficiency and poor accuracy, which reduces the production efficiency and acceptance accuracy of the mine.Thus, a method of automatically extracting open-pit mine step lines from the dense point cloud data of open-pit mines generated by sequence UAV images is proposed in this paper. This method uses progressive morphological filtering algorithm to preprocess the point cloud, and a three-dimensional edge detection and curvature index weighting method that takes into account the geometric properties of the neighborhood is proposed to extract the feature points of the step line, then, it uses the moving least squares method to accurately fit the step line. The experimental results show that the algorithm can automatically, efficiently and accurately extract the step line of the open-pit mine, and generate open-pit mining status map. It has important application value for open-pit mine production and safety.
Key words: open-pit minesstep line3D point cloudcurvatureautomatic extraction
露天矿在开采过程中, 为了便于从上至下逐层开采, 要把矿区内的矿物、岩石划分为具有一定厚度的水平分层, 这些分层的工作面就是台阶.露天矿台阶线对露天矿的采矿设计、计划编制、区块设计、爆破设计、矿区划分等业务均具有重要的作用.目前露天矿的工作实践中, 每个月需要大量的人力去计算采剥量，通过GPS获取有限的坐标点, 生成DSM (digital surface model, DSM) 进行人工绘制, 生成台阶线构建采场现状平面图, 这种方法工作量大、效率低、精度差, 降低了矿山的生产效率和验收精度^[1].本文的研究利用轻便的无人机, 获取露天矿序列影像, 进行三维建模, 自动提取采场台阶线, 生成采场现状三维模型, 对露天矿的智能开采和成本管理具有应用价值, 对加强矿山生产监控和安全管理具有重要意义^[2].
地物的几何属性可以由三维边缘参数表述出来, 现有的特征线提取方法主要有以下几种：
1) 通过计算采样点与其邻域点的法向量夹角或曲率变化来判定是否为边缘点, 如Demarsin等^[3]用主成分分析法求法向量, 然后采用区域生长算法对法向量变化较大的点进行分割, 并聚类分割数据, 从而得到特征线；Daniels等^[4]采用稳健移动最小二乘法(robust moving least squares, RMLS)对每个点的邻域进行曲面拟合, 然后计算点的投影残差值, 将该值较大的点判别为潜在的特征点；马聪聪等^[5]将点云数据中每个点的法向量与其邻域内各点的法向量点积的平均值作为特征度和点云密度结合组成判别参数提取特征点.
2) 根据采样点与其邻域点的特征相似性进行区域分割, 然后通过轮廓提取算法提取特征^[6], 如Sampath等^[7]将建筑与非建筑点云进行区域分割, 跟踪边界点形成规则的边界轮廓, 然后使用改进的凸包算法得到点云的凹轮廓; Xu等^[8]基于顶点法向量的聚类算法分割点云, 通过表面粗糙度区分原始面和断裂面, 最终基于表面分割提取边缘特征线.Bazazian等^[9]根据测地距离分割点云, 然后根据采样点的局部邻域属性定义多尺度算子, 可以利用该算子确定特征的连续性.
3) 计算采样点处拟合曲面的曲率进行区域分割, Kim等^[10]使用移动最小二乘法拟合曲面并估算曲率及其导数, 利用Delaunay三角剖分来计算邻域信息, 提取特征点, 最终根据曲率确定的线宽绘制山脊线和谷线.史皓良等^[11]使用主成分分析法和局部二次曲面拟合法估算平均法向量夹角和平均曲率, 根据双阈值检测方法进行特征提取.
综上所述, 本文针对现有方法中存在的边缘和特征线检测不全, 且当场景复杂时, 在非平面区域会产生现实中不存在的边缘或特征线等问题, 首先对露天矿点云数据使用渐进形态学滤波算法进行预处理, 提出了顾及邻域几何属性的三维边缘检测与曲率指数加权方法来提取台阶线特征点, 并采用移动最小二乘法高效、精准地拟合台阶线.
1 台阶线提取算法本文基于辽宁省鞍山市某铁矿的点云数据提取台阶线特征点, 然后使用移动最小二乘法拟合台阶线.算法流程图如图 1所示.
图 1(Fig. 1)

图 1 算法流程图Fig.1 Algorithm flow chart

1.1 渐进形态学滤波渐进形态学滤波是通过逐渐改变滤波器的窗口大小和高程差阈值, 移除点云数据中的车辆、植被、建筑等地物, 提取地面点.该方法可以根据研究区域的地形坡度确定高程差阈值, 假设坡度恒定, 则地形的最大高程差、窗口大小和地形坡度之间存在下列公式：

(1)

(2)

式中：dh₀表示初始高程差阈值; s表示斜率; c表示网格大小; dh_max表示最大高程差阈值; w_k表示运行到第k次窗口大小.
渐进形态学滤波在进行第一次迭代时, 将最小高程表面与初始过滤窗口大小作为输入.在后面的迭代中, 将从上一迭代中获得的已滤波表面和增加的窗口大小作为滤波器的输入, 根据不断增加的窗口对三维点云数据做开运算, 如果开运算前后的高程差小于阈值, 则程序将窗口大小赋值给标记数组, 然后继续进行开运算迭代直到窗口大小为最大窗口^[12].
1.2 曲率提取潜在特征点露天矿台阶特征线对采矿设计优化和采剥进度计划编制都具有重要作用, 台阶线的变化反映着露天矿台阶和边坡的变化, 且台阶线上的特征点均是具有曲率极值的点, 因此本文采用曲率加权的方法初步筛选潜在特征点.
对点云中某点P及其邻域内的点使用最小二乘法拟合曲面并对其进行参数化^[13]: r(x, y)=(x, y, ax²+bxy+cy²), 设r(x, y)的偏导数分别为r_x, r_y, r_xx, r_xy, r_yy, 然后根据式(3), 式(4)计算平均曲率H和高斯曲率K:

(3)

(4)

式中: L=r_x·n=0, E=r_x·r_x=1+(2ax+by)²; F=r_x·r_y=(2ax+by)(2cy+bx); G=r_y·r_y=1+(2cy+bx)²; ; .
最后, 通过下列公式计算主曲率:

(5)

(6)

曲率指数与坐标系无关, 具有旋转和位移不变性, 在三维特征点检测领域中被广泛使用.当拟合曲面弯曲程度较大时, 曲率指数较大; 反之, 曲率指数较小.曲率指数C_(p)由式(7)求得^[14]:

(7)

1.3 三维边缘检测算法传统的从点云中提取三维边缘方法只能提取出两个或多个平面相交形成的边缘, 且不能提取由小尺寸、狭长平面构成的三维边缘^[15].本文将倪欢等提出的AGPN (analysis of geometric properties of neighborhoods)算法^[16]应用在露天矿点云中, 细化特征点.使用该算法提取边缘时点密度对其检测三维边缘没有影响, 且可以从含有噪声的点云中检测出所有类型的三维边缘.
AGPN算法引入了RANSAC(random sample consensus)算法.露天矿地形复杂, 少有规则的曲面, 露天矿台阶线即为不规则曲面相交构成的边缘, 由于曲面的微小局部可以近似看成是平面, 平面模型可以有效地近似局部曲面的几何结构, 因此, 本文采用平面模型的RANSAC算法.如果潜在特征点中的一点P是折叠边缘, 则使用k-d树搜索邻域后, 空间向量(i=1，…，n)在点P邻域的拟合平面上的投影会有较大的夹角(如图 2a所示).该折叠边缘同时位于拟合平面和邻域内另一个相交表面上.当点P不是折叠边缘但也属于内点时, 那么该点位于内点的拟合平面的内部, 此时, 空间向量在拟合平面上的投影之间不存在较大的夹角, 是连续分布的(如图 2b所示).此外, 若内点集的点数过少, 则认为点P为噪声点.
图 2(Fig. 2)

图 2 两个相交表面局部邻域结构几何形态Fig.2 Geometry of local neighborhood structure of two intersecting surfaces (a)-点P位于折叠边缘；(b)-点P位于局部拟合平面内部.

1.3.1 法向量优化现阶段常用的对点云中某点P的法向量估计方法是对点P使用k-d树搜索邻域, 对其邻域采用最小二乘法拟合平面, 构造协方差矩阵, 求解协方差矩阵的特征值和特征向量, 然后选择最小特征值对应的特征向量, 并进行单位化, 最终将该向量作为点P的法向量.目前计算角度阶跃特征的方法, 大多数都采用该法向量构建空间坐标系, 并在垂直于该法向量的平面上计算角度阶跃值.这样将会导致此方法计算出来的角度阶跃值不能用于检测折叠边缘, 并且对于曲面轮廓, 也不能准确检测.
AGPN算法则是对点P的邻域使用平面模型的RANSAC算法, 将其分为内点和外点, 内点拟合平面的法向量作为点P优化后的法向量.当点P的邻域有多个平面时, 优化后的法向量将会垂直于其中一个平面.
1.3.2 角度阶跃计算根据优化后的法向量和拟合平面上一对相互垂直的向量建立空间直角坐标系, 设空间向量u在x轴上, 向量v在y轴上, 优化后的法向量为n, 则角度阶跃计算公式如下:

(8)

(9)

式中空间向量是由当前点P和其邻域内拟合平面提取的内点构成的.计算角度阶跃G_θ后, 将角度阶跃大于阈值的点作为提取的特征点.
1.4 移动最小二乘法传统的最小二乘法拟合曲线是在计算偏差后, 让平方和达到最小值, 然后求解线性方程组, 最终得到拟合曲线.该方法具有全局性, 但是不能满足定位的需要, 不适用于具有大量数据集的模型.与传统的最小二乘法相比, 移动最小二乘法在最终拟合曲线的定位和光滑性方面有较好的效果, 且拟合精度也比较高.移动最小二乘法有两个明显的优点: 一是不同于传统方法中的多项式, 拟合过程中使用的函数由基函数和系数向量组成; 二是引入了紧支的概念, 即点x处的值y只受x支撑域内节点的影响^[17].在要进行曲线拟合的部分中, 其目标子区域在拟合过程中使用的函数记作f(x):

(10)

其中: α(x)代表的为待定的系数矢量; p^T(x)代表的为m次多项式空间的一组基.

(11)

(12)

若计算点x的邻域内有n个节点, 近似函数f(x)在这些节点处的误差的加权平方和如下:

(13)

其中：y_i是x=x_i处的节点值; w_i(x)是具有紧支特性的权函数.
如果想要解出待定参数α(x), 就应让式(13)中的J取得最小值.对J取极值即对α(x)求导, 则

(14)

其中:

(15)

(16)

(17)

计算时需要用到形函数φ^k(x):

(18)

获得拟合函数:

(19)

2 实验分析本文使用的实验数据为辽宁省鞍山市某铁矿矿区的无人机数据.本文将无人机采集的铁矿矿区序列影像使用Smart 3D软件通过添加控制点, 空三解算等步骤生成密集点云数据, 再根据点云数据使用CloudCompare软件建立三维模型.图 3为建立的矿区三维模型.本文使用C++, Python和点云库PCL实现算法.
图 3(Fig. 3)

图 3 实验数据Fig.3 Experimental data

露天矿中挡土墙等人工地物会干扰台阶线特征点的提取, 本文使用渐进形态学滤波将点云分为地面点和非地面点, 然后将非地面点(挡土墙区域的点云)自适应地减去一定的高程, 最后将地面点和非地面点拼接, 提取特征点.滤波前后特征点提取结果对比如图 4所示.
图 4(Fig. 4)

图 4 渐进形态学滤波结果Fig.4 The results of progressive morphological filtering (a)-滤波前特征点提取结果；(b)-滤波后特征点提取结果.

为了保证露天矿上岩石或岩体的稳定性, 台阶设计的坡面角普遍为50°~75°, 但是在实际工程中, 坡面角通常会小于设计值.对露天矿点云建立DEM, 其三角网倾角的整体分布是多峰分布, 台阶面、道路和平台面坡度分界的临界值为30°, 该临界值可作为阈值过滤三角网, 这是露天矿坑、排土场所特有的.本文使用AGPN算法对点云数据进行法向量优化, 结果如图 5所示, 根据露天矿独有的特点, 使用不同的颜色表示具有不同方向法向量的点.其中法向量与水平面夹角为0°~50°使用紫色点表示, 与水平面夹角为50°~75°使用绿色点表示, 与水平面夹角为75°~90°使用红色点表示.根据法向量优化结果可以看出红色点覆盖了大部分的台阶面, 紫色点、绿色点和红色点有较明显的边界，台阶面、边坡面已经较好分离开.
图 5(Fig. 5)

图 5 法向量优化结果Fig.5 The results of normal vector optimization

在法向量优化后, 根据AGPN算法进行角度阶跃的计算, 本文采取的角度阶跃阈值为120°, 即将根据式(9)计算的G_θ>120°的点作为特征点. AGPN算法中RANSAC的距离阈值对最后的结果会有较大影响, 根据露天矿地形特征和实验总结, 本文采用的平面模型RANSAC算法的距离阈值为0.25 m, 最终提取结果如图 6所示.
图 6(Fig. 6)

图 6 特征点提取结果Fig.6 The results of feature point extraction

图 6中红色点为提取出来的特征点, 根据建立的露天矿三维模型和特征点可以得出: 本文使用的方法可以较好地提取露天矿台阶线特征点.最终采用区域增长法分割台阶线特征点, 对特征点采用移动最小二乘法拟合曲线, 拟合结果如图 7所示.
图 7(Fig. 7)

图 7 台阶线拟合结果Fig.7 The results of step line fitting

3 结语本文基于无人机影像生成的露天矿密集点云数据, 对露天矿台阶线自动提取算法进行了研究, 使用渐进形态学滤波处理点云数据, 根据曲率筛选潜在特征点, 提出对潜在特征点采用曲率指数加权的AGPN算法提取台阶线特征点, 最终采用移动最小二乘法拟合出台阶线.实验结果表明: 1) 渐进形态学滤波可以有效地解决露天矿中挡土墙等人工地物对实验结果的干扰; 2) 曲率指数结合AGPN算法可以较精确地提取出台阶线上的特征点.本文的研究成果对于提高露天矿测绘自动化程度, 加强露天矿生产与安全管理等方面, 具有重要的应用价值.
参考文献

[1]	Paraskevis N, Roumpos C, Stathopoulos N, et al. Spatial analysis and evaluation of a coal deposit by coupling AHP & GIS techniques[J]. International Journal of Mining Science and Technology, 2019, 29(6): 943-953. DOI:10.1016/j.ijmst.2019.04.002
[2]	Ulf K, Torsten G, Reimar S, et al. The Muruntau gold deposit(Uzbekistan)—a unique ancient hydrothermal system in the southern Tien Shan[J]. Geoscience Frontiers, 2016, 7(3): 495-528. DOI:10.1016/j.gsf.2015.09.005
[3]	Demarsin K, Vanderstraeten D, Volodine T, et al. Detection of closed sharp edges in point clouds using normal estimation and graph theory[J]. Computer-Aided Design, 2006, 39(4): 276-283.
[4]	Daniels J I, Ha L K, Ochotta T, et al. Robust smooth feature extraction from point clouds[C] //Proceedings of IEEE International Conference on Shape Modeling and Applications. Minneapolis, 2007: 123-136.
[5]	马聪聪, 李松, 曹菁菁, 等. 基于法向量和密度的点云特征点提取问题研究[J]. 计算机应用与软件, 2020, 37(5): 256-260,292. (Ma Cong-cong, Li Song, Cao Jing-jing, et al. Feature points extraction of point cloud based on normal vector and density[J]. Computer Applications and Software, 2020, 37(5): 256-260. DOI:10.3969/j.issn.1000-386x.2020.05.044)
[6]	Rouhani M, Lafarge F, Alliez P. Semantic segmentation of 3D textured meshes for urban scene analysis[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2017, 123: 124-139. DOI:10.1016/j.isprsjprs.2016.12.001
[7]	Sampath A, Shan J. Building boundary tracing and regularization from airborne lidar point clouds[J]. Photogrammetric Engineering & Remote Sensing, 2007, 73(7): 805-812.
[8]	Xu J Y, Zhou M Q, Wu Z K, et al. Robust surface segmentation and edge feature lines extraction from fractured fragments of relics[J]. Journal of Computational Design and Engineering, 2015, 2(2): 79-87. DOI:10.1016/j.jcde.2014.12.002
[9]	Bazazian D, Casas J R, Ruiz-Hidalgo J. Segmentation-based multi-scale edge extraction to measure the persistence of features in unorganized point clouds [C]// International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications. Porto, 2017: 317-325.
[10]	Soo-Kyun K. Extraction of ridge and valley lines from unorganized points[J]. Multimedia Tools and Applications, 2013, 63(1): 265-279. DOI:10.1007/s11042-012-0999-y
[11]	史皓良, 吴禄慎, 余喆琦, 等. 散乱点云数据特征信息提取算法[J]. 计算机工程, 2017, 43(8): 279-283. (Shi Hao-liang, Wu Lu-shen, Yu Zhe-qi, et al. Algorithm for feature information extraction from scattered point cloud data[J]. Computer Engineering, 2017, 43(8): 279-283. DOI:10.3969/j.issn.1000-3428.2017.08.047)
[12]	Zhang K, Chen S C, Whitman D, et al. A progressive morphological filter for removing nonground measurement from lidar data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(4): 872-882. DOI:10.1109/TGRS.2003.810682
[13]	Hawraa A, Yulia H, David M, et al. A 3D morphometric perspective for facial gender analysis and classification using geodesic path curvature features[J]. Computational Visual Media, 2018, 4(1): 17-32. DOI:10.1007/s41095-017-0097-1
[14]	李燕春. 基于特征点的三维人脸识别方法研究[D]. 南京: 东南大学, 2015. (Li Yan-chun. Research of 3D face recognition based on facial landmarks [D]. Nanjing: Southeast University, 2015. )
[15]	Lin Y B, Wang C, Cheng J, et al. Line segment extraction for large scale unorganized point clouds[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 102: 172-183. DOI:10.1016/j.isprsjprs.2014.12.027
[16]	倪欢. 机载激光雷达点云多实体多层次分类方法[D]. 武汉: 武汉大学, 2017. (Ni Huan. Multiple-entity based heriarchical classification method for airborne lidar point clouds[D]. Wuhan: Wuhan University, 2017. )
[17]	贺彤. 点云的特征线提取算法研究及应用[D]. 太原: 中北大学, 2018. (He Tong. Research and application of feature line extraction of point cloud[D]. Taiyuan: North University of China, 2018. )