Prediction of blasting mean fragment size using support vector regression combined with five optimization algorithms

Enming Li, Fnghao Yang, Mihng Rn, Xiliang Zhang, Jian Zhou,Manoj Khanlwal

a School of Resources and Safety Engineering, Central South University, Changsha, 410083, China

b Universidad Politécnica de Madrid - ETSI Minas y Energía, Ríos Rosas 21, Madrid, 28003, Spain

c College of Locomotive and Rolling Stock Engineering, Dalian Jiaotong University, Dalian,116028, China

d Department of Mathematics, University of California Santa Barbara, Santa Barbara, CA, 93106, USA

e State Key Laboratory of Safety and Health for Metal Mines, Maanshan, 243000, China

f School of Engineering, Information Technology and Physical Sciences, Federation University Australia, Ballarat, VIC, 3350, Australia

Keywords:Blasting mean fragment size e-support vector regression (e-SVR)V-support vector regression (v-SVR)Meta-heuristic algorithms Intelligent prediction

ABSTRACT The main purpose of blasting operation is to produce desired and optimum mean size rock fragments.Smaller or fine fragments cause the loss of ore during loading and transportation, whereas large or coarser fragments need to be further processed, which enhances production cost. Therefore, accurate prediction of rock fragmentation is crucial in blasting operations. Mean fragment size (MFS) is a crucial index that measures the goodness of blasting designs. Over the past decades, various models have been proposed to evaluate and predict blasting fragmentation.Among these models,artificial intelligence(AI)-based models are becoming more popular due to their outstanding prediction results for multiinfluential factors. In this study, support vector regression (SVR) techniques are adopted as the basic prediction tools, and five types of optimization algorithms, i.e. grid search (GS), grey wolf optimization(GWO), particle swarm optimization (PSO), genetic algorithm (GA) and salp swarm algorithm (SSA), are implemented to improve the prediction performance and optimize the hyper-parameters.The prediction model involves 19 influential factors that constitute a comprehensive blasting MFS evaluation system based on AI techniques. Among all the models, the GWO-v-SVR-based model shows the best comprehensive performance in predicting MFS in blasting operation. Three types of mathematical indices, i.e.mean square error(MSE),coefficient of determination(R2)and variance accounted for(VAF),are utilized for evaluating the performance of different prediction models. The R2, MSE and VAF values for the training set are 0.8355,0.00138 and 80.98,respectively,whereas 0.8353,0.00348 and 82.41,respectively for the testing set. Finally, sensitivity analysis is performed to understand the influence of input parameters on MFS. It shows that the most sensitive factor in blasting MFS is the uniaxial compressive strength.

1. Introduction

Blasting is widely used in many engineering operations including mining engineering, civil engineering and tunneling(Sanchidrián et al.,2007;Wang et al.,2018a,b;Hu et al.,2020;Zhou et al.,2020a,2021a).The primary purpose of the blasting operation is fragmentation and displacement of rock mass in mining engineering.Although significant developments have been achieved in explosive technology, only 20%-30% explosive energy is used for the actual fragmentation and displacement of rock mass(Khandelwal and Monjezi,2013;Ebrahimi et al.,2016),and the rest of the energy is dissipated in the ground or air and produces various hazardous effects,such as backbreak(Monjezi et al.,2012;Esmaeili et al., 2014), ground vibration (Zhou et al., 2020b), flyrock(Armaghani et al.,2014),and air blast.In mining operations,a blastpattern must be designed with due care to obtain optimum fragment size.Fragment size plays a crucial role in subsequent crushing and grinding operations. It is well-known that there is a tradeoff between fragment size and economic benefit.Large fragments need secondary blasting to reduce boulder size; however, finer fragments will increase the cost of mining due to higher explosive charge. Therefore, controlling and predicting the blast fragment size is important in mining operations.

Table 1 Previous work about blast fragmentation prediction using AI techniques.

A number of researchers have proposed and developed various empirical models to predict and estimate the blasting fragment size(e.g. Koulli and Rustan, 1993; Aler et al., 1996a,b; Zhang and Jin,1996; Chung and Katsabanis, 2000; Hamdi et al., 2002; Thornton et al., 2002; Esen et al., 2003; Moser, 2003; Ouchterlony, 2003;Spathis, 2004; Cunningham, 2005; Gheibie et al., 2009). These empirical models were created based on plenty of blast engineering cases, blasting mechanisms and mathematical statistics. For instance,Kuznetsov(1973)proposed a prediction function of mean fragment size (MFS), where rock mass and explosive properties were used as the main factors.Bergmann et al.(1973)developed a classical model named BRW model,in which a few new factors such as the velocity of detonation and the peak pressure in the gauge hole with filled water were associated. Later on, based on the research from Holmberg (1974) and Larsson (1974), a function of MFS in terms of the specific charge,explosive strength and spacing to burden ratio was proposed.Inspired by the BRW model,Aler and Du Mouza (1996a, b), Hamdi et al. (2002) and Ouchterlony (2003)drew up an equation, so that the BRW model can be utilized easily. Chung and Katsabanis (2000) created 80% passing size and MFS prediction equations based on rock mass factors and blasting design parameters. Gheibie et al. (2009) modified the Kuz Ram model and proposed a uniformity index method. Table A1 in Appendix A illustrates a detailed description of various empirical equations and the relative work of various researchers. These studies provide a valuable significance for future work, as they integrated different blast design parameters and rock mass parameters to analyze the relationship of various influential factors with blast fragmentation.

Blast fragmentation is influenced by numbers of blast design parameters. It is difficult to propose a fitting equation considering all the influential factors simultaneously. In such complex circumstances, artificial intelligence (AI)-based methods, which can develop a complicated relationship among various input and output variables, attract the attention of scholars worldwide. Over the past few decades, AI-based methods have brought satisfactory prediction results in mining engineering and various geoengineering fields (Hajihassani et al., 2015; Gordan et al., 2016;Zhou et al.,2016,2019b;Moayedi and Armaghani,2018;Armaghani et al., 2019; Sarir et al., 2019; Zhang et al., 2020b, 2021a). For instance, Shi et al. (2012) applied a support vector machine(SVM)to assessing rock fragment distribution and compared their prediction abilities with artificial neural network (ANN), multivariate regression analysis (MVRA) and Kuznetsov methods. Gao et al.(2018) used Gaussian process regression (GPR) with five kernel functions to propose various rock fragmentation prediction models.The GPR model performance was evaluated by the coefficient of determination(R2)and also compared with other AI-based models.It was found that the GPR-squared exponential model showed the best prediction performance. Monjezi et al. (2009) applied a fuzzy inference system (FIS) to predicting rock fragmentation and compared their results with regression analysis. They found that FIS-based models performed much better than the regression models. Hasanipanah et al. (2018) explored the feasibility of particle swarm optimization-adaptive network-based FIS(PSO-ANFIS)technique to estimate rock fragmentation and compared the prediction capability of PSO-ANFIS model with SVM, ANFIS and nonlinear multiple regression (MR) models. They found that the PSO-ANFIS model outperformed other models and its prediction capability was assessed with R2and root mean square error(RMSE).Asl et al. (2018) adopted ANN and firefly algorithm (FFA) as prediction tools to optimize the flyrock distance and rock fragment distribution and found that the obtained models can offer favorable evaluation results for flyrock and fragment size evaluation. More pertinent work about blast fragmentation prediction using AI methods is tabulated in Table 1.

A remarkable outcome could be achieved by implementing these models;however,there are still some shortcomings that need to be addressed. For instance, the selection of hyper-parameters is not determined by standard methods and probably ignores more effective hyper-parameters (Dimitraki et al., 2019). In the past studies, cross-validation was not implemented and thus the prediction accuracy may not be scientific(Ebrahimi et al.,2016).Apart from that, a limited number of influential parameters were taken into consideration in the past studies to develop the blast fragmentation prediction models. From the literature review,it can be found that various significant factors were ignored and thus the prediction models failed to provide convincing results.

In this study, 76 groups of blasting fragmentation datasets are collected from Sharma and Rai(2017) to predict blasting fragment size.Here,support vector regression(SVR)models including e-SVR and v-SVR are utilized as the main prediction tools combined with five optimization algorithms, i.e. GA, PSO, grey wolf optimization(GWO), salp swarm algorithm (SSA) and grid search (GS) method.Meanwhile, cross-validations are employed to examine the prediction capability of different models.Three mathematical indices,i.e.R2,mean square error(MSE)and variance accounted for(VAF),are used to assess the prediction performance. Finally, the sensitivity analysis is implemented to understand the sensitivity of each input parameter on blasting MFS.A general working framework of this study is demonstrated in Fig.1.

2. Data description

In this study,a total of 76 groups of blasting datasets are utilized.The datasets involve 19 influential factors and one output parameter (MFS (m)). These influential factors are categorized into three groups,i.e.blast design parameters,explosive parameters and rock mass parameters. Compared with the past research, more influential parameters are used for developing prediction models. The blast design parameters are blast hole diameter (D) (m), average bench height (H) (m), average sub-grade drilling (J) (m), average spacing(S)(m),average burden(B)(m),average stemming(ST)(m),average length(L)(m),average width(Wd)(m),S/B ratio(S.B),ST/B ratio(ST.B),stiffness ratio(H.B),J/B ratio(J.B),B/D ratio(B.D),length/width ratio (L.Wd), and number of holes (NH), whereas explosive parameters are total explosive amount (Qe) (t), linear explosive density(De)(kg/m)and powder factor(PF)(kg/m3).The rock mass parameter is reflected by the UCS (MPa).

General data distribution of each of the parameters is shown in Fig.2 based on the violin plot.A violin plot consists of a boxplot and a density plot,where a black spot is a median value.The black box represents the range from the lower quartile to the upper quartile.The black line indicates the 95% confidence interval. The outer shape of the box represents the density estimation of the data.The original datasets are divided into two parts at a ratio of 4:1.The 80%datasets are used to establish the training networks(Esmaeili et al.,2015).The remaining 20%of the datasets are not used in the model development but are employed to test the prediction performance of the network.

Fig.1. A framework of SVR-based model for blasting fragment size evaluation.

3. SVR and optimal algorithms

There are numbers of factors that influence MFS in blast operation. When the number of influential parameters is too high, the conventional empirical equations hardly provide a good prediction of MFS.Therefore,in this study,an SVR tool is applied to predicting the MFS.

3.1. V-SVR and e-SVR

SVM is a powerful machine learning technique developed by Vapnik (2013) to tackle classification problems based on mathematical statistics. SVM can solve regression problems by introducing the ε-insensitive loss function (Cherkassky and Ma, 2004;Shawe-Taylor and Cristianini, 2004). For the SVR, it can be described as follows.

Suppose that a group of datasets are represented by{(x1，y1)，(xa，ya)，…，(xAA，yAA)｝,where a represents the number of training data,AA denotes the total number of training data,x∊RNrepresents that the dimension of input variables is N,and y∊R1indicates an output.If the SVR is used for handling nonlinear regression problems,then it can be written as

where ω is the weight coefficient, θ is the coefficient used for transforming the nonlinear problem into a linear problem,and b is the model error.

The v-SVR is another branch of SVM (Schölkopf et al., 2000;Chang and Lin, 2002; Thomas et al., 2017). For the v-SVR, there is one important parameter, i.e. v with the range of (0, 1]. This parameter is used for balancing the number of support vectors and model errors. The primary task of v-SVR is as follows:

where‖w‖is the norm of w,ξaand ξ*arepresent two positive slack variables used for controlling the distance between the actual values and boundary values, and C is a regularization parameter used for balancing the model errors and model flatness. The ε-insensitive loss function is used for determining if the value of ω∙θ(xa)is located in the range of y±ε,then the calculation loss can be ignored.

Fig. 2. Data distribution of all parameters used for developing prediction models.

While in the e-SVR (Chang and Lin, 2002), the aforementioned formulation is written as follows and subject to Eq. (3):

where gamma defines the bandwidth of the RBF.

There are two parameters needed to be selected to control the model performance in the RBF kernel function,i.e.penalty factor Ct and RBF kernel deviation g.The penalty factor Ct is used to balance the model complexity and biasness. A higher Ct would probably bring better model fitting while being easier to get stuck into“overfitting”.A smaller Ct would decrease the model complexity but also weaken the model performance.The RBF kernel deviation g implies the space distribution of mapped data. Generally, larger g needs fewer support vectors,while smaller g would require more support vectors.The number of support vectors determines the speed of the model development. Two typical data redistributions by hyperplane in SVR are demonstrated in Figs.3 and 4.More details about SVR theory and its implementation can be found in other available literature (Drucker et al.,1997; Smola and Schölkopf, 2004; Awad and Khanna, 2015).

Manual selection of Ct and g costs a long time and thus it is not realistic to obtain the best parameters. Therefore, in this study,optimal algorithms are utilized to opt and optimize these two crucial parameters. Currently, many meta-heuristic algorithms,such as moth flame optimization (Zhou et al., 2021d), equilibrium optimizer algorithm (Yu et al., 2021a), Bayesian optimization algorithm(Zhou et al.,2021e)and whale optimization algorithm(Qiu et al.,2021),have been widely used to obtain satisfactory outcomes.In the current study, five types of optimization algorithms, i.e. GS,GA, PSO, SSA and GWO, are chosen to conduct parametric optimization. They are combined with SVR techniques and used to develop different prediction models.

Fig. 3. Data redistribution by hyperplane in two-dimensional space.

Fig. 4. Data redistribution by hyperplane in three-dimensional space.

3.2. PSO

PSO is a powerful method for solving optimization problems among many evolutionary search methods, prompted by simulating fish schools and bird schools (Eberhart and Kennedy,1995).In PSO, a herd of birds represents a group of particles, and a food source represents a functional purpose. By sharing and transmitting information about the distance between the bird herds and the food source,the location of the food source can be determined by bird schools.This cooperation allows the entire herd of birds to select the best information about the location of the food source and eventually gather around the food. By implementing these steps, the best food source can be found. In PSO, the particle numbers are initialized by PSO and each particle possesses the same probability of being selected as a candidate solution for the defined problem. In the next step, two crucial properties of each particle need to be determined carefully,i.e.updated speed(V)and iterative position (X) (Poli et al., 2007). The goodness of each particle is appraised by the fitness function and the positions of particle swarms are updated based on the evaluation results of the fitness function. By iterations, particle swarms reach the best position according to the pre-defined target function by users.

Some important parameters in PSO are updated as follows, by which the new position and velocity can be determined:

Fig. 5. A general process of PSO.

Fig. 6. A general process of GA.

where t represents the current iteration number; rand(A) and rand(B)signify the random number in the range of(0,1);Pbestand Gbestdenote the best position of the single-particle and the whole particle swarm, respectively; c1and c2are the constants which control the acceleration of particles (Zhou et al., 2013); and ω denotes an inertia weight which determines the balance between global optimization and local optimization(Shi and Eberhart,1998;Poli et al., 2007). Generally, the value of ω decreases with the iteration.It can be determined as follows:

where ωmaxis the maximum inertia weight,ωminis the minimum inertia weight, and Iterationmaxis the maximum iteration.

In SVR, the objective of the PSO is to optimize two important hyper-parameters Ct and g.By updating the velocity and position of the whole particle swarm, the particle swarm will eventually achieve global optimization. The general process of PSO optimization can be interpreted in Fig. 5.

3.3. GA

GA was first proposed and developed by Holland (1992).Inspired by Darwin’s theory of evolution (Zhou et al., 2021b), GA reckons that the strongest individual has more opportunities to reproduce,while undesirable individuals cannot obtain the right to spawn.

In GA,each chromosome is encoded as a candidate solution for a re-defined problem. GA can conduct global optimizations and by implementing the adaptive search process, a better optimal solution replaces the former one, as like other meta-heuristic algorithms. Generally, there are three steps in GA, i.e. chromosome initialization, ranking and selecting, and mating and mutation(David,1991).

Fig. 7. A general process of SSA.

In the beginning, GA randomly produces some chromosomes(individuals)and these individuals constitute a swarm.This step is so-called chromosome initialization. Next, the re-defined function judges the fitness of these individuals and gives corresponding scores.The individual with higher fitness will have the qualification to reproduce, while undesirable individuals will be eliminated. At the mating and mutation stage, those desirable chromosomes will change their genes to generate new chromosomes.Sometimes the mutation will occur to bring some fresh genes to join the reproduction process.The mutation probably is good for reproduction or may not, but it can prevent local convergence. By iteration, the evolutionary process tends to move to the best optimal solution.General working architecture of GA is shown in Fig. 6.

3.4. SSA

Like the above two mentioned algorithms, SSA belongs to a group of swarm intelligence-based optimization algorithms(Mirjalili et al., 2017; Faris et al., 2020). It was inspired by the hunting and navigation behaviors of salp swarms.A salp is a kind of marine mollusk with a similar appearance to a jellyfish.They move and hunt together by connecting like social animals. In the salp chains,there is a leader salp and a host of follower salps.The leader moves towards the direction of the food source and leads the movement of followers. This behavior can be regarded as global optimization and the followers explore the food in the local space.This behavior can be regarded as a local optimization. With the cooperation of leader salp and follower salps, falling into local optimum can be largely reduced.

In SSA, there are mainly three steps to optimize the hyperparameters, i.e. salp swarm initialization, update of leader, and update of followers. In the initialization process, the number and position of salp swarm will be defined. In addition, the best food source (the fitness function) will be determined according to the user requirement.In the next step, the optimization process starts in the search domain and this search domain can be represented by a matrix named Si:

In this search space,the leader salp will go into action at first so that it can provide the guide for follower salps.The update equation of the position of leader salp is written as

where gbjand hbjdenote the searching upper and lower boundaries in the j-dimension,respectively;fjrepresents the location of food in the j-dimension; and RBand RCindicate two random parameters with the value of[0,1].The parameter RBdetermines the length of salp movement,and it can control the direction of movement.The parameter RCplays a significant role in controlling the exploration and exploitation. It can be explained as follows:

Fig. 8. A diagrammatic sketch of GWO.

Fig. 9. A general process of GS optimization.

Fig.10. A schematic diagram of 5-fold cross-validation.

where kgrepresents the current iteration and Kgrepresents the total iteration number.From Eq.(13),it can be found that with the increase in iteration,parameter A will decrease.

For salp followers, they move under the guidance of a salp leader. The new location of salp followers obeys the following equation:

The aforementioned steps will be operated recursively before the optimization reaches the stopping criterion. The optimal process of SSA is shown in Fig. 7.

3.5. GWO

GWO is a kind of novel swarm-based algorithm introduced by Mirjalili et al. (2014). GWO forms a strict hierarchy structure inspired by the hunting and social behavior of wolves, where different wolves play different roles in a grey wolf group.Every wolf is dominated by the harsh social order. Basically, there are four types of wolves involving:

(1) The alpha wolf(α)may not be the strongest wolf but it plays the manager role in a whole wolf herd. It is responsible for commanding and all other wolves would submit the order from the alpha wolf.

(2) The beta wolves (β) mainly assist the alpha wolves to make decisions and provide advice. They make the action around the alpha wolves and the command from the alpha wolf would be transmitted from the beta wolves to the delta wolves and omega wolves.

Table 2 Parameter configurations of four meta-heuristic algorithms.

(3) The delta wolves(δ)would implement more trivial work like guarding,hunting and caretaking.They obey the order of the alpha wolves and beta wolves.

(4) The omega wolves (ω) are the lowest level of a wolf herd.They have to submit to all other wolves. However, they are still important because without them the wolf herd would have internal problems such as cannibalism.

In the implementation of GWO, the wolf herd updates its distribution and location according to the evaluation results of potential prey coordinates.The wolves make corresponding decisions and approach the prey step by step by iterations. Basically, the social behavior of grey wolves can be divided into three stages:social hierarchy, encircling and hunting:

(1) At the first stage, the hierarchy structure would be established. The grey wolves with the best fitness in the swarm would be labelled as “alpha”, “beta”, “delta” in order. While the rest wolves are regarded as “omega”wolves. During the iteration process, evolution would be directed by the aforementioned three wolves.

(2) The second stage is encircling prey. At this stage, wolves approach and encircle prey and this process can be formulated as

Fig.11. Optimization performance with different swarm sizes: (a) PSO-e-SVR and (b) PSO-v-SVR.

(3) Hunting stage. We assume that in each iteration, three wolves with the best fitness would be retained, i.e. “alpha”,“beta”and“delta”wolves.These three wolves have the best ability to recognize the position of prey and the other wolves would update their positions.This behavior can be shown as follows:

The best prey position can be procured by calculating the average of the prey position for “alpha”, “beta” and“delta”wolves as follows:

For global optimization, the “alpha”, “beta” and “delta” wolves collect the food information respectively and transmit the foodinformation to the whole wolf herd.By employing|A|>1,the wolf individual can conduct a global search. As aforementioned, C is a random vector with the value of [0, 2] which can avoid local optimization.A general diagrammatic sketch of GWO is shown in Fig.8.

Table 3 Performance and scores of different swarm sizes of PSO-e-SVR.

3.6. GS optimization

The GS is a kind of classical parametric optimization method and has been proved to be effective to optimize the model(Zhou et al.,2012). The main idea of this method is to test all combinations of the given parameters and find the most suitable one according to the fitness function. While implementing the GS, the search range and step need to be determined. A general process of GS optimization is demonstrated in Fig. 9.

4. Principal component analysis and cross-validation

A principal component analysis(PCA)is used for extracting and recognizing the significant information from a group of multidimensional variables (Cadima and Jolliffe, 1995; Croux and Haesbroeck, 2000; Higuchi and Eguchi, 2004). By matrix transformation, original data are mapped into a new dimension. The data with large dimensions are presented with fewer dimensions.However, data information is discarded by such operations. The original data are compressed and good for calculation. Byimplementing the PCA, 19 variables are transmitted into 7 new variables and the overfitting problems can be improved.

Table 4 Performance and scores of different swarm sizes of PSO-v-SVR.

Fig.12. Optimization performance with different swarm sizes: (a) GA-e-SVR and (b) GA-v-SVR.

Table 5 Performance and scores of different swarm sizes of GA-e-SVR.

Cross-validation is a useful method for assessing the model robustness and generalization.It can avoid over-fitting and underfitting to some extent. In this study, the k-fold cross-validation method is utilized (Kuhn and Johnson, 2013). Original training set is equally separated into k sets.k-1 sets are treated as new training sets and the remaining set is used for validation. This process is operated for k times.The average accuracy of these k models is used as the performance index of the regression model.As suggested by Kohavi (1995),10-fold or 5-fold cross-validation is recommended.In this study, 5-fold cross-validation is utilized according to the scale of the datasets.A general display of 5-fold cross-validation is shown in Fig.10. In this figure, P1, P2, P3, P4 and P5 represent the prediction results of the corresponding fold, respectively.

5. Evaluation metrics

To control and compare the performance of proposed algorithms, three mathematical evaluation metrics are adopted(Hasanipanah et al., 2016; Mehrdanesh et al., 2018; Zhang et al.,2020a, 2021a; Li et al., 2021b), i.e. R2(Eq. (23)), MSE (Eq. (24))and VAF (Eq. (25)). Generally, R2and VAF values equal to 100 and MSE values equal to zero indicate the best prediction performance.For different optimization algorithms, they would output different prediction results and thus these values can be used for determining the best optimization method. Meanwhile, for thesynthetical assessment performance of these five algorithms, a comprehensive evaluation system is employed(Zorlu et al.,2008).In this system,the best prediction performance is endowed with a higher score and inferior performance is given with a lower score.By adding all scores,each model obtains cumulative scores and the one with the highest scores is considered as the best.

Table 6 Performance and scores of different swarm sizes of GA-v-SVR.

6. Discussion

Fig.13. Optimization performance with different swarm sizes: (a) SSA-e-SVR and (b) SSA-v-SVR.

Table 7 Performance and scores of different swarm sizes of SSA-e-SVR.

Table 8 Performance and scores of different swarm sizes of SSA-v-SVR.

To develop e-SVM and v-SVM-based MFS prediction models,four types of meta-heuristic algorithms and GS method are tested and ten models are established and compared. The aforementioned 19 influential factors are utilized as the input parameters and the MFS is used as the output parameter.In metaheuristic algorithms, many parameters influence the optimal effect.Among these parameters,the swarm size and the number of iterations have significant impacts on optimization performance(Koopialipoor et al., 2019; Li et al., 2021a; Yu et al., 2021c; Zhou et al., 2021f). Therefore, in this study, these two parameters are discussed and compared. Detailed parameter configurations are listed in Table 2.

6.1. PSO-SVR optimization

To achieve the best performance of blasting MFS prediction,several parameters including velocity equation, number of iterations and swarm size,which have a significant influence on PSO,are carefully selected. Among these parameters, the number of iterations and swarm size bring a significant impact on optimization performance. Generally, with the increase of iteration, the optimization performance tends towards stability but the calculation time also increases.For judging the optimization performance,the MSE value is used in this section.

During testing, when the number of iterations is 500, the prediction performance is stable.Therefore,the number of iterations is determined to be 500 and the swarm sizes equal to 60, 70, 80, 90 and 100 are tested to choose the best optimal parameters.As Fig.11 depicts, the optimal process of each swarm size is different and with the increase of iteration, MSE value decreases. Each model obtains the lowest MSE value, when it reaches the end of optimization.

To compare the prediction performance of different swarm sizes, three evaluation performance metrics are referenced in this study and each performance is given with a corresponding score based on the comprehensive evaluation system(Zorlu et al.,2008).In this grading system, better performance is assigned a higher score and by adding all scores, each parameter configuration obtains comprehensive scores. The one with the highest score is regarded as the best parameter combination.

From Table 3,when the swarm size is equal to 80 for the PSO-e-SVR, the established model obtains the best prediction performance, where R2, MSE and VAF values are 0.8931, 0.00299 and 88.38, respectively for the training sets, and 0.8331, 0.01413 and 81.13,respectively for the testing sets.In addition,when the swarm size is 80, the evaluation performance is superior to the other models except for the MSE of the testing sets.

For the PSO-v-SVR, when the swarm size is 60, the developed model procures the best prediction performance, where R2, MSE and VAF values are 0.8743,0.00362 and 85.99,respectively for the training sets, and 0.8313, 0.01313 and 82.35, respectively for the testing sets. When the swarm size is 60, the prediction model shows strong fitness abilities in training sets, while for the testing sets, the model with the swarm size equal to 100 seems to have stronger fitness abilities.Detailed results of the training and testing sets can be found in Table 4.Then,the accumulation bar charts are shown in Figs. A1 and A2 in Appendix A to help to interpret the ranking results.

Fig.14. Optimization performance with different swarm sizes: (a) GWO-e-SVR and (b) GWO-v-SVR.

Table 9 Performance and scores of different swarm sizes of GWO-e-SVR.

Table 10 Performance and scores of different swarm sizes of GWO-v-SVR.

6.2. GA-SVR optimization

As mentioned earlier,GA can bring a positive effect on parameter selection and optimization. Before implementation, the most effective GA parameters that are utilized to govern GA-SVR prediction models should be opted. Similar to the development of PSO-SVR models, the impacts of the number of iterations and swarm sizeshould be discussed.For the sake of comparison,the same number of iterations and swarm size are employed in GA-SVR.The optimization process of GA-e-SVR and GA-v-SVR can be found in Fig. 12. By calculating the performance scores,it can be observed that when the swarm size is equal to 100 for GA-e-SVR, the proposed model can obtain the best prediction performance. At this time, R2, MSE and VAF values are 0.8474, 0.00496 and 81.05, respectively for training sets,and 0.8174,0.01488 and 77.97,respectively for testing sets.From Table 5,when the swarm size is equal to 100,all evaluation metrics obtain the best performance. Therefore, for the GA-e-SVR, it can be concluded that the model with a swarm size equal to 100 has the best robustness in this study.

Table 11 Performance and scores of different grid steps of GS-e-SVR.

Table 12 Performance and scores of different grid steps of GS-v-SVR.

Fig.15. GS-e-SVR optimization curve with grid step equal to 0.4.

Fig.16. GS-v-SVR optimization curve with grid step equal to 0.8.

For the GA-v-SVR, when the swarm size is 80, the procured model brings the best comprehensive prediction results, with R2,MSE and VAF values of 0.8347,0.005 and 80.84,respectively for the training sets, and 0.8368, 0.0125 and 82.44, respectively for the testing sets. From Table 6, although the model with swarm size equal to 80 does not obtain prominent performance in each single evaluation metric, its comprehensive performance is better than the other models.For the model with the swarm size equal to 70,it produces a good prediction effect in training sets but fails to bring desirable feedback in the testing sets. As for models with swarm sizes equal to 90 and 100,they obtain superior scores in R2and MSE in testing sets but show mediocre performance in the other metrics.The intuitive scoring results of GA-e-SVR and GA-v-SVR are shown in Figs. A3 and A4 in Appendix A, respectively.

6.3. SSA-SVR optimization

As a kind of novel swarm intelligence algorithm, SSA shows its superiority and feasibility in various optimal problems.SSA is easy to implement and adjust. For the sake of comparison, the same hyper-parameter values are also employed in the SSA-based models. The optimization process of SSA-based models with different swarm sizes is shown in Fig.13.

By comparing the model performance, it can be observed that when the swarm size is equal to 80 for SSA-e-SVR, the proposed model can obtain the best prediction performance according to Table 7. For the SSA-e-SVR, R2, MSE and VAF values are 0.8463,0.00139 and 80.9, respectively for the training sets, and 0.8157,0.00415 and 77.88, respectively for the testing sets. From Table 7,when the swarm size is equal to 80,the prediction model produces excellent prediction performance and each evaluation metric generates the best feedback.

From Table 8,when the swarm size is 100,the proposed model produces the best synthetical prediction results for the SSA-v-SVR with R2, MSE and VAF values of 0.8345, 0.00139 and 80.81,respectively for the training sets, and 0.8369, 0.00347 and 82.44,respectively for the testing sets.Although the model with a swarm size of 70 generates outstanding performance in training sets, it does not work well for the testing sets.Therefore,its generalization and robustness are poor under this situation. Finally, the intuitive ranking results of SSA-e-SVR and SSA-v-SVR are depicted in Figs.A5 and A6 in Appendix A, respectively.

6.4. GWO-SVR optimization

The GWO is easy to interpret as a kind of effective optimization strategy. In the GWO, three kinds of wolves control theoptimization process.Similarly,the number of iterations is set to be 500 and the swarm(wolf)sizes are set to be 60,70,80,90 and 100,respectively.From Fig.14,it can be seen that GWO optimization has a faster convergence velocity, which reflects its eminent optimization abilities.By comparing the model performance in Table 9,it can be observed that when the swarm size is equal to 90 for GWOe-SVR, the proposed model can obtain the best prediction performance. For the GWO-e-SVR, R2, MSE and VAF values are 0.8463,0.00139 and 80.89, respectively for the training sets, and 0.8156,0.00414 and 77.9,respectively for the testing sets.Although the R2value in testing sets is a little bit inferior to the one with the swarm size of 100, the model with the swarm size of 90 has the higher comprehensive score in any other evaluation metrics.

Table 13 Comparison of optimal algorithms of e-SVR.

Table 14 Comparison of meta-heuristic algorithms of v-SVR.

Fig.17. Predicted results using e-SVR-based models for (a) testing set and (b)training set.

For the GWO-v-SVR, the proposed model produces the best comprehensive prediction feedback when the swarm size is 100,as shown in Table 10,with R2,MSE and VAF values of 0.8355,0.00138 and 80.98, respectively for the training sets, and 0.8353, 0.00348 and 82.41, respectively for the testing sets. It is worth noting that although the model with swarm size equal to 100 has a superior performance in training sets,it fails to procure similar performance in testing sets. That is to say, the generalization of this model is worthwhile to be improved and discussed in the future.Perhaps,by changing and testing different swarm sizes and numbers of iterations,it produces more persuasive scenarios.Similarly,the ranking results of GWO-e-SVR and GWO-v-SVR are displayed in Figs.A7 and A8 in Appendix A, respectively.

Fig.18. Predicted results using v-SVR-based models for(a)testing set and(b)training set.

6.5. GS-SVR optimization

In the GS optimization, there are mainly two parameters that need to be tuned, i.e. grid step and grid bound. In this study,different grid bounds are changed and tested.It is found that grid bound almost does not have any influence on the optimization effect.Therefore,the search bound of Ct and g is set to(2-8,28).In the next step,five grid steps are tested,i.e.0.2,0.4,0.6,0.8 and 1 in the GS-e-SVR and GS-v-SVR models.From Table 11,the same grid steps (equal to 0.2, 0.4 and 0.6) can produce the same prediction performance. From Tables 11 and 12, it can also be said that finer grids always do not bring better prediction performance.It is due to the limitation of the GS method. It is almost impossible to examine the performance of all nodes in specified grid bound,unless the grid step is set to be ultra-small.Therefore,some more effective parameter combinations perhaps are missed. As per Figs.15 and 16, one grid node represents one kind of parameter combination and the longitudinal axis represents the corresponding MSE value. With the grid color changes from yellow to dark blue, a lower MSE value is obtained. When the grid step equals 0.2, 0.4 and 0.6, the GS-e-SVR method would obtain the best comprehensive performance, and for the GS-v-SVR method,the most suitable grid step is 0.8.

Fig.19. Sensitivity analysis of different factors on blasting MFS.

6.6. Comparison of optimal algorithms

Based on the previous discussions, the model with the highest comprehensive scores of each optimal algorithm is selected and compared. The detailed data are listed in Table 13. Among the e-SVR-based MFS prediction models, PSO-e-SVR obtains the best comprehensive scores.In addition,PSO-e-SVR also obtains the best performance in R2and VAF. Only considering the performance metric MSE,GWO-e-SVR outperforms the other algorithms for both training and testing sets. SSA-e-SVR has prominent performance with the metric MSE when being tested in training sets.

For the v-SVR-based MFS prediction, GWO-v-SVR procures the best comprehensive scores in Table 14.However,SSA-v-SVR is also competent because it has an excellent performance in the testing sets which means this algorithm has strong generalization ability based on the current data.Apart from that,the ability of PSO-v-SVR is worthwhile to be declared because it obtains the best scores for R2and VAF in the training set.

Finally, e-SVR and v-SVR-based models with the highest synthetical scores are screened and compared. Given the aforementioned discussions,it can be observed that PSO-e-SVR and GWO-v-SVR obtain the best scores and GWO-v-SVR shows better performance in four metrics,i.e.MSE of the training set,R2,MSE and VAF of testing set. While PSO-e-SVR only outperforms in two metrics,i.e. R2and VAF of the training set. The GWO-v-SVR model shows superior performance in the testing set which proves that it has better generalization and robustness capabilities. For AI-based models, subtle advantages can be magnified for large-scale data.Therefore, GWO-v-SVR can be regarded as the best method for blasting MFS prediction in this study. Corresponding swarm size and number of iterations are 100 and 500, respectively. Finally,predicted MFS values along with their actual values are demonstrated in Fig. 17 for e-SVR and v-SVR-based models,respectively, and the results are demonstrated in Fig.18.

7. Sensitivity analysis

For identifying and comparing the sensitivity of different influential factors on blasting MFS, the cosine amplitude method proposed by Yang and Zhang (1997) is utilized. Each input and output are considered as a column matrix.Then total of 20 column matrices are obtained as follows:

where the length of each matrix is equal to the number of datasets and then the sensitivity of different influential factors on blasting MFS can be analyzed by(Yang and Zhang,1997):

where xanrepresents one kind of input variable;xbnrepresents the output variable;and n represents the sequence number of the data,and for each variable, total of 76 groups of data are evaluated.

By implementing the cosine amplitude method, the relative strength of effects (RSE) of each variable can be evaluated. It is found that the most sensitive factor is UCS for blasting MFS, as shown in Fig.19.This result is reasonable because when the violent blasting shock wave occurs,the dynamic stress on the rock will also increase sharply and the strength of the rock itself counteracts this effect.Therefore,this indicates that UCS has a positive influence on the fragment size.The significance of UCS on the fragment size has also been reported in some previous researches(Bond,1952; Lilly,1986).

The sensitivity of each factor not less than 0.8 indicates that selected factors indeed contribute to the blasting MFS. Finally, the sensitivity of different parameters on blasting MFS based on available data can be sorted in ascending order as:Qe,J,J.B,L.Wd,H,L,ST,H.B,NH,De,ST.B,B,S,PF,B.D,Wd,D,S.B and UCS.It is noted that different data would bring different results. However, current analysis results could provide some references for blasting designs with similar conditions.

8. Limitation

By utilizing SVR as the predominant strategy to predict the blasting MFS, satisfactory prediction accuracies are procured.However,there are still some drawbacks and limitations that need to be improved in future work. Firstly, the scale of data used to establish the evaluation models is still small and only a several dozens of samples are collected. It tends to be easy to mine more useful information from larger datasets and thus the performance of prediction models remains to be strengthened.Secondly,deeply analyzing blasting mechanisms is significant for mining more related factors of blasting MFS. Thirdly, more advanced metaheuristic algorithms (Yu et al., 2020, 2021b) are worthwhile to be combined with SVR prediction models to improve the prediction accuracy. Apart from that, the potential of some other powerful regression tools such as random forest(Zhou et al.,2019a),extreme gradient boosting (Ding et al., 2020; Zhang et al., 2021b), and gradient boost machine(Zhou et al.,2019b;Zhang et al.,2020a)are not investigated and compared in this study.

9. Conclusions

The blasting MFS prediction and evaluation are always influenced by many factors and the relationship between these factors and fragment size is elusive.Therefore,it is hard to be evaluated by means of a certain function. AI-based techniques can simulate sophisticated relationships between influential factors and output targets compared to the conventional functions. In this study, a group of blasting fragment datasets which contain 19 types of influential factors are adopted.To control and estimate the blasting MFS, two types of SVR-based techniques, i.e. v-SVR and e-SVR,are employed. Then, five types of optimal algorithms are combined with SVR to optimize the hyper-parameters. Three types of mathematical indices combined with a comprehensive ranking system are utilized to evaluate the model performance. As a result, it is found that GWO-v-SVR obtains the most comprehensive prediction performance with R2, MSE and VAF values of 0.8355, 0.00138 and 80.98, respectively for the training set, and 0.8353, 0.00348 and 82.41, respectively for the testing set.

According to the sensitivity analysis results, the UCS plays the most important role in influencing the blasting MFS. Compared with other AI-based studies, this research takes more influential factors into consideration. On the one hand, the prediction results can provide significant guidance for complicated blasting design.On the other hand, the feasibility of AI-based techniques is validated again in open-pit mining fields. Finally, this study offers a more accurate prediction approach for multi-parameter coupling blasting MFS evaluation compared with the original paper(Sharma and Rai, 2017).

Declaration of competing interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication, and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

This research was funded by the National Natural Science Foundation of China (Grant No. 42177164) and the Innovation-Driven Project of Central South University(Grant No.2020CX040).The first author is supported by China Scholarship Council (Grant No.202006370006).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.07.013.