1. Introduction

Efficient and robust computational techniques for the inversion of geophysical datasets are essential for the development of meaningful density and velocity models of the Earth's substructure. Many global optimization methods have been developed and implemented for inverse problems in geophysics (e.g. genetic algorithms ^[1,2,3,4], simulated annealing ^[5,6], neighborhood algorithm ^[7,8], and Monte Carlo methods ^[9]) with the conventional formulation of the inverse problem posed as an unconstrained optimization formulation. Novel computational optimization techniques take advantage of the inclusion of equality and inequality constraints placed on the variables which restrict the model parameters to a smaller range of values shrinking the feasible space ^[10,11]. The use of a priori information reduces model ambiguity and allows us to focus only on those feasible models that provide a good fit to the data while satisfying the physical constraints ^[12]. The standard implementation of equality and inequality constraints uses a Lagrangian method with a powerful constrained optimization method known as Primal-Dual Interior-Point (PDIP) ^[11]. This framework has been an alternative to solve linear programming problems for single and joint inversion of different geophysical data sets ^[11,13,14] but not for the inversion of gravity anomalies in a constrained optimization framework.

Gravity and magnetic fields, also known as "potential fields", are widely used in the exploration of the Earth's substructure. Changes in the gravitational field constitute a passive source of data used to determine the structure of the Earth's subsurface by sensing the distribution of density in rocks. Gravity surveying constitutes a cheap, non-invasive, and non-destructive remote sensing method that helps associate the variations in the gravity field with subsurface density distributions and ultimately rock types. Through the interpretation of the changes observed in a gravitational survey, we can determine the distribution of stratum in a region of study. The juxtaposition of rocks of different densities can help determine the geological processes that have taken place at different locations and be used to locate faults, mineral or petroleum resources, and ground-water reservoirs ^[15].

The gravitational potential energy measures the force of attraction between two bodies (such as the Earth and an object); its strength is proportional to the mass responsible for the gravitational field and inversely proportional to the square of the distance between any part of that mass (usually the center of mass) and the observation point. Surveys using gravitational data focus on minuscule changes in gravity that occur from place to place caused by dissimilarities in the rock density of the subsurface (anomalies). Getting rid of unwanted gravity effects due to known sources, observations from gravity meters get corrected to isolate for regional or local structures (the unknown sources) ^[15]. This paper considers Complete Bouguer anomalies which focus on the effects in gravity corresponding directly to the density contrasts found in the subsurface.

Most methodologies used for the calculation of gravity anomalies can focus on 2-, 2.5-, and 3-D bodies. Available modeling techniques can be extended to 3-D using an approximation to irregularly shaped bodies with several smaller bodies represented by regular shapes whose gravitational field is easier to compute. Published works contain computations for gravity anomalies of basic geometrical forms such as spheres and cylinders ^{[16,17,18,19]} and polygonal prisms ^{[20,21,22,23,24]}, while additional publications examine the gravitational attraction originated by a right rectangular prism ^[25,26]. Now consider our problem: Given a single 2-D profile of a gravitational anomaly, what is the density distribution or shape of a two-dimensional mass with constant density $\rho$ that produces this anomaly? In general, the goal of inverse problems consists of finding the model parameters that allow the optimal reproduction of our set of observed measurements (the 2-D profile of gravity anomalies). Ideally, the exact theoretical framework prescribing how to transform the data in order to reproduce the model exists; in reality, an exact solution may not exist, but it may be sufficient to find the best approximate solution that produces a minimum misfit or residual using a predetermined norm ^[27,28]. A primary goal of gravity anomaly inversion consists of the detection and quantification of changes in the mass properties at different depths ^[29]. Here, even small changes in the model tend to greatly affect the outcome obtained from the inversion. Furthermore, the multi-dimensionality of the different geological density structures makes the problem at hand very complex from a mathematical point of view ^[30].

For many years, the optimization of this inverse problem -- varying the structure and physical properties of the subsurface model until the residual at each station on the surface is minimized -- was based on trial-and-error methods where the shape of the initial starting structure was perturbed. These algorithms, although useful and easy to implement, did not exhaust the possiblity of finding the most optimal model parameters within the solution space. Recently, many potential fields inversion algorithms have been developed ^{[14,29,31,32,33,34,35,36]}. Some of these algorithms use a priori information in the form of positivity constraints for the density values, empirical laws, and upper and lower density bounds; the inclusion of this information in the formulation of the optimization problem produces a reduction of the ambiguity originated from the inherent non-uniqueness. In this paper, we analyze the use of Interior-Point methods to improve the accuracy of the Earth's 2-D density distribution models through the inclusion of physical constraints coming from a priori information obtained from alternative geophysical surveys. We test our approach using different synthetic cases with varying complexity noise levels to show that our constrained optimization can be as accurate as commonly used unconstrained formulations of the problem while satisfying the physical constraints. We discuss the feasibility of our approach through the analysis of our results and conclude with the potential applications of this optimization scheme.

---

2. Methodology

Surveys using gravitational data focus on minuscule changes in gravity caused by dissimilarities in the rock density of the subsurface called gravity anomalies. These anomalies -- the differences between the corrected measured gravity and the theoretical gravity associated with a homogeneous ellipsoid -- originate from the vast multifarious structure of the Earth's interior and are the basis to understand the internal structure of the planet. Starting with the measured absolute gravity values (as recorded with gravimeters on the surface of the Earth), gravity datasets are corrected for all known sources of gravity; after all corrections have been implemented, the resulting dataset corresponds to the gravity anomalies associated with lithology changes based on density contrasts. The Bouguer gravity anomaly is defined as

$\Delta g_B= g_{obs}+ (g_{FA}- g_{BP} + g_{TC})- g_{n},$

(1)

where $g_{obs}$ is the measured absolute gravity value, $g_n$ is the theoretical (or normal) gravity value, and we have the free-air ($\Delta g_{FA}$), Bouguer plate ($\Delta g_{BP}$), and terrain ($\Delta g_{TC}$) corrections ^[37].

---

2.1. Interpretation of gravity anomalies

Gravity anomalies result from the irregular distribution of density within the Earth. The density contrast of a body is the difference between the density of rocks in an anomalous body $\rho$ and the density of the surrounding rocks: $\Delta \rho= \rho -\rho_0$ ^[37]. Anomalous bodies with higher density than host rock have positive density contrasts while bodies with lower density than host rock have negative density contrasts. A positive gravity anomaly is obtained over high-density bodies where the measured gravity is augmented; likewise, negative anomalies result over low-density bodies ^[37]. Therefore, analyzing the sign of an anomaly helps to determine the sign of the density contrast and whether the density of the geological body should be higher or lower than the surrounding rocks. The apparent "wavelength" of an anomaly refers to its horizontal extent and can be used to determine the depth of its anomalous mass ^[37]. Usually, a Bouguer gravity anomaly map contains superposed anomalies coming from several sources (in a 3-D sense). Long-wavelength anomalies (caused by deep density contrasts) called regional anomalies relate to large-scale structures of the Earth's crust such as mountain ranges, oceanic ridges, and subduction zones. Short-wavelength anomalies (due to shallow density constrasts) called residual anomalies often relate to shallow anomalous masses such as near-surface mineralized bodies ^[37]. Although this can be helpful to model the density distributions of the subsurface and determine the approximate location of the sources (large deep bodies often associate to broad long-wavelengths and low-amplitude anomalies while shallow bodies associate with narrow short-wavelength and sharp anomalies), a priori information from additional geophysical surveys is required to resolve ambiguities from non-uniqueness. After the removal of regional anomalies, residual gravity should be interpreted in terms of its approximate density distribution using iterative 2- and 3-D techniques.

---

2.2. 2-D modeling of gravity anomalies

We used a 2-D polygonal prisms algorithm ^[23,26,38] as our forward operator for the calculation of the gravity anomalies associated with 2-D density structure models (or profiles). Each density structure model contains a variety of $n$-sided polygons depicting the geometry of each body of constant density found in the subsurface. We used these polygons projected in the $y$ axis (in and out of the page) based on the geology of the deposit in order to reproduce the gravity anomalies recorded for each gravity station. The Bouguer anomaly associated with the density structure model is a function of the density contrasts, the geometrical shapes of the bodies ($x$ and $z$-coordinates of its vertices), the depths of the bodies, and the location of the gravity stations ^[39,40].

The total gravity anomaly associated with a gravity station is the sum of all contributions of all the vertices of all the polygons used to illustrate the substructure ^[23,26,38]. Placing each gravity station at the origin of an $xz$ coordinate system and following the geometrical convention shown in Figure 1 (in a strict clock-wise direction), we express the total vertical component of gravity anomaly at station $P$ as

$\Delta g_z = 2 \gamma \sum^l_{j=1}\bigg(\Delta \rho_j \sum_{i=1}^{n_j} B_i\bigg[(\theta_i - \theta_{i+1}) +\frac{(z_{i+1}-z_i)}{(x_{i+1}-x_i)}\mbox{ln}\bigg(\frac{\sqrt{{x_{i+1}}^2+{z_{i+1}}^2}}{\sqrt{{x_i}^2+{z_i}^2}}\bigg) \bigg] \bigg)$

(2)

with

$B_i= \frac{(x_{i+1}-x_i)(x_i z_{i+1}- x_{i+1} z_i)}{(x_{i+1}-x_i)^2+(z_{i+1}-z_i)^2} \mbox{ }\mbox{ and } \theta_i= \mbox{tan}^{-1}\bigg(\frac{z_i}{x_i}\bigg)$

where $\gamma$ is the gravitational constant $6.673848 \times {10^{ - 11}}{\rm{ }}{{\rm{m}}^3}/{\rm{kg }}{{\rm{s}}^2}$, $l$ is the number of bodies in the density structure model, $\Delta \rho_j$ is the density contrast of the $j$th polygon, and $n_j$ is the number of sides of the $j$th polygon.

It is important to distinguish between the origin of the actual density structure and the origin used for the calculation of gravity anomalies. As stated before, in order to calculate the gravity anomaly for a given station, it must be placed at point $(0, 0)$ such that all distances from this point to the vertices of all bodies can be easily determined following the scheme in Figure 1. On the other hand, the location of the origin for the actual density structure can be placed at the beginning of the profile on the surface or a base station where absolute gravity is known. For our purposes, the origin is placed at the far left side on the surface of our 2.5-D density profiles.

Using forward operator (2), we calculate the gravity anomaly associated with the given density structure model. The calculation of the density distribution in the Earth's substructure based on the measurements of its gravity field constitutes an example of an inverse problem in which we infer the causal factors that produce a set of observations.

---

2.3. Inversion of gravity anomalies

Most inverse methods approximate the values to the parameters of the source corresponding to the observed anomaly by iteratively solving equations related to the forward operator for the anomalies and inherent physical laws until $\Vert\Delta g_{obs}-\Delta g_{calc}\Vert \leq \varepsilon$ mGals. There are three important aspects that must be considered when solving inverse problems: solution existence, solution uniqueness, and instability of the solution process ^[27]. In theory, the calculated structure of the bodies associated with a gravity anomaly signature obtained from the inversion of surface observations would correspond to the "real" geological substructure of a region. However, the discrete nature of the observational data, the errors incurred in mathematical methods, equipment, and data corrections, and the assumptions made in geological and geophysical settings makes finding mathematically acceptable answers to inverse problems a non-trivial task. The ambiguity of geophysical data can lead to geologically unfeasible models that show a great level of agreement between observed and computed data, which translates to infinitely many models that fit the data in an adequate way ^[41,42]. Instability relates to ill-posedness where small alterations to the data result in large changes in the inferred models ^[27].

Using nonlinear integral equations, we assume the density contrasts associated with the 2-D density structure models and solve our system of equations for the geometrical characteristics (depths from the surface to each of the vertices in the polygonal prisms) associated with our synthetic multi-interfaced Earth structures. Assuming information regarding the mass distribution (source bodies' configurations), we assemble an initial structural model to start the inversion. Additional independent geophysical datasets can be used to constrain theoretical density models, reduce the solution space, and focus only on geologically feasible models. The inclusion of upper and lower bounds on the inversion parameters help us to find feasible models that will suit our needs and that conserve the accuracy, feasibility, and consistency associated with well-behaved inverse problems ^[27].

---

2.4. Using Primal-Dual Interior-Point methods

The adapted characterization of the constrained optimization framework solves a linearized version of the inverse problem by adding bound constraints over the model parameters. We define $Z$ as a density structure model with its elements specifying the location of the vertices of the geometrical shapes (depths or $z$-coordinates of each vertex) delineating the various bodies in the subsurface. Physical constraints preserve the relationships between the different bodies, the consistency of their boundaries, and the integrity of the final density structure models. These restrictions relate to the appropriate geometrical bounds pertaining the model parameters $Z$ in the inverse problem. Since each $Z$ in the parameter space relates to a geological model, each minimum point located within the solution space represents a possible structure.

From gravity prospecting, we create a hypothetical density structure model $Z\in \mathbb{R}^n$ and evaluate the non-linear forward operator $G \in \mathbb{R}^m$; $G(Z)$ represents the Earth's response or the calculated gravity anomaly. The relationship between the forward operator $G$ and the density structure dataset is given by

$G(Z)= [G_1(Z), G_2(Z), \ldots, G_m(Z)], \mbox{ for }Z=(z_1, z_2, \ldots, z_n),$

where $m$ is the number of observations and $n$ is the total number of vertices in the density distribution. For any dataset of gravity anomaly observations $A \in \mathbb{R}^m$, the inverse problem consists of finding the unknown density structure model $Z$, which allows the best approximation of $G(Z)$ to $A$, that is,

$\underset{Z}\min \mbox{ }\frac{1}{2}\Vert G(Z)-A\Vert^{2}= \underset{Z}\min \mbox{ }\frac{1}{2}\sum_{i=1}^{m}\left(G_{i}(Z)-A_{i}\right)^{2},$

(3)

which is usually posed as an unconstrained weighted non-linear least squares (NLSQ) problem.

The complexity of the non-linear operator $G$ makes the use of an iterative linearized least squares approach an alternative to avoid the computation of higher order derivatives in the minimization problem. In least squares problems, the objective function $F$ has the form

$F(Z)= \frac{1}{2} {\Vert R(Z) \Vert_2}^2,$

with the residual vector $R:\mathbb{R}^n \to \mathbb{R}^m$ defined as $R(Z)= G(Z)-A$.

The best match between model and data comes from minimizing function $F(Z)$ with respect to the desired parameters. We implement the constrained optimization framework proposed by ^{[11,43,44,45]} to solve the linearized version of the non-linear inverse problem. The constrained optimization problem is given by

$\underset{Z}\min \mbox{ } F(Z) \\ s.t. \mbox{ }h_E(Z)=0 \\ \mbox{ }h_I(Z) \geq 0,$

(4)

where $h_I(Z) \in \mathbb{R}^p$ are the inequality constraints and $h_E(Z)\in \mathbb{R}^q$ are the equality constraints associated to $Z$. The inequality constraints $h_I(Z)$ represent the physical bounds that correspond to the minimum and maximum depths of the vertices of the polygonal prisms and/or relationships between the different vertices used to conserve the integrity of the structural model. The equality constraints in $h_E(Z)$ represent the appropriate end-body conditions and the lower and upper boundaries for the top and bottom bodies in the cross-sections.

Problem (4) is redefined in the standard non-linear programming form

$\underset{Z}\min \mbox{ } F(Z) \\ s.t. \mbox{ }h_E(Z) =0 \\ \mbox{ }h_I(Z)-s =0 \\ s \geq 0,$

(5)

with $s \in \mathbb{R}^{p}$ so called slack variable. We define $S=\mbox{diag}(s_1, s_2, \ldots, s_{p})$ (a matrix with the elements of $s$ in the diagonal).

We solve (5) using Primal-Dual Interior-Point (PDIP) methods, starting with the Lagrangian function

$\mathcal{L}(Z, s, y_E, y_I)= F(Z) - {h_E}^T(Z) y_E-\big({h_I}(Z)-s\big)^Ty_I,$

(6)

where $y_E\in \mathbb{R}^{q}$ and $y_I\in \mathbb{R}^{p}$ are the Lagrange multipliers associated to the equality and inequality constraints, respectively, and $(s, y_I)>0$.

We define the Karush-Kuhn-Tucker (KKT) or necessary conditions for the optimization of this nonlinear programming problem as

$\begin{array}{r} \nabla F(Z) - \big(\nabla{h_E}(Z)\big)^Ty_E - \big(\nabla{h_I}(Z)\big)^Ty_I =0, \\ Sy_I=0, \\ h_E(Z)=0, \\ h_I(Z)-s=0. \end{array}$

The second condition, $Sy_I=0$, also known as the complementarity condition for the minimization problem, implies that one of the components of each product $\big(s_i\cdot({y_I})_i\big)$ must be equal to zero for each $i=1, 2, \ldots, p$ (and nonzeros of $s$ and $y_I$ should appear in complementary locations) ^[44,45].

We address the nonlinear programming problem with PDIP methods solving a sequence of approximate linear minimization problems iteratively ^[44]. As a first step we include a perturbation parameter $\mu$ to the complementarity condition (and use the notation $\nabla$ as the gradient operator) such that the perturbed Karush-Kuhn-Tucker (PKKT) conditions become

$\left[ \matrix{ \nabla z{\cal L}(Z, s, {y_E}, {y_I}) \hfill \cr S{y_I} - \mu e \hfill \cr \nabla {y_E}{\cal L}(Z, s, {y_E}, {y_I}) \hfill \cr \nabla {y_I}{\cal L}(Z, s, {y_E}, {y_I}) \hfill \cr} \right] = \left[ {\matrix{ {\nabla F(Z) - \nabla {h_E}^T(Z){y_E} - \nabla {h_I}^T(Z){y_I}} \cr {S{y_I} - \mu e} \cr {{h_E}(Z)} \cr {{h_I}(Z) - s} \cr } } \right] = \left[ {\matrix{ 0 \cr 0 \cr 0 \cr 0 \cr } } \right],$

(7)

where $e=(1, 1, \ldots, 1)\in \mathbb{R}^{p}$, and $(s, y_I)>0$.

Interior-point methods solve the PKKT conditions for a sequence of positive parameters ${\mu_k}$ that converges to zero while $(s, y_I)>0$ ^[44]. Letting sequence $\mu$ be strictly positive forces variables $s$ and $y_I$ to stay strictly positive, keeping the iterates away from the boundaries and in the interior of the constraints while converging to the optimal solution $(Z^*, s^*, y_E^*, y_I^*)$ as $\mu \to 0$ and satisfying the KKT conditions for (5) ^[44]. By respecting the bounds, interior-point methods avoid spurious solutions -- those minimizing the objective function but not meeting $(s, y_I)>0$ ^[45].

Applying Newton's Method, the primal-dual system associated with (7) becomes

$\left[ {\matrix{ {{\nabla ^2}{{\cal L}_{ZZ}}} & 0 & { - \nabla {h_E}^T(Z)} & { - \nabla {h_I}^T(Z)} \cr 0 & {{Y_I}} & 0 & S \cr {\nabla {h_E}(Z)} & 0 & 0 & 0 \cr {\nabla {h_I}(Z)} & { - I} & 0 & 0 \cr } } \right]\left[ {\matrix{ {\Delta Z} \cr {\Delta s} \cr {\Delta {y_E}} \cr {\Delta {y_I}} \cr } } \right] = - \left[ {\matrix{ {{\nabla _Z}{\cal L}(Z, s, {y_E}, {y_I})} \cr {S{y_I} - \mu e} \cr {{h_E}(Z)} \cr {{h_I}(Z) - s} \cr } } \right],$

(8)

in the variables $Z$, $s$, $y_E$, $y_I$ with $(s, y_I)>0$. Here, $Y_I=\mbox{diag}(y_1, y_2, \ldots, y_{p})$, $\mathcal{L}$ denotes the Lagrangian (6), and $\nabla ^2 \mathcal{L}_{ZZ}$ is the Hessian of the Lagrangian.

The primal-dual system can be rewritten in symmetric form

$\left[ {\matrix{ {{\nabla ^2}{{\cal L}_{ZZ}}} & 0 & {\nabla {h_E}^T(Z)} & {\nabla {h_I}^T(Z)} \cr 0 & \Sigma & 0 & { - I} \cr {\nabla {h_E}(Z)} & 0 & 0 & 0 \cr {\nabla {h_I}(Z)} & { - I} & 0 & 0 \cr } } \right]\left[ {\matrix{ {\Delta Z} \cr {\Delta s} \cr { - \Delta {y_E}} \cr { - \Delta {y_I}} \cr } } \right] = - \left[ {\matrix{ {{\nabla _Z}{\cal L}(Z, s, {y_E}, {y_I})} \cr {{y_I} - \mu {S^{ - 1}}e} \cr {{h_E}(Z)} \cr {{h_I}(Z) - s} \cr } } \right],$

(9)

where $\Sigma=S^{-1}Y_I$.

In order to ensure that the function value is decreasing at each step towards the solution $Z^*$, we can use a merit function to determine whether a step is productive and should be accepted ^[44]. Our merit function has the form

$\phi_{ {v}} = \bigg(F(Z)- \mu \sum_{i=1}^p\mbox{log}(s_i) \bigg)+ {v} \Vert h_E(Z) \Vert + {v} \Vert h_I(Z)-s \Vert,$

(10)

where ${v}$ is a penalty parameter used to force the solution towards feasibility (attempting a shorter step whenever the original step fails to decrease the merit function). Using $\mu>0$ and the logarithmic barrier term prevents the component of $s$ from getting too close to zero (since each $-\mu \mbox{log}(s_i)\to \infty$ as $s_i\to 0$) ^[44].

PDIP methods characterize by having a procedure to determine the step to take in an iteration and a measure of desirability of points in the search space ^[44]. At each iteration, the system advances to the next solution by taking the step $\Delta d= (\Delta Z, \Delta s, - \Delta y_E, - \Delta y_I)$ from the current point (also called the Newton direction); however, if the condition $(s, y_I)>0$ is violated, the full step is not a feasible or desired step ^[45]. A way around this problem consists in performing a line search along the Newton direction such that the new iteration becomes

$Z_{k+1}= Z_k+ \alpha \Delta Z, \\ s_{k+1}= s_k + \alpha \Delta s, \\ {y_E}_{k+1}= {y_E}_k+ \alpha(- \Delta y_E), \\ {y_I}_{k+1}= {y_I}_k+ \alpha(- \Delta y_I),$

(11)

for some $\alpha \in (0, 1]$ called the line search parameter.

We use parameter $\tau \in (0, 1)$ (usually equal to 0.995) to avoid moving to the boundaries too quickly ^[44] and define

${\alpha_{s}}^{\max}= \max \{\alpha \in (0, 1]: s+\alpha \Delta s \geq (1-\tau) s\}, \\ \mbox{ and } {\alpha_{y_I}}^{\max}= \max \{\alpha \in (0, 1]: y_I+\alpha \Delta y_I \geq (1-\tau) y_I\}.$

Using our merit function (10), we compute step lengths ${\alpha}_s \in (0, {{\alpha}_s}^{\max}]$ and ${\alpha}_{y_I}\in(0, {{\alpha}_{y_I}}^{\max}]$ that satisfy

$\phi_{ {v}}(Z_k +{\alpha}_s \Delta Z, s_k+\alpha_s \Delta s) \leq \phi_{ {v}}(Z_k, s_k) + \eta \alpha_s {D\phi}_{ {v}}(Z_k, s_k; \Delta Z, \Delta s)$

(12)

The implementation of the line search and (12) ensures finding the values of $(\alpha_s, \alpha_{y_I})$ that guarantee a sufficient decrease of the merit function ^[44]. The new iteration values are defined as

$Z_{k+1}= Z_k+ {\alpha_s}\Delta Z, \\ s_{k+1}= s_k+ {\alpha_s}\Delta s, \\ {y_E}_{k+1}= {y_E}_k+ {\alpha_{y_I}}\Delta y_E, \\ {y_I}_{k+1}= {y_I}_k+ {\alpha_{y_I}}\Delta y_I.$

(13)

Our algorithm contains optimized MATLAB functions to solve the approximated problem (10) using either a direct step (with Newton's method) or a Conjugate Gradient (CG) step. The direct step uses an LDL factorization of the matrix if the Hessian is positive definite; otherwise, the algorithm uses a CG step. In the CG approach, the problem is posed as minimizing a quadratic approximation to the barrier problem in a trust region subject to linearized constraints (more in ^[44]).

Gravity problems usually contain a large number of local minimum points each of them related to a geological structure (either well-defined or ill-defined) ^[46]. Using constraints for the parameter values and incorporating them into the optimization formulation allows us to focus only on those points within the feasible region that minimize our objective function.

---

3. Results

As previously stated, gravity data is mainly sensitive to the density distribution of anomalous bodies and their locations. Assuming constant values for the density contrasts located in the area, the forward operator $G$ used in the optimization algorithm depends non-linearly in the model parameters, $Z_i$ (for $i=1, \ldots, n$) that represent the depths from the surface (at sea-level) to the $n$ vertices of the polygonal prisms.

The inversion begins with an initial guess of the density structure model, $Z_0$, based on a priori geological information available for the area of study. Once the gravity anomaly has been calculated for the initial density structure model ($G(Z_0)$) and compared to the observation vector $A$, we obtain the residual vector $R(Z_0)$ and iterate the evaluation of the inverse problem to obtain different approximations of $Z_k$ using our constrained optimization method. The algorithm runs until it meets one of the stopping criteria: 1) the residual error is less than a tolerance $\varepsilon_1$, 2) the relative error (rms) is less than a tolerance $\varepsilon_2$, 3) the maximum number of iterations is reached, or 4) the difference between two consecutive iterations is less than a tolerance $\varepsilon_3$. The last iteration returns the latest updated model $Z_{k}$ which represent the best density structure model -- since the line search and sufficient decrease schemes included in the algorithm guarantee that the objective function $F(Z)$ decreases at each iteration. Our goal consists in finding the optimal $Z_k$ that resembles the hypothetical model as closely as possible (e.g. ${\rm{rms = }}\left\| {{{{Z_k} - {Z^*}} \over {{Z^*}}}} \right\|{\rm{ ＜ }}{\varepsilon _3}$) ^[11]. Moreover, at an optimal $Z_k$, the residual error between the predicted gravity anomaly (corresponding to the final density structure model from the inversion) and the observations (or gravity anomaly values associated to a hypothetical model) is less than a given tolerance (e.g. ${\rm{residual error = }}\left\| {{{G({Z_k}) - A} \over A}} \right\|{\rm{ ＜ }}{\varepsilon _1}$) and all the elements of $Z_k$ strictly meet the constraints (guaranteeing feasibility) ^[11,46].

---

3.1. Comparison with non-linear least squares

We compare the behavior of the proposed PDIP methods and the usual unconstrained non-linear least squares for the optimization of Bouguer gravity anomalies using the synthetic dataset for a simple structure. We used a total of 11 gravity observation points on the surface separated every kilometer. The structure consisted of two contiguous bodies: body B1 with density contrast 1000 Kg/m$^3$ and composed of tree vertices, and body B2 with density contrast 1870 Kg/m$^3$ and four vertices. However, we recorded each vertex only once and used these "unique" vertices as parameters for the optimization (adjusting the calculations of individual contributions of vertices to the total gravity anomaly and carrying out the calculations based on the unique vertices). Throughout the inversion the $x$-coordinates remained constant while the $z$-coordinates moved freely vertically. We had inequality conditions associated to $\{\forall z_i \in Z: z_i \in [-5, 0]\}$ and no equality constraints. Figure 2 shows the gravity profiles and 2-D density distributions of (a) the initial structural model $Z_0$, (b) the final structural model obtained using an unconstrained non-linear least squares approach, and (c) the final structural model resulting from our constrained PDIP method approach.

Figure 2 shows both methods seemingly recovering the actual gravity anomaly values recorded at each observation point. The residual error (a measure of the differences between final calculated and observed gravity values) is $1 \times {10^{ - 9}}$ for the structural model recovered by the PDIP algorithm and $4 \times {10^{ - 2}}$ for the NLSQ algorithm; however, only the PDIP method recovers the exact observations associated with the target density structure model. Comparing the final density structure models for both methods to the actual parameter values $Z^*$ (associated with the target geological structure), the relative error (rms) related to $Z_{PDIP}$ is $1 \times {10^{ - 9}}$, while that of the $Z_{NLSQ}$ is 1.7584. Moreover, from the bottom parts of Figure 2(b) and 2(c) we can see that some of parameters in $Z_{NLSQ}$ are outside the boundary conditions (a spurious solution) and far away from the values of $Z^*$ (dotted red lines included in xz plane at the bottom of Figure 2(b) and 2(c)) while $Z_{PDIP}$ coincide exactly with $Z^*$.

An additional advantage of the optimization with PDIP methods is its fast convergence to the final model. Figure 3 shows the number of iterations versus the value of the objective function for both the NLSQ and the PDIP methods associated with this synthetic example. While the NLSQ algorithm takes almost 600 iterations to converge to a final density structure model (since it looks at all possible solutions even if they lie outside the boundaries) with a final objective function value of 3.3205, the PDIP algorithm takes 20 iterations to converge to a solution with a final objective function value of $6.7190 \times {10^{ - 9}}$.

From our results, the inclusion of inequality constraints in our formulation forced the inversion using PDIP methods to stay in the feasible region (below the surface) while minimizing the residual between observed and calculated gravity anomalies -- a feature not accomplished by the NLSQ approach whose final structural model was a spurious solution outside the physical parameter constraints.

With respect to the equality and inequality constraints for the inversion parameters, special care must be placed on the formulations to avoid excessive severity that may exclude potential solutions. Moreover, for a more complex density structure model, it would be necessary to include constraints to preserve the structural composition of each one of its bodies and accurately represent the appropriate formation features. Figure 4 illustrates an example where the structural composition of a body has been compromised; this occurs when a vertex (in this case vertices 2 and 9) cross-cuts a side of a body beyond the limits imposed by its corresponding vertices (e.g., vertex 2 of model Final Zf is smaller than vertices 11 and 12 and cuts the line formed by them and vertex 9 is larger than vertices 4 and 5 and goes beyond the line formed by them). These conditions can easily be included in the formulation for PDIP methods using boundary constraints; however, the situation differs for NLSQ where such conditions should be excluded from the optimization and final density structure models for very complex areas suffer this kind of problem.

The representations used to depict different geological anomalous bodies may contain any finite number of vertices; using many vertices to represent a body can help to better model the structure while too small a number can compromise the final result or produce unacceptable solutions. However, total computation time increases with the number of vertices, hence, the trade-off should be considered when prioritizing desired outcomes from the optimization.

---

3.2. Synthetic examples

We implemented our constrained optimization using PDIP methods for three different synthetic models with varying complexity. We start with a simple geometrical representation for a sedimentary basin and move to two multi-interfaced continuous structures with higher intricacies in their density distribution and body inter-dependence; all synthetic datasets were created noise free.

---

3.2.1. Sedimentary basin

The sedimentary basin model has densities $2670$ kg/m$^3$ for granites and $2960$ kg/m$^3$ for greenstones. A total of 51 gravity observation points every kilometer helped us to determine the density distribution of the subsurface. We used 47 vertices to model bodies B1 and B2 and endbodies B3 and B4 (included in the density structure models to improve accuracy at the endpoints). 27 non-repeated vertices out of the total were used as parameters for the inversion. Throughout the inversion, only the $z$-coordinates varied freely. Bounds $\{\forall z_i \in Z: z_i \in [-2, 0]\}$ were used as inequality constraints; the equality constraints were related to the vertices on the surface and end-body structures and may have forced $z$-coordinates to remain constant during the inversion to ensure feasible boundary representations for all bodies. A base station with known absolute gravity was situated at $x=30$ Km and used to determine the gravity anomaly of the rest of the stations (the difference between the hypothetical and calculated gravity anomalies at this point is 0). In real applications base stations can be placed anywhere in the profile (by using well-known base stations or in-site base stations) and are important features used to ensure the accuracy of the calculated gravity anomalies.

Figure 5 includes: (a) the hypothetical model of the sedimentary basin and the gravity observations corresponding to each gravity station, (b) the initial 2.5-D model for the basin (usually based on a priori information but created randomly for our synthetic examples), its gravity signature, and hypothetical gravity observations, and (c) the final model obtained using our PDIP method as part of the optimization, its gravity signature, and the hypothetical gravity observations (densities are in kg/m$^3$).

The final density structure model obtained from the inversion, $Z$, recovers the hypothetical density structure model $Z^*$ almost exactly while also minimizing the objective function ($F(Z)=1.77 \times {10^{ - 3}}$). In fact, the relative error (rms) of $Z$ with respect to model $Z^*$ is $5.92 \times {10^{ - 2}}$ and the residual error between the final calculated gravity anomalies ($G(Z)$) and the observations (in this case $G(Z^*)$) equals $6 \times {10^{ - 4}}$.

From our results, we show that the algorithm converged to the optimal solution of the minimization problem when starting the inversion with a particular density structure model $Z_0$ (Figure 5(b) bottom). However, considering the non-uniqueness inherent to gravity inversion problems, starting the algorithm with a different $Z_0$ may result in different values for the final $Z_f$ that minimizes the objective function. Figure 6(a) shows the values of $Z_0$ for the vertices associated with ten different initial density structure models. Figure 6(b) contains the $Z_f$ values (and the resulting density distributions) derived from the optimizations. All inversions were tied to the base station at $x=30$ Km -- such that $|G(Z_0)-G(Z^*)|=0$ for all $Z_0$.

From Figure 6, we can observe that most inversions converged to the actual parameter values $Z^*$ (Figure 6(b)) even though they started at different initial models $Z_0$ (Figure 6(a)). Those final density structure models $Z_f$ that did not converge to $Z^*$ still had objective function values $F(Z_f) ＜1 \times {10^{ - 2}}$. All parameters associated with the final models met the boundary and the structural constraints, therefore, even though some of the initial values of the parameters of the ten $Z_0$ were outside the boundaries, all final values associated with all the inversions lie within the feasible region.

Figure 7 shows the ten initial calculated gravity profiles, $G(Z_0)$, and the final calculated gravity profiles, $G(Z_f)$, for the density structure models derived from the inversions; notice all $G(Z_0)= G(Z^*)$ at $x=30$ Km. Although none of the initial calculated gravity profiles associated with the ten initial density structure models coincide with the observed gravity anomaly at each station ($G(Z^*)$ or $A$) as shown in Figure 7(a), all the final calculated gravity anomaly profiles (corresponding to the final density distributions $Z_f$ in Figure 6(b)) are very close; therefore, all the final density structure models lie in the feasible space (a subspace of the solution space where all equality and inequality constraints hold), minimize the objective function, and recover the actual (or hypothetical) density structure model to some extent. A summary of the initial and final objective function values, and relative and residual errors of these implementations is included in Table 1 (only the inversions with highest final objective function values are shown).

---

3.2.2. Multiple-bodies continuous structure

The next structure consisted of seven inter-dependent bodies with density contrasts in the range $[-200, 700]$ Kg/m$^3$. We used 21 surface gravity observation points every 100 meters and modeled bodies B1 through B7 and endbodies B8 and B9 with a total of 43 vertices; 24 unique vertices served as parameters for the inversion where only the $z$-coordinates varied freely. We apply the following inequality constraints: $\{Z+1 \geq 0, -Z \geq 0\}$. Additional equality and inequality constraints were included to preserve structural integrity in all the bodies during the inversion. Figure 8 shows (a) the hypothetical model and its actual calculated gravity anomaly ($A^*$) and (b) the initial and (c) the final density models related to the optimization and their corresponding gravity anomalies. No base stations were used in this example.

The final density structure model from the inversion $Z_f$ did not recover the hypothetical density distribution model exactly (bottom of Figure. 8(c) and (a) respectively) even though the residuals between the observed and calculated gravity anomalies were minimal (top of Figure 8(c)). Although lateral bodies B1 and B7 have been completely recovered, the rest of the bodies were close but not quite the expected result.

We implemented our algorithm using 10 different initial models $Z_0$ for this synthetic example; Figure 9 illustrates the values of all vertices associated with the initial and final density structure models. Some of the parameters from the different $Z_0$ lie outside the constraints while all parameters of the final $Z$ meet the constraints and preserve the structural integrity of all the bodies (Figure 9(b)). Moreover, all final density structure models are guaranteed to lie in the feasible space (meeting the equality and inequality constraints), hence, no spurious solutions are included in our results. Although the exact density structure $Z^*$ was not recovered by the final inversion models $Z_f$, most of the parameter values ended up clustered together close to the actual parameter values.

Figure 10 contains the initial and final calculated gravity profiles ($G(Z_0)$ and $G(Z_f)$ respectively) for all inversions. Figure 10(a) shows that none of the initial calculated gravity profiles matched the gravity anomaly associated with the actual $Z^*$. Figure 10(b) shows that all inversions minimized the objective function value $F(Z)$ -- given that all $G(Z_f)$ coincide everywhere with the actual gravity profile ($G(Z^*)$ or $A$).

Table 2 contains a summary of the inversion results associated with the ten implementations; only the results for the five implementations with the highest final objective function values ($F(Z_f)$) are included.

---

3.2.3. Faulted structure

The last synthetic structure consisted of six inter-dependent bodies representing a 3 layer fault with densities in the range $[0, 2000]$ Kg/m$^3$. We used 21 surface gravity observation points separated every kilometer. Bodies B1 through B6 and seven endbodies were modeled using 112 vertices; the inversion had 55 unique vertices as variables. Only the $z$-coordinates varied freely within the feasible space given by $\{\forall z_i \in Z: z_i \in [-5, 0]\}$. We included additional equality and inequality constraints in the formulation to ensure structural integrity for all bodies in the density distribution model. Figure 11 shows the hypothetical, initial, and final density structure models and their corresponding gravity anomaly profiles for a layered fault structure. The initial model was obtained analyzing the structure of the hypothetical model and adding random noise. A base station was included at $x=10$ Km where $G(Z^*)=0$ mGals.

Analyzing the final 2-D density structure model recovered from the inversion, we can see that it differs slightly from the hypothetical density model (bottom parts of Figure 11(c) and Figure 11(a) respectively), even though the hypothetical and final calculated gravity anomalies ($G(Z^*)$ and $G(Z_f)$) are practically the same. The objective function decreased from $F(Z_0)= 2549.96$ to $F(Z_f)= 1.61 \times {10^{ - 3}}$, while the relative and residual errors associated with the final density structure model were equal to $4.18 \times {10^{ - 1}}$ and $2.0 \times {10^{ - 3}}$ respectively. Therefore, although the optimization algorithm did not recover $Z^*$ exactly, it minimized the objective function by converging to a feasible solution $Z$ close to it.

The results of the implementation of the optimization algorithm using ten different initial density structure models $Z_0$ are shown in Figure 12. The boundaries and their corresponding equality and inequality constraints remained the same for all the inversions. Note that some of the parameter values associated with the initial density structure models $Z_0$ lie outside the boundaries while the final parameter values associated with all $Z_f$ lie within the reduced solution space defined by the boundaries. Figure 12(b) shows the ten final $Z_f$ models obtained from the inversions. Analyzing the parameter values associated with these density structure models, we see that all the final values for $Z_f$ obtained from the inversions meet the constraints and are actually close to $Z^*$ (black dots in the figure) although none of them gets to the true parameter values. Therefore, most of the density structure models resulting from the inversion did not recover the synthetic density structure even though their final objective function values were close to zero. Looking carefully at Figure 12(b), we see that most of the final parameter values associated to $Z_f$ converge to the same solution and, in fact, form small clusters within the model of the structure; therefore, although we do not have convergence to the actual model, the majority of our inversions converged to the same values even when starting from different $Z_0$. Table 3 shows a summary of the inversion results.

Figure 13 illustrates the differences between the calculated gravity profiles associated to all $Z_0$ in Figure 12(a) and compares them to the observed dataset $G(Z^*)$ associated with the hypothetical density distribution. Note they are different at all points except at $x=10$ Km where a base station with known absolute gravity value is located. All the final calculated gravity profiles for the ten final $Z_f$ models coincide with $F(Z^*)$ at all points; hence, the objective function $F(Z)$ has been minimized by all the inversions.

Table 3 and Figure 12 and 13 show that although the objective function was minimized for all the different initial density structure models $Z_0$, the final structures obtained from the inversions did not recover the hypothetical (actual) density structure model $Z^*$ shown in Figure 11(a). However, the use of constraints in the formulation of our problem allowed the algorithm to converge to similar final density structure models as illustrated in Figure 12(b) and avoided convergence to unfeasible spurious solutions.

---

3.3. Noisy synthetic examples

The most common interferences in gravitational datasets are often caused by ambient, geologic, and cultural conditions (e.g. spatial variations in density, earthquakes, earth tides, extreme temperatures, vibration of vehicles, heavy equipment, etc.). In order to show how noise would affect our results, we added realistic noise (random Gaussian noise) to the vector of gravity anomaly observations $A$ and checked the robustness of the method. We used noise levels of 2.5 %, 5 %, and 10 % to the observed gravity anomalies $A$ (or $G(Z^*)$ for our synthetic models). Initial density structure models and all equality and inequality constraints remained the same as in part (b) in Figure 5,8, and 11. Figures 14 through 16 illustrate the numerical results for five tests performed for each noise level: parts (a), (c), and (e) show the final density structure models recovered with the resulting $Z$ from all inversions while parts (b), (d), and (f) show the final calculated gravity anomalies associated with the resulting models. Some of the results obtained from these implementations using 10 % noise are summarized in table 4.

Our results show that the density structures were well characterized by the PDIP method aided with the use of constraints in the inversion formulation. Out of the three structures, the sedimentary basin was the most affected by the addition of random noise; the range of gravity anomaly values for this cross-section was between $[-110, -70]$ mGals such that the noise represents a maximum of 11 mGals for each station (the maximum noise added to each station for the multi-bodied structure was around 1.8 mGals and around 3 mGals for the faulted structure).

From Table 4, we can see that even though the values of the final optimization function $F(Z_f)$ for all the implementations are very high, the ones for the sedimentary basin are the highest. The sedimentary basin representations recovered by the inversion algorithm for all noise levels are close to the actual density structure given by $Z^*$; moreover, their corresponding final calculated gravity anomalies are very close to the actual values $G(Z^*)$. However, for the sedimentary basin examples the final calculated gravity anomalies were not close to the noisy observations even though they were close to the original observations (without noise).

Even though the complexity associated with both the multi-bodied and faulted-layers models is higher with respect to the sedimentary basin model (since there are more bodies and they depend on each other) the use of physical and structural constraints in the parameter values of $Z$ helped us to stabilize the algorithm; keeping the final density structure models in the feasible region and allowing variation in the parameters only within the boundaries. Without these bounds, the parameters vary freely within the whole region and may converge to a spurious solution. Therefore, although noise in our gravity observations affects our results (since our final density structure models $Z_f$ are very different from each other), parameter constraints help us to ensure feasibility and structural composition at all times.

---

3.4. Additional Results

In order to determine the sensitivity of the algorithm to changes in the initial density structure models, we are including the mean and standard deviations associated with the final models obtained in sections 3.2 and 3.3. Parts (a)-(f) in Figure 17 through 19 include: (a) final Z values ($Z_f$) versus Initial Z values ($Z_0$) obtained from the inversion of 30 different initial density structure models -- '+' signs of the same color represent all the values associated with the same inversion parameter (variable) while vertical gray lines are the actual values for all those parameters -- and (b) the mean and standard deviation calculated from all 30 inversions compared to the actual $Z^*$ values.

In part (a) of these figures, the more spread out in the y-axis (vertically) a variable is the higher the change in the initial density structure experimented in the parameter. As an example, consider all the values associated with the final parameter of $Z_f$ at $-1.1$ Kms in Figure 17(a); the values of all used $Z_0$ at this particular point ranged between $[-2.8, 0.55]$ while the final values associated with each one of the inversions ($Z_f$) for this parameter range between $[-1.15, -1.05]$ (however, most of the values lie on the exact value or vertical gray line).

With respect to part (b) of these figures, we see that even when we start from different initial density structure models, the mean values (values associated with the mean density structure model $Z_{M}$ of all 30 final $Z_f$) of the parameters obtained from all inversion results lie very closely to the actual values associated with $Z^*$ (e.g. $Z_M \approx Z^*$)). Therefore, we can further show that our final density structure models converge to the most optimal solution by repeatedly implementing our algorithm and analyzing its behavior even in the presence of noise in the observations.

---

4. Discussion

We implemented Primal-Dual Interior-Point methods in our optimization scheme for the inversion of Bouguer gravity anomaly datasets. We used this approach in conjunction with a 2-D polygonal prisms' forward model for the imaging of the Earth's shallow subsurface. Using both equality and inequality constraints helped reduce the solution space by advancing towards an optimal solution while moving within its interior feasible region. The use of a perturbation parameter, a merit function based on the barrier problem, a line search sufficient decrease condition for the merit function, and additional techniques within the optimization of our constrained minimization problem allows us to ensure that our sequence of approximate linear minimization problems will iteratively and efficiently converge to an optimal solution. Unlike conventional non-linear least squares inversion methods, using a constrained optimization scheme allows us to ensure that the final 2-D density distribution models conserve their integrity and meet the boundary and structural conditions associated to the area. The formulation of the primal-dual system allows us to work with a symmetrical system (9) that is considerably better conditioned with respect to the original non-symmetric and usually highly indefinite system (8). Furthermore, additional modifications can be performed to our system of equations to reduce its size and improve its performance ^[44].

We applied our MATLAB algorithm to three synthetic examples with varying body inter-dependence and complexity levels. The use of a priori information by adding physical bounds over our body parameter $Z$ coming from the known structure (both seismic and non-seismic datasets in real life applications) makes the inversion more stable and reliable, reduces the model space, and avoids spurious solutions that although geologically unfeasible still satisfy the imposed mathematical conditions (e.g. minimization of the objective function). We used relative and residual error measurements comparing the final and initial structural models for different examples to determine the effectiveness demonstrated by our algorithm in terms of the recovery of the actual density distributions. We tested the nonuniqueness inherent to this type of dataset by running our algorithm using different initial density structure models $Z_0$ for each optimization and using different levels of noise in the observations. Our algorithm allows us to have starting density structure models $Z_0$ with elements outside of the boundaries; those unfeasible values are moved to the feasible region during the first iteration.

In order to validate our results, we compared the final structural models $Z$ and their corresponding calculated gravity profiles obtained from using our constrained optimization approach with those coming from a conventional unconstrained formulation using non-linear least squares. Under a variety of conditions, the inversions using PDIP methods were found to be more robust when compared to their NLSQ counterpart mainly due to the substantially improved conditions obtained from the a priori equality and inequality model constraints. The physical bounds included in the formulation of this problem helped us to bound the variability of the model parameter, $Z$, which reduced the model space to avoid spurious solutions that may minimize the objective function while being physically unfeasible. In order to improve the convergence to an optimal model that further minimizes the objective function and residuals, additional analysis to find the "ideal" perturbation parameters, merit function, and line search sufficient decrease conditions should be performed.

All our synthetic representations suffered from the non-uniqueness inherently associated with gravity anomaly inversions which affected our final results. The higher inter-dependence within the bodies in some of the density structure models may have played a role in the final results. Note that the bounds used for the values of $Z$ covered a fairly large range and the effect of small changes in the physical characteristics of a body may greatly affect its gravity signature; therefore, tighter individualized constraints for each vertex may help us converge closer to the hypothetical density structure model, $Z^*$, rather than going to a non-optimal feasible solution (as was the case for the NLSQ algorithm). However, placing too many constraints on the parameter values would further reduce our feasible region making it harder to locate a feasible solution which also minimizes the objective function.

Throughout the paper we mentioned the advantages and disadvantages observed in the implementation of the PDIP methods for the inversion of Bouguer gravity anomalies, which are summarized in Table 5.

In this paper, we focused on the use of PDIP methods for the inversion of 2-D Bouguer gravity anomalies, however, the technique may also be used for the optimization of 3-D Bouguer gravity anomalies for the location and delineation of ore bodies and additional anthropogenic features in archaeological prospecting through the implementation of a rectangular prisms algorithm as a forward model. Furthermore, the general technique can be extended to a wide variety of other geophysical data sets for which adequate a priori knowledge exists by applying the optimization algorithm to the joint inversion of distinct complementary datasets ^[11,47]. The current limits for the applicability of this method to big problems lay mainly on the reliability of gravity datasets and their resolution. Given the non-uniqueness associated with this type of datasets, it is important to employ other geological or geophysical information in order to draw definite conclusions about the area under study. Due to their reliability, resolution, and stability, different seismic datasets are often used in conjunction with gravity. As stated previously, the existence and uniqueness of a solution, and the stability of the solution process are essential elements of inverse problems; hence, we have to guarantee that these conditions exist even when integrating disparate data sets in the joint inversion (relating both types of datasets through a common structure ^[12] or similar structural variations of different medium properties ^[48]).

Future work will focus on the inclusion of this methodology into a joint inversion scheme for the inversion of gravity anomalies and compatible seismic datasets (i.e. receiver functions and surface wave dispersion data ^[49,50,51]). Adding explicit constraints of the physical structures within the Earth into our inverse problem formulation will help us to enforce structural similarity ^{[11,48,52,53]} to improve the stability of our joint inversion scheme.

---

5. Conclusion

We present a constrained formulation for the inversion of Bouguer gravity anomalies solved using Primal-Dual Interior-Point methods. We show that a priori geological and/or geophysical information can be added into the objective function through the inclusion of explicit equality and inequality physical constraints. This approach helps to reduce ambiguities raised by the inherent non-uniqueness of gravity datasets. Our work consists in the application of well-known constrained optimization techniques to this popular geophysical inverse problem often solved using alternative methods. Given the ease of use of our MATLAB algorithm and the behavior of our inversion method, we believe that our approach provides a good alternative for the optimization of Bouguer gravity anomalies.

---

Acknowledgments

This research project is based on work supported by the National Science Foundation under grants #HRD--0734825, #HRD--1242122, and #DGE--0947992, and the Program in Computational Science at the University of Texas at El Paso. The authors would like to thank our anonymous reviewers for their valuable comments and excellent suggestions to improve our manuscript.

---

Conflict of Interest

All authors declare no conflicts of interest in this paper.

---