岩体中广泛分布的不连续结构面使岩体具有强烈的不连续性、非均质性和各向异性.结构面对岩体力学性质的影响主要取决于结构面的发育状况.如岩性完全相同的两种岩体，由于结构面的空间方位、密度、规模、宽度、表面起伏形态及其组合关系不同，它们的力学性质会有很大的差异.因此，对岩体中的结构面各项参数进行实地测量和统计分析对于研究岩体的变形规律及破坏模式是至关重要的.

目前，对岩体结构面空间方向组合规律的研究主要是对结构面按照其产状进行归类，即岩体结构面产状优势分组.传统方法是利用玫瑰花图或极点密度等值线图进行目测分析，结合分析者的工程经验人为确定分组结果.显然，这种方法带有很强的主观性，缺少严格的理论依据，因此，有必要采用客观有效的数学方法对结构面产状进行优势分组.Shanley等^{[1, 2]}首次运用聚类算法分析结构面方向数据，该方法的结果受小球半径取值的影响很大，而如何确定合理的小球半径的问题一直没有得到有效解决.Hammah等^[3] 提出利用fuzzy C-means(FCM)的结构面组数划分方法，但此方法是一种局部搜索寻优法，且计算结果受初始聚类中心的影响较大.因此，蔡美峰等^[4]提出应用遗传算法来克服模糊C均值聚类算法的不足，但此方法需要编码和解码，操作较为复杂.卢波等^[5]采用小生境Pareto遗传算法，将结构面产状优势分组问题转化为多目标优化问题进行求解，由于遗传算法计算过程复杂，使得此方法也存在一定的缺陷.

随着岩体工程规模的逐渐扩大，需要分析的结构面数量急剧增加，寻找一种简单高效的结构面优势产状分析方法已经成为工程技术人员与相关学者亟须解决的问题.鉴于以上方法存在的不足，本文应用一种全新的智能算法——萤火虫算法来搜索最佳的聚类中心，进而获得合理的岩体结构面产状优势分组结果.萤火虫算法具有出色的全局搜索能力和较快的收敛速度，其结果要明显优于遗传算法和粒子群优化算法.最后将该算法应用于松塔电站坝址区平硐实测结构面产状数据的分组分析.

在现场调查统计结构面数据时，结构面的产状通常用倾向(α)和倾角(β)表示，Θ = (α,β)，其中，0°≤α＜360°，0°≤β≤90°.在结构面优势产状分析过程中，通常不考虑其厚度及起伏程度而将其简化为一空间平面，并用对应的单位法向量 X =(x,y,z) 来表示结构面的空间方位：

岩体结构面产状数据划分问题就是将空间方向相近的结构面分别归类分组，找出其优势方向，因此，首先需要选择合适的数学指标来度量结构面产状的相似性.在岩体结构面产状优势组数划分问题中存在一种特殊情况，即倾向相差约180°的两个高陡倾角结构面是相似的，应该划分为同一组.因此本文采用两结构面所夹锐角的正弦值来判断两结构面方向的接近程度.单位法向量分别为 X₁=(x₁,y₁,z₁)和 X ₂=(x₂,y₂,z₂)的两结构面所夹锐角为

给定N个结构面P_j(j=1,2,…,N)，划分组数C，每组的聚类中心为V_i(i=1,2,…,C ).定义u_ij为第j个结构面属于第i个聚类中心的隶属度^[6]：

萤火虫算法(firefly algorithm,FA)是由剑桥大学学者Yang^[7]在2008年提出的一种模拟萤火虫群体行为的随机优化算法.其计算原理可理解为：用萤火虫个体表示问题的解，根据其所处位置的亮度来度量解的优劣，利用亮度大的萤火虫吸引亮度小的萤火虫向其移动的过程来模拟问题的解向最优解逼近的过程，通过种群的不断更新最终获得问题的最优解^[8].

在萤火虫算法中，亮度与吸引度是两个非常重要的概念.特定位置s的萤火虫的亮度I(s)是由解空间对应位置的目标函数值决定的.萤火虫的吸引度η是随着个体之间的距离r的变化而变化的，同时也与萤火虫所处环境的介质密切相关.萤火虫的吸引度η由下式确定：

萤火虫i受到比其明亮的萤火虫k的吸引而向k的位置移动，萤火虫i移动后的位置由下式决定：

萤火虫算法具有操作简单、鲁棒性强、易于实现等特点.与遗传算法和粒子群算法相比，其搜索全局最优解的能力更强，收敛速度更快^[9]，目前已被应用到非线性问题优化^[10]、图像压缩处理^[11]等方面.

将萤火虫算法应用到岩体结构面优势方向的搜索，式(5)作为目标函数，约束条件为倾向α∈[0°,360°]，倾角β∈[0°,90°].通常意义上的萤火虫算法是要寻找亮度最大的个体，而结构面优势分组问题是要寻找当分组数给定时使目标函数最小时的聚类中心，为此，将s处萤火虫的亮度I(s)定义为

应用萤火虫算法进行岩体结构面产状数据划分的计算流程如图 1所示.

图 1(Figure 1)

Fig. 1

图1 算法流程图 Fig. 1 Flow chart of the algorithm 图中t，i和k均为循环变量.其中，t表示种群当前更新次数，i和k表示萤火虫个体

在寻找到对应于不同分组数C的一系列最优解后，通过聚类有效性函数来评价这些聚类结果的优劣，进而确定最佳聚类数.划分系数v_PC和Xie-Beni指标v_XB是常用的衡量聚类有效性的函数：v_PC越大，聚类结果越好；v_XB越小，聚类效果越好.计算公式如下：

怒江松塔水电站位于西藏自治区察隅县察瓦龙乡境内，坝址区属高山峡谷地貌，岩性主要为燕山晚期黑云二长花岗岩，节理裂隙比较发育.为了进行围岩稳定性评价，首先要对结构面的发育规律进行分析研究，确定结构面优势方向.本文选取坝址区右岸某高程相邻的四条平硐，实测1439条结构面产状数据，图 2为结构面的极点等密度图.应用文中提出的算法对这些结构面进行优势分组.算法参数设置为：萤火虫群的规模n=20，种群最大更新次数T_max=50，η₀=1，γ=1，ξ=0.1.计算不同的分组数对应的聚类有效性检验结果，见表 1.

图 2(Figure 2)

Fig. 2

图2 结构面极点等密度图 Fig. 2 Iso-density plot of discontinuity poles

表1

Table 1

表1(Table 1)

表1 聚类结果有效性检验 Table 1 Validity test of clustering results 指标 分组数C 2 3 4 5 6 7 8 9 10 vPC 0.62 0.64 0.56 0.52 0.48 0.45 0.43 0.41 0.37 vXB 0.21 0.11 0.32 0.29 0.34 0.31 0.33 0.32 0.29

表1 聚类结果有效性检验 Table 1 Validity test of clustering results

从表 1中可以看出，划分系数v_PC和Xie-Beni指标v_XB同时显示，当分组数为3时聚类效果最好，所以综合确定最优分组数为3，此时各组的中心产状及各组的结构面数目见表 2，详细分组结果见图 3a.

表2

Table 2

表2(Table 2)

表2 分组结果 Table 2 Clustering results of discontinuity data 组号 1 2 3 倾向/(°) 246.61 155.43 109.36 倾角/(°) 84.27 88.97 34.40 结构面数量 300 298 841

表2 分组结果 Table 2 Clustering results of discontinuity data

图 3(Figure 3)

Fig. 3

图3 分组极点图 Fig. 3 Pole map of clustering result (a)—新方法； (b)—文献[12]方法.

对比图 2与图 3a不难看出，新方法的分组结果是符合实际情况的，此方法不但准确搜索出各组的聚类中心，还给出了清晰的组间边界，这是传统方法做不到的.利用Xu等^[12]提出的变尺度混沌优化模糊聚类方法计算得出的分组结果如图 3b所示，可以看出，两种方法的分组结果是一致的，这也进一步说明了新方法的有效性.但由于后者采用混沌搜索具有一定的盲目性，计算时间很长(约90min)，而新方法的计算时间不足1min，说明新方法具有较高的计算效率，优势明显.图 4表示种群的更新过程，可以看出，算法的收敛速度是非常快的.

图 4(Figure 4)

Fig. 4

图4 群更新过程 Fig. 4 Population updating process

1) 基于萤火虫算法的岩体结构面产状数据优势组数划分的新方法能够自动识别各组的聚类中心以及清晰合理的组间边界，有效避免了传统的图形分析法主观性强的缺点.

2) 与文献^[12]方法的分组结果进行对比，进一步证明了新方法的正确性.但是，值得注意的是，新方法耗时远远低于文献^[12]方法，这一优点在工程应用中是非常重要的.因此，新方法具有较大的实际意义.