Prediction of blast-induced air overpressure using a hybrid machine learning model and gene expression programming (GEP): A case study from an iron ore mine

: Mine blasting can have a destructive effect on the environment. Among these effects, air overpressure (AOp) is a major concern. Therefore, a careful assessment of the AOp intensity should be conducted before any blasting operation in order to minimize the associated environmental detriment. Several empirical models have been established to predict and control AOp. However, the current empirical methods have many limitations, including low accuracy, poor generalizability, consideration only of linear relationships among influencing parameters, and investigation of only a few influencing parameters. Thus, the current research presents a hybrid model which combines an extreme gradient boosting algorithm (XGB) with grey wolf optimization (GWO) for accurately predicting AOp. Furthermore, an empirical model and gene expression programming (GEP) were used to assess the validity of the hybrid model (XGB-GWO). An analysis of 66 blastings with their corresponding AOp values and influential parameters was conducted to achieve the goals of this research. The efficiency of AOp prediction methods was evaluated in terms of mean absolute error (MAE), coefficient of determination (R 2 ), and root mean square error (RMSE). Based on the calculations, the XGB-GWO model has performed as well as the empirical and GEP models. Next, the most significant parameters for predicting AOp were determined using a sensitivity analysis. Based

by engineers.Uncontrollable parameters in the second group consist of the rock mass's RQD and compressive, shear and tensile strengths.These mainly depend on the complex geological and geotechnical conditions of the rock mass.
Remennikov and Rose [10] have proposed solutions such as improving the structure of buildings and glass doors or using barriers as shields to reduce the effects of AOp.However, their efficacy was not significantly proven [11].With similar objectives, experimental methods are also presented to calculate AOp.Also, Khandelwal and Kankar [12] and Armaghani et al. [13] showed that empirical equations have poor performance.In addition to low accuracy, conventional tools also have many limitations including studying specific areas, considering linear relationships between influence parameters and focusing on explosive charges per delay and monitoring distances [13].
Over the years, machine learning and artificial intelligence (AI) methods have been demonstrated to be acceptable and reliable for solving different engineering subjects, especially for prediction and optimization purposes [14][15][16][17][18][19].For example, AI can be used to analyze the results of a blast, using machine learning algorithms to identify patterns, and then this data can be used to adjust the blast design for future operations [20,21].Prediction of AOp has also been conducted using these methods.The artificial neural network (ANN) proposed by Khandelwal and Singh [22] was compared with regression methods to predict AOp.Results showed that ANN performed similarly to regression methods in predicting AOp with higher precision.Khandelwal and Kankar [12] for predicting AOp investigated and confirmed the efficacy of support vector machine (SVM).Nguyen et al. [23] evaluated and compared different ANN systems including ANN, BRNN and HYFIS to predict blastinduced AOp.Similarly, a heuristic algorithm, intended to optimize ANNs for predicting AOp, was used successfully by Armaghani et al. [24].Also, Hasanipanah et al. [25] developed a hybrid SVR model that was optimized with PSO to predict AOp.
Several soft computing models were presented by Bui et al. [4] for predicting AOp, including SVM, ANN, boosted regression trees, k-nearest neighbors and RF.In another study, Zhou et al. [26] developed a hybrid model which combines a fuzzy system (FS) and the firefly algorithm (FFA).Their study demonstrated that FFA-FS could be used to predict AOp efficiently.Tran et al. [27] have investigated the effect of meteorology on AOp, using AOp prediction models.Their results show that meteorological conditions, especially wind speed and air humidity, have a noteworthy impact on blastinduced AOp.Zeng et al. [28] have developed an efficient method based on the Levenberg-Marquardt (LM) algorithm and cascaded forward neural network (CFNN) to predict AOp.In addition, the accuracy level of this technique has been tested using the generalized regression neural network (GRNN) and extreme learning machine (ELM).Table 1 summarizes some of the relevant previous studies conducted by various researchers.
Blasting operations have a range of effects on the environment, including air overpressure (AOp) [29][30][31].Therefore, AOp prediction with high accuracy is essential to determining the safe zone around an operation site.In many studies, machine learning has been applied to evaluate and predict the adverse consequences of blasting.However, the studies have not addressed the evaluation and prediction of air overpressure via hybrid extreme gradient boosting (XGB).Thus, this article develops a hybrid model which combines extreme gradient boosting (XGB) with grey wolf optimization (GWO) to predict AOp in open pit iron mines accurately.Furthermore, gene expression programming (GEP) and an empirical model were used to assess the validity of the hybrid model (XGB-GWO).Additionally, unlike other methods, GEP can create a simple mathematical expression that can be applied to different mining conditions, which is another advantage of this article.

Methodology
Machine learning models require parameters based on datasets.Many studies have shown that heuristic algorithms improve machine learning accuracy and stability [40,41].Therefore, this study applied a hybrid predictive approach based on ideas from the extreme gradient boosting framework (XGB) as well as metaheuristic algorithms: namely, grey wolf optimization (GWO).Then, gene expression programming (GEP) and empirical models were employed to assess the validity of the optimized model (XGB-GWO).By using GWO we can find the optimal values of hyperparameters of the regression model.Therefore, the intelligent optimization algorithm by the adjustment of three key parameters of the XGB model (n_estimator, maximum_depth and learning_rate) achieves higher accuracy.N_estimator, maximum_depth and learning_rate represent the number of trees, the maximum depth of a tree and the shrinkage coefficient of the tree, respectively.

Extreme gradient boosting (XGB)
Chen et al. [42] proposed a method based on gradient boosting [43][44][45][46].Many engineering fields have implemented XGBoost for classifying and predicting problems [47].It has performed well because of the advantages of parallel processing, regularization and efficient tree pruning.XGBoost can solve a wide variety of data science problems quickly and accurately with parallel tree boosting.The core of this algorithm is optimizing the objective function [48].
In each iteration, XGBoost calibrates the previous predictor using the residual.Loss function optimization (LOF) is involved in this process.However, regularization is applied to the objective function during calibration to reduce overfitting.Using this description, Equation (1) describes the objective function as two parts: regularization and training loss.
In Equation ( 1), Θ and Ω represent the parameters trained from the data and related to regularization, respectively.Regularization controls the complexity of the model in order to avoid overfitting [49].L denotes the training loss function, which measures the model's fit to the training data.Complexity can be defined in a variety of ways.Equation (2), however, is often used to calculate the complexity of each tree.
where γ denotes each leaf's complexity, T represents the number of DT leaves, ω is the vector of scores on the leaves and λ scales the penalty.Next, XGBoost combines loss function (LOF) general gradient boosting with the second-order Taylor expansion.When the mean squared error is taken into account, equation ( 3) can obtain the objective function.
The MSE loss function's first and second derivatives are   and ℎ  , respectively.Also, the q function is used to assign data points to leaves.As a final step, we calculate the XGBoost objective function by referring to Equation (4).
Here   are independent of each other.The definitions of the two terms   and   are given in Equation (5).
In general, finding the minimum of a quadratic function can be applied to optimizing the objective function.In order to evaluate the performance of the model after splitting a particular node in DT, an objective function is used.If the model performs better than before, this division will be accepted.Otherwise, the division will come to a halt.Due to the inclusion of regularization phenomena, XGBoost can prevent over-installation more effectively [50].

Grey wolf optimization (GWO)
GWO is a metaheuristic optimization algorithm proposed by Mirjalili et al. [51] which reflects grey wolf family social systems in nature.Social hierarchy is very strict in the grey wolf community.Usually, wolf levels are divided into four categories: α, β, δ and ω. α is primarily responsible for making overall decisions at the first level, and the β wolf assists the α wolf at the second level.At the third level, δ must follow the decisions of α and β wolves.Furthermore, the lowest rank for wolves is ω.Also, this group must follow wolves with a higher rank.Mathematical models can be developed, and optimal solutions can be found based on wolves' hunting process and social hierarchy.
α represents the optimal solution, β denotes the second, and δ shows the third optimal solution.All other applicant solutions are denoted by ω.In the wolf pack, D indicates the distance between the prey and the individual (Equations ( 6) and ( 7)).
where C, D and XP(t) are the step length coefficient, the current number of iterations and prey location, respectively.Also, r1 and X are a random number ranging in (0,1) and the location of a grey wolf, respectively.
A pack of individuals tries to shorten the distance between their prey and themselves by constantly updating the following Equations ( 8) and ( 9): where r2 and A are a random number ranging in (0,1) and the convergence influence factor, respectively.A decreases linearly with the number of iterations from 2 to 0. Since α, β and δ possess high levels of intelligence, they can carry more information about where prey can be found.This will enable them to lead the group to the hunting grounds gradually.Gradually, the pack makes its way to the hunting grounds.As a result, three optimal solutions will now be available to ignore or find other solutions, and gradually they will find the global optimal solution based on the three optimal solutions, which are outlined in Equations ( 10) to (12).
A decision will be made as to which individuals to include in the remaining pack relying on a joint decision by α, β and δ.It is now necessary to move the position shown in Equation (13).
It is possible to summarize the grey wolf optimization algorithm in this way: during the optimization process, it constantly updates the solution area for the location search problem and finds the ideal solution at the end (Figure 1).

Gene expression programming (GEP)
The GEP method was introduced by Ferreira [52] and combines genetic programming (GP) and genetic algorithms (GA).As an AI developmental approach, GEP has corrected some of GA's and GP's weaknesses, such as tree systems.Both of these methods differ in the following ways: the concept of individuals or answers.In the GEP algorithm, individuals are defined as binary fixed-length chromosomes.Answers from the GP method can include tree systems of varying sizes.Because GEP combines GA and GP, it produces chromosomes of fixed length and tree systems with different sizes and shapes, identified as an expression tree (ET).The GEP algorithm's structure consists of several elements, including function, terminal, operator, fitness, and stop criteria [52].
Chromosomes are composed of two parts: a head and a tail, making them a fixed length.Functions and terminals are included in the head part, and terminals are also included in the end part.As a result of the complexity of the problems, a specific equation cannot be used to calculate the head part's length.This term is defined as the input of the GEP method.Trial and error are one of the solutions [53].Furthermore, the tail part's length can be calculated using Equation (14).
where   , h and t represent the number of arguments of the functions and the length of the head and tail.Solitary chromosomes in the initial population are evaluated according to a fitness function developed for gap problems.A number of genetic operators are used to adapt the considered chromosomes.Depending on the problem conditions, each chromosome may contain functions, constants and terminals [54,55].The following is a general description of the GEP algorithm: Step 1: Depending on the problem's conditions (size) being studied, a certain number of chromosomes should be determined (randomly).
Step 2: Expression trees and mathematical equations are used to select the initial population chromosomes.
Step 3: The chromosomes are fitted according to the overall fitness function (RMSE or R 2 ).Alternative methods, such as the roulette wheel method, are used if the stopping criterion is not met.
Step 4: During this step, the GEP algorithm's genetic operators, known as the core, must be linked to the rest of the chromosomes.
Step 5: Lastly, the process of creating the next generation begins, and the process is repeated until new structures are created.
In order to decode the programs on the chromosomes, Karva (K-Expression) was invented to express the codes.Inversions, mutations, triple recombination operators and triple transposition operators have been introduced so far as genetic operators that are used for chromosome modification [54].The structure of the GEP method is illustrated in Figure 2.
The equation recommended by USBM has been widely adopted by researchers and is expressed in the following manner [57]: By regression analysis, AOp in decibels (dB) is calculated, and  and  are site factors.

Model verification and evaluation
Evaluation and verification of models are essential steps during model development.In order to determine whether a model is of high quality and whether the results produced are adequate for the goals pursued, it is necessary to test it.In this study, training and testing sets are used to train and verify predictive models.Predicted and actual values are compared using relevant evaluation indicators.
In this study, hybrid models will be evaluated to ensure their reliability [58].R 2 , RMSE, and MAE are evaluation indicators.
As defined in Equation 17, between actual and predicted values, R 2 represents the square of the correlation.Further, the RMSE represents the standard deviation between predicted and actual values (Equation 18).Equation 19defines mean absolute error as the measure of error between paired observations describing the same phenomenon [59].

Case study and data collection
This study analyzes the case of Chadormalu, a large open-pit iron mine located 120 kilometers northeast of Yazd, Iran. Figure 3   In the current study, input and output parameters were used in developing the model.The influential input parameters used for air overpressure (AOp) prediction are given in Table 2. Burden (B), number of holes (NH), hole depth (H), spacing (S), powder factor (PF), distance (D), rock quality designation (RQD), stemming (ST) and the maximum charge per delay (MC) have been employed as the output parameters of the models to predict air overpressure (AOp).
Data were collected as follows: D was determined by handheld GPS, and AOp values from blasting operations were recorded with a sound level meter (SLM).Using blast design, the remaining parameters were collected.Table 2 shows the influential input and output parameters' details, symbols and statistical descriptions.
For model development, 66 data pairs were used to create a database related to air overpressure (AOp).From the organized database, 20% of the data set was selected to test the model to ensure consistency.

Results and discussion
Predicting air overpressure requires the preparation of a database.This database was divided into training and test sets using the most common division ratio of 80%-20% [60].The air overpressure prediction model was evaluated using several performance indicators, including R 2 , RMSE and MAE.All predictive models use the same training and test data sets.

Hybrid model result (XGB-GWO)
In order to avoid complexity in XGBoost modeling, three stopping criteria were considered, namely, n estimators, learning rate and maximum depth.If significant values are assessed for each parameter, overfitting can occur.Therefore, XGB parameters have been optimally determined using grey wolf optimization (GWO).Figure 1 illustrates the method used to develop models based on the XGB.It was time to set the XGB model's parameters.The optimization algorithm's parameters are shown in Table 3.Additionally, Table 3 shows the optimal parameters of the model based on the optimization process.Table 3.The optimization algorithm's parameter and optimal parameters for the hybrid model.

GEP model result
The flowchart shown in Figure 2 illustrates the gene expression programming (GEP) modeling process.The GEP was designed using the same testing and training datasets as the previous sections.The final relationship between the initial data and air overpressure was analyzed and determined by GeneXproTools (version 5.0).
To build an efficient model, it is imperative to consider the fitting parameters.The number of chromosomes determines how long it takes for the model to run, which is crucial to the model's performance.For a suitable architecture, the size of the head and the number of genes must also be considered.Each component's complexity and the number of related equations are determined by the head size and the number of genes.To achieve the most suitable GEP model, a trial and error mechanism was applied (Table 4).By using R 2 , RMSE and MAE indices, for the testing and training data sets, the performance prediction of GEP models was evaluated.Of the ten models stated in Table 4, the five which had the most accurate predicting of air overpressure were selected (Table 5).Therefore, based on the results in Table 5, model No. 6 was the most accurate model among all models made by the GEP method.Figure 5 illustrates the scatter plot of the predicted and measured air overpressure in the selected GEP model by training and test data.The number of genes and head size in model No. 6 are 8 and 4, respectively.The expression tree (sub-ETs of each gene) of this model is illustrated in Figure 6.In addition to the expression trees, functions were linked together to construct a large and complex ET.Expressions (20) to (23) are mathematical expressions of each of the genes.Finally, the developed GEP equation for predicting air overpressure is shown in Equation (24).

Empirical model result
In the empirical model, site factors  and  were computed from 52 blasting events (training data set).A multivariate regression analysis was performed using Microsoft Excel 2019.Therefore, the optimal values of  and  for the USBM model for predicting AOp are 173 and 0.1.Using the USBM model in this case, we can describe it as follows: Also, Figure 7 illustrates the scatter plot of the predicted and measured air overpressure in the empirical (USBM) model by training and test data.

Comparison of models and validation performance
The efficiency of the predictive models is evaluated in this section.In this study, R 2 , RMSE and MAE (Equations 15 to 17) were employed to assess predictive models' performance.Regarding Table 6, the above statistical criteria and performance comparison for techniques were determined for the testing and training data sets.Based on these results, the correctness rate and performance of the hybrid XGBoost (XGB-GWO) method are preferable to the GEP and empirical models.In the next step, the selected AOp prediction models' accuracies are compared, as shown in Figure 8.The accomplishment of the models in predicting the air overpressure in the testing and training data set is shown in Figure 8. Also, based on Figure 8, the GWO-XGB technique gives the most reliable and steady results in AOp prediction among the GEP and empirical models.As shown in Figure 9, the box plots show the distribution functions corresponding to the measured and predicted AOp values.Based on an exhaustive analysis of Figure 9, the XGB-GWO approach yielded the most promising performance relative to the GEP model due to its similar probability distribution to observational results.A Taylor diagram illustrating the performance of the predictive models is presented in this subsection.Taylor diagrams are used to assess the accuracy of models by showing them in two dimensions [61].Indicators of the relationship between the actual and predicted observations are R, RMSE and standard deviation.Each model is denoted by a term in the Taylor diagram.In an ideal model, the position of the point should coincide with the reference point.Figure 10 illustrates the predictive models developed in this study.While both models are highly accurate at predicting air overpressure, Figure 9 shows that the hybrid XGB-GWO model is better.

Decay of air overpressure
Understanding the decay of overpressures during the explosion of an explosive charge in an open pit mine is crucial in ensuring the safety of personnel and equipment.The pressure waves generated during an explosion can cause significant damage to structures and equipment, as well as pose a risk to the health and safety of personnel.By examining the decay of overpressures, it is possible to determine the safe distance for personnel and equipment from the explosion site, as well as the appropriate level of protective measures required.This information can also inform the design of blast patterns and other explosives-related procedures in order to minimize the risk of harm.Therefore, studying the decay of overpressures during explosions is an important aspect of ensuring the safe and effective use of explosives in open pit mines [62].
In this case, only two distances have been measured, and it is necessary to use empirical and GEP methods to estimate the decay of overpressures at different distances.The resulting plot can provide valuable insights into the behavior of air overpressure as it decays over distance, allowing us to identify any unexpected trends or anomalies that may require further investigation.In this regard, Figure 11   Analyzing the plot, we can see that the measured AOp values are highest at the closest distance to the explosion site (300 m) and decrease as the distance increases.Both the empirical and GEP predicted AOp values follow similar trends, with the GEP method predicting AOp values first that are generally higher than the empirical values.The plot also shows that the GEP method generally provides a reasonable estimate of the AOp values at different distances from the explosion site, although it tends to underestimate the AOp values at shorter distances.

Sensitivity analysis
During the final modeling stage, the output of the model is analyzed for its sensitivity to the input parameters.Analyzing the sensitivity of input parameters can provide insights into how they affect model output (objective function).For determining sensitivity analysis, the cosine amplitude method (CAM) is one method [63][64][65].Equation ( 26) describes the CAM method in terms of an n-dimensional space.
The m-dimensional length vector   , which is part of the array X, represents a variable in the given context.
= { 1 ,  2 ,  3 … ,   } (27) This implies that the different dimensions of   are interrelated with those of   , and the degree of correlation between them can be expressed mathematically as shown in the following equation.Input parameters have a greater influence on output as   approaches one.It has been shown that input parameters have a significant effect on the output when this parameter is above 0.9 [63].In Figure 12, the regression prediction parameters are sensitivity analyzed.According to Figure 10, stemming length and RQD had the most significant effect on the air overpressure among the input elements.

Conclusions
Blasting is a commonly used method in open-pit mines for breaking down rocks.However, it can result in adverse effects, such as air over-pressure (AOp), ground vibration, flyrock, backbreak, and dust, which can have a negative impact on the surrounding environment.To address this issue, it is essential to predict and control the effects of blasting.In this study, an efficient and practical hybrid machine learning model (XGB-GWO) was proposed for predicting AOp values, and its performance was evaluated using gene expression programming (GEP) and an empirical model.The accuracy of the predictive models was assessed using standard criteria such as R 2 , RMSE and MAE.Additionally, a sensitivity analysis was conducted using the cosine amplitude method (CAM) to determine the intensity of input parameters at AOp.Overall, this study highlights the significance of developing effective models to predict the impact of blasting and minimize its adverse effects on the environment.In conclusion, the hybrid XGB models proposed in this study demonstrate considerable potential in predicting AOp and can aid XGB in adjusting hyperparameters.The performance of the predictive models for test data falls within the following ranges: XGB-GWO (MAE: 0.69; RMSE: 1.42; R2: 0.983), GEP (MAE: 0.63; RMSE: 1.04; R2: 0.989) and empirical (MAE: 5.92; RMSE: 6.68; R2: 0.53).It is worth noting that the XGB-GWO hybrid model outperforms the other models in terms of overall performance.Furthermore, a sensitivity analysis technique called cosine amplitude was used to determine the importance of input variables, which revealed that the stemming length and RQD had the most significant impact on the penetration rate.This study demonstrates that the proposed XGB-GWO model is robust and effective in predicting blast-induced AOp, indicating its potential for practical applications in the mining industry.

Figure 1 .
Figure 1.Structure of hybrid method based on XGB and GWO for predicting air overpressure.

Figure 2 .
Figure 2. Structure of the GEP method for prediction of air overpressure.

2. 4 .
EmpiricalIn open-pit mines, the blast-induced AOp has been estimated using empirical methods.The technique proposed by the United States Bureau of Mines (USBM) is the most widely used empirical equation for predicting AOp among all available methods.Kuzu et al. [56] calculated the scaleddistance (SD) factor for AOp predictions.USBM is an empirical method for predicting AOp based on site factors, the maximum charge per delay (MC) and monitoring distance (D) in open-pit mines.SD values determine the relationship between MC and D [56]: depicts the location of the Chadormalu mine in Iran in relation to Tehran (the capital of Iran) and Yazd.It is estimated that there are about 400 million tons of mineable ore reserves.Magnetite and hematite are the main components of the deposit, according to mineralogical studies.Mine blasting is primarily carried out with ANFO explosives.

Figure 3 .
Figure 3.The location of the Chadormalu mine in Iran relative to Tehran (the capital) and Yazd city.
training data set, Figure 4 shows the correlation between the actual values and the predicted values.The training data are largely scattered near the closest fit line for the intelligent model, showing that training effects remain relatively favorable.XGB-based optimization techniques have been proposed that demonstrate high training effects with R 2 values of 0.96.After the hybrid intelligent model has been trained, it is verified and evaluated using the testing data set.Figure4shows that the test data set is mostly distributed near the perfectly fitted line.Because the predicted AOp values and the actual AOp values are within the correlation, they can be classified as correlated.Hybrid models achieve high levels of prediction accuracy, with R 2 values of 0.98.

Figure 4 .
Figure 4. Actual values and predicted values on the testing and training dataset by XGB-GWO Method.

Figure 5 .
Figure 5. Measured and predicted air overpressure on the training and test dataset by selected GEP Method.

1 (( 6 .Figure 6 .
Figure 6.Sub-ETs of each gene for the best GEP model with addition as a linking function.

Figure 7 .
Figure 7. Measured and predicted air overpressure on the training and test dataset by Empirical Method.

Figure 8 .
Figure 8. Predictive models with prediction values on training and testing datasets.

Figure 9 .
Figure 9.The box plot of predictive models (distribution function).

Figure 10 .
Figure 10.Performance comparison of predictive models with testing and training data in Taylor diagrams.
plots the measured, empirical and GEP predicted AOp values for different distances from the blasting site.

Figure 11 .
Figure 11.Comparison of measured, empirical, and GEP predicted AOp values at different distances from the blasting site.

Figure 12 .
Figure 12.The sensitivity analysis and effect of input data on air overpressure.

Table 1 .
Some recent work with machine learning techniques for air overpressure prediction.
Note: Symbols are explained in the Abbreviations section.

Table 2 .
Input and output data with details, symbols and statistical descriptions.

Table 5 .
Performance indices of GEP models.

Table 6 .
Statistics of predictive models and performance comparison.