Ground‐motion prediction equations (GMPEs), also called ground‐motion models and attenuation relationships, are empirical models widely used in probabilistic seismic hazard analysis (PSHA). They estimate the conditional distribution of ground shaking at a site given an earthquake of a certain magnitude occurring at a nearby location. In the past decade, the increasing interest in assessing earthquake risk and resilience of spatially distributed portfolios of buildings and infrastructure has motivated the modeling of ground‐motion spatial correlation. This introduces further challenges for researchers to develop statistically rigorous and computationally efficient algorithms to perform ground‐motion model estimation with spatial correlation. To this goal, we introduce a one‐stage ground‐motion estimation algorithm, called the scoring estimation approach, to fit ground‐motion models with spatial correlation. The scoring estimation approach is introduced theoretically and numerically, and it is proven to have desirable properties on convergence and computation. It is a statistically robust method, producing consistent and statistically efficient estimators of inter‐ and intraevent variances and parameters in spatial correlation functions. The performance of the scoring estimation approach is assessed through a comparison with the multistage algorithm proposed by Jayaram and Baker (2010) in a simulation‐based application. The results of the simulation study show that the proposed scoring estimation approach presents comparable or higher accuracy in estimating ground‐motion model parameters, especially when the spatial correlation becomes smoother. The simulation study also shows that ground‐motion models with spatial correlation built via the scoring estimation approach can be used for reliable ground‐shaking intensity predictions. The performance of the scoring estimation approach is further discussed under the ignorance of spatial correlation, and we find that neglecting spatial correlation in ground‐motion models may result in overestimation of interevent variance and underestimation of intraevent variance and thus inaccurate predictions.