Abstract

Based on the actual engineering geological survey, the mechanical response mechanism and sliding surface characteristics of the interbedded slope containing six layers of soft and hard rock during excavation under heavy rain conditions are studied by numerical simulation, and the composite sliding surface stability algorithm of the layered slope is established by using the penalty function. The research shows that, during the excavation, the distribution characteristics of stress, strain, plastic zone, and horizontal displacement along the bedding plane are obvious; the stress is concentrated in the hard rock layer; however, the strain and plastic zone are concentrated in the soft rock layer. When the ratio of the maximum vertical displacement to the maximum horizontal displacement is less than 1.85, the slope will slide. According to the relative position of the soft rock layer and the toe of the slope, the sliding surface of soft and hard rock interbedded slope can be divided into two types, namely, circular-straight sliding surface and circular-straight-circular sliding surface. The verification shows that the composite sliding surface algorithm objectively reflects the destruction of the slope by bedding sliding and can reasonably evaluate the safety of soft and hard rock interbedded slope.

1. Introduction

The geological structure and mechanical properties of soft and hard rock interbedded slope are very complex, and the stability problem during excavation is prominent. In order to prevent the instability of this kind of slope excavation, the mechanical response law and stability calculation method during the excavation process should be mastered. In recent years, scholars have paid more and more attention to layered slopes. In the aspect of bedding slope stability, Cao et al. revealed the influence of various human factors on the three-dimensional stability of soft rock slopes with weak bedding during the excavation process [1]. Xue-chi et al. studied the deformation law and stability of bedding rock slope under secondary excavation [2]. Li et al. revealed the control effect of double weak interlayers on slope stability [3]. Zhang et al. studied the deformation mechanism and stability of bedding soft rock slopes [4]. Liu et al. proposed that the stability of rock slopes can be evaluated by analyzing the plane properties of in situ bedding [5]. In terms of deformation and failure of bedding slopes, Fan et al. proposed that the failure modes of bedding slopes are mainly vertical tension cracks at the rear edge of the slope, bedding sliding of rock layers along the interlayer, and rock block caving at the top of the slope [6]. Zhang et al. analyzed the deformation and failure mechanism of bedding slopes with gentle dip angles [7]. Xi et al. studied the failure characteristics of bedding rock slopes with multiple weak interlayers [8]. Zhou and Tang studied the sudden slip failure mechanism of bedding slopes based on the energy principle [9]. Li and Cheng studied the creep characteristics of bedding landslides [10]. In terms of the sliding surface of bedding slopes, Zhou et al. found that the strongly weathered soft and hard rock interbedded slopes form a fold-line potential sliding surface from joints, fissures, and rock layers, and there is more than one potential sliding surface [11]. Ding et al. studied the deformation and failure process of typical bedding mining slopes and the evolution law of rock and soil properties on the sliding surface [12]. Liu et al. used a new method to predict the failure surface of bedding rock slopes [13]. In addition, some scholars have studied the failure mode of bedding rock landslides, such as Mao et al. found that the weak structure in bedding rock slopes plays a major role in the stability of the slope [14]. Ma and Zhang proposed a sliding-bending-pull-crack composite bedding landslide model and pointed out that the deformation of the slope is obviously controlled by the weak interlayer and rock mass structure [15]. Wang et al. proposed the sliding failure-shear layer failure-tensile failure mode of weak sandwich slopes [16]. Lai et al. proposed a sliding and traction-slip composite sliding mode [17].

At present, most of the relevant research studies focus on the stability, deformation and failure modes, dynamic tests, shear resistance of structural planes, and slope reinforcement of bedding slopes. The mechanical response mechanism and the stability algorithm of bedding slope excavation need further study. In order to explore the mechanical response and the stability algorithm during the excavation process of the soft and hard rock interbedded slope, the numerical simulation of the actual project is carried out in this study to study the stress, strain, plastic zone, displacement, and other mechanics of the excavation under heavy rain conditions. According to the characteristics of the sliding surface of the layered slope, the slope stability algorithm is established, and the algorithm is verified by an example.

2. Excavation Mechanical Response Analysis

2.1. Engineering Overview and Numerical Model

The deep excavation cutting slope of a highway in Guizhou, China, is a typical soft and hard rock interbedded slope. The slope has six soft and hard rock layers with an inclination angle of approximately 15°. The excavation is divided into seven levels, and the maximum excavation depth is 53.6 m. The geological conditions of the excavation area are complex, and the slope sliding trend along the bed is obvious. This area belongs to the tectonic erosion of the landform of low mountains and hills, and the cut slope to be excavated is located on the mountain slope. The local rainfall is abundant, with a total rainfall of 600–700 mm in summer, and floods often occur. The strata from top to bottom in the area discovered by engineering geological survey and drilling are as follows:① Silty clay (Q4dl+el): it belongs to grade II ordinary soil, and the maximum thickness revealed by drilling is 5.93 m② Strongly weathered sandstone (J1-2zl): it belongs to grade IV soft stone, and the maximum thickness revealed by drilling is 13.32 m③ Moderately weathered sandstone (J1-2zl): it belongs to V subhard stone, and the maximum thickness revealed by drilling is 7.98 m④ Strongly weathered mudstone (J1-2zl): it belongs to III hard soil, and the maximum thickness exposed by drilling is 8.6 m⑤ Medium-weathered sandstone (J1-2zl): it belongs to V subhard stone, and the maximum thickness revealed by drilling is 6.65 m⑥ Strongly weathered mudstone (J1-2zl): it belongs to grade III hard soil, and the maximum thickness exposed by drilling is 8.45 m⑦ Medium-weathered sandstone (J1-2zl): it belongs to grade V subhard stone, and the maximum thickness is 9.07 m revealed by drilling (not exposed)

The main unfavorable geology in the area is bedding slope containing strongly weathered rock stratum, which is prone to bedding sliding under excavation disturbance. This area belongs to the area below seismic intensity VI, the peak value of ground motion acceleration is 0.05 g, and the general characteristic period of ground motion response is 0.35 s. Based on the previous background analysis, it can be seen that the stability of the slope during the excavation process is easily affected by the heavy rain weather and the occurrence and lithology of the rock formation and is less affected by the earthquake factors, so the heavy rain conditions are mainly considered in the numerical simulation in this study. According to the field geological exploration data and experimental data, the selected geotechnical simulation parameters under heavy rain conditions are shown in Table 1.

According to the actual size of the slope, a 150 m × 70 m numerical model is established, as shown in Figure 1, and the CPE4 element type is used for mesh division. Abaqus is used for numerical simulation, the geotechnical parameters under heavy rain conditions in Table 1 are input into the model, the excavation data of slope under heavy rain conditions are obtained through numerical simulation, and the mechanical response analysis is carried out on this basis.

2.2. Stress Response Analysis

It can be seen from Figure 2 of the distribution of shear stress after excavation of slopes at all levels that the shear stress is always in a state of constant adjustment, and its distribution is closely related to the position of each rock layer and the strength of the rock mass and is mostly concentrated in the relatively high-strength moderately weathered sandstone layers. Therefore, the shear stress has obvious band distribution characteristics along the rock layer; in the later stage of excavation, the concentration of shear stress near the toe of the slope on both sides is gradually serious.

In Abaqus, Sxy and Syx are the shear stresses in the XY plane pointing to the Y-axis and X-axis directions, respectively, with the Y-axis direction as a positive unit, kPa. Figure 3 shows the maximum shear stress in the slope body during excavation of slopes at all levels. Comparing the Sxy and Syx curves, it is found that the average shear stress Sxy is about 1.85 times that of Syx during excavation of slopes at all levels, and the former is 95 kPa larger than the latter on average. During the excavation process, the sum of the absolute values of Sxy and Syx was in an increasing trend before the first-level excavation, with a cumulative increase of 51.31%. After the second-level excavation, a single increase of 31.53% was achieved, and it reached the maximum value of 463 kPa for excavation at all levels. However, after the first-level excavation, the shear stress dropped sharply by 22.03%, which indicated that the shear failure of the hard rock layer in the slope body occurred at this time, and most of the shear stress was suddenly released.

2.3. Strain and Plastic Zone Analysis

It can be seen from Figure 4 of the maximum strain cloud after excavation of the slopes at all levels that the strain during the excavation process is mainly distributed in the soft rock layer with low strength, forming a “strain zone” distributed along the soft rock layer, and the maximum value appears in the soft rock layer. The two adjacent “strain zones” have a tendency to connect up and down, and the two “strain zones” are completely integrated into one after the first-level excavation. The “strain zone” on the side adjacent to the slope surface gradually extended obliquely downward along the soft rock layer and extended to the slope surface after the fourth to the first stage of excavation.

At the same time, the distribution of plastic zone obtained by numerical simulation is basically the same as the distribution of strain; after excavation of slopes at all levels, the plastic zone and the maximum strain zone are mainly distributed along the soft rock layer. As shown in Figure 5, the maximum strain value shows a U-shaped change trend; after the excavation of the first three levels and the last four levels, it first decreased and then increased; especially, after the first level of excavation, the single increase rate was as high as 216.20%. The maximum strain value of the first three levels slope gradually decreased because there was a silty clay layer on the upper part of the slope, and the strain was large in the early stage of excavation, and then, the strain gradually decreased with the enhancement of unloading. The maximum strain value of the last four-level slope increases step by step because the strain gradually accumulates in the local soft rock layer under the combined action of structural change, unloading, and stress redistribution. After the first level of excavation, the toe of the slope is just located in the strongly weathered mudstone layer; under the combined action of the weight of the slope body, the free surface, and the soft rock layer, the strain rapidly accumulates and the maximum strain value increases by 216.20%, indicating that the plastic failure is developing rapidly, and the slope has serious stability problems.

2.4. Displacement Response Analysis

It can be seen from the cloud (Figure 6) of the horizontal and vertical displacement distribution that the horizontal displacement is mainly manifested as the movement of the layered slope to the outside of the slope along the free side, and the horizontal displacement of the intersection area between the slope toe and the soft rock layer is the most serious with the progress of excavation and forms a layered sliding area with the toe of the slope as the front end of the sliding.

Under the strong unloading action of the slope body, the vertical displacement is divided into the downward settlement of the slope top and the upward uplift of the slope bottom. The influence area of vertical displacement is very wide, and almost, the whole slope moves up or down in different degrees.

With the excavation at all levels, the subsidence area is gradually concentrated towards the top of the slope near the mountain body, and the uplift area is gradually concentrated towards the subgrade at the foot of the slope. The distribution of horizontal displacement along the rock layer is obvious and gradually decreases from the toe of the slope to the top of the slope, while the vertical displacement is mainly related to the depth of the rock mass and the progress of excavation.

U1 and U2, respectively, represent the maximum values of horizontal displacement and vertical displacement in the model in mm. It can be seen from the curve of the maximum displacement value after excavation of each level of slope shown in Figure 7 that the maximum displacement values in the horizontal and vertical directions both decrease step by step with the excavation before the fifth level and increase with the excavation after the fifth level. The excavation increased step by step; especially, after the first-level excavation, the increase of U1 and U2 reached 130.34% and 45.41%, respectively; the former was about 2.87 times the latter, indicating that, after the first-level excavation, the amount of displacement in the horizontal direction increases rapidly.

Comparing the U1 and U2 curves in Figure 7, it can be seen that the maximum displacement curves in the two directions are almost parallel, indicating that the two displacement changes are not isolated and have inherent coordination. The increase of the horizontal displacement is about 6.30 times that of the vertical displacement, but the vertical displacement is 5.29 mm larger than the horizontal displacement on average, indicating that the increase of the displacement in the slope is mainly the horizontal displacement, and the displacement value is mainly the vertical displacement. The ratio of U2/U1 is parabolic, and U2/U1 gradually decreases with the fifth stage of excavation, from a maximum of 5.99 to 1.85, indicating that the displacement values in the two directions are accelerating and approaching, and the combined displacement of the two must also be acceleration increases, predicting that the slope is about to slip.

In summary, the displacement increase of the layered slope is mainly horizontal displacement, the horizontal displacement has obvious distribution characteristics along the rock layer, and the displacement value is mainly vertical displacement; the change curves of vertical displacement and horizontal displacement are almost parallel. When the ratio of U2/U1 is less than 1.85, the risk of slippage of the slope is higher.

2.5. Sliding Surface Feature Analysis

Since the disturbance and bedding sliding trend caused by the excavation of the first level to fourth level slopes to the soft rock is obviously larger than that of the other levels, therefore, the sliding surface characteristics of the first level to fourth level slopes are mainly analyzed.

When the strength of slopes at all levels is reduced to the limit state, the sliding surface has obvious arc-shaped sliding and bedding sliding. When the slope toe is above the soft rock layer, the sliding surface will form an arc-shaped sliding surface between the bottom layer of the soft rock layer and the slope toe, as shown in Figures 8(a)–8(d). When the slope toe passes through or under the soft rock layer, the slope body will slide out directly along the bottom layer of the soft rock layer, as shown in Figures 8(b) and 8(c).

It can be seen from this that the sliding surface of the bedding slope between soft and hard rocks is complex. According to the relative position of the soft rock layer and the toe of the slope, the sliding surface of soft and hard rock interbedded slope can be divided into two types: arc-shaped—bedding sliding surface and arc-shaped—bedding—arc-shaped sliding surface. If the sliding surface is assumed to be curved surface, it is not consistent with the actual situation.

2.6. Mechanical Response Analysis of Strength Reduction Process

It can be seen from Figure 9 that, in the process of slope strength reduction at all levels, the responses of stress, strain, and displacement show a certain regularity.

In terms of maximum shear stress, Sxy is generally larger than Syx, and Sxy first rises and then falls. It can be seen that the maximum shear stress Sxy has an obvious process of first accumulation and then release, and the shear stress is continuously transferred from the soft rock layer to the hard layer. When the shear stress in the hard rock layer accumulates to its shear strength, shear failure occurs, resulting in the release of shear stress.

The change curves of the maximum strain E and the maximum displacements U1 and U2 are basically the same; in the early stage of reduction, they remain basically unchanged or slightly increase, and they suddenly increase rapidly when approaching the final value of reduction. The sudden accelerated development of strain and displacement indicates that a large mechanical change has occurred inside the slope body. Combined with the deformation cloud map of the numerical software, it can be seen that the slope is in a state of imminent sliding at this time. At the same time, comparing the sudden release of the maximum shear stress with the reduction time points of the maximum strain and displacement mutation, it can be seen that the former generally occurs earlier than the latter, so the sudden release of the maximum shear stress does not indicate that the slope is in the limit stable state, while the accelerated changes in maximum strain and displacement reveal that the slope is about to fail.

3. Research on Stability Algorithms

3.1. Slope Stability Algorithm for Conventional Arc-Shaped Sliding Surface

When the slope only bears gravity, its stress field can be obtained by using the elastic solution of elastic wedge. If there is no special geological structure in the slope, its sliding surface is generally arc-shaped.

As shown in Figure 10, if there is a sliding surface SAB passing through point A at the toe of the slope with a height of h, we take the coordinate system xoy, and the projection of SAB on the plane coordinate system xoy is LAB. Let the expression of LAB be y = f(x), and the rock and soil parameters of the slope are c, φ, and γ, respectively. The physical force component X = 0, Y = γ, the angle between the slope surface and the x-axis is θ, the coordinates of the toe of the slope are A(, h), and the coordinates of the intersection point between the sliding surface and the top of the slope are B(xb, 0), where , h, and θ are known quantities.

Referring to the method of Gu and Han [18], it is assumed that the slope stress function ψ is a pure cubic expression of x and y, namely,

From the elastic mechanics, it can be known that the expression of the stress component at a point in the slope is

We substitute equation (1) into equation (2) to obtain

The stress components satisfy the following boundary conditions on the boundary

Among them, there is no stress on the horizontal boundary (y = 0), i.e., , and the direction cosine of its outer normal is ; at the same time, there is no stress on the slope boundary , i.e., , and the direction cosine of its outer normal is , so, it can be obtained that

We substitute (5) into equation (3) to obtainwhere , , and are the normal stress and shear stress at a point in the slope, respectively.

Regarding the determination of the sliding surface, different methods have different assumptions, and the obtained sliding surface shapes are also different. For example, the simplified Bishop method assumes the sliding surface to be an arc surface. The mechanical analytical rule assumes that the expression of each sliding surface is a quadratic polynomial. The critical sliding surface obtained by the finite element method by connecting the failure points with local stability coefficient less than 1 is also close to the arc surface.

As shown in Figure 11, there is no layered structure between soft and hard rock in the slope. At this time, the curve segment LAB can be completely regarded as an arc segment, and the normal stress and tangential shear stress at any point on LAB can be expressed aswhere the expressions of l and m are

On the sliding plane, and satisfy , substitute equation (8) and into equation (7), and obtain

The Mohr−Coulomb constitutive model is used for slope rock mass, and the stability coefficient of slope mass can be expressed in the following equation:where and are the cohesion force and internal friction angle of the rock and soil, respectively, and the numerator and denominator represent the total antisliding force and total sliding force of the slope on the sliding surface SAB under the limit equilibrium state, respectively.

We substitute equation (9) into equations (10) and (11) to obtain the following:

We substitute equation (6) into (11) to obtain

We substitute into (12) to obtain

Equation (13) is the calculation formula of stability coefficient of slope with conventional arc-shaped sliding surface, where , , and are known, and is obtained. Then, the corresponding stability coefficient can be obtained.

3.2. Stability Algorithm of Soft and Hard Rock Interbedded Slope

Based on the analysis of sliding surface in practical engineering, it can be known that the sliding surface of soft and hard rock interbedded slope cannot simply be assumed as a curved surface. Therefore, combined with the engineering practice, the sliding surface of this kind of slope is assumed to be two combination types of arc-shaped—bedding sliding surface (I) and arc-shaped—bedding—arc-shaped sliding surface (II), as shown in Figure 11. The coordinates of the intersection point between the sliding plane and the rock surface or slope surface are set as . For the convenience of calculation, the coordinates of point A are set as , and the coordinates of the boundary points between the bedding plane and the arc plane are C and D.

The soft and hard rock interbedded slope is heterogeneous, and the parameters , , and in adjacent rock strata differ greatly. Therefore, equation (13) needs to be solved and accumulated according to lithology stratification. Therefore, after the transformation of equation (13), the sliding surfaces of the previous two combination types (I, II) can be expressed bywhere n is the number of soft and hard rock strata in the sliding zone, n = 2, 3, ... , , , and , and f(x) is the expression of composite sliding surface.

Equation (14) is the calculation formula of stability coefficient of soft and hard rock interbedded slope, and the key to solve the formula lies in the determination of f(x). If an expression of f(x) can be found so that Fs(x) in equation (14) achieves the minimum, then the minimum value is the slope stability coefficient, and f(x), in this case, is the composite sliding surface expression. If the piecewise function f(x) is calculated in strict accordance with the arc plane or along plane, it will be very tedious. Therefore, f(x) is approximately fitted into a multiple continuous function expression.

Taking the sliding surface after excavation of the first stage slope obtained by numerical simulation as an example, the curves of the sliding surface were fitted with second-, third-, fourth-, and fifth-order polynomial, respectively. By seeing Figure 12, it is found that the higher the order of fitting, the closer the curve is to the sliding surface projection obtained by numerical simulation. When the order of fitting is 4 or above, a relatively ideal effect can be achieved. However, if the order of fitting is too high, it is not conducive to calculation, so this study adopts the fourth-order polynomial for fitting. At the same time, it can be seen that when f(x) is a quadratic polynomial, the fitting curve differs greatly from the projection of the sliding plane, which also confirms that f(x) cannot be assumed as a quadratic polynomial only.

Based on the previous analysis results, can be assumed. The key of the problem is to obtain a set of (a, b, c, d, e) in f(x) so that the minimum value of Fs(x) can be obtained. In fact, it becomes a constrained optimization problem, for which the penalty function method is introduced to solve it.

Referring to the practice of Aizhong and Jia [19], the projection LAB of slope sliding surface SAB is isometrically divided into m segments along the vertical direction, as shown in Figure 13. Assuming that the expression of each segment is a quartic polynomial, the curve segment LAB can be expressed aswhere ai, bi, ci, di, ei (i = 1, 2, …, ) is the coefficient to be determined.

By ,

According to the continuity and smoothness of f(x), the following can be obtained:

According to equations (15) and (17), it can be concluded that

It can be obtained from the known quantities in coordinates of points A and B that

Combined with equations (16), (18), and (19), we can obtain the following:

In order to ensure that the sliding surface is always inside the slope body during the optimization process, constraint conditions are added to the ith subsection point as follows:where , xi, and are the abscissa values of the sliding surface projection LAB and slope surface projection LAO on the same segment line, respectively.

According to the types of the previous problems and characteristics of penalty functions, the mixed penalty function was selected for optimization, and the stability algorithm based on the composite sliding surface search model was established as follows:where is the objective function, s.t. is the constraint condition, is the mixed augmented objective function, is the penalty parameter, and . The calculation process can be realized by MATLAB programming.

4. Example Verification of Stability Algorithm

4.1. Stability Algorithm Calculation Based on the Composite Sliding Surface Search Model

The previous composite sliding surface model algorithm is used to focus on the calculation of the sliding surface of the slope from the first stage to the fourth stage with obvious bedding sliding, and the results are shown in Figure 14. The expressions and of potential sliding surface after slope excavation are obtained by calculation; the results are shown in Table 2.

4.2. Algorithm Verification Based on the Composite Sliding Surface Search Model

Firstly, we calculate the numerical simulation results of slope stability coefficient, and at the same time, we calculate the slope stability coefficient that only considers the sliding surface as a quadratic arc-shaped surface under the conventional algorithm ( is assumed to be a second order polynomial) to test the calculation effect of the slope stability algorithm based on the composite sliding surface search model.

In the Abaqus strength reduction process, when the strain does not converge and the calculation is forced to stop, it indicates that the slope is greatly deformed and unstable; at this time, the strength reduction coefficient is the slope stability coefficient.

Firstly, the calculated data in Abaqus are processed, and the relationship between strain and reduction coefficient during excavation is shown in Figure 15.

The ratio K of slope strain increment and reduction coefficient increment corresponding to time t in the calculation of strength reduction is shown in the following equation:where and are, respectively, strain variables and strain increment at time t and and are, respectively, reduction coefficient and reduction coefficient increment at time t. In (23), assuming that the maximum value is obtained when , the corresponding reduction coefficient is the slope stability coefficient, namely,where is slope stability coefficient. Equation (24) is used to calculate the stability coefficient on the curve in Figure 15, as shown in Figure 16.

In the calculation of slope stability coefficient under the arc-shaped sliding surface algorithm, is expressed in equation (13). In Figure 14, it is assumed that the expression of each segment is a quadratic polynomial; that is, and in equation (15). The slope stability coefficient under the quadratic arc-shaped sliding surface can be obtained by optimizing the mixed penalty function. The stability coefficients of the three methods are shown in Table 3.

As can be seen from Table 3 and the stability coefficient of slope at all levels, the result of the arc-shaped sliding surface algorithm is the largest, and the maximum deviation from the numerical simulation result is 14.42%. The maximum deviation between the composite sliding surface algorithm and the numerical simulation results is only 4.49%. Therefore, the composite sliding surface algorithm is closer to the numerical simulation results than the arc-shaped sliding surface algorithm.

Because the arc-shaped sliding surface algorithm only regards the sliding surface as the secondary arc surface and ignores the destruction of bedding sliding on the slope, the calculation results of slope stability coefficient are on the high side, which will overestimate its stability. However, the composite sliding surface algorithm has a relatively lower grade on slope stability, which is easy to draw people’s attention to slope safety. In summary, the composite sliding surface algorithm objectively reflects the destruction of bedding sliding on the slope and is more reasonable than the arc-shaped sliding surface algorithm in the safety assessment of such slopes.

5. Conclusion

(1)The stress, strain, plastic zone, and horizontal displacement of bedding slope have obvious distribution characteristics along the rock surface. The stress, strain, and plastic zone are concentrated in hard rock and soft rock, respectively. The increase of slope displacement is mainly horizontal displacement, and the larger displacement value is mainly vertical displacement.(2)The strain and displacement of bedding slope will accelerate when the slope is close to the ultimate stability state, and the sudden release of shear stress usually occurs earlier than the ultimate stability state.(3)According to the relative position of soft rock layer and slope foot, the sliding surface of bedding slope can be divided into two types, namely, arc-shaped—bedding sliding surface and arc-shaped—bedding—arc-shaped sliding surface.(4)The composite sliding surface algorithm objectively reflects the destruction of the slope by bedding sliding and can reasonably evaluate the safety of soft and hard rock interbedded slope.

Data Availability

The data used to support the findings of this study can be obtained from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (Grants 41472262 and 52009014), Innovation Research Group of Universities of Chongqing (Grant CXQT19021), key project of Chongqing Natural Science Foundation (Grant cstc2020jcyj-zdxmX0012), and the scientific personnel and technical staff at the State Key Laboratory of Mountain Bridge and Tunnel Engineering, Chongqing Jiaotong University, Chongqing.