一种基于激光雷达点云的矿山开采量计算方法

摘要

本发明属于激光雷达点云处理和矿山监测领域，具体技术方案为：一种基于激光雷达点云的矿山开采量计算方法，首先，基于激光雷达系统对矿山进行定期三维地形测绘，获取矿山原始点云数据；然后对激光雷达点云数据进行去噪、配准等预处理；其次检测点云局部空洞面并利用邻域信息内插三维数据点进行修补，前后两期点云数据形成某时段不规则疏密多变矿体的三维模型；最后，设置离散化步长对矿体三维模型进行块体离散化并提取块体体积，各单元块体体积累加得到整个矿山变化体积，即该时段的矿山开采量，本发明能够对矿山开挖方量进行快速估算，该方法具备测量高效安全、灵活便捷、即时性高等特点，可应用于三维变化监测、工程测管、智慧矿山建设等领域。

权利要求

1.一种基于激光雷达点云的矿山开采量计算方法，其特征在于，具体步骤如下：一、激光雷达点云获取：对矿山进行矿区勘探后确定监测范围，根据实际监测区域大小与高差选取激光雷达系统并制定数据采集方案，完成外业数据采集，经计算机软件处理后解算得到激光雷达点云数据；

二、点云预处理：对原始激光雷达点云数据进行去噪、配准预处理；

三、点云空洞面修补，包括以下两个步骤：

S31、空洞面检测：在三维空间中，基于原始三维点云数据构建初始三角网络，遍历计算所有三角网格的面积，设定最小空洞面积阈值P，自动检测面积大于阈值P的网格为需要修补的空洞面；

S32、空洞面修补：

利用邻域信息对空洞网格进行三维数据点内插，对需要修补的三角网格的每条边取中点，分别为点 ，连接三点后将空洞网格分割为四个新的三角形，分别取四个三角形的重心作为新的离散点内插至空洞网格中；

再次计算四个新的三角形的面积，仍大于阈值P的三角网格被划分为空洞网格，重复进行三角形分割与离散点内插，直到整体三角网格中的每一个网格面积均小于P，停止内插，得到精细矿山三维点云数据；

空洞面修补后，前后两期点云数据形成不规则疏密多变矿体的三维模型；

四、开采量计算：根据矿体三维模型的面积大小与点密度确定离散化步长，对其进行单元块体离散化，获取各块体的网格面积与块体相对高度，计算各单元块体的体积，将整个矿体所有单元块体的体积叠加获取整个矿山变化体积。

2.根据权利要求1所述的一种基于激光雷达点云的矿山开采量计算方法，其特征在于，在步骤一中，包括以下两个步骤：S11、激光雷达点云数据采集：对矿山的地形起伏、开采边界、典型地物进行勘探后确定矿山监测区域，根据监测区域的大小和高差选取续航时间与探测距离符合要求的激光雷达系统，结合系统具体的性能参数确定测量距离、扫描速度与扫描角度，对矿山实施外业数据采集；

S12、数据解算：基于原始惯导数据、原始激光测距数据、GPS数据通过计算机联合解算得到激光雷达点云数据。

3.根据权利要求2所述的一种基于激光雷达点云的矿山开采量计算方法，其特征在于，所述原始激光测距数据由激光扫描仪扫描采集得到。

4.根据权利要求1所述的一种基于激光雷达点云的矿山开采量计算方法，其特征在于，在步骤二中，包括以下两个步骤：S21、点云去噪：对原始激光雷达点云采用稀疏离群点移除方法进行去噪，该方法计算每个点到其所有临近点的平均距离，假设得到的结果是一个高斯分布，平均距离在标准范围之外的点，被定义为离群点并从数据集中去除掉；

S22、点云配准：对开采前后两期点云数据采用迭代最近点自动配准算法，ICP算法通过找到两组点云集合中距离最近的点对，根据估计的变换关系来计算距离最近点对经过变换之后的误差，不断的迭代直至使得设定的目标函数达到最小值，得到最优平移矩阵和旋转矩阵，以确定最终的变换关系，目标函数如下：式中，R表示旋转矩阵，T表示平移矩阵，k为待配准点云的点数，ɑi为参考点云；bi为待配准点云；

最终将待配准的点云与参考点云的数据对进行一一对应，即完成两期点云数据配准。

5.根据权利要求1所述的一种基于激光雷达点云的矿山开采量计算方法，其特征在于，在步骤S31中，初始三角网络构建的步骤为：

1）、以点云数据左下角任一点为起始点；

2）、找出与起始点最近的数据点相互连接形成三角形的一条边作为基线；

3）、在基线右侧进行扩展，搜索与基线距离最近的点构成起始三角形的第三点；

4）、基线的两个端点与第三点相连，成为新的基线；

5）、继续搜索与新的基线距离最近的第三点构成新的三角形，直至所有的点都被包含在三角网格中，停止搜索，初始三角网络构建完毕。

6.根据权利要求1所述的一种基于激光雷达点云的矿山开采量计算方法，其特征在于，在步骤四中，还包括以下两个步骤：S41、单元块体离散化：结合点云密度与测区大小确定离散化步长，将不规则疏密多变矿体三维模型划分成n个单元块体；

S42、矿山开采量计算：计算各单元块体的顶面积St和底面积Sd，取顶面积与底面积的平均值作为网格面积Sm，计算公式如下：其中，i为各单元块体的序号， ；

三维模型各单元块体对应的顶面与底面点云的高程值的平均高度作为单元块体高程值的最大值Zmax与最小值Zmin，对每个单元块体对应的顶面高程值与底面高程值做差得到块体相对高度H，计算公式如下：每个单元块体利用网格面积与块体相对高度相乘得到离散化块体的体积，经计算并累加求和得到矿体变化体积ΔV，计算公式如下：。

说明书

一种基于激光雷达点云的矿山开采量计算方法

技术领域

[0001] 本发明属于激光雷达点云处理和矿山监测领域，具体涉及一种基于激光雷达点云的矿山开采量计算方法。

背景技术

[0002] 矿产资源是人类赖以生存和经济发展的重要物资基础，开采量动态监测可为其科学高效的生产管理提供依据，对矿产资源的可持续发展与智慧矿山建设意义重大。矿山开采尤其是露天开采覆盖范围大且地形复杂多变，导致开采量测算难度大、效率低、成本高。如何快速准确地获取矿山地形变化并计算矿山开采量成为了一个亟待解决的问题。

[0003] 已有矿山开采量计算方法主要包括两种，一种是基于人工测绘的开采量计算，此类方法通过地测人员手动测量获取的开采区域工作面三维坐标生成地形图，勘测效率低且易受人为因素影响，导致开采量计算速度与精度不稳定；另一种是基于无人机可见光影像的开采量计算，此类方法通过获取时序无人机影像后重建三维模型，采集效率高于人工测绘方式，但由于露天矿地形起伏较大且边界持续移动，三维模型高程的真实性难以得到保障，依然无法在精度上满足矿山开采量计算的工程需求。

[0004] 激光雷达系统是一种集激光，全球定位系统(GPS)和惯性导航系统(INS)三种技术于一身的系统，可高效、精确、无损地获得地表三维信息。将激光雷达系统引入矿山开采量动态监测中，可通过定期获取矿山的精细三维点云数据，构建某时期不规则疏密多变矿体的三维模型，快速准确获取矿山开挖过程中的真实地形地貌变化，为开采量计算提供可靠的数据支撑。

发明内容

[0005] 本发明为了满足矿山开采量计算快速准确的需求，提供一种基于激光雷达点云的矿山开采量计算方法，通过对矿山进行定期三维地形测绘更新并基于精细三维点云数据计算矿山变化体积，从而实现矿山开采过程中的动态开采量计算。

[0006] 为实现上述目的，本发明所采用的技术方案为：一种基于激光雷达点云的矿山开采量计算方法，具体步骤如下：

[0007] 一、激光雷达点云获取，对矿山进行矿区勘探后确定监测范围，根据实际监测区域大小与高差选取激光雷达系统并制定数据采集方案，完成外业数据采集，经计算机软件处理后解算得到激光雷达点云数据。

[0008] 进一步是，步骤一包括以下步骤：

[0009] S11、激光雷达点云数据采集：对矿山的地形起伏、开采边界、典型地物进行勘探后确定矿山监测区域，根据监测区域的大小和高差选取续航时间与探测距离符合要求的激光雷达系统，结合系统具体的性能参数确定测量距离、扫描速度与扫描角度，对矿山实施外业数据采集。

[0010] S12、数据解算：基于原始惯导数据、原始激光测距数据（由激光扫描仪扫描采集得到）、GPS数据通过计算机联合解算得到激光雷达点云数据。

[0011] 二、点云预处理：对原始激光雷达点云数据进行去噪、配准等预处理。

[0012] 进一步是，步骤二包括以下步骤：

[0013] S21、点云去噪：对原始激光雷达点云采用稀疏离群点移除（Statistical Outlier Removal，SOR）方法进行去噪，该方法计算每个点到其所有临近点的平均距离，假设得到的结果是一个高斯分布，平均距离在标准范围（由全局距离平均值和方差定义）之外的点，被定义为离群点并从数据集中去除掉。

[0014] S22、点云配准：对开采前后两期点云数据采用迭代最近点（Iterative Closest Point, ICP）自动配准算法，ICP算法通过找到两组点云集合中距离最近的点对，根据估计的变换关系（旋转与平移）来计算距离最近点对经过变换之后的误差，不断的迭代直至使得设定的目标函数达到最小值，得到最优平移矩阵和旋转矩阵，以确定最终的变换关系，目标函数如下：

[0015]

[0016] 式中，R表示旋转矩阵，T表示平移矩阵，k为待配准点云的点数, ɑi与bi分别对应参考点云与待配准点云。

[0017] 最终将待配准的点云与参考点云的数据对进行一一对应，即完成两期点云数据配准。

[0018] 三、点云空洞面修补：在三维空间中根据点云数据构建初始三角网格，设定最小空洞面积阈值P，自动检测面积大于阈值P的点云空洞网格。利用邻域信息对空洞网格进行三维数据点内插，点云空洞面修补完成后得到矿山精细三维点云数据，前后两期点云数据形成不规则疏密多变矿体的三维模型。

[0019] S31、空洞面检测：矿山由于其本身反射率有限且存在部分水域，数据采集中往往存在三维点云局部空洞的问题，尤其是在露天矿的矿底部与侧面，开采量计算有赖于完整精细的点云数据，需要对空洞面进行检测与修补。基于原始三维点云数据构建初始三角网络，遍历计算所有三角网格的面积，设定最小空洞面积阈值P，自动检测面积大于阈值P的网格为需要修补的空洞面。

[0020] 其中，初始三角网络构建的基本步骤是：

[0021] 1）、以点云数据左下角任一点为起始点；

[0022] 2）、找出与起始点最近的数据点相互连接形成三角形的一条边作为基线；

[0023] 3）、在基线右侧进行扩展，搜索与基线距离最近的点构成起始三角形的第三点；

[0024] 4）、基线的两个端点与第三点相连，成为新的基线；

[0025] 5）、继续搜索与新的基线距离最近的第三点构成新的三角形，直至所有的点都被包含在三角网格中，停止搜索，初始三角网络构建完毕。

[0026] S32、空洞面修补：对需要修补的三角网格的每条边取中点，分别为点，连接三点后将空洞网格分割为四个新的三角形，分别取四个三角形的重心作为新的离散点内插至空洞网格中，重心坐标为
。

[0027] 再次计算四个新的三角形的面积，仍大于阈值P的三角网格被划分为空洞网格，重复进行三角形分割与离散点内插，直到整体三角网格中的每一个网格面积均小于P，停止内插，得到精细矿山三维点云数据。

[0028] 空洞面修补后，前后两期点云数据形成不规则疏密多变矿体的三维模型。

[0029] 四、开采量计算：根据矿体三维模型的面积大小与点密度确定离散化步长，对其进行单元块体离散化，获取各块体的网格面积与块体相对高度，计算各单元块体体积，将整个矿体所有单元块体的体积叠加获取整个矿山变化体积，即该时段的矿山开采量。

[0030] 进一步是，步骤四包括以下步骤：

[0031] S41、单元块体离散化：结合点云密度与测区大小确定离散化步长，将不规则疏密多变矿体三维模型划分成个单元块体。

[0032] S42、矿山开采量计算：计算各单元块体的顶面积St和底面积Sd，取顶面积与底面积的平均值作为网格面积Sm，计算公式如下：

[0033]

[0034] 其中，i为各单元块体的序号， ；

[0035] 三维模型各单元块体对应的顶面与底面点云的高程值的平均高度作为单元块体高程值的最大值Zmax与最小值Zmin，对每个单元块体对应的顶面高程值与底面高程值做差得到块体相对高度H，计算公式如下：

[0036]

[0037] 每个单元块体利用网格面积与块体相对高度相乘得到离散化块体的体积，经计算并累加求和得到矿体变化体积，即为该时段矿山的开采量，计算公式如下：

[0038] 。

[0039] 与现有技术相比，本发明的有益效果是：

[0040] 一、本发明将激光雷达探测手段引入到矿山开采量动态监测中，可以及时更新开采过程中的地形变化，基于激光雷达点云提出了一种快速准确的矿山开采量计算的方法。

[0041] 二、本发明通过对矿山点云数据进行去噪、配准与局部空洞面修补，保障了三维点云的质量，可形成完整的不规则疏密多变矿体三维模型，避免了由于点云缺失可能造成的算法失效，为矿山开采量计算提供可靠的数据支撑。

[0042] 三、本发明对不规则疏密多变矿体的三维模型进行单元块体离散化，通过确定最优步长可精确提取各单元块体的体积变化量，从而实现快速准确的矿山开采量计算，为矿山工程测管及产量监控提供依据。

附图说明

[0043] 图1为本发明基于激光雷达点云的矿山开采量计算方法流程图。

[0044] 图2为第一期露天煤矿矿山三维点云实景模型图。

[0045] 图3为第二期露天煤矿矿山三维点云实景模型图。

[0046] 图4为本发明得到的矿体相对变化高度图。

[0047] 图5为选取的矿山监测区域标注图。

[0048] 图6为第一期露天煤矿矿山三维点云数据图。

[0049] 图7为第二期露天煤矿矿山三维点云仿真数据图。

[0050] 图8为第一期、第二期露天煤矿矿山三维点云数据的数据剖面图高度对比图。

具体实施方式

[0051] 为了使本发明所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

[0052] 如图1所示，一种基于激光雷达点云的矿山开采量计算方法，采用无人机激光雷达系统对某露天煤矿进行两期矿山点云数据采集，具体计算方法主要包括以下步骤：

[0053] 一、激光雷达点云获取：对矿山进行矿区勘探后确定监测范围，根据实际监测区域选取激光雷达系统并制定数据采集方案，完成外业数据采集，经计算机软件处理后解算得到激光雷达点云数据。

[0054] 进一步是，步骤一包括以下步骤：

[0055] S11、激光雷达点云数据采集：对矿山的地形起伏、开采边界、典型地物进行勘探后确定矿山监测区域，根据监测区域的大小和高差选取续航时间与探测距离符合要求的激光2
雷达系统。此实施例中的监测区域约0.4km ，最大高差约85m，选用无人机激光雷达系统进行数据采集。结合该系统具体的性能参数确定测量距离、扫描速度与扫描角度，对矿山实施外业数据采集。此处设置飞行速度为5m/s，扫描重叠度为40%，选取90m的相对高度进行仿地飞行。

[0056] S12、数据解算：基于原始惯导数据、RTK基站数据、原始激光测距数据（由激光扫描仪扫描采集得到）、机载GPS数据通过计算机联合解算得到矿山激光雷达点云数据，此处数据解算选用CGCS2000大地参考坐标系完成点云数据解算。以上为数据解算的常规流程，均为本领域的技术人员熟知的。

[0057] 二、点云预处理：对原始激光雷达点云数据进行去噪、配准等预处理。

[0058] 进一步是，步骤二包括以下步骤：

[0059] S21、点云去噪：对原始矿山激光雷达点云采用稀疏离群点移除（Statistical Outlier Removal，SOR）方法进行去噪，该方法计算每个点到其所有临近点的平均距离，假设得到的结果是一个高斯分布，平均距离在标准范围（由全局距离平均值和方差定义）之外的点，被定义为离群点并从数据集中去除掉。本实施例中指定邻域点数量为20，以此计算平均距离，基于点云的平均距离的方差选取2。

[0060] S22、点云配准：对开采前后两期点云数据采用迭代最近点（Iterative Closest Point, ICP）自动配准算法，ICP算法通过找到两组点云集合中距离最近的点对，根据估计的变换关系（旋转与平移）来计算距离最近点对经过变换之后的误差，不断的迭代直至使得设定的目标函数达到最小值，得到最优平移矩阵和旋转矩阵，以确定最终的变换关系，目标函数如下：

[0061]

[0062] 式中，R表示旋转矩阵，T表示平移矩阵，k为待配准点云的点数，ɑi为参考点云；bi为待配准点云。

[0063] 最终将待配准的点云与参考点云的数据对进行一一对应，即完成两期点云数据配准。

[0064] 三、点云空洞面修补：在三维空间中根据点云数据构建初始三角网格，设定最小空洞面积阈值P，自动检测面积大于阈值P的点云空洞网格。利用邻域信息对空洞网格进行三维数据点内插，点云空洞面修补完成后得到矿山精细三维点云数据，前后两期点云数据形成不规则疏密多变矿体的三维模型。

[0065] S31、空洞面检测：矿山由于其本身反射率有限且存在部分水域，数据采集中往往存在三维点云局部空洞的问题，尤其是在露天矿的矿底部与侧面，开采量计算有赖于完整精细的点云数据，需要对空洞面进行检测与修补，基于原始三维点云数据构建初始三角网络。

[0066] 其中，初始三角网络构建的基本步骤是：

[0067] 1）、以点云数据左下角任一点为起始点；

[0068] 2）、找出与起始点最近的数据点相互连接形成三角形的一条边作为基线；

[0069] 3）、在基线右侧进行扩展，搜索与基线距离最近的点构成起始三角形的第三点；

[0070] 4）、基线的两个端点与第三点相连，成为新的基线；

[0071] 5）、继续搜索与新的基线距离最近的第三点构成新的三角形，直至所有的点都被包含在三角网格中，停止搜索，初始三角网络构建完毕。通过遍历计算所有三角网格的面积，设定最小空洞面积阈值P，自动检测面积大于阈值P的网格为需要修补的空洞面，此处P2
取0.0025m。

[0072] S32、空洞面修补：对需要修补的三角网格的每条边取中点，分别为点，连接三点后将空洞网格分割为四个新的三角形，分别取四个三角形的重心作为新的离散点内插至空洞网格中，重心坐标为
。

[0073] 再次，计算四个新的三角形的面积，面积仍大于0.0025 m2的三角网格被划分为空洞网格，重复进行三角形分割与离散点内插，直到整体三角网格中的每一个网格面积均小2
于0.0025 m，停止内插，得到精细矿山三维点云数据。

[0074] 空洞面修补后，前后两期点云数据形成矿体三维模型。

[0075] 四、开采量计算：根据矿体三维模型的面积大小与点密度确定离散化步长，对其进行单元块体离散化，获取各块体的网格面积与块体相对高度，计算各单元块体体积，将整个矿体所有单元块体的体积叠加获取整个矿山变化体积，即该时段的矿山开采量。

[0076] 进一步是，步骤四包括以下步骤：

[0077] S41、单元块体离散化：结合点云密度与测区大小确定离散化步长，矿山开采量监测中一般选取0.4m作为离散化步长，将某时段不规则疏密多变矿体的三维模型划分成n个单元块体。

[0078] S42、矿山开采量计算，计算各单元块体的顶面积St和底面积Sd，取顶面积与底面积的平均值作为网格面积Sm，计算公式如下：

[0079]

[0080] 其中，i为各单元块体的序号， ；

[0081] 矿体三维模型各单元块体对应的顶面与底面点云的高程值的平均高度作为单元块体高程值的最大值Zmax与最小值Zmin，对每个单元块体对应的顶面高程值与底面高程值做差得到块体相对高度H，计算公式如下：

[0082]

[0083] 每个单元块体利用网格面积与块体相对高度相乘得到离散化块体的体积，经计算并累加求和得到矿体变化体积ΔV，即为该时段矿山的开采量，计算公式如下：

[0084]

[0085] 具体实施时，所采用的无人机激光雷达系统的参数如下表：

[0086]

[0087] 利用该无人机激光雷达系统对露天煤矿矿山进行两期数据航摄，两期数据间隔3个月。如图2和图3给出了本实施例经预处理及局部点云修补后的两期露天煤矿矿山三维点2
云数据，该区域的两期点云数据均包含约1600万个三维坐标数据，点密度约为190pts/m 。
基于前后两期点云得到矿体三维模型，并进行高度变化值提取，可得到图4所示矿体的相对变化高度图。运用C++编程实现本发明基于激光雷达的矿山开采量计算方法，并按照不同的网格步长对矿体三维模型进行块体离散化后提取该时期的开采量。其计算平台为CPU：
Inter(R) Core(TM) i7‑9750H 2.60GHz，RAM 16G的Windows10计算机，对应实际矿山开采
3
量为335000m，选取不同步长得到的开采区变化体积计算结果如下表所示：

[0088]

[0089] 由上表的计算结果可知，步长选取较小时，单元匹配比例较低，无法正确计算开采量；步长选择较大时，尽管两期模型间各单元块体匹配度高，但选取的块体相对高度误差较大，开采量估算结果与实际开采量相比偏大。基于激光雷达点云对矿山开采量进行计算时，2
当矿山点云密度在100至200 pts/m的范围时步长设置为0.4m至0.6m，一般取0.4m即可，误差在5%以内。从运行时间方面来看，本方法计算开采量的运行速度较快，可用于工程量快速核算。

[0090] 基于某矿山于某时间航摄得到的原始三维点云数据的矿坑坑底及向上50m形成一个封闭区域，选择此区域作为矿山监测区域，如图5所示。将监测区域的点云高程通过仿真整体降低3m，得到一组变化后的模拟点云数据。变化前后两期点云数据如图6、图7所示，两3
期数据剖面图高度对比如图8所示。经计算，该矿山监测区域的开采量为1042239.8 m。仿
2 3
真该监测区域面积为0.3474km，变化体积为1042200 m，误差基本为0%。

[0091] 本发明将激光雷达探测手段引入到矿山开采量动态监测中，对不规则疏密多变矿体的三维模型进行单元块体离散化，通过确定最优步长可精确提取各单元块体的体积变化量，从而实现快速准确的矿山开采量计算，可以及时更新开采过程中的地形变化，为矿山工程测管及产量监控提供依据。

[0092] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包在本发明范围内。