The relationship between urban heat island and underlying surface was analyzed for the first time in Beijing，Shanghai，Shenyang and Wuhan in summer based on Landsat ETM+ data by combining mathematical statistics with spatial statistics. The urban underlying surface was represented by green vegetation and impervious surface generated by a linear spectral mixture analysis，while the relative urban heat island strength was represented by difference of radiant brightness temperature between surface and water. Radiant brightness temperature was retrieved from the thermal band. Results obtained in these cities suggest that a positive spatial autocorrelation exists in the distribution of green vegetation，impervious surface and heat island strength，with the existence of a similar autocorrelation range equivalent to the spatial pattern character of the combination of buildings and roads. As a result of spatial autocorrelation，the structural factor is more important than the random factor in contributing to spatial heterogeneity of urban green vegetation，impervious surface and heat island strength. Statistics show that the positive relationship between impervious surface and heat island strength exhibits different nonlinear characteristics in these four cities，which can be expressed by piecewise functions according to the quantity of impervious surface，and things are quite similar for the negative relationship between green vegetation and heat island strength. In these four cities，both the positive and the negative relationships demonstrate spatial heterogeneity and multi-scale character，along with an obvious spatial scale proximate to 550 meters detected by using cross-correlation. These conclusions may provide some ideas for civil planning，i.e.，through designing of reasonable intervals and proportions between green vegetation zones and buildings，the effect of urban heat island can be weakened.