Underground imaging by frequency-domain electromagnetic migration

A new method of the resistivity imaging based on frequency-domain electromagnetic migration is developed. Electromagnetic (EM) migration involves downward diffusion of observed EM fields whose time flow has been reversed. Unlike downward analytical continuation, migration is a stable procedure that accurately restores the phase of the upgoing field inside the Earth. This method is indented for the processing and interpretation of EM data collected for both TE and TM modes of plane-wave excitation. Until recently, the method could be applied only for determining the position of anomalous structures and for finding interfaces between layers of different conductivity. There were no well developed approaches to the resistivity imaging, which is the key problem in the inversion of EM data. We provide a novel approach to determining not only the position of anomalous structures but their resistivity as well. The main difficulty in the practical realization of this approach is determining the background resistivity distribution for migration. We discuss the method of the solution of this problem based on differential transformation of apparent resistivity curves. The final goal of migration is to provide a first order interpretation using a computational effort equivalent to a forward modeling calculation.

GEOPHYSICS, VOL. 61, NO. 3 (MAY-JUNE 1996); P. 666-682, 7 FIGS. Underground imaging by frequency-domain electromagnetic migration Michael S. Zhdanov*, Peter Traynin^, and John R. Booker** ABSTRACT A new method of the resistivity imaging based on frequency-domain electromagnetic migration is devel­oped. Electromagnetic (EM) migration involves down­ward diffusion of observed EM fields whose time flow has been reversed. Unlike downward analytical contin­uation, migration is a stable procedure that accurately restores the phase of the upgoing field inside the Earth. This method is indented for the processing and inter­pretation of EM data collected for both TE and TM modes of plane-wave excitation. Until recently, the method could be applied only for determining the position of anomalous structures and for finding inter­faces between layers of different conductivity. There were no well developed approaches to the resistivity imaging, which is the key problem in the inversion of EM data. We provide a novel approach to determining not only the position of anomalous structures but their resistivity as well. The main difficulty in the practical realization of this approach is determining the back­ground resistivity distribution for migration. We discuss the method of the solution of this problem based on differential transformation of apparent resistivity curves. The final goal of migration is to provide a first order interpretation using a computational effort equivalent to a forward modeling calculation. INTRODUCTION A problem of current practical interest is that of imaging inhomogeneous underground structures using surface or bore­hole electromagnetic data. The last decade has seen consider­able improvement in the ability to gather spatially dense, accurate EM induction data. Improved extraction of structural information from such data is important for many practical applications ranging from mineral exploration to waste and building site characterization. Recently there has been consid­erable progress in the direct inversion of magnetotelluric and other low-frequency induction data for multidimensional elec­trical conductivity structures (Berdichevsky and Zhdanov, 1984; Smith and Booker, 19881991; Wannamaker et al., 1989; Eaton, 1989; Madden and Mackie, 1989; deGroot-Hedlin and Constable, 1990; Xiong and Kirsh, 1992; Lee and Xie, 1993; Oristaglio et al., 1993; Pellerin et al., 1993; Torres-Verdin and Habashy, 1994; Tripp and Hohmann, 1993; Oldenburg et al., 1993; 1994; Zhdanov and Keller, 1994). However, these algo­rithms require repeated solution of large multidimensional forward problems. Several publications address simple and fast inversion tech­niques for transient electromagnetic data over inhomogeneous structures (Barnett, 1984; Macnae and Lamontagne, 1987; Eaton and Hohmann, 1989). The majority of these papers have been based on equating the transient response, measured at the surface of the Earth to the EM field of current filament images of the source. This approach originated in the pioneer­ing work of Nabighian (1979) describing the behavior of transient currents diffusing into the earth as a system of "smoke rings" blown by the transmitting loop into the earth. We will outline a different approach based on direct trans­formation of the observed wavefield into a resistivity image of the cross-section in a very rapid manner. We call this approach electromagnetic migration. The goal of migration is to provide a first-order interpretation using a computational effort equiva­lent to a forward model calculation. The basic ideas of EM migration were first formulated in the papers by Zhdanov and Frenkel(1983a, b), where the integral approach to the solution of the problem of downward extrap­olation in the reverse time has been exposed. The method described in those publications had two main limitations: first, the method is based on Stratton-Chu type integrals and therefore requires (in general) the observation of all six components of the EM field, which is difficult to realize in practice; second, it can be applied only to the models with the homogeneous background resistivity. The more advanced Manuscript received by the Editor April 5, 1994; revised manuscript received September 5, 1995. *University of Utah, Department of Geology and Geophysics, WBB 717, Salt Lake City, UT 84112. ^University of Utah, Department of Geology and Geophysics, 215 Mines, Salt Lake City, UT 84112. **University of Washington, Geophysics Program AK-50, Seattle, WA 98195. © 1996 Society of Exploration Geophysicists. All rights reserved. 666 Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ F re q u e n c y -d om a in EM M ig ra t io n 667 analysis of this approach has been presented in part in the more recent publications: Zhdanov, 1988; Zhdanov et al., 1988 (in Russian); Zhdanov and Booker, 1993; and Zhdanov and Keller, 1994. However, during recent years the method has been significantly developed and improved. Now we have a much clearer understanding of the physical principles and mathematical foundations of EM migration. Until recently, the method could be applied only for determining the position of anomalous structures and for finding interfaces between layers of different conductivity. There were no well developed ap­proaches to imaging the resistivity property itself, which is the key problem in the inversion of EM data. We present a new method of the resistivity imaging based on frequency-domain electromagnetic migration. We develop this method for the processing and interpretation of EM data collected for both transverse electric (TE) and transverse magnetic (TM) modes of plane-wave excitation. To make the presentation clearer, we feel it necessary to begin our paper with a short review of the concepts of the migration as applied to interpretation of EM data. So, for completeness, in the first sections of this paper we will outline the physical principles and mathematical foundations of EM migration. We illustrate the method with numerical examples, typical for geoelectrical exploration. PHYSICAL PRINCIPLES OF EM MIGRATION The physical principles of EM migration parallel those underlying laser holography and seismic migration. The EM field produced by a controlled or natural source and observed on the earth's surface is the combination of two fields: a primary field that propagates downwards into the earth, and a secondary field that propagates upwards after having been scattered back from internal structure. Both fields satisfy diffusion equations inside the earth. Their amplitudes decay and their phases are retarded in the direction of propagation that is downward for the source fields and upward for the back-scattered fields. Given measurements of the electric and magnetic fields at the Earth's surface it is possible to separate the downgoing and upgoing fields (Berdichevsky and Zhdanov, 1984). The measured surface expression of the upgoing field can then be used to reconstruct information about the field inside the earth and, hence, estimate the geoelectric structure (Zhdanov and Frenkel, 1983a, b; Lee et al., 1987; Zhdanov, 1988; Zhdanov et al., 1988). As in seismic migration, one could extrapolate the signal received back to reflectors or internal currents by using them as a boundary condition for a diffusion equation; this transformation is usually called an analytic continuation. In the presence of measurement noise, however, this is an unstable process. An alternate procedure is to reverse the time flow of the back-scattered signals received at each site and then diffuse these time-reversed signals downward using the ordinary diffusion equation. The diffused time-reversed fields are called migrated fields. Their amplitude will decay downward and will be very different from the original upgoing field whose amplitude increases downward. The usefulness of the migrated field arises because of the following facts: 1) The phase delay associated with diffusing the time- reversed signals corresponds to a downward phase ad­vance in ordinary time. Thus the downward phase behav­ior of the migrated fields is essentially the same as that of the original upgoing field and the migrated signal summed for sources at all the sites should exhibit the same constructive interference at internal scatterers as the original field. 2) The noise in the migrated field decays downwards be­cause the migrated field satisfies an ordinary diffusion equation. This contrasts with the explosion of error associated with downgoing analytic continuation by using the signal received itself as a boundary condition for a time-reversed diffusion equation. Thus the phase of the upgoing field inside the earth can be estimated more accurately from the migrated field. There are a number of algorithms developed for migration of ground-penetrating radar (GPR) reflection profile data (Hogan, 1988; Fisher et al., 1989; Fisher et al., 1992). All these algorithms are based on kinematic similarities between radar and seismic wave propagation and can be considered as the direct implementation of seismic migration techniques to radar data. However, at low frequencies or in conductive environ­ments where conduction currents are big enough, seismic type migration based on the wave equation is no longer appropriate for processing of EM data, because EM fields diffuse into the ground. There are also several publications concerning the possibility of transforming diffusive EM fields to wavefields (Lavrent'ev et al., 1980; Lee et al., 1989). The integral transform relates the diffusive field in time to a unique wavefield in a time-like domain weighted by an exponentially damped kernel. One can interpret these transformed wavefield data using the conven­tional techniques developed for the wavefields (say, using ray tomography: Lee and Xie, 1993). The main limitation of this approach is that the integral transform is ill-posed, so the explosion of noise could destroy the results of the transforma­tion if one doesn't apply a regularization. Here, we investigate a different approach to EM imaging. Instead of transforming the diffusive field into a wavefield, we transfer the principles of wavefield analysis to interpretation of the EM fields, which are governed by diffusion equations. Thus we develop the method of EM migration. First, we review some general ideas concerning seismic migration, or seismoholography. Suppose that we have a local source of seismic waves, located at some point on the earth's surface, and an array of receivers. Each receiver records oscillations at the earth's surface as a function of real time t. We introduce the reverse time t = (1) Now replace the receivers by auxiliary sources and make each of these sources operate in reverse time with a signal equal to the recording of the earth's surface oscillating in real time at the corresponding receiver. It is shown in the theory of seismic migration that this field is back propagating, that is, it goes from the observation surface into the earth (Claerbout, 1985). If we recalculate the back propagated, or migrated, field at any interior point of the medium at times corresponding to the arrivals of the direct waves from the actual source, the amplitude distribution of the migrated field will depict the positions of the reflectors and the diffraction points. Thus the restoration of the seismic image of the geological cross-section Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ 668 Z h d a n o v et al. is attained by assigning reverse time pseudo-sources to the receiver sites on the earth's surface. The analogous approach in principle may be applied to the interpretation of EM field data as well. Let us consider the situation where we have measured the total EM field, because of natural sources in the ionosphere or a controlled source. The system of synchronized receivers is located at the surface of the earth. We can replace the receivers by a system of artificial current or charge sources, which are determined by the observed EM field. When these artificial sources operate in reverse time, they produce a field that we will call the migrated EM field. As in the seismic case, this field in principle can "delineate" boundaries of the internal structure of the Earth and give us the "geoelectric image" of the earth interior. MATHEMATICAL CONCEPTS Now we will give a stricter definition for EM migration. This definition is much more general than one given in our previous publications. However, in the special case of a homogeneous background cross-section it is reduced to the original defini­tion. Consider a model in which the horizontal plane z = 0 separates the conductive Earth (z > 0) from a nonconducting atmosphere (z < 0). The conductivity of the Earth a(r) is an arbitrary function of the coordinates that can be represented as the sum of normal conductivity cr„ (r) and an anomalous conductivity Acr(r) so that a(r) = u„(r) + Aa(r). The EM field in the model is excited by arbitrary sources, located in the ionosphere or at the surface of the Earth. Field components observed on the Earth's surface are denoted by {E°(r, t), Ey(r, t), 0} and {H°x(r, t), Hy{r, t), H°z(r, t)}. We shall call the migrated field, Em (r, T)Hm (r, t), the field meeting the following conditions: I curl H m(r, t) = an{r)Em(r, t) for z > 0, (4) { Em( ^ _ <9Hm ^ ^ {Hm{r, t), Em(r, t)} 0 for |r| z > 0. (5) Thus we see that the migrated field Em, Hm is the EM field in the reverse time t. It is necessary to use the negative of the vertical component of the observed magnetic field at the right side of equation (3) to make the migrated field satisfy Maxwell's equations up to the surface of observation z = 0, because the observed field E , H satisfies Maxwell's equa­tions in the real time t. In real time t = - t, the migrated field satisfies the adjoint equations I curl Hm = <T„Em dHm (6) curl Em = (jlo • These equations imply that the migrated field is propagating in space from receivers to sources. That is to say, it is a back- propagating field. For the sake of simplicity, consider a situation when the background conductivity of the Earth an is constant. In this case, the electromagnetic field in the model will satisfy the diffusion equation V2H - aH - 0 d V2E - aE - 0 (7) everywhere outside the zones with anomalous conductivity, and we can discuss the problem of migration of any scalar component P(r, t) of the observed EM field. Note first of all that everywhere outside the zones with anomalous conductivity this component would satisfy the equation Let P° (r, t) stand for any of the components H®,Hy, H? , E®, or Ey measured at the Earth's surface. Then we shall call the migrated field P1" of the specified scalar component P° of the EM field, the field satisfying the conditions and V2Pm{r, t) - ix0g 9P T) = 0 for z > 0, ( 1 0 ) OT Notice that if we substitute the ordinary time t into equation (10) we shall have the adjoint diffusion equation, If the ordinary diffusion equation describes the process of field propagation from the sources to the receivers, then equation (12) describes the inverse process of field propagation from the receivers focusing to sources. Thus, the problem of establishing the migrated field reduces to a continuation of the field P0 from the Earth's surface to the lower half-space in the reverse time t. We call the solution of this problem EM migration. As is seen from the exposition above, the calculation of the migrated field is reduced to a boundary value problem de­scribed by the formulas (2)-(5) in general or by the formulas (9)-( 11) in the special case of uniform background. Now we can develop different techniques for solving these problems. In the following sections, we briefly describe solution techniques based on spectral representations of the field in the wavenum- ber-frequency domain and finite-difference approximations. Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ F re q u e n c y -d om a in EM M ig ra t io n 669 ANALYTIC CONTINUATION AND MIGRATION IN THE (K, co)-DOMAIN FOR A MODEL WITH CONSTANT BACKGROUND CONDUCTIVITY A specified component of the EM field P in the form of the Fourier integral with respect to spatial and time frequencies kx, ky, co is X exp [~i(kxx + kyy + coO] dkx dky du>, (13) where P(kx, ky, z, co) is the 3-D Fourier transform of the field component P. Let us rewrite expression (8) in the frequency domain for the Fourier transform P(kx, ky, z, co) d2 ^P(kx, ky, z, co) = v2P(kx, ky, z, co), (14) dz where = (kx + k2 - /co|x0a), Re (r;) > 0 is a wavenumber in the (k, co) domain and 0 <z< d with d being the distance from the earth's surface closest to the surface zone with anomalous conductivity. The general solution of the last equation has the form (15) P(kx, ky, z, co) = Pu(kx, ky, co) X exp (vz) + Pd(kx, ky, co) exp (-vz), where Pu (kx, ky , co) and Pd (kx, ky, co) are the spectrums of the upgoing and downgoing components of the field on the surface of the earth. Equation (15) solves the problem of EM field downward analytic continuation in the homogeneous layer 0 < z < d, if we know the upgoing and downgoing parts of the field. We can use the approach developed in Berdichevsky and Zhdanov (1984) for the field separation into the upgoing and downgoing parts. As an illustration, we present the simple technique for field separation in the 2-D case in Appendix A. Thus the analytic continuation of the downgoing and upgo- ing (or scattered) parts of the field is described by the formula Pcd(kx, ky, z, co) = Pd(kx, ky, 0, co) exp (-vz), (16) u(kx ryCd ky, z, co) = Pu(kx, ky, 0, co) exp (vz), w; - i \n,x, ivy, \j, vjj wwp \u*j, (17) where Pcd and Pcu are results of analytic continuation. Because the exponential in equation (17) is growing with depth, the downward continuation of the upgoing field is an unstable, ill-posed procedure, while the downward continua­tion of the downgoing field is stable. From the other side, the Fourier transform of the migrated upgoing field Pnu (kx,ky, z, co), according to equation (12) satisfies the equation dl dz2 Pmu*(kx, ky, z, co) = v2Pn *(kx,ky,z, co), (18) where v2 = (kx + k2 + ico |x0 a>, asterisk * denotes complex conjugate values, and 0 < z < + 00. Solving the last equation and taking into account condi­tion (11) we arrive at the expression for the complex conjugate spectrum of the migrated upgoing field at a depth z Pmu*(kx, ky, z, co) = Pu(kx, ky, 0, co) exp (-vz), (19) where we choose the branch where Re v = (kx + k2 + /cofiocr)1/2 > 0. Equation (19) gives us the frequency-domain algorithm for migration of the EM field components, which is an EM analog to the migration technique discussed in Gazdag (1978). Obviously, the function f(v, z) = exp (- vz) can be regarded as the frequency response of a low-pass, space-time filter. Therefore, the migration transformation of the EM field is a stable procedure. FINITE-DIFFERENCE MIGRATION IN A 2-D MODEL WITH THE SLOW HORIZONTAL VARIATION OF THE CONDUCTIVITY For the 1-D geoelectrical model, formulas (16), (17) and (19) will be reduced to the following: Pcd(z, co) = Pd(0, co) exp (ikn(z)z), Pcu(z, co) = Pu(0, co) exp (~ikn(z)z), Pmu*(z, co) = Pu(0, co) exp (~kn(z)z), (20) (21) (22) where kn (z) = VaojutQa~(z) is the wavenumber, Re kn > 0; P , Pcu are results of the analytic continuation of the downgoing and upgoing parts of the field, and Pmu is the migrated upgoing field. Let us now consider the 2-D geoelectrical model and analyze the TE mode first. Following Lee et al. (1987) by the analogy with equations (15), (20), and (21) we can represent the y component of electric field approximately by the formula: QdE(x, z, u>)eik"z + Que{x, z, u))e~ik"z, (23) Ey(x, z, co) where the coefficients depend slowly (i.e., continuous within each layer) on depth, and kn (x, z, co) = Vi(n^Q(jn(x, z) is a wavenumber and vn(x, z) is a back­ground conductivity. Here the term associated with the downward-decreasing exponential function corresponds to the downgoing part of the electric field, and the term associated with the downward increasing exponential function corresponds to the upgoing part of the electric field Ed(x, z, co) = QdE(x, z, u>)elknZ, Ey(x, z, co) , -iknz = Qe(x, z, co)e (24) The electric fields Ey,u (x, z, co) everywhere in the lower half-space satisfy the Helmholtz equations dx dz Ed,u(x, z, co) + kl(x, z, u>)Ed,u(x, z, co) = 0. (25) Magnetic fields can be expressed in terms of the electric fields by the second Maxwell equation 1 _ / CO JJL o dEfu dz ' d,u H d,u _ 1 BE /CO(JL o dX (26) Substituting the expressions for the electric fields from equation (24) into equation (25), neglecting the derivatives of kn (i.e., assuming that kn is locally constant), and omitting the exponential term we have Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ 670 Z h d a n o v et al. a2ef d2Qr ■ + --- ± 2ik n dx dz dQdEU dz o, (27) where "+" stands for downgoing field and stands for upgoing field. Differentiating the last expressions by z, we have a 3 n d'u x 2n dv 1 33Q£-h 2ik„ dx2dz d'QT ZzQe (28) Adding equations (27) and (28) with the proper sign and neglecting the third derivative of Qby z gives a3(2f d2Qt dx2dZ = +2 ikr dx 2~+ (2 ikn)2 dz (29) where "+" stands for downgoing field and "-" stands for upgoing field. Solving the last equations numerically and substituting the results in equation (24), we find the upgoing and downgoing components of the field inside the Earth. The corresponding numerical method is discussed in Appendix B. Now let us consider the TM mode. For this model by the analogy with equation (23) we can represent they component of magnetic field approximately by the following formula: , -ik„z Hy(x, z, co) = Qm(x, z, u>)elknZ + Qm(x, z, co)e (30) where Qfyu are the magnetic coefficients that slowly depend on depth. Here the term associated with the downgoing exponential function corresponds to the downgoing part of the field and the term associated with the upgoing exponential function corresponds to the upgoing part of the field Hy(x, z, co) = Qm(x, z, (o)e Hy(x, z, co) = Qm(x, z, co)e -iknz (31) Magnetic fields Hy,u (x, z, co) everywhere in the lower half-space satisfy the equation d I 1 d 'A dx 'kn(x, z, co) dx ) d( 1 aV J + dz \k%(x, z, co) dz) X H"'u(x, z, co) + Hfu{x, z, co) = 0. (32) Electric fields can be expressed in terms of the magnetic fields by the first Maxwell equation ?d,u , 1 dHd* Edu_ dz l dH: Tn dX d,u (33) Substituting the expressions for the magnetic field from equation (31) into equation (32) and repeating the transfor­mations described above, we obtain *3Qit *2Qit „ wit 9 - + 2ik n dx dz dx Y~+ (2 iknY dz (34) where "+" stands for downgoing field and "-" stands for upgoing field. In the case of migration taking into account equation (22), we use the following formulas for the downward extrapolation of the upgoing fields (35) In this case, the migrated upgoing electric field E™u (x, z, co) everywhere in the lower half-space satisfies the following Helmholtz type equation: (36) ET*(x,z, (o) - k2n(x, z, a)E™*(x, z, <o) = 0, (37) while the migrated upgoing magnetic field H™u (x, z, co) satisfi.es the equation X//ymM*(x, z, co) H™u*(x, z, co) = 0. (38) Combining equations (35), (36) and (37), (38) and repeat­ing the transformations described above we obtain dx dz °Qe,m dz ' (39) Equations (29), (34), and (39) contain only first-order ver­tical derivatives of Qe'$T • Therefore, by means of these equations it is possible to extend separately the downgoing and upgoing parts of the total electric or magnetic field known at the level z to the deeper level z + AZ, using the finite- difference approximation to equations (29), (34) and (39). A detailed description of this algorithm is presented in Appendix B. The important conclusion is that in the framework of the model with the slow conductivity variation, we can use the same numerical code for migration of both the electric field in the TE mode and the magnetic field in the TM mode. IMAGING GEOELECTRIC BOUNDARIES BY MIGRATION In this section, we formulate the principles of underground imaging, based on EM migration. Initially, let us consider a two-layered model with the slow variation of conductivity or t (x, z) and cr^ +1 (x, z) within each layer and sharp conduc­tivity contrast on the quasi-horizontal boundary S between two layers. First we analyze the behavior of the horizontal component of the electric field (TE mode) at the quasi-horizontal bound­ary {5: (x, ziS(jc))} between two layers. In the first layer, according to equation (24) we have CEy(x, z, crt) = q£(x, z, a>)elktz+QuE{x, zs, u)e~,ktz {E'y(x, z, o>) = ikeQi(x, z, u)eik,z-ikeQuE(x, z, w)e_ittz, (40) where prime denotes the vertical derivative of the electric field In the second layer we can write ik e+iz \Ey(x, z, co) = Qe(x, Z9 co)e (41) On the boundary S( z = z^ (x)) in the case of E -polariza­tion both components Ey and E'y are continuous. Therefore, the corresponding right-hand sides of equations (40) and (41) are equal. Solving this system of equations, we find Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ F re q u e n c y -d om a in EM M ig ra t io n 671 Que ^ ^ Qe ' 5 ' (42) where PW(x,ZS) = z_s) - ^J+lix, Z_s) yjvdx, Zs) + l(*> zs) (43) is the so-called reflectivity coefficient. Let us calculate the electric apparent reflectivity function as the ratio of the upgoing and downgoing electric fields: (B Ea(x, Z, Co) = E"(x, z, co) Qe(x, z, co) Ey{x, z, co) Qe(x, Z, co) According to equation (42), at the boundary S P Ea(x, ZS, CO) = P(€)(X, ZS). e-2iktz. (44) (45) So, at the geoelectrical boundary, the electric apparent reflec­tivity function is exactly equal to the true reflectivity coeffi­cient! In the case of the TM mode we can repeat all these calculations for the horizontal component of the magnetic field Hy. As the result, we obtain (46) Therefore we can introduce the magnetic apparent reflec­tivity function as the ratio of the upgoing and downgoing magnetic fields (note the minus sign) Z, Q>) = " Hy(x, Z, 0)) Hf(x, z, co) Qm(x, z, co) to) , -2ik ez (47) According to equation (46) at the boundary P Ma(x, ZS, CO) = PW(X, ZS) (48) So, at the geoelectrical boundary magnetic apparent reflectiv­ity function is also equal exactly to the true reflectivity coeffi­cient. Thus we see that although the phases of the downgoing and upgoing electric and magnetic fields cpJ,M (x, z, co) and (x, z, co) in general depend on the frequency co, close to the geoelectric boundaries their difference becomes approxi­mately independent of frequency, according to equations (45) and (48) (because the reflectivity coefficient (3^ is real). From the physical point of view, we have the effect analogous to the seismic case, when the primary (downgoing) and scattered (upgoing) components of the field have the same phase at the position of a reflector. This property of the wavefield is usually used as an imaging condition in a seismic phase migration algorithm (Claerbout, 1985). We can use this imaging condi­tion for electromagnetic fields as well, in which a geoelectric boundary plays the role of a reflector. We can express electric or magnetic apparent reflectivity functions Ma (r, co) as $Ma(x, z, co) = -Hy/Hy = \HuyIHdy\exp (i(<fuM -vdM± it)). (50) The normalized values of (3 E^Ma depend only on the phases differences PL = IP Ea I = exp (/(*£ - <p|)); (510 We have mentioned above that if the point of observation (x, z) approaches the geoelectrical boundary S, then where Acp (x, z) doesn't depend on frequency (it is equal to 0 or tt). Therefore we see that the phase of Ma is significantly frequency dependent except at an interface of high conductiv­ity gradient, where it will be approximately independent of frequency. Thus, stacking pjr^ for a spectrum of frequencies (coi < < co? < < co7) results in positive reinforce- ment at interfaces and destructive interference elsewhere: _ 1 p£,A/a(*> Z) ~ T 2 p£,Afa(*> Z> Mj) J=l and PE,Ma(X> PE,Ma(X' ►0, if(x, z) £ S; >1, if (x, z) G S. (53) (54) The same consideration can be applied not only to a two-layered model, but to a multilayered cross-section as well. Indeed, consider the N-layered model with the slow variation of conductivity (x, z) within each layer ( € = 1, 2, . . . N) and sharp conductivity contrast on the quasi-horizontal bound­aries S€ between the €th and (€ + l)-th layers. For any given boundary S € one can select a frequency oo^ so high that the skin depth S(co^) of the field penetration inside the Earth is less than the depth to the boundary S€ +1: 8(cow) < z s,,, ,if a> co (0 (55) In this case, one has exactly the same expressions for the field components inside the €th and (€ + l)-th layers as equations (40) and (41). Therefore, we can repeat all the mathematical analysis described above and obtain the same results. The only difference is that now the stacking of the normalized apparent reflectivity functions for different depths (z) corresponding to the skin depth §( coy ) has to be done for different frequency intervals: P EMa(X> Z ~ o )) - 1 J-j 0 + (56) where coyo is the lowest frequency of the stacking interval (coy0, co7). The procedure described above is analogous to the conventional vertical sounding of the geoelectrical cross sec­tion: when we migrate the field to shallow depth we use only high frequencies, while migrating deeper we involve more lower frequencies in the calculations, thus "scanning" the vertical cross section. It is important to emphasize that for this kind of imaging it is necessary to reconstruct only the phase of the upgoing field Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ 672 Zhdanov et al. inside the conductive earth. If we compare the phase frequency characteristics of the analytical continuation of the upgoing field [equation (17)] and of its complex conjugate migration [equation (19)], we see that they are equal! That means that the complex conjugate migrated upgoing field has the same phases as the upgoing field itself. Therefore for imaging we can use the migrated reflectivity function (x, z, w), equal to PEm = Euym'IEf„ or p lMm = -H™'IHdy, (57) where E *, Ht‘m * are the complex conjugate migrated upgoing fields. The inphase summation of the migrated reflec­tivity functions indicates the position of the boundaries be­tween layers with different conductivities. Now we can demonstrate the imaging principles for simple synthetic structures, thus providing evidence of the stability of EM migration. For practical calculations we have developed a simplified migration code, which is based on the assumption of a constant background conductivity. Note that in this paper we discuss mainly the results of numerical modeling for plane- wave excitation. In principle, the same approach is applicable for the controlled-source data as well, however, this problem remains to be analyzed and examined. Figure la depicts a mode1 of a 2-D step-wise structure. Figure 2a and Figure 2b show the corresponding apparent resistivity and phase pseudosections, computed for the TE mode. Migration in the frequency domain is realized by a finite-difference method. The migrated field Eum (x, z) corre­sponding to the surface values of the upgoing part of the observed field E"(x, 0) in the model was calculated in the lower half-space, using a finite-difference algorithm. The mi­grated reflectivity function p£m (x, z, ") was then calculated for each position (x, z) and each frequency u>. As we have shown above, stacking for a spectrum of frequencies results in positive reinforcement at interfaces and destructive interference elsewhere. Figure 2c shows a stacked normalized apparent reflectivity function (3Ea(x, z). The maximum in p£a(*, z) shown in Figure 2c almost coincides with the interface between two layers, which clearly demonstrates the phase coherence of the migrated field near the reflector. Thus the migration image produces the correct position of the interface. Figure 2f shows the same model derived from surface data to which 20% Gaussian noise was added as shown on the apparent resistivity and pseudo-sections in Figure 2d and Figure 2e. Nevertheless one can see that the result of the migration produced a rather clear image of the interface. Figure lb depicts a locally conductive rectangular-insert 2-D model. The resistivity of the inclusion is 0.5 ohm^m, and the resistivity of the host rocks is 50 ohm^m. Corresponding apparent resistivity and phase pseudosections are shown in Figure 3a and Figure 3b. The local maximum in (3£a(*, z) almost coincides with the top boundary of the anomalous structure (Figure 3c). Figure 3d and Figure 3e depict the response of the same model, but with 20% Gaussian noise added. The result of migration is shown in Figure 3f. The image becomes a little narrower, but still delineates the top boundary of the inclusion. Figure lc presents a 2-D model with three conductive rectangular inserts with resistivities 0.5 ohm^m within a 50 ohm^m background. These individual conductive bodies cannot be seen on the pseudosection of apparent resistivity (Figure 4a), and the phase pseudosection (Figure 4b), but they can easily be imaged by migration (Figure 4c). We have the same result even if we add 20% Gaussian noise to the observed data (Figure 4d, 4e, and 4f). Thus for these models EM migration produces stable images of the top geoelectrical interfaces. RESISTMTY IMAGING In practical applications, it is very important to be able to plot not only the geometry of the boundaries, but also the resistivity distribution. We discuss here the technique for the solution of this problem, based on the analysis of the vertical distribution of the reflectivity function for the same multilay­ered model introduced in the previous section. We begin our analysis with the two-layered model. We have shown above a. Distance ( m ) b. Distance (m) FIG. 1. Three 2-D resistivity models used to illustrate imaging geoelectric boundaries by EM migration. (a) Resistivity model of a step-wise structure. (b) Resistivity model with a locally conductive rectangular insert (resistivity of the inclusion is 0.5 ohm^m, resistivity of the host rock is 50 ohm • m). (c) Resistivity model with three conductive rectangular inserts (resistivity of the inclusions is 0.5 ohm • m, resistivity of the host rock is 50 ohm • m). Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ F re q u e n c y -d om a in EM M ig ra t io n 673 that the reflectivity function Pe.m« (x, z, t») at the depth of a geoelectrical boundary_ is equal to the true coefficient of reflectivity p t = (Vo-7 - Vu2)/(\'/cr7 + VcrT). Therefore the resistivity of the underlying layer p? can be calculated as where p 2 f[l + PfJJWaU, 2)]! : = \ [1 - Pt\A/aU, *)]j PI, I EMa ( *, ^) j S Pe,M<i(*^» wy)- (58) (59) On the other side, we know that stacking the migrated reflectivity function Pt.Mmt*, z, wy) for a spectrum of frequencies results in positive reinforcement at interfaces and destructive interferences elsewhere. We can introduce a nor­malized stacked migrated reflectivity function p ?) as _ 1 V P EMm(x,Z,tOj) PEMrn(x< z) M 2j In _ IV ■, ,., .M (60) N IP EMm(X, Z, U>;)| and calculate the migrated apparent resistivity pm (x, z) as Fig. 2. Apparent resistivity (a) and (d), EY phase pseudosections (b) and (e), and migration images (c) and (f) for the step-wise structure [model (a) in Figure l] without (left panel) and with (right panel) 20% Gaussian noise. Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ 674 Z h d a n o v et al. P m(x, Z) =' Pe,jZ) P ~ Pfc\Afm(*' 2)P EMAX, Z)\ (6 1) The migrated apparent resistivity pm (x, z) is equal to the resistiv­ity of the second layer at the interface and is equal to the re­sistivity of the background p[ elsewhere, because (3 E,Mm(x>z) -»■ 1 at the interfaces and 3£,*/«(*> z)-*0 elsewhere. Here, we use two reflectivity functions, apparent and migrated, to obtain a stable image of the geoelectric boundary and the resistivity contrast at the boundary simultaneously. The resistivity of the first layer p, can be found from the observed electromagnetic fields using the classical formulas. For example, in the case of plane-wave excitation we can use the admittance Y(W) estimated for a high enough frequency: p i = l/2o)n,0/?e2y(io) = l/2a)jjt0/w2y(<«>). (62) The algorithm described admits a simple generalization for the case of a multilayered background geoelectrical cross section. Indeed in this case the procedure of visualizing geoelectrical boundaries is realized successfully in a downward -40000 0 40000 Distance (m) -40000 0 40000 Distance (m) Fig. 3. Apparent resistivity (a) and (d), Ey phase seudosections (b and (e), and migration images (c) and f) for the model with a locally conductive rectangular insert [model (b in Figure l] without (left panel) and with (right panel) 20% Gaussian noise. Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ F re q u e n c y -d om a in EM M ig ra t io n 675 direction, and in each stage a band of frequencies is chosen at which the field penetrates the layers studied. In the first stage a specific electric resistivity of the first layer p, is evaluated by a standard formulas like equation (62). In the second stage, an analytic continuation (or migration) is made into the medium featuring the electric resistivity p,, and the reflectivity function fiE.Ma(x> z> w)> and the migrated reflectivity function zi w) are calculated. From the local maximum of the stacked reflectivity functions, the position of the top of the second layer and its specific electric resistivity p2 are evaluated. Then the migration and analytic continuation procedures are repeated but with a new parameter for the background me­dium p2. From the space-frequency distribution of the conven­tional and migrated reflectivity functions we find the top of the third layer and its resistivity, and so on. We can also generalize equation (61) for the model with the inhomogeneous (but slowly horizontally varying) background resistivity p„(x,z): -40000 0 40000 Distance (m) -40000 0 40000 Distance (m) Fig. 4. Apparent resistivity (a) and (d), Ey phase pseudosections (b) and (e), and migration images (c) and (f) for the model with three conductive rectangular inserts [model (c) in Figure l] without (left panel) and with (right panel) 20% Gaussian noise. Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ 676 Z h d a n o v et al. [1+ P£.M,n(*.z)PE.M»(*.*)1 [1 -P£JMm(Ar,z)P£JWa(-<C.2)]J 1 ■ 6 3 :■ Here, stacking of the reflectivity functions is done according to expression (56) at the stacking interval (t>>,-0, ,,) where <i>y0 is the frequency corresponding to the skin depth z = The last formula produces the resistivity of the underlying layer in the location of the interfaces and is equal to the background resistivity elsewhere. The main difficulty in realizing this approach is determin­ing the background resistivity distribution p„(jc, z). This problem can be solved, at least approximately, by various algebraic or differential transformations of the apparent resistivity curve. The most suitable transformation seems to be Niblett or Bostick transform (Niblett and Sayn- Wittgenstein, 1960; Berdichevsky and Zhdanov, 1984) which can be written as: I d I z \ 2 + m P„(*,z) = l/-^-j =Pa^, (64) where m = d log pa/d log VT, T = 2tt/oj and z is defined as effective depth of penetration: z = V,7'p(,/2'7T|AU. To obtain the background resistivity distribution mode1 with the slow horizontal variation, we can apply the electro­magnetic array profiling (EMAP) technique (Bostick, 1986, Torres-Verdin and Bostick, 1992) and apply some spatial filtering to the observed data. After that we calculate the conventional apparent resistivity and then recalculate it into the background resistivity using the Niblett transform (64). This is the first stage of our rapid imaging technique. In the second stage, we apply a migration transform and determine the migration apparent resistivity using formula (63). We will illustrate resistivity imaging for a model, containing one resistive and one conductive prismatic inclusions in a homogeneous background (Figure 5a) and in two layered background media (Figure 5b). The models are excited by a vertically propagating plane wave. The bottom parts of Figure 5 present the results of the resistivity imaging by migration (the plots of migrated apparent resistivity pm (x, z)). We can clearly see on these images the top boundaries of the geoelec­trical inhomogeneities, the correct values of the resistivities of the background and the inhomogeneities and the position of the top of the second layer on Figure 5d. We now examine the response of frequency-domain EM migration to host resistivity p/, errors. To make this analysis more clear, we consider again the simple model with one conductive body in the host medium with the resistivity yV: = 250 ohm • m (Figure 6a). The results of migration resistivity imaging are shown in Figure 6b. The migration image, calculated for correct background resistivity p„ = 250 a. b. Distance (m) Distance (m) Fig. 5. 2-D resistivity models used to illustrate resistivity imaging by EM migration, (a) Resistivity model with one resistive ( Pn = 2000 ohm • m) and one conductive (pi2 = 1 ohm • m) rectangular inserts in the host medium with the resistivity ph = 20 ohm • m. (b) One resistive ( pfl = 200 ohm • m) and one conductive ( pi2 = 2 ohm • m) rectangular inserts in the two layered host medium (resistivity of the first layer is pftJ = 20 ohm • m and of second layer is p*2 = 200 ohm • m). (c) Migration resistivity image for model (a), (d) Migration resistivity image for model (b). Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ F re q u e n c y -d om a in EM M ig ra t io n 6 7 7 . : •' •• ■- : :••••-: • ■ • ■ ® \ •> i „ &k & jgfcj; weWS* ■■ ' , ! ................■ ' : ■ ............,.„v... ■ % > . A% ' > -V s .....................■ : : - ■ , • Fig. 6. Study of the response of frequency domain EM migration to host resistivity p„ errors, (a) Resistivity model with a locally conductive rectangular insert (resistivity of the insert pj = 0.5 ohm • m, host resistivity p* = 250 ohm • m). (b) Migration resistivity image, calculated for the correct background resistivity p„ = 250 ohm • m. (c) Background resistivity p„ 60% higher (400 ohm • m) than the actual medium resistivity p/, causes "overmigration." (d) Background resistivity p„ 30% lower (180 ohm • m) than the actual medium resistivity p causes "undermigration.' (e), (f) Distortion of the migration resistivity image is stronger when, p„ is three times higher than pA [750 ohm • m (e)] or is three times lower than ph [85 ohm • m, (f)]. (g), (h) Migration resistivity image is completely destroyed when p„ is more than 10 times higher than py, (g) or is 10 times lower than pft (h). Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ 678 Z h d a n o v et al. ohm • m, clearly outlines the top boundary of the body and gives correct estimate of it's resistivity. Now we increase the background resistivity ( p„ ) to 400 ohm • m. The effect on the migration image is very clear-the left-hand and right-hand edges of the conductive structure go up (Figure 6c), while the central part of the image still describes well the location and the resistivity of the conductive body. if we decrease the background resistivity to p„ = 180 ohm • m the left-hand and right-hand edges of the conductive structure go down ((Figure 6d), but still the posi­tion of the central part of the image is quite correct. This result is very similar to one that takes place in seismic migration (Yilmaz, 1987). Actually in the first case we observe the "overmigration" effect when the resistivity of the background is 60% higher than the medium resistivity ph. In the second case, we observe the "undermigration" effect when the resis­tivity of the background is 30% lower than the medium resistivity p;,. As we amplify the errors in the background resistivity p„ (as it is shown in Figure 6e and 6f, where we used p„ of 750 ohm l m and 85 ohm*m correspondingly) the distor­tion level is increasing. However, only if we intensify the errors in the host resistivity determination significantly (to ten times or more), is the image significantly distorted, as it is shown in Figure 6g and 6h. The image is distorted faster in the case of decreasing the background resistivity than in the case of its increasing, because the highly resistive rocks are still transpar­ent to the EM field, while the conductive rocks absorb the EM field because of the skin effect. Thus, imaging of resistive targets in conductive background is hard. We can conclude that a correct estimation of the back­ground resistivity for migration is important, but the errors in p„ do not destroy the image dramatically as long as they are within reasonable limits (no more than one order of magnitude of the host resistivity: 0.2 < |p„/p*| < 5). IMAGING THE NORTH AMERICAN CENTRAL PLAINS CONDUCTMTY ANOMALY The North American Central Plains conductivity anomaly (known as NACP) was first discovered in the late 1960s (Reitzel et al., 1970) by a geomagnetic depth sounding (GDS) array. Jones and Craven (1989) conducted extensive studies of NACP using GDS, seismic and gravity data. We have calcu­lated the forward response of NACP using the finite-element code discussed in Wannamaker et al. (1987). The input model was slightly simplified with respect to the Jones and Craven original model to facilitate the assessment of the capability of the migration scheme to image a subsurface inhomogeneity such as that of NACP. The input model is illustrated in Figure 7a. Calculated field values were interpolated to create an equidistant set of data points with a spacing of 1.5 km. The E-polarization response was calculated at 68 periods over the four decades 0.l-l04 s and was interpolated for contouring. The resulting image of the migrated field is illustrated in Figure 7b. it is clear that separate bodies in the model are not resolved, but the general geometry of the model and its resistivity are imaged with good quality. Note, that this result corresponds very well to the result of the inversion of the synthetic NACP structures, obtained by the rapid relaxation inversion (RRI) method for TE data (Nong et al., 1993). We also studied the TM data for the same model. The TM field response was calculated at 68 periods over the four decades 0.l-l04 s as well. The migration image for TM data is presented in Figure 7c. It is possible to resolve all conductive bodies on the TM migration image, although the conductivity is underestimated by TM mode migration. Note, that in this example we use the simplified model of NACP structure without the surface conductor, which makes it possible to resolve the TM mode data. The real NACP response does not show an anomaly in the TM mode. One can expect that the joint TE and TM mode migration will produce the more accurate image, like in the case of RRI inversion (Nong et al., 1993). CONCLUSIONS Migration and analytic continuation make it possible to obtain a quick first image of the geoelectrical cross-section, provided one has available continuous profile electromagnetic observations on the surface of the earth, which are phase- referenced. It is important to notice that the computational efforts in this case are comparable to forward modeling. The results of imaging could be used either as the semiqualitative estimation of the geoelectrical model, or as a starting model for more comprehensive inversion algorithms. The practical application of electromagnetic migration re­quires addressing" several challenges: 1) Data must be spatially dense and must be collected in a manner that preserves intersite phase relations. This does not mean that all sites must be operated simultaneously, but it does require sufficient spatial overlap between separate deployments of instruments that intersite trans­fer functions can be calculated for all measured fields (see Egbert and Booker, 1989). 2) One must be able to separate the downgoing and upgoing (scattered) fields. In principle, this is accomplished for MT array data (Berdichevsky and Zhdanov, 1984; Zhdanov, 1988). For controlled sources, direct subtrac­tion of the primary field from the observed signal is appropriate to obtain the scattered (upgoing) field, but is subject to errors in measurement of the signals and source parameters, and prediction of the primary field. Despite these difficulties, EM migration clearly has potential advantages. For instance, it should be possible to quickly generate migrated images in the field and use them to optimize instrument deployments. In addition, migrated images for complicated source-receiver geometries and complex structure should be possible when other methods are not computation­ally feasible. Also, the migrated image can be used as the first approximation of the subsurface structure in an inversion scheme, based on more sophisticated forward modeling and inversion algorithms. ACKNOWLEDGMENTS Financial support for this work was provided by the National Science Foundation under grant No. EAR-9403925, and by the Office of Basic Energy Sciences of the United States Depart­ment of Energy. We also thank the Consortium for Electro­magnetic Modeling and Inversion at the Department of Geol­ogy and Geophysics, University of Utah, for the support of this work. We wish to thank P. Wannamaker for providing the code Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ F re q u e n c y -d om a in EM M ig ra t io n 679 -60000 -40000 -20000 O 20000 -40000 60000 Distance (m> Distance (m) 2.3 A3.8 83.0 126.3 167.3 208.8 230.0 Distance (m) FIG. 7. The North American Central Plains conductivity anomaly. (a) A simplified 2-D resistivity model of the NACP anomaly. (b) Migration resistivity image (TE mode). (c) Migration resistivity image (TM mode). for forward modeling and for fruitful discussions and com­ments. We are grateful to A. Tripp for correcting the manu­script and to 0. Portniaguine for help with the visualization code. We particularly thank the reviewers for their helpful and constructive review of this paper. REFERENCES Barnett, C. T., 1984, Simple inversion of time-domain electromagnetic data: Geophysics, 49 925-933. Berdichevsky, M. N., and Zhdanov, M. S., 1984, Advanced theory of dee geomagnetic sounding: Elsevier Science Publ. Bosttc , F. X., Jr., 1986, Electromagnetic array profiling: 56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 60-61. Claerbout, J. F., 1985, Imaging the Earth Interior: Blackwell Sci. Publ. deGroot-Hedlin, C., and Constable, S. C., 1990, Occam's inversion to generate smooth two-dimensional models from magnetotelluric data: Geophysics, 55, 1613-1624. Eaton, P. A., 1989, 3D Electromagnetic inversion using integral equations: Geophys. Prosp., 37, 407-426. Eaton, P. A., and Hohmann, G. W., 1989, A rapid inversion technique for transient electromagnetic soundings: PEPI, 53, 384-404. Egbert, G. D., and Booker, J. R., 1989, Multivariate analysis of geomagnetic array data I: The response space: J. Geophys. Res., 94, 14 227-14 248. Fisher, E., McMechan, G. A., Annan, A. P., and Cosway, S. W., 1992, Examples of reverse-time migration of single-channel, ground- penetrating radar profiles: Geophysics, 57,577-586. Fisher, E., McMechan, G. A., Gorman, M. R., Cooper, A. P. R., Aiken, C. L. V., Ander, M. E., and Zumberge, M. A., 1989, Determination of bedrock topography beneath the Greenland ice sheet by three-dimen­sional imaging of radar sounding data: J. Geophys. Res., 94,2874-2882. Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ Z h d a n o v et al. Gazdag, J., 1978, Wave equation migration with the phase shift method: Geophysics, 43, 1342-135 1. Jones, A. G., and Craven, J. A., 1990, The North American central plains conductivity anomaly and its correlation with gravity, mag­netic, seismic, and heat flow data in Saskatchewan, Canada: Phys. Earth Planet. Inter., 60, 169-194. Hogan, G., 1988, Migration of ground-penetrating radar data: A technique for locating subsurface targets: 58th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 345-347. Lavrent'ev, M. M., Romanov, V. G., and Shishatsky, S. P., 1980, Ill-posed problems of mathematical physics and analysis: Nauka (in Russian). Lee, K. H., Liu, G., and Morrison, H. F., 1989, A new approach to modeling the electromagnetic response of conductive media: Geo­physics, 54, 1180-l 192. Lee, K. H., and Xie, G., 1993, A new approach to imaging with low-frequency electromagnetic fields: Geophysics, 58, 780-796. Lee, S., McMechan, G. A., and Aiken, L. V., 1987, Phasefield imaging: The electromagnetic equivalent of seismic migration: Geophysics, 52, 679-693. Macnae, J., and Lamontagne, Y., 1987, Imaging quasi-layered conduc­tive structures by simple processing of transient electromagnetic data: Geophysics, 52, 545-554. Madden, T. R., and Mackie, R. L., 1989, Three-dimensional magne­totelluric modelling and inversion: Proc. IEEE, 77, 2, 318 -332. Nabighian, M. N., 1979, Quasi-static transient response of a conduct- i1n7g0 0h-a1l7f-0s5p.ace-An approximate representation: Geophysics, 44, Nekut, A., 1994, Electromagnetic ray-trace tomography: Geophysics, 59, 371-377. Niblett, E. R., and Sayn-Wittgenstein, C., 1960, Variation of electrical conductivity with depth by the magnetotelluric method: Geophysics, 25, 948-1008. Nong, Wr., Booker, J. R., and Smith, J. T., 1993, Rapid two­dimensional inversion of Coprod2 data: J. Geomag. Geoelectr., 45, 1073-1087. Oldenburg, D., and Li, Y., 1993, Inversion of induced polarization data: 63rd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 396-399. Oldenburg, D. W., Li, Y., and Eliis, R., 1994, Joint interpretation of DC, IP, magnetic, airborne EM, geological, and mineralization data: A case history from Mt. Milligan, 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 512-515. Oristaglio, M., Wang, T., Hohmann, G., Tripp, A., 1993, Resistivity imaging of transient electromagnetic data by conjugate-gradient method: 63rd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 347-350. Pellerin, L., Johnston, J., Hohmann G., 1993, Three-dimensional inversion of electromagnetic data: 63rd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 360-363. 680 Reitzel, J. S., Gough, D. I., Porath, D. I., and Anderson, III, C. W., 1970, Geomagnetic deep sounding and upper mantle structure in the western United States: Geophys. J. Roy. Astr. Soc., 19, 213-235. Samarsky, A. C., 1984, Theory of differencing schemes: Nauka (in Russian). Smith, J. T., and Booker J. R., 1988, Magnetotelluric inversion for minimum structure: Geophysics, 53, 1565-1576. ----------1991, Rapid inversion of two- and three-dimensional magne­totelluric data: J. Geophys. Res., 96, 3905-3922. Torres-Verdin, C., and Bostick, F. X., Jr., 1992, Principles of spatial surface electric field filtering in magnetotellurics: Electromagnetic array profiling (EMAP): Geophysics, 57, 603- 622. Torres-Verdin, C., and Habashy, T., 1994, Rapid 2.5-dimensional forward modeling and inversion via a new nonlinear scattering approximation: Radio Science, 29, 105 l-1079. Tripp, A. C., and Hohmann, G. W., 1993, Three-dimensional electro­magnetic crosswell inversion: IEEE Trans. on Geoscience and Remote Sensing, 31, 121-126. Wannamaker, P. E., Booker, J. R., Jones, A. G., Chave, A. D., Filloux, J. H., Waff, H. S., and Law, L. K., 1989, Resistivity cross-section through the Juan de Fuca subduction system and its tectonic implications: J. Geophys. Res., 94, 14 127-14 145. Wannamaker, P. E., Stodt, J. A., and Rijo, L., 1987, A stable finite-element solution for two-dimensional magnetotelluric model­ing: Geop. J. Roy. Astr. Soc., 88, 277-296. Wessel, P., and Smith, W. H. F., 1991, Free software helps map and display data: EOS Trans. AGU, 72, 441, 445-446. Xiong, Z., and Kirsch, A., Three-dimensional earth conductivity inversion: 1992 J. Comput. and Applied Math., 42, 109-121. Yilmaz, O., 1987, Seismic data processing, Series: Investigation in Geophysics, No. 2: Soc. Expl. Geophys. Zhang, J., Mackie, R., and Madden, T., 3-D resistivity forward modeling and inversion using conjugate gradients: 64th Ann. Inter­nat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 377-380. Zhdanov, M. S., 1988, Integral transforms in geophysics: Springer- Verlag. Zhdanov, M. S., and Booker, J., 1993, Underground imaging by electromagnetic migration: 63rd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 355-357. Zhdanov, M. S., and Frenkel, M. A., 1983a, The solution of the inverse problems on the basis of the analytical continuation of the transient electromagnetic field in reverse time: J. Geomag. Geoelectr., 35, 747-765. - 1983b, Electromagnetic migration: in Hjelt, S. E., Ed., 37-58, The development of the deep geoelectric model of the Baltic Shield, Part 2, Univ of Oulu, Oulu. Zhdanov, M. S., and Keller, G. V., 1994, The electromagnetic methods in geophysical exploration: Elsevier. Zhdanov, M. S., Matusevich, V. Yu., and Frenkel, M. A., 1988, Seismic and electromagnetic migration: Nauka (in Russian). APPENDIX A SEPARATION OF EM FIELD INTO DOWNGOING AND UPGOING COMPONENTS Here we present the simple technique for field separation in the 2-D case following the approach developed by Berdichev­sky and Zhdanov (1984). Consider a 2-D model with the constant background distri­bution of conductivity a n = const. Anomalous conductivity is concentrated in some closed domains or in some layer below the Earth's surface. In the layer with constant conductivity crn, the E-polarized electric field can be represented as a sum of downgoing and upgoing fields, ey(kx, co) exp (-vz) + e"(kx, co) exp (vz), (A-1) where ey (kx, z, oo) is the 2-D Fourier transform of the Ey ( x, z, co) component, ey (kx, co) and ey (A:x, co) are the spectra of the downgoing and upgoing components of the field on the surface of the earth, and v = V(kx - /co|jL0(jn), Re (v) > 0 is a wavenumber in (k, co) domain. From Maxwell's equations the 2-D Fourier transform hx (kx, z, co) of the Hx(x, z, co) component is X [e"{kx, a>) exp (vz) - ey(kx, a>) exp (-hz)]. (A-2) While equations (A-l) and (A-2) give KOfJLo ey(kx, 0, co) +------hx(kx, 0, co) zco|xo ey(kx, 0, co)---------hx{kx, 0, co; (A-3) Thus, applying an inverse Fourier transform to both sides of the equations (A-3) gives the downgoing and upgoing compo­nents of the electric field on the surface of the earth. Downloaded 26 May 2010 to 155.97.11.183. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/