ESTIMATION A N D CORRECTION OF WAVELET DISPERSION IN GROUND P E N E T R A T I N G R A D A R DATA by JAMES D. IRVING B.Sc, University of Waterloo, 1997 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department of Earth and Ocean Sciences) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April, 2000 © James D. Irving, 2000 In presenting this thesis in partial fulfillment 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 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 The University of British Columbia 2219 Main Mall Vancouver, BC, Canada V6T 1Z4 Date: 11 Abstrac t The attenuation of electromagnetic (EM) waves in many geological materials is strongly dependent upon frequency in the ground penetrating radar (GPR) range; high frequencies are attenuated much more quickly than lower ones during propagation. For this reason, the GPR wavelet often undergoes a significant change in shape as it travels through the subsurface, and reflections received at later times contain less high frequency information than those received at earlier times. This phenomenon is known as wavelet dispersion. In the GPR image, it is displayed as a characteristic "blurriness" that increases with depth. Correcting for wavelet dispersion in GPR data is an important signal processing step that should be performed before either qualitative interpretation or quantitative determination of subsurface electrical properties are attempted. Previous work by other researchers has shown that the EM wave attenuation parameter for many geological materials is approximately linear with frequency over the bandwidth of a GPR wavelet. Thus, the change in shape of a GPR pulse as it propagates can often be well described using one parameter, Q*, which is related to the slope of the linear region. In this thesis, we confirm and build on these results. Assuming that all subsurface materials can be characterized by some Q* value, the problem of estimating and correcting for wavelet dispersion in GPR data becomes one of determining Q* in the subsurface and deconvolving its effects through the use of an inverse Q filter. A method for the estimation of subsurface Q* from GPR data based on a technique developed for seismic attenuation tomography is presented. Essentially, Q* is determined from the downshift in the dominant frequency of the GPR wavelet with time down a trace. Once Q* has been obtained, an inverse Q filtering technique based on a causal, linear, model for constant Q wave propagation is proposed as a means of removing wavelet dispersion. Tests on field data collected near Langley, British Columbia indicate that these methods are very effective at enhancing the resolution of the GPR image. iii Table of Contents Abstract ii List of Tables v List of Figures vi Acknowledgements viii 1 Introduction 1 1.1 The GPR Method 1 1.2 The Motivating Problem 2 1.3 Thesis Overview 4 2 Electromagnetic Theory 5 2.1 Maxwell's Equations 5 2.2 Constitutive Relationships 7 2.3 Conduction and Displacement Currents 10 2.4 EM Wave Propagation 11 2.5 Reflection and Transmission Coefficients 15 2.6 Complex and Effective Parameters 17 3 Rock Physics 22 3.1 Reasons for Dielectric Dispersion at GPR Frequencies 22 3.2 Models for Fitting e(u) 25 3.3 Laboratory Measurements 29 3.4 Q and Q* Parameters 34 4 Estimation of Q* from GPR Data 3 7 4.1 Geometrical Assumptions 37 TABLE OF CONTENTS iv 4.2 Scattering Attenuation 38 4.3 Reflection and Transmission Effects 42 4.4 The Convolution Model 47 4.5 The Frequency Shift Method 51 4.6 The Wavelet Transform 57 4.7 Synthetic Tests 61 5 Inverse Q F i l t e r i n g 73 5.1 Constant Q Wave Propagation 73 5.2 Inverse Q Filter Design 77 5.3 Practical Implementation 80 5.4 Synthetic Examples : 83 6 A p p l i c a t i o n to F i e l d D a t a 90 6.1 The Langley Data Set 90 6.2 Estimation of Q* 92 6.3 Inverse Q Filtering 96 7 Conclusions 103 7.1 Summary 103 7.2 Future Work and Recommendations 104 References 105 V List of Tables 6.1 Determination of Q* from various traces in Figure 6.1 ; . . 99 VI List of Figures 3.1 Dipolar relaxation of pure water at 10 °C 24 3.2 Attenuation vs. frequency for a variety of sands 30 3.3 Attenuation vs. frequency for a variety of rocks 30 3.4 Attenuation and velocity vs. frequency for a variety of geological materials fitted using the Cole-Cole formula 32 3.5 Attenuation and velocity vs. frequency for a variety of rocks fitted using the Jonscher parameterization 33 4.1 Reflection coefficient modulus and phase vs. frequency for a number of subsurface interfaces 45 4.2 Two-way transmission coefficient modulus and phase vs. frequency for a number of subsurface interfaces 46 4.3 Synthetic GPR data used to test the frequency shift method 63 4.4 Wavelets used to create and analyze the synthetic trace in Figure 4.3c 64 4.5 Wavelet transform representation of the synthetic trace in Figure 4.3c 66 4.6 Centroid frequency versus time curve for the synthetic trace in Figure 4.3c . . . 68 4.7 Standard deviation versus time curve for the synthetic trace in Figure 4.3c . . . 69 4.8 Centroid frequency versus time curves for the synthetic trace in Figure 4.3c after inverse Q filtering to remove wavelet dispersion 71 5.1 Attenuation and velocity vs. frequency for various values of Q 76 5.2 Inverse Q filter applied to a noise-free, attenuated, reflectivity series 85 5.3 Inverse Q filter applied to the attenuated reflectivity series in Figure 5.2b with added Gaussian random noise 87 5.4 Inverse Q filter applied to the noisy, attenuated, reflectivity series in Figure 5.3b using incorrect values for Q 89 LIST OF FIGURES vii 6.1 Langley 100 MHz GPR data before correcting for wavelet dispersion 93 6.2 Centroid frequency versus time curve for trace number 300 in Figure 6.1 . . . . 97 6.3 Standard deviation versus time curve for trace number 300 in Figure 6.1 . . . . 98 6.4 Langley 100 MHz GPR data after correcting for wavelet dispersion by inverse Q filtering 101 6.5 Centroid frequency versus time curve for trace number 300 in Figure 6.4 . . . . 102 Acknowledgements Vlll Foremost, I would like to thank my supervisor, Dr. Rosemary Knight, for her tremendous guidance, generosity, encouragement, and confidence in my abilities over the past three years. Her helpful insight and constant enthusiasm for research have made my time at UBC both enjoyable and academically fruitful. I would also like to sincerely thank professor Tad Ulrych for introducing me to, and greatly inspiring me in, the areas of signal processing and time series analysis. I admire very much his philosophy of both science and life, which he has generously shared in many interesting discussions. Many thanks should also be extended to professors Rosemary Knight, Tad Ulrych, and Matt Yedlin for critically reviewing my thesis at such a busy time in their schedules. Finally, I would like to thank my family, and especially my wife, Melita, for their love, support, and encouragement during the preparation of this thesis. With very few complaints, Melita has spent many hours reading parts of this document when I was too tired to think, and has done just about everything around the apartment for the past six months. I owe you. This research was funded in part by an NSERC PGS-A scholarship, and by the U.S. Department of Energy (grant number DE-FG07-96ER14711). 1 1 In troduct ion 1.1 The G P R Method Ground penetrating radar (GPR) is a geophysical technique that employs radio waves, typically in the 10 to 1000 MHz frequency range, to map structure and features in the subsurface (Annan, 1992). Although the basic idea of using radio waves to investigate the subsurface is by no means a new one, the successful application of this idea to a wide range of geological and engineering problems is still in its infancy. In the 1960's, work with GPR began in the form of radio echo soundings (RES) to map the thickness of alpine glaciers and continental ice sheets in the Arctic and Antarctic (e.g., Waite, 1966; Rinker et al., 1966). The early 1970's saw the introduction of the technique into non-glacial environments. Initially, this work focused on applications in easily penetrated geological materials such as permafrost soil (Annan, 1976), dry rock salt (Unterberger, 1978), and coal (Coon et al., 1981). However, as the strengths and weaknesses of the method became better understood and as technology improved, GPR gained increasing popularity as an effective tool for high resolution imaging of the shallow subsurface in ordinary rock and soil. Today, GPR has a wide range of uses in many different disciplines including archaeology, civil engineering, Quaternary geology, hydrogeology, and glaciology. New applications for the technique are arising all the time. In its most common mode of operation, GPR is very similar in principle to seismic reflection profiling. First, a short electromagnetic (EM) pulse (referred to as a wavelet) is radiated vertically into the ground by a transmitter antenna. This wavelet propagates through the ground and is partially reflected back towards the surface at interfaces across which electrical properties change. The reflected wavelets are then picked up by a receiver antenna (usually oriented parallel to the transmitter, and separated by a minimal distance) 1 INTRODUCTION 2 and recorded as voltages in time by a central processing module. This process is repeated many times over as the transmitter and receiver antennas are moved along a survey line. The result is a collection of traces containing recorded voltages which can be put together and viewed just like a seismic reflection section. In contrast to seismic reflection profiling, however, GPR provides us with an image of the subsurface based on changes in its electrical, not its elastic, properties. 1.2 The Motivat ing Problem In the GPR frequency range, the attenuation of EM waves in many geological materials is strongly dependent upon frequency; that is, high frequency waves are attenuated much more quickly than lower frequency ones during propagation. To a much lesser extent, the velocity of EM waves is also slightly frequency-dependent in the GPR range. As a result, the GPR wavelet (which can be expressed as a linear combination of sinusoidal waves of various frequencies) often undergoes a significant change in shape as it propagates through the subsurface, and reflections received at later times contain noticeably less high frequency information than those received at earlier times. This phenomenon is known as wavelet dispersion. In the GPR image, it is displayed as a characteristic "blurriness" or lack of resolution that increases with depth. Correcting for wavelet dispersion in GPR data is important for a number of reasons. Qualitatively, it is desirable to have a high resolution, well-focused GPR image for the purpose of data interpretation. Quantitatively, the removal of wavelet dispersion is a necessary step before signal processing methods based on the assumption of a stationary GPR wavelet (i.e., a wavelet whose shape does not change with time down a GPR trace) can be successfully applied. Such methods include migration, whose purpose is to reposition reflection events in the GPR image to their true locations in space, and wavelet or "spiking" 1 INTRODUCTION 3 deconvolution, designed to remove the effects of the finite width of the GPR pulse. Indeed, past attempts at applying such techniques to GPR data without first correcting for wavelet dispersion have proven largely unsuccessful (e.g., Payan & Kunt, 1982; LaFleche et al., 1991; Powers, 1995). Further, since these signal processing methods lead to the recovery of the earth's reflectivity, it can be seen that correcting for wavelet dispersion is a necessary step before GPR data can be used for the quantitative determination of subsurface electrical properties. These issues will become even more important in the future as GPR systems possessing higher dynamic ranges and broader bandwidth antennas are introduced. Surprisingly, relatively little work has been done on the problem of correcting for wavelet dispersion in GPR data. The only paper devoted to this topic appears to be that of Turner (1994), who presents an algorithm for subsurface radar "propagation deconvolution" to remove the effects of frequency-dependent attenuation and velocity in the GPR range. There are three significant limitations to Turner's approach that we will attempt to overcome in this thesis. First, to correct for wavelet dispersion, Turner's algorithm requires detailed knowledge of the behaviour of the EM wave attenuation parameter with frequency in the earth. However, Turner provides little insight into what this behaviour should be and, more importantly, how it can be obtained in typical field situations. Secondly, the phase used in Turner's propagation deconvolution filter to correct for velocity dispersion is a digital minimum phase calculated from the discrete Hilbert transform of the logarithm of the filter's amplitude spectrum (given by the attenuation behaviour). Although the earth filter satisfies the minimum phase condition for continuous time, it is well known that significant differences can exist between a filter's true minimum phase spectrum and this digital minimum phase, which depends on the sampling interval of the data (Calvert et al., 1987; Varela et al., 1993). Finally, Turner's algorithm assumes that the behaviour of attenuation with frequency is constant with depth in the subsurface. It is more than likely, however, that different earth materials will possess significantly different attenuation-frequency relationships. 1 INTRODUCTION 4 1.3 Thesis Overview It is evident from the above discussion that wavelet dispersion is a common and significant problem in GPR data. The objective of this thesis is to (i) investigate in more detail the specific causes of wavelet dispersion in GPR data, (ii) develop a possible means of quantitatively estimating wavelet dispersion in GPR data, and (iii) introduce a method for removing the estimated wavelet dispersion. Chapter 2 of this thesis reviews some basic concepts in electromagnetic theory relevant to GPR and, specifically, the wavelet dispersion problem. Chapter 3 then examines some of the specific reasons for wavelet dispersion in GPR data from a rock physics perspective. Building on the results of other researchers, we will also see in this chapter that the change in shape that occurs in a GPR wavelet as it propagates through a homogeneous material can often be well described using a single parameter, Q*. Assuming that all subsurface materials can be characterized by some Q* value, Chapter 4 presents a technique for the estimation of wavelet dispersion in GPR data based on a method developed for seismic attenuation tomography. Chapter 5 then introduces a technique for removing the estimated wavelet dispersion known as inverse Q filtering. Finally, in Chapter 6 of this thesis, we present the application of these methods to a field GPR data set collected near Langley, British Columbia. 5 2 Electromagnetic Theory In order for us to understand and correct for wavelet dispersion in GPR data, it is first necessary to review some basic concepts in electromagnetic theory. This chapter begins by deriving the EM wave equations from Maxwell's equations and three constitutive relationships. From the plane wave solutions to the wave equations, we will then develop expressions for the propagation parameters that control wavelet dispersion: attenuation and velocity. Reflection and transmission coefficients will also be discussed. Finally, we will consider the complex and frequency-dependent nature of material properties, and introduce some effective parameters. 2.1 Maxwell 's Equations We begin our review of the necessary EM theory with Maxwell's equations, which describe the macroscopic behaviour of the electric and magnetic fields in terms of their sources. For an excellent treatment of classical electrodynamics leading to the derivation and study of Maxwell's equations, see Griffiths (1981). The four equations are: V - E = - (2.1) V B = 0 (2.2) V x E = - f (2.3) V x B = / i 0 J + Moe0-^ (2.4) where E and B are the electric and magnetic field vectors, respectively, p is the volume charge density, J is the current density vector, and e0 and /J0 are the dielectric permittivity and magnetic permeability of free space, respectively. The first equation (2.1) is the differential form of Gauss's law, which states that the total flux of the electric field through any closed 2 ELECTROMAGNETIC THEORY 6 surface is equal to the charge enclosed within that surface divided by CQ. Maxwell's second equation (2.2) states that the divergence of the magnetic field is always zero, or equivalently, that the total flux of the magnetic field through any closed surface must always be zero. This implies that, unlike with electric fields, there exist no stationary sources of magnetic fields; that is, "magnetic charges" do not exist. Maxwell's third equation (2.3) is the differential form of Faraday's law, which asserts that the electromotive force induced in a closed loop is equal to the negative time rate of change of magnetic flux through the loop. Finally, the fourth equation (2.4) is Ampere's law with a correction term added by Maxwell. It states that the magnetic field around a closed loop is dependent on both the current passing through the loop and the time rate of change of electric flux through the loop. It is important to note that, whereas (2.1) and (2.2) discuss the electric and magnetic fields separately, (2.3) and (2.4) show how they are fundamentally interdependent; a changing magnetic field induces an electric field, and vice versa. Since GPR involves the passage of EM waves through materials that are subject to charge separation in the presence of an electric field and slight magnetization in a magnetic field (induced electric and magnetic polarization, respectively), it is more convenient for us to work with the following versions of Maxwell's equations: V • D = P f V - B = 0 dB V x E V x H = J c + dt dB dt (2.5) (2.6) (2.7) (2.8) where D and H are the electric displacement and magnetic auxiliary fields, respectively, pf is the volume free charge density, and J c is the conduction or free current density. 2 ELECTROMAGNETIC THEORY 7 D and H are defined as follows: D = e 0 E + P H = — B - M (2.10) (2.9) where P is the electric polarization of the material or its electric dipole moment per unit volume, and M is the magnetic polarization of the material or its magnetic dipole moment per unit volume. Although equations (2.5) through (2.8) are just as general as (2.1) through (2.4), they are better suited to our purposes because they deal only with free charge and current densities. Materials that undergo electric and magnetic polarization also contain bound charge and current densities, the contributions from which are implicitly contained within the p and J terms in (2.1) and (2.4). Generally, it is much easier to know about and control the free charge and current density terms p/ and J C than it is to deal with the free and bound quantities together. 2.2 Constitutive Relationships In examining Maxwell's equations, note that they say nothing about material properties such as the electrical conductivity a, dielectric permittivity e, and magnetic permeability JJL. These parameters only take on meaning after three constitutive relationships have been introduced. Before this is done, however, we will make an assumption about linearity. In many materials, the amount of charge separation that occurs in the presence of an electric field is actually proportional to the field, provided that E is not too strong. Similarly, the amount of induced magnetization that occurs is often proportional to the magnetic field within a certain range. Materials that satisfy these relationships are called linear media. 2 ELECTROMAGNETIC THEORY 8 Following convention, the relationships are given by: P = eoXeE (2.11) M = — . Xm • B (2.12) A«o (Xm + 1) where Xe and Xm are constants called the electric susceptibility and magnetic susceptibility, respectively. In GPR, it is generally assumed that we deal with weak enough electric and magnetic fields to allow us to be within the range of linear behaviour for most geological materials (Powers, 1995; Knoll, 1996; Olhoeft, 1998). We can now introduce the three constitutive relationships as follows: J C = CTE (2.13) D = eE (2.14) B = pH (2.15) where e = e 0(l + Xe) (2.16) fi = fio(l + Xm) (2.17) The first relationship is Ohm's law, which asserts that the conduction or free current density is proportional to the electric field through the conductivity. From this we see that the conductivity measures the ease with which free charges flow through a material. The second relationship above was obtained by substituting (2.11) into (2.9). It states that the displace-ment field is proportional to the electric field through the dielectric permittivity, which is defined as the measure of a material's ability to support charge separation. Similarly, the 2 ELECTROMAGNETIC THEORY 9 third constitutive relationship was obtained by substituting (2.12) into (2.10). It states that the magnetic field is proportional to the auxiliary field through the magnetic permeability, which measures a material's ability to support induced magnetization. The three parameters cr, e, and fj, are called the constitutive parameters. We will see shortly that they control the nature of EM wave propagation within materials. Note that when the electric and magnetic susceptibilities of a material are zero, its dielectric permittivity and magnetic permeability parameters become those of free space. It should also be noted that the dielectric permittivity and magnetic permeability of a material are usually expressed relative to their free space values as follows: e r = - (2.18) lir = ± (2-19) where er and \xT are the relative parameters. The relative permittivity er is also commonly referred to as the dielectric constant, denoted by K. In general, the constitutive parameters a, e, and /J, are not real, constant values. They are complex, tensor quantities that vary with frequency and position within a medium. To simplify things, however, we will make two assumptions in this study. The first is that, in GPR, we are dealing with subsurface materials that are homogeneous at the measurement scale. In other words, we will assume that the wavelengths we are using to probe the ground are big enough to effectively "see" subsurface geological materials as homogeneous. This means that our constitutive parameters within a medium no longer vary with position. The second assumption that we will make is that the subsurface materials we deal with in GPR are isotropic; that is, that their electromagnetic properties do not vary with direction. When this is the case, the tensor nature of the quantities disappears. We are thus left with 2 ELECTROMAGNETIC THEORY 10 a, e, and [i being complex and frequency-dependent values. For the moment, however, we will ignore this and treat the constitutive parameters as real and constant. The lossy and dispersive nature of these quantities will be introduced and discussed in Section 2.6. 2.3 Conduction and Displacement Currents Looking at Maxwell's fourth equation (2.8), note that the curl of the auxiliary field H is determined by the conduction current density J C and the time derivative of the electrical displacement dD/dt. Since the latter term has the same units as JC, it is also known as the displacement current density 3d. We thus have: V x H = J C + Jd (2.20) where, using constitutive relationships (2.13) and (2.14), we get: J C = crE (2.21) <9E J D = e— (2.22) The above shows that, whereas the conduction current density is proportional to the electric field through the conductivity, the displacement current density is proportional to the time derivative of the electric field through the dielectric permittivity. Physically, conduction currents are those that exist in the presence of an electric field due to the movement of free charges such as ions and electrons. These currents are energy dissipation mechanisms. In moving through a material, free charges generate heat as they collide with other molecules; this results in an irrecoverable loss of energy. Displacement currents, on the other hand, only exist in the presence of a changing electric field and have no true physical significance. Part of the term in (2.20) accounts for contributions to 2 ELECTROMAGNETIC THEORY 11 the magnetic field that are generated directly by changes in the electric field. The rest accounts for contributions that are caused by the movement or adjustment of bound charges in response to a changing electric field (polarization currents). In contrast to conduction currents, displacement currents are energy storage mechanisms. Although work is done to move bound charges and to change the electric field, it is recoverable; the energy is stored in the material and in the field. The fundamental principle behind the use of GPR for high-resolution subsurface imaging is that the term must dominate the J C term in Maxwell's fourth equation. When this is the case, energy storage mechanisms dominate over energy loss mechanisms, and EM wave propagation can occur. When the J C term is dominant in (2.20), on the other hand, energy loss mechanisms prevail and the result is wave diffusion. 2.4 E M Wave Propagation As discussed previously, Maxwell's third and fourth equations tell us that time-varying magnetic fields produce electric fields, and that time-varying electric fields generate magnetic fields. This physical coupling between E and B has very far-reaching conclusions; in fact, it predicts the existence of electromagnetic waves. When we change the electric field at some point, for example, this produces a change in the magnetic field through the displacement current term in Maxwell's fourth equation. This change in the magnetic field then in turn generates a change in the electric field through Faraday's law, and the process repeats. In this way, energy can move through a material (or free space) as an EM wave by repeatedly changing form between the electric and magnetic fields. We will now show mathematically, from Maxwell's equations, that the electric and magnetic fields behave as waves. Substituting constitutive relations (2.13) to (2.15) into Maxwell's equations (2.5) to (2.8), and assuming that there is no net free charge density in our medium, we obtain the 2 ELECTROMAGNETIC THEORY 12 following: V - E = 0 (2.23) V • B = 0 (2.24) (9B V x E = - — (2.25) dE V x B = iioE + pe— (2.26) C/ V Taking the curl of the last two equations above, and making use of the identity V x V x A = V(V • A) - V 2 A (2.27) yields V(V • E) - V 2 E = -ix jt [oE + e^j (2.28) V(V • B) - V 2 B = (aB + e ^ ) (2.29) Equations (2.23) and (2.24) can now be used to simplify (2.28) and (2.29). We get: These are the vector wave equations for the electric and magnetic fields. There are an infinite number of solutions for E and B that satisfy (2.30) and (2.31) above. One solution, of course, is the rather complicated waveform that propagates through the subsurface in GPR. Fortunately, Fourier theory allows us to decompose any waveform that satisfies these equations into a linear combination of sinusoidal waves of various frequencies. 2 ELECTROMAGNETIC THEORY 13 Therefore, by studying the propagation of simple sinusoidal waves, we can learn about the propagation of more complicated waveforms, such as the GPR wavelet, in terms of their frequency components. We will thus begin by studying the propagation of monochromatic, EM plane waves. In complex notation, the solutions to (2.30) and (2.31) for linearly polarized plane waves traveling in the z-direction are: where E and B are the complex electric and magnetic field vectors as a function of position and time, E 0 and B 0 are the complex amplitudes containing the wave polarization and phase, u> is the angular frequency, and k is the wavenumber. The physical electric and magnetic fields E and B can be found simply by taking the real parts of (2.32) and (2.33). Substituting the solutions above into their respective wave equations, we arrive at: E(z, t) = E0ei{ut-kz) B(z,t) = B 0 e i { u t - k z ) (2.32) (2.33) k2 = ixeuJ1 — ifiau (2.34) Taking the square root of this yields the following expression for k: k — ft — ia (2.35) where (2.36) (2.37) 2 ELECTROMAGNETIC THEORY 14 Evidently, the wavenumber is a complex quantity that depends on the frequency of the waves and the constitutive parameters of the medium through which they travel. We will now use this expression for k to determine some important propagation parameters. Substituting (2.35) back into the plane wave solutions (2.32) and (2.33), we obtain: E(z,t) = E 0 e - a V ( < l r t - / J * ) (2.38) B(z,t) = B0e-azei^t~^ (2.39) This shows that the imaginary part of the wavenumber is responsible for an exponential decrease in the amplitude of the plane waves as they propagate. For this reason, a is called the attenuation. The real part of the wavenumber fl, on the other hand, can be seen to represent the spatial frequency of the waves. It is usually referred to as the phase parameter. The skin depth of the waves, defined as the distance to which they can travel within the medium before they are attenuated to 1/e of their original amplitude, is easily determined from a using: 6=- (2.40) a Similarly, the velocity v and wavelength A of the waves are obtained from ft as follows: v = ^ (2.41) A = ^ (2.42) In GPR, we are usually dealing with high enough frequencies and materials having low enough conductivities such that the low-loss condition a/cue < 1 is satisfied. This can be shown to correspond to the case where displacement currents (energy storage mechanisms) 2 ELECTROMAGNETIC THEORY 15 dominate over conduction currents (energy loss mechanisms) (Annan, 1992). When this relationship is valid, a binomial expansion of the inner square root term in (2.36) and (2.37) can be used to significantly simplify the propagation parameters above (Balanis, 1989). In particular, the expressions for the attenuation and velocity simplify to: (2.43) v ~ —— (2.44) It is important to note that these simplified expressions are not explicitly dependent upon frequency. However, we have already stated that a is very frequency-dependent in GPR (this is the main cause of wavelet dispersion). Clearly, this frequency-dependence must be largely due to the complex and frequency-dependent nature of the parameters cr, e, and \i. It should also be noted that the behaviours of the attenuation and velocity parameters with frequency are related through the Hilbert transform (Aki & Richards, 1980): ~ = — + n where Voo is the limit of v(u) as u —> oo, and % denotes the Hilbert transform operator. This equation is simply a statement of the Kramers-Kronig dispersion relations between the real and imaginary parts of a medium's index of refraction, and can be derived solely from the condition of causality without recourse to a specific wave equation. Knowing the attenuation behaviour of a medium, it allows us to calculate the velocity dispersion, and vice versa. 2.5 Reflection and Transmission Coefficients When an EM wave encounters an interface between two materials having different electric or magnetic properties, a portion of the wave is reflected and the other portion is transmitted. a(to) (2.45) 2 ELECTROMAGNETIC THEORY 16 The ratio of the complex amplitude of the reflected wave to that of the incident wave is called the reflection coefficient for the interface. Similarly, the ratio of the complex amplitude of the transmitted wave to that of the incident wave is called the transmission coefficient for the interface. We thus have: * = | (2.46) T = § (2.47) No-where R and T are the reflection and transmission coefficients, respectively. Note that these parameters are defined in terms of the amplitudes of the electric, not the magnetic, field. In general, solving for R and T for a particular interface is very complicated; their values depend on the geometry of the interface, the polarization, frequency, and angle of the incident wave, and the electric and magnetic properties of the two materials comprising the interface (Powers, 1995). This procedure can be simplified, however, if we make a few key assumptions. Assuming that (i) we are dealing with an interface that is smooth and planar at the wavelength scale, (ii) the incident wave is planar, uniform, and linearly polarized, and (iii) the incident wave's electric field is polarized perpendicular to the plane of incidence (referred to as T E mode), then R and T for the interface are given by: R _ iMth cos Bj - Hi y7k-i - kx2 sin2 Q{ ^ fi2ki cos di + p,\ y/k22 — k2 sin2 d{ 2^kx cos 9i fi2ki cos 6i + ii\ \Jk22 — k\ sin2 Q{ where Q{ is the angle of incidence measured from the interface normal, and the subscripts 1 and 2 denote the materials on either side of the interface containing the incident and transmitted waves, respectively. 2 ELECTROMAGNETIC THEORY 17 Equations (2.48) and (2.49) are known as the Fresnel equations, and can be rather lengthily derived from the boundary conditions for the electric and magnetic fields (e.g., Stratton, 1941; Griffiths, 1981). For normal incidence (t9j = 0), they simplify to: R = ^ ~ Mlf2 (2.50) T = , 2 M l , (2.51) The Fresnel equations demonstrate that any contrast in either a, e, or ii across an interface will generate a reflected wave (R ^ 0). Since the wavenumber A; is a complex quantity, they also show that the reflection and transmission coefficients for EM waves are actually complex; that is, R and T describe not only the relative amplitudes of the reflected and transmitted waves, but also any phase changes (other than 180°) that might occur upon reflection or transmission. Finally, it should be noted that R and T are, in general, frequency-dependent quantities. This means that each frequency component of a wavelet incident on an interface will be reflected and transmitted with a different amplitude and phase shift. The significance of this observation will be discussed in Sections 4.2 and 4.3. 2.6 Complex and Effective Parameters So far, we have treated the constitutive parameters a, e, and \i as real, constant values. However, in the context of EM waves, these parameters are actually dispersive; that is, they are frequency-dependent and complex. Following convention, we thus have: <T(W) = a'(u) + io"(u) (2.52) e(w) = e'{co) - ie"(u) (2.53) = - in"(u) (2.54) 2 ELECTROMAGNETIC THEORY 18 The dispersive nature of the conductivity can be briefly explained as follows. When we are dealing with low frequency EM waves and thus slowly varying electric fields in materials, inertial effects on free charges can essentially be ignored, and the conduction current density varies in phase with the electric field. At very high frequencies, however, these inertial effects can become significant and J c can actually lag in phase behind E . To account for this in Ohm's law, we must introduce a complex and frequency-dependent conductivity. In a similar fashion, we can also explain the dispersive nature of the dielectric permittivity. In materials, there exist many types of bound charges that contribute to the net polarization P in the presence of an electric field by shifting or rotating in response to the field. At low frequencies, all of the bound charge movements that contribute to the polarization occur in phase with the oscillating electric field, and thus P varies in phase with E . At higher frequencies, however, certain polarization mechanisms can start to lag behind the electric field and cause P to become out of phase with E . To account for this phenomenon, we must introduce a complex and frequency-dependent electric susceptibility, and thus a complex and frequency-dependent dielectric permittivity. Lastly, the dispersive nature of the magnetic permeability can be explained in the very same manner as the permittivity above, except that in this case we are dealing with the alignment of bound currents and a phase lag that can occur between the total magnetization M and the magnetic field B. It should be noted that the mechanisms by which J c , P and M lag behind the driving E and B fields are termed relaxation processes. For an excellent set of physical analogues that help to describe the relaxation processes mentioned above, the reader is referred to Griffiths (1981) and Balanis (1989). In general, it is a rather tedious process to accommodate the dispersive nature of all three constitutive parameters into the equations that we derived for EM wave propagation in Section 2.4. However, in GPR, we can make two important assumptions that simplify this procedure considerably. First, we can assume that conductivity relaxation is negligible 2 ELECTROMAGNETIC THEORY 19 at GPR frequencies, and thus that the conductivity is simply equal to its value at zero frequency adc (Keller, 1987; Olhoeft k Capron, 1994; Xu & McMechan, 1997). Secondly, we can assume in GPR that the effects of magnetic relaxation are also negligible, and that the magnetic permeability for subsurface materials is just equal to its free space value Ho (Topp et al, 1980; Turner & Siggins, 1994; Hollender & Tillard, 1998). Although it has been shown that this second assumption is not valid for materials containing large amounts of ferromagnetic substances such as iron, cobalt, and nickel (Olhoeft k. Capron, 1994), large quantities of these materials (which are highly conductive) are not likely to be found in the resistive environments suitable for EM wave propagation and thus GPR. We can now incorporate the complex and frequency-dependent nature of the dielectric permittivity into our existing framework for EM wave propagation as follows. Expressing the conduction and displacement current density terms in Maxwell's fourth equation (2.20) in terms of the electric field and our constitutive parameters, we have the following: V x H = crE + e— (2.55) Acknowledging the dispersive nature of the permittivity, our assumptions that a — adc and H = Ho, and also the fact that we are dealing with time-harmonic electric fields such as those described by the plane wave solution (2.32), we get: V x H = adcE + [e'(u)-e"(ij)]{iuB) = [ a d c W V ) ] E + e'(u;)^ (2-56) This shows that the imaginary or out-of-phase component of the dielectric permittivity contributes to energy loss through the conduction current density. For this reason, we can 2 ELECTROMAGNETIC THEORY 20 treat e"(uj) as a conductivity term and thus define an effective conductivity and dielectric permittivity as follows: V x H = aef(u) E + eef(u) ^ (2.57) where aef(co) = adc + coe"(co) (2.58) ee/(w) = e'(u) (2.59) Incorporating the complex and frequency-dependent nature of the permittivity into our existing framework for EM wave propagation is now straightforward. Any a and e terms in Section 2.4 should be replaced with the real, effective parameters above. The real and imaginary parts of the dielectric permittivity are not independent quantities. In fact, like the attenuation and velocity parameters, causality dictates that they be closely related through the Hilbert transform, once again satisfying the Kramers-Kronig dispersion relations (Jonscher, 1977; Bano, 1996): c ' M - C o o = ri[-e"(cj)} (2.60) e » = W ^ M - C o o ] (2.61) where denotes the high frequency limiting value of the real part of the permittivity. Also, a parameter commonly seen in the GPR literature is the loss tangent, which can now be defined as: tmS=^L (2.62) 2 ELECTROMAGNETIC THEORY 21 In order for GPR to serve as an effective tool for high-resolution subsurface imaging, the loss tangent must be less than 1. As mentioned previously, this corresponds to the case where displacement currents (energy storage mechanisms) dominate over conduction currents (energy loss mechanisms). Finally, it should be noted that, because we can generally ignore the complex and frequency-dependent nature of both the conductivity and magnetic permeability at GPR frequencies, the strong frequency dependence exhibited by the attenuation, and thus wavelet dispersion in GPR data, must largely result from the complex and frequency-dependent nature of the dielectric permittivity. This will be explored further in Chapter 3. 22 3 Rock Physics In the last chapter, we concluded that wavelet dispersion in GPR data is largely a result of the complex and frequency-dependent nature of the dielectric permittivity. This chapter begins by discussing some of the specific reasons for dielectric dispersion in geological materials at GPR frequencies. Then, two empirical models for fitting complex permittivity measurements will be discussed. Using laboratory data from different studies fitted with these two models, we will then demonstrate that the distortion that occurs in a GPR pulse as it travels through a homogeneous material can often be quantified using one parameter. 3.1 Reasons for Dielectric Dispersion at G P R Frequencies In geological materials, many different polarization mechanisms contribute to the complex and frequency-dependent nature of the dielectric permittivity. However, in the GPR fre-quency range from 10 to 1000 MHz, only two of these mechanisms are commonly cited as being important; these are interfacial polarization and the dipolar polarization of water molecules (Olhoeft & Capron, 1994; Powers, 1995). Interfacial polarization processes, otherwise known as Maxwell-Wagner mechanisms (Maxwell, 1891; Wagner, 1914), become possible when two materials having different conduc-tivities are combined. When such a composite is in the presence of an electric field, charged particles move easily through the more conductive material, but they become impeded and accumulate at boundaries with the more resistive medium. Consequently, such substances can undergo much higher degrees of polarization than their individual constituents. Indeed, "anomalously" high values for the measured dielectric permittivity at low frequencies are well known in rock physics research (Sen, 1981). As an example, consider the case of a resistive rock matrix saturated with a conductive pore fluid such as brine. When this rock is subjected to an electric field, ions in the pore fluid begin to move, but many are prevented 3 ROCK PHYSICS 23 from moving freely throughout the rock because of the geometry of the pore space. As a result, numerous ions become trapped at various locations within the rock and contribute to its polarization. This interfacial polarization mechanism resulting from ionic conductivity has been shown to be an important contributor to dielectric dispersion at GPR frequencies (Kenyon, 1984; Garrouch & Sharma, 1994). Other examples of interfacial processes that contribute to dielectric dispersion in the GPR frequency range are electrochemical double layer polarization, which occurs in moist, clay-bearing materials (Olhoeft & Capron, 1994), and surface effects seen in fully and partially-saturated sandstones whereby the presence of a thin, conductive surface layer coating the pore space has a large impact on the dielectric response (Knight & Nur, 1987). Note that, in order for the above mechanisms to contribute to dielectric dispersion at GPR frequencies, we must have relaxation; that is, the polarization processes must become unable to keep up with the oscillating electric field. It appears as if interfacial mechanisms relax at the low end of the GPR frequency range (Sherman, 1988). Olhoeft (1994) notes that frequency dependence due to interfacial polarization should only occur below 300 MHz. The relaxation of dipolar polarization is also commonly cited as an important contributor to dielectric dispersion at GPR frequencies (Hoekstra & Delaney, 1974; Carcione, 1996; Davis & Annan, 1989). This type of polarization occurs because water molecules, which are dipolar species, tend to undergo a common rotational alignment in the presence of an electric field. Most dipolar relaxation occurs from 1 to 100 GHz (Sherman, 1988), which is outside the GPR range, and thus most of the dielectric dispersion resulting from this mechanism occurs within this band. However, there is still enough dipolar relaxation from 100 to 1000 MHz to cause a significant variation in the attenuation with frequency over the GPR range. Figure 3.1 shows the dielectric permittivity, effective conductivity, and plane wave attenuation parameter versus frequency for pure water at 10 °C. Note that, although the real and imaginary parts of the permittivity change very little with frequency in the 3 ROCK PHYSICS 24 0 200 400 600 800 1000 "0 200 400 600 800 1000 frequency (MHz) frequency (MHz) Figure 3.1: Dipolar relaxation of pure water at 10 °C calculated using data from Hasted (1973). Relaxation mechanisms at higher frequencies have been ignored. The relative dielectric permittivity, effective conductivity, and attenuation are shown on a linear scale in the GPR frequency range. GPR range, just the presence of an imaginary part in this range due to water relaxation causes the effective conductivity to vary quite significantly with frequency above 100 MHz (this is clearly shown in the expression for the effective conductivity, aef = adc + toe"). As a result, the attenuation for pure water is extremely frequency-dependent in the upper part of the GPR band. 3 ROCK PHYSICS 25 Ground penetrating radar operates in a frequency range that is generally between the peaks of the interfacial and dipolar relaxation processes; the interfacial relaxation peak occurs at lower frequencies, whereas the dipolar peak occurs at higher frequencies. As a result, there is a region in the GPR frequency band between these two very dispersive, high-loss zones where the real and imaginary parts of the dielectric permittivity for geological materials exhibit very little variation with frequency (Olhoeft & Capron, 1994; Powers, 1995). This region is known as the high frequency window for electromagnetic wave propagation (Olhoeft, 1984). However, as shown above, just the presence of an imaginary component to the permittivity in this region is enough to cause the attenuation to be strongly dependent upon frequency. Therefore, it can be concluded that it is more the complex, rather than the frequency-dependent, nature of the dielectric permittivity that is the cause of wavelet dispersion in GPR. Finally, it should be noted that all of the polarization mechanisms mentioned here rely on the presence of water in subsurface materials. Indeed, most dry geological materials exhibit little dielectric dispersion at GPR frequencies (Hoekstra Sz Delaney, 1974; Sherman, 1988; Olhoeft & Capron, 1994). In reality, however, we are almost always dealing with materials containing some water in GPR, so dielectric dispersion (and therefore wavelet distortion during propagation) is usually present. 3.2 Models for Fi t t ing e(u) In order to conveniently describe the behaviour of the dielectric permittivity with frequency for various materials, researchers fit their laboratory measurements using a number of different empirical models. Two of the most commonly used and successful models in rock physics research are the Cole-Cole formula (Cole & Cole, 1941) and the power law model (Jonscher, 1977). The Cole-Cole formula is a generalization of the Debye relaxation equation 3 ROCK PHYSICS 26 (Debye, 1945), which describes the complex dielectric permittivity for a system containing a single polarization mechanism. The Debye equation states the following: e » - z e » = £ o o + (3.1) 1 + VJJT where e'(w) and e"(co) are the real and imaginary parts of the permittivity, respectively, Coo and es are the high and low frequency limiting values of the real part, ui is the angular frequency in radians, and r is a relaxation time constant for the system. Balanis (1989) derives this equation using a mechanical analogy for electric polarization consisting of a mass-spring-friction assemblage subject to an oscillating force field. In fact, the Debye equation applies to any physical system classified as an overdamped harmonic oscillator. A typical Debye response for the dielectric permittivity was seen in Figure 3.1, where we examined the dipolar relaxation of pure water. At low frequencies, the polarization varies in phase with the electric field, and thus there is no imaginary component to the permittivity. As frequencies increase, however, the polarization begins to lag in phase behind the field (relaxation begins), and the permittivity becomes complex. At some point, the imaginary component of the permittivity reaches a maximum value, which corresponds to the maximum amount of energy being lost from the system per cycle. The frequency at which this occurs is called the relaxation frequency. Eventually, at very high frequencies, the electric field oscillates far too rapidly for the polarization mechanism to operate at all. When this is the case, the real and imaginary contributions from the mechanism become zero. The Debye equation can be used to adequately describe the relaxation of any of the individual polarization mechanisms that were discussed in the last section (Carcione, 1996). In geological materials, however, more than one type of polarization process occurs and contributes to the measured permittivity. To account for this, we could sum over a number of Debye expressions having different relaxation time constants (Xu & McMechan, 1997), 3 ROCK PHYSICS 27 but fitting experimental data to a function of this form requires far too many parameters. Instead, the empirical Cole-Cole formula accounts for the effects of numerous polarization mechanisms through the addition of a distribution term over the time constants in the Debye equation (Cole k, Cole, 1941). That is, £ < M + T ^ - j ^ (3.2) where a is the time constant distribution parameter that varies between 0 and 1. The Cole-Cole formula has been used with great success in the GPR frequency range to fit dielectric permittivity measurements on geological materials within experimental errors (Taherian et al, 1990; Olhoeft & Capron, 1993). It can be seen that only four fitting parameters are required: e^ , es, r, and a. Also note that the Cole-Cole formula is a generalization of the Debye equation; when a = 1, the two expressions are the same. Another commonly used and successful model for fitting permittivity measurements versus frequency is the power law model. This empirical model stems from the work of Jonscher (1977), who noticed that the post-relaxation-peak behaviour of the imaginary component of the permittivity for most solid dielectrics can be described by the following power law relationship: e"(to) oc a / 1 - 1 where n < 1 (3.3) Jonscher also noted that, in many cases, n ~ 1; that is, the imaginary or lossy part of the permittivity is often very flat over a wide range of frequencies. An interesting property of (3.3) is that the Hilbert transform of such a power law function also has the same dependence upon frequency. Since, as discussed earlier, the real and imaginary parts of the dielectric permittivity are related through the Hilbert transform (see equations (2.60) and (2.61)), 3 ROCK PHYSICS 28 this led Jonscher to conclude that: e'(u) - oc un~l (3.4) In other words, away from a relaxation peak, the real and imaginary parts of the permittivity should have the same power law dependence upon frequency. Jonscher called this behaviour the "universal" dielectric response. Furthermore, he showed that the ratio between the real and imaginary parts is related to the exponent n as follows: „ , '— = cot — 3.5 Jonscher's results have been used in many different ways to fit dielectric permittivity measurements versus frequency. A particularly useful power law model that has been very successful in fitting data in the GPR frequency range is derived from the above equations by Hollender and Tillard (1998). They note that if we define eref — e'(u>) — C Q Q when u = uiref (3.6) where ujrej is an arbitrary reference frequency, then the constant of proportionality in (3.4) can be determined, and the relationship becomes: e'(w) - Coo = <W ( ) (3.7) \Wref/ Combining this with (3.5), the complex dielectric permittivity can now be expressed as: UJ n-l CJ, ref, e'(u;) - = eref [ — ) 1 - cot (^ — ) I + eoo (3.8) 3 ROCK PHYSICS 29 Hollender and Tillard call this formulation the Jonscher parameterization. It requires only three fitting parameters: ere/, n, and e^ . The Cole-Cole and power law models have both been shown to fit permittivity data in the GPR frequency range with a high degree of accuracy (Taherian et al, 1990). This is very interesting, considering that the forms of equations (3.2) and (3.8) are so markedly different. Knight and Nur (1987) explain this phenomenon by noting that, away from the relaxation frequency, the Cole-Cole formula reduces to an expression showing a power law dependence upon frequency of both e'(u) and e"(to). Since, as mentioned previously, the GPR frequency range tends to lie between the relaxation peaks of the interfacial and dipolar polarization mechanisms, it therefore makes sense that permittivity measurements in this range should be fitted well with either model. 3.3 Laboratory Measurements In order to estimate and correct for wavelet dispersion in GPR data, we must first learn from laboratory measurements how the wave propagation parameters that control wavelet dispersion behave with frequency in geological materials. As mentioned previously, these parameters are the attenuation and the velocity. Figures 3.2 and 3.3 show laboratory measurements by Turner and Siggins (1994) of attenuation versus frequency for a wide variety of sands and rocks in the GPR frequency range. The parameter Q* on the plots will be discussed in Section 3.4. Note that, although there is a wide variation in the behaviour of the attenuation between materials, the attenuation is generally very frequency-dependent. This confirms our previous statement that wavelet dispersion is a common and significant problem in GPR data. More important, however, is the nature of the attenuation-frequency relationship shown on these plots. Present commercial GPR systems radiate wavelets with a bandwidth of approximately two octaves (Hollender & Tillard, 1998). Over any two octave ROCK PHYSICS Figure 3.3: Attenuation properties of a variety of rocks as a function of frequency. From Turner and Siggins (1994). 3 ROCK PHYSICS 31 range, Figures 3.2 and 3.3 show that the attenuation for many geological materials can be approximated by a linear function of frequency. Noting that velocity dispersion can be calcu-lated from the attenuation behaviour using the Hilbert transform relationship (2.45), Turner and Siggins make an important conclusion from this result. They suggest that the change in shape that occurs in a radar pulse as it propagates through a homogeneous material can often be well described using one parameter. This parameter is related to the slope of the material's attenuation versus frequency curve over the bandwidth of the pulse. To investigate this further, the attenuation and velocity parameters have been calculated using equations (2.36), (2.37), and (2.41) for an additional number of geological materials whose dielectric properties have been measured in the GPR frequency range. Note that the effective parameters given by (2.58) and (2.59) were used for the calculations. Unfortunately, due to a significant lack of published data at GPR frequencies, the choice of materials was quite limited. Figure 3.4 shows the attenuation and velocity versus frequency for a collection of soils and rocks whose permittivity measurements have been fitted using the Cole-Cole formula. Figure 3.5 shows these parameters for a variety of rocks whose permittivities have been fitted using the Jonscher parameterization (Hollender & Tillard, 1998). Both figures are very informative. Indeed, over most two octave frequency ranges, it can again be seen that the attenuation is generally very frequency-dependent, and that this dependence can be approximated by a linear relationship. Thus, in accordance with Turner and Siggins, it seems reasonable that we should be able to describe using one parameter the change in shape that occurs in a GPR wavelet as it propagates through a homogeneous material. Figures 3.4 and 3.5 also show that the velocity for many geological materials is approximately constant (increasing very slowly with frequency) over typical wavelet bandwidths. This suggests that the change in wavelet shape due to velocity dispersion alone will often be minimal in GPR. 3 ROCK PHYSICS 32 o I 1 1 i i i i i i i I 0 100 200 300 400 500 600 700 800 900 1000 frequency (MHz) Figure 3.4: Attenuation and velocity calculated as a function of frequency for a variety of geological materials whose dielectric permittivities have been fitted using the Cole-Cole formula. (1) sandy soil (natural state), (2) sandy soil (saturated) (Olhoeft & Capron, 1993); (3) clay soil (dry), (4) clay soil (30.18 wt% water) (Olhoeft & Capron, 1994); (5) and (6) limestone saturated with brine, (7) and (8) sandstone saturated with brine (Taherian et al., 1990). 3 ROCK PHYSICS 33 o i i i i i i i i i i I 0 100 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 1 0 0 0 frequency (MHz) Figure 3.5: Attenuation and velocity calculated as a function of frequency for a variety of rocks whose dielectric permittivities have been fitted using the Jonscher parameteri-zation (Hollender & Tillard, 1998). (1) limestone (freshwater saturated) (Coutanceau, 1989); (2) andesite, (3) shale, (4) gabbro, (5) siltstone (Turner, 1993); (6) granite (freshwater saturated), (7) granite (dry), (8) schist (measured perpendicular to the schistocity) (Tillard, 1994). 3 ROCK PHYSICS 34 3.4 Q and Q * Parameters The quality factor, Q, is a measure of the maximum energy stored to the energy dissipated per cycle in a propagating wave. Q is given by the relation (Stacey et al, 1975): co where co is the angular frequency, v is the velocity, and a is the attenuation. In seismic studies, Q has been found to be largely frequency-independent over a wide frequency range (Strick, 1967; Stacey et al., 1975; Kjartansson, 1979). This is known as "constant Q", and it results because seismic wave attenuation is nearly linear with frequency and velocity dispersion is minimal over a broad bandwidth for many materials. Consequently, seismic Q is sometimes expressed as (Turner & Siggins, 1994; Quan & Harris, 1997): « = s b ( 3- 1 0 ) where VQ is the approximately constant value for the velocity, and s is the slope of a best-fit line through the attenuation-frequency curve, constrained to extrapolate through the origin. Q is a well known parameter in seismic studies. When combined with (2.45), the Hilbert transform causality relationship, it completely describes the change in shape of a seismic pulse as it propagates through a material. Values for seismic Q tend to lie between 50 and 300 (Sheriff, 1984), a lower value representing greater wavelet distortion during propagation. For electromagnetic waves, Q is closely related to the loss tangent, given by (2.62). This is most easily shown using equations (2.43) and (2.44), the low-loss approximations for the attenuation and velocity parameters, respectively. From these we have: wtan<5 a = -^r (3-n) 3 ROCK PHYSICS 35 and thus Q can be defined as (von Hippel, 1962): Q = l (3.12) tan S Unfortunately, constant Q has no real significance in GPR. A number of studies have shown that the loss tangent, and thus Q, are rarely constant over the bandwidth of a GPR wavelet, let alone over a wide frequency range like we have in seismology (Xu k, McMechan, 1997; Hollender & Tillard, 1998). The reason for this lies in the nature of the attenuation-frequency relationship in GPR. Although we have shown that, over the bandwidth of a GPR wavelet, attenuation is roughly linear with frequency and velocity dispersion is small, the line approximating the attenuation in GPR seldom passes directly through the origin (see Figures 3.4 and 3.5). According to (3.9), this is necessary for constant Q. It should be noted that Bano (1996) derives a constant Q relationship for GPR from the power law work of Jonscher (1977). However, in order to do this, he neglects the important contribution to the dielectric permittivity, and also the zero frequency conductivity o^. In order to describe the change in shape of a GPR wavelet as it propagates through a homogeneous material using one parameter, and also to conveniently compare wavelet dispersion in GPR with that in seismology, Turner and Siggins (1994) define a new parameter, Q*. This parameter is obtained from the slope of a best-fit line through the material's attenuation-frequency curve over the bandwidth of the wavelet as follows: (3.13) In the frequency region of interest, the attenuation can now be expressed as: (3.14) 3 ROCK PHYSICS 36 where a0 is the intercept of the best-fit line at to = 0. Note that Q* is a generalization of the Q parameter; when ao = 0, Q* = Q. Although Q* has no significance in terms of energy like Q, Turner and Siggins show that the change in shape of a pulse described by a particular value of Q* is the same as that described by the same value of Q. The only difference between a constant Q* and a constant Q response is in the total amplitude. Therefore, we can compare Q* directly with Q to evaluate wavelet dispersion. Using their laboratory results, Turner and Siggins also suggest that values for Q* in GPR tend to lie between 2 and 30 (see Figures 3.2 and 3.3). This range is significantly lower than the range of values given for seismic waves, which indicates that wavelet dispersion is far more severe in GPR than in seismology. It is important to note that, in the above discussion, we have expressed both constant Q and Q* using the slope of an attenuation-frequency curve and a constant velocity. In doing this, we may have implied that these parameters completely describe wavelet distortion during propagation by assuming that velocity dispersion does not occur, and that attenuation is strictly linear with frequency. This is not the case, and is actually a violation of causality (Futterman, 1962; Aki & Richards, 1980). Although both Q and Q* can be determined from the slope of a best-fit line to an attenuation-frequency curve and a reference velocity, they describe a much more complex scenario whereby attenuation is nearly linear with frequency and velocity is nearly constant. This will be explored further in Chapter 5. Assuming that the attenuation and velocity behaviour of all subsurface materials can be characterized using some Q* value, our problem of estimating wavelet dispersion in GPR data now becomes one of determining Q* in the subsurface. Chapter 4 is devoted to the development of a method for estimating general values for subsurface Q* from GPR data. Chapter 5 focuses on the removal of wavelet dispersion from GPR data using a technique known as inverse Q filtering. 37 4 Estimation of Q* from G P R Data In the last chapter, we concluded that the change in shape that occurs in a GPR wavelet as it travels through a homogeneous material can often be well described using one parameter, Q*. This parameter is related to the slope of the material's attenuation-frequency curve over the wavelet's bandwidth. Assuming that all subsurface materials can be characterized in this manner, this chapter develops a possible means of estimating general values for Q* in the subsurface from GPR data. First, we will make some geometrical assumptions. Then, the effects of other mechanisms that can produce distortion in a GPR wavelet (aside from frequency-dependent attenuation and velocity) will be evaluated. Using these results, a convolution model for a GPR trace will be introduced. Finally, we will develop from this model a possible method for estimating Q* using the wavelet transform. 4.1 Geometrical Assumptions In this chapter, we will make two important geometrical assumptions. First, we will assume that the traces that we record in GPR are effectively zero-offset (i.e., that the transmitter and receiver antennas can be considered to be coincident). This assumption is common in GPR research (e.g., Davis & Annan, 1989; Turner, 1994), and is especially valid at later times where wavelet dispersion is significant. This is because signals received at later times come from greater depths in the subsurface where the error in this geometrical approximation is negligible. In a 200 MHz GPR survey, for example, a commonly used antenna separation is 0.5 m (Annan, 1992). From personal experience, dispersion is rarely seen in 200 MHz data at times corresponding to depths above 4 m. At 4 m, the ratio of antenna separation to depth is 0.125, which is insignificant. Note that the assumption of zero-offset traces implies that all reflections reaching the receiver occur at normal incidence. This in turn implies that all reflected waves received have the same polarization as the incident wave. 4 ESTIMATION OF Q* FROM GPR DATA 38 The second geometrical assumption that we will make in this chapter is that we are dealing with linearly polarized, EM plane waves in GPR. In other words, we will assume that the propagating GPR wavelet can be expressed as a linear combination of monochromatic plane waves given by (2.32) and (2.33). This assumption is common to nearly all of GPR research. Away from the immediate vicinity of the transmitter antenna (in the radiating far-field region), it is well known that radar waves can be approximated by plane waves (Turner, 1994; Powers, 1995). Wavelet dispersion is significant only in reflections received from this region. Also, most GPR antennas are designed to create linearly polarized waves (Powers, 1995). Because the transmitter and receiver antennas are parallel in GPR data collection, the assumptions made above imply that the signal that we record in GPR represents the actual amplitude of the reflected electric wavefield, and not just the projection of this wavefield onto some direction. This is an important simplification that is necessary for the development of a convolution model. It allows us to ignore wave polarization and treat the electric field vector as a scalar quantity in terms of its amplitude. 4.2 Scattering Attenuation The attenuation that we observe on a GPR trace is actually a result of two mechanisms. Intrinsic attenuation, which is the only mechanism that we have discussed so far, results as EM wave energy is absorbed by subsurface materials and converted into heat during propagation. It clearly involves energy loss. Scattering attenuation, on the other hand, occurs as a result of subsurface geometry. It does not involve a loss of energy, but merely consists of the redistribution of energy to later arrival times and other propagation directions. There are many types of scattering that contribute to scattering attenuation. These include geometrical spreading, volume scattering, reflection and transmission, scattering at rough 4 ESTIMATION OF Q* FROM GPR DATA 39 interfaces, and multipathing. Since our goal is to estimate and correct for wavelet dispersion in GPR data, we are interested only in those types of scattering that are frequency-dependent; that is, we want to know about, and evaluate the effects of, scattering mechanisms that can cause a change in wavelet shape and therefore corrupt our estimation of subsurface Q*. A number of different types of scattering and their effects on the GPR wavelet are briefly discussed below. Although we have noted that radar waves can be approximated by plane waves in the far-field region of the transmitter antenna, we cannot ignore the fact that a great deal of amplitude loss down a GPR trace is caused by geometrical spreading. This type of "scattering" occurs because the energy leaving the transmitter is actually spread over the area of an ever-expanding surface (the antenna radiation pattern) during propagation. In a constant velocity medium, the amplitude of a GPR wavelet decays roughly as 1/r due to geometrical spreading, where r is the distance traveled from the transmitter. In realistic subsurface situations where velocity changes with depth, this spreading relationship is much more complicated (Powers, 1995). Because the velocity of EM waves in the GPR range is slightly dependent upon frequency, geometrical spreading is also slightly frequency-dependent, and can therefore, in theory, cause wavelet dispersion. In accordance with most GPR researchers, however, we will assume that this frequency dependence is negligible, and that geometrical spreading results in uniform amplitude losses during propagation. Volume scattering occurs when the wavelength of electromagnetic energy propagating through a material approaches the size of spatial heterogeneities within that material. When this is the case, energy is scattered randomly by the heterogeneities. Note that, in Section 2.2, we assumed that the wavelengths used in GPR are big enough to effectively "see" subsurface materials as homogeneous. Volume scattering results because this assumption is not entirely valid. Volume scattering causes a progressive decay in the amplitude of an EM wave as it propagates. More importantly, however, it is very frequency-dependent; high frequencies 4 ESTIMATION OF Q* FROM GPR DATA 40 are generally scattered much more than lower ones (Davis & Annan, 1989; Olhoeft, 1998). As a result, this scattering mechanism can potentially cause wavelet dispersion in GPR data identical to that caused by intrinsic attenuation. Fortunately, it appears as if the frequency-dependent effects of volume scattering are often negligible in the GPR range. Schaber et al. (1986) show that volume scattering only becomes significant at frequencies above 1 GHz for a dry sample of pebble channel alluvium (a very coarse material). Further, Olhoeft and Capron (1994) state that the intrinsic attenuation effects of only a few weight percent water within a material will usually dominate over the dispersive effects of volume scattering at GPR frequencies. For these reasons, we will assume that volume scattering has negligible frequency dependence in GPR. In the case where frequency-dependent attenuation due to volume scattering cannot be neglected, the general values for Q* that we estimate from our GPR data will include a contribution from this mechanism. Reflection and transmission of the GPR wavelet at subsurface interfaces make up another important contribution to scattering attenuation. Much of the decay in amplitude down a GPR trace is a result of transmission losses, which occur because signals reflected from deep interfaces encounter more transmission coefficients (and therefore more energy partitioning) along their path from transmitter to receiver than signals reflected from shallow interfaces. In addition to causing amplitude losses down a trace, reflection and transmission can also cause distortion of the GPR wavelet. As mentioned in Section 2.5, the reflection and transmission coefficients for EM waves are, in general, complex and frequency-dependent quantities. This means that each frequency component of an incident wavelet will be reflected and transmitted with a different amplitude and phase shift at an interface (i.e., the wavelet will undergo a change in shape). The significance of this distortion upon reflection and transmission will be evaluated in Section 4.3. In order to express the reflection and transmission coefficients for EM waves using the Fresnel equations (2.48) and (2.49), we assumed smooth and planar interfaces at the 4 ESTIMATION OF Q* FROM GPR DATA 41 wavelength scale. In GPR, however, this is often not the case, and amplitude losses due to random scattering at rough interfaces constitute yet another contribution to scattering attenuation. As could be expected, scattering at rough interfaces is frequency dependent; high frequencies are, in general, scattered more than lower ones (Annan & Davis, 1977). Thus, like volume scattering, this mechanism can potentially mask as frequency-dependent, intrinsic attenuation and cause a progressive broadening of received wavelets down a trace. We will assume in this study, however, that the frequency-dependent effects of rough inter-face scattering are insignificant compared to those of intrinsic attenuation. Once again, if this assumption is not valid, then our estimates of subsurface Q* will include a contribution from this scattering mechanism. The last type of scattering that we will consider is multipathing. Unlike the mechanisms discussed above, multipathing does not contribute to scattering attenuation by altering the shape or amplitude of a propagating GPR pulse. Instead, it causes apparent wavelet dis-tortion on a GPR trace through the interference of reflection events. For our purposes, multipathing can be defined as the process by which two or more reflected signals that have traveled via different paths in the subsurface reach the receiver antenna at approximately the same time. When this occurs, the signals combine constructively and destructively to produce unusual waveform shapes that can be mistaken for frequency-dependent material properties (Olhoeft, 1998). In seismic studies, multipathing has been shown to cause strong, progressive, signal broadening down a trace in geological environments that contain cyclic sequences of very thin layers (O'Doherty & Anstey, 1971; Schoenberger & Levin, 1974). This can result in estimates of Q that are much too low, and therefore inverse-Q deconvolu-tions that are unsuccessful. We will ignore such effects in this study, but it should be noted that our estimates of Q* from GPR data may be significantly affected by this phenomenon. We will also ignore the effects of isolated multipathing events in this study. Under the defi-nition of multipathing given above, this includes thin layer reflections. Although Hollender 4 ESTIMATION OF Q* FROM GPR DATA 42 and Tillard (1998) show that these events can result in significant signal distortion on a GPR trace, their effects will be spurious. In other words, isolated multipathing events will produce isolated regions of anomalous distortion on a GPR trace. These regions will occur on top of the general trend of pulse broadening caused by frequency-dependent intrinsic attenuation. In Sections 4.5 and 4.6, we will introduce a means of estimating Q* from GPR data that is somewhat insensitive to these anomalous regions. 4.3 Reflection and Transmission Effects As mentioned, reflection and transmission at subsurface interfaces can cause distortion of the GPR wavelet. This is because the reflection and transmission coefficients for EM waves are, in general, complex and frequency-dependent quantities. For the special case involving interfaces between constant Q materials, Turner and Siggins (1994) show that the reflection and transmission coefficients are independent of frequency. For interfaces between constant Q* materials, however, no such simplification can be made. Consequently, we cannot always obtain accurate reflected and transmitted waveforms in GPR simply by scaling our source wavelet by some factor. Instead, we must convolve it with time series representing the reflection and transmission coefficients. Assuming a 1-D layered earth and ignoring multiple reflections, a GPR wavelet reflected from the ith interface in the subsurface will encounter one reflection coefficient during its travels, and all of the transmission coefficients in both directions for the i — 1 overlying interfaces. Since convolution is a linear operation, we can thus think of our received wavelet (distorted only by reflection and transmission effects) as being the convolution product of the source wavelet with 1 reflection coefficient and i — 1 two-way transmission coefficients. Acknowledging our assumptions of zero-offset traces and that ix — no for geological materials, the Fresnel equations for normal incidence, (2.50) and (2.51), can be used to express the 4 ESTIMATION OF Q* FROM GPR DATA 43 reflection and two-way transmission coefficients for an interface in the frequency domain as follows: R { U J ) ~ M " ) + A*M { ] From these expressions, it is apparent that the following relationship exists between the coefficients R and T*: Tt(w) = 1 -i?(u;)2 (4.3) Note that, since we assumed that scattering from rough interfaces has negligible frequency dependence in GPR, our use of the Fresnel equations to obtain the above expressions is well justified. Power losses at rough interfaces will be accounted for in the next section through the introduction of a frequency-independent factor. To evaluate the significance of wavelet distortion caused by reflection and transmission, (4.1) and (4.3) were used to calculate the reflection and two-way transmission coefficients for a number of "typical" subsurface interfaces. Once again, due to a significant lack of published permittivity data at GPR frequencies, the choice of geological materials for these interfaces was quite limited. Figure 4.1 shows the modulus and phase of the reflection coefficients as a function of frequency. Figure 4.2 is a similar plot of the two-way transmission coefficients. It should be noted that, although the interfaces selected are typical in the sense that they are geologically reasonable, they by no means represent the majority of interfaces encountered in a GPR survey; that is, high-contrast interfaces were purposely chosen in order to get a feel for the maximum amount of distortion that can occur upon reflection and transmission. The vast majority of subsurface interfaces encountered by a GPR wavelet will have much smaller 4 ESTIMATION OF Q* FROM GPR DATA 44 R values than those shown, and values will be much closer to 1. More importantly, most reflection and two-way transmission coefficients will exhibit far less frequency dependence than the coefficients shown. Figures 4.1 and 4.2 illustrate a number of important points. First, it can be seen that the reflection and two-way transmission coefficients for the selected interfaces exhibit little frequency-dependence over most of the GPR range. Also, the phase of these coefficients tends to lie within tolerable limits of 0 and 180 degrees at most frequencies. This suggests that, except at low operating frequencies, wavelet distortion caused by reflection and transmission can essentially be ignored in GPR. Secondly, the figures show that two-way transmission coefficients are significantly less frequency-dependent than their corresponding reflection coefficients. In fact, for operating frequencies above 100 MHz, the frequency dependence of is negligible for the chosen interfaces. This is important because a received wavelet is filtered by many two-way transmission coefficients, but by only one reflection coefficient. It suggests that the combined effects of transmission through many interfaces will often result in minimal wavelet distortion. Finally, Figures 4.1 and 4.2 illustrate that, at low frequencies, the moduli of R and increase for some interfaces, but decrease for others; that is, some interfaces provide a low frequency "boost" and a high frequency "cut" upon reflection and two-way transmission, respectively, whereas other interfaces do the opposite. This suggests that the frequency-dependent effects of reflection and transmission are somewhat random in the subsurface, and will possibly cancel for a received wavelet. Based on these results, we will assume that the wavelet distortion introduced by reflection and transmission in GPR is negligible compared with that caused by frequency-dependent attenuation. In the case where this distortion cannot be ignored, it will most likely result in spurious effects that will not influence our estimation of Q*. In other words, reflection and transmission will not likely cause progressive wavelet broadening down a GPR trace like intrinsic attenuation. It should be noted that Hollender and Tillard (1998) assert that 4 ESTIMATION OF Q* FROM GPR DATA 45 100 200 300 400 500 600 700 800 900 1000 frequency (MHz) 300 400 500 600 frequency (MHz) 700 800 900 1000 Figure 4.1: Reflection coefficient modulus and phase versus frequency for a number of subsurface interfaces. Reflection coefficients were calculated using published Cole-Cole and Jonscher parameters. Numbers after material descriptions refer to weight percent water saturation. (1) sandy soil (1.23%) / sandy soil (10.53%), (2) sandy soil (1.23%) / clay soil (15.77%), (3) sandy soil (10.53%) / clay soil (30.18%), (4) clay soil (15.77%) / clay soil (30.18%) (Olhoeft and Capron, 1993, 1994); (5) granite / schist, (6) andesite / gabbro (Hollender & Tillard, 1998). Note that reflection coefficients for waves traveling in the opposite direction (i.e., from the second medium to the first) have the same modulus, but differ in phase by 180°. 4 ESTIMATION OF Q* FROM GPR DATA 46 Figure 4.2: Two-way transmission coefficient modulus and phase versus frequency for the interfaces in Figure 4.1. See Figure 4.1 for interface details. 4 ESTIMATION OF Q* FROM GPR DATA 47 reflection and transmission are the cause of significant wavelet distortion in GPR data. However, to make this point, they use an example involving an interface between a low-loss and extremely high-loss material. This type of interface is the exception, rather than the norm, in GPR surveys. Further, the high-loss material used in their example was found to not even permit radar wave propagation due to the attenuation being so high. 4.4 The Convolution Model Making use of the assumptions in the last three sections, we can now represent a recorded GPR trace as a series of convolutions. Specifically, it can be thought of as the convolution of the transmitter source wavelet with the impulse responses of the earth and the recording system (Turner, 1994). That is, where tr(t) is the GPR trace, s(t) is the source wavelet, y(t) is the recording system response, and e(t) is the earth's response. Note that we have ignored the effects of noise for now. The earth's response can also be expressed as a series of convolutions. For a subsurface containing N reflection coefficients, and ignoring multiple reflections, this yields: Here, c(t) is a filter that represents source coupling between the transmitter antenna and the ground, g(zj) accounts for frequency-independent amplitude losses in signals received from depth Zj, p{zj,t) is a propagation filter that accounts for all effects of velocity and attenuation, and rj(t) is a time series representing the jfth effective reflection coefficient. Note that g(zj) accounts for geometrical spreading, volume scattering, and scattering from rough tr(t) = s(t) * y(t) * e(t) (4.4) tr(t) = s(t) * y(t) * c(t) * V ^ ) \p(zj: t) * rj(t)] (4.5) 4 ESTIMATION OF Q* FROM GPR DATA 48 interfaces, all of which we assumed in Section 4.2 to have negligible frequency-dependence. Also note thatTj(t) represents all of the reflection and transmission coefficients encountered by the GPR wavelet as it travels from the surface to depth Zj and back; that is, rj(t) is the convolution of all of these time series. Equation (4.5) can also be expressed in terms of an "effective" source wavelet w(t) in the following manner: N tr(t) = w(t) * 5>(^ ) Wv*) * r M (4-6) where w(t) = s(t) * y(t) * c(t) (4.7) This formulation illustrates two very important points. First, it shows that a GPR trace is simply a sum of wavelets reflected from different interfaces in the subsurface. We will use this property shortly to develop a means of estimating general values for subsurface Q*. Secondly, (4.6) illustrates that we can think of a GPR trace as the convolution of an effective source wavelet w(t) with the term under the summation. This term is commonly referred to as a reflectivity series. In seismic studies where spectral contributions from reflecting interfaces are minimal and wavelet distortion during propagation can often be neglected, the reflectivity series is simply a series of spikes. In other words, reflections received from subsurface interfaces have the same spectral character as the source wavelet. In GPR, however, wavelet distortion due to frequency-dependent attenuation cannot be neglected. For this reason, the reflectivity series is a series of altered spikes that accommodate the non-stationarity. 4 ESTIMATION OF Q* FROM GPR DATA 49 For our purposes, it is much more useful to express the above convolution model for a GPR trace in the frequency domain. In order to do this, we take the Fourier transform of equation (4.6), which yields: N TR(u) = W{u,)Y,9(*i) [P{zj,u>)Rj{u>)] (4-8) If we assume, as discussed in the last section, that the wavelet distortion introduced by reflection and transmission at subsurface interfaces is minimal compared with that caused by frequency-dependent attenuation, then the effective reflection coefficients RJ{UJ) can be treated as real and frequency-independent. Lumping the RJ(OJ) and g(zj) together into the real, frequency-independent quantity Gj thus gives: N TR(u) = W(u>) GjP(zj, w) (4.9) 3=1 For propagation through a single material, the transfer function P(ZJ,UJ) can be determined from the plane wave solution for the electric field given by equation (2.32). Treating the propagating GPR wavelet as a sum of monochromatic plane waves of this form and taking the Fourier transform, we can obtain: P{zj,u) = e-2ik^ (4.10) where k(u) is the complex and frequency-dependent wavenumber of the material, and the factor of two arises because the wavelet is traveling to depth Zj and back. In reality, however, the GPR wavelet does not propagate through a single material having a constant wavenumber (if this were the case, we would not receive any reflected signals from the subsurface). We must therefore generalize the above equation to account for variations in k(u) with depth. 4 ESTIMATION OF Q* FROM GPR DATA 50 That is, P(zj,cu) = exp fZj —2i / k(uj,z) Jo dz (4.11) Using (4.9) and (4.11), the complex frequency spectrum of an individual wavelet reflected from depth Zj can now be expressed as follows: RWj(u)) = GjW{u)) exp Jo 2i / k(u, z) dz (4.12) Noting that k = f3 — ia, we can determine the amplitude spectrum of an individual reflection by taking the modulus of (4.12). This gives: ,RWj(w)\ = Gj\W(u>)\ exp —2 / a(u,z) Jo dz (4.13) A reflected wavelet cannot contain any frequencies outside the frequency band of the source wavelet w(t). If we now acknowledge our assumption that the attenuation behaviour of all subsurface materials can be characterized by some Q* value over the bandwidth of the source wavelet, then we can substitute (3.14) into (4.13) to get: RWj{u)\ = Gj\W{u)\ exp f(a»( 2v(z)Q*(z) dz (4.14) Note that the exponential resulting from the first term under the integral sign above is independent of frequency. This demonstrates that the intercept of the attenuation-frequency curve cvo, which distinguishes Q* from Q, has no effect on the form of the amplitude spectrum of a reflected wavelet, and thus the wavelet's shape in the time domain. It only affects the total amplitude, as was discussed in Section 3.4. Incorporating the frequency-independent 4 ESTIMATION OF Q* FROM GPR DATA 51 exponential in (4.14) into the term Gj, we thus have: RWj(u)\ = GJ\W(UJ)\ exp Jo v{z)Q u 1 i s * (4.15) In GPR, the traces that we record are time series; that is, we receive reflected wavelets as a function of time, not depth. However, (4.15) gives the amplitude spectrum of a reflected wavelet only as a function of the depth of the reflector Zj. Also note that it involves the velocity function v(z), which is generally unknown. Remembering that the goal of this chapter is to develop a means of estimating subsurface Q* from G P R data, it would be useful to have an equation that gives the amplitude spectrum of a reflected wavelet in terms of the time of its arrival tj, and without needing to know the velocity distribution. We can exploit the following relationship to accomplish this goal: dz = ^ - d t (4.16) Substituting (4.16) into (4.15), we arrive at the final result: RWj(u)\ = Gj\W(u)\ exp rh i 7o wr* -OJ / _ „ . . . . dt (4.17) This formula will be used in the next section. 4.5 The Frequency Shift Method In seismic studies, a broad variety of methods exist for the estimation of Q in the subsurface. Of these methods, the two most widely used and accepted are the rise time and spectral ratio techniques (Tonn, 1991). Since, as discussed previously, Q and Q* describe the same change in wavelet shape that occurs during propagation, both the rise time and spectral ratio 4 ESTIMATION OF Q* FROM GPR DATA 52 techniques can, in theory, be used to estimate Q* from GPR data. However, a number of practical limitations make the use of these two methods difficult, if not impossible, in GPR. First, we will discuss briefly the rise time and spectral ratio techniques, and the reasons why they are not particularly suitable for our purposes. Then, a variation of the frequency shift method will be developed from our results in Section 4.4 as a possible means of overcoming some of their limitations. Introduced by Gladwin and Stacey (1974), the rise time technique is a time domain method that involves the determination of subsurface Q (or Q*) from the rise times of received wavelets. The rise time r is most easily defined as the time difference between the 90% and 10% amplitude levels of a wavelet's first peak (Blair & Spathis, 1982). This technique relies on the following empirical relationship between r and Q: where ro is the rise time of the source wavelet, T is the travel time, and C is an empirical "constant". Equation (4.18) demonstrates that, as a pulse propagates and broadens due to frequency-dependent attenuation, its rise time increases at a rate inversely proportional to Q. The two main problems that result when attempting to apply this technique to GPR data are the following. First, it is difficult to obtain precise and robust measurements of r from field data because the rise time parameter is very sensitive to noise (Tarif & Bourbie, 1987; Quan & Harris, 1997). Since signal-to-noise ratios in reflection GPR data are often quite low, this suggests that rise time measurements from these data could be very inaccurate. Secondly, the parameter C in (4.18) is actually not constant, but depends on the nature of the source wavelet and on the value of Q when Q < 20 (Kjartansson, 1979; Tonn, 1991). Since our effective source wavelet in GPR is generally unknown due to the complexities of antenna/ground coupling (Hollender k. Tillard, 1998), and Q* values have been suggested to (4.18) 4 ESTIMATION OF Q* FROM GPR DATA 53 lie between 2 and 30 (Turner & Siggins, 1994), this means that C is an unknown quantity in GPR, and thus determination of Q* using the rise time technique will likely involve a great deal of speculation. In contrast to the above, the spectral ratio technique operates in the frequency domain and involves the determination of Q (or Q*) from the ratio of the amplitude spectra of received wavelets to that of a source or reference wavelet. This technique has been used for years by seismologists, and is given excellent treatment in Toksoz et al. (1979) and Sears and Bonner (1981). Taking the natural logarithm of (4.17), we can formulate the spectral ratio technique for GPR data from our results in Section 4.4 as follows: Since the term under the integral sign and Gj are both frequency-independent, (4.19) shows that the natural logarithm of the amplitude spectral ratio on the left-hand side is linear with frequency. The slope of this linear relationship can, in theory, be used to determine Q*. In practice, however, a number of problems arise when attempting to apply the spectral ratio technique to GPR data. First, the method requires a complete and accurate reference spectrum | W(u>)|, which is difficult to obtain in GPR. As mentioned, our effective source wavelet w(t) is generally unknown due to the complexities of antenna/ground coupling. Secondly, the technique requires the isolation of individual reflected arrivals in time in order to determine their amplitude spectra. This can be difficult with GPR data because the density of reflections down a trace is often quite high; that is, it can be hard to tell where one reflection begins and another ends. Finally, the spectral ratio technique is very sensitive to any extra signal present in a reflection as a result of wavelet interference. This extra signal causes scalloping and nulls in the amplitude spectrum, and thus tends to corrupt estimates of Q or Q* (Tarif & Bourbie, 1987; Dasgupta & Clark, 1998). (4.19) 4 ESTIMATION OF Q* FROM GPR DATA 54 In an attempt to overcome the limitations discussed above, we propose a variation of the frequency shift method as a means of determining general values for subsurface Q* from GPR data. Originally developed by Quan and Harris (1997) for seismic attenuation tomography, this method is based on the principle that, as a wavelet propagates and broadens due to frequency-dependent attenuation, the centroid of its amplitude spectrum undergoes a gradual downshift in frequency. Quan and Harris relate this downshift to Q for a number of different source spectra. Applying their results to equation (4.17) from Section 4.4 and converting LO (in rad/s) to / (in Hz) for practical convenience, the relationship between the centroid frequency of a received wavelet and that of our source wavelet in GPR can be expressed as follows: f™<=fw-cJ!'c7mjdt (4-20) where fRW. and fw are, respectively, the centroid frequencies of the received and source wavelets, and C is a constant that depends on the nature of the source wavelet's amplitude spectrum. The parameters JRWJ and fw are defined as follows: fw~ !?\w(i)\df ( 4 ' 2 2 ) For a boxcar source spectrum with bandwidth B, Quan and Harris show that C ~ B2/12. For a triangular source spectrum, C ~ B2/18. Finally, if the source spectrum is Gaussian in shape, then C is exactly equal to the spectrum variance given by: , _ io°°U-M'\w(f)\ df w ~ sr\w(S)\dj ( 4 2 3 ) 4 ESTIMATION OF Q* FROM GPR DATA 55 It is important to note that, like the rise time and spectral ratio techniques, the frequency shift method has no dependence on the frequency-independent quantity Gj, which contains the combined effects of a number of different scattering attenuation mechanisms, and is generally unknown. This is necessary for a robust measurement of Q*. To formulate the above results for the practical estimation of Q* from GPR data, consider the case of two interfaces in the subsurface separated by a layer of constant Q* material. In this case, the difference between the centroid frequencies of the two wavelets reflected from these interfaces can be expressed using (4.20) as: where t\ and i2 refer to the arrival times of the wavelets on a trace. This is an important result. If we assume that (i) we have a high density of reflections down our recorded GPR traces (this is often the case), and (ii) the subsurface can be described by one or a small number of general Q* values with depth, then (4.24) allows us to state the following: That is, we can obtain general or average values for subsurface Q* from the slope of a best-fit line to a centroid frequency versus time curve. This curve is determined by calculating the centroid frequency of the local amplitude spectrum at each time down a trace. In the situation where we have one general value for Q* with depth, the curve will exhibit only one slope. If two regions having considerably different Q* values are present in the subsurface, then the curve should exhibit two distinct slopes, and so on. At first glance, it appears as if the frequency shift method as formulated above is subject to the same practical limitations as the spectral ratio technique. Note that both methods IRW2 ~ IRWI — — 757 (h — * i ) (4.24) (4.25) 4 ESTIMATION OF Q* FROM GPR DATA 56 operate in the frequency domain to determine Q*. However, this is not the case. Although it seems as if this method still requires complete and accurate knowledge of our source spectrum | W ( / ) | for the determination of C, Quan and Harris (1997) show that estimates of attenuation obtained using the frequency shift method are relatively insensitive to small changes in source spectrum shape; that is, a certain value of C can be used with minimal error if the spectrum corresponding to that value of C is an approximate fit to the true source spectrum. We will assume that our source spectrum in G P R is roughly Gaussian in shape (Annan, 1992, p. 110), and therefore that C = a^y. Because Quan and Harris also show that a Gaussian input spectrum remains Gaussian with the same variance during constant Q propagation, we should be able to determine o^ relatively easily from the average of the variances of received wavelet spectra down a trace. It would also seem that the frequency shift method is subject to limitations associated with the isolation of individual reflections and wavelet interference like the spectral ratio technique. However, since this method involves the determination of amplitude spectra at all times down a trace, the isolation of individual reflections is not necessary. Further, we will see in Section 4.6 that windowing in time to determine these spectra is not needed using the wavelet transform. Also, unlike the spectral ratio technique, our variation of the frequency shift method uses the benefits of averaging to determine general values for subsurface Q*. Although some spectra will obviously be inaccurate due to wavelet interference, the fact that we are taking the slope of a best fit line through the centroids of many spectra down a trace means that the effects of any errors will likely be diminished. Moreover, we will see shortly that, by using the wavelet transform, we can reduce the problems associated with wavelet interference on the calculated spectra. Finally, it should be noted that the frequency shift method as formulated above is faster and easier to implement than the rise time and spectral ratio techniques. Since it is based on the statistics of the amplitude spectra down a trace, it may also be more robust in the presence of noise (Quan & Harris, 1997). 4 ESTIMATION OF Q* FROM GPR DATA 57 As with any Q estimation procedure, the frequency shift method is most reliable when the source spectrum is broad-band and absorption is large (White, 1992; Quan & Harris, 1997). This can be seen in equation (4.20), which shows that the centroid shift is greatest for large values of C and small values of Q*. Since frequency-dependent attenuation is far more severe in GPR than in seismology, it is thus likely that we can obtain reliable average estimates of subsurface Q* using this technique. This is especially the case for higher frequency antennas, which radiate wavelets having broader bandwidths. It should be stressed that our variation of the frequency shift method is useful only for determining general values for Q* in the subsurface; a high degree of resolution for Q* with depth cannot be obtained. For this reason, borehole radar attenuation tomography methods may be more suitable for estimating subsurface Q* in situations where Q* varies significantly with depth, or when subtle changes in Q* with depth are desired for quantitative purposes. 4.6 The Wavelet Transform As mentioned, we need to calculate the local amplitude spectrum at each time down a GPR trace in order to determine Q* using our variation of the frequency shift method. In other words, we must obtain the time-frequency representation of the trace. A critical issue in time-frequency analysis is the uncertainty principle. This principle states that we cannot know exactly what spectral components exist in a signal at a given instant in time; all we can know are the time intervals in which a certain band of frequencies are present. Another way of stating the uncertainty principle is that we cannot have arbitrarily good resolution in both time and frequency simultaneously. As our resolution (or certainty) in time increases, our resolution in frequency decreases, and vice versa. With this in mind, there are essentially two methods that have been developed to obtain the frequency content of a signal locally in time; these are the windowed or short-time Fourier transform and 4 ESTIMATION OF Q* FROM GPR DATA 58 the wavelet transform (Kumar &; Foufoula-Georgiou, 1994). We will discuss both of these techniques briefly, and the reasons why we have chosen to use the wavelet transform for our analysis of GPR data. The reader is referred to Kumar and Foufoula-Georgiou (1994) for further comparison between the two techniques in a geophysical context, as well as for an excellent list of references providing additional information. Until the mid 1980's, the standard method of obtaining a time-frequency representation of a signal was to isolate sections of the signal in time using a window function, and then take the Fourier transform of each of these sections. This technique has been formally named the short-time Fourier transform (STFT), and is expressed mathematically as follows: where f(u) is the time series to be analyzed, g(u) is the window function chosen for the analysis, and gWtt — g(u — t) e~lwu is the analyzing kernel. To compute the transform, it can be seen that the window function is moved along the time series. The result is the complex frequency spectrum around the time t, subject to the resolution limits of the uncertainty principle. It is important to note from (4.26) and (4.27) that the analyzing kernels for each frequency component in the STFT have constant support; that is, gw,t(u) has the same non-zero width for all to, which is determined by the width of the window function g(u). Because of this, the STFT has fixed resolution in time and frequency. In other words, all frequency components are resolved equally in time and in frequency with the STFT, as determined by the window function. Narrow windows provide good time resolution, but poor frequency resolution. Wide windows, on the other hand, provide good frequency resolution, but poor time resolution. When the Gaussian function is used as a window, the STFT is called the Gabor transform (Gabor, 1946). (4.26) (4.27) 4 ESTIMATION OF Q* FROM GPR DATA 59 For our purposes, the STFT possesses a number of limitations, all of which are associated with the windowing process. First, to compute the transform, we must choose a single window function that will be used for the entire time-frequency analysis. In our case, this is extremely difficult because the wavelets present on a dispersive GPR trace broaden with time. That is, we want to choose a window that is broad enough to span the widths of all reflections on a trace; however, we don't want this window to incorporate unwanted features outside of these reflections. Secondly, the frequency spectra obtained using the STFT often exhibit properties associated with the window function used, along with the signal being analyzed (Pike, 1994). This adds further complication to our choice of window function (i.e., which function gives the best spectrum for the estimation of Q*l). Finally, with the STFT, it is likely that spurious high frequency events, such as those resulting from wavelet interference and noise, will have a significant impact on many spectra because these events will be included in numerous windowed segments down a trace. It would serve useful to have a time-frequency analysis tool that is relatively "insensitive" to these anomalies. Like the STFT, the wavelet transform (WT) also examines the time-frequency character of a signal by integrating that signal with a series of analyzing kernels. However, in the WT case, the kernels are scaled and shifted versions of a prototype or "mother" wavelet function, rather than windowed exponentials. Developed in the early 1980's to overcome some resolution-related problems with the STFT, the wavelet transform can be expressed as follows: where f(u) is again the time series to be analyzed, ip(u) is the mother wavelet function a > 0 (4.28) (4.29) chosen for the analysis, a is a scale parameter, and iba>t — 1S the analyzing kernel. 4 ESTIMATION OF Q* FROM GPR DATA 60 Note that the wavelet transform does not involve frequency; instead, we obtain the analysis of our signal around the time t in terms of scale, which acts as the inverse of frequency. Scaling, as a mathematical operation, either dilates or compresses a signal. When a < 1, the mother wavelet is compressed for the analysis; when a > 1, it is dilated. The ^ term in (4.28) serves to ensure that the compressed and dilated wavelets have the same total energy as the mother wavelet function (Pike, 1994). It should be noted that the choice of ib(u) for the wavelet transform is neither unique nor arbitrary. The mother wavelet must be a function with unit energy, compact support (i.e., sufficiently fast decay), and zero mean (Kumar & Foufoula-Georgiou, 1994). Having said this, there are many well known mother wavelet functions to choose from. Our choice oiib(u) will be described in Section 4.7, where we apply our Q* estimation procedure to a synthetic G P R trace. In contrast to the STFT, (4.28) and (4.29) demonstrate that the analyzing kernels for each scale in the wavelet transform do not have the same support; that is, the non-zero width of ipaj is not constant, but varies inversely with a. For this reason, the wavelet transform does not resolve all spectral components equally like the STFT. With the W T , good resolution in time but poor resolution in frequency are obtained at high frequencies (low scales), and good resolution in frequency but poor resolution in time are obtained at low frequencies (high scales). As a result, the wavelet transform has the ability to examine features of a signal locally in time and in frequency with a detail matched to their scale (Kumar & Foufoula-Georgiou, 1994). This property is especially useful for the analysis of signals that are either noisy, non-stationary, or have important features of differing sizes (Liu & Oristaglio, 1998). A dispersive G P R trace possesses all of these qualities. With the wavelet transform, we can largely avoid the previously mentioned problems associated with the S T F T in our time-frequency analysis of G P R data. For example, because windowing is not necessary with the W T , the difficulties associated with choosing a window function are not an issue. Although we must still choose a mother wavelet function for 4 ESTIMATION OF Q* FROM GPR DATA 61 analysis with the wavelet transform, it seems reasonable that the wavelet used to analyze a G P R trace should somewhat resemble the G P R wavelet; that is, our choice of mother wavelet function should be relatively easy. Also, for our purposes, the wavelet transform is probably less sensitive than the S T F T to spurious high frequency events such as those resulting from wavelet interference and noise. Because the W T analyzes high frequencies with good resolution in time, these events will not likely affect as many spectra down a trace, and will thus have less influence on our estimates of Q*. Finally, it should be noted that, although the wavelet transform appears to be limited in the sense that it only provides us with a scale parameter, we can easily translate this parameter into an approximate frequency by determining the dominant frequency of the scaled, analyzing wavelet (e.g., Pike, 1994; Liu k Oristaglio, 1998). 4.7 Synthetic Tests To test our variation of the frequency shift method, the method was applied to a synthetic G P R trace possessing a rather conservative amount of constant Q* wavelet dispersion. As discussed previously, the reliability of the frequency shift method will improve with increasing wavelet dispersion (decreasing Q*). Figure 4.3a shows the reflectivity series used to create the synthetic data. The series contains 40 reflection coefficients randomly distributed in time between 20 and 350 ns. A sampling interval of 0.8 ns (typical for a G P R trace) was used, which implies a Nyquist frequency of 625 MHz. Note that the effects of transmission at interfaces have not been included in this simple model. Also note that, for a constant velocity of 0.1 m/ns, 400 ns would represent a depth of 20 m in the subsurface. To account for wavelet dispersion, the reflectivity series was convolved point-by-point with a series of forward Q propagation filters given by equation (5.13), which were used to simulate the effect of constant Q* = 30 with depth (see Chapter 5). Figure 4.3b shows the attenuated series. 4 ESTIMATION OF Q* FROM GPR DATA 62 Note the extreme broadening of the reflection spikes with time for this conservative value of Q*. This again confirms that wavelet dispersion is a significant problem in GPR data. Also note that, as time increases, many spikes "blur" together and become indistinguishable from one another; there is a significant decrease in resolution with time/depth. To create the synthetic GPR trace shown in Figure 4.3c, the attenuated reflectivity series was subsequently convolved with a Gaussian source wavelet having a center frequency of 200 MHz. The wavelet and its amplitude spectrum are shown in Figures 4.4a and 4.4b, respectively. The Gaussian spectrum was set with -20 dB points at 80 and 320 MHz to give the wavelet an effective bandwidth of about two octaves. It should be noted that attempting to estimate Q* from the trace in Figure 4.3c using either the rise time or spectral ratio techniques would be extremely difficult; individual reflections cannot be distinguished. To estimate Q* from the synthetic trace using our variation of the frequency shift method, the trace's time-frequency representation was first computed using the wavelet transform. For the mother wavelet function, the Morlet wavelet was used, which can be approximated by (Daubechies, 1992): 'if,(t) = e _' 2 / 2cos5i (4.30) Figures 4.4c and 4.4d show the Morlet wavelet and its amplitude spectrum, respectively. The wavelet shown was scaled to have a dominant frequency of 200 MHz for comparison purposes. The Morlet wavelet was chosen for our analysis for two reasons. First, this wavelet has already been used quite successfully for the time-frequency analysis of seismic and GPR data (e.g., Pike, 1994; Liu & Oristaglio, 1998). Secondly, it can be seen that the Morlet wavelet is smooth and symmetric in the time and frequency domains, and it closely resembles the Gaussian source wavelet in Figures 4.4a and 4.4b. As mentioned previously, it seems reasonable that the wavelet we use to analyze a dispersive GPR trace should closely 4 ESTIMATION OF Q* FROM GPR DATA 63 200 time (ns) 350 400 200 time (ns) 250 300 350 400 100 150 200 time (ns) 250 300 350 400 Figure 4.3: Synthetic GPR data used to test the frequency shift method, (a) random reflectivity series; (b) reflectivity series in (a) attenuated using constant Q* = 30; (c) the synthetic GPR trace, obtained by convolving the attenuated reflectivity series in (b) with the 200 MHz Gaussian wavelet shown in Figure 4.4a. 4 ESTIMATION OF Q* FROM GPR DATA 64 600 -20 -10 0 10 20 0 200 400 600 time (ns) frequency (MHz) Figure 4.4: Wavelets used to create and analyze the synthetic trace in Figure 4.3c. (a) 200 MHz Gaussian wavelet used to create the synthetic trace; (b) amplitude spectrum of (a); (c) Morlet wavelet used to analyze the synthetic trace with the wavelet transform (200 MHz Morlet wavelet is shown for comparison); (d) amplitude spectrum of (c). 4 ESTIMATION OF Q* FROM GPR DATA 65 resemble the GPR wavelet. When time is measured in nanoseconds, equation (4.30) yields a mother wavelet with a dominant frequency of 800 MHz. As a result, we can translate the wavelet transform scale parameter a into an approximate frequency (in MHz) for our analysis using the following relationship: Figure 4.5 shows the wavelet transform time-frequency representation of our synthetic GPR trace. The transform was calculated using wavelets having dominant frequencies from 1 to 400 MHz, in order to cover the full range of frequencies contained within the GPR signal (see Figure 4.4b). Individual amplitude spectra in Figure 4.5, which can be visualized as running vertically and out of the page, were normalized to their maximum values so as to better show the trend in the dominant frequency of the GPR signal with time. The wavelet transform plot illustrates two very important points. First and foremost, the plot clearly shows the gradual downshift in the dominant frequency of the GPR signal with time as a result of wavelet dispersion. At the beginning of the trace, it can be seen that the dominant frequency is around 200 MHz (the center frequency of the source wavelet); near the end of the trace, the dominant frequency is significantly lower (around 100 MHz). Secondly, the wavelet transform plot demonstrates the resolution limits of the uncertainty principle; high frequencies can be seen to be resolved well in time but poorly in frequency, whereas low frequencies are clearly resolved well in frequency but poorly in time. To determine a general value for subsurface Q*, the centroid and variance values for each local amplitude spectrum in Figure 4.5 were computed using equations (4.22) and (4.23), respectively. Figure 4.6 shows the centroid frequency versus time curve for our synthetic trace. Because the trace contains a high density of reflections, and because we have only one general value for Q* with depth, the centroid frequency can be seen to decrease 4 ESTIMATION OF Q* FROM GPR DATA 66 0 50 100 150 200 250 300 350 400 time (ns) Figure 4.5: Wavelet transform time-frequency representation of the synthetic trace in Figure 4.3c. Individual spectra at each time have been normalized to their maximum values in order to better show the downward trend in the dominant frequency of the GPR signal with time. 4 ESTIMATION OF Q* FROM GPR DATA 67 roughly linearly with time. This trend was fitted using least squares from 25 to 340 ns (the region of the trace containing reflections) with a line having slope -0.374 MHz/ns. It should be noted that, in some cases, fitting the centroid frequency versus time curve by eye may be more appropriate than using least squares. In regions where no reflections are present, for example, the centroid frequency tends to plunge towards zero. This can be seen at the beginning and ends of our trace, as well as briefly around 190 ns. These drops in the centroid frequency will greatly influence a best-fit line, and should therefore be ignored. It is also important to note that, based on the centroid frequency versus time curve for this simple example, fitting a large number of values for Q* with depth would be very difficult using our version of the frequency shift method; the curve clearly exhibits a large amount of variation, and it suggests that our estimation procedure will probably only work for a maximum of two or three values for Q* down a trace. Figure 4.7 shows the standard deviation versus time curve for our synthetic example. The standard deviation (square root of equation (4.23)) is plotted, rather than the variance, in order to better illustrate the spread of each local amplitude spectrum (in MHz) about the centroid frequency. As mentioned previously, Quan and Harris (1997) show that a Gaussian source spectrum remains Gaussian with the same variance during constant Q propagation. As a result, the standard deviation calculated down our synthetic trace should remain roughly constant with time. Indeed, although the curve in Figure 4.7 shows a great deal of variation, it can be seen to remain around a constant value, which was estimated using the average of values between 25 and 340 ns to be 62.4 MHz. It is important to note that this estimated value for aw is slightly higher than the actual standard deviation of our source wavelet, which was calculated to be 55.9 MHz. This observation, combined with the high variability shown in Figure 4.7 and the fact that Q* is dependent upon the square of the standard deviation measurement, suggests that the determination of C in equation (4.25) may be a significant source of error in our Q* estimation procedure. This is discussed below. 4 ESTIMATION OF Q* FROM GPR DATA 68 2501 1 1 r 0 50 100 150 200 250 300 350 400 time (ns) Figure 4.6: Centroid frequency versus time curve for the synthetic trace in Figure 4.3c. Curve was determined by calculating the centroid frequency of each local amplitude spectrum in Figure 4.5. Slope of best-fit line between 25 and 340 ns (the region of the trace containing signal) is -0.374 MHz/ns. 4 ESTIMATION OF Q* FROM GPR DATA 69 100 150 200 time (ns) 250 300 350 400 Figure 4.7: Standard deviation versus time curve for the synthetic trace in Figure 4.3c. Curve was determined by calculating the standard deviation of each local amplitude spectrum in Figure 4.5 around the centroid frequency value calculated for Figure 4.6. Average of values between 25 and 340 ns (the region of the trace containing signal) is 62.4 MHz. 4 ESTIMATION OF Q* FROM GPR DATA 70 Using the slope and standard deviation values determined above, equation (4.25) was used to yield an estimated value of Q* = 32.7 for our synthetic G P R trace. This estimate is very close to the correct value of Q* = 30, which indicates that our test of the frequency shift method was successful. To gain further insight into our Q* estimation procedure, the synthetic trace was subsequently inverse Q filtered to remove wavelet dispersion using three different values for Q. Chapter 5 discusses the design and implementation of the inverse Q filter used. Figure 4.8a shows again the centroid frequency versus time curve for the unprocessed synthetic trace. Figures 4.8b, 4.8c, and 4.8d show the curves determined for the inverse Q filtered traces using Q — 30, Q = 20, and Q = 40, respectively. Note that a running average filter could be applied to these plots to show more clearly the trend in the centroid frequency with time; this is done for our field data example in Chapter 6. The plots in Figure 4.8 illustrate an extremely important point. As mentioned above, the determination of C may be a significant source of error in our Q* estimation procedure. In fact, it seems reasonable that both our centroid frequency and variance measurements should be somehow affected by the resolution limits of the uncertainty principle, as well as by noise in the data (note that the latter was not considered in our synthetic example, and will be investigated more thoroughly in Chapter 6). Any errors in these parameters lead directly to errors in the estimation of Q*. However, Figure 4.8 demonstrates that the success of an inverse Q deconvolution, and thus the success of our estimation procedure, can be roughly verified using the frequency shift method. When the value of Q used in the inverse Q filter is too low, the input trace will be over-corrected for frequency-dependent attenuation, and the filtered output will show an increase in the centroid frequency with time, as illustrated in Figure 4.8c. On the other hand, when the value of Q used in the filter is too high, the input trace will be under-corrected for frequency-dependent attenuation, and the resulting output will show a decrease in the centroid frequency with time (see Figure 4.8d). Finally, a trace filtered using the correct value for subsurface Q* should show no general trend in the centroid 4 ESTIMATION OF Q* FROM GPR DATA 71 250 250 100 200 300 400 time (ns) 200 time (ns) 400 250 250 200 time (ns) 300 400 200 time (ns) 400 Figure 4.8: Centroid frequency versus time curves for the synthetic trace in Figure 4.3c after inverse Q filtering to remove wavelet dispersion, (a) before inverse Q filtering; (b) after inverse Q filtering using the correct value of Q — 30; (c) after inverse Q filtering using Q = 20; (d) after inverse Q filtering using Q = 40. 4 ESTIMATION OF Q* FROM GPR DATA 72 frequency with time, as illustrated in Figure 4.8b; that is, the centroid frequency versus time curve should have approximately zero slope. For this reason, it can be concluded that errors in our Q* estimation procedure are relatively unimportant. If the value determined for Q* using our variation of the frequency shift method yields an inverse Q filtered trace whose centroid frequency versus time curve is not approximately flat, then that value for Q* is simply incorrect and another one should be tried. In other words, a reasonably correct value for subsurface Q* can be found by iterating. In this case, our Q* estimation procedure is still useful in that it provides a good initial guess for the iterative process. 73 5 Inverse Q F i l t e r i n g Assuming that the attenuation behaviour of all subsurface materials can be characterized using some Q* value, we developed in the last chapter a means of estimating general values for Q* in the subsurface from the shift in the centroid frequency of local amplitude spectra down a GPR trace. Presuming that we have estimated Q* reasonably correctly using this method, this chapter presents a technique for the removal of wavelet dispersion from GPR data called inverse Q filtering. First, we will introduce a model for constant Q wave propagation. Then, the results of this model will be used to design a filter for the removal of constant Q (or Q*) wavelet dispersion. Finally, the practical implementation of this filter will be discussed. 5.1 Constant Q Wave Propagation As mentioned previously, Turner and Siggins (1994) have shown that Q and Q* describe the same change in wavelet shape that occurs during propagation; the only difference between a constant Q and constant Q* response is in the total amplitude. As a result, we can remove constant Q* wavelet dispersion from GPR data using a relatively common seismic technique known as inverse Q filtering, which removes the effects of propagation through constant Q materials. As discussed in Section 3.4, constant Q does not describe a situation whereby attenuation is strictly linear with frequency and velocity is constant. By invoking causality, it can been shown that frequency-dependent attenuation must be accompanied by velocity dispersion (e.g., Futterman, 1962; Aki & Richards, 1980). Consequently, constant Q describes a more complex scenario whereby attenuation is nearly linear with frequency and velocity is nearly constant over a wide frequency range, such that a and v satisfy the Hilbert transform causality relationship given by (2.45). Thus, to truly correct for constant Q (or Q*) wavelet dispersion, we must correct for both amplitude and phase spectrum changes in our source wavelet due to frequency-dependent attenuation and velocity, respectively. 5 INVERSE Q FILTERING 74 Many models have been developed to describe constant Q wave propagation, all of which yield similar results for the attenuation and velocity parameters as a function of frequency (Kjartansson, 1979; Varela et al, 1993). In this study, we will consider a causal, linear model based on a power law formula for the wavenumber k (Strick, 1967; Kjartansson, 1979; Bickel & Natarajan, 1985), which yields any positive value of Q exactly independent of frequency. With this model, we can express k as follows: where 7 is a constant that describes the medium, w0 is an arbitrary reference frequency, and a also describes the medium but depends on our choice of w0- Here, 0 < 7 < 1, with 7 = 0 corresponding to the case of no attenuation. Writing k above in terms of its real and imaginary parts, and considering only positive frequencies , we find that: = f3 — ia Note from equations (2.41) and (3.9) that Q is equal to one-half the ratio of the real to the imaginary parts of k (i.e., Q = L3/2CX). A S a result, we can obtain: (5.1) (5.2) Q = 1 (5.3) 2tan ( f ) which is clearly independent of frequency. We can now express 7 as follows: (5.4) 1 for relatively large Q (5.5) nQ 5 INVERSE Q FILTERING 75 Using the above results, the attenuation and velocity parameters for a constant Q material (according to the power law model) can be expressed in the following manner: where v0 = a/ cos(2^) is the velocity at the reference frequency co0. Figure 5.1 shows a(yj) and v(co) plotted over the GPR frequency range for various values of Q with v0 = 0.1 m/ns and WQ = 100 MHz. Note that, except in the highly dispersive case where Q = 2, the attenuation is nearly linear with frequency (when constrained to pass through the origin) and the velocity is nearly constant over a wide range. Also note that v(u) increases gradually with frequency, the rate of increase being fastest for low Q (highly frequency-dependent attenuation). This is consistent with Figures 3.4 and 3.5, in which we also observed a slight increase in the velocity with frequency. Incidentally, it should be mentioned that, when Q is relatively large and the frequency range is restricted to some degree around wQ, substitution of (5.5) into (5.6) and a subsequent Maclaurin series expansion of the exponential can be used to obtain the commonly seen Futterman (1962) dispersion relationship (Kjartansson, 1979): Futterman derived this relationship by assuming that attenuation was exactly linear with frequency over a specified range and invoking causality. Hence, in his results, Q is a function of frequency in order to satisfy (3.9). When Q is large and the frequency range is restricted, the Futterman results converge with the power law results because Q becomes essentially frequency-independent. (5.6) (5.7) (5.8) 5 INVERSE Q FILTERING 76 Figure 5.1: Attenuation and velocity calculated as a function of frequency for various values of Q. Equations (5.6) and (5.7), which are derived from a power law model for constant Q wave propagation, were used for the calculations. In all cases, a reference velocity of 0.1 m/ns at 100 MHz was used. 5 INVERSE Q FILTERING 77 Using equations (5.6) and (5.7), we can now express the wavenumber for a constant Q material in terms of Q and a reference velocity as follows: This formula will now be used to develop an inverse Q filter. 5.2 Inverse Q Filter Design In order to develop an inverse Q filter for the removal of constant Q (or Q*) wavelet dispersion from GPR data, we must first determine the forward Q filter operator; that is, we must first determine the operator that produces constant Q wavelet dispersion on a GPR trace. As discussed in Section 4.4, all effects of velocity and attenuation (and thus wavelet dispersion) are incorporated into our convolution model through the propagation filter P(ZJ, to). We will thus begin by examining this filter for the constant Q case. Assuming that we are dealing with a single constant Q medium, or equivalently, that the subsurface can be adequately described by one general value for Q with depth, equations (4.10) and (5.9) can be used to express P(ZJ,LO) as follows: It should be noted that, using a similar transfer function derived by Kjartansson (1979) from a power law formula for the wavenumber, Turner and Siggins (1994) successfully predict the change in shape of GPR wavelets that have propagated through the subsurface in a shallow, cross-hole, tomographic experiment (except for a difference in total amplitude). This indicates that the model for constant Q wave propagation introduced in Section 5.1 is a valid one for describing wavelet dispersion in GPR data. (5.9) (5.10) 5 INVERSE Q FILTERING 78 Equation (5.10) gives the constant Q propagation filter for an arrival received from depth Zj in the subsurface. Note, however, that this formula involves two parameters that are generally unknown in the context of a GPR trace: Zj and v0. Assuming that time on a trace represents the travel time experienced by the reference frequency UJ0, we can express the filter more conveniently as follows: Our constant Q propagation filter is now completely specified down a GPR trace by two parameters: Q and a reference frequency. If we further assume that time down a trace represents the travel time experienced by the center frequency of the transmitter antenna (a reasonable assumption), then u>o is determined, and the filter is completely specified by the parameter Q. It must be stressed that the constant Q propagation filter given by (5.11) is time-varying; reflections received at different times on a trace experience different delays and exhibit different amounts of wavelet dispersion. It should also be noted that, in many cases, we will have more than one general value for Q (or Q*) with depth. As discussed in Chapter 4, our variation of the frequency shift method has the ability to resolve a small number of general values for Q* in the subsurface. When this is the case, multiple Q values can be accounted for in (5.11) using the effective Q approach (Bickel & Natarajan, 1985; Varela et al., 1993), whereby an effective value for subsurface Q is calculated for each depth Zj, which represents the attenuation effect of the variable Q structure above that depth. In terms of time down a trace, this is expressed as follows: (5.11) Qef(tj) = (5.12) where the effective value Qef(tj) is simply substituted for Q in (5.11) above. 5 INVERSE Q FILTERING 79 As mentioned, because P(tj,u>) accounts for all effects of velocity and attenuation on a G P R trace, it accounts for wavelet dispersion. However, in accounting for all effects of velocity, our propagation filter is also responsible for placing a dispersed arrival around the time tj on a trace (the travel time of the reference frequency). We desire to know the filter that just produces wavelet dispersion about the time tj] that is, we wish to obtain the filter that, when convolved with a trace possessing no wavelet dispersion, produces constant Q wavelet dispersion in an arrival at time tj. This is easily accomplished by removing a linear phase shift term from (5.11). Multiplying by e1^, we arrive at: where U[tj,to) is our desired result, the forward Q filter operator. Note that this operator was used to produce constant Q wavelet dispersion in the random reflectivity series in Figure 4.3a; reflection coefficients at different times were subjected to different filters given by (5.13). Knowing the forward Q filter operator for each time down a G P R trace, our inverse Q filtering procedure becomes straightforward. To remove the effects of constant Q (or Q*) wavelet dispersion at the time tj on a trace, we simply convolve the trace with the time domain inverse of U(tj,u), which is easily seen to be (in the frequency domain): Thus, to deconvolve an entire trace of the effects of constant Q (or Q*) wavelet dispersion, we (in theory) just calculate H(tj,u)) for each time down the trace, take the inverse Fourier transform of each of these transfer functions, and then convolve the trace with each of the resulting time series. The j th time sample of the inverse Q filtered output trace is just equal to the jth. sample of the convolution product of the input trace with the jth time series. (5.13) (5.14) 5 INVERSE Q FILTERING 80 That is, tr?*={tr i»*r-1[H(tjlu,)]). (5.15) where trm and trout are the input and output traces, respectively, and T~x denotes the inverse Fourier transform. It should be noted that this point-by-point procedure not only corrects for the effects of frequency-dependent attenuation, but also for velocity dispersion such that all frequency components in the filtered trace travel at the reference velocity v0. Since higher frequencies travel more quickly than lower ones (see Figure 5.1), this means that frequencies less than UJQ in the input trace will experience positive delays, and frequencies greater than UJQ will experience negative delays (Varela et al., 1993). 5.3 Practical Implementation In practice, a number of important issues must be considered when applying the inverse Q filtering algorithm described above to G P R data. Most importantly, it should be noted that the true inverse Q filter operator given by (5.14) is an increasing exponential in time and frequency. To correct for wavelet dispersion, high frequencies must be boosted more than lower ones, and reflections at late times must be amplified more than those at earlier times. However, the fact that H(tj,u) is an increasing exponential makes it very unstable in the presence of noise, and thus unsuitable for the inverse Q filtering of real data. To see this more clearly, note that noise will tend to dominate signal on a G P R trace (i) at high frequencies, outside or near the edge of the frequency band of the G P R wavelet, and (ii) at late times, when the signal has become significantly attenuated. Since H(tj, to) increases with both time and frequency, it can thus be seen that the amplification performed by an inverse Q filter based on equation (5.14) will be greatest when a trace is "noisiest". As a result, filtering a noisy trace using the algorithm described above will yield an output trace which "blows up". 5 INVERSE Q FILTERING 81 Obviously, this must be avoided. To make our algorithm stable in the presence of noise, we will therefore use the least squares inverse of equation (5.13) for our inverse Q filter operator, which is given by (Berkhout, 1982): Hl2(tj,u) = —— .v. v r (5.16) where * here denotes the complex conjugate, and A(tj,u) is a regularization parameter that represents the inverse square of the input trace's signal-to-noise ratio as a function of time and frequency. It should be noted that, given an estimate of the initial signal-to-noise ratio at t = 0, and assuming that the noise present in our GPR data is both white and stationary (i.e., that it possesses uniform power at all frequencies and times), we can roughly track the signal-to-noise ratio in time and frequency down a trace using the forward Q filter operator. For this reason, we will express the regularization parameter in (5.16) as follows: A(ti,u) = =• (5.17) SNo(w)2|tf(ti,o;)|2 where SN0(co) denotes the input trace's initial signal-to-noise ratio as a function of frequency, for which a reasonable estimate can be obtained with knowledge of the amplitude spectrum of the GPR source wavelet. By using the least squares inverse of equation (5.13), noise is optimally taken into account in our inverse Q filtering procedure. In other words, we avoid excessively boosting regions of an input trace's time-frequency spectrum where noise is dominant, yet we still amplify in the most accurate manner those regions where signal-to-noise ratios are high. As mentioned previously, the amplification performed by an inverse Q filter based on equation (5.14) will be greatest when a trace is "noisiest". With the least squares inverse, our filter is prevented from becoming too large in the presence of significant noise by a large regularization parameter. 5 INVERSE Q FILTERING 82 On the other hand, when signal-to-noise ratios in the input trace are high, the regularization parameter in (5.16) will be small and the least squares inverse will approach the true inverse given by (5.14). It is important to note that, in the presence of noise, we can never obtain an inverse Q filtered trace that is completely devoid of wavelet dispersion. At frequencies and times where noise is dominant, signal cannot be boosted the appropriate amount without tremendous error. However, we will see in Chapter 6 that inverse Q filtering can still remove the majority of wavelet dispersion from field GPR data, and thus significantly enhance the resolution of the GPR image. In addition to stability in the presence of noise, a number of other practical issues must be considered when applying the inverse Q filtering procedure described in Section 5.2 to GPR data. First, since the algorithm involves the transformation of Hi2(tj,ui) from the frequency domain to the time domain using the inverse Fourier transform, we must be careful to avoid excessive ringing (Gibbs phenomena) in the resulting time domain filters. Gibbs phenomena occur after the inverse Fourier transform because of sharp cutoffs in a signal's amplitude spectrum, and can be significantly reduced by tapering the amplitude spectrum using a window function. Although tapering will obviously reduce the effectiveness of our inverse Q filter operator, it was found that the application of a Gaussian window (-10 dB point set at the Nyquist frequency) to the amplitude spectrum of Hi2(tj,co) before performing the inverse Fourier transform tended to produce the best inverse Q filtered results (i.e., those results that had the best trade-off between accuracy and lack of ringing). Another practical issue that must be considered when applying our inverse Q filter is the stability of the inverse Fourier transform of (5.16). Although the earth filter is a causal and minimum phase function (Aki Sz Richards, 1980; Varela et al, 1993), it was found that the windowed, least squares inverse of (5.13) is slightly acausal, and therefore "blows up" in a numerical inverse Fourier transform algorithm. To remedy this problem, a constant, linear phase shift must be added to the phase spectrum of Hi2(tj, ui) before transforming into the time domain, 5 INVERSE Q FILTERING 83 which makes the inverse Q filter operator causal and thus stable. The number of time samples corresponding to this phase shift must then be removed after the transformation. Finally, it should be noted that, in order to derive equations (5.6) and (5.7) which were used in the derivation of our forward Q filter operator, we considered only positive frequencies (see page 74). As a result, the expression for H~i2(tj,oj) given by equation (5.16) is valid only for OJ > 0. Since, for the inverse Fourier transform, we need to know the complex spectrum from —fpf to /AT, where is the Nyquist frequency, we must thus compute Hi2(tj,oj) at negative frequencies from the complex conjugate of its values at positive frequencies. This is a valid step because our time domain inverse Q filter operator is a real function. Using the above results, we can now summarize our inverse Q filtering algorithm for practical application to GPR data as follows: 1. Compute U(tj,u>) from 0 to Nyquist frequency using equation (5.13). 2. Determine regularized Hi2(tj,co) from U(tj,u) using equation (5.16). 3. Introduce constant, linear phase shift for stable inverse Fourier transform. 4. Window amplitude spectrum to reduce ringing in the time domain filter. 5. Set negative frequencies equal to complex conjugate of positives. 6. Inverse Fourier transform. 7. Convolve time domain filter with input trace. 8. Extract jth time sample. 9. Repeat for all times tj down the input trace. 10. Merge extracted time samples together to form inverse Q filtered trace. 5.4 Synthetic Examples To illustrate some of the points discussed above, the application of our inverse Q filter is now demonstrated on a few simple, synthetic examples. Figure 5.2a shows a reflectivity 5 INVERSE Q FILTERING 84 series containing 12 reflection coefficients. As in Figure 4.3, a sampling interval of 0.8 ns was used. Figure 5.2b shows this reflectivity series after accounting for constant Q = 30 wavelet dispersion using the forward Q filter operator given by equation (5.13). No noise was added to these data. Note again the significant broadening of the reflection spikes and lack of resolution that increases with time for this conservative value of Q. It is also now appropriate to point out the effects of velocity dispersion in the data. Since higher frequencies travel more quickly than lower ones, the attenuated spikes can be seen to be asymmetric, having a sharp onset with slowly decaying tails. Lastly, it should be noted that the first arrival on the attenuated trace exhibits slight ringing due to the Gibbs effect. This is a result of transforming U(tj,uj) into the time domain. To correct for wavelet dispersion, the attenuated reflectivity series was subsequently inverse Q filtered using the algorithm given on page 83. Figure 5.2c shows the inverse Q filtered output for Q = 30 and SN0(u) = 1.0 x 109. A constant value was used for the initial signal-to-noise ratio because the "GPR source wavelet" in this case is a spike, which contains equal amounts of all frequencies. The value chosen is exceptionally high because the input data contain no noise, except for numerical errors resulting from our forward modeling to produce wavelet dispersion. Note the extremely high resolution of the inverse Q filtered result. Comparing this plot with Figure 5.2a, it can be seen that we have successfully removed all wavelet dispersion (i.e., that caused by both frequency-dependent attenuation and velocity dispersion) from the input trace. Also note, however, the difference between the amplitudes of the reflection spikes in the initial and recovered spike series. Although not particularly important because relative amplitudes remain the same, it can be seen that the amplitudes of the spikes in Figure 5.2c are significantly less than those in Figure 5.2a. This is a result of our application of a Gaussian window to Hi2(tj,u)) in order to reduce Gibbs phenomena in the inverse Q filtered trace. Slight ringing can still be seen before and after each reflection; without windowing, the amplitude of this ringing would be significant. 5 INVERSE Q FILTERING 85 1.0 CD f o E to -1.0 (a) CD T3 Q. E to 0) T3 3 Q. E to 50 50 100 150 time (ns) time (ns) 100 150 time (ns) 200 200 250 250 250 Figure 5.2: Inverse Q filter applied to a noise-free, attenuated, reflectivity series, (a) noise-free reflectivity series with no attenuation; (b) reflectivity series in (a) attenuated using constant Q — 30; (c) attenuated series in (b) after inverse Q filtering, Q = 30, SN0{u) = 1.0 x 109. 5 INVERSE Q FILTERING 86 In order to demonstrate the performance of our inverse Q filter in the presence of noise, Gaussian random noise was next added to the attenuated reflectivity series in Figure 5.2b. The variance of this noise was set equal to 1% of the trace's maximum value. Figure 5.3b shows the resulting data. To illustrate the effects of too little regularization, the noisy trace was first inverse Q filtered using Q = 30 and SN0(co) = 1.0 x 109. In other words, the initial signal-to-noise ratio was set extremely high so that the regularization parameter in (5.16) would be much too small. Figure 5.3c shows the inverse Q filtered output. Note that, as discussed previously, the output "blows up" at late times because the amplification performed by the inverse Q filter is very high when signal-to-noise ratios are low; regions of the input trace's time-frequency spectrum where noise dominates are greatly boosted. Figure 5.3d, on the other hand, shows the inverse Q filtered result for the more realistic case of SNQ(LO) = 2.0 x 103 (proper regularization). In this case, the filter is stable in the presence of noise, and resolution is considerably improved by inverse Q filtering. Reflections down the output trace are reduced greatly in width, and all of the individual spikes in the input reflectivity series (with the exception of the one at 220 ns) can be easily distinguished. It is important to note that there is a noticeable broadening of the reflections with time in Figure 5.3d, even after inverse Q filtering. As mentioned previously, it is impossible to completely remove wavelet dispersion in the presence of noise. However, when properly applied, our algorithm clearly does a good job of removing most of the wavelet dispersion from the noisy input trace. As a final example, the noisy, attenuated, reflectivity series in Figure 5.3a was inverse Q filtered using two incorrect values for Q. Figures 5.4a and 5.4b show the output for Q = 20 and Q = 40, respectively, with SN0(co) = 2.0 x 103. In the first case, the input trace has been over-corrected for frequency-dependent attenuation. Although, at first glance, this result appears to be better than the one obtained in Figure 5.3d with Q = 30 in the sense that reflections seem sharper, note the significant number of artifacts created by the over-boosting 5 INVERSE Q FILTERING 87 CD T 3 E re _1 . 0 i 1 • 1 1 1 1 0 50 100 150 200 250 time (ns) _0 4 1 1 1 i i l ' 0 50 100 150 200 250 time (ns) Figure 5.3: Inverse Q filter applied to the attenuated reflectivity series in Figure 5.2b with added Gaussian random noise, (a) input reflectivity series; (b) attenuated series plus noise; (c) noisy, attenuated series in (b) after inverse Q filtering, Q = 30, SN0(u) = 1.0 x 109; (d) same as in (c) except SN0{u) = 2.0 x 103. 5 INVERSE Q FILTERING 88 of high frequencies in the data; new "reflections" appear that are not present in the original spike series. Obviously, this is an undesirable result. In the second case where Q = 40, our noisy input trace has been under-corrected for frequency-dependent attenuation. It should be noted that this trace is quite similar in appearance to Figure 5.3d. However, it can be seen that reflections are slightly more dispersed, and the four events between 140 and 160 ns are not as well separated in time; we have not removed all possible wavelet dispersion from the data. The above results demonstrate that it is quite important to use a reasonably correct value for Q when inverse Q filtering. When Q is too low, we get the creation of artificial reflections. When Q is too high, on the other hand, all possible wavelet dispersion is not removed from the input data. 5 INVERSE Q FILTERING 89 Figure 5.4: Inverse Q filter applied to the noisy, attenuated, reflectivity series in Figure 5.3b using incorrect values for Q. (a) series in Figure 5.3b after inverse Q filtering, Q = 20, SN0(to) = 2.0 x 103; (b) same as in (a) except Q = 40. 90 6 A p p l i c a t i o n to F i e l d D a t a In the last two chapters, we have presented a possible means of estimating and correcting for wavelet dispersion in GPR data using a variation of the frequency shift method and inverse Q filtering, respectively. We have seen these two techniques successfully applied to synthetic traces. This chapter presents the true test of these methods: their application to field GPR data. First, a 100 MHz GPR data set collected near Langley, British Columbia will be briefly introduced and discussed. Then, a general value for subsurface Q* will be estimated from these data using our variation of the frequency shift method. Finally, we will attempt to remove wavelet dispersion from the data set by inverse Q filtering. 6.1 The Langley Data Set As a final test, the techniques described in Chapters 4 and 5 for the estimation and correction of wavelet dispersion in GPR data were applied to a 100 MHz field data set collected near Langley, British Columbia. The field site consists of a sand and gravel aquifer underlain by a conductive clay, which was not penetrated with the radar. The data were gathered using a PulseEKKO IV GPR system with a trace spacing of 0.2 m and a sampling interval of 0.8 ns. Preprocessing included three basic steps. First, each trace in the data set was "dewowed" using a residual median filter (Gerlitz et al, 1993) in order to remove the slowly decaying transient or "wow" following the direct air and ground arrivals. This transient, present in all raw GPR data and characterized by a sharp onset and exponential decay, is a result of the low frequency component of the transmitter pulse diffusing, rather than propagating, into the ground. That is, the lowest end of the GPR signal spectrum sees an inductive (conduction current) type response as opposed to a propagating (displacement current) type response (see Section 2.3) (Annan, 1996). Removal of the wow in GPR data is necessary in order to obtain traces with approximately zero mean, and thus make small amplitude 6 APPLICATION TO FIELD DATA 91 reflections superimposed on the inductive response more distinct. Since dewowing is a form of high pass filtering, it should be performed with great care. Removing too much of the low frequency component on a GPR trace results in the removal of important reflection signal, and can often hide the presence of wavelet dispersion. Removing too little of the low frequency component, on the other hand, results in a trace possessing too much wow, on which reflections at early times are difficult to identify. For the Langley 100 MHz data set, it was found that a residual median filter with a window length of 51 points was optimal for the removal of wow and the preservation of signal. After dewowing, the Langley data were subsequently corrected for shifts in "time-zero". During a GPR survey, there may be a gradual drift or erratic jumps in the arrival time of the direct air wave (referred to as time-zero) from trace to trace. This can occur while system electronics warm up to their ambient operating temperature, and also as a result of changes in tension on system fiber optic cables (Annan, 1996). In order to correct for this problem, the Langley data were "flattened" on the first arrival by finding the first point on each trace whose amplitude exceeded a certain threshold value, and then setting this point as t = 0. In the case where a GPR survey is conducted in a region having variable topography, the data would first be corrected for time-zero shifting, and then moved up or down accordingly to account for the differences in elevation from trace to trace. The last preprocessing step that was applied to the Langley 100 MHz data set was an approximate correction for the geometrical spreading of energy. As discussed in Section 4.2, although radar waves can be approximated by plane waves in the far-field region of the transmitter antenna, a great deal of amplitude loss down a GPR trace occurs as a result of this phenomenon. In the case where velocity is constant with depth, geometrical spreading losses vary directly with the travel distance from transmitter to receiver. Thus, to roughly correct for geometrical spreading, the Langley data were subjected to a gain whose strength varied linearly with time. 6 APPLICATION TO FIELD DATA 92 Figure 6.1 shows a section of the Langley data set from 50 to 120 m after the preprocessing described above. In an attempt to compensate for the gradual decay in amplitude down the traces due to attenuation, a standard, frequency-independent, time-varying, exponential gain was also applied to the data for this figure. An average subsurface velocity of 0.1 m/ns, determined using the common mid-point (CMP) method (Annan, 1992), was used to convert time into an approximate depth. Note that the sand and gravel aquifer can be seen to increase in depth from left to right. The underlying clay layer is plainly indicated on the GPR image as the region containing no reflections; it is too conductive to permit effective radar wave propagation. More importantly, however, note that there is a significant amount of wavelet dispersion in these data, as indicated by the lack of resolution or "blurriness" in the image that increases with depth. Evidently, attenuation in the sand and gravel aquifer is quite dependent upon frequency. Further, the gradual broadening of reflections down the image indicates that processing techniques based on the assumption of a stationary GPR wavelet (such as migration and spiking deconvolution) would be ineffective without first correcting for wavelet dispersion. Finally, it should be noted that most traces across the Langley profile exhibit a high density of reflections. This suggests that the data set is a good candidate for Q* estimation using our variation of the frequency shift method. 6.2 Estimation of Q* To estimate a general value for Q* in the sand and gravel aquifer, eight traces across the Langley profile in Figure 6.1 were analyzed using our variation of the frequency shift method. As mentioned in Chapter 4, noise will obviously have some effects on the spectral centroid and variance values we calculate down a trace using equations (4.22) and (4.23), respectively. Since we are now dealing with real GPR data, it is appropriate here to briefly discuss these effects, along with some possible methods by which we can reduce the influence of noise on 6 APPLICATION TO FIELD DATA position (m) Figure 6.1: Section of Langley 100 MHz GPR data set from 50 to 120 m before correcting for wavelet dispersion. Notice the lack of resolution or "blurriness" in the image that increases with depth. 93 6 APPLICATION TO FIELD DATA 94 our estimates of subsurface Q*. In an ideal situation, the signal and noise regions of a trace's frequency spectrum would not overlap, and signal could be separated from noise relatively easily using simple filtering techniques. In the real world, however, a large part of the GPR signal bandwidth is always occupied by noise, and we must thus deal with the effects of noise in any quantitative analysis. At early times on a GPR trace when signal-to-noise ratios are high, the influence of noise on our calculated centroid and variance values will be minimal. In other words, the centroid and variance values determined from noisy traces at early times will be approximately equal to those values that would be obtained in the noise-free case. At late times, on the other hand, when the GPR wavelet has been greatly attenuated and signal-to-noise ratios are low, noise can result in significant bias in these two parameters. Because high frequencies in the GPR source wavelet are attenuated more quickly than lower ones, for example, centroid values obtained from noisy data at late.times will tend to be higher than those obtained in the noise-free case; that is, the calculated centroid values will be biased towards higher frequencies by the noise. Similarly, variance values obtained from noisy data at late times will also tend to be higher than those acquired from a clean signal; the noisy spectra will exhibit a greater spread about the centroid frequency. For these reasons, estimates of Q* determined from noisy data using the frequency shift method as described in Chapter 4 will probably be larger than they should be (see equation (4.25)). In other words, the presence of noise will most likely result in our underestimating the effects of frequency-dependent attenuation. To reduce the influence of noise on our estimates of subsurface Q*, three techniques were used in our analysis of the Langley 100 MHz data set. First, instead of computing the wavelet transform for all frequencies from zero to the Nyquist frequency, the WT was computed only for those frequencies estimated to lie within the bandwidth of the GPR source wavelet. This has the effect of reducing as much as possible the influence of noise outside the GPR signal range. For the Langley data set, this range was estimated to lie between 6 APPLICATION TO FIELD DATA 95 20 and 180 MHz from power spectra taken at the beginning of traces where signal-to-noise ratios were high. It should be noted that, for the reasons stated in Section 4.7, the Morlet wavelet was chosen for our analysis. Secondly, instead of using equation (4.22) to estimate the centroid of the GPR signal spectrum at all times down a trace, a weighted average between this value and the frequency fmax at which the maximum spectral amplitude occurs was employed. Although by no means a robust parameter, fmax was found to be a better estimator of the true centroid of the signal spectrum at late times when signal-to-noise ratios were low. Thus, a linear weighting in time was used, whereby the computed centroid value was weighted most strongly at early times, and fmax was weighted most strongly at late times. Finally, in order to reduce the effects of noise on our estimates of the source spectrum variance aw, an average of values at early times (where signal-to-noise ratios were highest) was employed. As mentioned in Chapter 4, assuming that we are dealing with an approximately Gaussian source spectrum, aw should remain roughly constant down a trace. It should be noted that the techniques described here for reducing the influence of noise on our estimates of subsurface Q* involve only our estimation procedure. That is, methods for suppressing noise in our GPR data before the estimation of Q* and inverse Q filtering were not considered in this study, and are a definite topic for future investigation. Figure 6.2 shows the centroid frequency versus time curve for trace number 300 in Figure 6.1, one of the eight traces analyzed across the Langley profile. The curve was obtained using the weighted average method described above, and then smoothed using a 13-point running average filter in order to show more clearly the general trend in the centroid frequency with time. Note that the curve exhibits a steady downshift as a result of frequency-dependent attenuation in the sand and gravel aquifer. The centroid frequency can be seen to change from approximately 100 MHz to 50 MHz over the course of the trace. This downshift was fitted using least squares from 20 to 360 ns (the region of the trace containing reflections, but not including the direct air and ground arrivals) with a single line having 6 APPLICATION TO FIELD DATA 96 slope -0.177 MHz/ns. Figure 6.3 shows the standard deviation versus time curve for the same trace. Note here that, at late times, the standard deviation can be seen to increase due to the effects of noise discussed above. To obtain an estimate of the standard deviation of our source spectrum, an average of values at early times (between 20 and 200 ns) was used, yielding aw = 50.9 MHz. This implies a source spectrum variance of a2^ — 2591 MHz 2. Using these results for the slope and variance, an average value of Q* = 45.9 was estimated for the sand and gravel aquifer from trace number 300. Table 6.1 summarizes the results of our analysis of the eight traces across the Langley profile. The centroid frequency and standard deviation versus time curves for each of these eight traces were very similar in appearance to Figures 6.2 and 6.3, respectively. Note that the Q* values obtained lie fairly close together (between 36.2 and 46.1). This attests to the robustness of our estimation method, and also indicates that the sand and gravel aquifer can be adequately characterized using a single Q* value. It should also be noted that these values are slightly above the upper limit of Q* = 30 suggested by Turner and Siggins (1994). This indicates that the effects of frequency-dependent attenuation in GPR data may be slightly less severe than predicted in their paper, or that noise has affected our results. 6.3 Inverse Q Filtering In order to remove wavelet dispersion, the Langley data in Figure 6.1 were inverse Q filtered using Q = 43, the average of the eight values obtained above for the sand and gravel aquifer. The algorithm given in Section 5.3 was employed. An initial signal-to-noise ratio of 1.0 x 103 was found to produce the best inverse Q filtered results. It should be noted that this initial ratio is independent of frequency; making SN0(UJ) dependent upon frequency based on an estimate of the amplitude spectrum of the GPR source wavelet was found not to produce significantly better results than using a constant value. Figure 6.4 shows the inverse Q 6 APPLICATION TO FIELD DATA 97 20 1 1 1 1 1 I I : I I 0 50 100 150 200 250 300 350 400 time (ns) Figure 6.2: Centroid frequency versus time curve for trace number 300 in Figure 6.1. Curve was smoothed using a 13-point running average filter to better show trend. Slope of best-fit line between 20 and 360 ns is -0.177 MHz/ns. 6 APPLICATION TO FIELD DATA 98 Figure 6.3: Standard deviation versus time curve for trace number 300 in Figure 6.1. Average of values between 20 and 200 ns is 50.9 MHz. 6 APPLICATION TO FIELD DATA 99 Table 6.1: Determination of Q* from various traces in Figure 6.1 Trace # Position (m) Slope (MHz/ns) Variance (MHz2) Q* 2 50.2 -0.237 3122 41.3 55 60.8 -0.202 2722 42.3 105 70.8 -0.188 2753 46.1 155 80.8 -0.231 2660 36.2 197 89.2 -0.187 2709 45.5 247 99.2 -0.201 2504 39.2 300 109.8 -0.177 2591 45.9 347 119.2 -0.197 2815 44.9 filtered data set. Notice the significant increase in resolution in this image compared with Figure 6.1. The GPR image is now well-focused and the widths of reflections at depth have been noticeably reduced. There are no visible signs of wavelet dispersion in the inverse Q filtered profile. By examining the image in more detail, it can also be seen that all reflections in Figure 6.4 correspond to events in Figure 6.1, and most can be followed laterally from trace to trace. This indicates that our inverse Q filtering has not produced any significant artifacts. Finally, it should be acknowledged that the data in Figure 6.4 are in severe need of migration, as indicated by the presence of numerous diffraction hyperbolae across the profile. Now that we have corrected for wavelet dispersion, migration would be the next logical processing step in order to further increase the resolution of the GPR image; this is another topic for future investigation. As mentioned in Chapter 4, after proper removal of wavelet dispersion, the slope of a trace's centroid frequency versus time curve should theoretically be reduced to zero. Figure 6.5 shows the centroid frequency versus time curve for trace number 300 in Figure 6.4 (the same trace analyzed previously to produce Figures 6.2 and 6.3). Once again, the curve 6 APPLICATION TO FIELD DATA 100 has been smoothed using a 13-point running average filter for clarity. Indeed, it can be seen that the centroid frequency now shows no general trend in time, except for a slight decrease at late times due to the fact that our inverse Q filter avoids boosting regions of a trace's frequency spectrum that are dominated by noise (see Chapter 5). In this case, it is obvious that our variation of the frequency shift method has worked extremely well, and there is no need to solve for a more accurate value for subsurface Q* by trial and error as described in Section 4.7. In conclusion, we have successfully estimated and corrected for wavelet dispersion in the Langley 100 MHz field data set. 6 APPLICATION TO FIELD DATA 101 position (m) Figure 6.4: Langley 100 MHz GPR data after correcting for wavelet dispersion by inverse Q filtering, Q = 43. Notice the significant increase in resolution compared with Figure 6.1. 6 APPLICATION TO FIELD DATA 102 1 6 0 | 1 1 1 r 0 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0 3 5 0 4 0 0 time (ns) Figure 6.5: Centroid frequency versus time curve for trace number 300 in Figure 6.4. Curve was smoothed using a 13-point running average filter. Notice the absence of a downshift in frequency with time after inverse Q filtering. 103 7 Conclusions 7.1 Summary In summary, wavelet dispersion is a common and significant problem in GPR data that primarily results from attenuation being strongly dependent upon frequency in the GPR range. This frequency dependence can be largely attributed to the presence of an imaginary or lossy component to the dielectric permittivity. Correcting for wavelet dispersion in GPR data is an important step that should be performed before either qualitative interpretation or quantitative determination of subsurface electrical properties are attempted. Over the bandwidth of a GPR wavelet, we have seen that the attenuation of EM waves in many geological materials is approximately linear with frequency. For this reason, the change in shape that occurs in a GPR pulse as it propagates through the subsurface can often be well described using a single parameter, Q*, related to the slope of the linear region. Assuming that the attenuation behaviour of all subsurface materials can be characterized using some Q* value, the problem of estimating and correcting for wavelet dispersion in GPR data becomes one of determining Q* in the subsurface and deconvolving its effects through the use of an inverse Q filter. In Chapter 4, a method for the estimation of general values for subsurface Q* from GPR data was developed. Essentially, this method involves the determination of Q* from the shift in the centroid frequency of the GPR signal spectrum with time down a trace. The method manages to largely overcome limitations associated with other means of estimating attenuation such as the rise time and spectral ratio techniques. Once Q* has been obtained, an inverse Q filtering technique based on a causal, linear, model for constant Q wave prop-agation was introduced in Chapter 5 as a means of removing wavelet dispersion. Tests on both synthetic and field GPR data from Langley, British Columbia showed that these two 7 CONCLUSIONS 104 techniques together provide a robust and effective means of estimating and correcting for wavelet dispersion in GPR data. More importantly, our Q* estimation procedure is especially useful in that it allows us to conveniently verify the success of an inverse Q deconvolution; a trace properly corrected for wavelet dispersion should show no general trend in the centroid frequency with time. 7.2 Future Work and Recommendations A number of topics for future investigation have presented themselves throughout the course of this thesis. First, since the methodology presented here provides an effective means of correcting for wavelet dispersion in GPR data, the successful application of migration and spiking deconvolution algorithms to further enhance the resolution of the GPR image may now be possible. This in turn may allow the determination of subsurface electrical properties from the earth's reflectivity. Secondly, since both our Q* estimation and inverse Q filtering algorithms are adversely affected by noise, it would serve useful to investigate methods of suppressing noise in GPR data before attempting to correct for wavelet dispersion. Finally, it should be stressed again that our variation of the frequency shift method is capable only of providing general estimates of Q* in the subsurface; a high degree of resolution for Q* with depth cannot be obtained. For this reason, borehole radar attenuation tomography methods may provide a more accurate means of determining subsurface Q*. 105 References Aki, K., & Richards, P.G. 1980. Quantitative Seismology: Theory and Methods. Vol. 1. W.H. Freeman and Company. Annan, A.P. 1976. Impulse radar sounding in permafrost. Radio Science, 11(4), 383-394. Annan, A.P. 1992. Ground Penetrating Radar Workshop Notes. Sensors Sz Software, Inc. Annan, A.P. 1996. pulseEKKO IV User's Manual, v. 4.2. Sensors Sz Software, Inc. Annan, A.P., Sz Davis, J.L. 1977. Radar range analysis for geological materials. In: Report of Activities, Part B. Geological Survey of Canada Paper 77-1B. Balanis, C A . 1989. Advanced Engineering Electromagnetics. John Wiley Sz Sons. Bano, M. 1996. Constant dielectric losses of ground-penetrating radar waves. Geophysical Journal International, 124(1), 279-288. Berkhout, A.J. 1982. Seismic Migration: Imaging of Acoustic Energy by Wave Field Ex-trapolation, Theoretical Aspects. Developments in Solid Earth Geophysics, vol. 14A. Elsevier Scientific Publishing Company. Bickel, S.H., Sz Natarajan, R.R. 1985. Plane-wave Q deconvolution. Geophysics, 50(9), 1426-1439. Blair, D.P., Sz Spathis, A.T. 1982. Attenuation of explosion-generated pulses in rock masses. Journal of Geophysical Research, 87(5), 3885-3892. Calvert, A.J., den Boer, L.D., Hirsche, W.K., Sz Hargreaves, N.D. 1987. Q attenuation and deconvolution. Can Soc. Expl. Geophys. Recorder. Carcione, J.M. 1996. Ground-penetrating radar: Wave theory and numerical simulation in lossy anisotropic media. Geophysics, 61(6), 1664-1677. Cole, K.S., & Cole, R.H. 1941. Dispersion and absorption in dielectrics. Journal of Chemical Physics, 9, 341-351. Coon, J.B., Fowler, J.C., Sz Schafers, C.J. 1981. Experimental uses of short pulse radar in coal seams. Geophysics, 46(8), 1163-1168. Coutanceau, N. 1989. Proprietes Dielectriques de Roches Non-Argileuses dans la Bande 20-1000 MHz: Methode de Mesure sur Eprouvette de Dimension Centimetrique. Ph.D. thesis, University of Paris. Dasgupta, R., Sz Clark, R.A. 1998. Estimation of Q from surface seismic reflection data. Geophysics, 63(6), 2120-2128. REFERENCES 106 Daubechies, I. 1992. Ten Lectures on Wavelets. SIAM. Davis, J.L., Sz Annan, A.P. 1989. Ground-penetrating radar for high-resolution mapping of soil and rock stratigraphy. Geophysical Prospecting, 37(5), 531-551. Debye, P.J.W. 1945. Polar Molecules. Dover Publishing. Futterman, W.I. 1962. Dispersive body waves. Journal of Geophysical Research, 67(13), 5279-5291. Gabor, D. 1946. Theory of communications. Journal of the Institute of Electrical Engineers, 93, 429-457. Garrouch, A.A., Sz Sharma, M.M. 1994. The influence of clay content, salinity, stress, and wettability on the dielectric properties of brine-saturated rocks: 10 Hz to 10 MHz. Geophysics, 59(6), 909-917. Gerlitz, K., Knoll, M.D., Cross, G.M., Luzitano, R.D., Sz Knight, R.J. 1993. Processing ground penetrating radar data to improve resolution of near-surface targets. In: Proc, Symposium on the Application of Geophysics to Engineering and Environmental Prob-lems. April 18-22, San Diego, CA, 561-574. Gladwin, M.T., Sz Stacey, F.D. 1974. Anelastic degradation of acoustic pulses in rock. Physics of the Earth and Planetary Interiors, 8(4), 332-336. Griffiths, D.J. 1981. Introduction to Electrodynamics, Second Edition. Prentice Hall. Hasted, J.B. 1973. Aqueous Dielectrics. Chapman and Hall. Hoekstra, P., Sz Delaney, A. 1974. Dielectric properties of soils at UHF and microwave frequencies. Journal of Geophysical Research, 79(11), 1699-1708. Hollender, F., Sz Tillard, S. 1998. Modeling ground-penetrating radar wave propagation and reflection with the Jonscher parameterization. Geophysics, 63(6), 1933-1942. Jonscher, A.K. 1977. The 'universal' dielectric response. Nature, 267, 673-679. Keller, G.V. 1987. Rock and mineral properties. In: Nabighian, M.N. (ed), Electromagnetic Methods in Applied Geophysics. Society of Exploration Geophysicists. Kenyon, W.E. 1984. Texture effects on megahertz dielectric properties of calcite rock samples. Journal of Applied Physics, 55(8), 3153-3159. Kjartansson, E. 1979. Constant Q-wave propagation and attenuation. Journal of Geophysical Research, 84(B9), 4737-4748. Knight, R.J., Sz Nur, A. 1987. The dielectric constant of sandstones, 60 kHz to 4 MHz. Geophysics, 52(5), 644-654. REFERENCES 107 Knoll, M.D. 1996. A Petrophysical Basis for Ground Penetrating Radar and Very Early Time Electromagnetics: Electrical Properties of Sand-Clay Mixtures. Ph.D. thesis, University of British Columbia. Kumar, P., Sz Foufoula-Georgiou, E. 1994. Wavelet analysis in geophysics: An introduction. In: Wavelets in Geophysics. Academic Press. LaFleche, P.T., Todoeschuck, J.P., Jensen, O.G., Sz Judge, A.S. 1991. Analysis of ground-probing radar data: Predictive deconvolution. Canadian Geotechnical Journal, 28, 134-139. Liu, L . , Sz Oristaglio, M . 1998. G P R signal analysis: Instantaneous parameter estimation using the wavelet transform. In: Proceedings of the Seventh International Conference on Ground Penetrating Radar. May 27-30, Lawrence, KS, 177-182. Maxwell, J.C. 1891. A Treatise on Electricity and Magnetism. Dover Publishing. O'Doherty, R.F. , Sz Anstey, N . A . 1971. Reflections on amplitudes. Geophysical Prospecting, 19(3), 430-458. Olhoeft, G.R. 1984. Applications and limitations of ground penetrating radar. In: SEG Extended Abstracts, 54th Annual International Meeting and Exposition. December 2-6, Atlanta, G A , 147-148. Olhoeft, G.R. 1998. Electrical, magnetic, and geometric properties that determine ground penetrating radar performance. In: Proceedings of the Seventh International Conference on Ground Penetrating Radar. May 27-30, Lawrence, KS, 177-182. Olhoeft, G.R., Sz Capron, D.E. 1993. Laboratory measurements of the radio frequency electrical and magnetic properties of soils from Yuma, Arizona. U.S.G.S. Open-File Report, 93-701. Olhoeft, G.R., Sz Capron, D.E. 1994. Petrophysical causes of electromagnetic dispersion. In: Proceedings of the Fifth International Conference on Ground Penetrating Radar. June 12-16, Waterloo, Ontario, 177-182. Payan, I., Sz Kunt, N . 1982. Subsurface radar signal deconvolution. Signal Processing, 4, 249-262. Pike, C.J . 1994. Analysis of high resolution marine seismic data using the wavelet transform. In: Wavelets in Geophysics. Academic Press. Powers, M . H . 1995. Dispersive Ground Penetrating Radar Modeling in 2-D. Ph.D. thesis, Colorado School of Mines. Quan, Y . , Sz Harris, J . M . 1997. Seismic attenuation tomography using the frequency shift method. Geophysics, 62(3), 895-905. REFERENCES 108 Rinker, J.N., Evans, S., k Robin, G.Q. 1966. Radio ice-sounding techniques. In: Proceedings of the Fourth Annual Symposium on Remote Sensing of Environment. 177-182. Schaber, G.G., McCauley, J.F., Breed, C.S., k Olhoeft, G.R. 1986. Shuttle imaging radar: Physical controls on signal penetration and subsurface scattering in the Eastern Sahara. IEEE Transactions on Geoscience and Remote Sensing, GE-24(4), 603-623. Schoenberger, M., k Levin, F.K. 1974. Apparent attenuation due to intrabed multiples. Geophysics, 39(3), 278-291. Sears, E.M., k Bonner, B.P. 1981. Ultrasonic attenuation measurement by spectral ratio utilizing signal processing techniques. IEEE Transactions on Geoscience and Remote Sensing, GE-19(2), 95-99. Sen, P.N. 1981. Relation of certain geometrical features to the dielectric anomaly of rocks. Geophysics, 46(12), 1714-1720. Sheriff, R.E. 1984. Encyclopedic Dictionary of Exploration Geophysics. Society of Exploration Geophysicists. Sherman, M.M. 1988. A model for the frequency dependence of the dielectric permittivity of reservoir rocks. The Log Analyst, 29(5), 358-369. Stacey, F.D., Gladwin, M.T., McKavanagh, B., Linde, A.T., k Hastie, L.M. 1975. Anelastic damping of acoustic and seismic pulses. Geophysical Surveys, 2(2), 133-151. Stratton, J.A. 1941. Electromagnetic Theory. McGraw-Hill Book Co. Strick, E. 1967. The determination of Q, dynamic viscosity, and creep curves from wave-propagation measurements. Journal of Geophysical Research, 13, 197-218. Taherian, M.R., Kenyon, W.E., k Safinya, K.A. 1990. Measurement of dielectric response of water-saturated rocks. Geophysics, 55(12), 1530-1541. Tarif, P., k Bourbie, T. 1987. Experimental comparison between spectral ratio and rise time techniques for attenuation measurement. Geophysical Prospecting, 35(6), 668-680. Tillard, S. 1994. Radar experiments in isotropic and anisotropic geological formations (gran-ites and schists). Geophysical Prospecting, 42(6), 615-636. Toksoz, M.N., Johnston, D.H., k Timur, A. 1979. Attenuation of seismic waves in dry and saturated rocks: Laboratory measurements. Geophysics, 44(4), 681-690. Tonn, R. 1991. The determination of the seismic quality factor Q from VSP data: A comparison of different computational methods. Geophysical Prospecting, 39(1), 1-27. REFERENCES 109 Topp, G.C., Davis, J.L., & Annan, A.P. 1980. Electromagnetic determination of soil water content: Measurements in coaxial transmission lines. Water Resources Research, 16(3), 574-582. Turner, G. 1993. The Influence of Subsurface Properties on Ground Penetrating Radar Pulses. Ph.D. thesis, University of Macquarie, Australia. Turner, G. 1994. Subsurface radar propagation deconvolution. Geophysics, 59(2), 215-223. Turner, G., & Siggins, A.F. 1994. Constant Q attenuation of subsurface radar pulses. Geo-physics, 59(8), 1192-1200. Unterberger, R.R. 1978. Radar propagation in rock salt. Geophysical Prospecting, 26(2), 312-328. Varela, C.L., Rosa, A.L.R., & Ulrych, T.J. 1993. Modeling of attenuation and dispersion. Geophysics, 58(8), 1167-1173. von Hippel, A.R. 1962. Dielectrics and Waves. John Wiley & Sons, Inc. Wagner, K.W. 1914. Erklarung der dielektrischen Nachwirkungsworgange auf Grund Maxwellscher vorrstellungen. Archiv. Electrotechnik, 2, 371. Waite, A.H., Jr. 1966. International experiments in glacier sounding, 1963 and 1964. Cana-dian Journal of Earth Sciences, 3, 887-892. White, R.E. 1992. The accuracy of estimating Q from seismic data. Geophysics, 57(11), 1508-1511. Xu, T., & McMechan, G.A. 1997. GPR attenuation and its numerical simulation in 2.5 O U N T dimensions. Geovhvsics. 62(1). 403-414.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Estimation and correction of wavelet dispersion in...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Estimation and correction of wavelet dispersion in ground penetrating radar data Irving, James D. 2000
pdf
Page Metadata
Item Metadata
Title | Estimation and correction of wavelet dispersion in ground penetrating radar data |
Creator |
Irving, James D. |
Date Issued | 2000 |
Description | The attenuation of electromagnetic (EM) waves in many geological materials is strongly dependent upon frequency in the ground penetrating radar (GPR) range; high frequencies are attenuated much more quickly than lower ones during propagation. For this reason, the GPR wavelet often undergoes a significant change in shape as it travels through the subsurface, and reflections received at later times contain less high frequency information than those received at earlier times. This phenomenon is known as wavelet dispersion. In the GPR image, it is displayed as a characteristic "blurriness" that increases with depth. Correcting for wavelet dispersion in GPR data is an important signal processing step that should be performed before either qualitative interpretation or quantitative determination of subsurface electrical properties are attempted. Previous work by other researchers has shown that the EM wave attenuation parameter for many geological materials is approximately linear with frequency over the bandwidth of a GPR wavelet. Thus, the change in shape of a GPR pulse as it propagates can often be well described using one parameter, Q*, which is related to the slope of the linear region. In this thesis, we confirm and build on these results. Assuming that all subsurface materials can be characterized by some Q* value, the problem of estimating and correcting for wavelet dispersion in GPR data becomes one of determining Q* in the subsurface and deconvolving its effects through the use of an inverse Q filter. A method for the estimation of subsurface Q* from GPR data based on a technique developed for seismic attenuation tomography is presented. Essentially, Q* is determined from the downshift in the dominant frequency of the GPR wavelet with time down a trace. Once Q* has been obtained, an inverse Q filtering technique based on a causal, linear, model for constant Q wave propagation is proposed as a means of removing wavelet dispersion. Tests on field data collected near Langley, British Columbia indicate that these methods are very effective at enhancing the resolution of the GPR image. |
Extent | 8471742 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-08 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0089439 |
URI | http://hdl.handle.net/2429/10414 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2000-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2000-0214.pdf [ 8.08MB ]
- Metadata
- JSON: 831-1.0089439.json
- JSON-LD: 831-1.0089439-ld.json
- RDF/XML (Pretty): 831-1.0089439-rdf.xml
- RDF/JSON: 831-1.0089439-rdf.json
- Turtle: 831-1.0089439-turtle.txt
- N-Triples: 831-1.0089439-rdf-ntriples.txt
- Original Record: 831-1.0089439-source.json
- Full Text
- 831-1.0089439-fulltext.txt
- Citation
- 831-1.0089439.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0089439/manifest