Stability evaluation remains a research hotspot in slope engineering, yet upper bound limit analysis (UBLA) of slopes simultaneously considering stratum distribution and slope morphology are scarce. Taking multi-layer and multi-stage soil slopes as research objects, this study constructs three kinematically admissible failure mechanisms: global rotational, local rotational, and rotational-translational failure along soil interface. Then, the deduced integral solutions of external work rate and internal energy dissipation rate are solved using the composite Simpson’s rule. The strength reduction method is applied to solve the safety factor (
Fs) in conjunction with a double exhaustion search strategy. The proposed method is verified by the Miaoling 750kV substation high slope, where the error between the calculated
Fs and the UBLA solution obtained from OptumG2 is less than 1%. Furthermore, the influence of geometric and shear strength parameters on the critical sliding surface and slope safety factor is discussed. The results show that multiple failure modes exist in multi-layer and multi-stage soil slopes, and both geometric and shear strength parameters affect the failure mechanism. A typical software for UBLA of two-layer and three-stage soil slopes is developed using MATLAB App Designer. This study provides valuable insights for understanding landslide mechanisms and optimizing slope design.