<p>For system reliability analysis of multi-layered slopes using the limit equilibrium method (LEM), the numerous highly nonlinear potential slip surfaces (PSSs) and the computationally intensive deterministic models pose significant challenges in accurately and efficiently estimating failure probability. Although the Bayesian compressive sensing (BCS)-based response surface method (RSM) has demonstrated its ability to construct high-accuracy response surfaces for highly nonlinear problems, directly constructing a response surface between the minimum factors of safety (FSs) corresponding to critical slip surfaces and input variables still necessitates a substantial amount of sampling data to guarantee construction accuracy. To specifically address this issue, this study develops an efficient multiple BCS (MBCS)-based RSM to enhance the computational efficiency for LEM-based multi-layered slope stability using sparse sampling data. In the proposed MBCS-based RSM, a few representative slip surfaces (RSSs) are identified from PSSs as key failure modes of the slope system, and a unique BCS-based response surface is constructed for each RSS, thereby circumventing the challenge posed by the high nonlinearity in constructing surrogate model between the minimum FSs and the input variables. Moreover, an innovative BCS-based RSM using sliced inverse regression is developed to address high-dimensional system reliability analysis problems that arise when representing the spatial variability of soil parameters using random field models. Investigations using three multi-layered slope reliability analysis problems indicate that the proposed MBCS-based RSM performs admirably.</p>