- © 2006 Society of Exploration Geophysicists. All rights reserved.
Gravity and magnetic data are inverted jointly in terms of a model consisting of an interface separating two layers having a constant density and magnetization contrast. A damped least-squares inversion is used to determine the topography of the interface. The inversion requires knowledge of the physical property contrasts across the interface and its average depth. Since the relationship between model parameters and data is weakly nonlinear, a constant damped least-squares inverse is used during the iterative solution search. The effect of this inverse is closely related to a downward continuation of the field to the average interface depth. The method is used to map the base of the Sept-Iles mafic intrusion, Quebec, Canada, and the shape of the central uplift at the Chicxulub impact crater, Yucatan, Mexico. At Sept-Iles, the intrusion reaches a thickness of , coincident with the maximum gravity anomaly, south of the intrusion center. At Chicxulub, the top of the central uplift is modeled to be deep and has a single peak form.
Multiple data sets offer the possibility of reducing ambiguity and nonuniqueness in geophysical data interpretation. The joint inversion of two data sets requires that each set of observations is caused by the same source, i.e., the spatial variation of each physical property is the same or, as commonly modeled, sources with constant physical properties have the same boundaries. If a magnetic field and gravity field are caused by the same source in which the magnetization direction and ratio of magnetization intensity to density are both constant, they are related by Poisson's relation. Hence, for noise-free data, one field can be calculated from the other. Thus, the utility of jointly inverting gravity and magnetic data would appear to be limited. However, where coincident gravity and magnetic anomalies appear to be caused by the same source, joint inversion provides a means to quantitatively assess whether this is indeed the case. If inversion results in good fits to both data sets, then a common source is indicated. However, poor fits to portions of either data set will pinpoint those parts of the source where physical property variations do not vary in concert and/or the boundaries of the sources causing the gravity and magnetic fields do not coincide. In these areas, the geological model must be adjusted or alternative interpretation methods could be used to satisfactorily model each data set independently. Bott and Ingles (1972) use an equivalent-layer approach to investigate the variation in magnetization/density ratio of sources without direct inversion of the data. Their method is independent of source type and is based on correlating coincident pseudogravity (calculated from the magnetic field) and gravity profiles. This approach does highlight variations in magnetization/density ratios and, consequently, implied variations in source boundaries; however, joint inversion provides additional information in the form of an acceptable geological model for the common source.
For many geological bodies, it is unlikely that magnetization (induced or remanent) will vary in exactly the same fashion as density. Density is primarily controlled by lithology, i.e., silicate mineralogy. In contrast, magnetization, which is predominantly determined by the occurrence of magnetite, is less directly related to lithology but is strongly influenced by metamorphism, alteration, and emplacement conditions. Variations in density are also likely to be much smaller in magnitude than those in magnetization, so the ratio between the two is expected to vary, e.g., between different phases of a single pluton. Consequently, a joint inversion of the gravity and magnetic data can delineate areas where the common source assumption is appropriate and where it is not. This provides useful geological information in terms of physical-property variation and source structure. The difference in resolution power between magnetic and gravity data also means that each data set responds differently to sources at various depths and/or levels of structural complexity. Magnetic data contain more high-frequency content than gravity data, so they are often more effective at resolving shallow or complex structures.
Several approaches to inverting gravity and magnetic data jointly are described in the literature. Menichetti and Guillen (1983) use a generalized inversion technique to determine body shapes for the 2.5D case with specified densities and magnetizations. A stable inversion is achieved by removing the effects of small eigenvalues, providing a suggested improvement over separate inversion of either data set. Serpa and Cook (1984) use a least-squares inversion method to solve for the vertices and physical properties of a 2.5D model. Their method does not include damping or regularization, so stabilizing the solution is effected by holding some parameters fixed while inverting for others. Zeyen and Pous (1993) use a 3D subsurface model comprising of an ensemble of vertical prisms and solved for physical properties and the tops and bottoms of each prism. They control solution stability by specifying a priori parameter and data covariances. Gallardo-Delgado et al. (2003) extend the 3D approach to include a continuous density variation with depth and the direction of magnetization as unknown parameters. Within a statistical framework, Monte Carlo methods have been used for joint gravity and magnetic inversion producing an a posteriori probability density function describing a suite of acceptable models rather than a single inverted model (Bosch, 1999; Bosch et al., 2001).
All of these inversion methods are based on solving a combined set of equations relating the model to both data sets. Depending on the model's definition, the problem can be posed with the number of data points either greater or smaller than the number of model parameters. In either case, the problem can be solved by minimizing (1) the difference between both observed and predicted data sets plus, for solution stability, (2) a measure of model length or roughness. This leads to methods that find a minimum norm solution (Menke, 1989). Tighter constraints on linking the different physical properties being solved for can be incorporated into the inversion by introducing some measure of structure into the model. Haber and Oldenburg (1997) define a measure of the model structure based on its curvature and then minimize the difference in structure between the two models subject to the constraint that both data sets are fit to within their expected errors. A similar approach is used by Gallardo and Meju (2003), who use a structure measure based on the vertical and horizontal gradients of the physical properties being solved for. Both of these approaches are again posed as minimization problems, with extra terms for structure added to the standard terms of data misfit and model length.
In this study, gravity and magnetic data are inverted jointly using a damped least-squares approach. The standard way of dealing with nonlinearity is by reevaluating matrix quantities at each model estimate during an iterative search for an acceptable solution. Since matrix evaluation and/or inversion is computationally expensive for large data sets and models, a fixed inverse with a simple form is used throughout the search. This approach is demonstrated on synthetic data and applied to coincident gravity and magnetic anomalies over the Sept-Iles mafic intrusion in Quebec, eastern Canada, and the Chicxulub impact crater, Mexico. For the two real data examples, joint inversion is used to determine how well the anomalies are described in terms of a common source and to focus on areas where this assumption breaks down.
The model used in the following consists of two homogeneous layers, each having a constant density and magnetization, separated by an interface whose topography is to be determined. The nonuniqueness associated with this kind of model is predominantly related to the depth and physical-property contrast of the interface (see Skeels, 1947). For a constant interface depth, varying the physical property contrasts between the two layers results in a range of models with an inverse relationship between the interface topography and the density or magnetization, all fitting the data with the same precision. For the gravity case, in general, the average interface depth need not be known a priori; but for the algorithm used below, it must be specified, since this parameter is not solved for. Similarly, for the magnetic case, the average interface depth must be known since the magnetic data do not provide information on its value. Consequently, the average interface depth and the density and magnetization contrast between the two layers are specified prior to inversion.
The expressions relating to gravity data g and magnetic data m to the model x are given by Parker (1973): (1) and (2) where j and t are unit vectors in the direction of the magnetization and the earth's magnetic field, ρ is the density contrast, J is magnetization contrast, G is the gravitational constant, F[ ] denotes Fourier transform, and k is the wavevector. Vector K is constructed from wavevector k, whose components in the coordinate system are , by . The average interface depth z is measured positive downward, and x is the interface topography measured (positive upward) with respect to z. Since equations 1 and 2 consist of a series of Fourier transforms, g, m, and x must be located on coincident grids and the appropriate treatment applied to minimize edge effects. Linearized versions of the above relations are, in partitioned matrix form, (3) or simply (4) with G and M can be written as (Pilkington and Crossley, 1986) (5) and (6) where E and are the matrix equivalents of Fourier and inverse Fourier transformation and the asterisk denotes complex conjugation. Matrices Γ and Λ are diagonal with elements (7) and (8) The factorization of G and M in equations 5 and 6 results from taking the first term in the series in equations 1 and 2 and setting x to zero. Since E and are unitary and Γ and Λ are diagonal, the factorizations in equations 5 and 6 are equivalent to the singular value decompositions (SVD) of G and M. Hence, Γ and Λ contain the eigenvalues of G and M, respectively.
Because the data sets have different units and magnitudes, one or both of them must be weighted to balance their contribution in the inversion. If, for example, numerical values of M and m are much greater than G and g in equation 3, the inversion will be dominated by the magnetic data and information from the gravity data will be lost. If both data sets are assumed to be contaminated by independent Gaussian errors, they can be transformed into sets of dimensionless values by dividing each observation by its expected standard deviation (Vozoff and Jupp, 1975; Serpa and Cook, 1984). If errors are dependent, the data transformation requires knowledge of the data covariances but is still implemented straightforwardly (Jackson, 1973). A similar weighting method consists of normalizing parameter and observed data sets by a priori model values (Hering et al., 1995). Finally, empirical weighting guided by estimates of data variances can be used. This is a flexible approach in that the weighting can be adjusted even if it requires that weights are not in agreement with estimated data variances (Lees and VanDecar, 1991; Zeyen and Pous, 1993).
In the statistical approach to inversion, weighting is quantified by the ratio of the standard deviation of the magnetic data errors to that of the gravity data. Determining reliable and usable estimates of data errors is difficult, to say the least, since the data used in inversion contain not only raw measurement errors but also the combined effects of error introduced by various data-reduction steps and, most serious of all, gridding. Consequently, estimates of data errors can be used to determine the weighting; but if it is clear that one data set is dominating the inversion solution, then the weighting should be adjusted accordingly. The empirical approach adopted here is implemented by weighting just one set of equations, so that equation 3 becomes (9) where w is the scalar weighting factor. An alternative approach that avoids the weighting issue is sequential inversion (Lines et al., 1988), in which one data set is inverted and the results are used as the starting model for inversion of the second data set. Coupling the two data sets may, however, prove weaker for this method.
Equation 9 is an overconstrained system solved by forming the damped least-squares solution (Menke, 1989) for the equivalent equation 4: (10) where I is the identity matrix, and θ is the damping factor. Substituting for G and M using equations 5 and 6 (11) where (12) and (13) Note that the contributions of each data set to the model are completely separable but are linked through the bracketed terms. Equations 11–13 can also be derived based on the SVD of A and forming its generalized inverse. Because , Λ, and are all diagonal matrices, determining x can be carried out efficiently in the wavenumber domain using gridded data without any need for explicit matrix inversion. Nonlinearity is taken care of by using the iterative scheme (14) where k is iteration, and are the calculated responses from model and , are the observed data. Equation 14 differs from the standard Gauss-Newton schemes used in nonlinear inversions by having a constant inverse that is not reevaluated at each model . The inverse used here is equivalent to that of a standard Gauss-Newton scheme with the inverse just evaluated once at the starting model of zero, i.e., a flat interface. The simple form of equations 5 and 6 is lost for nonzero and would require the inversion of in equation 9 either explicitly or through the solution of a large set of equations.
Use of a constant inverse in equation 14 is possible and effective because of the weak nonlinearity of the expressions in equations 1 and 2. Updating and at each iteration may result in better convergence than using a constant and , but the overall efficiency is reduced because of the computational burden of determining the derivatives and solving equation 9 at each iteration. Based on a similar problem using only profile data, Pilkington and Crossley (1986) compare the approximate and updated approaches and show the former is considerably faster (two orders of magnitude) computationally. A necessary requirement for convergence of equation 14 is that the inverse used will map small changes in the data residuals into small changes in the current model. The damping parameter θ, which is kept constant throughout the iterative procedure, also contributes to the stability of the iterations by suppressing the effects of small eigenvalues to the current model estimate. Any small eigenvalues in G and M lead to large numerical values of the terms in parentheses in equations 12 and 13, containing the inverse eigenvalues. If the latter become large, this can cause unacceptably high-amplitude, oscillatory behavior in the model and equation 14 diverges. Damping is always necessary because of the exponential terms in the eigenvalues (equations 7 and 8). The inversion is equivalent to a downward continuation of the field and conversion of the field at the average interface depth into topography through division by the nonexponential terms in equations 7 and 8. In the gravity case, this is equal to converting the continued field into topography using the Bouguer slab formula.
The effect of introducing damping through addition of θ is clarified when the damped inverse eigenvalues are compared with the undamped versions. Since the index i in equations 7 and 8 corresponds to wave vector, forming the ratio of damped to undamped inverse eigenvalues quantifies the frequency response of the smoothing introduced by the damping factor. Figure 1 shows an example case where θ is numerically equal to one of the eigenvalues ( in this case). Choosing θ from the range of or ensures that it is in the right numerical range to impose a useful amount of damping. If the solution is overdamped, the iterative search becomes more like the steepest descent method, and convergence slows. If there is no effective damping and oscillations in the model result from overamplified high-frequency components in the data causing divergence. The wavenumber response of the damping in Figure 1 is equivalent to a low-pass filter or smoothing operator. Both gravity and magnetic eigenvalues (equations 7 and 8) decay rapidly as wavenumber increases, so taking eigenvalues corresponding to the large wavenumbers results in less smoothing of the solution.
SYNTHETIC DATA EXAMPLE
The application of the iterative inversion scheme (equation 14) is demonstrated on a synthetic data set. Figure 2a shows a model topography with an average depth of , density contrast of , and induced magnetization of with inclination 70° and declination 30°. Figure 3a and b shows the gravity and magnetic fields produced by the model. These data sets were inverted jointly, resulting in the solution in Figure 2b after seven iterations. The contribution of each data set to the final solution is given in Figure 2c and d, which show, respectively, the solution components from the gravity and magnetic data. As expected, the magnetic contribution reflects greater power in the field at short wavelengths (see Figure 1), while the gravity contribution has a much smoother character. Figure 4 shows the maximum topographic variation for the gravity and magnetic contributions as a function of weighting factor w. The inversion in Figure 2 uses , corresponding to almost equal contributions to the solution, in terms of amplitudes, from each data type. For this synthetic example the gravity and magnetic sources are equivalent, so varying the weighting has little effect on the inversion results.
SEPT-ILES MAFIC INTRUSION, QUEBEC
The Sept-Iles mafic intrusion is located on the north shore of the St. Lawrence River in Quebec, Canada. Cambrian in age, it has a diameter of and is circular in shape with some elongation to the south-southeast. Only 5% of its extent is exposed, the remainder being covered by glacial deposits and water (Loncarevic et al., 1990). Geophysical data have therefore played a major role in delineating and characterizing the nature of the body. A combination of land, sea-bottom, and shipborne gravity measurements define a circular, positive, approximately anomaly (Figure 5a). The maximum is located south-southeast of the center of the main anomaly and is surrounded by a crescent of reduced values. The edge of the anomaly has its maximum gradients in the north, over land, while the southern portion shows a more gradual transition to higher-than-average background levels. Coincident with the gravity anomaly is a positive, magnetic anomaly (Figure 5b) with internal variations in character that match the changes in the gravity field. The magnetic data, however, reach maximum values of along the northern edge of the intrusion, while a secondary crescent-shaped high with values of occurs to the south and coincides with a similarly shaped gravity high. Again, gradients in the magnetic data are greatest along the northern edge of the body but are much reduced to the south.
So far, interpretation of geophysical data over Sept-Iles has consisted of 2D modeling of gravity data only. Goodacre and Nyland (1966) model the gravity anomaly as a vertical cylinder with a radius of , depth of , and density contrast of . Thomas (1974) interprets it in terms of a basin-like gabbroic body with a single density contrast of and a maximum depth extent of . Although Thomas (1974) uses only a single density for the body, he notes that the concentricity of the anomaly pattern suggests layering within the source. Based on mapping along the north Saint Lawrence shore, Loncarevic et al. (1990) divide the intrusion into a layered series ranging from ultramafics at the deepest levels through mixed gabbro, troctolite, and anorthosite at intermediate levels up to monzogabbro and mafic syenite near the surface. They assign density contrasts of 0.54, 0.25, and to the series, respectively, and carry out 2D modeling along two radial profiles. They interpret a funnel-shaped structure for the intrusion with a maximum thickness of .
The likelihood that both the gravity and magnetic fields over the Sept-Iles intrusion are predominantly the result of the same structure is demonstrated by the good coincidence of the two in Figure 5. Joint inversion would appear appropriate in this case. However, the observed layering of the intrusion and consequent variation of physical properties throughout the body means that the assumption of a single density/magnetization contrast for the inversion is probably too simplistic. Nevertheless, joint modeling of the observations using the two-layer model with just a single and constant physical property contrast can still provide useful information in terms of pinpointing areas where sources have more complex property variations.
Paleomagnetic studies of the Sept-Iles intrusive rocks indicate a stable remanent magnetization component (inclination , declination 333°) of Cambrian age (Tanczyk et al., 1987). Intensities of are also reported that, given the large amplitude of the anomalies over Sept-Iles, suggest induced magnetization predominates over the remanent component. Inspection of anomaly shapes and comparison to model effects calculated with the measured Cambrian magnetization clearly supports this position. The choice of susceptibility and density values to use in the inversion was guided by trial inversions of the magnetic and gravity data separately. The values (0.20 SI, ) that produced topographies similar in magnitude to the 2D values from Loncarevic et al. (1990) were adopted. These were later increased to 0.24 SI and after trial joint inversions and inspection of the residual fields. Several inversions were run, decreasing the damping factor to the smallest possible value (for greatest resolution) that did not cause divergence or any oscillatory behavior in the solution. An example inversion result is given in Figure 6a, which shows the topography of the base of the intrusion. The gravity and magnetic effects of the solution are given in Figures 5c and 5d. There is generally a good match between observed and calculated fields; the primary features in both fields over the intrusion can be accounted for by the computed model. The calculated fields are slightly smoother than the observed because of damping the solution.
Compared to the Loncarevic et al. (1990) 2D model, which shows a simple increase in depth toward the center of the intrusion, the model in Figure 6a has a more complex form, with its deepest point coinciding with the maximum gravity anomaly south-southeast of the intrusion center. Secondary features are the two crescent-shaped troughs to the north of this point, the inner coinciding with a magnetic and gravity high and the outer associated with only a magnetic high, on the northern edge of the body (see Figure 5). Both troughs are the result of the more dense and more magnetic late gabbro phase that forms ring dikes of varying thickness (Loncarevic et al., 1990). They show positive magnetic residuals but are well fit in the gravity case (Figure 6b and c), suggesting a higher magnetization value is appropriate for these parts of the intrusion. Over the gravity maximum ( north, east) magnetic residuals show a balance between positive and negative, whereas gravity residuals are generally positive, showing that model structure in this region is mainly determined by the magnetic data. This also appears the case for the shorter wavelength variations in the solution topography in the region just north of the gravity maximum. The gravity data here are unable to constrain the shorter wavelength portions of the topography as well as the magnetic data. To the northeast of the intrusion, the large gravity high ( north, east), although not completely fitted, dominates the solution, resulting in large negative magnetic residuals. Similarly, the extension of the gravity high over the intrusion to the southwest ( north, east) is not coincident with a magnetic high. As a result, the inversion solution shows deeper topography values in this region based on the gravity data but produces conspicuous negative magnetic residuals.
In the case of the Sept-Iles intrusion, assuming a constant source density and susceptibility leads to artifacts in the solution that would not appear in a model where physical properties are variable. Nevertheless, the joint inversion demonstrates that most of the observed gravity and magnetic data result from the same source. Inspection of the residuals for both data types delineates areas where property variations and/or source boundaries do not vary in concert and indicates how they deviate. This kind of information can form the basis for subsequent detailed modeling.
CHICXULUB IMPACT STRUCTURE, YUCATAN
The second data example is taken from the central portion of the Chicxulub impact structure located on the northwestern Yucatan coastline. This -diameter crater is buried approximately below the surface and was originally identified by its distinctive gravity and magnetic signatures (Penfield and Camargo, 1981). A characteristic of craters this size is a central uplift, which is caused by rebound of the crater floor immediately after excavation of the initial cavity. After impact, material comprising this central uplift can be many kilometers above its preimpact position. Therefore, denser and more magnetic rocks often make up the central uplift. This is certainly the case at Chicxulub, where the crater was formed in a 3 to -thick carbonate succession overlying crystalline basement. Previous interpretations of the gravity and magnetic data (Sharpton et al., 1993; Hildebrand et al., 1998; Pilkington and Hildebrand, 2000) indicate that gravity and magnetic highs at the center of the structure are caused by some form of central uplift, with a width of about and top at an estimated depth. However, in these studies, each data type was modeled in isolation. Espindola et al. (1995) model both data types in terms of a central uplift but only along a single profile. Here, residual gravity and magnetic data are inverted jointly over the central portion of the structure. The gravity effect of regional basement structures unrelated to the crater were modeled and removed (Hildebrand et al., 2003); then magnetic data were processed with a low-pass filter to remove the effects of short-wavelength magnetization changes within the crater's melt sheet that occur at depths of (Pilkington and Hildebrand, 2000).
Figure 7a and b shows the residual gravity and magnetic fields used in inversion with the corresponding calculated fields from the topography solution in Figure 8a, given in Figure 7c and d. Initial density and susceptibility values were set to values previously used in 2D modeling. These were adjusted to and 0.14 SI after some trial inversions.
Figure 9 shows the behavior of the gravity and magnetic data fits as a function of varying density and susceptibility values and carrying out inversions. The particular solution depicted in Figure 8a represents a compromise between the two data fits, with rms values of and . Comparing the observed and calculated data in Figure 7 shows that this solution does not represent an exemplary fit to the data. However, adjustment of the damping factor, average interface depth, and weighting factor did not improve the range of solutions shown in Figure 9. So these results are the best for the assumptions imposed. Clearly, the observed double gravity high (Figure 7a) will not produce a topography compatible with the single maximum seen in the magnetic data (Figure 7b). The solution in Figure 8a is a combination of both fields, with a widening of the uplift to the southwest, reflecting the two gravity maxima, while the single central magnetic high is responsible for the topographic maximum. The largest residuals for both the gravity and magnetic data occur over the southwest portion of the central uplift (Figure 8b and c), indicating that a common source cannot completely explain both data sets in this region. The southern lobe of the observed gravity maximum reaches , and the corresponding residual is approximately , so the assumption of a common source is still able to account for most of the power in the observations. The central and northeastern portions of the structure, with much smaller residuals for both data types, appear well described in terms of a common source consisting of a simple interface with constant physical-property contrasts. The residual maps are anticorrelated, suggesting that, for the weighting factor used, the inversion gives equal weight to both sets of data. If the inversion were recomputed with a higher weighting for the gravity data, the solution would take on more of a double-peaked form, reducing the gravity residuals but increasing the mismatch with the magnetic data.
Relaxing some of the prior assumptions, e.g., the constant density and magnetization, would likely produce a better fit to the data. Since there is no strong observational or theoretical evidence to support the existence of more than one peak at central uplifts in impact craters, the gravity data could be better modeled with a variable density structure. Preexisting basement density contrasts will persist through uplift formation, so a single constant density as assumed here may not be the optimal choice. Conversely, the magnetic data appear compatible with a constant magnetization throughout the uplift.
For simple two-layer models, gravity and magnetic data can be inverted jointly using a frequency-domain iterative scheme requiring no matrix inversion. The frequency response of the inverse eigenvalues used in the inversion algorithm (equation 14) show that the magnetic data contribute more short-wavelength information to the solution than the gravity data. The synthetic example also demonstrates this point.
The inversion of data sets from an intrusion and an impact crater central uplift shows that a simple common source assumption for these two bodies could explain most of the observed anomalies. The presence of residuals for either data type indicates areas that could not be satisfactorily explained in terms of the assumptions adopted. These regions are most likely caused by spatial variations in either or both density and magnetization. If such a subarea can be isolated, in which physical properties are constant (but not equal to the initial inversion values), then these property values can be adjusted and the data reinverted. Where the physical-property variations exist over small scales, alternative modeling methods must be used. The two data examples show that even the restrictive assumption of constant physical properties for the models can result in meaningful structural interpretations and provide a solid starting point for more detailed subsequent modeling if needed.
John Millar (University of Calgary) is thanked for providing the residual gravity data over Chicxulub. Two reviewers and the associate editor (J. Silva) provided detailed comments and suggestions that improved an earlier version of the manuscript. This is Geological Survey of Canada contribution 2004016.
- Received May 14, 2004.
- Revision received June 15, 2005.