Efficient estimation and mapping of soil heavy metal concentrations based on rapid and environmentally friendly multispectral data can provide a reference for the prevention and control of non-point source pollution and safety in ecosystems. In this study, the spatial variation of soil heavy metals (As, Cd, Cr, Cu, Ni, and Zn) in Wuqing was analyzed with 394 and 57 soil samples gathered in 2005 and 2019. Based on these samples, as well as multispectral variables (raw bands and their respective vegetation index and derivative transformation), random forest-ordinary kriging (RFOK) and artificial neural network-ordinary kriging (ANNOK) were combined to predict spatial variations of heavy metal content. Results indicated that the mean values of soil As, Cd, Cr, Cu, and Zn in 2019 exceeded their respective concentrations in 2005 by 1.40, 1.36, 1.12, 1.12, and 1.35 times, and exceeded the standard rate of background values (BGV) in 2019 for all soil heavy metals (except for Ni). Compared to RF, ANN, and RFOK, the ANNOK method exhibited the best estimation performance for As (R20052= 69.47% andR20192= 65.03%), Cd (R20052= 64.48% andR20192= 64.20%), Cr (R20052= 72.02% andR20192= 70.65%), Cu (R20052= 72.05% andR20192= 70.53%), Ni (R20052= 69.37% andR20192= 66.03%), and Zn (R20052= 71.16% andR20192= 67.11%) in the validation dataset. Mapping of heavy metals derived from ANNOK revealed a similar distribution in 2005 and 2019, where the high concentrations of all soil heavy metals were in the southern part of the study area. The effects of multispectral bands as indicated by vegetation index and derivative transformation exhibited sensitivity for predicting soil heavy metals based on hybrid geostatistical methods. Overall, greater efficiency of spatial estimation for the content of soil heavy metals is important for the prevention and control of non-point source pollution.