The geomechanical parameters for a particular site exhibit inherent uncertainties due to geological processes, and probabilistic back analysis incorporating field observation data can effectively reduce these uncertainties. Although the BUS (Bayesian Updating with Subset simulation) method can transform the high-dimensional probabilistic back analysis problem with the equality site information into an equivalent structural reliability problem, the value of the constructed likelihood function can become extremely small or even lower than the computer floating-point operation accuracy as the field observation data increase, which might seriously affect the computational efficiency and accuracy of probabilistic back analysis. To this end, this paper proposes an improved BUS method based on the parallel system reliability analysis. Starting from the Cholesky decomposition-based midpoint method, the total failure domain with a low acceptance rate is decomposed into several sub-failure domains with a high acceptance rate so as to avoid the “curse of dimensionality” arising from the integration of a large amount of field observation data, and to achieve accurate back analysis of the geomechanical parameters of slopes. Finally, the effectiveness of the proposed method is validated through a case study of an undrained saturated clay slope. The results show that the proposed method can integrate a large number of borehole data and the observation information of slope service state for efficient probabilistic back analysis of geomechanical parameters and slope reliability evaluation with reasonable accuracy. The proposed method provides an effective tool for high-dimensional probabilistic back analysis of spatially variable soil parameters and slope reliability evaluation.