Constraining the shape of a gravity anomalous body using reversible jump Markov chain Monte Carlo

Typical geophysical inversion problems are ill-posed, non-linear and non-unique. Sometimes the problem is trans-dimensional, where the number of unknown parameters is one of the unknowns, which makes the inverse problem even more challenging. Detecting the shape of a geophysical object underneath the earth surface from gravity anomaly is one of such complex problems, where the number of geometrical parameters is one of the unknowns. To deal with the difficulties of non-uniqueness, ill-conditioning and nonlinearity, a statistical Bayesian model inference approach is adopted. A reversible jump Markov chain Monte Carlo (RJMCMC) algorithm is proposed to overcome the difficulty of trans-dimensionality. Carefully designed within-model and between-model Markov chain moves are implemented to reduce the rate of generating inadmissible geometries, thus achieving good overall efficiency in the Monte Carlo sampler. Numerical experiments on a 2-D problem show that the proposed algorithm is capable of obtaining satisfactory solutions with quantifiable uncertainty to a challenging trans-dimensional geophysical inverse problem. Solutions from RJMCMC appear to be parsimonious for the given prior, in the sense that among the models satisfactorily represent the true model, models with higher posterior probabilities tend to have fewer number of parameters. The proposed numerical algorithm can be readily adapted to other similar trans-dimensional geophysical inverse applications. Keywords: trans-dimensional geophysical inversion, reversible jump Markov chain Monte Carlo; gravity anomaly.

💡 Research Summary

This paper tackles a notoriously difficult class of geophysical inverse problems: the estimation of the shape of a subsurface body from gravity anomaly data when the number of geometric parameters (i.e., the model dimension) is itself unknown. Such “trans‑dimensional” problems are ill‑posed, highly non‑linear, and non‑unique, making conventional deterministic inversion methods unreliable. To address these challenges, the authors adopt a fully Bayesian framework and develop a reversible‑jump Markov chain Monte Carlo (RJMCMC) algorithm that can explore both the parameter values and the model dimension within a single stochastic sampler.

The forward model assumes a two‑dimensional setting in which the anomalous body is represented by a polygon. The unknowns are the coordinates of the polygon’s vertices and the number of vertices. Prior distributions encode physical constraints (convexity, non‑self‑intersection) and a preference for parsimonious models, thereby penalising overly complex geometries. The likelihood is based on the misfit between observed gravity anomalies (contaminated with Gaussian noise) and the synthetic anomalies computed from the polygon via the standard Newtonian potential formulation. Because the posterior cannot be evaluated analytically, the authors resort to MCMC sampling.

RJMCMC is distinguished by two families of moves. Within‑model moves perturb the coordinates of existing vertices while respecting the geometric constraints; the proposal distribution is tuned to keep the acceptance rate high and to avoid generating inadmissible polygons. Between‑model moves change the dimensionality of the model: a “birth” move inserts a new vertex on an existing edge, and a “death” move removes a vertex. Both moves are carefully constructed to satisfy detailed balance: the Jacobian of the transformation is derived analytically, and the proposal probabilities are designed to reflect the prior on model complexity. By limiting the generation of illegal geometries, the algorithm achieves a high overall efficiency compared with naïve RJMCMC implementations.

The authors validate the method on synthetic data generated from a known five‑vertex polygon. After running the sampler for one million iterations, the posterior distribution concentrates on models that closely match the true shape while using the smallest number of vertices that still provide an adequate data fit. This demonstrates the algorithm’s ability to produce parsimonious solutions: among all models that explain the data, those with higher posterior probability tend to have fewer parameters, reflecting the Bayesian Occam’s razor effect. Moreover, the posterior samples provide full uncertainty quantification for both vertex locations and the number of vertices, something deterministic inversion cannot deliver.

Key contributions of the work include: (1) extending Bayesian trans‑dimensional inference to gravity‑based shape reconstruction; (2) designing geometry‑aware proposal mechanisms that dramatically reduce the rate of rejected, non‑physical models; (3) showing that the posterior naturally balances data fidelity against model complexity, yielding parsimonious yet accurate reconstructions; and (4) outlining how the same RJMCMC framework can be adapted to other geophysical modalities such as magnetics, electromagnetic, or seismic tomography.

In conclusion, the study provides a robust statistical tool for tackling trans‑dimensional geophysical inverse problems. By integrating model selection and parameter estimation into a single RJMCMC sampler, it delivers not only point estimates of subsurface geometry but also credible intervals that quantify epistemic uncertainty. The methodology is readily extensible to three‑dimensional settings and to real field data, opening avenues for more reliable and interpretable subsurface imaging in exploration geophysics.