Abstract. The spatiotemporal heterogeneity of b values has great potential for understanding the seismogenic process and assessing the seismic hazard. However, there is still much controversy about whether it exists or not, and an important reason is that the choice of subjective parameters has eroded the foundations of many researches. To overcome this problem, we used a recent developed non-parametric method based on the data-driven concept to calculate b values. The major steps of this method include: 1) perform a large number of Voronoi tessellation, Bayesian information criterion (BIC) value calculation and selection of the optimal models for the study area, and 2) use the ensemble median (Q2) and median absolute deviation (MAD) value to represent the final b value and its uncertainty. We investigated spatiotemporal variations of b values before and after the 2019 Changning MS 6.0 earthquake in Sichuan Basin, China. The results reveal a spatial volume with low pre-mainshock b values near the mainshock source region, and its size corresponds roughly with the rupture area of the mainshock. The anomalously high pre-mainshock b values distributed in the NE direction of the epicenter was interpreted to be related with fluid invasion or increased pore pressure. The decreases of b values during the aftershock sequence along with the occurrences of several strong aftershocks imply that b values could be an indicator of stress state. In addition, we found that although the distribution characteristics of b values obtained from different way of investigating are qualitatively consistent, they differ significantly in terms of their specific values, suggesting that the best way to study the heterogeneous pattern of b values is in the joint dimension of space-time rather than alone in time and space. Overall, our study emphasizes the importance of b value studies on assessing the earthquake hazards.