磷矿作为化工行业的关键原料，具有不可替代的战略意义，是保障化学工艺及产品生产的基石。在这一背景下，磷矿的安全开采显得尤为重要，合理的安全管理措施不仅能有效降低开采风险，还能提高生产效率，减少资源浪费。矿山巷道在磷矿开采过程中担负着矿井运输、通风、行人等重要任务。随着开采向深部延伸，主要巷道受压力增大，开始变形破坏[1-2]，继而发生冒顶、片帮乃至崩塌的危险，严重威胁着井下生产活动的开展与工作人员的生命安全。因此十分有必要对巷道的形变进行检测，以了解巷道围岩的稳定性[3-4]。现有的检测仪器有全站仪、手持式扫描仪等。然而，全站仪每次只能对单点进行检测；手持式扫描仪需要人工进行手持扫描，在存在风险的巷道内作业会危及扫描人员的生命安全。
近年来，随着激光雷达技术的国产化，激光雷达设备在提供更广泛扫描范围和更高精度的同时，其体积也实现了小型化[5-6]，因此其应用领域得到了进一步拓展。通过激光雷达对巷道进行扫描，能够采集到高精度并且包含丰富三维坐标的点云数据，为形变检测提供了基础数据支持[7-8]。
本文使用二维激光雷达作为获取主要数据的传感器，将其搭载在可移动载具平台上，并融合姿态传感器的位姿数据，构建了一种能够扫描并检测巷道形变的巷道形变检测系统。该系统不仅能够高效、准确地完成巷道形变的检测，还能在复杂的矿山环境中实现灵活操作，为矿山安全管理提供实时、可靠的数据支持。通过该系统的应用，能够采集完整的巷道三维点云数据，并且对巷道发生的形变位置与程度进行识别，为矿山的安全预警和风险评估提供科学依据，进而保障矿工的生命安全。
1 巷道形变检测系统硬件设计
巷道形变检测系统由硬件和软件两部分组成。硬件部分负责各个传感器的供电、数据采集和数据传输；软件部分负责将采集的传感器数据进行融合与处理，并完成形变分析。硬件部分由二维激光雷达、惯性测量单元（inertial measurement unit，IMU）、编码器、数据采集与处理主机、供电模块，以及搭载上述硬件设备的机器人小车组成，其硬件结构如图1所示。
<G:\武汉工程大学\2025\第3期\陈江-1.tif>[供电][数据传输][数据采集与
处理主机][供电][供电模块][机器人小车][二维激光雷达][编码器][IMU]
图1 巷道形变检测系统硬件结构图
Fig. 1 Hardware structure of tunnel deformation
detection system
激光雷达的选型为星秒科技的PAVO2-LS-50H二维激光雷达，近距离（0.1~15.0 m）测距精度为±10 mm，最小角分辨率达0.08°，激光雷达实物如图2（a）所示。IMU采用亚博智能CMP10A型号的IMU模块，内部集成了姿态测量、加速度测量、海拔测量和地磁信息等功能，内置的姿态解算算法能够快速求解出模块当前的实时运动姿态。IMU实物如图2（b）所示。编码器采用500线AB相高精度编码器，用两个编码器分别采集小车左右轮的里程数据，采用STM32芯片作为主控读取编码器的里程数据，再通过串口将数据发送给上位机，编码器模块如图2（c）所示。二维激光雷达垂直安装在机器人小车的前端，用于获取巷道截面的点云数据。编码器和车轮连接以获得设备的里程数据。IMU安装在机器人小车内部，获取设备的姿态数据，数据采集与处理主机则负责收集和处理来自各个传感器的数据，进行形变识别与体积计算，系统整体硬件如图2（d）所示。
2 数据采集和处理
巷道形变检测系统的软件算法主要分为数据采集和数据处理两个部分。在现场采集数据的过程中进行多传感器的数据融合，从而生成巷道三维点云数据；对获取的巷道三维点云数据依次进行点云数据预处理、形变识别与分割以及体积计算。
2.1 数据采集
在系统进行数据采集的过程中，为了使激光雷达采集的二维点云最终变换成三维点云，需要通过坐标变换和数据融合的方法，将左右编码器产生的车辆里程数据和IMU的姿态数据与二维点云进行关联，生成一帧带有位姿信息的二维点云，并将带有位姿信息的二维点云数据统一在同一公共的坐标系中。此过程中涉及的坐标系有：世界坐标系、小车运动学中心坐标系、IMU坐标系、编码器坐标系和激光雷达坐标系，上述坐标系均采用笛卡尔右手坐标系。世界坐标系的中心与小车初始位置的运动学中心坐标系重合，世界坐标系作为点云数据统一的公共坐标系，系统在数据采集过程中每得到一帧带有位姿信息的二维点云就会添加到世界坐标系中。以小车运动学中心为参考坐标系，两个编码器坐标系、IMU坐标系、激光雷达坐标系的坐标原点相对小车运动学中心的位置分别为[(-x1,y1,-z1)]、[(-x1,-y1,-z1)]、[(x2,0,z2)]、[(x3,0,z3)]，二维激光雷达的扫描范围是激光雷达坐标系的[YOZ]平面（激光雷达坐标系中向上的坐标轴与向右上的坐标轴组成的平面）。小车上各传感器坐标系相对位置如图3所示。
<G:\武汉工程大学\2025\第3期\陈江-3.tif>[IMU坐标系][激光雷达坐标系][编码器坐标系][小车运动学
中心坐标系]
图3 各传感器坐标系示意图
Fig. 3 Schematic diagram of each sensor coordinate system
二维激光雷达在进行扫描时，系统从二维激光雷达的驱动中获取当前时刻目标相对于二维激光雷达的旋转角[α]和极坐标半径[r]，并将其转换为笛卡尔坐标来获得二维点云数据[9-11]，二维激光雷达点云数据的笛卡尔坐标表示为[(xlidar,ylidar,zlidar)]，具体转换方式如式（1）所示。
[xlidar=0ylidar=r×cosα·π/180zlidar=r×sinα·π/180] (1)
设小车运动学中心的速度为[vωT]（其中[v]为运动学中心线速度，[ω]为运动学中心角速度），左右车轮的速度分别为[vl]、[vr]，小车的航向角为[θ]，小车沿世界坐标系的X轴正方向行驶时航向角[θ]被定义为0°。根据小车运动学中心与左右车辆的位置可以得出小车航向角随时间t变化的关系式，如式（2）所示。
[v=vl+vr2ω=vl-vr2y1θt=0tωdt] (2)
定义小车运动学中心在世界坐标系的位置变量为[X,Y,θT]，其中X、Y为小车运动学中心在世界坐标系中的坐标，由于世界坐标系的原点与小车的运动学中心坐标系在启动时完全重合，[θ]的变化直接反映了小车方向的改变，因此[θ]值可以作为小车运动学中心在世界坐标系中方向的表征。基于上述定义，[k]时刻与[k-1]时刻小车运动学中心在世界坐标系中具体位置分别为[Xk,Yk,θkT]，[Xk-1,Yk-1,θk-1T]。对[k-1]时刻小车运动的线速度和航向角进行分解，从而得出小车运动学中心[k-1]时刻到[k]时刻在世界坐标系中位置变量的关系矩阵，如式3所示。
[XkYkθk=Xk-1Yk-1θk-1+cosθk-1-sinθk-10sinθk-1cosθk-10001dXdYdθ]
(3)
根据式（3）可推导出[k]时刻的运动学增量[dXdYdθT]，如式（4）所示。
[dX=cosθk-1·vk-1dtdY=sinθk-1·vk-1dtdθ=ωk-1dt] (4)
其中：[vk-1]为k-1时刻小车运动学中心线速度，[ωk-1]为k-1时刻小车运动学中心线速度。
IMU使用四元数[a,b,c,e]表示其姿态数据，[k]时刻IMU的姿态数据为[qk=(ak,bk,ck,ek)]，[则k]时刻的增量[dq]如式（5）所示。
[dq=qk·q-1k-1] (5)
其中，[q-1k-1]是四元数[qk-1]的逆。
由于IMU和小车刚性连接，IMU的姿态变化可以直接映射到小车运动学中心，因此，如式（6）所示，小车运动学中心的姿态数据[qcar]与IMU的姿态数据[q]相同，小车运动学中心在[k]时刻的姿态增量[dqcar]与IMU在[k]时刻的姿态增量[dq]也相同。
[qcar=qdqcar=dq] (6)
由于车轮打滑和传感器噪声等因素的影响，位姿估计可能出现失真。为此，本文采用扩展卡尔曼滤波算法对车辆里程计与IMU数据进行融合，以提高位姿估计的准确性。扩展卡尔曼滤波算法分为预测和更新两个部分[12-13]。在此将k时刻小车运动学中心的完整状态变量定义为[xk=Xk,Yk,θk,qkT]。控制变量为[uk=vk,ωkT],小车的完整状态变量作为预测模型，车辆里程计与IMU数据作为观测模型。
系统的状态方程可由[k-1]时刻与[k]时刻小车运动学中心完整状态变量的关系表示：
[xk=[Xk-1+cosθk-1·vk-1dt,Yk-1+sinθk-1·vk-1dt,θk-1+ωk-1dt,qk-1?dqk-1]T]
(7)
其中：[dqk-1]是由IMU提供的姿态增量，[?]表示四元数的乘法。
系统的状态传递方程为：
[xk=fxk?1, uk+σk-1] (8)
其中：[fxk-1, uk]为[k-1]时刻系统的状态转移矩阵，[uk]为系统的控制输入，[σk-1]为系统在[k-1]时刻的高斯噪声。
系统的观测方程为[Zk]，其表达式如式（9）所示：
[Zk=hxk+gk] (9)
其中：[hxk]为观测到的模型传递矩阵，[gk]为观测到的高斯噪声。
扩展卡尔曼滤波算法包括预测和更新两部分。
（1）预测：计算的预测值为[xk=f(xk-1,uk-1)]，[k]时刻的协方差矩阵为[Pk=Fk-1Pk-1FTk-1+Qk-1]，其中[Fk-1]、 [Pk-1]和[Qk-1]分别为[k-1]时刻的状态传递雅可比矩阵、状态协方差矩阵和正态分布的方差。
（2）更新：在k时刻，卡尔曼增益[Kk]的计算、状态更新的表达式[xk]和协方差矩阵的更新矩阵[Pk]如式（10）所示，其中[Hk]为传感器函数[k]时刻的雅可比矩阵，[Rk]为[k]时刻正态分布的方差，I表示单位矩阵。
[Kk=PkHTkRk+HkPkHTk-1xk=xk-1+KkZk-hxk-1Pk=I-KkHkPk] (10)
以小车运动学中心为参考坐标系，二维激光雷达安装位置的坐标为[(x3,0,z3)]，与小车刚性连接。因此二维激光雷达产生的点云数据相对小车运动学中心具有固定的位移偏移量。点云数据在进行坐标转换的过程中需要考虑此偏移量，该位移偏移量可表示为[A=x3,0,z3,0T]。根据扩展卡尔曼滤波算法的推导，可以得出小车在运动的过程中，从小车坐标系中的点到世界坐标系的变换矩阵[B]。
[B=cosθk-sinθk0Xksinθkcosθk0Yk00100001] (11)
激光雷达在其自身的坐标系中采集的二维点云数据为[(0,ylidar,zlidar)]，转换为齐次坐标的表达式[Plidar=0,ylidar,zlidar,1T]，最终得到全局点云坐标转换的结果[Pworld]，表达式为：
[Pworld=BPlidar+A] (12)
2.2 数据处理
从采集到原始三维点云，到点云处理并识别形变得出结果的过程可分为点云处理、形变识别与分割和形变体积计算3个步骤，如图4所示。
<G:\武汉工程大学\2025\第3期\陈江-4.tif>[点云处理
点云滤波
点云降采样
点云配准][形变识别
与分割][形变体积计算]
图4 形变识别与体积计算算法流程图
Fig. 4 Flow chart of deformation recognition and volume calculation algorithm
2.2.1 点云处理 原始的三维点云并不能直接用作形变检测算法的输入源，因为激光雷达在实际采集数据的过程中会出现一些离群点与硬件噪声，因而需要对产生的噪声进行去除。首先是对点云进行滤波操作去除原始三维点云中的噪声和离群点，得到相对干净和均匀的点云数据。然后对采集到的点云数据进行降采样操作，寻找并提取点云图中关键点的同时降低整体的点云数量，有利于提高后期算法处理的速度。
上述操作完成后，两片点云还在其各自独立的坐标系下，因此要对两片点云进行配准。通过点云配准，使两片独立的点云能够实现最佳对齐，从而将其中一片点云准确地映射到另一片点云上，最终在同一坐标系下统一表示[14-16]。
本文采用迭代最近点算法（iterative closest point, ICP）作为配准算法。具体步骤为：（1）首先当两个点云初始位置差距较大时，对两片点云进行初始变换估计；（2）对于源点云中的每个点，查找其在目标点云中的最近点作为对应点对；（3）通过奇异值分解方法，最小化源点云和目标点云对应点对之间的欧氏距离，并计算出一个刚性变换（旋转矩阵和平移），使得源点云在变换后与目标点云更接近；（4）将计算得到的刚性变换应用到源点云上；（5）重复上述步骤直到满足某个收敛条件，如迭代次数达到上限，或者两次迭代之间的变换量小于某个阈值[17]。
2.2.2 形变识别与分割 对于已经配准的两期点云，以扫描时间靠前的一期点云P1为源点云，依次遍历扫描时间靠后的一期点云P2中的每一个点，设P2点云中当前正在被遍历的点为[Pi]，搜索P1点云中与[Pi]相距最小距离对应的点[Qi]和对应的最小距离[dmin]。理论上，当[Pi]与[Qi]对应现实世界同一位置或者相邻位置时，得到最小距离[dmin]，考虑到激光雷达硬件误差的存在，两期点云不可能精确对齐，因此要设置一个阈值距离用于判断点云是否发生明显的形变。
巷道点云数据在不同区域的局部点云密度不同，在形变识别的过程中如果使用固定阈值作为所有巷道点云是否发生形变的依据，容易导致误检或漏检。因此引入基于动态阈值的形变识别算法，根据局部环境和点云数据的变化，依次计算一片点云中每个点[Pi]的局部密度和噪声水平，动态调整当前点云中每个点的阈值[TD,i]，从而提高形变识别的准确性。具体包括以下步骤：
（1）计算点云中每个点[Pi]的局部密度。首先以点[Pi]为中心，定义半径为[rp]的邻域，然后，统计在该邻域内的点云数量，并计算[Pi]的局部密度[ρi]：
[ρi=j=1LIPi-Pj<rpVr] （13）
其中：L为当前点[Pi]所在点云的总点云数；[Pj]为当前点[Pi]所在点云中的另外一个点；[Vr]是半径为[rp]的邻域体积;[I(·)]为指示函数，当[Pi-Pj<rp]，即[Pi]与[Pj]的欧氏距离小于[rp]时，取值为1，反之为0。
（2）根据局部密度[ρi]和噪声水平，动态调整阈值[TD,i]：
[TD,i=T0?1-λ?ρi-ρρ] （14）
其中：[T0]为基准阈值，[λ]为调节系数，[ρ]为全局平均密度。
（3）在进行形变识别时，使用动态阈值[TD,i]进行判断：如果一片点云中的当前点[Pi]与另外一片点云的最小距离点[Qi]的距离关系满足[Pi-Qi<TD,i]，则说明此点属于未发生形变的点；反之，则说明此点云发生了形变，需要进行标记和点云数据的单独保存。
2.2.3 体积计算 形变识别与分割完毕后，使用切片法计算单独保存形变部分点云的体积：根据点云的具体形状，选择一个坐标轴作为切片方向，将点云沿该方向切分成N个厚度为Δh的点云切片截面。
对于每一个点云切片截面的点云数据，将其投影到与坐标轴垂直的二维平面上，形成一个由点云投影点围成的多边形；使用三角剖分法，将多边形分解为[M]个不重叠的三角形，假设剖分后的第i个三角形的顶点为[Ai(xi1,yi1)]、[Bi(xi2,yi2)]和[Ci(xi3,yi3)]，分别计算每个三角形的面积[St]，然后将所有三角形面积求和[18]，得到一个点云切片截面的面积[Ss,j]：
[St=12xi1yi2-yi3+xi2yi3-yi1+xi3yi1-yi2]
（15）
[Ss,j=i=1MSt,i] （16）
最后对每个点云切片截面的体积进行累加求和，得出形变部分整体的体积V：
[V=j=1NSs,j?Δh] （17）
3 试验数据分析
为了研究动态阈值形变识别算法在提升识别准确率方面的效果，选用真实矿山巷道的三维点云数据进行验证。通过控制变量法，逐个调试识别半径、基准阈值和调节系数，并观察识别准确率，以确定最佳的识别参数。然而在实际的参数调整过程中，每个参数的值需要缓慢调整，每次调试的结果提升并不显著。因此为了更加直观地展示动态阈值形变识别算法中参数的取值对点云形变识别的影响，本文选用了几组具有代表性的参数取值，并通过最终的误判点云数表示参数不同取值下算法的识别效果，具体结果如表1所示。所使用的巷道点云数据总点云数为5 947 017个，其中，误判点云数指在总点云中，巷道未发生形变但算法错误识别为形变的点云数量。误判点云占比的计算方法为：误判点云占比 = (误判点云数 / 巷道总点云数)×100%。
表1 动态阈值形变识别算法主要参数调试效果对比
Table 1 Comparison of debugging effects of main parameters in dynamic threshold deformation recognition algorithm
[序号 [rp] / m [T0] / m λ 误判点云数 /
个 误判点云占比 /
% 1 0.02 0.010 0.6 851 543 14.31 2 0.05 0.025 0.4 548 727 9.22 3 0.03 0.010 0.5 64 014 10.76 4 0.05 0.030 0.5 526 410 8.86 5 0.10 0.030 0.6 436 261 7.34 6 0.05 0.025 0.4 285 538 4.81 ]
通过表1可以看出，当识别半径[rp]= 0.05、基准阈值[T0] = 0.025、调节系数[λ] = 0.4时，动态阈值形变识别算法取得最佳效果，误判点云占比仅为4.81%。
为了验证系统是否拥有良好地对矿山巷道进行三维扫描、识别形变的能力，同时对比系统计算值与实际值的误差，进行了两次试验进行验证。一次模拟实验选用某校内某处长走廊模拟矿山巷道，使用沙堆模拟矿山巷道内发生的形变与脱落；另外一次现场试验的地点选在湖北省宜昌市兴山县树崆坪磷矿1172主平硐一号分段巷道内，作为一个已经开采完毕并且经过长时间沉降的矿山巷道，短期内发生形变的可能性很小，因此使用外加形变物体模拟巷道内部发生的形变。
具体的实验过程为：采集未发生形变的巷道作为一期点云数据；在模拟实验现场的不同位置放置模拟发生形变的物体，进行二期点云数据的采集。为了便于对比系统计算值与物体实际数值，模拟实验使用了体积相同但形状不同的沙堆模拟更加贴近真实的形变与脱落；现场试验则使用了纸团、圆锥、三棱锥和四棱锥模拟发生形变的部分，如图5、图6所示。
<G:\武汉工程大学\2025\第3期\陈江-5.tif>[（a）][（b）][（c）]
图5 模拟形变沙堆:(a)锥形沙堆, (b)异形沙堆1,
(c)异形沙堆2
Fig. 5 Simulated deformed sand pile:(a) conical sand pile, (b) special-shaped sand pile 1, (c) special-shaped sand pile 2
<G:\武汉工程大学\2025\第3期\陈江-6.tif>[（a）][（b）]
图6 模拟形变物体:(a)纸团, (b)圆锥、三棱锥、四棱锥
Fig. 6 Simulate deformable objects:(a)paper balls,
(b)cone, triangular pyramid, four-sided pyramid
数据处理时，将模拟实验与现场试验一期、二期的点云数据进行点云滤波和降采样处理，然后使用ICP算法对两期点云进行配准，现场试验的点云数据进行ICP配准前后整体点云的效果与部分点云的细节对比如图7所示。
<G:\武汉工程大学\2025\第3期\陈江-7.tif>[（a）][（b）][（c）][（d）][配准前][配准前][配准后][配准后]
图7 配准前后两期点云对比
Fig. 7 Comparison of point clouds before and after
registration
将两片点云配准以后，需要对二期点云以一期点云为参考进行形变识别与分割，具体操作为：遍历二期点云中的每个点，并在遍历过程中针对当前点，找到其在一期点云中的最近邻点。如果这两个点之间的距离超过当前遍历点的动态阈值[TD,i]，则将其进行标记和单独保存，形变部分点云在整体巷道点云中的标记效果如图8所示。
<G:\武汉工程大学\2025\第3期\陈江-8.tif>[（a）][（b）]
图8 动态阈值形变识别算法识别效果图:
(a) 形变纸团, (b) 形变椎体
Fig. 8 Recognition effect of the dynamic-threshold
deformation recognition algorithm： (a) paper balls,
(b) vertebral bodies
模拟实验和现场试验过程中识别到的形变部分点云，经过单独分割与保存后，其形状如图9所示。
<G:\武汉工程大学\2025\第3期\陈江-9.tif>[（a）][（b）][（c）][（d）][（e）][（f）][锥形沙堆][异形沙堆1][异形沙堆2][四棱锥][三棱锥][纸团]
图9 形变物体识别效果图
Fig. 9 Recognition effect of deformable objects
最后，在系统中测量形变物体的厚度，并利用切片法计算形变部分点云的体积，结果如表2、表3所示。其中误差率的计算公式为：误差率=|测量值-实际值|/实际值×100%。
从表2和表3可以看出，形变物体厚度的误差率平均为3%，形变物体体积的系统计算值与实际值的误差率平均为12%。
表2 形变物体厚度误差
Table 2 Thickness errors of deformed objects
[名称 形变物体厚度 / mm 误差率 / % 雷达测量值 实际值 纸团 23.8 23.2 2.8 三棱锥 51.7 50.1 3.2 四棱锥 48.9 50.4 2.9 锥形沙堆 85.7 83.0 3.2 异形沙堆1 51.3 49.8 3.0 异形沙堆2 48.5 47.2 2.8 ]
表3 形变物体体积误差
Table 3 Volume errors of deformed objects
[名称 形变物体体积 / mm3 误差率 / % 系统计算值 实际值 三棱锥 19 756.43 18 042.17 9.5 四棱锥 37 333.36 41 666.67 10.3 锥形沙堆 568 470.78 500 000.00 13.6 异形沙堆1 568 186.47 500 000.00 12.4 异形沙堆2 565 547.57 500 000.00 13.1 ]
4 结 论
本文设计并实现了一种基于激光雷达的巷道形变检测系统，介绍了系统的硬件组成和软件设计，实现了数据采集、数据融合、点云数据处理与形变检测等功能。在形变检测部分提出了动态阈值形变识别算法，将误判点云占比降至4.81%。现场实验结果表明，在测量形变物体厚度时误差的平均值为3%；通过对形变部分进行量化分析，在计算形变体积方面误差的平均值为12%。此装置能够较为精准地识别出巷道内形变发生的位置和形变部分的体积大小，满足实际巷道形变测量的需求。