Accurate spatial mapping of soil erodibility (K) is essential for assessing erosion risks and formulating conservation strategies. However, existing empirical models and spatial prediction face challenges, including underestimating spatial variability, static local environmental associations, and limited regional adaptability. This study proposed an innovative framework integrating empirical models, the second dimension of spatial association (SDA, incorporating multi-scale neighborhood features), and machine learning to select the optimal K values mapping. First, based on soil surveys and laboratory analyses in the Northeast China Black Soil (Mollisols) Region, three empirical models: the erosion-productivity impact calculator (K_EPIC), the Shirazi (K_Shirazi), and the Torri (K_Torri), were used to calculate K values. Second, SDA reconstructed environmental covariates by extracting quantile features (0-1) within radius-defined neighborhoods (100-3000 m), capturing multi-scale spatial dynamics. Third, Random Forest (RF) and Gradient Boosting Decision Tree (GBDT) were employed for digital mapping, while per configuration generated 40 ensemble models (10 random seeds & times; 4-fold cross-validation) to enhance model robustness. Results demonstrated that SDA-based models improved R 2 by 12-89 % compared to conventional static local association models. Considering data distribution, model accuracy, and spatial prediction, K_Shirazi demonstrated optimal regional representation (R 2 =0.4562, in SDA-GBDT). Moreover, climate and landscape are important driving factors for K_EPIC and K_Shirazi, while topography additionally influences K_Torri. The proposed framework offers a scientific and effective method to create the optimal pathway for soil erodibility mapping through multi-scale environmental feature extraction and integrated machine learning modeling, which could be transferred to other regions or fields.