# Geophysics

- Society of Exploration Geophysicists

## Abstract

We extend the cross-gradient methodology for joint inversion to three-dimensional environments and introduce a solution procedure based on a statistical formulation and equality constraints for structural similarity resemblance. We apply the proposed solution to the joint 3D inversion of gravity and magnetic data and gauge the advantages of this new formulation on test and field-data experiments. Combining singular-value decomposition (SVD) and other conventional regularizing constraints, we determine 3D distributions of the density and magnetization with enhanced structural similarity. The algorithm reduces some misleading features of the models, which are introduced commonly by conventional separate inversions of gravity and magnetic data, and facilitates an integrated interpretation of the models.

## INTRODUCTION

Our understanding of complex subsurface processes occurring in many geologic environments, from the combination of gravity, magnetic, and other geophysical data, demands a detailed analysis of the coupled distribution of their physical properties (e.g., Gallardo and Meju, 2003; Linde et al., 2006; Kowalsky et al. 2006). Although for some geologic environments this distribution can be represented by 2D or even 1D structures, a detailed analysis relies on an accurate determination of 3D models.

Unfortunately, the abundance of subsurface features inherent in a 3D model contrasts with the limited resolution capability of much geophysical data, such as gravity and magnetic data. This leads to inaccurate and nonunique models. To compensate for this deficiency in information, we can simplify our models by making them smooth or piecewise homogeneous, or by adding complementary information from, for instance, borehole logs or surface geology. Alternatively, we might expect that combined analysis of two apparently disparate geophysical data sets should be useful not only for lithology discrimination and classification (Bosch, 1999; Ezzedine et al., 1999; Chen and Rubin, 2003; Bedrosian et al., 2007) but also to increase accuracy and resolution of 3D geophysical models when they are estimated jointly (Zeyen and Pous, 1993; Afnimar et al., 2002; Gallardo-Delgado et al., 2003; Linde et al., 2006; Tryggvason and Linde, 2006; Linde et al., 2008).

By relying on direct or indirect parameter interdependence, joint inversion can restrict the model space successfully to only those models that satisfy some cross-linked characteristics. Link selection is a key issue that has led to diverse methodologies for joint inversion (e.g., Zhang and Morgan, 1996; Haber and Oldenburg, 1997; Bosch and McGaughey, 2001; Gallardo and Meju, 2003; Musil et al., 2003; Saunders et al., 2005; Chen et al., 2006; Pilkington, 2006). An emerging methodology relies on the idea that physical properties tend to change at the same location and focuses on the search for structural similarities (Zhang and Morgan, 1996; Haber and Oldenburg, 1997; Gallardo and Meju, 2003, 2004; Saunders et al., 2005).

Gallardo and Meju (2003, 2004) propose a technique to invert two data sets jointly based on cross-gradient constraints to determine 2D resistivity and velocity models where collocated changes in the parameters occur in the same direction regardless of their magnitudes. Tryggvason and Linde (2006) apply the cross-gradient function to 3D models to search for P- and S-wave velocity structure. Linde et al. (2006) apply it to improve hydrogeologic studies using crosshole electrical resistance and ground-penetrating-radar traveltime data. Gallardo and Meju (2003, 2004, 2007) and Gallardo (2007) incorporate the cross-gradient function as an equality constraint and solve the minimization problem using Lagrange multipliers, whereas Tryggvason and Linde (2006) and Linde et al. (2006, 2008) incorporate cross gradients in the objective function and minimize them.

Here we develop a 3D joint-inversion technique using cross-gradient constraint and apply it to gravity and magnetic data. We start by analyzing information provided by the full 3D cross-gradient function and its implications as a constraint in a joint-inversion formulation. Then, we extend the cross-gradient technique to 3D structures based on equality constraints and solve the 3D joint-inverse problem using singular-value decomposition (SVD) and the generalized nonlinear least-squares framework developed by Tarantola and Valette (1982). We test our formulation on illustrative synthetic experiments using gravity and magnetic data, and compare the results with those obtained through separate inversion of each data set. In addition, we apply our algorithm on a field-data example in Todos Santos Bay, Mexico.

## 3D CROSS-GRADIENT JOINT INVERSION FOMULATION AND SOLUTION

### 3D cross-gradient constraint

The 3D cross-gradient function, defined by Gallardo and Meju (2003), is given by

For 3D models, *x*-, *y*-, and *z*-components of the cross-gradient function are given by
*x*- and *z*-directions for both geophysical models.

To understand the information provided by each component of the cross-gradient function and how they can be combined to quantify structural similarity in 3D models, we constructed two illustrative heterogeneous models, an ellipsoid and a spheroid, illustrated in Figure 1. It is clear that, although the marked vertical slice suggests a similar 2D structure, the structure depicted in orthogonal planes is different. In this case, *x*-axis. Full 3D structural similarity should be observed when the three components are equal to zero as it occurs in the central vertical section of these models.

### Formulation of the generalized 3D cross-gradient joint inversion

As described in Gallardo et al. (2005) for the 2D case, incorporation of cross-gradient constraints does not solve ambiguities found in geophysical problems. Thus, formulation of a stable solution algorithm requires additional information in the form of other regularizing constraints. Although the formulation described by Gallardo and Meju (2003) seems extendable to a joint 3D inversion procedure, we prefer to base our inversion algorithm on the generalized least-squares formulation proposed by Tarantola and Valette (1982). Unlike the Lagrange multiplier solution developed by Gallardo and Meju (2003, 2004), we expect this formulation will provide the more robust statistical framework necessary to assess the joint 3D problem. In a broad sense, the approach of Tarantola and Valette (1982) can incorporate conventional regularizing constraints (e.g., Marquardt, 1970; de Groot-Hedlin and Constable, 1990) as a new random variable (i.e., additional nongeophysical information). However, incorporation of the cross-gradient constraint involves subtle differences.

Tarantola and Valette (1982) propose a nonlinear generalized least-squares problem that minimizes the objective function
**x** includes both parameters of subsurface **m** and their physical and geophysical characteristics **d** that form the physical system

Usually, nonlinear equation 5 represents a functional based on physical laws that relate the data to the parameters. It can be expressed generally as
**F** is the matrix of partial derivatives

To adapt this formulation to our joint-inversion problem, we redefine the vector of parameters as

Therefore, the a priori information described in equation 4 is given by the actual observed data

To construct the constraining equations 5, we considered geophysical relationships for both data types, smoothness operator **D** (applied to both models), and the cross-gradient constraint given by equation 2. This last constraint provides the correlation between the model parameters required to formulate a joint-inversion algorithm. All these relationships can be written in compact form as:

One difference from the original formulation of Tarantola and Valette (1982) is that **x**. This will lead to an alternative solution to those shown by Tarantola and Valette (1982). It constitutes an exact interdependence different from that followed by Tryggvason and Linde (2006) and Linde et al. (2006, 2008). To find this solution, we used equations 8 and 9 to develop an explicit expression for estimated model parameters **F**. These equations lead to the generalized nonlinear solution for cross-gradient joint inversion in Appendix A. To find an explicit solution to our application, we follow the fixed-point method suggested by Tarantola and Valette (1982) and propose the following iterative solution:
*n* is the number of the model parameters. A component of the

Matrices

Vector *k* iterations and

The Jacobian matrix associated with the cross-gradient function is given by

In compact notation, equation 14 can be written as

Analyzing the rank of matrices

The pseudoinverse of matrix **U** and **V** are orthonormal matrices with size *Λ* is the corresponding singular value matrix. These matrices can be partitioned as
*p* value such that

### Computational implementation

The subsurface is discretized into rectangular prisms of homogeneous properties which constitute model parameters

Corresponding expressions for discretized derivatives of the cross-gradient function are derived straightforwardly from equation 26. They provide the Jacobian matrix **B**. This matrix contains six nonzero elements on each row, given by
*n* rows of **B** in equation 21. Similar results can be obtained for the other components. In the same way, we define a centered seven-point discrete scheme of the Laplacian operator included in **G**. Geophysical forward-modeling solutions

Following previous experiences in cross-gradient-constrained problems in two dimensions (Gallardo and Meju, 2004; Gallardo et al., 2005), we select control parameters

Covariance matrices for a priori models

## 3D JOINT INVERSION OF GRAVITY AND MAGNETIC DATA

### Gravity and magnetic formulation

We implemented the developed algorithm to the joint inversion of gravity and magnetic data to search for the distribution of density and magnetization. We composed our subsurface model as an aggregate of rectangular prisms and computed its gravity and magnetic responses using equations developed by Bhaskara-Rao et al. (1990) and Bhattacharyya (1966), as written in Gallardo-Delgado et al. (2003). Following these formulations, the vertical component of the gravity attraction *γ* is the gravitational constant, *r* is the distance between the observation point *ρ* represents the mass density. Limits

Similarly, the magnetic field of a homogeneously magnetized 3D prism is given by
*J* is the total magnetization in the geomagnetic field direction *l*, *m*, and *n* are directive cosines of the geomagnetic field. The rest of the variables are the same as in equation 28. We assume there is no remanent magnetization and, therefore, that the magnetization vector is parallel to the geomagnetic field direction.

Following equations 28 and 29, we have a linear relationship between our model parameters *ρ* and *J* and our geophysical data

### Synthetic test model

To prove our 3D cross-gradient joint-inversion algorithm, we used *x*- and *y*-directions and

We calculated a total of 1681

### Separate versus joint inversion of test data

To gauge the main advantages of the joint 3D inversion formulation, we perform separate and joint inversion of gravity and magnetic data on two comparative inversion schemes:

#### Homogeneous smoothing

For this experiment, we selected a smoothing Laplacian operator with homogeneous weight for the whole model space. Thus, elements of the diagonal Laplacian covariance matrices are constant values, given by

#### Fixed bottom

Conditions imposed on this experiment are similar to homogenous smoothing but we incorporate some depth-weighting factor by fixing the bottom prisms of the models at

For all experiments, we followed similar procedures to select optimal factors. First, we set a priori and initial models as null homogeneous distributions of density and magnetization and assigned a priori covariance matrices with large values (*p*-value that coincides approximately with the number of parameters

Anomalies obtained after separate inversion for both experiments are shown in Figure 3b–f. The associated models obtained after separate and joint inversion are shown in Figure 4. It is noticeable that the models obtained after separate inversion using homogeneous smoothing are distorted with broad tails at depth (Figure 4a and b). This effect is common in the inversion of gravity and magnetic data based on different objective functions (cf., Li and Oldenburg, 1998; Portniaguine and Zhdanov, 2002). The selected large values of the covariance matrices of the a priori models and the noise added to the data yielded depth-extended models, which were controlled mainly by the smoothness term. However, the addition of other regularizing constraints and the cross-gradient constraint diminished this effect (Figure 4e through h). Note the poor approximation to hypothetical parameter values for all the experiments, which is less than 10% for the target. Although it is well known that this problem can be reduced by introducing a priori information into the inversion (for instance, through the a priori models), we assume that little information of this kind is usually available and concentrate on the effect of the cross-gradient constraint with some standard smoothness regularization.

By comparing density and magnetization models for separate and joint inversion for the two selected experiments shown in Figure 4, it is clear that all models obtained from joint inversion show several improvements. For instance, the density model is recovered significantly in its original position and the broad tails that appeared in separately estimated models using homogeneous smoothing (where there is no causative body) are reduced in their value and extent. Geometrical dispersion of the models is reduced, concentrating distribution of density and magnetization around the target.

Although data misfits account for geophysical support of the models, all models have rms values near one in this case. Structural similarity achieved by the models can be quantified by the norm of the cross-gradient vectors (see Figure 5). In the two separate inversion experiments (Figure 5a and b), regions with larger structural concordance correspond to homogeneous portions of either model or to regions with more parallel density and magnetization changes.

Magnitudes of the cross-gradient vectors for jointly inverted models (Figure 5c and d) are at least four orders smaller than in separately estimated models, indicating an increased structural similarity. In the illustrated sections, major differences in separately inverted models are concentrated above the cube, which is diminished after joint inversion.

Figure 6a shows the density and magnetization crossplot obtained from separate inversion for the homogeneous smoothing experiment. The largest cluster dispersion corresponds to cells on the bottom and covers the region of the target cube. However, this curve does not provide information about any trend toward real density value, which is largely underestimated (Figure 4a). Because of the good structural similarity achieved by joint inversion, we can observe a clearly aligned cluster in corresponding crossplots (Figure 6b) that suggests a direct proportionality between density and magnetization in all blocks. Under stronger data constraints, this might lead eventually to actual test values of

To show evidence of model validity and process stability, we plot data misfits (Figure 6c) and model convergence (Figure 6d) achieved for several iterative steps of the first experiment using Laplacian smoothness terms. The data misfit was calculated with the equations
*g* and *m* refer to gravity and magnetics, respectively, and

### Field example

To fully test the methodology in actual field data, we selected an area that covers the Todos Santos Bay, on the northwest coast of Mexico (

For the inversion process, we selected the region in Figure 7, rotated to optimize the block distribution and discretized with a grid of

For joint inversion of these field data, we selected

Figure 9 shows results obtained for the field experiment in the representative sections

## CONCLUSION

We have obtained an iterative formulation for cross-gradient joint inversion of two geophysical data sets associated with 3D models. The developed formulation gives solid support for model appraisal because it involves linear and nonlinear functional relationships between both data and parameters equally under a statistical framework. The computational solution found shows efficiency on regular-scale models, for which truncated singular values account for redundancy of 3D cross-gradient constraints. Although finding an optimal SVD cutoff level solved for this redundancy, it also allowed the exploration of structurally coupled 3D models, which strengthened our selection of geologically suitable models.

Results suggest that jointly inverted models show better definition of structural features and, in combination with conventional regularizing constraints, might partially overcome the inherent lack of depth resolution in gravity and magnetic data. Lateral resolution for potential data also facilitates the search for models with lateral geometrical similarity and the cross-gradient constraint seems to alleviate some problems with depth resolution. For our field experiment, the procedure could uncouple two basement units characterized by different magnetization and no significant density variation, conserving a cross-gradient measured similarity.

Although this work provides encouraging experiences in application of joint inversion of gravity and magnetic data on structurally constrained 3D models, we foresee that higher-resolution models still should demand additional efforts on large-scale inversion strategies.

## ACKNOWLEDGMENTS

The authors acknowledge a scholarship awarded to Emilia Fregoso by CONACYT. We appreciate the constructive comments made by X. Li, M. Pilkington, F. Guspí, C. Farquharson, J. Chen, and anonymous reviewers, which helped us to improve the readability of the manuscript and to enhance the methodology and its application to field data.

## APPENDIX A

### MATHEMATICAL FORMULATION FOR THE SOLUTION OF THE 3D CROSS-GRADIENT JOINT-INVERSION PROBLEM

To demonstrate equations 14 and 15, we begin by defining the partial derivatives matrix in equation 9 of relationships for the homogeneous equation system 13:**x**, given by equation 9. Matrices involved in A-1 are defined as
*n* is the number of model parameters.

Matrix **F** can be partitioned in the indicated submatrices. We label the block of identity matrices as **I** and the zero matrices in the inferior left-hand block as **O**. Matrix **G** is defined as
**G** contains sensitivity matrices for both models and matrices of smoothness operator derivatives.**t** with respect to both model parameters.

To find the corresponding expression for equation 8, we divide this equation into the following three elements:

For our application, the first element

To find an expression for

Solving

Identities A-4, A-5, and A-7 are substituted into equation 8 to obtain an expression for estimated parameters **d** and **m**, given by

The estimated parameters of interest from A-9 are

Substituting A-10 into A-11 and applying A20 to A22 from Tarantola and Valette (1982), we obtain equations 14 and 15.

- Received April 24, 2008.
- Revision received December 12, 2008.