<p>Because of the particularity of the geometric shape, the failure modes of multistage slope show diversity. In order to solve the problem of failure mode classification, a new multistage slope failure mode classification method based on system clustering method was proposed for the first time. Firstly, the multistage slope is discretized by noncommon node triangular elements, and on this basis, the mathematical programming model of the multistage slope with upper and lower bound method is established. By solving the model, the upper and lower bound solutions of the safety factor of the multistage slope, as well as the corresponding kinematically admissible velocity field and statically admissible stress field are obtained. Meanwhile, a method for calculating the similarity of slope failure modes based on the upper bound method is proposed, which realizes the quantitative analysis of failure modes. To verify the effectiveness of the proposed method, systematic analysis is carried out on two typical multistage slope engineering cases, and the analysis results are compared with those obtained by the finite element method and limit equilibrium analysis method. The results show that the stability analysis results of multistage slope obtained by the upper and lower bound method proposed in this paper are basically consistent with those obtained by the finite element method and limit equilibrium analysis method; the method can identify multiple failure modes in the tested multistage slope examples without predefining slip surfaces. The method proposed in this study has important theoretical reference value and practical engineering application significance for the research on failure modes of multistage slope.</p>