1 Introduction

As a common solution to eliminate the rock mass, blasting operations are used in some engineering works such as tunnel excavation, road construction, and hydraulic channels [1]. Most of the explosive operations have a lot of energy that can have impacts on the environment and surrounding areas [2,3,4,5,6]. The common environmental issues of blasting are flyrock, air overpressure, back-break and ground vibration [7,8,9,10,11]. Flyrock can cause the most important effects of damages among them according to several scholars [12].

In flyrock, the parameters of charge confinement, mechanical strength of the rock mass, explosive energy have an important relationships with each other [13]. Based on some researches, any mistakes in designing these parameters will result in flyrock [13, 14]. When flyrock phenomena has happened, a lot of fragmented rocks will be created and fly a distance from the blast face [15].

Three main categories of flyrock are included: cratering, rifling and face bursting. Cratering will occur because of the too small ratio of stemming length to diameter in blasting face. Rifling will happen when stemming material is incompetent or is insignificant. In the third case, which is named face bursting, flyrock may occur due to the production of high-pressure gases in weak rocky plates. Therefore, the explosion near the weak stone plates causes the face bursting state.

According to previous researches, controlled and uncontrolled factors can affect flyrock. The most important controllable factors are incompetent stemming, inappropriate burden and spacing, inaccurate drilling, too much explosive energy, inadequate delay timing and unwarranted powder factor [1, 16, 17]. In the case of uncontrollable factors, the most effective factors are related to the rock mass properties.

Several empirical relationships are presented for the prediction of flyrock from the blasting face [15, 18, 19]. These relationships used one or two influential factors, which lead to receive a low-performance prediction of this method. In addition, to increase the safety of the surrounding area, flyrock phenomena must be predicted with higher accuracy level before blasting [17, 20]. Therefore, to find the flyrock distance, more researches have to be done and models with higher performance prediction have to be presented.

The previous computational techniques developed to predict flyrock distance comprising of artificial neural network (ANN), fuzzy inference system (FIS), Monte Carlo simulation, multiple regression analyses, support vector machine (SVM), rock engineering systems (RES), and adaptive neuro-fuzzy inference system (ANFIS) models [14, 21,22,23,24]. According to the above-mentioned methods, providing new ways to predict this phenomenon is necessary.

In engineering sciences, the use of ANNs (as a branch of artificial intelligence) has been highlighted by many investigators [25,26,27,28,29,30,31]. Such networks are good tools for forecasting issues, however, they have several limitations such as low learning speed and falling into local minima [32,33,34]. As mentioned in literatures [20, 23, 35,36,37,38,39], using efficient optimization algorithms (OAs), these limitations can be overcome. Various OAs such as particle swarm optimization (PSO), imperialism competitive algorithm (ICA) and genetic algorithm (GA) can be applied to solve continuous and discontinuous problems. Based on powerful ability of global search of these OAs, weights and biases of an ANN network can be determined to improve its performance prediction. The mentioned hybrid models have been widely-utilized to solve non-linear and complicated engineering problems.

In this research, to predict flyrock phenomenon, three hybrid intelligent techniques, namely ICA–ANN, GA–ANN and PSO–ANN, are applied. These models were proposed based on the most important parameters influencing flyrock. In the following, after introducing flyrock prediction models and the applied models in this study, some explanations regarding the established database will be given. Then, modelling procedures of the applied techniques are described and finally, the best predictive model will be selected to predict flyrock.

2 Flyrock empirical methods

In recent years, many empirical studies have been delivered to predict flyrock in mining engineering. An empirical model was found by Lundborg et al. [19] according to two parameters as follows:

$${\text{Flyrock}}\,=\,260\,~ \times \,~{D^{2/3}},$$
(1)
$${T_{\text{b}}}\,=\,0.1\, \times \,{D^{2/3}},$$
(2)

where D is the hole diameter in inches, and Tb is the size of the fragmented rock in metre.

Raina et al. [40] performed a research according to selected parameters of rock mass and blast design for evaluating the horizontal (FSH) and vertical (FSV) safety factors of flyrock. In another research, the parameters including density, explosive density hole diameter, and confinement state were used by McKenzie [41] for prediction of flyrock and particle (rock) size. A new empirical equation to assess flyrock was used by Trivedi et al. [42]. For the development of this model, 95 explosive data sets were used and an equation based on specific charge, charge concentration, rock strength, burden, stemming length and rock quality designation was proposed to estimate flyrock.

Two power empirical equations were introduced in the study carried out by Marto et al. [23] who developed two high-performance empirical formulations for prediction of flyrock. These results were obtained from 113 operations where each of them contained charge per delay and powder factor. Furthermore, Jahed Armaghani et al. [14] presented an empirical technique for predicting flyrock. This technique was based on graph shown for different values of maximum charge per delay in a range of (75–550 kg) and also for various powder factor values in a range of (0.5–1.1 kg/m3).

3 Intelligent techniques

3.1 Artificial neural networks

Due to a structure of the human brain, artificial neural network (ANN) [43] can be created and developed to process information. The ANN structure consists of three main parts: input, hidden and output layers. In each network layer, the binding elements, namely neurons, transmit data from one layer to the next. Due to the strengthening or weakening of this transfer, the weights in each network control this transfer. To calculate the output of each layer’s neuron, an activation function of linear or sigmoid should be used. The number of used neurons in each layer is specified by the total number of inputs. Usually, the number of neurons can be obtained using a complex way or trial and error procedure. Among available network training algorithms, the back propagation (BP) training algorithm is more common in engineering sciences [44,45,46]. Briefly, two main parts of the ANN modelling are creation of a network construction and the relevant weights determination. Based on minimization error values, the network weights are adjusted by BP training algorithm. The values obtained at each stage are compared with the desired output values. If the errors are not desirable, the process should be continued to get desired values and reduce the system error [46,47,48].

3.2 Genetic algorithm

Genetic algorithm (GA), which was first developed by Holland [49], is one of the well-established optimization methods. This OA is inspired by the theory of natural selection. This method were expanded by Goldberg [50]. GA has been widely-performed to optimize various problems in engineering and science. One of the main advantages of this algorithm is its ability to solve complex and highly non-linear problems. For optimization purposes, such as linear or non-linear, static or dynamic (change with time), continuous or discontinuous or contain a random noise, GA can be used to solve. In addition, GA is considered as a problematic algorithm due to its limitations such as determining various parameters of the algorithm (population size and genetic operator rates) and creating the proper function. To determine these values, the designer should be very careful while they will affect the convergence of the algorithm and also its results [51, 52]. In GA, chromosomes have a fixed length that encodes issues to linear binary strings between 0 and 1. These chromosomes cause production of generation. As shown in Fig. 1, the chromosome is selected as a random characteristics and based on these characteristics, chromosomes are evaluated. Then, they are selected using genetic operators of the remaining chromosomes and start generating new generations. Crossover chooses between parents and mutation works in a range of 0–1. This process is repeated until creation of the best generations evaluated based on their performance [53, 54].

Fig. 1
figure 1

GA algorithm

Full size image

3.3 Particle swarm optimization

Another OA used in this study is particle swarm optimization (PSO) which was developed by Kennedy and Eberhart [55]. PSO is inspired by cumulative behaviour of particles. Among all advantages of PSO, a high learning speed and using less memory compared to GA should be noted. In PSO, to find the best position, a swarm of particles searches the best personal (pbest) and the best global (gbest) positions [35]. In other hands, in each system iteration, the particle moves toward finding the best positions (pbest and gbest). The velocity and position of particles are obtained as follows:

$${V_{{\text{new}}}}\,=\,w\, \times \,V\,+\,{C_1}\, \cdot \,{r_1}\left( {{p_{{\text{best}}}}\, - \,X} \right)\,+\,{C_2}.{r_2}({g_{{\text{best}}}}\, - \,X),$$
(3)
$${X_{{\text{new}}}}\,=\,X+{V_{{\text{new}}}},$$
(4)

where V and X denote current velocity and position of particle, respectively, C1 and C2 are two positive acceleration constants, Vnew and Xnew denote new velocity and position of particle, respectively, w denotes the inertial weight, and r1 and r2 represent the random numbers in (0, 1). More information about the PSO algorithm description/implementation can be found in the other researches [45, 56]. Furthermore, Fig. 2 shows details of a PSO algorithm.

Fig. 2
figure 2

PSO algorithm

Full size image

3.4 Imperialist competitive algorithm

Imperialist competitive algorithm (ICA) was first introduced and developed by Atashpaz-Gargari and Lucas [57] as a global search population-based on optimization problems. The beginning of ICA is with a production of randomly initial population called countries. The process continues to generate N number of countries (Ncountry), and then the number of imperialists, i.e. Nimp should be selected as a specific number of the countries which have the lowest costs. The remaining countries (Ncol) are used as special functions (functions of the imperialist normalized costs) among other empires. In this algorithm, imperialists are more powerful when they have more colonies. Three operators including assimilation, revolution and competition are the main parts of ICA [58, 59]. The part of ICA body is related to colonies that are equally absorbed by the imperialists. However, the revolution is causing many sudden changes. In the competition part, the imperialists are struggling to get more colonies, and in this competition, any empire that can achieve the desired criteria eventually wins. This process is repeated until the end of the desired benchmark. The number of decades in ICA has similar process of the number of generations in GA and the number of particles in PSO. To design them, evaluating in results of root mean square (RMSE) can be useful. More information/facts about ICA are available in several researches [33, 57, 60]. Figure 3 shows a structure of ICA algorithm.

Fig. 3
figure 3

ICA algorithm

Full size image

3.5 Hybrid algorithms

In engineering applications, many research have been carried out to enhance the ability of ANN models through OAs such as GA, PSO and ICA (e.g. [51, 61,62,63,64,65]). Due to the weakness of BP in finding the accurate global minimum, the ANN model may achieve undesirable results [66]. Nevertheless, the ANN model is more likely to be caught up in local minima, while OAs by setting weights and biases of ANN could solve the mentioned ANN problem. In this study, three methods of hybrid systems, i.e. ICA–ANN, GA–ANN and PSO–ANN are constructed to predict SF of the slopes under static and dynamic conditions. In the created systems, ICA, GA and PSO search for global minimum, and then ANN selected it for achieving the best system results.

4 Studied quarry sites and data collection

Data were collected to predict flyrock from six granite quarry mines in Johor state, Malaysia. Their names are Taman Bestari with latitude of 1°60′41″N and longitude of 103°78′32″E, Ulu Tiram with latitude of 1°36′41″N and longitude of 103°49′20″E, Trans Crete with latitude of 1°31ʹ21ʺN and longitude of 103°52′60″E, Putri Wangsa with latitude of 1°35ʹ32ʺN and longitude of 103°48′4″E, Ulu Choh with latitude of 1°31ʹ48ʺN and longitude of 103°32′41″E and Masai with latitude of 1°29ʹ42ʺN and longitude of 103°52′28″E. The target of explosions is to produce 8000–24,000 tons per month. The number of explosive operations in these sites is varied from 6 to 15. In the studied sections, the rock quality designation (RQD) and the values recorded by Schmidt hammer are 25–55% and 17–42, respectively. Figure 4 shows a view of Ulu Tiram-studied quarry.

Fig. 4
figure 4

A view of Ulu Tiram-studied quarry

Full size image

In the mentioned sites, flyrock phenomenon was considered as one of the most important environmental issues. 262 sets of blasting data were collected where, in each of them, data are included: burden to spacing ratio, blast-hole diameter, powder factor, stemming length, the maximum charge per delay, and blast-hole depth as inputs and flyrock distance as output. For blast operations, ammonium nitrate and fuel oil (ANFO) was used as explosive. In these operations, blast-hole diameters of 75, 89, 115 and 150 mm were utilized. The values of powder factor and stemming length have ranges of 0.44–1.14 kg/m3 and 1.4 and 4.5 m, respectively.

For recording the maximum flyrock, two video cameras were used. To measure the distance of flyrocks, blasting benches were coloured, and using the mentioned cameras, the flyrocks could be seen separately after the operations. Then, the maximum horizontal distance of fragments was considered as the maximum flyrock distance. It should be noted that the data used in this study have been previously utilized by Shirani et al. [67]. Methods of data collection were used similar to previous researchers [2, 5, 15, 33, 40]. So, more information regarding parameters used in this study can be found in the mentioned study.

5 Model development

In this section, descriptions of implementing hybrid models, namely ICA–ANN, PSO–ANN and GA–ANN, in predicting flyrock distance are presented. Effective parameters on ICA, PSO and GA are determined and used to receive higher accuracy level for flyrock prediction.

5.1 ICA–ANN

To obtain the best ICA–ANN model, its important factors/parameters should be investigated. Prior to investigation of ICA parameters, ANN architecture should be determined. This was accomplished by considering a trial and error process and it was found that an architecture of 6 \(\times\) 9\(\times\) 1 (or a model with nine hidden neurons) receives better results. Therefore, the mentioned architecture was used for all hybrid intelligent systems in this study. As mentioned earlier, Ncountry, Ndecade, and Nimp are considered as the most influential parameters on ICA. To determine Nimp, many models were designed using different values of Nimp, i.e., 5, 10, 15, 20, 25 and 30. In these models, Ncountry = 300 and Ndecade = 100 were utilized. Results of this parametric study showed that Nimp = 5 can obtain higher performance system capacity. To select the best value for Ndecade, as displayed in Fig. 5, various models with Ncountry values of 50, 100, 150, 200, 250, 300, 350 and 400 were constructed and evaluated based on their RMSE. As a result, RMSE results are not changed after Ndecade equal to 500. In the last step of modelling, using Nimp = 5 and Ndecade = 500, a various number of countries were considered and their ICA–ANN models were built. These models were evaluated based on performance indices (PIs), i.e., coefficient of determination (R2) and RMSE values as presented in Table 1. A ranking method introduced by Zorlu et al. [68] was used to choose the best hybrid models in this study.

Fig. 5
figure 5

ICA–ANN models with various Ncountry values

Full size image
Table 1 Different Ncountry effects in estimating flyrock distance
Full size table

A complete version of this technique can be found in Zorlu et al. [68] and according to it, a rank value was assigned for each PI in its group (training and testing). For example, values of 0.923, 0.938, 0.945, 0.938, 0.942, 0.943, 0.953 and 0.947 were achieved for R2 of training datasets of models 1–8, respectively, and values of 2, 3, 6, 3, 4, 5, 8 and 7 were assigned for their ranks, respectively. This has accomplished for RMSE results as well. Then, a summation value of rating of R2 train, RMSE train, R2 test and RMSE test was calculated and assigned to each model and, based on them, model 6 (Ncountry = 300) with total rank of 27 is the best ICA–ANN model. Evaluation of the best ICA–ANN result will be given later. It should be noted that all models were built using 80% of whole data as training and 20% of them as testing.

5.2 PSO–ANN

As pointed out before, several parameters such as coefficients of velocity equation, number of particle, number of iteration and inertia weight have a deep impact on PSO algorithm. According to literatures [55, 69], coefficients of velocity equation equal to 2 and inertia weight of 0.25 showed an acceptable results in other implemented PSO works. Therefore, these values were used in all PSO–ANN models. To select the best value for number of iteration, as displayed in Fig. 6, various models with swarm size values of 50, 100, 150, 200, 250, 300, 350 and 400 were built and evaluated based on their RMSE. As a result, RMSE results are not changed after swarm size of 500 for all models. To identify the optimum value for swarm size, a total number of eight PSO–ANN models were constructed to predict flyrock distance as tabulated in Table 2. Similar to previous section, ranking system proposed by Zorlu et al. [68] was performed and based on total rank values (Table 2), model 7 with swarm size of 350 and total rank of 29 shows the best system results. For this model, R2 of 0.958 and 0.959 were obtained for training and testing datasets, respectively. Evaluation of the selected PSO–ANN model will be discussed later.

Fig. 6
figure 6

PSO–ANN models with various swarm sizes

Full size image
Table 2 Different swarm size in estimating flyrock distance
Full size table

5.3 GA–ANN

As stated earlier, to design a GA–ANN model, influence of the effective GA factors should be investigated. Mutation probability values, percentage of recombination were set as 25, and 9%, respectively. As a cross-over operation, a single point with 70% possibility is used. A parametric study was conducted for determination of the maximum number of generation (Gmax) effects on network performance. To obtain the best Gmax, as shown in Fig. 7, a value of 1000 generation was assigned as stopping criteria considering RMSE values. As a result, like two previous hybrid models, after generation No. 500, the network performance is unchanged. Therefore, the optimum generation No. 500 was used in this study. In the final step, a series of hybrid GA–ANN models (Table 3) were created to determine the best population size (among size values of 50, 100, 150, 200, 250, 300, 350 and 400). Results showed that population size of 350 with total rank of 30 can provide higher performance prediction in terms of both R2 and RMSE indices. The obtained results of model number 5 (the best one) will be discussed in detail later.

Fig. 7
figure 7

GA–ANN models with various population sizes

Full size image
Table 3 Different population size in estimating flyrock distance
Full size table

6 Results and discussion

Prediction of flyrock distance due to blasting is the aim of the present study. Hence, the most influential parameters on flyrock were identified and used. Three hybrid intelligent systems, namely ICA–ANN, PSO–ANN and GA–ANN were applied to select the best predictive flyrock model among them. Many hybrid models were constructed for each predictive technique and the best of them was chosen. Results of the selected models of ICA–ANN, PSO–ANN and GA–ANN based on RMSE and R2 indices in predicting flyrock are presented in Table 4. Equations of RMSE and R2 can be found in other studies [70, 71]. High performances of the training datasets prove that the learning process of these predictive models is successful. While, a high accuracy level of testing datasets shows that the developed model is well-generalized. Although all predictive models are capable to predict flyrock, as a result, PSO–ANN predictive model can provide higher performance capacity in terms of R2 values of both training and testing phases. Additionally, RMSE values of (0.052, 0.045 and 0.057) and (0.045, 0.044 and 0.058) were achieved for training and testing of ICA–ANN, PSO–ANN and GA–ANN predictive models, respectively. These results indicated that lower system error can be obtained by developing PSO–ANN model among all implemented models. Predicted flyrock values together with their actual values for ICA–ANN, PSO–ANN and GA–ANN predictive models are displayed in Figs. 8, 9 and 10, respectively. In these figures, predicted results are presented for both training and testing datasets. According to these figures, although all models have acceptable prediction capacity in prediction flyrock distance, PSO–ANN model can introduce as a new hybrid model in this field.

Table 4 Results of the selected hybrid models in predicting flyrock distance
Full size table
Fig. 8
figure 8

Results of ICA–ANN model in estimating flyrock

Full size image
Fig. 9
figure 9

Results of PSO–ANN model in estimating flyrock

Full size image
Fig. 10
figure 10

Results of GA–ANN model in estimating flyrock

Full size image

It is worth mentioning that the same data used in this study have been utilized by Shirani et al. [67]. They developed and introduced a genetic programming (GP) model with R2 values of 0.908 and 0.819 for training and testing datasets, respectively. Comparing their results with the results obtained from this study revealed that all developed hybrid predictive models can apply better than GP model for the same database and could introduce as trustable models in the field of blasting operations.

7 Sensitivity analysis

A sensitivity analysis was performed on the data to determine the impact of each data on the output. Therefore, the method introduced by Yang and Zang in this study was used. All data pairs were utilized to construct a data array X as follows:

$$X=\{ {x_1},{x_2},{x_3}, \ldots {x_i}, \ldots ,{x_n}\} .$$
(5)

Variable xi in array X is a length vector of m as:

$${x_i}=\{ {x_{i1}},{x_{i2}},{x_{i3}}, \ldots {x_{im}}\} .$$
(6)

The strength of the relationship (rij) between datasets Xi and Xj can be expressed as follows:

$${r_{ij}}=\frac{{\mathop \sum \nolimits_{{k=1}}^{m} {x_{ik}}{x_{jk}}}}{{\sqrt {\mathop \sum \nolimits_{{k=1}}^{m} {x_{ik}}^{2}\mathop \sum \nolimits_{{k=1}}^{m} {x_{ik}}^{2}} }}.$$
(7)

Figure 11 shows the relationship between input data and output. As can be seen, among all data, the hole diameter and hole depth have the highest and lowest relation with the output, respectively.

Fig. 11
figure 11

Sensitivity analysis to determine the impact of each data on the output

Full size image

8 Conclusions

In this study, three hybrid models, i.e., ICA–ANN, PSO–ANN and GA–ANN were considered and developed to predict flyrock. To achieve this aim, the most influential parameters on OAs, i.e., ICA, PSO and GA, were identified based on background of these techniques and also available literatures. Then, these parameters were carefully designed using several rounds of parametric studies. At the end of each model designing, eight models were constructed and their related results were achieved based on RMSE and R2. In this step, evaluation of the obtained results was performed using the ranking technique and it was found that all hybrid models can offer a high level of accuracy in estimating flyrock distance. Nevertheless, a slightly higher performance prediction (lower error and higher coefficient of determination) was observed when a PSO–ANN model is developed. RMSE values of (0.052, 0.045 and 0.057) and (0.045, 0.044 and 0.058) were found for training and testing of ICA–ANN, PSO–ANN and GA–ANN predictive models, respectively, which show higher efficiency of the PSO–ANN model compared to other implemented models. Furthermore, in terms of R2, a similar trend was obtained. It can be concluded that if a predictive model with lowest error is needed, a hybrid PSO–ANN model is introduced as a superior one to predict flyrock distance. Additionally, results of sensitivity analysis showed that the effect of hole diameter on flyrock is slightly higher than the effect of others.