A method for slope system reliability analysis considering spatial variability of soil is developed to rationally evaluate the reliability of slope stability. The basic idea is that the slope domain is discretized into a series of “elements”, within each of which the soil property is fully correlated. The correlation coefficient among soil properties from different elements is obtained assuming autocorrelation function. Sets of random samples are generated by subset simulation, and harmony search algorithm is adopted to search the minimum factor of safety as the corresponding output. The method is applied to two examples to analyze the reliability and then compare it with existing results. Overall, the results show that the failure probability is lower than the case without considering the spatial variability and the potential failure modes are greater than that without taking into account the spatial variability.