Novel approach to the model appraisal and resolution analysis of regularized geophysical inversion

The existing techniques for appraisal of geophysical inverse images are based on calculating the model resolution and the model covariance matrices. In some applications, however, it becomes desirable to evaluate the upper bounds of the variations in the solution of the inverse problem. It is possible to use the Cauchy inequality for the regularized least-squares inversion to quantify the ability of an experiment to discriminate between two similar models in the presence of noise in the data. We present a new method for resolution analysis based on evaluating the spatial distribution of the upper bounds of the model variations and introduce a new characteristic of geophysical inversion, resolution density, as an inverse of these upper bounds.We derive an efficient numerical technique to compute the resolution density based on the spectral Lanczos decomposition method (SLDM). The methodology was tested on 3D synthetic linear and nonlinear electromagnetic (EM) data inversions, and also to interpret the helicopter-borne EM data collected by INCO Exploration in the Voisey's Bay area of Canada.

GEOPHYSICS, VOL. 71, NO. 6 (NOVEMBER-DECEMBER 2006); P. R79-R90,13 FIGS. 10.1190/1.2336347 A novel approach to the model appraisal and resolution analysis of regularized geophysical inversion Michael S. Zhdanov1 and Ekaterina Tolstaya1 ABSTRACT The existing techniques for appraisal of geophysical in­verse images are based on calculating the model resolution and the model covariance matrices. In some applications, however, it becomes desirable to evaluate the upper bounds of the variations in the solution of the inverse problem. It is possible to use the Cauchy inequality for the regularized least-squares inversion to quantify the ability of an experi­ment to discriminate between two similar models in the pres­ence of noise in the data. We present a new method for resolu­tion analysis based on evaluating the spatial distribution of the upper bounds of the model variations and introduce a new characteristic of geophysical inversion, resolution density, as an inverse of these upper bounds. We derive an efficient nu­merical technique to compute the resolution density based on the spectral Lanczos decomposition method (SLDM). The methodology was tested on 3D synthetic linear and nonlinear electromagnetic (EM) data inversions, and also to interpret the helicopter-borne EM data collected by INCO Exploration in the Voisey's Bay area of Canada. INTRODUCTION The results of geophysical data interpretation are usually present­ed in the form of a corresponding model of the earth's formations in the area of investigation. We determine this model by solving the in­verse problem for geophysical data, which are contaminated by noise and are acquired at a limited number of observation points. Be­cause of the ill-posed nature of inverse geophysical problems, the solutions are ambiguous and unstable. There are always many solu­tions that will tit the observed noisy data practically with the same data misfit. The variations of the inverse model parameters may be unreasonably large if we do not use regularization (Tikhonov and Arsenin. 1977). There are still significant uncertainties, however. even with the application of regularization to the inverse-problem solution. The question arises. What are the maximum possible varia­tions of the model parameters that would preserve the variation of the predicted data within the level of the noise in the observations? In other words, what is the practical resolution of the regularized inver­sion? This is one of the most important problems of exploration geo­physics. This problem arises in the initial stage of a geophysical in­vestigation when we design the geophysical survey. The same prob­lem appears at the final stage when we examine the results of the interpretation of the observed geophysical data. Actually, the ques­tion about sensitivity and resolution of the given geophysical meth­od is usually the first one asked by geologists working with geophys­ical data. The sensitivity of the geophysical method is determined as the ratio of the variation of the data to the variation of the model parame­ters. The sensitivity can be found by direct modeling of the theoreti­cal response for the given model perturbation, or by using a reciproc­ity principle (Rodi. 1976; McGillivray and Oldenburg. 1990; McGillivray et al.. 1994; Spies and Habashy. 1995; Zhdanov. 2002). The word resolution was introduced into geophysical inversion by Backus and Gilbert in their classic 1967 and 1968 papers about analysis of the general resolution power of the corresponding geo­physical method. In this paper, we are interested, instead, in the reso­lution study of a specific, regularized inversion of given geophysical data. In this sense, our approach provides the model appraisal of the regularized inversion. However, for this analysis, we use a mathe­matical technique, which is quite different from the ones discussed in previous publications (see. for example. Ramirez et al.. 1995; Alumbaugh and Newman. 2000). The existing techniques for appraisal of geophysical inverse im­ages are based primarily on calculating of the data and model resolu­tion and covariance matrices (Tarantola. 1987; Menke. 1989; Alum­baugh and Newman. 2000). These matrices make possible the a pos­teriori appraisal of the quality of the geophysical inversion by dis­playing a distribution of the variances of the model parameter m. which describes a standard deviation of the model parameters from Manuscript received by the Bditor August 31,2004: revised manuscript received April 9,2006: published online September 26,2006. 1 University of Utah, Department of Geology and Geophysics, Salt Lake City, Utah. B-mail: [e-mail]. © 2006 Society of Bxploration Geophysicists. All rights reserved. R79 Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ R80 Zhdanov and Tolstaya the inversion result. In geophysical applications, however, it may be useful to also estimate the upper bounds of the variations in the solu­tion of the inverse problem for the given errors in the observed data. These upper bounds of the model variations confine the actual reso­lution of the geophysical inversion. In this paper, we introduce a new characteristic of geophysical inversion, resolution density, which is determined as the inverse of the upper bounds of the model parame­ter variations, and develop a method for solving this problem. There have been previous examples of inversion procedures for generating the bounds on variables [see, for example, Parker (1975) on the theory of ideal bodies for gravity interpretation, Sabatier (1977a, 1977b, 1977c) on linear inverse problems with constraints, Oldenburg (1983) on funnel functions, and Stark et al. (1986) and Stark and Parker (1987) on tau-p inversions for seismic data inter­pretation]. However, the previous publications were focused mostly on constructing all kinds of extreme solutions for a specific geophys­ical problem, e.g., on finding the smallest envelope containing all ve­locity profiles consistent with the seismic data (Stark et al., 1986) or looking for an ideal body as one whose supremum is the smallest of all suprema of all solutions of the gravity inverse problem (Parker, 1975). Here, we consider the problem of evaluating the spatial distribu­tion of the upper bounds of the model parameter variations for the given inversion result. We introduce a novel approach to solving this problem using the Cauchy inequality for the regularized least- squares inversion. In the framework of this approach, we develop a method of resolution analysis for both the linear and nonlinear in­verse problems. We also develop a numerical method of resolution- density computation based on the spectral Lanzcos decomposition method (SLDM), which provides an efficient way of solving this problem for different values of the regularization parameter a (Zh­danov, 2002). The method is illustrated by the resolution study of 3D electro­magnetic (EM) inversions of airborne and magnetotelluric (MT) data. The case history includes interpretation of the helicopter-borne EM data collected by INCO Exploration in the Voisey's Bay area of Canada. We believe this new technique provides a useful tool for the analysis of the robustness of geophysical inversion. RESOLUTION OF GEOPHYSICAL INVERSION Astrict mathematical definition of the resolution of a geophysical method was introduced in Dmitriev et al. (1990); see also Zhdanov (2002, p. 31). According to this definition, the measure of the resolu­tion R of the given geophysical method is determined as the inverse of the norm of the inverse operator: where A is a linear forward-modeling operator for the given geo­physical problem. This definition comes from the equality 4nax = M<?=|, (2) K where 4milx is the maximum possible error in the solution of the in­verse problem for the given level of errors in the observed data S. Based on the last equations, one can say that two models, in, and m2, can be resolved if the following condition is satisfied: linii- m2|j > 4max = A The smaller the norm of the inverse operator, the larger the resolu­tion K and the closer to each other are the models that can be re­solved. If the inverse operator/^1 is not bounded, i.e., its norm goes to infinity, the resolution goes to zero, K = 0, and the maximum pos­sible variations in the determination of m are infinitely large. This situation appears in the case of ill-posed problems (Zhdanov, 2002). Note, however, that the aforementioned definition provides a global estimate of the resolution in the sense that we can estimate only a norm of the difference between two models that must be resolved. At the same time, it would be very important to be able to compute a lo­cal estimate of the resolution (resolution density), which would de­li ver a distribution of the upper bounds of the model parameter varia­tions in the regularized solution of the inverse problem for the given level of the errors in the observed data. Many papers analyze the effect of errors on the geophysical in­verse problem solution (e.g., Sabatier, 1977a, b, c; Stark et al., 1986; Stark and Parker, 1987; Menke, 1989). Generally speaking, there are two major points of view in addressing this problem: 1) The algebraic (deterministic) point of view [dating back to Lanczos (1961), Marquardt (1963, 1970), Backus and Gilbert (1967,1968), Backus (1970a, b, c), and Tikhonov and Arsenin (1977)]. 2) The stochastic (probabilistic) point of view [formulated in the pioneering papers of Foster (1961), Franklin (1970), Jackson (1972), Tarantola and Valette (1982), and Tarantola (1987)]. The stochastic point of view is widely used in geophysical litera­ture because it is closely associated with the statistical nature of noise in geophysical data (see Sambridge and Mosegaard, 2002). We would like to recall, however, Sabatier's remark (1977a, p. 125), "if one trusts a certain statistical interpretation of errors, and ergodicity, the solutions can be classed according to one's degree of confi­dence." At the same time, it has been demonstrated in many publica­tions [e.g., the classic work by Sabatier (1977a) that both points of viewresult in similarcomputational algorithms]. We analyze the solution of the geophysical inverse problem based on Tikhonov regularization, which corresponds to the algebraic (de­terministic) point of view (Tikhonov and Arsenin, 1977). RESOLUTION DENSITY Let us consider a linear matrix equation: d = Am. (3) Here, m is the vector of the model parameters of order N,„, d is the vector of the observed geophysical data of order Nj, and the matrix A is the Nj X N„, matrix of the linear forward-modeling operator. In the framework of Tikhonov regularization theory, the regular­ized solution of this inverse problem can be based on the parametric- functional minimization: Pa(m,d) = (WrfAm - Wf/d)*'(Wf/Am - Wf/d) + or(Wmm - Wmmaprf (Wmm - Wmmapr) = min, (4) Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ Resolution analysis o f geophysical inversion R81 where W,, and W,„ are some weighting matrices of the data and mod­el parameters; mllpr is some a priori model; * denotes the complex conjugate transpose matrix; and a is a regularization parameter. The detailed description of the optimal weighting matrices selection is given in Zhdanov (2002). A solution of the general least-squares problem of equation 4 is given by the following equations (Tikhonov and Arsenin, 1977): ma = (A*W;jA + Q-Wf;r,(A*W5d + «W,>upr). (5) Letus apply the variational operator iSto both sides of equation 5: <5ma. = (A*W5 A + aW,;)"'A*W^d. (6) We shall call the matrix R„ = (A*W;j A + aW-r1 A*W;j (7) a regularized inverse matrix. It is measured in the following units: [units of m] [units of RJ [units of d] The spatial variations of the resolution of the geophysical inver­sion can be found by individually analyzing the columns of matrix R„. Indeed, equation 6 in scalar notation can be written as Situ where KaiJ are the scalar components of R„, and S m, and Srf; are the components of vectors iSm„ and iSd, respectively. From the Cauchy inequality, we have mj where n, 2 and mj I -1/2 (8) (9) is a norm of the relative errors in the data. Note that the term represents a sum of the squares of the scalar components located in the /th column of R„. We can introduce a diagonal matrix 7?. formed by the elements 7i = 1,2, ..Nm. According to equa­tion 9, this matrix is related to the regularized inverse matrix R„ by the equation U = |idr>[diag(RX)]' -1/2 (10) We will call the diagonal matrix 7?. a resolution density matrix. Note that the resolution density matrix is computed based on the regularized inverse matrix R„. The last matrix, according to equation 7, is not a function of the data but of the operator of the forward prob­lem, the data and model parameters weights WJ and W2 used in the inversion, and the regularization parameter a. As a result, the resolu­tion density depends on the physics of the method used to collect data, source-receiverconfiguration, data components, coverage, etc. Therefore, similar to data and model resolution and covariance ma­trices (Menke, 1989), this matrix is affected by the geophysical methods under consideration and data acquisition parameters only. Now we can determine the upper bounds of the variations in the solution of the inverse problem for the given relative errors in the ob­served data, equal to e: sup \Snij I IlfSdM (11) Based on the last equations, we can determine the resolution den­sity of the inverse problem solution. Two models, mn) and m'2), in the vicinity of the point m0 can be resolved if the following condition is satisfied: ,<2>l irax ,p • (12) Thus, the upper bounds of the variations in the solution of the inverse problem are proportional to the relative errors in the data e. Note that the noise in the data affects the upper bounds of the model variations only; it does not affect the resolution density matrix introduced above. At the same time, equation 12 provides an appraisal of the in­verse problem solution by taking into account the errors (noise) in the data, the physics of the method used to collect the data, and the data acquisition parameters. The value 7?., is also related to the weights of the data and model and the regularization parameter a used in the inversion algorithm. The value 7?., is the measure of the resolution density for the given inverse problem solution. It is measured in the following units: Resolution density; [units of 7Zj\ 1 [units of m] (13) The larger the resolution density 7the closer to each other are the models that can be resolved. The low-resolution density 7?., corre­sponds to the area where even very different models cannot be re­solved. Note that both the maximum possible variations 4, milx and the resolution density 7?., depend on the cell number i. Thus, they de­scribe the spatial distribution of the variations in the inverse-prob­lem solution and of the resolution. By knowing the distribution of the resolution density in the area of inversion, we can identify the parts of the inverse model that are well resolved and the parts that are poorly resolved. Therefore, both the upper bound of the model varia­tions and the resolution density provide effective tools for a posterio­ri appraisal of the regularized inversion. It is important to notice that our model appraisal and resolution analysis method takes into ac­count the errors (noise) in the data, the physics of the method used to collect the data, the data acquisition parameters, and the inversion al­gorithm. Numerical calculation of the resolution density is a very challeng­ing computational problem. In Appendix A, we present an efficient algorithm for determining this parameter based on the SLDM (Lanc- zos, 1961; Druskin and Knizhnerman, 1994; Golub and Van Loan, 1996; Druskin et al., 1999; Zhdanov, 2002). The advantage of the SLDM method is that it allows us to compute the resolution density for all possible values of the regularization parameter a (Zhdanov, 2002). Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ R82 Zhdanov and Tolstaya RESOLUTION OF THE NONLINEAR INVERSE PROBLEM We presented a method for the resolution analysis of the linear in­verse problem. A similar technique can be introduced for nonlinear inverse-problem solution as well. Consider a nonlinear matrix equa­tion d = B(m), (14) where B is the nonlinear forward operator. Let us assume that m0 is a solution of the nonlinear inverse problem obtained by some inver­sion method. Let us perturb equation 14 in the vicinity of point m0. <?d = F0<?m, (15) where S m and S d are the perturbations of the model parameters and the data, respectively, and F0 is the Frcchet derivative matrix at m0. Our goal is to find out how the errors in the observed data (Sd will be transformed in the errors of the inverse-problem solution i5m. To solve this problem, we consider the regularized solution of equation 15. which can be expressed in the form (Zhdanov. 2002) Sm„ = (FqW^Fq + aW^r'FoW^d. We call matrix R„(m0). R„(m0) = (FXFo + «W;,)"1F;;Wfi (16) (17) a regularized inverse matrix of the nonlinear inverse problem at a point m0. Similar to the linear case, one can analyze, based on R„( m0). the maximum possible errors in model parameter distribution for the given level of errors in the observed data: I(5mJ" -z- 8" nl where ft 0; l|d||22 l^ra/(m0) ,i=i (18) (19) and e is a level of relative errors in the observed data: s = ||(5d||/||d||. The value 7?.0, is the measure of the resolution density for the given nonlinear inverse-problem solution m0. Figure 1. Schematic view of a conductive cubic model located with­in a resistive half-space and of a model HEM survey comprising five flight lines. Inequality 18 allows us to determine the maximum possible errors in the solution of the nonlinear inverse problem for the given errors in the observed data, using the equation similar to expression 11. The numerical technique for finding 7?.0, is similar to one devel­oped for a linear inverse problem in Appendix A. The only differ­ence is that we must substitute the linear forward-modeling operator by the Frcchet derivative in the corresponding equations of resolu­tion analysis. RESOLUTION STUDY OF THE LINEARIZED INVERSION OF A 3D SYNTHETIC HELICOPTER-BORNE EM SURVEY In this section, we illustrate the developed method of the resolu­tion analysis for the linearized airborne EM data inversion. Helicop­ter- borne EM (HEM) surveys are widely used in mineral explora­tion. The main difficulties in the modeling and interpretation of HEM data are related to the fact that for any new observation point, one must solve the forward problem anew for the corresponding position of the moving transmitter. In this situation, even forward modeling of HEM data over inhomogeneous structures requires an enormous number of computations. That is why. until recently, the interpretation of HEM data was restricted to simple 1D inversion. Zhdanov and Tartaras (2002) developed a new approach to the mod­eling and inversion of multisource array EM data based on the so- called localized quasi-linear (LQL) approximation. In the frame­work of this approach, forward modeling and inversion of multi­source data can be calculated at the same time for all different posi­tions of the transmitters. The LQL approximation also reduces the HEM data inversion to the solution of the linear inverse problem, which makes it possible to implement the linear resolution analysis developed in the previous sections of the paper. In the original paper by Zhdanov and Tartaras (2002). the linear EM inverse problem was solved using the conjugate gradient (CG) method. Zhdanov and Chernyavskiy (2004) introduced a new tech­nique for fast LQL inversion that employs the SLDM method (Druskin and Knizhnerman. 1994; Golub and Van Loan. 1996; Druskin et al.. 1999; Zhdanov. 2002). This technique helps to accel­erate HEM data inversion and provides a stable image of the geo­electrical target. We use a similar technique for the resolution analy­sis as well (see Appendix A). First, we consider a synthetic example of the resolution analysis of the HEM data inversion. We apply the integral equation software SYSEM (Xiong. 1992) to simulate such a survey over a relatively conductive (200ohm-m) cubic body located in a resistive (5000 ohm-m) half-space. Figure 1 depicts a 3D view of the model. Five lines were flown over the target at an altitude of 30 m and at a distance of 25 m apart. A schematic 3D view of the survey is shown in Figure 1. The moving transmitter-receiver system was a pair of vertical magnetic dipoles (simulating a horizontal coplanar coil pair) and a pair of horizontal magnetic dipoles (simulating a vertical coaxial coil pair) with 8 m of horizontal separation. The yy (coaxial) and zz (coplanar) components of the anomalous magnetic field were mea­sured every 15 m along the lines (50 observation points in each line). A 7.2-kHz frequency was used. We added 1 % random noise to the anomalous magnetic field and then inverted it using the SLDM method. The area of inversion, cen­tered around the body, was 150 m X 150 m X 150 m and was divid­ed into 12X 12X 12 cells. Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ Resolution analysis o f geophysical inversion R83 Figure 2 shows the vertical cross sections (along the.v-axis) of the 3D model obtained as a result of the regularized inversion with the minimum norm stabilizer (Zhdanov and Chernyavskiy. 2004). We should note that the inversion provides a correct position of the tar­get but underestimates the true conductivity of the body. Indeed, the recovered resistivity for the body is about 1000 ohm-m. while the true resistivity of the body is 200 ohm-m. This result comes without any surprise because it is well known that the linearized smooth in­version tends to underestimate the true physical parameters of the target (see. for example. Zhdanov. 2002. p. 46-49). To recover the correct conductivity, one should use the nonlinear inversion with the focusing stabilizer, which will be outlined in the next numerical ex­ample of the MT data inversion. This is. however, a very challenging problem in the case of HEM data collected with the moving trans- mitter-receiver pairs because any new position of the transmitter re­quires solving a different forward-modeling problem. The full 3D nonlinear inversion for the multitransmitter airborne data is still im­practical because of the huge computational time required in this case. At the same time, a linearized approach represents an effective solution of 3D inverse problem for the multitransmitter EM data, which can be widely used in practical interpretation of HEM data. We should also note that the limitations of the fast-forward mapping operators in inverse-problem solution are addressed in many publi­cations. including Zhdanov and Tartaras (2002) and Zhdanov and Chernyavskiy (2004). The errors in these approximations affect the uncertainty analysis in the same way that they affect the inversion re­sult itself. In this situation, it is especially important to evaluate how robust the linearized inversion is with respect to the noise in the data. The new method of model appraisal and resolution analysis provides the corresponding mathematical technique for solving this problem. Using the general resolution theory outlined above, we can find the resolution density 7?.f,T and the upper bounds of the conductivity variations, according to nr Aa ‘ Figure 2. The vertical cross sections (along the .v-axis) of the 3D model obtained as a result of regularized inversion for the synthetic observed HEM data for a conductive cubic model. The black out­lines show the true contour of the conductive body. The vertical cross sections of the upper bounds of the conductivity variations SAcr-, for the cubic model inversion are shown in Figure 3. while Figure 4 presents the 1Z.?'T distribution. We have assumed in these calculations that the relative error in the observed data is e = Wc. We can see in Figure 3 that the estimated variations in the HEM data inversion in the area of the cubic body are within 0.00001 S/m. while the estimated inverted conductivity of the body is about 0.001 S/m. Thus, the variations do not exceed Wc. which corresponds well to Wc level of noise in the data. This fact demon- Figure 3. The vertical cross sections of the upper bounds of the con­ductivity variations M^l r) computed for the cubic model HEM data inversion under the assumption that the relative error in the ob­served data is e = Wc. The white outlines show the true contour of the conductive body. Figure 4. The vertical cross sections of the resolution density 7Zf'T distribution for the cubic model HEM data inversion. The white out­lines show the true contour of the conductive body. Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ R84 Zhdanov and Tolstaya strates that our inversion algorithm does not amplify the noise in the data Based on this analysis, we can conclude that the LQL inversion outlined above is a very robust method. It provides a stable geomet­ric image of the target while underestimating the true conductivity. We should note, however, that the detailed analysis of Figure 3 shows that this robust solution is obtained only in the upper and cen­tral parts of the area of inversion. The estimated variations rapidly increase with the depth and to the sides of the inversion area, reach­ing almost 0,0001 S/m, The resolution density decreases for the bot­tom and edge parts of the inverse model, correspondingly, as shown in Figure 4, This simple numerical example shows that, in the practical inver­sion of geophysical data, it is notenough toplotjustthe inverse mod­el obtained by inversion. In principle, the data acquisition and inver­sion schemes, play rather important roles in the way the noise propa­gates into the inverse model. The newly developed method of model appraisal and resolution analysis makes itpossible to locate the parts of inverse image that experience the minimal effect of the noise in the data and the areas with the strongest distortions. As a result, the interpreter can identify the parameters of the inverse model that can be treated with the most confidence, as well as less reliable features. This is the major practical significance of our method of inverse im­age appraisal, RESOLUTION STUDY OF THE NONLINEAR 3D MAGNETOTELLURIC INVERSION ALGORITHM In this section, we demonstrate the application of the developed method of the resolution analysis to a nonlinear, 3D MT inverse problem. The foundations of the MT method were developed by Tikhonov (1950) andCagniard (1953). It is based on measurements of the nat­ural EM field at the surface of the earth. The interpretation of MT data is based on the calculation of the transfer functions between the horizontal components of the electric and magnetic fields, which form the so-called impedance tensor Z (Berdichevsky and Dmitriev, 2002): The components of the impedance tensor are determined from the horizontal components of the electric and magnetic fields at each ob­servation point. The corresponding technique for solving this prob­lem is outlined in Zhdanov and Keller (1994) and Berdichevsky and Dmitriev (2002), These data are inverted fora conductivity model of the earth. Thus, the MT inversion requires forward modeling of EM field components, the corresponding impedances, and the apparent resis­tivities and phases on each iteration step. This procedure is extreme­ly time consuming, resulting in enormous calculations to solve the inverse problem. To overcome this computational difficulty, Zh­danov and Golubev (2003) suggest using an approximate solution based on quasi-analytic (QA) approximation on the initial stage of the iterative inversion. The detailed description of the basic princi­ples of the QA approximation can be found in Zhdanov et al, (2000) and Zhdanov (2002), The approximate QA forward operators, intro­duced in the cited papers, can be used to compute the components of the impedance tensor Z. These operators significantly speed up the computations at each step of the inversion. In a general case, the corresponding equations of MT inversion can be expressed by an operator equation including the data vector d and the vector of model parameters m as d = B(m), (21) where B is the nonlinear forward operator representing the govern­ing equations of the MT impedance modeling problem, m is the vec­tor of the unknown conductivity distribution (model parameters), and d is the vector formed by the observed values of the components of the MT impedance tensor at the observation points. Inversion aims at estimating the model parameter vector m based on B and a known (observed) data vector d. This problem is usually ill posed, i.e., the solution can be nonunique and unstable. The con­ventional way of solving ill-posed inverse problems, according to regularization theory (Tikhonov and Arsenin, 1977; Zhdanov, 2002), is based on minimization of the Tikhonov parametric func­tional, similar to one shown in equation 4, To generate a focused image of the geoelectrical model, Zhdanov and Hursan (2000) and Mehanee and Zhdanov (2002) applied a min­imum support stabilizer, which is a nonquadratic functional of a form %s(m) = (m - mupr)7t(m - m.p,.)2 + f2i]_1(m - mupr), (22) where m and m.,F are N„, X N„, diagonal matrices of inverse-model parameters (current and a priori, respectively) m = diag(m1,m2> •••> mN)> Aupr = diag(m 1 upp MiapH • • • ■ ^Nmapr) ■ where e is the focusing parameter and I is an N„, X N„, identity ma­trix, Portniaguine and Zhdanov (1999) show that this functional minimizes an area of nonzero parameter distribution (minimizes the support of the inverse model) ife tends to zero: e^0. The principles of the minimum support inversion are discussed in detail in Zhdanov (2002). Recently, Zhdanov and Tolstaya (2004) suggest using a nonlinear parameterization to transform the nonquadratic, minimum-support stabilizing functional described by equation 22 into a quadratic one, described by m = [(m - m.p,.)2 + ?2I]-1/2(m - mupr), (23) and m - mupr = f*[i - m2]''"1/2m, (24) where m = {/«,■}. i = 1......N,„ is the original vector of the model pa­rameters; m = {/m,}, i'=l......N„, is a new vector of the nonlinear parameters; and m is a^„rXA,(I diagonal matrix with the diagonal formed by nonlinear model parameters, m = diagt/n^/nj......fnN ). We solve the minimization problem for the corresponding Tikhonov parametric functional by the regularized conjugate gradi­ent (RCG) method. The details of this algorithm are described in Zh­danov and Tolstaya (2004). The application of the QA approximation to forward modeling and Frechet derivative computations speeds up the calculation dra- Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ Resolution analysis o f geophysical inversion R85 matieally. However, to control the accuracy of the inversion, this method allows application of rigorous forward modeling in the final steps of the inversion procedure. We use an integral equation for­ward- modeling code based on the contraction integral equation method, which improves the convergence rate of the iterative solvers (Hursan and Zhdanov, 2002). Application of a few additional iterations with a rigorous forward-modeling solver improves the resolution of the inverse method and helps to gen­erate a more correct image of the target (Zhdanov an d Tol s tay a, 2004). We now present a numerical example of the MT data inversion and the resolution analysis. Consider a homogeneous half-space with a resis­tivity of 100 ohm-m, containing a conductive dike. The resistivity of the inhomogeneity is 3 ohm-m. The top of the dike is at a depth of 200 m, and its bottom is at a depth of 600 m be­neath the surface. This model is excited by a plane EM wave source. The .v- and y-components of the anomalous magnetic and electric fields for four different frequencies (1, 10, 100, and 1000 Hz) have been simulated at 225 receiver points arranged on a homogeneous grid, using inte­gral equation forward-modeling code INTEM3D (Hursan and Zhdanov, 2002). The coordinates of the receiver grid are .v and v from -700 to 700 ev­ery 100 m. The receiver system is located at the surface of the earth. The EM field components were recalculated into MT impedances, using the standard equations (Berdichevsky and Dmitriev, 2002). The area of inversion is covered by a ho­mogeneous mesh consisting of 16 X 25 X 8 cubic cells surrounding the anomalous structure to be inverted. Each cell has a dimension of 100 m in the a-, and .--directions. We select the focusing parameter as e = 0.016. The details of the tech­nique for the optimum e selection can be found in Zhdanov and Tol­staya (2004). Figure 5a shows the true model. For this model, we run inversion using different inverse methods. On the first stage, we ran 15 itera­tions of the RCG method with minimum norm stabilizer (Figure 5b); after that, we applied 60 iterations with a minimum-support stabiliz­er and QA forward modeling and an additional 20 iterations of the RCG method with minimum-support inversion and rigorous for­ward modeling (Figure 5c). The inversion curves, parametric func­tional /'[a], stabilizer ,y[m], misfit ^[m], and elapsed time versus it­eration number are shown in Figure 5d, as well. Figure 6 presents the vertical cross sections of the true model (panel a), the smooth inversion result with minimum-norm stabilizer (panel b), the intermediate result with minimum-support stabilizer and QAforward modeling (panel c), and the final focusing inversion result (panel d). One can see that the smooth minimum-norm result underestimates the true conductivity, while the focusing inversion reconstructs an image that is very close to the true model and with practically the same resistivity. We have analyzed the resolution of our nonlinear-inversion meth­od for the final model presented in Figure 5c and in Figure 6d. We computed the maximum possible variations in the solution of non­linear- inverse problems for a conductive dike for a given level of rel­ative errors in the observed data equal to 3%. The vertical sections of the upper bounds of the variations of conductivity distributions and the resolution density are shown in Figures 7 and 8. One can see that the resolution is higher in the central parts of the sections, and the Figure 5. (a) The true model of a conductive dike, (b) Inversion results after 15 iterations with the minimum-norm stabilizer, (c) Final resultafter60 iterations with the minimum- support stabilizer and QA forward modeling and 20 additional iterations of the RCG method with minimum-support inversion and rigorous full-forward modeling, (d) Inver­sion curves, parametric functional /'[«], stabilizer ,v[m], misfit and elapsed time versus iteration number. In this figure, we present 3D images of the resistivity distribution with volume rendering. The cutoff level of the resistivity for these images is shown in the corresponding panels. For example, the cutoff level p< 5 ohm-m means that only the cells with a value of resistivity less than 5 ohm-m are displayed (after Zhdanov and Tol­staya, 2004). Figure 6. Dike model. Vertical cross sections of the true model (a), the inversion result with the minimum-norm stabilizer (b), the inter­mediate result with the minimum-support stabilizer and QAforward modeling (c), and the final sharp inversion result (d) (after Zhdanov and Tolstaya, 2004). Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ R86 Zhdanov and Tolstaya variations increase with the depth and to the sides of the inversion area. The maximum possible variations in the area of the dike loca­tion do not exceed 0.03 S/m, which constitutes less than 10% of the anomalous conductivity. This example shows that by applying nonlinear regularized- focusing inversion, we obtain a very accurate image of the target. In this situation, the parameters of the true model are located within the upper bounds of the model parameter variations, provided by the resolution analysis. Figure 7. Model of the conductive dike. The vertical sections of the upper bounds of the conductivity variations &lcr,(r) computed for the dike-model MT data inversion under the assumption that the rel­ative error in the observed data is e = 3%. The white lines show the true position of the conductive dike. CASE HISTORY: RESOLUTION STUDY OF 3D INVERSION OF HEM DATA COLLECTED IN THE VOISEY'S BAY AREA We applied our method for the resolution analysis of the 3D inver­sion result to the real HEM data collected by INCO Exploration in the Voisey's Bay area in Canada to examine the resolution of the cor­responding inverse model. This area is characterized by high-con- ductivity nickel-copper sulfide deposits hosted by resistive troctolite dikes (Naldrett et al., 1996). A geologic map of the area with several identified deposits is shown in Figure 9. Zhdanov and Tartaras (2002) used localized quasilinear inversion for the 3D interpretation of the HEM data collected within an area outlined in Figure 9. This area corresponds to the location of the Ovoid deposit, which is a flat- lying deposit of very high conductance comprised of 70% massive sulfide (Balchetal., 1998). Here, we use the same data to demonstrate a new resolution-anal- ysis technique. In the first stage, following Zhdanov and Chern­yavskiy (2004), we apply an SLDM method for 3D inversion of the HEM data. Based on drilling information incorporated in Figure 9, we assume a 20-m deep, conductive overburden with resistivity of 10 ohm-m. We used the coplanar and coaxial components from the lowest frequency (900 Hz) because they are the least sensitive to the presence of the conductive overburden. The data were transformed from parts per million to anomalous field values in teslas, assuming a uniform background resistivity of 1900 ohm-m. The data comprise part of four flight lines at a distance of 200 m from each other (lines A, B, C, and D in Figure 9). The area of inversion was 700 m X 600 m X 160 m and was divided into 14 X 30 X 8 cells. Figure 10 presents the inversion result in the form of vertical slic­es through the model generated as a result of the inversion. Note that the vertical section at.v = 300 m corresponds to flight line B, while the vertical section at.v = 500 m corresponds to flight line C. Figure Figure 8. The vertical sections of the resolution density /$/r distribu­tion, computed for the model of the conductive dike. The white lines show the true position of the conductive dike. Figure 9. A map of the Voisey's Bay sulfide deposits, (a) Plan map. (b) North-facing longitudinal section. The data measured along four flight lines (A, B, C, and D) flown over the Ovoid deposit were used for inversion (afterZhdanov and Tartaras, 2002). Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ Resolution analysis o f geophysical inversion R87 11 shows the observed and predicted data along flight line. A, B, C, and D, shown in Figure 9. The observed data are shown by solid lines. The predicted data obtained by full SLDM comprised of 200 Lanczos steps are shown by the dashed lines, while the dotted curves correspond to the truncated SLDM with only 25 Lanczos steps (Zh­danov and Chernyavskiy, 2004). One can see that the agreement be­tween the three curves is very good. The results seem reasonable and in good agreement with the existing information about the Ovoid de­posit (Balch et al., 1998) and with the inversion result obtained by Zhdanov and Tartaras (2002) using the LQL method with the conju­gate gradient minimization. We have evaluated the resolution of this inversion using a tech­nique described in this paper. Figure 12 presents the vertical cross sections of the estimated variations of the conductivity Mcr/r) for the Voisey's Bay HEM data inversion. We can see in this figure that the estimated variations in the HEM data inversion in the area of the conductive target are within 0.03 S/m, while the maximum inverted conductivity of the Ovoid is about 1.5 S/m. We have assumed in these calculations that the relative error in the observed data is e = 1 %. Thus, variations in the inversion do not exceed 2% within the volume of the conductive target. If we assume that the relative errors of the data could be as high as 10%, the maximum variations in the inverse model will be increased up to 30%, correspondingly. Figure 13 presents the resolution density Ri"T distribution for the same in­verse model, shown in Figure 10. One can see that the resolution density decreases at the bottom and at the sides of the area of inver­sion. We observe the maximum resolution in the area of the conduc­tive Ovoid deposit. The resolution density distribution and the plots of the upper bounds of the conductivity variations, SAa-f, r), for the Voisey's Bay HEM data inversion provide a clear indication that the inversion result would not be distorted significantly by the noise in the data in the central area of the inversion domain, while the esti­mated variations increase with the depth and to the sides of the area of inversion. We should note, however, that the true conductivity of the ore­body in the area of investigation is much higher than 1.5 S/m ob­tained by the inversion. Therefore, in this case, the LQL inversion, as usual, underestimates the true conductivity, which is typical for a lin­earized smooth inversion, as in the case of our model study. At the same time, in airborne exploration for a mining target, similar to the HEM survey conducted in the Voisey's Bay area, it is very well known that the target has high conductivity. The problem is not to find the true conductivity of the mineralization zone but to correctly locate this target. In this situation, the LQL inversion provides a rea- Figure 11. The magnetic field values along four south-north lines A, B, C, and D above the Ovoid. The observed field is shown along with the predicted fields for the inverse model obtained by the full and truncated SLDM method. The labels CXQ and CXI refer to the quadrature and in-phase components of the HEM data, respectively. The units for the vertical axis are teslas (T). Figure 10. Vertical slices through the geoelectrical model generated as the result of the inversion for the Voisey's Bay HEM data. The ver­tical section at x = 300 m corresponds to flight line B; the section at x = 500 m corresponds to line C (Figure 9). The white outlines show the boundary of the conductive body with the conductivity more than 1 S/m, based on the inversion result. Figure 12. Vertical cross sections of the upper bounds of the conduc­tivity variations Mcr/r) for the Voisey's Bay HEM data inversion. The white outlines show the boundary of the conductive body with the conductivity more than 1 S/m, estimated based on the inversion result. Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ R88 Zhdanov and Tolstaya Figure 13. Vertical cross sections of the resolution density Wj^distri- bution for the Voisey's Bay HEM data inversion. The white outlines show the boundary of the conductive body with the conductivity more than 1 S/m. estimated based on the inversion result. sonable estimate of the location of the conductivity anomaly, while the resolution analysis confirms that this location would not change dramatically even if we consider the possible errors in the data and the inherent instability of the inverse problem. Therefore, we may conclude that we have a robust estimation of the location of the con­ductive body, even though the true conductivity is actually underes­timated. CONCLUSIONS The existing techniques for appraisal of geophysical inverse im­ages are based on calculating the data and model resolution and the covariance matrices. This technique provides, for example, a stan­dard deviation of the model parameters from the inversion result. In many practical applications, however, it may be useful to also esti­mate the maximum possible variations in the solution of the inverse problem for the given errors in the observed data. These upper bounds of the model variations determine the actual resolution of geophysical inversion. We have developed a novel approach to the resolution analysis of the regularized geophysical inversion, based on evaluating the upper bounds of the model variations using the Cauchy inequality for the regularized least-squares solution of the linear inverse problem. The inverse of the upper bounds of the model parameter variations determines a new characteristic of geophysical inversion - a resolution density. We have also developed an effi­cient numerical method for the resolution density calculation based on the SLDM. This new technique provides a quantitative evalua­tion of the stability of the inverse problem solution by generating a spatial distribution of the upper bounds of the variations in the in­verse model. At the same time, our method should be treated as a complement to the traditional methods of model appraisal and reso­lution analysis, which may provide additional insight in the inver­sion result. This approach was extended to nonlinear inverse problems as well. Our extension is based on studying the sensitivity and resolu­tion of the nonlinear inverse problem in the vicinity of the final solu­tion. In this case, we approximate the relationship between the errors in the observed data and the corresponding variations in inversion results by the linear Frcchet derivative operator computed for the fi­nal solution. This approximation allows us to conduct the nonlinear resolution analysis based on the same principles that have been de­veloped for the linear inverse problem. The method is illustrated by the resolution analysis of linearized 3D HEM data inversion and nonlinear 3D MT inversion. In the case of linearized smoothed LQL inversion, we underestimate the true conductivity of the target while providing a stable estimate of its lo­cation. The full nonlinear 3D MT inversion with focusing regular­ization makes it possible to accurately recover both the conductivity and the shape and the location of the target. We should note, howev­er. that the upper bounds of the model parameter variations may not represent the true errors of in version (the errors which reflect the dif­ference between the inverse model and the true model). The numeri­cal examples show that the interpretation of these bounds depends on the method used for the inversion. If we apply the nonlinear regu­larized focusing inversion based on rigorous forward modeling, the upper bounds provide the correct errors of the inversion result. In the case of smooth inversion based on the approximate forward solu­tions (like the LQL inversion), the resolution density analysis pro­vides useful information about the pattern of the noise transforma­tion into the inverse image. The important conclusion is that, using the new technique de­scribed in this paper, one can present the regularized inversion result in the form of both the model parameter and the estimated model variation distributions. Thus, the potential end user of the geophysi­cal survey may analyze not only the inverse model produced by a specific inversion algorithm but also a spatial distribution of the up­per bounds of the model parameter variations. This information may provide useful insight into the robustness of the inverse geophysical imaging. The resolution study in the case history of the 3D HEM survey, shows that this new technique provides a useful tool for the analysis of the robustness of a geophysical inversion, which may comple­ment the traditional methods of image appraisal. Although in this pa­per we present application examples in EM geophysical methods, it is possible to extend this theory to other geophysical methods as well. ACKNOWLEDGMENTS The authors acknowledge the support of the University of Utah Consortium for Electromagnetic Modeling and Inversion (CEMI). which includes BAE Systems. Baker Atlas Logging Services. BGP China National Petroleum Corporation. BHP Billiton World Explo­ration Inc.. Centre for Integrated Petroleum Research. EMGS. ENI. ExxonMobil Upstream Research Company. INCO Exploration. In­formation Systems Laboratories. MTEM. Newmont Mining Co.. Norsk Hydro. OHM. Petrobras. Rio Tinto-Kennecott. Rocksource. Schlumberger. Shell International Exploration and Production Inc.. Statoil. Sumitomo Metal Mining Co.. and Zonge Engineering and Research Organization. We also thank INCO Exploration, and particularly A. King, for providing us with the HEM data set. The authors are grateful to the anonymous reviewers for their comments and recommendations, which improved the paper. Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ Resolution analysis o f geophysical inversion R89 APPENDIX A APPLICATION OF THE SPECTRAL LANCZOS DECOMPOSITION METHOD (SLDM) FOR RESOLUTION DENSITY CALCULATION In this appendix, we construct a numerical algorithm for resolu- tion-density calculation. According to the definition, the resolution density is obtained by adding the squares of the column elements of the regularized inverse matrix R„. To find the ith column of this matrix, we can introduce a vector e, with unity in the ith position: [0,0, ..., 0,1,0, ..., 0]7 EnJ, where ENj is AT,,-dimensional Euclidean space, and recover the ith column of the regularized inverse matrix R'Q'* by simple multiplica­tion: R"> = (A'W^A + aW2 )-1A*'W;/ e, = W;(ASW‘A^ + aI)-‘A*W^ e, Introducing the notations (A-I) B = we finally obtain where A'WJAWj A'W^ e„ R (0 (A-2) /„(B) = (B + al)"1. Thus, we have arrived at the problem of computing a function of matrix B. This problem can be solved by the SLDM (Druskin and Knizhner- man, 1994; Golub and Van Loan, 1996; Zhdanov, 2002). First, we apply the Lanczos algorithm for QT decomposition of matrix B: A) = 1, q0 = 0, qi = c,-, while (3j =£ 0, q/V1 = - af Pi r/ = (B - IN)q,- - iq Pi (A-3a) (A-3b) (A-3c) j = I, 2, N - I. As a result, we find an orthogonal matrix Q, and the tridiagonal ma­trix T,, where L is an iteration step of the Lanczos algorithm. We can write expression A-2 as (i) R where i JD (A-4) JU [1,0, ..., 0,0, ..., o]7 el. The resolution density is computed now, according to equation 9, as Nj 2 k,7i L./=1 -1/2 (A-5) The advantage of the SLDM method is that we have to run the Lanczos algorithm only once for all different values of the regular­ization parameter a. After that we only have to invert a tridiagonal matrix (T, + al) for a different a, which is a much simpler opera­tion. The selection of the optimal regularization parameter a can be made using either Tikhonov's method or the /.-curve method. The detailed description of these methods can be found in Zhdanov (2002). REFERENCES Alumbaugh. D. L„ and G. N. Newman, 2(X)(), Image appraisal for 2D and 3D electromagnetic inversion: Geophysics, 65, 1455-1467. Backus, G. E., 1970a, Inference from inadequate and inaccurate data. I: Pro­ceedings of the National Academy of Sciences, 65, 1-7. -------, 1970b, Inference from inadequate and inaccurate data, II: Proceed­ings of the National Academy of Sciences, 65,281-287. 1970c, Inference from inadequate and inaccurate data. III: Proceed­ings of the National Academy of Sciences, 67,282-289. Backus, G. E., and T. I. Gilbert, 1967, Numerical applications of a formalism for geophysical inverse problems: Geophysical Journal of the Royal As­tronomical Society, 13,247-276. -------, 1968, The resolving power ofgross earth data: Geophysical Journal of the Royal Astronomical Society, 16, 169-205. Batch. S. J., T. J. Crebs, A. King, and M. Verbiski, 1998, Geophysics of the Voisev's Bav Ni-Cu-Co deposits: 68th Annual International Meeting, SEG, Expanded Abstracts, 784-787. Berdiche vsky, M. N„ and V. I. Dmitriev, 2002, Magnetotellurics in the con­text of the theory of ill-posed problems: SEG. Cagniard, L„ 1953, Theory of magnetotelluric method of geophysical pros­pecting: Geophysics, 18,605-635. Dmitriev, V. I., M. S. Zhdanov, V. A. Morozov, A. A. Nikitin, and H. P. Brus- netso v, eds., 1990, Computational mathematics and techniques in explora­tion geophysics (in Russian): Nedra Press. Druskin, V., and L. Knizhnerman, 1994, Spectral approach to solving three­dimensional Maxwell's diffusion equations in the time and frequency do­mains: Radio Science, 29,937-953. Druskin, V., L. Knizhnerman. and P. Lee. 1999. New spectral Lanczos de­composition method for induction modeling in arbitrary 3D geometry: Geophysics. 64.701-706. Foster. M„ 1961. An application of the Wiener-Kolmogorov smoothing the­ory to matrix inversion: Journal of the Society for Industrial and Applied Mathematics. 9.387-392. Franklin. J. N„ 1970. Well-posed stochastic extensions of ill-posed linear problems: Journal of Mathematical Analvsis and Applications. 31. 682­716. Golub. G. H., and C. F. Van Loan. 1996. Matrix computations. 3rd ed.: Johns Hopkins Univ. Press. Hursan. G., and M. S. Zhdanov. 2002. Contraction integral equation method in 3D electromagnetic modeling: Radio Science. 37. 1089 doi:10.1029/ 2001RS002513. Jackson. D. D„ 1972. Interpretation of inaccurate, insufficient and inconsis­tent data: Geophysical Journal of the Roval Astronomical Societv. 28.97­110. Lanczos. C., 1961. Linear differential operators: D. van Nostrand Co. Marquardt, D. W., 1963. An algorithm for least-squares estimation of nonlin­ear parameters: Journal of the Society for Industrial and Applied Mathe­matics. 11.431-441. -------. 1970. Generalized inverses, ridge regression, biased linear estima­tion. and nonlinear estimation: Technometrics, 12.591-612. McGillivray. P. R., and D. W. Oldenburg. 1990. Methods for calculating Frechet derivatives and sensitivities for the nonlinear inverse problem: A comparative studv: Geophysical Prospecting. 38.499-524. McGillivray. P. R„ D. W. Oldenburg. R. G. Eliis, and T. M. Habashy, 1994. Calculation of sensitivities for the frequency-domain electromagnetic problem: GeophysicalJournallnternational. 116. 1-4. Mehanee. S.. and M. S. Zhdanov. 2002, Two-dimensional magnetotelluric inversion of blockv geoelectrical structures: Journal of Geophysical Re­search. 107. B4.2065. doi: 10.1029/2001JB000191. Menke. W„ 1989. Geophysical data analysis: Discrete inverse theory: Aca­demic Press Inc. Naldrett. A. J.. H. Keats. K. Sparkes, and S. Moore. 1996. Geology of the Voi­sev's Bav Ni-Cu-Co deposit. Labrador. Canada: Exploration and Mining Geology' 5.169-179. Oldenburg. D. W., 1983. Funnel functions in linear and nonlinear appraisal: Jounuilof Geophysical Research. 88. B97387-7398. Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ R90 Zhdanov and Tolstaya Parker, R. I.., 1975, The theory of ideal bodies for gravity interpretation: Geophysical Journal of the Royal Astronomical Society, 42,315-334. Portniaguine, O., and M. S. Zhdanov, 1999, Focusing geophysical inversion images: Geophysics, 64,874-887. Ramirez, A. L„ W. D. Daily, and R.L. Newmark, 1995, Llectrical resistance tomography for steam injection monitoring and control: Journal of Bnvi- ronmentaland Bngineering Geophysics, 0,39-52. Rodi, W. L., 1976, A technique for improving the accuracy of tin ite element solutions for magnetotelluric data: Geophysical Journal of the Royal As­tronomical Society, 44,483-506. Sabatier, P. C., 1977a, On geophysical inverse problems and constraints: Journal of Geophysical Research, 43,115-137. ----, 1977b, Positivity constraints in linear inverse problems - I. General theory: Geophysical Journal of the Royal Astronomical Society, 48, 415-441. ----, 1977c, Positivity constraints in linear inverse problems - II. Appli­cations: Geophysical Journal of the Royal Astronomical Society, 48,443­459. Sambridge, M., and K. Mosegaard, 2002, Monte Carlo methods in geophysi­cal inverse problems: Reviews of Geophysics, 40,3,1-29. Spies, B. R., and T. M. Habashy, 1995, Sensitivity analysis of cross-well electromagnetics: Geophysics, 60,834-845. Stark, P. B., and R. L. Parker, 1987, Velocity bounds from statistical esti­mates of Hp) and X(/>): Journal of Geophysical Research, 92, B3, 2713-2719. Stark, P. B., R.L. Parker, G. Masters, and L. A. Orcutt, 1986, Strict bounds on seismic velocity in the spherical earth: Journal of Geophysical Research, 91,B14,13892-13902. Tarantola, A., 1987, Inverse problem theory: Blsevier Science Publ. Co. Inc.. Tarantola, A., and B. Valette, 1982, Generalized nonlinear inverse problem solved using the least squares criterion: Reviews of Geophysics and Space Physics, 20,219-232. Tikhonov, A. N., 1950, On determining electric characteristics of the deep layers ofthe earth's crust: DokladyAkademii NaukSSSR, 73,295-297. Tikhonov, A. N., and V. Y. Arsenin, 1977, Solution of ill-posed problems: W. H. Winston and Sons. Xiong, Z., 1992, BM modeling of three-dimensional structures by the meth­od of system iteration using integral equations: Geophysics, 57, 1556- 1561. Zhdanov, M. S., 2002, Geophysical inverse theory and regularization prob­lems: Blsevier Science Publ. Co., Inc. Zhdanov, M. S., and A.Chernyavskiy, 2004, Rapid three-dimensional inver­sion of multi-transmitter electromagnetic data using the spectral Lanczos decomposition method: Inverse Problems, 20, S233-S256. Zhdanov,M. S., V. I. Dmitriev, S.Fang, and G. Hursan, 2000, Quasi-analyti- cal approximations and series in electromagnetic modeling: Geophysics, 65,1746-1757. Zhdanov, M. S., and N. G. Golubev, 2003, Three-dimensional inversion of magnetotelluric data in complex geological structures, in J. Macnae and G. Liu, eds., Three-dimensional electromagnetics III, Australian Society of Bxploration Geophysics. Zhdanov, M. S., and G. Hursan, 2000, 3D electromagnetic inversion based on quasi-analytical approximation: Inverse Problems, 16,1297-1322. Zhdanov, M. S., and G. W. Keller, 1994, The geoelectrical methods in geo­physical exploration: Blsevier Science Publ. Co., Inc.. Zhdanov, M. S., and Ii. Tartaras, 2002, Inversion of multi-transmitter 3D electromagnetic data based on the localized quasi-linear approximation: Geophysical Journal International, 148,506-519. Zhdanov, M. S., and Ii. Tolstaya, 2004, Minimum support nonlinear parame­terization in the solution of 3D magnetotelluric inverse problem: Inverse Problems, 20,937-952. Downloaded 17 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/