The traditional slope reliability methods may not efficiently carry out slope reliability considering inherent spatial variability (ISV) of soil properties. This paper aims to develop a two-stage dimension reduction method for efficient slope reliability analysis considering ISV of soil properties. First, the framework of slope reliability analysis based on surrogate model is provided. Second, the two-stage dimension reduction method is proposed systematically. Thereafter, the implementation procedure for slope reliability analysis using the proposed method is summarized. The validity of two-stage dimension method is illustrated with a multiple-layered soil slope. The results indicate that the proposed two-stage dimension method can reduce the number of original model evaluations significantly and improve the efficiency of slope reliability analysis considerably. The surrogate model based on the proposed method is proved to be unbiased. The loss of accuracy in modeling ISV has a slight influence on the minimum factor of safety, sliding path and sliding volumes of slope. The proposed method can yield sufficiently accurate reliability results, which is more efficient than direct Monte Carlo simulation. For the multi-layered soil slope reliability considering ISV of soil properties examined in this study, the proposed method performs well in terms of accuracy and efficiency.