Absolute principal component score/multiple linear regression (APCS/MLR) and positive matrix factorization (PMF) were applied to a dataset consisting of 10 heavy metals in 300 surface soils samples. Robust geostatistics were used to delineate and compare the factors derived from these two receptor models. Both APCS/MLR and PMF afforded three similar source factors with comparable contributions, but APCS/MLR had some negative and unidentified contributions; thus, PMF, with its optimal non-negativity results, was adopted for source apportionment . Experimental variograms for each factor from two receptor models were built using classical Matheron's and three robust estimators. The best association of experimental variograms fitted to theoretical models differed between the corresponding APCS and PMF-factors. However, kriged interpolation indicated that the corresponding APCS and PMF-factor showed similar spatial variability. Based on PMF and robust geostatistics, three sources of 10 heavy metals in Guangrao were determined. As, Co, Cr, Cu, Mn, Ni, Zn, and partially Hg, Pb, Cd originated from natural source. The factor grouping these heavy metals showed consistent distribution with parent material map. 43.1% of Hg and 13.2% of Pb were related to atmosphere deposition of human inputs, with high values of their association patterns being located around urban areas. 29.6% concentration of Cd was associated with agricultural practice, and the hotspot coincided with the spatial distribution of vegetable-producing soils. Overall, natural source, atmosphere deposition of human emissions, and agricultural practices, explained 81.1%, 7.3%, and 11.6% of the total of 10 heavy metals concentrations , respectively. Receptor models coupled with robust geostatistics could successfully estimate the source apportionment of heavy metals in soils.