Practical Risk Assessment of Ground Vibrations Resulting from Blasting, Using Gene Expression Programming and Monte Carlo Simulation Techniques

Peak particle velocity (PPV) is a critical parameter for the evaluation of the impact of blasting operations on nearby structures and buildings. Accurate estimation of the amount of PPV resulting from a blasting operation and its comparison with the allowable ranges is an integral part of blasting design. In this study, four quarry sites in Malaysia were considered, and the PPV was simulated using gene expression programming (GEP) and Monte Carlo simulation techniques. Data from 149 blasting operations were gathered, and as a result of this study, a PPV predictive model was developed using GEP to be used in the simulation. In order to ensure that all of the combinations of input variables were considered, 10,000 iterations were performed, considering the correlations among the input variables. The simulation results demonstrate that the minimum and maximum PPV amounts were 1.13 mm/s and 34.58 mm/s, respectively. Two types of sensitivity analyses were performed to determine the sensitivity of the PPV results based on the effective variables. In addition, this study proposes a method specific to the four case studies, and presents an approach which could be readily applied to similar applications with different conditions.


Introduction
In mining and civil engineering projects, a common technique to remove rock mass is to blast the rock mass. This results in a problem of wasted blasting energy, which occurs for different reasons, and brings about a number of environmental issues, like ground vibrations, backbreak, dust and fumes, air overpressure and flyrock [1][2][3][4][5]. One of the most adversative effects is ground vibration (GV) which may cause damage to neighbouring structures [6,7]. In the literature, GV is defined generally as a wave motion that can be extended from the blast point to adjacent zones [8].
Thus, blast-induced GV needs to be predicted, and its risks should be analysed and assessed in a way to enable making the best decisions about the design of blasting operations and minimize associated environmental issues.

Case Study and Data Collection
In this study, in order to apply methods for PPV evaluation, four quarry blasting sites in Malaysia were selected. The names of these quarries and their latitudes and longitudes are presented in Table 1. In addition, geological map and locations of these quarries are shown in Figure 2. Generally in these quarries, the type of rock is granite, and blast operations are conducted in order to create aggregate to be used for construction purposes. In each site, roughly 30,000 tonnes of aggregate is produced per month. A value of 10 m was measured as the minimum bench height in the Kulai Quarry, while a value of 28 m was observed as maximum bench height in the Bukit Indah Quarry. Various weathering zones, ranging from moderately weathered to completely weathered, were observed in the studied sites based on the suggested method by ISRM [72], and as reported by Abad et al. [73]. Results of the Schmidt hammer test (based on ISRM [72]) were measured in the range of 19-37. In addition, in these sites, rock quality designation (RQD) as an indicator of rock mass fracturing was measured as a minimum of 22.5% to a maximum of 61.25%. According to previous studies [74,75], resistance parameters of the rock mass can have a deep impact on explosive operations.

Case Study and Data Collection
In this study, in order to apply methods for PPV evaluation, four quarry blasting sites in Malaysia were selected. The names of these quarries and their latitudes and longitudes are presented in Table 1. In addition, geological map and locations of these quarries are shown in Figure 2. Generally in these quarries, the type of rock is granite, and blast operations are conducted in order to create aggregate to be used for construction purposes. In each site, roughly 30,000 tonnes of aggregate is produced per month. A value of 10 m was measured as the minimum bench height in the Kulai Quarry, while a value of 28 m was observed as maximum bench height in the Bukit Indah Quarry. Various weathering zones, ranging from moderately weathered to completely weathered, were observed in the studied sites based on the suggested method by ISRM [72], and as reported by Abad et al. [73]. Results of the Schmidt hammer test (based on ISRM [72]) were measured in the range of 19-37. In addition, in these sites, rock quality designation (RQD) as an indicator of rock mass fracturing was measured as a minimum of 22.5% to a maximum of 61.25%. According to previous studies [74,75], resistance parameters of the rock mass can have a deep impact on explosive operations. The total number of conducted blasting operations which were investigated in these sites was 149 operations. In these operations, different blasting parameters, including burden (the shortest distance between the hole and the exposed bench face), spacing (the distance between blast-hole rows), stemming (the material which is placed on top of explosives in blast holes), maximum charge per delay (the maximum amount of explosive charge detonated on one delay within a blast with units defined as kg), powder factor (the amount of explosive needed to remove and a fragment of 1 m 3 of the rock mass with its unit defined as kg/m 3 ) and distance from the blast face, were recorded. In all operations, a value of 115 mm was utilized for the blast-hole diameter. It is of a high importance to measure, evaluate and simulate PPV in the investigated sites, because blasting operations are carried out at a close distance to adjacent residential buildings. The nearest building is located 350 m away from the investigated sites; therefore, an area in a distance ranging from 65 m to 640 m to the blast face was determined for PPV measuring stations. All of the obtained PPV values were recorded opposite the quarry's bench, and almost perpendicular to it, and the maximum PPV levels in a variety of directions were taken into consideration. The total number of conducted blasting operations which were investigated in these sites was 149 operations. In these operations, different blasting parameters, including burden (the shortest distance between the hole and the exposed bench face), spacing (the distance between blast-hole rows), stemming (the material which is placed on top of explosives in blast holes), maximum charge per delay (the maximum amount of explosive charge detonated on one delay within a blast with units defined as kg), powder factor (the amount of explosive needed to remove and a fragment of 1 m 3 of the rock mass with its unit defined as kg/m 3 ) and distance from the blast face, were recorded. In all operations, a value of 115 mm was utilized for the blast-hole diameter. It is of a high importance to measure, evaluate and simulate PPV in the investigated sites, because blasting operations are carried out at a close distance to adjacent residential buildings. The nearest building is located 350 m away from the investigated sites; therefore, an area in a distance ranging from 65 m to 640 m to the blast face was determined for PPV measuring stations. All of the obtained PPV values were recorded opposite the quarry's bench, and almost perpendicular to it, and the maximum PPV levels in a variety of directions were taken into consideration.  To effectively evaluate the PPV results, the parameters with the highest effects must be determined. As noted earlier, the two parameters of distance from the blast site and the maximum charge per delay have been widely taken into account by numerous scholars in the PPV prediction process. Furthermore, based on a number of studies available in the literature [7,11,22,60,76,77], powder factor and blasting pattern parameters play an important role in determining the PPV value. Therefore, in this study, the burden-to-spacing ratio (BS), stemming (ST), powder factor (PF), maximum charge per delay (MC) and the distance from the blast face (D) were considered as predictors or model inputs. The average, maximum, minimum, unit and symbol of the measured parameters are shown in Table 2. It is important to note that in the case of parameters which had several values in a blast (such as stemming length), their average values were considered and used in the established database.  To effectively evaluate the PPV results, the parameters with the highest effects must be determined. As noted earlier, the two parameters of distance from the blast site and the maximum charge per delay have been widely taken into account by numerous scholars in the PPV prediction process. Furthermore, based on a number of studies available in the literature [7,11,22,60,76,77], powder factor and blasting pattern parameters play an important role in determining the PPV value. Therefore, in this study, the burden-to-spacing ratio (BS), stemming (ST), powder factor (PF), maximum charge per delay (MC) and the distance from the blast face (D) were considered as predictors or model inputs. The average, maximum, minimum, unit and symbol of the measured parameters are shown in Table 2. It is important to note that in the case of parameters which had several values in a blast (such as stemming length), their average values were considered and used in the established database.

Gene Expression Programming (GEP)
To obtain a mathematical relation, the method of GEP has been developed for the prediction of PPV. This method was firstly introduced by Ferreira [78]. Recently, the application of this method has increased in civil and mining engineering [7,67,[79][80][81][82][83][84]. With the use of different mathematical operators and linear functions (e.g., +, −, sin, tan, cos, sqrt, x 2 , and x 3 ), this method forms a mathematical relation between the input and output variables. In the following, a summary of the GEP modelling process is presented: Step 1 On the basis of population number, a certain number of chromosomes are produced in a random way.
Step 2 The chromosomes of the initial population are denoted as mathematical equations.
Step 3 Each chromosome's fitness is measured based on the fitness function (coefficient of determination, R 2 , or root mean square error (RMSE)). If the stopping criteria are not met, the best of the first generation is chosen using the roulette wheel method.
Step 4 The genetic operators, which are known as the core of the GEP algorithm, are applied to the rest of chromosomes for the purpose of creating modified individuals.
Step 5 At this step, the chromosomes will start creating the next generation. This process is iterated for a certain number of generations. More information about this method can be found in the work of other researchers [78,80,83,84].
Various models of GEP were applied in which non-normal distributions of data were used. A total of 80% of the data was allocated to the training section and 20% to the testing section. Regarding the use of two statistical indicators of R 2 and RMSE [57,[85][86][87] for determining the best model, a scoring system was used, as introduced by Zorlu et al. [88]. The highest and lowest values of R 2 gain the highest and lowest score, respectively, while the situation is reversed for RMSE. For example, in the training section, values of R 2 for models 6 and 2 are 0.7283 and 0.6927, respectively, and the highest (10) and lowest (1) scores are obtained for them. However, for the RMSE of the training, the lowest error is related to model 4, and thus obtains the highest score (10), and model 10 also obtains the lowest score because of its maximum error value. This method of ranking is applied in different sections, and at the end, the total score of each row is placed into the last column. As a result, model number 6 obtains the highest score and is introduced as the best model for predicting PPV ( Table 3). As a result, R 2 of 0.7283 and 0.6823 and RMSE of 4.2814 and 4.0344 for training and testing, respectively, show suitable levels of prediction for developing the GEP model.
The best PPV equation obtained by the GEP model is presented in Equation (6). This equation is a summation of Equations (1) to (5), which are presented as follows: where, the variables of d(0), d(1), d(2), d(3), d (4), and Y are Burden to Spacing, Stemming, Powder Factor, Max charge Per Delay, Distance and PPV, respectively.   (6), for training and testing, respectively.      (6), for training and testing, respectively.

Background
MC simulation is a quantitative risk assessment technique which considers the effect of deterministic and probabilistic variables in a single model [70]. In other words, MC simulation is capable of analysing the effect of uncertain variables to predict or assess the risk. MC simulation has been adopted in various engineering fields due to its capabilities in predicting the effect of uncertain variables. In an MC risk analysis model, random sampling techniques are used based on the requirements of the model, and statistical analyses are conducted to ensure accurate results are obtained. MC simulation is developed for probabilistic circumstances in which deterministic values of variables are not provided, so that historical data can be used with various distribution functions.
According to Figure 5, which illustrates the operation principal of MC simulation, multiple variables are involved in an MC model with various ranges and distribution functions [89]. Although a random number is chosen for the variables at every iteration, the variables' distribution functions contribute in this random process. The number of required iterations is defined before conducting the analysis. Once the results of all iterations are obtained, their range of the outputs can be statistically analysed to predict the result. As a result, in an MC model, the inputs and outputs are ranges of values, rather than a deterministic value [90]. Notably, an MC simulation is capable of considering independent and dependent variables in a single model. Because there might be several variables involved in an MC model, in order to achieve more accurate results, the relationships among them should be investigated [91,92].

Background
MC simulation is a quantitative risk assessment technique which considers the effect of deterministic and probabilistic variables in a single model [70]. In other words, MC simulation is capable of analysing the effect of uncertain variables to predict or assess the risk. MC simulation has been adopted in various engineering fields due to its capabilities in predicting the effect of uncertain variables. In an MC risk analysis model, random sampling techniques are used based on the requirements of the model, and statistical analyses are conducted to ensure accurate results are obtained. MC simulation is developed for probabilistic circumstances in which deterministic values of variables are not provided, so that historical data can be used with various distribution functions.
According to Figure 5, which illustrates the operation principal of MC simulation, multiple variables are involved in an MC model with various ranges and distribution functions [89]. Although a random number is chosen for the variables at every iteration, the variables' distribution functions contribute in this random process. The number of required iterations is defined before conducting the analysis. Once the results of all iterations are obtained, their range of the outputs can be statistically analysed to predict the result. As a result, in an MC model, the inputs and outputs are ranges of values, rather than a deterministic value [90]. Notably, an MC simulation is capable of considering independent and dependent variables in a single model. Because there might be several variables involved in an MC model, in order to achieve more accurate results, the relationships among them should be investigated [91,92].

Application of MC Simulation in PPV Estimation
This section presents the modelling steps of MC techniques in the risk analysis of PPV. According to Figure 6, in this paper, six steps are considered for PPV estimation using MC simulation. First, the affecting variables (ST, BS, MC, PF and D), and their corresponding values which were gathered from quarry sites in Malaysia, were identified. Second, a predictive model which was carried out using GEP (Equations (1)-(6)) was developed. Third, the distribution function of each variable was considered. Because a large amount of data was gathered already, the distribution functions of all variables were calculated using the Risk Solver Platform [89]. The Risk Solver Platform is an MS-Excel plug-in feature which is used in this research for MC simulation. In this Risk Solver Platform, a 'best-fit' function is available to select the distribution function of the variables. The inputs for the MC simulation are presented in Table 4, together with their minimum and maximum amounts and the distribution functions. Notably, continuous probability distribution (CPD) was considered for all inputs during the use of this function.

Application of MC Simulation in PPV Estimation
This section presents the modelling steps of MC techniques in the risk analysis of PPV. According to Figure 6, in this paper, six steps are considered for PPV estimation using MC simulation. First, the affecting variables (ST, BS, MC, PF and D), and their corresponding values which were gathered from quarry sites in Malaysia, were identified. Second, a predictive model which was carried out using GEP (Equations (1)-(6)) was developed. Third, the distribution function of each variable was considered. Because a large amount of data was gathered already, the distribution functions of all variables were calculated using the Risk Solver Platform [89]. The Risk Solver Platform is an MS-Excel plug-in feature which is used in this research for MC simulation. In this Risk Solver Platform, a 'best-fit' function is available to select the distribution function of the variables. The inputs for the MC simulation are presented in Table 4, together with their minimum and maximum amounts and the distribution functions. Notably, continuous probability distribution (CPD) was considered for all inputs during the use of this function.  Fourth, because the current study is mainly aimed at the achievement of effective combinations in MC simulations, the input correlations were taken into consideration in the process of simulation. This program makes use of a correlation matrix between the model inputs and involves major rankorder correlations (Spearman's rho). In MC, the relations between model inputs have a deep impact on simulation results. Table 5 shows the mentioned relations obtained by Microsoft Excel Version 2013. Fifth, in developing the MC model, Latin hypercube sampling was used, due to its advantages compared with simple random sampling [94]. To ensure that every combination of variable values was analysed, 10,000 trials were performed. Then, the output ranges were generated as a sixth step.

Results and Discussion
In this study, in order to cover two different phases of prediction and simulation, GEP and MC techniques were proposed, respectively. A GEP model with R 2 values of 0.7283 and 0.6823 and RMSE values of 4.2814 and 4.0344 for training and testing, respectively, are able to provide a suitable accuracy level for PPV prediction.  Fourth, because the current study is mainly aimed at the achievement of effective combinations in MC simulations, the input correlations were taken into consideration in the process of simulation. This program makes use of a correlation matrix between the model inputs and involves major rank-order correlations (Spearman's rho). In MC, the relations between model inputs have a deep impact on simulation results. Table 5 shows the mentioned relations obtained by Microsoft Excel Version 2013. Fifth, in developing the MC model, Latin hypercube sampling was used, due to its advantages compared with simple random sampling [94]. To ensure that every combination of variable values was analysed, 10,000 trials were performed. Then, the output ranges were generated as a sixth step.

Results and Discussion
In this study, in order to cover two different phases of prediction and simulation, GEP and MC techniques were proposed, respectively. A GEP model with R 2 values of 0.7283 and 0.6823 and RMSE values of 4.2814 and 4.0344 for training and testing, respectively, are able to provide a suitable accuracy level for PPV prediction.
Then, the proposed GEP equation was used as an input for the MC technique for simulation purposes. The results of the MC simulation indicated that the minimum and maximum values of PPV are 1.13 mm/s and 34.58 mm/s, respectively, with a mean value of about 18.28 mm/s. According to the results, the distribution function of 10,000 outputs is 'logistic', which is illustrated in Figure 7. Based on the output values, with 90% confidence, the amount of PPV is less than 25.25 mm/s, and 99% of the results are less than 30.86 mm/s. In this study, the simulation of PPV was conducted after measurement on site and prediction using GEP. Because the means of PPV simulation and prediction are very close (18.28 mm/s and 16.23 mm/s, respectively), it can be stated that the accuracy of the MC model in the PPV simulation is acceptable. Figure 8 demonstrates the differences among the measured, predicted and simulated results of PPV. Furthermore, in order to have a better illustration regarding the contribution of the input values in the obtained output for each trial in MC simulation, Figure 9 plots the random amounts of all inputs, which have been taken by the MC simulator. This figure also shows their corresponding PPV results which were generated during MC simulation. It should be noted that Figure 9 was generated using the Risk Solver Platform, the MC simulator in our research, to show the involvement of each input in the PPV result. According to Duvall and Fogelson [95], as a main standard, PPV values of less than 50 mm/s can be considered as an acceptable value. In this paper, we can state with 100% confidence that the PPV for the blasting operations is in the allowable range.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 9 of 16 Then, the proposed GEP equation was used as an input for the MC technique for simulation purposes. The results of the MC simulation indicated that the minimum and maximum values of PPV are 1.13 mm/s and 34.58 mm/s, respectively, with a mean value of about 18.28 mm/s. According to the results, the distribution function of 10,000 outputs is 'logistic', which is illustrated in Figure 7. Based on the output values, with 90% confidence, the amount of PPV is less than 25.25 mm/s, and 99% of the results are less than 30.86 mm/s. In this study, the simulation of PPV was conducted after measurement on site and prediction using GEP. Because the means of PPV simulation and prediction are very close (18.28 mm/s and 16.23 mm/s, respectively), it can be stated that the accuracy of the MC model in the PPV simulation is acceptable. Figure 8 demonstrates the differences among the measured, predicted and simulated results of PPV. Furthermore, in order to have a better illustration regarding the contribution of the input values in the obtained output for each trial in MC simulation, Figure 9 plots the random amounts of all inputs, which have been taken by the MC simulator. This figure also shows their corresponding PPV results which were generated during MC simulation. It should be noted that Figure 9 was generated using the Risk Solver Platform, the MC simulator in our research, to show the involvement of each input in the PPV result. According to Duvall and Fogelson [95], as a main standard, PPV values of less than 50 mm/s can be considered as an acceptable value. In this paper, we can state with 100% confidence that the PPV for the blasting operations is in the allowable range.   Then, the proposed GEP equation was used as an input for the MC technique for simulation purposes. The results of the MC simulation indicated that the minimum and maximum values of PPV are 1.13 mm/s and 34.58 mm/s, respectively, with a mean value of about 18.28 mm/s. According to the results, the distribution function of 10,000 outputs is 'logistic', which is illustrated in Figure 7. Based on the output values, with 90% confidence, the amount of PPV is less than 25.25 mm/s, and 99% of the results are less than 30.86 mm/s. In this study, the simulation of PPV was conducted after measurement on site and prediction using GEP. Because the means of PPV simulation and prediction are very close (18.28 mm/s and 16.23 mm/s, respectively), it can be stated that the accuracy of the MC model in the PPV simulation is acceptable. Figure 8 demonstrates the differences among the measured, predicted and simulated results of PPV. Furthermore, in order to have a better illustration regarding the contribution of the input values in the obtained output for each trial in MC simulation, Figure 9 plots the random amounts of all inputs, which have been taken by the MC simulator. This figure also shows their corresponding PPV results which were generated during MC simulation. It should be noted that Figure 9 was generated using the Risk Solver Platform, the MC simulator in our research, to show the involvement of each input in the PPV result. According to Duvall and Fogelson [95], as a main standard, PPV values of less than 50 mm/s can be considered as an acceptable value. In this paper, we can state with 100% confidence that the PPV for the blasting operations is in the allowable range.

Sensitivity Analysis
In order to figure out the effect of each variable in the output values, correlation sensitivity analysis (CSA) and regression sensitivity analysis (RSA) were conducted for different purposes. It should be noted that both types of sensitivity analyses were performed using the Risk Solver Platform. In the CSA, the ranges of correlations were between −1 and +1. In CSA, MC simulation uses rank order correlation, rather than linear correlation, to calculate the relationship between two data sets by comparing the rank of each value in a data set. Rank order correlation has an acceptable performance if uncertain variables are involved in the simulation, or the distribution function of the variables are not specified precisely [96]. As shown in Table 6, the level of effectiveness of each variable in the final output of the simulation can be obtained from the results of CSA.

Sensitivity Analysis
In order to figure out the effect of each variable in the output values, correlation sensitivity analysis (CSA) and regression sensitivity analysis (RSA) were conducted for different purposes. It should be noted that both types of sensitivity analyses were performed using the Risk Solver Platform. In the CSA, the ranges of correlations were between −1 and +1. In CSA, MC simulation uses rank order correlation, rather than linear correlation, to calculate the relationship between two data sets by comparing the rank of each value in a data set. Rank order correlation has an acceptable performance if uncertain variables are involved in the simulation, or the distribution function of the variables are not specified precisely [96]. As shown in Table 6, the level of effectiveness of each variable in the final output of the simulation can be obtained from the results of CSA.
According to Table 6, D has the most impact on the results of PPV in an MC simulation, followed by ST, MC, PF and BS. It is worth mentioning that D and PF had negative impacts on the amount of PPV, while ST, MC and BS had positive impacts. In RSA, multiple regression analyses were conducted by varying one of the variable ranges (in this paper from −50% to +50%), while keeping the other variables constant [97]. It is noteworthy that these same processes were conducted for all variables to determine the influence of each variable on the output. The results of RSA indicated that the PPV result is highly sensitive to changes in D followed by ST and MC (see Figure 10). On the other hand, changes in the amount of BS have the lowest influence on the results. Results of the sensitivity analysis (especially regarding D, which is the most effective parameter on PPV) are in good agreement with studies carried out by [7,8,60]. For example, Monjezi et al. [60] proposed that the distance from the blast face is the most influential parameter on PPV among other used parameters (i.e., hole depth, stemming, maximum charge per delay). In addition, as discussed earlier, the distance from the blast face is one of the most effective parameters in empirical equations [15,17].
Appl. Sci. 2019, 9, x FOR PEER REVIEW 11 of 16 According to Table 6, D has the most impact on the results of PPV in an MC simulation, followed by ST, MC, PF and BS. It is worth mentioning that D and PF had negative impacts on the amount of PPV, while ST, MC and BS had positive impacts. In RSA, multiple regression analyses were conducted by varying one of the variable ranges (in this paper from −50% to +50%), while keeping the other variables constant [97]. It is noteworthy that these same processes were conducted for all variables to determine the influence of each variable on the output. The results of RSA indicated that the PPV result is highly sensitive to changes in D followed by ST and MC (see Figure 10). On the other hand, changes in the amount of BS have the lowest influence on the results. Results of the sensitivity analysis (especially regarding D, which is the most effective parameter on PPV) are in good agreement with studies carried out by [7,8,60]. For example, Monjezi et al. [60] proposed that the distance from the blast face is the most influential parameter on PPV among other used parameters (i.e., hole depth, stemming, maximum charge per delay). In addition, as discussed earlier, the distance from the blast face is one of the most effective parameters in empirical equations [15,17].

Conclusions
The PPV resulting from blasting operations-a critical environmental issue-needs to be predicted accurately for any blasting operations. The present research was centred upon the simulation of PPV in four quarry sites located in Malaysia. It took into account the input variables of stemming, powder factor, maximum charge per delay, burden-to-spacing ratio and the distance from the blast face. In addition, GEP was applied to the development of a predictive model, and the MC simulation was employed for the purpose of simulating all ranges of PPV results. All probable combinations of correlations and values of the variables were examined in the MC simulation process, and they were also analysed by carrying out 10,000 trials.

Conclusions
The PPV resulting from blasting operations-a critical environmental issue-needs to be predicted accurately for any blasting operations. The present research was centred upon the simulation of PPV in four quarry sites located in Malaysia. It took into account the input variables of stemming, powder factor, maximum charge per delay, burden-to-spacing ratio and the distance from the blast face. In addition, GEP was applied to the development of a predictive model, and the MC simulation was employed for the purpose of simulating all ranges of PPV results. All probable combinations of correlations and values of the variables were examined in the MC simulation process, and they were also analysed by carrying out 10,000 trials.
The results of the performance indices achieved for the developed GEP equation (training datasets: RMSE = 4.2814 and R 2 = 0.7283, and testing datasets: RMSE = 4.0344 and R 2 = 0.6823) confirmed its high capacity in terms of predicting the PPV value and also its applicability as an input equation into the MC simulation technique.
The range of PPV was simulated around 35 mm/s, and it starts from 1.13 mm/s up to 34.58 mm/s. Moreover, based on the results of the PPV, it can be stated with 100% confidence that the PPV is less than 50.8 mm/s, which is mentioned as the allowable maximum amount of PPV for surface mine blasting. Notably, based on the results of the RSA, the MC and ST have a direct relation with PPV, while any increase in the values of D and PF would lead to a reduction in the value of PPV. Moreover, based on the outcome of RSA and CSA, it was concluded that the result of the PPV is highly sensitive to the amount of D. These results are in good agreement with previous studies. The ST and PF are found to be the second and third most influential factors for PPV. Although BS was considered as an effective factor on PPV, its changes up to 50% were almost ineffective on PPV values. It is notable that the proposed equation is not applicable for other cases; however, the research framework can be adopted for other sites with different variables, values and functions for the purposes of predicting and simulating PPV.