Assessment of the Rock Elasticity Modulus Using Four Hybrid RF Models: A Combination of Data-Driven and Soft Techniques

Abstract

The determination of the rock elasticity modulus (EM) is an indispensable key step for the design of rock engineering problems. Traditional experimental analysis can accurately measure the rock EM, but it requires manpower and material resources, and it is time consuming. The EM estimation of new rocks using former published empirical formulas is also a possibility but can be attached of high uncertainties. In this paper, four types of metaheuristic optimization algorithms (MOA), named the backtracking search optimization algorithm (BSA), multi-verse optimizer (MVO), golden eagle optimizer (GEO) and poor and rich optimization algorithm (PRO), were utilized to optimize the random forest (RF) model for predicting the rock EM. A data-driven technology was used to generate an integrated database consisting of 120 rock samples from the literature. To verify the predictive performance of the proposed models, five common machine-learning models and one empirical formula were also developed to predict the rock EM. Four popular performance indices, including the root-mean-square error (RMSE), mean absolute error (MAE), the coefficient of determination (R²) and Willmott’s index (WI), were adopted to evaluate all models. The results showed that the PRO-RF model has obtained the most satisfactory prediction accuracy. The porosity (Pn) is the most important variable for predicting the rock EM based on the sensitive analysis. This paper compares the performance of the RF models optimized by using four MOA for the rock EM prediction. It provides a good example for the subsequent application of soft techniques on the EM and other important rock parameter estimations.

1. Introduction

In rock engineering, the rock elasticity modulus (EM) plays an important role for structure designs [1,2,3,4]. The EM is an important index for quantifying the rock behavior. It is also closely related to the rock’s durability, which can determine rock applications to a large extent [5]. Numerous experiments were developed according to the international society of rock mechanics (ISRM) to calculate the rock EM [6]. Nevertheless, the expensive sample costs and time-consuming laboratory operations have forced engineers to develop other methods for estimating the rock EM.

The empirical formula based on statistics is a popular method to estimate the rock EM in preliminary design phases [7,8,9,10,11,12]. The aim of the empirical formulas (simple-regression (SR) or multiple-regression (MR) formulas) is to establish a relationship between one or more rock properties and the EM. Numerous researchers proposed various empirical formulas to estimate the rock EM [3,4,13,14,15,16]. Beiki et al. [17] used the porosity (Pn) to predict the rock EM with a low prediction accuracy of the coefficient of determination (R²). Yasar and Erdogan [18] provided a linear SR formular using P-wave velocity (Vp) to estimate EM. Dincer et al. [19] established an SR formula between the Schmidt hammer rebound number (SHRN) and EM. The prediction accuracy of R² was equal to 0.85. Behzadafshar et al. [20] only utilized the point load index (PLI) to forecast the granite EM. The prediction accuracy was not satisfactory, resulting in a low R² of 0.58. Dehghan et al. [21] developed an MR formula consisting of the above four variables to predict the travertine EM. By reviewing the developed empirical formula, the Pn, Vp, SHRN and PLI are usually considered into these empirical equations, which can be obtained directly from tests [13]. However, the original data used to develop each empirical formula are fixed and different, and the formula performance decreases once it is used to predict the EM of a new rock [22]. From the prediction accuracy perspective, the intrinsic relationship of the SR and MR formulas cannot well describe the complex and nonlinear correlation between the rock properties and the EM [23].

In recent years, soft techniques have been widely used to solve prediction problems in rock engineering [24,25,26,27,28,29,30], especially machine-learning (ML) methods. For the rock EM prediction, Ocak and Seker [31] proposed an artificial neural network (ANN) model to predict the EM of intact rock. The results showed that the ANN model obtained a higher accuracy for the root-mean-square error (RMSE is equal to 0.191) than those of previous approaches. Pappalardo and Mineo [32] utilized ANN models to estimate the EM statis of rock samples. The results indicated that this artificial intelligence model is of more practical value in estimating the rock EM. Singh et al. [33] used an adaptive neurofuzzy inference system (ANFIS) to predict the EM values of 85 rock samples. The performance evaluation results illustrated that the ANFIS model has a better performance than the initial ANN model and fuzzy inference system (FIS) model by means of the lower RMSE value of 6.799. Umrao et al. [34] developed an ANFIS to estimate the EM of 45 heterogenous sedimentary rocks. This model has obtained a satisfactory prediction accuracy of R² = 0.935. Acar and Kaya [35] adopted a least-square support vector machine (LS-SVM) model to find the EM of weak rocks by considering the Vp, unit weight (r), PLI and tensile strength (Ts). The prediction results showed that the LS-SVM models can be a good substitute for experiments to measure the weak rock EM. Al-Anazi and Gates [36] used a support vector regression (SVR) model and a backpropagation neural network (BPNN) model to forecast the EM of reservoir rocks. The prediction results indicated that the former has a better performance than the latter. Matin et al. [37] used an integrated ML model, named the random forest (RF) model, to predict the EM and the uniaxial compressive strength (UCS) of various rocks. Based on the prediction results, the RF model not only achieves a satisfactory accuracy (R² = 0.91) for the EM prediction but can also accurately estimate the UCS (R² = 0.93). Other similar studies on the rock EM prediction by ML models can be referred to in the literature [38,39,40,41,42]. To improve the ML model’s performance, the metaheuristic optimization algorithms (MOA) are used to select the model hyperparameters. Tian et al. [43] used the imperialism competitive algorithm (ICA) and the particle-swarm optimization (PSO) to optimize an ANN model for predicting the EM of rock materials. The optimization results illustrated that the ICA-ANN has the best prediction accuracy for both the training and testing phases (R²: 0.952 and 0.955). Mokhtari and Behnia [44] combined the cuckoo optimization algorithm (COA) and ANN model to estimate the EM of limestone rocks. The results showed that the COA can obviously improve the prediction performance of the ANN model. Other optimized ML models for the EM prediction can be found in the literature [45,46,47,48,49,50].

Among common ML models, the RF model has unique advantages in resisting the overfitting phenomenon, and its combination with MOA can effectively solve the hyperparameter selection problem [51]. In general, MOA can be divided into four groups, i.e., the based evolutionary, based swarm intelligence, based human behavior and based physicochemical groups. Therefore, the aim of this paper is to generate four optimized RF models using four different MOA strategies for predicting the rock EM. An integrated rock database was established using the data-driven technology to train and test the proposed models. In addition, five widely used ML models and one empirical equation were also developed to compare the predictive performance with four hybrid RF models for predicting the rock EM. Finally, a sensitive analysis is carried out to calculate the variable importance.

2. Data Preparation

Reviewing the published studies on the rock EM prediction, the used rock samples are generally different, resulting in a loss of accuracy in predicting other rock sample properties from original empirical formulas or other models. To overpass this limitation, an integrated rock database consisting of 120 rock samples (e.g., granite and travertine) was established using the data-driven technology. The EM values of these rock samples was investigated by Armaghani et al. [3], Dehghan et al. [21] and Tuğrul and Zarif [52]. The reason for choosing these data is that the EM value of each rock is determined jointly by the porosity (Pn), P-wave velocity (Vp), Schmidt hammer rebound number (SHRN) and point load index (PLI). Before generating the prediction models, a correlation analysis needs to be conducted by outputting the correlation coefficient (CC) to determine the final variables used in the EM prediction. As illustrated in Figure 1, the scatter plot at the lower right shows the distribution of the four input variables and of one output variable. The histograms in the diagonal showed the data range of all variables, and the number at the upper left represents the CC between the two corresponding variables. If the number value is negative, the correlation between the two variables is negative. Otherwise, the correlation is positive. The absolute number value is used to evaluate the correlation between any of two variables. If the CC between two input parameters (or between each input and output parameter) is very high (or very low), one of them needs to be removed to increase the prediction efficiency. The results showed that the CC between the four variables is not high; Pn and EM have the highest CC value (−0.651). Especially, the low correlation between PLI and EM is caused by the diverse sources of databases used to predict EM in this paper, while the PLI is beneficial to accurately predict the rock EM [3]. Therefore, the Pn, Vp, SHRN and PLI are used as input variables to predict the EM (output variable). Their detailed information is shown in

Table 1. Detailed description of the input and output variables.

Variables	Statistical Information
	Sign	Unit	Min	Max	Mean	St. D
Point load index	PLI	MPa	0.890	12.530	4.365	2.839
Porosity	Pn	%	0.100	10.270	1.957	3.047
P-wave velocity	Vp	km/s	2.823	7.943	5.575	0.892
Schmidt hammer rebound number	SHRN	/	25.630	72.000	47.093	13.795
Elasticity modulus	EM	GPa	3.050	183.300	60.139	44.832

Note: Min—minimum value; Max—maximum value; St. D—standard deviation value.

.

3. Development of Hybrid RF Models for Predicting the Rock EM

3.1. Metaheuristic Optimization Algorithms

3.1.1. Backtracking Search Optimization Algorithm (BSA)

Civicioglu [53] proposed a MOA-based evolutionary method to solve optimization problems, namely, the backtracking search optimization algorithm (BSA). The algorithm is still inspired by an individual evolution, especially the traditional mutation, and crossover operators are uniformly replaced by breeding operators. The optimization process of the BSA can be described using five strategies: initialization, selection, mutation, crossover and selection. These strategies are defined by using Equations (1)–(5).

$$P_i=Low+rand(N)\cdot (Up-Low)$$

(1)

$$P_{oi}=\left\{\begin{array}{l}P\left[a,b~P_i\right), a<b\\ permuting(P_{oi}), a>b\end{array}\right.$$

(2)

$$Mutant=P_i+F\left(P_{oi}-P_i\right)$$

(3)

$$T_i=P_i+map\cdot F\left(P_{oi}-P_i\right)$$

(4)

$$P_{ui}=T_i, fitness(T_i)<fitness(P_i)$$

(5)

where $P_i$ and $P_{oi}$ represent the initial and historical position of the i-th individual, respectively. Low and Up indicate the lower and upper bounds of the dimension space. a and b are random numbers within the range of [0, 1]. F and map represent the control parameter of the mutation operator and a binary matrix consisting of crossover probability parameters, respectively. $T_i$ and $P_{ui}$ are the current position of the i-th tested individual and the updated position of the i-th individual, respectively.

3.1.2. Multi-Verse Optimizer (MVO)

The multi-verse optimizer (MVO) is a MOA-based physics algorithm proposed by Mirjalili et al. [54], which is inspired by the idea that the universe moves from white holes to black holes through wormholes to achieve a stable situation. The white holes are believed to be an important part of the original universe, the black holes have an irresistible pull on everything including light beams, and the wormhole is a bridge or passage connecting different universes. In the MVO algorithm, the birth of the universe is always related to the objects transfer. The universe with a low inflation rate is more likely to take in more objects. The optimization process of MVO can be described as follows:

(1)

Population—the initial population of the universes in the searching space is defined using the Equation (6).

$$U=\left[\begin{array}{cccc}{u}_1^1& u_1^2& \cdots & u_1^d\\ u_2^1& u_2^2& \cdots & u_2^d\\ ⋮& ⋮& ⋮& ⋮\\ u_n^1& u_n^2& \cdots & u_n^d\end{array}\right]$$

(6)

where $u_n^d$ indicates the parameter of the n-th universe in the d-dimension searching space.

(2)

Exploration and exploitation—the function of wormholes is to help objects move from one universe to another (see Figure 2). Thus, this mechanism by which objects are exchanged between universes through wormholes can be described as:

$$u_j=\left\{\begin{array}{l}{u}_j+TDR\cdot (\left(ub-lb\right)\cdot r_3+lb), r_2<0.5 \mathrm{and} r_1<WEP\\ u_j+TDR\cdot (\left(ub-lb\right)\cdot r_3+lb), r_2>0.5 \mathrm{and} r_1<WEP\\ u_j, r_1>WEP\end{array}\right.$$

(7)

where $u_j$ represents the j-th parameter of the best universe. ub and lb are the lower and upper bounds of the multi-universes space, respectively. $r_1$, $r_2$ and $r_3$ indicate three random numbers within the range of [0, 1]. TDR and WEP represent two coefficients, the former is the wormhole existence probability, and the latter is named the travelling-distance rate. These coefficients can be calculated using the Equations (8) and (9).

$$WEP=\min +t\cdot (\frac{\max -\min }{T})$$

(8)

$$TDR=1-\frac{t^{1/e}}{T^{1/e}}$$

(9)

where min and max represent the minimum maximum values, respectively. t and T indicate the current iteration and the maximum iteration, respectively. e is the exploitation accuracy.

3.1.3. Golden Eagle Optimizer (GEO)

Mohammadi-Balani et al. [55] developed a novel MOA-based swarm intelligence named the golden eagle optimizer (GEO) to provide an effective scheme for solving optimization problems. This algorithm is inspired by the hunting behavior of golden eagles, who can adjust their speed to hunt. The hunting behavior can be divided into three parts: (a) selecting the prey; (b) attacking the prey; and (c) cruising. It is worth noting that each hunting behavior of a golden eagle is carefully considered. It allows striking a balance between attack and cruise. Once an attack is launched, the golden eagle is unable to obtain food or replenish enough energy.

(I)

Selecting the prey—the selection can occur in a basic way, with each golden eagle randomly select a prey from the memory of any other group member to better explore the landscape. It is important to note that the chosen prey is not necessarily the nearest or furthest prey. Figure 3 shows how prey selection works.

(II)

Exploration and exploitation—after determining the prey, each golden eagle carries out the attacking and cruising behaviors. The attacking behavior can be expressed by the following mathematical formula:

$$Y_i^G=Y_l^P-A_i$$

(10)

where $Y_i^G$ and $Y_l^P$ represent the position of the i-th golden eagle and the prey determined by the l-th golden eagle, respectively. $A_i$ indicates the attacking distance between the prey and the i-th golden eagle. The cruising behavior is related to the attacking behavior, which can be expressed using the Equations (11) and (12).

$$h_1{y}_1+h_2{y}_2+\dots +h_m{y}_m=d\Rightarrow {\displaystyle \sum _{z=1}^m{h}_z{R}_z}$$

(11)

$${\displaystyle \sum _{z=1}^m{\alpha }_z{y}_z}={\displaystyle \sum _{z=1}^m{\alpha }_z^s{y}_z^{*}}$$

(12)

where $\left[h_1,h_2,\dots h_m\right]$, $\left[y_1,y_2,\dots y_m\right]$ and $\left[R_1,R_2,\dots R_m\right]$ represent the normal coefficients, variables and random points, respectively. s is the current iteration. $\left[y_1^{*},y_2^{*},\dots y_m^{*}\right]$ indicates the position of the selected prey. $\left[{\alpha }_1,{\alpha }_2,\dots {\alpha }_m\right]$ belongs to the $A_i$.

3.1.4. Poor and Rich Optimization Algorithm (PRO)

The poor and rich optimization algorithm (PRO) was proposed by Moosavi and Bardsiri [56]. It is inspired by people’s desires and attitudes towards money. The wealth accumulated by individuals can preliminarily divide society into two classes, namely, rich and poor. For the rich (i.e., wealth level is obviously higher than the average one), observing the behavior of the poor (i.e., wealth level is obviously lower than the average one) can help them to increase wealth and consolidate their class position. The poor tend to narrow the gap by learning from rich ideas about wealth and approaches to making money. In the PRO algorithm, the population distribution of rich and poor can be expressed using Equation (13).

$$PO{P}_{main}=PO{P}_{poor}+PO{P}_{rich}$$

(13)

where $PO{P}_{main}$ represents the main population, which is related to the $PO{P}_{poor}$ (poor population) and $PO{P}_{rich}$ (rich population). It should be noted that the position of the rich is better than the poor position. Their positions are calculated using the Equations (14) and (15), respectively.

$$P_{rich}^{*}=P_{rich}+c\cdot \left(P_{rich}-P_{poor}^{best}\right)$$

(14)

$$P_{poor}^{*}=P_{poor}+\left[c(\mathrm{Pattern})-P_{poor}\right]$$

(15)

where $P_{rich}^{*}$ and $P_{rich}$ represent the updated and old positions of the rich, respectively. $P_{poor}^{*}$ and $P_{poor}$ indicate the updated and old positions of the poor, respectively. $P_{poor}^{best}$ is the current position of the best people in the poor population. c is a random number within the range of 0 to 1. The Pattern value is calculated by using the Equation (16).

$$\mathrm{Pattern}=\frac{P_{rich}^{best}+P_{rich}^{mean}+P_{rich}^{worst}}{3}$$

(16)

where $P_{rich}^{best}$ and $P_{rich}^{worst}$ represent the current positions of the best and worst people in the rich population, respectively. $P_{rich}^{mean}$ indicates the average position of the people in the rich population. The position distribution of rich and poor is shown in Figure 4.

However, some sudden changes can occur for the rich and poor situations, such as the stock prices falling and rising, shortages of petroleum products and inflation [56]. Since the occurrence of the above situation is not predictable, $P_{rich}^{*}$ and $P_{poor}^{*}$ after mutation can be expressed as follows:

$$P_{rich}^{*}=P_{rich}^{*}+randn,if rand<Pmut$$

(17)

$$P_{poor}^{*}=P_{poor}^{*}+randn,if rand<Pmut$$

(18)

where rand is a random number within the range of 0 to 1. randn represents a value considering a normal distribution. $Pmut$ represents the mutation probability.

3.2. Hybrid RF Models

In this paper, the BSA, MVO, GEO and PRO are used to optimize the RF model for predicting the rock EM. The definition of the RF model is described in the literature [57,58,59,60]. The hyperparameter combination selection (i.e., number of trees (Nt) and random features (Maxdepth)) is a key step to tap into the RF model prediction potential. Therefore, four hybrid RF models were generated to find the optimal hyperparameter combination according to the following process:

(i)

Data preprocessing

A total of 120 rock samples with four input variables were used to predict the EM in this paper. All variables need to be extracted and normalized to [−1, 1]. The purpose of this step is to prevent a failure for establishing the accurate prediction relationships due to the parameter variability. After that, the train and test sets are separated from the initial database. The ratio of the train set to the test set is set equal to 4 to 1. It should be noted that the same train or test set is used to generate each hybrid RF model for predicting the rock EM and comparing their performance.

(ii)

Parameter settings

Although the Nt increase will not cause an overfitting of the RF model, a large parameter selection range can greatly increase the computation time. Therefore, the ranges of Nt and Maxdepth are set equal to [1, 100] and [1, 10], respectively. For the four MOA algorithms, the number of initial solutions (i.e., individuals of BSA, candidates of MVO, population of GEO and human of PRO) and the iteration time are the core factors that affect the optimization performance of these algorithms. To better activate the optimization performance, the solutions are set equal to 30, 60, 90, 120 and 150 during the 200 iterations.

(iii)

Optimization evaluation

The fitness function is utilized to evaluate the performance of each hybrid RF model with different solutions during the 200 iterations. The RMSE is adopted to represent the fitness values of all models in this paper. They do not need an absolute value to evaluate the model performance [51]. In other words, the best-optimized RF model has the lowest RMSE value among all hybrid models based on the same MOA. The flowchart for developing four hybrid RF models for predicting the rock EM is shown in Figure 5.

4. Performance Evaluation

The statistical indices widely utilized to evaluate the performance of the prediction models include the root-mean-square error (RMSE), the mean absolute error (MAE), the coefficient of determination (R²) and the Willmott’s index (WI). The RMSE and MAE are able to reflect the error between the predicted and the measured values. On the other hand, the R² and WI describe the fitting performance of the prediction models. In terms of values, the lowest values of RMSE and MAE and the highest values of R² and WI represent the best prediction model.

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(E_i-e_i\right)}^2}}$$

(19)

$${\mathrm{R}}^2=1-\frac{{\left[{\displaystyle \sum _{i=1}^n(E_i-e_i)}\right]}^2}{{\left[{\displaystyle \sum _{i=1}^n(E_i-\overline{E})}\right]}^2}$$

(20)

$$\mathrm{MAE}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|E_i-e_i\right|}$$

(21)

$$\mathrm{WI}=1-\left[\frac{{\displaystyle {\sum }_{i=1}^n{\left(E_i-e_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^n{\left(\left|e_i-\overline{E}\right|+\left|E_i-\overline{E}\right|\right)}^2}}\right]$$

(22)

where n is the number of the samples. $E_i$ and $e_i$ represent the measured and the predicted values of the rock EM, respectively. $\overline{E}$ is the average value of the measured rock EM.

5. Results and Discussion

5.1. Results of the Proposed Four Hybrid Models

To determine the optimal solution and corresponding hyperparameter combination (i.e., Nf and Maxdepth) of the RF model, all hybrid models were performed during 200 iterations. The iteration curves of each hybrid RF model with five solutions are shown in Figure 6. As it can be seen in Figure 6a, it is obvious that the BSA-RF with 60 solutions has a lower fitness value than the other four BSA-RF models during the 200 iterations. The solution of 90 is the most suitable for generating the MVO-RF model by means of the lowest fitness value (see Figure 6b). As illustrated in Figure 6c,d, the optimal solutions of the GEO-RF and PRO-RF models are equal to 120 and 90, respectively.

Table 2. Development results of all hybrid RF models.

Solutions	Fitness (RMSE)
	BSA-RF	MVO-RF	GEO-RF	PRO-RF
30	0.1941	0.1893	0.1987	0.1901
60	0.1868	0.1977	0.1974	0.1935
90	0.1947	0.1870	0.1934	0.1861
120	0.1940	0.1942	0.1925	0.1875
150	0.1927	0.1917	0.1940	0.1928
Optimal hyperparameter combination
Nf	19	21	20	17
MaxDepth	2	2	2	2

lists the results and the best hyperparameter combination of the four optimized RF models. The PRO-RF model has the lowest value of RMSE (0.1861) among all models. This optimized RF model shows that the optimal Nf and Maxdepth values are 17 and 2, respectively.

Four hybrid RF models with the optimal hyperparameter combinations were used to predict the rock EM in the training phase. The performance indices of each hybrid RF model are listed in

Table 3. Performance and ranking results of four hybrid RF models in the training phase.

Models	Performance Indices and Ranking Scores								Total
	R2	Score	RMSE	Score	MAE	Score	WI	Score
BSA-RF	0.9359	2	11.3203	2	7.8165	2	0.9824	2	8
MVO-RF	0.9407	3	10.8867	3	7.7471	3	0.9837	3	12
GEO-RF	0.9317	1	11.6807	1	8.0452	1	0.9809	1	4
PRO-RF	0.9423	4	10.7420	4	7.6514	4	0.9843	4	16

. As illustrated in this table, four models achieved a good prediction accuracy with high values of R² and WI and low values of RMSE and MAE. Compared with the other three models, the PRO-RF is the best prediction model for predicting the EM by means of the best performance indices, i.e., R² is equal to 0.9423, RMSE to 10.7420, MAE to 7.6514 and WI is equal to 0.9843. The ranking score results of all models indicated that the GEO-RF model with the lowest total score of 4 is the worst model among all hybrid models for predicting EM in the training phase. The score of the MVO-RF model is the second model in the rank just after the PRO-RF model.

As illustrated in Figure 7, the regression distribution of the predicted and the measured EM values represents the performance of the prediction models. The position of each rock EM data in the two-dimension regression diagram is determined by the predicted EM values (horizontal axis, x) and the measured EM values (vertical axis, y). If y = x, it means that the predicted EM value is equal to the measured EM value. The rock EM data are then located on the diagonal line. To this end, the PRO-RF model has obtained the most rock EM data on the diagonal or close to the line among the four hybrid RF models (Figure 7d). After the PRO-RF model, the MVO-RF model and the BSA-RF model have similar regression distributions for both the small and big rock EM data. The MVO-RF model has obtained more data close to the diagonal line than the BSA-RF model in the interval between 20 and 70.

To further determine the performance of all hybrid RF models, the performance indices were calculated again using the test set, as shown in

Table 4. Performance and ranking results of four hybrid RF models in the testing phase.

Models	Performance Indices and Ranking Scores								Total
	R2	Score	RMSE	Score	MAE	Score	WI	Score
BSA-RF	0.9322	3	10.8902	3	7.2155	3	0.9812	3	12
MVO-RF	0.9236	2	11.5567	2	7.5280	2	0.9785	2	8
GEO-RF	0.9123	1	12.3826	1	7.9739	1	0.9746	1	4
PRO-RF	0.9410	4	10.1548	4	6.0423	4	0.9840	4	16

. The best model is still the PRO-RF, which has the lowest values of RMSE and MAE (10.1548 and 6.0423) and the highest values of R² and WI (0.9410 and 0.9840). The predictive performance of the BSA-RF model is better than the MVO-RF model in the testing phase, the former not only has the better performance indices but also has a higher score (12) than the latter (8). In addition, the GEO-RF model still does not achieve better predictive performances than the other models using the test set.

Figure 8 illustrates the regression results of the four hybrid RF models in the testing phase. As it can be seen in Figure 8a–c, large errors between the predicted EM values by the BSA-RF, MVO-RF and GEO-RF models and the measured EM values indicate that the data point move away from the diagonal, especially when the EM values are in the range of [50, 80] and [0, 20]. The PRO-RF model has obtained the most EM data points close to the diagonal based on the best predictive performance in the testing phase. Therefore, the PRO-RF model is the best hybrid RF model for predicting the rock EM in this paper.

5.2. Performance Comparison between the Proposed and Other Models

To compare the predictive performance with the proposed hybrid models, five common ML models, named the ANN, SVR, extreme learning machine (ELM), kernel-extreme learning machine (KELM) and generalized regression neural network (GRNN), and one empirical formula, were also developed to predict the rock EM. The definition and the hyperparameter settings of the five ML models can be found in [61,62,63,64,65]. The multivariate- quadratic equation (MQE) of MR was used to generate an empirical formula as expressed using Equation (23).

$$\begin{array}{cc}\mathrm{EM} =\hfill & -136.609 + 7.864\mathrm{PLI} - 1.056{\mathrm{PLI}}^2 - 17.038{\mathrm{P}}_{\mathrm{n}} + 1.378{\mathrm{P}}_{\mathrm{n}}{}^2\hfill \\ & + 4.265{\mathrm{V}}_{\mathrm{p}} + 1.084{\mathrm{V}}_{\mathrm{p}}{}^2 + 5.297\mathrm{SHRN} - 0.042{\mathrm{SHRN}}^2\hfill \end{array}$$

(23)

The optimal hyperparameter combination of each model and the corresponding performance results are shown in

Table 5. Performance results of five common ML models using test set.

Models	Performance Indices				Hyperparameters
	R2	RMSE	MAE	WI
ANN	0.8683	15.1724	10.8323	0.9619	Nh = 2; Ne = 4,4
SVR	0.8592	15.6918	11.7625	0.9591	C = 128; Rk = 0.25
ELM	0.8795	14.5124	10.2086	0.9665	Nes = 65
KELM	0.8987	13.3074	8.4755	0.9716	Rc = 128; Rk = 1.0
GRNN	0.9010	13.1593	8.2674	0.9717	Sf = 0.3
MQE	0.8318	17.1476	13.5849	0.9497	Equation (7)

Note: N_h—the hidden layers number; N_e—the number of neurons in the corresponding hidden layer; C—penalty parameter; R_k—RBF kernel parameter; N_es—the number of neurons in a single hidden layer; R_c—regularization coefficient; S_f—smoothing factor.

. The GRNN obtained better performance indices (R²: 0.9010; RMSE: 13.1593; MAE: 8.2674 and WI: 0.9717) than the other models using the test set. After this model, the prediction accuracy of the KELM and the ELM models is superior to the ANN and the SVR models. The worst prediction model is the MQE; its index values of R², RMSE, MAE and WI are, respectively, equal to 0.8318, 17.1476, 13.5849 and 0.9497. Figure 9 shows the rank scores of all comparison models. It can be obviously observed that the performance ranking of the six models is GRNN (24), KELM (20), ELM (16), ANN (12), SVR (8) and MQE (4).

The regression diagrams of the six comparison models in the testing phase are presented in Figure 10. As it can be seen in these diagrams, all models have a finite number of rock EM data points near the diagonal line. Compared with the other five models, the GRNN model obtained a good EM regression distribution in the range of 70 to 90 (Figure 10e). For the MQE model, most of the rock EM data points are away from the diagonal, which also indicates a large error between the EM value predicted by this model and the measured EM value. As a result, the GRNN model is the best prediction model among the six other models by means of the best performance indices. It obtains the highest rank score value and the best regression distribution.

Thus, the GRNN model is used to compare the performance with the PRO-RF model for predicting the rock EM. The bar chart tool is utilized to make an intuitive comparison between the measured EM and the EM predicted by the GRNN model and PRO-RF model, as shown in Figure 11. In Figure 11a, the predicted EM values of the No. 1, No. 7, No. 11, No. 13 and No. 15 samples by the PRO-RF model showed deviations. On the other hand, the GRNN model also does not accurately predict the EM values of the No. 2, No. 3, No. 4 and No. 19 samples, in addition to the above samples. The error analysis results of the PRO-RF model and the GRNN model for predicting the rock EM are shown in

Table 6. Error analysis comparison between the PRO-RF model and the GRNN model.

Models	Statistical Indices
	Min	Max	Median	Mean	St. E	St. D	Sum
PRO-RF	0.134	31.524	2.678	6.042	1.702	8.337	145.014
GRNN	0.328	35.992	3.477	8.267	2.135	10.458	198.417

Note: St. E—standard error value. St. D—standard deviation value.

. It can be intuitively observed that the error in the statistical indices based on the PRO-RF model is better than the GRNN model ones, such as the error sum of the PRO-RF model (145.014) being lower than the GRNN model (198.417).

The relative deviation (RD) distribution is also an effective tool to evaluate the model performances [64]. The definition of RD is the ratio of the error between the predicted value and the measured value to the measured value. This also means that the models with a better performance have lower RD values. As illustrated in Figure 12a, the maximum RD of the PRO-RF model is 44.77%, and most of the RD values are lower than 40%. Especially, there are 11 RD values lower than 10%. It should be noted that the maximum RD of the GRNN model is higher than 70%, and only 9 RD values are lower than 10%. Therefore, the PRO-RF model has a better predictive performance than the GRNN model for calculating the rock EM.

5.3. Sensitive Analysis

After determining the best prediction model, the sensitive analysis was carried out to calculate the parameter importance of the rock EM prediction. The average impact on model output magnitude (mean) is one of indices in the Shapley additive explanations (SHAP), which is used to represent the parameter importance, as shown in Figure 13. As can be seen in this graph, the porosity (Pn) has the larger importance value (30.52) than the other input variables in the EM prediction based on the PRO-RF model. After the Pn, the importance ranking of other three variables are the SHRN (3.8), Vp (4.09) and PLI (2.08). In addition, another SHAP index is the impact on the model output (SHAP value), which is used to describe the correlation between the input and output variables. As illustrated in Figure 14, only the Pn has the greatest negative correlation with EM. The SHRN, Vp and PLI are positively correlated with the EM, and the correlation decreased successively.

6. Conclusions

In this study, four hybrid RF models, named BSA-RF, MVO-RF, GEO-RF and PRO-RF were developed to predict the elasticity modulus (EM) of 120 rock samples. Four rock properties, named the porosity (Pn), P-wave velocity (Vp), Schmidt hammer rebound number (SHRN) and point load index (PLI), were considered as the main factors for the EM prediction. In addition, five ML models (i.e., ANN, SVR, ELM, KELM and GRNN) and one empirical formula were also developed to predict the rock EM and compare the predictive performance with the proposed hybrid RF models. The main conclusion of this paper can be listed as follows:

i.

Four hybrid RF models have obtained a good prediction accuracy by means of four performance indices. In particular, the PRO-RF model is the best model among them.

ii.

The GRNN model has a better predictive performance than the other ML models and the empirical formula. It results in the higher values of R² (0.9010) and WI (0.9717) and the lower values of RMSE (13.1593) and MAE (8.2674). However, these four optimized RF models are superior to the GRNN model.

iii.

The porosity (Pn) is the most important variable by means of the highest average impact value of 30.52 for predicting the rock EM. Meanwhile, the Pn is also the only variable negatively correlated with EM.

This paper proposes four effective hybrid RF models to predict the rock EM. It shows a successful application of soft techniques for a rock parameter prediction. Nevertheless, more various rocks should be collected into the integrated database to increase the prediction model’s accuracy. Furthermore, other rock properties such as the density, water content and UCS could be considered in the rock EM prediction.