- Society of Exploration Geophysicists
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.
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 (1) where and are any two geophysical model parameters. Following the cross-gradient methodology, structural similarity is achieved when (2) i.e., when the collocated gradient vectors are parallel or one of them is null.
For 3D models, x-, y-, and z-components of the cross-gradient function are given by (3)and Each one of these components quantifies the structural similarities in planes normal to the component direction. For example, the component measures coupled spatial variations in the 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, vanishes over the whole model space but and still indicate the geometrical 3D differences of the 3D models. Figure 1c and d reflect these differences and provide clear evidence of the rotation origin of the models: radial cross-gradient components are zero and the remaining tangent component increases its magnitude as long as it spreads out from the rotation 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 (4) subject to the condition (5) where x includes both parameters of subsurface m and their physical and geophysical characteristics d that form the physical system (6) Variable is the a priori information with a covariance matrix . Quadratic function comprises all random variables, including data and prior model parameter information and contains physical and hypothetical exact constraining relationships.
Usually, nonlinear equation 5 represents a functional based on physical laws that relate the data to the parameters. It can be expressed generally as (7) The solution of the generalized minimization problem (equations 4 and 5), found by Tarantola and Valette (1982), under the Bayesian framework, is given by (8) where F is the matrix of partial derivatives (9) evaluated at point .
To adapt this formulation to our joint-inversion problem, we redefine the vector of parameters as (10) where and are the predicted observations for two geophysical data sets, and represent smoothness constraint values (based on either gradient or Laplacian operators) applied to the model parameters, and and are the two sets of model parameters.
Therefore, the a priori information described in equation 4 is given by the actual observed data and , a priori smoothness values and (in this case set to zero), and a priori model parameters and , which are all included in the vector (11) Accordingly, the covariance matrix is given by (12) Diagonal elements in correspond to covariance matrices and for each geophysical data set (which are assumed uncorrelated). and define weighting factors that regulate the level of smoothness required by the models. Finally, and are covariance matrices for a priori model parameters.
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: (13)In addition to the cross-gradient constraint, this formulation also permits incorporation of other physical or empirical relationships that might correlate both parameters directly. An example would be Poisson's relation for gravity and magnetic fields when they are produced by the same causative source, although we will not discuss it here.
One difference from the original formulation of Tarantola and Valette (1982) is that , , and are not included as random variables in 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 , by expanding the matrices involved in 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: (14) where (15) Furthermore, (16) includes sensitivity and smoothness matrices, for which a component of the matrix is given by , ; . The number of observed data related to model is , and n is the number of the model parameters. A component of the matrix is given by , ; .
Matrices and are given by (17) and (18)
Vector is the difference between observed and computed data. It includes the smoothness regularization terms (19) and contains components of the cross-gradient vector (20) Also, are the models after k iterations and are the a priori models.
The Jacobian matrix associated with the cross-gradient function is given by (21) where a component for the matrix is , ; .
In compact notation, equation 14 can be written as (22) where (23) In expressions 22 and 23, solves the least-squares minimization problem for separate estimation of both model parameters and agrees with equation 25 in Tarantola and Valette (1982), for separate inversion. The remaining terms are equality constraints given by the cross gradients, in addition to the original formulation of Tarantola and Valette (1982). Despite being different algebraically, equations 14 and 15 are equivalent to those obtained by Gallardo and Meju (2004) but are developed under a more robust statistical framework.
Analyzing the rank of matrices and involved in equations 14 and 15, we realize that although could be singular by itself, nonsingularity of the matrix can be assured by the covariance matrix of a priori models . Thus, can be calculated with any efficient inversion algorithm. Alternately, the matrix includes three equations for each couple of collocated parameters. Therefore, it represents a problem with more linearized equality constraints than parameters, which makes singular. Calculating the pseudoinverse of the matrix using SVD (Lines and Treitel, 1984; Strang, 1980) is a technique to solve this type of problem. Although it involves considerable computational effort, this is a good technique for inverting singular matrices and analyzing the condition of linear interdependency implicit in cross-gradient constraints.
The pseudoinverse of matrix is expressed by (24) for which U and V are orthonormal matrices with size , and Λ is the corresponding singular value matrix. These matrices can be partitioned as (25) for a given p value such that are considered as nonnull singular values.
The subsurface is discretized into rectangular prisms of homogeneous properties which constitute model parameters and . Based on these prisms, we used a forward-difference scheme analogous to that of Gallardo and Meju (2004) to discretize the cross-gradient vector and define (26)where relevant parameters are shown in Figure 2.
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 (27) for the case of , the first 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 are included according to the specific application.
Following previous experiences in cross-gradient-constrained problems in two dimensions (Gallardo and Meju, 2004; Gallardo et al., 2005), we select control parameters and of the same order of magnitude to achieve similar levels of smoothness in both models. Then we increase or decrease these regularization values in logarithmic steps until satisfactory data misfit and model convergence are achieved.
Covariance matrices for a priori models and can be selected according to geophysical parameter information provided by other sources. In our particular test experiments, we impose variance values about four orders of magnitude larger than parameter values to ensure that our choice of a priori model will not bias our estimated model parameters.
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 of a homogeneous rectangular prism at point is given by (28) where γ is the gravitational constant, r is the distance between the observation point and source point , and ρ represents the mass density. Limits , and are defined by the prism's edges, as shown in Gallardo-Delgado et al. (2003).
Similarly, the magnetic field of a homogeneously magnetized 3D prism is given by (29) where J is the total magnetization in the geomagnetic field direction , and 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 and T. Therefore, there is no need for a linearized approach for these functions. However, we still require an iterative optimal model search procedure because of the nonlinear cross-gradient constraint.
Synthetic test model
To prove our 3D cross-gradient joint-inversion algorithm, we used (2048) prisms in a volume that is in both horizontal directions and depth. The prisms are distributed on a grid of in x- and y-directions and in depth. Using this model, we set a cubic heterogeneity that is embedded in a homogeneous background of density and induced magnetization equal to 0. The cube is per side and it is buried at a depth of . Density contrast of the cube is and its induced magnetization is .
We calculated a total of 1681 -spaced data on the surface for both gravity and magnetic fields. Random noise was added to the data with standard deviation of and (2% of the maximum amplitude of their respective anomalies). The assumed direction of the inducing magnetic field is and . Figure 3a and d shows the resulting gravity and magnetic anomalies.
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:
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 and . Similarly, we set null a priori density and magnetization models with constant for the whole model space.
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 . We do this by setting , where For this case, and .
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 ( and A/m, correspondingly) to reduce the influence of the a priori model in the final estimates. Then, we applied several smoothing factors until a satisfactory data misfit and process convergence were achieved, in this case, . Separate inversion models were extracted during the first step of the process. Assisted by the corresponding singular value plots, we selected a p-value that coincides approximately with the number of parameters . We might expect that reducing the cutoff level of singular values should lead to less influence of the cross-gradient constraints, and therefore to models resembling those obtained after separate inversion. The reduction of the cutoff level of singular values will depend on incompatibility of the models and should accommodate a common structure when they involve large structural complexity.
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 and . Considering that the cross-gradient constraint acts exclusively on model structure and it does not constrain the geophysical values by itself, it seems reasonable that potential field data might require additional constraints to those used here, even when they are jointly inverted. Thus, we can expect that either reducing model complexity (by using homogeneity assumptions such as those implicit in Poisson's relation) or constraining selectively the ill-conditioned zones of the models might produce appealing density and magnetization values. However, we suggest keeping a close record of their biasing effect in every experiment.
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 (30) and (31) where subscripts g and m refer to gravity and magnetics, respectively, and and are the number of gravity and magnetic data, respectively. Convergence is obtained using (32) and (33) The fact that the rms values achieve a steady behavior indicates the process is selecting from equivalent models only those that are structurally similar without affecting the data misfit.
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 (, ). Gallardo-Delgado et al. (2003) and Pérez-Flores et al. (2004) describe the geologic settings in detail. The area consists of a sedimentary basin with an andesitic-granitic igneous basement, overlain by alluvial deposits. The basement is extensively exposed in the northern and southern bounds of the studied area and deepens up to below sea level (Gallardo-Delgado et al., 2003). Land, marine, and airborne surveys provided the 978 gravity and 812 magnetic data included in our field-inversion experiment (Figure 7). The assumed direction of the inducing magnetic field is and . The standard deviation considered was 0.1 and for gravity data and 50, 60, and for magnetic data.
For the inversion process, we selected the region in Figure 7, rotated to optimize the block distribution and discretized with a grid of blocks, which are wide, long, and thick. We excluded an area in the southern portion, located in a narrow mountain range, which could not be mapped accurately because of the cell size selected for the model. Following Gallardo-Delgado et al. (2003), we set a priori models, and , characterizing the basement with zero density and magnetization and the sedimentary filling with a density contrast of and a magnetization contrast of . Similar to the homogeneous smoothing experiment with the cubic model, covariance matrices of a priori models were selected large enough so as not to bias the models to predefined values with . Initial models were null homogeneous distributions of density and magnetization. We selected , taking into account the smoothness level of the models as well as the data misfit.
For joint inversion of these field data, we selected . Figure 8 shows rms and convergence values achieved throughout iterations in the experiment. These plots show that the data misfit achieved is about two times greater than expected data errors, which are close to unity. This is reasonable according to the grid size in our models. Even though the joint inversion fits the data as well as the separate inversion (Figure 8a), the convergence level greater than 50% (Figure 8b) shows differences in some parameter values. This means most of the models found through iterations adequately justify the data, and the cross-gradient constraint is indeed searching among equivalent models in terms of global data fit. The optimal model selected had a minimum convergence value for both processes simultaneously. It was the 18th iteration for the respective experiments.
Figure 9 shows results obtained for the field experiment in the representative sections and in Figure 7. Figure 9a–d corresponds to separate inversion and Figure 9e–h corresponds to the joint-inversion experiment. In addition, we show comparative profiles of calculated and observed data for each section, to give an account of model response and misfit level. Models found after joint inversion increase their geometrical resemblance as shown by the rms value of the cross-gradient function, which is one order of magnitude smaller than obtained by separately estimated models. As expected, the gravity field data characterize the sedimentary basin (with density contrast values of less than ), whereas the magnetic field data identify a magnetic basement (magnetization greater than ). Despite the fact that gravity and magnetic field data characterize different sources, the cross-gradient constraint yields a joint determination of the common boundary between the sedimentary basin and magnetic basement. Simultaneously, it permits the occurrence of a nonmagnetic basement, which is not distinguished individually by the gravity data.
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.
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.
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: (A-1) Here, derivatives are taken corresponding to random variables x, given by equation 9. Matrices involved in A-1 are defined as whereas , , , , , and are defined as , ; . For these matrices, and are the number of observed data related to the model and , respectively, although 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 (A-2) We observe that G contains sensitivity matrices for both models and matrices of smoothness operator derivatives. is a matrix that contains derivatives of cross-gradient vector-function t with respect to both model parameters.
To find the corresponding expression for equation 8, we divide this equation into the following three elements: , , and .
For our application, the first element (A-3) In compact notation, (A-4)
To find an expression for , we substitute from equation 12 and obtain (A-5)
Solving for our particular formulation, (A-6) Therefore, (A-7)
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 (A-8) We rewrite this expression as (A-9) where (A-10)
The estimated parameters of interest from A-9 are (A-11)
- Received April 24, 2008.
- Revision received December 12, 2008.