INTEGRATING ANN PREDICTION WITH HONEYBEE OPTIMISATION FOR FLYROCK MINIMISATION IN OPEN-PIT MINING

2. Field Investigation and Data Collection

To develop a predictive and optimization model for blast parameters aimed at minimizing flyrock, a total of 334 new datasets were collected from the Sungun copper mine, located in northwest Iran, approximately 35 km from Varzaghan in East Azerbaijan province. This mine is characterised by a hydrothermal ore deposit with mineralisation associated with the Cenozoic Sahand-Bazman orogenic belt, hosted within altered quartz-monzonite rocks. The primary ore minerals include chalcopyrite, pyrite, chalcocite, cuprite, and malachite, with copper being the principal extracted resource. Additionally, the deposit contains economically significant amounts of gold, silver, and molybdenite. Figure 1 illustrates the operational area of the Sungun copper mine.

Figure1. Drilling and Blasting Area in Sungun Copper Mine

Each dataset comprises key blasting parameters, including stemming, burden and spacing, blast hole diameter, blast hole length, sub-drilling length, specific charge, and flyrock distance. A descriptive statistical analysis was conducted on the collected datasets, with the results summarised in Table 1.

Table 1. Statistical Characteristics of the Data Used in the Network

The selected input variables include both fundamental blast design parameters and a derived parameter known as Specific Charge. Although Specific Charge is mathematically related to other inputs, its inclusion as a distinct variable is justified because it represents the overall concentration of explosive energy a critical factor in flyrock generation that geometric parameters alone do not fully capture. It is important to note that while uncontrollable geological factors can also contribute to flyrock, this study focuses on optimizing the controllable design parameters to minimize risk. The minimum recorded flyrock distance of 10 meters, although relatively short, was classified as a flyrock event according to site-specific safety protocols, as it landed beyond the designated safety perimeter, thereby necessitating its inclusion in the dataset.

3. Artificial Neural Network

Artificial neural networks (ANNs) are designed to mimic the structure and functionality of the human brain. Essentially, an ANN consists of interconnected neurons arranged in multiple layers, where each neuron acts as a simple processing unit that transmits information to others within the network. A large collection of these neurons forms a neural network. One of the most widely used learning algorithms in perceptron neural networks is the feedforward learning algorithm (Laguna & Martí, 2002), which operates based on the error correction learning law, a generalisation of the least-means algorithm.

Feedforward neural networks consist of three main layers: an input layer, a hidden (middle) layer, and an output layer. While the number of hidden layers is not restricted, a single hidden layer is often sufficient for solving complex nonlinear problems. The feedforward learning process is divided into two stages: the feedforward stage and the backward (error correction) stage. In the feedforward stage, inputs are passed sequentially through each layer, ultimately producing an output as the responses of the network. During the forward phase, synaptic weights are initialised. In the backward phase, these weights are adjusted based on error correction rules. The difference between the predicted response of the network and the desired (expected) response, known as the error signal, is propagated backward through the network, refining the synaptic weights to improve prediction accuracy.

To evaluate the network's performance, the coefficient of determination (R²) and root mean square error (RMSE) were used to measure the correlation between predicted and actual flyrocks. In this study, the neural network was implemented using the Lüneburg-Marquette learning algorithm, one of the most widely used optimisation techniques in ANN training. Based on prior research, the network structure was selected with three layers, an input layer with 7 input variables, a hidden layer and an output layer with a single output. A single hidden layer is sufficient to approximate any nonlinear function. Various studies have explored optimal methods for determining the number of neurons in the hidden layer, offering mathematical approaches to avoid the trial-and-error method. Equations 1 to 6 outline the established methodologies for selecting the appropriate number of hidden neurons (Hecht-Nielsen, 1987; Hornik et al., 1989; Masters, 1993; Paola, 1994; Ripley, 1993; Wang, 1994).

[TeX:] H \leq 2N_{i} + 1$H\le 2N_i+1$ (1)

[TeX:] H = \frac{N_{i} + N_{o}}{2}$H=\frac{N_i+N_o}{2}$ (2)

[TeX:] H = \frac{2 + N_{0} \times N_{i} + 0.5N_{0} \times (N_{0}^{2} + N_{i}) - 3}{N_{i} + N_{0}}$H=\frac{2+N_0\times N_i+0.5N_0\times (N_0^2+N_i)-3}{N_i+N_0}$ (3)

[TeX:] H = \frac{2N_{i}}{3}$H=\frac{2N_i}{3}$ (4)

[TeX:] H = \sqrt{N_{i} \times N_{0}}$H=\sqrt{N_i\times N_0}$ (5)

[TeX:] H = 2N_{i}$H=2N_i$ (6)

In this study, the neural network was implemented using the Levenberg-Marquardt learning algorithm. This algorithm was chosen due to its high efficiency and fast convergence rates for training small to medium-sized feedforward neural networks, making it well-suited for the dataset size in this research (Ampazis & Perantonis, 2000).

In this study, N_i represents the number of inputs, while N_o denotes the number of outputs in the model. Based on recommended parameter values, the Levenberg-Marquardt learning algorithm was implemented using a neural network with a hidden layer containing 2 to 15 neurons. To determine the most effective model, the coefficient of determination (R²) and root mean square error (RMSE) were evaluated for both the training and testing phases. A comparative analysis was conducted by assigning scores to each model configuration. The scoring process involved awarding the highest score to models with the largest R² values, while lower scores were assigned as R² decreased. Conversely, models with the lowest RMSE values received the highest scores, with scores decreasing as RMSE increased. At the end of the evaluation, the scores for each model configuration were summed to determine the overall performance ranking (Kaastra & Boyd, 1996).

The scoring process involved four criteria: R² for training, RMSE for training, R² for testing, and RMSE for testing. For each criterion, the 14 models were ranked from best (rank 14) to worst (rank 1). For R², a higher value received a higher rank, while for RMSE, a lower value received a higher rank. The final score for each model was the sum of its four ranks, allowing for a balanced assessment of both accuracy and generalisation ability.

Figure 2 presents a schematic diagram of how the neural network works.

Figure2. Schematic of a Basic Artificial Neural Network (ANN) with Input, Hidden, and Output Layers

4. Honeybee Optimisation Algorithm

The Honeybee Optimisation Algorithm is inspired by the collective behaviour of bee colonies, where individual bees, though simple on their own, work together to form a highly organised and efficient system for discovering and exploiting nectar resources. Within a colony, bees are divided into three main groups, each assigned a specific task in foraging.

The first group consists of scout bees, responsible for exploring the environment and new food resources. Once a scout bee finds a suitable resource, it returns to the hive and communicates its location through a movement known as the circular dance. The second group, the worker bees, focuses on exploiting these discovered food resources, ensuring efficient nectar collection. The third group, the onlooker bees, remains in the hive, observing the circular dance of the scout bees. Based on the quality of the reported resources, onlooker bees select the most promising sites for further extraction. The Honeybee Optimisation Algorithm was first introduced by Karaboga in 2005 (Karaboga, 2005). Since then, extensive research has been conducted on bee behaviour in nature, leading to the development of optimisation algorithms inspired by their social structure. While this algorithm has been applied in various fields, one of its notable implementations in mining is the prediction and optimisation of the backbreak caused by blasting operations (Ebrahimi et al., 2016; Kanellopoulos & Wilkinson, 1997; Rezaei et al., 2011; Sayadi et al., 2013; Zorlu et al., 2008). The algorithm follows a structured process consisting of four key stages. As illustrated in Figure 3, the workflow of the Bee Algorithm is structured into three core phases: Initialization, Iteration, and Termination Condition Evaluation, ensuring a systematic optimization process. The main steps of the Bee Algorithm are outlined in Karaboga & Basturk (2007):

Figure3. Flowchart of the Bee Algorithm illustrating its three primary phases: Initialisation, Iteration, and Termination Condition evaluation.

Step 1: In the initial phase of the algorithm, the bee population is equally divided into worker bees and non-worker bees. Each food source is assigned a single worker bee, which means the number of worker bees corresponds directly to the number of food sources around the hive. Consequently, within the defined solution space, the initial solution is generated based on the number of food sources. Once these initial solutions are established, their respective values are compared using problem-specific evaluation functions.

Step 2: In this step, for each of the answers to the problem, a new answer is created using the relationship:

[TeX:] v_{i.j} = x_{i.j} + \varphi_{i.j}(x_{i.j} - x_{k.j})$v_{i.j}=x_{i.j}+{\phi }_{i.j}(x_{i.j}-x_{k.j})$ (7)

[TeX:] i\text{ϵ\ }\left\{ 1.\ 2\ .\ .\ .\ \text{BN} \right\}$i\mathrm{\text{ϵ}}\left\{1.2...\mathrm{\text{BN}}\right\}$

[TeX:] j\text{ϵ\ }\left\{ 1.\ 2\ .\ .\ .\ D \right\}$j\mathrm{\text{ϵ}}\left\{1.2...D\right\}$

[TeX:] k\text{ϵ\ }\left\{ 1.\ 2\ .\ .\ .\ \text{BN} \right\}\ \&\ k \neq i$k\mathrm{\text{ϵ}}\left\{1.2...\mathrm{\text{BN}}\right\}\&k\ne i$

[TeX:] \text{φϵ\ }\left\lbrack - 1.1 \right\rbrack$\mathrm{\text{φϵ}}[-1.1]$

Where x_ij is the j^th parameter of the i^th solution of the problem, v_ij is the j^th parameter of the new solution, i is a number from one to the number of solutions to the problem, φ is a random number in the range of -1 to 1, k is a random number from one to the number of solutions to the problem, BN is the number of initial solutions to the problem, and D is the number of optimization parameters.

After creating a new solution, if the value of this solution is greater than the value of the previous solution, it will be replaced; otherwise, this solution will be forgotten.

Step 3: In this step, the probability of receiving a bee from each source is calculated using the following formula:

[TeX:] p_{i} = \frac{\text{fit}_{i}}{\sum_{n = 1}^{\text{SN}}\text{fit}_{n}}$p_i=\frac{{\mathrm{\text{fit}}}_i}{{\sum }_{n=1}^{\mathrm{\text{SN}}}{\mathrm{\text{fit}}}_n}$ (8)

Where fit_i is the fitness of source i (varies depending on the type and size of the problem and must be determined by the user) and p_i is the probability of selecting source i by observer bees. According to the fitness, a number of bees are assigned to each source. After calculating the value of each source using Equation 7, a new answer is created for the selected answers. If this answer has a higher value than the previous answer, this answer replaces the previous answer, and otherwise, a penalty is imposed. The purpose of the penalty is to create a counter for the number of non-improvement answers, and if the answer does not improve, one unit is added to its value.

Step 4: In this step, if the counter of the number of non-improvement answers reaches a predetermined limit (C_max), this answer will be replaced with a random answer. Also, in this step, the conditions for the end of the iterations are checked. If the termination conditions of the algorithm are met, the iterations will end, otherwise, it will return to step two.

The implementation of the Honeybee algorithm in this study involved setting several key parameters to ensure robust optimisation. The algorithm was run with a population size of 50 bees for 100 iterations. The fitness function was defined as the inverse of the flyrock distance predicted by the optimised ANN model, aiming to minimise this output. The termination condition was set to the maximum number of iterations. The search space for each blast parameter was constrained within the observed minimum and maximum values from the collected dataset presented in Table 1.

5. Conventional Predictor

Conventional empirical methods are used as a basis for predicting flyrock distance in surface blasting operations. These models are developed based on field data and simplified mathematical relationships, each considering specific parameters. Some of the most important conventional models are:

6. Results and Discussion

6.1 Prediction by ANN

The artificial neural network (ANN) model demonstrated robust predictive performance for flyrock distance. A total of 14 neural network configurations were evaluated, varying the number of neurons in the hidden layer (2 to 15). The best-performing model (Model 7) achieved a coefficient of determination (R²) of 0.8930 and 0.8874 for the training and testing phases, respectively, with root mean square error (RMSE) values of 0.2486 and 0.2512. These results highlight the model’s ability to generalise well to unseen data, indicating strong nonlinear approximation capabilities. The selected ANN structure consisted of 7 input parameters (blast hole diameter, length, spacing, burden, stemming, specific charge, and sub-drilling depth), a single hidden layer with 8 neurons, and one output (flyrock distance). Sensitivity analysis revealed that burden (B) and stemming length were the most influential parameters, contributing to 70.8% of the variance in flyrock predictions. Table 2 shows the results of this neural network analysis.

Table 2. Predicted Values for Flyrock

Figure4. Predicted values of Flyrocks for training Model No. 7

Figure5. Predicted values of Flyrocks for testing Model No. 7

6.2. Prediction by Conventional Predictor

The evaluation of traditional flyrock prediction methods, by Lundborg (1981), revealed significant limitations in their accuracy. For the Lundborg model, the root mean square error (RMSE) was calculated as 121.1 meters, with a coefficient of determination (R²) of 0.52. This indicates that the model explains only 52% of the variance in flyrock distances and suffers from oversimplification, as it ignores critical parameters like rock strength and explosive energy. For example, in sample predictions, the model estimated a flyrock distance of 110.1 meters for an actual value of 50 meters, highlighting its unreliability. Figure 6 presents a comparative analysis of the Lundborg method versus actual flyrock measurements.

Figure6. Comparison graph between the Lundborg method and actual flyrock

6.3. Flyrock Minimization using the Honeybee Algorithm

As previously mentioned, the search process continues until the minimum flyrock value is identified. Multiple iterations of the Honeybee algorithm were executed using different population sizes of bees. Figure 7 shows the execution of the algorithm in minimizing flyrock distance.

Figure7. Minimization of the amount of Flyrocks throw in meters

Based on the results, the minimized flyrock distance was determined to be 7.25 meters. As indicated in Table 1, the initial minimum flyrock distance was approximately 10 meters. By implementing the Honeybee optimisation algorithm, this value was reduced by 27.5%, demonstrating a significant improvement. Additionally, optimal values were obtained for key blasting parameters, including blast hole diameter and length, sub-drilling depth, burden, spacing, and stemming, as presented in Table 3.

Table 3. Blasting pattern optimisation values

6.4. Comparative Analysis

This paper compares flyrock prediction in mining blasting operations using two intelligent methods: Artificial Neural Networks (ANN) and Honeybee Optimization Algorithm, with traditional models. The ANN model demonstrated better performance compared to conventional models. Specifically, the ANN model achieved R² values of 0.8930 for training and 0.8874 for testing, indicating high predictive accuracy and effective simulation of flyrock behaviour. In contrast, the traditional Lundborg model, based on simplified mathematical relationships, only explained 52% of the variance in flyrock distance and showed poor prediction performance.

Furthermore, a statistical t-test on the prediction errors confirmed that the lower RMSE of the ANN model compared to the Lundborg model is statistically significant (p < 0.05), validating its superior performance.

Moreover, this model overlooks important parameters such as rock strength and explosive energy, which can lead to inaccurate results. The Honeybee Optimization Algorithm was then applied to optimize blast parameters and reduce flyrock. The algorithm optimized parameters like blast hole diameter, length, sub-drilling depth, spacing, and stemming length, reducing the flyrock distance from 10 meters to 7.25 meters, which represents a 27.5% improvement. This demonstrates that the Honeybee algorithm, using modern optimization methods, can significantly reduce flyrock distance. Sensitivity analysis in the study revealed that the burden parameter has the greatest impact on reducing flyrock, while stemming and specific charge have the most influence on increasing it.

Compared to traditional models, which rely on simplified empirical relationships, intelligent methods like ANN and Honeybee can provide more accurate predictions and more effective optimization of the blasting process in mines. These results show that intelligent techniques can significantly help in controlling flyrock and improving safety and productivity in mining operations.

Figure 8. Comparative Analysis of Flyrock Prediction Models

7. Evaluating the Influence of Blasting Parameters on Flyrock

Several advanced analytical methods, namely, Pearson correlation analysis, multiple regression modelling, random forest regression, and permutation importance analysis, were applied to evaluate the impact of various blasting parameters on flyrock. The key findings from these analyses are as follows:

7.1. Analysis with Pearson correlation

Figure 9 shows a heatmap illustrating the Pearson correlation between input parameters and flyrock distance. Among these parameters, burden (B) exhibits the most significant negative correlation with flyrock, indicating that an increase in burden leads to a considerable reduction in flyrock distance. Conversely, stemming length and specific charge show the strongest positive correlation, suggesting that an increase in these parameters results in a greater flyrock distance.

Figure 9. Pearson Correlation Heatmap

7.2. Multiple Regression Model

To quantify the contribution of each parameter, a multiple regression model was established (see Figure 10). The model achieved an R² value of 0.708, indicating that approximately 70.8% of the variations in flyrock can be explained by the selected variables. Among these parameters, burden exhibited the strongest negative correlation with flyrock, meaning an increase in burden leads to reduced flyrock distance. Conversely, stemming length and specific charge demonstrated a significant positive correlation, indicating that higher values of these parameters contribute to increased flyrock distance.

Figure 10. Multiple Regression Analysis: Coefficients with Confidence Intervals

7.3. Feature Importance Analysis (Random Forest Model)

To validate these results, a feature importance analysis was conducted using the Random Forest model (see Figure 11). The analysis confirmed that burden, stemming, and specific charge are indeed the most influential factors affecting flyrock behaviour. The results suggest that optimising Burden can significantly reduce excessive flyrock, while proper adjustments of stemming length and specific charge ensure a more controlled energy distribution.

Figure 11. Feature Importance (Random Forest)

7.4. Permutation Importance Analysis

In order to improve the sensitivity analysis, a permutation importance method was used (see Figure 12). Such a method systematically changes each variable with its value to see its significance on flyrock prediction. The findings are consistent with burden as the most impactful variable in reducing flyrock, with stemming and specific charge causing the distance of flyrock, respectively.

Figure 12. Permutation Feature Importance

The results of the various sensitivity analyses consistently show that burden, stemming, and specific charge are the three most crucial parameters influencing flyrock. However, although burden as a major factor helps fight the rock throw, stemming and specific charge increase flyrock and therefore require key optimisation.

These results have been statistically confirmed using a multiple regression model and correlation analysis, while Random Forest and Permutation Importance validated this systematic and well-established finding as machine-learning-based. The congruence of these methods further validates the results and allows for confirmation that changes to burden, stemming, and specific charge in the blast design can greatly enhance flyrock control.

To consolidate the findings from the different analytical techniques, Table 4 provides a summary of the parameter importance rankings. The consistency across all four methods strongly validates the identification of burden, stemming, and specific charge as the most critical parameters governing flyrock.

Table 4. Consolidated Summary of Sensitivity Analysis Results

8. Conclusions

This study successfully demonstrated the development and application of a hybrid intelligent system, combining an Artificial Neural Network (ANN) with a Honeybee Optimisation Algorithm, for the prediction and minimisation of flyrock in surface mining. The key novelty of this research lies in its integrated approach, moving beyond simple prediction to provide optimised, actionable blast design parameters for enhanced safety.

The developed ANN model, selected through a systematic scoring method, showed excellent predictive capabilities with high R² values (0.8930 for training, 0.8874 for testing) and low RMSE values, significantly outperforming traditional empirical methods. A comprehensive sensitivity analysis consistently identified burden, stemming, and specific charge as the most influential parameters affecting flyrock distance. The practical implication of this work is significant. By applying the Honeybee algorithm to the validated ANN model, blast parameters were optimised to achieve a minimum flyrock distance of 7.25 meters, representing a 27.5% reduction from the safest observed blasts in the field data. This result confirms that data-driven optimisation can lead to considerable improvements in operational safety and efficiency.

For future research, it is recommended to apply this hybrid model to other mine sites with different geological conditions to test its robustness. Furthermore, incorporating geomechanical variables (e.g. rock mass rating, joint properties) as inputs could further enhance the model's accuracy. Finally, the performance of the Honeybee algorithm could be benchmarked against other state-of-the-art metaheuristic algorithms to explore further potential for optimisation.

Author’s contribution

Mojtaba Rezakhah: conceptualization, data curation, formal analysis, investigation, resources, software, supervision, and writing – original draft. Erfan Nemati: formal analysis, investigation, methodology, software, validation, and writing – original draft. Amir Batarbiat: formal analysis, investigation, methodology, resources, software, validation, visualization, and writing – original draft. Manoj Khandelwal: data curation, methodology, resources, supervision, validation, visualization, and writing – review & editing

All authors have read and agreed to the published version of the manuscript.