露天矿数码电子雷管逐孔起爆条件下质点峰值振速预测

Blasting-Induced Peak Particle Velocity Prediction of Hole-by-Hole Blasting Operation Using Digital Electronic Detonator in Open-Pit Mine

摘要

针对目前露天矿爆破质点峰值振速预测研究存在模型可解释性不足、不适用于数码电子雷管逐孔起爆条件等问题，通过现场试验记录每孔爆破参数与测取爆破振动信号，结合轻型梯度提升机（LightGBM）算法与SHAP模型可解释性框架，建立了露天矿数码电子雷管逐孔起爆条件下的三轴质点峰值振速预测模型.从测试集均方根误差RMSE和拟合优度R ²而言，LightGBM总体RMSE相比于支持向量机与神经网络分别降低了25.9%和28.9%，总体R ²分别提高了12.7%和9.9%.LightGBM与萨道夫斯基经验公式相比，RMSE在径向X、切向Y和垂向Z上分别降低了63.4%、39.5%和68.3%， R ²分别提高了18.9%、27.7%和42.4%.除方向轴变量外，监测点距离、总药量、最小排距、平均装药长度、孔径与最大孔距为对质点峰值振速影响程度最大的6个变量，其中监测点距离与质点峰值振速为负相关关系，总药量、最小排距、平均装药高度与最大孔距则与质点峰值振速呈正相关关系.

关键词

采矿工程 / 逐孔起爆 / 爆破振速 / 机器学习 / 工程地质

Key words

mining engineering / hole-by-hole blasting operation / blasting vibration velocity / machine learning / engineering geology

0 引言

随着我国经济迅猛发展，各行各业对矿产资源的需求量与日俱增，露天开采作为重要的矿产资源获取方式，具有生产能力大、回采率高、作业安全、机械化及自动化程度高等优点，其中钻孔和爆破作业是露天矿山开采的主要工序和成本控制的重要因素.然而研究表明，只有30%的爆破能量直接用于岩石破碎，其余能量将以振动、飞石、噪声和空气冲击波等形式对周边环境造成危害（Trivedi et al.， 2014）.

在上述的爆破副作用中，振动因其危害范围广、程度大而一直被认为是首要危害形式.在现行的爆破安全规程GB 6722-2014中，以质点峰值振速（Peak Particle Velocity， PPV）作为地面振动幅度的评估指标并对其加以控制，大量的工程实践表明，毫秒延期爆破技术可有效降低爆破振速.毫秒延期爆破技术利用相邻炮孔先后起爆所产生的应力波叠加作用和新增自由面作用，在保证岩石的破碎效果下提高了能量利用率，同时减少了炸药使用量，从而达到降振效果（Tian et al.， 2019）.因此，毫秒延期爆破技术的关键在于爆破网络的设计以及延期时间的设定.理论上而言，短延时逐孔起爆可充分利用相邻应力波叠加作用，同时使爆破能量在空间中分散，从而使岩石进一步破碎，但导爆管雷管段别少、延期时间不精准的缺点制约了短延时逐孔起爆技术在露天矿中的应用.近年来，数码电子雷管因其精准的延期时间、优化简单的网络设计等优点得到快速发展，逐步取代了导爆管雷管并广泛应用于矿山开采和隧道建设，使得短延时逐孔起爆技术在露天矿中的应用逐渐增多.因此，对露天矿短延时逐孔起爆所引起的爆破振速进行预测并分析其影响因素具有重要的工程意义.

在目前对露天矿爆破振速影响因素和预测的研究中，研究者们通常使用三轴振动传感器采集径向、切向和垂向的爆破振速，采用经验公式、数值模拟、机器学习算法3种方法将质点峰值振速与爆破参数建立联系.经验公式因其简单易行而被广泛应用，其中最为著名的为萨道夫斯基公式和美国矿业局USBM公式，研究者们进一步通过量纲分析制定了多种修正经验公式（杨曦等， 2017； 何理等， 2021； 朱明等， 2021），这些经验公式对于初始爆破设计有借鉴意义，然而它们的拟合优度和精度较低；数值模拟在数值分析软件中构建现场岩体、炸药的模型并设定物理力学参数，建立动态本构方程，用数值算法对炸药爆炸全程进行模拟与求解（Zhou et al.， 2016； 张云鹏等， 2018； 张耿城等， 2020）.尽管数值算法的适用性广，可模拟岩体对爆炸冲击的非线性响应过程（Jiao et al.， 2021； Zhang et al.， 2021a），但数值模拟全流程所耗费的人力物力及时间成本过高，且所建立的模型难以及时修改参数应用于工程实践；机器学习算法相较于前两种方法而言，具有拟合优度高、精度高、鲁棒性强等优点（李文彬等， 2021； 郝慧珍等， 2021； Zhang et al.， 2021b； Zhang and Phoon， 2022），然而它的可解释性较差，需要一定量的数据建立模型（Hosseini et al.， 2019； Aladejare et al.， 2022）.上述研究多是选取单轴爆破振速数据进行建模，且在模型的可解释性上并未深挖.除此之外，上述研究多基于导爆管网络起爆进行试验，由于导爆管雷管的段别限制，单段起爆药量显著大于逐孔起爆的单孔药量，故并不适用于数码电子雷管逐孔起爆的情况.

本文依托于露天矿开采现场开展数码电子雷管逐孔起爆爆破试验，记录每孔爆破参数与爆破振动数据，结合构建变量、多次交叉验证选取模型超参数等方法，基于轻型梯度提升机（Light Gradient Boosting Machine，LightGBM）建立适用于露天矿数码电子雷管逐孔起爆条件下的三轴质点峰值振速预测模型.在此基础上，分别建立支持向量机（Support Vector Machine，SVM）和神经网络（Neural Network，NN）两种机器学习模型、并拟合常用经验公式，将其结果与LightGBM模型进行精度和拟合优度对比.进一步引入基于博弈论沙普利值（Shapley）概念的SHAP（Shapley Additive Explanations）模型可解释性框架，评估最优模型中爆破参数的总体重要性，分析爆破参数对质点峰值振速预测值的影响程度及方式.

1 建模方法及原理

本文主要基于LightGBM算法建立三轴质点峰值振速预测模型，为了与传统的经验公式进行比较，同时以回归方程的方式拟合萨道夫斯基公式和USBM公式.

1.1　经验公式拟合

露天矿中常用的爆破质点峰值振速预测经验公式主要为萨道夫斯基公式（简称萨氏公式）和美国矿业局USBM公式，其形式为：

v = K ( Q^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{3}$}\right.} {/R)}^a,(1)

v = K ( Q^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{2}$}\right.} {/R)}^a,(2)

其中，式（1）为萨氏公式，式（2）为USBM公式.式中：v为质点峰值振动速度，cm/s；Q为单段最大装药量，kg；R为监测点距离，m；K、α为与地形地貌相关的系数.

式（1）和式（2）两侧取对数，可得：

\mathrm{l} \mathrm{n} ( v ) = \mathrm{l} \mathrm{n} ( K ) + a \mathrm{l} \mathrm{n} ( Q^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{3}$}\right.} / R ),(3)

\mathrm{l} \mathrm{n} ( v ) = \mathrm{l} \mathrm{n} ( K ) + a \mathrm{l} \mathrm{n} ( Q^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{2}$}\right.} / R ),(4)

利用实测质点峰值振动速度，基于式（1）和式（2）进行最小二乘法拟合，可求出系数K与α.

1.2　轻型梯度提升机模型（LightGBM）

LightGBM是基于梯度提升决策树（Gradient Boosting Decision Tree， GBDT）的一种改进算法，它有更快的训练速度、较低的内存消耗和可处理海量数据的优点（Ke et al.， 2017）.

与常用的GBDT算法类似，LightGBM主要思想是利用决策树迭代训练以得到最优模型：

F_T ( x ) = {\displaystyle \sum _{i=\mathrm{1}}^T} f_i(x),(5)

式中： x为自变量，在本文中为爆破参数； f_i ( x )为第i个决策树； F_T ( x )为T个决策树的累加模型.

设本次迭代所获得的累加模型为 F_K ( x )，目标函数为 O b j_K ( x )，损失函数为 L ( y , F_K \left(x\right) )，对于质点峰值振速，可用L2损失函数，具体形式为：

L ( y , F_K ( x ) ) = {\displaystyle \sum _{i=\mathrm{1}}^N} (y_i-F_K(x_i{))}^{\mathrm{2}},(6)

O b j_K ( x ) = L ( y , F_K ( x ) ) + \mathrm{\Omega } ( f_K ( x ) ) + b,(7)

\mathrm{\Omega } ( f_K ( x ) ) = \gamma M + \mathrm{1}/\mathrm{2} \lambda {\displaystyle \sum _{j=\mathrm{1}}^M} w_j^{\mathrm{2}},(8)

式中： y为因变量，本文中为质点峰值振速； N为样本数量； \mathrm{\Omega } ( f_K ( x )为复杂度函数； M为决策树叶子节点数目；w_j 为第j片叶子的输出值； \gamma、 \lambda为约束系数.

在下一次迭代中，目标为寻找决策树 f_{K+\mathrm{1}} ( x )使得目标函数 O b j_{K+\mathrm{1}} ( x )最小：

F_{K+\mathrm{1}} ( x ) = F_K ( x ) + \epsilon f_{K+\mathrm{1}} ( x ),(9)

f_{K+\mathrm{1}} ( x ) = \mathrm{a} \mathrm{r} \mathrm{g} \mathrm{m} \mathrm{i} \mathrm{n} O b j_{K+\mathrm{1}},(10)

式中： \epsilon为学习率.

在优化求解的过程中，最耗时的步骤是在每一个叶子节点寻找合适的变量并获取最优划分点.假设在某一个叶子节点上二分裂为两个次节点，分别为左次节点（L）和右次节点（R），对于目标函数而言，分裂后的增益Gain为：

\begin{array}{l}Gain=(\mathrm{1}/\mathrm{2})[G_L^{\mathrm{2}}/(H_L+\lambda )+G_R^{\mathrm{2}}/(H_R+\lambda )-\\ (G_L+G_R{)}^{\mathrm{2}}/(H_L+H_R+\lambda )]-\gamma \end{array},(11)

式中：G、H分别为损失函数 L ( y , F_K \left(x\right) )一阶导和二阶导的样本累加值.

传统的GBDT和极端梯度提升（Extreme Gradient Boosting， XGBoost）算法需要遍历所有样本，在每一个叶子节点枚举所有可用变量的可能切分点，并采用按层生长（level-wise）策略，因此训练速度缓慢、内存消耗大.而LightGBM对原始的连续型变量采用直方图算法（Histogram-based Algorithm），将连续的变量值抽象为离散值并构建变量直方图，在此直方图基础上寻找可能切分点并计算分裂增益Gain.进一步地采用按叶生长（leaf-wise）策略，从当前所有叶子中寻找分裂增益Gain最大的一个叶子节点，然后继续分裂并以此循环.在分裂次数相同的情况下，按叶生长策略可获得更高的精度且计算速度更快.

除此之外，LightGBM修正了以往的梯度下降算法，提出单边梯度采样（Gradient-based One-side Sampling，GOSS）.GOSS算法通过削弱小梯度样本权重并提高大梯度样本权重以计算加权信息增益，其核心理念为大梯度样本对信息增益有更大影响.

1.3　模型评估指标

为评估不同模型与方法之间精度和拟合优度的差异，鉴于质点峰值振速为连续型变量，选取均方根误差RMSE（Root Mean Squared Error）和R ²作为评估指标.RMSE越接近于0则模型误差越小，R ²越接近于1则代表模型拟合程度越高，其定义如下：

\mathrm{R} \mathrm{M} \mathrm{S} \mathrm{E} = \sqrt[]{\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}(y_i-{\widehat{y}}_i{)}^{\mathrm{2}}},(12)

R^{\mathrm{2}} = \mathrm{1} - \frac{{\displaystyle \sum _{i=\mathrm{1}}^N}({\widehat{y}}_i-y_i{)}^{\mathrm{2}}}{{\displaystyle \sum _{i=\mathrm{1}}^N}(\overline{y}-y_i{)}^{\mathrm{2}}},(13)

式中：y_i 表示峰值振速真实值； \widehat{y_i}为模型预测值；N为样本数目.

2 模型可解释性

一般而言，模型的可解释性有两种含义.一为可判断模型中变量总体重要性；二为对于特定样本，可分析预测结果的输出过程.LightGBM模型尽管有着精度高、速度快等优点，可实质上是一个决策树累加模型，在可解释性上和神经网络都属于黑箱模型，而对于黑箱模型而言，要实现上述两种可解释性较为困难.为了打开黑箱模型、实现可解释性，Lundberg提出了一种基于博弈论Shapley概念的SHAP模型可解释性框架（Lundberg and Lee， 2017； Lundberg et al.， 2018）.

SHAP定义一个线性累加形式的解释性模型g（z）以解释原模型f（x）：

g ( z ) = {\phi }_{\mathrm{0}} + {\sum }_M {\phi }_j{z}_j,(14)

x = h_x ( x^{\text{'}} ),(15)

g ( z ) \approx f ( h_x ( z ) ) , z \approx x^{\text{'}},(16)

式中：x为输入变量向量；M为输入变量个数；z∈｛0，1｝ ^M，表示相应变量是否存在；x′为M维的简化变量；φ_j 为变量X_j 的Shapley值.由于LightGBM本质上是决策树累加模型，对于每一个样本所有变量均存在，故式（21）可简化为：

g ( z ) = {\phi }_{\mathrm{0}} + {\sum }_M {\phi }_j,(17)

若原模型f（x）是线性模型，则定义变量X_j 对预测值的贡献为：

{\phi }_j ( \widehat{f} ) = {\beta }_j x_j - E ( {\beta }_j X_j ) = {\beta }_j x_j - {\beta }_j E ( X_j ),(18)

式中：β_j 为变量X_j 的系数；E（β_jX_j ）为变量X_j 的平均影响估计.

将样本x的所有变量贡献相加，可得：

\begin{array}{l}{\displaystyle \sum _{j=\mathrm{1}}^p}{\phi }_j(\widehat{f})={\displaystyle \sum _{j=\mathrm{1}}^p}({\beta }_j{x}_j-E({\beta }_j{X}_j))\\ =\widehat{f}(x)-E(\widehat{f}(X)).\end{array},(19)

从式（26）可知，对于样本x，所有变量贡献之和为模型预测值减去模型平均预测值.但上述情况仅在模型为线性模型且变量之间相互独立的条件下成立，对于LightGBM这种决策树累加模型，需对每一个变量在考虑所有可能的变量组合条件下计算Shapley值，即：

\begin{array}{l}{\phi }_j=\sum _{S\subseteq F\backslash \{j\}}\frac{\left|S\right|!(\left|F\right|-\left|S\right|-\mathrm{1})!}{\left|F\right|!}[f_{S\cup \{j\}}(x_{S\cup \{j\}})-\\ f_S(x_S)]\end{array},(20)

式中：S为变量子集；F为变量全集；f_S∪｛j｝ （x_S∪｛j｝ ）表示将变量X_j 与子集S的变量结合并训练模型后所得预测值；f_s （x_s ）为子集S的变量训练模型后所得预测值.

对于任意样本，均可分别计算各个变量的Shapley值，通过Shapley值直接分析任意样本预测值的输出过程.进一步地，为了衡量变量总体重要性，可对各变量的样本Shapley绝对值取平均值，将其定义为变量总体重要度：

I_j = {\displaystyle \sum _{i=\mathrm{1}}^N} \left|{\phi }_j^{(i)}\right|,(21)

式中：N为样本数量；φ_j ^（i） 为第i个样本变量X_j 的Shapley值；I_j 为变量X_j 的总体重要度.

由此可见，借助SHAP模型可解释性框架，既可判断变量总体重要性，又可分析特定样本模型预测的输出过程，从而实现两种可解释性.

3 露天矿爆破振速数据采集

3.1　工程背景

本研究依托于陕西省某露天煤矿进行试验，此煤矿地表岩层以粉砂岩、泥岩为主，煤质属特低硫、低灰、特高发热量的优质长焰煤.破岩使用露天深孔台阶爆破方式，钻孔设备采用牙轮钻机和潜孔钻机，孔径有250 mm、165 mm和150 mm三种.台阶高度为5.6~18.6 m，设计孔网参数为17×9.5 m、14×9.5 m和8×4.5 m.炮孔内采用连续装药结构，使用混装乳化炸药.起爆网路以数码电子雷管连接，采用逐孔起爆形式，孔间延期36 ms，排间延期100 ms，起爆网路示意图见图1，数字为起爆时间，单位为ms.

3.2　现场试验

本研究选用中科测控TC-4850爆破测振仪搭配三维振动速度传感器进行试验以测量三轴爆破振速，测振仪的频响范围为5~500 Hz.本试验采集频率设定为2 000 Hz.为扩充样本量及测定保护建筑物附近的爆破振速，在宿舍楼、办公楼、锅炉房等建筑物周边和爆破现场内布置监测点.单次爆破试验最少布置5个监测点.本试验共记录了22次爆破试验、153个监测点、459条爆破振速数据，爆破质点峰值振速记录例表如表1所示.监测点距离范围为32.6~2 432 m，爆破质点峰值振速范围为0.05~12.83 cm/s.在钻孔及装药过程中，精确测量与记录每次爆破的每孔装药量、台阶高度、孔排距等爆破参数.表2为爆破参数记录例表.

4 模型验证与解释

4.1　数据处理

鉴于爆破参数是每孔记录，故需对单次爆破下的爆破参数进行处理.本研究采取以下措施：

（1） 定义当次爆破下的炮孔总数为孔数，并取孔径为变量.

（2） 对台阶高度、孔深、孔距、排距、装药高度和填塞高度，分别取当次爆破下的最大值、平均值和最小值.

（3） 对于药量，除上述3个值外，另计算当次爆破总药量.

（4） 对坐标轴变量进行独热编码，形成3个坐标轴指示变量.径向X为100，切向Y为010，垂向Z为001.

（5） 将连续性变量如台阶高度、孔深等做标准化处理，并将质点峰值振速进行对数变换.

将上述标准化爆破参数变量与爆破质点峰值振速记录相结合，形成以对数质点峰值振速为输出变量、标准化爆破参数为输入变量的数据集.此数据集的维度为459×29，其中459为样本数量，29为包含对数质点峰值振速的变量个数，表3为原始变量及模型变量列表.

为客观验证模型的精确度和拟合优度，从总体数据集中随机抽取20%作为测试集，余下80%数据作为训练集，训练集中包含367个样本，测试集中包含92个样本.LightGBM、SVM和NN模型调参及评估时，将预测的质点峰值振速转化为原尺度后再计算RMSE和R ².

4.2　经验公式拟合结果

式（1）和式（2）中所用药量为最大单段药量，在逐孔起爆情况下应进行修正.为进行比对，分别取单孔药量最大值和单次爆破总药量代入经验公式，并划分坐标轴数据分别进行拟合.经验公式拟合后测试集RMSE与R ²对比如图2所示.通过对比可知，在三方向上取单次爆破总药量代入萨氏公式可得较低的RMSE值和较高的拟合优度R ²：在径向X上，RMSE和R ²分别为1.72与0.74；在切向Y上分别为1.85与0.42；在垂向Z上分别为2.18与0.66.取单次爆破总药量代入的USBM公式仅在径向X上与萨氏公式相近，而两种经验公式均在切向Y上拟合优度较低，在垂向Z上误差较大.

4.3　LightGBM模型结果

LightGBM算法中包含多个超参数，如学习率、约束系数 \gamma和 \lambda、决策树最大深度等，这些超参数对最终模型的预测能力有很大影响.此外，为验证LightGBM模型算法的优越性，分别建立SVM和NN模型与其对比.为准确选取各模型的超参数并降低模型过拟合程度，采用多次交叉验证（He and Chalise， 2020）（Repeated K-fold Cross Validation）确定最终模型的超参数.表4为LightGBM一百次五折交叉验证后RMSE最低的超参数组，表5和表6分别为SVM和NN的超参数组.

基于上述最优超参数训练时，可分别计算3种机器学习模型所需计算时长.LightGBM、SVM和NN所需训练时长分别为0.062 s、 0.073 s和0.498 s，LightGBM的训练速度最快.各模型在训练集与测试集上的总体RMSE和R ²对比如图3所示.LightGBM在训练集与测试集上总体RMSE及R ²均显著优于SVM、NN和萨氏公式.在测试集上LightGBM总体RMSE较之于其他3个模型分别降低了25.9%、28.9%和55.2%，总体R ²分别提高了12.7%、9.9%和28.9%.图4为测试集模型真实值与预测值散点图，此图可直观地对比各模型的精度差异.总体而言，LightGBM模型精度显著优于其余3种模型，可随着真实峰值振速增大，其精度逐渐降低.对于高峰值振速样本，LightGBM模型同样优于其余模型.

由于模型变量中包含了坐标轴独热变量，且对三轴数据同时建模，为更细观地验证所建立模型的精度和拟合优度，将测试集预测结果按照坐标轴区分并分别计算RMSE和R ².测试集三轴RMSE与R ²对比见图5.LightGBM在3个坐标方向上同样优于SVM、NN和萨氏公式.LightGBM在径向X上，RMSE和R ²分别为0.63与0.88；在切向Y上分别为1.12与0.83；在垂向Z上分别为0.69与0.94.LightGBM与NN在径向X上表现相近，但在切向Y上LightGBM的RMSE降低了23.3%，R ²提高了15.3%；SVM则在三方向上均劣于LightGBM.与SVM和NN这2种机器学习模型相比，LightGBM在切向Y上的精度和拟合优度有显著提升，使得其总体精度和拟合优度优于其余2种模型.相对于萨氏公式而言，LightGBM的RMSE在3个坐标方向上分别降低了63.4%、39.5%和68.3%，R ²则分别提高了18.9%、27.7%和42.4%.

在实际爆破工程设计时，爆破设计技术人员依据爆破安全规程GB 6722-2014中各保护物的峰值振速安全阈值对装药量、孔排距等爆破参数进行调整，因此，评判模型对于高峰值振速的拟合效果有重要的工程意义.本文选取0.5 cm/s作为爆破质点峰值振速安全阈值，将真实峰值振速高于0.5 cm/s的样本定义为高峰值振速样本，图6为测试集上高峰值振速样本RMSE与R ²对比.与图3相比，各模型的RMSE均增大，R ²均下降，其中萨氏公式R ²下降幅度达52.2%.此结果表明萨氏公式对于高峰值振速样本的拟合效果较差，难以准确预测危险振速.LightGBM模型R ²虽下降了7.9%，可仍达到0.82并显著优于其余模型，故LightGBM模型可用于预测危险振速，从而为爆破参数调整提供可靠依据.

4.4　LightGBM模型解释

基于SHAP模型可解释性框架对建立的LightGBM模型计算训练集中每个样本的Shapley值，将结果按式（28）计算可得模型变量重要度，结果如图7所示.在28个变量中，监测点距离的重要性远超过了其余变量，而总药量的重要性处于第二位，此结果与经验公式所选取的变量相符.究其原因，从变量数值变化范围而言，药量、排距等爆破参数的变化范围较小，而监测点距离变化达上千米，在本数据集中其他变量变化不大的情况下，监测点距离决定了所测质点峰值振速的数量级.同时，本研究的数据采集形式是在单次爆破下布置多监测点，对于同一次爆破其爆破参数不变，而监测点距离变化，导致算法对监测点距离赋予了更大的权重.此外，孔径、最大孔距、最大台阶高度和孔数等爆破参数均与单孔装药量及总药量直接相关，LightGBM算法带有惩罚项，会削弱单变量的重要性，这在一定程度上分散了单孔装药量与总药量的重要度.

图8为模型变量密度散点图，该图中纵坐标为变量，横坐标为Shapley值，圆点代表样本，圆点的颜色由蓝到红表示变量本身数值由小变大.Shapley值为正，则代表此变量值对模型预测数值有正贡献，反之则为负贡献.由此图可直观看出，监测点距离越大，Shapley值越小，质点峰值振速越低；总药量、最小排距、平均装药高度与最大孔距则与质点峰值振速呈正相关关系.

图9为模型单样本变量影响图，选取训练集中第一个样本为例.图中纵坐标为变量，其对应的灰色数字为该变量标准化值，横坐标为Shapley值， E［f（X）］为训练集所有样本f（X）的期望，为-0.273.单样本变量影响图直接展示了各变量对输出预测值的影响，输出预测值f（x）为期望值E［f（X）］ 加上各变量Shapley值，在本例中为-0.503，即预测的质点峰值振速为e^-0.503=0.60 cm/s.

5 结论

本研究通过记录每孔爆破参数和测取质点峰值振速的现场试验，结合LightGBM算法与SHAP模型可解释性框架，建立了适用于露天矿数码电子雷管逐孔起爆条件下的三轴质点峰值振速预测模型，主要结论如下：

（1）与USBM经验公式相比，取单次爆破总药量代入萨道夫斯基公式可得相对较高的精度与拟合优度.径向X、切向Y和垂向Z上RMSE分别为1.72，1.85与2.18；R ²分别为0.74，0.42与0.66.经验公式在切向Y上R ²较低，在垂向Z上RMSE较大.

（2）LightGBM在测试集上总体RMSE为0.86，总体R ²为0.89，均显著优于支持向量机、神经网络和萨道夫斯基公式：总体RMSE分别降低了25.9%、28.9%和55.2%，总体R ²分别提高了12.7%、9.9%和28.9%.

（3）LightGBM在径向X、切向Y和垂向Z上RMSE分别为0.63，1.12与0.69；R ²分别为0.88，0.83与0.94，在3个坐标方向上均优于支持向量机、神经网络和萨道夫斯基公式.LightGBM相较于萨道夫斯基公式而言，RMSE在3个坐标方向上分别降低了63.4%、39.5%和68.3%，R ²分别提高了18.9%、27.7%和42.4%.

（4）基于SHAP模型可解释性框架，除去方向轴变量外，监测点距离、总药量、最小排距、平均装药长度、孔径与最大孔距为对质点峰值振速影响程度最大的6个变量，其中监测点距离的影响程度远高于其余变量.监测点距离与质点峰值振速为负相关关系，总药量、最小排距、平均装药高度与最大孔距则与质点峰值振速呈正相关关系.

（5）对高于质点峰值振速安全阈值的样本，LightGBM拟合优度可达到0.82并显著优于支持向量机、神经网络和萨道夫斯基公式.所建立的LightGBM模型可准确预测露天矿数码电子雷管逐孔起爆条件下质点峰值振速，可直接应用于其他类似爆破工程，为爆破工程技术人员调整爆破参数提供可靠依据.

参考文献

[1] Aladejare, A. E., Lawal, A. I., Onifade, M., 2022. Predicting the Peak Particle Velocity from Rock Blasting Operations Using Bayesian Approach. Acta Geophysica, 70(2): 581-591. https://doi.org/10.1007/s11600-022-00727-5

[2] Hao, H.Z., Gu, Q., Hu, X.M., 2021. Research Advances and Prospective in Mineral Intelligent Identification Based on Machine Learning. Earth Science, 46(9): 3091-3106 (in Chinese with English abstract).

[3] He, J., Chalise, P., 2020. Nested and Repeated Cross Validation for Classification Model with High-Dimensional Data. Revista Colombiana de Estadística, 43(1): 103-125.

[4] He, L., Yang, R. S., Zhong, D. W., et al., 2021. Calculation of Equivalent Charge Weight Per Delay and Vibration Velocity Prediction for Millisecond Delay Blasting. Explosion and Shock Waves, 41(9): 129-141 (in Chinese with English abstract).

[5] Hosseini, S. A., Tavana, A., Abdolahi, S. M., et al., 2019. Prediction of Blast-Induced Ground Vibrations in Quarry Sites: A Comparison of GP, RSM and MARS. Soil Dynamics and Earthquake Engineering, 119:118-129.

[6] Jiao, H., Du, X., Zhao, M., et al., 2021. Nonlinear Seismic Response of Rock Tunnels Crossing Inactive Fault under Obliquely Incident Seismic P Waves. Journal of Earth Science, 32(5):1174-1189.

[7] Ke, G., Meng, Q., Finley, T., et al., 2017. LightGBM: A Highly Efficient Gradient Boosting Decision Tree. Proceedings of the 31st International Conference on Neural Information Processing Systems. Long Beach, California.

[8] Li, W. B., Fan, X. M., Huang, F. M., et al., 2021. Uncertainties of Landslide Susceptibility Modeling under Different Environmental Factor Connections and Prediction Models. Earth Science, 46(10): 3777-3795 (in Chinese with English abstract).

[9] Lundberg, S. M., Erion, G. G., Lee, S.I., 2018. Consistent Individualized Feature Attribution for Tree Ensembles. Methods, 5(13):25.

[10] Lundberg, S. M., Lee, S. I., 2017. A Unified Approach to Interpreting Model Predictions. Proceedings of the 31st International Conference on Neural Information Processing Systems. Long Beach, California.

[11] Tian, X., Song, Z., Wang, J., 2019. Study on the Propagation Law of Tunnel Blasting Vibration in Stratum and Blasting Vibration Reduction Technology. Soil Dynamics and Earthquake Engineering, 126: 105813. https://doi.org/10.1016/j.soildyn.2019.105813

[12] Trivedi, R., Singh, T., Raina, A., 2014. Prediction of Blast-Induced Flyrock in Indian Limestone Mines Using Neural Networks. Journal of Rock Mechanics and Geotechnical Engineering, 6(5): 447-454.

[13] Yang, X., Zhang, Y.P., Li, Y., et al., 2017. Analysis of Seismic Wave Propagation Law of Multi-Step Terrain Blasting and Correction of Blasting Vibration Velocity Prediction Formula. Mining Research and Development, 37(5): 78-83 (in Chinese with English abstract).

[14] Zhang, G. C., Guo, L. J., Jia, J. J., et al., 2020. Response Characteristics of Blasting Vibration to Rock Slope in an Open-Pit Iron Mine. China Mining Magazine, 29(12): 165-169 (in Chinese with English abstract).

[15] Zhang, W., Li, H., Li, Y. et al, 2021a. Application of Deep Learning Algorithms in Geotechnical Engineering: A Short Critical Review. Artificial Intelligence Review, 54: 5633-5673. https://doi.org/10.1007/s10462-021-09967-1

[16] Zhang, X., Peng, X., Li, X., et al., 2021b. Three- Dimensional Seismic Response in Complex Site Conditions: A New Approach Based on an Auxiliary-Model Method. Journal of Earth Science,32 (5): 1152-1165.

[17] Zhang, W., Phoon, K. K., 2022. Editorial for Advances and Applications of Deep Learning and Soft Computing in Geotechnical Underground Engineering. Journal of Rock Mechanics and Geotechnical Engineering, 14(3): 671-673.

[18] Zhang, Y.P., Yang, X., Zhu, X.X., 2018. Critical Velocity and Safety Distance Calculation under Influence of Blasting Vibration. Industrial Minerals & Processing, 47(8): 28-30 (in Chinese with English abstract).

[19] Zhou, J. R., Lu, W. B., Yan, P., et al., 2016. Frequency-Dependent Attenuation of Blasting Vibration Waves. Rock Mechanics and Rock Engineering, 49(10): 4061-4072. https://doi.org/10.1007/s00603-016-1046-5

[20] Zhu, M., Hu, X.T., Feng, Z.W., et al., 2021. Optimization of Blasting Vibration Velocity Attenuation Model under the Condition of Soil-Rock Strata. Science Technology and Engineering, 21(10): 4211-4218 (in Chinese with English abstract).

[21] 郝慧珍, 顾庆, 胡修棉, 2021. 基于机器学习的矿物智能识别方法研究进展与展望. 地球科学, 46(9): 3091-3106.

[22] 何理, 杨仁树, 钟东望, 等, 2021. 毫秒延时爆破等效单响药量计算及振速预测. 爆炸与冲击, 41(9): 129-141.

[23] 李文彬, 范宣梅, 黄发明, 等, 2021. 不同环境因子联接和预测模型的滑坡易发性建模不确定性. 地球科学, 46(10): 3777-3795.

[24] 杨曦, 张云鹏, 李岩, 等, 2017. 多台阶地形爆破地震波传播规律分析及爆破振速预测公式修正. 矿业研究与开发, 37(5): 78-83.

[25] 张耿城, 郭连军, 贾建军, 等, 2020. 某露天铁矿爆破振动对边坡的动态响应特征研究. 中国矿业, 29(12): 165-169.

[26] 张云鹏, 杨曦, 朱晓玺, 2018. 爆破振动影响下露天边坡临界振速及安全距离计算. 化工矿物与加工, 47(8): 28-30.

[27] 朱明, 胡鑫涛, 冯志威, 等, 2021. 土岩复合地层条件下爆破振速衰减模型优选. 科学技术与工程, 21(10): 4211-4218.

基金资助

中央高校基本科研业务费专项资金项目(2017QL05)