Abstract: The failure probability of a soil slope system is strongly related to soil shear strength parameters and groundwater level. In the traditional slope reliability analysis, only the variability of soil shear strength parameters is generally considered, and the influence of random groundwater level is ignored. Finite element discrete technique, upper bound method theory, relevant non-Gaussian random field simulation method and, stochastic planning theory are combined to study the failure probability of the slope system under the effect of stochastic groundwater level considering the spatial variability of parameters. A stochastic programming model for slope reliability analysis is constructed upon the upper bound method for plastic limit analysis by using a finite element discrete slope soil body, taking the groundwater level of the slope as a random variable and considering the spatial variability of the soil shear parameters, and then the model is employed to obtain the stability safety coefficient and velocity field of the slope upon Monte Carlo simulation; the AP clustering analysis method is used for the first time to the system failure probability of the slope estimated based on the failure information of the elements in the slope, and the reliability analysis theory of the slope system is developed. A classical slope calculation case is systematically analyzed and compared with the results of Limit Equilibrium Analysis Method and Finite Element Method. The results show that: there are multiple failure modes of slopes under the action of random groundwater level considering the spatial variability of parameters, and the AP cluster analysis method can identify all the failure modes of slopes and the corresponding failure probabilities based on the location information of failed elements; an efficient parallel program of upper bound method is prepared upon Matlab, which greatly improved computing efficiency.