Abstract: The ground motion intensity parameters and structural response parameters at different sites of an earthquake event are spatially correlated, and the related seismic effects will lead to the sharp concentration and accumulation of earthquake losses, resulting in catastrophic events. However, the traditional regional seismic risk assessment does not consider the spatial correlation of ground motion intensity measures. In this paper, geo-statistical methods are used to characterize the spatial correlation of ground motion intensity measures in earthquake events. A geometrically robust estimation method of theoretical semi-variogram is proposed, and the weighted least squares method is used to accurately and efficiently evaluate the short-distance spatial correlation. Taking the Chuetsu-Oki earthquake as an example, based on the exponential semi-variogram model, the weighted least squares method is used to fit the geometrically robust estimation of semi-variogram value, and the spatial correlation functions of PGA, PGV and 0-10 s spectral acceleration are obtained. A prediction model of the range b value is proposed. The results show that the spatial correlation decay rate of the long-period spectral acceleration is smaller than that of the short-period spectral acceleration, and there is an inflection point at T=0.45 s. The developed spatial correlation model is applied to regional probabilistic seismic hazard analysis. The Monte Carlo method is used to simulate the spatially correlated random field of ground motion intensity measures of a hypothetic region. The traditional regional seismic hazard analysis method is refined. It is found that the annual exceedance probability considering spatial correlation is greater than that without considering spatial correlation in most cases, and it further shows that the consideration of spatial correlation has important guiding significance on risk assessment of portfolio buildings and infrastructures as well as the construction of resilient cities.