# Prediction of Final Displacement of Tunnels in Time-Dependent Rock Mass Based on the Nonequidistant Grey Verhulst Model

Authors: Un Chol Han; Chung Song Choe; Kun Ui Hong; Chol Pak
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Time-dependent deformations of rock are important factors for the design and construction of tunnels in the rock masses with time-dependent strength and deformation properties. Therefore, determination displacement and stable time of tunnel are an indispensable task in geotechnical engineering, especially during the design and construction of underground structures. However, in some types of rock masses, the tunnel displacement may increase for months or years after the excavation owing to the rheological behavior of the surrounding rock masses, which greatly influences the selection of the initial tunnel support system, the excavation layout, and the determination of its load-carrying capacity. In this paper, we present a method to quickly predict the final displacement and stable time of tunnels in time-dependent rock mass using the improved nonequidistant grey Verhulst model (INGVM) with Fourier series, which is used to reduce the periodic residual error. The proposed method is validated fairly via a test for a newly excavated roadway of −600 m level in a coal mine, Democratic People’s Republic of Korea. These results are useful references for similar projects in the future.


---

## Body

## 1. Introduction

The stability state of the underground structures such as road tunnels and mine roadways both during the excavation and afterwards can timely be known by monitoring the change of surrounding rock mass displacement. Therefore, the displacement prediction of rock mass surrounding the tunnel can be actually dealt with time series prediction problem.Heretofore, various methods have been proposed by several researchers for time series prediction of the tunnel displacement, such as traditional methods (i.e., analytical solutions or analytical functions) and artificial intelligence-based methods (i.e., artificial neural network (ANN), support vector machine (SVM), and grey models (GMs)), some of which are as follows.Sulem et al. [1] proposed a method based on analytical functions that the closure must be expressed as a function of the two parameters, namely, distance to the face and the time-dependent properties of the ground. They suggested that the closure must be expressed as a function of the above-mentioned parameters, which describes displacements in a plane perpendicular to the tunnel axis as a function of time and the advancing face. Zhao et al. [2] derived an analytical solution to the temporal-spatial displacement of the subsea tunnel lining based on the state equation of the generalized Kelvin constitutive model. Also, Grossauer et al. [3] calculated tunnel displacements from case histories using curve fitting techniques. Kim and Chung [4] suggested a method for predicting final displacement from displacement data measured in the initial stage by using statistical, exponential, and linear functions. Besides, Pan and Dong [5] considered the tunnel convergence as a function of the rheological properties and the stress state of the rock mass. Liu et al. [6] derived closed-form analytical solutions to predict the time-dependent displacement of deep tunnel in swelling rock under different water contents using the fractional viscoplastic model. However, the convergence confinement method, which can be used in circular tunnels in a hydrostatic stress field, is part of the rational approach and uses an analytical type calculation [7, 8]. Unlu and Gercek [9] presented expressions for the normalized radial displacements occurring around tunneling face in the linear elastic medium by using numerical analysis and nonlinear curve fitting, which is depending on Poisson’s ratio. Paraskevopoulou and Diederichs [10] analyzed the total displacements around a circular tunnel in a visco-elastic medium by performing an isotropic axisymmetric ﬁnite difference modeling, based on the convergence confinement method.Recently, along with the development of information technologies, a number of intelligent algorithms, such as regression analysis, grey theory, data smoothing processing, time series analysis, and genetic algorithm, have been applied to predict the displacement of rocks mass surrounding tunnel.Mahdevari and Torabi [11] developed a method based on ANN to identify dependencies between the convergence and the geotechnical and geological conditions encountered. They used a multilayer perceptron analysis to predict convergence during and after tunnel excavation and compared the results with ones of radial basis function and multivariable regression analyses. Xue et al. [12] predicted the nonuniform deformation of soft rock tunnels by using the back-propagation neural network model and Fuzzy Delphi-Rough Set with the deformation data of 30 excavation sections from the study area.Meanwhile, the SVM, as a new machine learning technology, is one of the methods of tackling this problem, which is a supervised machine learning method based on the statistical learning theory [13]. Yao et al. [14] and Wu et al. [15] predicted displacement of rock mass surrounding tunnel using SVM and ANN. Besides, Yao et al. [16] presented a hybrid prediction model based on SVM for predicting the tunnel surrounding rock displacement and genetic algorithm (GA) for optimizing the parameters in the SVM. Xue and Xiao [17] proposed a least-squares support vector machine (LSSVM) method based on particle swarm optimization (PSO) algorithm to predict the deformation of rock mass surrounding underground caverns in a hydropower station during excavation. From the above-mentioned methods, the researchers have noticed that SVM has a high-accuracy prediction than ANN or finite element method (FEM).Huang et al. [18] used time-series analysis termed as autoregressive moving average (ARMA (n, m)) model and grey model (GM(1,1)) to predict rock mass displacement at the key measuring points of the permanent shiplock in the Three Gorges Project, and Guo et al. [19] presented the nonequidistant grey model (NGM(1,1)) reduced origin error to predict final tunnel surrounding rock displacement based on the data of early 20 days. Han et al. [20] predicted the final displacement of tunnel in time-dependent rock mass by using the traditional nonequidistant grey Verhulst model (NGVM). Also, the Gaussian process (GP), SVM, the wavelet neural network (WNN), and GM(1,1) are analyzed comparatively for predicting the surrounding rock nonlinear deformation [21]. The conclusion from their findings states that the GP performs best, and the SVM performs better than the WNN and GM(1,1) methods.As stated in former studies, a traditional statistical method such as regression analysis needs a large amount of monitoring data and is confined to use the differential equations to build discrete stochastic model, and it is not easy to represent the nature and the inherent laws of the change progress. Also, the ANN model for time series prediction demands a great deal of training data. Furthermore, it is complex to obtain optimal model parameters (e.g., kernel function) through an effective method, because SVM is sensitive to the selection of the parameters. However, it is widely used for displacement prediction in geotechnical engineering practices, because GMs require only a limited amount of data, and has simple calculation process and higher forecasting accuracy. Nevertheless, the common grey prediction model only applies to strictly equal interval sequences, adding to the difficulty of engineering application. Especially, GM(1,1) is imperfect when the time series increases in the curve with S-type, or the increment of time series is in the saturation stage.Therefore, the present study in this paper attempts to use the improved nonequidistant grey Verhulst model (INGVM) with Fourier series to predict the displacement of rock mass surrounding tunnel, which is gradually stabilized in the shape of S-type saturated curve. The objectives of this study are to examine the feasibility of applying INGVM compared to traditional NGM(1,1) and NGVM, in order to predict the final displacement and stable time of tunnel excavated in time-dependent rock mass, with only monitored displacements during its excavation. The comparison shows that the prediction by the proposed model has high accuracy and strong adaptability.


## 2. Deformation Characteristics of Rock Mass Surrounding Tunnel

### 2.1. Time-Dependent Behavior of Surrounding Rock Mass during and after Tunnel Excavation

Time-dependent deformation of rock plays an important role in the design and construction of tunnels in the rock masses with time-dependent strength and deformation characteristics. Particularly, it is common in tunnels excavated in soft rock, heavily sheared weak rock masses, or rock masses with high in situ stress. It is generally believed that hard rock does not have obvious creep phenomena, but some studies show that hard rocks surrounding deep tunnel can exhibit creep behavior under high-stress conditions [22, 23]. In tunneling, creep behavior emerges as the ongoing increase of the radial displacements observed in the tunnel surface [9].The longitudinal deformation profile (LDP) or the deformation history for a point located on tunnel wall is similar to a sigmoid (S-type) curve [7, 9, 24]. Especially, for the circular tunnel in a hydrostatic in situ stress field, the radial displacement of a point located on the periphery ahead of the face occurs from distance approximately equal to 2 times of tunnel diameter as shown in Figure 1. Behind the face, the displacement increases rapidly and, finally, converges to the plane strain value (ur∞) at some distance, which is approximately 2–4 times of tunnel diameter for a linear elastic rock while it can be several diameters for elastoplastic rock [7, 25].Figure 1
Longitudinal deformation profile (LDP) of a point located on the surface of a circular tunnel excavated in a hydrostatic stress field [9].From the results of earlier research mentioned above, it is found that displacements of surrounding rock mass during tunnel excavation and operation are greatly influenced by face advance and duration time of the tunnel. Initially, when a tunnel is driven through a rock mass, which exhibits time-dependent behavior such as sandy shale and siltstone, tunnel displacements increase sharply at the early stage because of the disturbance in the in-situ stress state. Thereafter, tunnel displacement is gradually stabilized at a constant deformation rate in the shape of S-type saturated curve, which reflects the time-dependence of ground stress under the condition without additional disturbance in rock mass. Therefore, in this study, we choose the NGVM suitable to predict the displacement of nonequidistant time series as S-type curve, in order to predict the final displacement of rock mass surrounding the unsupported mine roadways.


### 2.2. Criteria for Stable Displacement of Rock Mass Surrounding Tunnel

In rock mechanics and rock engineering practice, criteria for stable displacement are of great importance for displacement back analysis of geotechnical parameters of rock mass, as well as stability analysis of underground structures.The total displacement (mm) and maximum stable deformation rate (mm/d) are often used to ensure the stability of tunnel during and after their construction. In many cases, the deformation rate is used as an important index to reflect the overall deformation and the stability of rock mass. For example, Jiang et al. [26] identified the stable regions in the three pilot tunnels, based on the deformation rate (less than 0.02 mm/d). Also, Yu et al. [27] proposed that if the deformation rate is less than 1.5 mm/d, the roadway tends to stabilize, based on the convergences measured in several roadways of deep mines under the conditions of high stress. Meanwhile, Yao et al. [16] and Wu et al. [28] have recognized that surrounding rock mass is in a steady state when the deformation rate of tunnel excavated in rock mass reached any specific value (generally less than 0.1 mm/d) at the monitoring stage.Based on practical experiences and in situ monitoring at northern coal fields of our country including Coal Mine “6.13” (see Section4.1 for a brief introduction about this mine), it could be found that rock mass surrounding the main roadways in different deep levels (more than −500 m below ground surface) was in a steady state with a gradual decrease in deformation rate, usually over a period of 6 months to 2 years after excavation, although deformation magnitude increased to some extent by various mining activities such as stope blasting and coal caving mining nearby the main roadways. As a result of the aforementioned considerations, we selected 0.05 mm/d as the maximum stable deformation rate to predict the final displacement of tunnel using grey prediction models, based on stability and deformation data measured in various roadways previously constructed in the upper and same levels of Coal Mine “6.13”.


## 2.1. Time-Dependent Behavior of Surrounding Rock Mass during and after Tunnel Excavation

Time-dependent deformation of rock plays an important role in the design and construction of tunnels in the rock masses with time-dependent strength and deformation characteristics. Particularly, it is common in tunnels excavated in soft rock, heavily sheared weak rock masses, or rock masses with high in situ stress. It is generally believed that hard rock does not have obvious creep phenomena, but some studies show that hard rocks surrounding deep tunnel can exhibit creep behavior under high-stress conditions [22, 23]. In tunneling, creep behavior emerges as the ongoing increase of the radial displacements observed in the tunnel surface [9].The longitudinal deformation profile (LDP) or the deformation history for a point located on tunnel wall is similar to a sigmoid (S-type) curve [7, 9, 24]. Especially, for the circular tunnel in a hydrostatic in situ stress field, the radial displacement of a point located on the periphery ahead of the face occurs from distance approximately equal to 2 times of tunnel diameter as shown in Figure 1. Behind the face, the displacement increases rapidly and, finally, converges to the plane strain value (ur∞) at some distance, which is approximately 2–4 times of tunnel diameter for a linear elastic rock while it can be several diameters for elastoplastic rock [7, 25].Figure 1
Longitudinal deformation profile (LDP) of a point located on the surface of a circular tunnel excavated in a hydrostatic stress field [9].From the results of earlier research mentioned above, it is found that displacements of surrounding rock mass during tunnel excavation and operation are greatly influenced by face advance and duration time of the tunnel. Initially, when a tunnel is driven through a rock mass, which exhibits time-dependent behavior such as sandy shale and siltstone, tunnel displacements increase sharply at the early stage because of the disturbance in the in-situ stress state. Thereafter, tunnel displacement is gradually stabilized at a constant deformation rate in the shape of S-type saturated curve, which reflects the time-dependence of ground stress under the condition without additional disturbance in rock mass. Therefore, in this study, we choose the NGVM suitable to predict the displacement of nonequidistant time series as S-type curve, in order to predict the final displacement of rock mass surrounding the unsupported mine roadways.


## 2.2. Criteria for Stable Displacement of Rock Mass Surrounding Tunnel

In rock mechanics and rock engineering practice, criteria for stable displacement are of great importance for displacement back analysis of geotechnical parameters of rock mass, as well as stability analysis of underground structures.The total displacement (mm) and maximum stable deformation rate (mm/d) are often used to ensure the stability of tunnel during and after their construction. In many cases, the deformation rate is used as an important index to reflect the overall deformation and the stability of rock mass. For example, Jiang et al. [26] identified the stable regions in the three pilot tunnels, based on the deformation rate (less than 0.02 mm/d). Also, Yu et al. [27] proposed that if the deformation rate is less than 1.5 mm/d, the roadway tends to stabilize, based on the convergences measured in several roadways of deep mines under the conditions of high stress. Meanwhile, Yao et al. [16] and Wu et al. [28] have recognized that surrounding rock mass is in a steady state when the deformation rate of tunnel excavated in rock mass reached any specific value (generally less than 0.1 mm/d) at the monitoring stage.Based on practical experiences and in situ monitoring at northern coal fields of our country including Coal Mine “6.13” (see Section4.1 for a brief introduction about this mine), it could be found that rock mass surrounding the main roadways in different deep levels (more than −500 m below ground surface) was in a steady state with a gradual decrease in deformation rate, usually over a period of 6 months to 2 years after excavation, although deformation magnitude increased to some extent by various mining activities such as stope blasting and coal caving mining nearby the main roadways. As a result of the aforementioned considerations, we selected 0.05 mm/d as the maximum stable deformation rate to predict the final displacement of tunnel using grey prediction models, based on stability and deformation data measured in various roadways previously constructed in the upper and same levels of Coal Mine “6.13”.


## 3. Methodology

Grey system theory is an interdisciplinary scientific area that was first introduced by Deng [29] in the early 1980s to study a system with incomplete information. Thereafter, the grey system has been broadly applied to many aspects of data processing. The displacement data monitored in situ can be considered as a grey system because they have characteristics of randomness and fuzziness. However, for most cases of tunneling engineering, tunnel displacement or convergence data obtained in situ is generally time sequence of unequal interval. In this paper, hence, we have constructed the nonequidistant grey system to predict data in an unequal interval sequence. In the following sections, the procedures in the construction of the nonequidistant grey models (NGMs, i.e., NGM(1,1) and NGVM) for the prediction of tunnel displacement are described.


### 3.1. Nonequidistant Grey Model (NGM(1,1))

Definition 1.
(see [30]). Suppose that X(0) = {x(0) (t1), x(0) (t2),…, x(0) (tn)} is a nonnegative sequence of original data, where △ti = ti − ti−1≠const, i = 2, 3, …, n, and then, X(0) is called a nonequidistant sequence. The first-order accumulated generation operation (1-AGO) sequence generated from nonequidistant sequence X(0) is defined as(1)X1=x1t1,x1t2,…,x1tn,where n is the number of data in the sequence, and(2)x1tk=x0t1,k=1,∑i=1kx0tiΔti,k=2,3,…,n..
The background value sequence generated from consecutive neighbors of the 1-AGO sequence is defined as(3)Z1=z1t2,z1t3,…,z1tn,where(4)z1tk=λx1tk+1−λx1tk−1,k=2,3,…,n.
In equation (4), λ is the adjustment coefficient of the background value sequence in the range of [0, 1], which is usually set to 0.5.Definition 2.
LetX(0), X(1), and Z(1) be the same as in Definition 1; then, the equation(5)x0tk+az1tk=b,is called the basic form of the NGM(1, 1), where a is the developmental coefficient and b is the grey action.
And, the differential equation(6)dx1tkdtk+ax1tk=b,is called the whitening equation of the NGM(1,1).
The parameter vector Φ = [a, b]T of NGM(1,1) can be obtained by the least-square method, namely,(7)Φ=a,bT=BTB−1BTY,where(8)B=−z1t21−z1t31⋮⋮−z1tn1,Y=x0t2x0t3⋮x0tn.
By substituting the initial condition into the solution of the whitening (6), the discrete-time response sequence of the grey model is formulated as(9)x^1tk=x0t1−b^a^e−a^tk−t1+b^a^,k=2,3,…,n.
Next, by using the first-order inverse accumulative generating operation (1-IAGO), the restored value of the grey model is calculated as(10)x^0tk=x^1tk−x^1tk−1Δtk,where x^0t1=x0t1.
Then, the sequence produced by the NGM(1,1) is represented as(11)x^0tk=1Δtkx0t1−b^a^1−ea^Δtke−a^tk1−t1,k=2,3,…,n.
By analyzing the above-stored sequence (11), it can be noticed that the NGM(1,1) has an ideal performance for the unequal interval time series with exponential characteristics. However, it is imperfect when the displacement increases in the curve with S-type or its increment is in the saturation stage. In this case, the predicting error of the grey system model will become larger, and the result is unaccepted in engineering practices.


### 3.2. Nonequidistant Grey Verhulst Model (NGVM)

GVM is a special kind of model within the grey system, which is first introduced by German biologist Pierre Franois Verhulst [31, 32]. The main purpose of GVM is to limit the whole development for a real system, and it is effective in describing some increasing processes, such as an S-curve that has a saturation region [20]. The NGVM has the same characteristics as the NGM(1,1). It primarily adds a restriction to the NGM(1,1).Definition 3.
(see [31]). Let X(0), X(1), and Z(1) be the same as in Definition 1; then,(12)dx1tkdtk+ax1tk=bx1tk2,is called the whitening equation of the NGVM.
Similar to NGM(1,1),(13)Φ=a,bT=BTB−1BTY,where(14)B=−z1t2z1t22−z1t3z1t32⋮⋮−z1tnz1tn2,Y=x0t2x0t3⋮x0tn.
The time response sequence of NGVM is defined as(15)x^1tk=ax0t1bx0t1+a−bx0t1eatk−t1,k=1,2,…,n.
Applying the IAGO, the sequence produced by the NGVM is obtained as(16)x^0tk=1Δtkax0t1a−bx0t1bx0t1+a−bx0t1eatk−t1∗1−eaΔtkeatk−1−t1bx0t1+a−bx0t1eatk−1−t1,k=2,3,…,n.
From the solution of the NGVM, it can be seen that whent⟶∞; if a > 0, then x(1) (tk)⟶0, and if a < 0, then x(1) (tk)⟶a/b [31]. That is, when t is sufficiently large, x(1) (tk) and x(1) (tk-1) will be sufficiently close. Therefore, the NGVM is applicable to describe a monotonously increasing curve with S-type.
Equations (11) and (16) can be used to make predictions for nonequidistant data sequence, such as displacement and convergence of rock mass surrounding tunnel.


### 3.3. Modification of NGVM Using Fourier Series of Error Residuals

The traditional NGMs (NGM(1,1) and NGVM) contain the model parameters (a and b). The solution curves of the NGM(1,1) and NGVM shown in equations (11) and (16) are influenced by the values of these parameters. Especially, NGVM is more suitable for S-type time series with saturation state trend. However, in the practical application process, the fitting situation of the original data always deviates from the ideal S-type, thus reducing the prediction accuracy of the model. At the same time, because of the internal of initial value x(1) (t1) and adjustment coefficient λ of the background value sequence for NGMs, the prediction error of the model will also increase. Therefore, in this study, Fourier series is introduced to reduce the periodic residual error in NGVM [31–33].Definition 4.
(see [32]). Consider the X(0) sequence in Definition 1 and the predicted values given by the NGMs, and then, the error sequence of X(0) can be defined as(17)ε0=ε0t2,ε0t3,…,ε0tn,where(18)ε0tk=x0tk−x^0tk,k=2,3,…,n.
The error residuals in equation (21) can be expressed in Fourier series as follows:(19)ε0tk≅12a0+∑i=1Faicos2πiTk+bisin2πiTk,k=2,3,…,n.
In equation (19), T=n−1 indicates the length of the residual series, and F=n−1/2−1 means the minimum deployment frequency of Fourier series. It is obvious that T will be an integer number, and F will be selected as an integer number [32].
Equation (22) can be rewritten as follows:(20)ε0≅PC.
P and C matrixes can be defined as follows:(21)P=1/2cos2πT2sin2πT2cos2π2T2sin2π2T2⋯cos2πFT2sin2πFT21/2cos2πT3sin2πT3cos2π2T3sin2π2T3⋯cos2πFT3sin2πFT3⋮⋮⋮⋮⋮⋮⋮⋮1/2cos2πTnsin2πTncos2π2Tnsin2π2Tn⋯cos2πFTnsin2πFTn,C=a0a1b1a2b2⋯anbnT.
One can use least-squares method to solve equation (20) and to calculate the matrix C:(22)C≅PTP−1PTε0.
Fourier series correction can be obtained as follows:(23)x^f0tk=x^0tk+ε^0tk,k=2,3,…,n+1.
By modifying the residual values of the traditional NGMs with Fourier series, the fluctuation of the grey prediction models can be reduced, and the adaptability and stability of the prediction model can be also improved.


### 3.4. Modeling Procedure

The detailed procedure of improved NGVM for predicting the final displacement of tunnel excavated in a time-dependent rock mass can be described as follows.Step 1.
Displacement data measured in situ during tunnel excavation and afterwards are collected, and the maximum stable deformation rate for rock mass surrounding tunnel is defined.Step 2.
Original sequenceX(0) and 1-AGO sequence X(1) are established, and then, background value sequence Z(1) is generated according to Definition 1. If measured displacement data are given as accumulative values of displacement and (or) convergence, the original sequence X(0) is obtained from 1-IAGO of in situ measured data.Step 3.
According to Definition3, the grey differential equations of NGVM are constructed, and then, the matrices B and Y are constructed, and the parameters a and b are estimated by the least-square method.Step 4.
The time response sequenceX^1 is calculated, and then, according to 1-IAGO, the predicted values X^0 are calculated using equations (11) and (16).Step 5.
The sequence of residual valuesε0tk is generated, and subsequently, the newly predicted residuals ε^0tk are generated with Fourier series, according to Definition 4.Step 6.
New prediction valuesX^f0 are generated using equation (23).Step 7.
The final displacementx^f0tk and stable time tk+1 of tunnel are predicted when all the difference between adjacent predicted values corresponding to the measurement components of tunnel displacement satisfied the stable deformation rate, namely, x^f0tk−x^f0tk−1<0.05 mm/d.
The flowchart for predicting tunnel displacement in time-dependent rock mass by using stable deformation rate and NGVM with Fourier series is shown in Figure2. The proposed NGVM and the contrastive grey models for case study are all realized by MATLAB platform 2019a.Figure 2
The flowchart for displacement prediction of the tunnel in time-dependent rock mass by INGVM.


### 3.5. Evaluation Criteria of the Predicting Performance

Prediction accuracy is an important criterion for evaluating predictive validity. In order to compare performances and reliability of proposed grey models, the two statistical indicators, that is, the mean absolute percentage error (MAPE) and the root mean square percentage error (RMSPE), are collected for this study. The MAPE and the RMSPE are deﬁned as [33, 34](24)MAPE=1n−1∑k=2nx^1k−x1kx1k×100%,RMSPE=∑k=2nx^1k−x1k/x1k2n−1×100%,where x1k is the actual value at time k, x^1k is its model value, and n is the number of data used for prediction.The values of MAPE and RMSPE sufficiently describe the goodness of the prediction effect. The lower the values of MAPE and RMSPE, the more accurate the prediction. In addition, the criteria of MAPE and RMSPE for prediction error are shown in Table1 [34, 35].Table 1
The MAPE and RMSPE criteria of prediction precision.
MAPE and RMSPE (%)<1010–2020–50>50Prediction precisionExcellentGoodReasonableUnacceptable


## 3.1. Nonequidistant Grey Model (NGM(1,1))

Definition 1.
(see [30]). Suppose that X(0) = {x(0) (t1), x(0) (t2),…, x(0) (tn)} is a nonnegative sequence of original data, where △ti = ti − ti−1≠const, i = 2, 3, …, n, and then, X(0) is called a nonequidistant sequence. The first-order accumulated generation operation (1-AGO) sequence generated from nonequidistant sequence X(0) is defined as(1)X1=x1t1,x1t2,…,x1tn,where n is the number of data in the sequence, and(2)x1tk=x0t1,k=1,∑i=1kx0tiΔti,k=2,3,…,n..
The background value sequence generated from consecutive neighbors of the 1-AGO sequence is defined as(3)Z1=z1t2,z1t3,…,z1tn,where(4)z1tk=λx1tk+1−λx1tk−1,k=2,3,…,n.
In equation (4), λ is the adjustment coefficient of the background value sequence in the range of [0, 1], which is usually set to 0.5.Definition 2.
LetX(0), X(1), and Z(1) be the same as in Definition 1; then, the equation(5)x0tk+az1tk=b,is called the basic form of the NGM(1, 1), where a is the developmental coefficient and b is the grey action.
And, the differential equation(6)dx1tkdtk+ax1tk=b,is called the whitening equation of the NGM(1,1).
The parameter vector Φ = [a, b]T of NGM(1,1) can be obtained by the least-square method, namely,(7)Φ=a,bT=BTB−1BTY,where(8)B=−z1t21−z1t31⋮⋮−z1tn1,Y=x0t2x0t3⋮x0tn.
By substituting the initial condition into the solution of the whitening (6), the discrete-time response sequence of the grey model is formulated as(9)x^1tk=x0t1−b^a^e−a^tk−t1+b^a^,k=2,3,…,n.
Next, by using the first-order inverse accumulative generating operation (1-IAGO), the restored value of the grey model is calculated as(10)x^0tk=x^1tk−x^1tk−1Δtk,where x^0t1=x0t1.
Then, the sequence produced by the NGM(1,1) is represented as(11)x^0tk=1Δtkx0t1−b^a^1−ea^Δtke−a^tk1−t1,k=2,3,…,n.
By analyzing the above-stored sequence (11), it can be noticed that the NGM(1,1) has an ideal performance for the unequal interval time series with exponential characteristics. However, it is imperfect when the displacement increases in the curve with S-type or its increment is in the saturation stage. In this case, the predicting error of the grey system model will become larger, and the result is unaccepted in engineering practices.


## 3.2. Nonequidistant Grey Verhulst Model (NGVM)

GVM is a special kind of model within the grey system, which is first introduced by German biologist Pierre Franois Verhulst [31, 32]. The main purpose of GVM is to limit the whole development for a real system, and it is effective in describing some increasing processes, such as an S-curve that has a saturation region [20]. The NGVM has the same characteristics as the NGM(1,1). It primarily adds a restriction to the NGM(1,1).Definition 3.
(see [31]). Let X(0), X(1), and Z(1) be the same as in Definition 1; then,(12)dx1tkdtk+ax1tk=bx1tk2,is called the whitening equation of the NGVM.
Similar to NGM(1,1),(13)Φ=a,bT=BTB−1BTY,where(14)B=−z1t2z1t22−z1t3z1t32⋮⋮−z1tnz1tn2,Y=x0t2x0t3⋮x0tn.
The time response sequence of NGVM is defined as(15)x^1tk=ax0t1bx0t1+a−bx0t1eatk−t1,k=1,2,…,n.
Applying the IAGO, the sequence produced by the NGVM is obtained as(16)x^0tk=1Δtkax0t1a−bx0t1bx0t1+a−bx0t1eatk−t1∗1−eaΔtkeatk−1−t1bx0t1+a−bx0t1eatk−1−t1,k=2,3,…,n.
From the solution of the NGVM, it can be seen that whent⟶∞; if a > 0, then x(1) (tk)⟶0, and if a < 0, then x(1) (tk)⟶a/b [31]. That is, when t is sufficiently large, x(1) (tk) and x(1) (tk-1) will be sufficiently close. Therefore, the NGVM is applicable to describe a monotonously increasing curve with S-type.
Equations (11) and (16) can be used to make predictions for nonequidistant data sequence, such as displacement and convergence of rock mass surrounding tunnel.


## 3.3. Modification of NGVM Using Fourier Series of Error Residuals

The traditional NGMs (NGM(1,1) and NGVM) contain the model parameters (a and b). The solution curves of the NGM(1,1) and NGVM shown in equations (11) and (16) are influenced by the values of these parameters. Especially, NGVM is more suitable for S-type time series with saturation state trend. However, in the practical application process, the fitting situation of the original data always deviates from the ideal S-type, thus reducing the prediction accuracy of the model. At the same time, because of the internal of initial value x(1) (t1) and adjustment coefficient λ of the background value sequence for NGMs, the prediction error of the model will also increase. Therefore, in this study, Fourier series is introduced to reduce the periodic residual error in NGVM [31–33].Definition 4.
(see [32]). Consider the X(0) sequence in Definition 1 and the predicted values given by the NGMs, and then, the error sequence of X(0) can be defined as(17)ε0=ε0t2,ε0t3,…,ε0tn,where(18)ε0tk=x0tk−x^0tk,k=2,3,…,n.
The error residuals in equation (21) can be expressed in Fourier series as follows:(19)ε0tk≅12a0+∑i=1Faicos2πiTk+bisin2πiTk,k=2,3,…,n.
In equation (19), T=n−1 indicates the length of the residual series, and F=n−1/2−1 means the minimum deployment frequency of Fourier series. It is obvious that T will be an integer number, and F will be selected as an integer number [32].
Equation (22) can be rewritten as follows:(20)ε0≅PC.
P and C matrixes can be defined as follows:(21)P=1/2cos2πT2sin2πT2cos2π2T2sin2π2T2⋯cos2πFT2sin2πFT21/2cos2πT3sin2πT3cos2π2T3sin2π2T3⋯cos2πFT3sin2πFT3⋮⋮⋮⋮⋮⋮⋮⋮1/2cos2πTnsin2πTncos2π2Tnsin2π2Tn⋯cos2πFTnsin2πFTn,C=a0a1b1a2b2⋯anbnT.
One can use least-squares method to solve equation (20) and to calculate the matrix C:(22)C≅PTP−1PTε0.
Fourier series correction can be obtained as follows:(23)x^f0tk=x^0tk+ε^0tk,k=2,3,…,n+1.
By modifying the residual values of the traditional NGMs with Fourier series, the fluctuation of the grey prediction models can be reduced, and the adaptability and stability of the prediction model can be also improved.


## 3.4. Modeling Procedure

The detailed procedure of improved NGVM for predicting the final displacement of tunnel excavated in a time-dependent rock mass can be described as follows.Step 1.
Displacement data measured in situ during tunnel excavation and afterwards are collected, and the maximum stable deformation rate for rock mass surrounding tunnel is defined.Step 2.
Original sequenceX(0) and 1-AGO sequence X(1) are established, and then, background value sequence Z(1) is generated according to Definition 1. If measured displacement data are given as accumulative values of displacement and (or) convergence, the original sequence X(0) is obtained from 1-IAGO of in situ measured data.Step 3.
According to Definition3, the grey differential equations of NGVM are constructed, and then, the matrices B and Y are constructed, and the parameters a and b are estimated by the least-square method.Step 4.
The time response sequenceX^1 is calculated, and then, according to 1-IAGO, the predicted values X^0 are calculated using equations (11) and (16).Step 5.
The sequence of residual valuesε0tk is generated, and subsequently, the newly predicted residuals ε^0tk are generated with Fourier series, according to Definition 4.Step 6.
New prediction valuesX^f0 are generated using equation (23).Step 7.
The final displacementx^f0tk and stable time tk+1 of tunnel are predicted when all the difference between adjacent predicted values corresponding to the measurement components of tunnel displacement satisfied the stable deformation rate, namely, x^f0tk−x^f0tk−1<0.05 mm/d.
The flowchart for predicting tunnel displacement in time-dependent rock mass by using stable deformation rate and NGVM with Fourier series is shown in Figure2. The proposed NGVM and the contrastive grey models for case study are all realized by MATLAB platform 2019a.Figure 2
The flowchart for displacement prediction of the tunnel in time-dependent rock mass by INGVM.


## 3.5. Evaluation Criteria of the Predicting Performance

Prediction accuracy is an important criterion for evaluating predictive validity. In order to compare performances and reliability of proposed grey models, the two statistical indicators, that is, the mean absolute percentage error (MAPE) and the root mean square percentage error (RMSPE), are collected for this study. The MAPE and the RMSPE are deﬁned as [33, 34](24)MAPE=1n−1∑k=2nx^1k−x1kx1k×100%,RMSPE=∑k=2nx^1k−x1k/x1k2n−1×100%,where x1k is the actual value at time k, x^1k is its model value, and n is the number of data used for prediction.The values of MAPE and RMSPE sufficiently describe the goodness of the prediction effect. The lower the values of MAPE and RMSPE, the more accurate the prediction. In addition, the criteria of MAPE and RMSPE for prediction error are shown in Table1 [34, 35].Table 1
The MAPE and RMSPE criteria of prediction precision.
MAPE and RMSPE (%)<1010–2020–50>50Prediction precisionExcellentGoodReasonableUnacceptable


## 4. Case Study

### 4.1. In Situ Measurements

The Coal Mine “6.13” in North Hamgyong Province of the Democratic People’s Republic of Korea is a bituminous coal mine with a designed annual production capacity of a hundred thousand tonnes. For developing the new №27 block in the coal mine, a main roadway of 80 m length with 3.6 m width and 3.0 m height is driven at a depth of about 650 m below ground surface, at a constant advance rate equal to 2-3 m/d with the full face using drill-and-blasting excavation method in sandy siltstone, which assumes relatively a continuous, homogeneous, and isotropic rock mass as a floor of 6–1# coal seam, as shown in Figure3(a) [20]. Dimensions of the main roadway and the layout of measuring points in a typical section are illustrated in Figure 3(b).Figure 3
(a) 3D model for development state of coal deposit and (b) layout of measuring points in a typical section of the main roadway.In order to monitor the deformation behavior of the rock mass, it is essential to set up measuring stations. Although there are many measuring stations in order to obtain the full displacement vector mostly for study in tunneling, used by several researchers in the past, 1–3 monitoring lines or 3–5 stations are required in practical engineering [20, 21, 36]. In this study, the two typical sections (i.e., Section I and II) for measuring the displacement of surrounding rock mass were, respectively, established at 17 m and 26 m away from start point of the roadway excavation along its axis. Also, the measuring stations are, respectively, installed at the crown (A) and both sidewalls (B and C) on the roadway wall in each monitoring section, resulting in that the measurement components for deformation of the roadway surface can be selected as the crown settlement (zA) and the sidewalls convergence (cBC), as shown in Figure 3(b).All measuring stations were arranged at a distance approximately 1 m from the working face, as soon as the face advanced away from the monitoring section. The roadway excavation started on 22 October 2017 and was finished in a period of 30 days. The displacement measurements with the total station instrument (model: DTM-352L) were performed once in a few days during the roadway excavation and once in 10–15 days after the excavation until all stations are fully stabilized with a stable deformation rate of less than 0.05 mm/d. The in situ monitored values for two monitoring sections are presented in Figure4.Figure 4
Displacement curves of roadway surface measured in each monitoring section.As indicated in Figure4, it could be found that rock mass surrounding roadway represents significant time-dependent deformation both during the excavation and afterwards. Moreover, the crown settlement and the sidewalls convergence at the end of the roadway excavation are, respectively, 14.2 mm and 37.6 mm in Section I and 12.4 mm and 34.2 mm in Section II. Afterwards, the surrounding rock mass was in a steady state with a stable deformation rate after approximately 8 months (230 days) from the start of roadway excavation. In detail, stable times of surrounding rock mass in Section I and Section II were surveyed as 210 days and 225 days, respectively. At the same time, the final values of the crown settlement and the sidewalls convergence were 22.2 mm and 65.8 mm in Section I and 21.8 mm and 63.3 mm in Section II, respectively.


### 4.2. Results and Discussion

It can be found from Figure4 that measured displacement data are given as accumulative values of the crown settlement and the sidewalls convergence. Therefore, we established the original displacement sequence X(0) by applying 1-IAGO on the displacement data measured in a period of the roadway excavation and constructed three NGMs (i.e., INGVM, NGM(1,1), and NGVM) for comparison, respectively. Then, the values of the parameters for different NGMs corresponding to measurement components of tunnel displacement in each monitoring section were estimated by the least-square method, which is indicated in Table 2.Table 2
The estimated values of parameters for NGMs by the least-square method.
Grey modelModel parameterSection ISection IIzAcBCzAcBCNGM(1, 1)A0.00850.00510.00670.0058B1.27433.21291.36744.0326NGVM and INGVMA−0.1505−0.1701−0.1096−0.1132B−2.8344e-3−1.1378e-03−1.0417e-3−3.0344e-4By using the estimated values of parameters in Table2, the fitting results of the prediction curves by the three compared models with the monitoring curves of roadway displacement at two monitoring sections are shown in Figure 5, and the values of two accuracy test criteria (MAPE and RMSPE) on fitting and prediction results of different prediction models are all listed in Table 3.Figure 5
Comparison curves between in situ monitored data and predicted results by three compared models. (a) Section I. (b) Section II.
(a)(b)Table 3
Accuracy check of fitting and prediction results for different prediction models.
Monitoring componentNGM(1,1)NGVMINGVMMAPERMSPEMAPERMSPEMAPERMSPEFittingSection IzA6.148.6515.5716.725.056.27cBC7.657.1018.0319.856.046.67Section IIzA7.418.1116.5918.025.377.67cBC8.7710.0318.8319.747.048.34Mean7.498.4717.4617.265.887.24PredictionSection IzA7.068.3514.2315.963.694.16cBC8.6510.1914.4015.673.644.13Section IIzA7.3510.5915.7217.154.766.09cBC8.1911.8013.5916.514.365.71Mean7.8110.2314.4916.324.115.02As shown in Figure5 and Table 3, it can be found that the MAPE and RMSPE between in situ measured data and predicted results by three compared models are all within 20%. Moreover, it can be seen from Figure 5 that, in prediction interval, INGVM and NGVM simulate a similar trend (i.e., S-type curve) of actual data, but the simulation of NGM(1,1) has an ascending trend. This may be the reason why the NGVM was applied in this study.From the viewpoint of error testing, it is obvious that the prediction accuracy of the INGVM is higher than that of the traditional NGMs in displacement prediction of the time-dependent rock mass. Especially, MAPE and RMSPE of INGVM were 5.88% and 7.24% in the fitting interval and 4.11% and 5.02% in prediction interval, respectively, indicating that prediction errors of INGVM are all smaller than those of the NGM(1,1) and NGVM. In addition, the mean values of MAPE and RMSPE of INGVM in the prediction interval are reduced from 14.49% to 16.32% to 4.11% and 5.02%, meaning the INGVM is more effective and applicable than the traditional NGVM. It is shown that Fourier series helped INGVM to further reduce the periodic residual error compared with traditional NGMs. Meanwhile, we can also see from Figure5 and Table 3 that the errors in fitting interval with relation to NGVM and INGVM are slightly lower than those in prediction interval because the fitting situation of displacement data measured in situ to be used at fitting interval deviates from the ideal S-type.At the same time, it could be also predicted that stable time, final crown settlement, and sidewalls convergence predicted by INGVM are 218 days, 22.43 mm and 66.37 mm in Section I, and 235 days, 22.01 mm and 62.83 mm in Section II, respectively, indicating that prediction results by the proposed model with criterion for the stable displacement of 0.05 mm/d are in good agreement with in situ monitored data. It can be seen that prediction results’ stable time and final displacement by using the proposed model tend to slightly increase than real values.From the above comparative analysis, we can conclude that the proposed INGVM method is a better fitting model compared with NGM(1,1) and NGVM to predict the final displacement and stable time of time-dependent rock mass surrounding roadways as a moderate accuracy, using only monitored in situ data during excavation. Unfortunately, it may not be suitable for the cases when the deformation of rock mass surrounding unsupported and supported underground structures is not just unstable but either continuously increased according to time changes or suddenly increased. However, if displacement histories of soil tunnel and highly jointed rock tunnel show the tendency gradually to be stabilized in S-type curve by rational excavation method, suitable support system, and ground reinforcement, proposed model may be also used for displacement prediction of such rock tunnels.


## 4.1. In Situ Measurements

The Coal Mine “6.13” in North Hamgyong Province of the Democratic People’s Republic of Korea is a bituminous coal mine with a designed annual production capacity of a hundred thousand tonnes. For developing the new №27 block in the coal mine, a main roadway of 80 m length with 3.6 m width and 3.0 m height is driven at a depth of about 650 m below ground surface, at a constant advance rate equal to 2-3 m/d with the full face using drill-and-blasting excavation method in sandy siltstone, which assumes relatively a continuous, homogeneous, and isotropic rock mass as a floor of 6–1# coal seam, as shown in Figure3(a) [20]. Dimensions of the main roadway and the layout of measuring points in a typical section are illustrated in Figure 3(b).Figure 3
(a) 3D model for development state of coal deposit and (b) layout of measuring points in a typical section of the main roadway.In order to monitor the deformation behavior of the rock mass, it is essential to set up measuring stations. Although there are many measuring stations in order to obtain the full displacement vector mostly for study in tunneling, used by several researchers in the past, 1–3 monitoring lines or 3–5 stations are required in practical engineering [20, 21, 36]. In this study, the two typical sections (i.e., Section I and II) for measuring the displacement of surrounding rock mass were, respectively, established at 17 m and 26 m away from start point of the roadway excavation along its axis. Also, the measuring stations are, respectively, installed at the crown (A) and both sidewalls (B and C) on the roadway wall in each monitoring section, resulting in that the measurement components for deformation of the roadway surface can be selected as the crown settlement (zA) and the sidewalls convergence (cBC), as shown in Figure 3(b).All measuring stations were arranged at a distance approximately 1 m from the working face, as soon as the face advanced away from the monitoring section. The roadway excavation started on 22 October 2017 and was finished in a period of 30 days. The displacement measurements with the total station instrument (model: DTM-352L) were performed once in a few days during the roadway excavation and once in 10–15 days after the excavation until all stations are fully stabilized with a stable deformation rate of less than 0.05 mm/d. The in situ monitored values for two monitoring sections are presented in Figure4.Figure 4
Displacement curves of roadway surface measured in each monitoring section.As indicated in Figure4, it could be found that rock mass surrounding roadway represents significant time-dependent deformation both during the excavation and afterwards. Moreover, the crown settlement and the sidewalls convergence at the end of the roadway excavation are, respectively, 14.2 mm and 37.6 mm in Section I and 12.4 mm and 34.2 mm in Section II. Afterwards, the surrounding rock mass was in a steady state with a stable deformation rate after approximately 8 months (230 days) from the start of roadway excavation. In detail, stable times of surrounding rock mass in Section I and Section II were surveyed as 210 days and 225 days, respectively. At the same time, the final values of the crown settlement and the sidewalls convergence were 22.2 mm and 65.8 mm in Section I and 21.8 mm and 63.3 mm in Section II, respectively.


## 4.2. Results and Discussion

It can be found from Figure4 that measured displacement data are given as accumulative values of the crown settlement and the sidewalls convergence. Therefore, we established the original displacement sequence X(0) by applying 1-IAGO on the displacement data measured in a period of the roadway excavation and constructed three NGMs (i.e., INGVM, NGM(1,1), and NGVM) for comparison, respectively. Then, the values of the parameters for different NGMs corresponding to measurement components of tunnel displacement in each monitoring section were estimated by the least-square method, which is indicated in Table 2.Table 2
The estimated values of parameters for NGMs by the least-square method.
Grey modelModel parameterSection ISection IIzAcBCzAcBCNGM(1, 1)A0.00850.00510.00670.0058B1.27433.21291.36744.0326NGVM and INGVMA−0.1505−0.1701−0.1096−0.1132B−2.8344e-3−1.1378e-03−1.0417e-3−3.0344e-4By using the estimated values of parameters in Table2, the fitting results of the prediction curves by the three compared models with the monitoring curves of roadway displacement at two monitoring sections are shown in Figure 5, and the values of two accuracy test criteria (MAPE and RMSPE) on fitting and prediction results of different prediction models are all listed in Table 3.Figure 5
Comparison curves between in situ monitored data and predicted results by three compared models. (a) Section I. (b) Section II.
(a)(b)Table 3
Accuracy check of fitting and prediction results for different prediction models.
Monitoring componentNGM(1,1)NGVMINGVMMAPERMSPEMAPERMSPEMAPERMSPEFittingSection IzA6.148.6515.5716.725.056.27cBC7.657.1018.0319.856.046.67Section IIzA7.418.1116.5918.025.377.67cBC8.7710.0318.8319.747.048.34Mean7.498.4717.4617.265.887.24PredictionSection IzA7.068.3514.2315.963.694.16cBC8.6510.1914.4015.673.644.13Section IIzA7.3510.5915.7217.154.766.09cBC8.1911.8013.5916.514.365.71Mean7.8110.2314.4916.324.115.02As shown in Figure5 and Table 3, it can be found that the MAPE and RMSPE between in situ measured data and predicted results by three compared models are all within 20%. Moreover, it can be seen from Figure 5 that, in prediction interval, INGVM and NGVM simulate a similar trend (i.e., S-type curve) of actual data, but the simulation of NGM(1,1) has an ascending trend. This may be the reason why the NGVM was applied in this study.From the viewpoint of error testing, it is obvious that the prediction accuracy of the INGVM is higher than that of the traditional NGMs in displacement prediction of the time-dependent rock mass. Especially, MAPE and RMSPE of INGVM were 5.88% and 7.24% in the fitting interval and 4.11% and 5.02% in prediction interval, respectively, indicating that prediction errors of INGVM are all smaller than those of the NGM(1,1) and NGVM. In addition, the mean values of MAPE and RMSPE of INGVM in the prediction interval are reduced from 14.49% to 16.32% to 4.11% and 5.02%, meaning the INGVM is more effective and applicable than the traditional NGVM. It is shown that Fourier series helped INGVM to further reduce the periodic residual error compared with traditional NGMs. Meanwhile, we can also see from Figure5 and Table 3 that the errors in fitting interval with relation to NGVM and INGVM are slightly lower than those in prediction interval because the fitting situation of displacement data measured in situ to be used at fitting interval deviates from the ideal S-type.At the same time, it could be also predicted that stable time, final crown settlement, and sidewalls convergence predicted by INGVM are 218 days, 22.43 mm and 66.37 mm in Section I, and 235 days, 22.01 mm and 62.83 mm in Section II, respectively, indicating that prediction results by the proposed model with criterion for the stable displacement of 0.05 mm/d are in good agreement with in situ monitored data. It can be seen that prediction results’ stable time and final displacement by using the proposed model tend to slightly increase than real values.From the above comparative analysis, we can conclude that the proposed INGVM method is a better fitting model compared with NGM(1,1) and NGVM to predict the final displacement and stable time of time-dependent rock mass surrounding roadways as a moderate accuracy, using only monitored in situ data during excavation. Unfortunately, it may not be suitable for the cases when the deformation of rock mass surrounding unsupported and supported underground structures is not just unstable but either continuously increased according to time changes or suddenly increased. However, if displacement histories of soil tunnel and highly jointed rock tunnel show the tendency gradually to be stabilized in S-type curve by rational excavation method, suitable support system, and ground reinforcement, proposed model may be also used for displacement prediction of such rock tunnels.


## 5. Conclusions

This study focused on predicting the final displacement of tunnel in time-dependent rock mass using improved nonequidistant grey Verhulst model (INGVM), in which Fourier series is used to reduce the periodic residual error.The maximum stable deformation rate was used as a criterion for the stable displacement of rock mass surrounding unsupported drift in a coal mine, which was equal to 0.05 mm/d because the deformation rate is an important index to reflect the overall deformation and the stability of the rock mass. However, this criterion can be set up as a lower value for more important underground structures (e.g., powerhouse caverns and mine primary roadways).Significantly, when the deformation of surrounding rock mass becomes gradually stable as a saturated curve with S-type, we found the fact that the tunnel displacement can be fully predicted by INGVM with proper stable criterion, based on only displacement data monitored during tunnel excavation, whereas NGM(1,1) is imperfect for such case discussed above.According to Figure5 and Table 3, the mean MAPE and RMSPE between displacements predicted by INGVM and in situ data in prediction interval are, respectively, 4.11% and 5.02%, which are higher than NGM(1,1) and NGVM. Therefore, it is found that the proposed grey model can be successfully applied to displacement prediction of time-dependent rock mass and to stability evaluation and support design of the underground structures, as well as mine openings and stopes, in which the displacement increases in the curve with S-type or its increment is in the saturation stage.Future researchers can be using different equations or different methodologies like the Markov Chain and PSO to improve the accuracy of NGVM. Furthermore, the proposed model can be applied in many other industries with the S-type saturated data.



---

*Source: 3241171-2022-03-26.xml*