Focusing geophysical inversion images

A critical problem in inversion of geophysical data is developing a stable inverse problem solution that can simultaneously resolve complicated geological structures. The traditional way to obtain a stable solution is based on maximum smoothness criteria. This approach, however, provides smoothed unfocused images of real geoelectrical structures. Recently, a new approach to reconstruction of images has been developed based on a total variational stabilizing functional. However, in geophysical applications it still produces distorted images. In this paper we develop a new technique to solve this problem which we call focusing inversion images. It is based on specially selected stabilizing functionals, called minimum gradient support (MGS) functionals, which minimize the area where strong model parameter variations and discontinuity occur. We demonstrate that the MGS functional, in combination with the penalization function, helps to generate clearer and more focused images for geological structures than conventional maximum smoothness or total variation functionals. The method has been successfully tested on synthetic models and applied to real gravity data.

GEOPHYSICS, VOI.. 64, NO. 3 (MAY-JUNE 1999); P. 874-887,15 FIGS. Focusing geophysical inversion images Oleg Portniaguine* and Michael S. Zhdanov* ABSTRACT A critical problem in inversion of geophysical data is developing a stable inverse problem solution that can si­multaneously resolve complicated geological structures. The traditional way to obtain a stable solution is based on maximum smoothness criteria. This approach, however, provides smoothed unfocused images of real geoelectri­cal structures. Recently, a new approach to reconstruc­tion of images has been developed based on a total varia­tional stabilizing functional. However, in geophysical ap­plications it still produces distorted images. In this paper we develop a new technique to solve this problem which we call focusing inversion images. It is based on specially selected stabilizing functionals, called minimum gradi­ent support (MGS) functionals, which minimize the area where strong model parameter variations and disconti­nuity occur. We demonstrate that the MGS functional, in combination with the penalization function, helps to generate clearer and more focused images for geologi­cal structures than conventional maximum smoothness or total variation functionals. The method has been suc­cessfully tested on synthetic models and applied to real gravity data. INTRODUCTION One of the critical problems in inversion of geophysical data is developing a stable inverse problem solution which at the same time can resolve complicated geological structures. Tra­ditional geophysical inversion methods are usually based on Tikhonov regularization theory, and they provide a stable so­lution of the inverse problem. This goal is reached, as a rule, by introducing a maximum smoothness stabilizing functional. The obtained solution provides a smooth image of real geo­electrical structures that sometimes makes it look geologically unrealistic. A new approach to reconstructing noisy images, developed in papers by Rudin et al. (1992) and by Vogel and Oman (1998), is based on a total variational stabilizing functional. The func­tional requires that the model parameter's distribution be of bounded variation. This requirement is much weaker than one of maximum smoothness because it can be applied even to dis­continuous functions. In this way, the total variation method produces better quality images for blocky structures. How­ever, it still decreases bounds of model parameter variation and therefore somehow distorts the real image. We study different ways of focusing geophysical images us­ing specially selected stabilizing functionals. In particular, we introduce a new stabilizing functional that minimizes the area where strong model parameter variations and discontinuity oc­cur. We call this new stabilizer a minimum gradient support (MGS) functional. We demonstrate that this MGS functional, in combination with the penalization function, helps to gener­ate a stable solution of the inverse problem for complex geo­logical structures and does not impose destructive restrictions on the bounds of model parameter variations. Thus, it helps to generate much more focused images than conventional meth­ods. We call this approach focusing inversion images. We also present a comparative analysis of inversion schemes based on different stabilizing functionals. This analysis shows that inversion codes based on the MGS stabilizing functional and the penalization function could be considered a good al­ternative to maximum smoothness or total variational-based inverse algorithms. TIKHONOV REGULARIZATION AND STABILIZING FUNCTIONALS Consider the inverse problem d = Am, (1) where A is the forward modeling operator; m = m(r), a scalar function describing geological model parameter distribution in some volume V in the earth (m e M, where M is a Hilbert space of models with L2 norm); and d = d(r), a geophysical data set (d e D, where D is a Hilbert space of data). Published on Geophysics Online, February 17,1999. Manuscript received by the Editor April 13,1998; revised manuscript received November 5, 1998. Published on line. *Dept. of Geology and Geophysics, University of Utah, Salt Lake City, Utah 84112-1183. E-mail: [e-mail]; mzhdanov@mines. utah.edu. © 1999 Society of Exploration Geophysicists. All rights reserved. 874 F o c u s in g In v e rs io n Im a g e s 875 Inverse problem (1) is usually ill posed, i.e., the solution can be nonunique and unstable. The conventional way of solving ill-posed inverse problems, according to regularization theory (Tikhonov and Arsenin, 1977; Zhdanov, 1993), is based on min­imization of the Tikhonov parametric functional: Pa(m) = 4>(m) + as(m), (2) where 4>(m) is a misfit functional determined as a norm of dif­ference between observed and predicted (theoretical) field, 4>(m) - ||Am - d\\2D = (Am - d, Am - d)i (3) Functional s(m) is a stabilizing functional (stabilizer). There are several common choices for a stabilizer. One is based on the least-squares criteria or, in other words, on an L2 norm for functions describing geoelectrical model parameters: sL2(m) = m = (m,m) = Jv dv - min. (4) The conventional argument in support of this norm comes from statistics and is based on the assumption that the least-squares image is the best over the entire ensemble of all possible images. Another stabilizer uses the minimum norm of the difference between the selected model and some a priori model mapr: SL2apr(m) = N ~ mapr\\2 = mill. (5) This criterion, as applied to the gradient of model parameters Vm, brings us to a maximum smoothness stabilizing functional: Sm*xsm(m) = ||Vm||2 = (Vm, Vm) = min. (6) It has been successfully used in many inversion schemes devel­oped for EM data interpretation (Berdichevsky and Zhdanov, 1984; Constable et al., 1987; Smith and Booker, 1991; Xiong and Kirsh, 1992; Zhdanov and Fang, 1996). This stabilizer pro­duces smooth geoelectrical models which in many practical situations do not describe properly the real blocky geological structures. It also can result in spurious oscillations when m is discontinuous. In a paper by Rudin et al. (1992), a total variation approach to reconstruction of noisy, blurred images is introduced. It uses a total variation stabilizing functional, which is essentially the Li norm of the gradient: sTv(m) = ||Vm||L L Vm I dv. 0) This criterion requires that the model parameter's distribution in some volume V be of bounded variation (Giusti, 1984). How­ever, this functional is not differentiable at zero. To avoid this difficulty, Acar and Vogel (1994) introduce a modified total variation stabilizing functional: sprv(m) JvV\v^2 ■ fi1 dv. (8) The advantage of this functional is that it does not require that the function m be continuous but that it be piecewise smooth (Vogel and Oman, 1998). The total variation norm does not penalize discontinuity in the model parameters, so we can re­move oscillations yet preserve sharp conductivity contrasts. At the same time it imposes a limit on the total variation of m and on the combined arc length of the curves along which m is discontinuous. That is why the functional produces a much better result than maximum smoothness functionals in imaging blocky structures. Total variation functionals sTV(m) and spTV(tn) tend to de­crease bounds of model parameter variation [equations (7) and (8)] and in this way still try to smooth the real image. However, this smoothness is much weaker than in the case of traditional stabilizers in equations (5) and (6). We can diminish this smoothness effect by introducing a new stabilizing functional that would minimize the area where sig­nificant variations of the model parameters and/or disconti­nuity occur. We call this new stabilizer a minimum gradient support (MGS) functional. For the sake of simplicity we first discuss a minimum support (MS) functional, which provides the model with the minimum area of the anomalous parame­ters distribution. MS AND MGS STABILIZING FUNCTIONALS The minimum support functional was considered first by Last and Kubik (1983), where the authors suggest seeking the source distribution with the minimum volume (compactness) to ex­plain the anomaly. This approach is modified in Guillen and Menichetti (1984) by introducing the functional minimizing moments of inertia with respect to the center of gravity or to a given axis. We consider an approach based on a minimum gradient support stabilizer which leads to the construction of models with sharp boundaries. Consider the following integral of model parameter distri­bution: n2 2T^2 dV' (9) v ml + p We introduce the support of m (denoted spt m) as the combined closed subdomains of V where m ^ 0. We call spt m a model parameter support. Then expression (9) can be modified as Mm) = [ Jv Jp(m) = f Jsx 1- P2 spt m = sptm + P2] dv P2 f J S[ ■ dv. spt m From expression (10) we can see that Jp(m) -> sptm, if fi -y 0. (10) (11) Thus, integral (m ) can be treated as a functional, proportional (for a small fi) to the model parameter support. We can use this integral to introduce a minimum support stabilizing functional spMs(m) as spMs(m) = Jp(m (m - mapr )2 lapr ^ Iv (m - mapr)2 ■ dv. (12) To justify this choice, we should prove that Sf)MS(m) can be considered as a stabilizer according to regularization theory. This fact is demonstrated in Appendix A. This functional has an important property: It minimizes the total area with the nonzero departure of the model parameters from the given a priori model. Thus, dispersed and smoothed distribution of the parameters with all values different from the a priori model mapr results in a big penalty function, while well-focused distribution with a small departure from mapr has876 P o rtn ia g u in e a n d Z h d a n o v a small penalty function.This approach was originally used by Last and Kubik (1983) for compact gravity inversion. We can use this property of a minimum support functional to increase the resolution of blocky structures. To do so we modify spMs(m) and introduce a minimum gradient support functional as f Vm•Vm s>"es{m) = J>pm] = <13) The value sptVm denotes the combined closed subdomains of V, where Vm ^ 0. We call sptVm a model parameter dis­tribution gradient support (or, simply, gradient support). Then expression (13) can be modified as SfiMGs(m) = sptVm - £2 f ---------- dv. (14) J sptVm Vm • Vm + pz From the last expression we can see that sPMGs(m) -+ sptVm, if 0. (15) Thus, we can see that s^Mas(m) can really be treated as a func­tional proportional (for a small fi) to the gradient support. A possible way to clarify the physical interpretation for the math­ematical form of equation (13) is to realize that the terms where gradient is nearing zero (or much less than fi) have zero con­tribution, while terms where any gradient exists (larger than ft) have contributions equal to one, even if the gradient is very large. Thus, solutions with sharp boundaries are promoted but the penalty for large gradients (discontinuous solutions) is not excessive. Repeating the considerations described in Appendix A for spMs(tn), one can demonstrate that the minimum gradient sup­port functional satisfies Tikhonov criteria for a stabilizer. PARAMETRIC FUNCTIONAL MINIMIZATION SCHEME Within the framework of the regularization theory, as dis­cussed above, the inverse problem solution is reduced to the minimization of the Tikhonov parametric functional Pa [equa­tion (2)], which can be written as Pa(m) = (Am - d, Am - d)r> + as(m). (16) Note that all stabilizing functionals introduced above can be written as the squared L2 norm of some function of the model parameters: s(m) = (f (m), f (m)). (17) For example, the maximum smoothness stabilizer appears if fmaxsm(m) - Vm. (18) In the case of the total variation stabilizing functional spTv(m), this function is equal to fpTv(m) = (|Vm|2 + £2)1/4. (19) In the case of the minimum support functional s,tMS(m), we obtain Finally, for the minimum gradient support functional spMGS(m), we find V Yt\ hMGs{m) = (Vm.Vm | /?-)•/-- (21) We can introduce a variable weighting function f(m) e(m) = ((m, m) + e2)1/2 ' (22) where s is a small number. Then the stabilizing functional in general cases can be written as the weighted least-squares norm of m: s(m) = (f(m), f (m)) « (we(m)m, we(m)m) = (m,m)We = I w^(m)m2 dv if s Jv 0. (23) f/3Ms(m) im ' + fi2')1!2 (20) The corresponding parametric functional can be written as Pa(m) - (Am - d, Am - d)D + a(m, m)We. (24) Therefore, the problem of the minimization of the paramet­ric functional introduced by equation (16) can be treated in a similar way as the minimization of the conventional Tikhonov functional with the L2 norm stabilizer sLiapr(m) [equation (5)]. The only difference is that now we introduce some a priori variable weighting functions we(m) for model parameters. This method is similar to the variable metric method; however, in our case the variable weighted metric is used to calculate the stabilizing functional only. The minimization problem for the parametric functional in­troduced by equation (24) can be solved using the ideas of traditional gradient type methods. The computational procedure to minimize the parametric functional (24) based on the reweighted conjugate gradient method is presented in Appendix B. PENALIZATION OF MATERIAL PROPERTY AND FOCUSING INVERSION In this section we discuss the possibility of using some ideas of the composite materials theory for solving the geophysi­cal inverse problem. Assume the geological model can be de­scribed as a composite of two materials with known physical properties (for example, density, magnetization, or electrical conductivity). One material corresponds to the background ho­mogeneous cross-section; the other one forms the anomalous body. In this situation, the values of the material property in the inversion image can be equal to the background value or to the anomalous value. However, the geometrical distribution of these values is unknown. We can force the inversion to produce a model which not only fits the data but which is also described by these known values, thus painting the geometry of the ob­ject. In the composite materials literature, this method is known as penalization. There is a simple and straightforward way of combining penalization and the MGS method. Numerical tests show that MGS generates a stable solution, but it tends to pro­duce the smallest possible anomalous domain. It also makes the image look unrealistically sharp. At the same time, the ma­terial property values m(r) outside this local domain tend to be equal to the background value which nicely reproduces first composite material, i.e., the background. We can impose the upper bound for the positive anomalous parameter valuesF o c u s in g In v e rs io n Im a g e s 877 mlb(T) (the; second material) and, during the iterative process, chop off all the values above this bound. This algorithm can be described as m (r) -mbg (r) = maub (r), if [m (r) - mbg (r)] > maub (r), . (25) m(r) - mbg(r) = 0, if [m(r) - mbg{r)] < 0. Thus, according to formula (25), the material property values m(r) are always distributed within the interval mhg{r) <m(r) < mbg(r) + maub(r). (26) A similar rule is applied in the case of negative anomalous parameter values. NUMERICAL COMPARISON We compare results of regularized inversion performed with the following stabilizers: maximum smoothness imaxJm(»i), the total variation functional sTV(m), the MS functional and the MGS functional spMGs(>Ti). We also consider a focus­ing inversion method that combines Tikhonov regularization with the MGS functional and penalization of material prop­erty. Minimization problems for all these cases were solved using a reweighted conjugate gradient optimization technique, discussed in Appendix B. We present synthetic examples of different geophysical data inversion. They include gravity field, stationary magnetic field, and EM field data. 2-D gravity data inversion Let us treat m as density distribution. In this case opera­tor A is a linear forward gravity operator. Figure la presents synthetic gravity data with 5% random noise (solid line) com­puted for a rectangular material body presented in Figure lb. Note that we have data in only ten observation points. The unknown densities in the grid shown in Figure lb form a large 20 x 15 matrix of unknown parameters. Thus, the inverse prob­lem is underdetermined, which can lead to multiple solutions. We have run four inversions with the different stabilizers and have obtained four different models, shown in Figure lc-f. The a) Data d) TV b) True model 0 50 100 c) Max smoothness 10 20 E 30 a D 50 60 70 11 0.8 0.6 0.4 0.2 50 100 e) Minimum Support 10 ?0 £ 30 -.40 Q> u 50 60 70 10 50 Distance, m 50 Distance, m Fig. 1. A 2-D gravity inversion for rectangular body. Grayscale shows normalized density. 878 P o rtn ia g u in e a n d Z h d a n o v theoretical data computed for these models fit the observed data practically with the same accuracy of 5% (all four pre­dicted data curves are shown by stars on Figure la). Figure lc shows the result of inversion with a maximum smoothness stabilizer. Figure Id shows the result obtained with a total variation stabilizer, which is better than the first one, but still the image is very dispersed. Figure le shows an inversion result with an MGS stabilizer. The image is oversharpened. Figure If presents the result of focusing inversion. For the fo­cusing inversion we assume we know the upper bound value of the anomalous density. The next set of inversions has been done for the model of two small bodies (Figure 2). Figure 2a depicts observed data with 5% random noise (solid line) and theoretical predicted fields for four inversion results (stars) shown in the other pan­els. Figure 2c shows the solution with the maximum smoothness functional, Figure 2d presents the bounded total variation so­lution, Figure 2e shows the solution with the minimum gradient support functional, and Figure 2f demonstrates the focused im­age. We again assume we know the upper bound value of the anomalous density. Figure 3 shows the set of equivalent solutions for steplike density distribution: (a) actual data with 5% random noise (solid line) and theoretical predicted data for four inversion results (stars) shown in the other panels, (b) actual model, (c) maximum smoothness solution, (d) the bounded total varia­tion solution, (e) the solution with a minimum gradient support functional, and (f) the result of focusing inversion. The focusing inversion produces the best image of the steplike structure. 2-D magnetic data inversion Now we assume that m is magnetic succeptibility and oper­ator A is the linear forward magnetic operator. We solve the stationary magnetic inverse problem. Figure 4 shows (a) syn­thetic observed magnetic data with 5% random noise and the­oretical predicted data for inversion results (stars) shown in Figure 4c, (b) the actual model, and (c) the bounded total variation inversion result. We now have two anomalous bod­ies with different susceptibilities. We first assume we know the anomalous property of both bodies. This knowledge is included in the algorithm as a priori information about the distribution a) Data d) TV 10 ■ 20 0.08 O CO UJ 0.06 Depth, (71 o o 0.04 60 0.02 70 0 50 100 10 . 20; e 30; f 40: 0) Q so; 6o; 70 ' b) True model e) Minimum Support 50 100 0.5 LJ 0 10 20 E 30 f 40 © Q 50 60 70 1.5 0.5 LJ 0 50 100 c) Max smoothness f) Focusing 10 20 E 30 f 40 a> ° 50 60 70 50 Distance, m 0.1 0.05 100 10 20 e 30 f 40 0 Q 50 60 70 50 Distance, m 0.5 LJ 0 100 Fig. 2. A 2-D gravity inversion for two small bodies. Grayscale shows normalized density. Focusing Inversion Images 879 of the constraints of anomalous susceptibility (Figure 5b). The resulting image is presented in Figure 5a. It resolves well the position and shape of both bodies. Figure 5d reflects the wrong assumption about the bodies' susceptibilities: one is two times bigger and the other is two times smaller. Figure 5c shows the focused image computed for this case. As one would expect, the sizes of the bodies correspondingly increase and decrease by two times. Figure 5f reflects the wrong a priori information about susceptibility: the susceptibility of the first body is two times smaller, and the susceptibility of the second body is two times larger than the true values. Figure 5e shows the corre­sponding focused image. This example suggests that even if we do not know the property exactly, focusing inversion still can be applied and can produce useful results. 3-D borehole induction data inversion We have applied different stabilizing functionals discussed above to solving the following EM inverse problem. Consider the model of two conductive bodies located at a depth of 1000 m (Figure 6). The bodies are prisms 20 x 20 m in the X and Y di­rections and 10 m in the Z direction. The observation array is formed by a vertical magnetic dipole transmitter and three- component magnetic field receivers located in the borehole with a vertical separation of 6 m. Resistivity of the bodies is 1 ohm-m, while the background resistivity is 1000 ohm-m. The theoretical frequency-domain EM field in this model was sim­ulated for frequencies of 16,32, 64, and 128 kHz using SYSEM integral equation forward-modeling code (Xiong, 1992). The transmitter-receiver installation was moving along the bore­hole from 950 m to 1050 m with observations every 10 m (stars in Figure 6). We use three-component data measured at the single obser­vation point in the borehole to obtain information about the location of the conducting bodies in the horizontal plane. The results of simulation are shown in Figure 7 (solid lines). The first experiment demonstrates the result of inversion with a total variation stabilizer: Prvim) = <p(m) + aspTv(m) = (Am - d, Am - d)o + f y/\Vm\2 + p2 dv - min. (27) Jv a) Data d) TV b) True model 10 20 £ 30 f 40 o ° 50 60 70 10 20 E 30 €"0 <D ° 50 60 70 50 100 e) Minimum Support I 40 20 -20 50 100 c) Max smoothness f) Focusing 10 20 E 30 f 40 (D Q 50 60 70 ■ 2 E 30 f 40 <D Q 50 1 60 70 1 ^ n 50 Distance, m 100 Fig. 3. A 2-D gravity inversion for step-like structure. Grayscale shows normalized density. 880 Portniaguine and Zhdanov The inversion image is shown in Figure 8. We cannot see two separate bodies in this figure. At the same time, the misfit between the observed and predicted data for this image is only 1.5%. The next numerical experiment demonstrates results of in­version using a minimum gradient support stabilizer, PMGs(m) = <Km) + aSgradsup(m) = (Am - d, Am - d)[) (Vm, Vm) L (Vm, Vm) + ft- dv, (28) and penalization, described above. Figure 9 presents the re­sults of inversion. Figure 7 shows the comparison between the a) Data 10 20 E 30 f 40 <D Q 50 60 70 10 20 £ 30 50 60 70 b) True model 1.5 0.5 50 100 150 Distance, m c) TV 200 0.4 0.3 0.2 0.1 0 50 100 150 Distance, m 200 Fig. 4. A 2-D magnetic field inversion. Grayscale shows normalized magnetization of the cells. observed (solid lines) and theoretical predicted data (stars) computed for the model shown in Figure 9. In this image, two bodies are obviously resolved. At the same time, the accuracy of fitting here is almost the same (1.5%) as that for total vari­ation inversion. The obtained results clearly demonstrate the advantages of the MGS plus penalization approach. FOCUSING INVERSION ON PENASQUITO GRAVITY DATA Gravity data for the Penasquito site, collected by Kennecott Exploration, are used as a test for inversion. The map of Bouguer anomalies for this site (terrain corrected for nearest 30 m) is shown in Figure 10a. The subsurface geology of the area is characterized by the presence of intrusions embedded into sedimentary formations. Core tests show lowest density for breccia and quartz porphyry samples (2.32-2.47 g/cm3). Density is 2.58-2.73 g/cm3 for both altered and unaltered background formations. Thus, negative gravity anomalies are possibly associated with breccia pipes. Most of the area is covered with alluvium; how­ever, one breccia pipe is outcropping at the central part of the map. Another breccia pipe was confirmed by drilling. In the inversion procedure, contrast for breccia and back­ground rock was taken as -0.3 g/cm3. However, there are many areas with positive gravity anomalies up to 1 mGal, which man­ifest formations with density higher than the background. Focusing inversion allows us to obtain a well-focused, sharp subsurface image, in contrast to widely known smooth in­version methods. It also requires application of the penaliza­tion technique, in which upper and lower limits of anomalous density variations are used to produce the density model. In this example we used values of -0.3 g/cm3 as lower limits which corresponded to the drilling core data about the brec­cia pipe's anomalous density. The positive constraints were taken as +0.3 g/cm3 to designate unknown high-density for­mations. The subsurface region under investigation was divided into cubic cells of 100 x 100 m horizontal size. Cell size increases with depth, starting with 50 m at the surface, then 50, 75,100, 150,200,300,400, and 500 m. There were 7800 data values and 20 000 model parameters. The density contrast within each cell was assumed constant, but it changed from cell to cell. Starting from the model with the zero anomalous density, the inversion procedure iteratively converged to the model that best fit the gravity data. Three separate experiments were done. First, focusing inver­sion was performed for the entire area. Second, minimum L2 norm (smooth) inversion was done for comparison, also for the entire area. Third, the inversion was applied to a local data set above one of known breccia pipes using a 3-D grid with small, uniform, cubic cells. The results of the third experiment were compared with drilling information available in the site. Results of inversion The focusing inversion was performed on the real grav­ity data from the Penasquito site. Only 10 minutes of com­puter time on a SPARC-20 were required to invert a relatively large 3-D model as demonstrated here. Figure 10b shows dataFocusing Inversion Images 881 predicted from the focusing inversion. The predicted data fit the observed data well. Figure 11a shows the residual field, which is the difference between observed and predicted grav­ity data. We can treat the residual field as the random noise which contaminates real data. Maximum errors are on the or­der of 0.1 mGal, but they occur only above one of the breccia pipes. Most of the errors are less than 0.02 mGal. A histogram of residuals is shown in Figure lib. The residuals form a Gaussian distribution, which is not surprising, given the fact that least- squares minimization of residuals was performed. Dispersion (square root of sum of error squares divided by number of samples) is 0.01 mGal for this plot. Most of the residual field is of short length, which means it represents random observation errors and near-surface features too small and shallow for the resolution of the method. The resulting inversion model is presented further as slices of anomalous density at different depths. Figure 12 presents slices at 200 and 325 m depth, respectively. The plot, corresponding to 200 m depth, is most informative and clear (Figure 12a). It shows two known breccia pipes at 0 N -0.5 E and -1 N 0.3 E. The prospective pipe at -0.3 N -2.2 E starts shifting to the north, and the prospective pipe at -1.5 N 2.2 E is shifting toward the south. There are also numerous positive-contrast density bodies. One of these features (at 0 N -1.2 E) is present on the deeper slice at 325 m depth (Figure 12b). a) 10 20 E 30 10 20 £ 30 f 40 <D Q 50 60 70 10 20 £ 30 f 40 <D Q 50 60 70 1.5 0.5 50 100 150 200 c) 1 50 100 e) 150 200 0.5 50 100 150 Distance, m 200 b) lIllllllllllllKlllll IIIIIIIIIIIIIIIIIIIII lllllillllllllilllll! IIIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIIII 1.5 0.5 50 100 150 200 10 20 £ 30 <D Q 50 60 70 1.5 0.5 50 100 150 200 Distance, m FlG. 5. A 2-D magnetic field inversion with various constraints. Grayscale shows normalized magnetization of the cells. Figure 5a,c, and e shows magnetic field inverse problem solutions for the model presented in Figure 4, with the different assumptions about the maximum anomalous susceptibility of each body (constraints). Figure 5b, d, and f shows distributions of the constraints applied in each case. Note that different constraints can be applied to different parts of the same model. For example, (b) shows that normalized susceptibility of the left body should not exceed 0.5 units and normalized susceptibility of the right body should not exceed 2. Figure 5a corresponds to a priori constraints shown in (b), (c) corresponds to (d), and (e) corresponds to (f).882 Portniaguine and Zhdanov Fig. 6. True model for 3-D borehole EM inversion. Grayscale shows anomalous conductivity in siemens/meter. Minimum L2 norm inversion results For comparison, the results of the minimum L2 norm inver­sion are also presented. This inversion produces smooth, murky images. However, it may also provide useful information. Smooth inversion produces the data which fit the observa­tion with almost the same accuracy as for focusing inversion. However, the inverse anomalous density model is different because it is a smooth model now. This result corresponds to the fact that the solution of gravity inverse problem is nonunique. By introducing a different stabilizing functional in the Tikhonov regularization scheme, we select different solu­tions from the class of possible inverse models. The resulting smooth inversion model is presented further as slices of anomalous density at different depths (Figures 13a,b). The slices look somewhat similar to the corresponding sharp pictures; however, the images here are more dispersed and unclear. It is hard to evaluate the shape of anomalous bodies from these pictures, and some bodies cannot be distinguished at all. Validation of results with drilling data There are no data so far to confirm or reject any hypothesis about deep structures; however, there are drilling data avail­able on the site to comparc with at depth up to 100-200 m. One known breccia pipe was selected for more careful inver­sion in the small window to better understand geometry of this particular pipe and to check the reliability of inversion results. A window 1.5 x 1.5 km was cut from the data, and inver­sion was performed for the data within this window. Cells were taken as cubes with sides of 100 m for all depths, up to 1.5 km. After focusing inversion, the cells with zero density were erased and a 3-D image of the body was generated, as shown in Fig­ure 14. Stars in Figure 14 show boreholes. The X axis is directed to the east, and the Y axis is directed to the north. Comparing these pictures with the images of the same body obtained from the entire map, we cannot see any significant difference, which demonstrates the robustness of the algorithm to the cell's size. Figure 15 shows the view of the body from the top and the contour of breccia pipe derived from drilling data. Inversion correctly predicts which wells are inside and which wells are outside of the breccia pipe. CONCLUSIONS The results of our work demonstrate that by choosing differ­ent types of stabilizing functionals we can generate inversion images resolving the anomalous bodies with different accu­racy. The maximum smoothness functional obviously produces a very diffuse image. The total variation functional generates a more focused image but still cannot resolve anomalous bod­ies well. Finally, the MGS functional in combination with pe­nalization produces the more resolved and focused image of anomalous structures. Thus, the MGS functional in combination with penalization helps to generate clearer and more focused images for geolog­ical structures than conventional maximum smoothness and total variation functionals. Focusing inversion code was performed on the real gravity data from the Penasquito site. The results of focusing inversion have been compared with conventional inversion and checked against drilling data. Comparison shows that focusing inversion produces a different kind of information than the conventional smooth method. The shape and size of the bodies are much better resolved, especially at smaller depths, and are confirmed by drilling data. As such, focusing inversion can be a useful tool in interpreting the data. ACKNOWLEDGMENTS We thank Mr. J. Inman from Rio Tinto for useful discussion and comments. We are grateful to Rio Tinto for permission to publish the Penasquito gravity data and associated results. The authors acknowledge the support of the University of Utah Consortium of Electromagnetic Modeling and Inversion (CEMI), which includes Advanced Power Technologies, BHP, INCO, Japan National Oil Corp., MIM Exploration, Min- deco, Naval Research Laboratory, Newmont Exploration, Rio Tinto, Shell International Exploration and Production, Schlumberger-Doll Research, Western Atlas, Western Mining, Unocal Geothermal Corp., and Zonge Engineering. REFERENCES Acar, R., and Vogel, C. R., 1994, Analysis of total variation penalty methods: Inverse Problems, 10,1217-1229. Berdichevsky, M. N., and Zhdanov, M. S., 1984, Advanced theory of deep geomagnetic sounding: Elsevier Science Publ. Co., Inc. Constable, S. C., Parker, R. C., and Constable, G. G., 1987, Occam's inversion: A practical algorithm for generating smooth models from EM sounding data: Geophysics, 52, 289-300. Eckhart, U., 1980, Weber's problem and Weiszfeld's algorithm in gen­eral spaces: Math. Program., 18,186-196. Giusti, E., 1984, Minimal surfaces and functions of bounded variations: Birkhauser-Verlag. Guillen, A., and Menichetti, V., 1984, Gravity and magnetic inver­sion with minimization of a specific functional: Geophysics, 49, 1354-1360. Last, B. I, and Kubik, K., 1983, Compact gravity inversion: Geophysics, 48,713-721. Rudin, L. I., Osher, S., and Fatemi, E., 1992, Nonlinear total variation based noise removal algorithms: Physica D, 60, 259-268. Smith, J. T., and Booker, J. R., 1991, Rapid inversion of two- and three­dimensional magnetotelluric data: J. Geophys. Res., 96, 3905-3922. Tikhonov, A. N., and Arsenin, V. Y., 1977, Solution of ill-posed problems: V. H. Winston and Sons.Focusing Inversion Images 333 Fig. 7. Predicted (stars) and observed (solid lines) normalized magnetic field components (Hx, Hy, Hz) for 3-D borehole EM inversion. 150s 4 200 v 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 Fig. 8. Total variation inversion result. Grayscale shows anomalous conductivity in siemens/meter. Fig. 9. Focusing inversion result. Grayscale shows anomalous conductivity in siemens/meter. Vogel, C. R., and Oman, M. E., 1998, Fast total variation based recon­struction of noisy, blurred images: IEEE Trans. Image Processing, 7, 813-824. Xiong, Z., 1992, EM modeling of three-dimensional structures by the method of system iteration using integral equations: Geophysics, 57,1556-1561. Xiong, Z., and Kirsch, A., 1992, Three-dimensional earth conductivity inversion: J. Comp. Appl. Math., 42 ,109-121. Zhdanov, M. S., 1993, Regularization in inversion theory: Tutorial, Colorado School of Mines. Zhdanov, M. S., and S. Fang, 1996, 3-D quasi linear electromagnetic inversion: Radio Science, 31, No. 4,741-754.Fig. 10. -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Easting, km Predicted (a) and observed data (b) for focusing inversion milligals. . Grayscale shows a) 1 0.5 i ; o L -1.5 1.5 c^oo,^•.»" Uf\- J*-.S -2 -1.5 -1 -0.5 0 0.5 Easting, km 1.5 b) 1500 </> 1000 500 0.1 0.05 -0.06 -0.04 -0.02 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Residual value. mGal Fig. 11. Residual (a) and residual distribution (b) for focusing inversion. Grayscale shows milligals. Portniaguine and Zhdanovb) 1.5 0.5 I 0 t i -0.5 -1 -1.5 ■ . - - - - -- ■ 1 ■ ■ ■ I 1 ■ - -1 0. ).5 1 1.5 2 Easting, km 0.3 0.2 0.1 -0.1 -0.2 -0.3 Fig. 12. Slices of focusing inversion result. Gra timeter. (a) Slice at 200 m suit. Grayscale shows grams pi depth; (b) slice at 325 m depth. er cubic cen- a) -2 -1.5 -1 -0.5 0 0.5 Easting, km 1.5 I0.2 0.1 -0.1 -0.2 -0.3 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Easting, km Fig. 13. Slices of smooth image. Grayscale shows grams per cubic centimeter, (a) Slice at 200 m depth; (b) Slice at 325 m depth. Focusing Inversion Images 885886 Portniaguine and Zhdanov Fig. 14. A 3-D image of breccia pipe, viewed from the south­west. Stars show borehole locations. Gray cubes show cells with anomalous density. * \ " I / m . =e!. .. 1 ■ ; . I / 4 / 0.2 0.4 X km Fig. 15. Top view of 3-D image of breccia pipe. Gray cubes show cells with anomalous density. Stars show borehole loca­tions. Dashed line shows contour of breccia pipe inferred from drilling data. APPENDIX A MINIMUM SUPPORT FUNCTIONAL AS A STABILIZER According to regularization theory (Tikhonov and Arsenin, 1977), a nonnegative functional s(m) on some Hilbert space M is called a stabilizing functional if for any real c > 0 the subset Mc of elements m e M for which s(m) < c is a compact. Consider the subset Mc elements of M satisfying the condition SfiMs(m) < c, (A-1) where spMs(>n) is a minimum support stabilizing functional de­termined by equation (14). From the other side, spMs is a mono- tonically increasing function of ||m - mapr ||2 : s/3Gs(mi) < spGS(m2) if \\m\ - mapr\\ < ||m2 - mapr\\. (A-2) To prove this, let us consider the first variation of the minimum support functional: Ss0Ms(m) = S L (m - mapr)2 (m ■ m apr )2 + P: Iv ((m h2s Jv p2 = I a S(m 'v r)2+P2)2 r)2dv, dv S(n lapr )2 dv (A-3) where a2 = ((m - mapr)2 + p2)2' (A-4) Using the theorem of the average value, we obtain ^ r Ss0MS(m) = a2 J S(m-mapr)2dv = a2S I (m - mapr)2 dv = a2S\\m - m Jv (A-5) apr II (A-6) 2a2\\m - mapr\\S\\m - mapr\\, /olurne ' apr ii > 0, we obtain equa- where a2 is an average value of a2 in volume V. Taking into account that a2 > 0 and ||m - m tion (A-2) from equation (A-6). Thus, from condition (A-1) we see that IIm - mapr || <q,me Mc, (A-7) where q > 0 is some constant, i.e., M, forms a sphere in the space M with a center at the point mapr. It is well known that the sphere is a compact in a Hilbert space. Therefore, we can use functional spMS(m) as a stabilizer in a Tikhonov regularization process.Focusing Inversion Images 887 APPENDIX B CONJUGATE GRADIENT REWEIGHTED OPTIMIZATION Consider a discrete inverse problem equivalent to the gen­eral inverse problem (1) in the case of the discrete model parameters and data. Suppose that N measurements are per­formed in some geophysical experiment. Then we can consider these values as the components of vector d of a length N. Sim­ilarly, some model parameters can be represented as the com­ponents of vector m of a length L. In this case, equation (1) can be rewritten in matrix notation: d = A(m), (B-l) where A is the matrix column of the operator A. The parametric functional (24) for a general nonlinear in­verse problem can be expressed using matrix notations: P"(rii) = (WrfA(m) - W,zd)*(WdA(m) - W,d) + a(Wcm)*Weih, (B-2) where YVj and W(, are weighting matrices of data and model parameters and the asterisk denotes a transposed complex con­jugate matrix. We = Wc(m) is the matrix of the weighting func­tion we(m) introduced above in equation (22). Thus, using as /(m) in equation (22) the corresponding expression (18), we obtain a maximum smoothness stabilizer. Determining /(rii) according to equation (19) yields a total variation stabilizer. Substituting corresponding formula (20) or (21) instead /(m) produces minimum gradient support or minimum support stabilizers. We use the conjugate gradient method to minimize the para­metric functional (B-2). It is based on the successive line search in the conjugate gradient direction 7(m„): ift„+1 =m„+Jm = m„- knl(mn). (B-3) The idea of the line search is that we present Pa (ih„ - fc7(m„)) as a function of one variable k and, evaluating it three times along direction 7(ih„), approximately fit it by parabola and then find its minimum and the value of kn, corresponding to this min­imum. The conjugate gradient directions 7(m„) are selected as follows. First, we use the gradient direction 7(mo) = l(nio) = F;W^(A(m0) - d) + «W;0m0, (B-4) where Wj0 = YV;'(riv,). Next, the conjugate gradient direction is the linear combina­tion of the gradient on this step and the direction 7(m0) on the previous step: 7(mi) = l(iiii) + PiK™o)- (B-5) On the n + 1th step, 7(ih„+i) = l(ihn+i) + j6„+i/(m„), (B-6) where i(m„) = F;W5(A(m„) - a) + (B-7) and W^„ = W2(m„). The coefficients /3„+1 are defined from the condition that the directions 7(m„+i) and 7(m„) are conjugate: Pn+l r(mn+1)i(mn+1) i*(mn)i(rii„) (B-8) This algorithm always converges to the local minimum because on every iteration we apply the parabolic line search. We call this algorithm conjugate gradient reweighted optimization be­cause the weighting matrix W2„ is updated on every iteration. One can find the formal proof of the convergence of this type of optimization technique in Eckhart, (1980). In the case of linear forward operator A, the parametric func­tional has only one local minimum, so the minimization of Pa is unique (Tikhonov and Arsenin, 1977). The advantage of the conjugate gradient reweighted opti­mization algorithm is that we do not have to know the gradient of /(m) for every iteration-only its value for corresponding model parameters, which is easy to calculate.