Optimized Random Forest Models for Rock Mass Classification in Tunnel Construction

Abstract

The accurate prediction of rock mass quality ahead of the tunnel face is crucial for optimizing tunnel construction strategies, enhancing safety, and reducing geological risks. This study developed three hybrid models using random forest (RF) optimized by moth-flame optimization (MFO), gray wolf optimizer (GWO), and Bayesian optimization (BO) algorithms to classify the surrounding rock in real time during tunnel boring machine (TBM) operations. A dataset with 544 TBM tunneling samples included key parameters such as thrust force per cutter (TFC), revolutions per minute (RPM), penetration rate (PR), advance rate (AR), penetration per revolution (PRev), and field penetration index (FPI), with rock classification based on the Rock Mass Rating (RMR) method. To address the class imbalance, the Borderline Synthetic Minority Over-Sampling Technique was applied. Performance assessments revealed the MFO-RF model’s superior performance, with training and testing accuracies of 0.992 and 0.927, respectively, and key predictors identified as PR, AR, and RPM. Additional validation using 91 data sets confirmed the reliability of the MFO-RF model on unseen data, achieving an accuracy of 0.879. A graphical user interface was also developed, enabling field engineers and technicians to make instant and reliable rock classification predictions, greatly supporting safe tunnel construction and operational efficiency. These models contribute valuable tools for real-time, data-driven decision-making in tunneling projects.

1. Introduction

The utilization of tunnel boring machines (TBMs) in rock tunneling has become widespread due to increased investment in transportation infrastructure, such as highways and railways, and continuous advancements in construction technology. This popularity is primarily attributed to their advantages, including high construction quality, minimal surrounding rock disturbance, low labor intensity, and cost-effectiveness [1]. However, TBMs are susceptible to variations in geological conditions, and unforeseen geological circumstances can impede tunneling efficiency and potentially lead to significant engineering incidents [2]. For example, a sudden increase in the strength of the rock mass without timely adjustment of the TBM tunneling parameters can result in the excessive wear and tear of machine components. This, in turn, may adversely affect the machine’s service life and prolong downtime for maintenance. Conversely, if the strength of the rock mass suddenly decreases and the excavation parameters are not adjusted in time, it can lead to the collapsing of the weak rock mass, causing TBM jamming accidents in serious cases [3]. Therefore, timely and accurate assessment of the surrounding rock quality ahead of the tunnel face is crucial for safe tunnel construction.

Currently, traditional methods for rock mass classification have undergone significant development, including the Rock Mass Rating (RMR), Q, Geological Strength Index (GSI), Rock Mass Index (RMi), Barton’s Q (BQ), and Hardness Coefficient (HC), among others. These methods find wide applications in mines, tunnels, and other underground projects [4,5,6]. However, these rock classification methods generally rely on the physical and mechanical properties of the surrounding rock, geological structure, and hydrological conditions. They often require the acquisition of numerous parameters, many of which necessitate on-site data collection and indoor testing. Consequently, the time required for rock mass classification is prolonged, significantly impacting tunnel construction efficiency. Moreover, these tests consume both labor and resources, contributing to increased construction costs.

Presently, various geological investigation techniques are employed in the geological assessment in front of the tunnel face, including drilling and geology analysis, nondestructive geophysical exploration, and integrated exploration and interpretation [7]. While these methods can obtain reliable, comprehensive, and accurate geological information, the exploration process requires the TBM to stop working, which fails to meet the demands of rapid TBM excavation. Additionally, the confined space between the tunnel face and the TBM cutter presents challenges for installing geological forecasting equipment. Common methods of surrounding rock classification and geological investigation often overlook the influence of engineering on the evaluation of rock quality. In TBM construction, the dynamic interplay between the TBM system and the surrounding geological context significantly affects the construction process. The disruptive impact of the TBM on the surrounding rock, coupled with delays in identifying surrounding rock classification, can lead to geohazards and safety accidents. Real-time prediction of surrounding rock classes using TBM tunneling parameters can help mitigate such problems. Ji et al. [8] attempted to establish regression equations for predicting rock mass classification solely based on TBM tunneling parameters but encountered poor prediction performance.

The integration of artificial intelligence (AI) and big data is gaining traction in tunnel engineering due to their rapid advancements [9,10,11]. Unlike conventional empirical formulations, AI methods excel at capturing complex nonlinear relationships and the effects of multivariate factors, yielding more accurate predictions, especially in complex geological and engineering environments. Therefore, the application of AI in surrounding rock classification during tunnel construction holds promising potential. For example, Zhang et al. [12] applied data compression techniques to process extensive TBM operational data, utilizing the K-means++ clustering algorithm to identify potential surrounding rock classes. Subsequently, a high-accuracy predictive model for classifying the surrounding rock was developed using a support vector machine (SVM). Liu et al. [13] constructed an ensemble learning model for surrounding rock classification prediction based on the Classification and Regression Tree (CART) and AdaBoost algorithms, utilizing TBM tunneling parameters and surrounding rock classifications acquired from the hydropower classification method. The model’s superiority was confirmed through comparisons with CART, SVM, artificial neural network (ANN), and k-nearest neighbor (KNN) models. Hou et al. [14] employed the stacking technique of ensemble learning to forecast surrounding rock classification using TBM operational data from the Songhua River diversion tunnel project in China. This method, when compared to other individual classifiers, exhibited superior adaptability and performance, especially for small and unbalanced datasets. Yang et al. [15] simplified rock mass classification into two categories, incompetent rock mass and competent rock mass, improving the prediction performance of random forest (RF) and convolutional neural networks using the Bayes boosting method. Sebbeh-Newton et al. [16] addressed imbalanced rock mass classes by employing oversampling techniques and established a predictive model for surrounding rock classification using TBM operating parameters and surrounding rock classes from the Japanese Highway Classification System. The ensemble learning algorithm, with RF being a popular choice, demonstrated enhanced stability and adaptability. However, the predictive ability of AI models is still limited by hyperparameters, and improper selection and adjustment may affect optimal model performance. For example, in the case of RF, if hyperparameters such as the number of trees or depth are set too low, the model may fail to capture the complex relationships in the data, leading to underfitting. On the other hand, if these parameters are set too high, the model could become overly complex and may start memorizing the noise in the training data, resulting in overfitting. Both underfitting and overfitting can negatively impact the accuracy of the model’s predictions. To enhance model performance, various optimization algorithms, particularly metaheuristic algorithms, have been proposed for hyperparameter tuning [17,18]. The versatility and robustness of metaheuristic algorithms have led to their widespread acceptance in the engineering field. For instance, He et al. [18] utilized the tunicate swarm algorithm (TSA), the whale optimization algorithm (WOA), and the gray wolf optimizer (GWO) to optimize the random forest (RF) model, achieving precise prediction of overbreak phenomena resulting from tunnel blasting. Similarly, Zhou et al. [19] enhanced the predictive performance of RF through the utilization of the moth-flame optimization (MFO), GWO, and the multi-verse optimizer (MVO) algorithm for stability assessment in underground entry-type excavations. While metaheuristic algorithms and RF are increasingly applied in geotechnical engineering, their application in the classification of surrounding rock for TBM construction has received limited research attention.

Based on the above analysis, this paper constructed a database using TBM tunneling parameters derived from the Pahang Selangor Raw Water Transfer (PSRWT) tunnel and the surrounding rock classification obtained through the Rock Mass Rating (RMR) method. Two metaheuristic algorithms were employed to optimize RF and create real-time predictive models for classifying rock mass in TBM excavation. To demonstrate the superior performance of the metaheuristic algorithms, we also compared them with the Bayesian optimization algorithm. It is worth noting that, at the current stage, there has not been an exhaustive study on hybrid model development for the real-time prediction of surrounding rock classes in TBM construction. Therefore, this study fills the gap in this field by developing three hybrid random forest models. The subsequent section outlines the paper’s structure: Section 2 offers a comprehensive explanation of the theoretical foundation of the models and delineates the construction process. Section 3 introduces the data sources used in this study and analyzes and preprocesses the data. In Section 4, the training process and evaluation metrics of the hybrid model are described. Section 5 showcases the outcomes of the model evaluation alongside the sensitivity analysis. Section 6 presents the engineering validation results of the models and shows the corresponding graphical user interfaces. Section 7 deliberates on the limitations of this study and proposes directions for future research. Lastly, Section 8 encapsulates the conclusions drawn from this study.

2. Methodology

2.1. Random Forest

Random Forest (RF) is an ensemble learning method that utilizes decision tree-based models to establish connections between features and sample categories through nonlinear fitting [20]. The decision tree, a nonparametric algorithm, can be visualized as a tree-like classifier with nodes and directed edges. RF applies the Bagging ensemble method [21], extracting data from the training set through random resampling techniques to create multiple independent training subsets. These subsets are then used to train independent classifiers. In contrast to the traditional Bagging algorithm, RF introduces a randomized feature selection mechanism during the training of the base learner.

In Random Forest (RF), during the creation of each decision tree node, a subset of size k is randomly selected from the set of d features, and then an optimal feature is chosen from that subset for classification. Throughout the process of building decision trees in RF, both sample perturbation and feature perturbation are incorporated into the initial training set, effectively mitigating overfitting without the need for pruning. The classification process of RF is outlined in Figure 1, encompassing the following steps:

(1) Prior to constructing the decision tree, a unique training sample subset is generated by randomly sampling samples from the original dataset with replacement.

(2) The dataset comprises d features, with k (where k ≤ d) features chosen randomly to form a feature subset.

(3) The decision tree model is trained using the newly formed sample subset and feature subset. Optimal features for splitting are selected during the splitting process.

(4) Steps (1), (2), and (3) are iteratively executed to build n decision trees, resulting in the formation of a random forest. Each decision tree is constructed without employing pruning operations.

(5) The classification outcomes of the n decision trees are aggregated through a voting mechanism to determine the final category.

2.2. Moth-Flame Optimization

The moth-flame optimization (MFO) algorithm is designed based on how moths navigate laterally during nighttime flight [22]. During the night, moths maintain a fixed flight angle relative to the moon, as depicted in Figure 2a. Due to the moth’s distance from the moon, moonlight can be considered parallel light, ensuring that the moth flies in a straight path by maintaining this angle. When faced with obstacles, the moth can adjust its flight path based on this angle without deviating from its original direction. However, artificial light sources can mislead moths, causing them to spiral around the source, ultimately resulting in exhaustion and the phenomenon known as “Moths to the Flame,” as shown in Figure 2b.

Inspired by this behavior, the MFO algorithm interprets the moths’ spiral flight around the artificial light source as a search for an optimal solution in the solution space. In this algorithm, moths represent individuals in constant search, with the flame symbolizing the optimal position attained thus far. Each moth updates its position according to a specific equiangular helix flying around the corresponding flame. To increase the likelihood of attaining the optimal solution, the current best solution is used as the reference position for the subsequent generation’s flame. After each iteration, flame positions are rearranged based on fitness values, resulting in an updated flame sequence. Consequently, subsequent generations of moths update their positions in accordance with the corresponding sequence of flames.

The optimization process of the MFO algorithm can be expressed as a ternary optimality-seeking problem:

$$\left\{\begin{array}{l}MFO=\left(I,P,T\right)\\ I:\varphi \to \left\{M,OM\right\}\\ P:M\to M^{\prime }\\ T:M\to \left\{true,false\right\}\end{array}\right.$$

(1)

• I is a function that randomly generates moth populations and corresponding fitness values.
• ϕ represents the fitness function.
• M denotes the position of moths, with OM indicating the corresponding fitness value of moths in M.
• P is a mechanism for moths to update their position with a helix trajectory around the flame, denoted by M′ for the updated positions of moths.
• T is a judgment iteration function. If it meets termination conditions, T returns true and stops iteration; otherwise, it continues to iterate to find the optimal solution.

Through continuous iterations, moths actively eliminate poorly adapted flames while searching for superior ones, gradually reducing the flame number. This process balances the algorithm’s global search capability and local development within the search space. The formula used to adaptively decrease the flame number is as follows:

$$flame\_no=round(N-l\times {\displaystyle \frac{N-1}{T}})$$

(2)

In this context, the term “flame_no” represents the current count of flames, where “N” represents the maximum number of flames, “l” denotes the current iteration count, and “T” signifies the maximum number of iterations.

2.3. Gray Wolf Optimization

The gray wolf optimization (GWO) algorithm is inspired by the sophisticated hunting strategies observed in gray wolves [23]. Renowned for its simplicity, rapid convergence, minimal parameter setup, and ease of implementation, this algorithm offers significant advantages. Within the GWO algorithm, the gray wolf pack is structured into four hierarchies, denoted as α, β, δ, and ω from highest to lowest rank, as depicted in Figure 3. α, as the leader, holds the highest decision-making authority, with other wolves following α’s command. β assists α in decision-making or other activities, and the best β becomes α’s successor. δ follows commands from α and β and is responsible for hunting, scouting, and pack safety. ω ranks lowest in the hierarchy, relying on higher-ranking members for guidance to maintain population balance.

The optimization process of the GWO algorithm involves several steps, including social hierarchy division, prey encirclement, hunting, attacking, and prey searching. Mathematical descriptions of these phases are outlined below:

Social Hierarchy:

In the GWO algorithm, each gray wolf represents a potential solution, with the top three performers designated as α, β, and δ, representing the optimal, superior, and suboptimal solutions, respectively. The remaining wolves (ω) constitute candidate solutions within the population. α, β, and δ guide the rest of the wolves in the population for positional updates.

Prey Encirclement:

Gray wolves approach prey gradually and then encircle their target. The mathematical representation of these behaviors is elucidated through the following equation:

$$\left\{\begin{array}{l}D=|C·X_p(t)-X(t)|\\ X(t+1)=X_p(t)-A·D\\ A=2a·r_1-a\\ C=2·r_2\end{array}\right.$$

(3)

In these equations, the variable “t” represents the current iteration number. “A” and “C” denote coefficient vectors, while “X_p” signifies the position vector of the prey, and “X” represents the position vector of the gray wolf. The coefficient “a” gradually decreases from 2 to 0 during the iteration process. Moreover, “r₁” and “r₂” refer to random vectors within the range [0, 1].

Hunting:

In the confined search space, the exact location of the prey (optimal solution) remains undisclosed. Drawing inspiration from the hunting behavior of gray wolves, it is suggested that α, β, and δ possess superior knowledge about the prey’s position, thereby guiding other gray wolves to adjust their positions accordingly. This procedural representation is captured by the following equation:

$$\left\{\begin{array}{l}{D}_{\alpha }=\left|C_1·X_{\alpha }-X\right|,D_{\beta }=\left|C_2·X_{\beta }-X\right|,D_{\delta }=\left|C_3·X_{\delta }-X\right|\\ X_1=X_{\alpha }-A_1·(D_{\alpha }),X_2=X_{\beta }-A_2·(D_{\beta }),X_3=X_{\delta }-A_3·(D_{\delta })\\ X(t+1)={\displaystyle \frac{X_1+X_2+X_3}{3}}\end{array}\right.$$

(4)

where D_α, D_β, and D_δ represent the distances between α, β, and δ and ω; X₁, X₂, and X₃ represent the distances between the distance between α, β, and δ and the prey.

Engaging with Prey:

Gray wolves initiate their assault once the prey halts movement. To emulate the gradual approach of gray wolves toward their target, the parameter “a” steadily diminishes from 2 to 0, thereby concurrently adjusting the value of “A”. When the random value of “A” falls within the interval [−1, 1], the subsequent position of the gray wolf may lie anywhere between its current position and that of the prey. When ∣A∣ < 1, the wolves will commence attacking the currently sought-after prey.

Pursuit of Prey:

Gray wolves primarily rely on information from α, β, and δ to search for prey. Additionally, wolves disperse during their search for prey and then focus on hunting upon locating it. To simulate this dispersion, ∣A∣ > 1 is employed to compel the wolves to maintain distance from the current prey. Another search coefficient in the algorithm is C, representing the random weight of the prey. C undergoes random changes throughout the iteration process, aiding the algorithm in avoiding local optima during the final iteration. The settings of parameters A and C balance the global and local search capabilities of the algorithm.

2.4. Bayesian Optimization Algorithm

The Bayesian optimization (BO) algorithm emerges as a potent global optimization technique, widely applied for tuning hyperparameters in machine learning models due to its minimal iteration count and swift pace [24]. BO maximizes the utilization of already explored spaces and adjusts hyperparameters accordingly, diminishing unnecessary sampling and enhancing efficiency in solving intricate problems. At the core of the BO algorithm lies the probabilistic agent model and the acquisition function. The probabilistic agent model comprises two primary elements: the prior probability distribution and the observation model. The prior probability distribution furnishes initial insights to guide the estimation of model parameters, while the observation model utilizes existing data to refine the prior probability distribution, yielding an enriched posterior probability distribution with additional information. In this study, the Gaussian process agent model is employed to depict the relationship between hyperparameter x and the objective function f(x), as depicted in Equation (5). The acquisition function utilizes agent model insights to determine the next point for evaluation, striking a balance between exploration and exploitation to achieve global optimization.

$$\left\{\begin{array}{l}f(x)\sim GP(\mu (x),k(x,x^{\prime }))\\ \mu (x)=E[f(x)]\\ k(x,x^{\prime })=E[f(x)-\mu (x)(f(x^{\prime })-\mu (x^{\prime }))]\end{array}\right.$$

(5)

The expression “GP (*)” signifies the Gaussian process, where “µ(x)” represents the mean function, and “k (x, x’)” denotes the covariance function.

2.5. Hybrid Models

The aim of this investigation is to develop a classification predictive model for surrounding rock in TBM operations. To achieve this, three hybrid models—GWO-RF, MFO-RF, and BO-RF—are constructed by adjusting the hyperparameters of the RF using two metaheuristic algorithms (MFO and GWO) and BO. Despite differences in optimization approaches between the metaheuristic algorithms and the BO algorithm, they share similar steps in optimizing the RF models, as illustrated in Figure 4. The specific optimization process involves the following:

• Data Analysis and Preprocessing: The dataset is analyzed and normalized. Subsequently, it is randomly divided into training and test sets, ensuring that all three hybrid models utilize the same dataset.
• Parameter Initialization: The hyperparameter range for the RF model is set. Additionally, the appropriate population size and number of iterations are determined for the metaheuristic algorithms. For the BO algorithm, the same number of iterations as the metaheuristic algorithms is set.
• Fitness Evaluation: A fitness function is defined, and fitness evaluation is performed on the initial model.
• Parameter Update: The metaheuristic algorithms update the hyperparameter combinations based on the results of the previous iteration to achieve improved optimization results. In contrast, the BO algorithm utilizes a probabilistic model to guide parameter selection in the subsequent step and dynamically updates the model to explore the parameter space more effectively.
• Stopping Condition Check: When the maximum number of iterations or convergence is reached, the optimization process is terminated, and the optimal hyperparameter combination of the hybrid model is obtained.

Figure 4 illustrates the overall construction process of the hybrid models.

3. Data

3.1. Data Source and Description

The data utilized in this study were gathered from the PSRWT tunnel located in Malaysia. Stretching over a length of 44.6 km, this tunnel serves the purpose of transporting raw water from the Semantan River in Pahang to the South Klang Valley region of Selangor state, addressing water scarcity challenges resulting from the area’s rapid population growth. The predominant lithology along the tunnel route consists of granite and metamorphic sedimentary rocks. Excavation for the PSRWT tunnel predominantly employed Tunnel Boring Machine (TBM) methods, covering 34.74 km, with the remaining section excavated using conventional drill and blast techniques and excavation overlay methods. The primary TBM utilized for excavation was a main girder type TBM provided by Robbins USA, featuring a cutter diameter of 5.23 m.

For this study, 544 datasets were collected from the PSRWT tunnel, forming the basis for constructing hybrid prediction models based on Random Forest (RF): GWO-RF, MFO-RF, and BO-RF, with the aim of predicting the class of surrounding rock in real-time during TBM excavation. These datasets comprised 258 sets from the fresh zone of the rock mass and 286 sets from the slightly weathered zone. During excavation, the TBM gathers substantial tunneling data, with average data recorded every 10 m. The condition of the rock mass at the tunnel face influences the TBM’s operational and performance parameters through interaction, leading to the selection of parameters such as thrust force per cutter (TFC), revolution per minute (RPM), penetration rate (PR), advance rate (AR), penetration per revolution (PRev), and field penetration index (FPI) as input parameters for predictive models. Specifically, thrust force per cutter (TFC) represents the force exerted by each cutter on the rock, indicating the interaction between the TBM cutter and the rock. Revolutions per minute (RPM) refers to the rotational speed of the cutterhead and is a key indicator of TBM performance, with RPM typically decreasing in harder rock to reduce tool wear. Penetration rate (PR) measures the depth the cutterhead advances per minute, reflecting the efficiency of the TBM; a higher PR is generally observed in softer rocks. Advance rate (AR) is the distance advanced by the TBM per unit of time, directly representing overall excavation efficiency. Penetration per revolution (PRev) indicates the depth the cutterhead advances per revolution, which reflects the cutting efficiency and the relationship with rock hardness. Finally, the field penetration index (FPI) measures the efficiency of cutting under a given thrust, providing a comprehensive view of the tool performance, thrust, and rock properties. These parameters are monitored and calculated in real time by the TBM’s multiple sensors.

Specific conditions of the rock mass at the tunnel face were determined through field investigations and indoor tests. The rock mass rating (RMR) system proposed by Bieniawski was employed for classifying the surrounding rock at the tunnel face. The RMR method involves six parameters: intact rock strength, rock quality index, joint spacing, joint status, and groundwater status. It classifies the rock mass into five classes ranging from class I to class V, representing very good, good, fair, poor, and extremely poor conditions, respectively. Widely recognized in geotechnical engineering, the RMR method is valued for its simplicity and comprehensive assessment of rock mass quality, making it widely applied in practical engineering design and construction.

3.2. Data Pre-Process

The dataset utilized in the study consists of three classes of surrounding rock samples: 167 class I samples (30.7%), 320 class II samples (58.8%), and 57 class III samples (10.5%), as depicted in Figure 5. Notably, there exists an imbalance in sample distribution among the classes, particularly evident in the differing number of class II and class III samples, with a ratio reaching 5.6:1. Such class imbalance can introduce biases in machine learning models, where the more frequent classes tend to dominate during training, potentially resulting in inferior predictions for less represented classes. Figure 6 presents the data distribution and correlation analysis among the six input parameters. It reveals a strong positive correlation (correlation coefficient of 0.73) between TFC and FPI, while the correlations among the remaining parameters are relatively weaker. The histograms in Figure 6 show the distribution of data for each feature, with the x-axis of the histograms indicating the range of feature values and the y-axis indicating the corresponding frequencies, and these histograms visualize the concentration trend and dispersion of the features. Figure 7 illustrates box plots for the TBM tunneling data, showcasing statistical metrics such as the median, upper and lower quartiles, and outliers. Outliers are determined based on the statistical principles of the boxplot. The normal range is defined by the quartiles (Q1 and Q3) and the interquartile range (IQR). Data outside the range from Q1 − 1.5 × IQR to Q3 + 1.5 × IQR are considered outliers, typically indicated as outliers in the box plot. Some variables exhibit outliers, notably TFC, which were retained in this study despite their unidentified origins, considering their potential informational value.

Randomized stratified sampling was employed to partition the 544 datasets into an 80% training set and a 20% test set, maintaining the original database’s data structure. The training set was utilized for the model to understand the relationship between TBM tunneling parameters and the surrounding rock class, while the test set was used to evaluate the model’s performance. Various strategies were implemented during data preprocessing to tackle class imbalance. While under-sampling and over-sampling are common approaches, under-sampling may result in the loss of crucial information by reducing sample numbers across most categories. In contrast, basic SMOTE randomly generates synthetic samples for the minority class but does not consider the distribution of samples near the decision boundary, which may lead to overlapping among classes and reduce classification precision. To mitigate this, the Borderline Synthetic Minority Over-sampling Technique (Borderline-SMOTE), an enhanced SMOTE oversampling method, was employed. Borderline-SMOTE enhances the model’s adaptability to imbalanced data by strategically generating synthetic samples near the decision boundary of minority categories [25]. This approach better preserves the minority class’s critical information, reducing the risk of overgeneralization and ensuring a more accurate representation of minority categories. This technique balanced the sample differences across surrounding rock classes, ensuring adequate representation of class III surrounding rock samples to enhance classification performance during training. Before model training, Z-score normalization was applied to eliminate magnitude differences and data biases among input parameters, facilitating model learning and optimization.

4. Modeling

4.1. Model Metrics

Evaluating machine learning model performance comprehensively and accurately is integral to the model development process. Proper selection of evaluation metrics is crucial for understanding the effectiveness and performance of the model [26]. In this study, three global evaluation metrics—accuracy, Kappa index (Kappa), and Matthews correlation coefficient (MCC)—were utilized to assess the model’s overall prediction capability for surrounding rock classification [27]. Accuracy refers to the proportion of correctly predicted samples relative to the total number of samples, providing an overall measure of the model’s predictive capability. Kappa assesses the agreement between the model’s predictions and the actual outcomes, adjusting for random chance. The Matthews Correlation Coefficient (MCC) incorporates true positives, false positives, true negatives, and false negatives, providing a comprehensive evaluation of the model’s performance. Precision measures the proportion of true positive samples among those predicted as positive, reflecting the model’s accuracy in predicting positive instances. The recall represents the proportion of true positives correctly identified by the model, indicating its ability to capture positive samples. The F1-score is the harmonic mean of precision and recall, offering a balanced evaluation of the model’s performance across both dimensions. Higher values of these metrics, approaching 1, indicate better model predictions, with 1 signifying perfect agreement between predicted and true values. Given the presence of unbalanced data in the database, relying solely on global evaluation metrics may not sufficiently reveal model performance across all categories. Hence, local evaluation metrics such as precision, recall, and F1-score were employed to gain insights into the model’s performance across different classes [28]. These global and local evaluation metrics can be computed from the confusion matrix, as depicted in Figure 8.

4.2. Model Training

This study utilized six TBM digging parameters (TFC, RPM, PR, AR, Prev, and FPI) as input parameters for predicting three rock mass classifications. Three optimization algorithms—MFO, GWO, and BO—were employed to adjust the hyperparameters of RF. The optimization ranges for the hyperparameters of RF were set to (10, 200) for n_estimators, (2, 10) for max_depth, (1, 10) for min_samples_leaf, and (2, 20) for min_samples_split. In metaheuristic algorithms like MFO and GWO, population size plays a crucial role in optimization ability. Larger populations enhance exploration but increase computational cost and slow convergence, while smaller populations converge faster but risk local optima. Therefore, selecting a reasonable population size for metaheuristic optimization algorithms is essential for performance and effectiveness. Population sizes were set to 20, 30, 40, 50, 60, and 70, and iterative calculations were performed accordingly. The effectiveness of the BO algorithm relies on both the probabilistic agent model and the acquisition function. In this study, the Gaussian process was chosen as the probabilistic agent model, while GP-Hedge served as the acquisition function. Accuracy was utilized as the fitness function for the iterations. After 100 iterations, the fitness values of the three optimization algorithms stabilized, as shown in Figure 9. Analysis of the iteration curves of the MFO algorithms (Figure 9a) reveals that results from different population sizes converge to the same value. Although prediction accuracies remain consistent across various populations, differences primarily occur in convergence speed. Notably, the model achieves maximum accuracy most rapidly with a population size of 50. Hence, the model corresponding to a population size of 50 was selected as the optimal MFO-RF model. Examination of the GWO algorithm’s iteration curve (Figure 9b) indicates that both GWO-RF models with populations of 30 and 50 achieve maximum accuracy. However, the model with a population of 50 converges faster, leading it to be considered the optimal GWO-RF model.

Table 1. Development results of three hybrid models.

Models	Swarm	Accuracy	Iteration
MFO-RF	20/30/40/50/60/70	0.992	100
GWO-RF	20/40/60/70	0.985	100
GWO-RF	30/50	0.987	100
BO-RF		0.985	100

presents the accuracy results for the three hybrid models. A comprehensive evaluation of the three RF-based hybrid models (MFO-RF with 50 swarms, GWO-RF with 50 swarms, and BO-RF) will be provided below. “Swarms” refers to the population size in metaheuristic optimization algorithms, representing the number of candidate solutions used in each iteration of the algorithm.

5. Results and Discussion

5.1. Model Evaluation

To provide a more intuitive assessment of classification accuracy for each model across individual classes, confusion matrices were constructed for each model using the training set, as presented in Figure 10. Normalized confusion matrices were utilized in this study, offering relative proportions for each category by normalizing at each row. This enhances visualization and comprehension, particularly in situations with data category imbalances. The figure illustrates classification results for the three surrounding rock classes (I, II, and III). Each row represents the actual surrounding rock class, and each column represents the predicted surrounding rock class. Percentages along the diagonal signify correctly classified instances, with off-diagonal percentages indicating misclassifications. Darker colors indicate higher percentage values. Analysis of the confusion matrix plot reveals that the unoptimized RF model demonstrates superior discrimination between classes I and II but weaker discrimination for class III, with 22.18% of class III samples misclassified as class II. The proportion of correctly classified samples in each category is significantly improved by the optimization algorithm, indicating that the classification performance of RF can be enhanced using different optimization algorithms.

Three global evaluation indices (accuracy, Kappa, and MCC) and three local evaluation indices (precision, recall, and F1-score) were used to comprehensively quantify the prediction ability of the surrounding rock classification for each model, as presented in

Table 2. The assessment results of each model in the training stage.

Model	Accuracy	Kappa	MCC	Precision	Recall	F1-Score
RF	0.801	0.702	0.708	0.823	0.801	0.806
MFO-RF	0.992	0.988	0.988	0.992	0.992	0.992
GWO-RF	0.987	0.980	0.980	0.987	0.987	0.987
BO-RF	0.985	0.978	0.978	0.986	0.985	0.985

. It is evident from the table that all optimized hybrid models exhibit high prediction accuracy on the training set, with accuracies exceeding 0.98. Among them, the MFO-RF model demonstrates optimal performance on both global and local evaluation metrics (0.992 for accuracy, 0.998 for Kappa, 0.988 for MCC, 0.992 for precision, 0.992 for recall, and 0.992 for F1-score). This indicates that the MFO-RF model possesses the strongest prediction ability for surrounding rock classes and can accurately distinguish between them, showing a high level of consistency between predicted and actual results. The GWO-RF and BO-RF models also exhibit comparable performance and outperform the unoptimized RF model. To determine the optimal model, evaluation metrics for the four models were scored on a scale from 1 to 4, where higher scores reflect superior predictive ability. This scoring system enables explicit comparison of performance between models and provides an objective basis for selecting the best model. Figure 11 displays the final scores of each model in the training phase. The MFO-RF model obtained full scores across all indicators, ultimately achieving the highest score of 24, indicating optimal performance. The prediction performance of each model in the training stage is ranked as MFO-RF > GWO-RF > BO-RF > RF.

To thoroughly assess the generalization ability of the constructed models, validation was conducted on a designated test set. Test set validation is a crucial step to verify whether the model can achieve good performance on unknown data and determine the reliability of the model in practical applications. Figure 12 illustrates the normalized confusion matrix of the four models on the test set. It is apparent that the unoptimized RF model exhibits a higher classification error rate compared to the hybrid models, particularly underperforming in class III where the sample size is limited. Notably, all three hybrid models demonstrate superior performance in class I, with accuracy above 90%. For class II, the hybrid model MFO-RF performs the best with 94.44% accuracy. GWO-RF shows 100% correct classification on class III.

Table 3. The assessment results of each model in the testing stage.

Model	ACC	Kappa	MCC	Precision	Recall	F1
RF	0.782	0.595	0.605	0.841	0.782	0.804
MFO-RF	0.927	0.852	0.853	0.933	0.927	0.929
GWO-RF	0.873	0.759	0.766	0.897	0.873	0.878
BO-RF	0.891	0.786	0.789	0.905	0.891	0.895

presents the global and local evaluation metrics of the four models on the test set. It is evident from the table that the hybrid models continue to exhibit satisfactory performance on the unknown dataset without experiencing overfitting issues. Moreover, the MFO-RF model still performs very well on the test set (0.927 for accuracy, 0.852 for Kappa, 0.853 for MCC, 0.933 for precision, 0.927 for recall, and 0.929 for F1-score). Figure 13 illustrates the scores of each evaluation metric for the four models on the test set. The MFO-RF model still obtains the highest score. Unlike the ranking observed during the training phase, the BO-RF model demonstrates superior performance compared to the GWO-RF model on the test set, confirming the noteworthy optimization capability of the BO algorithm.

The evaluation of performance on both training and test sets highlights the superior learning and generalization capabilities of the MFO-RF model. To further investigate its advantages, various commonly utilized machine learning models were developed for predicting surrounding rock classes. These models include decision tree (DT), SVM, KNN, ANN, logistic regression (LR), AdaBoost, and Naive Bayes (NB). Trained on the same dataset, these models were subsequently assessed using the same test set. Figure 14 illustrates the accuracy of the MFO-RF model compared to other machine learning models on both training and test datasets.

It is evident that while the DT model performs well on the training dataset, it exhibits poor performance on the test dataset, indicating severe overfitting problems. The other models also show inferior performance compared to MFO-RF on both the training and test datasets. This indicates that the developed hybrid model, MFO-RF, can achieve more accurate surrounding rock classification compared to other commonly used machine learning models.

5.2. Model Interpretation

The poor interpretability of machine learning models often hampers their application in practical engineering fields. Following the identification of the optimal hybrid model for predicting surrounding rock classification (MFO-RF), it is essential to analyze its interpretability. This involves extracting the key factors influencing model performance and revealing the internal relationship between surrounding rock classification and each excavation parameter of the TBM.

The constructed MFO-RF model was interpreted using the Shapley Additive Explanations (SHAP) interpretation method. SHAP is a nonparametric technique for interpreting machine learning models, rooted in the Shapley value from cooperative game theory [29]. It quantifies the contributions of features to model performance, revealing both local and global model behavior. Figure 15 presents a summary plot of SHAP, illustrating the contribution of each feature to different surrounding rock classes in the MFO-RF model. On the left side of the figure, features are ranked by importance to the prediction of surrounding rock classification, with their significance indicated by their position from top to bottom. Each dot in the figure represents a sample of surrounding rock classes, with the horizontal axis denoting the SHAP value, indicating the feature’s influence on prediction outcomes. A positive SHAP value increases the probability of prediction results, while a negative value decreases it. Sample point color, ranging from blue to red, represents feature values from small to large. Relative importance rankings of each feature across the three surrounding rock classes are depicted in Figure 15a–c.

According to Figure 15a, PR, AR, and RPM are significant predictors of surrounding rock class I. Lower values of these features correlate with a higher probability of the surrounding rock being classified as class I, affirming that more intact surrounding rock corresponds to a slower TBM digging speed. In Figure 15b, RPM, FPI, and TFC exert greater influence on predicting surrounding rock class II. For class III, shown in Figure 15c, TFC, PR, and FPI are the key features in model prediction, with higher feature values associated with a higher probability of predicting class III. For the intact and unweathered rock mass, the TBM needs to reduce the rotational speed of the tool to ensure safe and effective cutting of the rock mass, and the working efficiency of the TBM and the overall efficiency of the tunnel boring will also be reduced correspondingly. Conversely, relatively more fragmented rock allows for higher dig ability and TBM efficiency. Thus, employing TBM tunneling parameters enables accurate prediction of the corresponding surrounding rock class, guiding safe and efficient on-site construction.

Based on the analysis in Figure 15, Figure 16 demonstrates the relative importance of overall variables for the three surrounding rock classes. Notably, the contribution values of TFC, RPM, PR, AR, PRev, and FPI are 0.156, 0.214, 0.226, 0.215, 0.079, and 0.11, respectively. “Contribution values” refer to the relative importance or contribution of each variable in predicting the target class, indicating the extent to which each variable influences the model’s output. The contributions of PR, AR, and RPM are particularly prominent, with PR making the most substantial contribution to predicting surrounding rock classes.

6. Validation

To further validate the accuracy and reliability of the hybrid models developed for engineering applications, an additional 91 sets of data were collected, all derived from the PSRWT tunnel. The new dataset consists of 25 class I samples, 60 class II samples, and six class III samples, exhibiting a similar class imbalance phenomenon as the original dataset. Figure 17 displays categorized report plots of the three hybrid models on the validation dataset, offering a comprehensive assessment of their performance across each class and overall. Figure 17a,c indicate that the MFO-RF and BO-RF models outperform others in predicting surrounding rock class II, followed by class I, with relatively poorer performance in class III. In Figure 17b, the GWO-RF model performs similarly to the MFO-RF and BO-RF models in ranking the performance of each surrounding rock class, but the performance is slightly worse. Overall, the hybrid models demonstrate satisfactory prediction performance, notably the consistent performance of MFO-RF and BO-RF models, both achieving accuracies of 0.879. These results confirm the engineering applicability of the developed models.

Furthermore, a graphical user interface (GUI) was created to facilitate rapid prediction of rock mass classification, as presented in Figure 18. The interface consists of input parameters, a prediction model selection, and output results. Users simply input prediction index values and select the desired prediction model to obtain corresponding surrounding rock classifications.

7. Limitations and Future Studies

Although the developed hybrid models have shown promising results in predicting surrounding rock classes, this study has several limitations:

(1) The applicability of the developed models is limited to similar tunnel projects. Future efforts could focus on creating more universally adaptable models by collecting TBM tunneling data from various tunnel projects.

(2) The dataset comprises only three rock mass classes, with a relatively small number of samples for class III. This limited sample size may constrain the prediction performance of the model for class III. Expanding the training dataset could enhance the model’s predictive capabilities across all classes.

(3) The accuracy of the developed models may be influenced by outliers in the database. In this study, outliers were not addressed, as their source is unknown and they may contain valuable information. Future research could employ more robust outlier handling techniques to effectively utilize potentially valuable data insights.

8. Conclusions

TBMs are extensively utilized in underground engineering for efficient and safe tunnel construction. However, a persistent challenge lies in the timely perception of the geological environment ahead of the tunnel face, posing risks and uncertainties to construction. Addressing this issue, this study developed three novel hybrid models using TBM boring parameters, along with optimization algorithms like MFO, GWO, and BO, combined with RF models, to enable the real-time prediction of rock mass classification ahead of the tunnel face. The main conclusions of this study are as follows:

(1) Comparative analysis revealed that all developed hybrid models outperformed the unoptimized RF model in prediction accuracy. Among them, the MFO-RF model exhibited the highest prediction performance, achieving 0.992 training accuracy and 0.927 testing accuracy. Furthermore, compared to other commonly employed machine learning models, the MFO-RF model still showed superior learning and generalization performance in the prediction of surrounding rock classification.

(2) SHAP was introduced to interpret the MFO-RF model. PR, AR, and RPM were identified as the key input parameters for surrounding rock classification prediction, and the effects of these key parameters on predicting different rock mass classes were analyzed.

(3) Further data gathered from the PSRWT tunnel were utilized to verify the precision and consistency of the hybrid models. Validation results indicated that the hybrid models generally attained satisfactory outcomes, with the accuracy of the MFO-RF and BO-RF models recorded at 0.879, while the GWO-RF model exhibited an accuracy of 0.857.