- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Approximate sensitivities for the multi-dimensional...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Approximate sensitivities for the multi-dimensional electromagnetic inverse problem 1995
pdf
Page Metadata
Item Metadata
Title | Approximate sensitivities for the multi-dimensional electromagnetic inverse problem |
Creator |
Farquharson, Colin Glennie |
Date Created | 2009-04-16 |
Date Issued | 2009-04-16 |
Date | 1995 |
Description | Inversion of measurements from a geophysical electromagnetic survey to produce a two- or three-dimensional conductivity model of the Earth is computationally demanding. The classical approach is to linearise the inverse problem and iterate towards the solution. A typical, modern version of this approach is described in the initial part of this thesis, and applied to the one-dimensional inversion of time-domain electromagnetic data. One of the most time-consuming parts of a linearised, iterative inversion procedure is the generation of the Jacobian matrix of sensitivities. In this thesis, I have developed an approximate method for generating these sensitivities. The approximation is based on the adjoint-equation method in which the sensitivities are obtained by integrating, over an individual cell, the scalar product of an adjoint field (the adjoint Green’s function) with the electric field produced by the sources for the geophysical survey. Instead of calculating the adjoint field in the multi-dimensional conductivity model, an approximate adjoint field is computed in a homogeneous or layered halfspace, or using the Born approximation. This approximate adjoint field is significantly quicker to compute than the true adjoint field and leads to considerable reductions in the time required to generate the Jacobian matrix of sensitivities. The time-differences were found to be of one or two orders of magnitude for the small examples considered in this thesis, and this saving will further increase with the size of the problem. The approximate sensitivities were compared to the accurate sensitivities for two and three-dimensional conductivity models, and for both artificial sources and the plane wave source of magnetotellurics. In all the examples considered, the approximate sensi tivities appeared to be sufficiently accurate to allow an iterative inversion algorithm to converge to an acceptable model. This was emphasised by the successful inversion, using approximate sensitivities, of two sets of magnetotelluric data: a synthetic data-set and a sub-set of the COPROD2 data. |
Extent | 3975780 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-16 |
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.0052949 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of |
Degree Grantor | University of British Columbia |
Graduation Date | 1995-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/7282 |
Aggregated Source Repository | DSpace |
Digital Resource Original Record | https://open.library.ubc.ca/collections/831/items/1.0052949/source |
Download
- Media
- ubc_1995-059545.pdf
- ubc_1995-059545.pdf
- ubc_1995-059545.pdf [ 3.79MB ]
- ubc_1995-059545.pdf
- Metadata
- JSON: 1.0052949.json
- JSON-LD: 1.0052949+ld.json
- RDF/XML (Pretty): 1.0052949.xml
- RDF/JSON: 1.0052949+rdf.json
- Turtle: 1.0052949+rdf-turtle.txt
- N-Triples: 1.0052949+rdf-ntriples.txt
- Citation
- 1.0052949.ris
Full Text
APPROXIMATE SENSITIVITIES FOR THE MULTI-DIMENSIONAL ELECTROMAGNETIC INVERSE PROBLEM by Colin Glennie Farquharson B.Sc. (Geophysics), University of Edinburgh, 1990 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF GEOPHYSICS AND ASTRONOMY We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA August 1995 © Colin Glennie Farquharson, 1995 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. (Signature) a Department of o4’I-tYSuS T)2OiOry The University of British Columbia Vancouver, Canada Date i4LL-r ‘t’1’1 DE-6 (2)88) Abstract Inversion of measurements from a geophysical electromagnetic survey to produce a two- or three-dimensional conductivity model of the Earth is computationally demanding. The classical approach is to linearise the inverse problem and iterate towards the solution. A typical, modern version of this approach is described in the initial part of this thesis, and applied to the one-dimensional inversion of time-domain electromagnetic data. One of the most time-consuming parts of a linearised, iterative inversion procedure is the generation of the Jacobian matrix of sensitivities. In this thesis, I have developed an approximate method for generating these sensitivities. The approximation is based on the adjoint-equation method in which the sensitivities are obtained by integrating, over an individual cell, the scalar product of an adjoint field (the adjoint Green’s function) with the electric field produced by the sources for the geophysical survey. Instead of calculating the adjoint field in the multi-dimensional conductivity model, an approximate adjoint field is computed in a homogeneous or layered halfspace, or using the Born approximation. This approximate adjoint field is significantly quicker to compute than the true adjoint field and leads to considerable reductions in the time required to generate the Jacobian matrix of sensitivities. The time-differences were found to be of one or two orders of magnitude for the small examples considered in this thesis, and this saving will further increase with the size of the problem. The approximate sensitivities were compared to the accurate sensitivities for two and three-dimensional conductivity models, and for both artificial sources and the plane wave source of magnetotellurics. In all the examples considered, the approximate sensi tivities appeared to be sufficiently accurate to allow an iterative inversion algorithm to converge to an acceptable model. This was emphasised by the successful inversion, using approximate sensitivities, of two sets of magnetotelluric data: a synthetic data-set and a sub-set of the COPROD2 data. 11 Table of Contents Abstract.ii Table of Contents iii List of Tables vi List of Figures vii Acknowledgments ix Dedication x Chapter 1 Introduction 1 Chapter 2 One-Dimensional Inversion of Time-Domain Electromagnetic Data 7 2.1 Introduction 7 2.2 The forward problem 8 2.3 Calculation of the sensitivities 10 2.4 The inversion algorithm 13 2.5 Examples 17 2.5.1 Synthetic data 17 2.5.2 Field data 22 2.6 Conclusions 27 Chapter 3 Introduction to Approximate Sensitivities 30 3.1 Multi-dimensional problems 30 3.2 The need for approximate sensitivities 32 3.3 Previous forms of approximate sensitivities 34 3.3.1 The one-dimensional problem 34 3.3.2 The two-dimensional problem 35 111 3.3.3 The three-dimensional problem.37 3.4 Conclusions 37 Chapter 4 An Approximate Form of Sensitivity for the General Electromagnetic Inverse Problem 39 4.1 Introduction 39 4.2 The adjoint-equation method for accurately calculating the sensitivities ... 39 4.3 An approximate form of sensitivity 43 Chapter 5 Approximate Adjoint Fields for the Two-Dimensional Problem 45 5.1 Introduction 45 5.2 True adjoint field 45 5.3 Approximate adjoint field: homogeneous halfspace 47 5.4 Approximate adjoint field: layered halfspace 50 5.5 Approximate adjoint field: Born approximation 51 5.6 Computation times 52 5.7 Conclusions 53 Chapter 6 Approximate Sensitivities for the Two-Dimensional Magnetotelluric Inverse Problem 57 6.1 Introduction 57 6.2 Calculation of the sensitivities 57 6.2.1 Sensitivities for apparent resistivity and phase 57 6.2.2 Area integration 58 6.3 Comparison of approximate and accurate sensitivities 60 6.3.1 Homogeneous halfspace 60 6.3.2 Two-dimensional model 63 6.4 Computation times 66 6.5 Two-dimensional inversion of magnetotelluric data 68 6.5.1 Synthetic data 68 iv 6.5.2 COPROD2 data.73 6.6 Conclusions 79 Chapter 7 Approximate Sensitivities for the 25d and Three-Dimensional Inverse Problems 80 7.1 Introduction 80 7.2 Calculation of the sensitivities for the 2.5d problem 81 7.3 Comparison of approximate and accurate 2.5d sensitivities 83 7.3.1 Homogeneous halfspace 83 7.3.2 Two-dimensional model 84 7.4 Computation times for the 2.5d problem 86 7.5 Calculation of the sensitivities for the three-dimensional problem 88 7.6 Comparison of approximate and accurate three-dimensional sensitivities 89 7.6.1 Homogeneous halfspace 89 7.6.2 Three-dimensional model 92 7.7 Computation times for the three-dimensional problem 92 7.8 Conclusions 95 Chapter 8 Summary 97 References 100 Appendices 105 A Computation of the Electric and Magnetic Fields in a Layered Halfspace due to a Dipole Source 105 B Adjoint-equation method for the purely two-dimensional inverse problem .. 111 C Electric fields due to two-dimensional dipole sources on a homogeneous halfspace 112 D Sensitivities for the magnetotelluric apparent resistivity and phase 118 v List of Tables 5.1 Example computation times for the various types of two-dimensional adjoint fields 53 6.1 Example computation times for the various types of sensitivities for the two-dimensional magnetotelluric inverse problem 68 7.1 Example computation times for the accurate and approximate sensitivities for the 2.5d inverse problem 88 7.2 Example computation times for accurate and approximate sensitivities for the three-dimensional inverse problem 94 vi List of Figures 2.1 Linear transformations for the id TEM forward problem 10 2.2 Models produced by the id inversions of synthetic data 18 2.3 The synthetic data inverted to give the models in Fig. 2.2 19 2.4 Convergence curves for the id inversion of synthetic TEM data 21 2.5 A TEM sounding acquired during an environmental survey 23 2.6 Final models from the inversion of the data in Fig. 2.5 26 2.7 A second TEM sounding from the environmental survey 28 2.8 The model obtained from the inversion of the data in Fig. 2.8 29 5.1 The two-dimensional conductivity model used throughout this thesis 46 5.2 The mesh corresponding to the model in Fig. 5.1 47 5.3 The amplitudes of the true and approximate adjoint fields 48 5.4 The phases of the true and approximate adjoint fields 49 5.5 The differences between the logarithm of the amplitudes of the adjoint fields .. 55 5.6 The differences between the phases of the adjoint fields 56 6.1 Notation used for the bilinear interpolation 59 6.2 E-polarisation MT sensitivities for the homogeneous halfspace 62 6.3 11-polarisation MT sensitivities for the homogeneous halfspace 64 6.4 E-polarisation MT sensitivities for the model in Fig. 5.1 65 6.5 11-polarisation MT sensitivities for the model in Fig. 5.1 67 6.6 Synthetic MT data inverted using approximate sensitivities 70 6.7 The final model from the inversion of the data in Fig. 6.6 72 6.8 Convergence curves for the inversion of the synthetic MT data 73 6.9 The sub-set of the COPROD2 data inverted using approximate sensitivities 75 6.10 The mesh used in the inversion of the COPROD2 data 77 6.11 The final model produced by the inversion of the COPROD2 data 79 vu 7.1 2.5d sensitivities for the homogeneous halfspace 85 7.2 2.5d sensitivities for the model in Fig. 5.1 87 7.3 Geometry for the three-dimensional example in section 7.6 90 7.4 3d sensitivities for the homogeneous halfspace 91 7.5 3d sensitivities for the 3d conductivity model 93 A.1 Notation for the horizontally layered conductivity model 106 C.1 The coordinate system and geometry for purely two-dimensional problems ... 114 viii Acknowledgements I would like to thank my supervisor Dr. Doug Oldenburg for his advice and encour agement, and for his numerous shots of enthusiasm when my own was flagging. I am also grateful to Drs. Rob Ellis and Yaoguo Li for their advice and assistance on all matters concerning computers, inverse theory and electromagnetic geophysics. I would also like to thank the members of my committee, Drs. Garry Clarke, Richard Froese and Matt Yedlin, for their advice and support. I am grateful to Cliff Candy of Frontier Geosciences, Inc., for providing the field data, the well-log measurements and the output of the TEMEX-GL program that appear in Chapter 2. I would also like to thank Drs. Martyn Unsworth and Art Raiche for the use of their 2.5d and 3d modelling programs. I would like to acknowledge receipt of a Commonwealth Scholarship which enabled me to carry out the work for this thesis, and to do so in such a wonderful part of the world. Finally, thanks to everyone I have known in the Department of Geophysics and Astronomy for making my time here so very enjoyable. ix To my parents. x Chapter 1 Introduction There are many geophysical techniques which use the principles of electromagnetism to investigate the electrical conductivity of the Earth. The techniques involve a source, ei ther natural or artificial, that generates electric currents within the Earth. Measurements are then made of the electromagnetic fields associated with these currents. The behaviour of the currents is dependent upon the spatial variation of the Earth’s conductivity, and so information about the conductivity can be inferred from the measurements of the electro magnetic fields. A compendium of present-day techniques is given by Nabighian (1988). Geophysical electromagnetic techniques were pioneered by the mining industry and underwent considerable development in the period immediately following World War II (Brant 1966). Exploration for mineral deposits remains the main use of these techniques today. However, electromagnetic methods are also used to study both the deep and shallow structure of the Earth’s crust. Techniques such as the magnetotelluric method are used to investigate the large-scale features of the lithosphere (e.g., Jones 1993), and methods using small, controlled sources play a role in engineering and environmental studies of the shallow subsurface (e.g., Heiland 1940; Ward 1990). Before the general accessibility of computers, quantitative interpretation of the mea surements obtained from geophysical electromagnetic techniques was limited to the use of “type” or “master” curves (Grant and West 1965; Telford, Geldart & Sheriff 1990). Such curves showed the electromagnetic fields that would be observed over simple Earth models for certain sources and were obtained either from analytic formulae if such existed or from scale-model experiments. If a type curve could be found that matched to some extent the observations, then the corresponding model could be considered to resemble the region of the Earth under investigation. 1 CHAPTER 1: INTRODUCTION 2 With the coming of accessible computer resources, it became possible to compute the electromagnetic fields over more complex models than those used to produce type curves. Cagniard’s (1953) procedure for calculating the magnetotelluric response of a halfspace made up of uniform, horizontal layers was adapted by Wu (1968) for the nu merical computation of the response for an arbitrary number of layers. Algorithms were also developed for computing the electromagnetic fields generated in a layered halfspace by controlled sources (e.g., Ryu, Morrison & Ward 1970; Scriba 1974). As the speed and memory of computers increased, so did the complexity of the models and sources for which the electromagnetic fields could be computed. Programs were written to compute the magnetotelluric response over two-dimensional conductivity models (e.g., Madden & Swift 1969; Brewitt-Taylor & Weaver 1976; Rodi 1976; Jupp & Vozoff 1977) and to compute the electromagnetic fields generated in such models by controlled sources (e.g., Everett & Edwards 1992; Unsworth, Travis & Chave 1993). These models were described by large numbers of parameters, usually the conductivities of rectangular cells into which the Earth model was divided, and were general enough to represent arbitrary conductivity distributions. Computer programs have also been written to calculate the electromagnetic fields generated in arbitrary three-dimensional models, both for the plane wave source of magnetotellurics (e.g., Jones 1974; Mackie, Madden & Wannamaker 1993) and for con trolled sources (e.g., Gupta, Raiche & Sugeng 1989; Livelybrooks 1993). MDonald & Agarwal (1994) and Newman & Alumbaugh (1995) have developed three-dimensional forward-modelling programs specifically designed for the new generation of massively par allel computers. However, it is not yet possible, given generally available computer tech nology, to carry out forward modelling for the size and complexity of three-dimensional models that one would like. The development of automated interpretation schemes has followed closely behind at each stage in the evolution of forward-modelling programs. With the ability to compute what would be measured over an arbitrary conductivity structure, and to quickly compute the solution to matrix equations, it became possible to apply mathematical techniques to invert the nonlinear relationship between the parameters describing a model and the CHAPTER 1: INTRODUCTION 3 measured values of the electromagnetic fields. The first inversion techniques to be applied, and the ones which have proven to be the most enduring, are variations of the Gauss- Newton technique (Gill, Murray & Wright 1981). Wu (1968) and Jupp & Vozoff (1975) applied such techniques to the inversion of magnetotelluric data to obtain horizontally- layered conductivity models. The basis of this collection of techniques is to linearise the relationship between the model parameters, m, and the observations, d(0I): d(0 = d[mj + i[m] Sm + R. (1.1) The data that would be measured over the Earth if it were equivalent to the model defined by m are denoted by d[m]. The matrix J is the Jacobian matrix of “sensitivities”, or partial derivatives, of the data with respect to the model parameters: [mJ = i = 1,...,M, j = 1,...,N, (1.2) where M is the number of observations and N is the number of model parameters. Assuming the remainder R can be neglected, eq. (1.1) can be inverted, using some reg ularised procedure, to give a perturbation Sm to the model that will (hopefully) bring the predicted data closer to the observations. Successive iterations will then converge to a model that minimises the misfit between the observations and the predicted data, usually measured in terms of the 12 norm. This was the goal of the initial applications of Gauss-Newton techniques. Linearised, iterative inversion procedures for geophysical electromagnetic data have since undergone many refinements. The increased speed of computers now means for ward modeffing no longer constrains possible inversion techniques for the one-dimensional magnetotelluric problem. It is in the context of this problem that Gauss-Newton inver sion procedures for geophysical electromagnetic data have received much of their develop ment (Constable, Parker & Constable 1987; Smith & Booker 1988; Whittall & Oldenburg 1992), and are now concerned not only with convergence and robustness, but also with the nonuniqueness of the inverse problem. These sophisticated forms of the Gauss-Newton CHAPTER 1: INTRODUCTION 4 method involve many more model parameters than there are observations, and attempt to find the model that has the minimum value of some norm subject to the constraint that the observations are reproduced to an appropriate level of misfit. Sophisticated Gauss-Newton procedures have now been used to invert magnetotel luric data to give two- and three-dimensional conductivity models (e.g., deGroot-Hedlin & Constable 1990; Mackie & Madden 1993), and to invert data obtained in controlled-source surveys for two-dimensional models (Unsworth & Oldenburg 1995). However, computer resources restrict the size and complexity of the multi-dimensional problem that can presently be handled, both because of the time and memory requirements of the forward modeffing and because of the time required at each iteration to generate the Jacobian matrix of sensitivities and to invert the linear system of equations. Many other techniques have been applied to the inversion of geophysical electro magnetic data. Because the one-dimensional magnetotelluric problem is the simplest electromagnetic inverse problem it has received the most attention. An excellent dis cussion and comparison of the various techniques applied to this problem is given by Whittall & Oldenburg (1992). The majority of techniques consider the inverse problem as an optimisation problem, as do the Gauss-Newton methods, in which an objective function is to be minimised. The objective function is usually some combination of a data misfit and a model norm. Two such techniques are simulated annealing (Dosso & Oldenburg 1991) and genetic algorithms (Schultz et al. 1993). A considerable amount of work has also been done on the theoretical aspects of the one-dimensional magnetotel luric inverse problem, both to investigate the existence and uniqueness of solutions, and to devise means of constructing solutions (e.g., Bailey 1970; Weidelt 1972; Parker 1980). The alternatives to the linearised, iterative Gauss-Newton procedures are generally not as efficient and so have not yet been considered for the two- and three-dimensional electromagnetic inverse problems. And since even the Gauss-Newton algorithms are re stricted by current computing power, the only feasible approaches to multi-dimensional inversion of geophysical electromagnetic data are those involving some form of approxima tion. Most of the methods that have been developed are based on a linearised procedure CHAPTER 1: INTRODUCTION 5 but with one or more of the component calculations carried out using an approximation. For example, Torres-VerdIn & Habashy (1993) make use of an approximate forward modeffing algorithm, both for the forward modeffing within their inversion routine and to generate the Jacobian matrix of sensitivities, in their procedure for inverting cross- borehole measurements for a two-dimensional conductivity model. The generalised sub- space technique developed by Oldenburg, McGiffivray & Ellis (1993) enables the matrix equation at each iteration of a Gauss-Newton procedure to be inverted in much less time than the corresponding exact calculations. Mackie & Madden (1993) and Ellis & Old enburg (1994) use conjugate gradient techniques to solve the matrix equation at each iteration in their procedures for the multi-dimensional inversion of magnetotelluric data. These techniques have the additional advantage that the Jacobian matrix of sensitivi ties need not be explicitly generated, but rather that additional forward modellings are carried out to compute the product of the Jacobian matrix with the necessary vectors. Smith & Booker (1991) use a Jacobian matrix of approximate sensitivities in their inver sion procedure for the two-dimensional magnetotelluric problem. Until it is possible to carry out the necessary computations for multi-dimensional inverse problems in the same amount of time as is currently possible for the one-dimensional magnetotelluric prob lem, Gauss-Newton algorithms using accelerated, approximate ways of both generating the Jacobian matrix of sensitivities and solving the linear system of equations at each iteration offer the most practicable means of constructing solutions to multi-dimensional problems. This thesis is concerned with the development of a rapid, approximate method for generating the Jacobian matrix of sensitivities needed by any iterative, linearised in version procedure. In particular, an approximate form of sensitivity is sought that will be suitable for the multi-dimensional inversion of geophysical electromagnetic measure ments, and that will be sufficiently general to apply to any survey configuration. I begin in Chapter 2 by describing a Gauss-Newton procedure for inverting time-domain electro magnetic measurements to recover a one-dimensional conductivity model of the Earth. The purpose of this chapter is two-fold: (1) to illustrate the inversion philosophy and CHAPTER 1: INTRODUCTION 6 major components of a Gauss-Newton algorithm unrestricted by computing power, and (2) to develop an inversion program, albeit one-dimensional, for a common controlled- source survey configuration. The contents of Chapter 2 have been published (Farquhar son & Oldenburg 1993) and the inversion program is currently being used by a number of mineral exploration companies. Chapter 3 is an introduction to the subject of approximate sensitivities. In this chap ter, I discuss the need for a rapid method of approximating the sensitivities and describe the previous work on this subject. A statement of the approximate form of sensitivity that I present in this thesis is given in Chapter 4. The approximation is obtained from the adjoint-equation method (a derivation of which is included in Chapter 4) by replacing the true adjoint field with an approximate adjoint field that is considerably quicker to compute. In Chapter 5, I discuss the possible forms of this approximate adjoint field, and compare the main candidates with the true adjoint field for a purely two-dimensional example. The approximate sensitivities for the two-dimensional magnetotelluric problem obtained using the approximate adjoint fields described in Chapter 5 are investigated in Chapter 6. The approximate sensitivities are compared to accurate sensitivities for a simple conductivity model. A comparison of their respective computation times is also carried out. To conclude Chapter 6, the approximate sensitivities are tested in a Gauss- Newton inversion procedure analogous to that described in Chapter 2. Magnetotelluric data from a synthetic data-set and a field data-set (a sub-set of the COPROD2 data) are inverted to recover two-dimensional conductivity models. In Chapter 7, I compute approximate sensitivities for the 2.5d and three-dimensional controlled-source problems. Comparisons are made of the approximate and accurate sen sitivities for simple conductivity models. The comparisons for the three-dimensional problem are limited by computing power, but they are sufficient to illustrate the ma jor differences in the respective computation times for the approximate and accurate sensitivities. Finally, the major conclusions of this thesis are summarised in Chapter 8. Chapter 2 One-Dimensional Inversion of Time-Domain Electromagnetic Data 2.1 Introduction Gauss-Newton inversion procedures for geophysical electromagnetic data have un dergone much development in the context of the one-dimensional magnetotelluric prob lem. These procedures now concentrate on obtaining the model that has the minimum value of some norm subject to the constraint that the observations are reproduced to an acceptable level. Such a procedure has not yet been applied to the controlled-source sur vey configurations commonly used in the exploration industry, even for one-dimensional models. Fullagar & Oldenburg (1984) did construct an iterative, linearised procedure for the one-dimensional inversion of frequency-domain data collected using the horizon tal loop source-receiver configuration. However, Fullagar & Oldenburg minimised the norm of the model perturbation at each iteration rather than the norm of the model itself. In this chapter, a Gauss-Newton procedure is developed for inverting data from a controlled-source sounding. The procedure contains aid the features present in the mod ern approaches to solving the one-dimensional magnetotelluric inverse problem. As such, this chapter also serves to illustrate the ultimate capabilities of inversion procedures for multi-dimensional problems. The geophysical method for which the inversion procedure is developed in this chap ter is the time-domain electromagnetic (TEM) method. The TEM method is used exten sively in the exploration, geotechnical and environmental applications of geophysics. A review of the method and its uses is given by Nabighian & Macnae (1991). Typically, a step or ramp turn-off in the current flowing in a rectangular transmitter loop is used to induce currents in the Earth, and the variation with time of the vertical component of the 7 CHAPTER 2: 1 d INVERSION OF TEM DATA 8 magnetic field, or its time derivative, resulting from these induced currents is measured. These measurements can be at any point on the surface of the Earth, either inside or outside the transmitter ioop. The conductivity model for the inversion procedure is composed of a large number of horizontal layers of fixed thickness, and is terminated by a half-space. In general there are many more layers than observations, so the inverse problem is under-determined. This increases the nonuniqueness of the mathematical solution but allows the model that minimises a specific norm to be found from the infinity of models that adequately re produce the data. Suitable choice of the model norm enables the inversion procedure to produce the model that is the most plausible representation of the region of the Earth under investigation. Ideally, the model that is produced should exhibit the right “char acter” (that is smooth or blocky in accordance with the assumed geology), be as close as possible to the conductivity section obtained from a well-log or a neighbouring sounding (if such is available) and have a minimum amount of structure. This last point is partic ularly important since arbitrarily complicated structures, which would seem unlikely to resemble the Earth, can suffice as mathematical solutions to the inverse problem. The aim is to generate a model that contains just enough structure to fit the observations, but no more. The flexibility to generate models of a particular character, and hence produce the most plausible model for a particular geological setting, is a fundamental feature of the inversion algorithm developed in this chapter for the TEM method. 2.2 The forward problem Consider firstly the forward problem of calculating either the vertical component of the magnetic field, h(t), or its time derivative, 8h(t)/ot, induced at some point on the surface of a layered conductivity structure by a step turn-off in the current in a rectangular transmitter loop. To exploit the work that has been done on electromagnetic methods in the frequency domain, the analysis is carried out in the frequency domain and only at the very end are the results transformed to the time domain. The notation of CHAPTER 2: id INVERSION OF TEM DATA 9 Ward & Hohmann (1988) is used in which lower-case letters represent fields in the time domain and upper-case letters represent fields in the frequency domain. The magnetic field, H, due to a rectangular transmitter loop can be evaluated by in tegrating H due to a horizontal electric dipole around the transmitter ioop (Poddar 1982). The major portion of the forward problem therefore involves calculating H due to a hori zontal electric dipole. The computation of the magnetic and electric fields due to a dipole source on the surface of a homogeneous or layered halfspace is required throughout this thesis. The method used to compute these fields is described in Appendix A. The method given in Appendix A for computing the magnetic field in a layered halfspace due to a horizontal electric dipole source makes use of the Schelkunoff poten tials A = Ae and F = Fe (section 4, Ward & Hohmann 1988) where ê, is the unit vector in the z-direction. When computing just the vertical component of the magnetic field, only the potential F is required. In the th layer of the conductivity model (see Fig. A.1), F, the two-dimensional Fourier transform of F, satisfies the ordinary differen tial equation ( — = 0, (2.1) where u = k + k — ew2 + iwp.o-. and e are the magnetic permeability and electrical permittivity, respectively, of the layer: they are assumed fixed in the inversion process but can have different values in different layers. cr is the conductivity of the th layer. The tilde indicates the (k, ky,,) domain. The steps for solving this differential equation are given in Appendix A. Once F (kr, k, z,w) has been found for all j, the time-domain values of either h(t) or Oh(t)/Ot for the rectangular transmitter loop are obtained by the sequence of linear transformations shown in Fig. 2.1. The inverse Hankel transform of F that gives H as a function of space and frequency for the dipole is numerically evaluated using the digital filtering technique of Anderson (1979b). Integration of the dipole response around the rectangular transmitter ioop is done using a Romberg integration routine (e.g., Press et al. 1992). Finally, the time-domain values of h are computed from the frequency-domain CHAPTER 2: 1 d INVERSION OF TEM DATA 10 (km, k, z, w) dipole inverse Hankel transform H (x, y, z, Li)) di pole Integration around transmitter loop H (x, y, z, w) loop Inverse Fourier transform h (x, y, z, t) or 8h (x, y, z, t)/8t loop Figure 2.1 The sequence of integrations used to transform the values of the vector potential due to a horizontal electric dipole as a function of wavenum ber and frequency to values of the vertical component of the h field, or its time derivative, for the rectangular transmitter ioop as a function of space and time. This sequence of linear transformations is also used to obtain the sensi tivities 8h/9cr (z, y, z, t) for the rectangular transmitter loop from the values of OF/9o (k,k,z,w) for the dipole. values of H using the routine of Newman, Hohmann & Anderson (1986). This sequence of linear transformations is also used in the next section to obtain the sensitivity of h from the sensitivity of F. 2.3 Calculation of the sensitivities In order to construct an iterative, linearised algorithm for inverting values of h(t) or 8h(t)/8t for a horizontally layered conductivity model, a procedure is required for calculating the Jacobian matrix of sensitivities where the sensitivities are the partial CHAPTER 2: 1 d INVERSION OF TEM DATA 11 derivatives of the data with respect to the model parameters. For the problem in this chapter, = , (2.2) where h is the jth observation and o is the conductivity of the j’ layer. Here, and for the remainder of this chapter, h is used to represent either h or Oh/Ot, depending on which is being considered as the data. The method used here to calculate the sensitivities for the one-dimensional inverse problem is a special case of the general adjoint-equation method for calculating the sen sitivities for the multi-dimensional electromagnetic inverse problem. The general method will be described in Chapter 4, and will form the basis for the approximate form of sensitivities developed in this thesis. The layered conductivity structure that is the model for the one-dimensional inverse problem (see Fig. A.1) can be described in terms of a linear combination: (z) = (z), (2.3) where the basis function is equal to unity within the th layer and zero everywhere else. The coefficient is equal to the conductivity of the th layer. Using this definition for the layered conductivity structure, eq. (2.1) can be rewritten in a form that is valid for all zE(—oo,00): /2 N — u0 — iw crjb) F = 5, (2.4) where u = k + k — jt ew2 and S represents the dipole source in the appropriate layer. Differentiating with respect to 0k (McGifflvray & Oldenburg 1990), making use of the chain rule and realising that S is independent of o, eq. (2.4) becomes = iLL)/Tk.F, (2.5) CHAPTER 2: 1 d INVERSION OF TEM DATA 12 which is an inhomogeneous ordinary differential equation for the sensitivity OF/Ook. The boundary conditions are —* 0 as z — ±00. This boundary value problem can be solved using the adjoint Green’s function method (Lanczos 1961). Hence, = f iwk(z)F(z)Gt(z; dz, (2.6)- where the adjoint Green’s function Gt(z; c) satisfies N (_u+iwaj)Gt(z;) S(z-) (2.7) and G(z;) —* 0 as z —* +00. (2.8) To construct the adjoint Green’s function, two linearly independent solutions of the homogeneous form of eq. (2.7) are needed, one of which satisfies the boundary condition as z — —cc, and the other which satisfies the boundary condition as z —* +co. For the problem under discussion here, it is only the value of the sensitivity at the surface of the Earth (c = 0 in eq. 2.6) that is of interest. In the region —cc < z = 0, where the conductivity is zero, the adjoint Green’s function has the form exp(u0z). In the region = 0 z < cc, eq. (2.7) is the complex conjugate of eq. (2.4). And, since P —* 0 as z —* +oo, the adjoint Green’s function is proportional to F* in this region. Hence, Gt(z. — fce_, z_0, 29 “ ‘‘ — lc+F*(z), z<=0. At z = = 0, Gt(z; ) must be continuous, and its derivative with respect to z must be discontinuous by an amount equal to 1 (Roach 1982). This determines the two coefficients: *() c_ = I*(o) — uo*(0) 1 (2.10) c+ = - - FI*(O) — uoF*(0) (2.11) CHAPTER 2: 1 d INVERSION OF TEM DATA 13 where the prime denotes the derivative. So, using the explicit form of the adjoint Green’s function given above, and remem bering that is unity in the th layer and zero everywhere else, eq. (2.6) becomes Ui) = j [E(z)j2 dz. (2.12)F (0) —u0F(0) z—z This expression for the sensitivity links the change in F(k, k, z = 0,w) to the change in the conductivity of the j’ layer. F is known throughout the layers from the forward modeffing mentioned in the previous section. The desired sensitivity, Oh/Oo, can be obtained from 8F/Uu by the series of linear transformations shown in Fig. 2.1. 2.4 The inversion algorithm (obs) .Consider a set of M observations, h , i = 1, . . . , M, which can either be of the vertical component of the h field, or of its time derivative. The goal of the inversion algorithm is to find a set of layer conductivities that adequately reproduces these obser vations. Because the conductivities found in the Earth can vary by orders of magnitude, it is convenient to work with the logarithm of conductivity in the inverse problem. Also, working with logarithms ensures that o is positive. Let m = ln o, and let the vector m = (m1,. . . , m)” define the model for the inverse problem. The inverse problem is nonunique: if there is one model that adequately reproduces the finite set of observations, then there is an infinite number of such models. To find a specific model, a norm of the model is minimised which has the form 4’m = II Lm (m — m) j2 (2.13) where Wm is a weighting matrix, is a reference model and . represents the 12 norm. The model, m, that minimises 4m will have a character that depends on the particular choices for the weighting matrix and the reference model. The reference model CHAPTER 2: id INVERSION OF TEM DATA 14 can be used to include any a priori information that may be available about the possible conductivity structure. The appropriate numerical values of Wm for generating a model of a specific char acter can be obtained by considering a functional analogous to m for models that are continuous functions of depth, for example, = jw(z) T[m(z) _m(z)1I2 dz. (2.14) The operator T can be the identity operator, or the first- or second-order derivative with respect to z. The function w(z) is an additional weighting function which can be used to enhance or suppress structure over certain depth ranges. The weighting matrix Wm can be obtained by making m the discrete equivalent of m To determine whether or not the data produced by the model conductivity structure are sufficiently close to the observations, the following measure of misfit is used: = (h(0b — h[m]) 2 (2.15) where h[m] represents the data computed for the model m. W is a weighting matrix which is usually a diagonal matrix whose elements are the reciprocals of the error estimate of each datum. The objective in the inversion is to find a model which gives a misfit, d’ equal to a target misfit For the examples considered in this chapter (and throughout this thesis), it will be assumed that the measurement errors of the data are unbiased, independent and Gaussian. In such cases, q is equal to the x2 random variable. From the properties of the x2 random variable (Menke 1984; Parker 1994), the expected value of cbd is equal to M, the number of observations. Hence, the final target misfit in the inverse problem is usually taken as = M. The relationship between the observations and the model parameters is nonlinear and so an iterative method is required to solve the inverse problem. The data misfit after CHAPTER 2: id INVERSION OF TEM DATA 15 the (n + l)th iteration is given by (n+i) = 1 (h(0bs) — h[m’)]) 2 (2.16) To avoid excessive structure developing in the model during early iterations, the target th • (n) (tar)misfit for the (n + 1) iteration is chosen to be max (,Bq ) where typically 0.1 < 3 < 0.5. Also, it is the norm of the model produced by this iteration that is to be minimised: (n+i) = II Wm (m(’ — m) 2 (2.17) The inverse problem to be solved at the (n + l)th iteration is therefore: Minimise (2.18) subject to (n+1) max(,t) c, (2.19) The general approach to solving this optimisation problem is to construct an objec tive function (n+1) + ((fl+1) - q), (2.20) where -y is a Lagrange multiplier. This objective function is then linearised about the model using a Taylor series expansion, differentiated with respect to the model perturbation Sm = m(1) — and the resulting derivatives equated to zero. The set of simultaneous equations is then solved for the model perturbation. (See, for ex ample, Gill, Murray & Wright 1981; Oldenburg, McGiffivray & Ellis 1993.) It is this method of solving the linear inverse problem at each iteration that will be used for multi dimensional electromagnetic inverse problems elsewhere in this thesis. However, for the problem considered in this chapter, the method of singular value decomposition (SVD) was used. This method is particularly efficient when the matrix equation to be solved at each iteration is not too large (less than, say, 200x200), which is the case for the one-dimensional problem treated here. CHAPTER 2: 1 d INVERSION OF TEM DATA 16 Consider again the data misfit at the (ri+ i)tI iteration given in eq. (2.16). Consider a Taylor series expansion of h[m(’)] about the model h[m2419 = h[m] + J[m] Sm + R, (2.21) where = Oh/8m = uOh/ao (since the model is in terms of the logarithm of the layer conductivities) is the Jacobian matrix of sensitivities. Assuming the remainder term R can be neglected, substituting eq. (2.21) into eq. (2.16) gives a linearised estimate, çj’, of the true data misfit: = W4 (h(0b — h[m] — [Sm) 2 (n+1) (2.22) By writing Sm explicitly as the difference between the model parameters at two successive iterations (e.g., Oldenburg 1983; Constable, Parker & Constable 1987), and introducing the reference model, eq. (2.22) becomes = II W (h(0) — h[m] — 3 + — .j m’’ + 3 m) 2 d — (2.23) where d = W (h(0b — h[m] — Jm + Jm()) (2.24) = JW’, (2.25) — m9. (2.26) The linearised inverse problem to be solved at the (m + i)tI iteration can now be stated as: Minimise cbm = Ilth4’W2 (2.27) subject to q = — jth(n+1)II2 = (2.28) CHAPTER 2: id INVERSION OF TEM DATA 17 The SVD solution of an under-determined system of equations is the one with the smallest 12 norm (Wiggins 1972; Parker 1977; Menke 1984; Golub & Van Loan 1989). Writing the SVD of as J = UAVT the solution to the linear inverse problem in eqs. (2.27) and (2.28) is given by = VTA_1UT (2.29) where = s S/(-y s + 1), s is the th singular value of and is the Kronecker delta (Wiggins 1972). The Lagrange multiplier y (the same as that in eq. 2.20) could be chosen using a line search so that the constraint q q in eq. (2.28) is satisfied. However, it is a solution to the full nonlinear inverse problem that is required. A line (n+1) *search is therefore used to choose -y such that bd = cbd. The value of -y from the previous iteration is used as an initial estimate for y for the current iteration. The corresponding model is computed using eq. (2.29), and a forward modeffing carried out to determine the data misfit, bd, for this model. 7 is then halved, the corresponding model computed using eq. (2.29) and a forward modelling carried out to calculate the data misfit for this new model. 7 is then repeatedly altered, assuming a locally linear relationship between in cd and in , until either q5, or the smallest possible data misfit, is obtained. 2.5 Examples 2.5.1 Synthetic data Synthetic data were generated from the three-layer conductivity structure shown in Fig. 2.2. The transmitter loop was a square of side 50 m and the vertical component of the h field due to a step turn-off in a 1 A current was calculated 50 m from the centre of the loop. Both the transmitter and receiver were on the surface of the conductivity structure, and 20 values of h were calculated over the range of delay times shown in Fig. 2.3. Gaussian random noise with a standard deviation of 2.5% was added to the CHAPTER 2: id INVERSION OF TEM DATA 18 10—1 S cf b 102 Figure 2.2 The three-layer model (“true”) used to generate the synthetic data, and the three models produced by the inversion routine. For clarity, the individual layers in the smallest, flattest and smoothest models are not shown. values of h. The data, with error bars equal to the standard deviation of the added noise, are shown in Fig. 2.3. The inversion routine was used to obtain three different conductivity models, each with a distinct character, that reproduced the synthetic data to the same level of mis fit. These three models minimise, in turn, the difference between the model and a reference halfspace of 5 x 10—2 S m1, the gradient of the model, and the curvature. For convenience, these models are referred to respectively as the smallest, flattest and smoothest. The controlling factor for each model is the weighting matrix YLm To ob tain the appropriate form of Wm for the smallest model consider eq. (2.14) with w = 1 and T = I, where I is the identity operator. Consider also a description of the layered model analagous to that in eq. (2.3). Eq. (2.14) then becomes = J [m — m] (z)2dz (2.30)z0 j=1 100 101 102 Depth (m) S N I I I I I CHAPTER 2: id INVERSION OF TEM DATA 19 1 O 1 o io io Time (s) Figure 2.3 The synthetic data, and associated error bars, generated from the three-layer model in Fig. 2.2. The values of h were calculated 50 m from the centre of a 50 m x 50 m transmitter ioop. A step turn-off in a 1 A current was used as the transmitter current waveform and Gaussian random noise of standard deviation 2.5% was incorporated into the data. The continuous curves represent the data predicted by each of the three models shown in Fig. 2.2. N N J [m_m] [mk—mj(z)k(z)dz (2.31)z=O N = (m — m )2 f dz (2.32)j=1 z=zi_1 = (m — )2 , (2.33) because of the particular nature of the basis functions, t is the thickness of the th layer. Comparing eq. (2.33) with eq. (2.13), the weighting matrix required for the smallest model is therefore = diag (‘T (2.34) CHAPTER 2: 1 d INVERSION OF TEM DATA 20 The thickness of the basement halfspace is infinite, and so the final element of Wm is assigned the same value as the previous element. The value of the model in the basement halfspace will not therefore dominate the model norm. To obtain the weighting matrix for the flattest model consider eq. (2.14) with w = 1 and T = d/dz: d 2 = J -- [m(z) — m(z)] dz. (2.35)z=O uZ An analogous form of this expression for the discrete case that results in an appropriate model norm to minimise is N—i / (ref) (ref) 2i[m. -m ]-[m.--m. 1’ = 3 Lz. ‘ ) (2.36)j=i 3 Letting Lz equal the vertical distance between the centre of the th layer and the centre of the layer above, comparison of eqs. (2.36) and (2.13) implies that the weighting matrix for the flattest model is given by — _____ 1/2 (tl±t2) 1/2 0 . . . 0 0 —(t+t3)/2 (t2+t3)_h/2 0 0 _(tN_i)uh!2 (tN1)h/2 0 ... 0 -C (2.37) where the constant C was chosen as 103(tNi)/2. This value is effectively negligible compared with the other elements in the matrix, and yet is sufficiently large so that W’ can be computed. For the smoothest model, the weighting matrix is the square of the matrix used for the flattest model. The resulting three models are shown in Fig. 2.2. All have 100 layers (i.e., N=100), the thicknesses of which increase exponentially (to the base 1.05 with the thickness of the first layer equal to 0.5 m). During the inversion process, 3 (the desired reduction in misfit at each iteration) was kept fixed at 0.5. This gave a slow but steady convergence towards the final model. Fig. 2.4 shows the variation of q and m during the course CHAPTER 2: 1 d INVERSION OF TEM DATA 21 10 100 10-1 10-2 / b io—3 I 0 2 4 6 8 10 12 Itcration Figure 2.4 The variation of (a) the data misfit, c1d (b) the model norm, m’ and (c) the Lagrange multiplier, 7, during the inversion of the synthetic data in Fig. 2.3 to produce the flattest model shown in Fig. 2.2. of the inversion to produce the flattest model. The thirteen iterations shown in Fig. 2.4 took 90 minutes on a Sun SparclO workstation. The data calculated from the three models produced by the inversion routine are represented by the continuous lines in Fig. 2.3. The values of c5d for the smallest, flattest and smoothest models were all equal to 20 (M = 20). Although these three models reproduce the observations to the same level of misfit, they are noticeably different in character. These different characters are a direct consequence of the particular model norm that was minimised in the inversion. The most obvious differences appear in the depth ranges that are poorly constrained by the data: below approximately 200 m the smallest model returns to its reference halfspace of 5 x 10—2 5m’, the flattest model levels out to achieve zero gradient, and the smoothest model degenerates to a straight line when no longer influenced by the data. Similar behaviour can be seen at shallow depths above about 10 m. Even in the depth range to which the data are most sensitive CHAPTER 2: id INVERSION OF TEM DATA 22 (10— 200 m, approximately) there are differences in character between the three models. The smallest model manages to follow the block nature of the true model quite closely, whereas the smoothest model smears out the high conductivity layer as much as it can in order to produce the model with the minimum amount of curvature. From Fig. 2.2 it is obvious that the smallest model is the closest to the “true” Earth. However, the remarkable agreement above 10 m depth is somewhat fortuitous, since the model at these depths is very poorly constrained by the data. And the exact agreement below 200 m is only to be expected since the reference model was chosen to have the same conductivity as the basement halfspace in the true model. For many of the other data-sets that the inversion routine was tested upon, it was found that the smallest model often contained an unreasonable amount of structure. A model with the least amount of variability required to reproduce the observations is generally preferred. 2.5.2 Field data To test the inversion routine with more realistic data, two sets of field data acquired during an environmental study were inverted. Values of the time derivative of the vertical component of the h field were measured at the centre of a square (60m x 60m) transmitter loop using the Geonics Protem system (see Nabighian & Macnae 1991). The data, and their assumed error estimates, obtained from soundings at two different locations are presented in Figs. 2.5 and 2.7. The Protem instrument uses a linear ramp turn-off of length r instead of a pure step turn-off, which is impossible to generate in practice. This modification of the source current waveform can be taken into account when calculating the resultant time decay of the fields by convolving the time-derivative of the linear ramp turn-off with the values of Oh/Ot calculated for a pure step-off current source (eq. 1, Asten 1987): ahr 00 Oh(t) = j (u) r’(t - u) du. (2.38) (1 r() = —t/r t < —T, —T t 0, t>0. Its time derivative, ‘1’ is therefore a boxcar of length T and height 1/T. The response, Oh”/Ot, for the ramp turn-off was obtained from the response, Oh/Ot, for the step turn-off by numerically evaluating eq. (2.38) using a Romberg integration routine (e.g., Press et Model fi: 191. Model f2: ç = 179. CHAPTER 2: id INVERSION OF TEM DATA 23 1 3 1 2 10 10—1 > 1 —3 1 —4 101 Time (s) Figure 2.5 A TEM sounding acquired during an environmental survey. The transmitter ioop was 60 m x 60 m, and the receiver was at its centre. The data were acquired in three overlapping sweeps: 7 its —0.7 ms, 0.1 —2.8 ms and 0.8— 70 ms. The observations and assigned error estimates are represented by the error bars. The continuous curve indicates the data predicted from both models fi and f2 in Fig. 2.6. The linear ramp turn-off is represented by (2.39) i05 I 1111111 I 1111111 I I I II III 1 al. 1992). CHAPTER 2: id INVERSION OF TEM DATA 24 The observations shown in both Figs. 2.5 and 2.7 were acquired in three overlapping sweeps, the first sweep from 7 s to 0.7 ms, the second from 0.1 to 2.8 ms, and the third from 0.8 ms onwards. The length of the ramp, r, for these sweeps was 3.5 jis, 35 ,us and 45 ,us respectively. These different ramp times give the data a somewhat sawtooth appearance in the time-range for which the sweeps overlap. The linear ramp turn-off in the transmitter current must also be taken into account in the inversion process. Because the inversion is being carried out in the time domain, the sensitivities required for a ramp turn-off in the transmitter current can be obtained by applying the convolution in eq. (2.38) to the sensitivities calculated for the pure step turn-off. Estimates of the errors were not available for the data-sets shown in Figs. 2.5 and 2.7. Values that seemed reasonable were therefore assigned to the data. For the data in Fig. 2.5, errors of 1% were assigned to the measurements before 50 ,us, and errors of 5% were assigned to all other measurements, except over the interval 90 ,us to 3 ms. In the centre of this interval (0.3 ms to 0.8 ms) there is an obvious discrepancy between measurements from the first and second sweeps, suggesting larger errors in these data. Errors of 10% were therefore assigned to the whole interval (90 ,us to 3 ms) over which there was overlapping of the three sweeps, although the discrepancy between the first and second sweeps in the first part of the interval was subsequently found to be due to the different ramp times for the two sweeps. For the data in Fig. 2.7, errors of 1% were assigned to the measurements before 80is, errors of 5% to the measurements between 80is and 5 ms, and errors of 10% to those after 5 ms. Fig. 2.6 shows two inversion results for the data in Fig. 2.5. Model fl, represented by the solid line, was obtained using the same weighting matrix as for the flattest model in Section 5.2.1. Model f2, represented by the dotted line, was obtained using the weighting CHAPTER 2: id INVERSION OF TEM DATA 25 matrix —1 1 0 ... 0 0—11 0 = ... (2.40) —1 1 0 0-C where C = iO. Both models have 200 layers increasing exponentially in thickness (to the base 1.03, with the first layer of thickness 0.1 m). Because of this increasing layer thickness, the model norm constructed using Wm above is essentially a discretised version of = x: [:]2donz = jz [d(lnu)]2 (2.41) (Remember that m = ln u.) The model norm constructed using Wm in eq. (2.31) is a discretised version of = f° [d(ln)]2 (2.42) The dashed line in Fig. 2.6 corresponds to the result obtained from a parameter estimation program (TEMEX—GL), based on Anderson’s (1979a) routine, which was restricted to contain seven layers. The values of the misfit, q, for models fi and f2 were 191 and 179, respectively. The value of d for the seven-layer model was 983. Models fl and f2 produce nearly identical fits to the observations, as illustrated in Fig. 2.5. However, the two values of misfit quoted above are still significantly larger that the expected value of 56 (the number of observations). No model could be found which gave a smaller misfit than that for model f2. This suggests that the error estimates that were assigned to the data are smaller than the true errors in the measurements. This seems entirely plausible, especially for the late-time measurements in each data-set, where errors of 30 or 40% are more likely to be realistic. The similarity between the two models in Fig. 2.6, despite the different nature of the weighting matrices used to produce them, emphasises the robustness of an inversion Figure 2.6 fi and f2 are the two versions of the flattest model obtained from the inversion of the data in Fig. 2.5. The best-fitting seven-layer model produced by a parameter estimation routine (TEMEX—GL) is also shown. routine which looks for a minimum-structure model. Since models fi and f2 both contain just enough structure to fit the data, it is not surprising that they agree on the features that are required by the data. This provides confidence that the particular sequence of conductive and resistive layers to a depth of about 300 m is present in the real Earth. Below 300 m the models are less strongly constrained by the observations, and so the particular form of the weighting matrix begins to dominates their behaviour. The weighting term z that appears in eq. (2.41) but not in eq. (2.42) causes structure in model f2 to be suppressed at these depths compared with model fi. However, the fact that both models fi and f2 do show increases in conductivity around 500 m depth, rather than levelling off, might indicate the presence of a good conductor on the extreme limit of penetration of this sounding. When this increase in conductivity was removed from CHAPTER 2: id INVERSION OF TEM DATA 26 Solid — model fi, Dotted — model f2, Dashed — parameter estimation. 101 1 O 10—1 1 —2 b 101 100 liii I I 111111 I 1111111 I [____I__I_IIIIII 101 1 2 Depth (m) 1 3 CHAPTER 2: id INVERSION OF TEM DATA 27 model if, the misfit did increase to a value of 219, and the predicted data from this modified model did not fit the last two data points in Fig. 2.5 as well as those for model fi. Similarly, removing the final increase in conductivity from model f2 resulted in an increase in the misfit to 194, suggesting that this feature might be required by the data. At shallow depths, the different weighting matrices lead to structure being enhanced in model f2 relative to model fi. As a final example, the data from the second sounding (see Fig. 2.7) were inverted. The weighting matrix was the same as the one used to produce model f2 in Fig. 2.6. The resulting model is represented by the dotted line in Fig. 2.8. It is again made up of 200 layers whose thicknesses increase exponentially with depth. The solid line in Fig. 2.7 represents the data predicted from this model. The corresponding value of the misfit was cd = 1216. The assigned error estimates again seem to be too small and it was not possible to obtain a misfit close to the expected value of 40. The solid line in Fig. 2.8 represents well-log measurements of conductivity obtained from a borehole 70 m away from the location of the sounding. There is good agreement between the model and the well-log measurements in the depth range 10— 150 m, especially on the location of the upper regions of the two conductive zones around 35 and 90 m. The character of the well-log, comprising transition zones rather than a few layers of constant conductivity, is also reproduced by the model in Fig. 2.8. This type of model, constructed by minimising the 12 norm of its gradient, therefore seems particularly suited to this geological setting. 2.6 Conclusions In this chapter, I have presented an inversion algorithm that generates a layered conductivity structure from measurements, at any point on the surface of the Earth, of the time-variation of h, or Oh/8t, induced by a rectangular transmitter loop. The inversion algorithm is a Gauss-Newton algorithm that minimises a norm of the model subject to the constraint that the observations are reproduced to an appropriate level. The form of the model norm to be minimised is flexible, and can be constructed to enable the algorithm to find the most plausible model from the infinity of models that Model: ‘kd = 1216. I 1111111 I 1111111 I I io Time I CHAPTER 2: 1 d INVERSION OF TEM DATA 28 1 3 1 2 101 100 10—1 1 —2 1 3 10-2 (s) Figure 2.7 A second sounding from an environmental survey. The survey geometry and the three measurement sweeps are the same as for Fig. 2.5. The error bars indicate the observations and assigned error estimates, and the con tinuous curve represents the predicted data from the flattest model produced by the inversion and shown in Fig. 2.8. adequately reproduce the observations. Which model is the most plausible will depend on the geological setting and the amount of prior knowledge, and may be the one with the least amount of structure, or the one that is the closest to some reference model that represents the preconceived image of the region under investigation. The algorithm described in this chapter is robust: it does not become unstable in the presence of noisy data. In addition, the combination of (1) finding a model that has a minimum amount of structure, and (2) requiring that the desired misfit at each iteration is only gradually decreased, leads to a procedure that consistently finds a model Figure 2.8 The flattest model (dotted line) obtained from the inversion of the data in Fig. 2.7. The solid line represents the values of the conductivity measured in a borehole 70 m from the observation location. CHAPTER 2: id INVERSION OF TEM DATA 29 Dotted — model, Solid — well—log. 10° 10_i 1 —2 (I) b 10—1 100 NH I 1111111 I 11111111 I I I I I III 101 Depth (m) 1 2 that reproduces the observations and is a plausible representation of the Earth. The final model produced by the algorithm is also insensitive to the particular starting model used. Finally, the solution to the one-dimensional controlled-source inverse problem de scribed in this chapter is not restricted by computing power unlike the multi-dimensional problems considered later in this thesis. All the calculations are therefore exact, including the generation of the Jacobian matrix of sensitivities. Chapter 3 Introduction to Approximate Sensitivities 3.1 Multi-dimensional problems Before I start to discuss approximate sensitivities, I shall summarise some concepts relevent to multi-dimensional geophysical electromagnetic inverse problems. The source for the magnetotelluric method is assumed to be a plane wave impinging on the surface of the Earth. Because of this simple source, the magnetotelluric method is the one for which the forward and inverse problems can be most easily carried out. The measured quantities in a magnetotelluric experiment are the horizontal components of the electric and magnetic fields. Since the details of the plane-wave source are unknown, the ratio of orthogonal components of the electric and magnetic fields is calculated (in the frequency domain). An apparent resistivity and phase are obtained from this ratio (see eq. D.1), and it is these quantities that are considered as the data acquired in the experiment. To investigate the two- or three-dimensional structure of the Earth, these data have to be collected at many different locations on the surface of the Earth as well as at many frequencies. For a two-dimensional conductivity model, the electric and magnetic fields induced in the model by a plane-wave source divide into two decoupled modes. One mode involves the component of the electric field parallel to the strike direction (i.e., the direction in which the model is invariant) and the two components of the magnetic field perpendicular to the strike direction. This mode is known as the “E-polarisation” mode since the only component of the electric field within the model is parallel to the strike direction. It is also referred to as the “transverse electric” (TE) mode since the only component of the electric field in the model is perpendicular (or transverse) to the vertical direction (see section 4.1.1, Weaver 1994). The other mode involves the component of the magnetic field 30 CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 31 parallel to the strike direction and the components of the electric field perpendicular to the strike direction. It is therefore called the “H-polarisation” or “transverse magnetic” (TM) mode since the only component of the magnetic field in the model is parallel to the strike direction (and transverse to the vertical direction). This decoupling into two modes also applies to any situation in which the source is invariant in the strike direction of the model. In the exploration industry, controlled-source methods are used more often than the magnetotelluric method. A typical controlled-source method has either a current loop or a grounded line current as a source, and involves measurements of the frequency- or time-dependence of the electric or magnetic fields. When considering a two- or three- dimensional target region, these measurements are either collected at many different locations relative to a stationary source, or at a fixed distance from the source for many different source positions, or a combination of the two. In this thesis, the term “forward modeffing” is used to refer to the process of com puting the values of the measurements that would have been obtained if the geophysical experiment had been carried out over the mathematical conductivity model rather than over the Earth itself. One forward modeffing produces values for the measured quantities at every observation location for every source location. I shall further specify that one forward modeffing, as considered in this thesis, corresponds to only one frequency or de lay time: if measurements were collected at ten frequencies then ten forward mode]iings would be required to compute values for all the quantities measured in the experiment. I shall also assume that, in the process of computing the predicted data for a particular frequency or delay time, a forward modeffing produces the values of the electric field everywhere in the conductivity model. As a final point, it is worth noting the correspondence between “sensitivities” and “Fréchet derivatives”. Sensitivities are obtained by first parameterising the conductivity model and then differentiating the equations that define the forward problem with respect to the model parameters (see sections 2.3 and 4.2). An alternative approach is to consider CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 32 a general perturbation to the model before any parameterisation is introduced. The resulting effect on the observed quantity can be expressed as (Parker 1994): d [m(r) + Sm(r)] = d {m(r)] + (D(r),Sm(r)) + R, (3.1) where d represents the observed quantity and (.,.) represents the inner product. If the remainder R is such that jR/j6mj — 0 as — 0, then D is the Fréchet derivative of d. If the model parameterisation is now introduced such that m(r) = m&(r) and Sm(r) = Smib(r), eq. (3.1) becomes d [m+ Sm] = d [m] + (D(r),(r)) 5m + R, (3.2) where m = (m1,... , m)’. For the problems considered in this thesis, the inner prod uct (•,) is the volume integral over the whole of three-dimensional space. Hence, (D(r),(r)) = jDi(r)j(r)dv. (3.3) Comparing eq. (3.2) with eq. (1.1), it can be seen that the sensitivity is equivalent to the volume integral of the product of the Fréchet derivative and the basis function: I = I D(r) (r) dv. (3.4) um j 3.2 The need for approximate sensitivities The Gauss-Newton algorithm described in Chapter 2 is the inversion methodology that one would like to use to solve the two- and three-dimensional electromagnetic in verse problems. However, the computations required by some of the components of this algorithm become very time-consuming when applied to two- and especially three dimensional problems. The three most computationally intensive components are (1) the forward modeffing, (2) the generation of the sensitivity matrix and (3) the solution of CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 33 the linear system of equations. Each of these operations has to be performed at least once for each iteration. Highly sophisticated and efficient forward modeffing schemes are now being developed for two- and three-dimensional conductivity models (e.g., Everett & Edwards 1992; Unsworth, Travis & Chave 1993; Mackie, Madden & Wannamaker 1993; Livelybrooks 1993; MDona1d & Agarwal 1994; Newman & Alumbaugh 1995). The in version of the linear system of equations has also received attention, and techniques such as a generalised subspace method can be used to speed up the solution of the matrix equation without hindering the progress of the iterative inversion scheme (Oldenburg, WGiffivray & Ellis 1993; Oldenburg & Li 1994). The generation of the Jacobian matrix of sensitivities, however, has received little attention. The three main methods of calculating the sensitivities for the general electromag netic inverse problem are the brute-force or perturbation method, the sensitivity-equation method and the adjoint-equation method (McGiffivray & Oldenburg 1990). The compu tation times for these methods are roughly equivalent to N x Mf, N x Mf and M0 x Mf forward modeffings respectively where N is the number of model parameters, Mf is the number of frequencies (or delay times) and M0 is the number of observation locations (Mf x M0 = M). The actual time taken by each of these three methods will depend on their implementation for a particular problem. It is sometimes possible to exploit features of a forward modeffing algorithm to reduce the time needed to generate the sensitivities. For example, Oristaglio & Worthington (1980) make use of the already-factored matrix from their finite-difference forward-modeffing program when calculating the sensitivities using the sensitivity-equation method. However, the above estimates of the computation times for the three methods of calculating the sensitivities are reasonable when construct ing an inversion scheme around an existing forward-modeffing program, and are a good indication of how the computation times for the sensitivities are affected by the number of model parameters and data. For surveys over two- or three-dimensional geological structures, the typical number of measurements is of the order of one or ten thousand. The number of parameters re quired to adequately describe any model used in the inversion of such measurements will CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 34 also be of the order of one, ten or a hundred thousand. Given that forward-modelling programs for such two- or three-dimensional conductivity structures currently have com putation times of at least several minutes, it can be seen that calculating the sensitivities using any of the accurate methods described above will be very time-consuming, if not prohibitive. It therefore becomes necessary to consider approximate methods for gen erating the Jacobian matrix of sensitivities if multi-dimensional electromagnetic inverse problems are to be tackled. 3.3 Previous forms of approximate sensitivities 3.3.1 The one-dimensional problem There has been one investigation of approximate sensitivities for the one-dimensional electromagnetic inverse problem. Boerner & Holladay (1990) demonstrated that the sensitivities with respect to the layer conductivities do not depend strongly on the model when components of the electric or magnetic fields associated with only horizontally- flowing currents in a horizontally-layered model are considered. They illustrated this by considering a horizontal electric dipole source and measurements of the resulting vertical component of the magnetic field, with both the source and measurement location on the surface of the Earth. They computed the sensitivities of these measurements with respect to the conductivity of a thin test layer for three models: a homogeneous halfspace of conductivity 0.01 S m1, a halfspace of 0.01 S m1 in which a conductive layer of 0.05 S m1 was embedded, and another halfspace of 0.01 S m1 in which a resistive layer of 0.002 S m1 was embedded. The plots of the sensitivities as functions of the depth of the test layer illustrate the weak dependence of the sensitivities on the model: all three curves are smooth and continuous across the boundaries of the conductive or resistive layer, they all exhibit the same general characteristics, and differ at most by 10%. Boerner & Holladay call the sensitivity with respect to the conductivity the “incre mental” sensitivity. The sensitivity with respect to the logarithm of the conductivity, CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 35 what Boerner & Holladay call the “fractional” sensitivity, is obtained by scaling the in cremental sensitivity by the conductivity of the layer: c9/O(ln u) = o O/8o. This is the sensitivity that is applicable when the model is defined in terms of the logarithm of the conductivity, as is usually the case. Because of the weak dependence of the incremental sensitivities on the structure in the model, the dominant influence on the fractional sen sitivities is the scaling by the layer conductivities. Boerner & Holladay exploited this by inverting synthetic data using a simple, linearised procedure in which the fractional sensi tivities were approximated by the incremental sensitivities calculated for a homogeneous halfspace and scaled by the conductivity of the layers. They found that the convergence of the inversion procedure was not hindered by the use of these approximate sensitivities, and that in some instances the rate of convergence was increased. Boerner & Holladay state that the weak dependence of the incremental sensitivities on the conductivity model is not to be expected for measurements of the electric and magnetic fields associated with current flow across the boundaries of the layers. However, the plots of Boerner & West (1989) show that, for a situation in which currents flow across the layer boundaries, the fractional sensitivities are also strongly dependent on the scaling by the layer conductivities, and that the fractional sensitivities for a one-dimensional model can be approximated by scaling the incremental sensitivities for a homogeneous halfspace by the conductivities of the layers. (The situation Boerner & West considered was that of a grounded electric dipole source and measurements of the corresponding electric field, although calculated in the wavenumber domain, at zero frequency and with a geometrical factor removed.) 3.3.2 The two-dimensional problem Approximate sensitivities have been used to reduce the time required to invert mag netotelluric data for two-dimensional conductivity models. Oldenburg & Ellis (1993) and Effis, Farquharson & Oldenburg (1993) used the sensitivities from the one-dimensional magnetotelluric inverse problem as approximate sensitivities for the two-dimensional CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 36 problem. For each observation location, a horizontally-layered conductivity structure is imagined in which the conductivities of the layers are equal to the conductivities of the cells in the two-dimensional conductivity model immediately beneath the observation location. One-dimensional sensitivities are calculated with respect to the layer conductiv ities. These sensitivities are then used in the two-dimensional inversion as the sensitivities for the cells immediately below the observation location. It was found that iterative inversion procedures using these approximate sensitivities converged to acceptable models when the data were either for the E- or il-polarisation, or determinant averages of the data for the two modes. These approximate sensitivities were not as successful, however, when used to jointly invert observations from both the E- and il-polarisations, especially when there were significant differences between the data from the two modes. This is not surprising since the decoupling into two distinct modes does not exist in the one-dimensional magnetotelluric problem, and so the one- dimensional sensitivities used as approximate sensitivities in the joint inversion were the same for both modes. Smith & Booker (1991) developed a slightly better form of approximate sensitivities for the two-dimensional magnetotelluric problem. They assumed that the horizontal gradient of the electric field in the two-dimensional conductivity model was small relative to its vertical gradient. The resulting expressions for the approximate sensitivities are the same as the expressions for the one-dimensional sensitivities used by Oldenburg & Effis and Ellis et al., but the electric field in the two-dimensional conductivity model is used rather than the electric field in the layered model imagined to exist beneath each observation location. Smith & Booker’s approximate sensitivities are therefore different for the E- and H-polarisations, and through the electric field from the two-dimensional model, a small amount of information about the two-dimensional structure of the model is incorporated in the sensitivities. The advantages of the above two forms of approximate sensitivity are that they are quick to compute and are sufficiently accurate to enable, in many cases, iterative, linearised inversion schemes to converge to an acceptable model. The “Rapid Relaxation CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 37 Inversion” program of Smith & Booker (1991) is routinely used to invert data from regional magnetotelluric surveys. The disadvantage of the above forms of approximate sensitivity is that they lead to sparse Jacobian matrices: the approximate sensitivity for a measurement is only defined with respect to the conductivities of the cells immediately below the location of the mea surement. There is no information, therefore, in these approximate sensitivities about the dependence of a measurement on the cells that are close to the location of measurement but not immediately beneath it. Another consequence is that the discretisation of the model is restricted to have each column of cells below an observation location. 3.3.3 The three-dimensional problem In principle, the approximate sensitivities of Smith & Booker (1991) described in section 3.3.2 can be extended to the three-dimensional magnetotelluric problem. Also, attempts have been made to calculate approximate sensitivities using the Born approxi mation (see Hohmann & Raiche 1988). However, the amount of work that has been done on approximate sensitivities for three-dimensional problems is small. 3.4 Conclusions Sections 3.3.1 and 3.3.2 suggest that approximate sensitivities need not be extremely accurate to enable a linearised, iterative inversion procedure to converge to an acceptable model. The work of Boerner & Holladay (1990) shows that the fractional sensitivities are heavily dependent on the scaling by the conductivity structure of the model, rather than the specific features of the incremental sensitivities. A moderate approximation of the incremental sensitivities, when scaled by the conductivities of the model, should therefore lead to a good approximation of the fractional sensitivities. The relative success of the approximate sensitivities described in section 3.3.2, despite their somewhat poor level of approximation, strengthens the argument that approximate sensitivities need not be highly accurate to be useful. CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 38 It is also worth noting that in an iterative Gauss-Newton procedure for solving the non-linear inverse problem (see eq. 2.20 and the surrounding description), the Jacobian matrix of sensitivities is evaluated at the current model: J[m(’)]. Consider for the mo ment that the need for an iterative procedure has not yet been recognised. To find the model that minimises the objective function, therefore, the objective function would sim ply be differentiated with respect to the model parameters and the resulting expressions equated to zero. This would produce a matrix equation for the model parameters, m(s01), that minimise the objective function. However, this matrix equation involves the Ja cobian matrix evaluated for this unknown model: J[m(s01)]. The matrix equation is therefore non-linear and intractable. (See the section “The Nonlinear Problem” of Con stable, Parker & Constable 1987). Comparing this to the iterative approach, [[m(’)] can be thought of as an approximation to J[m(s01)], with only the dominant features of J[m(’)] influencing the direction in which the iterative procedure moves through model space, at least for the early iterations. This further suggests that approximate sensitivi ties might be effective in enabling an iterative, linearised inversion algorithm to converge to an acceptable solution. It is the aim of this thesis to develop a general form of approximate sensitivity that will be of use in the solution of any electromagnetic inverse problem, especially those involving controlled sources. Chapter 4 An Approximate Form of Sensitivity for the General Electromagnetic Inverse Problem 4.1 Introduction The three main methods of accurately calculating the sensitivities are the brute- force or perturbation method, the sensitivity-equation method and the adjoint-equation method (McGillivray & Oldenburg 1990). The computation times of these methods are roughly equal to N x Mf, N x Mf and M0 x Mf forward modellings respectively. (N is the number of model parameters, Mf is the number of frequencies and M0 is the number of observation locations). The inversion procedure envisioned for the multi-dimensional problem, just as for the procedure described in Chapter 2 for the one-dimensional prob lem, is one in which there are many more more model parameters than observations. The adjoint-equation method is therefore the most efficient method of calculating the sensitivities for such an under-determined problem. Because it is the natural method for an under-determined problem, and because it can be readily modified to give a use ful approximation, the adjoint-equation method is the basis for the approximate form of sensitivity presented in this thesis. An abridged derivation and description of the adjoint-equation method are given below. 4.2 The adjoint-equation method for accurately calculating the sensitivities The adjoint-equation method has previously been applied to electromagnetic inverse problems by, for example, Weidelt (1975) and Madden (1990). Also, the method used in Chapter 2 to calculate the sensitivities for the one-dimensional electromagnetic inverse problem considered in that chapter is a special case of the adjoint-equation method. 39 CHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 40 The derivation presented here follows that of McGiffivray et al. (1994). Consider a do main, D, whose conductivity, o, varies only as a function of position, and whose elec trical permittivity, e, and magnetic susceptibility, t, are constant. Consider also some general source which can be defined in terms of electric and magnetic current densities J and m = —iWILM. This description of the source terms is not that of McGifflvray et al., but that used by Ward & Hohmann (1988) where M is the magnetisation. This description is consistent with Appendices A and C in which a unit electric dipole source corresponds to eI = 1 and a unit magnetic dipole source corresponds to M = 1. The electric and magnetic fields, E and H respectively, that are generated in the domain are given by the frequency-domain Maxwell’s equations: V xE = iWUH+Jm (4.1) VXH(O+W6)E+Je. (4.2) The electric and magnetic fields satisfy the general boundary conditions o(flxU) + ,13(ñxñxVxU) = Q (4.3) on the boundary, ÔD, of the domain D. U represents either E or H, c and /3 are constants, ii is the unit normal to OD and Q is the appropriate electric or magnetic surface current density. For the purposes of the inversion, the conductivity is represented as a finite linear combination of suitable basis functions: (r) = (4.4) where are the basis functions and are the coefficients. Substituting this representa tion of the conductivity into eqs. (4.1) and (4.2), and differentiating these equations with CHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 41 respect to o, one obtains OE . 8H V x = —iw-— (4.5) 00k 811 . OE V x = (tr + zwe)--— + ct’kE. (4.6) The partial derivatives OE/ôok and OH/ôok are the sensitivities of the electric and mag netic fields with respect to the coefficient ok. These sensitivities satisfy the homogeneous boundary conditions a (ñx) + (ñxñxVx) = 0. (4.7) To arrive at a numerical method for calculating the sensitivities, a suitable solution to eqs. (4.5) and (4.6) has to be found. With this in mind, consider an auxiliary problem in which auxiliary electric and magnetic fields, Et and Ht, are generated in the conduc tivity, o, by as yet undefined electric and magnetic sources J and J. These auxiliary fields satisfy the following version of Maxwell’s equations: V x = —iw1uH + J (4.8) V x Ht = (o. + iw6)Et + J. (4.9) The auxiliary fields satisfy the homogeneous boundary conditions a(ñxU) + ,Bt(nxnxVxUt) = 0. (4.10) If the boundary 9D extends to infinity, at and ,B need not be the same as a and /3. If the domain D is finite, however, a = a and 13t = /3 for the analysis that follows. Using the vector identity V.(AxB)=B.(VxA)—A.(VxB) (4.11) CHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 42 eqs. (4.5) and (4.6) can be combined with eqs. (4.8) and (4.9) to give . (t x — x Hf) = J. + J . — Et . E&k. (4.12) Integrating this expression over the domain D and using the divergence theorem gives F / OH UE I tEx——---——xHinds J8D k 0k I = ID [Jt. 0(7k + 0o — E E k] dv. (4.13) It can be shown that, for the boundary conditions specified above, the left-hand side of eq. (4.13) is equal to zero. Hence, ID (Jt. + jt. -) dv = ID E k dv. (4.14) Eq. (4.14) is the main result from which the sensitivities for the electric and mag netic fields can be obtained for a particular problem. It requires appropriate choice of sources, J and J, for the auxiliary problem. For example, to obtain the sensitivity of the x-component of the electric field at an observation location r0, choose J = and J = 0. Substituting these expressions into eq. (4.14) gives = J E E bk dv. (4.15)0k r0 D For this example, the auxiliary electric field, E, is now defined as the electric field in the domain, D, due to an x-directed unit electric dipole at the location, r0, of the observations of E. E is the electric field due to the particular source (described by e and m) used in the geophysical survey. To obtain the sensitivity of the x-component of the magnetic field at r0, set J = 0 and J = S(r — r0) ê in eq. (4.14). The sensitivity is then given by = J E .Ek dv, (4.16)D CHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 43 where E is now the electric field in the domain due to an x-directed magnetic dipole of strength (—iwt) at r0. Eq. (4.14) can be obtained by a formal Green’s function solution of eqs. (4.5) and (4.6) (see the references given at the start of this section). The auxiliary electric field, E, in this section is equivalent to the complex conjugate of the appropriate adjoint Green’s function for eqs. (4.5) and (4.6). I shall therefore use the term “adjoint field” rather than “auxiliary electric field” to refer to E. 4.3 An approximate form of sensitivity To calculate the sensitivities using the adjoint-equation method outlined in the pre vious section, the following are needed: (1) the electric field, E, in the domain due to the source used in the geophysical experiment, (2) the adjoint field, E, in the domain due to the appropriate dipole source at the observation location, and (3) the evaluation of the volume integral (weighted by the basis function) of the scalar product of these two fields. The computation time required to numerically evaluate the volume integral is negligible compared to that for the other parts of the procedure. Also, in an iterative inversion procedure the electric field, E, will already have been computed in order to calculate the data misfit resulting from the previous iteration. No additional work therefore has to be done to fulfil the first requirement above. However, for the second requirement, the adjoint field has to be explicitly computed in the conductivity model for a dipole source at each observation location in turn. It is this that accounts for the vast majority of the computation time for this method, and hence the statement that the computation time for the adjoint-equation method is equivalent to that for M0 x Mf forward modellings. In this thesis, I present an approximate form of sensitivity that is much quicker to compute than the adjoint-equation method: instead of computing the adjoint field, E, in the multi-dimensional conductivity model, an approximate adjoint field is computed that is considerably quicker to obtain. The types of approximate adjoint fields given in this thesis are the electric fields in a homogeneous halfspace, or in a horizontally-layered CHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 44 halfspace, which are generated by the appropriate dipole source. The Born approximation of the adjoint field is also considered. Computing an approximate adjoint field will be quicker than computing the true adjoint field in the multi-dimensional conductivity model. Repeating this calculation for the dipole source at each observation location will lead to considerable time-savings in the generation of the Jacobian matrix. The approximate adjoint field will obviously not contain all the features of the true adjoint field in the multi-dimensional model. How ever, as I shall illustrate in the next chapter, the true adjoint field is dominated by its decay away from the dipole source. Suitable choice of conductivities for the homogeneous or layered halfspace result in an approximate adjoint field that has the same dominant behaviour. Such an approximate adjoint field, when combined in the volume integration with the electric field, E, from the forward modelling (which is exact) produces approxi mate sensitivities that are sufficiently accurate to allow an iterative inversion procedure to converge to the desired result. Furthermore, the adjoint-equation method of exactly calculating the sensitivities is completely general. The approximate form of sensitivities presented here can therefore be used in any electromagnetic inverse problem irrespective of source-receiver geometry. In the next Chapter, I shall investigate the forms of the approximate adjoint field for the two-dimensional electromagnetic inverse problem and emphasise the comparisons between the approximate adjoint fields and the true adjoint field for a two-dimensional conductivity structure. In Chapter 6, approximate sensitivities for the two-dimensional magnetotelluric problem are calculated for a simple model and compared with the ac curate sensitivities. The approximate sensitivities are then used in the inversion of a synthetic data-set and a field data-set. In Chapter 7, the approximate sensitivities for the 2.5d problem and the three-dimensional controlled-source problem are calculated and compared with the accurate sensitivities for simple conductivity models. Chapter 5 Approximate Adjoint Fields for the Two-Dimensional Problem 5.1 Introduction The approximate form of sensitivity presented in this thesis is based on approximat ing the adjoint field required in the adjoint-equation method of accurately calculating the sensitivities. The rationale for this was given in Chapter 4. Here I investigate vari ous possibilities for the approximate adjoint field in the context of the two-dimensional problem. The possibilities are the electric field computed in a homogeneous halfspace for the appropriate dipole source, the electric field computed in a horizontally-layered halfspace, and the Born approximation of the adjoint field. The approximate adjoint fields are computed and compared with the corresponding true adjoint field for a simple two-dimensional conductivity model. 5.2 True adjoint field Consider the two-dimensional conductivity model shown in Fig. 5.1 comprising a con ductive block of 0.1 S m1 in a more resistive background of 0.01 5 m1, and a conductive basement of 0.1 S m1. (The mesh used for all calculations relating to this conductivity model is shown in Fig. 5.2.) Suppose measurements are made of the electric and magnetic fields induced in this model by a source that is also invariant in the strike direction (for example, the plane-wave source assumed in magnetotelluric experiments). For this purely two-dimensional problem, the adjoint field is the electric field due to a two-dimensional dipole source at the measurement location (see Appendix B). This source is equivalent to an infinite line of point (in three-dimensional space) dipoles which is parallel to the strike direction and which passes through the measurement location. Suppose a measurement 45 CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 46 —080 —0.95 I,. -1.09 —1.24 —1 38 —1 52 —1.67 —1.81 —1.95 —2.10 log10 cr Figure 5.1 The two-dimensional conductivity model used here and in sub sequent chapters to illustrate the approximate adjoint fields and approximate adjoint sensitivities for the two-dimensional and 2.5d inverse problems. of the along-strike electric field, E, has been made at x0 = —3000 m on the surface of the model. From eq. (4.14) and Appendix B, the sensitivity of this measurement with respect to the model parameter is given by = fEt.Ejds. (5.1) ,=—3OOO A Here, the adjoint field, Et, is the electric field due to a unit two-dimensional y-directed electric dipole source (essentially an infinite line current in the along-strike direction) at the observation location x1, = —3000 m. The true adjoint field for the two-dimensional model in Fig. 5.1 and for the two dimensional p-directed electric dipole source at = —3000 m was computed using the N 100001 I 15000 (m) CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 47 U Figure 5.2 The mesh used for all in Fig. 5.1. calculations associated with the model shown finite-element program of Unsworth, Travis & Chave (1993). The adjoint field was com puted at each node in the mesh shown in Fig. 5.2. The amplitude of this adjoint field is shown in Fig. 5.3a. The phase is shown in Fig. 5.4a. The frequency of the source was 0.2 Hz. Note the distortion of the phase due to the conductive block, and the increased decay of the amplitude and the increase in the phase as the electric field penetrates the conductive basement. 5.3 Approximate adjoint field: homogeneous halfspace The first possible form of approximate adjoint field that I present is the electric field computed in a homogeneous halfspace for the appropriate dipole source. The method used for computing the electric field generated in a homogeneous halfspace by a two- dimensional dipole is briefly described in Appendix C. Fig. 5.3b shows the amplitude of the electric field computed in a homogeneous halfspace of conductivity 0.028 S m1 for a unit two-dimensional y-directed electric dipole source at x0 = —3000 m, and for a frequency of 0.2 Hz. The corresponding phase is shown in Fig. 5.4b. The above value of 5000 1 flflflflN 15000 20000 t’3 — (TI 0 (TI 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (11 I— — (‘3 0 0 Cl’ 0 0 0 0 0 0 0 0 0 0 0 0 x(m) z(m) z(m) Figure 5.3 The amplitudes of the true and approximate adjoint fields for a directed electric line source over the two-dimensional conductivity model shown in Fig. 5.1. The location of the line source is shown by the triangle. Panel (a) shows the true adjoint field, and panels (b) to (d) show respectively the approxi mate adjoint fields computed in a homogeneous halfspace, in a layered halfspace, and using the Born approximation. The grey-scale represents log10 Ej. CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 48 —5.80 —6.17 —6.53 —6.90 —7.27 —7.63 —8.00 —8.36 —8.73 —9.10 ct c) C C C C C C C C C C C C C U, C U C — — CU CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS —10000 —15000 410. 377. 344. 310. 277. 243. 210. 176. 143. 110. Figure 5.4 The phases of the true and approximate adjoint fields for the two-dimensional conductivity model shown in Fig. 5.1. The location of the line source is shown by the triangle. Panel (a) shows the true adjoint field, and panels (b) to (d) show respectively the approximate adjoint fields computed in a homogeneous halfspace, in a layered halfspace, and using the Born approxi mation. The grey-scale represents phase in degrees. 49 5000 0 15000 10000 —5000 —10000 —15000 —20000 20000 15000 10000 5000 0 —5000 —20000 C-) z(m) z(rn) CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 50 the conductivity was obtained by averaging the logarithm of the conductivities, weighted by across-strike area, of the model shown in Fig. 5.1. The approximate adjoint field shown in Figs. 5.3b and 5.4b obviously contains no information about the conductive block nor the conductive basement. However, the dom inant behaviour of the true adjoint field is the decay in amplitude and increase in phase away from the source, which is exactly the behaviour of the approximate adjoint field computed in the homogeneous halfspace. The effect of the two-dimensional structure in the conductivity model shown in Fig. 5.1 can almost be relegated to that of a perturbation to the field that would otherwise be present in the homogeneous halfspace. It is felt that this level of agreement between true and approximate adjoint fields should be sufficient for the subsequent approximate sensitivities to be useful in an iterative inversion scheme. 5.4 Approximate adjoint field: layered halfspace The second form of approximate adjoint field is that computed in a layered halfspace. An appropriate layered halfspace can be constructed by averaging the conductivities of the multi-dimensional model in each of a series of horizontal layers. The electric field gen erated in a layered model by a two-dimensional source can be computed using the method described in Appendix A. Doing this for a unit two-dimensional y-directed electric dipole at x = —3000 m and for a frequency of 0.2 Hz gives the electric field shown in Figs. 5.3c and 5.4c. The layered halfspace used for this example was obtained by averaging the logarithms of the conductivities, weighted by across-strike area, in each horizontal layer of the model in Fig. 5.1. Just as for the approximate adjoint field computed in a homo geneous halfspace (see Figs. 5.3b and 5.4b), there is, of course, no manifestation in the field shown in Figs. 5.3c and 5.4c of the two-dimensional effects of the conductive block. However, the rapid decrease in amplitude and increase in phase of the true adjoint field in the conductive basement is well reproduced in the approximate adjoint field computed in the layered halfspace. This is not surprising since this feature of the two-dimensional CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 51 model is replicated in the layered halfspace used for computing the approximate adjoint field. 5.5 Approximate adjoint field: Born approximation The third form of approximate adjoint field presented here results from the Born approximation. Consider the integral equation describing the true adjoint field, E, in the multi-dimensional conductivity model (eq. 27, Hohmann 1988): Et(r) = E(r) + j G(r; r’) Et(r’) (r’) dv’, (5.2) where is the electric field calculated in some background conductivity model and Lio is the difference between the background and true conductivity structures. G is the tensor Green’s function for the background conductivity model. Eq. (5.2) is exact. The Born approximation for E is obtained by replacing E in the integrand by E: E(r) E(r) + j G(r; r’) . E(r’) (r’) dv’. (5.3) The elements of the tensor Green’s function, /G G(r;r’) = G?,,z ) , (5.4) \ G J are such that is the x-component of the electric field generated in the background con ductivity model by a unit y-directed electric dipole. For the E-polarisation case, eq. (5.3) reduces to a scalar equation and the only relevent component of the tensor Green’s func tion, for a homogeneous halfspace background model is given by Hohmann (1988). For the il-polarisation case, the relevent components of the tensor Green’s function (G, G & G) for a homogeneous halfspace are given by Lee & Morrison (1984). CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 52 The Born approximation for the adjoint field was calculated for the conductivity model in Fig. 5.1. The background conductivity model in which and were caldu lated was a homogeneous halfspace of 0.01 S m1. This value is the same as that for the resistive background in the model in Fig. 5.1, and was chosen to minimise the number of cells over which the numerical integration in eq. (5.3) had to be carried out. The resulting approximate adjoint field generated by a two-dimensional y-directed electric dipole source is shown in Figs. 5.3d and 5.4d. As before, a frequency of 0.2 Hz was used. From Fig. 5.3 it can been seen that the Born approximation produces a poorer approximation of the amplitude of the adjoint field than that computed in the homogeneous halfspace. This is because of the difference between the background conductivity used to compute the Born approximation and the conductivity of the homogeneous halfspace used to produce Fig. 5.3b. However, the phase of the Born approximation (Fig. 5.4d) is somewhat better than that computed for the homogeneous halfspace above about 1000 m and does give an indication of the increase in phase that the true adjoint field exhibits in the conduc tive basement, although this is not nearly as well reproduced as in the approximation calculated in the layered model. 5.6 Computation times The true and approximate adjoint fields described in this chapter were computed on a Sun Sparclo workstation. The time required to compute each of these is listed in Table 5.1. For the given values, the adjoint field was computed at the 861 nodes in the mesh shown in Fig. 5.2 (21 vertically and 41 horizontally) for the two-dimensional y-directed electric dipole source. The approximate adjoint fields calculated in a homo geneous or layered halfspace have similar computation times which are significantly less than the computation time for the true adjoint field. The Born approximation takes considerably longer to compute than the true adjoint field. This is because the Green’s function needed in eq. (5.3), that is the electric field in the homogeneous halfspace due CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 53 True adjoint field 75 Approximate adjoint field: homogeneous halfspace 1 layered halfspace 2 Born approximation 25 Table 5.1 Example computation times for the various types of two-dimensional adjoint fields. to a two-dimensional electric dipole, has to be computed for dipoles at depths within the model, not just for one at the surface. 5.7 Conclusions In this chapter, I have presented three forms of approximate adjoint field. An exam ple for each was computed for a simple two-dimensional conductivity model and compared with the true adjoint field. The approximate adjoint field computed in the homogeneous halfspace replicates the dominant features of the true adjoint field: the decay in ampli tude and increase in phase away from the dipole source. The approximate adjoint field computed using the layered halfspace matches the behaviour of the true adjoint field in the layered features within the two-dimensional model. The Born approximation of the adjoint field exhibits some of the two-dimensional features present in the true adjoint field, although it does not match the amplitudes of these features. These points are further illustrated by Figs. 5.5 and 5.6: the differences between the logarithm of the amplitude of the true adjoint field and the logarithm of the amplitude of the three approximate adjoint fields are shown in Fig. 5.5, and the differences between the phases are shown in Fig. 5.6. It is further evident from these figures that all three approximate adjoint fields are good overall approximations of the true adjoint field. The approximate adjoint field computed in the homogeneous and layered halfspaces are significantly quicker to obtain than the true adjoint field. The Born approximation, however, is substantially slower to compute than the true adjoint field. CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 54 In the next chapter, I shall use the three forms of the approximate adjoint field introduced in this chapter to compute approximate sensitivities for the two-dimensional magnetotelluric problem. Figure 5.5 The differences between the logarithm of the amplitude of the true adjoint field shown in Fig. 5.3a and the logarithm of the amplitudes of the approximate adjoint fields shown in Figs. 5.3b to 5.3d. The panels show the differences between the true adjoint field and (a) the approximate adjoint field computed in a homogeneous halfspace, (b) that computed in a layered halfspace, and (c) using the Born approximation. CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 55 0 170 0.049 —0.072 —0.193 —0.314 —0.436 —0.557 —0.678 —0.799 —0.920 ——20000 z(m) S z(m) CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 56 Figure 5.6 The differences between the phase of the true adjoint field shown in Fig. 5.4a and the phase of the approximate adjoint fields shown in Figs. 5.4b to 5.4d. The panels show the differences between the true adjoint field and (a) the approximate adjoint field computed in a homogeneous halfspace, (b) that computed in a layered halfspace, and (c) using the Born approximation. C.) z (in) 90.0 73.3 56.6 40.0 23.3 6.7 —10.0 —26.7 —43.3 —60.0 S z(m) Chapter 6 Approximate Sensitivities for the Two-Dimensional Magnetotelluric Inverse Problem 6.1 Introduction In the previous chapter, I compared three forms of approximate adjoint field that could be used for a two-dimensional inverse problem. In this chapter, I use these adjoint fields to calculate approximate sensitivities for the two-dimensional magnetotelluric prob lem. These sensitivities are compared with sensitivities calculated using the brute-force (or perturbation) method. In the final two sections of this chapter, the approximate sensitivities are used within a minimum-structure, Gauss-Newton algorithm to invert synthetic and field data. 6.2 Calculation of the sensitivities 6.2.1 Sensitivities for apparent resistivity and phase The observed quantities in a magnetotelluric experiment are usually considered to be the frequency-domain values of apparent resistivity, pa, and phase, ç. For a two dimensional Earth, for which the magnetotelluric problem can be split into two decoupled modes, the apparent resistivity and phase are given by paQ”) = and b(w) = phase () (6.1) for the E-polarisation mode, and by pa(W) = and (w) = phase () (6.2) 57 CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 58 for the il-polarisation mode, where w is the angular frequency and t is the magnetic permeability. The strike direction is assumed to lie in the y-direction. The sensitivities for the apparent resistivity and phase can be expressed in terms of the sensitivities for the electric and magnetic fields (see Appendix D): = 2Pa {e () - e iL) } (6.3) and (1OE”\ (1OH’\ = sm - (6.4) These relationships therefore enable approximate sensitivities for the apparent resistivity and phase to be calculated from the approximate sensitivities for the electric and magnetic fields obtained using the method described in Chapter 4. 6.2.2 Area integration For the two-dimensional magnetotelluric inverse problem, the sensitivities for the electric and magnetic fields have the form (see section 5.2) = jEt.Ezi&jds (6.5) where F represents one component of the electric or magnetic field (either E, E, H or Hr). For the typical situation in which the model comprises rectangular cells of uniform conductivity, that is the basis function is equal to 1 within the j” cell and equal to zero everywhere else, the area integration in eq. (6.5) becomes OF = J J E . E dz dx, (6.6)9o-j 0 ZO CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 59 Z (x0,z) (x1,z0) (a0,z1) (xj,z1) Figure 6.1 The notation used for the bilinear interpolation and integration of the electric and adjoint fields over a cell in a two-dimensional conductivity model. where, for simplicity of notation, the th cell is assumed to extend from z0 to x and from z0 to z1 as shown in Fig. 6.1. The forward-modelled electric field E will usually only be known at the nodes in the model (that is, at the vertices of the cell shown in Fig. 6.1). To calculate the integral in eq. (6.6), the real and imaginary parts of the components of E are bilinearly interpolated over the rectangular cell: E(x,z) = c1xz + c2z + c3x + c4. (6.7) E represents the real or imaginary part of the component of the electric field. For a general cell, the 4 x 4 system of equations relating the values of E at the vertices of the cell to the coefficients c, i = 1,... , 4, in the above expansion was solved analytically. The resulting expressions are then used to calulate the coefficients for a particular cell from the values of E at the vertices of that cell. The adjoint field, E is bilinearly interpolated in an analogous manner. With the real and imaginary parts of the components of E and E in the form of eq. (6.7), the integration with respect to x and z in eq. (6.6) could also be carried out analytically for a general cell. The resulting expression is a function CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 60 of the coefficients in the bilinear expansions of E and E. For a particular cell, therefore, the integral in eq. (6.6) can be evaluated, using analytic formulae, from the values of the electric and adjoint fields at the vertices of the cell. It was found that the horizontal discretisation of the mesh (see Fig. 5.2) used to calculate the sensitivities in the next two sections was too coarse for the above technique to accurately evaluate the integral over the cells close to the dipole, especially for the H-polarisation mode. The adjoint field for this mode varies more rapidly near the source than the adjoint field for the E-polarisation mode shown in Figs. 5.3 and 5.4. The cells in the mesh were therefore sub-divided horizontally and the forward-modelled electric field (which does not vary rapidly horizontally) linearly interpolated over these sub-divisions. The adjoint field was explicitly computed at the vertices of the sub-divisions. This re quired an insignificant number of additional computations. The technique described in the previous paragraph was then used to evaluate the area integration over each sub division. 6.3 Comparison of approximate and accurate sensitivities 6.3.1 Homogeneous halfspace Programs were written to compute approximate sensitivities for the two-dimensional magnetotelluric problem using each of the three forms of approximate adjoint field de scribed in Chapter 5. The programs were first tested by computing the approximate sensitivities for a homogeneous halfspace. For this special case, the approximate adjoint fields are not in fact approximate but exact, and so the resulting sensitivities are exact. Consider a homogeneous halfspace of conductivity 0.01 S m1 with the same hori zontal and vertical dimensions as the model shown in Fig. 5.1. Suppose measurements of the E-polarisation apparent resistivity and phase had been made at x, = —3000m on the surface of the model at a frequency of 0.2 Hz. Sensitivities were calculated for these obser vations using the brute-force (or perturbation) method. The model was divided into the 800 cells (20 vertically and 40 horizontally) shown in Fig. 5.2. The conductivity of each CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 61 cell was perturbed in turn by 10% (this perturbation is large, but was required to obtain accurate values for the cells far from the observation location) and a forward modelling performed to calculate the resulting changes in the apparent resistivity and phase. The sensitivity with respect to the perturbed model parameter is then given by the change in the apparent resistivity or phase divided by the change in the model parameter. (The forward modelling was carried out using the transmission-surface modelling program of Madden 1972.) These sensitivities are shown in Figs. 6.2a and 6.2b. The sensitivity is plotted as a constant value over the area of the cell to which it refers. Approximate sensitivities were calculated for the homogeneous conductivity model described above. The forward-modelled electric field was computed using the program of Madden (1972). All three forms of approximate adjoint field were used, and the resulting sensitivities compared with the brute-force sensitivities shown in Figs. 6.2a and 6.2b. The E-polarisation sensitivities calculated using the adjoint field computed in a layered halfspace (all the layers of which had the same conductivity) are shown in Figs. 6.2c and 6.2d. There is very good agreement between the brute-force sensititivities and the sensitivities calculated using the adjoint field in the layered halfspace. The approximate sensitivities calculated using the adjoint field computed in a homogeneous halfspace and those calculated using the Born approximation of the adjoint were indistinguishable from those in Figs. 6.2c and 6.2d. A comparison was also made of the brute-force and approximate sensitivities for the H-polarisation mode. The brute-force sensitivities are shown in Figs. 6.3a and 6.3b. The approximate sensitivities calculated using the adjoint field in a layered halfspace are Figure 6.2 (Following page) Sensitivities for the E-polarisation apparent resis tivity and phase for a homogeneous halfspace of conductivity 0.01 S m1 and for a frequency of 0.2 Hz. The observation location is indicated by the tri angle. Panels (a) and (c), and the colour-bar on the left of the figure, show log10ãlnp/Olno. Panels (b) and (d), and the colour-bar on the right, show log8q5/6lncr. Panels (a) and (b) were produced by the brute-force method, and panels (c) and (d) were produced by the approximate-sensitivity program using an adjoint field computed in a layered halfspace. CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 62 o L) 0 0 C C o — — — — CV I I I CV N — tf) C CD 0 - CV CV CV C) 0 I I I I I 20000 15000 10000 5000 0 —5000 —10000 —15000 —20000 20000 16000 10000 —5000 —10000 —15000 —20000 I 1; 1 o N C N C N CC’ I C’- 0 CV N 0 -j CV CV CV C’) C’) C’) C’) C’) ‘ ‘ 0 —5000 —10000 —15000 —200000 L’ 0 0 o 0 o o z(m) CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 63 shown in Figs. 6.3c and 6.3d. Again there is good agreement between the brute-force and approximate sensitivities. 6.3.2 Two-dimensional model Approximate sensitivities were calculated for the two-dimensional model shown in Fig. 5.1. As in the section above, measurements of the apparent resistivity and phase made at Xc, = —3000 m on the surface of the model at a frequency of 0.2 Hz were considered. The E-polarisation brute-force sensitivities are shown in Figs. 6.4a and 6.4b, and the approximate sensitivities in Figs. 6.4c and 6.4d. The approximate sensitivities were calculated using the approximate adjoint field computed in a layered halfspace. (For this example, the approximate adjoint field is that shown in Figs. 5.3c and 5.4c.) There is good agreement between the brute-force and approximate sensitivities shown in Fig. 6.4. A comparable, although slightly poorer, level of agreement was found for the two other forms of approximate sensitivity. It should be noted that, because the sensitiv ities shown in Fig. 6.4 are with respect to ln a, they depend directly on the conductivity model as well as the behaviour of the electric and adjoint fields within the model. (This Figure 6.3 (Following page) The H-polarisation sensitivities for the homoge neous halfspace. Panels (a) and (c), and the colour-bar on the left of the figure, show log10jOlnp/Oln4 Panels (b) and (d), and the colour-bar on the right, show logOq/Olno. Panels (a) and (b) were produced by the brute-force method, and panels (c) and (d) were produced by the approximate-sensitivity program using an adjoint field computed in a layered halfspace. Figure 6.4 (Page 65) Brute-force and approximate sensitivities for the E polarisation apparent resistivity and phase for the conductivity model shown in Fig. 5.1. The observation location is indicated by the triangle. The fre quency was 0.2 Hz. Panels (a) and (c), and the colour-bar on the left of the figure, show log10Olnp/O1ncrj, and panels (b) and (d), and the colour-bar on the right, show log10Oq/Olnoj. Panels (a) and (b) show the sensitivities produced by the brute-force method, and panels (c) and (d) show the approxi mate sensitivities calculated using the approximate adjoint field computed in a layered halfspace. CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 64 0 If) C) — C) f) 0 C D C C- — If) C) C’) 0 - — CIf CIf CIf C’) C’) C’) I I I I I I I’ 0 0 0 0 o 0 If) C — C’) 20000 15000 10000 5000 0 —5000 —10000 —15000 —20000 -d C) z(m) z(m) ii 0 C’- 0 C’- 0 C’- C C’) CC C) C’) CC C) C’) CD C) C’) ,_i - — C C’) C’) C’) C’) C’) I I I I I I I CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT I Ii C LID OD ID 0cc Q C’) ID C C’) IL) C o———— C’2 CQ CIL C’ I I I I I I I I I 11 o II) 0) ID C’) C’- — 11) 0C C’) C C’- — ‘ — CQ CIL CL) C’) C’) C’) ‘ ‘ ‘ I I I I I I I I 65 0 0 ID z(m) z(m) CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 66 can be seen from the relation O/O(ln o) = o 6/0o.) The dramatically increased sensi tivity in the conductive block that results from this, and the enhanced sensitivity in the conductive basement, can be clearly seen in Fig 6.4. This is the two-dimensional mani festation of what Boerner & Holladay (1990) observed in their “fractional” sensitivities for the one-dimensional problem (see section 3.3.1). The H-polarisation sensitivities were also calculated for the two-dimensional con ductivity model. The brute-force sensitivities are shown in Figs. 6.5a and 6.5b, and the approximate sensitivities calculated using the adjoint field in a layered halfspace are shown in Figs. 6.5c and 6.5d. The agreement between the brute-force and approximate sensitivities is similar to that for the E-polarisation sensitivities. Again, the two other forms of approximate sensitivity gave a slightly poorer level of agreement. 6.4 Computation times Table 6.1 lists the computation times for the various methods of calculating the Jacobian matrix of sensitivities for a small two-dimensional magnetotelluric example. The sensitivities were calculated for 80 data (observations of apparent resistivity and phase for both E- and H-polarisations at four locations and at five frequencies) with respect to the conductivities of the 800 cells used to parameterise the model in Fig. 5.1. The approximate adjoint fields in both the homogeneous and layered halfspaces were only explicitly computed for one observation location and then translated horizontally to calculate the sensitivities for the three remaining observation locations. This results in a considerable time-saving that becomes more significant with increasing number of observation locations. In a similar manner, for the Born approximation of the adjoint Figure 6.5 (Following page) Brute-force and approximate sensitivities for the H-polarisation mode for the conductivity model shown in Fig. 5.1. Panels (a) and (c), and the colour-bar on the left of the figure, show log10Olnp/ôlncr, and panels (b) and (d), and the colour-bar on the right, show log10 q/Olnu. Panels (a) and (b) show the sensitivities produced by the brute-force method, and panels (c) and (d) show the approximate sensitivities. CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT o a a a a a a a a a c — q c’ CO a a a — — - C\ C’ 01 01 I I I I I I I ItJ *1 I 1 o 0) 0 CO U) C) 01 — 00) 0) 0 C’) 0) C — — CU CU C’) C LI) I I I I I I I I 67 I a a 0 U) z(m) a a a a a U) z(m) CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 68 field, both and G, which were calculated in a homogeneous halfspace, were explicitly computed for only the first observation location and then translated horizontally for the other observation locations. This reduces the computation time for this third form of approximate sensitivities, although it is still not less than that of the accurate sensitivities. It is not until there are tens of observation locations that the approximate sensitivities based on the Born approximation of the adjoint field require less time to calculate than the accurate sensitivities. Accurate sensitivities: brute-force 500 mm adjoint-equation 2.5 Approximate sensitivities: homogeneous halfspace 0.25 layered halfspace 0.5 Born approximation 5 Table 6.1 Example computation times for the various types of sensitivities for the two-dimensional magnetotelluric inverse problem. 6.5 Two-dimensional inversion of magnetotelluric data So far in this chapter, I have been concerned with the comparison of the values of the approximate and accurate sensitivities. Although this is important, the usefulness of the approximate sensitivities ultimately depends on their successful performance within an inversion routine. In this section, therefore, I test the behaviour of an iterative, linearised inversion procedure that uses the approximate sensitivities. The inversion strategy discussed in Chapter 2 is used to invert synthetic and field magnetotelluric data to recover two-dimensional conductivity models. The sensitivities are just one of several components in an inversion procedure, and so the effectiveness of the approximate sensitivities will depend on the details of the inversion strategy used: convergence of one algorithm using approximate sensitivities does not mean that all algorithms will converge. However, testing approximate sensitivities in all CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 69 possible algorithms is clearly impractical. Here I consider an algorithm that is typical of those that one would like to apply to the multi-dimensional electromagnetic inverse problem. 6.5.1 Synthetic data Approximate sensitivities, calculated using the approximate adjoint field in a layered halfspace, were used to invert a set of synthetic data generated from the conductivity model in Fig. 5.1. Values of the E- and H-polarisation apparent resistivity and phase were calculated at four observation locations (Xe —11000, —3000,5000 and 13000m) for five frequencies (1, 0.5, 0.2, 0.1 and 0.05 Hz). Gaussian random noise was added to these data. The standard deviation of the noise added to the apparent resistivity was 5%, and the standard deviation of the noise added to the phase was 2 degrees. The data are shown in Fig. 6.6. The 80 synthetic data were inverted using an iterative, linearised, minimum-structure inversion procedure in which the system of equations at each iteration was solved using a subspace technique. The details of the inversion algorithm are given in Oldenburg & Ellis (1993). The model norm that was minimised was the discrete equivalent of = fa3w(m — rn0)2 + JA 1 /O(m_mo)2 ______ cw Ox ) + wz ) ) ds, (6.7) where the model, m, was equal to ln o, and the reference model, m0, corresponded to a halfspace of conductivity 0.01 5 m1. The values of the coefficients were = 10—2, = 3 and c = 1. The weighting functions, w8, w and w were constant and equal to unity throughout the model, except for w, which was equal to 10 in the top four layers of the model. This was to suppress structure developing in the near-surface layers as a CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT Figure 6.6 The synthetic data (shown by the circles and error bars) generated from the model shown in Fig. 5.1, and the predicted data (shown by the lines) from the final model produced by the inversion of the synthetic data. The filled circles and solid lines correspond to the E-polarisation data, and the open circles and dashed lines to the H-polarisation data. The units of the apparent resistivity, Pa are fm and those of the phase, q, are degrees. 70 70 102 102 1 2 102 102 Ct’ q a : ,a E ‘. - 0.05 Hz —15000 0 15000 j Tz —15000 0 15000 rZ —15000 0 15000 0.5Hz —15000 0 15000 1Hz — I - 0.05 Hz E —15000 0 15000 70 0.1 Hz 60 - J —15000 0 15000 70 0.2 Hz 60- 5 - f---;_; 40 p —15000 0 15000 70 - 0.5 lIz —15000 0 15000 70 - 1 Hz : x(m) 40 —15000 0 15000 —15000 0 15000 x(m) CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 71 result of the large H-polarisation sensitivities for these layers (see Figs. 6.3 and 6.5). The subspace vectors were the steepest-descent vectors vk = (Wm)’ Vm (6.8) where YLm is the weighting matrix associated with the discrete equivalent of eq. (6.7), and is the data misfit for the kt partition of the data. The data were ordered according to their individual misfit and then divided into, in this case, 20 partitions. The approximate sensitivities were calculated using the approximate adjoint field in a layered halfspace. The conductivities of the layers were chosen by averaging the logarithms of the conductivities within each layer of the model. The final model is shown in Fig. 6.7. The conductive block has approximately the correct amplitude and is situated at the correct depth and horizontal location. The conductive basement has also been recovered. The predicted data produced by this model are shown in Fig. 6.6. The values of the data misfit, model norm and Lagrange multiplier during the inver sion are shown in Fig. 6.8. The data misfit decreases at the requested rate (/3 = 0.8) for the first seventeen iterations. The inversion algorithm can then only decrease the misfit more slowly than the requested rate. The target misfit of 80 is achieved after 31 iterations and the misfit then remains constant for the remaining iterations. The model norm grad ually increases before levelling off once the algorithm has achieved the target data misfit. The Lagrange multiplier steadily increases for the first twenty iterations, but then varies considerably for the next ten. This coincides with the iterations for which the algorithm cannot decrease the misfit at the requested rate. It appears as if the algorithm is having to work hard to lower the misfit to the target value. Once the target misfit has finally been achieved, the Lagrange multiplier levels off and shows only small variations as the algorithm attempts to further minimise the model norm. The results of the above inversion of synthetic data are what would be expected if accurate sensitivities had been used in place of the approximate sensitivities: the final CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT —0.80 —095 —1.09 —1.24 —1.38 —1.52 —1.67 —1.81 --1.95 —2.10 log10 Figure 6.7 The final model produced by the inversion of the synthetic data shown in Fig. 6.6. The model used to generate these data is the one shown in Fig. 5.1. The colour-bar is the same for both figures. The triangles indicate the observation locations. 72 model reproduces the data to the desired level of misfit and has the typical smeared- out appearance of a model produced by a minimum-structure inversion algorithm. The behaviour of the data misfit, model norm and Lagrange multiplier (with the possible exception of the oscillations between iterations 24 and 30) is also what would have been expected if accurate sensitivities had been used in the inversion procedure. As a final comment for this section, if the approximate sensitivities of Smith & Booker (1991) had been used to invert the synthetic magnetotelluric data considered in this section, the model could only have contained four columns of cells, one beneath each observation location. The horizontal resolution of the final model would therefore have been considerably less than that in Fig. 6.7. N 10000 I S— 15000 x(m) CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 73 a 102 - 0 10 20 30 40 i o - io 101 b 100 0 10 20 30 40 1 o3 100 0 0 10 20 30 40 L er a Li on Figure 6.8 The data misfit, q, model norm, 1m’ and Lagrange multiplier, , as functions of iteration during the inversion of synthetic data discussed in section 6.5.1. 6.5.2 COPROD2 data To conclude this chapter, the approximate sensitivities were used in the inversion of a set of field data. The data-set chosen was a sub-set of the COPROD2 data-set collected in southern Saskatchewan and Manitoba. A description of the COPROD2 data-set and their context is given by Jones (1993). The data are dominated by the effects of the North American Central Plains (“NACP”) anomaly: a linear zone of high conductivity at mid crustal depths that runs north from South Dakota through Saskatchewan before turning eastwards under Hudson Bay. The anomaly appears in the COPROD2 data as low values of the E-polarisation apparent resistivity and high values of the E-polarisation phase. In contrast, the H-polarisation measurements show little indication of the anomaly. CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 74 The COPROD2 data have been used by the magnetotelluric community as a test data-set for the comparison of two-dimensional inversion schemes (see Ingham, Jones & Honkura 1993). The considerable splitting of the E- and H-polarisation modes mentioned above makes the COPROD2 data-set a challenging test for any inversion procedure, es pecially considering the size of the problem is close to the limit of what can presently be contemplated. All the inversion programs thus far applied to the COPROD2 data-set have struggled to some extent, either because of the nature of the data or because of their considerable number. The sub-set of the COPROD2 data chosen for consideration here is the same as that used by Ellis, Farquharson & Oldenburg (1993), except that stations to the east of 118 km, and so relating to the Thompson belt (“TOBE”) anomaly, were omitted. The total number of data was therefore 896. The error estimates were the same as those of Ellis et al. The observed values of apparent resistivity and phase are shown in Fig. 6.9. For the inversion carried out here, the model was divided into 30 layers and 62 columns as shown in Fig. 6.10. The columns were arranged so there was one column beneath each observation location and an additional column between each pair of stations. There were therefore just over twice as many columns as there were stations. The resulting number of cells was 1860. The inversion algorithm was the same as that used in section 6.5.1. Approximate sensitivities were calculated using the approximate adjoint field computed in a layered halfspace. The conductivities of the layers of this halfspace were obtained by averaging the logarithms of the conductivities in each layer of the two-dimensional model. A one- dimensional inversion of the complete data-set was carried out to produce the best-fitting horizontally-layered model. This model, which has a conducting surface layer of 0.5 S m1 that ramps down to a resistive basement of 0.01 S m1 below about 7 km, was then used as the reference model in the two-dimensional inversion. This reference model was not included in the horizontal and vertical gradient terms (see eq. 6.7) in the weighting matrix Lm The subspace vectors were again obtained by partitioning the data in terms of misfit, but 200 partitions were used here rather than the 20 in section 6.5.1. CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 75 102 0.023 Hz 60 0.023 Hz I I ______ —200 —100 0 100 —200 —100 0 100 102 0.047 Hz 60 0.047 Hz I__I__I__I__I _ ____ _ _ —200 —100 0 100 —200 —100 0 100 102 0.093 Hz 60 0.093 Hz 101 20 100 0 I —200 —100 0 100 —200 -100 0 100 1 2 0.19 Hz 60 101 ° 100 _____________ 20 - 0.19 liz —200 —100 0 100 —200 —100 0 100 x (km) x (km) Figure 6.9 The sub-set of the COPROD2 data considered in section 6.5.2. The observations are represented by the circles, crosses and error bars. The lines indicate the predicted data for the final model produced by the inversion (see Fig. 6.11). The filled circles and solid lines correspond to the E-polarisation data, and the crosses and dashed lines to the il-polarisation data. The units of the apparent resistivity, Pa are m and those of the phase, q, are degrees. CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 76 1 2 100 102 60 40 20 —200 —100 0 100 60 40 20 102 CO 10 100 —200 1 2 CO LU 100 —100 —200 —100 0 100 0 100 x (km) 60 40 20 - 0.0059 Hz 0— 200 60 ‘10 - 20 0 100 —200 —100 0 100 0 100 x (km) CO hi —200 —100 0 100 0 0.0029 lIz iOo I I —200 —100 0 100 0.0029 1-Tz 0 —200 —1 00 0 100 0.0059 Hz I I I I 0.012 Hz - 0.012 Hz Figure 6.9 (contd.) CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 77 0 10 ______ 20 30 - ‘-‘ 40 50 60 I I I 0 — — — (31 0 0 o a a 0 o a a x (km) Figure 6.10 The mesh for the model used in the inversion of the COPROD2 data. The final model is shown in Fig. 6.11. The vertical exaggeration is 2.9:1 (this is different from Fig. 6.11). As mentioned at the start of this section, the COPROD2 data-set has proved to be challenging to many inversion programs. This was also the case here. The main difficulty was that both the E- and H-polarisation sensitivities were large for the uppermost layers in the model. The values of w3 and w for these layers were therefore increased to counteract the development of excessive structure that these large sensitivities would otherwise cause. The process of choosing the most effective form of this additional weighting was empirical. I therefore ran the inversion program numerous times, each using a different form of weighting in the near-surface layers. Once the misfit had stopped decreasing, I would alter the weighting to see if a further decrease in the misfit could be achieved. The final result of the inversion process was the model shown in Fig. 6.11. The predicted data for this model are shown in Fig. 6.9, and the corresponding value of the x2 misfit is 1370. This is greater than the target value of 896. However, the final misfit is close to this target misfit, and the comparison between the observed and predicted data in Fig. 6.9 is good. CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 78 In summary, therefore, the inversion of the CQPROD2 data using the approximate sensitivities has produced a model that reproduces the observations to almost the desired level of misfit. The final model also contains the multiple, distinct, mid-crustal conductors that have been produced by other inversions of the COPROD2 data. 6.6 Conclusions In this chapter, I have calculated approximate sensitivities for the two-dimensional magnetotelluric inverse problem. These sensitivities were compared with brute-force sen sitivities for a two-dimensional conductivity model. It was shown that the agreement be tween the approximate and brute-force sensitivities is good, especially when considering the sensitivity with respect to the logarithm of the conductivity. The level of agreement, as illustrated in Figs. 6.6 and 6.7, seems sufficient to enable iterative, linearised inver sions of magnetotelluric data which use these approximate sensitivities to converge to an acceptable model. This was emphasised by successful inversions of a synthetic data-set and of the COPROD2 data. It was also shown that the approximate sensitivities are considerably quicker to compute than even the most efficient method of accurately calculating the sensitivities. For the small example considered in section 6.4, a factor of five difference was obtained. This factor will increase with the size of the problem, since the approximate adjoint field is calculated for the dipole at only one observation location and then translated for all other observation locations. Calculating the accurate sensitivities using the adjoint equation method requires the adjoint field to be computed for a dipole at each observation location. The difference-factor for the inversion of the COPROD2 data in section 6.5.2 was estimated to be two orders of magnitude. Figure 6.11 (Following page) The final model produced by the inversion, using the approximate sensitivities, of the sub-set (shown in Fig. 6.9) of the COPROD2 data. The vertical exaggeration is approximately 2.3:1. () I 0. 90 0. 35 — 0. 19 — 0. 74 — 1. 28 — 1. 82 — 2. 37 — 2. 91 — 3. 45 — 4. 00 c (km ) lo g 1 0 a Chapter 7 Approximate Sensitivities for the 2.5d and Three-Dimensional Inverse Problems .1 Introduction In Chapter 6, I investigated approximate sensitivities in the context of the two- dimensional magnetotelluric problem. I presented a comparison of the approximate and accurate sensitivities for a simple two-dimensional conductivity model. This comparison covered all the observed quantities relevent to the two-dimensional magnetotelluric prob lem: measurements of apparent resistivity and phase for both the E- and il-polarisations. It was shown that there was good agreement between the accurate and approximate sen sitivities. It was also shown that the computation time for the approximate sensitivities was considerably less than for any of the accurate methods of calculating the sensitiv ities. Finally, in the last part of Chapter 6, I presented the results of two inversions using the approximate sensitivities, one of synthetic data and one of field data (a sub-set of the COPROD2 data). The approximate sensitivities performed successfully in both inversions. In this chapter, I consider the approximate sensitivities for the 2.5d problem and the three-dimensional controlled-source problem. The main goal is to construct com putational algorithms to calculate the approximate sensitivities for these more realistic and widespread problems, and to determine how much faster the approximate sensitiv ities are to compute than the accurate sensitivities. A complete investigation similar to that in Chapter 6, including test inversions using the approximate sensitivities, was not contemplated. It was assumed that if a similar level of agreement to that in Chapter 6 was obtained between the approximate and accurate sensitivities, then the approximate 80 CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 81 sensitivities would perform as well in a 2.5d or three-dimensional inversion as they did in the two-dimensional inversions of magnetotelluric data. In the first half of this chapter, I consider the approximate sensitivities for the 2.5d inverse problem in which the electric and magnetic fields vary in all three dimensions but the conductivity model varies in only two-dimensions (the vertical direction and one horizontal direction). The approximate sensitivities are compared with accurate sensitivities for a homogeneous halfspace and for a simple conductivity model, and the difference in computation time between the approximate and accurate sensitivities is demonstrated. In the second half of this chapter, I consider the approximate sensitivities for the three-dimensional controlled-source problem. This is the most general geophysical elec tromagnetic inverse problem and is the one we ultimately want to solve on a routine basis. However, the extent to which a comparison of approximate and accurate sensitivi ties could be made was severely restricted by the limitations placed on three-dimensional forward modelling by current computer technology. The only feasible comparison was therefore for a model in which a confined region had been parameterised into 5 x 5 x 5 cells. Using this parameterisation, the approximate and accurate sensitivities were com pared for a homogeneous model and for a simple three-dimensional conductivity structure. These comparisons also served to illustrate the considerable difference in computation times for the approximate and accurate sensitivities. 1.2 Calculation of the sensitivities for the 2.5d problem To calculate the sensitivities for the 2.5d problem, the general expression given by eq. (4.14) is first manipulated to take advantage of the particular nature of the problem (Unsworth & Oldenburg, 1995). As for the two-dimensional magnetotelluric problem discussed in Chapter 6, the model for the 2.5d problem is usually parameterised in terms of rectangular cells of uniform conductivity. These cells are infinite in extent in the strike CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 82 direction (assumed as before to lie in the y-direction). The expression for the sensitivity can therefore be written as = jfEt.Edsdy, (7.1) oo A3 where F represents the appropriate component of the electric or magnetic field (Es, E, E, H, H or H) and A is the area of the cell (in the across-strike plane). The forward modeffing for the 2.5d problem was carried out using the program of Unsworth, Travis & Chave (1993). This program essentially solves a two-dimensional problem for each in a sequence of along-strike wavenumbers k. The required electric or magnetic field is then obtained by performing the one-dimensional inverse Fourier transform in k. The forward-modelled electric field, E, in eq. (7.1) can therefore be written as E(x,y,z) = _j E(x,k,z) dk, (7.2) where y8 is the y-coordinate of the source. E(x, k, z) emerges naturally from the program of Unsworth et al. The approximate adjoint field, E, can be expressed in a similar manner: Et(x,y,z) = — j Et(x,1c,z) eio) dk, (7.3) where y0 is the y-coordinate of the observation location. Et(c, ks,, z) is easily obtained when it comes to considering an approximate adjoint field. Substituting the above two expressions for the forward-modelled and adjoint fields into eq. (7.1) gives = j-- J J f A(k,k) ei eio) dk dk dy, (7.4)3 —00 —00 —00 where the area integration, A, in the along-strike-wavenumber domain is defined as A(k,k) = f Et(x,k,z) .E(,k,z) ds. (7.5) CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 83 Using the Fourier transform definition of the Dirac delta function: S(k) = — j: e’ dy, (7.6) the integration with respect to y in eq. (7.4) can be performed to give = (7.7) assuming, for simplicity, that y = y0 = 0. This equation reduces to OF 1[ = (7.8) To calculate the sensitivities for the 2.5d problem, therefore, the area integration in the along-strike-wavenumber domain (eq. 7.5) is evaluated at each of the wavenumbers re quired by the forward modeffing. The area integration at each wavenumber is carried out using the same method as for the purely two-dimensional problem described in sec tion 6.2.2. The sensitivity for a particular cell is then obtained by the integration over the wavenumbers in eq. (7.8). 7.3 Comparison of approximate and accurate 2.5d sensitivities 7.3.1 Homogeneous halfspace As for the two-dimensional magnetotelluric problem, the program for computing the approximate sensitivities was first tested on a homogeneous halfspace conductivity model. For this simple model, the sensitivities calculated using this program should be equal to the accurate sensitivities. Consider a homogeneous halfspace of conductivity 0.01 Sm’ having the same di mensions as the model shown in Fig. 5.1 and parameterised using the mesh in Fig. 5.2. Consider a controlled-source electromagnetic experiment in which an along-strike electric dipole source (a point source in three-dimensional space) is located at x3 = —11000m and = = 0. Suppose measurements are made of E at x0 = l000m and y0 = = 0, and CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d S 3d 84 at a frequency of 0.2 Hz. The accurate sensitivities for this arrangement were calculated by both the brute-force (or perturbation) method using the forward-modelling program of Unsworth, Travis & Chave (1993) and the adjoint-equation method (Unsworth & Old enburg 1995). The sensitivities calculated using the latter method are shown in Figs. 7.la and 7.lb. The sensitivity is plotted as a constant value over the area of the cell to which it refers. The program written to compute approximate sensitivites for the 2.5d problem was used to calculate sensitivities for the homogeneous halfspace. The experimental configura tion was the same as that described above. The adjoint field was that due to a y-directed unit electric dipole and was computed in a layered halfspace using the method described in Appendix A. The resulting sensitivities are shown in Figs. 7.lc and 7.ld. There is very close agreement between the accurate sensitivities calculated using the adjoint-equation method and the sensitivities calculated using the approximate-sensitivity program (which should equal the accurate sensitivities for this special case of the homogeneous halfspace). 7.3.2 Two-dimensional model Approximate sensitivities were calculated for the two-dimensional model shown in Fig. 5.1. The controlled-source experimental configuration used in the previous section was also considered here. The accurate sensitivities calculated using the adjoint-equation Figure 7.1 (Following page) Sensitivities for the 2.5d problem described in sec tion 7.3.1. The conductivity model was a homogeneous halfspace of 0.01 5m’. The source and observation locations are indicated by the open and solid trian gles respectively. Panels (a) and (c) show log10Oln Ej/Oln °i, and panels (b) and (d) show log10 0 ç/0ln crj. IE and q are the amplitude and phase, respec tively, of the electric field. The phase is measured in radians. The colour-bar refers to all four panels. Panels (a) and (b) were produced by the adjoint equation method, and panels (c) and (d) were produced by the approximate sensitivity program using an adjoint field computed in a layered halfspace. CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d C (‘2 CD C — C) CD CD Ca iD C U C CD —. CD — 1’- (‘1 (‘2 C) ) CD CD CD CD I I I I I I —5000 —10000 —15000 85 Q C o o o o CD 0 — Cv —20000 0 0 0 0 0 o o 0 0 o o 0 0 CD 0 CD 0 — — Cv —20000 —5000 —10000 15000 c.) CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 86 method are shown in Figs. 7.2a and 7.2b. The approximate sensitivities calculated using an approximate adjoint field in a layered halfspace are shown in Figs. 7.2c and 7.2d. There is good agreement between the approximate and accurate sensitivities shown in Fig. 7.2. Because these sensitivities are with respect to the logarithm of the conductivity, the direct effect of the conductivity structure is clearly evident, just as for the two- dimensional magnetotelluric sensitivities shown in Chapter 6. The agreement between the approximate and accurate sensitivities for the 2.5d problem is better than that for the two-dimensional magnetotelluric sensitivities. This is because the adjoint field is now fully three-dimensional and so is even more dominated by its decay away from the dipole source compared to its decay away from the line source in two dimensions. 7.4 Computation times for the 2.5d problem The computation times for the accurate and approximate sensitivities for a small 2.5d example based on the model in Fig. 5.1 are listed in Table 7.1. These times are for observations of the amplitude and phase of E at six observation locations at a frequency of 0.2 Hz and for the 800 cells in the model. For this small example, the computation of the approximate sensitivities was quicker by a factor of five than the most efficient method of calculating accurate sensitivities. This factor will increase further as the size of the problem increases, just as for the magnetotelluric sensitivities discussed in Chapter 6 (see section 6.4). Figure 7.2 (Following page) 2.5d sensitivities for the conductivity model shown in Fig. 5.1. The experimental configuration was the same as that considered in section 7.3.1. The source and observation locations are indicated by the open and solid triangles respectively. Panels (a) and (c) show log1öln E/Olnoi, and panels (b) and (d) show log108q5/Olnu. The colour-bar refers to all four panels. Panels (a) and (b) were produced by the adjoint-equation method, and panels (c) and (d) were produced by the approximate-sensitivity program using an adjoint field computed in a layered halfspace. CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d C) D.? f) C’- o . CO CD U) C)C I! C IC) C CD CD ‘—u C’2 CQ CO Co U) IL) CD CD I I I I I I I I p20000 15000 10000 5000 —5000 —10000 —15000 —20000 87 ?jj WI z(m) z(m) CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 88 Accurate sensitivities: brute-force 5000 mm adjoint-equation 25 Approximate sensitivities: layered halfspace 5 Table 7.1 Example computation times for the accurate and approximate sen sitivities for the 2.5d inverse problem. 7.5 Calculation of the sensitivities for the three-dimensional problem The general expression for the sensitivity given by eq. (4.14) applies without simpli fication to the three-dimensional problem. It is usual for the model to be parameterised into cuboidal cells of uniform conductivity. For this situation, the sensitivity with respect to the conductivity of the th cell is given by = JET.Edv, (7.9) where F is any component of the electric or magnetic field, and V is the volume of the th cell. Most forward-modeffing programs produce the values of the electric field, E, at the vertices of the cells. For E in this form, the integral in eq. (7.9) can be evaluated by tn-linearly interpolating the real and imaginary parts of the components of B over the volume of the cell (an extension to three dimensions of the method described in section 6.2.2): E (at, y, z) = c1 xyz + c2 xy + c3 yz + c4 xz + c5 x + c6 y + c7 z + c8, (7.10) where B represents the real or imaginary part of a component of the electric field. The 8 x 8 system of equations relating the values of B at the vertices of the cell to the coefficients c, i = 1,... , 8, was solved for a general cell using an analytic software package (Wolfram Research, Inc. 1992). The adjoint field, E, is also tn-linearly interpolated over the cell in an analogous manner. With both B and Et in the form of eq. (7.10), the integration in eq. (7.9) can be performed analytically for a general cell to give an expression involving CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 89 the coefficients of the tn-linear expansions. These analytic formulae are then used to evaluate the integral in eq. (7.9) for a particular cell from the appropriate nodal values of E and E. 7.6 Comparison of approximate and accurate three-dimensional sensitivities 7.6.1 Homogeneous halfspace A program was written to compute the approximate sensitivities for the three- dimensional problem. The program was tested by comparing the sensitivities it pro duced for a homogeneous conductivity model with those calculated using the brute-force method. The forward-modeffing program “SAMAYA” (Gupta, Raiche & Sugeng, 1989) was used for the brute-force calculation of the sensitivities. A grounded wire source extending from (50, —5, 0) m to (50, 5,0) m, and measurements of the vertical compo nent of H at (0, 750,0) m were considered. The conductivity model was a homogeneous halfspace of 0.01 Sm’. The region extending from x = 400 to 650 m, from y = —100 to 150 m and from z 0 to 250 m was discretised into 5 x 5 x 5 cuboidal cells as shown in Fig. 7.3. The conductivity of each of these 125 cells was perturbed in turn by 1%. The resulting brute-force sensitivities for measurements of the amplitude and phase of H at a frequency of 100 Hz are illustrated in Figs. 7.4a and 7.4b. The sensitivities shown are those for the shaded cells in Fig. 7.3. The sensitivity is plotted as a constant value over the cell to which it refers. Sensitivities were calculated for the homogeneous halfspace using my approximate sensitivity program. The model parametenisation and experimental configuration were the same as described above. The forward-modelled electric field, E, was computed using SAMAYA (Gupta et al., 1989). The adjoint field was computed in a layered halfspace using the method described in Appendix A. Because observations of the vertical component of the magnetic field were being considered, the source for this adjoint field was a vertical magnetic dipole. The resulting sensitivities are shown in Figs. 7.4c and 7.4d. There is Figure 7.3 Geometry for the three-dimensional example considered in sec tions 7.6 and 7.7. Panel (a) is the side view and panel (b) is the plan view. S indicates the y-directed grounded wire source and R indicates the observa tion location. The shaded cells are the ones for which the sensitivities are shown in Fig. 7.4. The numbers correspond to distances in metres. good agreement between these sensitivities and those calculated using the brute-force method. Figure 7.4 (Following page) Sensitivities for the three-dimensional problem described in section 7.6.1. The conductivity model was a homogeneous halfspace of 0.01Sm1 Panels (a) and (c), and the colour-bar on the left of the figure, show log10 IOln IHI/ôlnoi, and panels (b) and (d), and the colour-bar on the right of the figure, show log10i8q/ôln. HI and 4’ are the amplitude and phase, respectively, of the vertical component of the H-field. The phase is measured in radians. The sensitivities shown are those for the shaded cells in Fig. 7.3. Panels (a) and (b) were produced by the brute-force method, and panels (c) and (d) were produced by the approximate-sensitivity program using an adjoint field computed in a layered halfspace. CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d 3d 90 a ____ z 650OS 50 400 b 750 250 Os -100 z x v4 400 R 150 650 a b N — 6. 20 — 6. 35 — 6. 49 — 6. 64 — 6. 78 — 6. 92 — 7. 07 — 7. 21 — 7. 35 — 7. 50 C a 0 LI) LI) 0 — 5. 40 — 5. 45 — 5. 49 — 5. 54 — 55 8 — 5. 62 — 5. 67 — 5. 71 — 5. 75 — 5. 80 y (m ) y (m ) CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 92 7.6.2 Three-dimensional model Brute-force and approximate sensitivities were also calculated for a simple three- dimensional conductivity model. The survey configuration and model parameterisation were the same as those in section 7.6.1. The conductivity model comprised a conductive block of 0.1 Sm’ in a background of 0.01 Sm’. The conductive block coincided with the central square of nine cells in the third plane of cells down from the surface and the central square of nine cells in the fourth plane of cells below the surface. The brute-force sensitivities are shown in Figs. 7.5a and 7.5b, and the approximate sensitivities calcu lated using the approximate adjoint field in a layered halfspace are shown in Figs. 7.5c and 7.5d. There is, in general, good agreement between the brute-force and approximate sensitivities shown in Fig. 7.5, although there are differences between the sensitivities of the phase for the cells outside the conductive block, especially in the second row (see Figs. 7.5b and 7.5d). However, the resolution of the comparison is poor because of the limited number of cells in the model, making it difficult to determine the significance or cause of these differences. 7.7 Computation times for the three-dimensional problem The computation times for the sensitivities calculated in sections 7.6.1 and 7.6.2 are given in Table 7.2. Six observation locations were considered at the single frequency of 100 Hz. The computations were carried out on a Sun Sparcl0 workstation. An estimated Figure 7.5 (Following page) Sensitivities for the three-dimensional problem described in section 7.6.2. The conductivity model comprised a homogeneous background of 0.015m1with a conductive block of 0.15m coinciding with the cells of enhanced sensitivity shown here. Panels (a) and (c), and the colour-bar on the left of the figure, show log10Oln jHI/ölno, and panels (b) and (a), and the colour-bar on the right of the figure, show log10 ô qf/öln oi. The sensitivities shown are those for the shaded cells in Fig. 7.3. Panels (a) and (b) were pro duced by the brute-force method, and panels (c) and (d) were produced by the approximate-sensitivity program using an adjoint field computed in a layered halfspace. a b — 4. 60 — 4. 74 — 4. 87 — 5. 00 — 5. 14 — 5. 27 — 5. 40 — 5. 54 — 5. 67 — 5. 80 C d I — 5. 10 — 5. 44 — 5. 77 — 6. 10 — 6. 44 — 6. 77 — 7. 10 — 7. 44 — 7. 77 — 8. 10 y (m ) y (m ) CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 94 time for the adjoint-equation method is also given. This estimate is equal to the time required to carry out six forward modeffings for this test-case using SAMAYA. Even for this small example the approximate sensitivities are almost two orders of magnitude quicker to compute than the accurate sensitivities. This time difference will increase as the number of observation locations increases for the same reason as that described in section 6.4. Accurate sensitivities: brute-force 260 mm adjoint-equation 12 mm Approximate sensitivities: layered halfspace 10 s Table 7.2 Example computation times for accurate and approximate sensitiv ities for the three-dimensional inverse problem. It is noteworthy that the approximate sensitivities for the three-dimensional problem are quicker to compute than those for the 2.5d problem for the same number of observation locations and frequencies, and for the same number of layers in the model. To understand this, consider the way in which the adjoint field is computed for the two problems (see Appendix A). For the three-dimensional problem, the two-dimensional inverse Fourier transform that gives the spatial dependence of the adjoint field is converted to a one- dimensional Hankel transform using eq. (2.10) of Ward & Hohmann (1988). This is then evaluated using the digital filtering program of Anderson (1979b). For the 2.5d problem, the k and k dependencies are considered separately, as described in section 7.2. The forward-modelled electric field and the adjoint field are combined in the (x, k) domain. The adjoint field is therefore required as a function of x and for a sequence (typically ten) of values of k?,,. To obtain the adjoint field as a function of x for each value of k, an inverse cosine/sine transform has to be carried out. This is done using the digital filtering technique of Newman, Hohmann & Anderson (1986). The time required to evaluate one Hankel transform and one cosine/sine transform is similar. Because of this, and because the evaluation of the Hankel and cosine/sine transforms accounts for the majority of the CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 95 computation times for the two methods, the sensitivities for the 2.5d problem are roughly ten times slower to compute than those for the three-dimensional problem. The number of cells in the example considered in sections 7.6.1 and 7.6.2 is tiny compared to the number that would be required in a three-dimensional inversion of a re alistic data-set. The number of cells required would be of the order of iOO. However, the number of cells used to produce the sensitivities shown in Figs. 7.4 and 7.5 was the largest for which the brute-force sensitivities could feasibly be calculated. A parameterisation involving 8 x 8 x 8 cells was investigated, but a single forward modeffing for this number of cells required 75 minutes. Carrying out 513 such forward modellings to calculate the sensitivities was considered impracticable. In addition to computation time, the amount of memory required by three-dimensional forward-modeffing programs severely constrains the number of cells that can be considered. For example, the model containing 8 cells required 70 Mbytes of memory. This is approaching the maximum amount of memory typically available on today’s computers. 7.8 Conclusions Programs were written to compute approximate sensitivities for the 2.5d and three- dimensional inverse problems. For the 2.5d problem, the approximate sensitivities for a simple conductivity model were compared with the accurate sensitivities. The agree ment was found to be better than that for the two-dimensional magnetotelluric sensitiv ities discussed in Chapter 6. The approximate sensitivities for the 2.5d problem should therefore enable an iterative, linearised inversion procedure to converge to an acceptable model. Also, for the small example considered in this chapter, the computation of the approximate sensitivities was quicker by a factor of five than the most efficient method of calculating the accurate sensitivities. This factor will increase further as the size of the problem increases, just as for the magnetotelluric sensitivities discussed in Chapter 6. For the three-dimensional problem, the approximate sensitivities were computed for a simple three-dimensional conductivity structure and compared to the brute-force CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 96 sensitivities. There was good agreement between the approximate and brute-force sensi tivities. However, the number of model parameters for which the brute-force sensitivities could be computed was severely restricted by the computation time required by the three-dimensional forward-modelling program. The comparison of the sensitivities was therefore of a much lower resolution than for the two-dimensional magnetotelluric and 2.5d problems. Although the example for the three-dimensional problem was severely restricted in size, it did show that the approximate sensitivities for this problem are considerably quicker to compute than the accurate sensitivities. The difference was roughly two orders of magnitude for the small example presented in this chapter. This time difference will increase significantly as the number of observation locations increases. Chapter 8 Summary The multi-dimensional inversion of geophysical electromagnetic data is too computa tionally demanding to be carried out on a routine basis, even using present-day computer technology. The aim of this thesis was to contribute a means of accelerating the multi dimensional inversion process. In Chapter 2, I described an iterative, linearised inversion procedure for inverting time-domain electromagnetic measurements for a one-dimensional conductivity model. This procedure is typical of the inversion strategy that one would ultimately like to apply to the multi-dimensional inversion of electromagnetic data. The most computa tionally intensive components of this inversion procedure are the forward modelling, the generation of the Jacobian matrix of sensitivities, and the solution of the linear system of equations. In this thesis, I concentrated on developing a quick, approximate method for calculating the sensitivities. The approximate form of sensitivity that I presented is based on the adjoint-equation method. This method is the most efficient way of accurately calculating the sensitivities when the inversion strategy is based on an under-determined problem. As described in Chapter 4, the vast majority of the computation time required by the adjoint-equation method is due to calculating the adjoint field, E, in the multi-dimensional conductivity model for dipole sources at every observation location. The approximate sensitivities are obtained by replacing the true adjoint field by an approximate adjoint field that is much quicker to compute. Three possible forms of the approximate adjoint field were investigated: the electric field computed in a homogeneous halfspace for the appropriate dipole source, the electric field computed in a horizontally-layered halfspace, and that computed using the Born approximation. 97 CHAPTER 8: SUMMARY 98 In Chapter 5, I investigated the three different forms of approximate adjoint field in the context of the purely two-dimensional problem. It was shown that all three forms replicated the dominant behaviour of the true adjoint field. The approximate adjoint field computed in a layered halfspace matched the features in the true adjoint field due to one-dimensional structures in the conductivity model. The approximate adjoint fields computed in the homogeneous and layered halfspaces were considerably quicker to com pute than the true adjoint field. The Born approximation of the adjoint field, however, was slower. Although the computation time of the sensitivities using the Born approxi mation could be made smaller than that for the accurate sensitivities, the computation times for the sensitivities calculated using the approximate adjoint fields in a homoge neous or layered halfspace were even quicker. The approximate adjoint field in the layered halfspace is therefore the most useful of the three possible forms investigated here. In Chapter 6, I investigated the approximate sensitivities for the two-dimensional magnetotelluric problem. I compared approximate and accurate sensitivities for a simple conductivity model. The level of agreement was good, and the approximate sensitivities using the approximate adjoint field in either a homogeneous or layered halfspace were considerably quicker to compute than the accurate sensitivities. The approximate sen sitivities were used within an iterative, linearised procedure to invert both a synthetic data-set and a field data-set. The inversion strategy was the same as that used for the one-dimensional inversion of time-domain electromagnetic data in Chapter 2, except that a subspace methodology was used to solve the linear system of equations. The approxi mate sensitivities performed successfully in the inversion of the synthetic data-set. The field data-set was a sub-set of the COPROD2 data-set which has proved challenging for the many inversion procedures that have been applied to it. The inversion using the approximate sensitivities did not achieve the final target misfit. However, the misfit that was achieved was close to the target value and the major features in the data were re produced. The final model also contained the multiple mid-crustal conductors that other investigations have produced. The inversion using the approximate sensitivities therefore performed as well as any of the other inversion procedures applied to the COPROD2 data. CHAPTER 8: SUMMARY 99 To conclude the work for this thesis, approximate sensitivities were computed for the 2.5d problem and three-dimensional controlled-source problem. These situations, especially the latter, are much more common in geophysics than the special case of two-dimensional magnetotellurics. These situations are also significantly more time- consuming. The approximate sensitivities for the 2.5d and three-dimensional problems were computed for simple conductivity models and compared to the accurate sensitivities. The agreement for the 2.5d problem was found to be better than for the two-dimensional magnetotelluric problem. The agreement for the three-dimensional problem was also found to be good, although the number of model parameters used for the comparison was severely restricted by the computation time required by the three-dimensional forward- modelling program. Since the level of agreement between the approximate and accurate sensitivities was found to be similar, if not better, than that for the two-dimensional magnetotelluric problem, it was assumed that these approximate sensitivities would per form just as well in an iterative, linearised inversion procedure. The more important result from the comparisons of the approximate and accurate sensitivities was the dif ference in the computation times. This difference was considerable, especially for the three-dimensional problem where a difference of two orders of magnitude was found, even for the small example that was considered. Also, the relative differences in compu tation times between the approximate and accurate sensitivities for the two-dimensional, 2.5d and three-dimensional problems all increase as the number of observation locations increases. In this thesis, therefore, I have presented an approximate form of sensitivity that is appropriate for the multi-dimensional inversion of any geophysical electromagnetic data set. The approximate sensitivities appear to be sufficiently good approximations of the accurate sensitivities to enable an iterative, linearised inversion procedure to converge to an acceptable model. And, as desired, the approximate sensitivities are considerably quicker to compute than the accurate sensitivities, especially for the three-dimensional controlled-source problem. References Anderson, W.L., 1979a. Programs TRANS_HCLOOP and TRANS_HZWIRE: calculation of the transient horizontal coplanar loop soundings and transient wire-loop soundings, USGS Open-File Report 79-590. Anderson, W.L., 1979b. Numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering, Geophysics, 44, 1287—1305. Asten, M.W., 1987. Full transmitter waveform transient electromagnetic modeling and inversion for soundings over coal measures, Geophysics, 52, 279—288. Bailey, R.C., 1970. Inversion of the geomagnetic induction problem, Proc. Roy. Soc. London, A315, 185—194. Boerner, D.E. & Holladay, J.S., 1990. Approximate Fréchet derivatives in inductive electromagnetic soundings, Geophysics, 55, 1589—1595. Boerner, D.E. & West, G.F., 1989. A spatial and spectral analysis of the electromagnetic sensitivity in a layered earch, Geophys. J. mt., 98, 11—21. Brant, A.A., 1966. Introduction, in Mining Geophysics, Vol. 1, pp. 3—6, Soc. of Expl. Geophys., Tulsa. Brewitt-Taylor, C.R. & Weaver, J.T., 1976. On the finite difference solution of two- dimensional induction problems, Geophys. J. R. astr. Soc., 47, 375—396. Cagniard, L., 1953. Basic theory of the magneto-telluric method of geophysical prospect ing, Geophysics, 18, 605—635. Constable, S.C., Parker, R.L. & Constable, C.G., 1987. Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data, Geo physics, 52, 289—300. deGroot-Hedlin, C. & Constable, S., 1990. Occam’s inversion to generate smooth, two- dimensional models from magnetotelluric data, Geophysics, 55, 1613—1624. Dosso, S.E. & Oldenburg, D.W., 1991. Magnetotelluric appraisal using simulated anneal ing, Geophys. J. mt., 106, 379—385. Ellis, R.G., Farquharson, C.G. & Oldenburg, D.W., 1993. Approximate Inverse Mapping inversion of the COPROD2 data, J. Geomag. Geoelectr., 45, 1001—1012. Ellis, R.G. & Oldenburg, D.W., 1994. The pole-pole 3-D DC-resistivity inverse problem: a conjugate-gradient approach, Geophys. J. mt., 119, 187—194. Everett, M.E. & Edwards, R.N., 1992. Transient marine electromagnetics: the 2.5-D forward problem, Geophys. J. mt., 113, 545—561. 100 REFERENCES 101 Farquharson, C.G. & Oldenburg, D.W., 1993. Inversion of time-domain electromagnetic data for a horizontally layered Earth, Geophys. J. mt., 114, 433—442. Fullager, P.K. & Oldenburg, D.W., 1984. Inversion of horizontal loop electromagnetic frequency soundings, Geophysics, 49, 150—164. Gradshteyn, I.S. & Ryzhik, I.M., 1994. Table of Integrals, Series, and Products, 5th ed., Academic Press, Boston. Gill, P.E., Murray, W. & Wright, M.H., 1981. Practical Optimization, Academic Press, London. Grant, F.S. & West, G.F., 1965. Interpretation Theory in Applied Geophysics, McGraw— Hill, New York. Golub, G.H. & Van Loan, C.F., 1989. Matrix Computations, The John Hopkins Univer sity Press, Baltimore. Gupta, P.K., Raiche, A.P. & Sugeng, F., 1989. Three-dimensional time-domain electro magnetic modeffing using a compact finite-element frequency-stepping method, Geo phys. J. mt., 96, 457—468. Heiland, C.A., 1940. Geophysical Exploration, Prentice—Hall, New York. Hohmann, G.W., 1988. Numerical Modeling for Electromagnetic Methods of Geophysics, in Electromagnetic Methods in Applied Geophysics, Vol. 1, pp. 313—363, ed. Nabighian, M.N., Soc. Expl. Geophys., Tulsa. Hohmann, G.W. & Raiche, A.P., 1988. Inversion of Controlled-Source Electromagnetic Data, in Electromagnetic Methods in Applied Geophysics, Vol. 1, pp. 469—503, ed. Nabighian, M.N., Soc. Expl. Geophys., Tulsa. Ingham, M.R., Jones, A.G. & Honkura, Y., (eds.) 1993. MT-DIW1 Special Section, J. Geomag. Geoelectr., 45, 931—1150. Jones, A.G., 1993. The COPROD2 Dataset: Tectonic Setting, Recorded MT Data, and Comparison of Models, J. Geomag. Geoelectr., 45, 933—955. Jones, F.W., 1974. The Perturbation of Geomagnetic Fields by Two-dimensional and Three-dimensional Conductivity Inhomogeneities, Fag eoph, 112, 793—800. Jupp, D.L.B. & Vozoff, K., 1975. Stable Iterative Methods for the Inversion of Geophys ical Data, Geophys. J. B. astr. Soc., 42, 957—976. Jupp, D.L.B. & Vozoff, K., 1977. Two-dimensional magnetotelluric inversion, Geophys. J. B. astro. Soc., 50, 333—352. Kaufman, A.A. & Keller, G.V., 1983. Frequency and Transient Soundings, Elsevier, Amsterdam. Lanczos, C., 1961. Linear Differential Operators, Van Nostrand, London. REFERENCES 102 Lee, K.H. & Morrison, H.F., 1984. A solution for TM-mode plane waves incident on a two-dimensional inhomogeneity, Lawrence Berkeley Lab., Rep. LBL-1 7857. Livelybrooks, D., 1993. Program 3Dfeem: a multidimensional electromagnetic finite element model, Geophys. J. mt., 114, 443—458. M8Donald, K.J. & Agarwal, A.K., 1994. Two and Three dimensional induction mod elling on massively parallel supercomputers, presented at the J2th Workshop on EM Induction in the Earth, Brest, France. McGillivray, P.R. & Oldenburg, D.W., 1990. Methods for calculating Fréchet derivatives and sensitivities for the non-linear inverse problem: a comparative study, Geophys. Prosp., 38, 499—524. McGiffivray, P.R., Oldenburg, D.W., Ellis, R.G. & Habashy, T.M., 1994. Calculation of sensitivities for the frequency domain electromagnetic problem, Geophys. J. mt., 116, 1—4. Mackie, R.L. & Madden, T.R., 1993. Three-dimensional magnetotelluric inversion using conjugate gradients, Geophys. J. mt., 115, 215—229. Mackie, R.L., Madden, T.R. & Wannamaker, P.E., 1993. Three-dimensional magne totelluric modeling using difference equations - Theory and comparisons to integral equation solutions, Geophysics, 58, 215—226. Madden, T.R., 1972. Transmission systems and network analogies to geophysical forward and inverse problems, Tech. Rep. 72—3, Dept. of Earth and Planetary Sciences, MIT. Madden, T.R., 1990. Inversion of low frequency electromagnetic data, in Oceanographic and Geophysical Tomography, pp. 379—408, eds. Desaubies, Y., Tarantola, A. & Vinn Justin, J., Elsevier, Amsterdam. Madden, T.R. & Swift, C.M., 1969. Magnetotelluric studies of the electrical conductivity structure of the crust and upper mantle, in The earth’s crust and upper mantle, Geo phys. Monogr. 13, pp. 469—479, ed. Hart, P.3., Am. Geophys. Union, Washington, D.C. Menke, W., 1984. Geophysical Data Analysis: Discrete Inverse Theory, Academic Press, Orlando. Nabighian, M.N., 1988. (ed.) Electromagnetic Methods in Applied Geophysics, Vol. 2, Soc. Expl. Geophys., Tulsa. Nabighian, M.N. & Macnae, J.C., 1991. Time-domain electromagnetic prospecting meth ods, in Electromagnetic Methods in Applied Geophysics, Vol 2., pp. 427—479, ed. Nabighian, M.N., Soc. Expl. Geophys., Tulsa. Newman, G.A. & Alumbaugh, D.L., 1995. Frequency-domain modeling of airborne elec tromagnetic resonses using staggered finite differences, Geophysical Prospecting (ac cepted). REFERENCES 103 Newman, G.A., Hohmann, G.W. & Anderson, W.L., 1986. Transient electromagnetic response of a three-dimensional body in a layered Earth, Geophysics, 51, 1608—1627. Oldenburg, D.W., 1983. Funnel functions in linear and non-linear appraisal, J. geophys. Res., 88, 7387—7398. Oldenburg, D.W. & Ellis, R.G., 1993. Efficient inversion of magnetotelluric data in two dimensions, Phys. Earth Planet. mt., 81, 177—200. Oldenburg, D.W. & Li, Y., 1994. Subspace linear inverse method, Inverse Problems, 10, 915—935. Oldenburg, D.W., McGifflvray, P.R. & Ellis, R.G., 1993. Generalised subspace methods for large scale inverse problems, Geophys. J. mt., 114, 12—20. Oristaglio, M.L. & Worthington, M.H., 1980. Inversion of surface and borehole electro magnetic data for two-dimensional electrical conductivity models, Geophys. Prosp., 28, 633—657. Parker, R.L., 1977. Understanding inverse theory, Ann. Rev. Earth planet. Sci., 5, 36—64. Parker, R.L., 1980. The inverse problem of electromagnetic induction: existence and construction of solutions based on incomplete data, J. geophys. Res., 85, 4421—4428. Parker, R.L., 1994. Geophysical Inverse Theory, Princeton University Press, Princeton. Poddar, M., 1982. A rectangular loop source of current on a two-layered Earth, Geophys. Prosp., 30, 101—114. Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., 1992. Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed., Cambridge University Press, Cambridge. Roach, G.F., 1982. Green’s Functions, Cambridge University Press, Cambridge. Rodi, W.L., 1976. A Technique for Improving the Accuracy of Finite Element Solutions for Magnetotelluric Data, Geophys. J. R. astr. Soc., 44, 483—506. Ryu, J., Morrison, H.P. & Ward, S.II., 1970. Electromagnetic Fields About a Loop Source of Current, Geophysics, 35, 862—896. Schultz, A., Kurtz, R.D., Chave, A.D. & Jones, A.G., 1993. Conductivity discontinuities in the upper mantle beneath a stable craton, G.R.L., 20, 2941—2944. Scriba, H., 1974. A Numerical Method to Calculate the Electromagnetic Field of a Horizontal Current Dipole, Pageoph, 112, 801—809. Smith, J.T. & Booker, J.R., 1988. Magnetotelluric inversion for minimum structure, Geophysics, 53, 1565—1576. REFERENCES 104 Smith, J.T. & Booker, J.R., 1991. Rapid Inversion of Two- and Three-Dimensional Magnetotelluric Data, J. geophys. Res., 96, 3905—3922. Telford, W.M., Geldart, L.P. & Sheriff, R.E., 1990. Applied Geophysics, 2’ ed., Cam bridge University Press, Cambridge. Torres-Verdin, C. & Habashy, T.M., 1993. Cross-Well Electromagnetic Tomography, Expanded Abstract, 3 International Congress of the Brazilian Geophysical Society, Rio de Janeiro. Unsworth, M.J. & Oldenburg, D.W., 1995. Subspace inversion of electromagnetic data - application to mid-ocean ridge exploration, Geophys. J. mt. (submitted). Unsworth, M.J., Travis, B.J. & Chave, A.D., 1993. Electromagnetic induction by a finite electric dipole source over a 2-D earth, Geophysics, 58, 198—214. Ward, S.H., 1990. (ed.) Geotechnical and Environmental Geophysics, Vols. 2 & 3, Soc. Expl. Geophys., Tulsa. Ward, S.H. & Hohmann, G.W., 1988. Electromagnetic theory for geophysical appli cations, in Electromagnetic Methods in Applied Geophysics, Vol. 1, pp. 131—311, ed. Nabighian, M.N., Soc. Expl. Geophys., Tulsa. Weaver, J.T., 1994. Mathematical Methods for Geo-Electromagnetic Induction, John Wiley, New York. Weidelt, P., 1972. The inverse problem of geomagnetic induction, Z. Geophys., 38, 257— 289. Weidelt, P., 1975. Inversion of two-dimensional conductivity structures, Phys. Earth Planet. mt., 10, 282—291. Whittall, K.P. & Oldenburg, D.W., 1992. Inversion of Magnetotelluric Data for a One- Dimensional Conductivity, Geophysical Monograph Series, 5, Soc. Expi. Geophys., Tulsa. Wiggins, R.A., 1972. The general linear inverse problem: Implication of surface waves and free oscillations for Earth structure, Rev. Geophys. Space Phys., 10, 251—285. Wolfram Research, Inc., 1992. Mathemaiica, version 2.1. Wu, F.T., 1968. The inverse problem of magnetotelluric sounding, Geophysics, 33, 972— 979. Appendices Appendix A: Computation of the Electric and Magnetic Fields in a Layered Halfspace due to a Dipole Source The electric or magnetic field generated in a layered halfspace by an electric or magnetic dipole source is required throughout this thesis: in Chapter 2 the vertical component of the magnetic field is required on the surface of a layered halfspace for a horizontal electric dipole source, and in Chapters 5 to 7 the approximate form of adjoint field that is found to be the most useful is the electric field generated in a layered halfspace. In this appendix, I outline the method that was used to calculate the electric and magnetic fields generated in a horizontally-layered halfspace by an electric or magnetic dipole source. The layered haifspace in which the electric or magnetic field is to be calculated is made up of uniform layers of constant conductivity as shown in Fig. Ad. The TE and TM mode Schelkunoff potentials of Ward & Hohmann (1988) are used: A = Ae, (A.1) F = Fe, (A.2) where ê is the unit vector in the z-direction. In the j” layer, A and F satisfy the following ordinary differential equations (section 4, Ward & Hohmann 1988): ( — ) A = o,( - = o, where the tilde represents the two-dimensional Fourier transform, = k + k — i€W2 + iw,uoj, w is the angular frequency, p and e are the magnetic permeability and electrical 105 y xcp 0• Figure A.1 Notation and coordinate system for the horizontally layered con ductivity model of the Earth used in this thesis. is the depth to the bottom of the th layer, and and are the conductivity and thickness, respectively, of the th layer. S represents the dipole source. permittivity, and is the conductivity of the th layer. The definition of the two- dimensional Fourier transform used by Ward & Hohmann is also used here: and aDO aDO (k,k,z) = / / F(x,y,z) e_i k) dxdy, (A.5) i-CO i-DO DO DO F(x,y,z) = U1 APPENDIX A: EM FIELDS IN A LAYERED HALFSPAGE 106 S z=-h o= 0 z=0 z’ z2 zi-’ ti zi ZN 2 UN1 tN-1 ZN1 P(k,k,z) e”’ d1cdk. (A.6) APPENDIX A: EM FIELDS IN A LAYERED HALFSPACE 107 The solutions for and are A(k,k,,z,w) = C(k,k,w) e z_1) + D(k,k,w) e_uj_zj_1), (A.7) = U(k,k,w) euj z_i) + Vj(k,k,w) e_Ui_Zj_). (A.8) In the region of free space above the layered halfspace, A and F satisfy A0(k,k,z,w) = Co(kz,ky,w)eu0z + Do(k,k,w)e_u0z, (A.9) Po(k,k,z,w) = U0(kz,ky,w)eu0z + V(kz,ky,w)e_u0z, (A.10) where ug = k + k. Firstly, consider the potential F. From Ward & Hohmann eqs. (1.152) and (1.153), the boundary conditions on F at z = are k, z=z_i,w) = k, z=z_i,w), (A.11) oP. ____ a2’ k, z=z_1,w) = ‘ (ku,, k, z_—z_1,w). (A.12) Substituting eq. (A.8) into the above boundary conditions gives + V = Ui_ie_1ti_1 + j_ieUj_huj_1, (A.13) uU — uVj = uj_iUi_ieui_lti1 — u_iVj_ieui_1ti_1 (A.14) = — is the thickness of the J’’ layer. These two boundary conditions can be combined in a matrix equation: (1 1 (U = ( e1ti1 e_Ui_lti_1 (A 15 u —Ui) ) U_i&L_1t_1 _u_ei1i1) a_i) This can be rewritten as (v’) = e’ M (), (A.16) APPENDIX A: EM FIELDS IN A LAYERED HALFSPA GE 108 where / (1 + ___)e_2u_1ti_1 (1 — M — (A 17 (‘—ti) (l+L) In eq. (A.16), the exponential term exp(u1t_has been factored out of the matrixM so that the propagation through the layers carried out below remains stable. Applying the boundary conditions at the surface of the halfspace and using eq. (A.lO) for the potential F in the free space above the layered halfspace leads to () = where /(1i!a’ ‘i—-’ _I\ u01 \ u0’1 ‘Al - (1-) (i+’)) uo Uo Eqs. (A.18) and (A.16) can be used to propagate the boundary conditions through the layers to produce an expression relating the coefficients of the potential P in the free space above the halfspace to those in the basement halfspace: N+1 N N+1(E.°) = () ex(>uiti) ]JM, (E11). (A.20) There can be no upward-decaying solution for F in the basement halfspace, so UN+l = 0. It is assumed, at the moment, that the source is at a height z = —h above the layered halfspace. If this source is a unit x-directed electric dipole then, from eq. (4.137) of Ward & Hohmann, v — _____ —uoh A2lo — 2k+k . ( . ) For a unit x-directed magnetic dipole, v—-- k —uoh2 (A.22) APPENDIX A: EM FIELDS IN A LAYERED HALFSPACE 109 using eq. (4.106) of Ward & Hohmann. Eq. (A.20) is therefore a pair of simultaneous equations in the two unknowns U0 and VN+l. Once U0 is known, eq. (A.18) can be used to obtain U1 and V1. Successive applications of eq. (A.16) then give the coefficients U and Vj, and hence the potential F, in every layer. A modified version of the above approach was used for the potential A. The boundary conditions on A at z = z_1 are, from Ward & Hohmann eqs. (1.182) and (1.183), A(k, k, z=z_1, w) = A1(kr, k, z=z_1, w), (A.23) 1 oA. 1 oA. = (k,,k,z=z_i,w). (A.24) Because of the awkwardness, numerically, of the second boundary condition at the surface of the layered halfspace (above which o = 0), the source is assumed to lie within layer 1 (at z = h). The final result for the source on the surface of the halfspace is obtained by letting the source approach the surface from below. In layer 1, therefore, eq. (A.7) is modified to contain a term representing the particular solution: A1 = G eU1Z + D1 e_U1Z + A e_U1_l. (A.25) Extending the analysis of Ward & Hohmann to determine explicit expressions for their quantity A, the appropriate expression for A eq. (A.25) for a unit c-directed electric dipole is = 2 k+k’ (A.26) and A — [LW — iW/2U1 ik A 27p — — 2u k+k’ for a unit x-directed magnetic dipole. The boundary conditions can now be propa gated through the layers in a similar manner to the potential F to obtain two simul taneous equations in the two unknowns C0 and DN+l (D0 is zero since there is now no downward-decaying part of the solution in the free space above the halfspace). To APPENDIX A: EM FIELDS IN A LAYERED HALFSPAGE 110 obtain the potentials for y-directed electric and magnetic dipoles, the transformation (x,y) —* (y,—x) [ —* (k,_k)] can be used. Once the potentials A and F are known, the components of the electric and magnetic fields can be calculated throughout the layered halfspace using eqs. (1.129) and (1.130) of Ward & Hohmann: E = a2A — (A.28)z o+iew 9xOz = 1 82A + , (A.29) cr + iew Oyôz 72 = . (- + k ) A, (A.30) c7+iEW \UZ J and 8A 1 L92F H = + -—-- , (A.31) Oy iwu 9x8z OA 1 02F (A.32) OX iW[L oyaz 72 H = (-- + k ) F, (A.33) ZWI \OZ / where k2 = — iwuo. The above equations are in fact transformed to the (k,k) domain and used to obtain B and H from A and F, and then the transformation back to the (x,y) domain that is appropriate for the particular problem is carried out. For the purely two-dimensional case discussed in Chapters 5 and 6, the electric field for the two-dimensional dipole sources is obtained by setting the along-strike wavenumber, k, equal to zero before carrying out the inverse cosine/sine transform to recover the x dependence of the electric field. For the 2.5d case discussed in Chapter 7, the inverse cosine/sine transform with respect to the across-strike wavenumber, k, is performed at different values of the along-strike wavenumber, k. The resulting (z, k )-dependence of the electric field is then used in the calculation of the sensitivities. Finally, for the three dimensional case also described in Chapter 7, the two-dimensional Fourier transform with respect to k and k is converted to a Hankel transform using eq. (2.10) of Ward & APPENDIX B: PURELY 2d SENSITIVITIES 111 Hohmann, and this Hankel transform performed to obtain the full spatial dependence of the electric field. The numerical evaluation of the cosine/sine transforms was carried out using the method of Newman, Hohmann & Anderson (1986), and the evaluation of the Hankel transform was carried out using the method of Anderson (1979b). Appendix B: Adjoint-equation method for the purely two-dimensional inverse problem For a purely two-dimensional problem, that is, one in which both the conductivity model and the source are invariant in the strike direction, the adjoint field required to calculate the sensitivities becomes the electric field due to a two-dimensional dipole source. This can be seen as follows. Consider the particular form of eq. (4.14) for a given measurement and assume that the model and ID are invariant in the strike direction which is chosen here as the y-direction. One obtains = f J J E(z,y,z) .E(z,z) (x,z) dxdydz (B.1)9o.j —00 —00 = jj(x,z) E(x,z).{f00Ef(x,Y,z) dy} dxdz. (B.2) where F represents the component of interest of the electric or magnetic field. With out loss of generality, the adjoint field can be expressed as a two-dimensional Fourier transform: E(x,y,z) = J J Et(k,k,z) dkdk. (B.3) -00 -00 APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 112 If this expression for the adjoint electric field is substituted into eq. (B.2), then the integration with respect to y can be carried out, since = S(k). (B.4) This reduces the adjoint field to only its zero along-strike-wavenumber component: = jfbj(xz) E(a,z).{_j dk} dxdz. (B.5) The term within the braces is the electric field due to a two-dimensional source, that is, one that is invariant in the strike direction. Eq. (B.5) can therefore be rewritten as = j E(x, z) . E(x, z) (x, z) ds, (B.6) where the adjoint field, Et, is now due to a two-dimensional dipole source at the observa tion location (x0,z0). This is the appropriate expression for the sensitivities for a purely two-dimensional inverse problem. Appendix C: Electric fields due to two-dimensional dipole sources on a homogeneous halfspace The first form of approximate adjoint field presented in Chapter 5 for the purely two dimensional inverse problem is the electric field generated in a homogeneous halfspace by a two-dimensional dipole source. To calculate the corresponding approximate sensitivities in Chapter 6 for the two-dimensional magnetotelluric inverse problem, this approximate adjoint field is required for both electric and magnetic dipole sources oriented in both the a- and y-directions on the surface of the halfspace. APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 113 Electric dipole source To develop an expression for the electric field generated in a homogeneous halfspace by an x- or y- directed two-dimensional unit electric dipole, initially follow the derivation of Kaufman & Keller (1983) for the electric field generated by a point (in three-dimensional space) electric dipole. Consider an electric dipole oriented in the y-direction and suppose that both the dipole and the point at which the electric field is to be calculated are situated within the conductive halfspace (see Fig. C.1). Assume that z is positive downwards. Kaufman & Keller make use of the vector potential A: A = (O,A,A) (C.1) such that E = iwA + V(V . A). (C.2) w is the angular frequency, and [I and cr are the magnetic permeability and conductivity of the halfspace. The two non-zero components of the vector potential are: = _ j {e’l + J0(Ap) dk (C.3) = j u2+A e_UJ1(p) dA, (C.4) where u2 = )2 + iw,ucr and p2 = x2 + y2. I have followed the convention that Ward & Hohmann (1988) use for the time- to frequency-domain Fourier transform. This accounts for the plus sign in the above expression for u. h is the z-coordinate of the dipole source. Using eq. (2.10) of Ward & Hohmann (1988) to convert the above Hankel transforms to two-dimensional Fourier transforms gives 1 °° °°1 — = J J — { Iz-hI + + e_u} dk dk (C 5) = — L L A(u+ A) e_U ik e ky) dk dk, (C.6) APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 114 o=O z=O S z=h 0 0• Figure C.1 The coordinate system and geometry for computing the electric field induced in a homogeneous halfspace by a two-dimensional dipole source, S. where now A2 = k + k. Substituting these two expressions into eq. (C.2) gives the three components of the electric field resulting from the y-directed finite electric dipole: 1 J kk (e’ u A 82 _ u+A _u(z+h)) + e + k1cu uA(u + A)e} dk dk, (C.7) 1 J J°° f (iwp — k (e_uIz_ + U —82 — u u) u+A 2ku — A(u + A)e} dkdk, (C.8) and 1 fofoIik / = — (u sgn(z — h) e_uIz_hI + u(u — A) e_u(z+h)) _cru u+A 2ik AV U + ( + A) e (z+h)} dk dk. (C.9) APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 115 Appendix B shows that the adjoint field required to calculate the sensitivities for a purely two-dimensional inverse problem can be obtained from the electric field for a finite dipole source by considering only the zero along-strike-wavenumber component. If this is done in the above equations the following expressions are obtained for the electric field generated in a homogeneous halfspace by a two-dimensional y-directed electric dipole source: = 0, (0.10) B = — f e_UZ dlc , (0.11)2r J_u+kj = 0. (0.12) where it has been assumed that the source is at the surface of the halfspace (h = 0) and now u2 = k + iwo. From the above equations it is clear that for this type of source there is only an along-strike component of the electric field. This is exactly what is required when calculating the sensitivities for the E-polarisation mode of the two- dimensional magnetotelluric inverse problem. Note that eq. (0.11) agrees, as it should, with eq. (4.208) of Ward & Hohmann for the electric field due to an infinite line current. To obtain expressions for the electric field produced by an x-directed two-dimensional electric dipole first consider the electric field produced by a finite x-directed electric dipole. This field is given by eqs. (0.7) to (0.9) after rotation of the coordinate axes corresponding to the transformation (x,y) —* (y, —x). This implies that (1c,k) —÷ (ks,, —kr) and (Br, B, E) —* (Er, —Er, Es). The electric field due to the two-dimensional dipole source is then obtained by considering only the components of the equations corresponding to zero along-strike wavenumber: ___jU(e’ + e_u) dk, (0.13) = 0, (0.14) = f k (sgn(z — h) eZ + e_u(z+h)) dk. (0.15) APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 116 The electric field for an x-directed two-dimensional electric dipole is therefore restricted to lie in the plane perpendicular to the strike direction. The integrals in the above equations can be evaluated (Lee & Morrison 1984) to give analytic expressions for the field components: = —- j- K(ik’r) + z) K1(ikr), (C.16) = K0(ikr) + --K1(ikr), (C.17) where k2 = —i wjio- and r2 = x2 + z2 for the electric dipole on the surface of the halfspace. K0 and K1 are the zeroth and first order Modified Bessel functions of the second kind. Magnetic dipole source To develop expressions for the electric field induced in a homogeneous halfspace by a two-dimensional magnetic dipole source initially follow the analysis of Ward & Hohmann (1988). Consider the Schelkunoff potentials A and F such that E — 1 ÔAZ OFZ C18 OxOz 182A OFz z C19 yuOyOz Ox’ = + k2) A, (C.20) where k2 = —iwiu. For a unit x-directed magnetic dipole = j°° / (e’ — e_u) ei(kx+kv) dk dk, (C.21) and F = i_: L (€_ulz_h + e_u) dk dk. (C.22) APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 117 Here, as before, u2 = k + k + iw,uo and )2 = k + k. Substituting these expressions into eqs. (0.18) to (0.20) gives J f (sgnz — h) — 1]e_uIz_ + 2A e u(z+h))8712 — j_ U + x dk dk, (0.23) / = j j { (sgn(z — + e_u) k2 U — e_u) } dk dk, (0.24)+ (eu1zh1 + + and j°° °° ik = 8 2 j ___(e_ul7_I — e_u(z+h)) dk dk. (0.25)71 j_ - U To obtain expressions for the electric field generated by a two-dimensional z-directed magnetic dipole on the surface of the halfspace consider the reduced forms of the above equations for k = 0: (0.26) iwI f U e_UZ dk, (0.27) = 2irju+k = 0. (0.28) For a y-directed two-dimensional magnetic dipole use the transformation (x, y) —* (y, —z) in eqs. (0.23) to (0.25) before considering the reduced form of these equations. This gives = f e — u z e z ai, (0.29)271 = 0, (0.30) = 0. (0.31) APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 118 Using eq. (3.914) of Gradshteyn & Ryzhik (1994), eq. (C.29) can be reduced to wpkz = — —K1(zkr) (C.32) where k2 = —iw,uo and r2 2 + The electric fields described above were computed by evaluation of the Bessel func tions or, if the expression for the electric field could not be reduced to one involving Bessel functions, by evaluation of the Fourier transform using the digital filtering code of Newman, Hohmann & Anderson (1986). Appendix D: Sensitivities for the magnetotelluric apparent resistivity and phase In the magnetotelluric inverse problem, the data are usually considered to be values of the apparent resistivity, pa, and phase, q, where p(w) = --- and (w) = phase (i). (D.1) E and H represent orthogonal horizontal components of the electric and magnetic fields. w is angular frequency and u is magnetic permeability. Since it is the apparent resistivity and phase that are the data in the inverse problem, the sensitivities are required for these data rather than for the electric and magnetic fields. However, the sensitivities for the apparent resistivity and phase can be obtained from those for the electric and magnetic fields as follows. Differentiating the apparent resistivity in eq. (D.1) with respect to the model parameter o gives (D28o wi H Ocr \jHj) APPENDIX D: SENSITIVITIES FOR MT 119 — 2 E fi OEI IE OH D 3 — w H UH H2 ( Note that all the quantities in the above expression, the electric and magnetic fields as well as the sensitivities, are to be evaluated at the observation location. Since the electric field in the frequency domain is a complex quantity it can be written as E = Ee4’E. (D.4) Treating both Ej and 4E as functions of the model parameters, differentiating eq. (D.4) with respect to o gives - = E —- (i) + e (D.5) Ocr ôcr Ocrj = iEea + EOjE{ (D.6)8cr E &r 1 ÔE ___ 1 ÔIEI = ä=ô (D.7) Equating the real and imaginary parts of eq. (D.7) gives ___ = () (D.8) and ___ = m (). (D.9) This argument can also be applied to the magnetic field resulting in ôIH = Hj Re (D.1O) 9uj \H 9c7jJ and = m (). (D.11) APPENDIX D: SENSITIVITIES FOR MT 120 Eqs. (D.8) and (D.10) can be substituted into eq. (D.3) to give a final expression for the sensitivity for the apparent resistivity: — - D12 — w H UHI E Ou) HI H o)j . = 2Pa {e () - e £Z) }. (D.13) This expression can be used to calculate the sensitivity for the apparent resistivity since Pa is known at the observation location, and OE/Ou and OH/ôo can be calculated using eq. (4.14). To determine the final expression for the sensitivity for the phase consider the ratio B — IEIei4’ — (D14HIHIe’H H By definition, = — (D.15) and so — __ — (D16) oo.j — oo = m () - m (E (D.17) using eqs. (D.9) and (D.11). The sensitivity for the phase can therefore be calculated using eq. (D.17) and eq. (4.14).
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
United States | 5 | 3 |
United Kingdom | 2 | 1 |
China | 2 | 0 |
Germany | 1 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 3 | 1 |
Mountain View | 3 | 3 |
Beijing | 2 | 0 |
Sunnyvale | 1 | 0 |
Ashburn | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Share to: