This study aimed to identify the spatial patterns of potentially toxic elements (PTEs), including the spatial distribution, spatial autocorrelation, and risk probability, and to quantify the sources of PTEs, to provide guidelines for soil management. Spatial distributions and probabilities of PTEs were determined by empirical Bayesian kriging (EBK), while spatial autocorrelation was estimated by Moran's I. Positive matrix factorization (PMF) was adopted for the quantitative source contributions of PTEs. More than 64.6% of Co, Cr, Mn, and Ni were derived from geogenic sources, with high regions and high-high clusters both correlated to sandstone. Thus, it can be deduced that parent materials dominated the spatial patterns of these PTEs. In addition, some hotspots were situated in urban areas, and the influence of human activities on these four PTEs should be considered. Industry-traffic discharge and parent materials both influenced As, Cu, Pb, and Zn. Nonetheless, the spatial patterns of As, Cu, Pb, and Zn were formed by anthropogenic emissions since hotspots and high-high clusters were contiguously situated in urban areas. 58.5% of Hg originated from atmospheric deposition related to industrial emissions, and 47.2% of Cd was controlled by the application of chemical fertilizers. High levels of Hg and Cd mainly corresponded with industrial sites and cultivated land, suggesting that industrial and geoponic production played major roles in the generation of spatial patterns for Hg and Cd, respectively. Furthermore, the Cd and Hg posed a severe risk to soils, with a high probability to surpass 1.5 times the backgrounds. The EBK, Moran's I, and PMF results showed that all ten PTEs were enriched to some degree due to natural or anthropogenic factors. The results of geostatistical analysis and the receptor model can be mutually verified, indicating the reliability of these methods.