Inverse Gravimetric Problem Solving via Ellipsoidal Parameterization and Particle Swarm Optimization

1. Introduction

Inverse problems in gravimetry aim to infer the geometry and physical properties of subsurface structures from gravity measurements collected at the surface. These problems are inherently ill-posed in the sense of Hadamard [1], often exhibiting non-uniqueness, instability, and sensitivity to noise due to the smoothing nature of gravity data [2,3]. As multiple models can equally explain the observations, the inversion process is typically ambiguous and requires additional constraints to stabilize the solutions [4].

Classical inversion techniques address these challenges by incorporating regularization strategies, typically following the frameworks established by Tikhonov and further formalized by Menke [5] and Aster et al. [6]. However, the choice of regularization parameters and reference models can significantly influence the final solution, potentially introducing biases or masking important geological features [7,8]. Furthermore, the topography of the misfit functional in both linear and nonlinear inverse problems is often characterized by multiple minima and complex structures [4,8], making the inversion process even more delicate.

Probabilistic approaches, particularly those grounded in Bayesian inference [8], have been proposed to better capture the inherent uncertainties of inverse problems. Techniques such as Monte Carlo sampling [9,10] and transdimensional methods [8,11] allow for a more rigorous exploration of the solution space but at the cost of increased computational demands. Parsimonious approaches to Bayesian inversion have also been explored to balance complexity and tractability [9].

In gravity inversion, model reduction strategies have emerged as an effective alternative to full-space discretizations. Instead of parameterizing the subsurface with finely meshed grids, model reduction uses simple geometric bodies, such as prisms, spheres, cylinders, or ellipsoids, to describe anomalies with a small number of parameters [7,12]. These methods reduce the dimensionality of the problem and mitigate the curse of dimensionality effects [13], while often preserving key geological realism. Among model reduction strategies, ellipsoidal parameterization provides an advantageous balance between simplicity and flexibility. By representing subsurface cavities or density anomalies as ellipsoids, it is possible to capture variations in depth, orientation, and elongation, while maintaining analytical or semi-analytical expressions for forward modeling [12,14].

Global optimization algorithms, particularly Particle Swarm Optimization (PSO) [15,16], have proven highly effective for solving inverse problems characterized by complex, multimodal objective functions. PSO-based approaches allow not only the search for optimal solutions but also the sampling of plausible models, facilitating uncertainty quantification without relying on explicit regularization [17,18,19]. Applications of PSO to gravity inversion problems have been successfully demonstrated in various geological contexts [17,20,21,22]. Additionally, recent advances have explored the integration of PSO with Bayesian ideas for model space exploration [8], transdimensional approaches [11], and hybrid strategies to improve the robustness of inversion results in noisy environments [23].

In this study, we address the detection and characterization of anomalous bodies by representing them as prolate ellipsoids embedded in a homogeneous background. Assuming the density contrast is known, the model is defined by seven parameters: the coordinates of the ellipsoid center, the lengths of two semi-axes, the azimuth, and the tilt. We solve the inverse problem using PSO, with an emphasis on simultaneously achieving optimization and uncertainty quantification. We try also to introduce the density contrast as an additional parameter to estimate.

The proposed methodology is validated through synthetic experiments contaminated with Gaussian noise. The results demonstrate that ellipsoidal model reduction combined with swarm-based global optimization provides an efficient, robust, and interpretable framework for gravity inversion problems where the assessment of uncertainty is critical.

The rest of this article is structured as follows: Section 2 introduces the mathematical and computational framework. Sections 3 and 4 presents inversion results from synthetic noisy datasets. Section 5 discusses the implications and performance of the method. Section 6 offers conclusions and future research directions.

2. Materials and Methods

We consider the forward problem of computing the gravity anomaly produced by an ellipsoidal body with known density contrast $\Delta \rho$ embedded in a homogeneous background. The corresponding inverse problem consists of estimating the model parameters—position, size, orientation, and density contrast—from a finite set of surface gravity measurements. Due to the fact that gravity derives from a potential field, the limited number of observations, and the presence of measurement noise, the problem is ill-posed: its solution is generally non-unique and unstable [1,3]. The objective function associated with this inverse problem exhibits a highly irregular, non-convex landscape characterized by multiple local minima and flat regions. Even simple gravimetric configurations can produce highly complex misfit surfaces, as demonstrated by [4], which limits the effectiveness of local gradient-based optimization methods. To address these challenges, we adopt a model reduction strategy by parameterizing each anomaly as a prolate ellipsoid. The model vector $\mathbf{m}$ includes eight parameters per body: center coordinates $(x_c,y_c,z_c)$, semi-major and semi-minor axis lengths $(A,B)$, azimuth $\left(\alpha \right)$, tilt angle $\left(\varphi \right)$, and density contrast $(\Delta \rho )$. This results in a compact parameterization of dimension $8n_a$, where $n_a$ denotes the number of anomalous bodies. The gravitational response of each ellipsoid is computed analytically using the formulation of MacMillan [24], ensuring both accuracy and computational efficiency in the forward modeling.

Given the limitations of deterministic and full Bayesian approaches in high-dimensional nonlinear problems, we use Particle Swarm Optimization (PSO), a population-based global optimizer that does not require gradient information [15,16]. Its stochastic nature and swarm dynamics make it well-suited for rugged misfit landscapes and allow direct uncertainty estimation via particle dispersion. Compared to methods such as Markov Chain Monte Carlo (MCMC) technique [9,10], PSO offers a more tractable trade-off between sampling depth and computational cost. Techniques like Transitional MCMC or Simulated Tempering [11,25,26] share similar goals but are more complex to implement. We employ a resilient variant of PSO (CP-PSO), designed to balance convergence and diversity, enabling robust sampling of the posterior distribution without additional regularization [15].

Rather than relying on post hoc resampling or linearized error analysis, uncertainty is estimated directly from the final distribution of particle positions. This approach allows for the identification of anisotropic uncertainty structures in parameter space, consistent with physical resolution limits of gravity data [8,27]. Confidence intervals and plausible model ensembles are computed without assuming Gaussianity or linearity.

2.1. Gravity Inversion and Model Reduction

In a three-dimensional domain, the forward problem is given by

$$\Delta g\left(\mathbf{x}\right)=G{\int }_V\rho \left({\mathbf{x}}^{\prime }\right){\displaystyle \frac{(z-z^{\prime })}{|\mathbf{x}-{\mathbf{x}}^{\prime }{|}^3}}d{V}^{\prime },$$

(1)

where G is the gravitational constant, $\mathbf{x}=(x,y,z)$ the observation point, and $\rho \left({\mathbf{x}}^{\prime }\right)$ the density distribution in the volume V that produces the signal. This equation follows from Newton’s law of universal gravitation [23].

The inverse problem seeks $\rho \left({\mathbf{x}}^{\prime }\right)$ from measured $\Delta g\left(\mathbf{x}\right)$, typically expressed as:

$$\Delta g\left(\mathbf{x}\right)={\int }_{V}K(\mathbf{x},{\mathbf{x}}^{\prime })\rho \left({\mathbf{x}}^{\prime }\right)d{V}^{\prime },$$

(2)

with kernel:

$$K(\mathbf{x},{\mathbf{x}}^{\prime })=G{\displaystyle \frac{(z-z^{\prime })}{|\mathbf{x}-{\mathbf{x}}^{\prime }{|}^3}}.$$

(3)

Discretizing $\rho \left({\mathbf{x}}^{\prime }\right)$ using basis functions ${\varphi }_k$,

$$\rho \left({\mathbf{x}}^{\prime }\right)=\sum _{k=1}^N{\rho }_k{\varphi }_k\left({\mathbf{x}}^{\prime }\right),$$

(4)

leads to an underdetermined system $N\gg m$ where m is the number of observations. Regularization methods (e.g., Tikhonov, Bayesian priors) are often necessary [6,8,28].

2.2. Model Reduction in 2D Using Elliptical Bodies

An effective strategy to reduce the indeterminacy of inverse problems is dimensionality reduction through model reparameterization, in which anomalous bodies are represented as ellipses rather than discretizing the domain into a large number of cells. This approach substantially decreases the number of model parameters, although it introduces nonlinearity, as the forward model depends nonlinearly on the parameters that define the shape, position, and orientation of the anomalous body. In a homogeneous medium, a simple parameterization of an elliptical anomaly involves the set $[{\rho }_a,{\rho }_b,x_i,x_f,z_i,z_f]$, where ${\rho }_b$ is the background density, ${\rho }_a$ the anomaly density, and $x_i,x_f,z_i,z_f$ are the horizontal and vertical bounds of the elliptical region. Solving the inverse problem in this setting typically requires a 2D grid with $n_x$ and $n_z$ cells along the x- and z-axes, respectively. The reparameterized model reduces the dimensionality from $n_x{n}_z$ to $5n_a+1$, where five parameters describe each of the $n_a$ anomalies, and one parameter is used for the background density. The forward response is computed using either the Barbosa and Silva formulation for gridded solutions or Talwani’s method for polygonal bodies, both yielding comparable results; the latter is preferred due to its computational efficiency. This methodology was successfully applied by Muñiz et al. [19].

2.3. Model Reduction in 3D Using Ellipsoidal Bodies

In 3D environments, an alternative to classical discretization the anomalous bodies can be represented by prolate ellipsoids of revolution, significantly reducing the number of model parameters and, consequently, the dimensionality of the inverse problem. While this parameterization introduces nonlinearity, since the forward model becomes a nonlinear function of the geometric parameters, it retains sufficient flexibility to capture essential features of subsurface anomalies.

Each ellipsoid is described by a set of eight parameters: the coordinates of its center $(x_c,y_c,z_c)$, the lengths of the semi-major and semi-minor $(A,B)$ (assuming a prolate geometry, i.e. $A>B$), the azimuth and tilt angles $(\alpha ,\varphi )$ defining its orientation in space, and the density contrast $(\Delta \rho )$ relative to the surrounding medium. This compact representation enables a substantial reduction in the number of degrees of freedom, from $n_x{n}_z$ in cell-based models to $8n_a$, where $n_a$ is the number of ellipsoidal bodies.

2.3.1. Gravitational Attraction of a Homogeneous Prolate Ellipsoid

We consider a homogeneous ellipsoid with constant density $\rho$, and aim to compute the vertical component of the gravitational attraction at a point located along its axis of symmetry (the z-axis), i.e., at $(0,0,z)$, with $z>A$. The vertical gravitational attraction at such a point is given by:

$$g_z\left(z\right)=2\pi G\rho \left[{\displaystyle \frac{2A}{e^2}}\left({\displaystyle \frac{z}{\sqrt{z^2-A^2}}}arctan\left({\displaystyle \frac{\sqrt{z^2-A^2}}{eB}}\right)-1\right)\right],$$

where G is the gravitational constant and e is the eccentricity defined by

$$e=\sqrt{1-{\displaystyle \frac{B^2}{A^2}}}(\mathrm{with}A>B).$$

This exact analytical expression results from solving the gravitational potential integral for a prolate spheroid (an axially symmetric ellipsoid) along its symmetry axis, followed by differentiation to obtain the gravitational field [24]. It offers computational efficiency for forward modeling due to its closed-form nature.

The formula assumes that the observation point lies directly above the ellipsoid along its symmetry axis. If the ellipsoid is inclined or rotated (i.e., has a non-zero inclination or azimuth), the coordinates of the observation point must be transformed into the ellipsoid’s local reference frame. After applying the formula in that frame, the resulting gravitational vector must then be rotated back into the global coordinate system.

The gravitational attraction $g_z$ decreases with increasing distance z, and its dependence on the eccentricity e illustrates how the body’s shape influences the gravitational field: more elongated ellipsoids (higher e) concentrate mass closer to the center, thereby altering the field observed at external points.

The associated inverse problem involves estimating the set of ellipsoidal parameters that best fit the observed gravity data. By reducing the dimensionality of the parameter space and leveraging the geometric compactness of ellipsoidal shapes, this modeling strategy enhances computational tractability while preserving physical interpretability and structural coherence.

2.4. Noise and Uncertainty

Although it may seem counterintuitive, the uncertainty analysis is performed within an error region that is broader than the error associated with the best-fit model. This is because, in the presence of noise, the inverted model is affected by the noise in the data and becomes the only model that cannot exactly reproduce the synthetic observations. In contrast, the true model—if it were known—would actually have a larger misfit than the best-fit model and would lie within the equivalence region.

To illustrate this concept, consider a linear regression problem of the form $y^{*}\left(x\right)=a_0+a_{1}x$. Given a true parameter vector ${\mathbf{a}}_{\mathrm{true}}=(a_0,a_1)=(0.5,30)$, we generate $m\gg 2$ observed data points as:

$$y^{\mathrm{obs}}\left(x_i\right)=y^{*}\left(x_i\right)+{\delta }_i,i=1,\dots ,m,{\delta }_i\sim \mathcal{N}(0,\sigma ).$$

This is a purely overdetermined linear inverse problem, and the least-squares solution is given by:

$${\mathbf{a}}_{\mathrm{LS}}={\left(F^{T}F\right)}^{-1}{F}^T{\mathbf{y}}^{\mathrm{obs}},$$

where

$$F=\left[\begin{array}{cc}{\mathbf{1}}_{{\mathbb{R}}^m}& \mathbf{x}\end{array}\right]\in {\mathbb{R}}^{m\times 2},{\mathbf{1}}_{{\mathbb{R}}^m}={(1,\dots ,1)}^T,\mathbf{x}={(x_1,\dots ,x_m)}^T,{\mathbf{y}}^{\mathrm{obs}}={(y_1,\dots ,y_m)}^T.$$

The effect of noise in both linear and nonlinear inverse problems has been thoroughly studied by Fernández Martínez et al. [29,30]. One of the main conclusions of that analysis is that noise distorts the solution obtained not only by linear and nonlinear optimization methods but also by sampling techniques. In fact, noise can act as a form of implicit regularization, altering the topology of the uncertainty region.

Figure 1 illustrates the core motivation for uncertainty quantification in inverse problems. Even in this simple linear regression example, the presence of noise causes the least-squares solution to deviate from the true model, emphasizing the importance of exploring a family of plausible models instead of relying solely on a single estimate.

The figure presents the results of sampling the uncertainty region using a bagging algorithm. The least-squares solution (in green), which lies at the center of the sampling hypercube, clearly deviates from the true model (in black), which has a larger misfit than the least-squares solution. This justifies performing the uncertainty analysis within a relaxed error region:

$$M_{tol}={\displaystyle \frac{\parallel {\mathbf{y}}^{obs}-{\mathbf{y}}^{pre}{\parallel }_2}{\parallel {\mathbf{y}}^{obs}{\parallel }_2}}<tol,\mathrm{with}\mathrm{tol}>{\displaystyle \frac{\parallel F{\mathbf{a}}_{\mathrm{LS}}-{\mathbf{y}}^{\mathrm{obs}}{\parallel }_2}{\parallel {\mathbf{y}}^{\mathrm{obs}}{\parallel }_2}},$$

(5)

where ${\mathbf{y}}^{\mathrm{pre}}={(y^{*}\left(x_1\right),\dots ,y^{*}\left(x_m\right))}^T$ is the prediction from a sampled model.

In high-dimensional nonlinear problems, uncertainty analysis becomes even more challenging due to the curse of dimensionality [13] and the complex topology of the objective function, which often results in non-Gaussian posterior distributions with multiple low-error basins that may be disconnected in parameter space. For this reason, in this work, we employ members of the Particle Swarm Optimization (PSO) family as a sampling technique to explore the posterior distribution of the model parameters.

Uncertainty is quantified directly from the ensemble of particle positions during the PSO search. For each parameter $p_i$, a $95\%$ confidence interval is computed using:

$${\mathrm{CI}}_{95\%}\left(p_i\right)=[p_i^{(min\%)},p_i^{(max\%)}],$$

(6)

where $p_i^{\left(q\right)}$ is the q-th percentile of sampled values. This approach does not assume Gaussianity or linearity, and maps of parameter uncertainty provide insight into model resolution.

2.5. Model Inversion Via the PSO Family

Particle Swarm Optimization (PSO) was introduced into geophysical inversion more than a decade ago [12,31,32], and has since demonstrated strong performance in gravity inversion problems involving complex geological structures, such as faults with variable density contrasts [33,34,35]. In particular, PSO has been widely applied in gravity inversion studies [21,22,36,37,38], showing its ability to deal with nonlinearity, discontinuities, and non-uniqueness in the inverse problem.

Although Monte Carlo-based algorithms can also be applied to such problems [39], PSO offers specific advantages for posterior uncertainty analysis, especially its capacity to explore equivalence regions more effectively. Nevertheless, both PSO and other global optimization algorithms share common limitations: they require a large number of forward model evaluations, which becomes computationally expensive when dealing with high-dimensional parameter spaces or forward models with high computational cost.

The PSO family member RR-PSO was used due to its high exploratory and exploitative capacity. PSO consists of various variants, each with distinct behaviors, some exhibit stronger exploratory characteristics, while others are more exploitative. Exploratory behavior means that the search space is sampled more uniformly across iterations, whereas exploitative behavior causes the particle swarm to concentrate in regions with better misfit values as iterations progress. The optimal PSO configuration balances exploration and exploitation to ensure both an accurate estimation of the cost function’s optimum and a well-distributed sampling of the search space. This prevents the algorithm from becoming trapped in local minima and allows for a robust uncertainty analysis.

Particle Swarm Optimization (PSO) offers several advantages in gravimetric inversion, including its simplicity of implementation, low computational cost, and robustness in exploring complex, high-dimensional search spaces ([15]). Unlike gradient-based methods, PSO does not require the computation of derivatives, making it suitable for non-linear and ill-posed inverse problems. Its population-based nature enhances the global search capability, reducing the risk of being trapped in local minima. When configured in its exploratory form, PSO also enables a meaningful uncertainty analysis, despite not being a Bayesian method, providing valuable insights into the stability and resolution of the inferred models. These features have led to its successful application in geophysical inversion problems, particularly in gravity data interpretation [22,36].

Figure 2 provides a schematic overview of the PSO-based inversion workflow. It outlines the key steps involved in the estimation of ellipsoidal model parameters and the quantification of uncertainty, from particle initialization and misfit evaluation to the final extraction of optimal and plausible solutions. This flowchart emphasizes the dual role of PSO in both optimization and probabilistic exploration of the parameter space.

3. Results

This section presents the inversion results of different synthetic data sets generated from a model composed of one, two, or three ellipsoids and contaminated with Gaussian noise. The true model is compared with the inverted model and the average model, along with its dispersion, obtained by sampling the nonlinear uncertainty region. The trade-offs between the model parameters are also analyzed through their sampled posterior distribution.

3.1. Case of one ellipsoid with $10\%$ Gaussian noise

The first case involves a model composed of a single ellipsoid (8 parameters), with data contaminated by Gaussian noise with zero mean and a standard deviation equal to $10\%$ of the average value of the generated anomaly

$$\delta \mathbf{g}⇝N(0,\sigma ),\sigma =0.1·mean\left({\mathbf{g}}^{obs}\right)$$

Table 1 summarizes the performance of the inversion method by comparing the true model (TM), the best-fitting model (BM) with the lowest prediction error ($7.28\%$), and the median of models with less than $10\%$ prediction error (MedM), along with uncertainty metrics. The BM closely reproduces the TM parameters, indicating high inversion accuracy, particularly for center position ($x_c$, $y_c$, $z_c$), shape (A, B), and density ($\rho$). In contrast, the MedM shows larger deviations, especially in orientation angles ($\alpha$, $\theta$), highlighting variability within the uncertainty region. This is further confirmed by the interquartile range (IQR) and its normalized version (IQRn)

$$IQRn=\left|{\displaystyle \frac{IQR}{MedM}}\right|,$$

which show that orientation parameters are less constrained ($IQRn>0.8$), while density and center coordinates are better resolved ($IQRn<0.6$). The parameter bounds (Min–Max) encompass the true values, ensuring an adequate search space. Overall, the method yields robust estimates for most parameters, though orientation remains more sensitive and uncertain.

By performing a sensitivity analysis of the model parameters over 500 independent simulations, using the same true model (TM) and identical search limits, the results summarized in Table 2 were obtained:

The coefficient of variation using standard deviation is given by:

$$cvs=\left|{\displaystyle \frac{\sigma }{\mu }}\right|,$$

where $\sigma$ is the standard deviation and $\mu$ is the mean of the accepted models (those with a prediction error below $10\%$). This coefficient quantifies the relative variability of each model parameter: lower values indicate more consistent and reliable estimates, whereas higher values reflect greater uncertainty or dispersion in the parameter space.

From Table 2, it is evident that the parameters $x_c$, $y_c$, $z_c$, and $\rho$ show lower variability (with IQRn and cvs values well below 0.6), suggesting stable estimations across the accepted models. In contrast, the geometrical parameters A, B, $\alpha$, and $\theta$ exhibit greater variability, particularly evident in their higher IQRn and cvs values, close to or exceeding $0.7$, indicating greater sensitivity or ambiguity in their estimation under the given inversion scheme. In this case, the lowest prediction error is $7.8\%$.

The figures presented in this subsection further illustrate both the accuracy and robustness of the inversion results. In the synthetic experiments, the recovered ellipsoidal models closely match the true anomalies in location, geometry, and density contrast, even in the presence of synthetic noise.

Uncertainty is explored through both visual inspection of model ensembles and statistical analysis of parameter distributions. The anisotropic nature of the uncertainty is consistent with the known resolution limitations of gravity data, where depth is better constrained than lateral extent or orientation. Figure 3 shows the marginal distributions and pairwise correlations of the model parameters obtained from the inversion of the one-ellipsoid synthetic case. The distributions are generally non-Gaussian, and the parameters appear uncorrelated, supporting the independence assumption in the PSO-based sampling.

Figure 4 displays the true ellipsoidal anomaly and the best-fitting model recovered by the inversion and their spatial positions using plan and elevation views. The strong spatial agreement between the two ellipsoids (true and inverted) confirms that the method accurately reconstructs the anomaly’s geometry and position even in the presence of noise.

Figure 5 illustrates the convergence behavior of the PSO algorithm during the inversion of the single-ellipsoid model. The misfit error stabilizes rapidly (around 100 iterations), reaching a minimum value of $7.28\%$, indicating that the algorithm efficiently locates a viable solution in the parameter space.

Finally, Figure 6 compares the observed gravity anomaly with the prediction from the best-fitting model in the 1-body case. The close match between the two confirms the fidelity of the forward model and the effectiveness of the inversion.

These results also highlight the importance of sensitivity analysis when evaluating the robustness of geophysical inversion models. Parameters with low IQRn and cvs values (such as the centroid coordinates and density contrast) can be considered well-constrained by the data, whereas parameters with higher variability (e.g., shape and orientation angles) may require additional regularization, prior information, or refined data acquisition strategies to improve their estimation. This analysis supports more informed model selection and uncertainty quantification, reinforcing the reliability of the inversion outcomes in practical scenarios.

3.2. Case of two ellipsoids with $10\%$ Gaussian noise

This second case involves a model composed of two ellipsoids (16 parameters), with data contaminated by Gaussian noise with zero mean and a standard deviation equal to $10\%$ of the average value of the generated anomaly, as in the one ellipsoid case.

Table 3 presents a comparison between the true model that generated the data and the inversion results: the model with the lowest prediction error ($4.62\%$), the median model obtained from the $10\%$ uncertainty region, the interquartile range, the model parameter coefficient of variation, and the maximum and minimum parameter values of the models sampled within that region.

Table 3 and Table 4 present the inversion results for two ellipsoidal bodies, including the true model (TM), the best-fitting model (BM), and the median model (MedM) among those with prediction errors below $10\%$. These results are complemented with the interquartile range (IQR), the normalized coefficient of variation (IQRn), and the minimum and maximum parameter values within the $10\%$ uncertainty region, providing a comprehensive picture of parameter recovery and associated uncertainty.

For both ellipsoids, the BM closely approximates the TM, with small deviations in most parameters. In Ellipsoid 1, parameters such as $x_c$, $z_c$, $\theta$, and $\rho$ exhibit errors around $6\%$, while $y_c$, A, and B show errors between $7.3\%$ and $8.6\%$. The azimuth $\alpha$ deviates by approximately $10\%$. For Ellipsoid 2, the location coordinates, azimuth, and density are recovered with errors below $4.1\%$, although errors increase for A ($8\%$), B ($12\%$), and $\theta$ ($21\%$).

The MedM shows larger deviations, particularly in parameters related to geometry and orientation, such as the semi-axes and angular variables, reflecting the non-uniqueness of the inverse solution. This is consistent with the IQR values, which reveal significant dispersion in these parameters, especially for $\alpha$, $\theta$, A, and B.

The IQRn metric underscores the varying degree of parameter identifiability. The dip angle exhibits the highest IQRn values in both ellipsoids, indicating its strong sensitivity to noise and nonlinearity. In contrast, the density contrast $\rho$ consistently presents the lowest IQRn, suggesting greater stability and data sensitivity. Parameters can be grouped by their IQRn behavior: (i) low variability: $\rho$; (ii) moderate variability: $x_c$ and $y_c$; and (iii) high variability: $z_c$, A, B, $\alpha$, and $\theta$.

The broad ranges between minimum and maximum parameter values further highlight the inherent ambiguity of the inversion process, particularly in geometrical and directional parameters. These findings emphasize the necessity of uncertainty quantification in inversion problems, as even well-fitting models may be accompanied by significant variability in parameter estimates, especially under complex or noisy conditions.

The spatial distribution of the true and recovered ellipsoids demonstrates that the inversion method is capable of capturing the essential features of both anomalies simultaneously. As shown in the plan and front views (Figure 7), the geometry and location of the best-fitting models resemble those of the true bodies, even under significant noise contamination. This highlights the method’s ability to handle increased model complexity without sacrificing resolution.

The convergence behavior of the PSO algorithm in this two-body case further supports its robustness, as it is shown in Figure 8. Although more complex than the single-anomaly scenario, the swarm efficiently reaches a stable region in the parameter space, as indicated by the stepwise reduction in the misfit error. The stabilization of the error curve reflects the successful balance between exploration and exploitation during the optimization process.

The match between the observed gravity anomaly and the prediction from the best-fitting model confirms the validity of the inversion results (see Figure 9). The predicted response effectively reproduces the main features of the observed data, indicating that the model not only fits globally but also honors local anomaly structures.

3.3. Exploratory capability of the PSO algorithm

The performance of Particle Swarm Optimization (PSO) strongly depends on its ability to effectively explore the search space. This exploratory behavior is influenced by both the swarm size and the dynamic evolution of the particle distribution throughout the optimization process. In this subsection, we analyze several aspects that provide insight into how exploration is achieved and maintained.

Figure 10(a) and 10(b) illustrate the internal behavior of the swarm over the course of the iterations. Figure 10(a) shows the percentage of swarm dispersion relative to its center of mass as a function of the number of PSO iterations. A brief initial exploratory phase is observed, during which the dispersion drops from approximately $24\%$ to $4\%$ within the first 10 iterations. Around iteration 30, the dispersion begins to increase again, reaching $8\%$ at approximately 100 iterations, and remains stable at that level through the full 200 iterations. This trend reflects a well-balanced exploration-to-exploitation transition, where initial diversity gives way to convergence, followed by a moderate recovery of diversity to avoid premature stagnation.

Figure 10(b) shows the percentage of particles that remain within the parameter search space over time. This percentage decreases sharply during the first three iterations, likely due to the random initialization of the swarm. However, by iteration 8, over $65\%$ of particles are already within valid bounds, and by iteration 15 this rises to over $75\%$, remaining steady thereafter. These results confirm that the PSO configuration effectively steers particles toward feasible regions early in the optimization.

In addition to the temporal evolution of swarm dynamics, the effect of swarm size on performance is analyzed. Figure 11(a) shows the dependence of CPU time on swarm size. This relationship is observed to be linear and increasing across most of the domain. Figure 11(b) displays the relative error of the best model as a function of swarm size. It can be seen that with a swarm of around 30 particles, the low-error equivalence region (less than $9\%$) is already reached, and with approximately 200 particles, the minimum error (below $8\%$) is achieved. This result suggests that it is more efficient to reduce the number of particles (between 100 and 200) and instead increase the number of iterations to more quickly obtain low-error models. However, given the stochastic nature of the algorithm, it is possible that in a specific run (as shown in the figure), the relative error may increase with swarm size.

Together, these results emphasize the importance of tuning both the swarm size and its dynamic behavior throughout the optimization process. Early exploration followed by moderate convergence and stable particle retention within the search region are key elements in ensuring both robust inversion and meaningful uncertainty quantification.

4. Inverse Gravimetric Detection of Water-Filled Cavities Using Ellipsoidal Model Reduction

The detection of subsurface cavities filled with water constitutes a significant inverse problem in applied geophysics, particularly in environmental, hydrogeological, and engineering applications. These cavities typically generate subtle negative gravity anomalies due to the lower density of water relative to the surrounding rock matrix. Accurately resolving such structures requires robust inversion methods capable of handling noisy, sparse, and non-unique data.

In this study, the cavity is modeled as a prolate ellipsoid, offering both analytical simplicity and geometric flexibility to represent realistic voids. The model is defined by a reduced set of parameters: the coordinates of the ellipsoid center, the lengths of the major and minor semi-axes, and its spatial orientation (azimuth and inclination). The density contrast with respect to the surrounding karstic host rock is considered known, as the densities of both water and typical karst limestones are well constrained. This assumption allows us to reduce the number of free parameters from eight to seven, transforming the ill-posed inversion problem into a lower-dimensional nonlinear optimization problem.

To solve the inverse problem, we employ Particle Swarm Optimization (PSO), a global, population-based optimization algorithm particularly suitable for nonlinear and multimodal objective functions. The misfit function is defined as the sum of squared differences between observed gravity values and those predicted by the forward model. Regularization is not explicitly applied, as the dimensionality reduction introduced by ellipsoidal modeling provides an implicit form of model stabilization.

The use of ellipsoidal geometries not only enables an interpretable and compact description of the anomaly but also facilitates the quantification of uncertainty through posterior sampling. Key attributes such as depth, volume, and orientation of the cavity can be reliably estimated, parameters that are critical for groundwater management, tunnel design, and karst risk assessment. Moreover, this modeling framework supports risk evaluation and improves communication with stakeholders in environmental and engineering contexts.

Water-filled karst cavities may vary substantially in shape and depth depending on the geological setting and the degree of karstification. Based on empirical data from geophysical surveys and cave studies, typical dimensions range from 5 to 50 meters in length (semi-axis A), 2 to 20 meters in width (semi-axis B), and 10 to 80 meters in depth to the cavity center ($z_c$), depending on lithological and hydrogeological conditions. A narrow water-filled cavity typically ranges from 10 to 30 meters in length, 2 to 6 meters in width, and lies at depths between 10 and 60 meters, with estimated volumes commonly between 100 and 800 cubic meters depending on geological conditions. In addition, the ratio between the two semi-axes ranges from $2.5$ to 5.

In our synthetic experiment, we simulate gravity data from a water-filled cavity embedded in a homogeneous karstic medium. The anomaly is modeled as a prolate ellipsoid with known density contrast, leading to a seven-parameter inversion problem. The results of the inversion are summarized in Table 5.

The prediction error is significantly smaller than in the general case of ellipsoids with eight free parameters, as the density contrast is assumed to be known. This simplification reduces the number of estimated parameters from eight to seven, thereby decreasing the dimensionality of the search space and improving inversion stability. Furthermore, when the cavity is thinner than average, it is necessary to constrain the parameter search bounds more tightly around the true values. This ensures that the algorithm can recover an ellipsoid with geometry and location close to the actual one, as illustrated in Table 6.

In the case of a standard-size cavity, the spatial agreement between the true and inverted ellipsoids, as shown in Figure 12, confirms the effectiveness of the reduced-parameter inversion strategy. The recovered geometry closely matches the actual one in both plan and elevation views, despite the presence of noise. These results suggest that assuming a known density contrast enhances the stability and reliability of the inversion, making it especially suitable for practical applications.

The convergence behavior of the algorithm, illustrated in Figure 13, reflects the increased stability introduced by the reduced parameter space. The PSO search rapidly identifies a region of low prediction error, benefiting from the lower dimensionality of the problem. This efficient convergence is particularly advantageous in environmental and engineering contexts, where reliable model estimates are required within limited computational time.

Table 5 and Table 6 summarize the inversion results for two synthetic test cases involving water-filled karstic cavities of different sizes. The statistical metrics reported include the best-fitting model (BM), the median of models with prediction error below $10\%$ (MedM), the mean model (MeanM), the interquartile range (IQR), and the corresponding normalized metrics (IQRn and coefficient of variation, cvs). The results show that the inversion is able to recover the true geometry and position of the cavity with reasonable accuracy and low uncertainty. In particular, the relative error is significantly lower than in full ellipsoidal inversions involving eight parameters, owing to the fact that the density contrast is assumed known in this scenario. This effectively reduces the dimensionality of the parameter space from eight to seven, simplifying the optimization landscape and enhancing the stability of the inversion. Furthermore, as illustrated in Table 6, for small-sized cavities it becomes necessary to define tighter search bounds to ensure convergence towards the true ellipsoidal configuration. These results highlight the advantages of incorporating prior geological knowledge to constrain the inversion and reduce non-uniqueness.

Figure 14 compares the gravity anomaly derived from noisy synthetic observations with the response predicted by the best-fitting model. The excellent agreement between the two graphics confirms the accuracy of the forward modeling and demonstrates the ability of the inversion scheme to reproduce realistic gravitational signatures of water-filled cavities.

In the narrow cavity case, accurate inversion results are achieved only when the parameter search bounds are appropriately tightened. These constraints are essential for guiding the algorithm toward the true solution, as the limited spatial extent of the anomaly makes it more difficult to detect under broad or uninformed priors. As shown in Figure 15, the recovered ellipsoid closely matches the true geometry, confirming the importance of incorporating prior knowledge when targeting small, high-contrast structures.

As shown in Figure 16, the convergence curve indicates that the relative error rapidly reaches a stable minimum, demonstrating the method’s reliability in detecting narrow anomalies, despite the reduced sensitivity of gravity data to such small-scale features.

Figure 17 shows that the predicted gravity response closely reproduces the observed anomaly, even in the presence of noise. This agreement confirms that the inversion method is capable of accurately recovering the gravitational signature of narrow, low-density cavities, provided that appropriate constraints are applied during the optimization process.

5. Discussion

The results presented in this study confirm the viability of the proposed model reduction strategy for gravity inversion, based on ellipsoidal parameterization and global optimization via Particle Swarm Optimization (PSO). By reducing the number of free parameters and leveraging the geometric simplicity of ellipsoids, the method offers a computationally efficient alternative to traditional discretized inversion schemes.

The uncertainty analysis revealed anisotropic confidence intervals, particularly in horizontal location and orientation, consistent with the known limitations of gravimetric sensitivity. In contrast, the estimation of the ellipsoid’s center and the density contrast proved to be more stable, benefiting from the higher resolution of gravity measurements with respect to depth variations near the surface. These results highlight the importance of uncertainty quantification as an integral part of the inversion process.

The application to noisy data further demonstrated the practical utility of the method. Notably, the inversion was performed without the need for explicit regularization, relying instead on the implicit regularizing effect of geometric model reduction and global sampling.

Compared to classical inversion techniques, the proposed approach offers several advantages. First, it eliminates the need for arbitrary discretization and mesh design, which can introduce bias or increase computational costs. Second, the PSO framework enables both optimization and exploration of the solution space, generating a distribution of plausible models rather than a single best-fit solution. This capability is particularly valuable in risk-sensitive applications such as cavity detection. In this paper, we have demonstrated the application to the detection of a water-filled karst cavity using gravimetric methods, considering both the standard case and that of a narrow cavity.

Nonetheless, the method has certain limitations. The use of simple ellipsoidal geometries may lead to an oversimplification of complex subsurface structures, potentially overlooking small-scale heterogeneities. Furthermore, the current implementation assumes a homogeneous background, which may not hold in geologically complex settings. These assumptions should be carefully evaluated when applying the method to new environments.

Future research could address these limitations by incorporating multi-ellipsoid models or extending the approach to heterogeneous backgrounds. Additionally, integrating gravimetric data with complementary geophysical measurements (e.g., seismic or electrical) through joint inversion could further constrain the solution space and enhance interpretability. Finally, embedding this approach within a Bayesian framework could provide a more rigorous statistical treatment of uncertainty while maintaining computational efficiency through PSO-based posterior sampling.

Overall, the results support the use of ellipsoidal model reduction and swarm-based global optimization as a robust and interpretable strategy for gravity inversion problems, particularly in scenarios where model simplicity, efficiency, and uncertainty assessment are critical.

6. Conclusions

This study introduces a reduced-parameter approach to three-dimensional gravity inversion that models subsurface anomalies as ellipsoidal bodies and employs Particle Swarm Optimization (PSO) for parameter estimation. The method significantly reduces the complexity of the inverse problem while providing a practical means of exploring the associated uncertainty. By avoiding traditional regularization and instead relying on the inherent diversity of the particle swarm, the approach yields stable and interpretable models without sacrificing accuracy.

The 3D inversion results from synthetic experiments demonstrate that the method is capable of accurately recovering the position and density contrast of ellipsoidal anomalies, even when the data are contaminated with $10\%$ Gaussian noise. Although parameters related to the shape and orientation of the anomalies, such as the azimuth, dip angle, and semi-axes, exhibit greater variability, this behavior is expected due to the limited resolving power of gravity data in those dimensions.

The analysis also highlights the method’s ability to provide meaningful estimates of parameter uncertainty through the distribution of accepted models within a defined prediction error threshold. This feature allows for the identification of well-constrained parameters and reveals the anisotropic nature of the inversion’s solution space. As such, the methodology offers a compelling alternative to more computationally intensive probabilistic techniques, while still supporting rigorous uncertainty assessment.

Overall, the proposed framework proves effective for detecting and characterizing subsurface gravity anomalies, with direct applications in geotechnical assessment, environmental studies, and mineral exploration. Future research may extend this approach to incorporate joint inversion with complementary datasets or explore adaptive schemes for parameter reduction to enhance model resolution.