随着我国对资源需求的日益增加和浅层矿产资源的逐渐枯竭，深部开采已经成为矿业发展的一种趋势[1]．岩爆是储存在围岩中的弹性能量突然释放导致岩石碎片猛烈弹出的动态失稳现象，具有突发性和极强的破坏性．岩爆的发生一般与围岩应力状态、岩性、地质构造、地形等因素相关．在开采作业期间，岩爆的发生会威胁开采人员的生命安全，导致支护结构的倒塌或变形、损坏设备，以及使运营单位承受经济损失[2-5]．岩爆灾害已成为我国深部金属矿山开采的重要制约因素，因此对岩爆灾害预测和预防的研究刻不容缓，这对安全高效地开采矿产资源具有重要的意义．岩爆预测方法总体可分为四类．第一类是基于指标判据的经验方法，如Barton判据[6]、Hoek判据[7]和二郎山判据[8]等．第二类是利用设备在现场监测的方法，如微震法[9]和声发射法[10]等．第三类是应用数值模拟的方法，应用的主要软件有FLAC 2D/3D[11]，PFC 2D/3D[12]和RFPA[13]等．第四类是应用数学的方法，又可分为两类：a. 基于指标的不确定性方法，如灰色系统理论[14]、多维云模型[15]、D-S证据理论[16]等；b. 基于数据驱动的机器学习方法，如支持向量机[17]、随机森林[18]和概率神经网络[19]等．目前，随着机器学习领域的快速发展和岩爆数据量的增多，岩爆预测的研究朝着更加智能化的方向发展，大量研究者已将不同类型的机器学习模型应用在岩爆预测上，并取得了较好的预测效果[20]．这些机器学习方法可大致分为两大类，即单一模型和集成模型(同质或异质)．其中单一模型存在鲁棒性低和不能对所有问题都得到最优解的固有缺点，而用装袋、提升和堆叠集成策略构建的集成模型可以克服单一模型存在的缺点．并且在这三种集成策略中，通过提升策略构建的集成树模型是当前机器学习领域非常流行的算法．尤其是在梯度提升树(GDBT)的基础上进一步优化后的极限梯度提升(XGBoost)算法，具有更高的精度和灵活性．文献[21]分析了XGBoost对缺失值容忍度，结果表明其不对缺失值填补也能保证较高的预测精度，且对缺失值容忍度较高．文献[22]利用交叉验证训练XGBoost岩爆预测模型，并对其适用性进行讨论．文献[23]对样本进行加权处理后建立CRITIC-XGBoost岩爆预测模型，与其他模型相比有更高的准确率．上述XGBoost模型虽然取得一定的预测效果，但是存在超参数的随机选择、训练数据量较少(40)和类不平衡的缺陷．超参数的随机选择通常会影响模型预测精度和收敛速度，数据量过少则会导致模型学习的信息量不足，此外类不平衡会导致模型只在多数类有较好的表现．本研究提出了一种基于非洲秃鹰优化算法(African vultures optimization algorithm，AVOA)的XGBoost岩爆预测模型，并引入合成少数类过采样技术(SMOTE)来解决以上模型存在的不足和缺陷．首先，收集国内外岩爆案例建立岩爆数据库，进行指标的筛选及对少数类数据的扩充；然后，建立AVOA-XGBoost岩爆预测模型，通过测试集验证AVOA-XGBoost模型的泛化能力；最后，将该模型与其他机器学习模型比较分析，并应用于三山岛金矿工程岩爆案例．1 数据收集和预处理选择合理有效的评价指标对岩爆预测至关重要，本研究初步选取最大切向应力σθ、单轴抗压强度σc、单轴抗拉强度σt、应力系数σθ/σc、脆性系数σc/σt和弹性能量指数Wet共6个指标．这些评价指标从不同角度反映岩爆发生的主要影响因素，且易于获取．岩爆烈度根据岩爆发生规模、破坏强度(潜在)和破坏模式，一般可分为无岩爆、轻微岩爆、中等岩爆和强烈岩爆．根据初步选取的评价指标和确定的分级标准，广泛收集国内外矿山、隧道及水电站等地下开挖工程岩爆案例．本研究参考文献[17，19，24]建立的岩爆案例数据库，剔除重复和含缺失值的样本后建立了包含326个岩爆案例的数据库．该数据集中包含53个无岩爆样本、99个轻微岩爆样本、119个中等岩爆样本和55个强烈岩爆样本．1.1　相关性分析及指标筛选当输入指标过多及指标之间相关性较强时，会增加模型构建的复杂性和模型训练的时长，甚至影响模型的预测效果；当输入指标过少时，则不能充分反映与岩爆发生相关的影响因素．因此，本研究通过计算皮尔逊(Pearson)相关系数来评估6个评价指标之间的相关性，计算公式为r=∑(x1i-x¯1)(x2i-x¯2)∑(x1i-x¯1)2∑(x2i-x¯2)2，式中：x¯1为指标x1的平均值；x¯2为指标x2的平均值；r为指标x1和x2之间的皮尔逊相关系数．由图1可知：只有σθ与σθ/σc和σt与σc/σt之间的皮尔逊相关性系数绝对值在0.6以上，才有较强的相关关系；这表明有些评价指标之间存在冗余信息，且须要剔除某些指标．因此，本研究采取Boruta算法对评价指标进一步筛选，Boruta算法计算原理参见文献[25]．在Python平台上利用10.13245/j.hust.231269.F001图16个评价指标间的相关性热图boruta_py库实现该算法，得到评价指标的筛选结果，其中σθ，σc，σt，σc/σt和Wet为“保留”，只有σθ/σc的状态是“剔除”，故将其从评价指标体系中剔除．1.2　不平衡数据集的处理数据集中明显存在类不平衡的问题，故模型在训练中会过多关注多数类，从而导致模型在少数类上的训练效果不佳．因此，本研究采用合成少数类过采样技术[26]对无岩爆和强烈岩爆两类进行扩充．无岩爆和强烈岩爆的采样倍率都为2，每个少数类样本都与选择的近邻样本利用下式合成新的样本，xnew=x+rand(0,1)x-xk，式中：xnew为新合成的少数类样本；x为少数类样本；rand(0，1)为(0，1)的随机数；xk为第k个近邻样本．最后，数据集中的无岩爆和强烈岩爆的样本数量分别被扩充到106个和110个．为了消除数据集量纲的差异和减少模型计算的成本，利用下式对扩充后的数据集进行标准化处理，Xij=(xij-xmini)/(xmaxi-xmini)，式中：Xij为第i个指标中第j个样本的标准化值；xij为第i个指标中第j个样本的实际值；xmini为第i个指标的最小值；xmaxi为第i个指标的最大值．2 模型和方法2.1　XGBoost算法XGBoost是一种基于树提升的可扩展机器学习系统，由一组训练的分类和回归树(classification and regression tree，CART)组成的集成模型[23]．除开始构建的分类和回归树外，后续的每个分类和回归树都会去拟合上一次预测的残差，具体原理如下．假设数据集D中由n个样本和m个特征组成(D∈Nn×m(R))，分类和回归树集成模型通过K个基模型函数求和来预测输出，y^i=ϕ(xi)=∑k=1Kfk(xi)  (fk∈F)，式中：y^i为输出结果；fk为第k个基模型分类器，每个fk都与一个独立的树结构q和叶子分数w相关；F为分类和回归树的空间．然后，利用对应的决策规则对叶子进行分类，并对相关叶子中的分数求和来计算最终的预测输出．基模型分类器fk通过最小化目标函数来学习，其中目标函数Lj*加入正则化项来避免过拟合问题，泰勒展开的二次项用来提高迭代效率，有Lj*=-12∑j=1T(∑i∈Ijgi)2/(∑i∈Ijhi+λ)+γT，式中：L为损失函数；T为叶子节点个数；Ij为结果中包含叶子节点j的实例集；gi和hi分别为损失函数的一阶和二阶梯度统计值；λ和γ为正则项系数．假设IL和IR为分割后左右节点的实例集(I=IL⋃IR)，那么叶子分裂的评分通过下式来计算，并作为决定最佳叶子分裂点的评判标准，有    Lsplit =∑i∈ILgi2/∑i∈ILhi+λ+∑i∈IRgi2/∑i∈IRhi+λ+∑i∈Igi2/∑i∈Ihi+λ /2-γ.2.2　非洲秃鹰优化算法非洲秃鹰优化算法是一种模仿非洲秃鹰觅食和导航行为的元启发式优化算法，有灵活性高、迭代速度快、能较快跳出局部最优解等优点[27]．该优化算法由确定种群中最好的秃鹰、计算秃鹰的饥饿率、探索阶段和开发阶段四个部分构成，具体流程如下．a. 设置相关参数和初始化种群空间．b. 计算种群个体的适应度，确定种群中排名第一和第二的秃鹰．c. 选择个体的移动方向(排名第一的秃鹰或第二的秃鹰)，R(i)=V1B  (Pi=L1);V2B  (Pi=L2),式中：R(i)为当前种群个体选择的最好的秃鹰；V1B为当前排名第一的秃鹰；V2B为当前排名第二的秃鹰；Pi为当前种群个体的选择概率(通过轮盘赌的方式确定)；L1和L2分别为选择排名第一和第二秃鹰的概率参数．d. 计算当前种群个体的饥饿度，F=[2×rand(0,1)+1]z(1-i/n)+t;t=hsinwπ2×in+cosπ2×in-1,式中：F为秃鹰的饥饿度；rand(0,1)为在(0，1)的随机数；z为在[-1，1]的随机数；i为当前的迭代次数；n为总的迭代次数；h为在[-2，2]的随机数；w为判断勘探和开采阶段中断的参数．e. 根据饥饿度来判别个体所进入的阶段(探索或开发)，然后分别采取不同的位置更新策略．f. 重复b～e，直至满足算法最大迭代次数．3 结果与讨论3.1　性能评价指标岩爆智能预测模型建立后，须要选取合适的评估指标来衡量模型的预测性能和泛化能力．本研究选择了4个类内分类模型评估指标，即准确率(ACC)、精确率(P)、召回率(R)和F1得分(F1)来衡量模型的性能．除此之外，还选择Kappa系数(κ)作为全局分类性能指标．3.2　模型构建XGBoost算法须要调整的超参数较多，手动调整是一个非常繁琐的过程，且不能保证模型的预测效果，故本研究用非洲秃鹰优化算法来快速有效地寻找超参数的最优值，具体过程如下．步骤1 将经过预处理的数据集划分为训练集和测试集．本研究采用分层抽样的方法，依次对每个岩爆类别的样本随机抽取20%，并将抽样结果放在一起作为测试集，而每个岩爆类别中剩余80%样本放在一起作为训练集．步骤2 确定XGBoost算法中须要调优的参数．由于XGBoost须要确定的超参数较多，因此本研究只对其中几个关键的超参数调优，分别为基分类器的数量(n_estimators)、学习率(learning_rate)、叶子节点中最小的样本权重和(min_child_weight)、CART的最大深度(max_depth)、最小损失减少(gamma)、二次采样(subsample)，其他参数取初始设定的默认值．在参考XGBoost算法对以上超参数的定义及取值范围后，最终确定的寻优区间见表1．10.13245/j.hust.231269.T001表1XGBoost算法超参数寻优范围超参数下限上限类型n_estimators100400离散learning_rate01连续gamma01连续max_depth015离散subsample01连续min_child_weight0.41连续步骤3 设置非洲秃鹰优化算法中的参数．该算法须要确定的参数主要有迭代次数N、种群的数量Pop、选择机制中的概率和确定勘探和开采阶段中断的参数．本研究对非洲秃鹰优化算法中存在参数的取值不进一步分析，取文献[27]中的参考值：N=100，Pop=20，L1=0.8，L2=0.2，w=2.5，P1=0.6，P2=0.4，P3=0.6．步骤4 计算适应度值．本研究将模型的准确率作为适应度值，并结合5折交叉验证进行计算．为了与其他群优化算法对比，利用粒子群优化(PSO)算法和非洲秃鹰优化算法一起对XGBoost算法的超参数寻优，迭代过程如图2所示．从图2中可以发现非洲秃鹰优化算法的寻优效果和速度都要优于粒子群优化算法．其中，非洲秃鹰优化算法确定的超参数组合为n_estimators(=212)，learning_rate(=0.12)，min_child_weight(=0.73)，gamma(=0.18)，max_depth(=4)，subsample(=0.7)．10.13245/j.hust.231269.F002图2适应度迭代曲线步骤5 建立AVOA-XGBoost岩爆烈度预测模型．在训练集上，利用步骤4确定的超参数组合训练XGBoost模型分类器．3.3　模型分类性能分析与比较将测试集输入到AVOA-XGBoost岩爆烈度预测分类器，得到该分类器在测试集上的预测结果(见图3)，其中无岩爆案例全部预测正确，其他三类的岩爆案例也只有少数被预测错误．计算得到准确率和Kappa系数分别为94.25%和0.92，从表2可以看出模型在各个类别的精确率和召回率都在0.8以上，表明AVOA-XGBoost岩爆烈度预测模型有着较优的预测性能，且泛化能力优异．10.13245/j.hust.231269.F003图3模型在测试集上的混淆矩阵10.13245/j.hust.231269.T002表2AVOA-XGBoost模型性能评估评估指标岩爆类别1234P0.951.00.920.91R1.000.900.920.95F10.980.950.920.93为了进一步验证本文模型的有效性，将AVOA-XGBoost分类器与经典的分类器算进行比较(见表3)，结果表明：本文模型的准确率最高，且比单独的XGBboost模型高3.45%，证明AVOA-XGBoost是一种有效的岩爆预测方法．10.13245/j.hust.231269.T003表3AVOA-XGBoost与经典模型比较预测模型准确率/%AVOA-XGBoost94.25XGBoost90.80支持向量机82.76随机森林86.21逻辑回归80.46注：经典模型的超参数利用交叉验证和网格搜索来确定．图4为AVOA-XGBoost和XGBoost模型的训练误差迭代曲线，从中可以看出AVOA-XGBoost的训练精度和收敛速度都有更优的表现．10.13245/j.hust.231269.F004图4模型训练迭代曲线3.4　特征变量的重要性分析本研究提出的AVOA-XGBoost模型是基于XGBoost算法构建的，可以统计每个特征在所有分类和回归树中作为“划分特征”的频数来衡量输入变量的重要性[23]．由图5可知：Wet和σc /σt的特征重要性最大，分别为112和98；这表明岩爆受Wet的影响最大，且通常发生在累积高弹性应变能的脆性岩体中，与已有研究成果[18，28]相符；而σc的特征重要性最小(69)，这表明σc对构建本模型的贡献是最小的．10.13245/j.hust.231269.F005图55个评价指标的重要性得分4 工程验证三山岛金矿是我国典型的海底深部开采金属矿山，其矿体主要储存于地下1 000～2 000 m之间．当该金矿开采深度达到915 m时，发生岩爆灾害，并在960 m处发生较严重的岩爆事件[29]．文献[28]通过对三山岛金矿的8个地区采样对岩爆灾害进行评估．本研究将该工程实例输入到AVOA-XGBoost模型分类器，得到预测结果见表4，由表4可知全部岩爆案例被预测正确，模型表现出优越的工程实用性．10.13245/j.hust.231269.T004表4三山岛金矿岩爆案例预测结果序号σθσcσtσc/σtWet真实值预测值128.4059.658.676.882.1222236.5890.4714.886.082.1122355.52115.6512.159.521.9022471.15131.4814.918.825.4633586.93165.2313.4412.299.10446113.56115.659.4612.233.9944736.46110.0112.948.504.5344897.0697.5111.938.174.9344为了方便现场工作人员使用本研究提出的AVOA-XGBoost预测模型，在建立岩爆烈度预测模型后，通过Matlab的App Designer模块设计模型预测平台系统的登录界面和预测界面，并将代码嵌入到可视化交互界面中．预测平台建立完成后，工作人员可在登录界面通过输入用户名和密码进入预测界面．然后，在预测界面中输入5个评价指标的数值，并点击运行按钮，则该平台预测结果对应的警示灯将变成红色．5 结论本研究提出了一种结合非洲秃鹰优化算法和XGBoost算法的新型混合机器学习模型，并将其用于深埋地下工程的岩爆烈度预测，得到的结论如下．a. 根据初步选取的6个评价指标建立了326个样本的初始数据集，由于个别指标之间存在较强的相关性，采用Boruta算法剔除了应力系数指标．最后，利用合成少数类过采样技术解决了数据库中存在的类不平衡问题，数据集被扩充到434个．b. 采用非洲秃鹰优化算法可以对XGBoost模型的超参数寻优，且效果优于粒子群优化算法．AVOA-XGBoost在测试集上的准确率为94.25%，Kappa系数为0.92，具有较优异的泛化能力．同时，也表明非洲秃鹰优化算法可以进一步应用于其他智能模型超参数的确定．c. 将AVOA-XGBoost与XGBoost、支持向量机、随机森林和逻辑回归比较，发现AVOA-XGBoost模型的准确率最高；且与XGBoost相比，其有更优的收敛速度，这也进一步说明本研究提出的模型有较高的优越性和可靠性．特征重要性分析结果表明Wet是岩爆烈度预测模型中贡献最大的特征．d. 将AVOA-XGBoost模型应用于三山岛金矿工程岩爆案例，其预测效果优异，具有较好的工程实用性；同时，可将该工程验证数据用于扩展岩爆数据库，以便后续研究使用．