*Article* **Performance of Statistical and Intelligent Methods in Estimating Rock Compressive Strength**

**Xuesong Zhang 1,\*, Farag M. A. Altalbawy 2,3, Tahani A. S. Gasmalla <sup>4</sup> , Ali Hussein Demin Al-Khafaji <sup>5</sup> , Amin Iraji <sup>6</sup> , Rahmad B. Y. Syah <sup>7</sup> and Moncef L. Nehdi 8,\***


**Abstract:** This research was conducted to forecast the uniaxial compressive strength (UCS) of rocks via the random forest, artificial neural network, Gaussian process regression, support vector machine, K-nearest neighbor, adaptive neuro-fuzzy inference system, simple regression, and multiple linear regression approaches. For this purpose, geo-mechanical and petrographic characteristics of sedimentary rocks in southern Iran were measured. The effect of petrography on geo-mechanical characteristics was assessed. The carbonate and sandstone samples were classified as mudstone to grainstone and calc-litharenite, respectively. Due to the shallow depth of the studied mines and the low amount of quartz minerals in the samples, the rock bursting phenomenon does not occur in these mines. To develop UCS predictor models, porosity, point load index, water absorption, P-wave velocity, and density were considered as inputs. Using variance accounted for, mean absolute percentage error, root-mean-square-error, determination coefficient (R<sup>2</sup> ), and performance index (PI), the efficiency of the methods was evaluated. Analysis of model criteria using multiple linear regression allowed for the development of a user-friendly equation, which proved to have adequate accuracy. All intelligent methods (with R<sup>2</sup> > 90%) had excellent accuracy for estimating UCS. The percentage difference of the average of all six intelligent methods with the measured value was equal to +0.28%. By comparing the methods, the accuracy of the support vector machine with radial basis function in predicting UCS was (R<sup>2</sup> = 0.99 and PI = 1.92) and outperformed all the other methods investigated.

**Keywords:** UCS; intelligent and statistical methods; prediction; sedimentary rocks

### **1. Introduction**

Stability of slopes, prediction of drilling rate, classification of rock masses, and modeling of foundations require knowledge of the uniaxial compressive strength (UCS) of the rocks for designing projects [1–3]. Indirect determination of the UCS in places where the preparation of standard samples is difficult requires lots of time and is expensive. Hence, various researchers have predicted the UCS of the limestones and sandstones using statistical and intelligent methods [4–7]. Aladejare et al. [8] collected empirical relationships and models between UCS and other rock characteristics from previous studies. Several models were developed to estimate the rock UCS using Gaussian process regression (GPR) [9–13], feedforward multilayer perceptron artificial neural network (FMP-ANN) [14–19], random forest algorithm (RFA) [20–23], adaptive neuro-fuzzy inference system (ANFIS) [24–28],

**Citation:** Zhang, X.; Altalbawy, F.M.A.; Gasmalla, T.A.S.; Al-Khafaji, A.H.D.; Iraji, A.; Syah, R.B.Y.; Nehdi, M.L. Performance of Statistical and Intelligent Methods in Estimating Rock Compressive Strength. *Sustainability* **2023**, *15*, 5642. https:// doi.org/10.3390/su15075642

Academic Editors: Mahdi Hasanipanah, Danial Jahed Armaghani, Jian Zhou and Jianjun Ma

Received: 27 December 2022 Revised: 26 February 2023 Accepted: 15 March 2023 Published: 23 March 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and multiple linear regression (MPLR) [3,4,7,8,25,29–31]. The results of Lawal et al.'s [9] study showed that the GPR method, with a correlation coefficient of almost 100%, is able to estimate the static and dynamic properties of sedimentary rocks. Moreover, a comparison of the RFA, MPLR, FMP-ANN methods in compressive strength estimation showed that FMP-ANN with the Levenberg–Marquardt algorithm has a higher accuracy than other methods [20]. The results of Matin et al.'s [22] study to select the effective variables using the random forest method showed that compression wave velocity is the most effective variable as an input for estimating compressive strength and the modulus of elasticity using predictive models. Hudaverdi [25] stated that the ANFIS method has a high efficiency in flyrock estimation with an average error of less than 8%. The results of the MPLR, ANFIS, and FMP-ANN methods in the UCS estimation showed that all three methods have a determination coefficient higher than 90%, while the ANFIS method has a better performance [26]. The comparison of the ANFIS, FMP-ANN and multiple regression methods by Yesiloglu-Gultekin and Gokceoglu [26] showed that the ANFIS method has higher accuracy for estimating compressive strength and the modulus of elasticity.

Mahmoodzadeh et al. [32] compared the K-nearest neighbor algorithm (KNNA), Gaussian process regression based on squared exponential kernel (GPR-SEK), support vector machine with radial basis function (SVR-RBF), and decision tree (DT) methods to forecast rock quality designation in a tunnel project and stated that the GPR-SEK method is more accurate than other methods. Xu et al. [33] forecasted the UCS of rock using intelligent technics. The SVM-RBF was used to predict UCS [34]. Rastegarnia et al. [19] used FMP-ANN and MPLR to predict the UCS of carbonates. They stated that FMP-ANN estimates the UCS more than the measured value. Trott et al. [35] used RFA to classify rock types. Barzegar et al. [36] predicted the UCS of travertine rocks using the SVM-RBF, FMP-ANN, and ANFIS methods and stated that the SVM-RBF showed higher accuracy than the other methods. Mohamad et al. [37] estimated the UCS of soft rocks using FMP-ANN and particle swarm optimization (PSO)-based ANN. Madhubabu et al. [6] used MPLR and FMP-ANN to estimate the UCS of the carbonate samples. Umrao et al. [24] used the ANFIS approach for estimating UCS based on density, porosity, and PWV. Moreover, using inteligent methods, Gül et al. [17] predicted the UCS of different rock types. Singh et al. [38] estimated the UCS of basalt samples via MPLR and ANFIS. Kaloop et al. [39] used GPR and multivariate adaptive regression splines (MARSs) to estimate rock UCS. They stated that the MARS showed higher acuracy than the GPR method. Some engineers and researchers are interested in simple empirical relationships using simple models such as simple and MPLR to estimate UCS. Therefore, simple empirical relationships are widely used to estimate rock UCS using statistical methods. Table 1 shows some of the relationships for estimating UCS by previous researchers.

This study was conducted to predict the UCS of sedimentary rocks based on porosity, point load index (PLI), density (D), water absorption by weight (WW), and P-wave velocity (PWV) using the FMP-ANN, GPR-SEK, KNNA, RFA, ANFIS, SVM, SR, and MPLR methods. Moreover, the types of kernel functions were investigated using the SVM method and the most accurate type of kernel function was introduced to estimate the UCS.

Sampling points, from 12 mines in the Bushehr province, south of Iran, are located between 50 and 52 degrees longitudes and 28 and 30 degrees latitudes. The mines are mainly travertine, limestone, and sandstone of the Aghajari and the Mishan formations.


**Table 1.** Relationships for estimating UCS by previous researchers.

#### **2. Methodology**

#### *2.1. Laboratory Tests*

Specimens with a diameter of 54 mm and a height to diameter ratio of 2 were prepared [46]. A wear device was used to parallel surfaces of specimens. Table 2 shows the methods used to measure geo-mechanical properties. Figure 1 displays some of the samples in laboratory tests.

**Table 2.** Methods used for performing tests.


**Figure 1.** Example samples for PWV, PLI, and UCS tests: (**a**) device for measuring PWV, (**b**) example samples after PLI test, and (**c**) sample under UCS test. **Figure 1.** Example samples for PWV, PLI, and UCS tests: (**a**) device for measuring PWV, (**b**) example samples after PLI test, and (**c**) sample under UCS test.

#### *2.2. Random Forest Algorithm (RFA) 2.2. Random Forest Algorithm (RFA)*

of the random forest [56].

The random forest method is one of the ensemble methods. In these methods, the model chosen for classification or regression is a combination of several models. Figure 2 The random forest method is one of the ensemble methods. In these methods, the model chosen for classification or regression is a combination of several models. Figure 2

shows the idea of the random forest algorithm. In this approach, each model issues its vote and the final result about the value is issued based on these votes [22,23,52]. The

combined methods is higher than the accuracy of its components. [53]. The RFA method has also been used in rock mechanics in recent years [20]. In the RFA method, the models used in the combined method, which are all of the decision tree type, form a forest. Each of the decision trees is made using a random selection of special attributes in each node to determine the branching. In other words, each tree is built based on the values of a random vector. These values have the identical scattering for all trees in the forest and are sampled independently. For classification, each tree issues its vote, and the final result is determined by the majority vote [54]. The number of trees and the number of chosen variables in each node are important parameters in the RFA [55]. In this method, by replacing the information every sampling time, some information is never sampled, and other data may be sampled several times. In other words, some input data for some trees will be out of the bag, that is, they will not participate in the creation of some trees. These data have the function of an internal validator for each tree, which is performed by estimating the outof-bag error. If the out-of-bag data itself is predicted through trees, there will be an error for these predictions, and the average of these errors is called the out-of-bag error. This error indicates the influence of the unselected samples on the error rate of the final result

shows the idea of the random forest algorithm. In this approach, each model issues its vote and the final result about the value is issued based on these votes [22,23,52]. The general principles of group training techniques are based on the assumption that their accuracy is higher than other training algorithms [53]. On the other hand, the accuracy of combined methods is higher than the accuracy of its components [53]. The RFA method has also been used in rock mechanics in recent years [20]. In the RFA method, the models used in the combined method, which are all of the decision tree type, form a forest. Each of the decision trees is made using a random selection of special attributes in each node to determine the branching. In other words, each tree is built based on the values of a random vector. These values have the identical scattering for all trees in the forest and are sampled independently. For classification, each tree issues its vote, and the final result is determined by the majority vote [54]. The number of trees and the number of chosen variables in each node are important parameters in the RFA [55]. In this method, by replacing the information every sampling time, some information is never sampled, and other data may be sampled several times. In other words, some input data for some trees will be out of the bag, that is, they will not participate in the creation of some trees. These data have the function of an internal validator for each tree, which is performed by estimating the out-of-bag error. If the out-of-bag data itself is predicted through trees, there will be an error for these predictions, and the average of these errors is called the out-of-bag error. This error indicates the influence of the unselected samples on the error rate of the final result of the random forest [56]. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 5 of 22

**Figure 2.** Concept of random forest algorithm. **Figure 2.** Concept of random forest algorithm.

#### *2.3. Gaussian Process Regression Based on Squared Exponential Kernel (GPR-SEK)*

*2.3. Gaussian Process Regression Based on Squared Exponential Kernel (GPR- SEK)*  Consider a d data set with n measurements: = ൛(, y )| = 1, . . . . . }, where xi is the input vector with D dimension and *yi* is the target output. This set, consisting of two components, input and output, will be denoted as measured points. To simplify the problem, the inputs of the collection are aggregated at = {, *x*ଶ,..., } matrix and the outputs are also combined at = ൛, yଶ,..., } matrix. Regression based on the data set *d* creates a new input *x\** to arrive at the predicted distribution for the corresponding values Consider a d data set with n measurements: *d* = {(*x<sup>i</sup>* , *y<sup>i</sup>* )|*i* = 1, . . . , *n*} , where *x<sup>i</sup>* is the input vector with D dimension and *y<sup>i</sup>* is the target output. This set, consisting of two components, input and output, will be denoted as measured points. To simplify the problem, the inputs of the collection are aggregated at *X* = {*x<sup>i</sup>* , *x*<sup>2</sup> , . . . , *xn*} matrix and the outputs are also combined at *Y* = {*y<sup>i</sup>* , *y*<sup>2</sup> , . . . , *yn*} matrix. Regression based on the data set *d* creates a new input *x\** to arrive at the predicted distribution for the corresponding values of the measured *y*\* data. The Gaussian process (GP) is a group of accidental parameters, a restricted number of which are combined with Gaussian distributions (GDs) [57]. The GP is a generalization of GD. The GD is actually scattered between accidental parameters,

of the measured *y*\* data. The Gaussian process (GP) is a group of accidental parameters, a restricted number of which are combined with Gaussian distributions (GDs) [57]. The GP is a generalization of GD. The GD is actually scattered between accidental parameters, while GP represents scattering between functions. The *f(x)* GP is described using the *m(x)*

) = (() − ())((<sup>ᇱ</sup>

points. The GP can be described as Equation (3).

() ∼ ((), c(x,x<sup>ᇱ</sup>

where () represents the arbitrary regression function and is the noise of the Gaussian function with zero mean and ଶvariance (i.e., ∼ (0, ଶ)). Furthermore, it is supposed that = ሾ(ଵ), f(ଶ), ..., f()ሿ் has a performance according to the GP (i.e., (|) = (0, C)). Here, C is the covariance matrix with the , = (, x) domains.

(, ) = ൮(ଵ, ଵ) c(ଵ, ଶ) ... <sup>c</sup>(ଵ, )

⋮ ⋮ ⋮ ⋮

(ଶ, ଵ) c(ଶ, ଶ) ... c(ଶ, )

Usually, for simplification, the value of the average function is considered equal to zero [58]. In the *GP*, the correlation between the target and the input vector is based on

() = (()) (1)

) is the covariance or kernel function, which is com-

= ()+ (4)

(, ଵ) c(, ଶ) ... c(, ) <sup>൲</sup> (5)

))) (2)

)) (3)

) − (<sup>ᇱ</sup>

(, <sup>ᇱ</sup>

In relationships 1 and 2, (, <sup>ᇱ</sup>

puted at the *x* and <sup>ᇱ</sup>

Equation (4).

while GP represents scattering between functions. The *f(x)* GP is described using the *m(x)* average and covariance functions according to Equations (1) and (2).

$$m(\mathbf{x}) = E(f(\mathbf{x})) \tag{1}$$

$$\mathcal{L}(\mathbf{x}, \mathbf{x}') = \mathbb{E}(f(\mathbf{x}) - m(\mathbf{x})) \left( f(\mathbf{x}') - m(\mathbf{x}') \right) \tag{2}$$

In relationships 1 and 2, *c*(*x*, *x* 0 ) is the covariance or kernel function, which is computed at the *x* and *x* 0 points. The GP can be described as Equation (3).

$$f(\mathbf{x}) \sim \text{GP}\left(\mathfrak{m}(\mathbf{x}), \mathfrak{c}(\mathbf{x}\mathbf{x}')\right) \tag{3}$$

Usually, for simplification, the value of the average function is considered equal to zero [58]. In the *GP*, the correlation between the target and the input vector is based on Equation (4).

$$y\_i = f(\mathbf{x}\_i) + \varepsilon \tag{4}$$

where *f*(*xi*) represents the arbitrary regression function and *ε* is the noise of the Gaussian function with zero mean and *σ* <sup>2</sup> variance (i.e., *<sup>ε</sup>* <sup>∼</sup> *<sup>N</sup>* 0, *σ* 2 ). Furthermore, it is supposed that *f* = [ *f*(*x*1), f(*x*2), . . . , f(*xn*)]*<sup>T</sup>* has a performance according to the GP (i.e., *p*(*f* |*X*) = *N*(0, C)). Here, C is the covariance matrix with the *ci*,*<sup>j</sup>* = *c xi* , x*<sup>j</sup>* domains.

$$\mathbf{C}(X,X) = \begin{pmatrix} c(\mathbf{x}\_1,\mathbf{x}\_1) \ c(\mathbf{x}\_1,\mathbf{x}\_2) \ \dots \ c(\mathbf{x}\_1,\mathbf{x}\_n) \\ c(\mathbf{x}\_2,\mathbf{x}\_1) \ c(\mathbf{x}\_2,\mathbf{x}\_2) \ \dots \ c(\mathbf{x}\_2,\mathbf{x}\_n) \\ \vdots \ \vdots \ \vdots \ \vdots \\ c(\mathbf{x}\_{\mathsf{n}},\mathbf{x}\_1) \ c(\mathbf{x}\_{\mathsf{n}},\mathbf{x}\_2) \ \dots \ c(\mathbf{x}\_{\mathsf{n}},\mathbf{x}\_n) \end{pmatrix} \tag{5}$$

The *ci*,*<sup>j</sup>* is the covariance between the latent function values of *f*(*xi*) and *f xj* . GP regression is used to calculate the predicted scattering for the *f\** function values in the test points of *X* ∗ = - *x* ∗ 1 , *x* ∗ 2 . . . *x* ∗ *m* . The distribution of y depends on the values of *f*, which is represented by an isotropic Gaussian as follows.

$$p\left(y \middle| f, \varkappa\right) = N(F, \sigma\_n^2 I) \tag{6}$$

In relation (6), *I* is the identity matrix. According to the characteristics of the Gaussian function, the marginal distribution of *y* is determined as follows.

$$p(y|X) = \int \left. p\left(y|f, X\right) p\left(f\middle|X\right) df = N(0, \mathbb{C} + \sigma\_n^2 I\right) \tag{7}$$

The integrated distribution of the observation data values, that is, the desired output, and the function values at the test points are written as follows [32].

$$
\begin{bmatrix}
\boldsymbol{y} \\
\boldsymbol{f}^\*
\end{bmatrix} \sim N\left(\mathbf{0}, \begin{bmatrix}
\mathbf{C}(\mathbf{X}, \mathbf{X}) + \sigma^2 I & \mathbf{C}(\mathbf{X}, \mathbf{X}\_\*) \\
\mathbf{C}(\mathbf{X}\_\*, \mathbf{X}) & \mathbf{C}(\mathbf{X}\_\*, \mathbf{X}\_\*)
\end{bmatrix}\right) \tag{8}
$$

According to relation (3), and using standard rules to limit Gaussian, the following conditional distribution can be obtained.

$$p(f\_\*|X, y, X\_\*) \sim N\left(\overline{f\_\*}, c(f\_\*)\right) \tag{9}$$

$$\bar{f}\_\* = \mathbb{C}(X\_\*, X) \left[ \mathbb{C}(\mathbb{C}, \mathbb{C}) + \sigma^2 I \right]^{-1} y \tag{10}$$

#### *2.4. The SVM-RBF*

To achieve the least error related to the test set, the SVM-RBF approach fits a linear line with epsilon (*ε*) thickness on the data [59]. In this method, a function such as *f(x)* = *m.x* + *b* is used for forecasting, where m is weight vector and *x* and *b* are weights.

For minimizing weight vector and test error, this technique utilizes error functions for ignoring errors that are at a determined range from the real errors [60]. Hence, some deviation (derived from Equation (11)) from *ε* must be overlooked by including Equation (11) in Equation (12), which considers the *ξ* − *i* and *ξ* + *i* deficiency parameters. Finally, the error values are optimized via Equation (12) using structural error minimization

$$|\mathfrak{f}|\_{\varepsilon} = \begin{bmatrix} 0 & \text{if } |\mathfrak{f}| \le \varepsilon \\ |\mathfrak{f}| - \varepsilon & \text{otherwise} \end{bmatrix} \tag{11}$$

$$\begin{array}{c} \text{Minimize}: \{ (\left\| m \right\| ^2) \* 1/2 \} + \{ (\sum\_{i=1}^N (\xi\_i^+ + \xi\_i^-)) \* \mathbb{C} \} \\ \text{ $\varepsilon$ Constrains:}: \begin{bmatrix} m.x\_i + b + \xi\_i^+ - y\_i \le \varepsilon & i = 1, 2, \dots, N \\ y\_i - (b + m.x\_i) \le \xi\_i^- + \varepsilon & i = 1, 2, \dots, N \\ \xi\_i^+ \ge 0 \quad \xi\_i^- \ge 0 & i = 1, 2, \dots, N \end{bmatrix} \end{array} \tag{12}$$

In Equation (12), {(k*m*k 2 )∗1/2} is the supervisory part, *N* is number of samples, *ε* is the allowable error, *C* is the complexity balance coefficient, and the *ε* values are the acceptable error range. As with the GPR method, various kernel functions are used in the SVM method [61]. Radial basis function (RBF), which is the most important kernel function, was used in the current research [62].

#### *2.5. K Nearest Neighbor Algorithm (KNNA)*

The KNNA is based on sample and performs classification based on K nearest neighbors. This method performs classification based on the similarity of the data. In fact, for each new test data, it calculates the K nearest neighbor distances and determines a label similar to the dominant label of this k neighbor for the desired point [63]. This method was introduced as a nonparametric method and does not make any assumption on the distribution of inputs. Therefore, it is extensively used in various fields [64].

In the KNNA classifier, an unknown value, is recognized by the similarity between known trained or labeled values based on the calculation of the distance between unknown values and labeled values. Then, K of the nearest values are selected as the basis for classification, and the unknown value (x test) is assigned to the class that has the most values among the closest values. For this purpose, three factors affect the KNNA classification: (1) the number of K of the neighbor and the changing of the value of K, which may amendment the classification results; (2) labeled dataset; therefore, adding or eliminating any value to the training samples affects the final results of the KNNA classifier; (3) the distance criterion. In KNNA, Euclidean distance is usually used as a distance criterion to measure the distance between two values [64,65]. This algorithm, as with the other algorithms used in this research, after examining the data in the program environment, divides the data into two parts, training data and test data, and builds the K nearest neighbor model and enters the training data into the model to train the model. Next, to determine the precision of the method, the test data is entered into the model for prediction and to evaluate the prediction accuracy in comparison with the labels of the test data [65,66].

#### *2.6. ANFIS and FMP-ANN*

The ANFIS and FMP-ANN methods have been widely introduced and described by previous researchers [18,25,67–70]. The transfer functions of neurons, membership functions, type of fuzzy inference system, and data training methods in these two methods are mentioned in the results section.

In the SVM-RBF, ANFIS, KNNA, *GPR-SEK*, and RFA methods, 30% and 70% of the whole data were used for the testing and training the models, respectively.

### *2.7. Performance Evaluation of Results*

To appraise the methods, the correlation coefficient, the MAPE %, the RMSE, VAF, and the PI are defined in the form of Equations (13)–(16).

$$\text{MAPE} = \frac{1}{n} \sum\_{i=1}^{n} \left| \frac{y - y'}{y} \right| \ast 100 \tag{13}$$

$$RMSE = \frac{1}{s^2 n} \sum\_{i=1}^{n} (y - y')^2 \tag{14}$$

$$VAF = 100\left[1 - \frac{var(y - y')}{var(y)}\right] \tag{15}$$

$$\text{PI} = \text{R}^2 + (\text{VAF}/100) - \text{RMSE} \tag{16}$$

In relationships 13 to 16, *y* is the value of the variable measured, *y* 0 is the predicted UCS, and *n* is the total data and *s* 2 is the sample variance. Equation (17) was used to normalize the data between −1 and 1.

$$X\_i = 2\left(\frac{X - X\_{\min}}{X\_{\max} - X\_{\min}}\right) - 1\tag{17}$$

In Equation (17), *x* is the measured variable, *Xmin* is the minimum of the data, and *Xmax* is the maximum of the data.

#### **3. Results and Discussion**

#### *3.1. Geomechanical Properties of Samples*

The maximum, minimum, and average engineering properties of 65 samples (37 samples of limestone, 11 samples of argillaceous limestone, and 17 samples of sandstone) are presented in Table 3. The average Es and UCS are 14.95 GPa and 37.54 MPa, respectively. Physical, mineralogical, and petrographic characteristics cause changes in the compressive strength of rocks [5,18,71]. Moreover, lithological properties such as the type of rock, the amount and type of minerals, the type of cement between the particles that comprise the rock and physical properties such as the amount of moisture, porosity, and water absorption have a significant effect on the compressive strength and, consequently, on the stability of mines [72,73]. As the amount of clay minerals increases, the resistance of the samples decreases [4,74]. Water absorption causes the swelling and instability of the mine wall in samples with a high percentage of clay minerals [75–77]. The number of joints changes the geomechanical properties and instability of the rocky slopes [78–80]. The engineering properties and stability of structures are affected by the amount of moisture [81].


**Table 3.** Laboratory results on sample.

In the sandstone samples of the present research, clay and gypsum were found. The cement of the samples is gypsum and calcite. The grains of these sandstones are semirounded to angular and have a moderate to poor grading. The study of the thin sections of the samples showed that the samples with higher clay content show lower resistance characteristics. Swelling clay minerals (such as montmorillonite) were not observed in the investigated samples. Generally, rock bursting occurs in deep mines and in quartz-rich rocks [12]. Because the depth of the studied mines is less than 50 m and the overburden stress is negligible, the risks of rock bursting have not been reported in them so far.

### *3.2. Petrographic Features*

Texture has a special effect on the engineering properties of sedimentary rocks [82]. In carbonates, the rock texture is very different, but their mineralogy is not much different [83,84]. According to microscopic studies, the most basic mineral of limestone rocks was calcite. Mishan formation limestone rocks, with an early Miocene age, based on the Dunham [51] classification, were classified in the range of mudstone to grainstone. Moreover, argillaceous limestone samples of this formation were classified in the mudstone to packstone categories. Sandstone samples of the Aghajari formation with an upper miocene age were classified as calc-litharenite according to the Folk [50] classification. These rocks consist of carbonaceous rock fragments (26 to 75), volcanic gravel (2 to 35%), meta-morphic fragments (2 to 18%), feldspar (1 to 10%), dark minerals (1 to 8%), quartz (0 to 22%), and chert (2 to 11%). *Sustainability* **2023**, *15*, x FOR PEER REVIEW 9 of 22 *3.3. Influence of Independent Variables on the UCS* 

#### *3.3. Influence of Independent Variables on the UCS* Figure 3 shows the effect of variables on the UCS. There is a reasonable tendency

Figure 3 shows the effect of variables on the UCS. There is a reasonable tendency among these characteristics. The UCS decreases with increasing WW and porosity. This Figure shows that porosity has the highest effect on the UCS. Numerous studies have reported linear relationships with high accuracy between the point load index (PLI) and UCS [8,30]. The results of the relationship between PVW and UCS show that PWV displays a high accuracy to estimate UCS (Figure 3). among these characteristics. The UCS decreases with increasing WW and porosity. This Figure shows that porosity has the highest effect on the UCS. Numerous studies have reported linear relationships with high accuracy between the point load index (PLI) and UCS [8,30]. The results of the relationship between PVW and UCS show that PWV displays a high accuracy to estimate UCS (Figure 3).

(**e**)

accuracy of the forecasted and actual UCS relationships were assessed.

**Figure 3.** Effect of (**a**) point load index (PLI), (**b**) water absorption by weight (WW), (**c**) %porosity,

For predicting the UCS of the rocks, some empirical relationships have been proposed (Table 1). In the current research, for each of the 65 samples of the present study, based on each of the empirical relationships in Table 1, UCS was predicted. Finally, the

Figure 4 displays the data scattering and the precision of correlation using PI and R2. The results revealed that there is good compatibility between actual UCS and the estimated one using previous studies (Figure 4). A performance index (PI) was introduced by Yagiz et al. [85] for evaluating empirical equations and models. The value of this index is equal to two in the best case, and the lower it is, the lower the relationship performance.

**Figure 3.** *Cont*.

(**d**) density (D), and (**e**) P−wave velocity (PWV) on the UCS.

*3.4. Evaluation of Previous Emperical Relationships* 

 (**c**) (**d**)

 (**a**) (**b**)

**Figure 3.** Effect of (**a**) point load index (PLI), (**b**) water absorption by weight (WW), (**c**) %porosity, (**d**) density (D), and (**e**) P−wave velocity (PWV) on the UCS. **Figure 3.** Effect of (**a**) point load index (PLI), (**b**) water absorption by weight (WW), (**c**) %porosity, (**d**) density (D), and (**e**) P−wave velocity (PWV) on the UCS.

#### *3.4. Evaluation of Previous Emperical Relationships 3.4. Evaluation of Previous Emperical Relationships Sustainability* **2023**, *15*, x FOR PEER REVIEW 10 of 22

*3.3. Influence of Independent Variables on the UCS* 

plays a high accuracy to estimate UCS (Figure 3).

Figure 3 shows the effect of variables on the UCS. There is a reasonable tendency among these characteristics. The UCS decreases with increasing WW and porosity. This Figure shows that porosity has the highest effect on the UCS. Numerous studies have reported linear relationships with high accuracy between the point load index (PLI) and UCS [8,30]. The results of the relationship between PVW and UCS show that PWV dis-

For predicting the UCS of the rocks, some empirical relationships have been proposed (Table 1). In the current research, for each of the 65 samples of the present study, based on each of the empirical relationships in Table 1, UCS was predicted. Finally, the For predicting the UCS of the rocks, some empirical relationships have been proposed (Table 1). In the current research, for each of the 65 samples of the present study, based on each of the empirical relationships in Table 1, UCS was predicted. Finally, the accuracy of the forecasted and actual UCS relationships were assessed. As can be seen, although the correlation coefficient is high, the performance index is negative, which indicates the poor performance of the previous researchers' relationships in estimating the UCS of the studied rocks (Figure 4). For this reason, various researchers

accuracy of the forecasted and actual UCS relationships were assessed. Figure 4 displays the data scattering and the precision of correlation using PI and R2. The results revealed that there is good compatibility between actual UCS and the estimated one using previous studies (Figure 4). A performance index (PI) was introduced by Yagiz et al. [85] for evaluating empirical equations and models. The value of this index is equal to two in the best case, and the lower it is, the lower the relationship performance. Figure 4 displays the data scattering and the precision of correlation using PI and R<sup>2</sup> . The results revealed that there is good compatibility between actual UCS and the estimated one using previous studies (Figure 4). A performance index (PI) was introduced by Yagiz et al. [85] for evaluating empirical equations and models. The value of this index is equal to two in the best case, and the lower it is, the lower the relationship performance. As can be seen, although the correlation coefficient is high, the performance index is negative, which indicates the poor performance of the previous researchers' relationships in estimating the UCS of the studied rocks (Figure 4). For this reason, various researchers have emphasized that empirical relationships should be determined for each region [85]. have emphasized that empirical relationships should be determined for each region [85]. The type of rock, strength amount, method of conducting experiments, the test conditions (such as loading rate), and the petrography of the specimens of a specific study reveal the applicability of the proposed relationships for forecasting the UCS of rocks in other regions. For example, the average UCS of the samples of Edet [3] study was 32.22 MPa, and the average UCS of the present study is 37.54 MPa, which shows that the resistance of the Edet [3] study samples is lower than the resistance of the current research samples. As a result, the UCS values estimated from this researcher's relationship are mostly below the diagonal line (Figure 4d). The sample breaks faster and shows more resistance when the loading rate is increased.

**Figure 4.** Measured UCS versus predicted UCS based on relationships of (**a**) Kahraman [43] and Teymen and Mengüç [40], (**b**) Aliyu et al. [30], Aldeeky and Al Hattamleh [42], (**c**) Salehin [41], Abdi and Khanlari [4], (**d**) Azimian [29], Edet [3], and Selcuk and Nar [31]. **Figure 4.** Measured UCS versus predicted UCS based on relationships of (**a**) Kahraman [43] and Teymen and Mengüç [40], (**b**) Aliyu et al. [30], Aldeeky and Al Hattamleh [42], (**c**) Salehin [41], Abdi and Khanlari [4], (**d**) Azimian [29], Edet [3], and Selcuk and Nar [31].

In the current work, MPLR analysis was performed using Minitab software (version 18). Equation (18) was developed to predict the UCS using this method. Various criteria

Various statistics (Tables 4 and 5) were used to evaluate relationship 18. The Durbin– Watson statistic DWS) and variance inflation factor (VIF) are used to evaluate the independence of errors and the correlation of independent variables, respectively [86]. The results showed that there is no problem in using relationships in terms of these two criteria because the DW is in the range of 1.5 to 2 and the VIF value is less than 10 (Table 5). Analysis of variance (ANOVA) results (Sig. <0.00) show that the model has been well developed

*3.5. Multiple Linear Regression (MPLR)* 

The type of rock, strength amount, method of conducting experiments, the test conditions (such as loading rate), and the petrography of the specimens of a specific study reveal the applicability of the proposed relationships for forecasting the UCS of rocks in other regions. For example, the average UCS of the samples of Edet [3] study was 32.22 MPa, and the average UCS of the present study is 37.54 MPa, which shows that the resistance of the Edet [3] study samples is lower than the resistance of the current research samples. As a result, the UCS values estimated from this researcher's relationship are mostly below the diagonal line (Figure 4d). The sample breaks faster and shows more resistance when the loading rate is increased.

### *3.5. Multiple Linear Regression (MPLR)*

In the current work, MPLR analysis was performed using Minitab software (version 18). Equation (18) was developed to predict the UCS using this method. Various criteria to evaluate this relationship are presented below.

Various statistics (Tables 4 and 5) were used to evaluate relationship 18. The Durbin –Watson statistic DWS) and variance inflation factor (VIF) are used to evaluate the independence of errors and the correlation of independent variables, respectively [86]. The results showed that there is no problem in using relationships in terms of these two criteria because the DW is in the range of 1.5 to 2 and the VIF value is less than 10 (Table 5). Analysis of variance (ANOVA) results (Sig. < 0.00) show that the model has been well developed using MPLR. Sig. values (related to T-test) in Table 5 indicate the presence of variables in the multivariate regression output relationship. The constant value, density, and water absorption were removed from Equation (18) because the sig. value is more than 0.05.

**Table 4.** Multiple linear regression results to estimate UCS and Es.



**Table 5.** Evaluation criteria of coefficients for relationship 18.

The normality of the error distribution is also one of the other criteria for evaluating empirical relationships. The normal distribution of errors related to the model provided by MPLR method shows that the proposed model can be used to estimate UCS (Figure 5).

UCS = 5.03PWV − 1.735

**Table 4.** Multiple linear regression results to estimate UCS and Es.

+ 2.667PLI 0.88 1.10 1.08 87.85 0.66 1.93 F-value = 79.37

**MAPE %** 

**(MPa)** 

**Table 5.** Evaluation criteria of coefficients for relationship 18.

**Equation R2 RMSE** 

**Figure 5.** Normal state of the residues of model 18 (**a**) histogram of residuals, and (**b**) normal probability plot. **Figure 5.** Normal state of the residues of model 18 (**a**) histogram of residuals, and (**b**) normal probability plot.

using MPLR. Sig. values (related to T-test) in Table 5 indicate the presence of variables in the multivariate regression output relationship. The constant value, density, and water absorption were removed from Equation (18) because the sig. value is more than 0.05.

**VAF** 

**Term Coefficients T-Value Significant Level (Sig.) VIF (Variance** 

The normality of the error distribution is also one of the other criteria for evaluating empirical relationships. The normal distribution of errors related to the model provided by MPLR method shows that the proposed model can be used to estimate UCS (Figure 5).

Constant −32.1 −1.34 0.187 -

PWV 5.03 2.44 0.018 7.58 D 21.4 1.82 0.074 3.02 WW 0.281 0.35 0.728 3.81 −1.735 −3.97 0.000 3.64 PLI 2.667 3.05 0.003 3.77

**% PI DWS ANOVA Results Eq. No.** 

*<sup>p</sup>*-value = 0.00 (18)

**Inflation Factor)** 

#### *3.6. The Results of Modeling Using RFA and GPR-SEK Methods*

*3.6. The Results of Modeling Using RFA and GPR-SEK Methods*  The RFA modeling was conducted using the R (version R4.2.1) software [54]. The GPR-SEK model was conducted using MATLAB software (MATLAB 2016b). In the RF method, the 10-fold cross-validation method was used to control the number of chosen parameters in each node of tree (m-try) and the number of trees (n-tree). According to this The RFA modeling was conducted using the R (version R4.2.1) software [54]. The GPR-SEK model was conducted using MATLAB software (MATLAB 2016b). In the RF method, the 10-fold cross-validation method was used to control the number of chosen parameters in each node of tree (m-try) and the number of trees (n-tree). According to this method, the number of 500 trees and five variables in each node has delivered the most satisfactory conditions for modeling. Therefore, these values were used for modeling purposes.

method, the number of 500 trees and five variables in each node has delivered the most satisfactory conditions for modeling. Therefore, these values were used for modeling purposes. The random forest method works well for large amounts of data and has high accuracy. In the random forest method, because the amount of error decreases with the increase of trees, 500 trees were used to develop the model. Upon model execution, the results were evaluated using an out-of-bag (OOB) error estimation. The model was appraised by the test data, the results of which are presented in Figure 6. One of the advantages of the random forest algorithm is that it can determine the importance of variables in a problem. In this research, the significance of the inputs was achieved using the The random forest method works well for large amounts of data and has high accuracy. In the random forest method, because the amount of error decreases with the increase of trees, 500 trees were used to develop the model. Upon model execution, the results were evaluated using an out-of-bag (OOB) error estimation. The model was appraised by the test data, the results of which are presented in Figure 6. One of the advantages of the random forest algorithm is that it can determine the importance of variables in a problem. In this research, the significance of the inputs was achieved using the Gini significance index [54]. The results showed that porosity has higher importance than other parameters. In Figure 6, the error histogram, the graph of the measured, the forecasted UCS using the RFA method, and the GPR-SEK model are drawn. The GPR-SEK model was implemented based on the squared exponential kernel function. As can be seen in the figure, the results are close to the bisector line, and it can be said that the values have been predicted with good accuracy. Theoretically, if R<sup>2</sup> equals 100%, all the observed values will be similar to the fitted values and all the data points will be on the fitted line [87]. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 12 of 22 Gini significance index [54]. The results showed that porosity has higher importance than other parameters. In Figure 6, the error histogram, the graph of the measured, the forecasted UCS using the RFA method, and the GPR-SEK model are drawn. The GPR-SEK model was implemented based on the squared exponential kernel function. As can be seen in the figure, the results are close to the bisector line, and it can be said that the values have been predicted with good accuracy. Theoretically, if R2 equals 100%, all the observed values will be similar to the fitted values and all the data points will be on the fitted line [87].

**Figure 6.** *Cont*.

GPR-SEK for all data.

forecasting UCS.

*3.7. The FMP-ANN Results*

(**c**) (**d**)

**Figure 6.** Results of (**a**) RFA for test data, (**b**) RFA for all data, (**c**) GPR-SEK for test data, and (**d**)

The FMP-ANN is widely used in engineering [88,89]. In the current study, for predicting UCS, various neurons in a hidden layer were investigated to develop optimal models. Based on equations proposed by previous researchers, the number of hidden layer neurons changes were determined (Table 6). The calculated number of hidden layer neurons changed from one to eleven according to Table 6. In this study, by checking this range using the FMP-ANN, this range was evaluated to achieve the ideal model architecture for [87].

(**a**) (**b**)

**Figure 6.** Results of (**a**) RFA for test data, (**b**) RFA for all data, (**c**) GPR-SEK for test data, and (**d**) GPR-SEK for all data. **Figure 6.** Results of (**a**) RFA for test data, (**b**) RFA for all data, (**c**) GPR-SEK for test data, and (**d**) GPR-SEK for all data.

Gini significance index [54]. The results showed that porosity has higher importance than other parameters. In Figure 6, the error histogram, the graph of the measured, the forecasted UCS using the RFA method, and the GPR-SEK model are drawn. The GPR-SEK model was implemented based on the squared exponential kernel function. As can be seen in the figure, the results are close to the bisector line, and it can be said that the values have been predicted with good accuracy. Theoretically, if R2 equals 100%, all the observed values will be similar to the fitted values and all the data points will be on the fitted line

#### *3.7. The FMP-ANN Results 3.7. The FMP-ANN Results*

The FMP-ANN is widely used in engineering [88,89]. In the current study, for predicting UCS, various neurons in a hidden layer were investigated to develop optimal models. Based on equations proposed by previous researchers, the number of hidden layer neurons changes were determined (Table 6). The calculated number of hidden layer neurons changed from one to eleven according to Table 6. In this study, by checking this range using the FMP-ANN, this range was evaluated to achieve the ideal model architecture for The FMP-ANN is widely used in engineering [88,89]. In the current study, for predicting UCS, various neurons in a hidden layer were investigated to develop optimal models. Based on equations proposed by previous researchers, the number of hidden layer neurons changes were determined (Table 6). The calculated number of hidden layer neurons changed from one to eleven according to Table 6. In this study, by checking this range using the FMP-ANN, this range was evaluated to achieve the ideal model architecture for forecasting UCS.

forecasting UCS. **Table 6.** Proposed equations by previous researchers to estimate the number of hidden layer neurons.


*N*<sup>0</sup> and *N<sup>i</sup>* are the numbers of input and output neurons, respectively.

The used FMP-ANN method has a hidden layer with five inputs (PWV, point load index (PLI), porosity, density, and water absorption) and one output (UCS). Using MATLAB software, the Levenberg Marquardt (LM) training algorithm was used to train the network. The neuron transfer functions were the selected Sigmoid between the input layers and hidden layers and the Purelin between the hidden layers and output layers. In FMP-ANN modeling, the percentages of the validation, test, and training data in the present study were randomly selected as 15%, 15%, and 70% of the total data, respectively. The validation set is used to prevent overfitting, the training group is used to determine weights, and the test group is used to evaluate the FMP-ANN results [97–100]. The results showed that the third neuron is the most accurate neuron for estimating UCS. Figure 7 displays the optimal FMP-ANN chart achieved in the current research.

rons.

**References** 

Kanellopoulos and Wilkinson

**Figure 7.** Used FMP-ANN structure. **Figure 7.** Used FMP-ANN structure.

Figure 8 shows the error variations in optimum results. The lowest error in epoch 4 was obtained for predicting the UCS (Figure 8). Moreover, in this research, the results of the FMP-ANN to estimate UCS have been compared with several methods. It was found Figure 8 shows the error variations in optimum results. The lowest error in epoch 4 was obtained for predicting the UCS (Figure 8). Moreover, in this research, the results of the FMP-ANN to estimate UCS have been compared with several methods. It was found that the accuracy of all methods was very high (the coefficient of determination is more than 97%). *Sustainability* **2023**, *15*, x FOR PEER REVIEW 14 of 22

**Table 6.** Proposed equations by previous researchers to estimate the number of hidden layer neu-

**Equations** 

( + )

<sup>ଶ</sup> + ) − 3

**Neuron Numbers Calculated for This Study** 

Hecht-Nielsen [90] ≤3 ≤2∗ + 1 Hush [91] 3 3

Ripley [92] 3 ( + )/2

Wang [94] 1 2/3

[96] 1 2

The used FMP-ANN method has a hidden layer with five inputs (PWV, point load

index (PLI), porosity, density, and water absorption) and one output (UCS). Using MATLAB software, the Levenberg Marquardt (LM) training algorithm was used to train the network. The neuron transfer functions were the selected Sigmoid between the input layers and hidden layers and the Purelin between the hidden layers and output layers. In FMP-ANN modeling, the percentages of the validation, test, and training data in the present study were randomly selected as 15%, 15%, and 70% of the total data, respectively. The validation set is used to prevent overfitting, the training group is used to determine weights, and the test group is used to evaluate the FMP-ANN results [97–100]. The results showed that the third neuron is the most accurate neuron for estimating UCS. Figure 7

Kaastra and Boyd [95] 2 ඥ ∗

N0 and Ni are the numbers of input and output neurons, respectively.

displays the optimal FMP-ANN chart achieved in the current research.

Paola [93] 11 2+ ∗ + 0.5 <sup>∗</sup> (

**Figure 8.** The FMP−ANN results for the optimum model (**a**) error reduction trend and (**b**) corelation cofficient between measured and predicted UCS. **Figure 8.** The FMP−ANN results for the optimum model (**a**) error reduction trend and (**b**) corelation cofficient between measured and predicted UCS.

#### *3.8. The KNNA Results 3.8. The KNNA Results*

10). Figure 10 shows the KNNA results for estimating the UCS.

To apply the KNNA method to the data and determine the best K value, the KNNA was written in the form of a program in MATLAB software, which was run 310 times for K values from 1 to 30 programs; moreover, the amount of error was measured (Figure 9). Of the total data, 70% and 30% were used to train and test the model. The results showed that the lowest estimation error of the UCS was obtained at K = 2 (Figure 9). The error of this network for estimating the UCS with respect to the K values is equal to 0.11 (Figure To apply the KNNA method to the data and determine the best K value, the KNNA was written in the form of a program in MATLAB software, which was run 310 times for K values from 1 to 30 programs; moreover, the amount of error was measured (Figure 9). Of the total data, 70% and 30% were used to train and test the model. The results showed that the lowest estimation error of the UCS was obtained at K = 2 (Figure 9). The error of this network for estimating the UCS with respect to the K values is equal to 0.11 (Figure 10). Figure 10 shows the KNNA results for estimating the UCS.

**Figure 9.** Obtained RMSE for forecasting UCS based on the KNNA for different values of K.

10). Figure 10 shows the KNNA results for estimating the UCS.

(**a**) (**b**) **Figure 8.** The FMP−ANN results for the optimum model (**a**) error reduction trend and (**b**) corelation

To apply the KNNA method to the data and determine the best K value, the KNNA was written in the form of a program in MATLAB software, which was run 310 times for K values from 1 to 30 programs; moreover, the amount of error was measured (Figure 9). Of the total data, 70% and 30% were used to train and test the model. The results showed that the lowest estimation error of the UCS was obtained at K = 2 (Figure 9). The error of this network for estimating the UCS with respect to the K values is equal to 0.11 (Figure

cofficient between measured and predicted UCS.

*3.8. The KNNA Results* 

**Figure 9.** Obtained RMSE for forecasting UCS based on the KNNA for different values of K. **Figure 9.** Obtained RMSE for forecasting UCS based on the KNNA for different values of K.

**Figure 10.** Accuracy of predicted UCS using KNNA: (**a**) correlation coefficient for test data and (**b**) error histogram for all data. **Figure 10.** Accuracy of predicted UCS using KNNA: (**a**) correlation coefficient for test data and (**b**) error histogram for all data.

#### *3.9. Results of SVM Method for Estimating UCS*

*3.9. Results of SVM Method for Estimating UCS* The SVM algorithm uses a set of mathematical functions that are named kernels [101]. The most important kernel functions for solving engineering problems are listed in Table 7. Normally, three radial basis kernel functions (RBFs), polynomial of degree, d, and linear, are used in the support vector machine, and the use of each of these functions with different parameters in the estimation of rock strength may lead to different results [60,101]. Therefore, it is necessary to evaluate the efficiency and accuracy of each of these functions and to choose the appropriate kernel function in predicting resistance. These three kernel functions were also used in this research. It should be mentioned that the calculation process of SVM was performed based on coding in a MATLAB environment and that the parameters of the kernel functions were optimized using a trial and error The SVM algorithm uses a set of mathematical functions that are named kernels [101]. The most important kernel functions for solving engineering problems are listed in Table 7. Normally, three radial basis kernel functions (RBFs), polynomial of degree, d, and linear, are used in the support vector machine, and the use of each of these functions with different parameters in the estimation of rock strength may lead to different results [60,101]. Therefore, it is necessary to evaluate the efficiency and accuracy of each of these functions and to choose the appropriate kernel function in predicting resistance. These three kernel functions were also used in this research. It should be mentioned that the calculation process of SVM was performed based on coding in a MATLAB environment and that the parameters of the kernel functions were optimized using a trial and error process. The results of these investigations are presented in Table 8. It can be observed that, based on the statistical criteria, the accuracy of the kernel functions is as RBF> PK > LK. In this regard, Nguyen [102] investigated the performance of various kernel functions using the support

Polynomial kernel (PK)

Radial basis function (RBF)

**Table 8.** Evaluation of SVM model performance in UCS estimation using various kernel functions.

and stated that the radial basis function has the highest performance.

**Function Description Kernel Function Type** 

This kernel is widely used in image processing, where d is the degree of the polynomial.

This kernel is used for general purposes. It is used when there is no prior knowledge about the data. In >0condition, = 1/2ଶparameter is used.

(, ) = . - Linear kernel (LK)

**Optimal Values of Parameters Test Period Train Period** 

PK 1.72 280.01 4 - 12.12 0.08 0.97 1.87 2.86 0.07 0.98 1.84 2.81 RBF 0.02 - - 1.10 27 0.06 0.99 1.90 2.82 0.06 0.99 1.90 2.80

LK 0.45 - - - 0.90 0.09 0.96 1.83 2.84 0.09 0.97 1.81

**t d c RMSE R2 PI MAPE RMSE R2 PI MAPE** 

(, ) = (. + 1)ௗ

(, ) = ( − ฮ − ฮ

**Kernel Function**  ଶ )

**Table 7.** The most important kernel functions for solving engineering problems [102].

process. The results of these investigations are presented in Table 8. It can be observed that, based on the statistical criteria, the accuracy of the kernel functions is as RBF> PK > LK. In this regard, Nguyen [102] investigated the performance of various kernel functions

next section.

did not occur.

**Table 9.** Modeling features using ANFIS**.** 

vector machine method in estimating blast-induced ground vibration and stated that the radial basis function has the highest performance.

**Table 7.** The most important kernel functions for solving engineering problems [102].


**Table 8.** Evaluation of SVM model performance in UCS estimation using various kernel functions.


The error histogram and predicted and estimated UCS relationship with the optimal function (RBF function) are presented in Figure 11. function (RBF function) are presented in Figure 11.

**Figure 11.** Accuracy of predicted UCS using SVM-RBF: (**a**) correlation coefficient and (**b**) error histogram for all data. **Figure 11.** Accuracy of predicted UCS using SVM-RBF: (**a**) correlation coefficient and (**b**) error histogram for all data.

#### *3.10. Results of ANFIS Method for Estimating UCS*

*3.10. Results of ANFIS Method for Estimating UCS* As with other intelligent methods, to test and train models using ANFIS, 30% and 70% of the whole data were used, respectively. The method of combining regression error propagation with least squares was used to train the model using the ANFIS. Table 9 shows the modeling features using the ANFIS method. A comparison of the performance of the methods for forecasting UCS based on different criteria has been reported in the As with other intelligent methods, to test and train models using ANFIS, 30% and 70% of the whole data were used, respectively. The method of combining regression error propagation with least squares was used to train the model using the ANFIS. Table 9 shows the modeling features using the ANFIS method. A comparison of the performance of the methods for forecasting UCS based on different criteria has been reported in the next section.

Number of membership functions (MFs) 6

Figure 12 shows the error histogram and correlation coefficient of the ANFIS model

in the test stage. This method, as with other used intelligent methods, has high accuracy in UCS estimation. The results of the intelligent models for estimating UCS from the test data performed better than the training data; therefore, it can be argued that overfitting

**FIS Generation Method GENFIS4**  Influence radius 0.60 Number of epochs 500

Error goal 0.00

Input MF type Gauss MF Output MF type Linear

Type Sugeno Rules 4

gram for all data.


**Table 9.** Modeling features using ANFIS.

Figure 12 shows the error histogram and correlation coefficient of the ANFIS model in the test stage. This method, as with other used intelligent methods, has high accuracy in UCS estimation. The results of the intelligent models for estimating UCS from the test data performed better than the training data; therefore, it can be argued that overfitting did not occur. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 17 of 22

**Figure 12.** Accuracy of predicted UCS using ANFIS: (**a**) correlation coefficient and (**b**) error histo-**Figure 12.** Accuracy of predicted UCS using ANFIS: (**a**) correlation coefficient and (**b**) error histogram for all data.

#### *3.11. Evaluation of the Used Methods*

*3.11. Evaluation of the Used Methods*  Table 10 and Figure 13 show the accuracy of the used methods for forecasting the UCS. According to the statistical criteria (i.e., R2, MAPE%, RMSE, VAF, and PI), the SVM-RBF model displays greater precision than other models because the SVM uses the minimizing structural risk theorem and adapts the ability of the model to existing training data [103]. The number of input variables, number of samples, and training algorithm type also Table 10 and Figure 13 show the accuracy of the used methods for forecasting the UCS. According to the statistical criteria (i.e., R<sup>2</sup> , MAPE%, RMSE, VAF, and PI), the SVM-RBF model displays greater precision than other models because the SVM uses the minimizing structural risk theorem and adapts the ability of the model to existing training data [103]. The number of input variables, number of samples, and training algorithm type also affect the accuracy of the methods [16,104,105]. Based on the correlation coefficient, all methods (R<sup>2</sup> > 90%) have excellent accuracy for estimating UCS.

affect the accuracy of the methods [16,104,105]. Based on the correlation coefficient, all methods (R2 > 90%) have excellent accuracy for estimating UCS. **Table 10.** Accuracy of approaches for predicting UCS. Considering that all six intelligent methods showed very high accuracy in UCS estimation, the percentage difference of the average of all six intelligent methods with the measured value in the laboratory is equal to +0.28%. This amount of difference is less than 1% and indicates the high capability of intelligent methods for forecasting the UCS.

 (**a**) (**b**)

**Figure 13.** Measured values versus predicted UCS using (**a**) SVR−RBF, ANFIS, and RFA methods,

and (**b**) FMP−ANN, KNN, and GPR−SEK methods.

KNNA 8.44 0.97 0.11 97.25 1.83 GPR-SEK 6.63 0.98 0.09 97.45 1.86 FMP-ANN 4.66 0.99 0.24 98.36 1.73

**Approaches MAPE% R2 RMSE VAF% PI** 

RFA 9.27 0.98 0.09 97.63 1.87

gram for all data.

*3.11. Evaluation of the Used Methods* 


**Table 10.** Accuracy of approaches for predicting UCS. SVM-RBF 2.83 0.99 0.06 98.96 1.92

methods (R2 > 90%) have excellent accuracy for estimating UCS.

**Table 10.** Accuracy of approaches for predicting UCS.

(**a**) (**b**)

Table 10 and Figure 13 show the accuracy of the used methods for forecasting the

UCS. According to the statistical criteria (i.e., R2, MAPE%, RMSE, VAF, and PI), the SVM-RBF model displays greater precision than other models because the SVM uses the minimizing structural risk theorem and adapts the ability of the model to existing training data [103]. The number of input variables, number of samples, and training algorithm type also affect the accuracy of the methods [16,104,105]. Based on the correlation coefficient, all

**Approaches MAPE% R2 RMSE VAF% PI** 

RFA 9.27 0.98 0.09 97.63 1.87

**Figure 12.** Accuracy of predicted UCS using ANFIS: (**a**) correlation coefficient and (**b**) error histo-

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 17 of 22

**Figure 13.** Measured values versus predicted UCS using (**a**) SVR−RBF, ANFIS, and RFA methods, and (**b**) FMP−ANN, KNN, and GPR−SEK methods. **Figure 13.** Measured values versus predicted UCS using (**a**) SVR−RBF, ANFIS, and RFA methods, and (**b**) FMP−ANN, KNN, and GPR−SEK methods.

#### **4. Conclusions**

The UCS of rocks is a basic parameter necessary for assessing the construction of civil and mining structures, such as the stability of the mines and the bearing capacity of foundations. UCS estimation using core specimens is costly, difficult, and, in some cases, impossible. After assessing the geo-mechanical features of 55 samples of sandstone, limestone, and argillaceous limestone specimens, predictive models for estimating UCS were developed via intelligent and statistical approaches. The results showed that the carbonate and sandstone samples were classified as mudstone to grainstone and calclitharenite, respectively. The PWV, WW, porosity, density, and PLI were considered as model inputs for predicting UCS. Statistical analysis allowed the development of equations with high accuracy to estimate UCS. Among the assessed linear, polynomial, and radial basis kernel functions, the accuracy of the other models was lower than that of SVM-RBF in forecasting UCS. The SVM-RBF model revealed that the R<sup>2</sup> and PI values were 0.99 and 1.92, respectively. The R<sup>2</sup> values of 98%, 98%, 97%, 98%, and 99% for forecasting the UCS were achieved using ANFIS, RFA, KNNA, GPR, and FMP-ANN, respectively. The number of samples and input variables had a significant impact on the performance of the methods. When the number of samples was small, the SVM method was more accurate. The percentage difference of the average of all six intelligent methods with the measured value was less than 1%, which indicates the superior capability of the intelligent methods in forecasting UCS.

**Author Contributions:** X.Z.: methodology, software, data curation. F.M.A.A.: Conceptualization, investigation. T.A.S.G.: writing—original draft preparation, methodology. A.H.D.A.-K.: analysis of results, validation. A.I.: performing field investigations and collecting samples R.B.Y.S.: writing —original draft preparation, resources. M.L.N.: supervision, project administration, funding acquisition. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** The authors are fully aware and satisfied with the contents of the article.

**Data Availability Statement:** The data used in this study has been appropriately described in the manuscript.

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
