*Article* **Estimating Flyrock Distance Induced Due to Mine Blasting by Extreme Learning Machine Coupled with an Equilibrium Optimizer**

**Ramesh Murlidhar Bhatawdekar <sup>1</sup> , Radhikesh Kumar 2,\*, Mohanad Muayad Sabri Sabri <sup>3</sup> , Bishwajit Roy <sup>4</sup> , Edy Tonnizam Mohamad <sup>1</sup> , Deepak Kumar <sup>5</sup> and Sangki Kwon <sup>6</sup>**


**Abstract:** Blasting is essential for breaking hard rock in opencast mines and tunneling projects. It creates an adverse impact on flyrock. Thus, it is essential to forecast flyrock to minimize the environmental effects. The objective of this study is to forecast/estimate the amount of flyrock produced during blasting by applying three creative composite intelligent models: equilibrium optimizer-coupled extreme learning machine (EO-ELM), particle swarm optimization-based extreme learning machine (PSO-ELM), and particle swarm optimization-artificial neural network (PSO-ANN). To obtain a successful conclusion, we considered 114 blasting data parameters consisting of eight inputs (hole diameter, burden, stemming length, rock density, charge-per-meter, powder factor (PF), blastability index (BI), and weathering index), and one output parameter (flyrock distance). We then compared the results of different models using seven different performance indices. Every predictive model accomplished the results comparable with the measured values of flyrock. To show the effectiveness of the developed EO-ELM, the result from each model run 10-times is compared. The average result shows that the EO-ELM model in testing (R2 = 0.97, RMSE = 32.14, MAE = 19.78, MAPE = 20.37, NSE = 0.93, VAF = 93.97, A20 = 0.57) achieved a better performance as compared to the PSO-ANN model (R2 = 0.87, RMSE = 64.44, MAE = 36.02, MAPE = 29.96, NSE = 0.72, VAF = 74.72, A20 = 0.33) and PSO-ELM model (R2 = 0.88, RMSE = 48.55, MAE = 26.97, MAPE = 26.71, NSE = 0.84, VAF = 84.84, A20 = 0.51). Further, a non-parametric test is performed to assess the performance of these three models developed. It shows that the EO-ELM performed better in the prediction of flyrock compared to PSO-ELM and PSO-ANN. We did sensitivity analysis by introducing a new parameter, WI. Input parameters, PF and BI, showed the highest sensitivity with 0.98 each.

**Keywords:** flyrock; weathering index (WI); equilibrium optimizer (EO); particle swarm optimization (PSO); extreme learning machine (ELM); artificial neural network (ANN)

### **1. Introduction**

Surface mining is basically for breaking in situ rock, construction activities, and excavation; blasting is the most popular method throughout the globe [1–3]. Blasting is a system consisting of interaction between explosives and rock [4]. Blast design, properties of explosives, and rock mass are primarily the existing critical parameters for blasting and its performance [5,6]. Standard operating procedures are adopted for blast execution [7,8]. Desired blast fragmentation throw and shape muck pile affect the efficiency

**Citation:** Bhatawdekar, R.M.; Kumar, R.; Sabri Sabri, M.M.; Roy, B.; Mohamad, E.T.; Kumar, D.; Kwon, S. Estimating Flyrock Distance Induced Due to Mine Blasting by Extreme Learning Machine Coupled with an Equilibrium Optimizer. *Sustainability* **2023**, *15*, 3265. https://doi.org/ 10.3390/su15043265

Academic Editors: Steve W. Lyon and Guang-Liang Feng

Received: 4 September 2022 Revised: 20 September 2022 Accepted: 17 October 2022 Published: 10 February 2023

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

of loading equipment [9]. Blast fragmentation is also crucial for loading equipment and downstream operation of hauling and crushing [10,11]. Hence, these parameters are called favorable parameters [12,13], whereas flyrock, ground vibration, air overpressure, and dust affect the environment [7,14]. Therefore, these parameters are known as unfavorable parameters [8,15].

Flyrock is a geotechnical issue based on various rock mass properties. Various studies have been carried out to resolve geotechnical issues. Failure in geotechnical structure was studied with a multiscale work analysis approach [16]. An impressive micromechanical modeling (MM) framework was proposed by utilizing the discrete element method (DEM) and the micro-mechanical (MM) model [17]. The impact of flow direction vis-a-vis gravity direction on suffusion in geotechnical structures or slopes was resolved using the pioneering computational fluid dynamic-discrete element method (CFD-DEM) [18]. Further, the CFD-DEM model was helpful in investigating the particle shape effect and various levels of transmission at the macro- and micro-behavior levels during suffusion [19,20]. The CFD-DEM method played a significant role in the investigation of seepage in the underwater tunnel face [21]. The novel multi-scale approach, by deploying the smoothed particle hydrodynamics (SPH) method, was found to be efficient in computational analysis to understand granular collapse [22]. A possible solution for several engineering and industrial processes was found by developing DEM for irregular 3D particle shapes [23]. An algorithm was developed to accomplish 3D realistic stones of irregular geometries at random for specified samples with quantitative adjustable control [24]. In the convention of blasting safety criteria, defining frequency characteristics of blasting is of practical significance. A computational method associated with the wavelet frequency domain parameter or a main frequency band was proposed [25]. The innovative liquid carbon dioxide rock-breaking technology was found to be safer than explosive blasting. This technology needs further investigation to show that this technology is more efficient than traditional non-explosive techniques [26].

Flyrock has also been considered as an incidental or haphazard, exorbitant throw of rock pieces originating from blasting operations. Rock fragments from blasting may be thrown beyond the expected normal distance. This may result in a serious hazard to people working around the mines or severe danger to the property and machinery near the blasting site [27–29]. There are several causes of flyrock: usage of types of explosives, improper blast design, or explicit or uncertain conditions of the rock mass. Flyrock accidents are caused by poor security or blast management practices [30–32]. During the last decade, several researchers have developed many computational methods concurrently to forecast flyrock due to blasting [33].

Various researchers have pointed out that ANN, which is a branch of artificial intelligence (AI), is suitable for forecasting engineering problems [34,35]. For decades, many researchers had used ANN as a prediction method for flyrock distance [36–38]. However, ANN has some limitations, including slow learning speed and it falls into local minima due to the use of a gradient-based optimizer [39]. Furthermore, a specific ANN model for the prediction of flyrock is not available. On the other hand, the boundary condition of ANN models depends upon the variation in data sets [31]. Further, despite significant learning cycles on the same data set, there may be a marginal improvement in prediction performance. ANN models can easily find out the significance and sensitivity of input parameters [31]. As stated in the literature, various metaheuristic optimization algorithms (MOAs) could be deployed to forecast flyrock created by blasting due to the MOA's higher efficiency and to avoid the limitation of ANN. Some of the researchers compared ANN with other models, such as the ICA-coupled ANN model, which provided better performance as compared to ANN for the prediction of flyrock [40]. Furthermore, ANFIS showed superior performance as compared to ANN [41,42]. The PSO algorithm provided a powerful equation to predict flyrock due to blasting [43]. Gene expression programming (GEP) and the firefly algorithm (FA) were used to compare the results of flyrock prediction [44]. Further, hybrid algorithms were developed using optimization algorithms and ANN. MOAs provide a powerful ability to search for the best local solution from global optima. Therefore, MOA's estimated biases and weights of ANN can improve the prediction task of flyrock. The hybrid model, PSO-ANN, was developed for the prediction of flyrock [45]. The results of two different hybrid models, genetic algorithm-based ANN (GA-ANN) and recurrent fuzzy neural network (RFNN-GA), were compared with ANN to predict flyrock [46]. The results of three hybrid models, ICA-ANN, PSO-ANN, and GA-ANN, were compared to predict flyrock [47].

In recent years, one of the most exciting areas of study is machine learning (ML) [48]. In general, ML describes the ways of making predictions about and learning from data. As a subcategory of artificial intelligence, ML mainly aims at proposing and developing algorithms that can learn automatically from data.

Specifically, AI seeks to identify the objects existing in the neighboring areas and predict how the environment behaves in ways to make informed decisions. As a result, the ML techniques have a higher tendency to predict instead of estimate. For instance, it discusses the way researchers can make use of data obtained from an interferometry experiment to predict the interference pattern that would be seen under a variety of conditions. Furthermore, ML-based methods are mostly used to address high-dimensional problems of a higher complexity compared to the problems that may generally arise during a conventional statistics course [49]. The capability of generating and analyzing large data sets has dramatically increased during the last three decades. Such "big data" revolution has been prompted due to an exponential upsurge in computation capacities and memory, which is recognized as Moore's law. The ML models SVM [50,51] and ORLEM [52] were used by researchers for the prediction of flyrock due to blasting. In this research work, three hybrid models: PSO-ANN, PSO-ELM, and EO-ELM, are developed for the prediction of flyrock and their results are compared.

#### **2. Models for the Prediction of Flyrock**

Numerous researchers have proposed various approaches for flyrock prediction, which includes empirical, semi-empirical, and mathematical models.

#### *2.1. Empirical Models for the Prediction of Flyrock*

Various empirical models were developed by several researchers, mainly for blast production, which depend upon blast design and/or rock mass properties. Figure 1 shows a schematic diagram of flyrock, which may result in face burst, cratering, and rifting. Face burst results when geological discontinuities or planes of weakness exist. Cratering happens due to the escape of gases in the stemming zone due to back breaks or weak rock. Furthermore, it may be caused due to incorrect delay sequence (back rows firing first as compared to front rows). Rifting is due to stemming release with air pulse and associated with air blast. Inadequate stemming length and inappropriate stemming material are the causes of rifting [53].

Various empirical equations were developed by several researchers for flyrocks, which occur in numerous sizes and shapes. The prediction of maximum flyrock was recommended based on a factor of safety and hole diameter [54]. Further, the relationship between the ratio of stemming length to burden was established to maximum flyrock distance [55]. *Sustainability* **2022**, *14*, x FOR PEER REVIEW 4 of 27

**Figure 1.** Schematic diagram of a mechanism of flyrock [56]. which is applied to determine the blast-hole pressure and effects. The impact time applied **Figure 1.** Schematic diagram of a mechanism of flyrock [56].

are the key parameters of cratering. Flyrocks increase due to cratering, with the increase in charge per m or decrease in stemming length. Similarly, rifting depends upon a charge per m, burden, and drill hole angle. Flyrock increase, with the increase in charge per m or decrease in burden and rifting, is minimum in cases of vertical drill hole angle.

may be suitable for a particular site only and are not accurate [40,59].

*2.3. Semi-Empirical Trajectory Physics-Based Models for the Prediction of Flyrock*

*2.2. Mathematical Models for the Prediction of Flyrock*

specific and having limited data input.

Flyrock prediction was developed based on blast design parameters (burden, stemming length, linear charge concentration, and specific charge) and rock mass properties (unconfined compressive strength and RQD) [58]. The empirical equation was established based on blast design parameters (stemming length, hole depth, burden, and spacing) and rock mass property (rock mass rating) for the prediction of flyrock [39]. These equations do not consider rock mass properties of tropically weathered rock. Various researchers have reported that empirical equations for the prediction of flyrock

Various researchers developed mathematical models for the prediction of flyrock. To estimate flyrock range, Lundborg [60] adopted a semi-empirical method to analyze the relationship between rock velocity and charge diameter. In terms of crater blasting in granite blocks, these authors introduced the relationship between the beginning velocity of the flyrock fragment, its size, and throw [61]. Two expressions were derived by Chiapetta (1983) in the case of distance that may be traveled by flyrock (Chiapetta RF, 1983). Furthermore, a relation was established by Roth (1979) to find out the flyrock travel range. According to this approach, all of the measurements were done on the flyrock range and the most important variable was to estimate the flyrock velocity at the beginning. Roth applied Gurney's proposed equation to compute the velocity at the beginning of the fragments thrown around through an explosion [61]. The limitation of mathematical equations is that their prediction is not accurate due to their being site-

In semi-empirical trajectory physics-based models, the focus is on the beginning velocity of flyrock; therefore, they are most desired. One of the models developed by St. George and Gibson was modified by Little and Blair [62]. These models generally suffer from inexplicity in defining the velocity of detonation and the density of the explosive,

Researchers developed equations for flyrock prediction for the calculation of initial velocity based on the scaled burden method [57,58]. As per the equation, a charge per m and burden are the key parameters of face burst, as face burst increases with the increase

Researchers developed equations for flyrock prediction for the calculation of initial velocity based on the scaled burden method [57,58]. As per the equation, a charge per m and burden are the key parameters of face burst, as face burst increases with the increase in charge per m or decrease in burden. Furthermore, stemming length or charging per m are the key parameters of cratering. Flyrocks increase due to cratering, with the increase in charge per m or decrease in stemming length. Similarly, rifting depends upon a charge per m, burden, and drill hole angle. Flyrock increase, with the increase in charge per m or decrease in burden and rifting, is minimum in cases of vertical drill hole angle.

Flyrock prediction was developed based on blast design parameters (burden, stemming length, linear charge concentration, and specific charge) and rock mass properties (unconfined compressive strength and RQD) [58]. The empirical equation was established based on blast design parameters (stemming length, hole depth, burden, and spacing) and rock mass property (rock mass rating) for the prediction of flyrock [39]. These equations do not consider rock mass properties of tropically weathered rock. Various researchers have reported that empirical equations for the prediction of flyrock may be suitable for a particular site only and are not accurate [40,59].

### *2.2. Mathematical Models for the Prediction of Flyrock*

Various researchers developed mathematical models for the prediction of flyrock. To estimate flyrock range, Lundborg [60] adopted a semi-empirical method to analyze the relationship between rock velocity and charge diameter. In terms of crater blasting in granite blocks, these authors introduced the relationship between the beginning velocity of the flyrock fragment, its size, and throw [61]. Two expressions were derived by Chiapetta (1983) in the case of distance that may be traveled by flyrock (Chiapetta RF, 1983). Furthermore, a relation was established by Roth (1979) to find out the flyrock travel range. According to this approach, all of the measurements were done on the flyrock range and the most important variable was to estimate the flyrock velocity at the beginning. Roth applied Gurney's proposed equation to compute the velocity at the beginning of the fragments thrown around through an explosion [61]. The limitation of mathematical equations is that their prediction is not accurate due to their being site-specific and having limited data input.

#### *2.3. Semi-Empirical Trajectory Physics-Based Models for the Prediction of Flyrock*

In semi-empirical trajectory physics-based models, the focus is on the beginning velocity of flyrock; therefore, they are most desired. One of the models developed by St. George and Gibson was modified by Little and Blair [62]. These models generally suffer from inexplicity in defining the velocity of detonation and the density of the explosive, which is applied to determine the blast-hole pressure and effects. The impact time applied to these equations is determined based on experimental observation instead of real monitoring.

#### *2.4. Artificial Intelligence Techniques*

Blasting is one of the major operations that causes several adverse environmental effects, such as generation of fines, ground vibrations, fumes, air blast, dust, and flyrock [44]. So, it is necessary to control these adverse effects while performing blasting operations [63]. During the last decade, various researchers have applied artificial intelligence (AI) techniques for the prediction, minimization, or optimization of these environmental effects. In recent years, ML is one of the trending methods. ML may be defined as computer algorithms that can improve automatically based on the nature of the signal or feedback provided to the learning system. It is divided into three categories, namely, reinforcement, unsupervised, and supervised learning. In "supervised learning", the model is trained based on provided inputs and their desired outputs and the goal is to learn the pattern by mapping input to output. In "Unsupervised learning", the model is trained based on provided inputs only where labels were not given so that it may find structure in the input on its own. In the case of "Reinforcement learning", the model gets across with a

dynamic environment to perform a certain goal and after performing this goal it obtains feedback that is equivalent to rewards, which the model attempts to maximize. ANN is one of the most significant algorithms of ML. ANN is a parallelly distributed system that epitomizes the neural network of the human brain to create information processing models by composing different networks and connections [64]. It has the advantages of self-organizing, adaptive, and real-time learning features that enable it to overcome the defects of traditional logic-based AI in handling unstructured information and intuition [65]. During 1992–1997, Vladimir Vapnik with colleagues developed support-vector machine (SVM) models at AT&T Bell Laboratories. SVM is one of the models of supervised learning that have associated learning algorithms to analyze data for regression or classification.

Despite intended flexibility in singleton AI models, the previous studies indicate that these algorithms may fail to deliver expected results and may experience poor generalization capability. The reason for this is that they may become stuck in a locally optimal solution due to the use of stochastic selection or gradient-based learning algorithms of learning parameters [66,67]. To get the precise result and to perform the tasks adequately, a hybrid combination of soft computing techniques may be used as it uses data pre-processing techniques or metaheuristic optimization approaches to solve those problems [65,67]. Furthermore, the metaheuristic approach coupled with the ML model may enhance the ML model's performance as it may adequately reach the local solution from the global best solution [64,67].

The ANFIS model is a hybrid model developed with ANN and fuzzy interface system (FIS). Various researchers have used the ANFIS model for the prediction of flyrock due to blasting [68]. The results are highly promising, and comparative analysis suggests that the proposed modeling approach outperforms ANNs and other traditional time series models in terms of computational speed, forecast errors, efficiency, peak flow estimation, etc. It was observed that the ANFIS model preserves the potential of the ANN approach fully, and eases the model building process.

Zhou [69] utilized PSO and ANN techniques to minimize flyrock due to blasting [68]. PSO-ANN was found better as compared to the ANN model. ICA, GA, extreme learning machine (ELM), and biography-based optimization (BBO) have been applied for the prediction of various geotechnical issues by many researchers. Murlidhar [70] reviewed the ANN-GA [69], ANN-PSO [71], and ANN- ICA models. Each of these hybrid models provided better accuracy as compared to single models [42]. Murlidhar [72] have applied PSO-ELM, BBO-ELM, and ELM models to predict flyrock due to blasting [73].

ELM has insufficient generalization ability in dealing with the samples due to the random initialization of parameters, and PSO uses the individuals sharing information in the group to move the whole group to evolve from disorder to obtain the optimal solution [74]. PSO-ELM takes advantage of PSO to search for global optimal solutions and ELM to quickly deal with the nonlinear relationship [70]. In other words, PSO-ELM uses the PSO algorithm to optimize the input weight matrix and the hidden layer bias in ELM to obtain an optimal network [75–78]. Therefore, the PSO-ELM model provides better performance as compared to the singleton ELM model. Similarly, the BBO-ELM model provides better performance as compared to the singleton ELM model because the hybrid model provides the advantage of both BBO and ELM models [72]. Hence, in this paper, the authors decided to compare three hybrid models: PSO-ANN, PSO-ELM, and EO-ELM.

#### **3. Background of Model**

#### *3.1. Extreme Learning Machine (ELM)*

The ELM was introduced by Guang-Bin Huang (2004) [79] for feedforward neural networks. It consists of one hidden layer with multiple hidden nodes and their parameters do not require the tuning of input weights. In ELM, output weights of hidden nodes are generally learned only in one step that result in the learning of a linear model. ELM has a very high generalization capability and is considerably faster than the feedforward neural network (back-propagation algorithm (BP). This is because of the dependency among parameters of the different layers in a feedforward neural network. Due to these dependencies, one is obliged to adjust all of the parameters (weights and biases). So, several iterative learning steps are required to improve the learning performance of the feedforward network. For these reasons, ELM has been widely applied for classification, regression, sparse approximation compression, clustering, future learning, and many more.

#### Essence of ELM

The hidden layer's parameters do not require tuning of input weights. "Randomness" is one of the ways to implement ELM other than semi-randomness, which is considered in many traditional methods. In ELM, the hidden layer's mapping follows the rules of ridge regress theory [80], universal approximation, and neural network generalization theory [81]. Figure 2 shows the mapping of input space and feature space. ELM has the ability that may bridge the current gaps among linear systems, neural networks, matrix theories, SVM, random projection, Fourier series, and others.

**Figure 2.** Mapping of input space and feature space.

The following steps need to be performed if a model consists of a particular number of hidden notes, node output function, and certain training set:


ELM is one of the ML algorithms which is free from tuning and works on the abovementioned steps. It consists of hidden nodes with high importance and has a high-speed learning process.

### *3.2. Artificial Neural Network*

The idea of an ANN is derived from the biological neural architecture where a very large number of biological neurons are interconnected through the links. It is an information processing model that is similar to the structure of neural networks present in the human brain in both structure and functions. It mimics the biological neural architecture in two ways: (i) the learning process is used by the network to acquire the knowledge from its surroundings and (ii) the acquired knowledge is stored by the synaptic weights (interconnection strength). An ANN consists of a network of interconnected processing units, which is capable of 'learning' to 'recognize' a complex input pattern, and predict the output pattern thereof. For this, the neural network is first 'trained' to analyze the input patterns and recognize the output that results from these inputs. The network is then able to recognize similarities in new input patterns and can predict the output. This property of a neural network makes it very useful for noisy (inexact) data to be interpolated and outputs predicted in terms of patterns that are already 'known' to it. This makes neural networks a ready replacement for older statistical techniques, such as linear regressing, multi-variable regression, and autocorrelation, etc. Even outputs that were previously not apparent to non-experts become recognizable, making the neural network a virtual expert. To solve different problems an ANN can be designed using the following three fundamental components:


### *3.3. Equilibrium Optimization*

The Equilibrium Optimizer (EO) was developed based on a generic mass balance equation [82]. The EO algorithm is designed with a high level of search capability from exploratory and exploitative systems to randomly change solutions and prevent local minima. In the Equilibrium Optimization (EO) algorithm, the equilibrium state (optimal result) is finally achieved through equilibrium candidates or best so-far solutions through search agents randomly [83]. If EO is compared with many metaheuristic algorithms, EO starts the optimization process based on the initial population. The construction of the initial concentrations is based on several dimensions and particles in the search space with consistent random initialization. Particles with their concentrations are known as search agents, corresponding to particles and positions in the PSO algorithm. During the process of optimization in the initial period, there is no awareness and understanding about the equilibrium state. In the beginning, equilibrium candidates are identified by updating concentrations randomly to fit solutions. During the whole optimization process, based on different experimentation and case studies four best-so-far particles are identified and another particle is the concentration of arithmetical mean of the aforesaid four particles. Exploration capability is based on best-so-far particles while exploitation capability is based on the average value. In several engineering problems, the generation rate can be expressed as a function of time [84]. In the case of EO, selecting an appropriate generation rate as well as updating concentrations randomly enhances EO's exploratory performance during the initial irritation and exploitation search in the final iterations. Thus, EO supports from beginning to end the complete optimization process and avoids local minima. Exploration and exploitation processes are balanced to obtain adaptive values for monitoring parameters resulting in a significant reduction in the movement of particles. To study the efficiency and effectiveness of EO, quantitative and qualitative metrics were endorsed. The EO algorithm showed higher efficiency (i.e., shorter computational time or limited iterations) in achieving optimal or close to optimal solutions with specific or most of the problems examined. The EO algorithm is undoubtedly a better algorithm as compared to the metaheuristic algorithms, such as GA and PSO, or the recently developed algorithms, such as GWO, GA, and GSA. The performance of EO is statistically comparable with SHADE and LSHADE-SPACM algorithms.

#### *3.4. Particle Swarm Optimization (PSO)*

The particle swarm optimization algorithm was first introduced by Kennedy and Eberhart (1995) [85]. To start operation, this algorithm distributes a set of entities (particles, each of which stands for a feasible solution) randomly in the search space. An objective function is considered as the determining factor of the swarm's goal. The fitness of every entity/particle is determined by the value that is correspondent to the objective function. Figure 3 shows a standard flow chart of a PSO.

algorithms.

*3.4. Particle Swarm Optimization (PSO)*

function. Figure 3 shows a standard flow chart of a PSO.

**Figure 3.** Standard flow chart of PSO [85]. **Figure 3.** Standard flow chart of PSO [85].

The PSO algorithm uses a random or stochastically selected population. The first step in the process is to select a population of particles or solutions and iterate the solutions until an optimum is reached. Each particle is assigned a value and then updated as per the 'gbest' fitness for two outstanding values. The solutions higher than 'gbest' fitness values are chosen and the 'gbest' fitness factors are updated. Subsequently, wherever the fitness of the particle is better or higher than that of 'pbest', the corresponding parameters The PSO algorithm uses a random or stochastically selected population. The first step in the process is to select a population of particles or solutions and iterate the solutions until an optimum is reached. Each particle is assigned a value and then updated as per the 'gbest' fitness for two outstanding values. The solutions higher than 'gbest' fitness values are chosen and the 'gbest' fitness factors are updated. Subsequently, wherever the fitness of the particle is better or higher than that of 'pbest', the corresponding parameters of 'pbest' are updated. The process would then enter the second phase, where the particles would be examined again.

computational time or limited iterations) in achieving optimal or close to optimal solutions with specific or most of the problems examined. The EO algorithm is undoubtedly a better algorithm as compared to the metaheuristic algorithms, such as GA and PSO, or the recently developed algorithms, such as GWO, GA, and GSA. The performance of EO is statistically comparable with SHADE and LSHADE-SPACM

The particle swarm optimization algorithm was first introduced by Kennedy and Eberhart (1995) [85]. To start operation, this algorithm distributes a set of entities (particles, each of which stands for a feasible solution) randomly in the search space. An objective function is considered as the determining factor of the swarm's goal. The fitness

#### of 'pbest' are updated. The process would then enter the second phase, where the particles would be examined again. *3.5. Case Study and Data Collection*

*3.5. Case Study and Data Collection* Basalt, granite, and limestone are common rocks for manufacturing aggregates in the construction industry of Thailand. Figure 4 shows potential rock resources of aggregates in different locations of Thailand. Small aggregate quarries produce less than 15,000 cubic meters per month. On the other hand, large aggregate quarries produce up to 150,000 cubic meters per month [86]. Vast quarries of limestone are used to supply limestone to manufacture Portland cement in factories situated 100 km from Thailand. The selected limestone quarry is an aggregate limestone quarry. Figure 5 shows photographs of an aggregate limestone quarry. Figure 6 shows the blasted muck pile of the face. During a blasting operation, flyrock is generated, which is concerning. In the selected quarry, the input parameters consisting of hole diameter, burden, stemming length, rock density, explosives charge per meter, powder factor, blast ability index, weathering index, and flyrock distance for 114 blasting events were collected. Figure 4 shows the location of the limestone quarry. Table 1 shows the details of all input parameters. The weathering index is site-specific.

*Sustainability* **2022**, *14*, x FOR PEER REVIEW 9 of 27

Based on weathering index, the blast input parameters are decided. Hence, sensitivity analysis of parameters is compared with weathering index. is site-specific. Based on weathering index, the blast input parameters are decided. Hence, sensitivity analysis of parameters is compared with weathering index.

Basalt, granite, and limestone are common rocks for manufacturing aggregates in the construction industry of Thailand. Figure 4 shows potential rock resources of aggregates in different locations of Thailand. Small aggregate quarries produce less than 15,000 cubic meters per month. On the other hand, large aggregate quarries produce up to 150,000 cubic meters per month [86]. Vast quarries of limestone are used to supply limestone to manufacture Portland cement in factories situated 100 kilometers from Thailand. The selected limestone quarry is an aggregate limestone quarry. Figure 5 shows photographs of an aggregate limestone quarry. Figure 6 shows the blasted muck pile of the face. During a blasting operation, flyrock is generated, which is concerning. In the selected quarry, the input parameters consisting of hole diameter, burden, stemming length, rock density, explosives charge per meter, powder factor, blast ability index, weathering index, and flyrock distance for 114 blasting events were collected. Figure 4 shows the location of the limestone quarry. Table 1 shows the details of all input parameters. The weathering index

**Figure 4.** Location map of Thailand quarry [86]. **Figure 4.** Location map of Thailand quarry [86].

**Figure 5.** Limestone quarry face at Thailand. **Figure 5.** Limestone quarry face at Thailand.

**Figure 6.** Blasted stone muck pile at quarry face.

**Rock Den-**

Symbol D B ST ρ CPM PF BI WI FR Unit mm m m Cum.t kg/m kg/cum % Ratio m Minimum 76 2.5 1.2 1.8 4.54 0.08 18.5 0.13 27 Quartile1 76 2.7 2 1.8 4.99 0.19 28.6 0.25 37 Average 90 3 2 2 7 0.30 43 0.76 81 Quartile3 102 3.6 2.95 2.5 8.99 0.40 54.6 0.88 82 Maximum 102 4.6 4 2.5 9.4 0.50 80.8 0.99 436

**sity**

**Table 1.** Input and output parameters.

**Burden Stemming** 

**Length**

**Parameters Hole Diam-**

**eter**

**Charge per M**

**Powder Factor**

Weathering index (WI) is a new parameter introduced based on rock mass properties,

such as water absorption (%), porosity (%), and point load index. Maximum values of water absorption and porosity are obtained for completely weathered granite. The maximum value of the point load index is obtained for fresh rock. At each blasting site,

**Blastability Index**

**Weathering Index** **Flyrock**

**Figure 6.** Blasted stone muck pile at quarry face. **Figure 6.** Blasted stone muck pile at quarry face.


**Figure 5.** Limestone quarry face at Thailand.


Symbol D B ST ρ CPM PF BI WI FR Unit mm m m Cum.t kg/m kg/cum % Ratio m Minimum 76 2.5 1.2 1.8 4.54 0.08 18.5 0.13 27 Quartile1 76 2.7 2 1.8 4.99 0.19 28.6 0.25 37 Weathering index (WI) is a new parameter introduced based on rock mass properties, such as water absorption (%), porosity (%), and point load index. Maximum values of water absorption and porosity are obtained for completely weathered granite. The maximum value of the point load index is obtained for fresh rock. At each blasting site, samples are collected, and each rock mass property is compared with the maximum value. The average of these ratios is known as the WI.

#### Average 90 3 2 2 7 0.30 43 0.76 81 **4. Model Development**

#### Quartile3 102 3.6 2.95 2.5 8.99 0.40 54.6 0.88 82 *4.1. Hybridization of PSO-ANN*

Maximum 102 4.6 4 2.5 9.4 0.50 80.8 0.99 436 Weathering index (WI) is a new parameter introduced based on rock mass properties, The PSO algorithm was developed as a bird swarm simulation by Kennedy and Eberhart. Swarm intelligence is the capability of an individual bird to deal with the previous experiences of the whole swarm. In PSO, the decision-making process is essential, and it can be made based on the following:


maximum value of the point load index is obtained for fresh rock. At each blasting site, Several researchers tried to enhance the generalization capabilities and performance of ANNs by using PSO algorithms because PSO is one of the robust global search methods that can enhance the performance capacity of ANN by adjusting its bias and weight of it. Furthermore, in the case of the local minimum, ANNs can increase the probability of convergence and, at that time, the PSO is likely to obtain the global minimum. Consequently, the developed PSO-ANN model acquires the search properties of both the ANN and PSO

models. So, in the case of the PSO-ANN model, the PSO searches for the global minimum and then exploits towards the local solution, which can be employed by the ANN to find the best results in the search space [87].

In a learning procedure for PSO-ANN, it starts with random assignment of weights and biases of a group of random particles. After this step the PSO-ANN model is trained based on the assigned weights and biases and then, at each iteration, the error is calculated between the predicted and actual value. After that, the calculated error is reduced by changing the particle position. By changing the particle position the best solution is selected and accordingly a new error is achieved. This complete learning process continues unless or until the termination criteria are fulfilled.

### *4.2. Hybridization of PSO-ELM*

Based on theory, Huang et al. [88] demonstrated the performance capability of ELM as a universal approximation, and the ability to engage many activation functions. Several researchers utilize ELM due to well-known features, such as fast learning capability and adequate ability to generalize, the same is deployed for prediction methods [79,89]. The generalizing ability of ELM is further enhanced by merging it with some other methods [90,91]. Several researchers from recent decades have successfully combined a naturebased/inspired algorithm to optimize the ELM model. Mohapatra et al. [92] designed a hybrid model consisting of the cuckoo search algorithm and ELM to classify medical data. The stability analysis of the photovoltaic interactive microgrid was carried out with a firefly algorithm by Satapathy et al. [93]. The evaluation of the aging degree of the insulated gate bipolar transistor was done with a whale optimization algorithm and ELM by Li et al. [94]. Figueiredo and Ludermir [95] studied the different topology of PSO—Global, Local, Von Neumann, Wheel, and Four Clusters—and showed, depending upon the problem, suitable topology, which need to be selected for the best PSO-ELM performance. Many researchers have used PSO-ELM for prediction in various engineering problems. PSO-ELM forecasting models were used to predict the regional groundwater depth [96]. The PSO-ELM approach was used for predicting landslide displacement interval [97]. The PSO-ELM model was deployed for predicting stabilized aggregate bases [98]. The PSO-ELM model was deployed to predict the vibration of the ground caused by the process of blasting [99]. Thus, from various research studies, an optimized version of ELM with other algorithms outperformed individual ELM accuracy levels in prediction jobs. ELM models generally get trapped in local minima because the initialization process is stochastic for the network input weights and hidden biases [100]. Various researchers have applied a combination of PSO and ELM to various areas reliably [101]. During the current study, to the best of our information/knowledge, the PSO-ELM model is developed for the first time to predict the flyrock caused by blasting. The flow chart of PSO-ELM is shown in Figure 7.

#### *4.3. Hybridization of EO-ELM*

This study proposes a new combination of hybrid ML models, called EO-ELM, where the EO optimizes the ELM learning parameters to find an optimal configuration of ELM for the prediction of flyrock. Here, the concentrations of EO are ELM learning parameters. The RMSE is considered an objective function for EO. The best equilibrium candidate found by EO is considered as the optimal configuration of ELM for prediction tasks.

chart of PSO-ELM is shown in Figure 7.

[101]. During the current study, to the best of our information/knowledge, the PSO-ELM model is developed for the first time to predict the flyrock caused by blasting. The flow

**Figure 7.** Flow chart of PSO–ELM. **Figure 7.** Flow chart of PSO–ELM.

*4.3. Hybridization of EO-ELM* This study proposes a new combination of hybrid ML models, called EO-ELM, where the EO optimizes the ELM learning parameters to find an optimal configuration of ELM for the prediction of flyrock. Here, the concentrations of EO are ELM learning parameters. The RMSE is considered an objective function for EO. The best equilibrium candidate found by EO is considered as the optimal configuration of ELM for prediction tasks. In EO-ELM, initially, all particles do not know the solution space. The collaboration In EO-ELM, initially, all particles do not know the solution space. The collaboration of five equilibrium candidates helps the concentration updating process of particles. At initial periods of iteration, the equilibrium candidates are diverse and the exponential term (F) produces large random numbers, which help particles to cover the entire solution space. Similarly, during the end period of iterations, particles are surrounded by equilibrium candidates which are in the optimum position with similar configurations. At these moments, the exponential term produces a lower value of random numbers, which helps fine-tune candidate solutions. The algorithm of EO-ELM is shown in the following Algorithm 1.

of five equilibrium candidates helps the concentration updating process of particles. At initial periods of iteration, the equilibrium candidates are diverse and the exponential term (F) produces large random numbers, which help particles to cover the entire solution space. Similarly, during the end period of iterations, particles are surrounded by

**Algorithm 1: The algorithm of EO-ELM. Exponential term (***F***),** *λ* **is a turnover rate and defined as a random vector in between 0 and 1,** *a***<sup>2</sup> is used to control the exploitation task.** *a***<sup>1</sup> is used to control the exploration task,** *sign*( <sup>→</sup>*<sup>r</sup>* <sup>−</sup>**0.5**) **component consequences the direction of intensification and diversification of particles,** *r* **is defined as a random vector in between 0 and 1, generation rate (G),** *r***1, and** *r***<sup>2</sup> denote the random values between 0 and 1.** *GCP* **is called generation rate control parameter** 1. Select training and testing dataset 2. Begin ELM training 3. Set hidden units of ELM 4. Obtain the number of input weights and hidden biases 5. Initialize the populations (P) 6. Initialize the fitness of four equilibrium candidates 7. Assignment of EO parameters value (*a*<sup>1</sup> = 2, *a*<sup>2</sup> = 1, *GP* = 0.5) 8. for it = 1 to maximum iteration number do 9. for i = 1 to P do 10. Estimate the fitness of the ith particle 11. if fitness ( → *Pi* ) < fitness ( → *Peq*[1] ) 12. Replace fitness ( → *Peq*[1] ) with fitness ( → *Pi* ) and → *Peq*[1] with → *Pi* 13. elseif fitness ( → *Pi* ) < fitness ( → *Peq*[1] ) & fitness ( → *Pi* ) < fitness ( → *Peq*[2] ) 14. Replace fitness ( → *Peq*[2] ) with fitness ( → *Pi* ) and → *Peq*[2] with → *Pi* 15. elseif fitness ( → *Pi* ) < fitness ( → *Peq*[1] ) & fitness ( → *Pi* ) < fitness ( → *Peq*[2] ) & fitness ( → *Pi* ) < fitness ( → *Peq*[3] ) 16. Replace fitness ( → *Peq*[3] ) with fitness ( → *Pi* ) and → *Peq*[3] with → *Pi* 17. elseif fitness ( → *Pi* ) < fitness ( → *Peq*[1] ) & fitness ( → *Pi* ) < fitness ( → *Peq*[2] ) & fitness ( → *Pi* ) < fitness ( → *Peq*[3] ) & fitness ( → *Pi* ) < fitness ( → *Peq*[4] ) 18. Replace fitness ( → *Peq*[4] ) with fitness ( → *Pi* ) and → *Peq*[4] with → *Pi* 19. end if 20. end for 21. → *Pmean* = ( → *Peq*[1] + → *Peq*[2] + → *Peq*[3] + → *Peq*[4] )/4 22. → *Peq*,*pool* = { → *Peq*[1] + → *Peq*[2] + → *Peq*[3] + → *Peq*[4] + → *Peq*[*mean*] } (Equilibrium pool) 23. Allocate *t* = <sup>1</sup> <sup>−</sup> *Iteration Maxiteration* (*a*2<sup>×</sup> *Iteration Maxiteration* ) 24. for i = 1 to P do 25. Random generation of vectors → *λ* and <sup>→</sup> *r* 26. Random selection of equilibrium candidate from equilibrium pool 27. Evaluate → *F* = *a*<sup>1</sup> × *sign*( → *r* − 0.5)[*e* − → *<sup>λ</sup><sup>t</sup>* <sup>−</sup> <sup>1</sup>] 28. Evaluate → *GCP* = { 0.5 ∗ *r*<sup>1</sup> *r*<sup>2</sup> ≥ *GP* 0*r*<sup>2</sup> < *GP* 29. Evaluate → *G*<sup>0</sup> = → *GCP* ∗ ( → *Peq* − → *λ* × → *P*) 30. Evaluate → *G* = → *G*<sup>0</sup> × → *F* 31. → *P* = → *Peq* + ( → *P* − → *Peq*)· → *F* + → *G* → *λV* ∗ (1 − → *F* ) (Concentration update) 32. end for 33. end for 34. Set ELM optimal input weights and hidden biases using → *Peq*[1] 35. Obtain output weights 36. ELM testing

### *4.4. Model Verification and Evaluation*

One of the most important aspects of the model development process is model verification and evaluation, as it is necessary to understand the behavior of the model to check the evolution of the model towards acquiring the accurate result and to know whether the

quality of the test model is excellent or not. To fulfil the desired, a training set is used to train the developed models and a different testing set is used to verify the model development. To evaluate the reliability of the developed model, seven different evaluation matrices, namely determination coefficient (R|2), root mean square error (*RMSE*), variance account factor (*VAF*), mean absolute error (*MAE*), Nash–Sutcliffe efficiency (*NSE*), mean absolute percentage error (*MAPE*), and a-20 index (A20), were used to define the relation between the actual and predicted value. Out of these evaluation matrices, the *RMSE* shows the standard deviation of the error between the actual and predicted values. The *MAPE* shows the error value percentage with the original data; having 0% *MAPE* shows the perfect model. The *NSE* value is a normalized statistic and is used to measure the goodness of fit of the model. Similarly, in the case of *MAE*, the goodness of the model increases with a decrease in the value of *MAE*. Further, R2 indicates the correlation between the actual and predicted values. The closer the value is to 1, the more perfect the model (1). Similarly, if the value of A20 is closer to 1, this shows a perfect prediction model (2). VAF shows the ratio of error variance to the measured data variance. The calculation formulas for different evaluation matrices are as follows:

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (y - \overline{y})^2} \tag{1}$$

$$R^2 = \left(\frac{\sum\_{i=1}^q \left(Y\_{E\_i} - Y\_{\overline{E}\_i}\right) \left(Y\_{O\_i} - Y\_{\overline{O}\_i}\right)}{\sqrt{\sum\_{i=1}^q \left(Y\_{E\_i} - Y\_{\overline{E}\_i}\right)^2 \sum\_{i=1}^n \left(Y\_{O\_i} - Y\_{\overline{O}^i}\right)^2}}\right)^2\tag{2}$$

$$MAE = \frac{1}{n} \sum\_{i=1}^{n} |(\hat{y}\_i - y\_i)| \tag{3}$$

$$MAPE = \frac{1}{n} \sum\_{i=1}^{n} \left| \frac{y\_i - \mathcal{Y}\_i}{y\_i} \right| \times 100\tag{4}$$

$$NSE = \left(1 - \frac{\sum\_{i=1}^{n} \left(R\_{O\_i} - R\_{E\_i}\right)^2}{\sum\_{i=1}^{n} \left(R\_{O\_i} - \overline{R\_{O\_i}}\right)^2}\right) \tag{5}$$

$$\text{VAF } (\%) = (1 - \frac{var(Y\_{E\_i} - Y\_{O\_i})}{var(Y\_{E\_i})}) \times 100\tag{6}$$

$$A20 = \frac{m20}{M} \tag{7}$$

#### **5. Results and Discussion**

The objective of this study was to predict the flyrock distance. Therefore, crucial blast design parameters were selected as input parameters (hole diameter, burden, rock density, stemming length, charge per meter, powder factor, blastability index, and weathering index). After that, the PSO-ELM, PSO-ANN, and EO-ELM models were developed and used for the prediction of flyrock. When data are split into train and test sets, it always becomes a challenging task to develop a generalized data-driven model. This work used 80% and 20% data split ratios for train and test data, respectively. These split data are used to test and compare the performance of the developed models.

To prove the model's effectiveness with the split dataset, a 10-times average run of three models is checked and compared. An average run of optimization-coupled ML models is useful to check the randomness in problem-solving for the optimal parameter set of optimization algorithms. Table 2 shows the optimal parameter set of metaheuristic algorithms, which are initially set by heat and trial method.


**Table 2.** Optimal parameter values of metaheuristic algorithms for eight hidden neurons in ELM.

During training, 500 iterations are set for each optimization-coupled ML model. The convergence plot of RMSE vs. iteration count is plotted and shown in Figure 8. Figure 8 shows faster (around 200 iterations) and better convergence ability of the EO-ELM compared to the PSO-ELM and PSO-ANN. The PSO-ELM and PSO-ANN become stuck in a local solution with premature convergence (Figure 8). Furthermore, the PSO-ANN may have an improper learning rate and overfitting issues due to the use of stochastic selection or gradient-based learning algorithms. As gradient-based learning algorithms intend to reach the minimum training error, but do not consider the magnitude of weights and only uses differentiable activation functions; due to this, they have less generalization performance [102]. In spite of this, other models use ELM, which is extremely fast and can train SLFNs using non-differentiable activation function to reach the solutions in a straightforward way without having issues, such as local minimum, improper learning rate, and overfitting. It tends to reach the smallest training error including the smallest norm of weights; due to this, it has the better generalization performance [102]. *Sustainability* **2022**, *14*, x FOR PEER REVIEW 16 of 27

The training phase scatter diagrams for the EO-ELM, PSO-ANN, and PSO-ELM are shown in Figure 9, Figure 10, and Figure 11, respectively. It is apparent from Figures 9–11 that the EO-ELM (Figure 9) predicts flyrock values more accurately compared to the PSO-

EO-ELM, PSO-ANN, and PSO-ELM are shown in Figure 12, Figure 13, and Figure 14, respectively. It is evident that the EO-ELM (Figure 12) predicts the test flyrock data better compared to the PSO-ANN (Figure 13) and PSO-ELM (Figure 14). Table 3 shows linear equations of the predicted and measured values for the EO-ELM, PSO-ANN, and PSO-

**Table 3.** Linear equations for predicted and measured values for the EO-ELM, PSO-ANN, and PSO-

**Model Training Data Testing Data** EO-ELM 67.10x + 73.93 96.66x + 100.59 PSO-ANN 58.06x + 74.58 88.9x + 92.43 PSO-ELM 64.58x + 73.92 99.10x + 98.85

**Figure 9.** Comparison of predicted flyrock with the EO-ELM vs. measured flyrock.

**Figure 8.** Convergence plot of three models (training period). **Figure 8.** Convergence plot of three models (training period).

ELM for training and testing data, respectively.

ELM for training and testing data.

The training phase scatter diagrams for the EO-ELM, PSO-ANN, and PSO-ELM are shown in Figures 9–11, respectively. It is apparent from Figures 9–11 that the EO-ELM (Figure 9) predicts flyrock values more accurately compared to the PSO-ANN (Figure 10) and PSO-ELM (Figure 11). The testing phase scatter diagrams for the EO-ELM, PSO-ANN, and PSO-ELM are shown in Figures 12–14, respectively. It is evident that the EO-ELM (Figure 12) predicts the test flyrock data better compared to the PSO-ANN (Figure 13) and PSO-ELM (Figure 14). Table 3 shows linear equations of the predicted and measured values for the EO-ELM, PSO-ANN, and PSO-ELM for training and testing data, respectively. ELM for training and testing data, respectively. **Table 3.** Linear equations for predicted and measured values for the EO-ELM, PSO-ANN, and PSO-ELM for training and testing data. **Model Training Data Testing Data** EO-ELM 67.10x + 73.93 96.66x + 100.59 PSO-ANN 58.06x + 74.58 88.9x + 92.43 PSO-ELM 64.58x + 73.92 99.10x + 98.85

The training phase scatter diagrams for the EO-ELM, PSO-ANN, and PSO-ELM are shown in Figure 9, Figure 10, and Figure 11, respectively. It is apparent from Figures 9–11 that the EO-ELM (Figure 9) predicts flyrock values more accurately compared to the PSO-ANN (Figure 10) and PSO-ELM (Figure 11). The testing phase scatter diagrams for the EO-ELM, PSO-ANN, and PSO-ELM are shown in Figure 12, Figure 13, and Figure 14, respectively. It is evident that the EO-ELM (Figure 12) predicts the test flyrock data better compared to the PSO-ANN (Figure 13) and PSO-ELM (Figure 14). Table 3 shows linear equations of the predicted and measured values for the EO-ELM, PSO-ANN, and PSO-

*Sustainability* **2022**, *14*, x FOR PEER REVIEW 16 of 27

**Figure 8.** Convergence plot of three models (training period).

**Figure 9. Figure 9.**  Comparison of predicted flyrock with the EO-ELM vs. measured flyrock. Comparison of predicted flyrock with the EO-ELM vs. measured flyrock.

**Figure 10.** Measured flyrock and predicted values of flyrock with the PSO-ANN for training data. **Figure 10.** Measured flyrock and predicted values of flyrock with the PSO-ANN for training data.

**Figure 11.** Measured flyrock and predicted values of flyrock with the PSO-ELM for training data.

**Figure 10.** Measured flyrock and predicted values of flyrock with the PSO-ANN for training data.

**Figure 11. Figure 11.**  Measured flyrock and predicted values of flyrock with the PSO-ELM for training data. Measured flyrock and predicted values of flyrock with the PSO-ELM for training data.

**Figure 12.** Comparison of measured flyrock vs. the EO-ELM prediction for testing data. **Figure 12.** Comparison of measured flyrock vs. the EO-ELM prediction for testing data. **Figure 12.** Comparison of measured flyrock vs. the EO-ELM prediction for testing data.

**Figure 13.** Comparison of measured flyrock vs. the PSO-ANN prediction for testing data. **Figure 13.** Comparison of measured flyrock vs. the PSO-ANN prediction for testing data. **Figure 13.** Comparison of measured flyrock vs. the PSO-ANN prediction for testing data.

**Figure 14.** Comparison of measured flyrock vs. the PSO-ELM prediction for testing data.

**Figure 14.** Comparison of measured flyrock vs. the PSO-ELM prediction for testing data.

**Figure 13.** Comparison of measured flyrock vs. the PSO-ANN prediction for testing data.

**Figure 12.** Comparison of measured flyrock vs. the EO-ELM prediction for testing data.

**Figure 14. Figure 14.** Comparison of measured flyrock vs. the PSO Comparison of measured flyrock vs. the PSO-ELM prediction for testing data. -ELM prediction for testing data. 17.60, NSE = 0.978, VAF = 97.88, A20 = 0.65) performed better compared to the PSO-ELM

**Table 3.** Linear equations for predicted and measured values for the EO-ELM, PSO-ANN, and PSO-ELM for training and testing data. (R2 = 0.959, RMSE = 35.7 MAE = 23.53, MAPE = 21.84, NSE = 0.96, VAF = 95.79, A20 = 0.56) and PSO-ANN (R2 = 0.924, RMSE = 48.12, MAE = 31.68, MAPE = 24.25, NSE = 0.93, VAF = 92.89, A20 = 0.35). Furthermore, for better representation in terms of model deviations, the


EO-ELM shows the minimum deviations followed by the PSO-ELM and PSO-ANN

Table 4 shows the better prediction efficiency of the EO-ELM in the training and testing period compared to the PSO-ANN and PSO-ELM in terms of seven matrices. In the testing period, the developed EO-ELM (R2 = 0.97, RMSE = 34.82, MAE = 20.3, MAPE = 17.60, NSE = 0.978, VAF = 97.88, A20 = 0.65) performed better compared to the PSO-ELM (R2 = 0.959, RMSE = 35.7 MAE = 23.53, MAPE = 21.84, NSE = 0.96, VAF = 95.79, A20 = 0.56) and PSO-ANN (R2 = 0.924, RMSE = 48.12, MAE = 31.68, MAPE = 24.25, NSE = 0.93, VAF = 92.89, A20 = 0.35). Furthermore, for better representation in terms of model deviations, the receiver operating characteristic (ROC) curve was drawn. It is evident that all of the models capture the good relationship in the prediction of flyrock during training (Figure 15) and minimum deviation was found for the EO-ELM followed by the PSO-ELM and PSO-ANN. Amongst the models, the EO-ELM shows comparatively lesser deviation during the training period. During the testing period, a similar pattern was observed, the EO-ELM shows the minimum deviations followed by the PSO-ELM and PSO-ANN (Figure 16). (Figure 16). **Table 4.** Comparison of three models in terms of seven matrices **Training Data Sets R<sup>2</sup> RMSE MAE MAPE NSE VAF A20** EO-ELM 0.942 17.02 11.26 21.20 0.946 94.62 0.53 PSO-ANN 0.827 29.5 21.07 27.04 0.821 82.19 0.43 PSO-ELM 0.907 21.56 15.64 24.27 0.900 90.08 0.45 **Testing Data Sets (Continued)** R<sup>2</sup> RMSE MAE MAPE NSE VAF A20 EO-ELM 0.973 34.82 20.3 17.60 0.978 97.88 0.65 PSO-ANN 0.924 48.12 31.68 24.25 0.93 92.89 0.35 PSO-ELM 0.959 35.7 23.53 21.84 0.96 95.79 0.56

**Figure 15.** Diagram of REC for training dataset. **Figure 15.** Diagram of REC for training dataset.

Table 4 shows the better prediction efficiency of the EO-ELM in the training and testing period compared to the PSO-ANN and PSO-ELM in terms of seven matrices. In the testing period, the developed EO-ELM (R2 = 0.97, RMSE = 34.82, MAE = 20.3, MAPE = 17.60, NSE = 0.978, VAF = 97.88, A20 = 0.65) performed better compared to the PSO-ELM (R2 = 0.959, RMSE = 35.7 MAE = 23.53, MAPE = 21.84, NSE = 0.96, VAF = 95.79, A20 = 0.56) and PSO-ANN (R2 = 0.924, RMSE = 48.12, MAE = 31.68, MAPE = 24.25, NSE = 0.93, VAF = 92.89, A20 = 0.35). Furthermore, for better representation in terms of model deviations, the receiver operating characteristic (ROC) curve was drawn. It is evident that all of the models capture the good relationship in the prediction of flyrock during training (Figure 15) and minimum deviation was found for the EO-ELM followed by the PSO-ELM and PSO-ANN. Amongst the models, the EO-ELM shows comparatively lesser deviation during the training period. During the testing period, a similar pattern was observed, the EO-ELM shows the minimum deviations followed by the PSO-ELM and PSO-ANN

**Figure 16.** Diagram of REC for testing dataset.

**Figure 15.** Diagram of REC for training dataset.

(Figure 16).

**Table 4.** Comparison of three models in terms of seven matrices.

**Table 4.** Comparison of three models in terms of seven matrices

EO-ELM 0.942 17.02 11.26 21.20 0.946 94.62 0.53 PSO-ANN 0.827 29.5 21.07 27.04 0.821 82.19 0.43 PSO-ELM 0.907 21.56 15.64 24.27 0.900 90.08 0.45

EO-ELM 0.973 34.82 20.3 17.60 0.978 97.88 0.65 PSO-ANN 0.924 48.12 31.68 24.25 0.93 92.89 0.35 PSO-ELM 0.959 35.7 23.53 21.84 0.96 95.79 0.56

**Training Data Sets R<sup>2</sup> RMSE MAE MAPE NSE VAF A20**

**Testing Data Sets (Continued)** R<sup>2</sup> RMSE MAE MAPE NSE VAF A20


#### *5.1. Average Performance of Models*

Table 5 shows the average results for the 10-times run of three models. The EO-ELM shows the best average prediction performance compared to the PSO-ELM and PSO-ANN in terms of all matrices (Table 5). Figure 17a shows the most generalized performance of the EO-ELM (training phase) at each of the runs (10-times) compared to the PSO-ELM and PSO-ANN. Figure 17b shows that the EO-ELM has the best average convergence rate compared to the PSO-ELM and PSO-ANN.

**Table 5.** Comparison of average results for 1-times run of models.


ANN.

Table 5 shows the average results for the 10-times run of three models. The EO-ELM shows the best average prediction performance compared to the PSO-ELM and PSO-ANN in terms of all matrices (Table 5). Figure 17a shows the most generalized performance of the EO-ELM (training phase) at each of the runs (10-times) compared to the PSO-ELM and PSO-ANN. Figure 17b shows that the EO-ELM has the best average convergence rate compared to the PSO-ELM and PSO-

**Table 5.** Comparison of average results for 1-times run of models.

**Training Data Sets R<sup>2</sup> RMSE MAE MAPE NSE VAF A20**

**Testing Data Sets (Continued) R<sup>2</sup> RMSE MAE MAPE NSE VAF A20**

EO-ELM 0.95 16.66 12.13 19.87 0.95 94.46 0.60 PSO-ANN 0.83 29.68 19.33 29.30 0.82 82.06 0.41 PSO-ELM 0.88 23.68 16.76 26.96 0.88 88.29 0.47

EO-ELM 0.97 32.14 19.78 20.37 0.93 93.97 0.57

**Figure 17.** Convergence of the models for ( **Figure 17.** Convergence of the models for ( **a**) 10 runs, and (**ab** ) 10 runs, and ( ) average of 10 **b**) average of 10 runs. runs.

### *5.2. Anderson–Darling (A–D) Test*

A non-parametric test called the A–D test was performed to assess the normality of all three models [68]. The *p*-values for the PSO-ELM, PSO-ANN, and EO-ELM models are less than the significance level of 0.05 (Table 6). Table 6 shows that the EO-ELM is the best performing model in estimating flyrock.

**Table 6.** Comparison of average results for 10-times run of models.


#### *5.3. Sensitivity Analysis* whereas, *zi* is a vector of length m in array *Z*, that as Equation (9):

*5.3. Sensitivity Analysis*

Sensitivity analysis of parameters, with respect to the weathering index, was carried out as shown in Figure 18. It showed the corresponding relationship of the parameter with respect to the measured flyrock distance based on the cosine amplitude. The application of this method was based on expressing all data pairs in a common *Z*-space. The following equation defines a data array *Z* based on data pairs of each input and output: = {1 , 2 , 3 … , }. 1 = {1 , 2 , 3 … , } (9) Each point required *m*-coordinates for describing completely by training each of these data pairs in m dimensional space. The results were achieved as all of the points were in the spaced pair. The following equation shows the strength of relation (*rij*) between the data set *Z<sup>i</sup>* and *Z<sup>j</sup>* as Equation (10):

2

A non-parametric test called the A–D test was performed to assess the normality of all three models [68]. The *p*-values for the PSO-ELM, PSO-ANN, and EO-ELM models are less than the significance level of 0.05 (Table 6). Table 6 shows that the EO-ELM is the best

**Count Mean Median SD AD** *p***-Value**

Sensitivity analysis of parameters, with respect to the weathering index, was carried out as shown in Figure 18. It showed the corresponding relationship of the parameter with respect to the measured flyrock distance based on the cosine amplitude. The application of this method was based on expressing all data pairs in a common *Z*-space. The following

, }. = {<sup>1</sup>

$$Z = \{z\_1, z\_2, z\_3, z\_{\dots}, z\_{\text{i}}, z\_{\text{n}}\}.\\Z = \{Z\_1, Z\_2, Z\_3, \dots, Z\_{\text{i}}, \dots, Z\_{\text{n}}\}\tag{8}$$

, <sup>2</sup>

, <sup>3</sup> … ,

, … , } (8)

whereas, *z<sup>i</sup>* is a vector of length m in array *Z*, that as Equation (9): √∑ =1 ∑ =1

= {<sup>1</sup>

and the measured value of flyrock distance.

, <sup>2</sup> , <sup>3</sup>

*Sustainability* **2022**, *14*, x FOR PEER REVIEW 21 of 26

**Table 6.** Comparison of average results for 10-times run of models.

Actual 114 81.307 50.5 85.927 0 1 PSO-ELM 114 78.951 54.632 75.877 2.843 0.03244 PSO-ANN 114 78.183 55.035 70.484 2.308 0.00619 EO-ELM 114 79.311 53.492 76.092 0.8886 0.004215

equation defines a data array *Z* based on data pairs of each input and output:

, … ,

*5.2. Anderson–Darling (A–D) Test*

performing model in estimating flyrock.

$$z\_{i} = \{z\_{i1}, z\_{i2}, z\_{i3}, \dots, z\_{l}\}.\\z\_{i1} = \{Z\_{i1}, Z\_{i2}, Z\_{i3}, \dots, Z\_{l}\} \tag{9}$$

2

**Figure 18. Figure 18.**  Sensitivity analysis of parameters with respect to weathering index. Sensitivity analysis of parameters with respect to weathering index.

Each point required *m*-coordinates for describing completely by training each of these data pairs in m dimensional space. The results were achieved as all of the points were in the spaced pair. The following equation shows the strength of relation (*rij*) between the data set *Z<sup>i</sup>* and *Z<sup>j</sup>* as Equation (10):

$$r\_{ij} = \frac{\sum\_{k=1}^{m} z\_{ik} z\_{jk}}{\sqrt{\sum\_{k=1}^{m} z\_{ik}^2 \sum\_{k=1}^{m} z\_k^2}}\tag{10}$$

The input parameters were selected, which were most sensitive, to apply to various prediction models and identify the best model suitable for comparing the predicted value and the measured value of flyrock distance.

#### **6. Conclusions**

This study uses three hybrid models: the EO–ELM, PSO–ANN, and PSO-ELM to predict flyrock. Out of these three hybrid models, the EO-ELM is proposed and the rest were used to validate the performance of the proposed model. The seven different matrices (*R* 2 , RMSE, MAPE, NSE, MAE, VAF, and A20) are used for comparing the efficacy of the developed model. The developed EO-ELM model performed better compared to PSOthe ELM and PSO-ANN in predicting flyrock (Table 4). Further, all models were run 10-times and average results are shown in Table 5. It was observed that the EO-ELM model outperformed the PSO-ELM and PSO-ANN in average results. Furthermore, the 10-times

run of the EO-ELM model showed better convergence capability (Figure 17) than the other two. Further, the A–D test showed that the EO-ELM model has better performance efficiency compared to the PSO-ELM and PSO-ANN. A sensitivity analysis was done introducing a new parameter, WI. The PF and BI showed the highest sensitivity with 0.98 each (Figure 18).

The limitation of this study is that this study was carried out for a particular limestone mine in Thailand. Therefore, the obtained results may not be suitable for other geological settings, i.e., in granite or other quarries. So, there is a need of further research by considering non controllable parameters of the blastability index, WI. On the other hand, by refining the controllable parameters, accuracy of the prediction of flyrock can be improved. There are a large number of mines near to the limestone mine under study and if large data sets are collected, further reliability of the prediction models can be improved through future research. The number of input parameters is eight. However, future studies of the prediction of flyrock can be done by limiting to five influential input parameters. The use of the latest technology, such as the video recording of flyrock with drones, will add value for future research. Furthermore, the wavelet frequency domain parameter or innovative liquid carbon dioxide rock-breaking technology are recent technologies used for blasting. Thus, there is a need to develop a new technology alternative to blasting.

**Author Contributions:** Conceptualization, E.T.M., R.M.B., D.K. and B.R.; methodology, R.M.B. and R.K.; software, R.K., D.K. and B.R.; formal analysis, R.M.B. and D.K.; resources, E.T.M. and R.M.B.; data curation, R.M.B.; writing—original draft, E.T.M., R.M.B., D.K., BR., M.M.S.S., S.K. and R.K.; writing—review and editing, R.MB., E.T.M., RK., M.M.S.S. and D.K.; supervision, E.T.M., D.K. and B.R.; funding acquisition, M.M.S.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research is partially funded by the Ministry of Science and Higher Education of the Russian Federation under the strategic academic leadership program 'Priority 2030' (Agreement 075-15-2021-1333 dated 30 September 2021).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** Authors are thankful to Edy Toonizam Mohamad, Dean of Department of Civil Engineering, Universiti Teknologi Malaysia for encouragement for this research. Authors are also thankful to Vitaly Sergeev, Centre of Peter the Great St. Petersburg Polytechnic University, 195251 St. Petersburg, Russia for guidance and supervision of this paper.

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

### **References**


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

**Yapeng Chen <sup>1</sup> , Tong Wu <sup>1</sup> , Xiaoshi Yan <sup>1</sup> , Shang Shi <sup>2</sup> , Jianyong Li <sup>2</sup> and Jinyu Dong 2,\***

<sup>1</sup> Beifang Investigation, Design & Research Co., Ltd., Tianjin 300222, China


**Abstract:** Based on a typical project in an altered rock area, this study carried out numerical simulations using the FLAC3D software to calculate the changes in the stress field, deformation field, and plastic zone of the surrounding rock during the unsupported and supported excavation of a water transfer tunnel. The degree of alteration of the surrounding rock was considered as the base point. The following results were obtained: in the unsupported state, the tunnel surrounding rock was affected by different degrees of alteration, and compressive stress concentration appeared within a certain range at the bottom of the chamber. The value of all-directional stress decreased with the deepening of the degree of alteration, while the opposite was the case for the depth of influence. The displacement changes at the bottom and side walls of the chamber were large and increased significantly with the deepening of the degree of alteration; the displacement monitoring points distributed around the tunnel exhibited the same deformation trend. The plastic zone of the surrounding rock obviously expanded as the degree of alteration deepened. The stress, deformation field, and plastic zone of the tunnel surrounding rock were effectively controlled after the adoption of support measures. The results obtained by this study can be used as a reference for similar projects in altered rock areas.

**Keywords:** alteration; tunnel; numerical simulation; stress field; deformation field; plastic zone

### **1. Introduction**

To solve the problems of the uneven distribution of water resources and the contradiction between supply and demand, many large-scale water conservancy construction projects have been built around the world, among which water transmission tunnels are a key component with high importance. Extreme engineering geological conditions, such as significant faults, seasonal thawing in permafrost regions, and high-pressure water action, frequently cause the issue of surrounding rock deformation and damage during tunnel excavation [1–4]. Similarly, the stability of the surrounding rock is also significantly impacted by the complicated stress changes, rock extrusion, and deformation brought on by tunnel excavation. Dong et al. [5] exposed how tectonic stress forces affected the rock surrounding the tunnel's stress and deformation damage pattern. After conducting an excavation simulation under high ground stress for the underground chamber complex of Jinping I hydropower facility, Qian and Zhou [6] discovered that the rock body will manifest as two-dimensional band disintegration phenomenon and suggested corresponding support measures. After examining the primary stress rotation mechanism and the rock extrusion and deformation law during the excavation of high and deep buried tunnels, Cai et al. [7,8] pointed out that the three-dimensional spatial effect is more significant for the analysis of the stability of the surrounding rock. They proposed a three-dimensional orthotropic analysis and rock strength based on the GZZ strength criterion that can weaken the three-dimensional spatial effect and exert the rock body's own strength.

**Citation:** Chen, Y.; Wu, T.; Yan, X.; Shi, S.; Li, J.; Dong, J. Stress and Deformation Characteristics of Tunnel Surrounding Rock under Alteration. *Sustainability* **2023**, *15*, 1161. https://doi.org/10.3390/ su15021161

Academic Editors: Marc A. Rosen, Mahdi Hasanipanah, Danial Jahed Armaghani and Jian Zhou

Received: 13 September 2022 Revised: 24 December 2022 Accepted: 28 December 2022 Published: 7 January 2023

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

Additionally, as a particular class of engineering geological problems, alteration rocks have been exposed in several projects around the world, such as the Sanjiang orogenic section of the Yunnan–Tibet Railway in China [9], the Kerman Tunnel in the Urmia–Dokhtar Magmatic Arc (UDMA) in central Iran [10], multiple geothermal power plants in the Kuril–Kamchatka island arc, Russia [11], and the Wheal Martyn china pit in Southwest England [12]. Alteration rocks are a new class of rock that forms after diagenesis by hydrothermal erosion, tectonic dynamics, and secondary weathering that alters the elemental composition and structural features of the original rock to various degrees [13,14]. Researchers have examined the composition and characteristics of altered rocks in various geological contexts and discovered that some feldspar and mica minerals are mostly changed into clay minerals, sericite, chlorite, and chlorite after alteration [15–17]. The destruction of the original internal structure of the rocks caused by the alteration of tiny minerals increases the internal porosity of the rocks. The development of internal porosity as well as microcracks in the rock will produce continuity fracture damage after being disturbed by excavation [18,19]. Along with the long-term geological tectonic evolution and the effect of ground stress in the region, the distribution of alteration rocks is irregular, the internal structure is highly fragmented and the physical and mechanical properties are poor [20–23], which have a direct impact on site selection and the design and construction of projects.

The previous research system to examine the mechanical behavior of deeply buried tunnels from the viewpoints of intrinsic model and geological structure is reasonably well developed; the theoretical research on micro mineral analysis, physical and mechanical property changes, and alteration degree classification of altered rocks is also reasonably mature [24,25], but the pertinent engineering case studies are slightly lacking. Specifically in deeply buried tunnels, where the damage to the mechanical properties of rocks due to alteration is frequently beyond our original prediction, it has been common to observe construction challenges, schedule delays, and even safety threats brought on by improper support and poor response in the field. This is because there is little consideration of the alteration effect of rocks during engineering construction in altered rock areas. Therefore, there is scientific and practical value in analyzing the stress and deformation damage characteristics of tunnel excavation surrounding rocks under the effect of alteration for the design and construction of projects in altered rock areas. This study considered a water transfer tunnel with different degrees of alteration as an example. The FLAC3D software was used to numerically simulate the excavation of a typical alteration tunnel section in its natural state and under two working conditions after support with the objective of analyzing the stress and deformation damage characteristics of the surrounding rock under different alteration degrees and elucidating the actual impact of alteration on the surrounding rock stability. The findings of this study provide the theoretical basis for optimizing the support scheme, and a reference for the design and construction of similar projects.

#### **2. Engineering Geological Conditions**

The water transmission tunnel is located in the central and western part of North Tianshan. This tunnel has a total length of 41.82 km, diameter of 5.3 m, and longitudinal slope of 1/564.8, which means that it is a deep and long buried tunnel. The study area has a complex geological structure with several northwest and northeast-trending compression–torsional faults, fold zones, and extrusion fracture zones. The tunnel's surrounding rocks belong to various lithologies, and mainly include Silurian, Devonian, Carboniferous sandstone, metamorphic sandstone, tuff, tuffaceous sandstone, and Hualixi-age granite. Among them, the granite section is 9.81 km long, dominated by diorite, granodiorite, and potassium granite, and was formed by the crustal movement of the Late Paleozoic and magmatic activity of the Hualixi period. After a long and complex tectonic–hydrothermal superposition modification, the granites in the study area have generally been affected by alteration,

with the chloritization of black mica and clayification of feldspar minerals as the main alteration types. original structure, and these rocks can be crushed by hand; the rebound value is less than 10, and the cavern is in a large area of debris collapse, exhibiting obvious plastic deformation when wet. The alteration of the surrounding rock in the tunnel is widely

more collapsed and larger in size, and cavities can easily form at the top arch. The structure of strongly altered rocks is completely destroyed, with visible traces of the

Based on a field survey and field test, the degree of rock alteration can be divided into three categories: slight alteration, moderate alteration, and strong alteration. Slightly altered rocks have a relatively intact structure, produce a brittle hammering sound, and their rebound values are between 30 and 50; the cave blocks are less collapsed and larger in size. Moderately altered rocks have a partially broken structure, produce a muffled hammering sound, and their rebound values are between 10 and 30; the cave blocks are more collapsed and larger in size, and cavities can easily form at the top arch. The structure of strongly altered rocks is completely destroyed, with visible traces of the original structure, and these rocks can be crushed by hand; the rebound value is less than 10, and the cavern is in a large area of debris collapse, exhibiting obvious plastic deformation when wet. The alteration of the surrounding rock in the tunnel is widely distributed, and the construction process often leads to the large deformation of the surrounding rock, collapse, and other geological problems. distributed, and the construction process often leads to the large deformation of the surrounding rock, collapse, and other geological problems. **3. Numerical Computational Model**  *3.1. Computational Model*  Based on site investigation and geological data, the geological model of the study area was generalized, and a three-dimensional (3D) numerical calculation tunnel model was established based on FLAC3D. As shown in Figure 1a, the water transmission tunnel model had an 8 m diameter, the height of the model was approximately 1313 m in the vertical direction (Z-axis), the calculated elevation at the bottom was 1000 m, and the highest elevation at the top surface was 2313 m. The model width was 200 m (X-axis direction), and 100 m was taken from each side of the tunnel's centerline. The extension

#### **3. Numerical Computational Model** of the model was 600 m (Y-axis direction), including 250 m for the slightly altered rock section (mileage K32 + 105**–**355), 200 m (mileage K32 + 355**–**555) for the moderately altered

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 4 of 18

#### *3.1. Computational Model* rock section, and 150 m (mileage K32 + 555**–**705) for the strongly altered rock section.

Based on site investigation and geological data, the geological model of the study area was generalized, and a three-dimensional (3D) numerical calculation tunnel model was established based on FLAC3D. As shown in Figure 1, the water transmission tunnel model had an 8 m diameter, the height of the model was approximately 1313 m in the vertical direction (Z-axis), the calculated elevation at the bottom was 1000 m, and the highest elevation at the top surface was 2313 m. The model width was 200 m (X-axis direction), and 100 m was taken from each side of the tunnel's centerline. The extension of the model was 600 m (Y-axis direction), including 250 m for the slightly altered rock section (mileage K32 + 105–355), 200 m (mileage K32 + 355–555) for the moderately altered rock section, and 150 m (mileage K32 + 555–705) for the strongly altered rock section. The geological model was simplified to four geotechnical materials, namely, the upper low resistance overburden and the lower three granites with different degrees of alteration. The tunnel passed through the lower altered rock layer, where the medium altered rock zone contained a fault. The geological model is shown in Figure 1. Since the model's overall size was too large, transitional meshing was employed to simplify the computation while guaranteeing that the mesh size close to the tunnel satisfied the requirements for calculation accuracy. The boundary conditions of the model were the Xdirectional displacement constraint along the tunnel extension boundary, the Ydirectional displacement constraint along the vertical tunnel boundary, and the fixed constraint at the bottom boundary.

**Figure 1.** Calculation model diagram. 1. Low-resistivity overburden; 2. Water tunnel; 3. Slightly altered rock zone; 4. Moderately altered rock zone; 5. Strongly altered rock zone. **Figure 1.** Calculation model diagram. 1. Low-resistivity overburden; 2. Water tunnel; 3. Slightly altered rock zone; 4. Moderately altered rock zone; 5. Strongly altered rock zone.

*3.2. Application of Ground Stress and Selection of Calculation Parameters*  To accurately obtain the tunnel ground stress field distribution, six sets of acoustic emission Kaiser effect ground stress tests and analyses were conducted on the borehole cores of the granite section, and the measured results are presented in Table 1. The test results reveal that the measured principal stress value of the borehole increased with the depth; the maximum principal compressive stress *σ*1 was approximately horizontal, the The geological model was simplified to four geotechnical materials, namely, the upper low resistance overburden and the lower three granites with different degrees of alteration. The tunnel passed through the lower altered rock layer, where the medium altered rock zone contained a fault. The geological model is shown in Figure 1. Since the model's overall size was too large, transitional meshing was employed to simplify the computation while guaranteeing that the mesh size close to the tunnel satisfied the requirements for calculation accuracy. The boundary conditions of the model were the X-directional displacement constraint along the tunnel extension boundary, the Y-directional displacement constraint along the vertical tunnel boundary, and the fixed constraint at the bottom boundary.

#### *3.2. Application of Ground Stress and Selection of Calculation Parameters*

To accurately obtain the tunnel ground stress field distribution, six sets of acoustic emission Kaiser effect ground stress tests and analyses were conducted on the borehole cores of the granite section, and the measured results are presented in Table 1. The test results reveal that the measured principal stress value of the borehole increased with the depth; the maximum principal compressive stress *σ*<sup>1</sup> was approximately horizontal, the dip angle *a*<sup>1</sup> was within ± 10◦ , and the dominant direction was north-northwest, which was more consistent with the direction of the regional tectonic stress field. The middle principal stress *σ*<sup>2</sup> had a gentle dip angle, and the dip angle *a*<sup>2</sup> varied within ± 28◦ ; the minimum principal stress dip angle *a*<sup>3</sup> was above 64◦ .


Note: *σ*1, *σ*2, and *σ*<sup>3</sup> are the maximum principal compressive stress, intermediate principal stress, and minimum principal stress, respectively; *a*1, *a*2, and *a*<sup>3</sup> represent the angle (inclination) between the maximum principal compressive stress, intermediate principal stress, and minimum principal stress, and horizontal plane, respectively, with positive values representing the elevation angle and negative values representing the pitch angle; *β*1, *β*2, and *β*<sup>3</sup> are the angle between the projection of the three principal stresses on the oxy plane and x-axis, respectively, with positive values indicating counterclockwise rotation and negative values indicating clockwise rotation.

Based on the data in Table 1, the ground stress in the granite cave section can be projected, and the linear correlation equations between the maximum horizontal principal stress (*σ*1), horizontal intermediate principal stress (*σ*2), and minimum horizontal principal stress (*σ*3) and the burial depth (*H*) can be derived as follows:

> *σ*<sup>1</sup> = 0.0196*H* + 16.653, *R* <sup>2</sup> = 0.9304 *σ*<sup>2</sup> = 0.0211*H* + 7.0488, *R* <sup>2</sup> = 0.9688 *σ*<sup>3</sup> = 0.0303*H* − 3.3618, *R* <sup>2</sup> = 0.9862

The principal stresses in the cavern line were projected according to the linear correlation equation, and the results are shown in Table 2. In the numerical calculation, *σ*1, *σ*2, and *σ*<sup>3</sup> were transformed along the X-direction (horizontal vertical tunnel axis direction), Y-direction (tunnel axis direction), and vertical direction to apply the Sxx, Syy, and Szz stresses to the model.

**Table 2.** Projected ground stress values in the calculated section of tunnel.


The elastic–plastic model and Mohr–Coulomb strength criterion were used in the calculation. The rock surrounding parameters were determined according to the standard of International Society of Rock Mechanics (ISRM, 2007) after testing at the Quality Inspection Center of Capital Construction Project of Haihe Water Conservancy Commission, Ministry of Water Resources of China and School of Earth Science and Engineering, North China University of Water Resources and Electric Power, with reference to the Engineering Rock Quality Grading Standard (GBT50218-2014) and after considering the degree of rock alteration, as presented in Table 3.


**Table 3.** Proposed values of the main geological parameters of the tunnel envelope.

According to the excavation design plan, after the tunnel excavation, the TBM shield was closed with synthetic coarse fiber concrete in time after the initial spraying. The slightly altered section adopted HW125 steel arch racks with a distance of 0.9 m; the moderately altered section adopted HW150 steel arch racks with a distance of 0.5 m; the strongly altered section adoptsedHW150 steel arch racks with a distance of 0.3 m. The longitudinal connection adopted Φ20 steel bars with a ring spacing of 1 m, and the top arch was equipped with Φ20 reinforcement rows within 150◦ . The support parameters are shown in Table 4, Table 5, and Table 6, respectively.

**Table 4.** Basic parameters of the anchor rods.


**Table 5.** Basic parameters of the steel arch.


**Table 6.** Basic parameters of the concrete primary lining.


#### **4. Analysis of the Calculation Results**

For the in-depth investigation of the distribution and changes of the stress field and deformation field of the surrounding rock, typical sections were selected for detailed analysis in the surrounding rock with different degrees of alteration. Section K32 + 230 was selected in the slightly altered surrounding rock with a burial depth of 905 m; section K32 + 455 was selected in the moderately altered surrounding rock with a burial depth of 860 m; section K32 + 630 was selected in the strongly altered surrounding rock with a burial depth of 930 m.

#### *4.1. Analysis of Stress Field*

After the tunnel excavation, the horizontal maximum principal stress was approximately in line with the tunnel axis, ignoring its influence on the tunnel surrounding rock. The horizontal intermediate principal stress was approximately perpendicular to the tunnel

axis, and the self-weight stress was the minimum principal stress. Provided that support measures were not used, under the influence of horizontal tectonic stress, a compressive stress concentration with a maximum value of 102.4 MPa appeared in the slightly altered surrounding rock section at the top and bottom 7.0 m of the chamber's depth, and gradually decreased to the initial stress state as the distance from the chamber increased, as shown in Figure 2a. Under the influence of chamber excavation and unloading, a stress reduction zone with a value of 54.6 MPa and influence depth of 16.0 m appeared in the surrounding rock on both sides and at the top of the chamber. Shear stress concentration appeared at the top and bottom corners of the cavern, and its maximum value was 29.5 MPa, as shown in Figure 2b. Compressive stress concentration with a maximum value of 87.36 MPa appeared in the medium alteration surrounding rock section 20 m at the top and bottom corners of the cavern, and gradually decreased to the initial stress state as the distance from the cavern chamber increased; the influence depth was 28.0 m, as shown in Figure 2c. A stress reduction zone appeared in the surrounding rock on both sides of the cavern chamber and at the top and bottom; its value was 42.3 MPa and its influence depth was 32.0 m. Shear stress concentration appeared at the top and bottom corner of the cavern, and its maximum value was 19.15 MPa, as shown in Figure 2d. In the strongly altered surrounding rock section, compressive stress concentration occurred at a depth of 13.0 m from the bottom of the cavern, and had the maximum value of 107.8 MPa, as shown in Figure 2e. In addition to the bottom of the cavern, a stress reduction zone with the minimum value of 49.04 MPa and small main stress influence depth of 35 m existed around the cavern. Shear stress concentration with an influence depth of 12.0 m existed at the top and bottom corners of the cavern. The maximum value was 19.7 MPa, as shown in Figure 2f. The depth of stress influence was significantly higher in the moderately altered and strongly altered surrounding rocks compared to the slightly altered surrounding rock.

After using support measures, the slightly altered surrounding rock section exhibited compressive stress concentration in the surrounding rock at the top and bottom corners of the cavern, and its maximum value was 102.1 MPa, as shown in Figure 3a. The stress reduction zone appeared at the bottom of the cavern, and its value was 59.9 MPa. Shear stress concentration appeared at the top and bottom corners of the cavern, and its maximum value was 31.5 MPa, as shown in Figure 3b. In the moderately altered surrounding rock section, compressive stress concentration appeared at the top and bottom corners of the cavern, and its maximum value was 101.6 MPa, as shown in Figure 3c. A stress reduction zone appeared at the side walls and bottom of the cavern, and its value was 61.4 MPa. Shear stress concentration appeared at the top and bottom corners of the cavern, and its maximum value was 24.9 MPa, as shown in Figure 3d. Compressive stress concentration appeared at the top and bottom corners of the cavern section, with a maximum value of 98.3 MPa, and gradually decreased to the initial stress state as the distance from the cavern increased, as shown in Figure 3e. A stress reduction zone with a minimum value of 61.06 MPa appeared in the side walls of the cavern. Shear stress concentration appeared at the top and bottom corners of the cavern, and its maximum value was 20.46 MPa, as shown in Figure 3f.

**Figure 2.** Main stress diagram of the tunnel excavation under unsupported conditions. (**a**) Middle main stress diagram of the slightly altered surrounding rock section. (**b**) Shear stress cloud diagram of the slightly altered surrounding rock section. (**c**) Middle main stress diagram of the moderately altered surrounding rock section. (**d**) Shear stress cloud diagram of the moderately altered surrounding rock section. (**e**) Middle main stress diagram of the strongly altered surrounding rock section. (**f**) Shear stress cloud diagram of the strongly altered surrounding rock section. **Figure 2.** Main stress diagram of the tunnel excavation under unsupported conditions. (**a**) Middle main stress diagram of the slightly altered surrounding rock section. (**b**) Shear stress cloud diagram of the slightly altered surrounding rock section. (**c**) Middle main stress diagram of the moderately altered surrounding rock section. (**d**) Shear stress cloud diagram of the moderately altered surrounding rock section. (**e**) Middle main stress diagram of the strongly altered surrounding rock section. (**f**) Shear stress cloud diagram of the strongly altered surrounding rock section.

After using support measures, the slightly altered surrounding rock section exhibited compressive stress concentration in the surrounding rock at the top and bottom corners of the cavern, and its maximum value was 102.1 MPa, as shown in Figure 3a. The stress reduction zone appeared at the bottom of the cavern, and its value was 59.9 MPa. Shear stress concentration appeared at the top and bottom corners of the cavern, and its maximum value was 31.5 MPa, as shown in Figure 3b. In the moderately altered surrounding rock section, compressive stress concentration appeared at the top and bottom corners of the cavern, and its maximum value was 101.6 MPa, as shown in Figure 3c. A stress reduction zone appeared at the side walls and bottom of the cavern, and its value was 61.4 MPa. Shear stress concentration appeared at the top and bottom corners of the cavern, and its maximum value was 24.9 MPa, as shown in Figure 3d. Compressive stress concentration appeared at the top and bottom corners of the cavern section, with a maximum value of 98.3 MPa, and gradually decreased to the initial stress state as the distance from the cavern increased, as shown in Figure 3e. A stress reduction zone with a minimum value of 61.06 in Figure 3f.

structure.

**Figure 3.** Main stress diagram of the excavation after the tunnel support. (**a**) Middle main stress diagram of the slightly altered surrounding rock section. (**b**) Shear stress cloud diagram of the slightly altered surrounding rock section. (**c**) Middle main stress diagram of the moderately altered surrounding rock section. (**d**) Shear stress cloud diagram of the moderately altered surrounding rock section. (**e**) Middle main stress diagram of the strongly altered surrounding rock section. (**f**) Shear stress cloud diagram of the strongly altered surrounding rock section. **Figure 3.** Main stress diagram of the excavation after the tunnel support. (**a**) Middle main stress diagram of the slightly altered surrounding rock section. (**b**) Shear stress cloud diagram of the slightly altered surrounding rock section. (**c**) Middle main stress diagram of the moderately altered surrounding rock section. (**d**) Shear stress cloud diagram of the moderately altered surrounding rock section. (**e**) Middle main stress diagram of the strongly altered surrounding rock section. (**f**) Shear stress cloud diagram of the strongly altered surrounding rock section.

MPa appeared in the side walls of the cavern. Shear stress concentration appeared at the top and bottom corners of the cavern, and its maximum value was 20.46 MPa, as shown

The comparison between the stresses in the unsupported and supported cavern chambers is shown in Figure 4. As can be seen from the support, except for the strong alteration surrounding the rock section, owing to the occurrence of surrounding rock damage deformation and X-direction compressive stress reduction, the other stresses were increased by the surrounding rock owing to the restraining effect of the support

The comparison between the stresses in the unsupported and supported cavern chambers is shown in Figure 4. As can be seen from the support, except for the strong alteration surrounding the rock section, owing to the occurrence of surrounding rock damage deformation and X-direction compressive stress reduction, the other stresses were increased by the surrounding rock owing to the restraining effect of the support structure.

**Figure 4.** Comparison of surrounding rock stress in each profile before and after excavation and support of alteration rock tunnel. **Figure 4.** Comparison of surrounding rock stress in each profile before and after excavation and support of alteration rock tunnel.

#### *4.2. Analysis of Surrounding Rock Displacement Field*

*4.2. Analysis of Surrounding Rock Displacement Field* In the natural state, the mountain displacement gradually decreases from top to bottom, with a stable trend in general, and the maximum displacement occurs in an area with large surface elevation and obvious surface undulation. The numerical simulation conducted by this study mainly considered the surrounding rock deformation caused by tunnel excavation. Therefore, the displacement in the natural state was considered to be zero, and the relative displacement of the surrounding rock after tunnel excavation was inves-In the natural state, the mountain displacement gradually decreases from top to bottom, with a stable trend in general, and the maximum displacement occurs in an area with large surface elevation and obvious surface undulation. The numerical simulation conducted by this study mainly considered the surrounding rock deformation caused by tunnel excavation. Therefore, the displacement in the natural state was considered to be zero, and the relative displacement of the surrounding rock after tunnel excavation was investigated on this basis.

tigated on this basis. In the tunnel excavation without support measures, the overall deformation of the surrounding rock after tunnel excavation was directed toward the cavern face; the top arch sank, the bottom slab rose, and the side walls moved inward. Figure 5 shows the displacement trend of each typical section of the surrounding rock under unsupported excavation. As can be seen, the maximum total displacement of the surrounding rock in the slightly altered section was 9.64 cm, and the displacement and impact depth of the top slab and sidewall were relatively large, while the displacement and impact depth of the bottom slab were relatively small. The deformation of the surrounding rock in the moderately altered section was larger than that in the slightly altered section, and the maximum total displacement was 166.4 cm. The maximum total displacement of the surrounding rock in the strongly altered section was 45.6 m, and the displacement at the base plate and sidewall was the most variable. In the displacement cloud distribution of each typical section (Figure 5), the area of the surrounding rock deformation caused by cavern excavation was pie shaped. Because of the large burial depth, the displacements of the surrounding rock in the vertical direction (Z-direction) and vertical tunnel direction (X-direction) were relatively large owing to the influence of self-weight stress and horizontal tectonic stress, while the displacement in the tunnel extension direction (Y-direction) was small. Owing to the influence of the degree of rock erosion, the displacement of the and moderately altered section (section 32 + 455, depth of 860 m) was much larger than that of In the tunnel excavation without support measures, the overall deformation of the surrounding rock after tunnel excavation was directed toward the cavern face; the top arch sank, the bottom slab rose, and the side walls moved inward. Figure 5 shows the displacement trend of each typical section of the surrounding rock under unsupported excavation. As can be seen, the maximum total displacement of the surrounding rock in the slightly altered section was 9.64 cm, and the displacement and impact depth of the top slab and sidewall were relatively large, while the displacement and impact depth of the bottom slab were relatively small. The deformation of the surrounding rock in the moderately altered section was larger than that in the slightly altered section, and the maximum total displacement was 166.4 cm. The maximum total displacement of the surrounding rock in the strongly altered section was 45.6 m, and the displacement at the base plate and sidewall was the most variable. In the displacement cloud distribution of each typical section (Figure 5), the area of the surrounding rock deformation caused by cavern excavation was pie shaped. Because of the large burial depth, the displacements of the surrounding rock in the vertical direction (Z-direction) and vertical tunnel direction (X-direction) were relatively large owing to the influence of self-weight stress and horizontal tectonic stress, while the displacement in the tunnel extension direction (Y-direction) was small. Owing to the influence of the degree of rock erosion, the displacement of the and moderately altered section (section 32 + 455, depth of 860 m) was much larger than that of the slightly altered section (section 32 + 230, depth of 905 m). Considering the discontinuity characteristics of

the rock, it is concluded that large deformation damage occurs frequently. The strongly altered section (section 32 + 630, depth of 930 m) had actually been completely destroyed by large deformation.

**Figure 5.** Displacement diagram of the tunnel excavation under unsupported conditions. (**a**) Total displacement diagram of the slightly altered surrounding rock section. (**b**) Total displacement diagram of the moderately altered surrounding rock section. (**c**) Total displacement diagram of the strongly altered surrounding rock section.

After adopting support measures, such as the anchors of a grouting reinforcement system, steel arch, shotcrete, and reinforcement mesh, the deformation area of the surrounding rock caused by tunnel excavation became butterfly shaped, and the deformation area of the surrounding rock was mainly concentrated at the side wall, whose displacement and influence depth were relatively large while the displacement and influence depth of the top arch and bottom slab were relatively small. Figure 6 shows the typical section displacement changes of each surrounding rock under supported excavation. As can been seen, the maximum total displacement of the slightly altered surrounding rock section was 1.02 cm, the maximum total displacement of the moderately altered surrounding rock section was 1.56 cm, and the maximum total displacement of the strongly altered surrounding rock section was 2.62 cm. The main reason for this is that the horizontal structural stress was the maximum main stress and the vertical stress was the minimum main stress in the tunnel. This led to the displacement of the surrounding rock, which was caused by the lateral expansion effect as the main contradiction, and the displacement of the top and bottom slab as the secondary contradiction, coupled with a high degree of fragmentation, extremely low mechanical strength, and the poor integrity of the strongly altered rock body, which led to the downward movement of the bottom slab under the influence of excavation disturbance and the constraint of the support effect during the excavation process (Figure 7a,c,e). *Sustainability* **2023**, *15*, x FOR PEER REVIEW 12 of 18

**Figure 6.** Comparison of the displacement in each profile after the excavation of the tunnel without support and after support. **Figure 6.** Comparison of the displacement in each profile after the excavation of the tunnel without support and after support.

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

support and after support.

**Figure 7.** Excavation displacement map after tunnel support. (**a**) Total displacement diagram of the slightly altered surrounding rock section. (**b**) Displacement monitoring map of the slightly altered surrounding rock section. (**c**) Total displacement diagram of the moderately altered surrounding rock section. (**d**) Displacement monitoring map of the moderately altered surrounding rock section. (**e**) Total displacement diagram of the strongly altered surrounding rock section. (**f**) Displacement monitoring map of the strongly altered surrounding rock section. **Figure 7.** Excavation displacement map after tunnel support. (**a**) Total displacement diagram of the slightly altered surrounding rock section. (**b**) Displacement monitoring map of the slightly altered surrounding rock section. (**c**) Total displacement diagram of the moderately altered surrounding rock section. (**d**) Displacement monitoring map of the moderately altered surrounding rock section. (**e**) Total displacement diagram of the strongly altered surrounding rock section. (**f**) Displacement monitoring map of the strongly altered surrounding rock section.

**Figure 6.** Comparison of the displacement in each profile after the excavation of the tunnel without

*4.3. Analysis of the Plastic Zone of Surrounding Rock* The size and distribution characteristics of the plastic zone reflect the mechanical properties of the surrounding rock. Additionally, they characterize the actual size of the loosening zone of the tunnel surrounding rock after excavation and unloading, and the degree of disturbance to the surrounding rock in each excavation phase. After excavation, because of the increase in the large principal stress and decrease in small principal stress in the surrounding rock, the surrounding rock unit was in the form of compression–shear damage; the distribution of its plastic zone is shown in Figure 8. In the natural state, the damage depth of the plastic zone was approximately 8.0–10.0 m for slightly altered surrounding rock (section 32 + 230, depth of 905 m), approximately 20.0 m for moderately altered surrounding rock (section 32 + 455, depth of 860 m), and the strongly altered section of the perimeter rock had been completely destroyed and the In the calculation process, to monitor the deformation evolution process of each type of surrounding rock after tunnel excavation, displacement monitoring points were set at the top plate, both walls, and bottom plate of the cavern; the monitoring results are shown in Figure 7b,d,f. It can be seen that after the use of support measures, the displacement monitoring amount of all types of surrounding rocks in the cavern was significantly reduced, and the surrounding rocks of the tunnel exhibited the largest displacement at the cavern side walls, followed by the bottom slab and smallest top slab. Among them, the displacement of the top plate of the slightly altered surrounding rock was 0.23 cm, the maximum displacement of both walls was 0.86 cm, and the displacement of the bottom plate was 0.48 cm. The maximum displacement of the top plate of the moderately altered surrounding rock was 0.32 cm, the maximum displacement of both walls was 1.30 cm, and the maximum displacement of the bottom plate was 0.88 cm. The displacement of the top

depth of impact was too great, after tunnel excavation, as shown in Figure 8a,c,e.

After the use of support measures, the extent of the plastic zone in the surrounding rock decreased significantly compared with the unsupported condition, and the depth of

in the moderately altered surrounding rock (section 32 + 455, depth of 860 m) decreased to approximately 1.0–3.0 m, and the depth of damage of the plastic zone in the strongly altered surrounding rock (section 32 + 630, depth of 930 m) decreased to approximately 0.6–1.0 m. The depth of damage in the plastic zone (section 32 + 630, depth of 930 m) de-

creased to approximately 1.0–3.0 m, as shown in Figure 8b,d,f.

plate of the strongly altered surrounding rock was 0.56 cm, the maximum displacement of both walls was 2.19 cm, and the displacement of the bottom slab was 1.62 cm.

### *4.3. Analysis of the Plastic Zone of Surrounding Rock*

The size and distribution characteristics of the plastic zone reflect the mechanical properties of the surrounding rock. Additionally, they characterize the actual size of the loosening zone of the tunnel surrounding rock after excavation and unloading, and the degree of disturbance to the surrounding rock in each excavation phase. After excavation, because of the increase in the large principal stress and decrease in small principal stress in the surrounding rock, the surrounding rock unit was in the form of compression–shear damage; the distribution of its plastic zone is shown in Figure 8. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 14 of 18

**Figure 8.** Cloud map of the plastic zone after tunnel excavation. (**a**) Cloud map of the plastic zone of the slightly altered rock section (without support). (**b**) Cloud map of the plastic zone of the slightly altered rock section (after support). (**c**) Cloud map of the plastic zone of the moderately altered rock section (without support). (**d**) Cloud map of the plastic zone of the moderately altered rock section (after support). (**e**) Cloud map of the plastic zone of the strongly altered rock section (without support). (**f**) Cloud map of the plastic zone of the strongly altered rock section (after support). **Figure 8.** Cloud map of the plastic zone after tunnel excavation. (**a**) Cloud map of the plastic zone of the slightly altered rock section (without support). (**b**) Cloud map of the plastic zone of the slightly altered rock section (after support). (**c**) Cloud map of the plastic zone of the moderately altered rock section (without support). (**d**) Cloud map of the plastic zone of the moderately altered rock section (after support). (**e**) Cloud map of the plastic zone of the strongly altered rock section (without support). (**f**) Cloud map of the plastic zone of the strongly altered rock section (after support).

During the tunnel boring process, it was found that the deformation of the steel arch

monitoring points with different degrees of alteration were randomly selected for each

*4.4. Recheck of the Surrounding Rock Deformation*

In the natural state, the damage depth of the plastic zone was approximately 8.0–10.0 m for slightly altered surrounding rock (section 32 + 230, depth of 905 m), approximately 20.0 m for moderately altered surrounding rock (section 32 + 455, depth of 860 m), and the strongly altered section of the perimeter rock had been completely destroyed and the depth of impact was too great, after tunnel excavation, as shown in Figure 8a,c,e.

After the use of support measures, the extent of the plastic zone in the surrounding rock decreased significantly compared with the unsupported condition, and the depth of damage of the plastic zone in the slightly altered surrounding rock (section 32 + 230, depth of 905 m) decreased to approximately 1.0–2.0 m. The depth of damage of the plastic zone in the moderately altered surrounding rock (section 32 + 455, depth of 860 m) decreased to approximately 1.0–3.0 m, and the depth of damage of the plastic zone in the strongly altered surrounding rock (section 32 + 630, depth of 930 m) decreased to approximately 0.6–1.0 m. The depth of damage in the plastic zone (section 32 + 630, depth of 930 m) decreased to approximately 1.0–3.0 m, as shown in Figure 8b,d,f.

#### *4.4. Recheck of the Surrounding Rock Deformation*

During the tunnel boring process, it was found that the deformation of the steel arch was mainly concentrated within 24 h after excavation, after which it gradually stabilized, and the final deformation was less different to the 24 h monitoring value. Therefore, 14 monitoring points with different degrees of alteration were randomly selected for each section of the excavation, and the above-mentioned support method was used to monitor and count the displacement changes of the steel arch after 24 h. The actual monitoring values are presented in Table 7 and Figure 9. Among them, the deformation of the steel arch in the slightly altered rock section ranged from 0.49 to 0.76 cm, and the average value was 0.62 cm. The deformation of the steel arch in the moderately altered rock section ranged from 0.57 to 1.29 cm, and the average value was 0.97 cm. The deformation of the steel arch in the strongly altered rock section ranged from 1.73 to 2.27 cm, and the average value was 1.93 cm.

After comparison with the numerical simulation results, it was found that the average value of the actual deformation of various altered rocks was basically consistent with the numerical simulation results, and the error was controlled within 1 cm, which validated the numerical simulation. The deformation of some moderately altered rocks was lower than the average value of deformation of slightly altered rocks, which indicates that moderately altered rocks still possessed a certain strength, and the rocks had good stability under supporting measures. The maximum deformation value of strongly altered rocks was 2.27 cm, which was significantly lower compared with unsupported excavation. Hence, the deformation of the surrounding rock was effectively controlled, which verified the reasonableness of the support measures.


**Table 7.** Comparison between the monitored and calculated values of the surrounding rock deformation in tunnel sections with different degrees of alteration.

section of the excavation, and the above-mentioned support method was used to monitor and count the displacement changes of the steel arch after 24 h. The actual monitoring values are presented in Table 7 and Figure 9. Among them, the deformation of the steel arch in the slightly altered rock section ranged from 0.49 to 0.76 cm, and the average value was 0.62 cm. The deformation of the steel arch in the moderately altered rock section ranged from 0.57 to 1.29 cm, and the average value was 0.97 cm. The deformation of the steel arch in the strongly altered rock section ranged from 1.73 to 2.27 cm, and the average

After comparison with the numerical simulation results, it was found that the average value of the actual deformation of various altered rocks was basically consistent with the numerical simulation results, and the error was controlled within 1 cm, which validated the numerical simulation. The deformation of some moderately altered rocks was lower than the average value of deformation of slightly altered rocks, which indicates that moderately altered rocks still possessed a certain strength, and the rocks had good stability under supporting measures. The maximum deformation value of strongly altered rocks was 2.27 cm, which was significantly lower compared with unsupported excavation. Hence, the deformation of the surrounding rock was effectively controlled, which verified

**Figure 9.** Actual monitoring value of tunnel surrounding rock deformation (Asl: calculated deformation value of the slightly altered surrounding rock; Amo: calculated deformation value of the moderately altered surrounding rock; Ast: calculated deformation value of the moderately altered surrounding rock. **Figure 9.** Actual monitoring value of tunnel surrounding rock deformation (Asl: calculated deformation value of the slightly altered surrounding rock; Amo: calculated deformation value of the moderately altered surrounding rock; Ast: calculated deformation value of the moderately altered surrounding rock.

### **5. Conclusions**

value was 1.93 cm.

the reasonableness of the support measures.

This study carried out numerical simulations of a typical tunnel section excavated in the study area under unsupported and supported conditions using the FLAC3D software. The stress field, displacement field, and plastic zone of a typical section of the surrounding rock were selected for calculation and comparison, and the results were analyzed. The following conclusions were drawn:


to the reduction in its stiffness and strength, and thus increased the deformation of the surrounding rock such that destabilization damage occurred. After implementing support measures, the deformation of the tunnel decreased significantly. The actual monitoring value of the surrounding rock displacement was consistent with the simulation results, and the support scheme was reasonable, which provides a theoretical basis for the design and construction of similar projects.

**Author Contributions:** Conceptualization, Y.C.; methodology, J.D.; software, J.L.; validation, T.W., Y.C. and J.D.; formal analysis, S.S.; investigation, X.Y.; resources, Y.C.; data curation, S.S.; writing original draft preparation, S.S.; writing—review and editing, J.D.; visualization, S.S.; supervision, J.D.; project administration, T.W.; funding acquisition, Y.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research work was sponsored by the National Key Research and Development Project of China (Grant No. 2019YFC1509704), the National Natural Science Foundation of China (Grant Nos. U1704243, 41741019, 41977249 and 42090052), Henan Province Science and technology research project (Grant No. 192102310006).

**Data Availability Statement:** The data that supports the findings of this study is available from the corresponding author upon reasonable request.

**Acknowledgments:** We thank all project team members for their contributions to this study.

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

### **References**


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