- 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 1997
pdf
Page Metadata
Item Metadata
Title | Reconstruction of conductivity and susceptibility from the inversion of EM data |
Creator |
Zhang, Zhiyi |
Date Created | 2009-04-17 |
Date Issued | 2009-04-17 |
Date | 1997 |
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 |
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 |
Degree |
Doctor of Philosophy - PhD |
Program |
Earth and Ocean Sciences |
Affiliation |
Science, Faculty of |
Degree Grantor | University of British Columbia |
Graduation Date | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/7309 |
Aggregated Source Repository | DSpace |
Digital Resource Original Record | https://open.library.ubc.ca/collections/831/items/1.0052317/source |
Download
- Media
- ubc_1997-251950.pdf [ 22.1MB ]
- 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
- Citation
- 1.0052317.ris
Full Text
RECONSTRUCTION OF CONDUCTIVITY AND SUSCEPTIBILITY FROM 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 i n 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 V6T 1Z4 Date: Abstract Both electrical conductivity and magnetic susceptibility are important physical para- meters 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 (FEM) data are taken in E M surveys by using different loop configurations, and the inversions of these data may provide inde- pendent 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 investig- ate 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 (TEM) data on the ground or in the borehole have found increased appli- cation 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, conductiv- ity, 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 geolo- gical 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 col- lected 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 1 2 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 Forward modeling for other coil configurations 16 2.3 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 2.4 The inverse algorithm 21 2.4.1 Joint inversion of the coplanar and coaxial data 23 2.5 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 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 4 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 Adjoint Green's function solution 98 4.4 Synthetic examples 104 4.5 Field example 109 4.6 Conclusions 114 5 1-D Simultaneous Inversion of E M data 117 5.1 Introduction 117 5.2 Methodology 120 5.2.1 The trade-off between conductivity and susceptibility 122 5.3 Numerical results 127 v 5.4 Field data example 139 5.5 Summary 148 6 Approximate inversion of 3-D E M data 150 6.1 Introduction 150 6.2 The methodology 153 6.3 Calculation of the sensitivities ! 156 6.4 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 7 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 Computation of the discrete sensitivities 199 vi F The trade-off between a and /x in the sensitivities 201 vi i List of Tables 5.1 D(a, fi0), D(a, / x 0 ) , A D ( c r , /x), and D(<x, fi) over a 0.01 S / m a n d 0.1 SI un i t half-space 5.2 T h e c o m p a r i s o n of the sensi t ivi t ies for the c o n d u c t i v i t y and suscep t ib i l i ty , f r o m a layer of 2 m th i ck at a d e p t h of 4 8 m , i n a 0.01 S / m half-space whose suscep t ib i l i ty is 0.1 SI un i t v m 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 25 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 27 2.6 The influence on the data from the resistive and the second conductive zones in the model used in the first example 28 2.7 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 38 2.14 The real component of the perpendicular data in ppm 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 hori- zontal 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 Mt . Milligan. 55 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 64 3.2 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 10 _ 1 seconds is used as the input data 66 3.4 The geometry for the inversion of borehole E M data 68 3.5 Frequency and transient responses over a half-space of 0.01 S/m for a central loop borehole system 73 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 74 3.7 The sensitivities for both surface and borehole central loop configurations in the time- and frequency domains 76 3.8 The sensitivities for borehole central loop configurations in the time- and frequency-domains in a half-space of 0.01 S/m 78 3.9 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 86 xi 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 inver- sion 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 in- version 106 4.6 D I G H E M data from Mt. MiUigan, at section Y9600 108 4.7 Inversion of D I G H E M data from Mt. Milligan, section Y9600 I l l 4.8 Recovered susceptibility models from inversions with different error schemes. 113 5.1 Individual inversions with accurate and inaccurate information about con- ductivity or susceptibility 117 5.2 The relationship among D(cr,fj,), D(a,fi0), D(<7,//0), and AD(<r,fi). . . . 123 5.3 The result of the simultaneous inversion with s = 6 126 5.4 The results from the inversion with s — 20 127 5.5 The inversion with s = 0.1 129 5.6 The 3-D model 130 xii 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 132 5.8 The Imaginary component of the predicated and observed data at 110, 7040 and 56320 Hertz for the 3-D model 133 5.9 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 AEM data 142 5.15 The x-y planview of the recovered susceptibility models from the simulta- neous inversion of AEM data and the inversion of aeromagnetic data. . . 144 5.16 The cross-section of the recovered susceptibility models from the simulta- neous inversion of AEM 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 160 xm 6.5 The real [panel (a)] and imaginary [panel (b)] components of the sensitiv- ities 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 Mt . 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 Mt . Milligan 169 6.11 The X - Y plan-view of the recovered conductivity from the 3-D approxim- ate inversion of D I G H E M data at Mt. 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 Mt . Milligan 171 x i v Acknowledgement I would like to thank my supervisor Prof. Doug Oldenburg for his advice and encourage- ment 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 x v i 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 (EM) 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 in- terpretation 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 advant- ages 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. How- ever, 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 (FEM) 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 (TEM) 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 meas- ure 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 meas- urements. For detection of very good conductors, B measurements will often give su- perior 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 Chapter 1. Introduction 4 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 con- figurations. Four commonly used loop-loop configurations are the horizontal coplanar (HC), vertical coplanar (VC), perpendicular (PP), and vertical coaxial (CA). In the PP 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 ter- minology. Theoretical results are often displayed using induction number. Typically the induction number B = {a\iLcj1yl2T is used in the frequency-domain, and the induction 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, depend- ing 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 two- layer 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), Ryu 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. Ryu 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 Ryu 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 pro- grams. 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 sup- posed to represent that portion of the earth under consideration. Normally the model- fitting 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, com- paring 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 appli- cations 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 non- uniqueness 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 Chapter 1. Introduction 7 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 para- meters m , and the matrix J is the Jacobian matrix of the data whose elements are given by j ( m ) = i = i , 2 , M , j = 1, 2 , N , (1.2) Orrij 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 inver- sion 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 hori- zontal 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 orienta- tion 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 sus- ceptibility 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 HC system with conjugate gradient method, quasi- Newton method, and Gaussian-Newton method. A l l techniques were very time consum- ing, and because of the limitation on computer memory, the model was severely under- parameterized. 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 simultane- ous 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 con- straint 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 pre- vious chapters. By 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 be- tween 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 approxim- ation 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 dif- ferent 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 per- centage 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 inver- sions 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 other p h y s i c a l const ra ints , to de te rmine the appropr i a t e er ror assign- men t . A s demons t r a t ed i n n u m e r i c a l examples , m y a l g o r i t h m is robus t w i t h different error ass ignments , a n d can ex t rac t i n f o r m a t i o n about the targets even w h e n the error i n the d a t a is cor re la ted . 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 orienta- tions 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 (VMD) 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 13 HC CA Tx Rx Tx R* A PP VC Tx A Tx Rx Figure 2.1: Four commonly used coil configurations: the horizontal coplanar (HC), the perpendicular (PP), the coaxial (CA), and the vertical coplanar (VC) systems. Tx 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. My 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 14 Tx a Z ( m ) j i | 2 o b s * ""i Rx R(m) earth's surface CT 1 K 2 Ki °1 Kj+1 ai+i K M a M K M + 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{. Z0bs is the vertical distance between Tx and Rx. 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 Iexut, where u> is the angular 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 = (2.1) oz 1 d iuHHz = — — (rEe), (2.2) r or and dHr dHz (iuje0 + o-j)Ee + J3, (2.3) dz dr where Js = I(uj)S(r - a)S(zobs), and i = z0bs is the distance between the source plane and the receiver plane. By ehminating the magnetic field Hr and Hz from equations (2.1), (2.2) and (2.3), I obtain an inhomogeneous scalar equation / d2 d2 Id 1 \ + + " a 9 + kj) Ee{r,z,u>) = iujfi0I(uj)6(r - a)S(zobs), (2.4) ydz2 ' dr2 r dr r2 3 , where k2- = u>2fij€0 — iu;[ij(Tj. Equation (2.4) can be converted into Hankel transform domain. The partial differential equation for the electric field, in the Hankel transform domain, is d 2 dz2 ^ 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 u2- = A 2 — k2. After inverse Hankel transformation, the secondary electrical field is given by r°° A 2 7^ Z E(r,zohs^-) = -iu>^aI(u) ^ _ - ^ _ ^ e « o ( * o . . - 2 f c ) j 1 ( A a ) J 1 ( A r ) A i A , (2.6) where the input impedance of the j th layer is given by (Wait, 1962) • _ Z ^ + -̂tanhK-fe,-) J Z j + ^ ' + 1 t a n h ( ^ j ) ' 1 ' Chapter 2. Recovering 1-D conductivity from the inversion of EM data 1 6 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 imped- ance is equal to the intrinsic impedance. That is, ZM = ZM- After inverse Hankel transformation, I obtain the secondary magnetic field in the frequency-domain zw/io Jo The induced secondary voltage measured in a receiver coil is the time derivative of the 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 PP 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 ( 2 . 9 ) magnetic flux and is expressed in the frequency domain as ( 2 . 1 0 ) Chapter 2. Recovering 1-D conductivity from the inversion of EM data 17 V C , C A , PP and HC systems are denoted as Hvc, HCA, HPP and Hue, then under the quasi-static assumption, the forward responses are Hca(T,LO,ZoU) = — - / RTEe-^2h-^XJ1(Xr)d\ 47rr Jo + ^ RTEe-X{2h-^h2J0(\r)d\, (2.11) 47T JO Hvc(r,u;,zoba) = — ^ e ^ - ^ ^ A J ^ A r ) ^ , (2.12) 47TT Jo HPP(r,u;,zobs) = — - / fee-^-^^A^^Ar)JA, (2.13) 772 /*oo HHcir^z*.) = — RTEe-^2h-z^X2J0{\r)dX, (2.14) 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, Z1, can be found by the recursive procedure outlined by Ryu et al. (1970). I note that the coil configurations do not provide independent information. In fact HCA = HHC — Hvc- The 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 18 Frequency ' Frequency F i g u r e 2.3: T h e secondary magne t i c f ield over a half-space of 0.01 S / m . Pane l s (a), (b) , (c) , a n d (d) show the absolute value of the secondary magne t i c f ie ld f r o m the coax ia l , pe rpend icu la r , v e r t i c a l coplanar and ho r i zon ta l coplanar co i l systems. S o l i d l ines denote the rea l components , and dashed lines denote the i m a g i n a r y components of the da ta . components of the d a t a f rom the other three configurat ions are larger t h a n the ampl i tudes of the real components at l ow frequencies. 2.3 Calculation of the sensitivities I n order to ca r ry out a r igorous invers ion , one needs to compu te the sensi t ivi t ies J , w h i c h connect the per tu rba t ions on the m o d e l m to the per tu rba t ions on the d a t a D(m). T h e elements of J are Chapter 2. Recovering 1-D conductivity from the inversion of EM data 19 ddi(m) (̂111) = - ^ ^ , i = l , 2 , . . . , M , j = l , 2 , . . . , iV , (2.16) orrij 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 d2 dE 2dE . Uj—= lujfijipjE. (2.18) dz2 do-j ~3 do-j The sensitivity problem can be solved by using adjoint Green's function solution (Zhang and Oldenburg, 1994): d E ^ Z ^ = i„N T + 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-U2:)Gj(\,Z,UJ) = 8 ( z - Z o b s ) Gj(\,zj,u,) = Gj+1(\,zj1u>) (2-20) G(\,z,u>) —> 0 when \z\ —> oo, Chapter 2. Recovering 1-D conductivity from the inversion of EM data 20 The sensitivities for the magnetic field are related to the sensitivities for the electric field through Faraday's law dH*f>z>") = - J - f dE^>Z>Uh0(\r)\>d\. (2.21) O0~j tUflj JO UO~j 2.3.2 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. By 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 em- bodied 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 dHXY{\oj,z) _ foo dR r ^ T E 0 X Y { X O J Z ) D K ( 2 .22) JO OCj do-j 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 _M2h_z)dE(\,u,z) ^ do-j 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 2 J 0 ( A r ) e - A ( 2 ^ 2 - ) 47T Chapter 2. Recovering 1-D conductivity from the inversion of EM data 21 Ovc = — A J 1 ( A r ) e - A ( 2 ' 1 - ^ ) OPP = - ^ A 2 J 1 ( A r ) e ~ A ( 2 ' 1 - ^ ) (2.24) 4TT OCA = OHC — Ovc One is listed here for completeness. 2.4 The inverse algorithm In nonlinear inverse problems, one is provided with observations d°bs,i = 1,N, and an associated error estimate £j for each datum. 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>m = a J wa(z)(m - m0)2dz + (1 - a) Jwf(z) penalizes vertical roughness and differences between the recovered model and a reference model m0. In equation (2.25) a is a parameter that controls the relative importance of the two terms. wa and Wf are weighting functions which can be prescribed by the user. When the earth is divided into layers of constant conductivity, equation (2.25) can be d(m — m0) dz. (2.25) Chapter 2. Recovering 1-D conductivity from the inversion of EM data 22 discretized and written as K =|| Wm{m-m0) | | 2 , (2.26) where m = (mi,m2, ...mj^) is a model parameter vector and Wm is an M x M weighting matrix. I chose a data misfit objective function fa =|| Wd(D°bs - D) || 2= £ ( - ? - ^ i ) 2 , (2.27) i=l Vi where rji is the error associated with the i th datum, and Doha and D are the observed and predicted data respectively. My goal is to find a model m that minimizes equation (2.26) subject to the constraint that (f)d in equation (2.27) is equal to a target misfit (j>*d. If the errors are Gaussian and independent then (f)d is a chi-squared variable and 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>d = <j>*d requires minimizing ^(m) = ^ m + i 9 - 1 ( ^ - ^ ) , (2.28) where 0 is a Lagrangian multiplier. The nonlinear nature of the problem requires lin- earization and iteration to a solution. Let m(n) be the model at the nth iteration and let 8m be a perturbation. The effect of the perturbation on the z'th datum is given through a Taylor's expansion M o 771 M di[m^ + 8m] ~ Fi[m^} + £ -^8mi = d\n) + E JH6mh (2-29) j = l omj 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 th layer. The problem becomes: minimize 4> =|| Wm[8m + rn'"' - m 0] | | 2 + /T 1 { | | Wd{D°hs - F[m^ + 8m}} | | 2 -fd(n+1)}, (2.30) where <^ n + 1 ^ is a target misfit at each iteration. Generally I choose (j)*}n+1^ = l/*f(j)*dn^ where 7 > 1. Writing F[m^n+1^] = F[m^] + J8m and setting the gradient Vsm4> = 0, I Chapter 2. Recovering 1-D conductivity from the inversion of EM data 23 obtain Sm = [f3W^Wm + fWjWaJ^iJTwJWrfDn + f3W^Wm[m0 - m W ] } , (2.31) where SDn is SDn = Dobs - F [ m W ] , (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) II2 + || WdCP{Dtp - DCP) ||2, (2.33) Chapter 2. Recovering 1-D conductivity from the inversion of EM data 24 where WdCA and WdCP a r e 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>m + /3 1(<j>dj ~ 4>*d)- Setting the gradient V$m(p(m) = 0 yields where J = 8D = and Wd dCA \ 0 Wdcp ) (2.34) 8m = [BWlWm + JTWjWdJ]-x{JTWjWd8Dn + (3W^Wm[m0 - m^]}, (2.35) (2.36) (2.37) (2.38) 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 gen- erated 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 in- versions 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 influence on the da t a f rom the resist ive and the second conduc t ive zones i n the m o d e l used i n the first example . T h e real and i m a g i n a r y components of the d a t a used i n the first example are denoted w i t h sol id and dashed l ines , and the rea l and i m a g i n a r y componen t s of the d a t a generated f r o m the same 1-D m o d e l , bu t w i t h o u t the resis t ive a n d the second conduc t ive zones, are denoted w i t h t r iangles a n d circles . Panels (a), (b ) , (c) , a n d (d) present the d a t a f rom the ve r t i ca l coplanar , coax ia l , pe rpend icu la r , and coplanar systems. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 31 I—i i 1 1 H I M — i i i nun—i i 1 1 it I i i 11 i i i M 10° 1 0 1 1 0 2 10° 1 0 1 1 0 2 Depth (m) Depth (m) 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. Recovering 1-D conductivity from the inversion of EM data 32 2.5.2 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 perpen- dicular 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 configura- tions 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 inver- sions 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) 700 500 ' 300 100 | 0 GO Ko 100 1—r" 300 X(m) X-Y Plan-view 500 700 100 0 15 30 ? N 50 100 X(m) 300 500 700 o n Ko X-Z Cross Section Figure 2.8: The 3-D model used to generate the synthetic data. The background con- ductivity cr0 = 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. Chapter 2. Recovering 1-D conductivity from the inversion of EM data 35 Predicted data Observed data o C o 600 400 200 f=880 — i 1 1 1 1 — 200 400 600 600 400 200 600 i 400 200 f=56k 0 200 400 600 Easting (m) 600 400 H 200 / / IB I , f=880 — I 1 1 — 200 400 600 600 400 200 i f f=7040 0 200 400 600 _ 1 1 1 1 1— 600 400 H 200 f=56k - 1 1 1 1 1 1— 0 200 400 600 Easting (m) F i g u r e 2.9: T h e rea l componen t of the v e r t i c a l coplanar d a t a i n p p m . Pane ls (a) , (b) , a n d (c) show the p red ic t ed da ta , and panels (d) , (e), and (f) show the observed da ta . 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 component of the p red ic t ed a n d observed v e r t i c a l cop lanar d a t a i n p p m . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 37 JZ -<—• a m • 50 100 150 i Conductivity (S/m) a i 100 200 300 400 5CC 600 700 E J= Q . CD Q 50 100 150 • >"*-~^__ \ Conductivity (S/m) | 100 200 300 400 500 600 700 EE a o • 50 H 100 150 Conductivity (S/m) 0 100 200 300 400 500 Easting (m) 600 700 F i g u r e 2.11: T h e x -z cross-sections of the recons t ruc ted 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 , f rom the invers ion of the v e r t i c a l coplanar da ta . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 38 S 50 i 100 150 Conductivity (S/m) 100 200 300 400 500 600 700 50 100 H 150 Conductivity (S/m) 100 200 300 400 500 600 700 ! - - - . - y) 50 - JZ 100 - CL 0 150 - Conductivity (S/m) [cj i r I I 100 200 300 400 500 Northing (m) 600 700 F i g u r e 2.12: T h e recovered m o d e l f rom the invers ion of the v e r t i c a l coplanar da t a , at (a) x = 2 5 0 m , (b) x = 3 5 0 m , and (c) z = 3 5 0 m . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 39 600 400 200 o d ( v • ! 1 i z=1 5m a 600 400 200 H z=40m 0 200 400 600 I I l 600 400 H 200 z=20j — i 1 1— 200 400 600 0 200 400 600 I ! I 600 400 200 z=50 200 400 600 600 - 600 - 400 - fi 400 - 200 - 1 1 W 200 - 1 : z=30m j z=60m if] o J • i i i 0 J 1 1 ! 0.195 0.139 0.100 0.071 0.051 0.036 0.026 0.019 0.013 0.010 0.195 0.139 0.100 0.071 0.051 0.036 0.026 0.019 0.013 0.010 0.195 0.139 0.100 0.071 0.051 0.036 0.026 0.019 0.013 0.010 0 200 400 600 Easting (m) 200 400 600 Easting (m) F i g u r e 2.13: T h e x - y p lan-v iews of the recovered conduc t iv i t y , f rom the inve r s ion of the v e r t i c a l cop lanar da ta , at s ix different depths . 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 inform- ation 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 over- fit 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 struc- tures 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 perpen- dicular 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 41 Predicted data Observed data o 600 g 400 o 200 600 400 200 0 200 400 600 200 400 600 Easting (m) 600 - 400 200 - 0 200 400 600 600 400 200 | f=7040 • el 0 200 400 600 600 400 200 0 200 400 600 Easting (m) F i g u r e 2.14: T h e rea l componen t of the pe rpend icu la r d a t a i n p p m . Pane ls (a) , (b ) , a n d (c) show the p red i c t ed da ta , and panels (d) , (e), and (f) show the observed da ta . Chapter 2. Recovering 1-D conductivity from the inversion of EM data 42 50 100 150 CA 1 i 1 1 : r~ 0 100 200 300 400 500 600 700 Easting (m) 0 100 200 300 400 500 600 Easting (m) 50 pp 0 100 200 300 400 500 600 Easting (m) 700 F i g u r e 2.15: T h e compar i son of recovered models f rom the invers ions of the coplanar , pe rpend icu l a r , a n d c o a x i a l d a t a at sect ion 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. At 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 overbur- den 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 rea l and i m a g i n a r y components of R over the sect ion at y = 2 5 0 m , at 880, 7040, and 56320 her tz . D a s h e d lines denote the i m a g i n a r y componen t , a n d the so l id l ines denote the rea l component of R. T h e grey bar denotes the h o r i z o n t a l p o s i t i o n of the conduc t ive 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 rea l ( the sol id l ine) and i m a g i n a r y (the dashed l ine) componen t s of R at y = 3 5 0 m . T h e ho r i z on t a l posi t ions of the two conduc t ive p r i sms are m a r k e d w i t h the 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 sig- nal 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. Recovering 1-D conductivity from the inversion of EM data 48 F i g u r e 2.19: T h e recovered m o d e l ( y = 4 0 0 m ) f rom the invers ion of the c o a x i a l da ta . T h e so l id and dashed lines p lo t the rea l and i m a g i n a r y components of the p r e d i c t e d da ta . T h e observed d a t a are p lo t t ed w i t h discrete dots. T h e two pr i sms a n d the overburdens are denoted w i t h w h i t e boxes. H s a n d H p are the secondary and the 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 over- burden 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 50 E Q. tx LO x 7500 6000 4500 H 3000 1500 H 56320 Hz -i 1 1 1 3500 m 3000 E 2500 Q. C L 2000 X 7040 Hz] j JL _ * 1 1 /*•*" • a_ _ -̂ — ^ — —» - — • - - ... -r-S b 1400 E 1200 C L C L 1000 C L 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 recovered m o d e l , at y = 4 0 0 m , f rom the inver s ion of the h o r i z o n t a l cop lana r da ta . T h e sol id and dashed lines p lo t the rea l and i m a g i n a r y componen t s of the p red i c t ed da ta . T h e observed d a t a are p lo t t ed w i t h discrete dots . T h e two pr i sms and the overburdens are denoted w i t h wh i t e boxes. H s and H p are the secondary 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 2C0 S/m VC 100 200 200 400 500 600 700 Y(m) 50 1 100 150 200 S/m HC o 100 200 300 400 Y(m) 500 600 700 F i g u r e 2.21: T h e recovered models a n d the da t a f rom the invers ions of the v e r t i c a l cop lanar and h o r i z o n t a l coplanar 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. How- ever, 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 HC 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 -\ Q . -9= 1200 C L X CO X E GL CL CL X «5 X C L Q. C L X c/j X C L CD Q 800 400 - 750 GOO 450 300 150 56320 Hz 7040 Hz I _ _ J> m 300 240 i80 120 H 60 880 Hz 50 100 150 200 S/m Id! i oo 200 300 400 X(m) 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 F i g u r e 2.22: T h e recovered m o d e l f rom the jo in t invers ion of the coarxial a n d the hor i - z o n t a l coplanar d a t a at y = 4 0 0 m . T h e sol id and dashed l ines plot the rea l a n d i m a g i n a r y componen t s of the p red i c t ed coax ia l da ta . 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 pr i sms and the overburdens are denoted w i t h w h i t e boxes. H s a n d H p are the secondary a n d the 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 observed and the p red i c t ed coplanar d a t a f r o m the j o i n t i n v e r s i o n of the c o a x i a l a n d coplanar da ta . T h e so l id a n d dashed l ines denote the rea l a n d i m a g i n a r y componen t s of the p red i c t ed da ta . 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 Mt . 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 DC 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 HC 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 DC resis- tivity 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 con- ductivity model from the inversion of 2-D DC 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 HC, C A , PP and V C systems. A Gauss-Newton algorithm was developed which minimises a norm of the model subject to the data con- straints. The model norm can be constructed to allow the inversion to find the most Chapter 2. Recovering 1-D conductivity from the inversion of EM data 57 12.4 12.6 12.8 13.0 13.2 13.4 50 •g 100 - ^ 150 § 200 - Q 250 12.4 12.6 12.8 13.0 Easting (km) 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 Figure 2.24: The recovered conductivity models over section Y9600 at Mt. 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 DC 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 ade- quately. 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, coax- ial, 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-to- noise 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 lm) 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 bod- ies 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 dif- ference 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 quad- rature 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 (TEM) form or in frequency domain (FEM) version are the down- hole 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 min- eralization 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. By 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 mat- ter 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 Mt. Morgan palaeochan- nel aquifer, which is located about 40 km southwest of Kambalda in the Goldfields region of Australia, was carried out by Western Mining Corporation Limited (WMC) in 1991. This aquifer is being developed as a source of saline process water for the St. Lves Gold Mi l l . 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 sur- face 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 frequency- domain. 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. An 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 IROTEM 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 64 3.2.1 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 r T „ , f f , ^ . , . ; , _ j dt 2 t°° 2 r°° — / Im[H(u>)} sinwtdu! = / Re[H(u>)] cosu/tdto, (3-1) 7T Jo TV Jo or = / —1 v n zosLotdu = h(0) - - / — l—^ J 1smut(ku, (3.2) 7T ./0 LO TT Jo LO where H(u>) and h(t) are the magnetic fields in the frequency- and time-domains respect- ively. 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 r°° dh(t) r°° Re[H(u>)} = — — — cos iotdt = -to [h{t) - h(0)] sin wtdt, (3.3) 27T Jo dt Jo and Im[H(w)] = f°° sin utdt = -to f°° h(t) cos totdt. (3.4) Jo dt Jo Equations (3.3) and (3.4) can be cast into a standard form: Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 65 poo F ( c ) = / f(t) Jo sin wt dt, (3.5) cosa>r. which can be easily calculated by using Anderson's digital filter. One of the major diffi- culties 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 frequency- domains over a half-space of 0.01 S/m. The current excitation in the frequency-domain is elwt 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~3 to 106 Hertz. Figure 3.1 shows the conversion. That 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 time- domain to the frequency-domain. A time span from 10~6 to 10 _ 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 108 Hertz. At frequencies lower 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~6 to 10 _ 1 seconds can provide satisfactory frequency data. For most T E M systems time spans are around the range of 1 0 - 4 to 10 _ 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~4 to 0.161 seconds) was used in the conversion. The converted Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 66 1 Q - 8 I I HII I I I I I I 111 I I I I I MM I ' I I I I i n | i i i i i n i 1 0 - 6 10~ s 10-" 10~3 10" 2 Time(s) Figure 3.1: The conversion of FEM data to TEM 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 - 3 to 106 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. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 67 Figure 3.2: The conversion of step turn-off TEM 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~6 to 10_ 1 seconds. The solid line and the stars denote the true FEM data and the converted FEM data respectively. Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 68 data deviate badly from the true data. 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 10 _ 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~6 seconds must be used in the conversion. In this paper, I calculate the forward responses and sensitivities in the frequency- domain, 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 Ielwt is placed inside an N-layer half-space at a depth of za. The receiver is placed at 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 1 0 - ' 1 0 - 2 10~ 3 10~ 4 1 0 " 5 ^ icr 6 <S 1 0 - 7 i _ io- 8 io- 9 1 0 - 1 0 io- 1' io-' 2 10- 1 3 1 0 - 3 1 0 - 2 1 0 _ 1 1 0 ° 10 1 10 2 10 3 1 0 4 1 0 5 10 8 10 7 1 0 8 F r e q u e n c y Figure 3.3: The conversion of TEM data to FEM data for a central loop system whose loop source has a radius of 50m. The impulse response of the magnetic field which spans 10~10 to 10 - 1 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. 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%ujt is placed at zs in the borehole, and the receiver takes the measurements at a depth of z and a radial distance of r. The governing equation for the electric field is given by Ej(\,z,u>) = iujfijaI(u>)Ji(\a)S(z — zs), ( 3 - 6 ) where A is the Hankel transform parameter, fij is the susceptibility in the jth. layer, and Uj = y A 2 + u>2fij60 — iu>Hj<jj. 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. dz2 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 71 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 com- ponents 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 3 Zj + Z^1 ianh(ujhj) y ' The secondary electric field in region C is given by Ej = A , - e u ' ' ( z _ z ' ' + l ) + Bje-U>(z-Zi+1\ (3.9) where and B ^ B ^ e ^ - ^ ^ + z / ^ (3-10) 7J+1 — Z- Aj = BjeUidi 7 - (3.11) 3 3 Zi+1 + Zj In region A the intrinsic and input impedances are and Zj(X,co) = 1-^, (3.12) z i + i = z ^ + ^ t a n l i ( u A - ) ( 3 . 1 3 ) 3 Z j + Z3 \an\\(ujhj) v ' Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 72 The solution for region A is Ej = Aje^"^ + Bje-u^z-^\ (3.14) where the coefficients are given by and where A*j+1 — Aj+1e~ui+lhi+1, and hj+x - zj+1 - Zj. 3.3.2 Solution in the source-containing layer In the source layer, care has to be taken because of the singularity existing at z = zs. When z is not equal to za, the solution for equation (3.6) is f A j e u z + Bje~uz z<zs Ej(r,z,co) = { (3.17) ( Cjeuz + D j e ' u z z> zs where the coefficients Aj,Bj, Cj and Dj can be determined in the following manner. At z = zs, E is continuous and its first order derivative is discontinuous E j \ z t = E i \ z 7 ' dEj dz + (3.18) '_ = iu>fiaJi(\a). From the definition of the input impedance at the j th and (j + l) th interfaces, I obtain (3.19) Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 73 and C = D e - 2 u ^ + 1 — '-. (3.20) Solving equation (3.18) along with equations (3.19) and (3.20), I obtain A=—e33 —T , (3.21) 2 e-2u>hiR1R2-l ' V ; anc UjZj-ujhj where and C UJZ, I _ Z).= — e"'^- ~ , (3.22) 2 e-2^h>R1R2 - 1 ' V ' Z3' — Z- Ri = 4 - r 4 , (3-23) Z> + Z 3 S is given by ZJ+1 - Z- R2 = fL f i . (3.24) S = ^ a J l ( A a ) . (3.25) 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 za and z are set to zero, the above results are reduced to the response from 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 ex- pressed with exponential functions with negative arguments: x _ -x i _ -2x tanh(z) = — — ^ = — — . (3.26) ex + e x 1 -f- e i 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: Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 74 Hz(r,u>,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. At 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 fila- ment 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~5 s, where 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 75 Frequency Time (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~6s. However, the surface data overpower the borehole data before 3 x 10~6s which indicates that the surface data carry more information about 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 obser- vation 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_2s. The radius of the Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 76 1e-05 0.0001 Time (s) 0.001 1e-05 0.0001 Time (s) 0.001 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 77 3.3.3 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 O0~j J ZJ 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 (& - uf) GJCM,<") = 6(z - Zob,) Gs(X1zj1u) = Gj+1(\1zj,u) (3-29) G(X,z,u>) —» 0 when \z\ —> oo, where z^a is the observation depth in the borehole. Note that G satisfies the same partial 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 ag'r-z-"> = -±-r a ^ ' " ' " ) J . ( A r ) A ^ . (3.30) 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 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. At 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 79 b to \ N X <o i c r 6 1 0 - 7 i c r 8 1 0 " 9 1 0 - 1 0 1 0 " 1 ' 1 0 ~ 1 2 1 0 - 1 3 1 0 - 1 4 1 0 ~ 1 5 _ a \ \ \ \ - / - \ \ \ 1 1 1 \ 1 0 _ 1 1 0 ° 1 0 1 1 0 2 1 0 3 10* 1 0 S 1C Frequency (Hertz) 10 10 10 10 10 eg 10- _ N 10" E 10" 1 0 - - - / / / \ 1 ~\l c 1 1 1 10" 10" 1 0 " 4 1 0 - 3 Time (s) 10" 10" 1 0 - 1 1 0 ° 1 0 1 1 0 2 1 0 3 1 0 4 1 0 S 1 0 6 Frequency (Hertz) 1 0 ° - C/) | 1 0 - 1 v E 1 0 ~ 2 / 1 0 - 3 1 0 " 4 \ 2 \ CD \ 1 0 " 5 \ N - C 1 0 " 6 £NJ CD i o - 7 1 0 " S 1 0 - 6 10~* 1 0 " 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% Gauss- ian 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 flat- test 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 rrSn\ then the line search is terminated, and the last 8 value is adopted. The inversion algorithm begins with an initial conductivity model and at each itera- tion 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 82 1 1 ' 1 1 1 1 I I 1 H H I—I I I I I IH 1 ••' • 0 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 I O - 3 1 0 " ! 1 0 - 1 Depth ( m ) T i m e (s) 0 1 0 0 2 0 0 3 0 0 iOO 5 0 0 I O " 1 1 0 " * Depth ( m ) T i m e (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 Chapter 3. 1-D Inversions of the Surface and Borehole Transient EM data 87 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 4 m 2 . The time span used in this survey is from 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' 100 200 300 400 Depth ( m ) V) 1 10° E < ^—* 10-' IT CO 10-2 CO \ N lO" 3 JZ CN CD 10"* 500 10"3 10~2 T ime (s ) 10-' 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 indic- ates the receiver position. 0 100 200 300 400 500 10"3 10"2 10"' Depth ( m ) T ime (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 89 , , 10' | , Depth ( m ) T ime (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 IROTEM 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 91 • i i i 11 I I i i i i i 1111 i i i i i n i l ' " I i ' ' 1 ' 1 * * i i 1 1 1 1 i c r 3 i o - 2 1 0 - ' 1 0 - 3 10-2 i o - 1 Time (s) Time (s) 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 IROTEM 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 IROTEM 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 92 "ICT3 1CT2 T i m e ( s ) 10 10"3 10~2 T i m e ( s ) io- 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 IROTEM 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 IROTEM 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. 1-D Inversions of the Surface and Borehole Transient EM data 93 3.6 Conclusions The inversion algorithm developed in this Chapter provides a new tool to obtain inform- ation 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. At later time, that relation is reversed, so the surface data have better signal-to- noise 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 suscep- tibihty 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. An 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. Recovering susceptibihty from 1-D Inversion of EM data 96 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 = mb for n < Kb. A n exponential function is used to represent susceptibihty values between KB and K \ . Figure 4.1 shows the mapping. The forward and inverse mappings, m = / ( « ) and K = f~1(m), are given by m ( « ) = K < KB K l [in +1 Kb < K < Ki K K > K i (4.1) and /c(m) m m i < m K\e K i mb < m < mi Kb m < m b (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 mh = K l [ln (^j + l j . (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 Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 98 are measured, the sensitivities for Hz are 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) itouHr(r,e,uj,z) = iupHz(r, 8,u;,z) = -1- {£[rEe(r, 6,u,, z)]} , (4.5) ajiA^A _ SEd^A = { i u j e Q + a ) E e ( r j d , w , z) + Is, where r, 0 and z are variables in the cylindrical coordinate system. Is is the current source given by _ I(u>)a6(r - a)S(zobs - 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*) = £ / M f c ( * ) , ( 4 -6 ) 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 th 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 9 9 given by tlx where the operator C is dm -[S(Z - Zi) - 6{z - z i + 1 ) \ + IU dz dz dfj-i c d*_ dz2 u ( 4 . 7 ) ( 4 . 8 ) 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 G d2 dE dG d dE r<*> dE J-oo dfii dz2dfn dzdzdmW^ ' - o ^ [dz roo G dE / iiotpi(z)(zuie0 + ai)GEdz — J— oo Pi OZ G - u2G dz ( 4 . 9 ) ( 4 . 1 0 ) 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) - u2G(\, to, z) = S(z - zohs) G ( A , w , z ) | , = , r =G{\,u,,z)\t=zt ( 4 - 1 1 ) G(\,w,z) —> 0 when \z\ —> oo, then the sensitivity for the electric field is dE(\,co, zohs) . , . dfii / iu>(iuje0 + o-i)G(X,u, z)E{\,w, z)dz - Jzi Pi dz Z i + l ( 4 . 1 2 ) 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 Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 100 a vertical magnetic dipole with unit strength at the observing point zobs. They are given by Ei(X,u,z) = Ai(X,co)eu^-^ + B i (A,a ; )e- u ' ^ - z ' ) Gi{\,u,z) = ai(\,u)e^z-^ + bi(X,oJ)e-u^-^ where the coefficients of the up-going wave and down-going wave are given by the fol- lowing formulae: Ai(\,u) = B i ( \ , u , ) e - ^ § ^ = < 4 ' 1 4 ) bi(\,w) = . Bifiw\ , ' V ' / v^mlaix (Aa) B0(\,io) - 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,zoba) 2u\ r rzi+i • ^ = / G(X, u>, z)E(\, to, z)dz - 2Aibihi dm m J +iu(iu>e0 +<Ti)Jli+1 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 Zi = Km{~)Zi = m, (4-16) w -»0 no Chapter 4. Recovering susceptibility from 1-D Inversion of EM data 101 be the normalized intrinsic impedance and i ^ h m ( - l ) ^ = / . / i + 1 + / X ' t a n y / l \ ) (4.17) <"-»° iu> pi + pi+1 tanh(A/j.j) be the normalized input impedance. The vertical magnetic field at zero frequency can then be expressed as H^wzob.) = al f°° X\Z! ~ / i ° ) e A ( ^ - 2 f c ) j 1 (Aa) J 0 (Aa)cJA. (4.18) • / - o o 2(Z1 + fi0) If pi = po, where i = 1,2, . . . , M , then Zx — /x0 and hence Hz(ui) will tend to zero. 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 = mp^po 8 f e 2 - r 2 2 47T /Xx +p0(4h2 - r - 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 sep- aration of 10 meters and at height 30 meters above a half-space of 10~2 S/m and a 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 mag- netic 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 103 hertz, the frequency-dependent term begins to dominate. For a given susceptibihty structure the transition frequency decreases as conductivity increases. The imaginary component, on the other hand, is completely frequency-dependent. At low frequencies, magnetic particles change orientation almost synchronously with the primary field, and therefore the amplitude of the imaginary component of the sensitiv- ities 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 con- tributes 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 observa- tion 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 sub- tracting 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. Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 104 frequency-dependent sensitivity is much smaller than the frequency-independent com- ponent; 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 struc- ture is assumed known and the mapping parameters Kb and Ki are fixed at 10~6 and 10~3 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 in- version 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 - 6 SI respectively. The true conductivity model, given in Figure 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 109 4.5 Field example As a field example I now invert airborne E M data acquired at Mt. 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 DC 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. At 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 DC 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 110 i i 1 2 . 4 1 2 . 8 1 3 . 2 X(km) 4 1 2 . 4 1 2 . 8 1 3 . 2 X(km) Figure 4.6: D I G H E M data from Mt. 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. Recovering susceptibihty from 1-D Inversion of EM data 111 beneath each station was used as a 1-D background conductivity for the susceptibihty in- version. 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 KB and K i were set to be 10~6 SI and 10~3 SI respectively. The starting model for inversions at 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 Mt. 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 non- uniqueness 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 Mt. 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 100 H Q- 200 CD Q 1000 (S/m) 12.4 12.6 12.8 13.0 13.2 13.4 12.4 12.6 12.8 13.0 13.2 13.4 8.258 7.395 6.533 5.670 4.808 3.945 3.083 2.220 1.358 0.495 101.10 90.95 80.83 70.71 60.59 50.48 40.36 30.24 20.12 10.00 Figure 4.7: Inversion of D I G H E M data from Mt. Milligan, section Y9600. (a) Recovered conductivity model from the inversion of 2-D DC data, (b) Susceptibihty model recon- structed from the 1-D inversion of D I G H E M data, (c) Magnetite content in 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. My 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 l n ( r v ) 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 susceptibili- ties. The resultant images had common features and the main feature coincided with a Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 115 400 500 0 100 £ 200 J 300 400 500 12.4 i 1 1 r 12.6 12.8 13.0 X (km) i i I 13.2 13.4 /a l to : IN b K x 1000 SI 12.4 12.6 12.8 13.0 X (km) 13.2 13.4 J 200 H f 300 4 0 0 500 K x 1000 Si 12.4 12.6 12.8 13.0 X (km) 13.2 13.4 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- 20.136 — 10.068 I—L 0.001 74.563 66.278 57.994 49 .709 41 .424 t- 33 .140 24 .855 16.570 8.286 0.001 F i g u r e 4.8: Recove red suscep t ib i l i ty models f rom inversions w i t h different er ror schemes: (a) the recovered m o d e l w h e n the s t anda rd devia t ions were 5 p p m plus 10 percent of the s t reng th of the data ; (b) the s tandard devia t ions were 1,4 a n d 10 p p m for d a t a at 900, 7200 and 56,000 her tz ; a n d (c) a constant s t andard dev ia t i on of 10 p p m was used for a l l the da ta . Chapter 4. Recovering susceptibihty from 1-D Inversion of EM data 116 region of high magnetite content inferred from visual estimates of borehole logs. Reason- ably 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 inform- ation 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 mag- netically 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 time- and 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 non- magnetic. 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 0, then the recovered conductivity model could be very 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 DC 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 DC 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 DC 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 con- ductivity or susceptibihty. Solid lines denote the recovered models, and dashed lines denote the true models, (a) The recovered susceptibihty from the inversion with inac- curate 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 Chapter 5. 1-D Simultaneous Inversion of EM data 120 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. One part, which is related to eddy currents through the term jumper, 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. My 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- aa) J w2(z) d(lna — ln<r0) 8~z dz, (5.1) and <f>* = OLK J w3(z)[m(K,) - m(n0)}2dz + (1 - ctK) j w4(z) d[m(n) — m(/€0)] dz dz, (5.2) Chapter 5. 1-D Simultaneous Inversion of EM data 121 where cr0 and Ko are the reference models for conductivity and susceptibihty. The pa- rameters aa and aK 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 conduct- ivity 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 re- covered 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 4><T = W^lnf—) and 4>* = \\WK[m(K) - rn(K0) (5.3) (5.4) where Wa and WK are M x M weighting matrixes. 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>a and (j>K into my cost function. Let the final cost function be 4>m = Q<I><T + 1<}>K, (5.5) where coefficients g and 7 are given by 1 1 + s 1 + s where 0 < s < 00 is the desired magnifying factor. When s —> 0, (j)m = <t>*, a n d when s —> oo,(j>m = 4>K. Now I solve the simultaneous inverse problem by minimizing the cost (5.6) Chapter 5. 1-D Simultaneous Inversion of EM data 122 function subject to the constraint that the data are adequately reproduced: minimize <f> = (gfa + >y<j>K) + 3 l((f>d - 4>tar), (5.7) where 3~x is a Lagrange multiplier and (f)tar is the target misfit level. In above equation, the data objective function <f>d is the same as used in Chapter 2 [equation (2.27)]. Let ma = ln(<r) and mK = m(n) be vectors with M components, and let Sm^ and 6mK denote perturbation at the nth iteration. The predicted data can be approximated as 5[m[ n ) + Sm^mW + SmK] « F[m<?\m^} + JJm, + JJmK, (5.8) where Ja and JK are the sensitivities whose elements are Jau = ddi/dmai and JRu = ddi/dmKi. Let J = ( J a , JK) be a global sensitivity matrix, m = (mC T,m r e) be a global model parameter, and ( gWa 0 \ W m = (5.9) V 0 1WK ) be the global weighting matrix. The linearized problem becomes the minimization of d> = 8 Wm[8m + -m0] ' + j Wd{Dohs - F[m^] + JSm}\2 - <j>tar} , (5.10) where cj)tar is the target level for data misfit at the (n + l) th iteration. Usually I reduce 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 = fi0(l + K) , and for convenience I use both K and fi in the following discussion. The governing partial Chapter 5. 1-D Simultaneous Inversion of EM data 123 differential equation for the electric field, in the Hankel transform domain, is given by equation (2.5): r d2 1 — - u 2 E(X,z,co) = iLOfi031(Xa)6(zobs - h), (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 Zi- Ei _ Z Z ^ + Zjtanh(ujhj) Hrj 3 Zj + ZJ+1 t a n h ( ^ ) ' ^ ' ; where Zj = —^L- Physically the secondary field results from both induced eddy cur- 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. Con- sider 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 Hz{X,co,zobs) = B o e ^ 2 h - ^ ^ ^ , (5.13) Xfi + ufi0 where u2 = A 2 + jtoficr, and X is the Hankel transform variable. B0 is given by 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 Hz(X,u,,zobs) = Boe-^2h-^\^. (5.15) X + u Chapter 5. 1-D Simultaneous Inversion of EM data 124 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(jJlioHra — ju/poa, where a — \ira. 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: AHz(\,u,zobs) = Boe-**-^ (^_^ _ ( 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 g ( . ,^) = - ^ f ^ ) 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,fi0), D(a,u,0), and AD(a,fi). Note that \iTa = ex. 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), D(cr, /x 0), D(a,fi0), and /L\D(a,fi) in Figure 5.2. Table 5.1 shows D(cr,/i), Z?(<T, i x 0 ) , D(a,/A0), and AZ) (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 fir > 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,0) by about 3% to 10%, and 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,fi0)\ is about 400% at 900 hertz, and about 14% at 56k hertz. These effective magnetic charges increase the amplitude of ImD(cr, fi0) from a Chapter 5. 1-D Simultaneous Inversion of EM data 126 Hz Real(Hz) (ppm) Imag(Hz) (ppm) Freq (hertz) 900 7200 56000 900 7200 56000 D(*,fi) -347 171 2362 220.1 970.9 2115 D((T, fXo) 52.9 528.1 2595 202.9 908.5 2021 D(a,ft,0) 59.5 578.3 2737 219.5 959.8 2052 AD(*,p) -406 -407 -375 0.6 11.1 63.1 Table 5.1: D(o~, po), D(a,^0), AD(a,fi), and D(o~,fi) over a 0.01 S/m and 0.1 SI unit half-space. 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 + JKSK, (5.19) where and JK are the sensitivities for conductivity and susceptibihty, whose elements are EtiWu 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 8DK((T, K) W [{J<r)liS<Ti + (JK)li6Ki\ 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 com- posed 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 in- dependent 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 128 Real Imag Freq (hertz) 900 7200 56000 900 7200 56000 dHz/d<r(xlO-r) -5.71 -43.6 -122 -14.4 -44.3 -26.8 dHz/dK(xlQ-3) 240 210 140 -1.0 -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 re- covered 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. Chapter 5. 1-D Simultaneous Inversion of EM data 129 mapping parameters K i and Kb are set to 10~3 and 10~6. The parameters aa and aK are 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 mS/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. By 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) suscep- tibihty 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 para- meter 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 Chapter 5. 1-D Simultaneous Inversion of EM data 132 700 5001 | >" 300| 1001 0 G O Ko T ~ T ~ 1 — I I I l I T " T ~ 1 ~ - 1 . 4 - - 1 I i - + - + . 1 _ _ r . I t • - t - i I 0 15 30 1 N 50 100 X(m) 100 300 500 700 100 300 500 X(m) X-Y Plan-view 700 ao Ko l 1 X-Z Cross Section Figure 5.6: The 3-D model. The background conductivity and susceptibility cr0 and K0 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, <r2 = 0.5S/m and K2 = 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. 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 neg- ative 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 im- possible 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 rep- resents 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. 1-D Simultaneous Inversion of EM data 134 Predicted data (ppm) Observed data (ppm) > 600 400 200 200 400 600 600 — 400 > 200 f=7040 |b 200 400 600 200 400 600 X(m) 600 400 200 600 - 400 H 200 H 200 400 600 600 400 200 1 - i ; J i f=56k| i i t 200 400 X(m) 600 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. Chapter 5. 1-D Simultaneous Inversion of EM data 135 Predicted data (ppm) Observed data (ppm) > 600 H 400 i 200 | f=110 -i r 0 200 400 600 600 - 400 i 200 200 400 600 0 200 400 600 X(m) 600 400 200 H f=110 200 400 POO 600 400 200 f=7040 S 200 400 i'Cf 600 400 200 200 400 X(m) I 126.800 116.114 105.429 94.743 84.058 73.372 62.687 52.001 41.316 30.630 1111.000 1068.089 1025.178 982.267 939.356 896.444 853.533 810.622 767.711 724.800 E 2064.00 2011.56 1959.11 1906.67 S - 1854.22 | - 1801.78 1749.33 | - 1696.89 1644.44 1592.00 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. Chapter 5. 1-D Simultaneous Inversion of EM data 136 m ) 50 ~ S — .c +-> CL 100 - CD Q 150 - Conductivity (S/m) 100 200 300 400 X(m) 500 600 700 0.195 0.139 0.100 0.071 0.051 0.036 0.026 0.019 0.013 0.010 F i g u r e 5.9: T h e x -z cross sections of the recovered c o n d u c t i v i t y f rom the inver s ion of 3-D da ta . W h i t e rectangles denote the posi t ions of those two pr i sms , (a) T h e x - z cross sect ion at y = 4 5 0 m ; (b) T h e x - z cross sect ion at y = 3 5 0 m ; (c) T h e cross sect ion at y = 2 5 0 m . Chapter 5. 1-D Simultaneous Inversion of EM data 137 S 50 CL 1 °0 CD • 150 o | 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 S 50 S~ CL 100 CD Q 150 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 50 Q . CD I I M K N X Susceptibility (SI) [c] 100 200 300 400 X(m) 500 600 700 ^ — r i 0.158 0.140 0.123 0.105 0.088 0.070 0.053 0.035 0.018 0.000 F i g u r e 5.10: T h e cross sections i n the x -d i r ec t i on of the recovered suscept ib ih ty . T h e pos i t ions o f the pr i sms are denoted w i t h wh i t e rectangles. Panels (a) , (b) a n d (c) show the cross sections at y=450 , 350, and 250m respect ively. Chapter 5. 1-D Simultaneous Inversion of EM data 138 600 - g - 400 > > E 200 z=20m 600 H 400 200 z=40m 600 400 i 200 1 1 1— 0 200 400 600 1 1 1 0 200 400 600 I I I J 1 ' 1 / z=60m 0 200 400 600 X(m) m- 600 -i 400 200 z=20m 600 400 200 w z=40m; 600 400 H 200 z=60m i d ! — i 1 r - 200 400 600 0 200 400 600 f 0 200 400 600 X(m) 0.172 0.153 0.134 0.115 0.096 0.077 0.058 0.038 0.019 0.000 0.172 0.153 0.134 0.115 0.096 0.077 0.058 0.038 0.019 0.000 0.172 0.153 0.134 0.115 0.096 0.077 0.058 0.038 0.019 0.000 F i g u r e 5.11: T h e x - y p lan-v iews of the recovered c o n d u c t i v i t y a n d suscept ib ih ty . Pane ls (a), (b) a n d (c) show the recovered c o n d u c t i v i t y at z=20 , 40, and 60 meters . Pane l s (d) , (e) a n d (f) present the suscep t ib ih ty models at the cor responding depths . Chapter 5. 1-D Simultaneous Inversion of EM data 139 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 id , 5.l ie, 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 con- ductivity 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 re- covered 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 Chapter 5. 1-D Simultaneous Inversion of EM data 140 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 aK were set to 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 mS/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. 1-D Simultaneous Inversion of EM data 141 Predicted data (ppm) Observed data (ppm) 10800 -gj gi 10600 - O 10400 V ~ — 10200 10000 i; 11000 10800 10200 10000 -+ <2U> 935 hertz ] 12200 12400 12600 12800 13000 13200 a; _ i 4600 hertz 12200 12400 12600 12800 13000 13200 11000 H ' ' 1 - r J 1 * 1 10800 S 10600 c 1 10400 o Z 10200 10000 -f 12200 12400 12600 12800 13000 13200 Easting (m) 12200 12400 12600 12800 13000 13200 60.000 52.111 44.222 36.333 28.444 20.556 12.667 4.778 -3.111 -11.000 12200 12400 12600 12800 13000 13200 Easting (m) F i g u r e 5.12: T h e rea l components of the p red ic t ed and observed da ta . D a r k region d o m i n a t i n g the p ic tu re is the M a i n Zone. Chapter 5. 1-D Simultaneous Inversion of EM data 142 Predicted data (ppm) Observed data (ppm) -|—1—r 12200 12400 12600 12800 13000 13200 11000 10800 — 10600 10200 10000 o S i 4600 hertz b "t 12200 12400 12600 12800 13000 13200 12200 12400 12600 12800 13000 13200 Easting (m) — I — ' —i— ' — I — ' — I — 1 12200 12400 12600 12800 13000 13200 Easting (m) F i g u r e 5.13: T h e i m a g i n a r y components of the p red ic t ed and observed da ta . D a r k region d o m i n a t i n g the p i c tu r e is the M a i n Zone. 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 aeromag- netic 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 anom- alies 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. 1-D Simultaneous Inversion of EM data 144 Conductivity Model (mS/m) 11000 1 0 8 0 0 1 0 6 0 0 1 0 4 0 0 1 0 2 0 0 1 0 0 0 0 Depth=30m 5 1 2 2 0 0 1 2 4 0 0 1 2 6 0 0 1 2 8 0 0 E a s t i n g (m) 1 3 0 0 0 1 3 2 0 0 1 0 0 0 0 1 0 2 0 0 1 0 4 0 0 1 0 6 0 0 E a s t i n g (m) 1 0 8 0 0 1 1 0 0 0 F i g u r e 5.14: Recove red c o n d u c t i v i t y models f rom the s imul taneous inve r s ion of H S S A E M da ta , (a) T h e x - y p l a n v i ew of the recovered c o n d u c t i v i t y at z = 3 0 m ; (b) T h e c o n d u c t i v i t y a long the sect ion denoted w i t h the whi t e hne i n pane l (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 con- ductivity. 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. Chapter 5. 1-D Simultaneous Inversion of EM data 146 e o 11000 10800 10600 10400 10200 10000 V 11000 10600 10400 3 9 3 — z=100m 12200 12400 12600 12800 13000 13200 12200 12400 12600 12800 13000 13200 E 11000 10800 10600 O 10400 H z 10200 - 10000 12200 12400 12600 12800 13000 13200 J I 1 1 1 1 1 1 12200 12400 12800 12800 13000 13200 ' 1 ' 1 1 1 1 1 ' 12200 12400 12600 12800 13000 13200 11000 10800 10600 10400 10200 10000 1 1 1 1 1 1 1 1 r 12200 12400 12600 12800 13000 13200 Eas t ing (m) Eas t i ng (m) F i g u r e 5.15: T h e x - y p l anv i ew of the recovered suscep t ib ih ty models f rom the s imu l t a - neous inver s ion of A E M da t a and the invers ion of aeromagnet ic da ta . Pane ls (a) , (b) , a n d (c) show the recovered suscep t ib ih ty f rom the inver s ion of aeromagnet ic da ta , a n d panels (d) , (e), and (f) plot the suscep t ib ih ty recovered f rom the s imul taneous inve r s ion . Chapter 5. 1-D Simultaneous Inversion of EM data 147 Line 12600 at HSS 50 4 g 100 £ 150 I" 200 250 300 1° Susceptibility (mSI) 10000 10200 10400 10600 50 100 a O 250 -\ 300 10000 10200 10400 10600 Northing (m) 10800 10800 11000 Susceptibility from mag3d (mSI) 11000 50 45 40 35 30 25 20 15 10 5 0 Figure 5.16: The cross-section of the recovered susceptibihty models from the simultane- ous 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. Chapter 5. 1-D Simultaneous Inversion of EM data 148 5.5 Summary Formulation of the simultaneous inverse problem requires minimization of an objective function including both conductivity and susceptibihty <J)=Q4><7 + 7</>K. By choosing dif- ferent weighting, I can obtain different solutions. Choice of the final relative weighting between (j)a and 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 sus- ceptibihty, 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 invert- ing for conductivity or susceptibility separately because the two model parameters are strongly coupled. Magnetic permeability \i affects the data in two ways. In the quasi- static 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 about suscep t ib i l i ty d i s t r i b u t i o n . 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 qualit- ative 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 Chapter 6. Approximate inversion of 3-D EM data 151 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 coopera- tive 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 approxi- mate, 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 approxi- mation (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 solv- ing 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 computa- tion 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. Chapter 6. Approximate inversion of 3-D EM data 153 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 th datum with respect to unit change in the jth. model cell. By neglecting the higher order term, and setting d[m^ + m) = doh3, where doba are the observations, I can introduce a data misfit objective function: <l)d=\\Wd(Jm-d)\\\ (6.2) where Wd is an N xN matrix whose diagonal elements are the reciprocals of the estimated error of each datum, and d = dohs — d(m^) + Jm^n\ The third step is to solve a linear inverse problem. The model objective function I want to minimize is <t>m{™>) = a i J ws(m - m0)2dv + a2 J W-, d{m — m 0 ) + « 3 d(m — m0) dy dv + CI4 J Wz dx d(m — m0) dz dv dv. (6.3) where m0 is the reference model, alya2,a3, and a 4 control the importance of each term. ws,wx,Wy and wz 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 Chapter 6. Approximate inversion of 3-D EM data 154 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 - m0)TW^Wm(m - m0), (6.4) and W^Wm is given by W*Wm = atWjWs + a2WlWx + a3WfWy + a4W?Wz. (6.5) In equation (6.5) Ws is a diagonal matrix with elements y^A^xAyAz, Wx = ^jAyAz/dxWn, Wy = AxAzjdyWD, and Wz = Ax Ay / dzWp, where Aa:, Ay, and Az are the lengths 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> = 4>m + B~\<f>d - 4>tar), (6.6) where (j)tar is the desired misfit level. Even for a small 3-D problem, the number of cells can be easily over 105, and number of data usually is greater than 103. Thus it is impractical to solve the linear inverse 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) Chapter 6. Approximate inversion of 3-D EM data 155 The linearized objective function in equation (6.6) can be rewritten as <f>m(™) = \\Wm[m^ + Va- m0}\\2+(3-l{\\Wd[d(m^) + JVa - d°bs]\\2 - 4>tar}. (6.8) Setting V(/)(a) — 0 yields Ba = c (6.9) where B and c are given by B = VT(JTWjWdJ + 3WlWm)V, (6.10) and c = -BVTWlWm[m^ - m 0] - VTJTWjWd[d(m^) - 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>d and (f>m, are used as basis vectors. The misfit objective function may be partitioned in a variety of ways. In my work, I partition <j)d according to the frequencies: 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)m corresponding to the four components in equation (6.6). So the next set of vectors is vk = {WlW^V^l, (6.13) where k = 1,2,3,4. The sum of all four terms in above equation equals the steepest descent vector for Vm<^>m. A constant vector is also included. So the total number of vectors used in the subspace inversion is NF + 4 + 1. Chapter 6. Approximate inversion of 3-D EM data 156 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 Gil l 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 = - j w i i H + J m (6.14) V x H = ((T + j W e ) E + J e ) (6.15) where J e is the electric current, and J m = — 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 M <r(R) = 5>*Mir), (6.17) fc=i where ipk is the basis function, then the sensitivity problem is „ <9E . <9H , v x 5S = -""aS ( 6 ' 1 8 ) V x g - = („ + „ * ) _ , (6.19) Chapter 6. Approximate inversion of 3-D EM data 157 and the corresponding boundary conditions are n{n x <9U ) + i{n x n x V x 8U ) = o. (6.20) d(Tk The auxiliary problem associated with above sensitivity problem is V x G = -jufiH* + 6(r - r 0) (6.21) V x H* = (<T + jue)G (6.22) where r 0 is the observation point, and G and H* are the auxiliary electric and magnetic fields. The sensitivities for the vertical magnetic field are 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 (6.23) where G is due to a vertical magnetic dipole of strength (—itou) 1 at r 0 . 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. My 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 = -iufi0Ia AJi(Aa)J!(Ar)dA, (6.24) JO U0 + Ui 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 . . . .A. i . 1 n 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 = Exnx + EyTiy-hEziix, and G = Gxnx-\-Gyny + Gznx, where n x , n y, and nz are the unit vectors in the x-, y-, and z-directions, then the real and imaginary components of the inner product are Re(E • G ) = ReEx ReGx - ImEx ImGx, +ReEy ReGy - ImEy ImGy ( 6 . 2 6 ) and Im(E • G ) = ReEx ImGx + ReGx ImEx + ReEy ImGy + ReGy ImEy. ( 6 . 2 7 ) Note that in 1-D Maxwell's system of equations Ez — 0. Since G satisfies the same partial differential equation (except a scaling factor) and boundary conditions as E , Gz also equals zero. If the amplitudes of E and G are denoted as E and G, then E& x<0,y<0 ( 6 . 2 8 ) £ f f x < 0, y > 0 E r Chapter 6. Approximate inversion of 3-D EM data 160 a G \y-R\ x < 0, y < R and Ey = E x < 0 Gy = G— x<0 (6.29) (6.30) where rp = y/x2 + y2 and ra = ^(y — R)2 + x2. So equations (6.26) and (6.28) have to be calculated accordingly in regions I, II, and III in Figure 6.2: Im(E • G ) p\y\\y-R\W y < Q p-y\y-R\W 0 < < R rpra — 3 py{y-R)W > R (6.31) rvTa and Re(E • G) = < Q \y\\y-R\+x2 y < 0 Q-y\y-R\+*2 o<y<R Qy(y-R)+x (6.32) 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) = E • G , then f(x,y,z) = (l-t)(l-u)(l-v)f1 + t(l-u)(l-v)f2 + tu(l-v)f3 +{l-t)u{l-v)fi + (l-t)(l-u)vf5 +t(l - u)vf6 + tuvf7 + (1 - t)uvfs, (6.33) where / i to f& are the values of f(x,y, z) at the eight nodes, and t = x — X\ x2 — X\ u y-yi , z-z! , and v = yi - yi Z2 - Zl (6.34) Chapter 6. Approximate inversion of 3-D EM data 161 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 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 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 8 (6.35) over regions of several hundred square-kilometers. When the earth is discretized into Chapter 6. Approximate inversion of 3-D EM data 162 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 z- direction. 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~5, a2 = Chapter 6. Approximate inversion of 3-D EM data 163 F i g u r e 6.4: T h e rea l a n d i m a g i n a r y components of the sensi t ivi t ies due to a layer be- tween z = 0 and 2 m d e p t h ins ide a 0.01 S / m half-space, a) T h e rea l componen t of the sensi t iv i t ies , b) T h e i m a g i n a r y component of the sensi t ivi t ies . Chapter 6. Approximate inversion of 3-D EM data 164 F i g u r e 6.5: T h e rea l [panel (a)] and i m a g i n a r y [panel (b)] components of the sensi t iv i t ies due to the layer ex t end ing f rom z = 45 to z = 5 0 m ins ide 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. a3 = a 4 = 1.0. The initial chi-square misfit was 5288. After 20 iterations, the misfit was 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 impos- sible 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. At 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 Mt. 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 Chapter 6. Approximate inversion of 3-D EM data 167 Apparent Conductivity (mS/m) Recovered Conductivity (mS/m) 6 0 0 4 0 0 2 0 0 0 - i- 7 0 . 0 - 6 3 . 1 - 5 6 . 2 - 4 9 . 3 - 4 2 . 4 - 3 5 . 5 - 2 8 . 6 - 2 1 . 7 - 1 4 . 8 - 7 . 9 _ n — L 1 .0 0 fill ) III z = 2 0 m [dj 2 0 0 4 0 0 6 0 0 2 0 0 4 0 0 6 0 0 6 0 0 4 0 0 2 0 0 6 0 0 4 0 0 H 2 0 0 2 0 0 4 0 0 E a s t i n g ( m ) 2 0 0 4 0 0 E a s t i n g ( m ) 9 8 . 0 8 8 . 7 7 9 . 4 7 0 . 1 6 0 . 8 5 1 . 5 4 2 . 2 3 2 . 9 2 3 . 6 1 4 . 3 5 . 0 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 Recovered Conductivity (mS/m) o 50 100 150 JS / 1 . f y=250m \ a 100 200 300 400 500 600 700 700 100 200 300 400 500 600 700 0 50 100 150 y=350m] 100 200 300 400 E a s t i n g (m) 500 600 I 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 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. Chapter 6. Approximate inversion of 3-D EM data 169 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 x = 2 x 10~5, and a 2 = a 3 = a 4 = 0.1. The 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 105. 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 104. After 20 iterations, the chi-square misfit 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 DC 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 geo- logical 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 conduc- tivity 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 104, to investigate the effect of different target 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 171 Real (ppm) Imag(ppm) 10.0 10.0 ~i 1 r 12.4 12.6 12.8 13.0 13.2 13.4 Easting (km) 12.4 12.6 12.8 13.0 13.2 13.4 Easting (km) 10.0 12.4 12.6 12.8 13.0 13.2 13.4 Easting (km) 12.4 12.6 12.8 13.0 13.2 13.4 Easting (km) 12.4 12.6 12.8 13.0 13.2 13.4 Easting (km) 12.4 12.6 12.8 13.0 13.2 13.4 Easting (km) Figure 6.9: D I G H E M data at Mt. 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. Chapter 6. Approximate inversion of 3-D EM data 172 Recovered Conductivity (mS/m) . — , _ . , R 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) 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) "1 1 ! 1 i 1 1 "1 1 r 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) 0 1 50 £ 100 J " 150 200 J L 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 13.5 Easting (km) 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 Mt. 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 DC resistivity data over section Y9600. Chapter 6. Approximate inversion of 3-D EM data 173 Recovered Conductivity (mS/m) -i r 1 2 . 4 1 2 . 6 1 2 . 8 1 3 . 0 1 3 . 2 1 3 . 4 I I 1 2 . 4 12 .6 i l l 1 2 . 8 1 3 . 0 1 3 . 2 1 3 . 4 7 E a s t i n g ( k m ) E a s t i n g ( k m ) 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 Mt. Milligan. Panels (a) to (f) plot the sliced model at z=10, 30, 50, 80, 100, and 120m. Chapter 6. Approximate inversion of 3-D EM data 174 Recovered Conductivity (mS/m) o g 50 T 100 CL G Q 200 ! . . . . . . . . . . } y J y=9.2km _y i a 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) 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?r r i 4 •K3? ) R a i n b o w y=9.6km M B X Stock / D y k e d 1 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) 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 Mt. 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 DC resistivity data over section Y9600. Chapter 6. Approximate inversion of 3-D EM data 175 T h e choice of target misf i t is i m p o r t a n t for th is inverse m a p p i n g . T o keep the l i n - ea r i za t ion v a l i d , the d a t a misfi t shou ld not be reduced too m u c h . T h e appropr i a t e target misf i t l eve l m a y be de t e rmined by 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 images f r o m several test runs . Because of the 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 needed w h e n a p p l y i n g this inverse m a p p i n g to d a t a col lec ted f rom regions w i t h large c o n d u c t i v i t y contras t . T h e o u t p u t f r o m this l inear m a p p i n g m a y be geological ly a n d geophys ica l ly in ter - pre tab le . I f the inverse m a p p i n g can be done reasonably fast t hen i t w i l l p rov ide p r o m p t feedback i n f ie ld surveys . T h i s technique allows the d a t a to be presented i n the r ight p h y s i c a l un i t s . I f there are a r t i f i c ia l s t ructures , the in te rpre te r w i l l h k e l y come to recog- nize a n d account for t h e m i n any p r e h m i n a r y in t e rp re ta t ion . T h e o u t p u t f r o m this hnear m a p p i n g m a y also be used as a weight ing func t ion i n the invers ion of another geophys ica l d a t a set. T h i s can be a first pass coopera t ive invers ion . T h e image o b t a i n e d f r o m the hnear m a p p i n g can also be used as a s ta r t ing m o d e l for subsequent invers ions . T h e results f rom the syn the t ic da t a test are encouraging . T h e resul tant images have p r o v i d e d m u c h more i n f o r m a t i o n about the dep th ex tens ion 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 conduc t iv i t y . T h e test on the f ie ld d a t a set also generated useful images w h i c h not on ly agree w i t h tha t ob t a ined f r o m 2 -D r igorous inver s ion of D C res i s t iv i ty da ta , bu t also are geological ly meaningfu l . T h e exis tence of the three ma jo r geological uni t s is obvious i n the recons t ruc ted c o n d u c t i v i t y image . 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. By 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 dif- ferent 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 sta- ble 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. Summary 179 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 as- sumption 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 sus- ceptibihty 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 sus- ceptibihty 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 situ- ations 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. 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 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 controlled- source 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 tha t can generate i n f o r m a t i o n about conduc- t i v i t y and suscep t ib ih ty f rom the rigorous s imul taneous invers ions of cont ro l led-source E M d a t a f rom different co i l systems for the first t ime . These techniques are po ten- t i a l l y useful to the m i n i n g i n d u s t r y and i n other places where magne t i c suscep t ib ih ty is a conce rn . I have also recons t ruc ted c o n d u c t i v i t y f rom the invers ions o f t rans ient E M borehole da ta , w h i c h can be s t i l l considered as a new technique . T h o s e 1-D inve r s ion me thods deve loped i n this thesis have p r o v i d e d necessary p h y s i c a l ins ights needed to a t t ack h igher d i m e n s i o n a l p roblems. T h e 3-D l inear m a p p i n g can be p o t e n t i a l l y useful i n the i n t e rp re t a t i on a n d invers ion of control led-source E M da ta . References [1] Anderson, W.L. , 1979, Computer Program: Numerical integration of related Han- kel 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 SEG 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 mag- netic 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, 279- 288. [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 1991- 1, 199-205. [9] DeMoully, G. T., and Becker, A. , 1984, Automated interpretation of airborne elec- tromagnetic 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, 384- 404. .References 184 [15] Elliot, C. L. , 1961, A electromagnetic drill hole instrument for the detection of conductive sulfide bodies: Presented at the 31th SEG Ann. Internat. Mtg. - 1966 Electromagnetic method and apparatus of geophysical exploration: Cana- dian patent 743,665. [16] Ellis, R. G., 1995, Joint 3D E M inversion: Abstract for 3D electromagnetics work- shop 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 electro- magnetic 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 electromag- netic 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 mini- mization of a special functional, Geophysics, 49, 1354-1360. [31] Gupta Sarma, D., Maru., U. M . , and Varadarajan, G. , 1976, An 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., Elec- tromagnetic 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 pro- specting, Pergamon Press, New York, NY, 519p. [36] Last, B. J. and Kubik, K., 1983, Compact gravity inversion, Geophysics, 48, 713- 721. [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 struc- tures by simple processing of transient electromagnetic data: Geophysics, 52, 545- 554. [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: Geoexplo- ration, 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, Vol- ume 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 electromag- netic 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 electro- magnetic 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: Geo- physics, 52, 1431-1435. [50] Noakes, J.E., 1951, An electromagnetic method of geophysical prospecting for ap- plication 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 hydroge- ological problems: Geophys. Prospect., 18, 236-254. [57] Qian, W., Gamey, J., Lo, B., and Holladay, J.S., 1996, AEM apparent resistiv- ity 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 exper- iment across Santa Clara Valley: Geophysics, 37, 351-374. [61] Sengpiel, K . P., 1988, Approximate inversion of airborne E M data from a multi- layered 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: Phys- ics, 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 . , Ed. , Electromagnetic Methods in Applied Geophysics: Inves- tigations 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 model- ing 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, Elec- tromagnetic 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, NY, 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 ho- mogeneous earth. Trans. Am. Geophys. Union, 34, 185-188. -1955, Mutual electromagnetic coupling of loops over a homogeneous ground: Geo- physics, 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, SEG publication, 131-311. [79] Ward, S.H., Ryu, J . , Glenn, W . E . , Hohmann, G.W. and Smith, B.D., 1974, Elec- tromagnetic 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, SEG 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 SEG annual meeting, 241-244. [85] Zhang, Z., Routh, P., and Oldenburg, D.W., 1996, Reconstruction of 1-D con- ductivity from coaxial, coplanar, vertical coplanar and perpendicular E M data: Expanded Abstract of 67th SEG 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 quasi- linear approximation: The abstracts for the international symposium on 3-D elec- tromagnetics 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 UmHz(r^zoba) = £ jT | _ ^ A > e ^ - » f c ) j 0 ( A r ) d A . (A.l) For a half space of fi = fii and a = a\, the input impedance equals the intrinsic impedance Z1 = Z L y / \ 2 — U>2€ofll - f ' i w/iiO -! When the frequency tends to zero, the normalized input impedance is A. -ILO Hi. Thus expression (A.l) can be further reduced to _ m (m - /xp) w—>0 l i m f l . C r . w , * ^ ) = ^-fL^El f*~ A V ^ - 2 h > J 0 ( A r ) d A . 47T (/i! + /i„) Jo From Gradshteyn and Ryzhik (1965, p712 ), poo J xm+1e-axJv(bx)dx = (-l)m+1b~v dm+l (Va2 + b2- a)v Va2 + b2 dam+1 Letting m = 1, b = r,v = 0, a = 2h and z o b s = 0 yields ; v IT r \ m _ /*o) 8h2 - r2 km HAr.u), zob.) = — ^ ^—r-, . „_<, ' ° w 4 f f ( ^ 1 + / i 0 ) ( 4 / l 2 + r 2)2 .5 ' (A.2) (A.3) (A.4) (A.5) (A.6) 193 Appendix B Adjoint Green's function solution for the sensitivities From the first two equations in equation (4.5) I obtain d2 8E . d dHr . dp dHr . dHr , . T T dtbi dz2 dpi dz dpi dz dpi dz dz and ( B . l ) -{-dr I r d V dpij d dH 9HZ (B.2) dr dpi " z dr ' 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 IUJ/J,, I obtain . d dHr . d dHz . . a . x dE t W t i ¥ z T p - ~ l^-dr-dp- = ^ + ^ d p - - (B.3) Solving equations (B.l) and (B.2) for iup-f^ and iup£^, and inserting these into equation (B.3) leads to d2 dE (. „di>i d n ; iu>Hr— 1- tu)ipi dz2 dpi dz dpi Above equation can be reorganized into d2 dE d ( 1 + —I- dr I r d_( dET dr \ r d p i j + tu-^—ifi or = (—iuj2 jj,ie0 + iuia) dE dpi' dz2 dpi dr d ( dEs dr \ dpit (—ia}2pe0 + iua) dE dpi (B.4) (B.5) 194 Appendix B. Adjoint Green's function solution for the sensitivities 195 . , fdHr dHz\ . dibi , N Since there is no artificial source in the ith layer, dHr dHz dz dr = (iioe0 + cr)E. (B.7) Hence the sensitivity problem can be expressed as d2 dE d [ 1 dz2 dm f9r I r d ( dE\]\ . 2 . dE = iwi>i(iwe0 + cx)E + iwHr^. (B.8) Since the box car basis function is actually the subtraction of two Heaviside functions, ipi = H(z — Zi) — H(z — zi+i), the derivative of the box car with respect to the depth is din dz = 8{z - z^ - S(z - zi+1). (B.9) Due to its symmetry this problem can be converted into the Hankel domain (Appendix B ) d2 dE 2dE . , N r , i dErc. , c . _ i n , where u2 = A 2 — o>2/zen + iwfio~. Compared with the partial differential equation for 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 L ~ d f1 d ( dE\ dr {r [dr \ dm}\ rJi(\r)dr = / Jo °° ( d 2 dE d BE IdE l r „ or2 dm dram rdm/ (Cl) Integrate by parts to obtain d 3E + /o rcV [ 7 , J l ( A r ) ] + ^ 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 2 ^ r ) flJ^Ar) J t ( A r ) r - r 5 r 2 J l ( A 7 , j + 2 _ 5 ; (9r r d2 _ . c9Ji(Ar) „ / % N , .(C.3) Let Xr - R, such that ^ = ^ f ^ r = A ^ . Carrying out this replacement in cooperation with the definition of the Bessel's function BR2 RdR R2 RX = -3l(R)RX, (CA) 196 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 d2 dE 20E . , . . _ 1 8 E r c f .•" f . . cfy<9tfr x i 2 ^ - = iuipi{iuje0 + <r)E + —K-[S{Z - «,-) - £(z - zi+1)} + w / - ^ . (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 _ 2utVo Hi dz First of all, the generic solution for the electric field and the Green's function in the ith layer take the following form = G(X,w,z)E(X,v,z)dz - 2Aibihi} (D.l) Therefore _ 1 fZi+1 2EQ JZi e 2 uf Ei = Aieu*z-Zi + Bie-^-^ d = a i e U i i z - X i + bie~Ui{z-Zi>> 2Aibihi = / 2Aibidz = — / 2A{Bidz Jzi EQ Jzi dz, (D.2) 1 (dE d. i 2 J zi dz = =r [' [EG - Zi 1 dGdE (D.S) u] dz dz J where E0 = i c ^ o / a J a ^ ^ e " 0 * 0 6 ' . Inserting equation (C-3) back into the right hand side of equation (D.3) results in 2 u | / G(\,uj,z)E(X,u;,z)dz-2Aibihi =^ Jzi J pii JZi GE 1 = 1 f + 1 (u 2 2 ^ dGdE\ , 2 O E + — — I <Zz 2 « 2 dz c9z 1 / d 2 £ a c ? ^ <iz <9z & = 1/ G dz Integrating by parts I obtain 1 /•«.-+! fgd^E dGdE\ dz_]L dz2 dz dz J m G dE dz JZi 5z <9z f Hi J'i "+1 5Ga/£ dz <9z" (D.5) 198 Appendix E Computation of the discrete sensitivities In the calculation of the sensitivities, I need to evaluate an integral / , defined by /(A,w) = P1 G(X,u,)E(XtU)dz = E\X,u)dz, (E.l) J Zi H/Q Jzi with E0 defined as in Appendix D. The partial differential equation for the electric field in the Hankel-transformed domain is d2E dz2 u2 = 0. (E.2) Hence . . . v 1 r*^Ed2E, 1 {f„dE\ti+1 ft+i (dE\2 } , v 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 rzi+i (8E\ , rzi+i , /. (&) d z = L (E.4) Inserting this equation back into equation (E.3) yields f(X,u,z) Since at the ith interface 1 EJ-z = 2E0u2 dEj _ m dEi_x dz fii-i dz ' + 2AiBihi. (E.5) (E.6) 199 Appendix E. Computation of the discrete sensitivities 200 equation (E.l) can be, for the convenience of computation, further written as 2E0u2 Ei dEj 8z 2.-+1 Ei-i dE^ + 2AiBihi, (E.7) where Ei dEj dz t+i 2BjZi+1 iumBi Z^ + ZiZ*1 + Zi -2u,hi (E.8) Appendix F The trade-off between u and fi in the sensitivities Introducing a new parameter r(pi, ô ) = jajm&i, I can rewrite the perturbation on the datum in equation (5.20) as M s n an M Since dr/dm = dr/dcri x <Ti/m, dDi dr , dDi dr 'dDi dr dDj,dr_s dA. ^ t dr dcri dr dm * dm \ + — —8 i = [Sen + ^8 • dr do~i dr dm * dr dcri \ 1 m / Noticing that dDi dr dr do~i I can further reduce equation (F.l) to (Ja)li, M / a- \ dD £ £ > ; ( < ? - , a ; , z o b 4) « Y\{J„)n [8^ + —8m j + Pi- i=i \ p* J dm (F.l) (F.2) (F.3) (F.4) This is the same as equation (5.20), if dDi/dm is defined as Qu. 201
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
Ukraine | 53 | 0 |
France | 10 | 0 |
China | 10 | 21 |
Australia | 2 | 0 |
Japan | 1 | 0 |
United Kingdom | 1 | 1 |
New Caledonia | 1 | 0 |
United States | 1 | 1 |
City | Views | Downloads |
---|---|---|
Unknown | 67 | 1 |
Wuhan | 5 | 0 |
Beijing | 4 | 0 |
Changchun | 1 | 0 |
Tokyo | 1 | 0 |
Ashburn | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Share to: