## 1 Introduction  
The analysis of slope stability continues to be a pivotal area of study in geotechnical engineering. While most current analyses are confined to two-dimensional (2D) models, real-world slope failures often exhibit pronounced 3D characteristics [1, 2]. Factors such as the failure width of landslide bodies, the intricate geometry of slopes, and heterogeneous geological conditions contribute to slope instability in 3D contexts. Simplified 2D analyses, due to their inherent limitations in precision, may lead to overly conservative slope designs, resulting in avoidable expenditures [3, 4]. Consequently, scholars are progressively transitioning from 2D to 3D in addressing slope stability challenges. In conventional slope stability assessments, the soil is typically treated as a homogeneous medium, with soil strength parameters averaged from statistical data to determine the safety factor, a measure of slope stability. However, soil parameters display spatial variability due to material composition, depositional history, and consolidation pressures [5; 6]. Recognizing the limitations of a single safety factor in comprehensively quantifying slope reliability, researchers have integrated reliability analysis techniques [7]. These methods offer a significant advantage by accounting for various factors, including the interdependence of geotechnical parameters, statistical uncertainties, and the unpredictability inherent in computational models. The failure probability derived from reliability analysis more accurately reflects the actual safety status of slopes, and its combination with safety factors enhances the precision of stability evaluations [8; 9]. Recent advances in reliability theory have led to the development of diverse methodologies, including the First-Order Reliability Method (FORM) [10; 11], Monte Carlo Simulation (MCS) [12; 13], Response Surface Method (RSM) [14; 15], Subset Simulation (SS) [16; 17], and Polynomial Chaos Expansions (PCE) [18; 19]. Among these, MCS has gained prominence due to its conceptual simplicity and compatibility with advancements in computational power. Despite its ability to yield accurate slope reliability assessments, the computational demands of direct MCS can be substantial. To mitigate this, hybrid approaches combining MCS with other techniques have emerged. For instance, [PERSON] et al. (2009) [7] integrated the elastoplastic finite element method and random field theory within the MCS framework to investigate the failure probabilities of both 2D and 3D slopes. [PERSON] et al. (2021) [14] employed a Stochastic Response Surface Method combined with variance reduction techniques (VRSRSM) to alleviate the computational load associated with implementing spatially correlated input parameters in 3D spaces. [PERSON] et al. (2017) [20] developed a Multiple Random Response Surface Method (MRRSM) to assess the systemic reliability of multi-layered soil slopes with spatially variable soil, streamlining the process by constructing multiple random response surfaces on the Representative Sliding Surface (RSS). These methods have reduced the simulation time for critical sliding surfaces while maintaining accuracy, demonstrating the potential for balancing computational efficiency and precision in slope reliability analysis. However, a balance must be struck between computational precision and efficiency in the modeling of critical sliding surfaces. To reconcile this dichotomy, heuristic optimization algorithms have been increasingly adopted [21; 22; 23; 24; 25]. Among these, the Particle Swarm Optimization (PSO) algorithm stands out for its user-friendly operation, high computational efficiency, and robust performance [26]. This algorithm has been successfully applied to determine the Critical Slip Surface (CSS) for homogeneous soil slopes [27], slopes under earthquake forces [28], and 3D asymmetric slopes [29]. Notably, [PERSON] et al.