煤矿巷道三维激光扫描关键技术及工程实践

王海军，刘再斌，雷晓荣，韩保山，陆自清

（中煤科工集团西安研究院有限公司，陕西 西安 710077）

近年来，随着工业互联网、大数据、云计算、人工智能、5G 等技术的飞速发展及在各行各业的广泛应用，传统的煤炭行业迎来了智能化、无人化发展的巨大机遇并取得了一系列成果。但是，必须清醒地认识到煤矿智能化开采距离真正的智能化、无人化还有很长的路要走，中国工程院王国法院士指出：智慧煤矿与智能化开采是一个复杂的巨系统，精准地质信息探测是当前智慧煤矿建设中的重点研发方向之一[1-2]。

巷道承载了矿井运输、通风、排水、供电、逃生路径等功能，巷道三维模型的快速获取成为地质透明化的重要一环，是煤矿智能化高效绿色开采的重要组成部分。传统巷道建模主要有3 种：基于巷道顶底板中心线和断面建模；基于巷道中心线和顶底板边界线建模；基于巷道中心线或者边界线和高程建模。一般步骤：采集巷道测量数据（实测数据、断面数据）、数据去重、规格化、模型构建、表面三角化、实体渲染[3]。其中最重要的是采集巷道测量数据，一种方法是利用井下巷道的导线点并结合巷道的地质写实图获取数据并建模，另一种方法是组织人力进行巷道尺寸、煤层起伏情况等数据测量。第一种方法的缺点是随着矿井采掘活动的进行，巷道会产生一定的形变，导致地质写实、图纸更新和建模需要定期重复；第二种方法的缺点是井下地质写实和测量费时费力，效率低下。

三维激光扫描技术可以快速获取被测物体的表面三维坐标，扫描得到的点云数据精度高、密度大、包含反射强度等丰富的语义信息。借助点云数据可以构建复杂巷道模型。一般通过三维激光扫描系统进行数据采集，该系统集成了全球卫星导航系统、惯性测量单元、一颗或多颗激光雷达、全景相机、里程计等传感器，安装在飞机、车辆、移动小车等不同的载体上，可实现不同场景三维点云数据的快速获取。受到卫星定位信号及精度的影响，在信号较差环境下（如地下停车场、隧道、煤矿井下等）点云数据的采集及拼接是业界的难题。

为此，很多学者针对这些特殊应用场景进行研究。江记洲等[4]采用圆柱面投影将三维点云转换为二维离散点，用分治算法进行三角剖分，结合二三维点云及三角网之间的拓扑关系重建三维巷道模型；石信肖等[5]针对三维激光扫描技术获取的点云数据，使用Delaunay 生长算法引入边和角约束条件、设置三角网边长阈值并构建了狭长形海量巷道点云模型；金卓等[6]针对单站点三维激光扫描点云，提出一种基于圆柱形投影面巷道建模方法；张君等[7]将移动式三维激光扫描设备的数据获取、数据处理、模型建立等引入井巷工程快速测绘中，提高了工作效率；郭良林等[8]提出两期单测站点云配准的方法，研究三维激光扫描在井下巷道变形监测中的应用；刘晓阳等[9]采用三维激光扫描方法对巷道顶板稳定性进行监测研究；赵小平等[10]基于三维激光扫描技术和三维GIS 构建三维巷道空间模型的数据库组合法则，实现了巷道数据模型构建。以上研究主要是基于点云数据在小段巷道模型构建、井巷测绘和巷道变形等方面取得了一定的应用效果，但是鲜有涉及大场景、长距离工作面巷道的多站数据校正、扫描、配准和重建等内容。

笔者在分析三维激光扫描方法煤矿井下环境面临的主要难题的基础上，重点研究了三维激光扫描的原理、关键技术和主要难题，设计了三维激光扫描系统的数据处理系统，为复杂巷道的快速扫描和重建提供了一种解决办法。

1 煤矿井下三维激光扫描面临的难题

在煤矿井下环境中，三维激光扫描系统面临的主要难题体现在以下几个方面：

（1）煤矿井下工况环境属于GB 3836.1-2010《爆炸性环境第1 部分：设备 通用要求》规定的爆炸性环境，要求在该场景下使用的电气设备必须具备矿用安全标志证书和防爆合格证。三维激光扫描系统的防爆改造是难点之一；

（2）在煤矿井下环境作业时，三维激光扫描系统自身的位姿必须通过惯性测量单元和里程计确定，由于组成惯性测量单元的陀螺仪存在漂移问题，长时间累积计算会产生较大误差[11]，所以单次点云采集时间不宜过长；

（3）分段测量的多个点云数据如何以较高的精度对准到统一的地理坐标系坐标，对准后大距离点云数据的处理和拼接难题；

（4）煤矿井下巷道距离长（需要分段采集）、煤岩特征相对单一（需要人工布设多个特征靶点）、巷道分支多（拼接困难）、粉尘较大（点云噪点增多）等特点，给点云数据的外业采集、内业处理和拼接带来很大的困难。

2 激光扫描原理及技术流程

2.1 激光扫描原理

三维激光扫描系统的测距原理为飞行时间测量法，即激光发射器发射出一束超短激光脉冲，激光投射到被测物体上，发生漫反射后返回光，激光接收器接收到漫反射光，通过记录激光束在空中的飞行时间，可准确计算被测物体到激光雷达传感器中心的距离d，表示为：

式中：d为距离；c为光速（一般为300 000 km/s）；t激光束飞行时间。

空间某数据点的坐标计算如图1 所示。

图1 空间点坐标计算Fig.1 Calculation of the space point coordinate

式中：X、Y、Z分别为某数据点在仪器坐标系下的坐标值；α为该激光束的垂直角分辨率；β为水平旋转间隔。

2.2 三维激光扫描技术流程

综合考虑煤矿井下巷道错综复杂的特点以及本次作业距离等因素，设计如图2 所示的三维激光扫描技术流程。

图2 三维激光扫描技术流程Fig.2 3D laser scanning technology process

①制定扫描采集方案，利用巷道现有的导线点坐标布设高反反光贴并测绘出相应的坐标点；

② 在每一个布设高反反光贴位置进行三维扫描系统的动态标定，而后开始采集原始点云数据，到下一个导线点位置处停止上一次数据采集并再次进行动态标定和原始点云数据采集，直到设定的采集路径全部采集后进入下一环节；

③对采集的每一站点云数据进行预处理（包含大尺度滤波和小尺度滤波），过滤掉大范围内的空间干扰点云；然后对点云数据的关键特征点进行提取和拼接，拼接处重复的点云采用滤波算法去重；

④ 点云高分重建、边界特征提取等专门应用；⑤ 最后输出标准格式的点云数据。

3 三维激光扫描重建关键技术

3.1 三维激光扫描系统动态标定

3.1.1 面临的难题

三维激光扫描系统获得的点云数据，是以扫描系统激光雷达中心为原点的仪器坐标系。三维激光扫描系统仪器坐标系的具体定义为：原点位于激光光源点（激光发射器），X轴正方向与激光雷达水平安装时0°激光线的发射方向一致，X轴顺时针旋转90°即为Y轴正方向，Z轴垂直于X轴与Y轴组成的平面上为正，如图3 所示。其中，灰色线条表示激光雷达接收和发射阵列各通道在垂直方向激光线束的分布。

图3 仪器坐标系Fig.3 Coordinate system of the instrument

煤矿井下使用的是地理坐标系（1980 西安坐标系或2000 国家大地坐标系），巷道扫描获取的点云也需要转换到地理坐标系。在煤矿井下环境作业时，三维激光扫描系统自身的位姿必须通过惯性测量单元和里程计确定，由于组成惯性测量单元的陀螺仪存在漂移问题，长时间累积计算会产生较大误差。为了获取准确的点云坐标，点云数据单次采集时间不宜过长，而且必须对三维激光扫描系统的坐标进行动态标定。

3.1.2 坐标变换

在煤矿井下具有已知导线点（测量坐标点）的位置作为采集的起始点。如图4b 所示，上部红色点为导线点，三维激光扫描系统放置在导线点的正下方，巷道两旁布置2 个高反标定物，并拿全站仪测量高反标定物中心的坐标。在激光扫描开始工作前，计算仪器采集的高反标志物的中心坐标与全站仪测量的高反标定物中心坐标，得出动态标定系数，其本质是两个三维空间坐标的3 个平移和3 个旋转运算，尺度变换参数取1，空间三维坐标变换原理如下：

式中：(x0,y0,z0)T为两坐标之间的平移量；A=AZ(φ)AY(θ)AX(ϕ)，为绕X轴旋转ϕ、绕Y轴旋转θ和绕Z轴旋转φ的3 个旋转矩阵的乘积，(xd′,yd′,zd′)T为初始坐标，(xd,yd,zd)T为目标坐标。

旋转角顺时针为正，逆时针为负[12]。

旋转矩阵A可改写为单位四元组元素表示，求解并最终得出平移量和旋转量[13]。

空间坐标转化为地理坐标的变换模型如下：

式中：L为大地经度；B为大地纬度；H为大地高程；a为椭球的长半轴，6 378.137 km；N为椭球的卯西圈曲率半径；e为椭球的第一偏心率，；b为椭球的短半轴，6 356.752 km；Φ、R为中间变量，其表达式如下。

利用式（7）计算时有交叉变量，先求出B的初值，带入式（7）求出H、N的初值，再次求出B的值。

3.2 点云预处理

1）大尺度噪声滤波

大尺度噪声指在主体点云周围，偏离主体点云且空间分布较稀疏的点或距离主体点云中心较远的点云，一般认为距离主体点云位置较远的点云所含的特征信息不丰富，可以考虑为噪点。采用统计滤波[14]去除大尺度噪声，主要步骤如下。

①对每个点云数据每个点的邻域进行统计分析，假设点云中所有点的距离构成高斯分布，其形状由均值μ和标准差ε决定。假设第n个点的坐标为Mn(Xn,Yn,Zn)，该点到任意一点Pm(Xm,Ym,Zm)的距离为：

② 遍历每个点到任意点Pm(Xm,Ym,Zm)之间距离平均值和标准差公式分别为：

③假设标准差的倍数为k，则计算过程中只需要输入点数n和k两个阈值即可，当某个点距离n个点的平均距离在标准差范围（μ-εk,μ+εk）内时保留，反之定义为噪点删除。

2）小尺度噪声滤波

小尺度噪声指混合在主体点云中的噪点，会影响点云模型的平整度，造成点云模型失真。移动最小二乘平滑滤波主要以重采样点偏差的平方和最小，对给定点集定义一个隐式全面，将采样点局部区域中的点投影到该曲面上，并用高阶多项式逼近以实现采样点的平滑和去噪[15]。移动最小二乘法计算原理如下：

已知节点X=[x1,x2,···,xu]T（xi∈Rr，u为节点数，r为节点空间维数），对应的节点值Y=[y1,y2,···,yu]T，求解拟合函数f(x)使得：

式中：w(·)为xi的权函数，x-xi即为节点间的欧式距离，权函数是一个减函数，随距离增加而减小，权函数具有紧支性，即在xi的影响域外权函数为零。推导后拟合函数可表示为：

式中：

3.3 点云关键点提取与特征描述

在三维点云空间关键点提取与特征描述时，除了考虑常规的坐标和特征外，需要考虑法线方向、曲率、纹理特征等描述子。点云关键点的数量相比于原始点云减少很多。与局部特征描述子结合在一起，可以完整地描述整个点云且不失代表性和可描述性。综合考虑到特征描述子应具备旋转不变性、较强表征力和不包含任何手工特征等因素，文中点云关键点提取选用尺度不变特征转换SIFT（Scale Invariant Feature Transform），特征描述选用快速点特征直方图算法FPFH（Fast Point Feature Histograms）。

1）SIFT 特征检测算法

SIFT 特征检测算法是D.G.Lowe[16]在1999 年提出的一种局部特征描述算法，2004 年完善后引入到三维点云中，实现了对点云特征的提取，具有对亮度变化、噪声、旋转和平移等因素保持较好的不变性。SIFT 算法的主要步骤如下。

①构建尺度空间。二维图像的尺度空间函数定义为：

其中，G(x,y,σ)是图像I(x,y)的可变高斯函数：

式中：(x,y)为图像I上某点的坐标；S表示尺度空间；σ为尺度空间因子，其值越小表示图像越平滑。

同一阶上两个相邻层的尺度函数相减得到高斯差分金字塔：

式中：k′为同一阶上两个相邻层之间的尺度比例。

② 尺度函数空间内检测极值点。将每一个采样点与同层、相邻层的所有相邻点进行比较，得到局部极值点（最大值或最小值）并作为下一个候选关键点。

③特征点过滤及定位。对局部极值点进行三维二次函数拟合可以得出特征点的位置和尺度，尺度函数进行泰勒展开并求偏导，去掉低对比度的特征点和不稳定的边缘点后得到极值点的位置：

④ 确定关键点方向值。利用梯度直方图统计领域像素的梯度方向，梯度直方图的主峰值为关键点的主方向。

⑤ 建立包含尺度、位置、方向等信息的特征描述子。

2）FPFH 特征描述算法

点特征直方图PFH（Point Feature Histograms）算法是R.B.Rusu 等[17]于2008 年提出的基于特征点与其邻域点的空间几何关系来编码的特征描述算法，n个点云计算PFH 特征的时间复杂度为O(nk2)，效率较低。在保留PFH 算法核心思想的基础上R.B.Rusu等[18]于2009 年又提出了快速点特征直方图（FPFH）特征描述子，计算时间复杂度降到O(nk)。该算法的主要步骤如下。

①构建点对坐标系。对于每一个组点对ps和pt，建立局部坐标系：

式中：×表示外积；h为点的法向。

② 计算特征算子。对于点云模型中所有点首先设定半径r邻域，计算α、η和γ三个特征算子：

式中：·为内积；d为ps和pt之间的欧式距离。

③特征编码。对于每个点pi计算包含两倍半径r内的其他点。查询点pq与周围c个点组成的点对，于新包含点的SPFH 进行加权并与点pq本身的SPFH求和，最后得到目标点pq的FPFH，公式如下：

式中：wi为点pq与近邻点pi的欧式距离；c为pq邻域内的近邻点个数。

3.4 点云配准

点云配准一般分为粗配准和精配准两部分。

1）粗配准

粗配准指在点云相对位姿未知的情况下对点云进行配准，可以为精配准提供较好的初始值。采用FPFH 特征描述算法，得到初始的变换矩阵，将源点云配准到目标点云上，具体步骤如下。

①预先设置阈值δ′，在源点云ps中选择几个特征点并确定这些点的空间距离大于阈值δ′，保证每个特征点的FPFH 特征不同。

② 依据特征点的FPFH 特征，在目标点云pt中找到相似的FPFH 特征一个或多个点，并从这些相似点中随机选取至少3 个点作为源点云在目标点云中的一一对应点。

③计算对应点间的刚体变换矩阵。先求出源点云与目标点云间的变换关系，然后依据该变换关系计算对应点变换后的距离误差和函数，将此函数作为评价配准性能指标。

式中：ml为预先设定的值；li为第i组对应点经过变换后的距离差。

粗配准后源点云与目标点云间的旋转和平移误差缩小，获得较好的初始位置。

2）精配准

迭代最近点ICP（Iterative Closest Point）算法[19-20]是点云精配准使用最多的方法之一。基本原理是：在待匹配的源点云集合Q 和目标点云集合P 中，按照一定的约束条件，找到最邻近点（pi,qi），然后计算出最优匹配参数R和T，使得误差函数E(R,T)最小，即：

主要计算步骤如下。

①计算最近点集。分别在源点云集合Q 和目标点云集合P 中，使用高维索引Kd-Tree（Kd 树）找出近邻点pi和qi，使得‖qi-pi‖最小。

② 分别计算旋转矩阵R和平移矩阵T。

③对源点云集合Q 中的点qi，应用旋转矩阵R和平移矩阵T求解出新的点云集：

④ 计算qi′与对应点集pi的平均距离：

如果d′′小于给定的阈值或大于设定的最大迭代次数，则停止迭代返回第②步。

⑤ 经过多次迭代计算，得到最优的旋转矩阵R和平移矩阵T，实现点云的精确配准。

4 工程实践

4.1 施工位置

唐家会煤矿地处准格尔煤田中部，井田面积大约28.573 km2，可采煤层4、5、6、9上和9下共5 层，主要开采6 煤，资源总储量7.66 亿t，核定生产能力900 万t/a。6 煤以暗煤、亮煤为主，大部含夹矸4～6 层，厚度13.723～22.716 m，平均厚度18.363 m。某工作面位于6 煤西南部，走向长度1 590 m，宽240 m。

4.2 扫描准备

扫描路径规划：从副立井开始，沿着6 煤辅运大巷、6 煤南辅运大巷、某工作面回风巷道、切眼、某工作面运输巷道、6 煤南回风大巷辅运联巷、6 煤南辅运大巷，到6 煤南辅运大巷与某工作面回风巷道交汇处终止，总扫描距离大约5 625 m。总导线点44 个,每个导线点布设一个高反反光贴组合（由两片单独的反光贴拼成一个直角组合），尺寸：30 mm×5 mm，进口强力黏胶。导线点和反光贴组合的布设如图4 所示。

图4 中顶部红色点表示导线点，巷道侧帮蓝色部分表示高反反光贴组合，反光贴组合顶部的红色位置表示2 个反光贴拼接后的直角的顶点位置。

图4 施工前准备Fig.4 Preparation before construction

4.3 巷道扫描

巷道激光扫描的难点在于井下巷道环境相对复杂。对比轮式激光扫描小车、无人机载三维扫描系统和移动式三维激光扫描系统的优缺点后，本文采用中煤科工集团西安研究院有限公司自主研发的移动式三维激光扫描系统进行数据采集。该系统由多线激光雷达、旋转云台、里程计、惯性单元等组成，主要技术参数见表1。

表1 移动式三维激光扫描系统主要技术参数Table 1 Main technical parameters of mobile 3D laser scanning

扫描、处理和拼接完成的点云模型如图5 所示。

图5 巷道点云成果Fig.5 Results of roadway points cloud

4.4 轮廓线提取

煤矿井下巷道轮廓线提取相对简单，主要步骤如下：

①把巷道点云模型整体做顶透视，然后把三维点云模型投影到水平面，设置内部和外部轮廓线宽度为单点大小，并批量导出内部和外部轮廓线的数据点坐标；

② 复制所有的轮廓线坐标数据，打开AutoCAD软件，点击多线段命令，在命令提示栏中粘贴所有坐标并空格确认输入即可。

提取的巷道点云轮廓（红色）和巷道（黑色）对比如图6 所示：

图6 巷道点云轮廓Fig.6 The contour of roadway points cloud

从图6 可以看出，点云模型提取的轮廓线与巷道平面图叠加后的整体吻合度较好，局部位置稍微有偏差（两导线点较远时）。

4.5 巷道扫描与工作面联合建模

采用4.4 节中点云模型投影到XOY平面提取的轮廓线的数据点坐标，过该坐标点作平行于Z轴的切片并采用软件系统配套软件测量出巷道纵向高度。在图6 中提取出巷道关键位置的坐标信息并与Auto-CAD 图纸构建的巷道模型进行局部修正，较精确地构建巷道模型。在此基础上，将巷道模型和扫描后的点云巷道模型融合可以实现待回采工作面巷道的快速测量（图7），有助于工作面地质模型的快速动态更新。

图7 巷道与工作面联合建模Fig.7 Joint modeling of the roadway and working face

图7 中，黑色部分为工作面地质模型，灰色部分为传统手段构建的巷道模型，彩色部分为巷道点云模型，红色部分为长度标尺。

4.6 精度分析

本次应用共44 个导线点（即44 站测量数据），44 站点云数据进行拼接时，每两站数据之间的重叠度为20%。点云拼接精度用局部纵向切割和整体偏差进行表示。

1）局部纵向切割精度

采用4.4 节中投影到XOY平面提取的轮廓线的数据点坐标，间隔5 m 抽取1 个样本点过该坐标进行纵向切割，测量点云数据巷道切割后剖面图中巷道顶底板之间的距离（单个点的纵向切割如图8 所示），共抽取5 个样本点（样本点越多，所取的巷道越长，对巷道中同一个样本点与点云定位越困难）与巷道顶底板测量距离进行对比分析，并给出偏差值（偏差值的计算方法为点云成果切片后巷道高的测量值与实际巷道高的测量值的差），结果见表2。

表2 纵向切割精度Table 2 Longitudinal cutting accuracy

图8 巷道切片Fig.8 Roadway slicing

从表2 可以看出，点云成果切片后巷道高度的测量值与实际巷道高度的测量值的差绝对值最大为0.027 m，最小为0.007 m，该偏差在移动式三维激光扫描系统中激光雷达的测量精度范围内。移动式三维激光扫描系统中激光雷达的极限测量距离为30 m，剔除煤矿井下粉尘、水雾等影响因素，实际测距10 m 处能较清晰地扫描到被测物体。因此，该系统适用于矿井所有巷道断面高度方向的扫描。

2）整体偏差

把点云模型投影到XOY平面后与巷道平面图叠合后进行对比分析，把二者的偏差距离按照测量起始点到5 625 m 对应长度进行分段统计，结果见表3。

表3 整体偏差Table 3 Overall deviation

从表3 可以看出，随着点云数据拼接长度的增加，点云巷道模型与实际巷道投影后的平面图重叠度越来越低，整体偏差越来越大。巷道距离超过5 400 m 时，偏差相对较大，超过了业界认可的0.3 m 以内的要求。

5 结 论

a.针对煤矿井下特殊的工况环境，设备防爆、陀螺仪惯性导航漂移、坐标对准和长距离拼接等是煤矿井下三维激光扫描面临的施工和技术难题。在此基础上，论述了三维激光扫描系统核心部件三维激光雷达的测距原理，包括飞行时间测量法的原理和空间坐标的转换公式；提出煤矿井下长距离巷道三维激光点云数据处理的技术实现流程。

b.煤矿井下长距离巷道三维激光点云扫描需要重点解决的关键技术，包括三维激光扫描系统动态标定、基于大小尺度噪声滤波的点云预处理技术、点云关键点提取与特征描述技术、基于粗配准和精配准的点云配准技术等。

c.采用自主研发的移动式三维激光扫描系统在准格尔煤田某矿工作面进行试验，验证了设备性能、施工流程和点云数据处理流程的可行性，为煤矿智能化开采巷道信息的精准探测和巷道三维模型的快速获取探索了一条可行的技术路径。

d.移动式三维激光扫描系统在扫描巷道距离5 400 m 内时整体拼接偏差能满足实际使用。因此，下一步研究的重点是提高大场景、长距离点云巷道模型的整体拼接精度。