2.5D regularized inversion for the interpretation of residual gravity data by a dipping thin sheet: numerical examples and case studies with an insight on sensitivity and non-uniqueness

Abstract

A new two-and-a-half dimensional (2.5D) regularized inversion scheme has been developed for the interpretation of residual gravity data by a dipping thin-sheet model. This scheme solves for the characteristic inverse parameters (depth to top z, dip angle θ, extension in depth L, strike length 2 Y, and amplitude coefficient A) of a model in the space of logarithms of these parameters (log(z), log(θ), log(L), log(Y), and log(|A|)). The developed method has been successfully verified on synthetic examples without noise. The method is found stable and can estimate the inverse parameters of the buried target with acceptable accuracy when applied to data contaminated with various noise levels. However, some of the inverse parameters encountered some inaccuracy when the method was applied to synthetic data distorted by significant neighboring gravity effects/interferences. The validity of this method for practical applications has been successfully illustrated on two field examples with diverse geologic settings from mineral exploration. The estimated inverse parameters of the real data investigated are found to generally conform well with those yielded from drilling. The method is shown to be highly applicable for mineral prospecting and reconnaissance studies. It is capable of extracting the various characteristic inverse parameters that are of geologic and economic significance, and is of particular value in cases where the residual gravity data set is due to an isolated thin-sheet type buried target. The sensitivity analysis carried out on the Jacobian matrices of the field examples investigated here has shown that the parameter that can be determined with the superior accuracy is θ (as confirmed from drilling information). The parameters z, L, Y, and A can be estimated with acceptable accuracy, especially the parameters z and A. This inverse problem is non-unique. The non-uniqueness analysis and the tabulated inverse results presented here have shown that the parameters most affected by the non-uniqueness are L and Y. It has also been shown that the new scheme developed here is advantageous in terms of computational efficiency, stability and convergence than the existing gravity data inversion schemes that solve for the characteristic inverse parameters of a sheet/dike.

Background

Gravity and magnetic methods have many successful applications in mineral prospecting, environmental applications, and crustal imaging (e.g., Abdelrahman et al. 1989; Ateya and Takemoto 2002; Batista-Rodríguez et al. 2013; Beiki and Pedersen 2011; Fedi 2007; Hinze 1990; Hinze et al. 2013; LaFehr and Nabighian 2012; Long and Kaufmann 2013; Mehanee 2014a; Nettleton 1976; Okubo et al. 2013; Paoletti et al. 2013; Pei et al. 2014).

Techniques based on some geometrically simple bodies (e.g., sheet) have been in use for the analysis of magnetic and gravity data acquired along profiles since the 1960s (e.g., Grant and West 1965). These techniques (e.g., Ateya and Takemoto 2002; Mushayandebvu et al. 2001; Sazhina and Grushinsky 1971) are applied to obtain the characteristic inverse parameters of the underlying target.

Three-dimensional (3D) thin-sheet model is a very good approximation to a prismatic structure (e.g., an ore vein), unless the thickness is somewhat greater than the depth to the top (Telford and Geldart 1976) (Fig. 1). 3D thin-sheet models can be used to describe and resemble, for example, dikes or veins in exploration geophysics (Grant and West 1965; Pirajno and González-Álvarez 2013; Santosh and Pirajno 2014). Furthermore, thin-sheet solutions are favored for their inexpensive computational cost (Holstein et al. 2010) compared to full parallelepiped (e.g., Sazhina and Grushinsky 1971) and polyhedral (e.g., Holstein and Ketteridge 1996; Holstein 2003) solutions.

Fig. 1

A sketch showing geometry and characteristic parameters (depth to the top z, extension in depth L, finite extension in the strike direction 2 Y, thickness t, and direction and amount of dip θ) of a 3D dipping thin-sheet-like target (Grant and West 1965)

Inverse interpretation of an isolated residual gravity anomaly by a 3D thin-sheet model remains of interest in geophysical prospecting (e.g., Grant and West 1965; Telford et al. 1976). The main objective of the interpretation in this case is to retrieve the characteristic inverse parameters of the 3D approximative model: depth to the top, thickness, extension in depth, extension in the strike direction, and direction and amount of dip.

Grant and West (1965) described a graphical method for the interpretation of a residual gravity anomaly measured along a profile by a 3D dipping thin-sheet model. The approach was successfully applied to sulfide prospecting. However, the drawback with this method is that it is based on sets of main and auxiliary characteristic curves, and it essentially demands interpolation between pairs of characteristic curves in order to estimate the magnitudes of the sheet parameters. Telford and Geldart (1976) presented and described the forward modeling formula of 3D thin-sheet model, which is the basis for the inversion scheme developed here. They also highlighted the limitations (which are reported above) of this formula. Holstein et al. (2010) used the concept of gravi-magnetic similarity, and extended the thin-sheet potential modeling formula to include the potential, field, and field gradient in both gravity and magnetic cases when the buried bodies have uniform density or magnetization. Holstein and Anastasiades (2010) derived an exact finite expansion for the forward modeling computation of the gravitational anomaly of a uniform thin-polygonal sheet. This expansion exhibits the required absolute numerical stability (Holstein and Anastasiades 2010). One can refer to Holstein and Ketteridge (1996) and Holstein et al. (2014), and the references therein, for a more detailed and comprehensive information on this particular thin-polygonal sheet subject. Beiki and Pedersen (2011) developed a data window constrained two-dimensional (2D) inversion technique for the interpretation of gravity gradient tensor data using dike/sheet and contact models. By definition, a data window contains a number of particular data points out of the entire measured gravity profile (Beiki and Pedersen 2011). These particular data points are measured with respect to the center (that is located at the maximum response of the profile) of this window. This means that the gravity data points located to the left and right flanks of a window were not considered in the inversion. In other words, the full measured gravity signature (that is, primarily influenced by the buried anomalous structure) was not considered in the interpretation. A series of successive data windows is constructed by increasing the number of data points used (out of the whole gravity profile) for estimating the target parameters (Beiki and Pedersen 2011). Beiki and Pedersen (2011) used the MATLAB nonlinear least squares optimization tool “lsqnonlin” to solve the minimization problem of their scheme. This minimization tool is based on the Levenberg-Marquardt method (Levenberg 1944; Marquardt 1963). The convergence of the Levenberg-Marquardt method and its corresponding inverse parameters (e.g., the depth, dip angle, and thickness of the buried target) obtained from each window were plotted in their paper. According to Beiki and Pedersen (2011), solution with the smallest data fit error is selected (out of the entire solutions retrieved from all attempted data windows) as the most reliable one. It appears from the numerical models analyzed by Beiki and Pedersen (2011) that the use of a few data points (out of the whole measured gravity profile) in inversion was essential so that the Levenberg-Marquardt method can converge. Note that the aforementioned limitations (that are pertinent to the number of data points used in inversion and to the convergence and stability of the minimizer) are satisfactorily treated and fully well illustrated here as will be seen. Note that the Levenberg-Marquardt method (Madsen et al. 2004) was used also in the seismic tomography inversions.

Holstein et al. (2014) derived the gravity potential, field, and field gradient for triangular sheet targets with continuous densities across the boundaries to make it possible to approximate general curved surfaces. They verified the obtained modeling formula via numerical testing and simulated a model described by a dike curvature and a density compaction.

From this review, it appears to us that Grant and West (1965) are the only authors to have developed a method for the quantitative inverse interpretation of a residual gravity profile measured over a 3D dipping structure by 3D dipping thin-sheet model. Consequently, to the best of our knowledge, a regularized scheme (based on 3D dipping thin-sheet type model) for the inversion of residual gravity anomaly measured along a profile over a 3D body was not developed before.

The objective of this paper is to develop an efficient and rigorous regularized inversion scheme (based on 3D dipping thin-sheet model) for the interpretation of a residual gravity anomaly measured along a profile traversing the center of a 3D buried anomalous structure and normal to the strike of this structure. This inversion scheme is called two-and-a-half dimensional (2.5D) inversion (see e.g., Hinze et al. 2013, p. 177). The developed scheme is capable of dealing with the nonlinearity and the ill-posedness of this 2.5D inverse problem and has many merits. First, it inverts the entire residual gravity data set rather than just a few characteristic points out of this data set. Second, it simultaneously recovers all the characteristic inverse parameters (depth z, direction and amount of dip θ, extension in depth L, strike length 2 Y, and amplitude coefficient A) of the model (Fig. 1). Note that the parameter A is given by Δ ρ t (A=Δ ρ t), where Δ ρ and t are the density contrast and thickness of the sheet (a discussion on this parameter is given in the “Forward modeling solution” section), respectively. Third, the algorithm uses the exact forward modeling formula and evaluates the Jacobian (Fréchet) matrix (sensitivity matrix) of the scheme analytically and accurately. Fourth, the developed scheme (referred here as the logarithmic-space algorithm) employs the Tikhonov regularization (Tikhonov and Arsenin 1977) and the steepest descent (SD) and Gauss-Newton (GN) methods (e.g., Zhdanov 2002) in the space of logarithms of the model parameters (log(z), log(θ), log(L), log(Y), and log(|A|)) rather than in the space of the parameters themselves (z,θ,L,Y, and A). The use of the logarithmed parameters is essential in order to impose the positivity property of these parameters and hence to maintain the convergence and stability of the inversion scheme. The main achievements of this paper are as follows: (1) the use of the logarithm in all inverse parameters of the anomalous body in order to maintain the positivity of these parameters and the convergence of the inversion scheme, (2) the investigation of the non-uniqueness (that turned out to be imperative as will be seen and discussed here) of this particular inverse problem, and (3) the insights pertinent to the model parameters sensitivities in relation to non- uniqueness.

The paper is structured as follows. First, we briefly describe the direct problem (“Forward modeling solution” subsection). Second, the formulation of this particular inverse problem and its solving is discussed. Third, before applying this method to real data, the accuracy of the scheme is assessed and verified to numerical models with and without noise. Finally, the applicability of the developed technique to real data is demonstrated on two published field data examples from mineral exploration.

Methods

Forward modeling solution

The gravity anomaly (g), due to a 3D dipping thin-sheet target at a point x along a profile traversing the center of the target and normal to the target’s strike direction (Fig. 1), has a closed-form solution and is given (Grant and West 1965; Telford and Geldart 1976) by

$$ {\fontsize{7.2}{12}{\begin{aligned} &\qquad\qquad\qquad g\left(x_{i}, \Delta \rho, t, \theta, z, Y, L\right) = 2 \; \gamma \; \Delta \rho \; t \;\\ &\left[\frac{1}{2} \; \sin \theta \; \log \; \left\{\frac{\left\{ (z + L \sin \theta)^{2} + (x_{i} + L \cos \theta)^{2} +Y^{2} \right\}^{\frac{1}{2}} - Y}{\left\{ (z + L \sin \theta)^{2} + (x_{i} + L \cos \theta)^{2} +Y^{2} \right\}^{\frac{1}{2}} + Y}\right.\right.\\ &\qquad\qquad\qquad\qquad\left. \times \frac{\left({x_{i}^{2}} + z^{2} + Y^{2}\right)^{\frac{1}{2}}+ Y}{\left({x_{i}^{2}} + z^{2} + Y^{2}\right)^{\frac{1}{2}}- Y} \right\}\\ &\left.- \; \cos \theta \; \tan^{-1} \; \left\{\frac{Y\; (z \sin \theta + L + x_{i} \cos \theta) }{\left\{ (z + L \sin \theta)^{2} + (x_{i} + L \cos \theta)^{2} + Y^{2} \right\}^{\frac{1}{2}} (x_{i} \sin \theta - z \cos \theta)} \right\}\right.\\ &\qquad\qquad\left.+ \; \cos \theta \; \tan^{-1} \; \left\{\frac{Y\; (z \sin \theta + x_{i} \cos \theta) }{\left\{ {x_{i}^{2}} + z^{2} + Y^{2} \right\}^{\frac{1}{2}} (x_{i} \sin \theta - z \cos \theta)} \right\} \right], \end{aligned}}} $$

((1))

where x _i is the coordinate (m) of the measurement station, Δ ρ is the density contrast (kg.m^-3), t is the thickness (m), θ is the dip angle (degrees) of the body measured anticlockwise, z is the depth (m) to the top, 2 Y is the strike length (m), L is the dipping extent (m) of the body, and γ is the gravitational constant (6.67384 × 10^-11 m ³ kg ⁻¹ s^-2).

As indicated in the “Background” section, the thin sheet described above is a good approximation to a 3D prismatic structure, unless the thickness is somewhat greater than the depth to the top (z) (Telford and Geldart 1976). The presence of the thickness (t) outside the main square bracket of formula (1) is a consequence of the thin-sheet assumption used in deriving this formula.

Note that formula (1) can be used to accurately simulate a two-dimensional (2D) target by substituting in this formula a large value for the parameter Y (Grant and West 1965). Analysis of formula (1) with respect to the 2D simulation accuracy is beyond the scope of this paper.

It is noted that the parameters Δ ρ and t can be combined into a single parameter (the so-called amplitude coefficient, A (kg.m^-2) =Δ ρ t) to help minimize the non-uniqueness nature of this particular inverse problem. Non-uniqueness means that many different models (approximative solutions) could fit the observed data with the same accuracy (e.g., Tarantola 1987). The matter pertinent to the aforementioned parameter combination will be discussed further in the “Discussion” section.

Formulation of this 2.5D inverse problem and its solving

In this paper, we seek to solve the discrete 2.5D nonlinear inverse problem described by the operator equation (e.g. Menke 2012)

$$ G(m) = g_{\circ}, $$

((2))

where G is the forward modeling operator acting on m to produce some predicted gravity data (g(x, A, z, Y, L, θ)) at a finite number (N) of observation points along a profile, m is a column vector of the model parameters (that is the vector of some approximative solution) we seek to retrieve from inversion (m=[A,z,Y,L,θ]^T), g _∘ is a column vector of a finite set of noisy gravity data measured along this profile, $\phantom {\dot {i}\!}{g}_{\circ }=\left [g_{{\circ }_{1}}, g_{{\circ }_{2}}, g_{{\circ }_{3}}, \ldots g_{{\circ }_{N}}\right ]^{T}$, and T is the transposition operator.

Equivalently, Eq. 2 can be written as

$$ {g}\left(x_{i}, A, z, Y, L, \theta\right) = {g}_{\circ_{i}}, \; \; 1 \leq i \leq N. $$

((3))

Recent advances to solve the inverse problem of gravimetric data are based on the use of deterministic approaches that utilize the regularized least squares techniques and the differentiability of the objective function subject to minimization (e.g., Zhdanov 2002). As indicated earlier, the inversion scheme described in this paper recovers the characteristic parameters of the buried target in a minimal time (a few seconds). Therefore, automatic deterministic approaches are much more efficient (see, for example, Mehanee 2015) than the approaches that are based on the trial-and-error modeling method.

Since the residual gravity data are usually corrupted with some noise, we do not seek to fit them exactly (Tarantola 2005). Rather, we seek (Ramlau 2005; Zhdanov 2002)

$$ \left\Vert {g}-{g}_{\circ}\right\Vert = \left [ \sum_{i=1}^{N} ({g}(x_{i}, A, z, Y, L, \theta) - {g}_{\circ_{i}})^{2} \right ]^{\frac{1}{2}} \; \leq \,\delta, $$

((4))

where g is a column vector of the predicted data calculated along the profile, using formula (1), from some approximative solution m, g=[g ₁,g ₂,g ₃,…g _N ]^T, and δ is the noise level embedded in the measured data (g _o ).

Solution of inverse problem (2) is non-unique (that is, different sets of approximative solution can equally fit the data subject to inversion) as it will be seen from the examples carefully analyzed here. Therefore, it is ill-posed (e.g., Gribok et al. 2002). The conventional method of solving ill-posed inverse problems is based on the minimization of the Tikhonov parametric functional (e.g., Mehanee and Zhdanov 2002; Tikhonov et al. 1998):

$$ \begin{aligned} &\qquad\qquad\qquad \phi \, (m, g_{\circ}, \alpha) \; = \; \left\Vert \; G(m) \; - \; g_{\circ} \; \right\Vert^{2} \; \; +\\ &\; \alpha \; \left\{ \; C_{z} \left\Vert \; z \; \right\Vert^{2} \; + \; C_{\theta} \left\Vert \; \theta \; \right\Vert^{2} \; + \; C_{L} \left\Vert \; L \; \right\Vert^{2} \; + \; C_{A} \left\Vert \; A \; \right\Vert^{2} \; \right.\\&\qquad\qquad\qquad\left.+ \; C_{Y} \left\Vert \; Y \; \right\Vert^{2} \; \right\} = min, \end{aligned} $$

((5))

where α (dimensionless) is the regularization parameter, the selection of which is addressed in the “Discussion” section, and $C_{z} \left (\frac {{\text {mGal}}^{\text {2}}}{{\mathrm {m}}^{\text {2}}}\right), C_{\theta } \left (\frac {\text {mGal}^{\text {2}}}{\text {degree}^{\text {2}}}\right), C_{L} \left (\frac {{\text {mGal}}^{\text {2}}}{{\mathrm {m}}^{\text {2}}}\right), C_{A} \left (\frac {{\text {mGal}}^{\text {2}} \mathrm {m}^{\text {4}}}{{\text {kg}}^{\text {2}}}\right)$, and $C_{Y} \left (\frac {{\text {mGal}}^{{\text {2}}}}{{\mathrm {m}}^{\text {2}}}\right)$ are some constants employed in (5) in order to adjust the inconsistency of units of the six terms of ϕ (m,g _∘,α) and in order to make the values of the C terms comparable (balanced). Given the aforementioned units of those constants, ϕ (m,g _∘,α) in this case (that is, in this particular choice of units) is given in the units of mGal². Equivalently, (5) can be readily re-written as:

$$ \phi \, (m, g_{\circ}, \alpha) \; = \; \left\Vert \; G(m) \; - \; g_{\circ} \; \right\Vert^{2} \; \; + \; \alpha \; \; \left\Vert \; C \; m \; \right\Vert^{2} \; \; = min, $$

((6))

where $C = \text {diag}\left [\sqrt {C_{z}}, \; \sqrt {C_{\theta }}, \; \sqrt {C_{L}}, \; \sqrt {C_{A}}, \; \sqrt {C_{Y}}\right ]$.

The first term of (6) is the data misfit functional, determined as the square norm of the difference between the observed and predicted data, and the second term is a stabilizing functional, the stabilizer.

The characteristic model parameters (A,z,Y,L, and θ) of a 3D thin sheet can have large variations in magnitude (e.g., z=100 m, A=10⁵ kg.m^-2) which can mean that the smaller parameters are ignored in the inversion as they have smaller sensitivities. Performing the inversion in the logarithmic space of the model parameters (log(|A|), log(z), log(Y), log(L), and log(θ)) makes these logarithmed parameters and their corresponding sensitivities comparable and balanced as will be seen. In other words, the logarithmed variants will not dominate one another as do the non-logarithmed values. Furthermore, the inversion in the logarithmic space of the model parameters forces and guarantees the positivity of the model parameters we seek. Inversion in the non-logarithmic (linear) space (that is, the space of the model parameters themselves: A, z, Y, L, and θ) can generate negative model parameters, and hence results in scheme divergence (e.g., Mehanee 2014a). Therefore, the new logarithmed-space objective functional takes the form:

$$ \Psi \, (\widetilde {m}, g_{\circ}, \alpha) = \tau \; \left\Vert \; \widetilde G(\widetilde {m}) \; - \; g_{\circ} \: \right\Vert^{2} \; + \; \alpha \; \left\Vert \; \widetilde {m} \; \right\Vert^{2} = min, $$

((7))

where τ is a scaling parameter used in order to make Ψ dimensionless and is set to 1 $\left (\tau = 1 \; \frac {1}{{\text {mGal}}^{\text {2}}}\right)$, and $\widetilde {m} = [\log (|A|/A_{\circ }), \log (z/z_{\circ }), \log (Y/Y_{\circ }), \log (L/L_{\circ }), \log (\theta /\theta _{\circ }) ]^{T}$. Each of the parameters A _∘, z _∘, Y _∘, L _∘, and θ _∘ (the so-called scaling parameters) is of a unit (positive) value and is introduced in order to make the quantities that are under the logs in $\widetilde {m}$ dimensionless. The introduction of these parameters is physically necessary because the logs do not allow dimensions to be defined. As mentioned above, α is a dimensionless quantity. Note that hereinafter the scaling parameters of $\widetilde {m}$ are dropped to simplify the notations.

All the numerical models and real data examples shown here are inverted and analyzed by the logarithmic minimization formulation presented in (7). We note that, in the restricted case of the noise-free and neighboring/interference effect-free numerical example (Model 1), the nonlinear minimization problem (7) is solved iteratively using a sequential hybrid technique that automatically combines the SD and GN methods (see Mehanee 2015, and the references therein) in order to verify and validate the developed inversion scheme prior to applying it to noisy and real data sets.

The use of the hybrid technique was essential. This is because the SD method converges very slowly or stagnates (as found from extensive noise-free numerical experiments) at the final stage of its iterative minimization process. It was also found from these experiments that the GN method essentially requires a very good initial guess in order to converge (Zhdanov 2002). That is why the SD method (which does not require a good initial guess to converge as will be seen in the “Numerical tests” subsection) is employed first in the developed scheme in order to generate and prepare a suitable initial guess for the GN method.

In the framework of this hybrid scheme, the SD method is used first until a normalized misfit (defined as $\frac {\left \Vert \widetilde G(\widetilde {m}) \; - \; {g}_{\circ }\right \Vert }{\left \Vert {g}_{\circ }\right \Vert } \times 100~\%$) of roughly around 5−10 % (depending on the model subject to inversion) is reached. After that, the GN method is applied to the observed data, where the inverse results produced from the SD method is used as the initial guess for the GN method. The GN method terminates when a normalized misfit below 10^-5 % is reached.

We note that when inverting the distorted (by noise and/or by interference due to nearby bodies) data sets (Model 2 and Model 3) shown in the “Numerical tests” subsection, and when inverting the two real data sets presented in the “Field data inversion” subsection, the minimization problem (7) was solely solved by the SD method. Unlike the noise-free data set (Model 1), we found in our experience for this particular inverse problem that the SD method can be sufficient for the inversion of the distorted and real field data examples. This is because in these examples (the distorted and real field data), we neither expect nor seek to exactly fit the observed data (see the figures pertinent to Model 2 and Model 3).

The entire computational steps of the GN and SD methods of the inverse scheme are presented and discussed in Appendices 1 and 2. The code of the scheme is implemented in MATLAB. We note that no built-in MATLAB optimization (minimization) functions were used in the 2.5D code of our paper.

Sensitivity (Fréchet) matrix calculation

The Jacobian (Fréchet) matrix ($\widetilde {{F}}$) used in this scheme is calculated analytically with respect to the inverse parameters we seek (log(|A|), log(z), log(Y), log(L), and log(θ)) as:

$$ \begin{aligned} \frac{\partial g}{\partial \log(|A|)} &= \frac{\partial g}{\partial A} \; \; \frac{\partial A}{\partial \log(|A|)} = A \; \; \frac{\partial g}{\partial A},\\ \frac{\partial g}{\partial \log(z)} &= \frac{\partial g}{\partial z} \; \; \frac{\partial z}{\partial \log(z)} = z \; \; \frac{\partial g}{\partial z},\\ \frac{\partial g}{\partial \log(Y)} &= \frac{\partial g}{\partial Y} \; \; \frac{\partial Y}{\partial \log(Y)} = Y \; \; \frac{\partial g}{\partial Y},\\ \frac{\partial g}{\partial \log(L)} &= \frac{\partial g}{\partial L} \; \; \frac{\partial L}{\partial \log(L)} = L \; \; \frac{\partial g}{\partial L}, \; \text{and}\\ \frac{\partial g}{\partial \log(\theta)} &= \frac{\partial g}{\partial \theta} \; \; \frac{\partial \theta}{\partial \log(\theta)} = \theta \; \; \frac{\partial g}{\partial \theta}. \end{aligned} $$

((8))

The use of the logarithmed model parameters has the benefit of making the aforementioned derivatives dimensionally the same, which is a good strategy for balancing the Jacobian terms. The quantities $\frac {\partial g}{\partial A}$, $\frac {\partial g}{\partial L}$, $\frac {\partial g}{\partial z}$, $\frac {\partial g}{\partial \theta }$, and $\frac {\partial g}{\partial Y}$ are evaluated analytically by differentiating the forward modeling formula (1).

Singular value decomposition

The singular value decomposition (SVD) of an m ×n matrix (in our case this matrix is $\widetilde {{F}} \; \epsilon \; R^{N \times 5}$, where N and 5, respectively, are the number of data points along the gravity profile and model parameters) is given (e.g., Jupp and Vozoff 1975) by:

$$ \widetilde {{F}}=U \, S \, V^{T}, $$

((9))

where U=[u ₁,u ₂,…,u _m ] ε R ^m×m and V= [v ₁,v ₂,…,v _n ] ε R ^n×n are both orthogonal matrices. S ε R ^m×n is a diagonal matrix with non-negative real values (the so-called singular values) arranged in decreasing order, s ₁≥s ₂……≥s _n ≥0. The sequence S=[s ₁,s ₂,…,s _n ] is referred to as the singular spectrum of $\widetilde {{F}}$. The columns of U=[u ₁,u ₂,…,u _m ] and V=[v ₁,v ₂,…,v _n ] are the left and right singular vectors in the input and output spaces of the transformation represented by $\widetilde {{F}}$, respectively. The magnitude of the singular values in S represents and determines the corresponding important singular vectors in the columns of U and V. The SVD is a useful tool for understanding the sensitivity analysis of the various model parameters of the 3D thin-sheet model as will be seen.

Results and discussion

Numerical tests

Prior to applying the developed inversion algorithm to real data examples, its accuracy is assessed and analyzed first on three numerical models. All inverse solutions are rounded to the first integer.

Model 1: noise-free example

The gravity response of a thin-sheet type model with z = 25 m, L = 50 m, Y = 500 m, θ=30^o, and A = 5700 kg/m ² is generated at the ground surface (Fig. 2).

Fig. 2

Model 1: noise-free data. Top panel: preliminary (initial) inversion results obtained from the steepest descent (SD) method. Bottom panel: final inversion results obtained from the Gauss-Newton (GN) method which used the inversion results produced from the SD (top panel) as initial guess; see the text for details. Bottom panel shows that the observed and predicted data are coincident with each other

The presented algorithm supplemented with the SD method was applied to this noise-free dataset for 6000 iterations (yielding to a normalized misfit of about 4.7 %) using the initial guess shown at the top panel of Fig. 2. The corresponding inversion results (the so-called the preliminary inverse solution) are shown at the top panel of Fig. 2. To avoid possible stagnation with the SD method, the same algorithm supplemented with the GN method was applied to the same residual gravity data. The aforementioned preliminary inverse solution evolved from the SD method was employed automatically, as initial guess, in the GN method. The bottom panel of Fig. 2 shows the final inverse solution obtained from the scheme.

It is seen that the algorithm has converged below 10⁻⁶ % (normalized misfit) and successfully recovered the true value of all inverse parameters of the buried body. Note that the observed data and predicted response are coincident with each other. The regularization parameter (α) used in the SD and GN methods is 10⁻⁵ and 10⁻¹², respectively. Figure 3 depicts the behavior of the misfit and objective functional corresponding to the SD and GN inversion results shown in Fig. 2.

Fig. 3

Model 1: noise-free data. Behavior of the misfit and objective functional corresponding to the SD (top panels) and GN (bottom panels) inversion results shown in Fig. 2

Model 2: noisy example

In order to assess and analyze the robustness and stability of the inverse algorithm, the developed scheme has been applied to a noisy profile. Subsequent noise levels of about 7, 11, and 20 % have corrupted the gravity data generated by a numerical model described by the following: z = 12 m, L = 35 m, Y = 100 m, θ=120°, and A = 12,000 kg/m ². Each noise level was separately generated and added to the noise-free data set $\widetilde G(\widetilde {m})$ using the MATLAB function “awgn” in order to produce the corresponding corrupted data set g _∘ subject to inversion. The aforementioned noise levels were separately calculated as $\frac {\left \Vert {g}_{\circ } \; - \; \widetilde G(\widetilde {m}) \right \Vert }{\left \Vert {g}_{\circ }\right \Vert } \times 100\,\%$.

The corresponding final inversion results produced by the SD method are rendered in Fig. 4. The scheme terminated when the embedded noise level is reached by the SD method. Figure 5 shows the behavior of the misfit and objective functional corresponding to the inversion results shown in Fig. 4. An α of 10⁻⁶ was used in all the inverse calculations of this model.

Fig. 4

Model 2: noisy data. Inversion results corresponding to various noise levels. Data corrupted with about 7 % (top panel), 11 % (middle panel), and 20 % (bottom panel) noise. The results shown in this figure and all subsequent figures are recovered by the SD method solely. See the text for details

Fig. 5

Model 2: noisy data. Behavior of the misfit and objective functional corresponding to the inversion results shown in Fig. 4. The top, middle and bottom panels correspond to those of Fig. 4

The obtained inversion results show that the developed inverse technique is stable with respect to the noise levels embedded. Tables 1 and 2 present the inverse solutions and the corresponding error (in percentage) of each model parameter measured with respect to the true value of this parameter. One can see that error range of about 4–39 % was encountered in the retrieved inverse solutions. It can be seen that the highest errors are associated with the data set contaminated with the highest noise level (20 %), which is not unexpected. Note that the model parameter θ encountered the smallest error in the case of the 20 % noise. Whereas the rest of the model parameters encountered the smallest error in the case of the 11 % noise. This could be attributed to the variation in the model parameter sensitivities and noise level that contaminated the data.

Table 1 Model 2 (noisy data). Inversion results of various noise levels. The true model parameters are z = 12 m, L = 35 m, Y = 100 m, A = 12,000 kg/m ², and θ= 120 °

Table 2 Model 2 (noisy data). Error of the inversion results shown in Table 1

Model 3: interference effect example

In some geologic settings, the measured gravity signature of a 3D body may sometimes be influenced by the neighboring/interference effect (that is, the effect of some nearby anomalous bodies) (e.g., LaFehr and Nabighian 2012). In order to examine and assess this effect on the accuracy of the inverse solution, we have computed, using formula (1), the forward modeling response of each of three spatially distributed bodies (Fig. 6, bottom panel). These bodies could resemble some intrusive dikes, which are simulated by thin sheets. The model parameters and forward modeling response of each body are illustrated at the top panel of Fig. 6. The profile of the composite gravity response of all three bodies is shown in the middle panel. This profile shows three anomalous signatures: the prominent anomaly of interest here (for example, the ore body) present in the middle of the profile and the two minor anomalies surrounding the prominent one.

Fig. 6

Model 3: neighboring/interference effect of three nearby bodies. Top panel: individual forward modeling responses computed using formula (1). Middle panel: composite gravity response subject to inversion (see Fig. 7). Bottom panel: a sketch showing the causative bodies whose characteristic parameters are shown in the top panel. See the text for details

Fig. 7

Model 3: neighboring/interference effect of three nearby bodies. Top panel: approximative inverse solution obtained from the inversion of the whole profile shown in the middle panel of Fig. 6. Bottom panel: approximative inverse solution obtained from the inversion of the extracted (truncated) response (marked by arrows in the middle panel of Fig. 6). See the text for details

Two interpretive scenarios are explored and investigated for analyzing and inverting the profile of the composite effect (Fig. 6, middle panel) by a single thin sheet (the subject of this paper). Scenario 1 (whole profile inversion) in which the whole profile is inverted and analyzed. Scenario 2 (extracted (truncated) profile inversion) in which only the distorted prominent anomaly (marked by arrows in Fig. 6, middle panel) out of the whole profile is inverted. Scenario 2, in this particular case, is probably more realistic (in terms of inverting its data by a single thin-sheet model) as its gravity response can resemble a single thin-sheet model. Scenario 1 is not so realistic (in terms of inverting its data by a single thin-sheet model) as its gravity response subject to inversion has three anomalous signatures that cannot be attributed to or represented by a single thin-sheet model. However, the inversion of scenario 1 has been carried out here solely for the sake of investigative purpose, better understanding and for the full illustration of the developed inversion scheme.

The top panel of Fig. 7 shows the true model parameters of the main body and the corresponding approximative inverse results of Scenario 1 for which a normalized misfit of about 10.85 % is reached using an α of 10⁻⁶. The inversion results of Scenario 2 are shown in the bottom panel of Fig. 7. A normalized misfit of 6.71 % is obtained, and an α of 10⁻⁹ was used in the inverse computations. This analysis shows that the parameters A, L, and Y of the main body exhibited significant inaccuracy.

On the basis of this investigation, we can conclude that the developed scheme produces more accurate inversion results for an isolated target, and thus it is more suitable for isolated anomalies or for anomalies with insignificant neighboring/interference effect. Isolated targets (such as ore veins and dikes) are frequent in many geologic settings and can be of a paramount economic interest as will be seen in the next two real data examples in which the gravity data were acquired for chromite, copper, and gold prospecting.

Field data inversion

In order to examine and assess the applicability of the developed inversion technique, two published field examples, from mineral exploration in Canada and Cuba, with various depth of burial, geologic complexity, and interference effects are analyzed. These particular data sets were selected in this research paper for a number of reasons. First, these data sets were generated by buried causative ore bodies, which can resemble thin-sheet models (Davis et al. 1957; Grant and West 1965; Siegel et al. 1957), the subject of this paper. Second, these data sets were measured from sites with known drilling information; hence, we can compare the numerical results yielded from the inversion against those confirmed from drilling. Third, these sites have core samples from which the density contrasts were accurately estimated in laboratory. As pointed out earlier, the knowledge of the density contrast (Δ ρ) is essential in order to compute the thickness (t) of the body from the amplitude coefficient (A) obtained from the inversion.

The Mobrun anomaly, Canada

Inverse results A prominent Bouguer gravity anomaly was observed over a massive sulfide ore body in the Noranda Mining District, Quebec, Canada (Grant and West 1965; Siegel et al. 1957). Figure 8 shows the SW-NE residual gravity profile (whose location is shown in Figure 10-1, page 272, Grant and West 1965) taken normal to the strike of the buried causative ore vein. Directional drilling intersected sulfides over a distance slightly greater than 30.5 m (Grant and West 1965, page 282, see also Figure 10-11 in page 281).

Fig. 8

The Mobrun anomaly, Canada. Inversion results obtained from two different regularization parameters (α) and same initial guess. The results shown at the top and bottom panels correspond to a normalized misfit of 2.93 and 2.48 %, respectively. The measured and predicted data are shown in solid and open circles, respectively

Grant and West (1965) reported that the average density of core samples of the mineral taken from drilling was 4600 ± 500 kg/m ³, and for the host rock, the density was about 2700 kg/m ³. Hence, the ore body has a density contrast (Δ ρ) of 1400–2400 kg/m ³. Thus, the aforementioned two density values are equal to the bulk densities of the ore and the host rock, respectively.

The inverse algorithm has been applied to the above mentioned residual gravity profile. Figure 8 shows the two approximative inverse solution sets obtained from two different regularization parameters (α) using an initial guess of (z = 2 m, L = 30 m, A = 2200 kg/m ², Y = 50 m, and θ=30°). The inversion results shown at the top panel were recovered using an α of 10⁻⁷ and correspond to a normalized misfit of 2.93 %. The bottom panel illustrates the results retrieved from the algorithm using an α of 10⁻⁴, for which a normalized misfit of 2.48 % was reached. It is noted that the density contrast range mentioned above was used to estimate the corresponding thickness (t) range of the ore body from the amplitude coefficient (A=t Δ ρ) evolved from the inversion.

In order to monitor and assess the non-uniqueness (equivalence) issue, which is not unexpected for this particular inverse problem, several initial guesses were also attempted in the scheme. Table 3 summarizes all obtained inversion results. It can be seen from the table that the regularization parameter indeed can affect the inverse solutions.

Table 3 The Mobrun anomaly, Canada. Tabulated inversion results. α is the regularization parameter

One can see that the two approximative solution sets (the so-called here “Equivalent Model 1” and “Equivalent Model 2”) shown at rows 3 and 5 of Table 3 are different, though each of them has the same normalized misfit value (2.7 %). This is attributed to the non-uniqueness nature of this particular inverse problem; that is, different inverse models can equally fit the observed data.

On the basis of the a priori information (which includes the contour form of the residual gravity map from which the gravity profile subject to inversion was taken), it can be suggested that the approximative inverse solutions presented in rows 4, 5, and 6 of Table 3 are not very realistic. On the other hand, these tabulated results suggest that the most common model (interpretive model to resemble the buried vein deposit) could be that shown in rows 1, 2, and 3. In light of this, it is reasonable to suggest that the corresponding inverse parameters of the interpretive model can be the averages of those shown in the aforementioned three rows. Therefore, the interpretive model is described by z = 12.4 m, L = 184.4 m, A = 55,791.2 kg/m ², Y = 118.6 m, θ=86°, and t = 23–40 m. Figure 9 shows the Mobrun sulfide ore body, delineated from drilling, at 50 m depth (from Grant and West 1965), and the interpretive model suggested here. The results confirmed from drilling (z, t, and Y), and those recovered from inversion are in reasonable match. We note that other gravity inversion methods (e.g., Li and Oldenburg 1998; Zhdanov et al. 2004) could produce a better fit.

Fig. 9

The Mobrun anomaly, Canada. Plan showing the Mobrun sulfide ore body (in gray), outlined from drilling, at 50 m depth (from Grant and West 1965). The inversion results are marked by the black box

The parameter A retrieved from inversion (A = 55,791.2 kg/m ²) suggests a thickness of 23–40 m (t = 23–40 m) which is greater than the depth recovered from inversion (z=12.4 m). Probably this is because the actual ore body is not a thin sheet. Supporting evidence is that the actual ore body has a thickness of 31 m as revealed from a slice (constructed from drilling) taken at 50 m depth (Fig. 9). And that the depth to the top of the target is less than 30.5 m as inferred from directional drilling that intersected the body at about 30.5 m. As pointed out in the “Forward modeling solution” subsection, the 3D thin sheet is a good approximation to a prismatic structure, unless the thickness (t) is somewhat greater than the depth to the top (z) (Telford et al. 1976). However, we sought in this paper to fit and interpret this residual gravity profile by a thin-sheet-like model for investigative purposes and the full illustration of the developed inverse scheme, and in order to gain some insights.

Prior to any drilling in this area, Siegel et al. (1957) analyzed this gravity profile based on a 2D trial and error modeling method and reported a corresponding solution of z=6 m, L=186 m, and t=31 m. It is relevant to note that the strike length of the deposit in this case is assumed to be infinite as the interpretation was essentially based on 2D assumption. Grant and West (1965) too interpreted this profile but as a 3D thin sheet, based on a graphical method, and the retrieved inverse solution is as follows: z=17 m, L=170 m, 2 Y (strike length) = 205 m, θ=83°, and t=36 m. These results, too, show that the approximative model is quasi-vertical. Using an Euler deconvolution technique, Roy et al. (2000) interpreted this profile by a 2D vertical ribbon model and obtained an inverse model of z = 22.7 m and L = 52 m. While their results are in some variations from the true results confirmed from drilling, these results are still useful as they could be used as a reasonable initial guess in deterministic (gradient-type) inversions. As noted above, other gravity inversion methods (e.g., Li and Oldenburg 1998; Zhdanov et al. 2004) could produce a more accurate inverse solution and fitting.

Sensitivity analysis

In order to get some insights on a possible mutual interrelation between the non-uniqueness of this particular inverse problem and the model parameter sensitivities (that is, the Fréchet matrix; $\widetilde {{F}} = \frac {\partial g}{\partial \widetilde {m}}$), the sensitivities of the aforementioned two equivalent solutions (“Equivalent Model 1” and “Equivalent Model 2”) and of the two initial guesses used in the inversion of these two equivalent solutions have been calculated.

The top panel of Fig. 10 depicts the sensitivities (absolute values) computed from the initial guess (shown at row 3 of Table 3) from which “Equivalent Model 1” was evolved. The sensitivities (absolute values, too) of the inverse parameters of “Equivalent Model 1” (shown, too, at row 3 of Table 3) are illustrated in the bottom panel of Fig. 10. The corresponding sensitivities of “Equivalent Model 2” and its initial guess (both are shown at row 5 of Table 3) are rendered in Fig. 11.

Fig. 10

The Mobrun anomaly, Canada. Sensitivity pertinent to this profile’s “Equivalent Model 1” and its initial guess (both are shown at row 3 of Table 3), see the text for details. Top panel: sensitivities calculated from the initial guess. Bottom panel: sensitivities calculated from “Equivalent Model 1”

Fig. 11

The Mobrun anomaly, Canada. Sensitivity pertinent to this profile’s “Equivalent Model 2” and its initial guess (both are shown at row 5 of Table 3), see the text for details. Top panel: sensitivities calculated from the initial guess. Bottom panel: sensitivities calculated from “Equivalent Model 2”

Note that the sensitivities of the inverse parameters log(L) and log(Y) (the two main parameters that appear to be most affected by the non-uniqueness of this particular inverse problem as will be strengthened further and seen in the sensitivity subsection of the field example of the Camaguey area, Cuba) of “Equivalent Model 1” and “Equivalent Model 2” are found to be positive real numbers. The bottom panel (corresponding to these two equivalent models) of Figs. 10 and 11 qualitatively reveals that the sensitivity curves of these two parameters (log(L) and log(Y)) have a similar behavior and form. The degree of conformity of these two sensitivity curves depends upon the values of the model parameters used in the sensitivity calculation. The top and bottom panels of Figs. 10 and 11 show that the parameter log(θ) has the dominant sensitivity. Pertinent quantitative analysis based on the singular value decomposition (e.g., Press et al. 1988) is provided in the next paragraph.

In order to understand better the mutual interrelation between the non-uniqueness of this particular inverse problem and the model parameters sensitivity, the SVD (e.g., Jupp and Vozoff 1975, Press et al. 1988) was carried out separately on the Fréchet matrix of “Equivalent Model 1” and of the initial guess of this equivalent model.

Table 4 shows the singular values and their variabilities (in percentage) of the sensitivity matrix shown at the top panel of Fig. 10. The model parameters used in the calculations of this matrix are also listed at this panel. It can be seen from Table 4 that the singular values encountered a very rapid decay, and that most of the variability (97.73 %) is captured by the first singular vector in the parameter space (right singular vectors). Note that the first three singular vectors in the parameter space collectively achieved a variability of more than 99.9 % (Table 4). Figure 12 depicts the various components of the three aforementioned singular vectors. The left top panel (in which the components corresponding to the first singular value (s ₁=19.7) are plotted in linear scale) shows that the parameter log(θ) has the most dominant effect. The right top panel (in which the absolute value of the same aforementioned components are plotted in log scale) shows that the order of the importance of parameters (for this particular model parameter set) is as follows: log(θ), log(|A|), log(L), log(Y), and log(z). This panel shows that the components pertinent to the parameters log(L) and log(Y) nearly have the same magnitude, this agrees well with the sensitivity curves (of these two parameters) rendered at the top panel of Fig. 10 as the first right singular vector has the most variability as indicated above (97.73 %). Furthermore, the middle and bottom panels (pertinent to the second and third right singular vectors, which collectively achieved a variability of 2.18 %) of Fig. 12 also show that the components pertinent to the parameters log(L) and log(Y) have a comparable magnitude.

Fig. 12

The Mobrun anomaly, Canada. Parameter eigenvectors corresponding to the first three singular values (denoted by s) (see the text for details) calculated from the singular value decomposition (SVD) analysis carried out on the Fr$\acute {\mathrm {e}}$chet matrix of the initial guess shown at the top panel of Fig. 10. Left panels: parameter eigenvectors in linear scale. Right panels: parameter eigenvectors’ absolute values in log scale

Table 4 The Mobrun anomaly, Canada. Singular values and their variabilities in percentage of the sensitivity matrix, which was calculated from the initial guess shown at the top panel of Fig. 10

Table 5 shows the singular values and their variabilities of the sensitivity matrix (shown at the bottom panel of Fig. 10) of “Equivalent Model 1”. This table indicates that most of the variability (96.7 %) is also ubiquitous to the first singular vector in the parameter space. Figure 13 depicts the various components of the first three right singular vectors. The top left panel shows that the parameter log(θ) has the most dominant effect. This agrees with that revealed from the top left panel of Fig. 12. The top right panel of Fig. 13 shows that the order of the importance of parameters this time (that is, for “Equivalent Model 1”) is as follows: log(θ), log(L), log(Y), log(|A|), and log(z). This panel shows that the components pertinent to the parameters log(L), log(Y), and log(|A|) have a comparable magnitude, which conforms with the sensitivity curves (of these three parameters) shown at the bottom panel of Fig. 10. Furthermore, the middle and bottom panels of Fig. 13 also reveal that the components pertinent to the parameters log(L) and log(Y) have nearly a similar magnitude.

Fig. 13

The Mobrun anomaly, Canada. Parameter eigenvectors corresponding to the first three singular values (see the text for details) calculated from the SVD of the Fr$\acute {\mathrm {e}}$chet matrix of “Equivalent Model 1” presented at the bottom panel of Fig. 10. Left panels: parameter eigenvectors in linear scale. Right panels: parameter eigenvectors’ absolute values in log scale

Table 5 The Mobrun anomaly, Canada. Singular values and their variabilities in percentage of the sensitivity matrix, which was obtained from “Equivalent Model 1” shown at the bottom panel of Fig. 10

The Camaguey Province anomaly, Cuba

Inverse results

Davis et al. (1957) carried out detailed gravity surveys in the Camaguey Province, Cuba, for chromite exploration. Residual gravity maps were collected over various spatially distributed ore bodies for determining approximately the quantity of chromite in any deposits found. In this paper, a profile from this area is inverted and interpreted. This profile is taken normal to the strike of the residual gravity anomaly shown in the middle part of Fig. 6 of Davis et al. (1957).

This gravity profile is associated with the largest chromite deposit found in the province (Davis et al. 1957). This prominent profile has a trend of S60 ^oW to N60 ^oE and overlies a chromite ore body which contains about 115,000 tons, dips steeply to the southwest, and comes within 3 m of the ground surface (Davis et al. 1957).

We have run the inversion using a number of various initial guesses to see the most common inverse solution. The top, middle, and bottom panels (which have identical observed data and point curves) of Fig. 14 show the three approximative solution sets yielded from three different initial guesses. A normalized misfit of about 7.56 % was reached in all cases. Set 1 has a solution of z = 5.5 m, L = 40 m, A = 6989 kg/m ², 2Y (strike length) = 60 m, θ=93°, and t = 4.7 m using an α of 10⁻⁴. Set 2 has a solution of z = 6.2 m, L = 29.5 m, A = 7647 kg/m ², 2Y= 82 m, θ=94°, and t = 5 m using an α of 10⁻⁶. Set 3 has a solution of z = 6.5 m, L = 29.8 m, A = 7979 kg/m ², 2Y=72 m, θ=94°, and t = 5.3 m using an α of 10⁻⁶. It is seen that the approximative solutions of the parametersL and Y of Set 1 and Set 2 are nearly swapped. This is readily attributed to the non-uniqueness nature (solution equivalence) of this inverse problem—these two solutions (Set 1 and Set 2) are so-called here “Equivalent Model 1” and “Equivalent Model 2”. It appears that the common solution is the one that is shown in Set 2 and Set 3. The depth shown in all sets is in a reasonable agreement with that obtained from drilling. Table 6 summarizes all inversion results.

Fig. 14

The Camaguey Province anomaly, Cuba. Residual gravity anomaly (redrawn from Figure 6 of Davis et al. 1957): inversion results retrieved from three various initial guesses

Table 6 The Camaguey Province anomaly, Cuba. Tabulated inversion results

The top panel of Fig. 15 illustrates the sensitivities (absolute values) calculated from the initial guess (shown at the top panel of Fig. 14) from which “Equivalent Model 1” was produced. The sensitivities (absolute values, too) of the inverse parameters of “Equivalent Model 1” (shown too at the top panel of Fig. 14) are illustrated in the bottom panel of Fig. 15. The corresponding sensitivities of “Equivalent Model 2”, and its initial guess (both are shown at the middle panel of Fig. 14) are rendered in Fig. 16.

Fig. 15

The Camaguey Province anomaly, Cuba. Sensitivity pertinent to this profile’s “Equivalent Model 1” and its initial guess (both are shown at the top panel of Fig. 14), see the text for more details. Top panel: sensitivities calculated from the initial guess. Bottom panel: sensitivities calculated from “Equivalent Model 1”

Fig. 16

The Camaguey Province anomaly, Cuba. Sensitivity pertinent to this profile’s “Equivalent Model 2” and its initial guess (both are shown at the middle panel of Fig. 14), see the text for more details. Top panel: sensitivities calculated from the initial guess. Bottom panel: sensitivities calculated from “Equivalent Model 2”

The sensitivity of the inverse parameters log(L), log(Y), and log(|A|) of “Equivalent Model 1” and “Equivalent Model 2” are found positive real number for this profile. The bottom panels of Figs. 15 and 16 show that these three parameters have quasi-similar behavior and form—this is consistent with the pertinent findings reported for the Mobrun anomaly, Canada (see the bottom panels of Figs. 10 and 11).

Sensitivity analysis

Table 7 shows the singular values and their variabilities (in percentage) of the sensitivity matrix shown at the top panel of Fig. 15. This table shows that most of the variability (79.1 %) is captured by the first right singular vector. Figure 17 illustrates the various components of the first three singular vectors. The top left panel shows that the parameter log(θ) has the most dominant effect as was observed, too, in “The Mobrun anomaly, Canada” subsection. The top right panel shows that the order of the importance of parameters (for this particular model parameter set) is as follows: log(θ), log(|A|), log(z), log(L), and log(Y). This panel as well as the middle and bottom panels show that the components of all parameters have diverse magnitudes.

Fig. 17

The Camaguey Province anomaly, Cuba. Parameter eigenvectors corresponding to the first three singular values (denoted by s) (see the text for details) calculated from the SVD analysis carried out on the Fréchet matrix of the initial guess shown at the top panel of Fig. 15. Left panels: parameter eigenvectors in linear scale. Right panels: parameter eigenvectors’ absolute values in log scale

Table 7 The Camaguey Province anomaly, Cuba. Singular values and their variabilities in percentage of the sensitivity matrix, which was calculated from the initial guess shown at the top panel of Fig. 15

Table 8 shows the singular values and their variabilities of the sensitivity matrix shown at the bottom panel of Fig. 15. This table shows that most of the variability (96.9 %) is held by the first right singular vector. Figure 18 illustrates the various components of the first three singular vectors. The top left panel shows that the parameter log(θ) has the most dominant effect as was observed earlier. The top right panel shows that the order of the importance of parameters (for this particular model parameter set) is as follows: log(θ), log(|A|), log(L), log(Y), and log(z). This panel and the middle and bottom panels show that the components of all parameters have comparable magnitudes.

Fig. 18

The Camaguey Province anomaly, Cuba. Parameter eigenvectors corresponding to the first three singular values (see the text for details) calculated from the SVD of the Fr$\acute {\mathrm {e}}$chet matrix of “Equivalent Model 1” presented at the bottom panel of Fig. 15. Left panels: parameter eigenvectors in linear scale. Right panels: parameter eigenvectors’ absolute values in log scale

Table 8 The Camaguey Province anomaly, Cuba. Singular values and their variabilities in percentage of the sensitivity matrix, which was calculated from “Equivalent Model 1” shown at the bottom panel of Fig. 15

The SVD analysis presented herein revealed that the parameter log(θ) has the highest importance in this inverse scheme. All field examples, carefully analyzed here, have shown that the parameter θ is determined very accurately (based on drilling information), and that this parameter is the least affected by the initial guess selection and the misfit stopping criteria. It is worthy to note that accurate determination of the amount and direction of dip (θ) of the buried target can be of a paramount importance for effective decision making on directional drilling.

On the basis of the inversion results of the two investigated field examples, it can be concluded that the model parameter which can be estimated from the developed method with the greatest accuracy is θ. The parameters z, A, L, and Y can be determined with an acceptable accuracy, especially when the gravity response is due to an isolated body.

Discussion

Gravity inversions based on regular models remain of interest in exploration geophysics (see for example, Biswas 2015, and the references therein). It is worthy to note that it is very rare to find a geologic target that is truly regular body. Nevertheless, these regular models often yield a first approximation solution, which is sometimes quite adequate (e.g., Abdelrahman et al. 1989; Mehanee 2014a).

As indicated in the “Background” section, although a geologic simplification, 3D thin-sheet model (the subject of this paper) is a good approximation to a prismatic structure (e.g., a dike or an ore vein) unless the thickness is somewhat greater than the depth to the top (Telford et al. 1976). Thin-sheet inverse schemes are favored for their quickness and inexpensive computational cost (Holstein et al. 2010) but have much reduced applicability compared to 3D rigorous inversions (e.g., Li and Oldenburg 1998; Zhdanov et al. 2004). The scheme developed in our paper is useful in exploration geophysics and applicable:

• For reconnaissance analysis and pilot studies (for first approximative interpretation) in geophysical prospecting prior to conducting large-scale gravity data surveys,

• For comparative and ambiguity studies and joint (e.g., gravity and magnetic data) inversion in order to minimize the non-uniqueness that is imperative here, and

• For integrated and unequivocal interpretation. For example, this scheme can determine the amount and direction of dip of the buried target with greatest accuracy (as seen in the investigated inversions of the noisy numerical data and the field examples), which was confirmed by the sensitivity analysis carried out in our paper.

This scheme is complimentary not contradictory to the existing interpretation methods. For example, the approximative inverse solution evolved from the thin-sheet inverse scheme can be used to build a reasonable initial guess for the rigorous 3D gravity inverse schemes. Note that a reasonable initial guess is a prerequisite for the rigorous 3D gravity inverse schemes. It is emphasized that we are not claiming that inversion schemes based on geometrically simple bodies replace the rigorous 3D gravity inverse schemes.

Rigorous 3D gravity inversion schemes (e.g., Li and Oldenburg 1998; Zhdanov et al. 2004) can simultaneously invert for several anomalous irregular bodies and normally require areal data coverage. However, these schemes take much longer computation time (than the thin-sheet schemes) in order to yield the approximative inverse solution (that is also non-unique), which comprises the 3D anomalous density distribution (Δ ρ) in the subsurface. As indicated above, a reasonable initial guess is a prerequisite in the rigorous 3D gravity inverse schemes. These rigorous 3D schemes can generate smooth or focused (compact) inverse images (depending primarily upon the particular objective of the interpretation) by incorporating the appropriate type of the stabilizer in the objective functional subject to minimization. It is relevant to note that the inversion of focused (compact) images requires accurate information about the lower and upper bounds of the anomalous density distribution (Δ ρ) of the buried targets.

Upon the completion of the gravity data acquisition and processing for an area, the following interpretation steps (e.g., Grant and West 1965; Hinze et al. 2013) are applied in order to estimate the characteristic parameters of a buried target:

First, contour the measured residual gravity data (regardless whether these data were acquired in a disordered or an ordered pattern). Second, examine the obtained contour map, and identify the strike of the buried target (Hinze et al. 2013). Third, take a profile traversing the center of the target and normal to the target’s strike direction (e.g., Hinze et al. 2013). Fourth, under the thin-sheet assumption, invert this profile by the 2.5D scheme developed here using several initial guesses. This will give first approximate interpretation about the underlying body. The value of the parameter θ (which determines the amount and direction of dip of the buried body) will be the most accurate (see column 12 of Tables 3 and 6) as confirmed from the numerical models, filed examples, and the sensitivity analysis. The obtained inverse solutions should be interpreted in an integrated manner with all available geological and geophysical information.

Grouping the density contrast (Δ ρ) and thickness (t) of a 3D dipping thin-sheet body into a single inverse parameter A (A=Δ ρ t) was essentially carried out because both of them appear only as a multiplicative combination. This grouping was done in order to help minimize the non-uniqueness nature of this particular inverse problem. Unfortunately, the rest of the inverse parameters (z,L,Y, and θ) we seek to recover from inversion are present frequently at various terms within formula (1) and are multiplicative (that is, they are coupled). Therefore, grouping any of these parameters (z,L,Y, and θ) is not really helpful in eliminating or mitigating the non-uniqueness of this inverse problem.

The particular 2.5D inverse problem developed here is non-unique as has been seen, thus it is ill-posed. A way of solving an ill-posed problem is via the use of regularization (Gribok et al. 2002; Tikhonov and Arsenin 1977). Regularization is required here for a number of reasons: First, in order to help the minimizer combat the entrapment in a local minimum. This can be accomplished by attempting various values for the regularization parameter in the scheme when inverting a gravity data set using an initial guess (see, for example, Fig. 8 and Table 3). Second, in order to find possible equivalent approximative solutions to understand better the non-uniqueness nature of this particular inverse problem, and thus to interpret these solutions in an integrated manner with the available geological, geophysical, and drilling data. Third, in order to make it possible to incorporate some a priori information, if available, in the stabilizer of the objective functional subject to minimization, which can help minimize the non-uniqueness nature of this inverse problem. The regularization parameter α is chosen such that it makes some balance between the misfit functional term and stabilizer (see, e.g., Li and Oldenburg 2003; Mehanee 2014b, and the references therein).

The inherent non-uniqueness of this inverse problem and the errors in the residual gravity data lead to the existence of equivalent solutions. If approximate prior knowledge on one of the target’s parameters (e.g., extension in depth L) is known, then we can incorporate this a priori information in the stabilizer of the objective functional. This could reduce the equivalence problem. It is not mere coincidence that the results obtained here for the two investigated field data examples have been found in some conformity with those confirmed by drilling. Therefore, the inverse scheme developed here can have great potential in exploration and mining geophysics for an isolated target described by a 3D thin-sheet type model.

Conclusions

We have developed an efficient regularized iterative 2.5D inversion scheme for the interpretation of a residual gravity profile measured over a dipping thin-sheet like target. The scheme determines the characteristic parameters (depth to top z, amount and direction of dip θ, extension in depth L, finite strike extension 2 Y, and amplitude coefficient A from which the thickness t is obtained) of a buried target. The algorithm of the scheme solves for the inverse parameters of a model in the space of their logarithms. This is advantageous for a number of reasons: First, in order to maintain the convergence of the objective functional subject to minimization. Second, in order to impose the positivity of the model parameters we seek, and hence, to produce realistic and meaningful inversion results. Third, in order to balance the Jacobian terms and to make the sensitivity derivatives dimensionally the same. It has been shown that this new scheme is advantageous in terms of computational efficiency, stability, and convergence than the existing schemes that solve for the characteristic inverse parameters of a sheet/dike model from gravity data inversions.

Before applying the method to real data examples, it has been successfully verified on noise-free numerical examples; it has recovered the actual model parameters. After that, the approach was assessed on noisy numerical data, and it is found stable and can estimate the parameters of the buried deposit with acceptable accuracy. However, some of the inverse parameters encountered some inaccuracy when the method was applied to residual data distorted, in terms of both magnitude and shape, by some significant neighboring gravity effects generated by nearby anomalous bodies.

The validity of the technique for practical applications has been successfully illustrated on two real field examples with various geologic settings and complexities from mineral exploration. The estimated inverse parameters of the investigated real data are found to generally conform well with those yielded from drilling.

The sensitivity analysis carried out on the Jacobian matrices of the field examples investigated here has shown that the parameter that can be determined with the greatest accuracy is θ (as confirmed from drilling information). Accurate determination of the amount and direction of dip (θ) of the buried target can be of a paramount importance for effective decision-making on directional drilling. The parameters z, L, Y, and A can be estimated with acceptable accuracy, especially the parameters z and A. Real data inversions have shown that the parameter θ is the least affected by the choice of the initial guess and the misfit stopping criteria.

The developed inversion scheme is useful in exploration programs and reconnaissance studies intended to delineate and map, for example, ore veins from residual gravity data. However, it can produce a non-unique inverse solution; a fact which should be kept in mind when interpreting the obtained approximative inverse solutions. Therefore, these equivalent solutions, as always, should be guided by geological information and other available geophysical results to help resolve any encountered non-uniqueness, which is not unusual in exploration geophysics. The non-uniqueness analysis and the tabulated inverse results presented here have shown that the parameters that are most affected by the non-uniqueness, the choice of the initial guess, and the misfit value of the stopping criteria are L and Y.

Simultaneous inversion of 3D dipping thin-sheet-like multi-bodies will be the subject of future research.

The gravity signature due to a buried target depends upon the density contrast among the other characteristic parameters (e.g., extension in depth and extension in the strike direction) of the buried target. The density contrast between the ore body and the country rock can be in some cases much smaller than the magnetic susceptibility contrast. Therefore, the buried target in this case may generate a more prominent magnetic signature than the gravimetric one. Thus, the magnetic method, in this case, could be of a more value than the gravity method. Therefore, the 2.5D inversion scheme developed here could also be extended and useful in the interpretation of magnetic data.

Appendix

Appendix 1 Gauss-Newton (GN) method

The computational steps of the GN method are summarized as (e.g., Menke 2012; Zhdanov 2002):

$$ \begin{aligned} {R}_{n} \, &= \, \widetilde G(\widetilde {{m}}_{n}) \, - \, {g}_{\circ} \, = \, G ({m}_{n}) \, - \, {g}_{\circ}\\ \widetilde {{l}}_{n} \, &= \, \widetilde {{F}}^{T}_{n} \; {R}_{n} \, + \, \alpha \; \widetilde {{m}}_{n}\\ \widetilde {{H}}_{n} \, &= \, (\widetilde {{F}}^{T}_{n} \; \widetilde{{F}}_{n} \, + \, \alpha \; {I})\\ \widetilde {\delta {m}}_{n} \, &= \, \widetilde {{H}}_{n}^{-1} \, \; \widetilde {{l}}_{n}\\ \widetilde {{m}}_{n+1} \, &= \, \widetilde {{m}}_{n} \, - \, \widetilde{\delta {m}}_{n}, \end{aligned} $$

((10))

where $\widetilde {{ m}}_{n}$ is the column vector of the model parameters at iteration n; $\widetilde {{ m}}_{n} = [\log (|A_{n}|)$, log(z _n ), log(Y _n ), log(L _n ), log(θ _n )]^T; T is the transpose operator; R _n is the column vector of the difference between the predicted (G(m _n )) and observed (g _∘) gravity data sets at iteration n; $\widetilde {{F}}_{n}$ is the Fréchet (Jacobian) matrix (Menke 2012; Tarantola 1987, 2005) computed at iteration n with respect to the log of the model parameters; α is the regularization parameter; $\widetilde {\delta {m}_{n}}$ is the model parameters update at iteration n; $\widetilde {{{l}}}_{n}$ is the direction of the steepest ascent computed at iteration n; and I and H are the identity and Hessian matrices.

Appendix 2: Steepest descent (SD) method

The computational steps of the SD method (Menke 2012; Zhdanov 2002) are summarized as:

$$ \begin{aligned} {R}_{n} \, &= \, \widetilde G(\widetilde {{m}}_{n}) \, - \, {g}_{\circ} \, = \, G ({m}_{n}) \, - \, {g}_{\circ}\\ \widetilde{{l}}_{n} \, &= \, \widetilde {{F}}^{T}_{n} \; {R}_{n} \, + \, \alpha \; \widetilde {{m}}_{n} \\ \widetilde {{m}}_{n+1} \, &= \, \widetilde{{ m}}_{n} \, + \, \widetilde {\delta {m}_{n}} \, = \, \widetilde{{ m}}_{n} \, - \, \zeta_{n} \; \widetilde{{l}}_{n}, \end{aligned} $$

((11))

where ζ _n is the step length which is given by:

$$ \zeta_{n} = \frac{\left\Vert \, \widetilde {{l}}_{n} \, \right\Vert^{2}}{\left\Vert \widetilde {{F}}_{n} \; \widetilde {{l}}_{n} \right\Vert^{2} + \alpha \; \left\Vert \, \widetilde {{l}}_{n} \, \right\Vert^{2}}. $$

((12))