Dynamic stability analysis method of anchored rocky slope considering seismic deterioration effect

The seismic deterioration effects of anchor cables and slope structural planes are often neglected in the dynamic stability analysis of anchored rocky slopes to the extent that the stability of slopes is overestimated. In this paper, a dynamic calculation method for anchored rocky slopes considering the seismic deterioration effect is established, and a stability evaluation method for anchored rocky slopes based on the Gaussian mixture model is proposed. The seismic deterioration effect on the stability of anchored rocky slopes is quantitatively analyzed with an engineering example, and the relationship between seismic intensity and the failure probability of slopes is clarified. The results show that compared with the calculation method without considering the seismic deterioration effect, the minimum safety factor and post-earthquake safety factor obtained by the proposed method in this paper are smaller. The number of seismic deteriorations of the slope is used as the number of components of the Gaussian mixture model to construct the failure probability model of the slope, which can accurately predict the failure probability of anchored rocky slopes. The research results significantly improve the accuracy of the stability calculation of anchored rocky slopes, which can be used to guide the seismic design and safety assessment of anchored rocky slopes.


Seismic deterioration effect
The safety factor of anchored rocky slopes is the ratio of anti-sliding force to sliding force.Under the action of an earthquake, the anti-sliding force provided by the anchor cable and the structural plane may be reduced to varying degrees due to the seismic deterioration of the slope.The main reasons are as follows.
(1) Slip deterioration effect of structural planes 11,16 : According to the Newmark displacement method, the slope slips when the seismic acceleration is greater than the yield acceleration of the slope.Slope slippage will inevitably cause abrasion of the structural plane, thereby reducing the shear strength and stiffness of the structural plane.
(2) Frictional attenuation effect of structural planes 5,15 : Under the seismic cyclic load, the sliding body and the bedrock will inevitably produce relative velocity during earthquakes.The relative velocity leads to a decrease in the friction coefficient of the structural plane, which reduces the strength of the structural plane.(3) Damage effect of prestressed anchor cables 17,20 : During earthquakes, damage effects may occur in anchor cables subjected to the tensile action of the sliding body, resulting in different degrees of deterioration of the axial force and stiffness of anchor cables.

Slip deterioration effect of structural planes
The slope slip process triggered by the earthquake can be regarded as the structural plane shear process under displacement-controlled loading, and the shear displacement is the slip distance of the slope, which can be obtained by the Newmark displacement method.Figure 1 is the shear stress path diagram of the structural plane under multiple slips of the slope.The stress path of the first slope slip process is O → A → B → C → D.Where OA is the rising section of shear stress, and point B is the peak point of shear stress.When the slip displacement exceeds the peak displacement u b , the shear stress begins to enter the BC descending section.Finally, the stress curve returns to point D when the displacement loading ends.On the whole, the shear stress path of the structural plane from the beginning to the end of the slip of the slope is basically the same as the stress path of the tangential loading-unloading of the structural plane.That is, there is an elastic deformation of the structural plane at the beginning of loading, and the tangential stress increases linearly with the displacement.As the loading displacement increases further, the shear stress decreases beyond the peak shear stress.At the end of the displacement loading, part of the elastic deformation of the structural plane is restored, and the shear stress is reduced to the static shear stress (shear stresses under static forces such as prestress and gravity).
The shear stress-shear displacement curves during the first slip of the slope usually show a significant peak.As the slip distance increases, the structural plane is continuously worn, and the shear stress-shear displacement curve of the structural plane no longer shows significant peaks during subsequent slips.
The peak shear stress of the structural plane during the first slip of the slope is: where ϕ b is the fundamental friction angle of a smooth structural plane, α k is the undulation angle of the struc- tural plane, and σ n is the normal stress exerted on the structural plane.
The peak shear stress of the structural plane that has been severely worn after several slips is: The degree of deterioration of the slope structural plane is closely related to the slope slip displacement.Gao et al. 16 proposed a method to fit the shear displacement-degradation coefficient with a negative exponential function, and the calculation formula of the shear strength degradation coefficient of the structural plane can be obtained by normalizing the shear displacement: where A , B and C are the coefficients to be determined, u rel is the one slip displacement calculated by the New- mark displacement method, and L is the length of the structural plane.
According to Plesha's study 6 , the calculation formula of the undulating angle degradation of the structural plane under cyclic shear action is as follows: (1) www.nature.com/scientificreports/where α k0 is the initial undulation angle of the structural plane, JRC is the roughness coefficient of the structural plane, JCS is the uniaxial compressive strength of rock, and W P is the plastic work.The formula for calculating the degradation of the undulation angle of the structural plane is organized as the product of the degradation coefficient of the undulation angle and the initial undulation angle.
The degradation coefficient is a function related to the plastic work, and its expression is as follows: According to Eq. ( 7), it is known that the degradation coefficient of the undulation angle of the structural plane has a negative exponential relationship with JRC(σ n /JCS) .Therefore, Eq. ( 3) can be further modified to obtain the Eq. ( 8) for the degradation coefficient of the undulation angle of the structural plane.
As a result, a degradation equation for the undulation angle of the structural plane related to the shear displacement can be obtained as follows: According to the theoretical derivation of Dong 46 , the correspondence between the undulation angle and the dilatancy angle of the structural plane is obtained as follows: Barton 47 conducted experiments on eight different rough structural planes and proposed a conversion formula for the dilatancy angle to the roughness coefficient of the structural plane.Equation ( 12) can be obtained by substituting Eqs. ( 9)-( 11) into Eq.(1): Therefore, the degradation coefficient of the undulation angle can be considered as the degradation coefficient of the roughness of the structural plane.
According to Wu's study 11 , it was found that the degradation law of tangential stiffness is similar to that of tangential strength of the structural plane.Therefore, the degradation formula of tangential stiffness of the structural plane can be expressed by negative exponential function as follows: Where k s0 is the initial tangential stiffness of the structural plane, and D , E and F are the coefficients to be determined.

Slip deterioration effect of structural planes
Wang 12 clarified that the dynamic friction coefficient is composed of the starting friction coefficient and the velocity-dependent function linearity through the sliding plane friction test of granite, and the expression is as follows: where f s is the starting friction coefficient, v s is the relative velocity of the upper block relative to the bedrock, and f µ (v s ) is a function that decreases as |v s | increases.
Based on Wang's research results, Ni 15 proposed the hypothesis that the peak shear strength of the structural plane has a negative exponential function relationship with the relative velocity, and gave the calculation formula of the relative velocity damage coefficient.
where γ r is the convergence value of the relative velocity attenuation, and a is a coefficient to be determined.Based on the calculation of Liu 5 , Ni 15 and Gao 16 , it's specified that γ r = 0.9 and a = 25 in this paper.
Therefore, considering the influence of the frictional attenuation effect and slip deterioration effect on the shear strength of the structural plane, the formula for calculating the shear strength of the structural plane can be obtained as follows: www.nature.com/scientificreports/Liu 5 transformed the formula for calculating the shear strength of the structural plane and deduced the formula for calculating the equivalent friction angle related to the undulation angle as follows: Substituting Eqs. ( 9)- (11) into Eq.( 17), the formula for calculating the equivalent friction angle with respect to the roughness coefficient of the structural plane is obtained as follows:

Damage effect of prestressed anchor cables
Anchor cables are flexible support structures that can only withstand tensile forces but cannot resist bending moments and shear forces.The failure modes of anchor cables are generally divided into brittle failures such as anchor head cracking and anchor pier collapse, and ductile failures such as anchor cable breakage.Among them, the maximum axial force is usually used as the failure index for brittle failure mode, while the elongation is used as the failure index for ductile failure due to the obvious yield stage of the anchor cable 17 .Based on the above considerations, the calculation model of the anchor cable for the brittle failure mode and the ductile failure mode is proposed in this paper.Figure 2 shows the p-s curve of the calculation model of the anchor cable.The brittle failure process of the anchor cable is composed of an elastic stage (AB) and a failure stage (BC), and the ductile failure process is composed of an elastic stage (Aʹ Bʹ), a plastic stage (Bʹ Cʹ) and a failure stage (Cʹ Dʹ).When the anchor cable is in the elastic phase (OA or OAʹ), the stiffness of the anchor cable is the initial stiffness, and the static axial force of the anchor cable is related to the slip displacement of the sliding body.If the slope does not slip, the static axial force is the prestress of the anchor cable.If the slope slips, the static axial force is the sum of the prestress of the anchor cable and the increase in the axial force of the anchor cable caused by the slope sliding.The dynamic axial force of the anchor cable is the sum of the static axial force and the increment of the axial force of the anchor cable caused by the earthquake.When the anchor cable is in the plastic stage (AʹBʹ), the stiffness of the anchor cable is 0, and the static axial force and dynamic axial force of the anchor cable are the axial force when the anchor cable yields.When the anchor cable is in the failure stage (BC or CʹDʹ), the stiffness of the anchor cable is 0, and the static axial force and dynamic axial force of the anchor cable are 0.
The calculation model of the anchor cable for the brittle failure mode is: where T p is the static axial force of the anchor cable, T 0 is the initial prestress of the anchor cable, T d is the dynamic axial force of the anchor cable, T ′ max is the maximum axial force of the brittle failure of the anchor cable, k f is the stiffness of the anchor cable, u p is the total slip displacement of the sliding body, u f is the total stretch of the anchor cable in the axial direction under the earthquake, and l is the length of free section of the anchor cable.
The calculation model of the anchor cable for the ductile failure mode is: where T max is the yield axial force of the anchor cable, and u max is the maximum stretch of the anchor cable.The total stretch of the anchor cable in the axial direction under the earthquake is where u s and u n are the displacement response of the sliding body in s and n directions during the earthquake.α is the inclination of the structural plane, and θ is the angle between the anchor cable and the horizontal plane.The stretches of the anchor cable in the s and n directions under the prestress are: The total slip displacement of the sliding body is the sum of each slip displacement.

Slope stability evaluation
In this section, the rocky slope reinforced by prestressed anchor cables is taken as the research object, considering the seismic deterioration effect of the structural plane and anchor cables during the earthquake, the time history of the safety factor of the slope during the earthquake is solved, and the stability evaluation method of the slope is established based on the Gaussian mixture model.

Dynamic calculation model of anchored rocky slope
The basic assumptions are inevitably made when the dynamic calculation model of anchored rocky slope is established as follows.
(1) The anchored rocky slope is simplified to a two-dimensional planar model in the calculation model; (2) The sliding body is an ideal rigid body that is homogeneous, continuous, and isotropic; (3) The influence of the self-weight of anchor cables is ignored; (4) Only the dynamic effect of horizontal seismic loads on the anchored rocky slope is considered.
According to the dynamic calculation model of rock mass proposed by Xue 48 , the structural plane of the slope can be regarded as a viscoelastic-plastic model.Anchor cables are usually considered as spring supports in the dynamic calculation model of the anchored slope 22 .The dynamic calculation model of the anchored rocky slope as shown in Fig. 3 can be established.
The Newmark displacement method is used as the calculation method for the plastic slip of slopes in the dynamic calculation model of the anchored rocky slopes 21 .The yield acceleration is the key to determining the plastic slip of the slope.The seismic force is assumed to be a horizontal static load, and the ultimate equilibrium method is used to calculate the seismic acceleration when the safety factor is 1, that is, the yield acceleration.It should be noted that with the real-time update of the deterioration parameters of the structural plane and anchor cables, the yield acceleration also changes dynamically.
Figure 4 shows the force diagram of the sliding body under earthquakes.Considering the force equilibrium condition of the sliding body, the normal force along the n direction is: where m is the mass of the sliding body, and a c is the yield acceleration.
The sliding force along the s direction is: The anti-sliding force and sliding force of the sliding body can be obtained by analyzing the force on the sliding body.According to the limit equilibrium method, the ratio of the anti-sliding force to the sliding force is the safety factor of the slope.( 22) where c is the cohesion of the structural plane.
The formula for calculating the safety factor of the slope can be obtained by substituting Eqs. ( 25) and ( 26) into Eq.( 27).
In the Newmark displacement method, the seismic acceleration with a safety factor of 1 is defined as the yield acceleration.Let F s = 1 , the expression for yield acceleration can be obtained by equation transformation as follows:

Seismic acceleration of the sliding body
The traditional Newmark displacement method takes the input seismic acceleration as the acceleration of the sliding body.However, according to Jia 21 , it was found that seismic waves passing through a structural plane change with the stiffness of the structural plane, but neglected the deterioration effect of the stiffness and the damping effect of the structural plane.In this paper, considering the stiffness degradation effect and the structural plane damping effect, the real seismic acceleration of the sliding body is obtained by establishing the dynamic balance equation in the horizontal direction and updating the parameters of the structural plane and anchor cables in real-time.
The dynamic calculation model of anchored rocky slope in the horizontal direction is shown in Fig. 5, and the dynamic equilibrium equation is established as follows: where ∂ 2 u g /∂t 2 is the seismic acceleration, c h is the damping of the structural plane in the horizontal direction, k h is the equivalent stiffness of the structural plane in the horizontal direction, u h is the displacement response in the horizontal direction, u hT is the static displacement in the horizontal direction under the prestress, and n is the number of anchor cables.
According to the principle of equal deformation energy, the equivalent spring stiffness of the horizontal direction of the structural plane can be obtained as follows:  where k n1 and k s1 are the stiffness of the structural plane in the n and s directions.
For single-degree-of-freedom systems, the Duhamel integral can be used to solve the dynamics equation 49 .The displacement and velocity at the initial moment are: where f (τ ) is the external load at time τ.
There are in the differential time interval where ξ and ω are the damping ratio and natural frequency without damping, and ω D is the natural frequency with damping.
The natural frequency with damping can be obtained as follows: The entire load history can be seen as consisting of a series of successive impulse loads, each of which produces the differential response shown in Eq. (33).The total response can be used as a superposition of all the differential responses resulting from the load time history.Integration of Eq. ( 33) can be obtained as follows: where There is some error in the trapezoidal method for solving the definite integral, especially in the early stages of the integration.The accuracy of the recursive method for solving the Duhamel integral is higher than that of the definite integral method, so the recursive method is used in this paper to solve the dynamic response of the slope.
The recursive method assumes that the response of the slope at time t is known and solves for the dynamic response of the slope at time t + Δt.When the time interval is Δt, the external load corresponding to time τ is The dynamic response of the slope at time t + Δt is a superposition of the following three cases: (1) Free vibration with initial condition u(0) = u t and v(0) = v t ; (2) Forced vibration with constant external load F t ; (3) Forced vibration with external load F(τ ) = F t + F t+�t −F t �t τ .The Duhamel integrals for the above three cases are superimposed to obtain the equation for the dynamic response of the slope at time t + Δt. (31)

The time history of the safety factor
When the seismic acceleration in the horizontal direction of the sliding body is less than the yield acceleration, the tangential plastic element is not triggered, and only the spring and damping play a role.The dynamic model of the anchored rocky slope is simplified by the centralized mass method 50 , and the dynamic response of the slope is decomposed in the s and n directions, as shown in Fig. 6.
The dynamic balance equations for the anchored rocky slope in the s and n directions are established as follows: where c s and c n are the damping of the structural plane in the s and n directions, k s and k n are the equivalent stiffness of the structural plane in the s and n directions, u s and u n are the static displacements in the s and n directions under the weight of the sliding body, u sT and u nT are the static displacements in the s and n direc- tions under the prestress.
The equivalent stiffness of the structural plane in the s and n directions can be obtained from the principle of equal deformation energy as follows: The tangential and normal forces acting on the structural plane during the earthquake can be obtained as follows: The seismic acceleration, velocity and displacement responses of the sliding body in the s and n directions at any t time can be obtained by iterative calculation of Eq. (36).Substituting the velocity response in the s direction at any time t into Eq.( 15) can obtain the corresponding relative velocity damage coefficient.Substituting the displacement response at any time t into Eq.(19) or Eq. ( 20) can obtain the corresponding dynamic axial force and stiffness of the anchor cable.Substituting the displacement response at any time t into Eqs.( 53) and ( 54), the corresponding shear force and normal force acting on the structural plane can be further obtained.
When the seismic acceleration in the horizontal direction of the sliding body is greater than the yield acceleration, the tangential plastic element of the structural plane is triggered, and the sliding body slides tangentially down along the structural plane with a slip acceleration.
The slip acceleration at time t can be defined as: (44) The slip acceleration at time t can be integrated to obtain the slip velocity of the sliding body at time t.
By integrating the slip velocity, the slip displacement of the sliding body at time t can be obtained as follows: The slip displacement of the sliding body is substituted into Eq.( 12) and Eq. ( 13) respectively to obtain the tangential strength and stiffness of the structural plane after slip deterioration.The slip displacement of the sliding body is substituted into Eq.(19) or Eq. ( 20) respectively to obtain the dynamic axial force and the stiffness of the anchor cable after damage.
Based on the limit equilibrium method [51][52][53] , the dynamic safety factor at time t can be obtained by substituting the tangential force and normal force acting on the structural plane at any time t, the degradation coefficient of the roughness at time t, and the relative velocity damage coefficient at time t into Eq.( 58).

Failure probability of anchored rocky slopes
Seismic loads are usually considered as random variables of time, so the dynamic safety factor of anchored rocky slopes can also be considered as a random function of time, and the stability of slopes during earthquakes can be evaluated by the probability evaluation method 5 .According to the central limit theorem, it is theoretically proved that the distribution of the dynamic safety factor of the slope obeys normal probability law 45 .However, as the seismic deterioration effect of the slope gradually appears, the nonlinearity of the Quantile-Quantum curve of the dynamic safety factor becomes more and more significant, indicating that the probability distribution of the dynamic safety factor deviates from the normal distribution.At this point, the Gaussian mixture model can better express the complex probability distribution of the dynamic safety factor.Gaussian mixture models are combinations of two or more Gaussian probability density functions and are very popular in density estimation and clustering [54][55][56] .Although in some cases the number of components in a Gaussian mixture model may be unlimited, in general, the number of components is limited to a finite 57 .The probability density function of a Gaussian mixture model can be written in the form of sum: where π is a parameter vector including the weight coefficients of the components of the Gaussian mixture model, µ is a parameter vector including the means of the components of the Gaussian mixture model, σ is a parameter vector including the variances of the components of the Gaussian mixture model, and k is the number of components of the Gaussian mixture model.
The seismic deterioration effect of the anchored rocky slope will cause shifts in the median of the safety factor time history of the slope, and the median value of the safety factor time history is closely related to the expectation of the components of the Gaussian mixture model.Therefore, when the Gaussian mixture model is used as the probability density function of the dynamic safety factor of the slope, the number of seismic deteriorations can be used as the number of components of the Gaussian mixture model.It should be noted that the number of seismic deteriorations is the sum of the number of slips of the sliding body and the number of anchor cable failures.
Maximum likelihood estimation is an efficient method to obtain probability distribution parameters by maximizing the logarithm of the likelihood function.The EM algorithm is commonly used to calculate maximum likelihood estimates in the presence of latent variables or missing data and can be used to determine the parameters of a Gaussian mixture model.The EM algorithm is an efficient iterative method where each iteration consists of two processes 57 : (1) E-step: Solve the expectation p Z|X, θ (t) .
(2) M-step: Solve the maximum and calculate the model parameter θ (t+1) for a new iteration.
In step E, given the current estimate of the parameter θ (t) , the conditional distribution of Z is determined as the proportional height of the normal density weighted by π: where θ k = (π k , µ k , σ k ) is a parameter vector including the weight, mean and variance of the Gaussian compo- nent of N ⊖ k .
In step M, the next estimate θ (t+1) can be determined by maximizing the conditional expectation Q in step E, and calculated by Eq. (64).
The expectation Q in Eq. (64) can be constructed as follows: The π k , µ k and σ k are separate linear terms, so they can be maximized independently.According to ∂Q ∂π k = 0 , ∂Q ∂µ k = 0 , ∂Q ∂σ k = 0 and s.t K k=1 π k = 1 , the next estimates can be obtained: Vol:.( 1234567890 www.nature.com/scientificreports/ The Technical Code for Building Slope Engineering (GB50330-2013) stipulates the minimum value of the stability safety factor of the slope of each safety grade under seismic conditions, and the required minimum values of the safety factor of the grade I, II and III slopes are 1.15, 1.1 and 1.05, respectively.In this evaluation method, the minimum safety factor value required by the code is used as the evaluation index of slope stability, and the probability of being less than the required minimum safety factor value is defined as the failure probability of the slope.

Stability evaluation steps for anchored rocky slopes
Based on the time history analysis method and probability analysis method, a dynamic stability analysis method of anchor rock slope considering the seismic degradation effect is proposed.The calculation steps are as follows.
(1) The quasi-static method is used to analyze the force of the sliding body, and the yield acceleration of the sliding body is obtained by Eq. ( 29). ( 2) The dynamic equation of the slope in the horizontal direction is established by Eq. (30).The seismic acceleration of the sliding body at the current time step is obtained by the Duhamel integral recursion method.By comparing the relationship between the yield acceleration and the seismic acceleration of the sliding body, it is determined whether the plastic slip of the slope is triggered.(3) When a c > a h , the sliding body does not slip.The dynamic equations of the slope in the s and n directions are established by Eqs. ( 49) and ( 50).The acceleration, velocity and displacement of the sliding body at the current time step are obtained by the Duhamel integral recursion method.Substituting the velocity into Eq.( 15) to obtain the relative velocity damage coefficient of the sliding body.Substituting the displacement into Eqs.(19) or (20) to obtain the static axial force, dynamic axial force and stiffness of the anchor cable.When a c < a h , the slip displacement of the sliding body at the current time step is calculated by Eq. (57).Substituting the slip displacement into Eqs.( 8) and ( 13) to obtain the degradation coefficient of the roughness and degradation stiffness of the structural plane.Substituting the slip displacement into Eqs.(19) or (20) to obtain the static axial force, dynamic axial force and stiffness of the anchor cable.(4) The results of the calculation in step (3) are used to update the current calculation parameters, and the updated calculation parameters are re-substituted into steps (1) to (3) to carry out the dynamic calculation for the next time step.This cycle is repeated until the end of the calculation.(5) The time history of the safety factor of the slope is obtained by substituting the tangential force and normal force of the structural plane, the relative velocity damage coefficient of the sliding body and the degradation coefficient of the roughness at each time step into Eq.( 58). ( 6) The number of seismic deteriorations of the anchored rocky slope is counted, and the number of seismic deteriorations is used as the number of components of the Gaussian mixture model to establish the probability density function of the Gaussian mixture model.(7) The appropriate slope stability evaluation index is selected, and the slope failure probability corresponding to this index is obtained from the cumulative distribution function of the anchored rocky slope.

The case study
The basic parameters of the slope In this paper, a specific engineering example is selected to quantitatively analyze the dynamic response and stability of anchored rocky slopes.As shown in Fig. 7, the slope height H = 14 m, the slope inclination angle β = 70°, and the structural plane inclination angle α = 45°.The rock mass density ρ = 2700 kg/m 3 , the structural plane initial roughness JRC = 12.3, the basic friction angle φ b = 25°, the cohesion c = 35 MPa, the rock strength JCS = 120 MPa, the damping ratio ξ = 0.15, the normal stiffness k n = 3 MPa, and the tangential stiffness k s = 1 MPa.The slope is arranged with six rows of prestressed anchor cables from top to bottom.The anchor cable type is 2φ15.2.The length of the free section l = 6 m, the anchor cable inclination angle θ = 10°, the initial prestress T f = 100 kN, and the maximum bearing tension force T max ʹ = 136.5 kN.The coefficients to be determined for the roughness and tangential stiffness degradation formulas for the structural plane obtained by least squares fitting are listed in Table 1.In this case, the input seismic wave is the ChiChi wave with a PGA of 0.8 g in the horizontal direction, and the time history curve and spectrum of the seismic wave after filtering and baseline correction are shown in Fig. 8. (66)

Influence of seismic deterioration effect
Figure 9 shows the comparison of the seismic acceleration time history of the sliding body with the input seismic acceleration time history.From Fig. 9, it can be seen that the waveform of the input seismic acceleration time history is basically similar to that of the seismic acceleration time history of the sliding body, but the fluctuation amplitude of the seismic acceleration time history of the sliding body is larger than that of the input seismic acceleration time history, which is 0.96 g.The difference in the amplitude between the seismic acceleration time history of the sliding body and the input seismic acceleration time history is mainly related to the frequency of the input seismic acceleration time history and the natural frequency of the slope.Failure of the anchor cable or slippage of the sliding body reduces the equivalent stiffness of the structural plane, which in turn reduces the natural frequency of the slope.This may result in the superior frequency of the input seismic acceleration time history being closer to the natural frequency of the slope, so that the slope produces a certain resonance effect, thereby increasing the amplitude of the seismic acceleration time history of the sliding body.In addition, due to the damping effect of the structural plane, the PGA of the sliding body is delayed with respect to the input PGA.
In the traditional Newmark displacement method, the yield acceleration is considered to be a constant value that does not change during the earthquake.In this paper, the slip deterioration of the structural plane and the damage of anchor cables are fully considered during the earthquake, and the yield acceleration of the sliding   body is updated in real-time to obtain a more reasonable slip displacement of the sliding body.As can be seen from Fig. 10, the yield acceleration of the sliding body decreases several times in different degrees.When the yield acceleration drops abruptly for the first time, the slip displacement of the sliding body is 0, indicating that the decrease in yield acceleration is not caused by the degradation of the structural plane but by the failure of anchor cables.Then, with the multiple slips of the sliding body, the yield acceleration gradually decreases, and the final yield acceleration tends to a constant.In general, the seismic deterioration effect will reduce the yield acceleration of the sliding body, and thus increase the slip displacement of the sliding body.
In order to clarify the influence of the relative velocity of the sliding body on the friction attenuation effect of the structural plane, the friction attenuation factor is used in this paper to describe the degree of friction attenuation of the structural plane.It is worth noting that the friction attenuation factor and the relative velocity damage coefficient are the sum of 1. Figure 11 shows the time history of the friction attenuation factor during the earthquake.The frictional attenuation effect occurs when the relative velocity of the sliding body and the bedrock is generated under the earthquake, which will cause a temporary reduction of the shear strength of the structural plane.It can be seen from Fig. 11 that the friction attenuation factor fluctuates in the range of 0-0.10, and its magnitude is closely related to the relative velocity.With the increase of the relative velocity, the relative friction attenuation factor increases, but its growth rate decreases.In addition, the frictional attenuation effect appears during the earthquake and disappears after the earthquake, and the friction attenuation factor after the earthquake is 0. This result is consistent with the results of Liu's research 5 on the relative velocity damage coefficient.
Figure 12 shows the time history of the equivalent friction angle during the earthquake.As can be seen from Fig. 12, the equivalent friction angle time history is mainly composed of two forms, one is a recoverable temporary reduction with seismic load fluctuations, and the other is a non-recoverable permanent reduction with slip displacement.The former is caused by the frictional attenuation effect associated with the relative velocity,  and the equivalent friction angle reduces non-permanently and recovers after the relative velocity disappears.The latter is caused by the slip deterioration effect associated with the slip displacement, the slip of the sliding body causes wear on the structural plane, resulting in a reduction in the roughness of the structural plane.Since the abrasion of the structural plane is permanent, the reduction of the equivalent friction angle due to the slip deterioration effect is irrecoverable, and the moment of sudden decrease of the equivalent friction angle corresponds to the moment of the slip of the sliding body.
Figure 13 shows the dynamic axial force and stiffness time history of the anchor cable during the earthquake.From Fig. 13, it can be seen that the dynamic axial force of the anchor cable fluctuates within a certain range in the early stage of the earthquake.When the dynamic axial force exceeds 136.5 kN, the brittle failure of the anchor cable occurs, and the dynamic axial force decreases abruptly to 0 kN.In this process, the anchor cable does not show a significant yield stage, indicating that the failure of the anchor cable is caused by the cracking of the anchorage device or the collapse of the anchor pier rather than the fracture of the anchor cable.Since the anchor cable is always in an elastic state before failure, the stiffness remains at the initial stiffness during this process and drops to 0 after the failure of the anchor cable.

The results obtained by the two calculation conditions
In order to analyze the seismic deterioration effect on the stability of the slope, two calculation conditions are used in this paper to calculate the anti-sliding force, sliding force and safety factor of the slope, respectively.Calculation condition A: The seismic deterioration effect is considered; Calculation condition B: The seismic deterioration effect is not considered.The time history of the anti-sliding force, the sliding force and the safety factor obtained under the two calculation conditions are shown in Fig. 14a-c, respectively.
As can be seen from Fig. 14a, the trend line of the anti-sliding force time history under calculation condition A shows both an asymptotical decrease and increase, as well as precipitous drops.In the early stage of the earthquake, compared with calculation condition B, the trend line of the anti-sliding force time history under calculation condition A gradually drifts downward, while the trend line of the anti-sliding force time history  www.nature.com/scientificreports/slowly rises with the decrease of the earthquake intensity in the later stage of the earthquake.This is because the friction attenuation effect is related to the relative velocity, the relative velocity of the sliding body that increases first and then decreases leads to the change of the trend line in the form of an asymptotically lower and then asymptotically higher under calculation condition A. The trend line of the anti-sliding force time history under calculation condition A first decreases abruptly at 26.29 s, and then there are several sudden drops of varying degrees with the progress of the earthquake.The reason for the sudden drops in the trend line is whether the anchor cable is broken or the structural plane deteriorates, which needs to be comprehensively analyzed in combination with the sliding force time history.From Fig. 14b, it can be found that compared with calculation condition B, the trend line of the time history of the sliding force under calculation condition A rises abruptly at 26.29 s, which is caused by the failure of anchor cables.There are two main reasons for the sudden increase in the trend line.On the one hand, the sliding body loses the normal force of anchor cables, which causes the friction force upward along the structural plane on the sliding body to decrease abruptly.On the other hand, the failure of anchor cables causes the tension upward along the structural plane on the sliding body to disappear suddenly.In addition, considering the variation trend of the time history of the anti-sliding force and the sliding force under calculation condition A in Fig. 14a,b, it can be seen that the first sudden drop of the trend line of the anti-sliding force time history is caused by the failure of anchor cables, and then the subsequent sudden drops of different degrees is caused by the slip deterioration effect of the structural plane.
As can be seen from Fig. 14c, the safety factor of the slope fluctuates to varying degrees during the earthquake.The minimum safety factor in the time history of the safety factor under calculation condition B is 0.7397, and the trend line of the time history of the safety factor remains horizontal and stabilizes at 1.3449.The minimum safety factor in the time history of the safety factor under calculation condition A is 0.6417, and the trend line of the time history of the safety factor decreases significantly and rises insignificantly, and finally stabilizes at 1.0756.Compared with calculation condition B, the minimum safety factor during the earthquake and the stable safety factor of the slope after the earthquake are reduced under calculation condition A. It can be seen that the seismic deterioration effect is considered in the calculation of the safety factor of the slope, which plays an important role in both the seismic design of the slope and the post-earthquake safety assessment of the slope.

The results obtained by the two calculation conditions
According to the time history of the safety factor under the two calculation conditions in Fig. 14c, the maximum and minimum safety factors of each curve are obtained, and the intervals of the maximum and minimum safety factors are discretized into multiple intervals with equal intervals.The safety factors of each interval are counted to obtain the number of safety factors in each discrete interval and the probability density histogram is plotted as shown in Fig. 15.As can be seen from Fig. 15, the probability density of the safety factor under calculation condition B is characterized by a unimodal distribution and obeys the normal distribution, which is consistent with Liu's research results 45 .However, the probability density of the safety factor under calculation condition A is characterized by a bimodal or multimodal and no longer obeys the normal distribution.The Gaussian mixture model can accurately fit the probability density of multimodal distribution, so the Gaussian mixture model is used for probability density estimation of the safety factor under calculation condition A.
The ChiChi waves with PGAs of 0.7-0.9g are used as the input seismic waves, and the probability distribution of the safety factor of the slope under the calculation condition A is fitted by using the normal distribution model and the Gaussian mixture model, respectively.The fitting effect is shown in Fig. 16.As can be seen from Fig. 16, the normal distribution model cannot accurately fit the probability distribution of the safety factor considering the seismic deterioration effect, while the Gaussian mixed model can accurately fit the probability distribution of the safety factor under this condition.Furthermore, the feasibility of using the number of seismic deteriorations of the slope as the number of components of the Gaussian mixture model is confirmed by accurately fitting the probability distribution of the safety factor of the slope under different seismic deterioration numbers.
In order to clarify the influence of seismic intensity on the failure probability of the slope, ChiChi waves with PGAs of 0.5-1.0g are used as input seismic waves.The cumulative distribution function curves shown in Fig. 17 are plotted and the failure probability of the slope shown in Table 2 is obtained by using different safety grades as evaluation indexes.As shown in Fig. 17, the variation trend of the cumulative distribution function curve under different seismic intensities is similar, and the main difference is in the safety factor range of 0.8-1.4.In addition, with the increase of seismic intensity, the failure probability of the slope under the same failure probability  www.nature.com/scientificreports/evaluation index increases.It can be seen from Table 2 that the failure probability of the slope increases with the increase of the safety grade of the slope under the same seismic intensity.The number of slips increases first and then remains unchanged with the increase of seismic intensity.It should be noted that the input maximum PGA in this paper is 1.0 g, but if the seismic intensity increases further, the sliding body may slip again.
The fitting results of the failure probability of the slope under different seismic intensities are shown in Fig. 18.As can be seen from Fig. 18, the failure probability increases linearly with the increase of seismic intensity under the same safety grade.In this paper, only the failure probability of the slope with input PGAs of 0.5-1.0g is calculated.In terms of the linear growth law of failure probability with seismic intensity, the failure probability of slope will further increase with the further increase of seismic intensity.

Discussion
In the example of this paper, anchor cables undergo brittle failure before the slip of the sliding body, and the static axial force of anchor cables suddenly decreases to 0. Therefore, the damage effect of anchor cables reduces the yield acceleration of the sliding body.However, it is worth noting that the static axial force of anchor cables may increase when the sliding body slides and anchor cables do not fail.In this case, the damage effect of anchor cables will increase the yield acceleration of the sliding body thereby reducing the slip displacement of the sliding body.No matter which of the above situations occurs, the important influence of the damage effect of anchor cables on the stability of the slope cannot be ignored.
Since anchor cables in this paper have the same specifications, they failed at the same time during the earthquake.However, in the actual project, the specifications of anchor cables at different positions in the seismic design of the slope may be different, and anchor cables may be damaged one by one during the earthquake.When determining the number of components of the Gaussian mixture model for the stability evaluation, it should be noted that the number of failures of anchor cables may be many.

Conclusion
(1) The seismic deterioration effect of the structural plane can be divided into the slip deterioration effect and the friction attenuation effect.The slip deterioration effect is related to the slip displacement of the sliding body, and the roughness of the structural plane decreases exponentially with the increase of slip displace- www.nature.com/scientificreports/ment and does not recover after the earthquake.The frictional attenuation effect is related to the relative velocity of the sliding body, which will cause temporary reductions in the shear strength of the structural plane and recovery after the earthquake.(2) Due to the resonance effect of the slope, the amplitude of the seismic acceleration time history of the sliding body may be greater than the amplitude of the input seismic acceleration time history.In addition, with the enhancement of the seismic deterioration effect, the yield acceleration of the sliding body gradually decreases.Therefore, the slip displacement calculated by the proposed method is larger than that calculated by the traditional Newmark displacement method, which is closer to the actual situation.(3) The minimum safety factor during the earthquake and the stable safety factor after the earthquake of the slope obtained by the calculation method in this paper are smaller than those of the calculation method without considering the seismic deterioration effect.It can be seen that the influence of seismic deterioration is considered in the calculation of the safety factor, which plays an important role in both the seismic design of the slope and the post-earthquake safety assessment of the slope.(4) For the seismic degradation effect of slopes, a slope stability evaluation method based on the Gaussian mixture model is proposed.The accuracy of the stability evaluation method and the feasibility of the number of seismic deteriorations as the number of components of the Gaussian mixture model are verified by an engineering example.The study on the slope failure probability under different seismic intensities shows that the slope failure probability is linearly correlated with the seismic intensity.

Figure 1 .
Figure 1.Tangential stress path diagram of structural plane.

Figure 3 .
Figure 3. Dynamic calculation model of anchored rocky slope.

Figure 4 .
Figure 4. Forces acting on the sliding body of the anchored rocky slope.

Figure 5 .
Figure 5. Dynamic calculation model of anchored rocky slope in the horizontal direction.

Figure 7 .
Figure 7. Schematic diagram of the anchored rocky slope.

Figure 8 .
Figure 8. ChiChi wave and the corresponding spectrum (a) Seismic acceleration time history (b) Spectrum.

Figure 9 .
Figure 9. Acceleration time history of the sliding body and input acceleration time history.

Figure 10 .
Figure10.Displacement calculated using the method in this paper.

Figure 13 .
Figure 13.Dynamic axial force and stiffness time history of the anchor cable.

Figure 14 .
Figure 14.Time history of anti-sliding force, sliding force and safety factor of the slope under different calculation conditions.(a) Anti-sliding force time history (b) Sliding force time history (c) Safety factor time history.

Figure 15 .Figure 16 .
Figure 15.Histogram of the safety factor of the slope under different calculation conditions.

Figure 17 .
Figure 17.Cumulative distribution function of the safety factor under different seismic intensities.
c is the damping of the structural plane, and k is the equivalent stiffness of the structural plane.The seismic acceleration of the sliding body at any time t can be obtained by iterative calculation of Eq. (36).

Table 1 .
The pending coefficient in the roughness and tangential stiffness degradation formulas.

Table 2 .
Failure probability of the slope under different seismic intensities.