Prediction of flyrock throw distance in quarries by variable selection procedures and ANFIS modelling technique

This study suggests application of variable reduction procedures for flyrock prediction. It was aimed to create robust and non-complex predictive models. Eleven operational blast parameters and rock mass properties were measured in an aggregate quarry. Dominant parameters for flyrock occurrence were determined by multivariate statistical methods. Two parallel ANFIS models were developed for flyrock prediction. The first ANFIS model was constructed based on the results of stepwise regression. Burden–hole diameter ratio, in-situ block size and specific charge are the input parameters of ANFIS 1. The second ANFIS model was created based on the results obtained by factor analysis. Burden–hole diameter ratio, bench height–burden ratio, number of holes and charge weight are used as input parameters for ANFIS 2. The calculated mean absolute percentage errors are lower than eight percent for the ANFIS predictions. The median absolute errors are lower than 5 m. The study also investigates alternative accuracy measures to evaluate forecasting performance. Standardized errors, normalized errors and Nash–Sutcliffe Efficiency were found to be useful for model validation. It is concluded that more than a single model can be created for a specific site. Pre-statistical analysis for variable reduction increases performance of the predictive models. Burden appeared to be a significant parameter for flyrock throw.


Introduction
Flyrock is one of the most important environmental effects of bench blasting. Flying rock particles may cause damage, injuries and fatal accidents. It is an important concern especially for quarries that operate near residential areas. On the other hand, flying rock particles may also be dangerous for quarry workers and structures. Several blast design parameters and rock mass properties influence flyrock throw. The first step of flyrock control is the prediction of throw distance using operational blast parameters.
Significant variations in rock mass may cause excessive and unexpected flyrock throw. According to Jimeno et al. (1995), intensely fissured and jointed rocks tend to create more flyrock than massive and homogenous rocks. International Society of Explosive Engineers (2011) indicates that unforeseen geological conditions may cause excessive flyrock. Geological structures such as clay seams, cavities and layers of unconsolidated materials may increase possibility of occurrence of flyrock. Unconsolidated layers, especially in the stemming region, should be observed during drilling operations. Olofsson (1990) and Bhandari (1997) suggest additional stemming for incompetent zones and voids. Charged weak zones may have weak resistance and cause excessive flyrock. Bhandari (1997) also points out large cavities. Large cavities should be passed by additional stemming. Wooden sticks may be used to fill up large cavities. In addition, variability of the strata should be observed. Competent strata, weathered parts and unconsolidated sediments should be recorded. Persson et al. (1994) discuss fissures, joints and weakness planes. They propose that, due to weathering, the surface rocks may be more inhomogeneous than deeper rock mass. A bench face inspection should be performed before starting drilling operations. Kricak et al. (2012) also suggest collecting base information from drill logs to determine incompetent rocks. Kecojevic and Radomsky (2005) focused on fissures, joints, voids and weakness zones. Incompetent zones may cause mismatch between detonation energy and the 1 3 281 Page 2 of 21 resistance of the rock. In addition to drill logging, they suggest stratigraphic modelling and discrete fracture network modelling for extreme conditions. Mishra and Mallick (2013) specifically indicate mud and clay seams. Explosive gases suddenly escape from weak mud and clay zones and cause severe flyrock. Karst features, weak zones with competent rocks, weak layers are defined as unfavorable conditions for flyrock by Raina et al. (2011). Especially weak beds dip into face should be monitored cautiously.
In recent years, researchers have applied different modeling techniques to predict flyrock throw distance. Faradonbeh et al. (2018) used gene expression programming and firefly algorithm to predict flyrock in a granite quarry. Rad et al. (2018) collected flyrock data from the Gole-E-Gohar iron mine in Kerman province of Iran. The researchers considered seven blast parameters and applied least squares support vector machines to forecast flyrock. Koopialipoor et al. (2019) compared genetic algorithms, particle swarm optimization and imperialist competitive algorithm. Burden to spacing ratio, blasthole diameter, powder factor, stemming length, the maximum charge per delay and blasthole depth were used as input parameters. Han et al. (2020) studied in the Malaysian quarries. The researchers applied random forest and Bayesian network techniques. Six different blast parameters were measured in the quarries to create a data set. Lu et al. (2020) used two machine learning models; extreme learning machine (ELM) and outlier robust extreme learning machine (ORELM). Spacing, burden, powder factor, stemming and rock density were considered as input parameters. Nikafshan Rad et al. (2020) applied recurrent fuzzy neural network (RFNN) combined with the genetic algorithm (GA), genetic algorithm-neural network (GA-ANN) and artificial neural network (ANN). Blasting operations were observed in 6-9 m benches in the two different quarries. Coefficient of determination (R 2 ) and root mean square error were used for model validation. Ye et al. (2021) applied two tree-based techniques, genetic programming (GP) and random forest (RF), considering six operational blast parameters. Monte Carlo simulation (MCS) was used with the developed GP equation to analyze the flyrock risk in the aggregate quarries. Guo et al. (2021) studied on deep neural network, artificial neural network and whale optimization algorithm for flyrock throw length prediction. Although different researches have been performed in recent years, a large number of studies were conducted in similar rock formation: Granite quarries in Johol Region, Malaysia (e.g., Armaghani et al. 2020;Han et al. 2020;Murlidhar et al. 2020;Zhou et al. 2020;Guo et al. 2021;Li et al. 2021). The diversity of blast parameters and rock formation is low.
Most of the aforementioned studies directly apply a modelling method to predict flyrock phenomenon. Researchers generally considered 6-8 input parameters for modelling Li et al. 2021;Ye et al. 2021). A variable selection procedure was not applied. Some studies do not consider any rock mass parameters as model input (e.g., Koopialipoor et al. 2019;Nikafshan Rad et al. 2020;Murlidhar et al. 2020). Selection of the input parameters by researchers is a subjective and sometimes misleading methodology. The main novelty of this research is application of multivariate statistics for variable selection. The flyrock phenomenon was discussed from a wide perspective. Eleven different variables, including a rock mass parameter, operational and dimensional parameters were considered. Stepwise regression and factor analysis methods were applied for variable reduction. These methods indicated different dominant variables for flyrock occurrence. By that way, two independent soft computing models were created to predict flyrock throw distance. One of the important finding of the research is that it is possible to create multiple effective models for a specific site. Variable selection procedures have enabled development of alternative predictor models. It is also concluded that burden is a significant parameter that affects flyrock throw distance.
The developed models were tested on the independent test blasts. Most flyrock studies apply insufficient number of error measures for validation, such as root mean square error, mean absolute error, and R 2 (e.g., Nikafshan Rad et al. 2020;Lu et al. 2020;Guo et al. 2021;Li et al. 2021). Model testing is an important phase of the model development. The other novel approach presented in this study is the investigation of alternative accuracy measures for ground vibration prediction. In addition to classical error measures, the standardized and normalized error measures were found to be useful for model validation. Nash-Sutcliffe Efficiency (NSE) (Nash and Sutcliffe 1970) was successfully applied as a goodness-of-fit index.
The general structure of the article can be explained into four sections: sites investigation-methodology, variable selection process, ANFIS modelling and model validation. In the upcoming pages, first, the site measurements are explained. The operational parameters are described and examined in detail. Application of stepwise regression and factor analysis are discussed in the variable selection part. After determination of the input variable, ANFIS modeling stage starts. In section four, first, a brief overview is given for the theoretical background of ANFIS. Then, the development of two independent ANFIS models is explained in detail. The model validation part starts with the introduction of accuracy metrics. The developed models were evaluated based on the calculated error values.

Site measurements
The site investigation was performed in an Istanbul Ayazaga region aggregate quarry. The formation dominantly contains lower carboniferous aged sandstones (greywackes), shale, siltstone and claystone. Andesite and diabase dikes are commonly observed. The thickness of the layers is between 10 cm and 2.5 m for especially sandstones and siltstones. The quarry generally produces coarse and medium-grained greywacke-sandstones. The average grain size is between 0.2 and 0.4 mm. The rock contains 60-75% quartz and lesser amount of feldspars, muscovite, chlorite and hematite minerals (Tugrul and Undul 2006;Ozgul 2012).
The rock formation is highly fractured and deformed in Akdaglar quarry. The rock quality can be described as poor to fair. In some parts, good rock quality is observed. The strike directions of the joints are mainly NW-SE and NE-SW. The aperture of joints ranges from 0.1 to 1.0 mm. The median in-situ block size is 1.21 m. The density of the rock is between 2.65 and 2.74 g/cm 3 . The compressive strength of the rock is between 61 and 77.5 MPa. The elasticity modulus is 15-16.9 GPa. The average porosity is 2.03%. The above mentioned rock properties were measured by the author (Hudaverdi and Akyildiz 2019;Akyildiz and Hudaverdi 2020).
The blasting operations are performed using ammonium nitrate based explosives and shock tube initiation systems. Figure 1 illustrates the operational parameters of bench blasting. All the measured design variables, a rock mass parameter and flyrock throw distance are presented in the Table 1. The columns 2-6 in Table 1 show the fundamental blast design parameters. These operational parameters are presented as ratios. Charge weight (CW) and specific charge (q) are the explosive parameters. They explain amount of explosive used per blasthole and per broken rock volume, respectively. Number of holes (N h ) and number of rows (N r ) were considered to represent blast dimensions ( Table 1). The selection of input parameters is in accordance with literature. All the recent studies, some of them are presented in the introduction, were reviewed. In addition, the fundamental references about rock blasting were examined (Konya and Walter 1990;Jimeno et al. 1995;Hustrulid 1999).
The tenth parameter is apparent in-situ block size (X B ) ( Table 1). Block size represents rock mass properties. The last parameter (AMF) is used to symbolize face conditions. Rock blasting is performed through a free face. Face condition directly affects efficiency and environmental effects. Face condition was classified from 1 to 4, based on amount of the material in front of the face. Code 1 represents a clean face. Code 4 indicates that amount of the material (unremoved muckpile) in front of the face exceeds half of the bench height. The last column of Table 1 shows the measured flyrock throw distance. Flyrock throw measurements were performed for 77 cases. Flyrock throw (FR) ranges from 28 to 142 m. The median FR is 74.0 m. The mean FR is 79.2 m with a standard deviation of 35.7 m.

Examination of the measured variables
Multicollinearity among input variables was investigated before the modeling stage. Figure 2 shows the multicollinearity matrix for the input variables. Several researchers indicated a correlation value of 0.9 or higher for the certain presence of multicollinearity (Dohoo et al. 1997;Gregorich et al. 2021;Darabi et al. 2021). In this research, all correlation values are lower than the threshold (0.9), and thus, there is no a multicollinearity problem. In addition, it should be noted that the study aims to create predictive ANFIS models. ANFIS is a non-linear technique. The books of Garson (2014) and Tabachnick and Fidell (2013) may be examined to review multicollinearity in linear models.
In the multicollinearity matrix, the last row and the last column show the relationship between flyrock throw distance and input variables. The highest correlations were obtained for specific charge (q) and B/D values. The calculated correlations are 0.85 and 0.76 for q and B/D ratio, respectively. All the other correlation values are lower than 0.5 (Fig. 2). Increase of specific charge means the increase of explosive usage per cubic meter broken rock. In our case, increase of specific charge results in increase of flyrock throw distance. B/D ratio is negatively correlated to the flyrock throw. It seems that increase of burden distance results in decrease of excessive flyrock throw. Burden is one of the most important parameters for flyrock throw. In general, results are in accordance with the fundamental blasting references (Konya and Walter 1990;Jimeno et al. 1995;Bhandari 1997 Faramarzi et al. (2014) investigated correlation between the measured parameters and flyrock. The highest observed correlation was 0.614 for H/B ratio. In this study, the correlation value between FR and H/B is 0.46 (Fig. 2). According to Dehghani et al.'s (2021) study, the correlation value between flyrock and specific charge was higher than 0.85. The other dominant parameters were the length of stemming, burden and spacing. Hasanipanah and Bakhshandeh Amnieh (2020) developed a fuzzy rule-based approach for prediction of flyrock. According to the sensitivity analysis, charge weight per hole was found to be most dominant parameters in the system. Burden, specific charge and H/B ratio were the other significant parameters. The present findings seem to be consistent with previous researches. On the other hand, it should be noted that blast researchers select the input parameters based on available data. Type and number of input parameters vary for each study. The researchers also apply different techniques to predict flyrock.
At this stage, some controversial debates on multicollinearity may be reviewed. In recent years, some researchers propose that multicollinearity is not a sufficient reason to remove variables from a model. According to Gregorich et al. (2021), the collinearity investigation may fail in determining suitable variables for a model. Multicollinearity results may misguide the researcher. An unnecessary variable may be kept in the model. Conversely, an essential variable, which is needed for proper adjustment in an explanatory model, may be omitted. Usage of correlated independent variables may increase standard errors. However, this possibility should be accepted to correctly answer a research problem. Some multivariate models should be constructed considering all the relevant parameters. Excluding relevant parameters may cause formation of incorrect and inconsistent forecasting models. O'Brien (2017) proposed that collinearity is not an acceptable reason to exclude an input parameter from a model. In some cases, removing a variable seems to be an easy solution. However, variable omission may cause incomplete evaluation of an engineering problem. The developed model may not support the researcher's asserted hypothesis. Gregorich et al. (2021) cited Hernan and Robins (2020) and suggested causal inference, such as conditional exchangeability and positivity for some particular cases.
A researcher should choose the variables based on the experience, the research topic and hypothesis. For example, in our case, the blast parameters were selected based on site observations and the researcher's experience. There is a logic behind the variable selection procedure. The input variables were determined considering the three categories: fundamental blast design parameters, parameters related to explosive consumption and parameters for blast dimension.

Research methodology
The measurement technique and research methodology were determined considering recent equipment and software. In addition, blasting and earth science literature was reviewed. The on-site observations revealed that maximum flyrock distance was due to face burst mechanism. Thus, the maximum horizontal distance between blast face and landed fragments was measured as flyrock throw distance. Similar measurement technique was also conducted by Ghasemi et al. (2012), Faradonbeh et al. (2018), Rad et al. (2018) and Nikafshan Rad et al. (2020). Theoretical modelling studies also expect maximum flyrock throw in front of the bench face (Richards and More 2004). Tape measure and hand-held GPS equipment were used to determine flyrock throw distance. Only the fragments, which had the capacity to cause damage, injury or fatality, were considered. Based on past investigations in the Akdaglar quarry and the literature related to flyrock (Little 2007;Ghasemi et al. 2012;Dehghani et al. 2021); these fragments have a diameter of about 10 cm.
Blasting is defined as a transformation from the state of in situ block size distribution to the state of blasted block size distribution . In this study, apparent in-situ block size was measured as a rock mass parameter. The benches were analyzed by the Wipjoint joint analysis software. The scaled images of the bench faces were captured before each blast. The images had details of all the small and large discontinues and apparent rock blocks. The images were delineated meticulously and average in situ block size was determined. The joint analysis software relies on the concept of stereology. The measured blocks are modeled matching equivalent sphere. In-situ block size is given as diameter of the equivalent sphere (Maerz 1996). Similar block size determination technique has been applied by several blast researchers (e.g., Hudaverdi et al. 2012;Azizi and Moomivand 2021). In general, rock mass contains several discontinuity sets. Spacing between joints, joint aperture, length of discontinuities, orientation and joint filling material are some of the properties of discontinuities. The formation of in-situ blocks of intact rock are determined by the mutual intersection of discontinuity sets with different spacing and orientation characteristics. The in-situ block sizes and shapes are governed mainly by the spacing, the persistence of discontinuities and the number of discontinuity sets Latham and Lu 1999). This research focuses on flyrock. On the other hand, blast fragmentation was also observed in the studied quarry. The observed blasts were quite successful. There were no excessive boulders in the muckpile. As a further reading, Hudaverdi et al.'s (2012) and Akyildiz and Hudaverdi's (2020) blast fragmentation researches may be examined. These investigations were also conducted in the studied Akdaglar quarry.
The flowchart of the research is given in Fig. 3. The first stage of the research is the measurement of input parameters in the quarry. After that, multivariate analysis procedures were applied for variable selection. IBM SPSS Statistics™ software was used to perform stepwise regression and factor analysis. Multivariate analysis helped to understand dominant parameters of flyrock occurrence. Stepwise regression and factor analysis indicated different parameters. The output of regression analysis is used to create ANFIS Model 1. The ANFIS Model 2 was created based on the results of factor analysis (Fig. 3). Both ANFIS models were developed by Matlab™ software. In the last stage, the developed ANFIS models were validated on the independent test blasts.
In this study, adaptive neuro fuzzy system (ANFIS) was selected to create predictive models. The ANFIS method combines the advantages of both artificial neural networks and fuzzy inference systems (Singh et al. 2012). ANFIS is a nonlinear modeling technique. It has ability to adapt quickly to the nonlinear mechanism of a process. It is a flexible method and has a rapid learning capacity (Mostafaei 2018). In addition, ANFIS is a well-proven modelling method. The background literature of ANFIS is relatively large in comparison to some other recent modeling techniques. such as particle swarm optimization (PSO), grey wolf optimizer (GWO), cuckoo search algorithms (CS) etc.

Stepwise regression
The first variable selection method was stepwise regression. In stepwise regression, model development is a sequential process. Independent variables are added to the model considering partial correlations. First the independent variable which is mostly correlated to dependent variable is included the model. In the second step, the other variables, with the highest partial correlation with the dependent, is added to model. This procedure is repeated until the addition of a remaining variable does not contribute to the equation (Garson 2014; Hudaverdi and Akyildiz 2019). The contribution of each variable is measured by considering the increase of R-squared value. Stepwise regression is generally applied to select a subset of independent variables in predicting the dependent variable. It eliminates the independent variables those do not contribute sufficiently to forecasting process (Tabachnick and Fidell 2013). Therefore, this technique may be considered as a useful tool for variable selection.
Eleven blast parameters given in Table 1 are the independent variables of the model. The dependent variable is flyrock throw distance (FR). Table 2 shows regression summary. The first model contains single parameter, B/D and an equation constant. Addition of X B increases the coefficient of determination (R 2 ) to 0.785. The standard error of the estimate decreases below 17. Addition of the specific charge (q) forms the final model. Equation (1) shows the obtained multiple regression equation. The adjusted R 2 values are smaller than R 2 values. In general, the increase rate of adjusted R 2 is slightly higher than R 2 in each step (Table 2).
(1) FR = −8.446 × (B∕D) − 7.527 × (X B ) + 94.520 × (q) + 239.740  The variable selection process of the stepwise regression is explained in Table 3. In the first phase, B/D ratio enters the model. Then, the remaining variables are examined. X B has the highest partial correlation (−0.307). X B is added to the equation. At this stage, there are nine variables remaining. Specific charge (q) has the highest correlation and it is included to the model ( Table 3). Addition of the other variables does not increase the coefficient of determination (R 2 ) significantly. Model development process is finalized (IBM SPSS Statistics Base 2016).

Factor analysis
The second variable selection method is factor analysis. This method can be applied to represent a large number of independent variables as a smaller number of factors. Variable subgroups can be created from a large data set. In addition, the method may be used to aggregate the cases or extreme values (Garson 2012). Factor analysis uses different factor extraction methods, such as principal component analysis, unweighted least squares, maximum likelihoods, alpha factoring and principal axis factoring. In this study, principal component analysis, which is most widely practiced technique in engineering literature, was selected as extraction method. Principal component analysis is appropriate for the purpose of data reduction to obtain the minimum number of factors to represent the studied data set (Ho 2014). During modelling, a rotation process is applied to create more understandable results. Rotation increases the strength of high correlations between factors (Tabachnick and Fidell 2013). Figure 4 illustrates results of the factor analysis. Eleven components (factors) were extracted. Figure 4a presents the eigenvalues calculated for each component. Towards the right of the graph, the eigenvalue decreases. In general, the components on the steep part of the slope are selected (IBM SPSS Statistics Base 2016). The components on the shallow slope contribute very weakly to target model. Four factors have an eigenvalue higher than 1. The factors with high eigenvalues represent the data more strongly. In our case, four factors were chosen. Figure 4b shows the variance explained by the factors cumulatively. For example, the cumulative variance for the third component is the sum of the percentage of variance for the first, second and third components. The first four factors represent 72.012% of the variability in the eleven blast parameters. Considering the eigenvalues and cumulative variance, four factors were selected to represent the database. Decision on number of extracted factors was largely discussed by Tinsley and Brown (2000). Table 4 is named as rotated component matrix table. It shows what the components (factors) explain. Italized cells indicate the selected variables. First component is highly correlated with B/D and q. B/D was selected from the first component to represent blast design parameters. From the second component, N h was selected to obtain an input variable that represents blast dimension. In the third component, two variables have a correlation value higher than 0.80 (Table 4). H/B was selected as a blast design parameters. CW was selected to represent explosive amount. There is no any parameter with a correlation higher than 0.800 in the fourth component. The variable selection process is finalized.

Fundamentals of ANFIS model
The ANFIS modelling process starts with definition of the problem. Input parameters and target parameters of forecasting process are determined. Then, fuzzification stage is started and fuzzy rules are defined. The next stage is the  (Ismail and Bendary 2020). Figure 5 shows the basic architecture of ANFIS model (Jang 1993). As shown in Fig. 5, there are two inputs and a single output parameter. The architecture consists of five different layers. Modelling process starts with the definition of two if-then rules: Rule 1: If x is A 1 and y is B 1 , then f 1 = p 1 x + q 1 y + r 1 , Rule 2: If x is A 2 and y is B 2 , then f 2 = p 2 x + q 2 y + r 2 , Where A 1 , A 2 , B 1 and B 2 are the lingustic variables for the parameters x and y; p 1 , q 1 , r 1 , p 2 , q 2 , r 2 are the output function variables (Jang 1993).
The composition of the layers in Fig. 5 can be described as follows: • Layer 1: A special node function is described for each node "i" in the layer 1: This layer is also named as fuzzification layer. Where x values are input nodes and A is a linguistic definition (small, medium, large, etc.). A i is the membership function of the linguistic variables. • Layer 2: In this multiplication layer, the shape of the nodes is circular and they are shown by the letter "β". Each node multiples the received data. The output of multiplication layer is called the firing strength: • Layer 3: Neurons are shown by symbol "ε" and identified by a circle. The main function of the third layer is normalization. ith rules firing strength is divided by summation of firing strength of all rules to get an output (Nezhad and Jafari 2020): • Layer 4: This layer is referred to defuzzyfing layer. All the nodes are adaptive type. The specific equation of the layer is given as: P i , q i and r i are the consequent parameters, w i is the normalized firing strength. • Layer 5: The last process is summation. Sum of the all incoming signals are calculated to reach an overall outcome. (Jang 1993;Rai et al. 2015):

ANFIS Model 1
The first ANFIS model was constructed based on the results of stepwise regression. The input parameters are B/D ratio, apparent in-situ block size (X B ) and specific charge (q). The general architecture of the proposed model is presented in Fig. 6. After input layer, there are three sequential sections. These are input membership function (MFs), rule and output layers, respectively. The last layer is referred to output layer. Number of rules is determined according to number of input membership functions (Mathworks 2018). Multiplication of input membership functions (2 × 2 × 5) creates twenty rules (Fig. 6). The type of output membership function is constant. Figure 7 shows membership functions for input parameters. Several membership functions, e.g., trapezoidal, generalized bell shape, gaussian, triangular, pi-shaped, are available for ANFIS modeling (Sihag et al. 2019). After the determination MFs, data are grouped based on linguistic variables, such as low, medium and high. All the possibilities were considered for determination of the optimum membership functions. The detailed trials were performed for different type of functions and for different subgroups. Generalized bell shape MFs was selected for input parameters (Fig. 7). Bell shape and trapezoidal MFs are widely used in the engineering literature (e.g., Mottahedi et al. 2018;Suthar 2020). B/D and X B graphs indicate two subgroups (Fig. 7a,b). MFs of specific charge (q) were divided into 5 subgroups from very low to very high (Fig. 7c). Figure 8 shows error graph during model training phase. Root mean square error is used as error metric. As seen in the figure, the error value constantly decreases until the fortieth epoch (repetition). The minimum obtained error value is 11.15.

ANFIS Model 2
The second ANFIS model was created based on the results of factor analysis. There are four input parameters: B/D, H/B, number of holes (N h ) and charge weight (CW). Increase in number of input parameters makes the model relatively more complex. 36 rules were created for fuzzification process. The membership function parameters of fuzzy inference system are controlled (adjusted) using a back propagation algorithm. Figure 9 presents general structure of the second ANFIS Model. Different type of membership functions were tested for input parameters. The output membership function is constant type. The membership functions of ANFIS Model 2 are presented in Fig. 10. MFs, sometimes called fuzzy reasoning, are the core of the fuzzy concept. Human thinking power and reasoning are imitated by membership function (Alipour and Ashtiani 2011). The selected MF for ANFIS 2 is trapezoidal. The trapezoidal membership function has a flat top and it looks like a truncated triangle. This membership function is in the form of straight lines and sometimes preferred for its simplicity (Mathworks 2018). B/D and N h parameters were divided into two subgroups (Fig. 10a, c). H/B and CW were classified into three subgroups; low, medium and high (Fig. 10b, d). Figure 11 shows the rule viewer window of the Matlab software. The first four columns are the input parameters. The last column shows the output; flyrock throw distance (FR). Each row indicates a different rule. The figure represents the first 25 rules. The rule viewer explains the mechanism of the fuzzy inference system (Güneri et al. 2011). If the red lines are moved by researcher, the system changes the value of input variables. Fuzzy inference system automatically calculates the corresponding output (Okwu and Tartibu 2020). In the figure, H/B, B/D, CW and N h values are 3.33, 24.2, 29.1 and 31.5, respectively. The calculated flyrock throw is 91.9 m. This value corresponds to a medium-high flyrock throw distance.

Selection of accuracy measures
Different accuracy measures were considered for model validation. Absolute errors and percentage errors are largely applied in engineering literature. Absolute errors are scale dependent and they cannot be used to compare different data series (Shcherbakov et al. 2013). Some researchers criticize percentage errors because of their unsymmetrical structure (Armstrong and Collopy 1992). Sometimes, error measures are largely affected by extreme values (outliers). Actually, each accuracy measure has some limitations. In this research, in addition to classical error metrics, recently suggested alternative measures were also applied. For example, in addition to regular absolute and percentage errors, standardized accuracy measures were also used (Casella and Berger 1990;Kusiak et al. 2010). A normalized root mean square error was calculated. NSE (Nash and Sutcliffe 1970) was proposed as an alternative goodness-of-fit index. All the suggested accuracy metrics are presented in Table 5.

Validation of the models on the test blasts
Thirty-two independent test blasts were used for model validation. The test blasts were also conducted in the studied quarry. All the test cases are presented in Table 6. Flyrock throw (FR) ranges from 31 to 140 m. The median FR are 84.0 m. Figure 12 shows the measured and predicted values  by the models. The arrows in the figure were replaced to help the readers to follow the details of the forecasting process. They show the highly accurate predictions performed by the models. Red and blue arrows, for ANFIS 1 and ANFIS 2, respectively, indicate the predictions achieved by an error lower than 2 m. As seen in the figure, models are quite successful. ANFIS 1 made successful predictions especially for the cases 7, 8, 10, 24 and 31. The largest deviation of ANFIS 1 is 11.94 m for Blast 20. ANFIS 2 made very accurate predictions for the cases 25-29 (Fig. 12). Predictions of regression largely deviate for the cases 11-15. Fourteen and sixteen cases were predicted with an error lower than 5 m by ANFIS 1 and ANFIS 2, respectively. Table 7 shows the calculated error values for the developed models. ANFIS 1 has the lowest mean absolute error values. The median absolute errors (MdAE) for ANFIS 1 and 2 are lower than 5 m. In general, median accuracy metrics are lower than mean accuracy metrics. Standardized absolute errors (Std_AE) and absolute percentage errors (Std_APE) clearly indicate superiority of the ANFIS models. The correlations values of ANFIS models are higher than the regression equation. However, the correlations are very close to each other. NSE indicates the superiority of ANFIS 1 more clearly. The NSE values are 0.970 and 0.959 for ANFIS 1 and 2, respectively (Table 7).

Discussion
The flyrock measurements were performed in a sandstone quarry and predictive models were created. The measured values may be considered normal. No excessive flyrock throw was observed during site measurements. Each study appearing in the literature has been performed in different The ANFIS models were validated on independent test blasts. The calculated mean absolute errors for the predictions are 5.40 and 5.73 for ANFIS 1 and 2, respectively. The median absolute errors are lower than 5 m. The NSE values of ANFIS models are higher than 0.95. The MAPE values are lower than eight percent. Evidently, the both ANFIS models perform better than the regression equation. Performance of ANFIS 1 is better than that of ANFIS 2 according to 10 out of 11 accuracy metrics.
At this stage, coefficient of determination (R 2 ) may be used to compare the developed ANFIS models and some recent models appearing in the literature. Blasting researchers mostly use coefficient of determination (R 2 ) to investigate prediction capability of the developed models. R 2 indicates the accord between measured and predicted flyrock values. Faradonbeh et al. (2018) andNikafshan Rad et al. (2020) were obtained a R 2 value of 0.924-0.944 by applying gene expression programming (GEP) and genetic algorithm ANN (GE-ANN) methods. Koopialipoor et al. (2019) achieved a R 2 value of 0.959 using particle swarm optimization ANN (PSO-ANN) model. Ye et al.'s (2021) genetic programing model (GE) has a R 2 value of 0.908. In the current study, the R 2 values are 0.972 and 0.964 for ANFIS 1 and ANFIS 2 models, respectively. It is clear that the R 2 values obtained by ANFIS models are quite satisfactory. The median absolute error of Hasanipanah et al.'s (2018) rock engineering system is 11.64 m. In this study, median absolute errors are 4.86 and 4.98 m. Lu et al. (2020) applied outlier robust extreme learning machine (ORELM) and obtained a mean absolute error value of 10.144 m. The mean absolute error values of the current study are quite low; 5.40 and 5.73 m. At this stage, it should be noted that comparison of the correlation values provides a basic evaluation. The mentioned studies consider similar blast variables. However, not all the variables are exactly same. In addition, each study appearing in the literature was performed in a different rock formation.
In that part, the geological condition of the blast benches may be briefly discussed. Bench faces were meticulously examined during blast measurements. Any unexpected geological structures such as cavities, voids, incompetent zones and clay seams were monitored. In addition, discontinues were examined. Discontinuities may play important role in flyrock occurrence. The blasts performed in massive or heavily fractured rocks may result in different level of flyrock throw. In this study, apparent in-situ block size was measured by image analysis software for each blast. Seventy four percent of the measured block size values are between 0.5 and 2 m. In our case, there is no a significant correlation Normalized root mean square error nRMSE = 1 Standardized absolute percentage error Pearson's correlation between in-situ block size and flyrock throw distance. The calculated correlation value is 0.30. Drilling operations were also observed in the quarry. The site measurements were started before explosive charging operation. No large cavities or in-hole collapses were observed.

Conclusions
Eleven different variables, including a rock mass parameter, were measured in a sandstone quarry. Operational parameters such as burden, spacing and stemming length were considered. Number of holes and number of rows were used as dimensional parameters. Charge weight and specific charge values were also calculated. Two alternative ANFIS models were created to predict flyrock throw distance. The input parameters of ANFIS 1 are B/D ratio, in-situ block size and specific charge. The input parameters of ANFIS 2 are B/D and H/B ratios, number of holes and charge weight per hole. The results revealed that ANFIS is a very efficient technique to predict flyrock throw distance. Although ANFIS 2 has four input parameters, ANFIS 1 with three inputs performed slightly better than ANFIS 2. The first model use in-situ block size as a rock mass parameter. ANFIS 2 mostly relies on operational blast parameters. The both models use burden/hole diameter ratio as an input parameter. It is believed that all forecasting models of flyrock should consider the burden. In addition, the other fundamental blast parameters should be included the models as much as possible. The developed ANFIS models are highly successful for the studied quarry. However, the models should be reconstructed for different rock formations. All the models that developed to predict blast fragmentation, flyrock or ground vibration are site specific (Kahriman 2004;Giraudi et al. 2009;Hudaverdi 2012, etc.). Factor analysis, regression analysis and ANFIS are non-rigid and very adaptive methods. The variable selection concept presented in this paper can be easily adapted to a particular situation. This study tries to develop a framework. Open pit mines can create different novel models based on the concept explained in this research. It is possible to repeat the proposed methodology by considering new blast parameters and using new blast databases in different geological conditions.
The study proved that application of pre-statistical analysis increases the performance of soft computing models. The developed ANFIS models have three and four inputs. Non-complex and robust models were created. In flyrock forecasting, predictive models should be constructed after a variable reduction process. The researcher should also examine the possibility of the creation of multiple models. This research created two parallel models for the studied sandstone quarry. It is possible to develop more than a single model for a specific site.
The determination of the membership functions plays an important role on the success of the developed ANFIS models. Increase of the input parameters and membership functions means increase in number of the rules. Large number of rules may make the models unnecessarily complex. Researchers tend to use a few accuracy metrics for model validation. In validation stage, eleven accuracy metrics were used. In addition to classical accuracy measures, standardized and normalized errors were found to be useful for model validation. NSE indicated the superior  of ANFIS models clearly. Application of NSE should be increased as an alternative goodness-of-fit index. In this research, flyrock cases were investigated in an aggregate quarry. Quarries use small diameter blast holes. Explosive charge per hole is generally low. The charge weight ranges from 17.97 to 40.17 kg in the studied blasts. Blasting operations were performed in low benches. The mean bench height is 7.08 m with a standard deviation of 0.90 m. In the future, flyrock studies should be performed in collieries or large open pits. Coal mines and metal mines use large blast holes. Size of the blasts is also larger than quarries. Higher and excessive flyrock throw cases may be observed and studied.