Surface subsidence caused by non-metal mining exhibits significant non-stationarity, and traditional Kriging methods struggle to effectively integrate multi-source geological information. This paper proposes a Bayesian Kriging with External Drift (BKED) model, which incorporates ore body depth, mined thickness and the comprehensive strength of the overlying strata as external drift variables, and employs the Markov Chain Monte Carlo (MCMC) method to quantify parameter and prediction uncertainties. Taking a phosphate mine in Yunnan as a case study, leave-one-out cross-validation was performed based on data from 27 surface monitoring points. The results indicate that the BKED model achieves a 90% prediction interval coverage of 92.6%, with an average CRPS reduced by 22.6% compared to conventional kriging, and a point prediction RMSE of 15.7 mm. Mined thickness is the dominant factor controlling subsidence, with a Sobol total effect index of 0.727. The BKED model provides reliable probabilistic prediction intervals whilst maintaining point prediction accuracy, making it suitable for mine subsidence risk assessment.