J. Cent. South Univ. (2016) 23: 1767-1777

DOI: 10.1007/s11771-016-3230-x

Joint inversion of gravity and multiple components of tensor gravity data

LU Guang-yin(鲁光银)¹, CAO Shu-jin(曹书锦)^{1, 2, 3}, ZHU Zi-qiang(朱自强)¹

1. School of Geosciences and Info-physics, Central South University, Changsha 410083, China;

2. School of Civil Engineering, Hunan University of Science and Technology, Xiangtan 411201, China;

3. Hunan Provincial Key Laboratory of Shale Gas Resource Utilization (Hunan University of

Science and Technology), Xiangtan 411201, China

Central South University Press and Springer-Verlag Berlin Heidelberg 2016

Abstract: Geological structures often exhibit smooth characteristics away from sharp discontinuities. One aim of geophysical inversion is to recover information about the smooth structures as well as about the sharp discontinuities. Because no specific operator can provide a perfect sparse representation of complicated geological models, hyper-parameter regularization inversion based on the iterative split Bregman method was used to recover the features of both smooth and sharp geological structures. A novel preconditioned matrix was proposed, which counteracted the natural decay of the sensitivity matrix and its inverse matrix was calculated easily. Application of the algorithm to synthetic data produces density models that are good representations of the designed models. The results show that the algorithm proposed is feasible and effective.

Key words: hyper-parameter regularization; full gravity gradient tensor; preconditioned matrix; Occam’s inversion; focusing inversion

1 Introduction

The gravity gradient tensor (GGT) is the spatial rate of change of gravitational acceleration. GGT data can be measured by partial or full tensor gravity gradiometer (which is called the FTG here) using land- or marine-based, or air- or space-borne platforms. It can be used to target oil, gas, and mineral deposits accurately, even for relatively deep objects [1-2]. In combination with pre-stack depth migration of seismic data, GGT data provide a potent mapping capability in pre-salt areas of the deepwater Gulf of Mexico. Traditional inversions of gravity data have been performed to recover a 3-D distribution of density contrast, for example, Occam’s inversion by LI and OLDENBURG [3], focusing inversion by ZHDANOV et al [4], and minimum- structure inversion by FARQUHARSON [5]. However, the interpretation of GGT data is not as straightforward as the more familiar gravity data. For a given source, GGT often produces a pattern of anomalies that are more complex than a simple single gravity anomaly. To overcome this difficulty, VASCO and TAYLOR [6] used gravity gradient tensors subset (g_xx, g_yy, and g_zz) to determine the basement topography in southwestern Oklahoma. LU and QIAN [7] developed a local level-set method for inverting 3D tensor gravity data based on a weighted L1-regularization method, but the convergence of the level-set algorithm relied heavily on initializing shape parameters of inversion. However, in practice, geological structures often exhibit smooth properties away from sharp discontinuities. Thus, the inversion of GGT data requires a suitable regularizer in two/ multiscale transformations. There is no specific sparse operator that can provide a perfect sparse representation of complicated geological models. The algorithm by O’BRIEN et al [8] combined seismic wave equation depth imaging and GGT data inversion, and was found to resolve the base of salt in the Gulf of Mexico much better than Kirchhoff depth imaging alone. The algorithm by ZHANG et al [9] integrated the ideas of the Occam’s inversion and sharp boundary inversion for 2-D magnetotelluric inversion using the minimum gradient support operator and conjugate gradient method. In fact, the convergence of the traditional sequential method is very slow when convex function p^a (m) contains hyper- parameter regularization terms. WRIGHT [10] proposed a primal-dual interior-point method for solving hyper-parameter optimization. However, the drawback is that this algorithm typically converges more slowly than traditional algorithms because the rate of convergence is linear. GOLDSTEIN and OSHER [11] proposed a very efficient split Bregman method, which can solve a very broad class of L1-regularized problems. JIA et al [12] produced a rigorous proof for the convergence of the split Bregman method. GHOLAMI and SIAHKOOHI [13] presented a novel and easy controllable joint sparsity regularization constraint algorithm for the inversion of seismic data, which benefited from two regularization terms. The inverse matrix of the preconditioned matrix in this algorithm must be calculated at each iteration or stored at the initialization of the inversion. However, it is possible that it could be so large that it may not be explicitly computed or stored.

In this work, joint inversion is analogous to the seismic inversion of GHOLAMI and SIAHKOOHI [13] to recover the properties of both smooth and sharp geological structures based on the split Bregman algorithm. A novel technique was proposed to avoid calculating or storing the preconditioned matrix. In all tests, the proposed algorithm successfully established more credible and meaningful geological models compared with other inverse methods.

2 Forward modelling of full gravity gradient tensor

In 1998, Bell Geospace took delivery of a full tensor gravity gradiometer built by Lockheed Martin for aiding navigation on a submarine used for oil and gas exploration. Since that time, the Bell Geospace’s full tensor gravity gradiometer has become commercially available. A full tensor gravity gradiometry survey has much higher bandwidth and delivers higher resolution images than conventional gravity surveys do. It can be used as a stand-alone service or in conjunction with existing seismic and other data, such as regional well logs, to provide cost-effective pictures of local geology that are much clearer. Current gravity gradiometers measure all gradients of three components of the gravity field forming a second order tensor, the elements of which are the second partial derivatives of the gravity potential f.

          (1)

The gravity potential is harmonic and well behaved outside source regions. Consequently, Γ is a full tensor gradient gravity matrix, which is symmetric and traceless, and thus only the five components indicated by the interior brackets are independent [14]. For all traditional synthetic tests, the designed model is divided into a set of prismatic cells using a 3-D orthogonal mesh and assuming a constant density within each cell (Fig. 1(a)). Given this discretization, taking g_zz as an example, which is produced by a prismatic cell with homogeneous density function at a given observation location (x, y, z), it can be written as [15]

          (2)

Here, the gravitational constant υ=6.672×10^{-8 m3}/(g·s²), and the scalars ξ_i, η_j and ζ_k are the coordinates of the opposing vertices of the prismatic cell.

where and

3 Joint inversion of gravity and multiple components of GGT data

According to Eq. (2), the sensitivity matrix for gravity and GGT data can be expressed in general matrix notation as

                 (3)

where G is the sensitivity matrix; the observation data d, such as g_α and g_αβ, consist of that unknown desired density model transformed by the sensitivity matrix G, where N is the number of observation points; M is the number of prismatic cells. For the inversion of gravity and single/multiple components of GGT data, the corresponding sensitivity matrix can be denoted by G=[G₁…G_i…G_n]^T. G_i can be represented by g_z, g_xx, g_xy, g_xz, g_yy, g_yz and g_zz, and d can be denoted by [d₁ … d_i … d_n 0]^T. The scalar n is the number of gravity gradient tensors subsets.

In general, the inverse problem (3) is regularized by the minimization of the Tikhonov functional to reduce the uncertainties associated with the estimation of geophysical variables. Then, this yields

                      (4)

where p^a (m) is a parametric objective function; is the data misfit functional; S(m) is the model objective functional; and μ is the regularization parameter. The Tikhonov regularizer isdefined as where D is the discrete derivative operator. The S^TK(m) attempts to penalize the coefficient of Dm due to the noise; thus, the large coefficients are associated with a strong penalty [13]. The L1 regularizer requires the solution to have a sparse representation with respect to a preselected basis or frame (e.g., wavelets) [16]. Geological structures often exhibit smooth properties away from sharp discontinuities (i.e., jumps in 1-D and edges in 2-D). The solutions of linear and nonlinear geophysical ill-posed problems require simultaneous sparse representation in two appropriate transformation domains. However, there is no specific operator available for this purpose [13]. To obtain a meaningful geophysical model, an attempt is made to recover the smooth and edge features using the hyper-regularization scheme [11, 13], which benefits from both the advantages of wavelets and L1 norm operators in the representation of the solution. Then, the hyper-regularization S(m) is redefined as

where 0<β<1 is a parameter that controls the relative weight of the two regularization terms, both γ and λ are regularization parameters,is a finite difference operator or a total variation operator, and W is a suitable sparse operator (e.g., Tikhonov regularization ||m||², weighted Tikhonov regularization ||Wm||₂, L1 wavelet). Equation (4) is reduced only with a wavelet-based sparsity regularizer (e.g., curvelet [16-17]) by assuming β=0. The wavelet transformation extracts multiscale information with too many small wavelet coefficients that might contain useful geological information. The sparsity regularizer applies a lesser penalty to the large coefficients, and vice versa. Then, the small wavelet coefficients are removed by the regularization parameter determined by noise level; however, this can lead easily to the fuzziness of the geological jumps caused by over-regularization. A piecewise-constant solution is obtained by taking β=1. To exploit the advantages of both the wavelet and finite difference regularizers, Eq. (4) can capture local jumps/discontinuities amongst the geological structures and variation of smoothness of the interior geological structures by setting β to a non-zero value.

4 Solution by split Bregman iterative algorithm

Following GOLDSTEIN and OSHER [11], SIMA et al [18], and PLONKA and MA [19], a novel inversion of gravity and single/multiple components of GGT data are introduced based on the split Bregman method, which separates the L2-norm and L1-norm terms from the cost function, to generate a sequence of unconstrained sub-problems that can be solved easily and more efficiently than by sequential inversion.

The hyper-regularization method is the minimizer of the cost function in Eq. (4) with multiple regularization terms. Therefore, Eq. (4) can be rewritten as

   (5)

where W is a class of multiscale L1 wavelet operator; is the total variation operator and its corresponding cost function is S^TV, which yields

According to the L1 norm [20], S^TV(m) becomes

where the subscripts of m_x, m_y and m_z denote the partial derivative with respect to the space variable x, y and z, respectively.

The local derivative of S^TV with respect to density m, yields

Following LORIS and VERHOEVEN [20], and RUDIN et al [21], a more stable form of is defined as

                          (6)

where

;

;

;

;

;

;

h_x, h_y and h_z are the length of side of prismatic cell. The subscripts of m_ijk denote the order number of prismatic cells with respect to the space variable x, y and z.

The preserving edge function f(a, b) is defined as

          (7)

where the sign function can be expressed as

Taking the partial derivative of m with respect to x_i (x₁=x, x₂=y, x₃=z), and letting the result equal zero, then if a is contrary to b, f(a, b) will be equal to zero and will be in contrast to Then, the discontinuities of the geological structure can be expressed mathematically as m_1jk=m_2jk, m_i1k=m_i2k, m_ij1=m_ij2, …, m_Mjk-m_(M-1)jk, m_iMk-m_i(M-1)k, m_ijM-m_ij(M-1).

To apply the split Bregman method, and Wm are replaced by p₁ and q₁ in Eq. (5), respectively [11]. The result is converted to the constrained problem:

                      (8)

where W is the L1 curvelet operator, which denotes a forward transform W^TW=I [19]; W^T is an inverse curvelet transform; I is an identity matrix.

We add penalty functional terms into Eq. (8) and convert it to an unconstrained problem. Therefore, Eq. (8) becomes:

To apply weakly enforced constraints for the split Bregman iteration, which are found by GOLDSTEIN and OSHER [11], the split formulation becomes:

       (9)

where is obtained by applying the difference operator to m and the temporary value b_W is obtained by the operator W.

Taking formula (9) with the derivative with respect to m and letting the result equal zero, the iterative equation can be expressed as

     (10)

where W_d is a data weighting matrix, in 3-D space, the subscripts of and denote the partial derivative with respect to the space variable x, y and z, the superscript T denotes the matrix transposition, and s^k+1 is denoted by

 (11)

where p₁ and q₁ are updated by the solved sub-programs at each iteration forand [22], w represents the multiscale wavelet coefficient.

The wavelet regularizer may yield many small wavelet coefficients of the geological model. Thus, this minimization is decomposed into sub-problems using the shrinkage technique with thresholding operation to update d_i and b_i. Simple iterative curvelet thresholding can follow MA [23]; sub-problems of the split Bregman method can follow GOLDSTEIN and OSHER [11]. In this case, the iterative form of the sub-problems can be written as

where the shrinkage operator is given by [19]

where subscripts ofrepresent frequency, orientation and position, respectively; T_σ is a soft shrinkage function defined by a fixed threshold σ>0:

Then, the temporary valuesand w are redefined as

5 Preconditioning technique

Because there are no specific gravity methods for geophysical exploration, except borehole gravity logging which can provide data for special depths, traditional gravity exploration only gathers 2-D surface data, meaning that there is no inherent depth resolution. Therefore, compared with other geophysical exploration methods, the inversion of surface gravity or GGT data is a more unstable and non-unique problem. For any kind of gravity inversion, regularization methods, efficient solvers, and physical constraints should be used to reduce the non-uniqueness and instability of the solutions in order to obtain geologically meaningful solutions.

The computational cost of the split Bregman method for joint inversion is largely dependent on Eq. (10). PILKINGTON [24] proposed a preconditioned matrix S by solving S^TSx=Ix based on the preconditioned conjugate gradient method to accelerate the 3-D magnetic inversion. Following this idea, the new preconditioned matrix can be defined as

where is a sparse matrix. In practice, Z₁ is a band matrix and its corresponding inverse matrix is too large to be constructed or stored; thus, it is only available for small-scale problems. For large-scale problems, the preconditioned matrix can be approximated with the diagonal matrix of G^TG. Thus, Z₁ is redefined as [25]:

where is a diagonal matrix. The inverse matrix of Z₂ must be calculated at each iteration of the preconditioned conjugate gradient method, which is a time-consuming problem with high demands on computational and physical memory usage because Z₂ is also a sparse band matrix. LI and OLDENBURG [3] introduced a depth weighting function into the objective function to counteract the natural decay of the sensitivity matrix, such that the inversion yields depth information, where z is the depth of the prismatic cell centroid, the decay index κ is dependent on the potential field, and the scalar z₀ is determined by the prismatic cell size and observation height. Thus, the preconditioner isThe natural decay of the sensitivity matrix can be approximated by adjusting the values of both z₀ and k. The computational cost of that inverts a sparse preconditioned matrix, which is much smaller than Z₁ and Z₂.

In practice, 1.5<κ<2.0 in the gravity inversion [3], 2.0

          (12)

Therefore, without consideration of the horizontal differences and boundary effects of discrete prismatic cells, Z_jj decays with depth. The scalar Z_jj attempts to penalize strong weighting on deep prismatic cells and weak weighting on shallow cells. This avoids the distribution that is an infinitesimally thin layer of mass just beneath the earth’s surface [27]. If γ and λ are assumed zero, then Z is equal to Z₂. The inverse matrix of Z and Z₃ needs to be calculated only once in the entire inversion process and the inverse matrix of Z₂ needs to be calculated at each iteration or stored in a full density matrix.

Integrating formula (12) into Eq. (5) as a preconditioned matrix yields:

       (13)

where and λ′=μ^-1λβ.

Similarly to formulas (8) and (10), the iterative form of Eq. (13) can be expressed as

       (14)

where

;

;

;

;

;

G_Z=W_dGZ, W_Z=ZW, and The inversion is initialized by setting the temporary values and to a zero value.

This split Bregman algorithm has several advantages over the traditional penalty function/ sequential methods for solving Eq. (14). First, the Split Bregman iteration converges very quickly when applied to certain types of objective functions, especially for problems where convex function P^a(m) contains an L1-regularization term. However, the sequential methods need several steps to solve each regularization term in Eq. (14). The second advantage is that the value of λ in Eq. (4) remains constant. Therefore, a value can be chosen for λ that minimizes the condition number of the sub-problems, resulting in faster convergence for iterative optimization methods. The split Bregman iteration also avoids the problems of numerical instabilities that occur as λ, which arises when using sequential methods [11].

For most wavelet functions in combination with the Haar operator, the relative weight β value ranges from 0.4 to 0.6; however, in the case of total variation operator, it ranges from 0.8 to 0.95 [13], and for the L1 wavelet operator, it ranges from 0.45 to 0.75.

6 Numerical experiments

The inversion codes are compiled using MATLAB pCode command under Windows 7. To verify the validity of the inversion algorithm, these three models have similar structures except their density distributions and depth. ZHDANOV et al [4] and ZHDANOV and MUKHERJEE [28] invert these three models by focusing inversion. The joint inversion based on the hyper-regularization scheme, Occam’s inversion, and focusing inversion are used to invert them in the following three tests.

6.1 Model of cubic bodies

Figure 1(a) is a representation of the designed model, which is divided into a set of prismatic cells using a 3-D orthogonal mesh and assuming a constant density within each prismatic cell. There are two cubic bodies buried at a depth of 300 m with the same density contrast (1.0 g/cm³) over the background, and the size of each body is 200 m×200 m×200 m. The designed model is discretized into 20×20×16=6400 prismatic cells in the x-, y-, and z-direction, respectively. The length of the prismatic cell edge is 50 m. The rectangular survey grid is divided into 400=20×20 observation points located on the earth’s surface. The observation height is 25 m and the space-sampling interval is 50 m in both the x- and y-directions. Then, we use our algorithm based on preconditioned matrix Z₃ inverted FTG data contaminated by 3% noise.

Fig. 1 Forward model of FTG data and comparison of effect of different inversion methods:

Figures 1(b)–(d) show the results of the Occam’s inversion, focusing inversion, and joint inversion of FTG data, respectively. The result of the Occam’s inversion (specifically, the flatness a_lx=1, a_ly=1 and a_lz=1 in the x-, y-, and z- directions, respectively and the amount of model smallness a_ls=0.0005. Readers are referred to LI and OLDENBURG [3] for further details) is smoothness or fuzziness. The focusing inversion is solved by reweighted regularized conjugate gradient method withminimum support operator where e=1.2×10^-6. It can be found a detailed explanation of by ZHDANOV et al [4]. It is easier to recover the sharp boundaries of the geophysical structures over the background with this method than with the Occam’s inversion. Figures 1(d) shows the result of joint inversion proposed in this work. Comparing it with the results of the Occam’s inversion published by LI and OLDENBURG [3] and shown in Fig. 1(b), the density boundaries can be represented more clearly and the positions of the cubic bodies are determined more accurately. Comparing the joint inversion proposed with the results of the focusing inversion published by ZHDANOV et al [4] and PORTNIAGUINE and ZHDANOV [29] and shown in Fig. 1(c), it reveals that the joint inversion has a relatively weak focusing effect and avoids inaccurate depth and obliqueness caused by over-regularization in the focusing inversion. The misfit of the results ranges from 0.03 to 0.1, which satisfies the accuracy of geophysical inversion (Fig. 1(e)). In particular, the parameter variation function reaches a minimum after the 8th iteration (Fig. 1(f)). Moreover, compared with other inverse methods, the joint inversion performs equally well in this case.

6.2 Influence of depth on joint inversion

Figure 2(a) shows the representation of the designed model. There are two cubic bodies, the densities of which are 2.0 and 1.0 g/cm³ for the left and right, respectively. A series of tests are set by moving Z_△ (moving step). The designed model is discretized into 40×40×10=16000 prismatic cells in the x-, y-, and z- direction, respectively. The length of the prismatic cell edge is 50 m. The rectangular survey grid is divided into 400=20×20 observation points located on the earth’s surface. The observation height is 25 m and the sampling interval is 50 m in the x- and y-direction. Then, we use our algorithm based on preconditioned matrix Z₃ inverted FTG data contaminated by 3% noise.

Fig. 2 Forward model of FTG data Z_△=200 (a) and comparison of effect of joint inversion with different depth sources Z_△=0 (b), Z_△=200 (c) and Z_△=400 (d)

Figures 2(b)-(d) shows the results of joint inversion based on preconditioned matrix Z₃ for Z_△=0, 200 and 400, respectively. Figure 2(b) shows that the final solutions for the two bodies have densities of 3.0 and 1.0 g/cm³, respectively. This is a good representation of the designed models. The numerical result in Fig. 2(c) shows that both geological bodies have almost equivalent anomalies because the left anomaly is much deeper than the other is; the potential value of the natural decay with depth is in accordance with Fig. 2(d).

Then, the proposed algorithm is used to invert the designed model (Z_△=400, contaminated by 3% noise) based on preconditioned matrices Z₂, Z₃ and Z, respectively.

Figure 3(a) shows the result of joint inversion based on preconditioned matrix Z₂, and compared with the result of joint inversion based on preconditioned matrix Z₃, shown in Fig. 2(d), the two results display the same effect. Figure 3(b) shows the result of joint inversion based on preconditioned matrix Z. Comparing this with the results of the same model shown in Figs. 2(d) and 3(a), it is illustrated that this method defines a position that is more accurate. Meanwhile, the densities of the shapes are in agreement with the models.

Figure 4(a) illustrates the NNZ curve of preconditioned matrix Z₂ as a dotted line and the other two curves are superpositioned. Figure 4(b) shows the efficiency of the calculation of the inverse matrix of the preconditioned matrices (NNZ represents number of non-zeros of the inverse matrix of preconditioned matrix). From these two figures, it can be determined that the inverse matrix of Z₂ could be so large that it may not be computed explicitly or stored, and that the processing of Z₃ is faster than that of the preconditioned matrix Z. However, it is difficult to choose between the joint inversion of gravity and GGT data. Compared with the others, the preconditioner Z is a more efficient preconditioned matrix, and the method proposed in this work is easier to calculate than the hyper-parameter regularization method of GHOLAMI and SIAHKOOHI [13].

6.3 Y-dyke model

A more complicated model consisting of two dykes dipping in opposite directions is presented, having different widths and lengths, but the same strike direction. This model is shown in one cross-section in Fig. 5(a). The two dykes have the same density contrast of 1.2 g/cm³ over the background. The designed model is discretized into 40×40×10=16000 prismatic cells in the x-, y-, and z-direction, respectively. The length of the prismatic cell edge is 50 m. The rectangular survey grid is divided into 400=20×20 observation points located at the earth’s surface. The observation height is 25 m and the sampling interval is 50 m in the x- and y-direction. Then, we use our algorithm based on formula (12), proposed in this work, to invert the FTG data contaminated by 3% noise.

Fig. 3 Comparison of effect of joint inversion of FTG data based on different preconditioned matrices:

Fig. 4 Comparison of computational cost of different preconditioned matrices with different designed model sizes:

Fig. 5 Forward model and comparison of effect of joint inversion on different combinations of gravity gradient tensor subsets:

Figure 5(b) shows the result of joint inversion of gravity gradient tensors subset (g_xx, g_yy and g_zz), the density boundaries are clearer and the right-hand anomaly is bigger than the Y-dyke in the designed model. However, the decay in depth of the result exhibits a slight difference over the designed model in the deeper geological structure. Figure 5(c) shows the result of joint inversion of FTG data. The deep anomaly is a good representation of the deep dyke in the designed model; however, the shallow anomaly with strictly vertical boundary caused by the deep anomaly is due to equivalence theory in Figs. 5(b) and 5(c). Figure 5(d) is the result of joint inversion of gravity and FTG data. Compared with the others, Fig. 5(d) shows very good agreement between our inversion results and the estimated depths to the mass centers obtained from the Occam’s inversion by LI and OLDENBURG [3] and the focusing inversion by ZHDANOV et al [4].

6.4 SEG/EAGE 3D salt model

FTG data are covered by the United States Munitions list 22.CFR121.1 and the export of FTG data must be licensed by the office of Defense Trade Controls, USA, Department of State [6]. Therefore, no FTG data are available. The SEG/EAGE salt model is designed as a velocity model for 3D seismic technology [30-31]. Figure 6(a) displays the perspective view of the density model by converting into density variations. To verify hyper-parameter regularization inversion, the salt body is assumed to have a constant density contrast of -0.2 g/cm³. The designed model is discretized into 676× 676×210=95964960 prismatic cells in the x-, y-, and z-directions, respectively. The length of the prismatic cell edge is 10 m. The rectangular survey grid is divided into 676×676=456976 observation points located on the earth’s surface. The observation height is 5 m and the space-sampling interval is 10 m in both the x- and y-directions. Then, we use our algorithm based on formula (12), proposed in this work, to invert the FTG data contaminated by 3% noise.

Compared with the designed model Fig. 6(a), the inversion results of gravity gradient tensor diagonal component g_xx and g_yy, as shown in Figs. 6(b) and (e), are presented as the deep structure of the model which has been successfully filled in with salt body, and the remainder of the model region has been reduced to zero.The density of inversion result ranges from -0.130 to 0.075. However, inversion results of gravity gradient tensor horizontal components g_xy, g_xz and g_yz, as shown in Figs. 6(c), (d), (f), exhibit a slight divergence over the designed model. Figure 6(g) shows the result of joint inversion of FTG data and a good representation of the designed model. The joint inversion of FTG data also avoids the distortion of density boundaries. Compared with the mass centers obtained from the binary’s inversion by KRAHENBUHL [31], Figs. 6(c), (d), (f) and (g) show very good agreement with the designed model.

Fig. 6 SEG/EAGE salt model (a), and comparison of effect of joint inversion on different combinations of gravity gradient tensor subsets g_xx (b), g_xy (c), g_xz (d), g_yy (e), g_yz (f), FTG data (g)

From the inversion results of these four examples above, we can determine that the joint inversion based on the hyper-regularization method represents not only the principal density boundaries more accurately, but it also avoids the distortion and over-regularization of density boundaries, Y-dyke and SEG/EAGE 3D salt model. This approach benefits from the advantages of both wavelet and finite difference operators in the representation of the solution. The accuracy of the results is also satisfactory.

7 Conclusions

A new regularization scheme for inverting gravity and GGT data is developed to recover a 3-D distributions of density contrast, and this scheme benefits from the advantages of two regularization terms in the representation of the desired model. By calculating three synthetic models and comparing them with previously published results from the Occam’s inversion and focusing inversion, the following conclusions are drawn.

1) The Occam’s inversion easily recovers smooth features and fuzzy geological edges.

2) The focusing inversion easily leads to very high density in small geological bodies caused by over- regularization.

3) The hyper-regularization scheme recovers the features of both smooth and sharp geological structures using the split Bregman algorithm. The numerical results present clearer density boundaries and positions that are more exact, and this inversion process is stable. Moreover, this new method could represent both the smoothness and sharpness of features. Crucial features of preconditioned matrix Z are that it counteracts the natural decay of the sensitivity matrix, its inverse matrix requires less memory, and its process is faster than the other methods.

References

[1] WANG Jing-bo, XIONG Sheng-qing, ZHOU Xi-hua, GUO Zhi-hong. The advances in the study of the airborne gravimetry system [J]. Geophysical and Geochemical Exploration, 2009, 33(4): 368-373. (in Chinese)

[2] O, GIMENEZ M, BRAITENBERG C, ANDRES F. Goce satellite derived gravity and gravity gradient corrected for topographic effect in the south central andes region [J]. Geophysical Journal International, 2012, 190(2): 941-959.

[3] LI Yao-guo, OLDENBURG D W. 3-D inversion of gravity data [J]. Geophysics, 1998, 63(1): 109-119.

[4] ZHDANOV M, ELLIS R, MUKHERJEE S. Three-dimensional regularized focusing inversion of gravity gradient tensor component data [J]. Geophysics, 2004, 69(4): 925-937.

[5] FARQUHARSON C G. Constructing piecewise-constant models in multi-dimensional minimum-structure inversions [J]. Geophysics, 2008, 73(1): K1-K9.

[6] VASCO D, TAYLOR C. Inversion of airborne gravity gradient data, southwestern Oklahoma [J]. Geophysics, 1991, 56(1): 90-101.

[7] LU Wang-tao, QIAN Jian-liang. A local level-set method for 3d inversion of gravity-gradient data [J]. Geophysics, 2015, 80(1): G35-G51.

[8] O’BRIEN J, RODRIGUEZ A, SIXTA D, DAVIES M, HOUGHTON P. Resolving the k-2 salt structure in the gulf of mexico an integrated approach using prestack depth imaging and full tensor gravity gradiometry [J]. The Leading Edge, 2005, 24(4): 404-409.

[9] ZHANG Luo-lei, YU Peng, WANG Jia-lin, WU Jian-sheng. Smoothest model and sharp boundary based two-dimensional magnetotelluric inversion [J]. Chinese Journal of Geophysics, 2009, 52(6): 1625-1632. (in Chinese)

[10] WRIGHT S J. Primal-dual interior-point methods[M]. Philadelphia: Society for Industrial and Applied Mathematics, 1997.

[11] GOLDSTEIN T, OSHER S. The split bregman method for L1-regularized problems [J]. SIAM Journal on Imaging Sciences, 2009, 2(2): 323-343.

[12] JIA Rong-qing, ZHAO Han-qing, ZHAO Wei. Convergence analysis of the bregman method for the variational model of image denoising [J]. Applied and Computational Harmonic Analysis, 2009, 27(3): 367-379.

[13] GHOLAMI A, SIAHKOOHI H. Regularization of linear and non-linear geophysical ill-posed problems with joint sparsity constraints [J]. Geophysical Journal International, 2010, 180(2): 871-882.

[14] KIRKENDALL B, LI Y G, OLDENBURG D. Imaging cargo containers using gravity gradiometry [J]. Geoscience and Remote Sensing, IEEE Transactions on, 2007, 45(6): 1786-1797.

[15] LI Xiong, CHOUTEAU M. Three-dimensional gravity modeling in all space [J]. Surveys in Geophysics, 1998, 19(4): 339-368.

[16] HERRMANN J, MOGHADDAM P, STOLK C C. Sparsity-and continuity-promoting seismic image recovery with curvelet frames [J]. Applied and Computational Harmonic Analysis, 2008, 24(2): 150-173.

[17] MA J, ANTONIADIS A, le DIMET F X. Curvelet-based snake for multiscale detection and tracking of geophysical fluids [J]. Geoscience and Remote Sensing, IEEE Transactions on, 2006, 44(12): 3626-3638.

[18] SIMA D M, van HUFFEL S, GOLUB G H. Regularized total least squares based on quadratic eigenvalue problem solvers [J]. BIT Numerical Mathematics, 2004, 44(4): 793-812.

[19] PLONKA G, MA Jian-wei. Curvelet-wavelet regularized split bregman iteration for compressed sensing [J]. International Journal of Wavelets, Multiresolution and Information Processing, 2011, 9(1): 79-110.

[20] LORIS I, VERHOEVEN C. Iterative algorithms for total variation-like reconstructions in seismic tomography [J]. GEM-International Journal on Geomathematics, 2012, 3(2): 179-208.

[21] RUDIN L I, OSHER S, FATEMI E. Nonlinear total variation based noise removal algorithms [J]. Physica D: Nonlinear Phenomena, 1992, 60(1): 259-268.

[22] DONOHO D L. Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition [J]. Applied and Computational Harmonic Analysis, 1995, 2(2): 101-126.

[23] MA Jian-wei. A single-pixel imaging system for remote sensing by two-step iterative curvelet thresholding [J]. Geoscience and Remote Sensing Letters, 2009, 6(4): 676-680.

[24] PILKINGTON M. 3-D magnetic imaging using conjugate gradients [J]. Geophysics, 1997, 62(4): 1132-1142.

[25] KOH K, KIM S J, BOYD S P. An interior-point method for large-scale L1-regularized logistic regression [J]. Journal of Machine Learning Research, 2007, 8(8): 1519-1555.

[26] PORTNIAGUINE O, ZHDANOV M S. 3-D magnetic inversion with data compression and image focusing [J]. Geophysics, 2002, 67(5): 1532-1541.

[27] PILKINGTON M. 3D magnetic data-space inversion with sparseness constraints [J]. Geophysics, 2008, 74(1): L7-L15.

[28] ZHDANOV M, MUKHERJEE S. Modeling and inversion of 3-d gravity tensor field[EB/OL].[2013-11-19]. http://www.cemi.utah. edu/appls/grav/tensgrav/tensor.html.

[29] PORTNIAGUINE O, ZHDANOV M. Focusing geophysical inversion images [J]. Geophysics, 1999, 64(3): 874-887.

[30] FERNANDO D, S. Adaptive learning 3d gravity inversion for salt-body imaging [J]. Geophysics, 2011, 76(3): I49-I57.

[31] KRAHENBUHL A. Binary inversion of gravity data for salt imaging [D]. Colorado, USA: Colorado School of Mines, 2005.

(Edited by FANG Jing-hua)

Foundation item: Projects(41174061, 41374120) supported by the National Natural Science Foundation of China

Received date: 2015-01-30; Accepted date: 2015-12-10

Corresponding author: LU Guang-yin, Professor, PhD; Tel: +86-13975894898; E-mail: luguangyin@csu.edu.cn