基于反赋权与MBCT-SR多维云模型算法岩爆预测研究*

宋英华，庞昭胜，李墨潇，江 晨，齐 石

(1.武汉理工大学 中国应急管理研究中心，湖北 武汉 430070；2.武汉理工大学 安全科学与应急管理学院，湖北 武汉 430070)

0 引言

岩爆是坚硬岩石开采和土木建筑中最常见的由连续岩体应力过大引起的破坏之一，其发生总是伴随着岩块开裂、剥落、抛掷、巨响等现象[1]。岩爆具有突发性、破坏性强、难控制等特点，极易对现场施工人员及工程设备造成损伤[2]。随着地下岩体开采深度和挖掘需求的增加，岩爆事故的发生频率与影响程度也越来越严重，岩爆已成为愈发被重视的工程问题。

由于岩爆机理的复杂性，且各指标数据测取具有较大的不稳定性，岩爆分类预测在世界范围内都一直是迫切需要解决的难题，岩爆烈度等级评价研究主要有3大类[3]：1)第1类是基于岩爆机理判据，如应力强度Russense判据、Barton判据，能量理论中能量比指标判据等直接进行评价；2)第2类则需要对岩爆现场进行实测，诸如微震法、声发射法等；3)第3类是基于岩爆影响因素的综合预测方法，是目前岩爆预测研究的热点。第3类方法又有2种不同的判别方式，一是基于岩爆工程实例数据，如采用XGBoost[4]、神经网络[5]、随机森林[6]等机器学习模型进行评判；二是基于岩爆指标判据的预测方法，主要运用模糊综合评判模型[7]，理想点模型[8]，云模型[9]等进行预测。

针对岩爆的随机性和模糊性特点，云模型的使用在岩爆预测中成为当下研究热点，已被证明具有一定的可靠性。Liu等[10]将云模型与粗糙集理论进行结合，根据冲击地压标准生成正态云模型，运用粗糙集理论确定权重，最后利用最大隶属度原理确定岩爆等级；周东良等[11]引进改进的AHP、熵权法、博弈论和模糊熵理论组合赋权法与二维云模型综合评判；刘晓悦等[12]将指标的权重融合到多维云模型中，生成多指标的等级综合云进行预测并验证了模型的可靠性及实用性。目前，云模型在岩爆预测方面的研究中，大多学者都采用正态云模型的正向云发生器算法，而正向云发生器算法的数字特征一般是经验值，往往具有较强的主观性。逆向云发生器算法是基于数据生成数字特征，具有较强的客观性，在岩爆预测研究中的使用较少。由于岩爆灾害等级确定受多指标综合影响，多维云模型也应运而生。多维云模型的权重分配研究目前大多是依靠主观赋权与客观赋权相结合的方式，但此类方式往往面临着主观权重的精度以及权重结合方法的科学性不足等问题，使问题复杂化。

本文采用逆向云发生器MBCT-SR算法[13]，解决确定云模型数字特征及权重时，由于主观干扰导致的预测结果与实际情况存在偏差等问题，该算法基于一定数量的样本实例，计算客观云数字特征，并结合多维云模型建立动态适应度函数，通过改进的遗传算法反求最优权重，从而建立完整的岩爆预测云模型，为岩爆预测研究提供更贴合实际案例数据的评价方法，使评价结果更加客观准确。

1 岩爆预测多维云模型

1.1 多维正态云模型

正态云模型[14](图1)可实现定性概念与定量表示的相互转换，反映定性概念的不确定性，其性质通过期望Ex、熵En以及超熵He3个数字特征表示。其中期望Ex为云的重心，最能代表定性概念；熵En表示定性概念的离散程度；超熵He则表示熵的离散程度，是熵的不确定性度量。

图1 云的数字特征Fig.1 Digital features of cloud model

在一维正态云模型的基础上，引入精确数值表示的n维集合U={X1,X2,X3,…,Xn}，其中U为n维定量论域，C为U上一定性概念。对于U中任意值X=(x1,x2,x3,…,xn)(X∈U)，都存在X对U有一定约束的随机数μ(X(x1,x2,x3,…,xn))，称为确定度。多维正态云模型确定度式[15]为式(1)：

(1)

式中：i表示n维数据中第i维数据(i=1,2,3,…,n)；yi为服从以熵Eni为均值，超熵的平方Hei2为方差的正态分布随机数；xi满足以期望Exi为均值；yi为方差的正态分布规律。

多维云模型包含2种云发生器，一种是正向云发生器(CGn)，另一种为逆向云发生器(CGn-1)，2种发生器可进行双向转换，如图2所示。

图2 多维云发生器Fig.2 Multi-dimensional cloud generator

n维正向云发生器是指将定性概念通过n组数字特征N(Exn,Enn,Hen)表示，并生成一定数量的云滴P(X(x1,x2,x3,…,xn))，μj的过程，公式(1)中用j来表示岩爆等级序列(j=1,2,3,4)；逆向云发生器是指基于一定数量的云滴，计算云数字特征的过程，在云滴较少的情况下数字特征一般为估计值。

1.2 MBCT-SR算法及综合确定度

通过式(1)可知数字特征决定确定度的取值及其波动范围，传统云模型数值特征一般采用经验式[16]进行计算，即式(2)：

(2)

式中：各指标上下限Cmax、Cmin和固定值K均为经验常数，这使得云模型规避了岩爆预测的模糊性特征，导致主观性较强。云模型的逆向云发生器算法基于大量云滴数据，客观计算云模型数字特征，减小经验常数的主观性。但传统的SBCT-1stM逆向算法常伴随计算结果漂移及不稳定等现象，因此，本文采用多步还原的逆向云变换MBCT-SR算法，提高结果的准确度和稳定性。具体算法步骤如下：

Step 1：输入m个样本点xk，k表示m个样本点中所取样本序列(k=1,2,3,…,m)。

Step 2：计算样本均值，得到期望Ex的估计值，即

(3)

(4)

Step 4：将Step 3中的结果代入式(5)计算出En2与He2的估计值，即式(5)：

(5)

多维模型的建立需结合各指标权重，由式(6)计算综合确定度Ω[17]，其中i为n维各项数据序列(i=1,2,3,…,n)，最终达到预测的目的。

(6)

2 指标权重反分析

依据文献收集[3,6,11,12,16-18]的真实样本数据192组，采用遗传算法对指标权重进行反分析。设指标权重为未知变量，以最大化满足样本数据为结果，通过优化算法求全局最优解。具体流程见图3。此方法降低权重确定主观性，其分析步骤为：①确定评价指标；②收集样本实例并进行数据预处理；③计算数字特征，构建多维云模型；④由式(6)，建立优化适应度函数；⑤利用遗传算法，寻求全局最优解；⑥回代，验证完整多维云模型的有效性及可行性。

图3 指标权重反分析方法流程Fig.3 Flow chart of index weight back analysis method

遗传算法只能从适应度函数中获取信息，故步骤④为反求权重的核心内容，在求取适应度函数时，由式(6)可得综合确定度中的参数yi为服从以熵Eni为均值，超熵的平方Hei2为方差的正态分布随机数(i=1,2,3,…,n)，故yi需由MBCT-SR算法计算获得其正态分布均值与方差，再进行随机取值，最终建立动态适应度函数fitness[18]如式(7)：

(7)

式中：j为岩爆等级序列；i为评价指标序列；k为样本序号；Ωj为第j个岩爆等级的综合确定度；n，p，m分别为指标总数、岩爆等级总数和样本总数；qk与Qk分别是第k个样本中的预测等级和实际等级。

由式(7)可知，优化算法以预测结果满足实例结果最大化为目标，规避主观干扰，建立适应度函数，求出最优权重。求得权重代入式(6)，根据综合确定度获得最终评价结果。

3 模型实例验证

3.1 评价指标选取及数据预处理

岩爆等级受多因素综合影响，本文基于真实数据分析，规避主观因素，对数据样本具有较高的要求，数据精度直接影响预测结果。结合以往研究经验，通过文献收集[3,6,11,12,16-21]的方法，最终确定数据样本较多、影响程度较大的3项指标：应力比Ts=σθ/σt、岩石脆性指数B=σc/σt以及弹性应变储能指数Wet作为岩爆等级预测评价指标，参考王元汉等[22]的相关研究，岩爆等级依托于各指标分为无岩爆(Ⅰ)、轻微岩爆(Ⅱ)、中等岩爆(Ⅲ)、强岩爆(Ⅳ)4个等级评价。

本文研究数据样本较全，不存在数据缺失的情况，但数据收集时无法避免存在噪点的情况出现。因此采用偏差分析结合箱线图进行预处理，采用aσ(a=2或3)规则[23]进行修剪，其中σ表示原始数据方差，a表示以a倍方差为临界值对离散数据进行截取。最终选取192 组国内外实际案例进行预测，其分布规律如图4所示。

图4 原始数据样本箱线Fig.4 Box line diagram of original data sample

从图4所知，由于各指标量纲不同，各指标数值差异较大，不利于进行科学计算及综合分析，故按式(8)对各指标进行归一化处理，即得

(8)

式中：xik为第i个指标中第k个原始样本数据；ximin，ximax分别是第i个指标中原始样本的最小值和最大值。

式(8)中原始数据归一化后其分布结果如图5所示。由图5可得，预测指标岩石脆性指数B=σc/σt依旧存在一定数量的离散点，数据中心偏移较为严重，而应力比Ts=σθ/σt以及弹性应变储能指数Wet分布较为均匀。指标上下四分位数跨度较小，其中Ts下四分位偏低，中位数靠近0.5，无离散点。弹性应变储能指数Wet上四分位数偏高，中位数稍低，但总体分布比较合理。由归一化数据分析初步猜测，本数据样本指标B稳定性较差，权重应最小，这也满足岩爆各指标权重分配的主观认知。

图5 归一化后样本箱线Fig.5 Normalized sample box plot

选取29组国内外实例数据做测试集，列举部分如表1所示，以剩余163组数据为训练样本，并与其他云模型预测方法对比，验证其准确性。

表1 国内外岩爆实例数据Table 1 Data of rock burst cases at home and abroad

3.2 计算数值特征

将163 组岩爆实例数据按级分类，并采用MBCT-SR算法对各级各指标岩爆数据分别计算，由于算法计算结果并不唯一，任选其一组将各数字特征示于表2 。

表2 数字特征计算结果Table 2 Digital feature calculation results

根据表2中各数字特征绘制多维云模型，如图6所示。从图6所知，“·”描绘数据代表无岩爆(Ⅰ级)，“+”描绘数据点代表弱岩爆(Ⅱ级)，“○”点代表中岩爆(Ⅲ级)，“*”点表示强岩爆(Ⅳ级)。指标Ts与实际岩爆顺序相同，相对于其他2个指标分布较为均匀，各等级界限较为明显。指标Wet的分布顺序与岩爆顺序也相同，但在无岩爆与弱岩爆中数据跨度较大，其分布特征相比于指标Ts较不稳定。B维度数据点分布顺序较乱，跨度较大，很明显相对于其他2指标规律性较差，最不适用于岩爆评价。此次MBCT-SR算法基于大量样本数据计算数字特征并生成三维云图的分析中，对第3.1节中假设进行了初步验证，并针对图6 做进一步假设，求得最终权重顺序应满足ωTs>ωWet>ωB。

图6 岩爆实例数据聚类云模型Fig.6 Clustering cloud model of rock burst case data

3.3 求取最优权重

将各数字特征代入式(7)构建适应度函数反求权重，在遗传算法中，适应度函数的每次使用将重新计算数字特征，确保结果真实性。

常规遗传算法迭代曲线保留原始样本最优解，寻求全局最优解，迭代曲线呈现单调上升趋势。本文由式(7)使用动态适应度函数，同一最优解会出现不同结果，满足岩爆评价实际情况。为降低时间成本，减少迭代次数，最快获取较为稳定全局最优解。本文增加种群数量，设置初始化种群数量N=1 000，迭代次数f=100，采用锦标赛选择和精英保留方法，增加图3中4个权重约束条件进行约束，运行Matlab最终获得迭代曲线如图7所示。

图7 迭代曲线Fig.7 Iteration curve

由图7得迭代48 次后满足条件样本数量保持稳定，此时各指标最优权重为ωTs=0.632 1、ωWet=0.257 7和ωB=0.110 2，最优权重与假设权重比例相符，满足样本实例。

3.4 实例验证与对比

将所求结果通过前文表1中未参与反求权重的29组实例样本进行评价。由于岩爆事故本身具有较强随机性，应将评价结果的合理波动考虑进来，增大预测范围，分析评价结果更多的可能性。本文对样本数据进行客观评价，采取保守预测方法，将确定度差值不足0.1的岩爆等级进行倾向性预测，其预测结果进行跨区间表示，例如：Ⅱ～Ⅲ级岩爆，此方法满足实际岩爆预测要求，使预测更加严谨，并与其他云模型评价方法进行结果对比，进一步验证本文方法的可靠性与有效性，具体情况见表3。

表3 岩爆倾向性评价结果Table 3 Evaluation index judgment interval

结果显示，预测结果与实际基本相符，与其他模型的预测结果基本吻合。本文与反赋权一维云模型权重确定方法相同。由表3得，反赋权一维云模型倾向性预测准确率为72%。组合赋权多维云模型预测结果准确率为86%，比一维正向云发生器的准确率更高。在逆向云发生器中，MBCT-SR算法与优化算法结合评价精准预测准确率可达89%，倾向性预测准确率可达100%。预测结果对比表明，本文评价模型更加合理且有效，算法充分考虑实际岩爆的随机性与模糊性，更贴合实际情况。

4 结论

1)结合当前机器学习的热门方向，依据国内外192 组岩爆实例数据，选取应力比Ts=σθ/σt、岩石脆性指数B=σc/σt以及弹性应变储能指数Wet作为评价指标，将逆向云发生器与优化算法相结合建立综合评价模型。本模型引进MBCT-SR算法，规避常规正向云发生器主观性过强、绝对性较大的缺点，建立更加客观的多维岩爆预测云模型。

2)指标权重优化结果来自于真实数据与多维云模型。本文依据图5、6中的数据可视化分析，对权重分布规律进行了满足ωTs>ωWet>ωB的初步假设，并由最终计算结果求得最优权重为ωTs=0.632 1、ωWet=0.257 7和ωB=0.110 2，对假设进行了科学验证。使得本指标权重确定方法更加准确严谨，最终结果与主观认知相符。

3)相比其他预测模型可得，本模型最大化降低主观因素干扰，在选取应力比Ts=σθ/σt、岩石脆性指数B=σc/σt以及弹性应变储能指数Wet作为评价指标进行预测时，精准预测率达89%，倾向性预测结果准确率可达到100%，其评价结果更加准确、真实。