Soil heavy metals (HMs) pollution has caused significant land degradation by affecting soil properties and functions. Clarifying sources and spatial distributions of soil HMs is necessary for the investigation of land degradation, but lacks accurate and efficient methods. This study proposes a combined method for improving source apportionment and spatial prediction of soil HMs. Finite mixture distribution modelling (FMDM) was used to explore the backgrounds and contamination thresholds, and can verify the source contribution calculated by positive matrix factorization (PMF). To improve the efficiency of traditional multivariate geostatistical simulation, independent component analysis and sequential Gaussian simulation were integrated into geostatistical independent simulation for revealing the spatial patterns and the polluted areas. The combined method was applied to an HMs dataset of Longkou City, eastern China as a case-study. Both PMF and FMDM produced three sources of HMs. Cr and Ni were dominated by natural sources. Agricultural inputs mainly affected Cd, Cu, and Zn concentrations. Hg, 66.8%, and Pb, 8.1%, originated from industrial, mining, and traffic emissions. Hazardous areas exceeding moderate pollution cutoffs comprised 203.88 km 2 , accounting for 23.2% of the study area. Less than 0.1% of the areas for all HMs were identified as hazardous zones exceeding the high pollution thresholds. The reasonable management of industrial and agricultural activities, traffic emissions, and mining could have an important influence on the protection of land resources. Combining FMDM, PMF, and geostatistical independent simulation proved valid for investigating the land degradation due to HMs pollution.