【Objective】The hilly region of western Chongqing is a transitional zone characterized by fragmented terrain and substantial elevation gradients, where soil thickness exhibits pronounced spatial variability driven by complex geomorphic processes. Soil thickness is a key indicator of soil resource status and plays a critical role in hydrological processes, ecosystem productivity, and geological hazard assessment. Traditional models used for soil thickness prediction often lack generalizability in hilly regions due to pronounced spatial heterogeneity, while pure machine learning approaches tend to overlook inherent spatial autocorrelation, potentially introducing local prediction biases. Therefore, accurate spatial prediction is essential for understanding soil resource distribution and supporting sustainable land management. 【Method】Based on 117 soil profiles (2020–2025) and 38 environmental covariates spanning topography, climate, remote sensing, parent material, and human activities, three benchmark models were compared: multiple linear regression (MLR), random forest (RF), and extreme gradient boosting (XGBoost) against a hybrid model, random forest residual Kriging (RFRK). A multi-stage variable selection procedure incorporating variance inflation factor, boruta algorithm, and recursive feature elimination addressed multicollinearity and overfitting. Climate variables (1 km) were resampled to 12.5 m via bilinear interpolation. After excluding five outliers (3σ/5σ criteria), 112 valid samples were partitioned into training (80%) and validation (20%) sets. Model performance was assessed using mean absolute error (MAE), root mean square error (RMSE), coefficient of determination (R2), and Lin’s concordance correlation coefficient (LCCC). For RFRK, variogram analysis characterized residual spatial structure, with model hyperparameters optimized through grid search. 【Result】(1) The screening strategy reduced 38 covariates to seven core variables, ranked by importance: total catchment area (TCA) > channel network baseline (CNBL) > soil enhancement index 2 (SER2) > channel network distance (CND) > soil enhancement index 1 (SER1) > aspect (ASP) > slope (SLP). This subset (five topographic, two remote sensing) underscores terrain-driven material redistribution as the primary control on soil thickness. (2) RFRK achieved the highest accuracy (R2 = 0.68, LCCC = 0.79, MAE = 14.59 cm, RMSE = 19.86 cm), significantly outperforming all benchmarks. Relative to the best benchmark (RF), RFRK improved R2 by 9.7%, reduced MAE by 12.3%, and reduced RMSE by 6.2%. MLR performed the poorest (R2 = 0.37), whereas XGBoost was intermediate (R2 = 0.56). (3) Variogram analysis revealed strong residual spatial autocorrelation (nugget/sill = 0.0017; range = 179.08 m), justifying kriging for local bias correction. The predicted map showed “thick northeast/southwest, thin central,” with thick layers (>80 cm) in valleys and gentle slopes, and thin layers (0–40 cm) on steep ridges and areas with intense anthropogenic disturbance. 【Conclusion】The RFRK hybrid model, integrating machine learning for nonlinear responses with geostatistical modeling of residual autocorrelation, characterizes soil thickness patterns in complex terrains more accurately than pure machine learning or linear models. This study validates the VIF-Boruta-RFE screening strategy and confirms RFRK’s superiority in capturing both deterministic environmental controls and stochastic spatial structures. The findings provide robust technical support for soil resource assessment and management and offer a generalizable methodological reference for digital soil mapping in comparable hilly and mountainous landscapes worldwide.