The spatial variability is an inherent uncertainty of soils. The random field theory is used to represent the spatial variability of soils, and the random field discretization is performed by the Karhunen-Loève (KL) expansion method. Using the slope upper bound analysis based on the discrete mechanism, the discretization results of the internal friction angle random field at each point in the space are considered when generating the velocity discontinuity surface. And the strength reduction technique, bisection searching, and sequential quadratic programming method are combined to solve the safety factor of slopes. The first-order reliability method (FORM) and subset simulation (SS) are employed for slope reliability analysis. Given the characteristics of SS and the shear strength reduction technique, an optimization algorithm coupling the two is proposed to improve computational efficiency. By calculating and analyzing an earth slope, the similarities and differences between FORM and SS based on the KL expansion method in solving the slope reliability index and failure consequence are clarified. The influence of the coefficient of variation of soil strength parameters on the slope reliability index and failure consequence is investigated, providing a theoretical basis for risk analysis and prevention of slopes.