基于花授粉算法的结构面自适应精细分组研究

2022-04-02王资平丁秀美黄星凯石广源

成都理工大学学报(自然科学版) 2022年2期

王资平, 林锋, 丁秀美, 黄星凯, 石广源, 王卫

(地质灾害防治与地质环境保护国家重点实验室(成都理工大学)，成都 610059)

岩体中结构面的发育特征是影响岩体力学性状的主要因素，精细描述岩体中结构面发育特征一直是工程地质和岩体力学领域的热点和难点问题[1-2]。岩体中结构面发育特征的精细描述，首先是进行结构面分组，完整的分组包括确定分组数、各组代表性产状和将测量的各条结构面归组。结构面是在岩体建造和改造过程中逐步形成的[1-3]，这是结构面分组的地质基础，也表明结构面在岩体中的分布并不是完全随机的。结构面分组通常依据结构面产状[1-4]，但也有一些学者建议分组时同时考虑间距、平面度、粗糙度、风化程度、长度和宽度等特征[5-8]。

对系统测量(如精测网法、测线法等)的大量结构面，常依据结构面产状采用极点等密度图、直方图和玫瑰花图进行分组[1-4]，这些方法简单、直观，但无法直接给出完整的分组结果，不便于进一步的统计分析。大量结构面的分组问题在数学上常概化为聚类分析问题。R.J.Shanley等[9]采用改进的一级模分析法依据产状进行聚类分析；R.E.Hammah等[5]采用模糊K均值聚类算法，分别依据“产状”和“产状+粗糙度”进行自动分组分析，但在寻优搜索时，较易陷入局部极小点[10]；王鹏等[11]引入遗传算法(GA)，与FCM算法一起依据产状分组，提高了算法的全局搜索能力；田景元等[12]采用综合模糊等价聚类方法分析优势产状；周玉新等[10]结合模糊等价聚类方法和模糊软划分聚类方法，形成综合模糊聚类方法，提高了分析效果；R.R.Jimenez等[13]采用谱聚类算法依据产状进行分组，该方法需指定分组数；卢波等[14]将分组问题概化为多目标优化问题，引入小生境技术和Pareto支配集理论进行求解，算法复杂；范雷等[15]采用改进的动态聚类方法和邓继辉等[16]采用改进的自组织聚类算法，均需要事先指定分组数；宋金龙等[17]将模糊聚类法与K均值动态聚类算法结合，改善了计算效率；徐黎明等[8]提出了基于变尺度混沌优化算法的多参数结构面优势组划分方法。综合分析这些研究成果可见，聚类分析可以用于结构面分组量化分析，但实现完全的精细分组还需要改进。

在水电开发等重大工程中，岩体特性现场调查会系统采集大量结构面数据。实践表明，中陡倾结构面在极点等密图上一般可以清晰分组，而缓倾角结构面难以分组[18]，为精细分析岩体中结构面发育特征，还需要改进聚类分析法。由于结构面产状与其成因密切相关，而成因决定了结构面其他特征，因此，本文将结构面分组问题抽象为以产状为基本变量的多目标组合优化求解问题，引入花授粉算法和Silhouette指标，建立一种新的聚类算法。花授粉算法(flower pollination algorithm, FPA)是英国剑桥大学Yang X.S.提出的一种模拟花朵授粉过程的种群智能优化算法[19]，随着该方法的进一步改进完善[20-23]和在函数优化、无线传感网、经济调度、图像分析等领域的推广应用[23-24]，可以解决多目标优化问题，具有全局搜索能力显著、收敛快等特点。

1 结构面分组问题的数学描述

1.1 结构面方位相似度及分组目标函数

结构面方位相似度d(Pi,Pj)定义为两结构面Pi和Pj的法线之间锐夹角φ正弦值的平方[8,14]，有

φ= arccos(Pi·Pj)

(1)

d(Pi,Pj) = sin2φ

(2)

假定待分组的N个结构面可划分为M组，其中任一组Gk(k= 1, 2, …,M)中有nk条结构面，Gk的中心记为Ck，称为分组中心。根据d(Pi,Ck)大小来判断Pi是否归属Gk组。用变量mik来描述归属关系

mik=1, (Pi∈Gk)

(3)

mik=0, (Pi∉Gk)

(4)

定义结构面分组目标函数(E)

(5)

由式(5)可知，当每一条结构面都单独成一组时，目标函数值为0；但实际分组数M是有限大小，最优分组时E值取最小值，这是一个组合优化问题。

1.2 最佳分组指标

Gk组中任一结构面Pj的轮廓宽度s(j,k)定义为

(6)

式中：a(j,k)为Pj与Gk组中所有其他结构面间方位相似度的均值；b(j) = min[a(j,i)]，(i≠k;i=1, 2, …,M)。

全局Silhouette指标的计算公式[25-26]为

(7)

Silhouette指标结合内聚度和分离度2种因素，既能评价聚类结果的优良程度，也能确定聚类个数。由式(7)可见，Gs取值范围为[-1, 1]。当Gs接近于1时，聚类效果较好；当Gs值接近于-1，则表示聚类结果错误。取Gs最大时对应的组数，作为最优聚类组数。

1.3 基于花授粉算法的最优分组求解

花朵授粉的过程分为2种，即异花授粉和自花授粉，前者对应目标函数最优求解的全局搜索，后者对应局部搜索[19-22]。花授粉算法在全局搜索与局部搜索之间的转换通过转换概率p控制[19]，p∈[0, 1]。

a.当算法执行全局搜索(异花授粉)时，由式(8)实现花粉位置更新。

(8)

(9)

其中，Г(λ)是标准伽马函数，并且该分布符合较大步长s>0的情况，常取λ=1.5。

b.当算法执行局部寻优(自花授粉)时，采用式(10)更新花粉位置。

(10)

根据上述原理，作者采用MATLAB平台编制了计算程序，主要计算步骤如下：

第一步：设定种群数Fn和最大迭代次数；初始化分组数M=2，转换概率p在文献[27]中建议取0.8，但作者通过多次试算发现p=0.4时收敛更快；删除输入数据中的相同数据，组成样本数据集P。

第二步：从P中任选M个样本作为初始聚类中心，据式(2)计算Pi与各聚类中心的方位相似度d，并对结构面归组，按式(5)计算目标函数值E。

第三步：寻优搜索，如果服从[0,1]均匀分布随机数rand>p，进行异花授粉，采用式(8)和(9)来计算更新解；如果rand

第四步：回到第二步，改变初始聚类中心。

第五步：记录最优解，计算Silhouette 指标；M=M+1，回到第二步，直至M>8。

第六步：比较Silhouette指标，找出Silhouette指标中最大的分组数，并输出最优解和所需图件。

2 结构面分组应用分析

西藏自治区芒康县境内澜沧江干流上某大型水电站坝址区基岩主要为英安岩和花岗岩，岩体中节理、裂隙十分发育，根据倾角大小，可分为中陡倾角结构面和缓倾角结构面，岩体结构类型以镶嵌结构、似互层状结构为主[18]。

在右岸平硐PDZ14中采用精测网法测量了640条结构面，现场判断发育6组，其极点等密度图见图1-A，得到的优势产状见表1，不能直接给出各组结构面条数。在左岸坝址区平硐PDZ11中全洞段测量了813条缓倾角结构面，现场判断可以分为4组，其极点等密图见图1-B，产状变化范围很大，无法精细分组，得到的优势产状为330°∠11°。

图1 平硐中实测结构面的极点等密图Fig.1 The pole figures showing measurement data from the precision mesh method in adits

采用花授粉算法，对PDZ14的分析成果如图2所示，当分组数取6时，全局Silhouette指标取得最大值，即最佳分组数为6组，优势产状及各组结构面数量见表1，优势产状与极点等密图法结果高度一致；对PDZ11中缓倾角结构面的分析成果见图3，Silhouette指标最大值对应的分组数为4，即最佳分组数为4，优势产状分别为191.8°∠14.6°、278.9°∠30.2°、346.0°∠14.7°、84.6°∠26.7°，这与现场判断确定的测量结果高度一致[18]。