0 引言

岩体结构控制理论认为，岩体的稳定性、物理力学性质及变形破坏机制均受到其结构的影响^[1]。在斜坡稳定性方面，岩体结构，尤其是软弱结构面，对基岩斜坡的变形与破坏具有显著的控制作用^[2]。众多学者根据多个岩体质量指标对不同岩体结构控制下的岩体进行了质量分级，方法包括BQ法^[3]、RMR法^[4]、水利水电HC分级法等。获取岩体结构面的几何特征参数(如产状、间距、迹长等)是岩体质量分级的关键环节。因此，结构面的几何特征参数对于斜坡岩体特性的评估、岩体质量与失稳机制的研究具有重要意义。

岩体结构特征的准确识别与测量是计算斜坡稳定性和预测岩体斜坡灾害的前提。斜坡岩体结构的调查和测量方法越来越趋向于远程化和智能化。传统手工测量方法效率低，难以应对高陡斜坡的测量需求。低空摄影测量和三维激光扫描等技术的发展，使得快速、高效获取高精度三维岩体结构数据成为可能^[5-7]。无人机摄影测量技术具有天然的优越性，不受环境限制，能够精确获取高陡斜坡的三维模型数据^[8-9]。Wang等^[10]使用无人机摄影技术对雅鲁藏布江断裂周边的斜坡岩体结构进行了识别与测量。刘宽富等^[11]利用无人机摄影技术建立了边坡的三维实景模型，以进行FLAC数值计算。Herrero等^[12]利用无人机航拍的照片和结构从动算法成功建立了西班牙某采石场的三维模型，以低廉、迅速、简便的手法获取岩体结构特征。Kong等^[13]结合无人机摄影测量与运动恢复结构(SfM)技术，建立了高分辨率数字露头模型。Caliò等^[14]利用无人机摄影与红外热成像技术，构建了数字化岩石悬崖模型，以研究意大利马里内洛湖自然保护区内不稳定悬崖的地质构造与形态特征。

获得斜坡岩体的三维结构模型后，对其数字化信息的提取与分析至关重要。在这方面的研究中，前人进行了大量探索和应用。Riquelme等^[15]采用K最近邻算法和核密度估计法，开发了半自动提取岩体不连续结构面的Matlab程序，进行优势结构面的聚类与提取。李华明等^[16]提出K-mean++的三维产状自动分组算法，实现了岩体结构面信息的自动提取。李水清等^[17]借助无人机摄影测量与三点法对结构面产状进行了半自动提取，并对优势结构面进行分类与参数获取，然而此过程仍需手动选定结构面的特征点。王俊杰等^[18]、Li等^[19]采用K-means聚类算法进行优势结构面的划分，准确获取了结构面产状。梁玉飞等^[20]运用色彩重建技术(HSV)和基于密度的带噪声空间聚类(DBSCAN)算法对滑坡区岩体的结构面进行分组与识别。Guo等^[21]、刘昌军等^[22]利用模糊C均值算法对岩体结构面进行优势结构面分类。

前人在研究中发现，与K-means等硬聚类算法相比，模糊C均值算法因引入隶属度指数，在优势结构面的聚类分组中表现出更佳的效果，被广泛认为是岩体结构面方向最好的数据聚类分析算法之一^[23]。然而，模糊C均值算法对噪声和离群点敏感，其关键参数(即模糊度值m)的取值直接影响聚类结果。当前该参数的取值多依赖经验^[24]，难以根据迭代过程动态更新，从而制约了聚类结果的精细化与准确性。在获取结构面产状后，仍需提取迹长等几何数据。当前，这一过程主要依赖手动测量或基于点云数据的提取^[25]。然而，结构面表面风化造成的凹凸差异、坡体表面植被干扰及数据获取过程中的不确定性等因素，导致许多非结构面点云数据难以被有效识别和剔除，从而影响几何数据提取的准确性，进而直接影响岩体稳定性预判的精度。

本研究在前人研究的基础上，改进模糊C均值算法，引入了新的参数以动态调整模糊度值，避免算法陷入局部最优解，增强其鲁棒性，并优化迭代过程，从而提升聚类效果。以四川省G108国道附近某露天边坡为案例，通过随机霍夫变换算法和改进的模糊C均值算法进行点云的聚类分析，基于Fast Marching算法提取闭合结构面及其产状，并使用Matlab编程结合ArcGIS软件，开展各组结构面间距和迹长的自动测量与统计研究。

1 研究方法

本研究利用无人机摄影测量技术，实现斜坡岩体结构面的自动识别与信息提取。首先，通过无人机快速获取斜坡岩体的三维点云模型，并采用随机霍夫变换算法计算点云数据的法向量。接着，利用改进的模糊C均值算法获得优势结构面分组。然后，基于优势结构面分组，使用Fast Marching算法识别单个闭合结构面并计算其产状。最后，结合Matlab编程与ArcGIS软件，对识别出的闭合结构面进行间距和迹长的自动计算，从而实现对斜坡岩体结构面信息的自动提取。

1.1 法向量计算

传统上，岩体三维模型的重建依赖于三角剖分原则，通过点云数据生成TIN网格，但该过程整体效率偏低。为了解决这个问题，本研究采用聚类分析技术直接从结构面点云数据中提取信息。获取点云数据的法向量信息是提取斜坡岩体结构特征的前提。法向量信息可用于研究岩体结构的空间倾向和倾角等特性，从而优化岩体结构几何特征信息的提取速率。

法向量信息的计算主要包括回归法、Delaunay三角分割法和邻域分割法这3种方法^[26-27]。回归法通常不易受数据噪声影响，但在存在异常值时，可能无法准确估计点云边缘的法向量；Delaunay三角分割法能够有效地保持点云边缘的尖锐特征，但在噪声条件下表现不佳；邻域分割法将点云分为若干邻域，再对每个邻域计算法向量。这3种方法在去噪、运算效率和精确度方面各有优劣，难以做到全面平衡。本研究采用Boulch等^[28]提出的随机霍夫变换算法计算法向量。该方法具有稳健的法向量评估能力，能够有效处理噪声、异常点及棱角锐利位置等情况，同时在精度、速度与鲁棒性之间实现了良好的平衡^[29]。

1.2 基于改进的模糊C均值算法进行结构面聚类

在获得结构面点云数据的法向量信息后，需进行聚类分析以提取岩体的优势结构面。本研究采用改进的模糊C均值算法来实现对结构面点云数据的聚类分组。作为一种常见的聚类算法，模糊C均值算法与其他聚类算法的不同之处在于它以隶属度作为衡量指标，为每个样本和每个聚类赋予权重，表示样本属于该聚类的程度，使其可以在[0, 1]内连续取值，从而提高了聚类结果的精准度。用U表示点云数据的隶属度矩阵，其表达如下：

(1)

(2)

(3)

(4)

式中，n为点云法向量数据的总数，即样本个数；c为聚类中心的总数，即选定的聚类组数；u_ij为第i个点云属于第j个聚类中心的隶属度，取值范围为[0, 1]；m为控制分类模糊度的参数，通常为大于1的实数；d_ij和d_ik为欧几里得距离，用于度量第i个样本x_i与聚类中心j或k之间的距离；v_j为第j个类别的聚类中心；x_i为第i个样本数据。

通过式(3)和式(4)反复修改聚类中心和数据隶属度并进行分类。当算法收敛时，理论上就得到了各类聚类中心及每个样本对于各个聚类的隶属度，从而完成模糊聚类的划分。模糊C均值聚类算法本质上就是寻找最佳分类，以使得该分类能够产生最小的目标函数值J_b。

(5)

(6)

式中，U为点云数据的隶属度矩阵；V为各类别的聚类中心组成的矩阵; v_i为第i个聚类中心。

模糊C均值算法对噪声和离群点敏感，且当模糊度参数m(以下简称为“模糊度”)的取值不当时，会影响最终结构面信息的提取效果。模糊度的取值对算法的效果影响很大。当m = 1时，模糊C均值算法退化为硬聚类方法。然而，限于取值方法的缺失，目前对于模糊度的取值仍主要依赖经验，如取m=2或3。当模糊度取较大值时，可以提高算法对数据的容忍度，提升算法的鲁棒性。对于一些数据分布模糊或者难以明确定量划分的情况(如边缘数据或存在异常值)，较大的模糊度有助于更合理地对数据进行分类，对那些不是明显属于某个类别的数据进行更加充分的考虑。相反，对于特征明显、易于划分的数据，过大的模糊度可能导致聚类结果模糊不清，簇的边界不够明确。模糊度取值较小时，在面对数据分布模糊或者难以明确定量划分的情况时，可能导致算法过早收敛至局部极小值，影响聚类结果的准确性。因此，合理选择模糊度对聚类结果至关重要。

鉴于不同情境下数据集的模糊性和类别显著性差异，需灵活调整模糊度以提高计算合理性。因此，本研究提出了一个新的参数α，用于控制模糊度m的更新，通过平滑调整的方式动态调整算法的模糊度，使其在连续范围内进行调整，避免突变或不连续，从而增强聚类结果的稳定性和连贯性。通过判断数据集的类别显著性，确定是否需要增大或减小模糊度。在算法迭代过程中，不断变化的模糊度通过式(3)不断调整隶属度值的计算，进而通过式(4)影响聚类中心的更新，优化算法迭代过程获得最佳结果。通过这种方式，改进模糊C均值算法，提升其稳定性和适应性。

根据隶属度矩阵U的分布情况动态更新模糊度。通过计算隶属度最大值和最小值的平方差，得到一个指标ex，以表征数据集的模糊性质及类别特性。

(7)

式中ex为根据隶属度矩阵计算出的指标，用于表征隶属度的模糊程度。

当指标ex超出预设阈值α时，表示类别明显，此时应降低模糊度m的值以增强隶属度的模糊性；若ex≤α，则表明类别不明显，需相应提高模糊度m的值以降低隶属度的模糊性。

(8)

通过这种动态调整模糊度的方式，使改进的模糊C均值算法能够更精细地调节隶属度的模糊性，从而优化算法的迭代过程，更好地适应各类数据集的特性和聚类任务的需求。利用改进的模糊C均值算法对结构面点云数据进行聚类分析后，可对各组优势结构面的点云数据进行结构面识别与信息提取。

1.3 结构面的识别与产状计算

在完成点云数据的聚类并获得优势结构面分组后，需要进行结构面的识别。其本质是根据共面度将点云划分为相邻点簇。本研究采用Cloud Compare软件的Fast Marching算法^[30]对分类的优势结构面点云数据进行封闭成面处理，完成结构面的识别与产状计算。该算法使用八叉树结构描述三维空间的树状数据要素，其主要流程是将初始空间用八叉树结构定义的规则网格进行细化，然后根据八叉树中相邻体元的均方根距离是否超越预设的最大距离来判定是否进行合并操作。如果相邻体元的均方根距离不大于预设值，则合并，并进行迭代运算。在识别过程中，需要根据结构面的大小和噪声因素确定适宜的参数，以实现最佳拟合效果，最终获得表面平滑无起伏的闭合结构面。通过这种方式提取的结构面较好地重现了实际结构面的几何形状与空间规律，避免了点云起伏及植被噪点的影响^[31]。

完成结构面的识别后，需要计算其产状。此过程同样由Fast Marching算法实现，可同时完成结构面的识别与产状的自动提取。该算法识别并提取结构面后，能够利用结构面的法向量信息n(n_x, n_y, n_z)进行产状的计算。所采用坐标系的z轴为竖直方向，与z轴垂直的xoy平面坐标轴的示意图如图 1所示，其中倾向的计算以北方向，即y轴正方向为起始点，顺时针旋转。产状计算公式如下：

(9)

(10)

式中，α为结构面倾向；β为结构面倾角；T为角度常量，取值为：

(11)

1.4 方法可靠性分析

为了进一步验证改进的点云聚类算法的可行性，对六面体点云数据进行聚类分析。六面体点云的示意图见图 2。在六面体点云所在的空间直角坐标系中，y坐标轴的正值方向代表北方向，x坐标轴的正值方向代表东方向，z坐标轴代表垂直方向。六面体点云的6组产状是已知的，如表 1所示。

完成六面体点云数据的法向量计算后，使用改进的模糊C均值算法对这些法向量进行聚类分析。完成聚类后，将原六面体的产状与各聚类中心的产状进行对比，结果如表 2所示。由表可知，聚类提取的产状与原六面体产状间的最大误差为0.08°，最小误差为0.06°，准确度较高。由于该点云数据具有高度规则性且表面无明显起伏，识别的精准度较高，所有误差均小于0.1°，满足要求。聚类识别效果良好，验证了改进的模糊C均值算法的有效性。六面体点云的聚类结果共有6种，其中，J1代表第1组和第4组聚类中心的目标点云；J2代表第2组和第6组；J3代表第3组和第5组。由图 3可知，六面体分类结果与其实际产状相一致，聚类效果较好。

1.5 结构面间距与迹长的自动化提取

在完成结构面的识别与产状计算后，可基于闭合结构面提取间距和迹长。结构面间距与迹长是评价岩体质量的关键指标。结构面迹长的提取主要依赖点云数据，通过计算点云中各点之间的距离来确定迹长。然而，由于结构面表面风化造成的凹凸差异、坡体表面的植被干扰及数据获取过程的不确定性等因素，许多结构面之外的点云数据难以有效识别和剔除，从而影响几何数据的提取，导致计算结果不准确。

为解决上述测量方法的不足，本研究采用Fast Marching算法对聚类后的点云数据进行封闭成面处理(该步骤已在1.3节的结构面识别中完成)，将闭合结构面内相距最远的2个边界转折点之间的最大距离作为迹长。通过Matlab软件编写算法，自动计算同一组优势结构面内所有结构面的迹长，实现自动计算，具体步骤如下。

(1) 结构面边界转折点的获取：首先，选取一组闭合结构面，将其导入到ArcGIS软件中；然后，使用要素折点转点工具，添加x，y坐标工具，提取该组结构面的所有边界转折点及其坐标信息。

(2) 边界转折点间距的计算：根据获得的单个结构面边界转折点坐标，利用式(12)直接计算任意两转折点之间的距离。

(12)

式中，i为结构面上某一转折点，i=1, 2, 3, …, n-1；j为结构面上除i以外的另一转折点，j=i+1, i+2, …, n。

(3) 单个结构面迹长的计算：基于式(12)，通过Matlab软件对同一结构面上任意两转折点间的距离进行循环迭代计算，直至完成所有点之间的间距计算，所得到的最大值即为单个结构面的迹线长度l。

(4) 结构面迹长的计算：利用Matlab编程重复步骤(2)和步骤(3)，即可提取一组结构面的迹长信息；然后，根据式(13)计算出该组所有结构面的平均迹长，即为该组优势结构面的迹长。

(13)

式中l_i为第i个结构迹长面的迹长。

对于结构面间距的计算，本研究基于Fast Marching算法提取的闭合结构面，并结合前人的方法^[25]，通过Matlab编写代码算法，自动计算同一组优势结构面内相邻结构面间的距离，即为2个结构面的间距。结构面间距计算的流程如下。

(1) 结构面的选取：从优势结构面组中选取一组闭合结构面，并对它们的平面参数a，b，c进行拟合，以获得最佳值；所选取结构面间的区别仅在于到原点的距离d。结构面方程为：

(14)

(2) 相邻结构面间距的计算：将各个结构面的距离d值按升序排列，通过式(15)计算相邻d值间的距离x_i，即为两相邻结构面间的间距。

(15)

(3) 结构面间距的计算：通过式(15)得到各相邻结构面间的距离后，根据式(16)计算出该优势组所有结构面间的平均间距，即为该组优势结构面的间距。

(16)

2 方法应用 2.1 研究区概况

研究区位于四川省G108国道北侧凉山彝族自治州会理县境内，坐标为102°8′34″E，26°34′18” N，海拔高程为2 245~2 455 m，高差为210 m。主要出露地层为震旦系观音崖组和灯影组并层的会理群，包括因民组、落雪组、黑山组。主要岩性为白云岩，下部砂岩、页岩、变质细碎屑岩、碳酸盐岩夹火山岩、火山碎屑岩和灰岩。根据现场走访，该斜坡在2014年10月发生过一次崩塌，其面积为2.4×10⁴ m²。崩塌岩体的主要岩性为灰岩。

2.2 无人机摄影测量

本研究采用无人机摄影测量技术获取岩质斜坡的多角度高清影像。随后，使用Metashape Professional软件对采集的影像进行三维模型重建，从而建立研究区的三维点云模型。所获得的三维点云模型具有清晰的斜坡结构面纹路与特征，有助于后续的分析与研究。

3 结果分析与验证 3.1 结构面产状的提取与算法对比

使用随机霍夫变换算法计算法向量，获得点云数据的法向量信息后，应用改进的模糊C均值算法对研究区的点云数据进行聚类分析。接着，利用Fast Marching算法获取结构面的产状信息，最终计算获得了6组优势结构面，如图 4所示。此外，还对结构面J1，J2，J4的典型剖面进行提取，如图 5所示。将提取结果与传统方法(即模糊度取2和3)获取的产状进行了对比，结果如表 3所示。

由表 3可见，改进后的算法在提取结构面产状方面的效果更佳。改进后的模糊C均值算法获得的结构面产状与模糊度取2和3的结果总体呈现出一致的数据特征，但前者的产状的标准差优于后者。相比之下，每组产状的标准差最高减小了2°和4°。这表明使用改进的模糊C均值算法所获得的结构面产状结果更为准确，验证了算法改进的有效性。

3.2 结构面间距与迹长的提取及精度验证

对识别出的结构面进行了间距和迹长的自动提取，并将提取结果与原始三维点云模型手动测量的值进行了对比，以验证自动提取算法的精确性。手动测量结构面间距的方法是根据结构面的产状与位置建立基准面，并测量相邻基准面之间的距离作为相邻结构面的间距。此外，使用Cloud Compare软件的距离测量工具，手动测量点云模型中结构面的迹长。由于结构面的数量较多，只将优势结构面J1和J4自动提取的结果与手动测量值进行比较，结果见表 4和表 5所示。

将结构面间距的手动测量与自动提取结果进行对比，发现其误差范围在厘米级，最大误差为1.1 cm，相对误差小于2%，误差较小，满足区域信息测量的误差容许要求和实际应用需求^{[25, 32]}。分析误差产生的主要原因是岩体出露结构面的表面不平，导致基准面位置的选取偏移。而间距的自动提取算法则是选取最佳基准面来计算间距，其结果比手动测量的结果更为合理。与手动测量结果相比，结构面迹长自动提取结果的最大误差为3.4 cm，相对误差小于2%，满足准确度要求。

通过对研究区岩体结构面的点云数据进行识别与信息自动提取，获取了岩体结构面的产状、间距、迹长及结构面的数量(见表 6)。

由图 4和表 6可知，研究区内主要发育6组结构面。结合剖面图(见图 5)与现场识别情况，分析该地区结构面类型结果如下：结构面J1主要表现为顺坡向的缓倾角结构面；结构面J2为斜向结构面；结构面J3和J4可能是因卸荷作用形成的卸荷结构面；结构面J5为顺坡向结构面；结构面J6可能是由于岩体风化作用形成的次生结构面。其中，结构面J3和J6发育较多，J1，J2，J5发育较少，结构面J2和J5组成X状结构面。由于该地区在2014年10月发生过一次崩塌，崩塌的上部呈现三角形，结构面J2和J5存在交叉现象，组成X状结构面。应用本研究的方法对研究区斜坡进行调查，结果表明，控制性结构面为J2和J5。需要注意的是，由于岩石表面的不均匀性，可能会导致算法识别的各组结构面之间存在部份交叉的现象，从而导致结构面显示不完全的情况。根据表 3的分析可知，该情况并不影响结构面的识别和产状的计算。

4 结论

本研究提出了一种基于点云数据的岩体结构面识别与信息自动提取的方法。通过引入新的参数改进模糊C均值算法，实现了模糊度的自适应调节，从而优化聚类结果并增强算法的鲁棒性。结合Fast Marching算法，实现了结构面间距与迹长的自动提取。将该方法应用于四川省G108国道北侧某边坡，结果表明，改进算法的产状标准差最高减少了4°，显著提高了识别精度。自动提取结构面间距与迹长的方法大幅提高了效率，最大测量误差仅为3.4 cm，相对误差不超过2%。这一结果不仅验证了所提方法的有效性，也为岩体结构面分析提供了一种高效、准确的技术方案。