基于密集匹配点云的露天矿台阶线提取方法研究

王继野

(中铁十九局集团矿业投资有限公司,北京 100161)

露天矿台阶线是矿山开采过程中产生的重要边坡迹线,是露天矿测量的主要对象之一。传统的露天矿台阶线提取主要是利用全站仪[1]、GPS[2]等测量方法。随着测绘科学技术的发展,出现了利用数字高程模型(Digital Elevation Model,DEM)[3]、三维激光点云数据、无人机倾斜摄影测量[4]等技术来提取露天矿台阶线的现代测量方法。其中,基于激光扫描与摄影测量技术的台阶线提取方法主要是利用现代化测绘手段构建三维点云模型,并在点云模型中提取台阶特征线。

目前主要的点云特征线提取方法可以归纳为3种:① 基于点云曲率或者法向量[5]。CHEN等[6]根据向量分布和聚类的性质提取点云特征点,采用三次B样条曲线拟合算法得到点云特征线。该方法适用于规则物体特征线的提取,对于一些地形起伏的场景点云,该方法难以提取出相应的特征点。② 基于邻近点云特征相似性。FU等[7]提出了一种基于点空间几何结构的特征线提取方法,提高了点云特征线提取的准确性和效率;WANG等[8]提出了一种利用激光线结构光获取三维数据的特征提取方法,并被应用于铁路特征线提取。总体来看,基于邻近点云特征相似性提取点云时,对点云颜色等信息要求较高,点云间相似性度量标准的选择难以把握,故而较难广泛使用。③基于点云分割的特征线提取。NI等[9]基于区域增长和模型拟合的混合方法实现点云特征线跟踪;GILANI等[10]通过点云分割的方式完成了建筑物特征线提取,解决了在遮挡和低采样点情况下的特征线提取难题。该类基于点云分割的特征线提取操作需要在分割的基础上实现,对于一些表面不规则的物体分割较为困难,难以保证特征线的提取精度,因而该方法应用具有一定的局限性。

上述特征线提取方法虽然可以较好地完成点云特征线提取,但是大部分主要是针对工业零部件、建筑物点云,该类点云模型形状较为规整,特征线以直线为主。目前研究中较少针对露天矿区不规则台阶线提取方法的研究,东北大学王植等[4]曾提出一种顾及邻域几何属性的三维边缘检测与曲率指数加权方法,实现了台阶线的自动化提取,但提取出的台阶线出现了少量不连续现象。由于露天矿台阶线具有不规则性且矿区地形复杂,因而利用点云数据进行露天矿台阶线提取具有一定的难度。鉴于此,本研究首先基于自适应邻域内高程变化特征提取台阶线特征点,然后采用三次B样条曲线拟合出露天矿台阶线,取得了较好的提取效果。

1 露天矿台阶线提取原理

1.1 点云局部特征判定

对于点云样本数据集而言,空间维度为3,此时,点云样本数据集S= { s1,s2,…,sj,…,sn}⊂R3,j=1,2,…,n,n为样本点数量。取数据集S中任意样本点s,其K近邻域点集Ns= { s1,s2,…,sb,…,sk},b=1,2,…,k,点s的协方差矩阵可表示为[4]

由于协方差矩阵M为对称正定阵,因此其特征值λ0、λ1、λ2具有非负性,同时可用协方差矩阵M对应的特征值对点周边邻域的维度特性进行描述,即:① 当 λ0＞＞λ1、λ2＞0 时,局部邻域呈线性特征分布;② 当 λ0、λ1＞＞λ2＞0 时,局部邻域呈平面特征分布;③ 当 λ0≈λ1≈λ2时,局部邻域呈三维曲面特征分布。

1.2 自适应邻域内高程变化的特征点提取

露天矿台阶线包括台阶坡顶线和台阶坡底线[11]。如图1所示,台阶上部平盘与下部平盘均是一个平面,平盘内各处的高程值大致相等,台阶坡面高程呈现渐变趋势。依据该特征可以进行台阶特征线提取。

图1 露天矿台阶组成示意Fig.1 Schematic of the composition of the open-pit steps

如图2所示,假设点P为台阶线中某一点,P1、P2,…,PN是以点P为中心、r为半径的邻近点。记邻近点数量为N。理论上,邻近点P1、P2,…,PN中高程(z分量)与中心点P高程相等的点的数量约占邻近点数量N的50%。但是平盘地形起伏以及台阶坡面隆起凹陷将会造成以下影响:① 位于平盘上的邻近点高程值并非处处相等。位于平盘上的邻近点P1、P2,…,Pn(图2中蓝色点)与中心点P的高程(z分量)之差小于某一高程差阈值Th,即 zPi-zP≤Th(zPi、zP分别表示邻近点Pi与中心点P的高程,i=1,2,…,n,n＜N),同时记满足该条件的点(平盘上的邻近点)的数量为n(n＜N)。② 平盘上邻近点的数量n与总邻近点数量N的比值在50%左右波动,即n/N∈ ( 0.5-l,0.5+l ),其中l为波动幅度。

图2 台阶线特征点示意Fig.2 Schematic of the characteristic points of step line

确定出高程差阈值Th以及波动幅度l后,即可对台阶线进行提取。如图3所示,红色点为中心点,r为搜索半径,中心点上下两条点划线表示到中心点高程差在Th范围的分界线,黄色点表示近邻点云中到中心点高程差小于Th的离散点。由图3可知:对于平坦区域(位置一)和陡峭区域(位置三)选用同样的阈值,会将位置一处的中心点误提取为特征点。

图3 台阶线特征点提取示意Fig.3 Schematic of the feature point extraction of step line

为此提出了一种基于自适应邻域内高程变化的特征点提取方法,依据前文1.1节中判定出的局部点云的特征维度选用不同的阈值。基本步骤如下:

(1)选定某一点P作为中心点,采用K近邻搜索(K Nearest Neighbor Search)法搜寻到中心点距离在半径r(本研究中半径r取2 m)以内的邻近点Pi(i=1,2,…,N,N为搜索到的邻近点数量)。

(2)依据式(1)判断中心点的邻近点云的局部邻域维度特性,按照中心点局部邻域特性选择不同的高程差阈值Th和波动幅度l。

(3)遍历点P的邻近点,比较邻近点云中各离散点Pi的z分量与中心点P的z分量之差,记满足zPi-zP≤Th的点的数量为n;计算n与N的比值,判定该中心点是否为台阶线特征点。

(4)以点云中各个离散点为中心点,循环重复步骤(1)、(2)、(3),直至执行完点云中所有离散点。

1.3 台阶线平滑

由于台阶平盘局部区域存在凹陷和突起现象,使得提取到的台阶线特征点中存在误提取现象。在特征点平滑连接之前需要剔除误提取的特征点。首先采用欧氏距离分割算法[12]对点云进行聚类,其次在完成聚类分割后计算各个点云簇中离散点的数量Ni,同时设定一个数量阈值Nthresh。若Ni＜Nthresh,则对该类进行剔除;反之,予以保留。本研究中数量阈值Nthresh经试验验证最终设为200。在试验中发现经剔除后的特征点中不仅包含有台阶线特征点,还含有矿区道路边线点云。

由于矿区车辆运输等实际需要,进出采坑的道路均具有一定的坡度;而露天矿平盘是一较为规则的水平面,地形变化较小。采用点云的高程值作为区分道路边线与台阶线特征点的依据。对于经欧氏距离分割得到的某一簇道路边线特征点Ω1而言,其高程值随坡度不断变化,该簇点云的高程值离散程度偏大,极差偏大;对于经欧氏距离分割得到的某一簇台阶线特征点Ω2而言,该簇点云的高程值离散程度较小,高程极差较小。研究中手动提取出了多簇道路边线特征点与台阶线特征点,经计算得到如下规律:当点云簇中各点高程值的方差小于0.75且极差值小于2 m时,判定该点云簇为台阶线特征点;反之,为道路边线特征点。

利用上述规律对提取到的特征点进行进一步筛选,将筛选后的特征点作为最终提取的露天矿台阶线特征点。

台阶线是由提取到的无数个特征点按序连接而成的。但是上述方法提取到的特征点是无序的,因此需要对特征点进行排序[13]。本研究主要依据离散点间距与偏转角进行点云排序。

如图4所示,以P1点作为初始中心点,依据最近点优先原则,P2点为P1点的连接点;以P2点作为中心点,P3点为P2点的连接点;以P3点作为中心点时,由于P3到P4的距离与P3到P5的距离相等,即P3的最近点为P4和P5点,此时最近点优先原则失效,改用偏转角来确定P3的连接点。从图4中可以看出,P1P2和 P2P3之间的夹角 θ1为顺时针转角;P2P3和P3P4之间的偏转角为θ2,且θ2为顺时针转角;P2P3和P3P5之间的偏转角为θ3,且θ3为逆时针转角。由图4可进一步看出,θ1与θ2转向一致,大小相近;而θ1与θ3转向相反,大小相差较大,故将P4点作为P3点的连接点。

图4 离散点排序规则[13]Fig.4 Sorting rules for discrete points

经排序得到有序点云后,采用三次B样条函数对特征点进行平滑连接处理[14]。三次B样条曲线的方程为

式中,Fi,k(t)为基函数,即:

因此,每4个离散点就可拟合一段曲线,例如P0、P1、P2、P3点可以拟合一段光滑曲线,P1、P2、P3、P4点可以拟合下一段,相邻两段曲线是平滑过渡的,以此类推,N个点可以拟合出(N-3)段平滑相接的曲线。

2 试验结果与分析

2.1 数据获取及预处理

本研究使用的数据是由大疆Phantom 4 RTK无人机采集的鞍山齐大山铁矿以及鞍千矿的影像数据。在Smart 3D软件中对获取到的影像进行空中三角测量、多视影像密集匹配、纹理映射等操作得到研究区的密集匹配点云模型[15],如图5所示。受到无人机拍照角度的影响,所构建的密集匹配点云中各处的点云密度存在差异性。利用CloudCompare软件计算出构建的点云平均点密度约为48个/m2,相邻点间最大高程差为0.048 9 m。

图5 研究区三维点云模型Fig.5 3D point cloud model of the study area

本研究使用C++语言与PCL点云库编写算法程序,计算机配置参数见表1。

表1 试验所用的计算机配置Table 1 Computer configuration for the test

2.2 参数确定

由前文1.2节可知,本研究方法中影响露天矿台阶线特征点提取精度的参数有高程差阈值Th和波动幅度l,且在不同地形条件下这两个参数的取值应是不同的。为了确定合理的参数取值,本研究基于鞍山齐大山铁矿数据集对参数合理取值进行试验分析。

2.2.1 高程差阈值Th选取

在保持区间波动幅度l恒定(波动幅度l=0.05)的情况下,确定高程差阈值Th的取值。由于所构建的密集匹配点云中相邻点间最大高程差为0.048 9 m,对该值进行取整,以0.05 m为间隔选取Th,分别取Th=0.05、0.10、0.15、0.20、0.25、0.30 m,不同高程差阈值下的露天矿台阶线特征点提取结果见图6。由图6可知:当Th取0.05 m时,研究区内地形平坦区域存在着大量误提取的点云;随着Th的增大,平坦区域内误提取的点云逐渐减少。当Th取0.10 m时,地形平坦区域中误提取的点较图5(a)中明显减少;当Th取0.30 m时,由图6(e)可见平坦区域内几乎没有误提取点。但是在Th增大的过程中,台阶上(下)部平盘中一些突起区域中逐渐出现了被提取出的点云,然而这些点云并非台阶线特征点,如图6(d)、图6(e)、图6(f)中椭圆标注的位置。从图6(c)中可以看出,台阶坡底线特征点提取结果优于Th=0.10 m和0.20 m的提取结果。

为了保证台阶线特征点的提取精度,依据前文1.1节中邻近点协方差矩阵的特征值进行了高程差阈值Th选取,地势平坦区域Th取0.30 m,地形复杂区域Th取0.15 m。

2.2.2 波动幅度l选取

利用图像构建出的密集匹配点云中相邻点间最大高程差为0.048 9m。随机选取了数据集中5个点作为中心点,同时分别搜索到与这5个中心点的距离为2m的邻近点,记录邻近点与中心点高程值之差小于(0.048 9t)m的离散点数量T,其中t=1,2,…,k且满足 0.048 9k≤2。然后,计算高程差每增加0.048 9T m的增长比例。

以选取的5个中心点为研究对象,试验发现高程差每增加0.048 9 m,满足条件的离散点数量T平均增长比例分别为 2.06%、3.64%、3.32%、2.51%、2.93%。为了保证提取的精度以及计算效率,本研究中将2%作为波动幅度l的调整幅度,即l=0.02k,且满足0.02k≤0.5。图7为部分波动幅度l取值的特征点提取结果。

图7 不同波动幅度l取值下的特征点提取结果Fig.7 Feature point extraction results under different fluctuation amplitude l values

由图7可知:随着波动幅度l的增大,提取到的特征点不断增多,特征点逐步从台阶坡线位置向台阶坡面上移动,开始远离台阶线所在的实际位置。当波动幅度l=0.12时,提取到的非台阶线特征点增多,如图7(d)所示;当波动幅度l=0.06时,提取到的特征点存在断续现象,漏提取的特征点较多。当波动幅度l=0.10时提取效果较为理想。

综上所述,高程差阈值Th及波动幅度l的取值如下:当离散点对应的邻域点云协方差矩阵M的特征值 λ0、λ1、λ2存在 λ0＞＞λ1、λ2＞0 或 λ0、λ1＞＞λ2＞0关系,即局部邻域呈线性或平面特征分布时,高程差阈值Th取0.30m,波动幅度 l=0.10;当 λ0、λ1、λ2存在λ0≈λ1≈λ2关系,即局部邻域呈三维曲面特征分布时,高程差阈值Th取0.15 m,波动幅度l=0.10。

2.3 台阶线提取结果

在确定出本研究算法中的两个重要参数后,对研究区进行了台阶线特征点提取,结果如图8所示。图8(a)是利用统一参数标准提取的特征点,该方法中未考虑地形起伏的影响,整个研究区内的参数保持一致,从提取结果中可以看出存在着较多误提取的离散点。图8(b)为自适应选用参数提取到的特征点,从中可以看出平坦位置或多或少还存在着误提取的特征点,但是较图8(a)相比误提取的点数量明显减少,提取效果得到了明显改善。

图8 特征点提取结果对比Fig.8 Comparison of the feature point extraction results

完成特征点提取后根据道路边线特征对特征点进行筛选,剔除道路边线点。最后采用三次B样条曲线对筛选的特征点进行曲线拟合得到台阶线。图9为最终提取到的两组数据集露天矿台阶线,从图中可以看出提取结果较好,可有效区分矿区的台阶线与道路边线。

图9 露天矿台阶线提取结果Fig.9 Extraction results of step line in open-pit mines

在提取效率方面,将该方法与传统的GPS点测量方法进行对比。对于1 km2作业区来说,使用一台GNSS(Global Navigation Satellite System)接收机每5 m采集一个碎部点的方式进行GPS点测量,其外业采集时间约15 h,内业数据处理耗时约2.5 h;利用本研究方法进行台阶线提取时,飞马D2000无人机航飞时间约22 min,内业处理时间约4.5 h。从时间上来看,本研究方法减轻了外业工作量,同时极大地提高了工作效率。

3 结 论

(1)基于无人机密集匹配点云数据,提出了自适应邻域内高程阈值的露天矿不规则台阶线提取方法。试验结果表明:相比采用统一参数提取到的特征点,该方法能有效减少误提取的特征点数量,并可区分出台阶线特征点与道路边线特征点,较好地提取出露天矿中不规则的台阶线。与传统台阶线提取方法相比,在减轻外业工作量的同时,显著提升了台阶线提取效率。

(2)目前该方法在部分区域由于阈值选取的偏差导致提取到的台阶线出现了错位现象,后期需要进一步进行方法优化。