Analysis and Warning Prediction of Tunnel Deformation Based on Multifractal Theory

: To better analyze the fluctuation characteristics and development law of tunnel deformation data, multifractal theory is applied to tunnel deformation analysis. That is, the multifractal detrended fluctuation analysis (MF-DFA) model is first utilized to carry out the multifractal characterization of tunnel deformation data. Further, Mann–Kendall (M–K) analysis is utilized to construct the dual criterion ( ∆ α indicator criterion and ∆ f ( α ) indicator criterion) for the tunnel deformation early warning study. In addition, the particle swarm optimization long-short-term memory (PSO-LSTM) prediction model is used for predicting tunnel settlement. The results show that, in reference to the tunnel warning level criteria and based on the Z-value results of the indicator criterion, the warning level of all four sections is class II. At the same time, through the analysis of tunnel settlement predictions, the PSO-LSTM model has a better prediction effect and stability for tunnel settlement. The predicted results show a slow increase in tunnel settlement over the next 5 days. Finally, the tunnel warning level and the predicted results of tunnel settlement are analyzed in a comprehensive manner. The deformation will increase slowly in the future. Therefore, monitoring and measurement should be strengthened, and disaster preparedness plans should be prepared.


Introduction
With the continuous development of transportation construction in China, the scale and technical level of tunnel construction have reached new heights [1].By the end of 2022, the national railroad operating mileage reached 155,000 km, boasting 17,873 operational railroad tunnels with a cumulative length of approximately 22,000 km.Of these, 259 extra-long railroad tunnels and 12 tunnels exceeding 20 km in length are in operation [2].Notably, over the past five years, railroad tunnel development has experienced remarkable acceleration.This surge is attributed to the rapid growth of the high-speed railroad industry and the consistent enhancement of tunnel trimming technology, resulting in the emergence of numerous tunnels traversing diverse types of special strata [3].However, concurrent with this progress, there exist design and construction irregularities in high-speed railroad tunnels that may lead to deformation, water leakage, cracking, and other issues [4][5][6][7][8].
Fractal theory, pioneered by the French mathematicians Mandelbrot et al. [27], describes chaotic, complex, and self-similar systems with regularity.Researchers, such as Zuo et al., have applied fractal theory to analyze and predict tunnel surface settlement, evaluating surface settlement stability [28].Ye et al. established a model based on fractal theory to study the interaction between microstructure evolution and tunnel leakage behavior, offering insight into the tunnel leakage mechanism [29].Additionally, MF-DFA, initially proposed by Grassberger [30], finely describes the volatility of tunnel deformation monitoring data at different levels.At present, the multifractal theory has been widely used.Lei et al. used multifractal theory and subcomponent combination prediction to synthesize the warning level of landslides [31].Mao et al. revealed early warning signals of rock slope deformation based on the multiple fractal time-varying response characteristics of micro-seismic waveforms associated with rock rupture [32].Zhou et al. applied multifractal theory to the deformation pattern analysis of dams and utilized MF-DFA, multivariate multifractal detrended fluctuation analysis (MV-MFDFA), and asymmetric multifractal detrended fluctuation analysis (A-MFDFA) to resolve the multiple fractal features of deformation states and their asymmetry [33].
Accurate prediction of tunnel deformation is crucial for disaster prevention and early warning.Currently, soil settlement prediction methods include empirical methods [34][35][36], theoretical methods [37][38][39][40], modeling methods [41][42][43][44], and machine learning methods [45][46][47][48].Long short-term memory (LSTM) is particularly effective in handing time-related information and is widely used in temporal tunnel settlement prediction [49].Li et al. employed three machine learning algorithms (LSTM, RF, and GRU) to forecast surface settlement.Their findings indicated that the LSTM algorithm demonstrated superior accuracy in predicting maximum surface settlement and effectively anticipated settlement progression in various strata [50].Cao et al. proposed the complete ensemble empirical mode decomposition with adaptive noise long short-term memory (CEEMDAN-LSTM) model, combining high prediction accuracy with acceptable computational efficiency [51].Duan et al. utilized the auto regressive integrated moving average (ARIMA) and LSTM models for predicting structural deformation trends during tunnel operation, determining that LSTM performs better with high data quality and sufficient samples [52].
From the above research results, it can be concluded that research on tunnel deformation analysis and early warning prediction is necessary.Therefore, this paper analyzes tunnel deformation data.Firstly, the deformation rate and the deformation of the tunnel are investigated by multiple fractal characterizations using MF-DFA.Secondly, the Mann-Kendall (M-K) analysis method [53,54] is used to evaluate the early warning classification of tunnel deformation.Furthermore, the prediction of tunnel deformation is realized by the particle swarm optimization LSTM (PSO-LSTM) prediction model.Finally, the results of tunnel deformation early warning classification and the tunnel deformation prediction results are jointly addressed.The tunnel deformation law is evaluated comprehensively to provide theoretical guidance for its monitoring, measurement, and prevention.This study calculates and analyzes the multifractal features of tunnel monitoring data by integrating the M-K method to determine the tunnel's warning level.Furthermore, the PSO-LSTM model is found to exhibit high accuracy in predicting tunnel deformation, offering a novel perspective for advancing research on tunnel early warning prediction.

Materials and Methods
The process of tunnel deformation analysis and early warning prediction, established in this paper based on multiple fractal theory, is illustrated in Figure 1.The steps can be summarized as follows: First, tunnel deformation monitoring and measurement are conducted to obtain tunnel deformation data.Second, the monitoring data are segmented into 12 groups, and multifractal eigenvalues are calculated and analyzed using the M-K method.The warning level is then determined in conjunction with the warning level criteria.Third, the parameters of the LSTM prediction model are optimized using the PSO optimization algorithm, and the PSO-LSTM prediction model is formulated for predicting tunnel settlement.Fourth, the warning level and prediction results are jointly analyzed to obtain the warning outcomes.Finally, corresponding preventive measures are implemented based on the results of the warning level.
in this paper based on multiple fractal theory, is illustrated in Figure 1.The steps can be summarized as follows: First, tunnel deformation monitoring and measurement are conducted to obtain tunnel deformation data.Second, the monitoring data are segmented into 12 groups, and multifractal eigenvalues are calculated and analyzed using the M-K method.The warning level is then determined in conjunction with the warning level criteria.Third, the parameters of the LSTM prediction model are optimized using the PSO optimization algorithm, and the PSO-LSTM prediction model is formulated for predicting tunnel settlement.Fourth, the warning level and prediction results are jointly analyzed to obtain the warning outcomes.Finally, corresponding preventive measures are implemented based on the results of the warning level.
This comprehensive approach provides a systematic framework for tunnel deformation analysis, early warning prediction, and subsequent preventive action.

Project Overview and Monitoring Data
The high-speed rail tunnel has a total length of 254.9 m, with the starting and ending mileage marked as DK115 + 133.0~DK115 + 387.9.It is positioned between the Sangujian No. 1 Tunnel and the Xikeng bridge.Geographically, it is situated on the western side of the mountain, in Tamkou Village, Xidi Town, Yixian County, Huangshan City.The surface vegetation on the mountain is characterized by bamboo forests, low shrubs, trees, and low weeds.The terrain is steep, with a natural slope ranging from 30° to approximately 55°.The tunnel is surrounded by dense mountain vegetation, and the entrance and exit are marked by ravine confluence, making traffic inconvenient.The perimeter rock grade is classified as Ⅳ, and the excavation employs the three-step method.
The tunnel is situated in the low hills, originating from Huangshan Mountain.It is characterized by undulating terrain with mountain elevations ranging from 250 m to 325 m and a maximum elevation difference of approximately 75 m.The mountain exhibits a monoclinic structure, with steeper south and north slopes.The tunnel primarily traverses carbonaceous mudstone of the Cambrian Hotang formation, with a maximum tunnel depth of about 65 m.(Refer to Figure 2 for a visualization of the tunnel cross-section.)This comprehensive approach provides a systematic framework for tunnel deformation analysis, early warning prediction, and subsequent preventive action.

Project Overview and Monitoring Data
The high-speed rail tunnel has a total length of 254.9 m, with the starting and ending mileage marked as DK115 + 133.0~DK115 + 387.9.It is positioned between the Sangujian No. 1 Tunnel and the Xikeng bridge.Geographically, it is situated on the western side of the mountain, in Tamkou Village, Xidi Town, Yixian County, Huangshan City.The surface vegetation on the mountain is characterized by bamboo forests, low shrubs, trees, and low weeds.The terrain is steep, with a natural slope ranging from 30 • to approximately 55 • .The tunnel is surrounded by dense mountain vegetation, and the entrance and exit are marked by ravine confluence, making traffic inconvenient.The perimeter rock grade is classified as IV, and the excavation employs the three-step method.
The tunnel is situated in the low hills, originating from Huangshan Mountain.It is characterized by undulating terrain with mountain elevations ranging from 250 m to 325 m and a maximum elevation difference of approximately 75 m.The mountain exhibits a monoclinic structure, with steeper south and north slopes.The tunnel primarily traverses carbonaceous mudstone of the Cambrian Hotang formation, with a maximum tunnel depth of about 65 m.(Refer to Figure 2 for a visualization of the tunnel cross-section.) To ensure the safety of tunnel construction operations, a total station is utilized for manual measurement of tunnel settlement and convergence.This paper involves monitoring the settlement of the tunnel vault and tunnel convergence utilizing the rear rendezvous method and the opposite side measurement method, respectively.For arch settlement monitoring, a single monitoring point is positioned at the arch of each section.In contrast, for tunnel convergence monitoring, two monitoring points are symmetrically arranged on both sides of the tunnel, employing the opposite side measurement method.Consequently, three monitoring points are allocated in each section, with a monitoring frequency of 2 times/day.Further details are depicted in Figure 3.To ensure the safety of tunnel construction operations, a total station is utilized for manual measurement of tunnel settlement and convergence.This paper involves monitoring the settlement of the tunnel vault and tunnel convergence utilizing the rear rendezvous method and the opposite side measurement method, respectively.For arch settlement monitoring, a single monitoring point is positioned at the arch of each section.In contrast, for tunnel convergence monitoring, two monitoring points are symmetrically arranged on both sides of the tunnel, employing the opposite side measurement method.Consequently, three monitoring points are allocated in each section, with a monitoring frequency of 2 times/day.Further details are depicted in Figure 3.As shown in Figures 4 and 5, this paper examines the monitoring data regarding tunnel settlement and convergence in the section from DK115 + 261~DK115 + 290.Additionally, the data from monitoring points are counted.The cumulative deformationtime diagram and deformation rate-time diagram of tunnel settlement and convergence are obtained, respectively.
In Figure 4, it is evident that both the cumulative tunnel settlement deformation and the cumulative tunnel convergence deformation exhibit an increasing trend over time.Notably, after 1000 h, the growth of these deformation measures reaches a turning point, transitioning into a slower growth trend.Throughout the monitoring period, the data display both significant and minor fluctuations, which can be effectively analyzed using multiple fractal theory.
Looking at Figure 5, it is apparent that both the tunnel settlement deformation rate and the tunnel convergence deformation rate demonstrate significant and minor fluctuations, displaying the characteristic traits of multiple fractals.These fluctuations can be effectively analyzed using multiple fractal theory.To ensure the safety of tunnel construction operations, a total station is utilized for manual measurement of tunnel settlement and convergence.This paper involves monitoring the settlement of the tunnel vault and tunnel convergence utilizing the rear rendezvous method and the opposite side measurement method, respectively.For arch settlement monitoring, a single monitoring point is positioned at the arch of each section.In contrast, for tunnel convergence monitoring, two monitoring points are symmetrically arranged on both sides of the tunnel, employing the opposite side measurement method.Consequently, three monitoring points are allocated in each section, with a monitoring frequency of 2 times/day.Further details are depicted in Figure 3.As shown in Figures 4 and 5, this paper examines the monitoring data regarding tunnel settlement and convergence in the section from DK115 + 261~DK115 + 290.Additionally, the data from monitoring points are counted.The cumulative deformationtime diagram and deformation rate-time diagram of tunnel settlement and convergence are obtained, respectively.
In Figure 4, it is evident that both the cumulative tunnel settlement deformation and the cumulative tunnel convergence deformation exhibit an increasing trend over time.Notably, after 1000 h, the growth of these deformation measures reaches a turning point, transitioning into a slower growth trend.Throughout the monitoring period, the data display both significant and minor fluctuations, which can be effectively analyzed using multiple fractal theory.
Looking at Figure 5, it is apparent that both the tunnel settlement deformation rate and the tunnel convergence deformation rate demonstrate significant and minor fluctuations, displaying the characteristic traits of multiple fractals.These fluctuations can be effectively analyzed using multiple fractal theory.As shown in Figures 4 and 5, this paper examines the monitoring data regarding tunnel settlement and convergence in the section from DK115 + 261~DK115 + 290.Additionally, the data from monitoring points are counted.The cumulative deformation-time diagram and deformation rate-time diagram of tunnel settlement and convergence are obtained, respectively.
In Figure 4, it is evident that both the cumulative tunnel settlement deformation and the cumulative tunnel convergence deformation exhibit an increasing trend over time.Notably, after 1000 h, the growth of these deformation measures reaches a turning point, transitioning into a slower growth trend.Throughout the monitoring period, the data display both significant and minor fluctuations, which can be effectively analyzed using multiple fractal theory.
Looking at Figure 5, it is apparent that both the tunnel settlement deformation rate and the tunnel convergence deformation rate demonstrate significant and minor fluctuations, displaying the characteristic traits of multiple fractals.These fluctuations can be effectively analyzed using multiple fractal theory.

MF-DFA
The MF-DFA model, which belongs to the multiple non-uniform fractal method, not only reveals the multiple fractal characteristics within the deformation sequence but also effectively evaluates their deformation trends [55].The calculation steps are as follows: Step 1: the data presented in Figures 4 and 5, respectively, are set as a time series x(t) with a series length of N, a series mean of x, and a series of cumulative deviations of x(t) with respect to x as y(t): Step 2: Set a time scale s.Equalize the sequence y(t) in terms of s, dividing it into a total of m equal-length consecutive and non-overlapping sub-intervals, m = int(N/s).
In practical arithmetic, N may not necessarily be divisible, resulting in potential tail data redundancy.Thus, alongside positive-order division, a reverse-order processing method is also employed simultaneously.The division operation is repeated from the end of the sequence to obtain 2m subintervals.
Step 3: Fit a trend to each subinterval and subtract the trend portion from the original.Obtain the corresponding residual series denoted as z v (t): where y v (t) represents the subinterval, and p k v (t) is a kth order fitting polynomial to the vth subinterval.v ranges from 1 to 2 m, and t ranges from 1 to s.
Step 4: calculate the mean square deviation F 2 (s, v) of the residual sequence z v (t): Step 5: Optimize the traditional MF-DFA division by employing a sliding window, which involves moving the values through a window of a specific length along the sequence at a defined step size.This helps to minimize pseudo-fluctuations in the data and maximize the utilization of data information.
The window length is denoted as s, the sequence length as N, and the sliding step is set to 1.The number of subintervals obtained in one run is N − S + 1, and the q-order fluctuation function is calculated according to Equation (4): Step 6: Repeat the previous steps to generate a series of point values for s-F q (s).If this time series exhibits a long-range correction, then F q (s) has a power-law relationship with s, as shown in Equation ( 5): By taking the logarithms of both sides of the previous equation, we obtain Equation ( 6): where F q (s) represents the q-order fluctuation function of the series, h(q) is the correspond- ing generalized Hurst exponent, and b is a constant coefficient.A plot of lgF q (s)-lgs was created and fitted to determine the generalized Hurst ex- ponent h(q).A fixed h(q) suggests a mono-fractal sequence without multifractal features.When h(q) < 0.5, the data sequence behaves as a memory process with inverse persistence.
A value of h(q) = 0.5 indicates uncorrelated stochastic behavior.For h(q) > 0.5, the data sequence behaves as a memory process with positive persistence.Values of h(q) > 1 indicate long-range positively correlated processes with strong non-stationarity.
Step 7: the multifractal spectrum f (α), which characterizes the fractal intensity and singularity of the time series, can typically be determined using Equation (7): where τ(q) is the Renyi index, also known as the scalar function.If it is a nonlinear upconvex function of q, the displacement sequence exhibits multifractal characteristics.If it is a linear function of q, the displacement sequence exhibits single fractal characteristics.
where the ∆α parameter is mainly used to evaluate the width of the multifractal spectrum of the data deformation sequence.As the value of ∆α increases, the intensity of multiple fractals becomes stronger, leading to more intense fluctuations.The ∆ f (α) parameter is mainly used to evaluate the proportion of large and small fluctuations in the waveform of the data deformation sequence.As the value of ∆ f (α) decreases, the proportion of large fluctuation waveforms increases.

M-K Test Method
The M-K test method is a non-parametric method.Non-parametric tests are also known as non-distributional tests.It has the advantage that the sample does not have to follow a particular distribution and is not affected by a few outliers.It is suitable for type and order variables, and it is relatively simple to calculate [56].The specific calculation steps are as follows: Step 1: The 12 sets of eigenvalues obtained from each set of data by the formula in 2.2 are set as a time series of 12 sample sizes {x 1 , x 2 , . . . . . . ,x 12 }.For all k(j ≤ n and k ̸ = j), the distributions of x k and x j are different and the difference function Sgn x j − x k is computed: Step 2: calculate the test statistic S: Step 3: S is normally distributed with mean 0. Calculate the variance Var(S): Step 4: calculate the standard normal statistical variable Z: , S < 0 (12) Step 5: The trend characteristics of the corresponding evaluation object can be judged by the size of Z.The Z α value is the critical value under the condition of corresponding significant level α.In this paper, the significance test with 99% confidence level is selected.That is, for a significance level α = 0.01, then Z 0.01 = 2.32.
If Z ≥ Z α , it means that the evaluation object has an increasing trend, and if it is larger, it indicates a stronger trend.
If −Z α < Z < Z α , it means that the evaluation object has a smooth trend.If Z ≤ −Z α , it means that the evaluation object has a decreasing trend, and if it is smaller, it indicates a stronger trend.

PSO-LSTM Prediction Modeling
2.4.1.LSTM LSTM, a special variant of recurrent neural networks, was introduced by Hochreiter in 1997 [57].It incorporates temporal memory units capable of learning dependency information across various time periods in a time series.This network is particularly effective at processing and predicting intervals and delayed events within time series data.The specific calculation steps are as follows: Step 1: tunnel monitoring statistics are normalized and calculated by the formula: where x(i, j) represents the original data, x ′ (i, j) denotes the normalized data, and x max (j) and x min (j) are the maximum and minimum values of the parameters of the jth model, respectively.
Step 2: LSTM consists of multiple cells, and the formulae for the forgetting gate, input gate, and output gate in each cell are: where h t−1 represents the hidden state from the previous moment, and σ denotes the sigmoid function.f t , i t , and O t are the outcomes of the state settlements for the oblivious gate, the input gate, and the output gate, respectively.W f , W i , and W o are the weight matrices for the oblivious gate, the input gate, and the output gate, respectively.b f , b i , and b o represent the biases for oblivious gate, the input gate, and the output gate, respectively.
Step 3: C t is a vector of candidate values, and the product of the input values and the vector of candidate values is used to update the cell state, calculated as: where W c represents the weight matrix for the input unit state, while b c is the bias term for the input unit state.The activation function is tan h.o t denotes the neuron output value.h t is the current moment's hidden state.

PSO
PSO is inspired by the behavior of birds foraging in nature.In this algorithm, each optimization problem is referred to as a "particle", and the process of particle swarm optimization is likened to the foraging behavior of a bird flock.These particles are equipped with a memory function to store optimization information, which is then shared within the flock.The optimal position information is subsequently selected as the optimal position information of the entire flock [58].
During each iteration, the position information of each particle is updated.The change in the position information x k id for the ith evolutionary example consists of a linear summation of the previous position information x k−1 id and the previous velocity information Furthermore, each particle possesses a velocity vector that dynamically changes in real time due to various factors.It includes the example velocity information V k id of the ith evolution, with the previous velocity WV k−1 id , the individual optimal position , and the population optimal position C 2 r 2 gbest id − x k−1 id .The formula is as follows: where C 1 and C 2 are acceleration constants, and r 1 and r 2 are random functions.

PSO-LSTM
LSTM demonstrates exceptional feature extraction capabilities from data exhibiting spatial and temporal correlation.This study applies the LSTM model for predicting tunnel settlement data.In order to further reduce the model error, PSO is employed to optimize the LSTM parameters and establish the PSO-LSTM model.
The LSTM model is constructed and trained using the training set.Subsequently, the prediction results and the actual values are compared and analyzed for error in the test set.The accuracy of the model is assessed using three indicators: the coefficient of determination (R 2 ), the mean absolute error (MAE), and the root mean square error (RMSE).The calculation formula is: where n represents the number of predicted outcomes, y o (i) represents the true outcome, y c (i) represents the predicted outcome, and y o represents the mean of the true values.

Characterization of Tunnel Deformation Rate Data
Before investigating the multifractal characteristics of the deformation rate of tunnel settlement and convergence, the basic statistical characteristics of the sample data should be understood.Here, the descriptive statistics of the monitoring data in Figure 2 are shown in Table 1.
According to the K-S criterion [59,60], the four evaluation indexes for average, standard deviation, skewness, and kurtosis are considered: 1.
From the average: the average tunnel deformation rate of each cross-section is greater than 0, which shows positive deformation.

2.
From the standard deviation: The standard deviation values of the tunnel deformation rate at each cross-section are close.In the deformation rate of tunnel settlement, the standard deviation at DK115 + 280 is the largest, which can be visualized in Figure 2a.
In the deformation rate of tunnel convergence, the standard deviation at DK115 + 270 is the largest, which can be visualized in Figure 2b.The standard deviations are all close to or greater than the mean, indicating a greater degree of dispersion in the data.

3.
From the skewness: The skewness value of each cross-section is greater than 0, indicating that the data distribution is right-skewed.That is, the dispersion on the right side of the mean is stronger than that on the left side, showing a certain degree of long tail on the right side.All skewness values are greater than the mean and standard deviation, indicating a greater degree of skewness in the data distribution.4.
From the kurtosis: the kurtosis value for each cross-section is greater than 0, indicating that the peaks of the data are steeper than the peaks of the normal distribution resulting in a spiky state.In summary, the combined Jarque-Bera (J-B) statistic results are all significantly greater than 0. Under the joint effect of kurtosis and skewness, the deformation rate of each crosssection does not follow a normal distribution.This means that the deformation rate of the tunnel has obvious fractal characteristics, which can be investigated by the relevant methods of fractal theory.

MF-DFA key parameter settings
The values of the parameters in multifractal theory affect the calculation results in different ways.The features of the non-stationary time series differ significantly in various scenarios, including the length of the signal time window and the fluctuation trend.Therefore, the key parameters need to be tried and predicted to estimate more reliable results.
In order to effectively use multiple fractal theory to analyze the fractal characteristics of the tunnel displacement deformation rate, a double logarithmic scatterplot of lgF q (s) − lgs is plotted by varying the fluctuation order q (q ∈ [−10, 10]).Least squares fitting was used with the slope of the generalized Hurst exponent h(q).Taking a set of multifractal features of tunnel settlement deformation rate and tunnel convergence deformation rate as an example, Figure 6 shows the trend of the q-order fluctuation function lgF q (s) − lgs through double logarithmic fitting.
The curve-fitting for the tunnel section settlement rate in Figure 6a shows better results when the value of lgs ranges from 1.5 to 1.65.At this time, the multifractal time length s is set as smin = 32 and smax = 45.
The curve-fitting for the tunnel section convergence rate in Figure 6b demonstrates improved performance when the value of lgs ranges from 1.1 to 1.4.At this time, the multifractal time length s is set as smin = 13 and smax = 25.The curve-fitting for the tunnel section settlement rate in Figure 6a shows better results when the value of lgs ranges from 1.5 to 1.65.At this time, the multifractal time length s is set as smin = 32 and smax = 45.
The curve-fitting for the tunnel section convergence rate in Figure 6b demonstrates improved performance when the value of lgs ranges from 1.1 to 1.4.At this time, the multifractal time length s is set as smin = 13 and smax = 25.

Multifractal characterization of tunnel deformation rates
According to the settings of the above parameters, the data series for tunnel settlement deformation rate and tunnel convergence deformation rate were subjected to multifractals using MATLAB software (R2018b), respectively.The variation of the generalized Hurst index of the tunnel deformation rate data series is shown in Figure 7.The variation of the scalar function   of the tunnel deformation rate data series is shown in Figure 8.The multifractal spectrum of the tunnel deformation rate data series is shown in Figure 9.As can be seen from Figure 7, as q ranges from −10 to 10, the generalized Hurst index of each section measured data is non-constant.Rather, it exhibits a nonlinear decreasing trend with the change of q, suggesting that the measured data of each cross-section

Multifractal characterization of tunnel deformation rates
According to the settings of the above parameters, the data series for tunnel settlement deformation rate and tunnel convergence deformation rate were subjected to multifractals using MATLAB software (R2018b), respectively.The variation of the generalized Hurst index of the tunnel deformation rate data series is shown in Figure 7.The variation of the scalar function τ(q) of the tunnel deformation rate data series is shown in Figure 8.The multifractal spectrum of the tunnel deformation rate data series is shown in Figure 9.The curve-fitting for the tunnel section settlement rate in Figure 6a shows be results when the value of lgs ranges from 1.5 to 1.65.At this time, the multifractal t length s is set as smin = 32 and smax = 45.
The curve-fitting for the tunnel section convergence rate in Figure 6b demonstr improved performance when the value of lgs ranges from 1.1 to 1.4.At this time, multifractal time length s is set as smin = 13 and smax = 25.

Multifractal characterization of tunnel deformation rates
According to the settings of the above parameters, the data series for tu settlement deformation rate and tunnel convergence deformation rate were subjecte multifractals using MATLAB software (R2018b), respectively.The variation of generalized Hurst index of the tunnel deformation rate data series is shown in Figu The variation of the scalar function   of the tunnel deformation rate data serie shown in Figure 8.The multifractal spectrum of the tunnel deformation rate data seri shown in Figure 9.As can be seen from Figure 7, as q ranges from −10 to 10, the generalized Hurst in of each section measured data is non-constant.Rather, it exhibits a nonlinear decrea trend with the change of q, suggesting that the measured data of each cross-sec As can be seen from Figure 7a, different fluctuation orders q, the generalized Hurst exponent curves of tunnel settlement rate at DK115 + 290 are concentrated in the lower fluctuations of the other sections.This indicates a weak multifractal nature.
As can be seen from Figure 7b, at different fluctuation orders q, the generalized Hurst exponent curves of the tunnel convergence rate of DK15 + 261 concentrate in the lower fluctuations of the other sections.This indicates a weak multifractal nature.On the basis of h(q), the Renyi index, i.e., the scaling function   , is calculated.As can be seen from Figure 8, the consistency of the scalar function of the monitoring data of each cross-section is good, and the central part is up-convex, which satisfies  0 1.Additionally, there is an overall nonlinear relationship, which further confirms that the monitoring data of each section have multifractal characteristics.As can be seen from Figure 9, each image of the multifractal spectrum shows a singlepeak convex distribution, which resembles a quadratic function curve.The local scales of the multiple fractals of the deformation rate of the tunnel in each section are not constant, reflecting the diversity of local variations at different moments.The singularity intensity  is mainly concentrated on the two sides of the image, reflecting the uneven distribution of the fractal structure of the data series.The uneven distribution of  also confirms the multifractal characteristics of the measured point series;   characterizes the fractal On the basis of h(q), the Renyi index, i.e., the scaling function   , is calculated.As can be seen from Figure 8, the consistency of the scalar function of the monitoring data of each cross-section is good, and the central part is up-convex, which satisfies  0 1.Additionally, there is an overall nonlinear relationship, which further confirms that the monitoring data of each section have multifractal characteristics.As can be seen from Figure 9, each image of the multifractal spectrum shows a singlepeak convex distribution, which resembles a quadratic function curve.The local scales of the multiple fractals of the deformation rate of the tunnel in each section are not constant, reflecting the diversity of local variations at different moments.The singularity intensity  is mainly concentrated on the two sides of the image, reflecting the uneven distribution of the fractal structure of the data series.The uneven distribution of  also confirms the multifractal characteristics of the measured point series;   characterizes the fractal As can be seen from Figure 7, as q ranges from −10 to 10, the generalized Hurst index of each section measured data is non-constant.Rather, it exhibits a nonlinear decreasing trend with the change of q, suggesting that the measured data of each cross-section display clear multifractal characteristics.When h(q) < 0.5, the data series exhibits a memory process with inverse persistence, while h(q) > 0.5 indicates a memory process with positive persistence.
As can be seen from Figure 7a, at different fluctuation orders q, the generalized Hurst exponent curves of tunnel settlement rate at DK115 + 290 are concentrated in the lower fluctuations of the other sections.This indicates a weak multifractal nature.
As can be seen from Figure 7b, at different fluctuation orders q, the generalized Hurst exponent curves of the tunnel convergence rate of DK15 + 261 concentrate in the lower fluctuations of the other sections.This indicates a weak multifractal nature.
On the basis of h(q), the Renyi index, i.e., the scaling function τ(q), is calculated.As can be seen from Figure 8, the consistency of the scalar function of the monitoring data of each cross-section is good, and the central part is up-convex, which satisfies τ(0) = −1.Additionally, there is an overall nonlinear relationship, which further confirms that the monitoring data of each section have multifractal characteristics.
As can be seen from Figure 9, each image of the multifractal spectrum shows a singlepeak convex distribution, which resembles a quadratic function curve.The local scales of the multiple fractals of the deformation rate of the tunnel in each section are not constant, reflecting the diversity of local variations at different moments.The singularity intensity α is mainly concentrated on the two sides of the image, reflecting the uneven distribution of the fractal structure of the data series.The uneven distribution of α also confirms the multifractal characteristics of the measured point series; f (α) characterizes the fractal dimension of the subintervals of the data series with the same singularity index α, which is correlated with the distributional characteristics and fractal intensity of the data series.The statistics of the multifractal features of the tunnel settlement rate and tunnel convergence rate measurement data sequences are calculated from Equation ( 8), as shown in Table 2.As can be seen from Figure 9a, the fractal spectrum of the tunnel settlement rate measurement data series for each section shows a clear right hook, indicating that the influence of small fluctuations is dominant.The results indicate that the tunnel settlement rate primarily exhibits minor fluctuations throughout the entire monitoring process.During the tunnel construction process, the settlement of the tunnel vault is minimally impacted by excavation and other external influences, remaining within the normal range of settlement changes.
As can be seen from Figure 9b, the sequence fractal spectrum of the tunnel convergence rate measurement data of DK115 + 290 is basically symmetrical, with good overall synergy and a stable development state.The sequence fractal spectrum of the tunnel convergence rate measurement data of DK115 + 261 shows an obvious right hook, indicating that the influence of small fluctuations is dominant.The fractal spectrum of the sequence of the tunnel convergence rate measurement data of DK115 + 270 and DK115 + 280 shows an obvious left hook, indicating that the influence of large fluctuations is dominant.
The fluctuation of the tunnel convergence rate in section DK115 + 290 remains stable throughout the entire monitoring process.During the tunnel construction process, the tunnel convergence in this section is minimally affected by excavation and other external influences, remaining within the normal range of convergence changes.Similarly, in section DK115 + 261, the tunnel convergence rate primarily exhibits minor fluctuations, with the tunnel convergence being less affected by excavation and other external influences, and remaining within the normal range of convergence changes.However, in sections DK115 + 270 and DK115 + 280, the tunnel convergence rate predominantly displays significant fluctuations.During the tunnel construction process in these sections, tunnel convergence is greatly affected by excavation and other external influences, leading to sudden changes.Consequently, it is imperative to reinforce monitoring in these areas.
Comparison of multifractal spectral width (∆α).For the tunnel settlement rate measurement data series, the multifractal spectral width ∆α = 0.5547 for DK115 + 290 is the maximum value.It shows that its multifractal intensity is larger, and the fluctuation is more intense and complex.The multifractal spectral width ∆α = 0.3946 of DK115 + 270 is the minimum value.It shows that its multifractal intensity is smaller, and the fluctuation is smoother.For the tunnel convergence rate measurement data series, the multifractal spectrum width ∆α = 0.5411 for DK115 + 270 is the maximum value.It shows that its multifractal intensity is larger, and the fluctuation is more intense and complex.The multifractal spectral width ∆α = 0.3674 of DK115 + 290 is the minimum value.It shows that its multiple fractal intensity is smaller, and the fluctuation is smoother.The results suggest that a larger ∆α corresponds to more dramatic and intricate fluctuations, while smaller data changes indicate that tunnel deformation is less affected by external influences such as excavation.Overall, the multiple fractal spectral widths of the measured data series of the tunnel settlement rate and tunnel convergence rate of the same section show an inverse change relationship, and the sum of the two is stable between 0.9 and 1, as shown in Figure 10a.
This indicates that the proportion of small fluctuations is small.For the tunnel convergence rate measurement data series, the Δ  for DK115 + 280 is Δ  0.3663, which is the maximum value.This indicates that the proportion of small fluctuations is larger.The Δ  for DK115 + 261 is Δ  0.4360, which is the minimum value.This indicates that the proportion of small fluctuations is small.The findings suggest that a larger Δ  corresponds to small data changes, indicating that tunnel deformation is less affected by external influences such as excavation.Overall, the Δ  in the measured data series of the tunnel settlement rate and tunnel convergence rate of the same section shows the same directional change relationship, as shown in Figure 10b.

Tunnel Deformation Warning Level Classification Criteria
Based on the research results in the literature [31,61,62], the tunnel displacement monitoring criterion is constructed through the ∆ and ∆  parameters derived from the tunnel displacement monitoring data in order to realize the tunnel displacement warning level classification.The specific criteria are set as shown in Table 3.  Compare the proportions of large and small fluctuations (∆ f (α)).For the tunnel settlement rate measurement data series, the ∆ f (α) in DK115 + 280 is ∆ f (α) = −0.2251,which is the maximum value.This indicates that the proportion of small fluctuations is larger.The ∆ f (α) for DK115 + 261 is ∆ f (α) = −0.4415,which is the minimum value.This indicates that the proportion of small fluctuations is small.For the tunnel convergence rate measurement data series, the ∆ f (α) for DK115 + 280 is ∆ f (α) = 0.3663, which is the maximum value.This indicates that the proportion of small fluctuations is larger.The ∆ f (α) for DK115 + 261 is ∆ f (α) = −0.4360,which is the minimum value.This indicates that the proportion of small fluctuations is small.The findings suggest that a larger ∆ f (α) corresponds to small data changes, indicating that tunnel deformation is less affected by external influences such as excavation.Overall, the ∆ f (α) in the measured data series of the tunnel settlement rate and tunnel convergence rate of the same section shows the same directional change relationship, as shown in Figure 10b.

Tunnel Deformation Warning Level Classification Criteria
Based on the research results in the literature [31,61,62], the tunnel displacement monitoring criterion is constructed through the ∆α and ∆ f (α) parameters derived from the tunnel displacement monitoring data in order to realize the tunnel displacement warning level classification.The specific criteria are set as shown in Table 3.The trends of the two indicators need to be satisfied at the same time.In case of inconsistency in warning levels, the final warning level is determined according to the most unfavorable principle.

Warning Classification of Tunnel Deformation
The tunnel displacement monitoring data in Figure 4 are used to obtain the required sets of ∆α and ∆ f (α) parameters.In order to realize the trend judgment of the two discriminant indicators, 100 groups of tunnel displacement monitoring data are divided into one group.Synthesizing the actual situation and monitoring data, a total of 12 groups are divided.Additionally, calculate ∆α and ∆ f (α) parameters for each group of tunnel displacement monitoring data, as shown in Figure 11.The trends of the two indicators need to be satisfied at the same time.In case of inconsistency in warning levels, the final warning level is determined according to the most unfavorable principle.

Warning Classification of Tunnel Deformation
The tunnel displacement monitoring data in Figure 4 are used to obtain the required sets of ∆ and ∆  parameters.In order to realize the trend judgment of the two discriminant indicators, 100 groups of tunnel displacement monitoring data are divided into one group.Synthesizing the actual situation and monitoring data, a total of 12 groups are divided.Additionally, calculate ∆ and ∆  parameters for each group of tunnel displacement monitoring data, as shown in Figure 11.The trends of the two discriminant indicators were judged using the M-K test to achieve the early warning grading of tunnel displacement.The results are analyzed as follows:  The analysis of ∆ index criterion results: Through calculations and statistical analysis, the results of the ∆ index criterion are obtained (see Table 4).In the crosssection settlement monitoring data, the Z values for DK115 + 261 and DK115 + 270 fall within the range [−2.32, 2.32], which is a steady trend.The Z value for DK115 + 280 is less than −2.32, which is a decreasing trend.The Z value for DK115 + 290 is more than 2.32, which is an increasing trend.In the cross-section convergence monitoring data, the Z value for DK115 + 280 falls within the range [−2.32, 2.32], The trends of the two discriminant indicators were judged using the M-K test to achieve the early warning grading of tunnel displacement.The results are analyzed as follows: • The analysis of ∆α index criterion results: Through calculations and statistical analysis, the results of the ∆α index criterion are obtained (see Table 4).In the cross-section settlement monitoring data, the Z values for DK115 + 261 and DK115 + 270 fall within the range [−2.32, 2.32], which is a steady trend.The Z value for DK115 + 280 is less than −2.32, which is a decreasing trend.The Z value for DK115 + 290 is more than 2.32, which is an increasing trend.In the cross-section convergence monitoring data, the Z value for DK115 + 280 falls within the range [−2.32, 2.32], which is a steady trend.The Z values for DK115 + 261 and DK115 + 290 are less than −2.32, which is a decreasing trend.The Z value for DK115 + 270 is more than 2.32, which is an increasing trend.• The analysis of ∆ f (α) index criterion results: Through statistical calculations, the results of ∆ f (α) index criterion are obtained (see Table 5).In the section settlement monitoring data, the Z value for DK115 + 261 falls between the ranges [−2.In order to select the tunnel settlement prediction model, the performance of the PSO-LSTM and Back Propagation Neural Network (BP-ANN) prediction models is compared and analyzed.The monitoring data of section DK115 + 261 is used as the sample data to establish the settlement prediction model of PSO-LSTM and BP-ANN.80% of the data is used as the sample training set, and 20% of the data is allocated to the test set.The delay step of prediction model is set to 6, and the prediction is carried out across 1 time point.The adaptation change curve of PSO-LSTM and BP-ANN are shown in Figure 12.The tunnel section settlement is predicted, and the prediction results of the two models are compared.The tunnel settlement prediction results for different models are shown in Figure 13.Further, the performance of each model is evaluated using the model evaluation index.The calculation results of the evaluation indexes for the model test set are shown in Table 7. PSO-LSTM predicts significantly better than BP-ANN.Consequently, PSO-LSTM was selected for predicting tunnel settlement data.To make the forecast results more concise, subsequent forecasts will be conducted daily.A program was written using Python software (3.9.5) to optimize the parameters of The tunnel section settlement is predicted, and the prediction results of the two models are compared.The tunnel settlement prediction results for different models are shown in Figure 13.Further, the performance of each model is evaluated using the model evaluation index.The calculation results of the evaluation indexes for the model test set are shown in Table 7. PSO-LSTM predicts significantly better than BP-ANN.Consequently, PSO-LSTM was selected for predicting tunnel settlement data.The tunnel section settlement is predicted, and the prediction results of the two models are compared.The tunnel settlement prediction results for different models are shown in Figure 13.Further, the performance of each model is evaluated using the model evaluation index.The calculation results of the evaluation indexes for the model test set are shown in Table 7. PSO-LSTM predicts significantly better than BP-ANN.Consequently, PSO-LSTM was selected for predicting tunnel settlement data.To make the forecast results more concise, subsequent forecasts will be conducted daily.A program was written using Python software (3.9.5) to optimize the parameters of

Optimization of PSO-LSTM Parameters
To make the forecast results more concise, subsequent forecasts will be conducted daily.A program was written using Python software (3.9.5) to optimize the parameters of the LSTM prediction model using the PSO algorithm program.The parameters in PSO are set and a range is delineated for the optimization parameters in LSTM: the number of neurons ranges from 10 to 200, the number of training rounds ranges from 10 to 100, and the batch size ranges from 1 to 10. Initialize each parameter of the PSO-LSTM prediction model.Since the total amount of monitoring data is different for each section, after each set of data is input into the PSO program and finishes running, the optimal number of neurons, the number of training rounds, and the batch size of the model are obtained for the 4 groups.The details are shown in Table 8.In this study, the PSO-LATM prediction model was constructed using the Python software platform to complete the prediction analysis of cumulative tunnel settlement.The parameters for the LSTM neural network were calculated by PSO.Eighty percent of the data from each section was used as a sample training set and 20% as a test set.R 2 , MAE, and RMSE were used as evaluation indexes for predicting accuracy.Tunnel settlement was predicted for the next 5 days, and the prediction results are shown in Table 9.For the test set results of settlement prediction for four cross-sections, the test set R 2 values for the prediction model are 0.96, 0.97, 0.98, and 0.98.The MAE values are 0.28, 0.10, 0.04, and 0.02.The RMSE values are 0.06, 0.01, 0.04, and 0.02.Among these, the closer the R 2 values for prediction result converge to 1, the smaller the MAE and RMSE become, and the higher the model prediction accuracy is.Therefore, the settlement prediction results of the PSO-LSTM settlement prediction model for the four sections are analyzed.The model's predictions for the tunnel section settlement are in good agreement with the measured data.
Analyzing the prediction results, the future settlements of the four sections generally show a slowly increasing trend.The accuracy of tunnel warning classification is further verified, and, based on the evaluation indexes, the PSO-LSTM model demonstrates high accuracy in predicting tunnel settlement.

Conclusions
In this study, the multifractal theory was utilized to the analyze tunnel deformation, and the sliding time window was employed to enhance the segmentation of traditional multifractal subintervals.Furthermore, the M-K analysis method was utilized to ascertain the multifractal feature trends and establish the warning level for the tunnel section.Subsequently, the accuracy of the warning level was validated by predicting tunnel settlement using the PSO-LSTM prediction model.The following conclusions were drawn: 1.
The tunnel settlement and convergence rates of the four sections exhibit distinct fractal sequence characteristics.The width of the multifractal spectra in the measured data series of tunnel settlement rate and tunnel convergence rate within the same section shows an inverse relationship, with the sum of the two remaining stable between 0.9 and 1.Additionally, the proportion of size fluctuations in the measured data series of the tunnel settlement rate and the tunnel convergence rate within the same section demonstrate a consistent trend.

2.
The analysis of tunnel settlement prediction indicates that the PSO-LSTM prediction model delivers superior predictive performance and stability in tunnel settlement forecasts.

3.
A comprehensive analysis of the tunnel warning level and tunnel settlement prediction results reveals a class II tunnel deformation warning level, which aligns with the actual tunnel conditions.This approach, leveraging quantitative data as a reference, enables a more precise determination of the tunnel warning level.

Figure 1 .
Figure 1.Flow chart of tunnel deformation analysis and early warning prediction.

Figure 1 .
Figure 1.Flow chart of tunnel deformation analysis and early warning prediction.

Figure 3 .
Figure 3. Tunnel deformation monitoring and measurement: (a) tunnel monitoring point layout, (b) manual monitoring and measurement.

Figure 3 .
Figure 3. Tunnel deformation monitoring and measurement: (a) tunnel monitoring point layout, (b) manual monitoring and measurement.

Figure 3 .
Figure 3. Tunnel deformation monitoring and measurement: (a) tunnel monitoring point layout, (b) manual monitoring and measurement.

Figure 6 .
Figure 6.Trend plot of q-order fluctuation function    double logarithmic fit settlement rate of section DK115 + 261, (b) convergence rate of section DK115 + 261.

Figure 11 .
Figure 11.Calculated values of ∆ and ∆  parameters: (a) ∆ parameter values for tunnel settlement, (b) ∆  parameter values for tunnel settlement, (c) ∆ parameter values for tunnel convergence, (d) ∆  parameter values for tunnel convergence.

Table 1 .
Descriptive statistics of monitoring data.

Table 3 .
Tunnel warning level classification criteria.

Table 3 .
Tunnel warning level classification criteria.
The analysis of final warning results: On the basis of the results for the ∆α indicator criterion and ∆ f (α) indicator criterion, the final warning results of the four monitoring cross-sections are analyzed (see Table6).In the section settlement, the warning level for DK115 + 261 is level III, and all other sections are at level II.In section convergence, the warning level of DK115 + 280 is grade III, and all other sections are grade II.Therefore, according to the most unfavorable principle of synthesis, the final warning level of the four sections are all level II.That is, monitoring and measurement should be strengthened, and corresponding engineering countermeasures should be taken if necessary.

Table 4 .
Results of the ∆α indicator criterion.

Table 6 .
Final warning results of tunnel displacement.

Table 7 .
Results of calculated performance evaluation for the model.

Table 7 .
Results of calculated performance evaluation for the model.

Table 7 .
Results of calculated performance evaluation for the model.

Table 9 .
Predicted cumulative tunnel settlement values.