R E C O N S T R U C T I O N OF CONDUCTIVITY AND SUSCEPTIBILITY F R O M T H E INVERSION OF E M DATA By Zhiyi Zhang B . Sc. (Geophysics) Changchun University of Earth Sciences M . Sc. (Geophysics) China University of Geosciences A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE. FACULTY OF GRADUATE STUDIES EARTH & OCEAN SCIENCES We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA March 1997 © Zhiyi Zhang, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Earth & Ocean Sciences The University of British Columbia 129-2219 Main Mall Vancouver, Canada V 6 T 1Z4 Date: Abstract Both electrical conductivity and magnetic susceptibility are important physical parameters in exploration geophysics, and information about their distributions can be used to determine subsurface structures, and to detect mineral deposits and other natural resources. The ultimate goal of this thesis is to recover conductivity and susceptibility from the inversions of electromagnetic (EM) data from various loop-loop systems. A large number of frequency-domain E M ( F E M ) data are taken in E M surveys by using different loop configurations, and the inversions of these data may provide independent information about the geological targets. Previous work on the inversions of E M data has involved only horizontal loop sources. Consequently, data measured with other coil systems have been constantly rejected from the inversions. In Chapter 2, I investigate the effect of coil configurations on the inversion and develop an inverse algorithm to invert E M data from different coil systems. E M data can also be measured in the boreholes. Large loop systems which measure transient E M ( T E M ) data on the ground or in the borehole have found increased application in exploration geophysics. However, the inversion of borehole T E M data has not been fully addressed. In Chapter 3, I investigate whether, and how, the use of borehole data in inversions enhances the recovered models. In geophysical explorations, E M responses are a function of the geometry, conductivity, and susceptibility. The influence from magnetic susceptibility on E M data has long been appreciated, but no existing literature has been found about the reconstruction of susceptibility through rigorous inversions. In most cases magnetic susceptibility has been treated as a source of "contamination" in the inversions, and people have been trying very n hard to eliminate that "contamination" by truncating or disregarding the inphase E M data in carrying out inversions. In doing so, useful information about the distribution of magnetic susceptibility is wasted, and the recovered conductivity models become less reliable. In Chapter 4, I study the effect of susceptibility on the data, and reconstruct susceptibility from the inversion when the conductivity distribution is specified. The problem with the individual inversions in Chapters 2, 3 and 4 is that accurate information about either conductivity or susceptibility is required in order to recover the other. Thus, it is necessary to explore the possibility of reconstructing 1-D conductivity and susceptibility simultaneously from the inversion of the E M data. In carrying out simultaneous inversions in Chapter 5, I minimize a global model objective function, which includes both conductivity and susceptibility, subject to the data constraints. The final conductivity and susceptibility models are obtained by adjusting the parameter that controls the relative weighting between the two terms in the model objective function. Much of the work in Chapter 2 to Chapter 5 is done in 1-D environment, while geological targets are usually 3-D. Idealy one would like to carry out 3-D inversion to obtain information about the 3-D targets. Currently, however, full 3-D rigorous inversions of E M data are computationally prohibitive, and approximate 3-D inversions are necessary. In Chapter 6, I develop a linear mapping that can be used to in interpret the data collected in a 3-D environment. The algorithm is applied to a field data set collected over Mt.Milligan. A l l algorithms have been tested on both synthetic and field data sets. Those tests show that the algorithms are robust to different error assignments, even when the error is correlated. The recovered conductivity and susceptibility models from the inversions of field data have provided useful geological information. in Table of Contents Abstract " List of Tables viii List of Figures ix Acknowledgement xvi 1 Introduction 2 1 Recovering 1-D conductivity from the inversion of E M data 12 2.1 Introduction 12 2.2 Forward modeling of the coplanar E M system 13 2.2.1 16 2.3 2.4 2.5 Forward modeling for other coil configurations Calculation of the sensitivities 18 2.3.1 Calculation of the sensitivities for the coplanar system 19 2.3.2 Calculating sensitivities for other coil systems 20 The inverse algorithm 21 2.4.1 23 Joint inversion of the coplanar and coaxial data Synthetic examples 24 2.5.1 1-D examples 25 2.5.2 3-D examples 32 2.6 Field data examples 55 2.7 Conclusion 56 iv 3 1-D Inversions of the Surface and Borehole Transient E M data 60 3.1 Introduction 60 3.2 Data conversion 63 3.2.1 Conversion of frequency-domain data to time-domain data . . . . 64 3.2.2 Conversion of transient data to frequency-domain data 64 3.3 4 5 Forward modeling 68 3.3.1 Solutions for the source-free regions 71 3.3.2 Solution in the source-containing layer 72 3.3.3 Calculation of the sensitivities 77 3.4 Synthetic data example 78 3.5 Field data example 86 3.6 Conclusions 93 Recovering susceptibility from 1-D Inversion of E M data 94 4.1 Introduction 94 4.2 Inversion algorithm 96 4.3 Calculation of the sensitivities 97 4.3.1 98 Adjoint Green's function solution 4.4 Synthetic examples 104 4.5 Field example 109 4.6 Conclusions 114 1-D Simultaneous Inversion of E M data 117 5.1 Introduction 117 5.2 Methodology 120 5.2.1 122 5.3 The trade-off between conductivity and susceptibility Numerical results 127 v 6 5.4 Field data example 139 5.5 Summary 148 Approximate inversion of 3-D E M data 150 6.1 Introduction 150 6.2 The methodology 153 6.3 6.4 7 Calculation of the sensitivities ! 156 Approximate sensitivities 157 6.4.1 The electric field and its adjoint field within a half-space 158 6.4.2 The calculation of the numerical integral in Equation (6.23) . . . 158 6.5 Synthetic examples 162 6.6 Field data example 166 6.7 Conclusions 170 Summary 176 References 182 Appendices 193 A Zero frequency solution for Hz over a half-space 193 B Adjoint Green's function solution for the sensitivities 194 C Hankel transformation for equation (B.8) 196 D Equivalence of equation (4.12) and (4.15) 198 E 199 Computation of the discrete sensitivities vi F The trade-off between a and /x in the sensitivities vii 201 List of Tables 5.1 D(a, fi ), 0 D(a, / x ) , A D ( c r , /x), a n d D(<x, fi) over a 0.01 S / m a n d 0.1 S I u n i t 0 half-space 5.2 T h e c o m p a r i s o n o f t h e s e n s i t i v i t i e s for t h e c o n d u c t i v i t y a n d s u s c e p t i b i l i t y , f r o m a l a y e r o f 2 m t h i c k at a d e p t h of 4 8 m , i n a 0.01 S / m half-space s u s c e p t i b i l i t y is 0.1 SI u n i t vm whose List of Figures 2.1 Four commonly used coil configurations 11 2.2 Geometry of the coplanar coil system 13 2.3 The secondary magnetic field over a half-space of 0.01 S/m 17 2.4 The results of the inversions of the 1-D data, at 10 frequencies, from different coil systems. The standard deviations are set to 10% of the data strength for all data 2.5 25 The results of the inversions of the 1-D data from different coil systems. The standard deviations are lppm for the coaxial data, and 2ppm for the data from the other three coil configurations 2.6 The influence on the data from the resistive and the second conductive zones in the model used in the first example 2.7 27 28 The recovered models from the inversions of the 1-D data at 110, 7040, and 56320 Hertz 30 2.8 The 3-D model used to generate the synthetic data 32 2.9 The real component of the vertical coplanar data in ppm 34 2.10 The imaginary component of the predicted and observed vertical coplanar data in ppm 35 2.11 The x-z cross-sections of the reconstructed conductivity at (a) y=250m, (b) y=350m, and (c) y=450m, from the inversion of the vertical coplanar data 36 ix 2.12 The recovered model from the inversion of the vertical coplanar data, at (a) x=250m, (b) x=350m, and (c) z=350m 37 2.13 The x-y plan-views of the recovered conductivity, from the inversion of the vertical coplanar data, at six different depths 2.14 The real component of the perpendicular data in ppm 38 39 2.15 The comparison of recovered models from the inversions of the coplanar, perpendicular, and coaxial data at section y=350m 40 2.16 The real and imaginary components of R over the section at y=250m, at 880, 7040, and 56320 Hertz 43 2.17 The real (the solid line) and imaginary (the dashed line) components of R at y=350m 44 2.18 The profiles of R at y=450m 45 2.19 The recovered model from the inversion of the coaxial data 47 2.20 The recovered model from the inversion of the horizontal coplanar data. 48 2.21 The recovered models and the data from the inversions of the vertical coplanar and horizontal coplanar data 49 2.22 The recovered model from the joint inversion of the coaxial and the horizontal coplanar data 51 2.23 The observed and the predicted coplanar data from the joint inversion of the coaxial and coplanar data 52 2.24 The recovered conductivity models over section Y9600 at M t . Milligan. 3.1 The conversion of F E M data to T E M data due to a step turn-off, for a central loop system over a half-space with 0.01 S/m conductivity 3.2 55 64 The conversion of step turn-off T E M data to F E M data for a central loop system over a 0.01 S/m half-space 65 x 3.3 The conversion of T E M data to F E M data for a central loop system whose loop source has a radius of 50m. The impulse response of the magnetic field which spans l O - 1 0 to 1 0 _1 seconds is used as the input data 3.4 The geometry for the inversion of borehole E M data 3.5 Frequency and transient responses over a half-space of 0.01 S/m for a central loop borehole system 3.6 66 68 73 The transient responses, as a function of the receiver position and time, due to a step turn-off over a half-space of 0.1 S/m for a central loop borehole system 3.7 74 The sensitivities for both surface and borehole central loop configurations in the time- and frequency domains 3.8 The sensitivities for borehole central loop configurations in the time- and frequency-domains in a half-space of 0.01 S/m 3.9 76 78 The models and the data from the inversions of the surface and borehole T E M data 80 3.10 The results from the inversion of the central loop data from the borehole configuration 82 3.11 The recovered model and the data from the joint inversion of the central loop data from both surface and borehole configuration. The receiver for the borehole configuration is placed at 150m along the axis of the borehole. 83 3.12 The recovered model and the data from the joint inversion of the central loop data from both surface and borehole configurations. The receiver for the borehole configuration is placed at 150m and 250m along the axis of the borehole 84 3.13 The results from the inversion of surface S I R O T E M data xi 86 3.14 The results from the inversion of borehole S I R O T E M data. The receiver is at 155m in the borehole 86 3.15 The recovered model and the observed and predicted data from the inversion of the S I R O T E M data collected at 230m in the borehole 87 3.16 The recovered model and the observed and predicted data from the joint inversion of the borehole (at 155m) and the surface S I R O T E M data. . . 88 3.17 Compatability of the surface and borehole data at 230m depth 89 3.18 Compatability of the surface and borehole data at 155m depth 90 4.1 Nonlinear mapping for positivity constraint 95 4.2 The absolute value of sensitivities (in logarithmic scale) over a half space. 101 4.3 Inversion of ground system data 103 4.4 Inversion with correct knowledge of conductivity structure 104 4.5 Effect of incorrect knowledge of the conductivity distribution on the inversion 106 4.6 D I G H E M data from M t . MiUigan, at section Y9600 108 4.7 Inversion of D I G H E M data from M t . Milligan, section Y9600 Ill 4.8 Recovered susceptibility models from inversions with different error schemes. 113 5.1 Individual inversions with accurate and inaccurate information about conductivity or susceptibility 117 5.2 The relationship among D(cr,fj,), D(a,fi ), 5.3 The result of the simultaneous inversion with s = 6 5.4 The results from the inversion with s — 20 5.5 The inversion with s = 0.1 129 5.6 The 3-D model 130 0 xii D(<7,// ), 0 and AD(<r,fi). . . . 123 126 127 5.7 The real components of the predicated and observed data for the 3-D model. Panels (a), (b) and (c) show the predicated data at 110, 7040 and 56320 Hertz. Panel (d), (e) and (f) plot the observed data at the same frequencies 5.8 132 The Imaginary component of the predicated and observed data at 110, 7040 and 56320 Hertz for the 3-D model 5.9 133 The x-z cross sections of the recovered conductivity from the inversion of 3-D data 135 5.11 The x-y plan-views of the recovered conductivity and susceptibility. . . . 136 5.12 The real components of the predicted and observed data 139 5.13 The imaginary components of the predicted and observed data 140 5.14 Recovered conductivity models from the simultaneous inversion of HSS A E M data 142 5.15 The x-y planview of the recovered susceptibility models from the simultaneous inversion of A E M data and the inversion of aeromagnetic data. . . 144 5.16 The cross-section of the recovered susceptibility models from the simultaneous inversion of A E M data and the inversion of aeromagnetic data over section 12600 at HSS 145 6.1 The geometry of the 3-D problem 149 6.2 The local coordinate system for the calculation of sensitivities 156 6.3 The order of the nodes for interpolation of E • G in each cell 158 6.4 The real and imaginary components of the sensitivities due to a layer between z = 0 and 2m depth inside a O.OlS/m half-space xm 160 6.5 The real [panel (a)] and imaginary [panel (b)] components of the sensitivities due to the layer extending from z — 45 to z = 50m inside a O.OlS/m half-space 161 6.6 (a) Misfit curve; (b) Model norm as a function of the number of iterations. 162 6.7 The x-y plan-views of the recovered model and apparent conductivity. . . 164 6.8 The x-z cross-sections of the recovered model at y=250, 350, and 450m. . 165 6.9 D I G H E M data at M t . Milligan 168 6.10 The X - Z cross section of the recovered conductivity, with overlayed geology (white lines), from the 3-D approximate inversion of D I G H E M data at M t . Milligan 169 6.11 The X - Y plan-view of the recovered conductivity from the 3-D approximate inversion of D I G H E M data at M t . Milligan 170 6.12 The X - Z cross section of the best-fit conductivity, with overlayed geology (white lines), from the 3-D approximate inversion of D I G H E M data at M t . Milligan 171 xiv Acknowledgement I would like to thank my supervisor Prof. Doug Oldenburg for his advice and encouragement in every aspect of my life in Canada and my research. Without Doug's moral and financial support, it would have been impossible for me to carry through my research here at U B C . I also thank Dr. Colin Farquharson, Dr. Yaoguo L i , and Roman Shekhtman for their advice and assistance on all matters concerning computers, inverse theory and electromagnetic geophysics. xv To my parents xvi Chapter 1 Introduction Exploration geophysics is applied to obtain information about the subsurface of the earth that is not available from surface geological observations. Because the electrical resistivity of different earth materials ranges over many orders of magnitude, electromagnetic ( E M ) methods can be used to map the subsurface resistivity structure. Early E M methods were largely designed by the Scandinavians and the Canadians for exploration in under-glaciated Precambrian shield conditions, where the resistivities of the host and the overburden are generally high. They did not work well in areas where the overburden or host rock was conductive. The lack of sophistication in data gathering and processing severely limited their exploration depth. Moreover, early E M systems were relatively heavy, cumbersome, and slow in operation. Modern E M methods are characterized by their emphasis on deeper exploration and by the need to acquire measurements of the response of the earth over a broad frequency range. In addition, to process the vast amount of digital data collected in the field and the inherently weaker signals originating from deeper targets, it has become necessary to rely heavily on modern computer technology. Nevertheless, most modern equipment is remarkably portable, considering its sophistication. Nowadays E M methods are widely used in exploration for resources such as petroleum, groundwater, mineral deposits and geothermal energy. They also find application in engineering studies to map structures and to help determine material properties of rocks and in hazardous waste site investigations to map conductive plumes and hydrological 1 Chapter 1. Introduction 2 features. E M methods are also used in archeology studies where layered rocks are usually encountered. E M methods include a variety of techniques, survey methods, applications, and interpretation procedures, which are further complicated by a bewildering array of trade names. Each technique, however, involves the measurement of one or more electric or magnetic field components by a receiver, from an E M source. There are two classes of sources: natural sources, and artificial or controlled sources. Each has its own advantages and disadvantages. For reconnaissance surveys and deep structure studies, natural sources are more economic to use. Due to power requirements, the investigation depth for controlled-source E M methods is limited to a few kilometers or less whereas natural field methods can be used to penetrate through the crust and into upper mantle. However, for studying shallow geological targets, the mobility and the resolution provided by controlled-source E M methods are superior to the natural source E M methods. The artificial sources can be grounded wires or isolated current-carrying loops. Grounded wire antennas are usually straight lengths of wire laid on the surface of the ground and connected to the earth through low resistance electrodes at either end. The transmitter can be connected to the wire at any convenient location. Large loops having dimensions of several hundred meters generally consist of a single turn of isolated wire laid on the ground in the form of a square or rectangle. Loops of intermediate size are generally formed from a multi-conductor cable and are laid on the ground or are suspended in the vertical plane by use of masts. Small loops generally consist of many turns of wire wound on a rigid form, and a ferrite or metal core may also be used. While wire sources can provide greater depth of penetration because the primary field from wires falls off less rapidly at large distance than loop sources, loops are popular because they can work in the presence of resistive overburdens, and can be deployed Chapter 1. Introduction 3 more quickly and easily. More important, since loop sources do not need to be connected to the ground, they can be used in airborne E M surveys, and that is a huge advantage over the wire sources. In most frequency-domain E M ( F E M ) systems, a current having approximately a sinusoidal or a square waveform is driven through the antenna by an amplifier or a switcher. If a square waveform is used, it is possible to simultaneously measure the base frequency and some of its harmonics. In any case the frequency is usually changed in discrete steps to make the measurements. In time-domain E M ( T E M ) systems, the most common waveform is a train of approximately square, bipolar pulses with an off-time between pulses. Other repetitive waveforms such as triangles are also used. The secondary electromagnetic fields induced by the primary fields are measured by different E M sensors. Three types of sensors are used for E M surveys: induction coils, magnetometers, and grounded wires. Wires are used to measure the electric field while loops are used to record the time rate of change of magnetic flux density, dB/dt. The time constant of the receiver coil ( L / R ) should be much less than the earliest time of measurement for T E M systems, and the self-resonant frequency of the coil should be higher than the frequencies of measurement for a F E M system. A variety of fluxgate, SQUID, feedback coil, and other types of magnetometers are sometimes used to measure B rather than dB/dt. For loop sensors, the sensitivity, the self-resonant frequency, and the noise level are the important parameters. The most important parameters of magnetometers are the sensitivity, the frequency response, and the noise level. There are several practical and theoretical differences between B and dB/dt measurements. For detection of very good conductors, B measurements will often give superior results (Mallick, 1978, Gupta Sarma et al., 1976). It is of interest to note that dB/dt devices effectively pre-whiten the noise (Spies and Frischknecht, 1988). In some frequency-domain instruments, the amplitude of one or more components or the ratio 4 Chapter 1. Introduction of two components is measured without reference to the transmitter. More often, some type of phase reference or time reference scheme is used. In principle a cable connection is the simplest. Different combinations of the source and receiver orientations form different coil configurations. Four commonly used loop-loop configurations are the horizontal coplanar (HC), vertical coplanar (VC), perpendicular (PP), and vertical coaxial (CA). In the P P configuration, one loop is oriented with its axis vertical and the other is oriented with its axis horizontal and pointing toward the axis of the first loop. The description of the other configurations are self-explanatory. A large loop is used primarily in the fixed-source mode. For a large loop it is practical to take measurements at any location inside or outside the loop, except in the immediate vicinity of the wire. A special case in which a vertical-axis receiver is put right at the center of the source loop is called central loop system, and it is very popular in T E M surveys. Coincident loop system is a special case where the source loop coincides with the receiver loop spatially. In the display and analysis of E M data, the induction number is an important terminology. Theoretical results are often displayed using induction number. Typically the induction number B = {a\iLcj1yl T is used in the frequency-domain, and the induction 2 number a = (afi/Aty^r or 3 = \/2ct is used in the time-domain. Here r is a geometrical factor which could be the radius of a sphere or a cylinder, or the coil separation, depending on the specific model under study. u> and t are the angular frequency and the time. <r and ft are the conductivity and susceptibility respectively. Theoretical studies on the effects of magnetic dipole sources placed over a conductive earth have been made by Belluigi (1949), Gorden (1951), Bhattacharyya (1955, 1963), Malqvist (1965), and in a series by Wait (1951, 1952, 1953, 1955, 1958, 1962). Slichter Chapter 1. Introduction 5 and Knopoff (1959) have presented some numerical results for homogeneous and twolayer earths. The theory of E M sounding with dipole sources has been later developed by Keller and Frischknecht (1966); Vanyan (1967), Dey and Ward (1970), R y u et al. (1970), Ward and Hohmann (1988). Keller and Frischknecht (1966) and Vanyan (1967) have presented extensive numerical results for the field and mutual impedances on, and above, a two-layer earth for different configurations of two-loop sounding. They mentioned an interpretation scheme with field examples. Dey and Ward (1970) computed the forward responses of a horizontal magnetic dipole (HMD) over a layered earth, by using the Hertz potential. Ward and Hohmann (1988) also studied the same problem by using Schelkunoff potentials and established a set of solutions under the quasi-static assumption. R y u et al. (1970) have developed the complete formulation of the E M response of a multilayered earth excited by a vertical magnetic dipole. Duckworth (1970) has suggested a simplified method to depth sounding in geometric mode for mining problems. Fuller and Wait (1972) have studied the response of a coplanar loop over ground where conductivity varies exponentially with depth. Ryu et al. (1972) discussed the field application of the two-loop sounding method for a multilayer earth. Sinha (1973) presented response curves over a multilayered earth for the case of E M two-loop sounding applicable to airborne surveys with a comparison of different loop configurations. Patra (1970) introduced a convenient method of frequency-sounding with a large loop known as central frequency sounding (CFS). He presented response curves for two-layer models for CFS based upon simplified assumptions and approximations. Later R y u et al. (1970) studied the central induction method and presented numerical results for a two-layer earth. Common techniques used in interpreting E M data include manual curve matching, interactive trial-and-error matching using a computer, and automated computer programs. A l l methods are still used today but automated computer inversion has become Chapter 1. Introduction 6 widespread and more popular than graphic methods, because of the advent of powerful, low-cost, digital computers. Those computer-based inversion methods can be divided into mainly two categories. The first one is model-fitting: computing the parameters of an assumed model supposed to represent that portion of the earth under consideration. Normally the modelfitting problem is overdetermined. Glenn et al. (1973) have discussed development of an interpretation technique of sounding data. In their paper they inverted the data from a two-layer model to recover the conductivity and thickness of the top layer. Later Glenn and Ward (1976) elaborated on the method, emphasizing experiment design: customizing a field survey based on an assumed model to maximize the information content of data. Ward et al. (1976) illustrated the use of the techniques in a ground water problem, comparing results obtained from the use of horizontal- and vertical-loop sources with those from electrical sounding (Schlumberger). Then Tripp et al. (1978), Gomez-Trevino and Edwards (1983), and Raiche et al. (1985) discussed model fitting techniques and applications further, in each case demonstrating the enhanced resolution engendered by joint inversion of electrical and E M soundings. Other publications illustrating the applications of 1-D model fitting for controlled-source E M are Frischknecht and Raab (1984), Spies and Frischknecht (Volume II), and Kaufman and Keller (1983). The main advantage for model-fitting methods is that it is computationally economic. However, because of the complexity of geological targets, usually it is very difficult to represent them adequately with models made of only a few parameters and cells. Further more, the primary concerns in model-fitting inversions is to fit the data, and the nonuniqueness of the inverse problems is not addressed. The other method is rigorous inversion which concerns not only with convergence and robustness, but also with the non-uniqueness of the inversion (Constable et al., 1987; Whittall and Oldenburg, 1992). Usually an underdetermined inverse problem is 7 Chapter 1. Introduction solved in rigorous inversions, and the prototype technique used are variations of the Gauss-Newton technique. Model norms, which can be 1-norm, 2-norm or other types of norms, are minimized subject to the data constraints. The nonlinear inverse problem is solved by linearizing the relationship between the model parameters m , and the observed data D o b s , D o b s = D ( m )+ J(m)£m + O, (1.1) where the predicted data D ( m ) are calculated from the earth defined by model parameters m , and the matrix J is the Jacobian matrix of the data whose elements are given by j ( m ) = Orrij i = i, 2 , M , j = 1, 2 , N , (1.2) where M and N are the numbers of observations and the model parameters respectively. If the remainder O can be neglected, then 8m, the perturbation on the model, can be obtained by using some regularized procedures. The model is updated by adding the perturbation to the starting model which is provided by the interpreter. Successive iterations will eventually reduce the data misfit to the desired tolerant level. Because of the restriction on computer resources, the application of rigourous inversion techniques to the inversion of controlled-source E M data is mainly restricted to 1-D problems. Fullagar and Oldenburg (1984) presented an inversion scheme to invert horizontal coplanar E M data to reconstruct 1-D conductivity. Farquharson and Oldenburg (1994) inverted transient E M data to recover 1-D conductivity. However, there are still many important questions, even for 1-D inverse problems, waiting to be addressed. Until now most work on the inversions of 1-D controlled source E M data involves horizontal coplanar data, and not enough attention has been put on the inversions of E M data Chapter 1. Introduction 8 obtained from other coil configurations. The earth responds to different source orientation differently. Hence E M data collected from different coil systems carry independent information about the targets. So the inversions of data from different coil systems need to be investigated. Another important issue is the effect of susceptibility on the E M data. The main physical property affecting E M measurements is conductivity where earth materials show many orders of magnitude variation, while for many rocks and materials, magnetic susceptibility are either zero or near zero. However, magnetic susceptibility must in many cases be considered in E M surveys. Large susceptibilities, substantially greater than 0, are often encountered in prospecting for mineral deposits containing high concentrations of magnetite or pyrrhotite. High susceptibility may also occur in flood basalts and other mafic rocks. The inphase data in F E M surveys can be severely affected by susceptibility. Negative inphase data in airborne E M surveys are the manifestation of the existence of magnetizations (Fraser, 1970). Those negative inphase data cannot be explained by any purely conductive model. Therefore inversion algorithms which can accommodate susceptibility are needed. Presently the rigorous inversions of 2-D or 3-D controlled-source E M data are still computationally prohibitive. Only a few authors have tried to attack the high-dimensional problems for controlled-source E M data by using rigorous inversion techniques. Ellis (1995) carried out 3-D inversion for a H C system with conjugate gradient method, quasiNewton method, and Gaussian-Newton method. A l l techniques were very time consuming, and because of the limitation on computer memory, the model was severely underparameterized. Newman and Alumbaugh (1995) used a conjugate gradient method to carry out a full 3-D inversion on massively parallel (MP) computer. Going to an M P platform is a necessity in that it allows large models to be reconstructed which are not under-parameterized in a reasonable amount of time. Chapter 1. Introduction 9 Since the main stream of the computers used by industry are PCs or workstations, the 2-D or 3-D inversions of controlled-source E M data are still computationally prohibitive, and cannot be used as practical tools in interpreting the data. As an alternative, low-cost techniques are needed to provide approximate 3-D images of the conductivity structures. Inverting 3-D data with 1-D algorithms can generate useful information about the 3-D structure, but the recovered images, could be biased by 3-D effects. Algorithms based on Born and extended Born (Harbashy, T . M . et al., 1993), quasi-linear (Zhdanov and Fang, 1994) approximations are developed to provide images less biased by 3-D effects than conventional 1-D inversions. This thesis is concerned with reconstructing the electrical conductivity and magnetic susceptibility structures from the inversions of E M data. The final goal is the simultaneous inversion in which the conductivity and susceptibility information is extracted from the E M data simultaneously. I begin in Chapter 2 by illustrating the inversion philosophy and major components of an inversion. In Chapter 2, I describe a Gauss-Newton procedure for a nonlinear problem. This procedure is applied to a frequency-domain E M problem. One dimensional conductivity is recovered from inversions of the horizontal coplanar, vertical coplanar, coaxial, and perpendicular coil systems. This inversion technique is also used to solve the nonlinear inverse problems in Chapters 3, 4, and 5. Chapter 3 is focused on recovering 1-D conductivity structures from inversions of the time-domain data from both surface and borehole configurations. The merits of joint inversions of borehole and surface data are explored. Synthetic and field data sets are inverted to recover 1-D conductivity models. In Chapter 4, I investigate the influence of magnetic susceptibility on E M data and develop an algorithm to recover 1-D susceptibility from the inversions of E M data when Chapter 1. Introduction 10 the conductivity is known. A nonlinear mapping is introduced to add a positivity constraint on the inversion. To conclude Chapter 4, the inverse algorithm is tested on both synthetic and field data sets. In Chapter 5 a 1-D simultaneous inversion code is built up based upon work in previous chapters. B y minimizing a global model objective function composed of individual model objective functions for conductivity and susceptibility, the algorithm recovers 1-D conductivity and susceptibility structures simultaneously from frequency E M data. In this chapter, I also discuss the ways to form feasible solutions, and the trade-off between conductivity and susceptibility. This simultaneous inversion algorithm is tested on synthetic and field data sets. In Chapter 6, I discuss approximate mapping for 3-D E M data. Instead of attacking the 3-D inverse problem directly, I introduce an 3-D mapping based on a Born approximation and solve a linear inverse problem. In order to deal with huge data sets, a subspace technique is used to solve the inverse matrices. This algorithm is tested on both 3-D synthetic and field data sets. The major conclusions of this thesis are summarized in Chapter 7. Finally, how to estimate the errors in the data is an important question because different error assignments can lead to different results. Common error assignments include uniform errors on all data, a percentage of the data, or a constant base value plus a percentage of the data. Since measurement error is usually given in fixed values, a constant error assignment can be used to account for it. 3-D effects are usually proportional to the amplitude of the data, so a percentage error can be used in 1-D inversions in 3-D environments. A small constant error can be added to the percentage error, to keep near-zero data from carrying too much weight in the inversion. Unlike synthetic data, the error in field data is usually unknown. So it is necessary to carry out several inversions with different error assignments, and compare the resultant models with geological Chapter 1. Introduction 11 i n f o r m a t i o n a n d o t h e r p h y s i c a l c o n s t r a i n t s , t o d e t e r m i n e t h e a p p r o p r i a t e e r r o r assignment. A s d e m o n s t r a t e d i n n u m e r i c a l e x a m p l e s , m y a l g o r i t h m is r o b u s t w i t h different error assignments, a n d c a n extract i n f o r m a t i o n about the targets even w h e n t h e error i n t h e d a t a is c o r r e l a t e d . Chapter 2 Recovering 1-D conductivity from the inversion of E M data 2.1 Introduction In electromagnetic surveys, data can be measured with different coil systems. Different orientations of the source probe the earth somewhat differently while different orientations of the receiver measure different components of the earth response. When only a limited number of frequencies are available, the data from different coil configurations provide complementary information about the conductivity. Figure 2.1 shows these four commonly used coil configurations. Most of the applications of those configurations are found in airborne E M surveys, where the transmitter is small enough compared to the coil separation so that it can be treated as a horizontal or vertical magnetic dipole. Ryu et al. (1970) presented a solution for the forward problem of a horizontal loop. Their solution can be easily adopted for a vertical magnetic dipole ( V M D ) source. Dey and Ward (1970) computed the forward responses of a horizontal magnetic dipole (HMD) over a layered earth, by using the hertz potential. They showed that, when the measurements are taken above the earth, the secondary field due to T M mode associated with a H M D can be ignored. Ward and Hohmann (1988) also studied the same problem by using Schelkunoff potentials and established a set of solutions under the quasi-static assumption. Although the problem of 1-D inversions of electromagnetic data has been studied extensively in the literature, most of these studies have centered on the inversion of the 12 Chapter 2. Recovering 1-D conductivity from the inversion of EM data Tx Rx 13 Tx A PP HC Tx CA A R* Tx Rx VC Figure 2.1: Four commonly used coil configurations: the horizontal coplanar (HC), the perpendicular (PP), the coaxial (CA), and the vertical coplanar ( V C ) systems. T x and Rx denote the transmitter and receiver respectively. horizontal coplanar data. In this Chapter I invert data obtained from the horizontal coplanar, vertical coplanar, perpendicular, and coaxial coil configurations to recover 1-D conductivity distribution. I first solve the forward problem for the different coil systems. Based upon Dey and Ward's work (1970), I established a linear relationship among the coaxial, vertical coplanar and horizontal coplanar data. This relationship can be used to detect the violation of the 1-D assumption, and be potentially used as a mapping tool to locate regions with strong variation in conductivity. I propose a simple way to obtain the sensitivities, which are the partial derivatives of the magnetic fields with respect to conductivity, for different coil systems. M y inverse algorithm is tested on both synthetic and field data. 2.2 Forward modeling of the coplanar E M system Horizontal coplanar system enjoys the most applications in geophysical exploration. The symmetry of the problem significantly reduces the complexity in the derivation of the theories for the forward modeling and the sensitivity computing. Solutions for the loop source can be readily reduced to those from a vertical magnetic dipole source, and the Chapter 2. Recovering 1-D conductivity from the inversion of EM data Tx R(m) i a | 14 2 obs * ""i Rx earth's surface CT 1 K 2 Ki °1 Kj+1 ai+i K a M M Z(m)j KM+I CJM+1, Figure 2.2: Geometry of the coplanar coil system. A horizontal transmitter loop (Tx) of radius a is located at height h above the surface of a 1-D earth. The source current has angular frequency to and amplitude I. The receiver (Rx) is situated at a radial distance r from the loop source. The earth is modeled as M layers. The conductivity and susceptibility in the ith layer is denoted as and K{. Z b is the vertical distance between T x and Rx. 0 s reflection coefficient can be used later to establish solutions for other coil configurations. Consider two horizontal coils separated by a distance r . The transmitter is at height h above the earth's surface and carries a harmonic current Ie , where u> is the angular xut frequency. The earth is characterized by a set of horizontal layers whose thickness, conductivity and susceptibility of the ith layer are given by (hi, ai, Ki). The geometry for the forward modeling is shown in Figure 2.2. In the following derivation, I use H and E to denote the magnetic and electric fields, and assume that the electric permittivity e takes on its value in free space 6Q. The circular symmetry allows me to employ a cylindrical coordinate system, reducing Chapter 2. Recovering 1-D conductivity from the inversion of EM data 15 Maxwell's equations to the following: iwfijHr = 1d = — — (rEe), r or iu H H (2.1) oz z (2.2) and dH dz dH dr r z (iuje + o-j)E 0 e + J, (2.3) 3 where J = I(uj)S(r - a)S(z ), s and i = obs z b is the distance between the source plane and the receiver plane. B y 0 s ehminating the magnetic field H and H from equations (2.1), (2.2) and (2.3), I obtain r z an inhomogeneous scalar equation / d d Id + + " a ydz ' dr r dr 2 2 2 2 1 \ 9 + j) r , Ee{r,z,u>) = iujfi I(uj)6(r - a)S(z ), k 2 3 0 obs (2.4) where k - = u> fij€ — iu;[ij(Tj. Equation (2.4) can be converted into Hankel transform 2 2 0 domain. The partial differential equation for the electric field, in the Hankel transform domain, is d2 dz 2 ^ E(\,z,w) — iio\LjaZ\{\a)S{z^, (2.5) where A is the Hankel transformation parameter, J i is the Bessel function of first order, and u - = A — k . After inverse Hankel transformation, the secondary electrical field is 2 2 2 given by r°° A 7^ Z = -iu>^aI(u) ^ _ - ^ _ ^ « o ( * o . . - 2 f c ) j ( A ) J ( A r ) A i A , 2 E(r,z ^-) ohs e 1 a 1 (2.6) where the input impedance of the j t h layer is given by (Wait, 1962) •_ J Z ^ + ^-tanhK-fe,-) Z + ^'+ tanh(^ )' 1 j j 1 ' Chapter 2. Recovering 1-D conductivity from the inversion of EM data 16 and the intrinsic impedance Zj is denned as (2.8) Uj In the half-space at the bottom there is no up-going wave and hence the input impedance is equal to the intrinsic impedance. That is, Z M = ZM- After inverse Hankel transformation, I obtain the secondary magnetic field in the frequency-domain (2.9) zw/io Jo The induced secondary voltage measured in a receiver coil is the time derivative of the magnetic flux and is expressed in the frequency domain as (2.10) where DS is the effective area of the receiver. Field data sets take on different forms. The responses can be the secondary magnetic fields or voltage, or they can be total magnetic fields or voltage; these latter require the inclusion of the primary field. When secondary fields or voltage are used, the data are usually normalized by the primary field and provided in parts per million (ppm) of the primary field. Responses in a field survey are acquired at a number of different frequencies and at each frequency both inphase and quadrature phase (or real and imaginary) data can be obtained. The phase determination is made with respect to the primary magnetic field. 2.2.1 Forward modeling for other coil configurations In this section I solve forward modeling for the V C , C A , and P P coil systems. Those coil configurations are mainly used in small loop systems where the coil separation is much bigger than the radius of the source loop, therefore I can use vertical or horizontal magnetic dipole to approximate the source loops. If the secondary magnetic fields for Chapter 2. 17 Recovering 1-D conductivity from the inversion of EM data V C , C A , P P and H C systems are denoted as Hvc, HCA, HPP and Hue, then under the quasi-static assumption, the forward responses are H (T,LO,Z ) ca oU = —- / + ^ 47rr 47T H (r,u;,z ) vc oba = R e-^ -^XJ (Xr)d\ 2h TE Jo 1 R e- -^h J (\r)d\, X{2h (2.11) 2 TE 0 JO — ^ e ^ - ^ ^ A J ^ A r ) ^ , (2.12) 47TT Jo H (r,u;,z ) = —- / HHcir^z*.) = — PP obs 772 fee-^-^^A^^Ar)JA, (2.13) /*oo R e-^ - ^X J {\r)dX, 2h z TE (2.14) 2 0 47T JO where RTE is the T E mode reflection coefficient, m is the moment of the source, z is the observing height. The above results are obtained by removing the primary field from the solutions given by Dey and Ward (1970), Ward and Hohmann (1988), and Ryu et al. (1970). I use intrinsic and input impedance to compute the reflection coefficient RTE'- =F T ! • ^ where Zo is the intrinsic impedance of free space, and input impedance at the first layer, Z , can be found by the recursive procedure outlined by Ryu et al. (1970). I note that the 1 coil configurations do not provide independent information. In fact HCA = HHC — HvcThe usage of this relation will be illustrated later on. A numerical example is given in Figure 2.3. The source and receiver are placed on the surface of a 0.01 S/m half-space. The magnetic susceptibility is set to zero. Panels (a),(b),(c) and (d) in Figure 2.3 show the data from coaxial, perpendicular, vertical coplanar and horizontal coplanar systems. Below 100 hertz, the amplitudes of the real and imaginary components of the coaxial data are about the same. The imaginary Chapter 2. Recovering 1-D conductivity from the inversion of EM data Frequency ' 18 Frequency F i g u r e 2.3: T h e s e c o n d a r y m a g n e t i c field over a half-space o f 0.01 S / m . P a n e l s ( a ) , ( b ) , (c), a n d ( d ) s h o w t h e a b s o l u t e value o f t h e s e c o n d a r y m a g n e t i c f i e l d f r o m t h e c o a x i a l , perpendicular, vertical coplanar a n d horizontal coplanar coil systems. Solid lines denote t h e r e a l c o m p o n e n t s , a n d d a s h e d lines d e n o t e t h e i m a g i n a r y c o m p o n e n t s o f t h e d a t a . c o m p o n e n t s o f t h e d a t a f r o m t h e o t h e r t h r e e c o n f i g u r a t i o n s are l a r g e r t h a n t h e a m p l i t u d e s of t h e r e a l c o m p o n e n t s at l o w frequencies. 2.3 Calculation of the sensitivities I n o r d e r t o c a r r y o u t a r i g o r o u s i n v e r s i o n , one needs t o c o m p u t e t h e s e n s i t i v i t i e s J , w h i c h c o n n e c t t h e p e r t u r b a t i o n s o n t h e m o d e l m t o t h e p e r t u r b a t i o n s o n t h e d a t a D(m). e l e m e n t s o f J are The Chapter 2. Recovering 1-D conductivity from the inversion of EM data ddi(m) ^(111) = - ^ ^ , i = l,2,...,M, orrij j = l,2,...,iV, 19 (2.16) where M and N are the numbers of observations and the model parameters respectively. There exist several different methods for calculating sensitivities. See McGillivray and Oldenburg (1990) for a review of some techniques. A general way of obtaining the Frechet derivative is the adjoint Green's function method, which is used here. 2.3.1 Calculation of the sensitivities for the coplanar system Let the earth be divided into a sequence of layers and let each layer have a constant but unknown conductivity that is to be found. Then <r(z) = jr (2.17) where N is the number of layers, and ipj is the box-car function. Inserting equation (2.17) into equation (2.6) and differentiating the electric field with respect to cr,- yields the sensitivity equation d dE 2 dE . lujfijipjE. 2 Uj—= ~ do-j dz do-j 2 3 (2.18) The sensitivity problem can be solved by using adjoint Green's function solution (Zhang and Oldenburg, 1994): d E ^ Z ^ = i„ T N + 1 Gj(X, z, u)Ej(\, z,u>)dz, (2.19). where E(\,z,u>) is the secondary electric field generated by the source loop in the j t h layer, and G(\,z,u>) is the Green's function which can be obtained by solving (£f-U :)G (\,Z,U ) = 2 j G (\,z ,u,) j j J 8(z-Z o b s = G (\,z u>) j+1 j1 G(\,z,u>) —> 0 when \z\ —> oo, ) (2-20) 20 Chapter 2. Recovering 1-D conductivity from the inversion of EM data The sensitivities for the magnetic field are related to the sensitivities for the electric field through Faraday's law = - J *f> >") dH z O0~j 2.3.2 - tUflj f ^> > h (\r)\>d\. dE Z U 0 JO (2.21) UO~j Calculating sensitivities for other coil systems While it is possible to obtain the sensitivities for other coil systems by using the vector potentials, and repeating the procedure outlined in the previous section, I compute the sensitivities through a simpler way. B y making use of the relationship governing the forward responses from different coil configurations and the sensitivities for the coplanar system, I build the solutions to the sensitivity problems of the other three systems. Equations (2.11) to (2.14) show that the effect of a on the data is exclusively embodied in the reflection coefficient RTE- The difference among the data of those four coil configurations is solely caused by geometry. Let HXY be any one of the secondary magnetic fields given in equations (2.11) to (2.14), and OXY be the geometry factor which can be decided from equations (2.11) to (2.14). Then the generic form of the sensitivities is dH {\oj,z) _ foo dR r XY do-j ^ T E OCj JO 0 X Y { X O J Z ) D K (2 .22) where the index XY denotes any one of VC, PP, CA and HC. The sensitivities of the electrical field for a horizontal coplanar system have been given in equation (2.19). Combining equations (2.22) and (2.19) yields 9RTE - 4 T T _ _ dE(\,u,z) = M2h do-j ^ z) iuifijm daj Substituting equation (2.23) back into equation (2.22) I obtain the sensitivities for the other three systems by applying the appropriate geometric factors: One = ^A J (Ar) - ( ^ -) 2 47T 0 e A 2 2 Chapter 2. Recovering 1-D conductivity from the inversion of EM data O = — OPP = -^A J (Ar)e~ ( ' -^) OCA = OHC — Ovc vc AJ (Ar)e- ( ' -^) A 2 1 1 2 4TT 21 A 2 (2.24) 1 1 One is listed here for completeness. 2.4 The inverse algorithm In nonlinear inverse problems, one is provided with observations d° ,i bs an associated error estimate £j for each datum. = 1,N, and A forward modeling is also supplied so that the ith datum can be written as di = Fi(cr) with the understanding that the magnetic susceptibility and dielectric permittivity are included in the forward mapping. The inverse problem is non-unique. I proceed by introducing a model objective function and then finding that model which minimizes the objective function subject to fitting the data. There are two requirements for the model parameter used in the inversion. Because conductivity is positive, a positivity constraint is required in the inversion. Conductivity can vary several orders in amplitude, so a nonlinear mapping is needed to prevent large values of conductivity from carrying too much weight in the inversion. A n obvious way to achieve those requirements would be the use of m = ln(cr) as the model. The model objective function, <f> = a J w (z)(m - m ) dz + (1 - a) 2 m a 0 d(m — m ) Jw (z) 0 f dz. (2.25) penalizes vertical roughness and differences between the recovered model and a reference model m . In equation (2.25) a is a parameter that controls the relative importance of 0 the two terms. w and Wf are weighting functions which can be prescribed by the user. a When the earth is divided into layers of constant conductivity, equation (2.25) can be 22 Chapter 2. Recovering 1-D conductivity from the inversion of EM data discretized and written as K =|| W {m-m ) m || , (2.26) 2 0 where m = (mi,m , ...mj^) is a model parameter vector and W is an M x M weighting 2 m matrix. I chose a data misfit objective function - D) || = £ ( - ? - ^ i ) 2 , fa =|| W (D° bs (2.27) 2 d i=l Vi where rji is the error associated with the i t h datum, and D oha and D are the observed and predicted data respectively. M y goal is to find a model m that minimizes equation (2.26) subject to the constraint that (f) in equation (2.27) is equal to a target misfit d (j>* . If the errors are Gaussian and independent then (f) is a chi-squared variable and d d its expected value is approximately equal to N for N > 5. Correspondingly, I often set (j>*n ~ N. The optimization problem of minimizing equation (2.26) subject to <f> = <j>* d d requires minimizing ^(m) = ^ + 9 - ( ^ - ^ ) , (2.28) 1 m i where 0 is a Lagrangian multiplier. The nonlinear nature of the problem requires linearization and iteration to a solution. Let m( ) be the model at the nth iteration and let n 8m be a perturbation. The effect of the perturbation on the z'th datum is given through a Taylor's expansion M M o 771 di[m^ + 8m] ~ Fi[m^} + £ -^8m i oj m j=l = d\ n) + E H h J (2-29) 6m j=i where J(j = ddi/dmj is the sensitivity which indicates how d{ is affected by a changing on the model parameter for the j t h layer. The problem becomes: minimize 4> =|| W [8m + rn'"' - m ] | | + / T { | | W {D° 2 m where < ^ n+1 0 1 hs d - F[m^ + 8m}} | | -f 2 }, (n+1) d (2.30) ^ is a target misfit at each iteration. Generally I choose (j)*} ^ = l/*f(j)* ^ n+1 n d where 7 > 1. Writing F[m^ ^] = F[m^] + J8m and setting the gradient Vsm4> = 0, I n+1 Chapter 2. Recovering 1-D conductivity from the inversion of EM data 23 obtain Sm = [f3W^W + fWjWaJ^iJTwJWrfDn m + f3W^W [m m 0 - m W]}, (2.31) where SD is n SD = D -F[mW], obs n (2.32) A nonlinear line search is carried out at each iteration to find that /3 which generates the desired target misfit, or if that is not achievable, to select that /3 which produces the minimum misfit. The inversion procedure continues until the initial desired misfit is achieved and until further iterations produce no significant reduction in the model objective function. The inversion algorithm begins with an initial conductivity model and in each iteration the model is upgraded until the model fits the data and recovers a smooth conductivity structure. The susceptibility distribution has to be specified at the very beginning of the inversion and is fixed throughout the inversion. 2.4.1 Joint inversion of the coplanar and coaxial data In many cases data from different coil configurations are collected at the same time. The most common example is the D I G H E M system, which takes measurements of coplanar and coaxial data simultaneously. The coplanar system has better coupling with horizontal targets while the coaxial configuration couples better with vertical structure. So a joint inversion could delineate the targets better in a 2-D or 3-D environment. However, when 1-D joint inversion algorithms are applied to 3-D data sets, the 3-D effects on data from different coil systems are different. Consequently it may be difficult to recover a model which can explain 3-D data from different coil configurations simultaneously. The model objective function is that given by equation (2.25). But the data objective function includes both coplanar and coaxial data: hi =|| WacA(D$l - DCA) II + || W {Dtp - D ) 2 dCP CP ||, 2 (2.33) Chapter 2. Recovering 1-D conductivity from the inversion of EM data where WdCA and WdCP a r e 24 the weighting matrixes for coaxial and coplanar data, D^p and Dcp are the observed and predicted coplanar data, and D^ and DCA are the observed and predicted coaxial data. The joint inversion can be solved by minimizing the model objective function subject to data constraints. That is, minimize <f>(m) = <t> + /3 (<j>dj ~ 4>* )1 m d (2.34) Setting the gradient V$ (p(m) = 0 yields m 8m = [BWlW + J WjW J]- {J WjW 8D T m x T d d n + (3W^W [m - m^]}, m 0 (2.35) where J = (2.36) 8D = (2.37) and dCA W (2.38) d \ 0 Wdcp ) where JCA and Jcp are the sensitivities for the coplanar and coaxial systems respectively. The value of ft is determined through a non-linear line search in the same manner as in the inversions of individual data sets. 2.5 Synthetic examples In this section I invert coplanar, coaxial, perpendicular and vertical coplanar data generated from both 1-D and 3-D synthetic models. Those synthetic models can be used to simulate sulfide deposits. Susceptibility is assumed to be zero and is fixed through out of the inversion. The data were normalized by the primary fields, and are given in Chapter 2. Recovering 1-D conductivity from the inversion of EM data 25 ppm. The primary field for the perpendicular system is zero, so I use the primary field associated with the coplanar system to normalize the perpendicular data. • 2.5.1 1-D examples In this section I test the algorithm on 1-D synthetic data. E M methods are good in picking up conductive zones but may have difficulty to delineate resistive zones. To make the test more realistic, I designed a 1-D model which consists of two conductive zones sandwiching a resistive layer. Through the test, I want to see whether the inversion can recover the resistive zone. The secondary magnetic fields for all four coil configurations were calculated at 10 frequencies ranging from 110 to 56320 hertz, with each frequency doubling the preceding one. The coil separation is 10 meters and the source and receiver are placed 30 meters above the earth. In performing the following inversions, the earth is divided into 44 layers to a depth of 500 meters. The parameter a is chosen as 0.02. The starting and reference models are 1 0 - 2 S / m half-spaces. The data were corrupted with 1% Gaussian noise. This error level is about what we can expect in field surveys. In the first inversion, I used a percentage error assignment in the inversions. The standard deviations were set to 1% of the amplitudes of the data. It takes 7 to 9 iterations for all four inversions to converge to the desired misfit levels. The results are presented in Figure 2.4. Figures 2.4a, 2.4c, 2.4e, and 2.4g show the recovered models from the inversions of the coplanar, vertical coplanar, perpendicular, and coaxial data. Figures 2.4b, 2.4d, 2.4f, and 2.4h plot the data. A l l inversions successfully recovered the two conductive zones. The resistive layer is also recovered from the inversions, but its position is shifted upwards slightly. It is noted that the conductivity values of the resistive layers recovered from the four inversions are much higher than the true model, and this is a common phenomenon in inverting E M data. The data from all four coil configurations Chapter 2. Recovering 1-D conductivity from the inversion of EM data 26 and at all frequencies are fit equally well, so the recovered models look similar. Different error assignments lead to different results from the inversions. Possible error assignments include a constant value or a certain percentage of the amplitudes of the data. A small constant number may be added to a percentage error scheme to prevent the standard deviations from being zero, and to keep small data, which have small signal-to-noise ratio, from carrying too much weight in the inversions. In the second example, I used the constant value error scheme to invert the same data sets used in the previous example. The standard deviations for the vertical coplanar, horizontal coplanar, and perpendicular data were set to 2 ppm at all frequencies. Since the primary field for the coaxial system is twice as big in amplitude as the primary field for the coplanar and vertical coplanar systems, the standard deviations for the coaxial data were set to 1 ppm at all frequencies, so that the actual errors after being converted back from ppm to A / m are the same for all data. Each of the four inversions converged to the desired misfit levels after 5 to 6 iterations. The recovered conductivity models and the data are presented in Figure 2.5. Panels (a), (c), (e) and (g) are the recovered models from the inversions of the horizontal coplanar, vertical coplanar, perpendicular, and coaxial data respectively. Panels (b), (d), (f) and (h) plot the corresponding predicted and observed data. The recovered model from the inversion of the coplanar data is the best among the four recovered models, followed by those obtained from the inversions of the vertical and coaxial data. The reconstructed conductivity from the inversion of the perpendicular data fails to pick up either the resistive layer or the second conductive zone in depth. This is because, given the error scheme, the perpendicular system has the smallest signal-to-noise ratio, and the signal pertaining to the resistive and the second conductive layers is washed out. By contrast, the signal related to those two layers in the coplanar system is well above the noise level, so the inversion 'of the coplanar data generated the best result. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 27 Figure 2.4: The results of the inversions of the 1-D data, at 10 frequencies, from different coil systems. The standard deviations are set to 1% of the data strength for all data. Panels (a), (c), (e), and (g) show the recovered (solid lines) and true models (dashed lines) from the inversions of the coplanar, vertical coplanar, perpendicular, and coaxial data respectively. Panels (b), (d), (f), and (h) plot the corresponding predicted (lines) and observed data (discrete points) from the inversions, for coplanar, vertical coplanar, perpendicular, and coaxial coil systems. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 28 Figure 2.5: The results of the inversions of the 1-D data from different coil systems. The standard deviations are lppm for the coaxial data, and 2ppm for the data from the other three coil configurations. Panels (a), (c), (e), and (g) show the recovered (solid lines) and true models (dashed lines) from the inversions of the HC, VC, PP, and CA data respectively. Panels (b), (d), (f), and (h) plot the predicted and observed data (discrete dots), for the HC, VC, PP, and CA data. The real and imaginary components of the predicted data are denoted with solid and dashed lines respectively. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 29 I explain this in Figure 2.6, by plotting the true data from previous inversions, along with the data from the same model but without the resistive and the second conductive zones. Panels (a), (b), (c), and (d) show the data for the vertical coplanar, coaxial, perpendicular, and coplanar systems. The coplanar data are the strongest in amplitude, followed by the vertical coplanar, coaxial, and perpendicular data. The information about the resistive and the second conductive layers is contained mainly in the real components of the data at the first three frequencies. By setting the standard deviations to lppm, that useful information in the perpendicular data is buried into the noise completely, therefore the inversion cannot bring out any information about those two layers in depth. The above two experiments indicate that the use of different error schemes generates different results from the inversions. When the signal-to-noise ratios are the same, the inversions of the data from different coil systems can produce the same results. In practice, however, different coil systems do have different signal-to-noise ratios, and hence produce different results. For a 1-D earth, the use of the coplanar system is highly recommended, and the use of the perpendicular system should be discouraged. In airborne E M surveys, data are measured typically at 3 frequencies. In the following example, I investigate the effect of the number of frequencies on the results of inversions. I repeat the inversions in the previous example, but this time only use 3 frequencies at 110, 7040, and 56320 hertz. After 8 to 9 iterations, the inversions converged to the desired chi-square misfit level. Figure 2.7 shows the results. Panels (a), (b), (c), and (d) show the recovered conductivity from the inversions of the coplanar, vertical coplanar, perpendicular, and coaxial data respectively. Compared with the corresponding models in Figure 5, the second conductive zone is shifted upwards slightly but otherwise the results in Figure 2.4 and 2.7 are similar. So it is possible to use fewer frequencies to obtain useful information from the inversions, as long as the range of the frequency is chosen appropriately, to provide the necessary depth of penetration. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 30 F i g u r e 2.6: T h e i n f l u e n c e o n t h e d a t a f r o m t h e r e s i s t i v e a n d t h e s e c o n d c o n d u c t i v e zones i n t h e m o d e l u s e d i n t h e first e x a m p l e . T h e r e a l a n d i m a g i n a r y c o m p o n e n t s o f t h e d a t a u s e d i n t h e first e x a m p l e are d e n o t e d w i t h s o l i d a n d d a s h e d l i n e s , a n d t h e r e a l a n d i m a g i n a r y c o m p o n e n t s o f t h e d a t a g e n e r a t e d f r o m t h e s a m e 1-D m o d e l , b u t w i t h o u t t h e r e s i s t i v e a n d t h e s e c o n d c o n d u c t i v e zones, are d e n o t e d w i t h t r i a n g l e s a n d c i r c l e s . P a n e l s (a), ( b ) , ( c ) , a n d ( d ) present t h e d a t a f r o m t h e v e r t i c a l c o p l a n a r , c o a x i a l , p e r p e n d i c u l a r , and coplanar systems. Chapter 2. 31 Recovering 1-D conductivity from the inversion of EM data I—i 10° i 11 i i nun—i i 10 10 Depth ( m ) H I M — i 1 2 11 it I 10° i i 11 10 10 Depth ( m ) 1 i ii M 2 Figure 2.7: The recovered models from the inversions of the 1-D data at 110, 7040, and 56320 hertz. The solid lines plot the recovered models, and the dashed lines denote the true model. The error is assumed to be 1% of the amplitude of the data. Panels (a), (b), (c), and (d) show the models from the inversions of the coplanar, vertical coplanar, perpendicular, and coaxial data. Chapter 2. 2.5.2 Recovering 1-D conductivity from the inversion of EM data 32 3-D examples Geological targets are in general 3-D, therefore it is desired to use 3-D algorithms in interpreting the data. However, due to the limitations on computer resources, 3-D solvers are still too expensive to use. As a result, it is quite often to apply 1-D algorithms on 3-D data sets. Hence it is necessary to investigate how well the 1-D algorithms can perform on 3-D models, and how reliable the images obtained from 1-D inversions of 3-D data sets are. In the next example, I invert 3-D data obtained from the model shown in Figure 2.8. The data were generated by Newman and Alumbaugh (1995) with a staggered finite differences method, and contaminated with 1% Gaussian noise. In order to examine the 3-D effect on the 1-D inversions, I start with this relatively simple model that consists of two conductive prisms sitting in a resistive half-space of 0.01 S/m. Those two prisms are separated from each other by 100m, so they can interact with each other yet the inversion could still recover separated anomaly bodies. The conductivity of the shallower prism is 0.1 S/m, and for the deeper prism it is 0.5 S/m. The susceptibility takes its value in free space. The frequencies used in this experiment are 880, 7040, and 56320 hertz. The coil separation and observation height are all the same as in the previous 1-D examples. The line spacing is 50 meters and the station interval is 25 meters. The standard deviations for the coplanar data were set to 5 ppm plus 10% of the data. For the coaxial, perpendicular, and vertical coplanar data, the standard deviations were set to 2 ppm plus 10% of the data strength. The error was determined by comparing 1-D forward modelings with the 3-D data. It was found that the discrepancies in those two data sets are about 5-10% over the prisms. In the inversions the data were all fit to the desired misfit level. Other parameters in the inversion were kept the same as in the 1-D inversions. In displaying the recovered models and the data, the positions of the two Chapter 2. Recovering 1-D conductivity from the inversion of EM data 33 prisms are denoted by white rectangles. Example 1 I first invert the vertical coplanar data. The real and imaginary components of the predicted and observed data from the inversion of the vertical coplanar data are presented in Figures 2.9 and 2.10. The data clearly indicate the presence of two anomalous bodies, but it is impossible to determine whether those two bodies are conductive or resistive, let alone the depths of burial. The imaginary component of the data forms two positive peaks at 880 hertz, one peak and one sink at 7040 hertz, and two toughs at 56320 hertz, over the tops of the two prisms. The real component of the data, on the other hand, shows two peaks and no toughs over the tops of the two prisms at all three frequencies. The y-z cross-sections of the recovered model at x=250, 350, and 450 meters are shown in Figure 2.11. The cross-sections at y=250, 350, and 450 meters are plotted in Figure 2.12. Figure 2.13 shows the recovered conductivity at z=15, 20, 30, 40, 50, and 60 meters depth. The recovered model is a good representation of the true model, even though it is a little wider and thinner than the true model. The existence of the two prisms is shown quite clearly in the reconstructed model. Because of the good coupling between the source and the vertical boundaries of the model in the y-z plan, the recovered model starts to see the deeper prism at x=250m. The coaxial and coplanar data look similar to the vertical coplanar data. The perpendicular data, however, seem to be more prone to 3-D effects, and are much rougher than those data from the other three coil systems. Figure 2.14 shows the real component of the perpendicular data. The data from coplanar, coaxial, and perpendicular coil configurations were inverted with the same parameters as in the previous inversion of the vertical coplanar data. The cross-sections at y=350m, of the recovered models from the inversions of the coaxial, coplanar, and perpendicular data, are presented in Figure 2.15. The Chapter 2. Recovering 1-D conductivity from the inversion of EM data 34 Y(m) Z(m) X(m) 700 500 0 1—r" 100 500 700 15 GO Ko 30 ? N ' 300 50 on 100 | 0 300 Ko 100 100 300 500 700 X(m) X-Y Plan-view X - Z Cross Section Figure 2.8: The 3-D model used to generate the synthetic data. The background conductivity cr = 0.01 S/m. The conductivity in the upper and lower prisms is 0.1 and 0.5 S/m respectively. Magnetic susceptibility is set to zero. Source and receivers are placed at 30 meters above the earth. The coil separation is 10 meters. Frequencies used in the forward modeling are 880, 7040, and 56320 hertz. The block volume which holds the two conductors are denoted with dashed white lines, and the white numbers are the coordinates of the left low corner of the block volume. 0 Chapter 2. Recovering 1-D conductivity from the inversion of EM data 35 O b s e r v e d data P r e d i c t e d data 600 600 400 400 H 200 200 //IBI, o f=880 f=880 —i 200 C 1 1 400 —I 200 1 1— 600 600 600 400 400 200 200 1 1— 400 600 400 600 i f o f=7040 0 200 _ 600 1 1 1— 600 i 400 400 H 200 200 f=56k f=56k - 0 1 1 200 400 Easting (m) 600 0 1 1 200 1 1 400 1 1— 600 Easting (m) F i g u r e 2.9: T h e r e a l c o m p o n e n t o f t h e v e r t i c a l c o p l a n a r d a t a i n p p m . Panels (a), (b), a n d ( c ) s h o w t h e p r e d i c t e d d a t a , a n d panels ( d ) , (e), a n d (f) show t h e o b s e r v e d d a t a . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 36 F i g u r e 2.10: T h e i m a g i n a r y c o m p o n e n t o f t h e p r e d i c t e d a n d o b s e r v e d v e r t i c a l c o p l a n a r data i n p p m . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 37 i 50 JZ < — - • 100 a m 150 Conductivity (S/m) • i 200 100 E J= Q. CD Q • "*-~^__ 50 > a o • 400 5CC 600 700 300 400 500 600 700 300 400 500 600 700 300 \ 100 150 Conductivity (S/m) | 100 EE a 200 50 H 100 150 Conductivity (S/m) 0 100 200 Easting (m) F i g u r e 2 . 1 1 : T h e x - z cross-sections o f t h e r e c o n s t r u c t e d c o n d u c t i v i t y at ( a ) y = 2 5 0 m , ( b ) y = 3 5 0 m , a n d (c) y = 4 5 0 m , from t h e inversion of the vertical coplanar data. Chapter 2. Recovering 1-D conductivity from the inversion of EM data S 38 50 i 100 150 Conductivity (S/m) 100 200 300 400 500 600 700 300 400 500 600 700 I I 50 100 H 150 Conductivity (S/m) 100 200 ! - - - . - y) 50 JZ 100 - CL 0 150 - [cj Conductivity (S/m) i 100 r 200 300 400 500 600 700 Northing (m) F i g u r e 2.12: T h e r e c o v e r e d m o d e l f r o m t h e i n v e r s i o n o f t h e v e r t i c a l c o p l a n a r d a t a , at (a) x = 2 5 0 m , ( b ) x = 3 5 0 m , a n d (c) z = 3 5 0 m . Chapter 2. 39 Recovering 1-D conductivity from the inversion of EM data 0.195 600 • 200 400 ( 1 ! 0.100 0.071 d v 400 0.139 600 o 0.051 0.036 i 200 0.019 0.013 z=40m a z=1 5m 0.026 H 0.010 0 200 400 I I 0 600 200 400 I ! l 600 I 0.139 600 600 0.195 0.100 0.071 400 400 H 0.051 0.036 0.026 200 200 0.019 0.013 z=50 z=20j —i 200 1 400 0.010 1— 600 400 200 600 0.195 600 - 600 - 400 - 400 - 0.139 0.100 0.071 200 - fi W 0.036 1 1 J 0 •i 200 0.026 1 200 : - 0.019 if] z=60m z=30m j o 0.051 i 400 i 600 Easting (m) 0 J 1 200 1 400 ! 600 0.013 0.010 Easting (m) F i g u r e 2.13: T h e x - y p l a n - v i e w s o f t h e r e c o v e r e d c o n d u c t i v i t y , f r o m t h e i n v e r s i o n o f t h e v e r t i c a l c o p l a n a r d a t a , at s i x different d e p t h s . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 40 recovered model from the inversion of the coaxial data has the most vertical structure, and beyond the vertical boundaries of the deeper prism, it bends downwards slightly. The model recovered from the inversion of the coplanar data is more localized and does not bend downward. The recovered model from the inversion of the perpendicular data is not a good representation of the true model. The roughness in the perpendicular data is caused by 3-D effects, and cannot be fit by 1-D inversions. Consequently, the recovered conductivity also shows strong roughness. The 1-D inversions of 3-D data in the above examples have generated useful information about the true model. In all recovered models, the existence of the two anomalous bodies is evident. Since I used relatively large standard deviations in order not to overfit the 3-D data, the amplitudes of the recovered conductivities are lower than the true value. Because of the good coupling between the source and the vertical boundaries, the recovered model from the inversion of the coaxial data has the sharpest vertical structures in the direction orthogonal to its source plane. The recovered conductivity from the inversion of the vertical coplanar data has the longest extension in the direction perpendicular to its source plane. From the 3-D example, it is found that the perpendicular coil system is most sensitive to 3-D effects, and consequently the recovered model, bear the most roughness. But at the same time, the model extracted from the perpendicular data shows two completely separated conductive bodies. A common feature in the inversions of all four coil systems is that, due to 3-D effects, the recovered models are thinner and wider in horizontal directions than the true models. So caution has to be exercised when interpreting the results from 1-D inversions of 3-D data, especially for the perpendicular data. Chapter 2. Recovering 1-D conductivity from the inversion of EM data O b s e r v e d data Predicted data 600 - 600 g 41 400 400 o 200 - 200 0 600 600 400 400 200 200 o 200 400 • el | f=7040 0 200 400 600 600 0 200 400 600 0 200 400 600 600 400 200 200 400 Easting (m) 600 Easting (m) F i g u r e 2.14: T h e r e a l c o m p o n e n t o f t h e p e r p e n d i c u l a r d a t a i n p p m . P a n e l s ( a ) , ( b ) , a n d (c) s h o w t h e p r e d i c t e d d a t a , a n d panels ( d ) , ( e ) , a n d (f) s h o w t h e o b s e r v e d d a t a . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 42 50 100 CA 150 0 1 i 1 100 200 300 1 400 r~ : 500 600 500 600 500 600 700 Easting (m) 0 100 200 300 400 Easting (m) 50 pp 0 100 200 300 400 700 Easting (m) F i g u r e 2.15: T h e c o m p a r i s o n o f r e c o v e r e d m o d e l s f r o m t h e i n v e r s i o n s o f t h e c o p l a n a r , p e r p e n d i c u l a r , a n d c o a x i a l d a t a at s e c t i o n y = 3 5 0 m . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 43 The linear relationship among the data As I mentioned previously, the relationship R = HCA — HHC + Hvc holds only over a 1-D earth. R should be zero over 1-D earth, and becomes non-zero when there exist 3-D effects. I plot R at all three frequencies over sections at y=250, 350, and 450 meters in Figures 2.16, 2.17, and 2.18 respectively. The horizontal positions of those two anomalous bodies are marked with grey bars. In all three cases, strong variations of R occur over the two conductive prisms. A t 56k hertz, the real component of R forms two local peaks at the two vertical boundaries of the anomalous targets, whereas the imaginary component of R generates troughs. At 7040 hertz, the real component of R creates local maximums over the top of the bodies and toughs over the edges of the conductive prisms. The imaginary component of R forms a minimum over the top of the deeper prism. Both components of R at 880 hertz register local maximums over the top of the deeper prism, even though the amplitude of the imaginary component is much greater than the real component's. So potentially I can use this relationship to locate the horizontal positions of underground geological targets. Tests on synthetic and field data are needed to make this technique a practical mapping tool. Example 2 In E M surveys, conductive overburdens reduce the depth of investigation, and make it more difficult to locate targets below the overburdens. As a last synthetic example, the effect of the overburden on the inversion is investigated. The inclusion of the conductive overburden brings the experiment closer to reality, especially for places like Australia, where E M surveys have been plagued by conductive overburden. In this test, an overburden with variable conductivity and thickness is added to the 3-D model shown in Figure 2.8. The overburden is separated into two units by a vertical fault at x=350m. The Chapter 2. Recovering 1-D conductivity from the inversion of EM data 44 F i g u r e 2.16: T h e r e a l a n d i m a g i n a r y c o m p o n e n t s o f R over t h e s e c t i o n at y = 2 5 0 m , at 880, 7040, a n d 56320 h e r t z . D a s h e d lines d e n o t e t h e i m a g i n a r y c o m p o n e n t , a n d t h e s o l i d lines d e n o t e t h e r e a l c o m p o n e n t o f R. T h e g r e y b a r denotes t h e h o r i z o n t a l p o s i t i o n o f the conductive p r i s m . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 45 F i g u r e 2.17: T h e r e a l ( t h e s o l i d l i n e ) a n d i m a g i n a r y ( t h e d a s h e d l i n e ) c o m p o n e n t s o f R at y = 3 5 0 m . T h e h o r i z o n t a l p o s i t i o n s o f t h e t w o c o n d u c t i v e p r i s m s a r e m a r k e d w i t h t h e grey bars. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 46 Figure 2.18: The profiles of R at y=450m. The solid line denotes the real, and the dashed line denotes the imaginary components. The grey bar indicates the horizontal position of the deeper conductive prism. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 47 conductivity and the thickness for the unit extending from x=0 to 350m are 0.02 and 30m, respectively. For the other unit, the conductivity is 0.1 S/m and the thickness is 25m. The background conductivity is lowered to 0.001 S/m from the original 0.01 S/m, so the conductivity contrast is enlarged to 100. The line spacing is increased from 50m to 100m, and the station interval is increased from 25m to 50m. Other parameters of the model remain the same. The data were contaminated with 1% Gaussian noise. I first invert the coaxial data at section y=400m. After 5 to 9 iterations at all stations, the chi-square misfit was reduced from the initial 50415 to 118. The standard deviation used in the inversion is lppm plus 2% of the data strength. Figure 2.19 shows the model and the data. The data clearly indicate the variation on conductivity of the overburden. The skin depth over a 0.1 S/m half-space is only about 53m at 880 hertz, and 19m at 7040 hertz, therefore only the data at 880 hertz may penetrate the overburden. The real component of the data at 880 hertz (panel (c)) forms a local peak over the top of the deeper prism. Because of the good coupling between the transmitter and the vertical edges of the prism, Re(Hx) drops sharply at x=400m and x=500m. That implies that the recovered model could define the horizontal extension of the deeper prism accurately. The reason why Im(Hx) in panel (c) does not see the deeper prism is that the imaginary component of the data decays faster than the real component. Panel (d) shows the recovered model. Because of the limited skin depth, the inversion only recovered the top part of the deeper prism. The shallower prism is picked up by the inversion, too. However, since the survey line passes the edge of the upper prism, the signal level related to this prism is low. Thus the amplitude of the recovered model is smaller than that of the true model. The position of the vertical fault is located accurately. The thickness and the conductivity of the overburden were recovered reasonably well. The horizontal coplanar data are inverted next. The error scheme used in the inversion is also 2ppm plus 2% of the data strength. The cumulative initial chi-square misfit over Chapter 2. 48 Recovering 1-D conductivity from the inversion of EM data F i g u r e 2.19: T h e r e c o v e r e d m o d e l ( y = 4 0 0 m ) f r o m t h e i n v e r s i o n o f t h e c o a x i a l d a t a . T h e s o l i d a n d d a s h e d lines p l o t t h e r e a l a n d i m a g i n a r y c o m p o n e n t s o f t h e p r e d i c t e d data. T h e o b s e r v e d d a t a are p l o t t e d w i t h d i s c r e t e dots. T h e t w o p r i s m s a n d t h e o v e r b u r d e n s are d e n o t e d w i t h w h i t e b o x e s . H s a n d H p a r e t h e s e c o n d a r y a n d t h e p r i m a r y fields. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 49 the whole section is 57469. After 4 to 6 iterations at each station, the cumulative misfit was reduced to the desired misfit level 90. Figure 2.20 shows the model and the data from the inversion. Similar to the coaxial data, the data at all frequencies indicate the variation on the conductivity of the overburden, and only the real component of the data at 880 hertz clearly indicates the existence of the lower prism. However, the peak of Re(Hz) over the vertical edges of the deeper prism changes more gradually, and that implies that the recovered model will have more smooth variation in the x-direction. The recovered model is presented in panel (d). The conductivity values for the overburden were accurately recovered. Vertical variations inside the overburden is reduced, and the horizontal extension of the prisms is blurred, compared to the recovered model from the inversion of the coaxial data. I also inverted the vertical coplanar data and horizontal coplanar data over section x=450m. The standard deviations are 2ppm plus 2% data strength. The accumulative starting and final Chi-square misfits for the vertical coplanar data over all 6 stations are 31825 and 53 respectively. For the horizontal coplanar data, the initial and final misfits are 35458 and 35. The recovered models are plotted in the same figure (Figure 2.21) for comparison. Because of its orientation, the vertical coplanar system has the best coupling with the vertical boundaries in the y-direction of the two prisms. The recovered model from the vertical coplanar data delineates the vertical boundaries of the prism much better, but the recovered conductivity overburden from the inversion of horizontal coplanar data is more vertically uniform. Joint inversion of the coaxial and coplanar data Since sources in different orientations couple with the 3-D geological targets differently, the inversions of the data from different coil systems can provide independent information. Therefore a 3-D joint inversion of those data sets is certainly beneficial. However, in 1-D Chapter 2. Recovering 1-D conductivity from the inversion of EM data 7500 E Q. tx LO 50 56320 Hz 6000 -i 1 1 JL _ * 1 1 1 4500 H 3000 x 1500 H 3500 3000 m E Q. j 7040 Hz] 2500 /*•*" • a_ _ ^- — ^ — —» - —• - - C L 2000 ... -r-S X b 1400 E 1200 CL C L 1000 CL X "35 X 880 Hz 800 600 400 200 70C 0.25100 0.18100 0.13100 0.09420 0.06790 0.04900 0.03530 0.02550 0.01840 0.01320 0.00955 F i g u r e 2.20: T h e r e c o v e r e d m o d e l , at y = 4 0 0 m , f r o m t h e i n v e r s i o n o f t h e h o r i z o n t a l coplanar data. T h e s o l i d a n d d a s h e d lines p l o t t h e r e a l a n d i m a g i n a r y c o m p o n e n t s o f the predicted data. T h e observed d a t a are p l o t t e d w i t h discrete dots. T h e two prisms a n d t h e o v e r b u r d e n s a r e d e n o t e d w i t h w h i t e b o x e s . H s a n d H p are t h e s e c o n d a r y a n d the p r i m a r y fields. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 51 50 100 150 H VC S/m 2C0 100 200 200 400 500 600 700 Y(m) 50 1 100 150 HC S/m 200 o 100 200 300 400 500 600 700 Y(m) F i g u r e 2.21: T h e recovered models a n d the data from the inversions of t h e vertical c o p l a n a r a n d h o r i z o n t a l c o p l a n a r d a t a at x = 4 5 0 m . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 52 inversions of 3-D data, the 3-D effects may not affect data in different coil configurations in the same manner. So a 1-D joint inversion of 3-D data may not necessarily generate better results than the inversions of individual data sets. It is of interest, therefore, to investigate whether a 1-D joint inversion can be helpful in interpreting 3-D data. Since the coaxial and coplanar data having been collected routinely in airborne E M surveys, I invert the coaxial and coplanar data in example 2 jointly. Figure 2.22 shows the model and the coaxial data from the joint inversion of the coaxial and coplanar data. Compared to the model recovered from the inversion of the coplanar data, the lateral extent of the model is better defined. The 0.1 S/m conductive overburden recovered from the joint inversion is more uniform than that recovered from the inversion of the coaxial data. The predicted and observed coplanar data from the joint inversion are shown in Figure 2.23. The standard deviations for both data sets were the same as in the separate inversions. Since the 3-D structure affects the coaxial and coplanar system differently, 1-D joint inversion usually cannot fit the data to the same degree as in the inversion of a single data set. The total starting Chi-square misfit over all 15 stations is 107883, and after 5 to 7 iterations, it was reduced to 407. In the above examples, all inversions indicate the existence of the lower prism. However, because of the conductive overburden, the depth of investigation is greatly reduced. For a 0.1 S/m half-space the skin depth is only 53m, therefore only the top part of the lower prism is recovered. Since the survey line at y=400 passes the edge of the upper prism, the signal level about this prism is relatively low. Further more, the overburden further blurs the signal needed to resolve the upper prism. As a result the inversions did not define the extension of the upper conductor clearly. The inversion of the H C data reproduces the 0.1 S/m conductive overburden the most accurately. The inversion of the C A data delineates the vertical fault and the vertical boundaries in the x-direction of the deeper prism the best, while the inversion of the V C data delineates the vertical Chapter 2. Recovering 1-D conductivity from the inversion of EM data 53 2000 E" 1600 - \ 56320 Hz Q . -9= 1200 C L X CO X E 800 400 - GL GOO CL 450 CL X «5 X 7040 Hz 750 300 I _ _ J> m 150 300 CL Q. CL 240 i80 X 120 X 60 c/j 880 Hz H 50 100 CL CD Q 150 Id! S/m 200 i oo 200 300 400 500 600 700 0.25100 0.18100 0.13100 0.09420 0.06790 0.04900 0.03530 0.02550 0.01840 0.01320 0.00955 X(m) F i g u r e 2.22: T h e r e c o v e r e d m o d e l f r o m t h e j o i n t i n v e r s i o n o f t h e coarxial a n d t h e h o r i z o n t a l c o p l a n a r d a t a at y = 4 0 0 m . T h e s o l i d a n d d a s h e d lines p l o t t h e r e a l a n d i m a g i n a r y components of the predicted coaxial data. T h e observed d a t a are p l o t t e d w i t h discrete d o t s . T h e t w o p r i s m s a n d t h e o v e r b u r d e n s a r e d e n o t e d w i t h w h i t e b o x e s . H s a n d H p are t h e s e c o n d a r y a n d t h e p r i m a r y fields. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 54 X(m) F i g u r e 2.23: T h e o b s e r v e d a n d t h e p r e d i c t e d c o p l a n a r d a t a f r o m t h e j o i n t i n v e r s i o n o f t h e c o a x i a l a n d c o p l a n a r d a t a . T h e s o l i d a n d d a s h e d lines d e n o t e t h e r e a l a n d i m a g i n a r y components of t h e predicted data. T h e observed d a t a are p l o t t e d w i t h discrete dots. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 55 boundaries of the lower prism in the y-direction the best. The joint inversion generated a model which exhibits characteristics of the models recovered from the inversions of the individual data sets. 2.6 Field data examples In the following example, I invert field data collected over the section Y9600 at M t . MiUigan which is Cu-Au porphyry deposit located in north central British Columbia. A D I G H E M airborne E M system was flown and inphase and quadrature phase data from the coaxial coils were collected at 900 hertz. Inphase and quadrature phase coplanar data at 900, 7200 and 56000 hertz were measured at the same time. The coil separation is 6.3 meters for the coplanar data at 56,000 hertz, and 8 meters for coaxial data and coplanar data at 900 and 7200 hertz. The inversion was carried out by parameterizing the earth into 44 layers. The model objective function given in equation (2.25) was minimized in the inversion. The parameter a is set to 0.02. The reference model was a conductive half-space of 1.7 mS/m, which was determined from the inversion of D C resistivity data over the same section. Because of the existence of strong magnetization in this region, the inphase data for both coaxial and coplanar systems at 900 hertz are negative at most of the stations. Since I do not know the susceptibility, I set /x = fio in the inversion as it is commonly done. Thus it is impossible to fit those negative inphase data, and consequently only the quadrature data were inverted. The error assigned to the coaxial data was 0.02 ppm plus 5% of the strength of the data. For the coplanar data, the error was 0.2 ppm plus 10% of the data. The reason for me to use a smaller error for the coaxial data is that there is only one coaxial datum at each station so I want to extract as much information as possible from it, even though I realize that I should not over fit the data when 3-D effects exist. After 8 iterations at each station, the total chi-square misfit level Chapter 2. Recovering 1-D conductivity from the inversion of EM data 56 at all the 13 stations was reduced to 14.6 from the accumulative initial misfit of 73756 for the coaxial data, and to 56.9 from 14777 for the coplanar data. The recovered models from the inversions of both coaxial and coplanar data at each station are combined to form the 2-D sections shown in Figure 2.6. Panel (a) shows the recovered conductivity structure from the inversion of the coaxial data. Panel (b) shows the recovered model from the inversion of the coplanar data over the same section. A joint inversion was carried out on the coaxial and coplanar data. The error used in the joint inversion was 1 ppm plus 10% of the strength of the data, and only the quadrature data were inverted because of the negative inphase data. The accumulative chi-squared misfit was reduced from 10888 to 1770 after 10 iterations at all stations. The result is shown in panel (c). Since 3-D structure affects C A and H C coil configuration in different ways, these two data sets may be incompatible with the 1-D assumption and the ascribed error. The recovered model is a compromise of the two individual inversions shown in panels (a) and (b). The two major conductive zones are recovered by the joint inversion, and vertical structure is increased, compared to the model in panel (b). Panel (d) shows the recovered conductivity model from the inversion of 2-D D C resistivity data (Oldenburg et al., 1995). The recovered model indicates two major conductive zones at both ends of the sections. These two conductive zones correspond to the conductivity model from the inversion of 2-D D C resistivity data. 2.7 Conclusion In this Chapter I have solved the inverse problem of the reconstruction of conductivity from the inversions of the data from H C , C A , P P and V C systems. A Gauss-Newton algorithm was developed which minimises a norm of the model subject to the data constraints. The model norm can be constructed to allow the inversion to find the most 57 Chapter 2. Recovering 1-D conductivity from the inversion of EM data 12.4 12.6 12.8 13.0 13.2 13.4 50 •g ^ § Q 100 150 200 250 12.4 12.6 12.8 13.0 13.2 13.4 63.096 30.824 15.058 7.356 3.594 1.756 0.858 0.419 0.205 0.100 63.096 30.824 15.058 7.356 3.594 1.756 0.858 0.419 0.205 0.100 63.096 30.824 15.058 7.356 3.594 1.756 0.858 0.419 0.205 0.100 63.096 30.824 15.058 7.356 3.594 1.756 0.858 0.419 0.205 0.100 Easting (km) Figure 2.24: The recovered conductivity models over section Y9600 at M t . Milligan. (a) The recovered conductivity model from the coaxial data; (b) recovered model from the coplanar data; (c) recovered model from joint inversion of the coaxial and coplanar data; and (d) reconstructed model from the inversion of 2-D D C resistivity data. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 58 plausible model from the infinite number of models that can reproduce the data adequately. The choice of the final model depends upon the geological setting and the prior knowledge. It could be the one which has the minimum structure, or the one with the smallest roughness, or the one that is closest to some reference model representing the preconceived image of the region under study. The following procedures have been taken to make the inversion stable and robust. First the misfit at each iteration is reduced gradually. Second the minimum norm and the flattest model norm solutions are sought. The final model recovered by this inversion algorithm is insensitive to the particular starting models used. The solution to the 1-D inverse problem described in this Chapter is not restricted by computer resources. So the algorithm developed in this Chapter can be used to invert larger amounts of data. But when this 1-D algorithm is applied to 3-D data, the recovered models will inevitably be affected by 3-D effects, and therefore caution is needed in interpreting the results. A linear relationship is established among the data from the vertical coplanar, coaxial, and horizontal coplanar systems, which can be used to check the validation of 1-D assumption, and can be potentially used as a structure mapping tool, as demonstrated in the synthetic example. The 1-D examples show that the horizontal coplanar data has the highest signal-tonoise ratio, followed by vertical coplanar and coaxial data. The perpendicular system has the lowest signal-to-noise ratio. This indicates order of preference for use in 1-D analysis. Different error assignments lead to different results in the inversions, but usually we do not know the errors. When a percentage error is used in the inversions of 1-D data from different coil systems, the data are fit to the same degree and hence the recovered models are the same. Some systems apply constant value for all data. With this fixed error the coplanar data might generate the best results, and the inversion of the perpendicular Chapter 2. Recovering 1-D conductivity from the inversion of EM data 59 data might lead to the worst result because of the low signal-to-noise ratio. From the particular example studied, the perpendicular data is most sensitive to 3-D effects. Therefore when 1-D inversion is carried out on those 3-D perpendicular data, the recovered model may be unreliable. But on the other hand, the 3-D inversion of those 3-D perpendicular data may delineate the targets better. This is yet to be verified when 3-D algorithms become available. The use of 1-D joint inversion in 3-D environment improves the quality of the recovered models obtained from the inversions of individual data sets, as shown in the synthetic and field data examples. In the 1-D joint inversion of 3-D data, the data cannot be fit to the same degree as in the inversions of individual data sets, because the 3-D effects may not affect data from different coil systems in the same manner. Whether the 1-D joint inversion of 3-D data is useful is model dependent, and difficulty may occur when the 3-D effects are too strong. Chapter 3 1-D Inversions of the Surface and Borehole Transient E M data 3.1 Introduction The course of development of E M methods applied in drill holes has been treated in some detail in Dyck (1975). E M methods were first applied in the well-logging industry. Doll (1949) described the basis for interpretation of induction logs as the "geometrical factor" concept in which interaction among eddy currents induced in the medium is neglected. The eddy currents are considered to be concentric with the drill hole for the simplest configuration in which both source and receiver are close together (of the order of l m ) and lie coaxial along the hole. The use of this rather large-diameter device has not widely spread to mining exploration. One of the large-scale prospecting systems designed to detect conducting sulfide bodies at some distance from the hole was introduced in Noakes (1951). The prototype system was similar to the surface Turam method in that the transmitter was a large, fixed loop on the surface and the receiver was a dual-coil device that measured the difference in the magnetic field at two points 15 meters apart. Concurrently, Ward and Harvey (1954) developed a tilt-angle method of down-hole E M surveying in which the configuration was minimum coupled. A third type of large-scale system employs the fixed-separation, coaxial-coil configuration akin to the induction logger. Elliot (1961) described such a system which operated at 1230 hertz and measured inphase and quadrature components of secondary field with a sensitivity of 100 ppm. Versions operated 60 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 61 at larger coil separation up to 150 m were subsequently developed but the only presently operating system is the somewhat smaller scale equipment, with coil separations up to 20m. At present, large-loop systems with a single component magnetic field receiver in either time-domain ( T E M ) form or in frequency domain ( F E M ) version are the downhole E M methods enjoying the greatest popularity in field application. These systems are made up by a large loop laid on the ground and a receiver down the borehole which can measure 3 components of the magnetic field. The geometric configuration is common to all three types: impulse-type, step-type, and multi-frequency systems. The loop is of a size comparable with the depth of the borehole. The transmitter loop is usually located at several different positions to get optimum coupling between the target and the system in high dimensional environment. Commercially available systems involve a single component induction-coil sensor, coaxial to the borehole, and linked to surface by a cable. A l l systems can be used on both surface and drill hole configurations. The borehole measurement device discussed in this section consists of two coils, one a transmitter and the other a receiver. The transmitter coil is energized with an alternating current at frequencies usually ranging from 10 hertz to several tens of kilohertz, or a series waveform, while the excited field is measured at the receiver. Application of drill-hole E M techniques used in exploring for bodies of conductive mineralization can detect and define the conductive target itself or other geological features which may lead indirectly to discovery of a mineral deposit. Drill-hole applications can be defined to include any configuration where either the source or the receiver (or both) is submersed in the lower half-space. Relatively large-scale prospecting systems which could be considered extensions of surface methods are emphasized here. B y comparison, inductive logging techniques employed in hydrocarbon exploration and evaluation have Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 62 seen relatively little use in mineral prospecting. The distinction is more than just a matter of scale: the latter devices, surrounded by the medium to which they are sensitive, yield a record of physical properties of the rocks in the immediate vicinity of the hole, i.e., an in situ rock-property measurement. E M prospecting systems, on the other hand, are sensitive to conductive bodies lying at some distance from the hole. Borehole E M techniques have also found wide applications in search for groundwater. As an example of such application, groundwater exploration at M t . Morgan palaeochannel aquifer, which is located about 40 km southwest of Kambalda in the Goldfields region of Australia, was carried out by Western Mining Corporation Limited ( W M C ) in 1991. This aquifer is being developed as a source of saline process water for the St. Lves Gold Mill. Surface and borehole transient E M data have been collected in this region. The inversions of surface T E M data have generated some information, but it is expected that the inclusion of the borehole T E M data in the inversion would considerably enhance the images of the conductivity distribution. However, the inversion of T E M borehole data has not been fully addressed, and more research needed. In this Chapter I investigate whether, and how, the joint inversion of borehole and surface T E M data improve the recovered models from the inversion of the surface T E M data. A n inverse algorithm is developed to recover conductivity from the inversions of transient E M data from a large loop system. The forward problem is solved in the frequencydomain. A two-way recursive procedure from both the surface and the bottom of the model is used to calculate the input impedance at each layer, and then the secondary magnetic field is expressed as a function of the input impedance. The frequency-domain data are then converted to time-domain data. A n adjoint Green's function solution is used to compute the sensitivities in the frequency-domain, which are then converted into the time-domain. The inverse problem is solved in the time-domain by minimizing the objective function subject to the data constraints. I invert synthetic and field data to Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 63 test the inversion algorithm. S I R O T E M data collected in the groundwater exploration from W M C is inverted, and the 3-D effects on the inversion are discussed. 3.2 Data conversion Time domain data can be inverted in two ways. One method is to transform data from the time-domain to the frequency-domain, and then invert the converted data in the frequency-domain. The other way is to keep the data in the time-domain. The forward modeling and calculation of Frechet derivatives are carried out in the frequency-domain, and then converted into the time-domain. The inversion is therefore carried out in the time-domain. With either of these two approaches I need to convert a time-domain function to the frequency-domain and vice versa. The advantage of solving the inverse problem in the time-domain is that the error in the data is clearly defined. The disadvantages include that the data have to be calculated at a wide range of frequencies and that there are fewer existing inversion codes in the time-domain. On the other hand, solving inverse problems in the frequency-domain can make use of many existing inversion software. However, the nature of the error in the converted data become unclear, and more study is needed to investigate how the original error in the time-domain data is altered during the transform. More important, the narrow time spans of the field data disallow accurate transforms. I address these important aspects in this Chapter. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 3.2.1 64 Conversion of frequency-domain data to time-domain data The Fourier transformation of the frequency response for a causal step turn-off can be calculated using a sine or cosine transform (Newman et al., 1986): dh(t)_2 2 t°° r „ , , ^ . , . ; , _ 2j r°° — / Im[H(u>)} sinwtdu! = / dt 7T Jo TV Jo T f f or = / — 1 v n zosLotdu = h(0) - - LO 7T ./0 / TT Re[H(u>)] cosu/tdto, — —^ smut(ku, l Jo J1 (3-1) (3.2) LO where H(u>) and h(t) are the magnetic fields in the frequency- and time-domains respectively. Numerical evaluation of these sine-cosine transforms is carried out with Anderson's (1979) digital filters. Even though either the real or imaginary component of the magnetic field can be used in the conversion, better results have been obtained by using Im[iJz(a;)] (Newman et al., 1986). Thus the imaginary component of the frequency-domain data is preferred in the conversion. 3.2.2 Conversion of transient data to frequency-domain data From the theory of sine and cosine transformations, I may obtain the magnetic field in the frequency-domain from equations (3.1) and (3.2) by means of inverse sine or cosine transformation, provided that the time-domain response is known. The real and imaginary components of the magnetic field are given by 1 Re[H(u>)} = — 27T r°° Jo dh(t) — — cos iotdt = -to dt r°° Jo [h{t) - h(0)] sin wtdt, (3.3) and Im[H(w)] = f°° sin utdt = -to f°° h(t) cos totdt. Jo dt Jo Equations (3.3) and (3.4) can be cast into a standard form: (3.4) Chapter 3. 65 1-D Inversions of the Surface and Borehole Transient EM data poo F(c) = / Jo f(t) sin wt (3.5) dt, cosa>r. which can be easily calculated by using Anderson's digital filter. One of the major difficulties here is that the narrow time spans of the field data disallow accurate transforms. I illustrate this in the following numerical example. I calculate the responses from a central loop system in both time- and frequencydomains over a half-space of 0.01 S/m. The current excitation in the frequency-domain is e lwt and in the time-domain it is a step turn-off with unit amplitude. Analytic solutions exist for this simple problem (Ward and Hohmann, 1988, p220-223). The transmitter is a circular loop with a radius of 50 meters, and it is placed on the surface of the earth. First I convert the data from the frequency- to time-domain. The frequency band used in the conversion is from 10~ to 10 Hertz. Figure 3.1 shows the conversion. That 3 6 the converted data agree with the true data well is attributed to the wide range of the frequency in the forward modeling. Next I convert the calculated time derivative of the magnetic field from the timedomain to the frequency-domain. A time span from 10~ to 1 0 6 _1 seconds is used in the conversion. Figure 3.2 shows the result. The transformed frequency-domain data agree with the theoretical data well between 10 and 10 Hertz. At frequencies lower 8 than 10 Hertz, the converted data become very unreliable. This means that for the frequency range usually used in inversions, a time span of 10~ to 1 0 6 _1 seconds can provide satisfactory frequency data. For most T E M systems time spans are around the range of 1 0 - 4 to 1 0 _1 seconds, and such a narrow span will lead to incorrect results in converting the data from the time- to frequency-domain. Figure 3.3 shows such an example. The standard S I R O T E M time span (4.87 x 10~ to 0.161 seconds) was used in the conversion. The converted 4 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 1 Q - 8 I I HII 10- I 6 I I I I I 111 10~ s I I I I I MM 10-" I ' I I 10~ I I 3 in | i i i ii 10" 66 ni 2 Time(s) Figure 3.1: The conversion of FEM data to T E M data due to a step turn-off, for a central loop system over a half-space with 0.01 S/m conductivity. The transmitter is a 50m loop in radius. The input FEM data have a frequency span from 10 to 10 Hertz. The solid line and the dashed line denote the true magnetic field and its time derivative respectively. The discrete symbols denote the converted transient data from the data in the frequency-domain. -3 6 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 67 Figure 3.2: The conversion of step turn-off T E M data to FEM data for a central loop system over a 0.01 S/m half-space. The source loop has a radius of 50m. The input T E M data have a time span from 10~ to 10 seconds. The solid line and the stars denote the true FEM data and the converted FEM data respectively. 6 _1 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data data deviate badly from the true data. 68 The error in the conversion is caused by the lack of the early-time data. This is because the data at earlier time channels are several orders of magnitude greater than the late-time data, and hence they carry much of the weight in the evaluation of equation (3.5). It is possible to improve the conversion by extrapolation. But the appropriate choice of extrapolation functions needs yet to be investigated. The truncation of data after 1 0 _1 second does not affect the conversion in a perceivable way as has been confirmed by numerical tests. In order to obtain accurate data within the frequency band used in inversions, the time-domain data measured as early as 10~ seconds must be used in the conversion. 6 In this paper, I calculate the forward responses and sensitivities in the frequencydomain, and then convert them into the time-domain. The inversion is carried out in the time-domain. 3.3 Forward modeling The T E M system considered in this Chapter consists of a large loop source on the surface and a receiver which is separated from the source by a radial distance r. The receiver can take measurements on either the surface or in a borehole. Forward responses for a central loop system can be obtained by setting r = 0. As I will show later, the solution for a loop source inside a layered half-space is needed in the calculation of the sensitivities with the adjoint Green's function method. So to be general, I consider the forward modeling problem for a loop source located within a 1-D earth. The geometry of the system is shown in Figure 3.4. A loop of radius a carrying current Ie lwt is placed inside an N-layer half-space at a depth of z . The receiver is placed at a a radial distance r from the source loop. The circular symmetry allows me to employ a cylindrical coordinate system, and solve the problem in the Hankel transform domain. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 69 10-' ^ <S 10- 2 10~ 3 10~ 4 10" 5 icr i_ 6 10- 7 ioio- 8 9 1 0 - 1 0 io- ' io-' 1 2 10- 13 10 - 3 10 - 2 10 _ 1 10° 10 1 10 2 10 3 10 4 10 5 10 8 10 10 7 8 Frequency Figure 3.3: The conversion of T E M data to F E M data for a central loop system whose loop source has a radius of 50m. The impulse response of the magnetic field which spans 10~ to 10 seconds is used as the input data. The solid line denotes the imaginary component of the true frequency-domain response, and the stars denote the converted frequency-domain response. 10 -1 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 70 Figure 3.4: The geometry for the inversion of borehole E M data. A source loop of radius a energized by a current Ie is placed at z in the borehole, and the receiver takes the measurements at a depth of z and a radial distance of r. %ujt s The governing equation for the electric field is given by dz 2 Ej(\,z,u>) = iujfijaI(u>)Ji(\a)S(z — z ), s (3-6) where A is the Hankel transform parameter, fij is the susceptibility in the jth. layer, and Uj = y A + u> fij6 — iu>Hj<jj. 2 2 0 For convenience, I define layers above and below the source-containing layer as regions A and C respectively, and define the source-containing layer as region B . I first derive the solutions for source-free regions A and C, and then discuss the solution for region B. 71 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 3.3.1 Solutions for the source-free regions In regions A and C, the source term on the right hand side of equation (3.6) disappears, and the general solution for this homogeneous equation is the combination of up- and down-going waves. The boundary conditions are the continuities of the tangential components of the electric and magnetic fields. I designed a two-way recursive procedure to match the boundary conditions. In region C , the intrinsic impedance is defined as Z i (\ t U ) = (3.7) Uj The input impedance for any layer in region C is given by Zj + Z^ 3 ianh(ujhj) 1 ' y The secondary electric field in region C is given by Ej = A , - e ' ' + Bje- > - i \ (3.9) B ^ B ^ e ^ - ^ ^ + z / ^ (3-10) u (z_z '' +l) U (z Z +1 where and 7J+ — Z1 Aj = Bje Uidi 3 3 Zi +1 7 + Zj - (3.11) In region A the intrinsic and input impedances are Zj(X,co) = -^, 1 (3.12) and z i + i = ^ z 3 Z j + ^tanli(uA-) + Z3 \an\\(ujhj) ( v 3 . 1 3 ' ) Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 72 The solution for region A is Ej = Aje^"^ + B - ^ -^\ u (3.14) z je where the coefficients are given by and where A* 3.3.2 — A e~ i+ i+ , u j+1 lh 1 j+1 and h - z j+x - Zj. j+1 Solution in the source-containing layer In the source layer, care has to be taken because of the singularity existing at z = z . s When z is not equal to z , the solution for equation (3.6) is a f A E (r,z,co) = { u z j e + B~ z<z uz je s j ( Cje + D uz j e ' u z z> (3.17) z s where the coefficients Aj,Bj, Cj and Dj can be determined in the following manner. A t z = z , E is continuous and its first order derivative is discontinuous s E j \ z t dEj dz = + '_ = E i \ z 7 ' (3.18) iu>fiaJi(\a). From the definition of the input impedance at the j t h and (j + l ) t h interfaces, I obtain (3.19) Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 73 and C = D - 2 ^ u e + — 1 '-. (3.20) Solving equation (3.18) along with equations (3.19) and (3.20), I obtain A=—e —T 33 2 anc , e- > iR R -l 2u 1 C UJZ, 2 ; _ ~ , e- ^ >R R 2 (3.21) V UjZj-ujhj I Z).= — e"'^- 2 ' h - 1' h 1 2 (3.22) V ' where Z ' — Z3 Ri = 4Z>- + r 4Z , (3-23) 3 and ZJ+ - Z1 R 2 fL = fi. (3.24) . (3.25) S is given by S = ^ a J l ( A a ) Once the amplitude of the electric field in the source-containing layer is determined, the electric field in regions A and C can be solved recursively, by using equations (3.7) to (3.16). If both z and z are set to zero, the above results are reduced to the response from a a surface configuration, and the two-way recursive procedure is reduced to the one-way recursive procedure outlined by Ryu et al. (1970). In numerical implementation of the forward modeling, the tanh(a;) function is expressed with exponential functions with negative arguments: x _ -x i _ -2x tanh(z) = — — ^ = — — . e + e 1 -f- e x x i (3.26) x Doing so can prevent the possible numerical overflow from happening. The secondary magnetic field in the frequency-domain is obtained by using Faraday's law: 74 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data H (r,u>,z) z -1 /•oo (3.27) ILOU.Q JO The final solution to the forward problem is obtained by converting the frequency response into the time-domain. As an example, frequency and transient responses are calculated over a half-space of 0.01 S/m. The receiver is placed alone the axis of the borehole at 50m depth. The radius of the transmitter is 50 meters. The left panel in Figure 3.5 shows the frequency response of the total magnetic field, and the right panel shows the transient response. Only the imaginary component of the frequency-domain data is plotted because the conversion of the data from the frequency- to time-domain only involves Im.(H Z(LO)) while at low frequencies, the responses measured at both the surface and in the borehole are almost identical. A t high frequencies, the borehole response lower because the high frequency component of the signal cannot reach the receiver at depth due to the skin effect. The corresponding transient data were shown in the right panel of Figure 3.5. To better understanding the behaviour of the transient data, I borrow the concept of "smoke ring" proposed by Nabighian (1979). According to his work the combined effect of all induced currents in the ground can be approximated by the effect of a single current filament of the same shape of the transmitter loop, and which moves downward and outward. The current intensity of the equivalent current filament varies inversely proportional with time. The transient magnetic field measured at the surface is similar to that measured in the borehole. The most visible difference of the two data sets occurs near 10~ s, where 5 the amplitude of the borehole data is slightly larger than that of the surface data. This overshoot occurs because, when the "smoke ring" passes through the plane of the receiver at depth, the receiver records a stronger signal than what it would read on the surface. In Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data Frequency Time 75 (s) Figure 3.5: Frequency and transient responses over a half-space of 0.01 S/m for a central loop borehole system. The radius of the transmitter is 50m, and the observing depth is 50 meters. The left panel shows the frequency response. The solid line denotes the borehole response, and the dashed line denotes the response on the surface. The right panel plots the transient response in the borehole due to a step turn-off. The lines denote the responses measured on the surface of the half-space. The discrete symbols denote the borehole responses. fact for this particular example the amplitude of the borehole data becomes greater than that of the surface data after 3 x 10~ s. However, the surface data overpower the borehole 6 data before 3 x 10~ s which indicates that the surface data carry more information about 6 the near surface structure. Therefore the surface data has higher signal-to-noise ratio at early times while the borehole data has higher signal-to-noise ratio at late times. The above characteristics of the magnetic field are even more evident in its time-derivative. Here 10 -5 s serves as a watershed with the borehole data exceeding surface data in amplitude after this time and the surface data having a higher signal level before this time. In the next example, I investigate how the borehole data vary as a function of observation depth by calculating the forward responses over a 0.1 S/m half-space, for receiver positions from 0 to 100m, and for a time span from 10 -5 to 10 s. The radius of the _2 76 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 1e-05 0.0001 Time (s) 0.001 1e-05 0.0001 0.001 Time (s) Figure 3.6: The transient responses, as a function of the receiver position and time, due to a step turn-off over a half-space of 0.1 S/m for a central loop borehole system. The radius of the transmitter is 50m, and the data are calculated at 34 locations between 0 and 100m. Panel (a) shows the magnetic field, and panel (b) plots the time derivative of the magnetic field. source loop is 50m. In Figure 3.6a I plot the amplitude of the magnetic field, and Figure 3.6b shows the amplitude of the time derivative of the magnetic field. For a given time channel, there exists an observation depth at which the strength of the data reaches the maximum value. That depth corresponds to the vertical distance between the surface and the plane of the "smoke ring" at that moment of time. When the observing depth is beyond that optimum depth, the amplitude of the borehole data drops, and it eventually becomes weaker than that of the surface data when the distance between the observation point and the plane of the "smoke ring" is larger than that from the surface to the plane of the "smoke ring". Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 3.3.3 77 Calculation of the sensitivities The sensitivity problem can be solved by essentially the same technique used in Chapter 2. Here however, the adjoint Green's function has to be calculated for receiver positions inside the borehole. The sensitivities for the electrical field are given by J ZJ O0~j where Ej(\, z,u>) is the secondary electric field generated by the source loop in the jth layer, and G is the Green's function which can be obtained by solving (& - f) u G (X z u) s 1 j1 GJCM,<") = ( 6 z - ob,) Z = G (\ z ,u) j+1 1 (3-29) j G(X,z,u>) —» 0 when \z\ —> oo, where z^ is the observation depth in the borehole. Note that G satisfies the same partial a differential equations as the electric field in equation (3.6) except a scaling factor on the right hand side. Thus the Green's function can be solved in the same manner as outlined in equations (3.6) to (3.25). The sensitivities for the magnetic field are given by ' - -"> = -±-r ag r z a ^ ' " ' " ) J . ( A r ) A ^ . As a numerical example, I calculate the sensitivities for the same model used in the example in Figure 3.5, for receiver locations at both surface and 50m in the borehole. Figures 3.7a and 3.7b show the real and imaginary components of the sensitivities in the frequency-domain, from the layer between 48 to 50 meters. The sensitivities for the borehole configuration exceed those for the surface configuration over all the frequencies in amplitude, and decay slower. Panels (c) and (d) in Figure 3.7 present the sensitivities for the magnetic field and its time-derivative due to a step turn-off. The sensitivities in (3.30) Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 78 the time domain for the borehole configuration are also larger than those for the surface configuration. Zero-crossings occur in the sensitivities at early time in the time-domain. To investigate how the sensitivities vary with respect to depth, in Figure 3.81 plot the sensitivities of the imaginary component of the magnetic field and the impulse response of the magnetic field from the same experiment, as a function of the time and depth. The left panel shows the imaginary component of the sensitivities in the frequency-domain. At a given frequency, the amplitude of the sensitivities of the magnetic field reaches a maximum value and then decreases as the depth further increases, and sign changes occur at about 70, 140, and 330m. A t a fixed depth, the amplitude of the sensitivities also reaches a maximum value at the optimum frequency, and decreases as the frequency increases further. The right panel shows the sensitivities for the impulse response. For a given depth, the amplitude of the sensitivities reaches a maximum at an optimum time channel. This optimum time occurs earlier for layers near the surface and later for deeper structures. This is consistent with the fact that the equivalent current filament travels from the surface to depth with increasing radius, therefore at early times the field is sensitive to the near surface structures and at late time sensitive to deep structures. 3.4 Synthetic data example The data from a central loop configuration over a synthetic model were calculated at 32 time channels, ranging from 0.487 ms to 161.4 ms, for receiver positions at 0, 150, 200, 250, and 300m. The model was designed based upon the resistivity logging in a ground water exploration in Australia, and consists of a conductive overburden and three conductive layers at depth, and are shown in Figure 3.9a (the dashed line). Through this example I want to examine whether the use of borehole data in the inversion can delineate the two conductive layers centered at 250 and 320m better. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data icr 10- 6 b to \ N X <o _ a 7 icr 10" 8 \ 9 -10 - 10" ' - 1 0 1 10~ 101 0 / \ \ \ \ \ 1 2 \ \ 1 3 -14 10~ 1 1 5 10 10° _ 1 10 1 10 1 1 2 Frequency 10 10* 3 10 S 1C 10- (Hertz) 10 eg _ N 10° - 10 10 10- C/) 10- E 10~ 2 10- 3 10" 4 10" 5 1 | -/ / / \ ~\l 2 CD \ N c CD 1 10" \ -C £NJ 10- 10" 10" 1 4 10- Time ( s ) 1 3 10° 10 1 10 2 10" 10" io10" 10" 10 10 3 4 10 S 10 6 (Hertz) - v / \ 1 10" 10" 1 Frequency 10 10 E 79 \ 6 7 S 1 -6 0 10~* 10" 3 Time (s) Figure 3.7: The sensitivities for both surface and borehole central loop configurations in the time- and frequency domains. These sensitivities are due to a layer with a thickness of 2m at 48 meters in a half-space of 0.01 S/m. The radius of the transmitter is 50m, and the observing depth for the borehole configuration is 50 meters below the ground. The solid lines denote the sensitivities for the borehole system, while the dashed lines denote the sensitivities for the surface configuration. Panel (a) plots the real components of the sensitivities in frequency-domain, and panel (b) plots the imaginary components of the sensitivities. Panels (c) and (d) show the sensitivities in time-domain, for the magnetic field and its time derivative. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 80 Figure 3.8: The sensitivities for borehole central loop configurations in the time- and frequency-domains in a half-space of 0.01 S / m . The radius of the transmitter is 50m, and the observing depth for the borehole configuration is 50 meters below the surface. The left panel plots the imaginary component of the sensitivities in the frequency-domain, and the right panel plots the sensitivities of the impulse response. Logarithmic scale is used to display the data. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 81 The data are denoted with triangles in Figure 3.9b, and are corrupted with 1% Gaussian noise. The source has an area of 900 square meters, and it is driven by a step turn off with unit amplitude. The same inversion procedure used in Chapter 2 is used to invert the data. The model objective function is the one given in equation (2.25). Only the flattest model is used here, so I need not to use a half-space reference model. This is a useful procedure when in practice the reference models are unknown. The earth is modeled as 101 layers up to a depth of 10000 meters. The starting models for the following synthetic examples are all 0.1 S/m half-spaces, and the standard deviations were set to 1% of the data strength. To save computational time a linear line search is carried out at each iteration to find that 8 which generates the desired target misfit. Such a procedure is valid as long as the perturbation on the model at each iteration is not too excessive. Two measures have been taken to ensure such an approximation valid. The first one is to use a small value for parameter 7 in reducing the misfit level at each iteration. It is found that 7 = 2 is a good choice for most cases, and it is not recommended to use any value greater than 7 = 5. The second measure is to control the size of the model perturbation on the model at each iteration. If the perturbation at the (n + l)th iteration for any layer exceeds 50% of rrS \ then the line search is terminated, and the last n 8 value is adopted. The inversion algorithm begins with an initial conductivity model and at each iteration the model is upgraded until the initial desired misfit is achieved and further iterations produce no significant reduction in the model objective function. I first invert the data acquired at the surface. After 8 iterations, the data were fit to the desired misfit level. Figures 3.9a and 3.9b show the models and data. The recovered model is a smoothed version of the true model. The inversion recovered the conductive overburden and the first resistive zone, but missed detail structures at depth. Next I inverted the data measured at 150m in the borehole. It took 9 iterations for Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 1 0 1 100 ' 200 Depth 0 100 200 Depth 1 300 1 400 1 500 I—I 1 I I 1 H H I O - (m) 300 (m) I I IIIH Time iOO 500 IO" 1 1 1 0 " 3 ••' ! 1 0 82 • - 1 (s) 10"* Time (s) Figure 3.9: The models and the data from the inversions of the surface and borehole T E M data. Panels (a), (c), (e), and (g) present the true (dashed line) and recovered models (solid line) from the inversion of the data measured at 0, 150, 250, and 300m. Panels (b), (d), (f), and (h) plot the corresponding observed (triangles) and predicted (solid line) data. The observation positions are indicated by arrows. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 83 the inversion to converge to the desired misfit level. Figures 3.9c and 3.9d show the results. The quality of the image of the conductive layers at depth has been dramatically improved. Now the existence of the two resistive zones and three conductive layers becomes evident. The price for such an improvement at depth is the degradation of the model near the surface. The borehole data are smaller at early time, and larger at late time in amplitude than the surface data. Figure 3.9e shows the recovered conductivity from the inversion of the borehole data at 250m depth. After 14 iterations, the data misfit level is reduced to the desired level. The inversion recovered the true amplitudes of the second and the third conductive layers. Delineation of the conductor near 280m is increased. However, the resolution for the conductive layer around 180m and the resistive layer around 190m is degraded. The predicted and observed data are plotted in Figure 3.9f. The amplitude of the data before 1.8 ms becomes smaller than that of the data measured at 150m. But after 1.8 ms, the amplitude of the data obtained at 250m is larger than that of the data measured at 150m. The above two examples indicate that the inversion of borehole data can provide more information about the structures at depth. However, for a given geometry and the time span in an experiment, the data contain information about structures only to a certain depth. If the observing depth is too great, then the signal-to-noise ratio will be too low, and hence the information generated from the inversions of those data is reduced. To explain this I invert the data measured at 300m depth. Figures 3.9g and 3.9h show the model and data. The inversion recovered a model with a single peak after 24 iterations. This is because the data do not have the high frequency components required to resolve detailed structure. Those high frequency components have been lost during the propagation of the E M field through the earth. The amplitude of the data is smaller at all time channels than that of the data at 250m. Thus the merits of the Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 84 Figure 3.10: The results from the inversion of the central loop data from the borehole configuration. The receiver is placed at 200m along the axis of the borehole. The solid and dashed lines in the left panel denote the recovered and true models, respectively. The right panel plots the observed (the triangles) and the predicted data (the solid line). The arrow indicates the observing point. borehole system no longer hold. One of the difficulties in the inversions of E M data is to pick up resistive targets in a conductive host. In the following example I investigate the possibility of improving the resolution for those resistive layers by taking measurements inside them. Figure 3.10 shows the models and the data from the inversion of the data observed at 200m, within the second resistive layer. The misfit was reduced to the desired level after 16 iterations. The recovered model is similar to that from the inversion of the data measured at the surface, except that the later has better resolution for the overburden and the first resistive layer. Compared to the recovered models from the inversion of data measured within conductive layers at 150 and 250m depth, this is not a good result. So it seems that, at least in this particular example, borehole data measured within resistive zones may not be necessarily helpful in delineating the structure at depth. The inversions of the surface data can resolve near surface structure better, while the inversions of the borehole data can delineate the structure at depth better. So joint Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 85 Figure 3.11: The recovered model and the data from the joint inversion of the central loop data from both surface and.borehole configuration. The receiver for the borehole configuration is placed at 150m along the axis of the borehole. The solid and dashed lines in the left panel denote the recovered and the true models, respectively. The right panel plots the observed and predicted data. The solid line denotes the predicted surface data, while the dashed line plots the predicted borehole data. The observed data are denoted with discrete points. The two arrows indicate the receiver positions. inversions of data measured at both surface and in the borehole is potentially beneficial. In the next two examples joint inversions are carried out. In these two inversions, all inverse parameters and the assigned error are kept the same as in previous inversions of individual data sets. I first invert the data measured at the surface and at 150m in the borehole. It took 11 iterations for the inversion to converge to the desired misfit. Figure 3.11 shows the recovered model and the data. The recovered model is a hybrid of the models obtained from the separate inversions of the surface and borehole data. It recovered a blurred conductive zone between 200 and 350m, and the existence of the two resistive zones is not obvious in the recovered model. To enhance the delineation of the resistive zones around 190 and 280m, I jointly invert the data measured at the surface, 150, and 250m. The chi-square misfit is reduced from 1933498 to the desired level of 96 after 14 iterations. Figure 3.12 shows the result. The Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 86 Figure 3.12: The recovered model and the data from the joint inversion of the central loop data from both surface and borehole configurations. The receiver for the borehole configuration is placed at 150m and 250m along the axis of the borehole. The solid and dashed lines in the left panel denote the recovered and the true models, respectively. The right panel plots the observed and predicted data. The solid line (curve 1) denotes the predicted surface data, while the dashed lines plots the predicted borehole data (curve 2 for 150m, and curve 3 for 230m). The observed data are denoted with discrete points. The three arrows indicate the receiver positions joint inversion recovered all conductive and resistive structures of the true model. The position of the second conductive layer around 250m is slightly shifted towards greater depth. The above two examples show that joint inversions of surface and borehole data can enhance the quality of the recovered model. Next I look at the application of this inverse technique to field data. 3.5 Field data example In this section I invert T E M data from the groundwater exploration mentioned earlier in the introduction of this Chapter. Transient E M data have been collected in this region with the S I R O T E M system. The transmitter is a square loop of 300m by 300m, and the 87 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data data were measured along the axis of the source loop, on the surface and in the borehole. The effective area of the receiver is 10 m . The time span used in this survey is from 4 2 0.487ms to 67.3ms, and the data are recorded at 26 time-channels. The borehole is 300m in depth, and its top 150 meters was cased with steel. The casing can certainly affect the data, and invalidated the 1-D assumption. The hole was logged below the casing and resistivity measurements from 155m to 195m were obtained. In carrying out the inversions, the earth is modeled with 101 layers to a depth of 5000m. A circular loop with radius of 169.3m is used to approximate the square loop. The radius of the loop is determined in such a way that the circular loop will have the same area as the 300m by 300m square loop used in the survey. For a loop of this size, this approximation is accurate enough for the inversions. The starting model is a 0.1 S/m half-space. The error is assumed to be 5% of the data strength plus 1 0 - 5 A / m . The constant component in the assigned error is intended to prevent late-time data, whose signal-to-noise ratio is very low, from carrying too much weight in the inversion. I first invert the surface data. After 8 iterations, the chi-square misfit level was reduced from 7190 to the desired misfit level. The recovered conductivity, along with the observed and predicted data, is shown in Figure 3.13. The apparent conductivity obtained from the down-hole resistivity logging is also plotted in the same Figure. The recovered conductivity picks up the conductive zone around 150 meters, but the resistive zone on the borehole conductivity at 270 meters does not show up in the recovered model. The data measured at 155m are inverted next. The inversion converged to the desired misfit level after 9 iterations. Figure 3.14 shows the results of the inversion. The recovered model shows a resistive zone sandwiched between two conductive layers at depth. But the position of the resistive zone is shifted upwards slightly, compared to the borehole conductivity obtained from the resistivity logging. The conductive overburden in Figure 3.13 does not show up in this recovered model. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 88 10' V) 1 E < ^—* 10° 10-' I T 10-2 CO CO \ JZ N CN CD 100 200 300 400 500 lO" 3 10"* 10" 3 Depth ( m ) 10~ 2 10-' Time (s) Figure 3.13: The restilts from the inversion of surface S I R O T E M data. The left panel plots the recovered conductivity (the solid line) and the borehole conductivity (the dashed line) obtained from down-hole resistivity logging. The right panel shows the predicted data (the continuous line) and the observed data (the discrete points). The arrow indicates the receiver position. 0 100 200 300 Depth ( m ) 400 500 10" 3 10" 2 10"' Time (s) Figure 3.14: The results from the inversion of borehole S I R O T E M data. The receiver is at 155m in the borehole. The left panel plots the recovered conductivity (the solid line) and the borehole conductivity (the dashed line) obtained from down-hole resistivity logging. The right panel shows the predicted data (the continuous line) and the observed data (the triangles). The arrow indicates the receiver position. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data , , Depth ( m ) 10' | 89 , Time (s) Figure 3.15: The recovered model and the observed and predicted data from the inversion of the S I R O T E M data collected at 230m in the borehole. The left panel plots the recovered conductivity (the solid line) and the borehole conductivity (the dashed line) obtained from down-hole resistivity logging. The right panel shows the predicted data (the continuous line) and the observed data (the triangles). The arrow indicates the receiver position. Data collected at 230m are also inverted. After 11 iterations, the chi-square misfit was reduced from the initial value of 10908 to the target level of 26. Figure 3.15 presents the models and the data from the inversion. The recovered model shows a very conductive zone around 250m, which contradicts the results from the resistivity logging. This makes the inversion suspect. Possible causes of this disagreement include the steel casing and other non 1-D effects. To enhance the information about structure at depth, a joint inversion is carried out to invert the surface data and borehole data at 155m depth. Knowing that the borehole data have better signal-to-noise ratio at late times and that the surface data have better signal level at early times, I assigned about 10% error to the borehole data before 3 ms, and 2% afterwards. For the surface data, the standard deviations were 1% before 20 ms, and 5% afterwards. After 20 iterations, the initial chi-square misfit was reduced from 7449 to 253. Figure 3.16 shows the recovered model and the data. Considerable extra Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 90 Figure 3.16: The recovered model and the observed and predicted data from the joint inversion of the borehole (at 155m) and the surface S I R O T E M data. The left panel shows the recovered model (the solid line) and the resistivity logging (the dashed line). The right panel plots the observed (discrete points) and the predicted data (the solid line). The receiver positions are indicated by arrows. structure, especially the conductor at 90m, has been recovered. The second conductive zone at 200m seems to have some connection with the .conductive zone in the resistivity logging but it is shifted towards shallower depth. The third conductive zone in the recovered model is consistent with the resistivity logging and is much improved from the recovered models from the inversions of individual data sets in Figures 3.13 and 3.14. The first conductive zone does not show up in the models recovered from the inversions of the individual data sets. The data at 230m were also inverted, but the inversion could not converge to the desired misfit level. , Possible causes responsible for the unsatisfactory results from the 1-D inversion of the field data are the casing and other non 1-D effects, which invalidate the 1-D assumption and make the surface and borehole data incompatible. To investigate the compatibility of the surface and borehole data, I first calculate the response at z=230m generated from the model recovered from the inversion of the surface data, and compare this forward Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data • i i i i 11 I I icr 3 i i ii i 1111 ioTime (s) 2 i i i i nil '" I i ' 10-' 10- 3 ' 1 ' 1 10-2 Time ( s ) * * i i 91 1111 io- 1 Figure 3.17: Compatability of the surface and borehole data at 230m depth. The solid lines denote the calculated data, while the triangles plot the field data. The left panel compares the borehole S I R O T E M data at 230m depth with the calculated borehole data of the model recovered from the inversion of the surface S I R O T E M data. The right panel shows the surface S I R O T E M data and the forward surface responses of the model recovered from the inversion of the borehole S I R O T E M data at 230m depth. response to the observed data at the same depth. Figure 3.17a shows the comparison. Those two data sets are vastly different from each other before I O - 2 seconds. This means that it is difficult to find a 1-D model which can reproduce the two field data sets at the same time. I also calculate the forward response at the surface from the model recovered from the inversion of the borehole data collected at 230m, and plot the calculated data in Figure 3.17b against the field data measured at the surface. These two data sets are also not even close at early time channels, therefore a 1-D joint inversion could not find a model which can explain these two data sets. i I also did the same analysis to the data measured at 155m. Figure 3.18 shows the comparison of the calculated and the field data. The left panel shows the predicted data at 155m from the model recovered from the inversion of the surface data, and the measured data at 155m. At early times the predicted data overshot the observed data. This is because the recovered model from the surface data is more conductive near the surface Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data "ICT 3 1CT 2 Time (s) 10 10" 3 10~ 2 92 io- Time (s) Figure 3.18: Compatability of the surface and borehole data at 155m depth. The solid lines denote the calculated data, while the triangles plot the field data. The left panel compares the borehole S I R O T E M data at 155m depth with the calculated borehole data of the model recovered from the inversion of the surface S I R O T E M data. The right panel shows the surface S I R O T E M data and the forward surface responses of the model recovered from the inversion of the borehole S I R O T E M data at 155m depth. than the recovered model from the borehole data. The right panel shows the field data measured at the surface and the predicted data at the surface from the model recovered from the inversion of the borehole data at 155m. The difference between these two data sets is also mainly at the early times. The overall difference between the borehole and surface data at 155m is much smaller than that for the data at 230m. Field data results are inconclusive but illustrate the importance of non 1-D data. The E M data are sensitive to a volume conductivity structure around the borehole, while the resistivity logging is sensitive to a much smaller region. The surveyed region is relatively 1-D but dipping structures do exist according to geological information. The non 1-D effects cause incompatibility between the surface and borehole data. Thus while it is possible to explain the surface data or the borehole data separately with distorted 1-D images of the earth, it is difficult to find a 1-D model which can reproduce the surface and the borehole data simultaneously. Chapter 3. 3.6 1-D Inversions of the Surface and Borehole Transient EM data 93 Conclusions The inversion algorithm developed in this Chapter provides a new tool to obtain information about the conductivity distribution over a 1-D earth from the inversion of transient E M data collected at both the surface and in the borehole. The amplitude of the borehole data is smaller than that of the surface data at early times. A t later time, that relation is reversed, so the surface data have better signal-tonoise ratio at early time channels while the borehole data enjoy higher signal-to-noise ratio at late times. The model recovered from the inversion of the surface data represents the near surface structure better, while the model obtained from the inversion of borehole data delineate the structure at depth better. Joint inversions of the data measured at the surface and in the borehole can improve the images obtained from the inversions of individual data sets. The field data example is inconclusive but shows the importance of non 1-D effects. Since geological structures are usually 3-D, it necessary to test this 1-D algorithm on 3-D synthetic data sets, to investigate how the 3-D effects alter the inversions. It is recognized at the outset that 3-D effects might present problems for the 1-D inversions of borehole T E M data. Nonetheless there are areas where this algorithm might be applicable. Chapter 4 Recovering susceptibility from 1-D Inversion of E M data 4.1 Introduction Electromagnetic data are sensitive to conductivity <r, magnetic permeability /x and to electrical permittivity e. The decay of the E M fields in the earth depends upon these parameters and the frequency of the source. It follows that a multi-frequency sounding contains information about all of these properties as a function of depth and, in principle, it is possible to simultaneously invert any set of data for <r, fi and e. The first step in this process of simultaneous inversion is to be able to invert for any one of the parameters when the others are specified. Electrical conductivity characteristically varies by orders of magnitude and there have been numerous papers devoted to recovering a when fi and e have been specified. Relative electrical permittivity varies from about 1 to 80 (Keller, 1990) and this parameter has received much attention for surveys carried out at high frequency. Although relative magnetic permeability may vary from 1 to 20 for various rocks and minerals, in practice it varies from 1 to less than 2.0. The effect of such small variations can often be ignored and E M data are commonly inverted after employing the assumption that fi = /io, the relative permeability of free space. Nevertheless, there are instances where permeability changes alter the data in a significant way. A well known example of this is the negative inphase data measured with a typical frequency domain airborne electromagnetic ( A E M ) system. Those negative data cannot be the response of a purely conductive model. Magnetic permeability greater than fia, or equivalently 94 Chapter 4. Recovering susceptibihty, from 1-D Inversion of EM data 95 positive magnetic susceptibility, must exist. Magnetic susceptibility is an important physical parameter in geophysical surveys, but the usual way to obtain information about the distribution of susceptibihty is through the inversion of static magnetic data obtained from magnetic surveys. Unfortunately these data can be reproduced by a layer of susceptible material at the earth's surface. This illustrates not only the extreme non-uniqueness inherent in the interpretation of magnetic data but also shows that there is no inherent information about the susceptibihty distribution with depth. Algorithms which obtain depth distributions do so by imposing parameterization on the model domain (Bhattacharyya, 1980; Zeyen and Pous, 1991; Wang and Hansen, 1990), by applying constraints to the solution (Last and Kubik, 1983; Guillen and Menichetti, 1984), or by introducing a depth weighting function to counteract the natural decay of the kernel functions. A n example of this last approach is given in L i and Oldenburg (1996). Rigorous inversion of E M data to estimate the distribution of magnetic susceptibihty over an arbitrary 1-D, 2-D or 3-D earth has not yet been fully investigated. Work has been done to estimate the physical and geometric parameters of some simple models. Ward (1959) described a method of determining the ratio of magnetic susceptibihty of a conducting magnetic sphere to the susceptibihty of the background rock. He used a uniform field for frequencies which span a large range and encompass the critical frequency at which the frequency independent magnetic field cancels the inphase component due to induced current. Fraser (1973) proposed a way to estimate the amount of magnetite contained in a vertical dike under the assumption that the body is nonconductive. He also developed (Fraser, 1981) a magnetite mapping technique for the horizontal coplanar coils of a closely coupled multi-coil airborne E M system. That technique yields contours of apparent weight percent magnetite under the assumption that the earth is a homogeneous half space. In this Chapter I attack the inverse problem of the reconstruction of susceptibility Chapter 4. 96 Recovering susceptibihty from 1-D Inversion of EM data by assuming that e = eo and that <r is variable, but known. I restrict myself to the 1-D problem. Since the forward modeling has been addressed in Chapter 2, I begin with the inversion procedure and derive expressions for the sensitivities. Synthetic and field data are then inverted and I present summarizing comments in a concluding section. 4.2 Inversion algorithm The same technique as in Chapter 2 is used to solve the inverse problem. Since most rocks have positive susceptibility, positivity constraint is needed. Positivity of the solution can be guaranteed in a number of ways. The simplest is to choose m = Inn as the model. Since Sinn = SK/K the sensitivities are easily obtained for this parameter. The difficulty with this mapping is that near-zero values carry too much weight in the model objective function, and large values of susceptibihty are over estimated due to the nature of the logarithm function. To overcome these difficulties I define m = f(n) as a three-piece mapping in which m = K for K greater than Ki and m = m for n < K . A n exponential b b function is used to represent susceptibihty values between K and K \ . Figure 4.1 shows B the mapping. The forward and inverse mappings, = / ( « ) and m K = f~ (m), 1 are given by K < m(«) = Kl [in +1 K B < b K K K < Ki (4.1) K > Ki and m i < m /c(m) K\e K b K i m b m m < m < <m b mi (4.2) Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 97 Figure 4.1: Nonlinear mapping for positivity constraint. Segment 1 is an exponential function and segment 2 is a straight line whose slope equals 1. Segment 3 is a straight line parallel to the m-axis. where mj is m = h K l [ln (^j + lj . (4.3) With this mapping the recovered susceptibility has a minimum value of «{, but this is chosen small enough so that its effect on the data is insignificant compared to the errors on the observations. 4.3 Calculation of the sensitivities The calculation of sensitivities J ; J = ddi/drrij is an important part of the inversion algorithm. Here I use the adjoint Green's function solution to calculate the sensitivities. There are a number of choices for the data and for the definition of m , but all of the sensitivities can be obtained once dE/dfij is known. For instance, if the secondary fields 98 Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data are measured, the sensitivities for H are z orrii itoui Jo orrii Because \i = Mo(l + K) the sensitivity for any m = / ( « ) can be easily generated. 4.3.1 Adjoint Green's function solution The geometry used in this chapter is the same as that in Chapter 2. The sensitivities can be computed using a modified adjoint Green's function solution. The Maxwell's equations for 1-D problems are (Ryu, 1970) itouH (r,e,uj,z) = r iupH (r, 8,u;,z) = - - {£[rE (r, 6,u,, z)]} , 1 z ajiA^A e _ SEd^A = { i u j e Q + a ) E e ( r j d (4.5) , , z) + I , w s where r, 0 and z are variables in the cylindrical coordinate system. I is the current s source given by _ I(u>)a6(r - a)S(z obs - h) r Due to the symmetry of the problem, electric and magnetic fields are no longer functions of 9. For simplicity, I use E to denote EQ. The permeability fi is a function of depth and for our layered earth it is represented as M M*) = £/Mfc(*), ( -6) 4 i=l where M is the number of layers, and tpi(z) is the box car function which is unity on the support of the i t h layer and zero elsewhere. Substituting equation (4.6) into equation( 4.5) and taking derivatives with respect to \i{ in each layer yields the partial differential equation for the sensitivity problem. In the Hankel transform domain, it is Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 99 given by dm dz tlx -[S(Z - Zi) - 6{z - z i + 1 ) \ + IU (4.7) dz dfj-i where the operator C is d*_ c (4.8) u dz 2 Detailed derivations of the above equations is given in Appendices A and B , where a Green's function is introduced into the problem. Multiplying both sides of equation ( 4 . 7 ) by the Green's function and integrating by parts we obtain d dE 2 G r<*> dE dG d dE dz dfn ' J-oo J-oo dfii d^ dzdzdmW^ 2 iiotpi(z)(zuie 0 uG 2 dz G dE — Pi OZ oo r / G - + ai)GEdz J— oo (4.9) [dz (4.10) The boundary term on the left hand side vanishes because the electric field for any finite source tends to zero at infinity and so do its derivatives. Thus if the Green's function satisfies the equation grG(A, u, z) - u G(\, to, z) = S(z - 2 G(A,w,z)|, , = G(\,w,z) r —> 0 z ) ohs =G{\,u,,z)\ (4-11) t=zt when \z\ —> oo, then the sensitivity for the electric field is dE(\,co, dfii z ) ohs / Jzi . ,. iu>(iuje + o-i)G(X,u, 0 z)E{\,w, z)dz - Z i + l Pi dz (4.12) where E is the primary field in the ith layer and G the corresponding Green's function. The primary field is that produced by the transmitter and the auxiliary field G is due to 100 Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data a vertical magnetic dipole with unit strength at the observing point z . They are given obs by Ei(X,u,z) = Ai(X,co)e ^-^ + B (A,a;)e- '^- ') Gi{\,u,z) = ai(\,u)e^ -^ + u u z i z b (X,o )e- ^-^ u i J where the coefficients of the up-going wave and down-going wave are given by the following formulae: A (\,u) = B (\,u,)e-^§^ i i = < ' ) 4 bi(\,w) = . fi \ Bi 'V ' / - B (\,io) 0 , w v^mlaix 1 4 (Aa) j^T ^ - The input impedance and intrinsic impedances can be calculated by using equations (2.7) an (2.8) in Chapter 2. For the convenience of computation, equation(4.12) can be reorganized to (Appendix C) dE(\,Lo,z ) ^ 2u\ r = oba dm +iu(iu>e +<Ti)Jl m i+1 0 rzi+i / • G(X, u>, z)E(\, to, z)dz - 2Aibihi J G(X,u,z)E(\,Lo,z)dz. (4.15) Appendix D outlines the detailed calculation of equation(4.15). Equation (4.15) shows the sensitivities consist of two terms. The first term, which represents the influence of sources on the two boundaries of each layer, will not tend to zero at zero frequency. The second term, on the other hand, will tend to zero when the frequency tends to zero. This is shown in Appendix E . To better understand the asymptotic behavior of the sensitivities, consider the value of the magnetic field when the frequency tends to zero. Let Z = Km{~)Z w-»0 no i i = m, (4-16) 101 Chapter 4. Recovering susceptibility from 1-D Inversion of EM data be the normalized intrinsic impedance and i ^ h m ( - l ) ^ <"-»° iu> = / ./ ' i + 1 + / X pi + pi y \ tanh(A/j.j) t a n +1 / l (4.17) ) be the normalized input impedance. The vertical magnetic field at zero frequency can then be expressed as H^wzob.) = al f°° \ ! X Z •/-oo ~ ° e (^- )j (Aa)J (Aa)cJA. /i ) A (4.18) 2fc 1 2(Z + fi ) 1 0 0 If pi = po, where i = 1,2, . . . , M , then Z — /x and hence H (ui) will tend to zero. x 0 z Conversely if any one of the layers has nonzero susceptibihty then the magnetic field will be nonzero. For a half-space with a conductivity a and a magnetic permeability ft illuminated by a dipole of moment m and with both source and receiver sitting at the same height h, the solution is reduced to = 2 8fe -r mp^po 47T/Xx 2 +p (4h 2 0 - r - 2 ) 2 2 - 5 ' V ' ; The derivation of the above result is given in Appendix F . Thus the induced static magnetic field due to a dipole located at h above the surface is equal to that of an image dipole of strength m(pi — po)/{pi ~ po) buried at a depth h beneath the surface. The effective magnetic charges are on the surface of the half space, where the susceptibihty is discontinuous. Similarly, the two boundary sources in the calculation of sensitivities can be viewed as layers of effective magnetic charges. Those surface magnetic charges are due to the sudden change of susceptibihty on boundaries and they remain as the frequency goes to zero. As frequency increases, eddy currents will be come stronger, and this adds a frequency dependent term to the sensitivities. A numerical example of the sensitivities for horizontal coplanar coils with a coil separation of 10 meters and at height 30 meters above a half-space of 10~ S/m and a 2 magnetic susceptibihty of 0.1 SI unit is given in Figure 4.2. Panels (a) and (b) show how Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 102 the amplitudes of real and imaginary components of the sensitivities vary with respect to frequency. At low frequency the real component remains at a constant value due to magnetic polarization. As frequency increases, the induced currents become stronger and the real component of the sensitivities become frequency-dependent. At a certain frequency, in this case around 10 hertz, the frequency-dependent term begins to dominate. For a 3 given susceptibihty structure the transition frequency decreases as conductivity increases. The imaginary component, on the other hand, is completely frequency-dependent. A t low frequencies, magnetic particles change orientation almost synchronously with the primary field, and therefore the amplitude of the imaginary component of the sensitivities is very small. It rises almost linearly with respect to frequency until 1000 hertz. As frequency arises further, the frequency-dependence is no longer linear, and a local maximum with respect to frequency is formed at about 10000 hertz. Both components of the sensitivities also behave differently with depth. The real component begins at a constant value and decreases with depth, but for the imaginary part, there is a depth at which the sensitivity is maximized. This characteristic contributes greatly to the depth resolution obtained in the inversion. Figure 4.2c presents the absolute value of the frequency independent part of the sensitivities as a function of depth and coil separation. When the coil separation is much smaller than the observation height, as in the case of this example, this term generally decreases with depth. The frequency-dependent term of the real component of the sensitivities is obtained by subtracting the frequency-independent term from Figure 4.2a. Figure 4.2d shows the result for a source-receiver separation of 10m. It is similar to the imaginary component of the sensitivities in that at frequencies higher than 1000 hertz it also reaches a maximum at depth. At frequencies lower than 1000 hertz, the amplitude of the frequency-independent component of the sensitivities is linear to frequency. However, the amplitude of the Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 103 Figure 4.2: The absolute value of sensitivities (on logarithmic scale) over a half space. The conductivity is 0.01 S/m and the susceptibility is 0.1 SI unit. The data were calculated for a coil separation of 10 meters and survey height of 30 meters, (a) Real component of the sensitivities; (b) Imaginary component of the sensitivities; (c) Frequency independent term in the real component of the sensitivities; (d) Frequency dependent term in the real component of the sensitivities. 104 Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data frequency-dependent sensitivity is much smaller than the frequency-independent component; thus, during an inversion, the frequency-independent component of the inphase part of the sensitivities may limit depth resolution. 4.4 Synthetic examples For all the synthetic data in this section I assume a coplanar system with coil separation 10m in which the transmitter has unit area and carries a harmonic electric current of 1 Ampere. The earth is divided into 44 layers and the thicknesses of the layers increases with depth to compensate for the loss of resolution with depth. The conductivity structure is assumed known and the mapping parameters K and Ki are fixed at 10~ and 10~ 6 3 b SI unit. As a first example, I invert data from a ground system that is 0.5m above the surface. The conductivity model, which is shown in Figure 4.3a, consists of a 40m conductive overburden of 0.1 S/m and a conductive zone at depth between 40m and 80m. The susceptibihty model consists of a single layer of 0.2 SI. This susceptible layer is 30m thick and straddles the boundary of the conductive overburden and the conductive zone at depth. The overburden and the offset of the conductive and susceptible zones make this a complicated but good example to test the reliability of the algorithm. The data are calculated at 10 frequencies: 110, 220, 440, 880, 1760, 3520, 7040, 14080, 28160, 56320 hertz and contaminated with 2% Gaussian noise. The model parameter used in this inversion is related to K through the nonlinear mapping given by equation (4.1). Parameter a in equation (2.25) is set to 0.02. The starting and reference susceptibility models were half-spaces of values 0.0 SI and 1 0 - 6 SI respectively. The parameter 7 was chosen as 5. The results are shown in Figure 4.3. After 7 iterations the inversion converged to the desired target misfit of 20. Figure 4.3b shows the reconstructed and true susceptibihty Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 105 Figure 4.3: Inversion of ground system data, (a) The conductivity structure used in the inversion. Its value below 80 meters is 0.01 S/m. (b) Recovered (solid line) and true susceptibihty models (dashed line), (c) Misfit curve for the inversion, (d) Model norm as a function of iteration. models. Figures 4.3c and 4.3d show the misfit curve and model norm as functions of the number of iterations. The inversion has recovered a very good representation of the true susceptibihty structure. In a second example, whose results are given in Figure 4.4,1 invert data from a typical airborne survey in which inphase and quadrature phase data at frequencies 900, 7200 and 56000 hertz were collected at a flight height of 30 meters. The data were calculated over the same model used in the previous example. It is more difficult to fit data with low noise level than with higher noise level. Hence I test the robustness of the algorithm by contaminating the data with 0.5% Gaussian noise. Starting and reference models in this example were 0.02 and 1 0 SI respectively. The true conductivity model, given in Figure - 6 4.4c, is assumed known. Two approaches to incorporate positivity will be illustrated. In Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 106 Figure 4.4: Inversion with correct knowledge of conductivity structure, (a) Recovered (solid line ) and true susceptibihty models (dashed line) from the inversion in which ln« was used as the model parameter (dashed line), (b) The true (dashed line) and recovered susceptibihty from the inversion in which the nonlinear mapping was used to incorporate positivity. (c) the true conductivity used! in the inversion, (d) The misfit curves corresponding to (a) (dashed line) and (b) (solid line). Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 107 the first I used logarithm of susceptibihty as the model parameter. The reconstructed model, obtained after seven iterations, fits the data to the desired level, but overshoots the true model. This is primarily the result of using ln(/c) as a model parameter. Using the nonlinear mapping [equation(4.1)] to guarantee positivity yields the model in Figure 4.4b. This model is a better representation of the true model and was obtained in four iterations. In both cases, the inversion converged to the desired target misfit 6. Plots of the data misfits for both inversions are provided in Figure 4.4d. In carrying out the inversion, parameters a and 7 were chosen as same as those in the previous example. The primary contribution to the E M responses is from eddy currents induced in the earth and the magnitude of the data is dependent upon the electrical conductivity structure. It follows that inversions for susceptibihty which are performed with incorrect knowledge of the electrical conductivity will suffer some deterioration. I illustrate this by repeating the last inversion but this time using an approximate conductivity model. The conductivity model is obtained by using a separate inversion algorithm which recovers a 1-D electrical structure from horizontal loop E M data by assuming that \i and e take their values in the free space. The algorithm was terminated after the second iteration when the misfit was 75.4, well above the desired value of 6. The true and the approximate conductivity models are shown in Figure 4.5a. Now I use the approximate conductivity and invert for K. The algorithm plateaued to a minimum misfit of about fa = 28 after 9 iterations. The recovered susceptibility in Figure 4.5b shows increased susceptibihty at about the right depth, but it overshoots the true model significantly and is not an accurate representation of the true model. This discrepancy increases as the conductivity model becomes a poorer representation of the true conductivity. Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 108 Figure 4.5: Effect of incorrect knowledge of the conductivity distribution on the inversion, (a) Solid line denotes the approximate conductivity model and dashed line denotes the true model, (b) The resultant susceptibihty model (solid line) and the true model (dashed line), (c) The misfit curve for the inversion, (d) The model norm as a function of iteration. Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 4.5 109 Field example As a field example I now invert airborne E M data acquired at M t . Milligan which is a Cu-Au porphyry deposit located in central British Columbia, Canada. The inphase and quadrature phase data at frequencies 900, 7200 and 56000 hertz were taken about every 10 meters along the flight line. The coil separation is 8.0 meters for 900 and 7200 hertz data and 6.3 meters for data at 56,000 hertz. The D I G H E M system was flown in north-south lines 100 meters apart. Even though the flight lines are north-south I invert data at 13 stations along an east-west line (Y9600). This is because D C resistivity data were collected on east-west lines and they have been inverted by Oldenburg, L i and Ellis (1996). I will use their 2-D conductivity model in our susceptibihty inversion. The D I G H E M data for the 13 stations are given in Figure 4.6. The real component of the vertical component of the magnetic field at 900 hertz is negative at most of the stations. There are also some negative inphase data at 7200 hertz. For airborne electromagnetic surveys, the negative inphase data is a direct result of magnetization. For such surveys, the flight height h is generally much greater than the coil separation r, and therefore the secondary magnetic field recorded at the receiver opposes the primary field. A t low induction numbers the magnetic field due to magnetic charges at boundaries of susceptibihty discontinuity can exceed the secondary fields generated by eddy currents in the earth, and since they oppose each other, it is possible for the inphase portion of the airborne E M response to be negative. At higher induction numbers, however, the effect of eddy currents will dominate. In performing the inversion the earth is divided into 22 layers which is the same as the number of layers used in the inversion of the 2-D D C data. A contoured representation of that conductivity model is shown in Figure 4.7a. The laterally averaged conductivity Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 12.4 12.8 i 13.2 110 i X(km) 4 12.4 12.8 13.2 X(km) Figure 4.6: D I G H E M data from M t . Milligan, at section Y9600. The real component is denoted by the solid line and imaginary component is denoted by the dashed line. The flight height varies between 25.2 and 48.8 meters. The coil separation is 7.98m at 900 and 7200 hertz, and 6.33m at 56,000 hertz. Chapter 4. 111 Recovering susceptibihty from 1-D Inversion of EM data beneath each station was used as a 1-D background conductivity for the susceptibihty inversion. The model objective function was that given in equation(lO) with the parameter a set to 0.02. The model parameter TTI(K) for the inversion is connected to susceptibihty through the nonlinear mapping given in equation(4.1). The mapping parameters K and B K i were set to be 10~ SI and 10~ SI respectively. The starting model for inversions at 6 3 all the stations was a half space of 0.02 SI. The reference model was a half space of 10~ 6 SI. The noise level in the data is assumed to be 10% of the amplitude of the data. This resulted in a minimum standard deviation of less than 1 ppm for some data and is likely to have been overly optimistic. The cumulative initial chi-squared misfit for the 13 stations was 417,206. The result of the inversion is shown in Figure 4.7b and the cumulative misfit has been reduced to 1795. Three regions of high susceptibihty are observed in the upper 200 meters and the maximum susceptibihty is 0.1 SI. This result can be compared with Figure 4.7c which shows magnetite concentration provided by DeLong et al. (1991) which was visually estimated from borehole samples over the same section. Topography information is incorporated in that Figure. The larger magnetic anomaly in the center is supported by 4 observations. The highest value in this anomaly is 8% and the magnetite content for the 3 other supporting points are 5%. The high susceptibihty at station 12.9 km in Figure 4.7b corresponds well, both vertically and horizontally, with the borehole information. There is an indication on Figure 4.7c of enhanced magnetite content near 13.1 km and another more elongated concentration near 12.9 km. These are not pronounced features, but they do correlate with the inversion result in Figure 4.7b. The data from a ground magnetic survey at M t . Milligan have been inverted to recover a 3-D model of susceptibihty (Li and Oldenburg, 1996). The cross-section from the recovered 3-D susceptibihty model is presented in Figure 4.7d. Three concentrations of susceptibihty are observed, with the largest amplitude of 0.047 SI occurring at 12.7 Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 112 km and at a depth of 200 meters. This is considerably deeper than the susceptibihty recovered by inverting the airborne E M data. Figures 4.7b and 4.7d both indicate high susceptibility at 12.7 km but there is a lateral difference of about 100m between the locations of the right most anomaly. With the exception of this lateral shift in the right hand anomaly, the greatest discrepancies between the models exist in the vertical direction. One possible explanation is that the depth of investigation, which is primarily controlled by skin depth and geometry, is less than 150 meters in this case. Therefore, the airborne E M data are primarily sensitive to structure in the top 150 meters and therefore structure with susceptibility lower than 0.1 SI at 200 meters will not greatly affect the data. This has been confirmed by forward modelling. Another possibility for the disagreement between Figs 4.7b and 4.7d is due to nonuniqueness in the inversion. The recovered model from the inversion of static magnetic data seems deeper and more spread out. Since the depth distribution in the 3-D model is a consequence of the depth weighting in the objective function, inappropriate design or use of that weighting function may affect the inversion. On the other hand, the quality of the results of the 1-D susceptibihty inversion can be affected by 3-D variations in conductivity and susceptibility, which surely exist at M t . Milligan, and by incorrect estimation of the background conductivity. Ideally I would like to incorporate the 3-D effects into the errors ascribed to the data, but I do not know how large these are. In the inversion in Figure 8b we assigned a constant percentage error. Other reasonable errors assignments are: (1) constant base level plus a percentage of the data; (2) a fixed but different value for each frequency and (3) uniform errors on all data. For a given data set, the inverted model depends upon the assigned errors and how well the data are misfit. To investigate this variability we carried out the inversions with different error assignments and additionally imposed a reasonable upper limit of 0.1 SI units on the recovered susceptibility. In Figure 4.8a I show the inversion result when the standard Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 113 8.258 7.395 6.533 5.670 4.808 100 H 3.945 Q- 200 3.083 1000 (S/m) CD Q 2.220 12.4 12.6 12.8 13.0 13.2 13.4 1.358 0.495 101.10 90.95 80.83 70.71 60.59 50.48 40.36 30.24 20.12 12.4 12.6 12.8 13.0 13.2 13.4 10.00 Figure 4.7: Inversion of D I G H E M data from M t . Milligan, section Y9600. (a) Recovered conductivity model from the inversion of 2-D D C data, (b) Susceptibihty model reconstructed from the 1-D inversion of D I G H E M data, (c) Magnetite content i n percentage from borehole information, (d) Susceptibihty model from the 3-D inversion of static magnetic data. Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 114 deviation for data at 900, 7200 and 56000 hertz is 5 ppm plus 10 percent of the strength of the data. For Figure 4.8b the standard deviations were 1, 4 and 10 ppm for data at the three respective frequencies and in Figure 4.8c the standard deviation for each datum was 10 ppm. There are differences between the three sections but they all identify anomaly highs at 12.6 km and 12.9 km. A l l susceptibility highs are concentrated within the top 150 meters. This provides confidence that the algorithm is producing meaningful results. 4.6 Conclusions The work presented here shows how electromagnetic data from a horizontal coplanar loop can be inverted to recover a 1-D susceptibihty structure under the assumption that the electrical conductivity is known. Since the strength of induced magnetization inside the earth depends upon the amplitude of the existing magnetic field, it follows that E M data at different frequencies are sensitive to susceptibilities at different depths. This is in contrast to static magnetic field data. M y algorithm follows traditional inversion methodologies for solving under determined nonlinear inverse problems and minimizes an objective function subject to fitting the data. Positivity is incorporated by using ln(rv) or using a nonlinear mapping of susceptibihty as model parameters. Synthetic inversions indicate that convergence with the nonlinear mapping usually requires fewer iterations to achieve the same misfit and generally produces a better representation of the true model. For field data, when using a 1-D inversion algorithm in complex environments, one is faced with the ubiquitous problem of specifying the observational errors and deciding how well the data should be fit. This remains problematic but in our example I used a variety of error assignments and imposed an upper limit on the constructed susceptibilities. The resultant images had common features and the main feature coincided with a Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 400 500 1 0 12.4 12.6 i i r 12.8 13.0 X (km) 1 I 13.2 13.4 i /alto: 100 £ 200 J 300 K b 400 IN x 1000 SI 500 12.4 J 200 f 300 12.6 12.8 13.0 X (km) 13.2 13.4 H K 400 x 1000 Si 500 12.4 12.6 12.8 13.0 X (km) 13.2 13.4 115 82.020 72.907 63.794 54.680 45.567 36.454 27.341 18.227 9.114 0.001 90.609 80.542 70.474 60.406 50.339 40.271 30.204 —j- 2 0 . 1 3 6 — 10.068 I—L 0.001 74.563 66.278 57.994 49.709 41.424 t- 3 3 . 1 4 0 24.855 16.570 8.286 0.001 F i g u r e 4.8: R e c o v e r e d s u s c e p t i b i l i t y m o d e l s f r o m i n v e r s i o n s w i t h different e r r o r schemes: (a) t h e r e c o v e r e d m o d e l w h e n t h e s t a n d a r d d e v i a t i o n s were 5 p p m p l u s 10 p e r c e n t o f t h e s t r e n g t h o f t h e d a t a ; ( b ) t h e s t a n d a r d d e v i a t i o n s were 1,4 a n d 10 p p m for d a t a at 9 0 0 , 7200 a n d 56,000 h e r t z ; a n d ( c ) a c o n s t a n t s t a n d a r d d e v i a t i o n o f 10 p p m was u s e d f o r a l l the data. Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 116 region of high magnetite content inferred from visual estimates of borehole logs. Reasonably accurate information about background conductivity is important for the inversion. If the true conductivity is known, the inversion can produce a good representation of the true magnetic susceptibihty. However, when the conductivity is not accurate, the recovered susceptibihty model will be distorted. This invites the challenge of carrying out simultaneous inversion of conductivity and susceptibihty. The method outlined in this Chapter is qualitatively useful. When accurate information about conductivity structure is available from other geophysical surveys such as a DC resistivity survey, our method may provide useful information about susceptibihty structure. It can also be used for depth mapping of fresh water lakes in shield terranes, or for airborne mapping where a non-magnetic overburden lies over resistive and magnetically permeable bedrock, so the conductivity structure is relatively separate from the susceptibihty structure. It can also be potentially useful for mapping magnetite based on magnetic susceptibihty, when these occur in resistive rocks. Chapter 5 1-D Simultaneous Inversion of E M data 5.1 Introduction Both conductivity and susceptibihty are important physical parameters. The traditional way to obtain information about susceptibihty is through inversions of static magnetic data. Since E M surveys are not affected by remanent magnetism, and the artificial sources used in E M surveys are highly localized compared to the relatively uniform geomagnetic field in magnetic surveys, E M surveys can provide complementary information about susceptibihty. Although the problem of 1-D inversion of electromagnetic data in both the timeand frequency-domains has been studied extensively in the literature, most of these studies have assumed knowledge of either the conductivity or the susceptibihty. A typical procedure in conductivity inversion is to assume that magnetic susceptibihty equals its value in free space. In many cases this assumption is valid since most rocks are nonmagnetic. However, quite often the geological targets are not only conductive but also magnetically permeable, so the data are affected by both conductivity and susceptibihty. A common example of the existence of strong magnetization is the negative inphase coplanar data in airborne E M surveys as discussed in Chapter 3. If I invert those data under the assumption that fi = /z , then the recovered conductivity model could be very 0 wrong, and valuable information about susceptibihty in the data is also wasted. On the other hand, when inverting E M data to recover susceptibility, incorrect knowledge 117 Chapter 5. 1-D Simultaneous Inversion of EM data 118 about conductivity may also cause severe distortion in recovered susceptibility models (Zhang and Oldenburg, 1996a). To illustrate this I invert a synthetic data set generated over a 1-D earth with variable conductivity and susceptibihty. The data were calculated at 900, 7200 and 56000 hertz. The coil separation is 10 meters and the observation height is 30 meters. Due to the presence of the susceptibihty structure, the inphase datum at 900 hertz has a negative value of -13.0 ppm. Gaussian noise with standard deviation of about 1% of the data strength was added to the data. In the first inversion I attempted to recover the susceptibihty distribution while specifying that the earth's conductivity is a half-space of 0.001 S/m. The recovered susceptibihty in Figure 5.1a is not a good representation of the true model. In the next inversion I used correct information about the conductivity. The recovered susceptibihty model represents the true model very well (Figure 5.1b). I carried out the same tests on conductivity inversion. Figure 5.1c shows the recovered conductivity model under the assumption that fi — /io- The recovered model is distorted severely at depth and overshoots the true model. However, when true information about susceptibihty is used, the recovered conductivity model, shown in Figure 5.Id, delineates the true model quite successfully. Those results show that in an individual inversion we need to have accurate information about either susceptibihty or conductivity in order to recover the other. In principal, I may invert D C resistivity data or static magnetic field data to provide information about conductivity or susceptibihty structure. I may then use the conductivity information from the inversion of D C resistivity data to carry out inversion of frequency- or time-domain E M data to recover susceptibility structure, or use the information about susceptibility from the inversion of static magnetic data to invert the E M data for conductivity structure. The problem is that I do not always have D C resistivity or static magnetic data along with frequency- or time-domain E M data. Therefore the ideal way to attack this problem is to recover both conductivity and susceptibihty at the same time through a simultaneous Chapter 5. 1-D Simultaneous Inversion of EM data 119 Figure 5.1: Individual inversions with accurate and inaccurate information about conductivity or susceptibihty. Solid lines denote the recovered models, and dashed lines denote the true models, (a) The recovered susceptibihty from the inversion with inaccurate information about conductivity. The conductivity model used in the inversion was a 0.001 S/m half-space, (b) Recovered susceptibihty from the inversion with accurate information about conductivity, (c) Recovered conductivity from the inversion under the assumption that the susceptibility is equal to its free-space value, (d) The recovered conductivity when accurate information about susceptibility was used in the inversion. inversion. In this Chapter I present a method to solve the simultaneous inverse problem in a layered earth for a horizontal coplanar E M system. The number of layers is chosen, based upon the estimation of the apparent conductivity, to be large enough to adequately represent the possible conductivity and susceptibihty structures. The thickness of each layer is fixed and increases with depth to compensate for the associated loss of resolution of the data due to the attenuation of the E M fields. In each of these homogeneous layers a pair of model parameters for the conductivity and susceptibility needs to be recovered. Methods similar to the one presented in this Chapter were previously used in the inversion 120 Chapter 5. 1-D Simultaneous Inversion of EM data of either the susceptibihty (Zhang and Oldenburg, 1995, 1996a) or conductivity structure (Fullagar and Oldenburg, 1984) in a 1-D environment. Here however, I simultaneously invert for two model parameters, <r and K. I show that the data can be decomposed into two parts to reflect the effect from susceptibihty. eddy currents through the term jumper, One part, which is related to is not resolvable without prior information. The other part is due to magnetic polarization and carries independent information about susceptibihty. I use a weighted sum of model objective functions of conductivity and susceptibihty to construct the cost function. I use synthetic and field data examples to show that simultaneous inversion can provide useful information. 5.2 Methodology Some aspects of techniques used in a simultaneous inversion have been addressed in previous Chapters. The forward modeling has been solved in Chapter 2. The sensitivities for conductivity and susceptibihty are given in equations (2.21) in Chapter 2 and (4.15) in Chapter 4. The question now is how to carry out the simultaneous inversion. The goal of the inversion is to find a model which reproduces the data and exhibits desired characteristics. M y choice for the objective function is guided by the desire to find a model that has minimum structure in the vertical direction and at the same time is close to a reference model. To accomplish this I set up the model objective functions for conductivity and susceptibihty as fa = oia J w\{z) hi {^j dz + (l- a) J a d(lna — w (z) ln<r ) 0 (5.1) dz, 8~z 2 and <f>* = OLK J w (z)[m(K,) 3 - m(n )} dz 2 0 + (1 - ct ) j K d[m(n) — m(/€ )] w (z) 4 0 dz dz, (5.2) Chapter 5. 121 1-D Simultaneous Inversion of EM data where cr and Ko are the reference models for conductivity and susceptibihty. The pa0 rameters a and a a K control the relative importance of smallest and flattest components in the model objective functions. The use of ln<r as the model parameter ensures the recovered conductivity is positive, and also accommodates the wide range of conductivity variations. The nonlinear mapping in equations (4.1) and (4.2) is used to project the susceptibihty into TO(K). This mapping can provide positivity constraints on the recovered susceptibihty model and prevent small values of susceptibihty from carrying too much weight in the inversion. For the discrete 1-D inversion, those two model objective functions can be rewritten as W^lnf—) (5.3) 4><T = and 4>* = \\W [m(K) K - (5.4) rn(K ) 0 where W and W are M x M weighting matrixes. a K One of the most distinguishing aspects of a simultaneous inversion problem is that there is more than one objective function that requires minimizing. Just as I combined two terms to make up the objective function in equations (5.1) and (5.2), here I can combine the two objective functions <f> and (j> into my cost function. Let the final cost a K function be 4>m = Q<I><T + 1<}>K, (5.5) where coefficients g and 7 are given by 1 1+ s (5.6) 1+ s where 0 < s < 00 is the desired magnifying factor. When s —> 0, (j) = <t>*, m a n d when s —> oo,(j> = 4> . Now I solve the simultaneous inverse problem by minimizing the cost m K Chapter 5. 122 1-D Simultaneous Inversion of EM data function subject to the constraint that the data are adequately reproduced: minimize <f> = (gfa + >y<j> ) + 3 ((f> l K where 3~ is a Lagrange multiplier and (f) x d 4> ), tar (5.7) is the target misfit level. In above equation, tar the data objective function <f> is the same as used in Chapter 2 [equation (2.27)]. Let d m a = ln(<r) and m = m(n) be vectors with M components, and let Sm^ and 6m denote K K perturbation at the nth iteration. The predicted data can be approximated as 5[m[ where J and J a K ddi/dm i. JJm , (5.8) are the sensitivities whose elements are J u = ddi/dm i and J u = n) + Sm^mW + Sm ] « F[m<?\m^} + JJm, K + K a a R Let J = ( J , J ) be a global sensitivity matrix, m = (m ,m ) be a global K a K CT re model parameter, and W m = 0 ( gW a V 0 W 1 K \ (5.9) ) be the global weighting matrix. The linearized problem becomes the minimization of d> = 8 W [8m + m -m ] 0 ' + j W {D ohs d - F[m^] + JSm}\ 2 - <j> } , tar (5.10) where cj) is the target level for data misfit at the (n + l ) t h iteration. Usually I reduce tar the target misfit level from one iteration to the next by a factor between 2 to 10. The above equation can be solved as in Chapter 2, and the solution has exactly the same expression as in equation (2.31). 5.2.1 The trade-off between conductivity and susceptibility Magnetic susceptibihty is related to magnetic permeability through p = fi (l + K), and 0 for convenience I use both K and fi in the following discussion. The governing partial 123 Chapter 5. 1-D Simultaneous Inversion of EM data differential equation for the electric field, in the Hankel transform domain, is given by equation (2.5): rd 1 2 — -u E(X,z,co) = iLOfi 3 (Xa)6(z 2 0 1 - h), obs (5-H) and the boundary conditions are the continuity of the tangential components of both electric and magnetic fields. The general solution for the above equation is the linear combination of exponential functions e . ±UjZ The independent information about the susceptibihty alone enters from the boundary terms related to the calculation of input impedances i- Z i _ E H Z ^ + Z 3 rj Z tanh(u h ) j j j Zj + ZJ+ t a n h ( ^ ) ' ^' 1 ; where Zj = —^ - Physically the secondary field results from both induced eddy curL rents and magnetic polarization. Eddy currents are related to the term jtoficr, and as long as the product of fi and a remains the same, the forward response is not affected. The simultaneous inversion is useful only when the secondary fields caused by magnetic polarization become sufficiently large. I illustrate this through a simple example. Consider a half-space of conductivity cr and permeability fi. In the Hankel transform domain, the vertical component of the secondary magnetic field for such a half-space, under the quasi-static assumption, is H {X,co,z ) z =B o e ^ obs 2 h - ^ ^ ^ , (5.13) Xfi + ufi 0 where u = A + jtoficr, and X is the Hankel transform variable. B is given by 2 2 0 Bo = (5.14) where J i is the Bessel function of the first kind. For a conductive half-space and a magnetic whole space with permeability u., the forward response is H (X,u,,z ) z obs = e-^ -^\^. (5.15) 2h Bo X+u 124 Chapter 5. 1-D Simultaneous Inversion of EM data The information about the susceptibihty comes in only through the term jwficr, which is contained within u. Since fi = poPr = /*o(l + «), the term juipa can be written as j(jJlioH a — ju/poa, r where a — \i a. r In another words, I can use a non-magnetic half- space with conductivity a to generate the same data as those that I would obtain from a magnetic whole space of \i and a conductive half-space of a. The difference between equations (5.13) and (5.15) tells me how strong the influence is from the discontinuity in the susceptibihty at the boundary: AH (\,u,z ) z = Boe-**-^ obs (^_^ _ ( 5 . 1 6 ) The data are functions of conductivity, susceptibihty, frequency, and geometric factors, but for simplicity I write the data as D(cr, //,). So in general I can define AD(a, /*) = D(cr, (i)- D(v, /*>). (5.17) This quantity can serve as a measure of the influence of discontinuity in susceptibility distribution on the data. Since the independent information about the susceptibihty structure is contained in AD((r,fi), I may obtain a useful result from the simultaneous inversion only if the amplitude of AD(cr, fi) is greater than the noise level in the data. For a half-space, AD(a, /i) can be obtained by carrying out an inverse Hankel transform to equation (5.16). Even for such a simple example, I still cannot find an analytic solution for above equation. Only when the frequency goes to zero, and when the source can be considered as a vertical magnetic dipole, can I obtain a closed expression for AD(cr, fi): UmA (.,^) = - ^ f ^ ) g J * ' . (5.18) The derivation of above equation is very similar to the work of Zhang and Oldenburg (1996a), who developed the asymptotic expression for the forward modeling. Note that Chapter 5. 1-D Simultaneous Inversion of EM data 125 Figure 5.2: The relationship among D(<7,/x), D(a,fi ), that \i a = ex. D(a,u, ), 0 0 and AD(a,fi). Note T when frequency tends to zero, there is no induction, so that conductivity disappears from the asymptotic expression of AD(a,fi). To obtain some insight into the magnitude of the trade-off between conductivity and susceptibihty, I plot the relationship among D(a,fi), in Figure 5.2. Table 5.1 shows D(cr,/i), Z?(<T, i x ) , 0 D(cr,/x ), D(a,fi ), 0 D(a,/A ), 0 0 and /L\D(a,fi) and A Z ) ( c , / x ) calculated for a half-space of 0.01 S/m and 0.1 SI. From basic E M theory, the angle <j> in Figure 5.2 varies from 0 for a very resistive earth to 7r/2 for a very conductive earth. Since fi r > 1 for a magnetically permeable earth, the effective conductivity a reinforces the eddy currents. In this example, at different frequencies a increases the real component of D(cr, /j, ) by about 3 % to 10%, and 0 the imaginary component by about 1.5% to 8%. The real component of the magnetic field associated with the effective magnetic charges is anti-phase to r\eD(<r, /^o), and the amplitude of jReAD(cr, fi)\/\ReD(<r,fi )\ 0 56k is about 4 0 0 % at 900 hertz, and about 1 4 % at hertz. These effective magnetic charges increase the amplitude of ImD(cr, fi ) from a 0 126 Chapter 5. 1-D Simultaneous Inversion of EM data Real(Hz) (ppm) 900 7200 56000 2362 -347 171 52.9 528.1 2595 59.5 578.3 2737 -375 -406 -407 Hz Freq (hertz) D(*,fi) D((T, fXo) D(a,ft, ) 0 AD(*,p) Table 5.1: D(o~, po), D(a,^ ), half-space. AD(a,fi), 0 Imag(Hz) (ppm) 7200 56000 900 220.1 970.9 2115 202.9 908.5 2021 219.5 959.8 2052 11.1 63.1 0.6 and D(o~,fi) over a 0.01 S/m and 0.1 SI unit negligible 0.6 ppm at 900 hertz to a maximum at 56k hertz of about 1.5%. This means that it would be difficult to recover susceptibihty distribution by inverting the imaginary component of the data alone. But on the other hand, this analysis suggests that when carrying out individual inversions without the knowledge of susceptibihty, the imaginary component of the data should be used to obtain information about conductivity, even though conductivity models obtained this way may be a little more conductive than they should be. Since the inverse problem is nonlinear, linearization and iteration are needed. At each iteration a perturbation is sought and the model is modified by this perturbation. The perturbation on the model will cause a change on the data 6D(a, n) m J^ba + J SK, (5.19) K where are and J are the sensitivities for conductivity and susceptibihty, whose elements EtiWu K a n d ZiLi(J*)u, 1 = 1,2, ,N, and i = 1 , 2 , . . . , M . N and M are the numbers of the data and the model cells respectively. In general, the perturbation of the fcth datum over a 1-D earth can be written as M 8D ((T, K) W K [{J<r)liS<Ti + (J )li6Ki\ K Chapter 5. 1-D Simultaneous Inversion of EM data 127 M £ (Ja)li f + —SKi ] + QuSKi (5.20) The proof of the above equation is given in the Appendix (F). Mathematically equation (5.20) indicates that the perturbation on the data is composed of three parts. The first two terms on the right hand side of the above equation are related to the term jcofia, and they are the same except a scaling factor. The independent information about susceptibihty is contained in the third term Q^i, which is the partial derivative of the data with respect to the susceptibihty in the boundary conditions. Further more, this boundary term retains a non-zero value as the frequency tends to zero, as proven by Zhang and Oldenburg (1996a). Physically equation (5.20) means that the perturbations on the data are caused by the perturbations on the conductivity, the perturbation associated with the eddy currents induced by susceptibihty, and the perturbation related to the effective magnetic charges on the boundaries. I calculate the sensitivities for both conductivity and susceptibihty over the same model used in Table 5.1. The results are given in Table 5.2. The real component of the sensitivities for susceptibihty is 4 to 6 orders of magnitude bigger, and the imaginary component of the sensitivities for susceptibility is 3 orders of magnitude greater, than those of the sensitivities for the conductivity. This means that the data are much more sensitive to the change of susceptibihty, and thus larger weighting should be applied to the model objective function for susceptibihty, to prevent the violation of linearization, and the uneven perturbation on conductivity and susceptibihty at each iteration. 5.3 Numerical results In the following synthetic examples, the earth is divided into 44 layers to a depth of 500 meters. In all inversions I use the mapped parameter m(/c) as the model parameter. The Chapter 5. 1-D Simultaneous Inversion of EM data Freq (hertz) dHz/d<r(xlO- ) r dHz/dK(xlQ- ) 3 900 -5.71 240 Real 7200 -43.6 210 56000 -122 140 128 900 -14.4 -1.0 Imag 7200 56000 -44.3 -26.8 -3.7 -5.4 Table 5.2: The comparison of the sensitivities for the conductivity and susceptibihty, from a layer of 2m thick at a depth of 48m, in a 0.01 S/m half-space whose susceptibihty is 0.1 SI unit. Figure 5.3: The result of the simultaneous inversion with 3 = 6. (a) The true and recovered conductivity models. The solid line denotes the recovered model, and the dashed line denotes the true model, (b) The true (dashed line) and reconstructed (solid line) susceptibihty models; (c) The data-misfit curve; (d) The real (solid line) and imaginary (dashed line) components of the predicted and observed data. Lines denote predicted data, and dots denote observed data. 129 Chapter 5. 1-D Simultaneous Inversion of EM data mapping parameters K i and K are set to 10~ and 10~ . The parameters a and a are 3 b 6 a K chosen as 0.02 for all of the synthetic examples. The starting and reference models for susceptibihty are 0.0 and 1 0 - 6 SI unit respectively. In the following 1-D examples, the starting and reference models for conductivity are all a 1 m S / m half-space. The same data set for the example in Figure 1 is re-inverted to simultaneously recover conductivity and susceptibihty. The use of such a relatively simple model makes it easier to investigate other aspects of the simultaneous inversion. In this inversion, I set s = 6. Figure 5.3 shows the results. The recovered conductivity model in Figure 5.3a is a good representation of the true model. The recovered susceptibihty model slightly undershoots the true model but it does recover the susceptibihty high at the right depth (Figure 5.3b). Figure 5.3c shows the misfit curve and Figure 5.3d plots the observed and predicted data. The appropriate choice of weighting parameters in equation (5.5) is important. B y increasing the parameter s I give conductivity more freedom to vary, and by reducing s, I give susceptibihty more room to vary. When s is too large, the susceptibihty will be over-depressed. On the other hand, if s is too small, large susceptibility is generated which may invalidate the linearization, and convergence difficulties may occur. I illustrate this by repeating the inversion with s = 20 and s = 0.1. Figure 5.4 shows the results of inversion with s = 20. The recovered conductivity in Figure 5.4a is a good representation of the true model. But the susceptibihty in Figure 5.4b is over-depressed into a near-surface layer, and it is not a good representation of the true model. The inversion converged after 6 iterations and fit the data to the desired misfit level. Figures 5.4c and 5.4d present the misfit curve and the data. When s = 0.1 is used in the inversion, the conductivity is recovered successfully (Figure 5.5a), but the recovered susceptibihty in Figure 5.5b overshoots the true model significantly. Beyond the ninth iteration, the inversion cannot reduce the misfit any Chapter 5. 1-D Simultaneous Inversion of EM data 130 Figure 5.4: The results from the inversion with s — 20. (a) The true and recovered conductivity models. The solid line denotes the recovered model, and the dashed line denotes the true model, (b) The true (dashed line) and reconstructed (solid line) susceptibihty models; (c) Data-misfit as a function of iterations; (d) The predicted and observed data: the solid line denotes the real component of the predicted data, and the dashed line denotes the imaginary component of the predicted data. Observed data are denoted by dots. Chapter 5. 1-D Simultaneous Inversion of EM data 131 Figure 5.5: The inversion with s = 0.1. In panels (a) and (b), the solid lines denote the recovered model, and the dashed lines denote the true model, (a) The true and recovered conductivity models; (b) The true (dashed line) and reconstructed (solid line) susceptibility models; (c) The data-misfit curve; (d) The real component of the predicted data (solid line) and the imaginary component of the predicated data (dashed line). Dots denote the observed data. further (Figure 5.5c), and the real component of the data at 900 hertz is fit to about 7%, instead of the desired level of 1%, of the amplitude of the data (Figure 5.5d). The 1-D inversions above, and other synthetic modelings, indicate that conductivity, compared to susceptibihty, is relatively insensitive to the choice of the value of parameter s. Because the data are more sensitive to the change on susceptibihty, more weight should be applied to the model objective function for susceptibihty. From my experience, s should be chosen between 2 to 20 in the simultaneous inversions. This choice is independent of the geometry of the experiments. Geological targets usually are 3-D. To determine how reliable the results are from the 1-D simultaneous inversions of 3-D field data, I tested my algorithm on 3-D synthetic data. In the following example, I inverted a 3-D synthetic data set generated from the 132 Chapter 5. 1-D Simultaneous Inversion of EM data X(m) 700 I I l I 300 500 30 5001 1 | >" 300| . 4 - - 1 I 1 _ _ I 100 r . 500 1 ao Ko •-t-i I t 300 l N 50 -+- + i . 1001 700 15 T"T~1~-1 Ko GO 0 0 T~T~1—I 100 100 700 X(m) X-Y Plan-view X-Z Cross Section Figure 5.6: The 3-D model. The background conductivity and susceptibility cr and K are 0.01 S/m and 0 SI unit respectively. The conductivity and susceptibihty in the upper prism are 0.1 S/m and 0.1 SI unit. For the lower prism, <r = 0.5S/m and K = 0.2 SI. Data are calculated at 110, 220, 440, 880, 1760, 3520, 7040, 14080, and 56320 hertz. The coil separation is 10m, and the flight height is 30m. The station interval is 25m and the line spacing is 50m. Flight lines are in the x-direction. 0 2 2 0 Chapter 5. 1-D Simultaneous Inversion of EM data 133 model shown in Figure 5.6. The background conductivity and susceptibihty are 0.01 S/m and 0 SI unit respectively. The conductivity in the upper prism is 0.1 S/m, and susceptibihty is 0.1 SI unit. In the lower prism, the conductivity is 0.5 S/m, and the susceptibihty is 0.2 SI unit. The observation height is 30 meters, and the coil separation is 10 meters. Data were calculated at 10 frequencies, ranging from 110 to 56320 hertz, with each frequency doubling the previous one. Line spacing is 50 meters and station interval is 25 meters. The weighting parameter s was set to 3. The standard deviations were 5ppm plus 10% of the data. Other parameters were kept the same as in previous 1-D examples. Figure 5.7 plots the real components of the predicted and observed data for this inversion. Panels (a), (b), and (c) show the predicted data at 110, 7040, and 56320 hertz, and panels (d), (e), and (f) plot the observed data at 110, 7040, and 56320 hertz. Figure 5.8 shows the imaginary components of the predicted and observed data at the same three frequencies. In the inversion the data were all fit into the desired misfit level, including the negative inphase data at 110 hertz. The positions of the two prisms are denoted by white rectangles. The data clearly indicate the presence of two anomaly bodies, but it is impossible to determine whether those two bodies are conductive, resistive, or permeable, and depth estimation is impossible. The simultaneous inversion recovered two conductive and permeable bodies at depth. Figures 5.9 and 5.10 show the recovered conductivity and susceptibihty models sliced in the x-direction. Figures 5.9a, 5.9b, and 5.9c show the cross-sections of the conductivity at y=250, 350, and 450 meters, and Figures 5.10a, 5.10b, and 5.10c present the corresponding cross-sections of susceptibihty. The white rectangles indicate the positions of those two prisms. The recovered conductivity represents the true model reasonably well, even though it is shallower and wider than the true model. The recovered susceptibihty is also a reasonable representation of the true Chapter 5. Predicted > > data (ppm) Observed 600 600 400 400 200 200 200 — 134 1-D Simultaneous Inversion of EM data 400 600 600 600 - 400 400 H 200 200 H f=7040 200 data (ppm) |b 400 200 600 400 600 1 - 600 ;J 400 i i 200 t f=56k| i i 200 400 X(m) 600 200 400 600 X(m) Figure 5.7: The real components of the predicated and observed data for the 3-D model. Panels (a), (b) and (c) show the predicated data at 110, 7040 and 56320 hertz. Panel (d), (e) and (f) plot the observed data at the same frequencies. 135 Chapter 5. 1-D Simultaneous Inversion of EM data Predicted Observed data (ppm) data (ppm) 126.800 600 H 600 400 i 400 116.114 105.429 94.743 84.058 > 73.372 62.687 200 H 200 |f=110 0 200 -i r 400 600 52.001 41.316 f=110 30.630 200 400 POO 1111.000 600 - 600 400 i 400 1068.089 I 1025.178 982.267 939.356 896.444 853.533 200 200 810.622 767.711 f=7040 S 724.800 200 400 200 600 400 600 400 i'Cf E 2064.00 2011.56 1959.11 1906.67 S - 1854.22 |- 1801.78 1749.33 200 |- 1696.89 1644.44 1592.00 0 200 400 X(m) 600 200 400 X(m) Figure 5.8: The Imaginary component of the predicated and observed data at 110, 7040 and 56320 hertz for the 3-D model. Panels (a), (b) and (c) show the predicated data. Panel (d), (e) and (f) plot the observed data. m) Chapter 5. 1-D Simultaneous Inversion of EM data S 136 50 ~ — .c +-> 100 - Q 150 - Conductivity (S/m) CL CD 100 200 300 400 500 600 700 0.195 0.139 0.100 0.071 0.051 0.036 0.026 0.019 0.013 0.010 X(m) F i g u r e 5.9: T h e x - z cross sections o f t h e r e c o v e r e d c o n d u c t i v i t y f r o m t h e i n v e r s i o n o f 3 - D d a t a . W h i t e r e c t a n g l e s d e n o t e t h e p o s i t i o n s o f those t w o p r i s m s , (a) T h e x - z cross s e c t i o n at y = 4 5 0 m ; ( b ) T h e x - z cross s e c t i o n at y = 3 5 0 m ; (c) T h e cross s e c t i o n at y = 2 5 0 m . Chapter 5. 1-D Simultaneous Inversion of EM data 137 o S CL CD • | 50 1 °0 150 Susceptibility (SI) 100 S - \'"'""~/ 200 300 400 500 700 600 0.158 0.140 0.123 0.105 0.088 0.070 0.053 0.035 0.018 0.000 50 S~ CL 100 CD Q 150 Susceptibility (SI) 100 200 300 400 500 700 600 ^ — IIMKNX 50 ri Q. CD [c] Susceptibility (SI) 100 200 300 400 500 600 700 0.158 0.140 0.123 0.105 0.088 0.070 0.053 0.035 0.018 0.000 0.158 0.140 0.123 0.105 0.088 0.070 0.053 0.035 0.018 0.000 X(m) F i g u r e 5.10: T h e cross sections i n t h e x - d i r e c t i o n o f t h e r e c o v e r e d s u s c e p t i b i h t y . T h e positions o f t h e prisms are denoted w i t h white rectangles. t h e cross sections at y = 4 5 0 , 350, a n d 2 5 0 m r e s p e c t i v e l y . P a n e l s ( a ) , ( b ) a n d ( c ) show Chapter 5. 1-D Simultaneous Inversion of EM data 138 0.172 0.153 600 -i 600 0.134 0.115 -g- > 400 400 0.096 0.077 0.058 200 200 0.038 z=20m z=20m 0 1 1 200 400 1— —i 600 200 1 400 id! 0.019 r 0.000 - 600 0.172 600 H 600 400 400 0.153 0.134 0.115 > w 0.096 0.077 0.058 200 200 0.038 z=40m 1 1 0 0.019 z=40m; 200 400 I I 0.000 1 0 600 200 400 600 0.172 I 0.153 600 600 0.134 0.115 E J' 1/ 400 i m- 1 0.096 400 H 0.077 0.058 200 200 0.038 f z=60m z=60m 0.019 0.000 0 200 400 X(m) 600 0 200 400 600 X(m) F i g u r e 5.11: T h e x - y p l a n - v i e w s o f t h e r e c o v e r e d c o n d u c t i v i t y a n d s u s c e p t i b i h t y . Panels (a), ( b ) a n d (c) s h o w t h e r e c o v e r e d c o n d u c t i v i t y at z = 2 0 , 4 0 , a n d 60 m e t e r s . P a n e l s ( d ) , (e) a n d (f) p r e s e n t t h e s u s c e p t i b i h t y m o d e l s at t h e c o r r e s p o n d i n g d e p t h s . Chapter 5. 139 1-D Simultaneous Inversion of EM data model. Two clearly separated susceptibility prisms are successfully recovered from the inversion. In contrast to conductivity, the recovered susceptibihty is a little deeper than the true model. Figures 5.11a, 5.11b and 5.11c show the x-y plan-views of the recovered conductivity at 20, 40, and 60 meters in depth, and Figures 5 . l i d , 5.lie, and 5.11f plot the corresponding x-y plan-views of the recovered susceptibility. At all corresponding depths, the susceptibihty is more localized than the conductivity. As a whole, the 1-D simultaneous inversion of the 3-D data set has generated conductivity and susceptibihty models which represent the true model reasonably well. Two conductivity highs and two separated susceptibility highs over the tops of the two prisms are recovered. The recovered susceptibihty has greater depth-extent than the true model, while the recovered conductivity is shallower and thinner than the true model. The recovered susceptibihty model also has better horizontal separation than the recovered conductivity. 5.4 Field data example In the following example, I invert field data collected over the Stratmat Main Zone which is located 40 km southwest of the City of Bathurst and 2 km north of the Heath Steele Mine site in Northern New Brunswick. The area is underlain by felsic to mafic and metasedimentary rocks of the Ordovician-aged Tetagouche Group, and is host to several major polymetalhc base metal sulphide deposits. A large meta-gabbro intrusion is adjacent to, and in some places, assimilates the sulphide deposits. The volcanic rocks are very resistive while the massive sulfide deposits are very conductive. Therefore E M surveys can register strong anomaly over the deposits. The magnetic minerals in the gabbroic Dyke can affect the E M data too. This Main Zone was discovered in 1957 as the result of follow-up surveys of an airborne 140 Chapter 5. 1-D Simultaneous Inversion of EM data E M anomaly. A n A E R O D A T airborne E M system was flown over this region. The coaxial data were collected at frequencies 935 and 4600 hertz, and the coplanar data were measured at 4175 hertz. The hne interval is about 200 meters, and the station spacing is about 8 meters. There are 696 stations all together. The flight height of the bird is between 20 to 50 meters. Since the coplanar data were measured at only one frequency, I needed to include the coaxial data into the inversion in order to obtain better depth resolution. Zhang et al. (1996) proposed a simple method to construct inverse algorithms for the coaxial, perpendicular, and vertical coplanar data, based upon existing inverse algorithms for the coplanar data. Following their work, I adapted my simultaneous inverse algorithm to invert the coplanar and coaxial data jointly. In carrying out the inversions, the earth was divided into 50 layers to a depth of 300 meters. The parameters ov and a were set to K 0.02. The reference model was a conductive and magnetic half-space whose conductivity and susceptibihty were 0.1 mS/m and 1 0 - 6 SI unit respectively. The starting model was a non-magnetic half-space with 1 m S / m conductivity. The error assigned to each datum was 1 ppm plus 10% of the datum strength. The parameter s was set to 3 and was fixed throughout the inversion. After 10 iterations at each station, the total chi-squared misfit level at all the stations was reduced to 4265 from the accumulative initial misfit of 187232. Figures 5.12 and 5.13 show the real and imaginary components of the predicted and observed data from the inversion. Both components of the data show a distinct anomalous high at the center of hne 12.8 km, where the Main Zone resides. The negative real component of the coplanar data is a manifestation of the existence of magnetization in this region, and justifies the necessity for a simultaneous inversion. The recovered models at each station are assembled to form a 3-D model. Figure 5.14 shows the recovered conductivity model. Panel (a) is the x-y plan-view of the Chapter 5. 141 1-D Simultaneous Inversion of EM data Predicted 10800 (ppm) Observed data (ppm) -gj g i 10600 O data V ~ — - 10400 <2U> 10200 935 hertz ] 10000 i; 12200 12400 12600 12800 13000 13200 a; _ 11000 10800 i 10200 10000 -+ 4600 hertz 12200 12400 12600 12800 13000 13200 11000 H ' ' 1 -r J 1 * 12200 12400 12600 12800 13000 13200 60.000 1 52.111 44.222 10800 36.333 S 10600 28.444 c 1 o 20.556 10400 12.667 Z 4.778 -3.111 10200 -11.000 10000 12200 -f 12400 12600 12800 13000 13200 Easting (m) F i g u r e 5.12: T h e real components 12200 12400 12600 12800 13000 13200 Easting (m) of the predicted a n d observed data. d o m i n a t i n g t h e p i c t u r e is t h e M a i n Z o n e . D a r k region Chapter 5. 1-D Simultaneous Inversion of EM data Predicted data Observed (ppm) 142 data (ppm) -|——r 1 12200 12400 12600 12800 13000 13200 11000 o 10800 — Si 10600 10200 4600 hertz 10000 b "t 12200 12400 12600 12800 13000 13200 —I—'—i—'—I—'—I— 1 12200 12400 12600 12800 13000 13200 12200 12400 12600 12800 13000 13200 Easting (m) Easting (m) F i g u r e 5.13: T h e i m a g i n a r y c o m p o n e n t s o f t h e p r e d i c t e d a n d o b s e r v e d d a t a . D a r k r e g i o n d o m i n a t i n g t h e p i c t u r e is t h e M a i n Z o n e . Chapter 5. 1-D Simultaneous Inversion of EM data 143 recovered conductivity at 30m. A conductivity high is recovered between stations at 10400m and 10500m of line 12800m. Panel (b) shows the y-z cross-section of the recovered conductivity at line 12800m. The white line in panel (a) indicates from where the section in panel (b) is taken. The recovered model suggests that the conductive body has a depth-extent of about 150 m. Other than that conductive high, the conductivity model is fairly uniform. The data from an aeromagnetic survey have been inverted to recover a 3-D model of susceptibihty (Li et al. 1996). I plot the recovered susceptibihty model from the 1-D simultaneous inversion along with that obtained from the 3-D inversion of the aeromagnetic data. Figure 5.15 shows the plan-sections of the recovered susceptibihty, at 20, 30, and 100 meters. Panels (d), (e) and (f) plot the susceptibihty recovered from the simultaneous inversion. The sources used in the E M survey are more localized than the geomagnetic field used in the magnetic surveys, so E M surveys can detect shallow and small structures, and magnetic surveys can detect deeper and larger structures. It is unclear which recovered susceptibihty model is more correct, but the major feature of the susceptibihty recovered from the simultaneous inversion is in general consistent with the result from the inversion of aeromagnetic data in panels (a), (b), and (c). The anomalies in panels (d), (e), and (f) appear as isolated peaks due to the sparse line spacing. Two major susceptibihty highs, one in the north, and one in the south, are registered by the simultaneous inversion. Those two anomaly highs coincide with the negative inphase data. But the anomaly centered at station 10450 of line 12800m in panels (a), (b), and (c) does not appear in the recovered model from the simultaneous inversion. This is because the signal level needed to resolve this "missing" susceptibihty structure is completely overwhelmed by the signal related to the eddy currents. In fact the E M responses over those neighboring susceptibihty units with the similar amplitudes are not very strong. For instance, over the susceptibihty high near station 10400 at line 12600 of the recovered Chapter 5. 144 1-D Simultaneous Inversion of EM data Conductivity Model (mS/m) 11000 10800 10600 10400 5 10200 Depth=30m 10000 12200 12400 12600 12800 13000 13200 10800 11000 Easting (m) 10000 10200 10400 10600 Easting (m) F i g u r e 5.14: AEM data, Recovered conductivity models from the simultaneous inversion of H S S (a) T h e x - y p l a n v i e w o f t h e r e c o v e r e d c o n d u c t i v i t y at z = 3 0 m ; ( b ) T h e conductivity along the section denoted w i t h the white hne i n panel (a). Chapter 5. 1-D Simultaneous Inversion of EM data 145 model from the inversion of the aeromagnetic data only causes weak responses of -1 ppm for the real component of the coaxial data at both frequencies, and -9 ppm for the real component of the coplanar data. Yet the maximum values for the real components of the coaxial data at 935 and 4600 hertz are 45 and 62 ppm respectively, and for the coplanar data the peak value is 287 ppm. In order not to overfit the data contaminated by the 3-D effects, I used 10% of the data strength as the error in the inversion, so the signal due to the susceptibihty is buried into the noise level. Figure 5.16 shows the cross sections of the recovered susceptibihty along the survey hne 12600m. The upper panel presents the results from the simultaneous inversion, and the lower panel plots the corresponding section of the recovered susceptibihty from the inversion of the aeromagnetic data. The two susceptibihty highs around stations 10200 m and 10350 m in panel (a) correspond to the anomalous highs in panel (b). The recovered model from the simultaneous inversion extents to a depth of about 150 meters, which is about the depth of investigation. Another anomalous high, which does not show in panel (b), is recovered at the right end of the section in panel (a). This anomaly is likely due to the influence of nearby 3-D susceptibihty structures on the data, and the 1-D algorithm treats those 3-D effects as if they are from a 1-D earth. The 1-D simultaneous inversion has generated 3-D images of susceptibihty and conductivity. The recovered conductivity model not only clearly indicates the existence of a conductor at the location of the Main Zone, but also suggests that the anomaly is very conductive and has a depth extent of about 150 to 200m. The recovered susceptibihty model from the 1-D simultaneous inversion is similar to that recovered from the 3-D inversion of the aeromagnetic data. 146 Chapter 5. 1-D Simultaneous Inversion of EM data 11000 11000 V 10800 e o 10600 10600 10400 10400 — 3 9 3 10200 z=100m 10000 12200 12400 12600 12800 13000 13200 12200 12400 12600 12800 13000 13200 12200 12400 12600 12800 13000 13200 12200 12400 12800 12800 13000 13200 11000 E O J I 1 1 1 1 1 1 11000 10800 10800 10600 10600 10400 H 10400 z 10200 10200 ' 1 ' 1 1 1 1 1 ' 1000012200 12400 12600 12800 13000 13200 E a s t i n g (m) 10000 1 1 1 1 1 1 1 1 r 12200 12400 12600 12800 13000 13200 E a s t i n g (m) F i g u r e 5.15: T h e x - y p l a n v i e w o f t h e r e c o v e r e d s u s c e p t i b i h t y m o d e l s f r o m t h e s i m u l t a neous inversion of A E M d a t a a n d the inversion of aeromagnetic data. Panels (a), (b), a n d (c) s h o w t h e r e c o v e r e d s u s c e p t i b i h t y f r o m t h e i n v e r s i o n o f a e r o m a g n e t i c d a t a , a n d p a n e l s ( d ) , (e), a n d (f) plot t h e s u s c e p t i b i h t y r e c o v e r e d f r o m t h e s i m u l t a n e o u s i n v e r s i o n . 147 Chapter 5. 1-D Simultaneous Inversion of EM data Line 12600 at HSS 50 1° 4 g 100 £ 150 I" 200 250 Susceptibility (mSI) 300 10000 10200 10400 10600 10800 11000 50 100 a O 250 -\ Susceptibility from mag3d (mSI) 300 10000 10200 10400 10600 10800 11000 50 45 40 35 30 25 20 15 10 5 0 Northing (m) Figure 5.16: The cross-section of the recovered susceptibihty models from the simultaneous inversion of A E M data and the inversion of aeromagnetic data over section 12600 at HSS. Panel (a) shows the recovered susceptibility from the simultaneous inversion, and panel (b) plots the susceptibihty recovered from the inversion of aeromagnetic data. 148 Chapter 5. 1-D Simultaneous Inversion of EM data 5.5 Summary Formulation of the simultaneous inverse problem requires minimization of an objective function including both conductivity and susceptibihty <J) Q4> + 7</>. B y choosing dif= <7 K ferent weighting, I can obtain different solutions. Choice of the final relative weighting between (j) and a requires additional input or knowledge on the part of the user. Since the data are much more sensitive to the changes of susceptibihty than the changes of conductivity, larger weight should be applied to the model objective function for susceptibihty, to avoid excessive perturbations on susceptibihty at any iteration. Of the two physical parameters, conductivity is more robust to the change of the weighting parameter. The numerical solution to the simultaneous inversion is more difficult than inverting for conductivity or susceptibility separately because the two model parameters are strongly coupled. Magnetic permeability \i affects the data in two ways. In the quasistatic assumption fi and a come together in the term ju>/u,a, and without prior information there is no possibility to separate them. Mathematically, independent information about susceptibihty arises only from the boundary conditions in the recursion formula for input impedances. This enlarges the scope of non-uniqueness compared to problems in which one parameter is sought. Simultaneous inversions can be useful only if the influence on data from the boundaries of the susceptibility discontinuity is greater than the noise level. Tests on both 1-D and 3-D synthetic data sets suggest that simultaneous inversion is a useful tool in providing information about conductivity and susceptibihty distributions. The results from the inversion of field data at Heath Steele Stratmat are encouraging. By inverting a single E M data set, I recovered not only the conductive deposit at the Main Zone, but also the susceptibihty structure associated with the gabbroic dyke. The negative inphase data are no longer a source of contamination - they have become an Chapter 5. 1-D Simultaneous Inversion of EM data i m p o r t a n t source for p r o v i d i n g i n f o r m a t i o n a b o u t s u s c e p t i b i l i t y distribution. Chapter 6 Approximate inversion of 3-D E M data 6.1 Introduction Geological targets are usually 3-D, therefore techniques for interpreting 3-D E M data are needed. Advances on systems of data acquisition in E M surveys make it possible to cover a large area within very short period time, so the data sets from E M surveys, especially in airborne E M surveys, could be huge. It is desirable to quickly obtain information about the buried geological targets. Quantitative information can be obtained through rigorous 3-D inversions of the E M data. However, at present, rigorous 3-D inversions are impractical because of the heavy computation involved in the forward modeling and the calculation of the sensitivities. As an alternative, various approximate solutions have been explored to glean qualitative information about the 3-D conductivity distribution. The simplest way to obtain qualitative information is to plot the data as apparent conductivities. Features seen in the apparent conductivity can sometimes be interpreted directly. However, due to the complexity of controlled-source E M problems, apparent conductivity can be a very poor representation of the true conductivity. The next step in sophistication is to construct conductivity by using an imaging technique, in which the observed data are mapped into conductivity through inverse mapping. Those inverse mappings demand much fewer computer resources, and can be used to interpret directly the geological structure during the survey. Studying and interpreting pictures presented in terms of the right physical 150 151 Chapter 6. Approximate inversion of 3-D EM data units in question (e.g. S/m) is far better than working on raw data in volt or volt/meter. Images generated from those simple mappings might be reasonable representations of the geological structure, but if nothing else, the data have been converted into a format that has the potential for offering insight about the changes of the physical property with spatial location. Those images may be used as a weighting function in a cooperative inversion of different geophysical data sets, when evidence indicating the existence of correlation among the anomaly sources. Also those images can serve as the starting model in, or even the imaging map itself may be adopted in the first several iterations of a rigorous inversion. This will save a great deal of computational time. It is noted that any algorithm which maps observed data into an element of model space is an approximate, or possibly complete, inverse mapping. The distinction between these is whether the mapped element in model space can adequately reproduce the data. Up to now, most imaging techniques solve a 1-D problem either in the frequency- Bergeron et al., 1987, Sengpiel, K . P., 1988) or the time-domain (Palacky and West, 1973; DeMoully and Becker, 1984; Zollinger et al., 1987; Nekut, 1987; Eaton and Hohmann, 1987; Macnae et al., 1991). In those methods, a nonlinear mapping is used to convert the data to a certain depth, and then a conductivity value is assigned to this depth. It is more difficult to find the counterpart to these mappings in 3-D problems. Recent work on a novel E M scattering approximation refered as the extended Born approximation (Habashy et al., 1993, and Torres-Verdin and Habashy, 1994) has provided the means to accurately simulate the electric field internal to the conductivity distribution without having to invert the large, often full, stiffness matrices that result from solving integral-equation or finite-difference simulation schemes. This approximation has recently been applied to the inversion of E M data and has demonstrated efficiency in solving geophysical exploration problems (Zhdanov and Fang, 1995). But the computation in this approximation is still so heavy that it is impractical to use this approximation Chapter 6. Approximate inversion of 3-D EM data 152 Figure 6.1: The geometry of the 3-D problem to solve large scale problems. Therefore its application has been confined to synthetic models typically consisting of relatively few model cells, and applications to field data sets have not published. In this Chapter, I use a linear inverse operator to map the data into conductivity images based on the Born approximation. The linear mapping is generated from the sensitivities calculated for a best-fit half-space. Then a 3-D model objective function is minimized subject to the constraints of a linear data objective function. This linear inverse problem is solved by using a subspace technique, which reduces the dimension of the linear inverse problem dramatically. The simplicity of this linear mapping makes it possible to solve problems involving tens of thousand model cells and hundreds of data on a workstation. With the redundancy of the data, I can still extract information about the geological targets from this first pass inversion of the data. Through this experiment I address questions such as what and how much information can be obtained by use of approximate inverse mappings on data measured at many stations, and how different the results are, compared to those from a rigorous 1-D inversion. 153 Chapter 6. Approximate inversion of 3-D EM data 6.2 The methodology The geometry for the linear mapping is shown in Figure 6.1. Both the source and the receiver are placed at a height of h above the surface of the earth. The receiver is separated from the source by a radial distance r. A regular grid system is used to divide the earth into a series regular cubic cells. This linear mapping consists of three major steps. First, a background conductivity model m^ is estimated through either 1-D inversion or apparent conductivity mapping. The perturbed data vector can then be written as d(m {n) + 8m) = d(mW) + J8m + H.O.T. (6.1) where J is the sensitivity matrix whose elements Jij quantify the change in the i t h datum with respect to unit change in the jth. model cell. By neglecting the higher order term, and setting d[m^ + m) = d , where d oh3 oba are the observations, I can introduce a data misfit objective function: <l) =\\W (Jm-d)\\\ d (6.2) d where W is an N xN matrix whose diagonal elements are the reciprocals of the estimated d error of each datum, and d = d ohs — d(m^) + Jm^ \ The third step is to solve a linear n inverse problem. The model objective function I want to minimize is <t> {™>) = a i J w (m - m ) dv + a JW-, 2 m s 0 d(m — m ) dv + 0 + «3 2 dy where m is the reference model, a a ,a , 0 w ,w ,Wy s x ly and w z 2 3 CI4 J d{m — m ) 0 dx d(m — m ) 0 W z dz dv dv. (6.3) and a control the importance of each term. 4 are the weighting functions to be supplied by the interpreter. The basis functions are chosen to be rectangular prisms of unit amplitude. Since the goal is 154 Chapter 6. Approximate inversion of 3-D EM data to find a model which minimizes a specific objective function, the inversion results should not depend upon the model parameterization. This means that a fine discretization is required and thus the number of cells will be large. Discretization of above equation yields </>m(m) = ( T O - m ) W^W (m - m ), T 0 m (6.4) 0 and W^Wm is given by W*W m = a WjW t + a WlW s 2 + a WfW x 3 y + a W?W . 4 (6.5) z In equation (6.5) W is a diagonal matrix with elements y^A^xAyAz, W = s W = y AxAzjdyW , D ^jAyAz/dxWn, x and W = Ax Ay / dzWp, where Aa:, Ay, and Az are the lengths z of the cell, and dx, dy, and dz are the distances between the central points of adjacent cells in the x-, y-, and z-directions. WD is the discrete differential operator. The inverse problem to be solved can be stated as minimize <f> = where (j) tar 4> m + B~\<f> d 4> ), tar (6.6) is the desired misfit level. Even for a small 3-D problem, the number of cells can be easily over 10 , and number 5 of data usually is greater than 10 . Thus it is impractical to solve the linear inverse 3 problem stated in above equation directly and hence the use of the generalized subspace method of Oldenburg et al. (1993) is necessary. This technique has been used to solve other linear inverse problems by L i and Oldenburg (1996) and Zhang and Oldenburg (1994). In the subspace method the model perturbation is restricted to be a linear combination of search vectors Vi, i = l,2,---,q, where q is much smaller than M, the number of model cells. Thus 8m = ^^ctiVi = Va. (6-7) 155 Chapter 6. Approximate inversion of 3-D EM data The linearized objective function in equation (6.6) can be rewritten as <f>m(™) = \\W [m^ + Va- m }\\ +(3- {\\W [d(m^) 2 m + JVa - d° ]\\ l 0 bs d 2 - 4> }. tar (6.8) Setting V(/)(a) — 0 yields Ba = c (6.9) where B and c are given by B = V (J WjW J T + 3WlW )V, T d (6.10) m and c = -BV WlW [m^ - m ] - V J WjW [d(m^) T T m T 0 d - d° }. bs (6.11) The matrix B is a q x q positive define and symmetric matrix and is easily inverted provided that q is relatively small. The steepest descent vectors associated with those two objective functions, 4> and d (f> , are used as basis vectors. The misfit objective function may be partitioned in a m variety of ways. In my work, I partition <j) according to the frequencies: d NF & = I>5, (6.12) where NF is the number of frequencies. There are also different ways to partition the model objective function. I choose to partition (j) corresponding to the four components m in equation (6.6). So the next set of vectors is v = {WlW^V^l, k (6.13) where k = 1,2,3,4. The sum of all four terms in above equation equals the steepest descent vector for V <^> . A constant vector is also included. So the total number of m m vectors used in the subspace inversion is NF + 4 + 1. 156 Chapter 6. Approximate inversion of 3-D EM data The success of a subspace approach hinges upon the spanning vectors for the activated subspace. In this Chapter, I focus on the use of the steepest descent vectors associated with the divided components of the misfit and model objective functions. Discussion on other possible choices of the spanning vectors can be found in the original papers by Gill et al. (1981) and Oldenburg et al. (1993). 6.3 Calculation of the sensitivities The key element in this linear mapping is the calculation of the sensitivities. McGillivray et al. (1994) reviewed several methods for the computation of the sensitivities. The frequency-domain Maxwell's equations are V x E = -jwiiH + J V x H = ((T + j e ) E + J W where J is the electric current, and J e m m (6.14) e ) (6.15) = — zu>/xM is the magnetic current, where M is the magnetisation. The general boundary conditions are 7](h x U ) + £(n x n x V x U ) = Q (6.16) where U stands for either E or H , 7/ and £ are constant, n is the unit normal to the boundary dD of the domain D , and Q is the appropriate electric or magnetic surface current density. If the discrete conductivity is 5>*Mir), M <r(R) = (6.17) fc=i where ipk is the basis function, then the sensitivity problem is „ v <9E x . <9H 5S = -""aS V x g - = („ + „ * ) _ , , (6 ' 18) (6.19) 157 Chapter 6. Approximate inversion of 3-D EM data and the corresponding boundary conditions are n{n x <9U ) + i{n x n x V x 8U d(T (6.20) ) = o. k The auxiliary problem associated with above sensitivity problem is 0 (6.21) V x H* = (<T + jue)G (6.22) V x G = -jufiH* + 6(r - r) where r is the observation point, and G and H* are the auxiliary electric and magnetic 0 fields. The sensitivities for the vertical magnetic field are (6.23) where G is due to a vertical magnetic dipole of strength (—itou) at r . 1 0 6.4 Approximate sensitivities The adjoint-equation method outlined in the previous section can be used to generate accurate sensitivities. However, the calculation of the electric field E and the auxiliary field G will incur considerable computational time, and this is exactly what I want to avoid. In this Chapter, I compute the electric field and the auxiliary field based upon the Born approximation. A similar technique has been used by a few authors to invert E M data. Oldenburg and Ellis (1993) used the sensitivities from the 1-D magnetotelluric inverse problem as approximate sensitivities for the 2-D problem, while Smith and Booker (1991) used 2-D electric field in their calculation of the sensitivities. In both cases, a small amount information about the 2-D structure is incorporated in the sensitivities. Farquharson (1995) used a similar technique to invert 2-D magnetotelluric data. Those approximations can lead to considerable saving on computational time, and in many cases enable the inversions to converge to an acceptable solution. However, the sensitivity Chapter 6. Approximate inversion of 3-D EM data 158 matrix does not include any information about the dependence of a measurement on the cells close to, but not immediately beneath, the station. Also, the discretization of the model is restricted to have each column of cells below a station. M y formulation does not have these restrictions. 6.4.1 The electric field and its adjoint field within a half-space In a Born approximation, the true electric field inside the scattering body is replaced with the background electric field. The choice of the background electric fields can be different but in this chapter, the electric field and its adjoint field within a half-space are used in evaluating equation (6.23). The electric field and its associated auxiliary field are in general functions of angular frequency w, coil separation r, and depth z, but for convenience I denote the electric and auxiliary fields as E and G. The electric and adjoint fields inside a layered earth have been given in Chapter 2. For a half-space, they are reduced to E = -iufi Ia 0 JO U 0 + Ui AJi(Aa)J!(Ar)dA, (6.24) and 6.4.2 The calculation of the numerical integral in Equation (6.23) In the calculation of the 3-D numerical integral in equation (6.23), both the electric and adjoint fields need to be evaluated at many locations inside each cell. To save computational time, the electric and its adjoint fields are approximated with a trilinear interpolation. Look-up tables with values of those two fields at all nodes in the model needs to be generated first. The data are put on a regular grid. Because of the symmetry of this problem, I only need to carry out the calculation on either the upper (x > 0) or Chapter 6. Approximate inversion of 3-D EM data 159 i. . ...A. n 1 in Figure 6.2: The local coordinate system for the calculation of sensitivities lower portion (a; < 0 ) of the model (Figure 6 . 2 ) . Let E = E n x x + EyTiy-hE ii , z and G = G n -\-G n x x x y + G n , where n , n , and n are y z x x y z the unit vectors in the x-, y-, and z-directions, then the real and imaginary components of the inner product are Re(E • G ) = ReE ReG - ImE ImG , +ReE ReG - ImE ImG (6.26) ReE ImG (6.27) x x x x y y y y y + ReG ImE . and Im(E •G ) = x + ReG ImE x x x + ReE ImG y Note that in 1-D Maxwell's system of equations E z y y — 0. Since G satisfies the same partial differential equation (except a scaling factor) and boundary conditions as E , G z also equals zero. If the amplitudes of E and G are denoted as E and G, then E E& x<0,y<0 (6.28) r £ f f x < 0, y > 0 160 Chapter 6. Approximate inversion of 3-D EM data G \y-R\ x < 0, y < R (6.29) a and Ey = E G y x < = G— 0 (6.30) x<0 where r = y/x + y and r = ^(y — R) + x . So equations (6.26) and (6.28) have to 2 2 2 p 2 a be calculated accordingly in regions I, II, and III in Figure 6.2: p\y\\y-R\W y -y\y-R\W Im(E • G ) 0 < p rr p Q < < — a y{y-R)W > p (6.31) R 3 R r Ta v and \y\\y-R\+x 2 Q Re(E • G ) = < Q-y\y- \+* y <0 R Qy(y-R)+x 2 o (6.32) <y<R y>R where Q = [ReE ReG - ImG ImE], and P = [Re£ ImG + ReG I m £ ] . A trilinear interpolation is used to approximate the dot product of the electric and the adjoint fields in equations (6.31) and (6.32). Let f(x,y,z) f(x,y,z) = (l-t)(l-u)(l-v)f + t(l-u)(l-v)f 1 +{l-t)u{l-v)f i = E • G , then + 2 + 7 where / i to f& are the values of f(x,y, t = x — X\ x 2 — X\ u 3 (l-t)(l-u)vf 5 +t(l - u)vf + tuvf + (1 - t)uvf , 6 tu(l-v)f (6.33) s z) at the eight nodes, and y-yi , , and v = yi - yi z-z! Z2 - Zl (6.34) 161 Chapter 6. Approximate inversion of 3-D EM data Y Z Figure 6.3: The order of the nodes for interpolation of E • G in each cell The order of the nodes in above equation is shown in Figure 6.3. The approximate sensitivities of the magnetic field are 8 (6.35) where Ax, Ay, and Az are the lengths in the x-, y-, and z-directions of the cell. In 3-D E M surveys, especially in airborne E M surveys, the data usually are collected over regions of several hundred square-kilometers. When the earth is discretized into a 3-D model, the number of the model cells can be potentially huge. Because of the small foot-print of loop-sources in E M surveys, only those cells close to the source can be illuminated. Thus there is no need to calculate the sensitivities for all cells of the model. In this thesis, I use a window to choose the model cells from which the sensitivities will be calculated. Once the sensitivity matrix g within the window is computed, the global sensitivity matrix J can be easily assembled from g. The size of the window depends upon the frequency and the geometry. The lower the frequency and the larger 162 Chapter 6. Approximate inversion of 3-D EM data the coil separation, the larger the window. For the common used frequency band and coil separations, a window of 400m by 400m is adequate. As an example, the sensitivities at 900 hertz over a O.OlS/m half-space are calculated. This half-space is divided into 80 cells in the x- and y-directions, and 12 cells in the zdirection. Each cell has the same length in the x- and y-directions (5m) and the thickness of the cells increases with respect to depth. The source is located at x = 200m, y = 200m, and the receiver is placed at x — 210m, y = 200m. The coil separation is 10m, and the observation height is 30m above the ground. Figure 6.4 shows the sensitivities for the layer between 0 and 2 meters depth. Panel (a) plots the real, and panel (b) plots the imaginary components of the sensitivities. Figures 6.6a and 6.6b present the real and imaginary components of the sensitivities due to a layer extending from 45 and 50 meters depth. The amplitude of the sensitivities drops quickly both horizontally and vertically. 6.5 Synthetic examples The data from the 2-prism model in Chapter 2 were inverted again with the approximate mapping algorithm. The model is divided into 18 by 36 cells in the x- and y-directions. Each cell has a dimension of 50m in the x-direction, and 25m in the y-direction. The thickness of the cells increases with respect to depth, to compensate for the loss in resolution. A total of 22 layers up to the depth of 200 meters was used in the inversion. This model discretization yields 15048 model cells. The error in the data is assumed to be Gaussian and independent, with standard deviation set to be 5% of the data strength. The conductivity of the best-fit half-space is 0.0105 S/m. The data objective function was partitioned according to the frequency, and the model objective function was partitioned by the four components in equation (6.3). A constant vector was also used in the inversion. The four parameters in equation (6.3) were set to oi = 2 x 10~ , a = 5 2 163 Chapter 6. Approximate inversion of 3-D EM data F i g u r e 6.4: T h e r e a l a n d i m a g i n a r y c o m p o n e n t s o f t h e sensitivities due t o a layer be- t w e e n z = 0 a n d 2 m d e p t h i n s i d e a 0.01 S / m half-space, a) T h e r e a l c o m p o n e n t o f t h e sensitivities, b) T h e i m a g i n a r y component o f the sensitivities. 164 Chapter 6. Approximate inversion of 3-D EM data F i g u r e 6.5: T h e r e a l [panel (a)] a n d i m a g i n a r y [panel (b)] c o m p o n e n t s o f t h e s e n s i t i v i t i e s d u e t o t h e l a y e r e x t e n d i n g f r o m z = 45 t o z = 5 0 m i n s i d e a O . O l S / m half-space. Figure 6.6: (a) Misfit curve; (b) Model norm as a function of the number of iterations. a = a = 1.0. The initial chi-square misfit was 5288. After 20 iterations, the misfit was 3 4 reduced to 2131. Figure 6.6a plots the misfit curve. Figure 6.6b shows the model norm as a function of iteration. Panels (a) to (c) in Figure 6.7 plot the apparent conductivity, which was obtained by fitting the data at each station to a half-space, at 900, 7200, and 56000 hertz. The apparent conductivity indicates the existence of two anomalous bodies. But it is impossible to tell whether those two bodies are conductive or resistive, let alone their buried depths. Figures 6.7d to 6.7f show the x-y plan-views of the model recovered from the inverse mapping, at 20, 45, and 80 meters depth. The white boxes indicate the positions of those two prisms. At z=20m, the inversion successfully recovered the upper conductive prism. The plan-view at z=45m indicates a second conductive prism. The recovered model at z=80m forms a single anomalous body centered at x=300m, which is not a good image of the true model. Figures 6.8a to 6.8c show the x-z cross sections of the recovered model at y=250, 350, and 450 meters. At y=250m, the upper prism is well defined by the inversion. A t section y=350, the inversion recovered an anomaly high, extending approximately from station Chapter 6. Approximate inversion of 3-D EM data 166 150m to station 500m. The recovered model at y=450m also indicates the existence of a conductive body, whose center is shifted towards the upper prism. Compared to the recovered model from 1-D inversion of the same data set (Figure 6.8d), the reconstructed conductivity from the 3-D linear mapping defines the thickness of the conductive bodies better. On the other hand, the recovered model from the 1-D inversion is more compact. The synthetic data example shows that, even though there is blurring and distortion of the conductivity blocks, there is much information about the 3-D conductivity structure that is evident in the images. 6.6 Field data example The most crucial test for any geophysical inverse algorithm is the inversion of field data sets. I use the linear mapping to image the D I G H E M data at M t . Milligan, over a 1.3 km by 1 km region shown in Figure 6.9. More information about the background of the data set can be found in Chapter 2, where a section of the data set (Y9600) was inverted. Figures 6.9a and 6.9b show the real and imaginary components of the coplanar data at 56000 hertz. Panels (c) and (d) plot the data at 7200 hertz, and panels (e) and (f) show the data at 900 hertz. Negative inphase data at 900 and 7200 hertz are observed at the center of the region. Those negative inphase data result from the strong magnetization in this region, and cannot be explained by any purely conductive structure. Only the quadrature phase of the data was therefore used in the inversion, since it is less affected by the magnetization. The hne spacing is about 100m and the station interval is about 25m. So the number of stations is 11x13 = 143, and the number of data is 143 x 3 = 429. A window of 400m by 400m was used to calculate the sensitivities. In order to apply the window in the inversion, the region is enlarged by 200m in each direction. The sizes of the model are chosen to be 50m in the x- and y-directions. The earth is divided into 25 layers 167 Chapter 6. Approximate inversion of 3-D EM data Recovered Conductivity (mS/m) Apparent Conductivity (mS/m) 600 i- 7 0 . 0 98.0 - 63.1 88.7 - 56.2 79.4 - 49.3 70.1 fill ) III - 42.4 400 - 35.5 - 28.6 200 0 - - 21.7 - 14.8 - 7.9 L 200 400 200 400 600 _ 1.0 z=20m 60.8 51.5 42.2 32.9 23.6 [dj n — 14.3 5.0 0 200 400 200 400 600 600 400 200 600 400 H 200 E a s t i n g (m) E a s t i n g (m) Figure 6.7: The x-y plan-views of the recovered model and apparent conductivity. Panels (a) to (c) plot the apparent conductivity obtained from data at 56000, 7.2k and 900 hertz. Panels (d) to (f) show the recovered model at z=20, 45, and 80 meters respectively. Chapter 6. Approximate inversion of 3-D EM data 168 R e c o v e r e d Conductivity (mS/m) o 50 100 150 JS y=250m \ I / 1 a . f 100 200 300 400 500 600 700 700 100 200 300 200 300 400 500 600 400 500 600 0 50 100 150 y=350m] 100 Easting 700 373.00 336.00 299.00 262.00 226.00 189.00 152.00 115.00 78.40 41.50 4.73 279.0 251.0 224.0 196.0 169.0 142.0 114.0 86.7 59.3 31.8 4.4 80.30 72.60 64.90 57.20 49.40 41.70 34.00 26.30 18.60 10.90 3.15 131.00 119.00 107.00 94.80 82.60 70.40 58.20 46.00 33.80 21.60 9.44 (m) Figure 6.8: The x-z cross-sections of the recovered model at y=250, 350, and 450m [panels (a) to (c)]. Panel (d) shows the recovered model from 1-D inversion. 169 Chapter 6. Approximate inversion of 3-D EM data in the vertical direction. So the total number of model cells is 32 x 28 x 25 = 22400. The parameters in equation (6.3) were chosen as a = 2 x 10~ , and a = a = a = 0.1. The 5 x 2 3 4 number of basis vectors was 8, and the vectors were the same as those in the previous inversion of the 3-D synthetic data. The error was assumed to be 10% of the data strength. The starting model is the best-fit half-space of 1.667 mS/m, which gives an initial chi-square misfit of 2.277 x 10 . 5 A crucial aspect in an inverse mapping is to decide how well the data should be fit. I first set the target misfit as 4 x 10 . After 20 iterations, the chi-square misfit 4 was converged to the specified target misfit. Figure 6.10 displays the recovered model at sections Y9200, Y9400, and Y9600. Two conductive highs are observed in panel (a), around stations at x=400 and 600m. Those two conductive highs become further and further away from each other at sections Y9400 and Y9600 (panels (b) and (c)). A resistive layer over the top of those two conductive anomalies is recovered from the inversion. Panel (d) presents the recovered conductivity from the 2-D inversion of D C resistivity data at section Y9600 (Oldenburg et al., 1996). The numerical range of the recovered model from the approximate 3-D inversion is wider than the model in panel (d). The conductive zones at the left upper corner in panels (c) and (d) have similar numerical values. In panel (c), the conductive anomaly around station at x=900m corresponds to the conductive zone at the right end of the section in panel (d), but dips towards left instead of right. The geological cross-section at section Y9600, which was provided by Delong et al. (1995), was incorporated into Figures 6.10c and 6.10d. The major features of the geological cross-section include: (1) the monzonite stock, which is truncated on the left by the Harris Fault, (2) the arcuate feature, refered as the Rainbow Dyke, extending from the stock to the surface on the right, and (3) trachyte Dyke between the Rainbow Dyke and the central monzonite unit. M B X Stock, and Rainbow Dyke. High conductivities Chapter 6. Approximate inversion of 3-D EM data 170 are observed at the left of the Harris Fault. To the east of the Harris Fault, low conductivity was recovered at the top 50m of the recovered model. The inverse mapping also recovered a conductive zone in the volcanics to the east where the Rainbow Dyke reaches the surface. Figure 6.11 shows the x-y plan-view of the recovered model at z=10, 30, 50, 80, 100, and 120m. A resistive zone is circled by a conductive ring (panels (a) to (c)). The resistive core disappears at depth (panels (d), (e) and (f)). Next I set the target misfit level to 2.7 x 10 , to investigate the effect of different target 4 misfit level on the inverse mapping. It took 24 iterations for the inversion to converge. The recovered conductivity image is shown in Figure 6.12. The images look similar to those in the previous figure, but the amplitude of the conductivity recovered from the inverse mapping becomes unreasonably high, and small artificial structures start to build up. Those are the indication that the data have been over-fit, and that the constructed model has become unreliable and should not be used in interpretation. The image obtained from the hnear inverse mapping has been useful in delineating fault regions. The size of the problem solved in the field data example indicates that this hnear inverse mapping may be used as a practical tool to extract information about 3-D conductivity distributions. 6.7 Conclusions A hnear inverse mapping based on Born approximation is developed in this Chapter, and it can be used to map E M data into 3-D images of conductivity distribution. This inverse mapping is obtained by linearizing a nonlinear objective function, and then solving a hnear inverse problem. To handle large scale problems, a subspace method is incorporated into this inverse mapping. Chapter 6. Approximate inversion of 3-D EM data Imag(ppm) Real (ppm) 10.0 171 10.0 ~i 1 r 12.4 12.6 12.8 13.0 13.2 13.4 12.4 12.6 12.8 13.0 13.2 13.4 Easting (km) Easting (km) 10.0 12.4 12.6 12.8 13.0 13.2 13.4 12.4 12.6 12.8 13.0 13.2 13.4 Easting (km) Easting (km) 12.4 12.6 12.8 13.0 13.2 13.4 12.4 12.6 12.8 13.0 13.2 13.4 Easting (km) Easting (km) Figure 6.9: D I G H E M data at M t . Milligan. Panels (a), (c), and (e) show the real components of the coplanar data at 56000, 7200, and 900 hertz respectively. Panels (b), (d), and (f) show the corresponding imaginary components. 172 Chapter 6. Approximate inversion of 3-D EM data Recovered Conductivity (mS/m) . — , 12.3 _ . , R 12.4 12.5 12.6 12.7 12.8 12.9 13.0 13.1 13.2 13.3 13.4 13.5 Easting (km) S" 150 12.3 12.4 12.5 12.6 12.7 12.8 12.9 13.0 13.1 13.2 13.3 13.4 13.5 Easting (km) 12.3 "1 1 ! 1 i 1 1 "1 1 r 12.4 12.5 12.6 12.7 12.8 12.9 13.0 13.1 13.2 13.3 13.4 13.5 Easting (km) J 0 L 1 50 £ 100 J " 150 200 MBX Stock Rainbow •• ' Dyke _d_ — 12.3 12.4 12.5 12.6 12.7 12.8 12.9 13.0 13.1 13.2 13.3 13.4 Easting (km) 13.5 100.000 50.100 25.100 12.600 6.310 3.160 1.580 0.794 0.398 0.200 0.100 100.000 50.100 25.100 12.600 6.310 3.160 1.580 0.794 0.398 0.200 0.100 100.000 50.100 25.100 12.600 6.310 3.160 1.580 0.794 0.398 0.200 0.100 15.800 11.200 7.940 5.620 3.980 2.820 2.000 1.410 1.000 0.708 0.501 Figure 6.10: The X - Z cross section of the recovered conductivity, with overlayed geology (white lines), from the 3-D approximate inversion of D I G H E M data at M t . Milligan. Panels (a), (b) and (c) plot the sliced model at Y9200, Y9400, and Y9600. Panel (d) shows the recovered conductivity from 2-D inversion of D C resistivity data over section Y9600. 173 Chapter 6. Approximate inversion of 3-D EM data Recovered Conductivity (mS/m) -i 12.4 r 12.6 12.8 13.0 13.2 13.4 I I 12.4 12.6 i 12.8 l l 13.0 13.2 13.4 7 Easting (km) Easting (km) Figure 6.11: The X - Y plan-view of the recovered conductivity from the 3-D approximate inversion of D I G H E M data at M t . Milligan. Panels (a) to (f) plot the sliced model at z=10, 30, 50, 80, 100, and 120m. 174 Chapter 6. Approximate inversion of 3-D EM data R e c o v e r e d Conductivity (mS/m) o g ! . . . . . . . . . . 50 } T 100 CL G Q J y=9.2km _y a i 200 12.3 12.4 12.5 12.6 12.7 12.8 12.9 13.0 13.1 13.2 13.3 y 13.4 13.5 Easting (km) 12.3 12.4 12.5 12.6 12.7 12.8 12.9 13.0 13.1 13.2 13.3 13.4 13.5 Easting (km) 12.3 12.4 12.5 12.6 12.7 12.8 12.9 13.0 13.1 13.2 13.3 13.4 13.5 Easting (km) - £ 100 8" 150 200 H?rri4 •K3? M B X Stock y=9.6km | -: Rainbow Dyke ) / d 1 12.3 12.4 12.5 12.6 12.7 12.8 12.9 13.0 Easting (km) 13.1 13.2 13.3 13.4 13.5 50100.0 13500.0 3630.0 977.0 263.0 70.8 19.1 5.1 1.4 0.4 0.1 50100.0 13500.0 3630.0 977.0 263.0 70.8 19.1 5.1 1.4 0.4 0.1 50100.0 13500.0 3630.0 977.0 263.0 70.8 19.1 5.1 1.4 0.4 0.1 15.800 11.200 7.940 5.620 3.980 2.820 2.000 1.410 1.000 0.708 0.501 Figure 6.12: The X - Z cross section of the best-fit conductivity, with overlayed geology (white lines), from the 3-D approximate inversion of D I G H E M data at M t . Milligan. Panels (a), (b) and (c) plot the cross-sections at Y9200, Y9400, and Y9600. Panel (d) shows the recovered conductivity from 2-D inversion of D C resistivity data over section Y9600. 175 Chapter 6. Approximate inversion of 3-D EM data The c h o i c e o f t a r g e t misfit is i m p o r t a n t for t h i s i n v e r s e m a p p i n g . T o keep the lin- e a r i z a t i o n v a l i d , t h e d a t a misfit s h o u l d n o t b e r e d u c e d t o o m u c h . T h e a p p r o p r i a t e t a r g e t misfit level m a y be determined b y e x a m i n i n g the a m p l i t u d e of the recovered c o n d u c t i v i t y i m a g e s f r o m s e v e r a l test r u n s . B e c a u s e o f t h e B o r n a p p r o x i m a t i o n , c a u t i o n is n e e d e d w h e n a p p l y i n g t h i s i n v e r s e m a p p i n g t o d a t a c o l l e c t e d f r o m regions w i t h l a r g e c o n d u c t i v i t y contrast. The output from this linear m a p p i n g m a y be geologically a n d geophysically inter- p r e t a b l e . I f t h e i n v e r s e m a p p i n g c a n b e done r e a s o n a b l y fast t h e n i t w i l l p r o v i d e p r o m p t f e e d b a c k i n field s u r v e y s . T h i s t e c h n i q u e allows t h e d a t a t o b e p r e s e n t e d i n t h e r i g h t p h y s i c a l u n i t s . I f t h e r e are a r t i f i c i a l s t r u c t u r e s , t h e i n t e r p r e t e r w i l l h k e l y c o m e t o r e c o g n i z e a n d a c c o u n t for t h e m i n a n y p r e h m i n a r y i n t e r p r e t a t i o n . T h e o u t p u t f r o m t h i s h n e a r m a p p i n g m a y also b e u s e d as a w e i g h t i n g f u n c t i o n i n t h e i n v e r s i o n o f a n o t h e r g e o p h y s i c a l d a t a set. T h i s c a n b e a first pass c o o p e r a t i v e i n v e r s i o n . T h e i m a g e o b t a i n e d f r o m t h e h n e a r m a p p i n g c a n also b e u s e d as a s t a r t i n g m o d e l for s u b s e q u e n t i n v e r s i o n s . The r e s u l t s f r o m t h e s y n t h e t i c d a t a test are e n c o u r a g i n g . T h e r e s u l t a n t i m a g e s h a v e provided m u c h more information about the depth extension a n d a m p l i t u d e of the 3-D c o n d u c t i v i t y d i s t r i b u t i o n t h a n apparent conductivity. T h e test o n t h e field d a t a set also g e n e r a t e d useful i m a g e s w h i c h n o t o n l y agree w i t h t h a t o b t a i n e d f r o m 2 - D r i g o r o u s i n v e r s i o n o f D C r e s i s t i v i t y d a t a , b u t also are g e o l o g i c a l l y m e a n i n g f u l . T h e e x i s t e n c e o f t h e t h r e e m a j o r g e o l o g i c a l u n i t s is o b v i o u s i n t h e r e c o n s t r u c t e d c o n d u c t i v i t y i m a g e . Chapter 7 Summary A E M source induces secondary currents in the earth, which in turn generate secondary magnetic fields which can be measured by E M receivers. The E M response depends upon the frequency, the conductivity, the magnetic susceptibihty, and the geometry. B y carrying out inversions of the measured data, information about subsurface structures can be obtained. Geological structure is usually 3-D, but the multi-dimensional rigorous inversion of E M data is currently too computationally demanding to.be carried out on a routine basis, even using the state-of-the-art computer technology. Therefore I focus on solving 1-D inverse problems. The aim of this thesis is to develop techniques for inverting controlled-source E M data to extract information about conductivity and susceptibihty. Different orientations of the source probe the earth somewhat differently while different orientations of the receiver measure different components of the earth response. When only a limited number of frequencies are available, the data from different coil configurations provide complementary information about the conductivity. E M data have long been routinely collected with different coil systems. Flexible, robust and stable algorithms to invert the data from the coaxial, perpendicular and vertical coplanar E M systems are long overdue. In Chapter 2, I described an iterative, linearized inversion procedure for inverting frequency-domain E M data collected with the horizontal coplanar, vertical coplanar, coaxial, and perpendicular coil configurations for a 1-D conductivity model. This procedure is typical of the inversion strategy that one would like to apply to the data in time-domain, and to multi-dimensional E M data. Up to now, the rigorous 176 Chapter 7. Summary 177 inversion algorithm developed in Chapter 2 is the only one in the literature which can invert data from different coil systems. This technique can be potentially extended into the time-domain, and can also be used to construct various approximate 3-D inversion algorithms. A hnear relationship among the horizontal coplanar, vertical coplanar, and coaxial data has been established. This relationship only holds over a 1-D earth, and therefore it can be used to signal the violation of the 1-D assumption. It can also be used as a potential mapping tool, to locate the regions of high conductivity variation. In a 1-D environment, inversions of data from different coil systems can generate similar results. But for those systems such as EM-31 and EM-38 which operate at only a very limited number of frequencies, joint inversion of 1-D data may be potentially beneficial. If, however, signal-to-noise ratio of the data is high and data for each coil configuration are measured at many frequencies, then the joint inversion of 1-D data may not provide extra information. In a 3-D enviroment, different coil systems couple with the targets differently, and thus E M data measured with different coil systems carry independent information about the anomalous bodies. Consequently it should be beneficial to use 3-D algorithms to jointly invert those E M data from different coil systems. However, because of the limitations of computer resources, 1-D algorithms in many cases are used to interpret 3-D data. The major difficulty with 1-D joint inversions of 3-D data is that 3-D structure affects data from different coil systems in different ways, so those data sets may be incompatible with the 1-D assumption and the ascribed error. Therefore it is unclear if 1-D joint inversion of 3-D data collected with different coil systems can generate a better result. Through synthetic and field data examples it is found that in 1-D joint inversions, the 3-D data cannot be fit to the same degree as in the inversions of individual data sets, but the joint inversion can be potentially useful because the recovered models carry the characteristics Chapter 7. Summary 178 of models recovered from the inversions of data from different coil configurations. According to inverse theory, model space can be divided into activated and null regions. The activated region is the part of the model space illuminated by the data. The larger the activated region, the more information about the model can be obtained. When only a limited number of frequencies or time channels are available, the activated region can be enlarged by measuring the E M response at as many positions as possible. Joint inversion of surface and borehole E M data is such an example. Surface data can illuminate one part of the model space, and borehole data can activate another part of the model space. These two activated subspaces in general do not completely overlap with each other, and the combination enlarges the activated portion of the model space. In Chapter 3, I extended the work in Chapter 2 to invert transient E M data from both surface and borehole configurations. The algorithm developed in this chapter is probably the only code which can carry out a rigorous inversion of large loop T E M borehole data. Surface T E M data have better signal-to-noise level at early times but the borehole data enjoy higher signal-to-noise ratio at late times. From the synthetic data examples it is concluded that inversion of surface data can delineate the near surface structure better, while inversion of borehole data can reconstruct the structures at depth better. Joint inversion of these surface and borehole data has enlarged the activated model space, and therefore the results from inversions of individual data sets are improved significantly. In practice, however, 3-D effects exist, which can cause difficulties for the 1-D joint inversion. The field data example is inconclusive but it did show that the use of borehole data in the inversions can enhance the structures at depth. The effect of magnetic susceptibihty on E M data has long been appreciated. A classic example is the negative inphase data measured with a typical frequency domain airborne electromagnetic system. These negative data can not be explained by any purely conductive model, and have been routinely excluded from the inversions in searching for Chapter 7. 179 Summary conductivity. In Chapter 4,1 developed an inversion algorithm for extracting information, when the conductivity is known, about the distribution of susceptibility from the 1-D inversion of E M data. The susceptibihty enters the problem through two avenues. In the quasi-static assumption fi and fj come together in the term ju>fia, and without prior information there is no possibility to separate them. Mathematically, independent information about susceptibihty arises only from the boundary conditions in the recursion formula for input impedances. The same technique as discussed in Chapter 2 was used to compute the sensitivities, but an extra boundary term was introduced into the solution. That boundary term represents the influence of effective magnetic charges accumulated on boundaries of susceptibihty discontinuity. This algorithm has generated very useful information about the distribution of susceptibihty from the inversion of a field data set. The technique is new, and can be useful when the conductivity is relatively uniform and resistive. The shortcoming is that conductivity has to be specified before hand, and in many field situations this is unknown. Therefore, the theoretical value of this technique is much greater than its practical value. A few authors (Qian et al., 1996; Beard et al., 1996) have been inspired by the work in this chapter, and have published their work subsequently. As a natural extension to the work in Chapter 4, I attacked simultaneous inversion in Chapter 5, which can recover 1-D susceptibihty and conductivity simultaneously from the inversion of a single E M data set. This algorithm eliminates the need to specify either conductivity or susceptibihty a priori, as required in an individual inversion. The major difficulty is the specification of an objective function that includes both conductivity and susceptibihty. I used the weighted sum of individual model objective functions for conductivity and susceptibihty to form a global model objective function, and minimized that global model objective function subject to fitting the single data set. The choice of Chapter 7. Summary 180 the appropriate weight is important. Conductivity is dominant in the forward modeling because it varies several orders of magnitude in amphtude, but the data are more sensitive to the changes of the susceptibihty. In order to avoid the violation of linearization, it is necessary to apply more weights to the model objective function for susceptibihty. I tested the algorithm on field data collected at Heath Steele Stratmat, and found that the results provide useful information about the geology of the region. B y inverting a single E M data set, I recovered not only the conductive deposit at the Main Zone, but also the susceptibihty structure associated with the gabbroic dyke. The simultaneous inversion technique to recover conductivity and susceptibihty from E M data is new in the literature. The impact of this inversion lies in that from now on, negative inphase A E M data will not be considered as contamination due to magnetic susceptibihty, and can be included in the inversion to provide useful information about the distribution of magnetic susceptibihty through inversions. Due to the hmitation of computer resources, most work on the inversions of controlledsource E M data is confined to the 1-D framework. Geological targets are usually 3-D; therefore 3-D inversion is highly desirable. Approximate 3-D inversion provides a possible solution to this stalemate. In Chapter 6, I used a Born approximation to construct a 3-D hnear mapping, which can generate a 3-D image of the conductivity from horizontal coplanar E M data. This mapping was posed as a hnear inverse problem, and it was solved with a subspace method. This algorithm was applied to synthetic and field data sets. It has successfully handled a relatively large scale problem that consisted of more than 22,000 cubic cells and more than 400 data. In 1993, when this study was finished, that was considered to be a large problem. I believe that this technique is still valuable. The resultant model can be used directly in the interpretation, and can be also used as a starting model for other linearized inversions. The hnear mapping itself can be potentially used in the first several iterations of a rigorous inversion to save computational time. Chapter 7. Summary 181 In s u m m a r y , I have presented techniques that c a n generate i n f o r m a t i o n about conduct i v i t y a n d susceptibihty from the rigorous simultaneous inversions of controlled-source EM d a t a f r o m different c o i l s y s t e m s for t h e first t i m e . T h e s e t e c h n i q u e s are p o t e n - t i a l l y useful t o t h e m i n i n g i n d u s t r y a n d i n o t h e r places w h e r e m a g n e t i c s u s c e p t i b i h t y is a c o n c e r n . I h a v e also r e c o n s t r u c t e d c o n d u c t i v i t y f r o m t h e i n v e r s i o n s o f t r a n s i e n t E M b o r e h o l e d a t a , w h i c h c a n b e s t i l l c o n s i d e r e d as a n e w t e c h n i q u e . T h o s e 1-D i n v e r s i o n m e t h o d s d e v e l o p e d i n t h i s thesis h a v e p r o v i d e d necessary p h y s i c a l i n s i g h t s n e e d e d t o a t t a c k h i g h e r d i m e n s i o n a l p r o b l e m s . T h e 3 - D l i n e a r m a p p i n g c a n b e p o t e n t i a l l y useful i n the interpretation a n d inversion of controlled-source E M data. References [1] Anderson, W . L . , 1979, Computer Program: Numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering: Geophysics, 44, 1287-1305. [2] Beard, L.P. and Nyquist, J.E., 1996, Inversion of airborne electromagnetic data for magnetic permeability: Expanded Abstract of the 67th S E G annual meeting, p940-943. [3] Belluigi, A . , 1949, Inductive coupling of a homogeneous ground with a vertical coil: Geophysics, 14, 501-507. [4] Bergeron, C. J . , Jr., Loup, J . W . , and Michel, G . A . II, 1987, Application of the modified image method to inversion of synthetic active electromagnetic bathymetric data: Geophysics, 52, 794-801. [5] Bhattacharyya, B . K . , 1980, A generalized multi-body model for inversion of magnetic anomalies: Geophysics, 45, 255-270. -1963, Electromagnetic fields of a vertical magnetic dipole placed above the earth's surface: Geophysics, 28, 408-425. -1955, Electromagnetic induction in a two-layer earth: J . Geophys. Res., 60, 279288. [6] Backus, G . E . and Gilbert, J.F., 1968, The resolving power of gross earth data: Geophys. J . Roy, Astr. Soc, 16, 247-276. 182 References 183 [7] Constable, S.C., Parker, R . L . and Constable, C . G . , 1987. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics, 52, 289-300. [8] DeLong, R. C , Godwin, E. I., Harris, M . H . , Caira, N . M . and Rebagliati, C. M . , 1991, Geology and Alteration at the Mountain Milligan Gold-Porphyry Deposit, Central British Columbia: In Geological Fieldwork 1990, British Columbia Ministry of Energy, Mines and Petroleum Resources - Geological Survey Branch, Paper 19911, 199-205. [9] DeMoully, G. T., and Becker, A., 1984, Automated interpretation of airborne electromagnetic data: Geophysics, 49, 1301-1302. [10] Dey, A . and Ward, S. H . , 1970, Inductive sounding of a layered earth with a horizontal magnetic dipole: Geophysics, 35, 660-703. [11] Doll, H . G . , 1949, Introduction to induction logging and application to logging of wells drilled with oil base mud: J.Pet. Tech., 1, 148-162. [12] Duckworth, K . , 1970, Electromagnetic depth sounding applied to mining problems: Geophysics, 35, 1086-1098. [13] Dyck A . V . , 1975, Electrical borehole methods applied to mineral prospecting: in Dyck A . V . , Ed., Borehole geophysics applied to metallic mineral prospecting: A review: Geol. Surv. Can. Paper 75-31, 13-29. [14] Eaton, P. A . and Hohmann, G.W., 1987, A rapid inversion technique for transient electromagnetic soundings: Physics of the Earth and Planetary Interiors, 53, 384404. .References 184 [15] Elliot, C. L . , 1961, A electromagnetic drill hole instrument for the detection of conductive sulfide bodies: Presented at the 31th S E G Ann. Internat. Mtg. - 1966 Electromagnetic method and apparatus of geophysical exploration: Canadian patent 743,665. [16] Ellis, R. G . , 1995, Joint 3D E M inversion: Abstract for 3D electromagnetics workshop at Schlumberger-Doll Research, 307-323. [17] Farquharson, C.G., 1995, Approximate sensitivities for E M inversion: PhD thesis, the University of British Columbia, Canada. [18] Farquharson, C . G . and Oldenburg, D.W., 1993, Inversion of time-domain electromagnetic data for a horizontally layered earth: Geophys. J . Internal., 114, 433-442. [19] Fraser, D . C , 1973, Magnetite ore tonnage estimates from an aerial electromagnetic survey: Geoexploration, 11, 97-105. - 1981, Magnetite with a multi-coil airborne electromagnetic system: Geophysics, 46, 1579-1593. [20] Frischknecht, F . C . , and Raab, P.V., 1984, Time-domain electromagnetic soundings at the Nevada Test Site, Nevada: Geophysics, 49, 981-992. [21] Fullagar, P . K . , 1984, A uniqueness theorem for horizontal loop electromagnetic frequency soundings: Geophys. J . Roy. Astr. Soc, 77, 559-566. [22] Fullagar, P.K. and Oldenburg, D.W., 1984, Inversion of horizontal loop electromagnetic frequency soundings: Geophysics, 49, 150-164. [23] Fuller, J . A . and Wait, J.R., 1972, High-frequency electromagnetic coupling between small coplanar loops over an inhomogeneous ground: Geophysics, 37, 997-1004. References 185 [24] Gill, P.E., Murray, W . , and Wright, M . H . , 1981, Practical optimization, Academic Press, New York, 401pp. [25] Glenn, W . E . , Ryu, J . , Ward, S.H., Peeples, W . J . and Phillipps, R . J . , 1973, The inversion of vertical magnetic dipole sounding data: Geophysics, 38, 1109-1129. [26] Glenn, W . E . and Ward, S.H., 1976, Statistical evaluation of electrical sounding methods: Geophysics (Part I), 41, 1207-1221. [27] Gomez-Trevino, E . , and Edwards, R . N . , 1983, Electromagnetic soundings in the sedimentary basin of southern Ontario-a case history: Geophysics, 48, 311-326. [28] Gorden, A . N . , 1951, The field induced by a oscillatory magnetic dipole outside a semi-infinite conductor, Q.J. Mech. Appl. Math., 4, 106-115. [29] Gradshteyn, S., Ryzhik, I.M., 1965, Table of Integrals, Series, and Products, 4th, 172, Academic Press Inc, New York and London. [30] Guillen, A. and Menichetti, V . , 1984, Gravity and magnetic inversion with minimization of a special functional, Geophysics, 49, 1354-1360. [31] Gupta Sarma, D., Maru., U . M . , and Varadarajan, G . , 1976, A n improved pulse transient electromagnetic system for locating good conductors: Geophysics, 41, 287-299. [32] Habashy, T . M . , Groom, R.W., and Spies, B.R., 1993, Beyond the Born and Rytov approximations: A non-linear approach to electromagnetic scattering: J . Geophys. Res., 98, 1795-1775. [33] Kaufman, A . A . , and Keller, G.V., 1983, Frequency and transient sounding: Elsevier scientific Publ. Co. References 186 [34] Keller, G. V., 1990, Rock and Mineral Properties, in Nabighian, M. N., Ed., Electromagnetic Methods in Applied Geophysics: Investigations in geophysics, Vol. 3, Soc. Expl. Geophys. publications, 13-51. [35] Keller, G. V. and Frischknecht, F.C., 1966, Electric methods on geophysical prospecting, Pergamon Press, New York, NY, 519p. [36] Last, B. J. and Kubik, K., 1983, Compact gravity inversion, Geophysics, 48, 713721. [37] Li, Y. and Oldenburg, D. W., 1996, 3-D inversion of magnetic data, Geophysics, 62, 394-408. [38] Li, Y., Oldenburg, D. W., Farquharson, C , and Shekhtman, R., 1996, Inversion of geological data at Heath Steele Stratmat, JACI annual report, Appendix A. [39] Macnae, J. C , and Lamontagne, Y., 1987, Imaging quasi-layered conductive structures by simple processing of transient electromagnetic data: Geophysics, 52, 545554. [40] Macnae, J. C , Smith, R., Polzer, B. D., Lamontagne, Y., and Khnkert, P. S., 1991, Conductivity-depth imaging of airborne electromagnetic step-response data: Geophysics, 56, 102-114. [41] Mallick, K., 1978, A note on decay pattern of magnetic field and voltage response of conducting bodies in electromagnetic time-domain system: Geoexploration, 16, 303-307. [42] Malmqvist, D., 1965, A numerical calculation of the electromagnetic field from a vertical and horizontal magnetic dipole above a homogeneous earth: Geoexploration, 3, 175-227. References 187 [43] McGillivray, P.R. and Oldenburg, D.W., 1990, Methods for calculating Frechet derivatives and sensitivities for the non-linear inverse problem: Geophys. Prospect., 38, 499-524. [44] Nabighian, M.N., Ed, 1988, Electromagnetic methods in applied geophysics, Volume 1, Theory: Soc. Expl. Geophys. - 1979, Quasi-static transient response of a conducting half-space: An approximate representation: Geophysics, 44, 1700-1705. [45] Nekut, A.G., 1987, Direct inversion of time-domain electromagnetic data (short note): Geophysics, 52, 1431-1435. [46] Newman, G. A. and Alumbaugh, D. A.,1995, 3D massively parallel electromagnetic inversion: Abstract for 3D electromagnetics workshop at Schlumberger-Doll Research, 287-306. [47] Newman, G. A. and Alumbaugh, D. A.,1995, Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences: Geophys. Prospct., 43, 1021-1042. [48] Newman, G. A., Hohmann G. W., and Anderson W. L., 1986, Transient electromagnetic response of a three-dimensional body in a layered earth: Geophysics, 51, 1608-1627. [49] Nekut, A. G., 1987, Direct inversion of time-domain electromagnetic data: Geophysics, 52, 1431-1435. [50] Noakes, J.E., 1951, An electromagnetic method of geophysical prospecting for application to drill holes: PhD thesis, Univ. of Toronto. References 188 [51] Oldenburg, D.W. and Ellis, R.G., 1993, Efficient inversion of magnetotelluric data in two dimensions, Phys. Earth Planet. Int., 81, 177-200. [52] Oldenburg, D. W., Li, Y., and Ellis, R.G., 1996, Inversion of geophysical data over a copper gold porphyry deposit: a case history for Mt. Milligan: Accepted by Geophysics. [53] Oldenburg, D.W., McGilhvray, P.R., and Ellis, R.G., 1993, Generalized subspace methods for large scale inverse problems: Geophys. J. Int., 114, pl2-20. [54] Palacky, G.J. and West, G.F., 1973, Quantitative interpretation of Input AEM measurements: Geophysics, 38, 1145-1158. [55] Parker, R. L., 1977, Understanding inverse theory: Ann. Rev. Earth Planet. Sci., 5, 35-64. [56] Patra, H.P., 1970, Central frequency sounding in shallow engineering and hydrogeological problems: Geophys. Prospect., 18, 236-254. [57] Qian, W., Gamey, J., Lo, B., and Holladay, J.S., 1996, AEM apparent resistivity and magnetic susceptibihty calculation: Expanded Abstract of the 67th SEG annual meeting, pl282-1285. [58] Raiche, A.P., Jupp, D.L.B., Rutter, H., and Vozoff, K., 1985, The joint use of coincident loop transient electromagnetic and Schlumberger sounding to resolve layered structures: Geophysics, 50, 1618-1627. [59] Ryu, J., Morrison, F. H., and Ward, S. H., 1970, Electromagnetic fields about a loop source of current: Geophysics, 35, 862-896. References 189 [60] Ryu, J . , Morrison, F. H . , and Ward, S. H . , 1972, Electromagnetic sounding experiment across Santa Clara Valley: Geophysics, 37, 351-374. [61] Sengpiel, K . P., 1988, Approximate inversion of airborne E M data from a multilayered ground: Geophys. Prosp., 36, 446-459. [62] Sinha, A . K . , 1968, Electromagnetic field of an oscillating magnetic dipole over an anisotropic earth: Geophysics, 33, 346-353. [63] Sinha, A . K . , 1969, Vertical electric dipole over an inhomogeneous and anisotropic earth: Geoexploration, 7, 9-28. [64] Sinha, A . K . , 1973, Comparison of E M coil systems placed over a multilayered conducting earth: Geophysics, 38, 894-919. [65] Sinha, A . K . and Bhattacharyya, P.K., 1967, Electric dipole over an anisotropic and inhomogeneous earth: Geophysics, 32, 652-667. [66] Slichter, L . B . , 1933, A n inverse boundary value problem in electrodynamics: Physics, 4, 411-418. [67] Slichter, L . B . and Knopoff, L . , 1959, Field of an oscillating magnetic dipole on the surface of a layered earth: Geophysics, 24, 77-88. [68] Smith, J.T. and Booker, J.R., 1991, Rapid inversion of two- and three-dimensional magnetotelluric data, J . Geophys. Res., 96, 3905-3922. [69] Speies, B . R. and Frischknecht, F. C , 1988, Electromagnetic sounding: in Nabighian, M . N . , E d . , Electromagnetic Methods in Applied Geophysics: Investigations in geophysics, Vol. 2, Soc. Expl. Geophys. publications, 289. .References 190 [70] Torres-Verdin, C and Habashy, T . M . , 1994, Rapid 2.5-dimensional forward modeling and inversion via a new nonlinear scattering approximation: Radio Sci., 29, 4, 1051-1079. [71] Tripp A . C , Ward, S.H., Sill, W.R., Swift, C M . , Jr., and Petric, W . R . , 1978, Electromagnetic and Schlumberger resistivity sounding in the Roosevelt Hot Springs K G R A : Geophysics, 43, 1450-1469. [72] Vanyan, L . L . , 1967, Electromagnetic depth soundings selected and translated by G.V. Keller, Consultant Bureau, New York, N Y , 312p. [73] Verma, R . K . , 1973, Topographic effects on electromagnetic depth sounding systems: Geophys. Prospect., 21, 1-25. [74] Wang, X . and Hansen, R. 0., 1990, Inversion for magnetic anomalies of arbitrary three-dimensional bodies, Geophysics, 55, 1321-1326. [75] Wait, J.R., 1951, A conducting sphere in a time varying magnetic field: Geophysics, 16, 666-672. -1952, Current carrying loop in a single inhomogeneous region. J . Appl. Phys., 23, 497. -1953, Induction by a horizontal oscillating magnetic dipole over a conducting homogeneous earth. Trans. Am. Geophys. Union, 34, 185-188. -1955, Mutual electromagnetic coupling of loops over a homogeneous ground: Geophysics, 20, 630-637. -1956, Mutual electromagnetic coupling of loops over a homogeneous ground - an additional note: Geophysics, 21, 479-484. References 191 -1958, Induction by an oscillating dipole over a two-layer ground: Appl. Sci. Res., B-7, 73-80. -1962, A note on the electromagnetic response of a stratified earth: Geophysics, 27, 382-385. [76] Ward, S.H., 1959, Unique determination of conductivity, susceptibihty, size, and depth in multi-frequency electromagnetic exploration, Geophysics, 24, 531-546. [77] Ward, S. H . , and Harvey, H . A . , 1954, Electromagnetic surveying of diamond drill holes: Canadian Mining Manual, National Business Publications, 19-30. [78] Ward, S. H . and Hohmann, G. W . , 1988, Electromagnetic theory for geophysical applications: Electromagnetic methods in applied geophysics, S E G publication, 131-311. [79] Ward, S.H., Ryu, J . , Glenn, W . E . , Hohmann, G . W . and Smith, B . D . , 1974, Electromagnetic methods in conductive terranes: Geoexploration, 12, 121-183. [80] Ward, S.H., Smith, B . D . , Glenn, W . E . , Rijo, L., and Inman, J.R., 1976, Statistical evaluation of electrical sounding methods: Geophysics (Part II), 41, 1222-1235. [81] Whittall, K . P., and Oldenburg, D. W . , 1992, Inversion of magnetotelluric data for a one-dimensional conductivity: Geophysics Monograph Series, 5, S E G publication. [82] Vidberg, H.J. and Riska, D.O., 1985, The inverse scattering problem for reflection of electromagnetic dipole radiation from Earth with vertical variation: J . Geophys., 56, 183-191. [83] Zeyen, H . and Pous, J . , 1991, A new 3-D inversion algorithm for magnetic total field anomalies, Geophys. J . Int. 104, 583-591. References 192 [84] Zhang, Z. and Oldenburg, D. W . , 1996a, Recovering magnetic susceptibihty from electromagnetic data over a 1-D earth: accepted for publication by Geophys. J . Internat.. -1996b, Simultaneous reconstruction of 1-D susceptibihty and conductivity from E M data: Expanded Abstract of 67th S E G annual meeting, 241-244. [85] Zhang, Z., Routh, P., and Oldenburg, D.W., 1996, Reconstruction of 1-D conductivity from coaxial, coplanar, vertical coplanar and perpendicular E M data: Expanded Abstract of 67th S E G annual meeting, 1294-1297. [86] -1994, Inversion of airborne data over a 1-D earth: Annual report for 1993, JACI consortium meeting, appendix H . [87] Zhdanov, M.S., and Fang, S., 1995, 3-D electromagnetic inversion based on quasilinear approximation: The abstracts for the international symposium on 3-D electromagnetics at Schlumberger-Doll Research center, 265-272. - 1994, Quasi-linear approximation in 3-D electromagnetic modeling: submitted to Geophysics. [88] Zollinger, R., Morrison, F . H . , Lazenby, P. G . , and Becker, A . , 1987, Airborne electromagnetic bathymetry: Geophysics, 52, 1127-1137. Appendix A Zero frequency solution for Hz over a half-space The asymptotic expression for the vertical component of the magnetic field , due to a dipole source of moment m, over a general 1-D earth at zero frequency is given by UmH (r^z ) z = £ oba jT | _ ^ A > ^ - » f c e )j (Ar)dA. (A.l) 0 For a half space of fi = fii and a = a\, the input impedance equals the intrinsic impedance Z 1 = Z (A.2) L y/\ — U> €ofll - f ' i w / i i O ! 2 2 - When the frequency tends to zero, the normalized input impedance is A. (A.3) Hi. -ILO Thus expression (A.l) can be further reduced to _ m (m - /xp) l i m f l . C r . w , * ^ ) = ^-fL^El f*~ A V ^ - > J ( A r ) d A . w—>0 47T (/i! + /i„) Jo 2 h 0 (A.4) From Gradshteyn and Ryzhik (1965, p712 ), poo J m+l (Va x e- J (bx)dx ax v = (-l) b~ da m+1 v IT r \ obs m „_<, ' ° w 4 f f v Va 2 + b 2 (A.5) = 0 yields ; km HAr.u), z .) = — ^ ob a) 2 v m+1 Letting m = 1, b = r,v = 0, a = 2h and z + b- 2 d m+1 (^ 193 _ 1 + / i 0 /*o) 8h - r ^—r-, . )(4/ 2 2)2.5' 2 l 2 + r (A.6) Appendix B Adjoint Green's function solution for the sensitivities From the first two equations in equation (4.5) I obtain d 8E 2 dz . d dH . r dpi 2 dp dH . r dz dpi dz dH , r . T dtbi T dz dpi (B.l) dz and -{dr 9H d dH Ir d V dpij Z dr dpi " z dr ' (B.2) dr Note that second terms on the right hand side of both equation (B.l) and (B.2) are zero because permeability in each layer remain constant and ipi is only a function of z. Multiplying both sides of the third equation in (4.5) with . t W d dH . r t i ¥ z T p - d dH . z ~ ^-dr-dpl = Solving equations (B.l) and (B.2) for iup-f^ . I obtain IUJ/J,, . a ^ dE x (B.3) ^ d p - - + and iup£^, and inserting these into equation (B.3) leads to d 2 dz 2 dE (. „di>i dpi + —I- dn; iu>H — dz d_( 1- tu)ipi dpi r dr = Ir dr dET \ dp + r i j dE (—iuj jj,ie + iuia) 2 0 dpi' tu-^—ifi or (B.4) Above equation can be reorganized into d 2 dE d dz dpi dr ( 1 d ( dE s dE (—ia} pe + iua) 2 2 dr \ 0 dp it 194 dpi (B.5) Appendix B. Adjoint Green's function solution for the sensitivities . , fdH dH \ r . z 195 dibi , N Since there is no artificial source in the ith layer, dH dH r z dz dr = (iioe + cr)E. (B.7) 0 Hence the sensitivity problem can be expressed as d dE 2 dz dm 2 [ 1 f9r I r d d ( dE\]\ . dE . 2 = iwi>i(iwe + cx)E + iwH ^. 0 (B.8) r Since the box car basis function is actually the subtraction of two Heaviside functions, ipi = H(z — Zi) — H(z — z ), the derivative of the box car with respect to the depth is i+i din dz = 8{z - z^ - S(z - z ). (B.9) i+1 Due to its symmetry this problem can be converted into the Hankel domain (Appendix B) d 2 dE dE 2 . , N r , i dE . rc , c . _ i n , where u = A — o> /zen + iwfio~. Compared with the partial differential equation for 2 2 2 the sensitivity of conductivity, this sensitivity problem has two extra terms located at the upper and lower boundaries of each layer. I solve equation (B.10) by introducing a Green's function G. Multiply both sides of equation (B.10) by G and integrate over the whole domain to obtain equation (4.15). The boundary term on the left hand side vanishes because the electric field for any finite source tends to zero at infinity and so do its derivatives. Thus if the Green's function satisfies equation (4.11) , then the sensitivity for the electric field is given by equation (4.12). ft Appendix C Hankel transformation for equation (B.8) The Hankel transform is defined as ~ d L f1 d ( dr {r [dr \ dE\ dm}\ °° ( = / Jo rJi(\r)dr d dE d BE or dm dram 2 2 IdE l „ r rdm/ (Cl) Integrate by parts to obtain d 3E + /o rcV [7,Jl(Ar)] + ^ r (C.2) Since the Bessel function of first order equals zero at r = 0, and the electric field and its derivatives diminish at infinity, the boundary terms in the above equation are equal to zero. The integrand can be further simplified by expanding and reorganizing aJx(Ar) Jx(Ar) 5 - 5r r r flJ^Ar) ^ r ) 2 2 J l ( A 7 , j d 2 + 2 _ _ 5 ; (9r . c9Ji(Ar) „ J (Ar) t / % N r , .(C.3) Let Xr - R, such that ^ = ^ f ^ r = A ^ . Carrying out this replacement in cooperation with the definition of the Bessel's function BR 2 RdR R 2 196 RX = -3 (R)RX, l (CA) Appendix C. Hankel transformation for equation (B.8) 197 (C.2) can be reduced to The complete expression of the Hankel transformation for the partial derivative equation of | ^ is d 2 dE 0E xi ^2 2 . ,.. _ 1 8E + <r)E + —K-[S{Z r c f = iuipi{iuje 0 .•" f . - «,-) - £ ( z - z )} i+1 . cfy<9tf + w / - ^ . r (C.6) Appendix D Equivalence of equation (4.12) and (4.15) To prove that equation (4.12) is equivalent to equation (4.15), I only need to show that G dE_ Zi+1 _ 2uVo t = Hi dz G(X,w,z)E(X,v,z)dz - 2Aibihi} (D.l) First of all, the generic solution for the electric field and the Green's function in the ith layer take the following form Ei = Aie * u d = a z + Zi Uiiz ie - Xi Bie-^-^ + bie~ - > Ui{z (D.2) Zi> Therefore 2Aibihi = / 2Aibidz = — / _ 1 2 f 1 Zi+1 2EQ J e Zi where E 0 = (dE uf ic^o/aJa^^e" * 0 0 6 '. i dz = d. 2A Bidz { EQ Jzi Jzi =r [' 1 [EG - 2 JJ zi Zi dGdE dz, u] dz dz J (D.S) Inserting equation (C-3) back into the right hand side of equation (D.3) results in 2u| / G(\,uj,z)E(X,u;,z)dz-2Aibihi J Jzi =1 f +1 2 ^ dGdE\ (u O E + — — 22 <9z & GE =^ pii I, <1=Zz 1 2 « dz c9z 2 J Zi 1 / / d £ 2 Gdz <iz ac?^ Integrating by parts I obtain 1 f /•«.-+! fgd^E Hi J'i G dE dz dz 2 "+1 5Ga/£ JZi 5z <9z dz <9z" 198 dGdE\ dz dz J _] dz L m (D.5) Appendix E Computation of the discrete sensitivities In the calculation of the sensitivities, I need to evaluate an integral / , defined by /(A,w) = P 1 G(X,u,)E(X )dz = tU E\X,u)dz, J Zi (E.l) H/Q Jzi with E defined as in Appendix D. The partial differential equation for the electric field 0 in the Hankel-transformed domain is dE 2 = 0. u 2 dz 2 (E.2) Hence . . . v 1 r*^Ed E, 1 2 ft+i {f„dE\ ti+1 (dE\ 2 } , The general solution of equation (E.2) is a Hnear combination of up-going and down-going waves, as given in equation (D.2). Thus 2 , r i+i (8E\ z /. (&) L d z = r i+i , z (E.4) Inserting this equation back into equation (E.3) yields f(X,u,z) 1 2E u 2 0 J-z E + = 2AiBihi. (E.5) Since at the ith interface dEj _ dz m fii-i 199 dEi_ x d z ' (E.6) v Appendix E. Computation of the discrete sensitivities 200 equation (E.l) can be, for the convenience of computation, further written as Ei 2E u 2 0 dEj 8z Ei-i dE^ + 2AiBihi, (E.7) 2.-+1 where Ei dEj dz 2BjZ i+1 iumBi Z^ + ZiZ* 1 t+i + Zi -2u,hi (E.8) Appendix F The trade-off between u and fi in the sensitivities Introducing a new parameter r(pi, o^) = jajm&i, I can rewrite the perturbation on the datum in equation (5.20) as M sn an M 'dDi dr t dDj,dr_ s dr dcri dr dm * dA. ^ dm \ (F.l) Since dr/dm = dr/dcri x <Ti/m, dr , dDi dr + — —8 dr do~i dr dm dDi i = * [Sen + ^8 • dr dcri \ m 1 (F.2) / Noticing that dDi dr dr do~i (F.3) (Ja)li, I can further reduce equation (F.l) to ££>;(<?-,a;, z M ob4 ) « Y\{J„)n i=i This is the same as equation (5.20), if dDi/dm 201 / a- \ [8^ + —8m j + \ p* J dD is defined as Qu. dm Pi- (F.4)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Reconstruction of conductivity and susceptibility from...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Reconstruction of conductivity and susceptibility from the inversion of EM data Zhang, Zhiyi 1997-12-31
pdf
Page Metadata
Item Metadata
Title | Reconstruction of conductivity and susceptibility from the inversion of EM data |
Creator |
Zhang, Zhiyi |
Date | 1997 |
Date Issued | 2009-04-17T18:19:20Z |
Description | Both electrical conductivity and magnetic susceptibility are important physical parameters in exploration geophysics, and information about their distributions can be used to determine subsurface structures, and to detect mineral deposits and other natural resources. The ultimate goal of this thesis is to recover conductivity and susceptibility from the inversions of electromagnetic (EM) data from various loop-loop systems. A large number of frequency-domain EM (FEM) data are taken in EM surveys by using different loop configurations, and the inversions of these data may provide independent information about the geological targets. Previous work on the inversions of EM data has involved only horizontal loop sources. Consequently, data measured with other coil systems have been constantly rejected from the inversions. In Chapter 2, I investigate the effect of coil configurations on the inversion and develop an inverse algorithm to invert EM data from different coil systems. EM data can also be measured in the boreholes. Large loop systems which measure transient EM (TEM) data on the ground or in the borehole have found increased application in exploration geophysics. However, the inversion of borehole TEM data has not been fully addressed. In Chapter 3, I investigate whether, and how, the use of borehole data in inversions enhances the recovered models. In geophysical explorations, EM responses are a function of the geometry, conductivity, and susceptibility. The influence from magnetic susceptibility on EM data has long been appreciated, but no existing literature has been found about the reconstruction of susceptibility through rigorous inversions. In most cases magnetic susceptibility has been treated as a source of "contamination" in the inversions, and people have been trying very hard to eliminate that "contamination" by truncating or disregarding the inphase EM data in carrying out inversions. In doing so, useful information about the distribution of magnetic susceptibility is wasted, and the recovered conductivity models become less reliable. In Chapter 4, I study the effect of susceptibility on the data, and reconstruct susceptibility from the inversion when the conductivity distribution is specified. The problem with the individual inversions in Chapters 2, 3 and 4 is that accurate information about either conductivity or susceptibility is required in order to recover the other. Thus, it is necessary to explore the possibility of reconstructing 1-D conductivity and susceptibility simultaneously from the inversion of the EM data. In carrying out simultaneous inversions in Chapter 5, I minimize a global model objective function, which includes both conductivity and susceptibility, subject to the data constraints. The final conductivity and susceptibility models are obtained by adjusting the parameter that controls the relative weighting between the two terms in the model objective function. Much of the work in Chapter 2 to Chapter 5 is done in 1-D environment, while geological targets are usually 3-D. Idealy one would like to carry out 3-D inversion to obtain information about the 3-D targets. Currently, however, full 3-D rigorous inversions of EM data are computationally prohibitive, and approximate 3-D inversions are necessary. In Chapter 6, I develop a linear mapping that can be used to in interpret the data collected in a 3-D environment. The algorithm is applied to a field data set collected over Mt. Milligan. All algorithms have been tested on both synthetic and field data sets. Those tests show that the algorithms are robust to different error assignments, even when the error is correlated. The recovered conductivity and susceptibility models from the inversions of field data have provided useful geological information. |
Extent | 23178684 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-04-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052317 |
URI | http://hdl.handle.net/2429/7309 |
Degree |
Doctor of Philosophy - PhD |
Program |
Earth and Ocean Sciences |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- ubc_1997-251950.pdf [ 22.1MB ]
- [if-you-see-this-DO-NOT-CLICK]
- Metadata
- JSON: 1.0052317.json
- JSON-LD: 1.0052317+ld.json
- RDF/XML (Pretty): 1.0052317.xml
- RDF/JSON: 1.0052317+rdf.json
- Turtle: 1.0052317+rdf-turtle.txt
- N-Triples: 1.0052317+rdf-ntriples.txt
- Original Record: 1.0052317 +original-record.json
- Full Text
- 1.0052317.txt
- Citation
- 1.0052317.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
Ukraine | 53 | 0 |
China | 17 | 21 |
United States | 16 | 2 |
France | 11 | 0 |
Russia | 6 | 0 |
Canada | 4 | 0 |
Philippines | 2 | 2 |
Australia | 2 | 0 |
United Kingdom | 2 | 2 |
South Africa | 1 | 2 |
Poland | 1 | 0 |
Germany | 1 | 0 |
New Caledonia | 1 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 73 | 11 |
Ashburn | 10 | 0 |
Saint Petersburg | 6 | 0 |
Wuhan | 5 | 0 |
Shenzhen | 5 | 21 |
Beijing | 4 | 0 |
Montreal | 2 | 0 |
Seattle | 2 | 0 |
Thornhill | 2 | 0 |
Kunming | 2 | 0 |
Fort Worth | 1 | 0 |
Bandung | 1 | 0 |
Julich | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052317/manifest