基于K-means聚类方法和I Index聚类有效性检验指标的岩体结构面自动分组及应用

，，

(1.黄河勘测规划设计有限公司 地质工程院,郑州 450003；2.中国地质大学(北京) 水资源与环境学院,北京 100083)

1 研究背景

岩体由岩块和分割它们的结构面组成。结构面是发育于岩体中，具有一定方向、延伸性和一定厚度的各种地质界面，对岩体的强度、变形和渗透性有着决定性影响[1]。在野外，常采用地表露头测量、探硐和钻孔岩芯裂隙编录等方法对结构面数据进行采集。工程地质、水文地质工作中，一个基础但十分重要的环节是对结构面产状数据进行统计分组研究，并确定每组结构面的优势产状。常用的方法主要是绘制倾向和走向玫瑰花图以及在吴氏、施密特网上作极点等密度图[2]。这些方法操作简单，结果直观，易于实现，但是由于地质工作者的专业素养不同，对划分结果的主观判断难免存在差异,当结构面的产状实测数据数量多且界线较为模糊时更是如此。因此，应借助数学手段对结构面的产状数据进行客观划分。

岩体结构面产状数据的计算机自动分组一直是国内外工程地质、水文地质界研究的热点问题,国外早在20世纪70年代便已经开始了相关研究。Shanley和Mahtab[3]首次提出了一套完整的岩体结构面产状数据自动划分的算法及其数值实现过程：通过对单位球面上的投影点进行密度分析完成优势组划分。Hammah和Curran[4]使用模糊C均值聚类方法进行了结构面优势组的自动划分。Klose等[5]基于随机逼近法，提出了一种无监督的机器自学习法，快速完成优势组划分。Jimenez-Rodriguez和Sitar[6]将谱分析方法和K-means聚类法相结合进行结构面优势组划分。Tokhmechi等[7]讨论了采用结构面的多个参数进行分组的重要性。我国在20世纪80年代也已开始了相关研究，并达到了国际领先水平[2, 8-9]。近些年来，国内一些学者基于不同数学原理，先后提出或改进了结构面自动划分的算法，并在实际工程中进行了应用，获得了不错的效果[10-14]。

本文以岩体结构面法向向量为样本，提出了基于K-means聚类方法和IIndex聚类有效性检验指标的岩体结构面自动分组方法，并开发了程序RDAP，服务于岩体变形、破坏和渗透性的分析与计算。

2 岩体结构面自动分组方法

2.1 结构面法向向量

野外测量得到的结构面产状数据通常使用倾向α和倾角β进行表示。根据空间解析几何的相关知识，基于右手法则，可将产状数据转化为单位向量[2],即

(1)

式中nx，ny，nz分别为结构面法向向量在x，y，z轴上的投影。

2.2 K-means聚类方法

K-means聚类方法是一种基于距离的聚类算法，采用距离作为相似性的评价指标，认为距离靠近的样本组成聚类组，聚类目标是使各聚类组总的距离平方和最小[15]。K-means算法的基本思想是初始随机给定K个聚类组中心，并根据距离按照最邻近原则把待分类的样本分配给各个聚类组，然后重新计算各个聚类组的中心，一直迭代，直到聚类组中心的移动距离小于某个给定值。相较于其他聚类方法，K-means聚类方法具有快速、简单的特点，适合大规模数据集的挖掘。

本文对传统K-means聚类方法稍加改进，沿袭田开铭和万力的方法[2]，使用结构面法向向量与聚类组最大概率方向(最大概率方向等价于聚类组结构面法向向量合矢量的方向)的夹角定义聚类距离，而不是K-means聚类方法通常使用的欧几里德距离。这样的改进突出了结构面产状数据的矢量性，具有更加显著的工程、水文地质意义。改进后的K-means聚类方法包括以下几个步骤：

(1)根据预定义聚类组数，随机指定若干聚类组中心，计算每个结构面产状样本数据与聚类组中心之间的角度。

(2)将每个结构面产状样本数据分配给与它夹角最小的聚类组中心，形成新的聚类组。

(3)将各组结构面产状样本数据的最大概率方向作为新的聚类组中心。

(4)重复步骤(2)和步骤(3)，直至所有聚类组的中心固定，结构面产状样本数据分配固定。

2.3 聚类有效性检验指标

K-means聚类方法均属于动态聚类法。动态聚类法通常无法直接确定聚类优势组数。实际应用中，除了人为根据经验预先设定聚类组数外，常根据聚类有效性指标确定最优聚类组数。聚类有效性指标较为常用的有：Davies-Bouldin (DB) Index，Dunn’s Index，Calinski Harabasz(CH) Index，Xie Beni(XB) Index[13]。近期，Maulik和Bandyopadhyay[16]提出了一种新型聚类有效性检验指标Iindex,即

(2)

(3)

式中：N为结构面聚类组数；v表示属于第i聚类组的结构面的法向单位向量；cNi表示第i聚类组结构面中心向量，即最大概率方向;d(v,cNi)表示v和cNi两向量之间的欧几里得距离。

IIndex聚类有效性检验指标使用各聚类组中心距离的最大值度量组间分离度，使用聚类组中各成员与组中心的距离和度量组内分离度。IIndex聚类有效性检验指标认为：IIndex取值愈大，其所对应的分组效果愈佳。此外，实例验证，当聚类距离定义为向量夹角时，相较于上述常用的聚类有效性检验指标，使用IIndex聚类有效性检验指标选择优势聚类组数的效果更好[13]。

2.4 岩体结构面自动分组程序RDAP

RDAP是岩体结构面自动分组程序英文名称Rock Discontinuities Auto-classification Program的简称。RDAP使用MatLab科学计算语言编制计算程序和GUI界面。原理上，RDAP首先根据Iindex进行聚类有效性检验，选取Iindex最大值所对应的聚类组数作为结构面优势组数，然后将其作为预定义聚类组数代入改进的K-means聚类算法实现岩体结构面产状数据的计算机批量自动分组(图1)。RDAP的自动分组结果可作为Dips,Unwedge,Steronet等商业软件的输入数据，服务于岩体变形、破坏和渗透性的分析与计算。

图1 I index聚类有效性检验指标随聚类优势组数的变化Fig.1 RDAP interface showing the plot of cluster validity index I versus cluster numbers

3 算法验证

Shanley和Mahtab[3]于1976年首次提出了一种结构面优势组自动划分的方法及其数值实现过程，并以算例的形式给出了286条结构面数据及自动分组结果。学术界对这组数据的分组形式已经形成了较为统一的认识，并常将其作为检验其他分组方法正确与否的试金石[5]。Shanley和Mahtab使用其算法，将这286条结构面数据分为了3个聚类组(图2)。

图2 Shanley和Mahtab[3] 分组结果与RDAP分组结果的对比Fig.2 Classing results by Shanley and Mahtab[3] and by RDAP

Iindex聚类有效性检验指标随聚类组数的变化图表明：Iindex聚类有效性检验指标在聚类组数为3时等于1.72，为最大值(图1)。因此，根据本文所提分组方法，286条结构面数据最优可以分为3个聚类组，这与Shanley和Mahtab所识别的最优聚类组数目一致。具体来讲，使用RDAP，第1组识别得到了114条结构面，第2组59条结构面，第3组113条结构面(图1)。为了直观地进行比较，我们将Shanley和Mahtab的分组结果和RDAP的分组结果分别绘制在极点图上(图2)。经比较，两者的分组结果基本一致，证明了本文所提分组方法的可靠性。

需要强调的是，除了可以给出分组结果，RDAP还能根据聚类组的最大概率方向给出每组结构面的优势、代表性产状(图1)，便于进行后续工程地质稳定性分析、水文地质的计算和防渗灌浆钻孔优势方位的选定。

4 实例应用

某地下工程位于广东江门市附近，坐落于花岗岩岩体内部。工程地下建筑物主要为斜井、竖井及实验厅3部分。目前，工程面临的主要问题是施工涌水。

4.1 涌水裂隙分组

竖井在施工过程中观测并记录了42条主要涌水裂隙(表1)。将这42条涌水裂隙的倾向和倾角数据作为输入数据，选择“程序自动分组方法”。RDAP将这42条主要涌水裂隙自动识别为5个裂隙组(表1)，并可给出这5个裂隙组的最大概率方向，即优势产状(表2)。由于第2和第3组涌水裂隙数据少(只有2条)，代表性差，可将其舍去，认为花岗岩岩体中主要发育3个裂隙组。在自动分组结果的基础上，根据组内每条涌水裂隙的出露高程(剔除异常点)、隙宽描述和每组裂隙的优势产状，可以计算得到每组裂隙的真张开度和真裂隙间距(表2)。

4.2 灌浆钻孔最佳方位的选取

涌水量的预测与防治是该工程施工中的重点与难点。灌浆防渗处理技术是工程涌水防治的有效、常规方法。其中，灌浆钻孔最佳方位的选取是影响灌浆效果的关键因素之一。

假设灌浆钻孔的倾向为a，倾角为b(与水平面之间的夹角)，灌浆钻孔的空间向量可以表示为

(4)

式中m，n，p为灌浆钻孔在x，y，z轴上的投影。

那么，灌浆钻孔揭露的视裂隙间距与真裂隙间距的关系可以表示为

(5)

表1 竖井涌水裂隙基本情况Table 1 Basic information of gushing fractures

表2 涌水裂隙统计分组结果Table 2 Classing result of gushing fractures

式中：sa为视裂隙间距；s为真裂隙间距；θ为灌浆钻孔与裂隙面法向向量之间的夹角。

灌浆钻孔的最佳方位是在同等压力下吃浆量最大的方位。换言之，应是灌浆钻孔穿过最多、最宽裂隙的方位[17]。若忽略裂隙宽度、粗糙度、充填情况等因素对灌浆效果的影响，灌浆钻孔最佳方位的选取等价于单位长度的灌浆钻孔尽可能多地与各组裂隙相交，即获得最大的灌浆钻孔综合视裂隙密度的方位。转化成数学语言，即为求解以下目标函数的最大值，其表达式为

(6)

图3 综合视裂隙密度三维渲染图Fig.3 Three-dimensional plot of the apparent fracture density

5 结 语

(1)以岩体结构面法向向量为样本，提出了基于K-means聚类方法和IIndex聚类有效性检验指标的岩体结构面自动分组方法，并开发了岩体结构面自动分组程序RDAP。

(2)通过与经典文献进行对比，本文所提结构面自动分组方法的可靠性得到了论证。

(3)使用本文提出的结构面自动分组方法对某地下工程的实测涌水裂隙资料进行了预处理，初步计算了灌浆钻孔的最佳方位，为工程涌水的防治提供了依据。