Reversible-jump Markov chain Monte Carlo and Shamos–Hoey algorithms for two-dimensional gravity inversion

SUMMARY

In this study, a modified 2-D gravimetric inversion algorithm is presented that is based on the reversible-jump Markov chain Monte Carlo method with the Talwani equation. To ensure the validity of the Talwani equation and accurate gravity anomaly calculations, the Shamos–Hoey algorithm is incorporated as an additional acceptance condition to prevent intersections in the model polygon. This improves upon the method proposed by Luo by iteratively refining a polygon model based on gravity anomalies while maintaining physical validity. Additionally, we suggest a revision of the prior density function to better test the proposal models. This method estimates the shape and location of subsurface intrusions, providing valuable insights into subsurface geological structures. This positions the algorithm as a valuable tool for geophysical research.

1 INTRODUCTION

Geophysical inversion involves estimating subsurface properties from observed geophysical data, such as seismic, electromagnetic, or gravity measurements (Menke 2012). It has numerous applications in comparing observed data with theoretical models. However, the complexity and inherent ambiguity of geophysical inversion require advanced methodologies capable of accurately representing uncertainties in both data and models.

Within geophysical inversion, gravity inversion is notably valuable as a non-invasive technique to characterize subsurface structures through variations in the gravitational field, reflecting density changes (Talwani et al. 1959). However, gravity inversion faces substantial challenges due to inherent non-uniqueness and nonlinearity (Tarantola 2005; Menke 2012). Bayesian inversion, particularly when implemented with methods based on ensembles and integrated geological insights, effectively addresses these challenges by generating ensembles of plausible subsurface models (Kaipio & Somersalo 2005; Gelman et al. 2013).

Gravitational anomalies in geodesy and geophysics are deviations from the expected gravitational acceleration caused by variations in Earth’s density. In geophysics, these are measured using the CGS system (centimetre-gramme-second), the galileo (Gal) and the eötvös (E). Although the SI system is the international standard (Taylor 2001; Kovalevsky & Quinn 2004), these non-SI remain common because they offer a practical scale for the variations encountered in geophysical surveys (Torge & Müller 2012).

The gal unit (Gal) represents acceleration due to gravity.

The eötvös (E) unit represents the spatial rate of change (gradient) in gravitational acceleration:

According to CODATA-2018 (Mohr et al. 2016; Tiesinga et al. 2021), the universal gravitational constant is given as:

Gravitational anomalies are of considerable importance due to their direct link to density contrasts within the Earth, a key to interpreting subsurface structures and modelling geological contexts (Camacho et al. 2001; Montesinos et al. 2003). Analysing variations in gravitational acceleration and its gradient allows geophysicists to infer the presence, shape and density of features such as ore bodies, sedimentary basins and faults. Thus, gravitational anomalies act as proxies for subsurface density variations, facilitating detailed geological modelling (Menke 2012).

Bayesian statistics provides a robust framework for quantifying uncertainty in subsurface investigations (Kaipio & Somersalo 2005). Using probability distributions to represent data and model parameters is a core concept. By systematically updating prior beliefs with empirical evidence, posterior distributions are produced that refine understanding and quantify uncertainty (Green 1995; Zapata & Restrepo 2021).

This iterative process allows for progressive model enhancement, particularly valuable in complex applications such as subsurface exploration (McCalman et al. 2014). As more data are incorporated, the posterior distribution increasingly reflects new evidence, balancing prior knowledge and improving model accuracy (Metropolis et al. 1953; Hastings 1970).

Bayesian inversion provides a rigorous framework for addressing these complexities (Gelman et al. 2013). By combining prior geological knowledge with observational data, Bayesian methods effectively quantify uncertainties, producing probabilistic solutions rather than a single-deterministic case (Menke 2012). Such approaches have increasingly been adopted in geophysical studies, given their strengths in managing nonlinearities, high-dimensional parameter spaces and incomplete data (Luo 2010; Izquierdo et al. 2020). Tricarico (2013) introduced a planetary gravity inversion method that accommodates bodies of arbitrarily shaped bodies by combining matrix inversion with Monte Carlo sampling.

The integration of gravity inversion with Bayesian methods provides a robust framework for subsurface investigations, such as seafloor mapping, crustal density modelling and aquifer delineation (McCalman et al. 2014). This approach not only determines the geometry and density of anomalous bodies, but also incorporates expert knowledge. By quantifying uncertainty through posterior probability distributions (Luo 2010), it leads to more informed geological and tectonic analyses. Using geological priors and uncertainty quantification, Bayesian geophysical inversion creates probabilistic models that enhance interpretation beyond deterministic methods, especially beneficial in data-limited areas for reliable solutions (Menke 2012).

Monte Carlo techniques are particularly suited to Bayesian gravity inversion. These include Markov chain Monte Carlo (MCMC) algorithms such as Metropolis–Hastings (Metropolis et al. 1953; Hastings 1970), Gibbs sampling (Gelman et al. 2013), Hamiltonian Monte Carlo (Betancourt 2017), Hessian-informed MCMC (Liang et al. 2023), and reversible-jump MCMC (RJMCMC) (Green 1995; Luo 2010; Izquierdo et al. 2020).

Unlike deterministic methods, these algorithms explore numerous potential models, enabling robust quantification of uncertainty (Kruschke 2010). Furthermore, hierarchical transdimensional Bayesian models improve flexibility by determining geological complexity, such as the number of layers or faults, directly within the inference process (Green 1995; Sampietro & Capponi 2023).

Several other iterative methods have been developed beyond Monte Carlo techniques for gravity inversion. Krahenbuhl and Li (Krahenbuhl & Li 2006) introduced a binary inversion algorithm using a genetic algorithm and Tikhonov regularization. Chen et al. (2006) proposed a hybrid-encoding genetic algorithm (HEGA) that combines binary and decimal coding. In contrast, Chen et al. (2006) developed a deterministic iterative method based on the iterative expansion of anomalous bodies within a 3-D prismatic partition to fit gravity data.

Interdisciplinary Bayesian inversion has emerged as a key trend, combining gravity data with complementary geophysical observations such as magnetic, seismic and electromagnetic data sets (McCalman et al. 2014). Object-based approaches, including techniques, for example, RJMCMC (Green 1995), dynamically adjust model complexity and offer enhanced uncertainty quantification. These methods are particularly valuable in scenarios where subsurface structures are not well known a priori (Luo 2010; Zapata & Restrepo 2021).

Data fusion through the integration of multimodal data sets with methods such as Bayesian inference (Luo 2010; Zapata & Restrepo 2021) and refined genetic algorithms (Chen et al. 2006), is set to transform geosciences. It does so by efficiently managing complex inverse problems. The future will likely see machine learning (Zhou et al. 2024) for model generation, while the increasing volume of data calls for broader approaches to algorithm optimization for effective processing (Kaipio & Somersalo 2005). Importantly, using geological knowledge to integrate geophysical and geological data sets (McCalman et al. 2014) will constrain models, lower non-uniqueness and improve interpretations when combined with classical inversion or iterative optimization.

2 OUR METHODOLOGY

To simulate gravity inversion, we employ the method proposed by Luo (2010), which is based on Bayesian inference and the RJMCMC method plus some computational geometric additions. It is important to note that the number of statistical parameters is not constant during each iteration. We used a dextral reference system, with the positive vertical axis pointing upward and the positive horizontal axis pointing to the right, as proposed by Luo (2010).

A 2-D gravity anomaly is represented by a polygon with a variable number of vertices in a clockwise sequence, defined by their coordinates in the x–z plane. The gravity anomaly is calculated using the Talwani equation and the model is iteratively updated using the Metropolis et al. (1953) and Hastings (1970) algorithm.

To ensure physical validity and improve the accuracy of the inversion, we incorporate an additional acceptance condition into the Metropolis–Hastings algorithm (Metropolis et al. 1953). Specifically, we propose modifying the acceptance criteria to omit test models defined by polygons with self-intersecting edges (loops), as these configurations represent unrealistic geological structures.

To achieve this, we use the Shamos–Hoey algorithm (Shamos & Hoey 2008), an efficient method for detecting the existence of intersections in a set of segments. This approach ensures that only physically plausible models have been accepted, complementing the Metropolis–Hastings framework and contributing to more reliable inversion outcomes.

In this example (Fig. 1 ), we sampled horizontally equidistant points 31 from |$\tilde{x}_{0} = 0$| to |$\tilde{x}_{25} = 500$|⁠, for |$j=0,1,\dots ,30$|⁠. The vertical component of height |$\tilde{z}_{0}$| varies sinusoidally. The gravity anomalies for the target model were then calculated at these points and compared with those of the test models.

Fig. 1 illustrates the inversion process, highlighting the evolution of the test model polygon and its corresponding gravity anomaly as the algorithm progresses. The final model closely approximates the target model, which underlines the effectiveness of the inversion method. The figure also depicts the trajectory of each 10 centroids of the accepted iteration models. The centroids of the initial and middle test polygons were indicated with dots ‘|$\cdot$|’ and the centroid of the target model with a cross ‘|$\times$|’.

For an additional comparison, Fig. 2 tracks the paths of both the mean positions of the polygon vertices and the centroids. This visualization provides insight into how the shape and position of the model evolve throughout the inversion process.

2.1 The Talwani equation

The direct method for gravity modelling, based on the Talwani equation (Talwani et al. 1959), calculates the vertical component of gravity anomalies |$\Delta {g}$| caused by 2-D geological structures. These structures were represented by polygons with k vertices, where each vertex |$\boldsymbol {v}_{j}$| is defined by its coordinates |$(x_{j}, z_{j})$| in the |$xz$| plane. The equation uses these vertex coordinates to compute the gravity anomaly.

This work uses the Talwani equation (Talwani et al. 1959) for forward gravity modelling, coupled with inversion methods inspired by Chen et al. (2006) and Luo (2010) to estimate subsurface structures from gravity anomalies.

Luo (2010) and Chen et al. (2006) introduced methodologies that take advantage of this equation to solve complex inverse problems in potential field data, simulating gravity anomalies without requiring precise initial models. Chen’s approach employed a HEGA, combining binary and decimal coding for enhanced genetic operations, while Luo focused on a more classical stochastic optimization approach. More recently, Srigutomo et al. (2019) applied a Very Fast Simulated Annealing method for gravity inversion based also on the Blakely (1995) modification of the Talwani equation.

Building upon established Bayesian gravity inversion methodologies, particularly the framework proposed by Luo (2010), our study introduces enhancements to algorithmic efficiency and robustness. Specifically, we integrate computational geometric techniques to improve the exploration of the model space and achieve a more refined characterization of subsurface geological structures. These techniques, including those related to geometric intersections (Shamos & Hoey 2008), polygon processing (centroid) (Bourke 1997; Nürnberg 2013), and ellipse fitting (Kanatani 1994; Fitzgibbon et al. 1996; Halir & Flusser 1998), enable a more precise representation and manipulation of subsurface geometries within the Bayesian inversion process.

This polygon, denoted as |$\boldsymbol {\theta}_k$|⁠, has clockwise ordered vertices as |$\boldsymbol {\theta}_{k} = (\boldsymbol {v}_{0}, \boldsymbol {v}_{1}, \dots , \boldsymbol {v}_{k-1})$|⁠. A test model |$\mathcal {M}_{k}$| uses this polygon to represent a single intrusion in a 2-D vertical section. The polygon is characterised by |$2k$| parameters, corresponding to the horizontal and vertical positions of its vertices. The number of vertices, k, which can vary between 3 and 20, determines the complexity of the polygon and is stored in a pseudo-vector |$\boldsymbol {\theta}_{k}$| of |$2k$| elements. Importantly, k can change with each iteration, influencing the geometry of the modelled intrusion.

We use forward modelling with the Talwani equation (Talwani et al. 1959; Chen et al. 2006) to calculate gravity anomalies for a polygon-defined intrusion. Random changes were introduced iteratively into a test model, evaluated using the Metropolis–Hastings algorithm (Metropolis et al. 1953; Hastings 1970), to achieve better data fit. Bayesian inference and RJMCMC (Luo 2010; Izquierdo et al. 2020) enable variations in statistical parameters such as polygon vertex count for target model optimization.

Initially, an iterator i is defined to iterate over the vertices of the polygon, and an iterator j is defined for the points at which gravity anomalies will be calculated. Consequently, the coordinates |$(x_j,z_j)$| of the j-th vertex of the polygon were related to the coordinates |$(\tilde{x}_j,\tilde{z}_j)$| of the j-the anomaly measurement point.

The coefficients |$A_{i,j}$|⁠, |$B_{i,j}$| and |$C_{i,j}$| are given by the following expressions:

where |$\xi _{i,j}$| and |$\eta _{i,j}$| represent the relative coordinates of the polygon’s vertices with respect to each point at which the gravity anomaly is calculated:

Density contrast values are crucial for interpreting subsurface geological structures, as they represent differences in density between anomalous geological features and their host rocks.

For example, volcanic islands such as the Canary Islands (Spain) show contrasts around |$\pm 440\, {\rm kg}\, {\rm m}^{-3}$| for features such as magma chambers (low density) and intrusive bodies (high density) (Camacho et al. 2001).

The São Sebastião depression (Azores archipielago) has used contrasts from |$-500\,$| to |$350\, {\rm kg}\, {\rm m}^{-3}$| (Montesinos et al. 2003), while Monowai caldera (New Zealand) studies identified contrasts of about |$450\, {\rm kg}\, {\rm m}^{-3}$| for dense intrusions (Paulatto et al. 2014).

Similarly, studies in the Cerro Machín volcano (Colombia) report contrasts of |$-370\, {\rm kg}\, {\rm m}^{-3}$| between volcanic units and the host rock, and |$-400\, {\rm kg}\, {\rm m}^{-3}$| between different geological formations (Pedraza et al. 2022).

These cases highlight the importance of selecting density contrasts based on geological context to ensure accurate density modelling. In the current simulation models, the density contrast |$\Delta {\rho }$| is a pre-defined constant hyperparameter established as follows for demonstration purposes:

Although |$2G\Delta {\rho }$| is currently constant in these models, future implementations of the algorithm could incorporate a new step to propose changes in density contrast |$\Delta {\rho }$|⁠, which introduces it as an additional unknown parameter to be estimated along with the position of the shape vertex of the model, complementing the movements already established in Section 4.1.

The resulting gravity anomaly curves are shown in the upper panel of Fig. 1. For illustrative purposes and visualization, these curves were initially generated using 501 points, calculated using the direct method based on the Talwani eq. (1). These points represent a finely sampled representation of the anomaly curve.

However, it is important to note that the inversion algorithm does not use all these 501 points. Instead, the algorithm operates on a limited set of pre-defined gravity anomaly values 31, which serve as only input for the inversion process. These values are computed using the direct method for both the initial and the proposed models.

This setup reflects the kind of data that is typically available in a real-world educational setting. Alternatively, without referring to a specific target polygon, one could use a list of gravity anomalies, such as the 31 used in this example. However, generating the complete set of anomaly curves solely from these data would not be viable as a pedagogical demonstration.

3 COMPUTATIONAL GEOMETRY

We present three modifications to the Luo algorithm (Luo 2010), informed by the principles of computational geometry. First, we introduce an additional criterion for acceptance, in addition to the Metropolis–Hastings acceptance, to ascertain the geometric and geophysical validity of the proposed model. This is achieved by employing the Shamos & Hoey (2008) algorithm to efficiently detect and reject polygon models that exhibit self-intersections among line segments. Furthermore, we review the computation method for the centroid of a polygon (Bourke 1997; Nürnberg 2013), moving away from the arithmetic mean of the vertices. Finally, we implement new methodologies for fitting an ellipsoid within a modified prior density function to a specified set of points for a test model (Fitzgibbon et al. 1996).

3.1 The Shamos–Hoey algorithm

In addition to the method described by Talwani et al. (1959), Chen et al. (2006) and Luo (2010), we used the Shamos–Hoey algorithm as a rejection condition for self-intersection test polygons, where the Talwani equation (Talwani et al. 1959; Chen et al. 2006) is not valid. During the iterative process, it is possible for the polygons to develop loops. Each loop turns in the opposite direction, causing the contributions of their gravity anomalies to cancel each other out.

We have noticed that in the Talwani equation (Talwani et al. 1959; Chen et al. 2006), the chirality of the vertex order (clockwise or anticlockwise) affects the sign of the total value of |$\Delta {g}$|⁠, that is, if the vertices of the polygon are anticlockwise, |$\Delta {g}$| takes negative values.

If the polygon has at least one loop, such as intersections on its sides, for example an infinite-shaped (horizontal case) with a positive density difference |$(\Delta {\rho }>0)$|⁠, the clockwise side will have positive anomalies and the anticlockwise side will have negative anomalies. In other words, the anticlockwise half behaves as if it had a density contrast of equal magnitude but opposite in sign |$(-\Delta {\rho })$| compared to the clockwise half (Fig. 3).

Instead, if the loop model is eight-shape (vertical case), the |$\Delta {g}$| values of each part will partially cancel each other out. These results are undesirable, so it is necessary to find a way to avoid them. We suggest modifying the Luo gravity inversion algorithm (Luo 2010) by using the Shamos–Hoey algorithm (Shamos & Hoey 2008) as an additional condition to reject test models with self-intersection polygons.

During certain iterations of the simulation, loops may form that invalidate the Talwani equation (Talwani et al. 1959; Chen et al. 2006), making it impossible to accurately compare gravity anomalies. We used the GitHub library lycantropos/bentley_ottmann (Bentley & Ottmann 1979; Ibrakov 2020) to detect and address these loop formations.

In this example, there are five segments, labelled A, B, C, D and E. The algorithm analyses consecutive pairs of segments, progressing from top to bottom and from left to right, shown in numbered columns). For each segment, only the starting black dot is considered; the ending white dot is ignored (Fig. 4).

Initially, at position (0), the segment A is encountered. As the scan progresses to the right (1), the segment B appears below A. No intersection is detected between these two.

Continuing the scan (2), segment C appears below B; again, no intersection is found. At position (3), the segment D begins between B and C, resulting in the order |$A-B-D-C$|⁠. No intersections are detected between the consecutive pairs |$B-D$| or |$D-C$|⁠.

Similarly, no intersections are found in the new consecutive pairs at position (4), with order |$E-A-B-D-C$|⁠, or at position (5), |$E-A-B-C$|⁠. However, at position (6), when the total order becomes |$E-A-C$|⁠, an intersection is detected between the segments A and C.

The algorithm terminates upon detecting this intersection. The key observation is that the relative vertical order of the starting points of A and C has been reversed compared to their initial appearance at position (8). The positions (7), and (9) mark the end points of the remaining segments, but are not relevant to the detection of intersections.

3.2 Centroid of a polygon

Since we limit our models to simple polygons, it is essential to consider the centroid |$(C_{x}, C_{z})$| of a polygon with k vertices, which is defined as the point where all lines that divide the polygon into two parts of equal area intersect. This concept is crucial in computational geometry and can be used in various applications, such as designing algorithms for image processing or simulating physical phenomena. (Bourke 1997; Nürnberg 2013):

In this expression, |$(x_{i}, z_{i})$| represent the coordinates of the vertices of the polygon, and A is the area of the polygon

The equations only apply to polygons without intersections between their sides (Bourke 1997; Nürnberg 2013). Using the Shamos–Hoey algorithm (Shamos & Hoey 2008), researchers can systematically determine the presence of intersecting sides within a polygon. This capability enables the prompt rejection of a test model exhibiting such geometrical anomalies. Please refer to Section 3.1 for exploring this topic further.

3.3 Ellipsoidal fitting for a set of points

The procedure involves conducting an ellipsoidal fit by employing a least-squares approach using the polygon vertices of the test model, and defining the ellipse with the polynomial definition of a conic in vector form, thereby minimizing the distance between the model points and the ellipsoidal surface to find the ellipse parameters that best fit the model data as described by Fitzgibbon et al. (1996), improved by Halir & Flusser (1998), and discussed by Kanatani (1994).

where

or in another form:

In this case, the ellipses satisfy the following condition: |$4ac-b^{2}=1$|⁠, or in its matrix form (Fitzgibbon et al. 1996):,

where the matrix |$\mathbb {C}$| is given by:

Our input matrix takes the form:

The scatter matrix is defined as (Fitzgibbon et al. 1996):

Next, we proceed to solve the following system:

The conditions for a viable solution are |$\mathbb {A}=\mathbb {S}^{-1}\mathbb {C}$| and |$\lambda _{max}>0$|⁠, where |$\lambda _{max}$| denotes the largest eigenvalue of |$\mathbb {A}$|⁠.

We adapt the EllipseModel class from the Scikit-image library (Pedregosa et al. 2011; Christian 2021), for the elliptical fitting. This function takes as input a matrix of points, where each row represents a point in the plane. Subsequently, it calculates the covariance matrix of the points and performs an eigenvalue decomposition with NumPy (Harris et al. 2020) to find the ellipse that best fits the data, selecting the maximum positive eigenvalue. The function returns a tuple that contains the centre, the semiaxes of the ellipse, and the inclination angle of the ellipse.

In cases where a circular fit is chosen in the code or when performing an elliptical fit becomes unfeasible due to certain conditions, a circular fitting process is executed. This circular fit involves using the average of the vertex radii from the centroid of the polygon, assigning this average value to both the minor and major semiaxes, with an orientation angle set to zero.

4 REVERSIBLE-JUMP MARKOV CHAIN MONTE CARLO SIMULATION

In Bayesian statistics, both the observed data and the model parameters, defined as the vertices defining a polygon |$\boldsymbol {\theta}$| in our case, are treated as random variables (Luo 2010). This means we assign probability distributions to both, reflecting the inherent uncertainty in their values. This approach fundamentally differs from frequentist methods, where parameters are considered fixed, unknown values.

The goal of Bayesian inference is to update our prior beliefs about the parameters |$\boldsymbol {\theta}$| with the observed data |$\boldsymbol {y}$| to obtain the posterior distribution |$\pi({\boldsymbol {\theta}}|{\boldsymbol {y}})$|⁠.

Taking into account the following probabilities (McCalman et al. 2014; Zapata & Restrepo 2021):

• Prior density|$\pi({\boldsymbol {y}})$|⁠: The probability distribution of observing |$\boldsymbol {y}$| before considering the statistical parameters.

• Likelihood|$\pi({\boldsymbol {y}}|{\boldsymbol {\theta}})$|⁠: The probability distribution of observing |$\boldsymbol {y}$|⁠, given that the proposition |$\boldsymbol {\theta}$| is true.

• Marginal likelihood:|$\int \pi({\boldsymbol {y},\boldsymbol {\theta}}){\rm d}\boldsymbol {\theta}$|⁠, is also called the “evidence” and acts as a normalizing constant.

The joint distribution of both data |$\boldsymbol {y}$| and statistical parameters |$\boldsymbol {\theta}$| can be represented as Luo (2010):

The probability distribution of θ, given that the proposition y is true, or the posterior density, is determined through Bayes' theorem Luo 2010:

The RJMCMC algorithm advances over conventional MCMC algorithms by allowing “jumps” between parameter subspaces of different dimensionality (Green 1995). This makes RJMCMC suitable for a wider range of applications, including those where the number of parameters is not known in advance, like geophysics inversion. In our context, this parameter vector represents the coordinates (vertices) of the polygon in a plane that describes the model.

Unlike conventional MCMC, which is restricted to problems with a fixed dimensionality of the parameter vector, RJMCMC can handle cases where the number of parameters changes. RJMCMC is an extension of the conventional Metropolis–Hastings algorithm that allows for jumps between parameter subspaces of different dimensions.

This RJMCMC methodology is founded on the seminal work of Luo constraining the shape of a gravity anomalous body using the reversible-jump Markov chain Monte Carlo (Luo 2010). An additional rejection criterion for test models is incorporated, as delineated in Section 3.1. For a more detailed explanation of the RJMCMC method, please refer to the original articles by Luo and Green (Green 1995; Luo 2010).

The algorithm iteratively proposes random modifications to the test models, including simple movements, birth movements, or death movements, with probabilities modulated by a constraint on the maximum number of vertices (see Table 1). The acceptance of these modifications is determined by an acceptance criterion based on the likelihood function, prior density and birth–death probabilities, thereby ensuring convergence toward the gravity anomaly profile of the target model.

4.1 Movement types

This subsection details the movement types (simple, birth, or death) within the RJMCMC algorithm (Green 1995). The algorithm iteratively proposes changes to a polygon model, defined by a set of vertices (Luo 2010). The number of vertices can change, via birth- or death-movements, making this a trans-dimensional MCMC method.

At each iteration, a vertex is selected uniformly at random. Then, a movement type is chosen. Each has a different probability, depending on the number of existing vertices, as shown in Table 1. It is worth noting that triangle death movements are not allowed and there is a maximum limit of 20 vertices, prohibiting any birth movement beyond that number.

The probabilities of model movements: simple |$s_{k}$|⁠, birth |$b_{k}$| and death |$d_{k}$| during the simulation are normalized to sum to 1 for each value of k based on the number of vertices |$(k)$| in the proposal model. For intermediate model sizes |$(4-19)$|⁠, all types of movements have equal probabilities |$(1/3)$| to encourage exploration. At the minimum |$(k=3)$| and maximum |$(k=20)$| vertex counts, probabilities |$(1/2)$| are adjusted to prevent invalid moves, death of a triangle vertex or birth beyond the maximum vertex limit have 0 probability.

• Simple (Within-model) movement: A selected vertex |$\boldsymbol {v_i}=(x_i, z_i)$| is moved to a new position |$\boldsymbol {v_i}^\prime =(x_i^\prime , z_i^\prime )$| according to:

  $$\begin{eqnarray} \begin{aligned} x_{i}^{\prime } &= x_{i}+r\cos\, (\vartheta ), \\ z_{i}^{\prime } &= z_{i}+r\sin\, (\vartheta ), \end{aligned} \end{eqnarray}$$

  (18)

  In this scenario, a random variable |$\vartheta \sim U(0, 2\pi )$| is drawn from a uniform distribution over |$[0, 2\pi )$|⁠, representing the angle of movement.

  Another random variable |$r \sim N(0, \sigma ^{2}_{a})$| is drawn from a normal distribution with a mean of 0 and a standard deviation of |$\sigma _{a}$|⁠. The value of |$2\sigma _{a}$| is set to half the minimum distance is half the minimum distance to the adjacent vertices (⁠|$\boldsymbol {v_{i-1}}$| and |$\boldsymbol {v_{i+1}}$|⁠).

  The grey point cloud in Fig. 5 visually represents the distribution of possible new vertex positions resulting from this process. Each grey point corresponds to a potential |$(x_i^\prime , z_i^\prime )$| coordinate generated by a specific pair of |$(r,\vartheta )$| values. The black vector shows a single instance of such a displacement, moving the vertex from |$\boldsymbol {v_{i}}$| to |$\boldsymbol {v_i}^\prime$|⁠.

• Birth movement: A new vertex |$\boldsymbol {v_{i+1}^\prime }$| is inserted at the midpoint between the randomly selected vertex |$\boldsymbol {v_{i}}$| and its subsequent neighbour, |$\boldsymbol {v_{i+1}}$|⁠. This new vertex is then subjected to a simple movement, as defined by eq. (18) and illustrated in Fig. 6. For this birth move, |$2\sigma _a$| is set to one-quarter of the distance between |$\boldsymbol {v_{i}}$| and |$\boldsymbol {v_{i+1}}$|⁠. This scaling is consistent with the simple move, when a new point is created halfway between two existing vertices, the resulting triangle has sides equal to each other and equivalent to half the distance between those vertices. Therefore, |$2\sigma _{a}$| is equivalent to the case of a simple movement.

• Death movement: The selected vertex |$\boldsymbol {v_{i}}$| is removed from the polygon (Fig. 7). This is the inverse of the birth move; the value of |$2\sigma _{a}$| used in a hypothetical reverse birth (to recreate the deleted vertex) is set to be one-quarter of the distance between the two neighbour vertices |$\boldsymbol {v_{i-1}}$| and |$\boldsymbol {v_{i+1}}$|⁠. This ensures reversibility.

4.2 Likelihood function

We assume that the gravity anomaly values between the test and target models are normally distributed with a constant standard deviation |$\sigma$| throughout the simulation (Luo 2010). Consequently, the likelihood function introduced there enables us to determine the probability that the observed values align with the proposed model, taking into account the uncertainty tied to |$\sigma$|⁠.

Eq. (1) by Talwani (Talwani et al. 1959; Chen et al. 2006) describes the relationship between the gravitational anomalies of the target model |$\Delta {g}_{j}$| and the gravity anomalies of the test model |$\hat{y}({\boldsymbol {\theta}_{k},\tilde{x}_{j}})$|⁠.

It is important to note that |$\sigma$| should not be confused with |$\sigma _{a}$|⁠, which corresponds to the standard deviation of the random movements of the selected vertex and changes for each iteration, depending on the shape of the polygon. Various values for standard deviation were tested to find an appropriate value for simulation purposes. The value proposed by Luo of |$\sigma =0.2$| (Luo 2010) was found to be the most suitable.

4.3 Prior density

We define the Prior density (Luo 2010) as:

where m = [1, 1, 1], each element can be 0 or 1 as a mask vector; and the vector P = [P_0, P_1, P_2] is defined by:

where γ ≥ 1 is another hyperparameter constant throughout the simulation.

where |$\omega _{k}=(k-2)\pi /k, \quad {k}\ge {3}$|⁠, is the interior angle of a regular polygon with k vertices and |$\phi _{i}({\boldsymbol {\theta}_{k}})$| is the internal angle of vertex i.

In this equation, |$r^{\mathrm{ poly}}_{i}$| and |$r^{\mathrm{ ell}}_{i}$| represent the radii of the polygon and ellipse, respectively, calculated at the i-vertex, from the centre of the fitted ellipse or the centroid of the test polygon when an ellipse cannot be fitted. |$r_{m}$| is the average of the radii of the ellipse or circle, while |$\beta =0.25$| is a new hyperparameter.

With the mask vector |$\boldsymbol {m}$|⁠, it is possible to control the functions |$P_{0}$| (eq. 21), |$P_{1}$| (eq. 22) and |$P_{2}$| (eq. 23). The exponential function |$P_{0}$| penalizes a large number of vertices based on the value of the exponent |$\gamma$|⁠. In contrast, the other two functions, |$P_{1}$| and |$P_{2}$|⁠, are non-normalized Gaussian functions that penalize the deformity of the test polygon.

The function |$P_{1}$| measures the internal angles of the polygon and compares them with those of a regular polygon. For this case, we assume that |$\gamma =1.6$| is constant throughout the simulation (Luo 2010). Luo provides a more detailed explanation of |$P_0$| and |$P_1$| in his original article (Luo 2010).

The function |$P_{2}$| is a new proposition that compares the radii of the polygon to those corresponding to an ellipse (or circle) fitted by least squares using the algorithm proposed by Fitzgibbon (Fitzgibbon et al. 1996), as explained in Section 3.3.

In other words, the prior probability assumes the shape of the polygon, and an amorphous test polygon with many sides is considered less probable than a round one with fewer vertices.

4.4 Birth and death probabilities

The acceptance probability, |$\alpha$|⁠, for birth-and-death proposals within the RJMCMC framework is defined as the product of the likelihood ratio, the prior ratio, the proposal ratio, and the Jacobian, consistent with the methodology described by Green (Green 1995) and Luo (2010).

where the proposal density (Luo 2010) is defined as:

Here, the parameter vector associated with these auxiliary variables is vector of the auxiliary variables for the birth move is |$\boldsymbol {u}_{k}=({r,\vartheta })$| (Luo 2010). And the normalized Gaussian function with mean |$\mu$| and standard deviation |$\sigma _{a}$| evaluated at r is given by:

Finally, the Jacobian of the transformation from polar to Cartesian coordinates and r:

Essentially, when choosing to accept a birth move (adding new a vertex), the algorithm evaluates the likelihood ratio of the new versus current model, adjusts based on the models’ prior probabilities (prior ratio), considers proposal asymmetry via proposal ratio, and finally, multiplies by the Jacobian to account for variable space changes from the coordinate transformation.

Specifically, we define the probabilities of the proposed birth–death moves, denoted by |$\Upsilon _{k}(\sigma _{a})$|⁠, as the product of the prior ratio and the Jacobian term:

Here, the product |$2\pi r$| reflects two elements: |$2\pi$| comes from the uniform distribution of the angle |$\vartheta \sim U(0, 2\pi )$|⁠. The probabilities for simple moves |$(s_{k})$|⁠, birth moves |$(b_{k})$| and death moves |$(d_{k+1})$| are set to be equally probable, with the additional constraint that the number of vertices vertices k reaches a minimum of 3 and a maximum value set at 20, as shown in Table 1.

4.5 Acceptance criteria

In Monte Carlo simulations, acceptance criteria are essential to ensure that the generated samples reflect the target search space accurately. As originally established by Metropolis et al. (Metropolis et al. 1953), the algorithm evaluates a proposed move by comparing the current model with a new candidate model. This evaluation involves calculating an acceptance probability, |$\alpha$|⁠, which is determined in our context by factors such as how well each model explains the observed gravity anomalies and the likelihood associated with proposing the move. A random number |$u \sim U(0,1)$| from a uniform distribution over [0,1) (Metropolis et al. 1953) and accepting the move if:

The acceptance criteria govern whether a proposed new configuration is accepted or rejected, and they are crucial to thoroughly exploring the parameter space. These criteria depend on the type of proposed change, such as parameter adjustments within-model or dimension changes (birth–death) between models. Moves that might not immediately improve model performance can still be accepted on the basis of these criteria. This approach prevents the algorithm from getting stuck in local optima and promotes a more comprehensive search of the parameter space.

We use the acceptance probabilities derived from the work of Luo (2010) and Green (Green 1995). The acceptance probabilities in the algorithm are dependent on the type of movement and the number of parameters is not fixed.

• Simple movement: A single vertex is moved without changing the total number of vertices. The acceptance probability is determined by the ratio of the new and old likelihoods and priors.

  $$\begin{eqnarray} \alpha _{{i}\rightarrow {j}} = \frac{ \pi({\boldsymbol {y}}|{\boldsymbol {\theta}^{\prime }_{j}, \mathcal {M}_{j})} }{ \pi({\boldsymbol {y}}|{\boldsymbol {\theta}_{i}, \mathcal {M}_{i})} } \times \frac{ \pi({\boldsymbol {\theta}^{\prime }_{k}, \mathcal {M}_{k})} }{ \pi({\boldsymbol {\theta}_{k}, \mathcal {M}_{k})} }. \hphantom{\Upsilon _{k}} \end{eqnarray}$$

  (28)
• Birth movement: A new vertex is added, increasing the vertex count. The acceptance probability compares the updated model with the original, including the probability of the birth and death movements |$\Upsilon _k(\sigma _a)$| to create the vertex.

  $$\begin{eqnarray} \alpha _{{k}\rightarrow {k+1}} = \frac{ \pi({\boldsymbol {y}}|{\boldsymbol {\theta}_{k+1}, \mathcal {M}_{k+1})} }{ \pi({\boldsymbol {y}}|{\boldsymbol {\theta}_{k}, \mathcal {M}_{k})} } \times \frac{ \pi({\boldsymbol {\theta}_{k+1}, \mathcal {M}_{k+1}}) }{ \pi({\boldsymbol {\theta}_{k}, \mathcal {M}_{k}}) } \times \Upsilon _{k}(\sigma _{a}). \end{eqnarray}$$

  (29)
• Death movement: An existing vertex is removed, reducing the vertex count. The acceptance probability is computed as the inverse of the birth move, adjusting for the removal of the vertex.

  $$\begin{eqnarray} \alpha _{{k+1}\rightarrow {k}} = \frac{ \pi({\boldsymbol {y}}|{\boldsymbol {\theta}_{k}, \mathcal {M}_{k})} }{ \pi({\boldsymbol {y}}|{\boldsymbol {\theta}_{k+1}, \mathcal {M}_{k+1})} } \!\times \!\frac{ \pi({\boldsymbol {\theta}_{k}, \mathcal {M}_{k}}) }{ \pi({\boldsymbol {\theta}_{k+1}, \mathcal {M}_{k+1}}) }\! \times \left[{\Upsilon _{k}(\sigma _{a})}\right]^{-1}\!. \end{eqnarray}$$

  (30)

Therefore, after each successful iteration, the gravity anomaly values of the test model will become slightly closer to those of the target model and the test polygon model, which is the intrusion that needs to be identified. This process corresponds to a gravity inversion using iterative methods.

5 ALGORITHM

Our geophysical inversion approach employs reversible-Jump Markov-Chain Monte Carlo (Green 1995) to determine subsurface geometries from gravity anomaly data, using single-precision float32 arithmetic for improved computational efficiency despite a higher risk of numerical overflow instead of double-precision float64. The method incorporates a Metropolis–Hastings acceptance criterion and the Shamos–Hoey algorithm to ensure valid polygon configurations, while half-precision float16 was rejected due to its greater incidence of floating-point invalid value errors.

A set of 31 gravity anomalies is measured at the coordinates |$(\tilde{x}_{j},\tilde{z}_{j})$| for |$i=0,1,\dots ,50$|⁠, and a target or test polygon model represents the shape of the subsurface body (see Fig. 1). Additionally, the methodology uses modulo arithmetic to maintain a cyclic list of vertices; after reaching the last vertex, the next index wraps around to the first, ensuring that the vertex list remains sequentially closed, as implemented using the modulo operator (\%) in Python.

Since the model can involve the deletion or creation of vertices, the number of statistical parameters is not fixed. In this formulation, only the vertex coordinates are treated as parameters, their number remaining unknown, while the density contrast |$\Delta {\rho }$| is assumed to be fixed and known value.

The gravity anomaly for a given polygon is computed using the direct method defined by the Talwani equation (Talwani et al. 1959; Chen et al. 2006) (eq. 1), which is then used to compare the gravity anomalies of the target model with those calculated for the trial model in each iteration.

An initial model is chosen, specifically, a diamond square, and the number of vertices k, in each model is determined. In each iteration, random changes to the polygon vertices are proposed, generating a new trial model based on a type of movement: death, simple or birth, chosen at random according to Table 1.

The core of the algorithm is an iterative process that proposes modifications and accepts or rejects them based on the Metropolis–Hastings acceptance criterion (Metropolis et al. 1953; Hastings 1970) and the Shamos–Hoey criteria (Shamos & Hoey 2008), with iterations being accepted only if no floating-point overflow errors occur in NumPy (Harris et al. 2020) calculations.

To determine an acceptance probability |$\alpha$|⁠, the algorithm evaluates the effectiveness of the current test model in explaining observed gravity anomalies. Then a uniform random number |$u \sim U[0, 1)$| is generated using random.rand in Python. The move to the new model is accepted if |$u < \min (1, \alpha )$|⁠; otherwise, the current model is retained. This adheres to the Metropolis criterion outlined by Metropolis et al. (Metropolis et al. 1953).

Once all conditions are satisfied, the new list of vertex coordinates for the polygon is stored, and the process is repeated. If any condition is not met, the new model is discarded, the current model remains unchanged, and random modifications are applied to the vertices of the current model until a new trial model is accepted. With a sufficient number of accepted iterations, the trial models are expected to converge towards the target model.

6 RESULTS

There are two different ways to measure the evolution of test models during simulation: one is through the midpoint of the vertices of the polygons (Luo 2010) (greyscale gradient Fig. 2), and the other is by calculating the centroid of the polygons (Bourke 1997; Nürnberg 2013) (rainbow/turbo line in Fig. 2).

In Fig. 2, it is observed that the centroid, represented by the rainbow path, serves as a more effective parameter for assessing the progression of the simulation when contrasted with the midpoint of the vertices, depicted by the greyscale path. The older metric generates a considerably more erratic random walk, which complicates the analysis of model progression during the simulation. Furthermore, the employment of the centroid is evidenced to be a reliable indicator of the convergence of the test models towards the target synthetic model.

The algorithm performed a general number of 180 818 iterations, during which 18 000 test models were successfully accepted. As the iterative process of the algorithm progressed, the models that generated increasingly resembled the target model. This observation underscores the algorithm’s ability to progressively enhance and refine the solution over a series of iterations. However, it is important to note that the improvement trajectory was not invariably smooth and systematic.

The negative logarithms of likelihood and prior density exhibited dynamic behaviour, accompanied by variations in the count of the vertices of the test model at each iteration, as shown in Figs 8 and 9. These variations suggest efficient exploration of the solution space, eventually achieving a stable range in the gravity anomaly profiles. The logarithmic scale demonstrates process sensitivity, even at low likelihood values. The test model converges robustly within a range, rather than to a single point value. This range-based convergence is a direct consequence of the birth–death movements, which, while causing small changes in polygon shape and area, allow for comprehensive sampling of diverse model geometries.

To study the overall patterns of behaviour, we examined our simulation’s progress (Fig. 10) by grouping each 30 successful iterations. The progression of the accepted test models can essentially be categorized into two distinct phases: initially, the phase where the centroid of these test models quickly converges towards the target model, exhibiting only gradual alterations in shape. Later, in the second phase, when the models had almost converged to the target model, the centroid converged more slowly, while the test models experienced faster shape changes.

The computation time |$\Delta {t}$| required in groups of successful iterations 30 during each stage is relatively consistent, being shorter and more constant for the first stage (Fig. 10). A similar behaviour can be seen in the number of attempts per 30 accepted models, since they were correlated metrics. The highest number of overflow errors occurs mostly during the first stage, with small sporadic ‘spikes’ during the second stage (Fig. 10).

However, these ‘spikes’ can lead to inadequate models. If a polygon vertex deviates considerably from the rest of the polygon, the resulting gravity anomaly contributions are excluded from the likelihood function because they lie too far from the 31 gravity anomaly points considered. Consequently, test models with intrusions beyond the polygon area might be mistakenly accepted, an outcome to avoid using others methods.

As seen in the simulation results, simpler movements dominate the process, an outcome that is not explicitly programmed but emerges naturally from the algorithm’s stochastic framework (Figs 11 and 12). This phenomenon reinforces the reliability of the RJMCMC approach in converging toward realistic models, even in cases where the number of vertices in the target model remains unknown.

The initial model’s vertex count is not critical; these generally collapse during iterations if it is far from the target systematic model; hence, triangles and squares were chosen as initial models. The model then rapidly converges towards the target model’s centroid, as illustrated in Fig. 8 and videos in supplementary data. This demonstrates the ability to precisely determine the target model’s centroid location, independent of its precise shape, solely from gravity anomaly values.

The simulation results show that the vertex count of the generated models converges to a narrow range (Figs 9 and 12), which may differ from the target model’s vertex count. As a result, the number of vertices tends to converge in a small range of values due to changes caused by birth–death movements. Despite this potential difference, the observed convergence (within the expected bounds) demonstrates the ability of the estimation method to reliably determine the location of the centroid of the target model. This remains true even when the final test model and the target model have differing vertex numbers.

7 DISCUSSION

This research enhances 2-D gravity inversion by integrating Bayesian inference with computational geometry. This approach ensures that the inversion favours physically realistic polygon shapes and boosting subsurface model reliability. This versatile method has broad applicability and significant potential for subsurface characterization in fields such as mineral exploration, hydrogeology and geotechnical engineering.

In particular, the use of the Shamos–Hoey algorithm (Shamos & Hoey 2008) is central to our method. It automatically rejects self-intersecting polygons, thereby eliminating physically unrealistic models from the iterative inversion process.

In addition, we employ the centroid (Bourke 1997; Nürnberg 2013) of the polygon as a robust metric to track the evolution of the model. Unlike simple vertex averaging, the centroid offers a stable and reliable measure of the overall geometry, smoothing out the erratic behaviour of individual vertex positions and providing clear indicators of convergence.

Finally, to further improve the geophysical inversion process, we introduce an ellipsoidal fitting (Kanatani 1994; Fitzgibbon et al. 1996; Halir & Flusser 1998) technique to sharpen the a priori density function. By comparing test model shapes to optimal ellipses or circles, we enforce a shape-based constraint that penalises deviations from geologically realistic rounded forms.

While |$P_{0}$| penalize models with an excessive number of vertices through an exponential term with a fixed |$\gamma = 1.6$| (Luo 2010), |$P_{1}$| evaluates the internal angles against those of a regular polygon.

The new hyperparameter |$\beta$| in |$P_{2}$| (eq. 23) enhances the flexibility of an RJMCMC algorithm for gravity inversion by controlling the penalization of non-circular or non-elliptical shapes, complementing existing functions |$P_{0}$| (eq. 21) and |$P_{1}$| (eq. 22) to incorporate higher-order shape features for more realistic models.

This provides an enhanced mechanism for regulating penalties associated with deviations from optimal shapes. The adjustment facilitates a more realistic modelling approach. Through the analysis of simulation results, it becomes evident that these hyperparameters effectively minimise the occurrence of amorphous shapes in favour of simpler or rounded ones.

The simulation results validate our approach, demonstrating convergence to the target intrusion’s centroid. The algorithm efficiently locates subsurface features within a Bayesian framework, achieving accurate results without unnecessary computational complexity. This confirms the method’s ability to provide reliable gravity-inversion solutions.

This method is adequate in using density changes for subsurface analysis, useful in geosciences. It is most effective in geologically stable areas suitable for 2-D approximations, such as elongated faults, regions of steady volcanic activity, and consistent sedimentary basins.

Future work will enhance the algorithm’s applicability by introducing density contrast variations and extending the parameter space to include additional subsurface features. We also plan to incorporate further geological constraints and verify performance using real-world data sets. Together, these efforts should improve both accuracy and interpretability.

We propose developing multipolygon models parametrized by the centroids of regions rather than their vertices, utilizing shared vertex structures derived from geological layering or Voronoi tessellation (Voronoi 1908; Shamos & Hoey 1975), also known as Dirichlet tessellation (Dirichlet 1850) or Thiessen polygons (Thiessen 1911). These modelling strategies offer a structured approach to subsurface representation, ensuring robustness and smoother convergence. In particular, integrating these techniques with the joint estimation of density contrast and the alignment of uncertainty within the inversion framework is expected to enhance both the accuracy of results and the interpretability of the model.

This study reinforces the validity of the proposed methodology by demonstrating its ability to generate physically plausible and geologically meaningful subsurface models. The integration of computational geometry with Bayesian inversion ensures stable convergence, reduces non-uniqueness, and enhances model reliability, even in data-limited scenarios.

8 CONCLUSIONS

In this study, we introduce an extensive 2-D gravity inversion methodology that combines Bayesian inference with computational geometric techniques. This approach uses the Talwani equation for forward modelling purposes and employs the RJMCMC algorithm. Our framework incorporates multiple substantial innovations aimed at improving not only the physical validity of the model, but also its convergence.

A central feature is the use of the Shamos–Hoey algorithm, which automatically rejects self-intersecting polygon models and prevents erroneous gravity anomaly calculations. In addition, our approach replaces vertex averaging with centroid computations to provide a more robust metric for assessing model progression, resulting in smoother convergence trajectories and more stable updates.

Furthermore, we introduce a novel function that evaluates polygon shape deformities by comparing them with fitted ellipses or circles. This strategy refines prior density estimation beyond conventional methods that rely on angles and vertex counts. The algorithm demonstrates rapid convergence toward the target model’s centroid, with early iterations effectively capturing the essential central features based solely on gravity anomaly data.

This early performance underscores the robustness of the Bayesian framework and its capacity to extract considerable geological information even from limited data. Although the overall convergence behaviour is encouraging, it suggests a reliable approximation of complex subsurface geometries within a stochastic framework. The persistent polygon ‘spikes’ and other evolving aspects of the shape of the polygons highlight the need for further refinement of the algorithm.

Finally, the computational performance of the algorithm reflects its adaptability. Although computation times generally stabilize as the test models approach the target, occasional peaks due to the stochastic nature of the simulations, along with sporadic numerical anomalies, highlight the need for further optimization, especially for real-time applications or scenarios with high computational demands.

Looking ahead, future updates to the inversion methodology will focus on further optimizing computational performance and incorporating additional geological constraints. Such enhancements promise to extend the applicability of our approach to real-time and data-limited scenarios, ultimately advancing geophysical inversion techniques. In general, the presented method offers a route for achieving physically grounded solutions in geophysical exploration and subsurface modelling.

ACKNOWLEDGMENTS

The authors acknowledge the use of Python, along with the NumPy library (Harris et al. 2020), and the Google Colab environment (Edwards et al. 2024), for computations essential to this work. We thank the Universidad de Antioquia and the Faculty of Physics for providing a research environment and valuable technical expertise.

DATA AVAILABILITY

The data created or examined in this investigation can be found in the published manuscript and its accompanying supplementary information documents. Additional details are available upon reasonable request by contacting the corresponding author.

REFERENCES