DETECTING UNEXPLODED ORDNANCE WITH TIME ELECTROMAGNETIC DOMAIN INDUCTION By Leonard Rodriguez Pasion B . Sc. (Honours) Mathematical Physics, Simon Fraser University, 1996 A THESIS SUBMITTED THE IN PARTIAL F U L F I L L M E N T OF REQUIREMENTS MASTER FOR T H E DEGREE OF OF SCIENCE in THE FACULTY OF GRADUATE DEPARTMENT STUDIES OF EARTH A N D OCEAN SCIENCES (GEOPHYSICS) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September 1999 © Leonard Rodriguez Pasion, 1999 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Earth and Ocean Sciences (Geophysics) The University of British Columbia 129-2219 M a i n M a l l Vancouver, Canada V 6 T 1Z4 Date: Abstract In this thesis I assume that the Time Domain Electromagnetic ( T E M ) response of a buried axisymmetric metallic object can be modelled as the sum of two dipoles centered at the midpoint of the body. The strength of the dipoles depends upon the relative orientation between the object and the source field, and also upon the shape and physical properties of the body. Upon termination of the source field, each dipole is assumed to decay as k (t + a)~^ e ' / . - 7 The parameters k, a, /3 and 7 depend upon the conductivity, permeability, size and shape of the object, and these can be extracted from field or laboratory measurements by using a nonlinear parametric inversion algorithm. A n investigation carried out using an analytic solution for a sphere and laboratory measurements of steel and aluminum rectangular prisms, suggest the following methodology. The value of f3 might be used as a diagnostic to assess whether the metallic object is non-magnetic or magnetic. If the object is thought to be magnetic, then the ratios of k\jk2 and /3i//3 are diagnostic indicators as to whether the geometry is plate-like 2 (uninteresting) or rod-like (a high candidate for being a U X O ) . ii Table of Contents Abstract ii List of Tables vii List of Figures viii Acknowledgement xiii Dedication xiv 1 Introduction 1 2 Development of an Approximate Forward Model 8 2.1 The Time Domain Electromagnetic Response of a Compact Metallic Object . . . 9 2.1.1 Maxwell's Equations 9 2.1.2 Constitutive Relationships 2.1.3 Boundary Conditions 2.1.4 Quasistatic Assumption 10 . 14 14 2.2 Electromagnetic Induction as a Diffusive Process 15 2.3 Approximating the Response of an Arbitrary Shape Target 19 2.4 The Time Domain Response of a Permeable, Conductive Sphere in a Uniform Field 2.5 20 The Decay Characteristics of the T D E M Response of a Sphere . 26 2.5.1 Decay Characteristics of B-field Response . . 26 2.5.2 Decay Characteristics of d3/dt Response 27 iii 2.6 3 Statement of Forward Model 28 2.6.1 Spatial Response Function M(r) For an Axi-symmetric Body 29 2.6.2 Time Decay Function L(t) 30 2.6.3 The Approximate Forward Model 33 2.6.4 Approximate Forward Modelling Examples 35 Non-linear Parameter Estimation Procedure 38 3.1 Defining The Objective Function 42 3.2 Variables of the Inversion Algorithm 44 3.2.1 Positivity Constraint 45 3.2.2 Scaling the Model Parameters 45 3.3 Newton's Method of Optimization 47 3.4 A n Algorithm For the Computation of Parameter Estimates . . . 49 3.5 Preliminary estimates of the model parameters 49 3.5.1 Analysis of Plotted D a t a . 50 3.5.2 Fitting A Best F i t Single Dipole at Each Time Channel 58 3.6 Computing a Search Direction 60 3.7 Computing a Step Length 64 3.8 Stopping Criteria 66 3.9 Error Bounds of the Parameter Estimates 67 3.10 Application to a Synthetic Data Set 69 3.11 Some Tests of the Parameter Estimation Procedure 73 3.12 Selection of Data i n a Field Survey 82 3.12.1 Time Channels 83 3.12.2 Spatial Coverage . 85 3.13 Summary . . iv 88 4 Relating Physical Parameters of a Metallic Target to Model Parameters 90 4.1 The Effect of Shape on the T D E M Response 90 4.1.1 Lab Setup 90 4.1.2 Computing B-field responses from dH/dt Data 92 4.1.3 T D E M Response of a Plate 93 4.1.4 T D E M Response of a Rod 98 4.2 Analysis of Laboratory Measurements of Metallic Prisms 4.3 Relating f3 to Magnetic Permeability 98 100 4.3.1 Variations of (3 with Magnetic Permeability when Using B-field Data . . . 101 4.3.2 Variations of (3 with Magnetic Permeability when Using dH/dt Data . . . 102 4.4 The Ratio *i/Jb 102 4.5 The Ratio / ? i / / 3 4.6 Analysis of Measured U X O 109 4.7 Proposed Algorithm For Determining Dimensionality of Target Ill 2 1 5 Field Application 5.1 5.2 6 115 Description Field Data Set 116 5.1.1 Target Information 116 5.1.2 Instrumentation 116 5.1.3 Survey Geometry 118 Characterization of EM61-3D Errors for Inversion 5.2.1 5.2.2 119 Method 1: Estimation of Noise B y Fitting Subsets of the Complete Data Set 5.3 0 2 123 Method 2: Estimation of Noise B y Using Measurements at a Distance as Noise Measurements 128 Preliminary Data Analysis 130 v 5.4 6 5.3.1 Estimates of Location and Orientation Obtained by Plotting the Data . . 130 5.3.2 Fitting A Single Dipole at Each Time Channel 5.3.3 Estimation of Time Decay Parameters 5.3.4 Summary of Preliminary Data Analysis 132 . 134 136 Inversion of the Field Data Set 137 5.4.1 Selection of Data 137 5.4.2 Inversion of Z- component oidTH/dt Data 5.4.3 Inversion of A l l Three (X, Y, Z) Components of Data . 138 142 5.5 Results and Discussion 145 5.6 Summary 155 Summary, Conclusions, and Future Work 158 Appendices 161 A Time Decay Parameters of Various Unexploded Ordnance 161 B 170 Time Decay Parameters of Metallic Prisms C Calculation of the primary field B p from a square loop References 173 176 vi List of Tables 2.1 Parameters for the dTl/dt response of a steel rod and a steel plate 3.1 Results from fitting a single dipole to some of the time channels of the synthetic 36 data set 59 3.2 Recovered parameters for the synthetic data set inversion 74 3.3 Errors in the recovery of parameters for a buried 75 mm anti-tank mortar . . . . 82 3.4 Recovered parameters when using different time channels 5.1 Results from fitting a single dipole to some of the time channels of the field data 84 set 135 5.2 A summary of the estimated decay parameters 136 5.3 Recovered model parameters when inverting 83/dt field data 142 5.4 Comparing expected parameters to recovered model parameters when inverting 5.5 dB/dt field data . 146 Results from applying the TJXO discrimination diagnostics 147 vii List of Figures 1.1 Typical U X O encountered in field surveys 2 1.2 The Geonics EM61-3D Time Domain Electromagnetic Sensor 5 2.1 The induction of currents in a buried confined conductor 9 2.2 The relationship between the magnetic field H and the magnetization M 2.3 Hysteresis curve for a ferromagnetic material 2.4 (a) A n arbitrary conductive body immersed i n a primary field B . (b) The field 15 2.5 The effect of permeability for a sphere in a uniform 2.6 The three stages of the T D E M response of a sphere 2.7 Set-up for a sphere in a uniform field B° 2.8 Graphical analysis to determine how many terms to keep in the infinite series 2.9 11 13 p step-off primary . . . . field 16 18 , 21 solution of a sphere 24 B-field dH/dt response of a steel and aluminum sphere 25 2.10 Comparison between the measured response of a metallic cube and the analytic response of a metallic sphere 26 2.11 Decay behavior changes due to variations in parameters k, a, (3, and 7 . 31 2.12 Fitting the time decay of a sphere using L(t) 32 2.13 Representation of the fields generated by induced fields i n a metallic target as a pair of orthogonal dipoles 33 2.14 Co-ordinate System for a buried target 34 viii 2.15 Synthetic 83/dt data for a steel rod and a steel plate using the approximate forward model '. 37 3.1 The 83/dt measurement taken over a 105 m m shell 39 3.2 Decay of the response at (X,Y) = ( l m , Ira) 40 3.3 The synthetic 83/dt response curves for a 75 m m anti tank mortar buried 58 cm beneath the surface 44 3.4 Values used for scaling the model vector m 3.5 The determination of strike (a) from plan view plots of the horizontal component t y p of the B-field data 3.6 46 51 The vertical component of secondary field for a thin rod and thin plate with various c 53 3.7 Estimating the depth and location of a thick steel rod 53 3.8 Definition of the "Full Width at Half Maximum" ( F W H M ) 54 3.9 Errors in the depth estimate of a buried cube when using the F W H M 55 3.10 Depth and location estimates for a steel rod and plate oriented at various dip . . 56 3.11 Approximating the response of an arbitrary compact object by that of the sphere 57 3.12 Results from fitting the Z-component of the response measured directly over steel rod (8 X 1 X 1) 58 3.13 Orientation estimate of a U X O by finding a best fit dipole for each time channel 3.14 The meaning of small and large eigenvalues in the Hessian 60 . 62 3.15 Adjusting the Newton Step using Directional Discrimination 63 3.16 X component of synthetically generated data for a 75 m m mortar 70 3.17 Y component of synthetically generated data for a 75 m m mortar 71 3.18 Z component of synthetically generated data for a 75 m m mortar 72 3.19 Time decay curves of the synthetic data set inversion 73 ix 3.20 Convergence results when maxstep = 1000 x ||m yp|| 75 3.21 Convergence results when maxstep — 10 x | | m t | | 76 3.22 The recovered parameters for different choices of e 77 3.23 Convergence behaviour when using various variable transformations 79 3.24 Convergence results when K = 2 80 3.25 Convergence results when K = 5 81 3.26 The error i n recovered parameters for a buried 75mm anti-tank mortar 83 3.27 Using only a portion of the time channels 85 t yp 3.28 Observed and predicted data for inversions using only a subset of the available time channels 86 3.29 Different survey configurations used to test how many stations are required for accurate model parameter recovery 87 3.30 Results from inverting only the zT-component of the dH/dt data of Line X = 0 m . 89 4.1 Geonics Laboratory Setup 91 4.2 Different target orientations for Geonics laboratory measurements 92 4.3 Eddy current induction in a non-permeable plate 93 4.4 Time Decay curves for 8x8x1 inch steel and aluminum plates 94 4.5 Eddy current induction in a magnetically permeable plate 95 4.6 Time Decay curves for 8x1x1 inch steel and aluminum rods 99 4.7 The behaviour of parameter f3 for various size spheres with varying permeability -101 /x 4.8 Relating the aspect ratio of a steel target with the ratio ki/k 4.9 Relating the aspect ratio of an aluminum target with the ratio ki/k 104 2 2 105 4.10 Relating the aspect ratio of a steel target with the ratio /3i//?2 107 4.11 Relating the aspect ratio of a aluminum target with the ratio ^1/^2 108 x 4.12 The T D E M response of an 81 m m mortar 4.13 The T D E M response of a 155 m m mortar 109 . 110 4.14 J3 averages for U X O 4.15 fci/&2 a n Ill d Pi 102 ratios when fitting to derivative of field dB/dt measurements of U X O 112 4.16 k and (3 ratios when applying parameterization to the B-field response . . . . . . 113 5.1 Information about a 105 m m shell .117 5.2 Setup for the York Field Data Set 118 5.3 Vertical component of dB/dt EM61-3D data taken over a 105 mm artillery shell at the York University Geophysical Test Site 5.4 120 Horizontal component of dB/dt EM61-3D data taken over a 105 mm artillery shell at the York University Geophysical Test Site. 5.5 121 Observed and predicted data when inverting synthetic dB/dt data from time channel <=1.788 msec 125 5.6 Noise model recovered in two synthetic test cases (Method 1) 126 5.7 Inversion of the t = 0.036maec time channel data acquired over a 105 m m U X O . 127 5.8 Noise model recovered for the field data by fitting a pair of dipoles to each time channel 5.9 128 Characterizing the noise of the EM61-3D data with e(i) = At B 130 5.10 Removal of data that fall beneath the noise level 131 5.11 Determining the strike of the buried 105 m m shell 132 5.12 Determining the depth of the buried 105 m m artillery shell 133 5.13 Determining the location of the buried 105 m m artillery shell 134 5.14 Estimation of time decay parameters by fitting the zT-component of the B-field and dB/dt response measured directly above the target xi 136 5.15 Surficial coverage of data used i n the inversions 138 5.16 Observed and predicted data when inverting the zT-component of the dB / dt field, plotted in plan view 139 5.17 Observed and predicted data when inverting the Z-component of the d3/dt field.140 5.18 Plan view of observed and predicted data when inverting all 3 components of dB/dt field, 143 5.19 Observed and predicted data when inverting all 3 components of dH/dt field . . 144 5.20 Comparison of forward modelling parameters obtained from lab measurements of a 105 m m U X O , to measurements made at a field site directly over a 105 m m 150 5.21 Forward modelling of parameters for a 105 m m shell • • • 152 5.22 Comparison of 63/dt responses of a 105 m m shell 153 C.l 173 A current carrying wire of finite length xii Acknowledgement First of all, I would like to thank my supervisor, Dr. Doug Oldenburg, for all his guidance and support throughout my work. His boundless energy, patience, and insight contributed immeasurably to this thesis. I would also like to thank Duncan McNeill and Miro Bosnar at Geonics L t d . for suggesting the problem and providing the invaluable lab and field T D E M data sets. It has been a pleasure interacting with the U B C - G I F Group members (Dr. Colin Farquarson, D r . Eldad Haber, Chrysophe Hyde, Francis Jones, D r . Yaoguo L i , Partha Routh, Roman Shektman, Sean Walker) on both a scientific and social level. We need to convince Doug to have more parties... Special thanks to D r . Cohn Farquarson for tolerating my more than daily visits to his office for advice; I am shocked he still hasn't hired security to protect himself from my constant badgering. Thanks to the fellow students in the Geophysics and Astronomy building for providing a wonderful and supportive environment for work. Special thanks to Kevin Kingdon and Stephane Rondenay for tolerating me as a cellmate i n the Mexican jail. Special thanks to my EOS Heavyweight Tag-Team partner, and second half of the heaviest office, Sean Walker for making everyday at office 227 a barrel of fun. I wouldn't have exchanged all the good times for saving the extra months we added on each others thesis. Finally, thank you Candice for always being there with your love and support. xiii Dedication o mom and dad. xiv Chapter 1 Introduction A n explosive ordnance is a munition that is either launched or fired with the intent of detonation at a specified target. Ordnance include small arms munitions, grenades, rockets, rifle grenades, and bombs (Figure 1.1). A n unexploded ordnance (TJXO) is an explosive ordnance that, due to some malfunction, remains undetonated. As a result, the ordnance can be found at the ground surface, partially buried, or buried at a depth of up to 8 m beneath the surface. Although they range i n length from less than 10 cm to several meters, TJXO have two identifying features. Firstly, they are generally axi-symmetric and rod-like in shape and typically a TJXO has a length to diameter aspect ratio of about 4. Secondly, most TJXO are made of steel. There exists 15 million acres of UXO-contaminated land in the United States (not including undersea U X O contamination) (Federal Advisory Committee, 1996). The majority of acreage consists of approximately 1500 different sites designated for weapons system testing and troop training by the Department of Defense (DoD). U X O are hazardous because they can explode when being struck or incinerated, or they can release toxic chemicals if the U X O casing has been compromised through impact or decay. The remediation of these sites has been made a high priority by the D o D in order to either maintain safe usage for continuing military operations or to permit land transfer to the private sector. Effective methods for clearing U X O contaminated land are crucial i n order to ensure public safety and to eliminate environmental concerns. The cleanup of UXO-contaminated lands have been hampered by labour, cost, and time intensive site remediation methods. The two main impediments to site remediation are the lack of an effective method of U X O characterization and the time and labour intensive methods 1 2 Chapter 1. Introduction Small Arms Munitions 62 cm are projectiles fired from a gun and are less than 10 c m long and less than 12 m m in diameter. \< 66 cm Mortars are made o f a metallic casing and range from 2.5 to 30 c m in diameter. H Projectiles H 47 cm range from 2.5 to 40 cm in diameter and from 5 to 120 c m in length. -142 cm -221 cm I Bombs 121 cm are metallic munition that range from 1 to 3 m in length. Rockets are self-propelled munition that range from 6 to more than 40 c m in diameter and from 30 c m to over 2 m in length. Rifle Grenades are fired from infantry rifles, and range from approximately 20 to 45 c m in length. Missiles are rockets except that they are guided by some controlling system. 99 mm 90 mm j cm 9.9 cm I JL Projection Grenades diameter. are most commonly 40 m m in 64 mm Hand Grenades 56 mm 64 mm are short-range munitions that are thrown. Grenades are made o f metal, plastic, cardboard, or rubber. F i g u r e 1.1: Typical U X O encountered in field surveys. U X O are generally projectiles that, depending on the fuse configuration, are designed either to detonate.on impact or after some specified delay time. Despite the wide range of U X O sizes, U X O are generally axi-symmetric, rod-like i n shape, and made of steel. The stabilizing fins found on some U X O are sometimes made of aluminum. Chapter 1. Introduction 3 required to excavate each suspect U X O . A t some sites the ratio of non-ordnance to ordnance items exceeds 100:1. Developing the ability to discriminate between ordnance and non-ordnance items is essential i n reducing the costs of cleanup. It has been estimated that a remediation program using current technologies to cleanup only 5% of the U X O contaminated land would require an investment of over 15 billion dollars and several decades of labour (United States Dept.of Defense, 1998). Several decades of research and millions of dollars have been invested towards developing U X O detection and identification technologies. This investment has yet to yield a system capable of efficiently and accurately finding and identifying U X O . It has been reported that approximately 70% of remediation costs are currently being used to excavate non-ordnance items (Butler et al., 1998). The traditional technique for U X O remediation is " M a g and Flag". This technique involves surveying the area with a metal detector and marking each location at which a buried object is detected with a flag. The term " M a g " is related to earlier times when a hand-held magnetometer would be the detection instrument. However, the term " M a g and Flag" is now used to describe the process of using any metal detector (e.g. electromagnetic induction sensors) without additional data processing that would discriminate between ordnance and non-ordnance items. In this case each signal must be excavated. The " M a g and Flag" method is thus severely hampered by an excessive false alarm rate. Current detection technologies for U X O detection include magnetometry, ground penetrating radar ( G P R ) , time and frequency domain electromagnetic induction, and infrared sensors. The effectiveness of each technique depends on the types of U X O , the depth at which the U X O is buried, and the examination area i n which the U X O is buried. These techniques are tested as part of the Unexploded Ordnance Technology Demonstration Program. This U . S . government funded program was established to set baselines of U X O detection technology by holding demonstrations by various geophysical and engineering contractors at the U . S . Army Jefferson Proving Ground i n Madison, Indiana. These 'competitions' consist of first detecting, Chapter 1. Introduction 4 then identifying, well characterized buried U X O within a set time frame. Prior to competition contractors are given a list of, as well as access to, all buried materials (both ordnance and non-ordnance). This situation has led to 'finger-printing' algorithms that match measured responses to a library of response type curves (references). For this thesis the method of time domain electromagnetic induction will be examined as a tool for detection and discrimination. Time domain electromagnetic ( T D E M ) induction surveys have been successful i n detecting subsurface ferrous and non-ferrous metallic objects, and these surveys are a mainstay amongst technologies currently utilized i n U X O clearance projects. In the T D E M method a time varying magnetic field is used to illuminate a conducting target. This primary field induces surface currents on the target, which then generate a secondary magnetic field that can be sensed above ground. The surface currents diffuse inwards, and the observed secondary field consequently decays. The rate of decay, and the spatial behaviour of the secondary field, is determined by the target's conductivity, magnetic permeability, shape, and size. Although the full T D E M decay curve can, i n principle, be measured, the primary T D E M instrument for U X O detection has been the E M 6 1 from Geonics L t d . , which measures only the vertical component of field at a single, late time. Recognizing that more information about a buried target would be obtained by measuring the full time domain behaviour of the three components of the secondary field, Geonics developed a prototype instrument EM61-3D (Figure 1.2). EM61-3D uses a square transmitter coil to produce a primary field. A t the center of the transmitter coil are three orthogonal receiver loops capable of measuring three components of the secondary field. The unit is mounted on a trailer, and pulled across the survey area by the operator. The three components of secondary field are measured between 0.037 milliseconds and 37 milliseconds after the termination of the primary field. The goal of this thesis is to investigate the possibility of using these data to discriminate between U X O and other buried objects. Also of interest to Geonics, is the importance, and need, of the horizontal component measurements Chapter 1, Introduction 5 Controller F i g u r e 1.2: EM61-3D Time Domain Electromagnetic Sensor. A l m x l m square transmitter coil produces a primary field. Three orthogonal receiver coils located at the center of the square transmitter measure three components of the secondary field. Chapter 1. Introduction of the field. 6 Building an instrument capable of simultaneous acquisition of horizontal and vertical component is difficult, and is preferably avoided. The goal of this research was to develop a method of interpreting T D E M data that has the capability of discriminating between U X O and non-ordnance objects. In Chapter 2 an approximate forward model is developed for describing the decay of a buried metallic object. This parameterization, which is based on the analytic solution of a sphere and the magnetostatic solution of a spheroid has a strong physical basis and is simple to compute. The time domain response of any metallic object is represented as a pair of perpendicular dipoles located at the center of the object. The strength of the dipoles decay with time, and the parameters that govern the time decay behaviour are related to the conductivity, permeability, shape, and size of the buried object. In Chapter 3 a non-linear parameter estimation procedure is developed to recover the 13 parameters of the model. Rough estimates of the model parameters can be obtained by applying some simple data preprocessing strategies. These initial estimates are then iteratively refined by solving a non-linear least squares problems. Newton's method is used to minimize the sum of the squared differences of the predicted and observed data. The simplicity of the approximate forward model allows the use of an exact Hessian matrix. The recovered model parameters will be a function of the noise i n the data, so error estimates of the recovered parameters are also calculated. Chapter 4 examines how magnetic permeability and object shape affect the T D E M response. These relationships are formed by observing forward modelled responses of spheres of different magnetic permeability and size, and by observing the measured time decay curves of aluminum and steel rectangular prisms of various aspect ratio. Each of the synthetic modelled and laboratory measured curves are then inverted for their decay parameters. Since these parameters encapsulate information about the physical attributes of the target, we can attempt to use them to determine i f the target is ferrous and i f the geometry is rod-like (most likely a Chapter 1. Introduction 7 U X O ) or plate-like (most likely a non-ordnance item). A n algorithm that might discriminate between U X O and non-UXO targets by considering values of the decay parameters is proposed. In Chapter 5 the discrimination algorithm is applied to an EM61-3D field data set. This data set was acquired over a 105-mm shell buried at the York University Geophysical Test Site. In order to invert this data i n a responsible manner, the nature of the noise i n the field data set is examined to generate a suitable noise model for the inversion algorithm. Chapter 6 concludes this thesis with a summary of the work performed and suggestions for further research. Chapter 2 Development of an Approximate Forward Model In this chapter the response of a metallic object in a T D E M survey is considered. In the T D E M method a transmitter current loop produces a primary magnetic field that illuminates a conductive target (Figure 2.1(a)). If the target is magnetic, as is the case with steel, the primary field will magnetize the target. When the primary field is terminated, the changing magnetic field induces eddy currents in the target according to Faraday's law. The finite conductivity of the target will cause these currents to decay, which then induces additional eddy currents. The rate of decay, and the spatial behaviour of the secondary field, is determined by the target's conductivity, magnetic permeability, shape, and size. These currents produce a secondary field that is measured by receiver loops at the surface (Figure 2.1(b)). Data are measured over several time channels to record the decay characteristics of the buried object. During a T D E M survey these measurements are repeated at several different locations on the surface to record the spatial behaviour of the secondary field. In order to invert measured T D E M data for the physical parameters of the target, it is necessary to have a forward model to describe the T D E M response for a buried metallic object. The T D E M response can be obtained by solving Maxwell's equations. Analytic solutions of Maxwell's equations for conductive and permeable confined bodies are limited to a sphere, and an electromagnetic numerical modelling code for numerical solutions of Maxwell's equations was not available to me. In this chapter, an approximate forward solution is constructed for a buried metallic target. I first present the relevant equations that govern electromagnetic induction, and solve for the T D E M response of a sphere in a uniform step-off magnetic field. 8 Chapter 2. Development of an Approximate Forward Model 9 F i g u r e 2.1: The induction of currents in a buried confined conductor, (a) Transmitter on-time. A transmitter, consisting of a square current-carrying coil, produces a primary magnetic field that illuminates the buried target, (b) Transmitter off-time. Once the transmitter current is terminated, the changing magnetic field induces currents in the buried conductive object. The induced currents produce a secondary field that is measured by receiver coils at the surface. T h i s s o l u t i o n provides insight regarding the spatial nature a n d t i m e decay behaviour of the secondary field. W i t h a d d i t i o n a l p h y s i c a l assumptions, it also leads to m y a p p r o x i m a t e forward model. 2.1 The Time Domain Electromagnetic Response of a Compact Metallic Object 2.1.1 Maxwell's Equations T h e physics of the t i m e d o m a i n electromagnetic m e t h o d is completely described by M a x w e l l ' s equations, d V x E = (2.1) B -5tT V x H = J , 5 D d (2.2) t = p (2.3) V -B = 0 (2.4) V D Chapter 2. Development of an Approximate Forward Model 10 The electromagnetic field is described by four field vectors: E, B, D, and H. The vector J is the electric current density and p is the free electric charge density. E is the electric field and B is the magnetic flux density. A stationary electric charge placed in an electric field will feel a force parallel to E and proportional to the strength of E. The direction of magnetic flux density B is defined to be the direction along which a dipole will align itself when placed in B. D is the electric displacement and its time derivative, dTt/dt is labelled the displacement current density. I will call H the magnetic field. The atomic circulating currents that are responsible for magnetization and the currents from the motion of bound charges in dielectrics are accounted for in the vector fields H and D, respectively. Consequently, J only describes the motion of free charge, such as eddy ctirrents. The electric displacement and magnetic field are introduced solely as a matter of convenience when considering polarizable and magnetizable materials, and Maxwell's equations are often expressed without H and D (for examples see Feynman (1964), Weaver (1994), and Lorrain and Corson (1970)). 2.1.2 Constitutive Relationships Constitutive relationships Unking H and D to B and E are required to solve the electromagnetic field equations. The electric displacement D is defined as D = e E+ P (2.5) 0 where e = (p. c ) 2 0 0 1 • In the absence of dielectric material, the polarization P = 0. For the following work, it will be assumed that all materials are non-polarizable, and thus the electric field and the electric displacement are directly proportional D = eE 0 (2.6) 11 Chapter 2. Development of an Approximate Forward Model The magnetic field H is defined to be related to B through the magnetization M (2.7) H = —B-M(H) Po In this equation M is explicitly written as a function of H . The magnetization vector M is defined as the average magnetic moment per unit volume in the material or, equivalently, the magnetic dipole polarization per unit volume. It is thus convenient to imagine the magnetization of a material as being from an assembly of magnetic dipoles. If these dipoles are distributed evenly throughout the material, the material is said to be miiformly magnetized. In real magnetic materials, a description of the magnetization process is, of course, more complicated. For a non-magnetic material, such as copper, there is no magnetization ( M = 0) and thus the magnetic flux density and the magnetic field are related simply by (2.8) The functional relationship of the magnetization with the magnetic field, M (H), helps classify the three main classes of magnetic materials: diamagnetic, paramagnetic, and ferromagnetic. The constitutive laws for the different magnetic materials are shown in figure 2.2. In diamagnetic a) M b) X<0 Diamagnetic M c) X>0 Paramagnetic Ferromagnetic Figure 2.2: The relationship between the magnetic field H and the magnetization M define the different types of magnetic material. Diamagnetic and paramagnetic materials are characterized by a linear relationship between H and M (Plots (a) and (b)). Plot (c) shows the non-single valued relation between H and M for ferromagnetic materials. and paramagnetic materials the magnetization M and the magnetic field H are linearly related M = xH (2.9) 12 Chapter 2. Development of an Approximate Forward Model where x 1 S the magnetic susceptibility (Figure 2.2(a) and (b)). The magnetic susceptibility of a diamagnetic material is negative, and has an order of magnitude of 10~ . When a diamagnetic 6 substance is placed in a magnetic field, the orbital electrons move in such a way to produce a moment that opposes the applied field (Feynmann, 1964). This is essentially Lenz's law at the atomic level. A l l materials are to some extent diamagnetic. Paramagnetic materials have positive susceptibilities between 10~ and 10~ . Aluminum 3 5 is an example of a paramagnetic metal. Paramagnetic substances are characterized by having permanent magnetic moments of atomic origin in the material. The potential energy of each magnetic moment in an applied field H a is given by — p m 0 • H . When H a a is applied to a paramagnetic material the moments align themselves with the applied field such the potential energy is minimized. The moments do not interact with each other, and when the applied field is terminated, the moments become randomly oriented due to thermal agitation. If the applied field is large enough, the material becomes saturated, and the linear relationship no longer applies. Metals derived from iron or steel are ferrous metals. In ferromagnetic materials it is more energetically favourable for the permanent magnetic moments throughout the material to be aligned (Feynmann, 1964). This is in contrast to a dipole interaction; it is not energetically favourable to have dipoles moments aligned. Throughout a ferromagnetic material there exist domains i n which all moments are aligned. The orientation of magnetization varies from domain to domain. When an external field is applied to the material, domain walls move such that regions that are magnetized opposite to the field are reduced in size. Unlike paramagnetic and diamagnetic materials, in ferromagnetic materials M may not be a single valued function of H. The hysteresis illustrated in Figure 2.2(c) is one of the defining characteristics of ferromagnetic materials. A consequence of hysteresis is that a target will remain magnetized even when the primary field is removed. A simplification in the mathematics Chapter 2. Development of an Approximate Forward Model 13 of the problem is obtained by assuming that the medium is isotropic and hysteresis is unimportant. This simplification is suitable when the H-field is small enough. Referring to figure 2.3, the region O A is reversible, and the linear relationship is very suitable, ff the H-field is small Figure 2.3: Hysteresis curve for a ferromagnetic material. The region OA is reversible. enough then the magnetization characteristic will remain along the O A region of the hysteresis curve. In this portion of the curve the magnetization will go to zero when the applied field is removed. For the analysis in this thesis, a linear relationship between M and H of equation 2.9 is assumed to also apply to ferromagnetic materials. Substitution of equation 2.9 into equation 2.7 gives B = M H + H) = M1 + ) H X X (2.10) This leads to the definition of the magnetic permeability /z = (1 + x)Mo = MrMo (2-11) where fi is the relative permeability of the material. The assumptions leading to the linear r relationship between B and H are considered valid in the analysis of eddy currents i n different types of steel (Stoll, 1974). Indeed, later in this chapter we will see that the measured time domain response of a steel cube is well modelled by an equivalent volume sphere that possesses Chapter 2. Development of an Approximate Forward Model 14 a linear relationship between B and H . Non-ferrous materials, such as aluminum, lead, brass, tin and copper, have unit relative permeability (fi = 1). It has been suggested that a suitable r value of relative permeability for steel targets is /x = 150 (Das et al., 1990). r 2.1.3 Boundary Conditions Maxwell's equations i n differential form describe the field at points where the divergence and curl of the fields exist. These requirements thus exclude surfaces where a and / i are discontinuous. Conditions on the field vectors at the surface are derived by using the integral form of Maxwell's equations. The derivation of boundary conditions can be found i n numerous texts (Jackson (1975), Ward and Hohmann (1991)), and only the final results will be repeated here. The normal component of the electrical field is a discontinuous function at an interface that separates regions of different conductivity and the normal component of magnetic flux density is a continuous at an interface that separates regions of different magnetic permeability n-(D -D )=^, 1 3 n-(B -B ) =0 1 2 (2.12) where p is surface charge density. The tangential components of electric and magnetic fields s are continuous functions n x ( E i - E ) = 0, 2 2.1.4 nx(Hi-H ) = 0 2 (2.13) Quasistatic A s s u m p t i o n A l l analysis i n this thesis will adopt the quasistatic assumption. That is, when using Maxwell's equations, the displacement current term dD/dt i n equation 2.2 will be omitted. To establish the validity of this assumption we can follow the steps i n Weaver (1994). The curl of equation 2.1 can be combined with equation 2.2 and Ohm's law J = crE to give a dimensionless equation i n E alone V ' x ( V ' x E ) + nidE/dt' + Tj d E/dt' 2 2 2 = 0 (2.14) 15 Chapter 2. Development of an Approximate Forward Model where primed variables indicate dimensionless variables, and 771 = pcrL jT and 772 = u-eL /T 2 2 2 are dimensionless coefficients. L is a characteristic length, and T is a characteristic time. The ratio of the coefficients is 771/7/2 = (rT/e indicates the relative importance of the second and third terms. The third term is the displacement current term. The first term of the equation is of order unity. B y choosing T to be the earliest time measurement for E M 6 1 (T = 0.01msec), letting e = e — 8.84 x 1 0 Q - 1 2 , and recognizing that the conductivity of metals is of the order 10 S / m , the value of the ratio is 771/772 = 1.13 x 1 0 . That is, the third term is not important 6 12 relative to the second term. Therefore the quasistatic assumption is valid for this problem. 2.2 Electromagnetic Induction as a Diffusive Process Figure 2.4 illustrates the physical situation we wish to consider. A n arbitrary shaped conductive (and possibly permeable) body is illuminated by a primary field B p (Figure 2.4(a)). A t t — 0, the primary field is terminated, and currents are induced i n the conductive target that give rise to a secondary B-field which is subsequently measured. B b ) x P i a = 0, u a, u. 0 Figure 2.4: (a) An arbitrary conductive body immersed in a primary field B . (b) The step-off primary p field. Prior to the termination of the primary field, the distribution of the B-field both exterior and interior to the metallic body can be determined by solving the appropriate magnetostatic equations. If the body is non-permeable, a uniform B - f i e l d will remain uniform everywhere. If p the body is permeable, the primary field will magnetize the body. The B-field inside the body will increase, and the field lines appear to be 'sucked' into the body (Figure 2.5). Chapter 2. Development of an Approximate Forward Model 16 A s a n example, let us consider the magnetostatic s o l u t i o n of a sphere. T h e whole space solution for the interior B-field of a sphere i l l u m i n a t e d by a constant, u n i f o r m field B° is 2+ ^ where u = xz/\i is the relative permeability. For a non-permeable target {fi — 1) the interior r 0 T field, B% w i l l be equal to the applied field. A s the p e r m e a b i l i t y increases, the interior field also increases, a n d asymptotes to B* = 3B Q for very permeable materials (/z > > 1). Therefore, the r interior B-field of a permeable sphere is greater t h a n the interior B-field of a non-permeable sphere. Subsequently, once the p r i m a r y field is terminated, the currents i n d u c e d on the surface of the permeable sphere w i l l be greater t h a n on the surface o f a non-permeable sphere. The field (and response) p r o d u c e d by currents on the permeable sphere w i l l be greater t h a n the field produced by a non-permeable sphere. F i g u r e 2.5: The effect of permeability in a uniform field B ° . The density of lines is proportional to the magnitude of the field B-field vector, and the lines are tangent to the direction of the field vector at each point. On the left, the permeability fi > /x and thus the interior B-field of the sphere, B is greater than B ° . O n the right, the sphere is non-permeable and thus there is no effect on B ° . 1 0 T h e extent to w h i c h the b o d y is magnetized (and the interior field is increased) is dependent on how m a g n e t i c a l l y permeable the target is, on the shape of the target, and o n the o r i e n t a t i o n of the target relative to the p r i m a r y field. T h i s effect, k n o w n as d e m a g n e t i z a t i o n , is an i m p o r t a n t factor i n the response of steel targets and w i l l be discussed i n C h a p t e r F o u r i n some detail. Chapter 2. Development of an Approximate Forward Model 17 At some time the primary field is terminated (at t = 0 in Figure 2.4(b)). The changing magnetic field then induces eddy currents on the surface of the conductor. The distribution of induced currents and B-field can then be determined by Maxwell's equations. B y assuming quasistatic fields and homogeneous media, Maxwell's equations can be rewritten to show that the electric field, magnetic flux density, and the induced currents behave according to the diffusion equation *v. = £ , * V B = f, *vj = | (2.i«) where the diffusivity k is k=— oyz (2.17) The process of electromagnetic induction can be described as a diffusive phenomenon. The time decay response is generally divided into three separate stages (Kaufman, 1985; McNeill and Bosnar, 1996; Nabighian and Macnae, 1991). The early time stage occurs immediately following the termination of the primary magnetic field. Once the primary field is terminated, currents are induced on the surface of the target according to Faraday's Law. The B-field outside the target will experience a step change upon termination of the primary field. The B-field inside the conductor, however, will be unchanged upon termination of the primary field. In accordance with Lenz's Law, currents are induced on the surface of the conductive body to maintain the B-field interior to the sphere prior to the termination of the primary field. A larger interior B-field will give rise to greater magnitude eddy currents being induced on the target surface. Once currents are induced on the target surface, the subsequent behaviour of the currents obeys the diffusion equation. Quantities that obey the diffusion (or heat) equation tend to move, or diffuse, from areas of high concentration to areas of low concentration. The induced surface currents therefore diffuse from the surface (a region of high current concentration), towards Chapter 2. Development of an Approximate Forward Model Early Time Intermediate Time 18 Late Time Figure 2.6: The three time decay stages of the T D E M response of a sphere. At the early time stage induced currents are restricted to the surface. At intermediate times currents diffuse towards the center of the target. At late times the currents achieve a final spatial distribution, at which point all currents decay exponentially. the center of the conductive body (a region of low current concentration). The diffusion of currents towards the center of the sphere characterizes the intermediate time stage. The physical picture used to describe the diffusion of the eddy currents towards the target center is based on Faraday's Law. Because the target has finite conductivity, the surface currents will decay, and the associated changes of B will induce more eddy currents towards the center of the target. This process continues until the spatial distribution of currents reaches a steady state. A t this late time stage the equilibrium distribution of currents simultaneously decay exponentially. The diffusivity k characterizes the diffusion rate of the fields and currents. The larger the diffusivity constant, the more diffusive the medium and thus the late time stage's steady state distribution will be reached sooner. Indeed, if the target were infinitely conductive then the material would be completely non-diffusive. Physically, currents would not decay and induce more currents towards the center of the target. In the case of an infinitely conductive target the currents would remain on the surface and not penetrate into the interior of the target. Two materials of interest are aluminum and steel. Aluminum has a smaller conductivity-permeability product (f/z), and therefore has a larger diffusivity than steel. The onset of the late time stage would occur sooner in an aluminum target than a steel target of the same shape and size. 19 Chapter 2. Development of an Approximate Forward Model 2.3 Approximating the Response of an Arbitrary Shape Target Kaufman (1994) derived a general form for the field caused by currents induced in a confined conductor. Consider an arbitrary shaped conductive body i n free space illuminated by a stepoff primary magnetic field (Figure 2.4). B y assuming quasistatic fields, the secondary field produced by currents i n the confined conductor can be written as oo B (t, )= l (BP P • i) £ d i n W~ e (- ) t/Tn 2 18 n=l where (t,p) is the secondary field in the 1 direction at a time t following the termination of a primary field. The field is observed at a point p. The value B P • 1 is the projection of the primary field along the 1 direction, d^ are coefficients that depend on the target location, size, and shape, and upon the geometry of the primary field. r n is also dependent on the size and shape of the target, but not the target location and geometry of the primary field. Each component has, in general, its own set of d ^ components. T\ is the largest value of r n n and is referred to as the diffusion time constant of the conductor. The late time decay of the field will be exponential as e a Tl = ^ a _ t / , T l . For a magnetic sphere this constant is ,2 (2.19) where a is the radius of the sphere. The value of T\ has been calculated for a number of different non-permeable targets i n a uniform field (Kaufman, 1994), and is used to identify the quality of conductors in mining exploration T D E M surveys (Nabighian and Macnae, 1991) and has been proposed as a diagnostic in discrimination and identifying U X O (McNeill and Bosnar). From equation 2.18 we see that the transient behaviour of a confined conductor can be expressed as an infinite sum of exponential functions. For a finite time window the response can be represented by only a finite number of terms. Attempts have been made to invert decay curves for a sum of exponentials. Such attempts have seen little success due to coupling among the coefficients (McFee and Das, 1995). 20 Chapter 2. Development of an Approximate Forward Model Rather than fitting decay curves to a sum of exponentials, I chose to develop an approximate forward model based on the analytic solution of a conductive and magnetically permeable sphere. The sphere solution was used in two ways. Firstly, I used the basic structure of the sphere solution as a blueprint for the approximate model. The solution of a sphere can be written as the product of two functions. The first function is dependent on the geometry of the primary field and observation location, and describes the secondary field as a dipolar field. The second function is dependent only on the physical characteristics of the target, such as size, conductivity, permeability, and shape. This second function determines the time decay characteristics of the measured response. The form of the approximate forward model retains this structure. A second way in which the solution for a sphere was utilized, was to provide insight on how the shape of decay curves vary with conductivity, permeability, and size. W i t h this information, a time decay law could be constructed such that responses predicted by the forward modelling can (at the least) reproduce the different features of the response decay of a sphere. 2.4 The Time Domain Response of a Permeable, Conductive Sphere in a Uniform Field Consider a sphere with a radius a, conductivity cr, and a magnetic permeability /j, i n an insulating whole space. When this sphere is illuminated by a uniform (primary) field B p that is terminated at t — 0, currents are induced on the sphere that give rise to a secondary field. The solution of the T D E M response is outlined here, with only the major points being mentioned. For all the details of the solution, see Kaufman's Inductive Mining Prospecting (1985). We use the quasistatic form of Maxwell's equations (2.20) V x H = <TE (2.21) Chapter 2. Development of an Approximate Forward 21 Model o- =0, u o 0 Figure 2.7: Set-up for a sphere in a uniform field B ° . V-E = 0 (2.22) V-B = 0 (2.23) The divergence of the electric field E is zero since the symmetry of the problem implies that the E-field does not intersect the sphere surface and thus no charges accumulate. B y introducing a vector potential A such that B = V x A (2.24) and by adopting a Gauge condition V • A = -au-V (2.25) we obtain a diffusion equation in A (9 A V A = <rp°— (2.26) 2 Once the diffusion equation is solved, the electric and magnetic fields can be obtained from B = V X A, dA dt E = —— v (2.27) ; The first step is to determine a solution to the diffusion equation. The diffusion equation is separable i n spherical co-ordinates. Once a suitable function is obtained for A , the boundary conditions of equations 2.13 requires that the tangential components of the field are continuous 22 Chapter 2. Development of an Approximate Forward Model at the sphere surface. Finally, the field interior to the sphere at the instant the primary field is terminated is given by the magnetostatic response of the sphere B int J ^ r _ = B P ( 2 + fir 2 2 g ) ' V Enforcing this as the initial condition allows us to determine the remaining unknown coefficients of A . Equations 2.27 can be used to determine the expression for the B-field from A . The expressions for the secondary field B are (Kaufman, 1985) s 3 B s 3 = B?^cos6L (t), B$ = B?~smeL (t), B R B$ = 0. B (2.29) The time decay behaviour of this secondary field is governed by L (t) B oo —ga* i (*) = 6/*rg B where r = ap,a , and 2 f t+ ( / i r fi r = e _ T pjfi 0 (2.30) 1 ) ( A | r +2 ) is the relative permeability. The values q are roots to the 3 transcendental equation t a n « « = ~2~r? TT ( 2 3 1 Is + (Mr - 1) This equation is solved by bracketing the interval where the root occurs, then using a bisection method to find the location of the zero i n the interval. The superscript B i n equation 2.30 indicates that this expression is part of the solution for the B-field response. The time derivative of the flux density 83/dt is of interest, since it is the quantity sensed by a receiver coil. Clearly, to obtain the time derivative of the field requires simply taking the derivative of L i n B equation 2.29. The time decay law for the derivative of the field is OO (t) = - 6 ^ Y —1st (2.32) In calculating the 83/dt and B-field responses, the infinite sum in equations 2.30 and 2.32, respectively, must be evaluated. Because each term of the series includes a solution to the Chapter 2. Development of an Approximate Forward Model 23 transcendental equation 2.31, it is difficult to perform standard convergence tests on either series. Therefore, to check convergence I am left with graphical techniques. The graphical analysis is also used to check how many terms in the sum must be retained in order to produce accurate results. Figure 2.8 shows how different terms i n the summation contribute to the calculation of dB/dt or the B-field of a steel sphere. In panel (a), individual terms of the summation used to calculate the B-field are plotted. Terms corresponding to s = 1, 5, 10, 50, 100, 500, 1000, and 2000 are plotted. Only the first few terms in the summation make a significant contribution to the response at later times. As s increases, the size of the corresponding terms decrease and only make a significant contribxition at earlier times. The contribution of the s = 500, 1000, and 2000 terms are not large enough to be included i n the plot. Figure 2.8(b) contains the B-field calculated by truncating the summation at different numbers of terms. The summation can be terminated earlier when the time window of interest is later in time. Indeed, the difference i n carrying 1000 or 2000 terms in the summation make a very small difference. Earlier times require more terms i n the summation to be calculation. It is clear that a finite number of terms of the summation is sufficient for the secondary response (either dB/dt or B-field) to be accurately calculated. Because the time channels at which EM61-3D makes measurements is between 0.01 and 100 milliseconds, my choice of carrying 5000 terms in the summation when calculating the response of a sphere, is more than adequate. This conclusion was reached after testing on both permeable and non-permeable spheres of different diameters. The calculated dB/dt and B-field responses for a steel and aluminum sphere are plotted in figure 2.9 The relative permeability of fi = 150 was chosen for modelling the steel sphere i nfigure2.9. r Because steel is ferromagnetic, a linear relationship B = / z H is generally not applicable. To ensure that the linear assumption, as well as the choice of \i is appropriate, the measured time T decay response of a cube and the analytic response of a sphere of equivalent volume and a 24 Chapter 2. Development of an Approximate Forward Model -q.t/1 e q + ( k-1)( k+2 ) (a) 10 10 10 0.01 100 1 1 x 10 1 0 0.01 Time (msec) 1 100 Time (msec) -qt/x 10 s = 2000^ f ! - 5 0 0 ^ s = 1000 S _ 0 U U s = 1 0 ° 10 0.01 1 100 Time (msec) F i g u r e 2.8: Examining the number of the terms that should be kept in the summation performed in calculating the B-field and dB/dt response of a 20 cm diameter steel sphere, (a) Individual terms of the summation used to calculate the B-field of a steel sphere are plotted for s = 1, 5, 10, 50, 100, 500, 1000, 2000. Terms corresponding to s = 500, 1000, and 2000 make a negligible contribution during the time window plotted above, (b) The B-field calculated by truncating the summation at a different number of terms. For earlier times, more terms of the summation must be retained, (c) Individual terms of the summation used to calculate dB/dt are plotted. Terms with large s only make a contribution for earlier times, (d) The dB/dt response calculated by truncating the summation at a different number of terms. For earlier times, more terms of the summation must be retained. When calculating the response for earlier times, more terms i n the summation must be included. 25 Chapter 2. Development of an Approximate Forward Model Time (msec) Time (msec) Figure 2.9: (a) The time decay behaviour of the magnetic flux density B . The B-field response is normalized by the strength of the primary field, (b) The time decay behaviour of the time derivative of the magnetic field dB/dt. relative permeability of \i = 150 are compared (Figure 2.10). The forward modelled steel and v aluminum sphere solutions match well with the curves of the measured responses for steel and aluminum cubes of equal volume, implying that the linear relationship between magnetization and H-field may be adequate in describing the steel targets. The measured response curves were provided by Geonics L t d . and measured i n a laboratory setting. Full details of the acquisition of these data are found in Chapter 4. In the preceding analysis the primary B-field was be assumed to be uniform i n the vicinity of the sphere. This assumption will be valid when the sphere is more than a radial length from an EM61 transmitter coil (Barrow, et.al. 1996). A major result for a sphere i n a uniform field is that the secondary field will be a dipolar field that decays as a function of size, shape, conductivity, and permeability. A non-uniform field will produce higher order moments that decay quickly with time (Grant and West, 1965). A second assumption, used above, is that the metallic sphere is sitting in free space. This greatly simplifies analysis and is justified since the conductivity of metals will be approximately 6 orders of magnitude greater than the ground (Das et al., 1990). Chapter 2. Development of an Approximate Forward Model 26 Steel Sphere (u = 150 n ) 0 10 -=i 1.0 10 10 Time (msec) 10 10 Figure 2.10: The responses of a metallic cube compared to a sphere of equal volume. The dB/dt responses of the cubes were normalized by the measured value of the steel cube at t = 0.01msec. The responses of the spheres were normalized by the calculated value for a steel sphere at t = 0.01msec. Normalization was required to counter any extra gain in the signal during measurement. 2.5 The Decay Characteristics of the T D E M Response of a Sphere Section 2.2 described the three stages of the time domain response i n terms of the formation and diffusion of induced eddy currents throughout a conductive target. This section describes in what manner these processes manifest themselves in the time domain response of a sphere. Each stage is characterized by a different behaviour of the currents. Figure 2.9 illustrates the time decay curve of both a non-permeable and permeable sphere. Let us consider the features of the dB/dt and B-field response separately. 2.5.1 Decay Characteristics of B-field Response Figure 2.9(a) contains the B-field response of both a steel and aluminum sphere. In the response of the steel sphere, there are three distinct sections of the B-field response. The relatively flat portion of the B-field response at the beginning of the response identifies the early time stage. The length of this first section is proportional to the conductivity, permeability, and radius of Chapter 2. Development of an Approximate Forward Model 27 the sphere. The intermediate time stage for the steel sphere is characterized by a straight line i n a log(B) versus log(i) plot. The response during the intermediate stage exhibits a t~@ decay behaviour, where /3 is approximately 0.5. This intermediate time behaviour is not seen i n the response of a non-permeable sphere, so this decay behaviour is clearly due to magnetic effects. This behaviour has been explained by treating the sphere as a perfect paramagnetic (Zhadnov, 1997; McNeill, 1997). That is, magnetic dipoles become aligned by the primary field and then relax to random orientations due to thermal agitation. At late times the B-field decays exponentially. This is not surprising since the T D E M solution is written as the sum of exponential functions exp{—t/r). The largest value of T, known as the diffusion time constant of the sphere, is r = a pa , where a is the sphere conductivity, p 2 i s the magnetic permeability, and a is the radius. The value of r determines the onset of the exponential, and final, time stage. The response of an aluminum sphere does not exhibit the three distinct sections observed for the steel sphere. The aluminum sphere has only a flat early time section, and an exponential, late time section. The intermediate time stage, and its i - 1 / 2 decay behaviour, is not evident. The onset of exponential behaviour occurs earlier due to a smaller diffusion time constant r . The early time B-field response of the aluminum sphere is weaker than the response of the steel sphere, as expected (see section 2.2). 2.5.2 D e c a y Characteristics o f 83/dt Response Figure 2.9(b) contains the 83/dt response of both a steel and aluminum sphere. As was seen in the B-field response of a steel sphere, the 83/dt response of the steel sphere has three distinct sections. In contrast to the B-field response, the 83/dt response early time behaviour is not flat. In the early time stage there is a t~$ decay behaviour with /? approximately 0.5. The intermediate time stage that follows decays again as t~@, but with /3 ss 1. The final decay stage 28 Chapter 2. Development of an Approximate Forward Model is again exponential. The dB/dt decay of the aluminum sphere only has two distinct sections. A t early time the response decays as r / . Kaufman (1994, pg. 229) shows that for early time the dB/dt - 1 2 response for a non-permeable sphere is „<; BR = 3Bf / a? 7^ c o s e 3B 6 = -T-F^ • , n a? P c B -7^ s i n „<j ^. B? = 0. , (2.33) At late times the decay is exponential. 2.6 Statement of F o r w a r d M o d e l The basic results of a permeable and conducting sphere provide a blueprint for an approximate forward model for the response of a metallic target. The analytic expression for both the secondary B-field and the time derivative dB/dt of the secondary field of a sphere consists of the product of two functions (also see McNeill & Bosnar 1996). Recall that the secondary B-field of a sphere is B S = B ^ cos9L*(t), B$ = B?^- smeL»(t), p R z B* = 0. 3 (2.34) These equations reveal that the B-field of a sphere in a uniform primary field is equivalent to the B-field of a single magnetic dipole located at the center of the sphere and oriented parallel to the primary field. The secondary B-field can be rewritten as B s = M(r)L(t) (2.35) where M (r) is the B-field due to a magnetic dipole m = 2 7 r B / / z . The decay of the dipole p 0 field is governed by L(t) = a L (t). 3 B L (t) is defined i n equation 2.30. Equation 2.35 also B shows that the time derivative of the B-field will have the same spatial behaviour (governed by M(r)) as the B-field. The first term M(r) is a function of the magnitude and geometry of the primary field, and includes no information on the material and shape properties of the target. This shows that Chapter 2. Development of an Approximate Forward Model 29 the spatial response of the sphere is identical to a dipole whose strength is proportional to the primary field. The second function L (t) contains all the information on the time decay of the sphere and it depends upon the material properties, shape, and size of the target. M y hypothesis is that more general metallic shapes can also be approximately modelled as the product of M(r) and L (t). However, choosing the right functional form of M(r) and L (t) will be crucial. 2.6.1 Spatial Response Function M(r) For an Axi-symmetric Body Most targets encountered i n a U X O survey can be adequately described as an axi-symmetric object. Unfortunately, analytic expressions for the time domain response of a permeable and conducting non-spherical axi-symmetric body are not available. However, recalling that for the time domain response of a sphere M(r) was identical to the response of a magnetostatic sphere (Equation 2.35, it is reasonable to consider the magnetostatic response for an axi-symmetric target as a possible candidate for M(r). Analytic solutions for the magnetostatic response of a magnetic prolate spheroid and the electromagnetic response of a perfectly conducting spheroid are available. As with the sphere, the spatial response is equivalent to a magnetic dipole, with moment m , 3ph induced at the spheroid center. This dipole moment is best understood as being the sum of two orthogonal dipole moments. Let us resolve the primary field into two orthogonal components, where one component is along the major axis of the spheroid. The first dipole moment is parallel to the main axis of the spheroid and its strength is proportional to the component primary field along that direction. The second dipole is perpendicular to the major axis, and its strength is proportional to the component of the primary field along that direction. The total induced dipole moment is (Das et al., 1990) m sph =k (all • # ) all + k p L (H T T - (all • #) p a") (2.36) Chapter 2. Development of an Approximate Forward Model 30 where a is the unit vector parallel to the main axis of the spheroid, k^ and kx, which premultiply the dipoles along the major and transverse axes respectively, are functions of the conductivity, permeability, shape, and size. The general expressions for and kx are well known for both the magnetic spheroid and perfectly conducting spheroid. It has been reported that this two dipole model represents the spatial response of U X O very well (Khadr et al., 1998). 2.6.2 T i m e Decay Function L(t) Recall that the time decay for a sphere is determined by the sum of exponentials. This result generalizes to the case of a conductive body of arbitrary size and shape in an insulating medium (equation 2.29). Thus we seek a form for L (t) that duplicates the time decay features as was observed for the sphere. More specifically, L (t) should be chosen such that early, intermediate, and late time stages can be realized. A n appropriate form of the decay law for the B-field response is * (t + a ) " e ^ . / 3 (2.37) The parameter k controls the magnitude of the modelled response. The three parameters a , (3, and 7 , control the duration and characteristics of the three different stages of the time decay curve. The duration of the relatively flat early time stage will be proportional to the parameter a. The linear decrease of response observed during the intermediate time stage is determined by t~&. The exponential decay characterizing the late time stage is controlled by the parameter 7 . The relationships between the decay parameters and the decay characteristics are illustrated in Figure 2.11. The 83/dt response of a sphere is plotted in Figure 2.12(b). In section 2.5.2, it was noted that the decay curve for a steel sphere coidd be divided into three stages. The decay nature of each stage differs slightly from the decay calculated for the B-field. The intermediate stage has a greater slope for 83/dt than B-field. The early time decay behaviour of 83/dt does 31 Chapter 2. Development of an Approximate Forward Model c) Changing p Changing y 10 Time (msec) F i g u r e 2.11: Decay behavior changes due to variations in parameters k, a, (3, and 7. (a) k shifts the curve up and down, (b) The length of the initial, flat stage of the response increases with increasing a. (c) The steepness of the intermediate stage of the response increases with increasing /?. (d) The onset of the exponential, final time stage is delayed with increasing 7. Chapter 2. Development of an Approximate Forward Model 32 Figure 2.12: (a) The time decay behaviour of the magnetic flux density B. The B-field response is normalized by the strength of the primary field, (b) The time decay behaviour of the time derivative of the magnetic field dB/dt. not possess the fiat section as is seen in the B-field response. However, during a typical time measurement window only the transition from the early to the intermediate time decay is seen. The time decay law 2.37 is a again a suitable description of the field behavior in such a case (Figure 2.12(b)). The termination of the early time stage will again be controlled by a. The difference in units between the dB/dt and B-field responses will be accounted for in the units oik. In the previous section I proposed that the spatial response of an axi-symmetric target be represented by a pair of dipoles located at the center of the target. It has been noted that the shape anomaly of the measured response changes with time (Grimm et al., 1997). The physical phenomena that gave rise to the temporal changes in shape anomaly was explained i n terms of the nature of the induced eddy currents. Eddy currents that circulate end-to-end in the U X O dominate at early time but decay away quickly, and eddy currents that circulate about the long axis extend later into time. This field behaviour can be duplicated by letting each of the two orthogonal dipoles of the previous section decay independently of each other. Let us assume that the currents circulating end-to-end and the currents circulating about the long axis generate secondary fields that can Chapter 2. Development of an Approximate Forward Model 33 Figure 2.13: Representation of the fields generated by induced fields in a metallic target as a pair of orthogonal dipoles. Both dipoles are located at the center of the target, with Dipole 1 being parallel to the axis of symmetry of the target and Dipole 2 being perpendicular to the axis of symmetry and in a plane defined by the primary field and Dipole 1. The moment of each dipole is proportional to the projection of the primary field onto a unit vector in the dipole direction. be approximated by the fields of a dipole perpendicular to the long axis and a dipole parallel to the long axis, respectively. B y assigning a different decay characteristic (governed by its decay parameters) to each dipole, the relative contribution by each dipole to the secondary field can vary with time. Indeed, the two sets of decay parameters can be chosen such that the dipole perpendicular to the symmetry axis (corresponding to currents circulating end-to-end) will be dominant at early times and the dipole parallel to the symmetry axis (corresponding to currents circulating about the long axis) dominant at late time. 2.6.3 The Approximate Forward Model Finally, we can write an approximate expression for the secondary field response of an axisymmetric target. The approximate forward model represents a buried metallic object as two orthogonal dipoles located at the center of the object (Figtire 2.13). Each dipole decays according to the decay law of equation 2.37. Let us consider a target whose center is located at R. The response measured at a receiver/transmitter location r and at a time t after the termination of the primary field, is then the sum of the responses of the two orthogonal dipoles £ (r, t) = B (r) L (t) + B (r) L (t), 1 x 2 2 (2.38) 34 Chapter 2. Development of an Approximate Forward Model Figure 2.14: Co-ordinate System for a buried target The response £(r,t) can represent either the B-field or the time derivative dB/dt. L is the x decay law of the i dipole moment th L (t) i = k (t + a y< 'e^ (2.39) 3 i i 1 and B j is the B-field due to the i th dipole moment mj (2.40) The induced dipole moment rti! is defined to be the projection of the primary field onto the axis of symmetry m i = (xx B )x (2.41) p x where X x is a unit vector parallel to the axis of symmetry, and B from a square loop source. The equations used to model B p p is the primary field generated are outlined in Appendix C . The unit vector parallel to the axis of symmetry of the target, i j , is simply the direction cosines of the axis of symmetry i i = cos(a) sin(c) x + sin(a) sin(c) y + cos(c) z (2.42) where the orientation angles a and c are defined in figure 2.14. The second dipole moment is chosen to be parallel with the inducing field and orthogonal to the main axis of the target 35 Chapter 2. Development of an Approximate Forward Model denned by X j . The unit vector defining the direction of the second dipole, x , is parallel to the 2 projection of the primary field onto a plane whose normal is X i . Two orthogonal unit vectors x a and Xb that lie in this plane are x = — cos(a) cos(c) x — sin(a) cos(c) y + sin(c) z a (2.43) and i b = sin(a) x — cos(a) y. (2.44) The projection of the primary field onto this plane, and therefore the magnitude of the second induced dipole moment ni2, is then m = (x • B ) x + (x • B ) x . p 2 a (2.45) p a b b In summary, the approximate response of buried metallic object can be generated from 13 parameters that describe the object. These model parameters are elements of the model vector m = [X, Y, Z, a, c, a , 0 , 71, & , a , f3 , 72] • t t 2 2 2 (2.46) The location the object is given by the parameters X, Y, and Z. X and Y is the location on the surface directly above the object, and Z is the depth of the object below the surface. The orientation of the target is described by the two angles a and c. The remaining parameters describe the decay characteristics of the two dipoles: fej, cti, /3i, and 71 describe the dipole parallel to the axis of symmetry ( m i ) , and k , a , (3 , and 72 describe the dipole perpendicular 2 to the axis of symmetry ( m ) . 2 2 2 Thus the inversion for the model m will immediately give estimates of target location and orientation. Information on the shape, size, and material parameters of the target may later be inferred from the remaining parameters. 2.6.4 Approximate Forward Modelling Examples I conclude this chapter with examples of the approximate forward modelling. Let us consider approximate dH/dt response of a steel rod and a steel plate. A rod and plate represent the 36 Chapter 2. Development of an Approximate Forward Model two end-members of the set of axi-symmetric targets, with a rod being geometrically similar to a U X O and a plate similar to scrap metal. These data sets are generated by substituting parameters recovered from lab measurements taken at Geonics L t d . into equation 2.38. The parameters and size of the rod and plate are listed in table 2.1. Further details of the different lab measurements and their subsequent inversion for parameters are covered i n Chapters Three and Four. Dipole 2 Dipole 1 Target Rod Plate Dimension (inch) 8x2x2 8x8x2 h 0.841 0.951 ai 0.0191 0.0204 Pi 1.01 1.30 7i 22.4 22.6 k 2 0.140 2.07 0.00702 0.0356 fa 1.32 1.11 72 2.92 22.2 Table 2.1: Parameters for the dB/dt response of a steel rod and a steel plate Figure 2.15 contains the forward-modelled response of the two targets using parameters listed in table 2.1. Each target is assumed to be located 0.5 m below the surface and at survey location (X, Y) — (2m, 2m). Each target's axis of symmetry with the axis of symmetry being oriented with a dip c = 4 5 ° and a strike of a = 4 5 ° . The top row of Figure 2.15(a) plots the plan view of the vertical component of the dB/dt response of a steel rod for three different time channels. The second row of Figure 2.15(a) plots the plan view of the horizontal component. Figure 2.15(b) contains responses for a steel plate. In both the case of a rod and a plate, there is a change in the anomaly shape in each time channel. The response of a rod changes from being due to a both dipoles at early times, to essentially the dipole parallel to the length of the rod at late times. The late time response of the plate is predominantly from the dipole perpendicular to the target and in the plane of the plate. Chapter 2. Development of an Approximate Forward Model (a) 37 8 x 2 x 2 inch rod t = 0.036 msec t = 0.29 msec x 1 o 3 t = 8.32 msec la Vertical Component Horizontal Component (b) 8 x 8 x 2 inch plate t = 0.036 msec t = 0.29 msec t = 8.32 msec Vertical Component l8B\ 2 Horizontal Component (t t (f i + F i g u r e 2.15: Synthetic dB/dt data for a steel rod and a steel plate using the approximate forward model, (a) The vertical and horizontal components of the secondary field for a steel rod are plotted respectively in the top two rows. The change in the observed anomaly is due to the relative contribution of the two dipoles changing with time, (b) The secondary field for a plate are plotted in the bottom two rows. Once again there is a change in the plotted anomaly pattern between early and late times. Chapter 3 Non-linear Parameter Estimation Procedure A n important step i n achieving the goal of U X O discrimination is to extract physical parameters of a buried metal object from time domain electromagnetic data. A T D E M survey involves generating a pulse primary magnetic field to induce currents within the target, and then measuring the secondary magnetic field that these currents produce. The measurement is made over a number of time channels, and at several locations on the surface above the target of interest. Figures 3.1 and 3.2 contain portions of the EM61-3D dB/dt data set acquired during a T D E M survey over a 105 m m U X O buried approximately 0.5 m below the surface. These data were collected over a 4 m x 4 m square survey area, along 9 lines (with l m line spacing) and with a station spacing of 50 cm along each line. The decaying secondary field was measured over 30 logarithmically spaced time gates. This data set include over 13 000 data. Figure 3.1 show the 3 components of the dB/dt response at t — 0.8025 msec and t = 27.92 msec. Figure 3.2 include both the dB/dt and B-field decay curves of the target measured at station (X,Y) = (lm, lm). The inversion of these observed data requires a forward modelling code to accurately predict the response of electromagnetic anomalies sensed by the EM61-3D detector. The ideal forward model would generate time domain electromagnetic responses for arbitrary shaped metallic objects by numerically modelling Maxwell's equations. Because such a code was unavailable during the time of this thesis, an approximate forward modelling that represents the response of a compact metallic object as two dipoles located at the center of the target was proposed. In 38 Chapter 3. Non-linear Parameter Estimation Procedure 39 t = 0.8025 msec 4ft* t = 27.92 msec -5 -10 lis Figure 3.1: The dB/dt measurement taken over a 105 mm shell. At earlier times the signal has the appearance of a dipole field. With later times the signal from the UXO is dominated by the noise. Chapter 2 either the dB/dt or B-field response, represented by £, of a buried metallic object was approximated by t(r,t) = B (r)L 1 1 2 2 where Li is the decay law of the i th Li(t) = *i(t + a i ) - e ^ , B i , the B-field due to the i th r ri dipole moment: » = 1,2. f t u„ / , B i = f 3 mi • ( 4x \ (3.1) (t) + B (r) L (t) (3.2) dipole moment mi, is - R (ri - R) mi ^ — - £ ri - R | |n - R | (3.3) The induced dipole moment m; is mi = (xi • B ) Xi p (3.4) For this thesis, I first assume that the response measured in a survey is due to a single body, and second, that the response of this single body can be accurately modelled with equation 2.38. Chapter 3. Non-linear Parameter Estimation Procedure 40 Figure 3.2: Decay of the response at [X, Y) — (lm, lm) (a) The early time channels for the dB/dt response is very noisy, and need not be fit too closely, (b) The process of integration smoothes the response for B. W i t h these hypotheses in place, an inversion procedure can be developed that utilizes the approximate forward model. The forward model can be expressed as d^F^xn], j = 1,2,3, ...N (3.5) This equation expresses the mapping of the model vector m to a datum dj by a functional Fj. The forward mapping Fj is denned by equation 2.38. The model vector m contains the 13 model parameters m = [X, Y , Z, a, c, fci, cti, /3i, 71, & , a , /3 , 72] • 2 2 2 (3.6) The goal of this inversion is to retrieve the 13 model parameters making up the model vector m from a vector of observed data d o b s . If there are L time channels and K locations where the EM61-3D data is collected, then there ' will N = 3KL total data points contained in the data vector d. It is clear from equations 3.1 to 3.4 that the N data will be non-linear in each of the M model parameters. Because data will be collected on several lines, with a number of stations per line, there will generally be far more data than model parameters (N >> M). Therefore the inversion for the model parameters of m involves solving an overdetermined system of non-linear equations. Chapter 3. Non-linear Parameter Estimation Procedure 41 For overdetermined inverse problems the goal is to find the model that produces the data that best fits the observed data. This type of problem is a departure from many of the standard geophysical inverse problems, where the subsurface is divided into many cells and the goal is to determine a physical parameter, such as conductivity in each cell. The parameterization of the subsurface generally leads to underdetermined problems, which are ill-posed, and possess a very large inherent non-uniqueness. In such problems a standard method of dealing with illposedness and non-uniqueness is to define a model objective function which forces the inversion to construct models with characteristics consistent with the inverter's notion of an appropriate model. Adopting the forward modelling of equation 2.38 essentially restricts the search of an appropriate model to very specific model type. The extreme non-uniqueness that characterizes most inverse problems, is not the major issue in solving the overdetermined inverse problem with which we are faced. The primary focus of the inversion, therefore, shifts from dealing with non-uniqueness to fitting the data. The basic model parameter estimation procedure is to first define an objective function which gives a measure of the agreement between the observed data and the data predicted by a set of model parameter estimates. Such an objective function is designed to decrease with improved agreement between observed and predicted data. Therefore, parameters are perturbed until a minimum in the objective function is achieved. Obtaining the best fit parameters can thus be expressed as a problem of multivariate optimization. The non-linearity of the forward modelling requires that the inversion procedure is an iterative process. This chapter outlines the implementation of a non-linear least squares algorithm to obtain an estimate of m. W i t h the use of synthetic data sets, this procedure is then evaluated to determine which of the 13 model parameters are most (and least) robust to variations i n data collection configurations and to noise. 42 Chapter 3. Non-linear Parameter Estimation Procedure 3.1 Denning The Objective Function The problem of determining the model parameters of a U X O from a data set can be solved as a non-linear least squares problem. In the method of least squares we try to find the model parameters that minimize the sum of the squared differences of the predicted and observed data. The non-hnear least squares problem is </> = ^\\W (F[UI] - d minimize d o b s ) || (3.7) 2 where F[m] is the forward modelling that produces the predicted data, d b is the observed c s data, and (j) is the least squares objective function that measures how closely our predicted data matches the observed data. Wd is the data weighting matrix. For ease of notation the least squares problem is rewritten as I minimize T (f> (m) — -R (m) 1 R (m) = - ^ N (m) (3-8) 2 i=l where R is the residual function R= W {F[m}-d ) d and T{ (m) is the i obs th component of the function R. The data weighting matrix Wd adjusts the relative contribution of each T{ to the objective function. Therefore Wd enables us to control how closely each datum is fit by the predicted data. There are two motivating factors to including a data weighting matrix. Firstly, the data weighting matrix allows me to deal the with large dynamic range exhibited i n T D E M data. The measurements of the fields i n a T D E M survey will take place from approximately 10~ to 5 0.1 seconds. Over such a large time range, the measured secondary field and induced voltage can vary by over 5 orders of magnitude for the EM-61. If we were to sum together squares of such disparate orders of magnitude, the sum would be dominated by the early time channel measurements, with negligible contribution from the late time measurements. As a result, data 43 Chapter 3. Non-linear Parameter Estimation Procedure points corresponding to the very small components of the residual vector will essentially be ignored when finding search directions with which to reduce the objective function. Secondly, some observations may be less 'reliable' than others. In T D E M measurements, the level of noise may begin to dominate the level of signal from the target at later times. Figure 3.1 contains measurements of the field at different time channels for data taken over a buried 105 m m shell. The late time channel is polluted with background noise. Figure 3.2 contains the time decay at a station in the survey. The early time behaviour of d~B/dt is very noisy. In this case we would not want to have the fitting algorithm work excessively hard to try and fit each of these noisy data points very closely. If the forward modelling is linear or if there is a large number of data with normally distributed errors, then the optimal choice for Wd would be such that Wj Wd is the inverse of the covariance matrix of the errors (Bard, 1974). I choose Wd to be TO« = ( where cr, is the standard error of the i th to be some percentage of the i th 3 ' 9 ) data point and e is a small positive constant. I set cr; datum cri = % error X di. The off-diagonal elements of Wd are zero. The denominator ai + e is the error assigned to the i th data point. The addition of a positive e ensures that small data points would have a reasonable error assigned to it, and thus we can avoid biasing the inversion to the small data points. It is also possible to increase or decrease how closely an individual datum is fit. This is done using the data weighting matrix to increase or decrease the relative importance of the datum of interest. That is, to increase how closely we fit the i th observed datum d° increase the value of (Wd)a- A decrease i n how closely we fit the i th achieved by decreasing (Wd)a- bs we need only observed datum d° bs is 44 Chapter 3. Non-linear Parameter Estimation Procedure Figure 3.3 contains the synthetic decay curve for a 75 m m anti-tank mortar located 0.58 m beneath the surface. Each datum is assigned 10% error and an e = 5 X 1 0 " . Error bars 16 9B 9B X dt 0.01 1 3B y 0.01 100 1 Time (msec) Z at dt 0.01 100 Time (msec) 1 i 100 Time (msec) Figure 3.3: The synthetic dB/dt response curves for a 75 mm anti tank mortar buried 58 cm beneath the surface. The circled data points indicate that the absolute value of the data point has been plotted. at early times are small enough to not extend past the plotting symbols. A t late times the addition of e increases the size of the error bars. 3.2 V a r i a b l e s o f the I n v e r s i o n A l g o r i t h m The model vector m contains the 13 model parameters required for the approximate forward modelling: m = [X, Y, Z, a, c, k , a i , /3i, fi,k , x 2 a , /3 , 72] • 2 2 (3.10) where parameters (X, Y, Z) represent the location of the target, (a, c) represent the orientation, and the remaining parameters (A?i, a l t p\, 71, k , a , (3 , 72) are positive values that determine 2 2 2 the relative strength and time decay behaviour of the two induced dipoles. The following two sections describe two simple adjustments to the model vector that improve the chances for a successful inversion. 45 Chapter 3. Non-linear Parameter Estimation Procedure 3.2.1 Positivity Constraint In the forward model described by equation 2.38 the time decay parameters a;, and 7^ are all positive. Recall that the a - parameters appear in equation 2.38 as part of the time decay law 2 Li = k (t + ai)' 0 { expi-t/ji). If ai is perturbed during an iteration such that it is negative, Li would be a complex quantity for early times t. To force positivity on any of the parameters, the parameter fitting problem could be written as one of constrained optimization. That is, minimize <f>, subject to where > 0 (3-11) is any of the parameters upon which positivity is to be forced. A very simple alternate method is to enforce the positivity of this parameter is to define a new parameter m ; = ^/ra^. For example, if we would like to ensure that a; is positive, Li could be rewritten as Li = ki(t + di ) 2 0 exp(-t/n) The recovered ai is then the squared value of di. G i l l , Murray, and Wright (Gill et al., 1981) refer to the optimization method that uses this simple transformation as a square-variable unconstrained problem minimize <f> (rh), where <^(rh) is just (j> with m; = mi . The squared-variable transformation must be imple2 mented with care because non-hnear transformations of the parameters may lead to an increase in the non-linearity of the objective function. 3.2.2 S c a l i n g the M o d e l P a r a m e t e r s At each iteration of the algorithm, the value of || ( m t — m k i ) || is monitored such that we + can keep track of stalling or convergence as discussed in section 3.8. The model parameters of 46 Chapter 3. Non-linear Parameter Estimation Procedure Location and Orientation c X Y Z a 45° 45° 1 1 1 Dipole 1 Oil fa 7i 10 0.5 1 1 Figure 3.4: Values used for scaling the model vector m t y p fc 1 2 Dipole 2 02 «2 0.5 1 72 10 . These values were chosen after fitting decay curves for several different metallic objects. this problem can vary by a couple of orders of magnitude. Without scaling the components of the model vector differences, the algorithm may fail to see relatively large changes i n small model parameters due to the presence of large model parameters. This situation would lead to a premature termination of the algorithm. To address this situation the model parameters are scaled using a diagonal matrix W . m us first define a vector m t y p whose components m t yp i Let are positive scalar values representing typical magnitudes of the model parameters. When the values of the diagonal are chosen to be the inverse of m t y p , i.e. 1 typ (3.12) m. each component of the model parameter will make an equal, contribution to the relative difference in models. A new model vector rh = W m m typical magnitudes chosen for m t y p is denned. The tables in Figure 3.4 list the . The Hessian, and gradient for the transformed objective function are then H = W- H W- , 1 1 and (3.13) = W 7 V</> 1 In my tests I have not noticed any affect on the convergence rates by including W . m This is not surprising because Newton type steps are not changed by a linear scaling of the parameters (Dennis and Schnabel, 1983). 47 Chapter 3. Non-linear Parameter Estimation Procedure 3.3 Newton's M e t h o d of Optimization We wish to minimize the weighted least squares objective function ^=i||^(F[m]-d o b s 2' )|| (3.14) 2 For ease of notation the objective function is written as NN 1 1 4> (m) = -R ( m f R (m) = - £ 2 n (m) (3.15) 2 where R is the residual function R= W (F[m]-d ) d and ri (m) is the i ohs th component of the function R. The objective function <)> c a n he expanded i n a Taylor series to obtain a quadratic model of 4> </>(m + Sm) ss </>(m ) + V</> ( m ) £ m + ^ m V V ( m ) ^ m . T k k T k k A necessary condition for a minimum is to have the gradient of <f> vanish. (3.16) Applying this requirement to the quadratic model leads to V</> (rn,) 6m = - V<p (m*), (3.17) 2 where m * minimizes the local quadratic model and minimum of <p. is the Newton direction towards the The algorithm that solves this equation as part of minimizing <fr is termed Newton's Method of unconstrained minimization. The derivatives of <j> are V<j>(m) = J ( m ) R ( m ) T where S (m) is N S (m) = ^ »=i r; (m) V r ; (m) 2 and V V (m) = J ( m ) J (m) + S ( m ) , T (3.18) 48 Chapter 3. Non-linear Parameter Estimation Procedure and the Jacobian matrix J is defined as Jij(m) = Or • —-, omj where i = 1,.., N and j = 1,.., M. M is the number of model parameters. V</> (m) is called 2 the Hessian matrix, which we write as H. The matrix S contains second order derivatives, and thus second order information of the objective function i n the vicinity of m . The Newton step Sm for the non-linear least squares problem is then given by H ( m ) Sm = - J ( m ) R ( m ) (3.19) T k k k Newton's method is locally convergent. That is, if V(/>(m*) — 0 for some m*, then there exists an open set S such that for every starting point m 0 in S the Newton iterates will remain in S and converge to m*. Within S, Newton's method produces a quadratically convergent sequence of model iterates m ||m k + k i - m*|| < /3||m - m * | | (3.20) 2 k where {3 > 0 (More and Sorensen, 1984). Therefore, not only is Newton's method locally convergent, but it is also very efficient within the neighbourhood 5. The problem with implementing Newton's method as a global method is that the region S may be small. Although we try to make reasonable initial estimates of the model parameters (Section 3.5), there is no guarantee that an initial guess m D would fall within the domain of attraction S. However, Newton's method's excellent local behaviour suggests that it woidd be a good basis for a global method. The parameter estimation procedure will force convergence by modifying Newton's method such that m 0 is iteratively brought into 5 by a series of modified Newton's steps. Once the current iterate is brought inside 5, we then switch back to Newton's method to take advantage of its superior local convergence behaviour. Chapter 3. Non-linear Parameter Estimation Procedure 49 A n Algorithm For the Computation of Parameter Estimates 3.4 The goal is to minimize the sum of the squares of the residuals, represented by the objective function. The approach taken here is to make an initial guess of the model parameters, then to repeatedly improve this guess until the data predicted by our guess matches as closely as possible to the actual observed data. The basic steps of the unconstrained minimization algorithm used to minimize the objective function are 1. Choose a starting model m 0 2. Compute a search direction 6m. The search direction indicates the direction i n which to perturb the current guess. The search direction is chosen to be the Newton Step that minimizes the local quadratic model about the current iterate mk3. Compute a step length A. The step length indicates how much the current gtiess should be perturbed i n order to decrease the objective function. Because the local quadratic model about the current iterate is only approximate, the minimum of this model will not necessarily reduce the actual objective function. A positive scalar A is chosen such that 0(mk 4- X8m) < ^>(mk). A is found by a line search. 4. Update estimate of model. Set ra^i = mt + ASm. 5. Test for convergence. If the updated guess m ^ ! is adequate, then the algorithm is terminated. Otherwise, return to step 2. The steps are discussed i n more detail in the following sections. 3.5 Preliminary estimates of the model parameters The first step of this inversion procedure is to make a starting guess for the 13 model parameters: m = [X, Y, Z, a, c, k , a , /3 , 71,fc ,a x x X 2 2 ) /3 , 7 ]. 2 2 (3.21) Chapter 3. Non-linear Parameter Estimation Procedure 50 The success of an inversion procedure, as well as the rate of convergence towards a solution is dependent on the quality of the initial guess (Bard, 1974). The following sections will discuss a number of simple data preprocessing strategies that enables us to make a reasonable initial guess of the above parameters. 3.5.1 Analysis of Plotted Data Figure 3.1 contains plan view plots of different components of a 83/dt measurement over a 105 m m artillery shell. The shape of the measured signal, when not swamped by the presence of noise, appears as a series of distinct peaks. B y plotting these peaks and noting the position, width, and shape of these peaks, simple approximate relationships can be observed between the location of the target on the survey grid (X, Y), the the distance of the target from the transmitter (Z), and the strike of the target (a). Because these techniques rely on plotting the data i n plan view, the precision of these techniques will be directly impacted by the density of station distribution during data acquisition. E s t i m a t e s o f Target O r i e n t a t i o n (a, c) The parameter a of equation 2.38 gives the strike of the target. This strike is defined as the angle the axis of symmetry of the target makes with the x-axis of our survey. The strike of a target can be seen by considering the horizontal component of the response. The horizontal component is defined as (3.22) Figure 3.5 contain plots of the horizontal component of synthetic data for a steel rod, steel plate, and aluminum rod. The plan view response over a 2 m X 2 m area is plotted at three time channels after the primary field has been terminated. Overlaying the plan view plots is an arrow indicating the strike of each target. The changes in the contour pattern with each Chapter 3. Non-linear Parameter Estimation 51 Procedure time channel are due to relative strengths of the two dipoles changing with time. In each plot the response contains a line of symmetry which coincides with the target strike. The early time responses of the aluminum rod and steel plate do not have a very distinguishable line of symmetry. Indeed, finding a distinct line of symmetry can be very difficult, making this method susceptible to multiple solutions. Time = 0.0453 ms Steel Rod ~ a = 30 deg -S-2 c = 70 deg Time = 0.8025 ms Time = 17.19 ms 7 Steel Plate a = 30 deg c = 45 deg A l u m i n u m 2.5 Rod £ a = 110 deg 2 c = 15 deg i.5 F i g u r e 3.5: The determination of strike (a) from plan view plots of the horizontal component of the B-field data. The black arrows indicate the actual strike. The parameter c indicates the dip of the target. I have not found a method to estimate this parameter i n any simple way directly from the plotted data. Chapter 3. Non-linear Parameter Estimation Estimating the Horizontal Location Procedure 52 (X,Y) The location (X, Y) is the place on the survey grid directly above the buried target. Because this is the point on the survey that is closest in proximity to the target, it is reasonable to expect that a plan view map of the anomaly would mark (X, Y) as a maximum in the plotted signal. The method of finding the location is very straightforward in the case of a sphere. The secondary field of a sphere in a uniform field is equivalent to the field of a single dipole induced at the sphere center. The plan-view response of the vertical (Z) component response possesses a single, distinct peak. A n accurate estimate of the (X, Y) location of a sphere is simply the location of this 03/dt signal peak. W i t h the 2-dipole model, the situation becomes more complicated. For non-spherical targets the vertical component of the response loses its circular symmetry and, in some cases, there will exist 2 distinct peaks. The added complexity is due to the response being dependent on the relative strength of dipoles ki and k as well as the orientation of the target. In addition, 2 because the relative strengths of the dipoles change with time, at some time channels there may be only 1 peak while at other times there may be 2 peaks. To estimate the location of a target, first take a profile along the strike. Figure 3.6 shows the profile of a thin steel rod and a thin steel plate for different dip c, buried at X = 2m and at a depth of 1 m . The response is normalized by the maximum response along the hne. Because the rod and plate are thin, only one of the two dipoles makes an appreciable contribution to the response. The plan view responses, for these cases, do not change appreciably with time. We see that as the value of c approaches 90 degrees for a plate and 0 degrees for a rod, the profile has a second peak. A n estimate of location can then be obtained by either taking the location of a single peak (when there is a single distinct maximum), or by taking the midpoint of two peaks. For targets that have two dipoles making an appreciable contribution, cases occur 53 Chapter 3. Non-linear Parameter Estimation Procedure Thin Plate Thin Rod -a- c c c -*- c 1.5 2 X (meters) 1.5 2.£ 2 X (meters) = 0 deg = 22.5 = 45 = 67.5 :90 2.5 Figure 3.6: The vertical component of secondary field for a thin rod and thin plate with various c. These profiles were are from the time channel at t — 1.998 msec. where there may be 1 peak on the profile at one time and two peaks at another time. This is because the relative strengths of the dipoles change with time. Figure 3.7 shows estimates of the location at several times. The accuracy of observing the peak(s) of the vertical component of the field changes with time. On average, however, the location is relatively close to the real location of 2 m . 1r ^wfV (fl «•<» c o Q. DaC> 0.6 rh SW* ^ / Tt\\ °— G "* * x x > V \\ \© ® ITV" A— —^ A 0.0723 ms 0.2900 ms 1.0030 ms 4.0550 ms 17.1900 ms Wf Estimated (msec) Estimated Location Time on Profile (m) Depth TJ 0.0723 1.7857 1.1303 N 0.2900 1.6633 0.9467 O =3 ° - 4 0.2 0.5 1.5 2 2.5 X (meters) 1.0030 2.1531 0.9796 4.0550 2.1837 1.1633 17.190 2.1837 1.1633 3.5 Figure 3.7: Estimating the depth and location of a thick steel rod with c - 70 at various times. The real "Location on Profile" is actually 2 m, and the real depth is 1 m. In summary, the method I choose to determine the location of a buried target is to first find the strike (a) of the target. Secondly, I plot the profile along the strike of the vertical co: mponent of the response at different time channels. At each time channel I either find the Chapter 3. Non-linear Parameter Estimation Procedure 54 midpoint of two distinct peaks of the signal, or I find the maximum of a single peak. Averaging the location estimates of the different time channels result in a good estimate of the location. The accuracy of this technique is hindered when there is very sparse data coverage, and also when it is difficult to determine the line of symmetry along which to take a profile for analysis. Estimating Depth from Receiver (Z) A standard method for finding the depth of a buried target is to take the width of the signal at half the maximum value (Das et.al., 1990 and Nabighian and Macnae, 1991). The "Full W i d t h at Half Maximum" ( F W H M ) is defined in Figure 3.8, and provides an estimate of the distance from the receiver loop (Z). X (meters) Figure 3.8: Definition of the "Full Width at Half Maximum" (FWHM). The width of a profile provides an estimate of the depth. The F W H M value in this example is 0.82 m. Extra care must be taken when there are multiple peaks. Figure 3.9 shows the error of this estimate for different distances of a steel cube from the transmitter loop and along the symmetry axis of the loop. The widths at half-maximum for this plot was taken on synthetic data using a primary field generated by a 1 m square transmitter loop and a dipole transmitter with the same moment as the square loop. For objects without the symmetry of a sphere or cube, we lose the single distinct peak. When obtaining an estimate of depth I consider a profile of the vertical component of either the B-field or dB/dt data taken along the strike of the target. In the case of a single distinct peak, the width at half-maximum Chapter 3. Non-linear Parameter Estimation Procedure 55 Depth from Loop (cm) Figure 3.9: Errors in the depth estimate of a buried cube when using the F W H M . The F W H M estimates the depth of the cube to be deeper than the real depth when the cube is buried within approximately 105 cm from a square transmitter loop. At depths greater than 105 cm the F W H M estimates places the cube closer to the surface than its real depth. is used. In the case of two distinct peaks, the distance between the two peaks is used as an estimate. This process is tested on synthetic data of a steel rod at a dip of c = 70 degrees (Figure 3.5.1). Once again the accuracy of this estimate changes with each different time slice. As was suggested i n obtaining the location of the target, we can then take an average of the depth estimates over the different time channels. Results from performing this estimate on synthetic data produced for a steel rod and steel plate at depths of 1 and 2 m , and located at X = 2.0m are foimd in tables 3.5.1 and 3.5.1. These results suggest that this method of depth estimate provides us with a good guess of depth for different dips of either a steel plate of rod target.- The accuracy of this method deteriorates when the density of the data decreases. This method requires an estimate of the strike (a), which introduces the possibility of inaccuracy if the strike is hard to determine. Estimates of Time Decay Parameters (&, a, / 3 , and 7 ) One method for obtaining an initial estimate of the time decay parameters ki, ai, and 7 ; (i = 1,2) would be to consider values of the decay parameters of typical U X O and scrap found 56 Chapter 3. Non-linear Parameter Estimation Procedure Steel R o d ( 8 x 1 x 1 Depth = inches) Depth=1.5 1.0 m c (degrees) Depth (m) Estimated Location Estimated O n Profile (m) m Estimated Estimated Depth (m) Location O n Profile 0 1.1916 2.3030 1.3294 2.1400 15 1.1242 2.2147 1.5688 2.2787 30 1.0742 2.1330 1.5214 2.1920 45 1.0485 2.0720 1.4685 2.0783 60 1.0476 2.0107 1.4631 1.9920 75 1.1079 1.9700 1.5536 1.957 90 1.0345 2.0000 1.3788 2.0000 Steel Plate ( 8 x 8 x 1 Depth = inches) Depth=1.5 1.0 m c Estimated Depth (m) Estimated Location O n Profile (m) m Estimated Estimated (degrees) Depth (m) Location O n Profile 0.9478 1.9800 1.3302 1.9800 15 1.0951 2.1860 1.5451 2.2680 30 0.9519 2.1460 1.3430 2.2180 45 0.9125 2.1060 1.2879 2.1560 2.0920 0 (m) 60 0.9015 2.0580 1.2727 75 0.9075 2.0140 1.2812 2.0280 90 0.9146 1.9800 1.2912 1.9800 (m) Figure 3.10: Depth and location estimates for a steel rod and plate oriented at various dip. In both cases the "Estimated Location On Profile" is at 2 m. i n the area. For example, the Jefferson Proving Ground is the site of an ongoing extensive U X O cleanup project. A t this location, there are records of the U X O tested at this site. Appendix A and B contain the decay parameters of U X O and lab prisms obtained from lab measxirements at Geonics. More details of these measurements are found i n Chapter 4. A second method would be to solve a simpler problem that woidd provide estimates of the problem parameters. Let us consider data taken only directly above the target and let us assume that the response is due to a sphere (Figure 3.11). The primary field will induce a dipole parallel that is parallel to the primary field. The response of the dipole is then described by t(r,t) = M(r)k(t + ay^e-*^, (3.23) 57 Chapter 3. Non-linear Parameter Estimation Procedure Square Transmitter Loop Figure 3.11: Approximating the response of an arbitrary compact object by that of the sphere. Directly above a sphere, the response looks like a single decay dipole located at the sphere center and located at a distance r from the square transmitter loop. where M is the magnetic flux due to the induced dipole moment m i n d , and r is the distance from the loop center to the dipole measured along the axis of the loop. For the geometry of figure 3.11, equations 2.40 and 2.41 give a simple expression for M : M(r) = (3.24) ^ B 4irr 6 where B p is the primary field. When the transmitter is approximated by a vertical dipole with unit moment (e.g. l m x l m square transmitter coil with 1 A m p current). The primary field at the target would be only in the vertical direction with a magnitude given by ~p Mo 2 4ir r (3.25) 3 where z is the vertical unit vector. Substitution into equation 3.24 gives M 4 x IO" 1 4 (3.26) -z. In the geometry of figure 3.11, there will only be a vertical (z) component of secondary field at the surface. Therefore equation 3.23 can be rewritten as „6 (4-fM= t(i+a) " where r has been replaced with r e3t v •t/~f (3.27) to indicate an estimate of the depth (e.g. by using the width of signal at half the maximum value (see section 3.5.1)). We can fit a single curve with 58 Chapter 3. Non-linear Parameter Estimation Procedure k(t + ot.yP e~ ^ for estimates of the decay parameters. Figure 3.12 contains the results for 4/ trying this technique on a steel 8 x 1 x 1 rod oriented at various dip c. I choose to plot k/2 k/2 a 2 20 1.4 3, 1.3 15 1.2 a, "0 50 c (degrees) 0 100 50 c (degrees) 1.1 1 100 0 50 c (degrees) p 100 2 10 "0 50 100 c (degrees) Figure 3.12: Results from fitting the ^-component of the response measured directly over steel rod (8x1x1). rather than k because we need to estimate the values of k for two dipoles. When the target is parallel to the surface (c = 90), only dipole 1 is excited and the estimated parameters correspond to those of dipole 1. When the target is perpendicular to the surface (c = 0), only dipole 2 is excited and the estimated parameters correspond to those of dipole 2. The behaviour of k/2 does not show this characteristic because of the inaccuracy of the estimated depth. These plots indicate that if the vertical component of the response measured directly over the target can be described by k (t + a ) _ / 3 e" ^, 1 then an appropriate starting guess for the decay parameters would be ki — k/2, oti = a, pi = (3, and 7 i = 0.75 x 7 (3.28) where i = 1,2. 3.5.2 Fitting A Best Fit Single Dipole at Each Time Channel A second method of obtaining estimates of depth and orientation is to fit a best single dipole to data from a single time channel. The dipole is denned to be oriented along the axis of the U X O , and the dipole strength is proportional to the projection of the primary field along the that direction. The data in each time slice is inverted for a location (X,Y), a depth (Z), orientation 59 Chapter 3. Non-linear Parameter Estimation Procedure (a and c), and dipole strength (k). These parameters make up a model vector mdip = [X,Y, Z,a,c,k], and the model vector my; producing the best fit is obtained by minimizing the least-squares p objective function. This technique is tested using the test synthetic data set. A t each time channel there are 15 locations, and 3 components of the secondary field are used, which makes 45 total data with which to recover the 5 parameter model vector. The recovered model parameters obtained for each of the time channels are listed in table 3.1. The real model is listed at the bottom of Time(ms) 0.01 0.01778 0.05623 0.100 0.1778 0.3162 0.5623 1.0000 1.778 3.162 5.623 17.78 31.62 56.23 RealModel X Y Z a c 1.03092 1.00006 0.99290 0.99101 0.99055 0.99072 0.99139 0.99287 0.99337 0.99452 0.99533 0.99420 0.99624 0.99055 1.00499 1.00028 1.00230 1.00218 1.00146 1.00072 0.99981 0.99815 0.99914 0.99627 0.99607 0.99814 0.99640 1.00146 1.01695 1.00388 1.00375 1.00352 1.00324 1.00306 1.00300 1.00300 1.00256 1.00337 1.00292 1.00190 1.00132 1.00324 99.9 23.5 31.2 29.5 28.2 27.3 26.6 25.9 25.8 25.2 25.1 25.3 155.2 28.2 80.1 37.3 147.0 33.4 33.7 34.0 34.1 34.4 34.5 34.8 34.9 35.1 35.6 33.7 1.00 1.00 1.00 25 35 Table 3.1: Results from fitting a single dipole to some of the time channels of the synthetic data set. All lengths are listed in meters, and all angles are in degrees. the table. There is an excellent match between the recovered location and depth for each time channel and the real location and depth. The recovered orientation, denned by a and b, are plotted in Figure 3.13. The orientation is well recovered for all but the first couple (tc=l,2) and last couple (tc=16,17) of time channels. 60 Chapter 3. Non-linear Parameter Estimation Procedure tc=4 y z Figure 3.13: Orientation estimate of a UXO by finding a best fit dipole for each time channel, 'tc' indicates the which time channel was inverted to obtain the orientation. For all but the early and late time channels, the orientation of the UXO is well estimated by the best fit single dipole. 3.6 Computing a Search Direction Newton's step is derived by firstly making a quadratic approximation for the objective function to be minimized, and secondly requiring that the gradient of the objective function is zero. This is only a necessary and not a sufficient condition for finding the minimum of a function (Dennis and Schnabel, 1983). Consequently, the Newton step will take us to the critical point, and not specifically to the minimum, of the local quadratic model. The nature of the critical point can be determined from the eigenvalues of the Hessian matrix. If all eigenvalues are negative then the Hessian is negative definite and the Newton step will take current iterate to a maximum of the quadratic model. If some eigenvalues are negative and some positive, then the Newton step will lead to a saddle point. A Newton step will only lead to a minimum of the local quadratic model if the Hessian is positive definite (all eigenvalues positive). In this final case, the Newton step is a descent direction. Therefore, at each iteration a non-positive definite Hessian matrix should be perturbed in 61 Chapter 3. Non-linear Parameter Estimation Procedure some way to make it positive definite before solving equation 3.17. This change in the Hessian effectively changes the local quadratic model of equation 3.16 in such a way that the new model has a unique minimum. The location of this minimum defines the modified Newton step. This step can then be paired with a line search technique to produce a convergent sequence of iterates. Theorem (4.8) in (More and Sorensen, 1984) gives the strongest convergence result for line search methods. This theorem basically states that if a search direction 8m is defined as 8m = -Bk" V<?!> (m ), (3.29) 1 k where {B } is sequence of positive definite matrices with bounded condition numbers, then k V^(m ) <5m 1 T k pmjj > -11^(^)11 (3.30) where K is a positive constant. In the limit as k approaches infinity, the sequence {V</>(m )} k converges to zero. If a positive definite matrix B k with a bounded condition number is chosen at each iteration, then the iterations will converge to a critical point of the objective function «/>. The method i n which the Hessian is changed has no effect on the final recovered parameters, but only on the iterative route to the final parameters. Of course, an irresponsible choice of a new Hessian may lead to very slow convergence. For example, i f H is defined to be the identity matrix / , then the step defined by equation 3.17 becomes the direction of steepest descent, which has a notoriously slow rate of convergence for many problems. The method I adopted for perturbing H uses information about the local quadratic model from the spectral decomposition of the Hessian (Bard, 1974; G i l l et al., 1981). The eigenvalue decomposition of H is H = UAU r = Vdiag{[Xi])V T (3.31) 62 Chapter 3. Non-linear Parameter Estimation Procedure where the columns of U are the eigenvectors of H , and A; is the i eigenvalue of H. The th presence of a small eigenvalue A; corresponds to a ridge or a trough. The ridge or trough is parallel to the eigenvector corresponding to the small eigenvalue (Figure 3.14). This ridge (or Small Negative Eigenvalue i,„, * ( l n ) = o \\ ' X 3 Small Positive Eigenvalue ____ x • Newton Direction ^ Direction when taking absolute \ <Km) = 3 value of eigenvalue <t>(m) = .2 4>(m) = 2 \ 4.(111) = 1 -A \ (j)(m) = 1 — Figure 3.14: Small eigenvalues indicate either a ridge or a trough in the quadratic model of <f>. A small positive eigenvalue corresponds to a direction down a ridge. A small negative eigenvalue corresponds to a direction up a ridge. The direction can be reversed to point down the ridge by taking the absolute value of the eigenvalue as in equation 3.32 (After Bard, 1974). trough) is concave down if the eigenvalue is negative, and the corresponding step will be towards a maximum (or saddle point) of the ridge (or trough). If the eigenvalue is positive, the step is towards a minimum (or saddle point) of the ridge/trough. The small eigenvalues correspond to long steps along the ridge themselves. Therefore ridges provide an efficient path along which to minimize the function. B y taking the absolute value of the small negative eigenvalue, the corresponding direction points down the ridge to minimize the function. Therefore a new Hessian is defined as H = X3diag(\\i\)X5 2 (3.32) If this choice of H has a large condition number then a small positive number fi, defined as A m a i - ( A „ X maxcond) m i maxcond — 1 (3.33) is added along the diagonal. This choice of p will produce a Hessian with a condition number of maxcond. Chapter 3. Non-linear Parameter Estimation Procedure 63 T o demonstrate h o w this change i n the Hessian affects the N e w t o n step, let us consider a simple two-parameter p r o b l e m . L e t us assume that the current iterate o f some inversion procedure is = (xk, Vk) = (1, 3), a n d let us also assume that the l o c a l quadratic m o d e l about this point is given b y the function f(x,y) = lOx — y . A surface p l o t o f this f u n c t i o n is contained i n figure 3.15(a). A saddle point of this function exists at (0,0). T h e Hessian m a t r i x F i g u r e 3.15: Adjusting the Newton Step using Directional Discrimination, (a) A surface plot of the quadratic function f(x, y) = 10a; — y . There is a saddle point at the origin. A contour plot is projected below, (b) The Newton step takes a current estimate nik = (1,3) to the saddle point. The Directional Discrimination step takes the current estimate to = (0,6). The Directional Discrimination step adjusts the Newton step such that the new step is in a direction that reduces the objective function. 2 2 at the current iterate mt is then H = ( ° _° ), a n d the gradient is Vf(xk, t/j.) = [20, — 6] . 2 T 0 2 T h e eigenvectors o f the Hessian m a t r i x are V ! = [1, 0] a n d v T = [0, 1 ] . T h e corresponding T 2 eigenvalues are A i = 20 a n d A = - 2 , a n d therefore the m a t r i x i n indefinite. 2 T h e small negative eigenvalue A corresponds to a long step along the t r o u g h t h a t increases f(x,y). T h e 2 large p o s i t i v e eigenvalue A i corresponds t o a shorter step perpendicular t o the t r o u g h t h a t Chapter 3. Non-linear Parameter Estimation decreases f(x,y). £m= - H _ 64 Procedure The Newton direction is then 1 V / = [-1, - 3 ] (3.34) r This direction, as indicated in figure 3.15(b), takes the current iterate to the saddle point of f(x,y). If a new Hessian is defined following equation 3.32, then H = ( ^ i j ) ' a n ( ^ the n e w direction is then 8m = - H _ 1 V / = [-1, 3 ] . (3.35) T The new iterate is then m k + i = (0, 6). This modification of the Newton step has led to a direction that decreases the local quadratic model (Figure 3.15b). 3.7 Computing a Step Length The goal of each iteration in the model algorithm is to reduce the least squares data objective function c6. The Newton step defined by equation 3.19 provides a descent direction in which to perturb the current estimate. The local quadratic model of equation 3.16 is only an approximation, and therefore, a full Newton step to the minimum of the quadratic model will not necessarily reduce <f>. Because the Newton direction is a descent direction, we can look along the direction m k + aSm such that <j> (m + a8m) < <f> (m ). The positive scalar a is the k k step-length. Once a descent direction 8m is calculated, I first check if the full Newton step (a = 1) decreases the data misfit. Making the full Newton step the first option i n the step length procedure ensures that full steps are taken when the current iterate is near to the solution (m k £ S). This way we can take advantage of the excellent local convergence behaviour of Newton's method near the solution. If a reduction of (f> is not realized with a full Newton step, then the step length a is reduced until </>(mk + a8m) < <^>(m) is achieved. Let us define k p(a) = c/>(m + ap). k Chapter 3. Non-linear Parameter Estimation 65 Procedure To find a , a mixed polynomial/cubic step length algorithm is implemented (Dennis and Schnabel, 1983). After checking if the full Newton step reduces the objective function, there are three pieces of information that we have of p(ct): p(0) = </>(mk), p' (0) = V</>(mk), and p (1) — (f> (mfc + p). Therefore p (a) can be fit by a polynomial Pquad (ct) = aa + ba + c 2 where a = {p(l) — p(0) — p' (0)), b — p'(O), and c = p(0). The minimum of this quadratic function is occurs at 2(p(l)-p(0)-p'(0)) If a i m is very small, then it is possible that the quadratic function does not model p (a) very n accurately. Therefore a lower limit / is set for a. That is, i f a ; „ < /, then we instead choose m a = / and check if p(l) < p (0). Dennis and Schnabel suggests setting I = 0.10. If this choice of a fails to reduce the objective function, then we again fit p (a) with a polynomial. However, there is now three pieces of information about p. In addition to p (0) and p' (0), we know the last two values of p (a). Let us designate the last two values of a*, as a _ i and a _ . We can fit p (a) with a cubic function 2 Pcubic = aa + ba 3 2 (3.36) + ca + d where a a_i - a_ -A a- l 2 a" -B b= a_i - a_ 2 a' a' d = p(0) c = y(0) and ^ = p(a_ )- t,(0)-p (0)a_i, , 1 / B = p(a: )-p{0)-p'(0)a- . 2 2 The minimum of the cubic is then -b + \/b 2 a r 3a - 3ac (3.37) 66 Chapter 3. Non-linear Parameter Estimation Procedure 3.8 Stopping Criteria We must determine when to terminate our algorithm. Recall that a necessary condition for m to be a minimum of the data objective function <f> is that the gradient of the objective function, Vc/>, equals zero. Clearly, in finite precision math this requirement is unsuitable. The two convergence criteria used in this inversion process were suggested i n Dennis and Schnabel (1983). The first convergence criteria adjusts the necessary condition that the gradient of the objective function is zero. The size of a relative gradient measure must be less than some tolerance level V(/>(m i)- max{|(m i)-| ,typ m j max{|(?!>(mA. i)-| ,typ <p} fe+ max fe+ i = 1, < ^gradient, (3.38) + where nik+i is the new model, typ m; is the i th component of a vector of typical magnitudes for the model parameters, and typ <f> is a typical magnitude of (j). If each component of the sum of squares is scaled properly, the value of (j> at solution will be approximately equal to the number of data. Therefore typ <f> is set to the number of data. Dennis and Schnabel (1983) suggest choosing (gradient = (eps)^/ \ 3 eps is the floating point relative accuracy, and characterizes the machine precision, eps is the smallest number such that the machine recognizes that (1 + eps) > 1. A l l tests reported in this thesis were performed on a Sun SparcStation5 with an eps = 2.22 x 10 - 1 6 , and therefore (eps^ / ) = 6.06 X 10~ . This value of 1 3 6 (gradient is smaller than necessary, and results i n extra, unnecessary iterations of the algorithm. Therefore I choose (gradient = 1 X 10 The second convergence criteria checks whether the algorithm has either stalled or converged by monitoring the relative change in successive models. When the size of the relative change in model is less than some predefined tolerance the algorithm is terminated. This criteria Chapter 3. Non-linear Parameter Estimation 67 Procedure is given by (m fe+1 )- - (m ). (3.39) fe max-lKmfe+ijJ , typ mi} I follow the suggestion from Dennis and Schnabel (1983) of selecting (model = 1 0 _ p where p is the number of significant digits desired at solution. 3.9 Error Bounds of the Parameter Estimates Once the model parameters m + which minimizes the objective function (p (m) have been ob- tained, we must examine the reliability and precision of the estimated parameters. Meas- urements are random and data are noisy. Thus it is not sufficient to obtain a set of model parameter estimates m + and claim that these parameters are the best estimates of the un- known parameters m*. The parameters that may best describe one measurement, may indeed be different than the parameters obtained from a second measurement on the same sample U X O . In Chapter 4 a U X O identification algorithm is defined which rely on the values of the recovered parameters. Therefore it is essential to understand how the parameters in our model respond to errors i n the data set i n order to establish how much confidence can be placed on the recovered model. The gradient of the objective is zero at a minimum The objective function <p(m) is written as <f> (m*,d ) to explicitly state the dependence of obs the objective function on the observed data. If the observed data were slightly perturbed, the location of the minimum would also be slightly perturbed. The gradient at this shifted minimum would then be zero: (3.40) Chapter 3. Non-linear Parameter Estimation 68 Procedure The gradient can then be expanded as a Taylor Series V<t> ( m , + Sm*, d° + Sd° ) bs = V</> ( m „ bs df ) s + V V (m*, < ) 4- b s _ 9^V</>(m*,dfs) dd where H.O.T. represents the higher order terms. When the objective function <j> is the sum of squares, the gradient of <j> is V4> = J T (V[m] - d obs ) . (3.42) Therefore the derivative of V(/> with respect to d obs is simply the Jacobian matrix J* . B y T subtracting equation 3.40 and retaining only first order terms we get Sm* K, H * " J* T 1 £d obs (3.43) where H * is the Hessian matrix evaluated at the minimum of <fi ( H * = V (j> (m*)). 2 The model covariance matrix V is defined as the expectation value of Sm* Sm* T m V = E (Sm* Sm* ) . T m Substitution of equation 3.43 gives V = E(H,- J / 1 m £d obs £d o b s T J , H *- l The Hessian and Jacobian i n the above expression are evaluated at m = m * , and are therefore constants. As a result they can be taken outside of the expectation value expression: V = R* 1 m Jj V J.Hr 1 d (3.44) where V is the covariance matrix of the data d V = E (Vd d obs Sd ^j obsT (3.45) Chapter 3. Non-linear Parameter Estimation 69 Procedure In the case when the observations all have normally distributed, uncorrelated errors the data covariance matrix reduces to a diagonal matrix Therefore, an estimate of the standard deviation of the i th (V )u = c (H* 2 m 1 model parameter m ; is then (3.46) J* J*H* T Model variance estimates applied to non-linear problems are not as reliable as when implemented in linear least squares problems, and should only be used as a very rough estimate (Bard, 1974; Dennis and Schnabel, 1983). 3.10 Application to a Synthetic Data Set The parameter estimation is now tested on a synthetically generated field data set. For this data set I will forward model the model parameters of a 75 m m anti-tank mortar buried at a depth of 67 cm (Z = l m ) , and located at (X,Y) — (2,2)m on the survey grid. The mortar is oriented such that a = 25 °, and c = 35 °. This survey consists of a 2m X 2m grid, containing 5 lines of data with 5 stations per line. At each station 3 components of the dB/dt response is generated for 17 logarithmically space time channels. The time channels start from 0.01 msec and end with 100 msec. The data set has 10% random Gaussian noise added to it. The residual in the objective function is multiplied by a data weighting matrix Wj = l / ( c - f e), where e = 1 X 1 0 ~ . Results 17 of the predicted and observed data are found on the following pages (Figures 3.16 to 3.19). F i g u r e 3.16: X component of synthetically generated dfl/dt data for a 75 m m mortar. The predicted and observed data are a very good match, with no correlated to noise. F i g u r e 3.17: Y component of synthetically generated dB/dt data for a 75 m m mortar. The predicted and observed data are a very good match, with no correlated noise. Chapter 3. Non-linear Parameter Estimation Procedure 72 F i g u r e 3.18: Z component of synthetically generated dB/dt data for a 75 m m mortar. The predicted and observed data are a very good match, with no correlated noise. Chapter 3. Non-linear Parameter Estimation Procedure 73 Time (msec) Figure 3.19: These three plots demonstrate that the predicted data matches well with the observed data. The final x 2 at solution is equal to 913.2 which is slightly less than the number of data (iV = 1275). The recovered parameters are indeed good. 3.11 Some Tests of the Parameter Estimation Procedure In this section several practical aspects of the least squares procedure are considered. These tests help to determine the exact form of the objective function, and to determine constants required by the algorithm. In particular, these tests should provide us with a feel as to when the parameter estimation might fail to converge. These tests will, in general, use the parameters Chapter 3. Non-linear Parameter Estimation m X Y Z a c h ai Pi 7i fc 2 «2 ft 72 true 2.00000 2.00000 1.00000 25 35 12.01800 0.00760 0.89310 17.65300 3.29960 0.00766 1.25230 11.53850 74 Procedure est. standard deviation 1.99974 0.00016 1.99998 0.00006 0.00226 1.00148 24.82 0.236 34.79 0.162 0.10108 11.75858 0.00748 0.00052 0.89157 0.00385 17.81922 0.14379 0.03450 3.26263 0.00712 0.00037 1.24198 0.00541 10.61993 0.24369 Table 3.2: Recovered parameters for the synthetic data set inversion. The true model m recovered model m are very close. The estimated standard deviation is very small. t T u e and the r e c for the (9B/dt response of a 75 mm anti-tank mortar. Test 1: Determining a Maximum Step Length The directional discrimination step guarantees a descent direction by forcing the Hessian to be positive-definite. This step may be too large and move the current iterate outside an area where the quadratic model is applicable. To avoid steps that are too large, a maximum step length maxstep is defined. Dennis and Schnabel (1983) suggest restricting the step length to be less than 1000 X ||mt ||, where m t yp y p is a vector of typical model parameters. I have found that the steps allowed by this size of maxstep are too large for this problem. Figure 3.20 contain result from inverting the 75 mm synthetic data set with maxstep = 1000 x 11m 11. Panels (a), typ (b), and (c) contain plots of the progress of the model parameters (scaled by the value of the real model parameter) as a function of iteration. The recovered parameters fail to converge to the correct model parameters. The recovery of a is particularly poor in this example. Panel 2 Chapter 3. Non-linear Parameter Estimation 75 Procedure (d) plots the values of the objective function with each iteration. The inversion was unable to reduce the objective function to the level that would be achieved if the real model parameters were recovered. Figure 3.21 contains results from inverting the 75 m m synthetic data set with maxstep = 10 X • ||m yp||• The parameters recovered from this inversion are correct. For the t inversion algorithm, the maximum step parameter is set to a conservative maxstep — 10 X 11mtyp11a) 8 10 Iterations Dipole 2 Parameters b) Dipole 1 Parameters 6 0 (3-47) O ^ O 6 O j 12 8 H O |) O j> 8 Objective Function Location (X,Y,Z) and Orientation (a,c) G—O |A—A B O ^ O j> O B X Y Z a e-o o o o aA — • .•—B—a—a—B—ffl ~~g"A~~A"A"A~A"& i—i—i—i—i—i—h 8 10 Iterations 12 14 16 Figure 3.20: Convergence results when maxstep - 1000 x ||m yp||. The recovered parameters fail to converge to the correct model parameters. The recovery of a is particularly poor in this example. Panel (d) plots the values of the objective function with each iteration. The inversion was unable to reduce the objective function to the level that would be achieved if the real model parameters were recovered. t 2 Chapter 3. Non-linear Parameter Estimation c) 76 Procedure Objective Function Location (X,Y,Z) and Orientation (a,c) Iterations F i g u r e 3.21: Convergence results when maxstep = 10 x ||m yp||. The correct model parameters are t recovered by the inversion. Panel (d) plots the values of the objective function with each iteration. The inversion achieved a data misfit that is less than the objective function to the level that would be achieved if the real model parameters were recovered. This indicates that some of the noise is being fit. Chapter 3. Non-linear Parameter Estimation 77 Procedure Test 2: Sensitivity of P a r a m e t e r s to e In Section 3.1 e was introduced into the data weighting matrix in order to remove any bias to fitting the small magnitude data points. Increasing the size of the error bars reduces how closely we must fit the smaller magnitude data points, which will clearly impact the accuracy of the recovered parameters. To examine how the choice of e affects the recovered parameters, I perform several inversions on synthetically generated data for a 75 m m anti-tank mortar. The data are contaminated with 10% Gaussian noise. The values of e for these tests ranges from 0 to 1 X 1 0 ~ . Recall that the minimum value of <9B/dt for this synthetic data set was set to 15 1 x 1 0 ~ . The results from these tests are reported in figure 3.22. The percent error in each 17 dipole 2 parameters dipole 1 parameters Location and Orientation 20 A k, * a, • P, Or, o ^> 15 10 5 CJ 0 1e-17 1e-16 1e-15 e Figure 3.22: The recovered parameters for different choices of e. Recall that for this data set that the minimum amplitude datum was set to 1 x 10~ V/m . Increases in e are reduces the accuracy with which we can recover the model parameters. The 7 parameter is the most sensitive to changes in e. 17 2 recovered parameter, defined as % error = \m; ,recovered I m- x 100%, \mis plotted as a function of e. W i t h increasing e there are accompanying increases in the percent error of the recovered model parameters. For parameters other than 71 and 72 there is no more than 4% error in parameter recovery (even for the unreasonably large e = 1 X 1 0 ~ ) . The 7 15 parameters are most sensitive to the e setting. This is not surprising because 7 is a measure of the onset of the late-time exponential behaviour of the dipole decay. The late-time data points are the smaller data points, and changes in e change how the inversion fits these small Chapter 3. Non-linear Parameter Estimation Procedure 78 data points to the greatest extent. The predicted standard deviation of the model parameters obtained from the model covariance matrix also increases with increasing e. Test 3: P o s i t i v i t y C o n s t r a i n t In Section 3.2.1 I discussed how a model parameter can be forced to be positive by defining a new parameter that is the positive square root of the original parameter. Any non-linear transformation on a model parameter would be expected to produce a more non-linear objective function. Tests must be done on specific problems to guarantee that the transformation does not cause problems with the chosen algorithm. In Figure 3.23, inversion results are presented using different transformations on the model parameters. The plots include convergence results when only the a parameters are squared, to when all the time decay parameters are squared. In Figure 3.23(a), each of the different transformations result i n recovering all the correct parameters. However, when the starting guess for the algorithm is located further from the location of the correct model, some of the different transformations result i n the algorithm converging to the incorrect model (Figure 3.23(b)). In the the final version of the algorithm, the transformation is applied only to the a; parameters. Test 4: S e n s i t i v i t y o f the S t a r t i n g Guess m 0 The success of an inversion procedure, as well as the rate of convergence towards a solution is dependent on the quality of the initial guess (Bard, 1974). Therefore it is appropriate to test how sensitive the success of the inversion procedure is to the initial guess. There is no systematic way to perform a test varying different parameters and different combinations of parameters due to the 13 degrees of freedom in the model. A test suggested in Dennis et al. (1981) is to try initial guesses that are further out on the ray from the solution m * to the standard starting point m . Figures 3.25and 3.24 contain results from a slightly modified D version of this test. Because preliminary estimates of the location and depth are qiute good, Chapter 3. 79 Non-linear Parameter Estimation Procedure a) ' 10 Objective Function Objective Function 7 * • |A G— * • A O 4 L(t) = k L(t) = k L(t) = k L(t) = k * • (t+a") exp(-t/y) (t+a ) " exp(-tM (t+a ) exp(-t/y ' (t+a )"" exp(-t/y ) 2 2 2 2 P 2 L(t) = k (t+a") * exp(-tfy) L(t) = k (t+a ) " exp(-t/y) L(t) = k (t+a ) exp(-t/y ) L(t) = k ( t + a ) ' exp(-t/y ) 2 2 2 2 2 2 -|i 2 2 2 2 * • 2 G 2 G 10 6 Iterations P 2 15 Iterations Figure 3.23: Convergence behaviour when using various variable transformations, (a) When using the standard starting guess, each of the different transformations on the model parameters converge to the correct parameters, (b) When the starting guess is taken further away from the solution (using equation 3.48, transformations to a more non-linear objective function lead to a cases where the inversion doesn't converge to the right answer. and because the orientation parameters are angles, only the time decay starting parameters will be varied i n these tests. That is, the starting guess is varied according to m i = 1,...,5, test (m where m est 0 oi - m»i)K - m*i (3.48) i = 6,...,13 is the new starting guess, m is the standard starting guess developed in the 0 previous section, m* is the real model, and K controls the distance along the ray from m* to Figure 3.24 shows convergence results when K = 2. Panels (a), (b), and (c) contain plots of the progress of the model parameters (scaled by the value of the real model parameter) as a function of iteration. The recovered parameters converge to the correct model parameters within 17 iterations. Panel (d) plots the values of the objective function with each iteration. The inversion achieved a final below the value of the objective function calculated with the true model. Therefore the parameterization is fitting some of the 10% Gaussian noise that was added to the model. Of course, the extent from which one can deviate from the standard Chapter 3. Non-linear Parameter Estimation a) Dipole 1 Parameters 80 Procedure b) Location (X,Y,Z) and Orientation (a,c) Dipole 2 Parameters Objective Function Iterations F i g u r e 3.24: Convergence results when K = 2. Panels (a), (b), and (c) contain plots of the progress of the model parameters (scaled by the value of the real model parameter) as a function of iteration. The recovered parameters converge to the correct model parameters within 17 iterations. Panel (d) plots the values of the objective function with each iteration. The inversion achieved a final data misfit less than would be achieved when calculated with the true model. Chapter 3. Non-linear Parameter Estimation 81 Procedure guess, while still successfully recovering parameters, is limited, and it is possible to make the algorithm fail. Figure 3.25 shows convergence results when K = 5. The recovered parameters fail to converge to the correct model parameters. The recovery of ct is particularly poor in this 2 example. Panel (d) indicates that the inversion was unable to reduce the objective function to the level that would be achieved if the real model parameters were recovered. Dipole 1 Parameters Dipole 2 Parameters h) 1800 1600 3000 1400 2500 1200 1000 m , 2000 m , 800 1500 600 1000 real real 400 500 200 ~15~ ~ 20 Iterations c) 15 25 Location (X,Y,Z) and Orientation (a,c) 20 Iterations 25 30 Objective Function Figure 3.25: Convergence results when K = 5. The recovered parameters fail to converge to the correct model parameters. The recovery of a is particularly poor. Panel (d) indicates that the inversion was unable to reduce the objective function to the level that would be achieved if the real model parameters were recovered. 2 82 Chapter 3. Non-linear Parameter Estimation Procedure Test 5: How well can we predict the error? To test how well we can predict the model parameter error for a solution m , the synthetic + data set generated using parameters for a 75 m m anti-tank mortar is inverted. The difference in the recovered and real parameters are then compared to the error predicted by the model covariance of equation 3.46. Inversions are performed with 5%, 10%, 15%, and 20% noise added to the data. The results are found in table 3.3 and figure 3.26. mi 5 % noise a error X Y Z a c 0.00002 0.00001 0.00023 0.00406 0.00019 0.00010 0.00004 0.00151 0.00291 0.00198 0.00005 0.00002 0.00010 0.00682 0.00052 0.00018 0.00007 0.00285 0.00505 0.00355 0.00007 0.00003 0.00008 0.00950 0.00088 0.00027 0.00009 0.00417 0.00699 0.00506 20 % noise error a 0.00009 0.00038 0.00005 0.00012 0.00014 0.00551 0.01282 0.00887 0.00120 0.00654 ki 0.05772 0.00003 0,00174 0.04067 0.06557 0.00035 0.00272 0.06751 0.22391 0.00006 0.00390 0.04177 0.11822 0.00067 0.00503 0.10890 0.50605 0.00021 0.00574 0.07263 0.16705 0.00101 0.00735 0.14853 0.91800 0.00062 0.00682 0.15344 0.01604 0.00013 0.00383 0.01987 0.02342 0.00025 0.00378 0.13031 0.06644 0.00029 0.00799 0.03611 0.04006 0.00048 0.00676 0.19509 0.15046 0.00044 0.01260 0.07983 0.05505 0.00071 0.00966 0.25418 0.27435 0.00055 0.01784 0.17695 Oil ft 7i k 2 OL 2 ?2 72 10 % noise CT error 15 % noise error a 0.21205 0.00137 0.00975 0.18869 0.06854 0.00094 0.01258 0.31472 Table 3.3: Errors in the recovery of parameters for a buried 75 mm anti-tank mortar, 'error' is the absolute value of the difference between the real and recovered parameter, a is the standard deviation of the recovered parameter predicted from the model covariance matrix (equation 3.46). 3.12 Selection of Data in a Field Survey The EM61-3D instrument can collect a large amount of data. A n EM61-3D data set collected over a U X O buried at the York University Geophysical Test Site contained over 13 000 different datum. It would be beneficial to identify which data are essential for parameter recovery, and which data provide redundant information. This would result in an immediate reduction in 83 Chapter 3. Non-linear Parameter Estimation Procedure 20 percent noise 5 percent noise 0.2 a 0.15 error 0.1 0.05 0 o o X Y Z a c k a P i i i Y i k a P v 2 2 2 Model Parameter X Y Z a c 2 k a p i i i Y i k a P r 2 2 2 2 Model Parameter Figure 3.26: The error in recovered parameters for a buried 75mm anti-tank mortar, 'error' is the absolute value of the difference between the real and recovered parameter, cr is the standard deviation of the recovered parameter predicted from the model covariance matrix (equation 3.46). The predicted error cr is a reasonable estimate for the actual difference in real and recovered parameter. the time required for the parameter estimation algorithm to converge, and also save time when collecting data from the field. 3.12.1 T i m e Channels The T D E M response is measured over a number of time channels. In the EM61-3D, the time decay response at each station is measured over 30 logarithmically spaced time channels. For an example of the measured values see figure 3.2. Thirty time channels may produce redundant information, and the data may be well represented by only a portion of the time channels. If the inversion could recover correct parameters using only half of these time channels, the number of data would be halved and the speed in which the computer code requires to complete an iteration would be increased. This section tests if parameters can accurately be recovered by using only some of the time channels. Figure 3.27 includes results from repeating inversions of synthetic 75 m m mortar dB/dt data using fewer time channels. In figure 3.27(a) results from using of the 17 time channels Chapter 3. Non-linear Parameter Estimation 84 Procedure are used. In figure 3.27(b) data from only the odd numbered time channels (time channels 1, 3, 17) are used. In figure 3.27(c) every third time channel (time channels 1, 4, 7, 10, 13, and 16) are used. For these three instances the inversion converged with the recovered model parameters being quite good (Table 3.4). In figure 3.27(d) (time channels 1, 5, 9, 13, 17) are used. In this final case the inversion did not converge and was therefore terminated at 100 iterations (of which the first 47 are plotted here). True Model 1.00000 1.00000 1.00000 25° 35° A l l Time Channels 1.00005 1.00002 1.00010 25.39° 34.97° Every 3rd Odd Time Channels 0.99979 1.00014 0.99602 24.97° 35.00° Every 3rd Time Channel (1,4, ...,16) 1.00054 0.99977 0.99951 25.45° 35.06° Time Channel (1, 5, 17) 0.99934 0.99976 0.99641 24.95° 35.41° 12.01800 0.00760 0.89310 17.65300 11.79409 0.00766 0.88920 17.69477 11.56766 0.00792 0.89354 17.93139 11.77030 0.00680 0.88834 17.62755 11.00768 0.00986 0.93181 21.74450 3.29960 Ot-2 0.00766 1.25230 02 11.53850 72 3.23316 0.00736 1.24431 11.57461 3.20491 0.00763 1.24714 11.47068 3.20602 0.00720 1.24941 11.41476 2.67732 0.01347 1.39148 -1122.20051 mi X Y Z a c Pi 7i k 2 Table 3.4: Recovered parameters when using different time channels. When enough of the time channels are omitted, the parameter estimation starts to fail. The results of fitting the 3 components of dB/dt at station ( X , Y ) = (0,0) are shown i n figure 3.28. As was expected, the data were best fit when all of the time channels are included, and the data was fit the poorest for the case of the fewest time channels being utilized. In particular, the early time channels could only be fit when the all the time channels were being used. In the case when only time channels 1, 5, 9, 13, and 17 were used, both the the early and late times were not fit by the predicted data. Chapter 3. Non-linear Parameter Estimation a) b) Model Parameters 85 Procedure Model Parameters G G X A—A Y Q • Z 1 1 c G O k, A A ai • * • * G--0 A--A B - -O X G A • 1 2 X O A o * * P, y, k G A a V G-O A--A • - -• ^ - -X- 1 / bifF m m 9 P a • • *• - -* 1z 2 3 G A— • — H X Y Z a c . k, «, P; y. k °? p 2 2 ft 4 5 Iterations Model Parameters c). 20 25 Iterations 30 F i g u r e 3.27: Using only a portion of the time channels, (a) A l l time channels are used, (b) Only odd time channels are used, (c) Time channels 1, 4, 7, 10, 13, and 16 are used, (d) Time channels 1, 5, 9, 13, and 17 are used. F r o m these results, 30 t i m e channels is more t h a n enough for accurate recovery of m o d e l parameters. 3.12.2 Spatial Coverage The spatial density of collected data during a T D E M survey is limited only by how accurately each station can be surveyed and the patience of the data collector. The synthetic field example of a 75 m m mortar that is used throughout this thesis represents a reasonable survey design. 86 Chapter 3. Non-linear Parameter Estimation Procedure b) X component of dB/dt o o ., Y component of dB/dt observed data All time channels Odd time channels time channels 1, 4, 7,10, 13,16 lime channels 1, 5, 9, 13, 17 C) ' Z component of dB/dt -„ 10 io"' 10"' 2 3 10"" co 10 10 H 6 10"" IO"' 8 observed data All time channels Odd time channels time channels 1, 4, 7, 10,13, 16 time channels 1, 5, 9, 13,17 10 10"' 10 10 Time (msec) Figure 3.28: Observed and predicted data for inversions using only a subset of the available time channels. Panels (a), (b), and (c) contain the fitted responses of the x, y, and z components, respectively, of the dB/dt response. The quality of fit decreases when fewer time channels are being used. A n increase i n both the field collection procedure and inversion computation time would be achieved i f an idea of some lower limit of required stations could be established. Figure 3.29 illustrates some of the different survey designs. In each of the panels, ' x ' marks the location of a station. In each case, an accurate recovery of the model parameters was achieved when using all three components of data. However, this was not the case when only the Z-component of 83/dt was considered. When only the Z-component is included as data, the spatial coverage of the survey becomes an issue. As an example, let us consider data taken 87 Chapter 3. Non-linear Parameter Estimation Procedure along Line X = 0 m . If Z-component time decay measurements are taken at 15 evenly space stations (every 14.3 cm) along this line and inverted, the inversion converges to the incorrect result. These results are plotted in figure 3.30. Panels (a) and (b) include the observed and predicted data from both early and late times, indicating a good data fit in both cases. Panels (c), (d), and (e) include the model parameters as a function of iteration. It is clear that the inversion has been unsuccessful in recovering the correct model parameters even though the data have been well fit. Panel (f) shows that the final data misfit is lower than when the data misfit is calculated using the real model. Correct convergence using only Z-component data occurs when stations surround the (X,Y) location of the target at the surface. Correct convergence for the inversion occurred for survey layout (a), (c), (e), and (f) of figure 3.29. In each of these layouts stations are distributed across the surface, as opposed to being restricted to one line. Convergence did not occur for layout (b) most hkely because there wasn't enough stations. a) X x = lm X - x x _ „ _ X -•x- : - X 1. - -X x -x X X 1 b) x = 2m X= Im 'i^ j.-X 7 d) c) * » X 5 X - - X "tm v— v X- > 1 . X - . x 0" X f) x = Om ^ x X „r * - -x • X—~—i » — L x Figure 3.29: Different survey configurations used to test how many stations are required for accurate model parameter recovery. The ' x' symbol marks the location of a station. Inversion attempts using only the Z-component of data resulted in correct convergence for survey configurations (a), (c), (e), and (f)- Chapter 3. Non-linear Parameter Estimation Procedure 3.13 88 Summary In this chapter a non-linear parameter was developed to recover the 13 parameters of the model. Rough estimates of the model parameters were obtained by introducing a number of simple data prepocessing strategies. These estimates were then refined by using a Newton's Method to minimize the least squares objective function. This process was shown to work well in the presence of random noise, and for a variety of data gathering configurations. In particular, the type and amount of data collected by EM61-3D (30 time channels, 3 components, 50 cm line spacing and 25 cm station spacing) is more than adequate to recover model parameters i n synthetic data inversions. Chapter 3. Non-linear Parameter Estimation 89 Procedure Time Channels 11, 12, 13, and 14 o t = 3.16 m s a t = 5.62 m s Station (m) Station (m) Figure 3.30: Results from inverting only the Z-component of the dB/dt data of Line X = 0 m . Panels ' (a) and (b) include the observed and predicted data from both early and late times, indicating a good data fit in both cases. Panels (c), (d), and (e) include the model parameters as a function of iteration. It is clear that the inversion has been unable to successfully recover the correct model parameters even though the data has been well fit. Panel (f) shows hat the final data misfit is lower than when the data misfit is calculated using the real model. Chapter 4 Relating Physical Parameters of a Metallic Target to Model Parameters U X O are typically rod-like rather than plate-like. This chapter examines how the time domain electromagnetic signal for a target varies with its shape. In the absence of analytic solutions for Maxwell's equations and the lack of suitable numerical electromagnetic forward modelling, I use lab measurements of the time decay curves of different metallic targets. Time domain curves for various sized steel and aluminum rectangular prisms compiled at Geonics L t d . were supplied for analysis. From this analysis I hope to use the different characteristics of the time domain curves to relate to them to the material and shape of the target. Ultimately, the goal of this analysis is to produce relationships between model parameters of the approximate forward modelling representation and actual, real parameters of the measured target. 4.1 The Effect of Shape on the T D E M Response 4.1.1 Lab Setup The experimental setup used by Geonics L t d . to obtain the decay response for different sized prisms and for U X O is illustrated in figure 4.1. A square transmitter loop with sides of 40 m was used to provide a relatively uniform field at the center of the loop. A i m diameter receiver coil was placed coaxial to the transmitter loop. Each target was placed at the center of the receiver loop. Geonics standard P R O T E M 47 time domain equipment was used for producing the transmitting field and for recording the time domain measurement. Of the targets measured using the experimental setup of Figure 4.1, included a series of metallic rectangular prisms of different dimensions. Each prism had at least one dimension of 90 Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 40 m ~H TEM47 Transmitter Figure 4.1: Geonics Laboratory Setup 8 inches and the targets ranged from a thin rod ( 8 x 1 / 4 x 1 / 4 inch) to a cube ( 8 x 8 x 8 inch) to a thin plate ( 8 x 8 x 1 / 4 inch). The suite of measurements included targets with no axis of symmetry (3 different dimensions). Measurements were made of steel and aluminum prisms. The majority of U X O are made of steel due to its strength. Steel is permeable (p = 150p ) o as well as conductive (cr = 0.67 X 10 5m). Aluminum is non-permeable (p = p ) and slightly 7 0 more conductive than steel (cr = 3.54 x 10 5m). 7 In addition to the sets of steel and aluminum prisms, a set of 24 sample U X O were measured. These targets included various ordnance items used by N A T O since World War II. The ordnance range in length from 18 to 85 cm, and in diameter from 6.05 to 15.92 cm. A diagram of all the ordnance, along with a table listing the dimensions of each ordnance, are included in appendix A . For the analysis of this thesis we consider only the measurements of axi-symmetric targets. The axi-symmetric targets were placed at the center of the receiver loop i n two orientations. The target was measured with the axis of symmetry perpendicular and parallel to the primary field (Figure 4.2). Measurements of the time decay response of these targets were recorded as plots of log(5B/(9t) vs. log(i). Since values were not recorded by a data logger, plots were digitized. Plots of the steel target responses were digitized by hand by J . D . McNeill at Geonics L t d . , and the aluminum targets were digitized by scanning the plots into Matlab. The digitization process i n both cases will have some error associated with it. Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters B Rod Orientation B Plate Orientation 1 1 B Plate Orientation Rod Orientation B 2 2 Figure 4.2: Different target orientations for Geonics laboratory measurements. In orientation 1, the axis of symmetry is parallel to the primary field B . In orientation 2, the axis of symmetry is perpendicular p to the primary field B . p 4.1.2 C o m p u t i n g B - f i e l d responses f r o m dB/dt D a t a The dB/dt data were also integrated to obtain the B-field response of the target using where dB/dt is the values measured by the lab equipment, and t is the time of the last meas0 urement. The integral from t to t is computed using a trapezoid ride. A spline interpolation 0 was applied to the data to increase the number of points used when applying the trapezoid integration. The second term of equation 4.1 is also calculated by integration B(t ) = jf ^-—jdt 0 The value of B ( t / a i e + B(t ) (4.2) late ) is not known. However, for large enough ti , B(ti ) ate When computing this integral I assume that ti ate Cl00t B(t )« D o dt J Jto is essentially zero. = 100£ , and approximate the integral as a / ate (4.3) This integration in carried out by first fitting the dB/dt data with (4.4) L(t) = k [t + a)~ fi e-'/T Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters The late time (t > t ) is then extrapolated using 4.4 Q T D E M Response of a Plate 4.1.3 Let us consider the response of a metal plate. Figure 4.4 contains the T D E M response of a steel plate and an aluminum plate. The two main physical processes that produce the decay responses plotted i n figure 4.4 are eddy current induction and magnetic induction. Aluminum is a conductive, non-ferrous metal, and therefore the response of an ahmiinum target will only exhibit the characteristics of eddy current induction. Steel is a conductive, ferrous metal, and therefore steel targets will experience both eddy current induction and magnetic induction. E d d y C u r r e n t Response The induction of eddy currents i n the presence of a changing magnetic field is predicted by Faraday's Law. When the primary field is on, magnetic flux passes through the plate. Referring to figure 4.3, the plate is said to be maximally coupled to the primary field when i n orientation 1. Plate Orientation 1 B Plate Orientation 2 B P ! \ Primary Field On Primary Field OffEddy Currents Induced Primary Field On Primary Field Off Eddy Currents Induced Figure 4.3: Eddy current induction in a non-permeable plate. When the direction of the primary field B is normal to the plane of the plate there is a maximum coupling between the field and the plate. The field and plate are minimally coupled when the primary field is parallel to the plane of the plate. p That is, the maximum amount of flux passes through the plate in this orientation. Conversely, the plate is minimally coupled to the primary field when in orientation 2. When the primary field is terminated, eddy currents are induced i n the plate which give rise to a secondary B-field. 93 Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 94 From Faraday's law, the magnitude to which the eddy currents are induced is proportional to the change i n flux. Therefore i n orientation 1, where there is a greater amount of flux passing through the plate than in orientation 2, there would be a greater eddy current response response. The measured B-field response of an 8 X 8 X 1 inch aluminum plate reflect the orientation dependence of the eddy current response (Figure 4.4(c)). When the primary field is perpendicular to the plane of the plate, the response is far greater than when the primary field is parallel to the plane of the plate. In both orientations 1 and 2, the B-field response of the aluminum plate is flat at earlier times, and decays exponentially at later times. The measured dB/dt A l u m i n u m Plate ( 8 x 8 x 1 inches) Steel Plate ( 8 x 8 x 1 inches) Figure 4.4: Time Decay curves for 8x8x1 inch steel and aluminum plates. Panels (a) and (b) contain the measured <9B/dt response. Panels (c) and (d) contain the B-field response integrated from the measured dB/dt response. Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 95 response of an 8 X 8 X 1 inch aluminum plate is plotted i n figure 4.4(a). When the primary field is perpendicular to the plane of the plate, the time derivative 83/dt is greater than when the primary field is parallel to the plane of the plate. In both orientations 1 and 2, the 83/dt response can be separated into 2 stages. A t early times 83/dt decays as a power law £ / . A t - 1 2 late times the decay is exponential. M a g n e t i c Induction Response The response of magnetic metals, such as steel, exhibits a different decay behavior than the response of non-magnetic metals. When the primary field is on, the magnetic target will become magnetized. The magnetization M is illustrated as as series of dipole moments i n figure 4.5. The extent to which a target will become magnetized is dependent on the shape of the target, Plate Orientation 1 Plate Orientation 2 B B » g gP i n t e n o r > r i m Primary Field On Plate Magnetized Primary Field OffEddy Currents Induced Primary Field On •» Plate Magnetized a y r B interior primary » B Primary Field OffEddy Currents Induced Figure 4.5: Eddy current induction in a magnetically permeable plate. When the direction of the primary field B is normal to the plane of the plate the demagnetization factor is large, and thus the magnetization is small. The demagnetization factor is small when the plane of the plate is parallel to the primary field, and thus the magnetization is large. p and the orientation of the target relative to the primary field. Exact expressions can be written for the magnetization of a uniformly magnetized ellipsoid in a uniform applied field B . The internal B and H fields are p int H jnt B = = P _ R P B + M D o . ( M 1 _ . ( ) 45 D . ) M . • ( 4 Q) Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters where the subscript i denotes the component along the ith principle axis (i — x,y, and z) and Di is the demagnetization factor along the i th constitutive relations B i n t = MoMr H l n t direction (page 33, Brown, 1962). If we assume the , we can solve the above equations for the magnetization M Bf M (Mr - 1) — = — Mi , ^ (4.7) Mo 1 + ^ ( ^ - 1 ) t ' 1 It is clear from equations 4.7 and 4.8 that larger demagnetizing factors will result in a smaller magnetization. Substition of Mi into equation 4.6 results in an expression for the internal field B int The shape dependence of magnetization is contained within the demagnetization factors Di. The demagnetization factors are geometrical factors determined by the axis ratios of the ellipsoid. The three demagnetization factors, corresponding to the three principal axes of the ellipsoid, are positive and constrained by D x +D +D = 1 y (4.9) z If the plate is approximated by an oblate spheroid then, due to the symmetry of the spheroid, there are only two relevant demagnetization factors: 7_?j| which is the demagnetization factor parallel to the axis of symmetry, and D± which is the demagnetization factor along the transverse axes. The two factors are related by D\\ + 2D± = 1. The demagnetization factor for a spheroid in the direction along the axis of symmetry (perpendicular to the plane of the plate) is Du = : v T ^ " " sin 1 f v l - r V* ' J\ ( 4 ' 1 0 ) where r is the ratio of the length along the axis of symmetry (thickness of oblate spheroid) and the width along the transverse axis (diameter of oblate spheroid) (page 89, Bertotti, 1998). Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters As an example, let's consider an oblate spheroid with a length to width ratio of r — 0.125 (e.g. a spheroid plate that is 1 inch thick and 8 inches across). The demagnetization factor along the axis of symmetry would be D\ = 0.831, and the demagnetization factor along the transverse axis would be much smaller; Z> = 0.0846. Substitution of these values into equation 4.7, reveal 2 that i f the plane of the plate is parallel to the primary field (orientation 2), then the induced dipole would be 9.17 times stronger than if the plane of the plate is perpendicular to the primary field (orientation 1). The presence of magnetization increases the B-field inside the target. Recall that when the primary field is terminated, currents are induced on the surface of the target that, by Lenz's Law, try to reproduce the interior field during the transmitter on time. When the interior B-field is larger, then the induced currents will be larger and the secondary field measured at the surface will be larger. Substitution of the demagnetization factors for the example into equation 4.8 reveals that the interior B-field is approximately 9.17 times stronger when the plane of the plate is parallel to the primary field, than when the plate is perpendicular to the primary field. Therefore I expect the response when the plane of the plate is parallel to the field to be larger than when the plate is perpendicular to the primary field. The measured B-field response of an 8 X 8 X 1 inch steel plate reflects the orientation dependence of the magnetic induction response (Figure 4.4(d)). When the primary field is perpendicular to the plane of the plate the B-field response is less than when the primary field is parallel to the plane of the plate. The measured dH/dt response of an 8 x 8 x 1 inch steel plate is plotted i n figure 4.4(b). When the primary field is perpendicular to the plane of the plate the time derivative dB/dt is less than when the primary field is parallel to the plane of the plate for all but the first few time channels. Chapter 4. Relating Physical Parameters of a Metallic Target to Model 4.1.4 98 Parameters T D E M Response of a Rod The measured B-field and dB/dt response curves of a rod are plotted i n figure 4.6. The arguments of the previous section can be used to predict the response of a rod. When the aluminum rod is parallel to the primary field (orientation 1), the rod is minimally coupled with the primary field and the response is less than when the rod is perpendicular to the primary field (Figures 4.6(a) and (c)). A steel rod response will also exhibit magnetic effects. A rod can be approximated by a prolate ellipsoid. The demagnetization factor of a prolate ellipsoid is 1 (4.11) where the aspect ratio r > 1 (page 89, Bertotti, 1998). When a magnetic permeable rod is parallel to the primary field (orientation 1), magnetic effects that increase the secondary response will be at a maximum. When the rod is perpendicular to the primary field (orientation 2), the response will be weaker. For example, let us consider a rod with an aspect ratio r = 8. The demagnetization factor along the symmetry axis of the rod would be D\ = 0.0284. The demagnetization factor along an axis transverse to the symmetry axis is then D 2 = 0.4858. Substitution of D\ and D into equations 4.7 and 4.8 reveal that the induced dipole moment, 2 as well as the interior B-field, i n orientation 1 is approximately 14 times stronger than when i n orientation 2. Figures 4.6(b) and (d) show measured dB/dt and B-field responses, respectively, of an 8 X 1X inch steel rod. For both the dB j dt and B-field the response is much greater when the primary field is parallel to the long axis of the rod. 4.2 Analysis of Laboratory Measurements of Metallic Prisms The decay curve of each axi-symmetric target measured using the experimental setup of Figure 4.1. These targets include both aluminum and steel targets that range from an 8 X 8 x 0.25 Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters Steel Rod (8x1x1 inches) Aluminum Rod (8x1x1 inches) a) 99 1 + 2 O 10 dB dt 10 10 10 10 ' 10 Time (msec) 10 10 10 10 10 10 Time (msec) 10 10 c) Figure 4.6: Time Decay curves for 8x1x1 inch steel and aluminum rods. Panels (a) and (b) contain the measured dB/dt response. Panels (c) and (d) contain the B-field response integrated from the measured dB/dt response. Each set of data is fit with the approximate forward model. inch, plate, to a 8 X 8 X 8 inch cube, to a 8 x 0.25 x 0.25 inch rod. B y measuring the decay curve of each rectangular prism in the two orientations illustrated in figure 4.2, the time decay constants of the approximate forward model can be determined. The approximate forward modelling describes the response of a target as a pair of dipoles that decay independently of each other, with each dipole having its own set of decay parameters. Denote the dipole along the symmetry axis of the target as dipole 1. This dipole is parallel to the major axis for a rod, and perpendicular to the plane of a plate. The transverse dipole is then be labelled dipole 2. Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 100 The two measurement orientations isolate the decay behaviour of each of the two dipoles. This is because the strength of each dipole is proportional to the projection of the primary field onto the dipole direction. Let us consider a plate. When the primary field is perpendicular to the plane of the plate the projection of the primary field onto dipole 2 is zero, thus the approximate forward model assumes the response can be modelled as a single dipole also perpendictdar to the plate. The decay parameters of dipole 1 [k\, a i , ft, and 71) can be estimated by fitting this curve to the decay law (equation 2.37). When the primary field is parallel to the plane of the plate, the response is due to dipole 2 and parameters fc , a , ft, and 7 can be recovered. 2 2 2 In the absence of a rigorous forward modelling code that solves Maxwell's equations for the response of a general metallic object, this set of measurements provided the only means of relating the shape of a target to the parameters of the forward model characterising the response. These relationships were established, firstly, by using a scaled down version of the non-linear least squares techniques outlined i n chapter 3 to obtain the decay parameters k, a, ft and 7 for both dipoles of each of the different targets. A complete list of the recovered parameters for each of the prisms are listed in appendix B . Secondly, I observed how recovered values of model parameters or combinations of parameters changed with the dimensions and magnetic properties of the measured prism. The patterns in the behaviour of the parameters, then led to the shape and permeabibty discrimination diagnostics that are proposed in the following sections. 4.3 Relating f3 to Magnetic Permeability U X O are generally made of steel, which is a ferrous material. Therefore, the magnetic permeability is most likely an identifying characteristic of U X O . To draw link between magnetic permeability and model parameters, I first performed forward modellings for a series of spheres varying in size and permeability. The analytic solution for the time domain response of a Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 101 a) Figure 4.7: The behaviour of parameter (3 for various size spheres with varying permeability p. Panel (a) contains results of fitting the B-field with the parameterization. Panel (b) contains results of fitting the dB/dt with the parameterization. Each set of data is fit with the approximate forward model. sphere in a uniform primary field was outlined in Chapter 2. The suite of data from these modellings were then fit with the time decay law to obtain best fit values of the decay parameters. Relationships between the time parameters and the magnetic permeabilities were then considered. Figures 4.7(a) and (b) suggest that the value of (3 obtained for a sphere may be diagnostic in determining whether the sphere is permeable or non-permeable. Both plots suggest setting a threshold value of /3, such that targets with recovered (3 values greater than the threshold are most likely permeable. This diagnostic is examined for both B-field and dH/dt data. 4.3.1 V a r i a t i o n s o f f3 w i t h M a g n e t i c P e r m e a b i l i t y w h e n U s i n g B - f i e l d D a t a B y applying the parameterization to the B-field of a steel sphere (p = 150^ ), we see that, o for spheres with radius between 5 to 15 cm, (3 falls between 0.4 and 0.5. We also see that a non-permeable sphere (// = p ) has a j3 value of approximately 0.1. Therefore, when applying Q the parameterization to the B-field response of a sphere, a threshold value of 0.3 is suitable. A recovered (3 greater than 0.3 would likely be permeable, while a recovered f3 less than 0.3 would Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 102 likely be non-permeable. The use of (3 as a diagnostic to determine permeability can be extended to non-spherical targets by looking at the recovered beta values for the aluminum and steel prisms. The model uses two values of (3, one for each of the excited dipoles, to describe a buried target. I suggest taking the average of the two recovered (3 values. I will label this average as (3. The 9 axisymmetric aluminum targets had average (3 value of 0.17 with a standard deviation of 0.03. The 11 steel targets had an average of 0.5 with a standard deviation of 0.2. These averages fall on either side of the 0.3 threshold obtained by fitting sphere B-field responses. 4.3.2 Variations of (3 with Magnetic Permeability when Using dB/dt Data B y applying the parameterization to dB/dt response of a steel sphere (/x = 150/i ), we see that o spheres with radii between 5 to 15 cm that (3 falls between 1.11 and 1.35. We also see that a non-permeable sphere (/x = p ) has a (3 value of approximately 0.5. This value of (3 corresponds Q to the early time i - 1 / behavior that Kaufman (1994) predicted for a non-permeable sphere (see 2 equation 2.33). Therefore, when applying the parameterization to the time derivative of the field, a value of (3 greater than about 0.5 indicated that the target was most likely permeable. 4.4 The Ratio k /k l 2 As we saw i n section 4.1.3, that for both the dB/dt and B-field response is greater for a steel plate when the primary field is parallel to the plane of the plate. A stronger induced dipole is reflected in a larger k value. Thus, for a steel plate the fc-ratio k\/k 2 < 1- For a steel bar the response is greater when the primary field is parallel to the bar (and thus along the axis of symmetry). Thus for a steel rod the fc-ratio k\/k 2 > 1. Figure 4.8 shows how the recovered fc-ratio varies for targets ranging from a steel plate to a steel rod. We also saw in section 4.1.3 that for both the B-field and dB/dt response of an aluminum plate is larger when dipole 1 is excited than when only dipole 2 is excited. Thus for an aluminum Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 103 plate the fc-ratio k\jk > 1- The opposite orientation effect was observed for an aluminiim rod, 2 and thus k /k 1 2 < 1 (Figure 4.9). Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 104 Steel -- Fitting B-field Block Dimensions (inches) c) Block Dimensions (inches) d) Steel - Fitting 5B Steel - Fitting at at k volume — © — Dipole 1 —A-- Dipole 2 plate ^ ' ffl ' i cube / rod r \ k 1 / k 2 i^--"^"?^^ i i l 1 1 1 1 1 Block Dimensions (inches) 6 £ Block Dimensions (inches) F i g u r e 4.8: Relating the aspect ratio of a steel target with the ratio ki/k2- " 2 Plot (a) contains the recovered k parameter when fitting the parameterization to the B-field response obtained by integrating the dB/dt response of a steel target. Plot (b) illustrates the relationship between the fci/&2 ratio and the shape of a steel target. Plot (c) contains the recovered k parameter from fitting the parameterization to the measured dB/dt response. shape of a steel target. Plot (d) illustrates the relationship between the ki/ko ratio and the 105 Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters a ) Aluminum - Fitting B-field b) Aluminum -- Fitting B-field n cub k,/k volume a CD _ • CO fi'-' CO ro 1 rod 2 10° A CO ' n ! 10' CO ^ CO *. PJ ro ^ Block Dimensions (inches) Aluminum - Fitting 1 1 1 1 . .. 1 •ftft i x s Block Dimensions (inches) d) Aluminum - Fitting aB at Block Dimensions (inches) Block Dimensions (inches) F i g u r e 4.9: Relating the aspect ratio of an aluminum target with the ratio ki/k2- Plot (a) contains the recovered k parameter when fitting the parameterization to the B-field response obtained by integrating the dB/dt response of an aluminum target. Plot (b) illustrates the relationship between the ki/k2 ratio and the shape of an aluminum target. parameterization to the measured Plot (c) contains the recovered k parameter from fitting the dB/dt response. ki/k2 ratio and the shape of a aluminum target. Plot (d) illustrates the relationship between the Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 4.5 The Ratio /V/3 2 In addition to the relative strength of the dipoles being shape dependent, the slope of the time decay response (either d~B/dt or B-field) during the intermediate time stage will be greatly dependent upon the target shape. This effect was seen in steel targets only. In Chapter 2, I showed that the steepness of the response during the intermediate time stage is reflected in the parameter {3 (Figure 2.11). A dipole that decays at a greater rate will have a larger f3. In section 4.1.3 it was seen that the rate of decay of both the B-field and dB/dt response is greater when the plane of a steel plate is perpendicular to the primary field (dipole 1), than when the plane of a steel plate is parallel to the primary field (dipole 2). Thus, for a steel plate the /3-ratiofli/p > 1- In the case of a rod, both the B-field and dB/dt response decay 2 faster (and thus /? is larger) when the main axis of the rod is perpendicular to the primary field (dipole 2). In the case of a steel rod the /3-ratio P\/02 < 1 (Figure 4.10). For aluminum targets the response shape look essentially the same for each of the targets. That is, the B-field response is essentially flat at early time and exponential at later times. The dB/dt response exhibits a power law decay o f t / - 1 2 and is exponential at later times. The decay curves for aluminum tar gets are essentially the same regardless of target shape. Therefore there is no relationship between the /3-ratio and the aspect ratio (Figure 4.11). 107 Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters a) Steel -- Fitting B-field Steel - Fitting B-field b) L plate /^ ' \^ -'^ i NSS \^n ' o i- rod _ 1 —-a - Block Dimensions (inches) Block Dimensions (inches) Steel - Fitting d) ISL B,/P cub p ' dB 8t 0 ! ' 1 rod " 2 - Block Dimensions (inches) Block Dimensions (inches) F i g u r e 4.10: Relating the aspect ratio of a steel target with the ratio Pi/fa. Plot (a) contains the recovered j3 parameter when fitting the parameterization to the B-field response obtained by integrating the dB/dt response of a steel target. Plot (b) illustrates the relationship between the /3i//?2 ratio and the shape of a steel target. Plot (c) contains the recovered (3 parameter from fitting the parameterization to the measured dB/dt response. Plot (d) illustrates the relationship between the /3i//?2 ratio and the shape of a steel target. 108 Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters a Aluminum -- Fitting B-field ^ b ) Aluminum - Fitting B-field ^ 0.5 0.45 — e — Dipole 1 — A — Dipole 2 ©-^ cub 'I B ' • rod f 0.4 0.35 iyp 0.3 2 0.25 - 0.2 / 0.15 0.1 0.05 Co CO Block Dimensions (inches) Block Dimensions (inches) ) Aluminum - Fitting £B dt Aluminum - Fitting 5B at d) plate 0.6 ^ "'^ 5 p ' B ! cub rod J 0.55 0.5 / ' A P /P 1 // 0.45 2 0.4 ' /' — e — 0.35 • / Dipole 1 Dipole 2 Block Dimensions (inches) Block Dimensions (inches) F i g u r e 4.11: Relating the aspect ratio of a aluminum target with the ratio Pi/fy- Plot ( ) contains the recovered (3 parameter when fitting the parameterization to the B-field response obtained by integrating the dB/dt response of an aluminum target. Plot (b) illustrates the relationship between the fii/fo ratio and the shape of an aluminum target. Plot (c) contains the recovered (3 parameter from fitting the parameterization to the measured dB/dt response. Plot (d) illustrates the relationship between the (3i/P2 ratio and the shape of a aluminum target. a Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 109 Figure 4.12: The T D E M response of an 81 m m mortar. Plot (a) contains the measured dB/dt response, and plot (b) contains the integrated B-field response. 4.6 Analysis o f Measured U X O The previous sections suggest values that may be diagnostic in determining the magnetic permeability of the target, as well as the basic shape of the target. However, these tests were performed on solid aluminum and steel rectangular prisms that may, or may not be, representative of real field UXO. As we saw in figure 1.1, UXO can be found that are longer or wider than 8 inches. Also UXO are generally made of steel, but may have some components (such as stabilizer fins) that consist of non-permeable material. In addition, UXO are clearly not solid objects and can possess complicated inner structure. It is not clear how these factors affect the measured T D E M response, and not obvious if we can apply the diagnostics of the previous sections. In the following analysis, the decay curves of each of the 24 UXO measured at Geonics will be inverted for the decay parameters. The decay curves of two of the UXO, an 81 mm mortar and a 155mm shell, are shown in figures 4.12 and 4.13. As with the prism data, the B-field response was obtained by integration. The basic characteristics of a permeable rod that we observed in the previous section are evident in these UXO. That is, dipole 1 is generally stronger and decays more slowly than dipole 2. Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters a) F i g u r e 4.13: The T D E M response of a 155 m m mortar. Plot (a) contains the measured dB/dt response, and plot (b) contains the integrated B-field response. The first diagnostic considered was to use (3 parameters to predict whether the target was permeable or non-permeable. Figure 4.14 contains the results from taking the average of /3\ and /?2- We see that for each U X O the recovered /3 value when applying the parameterization to the 83/dt response exceeds the threshold value of 0.8 that was suggested earlier. When using the B-field response, we observe that the appropriate threshold value is exceeded for 22 of the 24 U X O . We can then look at calculated k and (3 ratios for each of the U X O (Figures 4.15 and 4.16). Several inconsistencies are evident when comparing these results with those obtained from metallic prism measurements. Let us consider figure 4.15. There are two features that were not found in the measurements of solid metallic prisms. Firstly, we see that for U X O number 20 (diameter = 10.8 cm, length = 65 cm) and U X O 21 (diameter = 12.8 cm, length = 85 cm) that both k and (3 ratios are less than 1. Secondly, the k and j3 ratios U X O 12 are opposite to what is expected for a permeable rod-like target. The /? values for U X O 12 exceeds the (3 threshold when performing the analysis on either 83/dt or B-field. It is impossible to know what is the source of these inconsistencies without knowing more about the structure of the U X O , and without having a forward modelling code with which to generate the decay curves. Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters dB/dt p threshold B-field p threshold ! 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 UXO Sample # F i g u r e 4.14: /3 averages for U X O . When fitting the dB/dt response of a U X O with the parameterization, the average of the recovered /3's + ( 3 2 ) / 2 = (3) is greater than 0.8. When fitting the B-field response of a U X O with the parameterization, the average of the recovered /3's is greater than 0.3 (with the exception of U X O sample 21). F o r the r e m a i n i n g 21 U X O results are as expected for a permeable r o d - l i k e target. F i g u r e 4.16 contains results f r o m a p p l y i n g the parameterization to the B-field. M o r e inconsistencies occur when dealing w i t h the B-field. 4.7 Proposed Algorithm For Determining Dimensionality of Target Several diagnostics were i n t r o d u c e d i n this chapter that give hints to m a g n e t i c p e r m e a b i l i t y of a target a n d the r o u g h shape of a target. These can be c o m b i n e d as a n a l g o r i t h m for d e t e r m i n i n g the dimensionality. F i r s t l y a decision should be made as t o w h i c h d a t a , B-field or dB/dt, should be inverted for parameters. T h r o u g h o u t this chapter the decay l a w has been a p p l i e d to b o t h the B-field a n d dB/dt response of the data. For the m e t a l l i c p r i s m analysis b o t h forms of the d a t a p r o d u c e d essentially the same diagnostics. T h a t is, i f (3 is greater t h a n some threshold value then the target is l i k e l y permeable, and the k and j3 ratios follow the same patterns for Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 100.00 2.60 2.20 + 10.00 1.80 k\lki Pl/p2 1.40 1.00 1.00 0.60 L _ 0.10 0.20 UXO Sample # Figure 4.15: ki/k and (3i/(3 ratios when fitting to derivative offielddB/dt. The fli/(3 ratio is plotted with bars, and the line is the plotted ki/k ratio for each of the UXO. Recall that for a permeable rod-like prisms that ki/k > 1 and /3i/(3 < 1. 2 2 2 2 2 2 shape. Although this is the case, it is probably more robust to invert for the derivative of the field. Indeed, when examining the use of the diagnostics to sample U X O , fitting the dH/dt data produced more desirable results. There are errors introduced when integrating the finite number of dH/dt data, including having to extrapolate a function in order to approximate late time behavior. I conclude this chapter with a possible algorithm to determine the dimensionality of a target. The Algorithm For B-field Data 1. Recover model parameters by using a non-linear parameterization on the B-field data 2. Take average of ft and ft to obtain ft If (3 is greater than 0.30 then the target is most likely permeable. 3. Compute the ratios ft/ft and k\/k . There are two options: 2 • /3 > 0.3: If k\/k > 1 and ft/ft < 1 then a permeable rod-like target was measured. 2 112 Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters 113 100.00 2.60 2.20 10.00 1.80 hlki Pl/p2 1.40 1.00 1.00 0.60 + 0.20 T - ( N r t ^ l / ) ( D S C O O ) 0 ' - I N r t ' J i n ( D N C O O ) C -4- 0.10 UXO Sample # Figure 4.16: Recovered k and 0 ratios when applying parameterization to the B-field response. The Pi 102 ratio is plotted with bars, and the line is the plotted ki/k ratio for each of the U X O . Recall that for a permeable rod-like prisms that ki/k > 1 and 0i/f3 < 1. 2 2 2 If ki/k < 1 and 3i/8 > 1 then a permeable plate-like target was measured. 2 * P < 0.3: If k\/k 2 2 > 1 then non-permeable plate-like target was measured. If k\/k < 2 1 then the target is rod-like. Pi/P does not give supporting, or extra, information. 2 The same algorithm for the dB/dt field uses different threshold values for 0: The Algorithm For dB/dt Data 1. Recover model parameters by using a non-linear parameterization on the dB/dt data 2. Take average of Q\ and Q to obtain ji. If/3 is greater than 0.80 then the target is most 2 likely permeable. 3. Compute the ratios P1/P2 and ki/k - There are two options: 2 • P > 0.8: If k\/k2 > 1 and P1/P2 < 1 then a permeable rod-like target was measured. H ki/k < 1 and P1/P2 > 1 then a permeable plate-like target was measured. 2 Chapter 4. Relating Physical Parameters of a Metallic Target to Model Parameters • /3~ < 0.8: If k /k 1 2 114 > 1 then non-permeable plate-like target was measured. If k /k 1 2 1 then the target is rod-like. {3\jf3 does not give supporting, or extra, information. 2 < Chapter 5 Field Application The previous chapters of this thesis outlined the development of a data interpretation strategy for the location and identification of U X O . A n approximate forward modelling was developed that predicted a metallic object's secondary field response, due to a step off primary field, as a pair of decaying orthogonal magnetic dipoles. A model vector m was defined whose 13 elements included the target's location, orientation, and characteristic decay parameters. After making the assumption that the approximate forward modelling could adequately replicate the response of a metallic target, I developed a non-linear parametric inversion to determine the model m that minimizes the difference between the observed data set and a data set predicted by the approximate forward modelling. Finally, relationships between the relative values of the decay parameters and real physical parameters of a target were made based on lab T D E M measurements of a series of different sized metallic prisms. Ratios of the 3 and k parameters are possible indicators of the aspect ratio and magnetic permeability of the object. U X O are generally rod-like and magnetically permeable, so these relationships then can then form the basis of a U X O identification algorithm. To this point i n the thesis, this algorithm has been tested on synthetic data sets only. In this chapter, the inversion methodology and U X O identification diagnostics are applied to E M 6 1 3D secondary field measurements taken over a buried 105 m m shell. From this application, the algorithm's performance will be tested in the presence of real world natural and cultural noise sources, as well as i n the presence of errors resulting from adopting an inaccurate forward modelling. 115 Chapter 5. Field 5.1 116 Application Description Field Data Set 5.1.1 Target Information Geonics L t d . conducted a T D E M survey over a 105 m m dummy artillery shell buried at the York University Environmental Test Site. A 105 m m shell is approximately 40 cm long with a diameter of 10.5 cm, and therefore has an aspect ratio of 4 to 1. A photograph taken by Geonics Ltd. of the shell can be found in Figure 5.1(a). Although the dimensions of the buried shell were not measured prior to burial, I assumed that the buried shell was similar in shape and size to the photographed shell, because both artillery shells were supplied by the Canadian military. The target is both rod-like and magnetically permeable, and therefore possesses the features that I use to characterize U X O . The buried shell is part of the York University Environmental Test Site, and therefore the shell location, depth, and dip are relatively well characterized (Szeto, 1996). The shell was buried at a depth of 50 cm below the surface of the earth, and with a dip of 45 degrees (c = 45 °) (Figure 5.2). The strike (a) was not recorded upon burial of the shell. 5.1.2 Instrumentation The survey was carried out using the Geonics EM61-3D (Figure 1.2). The EM61-3D is a prototype system that records the 3 components of secondary field due to a pulse excitation. The primary field is produced by a 1 m square transmitter loop carrying a 50% duty cycle bipolar current. The transmitter loop generates a dipole moment of 64 A m . 2 Three orthogonal circular receiver loops located at the center of the transmitter loop measure the secondary field. Receiver coils with diameters of 60, 58, and 63 cm measure the x, y, and z components of the field. The fields sensed by the receiver coils are processed and stored in a P R O T E M 4 7 time domain receiver. The P R O T E M 4 7 receiver was configured to take measurements in 30 geometrically spaced time gates. The input is integrated over the length 117 Chapter 5. Field Application dB/dt B-field k a P 7 Orient. 1 Orient. 2 0.17 0.21 0.25 24.0 0.031 0.070 0.52 30.7 Orient. 1 36.1 0.010 0.76 28.7 Orient. 2 14.7 0.025 1.33 50.7 F i g u r e 5.1: (a) Photograph of a 105 m m artillery shell (courtesy of Geonics Ltd.). The shell measures 10.5 cm i n diameter and 40 cm long, and weighs 15kg. (b) Laboratory measurements of the secondary B-field of a 105 m m artillery shell, (c) Laboratory measurements of the secondary dB/dt of a 105 m m artillery shell. The responses in (b) and (c) are fit with the decay law L(t) = k(a +1)~^ exp-t/j. (d) The decay parameters from the fitted curves for the B-field and dB/dt data. 118 Chapter 5. Field Application a) 400 b) 100 cm 350 a E M 6 1 - Transmitter Loop Loop Height = 43 cm 3009 250 » X(cm) 2 0 0 & 150o Depth = 50 cm 100s 500 105 mmV V 0 50 100 150 200 250 300 350 400 Y(cm) Figure 5.2: (a) The 105 mm shell buried as part of the York University Test Site was buried at a depth of 50 cm. The transmitter loop is mounted on a trailer, and therefore sits parallel to the ground at a height of 43 cm. (b) Line and station layout of the the field data survey grid. The component data were collected on 9 lines, with a line spacing of 50 cm. A total of 17 stations were spotted on the line, with a station separation of 25 cm. The coordinate system was set up such that the Z direction points into the ground. of the time channel. The calculated mean value within the integration window is then recorded as the measurement at the center of the time window. The measurement is repeated a number of times and stacked to increase the signal to noise ratio. The system is mounted onto a trailer that is pulled along by the instrument operator. The transmitter loop is mounted parallel to the ground and at a height of 43 cm above the ground, and the P R O T E M 4 7 receiver is located in a backpack carried by the operator. 5.1.3 Survey Geometry EM61-3D measurements were taken over a 4mx4m square survey area centered on the ordnance location. The survey area was covered with total of 9 lines of data, with 50 cm spacing between the lines (Figure 5.2). Along each line there was a total of 17 stations, with a station separation of 25 cm. A t each station three components of dB/dt data were measured over 30 time channels. Figure 5.3 contain the vertical component, 0B /dt, of the time derivative of the secondary z B-field measured during the odd time channels. Figure 5.4 contain the horizontal component, 119 Chapter 5. Field Application (^JTT) + J 2 ' °f ti m e derivative of the secondary B-field measured during the odd time channels. The 'o' markers on the plot indicate locations of stations. 5.2 Characterization of EM61-3D Errors for Inversion In Chapter 3, the approximate forward model was expressed as d = F [m], j j = 1,2,3, ...N j (5.1) This equation expresses the mapping of the model vector m to a datum dj by a function Fj. The forward mapping Fj is defined by equation 2.38. This forward mapping is inexact because all data are noisy and the forward modelling is approximate. A more accurate representation of the situation encountered in field data is achieved by rewriting equation 5.1 as j = 1,2,3, ...N, dj = Fj[m] + ej, where €j is the error on the j th (5.2) datum. Because the goal of the inversion is to find a model that minimizes the misfit between the observed data and the data predicted by the forward modelling, the algorithm will try to fit the errors €j and therefore affect the recovered parameters. The source of the errors €j can generally be categorized as being either natural, cultural, or modelling errors. A discussion of the sources of error in electromagnetic measurements can be found in numerous papers (for example see (Nabighian and Macnae, 1991)). Sources of noise in T D E M measurements include ughtning discharges (sferics), power-line noise, amplitudemodulated ( A M ) long wave and very low frequency ( V L F ) transmitters. Modelling errors include any discrepancy between the approximate forward mapping F[m] and an exact forward model. Although this discrepancy may have a large effect on the values of recovered model parameters, these differences are not discussed here because their investigation requires access to an exact forward modelling code. 120 Chapter 5. Field Application Time = 0.036msec Time = 0.057msec Time = 0.148msec Time = 0.234msec Time = 0.092msec 11 Cr Time = 0.525msec Time = 0.8025msec Time = 1.998msec Time = 3.198msec Time = 13.49msec Figure 5.3: Vertical component of dB/dt EM61-3D data taken over a 105 mm artillery shell at the York University Geophysical Test Site. 121 Chapter 5. Field Application n Time = 0.057msec Time = 0.092msec Time = 0.148msec Time = 0.234msec Time = 0.3525msec Time = 0.525msec Time = 0.8025msec Time = 0.036msec 0 % 4] 2 0' ' ' ' 1raf -ml: 4 Time = 3.198msec Time = 1.998msec 0' Time = 13.49msec Time - 21.9msec Figure 5.4: Horizontal component of dB/dt EM61-3D data taken over a 105 mm artillery shell at the York University Geophysical Test Site. 122 Chapter 5. Field Application In Chapter 3, the noise was dealt with in the objective function. The objective function was denned to be <f> = \\W (F[m] - d° >) || b (5.3) 2 d where W d is denned such that WW d errors are independent, then W d d is the inverse of the covariance matrix of the data. If the is a diagonal matrix with elements equal to the inverse of the standard deviation of the error, and the objective function can be written i=i j where N is the number of data. When the errors are Gaussian with zero mean, then the sum in equation 5.4 becomes the x 2 measure. % is a statistical variable with an expected value of 2 approximately JV. As was reported in Chapter 3, this weighting of the objective function is optimal in the sense that the recovered model will be the one whose variance is smallest (Bard, 1974). In problems where the forward modelling is inaccurate, the presence of modelling error necessitate that the objective function be fit with x 2 > N- Although the above outlined method of dealing with data errors is reasonable (and in linear cases optimal in the sense of finding a least variance estimate), an estimate of the standard deviation of the error is required prior to the inversion. Effers0, et al. (1999) examined the effect of A M and V L F transmitters on T D E M measurements. They observed for their particular T D E M instrument that the standard deviation of the dB/dt signal exhibits a i _ 1 proportional- ity when A M transmitter noise is log-gated and stacked. Munkholm and Auken (1996) showed that log-gated and stacked white noise maps onto the T D E M response as errors with a standard deviation exhibiting a r - 1 / 2 decay. A decision must be made as to what error level should be assigned to the data in this inversion. The investigations into the nature of T D E M noise cited above suggests that for each time channel an error with constant standard deviation be assigned to all the data measured in 123 Chapter 5. Field Application that time channel. I will adopt this form of noise model for inversions of the field data. In this section I consider two methods of specifying the time variation of the noise. In the first method, inversions using data from single time channels are performed using a two-dipole model. Each inversion determines the minimum misfit for each time channel and that is used to estimate the noise for that time channel. In the second method, I assume that the measurements taken far from the target are of noise only and that these measurements can be used to characterize the noise. Based on these measurements a noise model is generated by fitting a time function e (t) = At~ B 5.2.1 to the measured noise. Method 1: Estimation of Noise By Fitting Subsets of the Complete Data Set In order to determine the standard deviation of the error for each time channel, the set of data during each time slice will be inverted separately to find the best fitting two dipole model. Because this is a subset' of the data, the fit to each data point will be better than when fitting the data from a number of time channels. Indeed, the misfit obtained by inverting the single time channel represents a minimum misfit of that time channel's data that could be hoped for when inverting all data from multiple time channels. The goal is to minimize an objective function N <p(U) = \\F[xn,U] ~ d obs 2 (r,)|| = 2 £ (^[m.fc] - df'{U)) (5.5) Here ti represents the time channel from which the data to be inverted are taken. There are N data in this inversion, where N = 3X (number of locations). The vector m is the vector of model parameters representing two dipoles: m = [X,Y, Z, a, c, ki, k ] 2 (5.6) where (X, Y) represents the location on the survey, Z is the depth, a and c specify the orientation, and k\ and k are the strengths of the two dipoles. As in the case with the approximate 2 124 Chapter 5. Field Application forward model, dipole 1 is parallel to the main axis of the target and has a strength that is proportional to the projection of the primary field in this direction. Dipole 2 is specified to be perpendicular to dipole 1. The model vector m is obtained by using the non-linear parameter estimation procedure of Chapter 3. A n estimate of the standard deviation of the error for time channel ti can then be obtained by considering the \ measure at the solution: 2 where £j (U) is the standard deviation of the j th data point. The expected value for x I 2 s approximately equal to the number of data N (when JV > 5). Since I assume that the error has a constant standard deviation in each time channel, the equation 5.7 can be rewritten ( - ) 5 ~ N 8 This inversion is repeated on the data from each time channel, and the estimated error in each time channel defines the time decay of the error, e(r). Test o n a S y n t h e t i c D a t a Set This method was first tested on a synthetic data set. Synthetic data from a 2 m x 2 m survey area containing 5 lines and 5 stations per lines were generated. A l l three components were used in the inversion, making a total of 75 data per time channel. A uniform error was set to be 1% of the maximum 5 B / d r . measurement in each time channels for the first test, and z then set to 5% for a second test. Each time channel was then fit with the best two dipole model, characterized by the 7 parameters X, Y, Z, a, c, k\, and k . Figure 5.5 shows how well 2 the data fitting worked for the time slice t = 1.788msec. Error with a standard deviation of 5% of the maximum value of the ^-component was added to the data. In the left column of Figure 5.5 are the synthetic observed data, in the middle column are the predicted data, and 125 Chapter 5. Field Application Observed Data Predicted Data d p r e d -d' n 12 _ Figure 5.5: Observed and predicted data when inverting dH/dt data from time channel 2=1.788 msec. The random appearance in the plot of the difference in the predicted and observed data reflect the random error assigned to the data of this time channel. in the furthest right column is the difference between the observed and predicted data. The random error assigned throughout the data is reflected in the random looking character of the difference plot. Figure 5.6 contains results of fitting each of the time channels in order to estimate the standard deviation of the error for that time channel. The solid plotted lines indicate the true standard deviation of the error for the time channel, and the dashed lines represent the error estimated by fitting each time channel and using equation 5.8. In both the 5% and 1% cases, the technique was able to estimate the error. 126 Chapter 5. Field Application Time (msec) Figure 5.6: Noise model recovered in two synthetic test cases. In the first case, the standard deviation of the assigned noise was set to be a 5% of the maximum value of the ^-component of the data of that time channel. In the second case the standard deviation of the assigned noise was set to 1% of the maximum. Solid lines represent the noise assigned in each time channel, and dashed lines represent the estimated noise by inverting each time channel for a pair of dipoles. The technique was successful in recovering the noise. Application to the Field Data Set The proposed error estimation technique can be applied to the field data set. For each of the 30 time channels, 3 components of dB/dt data at 153 stations for a total of 459 data, are inverted for the 7 parameter model vector of equation 5.6 representing the two dipole model. The major difference between this field data set and the synthetic data set to keep in mind is that an inaccurate forward model is being used in this study. Figure 5.7 contains the observed and predicted data, and their difference for the inversion of the first time channel data (r=0.036msec). The residual map (the difference between the predicted and observed data) indicates that the Z-component of the data is fitting the signal, but not the random noise. However, in the residual maps for the X and Y components there is correlated signal. This inversion was repeated for each of the remaining time channels, Chapter 5. Field Application Observed Data 127 Predicted Data *10~ ,! d pred_ obs d n5 x 1 0 - i—i Figure 5.7: Inversion of t = 0.036msec time channel dB/dt data acquired over a 105 mm UXO buried at the York University Environmental Test Site. and the predicted error from each time channel is plotted in Figure 5.8. Also included, for reference, in this plot are decay curves from measurements taken directly above the U X O ((X, Y) — (200,200) cm) and also along the line X — 0 cm located at the edge of the survey grid. The data along the X = 0 cm line are very noisy, fluctuating from negative ( V markers) to positive ('+' markers). The noise model (labelled 'Method 1 Noise Model' on the plot, is of the same order of magnitude as the data along the X = 0 cm line at early and late times. If this noise model were adopted for this data set, most of the measurements along X = 0 cm seem to be essentially noise. In the inversion analysis, these data should either be rejected, because any useful signal is being swamped by noise, or the measurements should be weighted in the data objective function such that they have close to zero influence on the inversion. Both options Chapter 5. Field 128 Application time (msec) Figure 5.8: Noise model recovered for the field data by fitting a pair of dipoles to each time channel. The recovered noise model is plotted, along with the Z-component dB/dt response measured directly above the UXO ((X, Y) = (200, 200) cm), and several (noisy) decay curves measured along the X = 0 cm line. will be tried later in this chapter. 5.2.2 Method 2: Estimation of Noise By Using Measurements at a Distance as Noise Measurements A second method of determining a suitable noise model is to assume that, when the transmitter/receiver is far away from the target, the only signal measured by the instrument is noise. Figure 5.8 contains measurements of the Z-component of the 63/dt field at different locations on the survey grid. In that plot, measurements at a few of the stations along the X = 0 cm line are plotted (thin plotting lines), as well as decay curves from station (X,Y) = (200, 200) cm (thick plotting lines). The decay curves measured along the line X = 0 cm fluctuate wildly, oscillating between positive values (represented by markers) and negative values (represented by ' o ' markers). It is reasonable to assume that these measurements can be used to characterize noise because the data has the erratic behaviour of noisy data, and because the vertical Chapter 5. Field Application 129 component of secondary field produced by a buried target does not have zero cross-overs. The bipolar nature of the data must therefore be created by some superimposed noise signal. I adopt a form of a noise model proposed by Munkholm and Auken (1996, equation (2)) Munkholm and Auken found that for log-gated, gated stacked data that white noise is seen as a r - 1 / 2 signal, and Effers0et. al. (1999) found that A M transmitter noise maps to a t~ l signal. In order to find the appropriate constants A and B for the field data survey, I first took measurements from stations at the extremities of the survey area and fit a best straight line (in a plot of log(|5B^/5t|) vs. log(r)) for the decay measurement at each station. Stations on lines X = 0 cm and X = 400 cm, were used for this purpose. The value of B is set to be the average slope of the fitted straight lines, and A is set such that all the noise falls beneath the line. The noise model resulting from this process is e(t) = (252.5) r - ° - . (5.10) 9 4 A l l data that fall below this line are considered to be essentially noise, with all meaningful signal from the target being swamped. Figure 5.9(a) plots the data used for determining the noise model, and the final noise model given is equation 5.10. From this figure it is clear that assigning noise to the data according to equation 5.10 results in a higher percent error being assigned to earlier time channels and late time channels, with a lower percent error being assigned to intermediate time channels. For some stations the response at late and early times is swamped by noise, but at intermediate time there is still some data above the noise line. The decay curve measured at (X,Y) = (350,150) cm is a clear example (Figure 5.10(a)). As expected, the percent error assigned to data further from the target will be higher since the signal to noise ratio will clearly be smaller. Figure 5.10 gives an example of how data away from the target falls below the noise level. 130 Chapter 5. Field Application Figure 5.9: Characterizing the noise of the EM61-3D data with e(t) = At . (a) Noise model obtained using method 2. The noise model, e(t) — 252.2f , was obtained by fitting the decay curves on lines X = 0 and X = 400 cm, which are plotted as 'o' and '+' symbols representing negative and positive values. The decay curve measured at (X, Y) = (350,150) cm is swamped by noise at early and late times, but appears to have signal at intermediate times, (b) The signals measured on lines X = 0 and X = 400 cm are plotted here to demonstrate that the mean of the noise is approximately zero, and that the standard deviation of the noise is bounded by the noise model. B -094 5.3 Preliminary Data Analysis In this section the methods of determining a reasonable starting guess for the inversion algorithm, as outlined in Chapter 3, are applied to the field data set. 5.3.1 Estimates of Location and Orientation Obtained by Plotting the Data The following estimates of strike, depth, and location are obtained using techniques outlined in section 3.5.1. Strike a The line of symmetry when plotting the horizontal component of the response generally indicates the strike of the target. Figures 5.11(a) and (b) indicate that it is not easy to discern a single, distinct line of symmetry for the data. The line of symmetry could be chosen as being either parallel to the lines or perpendicular to the lines. For the following inversions, the initial guess Chapter 5. Field Application ( ) a Time = 0.036msec Time = 0.525msec Time = 17.19msec Figure 5.10: Removal of data that fall beneath the noise level. In the top set of panels (a), are the Z-component of the dB/dt data, for three time channels. The bottom set of panels (b), have those data which fall beneath the noise level set by the e(t) = At~ B removed. for the strike will be set to a = 45 degrees. Depth Z In section 3.5, it was shown that the depth of the UXO could be estimated by taking the width of the z-component of the response at half the maximum of the signal. The width of the signal was to be taken along the strike of the target. However, the strike of the target is not easily determined from the plotted data. In figure 5.12(a) the z-component of the B-field is plotted for the tenth time channel. Two candidates for the strike of the U X O are indicated by white arrows. The profile of the signal along the line X = 200 cm is plotted in figure 5.12(b) and the profile of the signal along the stations Y = 200 cm is plotted in figure 5.12(c). Applying the half-width estimation technique to the profiles in (b) and (c) result in depth from surface estimates of 76.4 cm and 152.5 cm. Chapter 5. Field Application 132 Figure 5.11: The strike of the buried UXO are difficult to establish from plots of the horizontal component of the field. A line of symmetry for the above two plots can either be drawn parallel to the data lines, or perpendicular to the data lines. L o c a t i o n (X, Y) The location of a target is estimated to be directly below the peak i n the vertical (z) component of the response. Referring to figure 5.13, the maximum of the response is measured at (X,Y) = (200, 200) cm. The location of the maximum remains is static throughout the 30 time channels. 5.3.2 F i t t i n g A S i n g l e D i p o l e at E a c h T i m e C h a n n e l A second method of obtaining estimates of depth and orientation was to compute parameters for a best fitting single dipole to data from a single time channel (Section 3.5.2). The dipole is defined to be oriented along the axis of the U X O , and the dipole strength is proportional to the projection of the primary field along that direction. The data in each time slice was inverted for a location ( X , Y ) , a depth (Z), orientation (a and c), and dipole strength (k). The surficial coverage of data for this inversion was a subset of the complete survey grid. 5 lines (X = 100,150,200,250, and 300 cm) with 25 cm station spacing along each line were inverted for each time slice. This inversion was performed on both B-field and 83/dt data. Inversions using all three components of data and also with just the Z component of the data were carried out. The recovered model parameters obtained for some of the time channels by inverting all Chapter 5. Field Application 0 Y (cm) 133 1 2 Y(m) 3 4 X (cm) Figure 5.12: Determining the depth of the buried UXO. Panel (a) contains a plot of the vertical component of the B-field measured during the tenth time channel. The two white lines indicates where profiles of the signal have been taken for depth estimation. Panel (b) contains a profile taken along the Line X = 200 cm. The depth estimate (from the ground surface to the center of the target) using this profile is 76.4 cm. Panel (c) contains a profile taken perpendicular to the lines on stations Y = 200 cm. The depth estimate using this profile is 152.5 cm. Chapter 5. Field Application 134 Figure 5.13: Estimating the location of the buried UXO. B-field measured during (a) time channel 1 and (b) time channel 25. The maximum measured B-field response is located at (X, Y) = (200, 200)cm. three components of dB/dt data are listed in table 5.1. Listed on the row labelled 'Starting Model' are averages from the time channels, excluding the results from time channel 25, whose inversion produced an unsatisfactory fit. The results obtained for each time channel when using either the Z-component for the inversion, or using the B-field data were quite similar. Moreover, these calculated parameters for the starting model are quite close to the values of location, depth, and dip reported in the York University Test Site specifications (listed on the bottom line of the table). 5.3.3 Estimation of Time Decay Parameters Two methods that can be used to establish an initial guess of the time parameters were outlined at the end of Section 3.5. The first option is to obtain an inventory of the possible targets in the area, and use the time parameters associated with typical targets in the area. The second method involved fitting the time decay of the Z-component of the measured response directly above the target. This second method is a blind test, and more in the spirit of the thesis, so I will adopt that technique here. Figure 5.14 contains plots of the B-field and dB/dt response directly above the U X O . In 135 Chapter 5. Field Application Depth Time Channel 1 5 9 13 17 21 25 Starting Model Values in Test Site Specifications X 1.80 1.83 1.87 1.88 1.90 1.90 2.08 1.86 Y 2.26 2.12 2.15 2.18 2.18 2.18 2.22 2.18 (Z - loop height) 0.86 0.73 0.65 0.59 0.55 0.54 0.76 0.65 2.00 2.00 0.5 a 0.375 1.4324 0.178 0.0011 1.20 1.62 88.8 0.80 not reported c 37.544 38.723 34.2 30.82 26.0 26.9 91.7 32.5 45 Table 5.1: Results from fitting a single dipole to some of the time channels of the field data set. All lengths are listed in meters, and all angles are in degrees. 5 lines (X = 100,150, 200,250, and 300 cm) with 25 cm station spacing along each line were inverted for each time slice. The above tabulated results are for all 3 components of dB/dt data. Almost identical results were obtained for each time channel when using either the ^-component for the inversion, or using the B-field data. fitting the curve, the data was scaled to ensure that the k parameter will be similar to the k values recovered in the synthetic data tests of chapter 3. The Z-component of the B-field decay curve measured directly above the target is best fit with B^(<)^0.97(f + 0.64)- 00 8 9 e-'/ 512 (5.11) The initial guess for the time decay parameters based on this result is then: a = 0.64, 0 = 0.089, and 7 = 0.75 x 5.12 = 3.84. The dH/dt response is fit with ^ « 71.9(t + 0.43)- - e - ' / 0 29 503 This results in initial guesses of a = 0.43, Q = 0.29, and 7 = 0.75 X 5.03 = 3.75. (5.12) Chapter 5. Field 136 Application Figure 5.14: Fitting the Z-component of the B-field and dB/dt response measured directly above the target to find initial estimates of time decay parameters. Both set of data are plotted with a fitted f(t) - k(t + a)-? exp -t/i. 5.3.4 Summary of Preliminary Data Analysis I choose to use estimates of location and orientation by fitting a single dipole (Section 5.3.2). Time parameter estimates were determined in Section 5.3.3. The initial estimate of the model vector for both B-field and dB/dt data is Usted in Table 5.2. The predicted data are linear in (m ); 0 X Y Z a c fci, k 2 Cti, C t 2 7i. 72 B-field dB/dt 1.86 1.86 2.18 2.18 0.65 0.65 0.80 0.80 32.58 32.5 1 10 0.64 0.43 0.089 0.29 3.84 3.75 Table 5.2: A summary of the estimated decay parameters. parameters fci and k . I choose to scale the field data such that recovered values of ki and k 2 2 are of about the same order of magnitude as those used in the synthetic data sets. I accomplish Chapter 5. Field Application 137 this by making the maximum value of the Z-component of the field data to be about the same order of magnitude as the maximum value of the Z-component in synthetic data sets. 5.4 Inversion of the Field Data Set Up to this point in the chapter, a noise model has been developed which enables us to responsibly assign errors to the data and reject those data that are assumed to be noise. I have also obtained an initial estimate m of the model vector. With these prerequisite procedures completed, I can D now use the least squares algorithm laid out in Chapter 3 to recover the 13 model parameters. 5.4.1 Selection of Data The following inversions will use the 63/dt measurements, rather than the B-field response as data. This decision is, in part, motivated by the noise study presented earlier in this chapter. It was shown that the noise will swamp out the response from a target at early and late times for much of the 63/dt data set. The B-field response is obtained by integrating the measured 63/dt data curve at each station location (Section 4.1.2). If the 63/dt response is not known for late times, then this integration will be inaccurate. This inaccuracy motivates the use of 63/dt. The complete data set consists of 9 lines with 17 stations per line spanning a 4mx4m survey area. I will not use all these stations. For the following inversion I will use all the stations from within a 2mx2m square centered in the 4mx4m survey area. This area is shown in Figure 5.15. In this figure the circles represent locations where measurements were taken, and filled circles represent stations from which data for this inversion is used. Lines X = 100 cm, 150 cm, 200 cm, 250 cm, and 300 cm are used, with stations located at 25 cm intervals, for a total of 45 stations. A l l 30 time channels will be used. Therefore an inversion will have a total of 4050 total data if all 3 components of secondary field are used, and 1350 total data if only the vertical component is used. 138 Chapter 5. Field Application 400 e 350 o 300 e • Measured Decay Utilized 250 a- f o Measured Decay Not Utilized X(cm) 200 o150&100o- 50. 0 i 0 50 100 150 200 250 300 350 400 Y(cm) Figure 5.15: Surficial coverage of data used in the inversions of this chapter. The noise model of method 2 in section 5.2.2 is used. By throwing out data that falls beneath the estimated noise level, the total number of data when using all 3 components will be reduced from 4050 to 2753 data points. When inverting the Z-component of the data only, the total number of data is reduced from 1350 to 1166 data points. 5.4.2 Inversion of Z- component of dB/dt Data Figures 5.16 and 5.17 present the observed and predicted data from the inversion of the vertical component of secondary dB/dt field. The Z-component of observed data is appears oblong in plan view, with greater extent in the X direction, for each of the 30 time channels. As was observed in the noise analysis of this chapter, the early and late time channels start to be dominated by noise. In Figure 5.16, the predicted data reproduce the oblong character of the observed data throughout all the time channels. That is, the predicted data also appears as a single distinct peak elongated in the X direction. The predicted anomaly is slightly wider in the Y direction and narrower in the X direction than the observed anomaly. The difference between d pred and d * for 4 of the time channels are plotted in the third column of Figure 5.16. ob Chapter 5. Field Figure 5.16: Application Observed and predicted data when inverting the Z-component of the dB/dt field. 140 Chapter 5. Field Application Late Times Early Times ,-14 Line 1.0 m •1 Line 2.5 m x 10 x10 1 4 1.5 2 2.5 Station (m) x10 1 Line 1.0 m 1 5 1.5 2 2.5 Station (m) 3 Figure 5.17: Observed and predicted data when inverting the Z-component of the dB/dt field. The left column includes the data from some of the early time channels: t = 0.036 ms (o), 0.092 ms ( A ) , 0.234nw(D), and 0.525 ms(*). The right column of Figure 5.17 include the data from t - 1.003ms(o), 2.53 ms ( A ) , 6.54 ms (O), and 17.2 ms(*). The error bars represent the standard deviation of the error assigned to each data point. The predicted does not fit to within the errors specified by the error bars. Chapter 5. Field 141 Application A n area where the observed data are misfit substantially is on the X = 1 m line, reflecting the fact that the predicted anomaly is narrower in the X direction. Outside of this area, the largest misfits do not appear correlated. Figure 5.17 compares the data misfit hne-by-line. The left column of Figure 5.17 include the data from some of the early time channels: t = 0.036 ms (o), 0.092 ms ( A ) , 0.234 ms (•), and 0.525 ms (*). The right column of Figure 5.17 include the data from t = 1.003 ms (o), 2.53 ms ( A ) , 6.54ms (•), and 17.2 ms (*). The error bars represent the standard deviation of the error assigned to each data point. The predicted does not fit to within the errors specified by the error bars. The observed data on the line (X = 2 m) directly above the apparent location of the U X O are most closely fit. Data lines further away from the target are not being as closely fit. The x 2 measure of the misfit at the solution is x 2 = 4.46 X N. The result that x 2 > N value gives a numerical verification that not all the data were being fit to within the assessed error. This is not an entirely unexpected result, because the forward modelling is approximate, and its error has not been accounted for in the objective function. The recovered parameters and their estimated percent error are Usted in table 5.3. Recall that the U X O was reported to be located, according to our co-ordinate system, at (X, Y) — (200,200) cm and at a depth of 50 cm (Z = 92 cm). The recovered parameters estimate the location of the U X O to be located at (X,Y) = (214.3 ± 0.4, 219.5 ± 0.2) cm. The recovered depth is 127.5 ± 0.6 cm from the surface (164 cm from the transmitter loop). The dip of the U X O is reported to be 45 degrees in the test site specifications, but the recovered dip is 97.8 ± 0.4 °, which is 7 degrees from horizontal and 38 degrees from the specified dip. Again, depending on the precision with which the test site workers were able to bury the U X O , the dip may indeed be an accurate recovery of the true dip. Chapter 5. Field 142 Application Z- Component A l l Components standard standard deviation m c deviation TTli m X 1.84627 0.00723 2.14301 0.00362 Y 2.19203 0.00218 2.19510 0.00227 Z a 1.64292 3.0418 0.00378 0.1994 1.69517 3.2229 0.00556 0.2773 c 63.1 0.941 97.77 0.4263 h 0.00896 0.00008 0.04510 0.03031 0.48116 0.00019 0.07969 0.01412 0.41499 0.00908 4.35356 0.02261 Oil 0i r e c 0.48350 0.34064 r e 7i 1.26933 0.00988 k 0.52055 26.58743 0.44063 OL2 25.46438 0.14252 0.07716 0.17552 0.01625 02 72 0.22130 4.65282 0.01859 0.02247 0.22621 0.00433 0.00974 2 4.65889 Table 5.3: Recovered model parameters when inverting dB/dt field data. 5.4.3 Inversion of A l l Three (X,Y, Z) Components of Data Figures 5.18(a), (b), and (c), and Figure 5.19 present the observed and predicted data from the inversion of all three components of secondary dB/dt field. The data predicted from the recovered model has been able to reproduce the shape of the observed E M anomaly. Chapter 5. Field Application 143 145 Chapter 5. Field Application The recovered parameters and their estimated percent error are Usted in table 5.3. Recall that the U X O was reported to be located, according to our co-ordinate system, at (X, Y) = (200, 200) cm and at a depth of 50 cm (Z = 92 cm). The recovered parameters estimate the location of the U X O to be located at (X,Y) = (184.6 ± 0.7cm, 219.2 ± 0.2cm). The recovered depth is 122 ± 0.4 cm from the surface (164 cm from the transmitter loop). The dip of the U X O is reported to be 45 degrees in the test site specifications, but the recovered dip is 64 ± 4 ° . Again, depending on the precision with which the test site workers were able to bury the U X O , the dip may indeed be an accurate recovery of the true dip. 5.5 Results and Discussion Two different subsets of the York University Test Site field data were inverted in this chapter. In both cases, data acquired at stations within a 2 m x 2 m square centered on the apparent location of the U X O were inverted. In the first case the vertical (Z) component of the data were inverted, and in the second case all three (X, Y, and Z) components of data were inverted. The recovered parameters from the two inversions are Usted alongside the expected values of the parameters in Table 5.4 for comparison. The recovered parameters for both inversions are quite similar, with the exception of the X and c parameters. If we look carefully at the vertical component plot in plan view (see, for example, Figure 5.16), the contours appear to center about a point that is is sUghtly to the right and down from (X, Y) = (2,2)m. Recalling that the X-axis runs vertically along the page and y-axis runs horizontally in order to maintain right handed coordinate system with Z pointing into the page, I would expect that the recovered X is sUghtly less than 2, and Y would be slightly greater than 2. This result is found in the inversion of all the components. The recovered depth of the target is approximately 75 cm deeper than the reported target Chapter 5. Field 146 Application Expected All Components rrii (X,Y,Z) Z-Component Value X Y 1.84627 2.14301 2 2.19203 2.19510 2 1.64292 3.04° 1.69517 3.22° na 63.1° 97.8° 45° *i 0.00896 0.03031 36.1 <*\ 0.48350 0.48116 0.01 01 0.34064 0.41499 0.76 Ti 1.26933 4.35356 28.7 *2 «2 02 72 25.46438 0.14252 26.58743 14.7 0.17552 0.025 0.22621 1.33 50.7 Z (=Depth + Loop Height) a c 0.22130 4.65282 4.65889 0.92 Table 5.4: Comparing expected parameters to recovered model parameters when inverting 63/dt field data. depth of 50 cm. The recovered depth from the surface is 121.3 cm and 126.5 cm for the 3component and Z-component inversions, respectively. These values are actually close to the depth estimated by using the F W H M on the profile parallel to the X-axis (Figure 5.12). The recovered time decay parameters are compared to the time decay parameters obtained by fitting curves in the lab. The two sets of parameters (lab and field) are very different. The recovered value of &i is 3 orders of magnitude less than k from the field data set. The lab data 2 set resulted in k 2 being larger than k\. Since k\ is so much smaller than k , the first dipole 2 has almost no effect on the response. Also, since a.\, 0i, and 71 are all parts of functions that multiply Jbi, they can take on essentially any value without any effect on the predicted data. This fact means that 0i cannot be used as an indicator of whether the object is susceptible, not can it be used in any ratio needed by the discrimination algorithm. Chapter 5. Field All Components Conclusion Result Diagnostic ft ft/ft 147 Application 0.22 ± 0.02 1.54 ± 0.03 (3540 ± 2 ) X 10~ Y ki/k 2 non-permeable n.a. rod-like .Z-Component Result Conclusion 0.226 ± 0.004 1.836 ± 0.002 0.0011 ± 4 X I O - 1 0 non-permeable n.a. rod-like Table 5.5: Results from applying the UXO discrimination diagnostics. The value of ki is very small and therefore the value of 0i is unreliable because it can take on almost any value without affecting the data. Therefore only 0 is used in the diagnostic. 0 indicates that the target is non-permeable. The fc-ratio indicates that the target is rod-like. The diagnostics indicate a non-permeable rod. 2 2 From the inversions requiring all the components and only the Z-component, the ft values are 0.22 ± 0 . 0 2 and 0.226 ± 0.004, respectively. These recovered values are quite small, and not within the range of 0 parameters expected for metallic targets. The minimum value of 8 is that of a non-permeable target which is approximately 0.1 for B-field measurements and approximately 0.5 for dB/dt measurements. The fc-ratio test applied to these recovered parameters indicate that the target is rod-like. The diagnostics as applied to the recovered parameters indicate a non-permeable and rodlike target (Table 5.5). This result is in contradiction to the results established for a 105 m m shell measured in the Geonics lab. The detection diagnostics applied to those measurements indicated a magnetically permeable target. I will consider two possible explanations for this contradiction. The first possibility is that the buried target is indeed non-permeable and the discrimination diagnostics correctly recognized the target's features. This conclusion is reasonable because, although the target is known to be a 105 m m shell, there is no information about the metallic content of the target. In addition, the length of the buried 105 m m shell is 60 cm (Szeto, 1996), while the reported length of the shell measured in the Geonics lab was 40 cm. Clearly, the buried shell is a different type of 105 m m shell than the shell measured at Geonics. It is possible that the particular type of 105 m m shell buried at the York University Test Site is non-permeable. Chapter 5. Field Application 148 If this were the case, the algorithm would be considered successful in establishing the rod-like character of the target. A second possibility is that the buried 105 m m shell is indeed permeable and the recovered diagnostics are incorrect. In this case a possible reason for the failure of the diagnostics could be because the measurement apparatus in the lab is much different than that used in the field. Consequently the data sets and, subsequently, the representative decay parameters are different. Figure 5.1(b) and (c) contained the response of a 105 m m artillery shell measured by Geonics Ltd., using the same lab setup described in Chapter 4. The time decay parameters of the 105 m m shell were extracted from these measurements and listed in Figure 5.1(d). A fundamental assumption made in this thesis is that the parameters recovered in the lab would be the same as the parameters recovered in the field. This assumption allows the U X O identification diagnostics of Chapter 4 that were developed from the laboratory measurements of metallic prisms, to be applied to field data measurements. In the following two sections I assume that, under the same measurement conditions, the buried target that gave rise to the field data set would produce the same response as the 105 m m shell measured in the Geonics laboratory. The validity of this assumption can be examined with the York University data set. If the model parameters from the lab and the field were similar, it would be possible to take the parameters measured at the lab and forward model them to reproduce the measurements made in the field. Therefore, the following short study will examine whether the observed field data can indeed be modelled. If I assume that the location and depth could be adequately estimated by their listed values in the Test Site specifications, then the only real unknown element of the model vector is the dip c. Therefore I will try to reproduce the anomaly pattern of the field data by varying the dip c. This comparison will be carried out on the dB/dt and B-field decay curves measured directly above the U X O ((X, Y) — (200,200) cm), and also on the three component, 30 time channel, full dB/dt data set. Chapter 5. Field 149 Application If a reasonable fit of the field data set could be achieved for some value c, then it is possible that relationships between model parameters and physical parameters of the target from Chapter 4 are applicable to parameters recovered from field data. If it is not possible to attain a reasonable fit of the field data set, then I conclude that significant differences will exist between the response of an object in the lab and the response in the field. Should this be the case, I would surmise that the two data sets (lab and EM61-3D) are inconsistent, and that the interpretation of parameters developed from lab data would be inapplicable to the the field results. C o m p a r i s o n o f L a b a n d F i e l d T i m e D e c a y C u r v e s M e a s u r e d at a S i n g l e S t a t i o n Figure 5.20 contains the observed vertical component of the B-field and dB/dt response meas- ured directly above the U X O . Also included are forward modelled curves using model parameters obtained from fitting laboratory measurements. The model parameters are forward modelled for several different dip angles c. A l l curves are normalized such that the datum of the first time channel for all responses have the same value. There are several differences between the field measurements and the forward modelled curves. Let us first compare the dB/dt measurement to the forward modelled curves (Figure 5.20(b)). There are two major differences between the curves. Firstly, the basic character of the curves is different. The dominant feature of the forward modelled curves is the power law decay t~®, with a flattening of the response at early time, and an exponential decrease at the end. The measured dB/dt response is essentially exponential, and does not have section that appears to obey a power law. Secondly, due to the early onset of exponential decay in the field measured response, the forward modelled responses extend later into time. The differences in the B-field field measured response and forward modelled curves are less dramatic than the differences in the dB/dt response. However, the same differences exist. That is, the forward modelled curves extend later into time and exhibit a power law decay of t~&, and the measured curves look like Chapter 5. Field Application 150 Figure 5.20: Comparison of forward modelling parameters obtained from lab measurements of a 105 mm TJXO, to measurements made at a field site directly over a 105 mm. The model parameters obtained from lab measurements . In both the dB/dt and B-field responses, the forward modelled responses do not look like the EM61-3D field data. a single exponential. In both the dB/dt and B-field responses, the forward modelling was unable to use the lab measured 105 m m shell curves to reproduce the EM61-3D field data at a single point. Comparison of Spatial Behaviour Figures 5.3 and 5.4 plot, in plan view, the measured vertical and horizontal component of the dB/dt response for the odd time channels. The vertical component of the field appears elongated in plan view, and the horizontal component of the field has two lines of symmetry and two maxima in plan view over, essentially, all time channels. It should be possible to reproduce this measured response by forward modelling the decay parameters obtained from the laboratory measurement of a 105 m m U X O , provided that the buried target and the laboratory shell are similar and provided the forward model is accurate. Figure 5.21 contains plots where I try to duplicate the spatial features of the field data set when forward modelling the model parameters from the lab data. In this study I assumed that the U X O was buried 50 cm deep (as was Usted in the York University Field Site target Chapter 5. Field Application 151 specifications). The only model parameter to be varied was the dip c. I was unable to duplicate the features of the full time domain response measured in the field. Indeed, in only one instance was I able to match the features of the field response. When I forward model the decay parameters with the U X O having a dip of c = 90 ° , then the early time measurement of the vertical component of the field is elongated, and the horizontal component has the right symmetry and two peaks. This time slice is indicated by the black arrows in Figure 5.21. This test demonstrates that it is impossible to reproduce the field measurements from the parameters obtained from measurements of a 105 m m shell measured in the Geonics lab. This analysis also reveals the nature of the dipole that reproduces the measured electromagnetic anomaly of figures 5.3 and 5.4. In figure 5.21 the early time behaviour of a U X O with a dip of c = 9 0 ° is shown to have the same anomaly pattern as the field measurements. A t this early time slice the majority of the signal is produced by the dipole perpendicular to the long axis of the target (dipole 2). Therefore, the entire 30 channels of observed data could be reproduced with only dipole 2. The equivalent anomaly can be modelled by the two-dipole model when ki is much smaller than k , which was indeed a result from the inversions. 2 The tests from this section and the previous section were performed assuming that, under the same measurement conditions, the buried target that gave rise to the field data set would produce the same response as the 105 m m shell measured in the Geonics laboratory. It is clear from these tests that if this assumption is accurate, there is some doubt in the validity of transferring parameter interpretation results developed from the lab data analysis to field data. If the decay parameters obtained in the lab are inconsistent with those obtained in the field, then all the diagnostics developed to this point would be inapplicable to field data. Of course, these conclusions were drawn only after assuming that the buried target possessed the same response characteristics as the 105 m m shell in the lab. Further investigation has to be made into determining why the discrepancies exist, such that lab and field measurements can be confidently related. Chapter 5. Field Application 152 A JS 3 o -fl ed -Q -V <+3 ci Xi o 1) 3 .fl H 0< T60 l 1> fl <fl o T3 S V a cu 13 -fl JO o a fl o LO o 3 o CO II o CM II o CO II u II 1) fl s o 5 CO u 60 fl CO a. >1 LO "3 fl B O U o fl c 3 fa .fl 153 Chapter 5. Field Application Time (msec) Figure 5.22: Comparison of dB/dt response of a 105 mm shell. This plot includes the forward modelled response using parameters obtained from the Geonics Lab data ( A ) , the response measured above the UXO using the EM61-3D ( o ) , and the response measured using a modified EM61 with a 5 m square transmitter loop (*) (obtained from Naeva Geophysics). Comparison of the York University EM61-3D Field Data Set to Another Field Data Set These recovered parameters suggest, as was concluded in the short forward modelling study of the previous sections, that the York field data are inconsistent with the U X O discrimination algorithm. The differences in lab and field setup produce data that are too dissimilar to justify applying the interpretation technique, developed from lab data, to the field data. Figure 5.22 plots the Z-component of 83/dt data taken with the target co-axial to the receiver and trans- mitter loops, and again makes the comparison between the Geonics laboratory data and York field data for a 105 m m shell. Comparison of the experimental setup used to acquire the lab data (Figure 4.1) and the EM61-3D arrangement (Figure 1.2) reveal that lab data sets were acquired under different conditions. Two obvious differences between the lab and field setup are that the transmitter loop sizes are very different (a 40 m square loop vs. a i m square loop), and secondly that Chapter 5. Field Application 154 the lab measurements are taken in air while the field measurements obviously has the target buried. It is possible that the difference between the York field test site data and the lab data is, in part, due to the major difference in transmitter loop. A large transmitter loop will produce a more uniform field than a small transmitter loop. Let us consider a dipole source transmitter. The response of a sphere in a dipole field is written as an infinite series of multipoles (Grant and West, 1965; Nabighian, 1970). When the dipole source is far from the sphere, only the term corresponding to a uniform field is significant, and the secondary field of the sphere is then a dipolar field with a strength proportional to the primary field. When the sphere is closer to the dipole transmitter, the higher order multipoles are needed for accurate representation of the response, adding to the response at early times. The l m X l m EM61-3D loop can be approximated by a dipole for cases where U X O are more than 1 to 2 meters from the E M 6 1 transmitter coil (Barrow et al., 1996), however, the range of U X O size for which this rule is applicable was not clearly indicated in the cited paper. For distances less than 1 or 2 meters the response can still, in principle, be represented by a series of multipoles, and a small transmitter loop will induce more higher order multipoles in the target than induced by a large transmitter loop (Nabighian and Macnae, 1991). The effects of source/receiver/target geometries could not adequately be studied in this thesis without an electromagnetic forward modelling code flexible enough to reliably model a range of source/receiver/target geometries. However, there does exist a measurement of a 105 m m shell using a time domain system with a different loop size that provides an interesting comparison. Included in figure 5.22 is a measurement taken by Naeva Geophysics using a modified EM61 of a 105 m m shell. Plotted using '*' markers is the dB/dt decay curve measured as part of the Jefferson Proving Ground Phase IV competition. This is a field measurement directly above a 105mm shell buried at a depth of 0.59m and has a 4 5 ° declination. The shell has a length of 64.8cm, which is approximately 25cm longer than the 105mm shell measured in 155 Chapter 5. Field Application the lab. The measurement uses the same P R O T E M 4 7 receiver that is part of the EM61-3D, but uses a 5 m square transmitter loop. The Naeva Geophysics measurement is plotted alongside the forward modelled curve of a 105mm shell using parameters from the lab data and the York Field Data. The measurement setup for the two field measurements are almost identical. The target is buried at reportedly the same dip, and within 10 cm of the same depth. The major difference is that the transmitter loop size for the Naeva Geophysics configured system is 5 times wider. In comparing the responses of the measured and forward modelled responses plotted in Figure 5.22 it is clear that the response measured by Naeva Geophysics has a different character than the response measured by the EM61-3D. The major difference is in the slope of the dB/dt response. The EM61-3D measurements has an exponential decay, while the Naeva Geophysics measurements approximately follow a r _ 1 decay at earlier times. The Naeva Geophysics measurement is somewhat similar to the forward modelled response using parameters from lab measurements of a 105mm shell. In both these measurements a large transmitter loop, producing a more uniform field was used in the measurement. 5.6 Summary The previous chapters of this thesis outlined the development of a data interpretation strategy for the location and identification of TJXO. In this chapter, the inversion methodology and U X O identification diagnostics were applied to EM61-3D secondary field measurements taken over a buried 105 m m shell. Before applying the algorithm, two methods were suggested for characterizing the noise in the measured data. In the first method, inversions using data from single time channels were performed using a two-dipole model. The misfit for each time channel was used to estimate of the noise in that channel. In the second method, I assumed that field measurements taken Chapter 5. Field 156 Application far enough away from the target did not carry any useful signal from the target, but rather provided a measurement of the noise. These noise measurements could then be fit with a power law function A t . B The two methods provide similar and reasonable estimates of the noise. Inversions using each of the noise models generated very similar results. W i t h an acceptable noise model in place, the inversion procedure was then carried out. Two subsets of the complete York University Test Site field data set were inverted. In both cases only stations within a 2 m square region centered on the target were considered. In the first case only the vertical (Z) component of the 03/dt response was inverted. In the second case all three (X,Y,Z) components of 03/0t, collected within the same square region, were inverted. The two inversions produced similar results. In both cases the \ 2 misfit was approximately AN, where N is the number of data. I consider this misfit to be reasonable, considering the approximate forward modelling used to invert the data. The recovered parameters provide a good estimate of the location on the survey grid of the buried ordnance, but the recovered depth and dip of the U X O do not match the depth and dip Usted in the field site specifications. The algorithm, when appUed to either the 03 jOt or B-field data sets, placed the object at least 70 cm deeper than the reported depth. The recovered depth, however, is consistent with the estimation of depth using the Full Width Half Maximum. The recovered time decay parameters demonstrated that the data could be fit using a single dipole perpendicular to the long axis of the U X O . The diagnostic tests, developed in Chapter 4, indicated that the buried target was non-permeable and rod-Uke. The chapter concludes with a comparison of the field data and the lab measurements of a 105 m m shell. The major assumption in the final sections was that the buried 105 m m sheU and the shell measured in the Geonics laboratory, are similar in the sense that the two sheUs would produce the same response under the same measurement conditions. Given this assumption, the tests showed that the field data measured by the EM61-3D is inconsistent with the lab data used to develop the discrimination diagnostics. Fundamental to the accuracy of Chapter 5. Field Application 157 the forward modelling is the ability to model the response of a buried target as a pair of dipoles at the midpoint of the body. The results of these tests bring into question the validity of the dipole representation. This representation is exact for T D E M response of a sphere in a uniform, step off primary field and also for the magnetostatic solution of a spheroid in a uniform field. Indeed, if you are far enough away from any localized current distribution, the observed field will be that of a dipole (Jackson, 1975). There are a couple of key factors that weaken the assumption that the response can be represented as a single, point dipole for this field data set. Firstly, the primary field is produced by a l m square transmitter loop, which is likely too small to produce a uniform field. Secondly, we are observing a buried target with length equal to one-half the distance from the target center to the receiver. We may be located too close to the target to observe a dipole field. Given the different assumptions that may have been violated, or any other unaccounted errors, tests should be completed to determine under which source/receiver/target geometries the technique will be ineffective or successful. The possibility remains that the object buried at the York test site is not strongly magnetic and, in fact, is properly represented as a conducting, non-permeable rod. I suggest that the buried object be dug up and subjected to in-house analysis, or that the experiment of burying a U X O with measured laboratory responses be repeated. Chapter 6 Summary, Conclusions, and Future Work The goal of this thesis was to develop a method of discriminating between U X O and nonU X O items from time-domain electromagnetic data collected with the Geonics L t d . prototype instrument EM61-3D. In Chapter 2 an approximate forward model was developed for describing the time domain electromagnetic response of a buried metallic object. The model was based on the analytic solution of a sphere and the magnetostatic solution of a spheroid. It has a strong physical basis and is simple to compute. The time domain response of a metallic object is represented as a pair of perpendicular, decaying dipoles located at the center of the object. In Chapter 3 a non-linear parameter estimation procedure was developed to recover the 13 parameters of the model. Rough estimates of the model parameters can be obtained by applying some simple data preprocessing strategies. These initial estimates are then iteratively refined by solving a nonlinear least squares problems. These techniques were shown to work very well on synthetically generating data sets containing random, Gaussian noise. Chapter 4 outlined a proposed U X O discrimination method. In the absence of analytic solutions of Maxwell's equations and the lack of suitable numerical electromagnetic forward modelling, I use lab measurements of the time decay curves of difFerent metallic targets. The algorithm is based on two facts: (1) U X O are typically made of steel and (2) U X O are typically rod-like rather than plate-like. The discrimination method first requires that conductive targets , be distinguished from steel targets that are conductive and permeable. The 0 parameter was found to be an indicator as to whether the target is indeed permeable. 158 The second stage 159 Chapter 6. Summary, Conclusions, and Future Work then focuses upon determining if the steel target is rod-like or plate-like. The ratios of ki/k 2 and 0i/f$2 are diagnostic indicators as to whether the geometry is plate-like (uninteresting) or rod-like (a high candidate for being a TJXO). In Chapter 5 the parameter estimation procedure and discrimination algorithm is tested on a EM61-3D field data set. The data set was acquired over a 105-mm shell buried at the York University Geophysical Test Site. Techniques were developed in this chapter to estimate the noise level in the instrument during the survey. The standard deviation of the noise was found to follow a t - 1 decay. The parameter estimation algorithm was applied to a data set including all three (X, Y, and Z) components of dH/dt data, and to a data set containing only the vertical (Z) component of the data set. Both inversions returned approximately the same decay parameters, and the discrimination algorithm identified the buried target as a non-permeable rod-like target when applied to either data set. The chapter concludes with comparisons between a T D E M measurement of a 105 m m U X O in a laboratory setting and a field setting. The two electromagnetic anomalies are quite different, and therefore we must question how similar the 105 m m U X O buried at the York University Test Site is to the 105 m m U X O measured in the lab, and how does the differences the source/receiver/target geometries of the lab and field measurements affect the time domain response. Future Work Much work needs to be done before the discrimination algorithm can be used as a practical tool. Most importantly, the accuracy of the approximate forward modelling needs to be established. As was mentioned earlier, tests need to be completed to determine the conditions under which the approximate forward modelling is inadequate. Access to an accurate numerical solver that is flexible to source/receiver/target geometries will aid tremendously in this regard. Non-UXO buried metal objects will often be buried in close proximity to each other. A study of the T D E M response to multiple targets should be performed: It has been suggested Chapter 6. Summary, Conclusions, and Future Work 160 that the use of a magnetometer in concert with a T D E M sensor would provide the necessary information to identify 'clutter' (Barrow et al., 1996). Magnetometer surveys are performed on U X O remediation projects, and therefore a 'joint' inversion of magnetics and electromagnetic data sets would be useful. Appendix A Time Decay Parameters of Various Unexploded Ordnance In this appendix the time decay parameters k, a, 3, and 7 are Usted for the measurements of a 24 UXO provided by Geonics. Details of the measurement procedure are provided in Chapter 4. The 63/dt measurements were suppUed, and the B-field response was obtained by integrating dB/dt (see Chapter 4). The decay parameters are obtained by fitting the curve with least squares. A weighted least squares objective function: <p=\\W (f(k, ,3,r,t)-d° °)\\ b d (A.l) 2 a is minimized where (Wd)u = l/d° . bs This weighting is chosen to put equal importance on each data point in the fitting. The minimization is performed using Newton's method. 161 Appendix A. Time Decay Parameters of Various Unexploded Ordnance UXO 1 Length Diameter Weight 18 cm 6.05 cm 1 0 kg dB/dt k a 7 Orient. 1 8.4264 2.1278e-13 0.76425 4.4249 Orient. 2 1.034 0.0046863 1.3482 7.7877 B-field k a 0 7 Orient. 1 0.015929 0.14318 0.38228 4.0073 Orient. 2 0.0014048 0.020048 0.57791 5.0536 Orient. 1 0.03313 0.099158 0.28132 4.5931 Orient. 2 0.0048181 0.025894 0.53593 3.9331 Orient. 1 0.035966 0.13089 0.345 14.6366 Orient. 2 0.005295 0.037148 0.52799 8.2811 UXO 2 Length Diameter Weight 22 cm 7.64 2.9 kg dB/dt k a 0 7 Orient. 1 15.3405 0.00145 0.71778 5.34 Orient. 2 3.7056 0.00045603 1.1869 4.9646 B-field fe a 0 7 UXO 3 10 -i 10 Length Diameter Weight 23 cm 7.64 6.0 kg dB/dt fe a 0 7 to" 10° Time (msec) Orient. 1 12.018 0.0076046 0.89309 17.6534 10' 10 Orient. 2 3.2996 0.0076551 1.2523 11.5385 B-field fe a 0 7 162 Appendix A. Time Decay Parameters of Various Unexploded Ordnance dB/dt Length Diameter Weight 28 cm 8.28 cm 2.4 cm k a 0 1 Orient. 1 16.6154 0.0011347 0.72135 9.1392 Orient 2 3.1425 0.013708 1.3689 14.0355 Orient. 1 17.283 3.9309e-13 0.66067 7.0206 Orient. 2 4.7022 0.011705 1.307 9.0696 B-field k a & 7 163 Orient. 1 0.047564 0.11429 0.25381 7.4403 Orient. 2 0.0047484 0.039837 0.59659 9.2161 Orient. 1 0.050894 0.19161 0.24897 5.7971 Orient. 2 0.0069736 0.05303 0.60847 6.8361 UXO 5 dB/dt Length Diameter Weight 27 cm 7.96 cm 3.5 kg k a 0 7 B-field k a 0 7 UXO 6 10 Time (msec) Time (msec) dB/dt Length Diameter Weight 23 cm 7.64 3.0 kg k a 0 7 Orient. 1 13.4403 0.0015016 0.69969 9.2511 Orient. 2 6.8116 0.00044706 1.024 2.7785 B-field k a 0 7 Orient. 1 0.049723 0.5342 0.44479 8.9712 10 Orient. 2 0.0083152 0.017431 0.40575 2.2748 Appendix A. Time Decay Parameters of Various Unexploded Ordnance 164 UXO 8 10 Time (msec) Time {msec) Length Diameter Weight 31 cm 7.64 cm 5.5 kg dB/dt k a 7 Orient. 1 18.7578 0.021688 0.79351 22.2871 Orient. 2 5.4074 0.013519 1.2478 8.175 B-field k a 0 7 Orient. 1 0.10558 0.80482 0.48076 21.5926 10 Orient. 2 0.0088183 0.023736 0.44738 5.5922 UXO 9 10 Length Diameter Weight 32 cm 9.23 cm 10 kg dB/dt k a a 7 10 Time (msec) Orient. 1 18.5381 2.8288e-09 0.65171 27.1429 Time (msec) Orient. 2 8.3385 0.015513 1.2549 67.7316 B-field k a 7 Orient. 1 0.11009 0.29625 0.18691 0.18691 Orient. 2 0.019688 0.045709 0.42852 31.7268 Appendix A. Time Decay Parameters of Various Unexploded Ordnance 165 81 mm Mortar U X O 10 iu 10 Length Diameter Weight 43 cm 8.28 cm 3.7 kg dB/dt k a 0 7 ' 10"' 10° Time (msec) Orient. 1 19.6953 0.0038303 0.80704 3.5975 10' 10 Orient. 2 5.2157 0.0041505 1.1684 2.6491 B-field k a 0 7 Orient. 1 0.030839 0.1016 0.39821 3.224 Orient. 2 0.0054725 0.032929 0.56827 2.1543 81 mm Mortar U X O 11 10 Time (msec) Length Diameter Weight 43 cm 8.28 cm 3.5 kg dB/dt k a 0 7 Orient. 1 19.7178 2.8314e-05 0.67757 4.5845 Orient. 2 4.9078 0.015947 1.3145 4.5656 B-field k a 0 7 Orient. 1 0.04106 0.14366 0.30994 4.1266 10 Orient. 2 0.0060188 0.081855 0.73637 3.8707 Appendix A. Time Decay Parameters of Various Unexploded Ordnance 166 U X O 13 'V Length Diameter Weight 57 cm 7.96 cm 3.2 kg dB/at k a 0 7 10"' 10° Time (msec) Orient. 1 23.2215 5.5645e-10 0.5927 2.4582 io' 10 z Orient. 2 10.0824 0.006744 1.1617 2.0031 B-field k a 0 7 Orient. 1 0.037261 0.10843 0.24658 2.1857 Orient. 2 0.0092283 0.036822 0.59621 1.7036 Orient. 1 0.17694 0.52649 0.20864 20.6112 Orient. 2 0.034207 0.09479 0.56972 14.9911 Orient. 1 0.16956 0.21329 0.24762 23.9861 Orient. 2 0.031181 0.070159 0.51925 30.7427 81 mm Mortar U X O 14 10 ^ 10 Length Diameter Weight 42 13.37 12.5 dB/dt k a 0 7 10"' 10° Time (msec) Orient. 1 26.7977 1.3556e-08 0.65935 27.1799 10' 10 Orient. 2 17.6867 0.020031 1.2578 19.9042 s B-field k a 0 7 U X O 15 dB/fo' 10 Time (msec) Length Diameter Weight 40 cm 10.5 13.5 kg dB/at k a 0 7 Orient. 1 36.147 0.010285 0.75595 28.7217 Orient. 2 14.7409 0.025026 1.339 50.7483 B-field k ot 0 7 Appendix A. Time Decay Parameters of Various Unexploded Ordnance 167 U X O 16 10 dB/dt Length Diameter Weight 48 10.19 7.0 k dc 7 Orient. 1 36.7492 0.00061726 0.61302 6.1913 Orient. 2 12.1948 0.0080828 1.2124 3.7363 B-field k ot P 7 10 Time (msec) 10 Orient. 1 0.10094 0.22683 0.27539 5.648 Orient. 2 0.014044 0.036692 0.57446 3.0353 Orient. 1 0.273525 0.47898 0.33424 11.6425 Orient. 2 0.034819 0.037236 0.46327 7.8332 Orient. 1 0.19212 0.3148 0.1801 14.0868 Orient. 2 0.031204 0.039386 0.40618 6.9324 U X O 17 10"' 10° ,„ Time (msec) 10 dB/dt Orient. 1 Length Diameter Weight 40 10.82 12.0 k 48.5862 0.0011019 0.65187 12.5205 a 0 7 to' 10 B-field Orient. 2 20.5593 0.0068033 1.1529 9.9621 k Ot 0 7 U X O 18 10 10-! 10' - 10° 10' 10 2 Time (msec) dB/dt Length Diameter Weight 48 10.82 10.5 k a fi 7 Orient. 1 36.2589 5.7907e-ll 0.603 17.1082 Orient. 2 17.6064 0.010139 1.0819 8.7223 B-field k a & 7 Appendix A. Time Decay Parameters of Various Unexploded Ordnance 168 U X O 20 10 Time (msec) Time (msec) dB/dt Length Diameter Weight 65 cm 10.82 cm 10.0 kg k ot 0 7 Orient. 1 7.5588 1.3271e-09 0.71667 18.2876 Orient. 2 16.9015 0.025302 1.1765 12.5304 B-field k a 0 7 Orient. 1 0.030871 0.12357 0.20726 13.9887 10 Orient. 2 0.034962 0.19526 0.65286 0.11.1432 U X O 21 dB/ft< 10 Time (msec) dB/dt Length Diameter Weight 85 12.1 13.5 k a 0 7 Orient. 1 5.0162 0.00042785 0.74852 33.8514 Orient. 2 7.4084 0.023385 1.2394 2.3ell B-field k a 0 7 Orient. 1 0.02588 0.20311 0.22922 25.0263 10 Orient. 2 0.026974 0.036862 0.28924 388.9814 Appendix A. Time Decay Parameters of Various Unexploded Ordnance 169 UXO 23 10 Time (msec) Time {msec) Length Diameter Weight 58 cm 10.19 cm 10.0 kg dB/dt k a 0 7 Orient 1 19.4658 0.032818 1.0805 11.4199 Orient. 2 10.2121 0.013102 1.4264 6.7625 B-field k a 0 7 Orient. 1 0.052432 0.37856 0.73813 11.577 10 Orient. 2 0.012149 0.042451 0.71011 5.135 UXO 24 10 Time (msec) Length Diameter Weight 77 11.46 16 dB/dt k a 0 7 Orient. 1 19.4658 0.032818 1.0805 11.4199 Orient. 2 10.2121 0.013102 1.4264 6.7625 B-field k a 0 7 Orient. 1 0.052432 0.37856 0.73813 11.577 10 Orient. 2 0.012149 0.042451 0.71011 5.135 Appendix B Time Decay Parameters of Metallic Prisms In this appendix the time decay parameters k, a, fl, and 7 are listed for the measurements of steel and aluminum axi-symmetric targets provided by Geonics. These targets range from an 8 X 8 X 0.25 inch plate to a 8 X 0.25 X 0.25 inch rod. Details of the measurement procedure are provided in Chapter 4. The dB/dt measurements were supplied, and the B-field response was obtained by integrating dB/dt (see Chapter 4). The decay parameters are obtained by fitting the curve with least squares. A weighted least squares objective function: <j>=\\W (f(k,a,(3, ,t)-d° ')\\ b d (B.l) 2 r is minimized where (Wd)u = l/d° . bs This weighting is chosen to put equal importance on each data point in the fitting. The minimization is performed using Newton's method. 170 -»-> « a CO rH 00 00 X 00 X | 0.0049997 || | 0.087004 || |[ 0.38399 || | 14.9615 || | .95053 | 2.0745 | | 0.0017295 | 0.020377 | 0.035618 | | 0.048165 0.50746 | 1.3023 | 1.1128 | | 12.8818 | 22.5796 | 22.2348 | | CQ. | | CO cq 00 X X 00 CO co CO X X oo cq CJ <- -a a ca t- a ca c- CO 00 cq X 00 X Tf CO 8x CO rH X a ca 00 r H X X 00 rcq CO CO CO 8x0.5 xO.5 LO 00 X 00 coo CO rH CO 8x0.25 xO.25 Jc a ca 1.6748 |l| 0.00049227 | 0.0025385 || 0.30464 | || 0.0086219 | 0.0089813 ||| 0.018756 | 0.069077 || | 0.82183 || 0.49158 | 0.34433 || || 1.3212 3.3141 || 5.6646 | 2.7478 || || 11.6526 | •* O CQ. f- || 0.00069186 0.0030892 || 0.018414 | 0.09949 || 0.45907 | 0.38939 || || || 22.528 | 8.7827 || X e 1.3705 0.020107 0.9948 11.9242 X c- -K | | | | OO 0.34032 0.012776 1.3853 90.1963 t~ a ca | | ;| j| QO || 0.20299 | 0.0091345 || 0.72689 || 3.6951 | | | 0.014699 0.20098 0.45132 34.848 || || || || | 1.8608e-07 | 0.049106 || | 1.8329 || | 0.26928 || | 0.0013924 ||| 0.00035026 | 1.4054e-07 | | 0.013902 || | 8.9551e-05 ||| 0.10199 | 1.0787 || | 1.2619 ||| 0.31779 | 0.18193 ||| 3.2135 | 0.35199 || 0.001028 0.16903 0.3054 15.5135 | 0.03609 || 0.0014052 | 2.7009e-05 || | 0.011574 || | 0.0034494 || 0.14574 | 0.61933 || | 1.3037 || 0.35792 | 1.0041 || 13.3861 | 0.74721 || | | 0.0023002 | 0.0001559 || 0.10281 | 0.017102 || || || 0.36207 | 0.56658 || 1.8703 || || 15.5103 || | | 0.0049244 | 0.0018511 || 0.16051 | 0.060403 || 0.42977 | 0.5095 || || || 17.9691 | 11.0643 || | 0.30101 | 0.015712 || | 0.021459 | 2.1399e-05 || | 0.85908 1 1.0295 || | 20.6334 | 0.11215 || | 0.48849 | 0.019931 | 0.92726 | 17.4725 | | 0.84071 | 0.14034 0.019131 | 0.0070164 1.0078 1 1.3203 22.4194 | 2.9159 | | | 1.039 | 0.028859 | 1.2895 | 18.662 •4C r- 1.565 | | 0.00088423 | 0.0040879 || || 0.024692 || 0.029362 | 0.071046 || 0.51093 | 0.32517 || | 1.0288 || | 13.3727 || 5.9466 | 19.5259 || X -« a ca 0.56951 0.01287 1.3008 10.2103 QO 1.8320 0.03055 1.0482 23.5033 r- | | | | d | | | | rH X CO CO 0.0084561 | 0.18567 || | 0.43491 || | 18.6487 || -« 8 | | 0.0042864 | | 0.075231 0.45556 || 14.3573 || o 2.0124 | 2.9722 0.020758 | 0.028329 1.1724 1 1.0217 21.6217 | 23.4459 "Cf rH 0.20098 0.45132 34.848 to | | 4.8623 | 4.8623 || 0.040902 [ 0.040902 |i| 1.1097 | 1.1097 || || 46.867 | 46.867 Oi 1 | | | B-fi Orient. 1 | Orient. 2 || o> to 00 X 00 xO.25 CO | 0.014699 || | 0.20098 || | 0.45132 || | 34.848 || M 0.014699 0.20098 0.45132 34.848 o 4.8623 1 4.8623 || 0.040902 | 0.040902 || | 1.1097 | 1.1097 ||| | 46.867 | 46.867 ||| 1 <L> Orient. 1 | Orient. 2 || 3 x0.5 PH B-field || Orient. 1 | Orient. 2 || Orient. 1 | Orient. 2 || Appendix B. Time Decay Parameters of Metallic Prisms 171 a ca -tf a ca c- -s a ca c- HO CO 3 CO rS£ < rH || || || || || || || | || || || | 0.0063618 0.0038149 | 0.49342 | 3.9367e-13 | | 0.29248 | 0.40247 0.5537 | | 0.11297 | 0.15264 | 31.2522 | | 31.5899 | 24.1568 0.3568 | | 0.0050647 | 0.0016344 | 0.60208 | | 1.1304e-10 | 2.9783e-09 || 0.48497 | 0.22879 0.15008 | 0.52224 || 0.16367 | | 0.56652 9.0213 | 37.9015 | 11.277 || 30.4074 | 0.5726 | 0.3316 || 0.0043697 0.00082797 0.10958 | 5.3588e-09 | 2.2939e-ll || 0.53861 | 0.14249 | 0.53546 | 0.492 || 0.16692 | 3.1996 | 26.4055 | 4.0158 || 21.0603 | X oo OO X a oo X X oo <! J i < TC X X a <Q. CS oo 00 X Tf X CM rH Tf < 8x oo Tf a CN CN X -a a CO CO CS X 00 X rH < C^ < Tf oo a rS a 00 to CO 00 X <- X X rH ~K Cft rH < 8x0.5 a 0.0033937 0.00010717 | 0.52239 || 0.93802 0.87303 ]) || 4.0679e-10 | 1.7237e-ll |:| 0.52586 | 5.4038 || | 0.31702 ||| 0.20237 | || 0.37911 || 5.3196 | 1.3617 || | 0.21948 | j| 5.6039 OO 0.69953 6.118e-05 0.53195 T"H OO | | | || || || || < a 00 X xO.25 M 0.012039 1.0317 0.14904 74.4408 Je xl fl—' 1 | | | | || 0.012039 1.0317 0.14904 74.4408 0.11299 3.514e-10 0.56293 8.0072 0.19145 | | 0.00038641 | 0.00080864 || | 3.9402e-10 || 0.11184 j 0.17332 || | 0.15715 || | 0.55379 || 0.15123 5.8912 | 8.6427 || | 11.2565 || || || || | 0.0028411 | 0.40512 | 0.16213 | 21.4616 | | | | 0.040422 0.0038343 0.61374 2.25 | 0.094128 | 0.0032119 | 0.6008 | 2.1491 l 7S1 || 5.9427e-05 | 0.0001375 || | 0.04371 || || 0.05068 | 0.19534 || || 0.21149 || 1.8563 l ll I 0.028395 | 0.051813 || 1.2985e-05 | 3.015e-05 || l| 3.0485e-09 | 2.5515e-09 || 0.027442 | 0.053631 || | 0.25051 || || 0.45112 | 0.46347 ||| 0.21556 | 0.62544 || || 0.45513 | 0.64058 ||| 0.43618 | | | | | 0.28784 | 0.39413 1 | 0.0020248 0.34392 | 2.3614e-09 | 2.9448e-10 | | 0.15097 | 0.5406 | 0.537 I| 20.5874 | 25.81 | 26.2191 || || || || || || 0.012039 1.0317 0.14904 74.4408 | | | | | Orient. 2 B-fi | 0.8079 | 0.8079 || | 0.00043294 | 0.00043294 | | | 0.54751 | 0.54751 || || | 94.3956 | 94.3956 | Orient. 1 | - Orient. 1 | Orient. 2 |!| CO 0.8079 || 0.012039 | 0.8079 | 1.0317 | 0.00043294 | 0.00043294 || ||| 0.14904 | 0.54751 | 0.54751 | 94.3956 | 94.3956 |l| 74.4408 Orient. 2 ffl || Orient. 1 | «0 Orient. 2 o | .3 Orient. 1 Appendix B. Time Decay Parameters of Metallic Prisms a <a LO d X 172 Appendix C Calculation of the primary field B P from a square loop The B-field of square loop carrying a steady state current I can be determined by integrating the Biot-Savart law. The Biot-Savart law is simply fi I f dl X r' 0 ^ 47T To find the field of a loop I will determine an expression for the field from a finite segment of wire. The field from a square loop follows easily from this result. Let us consider a wire whose endpoints are located at r a and rj, (Figure C . l ) . Y w tfi ,. di = di i b u a b Figure C . l : A current carrying wire of finite length. Let I b be the unit vector pointing from r a a to rj,. We can first solve for the numerator of equa- tion C . l , and secondly the denominator, and then we can substitute the appropriate expressions in equation C . l . From figure C . l , r' = r 0 - / I b and dl = dllab. By Sine law, we can solve for a / and its derivative dl I = r (sin a cot (j> — cos a ) , dl — —r sin a esc <pd(f> 2 a a (C2) The numerator is then dl X r' = —r sin a esc <p (l b X r ^ d<f> 2 a a a 173 (C3) Appendix C. Calculation of the primary field B from a square loop p 174 The denominator is then * |r'| = r - / I a a 0 (CA) = r — r (sin a cot <fr — cos a) I j a a a = A — r sin a cot ^>Ib a a where A = r + r cos a l a a a 0 Substitution into Biot-Savart to gives us the field from a straight line pI [ ' 4TT J^ <t b=l3 0 >ab ~ ^ (* ) ^ | A - r sin a cot <f> I \ Ta s i n a c s c 2 afe x Ta d< (C.5) 3 A=A a ab To solve this integral make a change of variable by letting x = r sin a cot </>, and therefore a dx = — r sin a csc <f>d<f>: 2 a B A B ~ I F / '!Bo a (C.6) 3!l J,P — . r ) da: X (lab a This integral can be rewritten as da; B «» = ^ r ( ^ x r - ) / . [c + 26a; + aa; ] / 2 3 2 (C.l) aa; + b 47T ( V A C _62)( c + 26a; + A A ;2) / 3 2 where a = 1, 6 = - A • I , and c = A • A . Substitution to return to our original variables give a 0 <t>=/3, r'=r b {lab X r ) pI 0 / r q s i n a c o t ^ - (r a 47T r 2 (r -l y - a a •I ab - r cos a) a (C.8) r'=r ab a W i t h some manipulation we get /X J P (Iqb X r ) a •. /£a_£b\ ^ " * [r!-(r..U) ] 'U J r > U _ (C.9) where I used the fact that cos a lab •r a lab ' fb cos (3 = — r-b (CIO) Appendix C. Calculation of the primary held B p from a square loop 175 Equation C.9 therefore expresses the field of a straight wire of finite length in terms of the current in the wire J , and the endpoints of the line segment, given by r lab = ( r & - r ) / | r - r | ) . 0 6 0 v. 0 and r b (note that References Bard, Y . , 1974, Nonlinear Parameter Estimation: Academic Press. Barrow, B . J . , Khadr, N . , and Nelson, H . H . , 1996, Performance of Electromagnetic Incution Sensors for Detecting and Characterizing U X O : Proceedings of U X O Forum '96. Bertotti, G . , 1998, Hysteresis in Magnetism: Academic Press. Brown, W . F . , 1962, Magnetostatic Principles in Ferromagnetism: North-Holland Publishing Co. Butler, D . K . , 1997, Landmines and U X O : The Leading Edge, 16, no.10, 1460. Butler, D . , Cespedes, E . , O'Neill, K . , Arcone, S., Llopis, J . , Curtis, J . , Cullinane, J . , and Meyer, C , 1998, Overview of Science and Technology Program For J P G Phase IV: Proceedings of U X O Forum '98. Chikazumi, S., 1997, Physics of Ferromagnetism, 2nd ed.: Oxford University Press. Das, Y . , McFee, J . E . , Toews, J.and Stuart, G . C . , 1990, Analysis of an Electromagnetic Induction Detector For Real-Time Location of Buried Objects: I E E E Transactions on Geoscience and Remote Sensing, 28, no.3, 278-287. Dennis, J . E . , and Schnabel, R . B . , 1983, Numerical Methods for Unconstrained Optimization and Nonlinear Equations: Prentice-Hall. Effers0, F . , Auken, E . and S0rensen, K.I., 1999, Inversion of Band-Limited T E M Responses: Geophysical Prospecting, 47, 551-564. 176 177 References Federal Advisory Committee for the Development of Innovative Technologies, 1996, "Unex- ploded Ordnance (UXO): A n Overview." Feynman, R., 1964, The Feynman Lectures On Physics, volume 2: Addison-Wesley Publishing Company. Gill, P . E . , Murray, W . , and Wright, M . H . , 1981, Practical Optimization: Academic Press. Grant, F.S., and West, G . F . , 1965, Interpretation Theory in Applied Geophysics: McGraw-Hill, Inc. Grimm, R . E . , Blohm, M . W . , Lavely, E . M . , 1997, U X O Characterization Using Multicomponent, Multichannel Time-Domain Electromagnetic Induction: Proceedings of U X O Forum '97. Jackson, J . D . , 1975, Classical Electrodynamics: John Wiley and Sons. Kaufman, A . A . , 1994, Geophysical Field Theory and Method: Academic Press, Inc. Kaufman, A . A . , and Keller, G . V . , 1985, Inductive Mining Prospecting: Elsevier. Khadr, N . , Barrow, B . J . , and Bell, T . H . , 1998, Target Shape Classification Using Electromagnetic Induction Sensor Data: Proceedings of U X O Forum '98. Lorrain, P., and Corson, D . , 1970, Electromagnetic Fields and Waves, 2nd edition: W . H . Freeman and Company. McNeill, J . D . , and Bosnar, M . „ 1996, Application of Time Domain Electromagnetic Techniques to U X O Detection: Proceedings of U X O Forum '96. More, J . J . , and Sorensen, D . C . , 1984, Newton's Method, in Golub, G . H . , E d . , Studies in Numerical Analysis: Mathematical Association of America, 29-82. References 178 Munkholm, M . S . , and Auken, E . , 1996, Electromagnetic Noise Contamination on Transient Electromagnetic Soundings in Culturally Disturbed Environments: Journal of Environmental and Engineering Geophysics, 1, no.2, 119-127. Nabighian, M . N . , 1970, Quasistatic Transient Response of a Conducting Sphere in a Dipolar Field: Geophysics, 35, no.2, 303-309. Nabighian, M . N . , and Macnae, J . C , 1991, Time Domain Electromagnetic Prospecting Methods, in Nabighian, M . N . , E d . , Electromagnetic Methods In Applied Geophysics: Society of Exploration Geophysicists, Volume 2, Application, Part A , 427-520. Pavlov, D . , and Zhadnov, M . , 1998, Anomalous Conductivity and Magnetic Permeability Effects in Time Domain Electromagnetic Data: Proceedings from the 68th Annual Internat. M t g . , Society of Exploration Geophysicists, Volume 1, 440-443. Office of the Undersecretary of Defense (Acquisition & Technology), 1998, "Report of the Defense Science Board Task Force on Unexploded Ordnance ( U X O ) Clearance, Active Range U X O Clearance, and Explosive Ordnance Disposal ( E O D ) Programs." Pasion, L . , and Oldenburg, D . W . , 1997, The Inversion of Time Domain Electromagnetic Data for the Detection of Unexploded Ordnances: JACI1997 Ann.Rep., Dept.of Geophysics and Astronomy, Univ.of British Columbia. Szeto, A . M . K . , 1996, York University Environmental Test Site, Department of Earth and Atmospheric Sciences, York University. Varah, J . M . , 1990, Relative sizes of the Hessian terms in nonlinear parameter estimation: S I A M Journal on Scientific and Statistical Computing, 11, n o . l , 174-179. \ References 179 Ward, S.H.., and Hohmann, G . W . , 1991, Electromagnetic Theory for Geophysical Applications, in Nabighian, M . N . , E d . , Electromagnetic Methods In Applied Geophysics: Society of Exploration Geophysicists, Volume 1, Theory, 131-311. Weaver, J . T . , 1994, Mathematical Methods for Geo-electromagnetic Induction: John Wiley and Sons, Inc. \
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Detecting unexploded ordnance with time domain electromagnetic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Detecting unexploded ordnance with time domain electromagnetic induction Pasion, Leonard Rodriguez 1999
pdf
Page Metadata
Item Metadata
Title | Detecting unexploded ordnance with time domain electromagnetic induction |
Creator |
Pasion, Leonard Rodriguez |
Date Issued | 1999 |
Description | In this thesis I assume that the Time Domain Electromagnetic (TEM) response of a buried axisymmetric metallic object can be modelled as the sum of two dipoles centered at the midpoint of the body. The strength of the dipoles depends upon the relative orientation between the object and the source field, and also upon the shape and physical properties of the body. Upon termination of the source field, each dipole is assumed to decay as k (t + α)-β e-t/γ. The parameters k, α, β and γ depend upon the conductivity, permeability, size and shape of the object, and these can be extracted from field or laboratory measurements by using a nonlinear parametric inversion algorithm. An investigation carried out using an analytic solution for a sphere and laboratory measurements of steel and aluminum rectangular prisms, suggest the following methodology. The value of β might be used as a diagnostic to assess whether the metallic object is non-magnetic or magnetic. If the object is thought to be magnetic, then the ratios of k1/k2 and β1/β2 are diagnostic indicators as to whether the geometry is plate-like (uninteresting) or rod-like (a high candidate for being a UXO). [Scientific formulae used in this abstract could not be reproduced.] |
Extent | 9250282 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-07-07 |
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.0052881 |
URI | http://hdl.handle.net/2429/10329 |
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 | 2000-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2000-0135.pdf [ 8.82MB ]
- Metadata
- JSON: 831-1.0052881.json
- JSON-LD: 831-1.0052881-ld.json
- RDF/XML (Pretty): 831-1.0052881-rdf.xml
- RDF/JSON: 831-1.0052881-rdf.json
- Turtle: 831-1.0052881-turtle.txt
- N-Triples: 831-1.0052881-rdf-ntriples.txt
- Original Record: 831-1.0052881-source.json
- Full Text
- 831-1.0052881-fulltext.txt
- Citation
- 831-1.0052881.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-0052881/manifest