CN113362459B - 一种岩质滑坡破坏边界的特征提取方法 - Google Patents

一种岩质滑坡破坏边界的特征提取方法 Download PDF

Info

Publication number
CN113362459B
CN113362459B CN202110361027.3A CN202110361027A CN113362459B CN 113362459 B CN113362459 B CN 113362459B CN 202110361027 A CN202110361027 A CN 202110361027A CN 113362459 B CN113362459 B CN 113362459B
Authority
CN
China
Prior art keywords
structural surface
points
structural
landslide
boundary
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
CN202110361027.3A
Other languages
English (en)
Other versions
CN113362459A (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.)
Chengdu Univeristy of Technology
Original Assignee
Chengdu Univeristy of Technology
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 Chengdu Univeristy of Technology filed Critical Chengdu Univeristy of Technology
Priority to CN202110361027.3A priority Critical patent/CN113362459B/zh
Publication of CN113362459A publication Critical patent/CN113362459A/zh
Application granted granted Critical
Publication of CN113362459B publication Critical patent/CN113362459B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2200/00—Indexing scheme for image data processing or generation, in general
    • G06T2200/04—Indexing scheme for image data processing or generation, in general involving 3D image data
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/23—Dune restoration or creation; Cliff stabilisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Geometry (AREA)
  • Computer Graphics (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Remote Sensing (AREA)
  • Image Analysis (AREA)
  • Length Measuring Devices With Unspecified Measuring Means (AREA)

Abstract

本发明提供了一种岩质滑坡破坏边界的特征提取方法，属于岩质滑坡技术领域。包含采用地面三维激光扫描仪采集滑坡边界及影响区数据，通过Hough算法计算结构面法向量以确定其空间位置，将捕获的结构面进行HIS色彩重建并依据不同色彩种数明确结构面组数，基于彩色点云模型基础上密度带噪声空间聚类(DBSCAN)算法对结构面聚类分析，从而获得结构面参数进一步提取滑坡破坏边界特征。本发明相比传统法更加智能，测点区域覆盖广更具代表性，数据处理周期短。对岩质滑坡破坏边界特征提取在滑坡边界成因分析中具有极大的研究价值。

Description

一种岩质滑坡破坏边界的特征提取方法
技术领域
本发明属于岩质滑坡技术领域，尤其涉及一种岩质滑坡破坏边界的特征提 取方法。
背景技术
滑坡边界按传统方法由人工现场确定，考虑到安全性和实际工作量可能比 较大，因而人工法对滑坡所处地形条件要求较高。具体体现在高陡滑坡特别是 岩质滑坡的两侧边界及后缘，往往因其险要地势对野外调查者造成巨大阻碍； 同时在进行数据处理工作中，人工数据需要经历记录、输入等多个前期步骤， 使得工作量加大。为了解决这类问题，提出一套基于地面三维激光点云的岩质 滑坡破坏边界特征提取方法。
发明内容
针对现有技术中的上述不足，本发明提供的一种岩质滑坡破坏边界的特征 提取方法，解决了人工提取岩质滑坡破坏边界的特征工作量大、坡体盲区工作 者难以到达、全覆盖人工排查危险性高的问题。
为了达到以上目的，本发明采用的技术方案为：
本方案提供一种岩质滑坡破坏边界的特征提取方法，包括以下步骤：
S1、选择岩质滑坡不同部位布置人工测点，以及在滑区内设置若干台三维 激光扫描站点，构建三维点云模型；
S2、在所述人工测点和三维激光扫描站点处测量结构面产状数据；
S3、根据所述结构面产状数据对岩质滑坡边界进行点云数据分析；
S4、根据所述三维点云模型，利用随机霍夫变换计算得到所述点云数据的 法向量，并利用所述点云数据的法向量计算得到待检测结构面产状数据，将所 述待检测结构面产状数据与步骤S2中测量得到的结构面产状数据进行校核；
S5、根据校核后的结构面产状数据，利用HIS色彩重建方法确定结构面组 数；
S6、根据所述结构面组数，对结构面进行聚类分析；
S7、分别提取聚类后结构面的间距、延续性和粗糙度，完成岩质滑坡破坏 边界的特征提取。
本发明的有益效果是：本发明公开了一种岩质滑坡破坏边界特征提取方法。 包含采用地面三维激光扫描仪采集滑坡边界及影响区数据，通过Hough算法计 算结构面法向量以确定其空间位置，将捕获的结构面进行HIS色彩重建并依据 不同色彩种数明确结构面组数，基于彩色点云模型基础上密度带噪声空间聚类 (DBSCAN)算法对结构面聚类分析，从而获得结构面参数进一步提取滑坡破坏边 界特征。这种方法相比传统法更加智能，测点区域覆盖广更具代表性，数据处 理周期短，对岩质滑坡破坏边界特征提取在滑坡边界成因分析中具有极大的研 究价值。
进一步地，所述步骤S1包括以下步骤：
S101、在岩质滑坡边界上布置m个人工测点，在岩质滑坡影响区域内布置 3m个人工测点；
S102、在岩质滑坡区域内设置多个三维激光扫描站点；
S103、将三维激光扫描站点经数据拼接和坐标校准，建立三维点云模型。
上述进一步方案的有益效果是：在岩质滑坡边界上和影响区域内布设的人 工测点测结构面产状一方面用于三维激光扫描数据较核；另一方面用于与点云 数据对比，以判别这些结构面是由区域统一应力作用下形成的先期结构面，还 是与这该滑坡有关的新生断裂面。
再进一步地，所述步骤S4包括以下步骤：
S401、根据所述三维点云模型，在三维激光扫描数据的坐标系中设Y轴正 方向与北方向N相一致，设X轴正方向与东方向E相一致，设Z轴正方向为垂 直方向H；
S402、利用随机霍夫变换算法计算点云数据的法向量：
Ax+By+Cz+D＝0[A,B,C,D]∈R
其中，(x,y,z)表示点云数据坐标，(A,B,C)表示构成结构面平面法向量
Figure BDA0003005528820000031
D 表示平面的空间位置，R表示实数集。
S403、根据所述点云数据的法向量，分别计算得到结构面走向线与北方向N 夹角的关系以及结构面的倾角，并将所述结构面走向线与北方向N夹角的关系 以及结构面的倾角作为待检测结构面产状数据；
所述结构面走向线与北方向N夹角β的表达式如下：
Figure BDA0003005528820000032
所述结构面的倾角α的表达式如下：
Figure BDA0003005528820000033
其中，(A,B,C)构成结构面平面法向量
Figure BDA0003005528820000034
S404、将所述待检测结构面产状数据与步骤S2中测量得到的结构面产状数 据进行校核。
上述进一步方案的有益效果是：将经三维激光扫描数据展现在新的坐标系 下，以对点云的空间位置把控更加精准，利用随机霍夫变换算法计算点云数据 的法向量，通过法向量对结构面定位在清晰明了试图效果的基础上能极大简化 后续工作量。
再进一步地，所述步骤S6包括以下步骤：
S601、根据所述结构面组数，利用基于密度的带噪声空间聚类算法对岩质 滑坡破坏边界点云结构面进行分组，同时输入距离度量Eps和数量度量Min-pts 两个参数；
S602、由三维激光扫描接收器接收到的激光脉冲，利用时差原理计算得到 三维坐标点；
S603、根据分组结果，将所述三维坐标点按疏密程度分为若干团块；
S604、在某个团块之内点进行相邻点距计算，得到相邻点平均距离；
S605、将相邻点平均距离与标准差记为距离度量Eps；
S606、对所述数量度量Min-pts进行调整，并判断某点γ在距离度量Eps范 围内所包含的点个数是否达到数量度量Min-pts，且所有在距离度量Eps范围内 点个数总和达到阈值，若是，则数量度量Min-pts有效，并进入步骤S607，否 则，持续步骤S606；
S607、将所有点标记为核心点、边界点或噪声点，其中，所述核心点为在 距离度量Eps内含有超过数量度量Min-pts数量的点，所述边界点为在距离度量 Eps内点的数量小于数量度量Min-pts数量的点，所述噪声点为核心点和边界点 以外的点；
S608、删除所述噪声点，将距离不超过距离度量Eps的核心点相互连接， 并与纳入核心点领域内的点形成一个簇；
S609、删除点云数据小于100的簇，完成对结构面的聚类分析。
上述进一步方案的有益效果是：滑坡边界上出露的岩体结构面具有出露面 积差异大、几何形态不规则、空间分布随机性强等特点，基于密度的带噪声空 间聚类(DBSCAN)算法不需要预先设定形成的簇类的数量，可以对任意形状点云 进行聚类，对数据样本顺序不敏感，对噪声点具有很好的鲁棒性。
再进一步地，所述步骤S7中结构面间距提取的步骤如下：
A1、根据聚类后的结构面，利用最小二乘法计算得到各簇平面的结构面空 间位置Dij：
Figure BDA0003005528820000051
其中，Aij、Bij和Cij均表示k点的向量参数，n表示k的个数，
Figure BDA0003005528820000052
和
Figure BDA0003005528820000053
均表示一个簇中k点的坐标；
A2、将各所述簇平面的结构面空间位置以从小到大进行排序，并计算得到 相邻簇之间的间距；
A3、判断是否计算完所有相邻簇之间的间距，若是，则进入步骤A4，否则， 返回步骤A2；
A4、利用正态分布拟合得到所有相邻簇之间的间距，并以期望μ作为结构 面间距的代表值，完成结构面间距的提取。
上述进一步方案的有益效果是：结构面间距是衡量结构面发育程度的重要 指标，反映了岩体结构的完整性。通过将结构面空间位置量化为筛选工作提供 了很好的基础，由判断语句及相邻簇间距正态分布期望值作为结构面间距的提 取使得结果更具代表性。
再进一步地，所述步骤S7中结构面延续性提取的步骤如下：
B1、根据聚类后的结构面，利用迹长表征结构面的延续性；
B2、以倾向方向为x轴，走向方法为y轴，建立坐标系，其中，oxyz为原 始坐标系，o′x′y′z′为变换后的坐标系，点(x′，y′，z′)为局部坐标；
B3、沿o′x′提取结构面倾向方向延续性，沿o′y′方向提取结构面走向方向延 续性，并计算得到最大延伸长度，完成对聚类后结构面的延续性提取。
上述进一步方案的有益效果是：在变换后坐标系下以结构面走向和倾向方 向迹长提取延续性可以使获取步骤得以规范和程序化。
再进一步地，所述步骤B3中结构面倾向方向延续性Ldip的表达式如下：
Ldip＝max(x′(i，j))-min(x′(i，j))
其中，x′(i，j)表示结构面沿变换坐标轴o′x′倾向方向提取的迹线长度；
所述结构面走向方向延续性Lstrike的表达式如下：
Lstrike＝max(y'(i,j))-min(y'(i,j))
其中，y′(i，j)表示结构面沿变换坐标轴o′y′走向方向提取的迹线长度；
所述最大延伸长度Lmax的表达式如下：
Lmax＝max(Ch(X(i，j)))
其中，X(i，j)表示结构面沿倾向和走向提取的迹线长度，Ch(·)表示MATLAB 中convhull函数。
上述进一步方案的有益效果是：空间中迹线与结构面倾向和倾角的关系具 有不确定性，用结构面的最大延伸长度(Lmax)的平均值表征每组结构面延续性 能降低点云中真实迹线自动提取的困难。
再进一步地，所述步骤S7中粗糙度提取的步骤如下：
C1、根据聚类后的结构面，利用相对起伏度Ra和伸长率R反映结构面粗糙 度JRC；
C2、根据所述结构面粗糙度JRC的初略值，将结构面剖面轮廓分别划分为 平直状剖面曲线、波浪状剖面曲线以及锯齿状剖面曲线；
C3、选取以滑动方向所在各测窗上10m为直径的圆形区域，圆内从0°开始， 以倾向方向为0～180°，按顺时针方向每30°间隔设置剖面线，并分别计算得到 平直状剖面曲线的粗糙度、波浪状剖面曲线的粗糙度以及锯齿状剖面曲线的粗 糙度，完成对聚类后结构面的粗糙度提取。
上述进一步方案的有益效果是：宏观断面粗糙度目前还没有统一的描述方 法，大尺度粗糙度的测量对于了解粗糙度的尺寸效应，确定滑动方向具有重要 意义。
再进一步地，所述步骤C1中结构面粗糙度JRC的表达式如下；
JRC＝8.5R+0.2(100Ra)2-2
Figure BDA0003005528820000071
Figure BDA0003005528820000072
其中，L表示结构面剖面曲线长度，L0表示结构面剖面直线长度，RA表示 结构面剖面起伏度。
再进一步地，所述步骤C3中平直状剖面曲线的粗糙度JRC1的表达式如下：
JRC1＝12.862R+10.315arctan(100Ra)-4.017
所述波浪状剖面曲线的粗糙度JRC2的表达式如下：
JRC2＝2.179R+0.585(100Ra)2+6.169
所述锯齿状剖面曲线的粗糙度JRC3的表达式如下：
JRC3＝10.544R+0.327(100Ra)2+10.377
附图说明
图1为本发明的方法流程图。
图2为本实施例中照壁岩滑坡边界地质构造特征示意图。
图3为本实施例中南侧断壁h区法向量计算及HSI色彩渲染示意图。
图4为本实施例中DBSCAN算法点云划分示意图。
图5为本实施例中贯通结构面间距(d')与非贯通结构面间距(D')示意图。
图6为本实施例中结构面延续性示意图。
图7为本实施例中伸长率R与相对起伏度Ra计算参数示意图。
图8为本实施例中照壁岩滑坡破坏边界测窗示意图。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理 解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的 普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精 神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保 护之列。
实施例
如图1所示，本发明提供了一种岩质滑坡破坏边界的特征提取方法，其实 现方法如下：
S1、选择岩质滑坡不同部位布置人工测点，以及在滑区内设置若干台三维 激光扫描站点，并构建三维点云模型，其实现方法如下：
S101、在岩质滑坡边界上布置m个人工测点，在岩质滑坡影响区域内布置 3m个人工测点；
S102、在岩质滑坡区域内设置多个三维激光扫描站点；
S103、将三维激光扫描站点经数据拼接和坐标校准，建立三维点云模型。
本实施例中，选择滑坡不同部位布置人工测点，其中滑坡边界上布点个数m， 滑坡影响区域内布点3m个(m为正整数，依据滑坡规模大小取值)。在滑区内 布设多台ILRIS·LR地面三维激光扫描站点，一共获取N个测点。将站点经过数 据拼接和坐标校准，去除表面植被，从而建立滑坡三维点云模型。
S2、根据所述三维点云模型，在所述人工测点和三维激光扫描站点处测量 结构面产状数据；
本实施例中，人工测量结构面产状数据用于对三维激光扫描站点扫描数据 进行校核。
S3、根据所述结构面产状数据对岩质滑坡边界进行点云数据分析；
S4、根据所述三维点云模型，利用随机霍夫变换计算得到所述点云数据的 法向量，并利用所述点云数据的法向量计算得到待检测结构面产状数据，将所 述待检测结构面产状数据与步骤S2中测量得到的结构面产状数据进行校核，其 实现方法如下：
S401、根据所述三维点云模型，在三维激光扫描数据的坐标系中设Y轴正 方向与北方向N相一致，设X轴正方向与东方向E相一致，设Z轴正方向为垂 直方向H；
S402、利用随机霍夫变换算法计算点云数据的法向量：
Ax+By+Cz+D＝0[A,B,C,D]∈R
其中，(x,y,z)表示点云数据坐标，(A,B,C)表示构成结构面平面法向量
Figure BDA0003005528820000091
D 表示平面的空间位置，R表示实数集。
S403、根据所述点云数据的法向量，分别计算得到结构面走向线与北方向N 夹角的关系以及结构面的倾角，并将所述结构面走向线与北方向N夹角的关系 以及结构面的倾角作为待检测结构面产状数据；
所述结构面走向线与北方向N夹角β的表达式如下：
Figure BDA0003005528820000092
所述结构面的倾角α的表达式如下：
Figure BDA0003005528820000101
其中，(A,B,C)构成结构面平面法向量
Figure BDA0003005528820000102
S404、将所述待检测结构面产状数据与步骤S2中测量得到的结构面产状数 据进行校核。
本实施例中，在三维激光扫描数据处理软件的坐标系中，Y轴正方向与北(N) 方向相一致，X轴正方向与东(E)方向相同，Z轴正方向为垂直方向(H)。从而可 以根据坐标轴与方位的对应关系，应用几何知识求解结构面的产状。空间中某 一平面方程可用式(1)所示，采用随机霍夫(Hough)变换算法计算点云法向量。 则结构面产状(倾向/倾角)可按公式(2)、(3)推导计算：
Ax+By+Cz+D＝0[A,B,C,D]∈R (1)
其中，(x,y,z)表示点云数据坐标，(A,B,C)表示构成结构面平面法向量
Figure BDA0003005528820000103
D 表示平面的空间位置，R表示实数集。
结构面走向线与北方向N夹角β的表达式如下：
Figure BDA0003005528820000104
结构面的倾角α的表达式如下：
Figure BDA0003005528820000105
其中，(A,B,C)构成结构面平面法向量
Figure BDA0003005528820000106
S5、根据所述结构面产状，利用HIS色彩重建方法确定结构面组数；
本实施例中，如图2-图3所示，为使滑坡破坏边界上结构面分布情况能够 直观展示，采用HSI色彩重建技术渲染三维结构面，将施密特投影网与HSI色 彩相结合，H对应倾向，S对应倾角，强度设置为最大值1，从而为空间中每一 个结构面赋予特定的颜色，经HSI色彩渲染后形成的滑坡破坏边界彩图，在滑 坡边界上依据不同色彩种数明确结构面组数。
S6、根据所述结构面组数，对结构面进行聚类分析，其实现方法如下：
S601、根据所述结构面组数，利用基于密度的带噪声空间聚类算法对岩质 滑坡破坏边界点云结构面进行分组，同时输入距离度量Eps和数量度量Min-pts 两个参数；
S602、由三维激光扫描接收器接收到的激光脉冲，利用时差原理计算得到 三维坐标点；
S603、根据分组结果，将所述三维坐标点按疏密程度分为若干团块；
S604、在某个团块之内点进行相邻点距计算，得到相邻点平均距离；
S605、将相邻点平均距离与标准差记为距离度量Eps；
S606、对所述数量度量Min-pts进行调整，并判断某点γ在距离度量Eps范 围内所包含的点个数是否达到数量度量Min-pts，且所有在距离度量Eps范围内 点个数总和达到阈值，若是，则数量度量Min-pts有效，并进入步骤S607，否 则，持续步骤S606；
S607、将所有点标记为核心点、边界点或噪声点，其中，所述核心点为在 距离度量Eps内含有超过数量度量Min-pts数量的点，所述边界点为在距离度量 Eps内点的数量小于数量度量Min-pts数量的点，所述噪声点为核心点和边界点 以外的点；
S608、删除所述噪声点，将距离不超过距离度量Eps的核心点相互连接， 并与纳入核心点领域内的点形成一个簇；
S609、删除点云数据小于100的簇，完成对结构面的聚类分析。
本实施例中，基于密度的带噪声空间聚类(DBSCAN)算法用于岩质滑坡破坏 边界点云结构面的分组。首先输入两个参数：(a)Eps：距离度量，为核心点的邻 域半径；(b)Min-pts：数量度量，为聚类在一起点的最小数目。
本实施例中，参数设置：首先将所识别的点，按照疏密程度分为若干团块， 再在某个团块之内点进行相邻点距计算，将所得平均值加上标准差记为Eps。调 整Min-pts(记为ξ，为正整数)的值，当某点γ在Eps范围内所包含的点个数 达到ξ，且这样的γ点个数之和占团块总点数的百分之九十及以上则Min-pts值 有效，否则继续调整ξ使其达到要求。
本实施例中，如图4所示，DBSCAN聚类算法的基本步骤：将所有点标记 为核心点、边界点或噪声点，如图4所示，其中核心点为在Eps内含有超过Min-pts 数目的点；边界点为在Eps内点的数量小于Min-Pts，但是落在核心点的邻域内 的点；噪声点为既不是核心点也不是边界点的点；删除识别的噪声点；将距离 不超过Eps的核心点相互连接，并与纳入核心点邻域内的点，形成一个簇；删 除了点云数目小于100的簇，以剩余的簇作为结构面样本提取参数。
S7、分别提取聚类后结构面的间距、延续性和粗糙度，完成岩质滑坡破坏 边界的特征提取；
结构面间距提取的步骤如下：
A1、根据聚类后的结构面，利用最小二乘法计算得到各簇平面的结构面空 间位置Dij：
Figure BDA0003005528820000121
其中，Aij、Bij和Cij均表示k点的向量参数，n表示k的个数，
Figure BDA0003005528820000122
和
Figure BDA0003005528820000123
均表示一个簇中k点的坐标；
A2、将各所述簇平面的结构面空间位置以从小到大进行排序，并计算得到 相邻簇之间的间距；
A3、判断是否计算完所有相邻簇之间的间距，若是，则进入步骤A4，否则， 返回步骤A2；
A4、利用正态分布拟合得到所有相邻簇之间的间距，并以期望μ作为结构 面间距的代表值，完成结构面间距的提取。
结构面延续性提取的步骤如下：
B1、根据聚类后的结构面，利用迹长表征结构面的延续性；
B2、以倾向方向为x轴，走向方法为y轴，建立坐标系，其中，oxyz为原 始坐标系，o′x′y′z′为变换后的坐标系，点(x′，y′，z′)为局部坐标；
B3、沿o′x′提取结构面倾向方向延续性，沿o′y′方向提取结构面走向方向延 续性，并计算得到最大延伸长度，完成对聚类后结构面的延续性提取。
粗糙度提取的步骤如下：
C1、根据聚类后的结构面，利用相对起伏度Ra和伸长率R反映结构面粗糙 度JRC；
C2、根据所述结构面粗糙度JRC的初略值，将结构面剖面轮廓分别划分为 平直状剖面曲线、波浪状剖面曲线以及锯齿状剖面曲线；
C3、选取以滑动方向所在各测窗上10m为直径的圆形区域，圆内从0°开始， 以倾向方向为0～180°，按顺时针方向每30°间隔设置剖面线，并分别计算得到 平直状剖面曲线的粗糙度、波浪状剖面曲线的粗糙度以及锯齿状剖面曲线的粗 糙度，完成对聚类后结构面的粗糙度提取。
本实施例中，利用点云数据计算对滑坡破坏边界成因起重要作用的结构面 间距、延续性和粗糙度。公式(1)中的D值代表了结构面的空间位置，对使用DBSCAN算法聚类后的簇，采用最小二乘法，按公式(4)计算簇平面的D值。 把得到的D值从小到大排序，再计算相邻簇之间的间距，直到计算完所有簇。 用正态分布拟合得到的所有间距数据，以期望μ作为结构面间距的代表值，计 算中，如图5所示，区分了贯通结构面(d')和非贯通结构面(D')。
Figure BDA0003005528820000141
其中，Aij、Bij和Cij均表示k点的向量参数，n表示k的个数，
Figure BDA0003005528820000142
和
Figure BDA0003005528820000143
均表示一个簇中k点的坐标。
本实施例中，用迹长表征结构面延续性，迹长是岩体结构面与出露面相交 迹线的长度。对使用DBSCAN算法聚类后的簇，以倾向方向为x轴，走向方向 为y轴，建立坐标系，其中oxyz为原始坐标系，o′x′y′z′为变换后的坐标系，点 (x′，y′，z′)为局部坐标，如图6所示。倾向方向延续性Ldip沿o′x′方向提取，计算 公式为(5)；走向方向延续性沿o′y′方向提取，采用公式(6)计算；最大延伸 长度Lmax可以通过MATLAB中‘convhull’(简写为Ch)函数计算，见公式(7).
Ldip＝max(x′(i，j))-min(x′(i，j)) (5)
其中，x′(i，j)表示结构面沿变换坐标轴O′x′(倾向)方向提取的迹线长度。
Figure BDA0003005528820000144
其中，y′(i，j)表示结构面沿变换坐标轴O′y′(走向)方向提取的迹线长度。
Lmax＝max(Ch(X(i，j)))(7)
其中，X(i，j)表示结构面沿倾向和走向提取的迹线长度。Ch(·)表示MATLAB 中‘convhull’函数。
本实施例中，如图7所示，采用Barton提出的结构面粗糙度系数JRC来描 述结构面粗糙度。采用不同方向的剖面线定量描述断面的粗糙度。以Barton粗 糙度等级剖面曲线为基础，利用相对起伏度Ra和伸长率R反映结构面粗糙度(以 系数JRC表示)。JRC初略值按公式(10)计算，将结构面剖面轮廓按JRC值划 为平直状[0～8]、波浪状[8～16]、锯齿状[16～20]三类。再经过详细计算JRC值计 算公式见(11)～(13)。
伸长率R的表达式为：
Figure BDA0003005528820000151
相对伏度Ra的表达式为：
Figure BDA0003005528820000152
结构面粗糙度JRC初略值的表达式如下；
JRC＝8.5R+0.2(100Ra)2-2 (10)
其中，L表示结构面剖面曲线长度,L0表示结构面剖面直线长度,RA表示结构 面剖面起伏度。
平直状剖面曲线(0≤JRC≤8)，有
JRC1＝12.862R+10.315arctan(100Ra)-4.017 (11)
波浪状剖面曲线(8＜JRC≤16)，有
JRC2＝2.179R+0.585(100Ra)2+6.169 (12)
锯齿状剖面曲线(16＜JRC≤20)，有
JRC3＝10.544R+0.327(100Ra)2+10.377 (13)
其中，JRC结构面粗糙度系数，用来描述结构面粗糙度
本实施例中，粗糙度测量中剖面线的选取：采样间隔必须小于剖面长度的 2％，即当基于三维激光点云(有效点间距为S)时，粗糙度采样剖面需要大于50S。 通过数据检验，采样间距为0.2m时，可以得到较好的剖面线(例如在照壁岩滑 坡中剖面线长度设置为10m)，在滑坡破坏边界上选取10m为直径的圆形区域， 圆内从0°开始(以倾向方向为0～180°)，按顺时针方向每30°间隔设置剖面线，如 图8所示，计算每条剖面线的JRC值，从而反映滑坡破坏边界上不同方向的粗 糙度。

Claims (4)

1.一种岩质滑坡破坏边界的特征提取方法，其特征在于，包括以下步骤：
S1、选择岩质滑坡不同部位布置人工测点，以及在滑区内设置若干台三维激光扫描站点，构建三维点云模型；
S2、在所述人工测点和三维激光扫描站点处测量结构面产状数据；
S3、根据所述结构面产状数据对岩质滑坡边界进行点云数据分析；
S4、根据所述三维点云模型，利用随机霍夫变换计算得到所述点云数据的法向量，并利用所述点云数据的法向量计算得到待检测结构面产状数据，将所述待检测结构面产状数据与步骤S2中测量得到的结构面产状数据进行校核；
S5、根据校核后的结构面产状数据，利用HIS色彩重建方法确定结构面组数；
S6、根据所述结构面组数，对结构面进行聚类分析；
S7、分别提取聚类后结构面的间距、延续性和粗糙度，完成岩质滑坡破坏边界的特征提取；
所述步骤S4包括以下步骤：
S401、根据所述三维点云模型，在三维激光扫描数据的坐标系中设Y轴正方向与北方向N相一致，设X轴正方向与东方向E相一致，设Z轴正方向为垂直方向H；
S402、利用随机霍夫变换算法计算点云数据的法向量：
Ax+By+Cz+D＝0[A,B,C,D]∈R
其中，(x,y,z)表示点云数据坐标，(A,B,C)表示构成结构面平面法向量
Figure FDA0003800923740000011
D表示平面的空间位置，R表示实数集；
S403、根据所述点云数据的法向量，分别计算得到结构面走向线与北方向N夹角的关系以及结构面的倾角，并将所述结构面走向线与北方向N夹角的关系以及结构面的倾角作为待检测结构面产状数据；
所述结构面走向线与北方向N夹角β的表达式如下：
Figure FDA0003800923740000021
所述结构面的倾角α的表达式如下：
Figure FDA0003800923740000022
其中，(A,B,C)构成结构面平面法向量
Figure FDA0003800923740000023
S404、将所述待检测结构面产状数据与步骤S2中测量得到的结构面产状数据进行校核；
所述步骤S6包括以下步骤：
S601、根据所述结构面组数，利用基于密度的带噪声空间聚类算法对岩质滑坡破坏边界点云结构面进行分组，同时输入距离度量Eps和数量度量Min-pts两个参数；
S602、由三维激光扫描接收器接收到的激光脉冲，利用时差原理计算得到三维坐标点；
S603、根据分组结果，将所述三维坐标点按疏密程度分为若干团块；
S604、在某个团块之内点进行相邻点距计算，得到相邻点平均距离；
S605、将相邻点平均距离与标准差记为距离度量Eps；
S606、对所述数量度量Min-pts进行调整，并判断某点γ在距离度量Eps范围内所包含的点个数是否达到数量度量Min-pts，且所有在距离度量Eps范围内点个数总和达到阈值，若是，则数量度量Min-pts有效，并进入步骤S607，否则，持续步骤S606；
S607、将所有点标记为核心点、边界点或噪声点，其中，所述核心点为在距离度量Eps内含有超过数量度量Min-pts数量的点，所述边界点为在距离度量Eps内点的数量小于数量度量Min-pts数量的点，所述噪声点为核心点和边界点以外的点；
S608、删除所述噪声点，将距离不超过距离度量Eps的核心点相互连接，并与纳入核心点领域内的点形成一个簇；
S609、删除点云数据小于100的簇，完成对结构面的聚类分析；
所述步骤S7中结构面间距提取的步骤如下：
A1、根据聚类后的结构面，利用最小二乘法计算得到各簇平面的结构面空间位置Dij：
Figure FDA0003800923740000031
其中，Aij、Bij和Cij均表示k点的向量参数，n表示k的个数，
Figure FDA0003800923740000032
和
Figure FDA0003800923740000033
均表示一个簇中k点的坐标；
A2、将各所述簇平面的结构面空间位置以从小到大进行排序，并计算得到相邻簇之间的间距；
A3、判断是否计算完所有相邻簇之间的间距，若是，则进入步骤A4，否则，返回步骤A2；
A4、利用正态分布拟合得到所有相邻簇之间的间距，并以期望μ作为结构面间距的代表值，完成结构面间距的提取；
所述步骤S7中结构面延续性提取的步骤如下：
B1、根据聚类后的结构面，利用迹长表征结构面的延续性；
B2、以倾向方向为x轴，走向方法为y轴，建立坐标系，其中，oxyz为原始坐标系，o′x′y′z′为变换后的坐标系，点(x′，y′，z′)为局部坐标；
B3、沿o′x′提取结构面倾向方向延续性，沿o′y′方向提取结构面走向方向延续性，并计算得到最大延伸长度，完成对聚类后结构面的延续性提取；
所述步骤B3中结构面倾向方向延续性Ldip的表达式如下：
Ldip＝max(x′(i，j))-min(x′(i，j))
其中，x′(i，j)表示结构面沿变换坐标轴o′x′倾向方向提取的迹线长度；
所述结构面走向方向延续性Lstrike的表达式如下：
Lstrike＝max(y′(i，j))-min(y′(i，j))
其中，y′(i，j)表示结构面沿变换坐标轴o′y′走向方向提取的迹线长度；
所述最大延伸长度Lmax的表达式如下：
Lmax＝max(Ch(X(i，j)))
其中，X(i，j)表示结构面沿倾向和走向提取的迹线长度，Ch(·)表示MATLAB中convhull函数；
所述步骤S7中粗糙度提取的步骤如下：
C1、根据聚类后的结构面，利用相对起伏度Ra和伸长率R反映结构面粗糙度JRC；
C2、根据所述结构面粗糙度JRC的初略值，将结构面剖面轮廓分别划分为平直状剖面曲线、波浪状剖面曲线以及锯齿状剖面曲线；
C3、选取以滑动方向所在各测窗上10m为直径的圆形区域，圆内从0°开始，以倾向方向为0～180°，按顺时针方向每30°间隔设置剖面线，并分别计算得到平直状剖面曲线的粗糙度、波浪状剖面曲线的粗糙度以及锯齿状剖面曲线的粗糙度，完成对聚类后结构面的粗糙度提取。
2.根据权利要求1所述的岩质滑坡破坏边界的特征提取方法，其特征在于，所述步骤S1包括以下步骤：
S101、在岩质滑坡边界上布置m个人工测点，在岩质滑坡影响区域内布置3m个人工测点；
S102、在岩质滑坡区域内设置多个三维激光扫描站点；
S103、将三维激光扫描站点经数据拼接和坐标校准，建立三维点云模型。
3.根据权利要求1所述的岩质滑坡破坏边界的特征提取方法，其特征在于，所述步骤C1中结构面粗糙度JRC的表达式如下；
JRC＝8.5R+0.2(100Ra)2-2
Figure FDA0003800923740000051
Figure FDA0003800923740000052
其中，L表示结构面剖面曲线长度，L0表示结构面剖面直线长度，RA表示结构面剖面起伏度。
4.根据权利要求3所述的岩质滑坡破坏边界的特征提取方法，其特征在于，所述步骤C3中平直状剖面曲线的粗糙度JRC1的表达式如下：
JRC1＝12.862R+10.315arctan(100Ra)-4.017
所述波浪状剖面曲线的粗糙度JRC2的表达式如下：
JRC2＝2.179R+0.585(100Ra)2+6.169
所述锯齿状剖面曲线的粗糙度JRC3的表达式如下：
JRC3＝10.544R+0.327(100Ra)2+10.377。
CN202110361027.3A 2021-04-02 2021-04-02 一种岩质滑坡破坏边界的特征提取方法 Active CN113362459B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110361027.3A CN113362459B (zh) 2021-04-02 2021-04-02 一种岩质滑坡破坏边界的特征提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110361027.3A CN113362459B (zh) 2021-04-02 2021-04-02 一种岩质滑坡破坏边界的特征提取方法

Publications (2)

Publication Number Publication Date
CN113362459A CN113362459A (zh) 2021-09-07
CN113362459B true CN113362459B (zh) 2022-09-27

Family

ID=77525205

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110361027.3A Active CN113362459B (zh) 2021-04-02 2021-04-02 一种岩质滑坡破坏边界的特征提取方法

Country Status (1)

Country Link
CN (1) CN113362459B (zh)

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104183017A (zh) * 2014-07-29 2014-12-03 浙江大学 基于地面三维激光点云实现地质体产状自动提取的方法
CN106679579A (zh) * 2016-12-02 2017-05-17 中国电建集团昆明勘测设计研究院有限公司 一种移动式滑坡体变形监测装置及方法
WO2018201180A1 (en) * 2017-05-02 2018-11-08 PETRA Data Science Pty Ltd Automated, real time processing, analysis, mapping and reporting of data for the detection of geotechnical features
US10909650B2 (en) * 2017-06-23 2021-02-02 Cloud 9 Perception, LP System and method for sensing and computing of perceptual data in industrial environments
CN107860375A (zh) * 2017-11-03 2018-03-30 广州地理研究所 一种基于三维激光扫描技术的滑坡灾害体积快速提取方法
CN108489403B (zh) * 2018-06-08 2020-08-07 宁波大学 基于三维激光扫描的露天矿山边坡岩体节理产状快速精细取值方法
CN108489402B (zh) * 2018-06-08 2020-09-15 宁波大学 基于三维激光扫描的露天矿山边坡岩体节理规模快速精细取值方法
CN110288700A (zh) * 2019-06-26 2019-09-27 东北大学 一种岩质边坡结构面自动分组及位移预测方法
CN111553292B (zh) * 2020-04-30 2023-05-02 南京工业大学 一种基于点云数据的岩体结构面识别与产状分类方法
CN111340012B (zh) * 2020-05-19 2020-09-15 北京数字绿土科技有限公司 一种地质灾害解译方法、装置、终端设备
CN112529844B (zh) * 2020-11-24 2022-03-25 成都理工大学 一种基于三维激光扫描的岩体结构面识别与信息提取方法

Also Published As

Publication number Publication date
CN113362459A (zh) 2021-09-07

Similar Documents

Publication Publication Date Title
CN110030951B (zh) 一种基于三维激光扫描技术的引水竖井缺陷检测方法
CN111553292B (zh) 一种基于点云数据的岩体结构面识别与产状分类方法
Li et al. Terrestrial laser scanning assisted flatness quality assessment for two different types of concrete surfaces
CN110276732B (zh) 一种顾及地形特征线要素的山区点云空洞修复方法
CN109490072B (zh) 一种土木工程建筑用检测系统及其检测方法
CN106702995A (zh) 基于bim的岩土工程监测模型构建方法
CN107621231A (zh) 一种隧道二次衬砌厚度检测方法
CN108830317B (zh) 基于数字摄影测量的露天矿山边坡岩体节理产状快速精细取值方法
CN108801221B (zh) 基于数字摄影测量的露天矿山边坡岩体节理规模快速精细取值方法
CN109190593B (zh) 基于凹凸类别划分的山区道路沿线边坡稳定性初步判别方法
CN112989453B (zh) 一种基于bim的全息变形信息提取方法
CN114492984A (zh) 粉尘浓度的时空分布预测方法、装置、设备和存储介质
CN113362459B (zh) 一种岩质滑坡破坏边界的特征提取方法
CN111091049B (zh) 一种基于反向特征匹配的路面障碍物检测方法
CN110986815A (zh) 一种基于三维激光点云的隧道施工监控量测方法
CN114088015A (zh) 一种岩体三维裂隙网络模型快速智能生成和剖切方法
CN112967256B (zh) 一种基于空间分布的隧道椭圆化检测方法
CN113240637B (zh) 一种基于机器学习的墙面平整度信息化检测方法和系统
CN105043962B (zh) 一种定量测量砂岩质文物表面风化速度的方法
CN115455706A (zh) 考虑卸荷破裂效应的区域岩体质量评估方法及相关组件
CN115564725A (zh) 基于图像智能识别的边坡监测方法及装置
CN115035026A (zh) 一种基于三维点云信息的隧道爆破质量评估方法
CN109520461A (zh) 系列尺寸岩体结构面粗糙度试样的统计样本数确定方法
CN114001678B (zh) 基于车载激光雷达的路面平整度检测方法、装置和车辆
CN114781713B (zh) 基于卫星反演降水产品的雨量站网优化方法及装置

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