基于激光测量和FKM聚类算法的隧洞岩体结构面的模糊群聚分析

刘昌军，丁留谦，张顺福，廖井霞

中国水利水电科学研究院，北京 100038

0 引言

在隧道（洞）施工过程中，岩体结构面的存在对隧洞安全有着重要影响。快速、准确地获取岩体结构面的数量、几何信息和优势产状等参数是隧洞工程岩体稳定性评价和地质灾害防治等工作的前提和基础[1-2]。目前隧洞岩体结构面信息获取主要以人工现场接触测量为主，劳动强度大、效率低下，且难以满足快速施工要求。随着激光测量和摄影测量的快速发展，非接触测量已经应用在各行各业中。近几年来在岩土行业中得到快速发展，如高陡边坡地质勘察和测量、危岩体调查、岩体结构面统计、边坡变形监测等[3-6]。

近年来国内外学者也开始尝试应用摄影测量技术[3]和激光测量技术等非接触测量方法进行高陡边坡的岩体结构面的调查，取得了一定成果[4-10]。然而利用非接触测量方法研究隧洞岩体结构面的研究还较少，一般情况下，隧洞内可见光不足，采用摄影测量方法进行调查往往受到一定限制，而激光测量方法不受可见光影响，在隧洞调查中具有较强的应用优势。

针对隧洞岩体结构面人工现场量测和统计分析的难点，结合广泛使用的激光测量技术和激光点云数据处理的难点，提出了一整套应用激光测量技术进行隧洞岩体结构面的全自动统计分析方法，特别是针对隧洞表面激光点云数据的特点，应用柱面投影的Delaunay生长法对激光点云数据进行了三角网重构。基于激光测量技术和FKM（Fuzzy KMean）群聚方法进行隧洞岩体结构面优势产状统计分析的流程与方法如下：1）首先利用激光扫描仪获取隧洞表面岩体的激光点云数据，结合全站仪测定的反射片的坐标数据将各站激光点云数据进行坐标转换和拼接，使得拼接后的点云数据的y坐标与北方向重合。2）对处理后的激光点云数据进行柱面坐标投影，采用Delaunay生长算法对投影后的激光点云数据构建三角网，将三维点坐标数据按照平面建立的点云拓扑关系进行点号连接，得到隧洞岩体表面的空间TIN网格模型。3）利用TIN网格模型，求得的每个三角网平面方程的法向向量，以三角网平面方程的法向向量计算其产状。4）采用FKM模糊群聚方法和改进的方位矩阵方法对每个小三角网的产状进行全自动群聚分类，用三角网的产状聚类结果代表不同组的优势产状，并用不同颜色表示结构面的群聚分类结果。

基于以上原理，采用自主研发的海量激光点云数据的处理和岩体结构面统计分析软件FSS[4]对杭州某隧洞的岩体结构面进行了统计分析，得到了各组结构面的优势产状及其空间分布。工程实例的研究结果表明，该方法在隧洞内岩体结构面的调查中具有广阔的应用前景，为隧洞岩体结构面优势产状的快速调查和自动统计分析开辟了新的途径。

1 基于柱面投影的点云重构方法的实现

三维激光点云构网方法主要采用基于面剖分、体剖分和面投影的3种方法建立[11-14]。基于面投影是先将三维空间数据投影到平面上进行平面构网，再映射到三维空间最终达到三维构网目的。而基于面投影法可分为基于平面投影、基于临近区域切平面投影、圆柱面投影、球面投影和椭球面投影等方法。

结合隧洞激光点云数据的特点，采用基于柱面投影的Delaunay构网方法对隧洞表面激光点云进行了构网。基本思路如下：将转换成统一空间坐标后的激光点云数据以柱面为媒介投影到平面上，得到平面点坐标，在平面内构造平面三角网，最后将平面三角网映射到空间完成三维构网。

1.1 柱面坐标系的建立及投影

1.1.1 柱面坐标系

利用柱面投影坐标系将三维激光点云投影在二维平面内构三角网是一种快捷有效的构网方法。郑德华等[15]和徐秀川等[11]采用柱面投影法研究了三维激光点云的三角网构建方法，其核心是在保证封闭物体表面点之间的连接关系的前提下建立柱面坐标系统。

如图1a所示建立柱面坐标系，三维直角坐标系为左手坐标系，圆柱面底面圆心O为坐标中心，圆柱底面半径为R，圆柱面切割母线AA′。在上述坐标系中，设P（X，Y，Z）为任一点激光点的坐标数据，在圆柱面上投影坐标为P′（r，ω，h）。圆柱面长度为h。

图1 柱面坐标系和投影后的平面坐标系Fig.1 Cylindrical coordinate system and plane coordinate system after the projection

1.1.2 柱面坐标系统的转换

三维激光扫描获取的点云坐标系一般是以测站为原点的左手坐标系，为了建立三角形格网首先应将其转换到柱面坐标系统中。具体方法如下：三维激光扫描点云坐标中心点与柱面坐标系统的底面圆心点重合，将三维点坐标沿柱面半径向圆柱基准面投影[16]。三维空间直角坐标系和柱面坐标系的转换公式为

1.1.3 投影后的平面坐标转换为空间坐标

将柱面坐标系转换为平面二维坐标系，具体方法为文献[16]提出的沿圆柱基准面并沿切割母线AA′切开并展开的方法。展开后的平面坐标原点为图1b中的A′点，其中左侧的AA′轴为x轴，向上为正向，A′A′为y轴，向右为正向，见图1b。其转换坐标系矩阵为

式中，θ为切割母线在圆柱基准面上顺时针到XZ平面的角度。

1.2 基于柱面投影的构网方法的实现

激光点云数据柱面投影构网的基本思路是，首先确定圆柱面底面中心及半径以确定进行投影的圆柱面，然后将点云数据投影到圆柱面上，将直角坐标变换为柱面坐标。再将柱面坐标展开到平面坐标，在平面坐标系中进行构网。在IDL具体实现过程中，需构网的点云数据的数组为pointData（一个3×N的2维数组）。完成柱面投影构网算法包含以下几个步骤：

1）选择柱面底圆中心orign：根据图建立的柱面坐标系，确定激光点云数据坐标中心平移到柱面坐标系的底圆圆心坐标。

2）确定圆柱底圆半径radius：柱面底圆半径的赋值原则是让所确定的球与点云中每个点都尽可能的接近，这样能保证点云投影到柱面上尽可能地减少变形。这需要制定拟合优化策略，但在本算法中这并不是重点，因而在此作简单处理。在pointData中随机抽取sampleNum个点，将其投影到柱面底面上，将各点坐标平均即可得到柱面底圆半径，经实验这种方法效果也较好。

3）点云投影到柱面上：柱面确定后将点云数据以底圆圆心为原点按照公式（1）进行坐标变换，利用IDL中的函数CV＿Coord将激光点云从直角坐标变换为柱面坐标pointData1。

4）柱面坐标系转换为平面坐标：将柱面坐标按照图1b方法展开为平面二维坐标，在平面坐标内利用IDL函数QHull构建Delaunay三角网，得到点云数据的拓扑关系，并存放至triangle数组中。

5）平面坐标转换为空间坐标：利用公式（2）将平面坐标转换为空间坐标，将三维空间点按照平面三角构网确定的拓扑关系根据点号进行连接，形成空间Delaunay三角网。

2 FKM方法模糊聚类的基本理论

2.1 FKM模糊聚类方法

本研究将采用FKM计算方法进行结构面产状分类的分析[17]。FKM算法通过引入隶属度函数表示属于不同类别的程度对数据进行划分：首先估计聚类的中心，其次反复调整聚类中心，使数据集中的每个点距每个中心的距离之和达到最小或满足终止条件。

给定n个数据样本集合X＝（x1，x2，…，xn），观测样本xi的特性指标为xi＝（xi1，xi2，…，xin），也即每个样本具有m个属性，n个数据对象分成k（2≤k≤n）个不同类。数据集合X中的对象xi以一定的隶属度隶属某一类，每一类对应一个模糊矩阵U。如果X中含有k个类别，则有这k个聚类中心构成的向量矩阵为V，由此可建立FKM聚类算法的目标函数：

式中：u为法向向量对各群心法向量的隶属度矩阵；Xj≡（xj，yj，zj）为平面的单位法向向量，xj，yj，zj为Xj的分量，计算见式（4）；Vi≡（xi，yi，zi）为各群聚的群心法向向量；m为模糊度，由于m控制隶属度在各类之间共享的程度，所以m越大，模糊性就越大，d（Xj，Vi）为Xj与Vi的距离。

式中：β为岩体结构面的倾角；α为结构面的倾向。

R.E.Hammah和J.H.Curran[8]提出将结构面产状转换成在球体（半径为1）上的单位法向向量，然后计算单位法向向量与群聚的群心法向向量的正弦值，作为此两向量间的距离值：

FKM聚类算法是在满足隶属度uij值为

的条件下，求解目标函数J的最小值，即

且元素uij满足以下条件：

利用 R.E.Hammah和J.H.Curran[8]提出改进的方位矩阵S＊计算新的群心法向量[2，10]：

利用改进的方位矩阵可以改进新群心法向向量在椭圆形中不均匀分布情况。

取m＝2，利用反复迭代运算，寻找最终群心法向向量及相应隶属度矩阵。迭代收敛以两次隶属度矩阵的各元素差异值的绝对值小于允许值ε为准，即

式中：u为第一次隶属度矩阵；u′为第二次隶属度矩阵。

2.2 激光扫描和FKM聚类耦合方法的技术流程

基于激光扫描和FKM聚类耦合方法的岩体结构面优势产状统计的技术流程如图2所示。

2.3 激光扫描和FKM聚类耦合方法的适用条件

基于激光扫描技术和FKM聚类耦合方法研究岩体结构面优势产状比传统方法具有一定的优势，主要体现在以下几个方面：

1）采用激光扫描技术可以进行远距离岩体结构面快速扫描和统计，特别适合现场无法到达和存在危险区域的岩体结构面统计分析。

2）采用改进的方位矩阵和FKM模糊聚类方法可以模拟岩体结构面的发育分布规律，弥补了等密度图、玫瑰图等传统统计方法的不足，量化的给出了结构面自动分组及各组的优势产状，且可以分组进行三维空间显示。

3）由于三维激光扫描是自发光源，在一些无可见光或可见光不足情况下（如地下洞室、隧洞和巷道等）仍然可以使用，可以实时获取和分析施工过程岩体裂隙发育情况，用于指导施工和设计。

4）该方法是利用岩体表面三角网的群集特性进行的岩体结构面统计，岩体表面所有的结构面都参与了群集分类，三角网平面产状的群聚结果可以作为岩体结构面各组的优势产状。而传统的结构面优势产状统计都是只有局部或部分结构面产状参与了统计分析。

图2 基于激光点云的岩体结构面FKM聚类算法流程图Fig.2 Flow diagram of FKM clustering algorithm based on laser scanner

虽然该方法在岩体结构面自动统计中具有较大优势，但也存在一定的局限性：

1）该算法特别适合于裸露岩体，在一些植被茂密的岩体上，该方法有一定的局限性。

2）由于该方法是用三角网的产状替代裂隙产状，因此在一些岩体表面较为光滑，且不能区分结构面时，就会大大降低其统计结果的精度。

3）对于一些特别细小的结构面，在远距离扫描时，由于扫描精度问题，细小的结构面不能被识别出来，将影响统计结果的精度。

4）如果岩体表面特别潮湿或有地下水渗出的情况，该方法的使用会受到一定的限制。

2.4 激光测量和FKM聚类耦合结果的分析方法

使用基于激光测量和FKM聚类的耦合分析方法研究岩土结构面的优势产状与传统的统计方法不同。该方法得到的结果是岩体表面三角网的优势产状的聚类结果，用三角网产状的聚类结果来代替岩体结构面分组优势产状，而且随着不同群聚数的增加，三角网产状的聚类组数也相应增加。将该方法应用到实际工程时，需要利用人工现场判别结构面的分组将不同群聚数的分组结果加以区分，对人工开凿和非结构面的聚类结果进行剔除，进而确定合适的群聚数和分组结果。

图3 FSS软件系统框架和工作流程Fig.3 The system framework and workflow of FSS software

3 软件编制

采用IDL语言编制了基于激光点云数据的岩体结构面处理软件FSS，该软件主要包括数据处理模块、结构面统计分析模块和制图输出模块3个主要应用模块。其中：数据处理模块主要包括激光点云的处理功能（点云的局部删除、抽稀和坐标转换等）、点云的构网（基于平面、球面和柱面等投影算法的生长法构网）等功能；结构面统计分析模块包括单个结构面的空间拟合、绘制迹线和岩体结构面的模糊群聚分析等功能[4]。FSS软件的系统框架和工作流程见图3，软件界面见图4。

4 应用案例

4.1 研究区基本情况

以杭州将台山某人工开挖隧洞的岩体结构面统计分析为例。杭州将台山隧洞位于将台山东南麓，毗邻八卦田，遥望钱塘江。隧洞位于将台山东南麓边坡的西侧，为丁字形交叉洞室，洞口呈半圆形。该洞总长约200m，裂隙较为发育，多闭合，无充填，延伸长度较短。

4.2 现场扫描及数据处理

首先利用地面激光扫描仪对岩洞内的表面进行测量，获取岩洞内岩体的激光点云坐标约70万个，同时利用全站仪获取现场布置的15个圆形反射片中心点大地坐标。利用15个反射片中心点大地坐标将扫描获取的激光点云数据进行坐标转换，将转换成大地坐标的激光点云数据导入FSS软件，有关FSS软件功能及界面见文献[4]。本次调查选择了东西向的裂隙分布较为密集的一段山体隧洞，从扫描数据可查询该岩洞长约15m，高约5m，直径约6.5m。图5为洞口附近现场照片和洞口附近激光点云数据（点云赋值照片上的RGB色彩）。图6为隧洞研究段的点云图。

根据获取的激光点云数据，为了构网方便，将扫描获取的激光点云数据抽稀一半，对抽稀后的激光点云采用第2节介绍的方法进行空间三角网重构，构网后的隧洞三维网格图见图7。

图4 FSS软件界面图Fig.4 The interface diagram of FSS software

图5 洞口局部的现场照片和激光点云数据Fig.5 Local scene photos and laser point cloud data of the entrance of cave

图6 隧洞激光点云图Fig.6 Diagram of laser point cloud of the tunnel

图7 隧洞三维TIN网图Fig.7 3-D mesh figure of the tunnel

4.3 结构面模糊统计分析

从现场观察可知，该隧洞主要有5组岩体结构面较为发育。利用隧洞表面的空间TIN网，采用FKM的模糊群聚方法对隧洞存在的结构面优势产状进行统计分析。图8给出了不同群聚数下（k＝4～9）的隧洞岩体结构面的空间分布和优势产状，表1给出了相应群聚数下（k＝4～9）的结构面优势产状统计值和现场实测值。

由表1可知：第6组结构面在k＝6以后和第3组结构面倾向和倾角较为接近，第7组结构面在k＝7以后和第1组结构面的倾向和倾角较为接近；2组结构面在k＝6以后其优势产状的变化量变小，且在k＝6以后才能区分2组结构面。而第9组结构面在k＝9以后才出现，且均不和前8组结构面相近。因此，对于该隧洞的岩体结构面取k＝6的分类结果较为合适。从分类的结果图（图8）上可看出，第8组结构面为隧洞的地面和顶面，且不同群聚数下的该组结构面的分类结果呈现离散性，不作为单独一组。因此从k＝6～9的分类结果看，隧洞内的岩体结构面的优势产状共5组，且和实际测量值接近。另外，由以上模糊群聚的自动分析结果可知，计算的群心法向向量（即优势产状）随群聚数的增加，均出现变小的结果，即优势产状随着群聚数增加，其值变化不大，且能将产状接近的不同结构面加以辨识。

由上面的应用案例的统计结果可看出，本文提出的基于激光点云数据和FKM的岩体结构面群聚分析方法可用于复杂洞室的结构面优势产状统计和分类，计算结果可靠，计算速度快，能够适应于隧洞内特殊条件的岩体稳定性分析。另外，由激光点云数据可以对隧洞裂隙的几何信息（如迹长、间距和缝宽等）进行分析，也可给出其地质编录图，详细研究成果见文献[4]。

图8 不同群聚数下的岩体结构面模糊群聚结果Fig.8 Clustering results of rock mass discontinuity occurrence under different clustering numbers

表1 不同群聚数下的岩体结构面优势产状统计Table1 Statistics of different number of the attitude of rock mass discontinuity in advantage

5 结论

1）将激光扫描技术和FKM的模糊群聚方法应用于隧洞等复杂条件下的岩体结构面统计分析中，该方法可以实现非接触测量的岩体结构面优势产状的自动聚类分析。群聚结果表明，部分群聚的几何中心（群心值）随着群聚数目增加，其变化量将会减小，该群心值类别的颜色可反映其结构面的空间分布，且群心值可作为该结构面优势产状的代表值。工程实例应用表明，该方法计算结果可靠、合理，具有较强的通用性和推广性。

2）采用的基于柱面投影的Delaunay生长法进行复杂激光点云的构网方法速度快、构网质量高、重构得到的三角网表面与被采样的岩体表面拓扑差别最小。同时还提供了基于球面投影和平面投影的激光点云的生长法的构网算法，解决了复杂条件下的激光点云的快速构网问题。

3）利用自主研发的基于海量激光点云数据的岩体结构面统计分析软件FSS，可以快速对扫描获取的岩体表面的激光点云数据进行岩体结构面统计分析和编录工作，减少了高陡边坡人为攀爬的危险，减少了现场工作和数据后处理工作。

（References）：

[1]王述红，张靖杰，李云龙，等.工程岩体开挖过程全空间块体搜索及系统研制[J].东北大学学报：自然科学版，2010，31（7）：1026-1029.Wang Shuhong，Zhang Jingjie，Li Yunlong，et al.Quickly Spatial Identification Process of Blocks During Rock Excavation and Its System Modeling[J].Journal of Northeastern University：Natural Science，2010，31（7）：1026-1029.

[2]王述红，杨勇，王洋，等.基于数字摄影测量的开挖空间模型及不稳定块体的快速识别[J].岩石力学与工程学报，2010，29（增刊1）：3432-3437.Wang Shuhong，Yang Yong，Wang Yang，et al.Spatial Modeling and Quickly Identification of Unsteady Rock Blocks Based on Digital Photogrammetry [J].Chinese Journal of Rock Mechanics and Engineering，2010，29（Sup.1）：3432-3437.

[3]王凤艳，陈剑平，杨国东，等.基于数字近景摄影测量的岩体结构面几何信息解算模型[J].吉林大学学报：地球科学版，2012，42（6）：1839-1846.Wang Fengyan，Chen Jianping，Yang Guodong，et al.Solution Models of Geometrical Information of Rock Mass Discontinuity Based on Digital Close Range Photogrammetry[J].Journal of Jilin University：Earth Science Edition，2012，42（6）：1839-1846.

[4]刘昌军，丁留谦，孙东亚.基于激光点云数据的岩体结构面全自动模糊群聚分析及几何信息获取[J].岩石力学与工程学报，2011，30（2）：358-365.Liu Changjun，Ding Liuqian，Sun Dongya.Automatic Fuzzy Clustering Analysis and Geometric Information Acquisition of Rock Mass Discontinuity Based on Laser Point Cloud Data[J].Chinese Journal of Rock Mechanics，2011，30（2）：358-365.

[5]刘昌军，高立东，丁留谦，等.应用激光扫描技术进行岩体结构面的半自动统计研究[J].水文地质工程地质，2011（2）：52-57.Liu Changjun，Gao Lidong，Ding Liuqian，et al.Research of Semi-Automatic Statistics of Rock Mass Discontinuity Applying Laser Scanning Technology[J].Hydrogeology and Engineering Geology，2011（2）：52-57.

[6]刘昌军，丁留谦，孙东亚.三星堆月亮湾城墙遗址覆土方量计算[J].水利水电科技进展，2011（2）：81-84.Liu Changjun，Ding Liuqian，Sun Dongya.Calculation of Soil Quantity in Moon Bay Wall Site of Sanxingdui[J].Advances in Science and Technology of Water Resources，2011（2）：81-84.

[7]Feng Q，Sjögren P，Stephansson O，et al.Measuring Fracture Orientation at Exposed Rock Faces by Using a Non-Reflector Total Station [J].Engineering Geology，2001，59（1／2）：133-146.

[8]Hammah R E，Curran J H.Fuzzy Cluster Algorithm for the Automatic Identification of Joints Sets[J].International Journal of Rock Mechanics and Mining Sciences，1998，35（7）：889-905.

[9]王凤艳，陈剑平，付学惠，等.基于VirtuoZo的岩体结构面几何信息获取研究[J].岩石力学与工程学报，2008，27（1）：169-175.Wang Fengyan，Chen Jianping，Fu Xuehui，et al.Study on Geometrical Information of Obtaining Rock Mass Discontinuites Based on VirtuoZo[J].Chinese Journal of Rock Mechanics and Engineering，2008，27（1）：169-175.

[10]洪子恩，冯正一，吴宗江.应用三维激光扫描于岩坡露头位能的量测[J].水土保持学报，2007，39（3）：247-267.Hong Zien， Feng Zhengyi， Wu Zongjiang.Application of 3DLaser Scanning on Measuring Attitudes of Rock Slop Outcrops[J].Journal of Soil and Water Conversion，2007，39（3）：247-267.

[11]徐秀川，花向红，龚子桢，等.基于球面投影的散乱点云三维建模算法实现与效果分析[J].测绘工程，2011，20（3）：5-8.Xu Xiuchuan，Hua Xianghong，Gong Zizhen，et al.The Implementation Algorithm of Three-Dimensional Modeling to Point Cloud Based on Spherical Surface Projection and Effect Analysis[J].Engineering of Survey and Mapping，2011，20（3）：5-8.

[12]武晓波，王世新，肖春生.Delaunay三角网的生成算法研究[J].测绘学报，1999，28（1）：28-35.Wu Xiaobo，Wang Shixin，Xiao Chunsheng.A New Study of Delaunay Triangulation Creation[J].Acta Geodaetica et Cartographic Sinica，1999，28（1）：28-35.

[13]张帆，黄先锋，李德仁.基于球面投影的单站地面激光扫描点云构网方法[J].测绘学报，2009，38（1）：48-54.Zhang Fan，Huang Xianfeng，Li Deren.Spherical Projection Based Triangulation for One Station Terrestrial Laser Scanning Point Cloud[J].Acta Geodaetica et Cartographic Sinica，2009，38（1）：48-54.

[14]李逢春，龚俊，王青.基于三维TIN的精细表面建模方法[J].计算机应用研究，2006，23（8）：159-161.Li Fengchun，Gong Jun，Wang Qing.Highly Refined Surface Reconstruction Method Based on 3DTIN[J].Application Research of Computers，2006，23（8）：159-161.

[15]郑德华，张云涛.基于物体表面散乱三维激光扫描点的三角形格网建立[J].测绘工程，2004，13（4）：62-65.Zheng Dehua，Zhang Yuntao.Triangular Establishment Based on Disordered 3DLaser Scanning Point of Object Surface[J].Engineering of Surveying and Mapping，2004，13（4）：62-65.

[16]潘华伟，邹北骥.一种圆柱形全景图生成新算法及其实现[J].计算工程与科学，2003，25（6）：13-16.Pan Huawei，Zou Beiji.A New Algorithm for Generating Cylinder Panorama and Its Implementation[J].Computer Engineering and Science，2003，25（6）：13-16.

[17]田赤英，张旭男，宋士吉，等.改进的模糊C均值聚类算法及其在海底热液硫化物组分分析中的应用[J].海洋学研究，2010，28（4）：22-28.Tian Chiying，Zhang Xunan，Song Shiji，et al.Algorithm of Improved FCM Clustering and Its Application to Seafloor Hydrothermal Sulfide Composition Analysis[J].Journal of Marine Science，2010，28（4）：22-28.