基于HDBSCAN算法的岩体结构面产状识别及分组方法

2022-07-08王述红

东北大学学报（自然科学版） 2022年6期

王述红, 魏崴, 陈浩, 尹宏

(东北大学资源与土木工程学院, 辽宁沈阳 110819)

岩体是一种整体表现出不连续、非均质及各向异性等特征的复杂地质结构,主要由岩块和结构面两部分组成.一般认为,岩块由以节理、断层、软弱夹层为主要形式的结构面切割形成,内部物质相对连续密实,大多数情况下不会发生破坏;而结构面作为岩体内部力学属性最薄弱的地质界面,其空间姿态、尺寸、分布情况等因素在很大程度上决定了岩体的力学性质和稳定性状态[1].

现阶段,实际工程中常用赤平投影法[2]等方法对岩体稳定性进行初步评估,结构面的倾向、倾角等产状信息决定了投影线在赤平面上的位置和组合关系,对判断可能失稳岩体的滑动形式、规模、形态等产生直接影响.主流的岩体稳定性模拟方法诸如离散元法[3](DEM)、数值流形方法[4](NMM)、非连续变形分析[5](DDA)等,为了更好地模拟工程岩体由于受到非连续界面的切割而表现出的不连续、非均质及各向异性等力学特征,往往需在建模过程中对岩体结构面进行精细表征.一些块体稳定性分析软件例如GeoSMA-3D[6-7]、KBTA[8-9]等同样需要准确的结构面产状信息来进行地质统计分析,以期获得更加真实的结构面随机分布规律和统计函数,在此基础上进一步通过计算机模拟技术生成确定和随机结构面的裂隙网络.因此,如何准确、快速、高效地获取结构面产状信息,是开展岩体稳定性分析工作亟待解决的首要问题.

传统上多采用人工测量结构面产状的方法,借助罗盘、测绳、皮尺等工具在现场对工程岩体进行接触测量.该方法效率低、工作量大,并且存在着较高的安全隐患,无法满足现阶段对大规模结构面信息采集的需要.同时,由于测量空间受限,对于位置较高、跨度较大的岩体只能凭借经验和目测推断结构面产状,存在较强的主观性.随着测量技术的发展和进步,各种新型的非接触测量手段,例如数字近景摄影测量技术、三维激光扫描技术、无人机测量技术等已被成功应用于岩土工程的多个领域[10-12].近年来,国内外学者依托上述测量手段构建岩体的真实三维点云模型,进一步探索岩体结构面产状提取和表征的新方法[13-16].

为了更加高效、准确地获取结构面产状信息,本文采用无人机倾斜摄影技术构建研究区域的三维点云模型.针对基本DBSCAN(density-based spatial clustering of applications with noise)算法在处理点云数据时聚类指标不够全面的情况,提出了基于点云附加属性扩展判据的HDBSCAN(hybrid density-based spatial clustering of applications with noise)算法用于岩体结构面点云的识别和分组.同时,对于非全裸露的自然岩体(表面覆土、覆雪和覆盖植被等),基于点云的颜色及密度属性对岩体表面的非结构面点云进行降噪处理,在较大限度保留原始数据的情况下避免非目标点云对结构面产状提取的影响.最后,将上述方法应用于云南玉溪大石洞灰岩矿某一典型覆土岩坡结构面产状的提取和分组工作,用于验证本文所提方法的有效性及可靠性.

1 点云法向量计算

图1 点云法向量计算Fig.1 The normal vector calculationof the point cloud(a)—点云局部平面拟合； (b)—法向量方向调整； (c)—原始点云数据； (d)—点云法向量求解.

Apix+Bpiy+Cpiz+Dpi=0 .

(1)

可以得到点pi的协方差矩阵为

(2)

求解Σi矩阵的特征值及其对应的特征向量:

Σivj=λjvj(j=0,1,2).

(3)

式中λj为vj所对应的特征值.设λ0<λ1<λ2,则v0即为该测点k邻近平面的法向量npi=(Api,Bpi,Cpi),即点pi的法向量.通过文献[18]方法调整法向量指向物体表面外部.

2 基于HDBSCAN算法的结构面提取

2.1 改进DBSCAN算法

DBSCAN算法是一种基于密度的计算机迭代算法[19].该算法无需提前设置聚类簇的数目,即可拟合任意几何形状的点云簇,并且有效识别并剔除异常数据,多被应用于结构面点云簇的提取和分割[12-13,15].

算法的核心思想是借助Eps和Minpts两个输入参数(Eps表示当前数据点邻域的最小范围,Minpts表示当前邻域内点的总数量)构建单位超球区域,通过密度直达、密度可达等判定剔除背景点,并不断将核心点和边界点合成为新簇,最终实现求解当前数据集密度相连点最大集合的目的.简单来说,Eps参数用于表征数据点之间的距离,Minpts用于描述Eps范围内数据的“密度”,这使得DBSCAN算法存在较强的参数敏感性.对于密度分布不均匀、聚类间距差异相对较大数据集,聚类效果较大程度上取决于两个参数的联合调参,全局单一的参数设置容易导致过分割和欠分割情况的发生,即将本应划分为同一簇的数据点划分为多簇或将其识别为噪声点剔除.

需要注意,传统的DBSCAN算法常应用于二维或三维点数据集的聚类,因此采用欧氏距离度量各数据点之间的距离,聚类的评价指标相对单一.与普通空间离散点不同,点云不仅具备最基本的三维坐标信息,而且还具备能反映真实物体曲面在某点处的RGB颜色、局部形状(法向量、曲率)等重要属性.忽略这些关键信息,仅使用欧氏距离判据对点云数据进行聚类分析显然不具备充分的合理性与科学性.针对这一问题,本文尝试在聚类过程中引入“候选点”概念,并且在原算法密度聚类的基础上混合加入方向相似性判据、HSI颜色距离判据及马氏距离判据,旨在提高算法对点云数据聚类的鲁棒性和准确性.

算法改进的主要思路为:1)将密度聚类过程中识别的背景点或低于预设数量的点云簇标记为候选点;2)搜索各个候选点与其距离最近的核心点,分别计算两者间的HSI颜色距离、方向相关系数及马氏距离;3)将同时满足以上三种判据的候选点合并至其最近的核心点所属的簇,将最后的无分簇的数据点作为噪声点剔除.

如图2所示HDBSCAN算法的基本步骤,将包含空间坐标信息、RGB颜色属性及点云法向量属性的点云集作为输入数据,采用HDBSCAN算法对其进行聚类和分组.首先进行基本的密度聚类,此时的点云数据点可被划分为如图所示的核心点、边界点及背景点,同时该数据集被划分为红、蓝两部分点云簇.为了方便介绍原理,称红、蓝两部分点云簇分别为密度分簇1和密度分簇2,并且假设密度分簇2中的数据点数量低于预设数量,则图中的背景点及密度分簇2中的数据点作为候选点继续进行扩展判据的聚类纠偏与点云簇的合并步骤,直至遍历每一个点云簇,聚类结束并将最后无分组的数据点作为噪声点进行剔除.应说明的是,由使用者自定义设置的点云簇数量阈值并非是候选点判别的唯一准则,在实际的聚类过程中,可以根据点云密度及聚类效果进行调整,进一步确保算法的有效性和灵活性.

图2 HDBSCAN算法基本步骤Fig.2 Basic calculation steps of the HDBSCAN algorithm

设样本集中任意一个点云为pi=(xi,yi,zi),该点处的法向量为ni=(ni,x,ni,y,ni,z),其中元素分别表示法向量在X、Y、Z坐标轴上的分量,所包含的RGB颜色属性向量为Ci=(Ri,Gi,Bi),分别对应R(红)、G(绿)、B(蓝)三个颜色通道的亮度值.

1) 方向相似性判据.本文以法向量夹角θ作为衡量两点云方向相似性的依据.对于任意两点云pi和pj,其方向相似性可用法向量夹角θi,j的大小来衡量[20],具体表达式如下:

(4)

式中,ni及nj已做无量纲处理,且各元素均为正值.此时,θi,j取值范围为[0,π/2],取值越小说明两法向量越接近,即两点云的方向越接近.对于给定的夹角阈值ηθ,若θi,j<ηθ,则称点pi和点pj具有方向相似性.

2) HSI颜色距离判据.值得注意的是,RGB颜色坐标系统不易直观地呈现颜色的感知属性,两颜色之间的差异程度不能简单地用颜色空间内对应两色点欧氏距离的远近来评判.为了更好地表征两点云之间的颜色相似性,本文将点云原始色彩信息由RGB色彩空间转换到HSI色彩空间(图3).由于HSI的三个分量中表示色调的H(hue)分量对色彩的描述能力与人眼的感知能力最为接近,为了简化运算和提高计算效率,不考虑饱和度(saturation)和亮度(intensity)分量的作用,仅通过H分量的差异描述两点云之间颜色的相似程度.

图3 RGB、HSI色彩空间Fig.3 The RGB and HSI color space

将点云颜色由RGB格式转换为HSI格式,计算公式如下:

(5)

在HSI色彩空间中,两种色彩Ci,Cj之间的度量定义为

DH(Ci,j)=[(Ii+Ij)2+(Si(cosHi-SjcotHj)2+
(SisinHi-SjsinHj)2]1/2.

(6)

(7)

考虑到岩体的同一表面本身可能存在一定的颜色差异,本文适当放宽颜色距离阈值ηDH,取ηDH=0.15.若DH(Ci,j)<ηDH,则称点pi和点pj满足颜色相似性条件.

3) 马氏距离判据.马氏距离由Majalanobis提出,用于描述数据之间的协方差距离.由于考虑了数据间各种特性的关联性,对点云数据的聚类效果较好.本文选取马氏距离计算点pi和点pj之间的绝对距离,具体表达式如下:

(8)

DM(pi,pj)<ηDM.

(9)

式中:Σ为点pi和点pj的协方差矩阵;ηDM为马氏距离判据,本文取ηDM=1.5 Eps.若满足式(9),则称点pi和点pj满足马氏距离判据.

2.2 点云结构面产状计算及分组

进一步对HDBSCAN算法聚类分割得到的任意结构面点云簇Si=(p1,,p2,…,pi)i≥3的产状进行提取.与求解点云法向量的过程类似,计算协方差矩阵和特征值分解步骤不再赘述.通过最小二乘法得到该结构面点云簇的最佳拟合平面:

Aix+Biy+Ciz+Di=0 .

(10)

该平面的单位法向量为

(11)

则该结构面的倾向即可表示为

(12)

倾角可以表示为

(13)

同样,以任意两结构面法向量所夹的锐角作为结构面产状优势分组的依据,通过设置结构面法向夹角阈值对提取的结构面产状进行优势分组.但与点云方向的相似性判断有所不同,结构面分组需将倾向相差约180°且倾角较陡的结构面合并划分为同一组类.为了避免这一特殊情况造成的分组失误,本文采用结构面点云簇法向量夹角正弦值的平方来度量结构面优势产状的相似性,具体表达式为

(14)

式中:γ为Ni和Nj法向量所夹锐角；Ni和Nj分别为结构面点云簇Si及Sj的法向量.若γ<ηγ,则认为两结构面隶属于同一优势分组.

图4 HDBSCAN算法流程Fig.4 Flow chart of HDBSCAN algorithm

2.3 算例验证

选取文献[21]中的标准几何体数据集中的二十面体和立方体点云数据验证本文所提改进HDBSCAN算法的有效性.该数据集由瑞士洛桑大学使用三维扫描仪(Konica Minolta Vivid 9i)采集,所得点云数据来自不同角度的扫描数据的重叠和累加.该设备扫描视线与顶点的倾斜角度约为301°,不可避免地导致部分点云数据存在密度分布不均匀的情况发生.由于多面体本身不具备“产状分组”的性质和需求,仅对算法的聚类效果进行展示.这里取DBSCAN算法及HDBSCAN算法的基本参数Minpts=4,Eps设置为文献[22]推荐值,即每个点云与其最邻近4个点距离矩阵的平均值与两倍标准差的和,其余参数取为默认参数值.

详细的聚类结果及可视化效果如图5所示.其中,图5a,5e为原始的二十面体及立方体点云数据(无颜色属性),图5b,5f为基本DBSCAN算法的聚类结果,图中不同颜色代表不同的点云簇族.可以清楚地看到,对于密度分布不均匀的点云数据,直接采用固定参数的DBSCAN算法容易导致明显的欠分割和过分割的情况发生.对于二十面体来说,欠分割的情况较为严重,算法将大部分三角面归为一类,分割效果相对较差.而立方体由于上表面点云密度较大,文献[22]的参数设置仅对上表面的适应效果较好,对侧面的适应性较差,导致较为严重的过分割现象发生.将各个数据集拟聚类的平面用不同的颜色区分,用于模拟真实点云数据集的颜色属性如图5c,5g所示,采用本文提出的HDBSCAN算法进行聚类分析,聚类效果如图5d,5h所示.可以看到,由于在基本密度聚类的基础上引入点云颜色、方向和马氏距离扩展判据,本文算法能够准确地分割出点云数据集的平面,并且能在一定程度上避免DBSCAN算法的参数敏感性问题,聚类效果令人满意.

图5 标准多面体聚类结果Fig.5 Standard polyhedral clustering results(a)，(e)—二十面体及立方体初始点云数据; (b)，(f)—DBSCAN聚类结果; (c),(g)—点云颜色划分; (d),(h)—HDBSCAN聚类结果.

2.4 工程实例

本文依托于云南省玉溪市大石洞灰岩矿项目,以实际矿山岩坡点云数据作为研究对象,用于验证本文所提HDBSCAN算法的有效性.该矿区位于距玉溪市中心城区直平距约为6.3 km的玉溪盆地西侧边缘山区,地处低纬度高原地带,属亚热带半湿润高原季风气候,5月至10月为雨季.选取矿山北帮第四台阶某一典型岩坡作为研究区域,该区域岩体以石灰岩为主,出露结构面及裂隙发育,岩体破碎程度较高.岩体左侧坡脚处坡面相对平整,右侧坡面岩体破碎程度较大,揭露情况较为复杂,存在较多随机结构面.由于现场考察时间正值雨季,受到雨水裹挟和开挖扰动作用,坡顶岩土体存在一定的表层滑移,导致岩体表面呈现红褐色,并伴随泥土和碎石滑落,传统人工接触测量的方法存在较大的安全隐患.

本文采用无人机近景测量的方法建立真实边坡的三维点云模型(图6),主要步骤为:1)SFM算法[23]检测、匹配无人机图像特征生成稀疏点云;2)光束法平差[24]精校稀疏点云在物方坐标系中的坐标参数；3)多视觉立体视觉算法PMVS[25]三维重建,生成岩体表面稠密点云.该区域共获取点云数据220 241个,测量区域内结构面点云未出现明显缺损的情况.岩体表面存在覆土现象,其中坡顶和坡脚处覆土和岩土体堆积情况相对严重,中部岩体表面轻微覆土.为避免对后续结构面的识别和提取工作产生影响,本文基于点云的颜色及密度属性对边坡非岩体结构面点云数据进行过滤.考虑到岩体表面轻微覆土对结构面产状提取的影响不大,为了最大限度地保留结构面点云信息,仅针对岩体表面覆土情况较为严重的区域(坡顶覆土区域及坡底岩土体堆积区域)进行降噪处理,并进一步通过剔除离群和非结构面点云等操作,得到最终的结构面识别区域如图7所示.

图6 研究区域及三维点云模型建立Fig.6 The studied area and the 3D point cloud modeling(a)—矿山地理位置及全貌; (b)—无人机三维点云建模; (c)—研究区域实拍.

图7 非结构面点云及噪点数据剔除Fig.7 Unstructured surface cloud and noise data elimination(a)—原始边坡点云数据; (b)—点云颜色及密度属性叠加云图; (c)—非结构面点云数据剔除; (d)—结构面识别区域.

应用本文所提出的HDBSCAN算法对研究区域内结构面进行聚类分析,本文取DBSCAN算法基本参数Minpts=4,Eps设置为文献[22]推荐值.其余参数分别设置为:θi,j=5°,ηDH=0.15,ηDM=1.5 Eps,ηγ=3°.结果如表1和图8所示.与人工测量结构面产状的方式相比,改进后的DBSCAN算法得益于点云附加属性判据的补充和扩展,不仅可以准确地分割出大块完整的结构面,而且对于破碎结构面也有较为准确的识别效果,能够较大程度地还原真实岩体结构面的几何形貌.所提取结构面倾向、倾角产状与人工测量方法吻合度较高,最大相对误差不超过0.59%,平均绝对误差分别为0.32%和0.37%,表明该方法具有较高的可靠性,能够满足实际的工程需求.

表1 部分结构面产状识别和提取结果Table 1 Partial identification and extraction results of the structural plane strike

图8 HDBSCAN点云结构面聚类结果Fig.8 HDBSCAN point clouds clustering of the structure planes(a)—结构面聚类结果; (b)—结构面产状密度等值线云图.

进一步基于产状信息对结构面产状进行优势分组.首先对点云簇识别的结构面产状进行统计处理,得到结构面产状的极点和等密度云图.从图8b可以看到,结构面产状等密度图存在多个局部最小值,各极点产状的相似性较高,优势产状的分组界限相对模糊,仅从极点及等密度云图难以准确地判断结构面的平均产状.

采用2.2节所述方法对结构面进行优势分组,将不同分组的结构面用相同的颜色进行表示,分组结果如图9所示.研究区域内的结构面共划分为J1,J2,J3,J4及J5五组优势产状(倾向/倾角分别为:263.64°/53.38°,90.79°/72.78°,125.89°/76.87°,55.71°/87.12°,315.43°/33.92°).由图9a 和9b可以直观清楚地看到结构面的分组情况,分组结果也较为符合密度等值线的分布规律,并且能够将倾向产状相差约180°的结构面准确地划分为同一组别.