This paper presents a new algorithm for assessing the reliability of three-dimensional (3D) slope stability considering the spatial variability of soil based on the Particle Swarm Optimization (PSO) algorithm. First, a 3D random field is generated using the Karhunen–Loève (K-L) expansion method. Then, the simplified Bishop method of limit equilibrium is coupled with the PSO algorithm to calculate safety factors of the slope. Finally, the failure probability of the slope is determined using the Monte Carlo Simulation method. After validating the rationality of the proposed method through a typical case study, this paper offers an in-depth examination of how soil spatial variability affects the stability of 3D slopes. It is observed that, given identical soil correlation lengths, slope geometric parameters, and failure surface widths, the failure probability is positively correlated with soil spatial variability parameters, while the mean safety factor demonstrates an inverse relationship with these variability parameters. Additionally, the failure probability tends to increase as the soil correlation lengths increase, and it also escalates with the expansion of the failure surface width. In contrast, the mean safety factor exhibits an upward trend with the augmentation of the horizontal correlation length, while it diminishes progressively as the vertical correlation length grows, and it also shows a decline with the widening of the failure surface width. The proposed algorithm significantly improves computational efficiency while ensuring accuracy, making it suitable for the reliability analysis of three-dimensional slopes.