A GIS-based geostatistical method was put forward for analyzing spatial variation of combined heavy metal pollution in soils. As a case study, an area, about 7.7 km², around the Jiangxi Copper Smeltery in a red soil hill region in Guixi County, Jiangxi Province was selected. The sequence of the analysis method goes as follows: a) collect soil, plant or water samples by random grid from 50m to 200m with DGPS, and analyze the samples for heavy metal contents and other properties; b) explore distribution pattern and spatial variation of the data, then transform the data to get rid of the trend; c) calculate semivariogram, and select the best model and appropriate parameters; d) select the best kriging method based on cross validation; e) interpolate soil heavy metal contents and draw a heavy metal pollution map, then analyze spatial distribution characteristics and their controlling factors. The results show that there is combined Copper and Cadmium pollution in the study area, and the above method based on the geostatistical module in ArcGIS can exactly predict the spatial distribution pattern of combined heavy metal pollution in the hill region.