Abstract

Reliable prediction of blast-induced rock fragmentation is essential for optimizing blasting performance, reducing operational costs, and improving the efficiency of downstream mining operations. This study develops and evaluates hybrid machine learning models by integrating Backpropagation Neural Network (BPNN) and Support Vector Regression (SVR) with three optimization algorithms, namely the Zebra Optimization Algorithm (ZOA), Grey Wolf Optimizer (GWO), and Bayesian Optimization (BO). The resulting hybrid models, ZOA-BPNN, GWO-BPNN, BO-BPNN, ZOA-SVR, GWO-SVR, and BO-SVR, were compared with their standalone counterparts using the Mean Absolute Error (MAE), Mean Squared Error (MSE), Correlation Coefficient (R), and Coefficient of Determination (R²). Among the developed models, ZOA-BPNN achieved the highest predictive performance, with a validation MAE of 0.021, MSE of 0.046, R of 0.992, and R² of 0.984. Model interpretation using SHapley Additive exPlanations (SHAP) identified powder factor as the most influential variable governing fragment size, followed by charge per hole, charge length, and bench height, thereby providing valuable insights for blast design optimization. Overall, the results demonstrate that integrating metaheuristic optimization with machine learning substantially enhances prediction accuracy and robustness, with ZOA-BPNN emerging as the most effective framework for intelligent prediction and decision support in surface mine blasting.

Keywords

Rock Fragmentation, Metaheuristic Optimization, Zebra Optimization Algorithm, Back Propagation Neural Network, Support Vector Regression

Share and Cite:

1. Introduction

Blasting remains central to mining, tunneling, and civil construction, since it is the primary means of breaking rock masses for excavation, loading, and processing [1]-[3]. Despite advances in explosives and detonation technology, only about 20% - 30% of the released energy actually fragments and displaces rock [4] [5]; the rest escapes as seismic waves and air shock, producing flyrock, ground vibration, air overpressure, and backbreak [6]-[10]. These side effects raise safety and environmental concerns and add cost through operational delays and equipment damage.

Fragment size governs almost everything downstream. Well-sized fragments load and haul faster, protect crushers, and cut grinding energy, so blast quality is a direct lever on productivity and cost [11] [12]. Oversized boulders demand secondary breaking, while excess fines waste processing capacity and wear out equipment faster. Accurately predicting and controlling fragment size distribution is therefore central to running a cost-effective operation.

Fragmentation prediction has moved through several generations of models built on field data, expert judgment, and statistics [13]-[17]. Cunningham [18] [19] introduced the Kuz-Ram model, which combines the Rosin-Rammler distribution with empirical relationships to estimate mean fragment size (X₅₀) and uniformity (n) from explosive and rock-mass parameters. It remains widely used for its simplicity but has long been criticized for inaccuracy. Lawal [20] refined it using image analysis, striking a better balance between accuracy and practicality. The Crush Zone Model and Two-Component Model, developed at the Julius Kruttschnitt Mineral Research Centre [21] [22], added mechanisms for the fines generated near the borehole, and the Swebrec function [23] [24] later gave a better fit across the full-size range, feeding into the Kuznetsov-Cunningham-Ouchterlony (KCO) model [25]. Fragmentation-energy fan and percentile-based models extended this further [26] [27].

Even so, empirical models rest on fixed assumptions and site-specific calibration, so they generalize poorly to new conditions [28]. Traditional models such as the Kuznetsov equation [29] and Kuz-Ram [18] [19] were fitted to limited datasets and assume largely linear relationships, which cannot capture the non-linear interactions that actually govern fragmentation [30]-[32], and they were calibrated over narrow ranges of blast parameters [33]. The Swebrec function corrected some of Rosin-Rammler’s known weaknesses in underestimating fines and coarse fragments [23] [26] and gave a better S-shaped fit across the size spectrum [28], but it still depends on predetermined parameters that cannot flex with changing blast and geological conditions. This limitation, common to every empirical approach, is what motivates the machine learning approach used here.

Machine learning has become an attractive alternative for the same reason empirical models fall short: machine learning models can model complex, non-linear systems with many interacting variables, and larger datasets have made this practical [34]-[37]. In mining, AI already supports air overpressure and ground vibration prediction and production scheduling [38]-[40], as well as equipment condition monitoring [41]-[43]. Among these techniques, Back Propagation Neural Network (BPNN) and Support Vector Regression (SVR) are widely used for fragmentation prediction because they handle non-linear, multi-dimensional data well [44] [45] and have proven flexible and robust across mining and geotechnical applications [46]-[49]. Each still has a weakness, however: BPNN training is sensitive to weight initialization and network architecture, which risks overfitting, slow convergence, or entrapment in local minima [50], while SVR generalizes well on small or high-dimensional data but depends heavily on kernel and regularization choices that are usually tuned by trial and error [51].

This study addresses both weaknesses by pairing BPNN and SVR with metaheuristic optimizers, in two stages. First, the best baseline configuration is identified for each learner: BPNN is tested with Levenberg-Marquardt, Scaled Conjugate Gradient, and Bayesian Regularization training algorithms, while SVR is tested with linear, Gaussian, and polynomial kernels. Second, the Zebra Optimization Algorithm (ZOA), Grey Wolf Optimization (GWO), and Bayesian Optimization (BO) are applied to search for the best internal parameters of each baseline model, removing the guesswork of manual tuning and improving convergence and generalization.

The three optimizers were chosen because they search in different ways. ZOA mimics zebra herd behavior, specifically how zebras move collectively and evade predators, balancing broad exploration with local refinement [52]. GWO mimics the hunting hierarchy of grey wolves, using a structured pursue-and-encircle strategy well suited to fine-tuning a model’s core parameters [53]. BO takes a different, model-based approach by constructing a probabilistic surrogate of the objective function and using it to identify promising parameter combinations with relatively few computationally expensive evaluations [54].

This study develops and benchmarks, for the first time, six hybrid models, namely ZOA-BPNN, GWO-BPNN, BO-BPNN, ZOA-SVR, GWO-SVR, and BO-SVR, against their unoptimized counterparts, all implemented in MATLAB. To enhance model transparency and support engineering decision-making, the study further employs SHapley Additive exPlanations (SHAP), a model-agnostic interpretability technique derived from cooperative game theory that quantifies the contribution of each input variable to model predictions. SHAP provides both global feature importance and insight into how individual blasting parameters influence predicted fragmentation outcomes, thereby improving the interpretability and practical applicability of the developed models. Finally, the comparative evaluation provides practical insights into optimizer-model compatibility, demonstrating that the effectiveness of a metaheuristic algorithm depends not only on its optimization capability but also on its interaction with the underlying machine learning model. These contributions provide a more reliable, transparent, and interpretable framework for intelligent blast-induced rock fragmentation prediction.

The rest of the paper is organized as follows. Section 2 reviews related hybrid modelling studies. Section 3 describes the study site, and Section 4 the dataset. Sections 5 and 6 present the modelling methodology and model development process. Section 7 reports and discusses the results, including SHAP-based interpretability, and Section 8 concludes with implications and directions for future work.

2. Related Studies on Hybrid Machine Learning Fragmentation Models

This section reviews hybrid machine learning approaches used to predict blast-induced rock fragmentation, to position the present study within current research. A systematic search of ScienceDirect and Google Scholar used terms such as “rock fragmentation,” “metaheuristic optimization,” “hybrid machine learning,” and “blasting prediction.” Studies were retained only if they combined a machine learning model with a metaheuristic optimizer for fragmentation prediction, and only if published between 2015 and 2025.

Table 1 summarizes the resulting studies, their hybrid techniques, base models, optimizers, and reported performance. Notably, few of these studies first benchmark different learning functions within the standalone model before applying metaheuristic optimization—a gap this study addresses directly.

Table 1. A review of hybrid models for predicting blast-induced rock fragmentation.

JSO: Jellyfish Search Optimizer; LSSVM: Least Squares Support Vector Machine; AFSA: Artificial Fish Swarm Algorithm; APSO: Adaptive Particle Swarm Optimization; RF: Random Forest; ERT: Extremely Randomized Trees; AHA: Artificial Hummingbird Algorithm; GBoost: Gradient Boosting; AdaBoost: Adaptive Boosting; ANN: Artificial Neural Network; ABC: Artificial Bee Colony; ANFIS: Adaptive Neuro-Fuzzy Inference System; PSO: Particle Swarm Optimization; GA: Genetic Algorithm; FFA: Firefly Algorithm; BGAM: Boosted Generalized Additive Model; CSO: Cat Swarm Optimization; GPR: Gaussian Process Regression; LightGBM: Light Gradient Boosting Machine; BFO: Bacterial Foraging Optimization; BOA: Bayesian Optimization Algorithm; ICA: Imperialist Competitive Algorithm.

3. Site Description

The mine lies in one of Ghana’s most productive gold belts, about 10 km southwest of Tarkwa, 70 km from Sekondi, and 320 km from Accra. Geologically, it sits within the Tarkwaian Group of the West African Craton, made up mainly of metavolcanic and metasedimentary rocks from the Birimian Supergroup [65] (Figure 1). Operations follow a standard open-pit sequence of drilling, blasting, loading, and haulage, with blasting typically conducted daily unless delayed by weather or other operational factors. The mine operates several production drill rigs, including Rock Commander, Epiroc SmartROC T45, Epiroc DML, and Sandvik drill rigs. However, the dataset analyzed in this study was obtained from production blasts drilled primarily using Sandvik drill rigs with hole diameters of 115 mm and 127 mm, depending on the ground conditions. Accordingly, the descriptive statistics presented in Table 2 reflect only the hole diameters represented in the analyzed dataset.

Figure 1. Study area.

4. Data Description

A total of 498 blasts were monitored at the mine, with data recorded across three parameter groups: blast design, explosive characteristics, and rock properties. Blast design parameters—spacing (S), burden (B), initial and final stemming height (IS, FS), bench height (H), and charge length (L), all in meters, plus hole diameter (D) in millimeters—were measured alongside explosive parameters: charge per hole (C, kg) and powder factor (PF, kg/m³). Powder factor was calculated as charge per hole divided by blast volume, where blast volume is the product of burden, spacing, and bench height (PF = C/(B × S × H), kg/m³). Rock properties were represented by uniaxial compressive strength (UCS, MPa). These ten parameters served as the independent variables, with the actual mean fragment size (X₅₀, cm) as the dependent variable. Mean fragment size (X₅₀) for each blast was determined using WipFrag image-analysis software applied to digital images of the resulting muckpile, captured immediately after each blast across the full pile to ensure representative coverage. Oversized fragments that could not be reliably delineated within the image frame were handled using WipFrag’s boulder-tracing function to correct the size distribution. Table 2 summarizes the dataset statistically.

Table 2. Statistical overview of the rock fragmentation dataset.

Figure 2 shows the Pearson correlation matrix, used to quantify linear relationships among the input variables and X50. Bench height correlates strongly with charge length (R ≈ 0.98) and charge per hole (R ≈ 0.96), and initial and final stemming height are perfectly correlated (R = 1.00); hole diameter also correlates with both stemming heights (R ≈ 0.71), and charge length with charge per hole (R ≈ 0.94). Against X₅₀, bench height, stemming heights, charge length, and charge per hole all show moderate positive correlation—coarser blast designs tend to produce coarser fragments, while powder factor shows a strong negative correlation (R ≈ −0.85), confirming its dominant role in controlling fragment size. This matches expected blasting mechanics: more explosive energy per unit volume breaks rock finer.

Figure 2. Illustration of Pearson correlation plot.

Figure 3 shows the Spearman rank correlation, which captures monotonic relationships without assuming linearity—a useful check given the non-linear interactions typical of blasting data. The Spearman results largely confirm the Pearson findings, with initial and final stemming height again perfectly correlated (ρ = 1.00) and both strongly linked to charge per hole (ρ ≈ 0.85); bench height also correlates strongly with IS, FS, and C (ρ ≈ 0.76). The inverse relationship between powder factor and X50 is even stronger here (ρ ≈ −0.95), confirming that this effect holds consistently across both linear and non-linear analysis.

Figure 3. Illustration of Spearman correlation plot.

5. Methodology

5.1. Backpropagation Neural Network

Backpropagation Neural Networks (BPNN) are a simple, effective type of feedforward network [66] [67], made up of an input layer, one or more hidden layers, and an output layer. Each neuron computes a weighted sum of its inputs, adds a bias, and passes the result through an activation function, as shown in Equation (1) [68].

$y_j=f\left({\displaystyle \sum _{i=1}^n\left(w_{ij}\times x_i\right)+b_j}\right)$ (1)

where the values at the nodes in the previous layer (i) and the layer currently active (j) are illustrated by $x_i$ and$y_j$ , respectively. The total nodal values obtained from the previous layer are depicted by n. The network’s biases and weights are denoted by $b_j$ and $w_{ij}$ .

For this study, a tangential hyperbolic function was utilized based on its ability to produce more precise outcomes. The function’s range is bounded by −1 and 1, and its mathematical expression is given in Equation (2) as:

$y_j=f\left(net\right)=2\times \left(\frac{1}{e^{-2.net}+1}-0.5\right)$ (2)

The network’s final output ($y_j$ ) is calculated by feeding the computed internal network value (net) through an activation function (f). During the training phase, weights and biases values are repeatedly adjusted to minimize the error between the predicted network values and the actual target values This mechanism is essentially a minimization problem [69]. For this study, the BPNN model was built and executed using MATLAB 2023a software. For a deeper look at how the BPNN works, comprehensive details on the BPNN itself are available in resources like the work presented by Jain et al. [70].

5.2. Support Vector Regression

Support Vector Regression (SVR) is a variant of the Support Vector Machine adapted for regression tasks. As Zhang and O’Donnell [71] note, its structure is not fixed in advance—the input vectors that define it are determined during training. For a training set with input vector x and target f, the SVR estimator f(x) is given in Equation (3), where w is the weight vector, b the bias, and K(x) a non-linear basis function mapping the input into a higher-dimensional space.

$f\left(x\right)=b+w\cdot \left(\kappa \left(x\right)\right)$ (3)

where weight vector $w=\left(w_1,w_2,w_3,\cdots ,w_z\right)$ , the bias term is b, and the basis function vector $\kappa \left(x\right)=\left({\kappa }_1\left(x_1\right),\cdots ,\kappa \left(x_z\right)\right)$ referring to a series of non-linear changes. The non-linear basis function is used to convert the R-dimensional input vector from a low-dimensional phase to a high-dimensional space $\kappa \left(x\right)$ . Support Vector Regression (SVR) aims to identify the optimal separating hyperplane that enhances the model’s capacity to learn from data and effectively address a given problem [72]. Solving Equation (3) necessitates the use of Equation (4), while adhering to the constraints outlined in Equation (5). It is important to highlight that the structure of an SVR is not fixed in advance; rather, the input variables that contribute to the model are typically selected throughout the training process.

$\underset{w,b,\xi ,{\xi }^{\ast }}{\text{minimize}}\frac{1}{2}{〚w〛}^2+M\cdot {\displaystyle \sum _{J=1}^Z\left({\xi }_j\text{+ }{\xi }_j{}^{\ast }\right)}$ (4)

$\text{s}\text{.t}\text{.}\{\begin{array}{c}\left(b+w^{\text{T}}\kappa \left(x_j\right)\right)-y_j\le \epsilon -{\xi }_j\\ b+w^{\text{T}}\kappa \left(x_j\right)-y_j\le {\xi }_j{}^{\ast }+\epsilon \\ {\xi }_j{}^{\ast },{\xi }_j\ge 0\end{array}\}j=1,2,\cdots ,Z$ (5)

where the positive trade-off parameters, ɛ and M, regulate the measured deviation scale in maximization, while the slack factors $\xi$ and ${\xi }^{\ast }$ resolve the errors using the expression in Equation (5).

Support Vector Regression (SVR) relies on kernel functions, which function similarly to activation functions by implicitly projecting input data into a higher-dimensional space to enable better separation of classes. To construct accurate and efficient models, commonly used SVR kernels were evaluated, including the linear kernel, Gaussian kernel, and polynomial kernel functions [73]. These are represented by the following mathematical expressions, respectively.

$K\left(x,\overline{x}\right)=x^T\cdot \overline{x}$ (6)

$G\left(x,\overline{x}\right)=\exp \left(\gamma {\Vert -x+\overline{x}\Vert }^2\right)$ (7)

$P\left(x,\overline{x}\right)={\left(c+\left(x^T\cdot \overline{x}\right)\right)}^q$ (8)

where x, $\overline{x}$ represents the input feature vectors, $x^T\cdot \overline{x}$ is the product of the vectors γ is the scaling factor, c is the coefficient term, and $-x+\overline{x}$ refers to the Euclidean distance between two vectors.

5.3. Zebra Optimization Algorithm

Metaheuristic algorithms draw on natural behaviors to solve complex optimization problems [52]. The Zebra Optimization Algorithm (ZOA) used here is a swarm-based method known for fast convergence, modelled on how zebra herds forage and evade predators together, balancing broad exploration with local refinement.

5.4. Grey Wolf Optimization

The Grey Wolf Optimizer (GWO) mimics the social hierarchy and hunting behavior of grey wolves [53] and is widely used across engineering and machine learning [74]. It identifies the three best solutions (α, β, δ) as leaders that guide the remaining wolves toward the global optimum, through three phases: encircling, hunting, and attacking the prey.

5.5. Bayesian Optimization

Bayesian Optimization (BO) uses a probabilistic surrogate model, typically a Gaussian Process, to approximate an expensive or complex objective function [75], rather than evaluating it directly many times. The surrogate gives a mean estimate and an uncertainty at each point, and BO uses this to choose where to sample next by maximizing an acquisition function, balancing exploration of new regions against exploitation of known good ones [75] [76].

6. Model Development

This section outlines the data preprocessing, formulation of the baseline models, the metaheuristic optimization procedures used to develop the hybrid predictive models, and the SHAP-based interpretation approach adopted to explain the influence of the input variables on model predictions.

6.1. Data Preprocessing

The complete dataset was randomly divided into a training dataset (80%) and a validation dataset (20%). The training dataset was used for model development and hyperparameter optimization. During model development, five-fold cross-validation was performed exclusively within the training dataset to identify the optimal model configurations and hyperparameters. In each fold, normalization parameters were computed (using Equation 9) only from the corresponding training subset and subsequently applied to the associated validation fold without refitting, thereby preventing information leakage. After model development, the final models were evaluated using the independent validation dataset, which was not involved in the training or hyperparameter optimization process [77].

$N_i=\left\{\frac{2}{(s_{\max }-s_{\min })}\times \left(s_i-s_{\min }\right)\right\}-1$ (9)

where, the normalized data point of variable $s_i$ is represented by $N_i$ , the maximum scores are defined by $s_{\max }$ , and the minimum scores of the measured samples are indicated by $s_{\min }$ .

6.2. Baseline Model Formulation

The baseline BPNN and SVR models were developed by evaluating multiple candidate learning strategies and selecting the best-performing configurations. For the BPNN, three training algorithms, namely Levenberg-Marquardt (LM), Scaled Conjugate Gradient (SCG), and Bayesian Regularization (BR), were investigated. The number of hidden neurons was varied from 1 to 150 to determine the optimal network architecture. Model performance was evaluated using the Mean Absolute Error (MAE), Mean Square Error (MSE), correlation coefficient (R), and coefficient of determination (R²). The training algorithm and network architecture that produced the lowest prediction errors and highest goodness-of-fit metrics on the validation dataset were selected as the baseline BPNN model.

For the SVR, three kernel functions, namely Linear, Polynomial, and Gaussian, were evaluated. The kernel yielding the best predictive performance was selected as the baseline model. The principal SVR hyperparameters, including the regularization parameter (C), kernel parameter (σ), and insensitive loss parameter (ε), were subsequently considered for optimization to further improve predictive performance.

6.3. Hybrid Model Development and Hyperparameter Optimization

The selected baseline BPNN and SVR models were subsequently integrated with the Zebra Optimization Algorithm (ZOA), Grey Wolf Optimizer (GWO), and Bayesian Optimization (BO) to improve predictive accuracy through optimal hyperparameter selection. For the BPNN, the optimization algorithms searched for the optimal network weights and biases, whereas for the SVR they optimized the regularization parameter (C), kernel parameter (σ), and insensitive loss parameter (ε). The objective function was defined as the mean squared error (MSE) computed over the training folds.

The optimization variables were normalized within the interval [−1, 1] to maintain a compact search space, improve computational efficiency, and achieve a balanced trade-off between exploration and exploitation. The preliminary experiments were conducted to identify ranges that consistently produced stable convergence while avoiding unnecessary computational cost. For ZOA, the population size varied from 10 to 120 zebras with iteration numbers ranging from 100 to 1500. The swapping probability (P_s) and exploration constant (C) were fixed at 0.5. For GWO, the population size ranged from 1 to 50 wolves, while the maximum number of iterations varied between 10 and 150 for the BPNN and between 10 and 200 for the SVR. Bayesian Optimization employed a Gaussian Process surrogate model with a Matérn kernel, using 5 to 50 initial sampling points and 50 to 400 optimization iterations. Expected Improvement and Upper Confidence Bound acquisition functions guided the search for the optimal solution.

Each optimization algorithm iteratively updated candidate solutions until the maximum number of iterations or the convergence criterion was satisfied. The optimal parameter set identified by each optimization algorithm was then used to retrain the corresponding BPNN or SVR model. The optimized models were subsequently evaluated using the validation dataset, while the training dataset was used for model development and optimization. Figure 4 presents the overall model development and optimization framework. All simulations were performed on an HP laptop equipped with an AMD Ryzen 5 processor and 12 GB of RAM.

6.4. SHAP-Based Model Formulation

To enhance the interpretability of the developed machine learning framework, SHapley Additive exPlanations (SHAP) [78] were employed to quantify the contribution of each input variable to model predictions. SHAP is a model-agnostic interpretability technique derived from cooperative game theory that assigns each feature a Shapley value, representing its average marginal contribution to the prediction across all possible combinations of input variables. Unlike conventional feature importance measures, SHAP provides both the magnitude and direction of each feature’s influence, thereby offering consistent and locally accurate explanations of model behavior.

Following model development and performance evaluation, SHAP analysis was performed on the model that demonstrated the highest predictive performance. Global feature importance was quantified using the mean absolute SHAP values computed from all observations in the validation dataset, allowing the relative influence of each blasting parameter on fragment size prediction to be ranked. In addition, SHAP summary plots were generated to illustrate how variations in individual input variables influenced the predicted fragmentation outcomes.

Figure 4. Conceptual framework of BPNN and SVR models integrated with optimization algorithms for blast-induced rock fragmentation prediction.

6.5. Model Performance Evaluation

Model performance was assessed using MAE, MSE, correlation coefficient (R), and coefficient of determination (R²), which together capture prediction accuracy, error magnitude, and fit. Lower MAE and MSE indicate better accuracy, while higher R and R² indicate stronger agreement between predicted and observed values. These metrics are defined in Equations. (10)-(13) [79]-[81], where OF and PF denote the observed and predicted fragment sizes, their subscripted averages denote the dataset means, and n is the total number of samples.

$\text{MAE}=\frac{1}{n}\times {\displaystyle \sum _{i=1}^n\left|O{F}_i-P{F}_i\right|}$ (10)

$\text{MSE}=\frac{1}{n}\times {\displaystyle \sum _{i=1}^n{\left(O{F}_i-P{F}_i\right)}^2}$ (11)

$\text{R}=\frac{{\displaystyle \sum _{i=1}^n\left(O{F}_i-O{F}_{av}\right)\left(P{F}_i-P{F}_{av}\right)}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(O{F}_i-P{F}_{av}\right)}^2}}\times \sqrt{{\displaystyle \sum _{i=1}^n{\left(P{F}_i-P{F}_{av}\right)}^2}}}$ (12)

${\text{R}}^2=1-\left(\frac{{\displaystyle \sum _{i=1}^n{\left(O{F}_i-P{F}_i\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(O{F}_i-O{F}_{av}\right)}^2}}\right)$ (13)

In the statistical equations, OF_i and PF_i denote the observed and predicted fragment sizes for the i-th sample, respectively. The terms OF_av and PF_av represent the mean values of the observed and predicted fragment sizes across the entire dataset. The total number of samples used for model development is designated by n.

7. Results and Discussion

This section presents the results of the modelling process: first the baseline training algorithms and kernel functions with their hyperparameter settings, then the comparative performance of the optimized hybrid models, followed by feature importance and the broader implications of the findings.

7.1. Developed Models

Table 3 shows that all three BPNN training algorithms exhibited good generalization, with training performance consistently outperforming validation performance by a small margin. This pattern indicates effective learning while maintaining good generalization to unseen data, with no evidence of overfitting. Among the evaluated algorithms, Bayesian Regularization achieved the best overall performance, recording the lowest training MAE (0.362) and MSE (0.781), together with the highest R (0.836) and R² (0.698). Its validation performance remained consistently strong (MAE = 0.389, MSE = 0.853, R = 0.812, and R² = 0.659), demonstrating superior predictive capability and robust generalization compared with the Levenberg-Marquardt and Scaled Conjugate Gradient algorithms. Consequently, Bayesian Regularization was selected as the baseline BPNN training algorithm for subsequent hybrid optimization.

Table 3. Statistical analysis results of the training algorithms modelled in the BPNN.

Table 4 shows a similar pattern for the SVR models, with all three kernel functions exhibiting slightly better performance during training than during validation. The small differences between the two phases indicate good generalization ability while suggesting that the models did not overfit the training data. Among the evaluated kernels, the Gaussian kernel demonstrated the best overall performance, achieving the lowest training MAE (2.035) and MSE (4.612), together with the highest R (0.615) and R² (0.378). Its validation performance remained comparably strong (MAE = 2.127, MSE = 4.831, R = 0.592, and R² = 0.350), confirming its superior predictive capability and robustness compared with the Linear and Polynomial kernels.

Table 4. Statistical performance of candidate SVR kernel functions.

Based on these results, Bayesian Regularization and the Gaussian kernel were selected as the baseline BPNN training algorithm and SVR kernel, respectively. These baseline models were subsequently optimized using the Zebra Optimization Algorithm (ZOA), Grey Wolf Optimizer (GWO), and Bayesian Optimization (BO). The optimization variables were bounded within the interval [−1, 1] to maintain a compact search space, improve computational efficiency, and achieve a balanced trade-off between exploration and exploitation during the optimization process. The reported values correspond to the optimal solution returned by each optimization algorithm. Table 5 presents the baseline model architectures together with the optimal hyperparameter settings for each optimization algorithm.

Table 5. Baseline model architecture and optimal hyperparameter configurations for the optimization algorithms integrated with the BPNN and SVR models.

7.2. Statistical Comparison of Optimized Models

Table 6 compares the training and validation performance of the baseline models and their optimized hybrid counterparts.

Table 6. Training and validation performance of the baseline and optimized models.

Among all the evaluated models, ZOA-BPNN demonstrated the best predictive performance, achieving a validation MAE of 0.021, MSE of 0.046, R of 0.992, and R² of 0.984. These results indicate that the model explained approximately 98.4% of the variability in rock fragmentation, highlighting its excellent predictive accuracy and strong generalization capability. The superior performance of ZOA-BPNN can be attributed to the complementary strengths of the Zebra Optimization Algorithm and the BPNN. Conventional BPNN training relies on gradient-based optimization, which is susceptible to premature convergence and local optima. In contrast, ZOA combines effective global exploration with local exploitation, enabling a more comprehensive search of the network’s weight and bias space. By optimizing the initial weights and biases before training, ZOA provided a more favorable starting point for the learning process, thereby reducing convergence to poor local solutions, improving model stability, and enhancing prediction accuracy. The validation coefficient of determination (R²) results presented in Figure 5 further support these findings, with ZOA-BPNN exhibiting the highest R² value among all the developed models. GWO-BPNN and BO-BPNN also outperformed the standalone BPNN, increasing the coefficient of determination from 0.659 to 0.776 and 0.849, respectively. These improvements demonstrate the effectiveness of metaheuristic optimization in enhancing the predictive capability of neural network models, although neither algorithm matched the optimization efficiency achieved by ZOA. Among the support vector regression models, the standalone SVR produced the lowest predictive performance, with a validation MAE of 2.127 and an R² of 0.350. Nevertheless, integrating SVR with metaheuristic algorithms substantially improved its performance. ZOA-SVR emerged as the strongest SVR-based model, achieving a validation MAE of 0.033 and an R² of 0.837. This result further demonstrates the robustness and adaptability of ZOA, indicating that its optimization capability extends beyond neural networks to kernel-based learning algorithms.

Figure 5. Comparison of coefficient of determination results of the applied models.

7.3. Model Ranking

To provide a comprehensive comparison of the developed models, the validation results were ranked using the four evaluation metrics, namely MAE, MSE, R, and R². For MAE and MSE, lower values indicate superior predictive performance and were therefore assigned higher rankings. Conversely, for R and R², higher values indicate stronger predictive capability and received higher rankings. Each model was ranked from 1 (best) to 8 (worst) for each metric, and the individual rankings were summed to obtain an overall score. Consequently, a lower total score represents better overall model performance. This ranking approach enables simultaneous consideration of prediction accuracy, goodness-of-fit, and model reliability, thereby facilitating an objective comparison of the standalone and hybrid intelligent models. Table 7 ranks all eight models across the four metrics.

Table 7. Overall ranking of the developed intelligent models based on validation performance.

The ranking confirms that ZOA-BPNN was the best-performing model, achieving first place across all four evaluation metrics with the lowest overall score of 4. ZOA-SVR ranked second, followed by BO-BPNN, demonstrating that the Zebra Optimization Algorithm consistently produced the greatest performance improvements irrespective of the underlying learning algorithm. GWO-BPNN also outperformed the standalone BPNN, while BO-SVR and GWO-SVR showed substantial improvements over the conventional SVR model. Among the baseline models, BPNN ranked fifth and considerably outperformed the standalone SVR, which placed last. These findings demonstrate that integrating metaheuristic optimization with machine learning models substantially enhances predictive performance, although the magnitude of improvement depends on the compatibility between the optimization algorithm and the underlying predictive model.

7.4. SHAP-Based Interpretation of Feature Influence on Fragmentation Prediction

Following the comparative evaluation of all developed models, the ZOA-BPNN model demonstrated the highest predictive performance and was therefore selected for SHAP-based interpretation. Figure 6 presents the SHAP summary plot, which ranks the blasting parameters according to their mean absolute SHAP values and illustrates both the magnitude and direction of each variable’s influence on the predicted fragment size.

Powder factor emerged as the most influential predictor by a substantial margin, exhibiting the widest range of SHAP values (approximately −15 to +25). Higher powder factor values consistently shifted the predictions towards smaller fragment sizes, confirming its dominant role in controlling blast fragmentation. This finding is consistent with established blasting theory, where increasing the explosive energy per unit rock volume generally enhances rock breakage and produces finer fragmentation.

The SHAP feature importance plot (Figure 7) further corroborates these findings, ranking charge per hole and charge length as the second and third most influential variables, with mean absolute SHAP values of approximately 0.325 and 0.230, respectively. Bench height followed with a mean absolute SHAP value of approximately 0.194, while the remaining variables contributed comparatively less to the model predictions. Collectively, these results indicate that explosive energy distribution and blast geometry are the principal factors governing fragmentation performance.

It is important to note that, in this study, the blast design variables exhibit strong correlations, as demonstrated by the Pearson and Spearman analyses. In particular, powder factor, charge per hole, charge length, bench height, and the stemming variables are operationally related and jointly define the distribution of explosive energy within a blast. Consequently, their SHAP values should be interpreted within the context of these interdependent blasting parameters rather than as completely independent effects. Despite these correlations, the consistent identification of powder factor as the dominant predictor reinforces its critical influence on fragment size prediction and supports the robustness and engineering reliability of the developed ZOA-BPNN model.

Figure 6. SHAP summary plot showing global feature influence.

Figure 7. Mean absolute SHAP value bar chart for feature importance ranking.

8. Conclusions and Future Proposals

This study developed and evaluated hybrid machine learning models for predicting blast-induced rock fragmentation by integrating BPNN and SVR with the Zebra Optimization Algorithm (ZOA), Grey Wolf Optimizer (GWO), and Bayesian Optimization (BO). The results demonstrate that metaheuristic optimization significantly improved the predictive performance of both machine learning models compared with their standalone counterparts. Among the developed models, ZOA-BPNN achieved the best performance, attaining a validation MAE of 0.021, MSE of 0.046, R of 0.992, and R² of 0.984, thereby explaining approximately 98.4% of the variability in fragment size. ZOA-SVR also performed well, while BO-BPNN and GWO-BPNN provided competitive alternatives.

SHAP analysis identified powder factor as the most influential predictor, followed by charge per hole, charge length, and bench height, confirming that explosive energy distribution and blast geometry are the primary factors controlling fragmentation. These findings enhance the interpretability of the developed models and provide practical guidance for blast design optimization.

Although the proposed framework demonstrated excellent predictive capability, the study was based on data from a single mine and did not explicitly consider additional rock mass characteristics and blast sequencing variables. Future research should validate the models using multi-mine datasets, incorporate additional geological and operational parameters, and extend the framework into intelligent real-time decision-support systems for surface mine blasting.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References