<p>Gravity anomalies reflect the geophysical response to subsurface density structures. Traditionally, the terrain density is assumed to be a constant when calculating Bouguer gravity anomaly. But deviations from this assumption may induce high-frequency signals to Bouguer gravity anomaly. This study introduces a Bayesian method for computing Bouguer gravity anomaly. It incorporates a smoothness prior for the Bouguer gravity anomaly and estimates near-surface density parameters to minimize the Akaike’s Bayesian Information Criterion (ABIC) value. The effectiveness of this method is validated through theoretical model tests and calculations on two observed gravity profiles in Yunnan. The results indicate that the Bouguer gravity anomaly profiles estimated using the Bayesian approach need no extra filtering, exhibit correlations with the profile’s crustal structure, and effectively reveal subsurface crustal density variations. Moreover, the obtained density variations offer insights into the near-surface rock density of geological time. Specifically, Cenozoic formations have a density of roughly 2.65 to 2.90 g·cm⁻³, Mesozoic formations 2.61–2.91 g·cm⁻³, and Paleozoic formations 2.61–2.92 g·cm⁻³. Magmatic rock regions generally show higher density values. Additionally, these estimated densities show a positive correlation with the global <i>V</i>_<i>S</i>30 seismic velocity estimates, suggesting a new geophysical approach for seismic site classification. The findings of this study are significantly valuable for near-surface density estimation and Bouguer gravity anomaly calculations.</p>