Research on Slope Stability Based on Bayesian Gaussian Mixture Model and Random Reduction Method

Abstract

Slope stability analysis is conventionally performed using the strength reduction method with the proportional reduction in shear strength parameters. However, during actual slope failure processes, the attenuation characteristics of rock mass cohesion ($c$) and internal friction angle ($\phi$) are often inconsistent, and their reduction paths exhibit clear nonlinearity. Relying solely on proportional reduction paths to calculate safety factors may therefore lack scientific rigor and fail to reflect true slope behavior. To address this limitation, this study proposes a novel approach that considers the non-proportional reduction of $c$ and $\phi$, without dependence on predefined reduction paths. The method begins with an analysis of slope stability states based on energy dissipation theory. A Bayesian Gaussian Mixture Model (BGMM) is employed for intelligent interpretation of the dissipated energy data, and, combined with energy mutation theory, is used to identify instability states under various reduction parameter combinations. To compute the safety factor, the concept of a “reference slope” is introduced. This reference slope represents the state at which the slope reaches limit equilibrium under strength reduction. The safety factor is then defined as the ratio of the shear strength of the target analyzed slope to that of the reference slope, providing a physically meaningful and interpretable safety index. Compared with traditional proportional reduction methods, the proposed approach offers more accurate estimation of safety factors, demonstrates superior sensitivity in identifying critical slopes, and significantly improves the reliability and precision of slope stability assessments. These advantages contribute to enhanced safety management and risk control in slope engineering practice.

1. Introduction

Slope stability has long been a critical concern in geotechnical engineering [1]. Unstable slopes may trigger landslides, collapses, and other geological hazards, posing severe threats to mining safety and economics. Analyzing slope stability is therefore of significant importance [2]. Currently, the primary metric for evaluating slope stability is the safety factor [3], which represents the safety margin and is typically derived through the strength reduction method combined with numerical simulation.

Traditional strength reduction methods [4] proportionally reduce cohesion ($c$) and internal friction angle ($\phi$) until the slope reaches its limit state, using the reduction coefficient as the safety factor. As early as 2002, Zhao Shangyi et al. [5] combined finite element methods with strength reduction to calculate slope safety factors. Subsequently, Peng Xiaogang et al. [6] compared three numerical criteria for slope instability—numerical non-convergence, sudden displacement changes at characteristic points, and generalized plastic strain penetration—and proposed using computational convergence as the primary criterion. Sun Weijian et al. [7] introduced residual displacement increments based on characteristic point displacement mutation criteria to eliminate subjective errors in manual judgment of displacement mutations.

Beyond numerical methods, machine learning has been applied to predict slope stability. Singh, S.K. et al. [8] utilized multiple classification algorithms to predict slope stability under different joint types (parallel deterministic, cross-jointed, Baecher, Veneziano, and Voronoi). Malekian, M. et al. [9] employed the ARIMA model to forecast slope failure timing. Furthermore, considering the inherent uncertainties and variability in geotechnical parameters, the traditional deterministic safety factor framework exhibits certain limitations. Consequently, probabilistic methods have gained significant attention in slope stability assessment. For instance, Chakraborty Rubi et al. [10]. comprehensively reviewed currently available probabilistic approaches, providing detailed descriptions of various methodologies whilst elucidating the efficiency and drawbacks of each technique and their evolution in handling probabilistic slope analysis. Although these methods hold considerable potential in slope stability research, the field remains nascent and requires further refinement. Consequently, quantitatively determined deterministic safety factors typically remain the preferred metric for engineering practice.

As research on rock mass failure mechanisms deepens, studies reveal that the decay rates of $c$ and $\phi$ during progressive failure are inconsistent [11]. Traditional proportional reduction assumptions contradict this mechanism, prompting the development of dual-coefficient reduction methods. Yuan Wei et al. [12] established curves for critical slopes under varying geometries and rock densities, proposing a dual-parameter reduction method based on the shortest path to the critical curve. However, this method assumes specific reduction paths and lacks theoretical justification for its composite safety factor. Wu Shunchuan et al. [13] derived reduction ratios by assuming linear decay between peak and residual strength parameters, defining safety factors via reference slope probabilities. Although physically meaningful, this approach relies on linear decay assumptions. Lu Feng et al. [14] proposed a dual-coefficient method weighted by energy dissipation contributions during parameter reduction, yet it still presumes reduction paths. However, the approach of integrating energy evolution in this paper is noteworthy.

In summary, whilst current research addresses non-proportional reduction $c$ of and $\phi$, critical gaps remain as composite safety factors lack physical meaning and reduction paths are often arbitrarily assumed. This study addresses these issues by proposing a Bayesian Gaussian Mixture Model to classify slope stability under random reduction paths and define safety factors as shear strength ratios. Unlike prior machine learning approaches that discard slope model parameters, our method integrates numerical simulations of dissipated energy, ensuring rigorous parameter retention. By training unsupervised models on dissipated energy and reduction coefficients, the method accurately identifies critical instability thresholds, eliminating subjective thresholds in numerical simulations.

2. Construction of Slope Instability Criteria Using Bayesian Gaussian Mixture Model

2.1. Energy Dissipation Mutation Criterion

During the progressive failure of slopes, the energy field within the slope system and its surrounding environment remains in a state of dynamic equilibrium [15]. Studying slope failure from the perspective of energy evolution allows us to bypass complex stress–strain analyses and avoid tracking detailed damage states of geomaterials at each intermediate step. Instead, the focus shifts to the pathways of energy input, storage, transformation, and dissipation. In slope systems, part of the externally applied energy is stored as elastic strain energy within the material, whilst the remaining energy is gradually dissipated as damage-induced energy through deformation and damage processes [16,17]. According to the first law of thermodynamics, this process can be expressed as follows:

(1)$U=U^d+U^e$

Under static equilibrium, the energy input is primarily driven by the self-weight of the slope. Thus, Equation (1) can be rewritten as follows:

(2)$U^g=U^d+U^e$

where $U$ represents the total energy input; $U^d$ denotes the dissipated energy; $U^e$ is the elastic strain energy; and $U^g$ stands for the gravitational potential energy.

Both the storage of elastic strain energy and the dissipation of energy through deformation in a slope have inherent capacity limits. When the total energy input from external sources exceeds the slope’s energy-bearing threshold, the surplus energy that cannot be stored or dissipated rapidly transforms into kinetic energy. At this point, the slope reaches a state of limit equilibrium, triggering sliding and leading to failure. By simplifying slope stability analysis to the judgment of energy evolution patterns, the critical conditions for slope instability can be identified more efficiently and accurately. The energy evolution process of the slope can then be expressed by the following equation:

(3)$U^g=U^d+U^e+U^k$

where $U^k$ is the kinetic energy.

Taking any volumetric element within the slope system as the research object, the gravitational potential energy $U^g$, dissipated energy $U^d$ elastic strain energy, $U^e$ and kinetic energy $U^k$ can be calculated using the following formulas:

(4)$U^g={\displaystyle {\int }_V\rho ghdV}$

(5)$U^d={\displaystyle {\int }_v\sigma }:d{\epsilon }^p$

(6)$$\begin{gathered} U^e={\displaystyle {\int }_V\frac{1}{2}}\left({\sigma }_x{\epsilon }_x+{\sigma }_y{\epsilon }_y+{\sigma }_z{\epsilon }_z\right)dV= \\ {\displaystyle {\int }_V\frac{1}{2E}}\left[{\sigma }_x^2+{\sigma }_y^2+{\sigma }_z^2-2u\left({\sigma }_x{\sigma }_y+{\sigma }_y{\sigma }_z+{\sigma }_z{\sigma }_x\right)\right]dV \end{gathered}$$

(7)$U^k={\displaystyle {\int }_V\frac{1}{2}}\rho v^{2}dV$

where $\rho$ represents the density; $g$ is the gravitational acceleration; $h$ is the height of the element’s centroid; $V$ is the volume of the element; $\sigma$ is the stress tensor; “:” is the inner product operation; ${\epsilon }^p$ is the plastic strain tensor; ${\sigma }_x$, ${\sigma }_y$, and ${\sigma }_z$ are the principal stresses in the x, y, and z directions, respectively; ${\epsilon }_x$, ${\epsilon }_y$, and ${\epsilon }_z$ are the principal strains in the x, y, and z directions, respectively; $E$ is the elastic modulus of the material; $\upsilon$ is the Poisson’s ratio of the material; and $v$ is the velocity at the element’s centroid.

The energy dissipation process of slopes exhibits significant differences between unstable and stable states. During the stable phase, the dissipated energy typically follows a relatively smooth trend, with fluctuations within a certain range when influenced by external factors. However, as the slope approaches the critical state of instability, the internal structure of the rock and soil undergoes drastic changes, such as rapid crack propagation and imbalances in particle interactions, causing the energy dissipation to deviate from its original smooth pattern and exhibit abrupt changes [18]. Although the theoretical understanding of the distinct energy dissipation behaviors in different slope states is clear, accurately identifying the mutation point is crucial for the precision of subsequent safety factor calculations. Previous studies have attempted to locate the mutation point mathematically by constructing regression equations between reduction coefficients and dissipated energy [19]. However, these methods struggle to achieve precise identification in complex scenarios involving dual reduction coefficients due to the inherent limitations of the equations and the complexity of multi-variable interactions.

Unsupervised clustering algorithms [20], as important techniques in the field of data mining, can achieve automatic classification based on similarity metrics among data points under the condition of unlabeled datasets. This approach does not require pre-setting target variables and mines the inherent data structure with high adaptability. Common unsupervised clustering algorithms have distinct characteristics. The K-means algorithm realizes data partitioning by iteratively updating cluster centers, featuring high computational efficiency and rapid convergence, which is suitable for datasets with uniform sample distribution and approximately spherical cluster shapes; hierarchical clustering algorithms intuitively present the hierarchical relationships of data through a dendrogram structure, allowing flexible adjustment of the number of clusters, especially suitable for exploratory data analysis and scenarios requiring the revelation of clustering hierarchical relationships; the Bayesian Gaussian Mixture Model [21] (BGMM), based on probability theory, completes clustering by estimating the probability that data belongs to each Gaussian distribution. It can not only output the probability values of samples belonging to different categories but also effectively addresses data uncertainty and fuzziness, which has been widely applied in fields such as image recognition and speech processing.

In slope stability analysis, compared with the method of constructing regression equations between reduction coefficients and dissipated energy to find mathematical mutation points, unsupervised clustering methods in the field of machine learning show unique potential in identifying dissipated energy mutations under dual-parameter reduction. The stability (or instability) of slopes is essentially a probabilistic problem—with minor changes in reduction coefficients, their mechanical responses (such as dissipated energy) exhibit continuous gradual changes. Near the critical region, dissipated energy data points do not abruptly jump from “stable” to “unstable” but exist in a fuzzy region with overlapping states, reflecting the uncertainty of state evolution in slope systems. The core advantage of the Bayesian Gaussian Mixture Model lies in modeling data as a weighted mixture of several probability distributions, which can accurately describe the overlapping data distributions generated by different potential states. In this study, “stable” and “unstable” are regarded as two potential states. By assigning each data point the probability (membership degree) of belonging to the “stable class” and “unstable class”, BGMM perfectly matches the physical essence of state uncertainty in the critical region. Therefore, this study uses this model to perform clustering analysis on slope dissipated energy data, achieving a scientific classification of slope states into two categories: “stable” and “unstable”. To objectively demonstrate the applicability of BGMM, Table 1 systematically compares the core characteristics of each algorithm.

2.2. Gaussian Mixture Model (GMM)

The Gaussian Mixture Model (GMM) is a probabilistic model and an unsupervised clustering machine learning algorithm [22]. It assumes that all data points $X=\left\{x_1,x_2,\dots ,x_n\right\}$ are independently sampled from the same probability distribution. Its probability density function can be expressed as a linear combination of multiple Gaussian functions, which are as follows:

(8)$p(x)={\displaystyle \sum _{k=1}^K{\pi }_k}\mathcal{N}\left(x;{\mu }_k,{\displaystyle {\sum }_k }\right)$

(9)$s.t.{\displaystyle \sum _{k=1}^K{\pi }_k}=1,0\le {\pi }_k\le 1,\forall k=1,\dots ,K$

where $K$ represents the number of Gaussian distributions, i.e., the number of mixture components; ${\pi }_k$ is the prior probability, i.e., the mixing coefficient; $\mathcal{N}\left(x;{\mu }_k,{\displaystyle {\sum }_k}\right)$ is the probability density function of the k-th Gaussian distribution, where ${\mu }_k$ is the mean vector of the k-th Gaussian distribution; and ${\displaystyle {\sum }_k }$ is the covariance matrix of the k-th Gaussian distribution.

The specific form of the probability density function $\mathcal{N}\left(x;{\mu }_k,{\displaystyle {\sum }_k}\right)$ for the Gaussian distribution is as follows:

(10)$\mathcal{N}\left(x;{\mu }_k,{\displaystyle {\sum }_k }\right)={\displaystyle \frac{1}{{(2\pi )}^{\frac{d}{2}}{\left|{\displaystyle {\sum }_k }\right|}^{\frac{1}{2}}}}{e}^{-{\displaystyle \frac{1}{2}}{\left(x-{\mu }_k\right)}^T{\displaystyle \sum _k^{-1}\left(x-{\mu }_k\right)}}$

where $d$ is the dimensionality of the data point x; $\left|{\displaystyle {\sum }_k}\right|$ is the determinant of the covariance matrix ${\displaystyle {\sum }_k}$; and ${\displaystyle \sum _k^{-1} }$ is the inverse matrix of the covariance matrix ${\displaystyle {\sum }_k}$.

The core idea of the Gaussian Mixture Model lies in the assumption that the data is generated from a mixture of multiple Gaussian distributions, with each Gaussian distribution corresponding to a cluster. The model fits the probability distribution of the data by estimating the mean, covariance matrix, and mixing coefficient for each cluster. Its operational mechanism involves first identifying the Gaussian function center corresponding to each sample point, then calculating the posterior probability of the sample point belonging to each cluster, and finally assigning the sample point to the appropriate cluster based on these probabilities [23]. In practical implementation, within the Gaussian Mixture Model the mean, covariance matrix, and mixing coefficient for each Gaussian distribution can be estimated using either the maximum likelihood estimation method or the Expectation Maximization (EM) algorithm. Once these parameters have been estimated, the sample points can be assigned to the cluster corresponding to the Gaussian function with the highest probability density by comparing the probability density of each sample point across the Gaussian distributions.

The likelihood estimation method is a commonly used statistical inference approach. Its fundamental idea is that, given a probability model and a set of observed data, the goal of likelihood estimation is to find the model parameter values that are most likely to have generated the observed data. The likelihood function is the core tool of likelihood estimation, representing the probability of observing the data given the parameter values. It is typically denoted as $L\left(\theta ;x\right)$, where $\theta =\left\{{\theta }_1,{\theta }_2,\dots ,{\theta }_n\right\}$ are the unknown parameters of the model, and the observed data is $X=\left\{x_1,x_2,\dots ,x_n\right\}$, with the probability density function $f\left(x_i;\theta \right)$. The expression for the likelihood function is as follows:

(11)$L(\theta ;x)=f(x_1;\theta )f(x_2;\theta )\dots f(x_n;\theta )={\displaystyle \prod _{i=1}^{n}f}(x_i;\theta )$

By replacing the probability density function $f\left(x_i;\theta \right)$ with the probability density function of the Gaussian Mixture Model, the likelihood function of the Gaussian Mixture Model can be derived as follows:

(12)$L\left(\theta ;x\right)={\displaystyle \prod _{i=1}^N{\displaystyle \sum _{k=1}^K{\pi }_k}}\mathcal{N}\left(x;{\mu }_k,{\displaystyle {\sum }_k }\right)$

2.3. Bayesian Estimation Algorithm

The likelihood equation of the Gaussian Mixture Model typically lacks an analytical solution. Combining the Bayesian algorithm for iterative computation is one method for solving the Gaussian Mixture Model. Bayesian estimation [24], as a statistical inference method based on Bayes’ theorem, is fundamentally different from the traditional frequentist approach. In the Bayesian estimation framework, the unknown parameters of the model are treated as random variables, unlike the traditional frequentist approach, which considers them as fixed but unknown constants. The key point lies in accurately calculating the posterior distribution of the parameters based on prior knowledge and observed data. Here, the prior distribution carries the accumulated knowledge and assumptions about the parameters before observing the data. Subsequently, using Bayes’ formula the prior distribution and the likelihood function of the observed data are integrated to derive the posterior distribution. Finally, parameter estimation and inference are conducted based on the obtained posterior distribution. The specific steps for estimating the parameters of the Gaussian Mixture Model using the Bayesian algorithm are as follows.

In the Gaussian Mixture Model, the parameters to be estimated are the mixing coefficient ${\pi }_k$, the mean vector ${\mu }_k$, and the covariance matrix ${\displaystyle {\sum }_k }$. The probability density functions of their prior distributions can be expressed as follows:

(13)$p\left({\pi }_k|\alpha \right)={\displaystyle \frac{\mathrm{\Gamma }({\displaystyle \sum _{k=1}^K{\alpha }_k})}{{\displaystyle \prod _{k=1}^K\mathrm{\Gamma }}({\alpha }_k)}}{\displaystyle \prod _{k=1}^K{\pi }_k^{{\alpha }_k-1}}$

(14)$p({\mu }_k)={\displaystyle \frac{1}{{\left(2\pi \right)}^{\frac{d}{2}}{\left|V_0\right|}^{\frac{1}{2}}}}{\mathrm{e}}^{-{\displaystyle \frac{1}{2}}{({\mu }_k-m_0)}^T{V}_0^{-1}({\mu }_k-m_0)}$

(15)$p({\mathrm{\Sigma }}_k)={\displaystyle \frac{|S_0{|}^{\frac{{\nu }_0}{2}}}{2^{\frac{{\nu }_{0}d}{2}}{\mathrm{\Gamma }}_d(\frac{{\nu }_0}{2})}}|{\mathrm{\Sigma }}_k{|}^{-{\displaystyle \frac{{\nu }_0+d+1}{2}}}\exp \left(-{\displaystyle \frac{1}{2}}tr(S_0{\mathrm{\Sigma }}_k^{-1})\right)$

where the mixing coefficient follows a Dirichlet distribution, ${\pi }_k~Dir(\alpha )$, where $\alpha =({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_k)$ is the parameter vector of the Dirichlet distribution; the mean vector follows a normal distribution, ${\mu }_k~\mathcal{N}(m_0,V_0)$, where $m_0$ is the prior mean and $V_0$ is the prior covariance; the covariance matrix follows an inverse Wishart distribution, ${\displaystyle {\sum }_k }~IW(v_0,S_0)$, where $v_0$ and $S_0$ are the parameters of the inverse Wishart distribution; and $d$ is the data dimensionality.

Assuming a set of observed vectors as $X={\left[x_1,x_2,\dots ,x_n\right]}^T$, and letting the unknown parameter be $\theta$, the posterior probability of $\theta$ can be expressed as follows:

(16)$p(\theta |X)={\displaystyle \frac{p(X|\theta )p(\theta )}{p(X)}}$

where $p(\theta |X)$ is the posterior probability of $\theta$; $p(X|\theta )$ is the probability density function of $X$; $p(\theta )$ is the prior estimate of $\theta$; and $p(X)$ is the observed probability of $X$.

Under the given observed event $\theta$, $p(X)$ is a measurable constant, and $p(X|\theta )$ can be considered the likelihood function of event $\theta$. Thus, the posterior probability $p(\theta |X)$ is proportional to $p(X|\theta )p(\theta )$, which is shown as follows:

(17)$p(\theta |X)\propto p(X|\theta )p(\theta )$

From Equation (17), the posterior distribution of the mixing coefficient ${\pi }_k$ is ${\pi }_k~Dir({\alpha }^{*})$, where $\alpha *=({\alpha }_1+n_1,{\alpha }_2+n_2,\dots ,{\alpha }_k+n_k)$, and $n_k$ is the number of samples belonging to the k-th Gaussian component. The posterior distribution of the mean vector ${\mu }_k$ is ${\mu }_k~\mathcal{N}(m^{*},V^{*})$ where $V^{*}={(nk{V}_0^{-1}+V_0^{-1})}^{-1}$ and $m^{*}={\displaystyle \frac{V^{*}(n_k{\overline{x}}_k{V}_0^{-1}{m}_0)}{n_k{V}_0^{-1}+V^{*}}}$. The posterior distribution of the covariance matrix ${\displaystyle {\sum }_k }$ is ${\displaystyle {\sum }_k }~IW(v^{*},S^{*})$, where ${\nu }^{*}={\nu }_0+n_k$ and $S^{*}=S_0+S_k$.

By calculating the expected value of the posterior distribution, the parameter estimates can be obtained. The specific principle of the Bayesian estimation algorithm for computing the parameter values of the Gaussian Mixture Model is illustrated in Figure 1:

3. Comprehensive Safety Factor

3.1. Definition of Comprehensive Safety Factor

The traditional strength reduction method defines the safety factor by continuously reducing the strength parameters of the slope until it reaches the critical state, using the reduction coefficient at that point as the safety factor. In the dual-parameter reduction method, since there are two reduction coefficients for cohesion and internal friction angle defining the comprehensive safety factor is the core issue of the dual-parameter reduction method [25]. If the comprehensive safety factor is simply defined as the average of the reduction coefficients for cohesion and internal friction angle from a mathematical perspective, it will lack physical meaning and will be difficult to gain recognition for guiding real-world engineering problems. In traditional safety factor definitions, the ratio of the slope’s anti-sliding force to the sliding force is also used, which is shown as follows:

(18)$\left\{\begin{array}{l}{F}_s={\displaystyle \frac{f_{Anti-sliding}}{f_{sliding}}}\\ {\displaystyle \frac{{f_{Anti-sliding}}^{\prime }}{{f_{sliding}}^{\prime }}}=1\end{array}\right.$

where $F_s$ is the slope safety factor; $f_{Anti-sliding}$ and $f_{sliding}$ are the anti-sliding force and sliding force of the initial slope, respectively; and ${f_{Anti-sliding}}^{\prime }$ and ${f_{sliding}}^{\prime }$ are the anti-sliding force and sliding force of the critical slope, respectively, at which point the slope safety factor is 1.

Since the reduction in cohesion and internal friction angle does not change the sliding force, i.e., $f_{sliding}={f_{sliding}}^{\prime }$, it follows that:

(19)$F_s={\displaystyle \frac{f_{Anti-sliding}}{{f_{Anti-sliding}}^{\prime }}}={\displaystyle \frac{cA+{\sigma }_{\mathrm{n}}A\tan \phi }{cA/k_c+{\sigma }_{\mathrm{n}}A\tan \phi /k_{\phi }}}$

where $c$ and $\phi$ are the initial cohesion and internal friction angle, respectively; ${\sigma }_{\mathrm{n}}$ is the normal stress; $A$ is the sliding surface area; and $k_c$ and $k_{\phi }$ are the reduction coefficients for cohesion and internal friction angle, respectively.

In the traditional strength reduction method, $k_c=k_{\phi }$. If it is assumed that the normal stress on the slope element remains approximately equal after reducing only cohesion and the internal friction angle, the definition of the safety factor in the traditional strength reduction method can be obtained as follows:

(20)$F_s=k_c=k_{\phi }$

Thus, it can be considered that the comprehensive safety factor defined by Equation (19) holds equivalent reference significance to the safety factor defined by the traditional strength reduction method. Simplifying this equation yields the following:

(21)$F_s={\displaystyle \frac{cA+{\sigma }_{\mathrm{n}}A\tan \phi }{cA/k_c+{\sigma }_{\mathrm{n}}A\tan \phi /k_{\phi }}}={\displaystyle \frac{c+{\sigma }_{\mathrm{n}}\tan \phi }{c/k_c+{\sigma }_{\mathrm{n}}\tan \phi /k_{\phi }}}$

Introducing the following shear strength formula:

(22)$\tau =c+{\sigma }_{\mathrm{n}}\tan \phi$

The expression for the safety factor can be defined as follows:

(23)$F_s={\displaystyle \frac{\tau }{{\tau }^{\prime }}}$

where $\tau$ and ${\tau }^{\prime }$ are the shear strengths of the rock mass in the initial state and critical state, respectively.

3.2. Stress Analysis of the Element

Figure 2 illustrates the stress state of an arbitrary infinitesimal element within the slope. On this element, a pair of normal stresses ${\sigma }_x$ and ${\sigma }_y$, as well as shear stresses ${\tau }_x$ and ${\tau }_y$, are applied. $\beta$ is the angle formed between the principal stress direction and the positive direction of the x-axis.

For the signs of stresses $\sigma$, $\tau$, and angle $\beta$, the following rules apply:

(1). Normal stresses and principal stresses: When they exhibit a compressive state, their values are positive. If they are in a tensile state, their values are negative.

(2). Shear stresses: If a shear stress causes the element to rotate counterclockwise, its value is positive. If it causes clockwise rotation, its value is negative.

(3). Angle $\beta$: Its value ranges between 0° and 90°. When rotating from the positive direction of the x-axis toward the outward normal direction of the acting surface, if the rotation is counterclockwise, $\beta$ is positive. If it is clockwise, $\beta$ is negative.

Based on the above rules, analyzing Figure 3 reveals that the normal stress in the figure is positive, the angle is also positive, and the shear stress is negative.

According to the definition of principal stress, the calculation formulas for the maximum principal stress ${\sigma }_1$, the minimum principal stress ${\sigma }_3$, and the included angle $\beta$ on the principal plane are as follows:

(24)$\left.\begin{array}{l}{\sigma }_1\\ {\sigma }_3\end{array}\right\}={\displaystyle \frac{{\sigma }_x+{\sigma }_y}{2}}\pm \sqrt{{\left({\displaystyle \frac{{\sigma }_x-{\sigma }_y}{2}}\right)}^2+{\tau }_x^2}$

(25)$\beta =-{\displaystyle \frac{1}{2}}\arctan {\displaystyle \frac{2{\tau }_x}{{\sigma }_x-{\sigma }_y}}$

Figure 3 is a schematic diagram for calculating the stress on any inclined plane within the infinitesimal element.

Within the infinitesimal element, an arbitrary inclined plane is selected. Based on the static equilibrium conditions, it can be concluded that if the values of ${\sigma }_1$ and ${\sigma }_3$ are known then on an inclined plane at an angle $\alpha$ (where the angle $\alpha \in \left[0°,90°\right]$) relative to the acting surface of ${\sigma }_1$, the normal force ${\sigma }_n$ and shear stress $\tau$ will exhibit specific numerical relationships, as described below:

(26)${\sigma }_{\mathrm{n}}={\displaystyle \frac{{\sigma }_1+{\sigma }_3}{2}}+{\displaystyle \frac{{\sigma }_1-{\sigma }_3}{2}}\cos 2\alpha$

(27)$\tau ={\displaystyle \frac{{\sigma }_1-{\sigma }_3}{2}}\sin 2\alpha$

(28)$\alpha ={45}^{°}+\phi /2$

3.3. Calculation Process of Comprehensive Safety Factor

Using FLAC3D version 6.00.69 software (a three-dimensional finite difference program), the parameters ${\sigma }_1$ and ${\sigma }_3$ can be directly obtained. For calculating the shear strength on the potential sliding surface of the initial slope, the angle ${\alpha }_0$ between the failure plane ${\sigma }_1$ and the acting surface is determined based on Equation (28). Substituting ${\alpha }_0$ into Equation (26) and combining it with Equation (22), the shear strength ${\tau }_0$ of each element on the potential sliding surface in the initial state can be calculated. Similarly, for the shear strength on the potential sliding surface of the critical slope, following the same steps, the shear strength ${\tau }_{critical}$ of each element on the potential sliding surface in the critical state can be calculated. Subsequently, using a random sampling algorithm, elements are selected from the potential sliding surface, and the mean shear strengths of the slope elements in the critical and initial states are calculated, denoted as ${\tau }^{\prime }$ and $\tau$, respectively. Finally, substituting ${\tau }^{\prime }$ and $\tau$ into Equation (23) yields the final safety factor.

4. Case Verification

4.1. Numerical Model

A three-dimensional model, as shown in Figure 4, is established. The height of the slope bench is 15 m, the width of the bench is 5 m, and the overall dimensions of the model are 55 m in height, 60 m in length, and 30 m in width. The material parameters are listed in Table 2. A grid size of 1 m-by-1 m was adopted (validated by preliminary tests conducted prior to the numerical simulation, confirming that this size meets the requirements and ensures the accuracy and reliability of the simulation results).

4.2. Establishment of Training Sets for Reduction Coefficients and Dissipated Energy

To build a training set, FLAC3D was utilized for conducting numerical simulations of a slope model employing the Mohr–Coulomb constitutive model. A computational program was developed based on Equation (5) to calculate the dissipated energy under independent reductions in cohesion ($c$) and internal friction angle ($\phi$). The computational program utilizes FLAC3D’s built-in FISH programming environment. The range of reduction coefficients for $c$ and $\phi$ was set to [1,1.5] with an initial step size of 0.1. To optimize computational efficiency, data requiring over 100,000 convergence steps in FLAC3D were eliminated. The local refinement of step sizes (0.01 and 0.001) was applied to specific intervals, resulting in a final dataset of 609 entries. Partial results are summarized in Table 3, while the complete dataset is illustrated in Figure 5.

In geotechnical engineering, to solve the system of differential equations governing the medium’s response (equilibrium, geometric compatibility, and constitutive equations) along with boundary conditions, the Finite Difference Method first discretizes the problem domain and establishes a computational grid. Subsequently, at each grid node, it utilizes Taylor series expansions to approximate the spatial derivatives in the governing equations as finite differences using function values at adjacent nodes. These finite difference approximations are then substituted back into the original differential equations, thereby transforming them over the entire solution domain into a system of linear or nonlinear algebraic equations with respect to nodal unknowns (such as displacements or stresses). Finally, the approximate solution is obtained by numerically solving this algebraic equation system.

A visualization of dissipated energy contours for selected slope cases is shown in Figure 6.

When cohesion ($c$) and internal friction angle ($\phi$) are reduced proportionally, this scenario can be considered a special case within random reduction paths. The corresponding variations in dissipated energy are illustrated in Figure 7.

4.3. Construction of Bayesian Gaussian Mixture Model

For the construction of the Bayesian Gaussian Mixture Model (BGMM), Python code was developed for the training phase. The pre-obtained training dataset was saved in an Excel file. Subsequently, the implemented program read this file to train the BGMM and produced the requisite charts and evaluation metrics. Figure 8 illustrates the Gaussian component decomposition diagram of the trained model. The data histogram shown in the “Data Histogram” figure displays the distribution histogram of the original data. The data distribution exhibits significant bimodal characteristics, indicating the presence of two different scales of energy dissipation behavioral patterns in the system. Component 1 and Component 2 confirm that subsets of the original data can be described by Gaussian distributions.

The model training outcomes are summarized in Table 4, which presents the classification results of the raw data training (partial data). For a comprehensive visualization of all data points, the complete results are illustrated in Figure 9,

4.4. Performance Evaluation and Metric Analysis of the Bayesian Gaussian Mixture Model

To comprehensively evaluate the classification performance of the Bayesian Gaussian Mixture Model (BGMM) for slope energy dissipation patterns, this study uses the comparison of data distributions and four categories of complementary indicators to comprehensively verify the model’s effectiveness. By programming in Python (v.3.11.) to retrieve the model output results (see Figure 10 and Table 5), the evaluation covers the following: internal compactness indicators (silhouette score, quantifying intra-cluster aggregation, and inter-cluster separation); statistical fitting indicators (log-likelihood, testing the probability fitting goodness of the model to energy dissipation data distributions); clustering validity indices (Calinski–Harabasz Index and Davies–Bouldin Index, providing standardized assessments of clustering distinguishability and compactness); and stability indicators (Weight Variability and Mean Separation, verifying the model’s anti-interference ability and robustness in identifying slope critical states).

These four categories of indicators complement each other, constructing a cross-validation system from the following four dimensions: clustering structure rationality, probability distribution matching degree, statistical validity, and engineering application stability. This ensures that the model maintains both theoretical rigor and engineering practicality in capturing energy mutation characteristics.

This figure illustrates the comparison between synthetic data (generated by sampling from the trained BGMM) and original data. The peaks of both the original and synthetic data align in the same position, indicating that the model successfully captured the primary modal characteristics of slope dissipated energy. The left-skewed extensions of both datasets exhibit a steep decline, reflecting the rapid convergence of low energy dissipation. The right tails of the original and synthetic data show an extremely high degree of overlap, demonstrating that the synthetic data can accurately replicate the features of the original data. To quantitatively evaluate the similarity between the distributions of the original and synthetic data, the following three statistical parameters are presented in Figure 10: the Overlap Coefficient (0.8180), Bhattacharyya Coefficient (0.9386), and Jensen–Shannon Divergence (0.0490). Specifically, the Overlap Coefficient (0.8180) indicates a high spatial consistency in energy dissipation ranges. The Bhattacharyya Coefficient (0.9386), approaching 1, demonstrates an accurate replication of the distribution morphology of the original data. The Jensen–Shannon Divergence is below 0.05 (0.0490), implying that the information loss during distribution reconstruction is less than 5%. Collectively, these metrics verify that the BGMM accurately captures the statistical patterns of dissipated energy data. This enables the trained model to replicate the transition of slope energy dissipation from stable to unstable states, thereby achieving the reliable classification of slope models into two distinct categories: stable and unstable.

(1) Log-likelihood value: This value measures the model’s fit to the data, representing the logarithm of the probability of observing the data given the model parameters. A positive value indicates the model’s explanatory power, with higher values signifying better performance [26]. A negative value reflects model bias, with smaller absolute values being preferable. The log-likelihood value of the model in this study is −0.69, with an extremely small absolute value, indicating that the model effectively captures the distribution characteristics of the data and can accurately identify the critical state of slopes.

(2) Silhouette coefficient: This metric evaluates the clarity of sample assignments in clustering results, reflecting the balance between intra-cluster cohesion and inter-cluster separation. Values closer to 1 indicate higher clustering quality [27]. The silhouette coefficient of the model in this study is 0.88, approaching the theoretical optimum (1.0), demonstrating significant separability between the energy dissipation patterns of stable and unstable slopes and low error rates in critical state identification.

(3) Calinski–Harabasz index: This index assesses clustering effectiveness by the ratio of inter-cluster variance to intra-cluster variance, with higher values indicating greater inter-cluster differences [28]. The model in this study achieves an exceptionally high value of 916.52 (typical model range: 200–500), revealing substantial differences in energy distribution between stable and unstable states, which aligns well with the nonlinear instability mechanisms of slopes.

(4) Davies–Bouldin index: This metric quantifies inter-cluster similarity, with smaller values indicating lower risk of cluster overlap [29]. The index for the model in this study is 0.64 (significantly better than the threshold of 1.0), proving that the overlap region between stable and unstable modes accounts for less than 5%, with minimal misjudgment risk and robust performance even in noisy monitoring data.

(5) Weight volatility: This parameter characterizes the stability of the model’s component weights, with smaller values indicating stronger consistency in identifying primary and secondary states. The volatility of the model in this study is 0.36, demonstrating highly stable weight assignments for distinguishing between stable and unstable slope states, making it widely applicable in engineering contexts.

(6) Mean difference degree: This value quantifies the degree of difference between two probability distributions, model components, or system states. Its core objective is to transform abstract differences into computable numerical values to support decision-making. The mean difference degree in this study is 2.02 (in normalized space), reflecting the significance of changes in dissipated energy, which aligns with the abrupt energy dissipation observed during slope instability.

4.5. Calculation of Safety Factor

The constructed Bayesian Gaussian Mixture Model successfully classifies slopes under different reduction coefficients into stable and unstable categories. For data in the unstable state, if no other unstable state data exists where the reduction coefficients of cohesion and internal friction angle are both smaller than those of this dataset, the dataset is defined as critical slope data.

To address potential deviations caused by insufficient precision in the reduction step size, further processing is applied to the critical slope data: data points with dissipated energy exceeding 20% of the average value are excluded. The critical slope data refined according to this rule are detailed in Table 6.

Taking the slip surface of the critical slope in Group 1 as a case study, the initial slope is labeled as Group 0. Seven elements were randomly and uniformly sampled, as illustrated in Figure 11. The stress parameters (${\sigma }_1$, ${\sigma }_3$) of the selected elements under both initial and critical state (Group 1) were extracted using FLAC3D software. The average shear strength of the selected elements was calculated as follows:

Shear strength on the potential slip surface of the initial slope: The average angle (${\alpha }_0=60.02$) between the ${\sigma }_1$ action surface and the failure surface of the selected elements under the initial state was calculated using Equation (28). Substituting ${\alpha }_0$ into Equation (26) ${\alpha }_0$ and combining it with Equation (22), the average shear strength of the selected elements on the potential slip surface was determined as ${\tau }_0=193235.64Pa$.

Shear strength on the potential slip surface of critical slope 1: The average angle (${\alpha }_1=56.86°$) between the failure surface ${\sigma }_1$ and the action surface of the selected elements was calculated using Equation (28). Substituting ${\alpha }_1$ into Equation (26) and combining it with Equation (22), the average shear strength of the selected elements on the slip surface of the critical slope was determined as ${\tau }_1=141124.40Pa$.

Similarly, the shear strengths for all critical slopes were calculated and are listed in Table 7.

From this, it can be determined that $\tau =193.24kPa$ and ${\tau }^{\prime }=144.43kPa$, resulting in a final safety factor of 1.338.

4.6. Comparison of Results

In the traditional strength reduction method, the same reduction coefficient is applied synchronously, and the reduction coefficient of the critical slope is taken as the safety factor. The safety factor $F_1$ can be obtained iteratively through numerical simulation software, as shown in Figure 12.

When the dual-coefficient reduction method is applied, two safety factors are obtained for cohesion and internal friction angle, respectively. Some studies propose using the average of the two reduction coefficients as the comprehensive safety factor $F_2$, which is shown as follows:

(29)$F_2={\displaystyle \frac{k_c+k_{\phi }}{2}}$

Alternatively, the arithmetic square root of the product of the two reduction coefficients can be used as the safety factor $F_3$, which is shown as follows:

(30)$F_3=\sqrt{k_c·k_{\phi }}$

Other studies have proposed integrating statistical principles by analyzing large datasets of slopes under diverse conditions to derive the comprehensive safety factor $F_4$ through the following mathematical fitting:

(31)$F_4={\displaystyle \frac{\sqrt{2}{k}_c{k}_{\phi }}{\sqrt{k_c^2+k_{\phi }^2}}}$

Other studies define the comprehensive safety factor $F_5$ based on the shortest path of strength reduction, which is shown as follows:

(32)$\begin{array}{l}{F}_5={\displaystyle \frac{1}{1-R/\sqrt{2}}}\\ R=\sqrt{{\left(1-{\displaystyle \frac{1}{k_c}}\right)}^2+{\left(1-{\displaystyle \frac{1}{k_{\phi }}}\right)}^2}\end{array}$

The safety factors calculated using the methods described above are compiled in Table 8.

The results in Table 8 indicate that the safety factor derived from the proposed method is 2.5% lower than that obtained using traditional methods. This method achieves higher precision in classifying slopes into stable and unstable states. By employing a Bayesian Gaussian Mixture Model to process dissipated energy data it demonstrates significant advantages in complex geological conditions, particularly for slopes near critical states where it provides more accurate stability assessments.

This method pioneers a new paradigm in slope stability analysis. On the one hand, scholars have conducted extensive research on the reduction in cohesion and internal friction angle under various slope conditions, while the random reduction method proposed in this study can be practically applied by integrating specific reduction paths derived from these investigations. Furthermore, it establishes a novel direction for slope stability research. With continuous refinement and optimization of algorithms, this method will further enhance its ability to distinguish between stable and unstable slopes with greater precision.

5. Conclusions

Integrating energy dissipation theory and the Bayesian Gaussian Mixture Model (BGMM), this study innovatively proposes a slope stability analysis method that incorporates random reduction paths. By numerically simulating dissipated energy under varying reduction coefficients and using energy evolution characteristics as criteria, a comprehensive safety factor with clear physical significance is defined. The main conclusions are as follows:

(1) The limitations of conventional non-proportional reduction methods (e.g., fixed reduction paths and ambiguous physical meaning of safety factors) are addressed through the proposed random-path strength reduction method. This approach resolves the subjectivity in path selection and clarifies the physical interpretation of safety factors.

(2) A BGMM–dissipated energy joint criterion is established to achieve precise classification of slope states. Compared to traditional displacement mutation criteria, this model demonstrates higher accuracy in distinguishing stable/unstable boundaries for slopes under complex geological conditions (e.g., jointed or multilayered soils). In practice, the method incurs higher computational costs. However, unlike the iterative calculations of traditional methods, it enables parallel computing. When multiple devices are employed for simultaneous calculations, the required time may be further reduced.

(3) A novel definition of the comprehensive safety factor based on the shear strength ratio is proposed, explicitly reflecting the dynamic equilibrium between material strength reserves and stress states on potential slip surfaces. This provides an intuitive quantitative metric for engineering decision-making.

(4) Comparative validation with existing methods demonstrates the superiority of the proposed approach. For slopes under complex geological/engineering conditions, this method achieves more accurate classification of stable and unstable states, enhancing robustness while maintaining independence from specific reduction paths.

Future Prospects:

This study has preliminarily validated the effectiveness of the proposed method under relatively idealized geological conditions. It should be noted that the characteristics of the BGMM render it not amenable to conventional cross-validation methods. Meanwhile, the model has not been validated with actual data (e.g., field monitoring measurements). Therefore, the results are only applicable to conditions similar to the rock types and configurations considered in this study. To showcase its reliability and generalizability in real-world, complex geological scenarios, future research will focus on applying the method to a wider range of geologically more complex environments (e.g., those exhibiting strong heterogeneities, fault zones, or abrupt lithological changes). By conducting extensive case studies in these diverse and challenging geological settings, we will further elucidate the method’s robustness and potential for adaptation in addressing practical geological exploration challenges.

Furthermore, regarding the model analysis, while adopting the non-proportional reduction in dual parameters represents an improvement over the traditional proportional reduction method by better approximating the progressive failure mechanism of actual slopes (e.g., the inconsistency in the decay paths of cohesion and internal friction angle), it still assumes a uniform reduction in parameters across the entire slope mass. Real slope failures often involve the strain-softening behavior of the material and the development of localized failure zones, exhibiting strong spatial heterogeneity. Future research will incorporate softening constitutive models for the material and, based on this, explore non-uniform strength reduction states throughout the entire slope. This will be combined with the energy evolution analysis and Bayesian Mixture Model approach proposed in this study to establish a physically more realistic and more predictive framework for slope stability analysis.

Author Contributions

All authors contributed to the study conception and design. Conceptualization and methodology were developed by J.H. and T.D. Formal analysis and investigation were conducted by J.H., T.D. and S.P. The original draft was written by J.H., with critical revisions and editing by T.D. Funding acquisition was secured by X.Z. and resources were provided by X.P., D.W. and S.Z. The research was supervised by T.D. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Some or all data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

Thanks are extended to the editors and reviewers.

Conflicts of Interest

Authors Shouxing Peng, Xing Pang, Daochun Wan and Shaojun Zhang were employed by the company Pangang Group Mining. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

BGMM	Bayesian Gaussian Mixture Model
FOS	Factor of Safety

Footnotes

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

Figures and Tables

Figure 1 Schematic diagram of the Bayesian Gaussian Mixture Model.

Figure 2 Force diagram of the infinitesimal element.

Figure 3 Stress diagram on inclined plane in differential element.

Figure 4 Slope calculation model.

Figure 5 Three-dimensional relationship between reductions and dissipated energy.

Figure 6 Dissipated energy contours under varying reduction coefficients: (a) c:1, φ:1; (b) c:1, φ:1.2; (c) c:1.366, φ:1.368; (d) c:1.368, φ:1.366; (e) c:1.366, φ:1.37; (f) c:1.5, φ:1; (g) c:1.5, φ:1.1; and (h) c:1.5, φ:1.2.

Figure 7 Variation in dissipated energy under proportional reduction.

Figure 8 Gaussian component decomposition diagram.

Figure 9 Visualization of training results for the Bayesian Gaussian Mixture Model.

Figure 10 Comparison of data distributions.

Figure 11 Displacement diagram of critical slope No. 1.

Figure 12 Contour map of plastic zone distribution in the critical slope.

Table 1

Comparative analysis of clustering algorithms for slope energy data.

Feature	BGMM	K-Means	Hierarchical	DBSCAN
Core principle	Probabilistic density modeling	Distance minimization	Tree-based merging	Density expansion
Critical state identification	Handles overlapping states	Ignores transition fuzziness	Depend on subjective thresholds	Density-sensitive
Computation	High (Bayesian iteration)	Low (linear complexity)	Medium (matrix calculation)	Medium (Neighborhood search)
Slope data fit	Optimal for energy mutation	Limited for critical zones	Moderate for hierarchies	Moderate for uniform density

Table 2

Material parameters.

Cohesion/kPa	Internal Friction Angle/(°)	Elastic Modulus/MPa	Poisson’s Ratio	Density/(kg•m−3)
102	30.04	1850	0.27	2600

Table 3

Reduction coefficients and dissipated energy.

$\mathit{c}$-Reduction	$\mathit{\phi }$-Reduction	Dissipated Energy/J
1	1	241.4
1	1.1	241.4
1	1.2	241.5
…	…	…
1.368	1.366	3387.7
1.366	1.368	3572.2
1.366	1.37	3783.9
…	…	…
1.5	1	408
1.5	1.1	674.4
1.5	1.2	1272.2

Table 4

BGMM clustering outcomes.

$\mathit{c}$-Reduction	$\mathit{\phi }$-Reduction	Dissipated Energy/J	Membership Probability		Cluster Label
			Cluster 1	Cluster 2
1	1	241.4	0.980	0.020	Stable
1	1.1	241.4	0.980	0.020	Stable
1	1.2	241.5	0.980	0.020	Stable
…	…	…	…	…	…
1.368	1.366	3522	0.743	0.257	Stable
1.368	1.368	3801.5	0.564	0.436	Stable
1.368	1.37	4107	0.327	0.673	Unstable
…	…	…	…	…	…
1.5	1	408	0.983	0.017	Stable
1.5	1.1	674.4	0.986	0.014	Stable
1.5	1.2	1272.2	0.987	0.013	Stable

Table 5

Model evaluation metrics.

Log-Likelihood	Silhouette Score	Calinski–Harabasz Index	Davies–Bouldin Index	Weight Variability	Mean Separation
−0.69	0.88	916.52	0.64	0.36	2.02

Table 6

Critical slope data.

Group No.	$\mathit{c}$-Reduction	$\mathit{\phi }$-Reduction	Dissipated Energy
1	1.43	1.32	3890
2	1.372	1.366	3958.9
3	1.36	1.376	3968.9
4	1.37	1.368	3974.3
5	1.33	1.4	4037.5
6	1.38	1.36	4045.5
7	1.378	1.362	4077
8	1.368	1.37	4107
9	1.366	1.372	4145.2
10	1.376	1.364	4262.3
11	1.364	1.374	4387.7

Table 7

Shear strength calculation.

Group No.	$\mathit{c}$/kPa	$\mathit{\phi }$/(°)	${\mathit{\sigma }}_1$/kPa	${\mathit{\sigma }}_3$/kPa	$\mathit{\alpha }$/(°)	${\mathit{\sigma }}_{\mathit{n}}$/kPa	$\mathit{\tau }$/kPa
Initial Slope No. 0	102.00	30.04	359.60	90.79	60.02	157.26	193.24
Critical Slope No. 1	71.33	23.73	369.24	69.48	56.86	158.33	141.12
Critical Slope No. 2	74.34	23.01	377.09	72.03	56.51	164.19	144.29
Critical Slope No. 3	75.00	22.86	379.03	72.69	56.43	165.61	145.03
Critical Slope No. 4	74.45	22.98	377.54	72.24	56.49	164.55	144.45
Critical Slope No. 5	76.69	22.51	383.87	74.34	56.25	169.11	146.98
Critical Slope No. 6	73.91	23.10	376.26	71.76	56.55	163.53	143.88
Critical Slope No. 7	74.02	23.07	376.49	71.91	56.54	163.78	143.99
Critical Slope No. 8	74.56	22.95	378.36	72.51	56.48	165.06	144.67
Critical Slope No. 9	74.67	22.92	378.73	72.71	56.46	165.39	144.81
Critical Slope No. 10	74.13	23.04	377.82	72.47	56.52	164.65	144.37
Critical Slope No. 11	74.78	22.89	380.02	73.26	56.45	166.23	145.18

Table 8

Comparison of safety factors using different methods.

Method	FOS
Proposed method	1.338
Traditional SRM	1.371
Mean value method	1.369
Arithmetic square root	1.369
Mathematical fitting	1.369
Shortest path SRM	1.587

References

1. Dawson, E.M.; Roth, W.H.; Drescher, A. Slope stability analysis by strength reduction. Geotechnique; 1999; 49, pp. 835-840. [DOI: https://dx.doi.org/10.1680/geot.1999.49.6.835]

2. Guo, P.; Yuan, Y.; Peng, Y.; Zhang, H.; Zhao, Y. Analysis of Slope Stability and Disaster Law under Heavy Rainfall. Geofluids; 2021; 2021, 5520686. [DOI: https://dx.doi.org/10.1155/2021/5520686]

3. Xiao, S.; Dai, T.; Li, S. Review and Comparative Analysis of Factor of Safety Definitions in Slope Stability. Geotech. Geol. Eng.; 2024; 42, pp. 4263-4283. [DOI: https://dx.doi.org/10.1007/s10706-024-02793-6]

4. Nie, Z.; Zhang, Z.; Zheng, H. Slope Stability Analysis Using Convergent Strength Reduction Method. Eng. Anal. Bound. Elem.; 2019; 108, pp. 402-410. [DOI: https://dx.doi.org/10.1016/j.enganabound.2019.09.003]

5. Zhao, S.; Zheng, Y.; Shi, W. Analysis on Safety Factor of Slope by Strength Reduction FEM. Chin. J. Geotech. Eng.; 2002; 3, pp. 343-346. (In Chinese)

6. Peng, X.; Li, J.; Li, Y.; Jiang, J. Slope Instability Judgment Criteria in FEM Based on Strength Reduction Method. IOP Conf. Ser. Earth Environ. Sci.; 2020; 571, 012104. [DOI: https://dx.doi.org/10.1088/1755-1315/571/1/012104]

7. Sun, W.; Wang, G.; Zhang, L. Slope Stability Analysis by Strength Reduction Method Based on Average Residual Displacement Increment Criterion. Bull. Eng. Geol. Environ.; 2021; 80, pp. 4367-4378. [DOI: https://dx.doi.org/10.1007/s10064-021-02237-y]

8. Singh, S.K.; Kumar, S.; Chakravarty, D. Predicting the Stability of Rock Slopes in the Presence of Diverse Joint Networks and External Factors Using Machine Learning Algorithms. Min. Metall. Explor.; 2024; 41, pp. 2421-2440. [DOI: https://dx.doi.org/10.1007/s42461-024-01060-9]

9. Malekian, M.; Momayez, M.; Bellett, P.; Carra, F.; Tennakoon, E. Machine Learning for Slope Failure Prediction Based on Inverse Velocity and Dimensionless Inverse Velocity. Min. Metall. Explor.; 2023; 40, pp. 1557-1566. [DOI: https://dx.doi.org/10.1007/s42461-023-00781-7]

10. Chakraborty, R.; Dey, A. Probabilistic slope stability analysis: State-of-the-art review and future prospects. Innov. Infrastruct. Solut.; 2022; 7, 177. [DOI: https://dx.doi.org/10.1007/s41062-022-00784-1]

11. Guo, J. Stability Analysis of Bedding Rock Slope Based on Double Strength Reduction Method. J. Chongqing Jiaotong Univ. (Nat. Sci.); 2018; 37, pp. 61–66+72. (In Chinese)

12. Yuan, W.; Li, X.; Wang, W.; Bai, B.; Qi-Zhi, W.; Xiang-Jun, C. Study on Strength Reduction Method Based on Double Reduction Parameters. Rock Soil Mech.; 2016; 37, pp. 2222-2230. (In Chinese)

13. Wu, S.C.; Han, L.Q.; Li, Z.P.; Guo, C.; Zhou, J.X. Discussion on the Methods for Determining Slope Safety Factor Based on Stress State of the Sliding Surface. J. China Univ. Min. Technol.; 2018; 47, pp. 719-726. (In Chinese)

14. Lu, F.; Qiu, W. A Multiparameter Non-Proportional Shear Strength Reduction Method for Slope Stability Analysis Based on Energy Evolution Theory. Rock Soil Mech.; 2021; 42, 8.(In Chinese) [DOI: https://dx.doi.org/10.16285/j.rsm.2020.5435]

15. Su, X.; Wei, W.; Ye, W.; Meng, X.; Wu, W. Predicting Landslide Sliding Distance Based on Energy Dissipation and Mass Point Kinematics. Nat. Hazards; 2019; 96, pp. 1367-1385. [DOI: https://dx.doi.org/10.1007/s11069-019-03618-z]

16. Huang, L.; Huang, S.; Lai, Z. On an Energy-Based Criterion for Defining Slope Failure Considering Spatially Varying Soil Properties. Eng. Geol.; 2020; 264, 105323. [DOI: https://dx.doi.org/10.1016/j.enggeo.2019.105323]

17. Guo, Q.; Pan, J.; Cai, M.; Zhang, Y. Analysis of Progressive Failure Mechanism of Rock Slope with Locked Section Based on Energy Theory. Energies; 2020; 13, 1128. [DOI: https://dx.doi.org/10.3390/en13051128]

18. Zhang, E.; Zhou, B.; Yang, L.; Li, C.; Li, P. Failure Characteristic Analysis and Warning Method of Sandstone Based on Catastrophe Theory. Min. Metall. Explor.; 2023; 40, pp. 1865-1877. [DOI: https://dx.doi.org/10.1007/s42461-023-00830-1]

19. Yin, X.; Lin, H.; Chen, Y.; Tang, Y.; Wang, Y.; Zhao, Y.; Yong, W. Stability Analysis of Slope Considering the Energy Evolution of Locked Segment. Geotech. Geol. Eng.; 2022; 40, pp. 1729-1738. [DOI: https://dx.doi.org/10.1007/s10706-021-01989-4]

20. Gao, W. Stability Analysis of Rock Slope Based on an Abstraction Ant Colony Clustering Algorithm. Environ. Earth Sci.; 2015; 73, pp. 7969-7982. [DOI: https://dx.doi.org/10.1007/s12665-014-3956-4]

21. Zhang, W.; Liu, S.; Wang, L.; Sun, W.; He, Y.; Wang, Y.; Sun, G. The Overall Stability of a Partially Unstable Reservoir Bank Slope to Water Fluctuation and Rainfall Based on Bayesian Theory. Landslides; 2024; 21, pp. 2021-2032. [DOI: https://dx.doi.org/10.1007/s10346-024-02250-8]

22. Fan, C.; Zhang, N.; Jiang, B.; Liu, W.V. Weighted Ensembles of Artificial Neural Networks Based on Gaussian Mixture Modeling for Truck Productivity Prediction at Open-Pit Mines. Min. Metall. Explor.; 2023; 40, pp. 583-598. [DOI: https://dx.doi.org/10.1007/s42461-023-00747-9]

23. Wan, H.; Wang, H.; Scotney, B.; Liu, J. A Novel Gaussian Mixture Model for Classification. Proceedings of the 2019 IEEE International Conference on Systems, Man and Cybernetics (SMC); Bari, Italy, 6–9 October 2019; IEEE: Piscataway, NJ, USA, 2019; pp. 3298-3303. [DOI: https://dx.doi.org/10.1109/SMC.2019.8914215]

24. Rukhaiyar, S.; Ramulu, M.; Choudhury, P.B.; Pradeep, G.; Singh, P.K. Uncertainty Estimation in Blast Vibration Attenuation Model Using Bayesian Probabilistic Approach. Min. Metall. Explor.; 2021; 38, pp. 2405-2415. [DOI: https://dx.doi.org/10.1007/s42461-021-00492-x]

25. Deng, D.; Li, L.; Zhao, L. Limit Equilibrium Method (LEM) of Slope Stability and Calculation of Comprehensive Factor of Safety with Double Strength-Reduction Technique. J. Mt. Sci.; 2017; 14, pp. 2311-2324. [DOI: https://dx.doi.org/10.1007/s11629-017-4537-2]

26. Pinheiro, J.C.; Bates, D.M. Approximations to the Log-Likelihood Function in the Nonlinear Mixed-Effects Model. J. Comput. Graph. Stat.; 1995; 4, pp. 12-35. [DOI: https://dx.doi.org/10.1080/10618600.1995.10474663]

27. Layton, R.; Watters, P.; Dazeley, R. Evaluating Authorship Distance Methods Using the Positive Silhouette Coefficient. Nat. Lang. Eng.; 2013; 19, pp. 517-535. [DOI: https://dx.doi.org/10.1017/S1351324912000241]

28. Lima, S.P.; Cruz, M.D. A Genetic Algorithm Using Calinski-Harabasz Index for Automatic Clustering Problem. Rev. Bras. Comput. Apl.; 2020; 12, pp. 97-106. [DOI: https://dx.doi.org/10.5335/rbca.v12i3.11117]

29. Ros, F.; Riad, R.; Guillaume, S. PDBI: A Partitioning Davies-Bouldin Index for Clustering Evaluation. Neurocomputing; 2023; 528, pp. 178-199. [DOI: https://dx.doi.org/10.1016/j.neucom.2023.01.043]