Machine Learning-Based Intelligent Prediction of Elastic Modulus of Rocks at Thar Coalfield

Abstract

Elastic modulus (E) is a key parameter in predicting the ability of a material to withstand pressure and plays a critical role in the design of rock engineering projects. E has broad applications in the stability of structures in mining, petroleum, geotechnical engineering, etc. E can be determined directly by conducting laboratory tests, which are time consuming, and require high-quality core samples and costly modern instruments. Thus, devising an indirect estimation method of E has promising prospects. In this study, six novel machine learning (ML)-based intelligent regression models, namely, light gradient boosting machine (LightGBM), support vector machine (SVM), Catboost, gradient boosted tree regressor (GBRT), random forest (RF), and extreme gradient boosting (XGBoost), were developed to predict the impacts of four input parameters, namely, wet density (ρ_wet) in gm/cm³, moisture (%), dry density (ρ_d) in gm/cm³, and Brazilian tensile strength (BTS) in MPa on output E (GPa). The associated strengths of every input and output were systematically measured employing a series of fundamental statistical investigation tools to categorize the most dominant and important input parameters. The actual dataset of E was split as 70% for the training and 30% for the testing for each model. In order to enhance the performance of each developed model, an iterative 5-fold cross-validation method was used. Therefore, based on the results of the study, the XGBoost model outperformed the other developed models with a higher accuracy, coefficient of determination (R² = 0.999), mean absolute error (MAE = 0.0015), mean square error (MSE = 0.0008), root mean square error (RMSE = 0.0089), and a20-index = 0.996 of the test data. In addition, GBRT and RF have also shown high accuracy in predicting E with R² values of 0.988 and 0.989, respectively, but they can be used conditionally. Based on sensitivity analysis, all parameters were positively correlated, while BTS was the most influential parameter in predicting E. Using an ML-based intelligent approach, this study was able to provide alternative elucidations for predicting E with appropriate accuracy and run time at Thar coalfield, Pakistan.

1. Introduction

Elastic modulus (E) is a key parameter in predicting the ability of a material to withstand pressure and plays a critical role in the design process of rock-related projects. E has broad applications in the stability of structures in mining, petroleum, geotechnical engineering, etc. Accurate estimation of deformation properties of rocks, such as E, is very important for the design process of any underground rock excavation project. Intelligent indirect techniques for designing and excavating underground structures make use of a limited amount of data for design, saving time and money while ensuring the stability of the structures. This study has economic and even social implications, which are integral elements of sustainability. Moreover, this paper aims to determine the stability of underground mine excavation, which may otherwise result in a disturbed overlying aquifer and earth surface profile, adversely affecting the environment. E provides insight into the magnitude and characteristics of the rock mass deformation due to changes in the stress field. Deformation and behavior of different types of rocks have been examined by different scholars [1,2,3,4]. Usually, there are two common methods, namely, direct (destructive) and indirect (non-destructive), to calculate the strength and deformation of rocks. Based on the principles suggested by ISRM (International Society for Rock Mechanics) and the ASTM (American Society for Testing Materials), direct evaluation of E in the laboratory is a complex, laborious, and costly process. Simultaneously, in the case of fragile, internally broken, thin, and highly foliated rocks, the preparation of a sample is very challenging [5]. Therefore, attention should be given to evaluate E indirectly by the use of rock index tests.

Several authors have developed prediction frameworks to overcome these limitations by using machine learning (ML)-based intelligent approaches such as multiple regression analysis (MRA), artificial neural network (ANN), and other ML methods [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]. Advances in ML have so far been driven by the development of new learning algorithms and theories, as well as by the continued explosion of online data and inexpensive computing [22]. Similarly, Waqas et al. used linear and nonlinear regression, regularization and ANFIS (using a neuro-fuzzy inference system) to predict the dynamic E of thermally treated sedimentary rocks [23]. Abdi et al. developed ANN and MRA (linear) models, including porosity (%), dry density (γd) (g/cm³), P-wave velocity (Vp) (km/s), and water absorption (Ab) (%) as input features to predict the rock E. According to their results, the ANN model showed high accuracy in predicting E compared to the MRA [10]. Ghasemi et al. evaluated the UCS and E of carbonate rocks by developing a model tree-based approach. According to their findings, the applied method revealed highly accurate results [24]. Shahani et al. developed a first-time XGBoost regression model in combination with MLR and ANN for predicting E of intact sedimentary rock and achieved high accuracy in their results [25]. Ceryan applied the minimax probability machine regression (MPMR), relevance vector machine (RVM), and generalized regression neural network (GRNN) models to predict the E of weathered igneous rocks [26]. Umrao et al. determined strength and E of heterogeneous sedimentary rocks using ANFIS based on porosity, Vp, and density. Thus, the proposed ANFIS models showed superb predictability [27]. Davarpanah et al. established robust correlations between static and dynamic deformation properties of different rock types by proposing linear and nonlinear relationships [28]. Aboutaleb et al. conducted non-destructive experiments with SRA (simple regression analysis), MRA, ANN, and SVR (support vector regression) and found that ANN and SVR models were more accurate in predicting dynamic E [29]. Mahmoud et al. employed an ANN model for predicting sandstone E. In that study, 409 datasets were used for training and 183 datasets were used for model testing. The established ANN model exposed highly accurate results (coefficient of determination (R²) = 0.999) and the lowest mean absolute percentage error ((AAPE) = 0.98) in predicting E [30]. Roy et al. used ANN, ANFIS, and multiple regression (MR) to predict the E of CO₂ saturated coals. Thus, ANN and ANFIS outperformed the MR models [31]. Armaghani et al. predicted E of 45 main range granite samples by applying the ANFIS model in comparison with MRA and ANN. Based on their results, ANFIS proved to be an ideal model against MRA and ANN [32]. Singh et al. proposed an ANFIS framework for predicting E of rocks [33]. Köken predicted the deformation properties of rocks, i.e., tangential E (E_ti) and tangential Poisson’s ratio (v_ti) of coal-bedded sandstones located in the Zonguldak Hard Coal Basin (ZHB), northwestern Turkey, using various statistical and soft computing methods such as different regression and ANN evaluations including the physicomechanical, mineralogical, and textural properties of the rocks. According to this analysis, the remarkable results were that the mineralogical characteristics of the rock have a significant influence on the deformation properties. In addition to comparative analysis, ANN was considered as a more effective tool than regression analysis in predicting E_ti and v_ti of coal-bed sandstones [34]. Yesiloglu-Gultekin et al. used the different ML-based regression models such as NLMR, ANN, and ANFIS, and 137 datasets using unit weight, porosity, and sonic velocity to indirectly determine E of basalt. Based on the results and comparisons of various performance matrices such as R², RMSE, VAF, and a20-index, ANN was successful in predicting E over NLMR and ANFIS [35]. Rashid et al. used non-destructive tests, i.e., MLR and ANN, to estimate the Q-factor and E for intact sandstone samples collected from the Salt Range region of Pakistan. The ANN model predicted Q-factor (R² = 0.86) and E (R² = 0.91) more accurately than MLR regression for Q-factor (R² = 0.30) and E (R² = 0.36) [36]. E was predicted using RF by Matin et al. For comparison, multivariate regression (MVR) and generalized regression neural network (GRNN) were used for the prediction of E. The input V_p-R_n was used for E. According to their results, RF yielded more satisfactory conclusions than MVR and GRNN [37]. Cao et al. used an extreme gradient boosting (XGBoost) integrated with the firefly algorithm (FA) model for predicting E. consequently, the proposed model was appropriate for predicting E [17]. Yang et al. developed the Bayesian model to predict the E of intact granite rocks; thus, the model performed with satisfactory predicted results [38]. Ren et al. developed several ML algorithms, namely, k-nearest neighbors (KNN), naive Bayes, RF, ANN, and SVM, to predict rock compressive strength by ANN and SVM with high accuracy [39]. Ge et al. determined rock joint shear failures using scanning and AI techniques. Thus, the developed SVM and BPNN were considered as sound determination methods [40]. Xu et al. developed several ML algorithms, namely, SVR, nearest neighbor regression (NNR), Bayesian ridge regression (BRR), RF, and gradient tree boosting regression (GTBR), to predict microparameters of rocks by RF with high accuracy [41].

Based on the above literature and the limitations of the conventional predictive methods, a single model has low robustness, cannot achieve ideal solutions for all complex situations, and its performance varies with the input features. Therefore, authors have endeavored to use ML-based intelligent models that integrate multiple models to overcome the drawbacks of individual models and play a key role in determining the accuracy of the corresponding data for tests performed in the laboratory. However, there are few studies in predicting E. In addition, there are no comprehensive studies on the selection and application of such models in E prediction. To address this gap, this study developed six models based on an intelligent prediction approach, namely, light gradient boosting machine (LightGBM), support vector machine (SVM), Catboost, gradient boosted tree regressor (GBRT), random forest (RF), and extreme gradient boosting (XGBoost) to predict E, including wet density (ρ_wet) in gm/cm³, moisture in %, dry density (ρ_d) in gm/cm³, and Brazilian tensile strength (BTS) in MPa as input features under intricate and unsteady engineering situations. Next, 70% of the actual dataset of 106 is used for training and 30% for testing each model. To enhance the performance of the developed models, a repetitive 5-fold cross-validation approach is used. Intelligent prediction of E of sedimentary rocks from Block-IX of Thar coalfield has been applied for the first time. To the best of the author’s knowledge, application of intelligent prediction techniques in this scenario is lacking. Figure 1 depicts a systematic ML-based intelligent approach for predicting E.

2. A Brief Summary of the Study Area

The Thar coalfield is located in Sindh Province of Pakistan and is the seventh largest coal mine around the world in terms of coal potential [42]. Thar coal is classified as 175.5 billion tons of lignite, which can be used for fuel and power generation. The Thar coalfield is distributed in twelve different blocks as shown in Figure 2. The Thar coalfield is enclosed by dune sand that spreads to a normal distance of 80 m and rests upon an essential stand in the eastern portion of the desert. The general stratigraphic arrangement in the Thar coalfield encompasses the Basement Compound, coal posture Bara Formation, alluvial deposits, and dune sand. For coal mining in the region, both open-pit and underground mining methods can be preferred. Particularly, Sindh Engro Coal Mining Company (SECMC) has fully developed Block-II of the twelve blocks using an open-pit mining method, whereas Block-1 is under a development stage by Sino-Sindh Resources Ltd. (SSRL) in partnership with China and the Sindh government of Pakistan. Block-IX has been recommended for the underground mining method. The thickness of the coal seam of Block-IX of the Thar coalfield is approximately 12 m, the dip angle is 0° to 7°, and the top-bottom plate is siltstone–claystone to claystone. Shahani et al. proposed the use of the mechanized longwall caving (LTCC) method at Block-IX of the Thar coalfield in Pakistan for the first time [42,43]. In addition, Shahani et al. developed various gradient boosting machine learning algorithms to predict UCS of sedimentary rocks of Block-IX of the Thar coalfield [44]. Similarly, correct determination of the mechanical properties of Block-IX of the Thar coalfield, particularly E, plays an important role in fully understanding the behavior of the roof and ground prior to mining operations.

3. Data Curation

In this research, 106 samples of soft sedimentary rocks, i.e., siltstone, claystone, and sandstone were collected from Block-IX of the Thar coalfield, as shown Figure 2, with the location map in the green. Then, the rock samples were prepared and partitioned according to the principles suggested by ISRM [45] and the ASTM [46] to maintain the same core size, and geological and geometric characteristics. In the laboratory of the Mining Engineering Department of Mehran University of Engineering and Technology (MUET), the experimental work was conducted on the studied rock samples to determine the physical and mechanical properties such as wet density (ρ_wet) in g/cm³, moisture (%), dry density (ρ_d) in g/cm³, Brazilian tensile strength (BTS) in (MPa), and elastic modulus (E) in (GPa). Figure 3 shows (a) collected core samples, (b) universal testing machine (UTM), (c) deformed core sample under compression for E test, and (d) deformed core sample for BTS test. The purpose of the UCS test was conducted on the standard core samples of NX size 54 mm in diameter with an applied load of 0.5 MPa/s using UTM according to the recommended ISRM standard to find the E of the rocks. Similarly, in order to find the tensile strength of the rock samples indirectly, we performed the Brazilian test using UTM. Figure 4 illustrates the statistical distribution of the input features and output in the original dataset used in this study. In Figure 4, the legend of boxplots can be explained as: ▭ 25~75%, ⌶ Range within 1.5 IQR, ─ Median line, and ○ Outliers.

In order to visualize the original dataset of E, the seaborn module in Python was employed in this study, and Figure 5 demonstrates the pairwise correlation matrix and distribution of different input features and output E. It can be seen that BTS is moderately correlated to the E, whereas ρ_wet and ρ_d are negatively correlated to the E. Moisture representation does not correlate with E. It is worth mentioning that each feature cannot be well correlated with E independently, so all features are evaluated together to predict E.

4. Developing ML-Based Intelligent Prediction Models

4.1. Light Gradient Boosting Machine

Light gradient boosting machine abbreviated as LightGBM, an open-source gradient boosting ML model from Microsoft, uses decision trees as the base training algorithm [47]. LightGBM puts continuous buckets of elemental values into separate bins with greater adeptness and a fast speed of training. It uses a histogram-based algorithm [48,49] to improve the learning phase, reduce consumption of memory, and integrate updated communication networks to enhance the regularity of training and is known as a parallel voting decision tree ML algorithm. The data for learning were partitioned into several trees, and local voting techniques were executed in each iteration to select top-k elements and gain globing voting techniques. As shown in Figure 6, LightGBM operates the leaf-wise approach to identify the leaf with the maximum splitter gain. LightGBM is best adopted for regression, classification, sorting, and several ML schemes. It builds a more complex tree than the level-wise distribution method through the leaf-wise distribution method, which can be considered as the main component of the execution algorithm with greater effectiveness. For all that, it can cause overfitting; however, by using the maximum depth element in LightGBM, it can be disabled.

LightGBM [47] is a widespread library for performing gradient boosting, with some modifications intended. The implementation of gradient boosting is mainly focused on algorithms for building a computational system. The library includes tenfold training hyperparameters to validate the implementation of the framework in different scenarios. The implementation of LightGBM also demonstrates advanced capabilities on CPUs and GPUs, which can work like gradient boosting with multifold integrations, comprising column randomization, bootstrap subsampling, and so on. The main features of LightGBM are gradient-based one-sided sampling and unique attribute bundling. Gradient-based one-sided sampling is a sub-sampling technique used to construct the base tree of learning data as an ensemble. In the AdaBoost ML algorithm, the purpose of this technique is to increase the significance of samples with greater likelihood that are connected with samples with higher gradients. When gradient-based one-sided sampling is executed, the base learner’s learning data are articulated based on the top portion of samples with greater gradients (a) plus the portion of arbitrary orders (b) recouped from samples with lower gradients. To compensate for changes in measurement propagation, samples from the lesser gradient class are organized together and weighted by (1 − x)/y, and at the same time, computing the data gain. In contrast, the unique attribute bundling technique accrues meager elements into an individual element. This can be ended in the absence of impeding any information when these elements do not contain a non-zero number of coincidences. Both mechanisms predict a gain in the complementary learning rate.

4.2. Support Vector Machine

In 1997, Vapnik et al. originally proposed support vector machines (SVMs), which are a type of supervised learning [50]. SVMs can be widely used for regression analysis and for classification using hyperplane classifiers. The ideal hyperplane enhances the boundary between the two classes in which the support vector is positioned [51]. The SVM utilizes a high-extent feature space to develop the forecast function by proposing kernel functions and Vapnik’s 𝜀-insensitive loss function [52].

For a dataset P = {($x_1$, $y_1$), ($x_2$, $y_2$)…($x_n$, $y_n$)}, where $x_i$ ∈ $R^n$ is the input and $y_i$ ∈ $R^n$ is the output, the SVM employs a kernel function to plot the nonlinear input data in a high-extent feature space, and attempts to discover the best hyperplane to disperse them. This permits the narration of the original input to the output by a linear regression function [53,54,55] characterized as follows in Equation (1).

$$f\left(x\right)=M_v·\phi \left(x\right)+l_b$$

(1)

where $\phi \left(x\right)$ shows the kernel function, $M_v$ and $l_b$ show the weight vector and bias term, respectively. In order to obtain $M_v$ and $l_b$, the cost function proposed by Cortes and Vapnik [56] is required to be reduced as follows in Equation (2).

$$cost function=\frac{1}{2}{M}_v^2+C {{\displaystyle \sum }}_{i=1}^k\left({\xi }_i^{-}+{\xi }_i^{+}\right)$$

(2)

$$\mathrm{Subject} \mathrm{to}: \left\{\begin{array}{c}{y}_i-\left(M_v·\phi \left(x_1\right)+l_b\right) \le {\epsilon }_0+{\xi }_i^{+}\\ \left(M_v·\phi \left(x_1\right)+l_b\right)-y_i \le {\epsilon }_0+{\xi }_i^{-} \\ {\xi }_i^{-},{\xi }_i^{+}\ge 0, i=1, 2, \dots .,n \end{array}\right.$$

When converted to the dual space using the Lagrange multiplier method, Equation (2) can be reduced to obtain the following solution in Equation (3).

$$f\left(x\right)={\displaystyle \sum }_{i=1}^n\left({\infty }_i-{\infty }_i^{\prime }\right)\phi \left(x_i,x_j\right)+l_b$$

(3)

where ${\infty }_i$ and ${\infty }_i^{\prime }$ are Lagrange multipliers with 0 ≤ ${\infty }_i$ and ${\infty }_i^{\prime }$≤ C, though $\phi$ $\left(x_i,x_j\right)$ is the kernel function. The choice of the latter is important to the accomplishment of SVR. A large number of kernel functions have been studied in SVM, such as linear, polynomial, sigmoid, gaussian, radial basis, and exponential radial basis [54]. Figure 7 illustrates the basic structure of the SVM model.

4.3. Catboost

Catboost is a gradient boosting algorithm recently inherited by Dorogush et al. [57]. Catboost solves the complex problem of regression and classification simultaneously and is publicly available in an open-source multi-platform gradient boosting library [57,58]. In the Catboost algorithm, the decision tree is used as the underlying weak learner and the gradient boosting is successively fitted to the decision tree. To improve the implementation of the Catboost algorithm and to avoid overfitting, an inconsistent arrangement of gradient learning information is used [57].

The purpose of the Catboost algorithm is to reduce the predictive movement that occurs in the learning phase. Propagation movement is the deletion of F(y)|(yi) in the case that yi is the learning sample, related to F(y)|(yi) of the test sample y. In the learning phase, gradient boosting uses the same samples to compute the gradient and the model for lowering that gradient. The Catboost’s concept is to initiate j... n, the framework underlying the repetition of separate P enhancing. The mth recurrence’s ith framework is learned from the permutation of the initial ith sample and is suitable for computing the j + 1 sample’s gradient of the p + 1 recurrence. Subsequently, in order not to be limited by the starting arbitrary permutation, the technique employs an arbitrary permutation of s reciprocals. Each repetition constructs a distinguished framework that achieves all combinations and frameworks. Symmetric trees are used as the basis of the framework. These trees are prolonged by using the same partitioning criterion so that all leaf nodes grow level-wise.

In the Catboost algorithm, the mechanism proposed is to compute the identical up-to-date characters as the ones imitated when the network was built. Thus, for any specified samples’ permutation, data samples <i are utilized to computing the character values for each sample i. Then, different combinations are implemented, and the character values obtained for each sample are averaged. Catboost is a large-scale comprehensive library consisting of several elements such as GPU learning, standard boosting, and including fivefold hyperparameter optimization to amend to various practical examination situations. Standard gradient boosting is to be considered as part of the Catboost algorithm also. Figure 8 shows an explanation of the Catboost algorithm.

It is very important to note that the Catboost algorithm’s training ability is managed by its framework hyperparameters, i.e., iterations number, rate of learning, maximum depth, etc. Determining the optimal hyperparameters of a model is a challenging, laborious, and tedious task that depends on the user’s skills and expertise.

4.4. Gradient Boosted Regressor Tree

The gradient boosted regressor tree (GBRT) regression integrates the weak learner, i.e., the learner algorithms with average performance compared to random algorithms into a robust learner with an iterative method [59]. Contradictory of the bagging method, the boosted algorithm continuously generates the underlying framework. The soundness of the predictive framework is improved by prioritizing this hard-to-evaluate learning information to generate several frameworks in a series. In the boosting algorithm, underlying frameworks that were previously not suitable for estimation are frequently established in the training dataset compared to those frameworks that have been accurately evaluated. Each complementary underlying framework is designed to correct inaccuracies arising from its prior underlying framework. The occurrence of the boosting mechanism comes from reaction of Schapire’s feedback to Kearns’ investigation [60,61] (Kearns): Is the aggregation of weak learners a substitute for distinguishing strong learners? Weak learners are explained as the algorithms that work well compared to random approximation; strong underlying frameworks are a more realistic classification or regression algorithms that are incoherent with their effective counterparts to the problem. The reaction to such inquiry is extremely noteworthy. Assessments of weak frameworks tend to be less challenging than strong frameworks, and Schapire’s establishment of a “yes” response to Kearns’ inquiry, as evidenced by the combination of several weak frameworks into an upgraded and independent sound framework. The key dissimilarity between the boosting and bagging mechanism is that in the boosting approach, the training dataset is analytically investigated in order to predict the most appropriate instructions for each subsequent framework. In every phase of training, the modified propagation is dependent on the inaccuracies raised by the previous frameworks. On the contrary, in the bagging mechanism, every trial is constantly specified to produce a training dataset, and for the boosting mechanism, the vagueness of specifying an independent trial is conflicting. Trials which were erroneously assessed were more likely to set higher weights. Accordingly, each newly evolved framework underscores trials that are inaccurately assessed by subsequent frameworks.

Boosting assembles the auxiliary frameworks that decrease a specific loss function averaged over the learning dataset, i.e., the MAE or the MSE. The loss function computes the total number of predicted values that differ from the investigated values. The advanced staged modeling method is one of the assessed elucidations to the problem. This modeling method continuously attaches the new underlying framework without substituting the coefficients and specifications of the previously connected model. Referring to the regression model, the boosting mechanism is a “function gradient descent” configuration. Functional gradient descent is an optimization mechanism that minimizes the loss function by connecting the underlying framework to each stage to reduce the loss function by a certain amount. Figure 9 demonstrates the schematic diagram of GBTR (after [62,63]).

Friedman recommended improving the gradient boosting regression model by using pre-established regression trees for the underlying framework. The improved framework amplifies the performance of Friedman’s model [64]. For predicting E, the improved version of gradient boosted regression was used. Considering that the number of leaves is l, each tree divides the input space into l independent territory $T_{1p}$, $T_{2p}$……… $T_{lp}$ and a perpetual value $k_{lp}$ is predicted for territory $T_{lp}$. Equation (4) represents the gradient boosting regression tree as follows:

$$f_p\left(a\right)={{\displaystyle \sum }}_{l=1}^L{k}_{lp}F (a \in T_{lp})$$

(4)

where $F\left(a \in T_{lp}\right)$ = $\left\{\begin{array}{c}1, if a \in T_p \\ 0, otherwise \end{array}\right.$.

By using a regression tree to recover $f_p \left(a\right)$ in the generic gradient boosting mechanism, the framework gradient descent stage size and updating equation are given by Equations (5) and (6), respectively.

$$f_p\left(a\right)=f_{p-1}\left(a\right)+{\rho }_p{g}_p \left(a\right)$$

(5)

$${\rho }_p=argmi{n}_{\rho } {{\displaystyle \sum }}_{l=1}^{L}M(b_i,f_{p-1}\left(a_i\right)+\rho g_p\left(a_i\right))$$

(6)

Hence, Equations (5) and (6) fit as Equations (7) and (8).

$$f_p \left(a\right)=f_{p-1}\left(a\right)+{{\displaystyle \sum }}_{l=1}^L{\rho }_p{k}_{lp}F( a \in T_{lp})$$

(7)

$${\rho }_p=argmi{n}_{\rho } {{\displaystyle \sum }}_{l=1}^{L}M(b_i,f_{p-1}\left(a_i\right)+{{\displaystyle \sum }}_{l=1}^L{\rho }_p{k}_{lp}F( a \in T_{lp}))$$

(8)

By applying a discrete ideal ${\rho }_{lp}$ for each territory $T_{lp}$ $k_{lp}$ should be separated. The simplified framework Equations (9) and (10) are given by

$$f_p \left(a\right)=f_{p-1}\left(a\right)+{{\displaystyle \sum }}_{l=1}^L{\rho }_{p}F( a \in T_{lp})$$

(9)

$${\rho }_p=argmi{n}_{\rho } {{\displaystyle \sum }}_{l=1}^{L}M(b_i,f_{p-1}\left(a_i\right)+{{\displaystyle \sum }}_{l=1}^L{\rho }_{p}F( a \in T_{lp}))$$

(10)

The overfitting of the framework can be limited by managing the number of iterations of gradient boosting or more competently by evaluating the degree of benefit of each tree by J ∈ (0, 1). Thus, the simplified model is given by Equation (11).

$$f_p \left(a\right)=f_{p-1}\left(a\right)+L{{\displaystyle \sum }}_{l=1}^L{\rho }_{p}F( a \in T_{lp})$$

(11)

4.5. Random Forest

In 2001, Breiman originally proposed random forest (RF), a type of ensemble machine learning algorithm [65]. RF can be widely used for regression analysis and classification. RF is a state-of-the-art approach to bootstrap aggregating or bagging. RF perceives a one-of-a-kind relationship of model embodiment and predictive accuracy among alternative recognized AI computing [66].

To calculate the performance of the model, RF of 100 trees with a range of default settings was chosen for the study. Figure 10 shows the basic structure of the RF model.

4.6. Extreme Gradient Boosting

Extreme gradient boosting (XGBoost) is an important type of ensemble learning algorithm in ML approaches [67]. XGBoost consists of usual regression and classification trees with the addition of analytical boosting methods. The boosting method improves the accurate framework assessment by constructing different trees as alternatives to develop an addressed tree and then connecting them to estimate a systematic predictive algorithm [68]. It instigates the tree by consecutively holding the residuals of the historical trees as the effect of the resultant tree. Because of this, the resultant tree builds a full prediction by generating errors in the past tree. In the loss function reduction stage, the consecutive framework structural relationship can be subdivided into gradient descent types, which advances the forecast by connecting a supplementary tree at every stage to lessen the depletion [69]. Tree development stops at the time of the most unprecedented tree’s predetermined number is obtained, or at the time of the training stage error when it cannot be amplified to the predicted number of sequential trees. By attaching an arbitrary survey, the performance timeliness and estimation accuracy of gradient boosting can be greatly improved. In particular, for all symmetric trees, a random subsample of training information from the whole training dataset is considered, without substitutions. This arbitrarily described subsample replaces the whole sample, which is later used to adapt the tree and is identified as an improved framework. XGBoost is a state-of-the-art rearranged gradient boosting ML algorithm that manages and implements the latest prediction demonstrations [49]. The loss function’s subsequent assessment is used in the XGBoost and is fast and rapidly matched to the usual gradient boosting algorithms. XGBoost has widely been used to mine the features of gene coupling. Figure 11 shows the general structure of XGBoost models.

Consider $\overline{u_i}$ is the predicted outcome of the $i$th data, where the feature vector is $V_i$; E denotes the number of estimators and for every estimator $f_k$ (with k from 1 to E) analogous to the analysis of a single tree. $u_i^0$ describes the initial hypothesis and is the mean of the examined features in the data for learning. Equation (12) executes different extension functions to predict the results.

$$\overline{u_i}=u_i^0+\eta {{\displaystyle \sum }}_{K=1}^{\mathrm{E}}{f}_k(V_i)$$

(12)

Additionally, the η parameter is the learning rate, which is contiguous to the implementation of the improved model to enhance the model, perform rhythmically when connecting the latest trees, and combat overfitting.

With respect to Equation (12), the kth character is attached to the model in the kth state, the kth prediction $u_i^{-k}$ is realized by the prediction $u_i^{-\left(k-1\right)}$ of the previous state, and the additional kth character augmentation $f_k$ is described as in Equation (13).

$$u_i^{-k}=u_i^{-\left(k-1\right)}+\eta f_k$$

(13)

where $f_k$ denotes the weight of the leaves established by decreasing the kth tree’s objective function signified by Equation (14).

$$obj=\gamma N+{{\displaystyle \sum }}_{a=1}^N[ U_a{\omega }_a+\frac{1}{2} \left(V_a+\lambda \right) {\omega }_a^2$$

(14)

where N represents the kth tree’s leaves, ${\omega }_a$ denotes the leaves’ weights from 1 to N. $\gamma$ and $\lambda$ are the regularity attributes to achieve anatomical consistency to avoid the model’s overfitting. The $V_a$ and $U_a$ parameters are the sets of whole data connected with the prior data’s leaf and the gradient of the posterior loss function, respectively.

In the process of building the kth tree, individual leaves are partitioned into a different number of leaves. Equation (15) represents the dissection using the gain parameter. Consider that $U_R$ and $V_R$ describe inter-reliant right leaves and $U_L$ and $V_L$ are inter-reliant left leaves for divergence. At this point, the gain parameter is near to zero, which is traditionally considered as the benchmark for divergence. γ and λ are uniform features that affect the gain features, such as the gain parameter being reduced by a higher regularization parameter and thus avoiding leaf convolution. However, it reduces the adaptability of the framework to the training dataset.

$$gain=\frac{1}{2} \left[\frac{U_L^2}{V_L+\lambda }+\frac{U_R^2}{V_R+\lambda }+\frac{ {\left(U_L+U_R\right)}^2}{V_L+V_R+\lambda }\right]$$

(15)

XGBoost is a broadly adopted ML algorithm that brings together articulation and logical achievements of gradient boosting ML algorithms. A numerical value causes the problem with the prediction regression model. XGBoost can be accomplished in a timely manner in a probabilistic regression framework. The ensemble is built from a decision tree model. Ensembles are continuously connected trees that can be adjusted to predict imprecise models. These ensemble-type ML methods are called boosting. These frameworks are built by executing any random gradient descent optimization method with a unique loss function. When the model is executed, the gradient loss function is reduced, and so this technique is known as “gradient boosting”. Compared with LightGBM, SVM, Catboost, GBRT, and RF, the XGBoost model performed well on the E dataset with the identical parameters n_splits = 5, n_repeats = 3, and random_state = 1 (all remaining parameters were used as default parameters in Python). Although in this study, gbm_param_grid was further implemented in order to further improve the XGBoost model’s performance.

4.7. K-Fold Cross-Validation

The K-fold cross-validation is a technique employed to regulate the hyperparameters [70]. The technique accepts a search within a demarcated hyperparameters’ range and describes the predictable outcomes leading to the best outcome for calculation criteria such as R², MAE, MSE and RMSE. In the scikit-learn Python programing language, K-fold Cross-Validation has been implemented to handle this approach. This method simply calculates the score of CV for all hyperparameters integrated with a specific range. In this study, a 5-fold iterated arbitrary arrangement practice was integrated into the CV command as illustrated in Figure 12. GridSearchCV() permits not only the calculation of the anticipated hyperparameters, but also the evaluation of the metric values to their anticipated results.

4.8. Models Performance Evaluation

To accurately and approximately evaluate the performance of ML-based intelligent models, different authors have used different estimation criteria, namely, coefficient of determination (R²) [71], mean absolute error (MAE), mean square error (MSE), root mean square error (RMSE) [72], and a20-index [71]. Performance criteria are the main metrics used to assist in the highly accurate model evaluation, with the highest R², minimum MAE, MSE, RMSE, and appropriate a20-index. The following performance indices are employed to evaluate the performance of each model in E prediction.

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n\left({\mathrm{E}}_i-{\widehat{\mathrm{E}}}_i\right)}{{\mathrm{E}}_i-{\overline{\mathrm{E}}}_i}$$

(16)

$$\mathrm{MAE}=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^n|{\mathrm{E}}_i-{\widehat{\mathrm{E}}}_i|$$

(17)

$$\mathrm{MSE}=\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\mathrm{E}}_i-{\widehat{\mathrm{E}}}_i\right)}^2}{N}$$

(18)

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\mathrm{E}}_i-{\widehat{\mathrm{E}}}_i\right)}^2}{N}}$$

(19)

$$a20-index=\frac{m20}{N}$$

(20)

where ${\overline{\mathrm{E}}}_i$ and $\widehat{{\mathrm{E}}_i}$ are the mean values of the measured and predicted values of E, and ${\mathrm{E}}_i$ are measured values of E, respectively. m20 represents the datasets with a value of rate original/estimated values between 0.80 and 1.20 and N denotes the number of datasets.

5. Analysis of Results and Discussion

This study aims to examine the capability of various ML-based intelligent prediction models, namely, LightGBM, SVM, Catboost, GBRT, RF, and XGBoost, for predicting a substantial E using Python programming. In order to propose the most suitable prediction model to predict E, the selection of appropriate input features can be considered as one of the most important tasks. In this study, wet density (ρ_wet) in gm/cm³, moisture (%), dry density (ρ_d) in gm/cm³, and Brazilian tensile strength (BTS) in (MPa) were taken as the input features for all developed models.

Later, the measured and predicted output values were organized and plotted to facilitate the performance analysis and correlation of the developed models. The final output was examined using various analytical indices such as R², MAE, MSE, RMSE, and a20-index as performance criteria to analyze and compare the anticipated models and to evaluate the ideal model in terms of data prediction. The106 data points of the overall dataset were allocated as 70% (74 data points) for training and 30% (32 data points) for testing the model.

Figure 13 illustrates the scatter plots of predicted E of the test data by LightGBM, SVM, Catboost, GBRT, RF, and XGBoost models. The R² value of each model is determined according to the test prediction. The R² value of LightGBM, SVM, Catboost, GBRT, RF, and XGBoost models is 0.281, 0.32, 0.577, 0.988, 0.989, and 0.999, respectively.

Furthermore, to further understand the performance of the predicted E, it will be interesting to study the prediction rules of six developed ML-based intelligent models due to the wide dispersion of the range of values of E in the established dataset of the test data. The residuals (GPa) and percentage errors (%) of six models were utilized to view the predicting results. The residuals allow for the observation of the contrast between the predicted E and the measured E for each data point, and the percentage error shows the percentage by which the predicted E surpasses the measured E. They are expressed as Equations (21) and (22).

$$r={\mathrm{E}}_m-{\mathrm{E}}_p$$

(21)

$$p_{error}=\raisebox{1ex}{$r$}\!\left/ \!\raisebox{-1ex}{${\mathrm{E}}_m$}\right.\ast 100$$

(22)

where r = residual in GPa; E_m and E_p are the measured and predicted E, respectively; and p_error is the percentage error in %.

In Figure 14, the residuals indicate a direct relationship with the E, since the corresponding residuals can increase as the E increases. In contrast, in Figure 15, the percentage error shows an inverse relationship with the E, because it decreases as the E increases. Some models show negative residuals and percentage errors for smaller E measures and positive values for larger E measures. It revealed that these ML-based intelligent models seem to tend to have a predicted E higher than the measured E when the measured E was small, and tend to have a predicted E smaller than the measured E when the measured E is higher.

Table 1. Performance indices of the developed ML-based intelligent models in this study.

Model	Training					Testing
	R2	MAE	MSE	RMSE	a20-Index	R2	MAE	MSE	RMSE	a20-Index
LightGBM	0.496	0.1272	0.0470	0.2168	0.836	0.281	0.1340	0.0269	0.1640	1.012
SVM	0.324	0.1461	0.0805	0.2837	1.07	0.32	0.1031	0.0259	0.1609	1.22
Catboost	0.891	0.1091	0.0113	0.1069	1.04	0.577	0.218	0.0948	0.3101	0.86
GBRT	0.995	0.0162	0.0004	0.0200	0.96	0.988	0.0147	0.0003	0.0173	0.962
RF	0.991	0.0102	0.0018	0.0424	0.99	0.989	0.0284	0.0016	0.0400	0.943
XGBoost	0.999	0.0008	0.0004	0.0089	0.914	0.999	0.0015	0.0008	0.0089	0.996

exhibits the performance indices of the developed LightGBM, SVM, Catboost, GBRT, RF, and XGBoost models computed by Equation (16) to Equation (20). In this study, according to the proposed LightGBM, SVM, Catboost, GBRT, RF, and XGBoost models, the XGBoost outperformed when at the test data with R² of 0.999, MAE of 0.0015, MSE of 0.0008, RMSE of 0.0089, and a20-index of 0.996, for E prediction. In addition, GBRT and RF have also shown high accuracy and achieved second place next to XGBoost in predicting E, but they can be used conditionally. Therefore, XGBoost is an applicable ML-based intelligent approach that can be applied to accurately predict E, as shown in Figure 16.

The Taylor diagram explains a brief qualitative depiction of the best fit of the model to standard deviations and correlations. The expression for the Taylor diagram is given in Equation (23) [73].

$$R=\frac{\frac{1}{P} {{\displaystyle \sum }}_p^P\left(r_n-\overline{r}\right)\left(f_n-\overline{f}\right)}{{\sigma }_r{\sigma }_f}$$

(23)

where R denotes a correlation, P shows the discrete point number, $r_n$ and $f_n$ show two vectors, ${\sigma }_r$ and ${\sigma }_f$ show r and f standard deviation, and $r_n$ and $f_n$ illustrate the mean value of vectors $r_n$ and $f_n$, respectively.

Figure 17 represents the correlation between the predicted E and the measured E for the LightGBM, SVM, Catboost, GBRT, RF, and XGBoost models from Figure 16 in terms of standard deviation (STD), RMSE, and R². Based on the consequences, the XGBoost model was highly correlated with measured E than the other models developed in this study in predicting E.

Furthermore, the standard deviation (STD) of XGBoost was closest to the measured STD. Thus, compared to the existing published literature [8,74,75,76], XGBoost exhibits high accuracy and proved to be a highly accurate model for predicting E. The STD of GBTR and RF was also close to the measured STD but indicates the lowest R² values. Meanwhile, LightGBM, SVM, and Catboost showed the least correlation and were far from the measured STD.

6. Sensitivity Analysis

It is very important to correctly evaluate the essential parameters that have a large impact on the E of rock, which is undoubtedly a challenge in the design of rock structures. Thus, in this study, the cosine amplitude method [77,78] was adopted to investigate the relative impact of the inputs over the output. The general formulation of the adopted method is shown in Equation (24).

$$XGBoos{t}_{ij}=\frac{{{\displaystyle \sum }}_{k=1}^n({\mathrm{E}}_{ik}{\mathrm{E}}_{jk})}{\sqrt{{{\displaystyle \sum }}_{k=1}^n{\mathrm{E}}_{ik}^2{{\displaystyle \sum }}_{k=1}^n{\mathrm{E}}_{jk}^2}}$$

(24)

where ${\mathrm{E}}_i$ and ${\mathrm{E}}_j$ are the input and output values, respectively, and n is the number of datasets in the test phase. Finally, the range of $XGBoos{t}_{ij}$ is between 0 and 1, additionally proving the precision between each variable and the target. According to Equation (24), if $XGBoos{t}_{ij}$ of any parameter has a value of 0, it shows that there is no significant relationship between that parameter and the target. On the contrary, when $XGBoos{t}_{ij}$ is equal to 1 or nearly 1, it can be considered a significant relationship that has a large effect on the E of the rock.

Because of the high accuracy of the XGBoost model in predicting E, only a sensitivity analysis was performed on it at the testing level. Figure 18 shows the relationship between each input parameter of the developed model and output. Therefore, it can be seen from the figure that all parameters are positively correlated, while BTS is the most influential parameter in predicting E. The feature importance of each input parameter is given as ρ_wet = 0.0321, moisture = 0.0293, ρ_d = 0.0326, and BTS = 0.0334.

7. Conclusions

Elastic modulus (E) plays a key role in the designing of any rock engineering project. Therefore, an accurate determination of E is a prerequisite. In this study, six novel ML-based intelligent models, namely, LightGBM, SVM, Catboost, GBRT, RF, and XGBoost, were developed to predict E, including four input features, namely, ρ_wet, moisture, ρ_d, and BTS. To avoid overfitting of these models, the original dataset was distributed into 70% for the training and 30% for the testing of 106 data points. The study concludes that the XGBoost model performed more accurately than the other developed models, such as LightGBM, SVM, Catboost, GBRT, and RF, in predicting E with R², MAE, MSE, RMSE, and a20-index values of 0.999, 0.0015, 0.0008, 0.0089, and 0.996 of the test data, respectively. By employing the ML-based intelligent approach, this study was able to provide alternative elucidations for predicting E with appropriate accuracy and run time.

In future rock engineering projects, it is highly recommended to undertake proper field investigations prior to decision making. The XGBoost ML-based intelligent model performed well to predict the E. The conclusions of GBRT and RF are also applicable for the prediction of E; however, these methods can be used conditionally. Thus, for a large-scale study, this study recommends an adequate dataset to overcome the above limitation. In order to undertake other projects, the model proposed in this study should be considered as a foundation and the result should be reanalyzed, reevaluated, and even re-addressed.