We propose a new Bayesian approach to estimate continuous crustal strain-rate fields from spatially discrete displacement-rate data, based on Global Navigation Satellite System (GNSS) observations, under the prior constraint on spatial flatness of the strain-rate fields. The optimal values of the hyperparameters in the model of strain-rate fields are determined by using Akaike's Bayesian Information Criterion. A methodological merit of this approach is that, by introducing a two-layer Delaunay tessellation technique, the time-consuming computation of strain rates can be omitted through the model estimation process. We apply the Bayesian approach to GNSS displacement-rate data in Mainland China to obtain a unique GNSS velocity field for sub-regions with different tectonic backgrounds. We also examine the correlation between the estimated strain-rate fields and seismic activity by using the Molchan Error Diagram. The results show that the rate of maximum shear strain is positively correlated with the occurrence of earthquakes, indicating the strain rate can be used to augment probabilistic earthquake models for background seismicity forecasting.