After an earthquake occurs, it is a key technology to quickly generate a reasonable isoseismal map using strong motion observation data. The intensity isoseismal map generated based on source parameters and intensity attenuation law is too regular to reflect the local intensity distribution characteristics. However, the isoseismal map directly generated from the observed data by interpolation cannot reasonably reflect the source effect of large earthquakes. This article proposes an intensity field fusion Kriging interpolation method that considers the spatial correlation of seismic motion and the source effect of large earthquakes. The method utilizes the distance between points to consider the spatial correlation of intensity, and the distance between points to source to consider the source effect of large earthquakes. First, in the framework of Kriging interpolation method, Particle Swarm Optimization algorithm was used to realize the fast automatic fitting of semivariogram. Secondly, the consideration of spatial correlation and source effects was achieved through fitting two semivariogram functions with different independent variables. Finally, the interpolation results of earthquake data from Japan and China were used to verify the rationality of the method. The results show that the shape of the intensity isoseismal map obtained by the fusion Kriging interpolation method is regular, which can reflect the source effect of an earthquake, and has a good adaptability to the sparse situation of stations, meeting the requirements of rapid seismic intensity reporting. This method can provide technical support for rapid seismic intensity reporting and the drawing of seismic intensity isoseismal map.