Inversions of geodetic data are regularly used to estimate interseismic locking in subduction zones. However, the ill‐posed nature of these problems motivates us to include prior information, physically consistent with processes of the subduction seismic cycle. To deal with model instabilities, we present an inversion method to estimate both plate‐locking and model uncertainties by inverting Global Navigation Satellite System derived velocities based on a Bayesian model selection scheme. Our method allows us to impose positivity constraints via a multivariate folded‐normal distribution, with a specified covariance matrix. Model spatial correlations are explored and ranked to find models that best explain the observed data and for a better understanding of locking models. This approach searches for hyperparameters of the prior joint multivariate probability density function (PDF) of model parameters that minimize the Akaike Bayesian Information Criterion (ABIC). To validate our approach, we invert synthetic displacements from analytic models, yielding satisfactory results. We then apply the method to estimate the plate‐locking in Central Chile (28°–39°S) and its relation to the coseismic slip distribution of earthquakes with magnitudes Mw > 8.0, on the subduction zone since 2010. We also search among different prior PDFs for a single ductile‐fragile limit depth. Our results confirms a spatial correlation between locked asperities and the 2010 Mw 8.8 Maule and 2015 Mw 8.3 Illapel earthquake rupture zones. The robustness of our locking model shows potential to improve future seismic and tsunami hazard estimations.