Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Detecting unexploded ordnance with time domain electromagnetic induction Pasion, Leonard Rodriguez 1999

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

Item Metadata

Download

Media
831-ubc_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

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.  \  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items