The accurate identification of spatial drivers is crucial for effectively managing soil heavy metals (SHM). However, understanding the complex and diverse spatial drivers of SHM and their interactive effects remains a significant challenge. In this study, we present a comprehensive analysis framework that integrates Geodetector, CatBoost, and SHapley Additive exPlanations (SHAP) techniques to identify and elucidate the interactive effects of spatial drivers in SHM within the Pearl River Delta (PRD) region of China. Our investigation incorporated fourteen environmental factors and focused on the pollution levels of three prominent heavy metals: Hg, Cd, and Zn. These findings provide several key insights: (1) The distribution of SHM is influenced by the combined effects of various individual factors and interactions within the source–flow–sink process. (2) Compared with the spatial interpretation of individual factors, the interaction between Hg and Cd exhibited enhanced spatial explanatory power. Similarly, interactions involving Zn mainly demonstrated increased spatial explanatory power, but there was one exception in which a weakening was observed. (3) Spatial heterogeneity plays a crucial role in determining the contributions of environmental factors to soil heavy metal concentrations. Although individual factors generally promote metal accumulation, their effects fluctuate when interactions are considered. (4) The SHAP interpretable method effectively addresses the limitations associated with machine-learning models by providing understandable insights into heavy metal pollution. This enables a comparison of the importance of environmental factors and elucidates their directional impacts, thereby aiding in the understanding of interaction mechanisms. The methods and findings presented in this study offer valuable insights into the spatial heterogeneity of heavy metal pollution in soil. By focusing on the effects of interactive factors, we aimed to develop more accurate strategies for managing SHM pollution.