Time domain electromagnetic migration

This thesis develops a new approach to processing transient electromagnetic data. This approach is called time domain electromagnetic migration (TDEMM). The migration method is specifically designed tp process multidimensional time-domain electromagnetic data in order to obtain a resistivity image of subsurface structure which clearly indicates the spatial location of the targets and estimates their anomalous conductivity. The theoretical foundations of migration, a model study illustrating the method's stability and resolution, and an application of the migration procedure to real data interpretation are subsequently considered in the work.

T I M E D O M A I N E L E C T R O M A G N E T IC M I G R A T I O N by Portniaguine The University of Utah in partial fulfillment of t h e requirements for the degree of in Geophysics The University of Utah August 1995 TIME DOMAIN ELECTROMAGNETIC MIGRATION by Oleg Nikolaevich Portniagui ne A thesis submitted to the faculty of the Master of Science 1D Geophysics Department of Geology and Geophysics T he © Copyright © Oleg Nikolaevich Portniaguine 1995 All Rights Reserved THE U N I V E R S I T Y OF UTAH G R A D U A T E SCHOOL S U P E R V I S O R Y C O M M I T T E E A P P R O V A L supervisory satisfactory. UNIVERSITY GRADUATE SCHOOL SUPERVISORY COMMITTEE APPROVAL of a thesis submitted by Oleg Portniaguine This thesis has been read by each member of the following supeIVisory committee and by majority vote has been found to be satisfactory. Chair: Michael S. Zhdanov Gerard T. Schuster L CI~b~ Alan C. Tripp , , THE U N I V E R S I T Y OF U T A H G R A D U A T E S C H O OL F I N A L R E A D I N G A P P R O V A L form final UNIVERSITY UTAH GRADUATE SCHOOL FINAL READING APPROVAL To the Graduate Council of the University of Utah: I have read the thesis of Oleg Nikolaevich Portniaguine in its final fonn and have found that (1) its format, citations, and bibliographic style are consistent and acceptable; (2) its illustrative materials including figures, tables, and charts are in place; and (3) the fmal manuscript is satisfactory to the supervisory committee and is ready for submission to The Graduate School. Date Michael S. Zhdanov Chair: Supervisory Committee Approved for the Major Department Approved for the Graduate Council Ann W. Hart Dean of The Graduate School A B S T R A C T thesis called migration method is specifically designed to process multidimensional time-domain elec­tromagnetic data in order to obtain a resistivity image of subsurface structure which clearly indicates the spatial location of the targets and estimates their anomalous conductivity. The theoretical foundations of migration, a model study illustrating t h e method's stability and resolution, and an application of the migration procedure to real data interpretation are subsequently considered in this work. ABSTRACT This t hesis develops a new approach to processing transient electromagnetic data. This approach is caJ1ed time domain electromagnetic migmlion (TDEMM). The mi­gration structUl'c spat ial theo retical the method 's procedUl'c data interpretat ion a re t his C O N T E N T S A B S T R A C T iv L I S T O F F I G U R E S vi A C K N O W L E D G E M E N T S viii I N T R O D U C T I O N 1 1.1 characterization T H E O R Y O F E M M I G R A T I O N 6 realization algorithm M O D E L S T U D Y 3.1 method 3.2 3.3 field 3.4 39 3.5 Resolution of migration method 42 3.6 42 49 E M M I G R A T I O N I N W A S T E S I T E C H A R A C T E R I Z A T I O N 4.1 4.2 Survey design and data description C O N C L U S I O N A P P E N D I C E S A . A P P A R E N T R E F L E C T I V I T Y F U N C T I O N 71 M I G R A T I O N A P P A R E N T R E S I S T I V I T Y R E F E R E N C E S CONTENTS ABSTRACT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. IV LIST OF FIGURES.. . . . . . . .. . . . . . .. . . . . . . . . . . . .. . . . . . .. . . . . . . . . . . . . . . . . . . VI ACKNOWLEDGEMENTS ............................................... Vlll INTRODUCTION. . .. . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . .. . . . . . .. . . . . . . . . . . . . 1.] Time domain electromagnetic soundings ...... 1 1.2 Introducing time domain electromagnetic migration . . . . . . . . .. 2 1.3 Waste site characteri7-ation . . . . . . . . . . . . . . . . . . . . . . .. 5 THEORY OF EM MIGRATION. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Downward continuation in reverse time . . . . . . 6 2.2 Resistivity imaging by TDEM migration . . . . . 11 2.3 Numerical reali7-ation of the migration algorithm. 20 2.4 Testing of the algorithm . . . . . . . . . . . . . . 22 MODEL STUDY. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 27 3.7 Modeling strategy and method. . . . . . . . . Model with one local conductive body The separation of primary and secondary field Migration of noisy data . . . . . . . . . . . . . . . . . . . . Resolution of conductive and resistive bodies . Background resistivity influence on the EM migration image 27 28 28 EM MIGRATION IN WASTE SITE CHARACTERIZATION. ....... 56 4.3 INEL Waste Complex Characterization Project .................. . Application of migration for data interpretation ........... . 56 57 60 CONCLUSION. . . . . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . .. . . . . . .. . . . . .. . . . . . . . . .. 70 APPENDICES 71 A.APPARENT REFLECTIVITY FUNCTION ........................ 71 B. MIGRATION APPARENT RESISTIVITY ......................... 75 REFERENCES ............................................................ 78 L I S T O F F I G U R E S t h e field t h e wave, a) zone, slingram mode . . t h e t h e t h e cable field. Aperture dots) field. Grid dots) 7 Local conductive body in a homogeneous media 29 7 7 7. . 7 7 7 0.2/is separation, sec­ondary field is directed differently at symmetric points 38 field's maximum data separation separation separation separation separation LIST OF FIGURES 1 Approximation of the primary field propagation by the plane wave. a) Primary field propagation in the far zone. b) Primary field in slingram mode. ................................. 14 2 Surface response of an electrical field simulated for a buried cable. .. 15 3 Migrated field of the electrical cable response. Cross shows the location of the cable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 16 4 Secondary field for the model with one horizontal layer (solid line) and migration kernel (dashed line) . . . . . . . . . . . . . . . . . . . . .. 24 5 Analytic solution (solid line) and numerically migrated field. Aperture is 1000 m (dashed line) and 500 m (dots). ............... 25 6 Analytic solution (solid line) and numerically migrated field. Grid spacing is 40 m (dashed line) and 100 m (dots). 26 media. . . . . . . 8 Time derivative of primary field for the model of Figure 7. 30 9 Time derivative of secondary field for the model of Figure 7. 31 10 Migrated time derivative of magnetic field for the model of Figure 7.. 32 11 Migrated electric field for the model of Figure 7. . . . . . 33 12 Apparent reflectivity function for the model of Figure 7. 34 13 Migration apparent resistivity for the model of Figure 7. 35 14 Secondary field at 1s .. . . . . . . . . . . . . . . . . . 37 15 Primary/secondary field separation. a) The primary field does not change along the profile above the horizontally layered media, b) sec-ondary points. . . 16 Secondary field with 5% random noise to field's maximum. 40 17 Migration image produced from noisy data. . 41 18 Conductive bodies with 100 m separation. . 43 19 Migration image for conductive bodies with 100 m separation. 44 20 Conductive bodies with 200 m separation. ........... 45 21 Migration image for conductive bodies with 200 m separation. 46 22 Conductive bodies with 300 m separation. ........... 47 Migration conductive separation separation bodies sepa­ration resistive depths conductive at different depths 53 28 Migration with incorrect background resistivity (80 Ohm-m) 54 29 Migration with incorrect background resistivity (300 Ohm-m) 55 30 Location of anomalous objects in Cold Test Pit 58 31 Slingram mode TDEM survey design at Cold Test Pit 59 32 Cold Test Pit TDEM data (profile 0s) 61 33 Cold Test Pit TDEM data (profile 35s) 62 34 "Smoke ring" image of 0s profile 63 35 Migration apparent resistivity for Cold Test Pit (profile 0s) 64 36 Migration apparent resistivity for Cold Test Pit (profile 35s) 65 37 Apparent resistivity map of Cold Test Pit at the depth 3 m 67 38 Apparent resistivity map of Cold Test Pit at the depth 5 m 68 39 Horizontal slices of Cold Pit resistivity image 69 vii 23 rvLigratioll image for conducti ve bodies with 300 m separation. 48 24 Resistive and conductive bodies with 300 m separation. . . . . 50 25 Migration image for resistive and conductive bodjes with 300 m sepa-ration. . . . . . . . . . . . . . . . . .. ....... 51 26 Conductive and resist'ive body located at different depths. 52 27 Migration image for model with condudive and resistive bodies located at different depths. . . . . . . . . . . . . . . . . . 53 28 Migration with incorrect background res istivity (80 Ohm-m). . 54 29 Migration with incorrect background res istivity (300 Ohm-m). 55 30 Location of anomalous objects in Cold Test Pit. . . . 58 31 Slingram mode TDEM survey design at Cold Test Pit. 59 32 Cold Test Pit TDEM data (profile Os). . 61 33 Cold Test Pit TDEM data (profile 35s). . 62 34 "Smoke ringn image of Os profile. .... 63 35 Migration apparent resistivity for Cold Test Pit (profile Os). 64 36 Migration apparent resistivity for Cold Test Pit (profile 355). 65 37 Apparent resistivity map of Cold Test Pi t at the depth 3 m. 67 38 Apparent resistivit.y map of Cold Test Pit at the depth 5 m. 68 39 Horizontal slices of Cold Pi t resistivity image. . . . . . . . . 69 VII A C K N O W L E D G E M E N T S I the work. I Dr. the directing thanks with Tartaras am fruitful t h e migration discussions, I of this work. ACKNOWLEDGEMENTS T thank t he Consortium on Electromagnetic Modeling and Inversion at the De­partment of Geology and Geophysics, University of Utah, for providing support for this work . 1 also thank my adviser, DT. M.S. Zhdanov, for sharing t he ideas and directing this research. Special t hanks to my colleagues, Dr. P. Traynin for his help wilh handling and visualizing the data, and E. Tal'taras for his contribution to model calculations. I grateful to Dr. H. D. MacLean for providing the TDEM data and fru itful discussions of the results of mi~;ration and resistivity imaging, and Dr. A.C.Tripp for di scussions, comments, and corrections. 1 appreciate Dr. G.T. Schuster for the ideas and suggestions about future devel­opment 1 . I N T R O D U C T I O N 1 . 1 T i m e d o m a i n e l e c t r o m a g n e t i c s o u n d i n gs t r a n s m i t t e r design. The transmitter is usually a square or circular loop or coil, excited by an abruptly terminated current. Surface and borehole surveys can be done in many different ways. In the fixed transmitter technique the induced magnetic field is mea­sured at many receiving points, located at a variable distance from the transmitter. A second survey possibility, the slingram mode, requires the receiver and transmitter to be a fixed distance apart, moving together along the profile. Transient electromagnetic d a t a are most often interpreted using conventional one-dimensional ID) inversion techniques, or approximate, fast imaging techniques. ID fitting ID structure to t h e real data in a least-square sense (Farquharson and Oldenburg, 1993). It is conventionally applied to every survey point. Combining ID images together forms the resistivity-cross section. Fast resistivity imaging techniques have been developed in (Eaton and Hohmann, 1989; Macnae and Lamontagne, 1987; Barnett, 1981). The majority of these papers have been based on equating the transient response, measured at the surface of the earth, to t h e electromagnetic (EM) field of current filament images of t h e source. This approach was originated in t h e pioneering work of Nabighian (1979), which describes t h e transient currents diffusing into the earth as a system of "smoke rings" blown by 1. INTRODUCTION 1.1 Time domain electromagnetic soundings Transient electromagnetic methods are widely used in exploration geophysics and environmental studies. These are inexpensive and reliable nondestructive methods of subsurface and borehole investigations. Transient electromagnetic soundings can indicate the location of structures with anomalous conductivities with respect to the surrounding medium. Classification of transient sounding methods is based on transmitter type and sur­vey transmitter. transmitter fixed data one­dimensional (lD) The theory of 1D inversion is based on fitting the theoretical response from a 1D the 1D Hohmann, the filament the the the t h e t r a n s m i t t i n g 2 1 . 2 I n t r o d u c i n g t i m e d o m a i n e l e c t r o m a g n e t ic m i g r a t i o n This time domain electromagnetic migration. formulated originally in the pioneering works of Zhdanov and Frenkel (1983) and they have been developed in several publications by Zhdanov and coauthors (Zhdanov and Frenkel, 1983a,b; Zhdanov and Podchufarov, 1988; Zhdanov et al., 1988; Zhdanov and Keller, 1994). t h e major advantage in comparison with the conventional ID inversion and the "smoke ring" approaches. In this thesis a 2D migration algorithm is considered, but it can be readily generalized to the 3D case. ID t h a t the "smoke ring" imaging could be used for an estimation of background resis­tivity. migration resistivity imaging. We then tested the migration method on a set of numerical models and real TDEM data. et al., 1988; Claerbout, 1985) but differs in that electromagnetic migration is done on the basis of a diffusion equation rather than a wave equation. Seismic migration is based on downward continuation of the source field and the scattered field on the basis of a wave equation. At zero time the phases of the incident and t h e scattered fields are typically coincident at the boundary of the scatterer, and t h e ratio of incident and scattered field amplitudes is proportional to the reflection the transmitting loop into the earth. 1.2 Introducing time domain electromagnetic migration T'his thesis focuses on an approach to processing transient EM data, based on downward data extrapolation in reverse time. This method is called domain The principles of EM migration have been formulated The migration method is a multidimensional imaging method, which is the major 1D Comparison of migration with 1D imaging is done by Traynin (1995). This work concludes that the migration method has better resolution of localized bodies. How­ever, resis-tivity. In this work we modified a preceding time domain electromagnetic migration algorithm, and developed an effective numerical method and a computer code for Electromagnetic migration is similar in many aspects to seismic migration (Zh­danov the the reflection coefficient for this boundary. The wave field backward continuation problem could be formulated as a boundary value problem for the same wave equation, but with a changed boundary condition for the vertical derivative of the seismic field. Reversing t h e time affects only the boundary condition, but not the form of t h e wave equation. Thus, upward and downward continuation of the wave field do not differ much from one another. Physically, propagation of seismic waves in rocks is characterized by some attenuation, with a large amplitude decrease caused by geometrical spreading and losses from reflections. fact, the wave equation which is successfully used in most cases, does not contain an attenuative factor at all. That is why it is relatively easy to backward propagate the seismic field. was demonstrated in the literature t h a t seismic migration can be t r e a t e d as the first iteration in t h e process of minimizing the waveform misfit functional (Tarantola, 1987). This also can be shown for the electromagnetic high-frequency field (Cai et al., 1995). We expect that the same type of generalization can be developed for t h e quasi-stationary EM field as well. However, this consideration lies outside the framework of this thesis. Electromagnetic field propagation is characterized by strong attenuation. For low frequencies the attenuation dominates over the wave propagation, which provides the basis for a quasi-stationary approximation that describes low-frequency electromag­netic field propagation as a pure diffusion process. Although the same downward continuation could be done on the basis of the diffusion equation for the electromag­netic field, such a transformation would be unstable, for it dramatically amplifies high frequencies. To avoid instability, electromagnetic migration is defined as down­ward continuation of t h e diffusive field in reverse time. This transformation correctly reconstructs only the phase of the backward propagating field, but not the ampli­tude. Migration selectively suppresses frequencies too high for penetration at the given depth, and preserves the frequencies that can penetrate to this depth. The sup­pression makes migration stable in the presense of high frequency noise, but limits t h e resolution. Analogous to the seismic case, the ratio of the backward propagated scattered electromagnetic field and the primary field is proportional to t h e reflectivity 3 field. the the from cansed In field. It that treated the (Tarantola, aI., the downward field, amplifies the the propagated the reflectivity coefficient at a geoelectrical boundary. This property is used to image geoelectrical of t h e diffusive field in reverse time, we need to solve numerically a boundary value problem for the diffusion equation. This could be done, for example, using a finite-difference approach. This approach is attractive because it allows us t o back­ward propagate the field through inhomogeneous media. We, however, have chosen a simpler approach, which is based on Green's theorem. According to Green's theorem, t h e field at any point inside some domain could be expressed via t h e values of t h e field at the domain's boundaries. Approximating the earth surface by the plane and con­sidering the lower half-space as a domain for analytic continuation, one can express t h e value of the backward propagating field at any time moment at any particular depth via an electromagnetic analog of the Rayleigh integral. This integral has the same form as t h e well-known seismic Rayleigh integral, the only difference being that instead of t h e wave equation Green's function we use the diffusion equation Green's function, expressed in reverse time. Thus, the backward propagated field focuses the secondary field on the boundary of a scatterer. At this point, the migrated field should be transformed into a resistivity image of the resistivity ingomogeneity. The key step in such a transformation is estimation of the incident field, because the ratio of the incident and the scattered field gives t h e reflectivity function, which could be recomputed to give t h e anomalous resistivity. t h e surface to the lower half-space. However, for the sake of simplicity, in this work we approx­imate the primary field by the plane wave propagating in a homogeneous medium. This approach allows us to analytically estimate t h e primary field at every point in the lower half-space at the zero time moment, thus providing a simple algorithm which transforms the amplitudes of the migrated field to the migration apparent resistivity. 4 geoclectrical structures. Thus, formulating the electromagnetic migration problem as a backward propa­gation the difference to the the the field the the the boundary scattered the the anomalous resisti vi ty. It is possible numerically to downward continue the primary field from the surface the resistivity. Based on this theory, the computer code was created. The code consists of two major parts. The first part of the code takes the values of the field component on 5 t h e space, using a Rayleigh-type integral. This part of the code was tested by comparing the analytic solution for the plane wave in a two-layered model with the numerical migration code is tested on 2D numerical models. The stability and resolution of the EM migration procedure are illustrated by the results of numerical modeling. 1 . 3 W a s t e s i t e c h a r a c t e r i z a t i on t h e (RWMC) at the Idaho National Engineering Laboratory (INEL) (MacLean, 1993). TDEM t h e obtained set of profiles, intersecting the Cold Test Pit. As a result a large volume of TDEM d a t a were acquired. TDEM migration was employed to process and interpret the 3D TDEM data. TDEM migration and resistivity imaging made it possible to outline the boundaries of t h e Cold Test Pit and characterize the distribution of the anomalous conductivity, thus directly pointing to the location of the drums. Migration results were compared with "smoke ring" imaging, which was done on the same data. "Smoke ring" imaging failed to resolve separate conductive bodies, but provided a reasonable estimation of t h e background resistivity. In contrast, migration apparent resistivity maps indicate not only the thickness of capping but also t h e depth extent of t h e Pit. the surface and computes the migrated field at zero time moment at every point of comparing numerical solution. The second part of the code takes values of the horizontal component of electric field at the zero time moment in the lower half-space and computes the migration apparent resistivity, approximating the primary field by the plane wave. The whole EM 1.3 Waste site characterization EM migration was applied to the interpretation of a TDEM data set acquired at the Cold Test Pit site within the Radioactive Waste Management Complex (RWMC) The Cold Test Pit was specifically constructed to investigate the ability of TDEM systems to characterize DOE nuclear waste sites. The prime goal of the geophysical survey at the Pit was to locate conductive drums, simulating nuclear waste containers. The containers were buried at known places in the Pit. The data were obtained as a result of a high density TDEM profiling survey with the Geonics EM47 along the TDEM data the conductivity, compared the the the 2 . T H E O R Y O F E M M I G R A T I ON t h e conductivity using t h e migrated field. The last two subsections outline important theoretical aspects of algorithm realization and describe the analytic solution which was used to test the code. 2 . 1 D o w n w a r d c o n t i n u a t i o n i n r e v e r s e t i m e z = conductive [> < t h e cr(an = const, Aa(r) (Berdichevsky c-(r) = an + Ao-(r). field -* dH -> d E AH-fi0crn- = 0, AE-fi0o-n- = 0. component, P ( r , t ) , field. field, P™, component, P ° , t h e field 2. THEORY OF EM MIGRATION In this section the theory of the EM migration method is considered. The fol­lowing three subsections deal with the mathematical formulation of the migration problem, the extremal properties of the migrated field and estimation of the anoma­lous the field. important 2.1 Downward continuation in reverse time Consider a model in which the horizontal plane, 0, separates the conductive earth (z > 0) from the insulating atmosphere (z < 0). The conductivity in the earth, O"(r), is represented as the superposition of a constant normal (background) conduc­tivity, O"n canst, and an anomalous conductivity function, ~O"( r) (Berdichevsky and Zhdanov, 1984): (1) Everywhere in the earth outside the anomalous region, the electromagnetic field sat­isfies the diffusion equation: (2) For this model, we can discuss the problem of migration of any scalar component, P(f', t), of an observed electromagnetic field. Let us define the migrated field) pm, obtained from a specific scalar component, po, of the electromagnetic field observed at the earth's surface as being the field that satisfies the following conditions: 7 p m ( ? T ) l _( P°(r,T-r)z=0 forO<r<T, P (r) \ z = 0 - < Q forr<0, r >T [6) dPm(r T) APm ( f , T) fi0a = 0 for z > 0, (4) P r a ( r , r ) ^ 0 , | r \-* oo, z r - is r , 8Pm(r t) APm(r,t) + ti0<7 g\ ' ' 0 . (6) receivers, t h e n eq. (6) describes the inverse process of propagation from receivers to sources. P° time, r . This procedure is called field (Zhdanov and Frenkel, 1983; Zhdanov and Keller, 1994). field is reduced to a boundary value problem described by eqs. (3) through (5). This boundary value problem can be solved using Green's theorem and t h e Green's function for the 3D diffusion equation (Morse and Feshbach, 1953): G(r',t' \r,t) = - J y ) ^ ) 3 / a e - w ' | f ' - r | ' / 4 ( ' ' - t ) g ( « ' - 0 , (?) equation: Pm(--> ) 1 { pO(r,T - T)z=O for 0:::; T :::; T, r, T z=O= 0 for T < 0, T > T (3) Apm(-) apm(r,T) L..l. r,T - f-l00" aT =Oforz>O, (4) and (5) for 1 r 1----+ <Xl, Z > 0, where T = T - t is reverse time, and T IS the interval of electromagnetic signal recording. Note that if we exchange the reverse time, T, for ordinary time, t, in eq. (4), we have an equation which is the adjoint to the diffusion equation: A m(--> ) apm(r, t) L..l.P r,t f-loO" at = 0. (6) If the ordinary diffusion equation describes field propagation from the sources to then The problem of establishing the migrated field reduces to the downward extrap­olation of the field pO from the earth's surface into the lower half-space in the re­verse T. electromagnetic field migration It can be seen from these considerations that the calculation of a migrated field This the function G ( ) 1/2 G( , t' 1 r t) = - f-loO" -/toalr l -rI2 /4(t '-t) H(t' - t) r " 87r3 / 2 ( t' _ t )3/2 e , (7) which satisfies the equation: A C + fiQ<T- = S(f- ?)6(t - t'). (8) H(t' - (step-function): H(t - t') = 0; t ' - t <0 f - > 0 . surface) can be solved using equation (6) (Zhdanov and Keller, 1994). Let us denote t h e lower half-space as and the surface of t h e earth as Introducing a hemisphere O of large radius R we can construct a closed contour S + O around the domain D. Multiply the the left and right sides of equation by and integrate both sides over the infinite time range and over the space domain D. We then transform the volume integral over D to a surface integral over S + O using Green's formula and taking equation into consideration. Thus we have r' + g equation g in D g g = -0(7" 6G + JioO" ~~ 8(r' - t)8(t t'). Here H (t' - t) is a Heaviside excitation function (step-function): H(t _ t') = {o; t' - t < 0 1; t' - t > o. 8 (8) The problem of downward continuation of the field P (the field observed on the the D, the S. hemisphere o + 0 D. (6) G + 0 (6) P m(--<r / ,t' ) = 100 1 (PaG- ap) a - G-a dsdt, -00 5+0 n n (9) for all inside domain D. The same equation remains valid if we replace G with the sum (G + g), where 9 is an arbitrary solution to the homogeneous diffusion equation D pm(t, t') = 100 r (pa(G + g) _ (G + g) ap) dsdt. -ooJ5+0 an an Assume that the auxilary field 9 goes to zero as R goes to infinity. The S is the surface of observations, thus the normal derivative is equal to vertical derivative, so (10) If we take 9 as 9 = -G(t - t ' , t), r" r' z 9 = G + g = 0 d(G + g) _ dG dz dz' for z = 0. ( i i ) t h e field: G G Morse and Feshbach, 1953; Zhdanov et al., 1988). The vector r represents the image point of t h e migration field, and the vector r represents the integration point on t h e surface. der, G*, G t h e seismic case, a migration transformation of the electromagnetic field yields the upgoing field. t h e t h e field t h a t Y In this case the expressions for the different components of the migrated magnetic field, H®z(x, 0, £), observed along the profile X on the surface of t h e earth, will have t h e form: t h e Schnei function H™z(x',z',T-t') 47T X where the point is located symmetrically to with respect to the plane 0. Then G+g=o a(G + g) aG az = 2 az' for z 0. (11) Substituting equations (11) in equation (10), and after some algebraic operations we have the following representation for the migrated field: (12) where 0 is the adjoint to the Green's function for the diffusion equation (Morse i ' the field, i the surface. Equation (12) is the electromagnetic counterpart to the Rayleigh integral (Schnei­der, 1978; Claerbout, 1985). Just as in the seismic application, eq. (12) defines in space and normal time a field propagating towards the surface of observation (that is, upgoing waves), as can be seen from the fact that eq. (12) contains the function 0, adjoint to the Green's function of the diffusion equation. Hence, just as in the field. Let us consider the special case of the 2D model of the electromagnetic field (for example the E polarized mode) and a profile observation. We assume that the axis X coincides with the profile, and axis is orthogonal to the profile of observations. magnetic H~ ,z (x, t), the the form: m (" ') /-loa nHxz x,z,T-t = ---Z' iT ~ Hx0 z ( x,O,t ) 1 2 , 47r t' x' (t - t') (13) and analogously for the electric field: x exp fi0crr 4 (t - f) (x - x)2 + (z - z)2^ I dxdt. of t h e magnetic field. Therefore, the actual data are proportional to the time-derivative of the magnetic field variations j^H®z(xy 0, Thus we have to modify -H m ( x ' z' T-t') / " " T " * ' r f -H° (x 0 t) 1 x x exp 4 ( « - * ' ) (a;' - a;)2 + (z' - z ) 2 f cfoofa. (15) The last expression gives a possibility in the 2D case of calculating the migrated electric field from the observed vertical component of the magnetic field, since from dE™ _ OH™ and rx Q E; (x, z\ T-t') fi0 J* QJH?{X\ z\ - dx', (16) 'XL (x, z') XL coordinate of t h e left end point on the profile. Thus, from the observed vertical component of t h e magnetic field we can calculate t h e z' field E™ at the same level. 10 field: Eym ( x I ,zI ,T - t ') = -/10C-Jn-Z' iT] Ey0 (x,O,t ) 1 2 x 471" tf X (t - t') (14) Note, that typical EM equipment uses receiver loops for measuring the compo­nents the field. time­derivative ~tH~,z(x, t). modify (13): ~uati Hxm , z (x', z I ,T - tI) = /10C47J1n" Z' l'tfT ]x u~a Htx,0 z ( x, 0, t ) (t -1 t' ) x (15) field, from the second Maxwell's equation for the migrated field we have: oEm oHm y z -ax- = -/10aT- -' E;;Ci, Zl, T - t') = /10 r;; :,H';(x', Zl, T - t')dx', iXL ut (16) where (x, Zl) is a current point of electrical field calculations, and XL is a horizontal the the the migrated magnetic field at any level Zl and then determine the migrated electric E;: In a TDEM sounding one typically measures the electromagnetic response of a t h a t it is important to have a capability of converting these data into the delta-pulse current response. The conversion of electromagnetic measurements due to a step-function current to that due to delta-pulse current can be done in the process of migration. The delta-pulse response can be obtained from t h e step-pulse response by differen­tiation with respect to time. Therefore we can use equation (15) and after integration by parts obtain the time derivative of the migration field, which is equal to the mi­grated step response electrical field: H£(x', Z', T-t') ^H^"(x', z', T-t') = ff h o, t) x (17) x exp { - ^ r y [(&' x)2 + z)2] } dxdt, where H^p are step-function responses and H™* are migrated delta-pulse responses. 2 . 2 R e s i s t i v i t y i m a g i n g b y T D E M m i g r a t i on At this point we will give a description of the basic concepts in which migration imaging of geoelectric structures is based. To illustrate the "EM radiating-inhomogeneities" concept, we represent the ob­served electromagnetic fields E, H as being the sum of a primary field, EP,HP, and field, Es, Hs : E = EP + E% H = HP + H% source in the earth with the normal distribution of the electrical conductivity an(r) (say, in a homogeneous half-space) and the secondary field is due to t h e anoma­lous conductivity distribution. In other words we can treat the secondary field as 11 step-function current waveform. However we will see below that for migration imaging step­function the differen­tiation integration field: Hm8 (x' z' T - t') = JL Hmstep (x' z' T - t') = x,z , , ot' x,z " (17) X {- 4(~~;) [(x' - X)2 + (z' - Z)2]} dxdt, H~t,; p H;:,; 2.2 Resistivity imaging by TDEM migration E, jj field, EP, jjP, a secondary field, jjjs, jjs : (18) where the primary field is defined as being that field which is generated by a spe­cific conductivity O"n(r) the anoma­lous t h e t h e diffusion equation, this field does not feature such simple and effective geometrical properties as the Claerbout upgoing wave, which ensures a direct reconstruction of reflection boundaries by the location of backpropageted wavefront at the instant t' = 0. We will see nevertheless that the migration electromagnetic field also can be used for detecting local inhomogeneities in a conducting medium, and for determination of geoelectrical boundaries. t h e different electromagnetic t h e time. When the time reaches zero, the extrapolated field is at the location of the scatterer. Thus, in the case of the secondary field this downward extrapolation in reverse time will result in the focusing of the migrated field at zero time at the location of t h e source. complicated than the delta pulse. Since the propagation is dispersive, these sources "switch on" instantaneously, and the field builds up to t h e maximum with some delay, related to the time of propagating the primary field from the earth's surface to the geoelectrical inhomogeneties. That is why we have two options: 1) to reconstruct the migration field in the lower halfspace not at zero time, but at a retarded time, or 2) to use for the downward extrapolation not the actual background conductivity, but instead some effective migration conductivity o~m proportional in the general case to t h e normal conductivity am = jcrn, where 7 is called the migration constant. De­tailed analysis of geoelectrical models shows t h a t in t h e second case we can select the 12 the field generated by the extrinsic (anomalous) currents concentrated in localized inhomogeneous domains. The fact that the integral formula (12) contains the adjoint Green's function for the diffusion equation means that, as in the wave case, a field determined in the space by these integrals represents an assembly of electromagnetic "waves" moving towards the observation surface. However, in view of the specifics of the diffusion reflection tf O. To simplify our discussion we assume that the secondary field is generated by the delta-pulse currents switching on at zero time. This field propagates in different directions, and is observed on the surface of the earth. The idea of electromagnetic migration consists of reverse downward extrapolation of the observed field in reverse the Note that in reality the behavior of the sources of the secondary field is more the am the am "'ran, I that the am t h e t h e instead of the background conductivity allows us to keep t h e traditional imaging conditions, namely to construct the image as a function of the migration secondary field at zero time moment. is shown in the Appendix A that such an imaging condition is correct for the earth with slow horizontal variations of conductivity and plane wave excitation. the Appendix B the migration apparent resistivity concept is discussed. The migration procedure is used to focus the field at the anomalous structure, then the analytic formula for the two-layered model and the quasi-plane wave exitation is used for an arbitrary earth model and an arbitrary source to image t h e conductivity structure. t h e t h e first t h e t h e far field consists of the part propagating through the air and the part propagating through a conductive earth. The air field can be approximated by the plane wave because it arrives at every point simultaneously, and the part propagating through t h e ground can be neglected because it attenuates in the conductive e a r t h (see Figure l a .) use an analogy with seismic migration. As it is shown by Claerbout (1972), suming up records from many receivers (homogeneusly distributed along the profile) one can obtain a single record corresponding to a plane wave incident field. The primary field propagation for this case is shown in Figure l b. t h e field t h e surface by the buried electrical cable could be computed using the Green's function of the diffusion equation C, which is given by formula (7). The field of an electrical cable on the surface is presented in Figure 2. The background resistivity is 100 ftm and t h e cable depth is 100 m. 13 migration conductivity O'm in such a way that the extremum of the migration field at zero time coincides with the geoelectrical boundaries. Using the migration conductiv­ity the It horiwntal In the Approximation of the source field by the plane wave is valid for two cases. The first case is the far-field of the fixed transmitter (Zhdanov and Keller, 1994). The source's the earth lao ) The second case is a slingram mode survey design. For the slingram mode we can field. field lb. To illustrate the focusing properties of the migration, we migrate the field gener­ated by an infinitely long conductive cable. The electrical field excited on the surface G, Om the The result of migration is shown in Figure 3. It indicates that the migrated field 14 S o u r c e A i r ( i n s u l a t o r) P r i m a r y f i e l d p r o p a g a t i on O b s e r v a t i o n profile a ) E a r t h ' s surface C o n d u c t i v e e a r th A i r ( i n s u l a t o r) A n o m a l o u s b o dy S o u r c e s E a r t h ' s s u r f a ce C o n d u c t i v e e a r th A n o m a l o u s b o dy wave, zone, Air (insulator) Source Primary field propagation Observation profile - - - -~ - - - - - -::::--- - - - - - - ~ - - -,- - - -===- r -.- - _._, I I Earth's surface i I ,i j I I : i i I I j I i j V i , I Conducti ve earth ~ VI - Anolllalous body a) Air (insulator) Sources Earth's surf-ace : i i I ! I b) I i , I -.!t ! I i i i ~ Conductive earth I I VI - I VI Anomalous body Figure 1: Approximation of the primary field propagation by the plane wave. a) Primary field propagation in the far zone. b) Primary field in slingram mode. ]_ : : 2 - illlHIPWill^Bllillll^lllllMlliiHI^MiilllMHHIIHHIilllllll 0.0 0.0 0.0 1.2 13.3 42.6 931.7 Amplitude 1 5 10 20 ~ 50 g 100 t=1§: 200 500 1000 2000 5000 10000 -1000 -800 -600 -400 -200 0 200 400 Distance (tn) Field Atnplitude 600 800 1000 931.7 Figure 2: Surface response of an electrical field simulated for a buried cable. 15 -1000 -800 -600 -400 -200 0 200 400 600 800 1000 Distance (m) 40.0 70.0 100.0 130.0 160.0 190.0 220.0 250.0 280.0 t h e electrical t h e t h e 16 0 50 >0O :§: -s >50 !!; 200 250 - o O iSluncc (m) I Field Amplitude Figure 3: Migrated field of the elect rical cable response. Cross shows the location of the cable. migrated point at the source location, whereas the migration field has a maximum at t h e source location. In the real case where the subsurface inhomogeneous structure is excited by the field source located on the surface, the migration field still has a maximum at t h e point of an anomalous structure, which acts like a secondary source. fields t h e was shown in (Zhdanov and Booker, 1993; Zhdanov and Keller, 1994) that in t h e frequency domain for plane wave secondary (upgoing) and primary (downgoing) fields everywhere inside t h e conducting layer are characterized by different amplitudes and phases. On the geoelectric boundaries their phases coincide (or shifted by while the amplitude of the upgoing wave is proportional to the amplitude of the downgoing wave and the proportionality factor is equal to the reflectivity coefficient j3. When passing into the time domain the above regularities manifest themselves in t h e fact that the time pulses of the upgoing and downgoing waves differ everywhere inside t h e conducting layer and mutually proportional to t h e coefficient (3 only at the geoelectric boundary (see Appendix A). We can introduce a time domain apparent reflectivity function, (3ta(r,t), defined as the ratio of t h e secondary E*(r,t) and primary Ep(r,t) fields, as the following: t3ta(r,t) = Es y(r,t)/E*(r,t). (19) downward extrapolating upgoing and downgoing fields in t h e lower half space. How­ever, it can be generalized to the results of time-domain electromagnetic migration. In Appendix A we show that the migration apparent reflectivity function at zero time, fi™r), can be determined also from the values of t h e migrated secondary field, calculated at zero time and normalized by some function D(f), 17 has a maximum at the point of the cable location. This illustrates how the migrated field represents the location of a subsurface target. The actual field has a singular the the Furthermore, we develop the technique to transform the EM migration fields and their different components into resistivity images of the vertical cross-section. It the the amplitudes 7r), up going coefficient (3. the everywhere the the f3ta(i', t), defined the E;(i', t) E~(i', t) fields, The determination of this function, obviously, depends upon the procedure for the 13;: (i'), the field, D(i'), « * ( ? ) Er(P,0)/D(f), (20) D(f) t h e <p(z,P) f+oo f) = x, = / H?*(x, t)tf (t) dt, Jo (p(z,t) = a- exp T equation: 27/27T5/2 and (21) (22) (23) fi™ reflectivity j3 the fi. (24) on proportional to the reflectivity of the boundary and to the magnitude of t h e downgoing (primary) field at t h e same point. The reflectivity coefficient is connected with the conductivity contrast ACT at the geoelectrical boundary by the simple formula ( see Appendix A): - y/(Jn + ACT ro\l\- yj(jn + ACT t h e determine t h e resistivity below t h e (within the anomalous domain or layer): P = i-fi 18 1I:;'(i') = E;"'(f, 0)/ D(i'), (20) where D(i) is the convolution of the migrated primary field at the same depth and the analytical function tp(z , I) ( +ro D(i') D(x, z) 10 ll';"( x,z, l)<p(z, l)dl, (21) where (22) In the last equation: (23) The function {1': is equal to the actual rcOec.tivity coefficient (3 exactly at the position of t he boundary: 11:;'( i') = II (24) This result is based 011 the fact that the migrated secondary (upgoing) field has a local extremum exactly at the geoelectric boundaries. This extremum is proportionaJ the (he tJ.a 11- ,;a;; - Vad 6.a - ,;a;; + van + 6.a· From the last formula we can determi[]e the resist ivity beJow the geoelectrical bound­ary [ 1+ II]' p= 1-11 po, (25) pn t h e t h e migration resistivity. t h e resistivity pm t h e t h e pm function We can compare this situation with t h e time domain field, where the depth of the field penetration inside the earth is proportional to the square root of t h e time after t h e current pulse in the transmitter: for early times the field is concentrated in the near surface layers; for later times t h e field penetrates deeper in the earth. In the case of the migrated field at early times this field is in the vicinity of its sources - the geoelectrical boundaries. At later times this field propagates upward to the e a r t h ' s surface. So at early time the migrated field behavior is determined only by t h e properties of the medium near t h e boundary between two layers. Therefore it is plausible that eq. (19) could be used for determining the conductivity contrast between these layers. determine the resistivity of the second layer, then find the conductivity contrast between t h e second and t h i rd layer and determine t h e resistivity of t h e t h i rd layer, etc. This layer by layer process opens the way for direct imaging of complex geoelectrical structures. In a model with a complex multilayered structure we can take as pn an apparent resistivity of the earth pa{t) obtained by averaging along the observation profile. The formula connecting the t i m e with depth d for the apparent resistivity calculation is provided by the following approximate expression for the skin-depth (Zhdanov and Keller, 1994): 19 where Pn is the background (normal) resistivity. The last expression gives the exact value of the resistivity of the second layer in a two-layered model (see Appendix A). However, we can use this expression for an arbitrary model to define the migration apparent resistivity. So, the migration apparent resistivi ty pm can be calculated using the formula (40) at any point of the profile. Since the migration apparent resistivity pm is a function of depth, we obtain immediately the geoelectrical cross section. The migrated secondary field in expression (48) (Appendix A) is calculated at zero time, which corresponds to the secondary field in the vicinity of the source. the the after the the earth's the the In the case of a multilayered cross section we can start from the first layer and the third the the third Pn Pa (t) time t (26) Thus an inhomogeneous profile of the earth's resistivity can be approximated by an homogeneous model which is used as the background (normal) resistivity for the migration. We call this normal resistivity profile the mean resistivity model. velocity distribution is replaced by the mean velocity (Claerbout, 1985) Also, as in seismic prospecting, we can use recursive migration algorithms which are based on successive downward extrapolation from level to level with different background resistivities. The algorithm described above has been coded. Sample computations using this code are presented in the next section. 2 . 3 N u m e r i c a l r e a l i z a t i o n o f t h e m i g r a t i on a l g o r i t h m For numerical realization of the EM migration algorithm, equations (15) and (17) must be reformulated in a discrete sense. These equations are essentially convolutions of the field observed on the surface with the Green's function or a derivative of the Green's function. Implementing these convolutions one can expect certain problems connected with the Green's function behavior. These problems are related to the fact that depending on the background conductivity value and migration depth, the Green's function and especially its time derivative can vary much faster than the observed field. Besides, in practice electromagnetic data are usually represented in logarithmic time, and the observed field is approximated by t h e average value within each specific interval or time gate. So, rewriting time and spatial integrals in (15) over t h e whole integration domain as a sum of integrals over subdomains (time gates At(i) and x-intervals Ax(k) we have: d_ dV lJ'0o~mz' 47T ,/ Nx Nt 1 H™(x',z\T-t') BZJx&t) (t t>) 20 This approach is similar to the one used in seismic prospecting where the real resistivi ties. 2.3 Numerical realization of the migration algorithm convolutions problems field. the time gate. the domain sub domains llt(i) llx(k) ) -aa Hxmz (x" ,z,T-t' ) = f.lofJmZ' ~L L~ 1 1 -aa Hx0 z ( X,O,t ) -'-(---1--- :-)-2 X t' , 47r i=l k=l 6.t(i) 6.x(k) t, - t' x exp /i0ar 4 (t t>) (a;' - a:)2 -f (z' - z)' dxdt. t e rm derivative, that H® z(x,0,is a step-response, to get f) .. „ yi NX Nt 4 7 F fcl k=\ J A t ^ JAJ*x*((kk)) d x (j,0cr7 d t { ( t - t > y e x p \ 4{t-t>) (x' - x)2 [z - z)' dxdt. Assuming that H®z(x,0,t) is constant within each time gate At(i) and each x-interval Ax(k) t h e h e t h e derivative a discrete form of the continuous convolution integral (15) as \Hl,z{Xi^^i) G™{x',Xk,Xk+1,tk,tk+l,t',Z ~ Z)) , Nx Nt E E i=i k=i where xk+i and xk are x-interval boundaries Ax(k) x k + 1 - Xi, and ti+i and t{ are time gate boundaries At(i) = U+i - U The function G™(xkl Xk+i, t k , t^+i, - migration kernel G^(x\xk}xk+Utk,tk+Ut\z' - z) = - m X 47T x exp { _~oum [(x' - x)' + (z' - z)'] } dxdl. 4(1 - 1') 21 We can express the second term as a time deri vati ve, assuming t hat H2,z (x, O, t) a m I I r I J10amz I N~ N r ;, 0 al ,H.,(x ,Z ,T -I) = 4 I: I: [H.,(x, 0, t) x , 7r ;=lk=1 llt(i) JlJ. x(k) , a (I r al (I _ I')' exp 1. - 4 (/LI o_u mI' ) [(x ' - x )' + (.z , - z) ']}) dxdl. H2,zCx, O,t) .6.t(i) x­interval ~x( k ) we can take it out of the integral and after that calculate t he integral over t of the time deriva.tive of the Green's function. As a result we obtain a. di screte cont inuolls convolut ion ·a-aI' HZm,'z ( x, ,z, , T - I') = Ncr Nr 2: L: (H~,z(Xj, 0, td X G~I(XI, Xk, XHI , tk, tk+l, t', z' - z)) , (27) ;=1 k=l XHI Xk l>x(k) = X'+1 - X;, lit1 ti l>1 (i) 1;+. I;. T he fundion G;;(x', Xk, Xk+ I , k , lHI! t', z' - z) is called the migmtion kernel x exp 4 ( t , + i - t ' ) (x' - xf + (z - z)' dx. (28) kernel (28) in the formula (27) by its time integral, (which could be calculated ana­lytically) G 5 ( * ' , **, xk+1, tk , t k + 1 , t\ z'-z) = £ + 1 - _ - _ - X x ( e x p { - 4 ( C T [ ^ - ^ ) 2 + ^ - ^ ) 2 ] } - )<fe. precomputes values of t h e migration kernel G™xk, xk+i,tk, tk+i,t', z') for all sampling intervals Ax(k) and At(i) at particular depths z'. These values are stored as a matrix G™Sampling over time and space should be chosen according to survey parameters. Then t h e migration program convolves matrices G™with t h e matrix of magnetic field values H® (i,k), thus obtaining a migrated field at the depth z'. Repeating this procedure over all values z' one obtains a migrated field at each depth. 2 . 4 T e s t i n g o f t h e a l g o r i t hm Ey 22 X rXk +1 1 {/-lOam [( I )2 (' )2]} JXk (ii+l- t')2exp -4(ii+l- tl ) X -X z -z - 1 {/-lOam [( I )2 (' )2]} - (ti_tl)2exp -4(ti- t' ) X -X + z -z dx. To obtain the formula for the delta-pulse excitation it is enough to replace the X (exp {- (/-lOam ') [(x' _ X)2 + (z' _ Z)2]}_ 4 ti+l - t - exp { - 4 ~~~mil) [(x' - X)2 + (Z' - Z)2]} )dx. (29) The migration computer code is based on formulas (27) and (28). The code the G;;o (Xk' Xk+l, tk, tk+l, t', Zl) 6.x( k) 6.t( i) Z'. G;;o (i, k) . the G;;o (i, k) the H~z,( i, k), Z'. Zl 2.4 Testing of the algorithm To check the accuracy of the migration algorithm we compared the secondary field in a two-layer earth model excited by a plane wave, migrated analytically and numerically. The analytical solution for the migration problem for a plane wave in a layer covering a homogeneous half-space is in Appendix A. The Ey component of the secondary electric field for the model with one layer of a\ <J2, Q E a r t h ' s E;(z = 0,t) = Q/3a^exp(-2T2(yf), (30) a r expressions 27/2^5/2 a T = 27Tyj2t/({10-1), t h e formula moment £ ' w - 0 ) - ^ + ^ ' (31) <Ti 5 / m , o-2 m, am = 0.01333 m. $ = 0 . 0 5 V/m. z=0 A£ Ifis) used. We can see t h a t the migration kernel changes with time much faster than the secondary field. This is t h e reason why we cannot directly discretize formula (15) for its numerical realization, applying instead the formula (27). Formula (27) uses the discrete migration kernel (27), which is an average value of t h e Green's function over t h e computed 23 thickness h and conductivity 0"1 covering an infinite half-space with conductivity 0"2, excited by a plane wave with amplitude Q is expressed on the Earth's surface at any time moment t using the following formulas: where variables and T are given by the expressions a=--- and (3 is the reflectivity coefficient, which is given by the formula fJ = 0il- 0i2. 0i1 + 0i2 The migrated electrical field amplitude at any depth and at zero time moment can be expressed as (31) The parameters for the test are 0"1 = 0.01 S/m, 0"2 =0.02 S/m, and O"m 0.01333 S/m. The first layer thickness is 100 m, and Q=0.05 Vim. The secondary field, given by formula (30) and the migration kernel, both at z=O are presented in Figure 4. Here the constant time discretization (t::.t = IllS) was that t~an the the each particular time gate. The comparison of the analytical formula (31) with the migrated field computed 60 Time (|J,sec) the 24 W t------=~------~------~--------~------T ' 0 - 40 ·..'. . . . ~ :· .. ~<• 30 ····· ···· 20 ····· , ' 0 , O.""<;e(:) Figure 4: Secondary field for t.he model with one horizontal layer (solid line) and migration kernel (dashed line) 30 O O-I SO1 lOO1 Dep1th5,01 ( in) 2001 2SO, 300 field. 1000 m (dashed line) and 500 m (dots). algo­r i t hm numerical realization the time gates were chosen with 12 points per decade. Spatial sampling interval was equal to 40 m. The study has shown that the algorithm is sensitive to the spatial limits (aperture) of integration. To produce a reliable picture t h e limits of integration must be approximately 10 times wider t h a n the depth of the target. grid spacings of 40 and 100 m. We can see that error decreases with decreasing spacing. An empirical conclusion from this study is t h a t spatial sampling interval must be less t h a n half of t h e depth of the target to avoid significant integration errors. E -< 3°T-------~----~------~------~----~------T 20 , " . o+-------~----~------__ ------~----__ ------+ o 50 100 1SU Depl.lI .. (In.) 200 250 300 25 Figure 5: Analytic solution (solid line) and numerically migrated field. Aperture is numerically for different apertures (see Figure 5) shows that the error in the algo­rithm is small and significantly decreases with increasing integration limits. For the the than Figure 6 compares the migration results obtained for an aperture of 2000 m and spacmg. that than the 26 26 ,0 r---- - - ----- -------, ,0 .-- Figure 6: Analytic solution (solid line) and numerically migrated field. Grid spacing is 40 m (dashed line) and 100 m (dots). 3 . M O D E L S T U DY 3 . 1 M o d e l i n g s t r a t e g y a n d m e t h o d t h e t h e bodies; and 2) models with local resistive and conductive bodies. These models simulate anomalous zones which are targets for geophysical surveys. t r a n s m i t t e r offset mode with fixed transmitter - receiver separation. Time domain electromagnetic data are obtained as a result of forward model­ing using computer code developed by I. Adhidjaja, G.W. Hohmann, and M.L. Oristaglio (1985). This code is based on t h e finite-difference method which gives a forward time domain solution for a 2D body in a conductive earth. This code is a reformulation of t h e TEM modeling method of Oristaglio and Hohmann (1984) in terms of t h e secondary field. The EM field in t h e models is generated by an infinitely long cable. The observed field is dHz/dt and t h e transmitter-receiver separation (off­set) is equal to 20 m. The cross-sectional dimensions of all anomalous bodies in t he models (representing local inhomogeneities) are 20 x 20 m 2 . and stability for imaging different geophysical targets. We will demonstrate t h a t TDEM migration is a fast and stable method of geoelectrical imaging, satisfac­torily resolves the local inhomogeneities and gives a good estimate of the location, t h e t h e t h e t h e 3. MODEL STUDY 3.1 Modeling strategy and method In this section the effectiveness of the TDEM migration method is analyzed for dif­ferent geoelectrical models that can be considered typical for geophysical exploration and environmental studies. The set of models under investigation includes: 1) models with local conductive The theoretical survey in these models was conducted in the transmitter offset J.1. the the the the infinitely dHz/dt the the 2 • The results of numerical modeling and TDEM migration demonstrate the effec­tiveness demonstrate that dimension and resistivity of the anomaly. The method, however, requires a good initial estimation of the background distribution of the conductivity. It is also based on the availability of the spatially dense measurements along the profile (2D) or over 3 . 2 M o d e l w i t h o n e l o c a l c o n d u c t i v e b o d y t h e - m is embedded in an homogeneous half-space - t h e t h e modelling transmitter/separation, where t h e transmitter and receiver were moving together along t h e profile. is assumed that the transmitter is excited by a step pulse and the receiver measures t h e time derivative of the magnetic field. t h e procedure performed with a known background resistivity using formula (27) produces the time derivative of the migrated magnetic field, which is shown in Figure 10. Applying the transformation which is given by formula (16) one obtains t h e migrated electrical field, which is shown in Figure 11. This field indicates the depth of t h e anomalous body, as follows from the extremal properties of the migrated field. The estimated reflectivity function is shown in Figure 12. Using relationship (40), Appendix A, one finally arrives at t h e migration apparent resistivity cross-section (Figure 13), which could be viewed as a final result of the migration procedure. One can see t h a t the migration apparent resistivity indicates the position and the resistivity of t h e conductive body. 3 . 3 T h e s e p a r a t i o n o f p r i m a r y a n d s e c o n d a ry field algo­r i t hm field associated field. 28 a surface area (3D). 3.2 Model with one local conductive body The steps of the migration can be illustrated using the following model. One local conductive body with resistivity 10 Ohm - Tn is embedded in an homogeneous half­space with resistivity 100 Ohm - m. The vertical cross-section of the body is 20 x 20 m, while the depth to the body's top is 100 m (Figure 7). The numerical modelling code simulated a survey in slingram mode with a fixed 20 m transmitter/receiver the the It the The modeling program produces a primary field which does not change along the profile (Figure 8). The time derivative of the secondary field, however, clearly indicates the spatial location of the anomaly (Figure 9). The migration procedure time field, the field, the field. reflectivity finally the section that migration the 3.3 The separation of primary and secondary field An important problem arising in the practical realization of the migration algo­rithm is the separation of the secondary (scattered) field, i.e., the field associated with targets (scatterers), from the primary field. EM receivers measure the total field, which could be represented as a sum of the primary field (the field which is Figure 7: Local conductive body in a homogeneous media. 0 20 40 60 8' 80 .s 100 "'"- <l.> Q 120 140 160 180 200 200 400 600 20 800 1000 Distance (In) 100 OhIn-In 1200 1400 1600 1800 100 29 8000 10000 200 400 600 800 1000 1200 Distance (m) 1400 1600 1800 •mmm i 3.0e-05 4.0e-05 6.0e-05 8.0e-05 1.2e-04 1.8e-04 3.3e-04 7.4e-04 3.0e-03 5.8e+00 Amplitude ~ ~ ~ 0 2000 4000 6000 800 1000 1200 Distance (Ill) 3.0e-05 4.0e-05 6.0e-05 8.0e-05 1.2e-04 1.8e-04 3.3e-04 7.4e-04 3.0e-03 5.8e+00 Alllplitude Figure 8: Time derivative of primary field for the model of Figure 7. 30 31 Fi gure ~ 5 ~ 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 200 400 600 800 1000 1200 1400 1600 1800 Distance (m) -2e+00 -3e-04 -1e-04 -1e-04 -Oe+OO Oe+OO 1e-04 1e-04 3e-04 2e+00 AIllplitude Figure 9: Time derivative of secondary field for the model of Figure 7. t h e 32 0 20 40 60 80 100 g 120 ..s 140 80. . 160 180 200 220 240 260 280 200 400 600 800 lOOO 1200 1400 1600 1800 Distance (m) -1 e-02 -1 e-03 -3e-04 -le-04 -Oe+OO Oe+OO le-04 3e-04 le-03 le-02 Amplitude Figure 10: Migrated time derivative of magnetic field for the model of Figure 7. -3e-06 -7e-07 -2e-07 -le-07 -5e-08 -3e-08 -le-08 -7e-09 -4e-09 -9e-ll Amplitude :g -= ~ 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 200 400 600 800 1000 1200 1400 1600 1800 Distance (Ill) ge-AIllplitude Figure 11: Migrated electric field for the model of Figure 7. 33 _4e-01 -6e-02 -2e-02 -le-02 -6e-03 -3e-03 -2e-03 -8e-04 -3e-04 -Oe+00 -Oe+00 Amplitude g -S c.. <U 0 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 200 400 600 800 1000 1200 1400 1600 1800 Distance (m) -OO OO Figure 12: Apparent reflectivity function for the model of Figure 7. 34 200 400 600 800 1000 1200 1400 1600 1800 m) 17 20 25 30 37 45 Ohm-m 55 67 82 100 Fi gure 35 0 20 40 60 S 80 i 100 0 120 140 160 180 200 1000 1200 1400 1600 1800 Distance (rn) 37 45 Olun-m Figure 13: Migration apparent resistivity for the model of Figure 7. 36 of slingram mode survey design and homogeneous media the primary and secondary field can be separated easily. We can see in Figure 8 that the primary field does not change along t h e profile. Figure 9 shows that the secondary field varies in both space and time. A plot of t h e secondary field values at 0.2 jisec versus X is shown in Figure 14. We can see from the plot that the area bounded by positive portions of t h e curve and the X-axis is equal to that bounded by the negative curve segments and the X axis. In Figure 14 we can see that the curve has mirror symmetry with respect to t h e point right above the anomalous body (X=0). We can understand the behavior of the primary and secondary fields by analyzing the model presented in Figure 15. For a constant separation between transmitter and receiver, t h e primary field does not change as the system moves along the profile, since the geometric relationship between the source, t h e receiver and the horizontal layers do not change with respect However, t h e secondary field changes along t h e profile. Let us assume, for example, that the separation between the transmitter and receiver is small in comparison with t h e target depth. In this case, considering secondary fields we can assume the offset is zero. Now, consider two points Xi and 2 , symmetrically located about the center of the body. The transmitter (in this case the infinite cable) located at the point X\ produces an anomalous current in t h e body, with exactly the same direction and amplitude as if we place it at the symmetrical point X 2. The schematic vector lines of magnetic induction vector B are shown as dashed lines in Figure 15. The receiver (a small horizontal loop in this case) measures t h e inductive current which is proportional to the vertical component of magnetic field. We can see t h a t vector lines at point X\ are oriented upward, whereas at the point X 2 they are oriented downward. This means that the Bz components at the points X\ and X 2 have opposite signs. But, because the geometry is symmetric, they have exactly the same amplitudes. Looking at Figure 15 we can conclude that the spatial spectrum of the anomalous magnetic field does not contain an harmonic with zero spatial frequency. In contrast, the spatial spectrum of the primary field contains only the zero frequency component. observed in the absence of inhomogeneities), and the secondary field. In the case the the flsec the the O). the the to the source-receiver pair position. the the the Xl X 2 , Xl the 2 • the proportional field. that Xl 2 Bz Xl 2 component. I 1 1 1 1 I 1 1 1 200 400 600 800 1000 1200 1400 1600 1800 Distance, m 14: 0.2fis. 0 . ' i., 0 .0 -0. 1 200 400 /\ 600 800 V ' 000 Ois tance. m ' 400 Figure 14: Secondary field at O.2JlS. 37 '600 '800 Transmitter and receiver pair V V V Layer 1 Layer 2 Layer a ) Magnetic vector anomalous body Magnetic vector ~- ^ Magnetic field vector - - line _ - - ' b ) ~ Figure 15: Primary/secondary field separation, a) The primary field does not change horizontally media, directed differently at symmetric points. a) 38 Trans",i!l"r I Lay c .. 3 Trans mitter / _____ - - - - - - _ _ _ _ _ Trans millcr R:~~vh "",,"" -------___ ' n R~~~V"' < B.'lrth s urfac e , b) ,, , , -, anom a lo u s h<xJy D Mag netic field v ect o r ____ ~inc ____ _ " ,,, " M"g n ctic vecto r secondary separation. along the profile above the hori zontally layered medi a, b) secondary field is directed different ly 39 t h e dBgj^ dB'(t) = 1 ^dBt(t,Xi) dt Nx2-! dt X{ formula dt ~ dt dt • (Si> gives complex, application of spatial filtering for primary/secondary field separation has produced a reasonable result for the field data. 3 . 4 M i g r a t i o n o f n o i s y d a ta t h e method, dHz/dt E a r t h ' s t h e between - 5% and + 5% of secondary field maximum. model without adding noise is shown in Figure 13. As we can see, comparing the two pictures, t h e migration apparent resistivity calculated from the noisy d a t a is very close to t h a t calculated from the noiseless data. Similar results are obtained from the migration of noisy data for other models, discussed in the following sections. high-frequency random noise. This type of noise is typical for geophysical data. However, it does Thus, using spatial filtering, we can separate primary and secondary fields. This filtering can be realized by calculating the average value of the total field B~~P) along every time interval aBf(t) = _1_ ~ aB!(t, Xi) at Nx i=1 at ' BBP(t) where at represents the primary field. The secondary field can be obtained for every point Xi using the formula aB;(t, Xi) at aB!(t, Xi) at aBf( t, Xi) at (32) Formula (32) gIves an exact solution for horizontally-layered media containing local inhomogeneous bodies. In spite of the fact that the real situation is more 3.4 Migration of noisy data In order to examine the effect of noise on the resolution of the migration method, we have added 5% random noise to the secondary field dHz/dt computed on the Earth's surface. Figure 16 shows the same secondary field as in Figure 9, but with random noise added. The amplitude of the noise was distributed with constant prob­ability maximum. Figure 17 represents the apparent resistivity image produced after the migration of this noisy secondary field. The migration resistivity image produced for the same l:t the data that This model illustrates the stability of the migration procedure for high-frequency field's maximum. ~ ~ ." .§ E-< 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 600 800 1000 Distance (m) 1200 1400 i· .. ~;m"'''''f''~''%''m,',;"" 'I· i -9.860 -00418 -0.296 -0.177 -0.056 0.056 0.175 0.294 00419 9.579 Amplitude Figure 16: Secondary field with 5% random noise to field's maximum. 40 200 400 600 800 1000 1200 1400 1600 1800 Distance (m) 17 20 25 30 37 45 55 67 82 100 Ohm-m data. 0 20 40 60 § 80 -S 100 e- O 120 140 160 180 200 0 17 200 400 600 20 25 30 800 1000 1200 1400 1600 1800 Distance (Ill) Ohlll-rn Figure 17: Migration image produced from noisy data. 41 42 3 . 5 R e s o l u t i o n o f m i g r a t i o n m e t h od t h e Ohm-and separation 100 m, embedded in an homogeneous medium of resistivity 100 Ohm-m at a depth of 100 m. Figure 18 shows the geoelectrical cross section of this model. Figure 19 shows that these two conductive bodies produce one single spot on the migration apparent resistivity image. Thus, for this separation, migration does not allow us to distinguish the two bodies separately. The next model is similar to the previous one with the difference that the separa­tion between the two conductive bodies is now increased to 200 m (Figure 20). The migration apparent resistivity image on Figure 21 shows that the two bodies can now be distinguished. demonstrates that for this separation the two bodies are very clearly distinguished and migration gives a good estimation of their location, resistivity and dimensions. the distance between them should be at least two times greater than their depth. 3 . 6 R e s o l u t i o n o f c o n d u c t i v e a n d r e s i s t i v e b o d i es is known t h a t a resistive body generally produces a much smaller anomaly than a conductive one. In order to demonstrate the capability of migration to distinguish a resistive body in t h e presence of a conductive body, we consider a model consisting of one conductive (20 Ohm-m) and one resistive (1000 Ohm-m) body embedded in a not truly represent geological noise, produced for example, by near-surface inhomo­geneities. 3.5 Resolution of migration method In order to investigate the capability of migration method to resolve targets located close to each other, we compute a set of models with two conductive bodies, gradually changing the distance between the bodies. The starting model consists of two conductive bodies with resistivity 20 Ohm-m Ohm­m distinguished. In the next model the separation between the two conductive bodies is increased to 300 m, as Figure 22 shows. The migration apparent resistivity image on Figure 23 distinguished The analysis of these models suggests that for resolution of two conductive struc­tures 3.6 Resolution of conductive and resistive bodies It that distinguish the consisting Q I I O 20 40 60 80 100 120 140 160 180 200 200 400 600 800 1000 1200 1400 Distance (m) 20 100 100 Ohm-m 0 S -5 fr Q 0 200 400 600 800 Distance (m) Ohm-m 1000 1200 1400 Figure 18: Conductive bodies with 100 m separation. 43 3 5 7 10 15 22 32 46 68 100 Ohm*m separation. 44 0 20 40 60 § 80 .,s 100 c.. ~ 120 140 160 180 200 0 200 400 600 800 1000 1200 1400 Distance (rn) 3 Ohrn*rn Figure 19: Migration image for conductive bodies with 100 m separation. 600 800 Distance (m) 1000 1200 140O 20 100 Ohm-m m separation. 0 20 40 60 g 80 .,s 100 c.. ~ 120 140 160 180 200 0 200 400 Distance (Ill) Ohm-Ill 1400 Figure 20: Conductive bodies with 200 ill separation. 45 100 46 Ohm-m separation. 0 20 40 60 E- 80 .,s 100 0.- ~ 120 140 160 180 200 0 200 400 600 800 1000 1200 1400 Distance (rn) 12 14 17 20 24 28 34 41 49 58 70 84 100 Olun-rn Figure 21: Migration image for conductive bodies with 200 m separation. 200 0 200 400 600 800 1000 1200 1400 Distance (m) 20 100 100 Ohm-m 0 20 40 60 :§: 80 .,s 100 fr ~ 120 140 160 180 0 200 400 600 800 Distance (Ill) Ohm-m 1000 1200 1400 Figure 22: Conductive bodies with 300 m separation. 47 48 separation. - ~, ... ,I -:: _; '\. " ;{ I 0 20 40 60 :§: 80 .,s 100 e- O 120 140 160 180 200 0 200 400 600 800 1000 1200 1400 Distance (rn) 15 18 23 28 35 43 53 66 81 100 Obrn.-rn Figure 23: Migration image for conductive bodies with 300 m separation. medium of resistivity t h e separation t h e t h e estimation of the bodies' location, dimensions and resistivity. Of course it is more difficult to see t h e resistive body on t h e migration image, because t h e conductive body is "illuminated" much more brightly. The model presented in Figure 26 is similar to t h e previous one, but now the resistive body is located at a depth of 50m from the surface, whereas the conductive body is located at 110 m under the surface. The migration apparent resistivity image in Figure 27 gives again a good estimation of t h e bodies' location, dimensions and resistivity. Thus, migration can be used not only for detecting a highly conductive structure but a resistive structure as well. 3 . 7 B a c k g r o u n d r e s i s t i v i t y i n f l u e n c e o n t he E M m i g r a t i o n i m a ge t h e t h e t h e background t h e resistivity. However, in real cases it is not usually known well and can only be approximately estimated. In order to examine the dependence of the migration results on t h e choice of background resistivity, we have applied the migration with a different background resistivity. The model used in this test is t h a t of one conductive and one resistive body (Figure 24). The migration apparent resistivity image ob­tained for this model using the "true" resistivity (100 Ohm-m) is shown in Figure 25. The migration images for the same model computed with the background of 80 and 300 Ohm-m are shown in Figure 28 and 29 respectively. In these cases, t h e position of t h e two bodies is shifted upwards or downwards. This effect is very similar to the one known in seismic migration. This is why velocity analysis is an important part of seismic d a t a processing before migration. In EM migration, it is also important to correctly determine the background resistivity in order to obtain best results. One approach to this problem is discussed in the paper by Traynin (1995). 49 homogeneous mediumofresistivity 100 Ohm-m at the depth of 100 m. The separation between the two bodies is 300 m (the resolution distance for two bodies as shown by the previous set of models). The geoelectrical cross section of the model is shown in Figure 24. The migration apparent resistivity image in Figure 25 gives a good the the the the the 3.7 Background resistivity influence on the EM migration image One of the requirements of the migration method is a knowledge of the background resistivities. Until now, we have used for the migration of our models the true back­ground the that (loa the the data 50 o 20 L ,1 1 _ J _ I, „ I „_ 40 - 60 - 80 100 120 140 - 160 180 - 200 200 400 600 800 Distance (m) 1000 1200 1 1400 [ ' 20 100 1000 1000 separation. 0 g -5 0.. ~ 0 1200 1400 Distance (m) 1000 Ohm-m Figure 24: Resistive and conductive bodies with 300 m separation. Figure 25: Migration image for resistive and conductive bodies with 300 m separation. 51 0 20 40 60 § 80 oS 100 ~ 120 140 160 180 200 0 200 400 600 800 1000 1200 1400 Distance (lTI) 17 21 27 34 42 53 67 85 106 134 OhlTI-lTI 20 - 40 - 60 80 - -fiiiiiiiiiiiit^ IOO - ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ p ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ H - 140 - 160 - - 180 - 2oo I i | | i 1000 1200 1400 Distance (m) | | •mmmm^^mmmammmamammmmmmm^ 20 100 1000 1000 Ohm-m 52 0 20 40 60 § 80 -5 100 0.. 8 120 140 160 180 200 0 200 400 600 800 Distance (m) 20 100 1000 1000 Olnn-Ill Figure 26: Conductive and resistive body located at different depths. 0 20 40 60 § 80 t 100 0 120 140 160 180 200 0 200 400 600 800 1000 1200 1400 Distance (ITl) ":::·:·:::·:::::::f.;:::::,r~ j 5. .". j 17 32 47 61 76 91 106 120 135 150 Olnn-ITl Figure 27: Migration image for model with conductive and resistive bodies located at different depths. 54 Ohm-m 0 20 40 60 E- 80 -= 100 e- O 120 140 160 180 200 0 200 400 600 800 1000 1200 1400 Distance (Ul) 10 13 17 22 29 38 50 65 85 111 OhlTl-1Tl Figure 28: Migration with incorrect background resistivity (80 Ohm-m). 55 0 20 40 60 § 80 -a 100 <U 0 120 140 160 180 200 0 200 400 600 800 1000 1200 1400 Distance (m) ,.,.,.,. .•., .. ....., .• . ,.. .. 13 64 115 166 217 268 319 370 421 472 Ohm-m Figure 29: Migration with incorrect background resistivity (300 Ohm-m). 4 . E M M I G R A T I O N I N W A S T E S I TE C H A R A C T E R I Z A T I O N to characterization existing DOE buried waste facility, the Cold Test Pit, where the location of targets was specified. Transmitting loops that were comparable in size to the targets were used for t h e TDEM survey. Because of this, ID interpretation fails to resolve separate objects, but it does provide a reasonable estimation of t h e background resistivity. The data, however, reflect the lateral location of the targets, thus indicating that more information could be potentially extracted. Since t h e n a t u r e of migration is essentially multidimensional, it was hoped that its application could better resolve t h e targets. Each profile was migrated separately. After that a composed 3D image of t h e Pit was constructed. Migration not only resolved separate objects laterally, but also provided reasonable vertical resolution and estimation of anomalous resistivity. 4 . 1 I N E L W a s t e C o m p l e x C h a r a c t e r i z a t i o n P r o j e ct of Energy (DOE) operational areas were investigated through a subtask of Technical Task Plan (TTP) AL911201, " Nonintrusive Sensing of Environmentally Sig­nificant Objects and Sites" (Chem-Nuclear Geotech 1992a). The subtask (AL911201- G2, "Three-Dimensional Site Characterization Using Broadband Electromagnetics") was developed to investigate various broadband electromagnetic (BBEM) systems t h a t demonstrate promising systems under actual field conditions. The BBEM methods investigated 4. EM MIGRATION IN WASTE SITE CHARACTERIZATION This section describes the application of TDEM migration Lo the characterization of a US DOE waste site. A large volume of TDEM data was gathered over an TesL PiL, loca.tion ta.rgets the lD the daLa, the nature essenLially the the 4.1 INEL Waste Complex Characterization Project Nonintrusive methods of characterizing waste sites and waste forms at U.S. De­partment Chern-Electromagnetics") that potentially could characterize DOE waste sites in a 3D sense and to demonstrate were intended not only to locate conductive objects and zones but also to locate accurately the lateral and vertical boundaries of the conductive zones, measure the thickness of any capping material, and provide a qualitative estimate of the type of DOE (BWID) 4 . 2 S u r v e y d e s i g n a n d d a t a d e s c r i p t i on oppor­t u n i t y simulating two important features, namely conductive objects located near the 0s (South) profile and a conductive zone near the 35s profile. Objects (stacked conductive drums) were buried at depths from 3-5 m. offset or slingram mode as described in MacLean (1993) and as illustrated in Figure 31. this mode of operation, the transmitter and receiver coils are offset by a fixed distance. The receiver is placed outside of t h e transmitter loop and is separated from t h e nearest side of t h e loop by a distance of 10 m. The reading point is taken at the midpoint of t h e two coils, because of the reciprocal relation between transmitter and receiver. Unlike FDEM systems, the coil separation is not critical to the depth of investi­gation because the primary or inducing current is off during the measuring period. For t h e same reason, the TDEM system is not sensitive to noise and signal variations caused by small changes in t h e system geometry. TDEM readings were taken over a broad range of time delays on these lines at 2 m intervals. The initial set of readings was recorded over the earliest available series of 20 time delay windows. This initial data set was reviewed in the field; additional d a t a 57 waste contained in a particular pit. A full-scale field demonstration of the PROTEM 47 system at an existing DOE waste site facility was proposed. The Buried Waste Integrated Demonstration (BWID) supported this demonstration at three waste sites in the Radioactive Waste Manage­ment Complex (RWMC) at the Idaho National Engineering Laboratory (INEL). 4.2 Survey design and data description The locations of objects in the Pit were known a priori, which provided an oppor­tunity to check interpretation results. The Pit contains conductive objects simulating waste. The location of the objects is shown in Figure 30. Overall, the Pit possesses Os BBEM survey traverses were conducted over a grid consisting of 5 m spaced traverse lines crossing the Pit. The grid was extended over a sufficiently large area to cover the region of interest. The survey was conducted in the transmitter offset 3l. In fixed the from the the the receIver. the the TD EM data sets at longer time delays were recorded if a measurable signal was observed. D i s t a n c e ( m e t e r s) 0 10 20 30 40 50 60 70 80 90 N o r t h I n f e r r e d p i t b o u n d a r y S u r v e y m a r k er 2 3 2 3 f 1 )• a I n f e r r e d p i t b o u n d a r y - 1 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 O o CD B CD «-K CD E a s t 1 ) C o n d u c t i v e o b j e c ts 2 ) E a r t h b e r m s 3 ) N o n c o n d u c t i v e o b j e c ts 58 Distance (meters) o 10 20 30 40 50 60 70 80 90 -10 -::::::::::::::::::::::::::::liiferred----i~:~:~~:~~:~~y::::::::::::::::: 0 __~ _':l_~y~y._~~~~~~_____ _ ___ _____________________________________ _ - - - - - -- - - - - -- -- - -- -- - - - - -- - ---1------:::---+ -- -- -- - - - -- -- - -- - - - - -- -- - - -- - - - - - -- - - -- 10 ------------------------------ ---------3-------- --------------------------------------- North ------------------------------ 20 U i-" • [J) f'""'t- ------------------------------ ---------3-------- --------------------------------------- ---------_._---------------------------.-.----------------_._--------------------_ ... _--- ~ 30 ::::s (j (l) ,,-... 40 a (l) --------------------------------------3----------------------------------------------- f'""'t- (l) _::::::::::::::::::::::::::::: ~#f~ij~~: p.i~ :~:~~:~~:~~y::::::::::::::::: 50 ~ [J) "-' 60 70 East 1) Conductive objects 2) Earth berms 3) Nonconductive objects Figure 30: Location of anomalous objects in Cold Test Pit. S turn 5m Grid spacing 2m St-' Pit. Transmitters, 8-tum Sm square loops 1.9 A current " " s .~ N .... ".~.: receIver '.- --- --/-g-----@------ -- --- -- ----- 59 . _____ @ ____ @_____ __ ________ ____ _ ~Ofile ~ spacing Sm ---@ -----@------------------------ . . 'E • Transmitter/receiver separation 12.5 m Figure 31: Slingram mode TDEM survey design at Cold Test Pi t. 60 collected data are represented in Figure 32 and Figure 33. We can see from these figures the clear anomaly within the pit borders. Figure 32 presents the total field while Figure 33 presents the secondary field, which was separated using t h e algorithm described in subsection 3.3. using in Figure 34, but no vertical resolution is achieved. The image starts below 4 m depth, because according to the ID assumption which is the basis of "smoke ring" concept, there is no information in the data above 4 m. The upper part of t h e image is distorted for the same reason. On the part of the profile without the anomaly, at a depth of 6 m the resistivity is about 100 Ohm-m. This value was taken as the background resistivity for the migration procedure. ID was achieved. d a t a sufficient quality to be migrated. 4 . 3 A p p l i c a t i o n o f m i g r a t i o n f o r d a t a i n t e r p r e t a t i on Figure 30): 15n, 0s, 5s, 10s, 15s, 20s, 25s, 30s, 35s, 40s, 45s, 50s, 60s. An effective background resistivity of pn = 100 Ohm - m was used. As a result of processing the TDEM data by t h e migration method a set of vertical cross sections in t h e Cold Test Pit was obtained. Cross sections of migration apparent resistivity along the profiles 0s and 35s (see Figure 30) are presented in Figure 35 and Figure 36. We can clearly see not only the lateral extent of t h e anomalous targets, but also the vertical extent of t h e bodies, which is better than the results obtained by ID interpretation. The profile images have been interpolated horizontally to obtain the 3D resistivity volume image of t h e Pit. As an example, the maps of resistivity at the depth 3 and 5 m are The BWID Cold Waste Pit was covered by the profiles located 5 m from each other. The grid for the Cold Test Pit survey is shown in Figure 30. Two samples of \Ve the algorithm The "smoke ring" image, obtained usmg the program written by Eaton and Hohmann, (1989) is presented in Figure 34. The conductive zone is outlined lat­erally 1D the 1D inversion for these data was also done by MacLean, (1993). No vertical reso­lution Having a known background resistivity, in spite of random noise, these data have 4.3 Application of migration for data interpretation Thirteen profiles were processed, numbered from the North to the South (see Os, lOs, effective pn - m was used. As a result of processing the the the Os the the 1D the III 10 20 '[ E,- o 30 E E::: 40 50 0 10 20 30 -990 -654 -76 118 40 50 60 70 80 Distance (Ill) ..... i .' ~-T' .Cw ... w .• "" i tim",to.: .. ::.· 166 220 287 370 482 672 995 AIllplitude Figure 32: Cold Test Pit TDEM data (profile Os). 61 -10 0 10 20 30 40 50 60 70 80 Distance (m) -988.1 -234.8 10.0 15.9 24.8 37.3 Amplitude 54.5 86.0 148.0 290.2 974.7 TDEM -10 o 10 30 Distance (m) 24.8 37.3 Amplitude Figure 33: Cold Test Pit TD EM data (profile 35s). 62 8 H & 10 12 14 10 20 30 40 50 Distance (m) 60 70 80 90 50 100 150 250 Ohm-m 300 425 500 600 4 6 § -c<=1.) . 10 0 0 o 40 50 Distance (m) 250 Obm-m Figure 34: "Smoke ring" image of Os profile. 0 2 4 § ~ 6 0 8 10 12 lO 20 30 40 50 60 70 Distance (Ul) 10 18 27 35 43 52 60 68 77 85 93 102 110 Ohm-m Figure 35: Migration apparent resistivity for Cold Test Pit (profile aS). 64 Distance (m) 10 18 27 35 43 52 60 Ohm-m 68 77 85 93 102 110 0 2 4 S ~ 6 0 8 10 12 10 20 30 40 50 60 70 52 60 OllIn-m Figure 36: Migration apparent resistivity for Cold Test Pit (profile 35s). 65 p i t ' s sections filled with the drums and boxes. Note that elongation of t h e conductive anomalous zone in the northern section of the pit outside the formal boundary of the pit (shown as a solid line) could be explained as an interpolation effect, since there is no profile in between profiles 15n and Os. extrapolated t h a t characteristic 15s are reasonably well represented on this map. t h e conductive The t h e t h e t h e about 5 m. These pictures provide a reasonable volume image of t h e bulk Pit conductivity. application. The principal assumption is t h a t TDEM migration theory was developed for the 2D case, yet it was applied to 3D data. We assumed a known background resistivity, while there were no reliable data available for shallow depth range (0-4 m), because t h e useful image produced by "smoke rings" starts from the 4 m depth. However, comparison of t h e results of t h e time domain electromagnetic migration and resistivity imaging with the known object locations demonstrate that this method could be used to determine t h e structure of anomalous resistivity distributions in the cases where the conductivity distribution is 3D. 66 presented in Figure 37 and Figure 38. The heavy lines on these maps show the known boundary of the pit. Let us analyze, for example, the migration resistivity map at the 5 m depth (Figure 38). We observe on this map several conductivity anomalies that correlate very well with the pit's the conductive A 3D image was obtained from the set of 2D profiles using simple linear interpo­lation. Thus, the conductive structure from profile Os could be formally extrapolated outside the pit by the interpolation code. However, comparing the map on Figure 38 with the known overall Pit scheme (Figure 30) we can see that in general the prin­cipal boundaries of the pit and its internal structure, that is the two characteristic conductive zones along the profiles Os and 35s and the slightly conductive zone along Also, 3D resistivity images were constructed on the basis of resistivity maps at the depth 3, 5, 7 and 9 m (see Figure 39). One can see on this map that conductive structures appear near 3 m deep, and disappear between 7 and 9 m deep. The maximum of the conductive structures on the resistivity images lies at the depth about the conductivity. At this point it will be appropriate to outline all the assumptions made in this that developed background the the the method the 67 -10 0 10 20 40 50 60 70 80 90 Distance (In) "." .,·:<:~::g~e;LI·~·w .. ~."",~", I 9 21 33 45 57 70 82 94 106 118 130 Ohm-m Figure 37: Apparent resistivity map of Cold Test Pit at the depth 3 m. 68 -10 0 10 :§: 20 8 i~s 30 -10 o 10 20 30 40 50 70 80 Distance (rn) 21 33 45 57 70 82 94 106 118 130 OhlTl-rn Figure 38: Apparent resistivity map of Cold Test Pit at the depth 5 m. 69 depth=3m depth=5m depth=7m depth=9m Ohm-m image. I .. / ·"'·'·"·"""'''''''''·1· ; i i i i i i i i 20 28 36 44 52 60 68 76 84 92 100 Figure 39: Horizontal slices of Cold Pit resistivity image. 5 . C O N C L U S I ON migration to synthetic and field data. Results suggest that TDEM migration and resistivity imaging make it possible to locate anomalous geoelectric zones and estimate their resistivities. The resistivity images can be obtained from the vertical component of magnetic field Hz, which is the component routinely measured in practice. The measurements t h e of estimating t h e background resistivity affects the accuracy of t h e migration method. of testing t h e t h e Radioactive demon­s t r a t ed d a t a field conditions. 5. CONCLUSION The main goal of this thesis research was to apply the method of TDEM migration Hz, The TDEM migration is a stable and fast method of geoelectrical imaging. It provides better resolution of local objects than conventional ID imaging. It does not require repetitive costly forward modeling as do 2D inversion techniques. The application of migration is limited by the necessity for spatially dense measurements of EM responses along the profile or over an area on the surface of the earth. Errors the the method. The results oftesting TDEM migration at the Cold Test Pit within the Radioactive Waste Management Complex at the Idaho National Engineering Laboratory demon­strated the method's effectiveness for interpretation of data acquired under real field A P P E N D I X A A P P A R E N T R E F L E C T I V I T Y F U N C T I ON t h e coefficient field find field, we show here. crn(x,z) contrast on the horizontal boundary between two layers. Consider a 2D model of the electromagnetic field (E-polarization). a model with a slow horizontal variation of both the conductivity and the field we can represent the different components of t h e electromagnetic field in the frequency domain approximately by the following formulas (Smith and Booker, 1991; Zhdanov et al., 1994) Z(X,Z,UJ) = Qd z(x,z,u)eik»z + Qu z{x, z,Lo)e~^z, x(x, z, J) = Qd x(x, z, uj)+ Ql(x, z, u)e~ik"z, y(x,z,u) = x,z,u>)eiknZ + Q^(x,z,u)e~iknZ, Qx\y,z v a i T with depth, kn x, z, LJ) = yJiu)/ioan crn conductivity. function field: HZ(X,Z,LJ) = Qd z(x,z,uj)eiknZ, H^{X,Z,OJ) = Qu z(x,z,to)e~iknZ, Hd x{x,z,u>) = Qd x{x,z^)eik»z, Hu x{x,z,u;) = Qu x(x,z,u)e-ik»z, E*(x,z,u;) = QdJx,z,u)eik»z, Eu(x,z^) = Q"(x,z,u>)e-ik"z, APPENDIX A APPARENT REFLECTIVITY FUNCTION The ratio of the primary and secondary fields is equal to the reflection coefficient at a geoelectrical boundary. In the migration concept the actual secondary field is replaced by the field in the reverse time. However, it is still possible to find a reflectivity coefficient using a ratio of migrated secondary field and primary field, as Let us consider, for the sake of simplicity the two-layered model with a slowly varying conductivity an (x,z) within each layer and a sharp conductivity contrast S field In the following H (x z w) Qd(x z w)eiknZ + QU(x z w)e-iknZ z , , z' , z , , , H (x z w) = Qd(x z w)eiknZ + QU(x z w)e-iknZ x , , x, , x , , , (33) E (x z w) = Qd(x z w)eiknZ + QU(x z w)e-iknZ Y , , Y' , Y , , , where the coefficients Q~:~,z vary relatively slowly with the depth, kn (x, z, w) ViWlloan (x, z) is a wave number and an (x, z) is a background conductivity. Here the terms associated with the downgoing exponential function correspond to the primary field and the terms associated with the upgoing exponential function correspond to the secondary field: Hd(x Z w) Qd(x z w)eiknZ z , , z' , , HzU ( x,z,w ) -- QUz ( x,z,w ) e -iknz , Hd( ) _ Qd( ) iknz xx,z,w-~xx,z,we , HxU(X 'Z'W) _- QUx(X ,z,w)e -iknz , (34) Ed(x '" w) = Qd(x z w)eiknZ EU(x z w) = QU(x z w)e-iknZ y ,..." y' , 'y' , y' , , component t h e horizontal x,d(x))} two first-layer t h e Ey(x, d, UJ) = Qd(x, d, to)e^d + Qfa, d, uj)e~^d, y(x, d,u) = nQd(x, d,u;yk"d - nQu y(x, d,u)e~ik"d, Ey(x:d,Lo) = Qy(x,d,uj)elkn+ld, E'{x,d,u) ikn+1Qz(x,d,u)elkn+ld 1 t h e Ey E'y the right-equal. system < ^ = P(x,d)e2^d J Qy where yjan {x,d)-yfan+1 (x,d) (3 (x, d) = (38) ^jo~n (x,d) + yj<rn+i (x, d) t h e reflectivity coefficient. frequency function h e t h e E^(x.z.u) Q™(x,z,uA , . E;(x,z,u) Qd y{x,z,u) pua(x,d,w)=P{xtd). (40) So, at the geoelectrical boundary the apparent reflectivity function is equal to the r u e reflectivity coefficient. 72 We will analyze the behavior of the horizontal camp,onent of the electric field at the quasi-hori zontal boundary {S : (x ,d(x» } between L,WQ layers. On the first.-Iayer side of the boundary we have: (35) E' (x d w) ik Q'(x d w)eik.' - ik Q"(x d w)C;k.' II" ny" ny " , while on the second-layer side we have: (36) E1'1 (x ' d, w) = ikn +,Q ;(x, d , w)e;k"+" , where " prime" denotes the vertical derivative of the electric field. On the boundary S in the case of E-polarization both components Ey and E~ are continuous. Therefore, t.he corresponding right.-hand sides of (35) and (36) are equaL Solving this sys tem of equations, we find: (37) (38) is the reflect ivi.ty The [requeney domain apparent reflectivity fun ction as t he ratio of the secondary and primary fields is: (39) According to (37) at the boundary: /3wo(x, d,w) =/3(x, d) . (40) SO, t he geoelect rical reflecti vity t rue refl ectivi ty coeffic ient. reflectivity al., is to: EF(z,z,u>) = Q«(x,z1u>)e-k°>*, 41) E?*(x,z,U) = Qi y(x,z,w)e-k"", where km u>) = \jiu[iQUm is a migration wave number. function as x E?'(x,z,w) Q*(x,z.w) Eyv(x,z,u) Qd y{x,z,uj) Calculating function apparent reflectivity function using and (43) gives: (3Z(x,z^) = P„a(x,z,u>)e2ik»z. (45) S : PZ(x,d,u) f3wa(x,d,u)e2ik"d p(x,d)e2ik»d. proportional t h e boundary: 73 Now let us find the way to calculate the apparent reHectivity function from the migrated fields. According to the definition (Zhdanov et aI. , 1994a), the migrated secondary field i$ equal La: em' (x ' w) - Q"(x ' w)e-km' .Jy ,~, - lJ ,-, , (41 ) and the migrated primary field is equal to: (42) km (x, z, w) iiwJ.!oum (x, z) We can introduce the frequency domain migration apparent reflectivity function as : (43) Ca.\cu JaLing the ratio of the apparent reflectivity funelion and migration apparent fun ction (39) Thus we obtain: '(;3:..".,a+(X',- z'-,-w+) = e -21k .. ~• (3~~(X, z,w) According to eq. (40) at the boundary S: i•J•Wn ((xJ" d w) = (iW G(x" d w)e2" "' = p• (x, d) e'ikod . (44) (45) ( 46) This means that the migrated secondary and primary fields are proport ional at lhe boun dary: 74 E™s(x, d,w) 0 x, E™p(x, d,to)e Equation 47) reflec­t i v i ty time : ff?{z,z) = l%(z,z,t = 0) = E™(x,z,t = 0) (48) D(x, z) P^(x, t) t h e migration reflectivity obtained from ft™a(x, d,cj) inverse Fourier transform from the frequency domain to the time where x, z) +oo E™*>(t)<p(z,t)dt, ip(z,t) - a-exp (49) (50) ft™ z) reflectivity (3 t h e = /»(*,<Q- (52) Em,(x d w) = f3 (x d) Em,(x d w)e'''"' 11 1 1 'Y 1 1 • (47) Equat ion (4 7) can be used for calculating the time domain migration apparent reflec­tivity function at zero time: (48) Here ,8:;; (x ,z, t) is the time domain mi gra~ ion apparent renectivity function obtallled P:::u( .'t,d,w) by Pouder tran sform domail1 t.he domain, and 1.+00 D(x,z) = 0 E;;'P(x,z, t)<p(z, t)dt, (49) (50) and (51) The function fJ;: (x, z) is equal to the actual reAectiyjty coefficient (3 exactly at the position of the boundary: f3::' (x, d) = f3 [x, d) . (52) A P P E N D I X B M I G R A T I O N A P P A R E N T R E S I S T I V I TY primary field and considering a two-layered model, one can have an expression con­necting the reflectivity function with the amplitude of t h e migrated field. Using this relationship one can express the anomalous resistivity via the migrated field ampli­tude. Although derived for the plane wave and two-layered model, this expression can be used in the general case, thus leading to a migration apparent resistivity. coefficient Qd u>) in equation (36) does not depend on to and very slowly changes with Qd y(x,z,u>)nQ0(x). 6-pulse Ed y(x,z,uj) = QQ{x)eiknZ. E%{x,z,u) Qo(x)(312(x,d)e2lk"de^z. t h e to E™s(x, z, to) = Q0 (x) (3{x,d) e 2 ^ d e - k - z . APPENDIX B MIGRATION APPARENT RESISTIVITY As shown in Appendix A, the apparent reflectivity function is equal to the re­flectivity coefficient at a geoelectrical boundary. Substituting a plane wave for the the field. Let us consider the model that is based on the assumption that the coefficient Q~ (x, z,w) wand z : Q~(x,z,w) ~ Qo(x). (53) These conditions for example take place in the case of constant conductivity of the first layer and if the primary field can be approximated by the plane wave with a b-pulse waveform on the surface of the earth. Under these conditions we can write approximately: (54) Therefore from (37) and (36) we have: E; (x, z, w) = Qo (x) /3l2 (x, d) e2ikndeiknz. (55) According to definition, the migrated secondary field in the frequency domain is equal (56) 76 Fourier transform t h e time-domain, find the I r+oo E™s(x, z, t) = -Qo (x) p (x, d) / e2tKde-kmZe-iujtdw. (57) integral tabulated = 0. calcula­tions, where 7 = o-m/crn constant. us depth extremum t h e field. t h e vertical migrated t h e = if 7 = value field the extremum is equal to B^(X,d,t = 0) = - ^ ^ . From the last formula we can find 3 (x,d) P(x,d) = ~^°mf.E™{x,d,t = 0), V0Q0 (x) calculate the resistivity Pn+l (x) i + (3(x,dyl 1 -0(x,d) Expression (62) gives the exact solution for determining pn+1 only for the simple under introduce resistivity, which If we apply the inverse Fou.rier t ransform from the frequency domain to the time­domain , we nod the migrated secondary field in t he time domain: (57) The last in tegral can be reduced to a ta,hulated one, if t O. Omitting long calcul a- Lions, we write: (58) where., u.,../u" is the migration constant. Let find now the deptu of the extremulll of the migration fi eld . Since the verticaJ derivative of the migra.ted field is equal to zero exactly at the depth of the boundary z d, jf the migration constant i 4/3. The vaJue of the migration Held in cao (J(x,d) f3 ( d) 47rI'Um", Em,( d ) x, = - r.:; ()' 11 X, ,t = 0 , v3Qo x and calcu late t.he resis tivi ty of the second layer: [ 1 + f3(x, d)]' po+1 (x) = l-f3(x,d) pno (60) (61) (62) P,,+l (x) tbe model u_nder consideration. However, following a traditional approach to electro­magnetic sounding, we can int roduce the migration apparent resistivity) wbich is 77 formula: 1 2 t h e (3™ (x, z) is given by eq. (48). h e expression can equations 77 determined by the same formul a.: ( ) [ 1+,8:;' (x,Z) ] , Pm X,Z = l -{3::(x,z) p" . (63) In the last expression f3: (x, z) is gi ven by eq. (48) . Note in t he conclusion that ex pression (61) ca.n be obtained directly from equations (48) and (49). R E F E R E N C E S L, transient 2861. geomagnetic Amsterdam. J . F . , moveout-J . F . , Oxford,Eng. d a t a electromag­netic T T P Chem-two-dimensional Geophysics, 4 3 , 49-76. J.J.three-dimensional Elsevier,Amsterdam. REFERENCES Adhidjaja, J.I., Hohmann, G.W. and Oristaglio, M.L., 1985, Two-dimensional transient electromagnetic responses, Geophysics, 50, 2849-286l. 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 deep geomagnetic sounding, Elsevier, Amsterdam. Cai, W., Qin, F., and Shuster, G.T., 1995, Electromagnetic velocity inversion using two-dimensional Maxwell's equations, submitted to Geophysics. Claerbout, J.F., and Doherty, S.M., 1972, Downward continuation of moveout-corrected seismograms, Geophysics, 37, 741-768. Claerbout, J.F., 1985, Imaging the earth's interior, Blackwell Sci. Publ., Oxford,Eng. Eaton, P.A., and Hohmann, G.W., 1989, A rapid inversion technique for transient elec­tromagnetic soundings, PEPI, 53, 384-404. Farquharson, C.G., and Oldenburg, D.W., 1993, Inversion of time-domain electromag­netic data for a horizontally layered earth, Geophys. J. Int., 144, 433-442. Lee, S., McMechan, G.A., and Aiken, L.V., 1987, Phasefield imaging: The electromag­netic equivalent of seismic migration, Geophysics, 52, 679-693. Macnae, J., and Lamontagne, Y., 1987, Imaging quasi-layered conductive structures by simple processing of transient electromagnetic data, Geophysics, 52, 545-554. MacLean, H.D., 1993, Time domain electromagnetic survey of three waste burial pits at INEL radioactive waste management complex, TTP AL 911201-G2, Chern-Nuclear Geotech, Inc. Morse, P.M., and Feshbach, U., 1953, Methods of theoretical physics, McGrawHill Book Co., New York. Nabighian, M.N., 1979, Quasi-static transient response of a conducting half-space - an approximate representation, Geophysics, 44, 1700-1705. Oristaglio, M., and Hohmann, G., 1984, Diffusion of electromagnetic fields into a two­dimensional earth: A finite-difference approach, Geophysics, 49, 870-894. Schneider, W.A., 1978, Integral formulation for migration in two and three dimensions, 43, Smith, J .T., and Booker, J .R., 1991, Rapid inversion of two- and three-dimensional inversion of magnetotelluric data, J. Geoph. Res.,96, 3905-3922. Tarantola, A., 1987, Inverse problem theory. Methods for data fitting and model pa­rameter estimation, Elsevier,Amstcrdam. Traynin, P.N., 1995, Analysis of rapid inversion techniques and time-domain electro- t h e method 2 7 1, development t h e Geoelec, 35 Frenkel, electromagnetic Nauk 299, Amsterdam. electromagnetic 0 . , 0 . , tim e 79 magnetic migration results for transient electromagnetic soundings, In preparation. Xiong, Z., 1992, Electromagnetic modelling of three-dimensional structures by the method of system iterations using integral equations, Geophysics, 57, 1556-1561. Zhdanov, M.S., and Frenkel, M.A. 1983, Migration of electromagnetic fields in the so­lution of geoelectrical inverse problems, Doklady Akademii Nauk SSSR, N3, 271, 589-594. Zhdanov, M.S., and Frenkel, M.A., 1983a, Electromagnetic migration. The development of the deep geoelectric model of the Baltic Shield, Oulu (Finl.): University of Oulu, part 2, 37-58. Zhdanov, M.S., and Frenkel, M.A., 1983b, The solution of the inverse problems on the basis of the analytical continuation of the transient electromagnetic field in reverse time, J.Geomagn. and Geoelec., 747-765. Zhdanov, M.S., 1988, Integral transforms in geophysics, Springer-Verlag, New York. Zhdanov, M.S., Matusevich, V.U., and Frenkel; M.A., 1988, Seismic and electromagnetic migration, (in Russian) Nauka, Moscow. Zhdanov, M.S. and Podchufarov, V.V., 1988, Migration of scalar components of electro­magnetic fields for geoelectrical problems, Dokladu Akademii N auk SSSR, N3, 331-335. Zhdanov, M.S., and Booker, J.R., 1993, Underground imaging by electromagnetic mi­gration, SEG Expanded Abstracts 63rd Annual meeting, Washington D.C., 355-357. Zhdanov, M.S., and Keller, G.V., 1994, The geoelectrical methods in geophysical explo­ration, Elsevier, Amsterdam. Zhdanov, M.S., Traynin, P., and Booker, J.R. , 1994a, Principles of electromagnetic migration, (submitted to Geophysics). Zhdanov, M.S., Traynin, P., and Portniaguine, 0., 1994b, Migration and analytic con­tinuation in geoelectric imaging, SEG 64th Annual Meeting, Expanded Abstracts, Los Angeles, 357-360 Zhdanov, M.S., Traynin, P., and Portniaguine, 0., 1995a, Resistivity imaging by time do­main electromagnetic migration (TDEMM), (submitted to Exploration Geophysics). Zhdanov, M.S., Traynin, P., Portniaguine, O.N., and MacLean, D. 1995b, Time domain electromagnetic migration in INEL RWMC Cold Test Pit characterization, Proceed­ings of the Symposium on the Applications of Geophysics to Engineering and Envi­ronmental Problems, SAGEEP 95, Orlando, Fl., 919-924.