Advanced Prediction of Roadway Broken Rock Zone Based on a Novel Hybrid Soft Computing Model Using Gaussian Process and Particle Swarm Optimization

Abstract

A simple and accurate evaluation method of broken rock zone thickness (BRZT), which is usually used to describe the broken rock zone (BRZ), is meaningful, due to its ability to provide a reference for the roadway stability evaluation and support design. To create a relationship between various geological variables and the broken rock zone thickness (BRZT), the multiple linear regression (MLR), artificial neural network (ANN), Gaussian process (GP) and particle swarm optimization algorithm (PSO)-GP method were utilized, and the corresponding intelligence models were developed based on the database collected from various mines in China. Four variables including embedding depth (ED), drift span (DS), surrounding rock mass strength (RMS) and joint index (JI) were selected to train the intelligence model, while broken rock zone thickness (BRZT) is chosen as the output variable, and the k-fold cross-validation method was applied in the training process. After training, three validation metrics including variance account for (VAF), determination coefficient (R²) and root mean squared error (RMSE) were applied to describe the predictive performance of these developed models. After comparing performance based on a ranking method, the obtained results show that the PSO-GP model provides the best predictive performance in estimating broken rock zone thickness (BRZT). In addition, the sensitive effect of collected variables on broken rock zone thickness (BRZT) can be listed as JI, ED, DS and RMS, and JI was found to be the most sensitive factor.

1. Introduction

After the excavation of the tunnel, the original stress balance of the rock mass around the roadway was destroyed, which leads to the redistribution of original stress and has an impact on the surrounding underground structures including surrounding rock mass and surrounding filling bodies [1,2,3]. Meanwhile, the reduction of surrounding rock mass strength and the stress concentration will appear. The rock mass will be ruptured while the concentrated stress is greater than the surrounding rock mass strength, until the new stress balance created, and a broken rock zone (BRZ) will be developed during that process (see Figure 1). According to the study of Zhao and Wu [4] and Xu and Jing [5], the broken rock zone thickness (BRZT) can be used to describe the broken rock zone (BRZ), analyze the roadway surrounding rock stability and guide the roadway support design. Therefore, a simple and high-precision broken rock zone (BRZ) predictive method is important and meaningful.

For the accurate broken rock zone thickness (BRZT) determination, three methods including the empirical formula method, the on-site measurement method and the machine learning method were proposed.

For the empirical formula method, Yan [6] proposed an empirical formula for the thickness prediction based on the wave velocity of rock mass and rock. Shemyakin [7,8,9,10] proposed the discontinuous zone concept and the empirical formula of the discontinuous zone thickness. Dong [11] proposed a broken rock zone predictive formula based on the indoor test and filed experience. Although these formulas can be used to predict the broken rock zone thickness (BRZT), only a few factors are considered, a few engineering sites have been verified, and it is hard to be applied in the complex engineering environment.

For the on-site measurement method, the geological radar method, the borehole camera method and the ultrasonic method were considered in the broken rock zone (BRZ) measurement, and were applied in many successful cases, such as Tongsheng Tunnel [12] and Jiaojia Gold Mine [13]. These on-site measurement methods are practical but time- and labor-consuming.

For the machine learning method, several studies were conducted by worldwide researchers. For example, Zhou and Li [14] tested the feasibility of using a support vector machine to evaluate the broken rock zone thickness (BRZT) based on 132 datasets collected from many coal roadways. Zhao and Wu [4] also proposed a support vector machine (SVM) model to determine the broken rock zone thickness (BRZT), and the average relative error of that model is 2.66%, which means that the obtained results agree well with the actual BRZT values. Ma [15] combined the PSO algorithm and least square support vector machines (LSSVM), based on a 32 datasets database, and the developed PSO-LSSVM was verified to have the characteristic of fast convergence speed and high calculation accuracy. Zhu [16] utilized the artificial neural network (ANN) model to predict the broken rock zone thickness (BRZT), and the relative error of 9% can powerfully prove that the developed ANN model has high precision. Peng [17] proposed an RBF neural network model based on Matlab 6.5, and the relative error of two testing datasets is respectively 0.4% and 0.92%. Jing [18] and Xu et al. [19] applied the adaptive neuro-fuzzy network to predict the broken rock zone thickness (BRZT) and achieved nice predictive performance. Xue [20] combined the genetic algorithm (GA) and backpropagation (BP) neural network to analyze an 18 datasets database and found that the predictive performance of the new proposed GA-BP is much better than the original BP model. Liu et al. [21] utilized a wavelet relevance vector machine in broken rock zone (BRZ) prediction, and the proposed model is showing a square correlation coefficient of 0.95.

In recent years, artificial intelligence (AI) technology has developed very rapidly [22,23]. Some new methods such as Gaussian process (GP) were proposed, and the performance of that method was verified by many successful applications [24,25,26,27,28]. Although some artificial intelligence models were applied in the roadways broken rock zone (BRZ) prediction, the application of the Gaussian process (GP) was not analyzed, the effect of various factors on the broken rock zone thickness (BRZT) was also not conducted. Hence, the Gaussian process (GP) method was selected to predict the broken rock zone thickness (BRZT) in this study, the PSO algorithm was utilized to search the optimal hyper-parameter combination for optimal predictive performance, and a sensitivity analysis was carried out to find the most sensitive factor in the present paper.

2. Theoretical Background

2.1. Gaussian Process (GP)

GP is a nonparametric model and consists of a set of random variables with the Gaussian distribution. For a GP model, it can be determined after selecting the mean function $h\left(x\right)$ and covariance function $c\left(x,x^{\prime }\right)$. So a GP model can be expressed by using the following equation [29,30].

$$f\left(x\right)~GP\left(h\left(x\right),c\left(x,x^{\prime }\right)\right)$$

(1)

When used for regression, GP will encrypt the uncertainty before executing the training. Then, Bayes rule will be used to update the beliefs through the function and calculate the posterior distribution [31]. To the authors’ best knowledge, the GP model has not been used in broken rock zone (BRZ) prediction, so a GP model associated with the rational quadratic covariance function was developed. For a high-performance GP development, it is necessary to select the optimal hyper-parameter combination, due to its huge impact on model performance. That combination can be obtained by using an optimization method or test-and-trial method, and the most widely used approach is to maximize the log marginal likelihood, and the conjugate gradient method is an effective method for these aims.

It should be noted that only a brief introduction was described here, more theory of GP can be found in the study of Rasmussen [29].

2.2. Particle Swarm Optimization Algorithm (PSO)

PSO is a popular meta-heuristic optimization algorithm, and the successful application of PSO can be found in many cases, including flyrock [32,33,34,35,36], air-overpressure [37] and rock strength [38,39]. In PSO, various particles will be created in a search space with randomly set location and velocity, then, the fitness value of each particle will be calculated by using the specified fitness function. In practice, particles in search space represent various solutions, and the main task of PSO is to find the optimal solution to the selected problem. For this aim, these particles move according to the velocity and location update equations with the guidance of the fitness value, and that procedure will be ended when the optimal solution was found [40].

To find the best position (g_best) of the particle, the velocity and position of particles were updated from the current position (X_c) and velocity value (V_c) to a new position (X_n) and velocity value (V_n), using the following equations [41,42].

$$\{\begin{array}{c}{V}_n=\omega V_c+c_1{p}_1\left(P_{best}-X_c\right)+c_2{p}_2\left(g_{best}-X_c\right)\\ X_n=X_c+V_n\hfill \end{array}$$

(2)

where the optimal position of a particle is expressed by using P_best. p₁, p₂, c₁ and c₂ are two random values and two positive acceleration constants.

2.3. PSO-GP Model Process

Using the PSO to optimize the GP model includes the following 5 steps (see Figure 2):

(1)

Randomly divide the database into training datasets and testing datasets with a ratio of 80% and 20% [43,44,45].

(2)

Scale the training datasets and testing datasets into the range of 0 to 1 for reducing the calculation difficulties.

(3)

Find the optimal hyper-parameter combination of the GP model.

(4)

Develop the optimal training GP model.

(5)

Check the predictive performance of the trained GP model using testing datasets and verification metrics.

2.4. Artificial Neural Network (ANN)

ANN has been utilized in the science and engineering area for more than 80 years, and its performance in creating a high-level and non-linear relationship for the selected variables have been verified in many cases [46,47,48,49,50,51,52]. ANN is inspired by the neural structure of the human brain, and various network architecture, learning algorithms and transfer functions were proposed after long-term development [53,54]. Among ANN, the multi-layer perceptron (MLP) was considered to be the best ANN architecture and contains five parts, including the input layer, the hidden layer, the output layer, neurons in the layer and weight among these neurons. The back-propagation (BP) algorithm is the most popular MLP, and the knowledge form the training datasets will be forwarded, while the errors will be passed back to guide the updating of weights value [55,56,57]. During BP training, the weight updated according to the cost function for achieving better training performance, and that process will be ended when the cost function equal to or less than the setting level [58].

2.5. Multiple Linear Regression (MLR)

MLR is a simple method for creating a linear relationship between the independent variables and dependent variables, and that method can be carried out by many software such as Matlab, SPSS and Python [59]. In many published papers, MLR was often utilized as a robust regression method to summarize rules and predict engineering phenomena. Generally, the empirical formula obtained by using MLR can be expressed as [60,61]:

$$Y=a_1{X}_1+a_2{X}_2+a_3{X}_3+\cdots +a_n{X}_n+b$$

(3)

In Equation (3), a and b represent the obtained constant of the developed linear relationship, X and Y are respectively the selected input and output variables.

2.6. K-Fold Cross-Validation and Verification Metrics

K-fold cross-validation is an obvious method for generalization performance evaluation. Generally, the k-fold cross-validation is inserted in the model training process, the results of k-fold cross-validation can guide the hyper-parameters updating [62,63]. After reviewing, 10-fold cross-validation was accepted by most researchers, and considered to be a suitable model verification method [64]. For carrying out that procedure, the training dataset will be divided into 10 parts, and 9 parts will be trained with current hyper-parameters. After this, the developed predictive model will be checked using the remaining part, and each dataset has the opportunity to be tested during the 10-fold cross-validation process [65]. After 10 runs, the average performance, which is normally evaluated by RMSE, will be used to represent the generalization performance of the trained model.

After model training, the trained model should be evaluated to check the applicability of unknown data prediction by using testing datasets, so some verification metrics should be selected. In this study, VAF (variance account for), R² (determination coefficient) and RMSE (root mean squared error) were employed to evaluate the relationship between predicted BRZT (BRZT^PRE) and actual BRZT (BRZT^ACT) [66,67,68,69,70,71].

$$VAF=\left[1-\frac{\mathrm{var}\left(BRZ{T}_i^{ACT}-BRZ{T}_i^{PRE}\right)}{\mathrm{var}\left(BRZ{T}_i^{ACT}\right)}\right]\times 100$$

(4)

$$R^2=\frac{{\left(n{\displaystyle \sum _{i=1}^{n}BRZ{T}_i^{PRE}BRZ{T}_i^{ACT}}-{\displaystyle \sum _{i=1}^{n}BRZ{T}_i^{PRE}}{\displaystyle \sum _{i=1}^{n}BRZ{T}_i^{ACT}}\right)}^2}{\left(n{\displaystyle \sum _{i=1}^n{\left(BRZ{T}_i^{PRE}\right)}^2}-{\left({\displaystyle \sum _{i=1}^{n}BRZ{T}_i^{PRE}}\right)}^2\right)}\cdot \frac{1}{\left(n{\displaystyle \sum _{i=1}^n{\left(BRZ{T}_i^{ACT}\right)}^2}-{\left({\displaystyle \sum _{i=1}^{n}BRZ{T}_i^{ACT}}\right)}^2\right)}$$

(5)

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(BRZ{T}_i^{PRE}-BRZ{T}_i^{ACT}\right)}^2}}$$

(6)

Here, n means the number of checked datasets, the best values of VAF, R², and RMSE are respectively 100, 1 and 0.

3. Application

3.1. Broken Rock Zone Database

After reviewing the previous literature [14,18,72,73], a roadways broken rock zone (BRZ) database with 181 datasets consists of embedding depth (ED), drift span (DS), surrounding rock mass strength (RMS), joint index (JI) and broken rock zone thickness (BRZT) was collected. For the deamination of JI values, the overall, block, layered, fragmented and loose rock mass were assigned a value of 1 to 5, respectively. For the determination of BRZT values, some mines in this database use the ultrasonic detection method, and some mines use the geophysical detecting radar method. Among these 181 datasets, there are 154 datasets were collected from the roadway with depth of less than 800 m. Considering that there is a zoning rupture phenomenon when the excavation works enter deep mining, the concept of broken rock zone (BRZ) cannot fully describe the rapture of rock around the roadway [4]. In China, most scientists think that deep than 800 m can be judged as deep mining [74,75], so only the datasets (see Figure 3) with depth less than 800 m were analyzed in this paper. In Figure 3, the asterisks means the significance levels of the correlations, the more asterisks, the smaller the p value.

As shown, the lower and upper triangle respectively represent the pairwise relationships and correlation coefficients, and the marginal distribution was shown by using the diagonal of Figure 3. It can be seen that there is an important relationship between these input variables. For better understanding, the maximum, minimum, standard deviation values of these five collected variables were tabulated by using

Table 1. Descriptive statistical values of the five collected variables in this paper.

Variable	Maximum	Minimum	Mean	Median	Standard Deviation	Skew	Kurtosis
ED	780.00	97.00	455.78	424.00	180.14	0.08	−1.19
DS	10.00	2.40	3.78	3.40	1.29	2.93	10.59
RMS	110.20	1.10	27.03	16.80	21.76	1.71	2.78
JI	5.00	1.00	3.25	3.00	0.92	0.07	−0.39
BRZT	3.45	0.30	1.43	1.40	0.56	0.63	0.71

.

3.2. Artificial Intelligence Model Results

3.2.1. Multiple Linear Regression Model

For broken rock zone prediction, the MLR method was utilized and the following empirical formula was created after using SPSS software.

$$BRT=0.001\times ED+0.113\times DS-0.003\times RS+0.384\times JI-0.697$$

(7)

By using the above formula, the predicted BRZT can be calculated and the predicting results with VAF of 76.2155, R² of 0.7639 and RMSE of 0.2973 for training datasets and VAF of 88.6696, R² of 0.8899 and RMSE of 0.1678 for testing datasets can be found in Figure 4.

3.2.2. Artificial Neural Network Results

For an accurate ANN model development, the structure of the hidden layer including the number of hidden layers and nodes in hidden layers should be determined first. Given the study of Mohamad et al. [38] and Hasanipahah et al. [40], the application of only a hidden layer also has a nice performance in solving the regression problems. Meanwhile, studies of Nguyen [76,77] also show that the ANN model with two layers is often used for more complex problems. Hence, a parameter investigation was carried out to find out the optimal number of nodes and the optimal number of hidden layers for BRZ prediction. In this investigation, various ANN models with the number of nodes from 1 to 9 and the number of hidden layers for 1 to 3 were developed, the 10-fold cross-validation method was applied, each model was run five times, and the optimal ANN models with different structure will be compared. After comparing, 2 was found to be the optimal number of hidden layers, 7 and 4 were respectively found to be the optimal number of nodes in the first hidden layer and the second hidden layer. Then, the predictive performance with VAF of 90.7759, R² of 0.9086 and RMSE of 0.1779 for training datasets and VAF of 88.4466, R² of 0.8993 and RMSE of 0.1630 for testing datasets was achieved (see Figure 5).

3.2.3. Gaussian Process Results

Using the Gaussian process method, the hyper-parameters were searched, and the prepared training datasets were trained. After analyzing, the results show that the Gaussian process model can provide the predictive performance with VAF of 90.4209, R² of 0.9051 and RMSE of 0.1798 for training datasets and VAF of 86.7895, R² of 0.8900 and RMSE of 0.1733 for testing datasets, and the actual BRZT and predicted BRZT were plotted in Figure 6.

3.2.4. PSO-GP Results

Before using the PSO to find the hyper-parameter of the GP model for achieving the best predictive performance, several parameters including the number of particles and the maximum number of iterations of PSO should be defined by the user. In the study of Fang et al. [78], Koopialipoor et al. [79] and Nguyen et al. [80], it can be found that the optimization algorithm such as GA, PSO and ABC can usually get stable results with the maximum number of iterations of 1000. Therefore, a parameter investigation for 9 PSO-GP models with the various number of particles from 20 to 300, the maximum number of iterations of 1000 and a search range of −10 to 10 was carried out, and it was found that the PSO model marked by 300 particles and 1000 iterations yield the lowest fitness curve. Moreover, no change of the fitness curve can be found after 200 iterations, so 300 and 200 were determined to be optimal parameters combination of the PSO-GP model. After calculation, the optimal PSO-GP model is showing the VAF of 94.3008, R² of 0.9438 and RMSE of 0.1387 for training datasets, and the VAF of 89.2192, R² of 0.9153 and RMSE of 0.1591 for testing datasets, and the actual BRZT and predicted BRZT were shown in Figure 7.

3.3. Comparisons and Discussions

To compare the predictive performance of developed four predictive models and verify the optimization capability of PSO, the predictive performance of these models was compared by using Figure 8 and a ranking method provided by Zorlu [81] (see Figure 9).

As shown in Figure 8, predicted results provided by PSO-GP are most consistent with the actual measurement results, and the best predictive performance was also proved by the highest total ranking values of 24 by using Zorlu’s method.

For a better understanding of roadways broken rock zone (BRZ), a sensitivity analysis using cosine amplitude method [82] was carried out here to determine the influenced factors of these input variables on broken rock thickness. The sensitivity factor was obtained by using the following equations:

$$ED=\left\{97,125,125,\cdots ,780\right\}$$

(8)

$$DS=\left\{2.6,2.8,2.8,\cdots ,3.0\right\}$$

(9)

$$RS=\left\{11.2,13.3,13.3,\cdots ,70.6\right\}$$

(10)

$$JI=\left\{3,3,2,\cdots ,1\right\}$$

(11)

$$BRT=\left\{1.20,1.30,0.70,\cdots ,0.65\right\}$$

(12)

$$S_{ED\text{-}BRZT}=\frac{{\displaystyle {\sum }_{i=1}^{154}E{D}_i\cdot BRZ{T}_i}}{\sqrt{{\displaystyle {\sum }_{i=1}^{154}E{D}_i^2}}\cdot \sqrt{{\displaystyle {\sum }_{i=1}^{154}BRZ{T}_i^2}}}$$

(13)

$$S_{DS\text{-}BRZT}=\frac{{\displaystyle {\sum }_{i=1}^{154}D{S}_i\cdot BRZ{T}_i}}{\sqrt{{\displaystyle {\sum }_{i=1}^{154}D{S}_i^2}}\cdot \sqrt{{\displaystyle {\sum }_{i=1}^{154}BRZ{T}_i^2}}}$$

(14)

$$S_{RS\text{-}BRZT}=\frac{{\displaystyle {\sum }_{i=1}^{154}R{S}_i\cdot BRZ{T}_i}}{\sqrt{{\displaystyle {\sum }_{i=1}^{154}R{S}_i^2}}\cdot \sqrt{{\displaystyle {\sum }_{i=1}^{154}BRZ{T}_i^2}}}$$

(15)

$$S_{JI\text{-}BRZT}=\frac{{\displaystyle {\sum }_{i=1}^{154}J{I}_i\cdot BRZ{T}_i}}{\sqrt{{\displaystyle {\sum }_{i=1}^{154}J{I}_i^2}}\cdot \sqrt{{\displaystyle {\sum }_{i=1}^{154}BRZ{T}_i^2}}}$$

(16)

After calculating, the influenced factors of ED, DS, RMS and JI to BRZT are, respectively, 0.9384, 0.9162, 0.6651 and 0.9711, and it can be concluded that the joint index (JI) has the most influence on broken rock zone thickness (BRZT) prediction under the collected database.

In this study, four predictive models for roadway the broken rock zone (BRZ) prediction, including MLR model, ANN model, GP model and PSO-GP model were proposed and can provide the nice predictive performance of that engineering problem. However, some limitations can be found as follows:

• Due to the difference in mechanical behavior between shallow and deep mining, this paper focuses on shallow mining and can just be used in the prediction of roadway broken rock zone (BRZ) with depth less than 800 m.
• A large-scale and comprehensive comparison with various artificial intelligence models for a specified engineering problem and the same database is meaningful and can provide a guide for the on-site application. Although MLR, ANN, GP and PSO-GP were utilized in this paper, other methods, such as random forest (RF), support vector machine (SVM), extreme learning machine (ELM) and gene expression programming (GEP) were not analyzed, and other optimization algorithms, such as artificial bee colony algorithm (ABC), Salp Swarm algorithm (SSA), Harris Hawks optimization algorithm (HHO) and cuckoo search algorithm (CS), were not utilized.
• Although it is now possible to use some methods, such as borehole camera, geophysical detecting radar and ultrasonic detection, to measure the broken rock zone, it is still hard to develop a big enough database with enough datasets and enough factors. Though the authors tried the best to collect datasets, only 181 datasets were found, and that database is far from the big data level. Some studies [83,84] show that the broken rock zone is quite complex issue and is affected by many factors, such as rock beds gradient, drivage in the fold, seismic activity, goafs, drivage with blasting, etc. A qualified database should contain as many mine data as possible, due to the different conditions of each mine, we failed for developing a big enough database with various factors, so only four factors were utilized in this paper. In addition, the measurement method of broken rock zone thickness for some mines are different, for example, the BRZT was determined by using ultrasonic detection in Baizigou Coal Mine, Shuanggou Coal Mine and Jianchang Coal Mine, while it was determined by using geophysical detecting radar in Huafeng Coal Mine. Although these methods were utilized in many mines for BRZT measurement, the BRZT values in the collected database may not maintain uniform precision, due to the influence of the measurement method and measuring instruments.

As a data-driven method, the proposed models in this study are developed based on the collected BRZT database, so these models may provide nice performance when predicting BRZT of the mines contained in this database, but the performance cannot be guaranteed for the mines that have never been trained. Meanwhile, the limitations analyzed above may affect the analysis results, including intelligence predictive performance and sensitivity analysis results, and these limitations should be considered in the future investigation.

Nevertheless, this study is a powerful supplement to the roadway broken rock zone (BRZ) prediction and is a pioneer work, due to the application of GP and PSO-GP having not been introduced for that engineering problem.

4. Conclusions

For the simple and accurate prediction of roadways broken rock zone (BRZ), a Gaussian process model optimized by PSO based on a 181 datasets database was proposed. During the training process, 10-k-fold cross-validation was utilized for evaluating the generalization performance and guide the model training, three validation metrics including VAF, R² and RMSE were used to verify the predictive performance of developed intelligence models. After comparing with the prediction performance of MLR (VAF of 76.2155, R² of 0.7639 and RMSE of 0.2973 for training datasets and VAF of 88.6696, R² of 0.8899 and RMSE of 0.1678 for testing datasets), ANN (VAF of 90.7759, R² of 0.9086 and RMSE of 0.1779 for training datasets and VAF of 88.4466, R² of 0.8993 and RMSE of 0.1630 for testing datasets) and GP (VAF of 90.4209, R² of 0.9051 and RMSE of 0.1798 for training datasets and VAF of 86.7895, R² of 0.8900 and RMSE of 0.1733 for testing datasets), the obtained results show that the PSO-GP shows the best performance with VAF of 94.3008, R² of 0.9438 and RMSE of 0.1387 for training datasets and the VAF of 89.2192, R² of 0.9153 and RMSE of 0.1591 for testing datasets. Meanwhile, the obvious optimization performance of the PSO algorithm was also verified. Last, the sensitivity analysis shows that there is no big difference between the effect of ED, RMS and JI on BRZT, and JI was found to be the most sensitive factor. The developed predictive model in this paper provides a simple and accurate method for broken rock zone thickness (BRZT) prediction, and can be referenced when solving similar complex engineering problems.