Download Estimating Complex Gravity Sources with Varying Densities: Interactive Inversion - Prof. X and more Papers Engineering in PDF only on Docsity!
GEOPHYSICS, VOL. 71, NO. 1 (JANUARY-FEBRUARY 2006); P. J1–J9, 15 FIGS.
10.1190/1.
Interactive gravity inversion
Jo ˜ao B. C. Silva 1 and Val ´eria C. F. Barbosa 2
ABSTRACT
We have developed a new approach for estimating the location and geometry of several density anoma- lies that give rise to a complex, interfering gravity field. The user interactively defines the assumed outline of the true gravity sources in terms of points and line seg- ments, and the method estimates sources closest to the specified outline to achieve a match between the pre- dicted and observed gravity fields. Each gravity source is assumed to be a homogeneous body with a known density contrast; different density contrasts may be as- signed to each source. Tests with synthetic data show that the method can be of use in estimating (1) multiple laterally adjacent and closely situated gravity sources, (2) single gravity sources consisting of several homo- geneous compartments with different density contrasts, and (3) two gravity sources with different density con- trasts of the same sign, one totally enclosed by the other. The method is also applied to three different sets of field data where the gravity sources belong to the same cat- egories established in the tests with synthetic data. The method produces solutions consistent with the known geologic attributes of the gravity sources, illustrating its potential practicality.
INTRODUCTION
Potential field-data interpretations aimed at locating the horizontal projection of source boundaries and trends are usually performed via linear transformation techniques such as derivatives, shaded relief maps, upward and downward continuations, and apparent density and susceptibility map- ping representations (e.g., Kowalik and Glenn, 1987; Arkani- Hamed and Urquhart, 1990; Broome, 1990; Grauch et al., 2001). These techniques are commonly implemented via user- friendly software for image processing and interactive graphi- cal interpretation. Complex multiple laterally adjacent and
Manuscript received by the Editor December 18, 2003; revised manuscript received May 10, 2005; published online January 12, 2006. (^1) Federal University of Par ´a, Dept. of Geophysics, Caixa Postal 1611, Bel ´em, Par ´a 66017-900, Brazil. E-mail: joaobcs@directbr.com.br. (^2) LNCC, Av. Get ´ulio Vargas, 333, Quitandinha, Petr ´opolis, Rio de Janeiro, 25651-075, Brazil. E-mail: valcris@lncc.br. © c 2006 Society of Exploration Geophysicists. All rights reserved.
closely situated gravity sources can be handled easily by these methods. On the other hand, interpretations involving location and delineation of the sources themselves instead of their horizon- tal projections are usually performed in two different ways. The first approach comprises all inversion methods, automat- ically determining the position and geometry of a causative body, provided that sufficient a priori information about the source is incorporated by the method (e.g., Last and Kubik, 1983; Guillen and Menichetti, 1984; Barbosa and Silva, 1994; Li and Oldenburg, 1998; Silva et al., 2000). The advantage of this approach is its efficiency in finding a solution that not only fits the observations within the measurement errors but also possesses the desirable physical attributes specified by the interpreter and incorporated by the method. However, the application of this approach to complex geological settings is severely limited; the anomalous sources must be isolated or must display a relatively simple geometry — or both. In ad- dition, all geologic information must be mathematically trans- lated and automatically incorporated by the method, leaving little or no room for interactive supervision by the user. The second approach, frequently adopted when true gravity sources are close to each other (either vertically and laterally) and possess complex shapes, is the interactive 2D modeling method, which imposes virtually no limitations on the com- plexity of the interpreted source (e.g., Paul and Bain, 1998; Abbott and Louie, 2000; Grauch et al., 2001). A drawback of this approach is the tremendous difficulty in obtaining, in some cases, a reasonable fit of the observations. Even when the interactive 2D modeling produces correct fits, the solution may not match the geology unless enough other constraints are available to limit the choices. We present a new approach for interpreting 2D gravity anomalies produced by multiple and complex gravity sources that are separated (vertically and/or laterally) from each other by short distances. This is a step forward in combining the best features of automatic inversion and interactive model- ing. The assumed interpretation model is a grid of juxta- posed 2D prisms whose density contrasts are the parameters to be determined. The interpreter specifies, in a user-friendly
J
J2 Silva and Barbosa
environment, the outlines of the gravity sources in terms of geometric elements (line segments and points) and the den- sity contrast associated with the geometric elements defining each gravity source framework (this amounts to specifying the assumed density contrast for each source). The method then estimates the density-contrast distribution that fits the observed anomaly within the measurement errors and repre- sents compact gravity sources closest to the specified geomet- ric elements. The user can either accept the interpretation or modify the gravity-source framework, changing the position of geometric elements and/or the density contrast associated with each of the elements and restart the inversion. We use Guillen and Menichetti’s (1984) inversion method, modified to permit the interpretation of anomalies caused by complex, closely situated gravity sources. The first mod- ification consists of allowing different density contrasts to be assigned to different tentative gravity sources (we use the term tentative sources to describe the sources the interpreter pre- sumes exist). With this facility, multiple, complex, and closely separated gravity sources with different density contrasts can be delineated. The second modification consists of combin- ing Guillen and Menichetti’s (1984) methodology with a ro- bust procedure to interpret anomalies caused by small gravity sources embedded in larger ones. Our method’s potential value in producing stable and ge- ologically meaningful results is illustrated by inverting syn- thetic data produced by complex simulated geological settings. Three real-data profiles are interpreted with our proposed method. The first profile consists of several positive anoma- lies produced by metabasalts and metagabbros from a green- stone belt located in the Rio Maria region in the state of Par ´a, Brazil. The second one is a negative gravity anomaly produced by the Bodmin Moor Granite, which is part of the Cornu- bian batholith in southwestern England. The third anomaly is produced by the layered East Bull Lake gabbro-anorthosite intrusion in northern Ontario, Canada. In all cases, the struc- tures obtained by the interactive gravity inversion are consis- tent with the known geological attributes of the true gravity sources. When compared with interactive modeling, our re- sults are easier and quicker to obtain, and our method may produce a better anomaly fit.
Figure 1. Interpretation model consisting of a set of 2D ver- tical juxtaposed prisms whose density contrasts are the pa- rameters to be estimated. The outlines of anomalous grav- ity sources Sr , r = 1 ,... , R, are defined from a set of L presumably known geometric elements (axes and points) ei , i = 1 , 2... , L.
METHODOLOGY
Let S (^) r , r = 1 , 2 ,... , R, be a set of 2D gravity sources having arbitrary shapes and arbitrary density-contrast distributions, and assume that outlines of these sources may be constructed by a combination of axes and points totaling L geometric el- ements (Figure 1). Let T be the set of all geometric elements ei , i = 1 , 2 ,... , L (points and axes) ordered in an arbitrary way. Each element ei of T is assigned a target density con- trast. Additionally, we assign to the j th gravity source a subset t (^) j of T containing Kj geometric elements andKj target density contrasts. By combining (1) the presumably known sources outline (axes and points), (2) the corresponding target density contrasts, and (3) the measurements of the gravity anomaly produced by the R gravity sources, we can improve the source delineation. To this end, we first assume an interpretation model con- sisting of an Nx × Nz grid of 2D vertical juxtaposed prisms (Figure 1) whose density contrasts are the only unknown pa- rameters. The gravity anomaly gi ≡ g(xi ) produced by such an interpretation model at x = xi is given by
g i =
∑^ M
j = 1
p j Aij , i = 1 , 2 ,... , N, (1)
where Aij is numerically equal to the gravity anomaly pro- duced at xi by the j th prism with unit density contrast, N is the number of observations, M is the total number of prisms, and p (in equation 2) is an M ×1 vector of unknowns whose el- ement p (^) j is the density contrast of the j th prism. Using matrix notation, equation 1 becomes
g = Ap. (2)
Here, g is an N × 1 vector whose ith element is gi and A is an N × M matrix whose generic element is given by Aij. To obtain a stable solution of the linear system given in equation 2, we look for the solution satisfying the gravity anomaly and presenting most of its mass excess (or defi- ciency) concentrated about the specified geometric elements. This is accomplished by generalizing the iterative approach of Guillen and Menichetti (1984), which consists of the following steps. First, we obtain a standard minimum-norm solution,
p
o
= A T^ ( AA T^ + μ I )−^1 g , (3)
where μ is a nonnegative scalar, T is a transposition opera- tor, and I is the identity matrix. The larger the value of μ, the smaller the Euclidean norm of p
o
. Then, we update the param- eter estimate iteratively by
p
(k+1)
= pˆ
(k)
F +^ ∆ p
(k), (4)
where
p (k)^ = W −(k^1 ) A T^
AW −(k^1 ) A T^ + μ I
g − A ˆp ( Fk)
The symbol W (k) represents a diagonal matrix whose nonzero elements are given by
w jj =
d j^2
∣ pˆ( jk)
J4 Silva and Barbosa
to the anomaly produced by the larger source than to the ob- served anomaly. Then, the observation producing the maxi- mum residual is replaced by the fitted anomaly at this position. This new data set is then used as input for the next inversion, and the rms of the difference between the input and the fitted data is computed. The process is repeated until successive rms values along the iterations do not change appreciably. Upon convergence, the estimated density-contrast distribution will approximate the larger source; hence, the fitted anomaly is subtracted from the original observations to produce an estimate of the anomaly produced by the smaller source, which is then in- verted by the same technique described in the previous section.
TESTS WITH SYNTHETIC DATA
In this section we apply our proposed method to interpret- ing data from three simulated environments. We also analyze our method’s sensitivity to introducing incorrect a priori infor- mation.
Figure 2. Anomalous gravity source S 1 entirely surrounded by anomalous source S 2.
Figure 3. Synthetic example simulating multiple, complex, and laterally adjacent gravity sources. (a) Noise-corrupted Bouguer (crosses) and fitted (blue solid line) anomalies. (b) Simulated salt canopies (black, solid thick lines) having a den- sity contrast of −0.2 g/cm^3 and the inversion result (colored cells) using four presumably known axes (green lines).
Complex, multiple, and closely separated gravity sources
To illustrate the utility of our approach in interactive in- version of a gravity profile caused by multiple gravity sources with complex shapes and closely situated to each other, we model three salt canopies with roots having a density con- trast of −0.2 g/cm^3. The interpretation model consists of a grid of 160 × 64 cells with dimensions of 0.125 km in both the x - and z -directions and inversion parameters of μ = 0.2, f = 50 000, and τ = 0.01. The anomaly is contaminated with pseudorandom Gaussian noise with zero mean and a standard deviation of 0.1 mGal. All steps — design of the simulated gravity source, specification of the parameters related to the noise contamination of the anomaly, input of the geometric elements, and choice of the inversion parameters — are done interactively in a user-friendly environment. Figure 3 shows the output screen for this test. Figure 3a shows the observed (black crosses) and fitted (solid blue line) anomalies. Figure 3b shows the true gravity sources (black solid lines), the axes defining the sources (outlined in green), and the inverted density contrasts, mapped according to the color bar. The target density contrasts assigned to all axes are equal to −0.2 g/cm^3. We note that the observations are fit- ted within the measurement errors and the estimated gravity sources are close to the true ones because the specified geo- metric elements reflect factual geometric attributes of the true sources.
Single source with variable density contrasts
Figure 4a shows the gravity anomaly (dots) contaminated with pseudorandom Gaussian noise with zero mean and a standard deviation of 0.2 mGal produced by a simulated batholithic intrusion presenting two homogeneous sectors with density contrasts of 0.25 and 0.41 g/cm 3 for the leftmost and rightmost sectors, respectively (solid lines in Figure 4b). We inverted this anomaly using an interpretation model con- sisting of a grid of 80 × 64 cells with dimensions of 0.25 and 0.0625 km in the x- and z-directions, respectively, and setting μ = 0.25, f = 5000, and τ = 0.01. The sectors were assumed to be delineated by two axes each (labeled e 1 and e 2 for the leftmost sector and e 3 and e 4 for the rightmost sector in Fig- ure 4b). The target density contrasts assigned to the axes de- lineating the leftmost and rightmost sectors were 0.25 and 0.41 g/cm^3 , respectively. Figure 4c shows the spatial distribu- tion of the target density contrast assigned to each cell. The inversion results are shown in Figure 4b, indicating both sec- tors are well delineated.
Small source embedded in larger source
Figure 5a shows the gravity anomaly (dots) contaminated with pseudorandom Gaussian noise with zero mean and a standard deviation of 0.1 mGal produced by a simulated gabbro pendant in granite with density contrasts of 0.6 and 0.2 g/cm^3 , respectively (solid lines in Figure 5b). First, we in- verted this anomaly using the robust procedure described in the Methodology section by assuming an interpretation model consisting of a 176 × 64 grid with dimensions of 0.0625 km in both the x- and z-directions. We set μ = 0.2, f = 50 000, and τ = 0.01 and assumed that a single point outlines the larger
Interactive Gravity Inversion J
gravity source ( e 1 in Figure 5b) that was assigned a target den- sity contrast of 0.2 g/cm^3. The inversion result is shown in Fig- ure 5b, indicating the larger source has been well delineated. Next, we subtracted the fitted anomaly shown in Figure 5a (solid line) from the gravity anomaly shown in Figure 5a (dots) to obtain an estimate of the anomaly produced by the smaller gravity source with a density contrast of 0.4 g/cm^3 (dots in Fig- ure 6a). This residual anomaly was then inverted by the stan- dard, nonrobust procedure by setting μ = 0.75, f = 50 000, and τ = 0.01 and assuming that a single point outlines the smaller source ( e 2 in Figure 6b), which was assigned a target density contrast of 0.4 g/cm^3. The inversion result is shown in Figure 6b, and the fitted anomaly is shown in Figure 6a (solid line). This result indicates that the smaller gravity source has been delineated correctly. Figure 7a shows the fitted anomaly resulting from the sum of the anomalies displayed in Figure 5a and Figure 6a, pro- duced by the sum of the inverted density-contrast distributions shown in Figure 5b and Figure 6b.
Figure 4. Synthetic example simulating a single gravity source with variable density contrasts. (a) Noise-corrupted Bouguer (dots) and fitted (solid line) anomalies. (b) Simu- lated batholithic intrusion presenting two homogeneous sec- tors (solid lines) having density contrasts of 0.25 g/cm 3 (left- most sector) and 0.41 g/cm 3 (rightmost sector) and inversion result (grayscale cells), assuming axes e 1 and e 2 (dashed lines) for outlining the leftmost sector with a target density contrast of 0.25 g/cm 3 and axes e 3 and e 4 (dashed lines) for outlining the rightmost sector with a target density contrast of 0.41 g/cm^3. (c) Target density contrasts assigned to the cells of the inter- pretation model. Dark gray cells have a target density contrast of 0.41 g/cm 3 , light gray cells have a target density contrast of 0.25 g/cm^3 , and white cells have a target density contrast of 0 g/cm 3.
Sensitivity to incorrect a priori information
Figure 8 shows the gravity anomaly contaminated with pseudorandom Gaussian noise with zero mean and a standard deviation of 0.1 mGal produced by the simulated rootless salt canopy with density contrasts of −0.4 g/cm^3 , shown as a solid line in Figure 9. We inverted this anomaly using an interpre- tation model consisting of a 120 × 32 grid measuring 0.25 km in both the x- and z-directions. We set μ = 0.1, f = 500, and τ = 0.05 and assumed that three axes outline the gravity source ( e 1 , e 2 , and e 3 in Figure 9). The inversion results are shown in Figure 9.
Figure 5. Synthetic example simulating a small gravity source embedded in a larger source. (a) Noise-corrupted Bouguer (dots) and fitted (solid line) anomalies. (b) Simulated gabbro (inner body in solid line) enclosed in granite (outer body in solid line) having density contrasts of 0.6 and 0.2 g/cm 3 , re- spectively, and robust inversion result (grayscale cells) using point e 1 (dot) for outlining the outer body with a target den- sity contrast of 0.2 g/cm^3.
Figure 6. Synthetic example simulating a small gravity source embedded in a larger source. (a) Residual anomaly (dots, obtained by subtracting the gravity anomaly from the fitted anomaly in Figure 5a) and fitted anomaly (solid line). (b) In- version result (grayscale cells) using point e 2 (dot) for outlin- ing the inner body and a target density-contrast of 0.4 g/cm^3.
Interactive Gravity Inversion J
TESTS WITH REAL DATA
In this section we illustrate our method’s practical applica- tions by applying it to three sets of gravity data from different geologic settings.
Rio Maria greenstone belt
Figure 10a shows a gravity profile (solid line) across a green- stone belt consisting of metavolcano sedimentary rocks in the Rio Maria region of the state of Par ´a, Brazil. This unit was compressed by two blocks of granitoid rocks in a dextral trans- pression regimen (Souza et al., 1992). Density measurements of rock samples collected from outcrops indicate that the den- sity contrast between the metavolcano sedimentary and grani- toid rocks is about 0.3 g/cm^3. Figure 10a shows the fitted grav- ity anomaly (dashed line) produced by Souza et al. (1992) interpretation (Figure 10b) using interactive gravity modeling based on the expected synformal geometry for the greenstone belt unit and assigning uniform density contrasts of 0.3, 0.32, and 0.32 g/cm 3 to gravity sources A, B, and C, respectively. We inverted the same anomaly assuming an interpretation model consisting of a 72 × 64 grid with dimensions of 0.5 and 0.125 km in the x- and z-directions, respectively. We set μ = 0.05, f = 50 000, and τ = 0.05. The geometric elements ( e 1 – e 6 in Figure 11b) were defined so as to produce estimated gravity sources as close as possible to the interpretation given in Fig- ure 10b. All geometric elements were assigned a target den- sity contrast of 0.3 g/cm^3. The result (Figure 11b) shows that the proposed approach may lead to interpretations of multi- ple and complex gravity sources equivalent to the ones ob- tained by interactive modeling, but in a much easier and faster way and with the certainty of obtaining an acceptable fit to the data, as shown in Figure 11a (solid line).
Figure 10. Rio Maria greenstone belt. (a) Observed (solid line) and fitted (dashed line) Bouguer anomalies. (b) Interactive modeling (stippled polygons) according to Souza et al. (1992), assigning uniform density contrasts of 0.3, 0.32, and 0.32 g/cm^3 to gravity sources A, B, and C, respectively.
Cornubian batholith
Figure 12a shows the gravity anomaly (dots) produced by the Bodmin Moor pluton, part of the Cornubian batholith located in the county of Cornwall, England. The batholith has a granitic composition and intrudes low-grade, region- ally metamorphosed sediments and igneous rocks (Edmonds et al., 1975). Bott and Scott (1964) have modeled this pluton
Figure 11. Rio Maria greenstone belt. (a) Observed (dots) and fitted (solid line) Bouguer anomalies. (b) Inversion result us- ing our method (grayscale cells) and assigning a single target density contrast of 0.3 g/cm^3 to axes e 1 , e 4 , and e 6 (dashed lines) and to points e 2 , e 3 , and e 5 (dots).
Figure 12. Cornubian batholith. (a) Observed Bouguer anomaly (dots) and fitted anomalies; after Bott and Scott (1964) (dashed line) and using our approach (solid line). (b) Interactive modeling after Bott and Scott (1964), us- ing three polygons (solid line) with density contrasts of −0.16 g/cm^3 (southernmost), −0.13 g/cm^3 (intermediate), and −0.10 g/cm^3 (northernmost). The gray cells show the inversion result using our method, with geometric elements e 1 – e 5 indi- cated by dashed lines (axes) and dots (points). Elements e 1 – e 3 were assigned a target density contrast of −0.16 g/cm^3 , and el- ements e 4 – e 5 were assigned target density contrasts of −0. and −0.10 g/cm 3 , respectively.
J8 Silva and Barbosa
by incorporating the assumption that its density contrast in- creases to the north via a model consisting of three homo- geneous compartments with density contrasts of −0.16 g/cm^3 (the southernmost), −0.13 g/cm^3 (the intermediate), and −0.10 g/cm^3 (the northernmost). Bott and Scott’s (1964) in- terpretation and the corresponding fitted gravity anomaly are shown in Figures 12b (solid line) and Figure 12a (dashed line), respectively. We inverted this anomaly by assuming an interpretation model consisting of 58 × 24 cells with dimensions of 1.0 and 0.5 km in the x- and z-directions, respectively. We set μ = 0.5, f = 500 000, and τ = 0.01. The geometric elements ( e 1 to e 5 in Figure 12b) were introduced to produce an estimated gravity source close to Bott and Scott’s (1964) interpretation. Elements e 1 – e 3 have a target density contrast of −0.16 g/cm^3 , and elements e 4 and e 5 were assigned target density contrasts of −0.13 and −0.10 g/cm^3 , respectively. The result is shown in Figure 12b, which is very close to Bott and Scott’s (1964) inter- pretation but which displays a better anomaly fit (Figure 12a, solid line).
East Bull gabbro-anorthosite complex
Figure 13a displays the gravity anomaly (dots) over East Bull Lake, Ontario, Canada (Paterson and Reeves, 1985). It is caused by a gabbroic-anorthositic lopolith with the anorthositic rocks underlying the gabbroic rocks. Both units have a higher density contrast than the tonalitic country rocks. Automatic inversion methods applied to this anomaly present a severe limitation, caused by the strong vertical interference of the sources. We applied to this anomaly the robust procedure described in the Methodology section. First, we estimated the geome- try of the anorthositic unit by establishing an interpretation model consisting of 60 × 40 cells with dimensions of 0.1 km in both the x- and z-directions. We set μ = 1, f = 5000, and τ =
Figure 13. East Bull gabbro-anorthosite complex. (a) Ob- served (dots) and fitted (solid line) Bouguer anomalies. (b) Robust inversion result (gray cells) using a single axis ( e 1 in dashed line) for outlining the anorthositic body with a target density contrast of 0.2 g/cm^3.
0.01. The geometric element consists of a single horizontal axis at the surface, extending over the known anorthosite outcrops ( e 1 in Figure 13b). The assigned target density contrast to this axis was 0.2 g/cm^3. The value of μ was selected in such a way as to produce an interpreted geometry of the gabbroic rocks (see below), displaying a surface extent close to the known ex- tent of the gabbroic rocks. In this way, just a single inversion was executed in the robust procedure, and the anomaly pro- duced by the anorthosite was taken as the fitted anomaly. The rejected data points were completely determined in this case, just by the choice of μ. The estimated anorthositic unit and the
Figure 14. East Bull gabbro-anorthosite complex. (a) Residual anomaly (dots, obtained by subtracting the Bouguer anomaly from the fitted anomaly in Figure 13) and fitted anomaly (solid line). (b) Inversion result (gray cells) using a single axis ( e 2 in dashed line) for outlining the gabbroic body with a target density contrast of 0.08 g/cm^3.
Figure 15. East Bull gabbro-anorthosite complex. (a) Bouguer anomaly (dots) and fitted anomaly (solid line) resulting from the sum of the fitted anomalies displayed in Figure 13a and Figure 14a. (b) Inversion result (gray cells) obtained by sum- ming the inverted density contrast distributions displayed in Figure 13b and Figure 14b.
