INVERSION OF E M DATA T O R E C O V E R 1-D C O N D U C T I V I T Y A N D A G E O M E T R I C SURVEY P A R A M E T E R By Sean Eugene Walker B. Sc. (Honours), Geology & Physics, McMaster University, 1996 A T H E S I S 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 REQUIREMENTS FOR T H E DEGREE OF M A S T E R OF SCIENCE 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 SCIENCES (GEOPHYSICS) 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 degree at the thesis in University of partial fulfilment of of department this thesis for or by his or scholarly purposes may be her representatives. permission. Department of The University of British Columbia Vancouver, Canada DE-6 (2/88) 18, 1779 for an advanced Library shall make it agree that permission for extensive It publication of this thesis for financial gain shall not Date requirements British Columbia, I agree that the freely available for reference and study. I further copying the is granted by the understood that head of copying my or be allowed without my written 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 E M (AEM) and ground based horizontal loop E M (HLEM). When dealing with A E M 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 2 Introduction 1 1.1 Motivation for the Thesis 2 1.2 Outline of the Thesis 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 3 10 Details of the Inverse Problem 24 2.2.1 The Forward Problem 24 2.2.2 The Inverse Problem 29 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 3.2 3.3 How do measurement height errors affect the inversion results? . . 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 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 4.1 4.2 4.3 4.4 44 62 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 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 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 Field Examples 108 4.4.1 Mt. Milligan 110 4.4.2 Sullivan 114 v 4.5 5 114 Application of G C V 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 5.3 6 Summary 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 Summary 142 Summary 143 References 145 vi List of Tables 3.1 Results from the inversion of AEM data with the incorrect h 3.2 Results from the inversion of small separation HLEM data with the incor- eat values. . . rect s value 3.3 54 Results from the inversion of small separation HLEM data with the incorrect s value and noise estimated to account for inphase modelling errors. 3.4 44 56 Results from the inversion of large separation HLEM data with the incorrect 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 5.1 107 Comparison of parameter values from the AEM discrepancy principle and GCV inversions with a fixed value of h 5.2 Comparison of results from AEM example discrepancy principle and GCV inversions with h included in the model 5.3 125 132 Comparison of results from the four frequency AEM example discrepancy principle and GCV inversions with a fixed value of h vii 137 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. 1.4 Comparison of the inphase and quadrature data from two soundings along . . 5 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 looploop E M 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 A E M 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 A E M soundings with h 43 true Vlll = 24 m, 30 m, and 36 m. 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 3.7 Synthetic data from small separation HLEM soundings with s 48 eat = 10 m and Strue = 9 m, 10 m, and 11 m 3.8 50 Synthetic data from large separation HLEM soundings with s e3t = 50 m and s^ue = 49 m, 50 m, and 51 m 3.9 52 Results from the inversion of small separation HLEM data with the incorrect 5 value 55 3.10 Results from the inversion of small separation HLEM data with the incorrect 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 incorrect 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 fixedfiinversion with afixedvalue of h 4.3 74 Convergence curves from AEM example fixed j3 inversion with afixedvalue of h 4.4 75 Flowchart for the discrepancy principle inversion algorithm with a fixed value of p 4.5 78 Recovered models and predicted data from AEM example discrepancy principle inversion with a fixed value of h 4.6 80 Convergence curves from AEM example discrepancy principle inversion with afixedvalue 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 4.8 Results from eleven inversions of the small separation HLEM sounding data, 4.9 83 Stme = 11 na, with fixed values of s ranging from 12.5 m to 7.5 m. . 84 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 a = oo and 0 88 p 4.11 Values calculated during the line search for (3 during the first iteration of the inversion of the AEM data with a = oo, 10 , and 10 -6 -3 p 4.12 Results from the search for a best fit p f re value for the AEM example. . . 4.13 Results from the search for a best fit p f value for the HLEM examples. Te 90 92 93 4.14 Contour plots of (f>d for a range of cr f and p values for the AEM example Te and both HLEM examples 95 4.15 Flowchart for thefixed{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 example 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 example 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 Ill 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. Milligan 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 5.3 Convergence curves from AEM example discrepancy principle inversion with a fixed value of h 5.4 126 Convergence curves from AEM example GCV inversion with a fixed value 127 of h 5.5 125 GCV Curves and recovered models from iterations 1, 3, and 6 of the AEM example GCV inversion with a fixed value of h xi 129 5.6 Flowchart for the GCV inversion algorithm with p included in the model. 5.7 Comparison of recovered models and predicted data from A E M example discrepancy principle and GCV inversions with h included in the model. . 5.8 132 Convergence curves from AEM example discrepancy principle inversion with h included in the model 5.9 131 134 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 frequency A E M example discrepancy principle and GCV inversions with a fixed value of h 138 5.12 Convergence curves from four frequency AEM example discrepancy principle 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 xii 141 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 enjoyable. 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 estimation 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 E M 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 configuration 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. 3 Introduction o -40 • o 7040 Hz ° 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, 4 Introduction 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. (a) Co 0 N X Introduction 5 20 o s 1^ - o 2 0 f ° Voo° p c 9=-40. 0 0 oOOOCe 50 o °°o 00 100 00 oPWP^W^co^eP^oo 150 200 250 « 'o°ooo 300 c IP Obs IP Pre -- Q Obs QPre 350 Easting as 20 XL --20f°°^o o .oOo^oP°coooo° o ^ o c o ^ c o o a ^-40 50 100 150 200 0 0 0 250 o„cP 300 350 Easting 20 (c) N X oo pod O £-40 50 100 150 200 250 300 350 250 300 350 Easting £ 20 (d) £ 0 8-20?° V x ^ V * o Q.-4Q 50 0 0 100 150 200 Easting Figure 1.3: Observed and predicted data from the inversion of the Sullivan data. 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 1. 6 Introduction 20 20 G O IP 10 • X X Q OOIP 10 X X Q •o X x X 00 X I d CO X X X o "D p !-io "D 0) (0 !-io in CB |-20 "I-20 -30 10 O o 10* Frequency (Hz) 10' -3010 o O 10' Frequency (Hz) 10 a 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 o o o- — - o -B a 250 200 900 Hz IP x .1 900HzQx.1 7200 Hz IP x.1+25 7200 Hz Q x .1+25 56 kHz IP x .1+75 56 kHz Q x . 1+75 -0 0-—0 C/5 co150 cd Q 100 50 4 o - G . _ _ ^ . o -o - © - e - - e - o- -a - o - - e - e - e — e — e — e — o o o — o o 'o—e—e—e—o o—e—e—e—o' o n n n =& 9D00 e ) 1 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 8 Chapter 1. Introduction 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 similar 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. 1.2 Introduction 10 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 E M 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 E M 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 looploop E M 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 looploop surveys. The four most common geometries are horizontal coplanar (HC), vertical coplanar (VC), coaxial (CA), and perpendicular (PP). Examples of these transmitterreceiver 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 14 Chapter 2. Background Information S M ™ "i| Rx Tx y Z ; X —' h • \ 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 I is the amplitude of the current, o> is the angular frequency, and t is time. This 0 current produces the primary magnetic field H . p 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 V , and it has the form 0 V, = - (2.2) w where $ is the magnetic flux through the coil which is defined as $ = NJ B-dS, R (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 = -iufioN J H p R • dS, (2.4) where / x is the magnetic permeability of free space. Since the radius of the receiver coil 0 is assumed to be much smaller than the separation of the transmitter and receiver coil I can assume that H p is constant across the surface S. This simplifies the expression for the primary voltage to Vo = -UJ^N A H (S), (2.5) p R R where AR is the area enclosed by the receiver coil and H P 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 transmitterreceiver pair, shown in Figure 2.3(b). When H P 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 H . S a combination of both H and H . P S The field sensed at the receiver coil H R will be 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 (H (s) + fff (u;, cr, h, s)) . P (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 E M 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 18 Chapter 2. Background Information This is the value that is recorded by the EM system. However, it is important to note that H S lags H , P due to the inductive interaction of H P 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~\ H z \ S ( ) x 211 ) and Quadrature = Im (H*(u;,(T,h,s)) ' ' ' x f. V (2.12) v H z { S ) where £ is a multiplicative factor. This factor is usually set equal to 100 or 10 such that 6 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 m T = IA N , T T (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 = iwp N A -^-. 0 R R Airs 6 (2.15) 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 A E M 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 isflown,it is determined using a laser altimeter mounted Chapter 2. Background Information 20 S SB Tx Rx Bx VB V gain control 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 A E M systems as described in Fitterman (1998). The bucking coil, shown in Figure 2.5, is located between the transmitter and the receiver. Since the goal is to remove the primary field from the field sensed at the receiver coil, I want the voltage VB in the bucking coil to be entirely due to the primary 21 Chapter 2. Background Information 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 A E M 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 A E M 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 T ~ = -^r ' Z 1 Zo - V - aVn1 (2 16) OLVB and since I am assuming that Z OLVB = Vo I get Hf((jj,a,h,s) / n * - = Ww • (217) 1 Due to the fact that the secondary field values are very small in A E M 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 ) 22 Chapter 2. Background Information s Figure 2.6: Cartoon of a typical HLEM system, h is the measurement height and s is the coil separation. and Quadrature(PPM) = Im (Hf(u,(T,h,s)) ^ , } H „ x 10 . 6 (2.19) 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 amplitude 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 parametric 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 23 Chapter 2. Background Information 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 H z \ S ' ) } x 100, (2.22) ) anc / < w x Quadrature(%) = Im(H?(w,a,h,s)) ' ' x 100. V H z K S (2.23) ) Comparison of A E M and H L E M While both A E M and HLEM surveys can be carried out with horizontal coplanar coils there is an important difference between them which should be reinforced. As was mentioned 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 practice, accurately positioning the coils at a fixed distance s can be difficult. Therefore, both A E M 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 important 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 i th Re(di) = measurement di will have the form Re ( f f f K c r , H (si) P (2.24) Chapter 2. Background Information 25 and Im(di) lm (iff (u>,,<7, hi, Si)) (2.25) H (si) p where each of the subscripted parameters represent the parameter value associated with the i th 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 instead 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 represented 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 inhomogeneity I am dealing with. A E M 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 structures. 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 necessary. 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 26 Chapter 2. Background Information Rx surface Gi hi di hi CJM flM i = 2,... M - l 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 i layer, and <7; is the conductivity of the :th layer. th 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 0 th layer and the thickness of this layer is equal to the measurement height such that h = h. The 0 conductivity of the air is <r = 0. 0 As with any E M problem I start with Maxwell's equations. I will make use of the 27 Chapter 2. Background Information 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 E M 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 , (2.27) s 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. H = - j j - * oz 1d u>n H = — — {rEj,), r or 0 0 ^ oz - ^ oz (2.28) T (2.29) z = „J5, + J.. (2.30) 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 H from Equation 2.28, and H from r z 28 Chapter 2. Background Information Equation 2.29 into Equation 2.30, I get the following partial differential equation for the tangential electric field (F \ 9 2 where A; = 2 1d 9 dr (rE*) + k Ej, = (2.31) iup Js, 0 Following Ryu et al. (1970) I can transform Equation 2.31 into the —IUJLIQCT. Hankel domain to get an ordinary differential equation for the tangential component of the electric field 0 s where u = A — k , and 2 (2.32) E,p = iu)p, j , u 2 and J, are the Hankel transforms of Ej, and J . 2 s As mentioned previously, these equations only hold for regions of constant conductivity. 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 -iu)Lio 47T e -2uo'»o -R Ji{\r)\ d\, 2 TE (2.33) 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 Z l R t E ~ -Z (2.34) 0 Zi + Zo' where the intrinsic impedance of the i th Zi interface Z{ is defined as 29 Chapter 2. Background Information and the input impedance of the i layer Z is denned by the recursion relation th l -Uj2/ij I Z i+1 + Z (1 + e" r Z = Zi l -U{2hi i = M,...,l, (2.35) -Ui2hi Zi + z^ (1 + e--Uj2hi where 17M+I ry (2.36) — Aftf+1- 6 However, the value I want to calculate is H (r,u>, z). Using Equation 2.29, and applying z the inverse Hankel transform, I get H (r,u;,h ) z 0 TUT t°° = / 4TT J 0 e~ ° 2u ho R J {Xr)X d\. (2.37) 3 TE U 0 0 The expression for H (r,w, ho) in Equation 2.37 is equal to the secondary magnetic field z 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 cr (z). The observed data d° " is a vector containing h true 7Y data which are the result of an EM experiment overtop of o~ (z). The data include true some unknown amount of measurement noise. The goal of the inversion is to recover the function o~ (z) however, this problem is seriously underdetermined since I have only true 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 <p of the form m d(M-M y w (z) dz 2 <Pm{M) = a, J w,(z)(M - M ffdz re +a J z ref dz, z (2.39) where M. f is some reference model. The first term of (p is a measure of how close M. re m is to the reference model M. f- This term is referred to as the smallest model norm. The Te 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 a determine the z relative importance of the smallest andflattestcomponents of the model norm and the functions w,(z) and w (z) are spatial weighting functions which can be used to include z further information about the model. This formulation provides the ability to recover the smallest model when a, = 1 and a = 0, theflattestmodel a, = 0 and a — 1, or any z z combination of the two. Usually I want to find a model that is as featureless as possible and hence a, and a are selected such that the second term in <f) is dominant over the z m Chapter 2. Background Information 31 first. When no extra information is available it is common to set w (z) = w (z) = 1. s z This is the general definition of <f) that I will adopt throughout the thesis. m 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{<7 ),...,log{a )] . (2.40) T 1 M It is important to note that M » N and therefore, the problem remains underdeter- mined. The discrete form of <£ in Equation 2.39 will be m cj> {m) = {m- m f m [a WjW Tef a +aW W] T a z z (m - m ), (2.41) ref where W is the smallest model weighting matrix and W is the flattest model weighting 3 z matrix. These terms can be combined such that 4> {m) = {m- m ) WlW (m T m where W m ref m - m ), (2.42) Tef 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™) = \\W (m - m )\\ , 2 m ref (2.43) where || • || is the Euclidean 2-norm. The design of the model objective function is such that if I simply minimize <p , and m W m 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° ° = F[m ] + e, b true where m (2.44) is the discrete representation of the true model, F is the forward operator t r u e 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> is a measure of how well the data predicted from a given model m fits the d 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 i th datum. I can rewrite 4>d in matrix notation as follows fa=\\W (F[m}-d° °)\\ , b (2.46) 2 d where Wd is the diagonal data weighting matrix which is equal to W = diag(-,...,—). (2.47) d \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 offafrom Equation 2.46 I see that it is a random variable with a % distribution 2 and therefore, it will have an expected value approximately equal to N. Thus the target misfit <j> should also be approximately equal to N. d 33 Chapter 2. Background Information 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> such that fa = <p . m (2.48) d 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 ( F H - ^ ) | | d 2 + ^||W (m-m m )|| , (2.49) 2 r e / where (3 is the regularization or trade-off parameter. The problem will then become: Minimize <f> = fa + {3<j> such that fa = <f> , ^ m (2.50) d 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™) = \\W (Gm-d° ')\\ b d 2 + (3\\W (m-m )\\ . (2.52) 2 m ref 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) = 2G WjW {Gm - (f**) + 2f3WZW {m - m ). T d m (2.53) ref Setting g(m) = 0 leads to the matrix system of equations [G WjW G T d + L3WlW )m = G WjWdd° T m b> + f3WZW m , m ref (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) is the target misfit, <f> is the model objective function, and f3 is the trade off parameter. d m 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) ((3). This plot is called the Tikhonov curve. A cartoon of a typical Tikhonov m curve is shown in Figure 2.8. This curve illustrates how the choice of /3 dictates the values of $d and (j) and therefore, the character of the recovered model. When /3 gets m large, 4>(m) ~ /3(f> and the minimization is strictly searching for the minimum structure m 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> and m 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 offindingj3 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> . If such a /3 does exist, as shown in Figure 2.8, it is a simple d 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 particular /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) = 2J WjW (F[m] - d° ') + 2/3W^W (m - m T b d 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 m . The linearized version of g is k expressed as g(m + 8m) = g{m ) + H(m )6m, k k where H(m ) k is the Hessian matrix, evaluated at m , which is defined as k H(m ) = k (2.57) k Jp^ J _,*) + ... 9 = 2 om W om k Wd{F[mk] k 2J{m ) WjW J(m ) T k d k + 2/3W£W . m (2.58) In order to avoid evaluating the derivative of the Jacobian, the Hessian is approximated as H(m ) k * 2J(m ) WjW J{m ) T k d k + 2/3W^W . m (2.59) Chapter 2. Background Information 36 + 8m) — 0 yields the matrix system of equations Setting g(m k [J(m ) WjW J(m ) + L3WlW )8m T k d k m J(m ) WjW (d° > T k b d = ... - F[m }) k fiW^W^m - m ), (2.60) ref which can be solved to give 8m such that the model for the next iteration will be equal to m i k+ = m k + 8m. T h e process of linearization and iteration is continued until the solution converges to the m i n i m u m of Equation 2.49. This process can be repeated for a number of j3 values, as i n the linear case, to find the f3 which satisfies (f> — cj> . d 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 incomprehensible. For this reason the pertinent survey parameter values are recorded with the data output from the E M 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 interpretation 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 geometric 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 E M system records in the data output is referred to as the expected or 37 38 Chapter 3. Geometric Survey Parameter Errors estimated value of p and it is represented with the variable p . eat The true value of p is represented with the variable p e- The parameter error is represented with the variable trU 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. Thefirstsurvey simulates a typical A E M 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 A E M survey consists of measurements of the inphase and quadrature components 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 A E M sounding. The conductivity model over which the A E M 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 A E M sounding it is necessary to define the true measurement height value. Both of the HLEM surveys consist of measurements of the inphase and quadrature 39 Chapter 3. Geometric Survey Parameter Errors 2500i -0.5i (a) AEM O O Inphase x x Quad. E w w (b) S2000| Q. pO x §1500| > o 5-1.5 c o O TJ 9 co g- 500 o -2.5 O x C01000 50 100 Depth (m) 150 10 a- X O x O 10 10 10' Frequency (Hz) J -0.5i O O Inphase (c) g, (d) -1 x x Quad. Small HLEM TJ -1.5 c o O O -2.5 50 100 Depth (m) 6 6 o ° oi^ 10 150 3 10 10" Frequency (Hz) 30 -0.5i o 20 (e) Large HLEM E —» cn (f) g I3 o O o 10 o5 6 o •o-10| CO -1.5 o O 8-20I CO cL-30 -2.5 50 100 Depth (m) 150 O O Inphase -40 -50' 10' x x Quad. 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 h true 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 A E M system. Since this value may not be equal to h tTue I will refer 41 Chapter 3. Geometric Survey Parameter Errors 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 h . est Using Equation 3.1 the measurement height error 8h can be expressed as Sh = h true A cartoon of the A E M system and the values used to calculate h (3.2) eat 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 Equation 3.2 I get the expression Sh = 6a — 6b, (3.4) 42 Chapter 3. Geometric Survey Parameter Errors Case 2. Case 1. 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 h true denote the true height, the AEM data from Equations 2.18 and 2.19 are Inphase(PPM) = Re H^(UJ,(T, h , S)" true (3.5) x 10 e and Quadrature(PPM) = Im [ ^i"^*™,*) H {s) } ^ (3.6) 1 Q 0 p By generating synthetic data from A E M soundings with different h true 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 A E M sounding data 43 Chapter 3. Geometric Survey Parameter Errors 3500 4000, (a) < 00 10 10 Frequency (Hz) Frequency (Hz) Figure 3.4: Synthetic data from AEM soundings with h = 24 m, 30 m, and 36 m. (a) The inphase component of the data ( h — 24 m - diamonds, h = 30 m - circles, and htme — 36 m - squares) (b) The quadrature component of the data ( h = 24 m diamonds, h = 30 m - crosses, and h e — 36 m - squares). t r u e t r u e t r u e t r u e t r u e t r U 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 htrue t r U e = 24 m data are represented by diamonds and the = 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. 44 Chapter 3. Geometric Survey Parameter Errors h »t 4>m fa 24 m 0.136 20.0 30 m 0.127 20.0 36 m 0.353 25.1 e 20.0 20.0 20.0 Table 3.1: Results from the inversion of AEM data with h values of 24 m, 30 m, and 36 m when h e is equal to 30 m. <f) is the recovered model norm, fa is the data misfit achieved by the inversion and <f> is the target misfit. These values correspond to results plotted in Figure 3.5 eat trU m d 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 h eat 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 Figures 3.5(a), 3.5(c), and 3.5(e) and the observed and predicted data are plotted in Figures 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 % misfit level of 20 for both h 2 eat fairly well (fa « 25) for h eat — 24 and 30 m and has done = 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 h eat is different. Since the inversion results from h eat = 30 m (8h = 0) are free of measurement height errors I will adopt them as the results against which I compare the other inversions. 45 Chapter 3. Geometric Survey Parameter Errors 3000. - True Model — Recovered Model (a) (b) 55 24 m jr-J 50 Depth (m) —-0.5 (c) h = 30 m 150 — True Model — Recovered Model I (d) I > 1-1.5| •o c o O -21 -2-2.5] 50 100 Depth (m) 150 10' 10' Frequency (Hz) 3000, -•- True Model — Recovered Model (e) (f) h = 36 m -31 50 100 Depth (m) 150 10* 10" Frequency (Hz) Figure 3.5: Results from the inversion of A E M data with h values of 24 m, 30 m, and 36 m when h 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. tTUe 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 h est = 24 m is slightly pushed down in comparison to the structure recovered with h eBt = 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 variations 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 s eat 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 s etru Using Equation 3.1, I can express the coil separation error Ss as Ss — s true - s . eat (3.7) When a survey is performed on level ground separation errors are likely due to incorrect 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. 48 Chapter 3. Geometric Survey Parameter Errors Case 2. \Ss > 0 Case 1. \6s < 0 "si.V 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 s tTUe and s eat 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 ( i f f (Strue) + H^(strue)) V = -iup m, 0 (3.8) However, the estimate of the primary voltage calculated by the HLEM system depends upon s . Therefore, the voltage Vc is denned as eat V = c (3.9) -iu>p m H (s t). P 0 a e!> 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 H JU, *, h, Strue) + H {s ) Z_ _ Z 0 S ~ - P Z true H?(s ) eat H?{s ) (3.10) eat From this expression it is clear that when s eat ^ 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 s eat Inphase(%) = Re f H!(^ h^ s^e)\ as H (strue,seat), P x l f j 0 + S (3.11) 49 Chapter 3. Geometric Survey Parameter Errors and (3.12) x 100 Quadrature(%) = lm where 8H is the primary field error and is denned as P (3.13) These equations show that both s e trU and s est 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 8H in the inphase component. P In order to get more insight into the behavior of 8H I substitute the definition of P the primary field from Equation 2.13 into Equation 3.13 and end up with (3.14) \\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 8H term. From Equation 3.14 I see P that when 8s > 0 (s > s ) the value of 8H will be negative and therefore the inphase P true est data will be shifted in the negative direction. Similarly when 8s < 0 (a true < s , ) the e t value of 8H will be positive and the inphase data will be shifted in the positive direction. P 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 s t value for all of ea the small separation soundings is set to 10 m and the s true value for each sounding is 11 m, 10 m, and 9 m respectively. These values correspond to coil separation errors of Chapter 3. Geometric Survey Parameter Errors 50 (a) 40 >--o--o--^-o-o-o-" "' 30 0-0 <> s 50 ,-4 <3> (b) - 9m £ 20 G - O s f ^ l O m B-B s. = 11m true cu <2 C O 10 0--Q ob>--G-e-0-e--o-©- " e -10 -20 -B' - a L.. .. ..^... -e-a-B-' -30' 10' 10* 10* 10" Frequency (Hz) Q Q & 10 10 Frequency (Hz) i>~o~o~^~o--o~o~^~o--o 00 (c) s, = 9m O-O s = 10m „ „ ., ., Q-Q s. = 11 m true r U 8 (d) TRUE Q. X ra s oE>-o--e-e-e--e-e-e-e--o Q &--e--Q—a—B--Q--B—B~-B--0 -30' 10' 10 10 10 Frequency (Hz) 10 10 Frequency (Hz) Figure 3.7: Synthetic data from small separation HLEM soundings with s = 10 m and •Strue = 9 m, 10 m, and 11 m. (a) Inphase component of the data (s = 9 m - diamonds, Stme = 10 m - circles, s = 11 m - squares), (b) quadrature component of the data (•Strue = 9 m - diamonds, s = 10 m - crosses, s = 11 m - squares), (c) real part of normalized secondary field (as in panel (a)), and (d) primaryfielderror (as in panel (a)). est true true true true + 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 51 Chapter 3. Geometric Survey Parameter Errors the data respectively. The shift due to 8H seen in the inphase component, is much P 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 A E M 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 s eat = 10 m is equivalent to the shift caused by 8s = 5 m in a survey with s at — 50 m. However, it is unlikely to encounter coil separation errors larger than a e 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 s eat 50 m and the value of s true is equal to 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 s true = 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 35 30 (a) 25 (b) „ 20 5? ~Z o) 1 ^ / / / 5 a-. „ © \ E3 \ V j= 10| CL 51 £ o^-0-^ y , „ . J 3 -5b-Q 0-0 s G>-Os Q-Q . s -10' 10 t r u e , r u e true = 49m'^ =50m = 51m 10 10 Frequency (Hz) 10" 10 10 Frequency (Hz) 30 0-0 s = 49m G-O s / = 50m ., —.25 _Q--Q„ s.t r u e = _51m ( r u e u e (c) true /v "0 {>-o-0"^-^>~o-o--^--o--o 6 ::<3\ 0-0 s = 49m Q-O s , = 50m „ „ true Q-Q s. = 51m (d) 120 r u e true CL X. ra BCD o&-o-o-e-e-0-e-e-e-G ilO| CO 6 cr Q vO • Jf 10" 10" 10' Frequency (Hz) 10" -6' 10' 10° 10" Frequency (Hz) 10" Figure 3.8: Synthetic data from large separation HLEM soundings with s = 50 m and Stme = 49 m, 50 m, and 51 m. (a) Inphase component of the data (strue = 49 m diamonds, s ue = 50 m - circles, s = 51 m - squares), (b) quadrature component of the data (s e = 49 m - diamonds, s e — 50 m - crosses, s ue = 51 m - squares), (c) real part of the normalized secondary field (as in panel (a)), and (d) primary field error (as in panel (a)). e a t tr trU tPue trU tr 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 algorithms 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 54 Chapter 3. Geometric Survey Parameter Errors Strue 4>m 11 m 3.46 xlO" 10 m 0.116 0.377 9m 3 <f>d 2.20 xlO 20.0 3.58 xlO 20.0 20.0 20.0 4 4 Table 3.2: Results from the inversion of small separation HLEM data with an s value of 10 m when s is equal to 11 m, 10 m, and 9 m. <f> is the recovered model norm, 4>d is the data misfit achieved by the inversion and <f> is the target misfit. These values correspond to results plotted in Figure 3.9 true m d 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 s tTUe s est values of 11 m, 10 m, and 9 m, and = 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 tofitthe inphase data. This suggests that the Chapter 3. Geometric Survey Parameter Errors 55 —- True Model — Recovered Model -0.5 (a) (b) g 11 m. o o IP Obs. — IPPred. x x QObs. 1-10] — QPred. o ra 0] w •15 o-20 50 100 Depth (m) -25^ 150 10 o n 10' 10' Frequency (Hz) 10 7 True Model Recovered Model -0.5 (c) Sir (d) = 10 m. g o o IP Obs. 6 — IPPred. x x QObs. QPred. 5 § 4 o Y (0 31 S 2l (0 £.1 -2.5 50 100 Depth (m) 1r —- True Model — Recovered Model 0.5 (e) f Pi, o Str = 9 m. Io' 150 10 10 Frequency (Hz) 10 50i ^ b o o o o °o 40 (f) g J30 o o ip Obs. — w IPPred. :&-0.5T O 20 x x QObs. 1 o " « 10| 1 <1) O) (0 Q. O) o -2 -4 QPred. TJ "-10I 50 100 Depth (m) 150 -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 s 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. true Chapter 3. Geometric Survey Parameter Errors ^ true 11 m 10 m 9m K 0.086 0.116 0.061 56 (j>d 20.0 20.0 20.0 20.0 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 is equal to 11 m, 10 m, and 9 m and noise estimated to account for inphase modelling errors. (p is the recovered model norm, (pj is the data misfit achieved by the inversion and (p* is the target misfit. These values correspond to results plotted in Figure 3.10 true m d 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 true = 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 Figures 3.10(b) and 3.10(f) show that it is not possible to adequately reproduce the inphase observations for the s ue tr = 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 True Model Recovered Model -0.5 (b) g •Strue = 11 HI. o o IP Obs. — IPPred. x x QObs. -—• QPred. o I-10I cs CD cc3-15| Q. C - 50 100 Depth (m) (c) 10 (d) m. o -25^ 10 150 True Model Recovered Model -0.5 Strue = -20] 10 o 10 Frequency (Hz) o o |p Obs. — IPPred. x x QObs. - - Q Pred. g o 1 3| co CD ol w 'I IS f 1 -2.5 50 100 Depth (m) 150 10 10 Frequency (Hz) 10 45 True Model Recovered Model -0.5 40 (e) •Sfr -1 10 S35 ) ° o O c o o o o IP Obs. — IPPred. x x QObs. 0 5| - - Q Pred. TJ 9 m. 2 c320| CD S3is| -2.5 50 100 Depth (m) 150 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 s e 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. trU 58 Chapter 3. Geometric Survey Parameter Errors Strue <j>m fa 51 m 0.203 20.0 50 m 0.171 20.0 49 m 0.152 39.7 ti 20.0 20.0 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> is the recovered model norm, fa is the data misfit achieved by the inversion and (f> is the target misfit. These values correspond to results plotted in Figure 3.11 m d 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 59 Chapter 3. Geometric Survey Parameter Errors -0.5i (a) Sir 30 20 10 True Model Recovered Model 00 !»' o o = 51 m. ?-10| CD O < (/)• co 20 — 9-301 -2.5 -0.5. C •Strue = 50 IPPred. x x QObs. -40 -50-; 10' -- QPred. 1 50 100 Depth (m) 150 10 10 10" 10' 10" Freguencv (Hz) 30 20 10 True Model Recovered Model () IP Obs. (d) « iE-101 «oo m. IP Obs. co -2.5! 50 100 Depth (m) 150 S -201 — IPPred. cd g-30 x x Q O b s . -40 — Q Pred. -50' 10' 10" Frequency (Hz) -0.5i True Model Recovered Model (e) Sir 0) £ = 49 m. -2.5 50 100 Depth (m) 150 10 10 Frequency (Hz) Figure 3.11: Results from the inversion of large separation HLEM data with an s value of 50 m when s e 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. trU Chapter 3. Geometric Survey Parameter Errors that the target misfit was achieved for s true 60 = 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 h est is smaller than h etTU 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 3.3.2 61 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 parameters 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 attempting 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 4.1 63 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 specifically to details that relate to either AEM or HLEM data I will use a superscript A or a superscript H respectively (e.g. d will represent the forward modeled AEM data). HowA ever, outside of the discussion in this section, the methodology will not specify between the inversion of A E M 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(a ),p] , (4.1) T M 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 m A = [log{<Ti),log(<r ), M h] , (4.2) T which is used when inverting AEM data, and m H = [log(<Ti),log((T ), M s] . T which is used when inverting HLEM data. (4.3) Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 4.1.2 64 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(d ), Im(d )] , (4.4) T N N 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 i th measurement in a given data set to be Re(di) = Re(Fi[m}), (4.5) Im(di) = Im(Fi[m]), (4.6) and where the details of the forward modelling F{ depend upon the survey being considered. In the A E M case I denormalize the data from Equations 3.5 and 3.6 and substitute h for h to get the following expression for the data from the i th true A E M measurement Re{di)^Re{H {^ ;a,h)\ (4.7) s z Si and Im(df) = Im{H { , s z Ui S i ] a, h)). (4.8) Since the estimated measurement height h est inclusion of h true has no direct effect upon AEM data, the does not change the forward modelling calculations. Similarly, by denormalizing Equations 3.11 and 3.12 substituting s for from the i th s e trU the data HLEM measurement will have the form Re(d?) = Re{Hl{^ h,-a, s)) + H (s) - H (s ), P P ett (4.9) Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 65 and Im(df) = Im(H (ui,hi] s z 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 ij th sensitivity of the i th datum to a change in the j th entry in J represents the model parameter as defined by As a result the generic sensitivity matrix will have the form J=[J„J ], (4.12) P where is the JV x M conductivity sensitivity matrix and J is the JV x 1 parameter p sensitivity matrix. From the expressions for both d and d it is clear that the subsurface conductivity A H 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 E M surveys I am dealing with have been determined in Zhang & Oldenburg (1994) using the adjoint Green's function method. The sensitivity of the i th datum to a change in the conductivity of the j th ddi d(Tj iuipomx layer can be expressed as J™ 0T ^ )°( ) ' - 2 T T +1 dz j complex A5 A3dA (-) 4 13 Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter where Zj is the depth to the top of the j th layer, 66 is the tangential component of the electric field in the Hankel transform domain, and J is a zeroth order Bessel function. 0 However, the model parameter I am concerned with is log(o-j) therefore, the sensitivity value I want is denned by (4-14) J« - - jr*. o-j do-j These values will fill J which makes up the first M columns of the sensitivity matrix. a 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 A E M case I must determine the sensitivity of the data to a change in measurement 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 d in complex form as A d? = -^ R \ J {\ )d\. (4.15) 3 4TT TE J 0 Si U 0 0 By differentiating Equation 4.15 with respect to h I can calculate the sensitivity of the i th datum to a change in measurement height. The complex sensitivity is equal to ^ = ^1°° ~ e 2U0h W J (A«)«tt. (4-16) 0 These values will fill J when inverting AEM data. p 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 . ) « The sensitivity of the i th - + ( ^ ) . (4,T) 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? ds m T Air r e- 2uohi n ,3 d , „ . T xx „ m T d ( 1 Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 67 Expanding the partial derivative in the integrand I get = £[J-i(A«) - (4.19) MXs)}. Using the identity J-n(As) (-1)" (4.20) J„(A«), I am left with dJo(Xs) ds (4.21) AJi(Aa), which can be substituted back into the integrand to give the complex sensitivity value (4.22) + These values will fill J when inverting HLEM data. p 4.1.4 The Objective Function Since I have modified the model it is necessary to modify the objective function accordingly. 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 <P = <Pd + (4.23) Pfa where <pd is the data misfit, <p is the model objective function, and /3 is the regularization m 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 <p as m (4.24) Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 68 where 4> is the conductivity model norm, (p is the parameter norm, and a is weighting a p p parameter which defines the relative importance of <p within the model norm. The conp ductivity norm determines the character of the conductivity structure I want to recover, and is defined as 4>« = \\W {log(a) - log(a ))\\\ v (4.25) Tef where W is the conductivity weighting matrix and <r / is the reference conductivity a re model. Wg- is chosen to be identical to the model weighting matrix defined in Equation 2.43. The parameter norm determines the closeness of the recovered parameter to some reference value p f, a n re d is defined as ^={P-Preff. (4.26) Using the model m from Equation 4.1 it is possible to express <6 as m (p = \\W {m-m )\\ , 2 m m where W m ref (4.27) is equal to W a 0 0 (4.28) ,/a 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) and how this choice p 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> - In practice however, this is computationally expensive. d 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) . These two algorithms will be referred to d 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>(m ) = cj> +f3<f>^ , k +1 k+1 , 1 d (4.29) and the objective function corresponding to the model from the previous iteration <Km ) = 4> + 3<ft. (4.30) h k d r If <j)(m i) < </>(m ) then 8m is accepted since it decreases the global objective function. k+ k However, when </)(m i) > (f>(m ) the step 8m should not be taken. Therefore, the step k+ k 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 m i is updated and <f>(m i) is recalculated and compared to <f){m ). This k+ k+ k process is repeated until a 8m, that decreases the objective function, is found. When an acceptable m +i is found the algorithm must decide whether to continue on to the next k iteration or to accept the model m i as the solution to the minimization problem. The k+ stopping criteria used in the algorithm are (1) the relative gradient of <j>, i.e. is m i at a k+ minimum, and (2) the relative change between m i and m , i.e. is the model stationary. k+ k The relative gradient of <f> is considered small when \gW(m ) k+1 max m^ \ +1 i = l,...,M <£("ife+i) where g(m i) <e u (4.32) is the gradient of (j> and t\ is a tolerance level. The model is considered k+ to be stationary when ||m - m \\ max{\\m i\\, ||7nfe||) fc+1 k+ k < e , 2 Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 71 where e is a tolerance level. The algorithm is said to have converged when either of 2 these conditions are satisfied. In my algorithm both of the tolerance levels are set equal to 10 . With all of these pieces explained in detail it is now possible to map out the -3 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,t , and m f, and (3) parameters needed to define the objective aPt re functions such as a and a , and error estimates for the data. The value of fi must 3 z also be provided. 2. The initial values of d, <pd, <p , and J are calculated using m m, . tart 3. Iteration counter k is set equal to zero and the variables calculated in Step 2. are labeled to belong to the k th iteration as m.fe, d , <p , (p , and J(mfe). h d m 4. Begin k + 1 iteration. th 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 User Input including P v l (2) " I Calc initial values I Set k = 0 Initialize Variables) L3J • ( (4) Begin iter, k+1 ) i ' ^ 5 . U s e fixed value of P I ;, ^ Calc G N step 8m Setm =m +8m (k+1) 00 I Calc k+1 values ( | A ) Use D G N <;8m = co 8m and m =m + 8m > ,k+n °A N1 ( T ) Has § decreased? - (k, Set k=k+1 Update variables N 9 Convergence ? 10 Exit Figure 4.1: Flowchart for the fixed (3 inversion algorithm with a fixed value of 73 Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 7. The values that were calculated in Step 2. are calculated using m +\, to be k <f> , 4> +\ and k+1 k d , k+1 J(m ). k+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 A E M sounding data with a fixed value of (3. The value of /? was selected, through trial and error, such that the goal fa = <f> was d 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) , (3 and <f>, as a function of iteration. It is shown that at each iteration m 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>* . The best plan is d 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 3000 True Model ^2500 Recovered Model (a) (b) o o IP Obs. CL CL — T3 X t'2000 CO IP Pred. X o 1500 X QObs. Q Pred. -o c to gioool CQ S SI Q. £ / // X 500 X' -3 50 100 150 Depth (m) 10' 10" 10' Frequency (Hz) 10° Figure 4.2: Recovered models and predicted data from A E M example fixed (3 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 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 (3 k+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>*} ^ is equal to the target for the inversion (j) . Therefore, at k+1 d each iteration the target <f\i ^ is defined as k+1 (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>*} \ the curve of k+1 fa versus (3 usually has a parabolic shape. Therefore, there may well exist two values of (3 that give fa = (j>*} ^• In this case the smaller of the two (3 values should be chosen k+ 75 Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 10 *l 1 . 2 3 . 4 . Iteration 5 . 6 1 1 0 'l 1 2 • 3 • Iteration 4 • 5 1 6 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) (p (circles) and target 4>d (crosses) vs. iteration, (b) c6 vs. iteration, (c) fixed f3 vs. iteration, and (d) (p(mk+i) (circles) and <p(mk) (crosses) vs. iteration. d m Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter as fi k+1 76 such that the least amount of structure is added to the model. In the case that it is not possible to find afivalue which satisfies (pd = <p*j ^ the value offiwhich gives k+1 the minimum possible <j> value is chosen as fi . k+1 d After having found fi k+1 the Gauss-Newton step Sm is calculated using Equation 2.60 and the model is updated to be m +i = m-k + Sm. At this point it is necessary to check if k the step Sm is beneficial. Sincefiis changing at each iteration, the <p values from iteration k, evaluated with fi , can not be compared with those from iteration k + 1, evaluated k with fi . Therefore, in order to evaluate the changes taking place in cp at successive k+1 iterations both (p(mk) and <p(mk+i) are calculated with fi = fi k+1 such that the values are defined as 4>(m ) = cMm ) + / 3 V m ( m ) , (4.34) fc+ k fc fc and <j>{m ) = (p {m ) + k+l d k+1 fi <p (m ). k+1 m k+1 (4.35) In practice, it has been shown that the choice of fi , corresponding to a small k+1 decrease in <pd, can result in a solution that both minimizes the global objective function and satisfies (pd = (p* (Parker, 1994). However, this is not always the case. When (p* d d is set too small then the target misfit is not attainable and the algorithm can reach a stage at which a decrease in (p can only be achieved by choosing small fi values that add d 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 m from the previous iteration. k Requiring c6(mfe ) < 4 ( k) is a consistency check on the algorithm. It at least > m +1 ensures that the previous model is not a better model than the updated model m +i — m + Sm. If <p(mk+i) > (p(mk), then a number of actions are possible. First the k k 77 Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 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) = 4> ). Empirically, the case that (j)(mk+i) > 4>{mk) occured when the misfit { k) 1 d was beginning to plateau out and could no longer be significantly reduced, even with a large perturbation 8m. The achieved misfit </> was inferred to be a representation of mm the minimum misfit. Accordingly (j) should Ue above this value. Therefore, when the d algorithm detects that <^(mj. i) > <f>(mk), the value of <j) is increased and the inversion d + is restarted. In my algorithm the restart value of <f> is set equal to 1.1 x <^> . mm d 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 m f, and (3) parameters needed to define the objective re functions such as a, and ct , and error estimates for the data. The target misfit <j> z d must be provided. 2. The initial values of d, fa, (f> , and J are calculated using m m, . taTt 3. Iteration counter k is set equal to zero and the variables calculated in Step 2. are labeled to belong to the k th iteration as m , d , <j> , (fa, and J(mfe). k k d Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter Discrepancy Principle Inversion Algorithm for fixed value of p 78 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 J (k+,) I r~T\ Calc GN step 8m ' Set m = m + 8m w (k+1) (k) I ! (Z) J5A^ Choose larger <()* and restart 9A Setk=k+1 Update variables Calc k+1 values ) <jjf) Has <)| decreased? i>[ C +P C <^ +B +1) (k+1) ) ) (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 4. Begin k + 1 iteration. Calculate (f> th ' using Equation 4.33. +1 d 5. The f3 value which satisfies fa = <j)*} ^ is selected as /3 k+1 6. The Gauss-Newton step 8m, corresponding to /3 fe+1 fe+1 . , is calculated using Equa- tion 2.60. The model for this iteration is defined as m k+i = m + 8m. k 7. The values that were calculated in Step 2. are calculated using m i, to be k+ <f> , <fft\ and J(m k +1 d fc+1 79 d , k+1 ). 8. The values of 0(mfc ) and <j>(m ) are calculated using Equations 4.34 and 4.35. If +1 k 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> . Then the inversion is restarted d 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 A E M 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} , /3 and <j>, as a function of iteration. m 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 — True Model — Recovered Model E ( ) ^ a > o (b) F T 3- n c o O o Depth (m) 150 10" 10 Frequency (Hz) Figure 4.5: Recovered models and predicted data from AEM example discrepancy principle 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 ct and p f- A discussion of the ways in which each of these values affect p re 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 ct and p f- The discussions in this p re section will deal strictly with the discrepancy principle inversion algorithm. 80 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/> vs. iteration, (c) /3 vs. iteration, and (d) </>(mfei) (circles) and (f>(mk) (crosses) vs. iteration. m + Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 82 Choosing <x p In order to choose a good a value it is necessary to understand how this choice affects p the problem. The value of ot determines the relative importance of the parameter norm p term. Therefore, it controls how strictly the requirement that the recovered p is close to p f will be enforced. To get a feeling for the effects I set a equal to its two extreme re p values: 0 and oo. As a approaches infinity, p is fixed at the value of p f. However, p re when a is equal to zero, p is not required to be close to the p f value at all. Looking at p re how the problem behaves at these two extremes provides insight into choosing the value of a . p I begin by looking at the case where a approaches infinity. Since it is not mathemp atically possible to set a = oo I simply fix p = p f for each inversion. For each of the p re synthetic examples I carried out a series of inversions, each with a different p f value. re 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> value, m 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 h true = 30 m. The data were inverted 11 times with h fixed equal to different values of h f ranging from 45 m re 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 <6 values from each inversion. These values are m equal to <$> since <6 is equal to zero as a result of p being set equal to p f. The middle a P re 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 l o g o Measurement Height (m) 1Q 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/> value, middle panel is final <f>d value, and the bottom panel is the recovered conductivity models plotted next to one another. m These results are identical to those described in Section 3.1. When the h value used by the inversion is smaller than h e, trU 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 h ue, tr the recovered model has a high conductivity structure at the surface, and, as h increases, it becomes impossible tofitthe 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) values that correspond to this subset, m h — 33 m to 15 m, it can be seen that there exists a minimum (f> value that corresponds m to a particular value of h. For this example, this value of h is approximately equal to 27 m which is close to h etrU 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) value, middle panel is final fa value, and the bottom panel is the recovered conductivity models plotted next to one another. m case I chose to invert the data that had an s e value of 11 m. Eleven inversions were trU carried out with s f values ranging from 12.5 m to 7.5 m. For each of the inversions re the data were assigned an error of 0.5% of the primaryfield,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 s true = 51 m were inverted eleven times with s f values ranging Te 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> value, middle panel is final value, and the bottom panel is the recovered conductivity models plotted next to one another. m 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 s . true 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 = 0 is a good choice. p Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 86 Instabilities when oc = 0 p While it seems that setting a = 0 is a good idea it can cause problems during the first p 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> return to their starting values. By m decreasing fi slowly the problem is regularized and any poor steps that may have been taken are avoided. However, if a is set equal to zero, 4> will not be regularized by fi p p 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 a = 0. In order to do this I will redefine the terms in Equation 2.60 in p 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 J WjW J T d JJ T becomes J J T r L which can denned explicitly, using Equation 4.12, as JTJ o- °~ J J </<7 y — p . J p °<r J o~ J u p jTj P J (4.37) °p 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 WjW a 0 0 (4.38) 87 Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter Substituting all of these details into Equation 2.60 we can write out the Gauss-Newton step in matrix form as jjj. (w?w , JJJ \ +/3 P jTj jTj 0 °p °<r p °pj u 'jf a WfW a Sd-j3 Slog(<r) | 0 < 0 Sp E 0 (4.39) log(a) - a„ log(a ) ref V - Vref where Sd is equal to d° * — F[m]. b For a = 0 and /3 —> oo Equation 4.39 becomes p WjW, 0 0 a Sp% P J -W?W Slog(a) {log(o-) - log(a )) ref (4.40) J Sd T such that the conductivity and parameter perturbations are decoupled and the solution has the form Slog(o-) = - {log(a) - h = (J T P Jp)" 1 (4.41) log(a )) ref (4.42) Jjw. From Equation 4.41 it is clear that Sa is equal to zero when log(cr) = log(o~ f). The re 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 a = 0 can be shown. p For this example I chose an h f value of 42 m and denned the rest of the details re 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 a = oo and the a = 0 cases. Figure 4.10 shows plots of the p p 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 A E M data. The starting height was set equal to 42 m. The two fines in the plot correspond to values of a = oo (circles) and 0 (plusses). (a) <f>& vs. fi, (b) vs. (3, (c) h vs. f3, and (d) (f> vs. fi. p Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 89 The solid lines with circles refer to the case where a = oo and the dashed lines with D plusses represent ct = 0. p When ct = oo the problem behaves as expected. At large (3 values the fa curve p 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 m f and therefore, fa is approaching its starting value. However, setting re a = 0 causes a change in the behavior. The fa versus (3 curve still approaches zero as p (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 ct = 0 curves). Therefore, the linearized Gauss Newton step Sp is p too large and it is necessary to introduce some additional regularization to ensure that the step is acceptable. A simple solution is to set a equal to some non zero value. p When a is non zero and /3 — > oo the Gauss-Newton step will have the form p WjW, 0 Slog(a) h Sp '-W^W (log{a)-log(<r )) a ref (4.43) 0 -(P-Pref) and the solution will be Slog(cr) = - (log(o-) Sp - log(a )) (4.44) ref -(p-pr f), (4.45) e such that when log(o~) = log(o~ f) and p = p f, Sp will be equal to zero. It is important re re 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 f3a . Therefore, any non-zero value of a will result in Sp —> 0 p p as f3 gets sufficiently large, but the point at which this occurs depends upon a . This p 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 A E M data. The starting height was set equal to 42 m. The three lines in the plot correspond to values of a = oo (circles), 10 (plusses), and 10 (crosses), -6 -3 p (a) 4>d vs. fi, (b) 4><r v s - fit ( ) h vs. fi, and (d) <f> vs. fi. c where ot = oo, the dashed lines with plusses represent a — 10~ , and the dotted lines 6 p p with crosses represent a = 10 . From the curves in Figure 4.11 it is clear that as a -3 p p increases the value of fi at which 8p goes to zero decreases. For the algorithm to be stable it is necessary for a to be non zero at early iterations. p However, it has been shown that a should be equal to zero in the late stages of the p inversion. Therefore, we propose a cooling schedule for a . This entails decreasing a p p 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 a provides the algorithm with the stability p 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 PrefChoosing p ref During the first few iterations the inversion algorithm will demand that p be close to p f therefore, it is a good idea to have a realistic estimate for its value. The problem re of selecting a value for p f has some inherent difficulties associated with it. While, re both the A E M 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 p f- The fact that an inversion re algorithm is being used to recover the true parameter value indicates that it has already been assumed that p may be incorrect. Therefore, it may be useful to find a better eat estimate of p fre A simple choice of p f is to find the value that minimizes (f>d for the reference conre ductivity model a j provided by the user. Thefirststep infindingthis value of p is to Te fix the conductivity model equal to <r f and forward model data for a range of p values. re 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 /. This estimate of p f will depend upon the definition of cr f. In the case that the re Te re reference conductivity model is a good first order approximation to the true conductivity model, the choice of p f should be accurate. However, when the guess of a f is not as re Te good the results may vary. I investigated the influence of cr f on the choice of p f for re each of the three synthetic examples. re Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter -0.5i 10 — — (a) 92 * a = 0.1 S / m a = 0.01 S / m a = 0.001 S / m () I b AEM -2 3 TJ c o True o — — 0-2.5 tj = 0.1 S/m o = 0.01 S/m o = 0.001 S/m •2 -3 0 20 40 Measurement Height (m) -3.5, 60 50 100 Depth (m) 150 Figure 4.12: Results from the search for a best fit p t value for the AEM sounding data where h = 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) . re tTUe 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~ f models (0.1 S/m, 0.01 S/m, and 0.001 S/m) Te are plotted in Figure 4.12(b). It can be seen that when a f re = 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- j is varied by one order of magnitude in either direction the results are not as good. re In fact, the o~ f — 0.1 S/m (solid line) and the o~ f — 0.001 S/m (dotted line) have re re 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 93 Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter — o = 0.1 S/m —- o = 0.01 S/m o = 0.001 S/m (a) -0.5, (b) I£-1.51 Small HLEM TJ C — — — o 0-2.5 -31 3 8 9 10 11 Coil Separation (m) -3.5 12 (c) Large HLEM o = 0.1 S/m o = 0.01 S/m o = 0.001 S/m (d) 150 ? " 1 55 J-1.5. > o >N~ 1 — — * * True o — o = 0.1 S/m — o = 0.01 S/m o = 0.001 S/m o 1 0 0-2.5 CD 2 10 100 -0.5i 10' 2 ra ra Q 50 Depth (m) 10 — — • True o o = 0.1 S/m o = 0.01 S/m o = 0.001 S/m 48 49 50 51 Coil Separation (m) -31 -3.5, 52 50 100 150 Depth (m) Figure 4.13: Results from the search for a best fit p j value for the HLEM examples. Small HLEM sounding data where s e = 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 3 = 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). re trU true sounding I attempted to estimate s using the noisy data generated with s tTue = 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 frue 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 a f = 0.01 S/m. However, each of the minimum values he re 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 individually, this method of estimation appears to work best when the parameter dominates the response. Otherwise it seems that a poor choice of <r f could result in an incorrect re 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 cr f choice, but re also to the minimum with the smallest to get an accurate estimate of pt correspond to the minimum rue value. This suggests that it may be possible by searching for the pair of p / and a f values that re re value. In order to investigate this possibility the data misfits for a range of <r f halfspace values and p can be calculated and contoured over re the <r f and p values. re 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 a f halfspaces from 0.001 S/m to 0.1 S/m, and a range of h values from 15 re 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 tr / 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 re 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 re / 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 <r f value. re 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 <r / halfspaces from re 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 <r / and p f, can be Pe re 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 p f. The decision re Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 97 to discard the calculated a f value is based on the fact that while the estimate of <r / re re may be a good one, since the a f value provided by the user was most likely chosen for re 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 p f value could easily be included in an inversion algorithm, for my purposes I re have decided to always use the user supplied value of cr /. re This section has provided a reliable way in which to choose p f and a stable way Te in which to define a , and hence it is now possible to define the inversion algorithm to p 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 oc that affects Steps p 1 and 9A, and the inclusion of the step in which p f is determined as described in re 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 ct . p 2. The starting p f value is calculated as discussed in Section 4.2.2. re Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter Fixed (3 Inversion Algorithm with p included in the model 98 ( • C D User Input including ot J P Calc initial values including p, r. 2 e (3) Setk = 0 Initialize Variables "•I't 4 ; Begin iter, k+1 ' i ( ^5) Use fixed value of P I Calc GN step 8m ' Setm =m +8m w (k+1) 7 8A Use DGN 5m = co 8m and m = m + 8m (U1, N ] ($) Has < > | decreased?: (k, Set k=k+l Update variables oc = 0 a P Calc k+1 values I T L 9A (k) N / ® Convergence ?!;j) P 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 3000 True Model (a) ^2500 Recovered Model (b) o o IP Obs. — CL CL £ '2000 TJ CO 3 o1500 IPPred. / x QObs. x x //' —- Q Pred. A* TJ C CO 8>1000| CO £ 500] •X' 50 100 Depth (m) 150 j?; 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 a is decreased by a factor of 10 after each iteration. p This algorithm was used to invert the AEM sounding data with afixedvalue of j3. The value of /3 was selected such that the goal fa = <f> was achieved. The value of a was d p 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 thefixedp 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 h , it true does recover a model that can fit the data, has minimum structure, and reproduces the true model quite well. Therefore, the inversion usingfixedf3 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 o-e if new (f) *--* 0 old o c •10' CD * > CD '5 ca W -Q I. A A A . o CD 2 10 4 2 Iteration 4 Iteration Figure 4.17: Convergence curves from A E M 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 h (crosses), (d) fa - iteration, (e) fixed (3, and (f) 4>(mk+i) (circles) and ^(mjt) (crosses) vs. iteration. v s trtie Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 101 Details of discrepancy principle A l g o r i t h m 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 ct that affects Steps 1 and 9A, and the inclusion of the step in p which p f re 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 E x a m p l e s 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. equal to 1.0. The value of a p was set 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 has ensured that p is not able to wander too far from p f during p re Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter Discrepancy Principle Inversion Algorithm with p included in the model Ji) User Input including a 102 P -('(2) Calc, initial values . including p r. J re ((Jy I, Setk = 0 ^ 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 =m +8m ^ (k+l) (k) I 'JD (8A)Choose larger $ and restart v9AJ Setk=k+1 Update variables Op = 0 cx ^ j Calc k+1 values j c "+r" c° < $ / Q ( (jT) Has (() decreased? 7\ X Convergence ? P (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 3000i True Model Recovered Model (a) (b) o o IP Obs. — IPPred. ^2500| x x QObs. Q. -- QPred. Q_ £ 1 X ^2000| CO 3 O 1500 / "O c X (0 $1000| co Q. £ 50 100 Depth (m) 150 500 10' . X / 10' 10" Frequency (Hz) 10* Figure 4.19: Recovered models and predicted data from AEM example discrepancy principle 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 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). 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 A E M 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 conductivity model the inversion results from the fixed p and variable p, inversions of the A E M 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 (b) *--» Target Misfit 5 Iteration 10 15 o-e Recdh * - - * True h (d) 10 O-o ) new (f) o * - • * tiold <2io! 3 X . * o «10 ! o 8 8 < l -®-® O 5 Iteration 10 10 5 Iteration 10 15 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 h (crosses), (d) fa vs. iteration, (e) (3 vs. iteration, and (f) (j)(mk+i) (circles) and (f>(m,k) (crosses) vs. iteration. true Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 0 m :p. = 30m m :p = 27.17m rec true (a) E -0.5 CO r rec r e c r m, -1 105 :p. true true = 30m (b) r > •3-1.5 o c o O -2 o —- m :p, = 30m rec true — m :p = 29.88m r -2.5 rec r e c r m, -3. 50 100 Depth (m) 150 :p, true Hrue -3 =30m 50 100 Depth (m) 150 Figure 4.21: Comparison of models recovered from the inversion of A E M data with fixed h (dashed line) and variable h (solid line). The true model is also plotted (dotted line), (a) Models from the synthetic A E M data used throughout the chapter, (b) Models from a different set of A E M in which the model was more complex. Fix h Var. h h 30.0 27.2 4>m 2.50 x IO" 2.17 x IO" 2 2 4> d 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, ci is the recovered model norm, fa is the data misfit achieved by the inversion, and (j)* is the target misfit. Values are tabulated for both the fixed h (Fix h) and variable h (Var. h) inversions. m d 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 h . true 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 confirmed by comparing the <fi m 106 values from Table 4.1. These types of distortions should be expected in the A E M 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 A E M 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/» value. m Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter h Fix h Var. h 30.0 29.9 (f>m 2.50 2.49 x x fa 10" 10~ 1 2 107 ti 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) is the recovered model norm, fa is the data misfit achieved by the inversion, and <j> is the target misfit. Values are tabulated for both the fixed h (Fix h) and variable h (Var. h) inversions. m d 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 a was p 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 Depth (m) 108 Frequency (Hz) Figure 4.22: Recovered models and predicted data from small separation HLEM example 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 10 o-o (a) Recovered Misfit (b) * - - * Target Misfit co Q 10 5 Iteration 15 10 11.01 1 (c) 10.99! .10.98 1 * ] f — - * « * o-o Reed s *~* True s (d) ooo 10 Iteration 15 Iteration Figure 4.23: Convergence curves from small separation HLEM example discrepancy principle 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 s e (crosses), (d) fa vs. iteration, (e) (5 vs. iteration, and (f) <f>(mk+i) (circles) and 4>{mk) (crosses) vs. iteration. tTU Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 0 50 100 Depth (m) 150 10 2 10 3 10" Frequency (Hz) 110 10 s 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 separation 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, quadrature - 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 wasfixed(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 Parameter111 Iteration Iteration Figure 4.25: Convergence curves from large separation HLEM example discrepancy principle 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 (crosses), (d) fa vs. iteration, (e) fi vs. iteration, and (f) <£(mfc i) (circles) and 4>(m-k) (crosses) vs. iteration. P u e + Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 112 ij5 2 0 0 I \ -100 §"200 • ••••• I I • -3 : • 00 9200 9400 9600 N o r t h i n g (m) 9800 iog a 3.5 10 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 - h , . e 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 h est 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 methodology 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 4.4.2 114 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 distribution 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 (a) W ^ 20 N 0 ^ IP Obs IP Pre i**xxxx-xxx XL QObs r-_20< O Q Pre 9=-40. (b) 50 100 150 200 Easting 250 300 350 50 100 150 200 Easting 250 300 350 50 100 150 200 Easting 250 300 350 50 100 150 200 Easting 250 300 350 d 20 ? x 0 ^-20' 0 O £-40 5? w 115 1 20 0 X CVJ-20 O Sl-40 (d) °S 0 20 N X XL CO io_20 o l 0 ir-40. 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 CO o 110 • §10° >' 2 1*1 1 t ;, \A 1 , - - CO 0 Figure 4.29: Inversion top panel: D a t a misfit recovered (dashed line) models plotted next to 100 _ 200 Easting (m) • Vv_y'- 300 , ,- 116 ii i ij l o 9io° results from Sullivan D a t a w i t h s included i n the model. T h e values (dashed). M i d d l e panel: T h e estimated (solid line) and measurement height. B o t t o m panel: T h e recovered conductivity one another. Chapter 4. Recovering 1-D Conductivity and a Geometric Survey Parameter 117 for both A E M 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 A E M 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 representation 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 G C V 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 parameter 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 geophysical 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 119 Chapter 5. Application of GCV to 1-D EM Inversion 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 o6i , 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 d obs to the correct amount. I will begin as usual by denning a global objective function, <f>(m) = fa +fi<j> = \\Gm - d° '\\ b m 2 + fi\\W m\\ 2 m (5.1) where fa is the data misfit, fi is the regularization parameter, <j> is the model norm, G m is a linear forward operator, and W m 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, GCV{j3) [trace(I - G(G G T It has been shown that the + (5.2) 3W^W )- (P')] ' 1 2 m r 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 d '. ob When each datum within d oha 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 nonlinear objective function is denned as <t>{m) = cf> + (3(j> = \\F[m] - d° \\ b> d m 2 + f3\\W m\\ , 2 m (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 knowingfi*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 m . My goal k now, is to find the perturbation 8m to the model such that m \ = m -\-8m will minimize k+ k the objective function <f>. By expanding the operator F in a Taylor series I am left with F(m ) = F(m + Sm) = F(m ) + J(m )8m + Af(0{Sm )), (5.4) 2 k+1 k k k where J(m ) is the Jacobian defined at m and N represents the higher order terms of k k the expansion. Neglecting higher order terms, the data misfit fa can be approximated by the linearized misfitft™which has the form fa = \\F(m ) + J(m )m lin k k - J(m )m k+1 k - d \\ . oh> k (5.5) 2 This results in a linearized objective function of the form n (m ) = fa lin fc+1 where r = d oba k +ficl> = \\J(m )m + J(m )m m k k k — F(m ). k - r \\ + 2 k+1 k fi\\Wm \\ , (5.6) 2 k+1 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 tofitthe noise associated with d obs properly. After selecting fi using GCV, the perturbation to the model is calculated using Equation 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(6m ) they may be sufficiently large such 2 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 A E M 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 E M 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, -•' \A•) Begin iter, k+1 i —i v CD Estimate p ' using GCV >(k+ ) J ^ Calc GN step 8m Setm =m + 8m (k+,) (k, I *^CZ) <8A; Use DGN 5m = o) 8m and m = m + 8m ,k+1, 9A lk, Setk=k+1 Update variables N Calc k+1 -values •'"sT) +Has (j) decreased? i) (k+,) i) c +p c<c+rc v .N/® i) T Convergence? 1 10) Exit Figure 5.1: Flowchart for the GCV inversion algorithm with a fixed value of p. 124 Chapter 5. Application of GCV to 1-D EM Inversion 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 f3 . k+1 8. The values of <p(m,k+i) and (j)(mk) are evaluated using /3 fc+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 A E M 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. Thefinal(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 125 Chapter 5. Application of GCV to 1-D EM Inversion True Model (a) Recovered Model o o (b) — x IP O b s . 0 IPPred. /* x Q Obs. QPred. ^15001 C CO SJ1000 « CL S 500 JZ 50 100 D e p t h (m) GCV True Model (c) Recovered Model 10' 10' F r e q u e n c y (Hz) 10 o o (d) Jo — x IP O b s . IPPred. x Q 50 100 D e p t h (m) / QObs. x Pred. 10 10 F r e q u e n c y (Hz) Figure 5.2: Comparison of recovered models and predicted data from A E M example discrepancy principle and GCV inversions with a fixed value of h. The measurement height was set equal to 30 m. Discrepancy principle (x ) 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)) 2 x 2 GCV 4>m 1020.4 257.5 fa 1.87 x IO" 2.35 x IO" 2 2 20.0 17.4 20.0 N/A Table 5.1: Comparison of parameter values from the A E M discrepancy principle and GCV inversions with a fixed value of h. fi is the final regularization parameter value, <£ 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 (x ) and GCV inversions. m 2 126 Chapter 5. Application of GCV to 1-D EM Inversion o—o Recovered Misfit (b) --* Target Misfit 10 4 6 4 10 6 8 Iteration Iteration 10' »,• 10 — G - G 0 new * - - » tj) old (d) O10 O CO n o — ! O 2 4 10 6 Iteration 2 4 6 10 Iteration 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) <p vs. iteration, (c) fixed (3 vs. iteration, and (d) c/>(mfci) (circles) and <£(mfc) (crosses) vs. iteration. m + Chapter 5. Application of GCV to 1-D EM Inversion Iteration 127 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> vs. iteration, (c) fixed /3 vs. iteration, and (d) <j>(mk+i) (circles) and 4>{m,k) (crosses) vs. iteration. m 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> and ft curves in the GCV results are caused by the fact that the m 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 afixedvalue 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 5.2.2 130 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 . These p 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 p f value is calculated as discussed in Section 4.2.2. re 9A. The value of a is decreased by a factor of 10 after each iteration. p Now having designed the algorithm it can be tested on a synthetic example. Synthetic A E M Example The same synthetic A E M data used in the fixed p example was inverted with p as part of the model. For this inversion the starting value of a was set equal to 1. Figure 5.7 shows p 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 GGV Inversion Algorithm with p included,in , • , the model 131 S) User Input including q J { s P 2 C a l c initial values including p f. J re I I (3) Set k = 0 - .1 Initialize Variablesy ' . ^ •I C4j Begin iter, k+1 ; I ; .5) Estimate p k l ) using GCV i I ^ w "n v L ' ' 8A Use DGN 8m = u) 8m and m = m + om 9A Setk=k+1 Update variables 0Cp = 9 oc N] N Calc GN step 8m j ^Setm =m +8m J (k+1) (k) Calc k+1 values j r+rc<c+r c (%) Has <}) decreased? i, 9 : ,) Convergence ? P 10 Exit Figure 5.6: Flowchart for the GCV inversion algorithm with p included in the model. 132 Chapter 5. Application of GCV to 1-D EM Inversion True Model (a) 0 0 IP O b s . (b) Recovered Model 0. — 052000 x CQ 3 x - - ° 1 5 0 0 0 IPPred. /r QObs. Q Pred. C CC $1000 cc sz CL £ 500 Of 50 100 10 D e p t h (m) 10 10 F r e q u e n c y (Hz) 3000i GCV True Model (c) 0 0 IP O b s . (d) Recovered Model — x ^1500| 0 IPPred. x - - /r QObs. /*/ Q Pred. C CO gjiooo CO JZ Q. S 50 500 100 10 D e p t h (m) 10 10 F r e q u e n c y (Hz) Figure 5.7: Comparison of recovered models and predicted data from A E M example discrepancy 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 (%) 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)) 2 x 2 GCV h 29.1 28.9 <f>m 1235.3 259.9 fa 1.84 x 10~ 2.26 x IO" 2 2 20.0 17.4 20.0 N/A Table 5.2: Comparison of results from A E M 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} 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 (%) and GCV inversions. m 2 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 G C V 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 134 Chapter 5. Application of GCV to 1-D EM Inversion io l J 0 • 5 • 1 no' 10 15 0 1 Iteration • • 5 10 I 15 Iteration Figure 5.8: Convergence curves from A E M 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 h (crosses), (d) fa vs. iteration, (e) /3 vs. iteration, and (f) (j)(mk+i) (circles) and <f>(mk) (crosses) vs. iteration. true Chapter 5. Application 135 of GCV to 1-D EM Inversion 10 o—o Recovered Misfit Target Misfit 2 Iteration (b) 4 10" o-o Reed h *---* True h 10 (d) o Z10 . Q. 10 10 10 o-o (j) new (e) *-•* < > j old (f) h (1 O «10' '* 1 .Q o O 2 3 ' & 10 4 Iteration 3 ® o 4 Iteration Figure 5.9: Convergence curves from A E M 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 h (crosses), (d) fa vs. iteration, (e) /3 vs. iteration, and (f) c/>(mfci) (circles) and c/>(m.fe) (crosses) vs. iteration. true + 136 Chapter 5. Application of GCV to 1-D EM Inversion , (b) it=l r J " " ' - A j 1 r-l' d- i — GCV Model DGN ndamp = 1 50 100 Depth (m) 150 (d) it=3 — GCV Model — DGN ndamp = 0 50 100 Depth (m) 150 (f) it=6 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 4>m P GCV 2231.7 1.44 x IO" 3 137 fa 3.10 x 10~ 2.61 2 8.0 12.5 ti 8.0 N/A Table 5.3: Comparison of results from the four frequency A E M example discrepancy principle and GCV inversions with a fixed value of A. /3 is the final regularization parameter value, (j) is the recovered model norm, fa is the data misfit achieved by the inversion, and 4>* is the target misfit. Values are tabulated for both the discrepancy principle (%) and GCV inversions. m 2 d 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 A E M 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 138 Chapter 5. Application of GCV to 1-D EM Inversion Figure 5.11: Comparison of recovered models and predicted data from the four frequency A E M example discrepancy principle and GCV inversions with a fixed value of h. The measurement height was set equal to 30 m. Discrepancy principle (%) 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)) 2 Chapter 5. Application of GCV to 1-D EM Inversion 139 10' o-o Recovered Misfit *--* Target Misfit (a) 10 10 Iteration 20 (b) 30 10 Iteration 20 10' § new * - - * ij> old o-o (c) * (d) 10 10 20 30 Figure 5.12: Convergence curves from four frequency A E M 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) <£ vs. iteration, (c) fixed fi vs. iteration, and (d) 0(mfc ) (circles) and <j>(mk) (crosses) vs. iteration. m +1 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> , fl, m 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 function 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 thefirstfew 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 141 Chapter 5. Application of GCV to 1-D EM Inversion it=l — GCV Model — DGN ndamp = 1 50 100 Depth (m) 150 GCV Model it=2 DGN ndamp = 2 50 100 Depth (m) GCV Model it=3 150 (f) DGN ndamp = 0 50 100 Depth (m) 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 GaussNewton approach was developed. The algorithm was shown to work well when inverting frequency domain EM soundings which consist of twenty data. The GCV inversion algorithm 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 developed 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 estimating 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, measurement 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 A E M 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 conductivity 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 nonlinear 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 environmental 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 Optimization and Nonlinear Equations: SIAM Philadelphia. [4] Ellis, R.G., and Shehktman, R., 1994, ABFOR1D k ABINV1D: Programs for forward and inverse modelling of airborne E M 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 mapping: 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 electromagnetic 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 theoretical and field estimates: Geophysics, 60, 374-380. [11] Norton, S.J., San Filipo, W.A., and Won, I.J., 1999, Reciprocity formulas for electromagnetic induction systems: Proc. of the Symposium on the Application of Geophysics 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 problem: 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 Symposium on the Application of Geophysics to Environmental and Engineering Problems, 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 Geophysics: 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 susceptibility 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 E M data: Geophysics.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Inversion of em data to recover 1-d conductivity and...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Inversion of em data to recover 1-d conductivity and geometric survey parameter Walker, Sean Eugene 1999
pdf
Page Metadata
Item Metadata
Title | Inversion of em data to recover 1-d conductivity and geometric survey parameter |
Creator |
Walker, Sean Eugene |
Date Issued | 1999 |
Description | 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 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. |
Extent | 10867030 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-06-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0053001 |
URI | http://hdl.handle.net/2429/9271 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1999-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
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
Cite
Citation Scheme:
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>
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