## 1 Introduction  
The instability of slope can not only endanger people's lives and property but also lead to great loss of national economy [1; 2; 3]. Under the long-term geological effects of deforestation, deposition and post-deposition, the shear strength parameters of slope soil mass show spatial variability. Research shows that the spatial variability of soil parameters is a major source of uncertainty in slope engineering. It not only affects bearing capacity and safety performance, but also has a significant influence on the slope failure mode [4; 5; 6]. Different failure modes correspond to different failure consequences. The traditional integrated failure probability (IFP) method based on the slope stability coefficient for failure risk analysis needs to count the failure modes individually. In addition, the calculation of the slope stability coefficient involves complex processes such as limit state definition and strength reduction. Therefore, it is urgent to develop new methods for slope failure risk analysis. At present, the three most commonly used methods for slope stability analysis are the rigid body limit equilibrium method (LEM) [7; 8; 9], the finite element method (FEM) [10; 11; 12] and the limit analysis method (LAM) [13; 14; 15; 16; 17; 18; 19]. The LEM first assumes the sliding surface, then divides the soil on the sliding surface into strips and turns a hyper-static problem into a static problem by assuming the magnitude and direction of the inter-strip force. Based on different assumptions, a variety of LEMs have been formed, such as the Ordinary/Fellenius method, the Bishop Simplified method, the Morgenstern-Price method and the Spencer method. The LEM is widely used in engineering because of its clear concept and simple calculation. However, it possesses the following drawbacks when using this method for slope stability analysis considering spatial variability of parameters: (1) only considering the spatial variability of soil material in one dimension on the sliding surface, ignoring the influence outside the sliding surface [20]. Due to the deficiencies of LEM, the stress-strain relationship of soil cannot be accurately obtained [18; 20]. Based on the principle of minimum potential energy, the FEM turns differential equations into linear algebraic equations to solve and simultaneously satisfies the constraints of equilibrium equations, geometric equations and constitutive relations of discrete triangular elements. LAM applies the upper and lower limit theorems generally observed by an ideal elastoplastic body and obtains the upper and lower limit solutions when the slope is in the limit state by constructing kinematically-admissible velocity fields and static admissible stress fields. It is known from the upper and lower limit theorems that the true solution must lie between the upper and lower limit solutions. This method can efficiently and accurately obtain the ultimate load [18; 19], stability coefficient [13; 14; 15; 16; 17; 18] and failure mode [5; 18] of a slope in the limit state without considering the complex constitutive relationship of rock and soil mass. [PERSON] et al. [18] and [PERSON] et al. [19] set the shear strength parameters of soil as random variables and solved the upper and lower limit solutions. [PERSON] et al. [21] studied the failure probability of slopes under the action of random groundwater levels from the perspective of element failure probability (EFP). [PERSON] et al. [20] found that the spatial variability of soil shear strength parameters has a significant influence on the slope failure mode. The above studies mostly focus on the slope failure probability and pay less attention to the failure consequences, which is not necessarily. For example, when a slope with multi-failure modes fails, different failure modes (shallow landslide/deep landslide) correspond to different failure consequence coefficients (sliding mass volume).