CN112945196B - 一种基于点云数据的露天矿台阶线提取和边坡监测的方法 - Google Patents

一种基于点云数据的露天矿台阶线提取和边坡监测的方法 Download PDF

Info

Publication number
CN112945196B
CN112945196B CN202110109374.7A CN202110109374A CN112945196B CN 112945196 B CN112945196 B CN 112945196B CN 202110109374 A CN202110109374 A CN 202110109374A CN 112945196 B CN112945196 B CN 112945196B
Authority
CN
China
Prior art keywords
point
points
curvature
point cloud
fitting
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110109374.7A
Other languages
English (en)
Other versions
CN112945196A (zh
Inventor
王植
安世缘
邹俊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northeastern University China
Original Assignee
Northeastern University China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northeastern University China filed Critical Northeastern University China
Priority to CN202110109374.7A priority Critical patent/CN112945196B/zh
Publication of CN112945196A publication Critical patent/CN112945196A/zh
Application granted granted Critical
Publication of CN112945196B publication Critical patent/CN112945196B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C9/00—Measuring inclination, e.g. by clinometers, by levels
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B21/00—Measuring arrangements or details thereof, where the measuring technique is not covered by the other groups of this subclass, unspecified or not relevant
    • G01B21/02—Measuring arrangements or details thereof, where the measuring technique is not covered by the other groups of this subclass, unspecified or not relevant for measuring length, width, or thickness
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C5/00—Measuring height; Measuring distances transverse to line of sight; Levelling between separated points; Surveyors' levels

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Image Processing (AREA)

Abstract

本发明公开一种基于点云数据的露天矿台阶线提取和边坡监测的方法，该方法通过渐进形态学滤波来消除挡土墙等人工地物对露天矿台阶线提取的干扰，提高了精度；根据台阶线特征点的几何属性在提取特征点前先计算每点曲率，将曲率大于阈值的点作为潜在特征点，再对潜在特征点进行细化，大大减少了计算复杂度，提高了整个运行速度；将曲率指数结合AGPN算法提高了提取台阶线上的特征点的精度，本发明的方法从整体上提高了露天矿台阶线提取的自动化程度。此外通过对两期台阶线提取结果分析可以得出露天矿边坡坡度变化和边坡位移量，这对露天矿的安全生产具有重要意义。

Description

一种基于点云数据的露天矿台阶线提取和边坡监测的方法
技术领域
本发明涉及边坡数据采集与监测技术领域，尤其涉及一种基于点云数据的露天矿台阶线提取和边坡监测的方法。
背景技术
露天矿在进行开采工作时，为了便于从上至下逐层开采，要把矿区内的矿物、岩石划分为具有一定厚度的水平分层，这些分层的工作面就是台阶。露天矿台阶线对露天矿的勘测、采矿设计、计划编制、区块设计、爆破设计、矿区划分等业务均具有重要的作用。
现阶段在实际生产过程中提取露天矿上台阶线等特征线使用的主要方法是根据DSM进行人工绘制，这种方法工作量大、效率低。利用轻便的无人机，获取露天矿序列影像，进行三维建模，检测采场台阶线，自动精确解算剥离量，对加强矿山生产监控和成本管理具有重要意义。史红霞等在中国机械工程中发表的《基于法向量区域聚类分割的点云特征线提取》中提出采用自适应邻域的主成分分析法(PCA)估算模型法向，利用萤火虫算法(FA)优化的模糊C均值(FCM)算法对法向进行聚类实现模型的有效分割。构造点集剔除与合并准则从各分割块边界点集中析取候选特征点，再以局部邻域主轴方向为基准提取特征点的方法。柏宏强等在激光与光电子学进展中发表的《基于三维激光点云特征线提取的溶洞多分辨率三维重建方法研究》中提出采用改进邻近点几何特征提取溶洞特征值，增加法向量角作为检测特征点的参数；用社会粒子群(SPSO)算法与模糊C-均值(FCM)聚类算法实现点云分类；采用折线生长方法将特征点连接成特征线的方法。上述方法中，存在边缘和特征线检测不全，受噪声点干扰大，且当场景复杂时，在非平面区域会产生现实中不存在的边缘或特征线等问题。
发明内容
针对上述现有技术的不足，本发明提供一种基于点云数据的露天矿台阶线提取和边坡监测的方法。
为解决上述技术问题，本发明所采取的技术方案是：一种基于点云数据的露天矿台阶线提取方法，包括如下步骤：
步骤1：对露天矿点云数据使用渐进形态学滤波算法进行预处理，通过逐渐增加滤波器的窗口大小和不断增大的高程差阈值移除点云数据中的地物数据，提取地面点，过程如下：
步骤1.1：渐进形态学滤波算法在进行第一次迭代时，将最小高程表面与初始滤波器的窗口大小作为输入，对三维激光点云数据作滤波运算，得到已滤波表面；
步骤1.2：比较已滤波表面在滤波前后的高程差是否小于高程差阈值，若小于则将窗口大小赋值给标记数组，来增加窗口大小；
所述高程差阈值根据研究区域的地形坡度确定，方法如下：
假设坡度恒定，则地形的最大高程差，窗口大小和地形坡度之间存在下列公式：
Figure BDA0002918689120000021
Figure BDA0002918689120000022
其中，dh0表示初始高程差阈值；s表示斜率；c表示格网大小；dhmax表示最大高程差阈值；wk表示运行到第k次窗口大小，dhT，K表示高程阈值，dhmax(t)，K表示地形的最大高程差。
步骤1.3：将上一迭代中获得的已滤波表面和增加的窗口大小作为滤波器的输入，继续进行滤波运算迭代，直到窗口大小为最大窗口。
步骤2：基于顾及邻域几何属性的三维边缘检测与曲率指数加权方法提取台阶线特征点，过程如下：
步骤2.1：点P为点云中任意一点，对点云中点P及其邻域内的点使用最小二乘法拟合曲面，并对其进行参数化：r(xy)＝(x，y，ax2+bxy+cy2)，其中，x为点P的横坐标，y为点P的纵坐标，a、b、c为对曲面参数化后的参数；
步骤2.2：设r(x，y)的偏导数分别为rx，ry，rxx，rxy，ryy，然后根据以下公式计算平均曲率H和高斯曲率K：
Figure BDA0002918689120000023
Figure BDA0002918689120000024
其中，
Figure BDA0002918689120000025
为曲面的单位法向量；E＝rx·rx＝1+(2ax+by)2；
F＝rx·ry＝(2αx+by)(2cy+bx)，G＝ry·ry＝1+(2cy+bx)2，
Figure BDA0002918689120000026
步骤2.3：最后通过下列公式计算主曲率：
Figure BDA0002918689120000027
Figure BDA0002918689120000028
其中，kmax为最大主曲率，kmin为最小主曲率；
步骤2.4：计算点云中点P的曲率指数：
Figure BDA0002918689120000031
步骤2.5：将曲率指数大于曲率指数阈值的点定义为潜在特征点，对潜在特征点使用顾及邻域几何属性的三维边缘检测算法AGPN细化特征点，将潜在特征点中的内点集进行法向量优化；将潜在特征点中的外点集直接定义为非特征点，过程如下：
步骤2.5.1：对曲率指数大于曲率指数阈值的点其中任意一点P使用k-d tree算法检索其邻域点集S；
步骤2.5.2：对邻域点集S使用RANSAC算法拟合平面，将点集S分为内点集和外点集；
步骤2.5.3：若点P在内点集中，则计算拟合平面的法向量，将其作为点P的法向量进行法向量优化，将潜在特征点中的外点集直接定义为非特征点。
步骤2.6：根据优化后的法向量和拟合平面上一对相互垂直的向量建立空间直角坐标系，设空间向量
Figure BDA0002918689120000032
在x轴上，向量
Figure BDA0002918689120000033
在y轴上，优化后的法向量为
Figure BDA0002918689120000034
计算角度阶跃Gθ，公式如下：
Figure BDA0002918689120000035
Gθ＝max(θi+1-θi)(i＝1…N-1) (9)
式中，空间向量
Figure BDA0002918689120000036
是由当前点P和其邻域内拟合平面提取的内点构成的，θi为向量
Figure BDA0002918689120000037
与向量
Figure BDA0002918689120000038
夹角；
步骤2.7：将角度阶跃大于阈值的点作为提取的台阶线特征点。
步骤3：采用移动最小二乘法拟合台阶线，得到提取后的露天矿台阶线，过程如下：
步骤3.1：对提取出的台阶线特征点进行曲线拟合，在要进行曲线拟合的部分中，其目标子区域在拟合过程中使用的函数记作f(x)：
Figure BDA0002918689120000039
α(x)＝[a1(x)，ax(x)，……，am(x)]T (11)
pT(x)＝[p1(x)，p2(x)，......，pm(x)] (12)
其中，α(x)代表的为待定的系数矢量；pT(x)代表的为m次多项式空间的一组基，ai(x)为待定的系数矢量，pi(x)为基函数；
步骤3.2：若计算点x的邻域内有肝个节点，近似函数f(x)在这些节点处的误差的加权平方和如下：
Figure BDA0002918689120000041
其中，yi是x＝xi处的节点值，wi(x)是具有紧支特性的权函数；
步骤3.3：解出待定参数α(x)，让公式(15)中的J取得最小值；
对J取极值即对α(x)求导得
Figure BDA0002918689120000042
α(x)＝A-1(x)B(x)y (15)
Figure BDA0002918689120000043
B(x)＝[w(x1)p(x1)，w(x2)p(x2)，......，w(xn)p(xn)] (17)
计算时需要用到形函数
Figure BDA0002918689120000044
Figure BDA0002918689120000045
步骤3.4：获得最终的拟合函数：
Figure BDA0002918689120000046
另一方面，本发明提供一种基于台阶线提取的边坡监测方法，包括如下步骤：
步骤1：上述基于点云数据的露天矿台阶线提取方法进行台阶线提取；
步骤2：通过对两期露天矿台阶线提取结果进行分析，计算边坡曲面上每一部分的坡度值和位移变化量。
使用移动最小二乘法拟合台阶线后，将边坡曲面上每一部分对应的台阶线特征点坐标带入公式中计算坡度与位移变化量；
坡度，即两点的高程差与其水平距离的百分比，其计算公式如下：
h＝zi+1-zi (20)
Figure BDA0002918689120000047
Figure BDA0002918689120000048
式中：h为高程差，l为水平距离，x，y，z为台阶线特征点坐标；
位移变化量计算公式如下：
Figure BDA0002918689120000051
式中，x，y，z为台阶线特征点坐标。
采用上述技术方案所产生的有益效果在于：
1、本发明提供的方法通过渐进形态学滤波来消除挡土墙等人工地物对露天矿台阶线提取的干扰，提高了精度。在未经过滤波处理时提取露天矿台阶线，挡土墙的脊点会被提取出来，该点并不是台阶线特征点。所以使用渐进形态学滤波消除干扰。
2、本发明根据台阶线特征点的几何属性在提取特征点前先计算每点曲率，将曲率大于阈值的点作为潜在特征点，再对潜在特征点进行细化。这样大大减少了计算复杂度，提高了整个运行速度。
3、本发明通过对两期台阶线提取结果分析可以得出露天矿边坡坡度变化和边坡位移量，这对露天矿的安全生产具有重要意义。
附图说明
图1为本发明实施例中基于点云数据的露天矿台阶线提取和边坡监测的方法流程图；
图2(a)为本发明实施例中点P位于折叠边缘时的空间向量在拟合平面上的投影情况示意图；
图2(b)为本发明实施例中点P位于局部拟合平面内部时的空间向量在拟合平面上的投影情况示意图。
图3为本发明实施例中台阶线提取效果图。
具体实施方式
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
本实施例中基于点云数据的露天矿台阶线提取和边坡监测的方法如图1所示，其中基于点云数据的露天矿台阶线提取方法，包括如下步骤：
步骤1：对露天矿点云数据使用渐进形态学滤波算法进行预处理，通过逐渐增加滤波器的窗口大小和不断增大的高程差阈值移除点云数据中的地物数据，提取地面点，过程如下：
步骤1.1：渐进形态学滤波算法在进行第一次迭代时，将最小高程表面与初始滤波器的窗口大小作为输入，对三维激光点云数据作滤波运算，得到已滤波表面；
步骤1.2：比较已滤波表面在滤波前后的高程差是否小于高程差阈值，若小于则将窗口大小赋值给标记数组，来增加窗口大小；
所述高程差阈值根据研究区域的地形坡度确定，方法如下：
假设坡度恒定，则地形的最大高程差，窗口大小和地形坡度之间存在下列公式：
Figure BDA0002918689120000061
Figure BDA0002918689120000062
其中，dh0表示初始高程差阈值；s表示斜率；c表示格网大小；dhmax表示最大高程差阈值；wk表示运行到第k次窗口大小，dhT，K表示高程阈值，dhmax(t)，K表示地形的最大高程差。
步骤1.3：将上一迭代中获得的已滤波表面和增加的窗口大小作为滤波器的输入，继续进行滤波运算迭代，直到窗口大小为最大窗口。
所述地物数据包括：车辆、植被、建筑、挡土墙等。消除地物对露天矿台阶线提取的干扰，提高了精度。
步骤2：基于顾及邻域几何属性的三维边缘检测与曲率指数加权方法提取台阶线特征点，过程如下：
台阶线上的特征点均是具有曲率极值的点，因此本发明根据曲率初步筛选潜在特征点。
曲率是描述几何体弯曲程度的量，例如二维曲线偏离直线的程度，或者三维曲面偏离平面的程度，也可以用来确定曲面类型。
平均曲率、主曲率和高斯曲率是曲率的三个基本要素，(1)平均曲率：是指曲面上过一点P有无数个不同方向的曲线，不同方向具有不同的曲率，任意取其中一对相互垂直的正交曲率，取平均值。其计算公式为：甘＝(K1+K2)/2，式中K1，K2为点P的一对正交曲率。(2)主曲率：过曲面上一点P上具有无穷个正交曲率，其中存在一条曲线使得该曲线的曲率为极大，这个曲率为极大值Kmax，垂直于极大曲率面的曲率为极小值Kmin，极值方向为主方向。(3)高斯曲率：其计算公式为K＝Kmax·Kmin，表示的是点P的弯曲程度。
步骤2.1：点P为点云中任意一点，对点云中点P及其邻域内的点使用最小二乘法拟合曲面，并对其进行参数化：r(x，y)＝(x，y，ax2+bxy+cy2)，其中，x为点P的横坐标，y为点P的纵坐标，a、b、c为对曲面参数化后的参数；
步骤2.2：设r(x，y)的偏导数分别为rx，ry，rxx，rxy，ryy，然后根据以下公式计算平均曲率H和高斯曲率K：
Figure BDA0002918689120000063
Figure BDA0002918689120000071
其中，
Figure BDA0002918689120000072
为曲面的单位法向量；E＝rx·rx＝1+(2ax+by)2；
F＝rx·ry＝(2ax+by)(2cy+bx)，G＝ry·ry＝1+(2cy+bx)2，
Figure BDA0002918689120000073
步骤2.3：最后通过下列公式计算主曲率：
Figure BDA0002918689120000074
Figure BDA0002918689120000075
其中，kmax为最大主曲率，kmin为最小主曲率；
步骤2.4：计算点云中点P的曲率指数：
Figure BDA0002918689120000076
曲率指数与坐标系无关，具有旋转和位移不变性，在三维特征点检测领域中被广泛使用。
步骤2.5：将曲率指数大于曲率指数阈值的点定义为潜在特征点，对潜在特征点使用顾及邻域几何属性的三维边缘检测算法AGPN细化特征点，将潜在特征点中的内点集进行法向量优化；将潜在特征点中的外点集直接定义为非特征点；
AGPN算法引入了平面模型的RANSAC算法，RANSAC算法可以从一组包含“局外点”的观测数据集中通过迭代方式估计数学模型的参数，并且该数学模型能够适用于所有的“局内点”，它能从包含大量局外点的数据集中估计出高精度的参数，具有很好的鲁棒性。
步骤2.5的过程如下：
步骤2.5.1：对曲率指数大于曲率指数阈值的点其中任意一点P使用k-d tree算法检索其邻域点集S；
步骤2.5.2：对邻域点集S使用RANSAC算法拟合平面，将点集S分为内点集和外点集；
露天矿地形复杂，少有规则的曲面，当检测不规则曲面相交构成的边缘时，由于曲面的微小局部可以看成是平面，平面模型可以有效地近似局部曲面的几何结构。如果潜在特征点中一点P是折叠边缘，则使用k-d树搜索邻域后，空间向量
Figure BDA0002918689120000077
在点P邻域拟合平面上的投影会有较大的夹角，如图2(a)所示，该折叠边缘同时位于拟合平面和邻域内另一个相交表面上。当点P不是折叠边缘但也属于内点时，那么该点位于内点的拟合平面的内部，此时，空间向量
Figure BDA0002918689120000078
在拟合平面上的投影之间不存在较大的夹角，是连续分布的，如图2(b)所示。此外，若内点集的点数过少，则认为点P为噪声点。
步骤2.5.3：若点P在内点集中，则计算拟合平面的法向量，将其作为点P的法向量进行法向量优化；并将潜在特征点中的外点集直接定义为非特征点。
法向量优化过程为：对点云中某点P的法向量估计方法是对点P使用k-d树搜索邻域，对其邻域采用最小二乘法拟合平面，构造协方差矩阵，求解协方差矩阵的特征值和特征向量，然后选择最小特征值对应的特征向量，并进行单位化，则该向量为点P法向量。当点P的邻域有多个平面时，优化后的法向量将会垂直于其中一个平面。
步骤2.6：根据优化后的法向量和拟合平面上一对相互垂直的向量建立空间直角坐标系，设空间向量
Figure BDA0002918689120000081
在x轴上，向量
Figure BDA0002918689120000082
在y轴上，优化后的法向量为
Figure BDA0002918689120000083
计算角度阶跃Gθ，公式如下：
Figure BDA0002918689120000084
Gθ＝max(θi+1-θi)(i＝1…N-1) (9)
式中，空间向量
Figure BDA0002918689120000085
是由当前点P和其邻域内拟合平面提取的内点构成的，θi为向量
Figure BDA0002918689120000086
与向量
Figure BDA0002918689120000087
夹角；
步骤2.7：将角度阶跃大于阈值的点作为提取的台阶线特征点。
步骤3：采用移动最小二乘法拟合台阶线，得到提取后的露天矿台阶线；
传统的最小二乘法拟合曲线是在计算偏差后，让平方和达到最小值，从而将问题转化为求解线性方程组，最终得到拟合曲线。该方法具有全局性，但是不能满足定位的需要，不适用于具有大量数据集的模型。与传统的最小二乘法相比，移动最小二乘法在最终拟合曲线的定位和光滑性方面有较好的效果，且拟合精度也比较高。移动最小二乘法有两个明显的优点：第一点是不同于传统方法中的多项式，拟合过程中使用的函数由基函数和系数向量组成；第二点是引入了紧支的概念，即点x处的值y只受x支撑域内节点的影响。
步骤3的过程如下：
步骤3.1：对提取出的台阶线特征点进行曲线拟合，在要进行曲线拟合的部分中，其目标子区域在拟合过程中使用的函数记作f(x)：
Figure BDA0002918689120000088
α(x)＝[a1(x)，ax(x)，......，am(x)]T (11)
pT(x)＝[p1(x)，p2(x)，......，pm(x)] (12)
其中，α(x)代表的为待定的系数矢量；pT(x)代表的为m次多项式空间的一组基，ai(x)为待定的系数矢量，pi(x)为基函数；
步骤3.2：若计算点x的邻域内有n个节点，近似函数f(x)在这些节点处的误差的加权平方和如下：
Figure BDA0002918689120000091
其中，yi是x＝xi处的节点值，wi(x)是具有紧支特性的权函数；
步骤3.3：解出待定参数α(x)，让公式(15)中的J取得最小值；
对J取极值即对α(x)求导得：
Figure BDA0002918689120000092
α(x)＝A-1(x)B(x)y (15)
Figure BDA0002918689120000093
B(x)＝[w(x1)p(x1)，w(x2)p(x2)，......，w(xn)p(xn)] (17)
计算时需要用到形函数
Figure BDA0002918689120000094
Figure BDA0002918689120000095
步骤3.4：获得最终的拟合函数：
Figure BDA0002918689120000096
在本实施例中，某露天矿的台阶线提取效果如图3所示。
另一方面，本发明提供一种基于台阶线提取的边坡监测方法，包括如下步骤：
步骤1：上述基于点云数据的露天矿台阶线提取方法进行台阶线提取；
步骤2：通过对两期露天矿台阶线提取结果进行分析，计算边坡曲面上每一部分的坡度值和位移变化量。
使用移动最小二乘法拟合台阶线后，将边坡曲面上每一部分对应的台阶线特征点坐标带入公式中计算坡度与位移变化量；
坡度是地形描述中比较常用的参数，是一个具有大小与方向的矢量，坡度是指曲面上某一点的垂直方向Z与法线方向N之间的夹角值。
坡度，即两点的高程差与其水平距离的百分比，其计算公式如下：
h＝zi+1-zi (20)
Figure BDA0002918689120000101
Figure BDA0002918689120000102
式中：h为高程差，l为水平距离，x，y，z为台阶线特征点坐标；
位移变化量计算公式如下：
Figure BDA0002918689120000103
式中，x，y，z为台阶线特征点坐标。

Claims (2)

1.一种基于点云数据的露天矿台阶线提取方法，其特征在于，包括如下步骤：
步骤1：对露天矿点云数据使用渐进形态学滤波算法进行预处理，通过逐渐增加滤波器的窗口大小和不断增大的高程差阈值移除点云数据中的地物数据，提取地面点，过程如下：
步骤1.1：渐进形态学滤波算法在进行第一次迭代时，将最小高程表面与初始滤波器的窗口大小作为输入，对三维激光点云数据作滤波运算，得到已滤波表面；
步骤1.2：比较已滤波表面在滤波前后的高程差是否小于高程差阈值，若小于则将窗口大小赋值给标记数组，来增加窗口大小；
所述高程差阈值根据研究区域的地形坡度确定，方法如下：
假设坡度恒定，则地形的最大高程差，窗口大小和地形坡度之间存在下列公式：
Figure FDA0003842392910000011
Figure FDA0003842392910000012
其中，dh0表示初始高程差阈值；s表示斜率；c表示格网大小；dhmax表示最大高程差阈值；wk表示运行到第k次窗口大小，dhT,K表示高程差阈值，dhmax(t),k表示地形的最大高程差；
步骤1.3：将上一迭代中获得的已滤波表面和增加的窗口大小作为滤波器的输入，继续进行滤波运算迭代，直到窗口大小为最大窗口；
步骤2：基于顾及邻域几何属性的三维边缘检测与曲率指数加权方法提取台阶线特征点，过程如下：
步骤2.1：点P为点云中任意一点，对点云中点P及其邻域内的点使用最小二乘法拟合曲面，并对其进行参数化：r(x,y)＝(x,y,ax2+bxy+cy2)，其中，x为点P的横坐标，y为点P的纵坐标，a、b、c为对曲面参数化后的参数；
步骤2.2：设r(x,y)的偏导数分别为rx，ry，rxx，rxy，ryy，然后根据以下公式计算平均曲率H和高斯曲率K：
Figure FDA0003842392910000013
Figure FDA0003842392910000014
其中，
Figure FDA0003842392910000015
Figure FDA0003842392910000016
为曲面的单位法向量；E＝rx·rx＝1+(2ax+by)2；
F＝rx·ry＝(2ax+by)(2cy+bx)，G＝ry·ry＝1+(2cy+bx)2，
Figure FDA0003842392910000017
步骤2.3：最后通过下列公式计算主曲率：
Figure FDA0003842392910000021
Figure FDA0003842392910000022
其中，kmax为最大主曲率，kmin为最小主曲率；
步骤2.4：计算点云中点P的曲率指数：
Figure FDA0003842392910000023
步骤2.5：将曲率指数大于曲率指数阈值的点定义为潜在特征点，对潜在特征点使用顾及邻域几何属性的三维边缘检测算法AGPN细化特征点，将潜在特征点中的内点集进行法向量优化；将潜在特征点中的外点集直接定义为非特征点，过程如下：
步骤2.5.1：对曲率指数大于曲率指数阈值的点其中任意一点P使用k-d tree算法检索其邻域点集S；
步骤2.5.2：对邻域点集S使用RANSAC算法拟合平面，将点集S分为内点集和外点集；
步骤2.5.3：若点P在内点集中，则计算拟合平面的法向量，将其作为点P的法向量进行法向量优化；并将潜在特征点中的外点集直接定义为非特征点；
步骤2.6：根据优化后的法向量和拟合平面上一对相互垂直的向量建立空间直角坐标系，设空间向量
Figure FDA0003842392910000024
在x轴上，向量
Figure FDA0003842392910000025
在y轴上，优化后的法向量为
Figure FDA0003842392910000026
计算角度阶跃Gθ，公式如下：
Figure FDA0003842392910000027
Gθ＝max(θi+1-θi)i＝1,…N-1 (9)
式中，空间向量
Figure FDA0003842392910000028
是由当前点P和其邻域内拟合平面提取的内点构成的，θi为向量
Figure FDA0003842392910000029
与x轴的夹角；
步骤2.7：将角度阶跃大于阈值的点作为提取的台阶线特征点；
步骤3：采用移动最小二乘法拟合台阶线，得到提取后的露天矿台阶线，过程如下：
步骤3.1：对提取出的台阶线特征点进行曲线拟合，在要进行曲线拟合的部分中，其目标子区域在拟合过程中使用的函数记作f(x)：
Figure FDA00038423929100000210
α(x)＝[a(x1),a(x2),……,a(xm)]T (11)
pT(x)＝[p(x1),p(c2),……,p(xm)] (12)
其中，α(x)代表的为拟合曲线的一组系数；pT(x)代表的为m次多项式空间的一组基，a(xi)为待定的第i个系数，p(xi)为基函数；
步骤3.2：若计算点x的邻域内有n个节点，近似函数f(x)在这些节点处的误差的加权平方和如下：
Figure FDA0003842392910000031
其中，yi是x＝xi处的节点值，w(xi)是具有紧支特性的权函数；
步骤3.3：解出待定参数α(x)，让公式(13)中的J取得最小值；
对J取极值即对α(x)求导得：
Figure FDA0003842392910000032
α(x)＝A-1(x)B(x)y (15)
Figure FDA0003842392910000033
B(x)＝[w(x1)p(x1),w(x2)p(x2),……,w(xn)p(xn)] (17)
计算时需要用到形函数
Figure FDA0003842392910000034
Figure FDA0003842392910000035
步骤3.4：获得最终的拟合函数：
Figure FDA0003842392910000036
2.一种基于台阶线提取的边坡监测方法，其特征在于，包括如下步骤：
步骤1：采用权利要求1所述的基于点云数据的露天矿台阶线提取方法进行台阶线提取；
步骤2：通过对两个不同时期的露天矿台阶线提取结果进行分析，计算边坡曲面上每一部分的坡度值和位移变化量；
使用移动最小二乘法拟合台阶线后，将边坡曲面上每一部分对应的台阶线特征点坐标带入公式中计算坡度与位移变化量；
坡度i，即两点的高程差与其水平距离的百分比，其计算公式如下：
h＝zi+1-zi (20)
Figure FDA0003842392910000037
Figure FDA0003842392910000038
式中：h为高程差，l为水平距离，x，y，z为台阶线特征点坐标；
位移变化量计算公式如下：
Figure FDA0003842392910000039
式中，x，y，z为台阶线特征点坐标。
CN202110109374.7A 2021-01-27 2021-01-27 一种基于点云数据的露天矿台阶线提取和边坡监测的方法 Active CN112945196B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110109374.7A CN112945196B (zh) 2021-01-27 2021-01-27 一种基于点云数据的露天矿台阶线提取和边坡监测的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110109374.7A CN112945196B (zh) 2021-01-27 2021-01-27 一种基于点云数据的露天矿台阶线提取和边坡监测的方法

Publications (2)

Publication Number Publication Date
CN112945196A CN112945196A (zh) 2021-06-11
CN112945196B true CN112945196B (zh) 2022-11-04

Family

ID=76237635

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110109374.7A Active CN112945196B (zh) 2021-01-27 2021-01-27 一种基于点云数据的露天矿台阶线提取和边坡监测的方法

Country Status (1)

Country Link
CN (1) CN112945196B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115861296B (zh) * 2023-02-14 2023-05-16 湖北工业大学 基于无人机点云的高陡边坡危岩体自动识别方法及系统
CN117036621B (zh) * 2023-10-08 2024-02-09 山东瑞鑫时空信息科技有限公司 基于物联网的地理信息测绘仪数据管理方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105275471B (zh) * 2015-11-09 2017-11-21 安徽马钢工程技术集团有限公司 一种确定露天矿山开采境界的方法
CN107037496A (zh) * 2017-04-01 2017-08-11 重庆地质矿产研究院 一种露天矿山的实地动态检测方法
US10839264B2 (en) * 2018-11-09 2020-11-17 International Business Machines Corporation Scalable feature classification for laser scanning data and digital elevation models
CN109508508B (zh) * 2018-12-08 2024-03-22 河北省地质矿产勘查开发局国土资源勘查中心(河北省矿山和地质灾害应急救援中心) 一种露天矿山治理勘查设计方法
CN110047141B (zh) * 2019-03-02 2023-07-14 长沙迪迈数码科技股份有限公司 露天矿山开采现状三维模型更新方法
CN110751315B (zh) * 2019-09-23 2023-09-19 中南大学 露天矿道路系统人机交互式选线方法、系统及控制器

Also Published As

Publication number Publication date
CN112945196A (zh) 2021-06-11

Similar Documents

Publication Publication Date Title
CN106780524B (zh) 一种三维点云道路边界自动提取方法
CN107526360B (zh) 一种未知环境下排爆机器人多阶自主导航探测系统及方法
CN108958282B (zh) 基于动态球形窗口的三维空间路径规划方法
CN112184736B (zh) 一种基于欧式聚类的多平面提取方法
CN112945196B (zh) 一种基于点云数据的露天矿台阶线提取和边坡监测的方法
CN113269094B (zh) 基于特征提取算法和关键帧的激光slam系统及方法
CN111340723B (zh) 一种地形自适应的机载LiDAR点云正则化薄板样条插值滤波方法
CN102169581A (zh) 一种基于特征向量的快速高精度鲁棒性匹配方法
CN110335352B (zh) 一种机载激光雷达点云的双基元多分辨率层次滤波方法
CN114463338B (zh) 一种基于图割与后处理的建筑物激光脚点自动提取方法
CN114549549B (zh) 一种动态环境下基于实例分割的动态目标建模跟踪方法
CN109242786B (zh) 一种适用于城市区域的自动化形态学滤波方法
Sun et al. Automated segmentation of LiDAR point clouds for building rooftop extraction
CN117029870A (zh) 一种基于路面点云的激光里程计
CN116079713A (zh) 一种钻锚机器人多钻臂协同控制方法、系统、设备及介质
CN115131571A (zh) 一种基于点云预处理六领域的建筑物局部特征点识别方法
CN113554705A (zh) 一种变化场景下的激光雷达鲁棒定位方法
CN114897967A (zh) 一种面向挖掘装备自主作业的物料形态识别方法
CN111127474B (zh) 机载LiDAR点云辅助的正射影像镶嵌线自动选取方法及系统
CN112815944A (zh) 一种基于边角联合特征结构体的激光反射板定位方法
CN112348950A (zh) 一种基于激光点云分布特性的拓扑地图节点生成方法
Kaleci et al. Plane segmentation of point cloud data using split and merge based method
Luo et al. Comparison of an L1-regression-based and a RANSAC-based planar segmentation procedure for urban terrain data with many outliers
Skinner et al. Detection and segmentation of underwater archaeological sites surveyed with stereo-vision platforms
Zou et al. An adaptive strips method for extraction buildings from light detection and ranging data

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant