UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Inversion of em data to recover 1-d conductivity and geometric survey parameter Walker, Sean Eugene 1999

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1999-0390.pdf [ 10.36MB ]
Metadata
JSON: 831-1.0053001.json
JSON-LD: 831-1.0053001-ld.json
RDF/XML (Pretty): 831-1.0053001-rdf.xml
RDF/JSON: 831-1.0053001-rdf.json
Turtle: 831-1.0053001-turtle.txt
N-Triples: 831-1.0053001-rdf-ntriples.txt
Original Record: 831-1.0053001-source.json
Full Text
831-1.0053001-fulltext.txt
Citation
831-1.0053001.ris

Full Text

INVERSION OF E M DATA TO RECOVER 1-D CONDUCTIVITY AND A GEOMETRIC SURVEY PARAMETER By Sean Eugene Walker B. Sc. (Honours), Geology & Physics, McMaster University, 1996 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E STUDIES D E P A R T M E N T O F E A R T H A N D O C E A N S C I E N C E S ( G E O P H Y S I C S ) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A June, 1999 © Sean Eugene Walker, 1999 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. Department of The University of British Columbia Vancouver, Canada Date 18 , 1 7 7 9 DE-6 (2/88) Abstract The presence of geometrical survey parameter errors can cause problems when attempting to invert electromagnetic (EM) data. There are two types of data which are of particular interest: airborne EM (AEM) and ground based horizontal loop EM (HLEM). When dealing with AEM data there is a potential for errors in the measurement height. The presence of measurement height errors can result in distortions in the conductivity models recovered via inversion. When dealing with HLEM there is a potential for errors in the coil separation. This can cause the inphase component of the data to be distorted. Distortions such as these can make it impossible for an inversion algorithm to predict the inphase data. Examples of these types of errors can be found in the Mt. Milhgan and Sullivan data sets. The Mt. Milligan data are contaminated with measurement height errors and the Sullivan data are contaminated with coil separation errors. In order to ameliorate the problems associated with geophysical survey parameter errors a regularized inversion methodology is developed through which it is possible to recover both a function and a parameter. This methodology is applied to the 1-D EM inverse problem in order to recover both 1-D conductivity structure and a geometrical survey parameter. The algorithm is tested on synthetic data and is then applied to the field data sets. Another problem which is commonly encountered when inverting geophysical data is the problem of noise estimation. When solving an inverse problem it is necessary to fit the data to the level of noise present in the data. The common practice is to assign noise estimates to the data a priori. However, it is difficult to estimate noise by observation alone and therefore, the assigned errors may be incorrect. Generalized cross validation is ii a statistical method which can be used to estimate the noise level of a given data set. A non-linear inversion methodology which utilizes GCV to estimate the noise level within the data is developed. The methodology is applied to 1-D EM inverse problem. The algorithm is tested on synthetic examples in order to recover 1-D conductivity and also to recover 1-D conductivity as well as a geometrical survey parameter. The strengths and limitations of the algorithm are discussed. iii Table of Contents Abstract ii List of Tables vii List of Figures viii Acknowledgement xiii 1 Introduction 1 1.1 Motivation for the Thesis 2 1.2 Outline of the Thesis 10 2 Background Information 12 2.1 Denning the EM Data 12 2.1.1 Generic Horizontal Coplanar Systems 14 2.1.2 Specific Horizontal Coplanar Systems 18 2.2 Details of the Inverse Problem 24 2.2.1 The Forward Problem 24 2.2.2 The Inverse Problem 29 3 Geometric Survey Parameter Errors 37 3.1 Measurement Height Errors in AEM 40 3.1.1 How do measurement height errors occur? 40 3.1.2 How do measurement height errors affect the data? 42 iv 3.1.3 How do measurement height errors affect the inversion results? . . 44 3.2 Coil Separation Errors in HLEM data 47 3.2.1 How do coil separation errors occur? 47 3.2.2 How do coil separation errors affect the data? 48 3.2.3 How do coil separation errors affect the inversion results? 53 3.3 Parameter Errors in the Field Data Sets 60 3.3.1 Mt. Milligan 60 3.3.2 Sullivan Data 61 4 Recovering 1-D Conductivity and a Geometric Survey Parameter 62 4.1 Changes due to the extra parameter 63 4.1.1 The Model 63 4.1.2 The Data 64 4.1.3 The Sensitivities 65 4.1.4 The Objective Function 67 4.2 Inversion Algorithms 68 4.2.1 Details of algorithm to recover 1-D conductivity with fixed p . . . 69 4.2.2 Defining the Parameter Norm 80 4.2.3 Inversion algorithm to recover 1-D conductivity and p 97 4.3 Synthetic Examples 101 4.3.1 AEM Sounding 101 4.3.2 Trade off between a and h 103 4.3.3 HLEM sounding 107 4.4 Field Examples 108 4.4.1 Mt. Milligan 110 4.4.2 Sullivan 114 v 4.5 Summary 114 5 Application of GCV to 1-D E M Inversion 118 5.1 Estimating Noise using GCV 119 5.1.1 GCV in Linear Inverse Problems 119 5.1.2 GCV in Non-linear Problems 120 5.2 GCV and the 1-D EM inverse problem 122 5.2.1 GCV when p is fixed 122 5.2.2 GCV when p is included in the model 130 5.2.3 Application of GCV to small data sets 133 5.3 Summary 142 6 Summary 143 References 145 vi List of Tables 3.1 Results from the inversion of AEM data with the incorrect heat values. . . 44 3.2 Results from the inversion of small separation HLEM data with the incor-rect s value 54 3.3 Results from the inversion of small separation HLEM data with the incor-rect s value and noise estimated to account for inphase modelling errors. 56 3.4 Results from the inversion of large separation HLEM data with the incor-rect s value and noise estimated to account for inphase modelling errors. 58 4.1 Comparison of results from the AEM inversions with fixed h and variable h.lOb 4.2 Comparison of results from the overburden AEM inversions with fixed h and variable h 107 5.1 Comparison of parameter values from the AEM discrepancy principle and GCV inversions with a fixed value of h 125 5.2 Comparison of results from AEM example discrepancy principle and GCV inversions with h included in the model 132 5.3 Comparison of results from the four frequency AEM example discrepancy principle and GCV inversions with a fixed value of h 137 vii List of Figures 1.1 The inphase and quadrature components of the Sullivan data 3 1.2 Inversion results from the Sullivan data 4 1.3 Observed and predicted data from the inversion of the Sullivan data. . . 5 1.4 Comparison of the inphase and quadrature data from two soundings along Line 1 of the Sullivan data 6 1.5 The inphase and quadrature components of the Mt. Milligan data 7 1.6 Inversion results from the Mt. Milligan data 8 2.1 The four common types of transmitter-receiver geometries used in loop-loop EM surveys 13 2.2 Geometry of the generic horizontal coplanar EM system 14 2.3 Cartoon showing the behavior of a horizontal coplanar EM system in free space and in the presence of a conductive body 15 2.4 Cartoon of a typical airborne EM system 19 2.5 The bucking coil system used in AEM 20 2.6 Cartoon of a typical HLEM system 22 2.7 Discretization of the problem domain for the forward problem 26 2.8 Cartoon of a typical Tikhonov curve 34 3.1 Conductivity models and data from the synthetic examples 39 3.2 Values used to determine measurement height 41 3.3 The two situations that will result in measurement height errors 42 3.4 Synthetic data from AEM soundings with htrue = 24 m, 30 m, and 36 m. 43 V l l l 3.5 Results from the inversion of AEM data with the incorrect h values. . . . 45 3.6 The two situations that will result in coil separation errors 48 3.7 Synthetic data from small separation HLEM soundings with seat = 10 m and Strue = 9 m, 10 m, and 11 m 50 3.8 Synthetic data from large separation HLEM soundings with se3t = 50 m and s^ue = 49 m, 50 m, and 51 m 52 3.9 Results from the inversion of small separation HLEM data with the incor-rect 5 value 55 3.10 Results from the inversion of small separation HLEM data with the incor-rect s value and the noise was estimated to account for inphase modelling errors 57 3.11 Results from the inversion of large separation HLEM data with the incor-rect s value and the noise was estimated to account for inphase modelling errors 59 4.1 Flowchart for the fixed inversion algorithm with a fixed value of p. . . . 72 4.2 Recovered models and predicted data from AEM example fixed fi inversion with a fixed value of h 74 4.3 Convergence curves from AEM example fixed j3 inversion with a fixed value of h 75 4.4 Flowchart for the discrepancy principle inversion algorithm with a fixed value of p 78 4.5 Recovered models and predicted data from AEM example discrepancy principle inversion with a fixed value of h 80 4.6 Convergence curves from AEM example discrepancy principle inversion with a fixed value of h 81 ix 4.7 Results from eleven inversions of the AEM sounding data with fixed values of h ranging from 45 m to 15 m 83 4.8 Results from eleven inversions of the small separation HLEM sounding data, Stme = 11 na, with fixed values of s ranging from 12.5 m to 7.5 m. . 84 4.9 Results from eleven inversions of the large separation HLEM sounding data, Stme = 51 m, with fixed values of s ranging from 52.5 m to 47.5 m. 85 4.10 Values calculated during the line search for j3 during the first iteration of the inversion of the AEM data with ap = oo and 0 88 4.11 Values calculated during the line search for (3 during the first iteration of the inversion of the AEM data with ap = oo, 10 - 6, and 10 - 3 90 4.12 Results from the search for a best fit pref value for the AEM example. . . 92 4.13 Results from the search for a best fit pTef value for the HLEM examples. 93 4.14 Contour plots of (f>d for a range of crTef and p values for the AEM example and both HLEM examples 95 4.15 Flowchart for the fixed {3 inversion algorithm with p included in the model. 98 4.16 Recovered models and predicted data from AEM example fixed /3 inversion with h included in the model 99 4.17 Convergence curves from AEM example fixed (3 inversion with h included in the model 100 4.18 Flowchart for the discrepancy principle inversion algorithm with p included in the model 102 4.19 Recovered models and predicted data from AEM example discrepancy principle inversion with h included in the model 103 4.20 Convergence curves from AEM example discrepancy principle inversion with h included in the model 104 4.21 Comparison of models recovered from the inversion of AEM data with fixed h and variable h 105 4.22 Recovered models and predicted data from small separation HLEM ex-ample discrepancy principle inversion with s included in the model. . . . 108 4.23 Convergence curves from small separation HLEM example discrepancy principle inversion with s included in the model 109 4.24 Recovered models and predicted data from large separation HLEM ex-ample discrepancy principle inversion with s included in the model. . . . 110 4.25 Convergence curves from large separation HLEM example discrepancy principle inversion with s included in the model I l l 4.26 Inversion results from Mt. Milligan Data with h included in the model. . 112 4.27 Recovered 8h superimposed on the fixed h inversion results from Mt. Mil-ligan 113 4.28 Observed and predicted data from the inversion of the Sullivan data with s included in the model 115 4.29 Inversion results from Sullivan Data with s included in the model 116 5.1 Flowchart for the GCV inversion algorithm with a fixed value of p 123 5.2 Comparison of recovered models and predicted data from AEM example discrepancy principle and GCV inversions with a fixed value of h 125 5.3 Convergence curves from AEM example discrepancy principle inversion with a fixed value of h 126 5.4 Convergence curves from AEM example GCV inversion with a fixed value of h 127 5.5 GCV Curves and recovered models from iterations 1, 3, and 6 of the AEM example GCV inversion with a fixed value of h 129 xi 5.6 Flowchart for the GCV inversion algorithm with p included in the model. 131 5.7 Comparison of recovered models and predicted data from AEM example discrepancy principle and GCV inversions with h included in the model. . 132 5.8 Convergence curves from AEM example discrepancy principle inversion with h included in the model 134 5.9 Convergence curves from AEM example GCV inversion with h included in the model 135 5.10 GCV Curves and recovered models from iterations 1,3, and 6 of the AEM example GCV inversion with h included in the model 136 5.11 Comparison of recovered models and predicted data from the four fre-quency AEM example discrepancy principle and GCV inversions with a fixed value of h 138 5.12 Convergence curves from four frequency AEM example discrepancy prin-ciple inversion with a fixed value of h 139 5.13 GCV Curves and recovered models from iterations 1, 2, and 3 of the AEM example GCV inversion with a fixed value of h 141 xii Acknowledgement First of all I would like to thank Doug Oldenburg for opening my eyes to the wonderful world of inversion. The amount that I have learned by being a part of his research group over the past three years is startling. I would also like to thank him for his support throughout my degree. I would be hard pressed to find a more enthusiastic, positive, and understanding person anywhere. I would like to thank Tanya for her love and support. Special thanks go to Len for making up the second half of the 400 lb. office. The research might have gone faster without him, but it wouldn't have been nearly as enjoy-able. Thanks also go to all of my friends in the Geophysics building for providing the opportunity for many not-so-scientific experiences. I must also thank Colin Paton who was omitted from thank yous in a previous work of mine. Thanks for giving me an excuse to go to Edmonton twice a year. xiii Chapter 1 Introduction Both airborne and ground based frequency domain electromagnetic (EM) geophysical methods were initially developed as mineral exploration tools. The early EM systems were rather crude and their applications were limited to locating highly conductive ore bodies within resistive host rocks. Further research and development, coupled with the use of digital technology, resulted in modern EM instruments with greater measurement precision and a wider range of measurement frequencies. These advancements have made it possible to apply EM methods to a wide variety of problems which include the esti-mation of sea ice thickness (Kovacs et ai, 1995), the detection of contaminant plumes (Sauck et al, 1998), and the mapping of hydrogeologic structures (Wynn & Gettings, 1998). The widespread use of geophysical EM methods has generated a need for new data interpretation and processing tools. While traditional curve matching techniques are still used, the importance of EM inversion is now being realized. The development of fast inversion techniques and the affordability of high powered personal computers has helped to make EM inversion a feasible processing option for most practicing geophysicists. The ideal inversion scheme would have the ability to invert multiple lines of EM data in order to recover a 3-D distribution of conductivity. However, fast and reliable 3-D forward modelling codes are still in the development stages. Therefore, most of the available codes are restricted to inverting data from individual soundings to recoverl-D conductivity structures. The recovered 1-D models are then plotted one next to the other 1 Chapter 1. Introduction 2 to produce a pseudo 2- or 3-D conductivity structure. The theoretical details of solving the 1-D EM inverse problem are well established (Fullagar & Oldenburg, 1984 and Zhang & Oldenburg, 1999) and the algorithms have been proven to be successful when applied to field data sets. However, it is always important to be aware of the problems that can arise when attempting to invert field data. 1.1 Motivation for the Thesis The motivation for this thesis comes from problems encountered when attempting to invert two field data sets. The first data set was from a Max-Min ground based EM survey performed over a tailings pond at Cominco's Sullivan Mine in Kimberley, British Columbia. The goal of the survey was to investigate the shallow subsurface conductivity structure around the pond in order to determine the effects of groundwater flow on the tailings (Jones, 1996). The data were collected using the horizontal coplanar configura-tion of the Max-Min system. Measurements of the inphase and quadrature component of the secondary fields were taken at 4 frequencies ranging from 7040 Hz to 56320 Hz. The coil separation was preset to be 5 m for each measurement. I will concentrate on Line 1 which was 350 m long with stations every 5 m. The inphase and quadrature data are shown in Figure 1.1. Upon inspection the Sullivan data appear to defy interpretation. The measurements oscillate from one station to the next and the inphase component is almost entirely negative. The inphase component of the data generated from a horizontal coplanar EM survey over top of a 1-D conductivity structure should asymptote to zero as the measurement frequency decreases however, the Sullivan data do not exhibit this behavior. Therefore, in order to invert the Sullivan data with existing algorithms it is necessary to Chapter 1. Introduction 3 o o 7040 Hz -40 • ° 14080 Hz 0 28160 Hz 56320 Hz 0 50 100 150 200 Easting (m) 250 300 350 Figure 1.1: The inphase and quadrature components of the Sullivan data. 7040 Hz (circles), 14080 Hz (squares), 28160 Hz (diamonds), and 56320 Hz (stars). assign error estimates to the inphase measurements that are large enough to counteract the problems caused by the negative inphase values. This is equivalent to discarding the negative inphase values. When inverting the Sullivan data I assigned standard deviations of 20% and 1% of the primary field to the inphase and quadrature components of the data respectively. The target data misfit for each sounding was equal to 8. The results from inverting each sounding in the Sullivan data are shown in Figure 1.2. The top panel of Figure 1.2 shows that it was possible to achieve the target misfit at most stations. The recovered conductivity models plotted in the bottom panel of Figure 1.2 show a trend within the subsurface from a resistive zone to a more conductive zone as the stations move from east to west along the Une. The survey information provided with the Sullivan data indicated that the eastern end of Line 1 (OE to 40E) lies on top of a bedrock outcrop while the remainder of the Une extends over the taihngs pond itself. Therefore, the recovered model seems to be in agreement with the available a priori information. These results suggest that the SulUvan data have been successfully inverted. However, by plotting the observed data and the data predicted by the recovered model (Figure 1.3) it is clear that I have Chapter 1, Introduction 4 Easting (m) Figure 1.2: Inversion results from the Sullivan data. Top panel: the final data misfit values from each inversion plotted by station location. Bottom panel: recovered 1-D conductivity models from each inversion plotted by station location. only predicted the quadrature component of the data. While in some sense the inversion has been successful I am unsure as to whether I should believe the results because such large errors have been assigned to the data. I am also left to wonder whether the information contained in the inphase component could be important. These problems prompted me to determine the cause of the negative inphase values. A potential cause could have been the presence of magnetizable materials within the tailings however, the negative inphase data is seen at all station locations along Line 1. The data from soundings at two stations are plotted for comparison in Figure 1.4. The sounding in Figure 1.4(a) is from station 5E at the east end of Line 1 on a bedrock outcrop on the edge of the tailings pond. It is fairly safe to assume that this rock will not contain large amounts of magnetite. Figure 1.4(b) shows a sounding from station 150E which is located in the center of the tailings pond. Both soundings show comparable negative shifts of the inphase data. This suggests that these shifts are due to something other than magnetic susceptibility. Chapter 1. Introduction 5 (a) Co 20 0 s N X 1^ o 9=-40. - 2 0 f ° Voo° cp o O O O C e 0o 0 0°°o 0 0 o P W P ^ W ^ c o ^ e P ^ o o 'o°oooc 0 50 100 150 200 Easting 250 300 350 o IP Obs - IP Pre « Q Obs -- QPre as 20 XL a ^ - 4 0 - - 2 0 f ° ° ^ o o .oOo^oP°coooo° o ^ o c o ^ c o o 0 0 0 o„cP 50 100 150 200 250 300 350 Easting (c) N X oo 20 O £ - 4 0 pod 50 100 150 200 250 300 350 Easting (d) £ 20 £ 0 8-20?° V x ^ V * 0 0 o Q . - 4Q 50 100 150 200 Easting 250 300 350 Figure 1.3: Observed and predicted data from the inversion of the Sullivan data. Ob-served inphase (circles), predicted inphase (solid line), observed quadrature (crosses), and predicted quadrature (dashed line), (a) Inphase (IP) and quadrature (Q) 7 kHz, (b) IP and Q 14 kHz, (c) IP and Q 28 kHz, and (d) IP and Q 56 kHz. Chapter 1. Introduction 6 20 10 •o CO "D !-io in CB "I-20 -30 10 G O IP • X X Q X x X X o 00 I 20 10 d o "D !-io 0) (0 | - 2 0 -30 O O I P X X Q X X X p O o O 10* Frequency (Hz) 10' 10 10' Frequency (Hz) 10a Figure 1.4: Comparison of the inphase and quadrature data from two soundings along Line 1 of the Sullivan data, (a) Data from Station 5E and (b) Data from Station 150E. Inphase (circles) and quadrature (crosses). The Sullivan data was collected using a 5 m coil separation which is the smallest coil separation at which it is possible to perform a Max-Min survey. When performing such a survey it is not uncommon for the inphase component of the data to be affected by coil positioning errors (Alumbaugh & Newman, 1997). This suggested the possibility of coil separation errors as the cause of the shifts in the inphase data. Examination of the reference cable used to collect the Sullivan data revealed that it was chained correctly to be 5 m in length however, due to the design of the Max-Min system, when the cable was pulled tight, the coils were in fact 5.5 m apart. If this were to happen during a survey it would result in a coil separation error of 0.5 m. Such an error is negligible when considering a survey using a 100 m coil separation however, at a separation of 5 m it represents a 10% error. The negative inphase measurements in the Sullivan data are therefore, likely the result of coil separation errors. The second set of field data that prompted this research was from a Dighem airborne EM survey performed over the Mt. Milligan deposit in north central British Columbia. Mt. Milligan is a Cu-Au porphyry deposit and the goal of the survey was to delineate Chapter 1. Introduction 7 300 250 200 C/5 co cd Q 150 100 50 o o 900 Hz IP x .1 o- — - o 9 0 0 H z Q x . 1 -B 7200 Hz IP x.1+25 a 7200 Hz Q x .1+25 -0 56 kHz IP x .1+75 0 - — 0 56 kHz Q x . 1+75 4 o - G . _ e _ ^ ) . o -o - © - e - - e - o- -a - o - - e - e - e — e — e — e — o o o — o o 'o—e—e—e—o 1 o—e—e—e—o' o n n n =& 9D00 9200 9400 9600 Northing 9800 10000 Figure 1.5: The inphase and quadrature components of the Mt. Milligan data. 900 Hz (circles), 7200 Hz (squares), and 56000 Hz (diamonds). Inphase (solid Une) and quadrature (dashed line). the intrusive stock and the surrounding zone of mineralization (Oldenburg, et al., 1997). The data were collected using three horizontal coplanar transmitter and receiver coil pairs with measurement frequencies of 900 Hz, 7200 Hz, and 56000 Hz and coil separations of 8 m, 8 m, and 6.4 m respectively. I will concentrate on Line 12625E which is 1 km long with measurements approximately every 10 m. The inphase and quadrature components of the data are shown in Figure 1.5. In contrast to the Sullivan data the Mt. Milligan data appear to be interpretable and free of any serious errors. The inphase component of the Mt. Milligan data asymptotes to zero at low frequency and both components of the data vary smoothly from station to station. When inverting the Mt. Milligan data I assigned a standard deviation of 25 Dighem units (5 PPM) to the 900 Hz measurements and a standard deviation of 50 Dighem units Chapter 1. Introduction 8 Northing (m) Figure 1.6: Inversion results from the Mt. Milligan data. Top panel: the final data misfit values from each inversion plotted by station location. Bottom panel: the recovered 1-D conductivity models from each inversion plotted by station location. (10 PPM) to the 7200 Hz and 56000 Hz measurements. The target data misfit for each sounding was equal to 6. I can see from the inversion results in the top panel of Figure 1.6 that it has been possible to achieve the desired misfit at most stations. However, the recovered conductivity model shown in the bottom panel of Figure 1.6 shows some regions of anomalously low conductivity near the surface. This has been described by Ellis & Shehktman (1994) as an "air layer" effect and it is commonly seen in inversion results when the incorrect measurement height has been used. The measurement heights used in the inversions are the values that were measured by the Dighem system and it is well known that these values can be incorrect in certain situations. In Fraser (1986) it was suggested that both extreme terrain and tree cover can lead to measurement height errors. I know that the region of the Mt. Milligan deposit has topographical variation and it is covered by large trees in some areas. Therefore, it is likely that the presence of measurement height errors in the Mt. Milhgan data is the cause of the distortions in the recovered conductivity models. Chapter 1. Introduction 9 From a mathematical viewpoint the Sullivan and Mt. Milligan data sets are simi-lar in that each is difficult to interpret because a survey parameter (coil separation or measurement height) was not accurately provided. Therefore, these two field data sets invite the solution of the general problem of inverting EM data when geometrical survey parameters are inaccurate. Thus an inversion methodology must be generated to find both the electrical conductivity and the erroneous survey parameter. The major part of the work contained in this thesis is concerned with solving this problem. A second problem which will be investigated concerns the estimation of noise when inverting geophysical data. Correctly estimating the amount of random noise associated with a given data set can increase the amount of information that can be extracted via inversion. Noise estimates for inverse problems are usually made by inspecting the data and assigning the level to which the data should be fit. By examining the erratic nature of the Sullivan data shown in Figure 1.1 it is clear there is a fair bit of noise present however, the task of estimating the amount of noise is an extremely difficult one. While information about the accuracy of Max-Min measurements is available, it can not provide me with details about the amount of random error associated with a given measurement. A more desirable method of estimating noise would make use of the information contained in the data itself. Generalized cross validation (GCV) is such a method. GCV is a statistical method which can be used to estimate the amount of random noise associated with a given data set. It has been applied to widely used to estimate noise in linear problems (Wahba, 1990) and has recently been applied to geophysical inverse problems (Haber, 1997). The application of GCV to the 1-D EM inverse problem would help to ameliorate some of the problems caused by estimating errors by inspection. Chapter 1. Introduction 10 1.2 Outline of the Thesis The goal of this thesis is to develop an inversion methodology that will be able to recover a 1-D conductivity distribution as well as a geometrical survey parameter from frequency domain EM data. In order to address this problem effectively some background information is needed. The necessary information is provided in Chapter 2. This section includes a description of frequency domain EM data and how it is generated using airborne and ground based EM systems. It also includes a discussion of how the 1-D EM inverse problem is currently solved. This review should provide the tools needed to effectively attack the problem. Chapter 3 presents a discussion of the problems associated with parameter errors in EM data. It describes the way in which parameter errors occur and how these errors affect the EM data. A description of how these errors affect inversion results is also provided. Each of these topics is covered in terms of both AEM and HLEM surveys. The final section of this chapter attempts to confirm the presence of geometrical survey parameter errors in the two field data sets. Chapter 4 deals with the development of an inversion methodology which has the ability to recover a function and a parameter. The changes that must be made to the inversion methodology in order to include an extra parameter are addressed first. The most important of these changes is made to the model objective function. The way in which the model objective function is defined greatly affects the stability of the algorithm therefore, it must be defined carefully. Once all of the changes are completed the new algorithm is applied to both synthetic and field data sets. Chapter 5 is concerned with the application of GCV to the 1-D EM inverse problem. A brief description of how GCV is used to estimate random noise is presented. A discussion of how it has been applied in both linear and non-linear inverse problems is also included. Chapter 1. Introduction 11 An inversion algorithm is developed that uses GCV to estimate noise. The algorithm is developed for the 1-D EM inverse problem to recover a conductivity structure, as well as to recover both conductivity and a geometric survey parameter. The algorithm is then tested on synthetic data. A discussion of the applicability of GCV to small data sets is also presented. The final chapter will summarize the results of the entire thesis. Chapter 2 Background Information As was stated in Chapter 1 this thesis will address a number of problems that can arise when attempting to invert problematic field data sets. In order to do this it is necessary to have a good understanding of all aspects of the inversion process. The purpose of this chapter is to provide the required background information. I will begin by denning the EM data used here. This will be followed by a discussion of how the 1-D EM inverse problem is currently solved. This will provide the tools needed to deal with the proposed inverse problems. 2.1 Denning the E M Data This discussion of frequency domain geophysical EM methods will be focused toward small loop EM surveys. These surveys employ a small loop of wire as a transmitter. In order for a transmitter loop to be considered small it must have a radius that is much smaller than the distance that separates the transmitter and the receiver. The receiver is typically a small coil that measures the time rate of change of the magnetic flux density. Methods which employ this type of transmitter receiver pair are commonly called loop-loop EM surveys. While it is possible to measure the magnetic flux density itself using a flux-gate magnetometer as a receiver, both the airborne and ground based surveys I am concerned with use coils as receivers and therefore, I will not discuss the details of using magnetometers, but the alteration of my methodology to work with flux is simple. 12 Chapter 2. Background Information 13 Figure 2.1: The four common types of transmitter-receiver geometries used in loop-loop EM surveys. The figure is labeled as follows : Tx - the transmitter, Rx - the receiver, HC - horizontal coplanar, VC - vertical coplanar, CA - coaxial, and PP - perpendicular. There are many possible configurations of transmitter and receiver coils in loop-loop surveys. The four most common geometries are horizontal coplanar (HC), vertical coplanar (VC), coaxial (CA), and perpendicular (PP). Examples of these transmitter-receiver pairs are shown in Figure 2.1. In the HC geometry the axes of both coils are perpendicular to the surface of the earth and therefore, the coils he in a common plane parallel to the surface of the earth. The VC geometry can be thought of as the HC geometry shown in Figure 2.1 rotated 90° out of the page. The VC coils He in a common plane perpendicular to the surface of the earth and the axes of the coils are parallel to the surface. The coils in the CA geometry he on a common axis which is parallel to the surface of the earth. Finally, in the PP geometry, the axes of the two coils are perpendicular to one another and both of these axes he in a plane perpendicular to the surface of the earth. Each of these coil geometries enables the EM system to sample the earth in a different way and therefore they will each generate different data. Since the field data sets I am investigating contain data collected using only the HC geometry I will not discuss any of the other geometries within the thesis. However, it is important to remember that the principles of loop-loop EM measurements remain the Chapter 2. Background Information 14 S Tx M | Rx "i™ X Z ; h • y — ' \ Figure 2.2: Geometry of the generic horizontal coplanar EM system. Tx and Rx denote the transmitter and receiver coils respectively, s is the coil separation, and h is the measurement height above the surface of the earth. same for any geometry, and all of the methods that are developed within the thesis could easily be applied to other transmitter-receiver geometries (Zhang, et al., submitted for publication). 2.1.1 Generic Horizontal Coplanar Systems A diagram of the generic horizontal coplanar EM system is shown in Figure 2.2. In this configuration the transmitter coil is located a vertical distance h above the surface of the earth and is oriented with its axis perpendicular to the surface of the earth. The receiver and transmitter coil he in the same horizontal plane and are separated by a horizontal distance s. While I want to understand what goes on in a measurement taken over a given conductivity structure it is always good practice to start simple. Therefore, I will begin by considering the case when the transmitter and receiver pair are in free space, shown in Figure 2.3(a). A time varying current I is set up in the transmitter loop. This current has the form /(«) = Joe**, (2.1) where I0 is the amplitude of the current, o> is the angular frequency, and t is time. This current produces the primary magnetic field Hp. It is known from Faraday's Law that Chapter 2. Background Information 15 Figure 2.3: Cartoon showing the behavior of a horizontal coplanar EM system (a) in free space and (b) in the presence of a conductive body. Hp is the primary field, Hs is the secondary field, J is the source current, and Is are the eddy currents. Chapter 2. Background Information 16 a loop in the presence of a time varying magnetic field will have an electromotive force (emf) associated with it. The emf measured in the receiver coil is referred to as the primary voltage V0, and it has the form V, = - w (2.2) where $ is the magnetic flux through the coil which is defined as $ = NRJ B-dS, (2.3) where NR is the number of turns in the receiver coil, S is the surface defined by the receiver coil, and B is the magnetic flux density. By substituting Equation 2.3 into Equation 2.2, and by making use of both the constitutive relation B = UQH and the time dependence of the magnetic field, I am left with Vo = -iufioNR J Hp • dS, (2.4) where / x 0 is the magnetic permeability of free space. Since the radius of the receiver coil is assumed to be much smaller than the separation of the transmitter and receiver coil I can assume that Hp is constant across the surface S. This simplifies the expression for the primary voltage to Vo = -UJ^NRARHp(S), (2.5) where AR is the area enclosed by the receiver coil and HP is the vertical component of the primary field at the center of the receiver coil. The primary voltage is commonly normalized by the amplitude of the current in the transmitter. This gives me the mutual impedance of the two coils in free space Chapter 2. Background Information 17 Zo is often referred to simply as the free space impedance. The next step is to introduce the effect of a conductive body near by the transmitter-receiver pair, shown in Figure 2.3(b). When HP interacts with a conductive body cr the primary field will induce eddy currents 7s in the conductor. These currents will in turn produce a secondary magnetic field HS. The field sensed at the receiver coil HR will be a combination of both HP and HS. Following the same steps I used to calculate VQ I see that the voltage V measured in the receiver coil will have the form V = -iufioNnAR (HP(s) + fff (u;, cr, h, s)) . (2.7) In Equation 2.7 the primary field at the receiver coil only depends upon the coil separation while the secondary field is a function of the frequency, the conductivity structure, the measurement height, and the coil separation. The voltage V is normalized by Io to give an expression for the mutual impedance Z of the two coils in the presence of a conductive body Z = y- (2.8) The quantity which the EM system measures is V. However, a plot of V as a function of frequency is not very informative. In order to glean some information about subsurface conductivity structure from EM measurements the data are presented as the mutual impedance ratio Z/Zo or more commonly the relative change in the impedance ZQ ZQ where AZ = Z — Zo- By substituting Equations 2.5 and 2.7 into the definitions of Zo and Z, I get Chapter 2. Background Information 18 This is the value that is recorded by the EM system. However, it is important to note that HS lags HP, due to the inductive interaction of HP and the conductor, and therefore each datum is a complex quantity. A single datum from Equation 2.10 can be split into a real and imaginary part. The real part is called the inphase component and the imaginary part is called the quadrature component. The final form of the data from the generic horizontal coplanar EM system is Inphase = T7p7~\ x (211) H z \ S ) and Im (H*(u;,(T,h,s)) Quadrature = V ' v ' ' x f. (2.12) H z { S ) where £ is a multiplicative factor. This factor is usually set equal to 100 or 106 such that the data units are either percent or parts per million of the primary field. The units and multiplicative factors have not been defined in the above equations so as to make them as generic as possible. 2.1.2 Specific Horizontal Coplanar Systems As mentioned previously, the voltage V is actually measured by the system. Therefore in order to generate the data in the form shown in Equations 2.11 and 2.12 the primary voltage Vo must be calculated. When using the horizontal coplanar geometry the primary field value at the receiver is easy to calculate. It is equal to the vertical component of the magnetic field from a dipole source. where mx is the dipole moment of the transmitter coil which is denned as mT = IATNT, (2.14) Chapter 2. Background Information 19 Figure 2.4: Cartoon of a typical airborne EM system, h is the measurement height and s is the coil separation. where AT is the area enclosed by the transmitter coil, and NT is the number of turns in the transmitter coil. By substituting Equation 2.13 into Equation 2.5 I get Vo = iwp0NRAR-^-. (2.15) Airs6 While the horizontal coplanar geometry is used widely in both airborne and ground based surveys there are slight differences in the way the data are generated. The difference lies in how the value for the primary voltage is attained. Therefore I will present a brief discussion of the details of airborne and ground based EM surveys. Airborne E M There are many types of airborne EM (AEM) systems in use. I will concentrate on a Dighem type of AEM system in which the coil pairs are housed in a "bird" which is towed behind a helicopter. A cartoon of a typical AEM system is shown in Figure 2.4. The bird contains a number of coils at fixed separations, each of which can make measurements at a particular frequency. Since it is not possible for the measurement height h to be fixed at a known value while the survey is flown, it is determined using a laser altimeter mounted Chapter 2. Background Information 20 S SB Tx Bx Rx gain control VB V Figure 2.5: The bucking coil system used in AEM. The coils are labeled as follows: Tx is the transmitter, Rx is the receiver, and Bx is the bucking coil. The separations are labeled as follows: s the separation between Tx and Rx and SB is the separation between Tx and Bx. V is the voltage measured at the receiver coil and VB is the voltage measured in the bucking coil. The gain control a is adjusted such that OVB is equal to the primary voltage at the receiver Vo. The output of the system V — OLVB is equal to the voltage induced in the receiver coil by the secondary fields. on the helicopter. The measurement height is calculated by subtracting the estimated distance between the bird and the helicopter from the reading on the altimeter. The calculated value of h is included in the data output from the system. The fact that the coils are rigidly mounted within the bird ensures that the coil separation is constant throughout the survey. This enables the use of a bucking coil to remove the effect of the primary field at the receiver. A discussion of the types of bucking coils commonly used in EM systems is presented in Norton, et al. (1999). My discussion will concentrate on the type used in Dighem AEM systems as described in Fitterman The bucking coil, shown in Figure 2.5, is located between the transmitter and the receiver coil, I want the voltage VB in the bucking coil to be entirely due to the primary (1998). receiver. Since the goal is to remove the primary field from the field sensed at the Chapter 2. Background Information 21 field. Therefore, the distance SB between the transmitter and the bucking coil is chosen such that the amplitude of the primary field at the bucking coil is much larger than the amplitude of the secondary field. In order to satisfy this condition the bucking coil is placed in close proximity to the transmitter coil. The bucking coil voltage is then adjusted using a gain control a (shown in Figure 2.5) such that Vo — QVB- This gain calibration is done prior to the survey on the ground. It is assumed that the ground over which the AEM system is calibrated is highly resistive and therefore, the effect of the secondary fields on the system is negligible. Fitterman (1998) discusses problems that can be encountered when the subsurface is in fact conductive. However, for this investigation I will assume that the calibration procedure is successful. Therefore, the value OLVB will provide an accurate estimate of Vo-Figure 2.5 shows that the bucking coil and receiver coil are wired together. This allows the estimated primary voltage to be removed from the measured voltage to leave the secondary voltage. The AEM system also normalizes the secondary voltage by OLVB-This is equivalent to substituting the gain adjusted bucking coil voltage into Equation 2.10 such that Z V - aVn T ~1 = -^r1' (2-16) Zo OLVB and since I am assuming that OLVB = Vo I get Z Hf((jj,a,h,s) / n * -1 = Ww • (217) Due to the fact that the secondary field values are very small in AEM surveys, the data are expressed in parts per million (PPM) of the primary field. Therefore, a single AEM measurement will result in data of the following form Inphase(PPM) = M V ^ ' M ) x 1 0 > ( 2 . 1 8 ) Chapter 2. Background Information 22 s Figure 2.6: Cartoon of a typical HLEM system, h is the measurement height and s is the coil separation. and Im (Hf(u,(T,h,s)) „ Quadrature(PPM) = ^ , } x 106. (2.19) H z \ S ) Ground based E M There are many types of ground based horizontal loop EM (HLEM) systems in use. I will concentrate on a Max-Min type of system in which the transmitter and receiver coils are independent and are carried by two people. This makes it possible to take measurements at various preset coil separations as well as at a range of frequencies. A cartoon of a typical HLEM system is shown in Figure 2.6. While the transmitter and receiver are not rigidly connected they are linked via the reference cable.The reference cable is used to pass information about the phase and amp-litude of the primary field to the receiver, and to provide the operators with an estimate of the coil separation. The HLEM system can be used for both parametric soundings (when the frequency is varied at a fixed coil separation) and geometric soundings (when the coil separation is varied at a fixed frequency). While analysis will focus on paramet-ric sounding results all of the methods developed in the thesis are equally applicable to geometric sounding data. Since the coils are being carried by people the height of the transmitter coil may not be equal to the height of receiver coil. It is also possible for the coils to be tilted Chapter 2. Background Information 23 incorrectly with respect to one another such that they are not coplanar. However, for my investigation, I assume that the effects of such errors are negligible and hence it is assumed that the measurement height is constant and the coils are oriented properly. When taking a measurement, the transmitter and receiver coils must be positioned a preset distance apart from one another. Their separation is commonly estimated using the length of the reference cable or by using station location markers. The preset separation value is used to calculate the voltage Vc that would be induced in the receiver by the primary field alone. This gives me and since Vc is calculated using Equation 2.5 I can assume Vo = Vc such that f - l = % p M . (2.21) Due to the fact that both the transmitter and receiver coils are positioned close to the surface of the earth, the amplitude of HLEM data is much larger than AEM data. Therefore, HLEM data are usually expressed in units of percent of the primary field. A single HLEM measurement will result in data of the following form Inphase(% = V „p / ' ) } x 100, (2.22) H z \ S ) anc / < w x Im(H?(w,a,h,s)) Quadrature(%) = V ' ' x 100. (2.23) H z K S ) Comparison of A E M and H L E M While both AEM and HLEM surveys can be carried out with horizontal coplanar coils there is an important difference between them which should be reinforced. As was men-tioned AEM surveys use coils which are rigidly mounted such that the coil separation s Chapter 2. Background Information 24 is fixed. The measurement height h in these surveys is determined by the AEM system. It is important to note that this value of h may not be correct. HLEM surveys differ in the fact that the transmitter and receiver coils are suspended at a constant height h above the earth. Another difference is that the coils are independent and therefore, they must be positioned at the expected separation s for each measurement. In prac-tice, accurately positioning the coils at a fixed distance s can be difficult. Therefore, both AEM and HLEM each have an important geometrical survey parameter that is not accurately known at the time of data collection. While the presence of these errors is important, most inversion algorithms do not recognize their existence and therefore, AEM and HLEM data are treated identically. 2.2 Details of the Inverse Problem In order to understand the problems that can occur when inverting EM data it is impor-tant to understand the details of the inversion process. I will treat this in two sections: the forward problem and the inverse problem. 2.2.1 The Forward Problem Having discussed the details of how the data from the two EM systems are generated I would like to be able to simulate the response of an EM experiment numerically. The data I collect has the generic form shown in Equations 2.11 and 2.12. The set of measurements from a given sounding which I wish to invert will be stored in a data vector d. The entries of this vector will be the real and imaginary parts of the data values. The total number of data N will be equal to twice the number of measurements and the two parts of the data from the ith measurement di will have the form Re(di) = Re ( f f f K c r , HP(si) (2.24) Chapter 2. Background Information 25 and Im(di) lm (iff (u>,,<7, hi, Si)) Hp(si) (2.25) where each of the subscripted parameters represent the parameter value associated with the ith measurement and £ is some multiplicative factor. It has been shown in Equa-tion 2.13 that the primary field value at the receiver can easily be calculated using the separation of the transmitter and receiver. However, the task of calculating the resultant secondary magnetic fields for a given coil separation, measurement height, frequency, and conductivity distribution is a little more involved. While I know that the subsurface conductivity structure is usually 3-D I choose in-stead to model the response of a 1-D structure. The reason for this is that in many cases the volume of the earth that is being sampled by the EM system can be adequately rep-resented by a 1-D conductivity structure. The validity of this assumption depends upon the size of the "footprint" of the system relative to the scale of the lateral inhomogene-ity I am dealing with. AEM systems have quite a large "footprint"; however, they are commonly used for mapping large scale features. The "footprint" of an HLEM system depends upon the coil separation being used and the separation is chosen according to the size of the feature that is being mapped. Surveys with small coil separations are used to map small scale structures and larger separations are used to map larger struc-tures. In many of these situations the 1-D assumption is acceptable and therefore, I will use a forward modelling code that generates the secondary fields from an arbitrary 1-D conductivity structure. The choice of a 1-D model instead of a 3-D model also provides me with a faster forward modelling routine which is important since I am planning to incorporate it into an inversion scheme in which multiple forward modellings will be ne-cessary. The general solution of the 1-D frequency domain EM problem can be found in Ward & Hohmann (1988) and I have based my solution on their work. I will present an Chapter 2. Background Information 26 Rx surface Gi hi di hi i = 2,... M - l CJM flM OM+I Figure 2.7: Discretization of the problem domain for the forward problem. Tx and Rx denote the transmitter and receiver coils respectively, s is the coil separation, h is the measurement height, hi is the thickness of the ith layer, and <7; is the conductivity of the :th layer. abbreviated description of the specific method used to calculate the secondary fields. The details of the experiment are the same as in Section 2.1. The horizontal coplanar coil geometry is shown in Figure 2.7 above a layered subsurface. The transmitter and receiver coils are separated by a distance s and he in a common horizontal plane at a height h above the surface of the earth. The current I in the transmitter is defined in Equation 2.1. In this coordinate system the z direction is downwards. The subsurface is discretized into M layers. Each of these layers is assigned a thickness, hi, and a conductivity, o~i. These M layers overlie a half-space of conductivity O~M+I- I have assumed that the magnetic permeability of each layer is equal to that of free space. The region between the surface and the coils will be referred to as the 0th layer and the thickness of this layer is equal to the measurement height such that h0 = h. The conductivity of the air is <r0 = 0. As with any EM problem I start with Maxwell's equations. I will make use of the Chapter 2. Background Information 27 frequency domain quasi-static Maxwell's equations in which the effects of displacement currents are neglected. The assumption that the effects of displacement currents are negligible is described in Weaver (1994) to be equivalent to the assumption that the time taken by an EM wave to travel over the region of interest (i.e. the region of the sounding) is much less than the scale time scale over which the EM fields are changing (i.e. a time scale of t = 2TZ/UJ ). For the length scales and measurement frequencies considered in my EM experiments the quasi-static assumption holds for any region within the problem domain. Therefore, the equations that describe the fields in a region of constant conductivity cr are V x E + ifiouH = 0 (2.26) V x H - <TE = J s , (2.27) where E is the electric field, H is the magnetic field, and Js is an electric current source. The symmetry of the loop source allows me to pose the problem in cylindrical coordinates. Since the current source is in the <p direction, and all of the conductivity boundaries are parallel to the f — z plane, all of the induced current is forced to flow in the <p direction. The vertical and radial components of the electric field will thus be equal to zero and therefore I see from Equation 2.26, that the tangential component of the magnetic field is also equal to zero. This leaves me with the following set of equations u;u.0HT = - j j - * (2.28) oz 1 d u>n0Hz = — — {rEj,), (2.29) r or ^ - ^ = „J5, + J.. (2.30) oz oz Each of the E and H field components is a function of r, z, and u> even though it has not been explicitly stated. By substituting Hr from Equation 2.28, and Hz from Chapter 2. Background Information 28 Equation 2.29 into Equation 2.30, I get the following partial differential equation for the tangential electric field 1 d 9 2 (F \ 9 dr (rE*) + k Ej, = iup0Js, (2.31) where A;2 = —IUJLIQCT. Following Ryu et al. (1970) I can transform Equation 2.31 into the Hankel domain to get an ordinary differential equation for the tangential component of the electric field u E,p = iu)p,0js, (2.32) where u2 = A2 — k2, and and J, are the Hankel transforms of Ej, and Js. As mentioned previously, these equations only hold for regions of constant conductiv-ity. Therefore, I must solve Equation 2.32 in each of the M + 2 regions that were denned in Figure 2.7. Then by matching the interface conditions I will be able to propagate the solution from the lower half-space up to the surface. The expression for the secondary electric field at the receiver loop, in the Hankel domain is -2uo'»o -iu)Lio e -RTEJi{\r)\2d\, (2.33) 47T UQ where J\ is a 1st order Bessel function and RTE is the transverse electric (TE) field mode reflection coefficient. The reflection coefficient is defined as Zl -Z0 R t E ~ Zi + Zo' where the intrinsic impedance of the ith interface Z{ is defined as (2.34) Zi Chapter 2. Background Information 29 and the input impedance of the ith layer Zl is denned by the recursion relation I Zi+1 + Zr Zl = Zi Zi + z^ -Uj2/ij (1 + e" -U{2hi -Ui2hi (1 + e--Uj2hi i = M , . . . , l , (2.35) where 17M+I ry 6 — Aftf+1-(2.36) However, the value I want to calculate is Hz(r,u>, z). Using Equation 2.29, and applying the inverse Hankel transform, I get TUT t°° e~2u°ho Hz(r,u;,h0) = / RTEJ0{Xr)X3d\. 4 T T J0 U0 (2.37) The expression for Hz(r,w, ho) in Equation 2.37 is equal to the secondary magnetic field Hf. In order to convert these values into either percent or PPM of the primary field format it is just a matter of normalizing the values to the primary field and multiplying by the appropriate factor. Therefore, I now have a method of calculating the data for use in the inversion algorithm. 2.2.2 The Inverse Problem I have assumed that within the region of the survey the conductivity can be represented as a one-dimensional function of depth crtrue(z). The observed data d°h" is a vector containing 7Y data which are the result of an EM experiment overtop of o~true(z). The data include some unknown amount of measurement noise. The goal of the inversion is to recover the function o~true(z) however, this problem is seriously underdetermined since I have only 7Y data constraints and a function has infinitely many degrees of freedom. Therefore, the solution to the inverse problem has a non-unique solution. This means that if one Chapter 2. Background Information 30 <r(z) can be found which reproduces the data there exist an infinite number of other (T(Z) functions which can also fit the data. In order to deal with the problem non-uniqueness it is necessary to provide information about the specific type of conductivity model I want to recover. This is where it is possible to incorporate any a priori information that is available about the model. The first pieces of information that I can include in the model are that conductivity is always positive and that values for earth materials can vary over several orders of magnitude. Therefore, I will define the continuous model M. to be equal to M = log(a(z)). (2.38) This form of the conductivity is also a good choice because it treats relative conductivity contrasts well. Now that I have chosen a model, more information can be added by introducing a model objective function <pm of the form <Pm{M) = a, J w,(z)(M - Mreffdz + az J wz(z) d(M-Mrefy2 dz, (2.39) dz where M.ref is some reference model. The first term of (pm is a measure of how close M. is to the reference model M.Tef- This term is referred to as the smallest model norm. The second term provides a measure of the derivative of M. in the vertical direction. This term is referred to as the flattest model norm. The parameters a, and az determine the relative importance of the smallest and flattest components of the model norm and the functions w,(z) and wz(z) are spatial weighting functions which can be used to include further information about the model. This formulation provides the ability to recover the smallest model when a, = 1 and az = 0, the flattest model a, = 0 and az — 1, or any combination of the two. Usually I want to find a model that is as featureless as possible and hence a, and az are selected such that the second term in <f)m is dominant over the Chapter 2. Background Information 31 first. When no extra information is available it is common to set ws(z) = wz(z) = 1. This is the general definition of <f)m that I will adopt throughout the thesis. While I would like to be able to recover the continuous function , it is necessary to discretize my model in order to both solve the inverse problem and calculate the forward modelling. Thus the model is divided into M layers of constant conductivity (as shown in Figure 2.7 ) such that m, the discrete representation of Ai, will have the form m=[log{<71),...,log{aM)]T. (2.40) It is important to note that M » N and therefore, the problem remains underdeter-mined. The discrete form of <£m in Equation 2.39 will be cj>m{m) = {m- mTeff [aaWjWa + azWTWz] (m - mref), (2.41) where W3 is the smallest model weighting matrix and Wz is the flattest model weighting matrix. These terms can be combined such that 4>m{m) = {m- mref)TWlWm(m - mTef), (2.42) where Wm is the model weighting matrix which encompasses all of the details incor-porated into my model objective function. The model objective function can also be expressed as the following <M™) = \\Wm(m - mref)\\2, (2.43) where || • || is the Euclidean 2-norm. The design of the model objective function is such that if I simply minimize <pm, and Wm is invertible, the model I recover will be equal to the reference model. Now that I have defined the character of the model I want to recover I can use the observed data to refine the model further. The observed data can be expressed Chapter 2. Background Information 32 mathematically as d°b° = F[mtrue] + e, (2.44) where m t r u e is the discrete representation of the true model, F is the forward operator which generates the data, and e is a vector of length N which contains the noise associated with the data. I will assume that e is Gaussian random noise. The fact that the observed data are contaminated with noise suggests that fitting the data exactly is a bad idea. Therefore, I introduce a data misfit function fa. As its name indicates 4>d is a measure of how well the data predicted from a given model m fits the observed data. Since I have assumed that the errors are Gaussian, a sensible misfit term would be where rji is an estimate of the standard deviation of the noise on the ith datum. I can rewrite 4>d in matrix notation as follows fa=\\Wd(F[m}-d°b°)\\2, (2.46) where Wd is the diagonal data weighting matrix which is equal to Wd = diag(-,...,—). (2.47) \m VNJ Different misfit criteria may be required in some problems and they can be easily defined by choosing a different Wd, a different norm, or by completely redefining fa. Using the definition of fa from Equation 2.46 I see that it is a random variable with a %2 distribution and therefore, it will have an expected value approximately equal to N. Thus the target misfit <j>d should also be approximately equal to N. Chapter 2. Background Information 33 Now that I have dealt with the noise in the data and the problem of non-uniqueness I can express the goal of the inversion as: Find m which minimizes <j>m such that fa = <pd. (2.48) This can also be expressed in the form of an optimization problem where I want to minimize a global objective function <J> which is equal to <f>(m)= | | V 7 d ( F H - ^ ) | | 2 + ^ | | W m ( m - m r e / ) | | 2 , (2.49) where (3 is the regularization or trade-off parameter. The problem will then become: Minimize <f> = fa + {3<j>m such that fa = <f>d, ^ (2.50) The methods used to solve this problem depend upon the details of the forward modelling. When the forward modelling F[m] is linear the data can be expressed in general as F[m] = Gm, (2.51) where G is an N x M matrix. Substituting Equation 2.51 into the global objective function from Equation 2.49 results in the objective function for the linear inverse problem <K™) = \\Wd(Gm-d°b')\\2 + (3\\Wm(m-mref)\\2. (2.52) The minimization of the linear functional in Equation 2.52 for a particular /3 is solved by calculating the gradient of <f> and setting it equal to zero. The gradient g is obtained by differentiating <j) with respect to the model m g(m) = 2GTWjWd{Gm - (f**) + 2f3WZWm{m - mref). (2.53) Setting g(m) = 0 leads to the matrix system of equations [GTWjWdG + L3WlWm)m = GTWjWdd°b> + f3WZWmmref, (2.54) Chapter 2. Background Information 34 <f>d 4>d /? -> o Figure 2.8: Cartoon of a typical Tikhonov curve, fa is the data misfit, <j)d is the target misfit, <f>m is the model objective function, and f3 is the trade off parameter. which can be solved to provide m which minimizes Equation 2.52. By carrying out the minimization at a number of (3 values it is possible to generate a plot of fa({3) versus <f)m((3). This plot is called the Tikhonov curve. A cartoon of a typical Tikhonov curve is shown in Figure 2.8. This curve illustrates how the choice of /3 dictates the values of $d and (j)m and therefore, the character of the recovered model. When /3 gets large, 4>(m) ~ /3(f>m and the minimization is strictly searching for the minimum structure model. This can result in models that do not fit the data very well. Whereas, when (3 approaches zero, <j)(m) m fa. In this case the minimization will attempt to find the model that provides the best fit to the data however, this can result in models with large amounts of structure. In between these two extremes the relative importance of <f>m and fa is determined by the value of j3. The Tikhonov curve quantifies the way in which the choice of /3 regularizes the solution. In linear inverse problems the process of finding j3 which provides the desired fa value is straight forward. By plotting the Tikhonov curve it is possible to determine whether a (3 exists such that fa = 4>d. If such a /3 does exist, as shown in Figure 2.8, it is a simple Chapter 2. Background Information 35 matter to find it. In the case that the forward modelling F[m] is non-linear however, the task of finding m which minimizes the objective function from Equation 2.49 for a given (3 is more involved. It is necessary to linearize the problem and solve it iteratively. This is done using a standard Gauss-Newton methodology. As in the linear case the minimization of the functional in Equation 2.49 for a partic-ular /3 is solved by calculating the gradient of <p and setting it equal to zero. The gradient g of (f> is equal to g(m) = 2JTWjWd(F[m] - d°b') + 2/3W^Wm(m - m r e / ) , (2.55) where J is the Jacobian matrix defined as dm The gradient from Equation 2.55 is a non-linear equation in m and as a result, it is necessary to linearize g about some current model mk. The linearized version of g is expressed as g(mk + 8m) = g{mk) + H(mk)6m, (2.57) where H(mk) is the Hessian matrix, evaluated at mk, which is defined as H(mk) = = 29Jp^WJWd{F[mk] _,*) + ... omk omk 2J{mk)TWjWdJ(mk) + 2/3W£Wm. (2.58) In order to avoid evaluating the derivative of the Jacobian, the Hessian is approximated as H(mk) * 2J(mk)TWjWdJ{mk) + 2/3W^W m . (2.59) Chapter 2. Background Information 36 Setting g(mk + 8m) — 0 yields the matrix system of equations [J(mk)TWjWdJ(mk) + L3WlWm)8m = ... J(mk)TWjWd(d°b> - F[mk}) - fiW^W^m - mref), (2.60) which can be solved to give 8m such that the model for the next iteration wi l l be equal to mk+i = mk + 8m. The process of linearization and iteration is continued unti l the solution converges to the minimum of Equation 2.49. This process can be repeated for a number of j3 values, as in the linear case, to find the f3 which satisfies (f>d — cj>d. Chapter 3 Geometric Survey Parameter Errors In order to interpret EM data it is necessary to know the parameter values that define the survey from which the data were collected. Without knowledge of these parameters the data are merely a set of numbers whose relationship to one another is incomprehens-ible. For this reason the pertinent survey parameter values are recorded with the data output from the EM system. However, both AEM and HLEM surveys are susceptible to geometric survey parameter errors that can cause problems in the interpretation and/or inversion of the data. In Section 1.1 I conjectured that some of the problems associated with the interpreta-tion and inversion of the field data sets may be the result of geometric survey parameter errors. The goal of this chapter is to show that these problems are indeed due to geomet-ric survey parameter errors. I will discuss the details of how geometric survey parameter errors occur and how they affect both EM data and the inversion results. This will provide me with the ability to identify the presence of geometric survey parameter errors in the field data sets and inversion results. However, before I begin it is necessary to introduce the notation I will use throughout the thesis as well as the synthetic models with which I will be concerned. Notation In order to introduce notation I will consider a generic survey parameter p. The value of p that the EM system records in the data output is referred to as the expected or 37 Chapter 3. Geometric Survey Parameter Errors 38 estimated value of p and it is represented with the variable peat. The true value of p is represented with the variable ptrUe- The parameter error is represented with the variable 8p and is defined as 8p = Ptrue ~ Pe,t- (3-1) This notation will be applied to all of the parameters discussed in the thesis. Synthetic Examples Throughout this chapter I will use synthetic data from three EM surveys. The first survey simulates a typical AEM survey, the second simulates an HLEM survey performed at a small coil separation, and the third simulates an HLEM survey performed at a larger coil separation. A description of the details of each survey and the conductivity models are provided along with an example of synthetic data from each survey. The data shown are free of parameter errors and no noise has been added so that it is possible to become familiar with the uncontaminated form of the data before I begin the investigation of parameter errors. The AEM survey consists of measurements of the inphase and quadrature compo-nents of the secondary fields at 10 frequencies ranging from 110 Hz to 56320 Hz. Each measurement is taken at a precise coil separation of 10 m and at some height above the earth. I will refer to this type of measurement as the AEM sounding. The conductivity model over which the AEM soundings will be simulated is shown in Figure 3.1(a). It consists of a 20 m thick 0.1 S/m layer buried 30 m deep in a 0.01 S/m halfspace. The data shown in Figure 3.1(b) are acquired at a height of 30 m. It is important to note that to simulate an AEM sounding it is necessary to define the true measurement height value. Both of the HLEM surveys consist of measurements of the inphase and quadrature Chapter 3. Geometric Survey Parameter Errors 39 AEM Small HLEM Large HLEM (c) (e) -0.5i (a) w E > 5-1.5 T J c o O o -2.5 -0.5i -1 T J c o O -1.5 -2.5 -0.5i E —» cn -1.5 o O -2.5 50 100 Depth (m) 150 50 100 Depth (m) 150 50 100 Depth (m) 150 2500i (b) S2000| Q. p-§ 1 5 0 0 | o C01000 co g- 500 O O Inphase x x Quad. O x O x 9 a-X x O O 10 10J 10' Frequency (Hz) 10 (d) g, O O Inphase x x Quad. O oi^  6 6 o 3 ° 10 30 20 10 10" Frequency (Hz) (f) g 10 I o 5 6 o 3 o •o-10| CO 8-20I CO cL-30 -40 o O o O O Inphase x x Quad. -50' 10' 10' 10' Frequency (Hz) 10" Figure 3.1: Conductivity models and data from the synthetic examples. Panels (a),(c), and (e) show the conductivity models used for the AEM, small separation HLEM, and large separation HLEM examples. Panels (b), (d), and (f) show synthetic data (inphase - circles, quadrature - crosses) used for the AEM, small separation HLEM, and large separation HLEM examples. Chapter 3. Geometric Survey Parameter Errors 40 components of the secondary fields at 10 frequencies ranging from 110 Hz to 56320 Hz, each with a measurement height of 1 m. The difference between the two synthetic HLEM examples is the coil separation; the first survey is performed with a coil separation of 10 m and the second has a coil separation of 50 m. The conductivity model used for the small separation HLEM soundings is shown in Figure 3.1(c). It consists of a 20 m thick 0.1 S/m layer buried 10 m deep in a 0.01 S/m halfspace. Figure 3.1(d) shows the data from a measurement with a coil separation of 10 m. Since the depth of penetration increases with coil separation it is possible to look at deeper structures therefore, I chose a different model for the large coil separation example. The conductivity model I use for the large separation HLEM soundings is shown in Figure 3.1(e). It consists of a 30 m thick 0.1 S/m layer buried 20 m deep in a 0.01 S/m halfspace. The synthetic data generated with a coil separation of 50 m is shown in Figure 3.1(f). Now that I have introduced the notation and synthetic examples it is possible to begin the discussion of parameter errors. 3.1 Measurement Height Errors in A E M As mentioned in Section 2.1.2 the transmitter-receiver pairs in AEM surveys are rigidly mounted within the bird and therefore, the coil separation is assumed to be fixed. Because the bird is towed below the helicopter variations in height are unavoidable and therefore, the geometric survey parameter of interest in AEM surveys is the measurement height. 3.1.1 How do measurement height errors occur? The true measurement height htrue is denned as the distance from the bird to the surface of the earth. However, the measurement height which is recorded in the data file is a value calculated by the AEM system. Since this value may not be equal to htTue I will refer Chapter 3. Geometric Survey Parameter Errors 41 Figure 3.2: Cartoon of an AEM system showing the values used to calculate measurement height h. a is the altitude of the helicopter and b is the distance from the helicopter to the bird. to it as the estimated measurement height hest. Using Equation 3.1 the measurement height error 8h can be expressed as Sh = h true (3.2) A cartoon of the AEM system and the values used to calculate heat are shown in Figure 3.2. The altitude of the helicopter a is determined using a laser altimeter and the distance from the helicopter to the bird b is assumed to be constant. The assumption about the fixed value of b is based on the fact that the bird, and the cable that connects it to the helicopter, are designed so that b is equal to a known value when the helicopter is traveling at a constant speed. Using the values of a and b, h is calculated to be h = a — b. (3-3) By substituting the subscript est and subscript true versions of Equation 3.3 into Equa-tion 3.2 I get the expression Sh = 6a — 6b, (3.4) Chapter 3. Geometric Survey Parameter Errors 42 Case 1. Case 2. Figure 3.3: The two situations which that result in measurement height errors. Case 1. Sh > 0, measurement height underestimated because of tree cover. Case 2. Sh < 0, measurement height overestimated due to extreme terrain. which defines the measurement height error in terms of errors in a and b. Errors possibly exist in a and b however, for my purposes the bottom line is the value of Sh. When Sh > 0 the bird is farther away from the surface than expected and when Sh < 0 the bird is closer to the surface than expected. Figure 3.3 illustrates two possible cases in which measurement height errors can exist. 3.1.2 How do measurement height errors affect the data? In order to assess the effect of measurement height errors I must first investigate how the data are affected by measurement height. Letting htrue denote the true height, the AEM data from Equations 2.18 and 2.19 are H^(UJ,(T, htrue, S)" Inphase(PPM) = Re x 10e (3.5) and Quadrature(PPM) = Im [ ^i"^*™,*) } ^ 1 Q 0 (3.6) Hp{s) By generating synthetic data from AEM soundings with different htrue values it will be possible to investigate the interpretation problems caused by measurement height errors. Figure 3.4 shows the inphase and quadrature components of AEM sounding data Chapter 3. Geometric Survey Parameter Errors 43 4000, (a) < 3500 00 Frequency (Hz) 10 10 Frequency (Hz) Figure 3.4: Synthetic data from AEM soundings with h t r u e = 24 m, 30 m, and 36 m. (a) The inphase component of the data ( h t r u e — 24 m - diamonds, h t r u e = 30 m - circles, and htme — 36 m - squares) (b) The quadrature component of the data ( h t r u e = 24 m -diamonds, h t r u e = 30 m - crosses, and h t r U e — 36 m - squares). with h^ue values of 24 m, 30 m, and 36 m. The inphase and quadrature component of the h^ue = 30 m data are represented by circles and crosses respectively in order to show that they are identical to the data plotted in Figure 3.1(b). Both the inphase and quadrature component of the h t r U e = 24 m data are represented by diamonds and the htrue = 36 m data are represented by squares. The curves in Figure 3.4 clearly illustrate the effect of measurement height on AEM data. In general, the closer the bird is to the surface, the greater the signal, and vice versa. However, this effect does exhibit some frequency dependence. At lower frequencies the difference in the data is not large while at the highest frequency a ±6 m change in measurement height can lead to a change of greater than ±30% of the data amplitude. It follows that if h t r U e = 30 m, and if the AEM data were interpreted with an incorrect measurement height (either 24 m or 36 m), the high frequency data would need to be assigned an error of at least 30% of the data amplitude. Otherwise, any interpretation would be the result of fitting "noise" caused by the presence of measurement height errors. Chapter 3. Geometric Survey Parameter Errors 44 he»t 4>m fa 24 m 0.136 20.0 20.0 30 m 0.127 20.0 20.0 36 m 0.353 25.1 20.0 Table 3.1: Results from the inversion of AEM data with heat values of 24 m, 30 m, and 36 m when htrUe is equal to 30 m. <f)m is the recovered model norm, fa is the data misfit achieved by the inversion and <f>d is the target misfit. These values correspond to results plotted in Figure 3.5 3.1.3 How do measurement height errors affect the inversion results? In order to obtain some insight about the effects of assuming an incorrect measurement height I will invert a single data set, acquired at a true height of 30 m, three different times. The results of an actual sounding were simulated by adding Gaussian random noise with a standard deviation equal to 5% of the data amplitude + 10 PPM was added to each synthetic datum. The noisy data were then inverted three times using heat values equal to 24 m, 30 m and 36 m respectively. The reference model for each inversion was a 0.01 S/m halfspace and the target misfit was equal to 20. The conductivity models recovered from the three inversions are plotted in Fig-ures 3.5(a), 3.5(c), and 3.5(e) and the observed and predicted data are plotted in Fig-ures 3.5(b), 3.5(d), and 3.5(f). The final data misfit and model norm values from each inversion are shown in Table 3.1. From the data misfit values, and the data plotted in Figure 3.5, I can see that the inversion algorithm has been able to recover models that fit the data to the desired %2 misfit level of 20 for both heat — 24 and 30 m and has done fairly well (fa « 25) for heat = 36 m. However, I can also see from both the model norm values and the recovered models plotted in Figure 3.5 that the conductivity structure recovered for each heat is different. Since the inversion results from heat = 30 m (8h = 0) are free of measurement height errors I will adopt them as the results against which I compare the other inversions. Chapter 3. Geometric Survey Parameter Errors 45 24 m h = 30 m h = 36 m (a) (e) 55 jr-J - True Model — Recovered Model (c) I —-0.5 > 1-1.5| •o c o O -21 -2-2.5] -31 3000. 50 Depth (m) 150 — True Model — Recovered Model 50 100 Depth (m) 150 -•- True Model — Recovered Model 50 100 Depth (m) 150 (b) (d) I (f) 10' 10' Frequency (Hz) 3000, 10* 10" Frequency (Hz) Figure 3.5: Results from the inversion of AEM data with h values of 24 m, 30 m, and 36 m when htTUe is equal to 30 m. Panels (a), (c), and (e) show the true (dashed line) and recovered (solid line) conductivity models and panels (b), (d), and (f) show the observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature - dashed line) from the three inversions. Chapter 3. Geometric Survey Parameter Errors 46 The model recovered with an estimated measurement height of 24 m (Sh > 0) plotted in Figure 3.5(b) has a region of low conductivity at the surface which is not present in Figure 3.5(d). This has been described as the "air layer" effect by Ellis & Shekhtman (1994) and is a result of the estimated height being smaller than the true value. I also notice that, as a result of the "air layer", the conductivity structure recovered using hest = 24 m is slightly pushed down in comparison to the structure recovered with heBt = 30 m. This distortion is an attempt to decrease the amplitude of the predicted data. The "air layer" is trying to provide the extra distance between the surface and the bird that is needed recreate the survey that generated the observed data. The distance by which the recovered conductivity structure is pushed downwards is approximately equal to Sh. I can see from the recovered model plotted in Figure 3.5(f) that the opposite effect is observed when the data is inverted with an estimated measurement height of 36 m (Sh < 0). In this case there is a region of high conductivity at the surface which is not present in Figure 3.5(d). It also appears that in the model recovered using h = 36 m the conductivity structure below the high conductivity region is distorted and pulled up towards the surface in comparison to the structure recovered using h = 30 m. In this case the algorithm has a harder time correcting the problems caused by Sh. Since it is not possible to decrease the distance between the bird and the surface it is necessary to find another way to increase the amplitude of the data. This is achieved by increasing the conductivity of the model near the surface. The extra structure at the surface causes the rest of the recovered conductivity structure to be pulled up towards the surface. The presence of these features are due to the estimated measurement height being larger than the true value. It has been shown that the presence of measurement height errors in AEM data can have a definite effect upon inversion results. From the above examples I see that while it is still possible to achieve the target misfit the recovered models are distorted by Chapter 3. Geometric Survey Parameter Errors 47 the incorrect measurement height estimates. While this makes it difficult to assess the presence of measurement height errors when the true conductivity model is not known at least I have an idea of the ways in which the inversion results can be affected. 3.2 Coil Separation Errors in H L E M data As mentioned in Section 2.1.2 the transmitter and receiver coils in HLEM surveys are independent of each other. The coils are carried separately by two individuals and varia-tions in measurement height are considered to be negligible and as such the measurement height is assumed to be fixed. However, the separation of the coils is variable. Therefore, the geometric survey parameter of interest in HLEM surveys is the coil separation. 3.2.1 How do coil separation errors occur? Prior to collecting data the HLEM system must be set to operate at a particular coil separation. This preset value is the expected coil separation seat and this is the coil separation value recorded in the data file. The distance which actually separates the transmitter and receiver is the true coil separation strue- Using Equation 3.1, I can express the coil separation error Ss as Ss — strue - seat. (3.7) When a survey is performed on level ground separation errors are likely due to in-correct chaining of the station markers or reference cable. When a survey is performed in extreme terrain separation errors can occur even when the reference cable is chained correctly. These types of errors can occur when passing over topographic features such as gullies, cliffs, and valleys. Figure 3.6 shows the two situations that will result in coil separation errors. Chapter 3. Geometric Survey Parameter Errors 48 Case 1. \Ss > 0 "si.V Case 2. \6s < 0 Seat Figure 3.6: The two situations that will result in coil separation errors. Case 1. 8s > 0, the coils are too far apart, and Case 2. 8s < 0 the coils are too close together. 3.2.2 How do coil separation errors affect the data? In order to assess the effect of coil separation errors I must determine the roles that both stTUe and seat play in the HLEM data. The magnetic fields sensed at the receiver coil depend upon stTue- Therefore, the voltage V that is induced in the receiver is equal to V = -iup0m, ( i f f (Strue) + H^(strue)) (3.8) However, the estimate of the primary voltage calculated by the HLEM system depends upon seat. Therefore, the voltage Vc is denned as Vc = -iu>p0maHP(se!>t). (3.9) By substituting V, from Equation 3.8, and Vc, from Equation 3.9, into the expression for the HLEM measurement, from Equation 2.20, I am left with Z_ _ HSZ JU, *, h, Strue) + HP{strue) - H?(seat) Z0 ~ H?{seat) (3.10) From this expression it is clear that when seat ^ strue the HLEM measurements can not be represented by the data from Equations 2.22 and 2.23. Therefore, the HLEM data must be rewritten in terms of strUe and seat as Inphase(%) = Re fH!(^h^ s^e)\ x l f j 0 + SHP(strue,seat), (3.11) Chapter 3. Geometric Survey Parameter Errors 49 and Quadrature(%) = lm x 100 (3.12) where 8HP is the primary field error and is denned as (3.13) These equations show that both strUe and sest have a definite effect upon HLEM data. This effect can be subdivided into the effects on the normalized secondary fields and the effects on the primary field error. I can see from Equations 3.11 and 3.12 that the real and imaginary parts of the normalized secondary fields behave in a fashion similar to the inphase and quadrature components of AEM data in Equations 3.5 and 3.6. However, this problem is compounded by the presence of 8HP in the inphase component. In order to get more insight into the behavior of 8HP I substitute the definition of the primary field from Equation 2.13 into Equation 3.13 and end up with \\StrueJ J This shows that the presence of the coil separation errors has a direct effect upon the inphase component of HLEM data by way of the 8HP term. From Equation 3.14 I see data will be shifted in the negative direction. Similarly when 8s < 0 (a t r u e < se,t) the All of the effects on the data can be illustrated using the two HLEM examples. I will begin by considering a survey in which three synthetic small separation HLEM soundings are taken over the conductivity model shown in Figure 3.1(c). The seat value for all of the small separation soundings is set to 10 m and the strue value for each sounding is 11 m, 10 m, and 9 m respectively. These values correspond to coil separation errors of (3.14) that when 8s > 0 (strue > sest) the value of 8HP will be negative and therefore the inphase value of 8HP will be positive and the inphase data will be shifted in the positive direction. Chapter 3. Geometric Survey Parameter Errors 50 (a) 50 40 30 £ 20 cu <2 10 >--o--o--^ -o-o-o-"<>"'<3> 0-0 s CO -10 -20 ,-4 - 9m G - O s f ^ l O m B-B s. = 11m true 0--Q ob>--G-e-0-e--o-©-e" L . . Q . . Q . . ^ . . . & - e - a - B - ' - B ' - a -30' 10' (b) 10* 10* Frequency (Hz) 10" 10 10 Frequency (Hz) (c) (d) Q. X ra i>~o~o~^ ~o--o~o~^ ~o--o 0 0 s, = 9m O-O s r U 8 = 10m „ „ T R U E ., ., Q-Q s. = 11 m true Q s oE>-o--e-e-e--e-e-e-e--o &--e--Q—a—B--Q--B—B~-B--0 10 10 Frequency (Hz) -30' 10' 10 10 Frequency (Hz) 10 Figure 3.7: Synthetic data from small separation HLEM soundings with sest = 10 m and •Strue = 9 m, 10 m, and 11 m. (a) Inphase component of the data (strue = 9 m - diamonds, Stme = 10 m - circles, strue = 11 m - squares), (b) quadrature component of the data (•Strue = 9 m - diamonds, strue = 10 m - crosses, strue = 11 m - squares), (c) real part of normalized secondary field (as in panel (a)), and (d) primary field error (as in panel (a)). + 1 m, 0 m, and -1 m. The synthetic data from each of the soundings, with no noise added, are plotted in Figure 3.7. The inphase and quadrature components of the strUe = 10 m data are represented by circles and crosses respectively in order to show that they are identical to the data plotted in Figure 3.1(d). This data will be referred to as the expected data. Both the inphase and quadrature component of the strUe = 11 m data are represented by squares and the 3 t r u e = 9 m data are represented by diamonds. The plots in Figures 3.7(a) and (b) show the inphase and quadrature components of Chapter 3. Geometric Survey Parameter Errors 51 the data respectively. The shift due to 8HP seen in the inphase component, is much larger than the variation of the response with frequency. The two portions of the inphase response, as per Equation 3.11, are plotted separately in Figures 3.7(c) and (d). Figures 3.7(c) and (b) illustrate that the closer the coils are to one another the larger the amplitude normalized secondary fields. These errors are similar to those seen in the AEM case. However, the effect of the primary field errors (Figure 3.7(d)) has a direct effect upon the data. These shifts are independent of frequency and therefore, affect the entire inphase component of each sounding. This type of error can render the inphase component of HLEM data (Figure 3.7(a)) uninterpretable. From Equation 3.14 it can be seen that the magnitude of the shift depends upon the relative size of the coil separation error. Therefore, the shift caused by 8s = 1 m in a survey with seat = 10 m is equivalent to the shift caused by 8s = 5 m in a survey with seat — 50 m. However, it is unlikely to encounter coil separation errors larger than a few metres, even when considering surveys with large coil separation values. Therefore, surveys performed with small coil separations are more likely to be detrimentally affected by coil separation errors. This difference can be illustrated with the synthetic data from the large separation HLEM example. Consider a survey in which three large separation HLEM soundings are taken over the conductivity model shown in Figure 3.1(e). For each of the soundings seat is equal to 50 m and the value of strue at each sounding is 51 m, 50 m, and 49 m respectively. These separation values correspond to relative errors of +2%, 0%, and -2%. The synthetic data from each of the soundings, with no noise added, are plotted in Figure 3.8. The inphase and quadrature components of the strue = 50 m data are represented by circles and crosses respectively in order to show that they are identical to the data plotted in Figure 3.1(f). This data will be referred to as the expected data. Both the inphase and quadrature components of the strUe = 51 m data are represented by squares and the strUe Chapter 3. Geometric Survey Parameter Errors 52 (a) (c) 35 30 25 „ 20 5? ~Z 1 5 o) j= 10| CL £ 51 ^ / a-. „ © \ / / E3 \ V o ^ - 0 - ^ y 0-0 s t r u e = 49m'^ , „ . J 3 G>-Os =50m - 5 b - Q , r u e -10' 10 30 Q-Q s . = 51m t r u e 10 10 Frequency (Hz) 10" —.25 120 i l O | CO cr 0-0 s ( r u e = 49m G-O s / u e = 50m _ „ t r u e _ . , Q--Q s. = 51m t r u e / v "0 ::<3\ 6 Jf vO • 10" 10" 10' Frequency (Hz) 10" (b) (d) 10 10 Frequency (Hz) 6{>-o-0"^ -^ >~o-o--^ --o--o 0-0 s = 49m Q-O s , r u e = 50m „ „ t r u e Q-Q s. = 51m t r u e CL X. ra CD Q B o&-o-o-e-e-0-e-e-e-G -6' 10' 10° 10" Frequency (Hz) 10" Figure 3.8: Synthetic data from large separation HLEM soundings with s e a t = 50 m and Stme = 49 m, 50 m, and 51 m. (a) Inphase component of the data (strue = 49 m -diamonds, strue = 50 m - circles, s t P u e = 51 m - squares), (b) quadrature component of the data (strUe = 49 m - diamonds, strUe — 50 m - crosses, strue = 51 m - squares), (c) real part of the normalized secondary field (as in panel (a)), and (d) primary field error (as in panel (a)). Chapter 3. Geometric Survey Parameter Errors 53 = 49 m data are represented by diamonds. The inphase and quadrature components of the data are plotted in Figures 3.8(a) and (b). The plots in Figures 3.8(c) and (d) are the real part of the normalized secondary field and the primary field errors respectively. By comparing the primary field error plotted in Figure 3.8(d) to the error plotted in Figure 3.7(d) I see that the shift due to a 2% coil separation error is a factor of five smaller than the shift due to a 10% coil separation error. Although the shift errors in the large separation HLEM example are comparatively small they still dominate the inphase component of the data at low frequency. Using these two HLEM examples I have shown that the presence of coil separation errors will have a direct effect upon HLEM data. These errors manifests themselves in a shift of the inphase data. While it is more likely that the inphase component of the data from surveys performed at small coil separations will be rendered uninterpretable by coil separation errors, these problems can also occur in large separation HLEM. I have shown that coil separation errors can also cause discrepancies between the measured and expected secondary field values. This result is similar to the AEM case in that it creates problems with the interpretation, and, while it is present in both the inphase and quadrature components of the data, the inphase portion is usually overshadowed by primary field errors. This means that even if someone were to attempt to deal with the shift in the inphase by only interpreting the quadrature part of an HLEM data set the results may still be incorrect. 3.2.3 How do coil separation errors affect the inversion results? The previous results show that substantial differences exist between the expected data and those obtained when the coil separation is incorrect. Since many inversion algo-rithms can not account for these differences the inversion of data contaminated with coil separation errors may be difficult. In order to explore these difficulties I will attempt to Chapter 3. Geometric Survey Parameter Errors 54 Strue 4>m <f>d 11 m 3.46 xlO" 3 2.20 xlO 4 20.0 10 m 0.116 20.0 20.0 9 m 0.377 3.58 xlO 4 20.0 Table 3.2: Results from the inversion of small separation HLEM data with an s value of 10 m when strue is equal to 11 m, 10 m, and 9 m. <f>m is the recovered model norm, 4>d is the data misfit achieved by the inversion and <f>d is the target misfit. These values correspond to results plotted in Figure 3.9 invert the data from the two HLEM examples. I will start by considering the data from the three small separation HLEM soundings. In order to simulate the data from an actual sounding Gaussian noise with a constant standard deviation of 0.5% of the estimated primary field at the receiver was added to each synthetic datum. Noise with a standard deviation which is a percentage of the data amplitude was not added in this case due to the fact that the soundings contaminated with separation errors would have been assigned unrealistically large amounts of noise. The noisy data are the result of soundings with stTUe values of 11 m, 10 m, and 9 m, and sest = 10 m. Therefore, they were each inverted using a coil separation of s — 10 m. The reference model for each inversion was 0.01 S/m. Each datum was assigned an error estimate of 0.5% of the primary field and the target misfit was set equal to 20. The final data misfit and model norm values from the three inversion are shown in Table 3.2 and the recovered conductivity models are plotted in Figures 3.9(a), 3.9(c), and 3.9(e) and the observed and predicted data are plotted in Figures 3.9(b), 3.9(d), and 3.9(f). From the data misfit values in Table 3.2 it can be seen that it was not possible to achieve the target misfit when strUe is equal to either 9 m or 11 m. In both cases the recovered data misfit was much larger than the target of 20. Figures 3.9(b) and 3.9(f) indicate that the algorithm was not able to fit the inphase data. This suggests that the Chapter 3. Geometric Survey Parameter Errors 55 11 m. Sir = 10 m. Str = 9 m. (a) (c) -0.5 -0.5 -2.5 1r 0.5 (e) f Pi, o :&-0.5T 1 "1 o O) o -2 -4 —- True Model — Recovered Model 50 100 Depth (m) 150 True Model Recovered Model 50 100 Depth (m) 150 —- True Model — Recovered Model 50 100 Depth (m) 150 (b) g (d) o 1-10] ra 0] w o-•15 -20 o o IP Obs. — IPPred. x x QObs. — QPred. -25^ 10 o n 10' 10' Frequency (Hz) 10 7 6 g 5 § 4 o Y 31 (0 S 2l (0 £ . 1 o o IP Obs. — IPPred. x x QObs. QPred. Io' 10 10 Frequency (Hz) 10 50i (f) g ^ 4 0b o o o o ° o J30 w O 20 TJ « 10| <1) O) (0 Q. "-10I o o ip Obs. — IPPred. x x QObs. QPred. -20h 10 10" 10' Frequency (Hz) 10" Figure 3.9: Results from the inversion of small separation HLEM data with an s value of 10 m when strue is equal to 11 m, 10 m, and 9 m. Panels (a), (c), and (e) show the true (dashed line) and recovered (solid line) conductivity models and panels (b), (d), and (f) show the observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature - dashed line) from the three inversions. Chapter 3. Geometric Survey Parameter Errors 56 ^ true K (j>d 11 m 0.086 20.0 20.0 10 m 0.116 20.0 20.0 9 m 0.061 20.0 20.0 Table 3.3: Results from the inversion of small separation HLEM data with an s values of 10 m when s t r u e is equal to 11 m, 10 m, and 9 m and noise estimated to account for inphase modelling errors. (pm is the recovered model norm, (pj is the data misfit achieved by the inversion and (p*d is the target misfit. These values correspond to results plotted in Figure 3.10 differences between the expected data and the erroneous data, which are in fact modelling errors, must be dealt with. One way to handle them in the inversion is to assign to each datum a standard deviation that is about the same size as the modelling error. For the small HLEM sounding this means assigning errors of 25%, 0.5%, and 38% of the primary field to the inphase components of the s t r u e = 11 m, 10 m, and 9 m data respectively. An error of 0.5% of the primary field was assigned to all of the quadrature data. Each of the inversions was carried out again with a reference model of 0.01 S/m and the target misfit set equal to 20. The final data misfit and model norm values from the three inversion are shown in Table 3.3 and the recovered conductivity models are plotted in Figures 3.10(a), 3.10(c), and 3.10(e) and the observed and predicted data are plotted in Figures 3.10(b), 3.10(d), and 3.10(f). I see from the data misfit values in Table 3.3 that the algorithm has been able to achieve the target misfit of 20 for each sounding. However,the data plotted in Fig-ures 3.10(b) and 3.10(f) show that it is not possible to adequately reproduce the inphase observations for the strue = 11 m and 9 m soundings. Only when Ss = 0 is it possible to predict both the inphase and quadrature data(shown in Figure 3.10(d)). The recovered models are quite similar. This is good since it shows that the inversion result is not extremely sensitive to bias errors in the inphase data. Nevertheless, the fact that I have Chapter 3. Geometric Survey Parameter Errors 57 •Strue = 11 HI. Strue = 10 m. •Sfr 9 m. (c) (e) -0.5 -0.5 -2.5 -0.5 -2.5 True Model Recovered Model 50 100 Depth (m) 150 True Model Recovered Model 50 100 Depth (m) 150 True Model Recovered Model 50 100 Depth (m) 150 (b) g o I-10I cs CD c3-15| c Q. C --20] o o IP Obs. — IPPred. x x QObs. -—• QPred. - 2 5 ^ 10 o o 10 10 Frequency (Hz) (d) g o 1 3| co CD ol w 'I IS f 1 -1 10 o o |p Obs. — IPPred. x x QObs. - - Q Pred. 10 10 Frequency (Hz) 10 45 40 S35 0 2 5 | TJ c320| CD S3is| ) ° o O c o o IP Obs. — IPPred. x x QObs. - - Q Pred. o o 10" 10' Frequency (Hz) 10 Figure 3.10: Results from the inversion of small separation HLEM data with an s value of 10 m when strUe is equal to 11 m, 10 m, and 9 m and the noise was estimated to account for inphase modelling errors. Panels (a), (c), and (e) show the true (dashed line) and recovered (solid line) conductivity models and panels (b), (d), and (f) show the observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature - dashed line) from the three inversions. Chapter 3. Geometric Survey Parameter Errors 58 Strue <j>m fa ti 51 m 0.203 20.0 20.0 50 m 0.171 20.0 20.0 49 m 0.152 39.7 20.0 Table 3.4: Results from the inversion of large separation HLEM data with an s values of 50 m when strUe is equal to 51 m, 50 m, and 49 m and noise estimated to account for inphase modelling errors. <f>m is the recovered model norm, fa is the data misfit achieved by the inversion and (f>d is the target misfit. These values correspond to results plotted in Figure 3.11 not predicted the observations more closely is unsettling. It leaves me to wonder if there are artifacts in the recovered models and if there is possibly more information about the earth conductivity that could be extracted if the modelling errors were accounted for. Similar problems can be found in the large separation HLEM example. In this case the shifts caused by the primary field errors are much smaller however, in order to achieve the target misfit it is still necessary to account for modelling errors in the inphase component. In order to simulate actual HLEM measurements Gaussian noise with a constant standard deviation of 0.5% of the primary field was added to the data from each of the large separation HLEM soundings. Then the noisy data from the three soundings were each inverted using a coil separation of s = 50 m even though the true coil separations for each of the soundings were actually 51 m, 50 m, and 49 m. The reference model for each inversion was 0.01 S/m. The error estimates assigned to the inphase data were 8% of the primary field for both strUe = 51 m and 49 m and the target misfit was set to 20. The final data misfit and model norm values from the three inversion are shown in Table 3.4 and the recovered conductivity models are plotted in Figures 3.11(a),(c), and(e) and the observed and predicted data are plotted in Figures 3.11(b), (d), and (f). In all three inversion results the quadrature component of the data appears to have been predicted quite well and while the shape of the inphase response has been predicted well, only the 5 t r u e data looks to have been reproduced completely. I notice from Table 3.4 Chapter 3. Geometric Survey Parameter Errors 59 -0.5i Sir = 51 m. •Strue = 50 m. Sir = 49 m. (a) (e) -2.5 -0.5. (C) « -2.5! -0.5i -2.5 True Model Recovered Model 50 100 Depth (m) 150 True Model Recovered Model 50 100 Depth (m) 150 True Model Recovered Model 50 100 Depth (m) 150 00 (d) 30 20 10 !»' ?-10| CO 20 <D (/)• co 9-301 -40 o o IP Obs. — IPPred. x x QObs. -- QPred. -501-; 10' 30 20 10 i« E-101 co S-201 cd g-30 -40 10 10 Freguencv (Hz) o o IP Obs. — IPPred. x x QObs. — Q Pred. -50' 10' 10" 10' Frequency (Hz) 0) £ 10 10 Frequency (Hz) 10" 10" Figure 3.11: Results from the inversion of large separation HLEM data with an s value of 50 m when strUe is equal to 51 m, 50 m, and 49 m and the noise was estimated to account for inphase modelling errors. Panels (a), (c), and (e) show the true (dashed line) and recovered (solid line) conductivity models and panels (b), (d), and (f) show the observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature - dashed line) from the three inversions. Chapter 3. Geometric Survey Parameter Errors 60 that the target misfit was achieved for strue = 51 m and 50 m, and the value recovered for •strue = 49 m was close (fa « 40). All of the recovered models shown in Figures 3.11(a), (b), and (c) are good representations of the true model. I have shown that it is possible to achieve the target misfit by assigning large errors to the inphase data. However, when this is done the modelling errors in the inphase data are not being adequately dealt with. This suggests that there may be a better way to deal with coil separation errors. 3.3 Parameter Errors in the Field Data Sets Now that I have discussed the details of measurement height and coil separation errors I should be able to confirm the role, if any, that geometric survey parameter errors play in the field data sets. For both the Mt. Milligan and Sullivan field data sets I will address the possible sources of geometrical survey parameter errors, how the errors affect the data, and how they affect the inversion results. 3.3.1 Mt. Milligan Mt. Milligan is an area with topographical variation and with variable tree cover, both of which can introduce measurement height errors. The conductivity model recovered by inverting the data (Figure 1.6) shows regions of very low conductivity near the surface. These were identified as "air layers" which can arise when hest is smaller than htTUe- It would appear, that at least in some areas, measurement height errors caused by tree cover could be a problem. Chapter 3. Geometric Survey Parameter Errors 61 3.3.2 Sullivan Data The possibility of coil separation errors due to topography at Sullivan was discarded due to the fact that the tailings pond is relatively flat. However, it was discovered that the design of the Max-Min system used to collect the data presented the possibility of a maximum coil separation error of +0.5 m. Since the expected coil separation in the Sullivan survey was equal to 5 m, this represented a possible +10% coil separation error. In Section 3.2 the small separation HLEM example was used to model synthetic data with a +10% relative coil separation error. This error resulted in a shift in the inphase component of the data by approximately -25% of the primary field while there was no evidence of distortion in the quadrature component. When I re-examine the Sullivan data shown in Figure 1.1 I see that the 7040Hz measurements of the inphase component have values between -20% and -30% of the primary field. The model that was recovered by including large error estimates on the inphase component is shown in Figure 1.2. From this plot it is seen that the target misfit of 8 has been achieved in most cases, and the subsurface conductivity agrees with the available a priori information. However, the fact that it is not possible to predict the inphase component of the data (Figure 1.3) leaves something to be desired. It appears that coil separation errors are causing problems with the inversion of the Sullivan data. Chapter 4 Recovering 1-D Conductivity and a Geometric Survey Parameter The previous chapter contained a discussion of the problems that erroneous survey para-meters can cause when inverting both AEM and HLEM data. In this chapter a solution to these problems is proposed in the form of an inversion methodology which can recover both a 1-D conductivity distribution and a geometric survey parameter p. The idea of inverting geophysical data to recover a property distribution as well as information about the experiment has been done previously by Pavlis & Booker (1980) and deGroot-Hedlin (1991). The approach I choose is to include the survey parameter in the model vector and solve the inverse problem using a methodology similar to that used in the absence of the extra parameter. There is no reason that p could not be a vector of parameter values. The formulation and numerical solution that I will present is designed in such a way that this could be accommodated. However, in my case a single parameter is all that is needed to correct the field data sets and hence I will treat pas a single parameter from now on. I will begin with a discussion of the changes that are necessary when dealing with an extra parameter. These include changes to the model, the data, and the sensitivity calculations. This will be followed by the development of a modified objective function, and a discussion of some of the problems and possible solutions encountered when at-tempting to solve the inverse problem. I will show synthetic results from the AEM and HLEM examples and apply the new inversion methodology to the two field data sets. 62 Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 63 4.1 Changes due to the extra parameter In this section I will discuss how the model, the forward modelling and the sensitivity calculations are affected by the inclusion of an extra parameter. When referring specif-ically to details that relate to either AEM or HLEM data I will use a superscript A or a superscript H respectively (e.g. dA will represent the forward modeled AEM data). How-ever, outside of the discussion in this section, the methodology will not specify between the inversion of AEM and HLEM data. Each inversion will be carried out in terms of generic models, data, and sensitivities. 4.1.1 The Model To deal with survey parameter errors I have included an erroneous parameter p in the model vector I am attempting to recover. The generic model m now has the form m = [log^x),... ,log(aM),p]T, (4.1) where M is equal to the number of layers in the conductivity model. The vector m has a length of M + 1. The parameter is either measurement height or coil separation depending upon whether the data are AEM or HLEM. As a result the two specific models which must be denned in order to deal with parameter errors are mA = [log{<Ti),log(<rM), h]T, (4.2) which is used when inverting AEM data, and mH = [log(<Ti),log((TM), s]T. which is used when inverting HLEM data. (4.3) Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 64 4.1.2 The Data As in Section 2.2.2 the data values are stored in the vector d which has the form d = [iZe(di), Imidr),Re(dN), Im(dN)]T, (4.4) where 7Y is equal to the number of data. Both AEM and HLEM data, as usually denned, are normalized by an estimate of the primary field value at the receiver and multiplied by a constant. For my purpose it is simpler to work with the denormalized data. Therefore, I define the generic form of the ith measurement in a given data set to be Re(di) = Re(Fi[m}), (4.5) and Im(di) = Im(Fi[m]), (4.6) where the details of the forward modelling F{ depend upon the survey being considered. In the AEM case I denormalize the data from Equations 3.5 and 3.6 and substitute h for htrue to get the following expression for the data from the ith AEM measurement Re{di)^Re{Hsz{^Si;a,h)\ (4.7) and Im(df) = Im{Hsz{Ui, S i ] a, h)). (4.8) Since the estimated measurement height hest has no direct effect upon AEM data, the inclusion of htrue does not change the forward modelling calculations. Similarly, by denormalizing Equations 3.11 and 3.12 substituting s for strUe the data from the ith HLEM measurement will have the form Re(d?) = Re{Hl{^ h,-a, s)) + HP(s) - HP(sett), (4.9) Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 65 and Im(df) = Im(Hsz(ui,hi] cr, s)). (4.10) From Equation 4.9 it is clear that the presence of a shift in the inphase data is due to a coil separation error. 4.1.3 The Sensitivities In order to solve the inverse problem it is necessary to be able to calculate the sensitivities of the data to changes in the model parameters. These quantities are expressed in the form of the Jacobian, or sensitivity, matrix J. The ijth entry in J represents the sensitivity of the ith datum to a change in the jth model parameter as defined by As a result the generic sensitivity matrix will have the form J=[J„JP], (4.12) where is the JV x M conductivity sensitivity matrix and Jp is the JV x 1 parameter sensitivity matrix. From the expressions for both dA and dH it is clear that the subsurface conductivity structure only affects the secondary field terms in AEM and HLEM data. Therefore, regardless of type of data being inverted it is necessary to calculate the sensitivities of the secondary fields to changes in conductivity. The sensitivities for the type of frequency domain EM surveys I am dealing with have been determined in Zhang & Oldenburg (1994) using the adjoint Green's function method. The sensitivity of the ith complex datum to a change in the conductivity of the jth layer can be expressed as ddi - 2 T T d(Tj iuipomx J™ 0T+1 ^ dz) j°(A5)A3dA' (4-13) Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 66 where Zj is the depth to the top of the jth layer, is the tangential component of the electric field in the Hankel transform domain, and J 0 is a zeroth order Bessel function. However, the model parameter I am concerned with is log(o-j) therefore, the sensitivity value I want is denned by J« - - jr*. (4-14) o-j do-j These values will fill Ja which makes up the first M columns of the sensitivity matrix. The final column of the matrix will contain the sensitivities of the data to the survey parameter p. The details of these sensitivities depend upon the survey being considered. In the AEM case I must determine the sensitivity of the data to a change in measure-ment height h. Since the AEM data in Equations 4.7 and 4.8 are simply the secondary magnetic field we can use Equation 2.37 to express dA in complex form as d? = -^ RTE\3J0{\Si)d\. (4.15) 4 T T J0 U0 By differentiating Equation 4.15 with respect to h I can calculate the sensitivity of the ith datum to a change in measurement height. The complex sensitivity is equal to ^ = ^ 1°° e~2U0h W J0(A«)«tt. (4-16) These values will fill Jp when inverting AEM data. When dealing with HLEM data I want to calculate the sensitivity to changes in the coil separation s. I can expand the HLEM data in Equations 4.9 and 4.10 in complex form as * - = f ^ W * ( A . ) « + - ( ^ ) . (4,T) The sensitivity of the ith datum due to a change in the true coil separation can be calculated by differentiating Equation 4.17 with respect to s as follows dd? mT r e-2uohi n , 3 d , T „ .xx „ mT d ( 1 ds Air Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 67 Expanding the partial derivative in the integrand I get = £ [ J - i ( A « ) - MXs)}. (4.19) Using the identity J-n(As) (-1)" J„(A«), (4.20) I am left with dJo(Xs) ds AJi(Aa), (4.21) which can be substituted back into the integrand to give the complex sensitivity value These values will fill Jp when inverting HLEM data. 4.1.4 The Objective Function Since I have modified the model it is necessary to modify the objective function accord-ingly. I want the solution of the modified problem to remain as close to the parameter error free version as possible. Therefore, I choose a global objective function (f> similar to the one defined in Section 2.2.2. The objective function has the form where <pd is the data misfit, <pm is the model objective function, and /3 is the regularization parameter. I will choose c6j to be the same as the misfit term denned in Equation 2.46. However, since the model has changed I have to redefine <pm as + (4.22) <P = <Pd + Pfa (4.23) (4.24) Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 68 where 4>a is the conductivity model norm, (pp is the parameter norm, and ap is weighting parameter which defines the relative importance of <pp within the model norm. The con-ductivity norm determines the character of the conductivity structure I want to recover, and is defined as 4>« = \\Wv{log(a) - log(aTef))\\\ (4.25) where Wa is the conductivity weighting matrix and <rre/ is the reference conductivity model. Wg- is chosen to be identical to the model weighting matrix defined in Equa-tion 2.43. The parameter norm determines the closeness of the recovered parameter to some reference value pref, a n d is defined as ^={P-Preff. (4.26) Using the model m from Equation 4.1 it is possible to express <6m as (pm = \\Wm{m-mref)\\2, (4.27) where Wm is equal to Wa 0 0 , / a (4.28) Having defined the new objective function it is now possible to solve the inverse problem. 4.2 Inversion Algorithms Now that I have made all of the changes needed to deal with the extra parameter it is possible to design my inversion algorithm. This thesis, so far, has avoided discussing the Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 69 details of my inversion algorithm. Therefore, I begin this section with a description of the algorithm used to invert EM data when p is fixed. The modifications to the algorithm that are necessary in order to accommodate the extra parameter will depend upon my definition of fa. Hence, I will examine the best possible choice of <j)p and how this choice affects the performance of the current inversion algorithm. The insight gained through these investigations is used to design the algorithm which can invert EM data to recover both conductivity and a geometric survey parameter. 4.2.1 Details of algorithm to recover 1-D conductivity with fixed p In order to provide as much detail as possible I will present a description of the two inversion algorithms that can be used to recover a 1-D conductivity structure with a fixed value of p. The first algorithm is a detailed description of the method provided in Section 2.2.2 to solve the non-linear inverse problem for a particular f3 value. It was stated that this process could be repeated for a number of f3 values in order to find the one which satisfies fa — 4>d- In practice however, this is computationally expensive. Therefore, I will describe a second algorithm in which a j3 is varied at each iteration, such that the eventual solution will satisfy fa = <f)d. These two algorithms will be referred to respectively as the fixed /3 and discrepancy principle algorithms. Details of fixed (3 Algorithm This section will present the details of the algorithm I use to invert EM data with a fixed value of (3. As mentioned in Section 2.2.2 the inverse problem can be solved by finding the model m that minimizes (j>. This is done using a damped Gauss-Newton (DGN) approach. The theory of this approach can be found in Dennis and Schnabel (1996) and some background on its application to non-linear inverse problems can be found in Haber (1997). In my algorithm the DGN approach proceeds as follows: at each iteration the Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 70 Gauss-Newton step 8m is checked to see whether or not it has actually decreased (j). This done by calculating the objective function corresponding to the updated model <f>(mk+1) = cj>kd+1+f3<f>^1, , (4.29) and the objective function corresponding to the model from the previous iteration <Kmk) = 4>hd + r3<ft. (4.30) If <j)(mk+i) < </>(mk) then 8m is accepted since it decreases the global objective function. However, when </)(mk+i) > (f>(mk) the step 8m should not be taken. Therefore, the step is corrected by damping it as follows 8m = u>8m, (4-31) where a; is a constant between 0 and 1. In my algorithm u> = 0.5. After 8m is damped, the model mk+i is updated and <f>(mk+i) is recalculated and compared to <f){mk). This process is repeated until a 8m, that decreases the objective function, is found. When an acceptable mk+i is found the algorithm must decide whether to continue on to the next iteration or to accept the model mk+i as the solution to the minimization problem. The stopping criteria used in the algorithm are (1) the relative gradient of <j>, i.e. is mk+i at a minimum, and (2) the relative change between mk+i and mk, i.e. is the model stationary. The relative gradient of <f> is considered small when \gW(mk+1) m^+1\ max i = l,...,M <eu (4.32) <£("ife+i) where g(mk+i) is the gradient of (j> and t\ is a tolerance level. The model is considered to be stationary when | |m f c + 1 - mk\\ max{\\mk+i\\, ||7nfe||) < e 2 , Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 71 where e2 is a tolerance level. The algorithm is said to have converged when either of these conditions are satisfied. In my algorithm both of the tolerance levels are set equal to 10 - 3. With all of these pieces explained in detail it is now possible to map out the entire algorithm. Each of the steps taken during the inversion is summarized in the flow chart shown in Figure 4.1 and each entry in the following list corresponds to a step in the flow chart. While some of the steps may seem trivial for the fixed fi algorithm, they are included such that the numbering scheme will be consistent for subsequent algorithms. 1. Input pertinent information for the inverse problem. This includes: (1) parameters that concern the data such as forward modelling details (measurement frequencies, measurement heights, coil separations, etc.) and the observed data, (2) parameters which concern the model such as discretization details (number of layers, depth to layer interfaces), m,taPt, and mref, and (3) parameters needed to define the objective functions such as a3 and az, and error estimates for the data. The value of fi must also be provided. 2. The initial values of d, <pd, <pm, and J are calculated using m,tart. 3. Iteration counter k is set equal to zero and the variables calculated in Step 2. are labeled to belong to the kth iteration as m.fe, dh, <pd, (pm, and J(mfe). 4. Begin k + 1th iteration. 5. The value of fi supplied in Step 1. is selected as the regularization parameter. 6. The Gauss-Newton step 8m is calculated using Equation 2.60. The model for this iteration is defined as m-k+i = m-k + 8m. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter Fixed P Inversion Algorithm for fixed value of p ( | A ) Use D G N <;8m = co 8m and m , k+n=m (k ,+ 8m > °A Set k=k+1 Update variables User Input v including P l " I (2) Calc initial values I L3J Set k = 0 Initialize Variables) • ( (4) Begin iter, k+1 ) i ' ^5 .Use fixed value of P I ;, ^  Calc G N step 8m Setm ( k + 1 )=m0 0+8m I Calc k+1 values N1 (T) Has § decreased? -N 9 Convergence ? 10 Ex i t Figure 4.1: Flowchart for the fixed (3 inversion algorithm with a fixed value of Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 73 7. The values that were calculated in Step 2. are calculated using mk+\, to be dk+1, <f>k+1, 4>k+\ and J(mk+l). 8. The values of (f>(mk+i) and <f>(mk) are calculated using Equations 4.29 and 4.30. If 8m decreases <f>, then go to Step 9. However, if 8m results in an increase in <f> go to Step 8A and damp 8m, then go back to Step 7. 9. Check for convergence. If solution has converged go to Step 10. Otherwise go to Step 9A and set k = k + 1 and switch the index of all values then go to Step 4. 10. Exit algorithm. Return the recovered model, predicted data, and the final data misfit and model norm values. This algorithm was used to invert the AEM sounding data with a fixed value of (3. The value of /? was selected, through trial and error, such that the goal fa = <f>d was achieved. The measurement height was set equal to 30 m and the rest of the details of the inversion were the same as those used in previous AEM inversions. Figure 4.2 shows the recovered model and the predicted data from the inversion and Figure 4.3 shows the values of fa, <f)m, (3 and <f>, as a function of iteration. It is shown that at each iteration the objective function was decreased until it reached a minimum value. Since /3 had been properly selected the solution achieves the target misfit of 20. However, when the proper {3 is not known it can be time consuming to solve a number of problems in order to find it. Therefore, next step is to construct an algorithm that will find the desired (3 value. Details of discrepancy principle Algorithm The desired algorithm will find the f3 value that satisfies fa = <j>*d. The best plan is to decrease (3 gradually such that the addition of unwanted structure to the model is avoided. Therefore, the algorithm will start with a large /3 value and decrease it slightly Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 74 (a) -3 True Model Recovered Model (b) t 3000 ^2500 C L C L '2000 1500 T3 CO o -o c to gioool CQ SI Q. £ 500 o o IP Obs. — IP Pred. / X X X QObs. / / Q Pred. S  X X' 50 100 Depth (m) 150 10' 10" 10' Frequency (Hz) 10° Figure 4.2: Recovered models and predicted data from AEM example fixed (3 inversion with a fixed value of h. The measurement height was set equal to 30 m. (a) True conduct-ivity (dashed line) and recovered conductivity (solid line) and (b) observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature -dashed line). at each iteration until the desired (3 value is reached. The next step is to decide how to decrease f3. The method that is adopted is referred to as the discrepancy principle. The following is the inversion algorithm I use to invert EM data using the discrepancy principle for a fixed value of p. The algorithm proceeds such that at iteration k + 1 the value (3k+1 is chosen such that it will result in decreasing fa by a fixed amount. This is continued until the target for the current iteration (f>*}k+1^ is equal to the target for the inversion (j)d. Therefore, at each iteration the target <f\ik+1^ is defined as (4.33) where 7 is a constant which should have a value between 0.5 and 1.0. In my algorithm 7 is set equal to 0.5. It is important to note that when searching for <f>*}k+1\ the curve of fa versus (3 usually has a parabolic shape. Therefore, there may well exist two values of (3 that give fa = (j>*}k+^• In this case the smaller of the two (3 values should be chosen Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 75 1 0 * l . . . . 1 1 0 ' l • • • 1 1 2 3 4 5 6 1 2 3 4 5 6 Iteration Iteration Figure 4.3: Convergence curves from AEM example fixed /3 inversion with a fixed value of h. The measurement height was set equal to 30 m. (a) (pd (circles) and target 4>d (crosses) vs. iteration, (b) c6m vs. iteration, (c) fixed f3 vs. iteration, and (d) (p(mk+i) (circles) and <p(mk) (crosses) vs. iteration. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 76 as fik+1 such that the least amount of structure is added to the model. In the case that it is not possible to find a fi value which satisfies (pd = <p*jk+1^ the value of fi which gives the minimum possible <j>d value is chosen as fik+1. After having found fik+1 the Gauss-Newton step Sm is calculated using Equation 2.60 and the model is updated to be mk+i = m-k + Sm. At this point it is necessary to check if the step Sm is beneficial. Since fi is changing at each iteration, the <p values from iteration k, evaluated with fik, can not be compared with those from iteration k + 1, evaluated with fik+1. Therefore, in order to evaluate the changes taking place in cp at successive iterations both (p(mk) and <p(mk+i) are calculated with fi = fik+1 such that the values are defined as 4>(mk) = cMmfc) + /3 f c + Vm (m f c ) , (4.34) and <j>{mk+l) = (pd{mk+1) + fik+1<pm(mk+1). (4.35) In practice, it has been shown that the choice of fik+1, corresponding to a small decrease in <pd, can result in a solution that both minimizes the global objective function and satisfies (pd = (p*d (Parker, 1994). However, this is not always the case. When (p*d is set too small then the target misfit is not attainable and the algorithm can reach a stage at which a decrease in (pd can only be achieved by choosing small fi values that add large amounts of structure to the model. This can cause the recovered model to be quite poor. The problem comes from the fact that there is no guarantee that the step Sm will provide a new model mk+i that is better than the model mk from the previous iteration. Requiring c6(mfe+1) < 4>(mk) is a consistency check on the algorithm. It at least ensures that the previous model is not a better model than the updated model mk+i — mk + Sm. If <p(mk+i) > (p(mk), then a number of actions are possible. First the Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 77 size of the perturbation step may be damped to see if <f)(mk + u)8m) < 4>{mk) for some factor UJ. This did not seem to be necessary here since the step sizes have already been limited by carrying out successive iterations with the goal of slowly reducing the misfit (<£*(fc+1) = 14>{dk)). Empirically, the case that (j)(mk+i) > 4>{mk) occured when the misfit was beginning to plateau out and could no longer be significantly reduced, even with a large perturbation 8m. The achieved misfit </>mm was inferred to be a representation of the minimum misfit. Accordingly (j)d should Ue above this value. Therefore, when the algorithm detects that <^(mj.+i) > <f>(mk), the value of <j)d is increased and the inversion is restarted. In my algorithm the restart value of <f>d is set equal to 1.1 x <^ >mm. The convergence criteria for the discrepancy principle algorithm are the same as for the fixed f3 case. With all of these pieces explained in detail it is now possible to map out the entire algorithm. Each of the steps taken during the inversion is summarized in the flow chart shown in Figure 4.4 and each entry in the following list corresponds to a step in the flow chart. 1. Input pertinent information for the inverse problem. This includes: (1) parameters that concern the data such as forward modelling details (measurement frequencies, measurement heights, coil separations, etc.) and the observed data, (2) parameters which concern the model such as discretization details (number of layers, depth to layer interfaces), mBtart, and mref, and (3) parameters needed to define the objective functions such as a, and ctz, and error estimates for the data. The target misfit <j>d must be provided. 2. The initial values of d, fa, (f>m, and J are calculated using m,taTt. 3. Iteration counter k is set equal to zero and the variables calculated in Step 2. are labeled to belong to the kth iteration as mk, dk, <j>d, (fa, and J(mfe). Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 78 Discrepancy Principle Inversion Algorithm for fixed value of p J5A^  Choose larger <()* and restart 9 A Setk=k+1 Update variables LL' User Input I Calc initial values J JT; Set k = 0 Initialize Variables 5} Begin iter, k+1 Set target cj)d for this iter. (|>r+,)=max(Y <$>,((>*) ® Line Search for fi (k+, ) J I r~T\ Calc GN step 8m w ' Set m(k+1)= m(k)+ 8m I ! (Z) Calc k+1 values ) <jjf) Has <|) decreased? i>[ C+1)+P(k+1)C)<^)+B(k+1) .00 Convergence ? 1 10 Exit Figure 4.4: Flowchart for the discrepancy principle inversion algorithm with a fixed value of p. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 79 4. Begin k + 1th iteration. Calculate (f>d +1' using Equation 4.33. 5. The f3 value which satisfies fa = <j)*}k+1^ is selected as /3 f e + 1. 6. The Gauss-Newton step 8m, corresponding to /3 f e + 1, is calculated using Equa-tion 2.60. The model for this iteration is defined as mk+i = mk + 8m. 7. The values that were calculated in Step 2. are calculated using mk+i, to be dk+1, <f>kd+1, <fft\ and J(m f c + 1). 8. The values of 0(mfc+1) and <j>(mk) are calculated using Equations 4.34 and 4.35. If the step 8m decreases </> go to Step 9. If, however, 8m results in an increase in (f> go to Step 8A and increase the global target misfit <f>d. Then the inversion is restarted by going to step 3. 9. Check for convergence. If solution has converged go to Step 10. Otherwise go to Step 9A and set k — k + 1 and switch the index of all values then go to Step 4. 10. Exit algorithm. Return the recovered model, predicted data, and final data misfit and model norm values. This algorithm was used to invert the AEM sounding data using the discrepancy principle. The target misfit was set equal to 20. The measurement height was set equal to 30 m and the rest of the details of the inversion were the same as those used in previous AEM inversions. Figure 4.5 shows the recovered model and the predicted data from the inversion and Figure 4.6 shows the values of fa, (j}m, /3 and <j>, as a function of iteration. The inversion results are identical to those presented in Section 3.1 when the AEM data were inverted with the correct h value. Figure 4.6(a) shows that after each iteration the data misfit was decreased by a constant amount until reaching the target value of 20. It Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 80 ( a) ^ E > o 3 -n c o O o F T — True Model — Recovered Model (b) Depth (m) 150 10" 10 Frequency (Hz) Figure 4.5: Recovered models and predicted data from AEM example discrepancy prin-ciple inversion with a fixed value of h. The measurement height was set equal to 30 m. (a) True conductivity (dashed line) and recovered conductivity (solid fine) and (b) observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature - dashed line). can also be seen in Figure 4.6(b) that when the target level of the misfit was attained at iteration 8, the model norm was then decreased slightly until the solution converged. This results in a model which both fits the data and has minimum structure. The plot of (f> in Figure 4.6(d) shows that at each iteration the objective function was decreased. In order to define the details of the modified inversion algorithms it is necessary to define the parameter norm. 4.2.2 Defining the Parameter Norm The new objective function, which includes a geometric survey parameter, requires the specification of ctp and pref- A discussion of the ways in which each of these values affect the problem is presented in this section. The goal of this discussion is to determine an educated way in which to assign values to both ctp and pref- The discussions in this section will deal strictly with the discrepancy principle inversion algorithm. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 81 Figure 4.6: Convergence curves from AEM example discrepancy principle inversion with a fixed value of h. The measurement height was set equal to 30 m. (a) fa (circles) and target <j>d (crosses) vs. iteration, (b) c/>m vs. iteration, (c) /3 vs. iteration, and (d) </>(mfe+i) (circles) and (f>(mk) (crosses) vs. iteration. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 82 Choosing <xp In order to choose a good ap value it is necessary to understand how this choice affects the problem. The value of otp determines the relative importance of the parameter norm term. Therefore, it controls how strictly the requirement that the recovered p is close to pref will be enforced. To get a feeling for the effects I set ap equal to its two extreme values: 0 and oo. As ap approaches infinity, p is fixed at the value of pref. However, when ap is equal to zero, p is not required to be close to the pref value at all. Looking at how the problem behaves at these two extremes provides insight into choosing the value of ap. I begin by looking at the case where ap approaches infinity. Since it is not mathem-atically possible to set ap = oo I simply fix p = pref for each inversion. For each of the synthetic examples I carried out a series of inversions, each with a different pref value. For each inversion I attempted to fit the data to a specified target misfit. The ability of the algorithm to achieve the desired misfit level, as well as the resultant (f>m value, provided insight into the importance of p in both AEM and HLEM inversions. The first example I looked at was the noisy AEM data with htrue = 30 m. The data were inverted 11 times with h fixed equal to different values of href ranging from 45 m to 15 m. For each of the inversions the data were assigned an error of 5% of the data amplitude + 10 PPM, the target misfit was equal to 20, and a reference model of 0.01 S/m was used. The results from all of the inversions are plotted together in Figure 4.7. The values plotted in the top panel are the final <6m values from each inversion. These values are equal to <$>a since <6P is equal to zero as a result of p being set equal to pref. The middle panel shows the final c/>d values from each inversion and the bottom panel shows the recovered conductivity models plotted next to one another. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 83 45 42 39 36 33 30 27 24 21 18 15 log 1 Q o Measurement Height (m) Figure 4.7: Results from eleven inversions of the AEM sounding data with fixed values of h ranging from 45 m to 15 m. The top panel is final c/>m value, middle panel is final <f>d value, and the bottom panel is the recovered conductivity models plotted next to one another. These results are identical to those described in Section 3.1. When the h value used by the inversion is smaller than htrUe, the recovered model has a low conductivity structure at the surface and it is possible to achieve the target misfit level. When the h value used by the inversion is larger than htrue, the recovered model has a high conductivity structure at the surface, and, as h increases, it becomes impossible to fit the data. The success of any inversion is gauged, in part, on the ability of the recovered model to predict the observed data. Figure 4.7 shows that there are a number of h values that satisfy this condition. By considering only the (f)m values that correspond to this subset, h — 33 m to 15 m, it can be seen that there exists a minimum (f>m value that corresponds to a particular value of h. For this example, this value of h is approximately equal to 27 m which is close to htrUe-I also inverted the noisy data from the HLEM examples. From the small separation Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 84 §"100 150 12.5 12 11.5 11 10.5 10 9.5 9 8.5 8 7.5 logger Coil Separation (m) Figure 4.8: Results from eleven inversions of the small separation HLEM sounding data, •Strue = 11 rn, with fixed values of s ranging from 12.5 m to 7.5 m. The top panel is final (j)m value, middle panel is final fa value, and the bottom panel is the recovered conductivity models plotted next to one another. case I chose to invert the data that had an strUe value of 11 m. Eleven inversions were carried out with sref values ranging from 12.5 m to 7.5 m. For each of the inversions the data were assigned an error of 0.5% of the primary field, the target misfit was equal to 20, and a reference model of 0.01 S/m was used. The results are plotted next to one another in Figure 4.8. The middle panel shows that the data could only be fit to the desired misfit level when s was equal to 11 m. This value is near a distinct minimum in the versus s curve. Both the model norm values and the recovered conductivity structures from the other s values provide little information since it was not possible to fit the data. The inversion results from the large separation HLEM example show the same effect. The noisy data with strue = 51 m were inverted eleven times with sTef values ranging from 52.5 m to 47.5 m. Again the data were assigned an error of 0.5% of the primary Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 85 Figure 4.9: Results from eleven inversions of the large separation HLEM sounding data, •Strue = 51 m, with fixed values of 5 ranging from 52.5 m to 47.5 m. The top panel is final 4>m value, middle panel is final value, and the bottom panel is the recovered conductivity models plotted next to one another. field, the target misfit was equal to 20, and a reference model of 0.01 S/m was used. Figure 4.9 shows that the results are similar to the small separation case. The target data misfit was only achieved when s was set equal to 51 m and the inversion results from every other s are inconclusive. These HELM examples show that there exists a narrow range of p values for which a conductivity model can be found that will satisfy the data. In both of these examples the values of s were close to strue. Each of the three examples suggest that if p were allowed to vary freely it would choose a satisfactory value. In the AEM example, when there is a wide range of h values that can achieve the target misfit, the selection of that p that provides a minimum model norm value is a good estimate of the true value. In both of the HLEM examples a narrow range of s values that could achieve the target misfit were found. Therefore, it appears that setting a p = 0 is a good choice. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 86 Instabilities when ocp = 0 While it seems that setting ap = 0 is a good idea it can cause problems during the first few iterations. As described in Section 4.2.1 the goal of the inversion algorithm is to add structure to the model slowly. Therefore, the starting model is usually set equal to the reference model. This means that in the first iteration, as fi —> oo, the recovered model reverts to the reference model. Hence, both <f>d and 4>m return to their starting values. By decreasing fi slowly the problem is regularized and any poor steps that may have been taken are avoided. However, if ap is set equal to zero, 4>p will not be regularized by fi and there is no guarantee that the step 8p will be beneficial. By examining the system of equations that are solved to give the Gauss-Newton step it is possible to see what happens when ap = 0. In order to do this I will redefine the terms in Equation 2.60 in accordance with changes made in Section 4.1. The model perturbation Sm becomes, Sm = [Slog(a), Sp]. (4.36) For ease of presentation I will assume that Wd is equal to the identity matrix. As a result JTWjWdJ becomes JTJ which can denned explicitly, using Equation 4.12, as JTJ r JTJ jTj Jo- J°~ Jo~  JP </<7 Jp — L y . Jp °<r u p °p (4.37) where the matrix has a block format, and the upper left hand block of the matrix is M x M, the upper right hand block is 1 x M , lower left hand block is M x 1, and the lower right hand block is 1 x 1. From Equation 4.1.4 the matrix W^Wm, is defined as WjWa 0 0 (4.38) Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 87 Substituting all of these details into Equation 2.60 we can write out the Gauss-Newton step in matrix form as Slog(<r) | Sp jjj. J J J P \ +/3(w?w<, 0 jTj jTj °p °<r up °pj 'jf Sd-j3 0 aE WfWa 0 0 a„ log(a) - log(aref) V - Vref (4.39) where Sd is equal to d°b* — F[m]. For ap = 0 and /3 —> oo Equation 4.39 becomes WjW, 0 0 J P % Slog(a) Sp -W?Wa {log(o-) - log(aref)) JTSd (4.40) such that the conductivity and parameter perturbations are decoupled and the solution has the form Slog(o-) = - {log(a) - log(aref)) h = ( J P T J p ) " 1 J j w . (4.41) (4.42) From Equation 4.41 it is clear that Sa is equal to zero when log(cr) = log(o~ref). The parameter perturbation however, is equal to the Gauss-Newton minimization of fa with respect to p which is independent of (3. This clearly shows that fa is not being regularized at all. Using the AEM sounding data the possible problems associated with settling ap = 0 can be shown. For this example I chose an href value of 42 m and denned the rest of the details governing the inversion such that they were identical to those used when the data was inverted with a fixed h value in Section 3.1. The first iteration of the inversion was then performed for both the ap = oo and the ap = 0 cases. Figure 4.10 shows plots of the values of fa, fa, <f>, and h evaluated at a series of (3 values during the search for c^ 1 .^ Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 88 Figure 4.10: Values calculated during the line search for fi during the first iteration of the inversion of the AEM data. The starting height was set equal to 42 m. The two fines in the plot correspond to values of a p = oo (circles) and 0 (plusses). (a) <f>& vs. fi, (b) vs. (3, (c) h vs. f3, and (d) (f> vs. fi. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 89 The solid lines with circles refer to the case where a D = oo and the dashed lines with plusses represent ctp = 0. When ctp = oo the problem behaves as expected. At large (3 values the fa curve asymptotes to a constant value. By taking note of the fa versus (3 curve, which goes to zero as f3 —> oo, and the fact that h is fixed, it is evident that as /3 gets large m is reverting to mref and therefore, fa is approaching its starting value. However, setting ap = 0 causes a change in the behavior. The fa versus (3 curve still approaches zero as (3 increases however, h plateaus to a value approximately equal to 20 m. This causes the <j>d values, and hence the <f> values, at large (3, to be larger than their initial values (which are known from the ctp = 0 curves). Therefore, the linearized Gauss Newton step Sp is too large and it is necessary to introduce some additional regularization to ensure that the step is acceptable. A simple solution is to set ap equal to some non zero value. When ap is non zero and /3 —> oo the Gauss-Newton step will have the form WjW, 0 0 h and the solution will be Slog(a) '-W^Wa(log{a)-log(<rref)) Sp -(P-Pref) (4.43) Slog(cr) = - (log(o-) - log(aref)) Sp - -(p-pref), (4.44) (4.45) such that when log(o~) = log(o~ref) and p = pref, Sp will be equal to zero. It is important to note however, that while Slog(a) goes to zero as a function of (3, Sp goes to zero as a function of the product f3ap. Therefore, any non-zero value of ap will result in Sp —> 0 as f3 gets sufficiently large, but the point at which this occurs depends upon ap. This is shown by the plots of fa, fa, <j), and h versus /3 from the first iteration of the AEM sounding inversion shown in Figure 4.11. The solid lines with circles refer to the case Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 90 Figure 4.11: Values calculated during the line search for fi during the first iteration of the inversion of the AEM data. The starting height was set equal to 42 m. The three lines in the plot correspond to values of a p = oo (circles), 10 - 6 (plusses), and 10 - 3 (crosses), (a) 4>d vs. fi, (b) 4><r v s - fit (c) h vs. fi, and (d) <f> vs. fi. where otp = oo, the dashed lines with plusses represent ap — 10~6, and the dotted lines with crosses represent ap = 10 - 3. From the curves in Figure 4.11 it is clear that as ap increases the value of fi at which 8p goes to zero decreases. For the algorithm to be stable it is necessary for ap to be non zero at early iterations. However, it has been shown that ap should be equal to zero in the late stages of the inversion. Therefore, we propose a cooling schedule for ap. This entails decreasing ap gradually at each iteration such that in the final iterations it is equal to zero. This is, in Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 91 fact, a simplistic application of the Levenberg-Marquardt method (Dennis & Schnabel, 1996). Implementing a cooling schedule for ap provides the algorithm with the stability it needs at early iterations as well as permitting the parameter to roam freely at the end of the inversion. The final step needed to finalize the algorithm is to determine a method to choose Pref-Choosing pref During the first few iterations the inversion algorithm will demand that p be close to pref therefore, it is a good idea to have a realistic estimate for its value. The problem of selecting a value for pref has some inherent difficulties associated with it. While, both the AEM and HELM systems provide an estimate of the parameter of interest it is unclear as to whether or not this value should be used as pref- The fact that an inversion algorithm is being used to recover the true parameter value indicates that it has already been assumed that peat may be incorrect. Therefore, it may be useful to find a better estimate of pref-A simple choice of pref is to find the value that minimizes (f>d for the reference con-ductivity model aTej provided by the user. The first step in finding this value of p is to fix the conductivity model equal to <rref and forward model data for a range of p values. Using the data generated from each p value, the data misfit is plotted as a function of p and the value that corresponds to the minimum in the <i>j versus p curve is chosen as p r e/ . This estimate of pTef will depend upon the definition of crref. In the case that the reference conductivity model is a good first order approximation to the true conductivity model, the choice of pref should be accurate. However, when the guess of aTef is not as good the results may vary. I investigated the influence of crref on the choice of pref for each of the three synthetic examples. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 92 1 0 (a) * AEM — a = 0.1 S / m — a = 0.01 S / m a = 0.001 S / m -0.5i (b) I 3 -2 TJ c o 0-2.5 •2 -3 0 20 40 60 Measurement Height (m) -3.5, True o — tj = 0.1 S/m — o = 0.01 S/m o = 0.001 S/m 50 100 Depth (m) 150 Figure 4.12: Results from the search for a best fit pret value for the AEM sounding data where htTUe = 30 m. (a) fa versus h curves for the three halfspaces: 0.1 S/m (solid line), 0.01 S/m (dashed line), and 0.001 S/m (dotted line), and (b) the halfspace models (as in (a)) and the true conductivity structure (thick gray line) . The data misfit values, calculated for the AEM sounding data at a series of h values ranging from 0 m to 60 m, are plotted in Figure 4.12(a). The true conductivity model as well as the corresponding halfspace o~Tef models (0.1 S/m, 0.01 S/m, and 0.001 S/m) are plotted in Figure 4.12(b). It can be seen that when aref = 0.01 S/m (dashed line) the (j)d versus h curve has a broad minimum. However, the minimum does appear to be centered close to the true measurement height of 30 m. On the other hand, when o-rej is varied by one order of magnitude in either direction the results are not as good. In fact, the o~ref — 0.1 S/m (solid line) and the o~ref — 0.001 S/m (dotted line) have minima at approximately 50 m and 11m respectively. This leads to the assumption that the influence of the measurement height on the data can easily be overpowered by the conductivity structure. Therefore, attempting to estimate a good reference value for h in this manner is unreliable. Before addressing a solution to this problem I will present the results from the HLEM examples. Since the influence of s on HLEM data is very strong I would expect the fa versus s curves to be much different than those shown in Figure 4.12. For the small HLEM Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 93 Small HLEM Large HLEM (a) (c) 1 0 10' 2 ra ra Q 1 0 1 0 — o = 0.1 S/m —- o = 0.01 S/m o = 0.001 S/m -0.5, 8 9 10 11 12 Coil Separation (m) — o = 0.1 S/m — o = 0.01 S/m o = 0.001 S/m > N ~ —— 48 49 50 51 52 Coil Separation (m) (b) (d) £-1.51 I -TJ C o 0-2.5 3 -31 -3.5 -0.5i ? "1 55 J-1.5. > o 1 o 0-2.5 CD 2 -31 -3.5, • True o — o = 0.1 S/m — o = 0.01 S/m — o = 0.001 S/m 50 100 Depth (m) 150 * * True o — o = 0.1 S/m — o = 0.01 S/m o = 0.001 S/m 50 100 Depth (m) 150 Figure 4.13: Results from the search for a best fit prej value for the HLEM examples. Small HLEM sounding data where strUe = 11.0 m : (a) fa versus h curves for the three halfspaces: 0.1 S/m (solid line), 0.01 S/m (dashed line), and 0.001 S/m (dotted line), and (b) the halfspace models (as in (a)) and the true conductivity structure (thick gray line). Large HLEM sounding data where 3true = 51.0 m : (c) fa versus h curves for the three halfspaces: (as in (a)), and (b) the halfspace models (as in (a)) and the true conductivity structure (thick gray line). sounding I attempted to estimate s using the noisy data generated with stTue = 11 m. Values of the data misfit were calculated at s values ranging from 7.5 m to 12.5 for conductivity values of 0.1 S/m, 0.01 S/m, and 0.001 S/m. The resulting fa versus s curves are shown in Figure 4.13(a). Each of the fa curves have sharp minima and they are all located at about 11 m. Choosing s which corresponds to the minimum fa value has provided a good estimate of s f r u e in this case due to the fact that the effect of coil separation errors on the data greatly overpowers that of the conductivity structure. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 94 However, as the coil separation increases, and the importance of primary field errors decreases, this method of estimating s becomes less effective. When trying to estimate the value of s for the large separation HLEM data the results are similar. The minima for each curve shown in Figure 4.13(c) are quite broad and the best estimate is made when aref = 0.01 S/m. However, each of the minimum values he within 2 m of the correct separation. By looking at the curves on the fa versus p plots shown in Figures 4.12 and 4.13 indi-vidually, this method of estimation appears to work best when the parameter dominates the response. Otherwise it seems that a poor choice of <rref could result in an incorrect estimate. However, by examining the three curves from each plot together it is possible to see that the best estimate of ptrue not only corresponds to the best crref choice, but also to the minimum with the smallest value. This suggests that it may be possible to get an accurate estimate of pt r u e by searching for the pair of p r e / and aref values that correspond to the minimum value. In order to investigate this possibility the data misfits for a range of <rref halfspace values and p can be calculated and contoured over the <rref and p values. From the preliminary investigation, the AEM data appears to be more sensitive to variation in ground conductivity than the HLEM data. Since the goal is to estimate h the fa value used for the contour plot was calculated using the data from the two frequencies alone. It is hoped that this will limit the influence of the conductivity structure at depth on the data. The contour plot for the AEM example in Figure 4.14 shows fa values for a range of aref halfspaces from 0.001 S/m to 0.1 S/m, and a range of h values from 15 m to 45 m. It can be seen that a minimum fa value exists at approximately 0.01 S/m and 30 m. This minimum lies within a relatively shallow bowl that trends diagonally from low conductivity and small height values to high conductivity and large height values. This Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 95 Coil Separation (m) Figure 4.14: Contour plots of <j>d for a range of t r r e / and p values for the A E M example and both H L E M examples, (a) A E M example in which only two highest frequencies were used to calculate fa, (b) small H L E M example and (c) large H L E M in which every frequency was used to calculate fa Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 96 is in general agreement with the inversion results when p is incorrect from Section 3.1. When the bird is too close to the earth the data can best be predicted by a less conductive structure, and vice versa. This behavior explains the shape and depth of the bowl in fa contour plot. For each of the HLEM examples, all of the data were included in the parametric inversion. This was done because it appears to be impossible to overpower the effects of coil separation errors with incorrect conductivity estimates. The contour plot for the small HLEM example in Figure 4.14(b) shows fa values for a range of <7 r e/ halfspaces from 0.001 S/m to 0.1 S/m, and a range of s values from 7.5 m to 12.5 m. The minimum fa value is located near 0.01 S/m and 11 m. This minimum lies in a deep trough which runs parallel along the s = 11 m value for each conductivity value. This is a clear indication that the true value of s can not be masked by an incorrect <rref value. The results from the large HLEM example are similar. The contour plot for the small HLEM example in Figure 4.14(c) shows fa values for a range of <rre/ halfspaces from 0.001 S/m to 0.1 S/m, and a range of s values from 47.5 m to 52.5 m. The minimum fa value is located near 0.01 S/m and 51 m. The pattern of the contours about the minimum appear more spread out than those from the small HLEM example. This is due to the fact that the distortions caused by coil separation errors are a function of the relative coil separation error. Since the small separation example represents a 10% error, and the large example only represents a 2% error, the difference in the patterns make sense. These examples have shown that accurate estimates of both <rPe/ and pref, can be obtained from either AEM or HLEM data. Therefore, the inclusion of a pre-processing step in which a best-fit 2-parameter model is determined would be beneficial to the inversion algorithm. The inversion could then be carried out using the calculated value of pref. The decision Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 97 to discard the calculated aref value is based on the fact that while the estimate of <rre/ may be a good one, since the aref value provided by the user was most likely chosen for a particular reason. One reason could be that the sounding is part of line of inversions being performed, and for consistency it is best to use the same reference model for each inversion along the line. While the option to choose between using the estimated and supplied pref value could easily be included in an inversion algorithm, for my purposes I have decided to always use the user supplied value of crre/. This section has provided a reliable way in which to choose pTef and a stable way in which to define a p , and hence it is now possible to define the inversion algorithm to recover both conductivity and a parameter. 4.2.3 Inversion algorithm to recover 1-D conductivity and p In order to provide as much detail as possible I will present a description of the two inversion algorithms that can be used to recover a 1-D conductivity structure and a survey parameter p. The first algorithm solves the problem for a fixed /3 value and the second algorithm uses the discrepancy principle. Details of fixed (3 Algorithm This algorithm is very similar to the one presented Section 4.2.1 for the fixed p solution. The only modifications are the inclusion of a cooling schedule for ocp that affects Steps 1 and 9A, and the inclusion of the step in which pref is determined as described in Section 4.2.2. A flow chart of the algorithm is shown in Figure 5.6. The changes in the algorithm are as follows: 1. The user must provide a starting value of ctp. 2. The starting pref value is calculated as discussed in Section 4.2.2. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 98 Fixed (3 Inversion Algorithm with p included in the model ( • C D User Input 2 including otP J Calc initial values including p,er. (3) Setk = 0 Initialize Variables 8A Use DGN 5m = co 8m and m ( U 1 ,= m(k,+ 8m L 9 A Set k=k+l Update variables ocP = 0 a P "•I't 4 ; Begin iter, k+1 ' i ( ^5) Use fixed value of P I Calc GN step 8m w ' Setm ( k + 1 )=m ( k )+8m 7 Calc k+1 values I N ] ($) Has <|> decreased?: T N / ® Convergence ?!;j) 10 Exit Figure 4.15: Flowchart for the fixed (3 inversion algorithm with p included in the model. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 99 (a) True Model Recovered Model (b) £ 3000 ^2500 CL CL '2000 1500 T J CO 3 o T J C CO 8>1000| CO £ 500] o o IP Obs. — IPPred. / x x x QObs. / / ' —- Q Pred. A* •X' j?; 50 100 Depth (m) 150 10 10 Frequency (Hz) 10 Figure 4.16: Recovered models and predicted data from AEM example fixed j3 inversion with h included in the model. The recovered measurement height was equal to 27.2 m. (a) True conductivity (dashed line) and recovered conductivity (solid line) and (b) observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature - dashed line). 9A. The value of ap is decreased by a factor of 10 after each iteration. This algorithm was used to invert the AEM sounding data with a fixed value of j3. The value of /3 was selected such that the goal fa = <f>d was achieved. The value of ap was set equal to 1.0. The rest of the details of the inversion were the same as those used in previous AEM inversions. Figure 4.16 shows the recovered model and the predicted data from the inversion and Figure 4.17 shows the values of <j>d, fa,p, fa, f3 and c/>, as a function of iteration. The convergence curves <j>d, fa, and (f> curves are very similar to those from the fixed p example. However, the important curve is the plot of p versus iteration in Figure 4.17(c). The value of p decreases gradually from it starting value of 30.5 m to eventually plateau at a final value of 27.2 m. While this value is not equal to htrue, it does recover a model that can fit the data, has minimum structure, and reproduces the true model quite well. Therefore, the inversion using fixed f3 has been successful. The next step is to test the methodology using the discrepancy principle. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 100 10 ©-© Recovered Misfit *--* Target Misfit (b) (d) 10 2 4 Iteration (f) o c CD •10' > CD '5 W -Q o CD 10 o-e if new *--* 0 old * ca I. A A A . 2 4 Iteration Figure 4.17: Convergence curves from AEM example fixed (3 inversion with h included in the model, (a) fa (circles) vs. iteration and target fa (crosses), (b) fa vs. iteration, (c) recovered h (circles) vs. iteration and htrtie (crosses), (d) fa v s - iteration, (e) fixed (3, and (f) 4>(mk+i) (circles) and (^mjt) (crosses) vs. iteration. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 101 Details of discrepancy principle Algorithm This algorithm is very similar to the one presented in Section 4.2.1 for the discrepancy principle solution. The only modifications that had to be made were the inclusion of a cooling schedule for ctp that affects Steps 1 and 9A, and the inclusion of the step in which pref is determined as described in Section 4.2.2, which affects Step 2. Since these changes were outlined with respect to the fixed /? algorithm they will not be listed again here. However, a flow chart of the algorithm is shown in Figure 4.18. The discrepancy principle algorithm will be used for the rest of the chapter to invert E M data, henceforth it will simply be referred to as the inversion algorithm. The next step is to test it on the synthetic examples. 4.3 Synthetic Examples In this section the algorithm is tested on each of the synthetic A E M and H L E M examples that have been used throughout the thesis. 4.3.1 A E M Sounding The algorithm was used to invert the A E M sounding data. The value of a p was set equal to 1.0. The rest of the details of the inversion were the same as those used in previous A E M inversions. Figure 4.19 shows the recovered model and the predicted data from the inversion and Figure 4.20 shows the values of fa, fa,p, fa, ft and (f>, as a function of iteration. Figure 4.19 shows that the algorithm has also been able to recover a conductivity structure that is a good representation of the true model. The measurement height recovered using the discrepancy principle has a value of 27.2 m, which is close to the true height of 30 m. The plot of p in Figure 4.20 illustrates the way in which the cooling schedule for a p has ensured that p is not able to wander too far from pref during Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 102 Discrepancy Principle Inversion Algorithm with p included in the model (8A)Choose larger $ and restart v9AJ Setk=k+1 Update variables Op = 0 cxP Ji) User Input including aP -('(2) Calc, initial values . including prer. J ((Jy Setk = 0 ^ I, Initialize Variables ) Begin iter, k+1 Set target ())d for this iter. <j>r+,)=max(YC><iO ^ i f CD Line Search for P(k+1) { ensure 8p not too large I Calc GN step 8m ^ Setm(k+l)=m(k)+8m I 'JD Calc k+1 values j ^ / Q ( (jT) Has (() decreased? j c "+r" c° < $ 7\ X Convergence ? (lOt) Exit Figure 4.18: Flowchart for the discrepancy principle inversion algorithm with p included in the model. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 103 (a) True Model Recovered Model 3000i (b) £ ^2500| Q. Q_ ^2000| CO 3 O "O c (0 $1000| co Q. £ 500 1500 o o IP Obs. — IPPred. x x QObs. -- QPred. 1 X /X . X / 50 100 Depth (m) 150 10' 10' 10" Frequency (Hz) 10* Figure 4.19: Recovered models and predicted data from AEM example discrepancy prin-ciple inversion with h included in the model. The recovered measurement height was equal to 27.2 m and the true measurement height was equal to 30 m. (a) True conductiv-ity (dashed line) and recovered conductivity (solid line) and (b) observed data (inphase -circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature - dashed line). the early iterations when (3 is large. The remaining curves show a character similar to those seen in the discrepancy principle inversions when p was fixed (Figure 4.6). It is satisfying to see the inversion recover both measurement height and conductivity from AEM data has been successful. However, the fact that the recovered height is different from the true height by an amount of 3 m opens up the results to further investigation. In particular, the question of the nature of the trade off which takes place between measurement height and conductivity structure when p is allowed to vary should be explored. 4.3.2 Trade off between a and h In order to determine the nature of the trade off between measurement height and con-ductivity model the inversion results from the fixed p and variable p, inversions of the AEM data will be compared. Figure 4.21(a) shows the recovered models from the fixed Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 104 10 G - Q Recovered Misfit *--» Target Misfit (b) 5 10 Iteration 15 o-e Recdh * - - * True h (d) (f) 5 10 Iteration 10 o <2io3! o « 1 0 ! o O 10 O-o ) new * - • * tiold . * X 8 8 < l -®-® 5 10 15 Iteration Figure 4.20: Convergence curves from AEM example discrepancy principle inversion with h included in the model, (a) <j>d (circles) vs. iteration and target fa (crosses), (b) fa vs. iteration, (c) recovered h (circles) vs. iteration and htrue (crosses), (d) fa vs. iteration, (e) (3 vs. iteration, and (f) (j)(mk+i) (circles) and (f>(m,k) (crosses) vs. iteration. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 105 0 -0.5 (a) -1 E CO > 3-1.5 •o c o O -2 o -2.5 -3. m :p. = 30m rec rtrue m :p = 27.17m rec r rec m, :p. = 30m true rtrue (b) 50 100 Depth (m) 150 -3 —- m :p, = 30m rec rtrue — m :p = 29.88m rec r rec m, :p, =30m true Hrue 50 100 Depth (m) 150 Figure 4.21: Comparison of models recovered from the inversion of AEM data with fixed h (dashed line) and variable h (solid line). The true model is also plotted (dotted line), (a) Models from the synthetic AEM data used throughout the chapter, (b) Models from a different set of AEM in which the model was more complex. h 4>m 4>d Fix h Var. h 30.0 27.2 2.50 x IO"2 2.17 x IO"2 20.0 20.0 20.0 20.0 Table 4.1: Comparison of results from the AEM inversions with fixed h and variable h. h is the recovered measurement height, cim is the recovered model norm, fa is the data misfit achieved by the inversion, and (j)*d is the target misfit. Values are tabulated for both the fixed h (Fix h) and variable h (Var. h) inversions. p (dashed fine) and variable p (solid fine) inversions. The true model is also plotted (dotted fine). The final measurement height, model norm, data misfit, and target misfit values for both of the inversions are shown in Table 4.1. As was mentioned previously, the recovered h value is close but not equal to htrue. By comparing the recovered models plotted in Figure 4.21 it appears that the difference in measurement height is accounted for by differences in the recovered conductivity structure. While the amplitude of the recovered conductive layer at 40 m is similar for both algorithms, the other regions of the structure are different. The model recovered from the variable p inversion appears to have been smoothed out both above and below the peak conductivity. This can be Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 106 confirmed by comparing the <fim values from Table 4.1. These types of distortions should be expected in the AEM case. In Figure 4.7 it was shown that there are a wide range of h values that can fit the data. The algorithm selected the one with the minimum model norm value, which, in this case, corresponded to a height of 27.2 m. The size of the difference between the recovered height and the true height will depend a great deal on the model being considered. In the AEM example shown above, the structure above and below the conductive layer had a secondary influence on the data. Therefore, perturbations to the data due to changes in flight height were compensated by changes to the recovered conductivity structure. In this example the conductivity model was smoothed out and a smaller measurement height was recovered however, a different synthetic model will produce a different result. I will now attempt to invert data collected above a slightly more complex structure to examine the trade off between conductivity and h. In this case a 5 m overburden with a conductivity of 0.1 S/m was placed at the surface and the conductivity of the buried layer was increased to 0.5 S/m. The data were simulated at ten frequencies ranging from 110 Hz to 56320 Hz with a coil separation of 10 m and measurement height of 30 m. This example will be referred to as the overburden AEM example. In order to simulate actual measurements, random noise with a standard deviation equal to 5% of the data amplitude + 10 PPM was added to each datum. Figure 4.21(b) shows the models recovered from the fixed p (dashed line) and variable p (solid line) inversions of the overburden AEM data. The final measurement height, model norm, data misfit, and target misfit values for both of the inversions are shown in Table 4.2. It can be clearly seen that there is very little trade off possible in this example. This is a direct result of the presence of the overburden at surface. In order to fit the data, it is necessary to recover the conductive layer at the surface accurately. Therefore, it is not possible to vary the flight height in order to recover a model with a smaller c/»m value. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 107 h (f>m fa ti Fix h Var. h 30.0 29.9 2.50 x 10"1 2.49 x 10~2 20.0 20.0 20.0 20.0 Table 4.2: Comparison of results from the overburden AEM inversions with fixed h and variable h. h is the recovered measurement height, <f)m is the recovered model norm, fa is the data misfit achieved by the inversion, and <j>d is the target misfit. Values are tabulated for both the fixed h (Fix h) and variable h (Var. h) inversions. The insight provided by these two examples is encouraging. It has been shown that the trade off that can occur between the conductivity and the measurement height will only affect parts of the conductivity structure that are of secondary importance to the data. Therefore, while these trade offs can result in incorrect recovery of the measurement height, they will not distort the important features of the conductivity model. 4.3.3 H L E M sounding This section will present the inversion results from both of the HLEM examples. First of all, the algorithm was used to invert the small HLEM data. The value of ap was set equal to 1.0. The rest of the details of the inversion were the same as those used in previous small HLEM inversions. Figure 4.22 shows the recovered model and the predicted data from the inversion and Figure 4.23 shows the values of fa, fa,p, fa, f3 and <f>, as a function of iteration. Figure 4.22 shows that the algorithm has also been able to recover a conductivity structure that is a good representation of the true model. The coil separation recovered using the discrepancy principle has a value of 10.98 m, which is very close to the true separation of 11 m. The plot of p in Figure 4.23 shows that p never wandered very far from the value determined by the best fit half space. This type of behavior is expected since the HLEM is dominated by the true coil separation value. Similar results are seen when the large HLEM data is inverted. Figure 4.24 shows the recovered model and the predicted data from the inversion and Figure 4.25 shows the Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 108 Depth (m) Frequency (Hz) Figure 4.22: Recovered models and predicted data from small separation HLEM ex-ample discrepancy principle inversion with s included in the model. The recovered coil separation was equal to 10.98 m and the true coil separation was equal to 11 m. (a) True conductivity (dashed line) and recovered conductivity (solid line) and (b) observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature - dashed line). values of (pd, <Pa,P, $P, fi a n d 4>, as a function of iteration. The algorithm has been successful in recovering both coil separation and conductivity for HLEM data. Due to the fact that the data is dominated by the effects of coil separation there is very little trade off between the recovered conductivity structure and the coil separation value. This is as expected since it was shown previously that there exists a narrow range of separation values that will be able to fit the data (Figures 4.8 and 4.9). Now that the algorithm has been tested on synthetic examples it is time to apply it to the field data sets. 4.4 Field Examples I have already determined that the field data sets may be contaminated with geometrical survey parameter errors. Therefore, I will now apply the modified inversion methodology Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 109 (a) 10 co Q 10 o-o Recovered Misfit * - - * Target Misfit 5 10 Iteration 15 (b) 11.01 (c) 1 1 * ] f — - * « * o-o Reed s 10.99! .10.98 *~* True s o o o (d) 10 15 Iteration Iteration Figure 4.23: Convergence curves from small separation HLEM example discrepancy prin-ciple inversion with s included in the model, (a) fa (circles) vs. iteration and target fa (crosses), (b) fa vs. iteration, (c) recovered s (circles) vs. iteration and stTUe (crosses), (d) fa vs. iteration, (e) (5 vs. iteration, and (f) <f>(mk+i) (circles) and 4>{mk) (crosses) vs. iteration. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 110 0 50 100 150 102 103 10" 10s Depth (m) Frequency (Hz) Figure 4.24: Recovered models and predicted data from large separation HLEM example discrepancy principle inversion with s included in the model. The recovered coil sepa-ration was equal to 50.93 m and the true coil separation was equal to 51 m. (a) True conductivity (dashed line) and recovered conductivity (solid line) and (b) observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid Une, quadra-ture - dashed line). to each of the data sets. When inverting the Mt. Milligan data I will attempt to recover the true measurement height and when inverting the Sullivan data I will attempt to recover the true coil separation. 4.4.1 Mt. Milligan When inverting the Milligan data I assigned a standard deviation of 25 Dighem units (5 PPM) to the 900 Hz measurements and a standard deviation of 50 Dighem units (10 PPM) to the 7200 Hz and 56000 Hz measurements. The target data misfit for each sounding was equal to 6. The results from the inversion are plotted in Figure 4.26. The top panel shows the final fa values achieved when h was fixed (solid fine), and the final fa values achieved when h was included in the model (dashed line). The modified inversion algorithm has fit the data to approximately the same degree at each station. This is expected in the case where the bird is too high, which is what was assumed in Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 111 Iteration Iteration Figure 4.25: Convergence curves from large separation HLEM example discrepancy prin-ciple inversion with s included in the model, (a) <pd (circles) vs. iteration and target fa (crosses), (b) fa vs. iteration, (c) recovered s (circles) vs. iteration and 5 t P u e (crosses), (d) fa vs. iteration, (e) fi vs. iteration, and (f) <£(mfc+i) (circles) and 4>(m-k) (crosses) vs. iteration. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 112 ij5 20 0 -100 §"200 I \ • • • • • • • 00 9200 9400 9600 9800 I I •: -3 N o r t h i n g ( m ) iog 1 0 a 3.5 Figure 4.26: Inversion results from Mt. Milligan Data with h included in the model. The top panel: Data misfit values from the fixed h inversion (solid line) and the variable h inversion (dashed). Middle panel: The estimated (solid line) and recovered (dashed line) measurement height. Bottom panel: The recovered conductivity models plotted next to one another. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 113 Northing (m) Figure 4.27: Plot of the recovered Sh superimposed on the fixed h inversion results from Mt. Milligan. The solid fine is Sh = h - he,t. the Mt. Milligan area. The plots of the recovered (dashed fine) and estimated (solid fine) measurement heights in the middle panel of Figure 4.26 indicate the presence of measurement height errors as large as +20 m. However, errors of this size should be expected in an area of dense tree cover such as Mt. Milligan. The conductivity models plotted in the bottom panel of Figure 4.26 also show a big difference from those recovered using the estimated measurement. The "air layers" that were present in Figure 1.6 have now been removed. There is a direct correspondence between the removal of the air layers and the recovery of a measurement height error. The difference between the recovered h values and the measured hest values are, in fact, a good approximation of the thickness of the air layers. In order to illustrate this point I have plotted Sh along the survey fine and have superimposed it upon the conductivity section recovered using the estimated measurement height. This plot is shown in Figure 4.27. Overall, the inversion of the Mt. Milligan data with the modified inversion method-ology has been successful. While the algorithm has not been able to recover any new structures at depth, it has been able to remove the distortions in the model that were caused by the presence of measurement height errors. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 114 4.4.2 Sullivan When inverting the Sullivan data I assigned standard deviations of 1 percent of the primary field to both the inphase and quadrature components of the data. The target data misfit for each sounding was equal to 8. The observed data are plotted with the data predicted by the modified inversion algorithm in Figure 4.28. The improvement, when compared to the predicted data from the fixed s inversion (Figure 1.3), is evident immediately. By allowing s to vary we have been able to predict the inphase component of the data and we have also been able to fit the quadrature data better. The rest of the results from the inversion are plotted in Figure 4.29. The data misfit values plotted in the top panel indicate that the algorithm has been able to fit the data to the desired level at most stations. The middle panel shows the estimated (solid line) and recovered (dashed line) coil separation value. It appears that the coil separation is approximately equal to +0.5 m at most stations along the line. This is in accordance with the error that was estimated by examining the Max-Min system used to collect the data. It is reassuring to see that the results match well with the expectations. The recovered conductivity structure varies slightly at depth from the models recovered from the inversion with fixed s (Figure 1.2). However, the fact that the model can now predict the inphase data improves my confidence in the recovered model. The use of the modified inversion algorithm has been able to improve the results from both the Mt. Milligan and Sullivan data. 4.5 Summary In this Chapter a general methodology was developed to recover both a function and a parameter. It was then applied to the problem of recovering a 1-D conductivity distrib-ution and a survey parameter from EM data. The details of the problem were developed Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 115 ^ 20 (a) ^ W N 0 XL r-_20< O 9=-40. i**xxxx-xxx 50 100 150 200 250 Easting 300 350 IP Obs IP Pre QObs Q Pre (b) d 20 ? 0 x ^ - 2 0 ' O £ - 4 0 0 50 100 150 200 Easting 250 300 350 5? 20 w 1 0 X CVJ-20 O S l - 4 0 0 50 100 150 200 250 300 350 Easting (d) °S 20 N X XL CO io_20 l o i r - 4 0 . 0 50 100 150 200 Easting 250 300 350 Figure 4.28: Observed and predicted data from the inversion of Sullivan data with s included in the model. Observed inphase (circles), predicted inphase (solid line), observed quadrature (crosses), and predicted quadrature (dashed line), (a) Inphase (IP) and quadrature (Q) 7 kHz, (b) IP and Q 14 kHz, (c) IP and Q 28 kHz, and (d) IP and Q 56 kHz. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 116 CO o 1102 CO § 1 0 ° • 1*1 1 t , 1 , - i i , i , - - ; \A _ • Vv_y'- ij >' 0 100 200 Easting (m) 300 l o 9 i o ° Figure 4.29: Inversion results from Sullivan Data with s included in the model. The top panel: Data misfit values (dashed). Middle panel: The estimated (solid line) and recovered (dashed line) measurement height. Bot tom panel: The recovered conductivity models plotted next to one another. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 117 for both AEM and HLEM data. An investigation into effects of adding the parameter to the inversion resulted in a modified inversion algorithm. The algorithm was tested on synthetic AEM and HLEM data sets. The results have shown that the algorithm is stable and provides good results. In the AEM case a good estimate of the measurement height was recovered. The accuracy of the estimate depends in part upon the model being recovered. However, the height that is recovered corresponds to the minimum structure model that fits the data. For the examples shown the conductivity structure that is recovered is a good repres-entation of the true model. When inverting HLEM data the algorithm has recovered an accurate estimate of the coil separation and a good representation of the subsurface conductivity structure. The use of the modified inversion algorithm has been able to improve the results from both the Mt. Milligan and Sullivan data. The inversion results from the Mt. Milligan data show that the distortions due to the "air layers" have been removed by recovering a better estimate of the flight height. This has generated a conductivity model which is much easier to interpret. The inversion of the Sullivan data made it possible to predict the inphase component of the data by recovering an estimate of the true coil separation. While the inclusion of the coil separation in the model does not drastically change the recovered conductivity structure, the plausibility of the inversion result has been increased by the fact that the algorithm has been able predict the inphase component of the data. Chapter 5 Application of GCV to 1-D E M Inversion When inverting geophysical data the presence of noise is always an issue. The amount of information that can be extracted from a noisy set of data is dependent upon our knowledge of the noise. The standard way to deal with noise is to provide an estimate of the noise and attempt to fit the data accordingly. Since it is difficult to estimate the amount of noise associated with the data simply by inspection, this method can cause problems. Generalized cross validation (GCV) is a method of estimating the regularization pa-rameter that will account for the noise associated with a given data set. This technique has been widely used in statistics (Wahba, 1990) and has recently been applied to geo-physical inverse problems (Haber, 1997). This chapter will present a brief review of GCV in both the linear and non-linear inverse problems. Using this methodology, a non-linear GCV inversion algorithm will be developed for the 1-D frequency domain EM inverse problem. Both the problem of when p is fixed and when it is included in the model will be dealt with. Each of the algorithms will be tested on synthetic data sets. Since field data sets, such as the Sullivan and Mt. Milligan data sets, consist of limited numbers of data, the performance of a statistical method such as GCV is uncertain. Therefore, the applicability of GCV to data sets with small numbers of data will also be investigated. 118 Chapter 5. Application of GCV to 1-D EM Inversion 119 5.1 Estimating Noise using GCV This section will begin with a discussion of G C V in a linear inverse problem. After the basics have been covered, a methodology by which G C V can be applied to nonlinear inverse problems will be presented. 5.1.1 GCV in Linear Inverse Problems The idea behind G C V is that a good model should be able to predict measurements that have yet to be taken. However, it is not feasible to take new measurements after each inversion in order to test the recovered model. A simpler approach would be to leave one datum out of the data vector and search for the model that is best able to predict the missing datum. This is the method used by G C V to determine the model which is affected least by any given datum. The theoretical background of G C V can be found in Wahba (1990) and will not be covered here. However, an overview of how G C V can be used in inverse problems will be discussed. I will begin by considering a generic linear inverse problem in which I am attempting to invert the observed data d o 6 i , that are contaminated by random noise with a uniform, but unknown, standard deviation e. I want to recover the model m with the minimum amount of structure that fits dobs to the correct amount. I will begin as usual by denning a global objective function, <f>(m) = fa + fi<j>m = \\Gm - d°b'\\2 + fi\\Wmm\\2 (5.1) where fa is the data misfit, fi is the regularization parameter, <j>m is the model norm, G is a linear forward operator, and Wm is a weighting matrix which penalizes a certain type of structure in the model. As mentioned in Section 2.2.2 the solution of this problem can be found by minimizing <£ at the correct value of fi, such that our expectations for the Chapter 5. Application of GCV to 1-D EM Inversion 120 data misfit are satisfied. GCV is a method which can be used to get an estimate of the correct value of j3. The GCV function is denned as, It has been shown that the value corresponding to the minimum of the GCV function will provide a good estimate for the optimal regularization parameter in problems with the form shown in Equation 5.1. The model recovered using the GCV estimate of /3 sometimes over-fits the data, but generally provides a reliable estimate of the noise associated with When each datum within doha is contaminated by Gaussian noise with the same standard deviation, the GCV estimate of f3 requires no information about the noise from the user. However, in the case where the standard deviation of the noise varies with each datum it is necessary to scale the data such that the relative noise levels are equal for each datum. This is due to the fact that the mathematics of GCV are done in terms of noise with a single standard deviation. However, it is important to note that the user need only provide the relative noise estimates, the GCV estimate of /3 will determine the absolute noise level associated with the data. It is also important to note that the GCV estimate of 0 can not account for correlated noise. This fact is what suggests that the application of GCV to nonlinear problems could be useful. 5.1.2 G C V in Non-linear Problems The non-linear problem is approached in a manner similar to the linear case. The non-linear objective function is denned as GCV{j3) (5.2) [trace(I - G(GTG + r3W^Wm)-1(P')]2' dob'. <t>{m) = cf>d + (3(j>m = \\F[m] - d°b>\\2 + f3\\Wmm\\2, (5.3) Chapter 5. Application of GCV to 1-D EM Inversion 121 where F[m] is a non-linear forward operator. The problem of not knowing fi* still exists, but now the problem is compounded by the fact that the function being minimized is nonlinear. A common way to proceed in a case like this is to linearize the system and solve it iteratively. Therefore, I assume that the current model is equal to mk. My goal now, is to find the perturbation 8m to the model such that mk+\ = mk -\-8m will minimize the objective function <f>. By expanding the operator F in a Taylor series I am left with F(mk+1) = F(mk + Sm) = F(mk) + J(mk)8m + Af(0{Sm2)), (5.4) where J(mk) is the Jacobian defined at mk and N represents the higher order terms of the expansion. Neglecting higher order terms, the data misfit fa can be approximated by the linearized misfit ft™ which has the form falin = \\F(mk) + J(mk)mk+1 - J(mk)mk - doh>\\2. (5.5) This results in a linearized objective function of the form n(m f c + 1) = falin + ficl>m = \\J(mk)mk+1 - rk\\2 + fi\\Wmk+1\\2, (5.6) where rk = doba + J(mk)mk — F(mk). The fact that the linear objective function 4>Un has a form similar to the objective function shown in Equation 5.1 suggests that using GCV to estimate fi might be a good idea. The usual problem of choosing a regularization parameter is made more difficult in the nonlinear case by the fact that the equations being minimized are linear approximations to the real ones. It is not unreasonable to think of the inaccuracies that are introduced by linearization as correlated noise, and if so the GCV estimate of fi might not be affected by these terms. Therefore, GCV should provide an estimate of fi which will be able to fit the noise associated with dobs properly. After selecting fi using GCV, the perturbation to the model is calculated using Equa-tion 2.60. However, the errors due to the linearization may cause problems with the Chapter 5. Application of GCV to 1-D EM Inversion 122 Gauss-Newton step. Since these errors are of 0(6m2) they may be sufficiently large such that the step 8m will not result in a decrease in the non-linear objective function. This problem is dealt with by taking a damped Gauss-Newton (DGN) approach, as was done in Section 4.2.1. By incorporating these pieces into the fixed /3 inversion algorithm it will be straightforward to design the GCV algorithm. 5.2 G C V and the 1-D E M inverse problem Haber (1997) successfully applied a GCV based inversion methodology to a non-linear problem in EM {viz magnetotellurics) using a DGN approach. This methodology will be applied to the 1-D frequency domain EM problem. It will first be applied to the case in which the parameter p is fixed during the inversion, and then the algorithm will be modified to accommodate the parameter p as part of the model. Each of these algorithms will be tested on synthetic AEM data. The final part of the section will deal with the application of GCV to small data sets. 5.2.1 G C V when p is fixed When investigating a new inversion methodology the best starting point is to solve the simplest problem possible. Therefore, I will begin by applying GCV to the inversion of EM data when p is fixed. Inversion Algorithm The inversion algorithm for GCV with a fixed value of p is very similar to the algorithm using a constant f3 for fixed p that was presented in Section 4.2.1. Each of the steps taken during the inversion are summarized in the flow chart shown in Figure 5.1. The differences between this algorithm and the constant f3 algorithm are in Steps 5 and 8, Chapter 5. Application of GCV to 1-D EM Inversion GCV Inversion Algorithm with a fixed value of p CO User Input (T)„ Calc initial values j : ' ( (3) Setk = 0 Initialize Variables, <8A; Use DGN 5m = o) 8m and m,k+1,= mlk,+ 8m 9A Setk=k+1 Update variables -•' \A•) Begin iter, k+1 i v — i CD Estimate p>(k+') using GCV J Calc GN step 8m ^ Setm(k+,)=m(k,+ 8m I *^CZ) Calc k+1 -values N •'"sT) Has (j) decreased? v c+i)+p(k+,)ci)<c+ri)c T . N / ® Convergence? 1 10) Exit Figure 5.1: Flowchart for the GCV inversion algorithm with a fixed value of p. Chapter 5. Application of GCV to 1-D EM Inversion 124 and the modifications are summarized in the following list: 5. The GCV estimate of f3 that corresponds to the minimum of the GCV function is selected as f3k+1. 8. The values of <p(m,k+i) and (j)(mk) are evaluated using /3 f c + 1, as in the discrepancy principle algorithm in Section 4.2.1. Now, having designed the algorithm, it can be tested on a synthetic example. Synthetic A E M Example A set of synthetic AEM data were generated over the conductivity model shown in Figure 3.1(a). The data were simulated at ten frequencies ranging from 110 Hz to 56320 Hz with a coil separation of 10 m and measurement height of 30 m. In order to simulate actual measurements, random noise with a standard deviation equal to 5% of the data amplitude + 10 PPM was added to each datum. The noisy data were then inverted using the GCV algorithm. The noise estimates used for the scaling matrix were 5% of the data amplitude + 10PPM and the reference model was 0.01 S/m. Figure 5.2 shows a comparison of the recovered models and predicted data from the discrepancy principle inversion and the GCV inversion. The final (3, model norm, data misfit, and target misfit values for both of the inversions are shown in Table 5.1. All of these results indicate that the two inversion results are both qualitatively and quantitatively similar. In comparison with the discrepancy principle however, the GCV algorithm has over-fit the data slightly. The convergence curves from both inversions are plotted in Figures 5.3 and 5.4. It is clear from the curves in these Figures that the paths taken by the two algorithms to reach the solution were quite different. While the data misfit is decreased at each iteration in both inversions, the decreases in the GCV case are not at a constant rate. The model norm does not begin small and increase gradually, as Chapter 5. Application of GCV to 1-D EM Inversion 125 (a) GCV (c) T r u e M o d e l R e c o v e r e d M o d e l (b) 50 100 D e p t h (m) T r u e M o d e l R e c o v e r e d M o d e l (d) 50 100 D e p t h (m) 1^5001 C CO SJ1000 « CL S 500 o o I P O b s . 0 — I P P r e d . x x Q O b s . / * Q P r e d . Jo JZ 10 10' 10' F r e q u e n c y ( H z ) o o I P O b s . — I P P r e d . x x Q O b s . / x Q P r e d . 10 10 F r e q u e n c y ( H z ) Figure 5.2: Comparison of recovered models and predicted data from AEM example discrepancy principle and GCV inversions with a fixed value of h. The measurement height was set equal to 30 m. Discrepancy principle (x2) results: (a) True conductivity (dashed Une) and recovered conductivity (sohd Une) and (b) observed data (inphase -circles, quadrature - crosses) and predicted data (inphase - sohd Une, quadrature - dashed Une). GCV results: (c) conductivity models (labeled as in (a)) and (d) data (labeled as in (b)) 4>m fa x2 GCV 1020.4 257.5 1.87 x IO"2 2.35 x IO"2 20.0 17.4 20.0 N/A Table 5.1: Comparison of parameter values from the AEM discrepancy principle and GCV inversions with a fixed value of h. fi is the final regularization parameter value, <£m is the recovered model norm, fa is the data misfit achieved by the inversion, and is the target misfit. Values are tabulated for both the discrepancy principle (x2) and GCV inversions. Chapter 5. Application of GCV to 1-D EM Inversion 126 10 o—o Recovered Misfit --* Target Misfit (b) 4 6 Iteration 10 (d) 2 4 6 Iteration 10' O 1 0 O CO n o O 10 4 6 8 10 Iteration » , • — G - G 0 new * - - » tj) old — ! 2 4 6 Iteration 10 Figure 5.3: Convergence curves from AEM example discrepancy principle inversion with a fixed value of h. The measurement height was set equal to 30 m. (a) fa (circles) and target fa (crosses) vs. iteration, (b) <pm vs. iteration, (c) fixed (3 vs. iteration, and (d) c/>(mfc+i) (circles) and <£(mfc) (crosses) vs. iteration. Chapter 5. Application of GCV to 1-D EM Inversion 127 Iteration Iteration Figure 5.4: Convergence curves from AEM example GCV inversion with a fixed value of h. The measurement height was set equal to 30 m. (a) fa (circles) and target fa (crosses) vs. iteration, (b) <f>m vs. iteration, (c) fixed /3 vs. iteration, and (d) <j>(mk+i) (circles) and 4>{m,k) (crosses) vs. iteration. Chapter 5. Application of GCV to 1-D EM Inversion 128 it does in the discrepancy principle example, instead it remains close to the same value throughout the inversion. The way in which /3 varies through the inversion is opposite in the two algorithm. When the discrepancy principle algorithm is used, (3 generally starts out large and decreases at each iteration. Conversely the GCV algorithm usually starts out by estimating small values of (3 which gradually increase as the inversion proceeds. The differences in the <j>m and ft curves in the GCV results are caused by the fact that the GCV tends to over-fit the noise associated with the linearized problem at each iteration. This can be seen by plotting the GCV function values as a function of /3 for a particular iteration along with the model that was calculated using the (3 chosen by the GCV. The plots were generated for iterations 1, 3, and 6 and are shown in Figure 5.5. In the first iteration the minimum in the GCV function lies at a small (3 value. The perturbation 8m calculated using this (3 was quite large and had to be damped once. While the damping ensured that the step decreased <j>, the model recovered at this iteration still has a lot of structure. By the third iteration however, the GCV curve has selected a larger j3 value and is calculating good perturbations. In fact, by the sixth and final iteration the minimum of the GCV curve has settled to a relatively constant value. As it stands, the GCV algorithm must work to remove the structure added to the model due through the selection of small (3 values during early iterations. While this is not the path the discrepancy principle algorithm takes, it leads to an acceptable result. The recovered conductivity structure is a good representation of the true conductivity and the data is fit to an appropriate level. Since the GCV algorithm has be successful for a fixed value of p there is no reason why it should not work when p is included in the model. Chapter 5. Application of GCV to 1-D EM Inversion 129 Figure 5.5: GCV Curves and recovered models from iterations 1, 3, and 6 of the AEM example GCV inversion with a fixed value of h. The measurement height was set equal to 30 m. Panels (a), (c), and (e) GCV function vs. /3 (solid line), f3 selected by GCV algorithm (circle). Panels (b), (d), and (f) model calculated using GCV estimate of y3 (solid fine), mode after it has been damped according to DGN approach (dashed fine), true model (dotted line), and ndamp is the number of dampings performed. Chapter 5. Application of GCV to 1-D EM Inversion 130 5.2.2 G C V when p is included in the model The next step is to modify the algorithm to include p in the model, such that it is possible recover both conductivity and a parameter value. Inversion Algorithm This algorithm is almost identical to the one presented Section 5.2.1 for the fixed p solution. The only modification is the inclusion of a cooling schedule for a p. These changes He in Steps 1,2, and 9A. A flow chart of the algorithm is shown in Figure 5.6. The changes in the algorithm are as follows: 1. The user must provide a starting value of a p . 2. The starting pref value is calculated as discussed in Section 4.2.2. 9A. The value of ap is decreased by a factor of 10 after each iteration. Now having designed the algorithm it can be tested on a synthetic example. Synthetic A E M Example The same synthetic AEM data used in the fixed p example was inverted with p as part of the model. For this inversion the starting value of ap was set equal to 1. Figure 5.7 shows a comparison of the recovered models and predicted data from the discrepancy principle inversion and the GCV inversion. The final measurement height, (3, model norm, data misfit, and target misfit values for both of the inversions are shown in Table 5.2. These results show that the GCV algorithm behaves similarly when p is included in the model. The discrepancy principle inversion recovered a measurement height of 29.1 m and the GCV recovered a value of 28.9 m, both of which are very good estimations of the true height of 30 m. In both cases the recovered conductivity models are good representations Chapter 5. Application of GCV to 1-D EM Inversion 131 GGV Inversion Algorithm with p included,in , • s, the model 8A Use DGN 8m = u) 8m and m = m + om 9 A Setk=k+1 Update variables 0Cp = 9 ocP {S) User Input including qP J 2 C a l c initial values including pref. J I (3) Set k = 0 - .1 I Initialize Variablesy ' . ^ •I C4j Begin iter, k+1 ; I ; .5) Estimate pk l )using GCV i I Calc GN step 8m j ^ w ^Setm(k+1)=m(k)+8m J "n v L ' ' Calc k+1 values j (%) Has <}) decreased? : N] r+rci,<c+r,)c N 9 Convergence ? 10 Exit Figure 5.6: Flowchart for the GCV inversion algorithm with p included in the model. Chapter 5. Application of GCV to 1-D EM Inversion 132 (a) GCV (c) T r u e M o d e l R e c o v e r e d M o d e l 50 100 D e p t h (m) T r u e M o d e l R e c o v e r e d M o d e l 50 100 D e p t h (m) (b) (d) 0. 0-52000 CQ 3 ° 1 5 0 0 C CC $ 1 0 0 0 cc sz CL £ 500 0 0 I P O b s . 0 — I P P r e d . x x Q O b s . / r - - Q P r e d . Of 10 10 10 F r e q u e n c y ( H z ) 3000i ^ 1 5 0 0 | C CO gjiooo CO JZ Q. S 500 0 0 I P O b s . 0 — I P P r e d . x x Q O b s . / r - - Q P r e d . /*/ 10 10 10 F r e q u e n c y ( H z ) Figure 5.7: Comparison of recovered models and predicted data from AEM example dis-crepancy principle and GCV inversions with h included in the model. The measurement height recovered with discrepancy principle was 29.1 m, with GCV was 28.9 m, and the true value was 30 m. Discrepancy principle (%2) results: (a) True conductivity (dashed line) and recovered conductivity (sohd line) and (b) observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - sohd line, quadrature - dashed line). GCV results: (c) conductivity models (labeled as in (a)) and (d) data (labeled as in (b)) h <f>m fa x2 29.1 1235.3 1.84 x 10~2 20.0 20.0 GCV 28.9 259.9 2.26 x IO"2 17.4 N/A Table 5.2: Comparison of results from AEM example discrepancy principle and GCV inversions with h included in the model, h is the recovered measurement height, /3 is the final regularization parameter value, <f}m is the recovered model norm, fa is the data misfit achieved by the inversion, and </>*, is the target misfit. Values are tabulated for both the discrepancy principle (%2) and GCV inversions. Chapter 5. Application of GCV to 1-D EM Inversion 133 of the true model. However, once again the GCV algorithm has over-fit the data slightly. The convergence curves from both of the inversions are plotted in Figures 5.8 and 5.9. Upon comparison, the convergence curves from these inversion are quite similar to those from the fixed p inversions. The plots of fa, fa, and fi versus iteration from the GCV inversion exhibit a character similar to that seen when the height was held fixed at 30 m. In the GCV inversion the p versus iteration plot follows a much more direct route to the final parameter value than that of the discrepancy principle solution. This is a function of the small values of fi values being selected by the algorithm. As discussed in Section 4.2.2, when fi is large the equations that govern Slog(cr) and Sp are decoupled and the steps taken by p. This is the reason why in the early iterations of the discrepancy principle inversions, p tends to wander slightly. However, the small f3 values selected by the GCV algorithm ensure that both parts of the step Sm are influenced by one another, and as a result p approaches its final value in a direct manner. Plots of the GCV function versus fi and the resultant models were generated for iterations 1,3, and 6 and are shown in Figure 5.10. The GCV curves and the recovered models from this inversion are quite similar to those from the fixed p inversion. The algorithm introduces unwanted structure in the first iteration by selecting a small value of fi, and the remaining iterations are used to remove the structure until a model which fits the data is found. For this example the GCV algorithm has worked well to recover good representations of both conductivity as well as a parameter value. 5.2.3 Application of GCV to small data sets The theory of GCV requires a large number of data. The previous examples have shown that GCV can work well with ten frequencies, that is twenty data. However, both of the field examples which have motivated this thesis consist of soundings with a small number of measurement frequencies. A sounding from the Sullivan data has four frequencies and Chapter 5. Application of GCV to 1-D EM Inversion 134 ioJl • • 1 no1' • • I 0 5 1 0 1 5 0 5 1 0 1 5 Iteration Iteration Figure 5.8: Convergence curves from AEM example discrepancy principle inversion with h included in the model, (a) fa (circles) and target fa (crosses) vs. iteration, (b) fa vs. iteration, (c) h (circles) vs. iteration and htrue (crosses), (d) fa vs. iteration, (e) /3 vs. iteration, and (f) (j)(mk+i) (circles) and <f>(mk) (crosses) vs. iteration. Chapter 5. Application of GCV to 1-D EM Inversion 135 10 o—o Recovered Misfit Target Misfit 2 4 Iteration (b) o-o Reed h *---* True h (d) 10" 10 o . Z 1 0 Q. 10 10 (e) 10 2 3 4 Iteration (f) O « 1 0 ' .Q o O 10 o-o (j) new h *-•* <j> old (1 '* 1 ' & ® o 3 4 Iteration Figure 5.9: Convergence curves from AEM example GCV inversion with h included in the model, (a) fa (circles) and target fa (crosses) vs. iteration, (b) fa vs. iteration, (c) h (circles) vs. iteration and htrue (crosses), (d) fa vs. iteration, (e) /3 vs. iteration, and (f) c/>(mfc+i) (circles) and c/>(m.fe) (crosses) vs. iteration. Chapter 5. Application of GCV to 1-D EM Inversion 136 i t = l A j r-l' d- i , (b) r J " " ' - 1 — GCV Model DGN ndamp = 1 50 100 Depth (m) 150 it=3 (d) — GCV Model — DGN ndamp = 0 50 100 Depth (m) 150 it=6 (f) GCV Model DGN ndamp = 0 50 100 Depth (m) 150 Figure 5.10: GCV Curves and recovered models from iterations 1, 3, and 6 of the AEM example GCV inversion with h included in the model. The measurement height was set equal to 30 m. Panels (a), (c), and (e) GCV function vs. fi (sohd line), fi selected by GCV algorithm (circle). Panels (b), (d), and (f) model calculated using GCV estimate of fi (solid line), mode after it has been damped according to DGN approach (dashed line), true model (dotted Une), and ndamp is the number of dampings performed. Chapter 5. Application of GCV to 1-D EM Inversion 137 P 4>m fa ti GCV 2231.7 1.44 x IO"3 3.10 x 10~2 2.61 8.0 12.5 8.0 N/A Table 5.3: Comparison of results from the four frequency AEM example discrepancy prin-ciple and GCV inversions with a fixed value of A. /3 is the final regularization parameter value, (j)m is the recovered model norm, fa is the data misfit achieved by the inversion, and 4>*d is the target misfit. Values are tabulated for both the discrepancy principle (%2) and GCV inversions. the Mt. Milligan data has three. In this Section I will apply the GCV algorithm for fixed p to a set of data collected using four frequencies. This will shed some light onto the applicability of my GCV algorithm to the field data sets. Four Frequency A E M Example In order to investigate how GCV behaves for a small number of data it is applied to a synthetic AEM data set. The data were generated over the conductivity model shown in Figure 3.1(a). The data were simulated at four frequencies ranging from 110 Hz to 56320 Hz with a coil separation of 10 m and measurement height of 30 m. In order to simulate actual measurements random noise with a standard deviation equal to 5% of the data amplitude + 10 PPM was added to each datum. In order to invert the data using GCV the data were scaled by the standard deviations of the errors which were equal to 5% of the data amplitude + 10PPM. The data were inverted with a reference model equal to 0.01 S/m. Figure 5.11 shows a comparison of the recovered models and predicted data from the discrepancy principle inversion and the GCV inversion of the four frequency data. The final /3, model norm, data misfit, and target misfit values for both of the inversions are shown in Table 5.3. These results show that the GCV algorithm has not been successful. The plot of the recovered model in Figure 5.11(c) does not resemble the true model, and the values from Table 5.3 show Chapter 5. Application of GCV to 1-D EM Inversion 138 Figure 5.11: Comparison of recovered models and predicted data from the four frequency AEM example discrepancy principle and GCV inversions with a fixed value of h. The measurement height was set equal to 30 m. Discrepancy principle (%2) results: (a) True conductivity (dashed line) and recovered conductivity (solid line) and (b) observed data (inphase - circles, quadrature - crosses) and predicted data (inphase - solid line, quadrature - dashed line). GCV results: (c) conductivity models (labeled as in (a)) and (d) data (labeled as in (b)) Chapter 5. Application of GCV to 1-D EM Inversion 139 10' (a) (c) 10 o-o Recovered Misfit *--* Target Misfit 10 20 Iteration 30 (b) (d) 10' 10 10 20 Iteration o-o § new * - - * ij> old * 10 20 30 Figure 5.12: Convergence curves from four frequency AEM example discrepancy principle inversion with a fixed value of h. The measurement height was set equal to 30 m. (a) (f)d (circles) and target fa (crosses) vs. iteration, (b) <£m vs. iteration, (c) fixed fi vs. iteration, and (d) 0(mfc+1) (circles) and <j>(mk) (crosses) vs. iteration. that the GCV solution in no way resembles the solution reached using the discrepancy principle algorithm. It is clear from this that my GCV algorithm is unstable when attempting to invert the four frequency data. However, determining the cause of the instability may help in developing a better GCV algorithm. In order to gain some insight into why the algorithm failed, it is useful to examine the convergence curves from the inversion. The plots of fa, <j>m, fl, and <f> versus iteration shown in Figure 5.12 are jagged. The most important curve, which is also the most erratic, is the regularization parameter. The GCV function was never able to reach a constant minimum value. A possible reason for this can be seen by Chapter 5. Application of GCV to 1-D EM Inversion 140 plotting the GCV curves from various stages during the inversion. GCV function versus fi and the resultant models, for the first three iterations, are plotted in Figure 5.10. The GCV curve from the first iteration has a well denned minimum and appears to be normal; however, the perturbation 8m which it generates is extremely large. The size of the step was damped once, such that the objective function decreased. The damped model shown in Figure 5.13(b) was physically unacceptable; however, the algorithm accepted the step and carried on to the next iteration. It is at this point that the problems continued. The GCV curve from the second iteration, which is plotted in Figure 5.13(c), does not have a unique minimum. As a result the algorithm chooses an even smaller (3 value than it chose in the first iteration. As would be expected, the step taken by the model in this iteration is another unphysical one. This type of behavior continued until the algorithm stopped after the thirtieth iteration. The behavior exhibited in this example could be construed as the failure of GCV however, this is not the case. When working with a linearized problem the solution that is attained can only be as good as the linearization itself. Therefore, when evaluating the GCV function during a given iteration, the chance of success depends in part on the model about which the problem is linearized. Therefore, it should not come as a surprise that the GCV func-tion calculated during the second iteration was problematic. If the model about which the problem was linearized was physically unreasonable the GCV had little chance of accurately estimating the noise. These results suggests that the instabilities of my algorithm could be ameliorated by adjusting the choice of fi during the first few iterations. While choosing the GCV estimate of fi in the first iteration may well predict the noise level correctly, it may also result in an unacceptable amount of structure being added to the model. A better choice of fi for the first iteration would be a large value which decreased the data misfit slightly without adding too much structure, such as the step taken by the discrepancy principle. This can Chapter 5. Application of GCV to 1-D EM Inversion 141 it=l — GCV Model — DGN ndamp = 1 50 100 Depth (m) 150 it=2 GCV Model DGN ndamp = 2 50 100 Depth (m) 150 it=3 GCV Model DGN ndamp = 0 50 100 Depth (m) (f) 150 Figure 5.13: GCV Curves and recovered models from iterations 1, 2, and 3 of the AEM example GCV inversion with a fixed value of h. The measurement height was set equal to 30 m. Panels (a), (c), and (e) GCV function vs. fl (solid line), fl selected by GCV algorithm (circle). Panels (b), (d), and (f) model calculated using GCV estimate of fl (solid line), mode after it has been damped according to DGN approach (dashed line), true model (dotted line), and ndamp is the number of dampings performed. Chapter 5. Application of GCV to 1-D EM Inversion 142 help to ensure that the recovered model will remain sensible and that the GCV estimate of fl will be reliable. An algorithm in which the value of (3 was gradually cooled with the eventual goal of using the GCV estimate of fl could help deal with the instabilities of the algorithm as it is currently denned. The implementation of such an algorithm is necessary in order to effectively assess the applicability of GCV to small data sets. 5.3 Summary In this section a GCV based inversion algorithm which makes use of the damped Gauss-Newton approach was developed. The algorithm was shown to work well when inverting frequency domain EM soundings which consist of twenty data. The GCV inversion algo-rithm was applied to both the problem of recovering conductivity as well as the problem of recovering conductivity and a survey parameter. The results from the synthetic AEM compared well with those attained using the the discrepancy principle algorithm devel-oped in Chapter 4. When applied to data sets with smaller numbers of data the GCV algorithm, as it has been designed, was unstable. This suggests the need for a new algorithm that combines the stability of the discrepancy principle and the noise estimat-ing properties of GCV. Overall, the use of GCV to estimate noise in non-linear inverse problems looks promising. Chapter 6 Summary This thesis seeks a method to improve the inversion results of the Mt. Milligan and Sullivan field data sets. It was conjectured that the problems associated with the data sets were due to geometric survey parameter errors. In the Mt. Milligan data, meas-urement height errors were causing distortions in the conductivity structure and in the Sullivan data, coil separation errors were making it impossible to predict the inphase data. A solution was proposed in the form of an algorithm that would recover both 1-D conductivity and a geometric survey parameter. Making use of background information on electromagnetic (EM) methods and the 1-D EM inverse problem, presented in Chapter 2, as well as the insight into the nature of geometric parameter errors, gained in Chapter 3, the stage was set to solve the problem. In Chapter 4 a modified inversion methodology was proposed to incorporate the extra parameter into the model recovered during the inversion. The necessary changes were made to the current algorithm and the new algorithm was designed. The algorithm was tested on synthetic AEM and HLEM data sets. The results from these tests showed that it was possible to recover both good estimates of parameter values as well as realistic representations of the subsurface conductivity structure. When applied to the field data sets the problems that were encountered before were corrected. The recovery of correct measurement heights in the Mt. Milligan case removed the distortions from the conduc-tivity models and in the Sullivan case the recovery of correct coil separations made it possible to predict the inphase data. 143 Chapter 6. Summary 144 The problem of noise estimation was also addressed in the thesis. Chapter 5 presented a discussion of the use of generalized cross validation (GCV) to estimate noise in non-linear inverse problems. A GCV based inversion algorithm was developed for the 1-D EM inverse problem. The algorithm worked well when applied to synthetic AEM data set with 20 data. However, some of the algorithms instabilities were revealed when it was applied to smaller data sets. While it appears that GCV worked well to estimate the noise associated with the data, the algorithm needs modification in order to fully exploit the potential of GCV. The result of these investigation was a better understanding of how GCV should be used within the context of non-linear inverse problems. References [1] Alumbaugh, D.L., and Newman, G.A., 1997, 3-D electromagnetic inversion for envi-ronmental site characterization: Proceedings of the Symposium on the Application of Geophysics to Environmental and Engineering Problems, 407-416. [2] deGroot-Hedlin, C., 1991, Removal of static shift in two dimensions by regularized inversion: Geophysics, 56, 2102-2106. [3] Dennis, J.E., and Schnabel, R.B., 1996, Numerical Methods for Unconstrained Op-timization and Nonlinear Equations: SIAM Philadelphia. [4] Ellis, R.G., and Shehktman, R., 1994, ABFOR1D k ABINV1D: Programs for for-ward and inverse modelling of airborne EM data, in JACI 1994 Annual Report: Dept. of Geophysics & Astronomy, University of British Columbia. [5] Fitterman, D.V., 1998, Sources of Calibration Errors in Helicopter EM Data: Proc. of the International Conference on Airborne Electromagnetics, 1.8.1-1.8.10. [6] Fraser, D . C , 1986, Dighem resistivity techniques in airborne electromagnetic map-ping: in Palacky, G.J., Ed., Airborne Resistivity Mapping: Geol. Survey of Canada, Paper 86-22, 49-54. [7] Fullagar, P.K., and Oldenburg, D.W., 1984, Inversion of horizontal loop electromag-netic frequency soundings: Geophysics, 49, 150-164. [8] Haber, E., 1997, Numerical Strategies For the Solution of Inverse Problems: Ph. D. Thesis, Univ. of British Columbia. 145 References 146 [9] Jones, F.H.M., 1996, Preliminary Report: Geophysical Surveys at Sullivan Mine Tailings: Dept. of Earth & Ocean Sciences, Univ. of British Columbia. [10] Kovacs, A., Holladay, J.S., and Bergeron, C.J.Jr., 1995, The footprint/altitude ratio for helicopter electromagnetic soundings of sea-ice thickness: Comparison of theo-retical and field estimates: Geophysics, 60, 374-380. [11] Norton, S.J., San Filipo, W.A., and Won, I.J., 1999, Reciprocity formulas for elec-tromagnetic induction systems: Proc. of the Symposium on the Application of Geo-physics to Environmental and Engineering Problems, 961-970. [12] Oldenburg, D.W., Li, Y., and Ellis, R.G., 1997, Inversion of geophysical data over a copper gold porphyry deposit: A case history for Mt. Milligan: Geophysics, 62, 1419-1431. [13] Parker, R.L., 1994, Geophysical Inverse Theory: Princeton University Press. [14] Pavlis, G.L., and Booker, J.R., 1980, The mixed discrete-continuous inverse prob-lem: application to the simultaneous determination of earthquake hypocenters and velocity structures: Journal of Geophysical Research, 85, 4801-4810. [15] Ryu, J., Morrison, H.F., and Ward, S.H., 1970, Electromagnetic fields about a loop source of current: Geophysics, 35, 862-896. [16] Sauck, W.A., Atekwana, E.A., and Bermejo, J.L., 1998, Characterization of a newly discovered LNAPL plume at Wurtsmith AFB, Oscoda, Michigan: Proc. of the Sym-posium on the Application of Geophysics to Environmental and Engineering Prob-lems, 399-407. [17] Wahba, G., 1990, Spline models for observational data: SIAM Philadelphia. References 147 [18] Weaver, J.T., 1994, Mathematical methods for geo-electromagnetic induction: John Wiley & Sons Inc. [19] Ward, S.H., and Hohmann, G.W., 1988, Electromagnetic theory for geophysical applications: in Nabighian, M.N., Ed., Electromagnetic Methods in Applied Geo-physics: Soc. Expi. Geophys., Vol 1, 131-311. [20] Wynn, J., and Gettings, M., 1998, An airborne electromagnetic (EM) survey used to map the upper San Pedro River aquifer, near Fort Huachuca in southeastern Arizona: Proc. of the Symposium on the Application of Geophysics to Environmental and Engineering Problems, 993-998. [21] Zhang, Z., and Oldenburg, D.W., 1994, Inversion of airborne data over a 1-D earth, in JACI 1994 Annual Report: Dept. of Geophysics & Astronomy, University of British Columbia. [22] Zhang, Z., and Oldenburg, D.W., 1999, Simultaneous reconstruction of 1-D suscep-tibility and conductivity from electromagnetic data: Geophysics, 64, 33-47. [23] Zhang, Z., Routh, P.S., Oldenburg, D.W., Alumbaugh, D.L., and Newman, G.A., submitted for publication, Reconstruction of ID conductivity from dual-loop EM data: Geophysics. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0053001/manifest

Comment

Related Items