- Email: [email protected]
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
CHAPTER 11
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
Chapter Outline 11.1 Gravity Gradiometry Data 262 11.2 Migration of 3-D Gravity and Gravity Gradiometry Data 264 11.2.1 Adjoint Operators for Gravity and Gravity Gradiometry Inversion 264 11.2.2 Adjoint Operator for 3-D Gravity Fields 265 11.2.3 Adjoint Operator for 3-D Gravity Tensor Fields 266
11.3 Fast Density Imaging Based on Migration 267 11.3.1 11.3.2 11.3.3 11.3.4
Principles of Fast Inverse Imaging 267 Migration of Gravity and Gravity Tensor Fields and 3-D Density Imaging 268 Integrated Sensitivity of 3-D Gravity Fields 270 Integrated Sensitivity of 3-D Gravity Tensor Fields 272
11.4 Migration of Total Magnetic Intensity Data 275 11.4.1 11.4.2 11.4.3 11.4.4
Adjoint Operator for the Total Magnetic Intensity 275 Migration of the Total Magnetic Intensity 277 Integrated Sensitivity of the Total Magnetic Intensity 278 Model Study 279
References 281
Gravity gradiometry has become widely used in geophysical exploration since it can provide an independent measure of the subsurface density distribution. The advantage of gravity gradiometry over other gravity methods is that the data are extremely sensitive to local density anomalies within regional geological formations. High quality data can be acquired from either airborne or marine platforms over very large areas for relatively short time. A number of publications have discussed the use of the regularized inversion with both smooth and focusing stabilizers for the interpretation of gravity gradiometry data (e.g., Li, 2001; Zhdanov ˇ ˇ et al., 2004; Cuma et al., 2012; Cuma and Zhdanov, 2014). A variety of fast imaging techniques related to Euler decomposition have also been developed. Most of these are based on the superposition of analytical responses from specific sources (e.g., Fedi, 2007). These imaging methods typically estimate the positions and some parameters of the sources based on Inverse Theory and Applications in Geophysics. http://dx.doi.org/10.1016/B978-0-444-62674-5.00011-6 Copyright © 2015 Elsevier B.V. All rights reserved.
261
262
Chapter 11
field attenuation characteristics. In this chapter, I present a different approach to imaging based on the ideas of potential field migration as originally introduced by Zhdanov (2002). Migration can be mathematically described as the action of an adjoint operator on observed data. This concept has been long developed for seismic wavefields (e.g., Schneider, 1978; Berkhout, 1980; Claerbout, 1985; Tarantola, 1987) and was also developed for electromagnetic fields (Zhdanov, 1988, 2002, 2009), where the adjoint operators manifest themselves as the (backward) propagation of seismic or electromagnetic fields in reverse time. As applied to potential fields, such as gravity and magnetic fields, migration manifests itself as a special form of downward continuation of the potential field and/or its gradients. This downward continuation is applied to the auxiliary field obtained by moving the sources of the true observed field into the upper half-space as the mirror images of the true sources. This transformation results in extrapolation of the field downward and, contrary to conventional downward continuation, away from the mirror images of the sources. Thus migration is a stable transformation similar to conventional upward continuation. As I will demonstrate below, the migration field does contain remnant information about the original source distribution, which is why it can be used for subsurface imaging.
11.1 Gravity Gradiometry Data Recent technological developments make it possible to accurately measure all the independent tensor components of the gravity gradient field from a moving platform. It has been demonstrated that the use of gravity gradiometry data can significantly improve inversion results for mineral (e.g., Zhdanov et al., 2004) and hydrocarbon (Wan and Zhdanov, 2008) exploration. The technology that enables such rapid and accurate data acquisition has motivated the research to develop other methods for processing and interpreting gradiometry data. Moreover, these advancements continue to stimulate a growing interest in the application of gravity gradient data in geophysical exploration. Here, I provide a brief description of the gravity tensor components measured by gravity gradiometers. First, we know that the gravity field, g, must satisfy the following equations (Zhdanov, 1988): ∇ · g = −4πγρ,
∇ × g = 0,
(11.1)
where γ is the universal gravitational constant and ρ is the anomalous density distribution within a domain D. The solution of these equations is given by the equation: r − r g(r) = γ ρ(r ) dv , (11.2) 3 |r − r| D where integration is conducted over the variable r . The gravity field can be expressed by the gravity potential U(r) as g(r) = ∇U(r),
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
263
where U(r) = γ D
ρ(r ) dv . |r − r|3
(11.3)
The second spatial derivatives of the gravity potential U(r), gαβ (r) =
∂2 U(r), ∂α∂β
α, β = x, y, z,
(11.4)
form a symmetric gravity tensor: ⎤ gxx gxy gxz gˆ = ⎣ gyx gyy gyz ⎦ , gzx gzy gzz ⎡
where gαβ =
∂gα , ∂β
α, β = x, y, z.
(11.5)
The expressions for the gravity tensor components can be calculated based on (11.4) and (11.3): ρ(r ) K (r − r) dv , (11.6) gαβ (r) = γ − r|3 αβ |r D where the kernels, Kαβ , are equal to Kαβ (r − r) =
⎧ (α−α )(β−β ) ⎨ 3 |r −r|2 , α = β ⎩
)2 3 (α−α |r −r|2
− 1, α = β
,
α, β = x, y, z.
(11.7)
In addition to the gravity tensor components described by (11.6) and (11.7), the gravity gradiometers also measure the difference between the gradients: g =
1
gxx − gyy , 2
(11.8)
ρ(r ) K (r − r) dv , 3 |r − r|
(11.9)
3 (x − x)2 − (y − y)2 . 2 |r − r|2
(11.10)
which can be expressed as g = γ
D
where K (r − r) =
264
Chapter 11
11.2 Migration of 3-D Gravity and Gravity Gradiometry Data 11.2.1 Adjoint Operators for Gravity and Gravity Gradiometry Inversion Let us consider a problem of the inversion of gravity and gravity gradiometry data using gradient-type methods. We begin our discussion with gravity field inversion. According to (11.2), we have the following expression for the gravity field: ρ(r ) g Kα (r − r) dv , r ∈ / D, (11.11) gα (r) = Aα (ρ) = γ 3 D |r − r| g
where Aα (ρ), α = x, y, z, denotes the forward modeling operator for different gravity field components, and the kernel Kα (r − r) is equal to: Kα (r − r) = α − α,
α = x, y, z.
(11.12)
Let us assume that we have observed some component of the gravity field gobs α (r) on the observational surface S, and domain D is located in the lower half-space. The problem is to determine the density distribution, ρ(r ). We introduce a real Hilbert space G of gravity data with this metric: f (r)g(r) ds, f , g ∈ G, (11.13) (f , g)G = S
and a real Hilbert space M of models (density distribution, ρ(r )) with this metric: (σ , ρ)M =
σ (r )ρ(r ) dv ,
σ , ρ ∈ M.
(11.14)
D
For simplicity, we first ignore the ill-posedness of gravity inversion and reduce the inverse problem to the minimization of the misfit functional between the observed and predicted data: 2
obs obs = min . (11.15) φ(ρ) = gα − gα = Agα (ρ) − gobs α , Aα (ρ) − gα G
G
To solve the minimization problem (11.15) we calculate, as usual, the first variation of the misfit functional:
obs − g , A δφ(ρ) = δ Aα (ρ) − gobs (ρ) α α α G
= 2 δρ, A α Ag (ρ) − gobs = 2 (δρ, l (ρ))M , (11.16) α M
where the star denotes the adjoint operator, and l (ρ) is the direction of the steepest ascent at the point ρ of the space M of the model parameters: obs − g = A g . l (ρ) = A α Agα (ρ) − gobs α α α α
(11.17)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data The direction of the steepest ascent in the initial model with zero density is equal to:
. l (ρ = 0) = −A α gobs α
265
(11.18)
11.2.2 Adjoint Operator for 3-D Gravity Fields Let us find an explicit form of the adjoint operator Agα for the gravity fields:
(Aα (ρ) , f )G = ρ, A α (f ) M .
(11.19)
After some algebra, we obtain: Aα (ρ) f (r) ds = γ (Aα (ρ) , f )G =
ρ(r ) K (r − r) dv f (r) ds, − r|3 α |r S S D ρ(r ) ρ(r ) γ K (r − r)f (r) ds dv = − r|3 α |r D S f (r) = ρ, γ K (r − r) ds . (11.20) 3 α S |r − r| M
By comparing (11.19) and (11.20), we can see that f (r) ρ, Aα (f ) − γ K (r − r) ds = 0, 3 α S |r − r| M
(11.21)
where the star means adjoint operator. Equation (11.27) holds for any ρ, for example: f (r) ρ = Aα (f ) − γ K (r − r) ds. (11.22) 3 α S |r − r| Substituting (11.22) into (11.21), we obtain 2 f (r) A (f ) − γ Kα (r − r) ds = 0. α 3 S |r − r| M From the last equation, we find that adjoint operator, Agα , applied to some function f (r), is given by f (r) Kα (r − r) ds. (11.23) Aα (f ) = γ 3 S |r − r| It is shown above that the adjoint operator Agα for the gravity problem is equal to: f (r) K (r − r) ds. Aα (f ) = γ 3 α |r − r| S
(11.24)
266
Chapter 11
Therefore, according to (11.17), the direction of steepest ascent is equal to: gα (r) − gobs α (r) l (ρ) = γ Kα (r − r) ds, |r − r|3 S
(11.25)
where gα (r) is the predicted gravity field on the observation surface.
11.2.3 Adjoint Operator for 3-D Gravity Tensor Fields Let us consider now an adjoint operator for gravity tensor fields. Using the definitions (11.13) and (11.14) of inner products, and (11.6) for forward operator, we can rewrite (11.19) as follows:
ρ(r ) gαβ Aαβ (ρ) , f G = A (ρ) f (r) ds = γ K (r − r) dv f (r) ds − r|3 αβ |r S S D ρ(r ) ρ(r ) γ Kαβ (r − r)f (r) ds dv = 3 D S |r − r| f (r) = ρ, γ K (r − r) ds . (11.26) 3 αβ S |r − r| M By comparing (11.19) and (11.26), we can see that f (r) K (r − r) ds = 0, ρ, Aαβ (f ) − γ 3 αβ S |r − r| M
(11.27)
where the star means adjoint operator. Equation (11.27) holds for any ρ, for example: f (r) ρ = Aαβ (f ) − γ K (r − r) ds. (11.28) 3 αβ |r − r| S Substituting (11.28) into (11.27), we obtain 2 f (r) = 0. A (f ) − γ K (r − r) ds αβ αβ 3 S |r − r| M From the last equation, we find that the adjoint operator, A αβ , applied to some function f (r), is given by f (r) Aαβ (f ) = γ K (r − r) ds. (11.29) 3 αβ S |r − r| Similarly, in the case of the g component, the adjoint operator, A , applied to some function f (r), is given by f (r) K (r − r) ds. (11.30) A (f ) = γ 3 |r − r| S
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
267
Let us now consider gravity gradiometry inversion. We assume that we have observed some gravity gradients gobs αβ (r) on the observational surface S, and the domain D is located in the lower half-space. Again, the problem is to determine the density distribution, ρ(r ). We have the following expression for the gravity tensor field: ρ(r ) gαβ (r) = Aαβ (ρ) = γ K (r − r) dv , r ∈ / D. (11.31) − r|3 αβ |r D The adjoint operator for gravity gradients is given by (11.29). Therefore, according to (11.17), the direction of the steepest ascent is equal to: g (r) − gobs (r) αβ αβ l (ρ) = γ Kαβ (r − r) ds, (11.32) − r|3 |r S where gαβ (r) is the predicted gravity field on the observation surface. Note that, according to (11.17), the direction of the steepest ascent in the initial model with zero density is equal to: l (ρ = 0) = −A αβ gobs αβ .
(11.33)
We will use (11.18) and (11.33) for introducing the 3D gravity and gravity gradiometry migration fields.
11.3 Fast Density Imaging Based on Migration 11.3.1 Principles of Fast Inverse Imaging It is well known that the gravity and gravity gradiometry inversion is an ill-posed problem. We have to use regularization to generate a stable solution. In the framework of regularization theory, the solution of the inverse problem is reduced to the minimization of the Tikhonov parametric functional: Pλ (ρ) = ϕ (ρ) + λs (ρ) = min,
(11.34)
where ϕ (ρ) is a misfit functional between the theoretical predicted A (ρ) and the observed data d, and s (ρ) is a stabilizing functional. The traditional approach to implementing regularization is based on the class of solutions with smooth distributions of the model parameters. This approach is widely used in inversion. In particular, we can write the parametric functional (11.34) with a minimum norm stabilizer as follows:
2 Pλ (σ ) = Wd (A(ρ) − d)2G + λ Wm ρ − ρ apr G = min, (11.35) where Wd and Wm are positively determined linear data and model parameter weighting operators, respectively; λ is the regularization parameter; and ρ apr is an a priori model. Note
268
Chapter 11
that, below, we will use the weighting operators Wd and Wm described by multiplication of the observed data d by a positive scalar function wd (ri ), and by multiplication of the anomalous density ρ(r) by a positive scalar function, wm (r). In the framework of the Newton method one finds a solution to minimization problem (11.35) in one step: ρ 1 = ρ 0 + δρ.
(11.36)
The optimum step, δρ, satisfies the following normal equation (see Section 5.2.2): (Hm0 + λWm Wm )δρ = −l0 − λWm Wm (ρ 0 − ρ apr ),
(11.37)
where Hm0 is a Hessian operator and l0 is a gradient vector at the initial iteration σ 0 . Following Section 5.2.3, we can find an approximate solution to the normal (11.37) if we
assume that the regularization parameter λ is large enough to neglect the term Hm0 δρ with respect to the term λWm Wm δρ in (11.37). Then we obtain λWm Wm δρ ≈ −l0 − λWm Wm (ρ 0 − ρσ apr ). −1
Applying the inverse weighting operator Wm Wm to both sides of the last equation, we find:
−1 l0 − (ρ 0 − ρ apr ). δσ ≈ −λ−1 Wm Wm
(11.38)
Note that the coefficient λ−1 can be treated as a scaling factor that can be found by minimizing the misfit between the observed and predicted data. Equation (11.38) plays an important role in rapid imaging since it provides an approximate solution to the inverse problem:
−1 ρ 1 = ρ 0 + δρ ≈ ρ apr − λ−1 Wm Wm l0 , (11.39) which requires only gradient direction calculations. In the special case where the initial and a priori models of anomalous density are equal to zero (i.e., we start the inversion with no background density model), we arrive at the following extremely simple and important equation:
−1 ρ 1 ≈ −λ−1 Wm Wm l0 , (11.40) which serves as a basis for migration imaging.
11.3.2 Migration of Gravity and Gravity Tensor Fields and 3-D Density Imaging Following the principles introduced in Chapter 10 for 2-D migration, the migration gravity field, gm α (r), is introduced as a result of application of the adjoint gravity operator, Aα , to the observed component of the gravity field: gm α (r) = Aα gα ,
(11.41)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data where the adjoint operator Agα for the gravity problem is equal to: f (r) K (r − r) ds. A α (f ) = − r|3 α |r S
269
(11.42)
From the physical point of view, the migration field is obtained by moving the sources of the observed gravity field above the observational surface. Nevertheless, the migration field contains some remnant information about the original sources of the gravity anomaly. That is why it can be used in imaging the sources of the gravity field. In a similar way, we can introduce a migration gravity tensor field gm αβ (r) and use the following notations for the components of this tensor field: gm αβ (r) = Aαβ gαβ ,
(11.43)
gm (r) = A g ,
(11.44)
where the adjoint operators, A αβ and A (f ), applied to some function f (r), are given by the formulas f (r) K (r − r) ds, (11.45) Aαβ (f ) = 3 αβ S |r − r| f (r) K (r − r) ds. (11.46) A (f ) = 3 |r − r| S We should note, however, that the direct migration of the observed gravity and/or gravity tensor fields does not produce an adequate image of the subsurface density distribution because the migration fields rapidly attenuate with the depth, as one can see from expressions (11.42), (11.45), and (11.46). In order to image the sources of the gravity fields at their correct location, one should apply an appropriate spatial weighting operator to the migration fields. This weighting operator is constructed based on the integrated sensitivity of the data to the density. We shall now apply (11.40) for rapid imaging. Substituting (11.18) for the direction of steepest ascent into (11.40), and taking into account (11.41), one can find a distribution of the density of the gravity field sources, described by the following expression:
−1 l0 , (11.47) ρ 1 ≈ −λ−1 Wm Wm
−1 m Aα gα = kα w−2 (11.48) ραm (r) = λ−1 Wm Wm α (z) gα (r), where unknown coefficient kα = λ−1 can be determined by a linear line search (Zhdanov, 2002) according to the following: w 2 A gα α M kα = (11.49) , Aw Aw gα 2 α α D −1 Aw α = Aα Wα ,
(11.50)
270
Chapter 11
and the linear weighting operator Wm = Wα is selected as a linear operator of the multiplication of the density ρ by a function, wα , equal to the square root of the integrated sensitivity of the complex intensity of the gravity field, Sα : wα = Sα . (11.51)
11.3.3 Integrated Sensitivity of 3-D Gravity Fields The integrated sensitivity, in accordance with the definition, is calculated as: obs δg α G , Sα = δρ
(11.52)
where δgobs α is the perturbation of the gravity field resulting from a local perturbation of the density, δρ(r) = ρ(r) dv, within a differential element of volume dv, located at the point r = (x, y, z) of the lower half-plane (z > 0): obs δgobs α = δgα (r ) = γ
ρ(r) Kα (r − r ) dv. |r − r |3
(11.53)
Substituting (11.53) into (11.52), we find 1 Sα = ρ(r) dv
=γ S
S
2 ds δgobs α (r )
Kα2 (r − r ) ds , |r − r |6
(11.54)
where S is the surface of observations of the gravity field. In particular, if the profile of observations coincides with the horizontal plane, z = 0, the definite integral in (11.54) can be evaluated as: ∞ K 2 (r) α Sα (z) = γ dx dy. (11.55) 6 |r| −∞ Using a polar system of coordinates, then x = R cos ϕ,
y = R sin ϕ,
ds = dx dy = R dR dϕ,
3 3 |r|6 = x2 + y2 + z2 = R2 + z2 , Kα (r − r) = α − α,
α = x, y, z,
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data and we can express (11.54) for α = x, y, z, as follows: ∞ 4π 2 R cos2 ϕ Sx (z) = γ 3 R dR dϕ, 0 0 R2 + z2 Sy (z) = γ
∞ 4π
0
Sz (z) = γ |z|
∞ 4π 0
Taking into account that: 4π 2 cos ϕ dϕ = 0
4π
(11.56)
R2 sin2 ϕ 3 R dR dϕ, R2 + z2
0
0
(11.57)
R dR dϕ 3 . R2 + z2
(11.58)
sin ϕ dϕ = 2π, 2
4π
and
0
271
dϕ = 4π,
(11.59)
0
we obtain the following: Sx (z) = Sy (z) = γ 2π
∞
0
Sz (z) = γ |z| 4π
4π
0
R2 3 R dR, R2 + z2
R dR R2 + z2
3 .
Calculating the tabulated integrals in (11.60) and (11.61), we finally find: √ 1 π 1 , Sz (z) = γ π . Sx (z) = Sy (z) = γ |z| 2 |z|
(11.60)
(11.61)
(11.62)
Thus, the integrated sensitivity, Sα , can be computed from the following equations: Sα = c α
1 , |z|
z < 0, α = x, y, z,
(11.63)
and cα are the corresponding constants for different components equal to: √ π cx = cy = γ , cz = γ π. 2 g
Equation (11.48) is called a migration density, ρm (ζ ): m ραm (r) = kα w−2 αg (z) gα (r),
(11.64)
272
Chapter 11
and is proportional to the weighted migration field, gm α: kα |z| gm ραm (r) = α (r), cα where
gm α (r) =
S
gα (r ) Kα (r − r ) ds . |r − r |3
(11.65)
(11.66)
Thus, the migration transform with spatial weighting provides a stable algorithm for calculating the migration density. In a similar way, we can introduce a migration density based on the gravity tensor migration: m m (r) = kαβ w−2 ραβ αβ (z) gαβ (r),
where T kαβ
2 w Aαβ gαβ M = , w W 2 Aαβ Aαβ gαβ D
m m ρ (r) = γ k w−2 (z) g (r),
T k
w 2 A g M = , Aw AW g 2 D
(11.67)
(11.68)
where functions wαβ and w are equal to the square root of the integrated sensitivity of the gravity tensor fields, Sαβ and S , respectively: wαβ = Sαβ , w = S . (11.69)
11.3.4 Integrated Sensitivity of 3-D Gravity Tensor Fields The integrated sensitivity of the tensor field components is calculated as: obs δgαβ G , S= δρ
(11.70)
where δgobs αβ is the perturbation of the corresponding component of the gravity tensor field resulting from a local perturbation of the density, δρ(r) = ρ(r) dv, within a differential element of volume dv, located at the point r = (x, y, z) of the lower half-plane (z > 0): ρ(r) obs Kαβ (r − r ) dv. (11.71) δgobs αβ = δgαβ (r ) = γ |r − r |3 Substituting (11.71) into (11.70), we find the following: 2
1 ) (r ds δgobs Sαβ = αβ ρ(r) dv S K 2 (r − r ) αβ =γ ds , 6 S |r − r |
(11.72)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
273
where S is the surface of observations of the gravity tensor field. In particular, if the profile of observations coincides with the horizontal plane, z = 0, the definite integral in (11.72) can be evaluated as: ∞ K 2 (r) αβ Sαβ (z) = γ dx dy. (11.73) 6 −∞ |r| Using a polar system of coordinates, then: x = R cos ϕ,
y = R sin ϕ,
ds = dx dy = R dR dϕ,
3 3 |r|6 = x2 + y2 + z2 = R2 + z2 , Kαβ (r) =
αβ 3 |r| 2 , α = β
3 (α) − 1, α = β |r|2 2
,
α, β = x, y, z.
We can express (11.72) for α, β = x, y, z, as: 2 ∞
4π 2z2 − R2 R dR dϕ, Szz (z) = γ 5 0 0 R2 + z2 ∞ 4π R2 cos2 ϕ dϕ, Sxz (z) = Syz (z) = γ 3 |z| 5 R dR 2 2 0 0 R +z 2 ∞ 4π 1 R2 cos2 ϕ −1 Sxx (z) = Sxx (z) = γ 3 2 3 R dR dϕ. 2 2 R +z 0 0 R + z2 Taking into account that 4π 3 cos4 ϕ dϕ = π, 2 0
4π 0
cos ϕ dϕ = 2π, 2
4π
dϕ = 4π,
(11.74)
(11.75)
(11.76)
(11.77)
0
we obtain the following:
Szz (z) = γ 4π
2 2z2 − R2 5 R dR, 0 R2 + z2 ∞ R2 Sxz (z) = Syz (z) = γ 3 |z| 2π 5 R dR, 0 R2 + z2 ∞ R4 1 Sxx (z) = Syy (z) = 3γ 4π 2 3 R dR. 0 R2 + z2 R2 + z2 ∞
(11.78)
(11.79)
(11.80)
274
Chapter 11
Calculating the tabulated integrals in (11.78) to (11.80), we finally find: √ 3π Szz (z) = γ , z < 0, 2z2 √ 3π , z < 0, Sxz (z) = Syz (z) = γ 2z2 √ 3 π , z < 0. Sxx (z) = Syy (z) = γ 4z2
(11.81) (11.82) (11.83)
Thus, the integrated sensitivity of the gravity tensor field is calculated from the following equations: 1 1 Sαβ = cαβ 2 , S = c 2 , (11.84) z z where cαβ are the corresponding constants for different components equal to: √ √ 3π 3 π czz = czx = czy = γ , cxx = cyy = γ . 2 4 Expression (11.67) is called a tensor field migration density. It is proportional to the magnitude of the weighted tensor migration field gm αβ . Thus, migration transformation provides a stable algorithm for calculating migration density. Substituting (11.84) for the weighting function wT back into (11.69) and (11.67), we find that: kαβ 2 m k 2 m m m (r) = z gαβ (r), ρ (r) = z g (r), (11.85) ραβ cαβ c where
gαβ (r ) K (r − r ) ds , |3 αβ |r − r S g (r ) gm (r) = K (r − r ) ds . |3 |r − r S
gm αβ (r)
=
(11.86) (11.87)
Finally, we can consider joint migration of several components of the gravity tensor. For example, we can jointly migrate gzz , gxx , gzx , and the g components and find the corresponding migration density as per the following: kxx m kzx m k m m 2 kzz m ραβ (r) = z g (r) + g (r) + g (r) + g (r) . (11.88) czz zz cxx xx czx zx c Note that (11.88) can be simplified to: m m m m ραβ (r) = z2 azz gm zz (r) + axx gxx (r) + azx gzx (r) + a g (r) ,
(11.89)
where azz , axx , azx , and a can be treated as the weights of the corresponding migration fields in the density model, which can be empirically determined from model studies.
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
275
Thus, for gravity fields and their gradients, we have shown that potential field migration is an integral transformation of the gravity field and/or gradients into 3-D density distributions. Potential field migration is very fast and stable, and the results are effectively equivalent to those obtained from 3-D regularized inversion with smooth stabilizers (Zhdanov et al., 2010, 2011).
11.4 Migration of Total Magnetic Intensity Data The earth’s magnetic field is the vector sum of contributions from two main sources—a background field due to the dynamo effect of the earth’s liquid core, and anomalous fields due to magnetic rocks and minerals above the Curie isotherm. Magnetic vector data, measured by orthogonal fluxgate magnetometers, are dominated by the earth’s background field, and are thus very sensitive to instrument orientation. The development of reliable and low cost optically pumped magnetometers in the 1960s enabled direct measurement of the total magnetic intensity (TMI) regardless of instrument orientation. It is now routine practice that every airborne geophysical survey produces TMI data as a standard deliverable. Relative to the millions of line kilometers of TMI data acquired each year, 3-D inversions are rarely performed. This is a reflection of the limited capability of existing 3-D inversion software to invert entire surveys to 3-D earth models with sufficient resolution in sufficient time so as to affect exploration decisions. Most interpretations are simply based on picking lineaments from the maps of the TMI’s first vertical derivative. Structural interpretations are usually based on some kind of Euler deconvolution, eigenvector, wavelet, analytic signal, or depth-from-extreme-points methods. While such methods may provide information about the sources of the TMI field, it is not immediately obvious how that information can be quantified in terms of a 3-D susceptibility model. To that end, Zhdanov et al. (2012) developed a method of rapid 3-D imaging of the TMI field data based on the principles of potential field migration discussed above.
11.4.1 Adjoint Operator for the Total Magnetic Intensity We measure TMI data on a surface S above a domain V that is filled by magnetic sources with the intensity of magnetization I(r). The problem is to determine the magnetic susceptibility, χ(r). In what follows, we adopt the common assumptions that there is no remanent magnetization, that the self-demagnetization effect is negligible, and that the magnetic susceptibility is isotropic. Under such assumptions, the intensity of magnetization is linearly related to an inducing magnetic field, H0 (r), through the magnetic susceptibility: I(r) = χ(r)H0 (r), where r is the radius vector of a point within the volume V.
(11.90)
276
Chapter 11
It is well known that the anomalous scalar TMI data T generated by the magnetic sources within the volume V can be represented by the linear operator equation (Li and Oldenburg, 1996; Portniaguine and Zhdanov, 2002): χ(r) T(r ) = A (χ) = H0 K(r − r)dv, r ∈ / V, (11.91) 3 |r | − r V where A (χ) denotes the forward modeling operator, H 0 is the magnitude of the inducing field, and l is a unit vector in the direction of magnetization: H0 (r) = H 0 l(r), and K is the TMI kernel:
2 3 l · r −r K(r − r) = − 1. (11.92) |r −r|2 We introduce a real Hilbert space D of the TMI data with the following metric: f (r )g(r ) ds , f , g ∈ D, (f , g)D =
(11.93)
S
and also introduce a real Hilbert space M of models (i.e., magnetic susceptibility, χ) with the following metric: χ(r)η(r) dv, χ, η ∈ M. (11.94) (χ, η)M = V
According to definition, the TMI adjoint operator, A , satisfies to the following equation:
(11.95) (A (χ) , f )D = χ, A (f ) M . Let us write the left-hand side of this equation in explicit form: A (χ) f (r ) ds (A (χ) , f )D = S
χ(r) K(r − r) dv f (r) ds 3 − r| S V f (r ) χ(r ) H0 K(r − r)f (r) ds dv = 3 D S |r − r| f (r ) = χ, H0 K(r − r) ds . 3 S |r − r| M = H0
|r
By comparing (11.95) and (11.96), we can see that: f (r ) K(r − r) ds = 0, χ, A (f ) − H0 3 S |r − r| M
(11.96)
(11.97)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data where “ ” denotes the adjoint operator. Equation (11.97) holds for any χ, e.g.: f (r ) χ = A (f ) − H0 K(r − r) ds. 3 S |r − r| Substituting (11.98) into (11.97), we obtain: 2 f (r ) A (f ) − H0 = 0. K(r − r) ds 3 S |r − r| M
277
(11.98)
(11.99)
From (11.99), we find that the adjoint operator, A , applied to any function f (r ), is given by the following expression: f (r ) A (f ) = H0 K(r − r) ds. (11.100) 3 S |r − r| We note that, the TMI adjoint operator produces a function that is analytical everywhere in the subsurface. The significance of this will become apparent momentarily.
11.4.2 Migration of the Total Magnetic Intensity We have established in the previous sections that migration is the action of the adjoint operator on the observed data. It follows that the migration TMI field, T m (r), is introduced as the action of the TMI adjoint operator, A , on the measured TMI field, T: T m (r) = A T(r ).
(11.101)
From the physical point of view, the migration TMI field is obtained by moving the sources of the observed TMI field above the surface, S. Nevertheless, the migration TMI field contains information about the original sources of the TMI field, which is why it can be used for imaging. We note, however, that direct migration of the measured TMI field does not produce an adequate image of the susceptibility distribution because the migration fields rapidly attenuate with the depth. In order to image the sources of the TMI field at their correct locations, one should apply an appropriate spatial weighting operator to the migration TMI field. Similar to 3-D gravity migration, the migration magnetic susceptibility can be computed from:
−1 χ m (r) = k W W A T = kw−2 (z) T m (r). (11.102) In the last formula, the unknown coefficient k can be determined by a linear line search according to the following: k=
Aw T2M Aw Aw T2D
,
(11.103)
where Aw =AW −1 ,
(11.104)
278
Chapter 11
and the linear weighting operator, W, is selected as an operator of multiplication of the susceptibility χ by a depth weighting function, w, which is equal to the square root of the integrated sensitivity, S (z), of the TMI data: δT obs D w(z) = (11.105) = S (z), δχ
11.4.3 Integrated Sensitivity of the Total Magnetic Intensity The integrated sensitivity of the TMI is calculated as: δT obs D S= , δχ
(11.106)
where δT obs is the perturbation of the corresponding component of the TMI resulting from a local perturbation of the magnetic susceptibility, δχ (r) = χ(r) dv, within a differential element of volume dv, located at the point r = (x, y, z) of the lower half-space (z > 0): δT obs = δT obs (r ) = H0 Substituting (11.107) into (11.106), we obtain: S = H0 S
χ(r) K(r − r ) dv. |r − r |3
K 2 (r − r ) ds , |r − r |6
(11.107)
(11.108)
where S is the surface of observations of the TMI. In particular, if the surface of observations coincides with the horizontal plane, z = 0, the definite integral in (11.108) can be evaluated as follows: 2 2 ∞ 4π
lx R cos ϕ + ly R sin ϕ + lz z R dR dϕ S (z) = H0 (11.109) 3 −1 3 , 2 2 R +z 0 0 R2 + z2 where we use a polar system of coordinates (ϕ, R, z): x = R cos ϕ, Taking into account that 4π 3 cos4 ϕ dϕ = π, 2 0
y = R sin ϕ, 0
4π
ds = dx dy = R dR dϕ.
cos ϕ dϕ = 2π, 2
0
4π
dϕ = 4π,
(11.110)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
279
we obtain the following expression:
S (z) = 6H0 π 0
∞
R5
5 dR. R2 + z2
(11.111)
Calculating the tabulated integral in (11.111), we finally find a simple analytical formula for integrated sensitivity of the TMI field: √ 3 π (11.112) S (z) = H0 2 , z < 0. 4z Substituting (11.101), (11.105), and (11.112) into (11.102), we find that: T(r ) 4z2 k m K(r − r ) ds . (11.113) χ (r) = √ − r|3 |r 3H02 π S The migration magnetic susceptibility (11.113) is proportional to the magnitude of the weighted migration field, T m (r), which is analytical everywhere in the subsurface, implying that migration is a well posed and stable transform. It can be shown that, if the surface S coincides with the horizontal axis, z = 0, migration is mathematically equivalent to a form of downward continuation of a function analytical everywhere in the lower half-space. It should be noted that the downward continuation of the measured TMI field and the migration TMI field are significantly different. The migration TMI field has singular points in the upper half-space; therefore its downward continuation is a well posed and stable transform. The downward continuation of the measured TMI field has singular points in the lower half-space associated with its sources, so its downward continuation can only be extended down to these singularities, making it an ill-posed and unstable transform (e.g., Strakhov, 1970; Zhdanov, 1988).
11.4.4 Model Study To investigate the performance of 3-D TMI migration, we considered two discrete bodies with a susceptibility of 0.05 buried about 110 m below variable topography (Figure 11.1a). The inducing field had an inclination of 75 degrees and a declination of 25 degrees. The synthetic TMI data were computed on a 20 m regular grid draped over a curvilinear surface covering an area of one square kilometer. Figure 11.1b shows an example of a profile of the TMI data over both bodies. We applied 3-D migration to the entire synthetic TMI dataset, and calculated the migration magnetic susceptibility according to (11.112). As an example, Figure 11.1c shows a vertical cross section through the 3-D migration magnetic susceptibility model. The correct locations of the two bodies can be clearly determined. In order to demonstrate the robustness of migration to noise, we contaminated the data with 20% random Gaussian noise, as shown in Figure 11.1d. Figure 11.1e shows the same vertical cross section of results for the 3D migration of the noisy data. As expected, 3-D migration produced a very robust image of the susceptibility distribution.
(a)
(b)
(c)
(d)
(e) Figure 11.1 (a) Perspective of 3-D model consisting of two blocks of 0.05 susceptibility. Observation sites (dark gray dots) are located across a variable topography. (b) Synthetic TMI data along the profile y = 0 m. (c) Vertical cross section of susceptibility along y = 0 from 3-D migration of synthetic TMI data. (d) Synthetic TMI data with 20% random Gaussian noise along the profile y = 0 m. (e) Vertical cross section of susceptibility along y = 0 from 3-D migration of noisy synthetic TMI data.
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
281
References Berkhout, A.J., 1980. Seismic Migration. Elsevier, Amsterdam. Claerbout, J.F., 1985. Imaging the Earth’s Interior. Blackwell Scientific Publications, Oxford. ˇ Cuma, M., Wilson, G.A., Zhdanov, M.S., 2012. Large-scale 3D inversion of potential field data. Geophys. Prospect. 60 (6), 1186-1199. ˇ Cuma, M., Zhdanov, M.S., 2014. Massively parallel regularized 3D inversion of potential fields on CPUs and GPUs. Comput. Geosci. 62, 80-87. Fedi, M., 2007. DEXP: a fast method to determine the depth to the sources of potential fields. Geophysics 72, L1-L11. Li, Y., Oldenburg, D.W., 1996. 3-D inversion of magnetic data. Geophysics 61, 394-408. Li, Y., 2001. 3D inversion of gravity gradiometer data. In: 71st Annual International Meeting, SEG, Expanded Abstracts, pp. 1470-1473. Portniaguine, O., Zhdanov, M.S., 2002. 3-D magnetic inversion with data compression and image focusing. Geophysics 67, 1532-1541. Schneider, W.A., 1978. Integral formulation for migration in two and three dimensions. Geophysics 43 (2), 49-76. Strakhov, V.N., 1970. Some aspects of the plane inverse problem of magnetic potential. Izv. Acad. Nauk SSSR Fiz. Zemli 9, 31-41 (in Russian). Tarantola, A., 1987. Inverse Problem Theory. Elsevier, Amsterdam. Wan, L., Zhdanov, M.S., 2008. Focusing inversion of marine full-tensor gradiometry data in offshore geophysical exploration. In: 76th Annual International Meeting, SEG, Expanded Abstracts, pp. 751-754. Zhdanov, M.S., 1988. Integral Transforms in Geophysics. Springer-Verlag, Berlin. Zhdanov, M.S., 2002. Geophysical Inverse Theory and Regularization Problems. Elsevier, Amsterdam. Zhdanov, M.S., Ellis, R.G., Mukherjee, S., 2004. Regularized focusing inversion of 3D gravity tensor data. Geophysics 69, 925-937. Zhdanov, M.S., 2009. Geophysical Electromagnetic Theory and Methods. Elsevier, Amsterdam. Zhdanov, M.S., Liu, X., Wilson, G., 2010. Potential field migration for rapid 3D imaging of gravity gradiometry surveys. First Break 28 (11), 47-51. Zhdanov, M.S., Liu, X., Wilson, G.A., Wan, L., 2011. Potential field migration for rapid imaging of gravity gradiometry data. Geophys. Prospect. 59, 1052-1071. Zhdanov, M.S., Liu, X., Wilson, G.A., Wan, L., 2012. 3D migration for rapid imaging of total-magnetic-intensity data. Geophysics 77 (2), J1-J5.
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
Chapter Outline 11.1 Gravity Gradiometry Data 262 11.2 Migration of 3-D Gravity and Gravity Gradiometry Data 264 11.2.1 Adjoint Operators for Gravity and Gravity Gradiometry Inversion 264 11.2.2 Adjoint Operator for 3-D Gravity Fields 265 11.2.3 Adjoint Operator for 3-D Gravity Tensor Fields 266
11.3 Fast Density Imaging Based on Migration 267 11.3.1 11.3.2 11.3.3 11.3.4
Principles of Fast Inverse Imaging 267 Migration of Gravity and Gravity Tensor Fields and 3-D Density Imaging 268 Integrated Sensitivity of 3-D Gravity Fields 270 Integrated Sensitivity of 3-D Gravity Tensor Fields 272
11.4 Migration of Total Magnetic Intensity Data 275 11.4.1 11.4.2 11.4.3 11.4.4
Adjoint Operator for the Total Magnetic Intensity 275 Migration of the Total Magnetic Intensity 277 Integrated Sensitivity of the Total Magnetic Intensity 278 Model Study 279
References 281
Gravity gradiometry has become widely used in geophysical exploration since it can provide an independent measure of the subsurface density distribution. The advantage of gravity gradiometry over other gravity methods is that the data are extremely sensitive to local density anomalies within regional geological formations. High quality data can be acquired from either airborne or marine platforms over very large areas for relatively short time. A number of publications have discussed the use of the regularized inversion with both smooth and focusing stabilizers for the interpretation of gravity gradiometry data (e.g., Li, 2001; Zhdanov ˇ ˇ et al., 2004; Cuma et al., 2012; Cuma and Zhdanov, 2014). A variety of fast imaging techniques related to Euler decomposition have also been developed. Most of these are based on the superposition of analytical responses from specific sources (e.g., Fedi, 2007). These imaging methods typically estimate the positions and some parameters of the sources based on Inverse Theory and Applications in Geophysics. http://dx.doi.org/10.1016/B978-0-444-62674-5.00011-6 Copyright © 2015 Elsevier B.V. All rights reserved.
261
262
Chapter 11
field attenuation characteristics. In this chapter, I present a different approach to imaging based on the ideas of potential field migration as originally introduced by Zhdanov (2002). Migration can be mathematically described as the action of an adjoint operator on observed data. This concept has been long developed for seismic wavefields (e.g., Schneider, 1978; Berkhout, 1980; Claerbout, 1985; Tarantola, 1987) and was also developed for electromagnetic fields (Zhdanov, 1988, 2002, 2009), where the adjoint operators manifest themselves as the (backward) propagation of seismic or electromagnetic fields in reverse time. As applied to potential fields, such as gravity and magnetic fields, migration manifests itself as a special form of downward continuation of the potential field and/or its gradients. This downward continuation is applied to the auxiliary field obtained by moving the sources of the true observed field into the upper half-space as the mirror images of the true sources. This transformation results in extrapolation of the field downward and, contrary to conventional downward continuation, away from the mirror images of the sources. Thus migration is a stable transformation similar to conventional upward continuation. As I will demonstrate below, the migration field does contain remnant information about the original source distribution, which is why it can be used for subsurface imaging.
11.1 Gravity Gradiometry Data Recent technological developments make it possible to accurately measure all the independent tensor components of the gravity gradient field from a moving platform. It has been demonstrated that the use of gravity gradiometry data can significantly improve inversion results for mineral (e.g., Zhdanov et al., 2004) and hydrocarbon (Wan and Zhdanov, 2008) exploration. The technology that enables such rapid and accurate data acquisition has motivated the research to develop other methods for processing and interpreting gradiometry data. Moreover, these advancements continue to stimulate a growing interest in the application of gravity gradient data in geophysical exploration. Here, I provide a brief description of the gravity tensor components measured by gravity gradiometers. First, we know that the gravity field, g, must satisfy the following equations (Zhdanov, 1988): ∇ · g = −4πγρ,
∇ × g = 0,
(11.1)
where γ is the universal gravitational constant and ρ is the anomalous density distribution within a domain D. The solution of these equations is given by the equation: r − r g(r) = γ ρ(r ) dv , (11.2) 3 |r − r| D where integration is conducted over the variable r . The gravity field can be expressed by the gravity potential U(r) as g(r) = ∇U(r),
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
263
where U(r) = γ D
ρ(r ) dv . |r − r|3
(11.3)
The second spatial derivatives of the gravity potential U(r), gαβ (r) =
∂2 U(r), ∂α∂β
α, β = x, y, z,
(11.4)
form a symmetric gravity tensor: ⎤ gxx gxy gxz gˆ = ⎣ gyx gyy gyz ⎦ , gzx gzy gzz ⎡
where gαβ =
∂gα , ∂β
α, β = x, y, z.
(11.5)
The expressions for the gravity tensor components can be calculated based on (11.4) and (11.3): ρ(r ) K (r − r) dv , (11.6) gαβ (r) = γ − r|3 αβ |r D where the kernels, Kαβ , are equal to Kαβ (r − r) =
⎧ (α−α )(β−β ) ⎨ 3 |r −r|2 , α = β ⎩
)2 3 (α−α |r −r|2
− 1, α = β
,
α, β = x, y, z.
(11.7)
In addition to the gravity tensor components described by (11.6) and (11.7), the gravity gradiometers also measure the difference between the gradients: g =
1
gxx − gyy , 2
(11.8)
ρ(r ) K (r − r) dv , 3 |r − r|
(11.9)
3 (x − x)2 − (y − y)2 . 2 |r − r|2
(11.10)
which can be expressed as g = γ
D
where K (r − r) =
264
Chapter 11
11.2 Migration of 3-D Gravity and Gravity Gradiometry Data 11.2.1 Adjoint Operators for Gravity and Gravity Gradiometry Inversion Let us consider a problem of the inversion of gravity and gravity gradiometry data using gradient-type methods. We begin our discussion with gravity field inversion. According to (11.2), we have the following expression for the gravity field: ρ(r ) g Kα (r − r) dv , r ∈ / D, (11.11) gα (r) = Aα (ρ) = γ 3 D |r − r| g
where Aα (ρ), α = x, y, z, denotes the forward modeling operator for different gravity field components, and the kernel Kα (r − r) is equal to: Kα (r − r) = α − α,
α = x, y, z.
(11.12)
Let us assume that we have observed some component of the gravity field gobs α (r) on the observational surface S, and domain D is located in the lower half-space. The problem is to determine the density distribution, ρ(r ). We introduce a real Hilbert space G of gravity data with this metric: f (r)g(r) ds, f , g ∈ G, (11.13) (f , g)G = S
and a real Hilbert space M of models (density distribution, ρ(r )) with this metric: (σ , ρ)M =
σ (r )ρ(r ) dv ,
σ , ρ ∈ M.
(11.14)
D
For simplicity, we first ignore the ill-posedness of gravity inversion and reduce the inverse problem to the minimization of the misfit functional between the observed and predicted data: 2
obs obs = min . (11.15) φ(ρ) = gα − gα = Agα (ρ) − gobs α , Aα (ρ) − gα G
G
To solve the minimization problem (11.15) we calculate, as usual, the first variation of the misfit functional:
obs − g , A δφ(ρ) = δ Aα (ρ) − gobs (ρ) α α α G
= 2 δρ, A α Ag (ρ) − gobs = 2 (δρ, l (ρ))M , (11.16) α M
where the star denotes the adjoint operator, and l (ρ) is the direction of the steepest ascent at the point ρ of the space M of the model parameters: obs − g = A g . l (ρ) = A α Agα (ρ) − gobs α α α α
(11.17)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data The direction of the steepest ascent in the initial model with zero density is equal to:
. l (ρ = 0) = −A α gobs α
265
(11.18)
11.2.2 Adjoint Operator for 3-D Gravity Fields Let us find an explicit form of the adjoint operator Agα for the gravity fields:
(Aα (ρ) , f )G = ρ, A α (f ) M .
(11.19)
After some algebra, we obtain: Aα (ρ) f (r) ds = γ (Aα (ρ) , f )G =
ρ(r ) K (r − r) dv f (r) ds, − r|3 α |r S S D ρ(r ) ρ(r ) γ K (r − r)f (r) ds dv = − r|3 α |r D S f (r) = ρ, γ K (r − r) ds . (11.20) 3 α S |r − r| M
By comparing (11.19) and (11.20), we can see that f (r) ρ, Aα (f ) − γ K (r − r) ds = 0, 3 α S |r − r| M
(11.21)
where the star means adjoint operator. Equation (11.27) holds for any ρ, for example: f (r) ρ = Aα (f ) − γ K (r − r) ds. (11.22) 3 α S |r − r| Substituting (11.22) into (11.21), we obtain 2 f (r) A (f ) − γ Kα (r − r) ds = 0. α 3 S |r − r| M From the last equation, we find that adjoint operator, Agα , applied to some function f (r), is given by f (r) Kα (r − r) ds. (11.23) Aα (f ) = γ 3 S |r − r| It is shown above that the adjoint operator Agα for the gravity problem is equal to: f (r) K (r − r) ds. Aα (f ) = γ 3 α |r − r| S
(11.24)
266
Chapter 11
Therefore, according to (11.17), the direction of steepest ascent is equal to: gα (r) − gobs α (r) l (ρ) = γ Kα (r − r) ds, |r − r|3 S
(11.25)
where gα (r) is the predicted gravity field on the observation surface.
11.2.3 Adjoint Operator for 3-D Gravity Tensor Fields Let us consider now an adjoint operator for gravity tensor fields. Using the definitions (11.13) and (11.14) of inner products, and (11.6) for forward operator, we can rewrite (11.19) as follows:
ρ(r ) gαβ Aαβ (ρ) , f G = A (ρ) f (r) ds = γ K (r − r) dv f (r) ds − r|3 αβ |r S S D ρ(r ) ρ(r ) γ Kαβ (r − r)f (r) ds dv = 3 D S |r − r| f (r) = ρ, γ K (r − r) ds . (11.26) 3 αβ S |r − r| M By comparing (11.19) and (11.26), we can see that f (r) K (r − r) ds = 0, ρ, Aαβ (f ) − γ 3 αβ S |r − r| M
(11.27)
where the star means adjoint operator. Equation (11.27) holds for any ρ, for example: f (r) ρ = Aαβ (f ) − γ K (r − r) ds. (11.28) 3 αβ |r − r| S Substituting (11.28) into (11.27), we obtain 2 f (r) = 0. A (f ) − γ K (r − r) ds αβ αβ 3 S |r − r| M From the last equation, we find that the adjoint operator, A αβ , applied to some function f (r), is given by f (r) Aαβ (f ) = γ K (r − r) ds. (11.29) 3 αβ S |r − r| Similarly, in the case of the g component, the adjoint operator, A , applied to some function f (r), is given by f (r) K (r − r) ds. (11.30) A (f ) = γ 3 |r − r| S
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
267
Let us now consider gravity gradiometry inversion. We assume that we have observed some gravity gradients gobs αβ (r) on the observational surface S, and the domain D is located in the lower half-space. Again, the problem is to determine the density distribution, ρ(r ). We have the following expression for the gravity tensor field: ρ(r ) gαβ (r) = Aαβ (ρ) = γ K (r − r) dv , r ∈ / D. (11.31) − r|3 αβ |r D The adjoint operator for gravity gradients is given by (11.29). Therefore, according to (11.17), the direction of the steepest ascent is equal to: g (r) − gobs (r) αβ αβ l (ρ) = γ Kαβ (r − r) ds, (11.32) − r|3 |r S where gαβ (r) is the predicted gravity field on the observation surface. Note that, according to (11.17), the direction of the steepest ascent in the initial model with zero density is equal to: l (ρ = 0) = −A αβ gobs αβ .
(11.33)
We will use (11.18) and (11.33) for introducing the 3D gravity and gravity gradiometry migration fields.
11.3 Fast Density Imaging Based on Migration 11.3.1 Principles of Fast Inverse Imaging It is well known that the gravity and gravity gradiometry inversion is an ill-posed problem. We have to use regularization to generate a stable solution. In the framework of regularization theory, the solution of the inverse problem is reduced to the minimization of the Tikhonov parametric functional: Pλ (ρ) = ϕ (ρ) + λs (ρ) = min,
(11.34)
where ϕ (ρ) is a misfit functional between the theoretical predicted A (ρ) and the observed data d, and s (ρ) is a stabilizing functional. The traditional approach to implementing regularization is based on the class of solutions with smooth distributions of the model parameters. This approach is widely used in inversion. In particular, we can write the parametric functional (11.34) with a minimum norm stabilizer as follows:
2 Pλ (σ ) = Wd (A(ρ) − d)2G + λ Wm ρ − ρ apr G = min, (11.35) where Wd and Wm are positively determined linear data and model parameter weighting operators, respectively; λ is the regularization parameter; and ρ apr is an a priori model. Note
268
Chapter 11
that, below, we will use the weighting operators Wd and Wm described by multiplication of the observed data d by a positive scalar function wd (ri ), and by multiplication of the anomalous density ρ(r) by a positive scalar function, wm (r). In the framework of the Newton method one finds a solution to minimization problem (11.35) in one step: ρ 1 = ρ 0 + δρ.
(11.36)
The optimum step, δρ, satisfies the following normal equation (see Section 5.2.2): (Hm0 + λWm Wm )δρ = −l0 − λWm Wm (ρ 0 − ρ apr ),
(11.37)
where Hm0 is a Hessian operator and l0 is a gradient vector at the initial iteration σ 0 . Following Section 5.2.3, we can find an approximate solution to the normal (11.37) if we
assume that the regularization parameter λ is large enough to neglect the term Hm0 δρ with respect to the term λWm Wm δρ in (11.37). Then we obtain λWm Wm δρ ≈ −l0 − λWm Wm (ρ 0 − ρσ apr ). −1
Applying the inverse weighting operator Wm Wm to both sides of the last equation, we find:
−1 l0 − (ρ 0 − ρ apr ). δσ ≈ −λ−1 Wm Wm
(11.38)
Note that the coefficient λ−1 can be treated as a scaling factor that can be found by minimizing the misfit between the observed and predicted data. Equation (11.38) plays an important role in rapid imaging since it provides an approximate solution to the inverse problem:
−1 ρ 1 = ρ 0 + δρ ≈ ρ apr − λ−1 Wm Wm l0 , (11.39) which requires only gradient direction calculations. In the special case where the initial and a priori models of anomalous density are equal to zero (i.e., we start the inversion with no background density model), we arrive at the following extremely simple and important equation:
−1 ρ 1 ≈ −λ−1 Wm Wm l0 , (11.40) which serves as a basis for migration imaging.
11.3.2 Migration of Gravity and Gravity Tensor Fields and 3-D Density Imaging Following the principles introduced in Chapter 10 for 2-D migration, the migration gravity field, gm α (r), is introduced as a result of application of the adjoint gravity operator, Aα , to the observed component of the gravity field: gm α (r) = Aα gα ,
(11.41)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data where the adjoint operator Agα for the gravity problem is equal to: f (r) K (r − r) ds. A α (f ) = − r|3 α |r S
269
(11.42)
From the physical point of view, the migration field is obtained by moving the sources of the observed gravity field above the observational surface. Nevertheless, the migration field contains some remnant information about the original sources of the gravity anomaly. That is why it can be used in imaging the sources of the gravity field. In a similar way, we can introduce a migration gravity tensor field gm αβ (r) and use the following notations for the components of this tensor field: gm αβ (r) = Aαβ gαβ ,
(11.43)
gm (r) = A g ,
(11.44)
where the adjoint operators, A αβ and A (f ), applied to some function f (r), are given by the formulas f (r) K (r − r) ds, (11.45) Aαβ (f ) = 3 αβ S |r − r| f (r) K (r − r) ds. (11.46) A (f ) = 3 |r − r| S We should note, however, that the direct migration of the observed gravity and/or gravity tensor fields does not produce an adequate image of the subsurface density distribution because the migration fields rapidly attenuate with the depth, as one can see from expressions (11.42), (11.45), and (11.46). In order to image the sources of the gravity fields at their correct location, one should apply an appropriate spatial weighting operator to the migration fields. This weighting operator is constructed based on the integrated sensitivity of the data to the density. We shall now apply (11.40) for rapid imaging. Substituting (11.18) for the direction of steepest ascent into (11.40), and taking into account (11.41), one can find a distribution of the density of the gravity field sources, described by the following expression:
−1 l0 , (11.47) ρ 1 ≈ −λ−1 Wm Wm
−1 m Aα gα = kα w−2 (11.48) ραm (r) = λ−1 Wm Wm α (z) gα (r), where unknown coefficient kα = λ−1 can be determined by a linear line search (Zhdanov, 2002) according to the following: w 2 A gα α M kα = (11.49) , Aw Aw gα 2 α α D −1 Aw α = Aα Wα ,
(11.50)
270
Chapter 11
and the linear weighting operator Wm = Wα is selected as a linear operator of the multiplication of the density ρ by a function, wα , equal to the square root of the integrated sensitivity of the complex intensity of the gravity field, Sα : wα = Sα . (11.51)
11.3.3 Integrated Sensitivity of 3-D Gravity Fields The integrated sensitivity, in accordance with the definition, is calculated as: obs δg α G , Sα = δρ
(11.52)
where δgobs α is the perturbation of the gravity field resulting from a local perturbation of the density, δρ(r) = ρ(r) dv, within a differential element of volume dv, located at the point r = (x, y, z) of the lower half-plane (z > 0): obs δgobs α = δgα (r ) = γ
ρ(r) Kα (r − r ) dv. |r − r |3
(11.53)
Substituting (11.53) into (11.52), we find 1 Sα = ρ(r) dv
=γ S
S
2 ds δgobs α (r )
Kα2 (r − r ) ds , |r − r |6
(11.54)
where S is the surface of observations of the gravity field. In particular, if the profile of observations coincides with the horizontal plane, z = 0, the definite integral in (11.54) can be evaluated as: ∞ K 2 (r) α Sα (z) = γ dx dy. (11.55) 6 |r| −∞ Using a polar system of coordinates, then x = R cos ϕ,
y = R sin ϕ,
ds = dx dy = R dR dϕ,
3 3 |r|6 = x2 + y2 + z2 = R2 + z2 , Kα (r − r) = α − α,
α = x, y, z,
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data and we can express (11.54) for α = x, y, z, as follows: ∞ 4π 2 R cos2 ϕ Sx (z) = γ 3 R dR dϕ, 0 0 R2 + z2 Sy (z) = γ
∞ 4π
0
Sz (z) = γ |z|
∞ 4π 0
Taking into account that: 4π 2 cos ϕ dϕ = 0
4π
(11.56)
R2 sin2 ϕ 3 R dR dϕ, R2 + z2
0
0
(11.57)
R dR dϕ 3 . R2 + z2
(11.58)
sin ϕ dϕ = 2π, 2
4π
and
0
271
dϕ = 4π,
(11.59)
0
we obtain the following: Sx (z) = Sy (z) = γ 2π
∞
0
Sz (z) = γ |z| 4π
4π
0
R2 3 R dR, R2 + z2
R dR R2 + z2
3 .
Calculating the tabulated integrals in (11.60) and (11.61), we finally find: √ 1 π 1 , Sz (z) = γ π . Sx (z) = Sy (z) = γ |z| 2 |z|
(11.60)
(11.61)
(11.62)
Thus, the integrated sensitivity, Sα , can be computed from the following equations: Sα = c α
1 , |z|
z < 0, α = x, y, z,
(11.63)
and cα are the corresponding constants for different components equal to: √ π cx = cy = γ , cz = γ π. 2 g
Equation (11.48) is called a migration density, ρm (ζ ): m ραm (r) = kα w−2 αg (z) gα (r),
(11.64)
272
Chapter 11
and is proportional to the weighted migration field, gm α: kα |z| gm ραm (r) = α (r), cα where
gm α (r) =
S
gα (r ) Kα (r − r ) ds . |r − r |3
(11.65)
(11.66)
Thus, the migration transform with spatial weighting provides a stable algorithm for calculating the migration density. In a similar way, we can introduce a migration density based on the gravity tensor migration: m m (r) = kαβ w−2 ραβ αβ (z) gαβ (r),
where T kαβ
2 w Aαβ gαβ M = , w W 2 Aαβ Aαβ gαβ D
m m ρ (r) = γ k w−2 (z) g (r),
T k
w 2 A g M = , Aw AW g 2 D
(11.67)
(11.68)
where functions wαβ and w are equal to the square root of the integrated sensitivity of the gravity tensor fields, Sαβ and S , respectively: wαβ = Sαβ , w = S . (11.69)
11.3.4 Integrated Sensitivity of 3-D Gravity Tensor Fields The integrated sensitivity of the tensor field components is calculated as: obs δgαβ G , S= δρ
(11.70)
where δgobs αβ is the perturbation of the corresponding component of the gravity tensor field resulting from a local perturbation of the density, δρ(r) = ρ(r) dv, within a differential element of volume dv, located at the point r = (x, y, z) of the lower half-plane (z > 0): ρ(r) obs Kαβ (r − r ) dv. (11.71) δgobs αβ = δgαβ (r ) = γ |r − r |3 Substituting (11.71) into (11.70), we find the following: 2
1 ) (r ds δgobs Sαβ = αβ ρ(r) dv S K 2 (r − r ) αβ =γ ds , 6 S |r − r |
(11.72)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
273
where S is the surface of observations of the gravity tensor field. In particular, if the profile of observations coincides with the horizontal plane, z = 0, the definite integral in (11.72) can be evaluated as: ∞ K 2 (r) αβ Sαβ (z) = γ dx dy. (11.73) 6 −∞ |r| Using a polar system of coordinates, then: x = R cos ϕ,
y = R sin ϕ,
ds = dx dy = R dR dϕ,
3 3 |r|6 = x2 + y2 + z2 = R2 + z2 , Kαβ (r) =
αβ 3 |r| 2 , α = β
3 (α) − 1, α = β |r|2 2
,
α, β = x, y, z.
We can express (11.72) for α, β = x, y, z, as: 2 ∞
4π 2z2 − R2 R dR dϕ, Szz (z) = γ 5 0 0 R2 + z2 ∞ 4π R2 cos2 ϕ dϕ, Sxz (z) = Syz (z) = γ 3 |z| 5 R dR 2 2 0 0 R +z 2 ∞ 4π 1 R2 cos2 ϕ −1 Sxx (z) = Sxx (z) = γ 3 2 3 R dR dϕ. 2 2 R +z 0 0 R + z2 Taking into account that 4π 3 cos4 ϕ dϕ = π, 2 0
4π 0
cos ϕ dϕ = 2π, 2
4π
dϕ = 4π,
(11.74)
(11.75)
(11.76)
(11.77)
0
we obtain the following:
Szz (z) = γ 4π
2 2z2 − R2 5 R dR, 0 R2 + z2 ∞ R2 Sxz (z) = Syz (z) = γ 3 |z| 2π 5 R dR, 0 R2 + z2 ∞ R4 1 Sxx (z) = Syy (z) = 3γ 4π 2 3 R dR. 0 R2 + z2 R2 + z2 ∞
(11.78)
(11.79)
(11.80)
274
Chapter 11
Calculating the tabulated integrals in (11.78) to (11.80), we finally find: √ 3π Szz (z) = γ , z < 0, 2z2 √ 3π , z < 0, Sxz (z) = Syz (z) = γ 2z2 √ 3 π , z < 0. Sxx (z) = Syy (z) = γ 4z2
(11.81) (11.82) (11.83)
Thus, the integrated sensitivity of the gravity tensor field is calculated from the following equations: 1 1 Sαβ = cαβ 2 , S = c 2 , (11.84) z z where cαβ are the corresponding constants for different components equal to: √ √ 3π 3 π czz = czx = czy = γ , cxx = cyy = γ . 2 4 Expression (11.67) is called a tensor field migration density. It is proportional to the magnitude of the weighted tensor migration field gm αβ . Thus, migration transformation provides a stable algorithm for calculating migration density. Substituting (11.84) for the weighting function wT back into (11.69) and (11.67), we find that: kαβ 2 m k 2 m m m (r) = z gαβ (r), ρ (r) = z g (r), (11.85) ραβ cαβ c where
gαβ (r ) K (r − r ) ds , |3 αβ |r − r S g (r ) gm (r) = K (r − r ) ds . |3 |r − r S
gm αβ (r)
=
(11.86) (11.87)
Finally, we can consider joint migration of several components of the gravity tensor. For example, we can jointly migrate gzz , gxx , gzx , and the g components and find the corresponding migration density as per the following: kxx m kzx m k m m 2 kzz m ραβ (r) = z g (r) + g (r) + g (r) + g (r) . (11.88) czz zz cxx xx czx zx c Note that (11.88) can be simplified to: m m m m ραβ (r) = z2 azz gm zz (r) + axx gxx (r) + azx gzx (r) + a g (r) ,
(11.89)
where azz , axx , azx , and a can be treated as the weights of the corresponding migration fields in the density model, which can be empirically determined from model studies.
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
275
Thus, for gravity fields and their gradients, we have shown that potential field migration is an integral transformation of the gravity field and/or gradients into 3-D density distributions. Potential field migration is very fast and stable, and the results are effectively equivalent to those obtained from 3-D regularized inversion with smooth stabilizers (Zhdanov et al., 2010, 2011).
11.4 Migration of Total Magnetic Intensity Data The earth’s magnetic field is the vector sum of contributions from two main sources—a background field due to the dynamo effect of the earth’s liquid core, and anomalous fields due to magnetic rocks and minerals above the Curie isotherm. Magnetic vector data, measured by orthogonal fluxgate magnetometers, are dominated by the earth’s background field, and are thus very sensitive to instrument orientation. The development of reliable and low cost optically pumped magnetometers in the 1960s enabled direct measurement of the total magnetic intensity (TMI) regardless of instrument orientation. It is now routine practice that every airborne geophysical survey produces TMI data as a standard deliverable. Relative to the millions of line kilometers of TMI data acquired each year, 3-D inversions are rarely performed. This is a reflection of the limited capability of existing 3-D inversion software to invert entire surveys to 3-D earth models with sufficient resolution in sufficient time so as to affect exploration decisions. Most interpretations are simply based on picking lineaments from the maps of the TMI’s first vertical derivative. Structural interpretations are usually based on some kind of Euler deconvolution, eigenvector, wavelet, analytic signal, or depth-from-extreme-points methods. While such methods may provide information about the sources of the TMI field, it is not immediately obvious how that information can be quantified in terms of a 3-D susceptibility model. To that end, Zhdanov et al. (2012) developed a method of rapid 3-D imaging of the TMI field data based on the principles of potential field migration discussed above.
11.4.1 Adjoint Operator for the Total Magnetic Intensity We measure TMI data on a surface S above a domain V that is filled by magnetic sources with the intensity of magnetization I(r). The problem is to determine the magnetic susceptibility, χ(r). In what follows, we adopt the common assumptions that there is no remanent magnetization, that the self-demagnetization effect is negligible, and that the magnetic susceptibility is isotropic. Under such assumptions, the intensity of magnetization is linearly related to an inducing magnetic field, H0 (r), through the magnetic susceptibility: I(r) = χ(r)H0 (r), where r is the radius vector of a point within the volume V.
(11.90)
276
Chapter 11
It is well known that the anomalous scalar TMI data T generated by the magnetic sources within the volume V can be represented by the linear operator equation (Li and Oldenburg, 1996; Portniaguine and Zhdanov, 2002): χ(r) T(r ) = A (χ) = H0 K(r − r)dv, r ∈ / V, (11.91) 3 |r | − r V where A (χ) denotes the forward modeling operator, H 0 is the magnitude of the inducing field, and l is a unit vector in the direction of magnetization: H0 (r) = H 0 l(r), and K is the TMI kernel:
2 3 l · r −r K(r − r) = − 1. (11.92) |r −r|2 We introduce a real Hilbert space D of the TMI data with the following metric: f (r )g(r ) ds , f , g ∈ D, (f , g)D =
(11.93)
S
and also introduce a real Hilbert space M of models (i.e., magnetic susceptibility, χ) with the following metric: χ(r)η(r) dv, χ, η ∈ M. (11.94) (χ, η)M = V
According to definition, the TMI adjoint operator, A , satisfies to the following equation:
(11.95) (A (χ) , f )D = χ, A (f ) M . Let us write the left-hand side of this equation in explicit form: A (χ) f (r ) ds (A (χ) , f )D = S
χ(r) K(r − r) dv f (r) ds 3 − r| S V f (r ) χ(r ) H0 K(r − r)f (r) ds dv = 3 D S |r − r| f (r ) = χ, H0 K(r − r) ds . 3 S |r − r| M = H0
|r
By comparing (11.95) and (11.96), we can see that: f (r ) K(r − r) ds = 0, χ, A (f ) − H0 3 S |r − r| M
(11.96)
(11.97)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data where “ ” denotes the adjoint operator. Equation (11.97) holds for any χ, e.g.: f (r ) χ = A (f ) − H0 K(r − r) ds. 3 S |r − r| Substituting (11.98) into (11.97), we obtain: 2 f (r ) A (f ) − H0 = 0. K(r − r) ds 3 S |r − r| M
277
(11.98)
(11.99)
From (11.99), we find that the adjoint operator, A , applied to any function f (r ), is given by the following expression: f (r ) A (f ) = H0 K(r − r) ds. (11.100) 3 S |r − r| We note that, the TMI adjoint operator produces a function that is analytical everywhere in the subsurface. The significance of this will become apparent momentarily.
11.4.2 Migration of the Total Magnetic Intensity We have established in the previous sections that migration is the action of the adjoint operator on the observed data. It follows that the migration TMI field, T m (r), is introduced as the action of the TMI adjoint operator, A , on the measured TMI field, T: T m (r) = A T(r ).
(11.101)
From the physical point of view, the migration TMI field is obtained by moving the sources of the observed TMI field above the surface, S. Nevertheless, the migration TMI field contains information about the original sources of the TMI field, which is why it can be used for imaging. We note, however, that direct migration of the measured TMI field does not produce an adequate image of the susceptibility distribution because the migration fields rapidly attenuate with the depth. In order to image the sources of the TMI field at their correct locations, one should apply an appropriate spatial weighting operator to the migration TMI field. Similar to 3-D gravity migration, the migration magnetic susceptibility can be computed from:
−1 χ m (r) = k W W A T = kw−2 (z) T m (r). (11.102) In the last formula, the unknown coefficient k can be determined by a linear line search according to the following: k=
Aw T2M Aw Aw T2D
,
(11.103)
where Aw =AW −1 ,
(11.104)
278
Chapter 11
and the linear weighting operator, W, is selected as an operator of multiplication of the susceptibility χ by a depth weighting function, w, which is equal to the square root of the integrated sensitivity, S (z), of the TMI data: δT obs D w(z) = (11.105) = S (z), δχ
11.4.3 Integrated Sensitivity of the Total Magnetic Intensity The integrated sensitivity of the TMI is calculated as: δT obs D S= , δχ
(11.106)
where δT obs is the perturbation of the corresponding component of the TMI resulting from a local perturbation of the magnetic susceptibility, δχ (r) = χ(r) dv, within a differential element of volume dv, located at the point r = (x, y, z) of the lower half-space (z > 0): δT obs = δT obs (r ) = H0 Substituting (11.107) into (11.106), we obtain: S = H0 S
χ(r) K(r − r ) dv. |r − r |3
K 2 (r − r ) ds , |r − r |6
(11.107)
(11.108)
where S is the surface of observations of the TMI. In particular, if the surface of observations coincides with the horizontal plane, z = 0, the definite integral in (11.108) can be evaluated as follows: 2 2 ∞ 4π
lx R cos ϕ + ly R sin ϕ + lz z R dR dϕ S (z) = H0 (11.109) 3 −1 3 , 2 2 R +z 0 0 R2 + z2 where we use a polar system of coordinates (ϕ, R, z): x = R cos ϕ, Taking into account that 4π 3 cos4 ϕ dϕ = π, 2 0
y = R sin ϕ, 0
4π
ds = dx dy = R dR dϕ.
cos ϕ dϕ = 2π, 2
0
4π
dϕ = 4π,
(11.110)
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
279
we obtain the following expression:
S (z) = 6H0 π 0
∞
R5
5 dR. R2 + z2
(11.111)
Calculating the tabulated integral in (11.111), we finally find a simple analytical formula for integrated sensitivity of the TMI field: √ 3 π (11.112) S (z) = H0 2 , z < 0. 4z Substituting (11.101), (11.105), and (11.112) into (11.102), we find that: T(r ) 4z2 k m K(r − r ) ds . (11.113) χ (r) = √ − r|3 |r 3H02 π S The migration magnetic susceptibility (11.113) is proportional to the magnitude of the weighted migration field, T m (r), which is analytical everywhere in the subsurface, implying that migration is a well posed and stable transform. It can be shown that, if the surface S coincides with the horizontal axis, z = 0, migration is mathematically equivalent to a form of downward continuation of a function analytical everywhere in the lower half-space. It should be noted that the downward continuation of the measured TMI field and the migration TMI field are significantly different. The migration TMI field has singular points in the upper half-space; therefore its downward continuation is a well posed and stable transform. The downward continuation of the measured TMI field has singular points in the lower half-space associated with its sources, so its downward continuation can only be extended down to these singularities, making it an ill-posed and unstable transform (e.g., Strakhov, 1970; Zhdanov, 1988).
11.4.4 Model Study To investigate the performance of 3-D TMI migration, we considered two discrete bodies with a susceptibility of 0.05 buried about 110 m below variable topography (Figure 11.1a). The inducing field had an inclination of 75 degrees and a declination of 25 degrees. The synthetic TMI data were computed on a 20 m regular grid draped over a curvilinear surface covering an area of one square kilometer. Figure 11.1b shows an example of a profile of the TMI data over both bodies. We applied 3-D migration to the entire synthetic TMI dataset, and calculated the migration magnetic susceptibility according to (11.112). As an example, Figure 11.1c shows a vertical cross section through the 3-D migration magnetic susceptibility model. The correct locations of the two bodies can be clearly determined. In order to demonstrate the robustness of migration to noise, we contaminated the data with 20% random Gaussian noise, as shown in Figure 11.1d. Figure 11.1e shows the same vertical cross section of results for the 3D migration of the noisy data. As expected, 3-D migration produced a very robust image of the susceptibility distribution.
(a)
(b)
(c)
(d)
(e) Figure 11.1 (a) Perspective of 3-D model consisting of two blocks of 0.05 susceptibility. Observation sites (dark gray dots) are located across a variable topography. (b) Synthetic TMI data along the profile y = 0 m. (c) Vertical cross section of susceptibility along y = 0 from 3-D migration of synthetic TMI data. (d) Synthetic TMI data with 20% random Gaussian noise along the profile y = 0 m. (e) Vertical cross section of susceptibility along y = 0 from 3-D migration of noisy synthetic TMI data.
Migration of 3-D Gravity, Gravity Tensor, and Total Magnetic Intensity Data
281
References Berkhout, A.J., 1980. Seismic Migration. Elsevier, Amsterdam. Claerbout, J.F., 1985. Imaging the Earth’s Interior. Blackwell Scientific Publications, Oxford. ˇ Cuma, M., Wilson, G.A., Zhdanov, M.S., 2012. Large-scale 3D inversion of potential field data. Geophys. Prospect. 60 (6), 1186-1199. ˇ Cuma, M., Zhdanov, M.S., 2014. Massively parallel regularized 3D inversion of potential fields on CPUs and GPUs. Comput. Geosci. 62, 80-87. Fedi, M., 2007. DEXP: a fast method to determine the depth to the sources of potential fields. Geophysics 72, L1-L11. Li, Y., Oldenburg, D.W., 1996. 3-D inversion of magnetic data. Geophysics 61, 394-408. Li, Y., 2001. 3D inversion of gravity gradiometer data. In: 71st Annual International Meeting, SEG, Expanded Abstracts, pp. 1470-1473. Portniaguine, O., Zhdanov, M.S., 2002. 3-D magnetic inversion with data compression and image focusing. Geophysics 67, 1532-1541. Schneider, W.A., 1978. Integral formulation for migration in two and three dimensions. Geophysics 43 (2), 49-76. Strakhov, V.N., 1970. Some aspects of the plane inverse problem of magnetic potential. Izv. Acad. Nauk SSSR Fiz. Zemli 9, 31-41 (in Russian). Tarantola, A., 1987. Inverse Problem Theory. Elsevier, Amsterdam. Wan, L., Zhdanov, M.S., 2008. Focusing inversion of marine full-tensor gradiometry data in offshore geophysical exploration. In: 76th Annual International Meeting, SEG, Expanded Abstracts, pp. 751-754. Zhdanov, M.S., 1988. Integral Transforms in Geophysics. Springer-Verlag, Berlin. Zhdanov, M.S., 2002. Geophysical Inverse Theory and Regularization Problems. Elsevier, Amsterdam. Zhdanov, M.S., Ellis, R.G., Mukherjee, S., 2004. Regularized focusing inversion of 3D gravity tensor data. Geophysics 69, 925-937. Zhdanov, M.S., 2009. Geophysical Electromagnetic Theory and Methods. Elsevier, Amsterdam. Zhdanov, M.S., Liu, X., Wilson, G., 2010. Potential field migration for rapid 3D imaging of gravity gradiometry surveys. First Break 28 (11), 47-51. Zhdanov, M.S., Liu, X., Wilson, G.A., Wan, L., 2011. Potential field migration for rapid imaging of gravity gradiometry data. Geophys. Prospect. 59, 1052-1071. Zhdanov, M.S., Liu, X., Wilson, G.A., Wan, L., 2012. 3D migration for rapid imaging of total-magnetic-intensity data. Geophysics 77 (2), J1-J5.