Journal of Applied Geophysics 121 (2015) 54–62

Contents lists available at ScienceDirect

Journal of Applied Geophysics

journal homepage: www.elsevier.com/locate/jappgeo

Interpretation of gravity data using 2-D continuous wavelet

transformation and 3-D inverse modeling
Amin Roshandel Kahoo a,⁎, Ali Nejati Kalateh a, Farshad Salajegheh b
a
School of Mining, Petroleum & Geophysics Engineering, University of Shahrood, Iran
b
CCFS-GEMOC, Department of Earth and Planetary Sciences, Macquarie University, Australia

a r t i c l e i n f o a b s t r a c t

Article history: Recently the continuous wavelet transform has been proposed for interpretation of potential field anomalies. In
Received 31 December 2014 this paper, we introduced a 2D wavelet based method that uses a new mother wavelet for determination of the
Received in revised form 17 June 2015 location and the depth to the top and base of gravity anomaly. The new wavelet is the first horizontal derivatives
Accepted 4 July 2015
of gravity anomaly of a buried cube with unit dimensions. The effectiveness of the proposed method is compared
Available online 16 July 2015
with Li and Oldenburg inversion algorithm and is demonstrated with synthetics and real gravity data. The real
Keywords:
gravity data is taken over the Mobrun massive sulfide ore body in Noranda, Quebec, Canada. The obtained results
2D continuous wavelet transform of the 2D wavelet based algorithm and Li and Oldenburg inversion on the Mobrun ore body had desired similar-
Li and Oldenburg inversion algorithm ities to the drill-hole depth information. In all of the inversion algorithms the model non-uniqueness is the chal-
Mobrun massive sulfide ore body lenging problem. Proposed method is based on a simple theory and there is no model non-uniqueness on it.
Gravity data © 2015 Elsevier B.V. All rights reserved.
Depth estimation

1. Introduction densities for inverting the data of the iron ore province of Quadrilátero
Ferrífero, southeastern Brazil.
Gravity surveys have been widely used in the hydrocarbon and min- Wavelet transform (Mallat, 2008) is a powerful tool for non-
eral explorations. Investigation of the earth's heterogeneities is being in- stationary geophysical data analysis. Fedi and Quarta (1998) applied
creasingly referred by various methods such as inversion and wavelet the discrete wavelet transform on potential field data to separate the re-
transform of gravity data (Paterson and Reeves, 1985). gional and residual. Ridsdill-Smith and Dentith (1999) have developed
The inversion is the most essential step in the quantitative interpre- the application of wavelet transform for processing of aeromagnetic
tation of potential field data. This step of interpretation of potential field data. Sailhac et al. (2000) used the wavelet transform modulus maxima
data suffers from the non-unique determination of the source parame- (WTMM) method to estimate the source depth of magnetic data. The
ters (Fedi and Rapolla, 1999). Some researchers proposed the density discrete version of wavelet transform can be used to denoise aeromag-
variation and tried to invert the geometrical parameters (Guspi, 1992; netic and gravity gradiometry data (de Oliveira Lyrio et al., 2004;
Pedersen, 1977). Another group of researchers used more qualitative Leblanc and Morris, 2001). The estimation of source depth can be per-
prior information for inversion of potential field data. Last and Kubik formed by ratio of the first- and second order coefficients of complex
(1983) minimized the total volume of the source of anomaly and wavelet transform (Valee et al., 2004). Cooper (2006) introduced an-
Guillen and Menichetti (1984) used the moment of inertia of causative other CWT based method for interpretation of gravity and magnetic
body. Barbosa and Silva (1994) introduced the method based on gener- data acquired along a profile. His method could estimate the depth
alized compact gravity inversion technique to allow compactness along and position of potential source anomaly using continuous wavelet
the several axis using Tikhonov's regularization (Zhdanov, 2002) meth- transform. He used the first and second horizontal derivatives of buried
od. Li and Oldenburg (2003) used the logarithmic barrier method to in- cylinder and vertical derivative of them as a mother wavelet in his
verse the large-scale magnetic data for recovering the 3D distribution of study.
magnetic susceptibility. Most recently, another approach for the solu- We propose a new method for interpretation of gravity data which
tion of large-scale 3D inversion has been used under the name of “mov- was acquired on a regular grid. This method uses the 2D continuous
ing footprint” (Cox et al., 2010). Uieda and Barbosa (2012) have wavelet transform for position and depth estimation of gravity field
developed a new robust iterative algorithm by planting anomalous sources. This method is an extension of Cooper (2006) method. We
use the first horizontal derivatives of the anomaly of a buried cube
with unit dimensions as the mother wavelet. We illustrate our method
⁎ Corresponding author. with synthetics and real data and compare the results by 3D inversion of
E-mail address: [email protected] (A. Roshandel Kahoo). gravity data using Li and Oldenburg algorithm (Li and Oldenburg, 2003).

http://dx.doi.org/10.1016/j.jappgeo.2015.07.008
0926-9851/© 2015 Elsevier B.V. All rights reserved.
 A. Roshandel Kahoo et al. / Journal of Applied Geophysics 121 (2015) 54–62 55

Fig. 1. (a) Unit buried cube and (b) its gravity response for 10 m grid spacing in X and Y directions.

Fig. 2. First horizontal derivative of gravity anomaly over the buried unit cube in (a) X direction and (b) Y direction.

The synthetic data are contaminated by different levels of Gaussian 2. Methodology

noise. The real gravity data were taken over the Mobrun sulfide body
in Noranda, Quebec, Canada (Grant and West, 1965). 2.1. 2D continuous wavelet transform (2DCWT)

The one dimensional continuous wavelet transform of a continu-

ous signal f(x) at a scale s and translational value u is expressed by
Eq. (1):

Z þ∞ x−u
1
W f ðu; sÞ ¼ f ðxÞ pffiffi ψ dx ð1Þ
−∞ s s

where ψ(x) is a continuous function in spatial domain which is

called the mother wavelet and “*” denotes the conjugate complex
(Mallat, 2008). Thus, the wavelet transform is calculated as the
inner product of f(x) and translated and scaled versions of a mother
wavelet ψ(x).

Table 1
Depth and lateral position of three 3D prism.

Body no. xmin xmax ymin ymax zmin zmax

1 825 875 825 875 50 100

2 800 900 800 900 100 300
3 825 875 825 875 300 350
Fig. 3. 3D view of the model geometry of synthetic example.
56 A. Roshandel Kahoo et al. / Journal of Applied Geophysics 121 (2015) 54–62

Fig. 4. (a) Noise free synthetic gravity anomaly for model geometry shown in Fig. 3 and noise contaminated synthetic anomaly with (b) 5%, (c) 10% and (d) 15% Gaussian noise. The color
scale bar shows the gravity field in mgal.

In two dimensional cases, the continuous wavelet transform of a 2D x1 ≤ x ≤ x2, y1 ≤ y ≤ y2 and z1 ≤ z ≤ z2 has a vertical attraction at the origin
continuous signal f(x, y) is defined as Eq. (2): given by Eq. (3).
Z þ∞ Z þ∞  
1 x−a y−b X 2 X 2 X    
W f ða; b; s; θÞ ¼ f ðx; yÞψθ ; dxdy ð2Þ 2
xi y j  
jsj −∞ −∞ s s g ¼ γρ μ i jk zk arctan −xi log Ri jk þ y j −y j log Ri jk þ xi
i¼1 j¼1 k¼1
z R
k i jk

where, ψ(x, y) is a 2D mother wavelet, a and b are the spatial translation ð3Þ
parameters, s is the dilation parameter and θ is the orientation parame- qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ter of the mother wavelet which set to zero in our proposed method (Tai where, Ri jk ¼ x2i þ y2j þ z2k and μ ijk = (−1)i(−1) j(−1)k.
Sing, 1996). When g is used as the basis for a family of wavelets, γ and ρ are
The anomaly of a buried unit cube can be useful in applying 2D con- not needed and depth of unit cube is set to one. Because the
tinuous wavelet transform on potential field data, due to any source 2DCWT resolution is directly related to the size of the cube, it
body that can be represented as a linear combination of unit cubes should not be chosen bigger than grid spacing. For smaller quanti-
based on grid spacing. Fig. 1 shows the unit cube and its gravity re- ties, we obtained a small increase in resolution with higher compu-
sponse for 10 m grid spacing in X and Y directions which were calculat- tational cost. As can be seen, the anomaly of a buried unit cube does
ed based on Plouff (1976). Blakely (1996) shows that a rectangle cube not have a zero mean therefore, it does not meet the requirements
with uniform density ρ and with dimensions described by the limits of a mother wavelet (Mallat, 2008). So, we used its first horizontal
 A. Roshandel Kahoo et al. / Journal of Applied Geophysics 121 (2015) 54–62 57

Fig. 5. Three slices in different directions of total 2D wavelet coefficient of horizontal derivatives of (a) noise free synthetic data, (b) 5% noise contaminated data, (c) 10% noise contam-
inated data and (d) 15% noise contaminated data.

derivatives of unit cube gravity response as a mother wavelet in X 2.2. Li–Oldenburg inversion algorithm
and Y directions which are shown in Fig. 2. The first horizontal
derivatives in X and Y directions were calculated by numerical gra- The forward solution of gravity data can be defined as the follow-
dient. It is important to know that the 2D continuous wavelet ing matrix equation by dividing the source region into the 3D prism
transform analysis using the abovementioned mother wavelets is (Plouff, 1976).
applied to the first horizontal derivatives of gravity anomaly in X
and Y directions. Gρ ¼ d ð4Þ
58 A. Roshandel Kahoo et al. / Journal of Applied Geophysics 121 (2015) 54–62

Fig. 6. Illustration of the effect of the density bound constraint. The panels on the left are models obtained without positivity for (a) 5%, (c) 10% and (e) 15% levels of noise. The panels on the
right are the models obtained with positivity imposed for (a) 5%, (c) 10% and (e) 15% levels of noise.

where ρ is the vector of model parameter densities, d is the vector of point, and (z − zi) is the vertical distance between jth cell and observa-
observed gravity data and matrix G has as elements Gij that quantify the tion point (Li and Oldenburg, 1998).
contribution to the ith datum of a unit density in the jth cell as Eq. (5). Gravity inverse problem is formulated as an optimization problem
Z where an objective function of the model and data is minimized. The ob-
z−zi jective function of the Li–Oldenburg algorithm is defined as Eq. (6).
Gi j ¼ γ dv ð5Þ
ΔV j jr−r i j3
X
M  
where ΔVj is the cuboidal volume within jth cell, γ is the gravitational ϕðλÞ ¼ ϕd þ μϕρ −2λ ln ρ j ð6Þ
constant, |r − ri| is the distance between jth cell and observation j¼1
 A. Roshandel Kahoo et al. / Journal of Applied Geophysics 121 (2015) 54–62 59

where λ is the barrier parameter, μ is the regularization parameter, M is

the number of model parameter, ϕd and ϕρ are the data misfit and the
model objective functions, respectively (Li and Oldenburg, 2003). ϕd
and ϕρ can be computed as Eqs. (7) and (8).

T 2
ϕd ¼ ðGρ−dobs Þ W d ðGρ−dobs Þ ¼ kW d ðGρ−dobs Þk ð7Þ

2
ϕm ðρÞ ¼ ðρ−ρ0 ÞT W ρ ðρ−ρ0 Þ ¼ W ρ ðρ−ρ0 Þ ð8Þ

where, dobs is the field data vector, Wd = diag{1/σ1, …, 1/σN}, σi is the

standard deviation error related with the ith data point, ρ0 is the refer-
ence model and Wρ is the total weight of model parameters presented
by Li and Oldenburg (2003).

3. Synthetic examples

We now apply the 2DCWT and Li–Oldenburg inversion algorithms

to a synthetic example. The model consists of three 3D prismatic bodies
with a density contrast equal to 1.0 g/cm3. Fig. 3 shows a 3D view of the
model geometry which is summarized in Table 1. We have calculated
gravity anomaly data above the plain surface on 29 by 29 of 25 m spac-
ing. The noise free gravity anomaly of synthetic model is shown in
Fig. 7. Gravity anomaly overlayed on 2D view of the model geometry which is summarized Fig. 4a. Fig. 4b–d shows the noise-contaminated data at different levels
in Table 2. The color scale bar shows the gravity field in mgal. of signal to noise ratio.
We apply the 2DCWT on both noise free and noisy synthetic gravity
data with two mother wavelets shown in Fig. 2. The first horizontal de-
rivatives of the data were calculated in both X and Y directions. Then, we
apply the 2DCWT on the first horizontal derivatives of the synthetic data
with mother wavelets in the same direction. In Fig. 5, we show the ab-
Table 2 solute value of 2D wavelet coefficient in X and Y directions based on
Geometrical properties of two 3D prism.
Eq. (9).
Body Width Length Height Top center Rotation Density
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
no. (m) (m) (m) position (x,y,z) contrast g/cm3
TWCHD ¼ 2DWT 2THDX þ 2DWT 2THDY ð9Þ
1 100 100 100 (300,700,50) 45 1
2 140 140 150 (700,300,100) 45 −1
where TWCHD is the total 2D wavelet coefficient of horizontal deriva-
tives, 2DWTTHDX is the 2D wavelet coefficient first horizontal derivative

Fig. 8. The result of applying the 2DCWT on the model with sides that were not aligned EW and NS directions with different density contrasts. (a) A vertical contour section of the TWCHD
overlaid by isosurface of the wavelet coefficient value for a cutoff of 0.9 along the PP′ profile. (b) A horizontal contour slice of the TWCHD overlaid by isosurface of the wavelet coefficient
values for a cutoff of 0.9 at 120 (m) depth (c) and (d) are the XZ and YZ views of oblique slice, respectively.
60 A. Roshandel Kahoo et al. / Journal of Applied Geophysics 121 (2015) 54–62

of data in X direction based on mother wavelet shown in Fig. 2a and

2DWTTHDY is the 2D wavelet coefficient first horizontal derivative of
data in Y direction based on mother wavelet shown in Fig. 2b.
In Fig. 5a–d, we show the three slices in different directions of total
2D wavelet coefficient of horizontal derivatives of synthetic data with
different levels of noise. Overlain on all the 2DCWT slices is the synthetic
model geometry. The results show an appropriate correlation between
the synthetic model location and the maxima of the 2DCWT. Based on
the obtained results for the noisy data, we find out that the proposed
method have stability results against the noise.
To perform the inversion using Li–Oldenburg inversion algorithm,
we chose two different strategies to evaluate the effect of density con-
strains for synthetic model recovering. In the first strategy, we consider
a wide density model bound between −3.0 g/cm3 and 3.0 g/cm3. In the
second strategy, we chose a positivity density constraint with 1.0 g/cm3
as maximum density contrast of model. It's useful when a reliable esti-
mate of maximum density contrast is available.
Fig. 6 shows the inverted model of synthetic data with different
levels of noise by two mentioned strategies. Superimposed on the inver-
sion results are the true model (shown as red rectangle blocks). Fig. 6a, c
and e is the results of inversion algorithm by choosing the broad bound
density. They are characterized by broad tail at depth and wide density
distribution around the true model. Fig. 6b, d and f shows the inversion
results when positivity constraint is imposed. By considering the posi-
tivity constraint, the density distribution has been concentrated at the
true model location.
In the next synthetic example, we used a model with sides which
were not aligned EW and NS directions. To make the model more real-
istic, two separated bodies with different density contrasts were used
for evaluation of the model alignment and density contrast on the re-
Fig. 9. (a) Bouguer anomaly field of a Mobrun ore body near Noranda, Quebec. (b) Bouguer
gravity data over the AB profile.
sults of 2DCWT method. Fig. 7 shows the gravity anomaly overlayed

Fig. 10. The result of applying the 2DCWT on the Mobrun gravity anomaly. (a) A contour section of the total 2D wavelet coefficient overlaid by isosurface of the wavelet coefficient value for
a cutoff of 2.7. (b) Three cross sections from the total 2D wavelet coefficient at different directions. (c) Oblique cross sections from the total 2D wavelet coefficient along the AB profile.
 A. Roshandel Kahoo et al. / Journal of Applied Geophysics 121 (2015) 54–62 61

on 2D view of the model geometry which is summarized in Table 2. The

gravity anomaly has calculated above the plain surface on 40 by 40 of
25 m spacing.
Fig. 8a and b shows two vertical and horizontal slices from total 2D
wavelet coefficient of horizontal derivatives of synthetic data, respec-
tively. Overlain on the slice is the isosurface of the wavelet coefficient
value for a cutoff of 0.9. The XZ and YZ views of oblique slice are
shown in Fig. 8b and c. The results show an appropriate correlation be-
tween the position of the synthetic model and the maxima of the
TWCHD in the model with sides that were not aligned in EW and NS di-
rections with different density contrasts.

4. Real example

Using the new proposed algorithm, we analyzed field gravity data

(Fig. 9a) taken over a sulfide deposit, near Noranda, Quebec, which led
to the discovery of Mobrun ore body. Fig. 9b shows the Bouguer gravity
data on AB profile. This is a massive body of base metal sulfide (mainly
Fig. 11. 3D regularized massive sulfide density distribution obtained by the Li–Oldenburg pyretic) which has displaced volcanic rocks of middle Precambrian age
inversion algorithm. (Grant and West, 1965). The dataset consists of a regular grid of 60 × 57

Fig. 12. Two sections in two different directions from the 3D obtained massive sulfide density distribution.

Fig. 13. Center section of the Mobrun sulfide body (after Grant and West, 1965) with geophysical interpretation.
62 A. Roshandel Kahoo et al. / Journal of Applied Geophysics 121 (2015) 54–62

data in x and y directions respectively, 10 m spaced and digitized from stability against the random noise. The obtained results of 2DCWT anal-
Fig. 10.1 in Grant and West (1965). The average density of core samples ysis and 3D inversion by Li and Oldenburg algorithm on the Mubrun
of the mineral was about 4.6 g/cm3 and for the host rocks the density massive sulfide ore body had favorable similarity to the drill-hole
was about 2.7 g/cm3. Thus the density contrast of the sulfide body was depth information.
about 1.9 g/cm3 (Grant and West, 1965).
The result of applying the 2DCWT on the Mobrun gravity anomaly References
is shown in Fig. 10. In Fig. 10a–b, we show the cross sections of the
Barbosa, V.C.F., Silva, J.B.C., 1994. Generalized compact gravity inversion. Geophysics 59
total 2D wavelet coefficient. Overlain in Fig. 10a is the isosurface of (1), 57–68.
massive sulfide body for a cutoff of 2.7 of the wavelet coefficient Blakely, R.J., 1996. Potential Theory in Gravity and Magnetic Applications. Cambridge Uni-
value. The obtained results demonstrate the depth-to-top of the versity Press.
Cooper, G.R.J., 2006. Interpreting potential field data using continuous wavelet transforms
ore body is less than 25 m and the depth-to-base of the ore body is of their horizontal derivatives. Comput. Geosci. 32 (7), 984–992.
more than 220 m. Fig. 10c shows the oblique slice along the AB pro- Cox, L.H., Wilson, G.A., Zhdanov, M.S., 2010. 3D inversion of airborne electromagnetic data
file from TWCHD in which a black bold line indicates the Bouguer using a moving footprint. Explor. Geophys. 41, 250–259.
de Oliveira Lyrio, J.C.S., Tenorio, L., Li, Y., 2004. Efficient automatic denoising of gravity
gravity anomaly. gradiometry data. Geophysics 69 (3), 772–782.
The 3D regularized massive sulfide density distribution obtained by Fedi, M., Quarta, T., 1998. Wavelet analysis for the regional–residual and local separation
the Li–Oldenburg inversion algorithm for a cutoff of 0.85 g/cm3 is shown of potential field anomalies. Geophys. Prospect. 46 (5), 507–525.
Fedi, M., Rapolla, A., 1999. 3-D inversion of gravity and magnetic data with depth resolu-
in Fig. 11. For better illustration of density distribution inside the ore
tion. Geophysics 64 (2), 452–460.
body, we show the two sections in two different directions from the Grant, F.S., West, G.F., 1965. Interpretation Theory in Applied Geophysics. McGraw-Hill.
3D obtained model in Fig. 12a and b. The recovered model shows that Guillen, A., Menichetti, V., 1984. Gravity and magnetic inversion with minimization of a
the depth-to-top of the ore body is less than 20 m and the depth-to- specific functional. Geophysics 49 (8), 1354–1360.
Guspi, F., 1992. Three-dimensional Fourier gravity inversion with arbitrary density con-
base of the ore body is more than 200 m. trast. Geophysics 57 (1), 131–135.
Several boreholes are drilled near Noranda, Quebec on the Mubrun Last, B.J., Kubik, K., 1983. Compact gravity inversion. Geophysics 48 (6), 713–721.
massive sulfide ore body. BH.1, BH.2, BH.15 and BH.17 are four bore- Leblanc, G.E., Morris, W.A., 2001. Denoising of aeromagnetic data via the wavelet trans-
form. Geophysics 66 (6), 1793–1804.
holes located along the A–B profile (Fig. 9) which are shown in Fig. 13. Li, Y., Oldenburg, D.W., 1998. 3-D inversion of gravity data. Geophysics 63 (1), 109–119.
The depths to top and base of Mobrun sulfide body obtained by the Li, Y., Oldenburg, D.W., 2003. Fast inversion of large-scale magnetic data using wavelet
abovementioned methods approve very well drilling information transforms and a logarithmic barrier method. Geophys. J. Int. 152 (2), 251–265.
Mallat, S., 2008. A Wavelet Tour of Signal Processing: The Sparse Way. Elsevier Science.
which showed that the depth-to-top was 17 m and the depth-to-base Paterson, N.R., Reeves, C.V., 1985. Applications of gravity and magnetic surveys: the state‐
was 187 m (Roy et al., 2000). of‐the‐art in 1985. Geophysics 50 (12), 2558–2594.
Pedersen, L.B., 1977. Interpretation of potential field data a generalized inverse approach.
Geophys. Prospect. 25 (2), 199–230.
5. Conclusion Plouff, D., 1976. Gravity and magnetic fields of polygonal prisms and application to mag-
netic terrain corrections. Geophysics 41 (4), 727–741.
Ridsdill-Smith, T.A., Dentith, M.C., 1999. The wavelet transform in aeromagnetic process-
In this paper, we developed a new two dimensional continuous ing. Geophysics 64 (4), 1003–1013.
wavelet based method for determination of the location and the depth Roy, L., Agarwal, B.N.P., Shaw, R.K., 2000. A new concept in Euler deconvolution of isolated
gravity anomalies. Geophys. Prospect. 48 (3), 559–575.
to the top and base of gravity field sources. The method is compared
Sailhac, P., Galdeano, A., Gibert, D., Moreau, F., Delor, C., 2000. Identification of sources of
with Li and Oldenburg 3D inversion algorithm and is demonstrated potential fields with the continuous wavelet transform: complex wavelets and appli-
with a synthetic example and on real gravity data from near Noranda, cation to aeromagnetic profiles in French Guiana. J. Geophys. Res. Solid Earth 105
Quebec. To perform the 2DCWT method, we used the first horizontal (B8), 19455–19475.
Tai Sing, L., 1996. Image representation using 2D Gabor wavelets: pattern analysis and
derivatives of the anomaly of a buried cube with unit dimensions as machine intelligence. IEEE Trans. 18 (10), 959–971.
the mother wavelet. 2DCWT analysis on synthetic example can give Uieda, L., Barbosa, V.C.F., 2012. Robust 3D gravity gradient inversion by planting anoma-
useful information on the positions and depths of sources which have lous densities. Geophysics 77 (4), G55–G66.
Valee, M.A., Keating, P., Smith, R.S., St-Hilaire, C., 2004. Estimating depth and model type
a good agreement with 3D inversion density distribution. The method using the continuous wavelet transform of magnetic data. Geophysics 69 (1), 191–199.
is also applied on various levels of noise-contaminated synthetic data Zhdanov, M.S., 2002. Geophysical Inverse Theory and Regularization Problems. Elsevier
and the obtained results show that the proposed method has a good Science.