Published online 25 February 2008
Published in Vadose Zone J 7:171-183 (2008)
DOI: 10.2136/vzj2006.0128
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: GROUND PENETRATING RADAR IN HYDROGEOPHYSICS
Measuring the Electrical Properties of Soil Using a Calibrated Ground-Coupled GPR System
C. P. Odena,*,
G. R. Olhoefta,
D. L. Wrightb and
M. H. Powersb
a Earth Science Systems, LLC, Golden, CO 80402
b U.S. Geological Survey, Denver, CO
* Corresponding author (coden{at}earthsciencesystems.com).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 1 September 2006.
 |
ABSTRACT
|
|---|
Traditional methods for estimating vadose zone soil properties using ground penetrating radar (GPR) include measuring travel time, fitting diffraction hyperbolae, and other methods exploiting geometry. Additional processing techniques for estimating soil properties are possible with properly calibrated GPR systems. Such calibration using ground-coupled antennas must account for the effects of the shallow soil on the antenna's response, because changing soil properties result in a changing antenna response. A prototype GPR system using ground-coupled antennas was calibrated using laboratory measurements and numerical simulations of the GPR components. Two methods for estimating subsurface properties that utilize the calibrated response were developed. First, a new nonlinear inversion algorithm to estimate shallow soil properties under ground-coupled antennas was evaluated. Tests with synthetic data showed that the inversion algorithm is well behaved across the allowed range of soil properties. A preliminary field test gave encouraging results, with estimated soil property uncertainties (±
) of ±1.9 and ±4.4 mS/m for the relative dielectric permittivity and the electrical conductivity, respectively. Next, a deconvolution method for estimating the properties of subsurface reflectors with known shapes (e.g., pipes or planar interfaces) was developed. This method uses scattering matrices to account for the response of subsurface reflectors. The deconvolution method was evaluated for use with noisy data using synthetic data. Results indicate that the deconvolution method requires reflected waves with a signal/noise ratio of about 10:1 or greater. When applied to field data with a signal/noise ratio of 2:1, the method was able to estimate the reflection coefficient and relative permittivity, but the large uncertainty in this estimate precluded inversion for conductivity.
Abbreviations: FDTD, finite difference time domain • GPR, ground penetrating radar • IMSP, inverse model for soil properties • TDR, time domain reflectometry
|
INTRODUCTION
|
|---|
Ground penetrating radar (GPR) is a valuable technique in vadose zone studies for many reasons. First, GPR surveys produce information with higher spatial resolution than other surface geophysical methods, which is why it is used to study the scale of heterogeneity in the vadose zone (Bowling et al., 2005; Hendrickx et al., 2001; Rea and Knight, 1998). Furthermore, GPR is used to estimate moisture content within a few wavelengths of the surface (Lambot et al., 2006, 2004a,b; Schmalholz et al., 2004; Grote et al., 2003; Huisman et al., 2003; Olhoeft, 2000) and at depth (Loeffler and Bano, 2004; Grote et al., 2002). Ground penetrating radar has also been used to delineate transport paths for water and contaminants (Wollschläger and Roth, 2004; Gish et al., 2002; Brewster et al., 1995; Sander and Olhoeft, 1994; Greenhouse et al., 1993). The electrical properties measured by GPR can be related to hydrogeologic properties such as moisture content and porosity using mixing models (Johnson and Poeter, 2005; Sneddon, 2003; Wtorek, 2003; Sihvola, 1999). The ability to measure electrical properties in three dimensions enables generation of three-dimensional hydrogeologic property maps. Ground penetrating radar is also a good tool for soil property measurements because surveys can be conducted continuously, rapidly, and inexpensively. These favorable economics make it feasible to investigate changes in the subsurface as a function of time (Kowalsky et al., 2005; Binley et al., 2002; Brewster et al., 1995). The difficulties with modeling fluid transport in the vadose zone can be addressed by conducting frequent GPR surveys to determine the location of fluids with time. Finally, GPR makes measurements noninvasively without the need to push sensors into the soil, as with time domain reflectometry (TDR) or neutron measurements.
Ground-coupled GPR antennas are often used in surveys over soils to increase both energy transfer to the ground and penetration depth over that of air-launched antennas. This study investigated calibration and data processing techniques specific to small-offset ground-coupled antennas that enable additional methods for measuring the electrical properties of the subsurface. These methods require a calibrated GPR system (including the antennas) so that the amplitude and spectral content of waves transmitted into the ground are known. When ground-coupled antennas are used, the response of the antennas depends on the properties of the soil directly beneath the antennas and on the height of the antennas above the ground. Since these properties often vary with location along a survey line, the response of the antennas changes accordingly. Therefore, a method is needed to continually estimate the changing response of the GPR system during a survey. We first discuss methods for calibrating GPRs with ground-coupled antennas. Next, a method is presented to estimate the soil properties beneath ground-coupled antennas, and subsequently estimate the GPR system response. Last, we present a procedure for estimating subsurface reflector properties using deconvolution. This can be done when the system response is known, when the shape of the reflector is known, and when the reflector is contained in an effectively homogeneous host material (described below). With this deconvolution method, GPR surveys are interpreted in a manner similar to TDR surveys to provide estimates of the electrical properties of subsurface. These methods were evaluated with synthetic data, and examples are given using field data. These techniques are tailored for use with impulse radars, but could be adapted for use with other radar types as well.
In this study, we were concerned with the dielectric permittivity and electrical conductivity of the subsurface, and assumed that magnetic media were not present (though they do exist; Olhoeft, 2000). We also assumed that the subsurface can be separated into an arbitrary number of sections with arbitrary shapes and sizes, and that each section is an effectively homogeneous medium. In such a material, any reflections returned from the medium will have very small (if not undetectable) amplitudes. Scatterers in these media are generally less than about one-third wavelength in diameter (Rossiter, 1977). The effective material properties can be determined from the properties of the constituents using mixing rules (Sihvola, 1999; Wtorek, 2003). Note, however, that when soil properties gradually change, the medium cannot be described by effective properties. The methods described here may give erroneous results when used with graded media.
|
Network Analysis of a GPR System
|
|---|
Heimovaara et al. (2004), Robinson et al. (2003), Friel and Or (1999), and Feng et al. (1999) described the processing steps necessary to estimate the electrical properties of horizontally layered media from TDR measurements using deconvolution and inversion. In these measurements, guided electromagnetic waves travel along the transmission line formed by the TDR rods inserted into the soil. These waves propagate along the transmission line until a change in electromagnetic properties causes some of the energy to be scattered (reflected, refracted, or diffracted) while some continues (is transmitted) through the change. The change commonly occurs at the contact interface between two horizontal layers, but changes that are more complicated can occur (e.g., capillary fringe above the water table, rough surface, or volume scattering). The amount of energy transmitted and reflected depends on the contrast in the material properties across the change, the geometry of the change, and the polarization of the incident waves. The received waveform can have a complicated shape because the reflections and transmissions that occur at multiple interfaces or multiple propagation paths are separated in time according to their propagation paths, and they recombine with constructive and destructive interference. Waveforms also change shape with propagation through lossy materials that act like frequency-dependent filters. The response of this transmission line system can be described by a time-invariant linear network. In such a network (McGillem and Cooper, 1984), waves enter and exit the network on ports (ports are a portion of an arbitrary but defined bounding surface of the network), and these waves are related by relatively simple equations. For example, the response of a vertical transmission line can be modeled as a two-port network, with a port on the top and another on the bottom of the line. A wave exiting a port due to a wave incident on a port can be calculated using
 | [1] |
where V+ is the waveform that is incident on a port, V– is the waveform leaving a port due to the incident waveform V+, S is the impulse response of the network between the two ports, t is time, and 
is the convolution operator. Equation [1] can be used to calculate the waveform exiting any port due to a waveform incident on any port, and the complete characterization of a network involves finding the impulse response for all port combinations. For a given pair of network ports, if V+ is known and V– is measured, then the impulse response can be found by deconvolution.
A GPR system can be divided into subsystems where the response of each subsystem is described by an N-port network (N = 2 for a two-port network). An N-port network is a linear time-invariant system, where the waves exiting from each port of the network can be described in terms of the waves incident on each port according to
 | [2] |
where V+ is a vector containing the complex amplitudes of time-harmonic waves incident on each port of the network, V– is a vector containing the complex amplitudes of the waves leaving each port of the network, and S is the complex scattering matrix that describes the network response. The amplitudes V+ and V–, and the scattering matrix S are all functions of frequency (
). Note that Eq. [1] was written in the time domain, and Eq. [2] is the multidimensional equivalent (for N ports) to Eq. [1] in the frequency domain. Scattering matrices contain everything necessary to describe propagating waves within the GPR system and in the ground (as long as they are linear systems). Consult Weber (2001), Collin (1992), Kerns (1981), and Simpson (1981) for discussions on the theory of N-port scattering matrices.
An impulse GPR system using ground-coupled antennas can be represented as a series of N-port networks, as shown in Fig. 1
. Although the division between individual networks is somewhat arbitrary, the ports of each network must be defined so that together they form a closed surface containing the network. Ports do not need to be defined on the bounding surfaces of a network that are electrically shielded or are at an infinite distance from the other ports, because no waves can pass through these bounding surfaces. Multiple networks can be combined and represented by a single network. For example, a network representing a clay lens and another representing the surrounding soil can be described by a single network for simplicity. The waves entering and leaving the two-port networks in Fig. 1 are fully trapped, guided, single-mode waves in transmission lines or within signal-processing circuitry. The ports that couple the antenna systems to the subsurface pass plane waves with specific propagation directions and polarizations. In this study, the shallow soil under ground-coupled antennas (to a depth of about a half wavelength) is included in the antenna N-port networks because the properties of this soil affect the antenna response. Equation [3] specifies the output V– of the GPR system (i.e., waves incident on the mass storage unit and on the impulse generator) in terms of the input V+ (i.e., waves produced by the impulse generator and the mass storage unit):
 | [3] |
The quantities V+ and V– are the complex amplitudes of the input and output waves at the given frequency, 
is the cascade operator for scattering matrices (a linear matrix operation, see Simpson, 1981), and the matrix in parentheses is the total response function for the radar system placed in a particular geometric orientation and position over a particular subsurface. The antenna scattering matrices (Sr and St) are specified for specific antenna locations, attitudes, polarizations, height above the ground, and the specific near-field soil properties. The subsurface scattering matrix (Sss) describes a specific subsurface containing specific scatterers at specific attitudes and locations. The N-port network described by Sss is bounded on the bottom by a hemisphere of infinite radius, and on top by a horizontal plane about a half wavelength deep. Calibrating an impulse radar system (i.e., determining the system response) amounts to determining each scattering matrix listed in Eq. [3] (with the exception of Sss) for each frequency in a discretized frequency spectrum. The scattering matrix approach utilizes a steady-state frequency-domain representation that accounts for all multiples.
View larger version (13K):
[in this window]
[in a new window]
|
FIG. 1. Depiction of a ground penetrating radar system as a series of N-port networks. The Pt,t port and the Pr,r port of the antenna system refer to the transmitting antenna feed port and the receiving antenna feed port, respectively. Note that the antenna networks include the soil that lies within the reactive near-field region of the antennas (i.e., the zone up to about a half wavelength in depth beneath the antennas).
|
|
|
GPR System Calibration
|
|---|
Calibrating a ground-coupled GPR system is not a trivial task, and no commercial GPR manufacturer currently provides calibration information. The calibration of a GPR system involves determining the scattering matrix for each subsystem (assuming that each subsystem can be represented by a time-invariant linear N-port network). Accurately measuring the response of the GPR system elements described above requires carefully conducted experiments made such that the terminating impedance at network ports is not changed by the measurement equipment (as would happen by simply connecting an oscilloscope probe). Oden (2006) discussed a variety of techniques that can be used to calibrate impulse GPR systems. The scattering matrix for the transmission lines and the receiver electronics can be measured using TDR, or using a frequency-domain vector network analyzer. The output of the impulse generator can be measured by sampling the signal while the generator is presented with an appropriate load. Calibrating the response of ground-coupled GPR antennas is more difficult than calibrating other system components since the waves that couple to the antenna ports are traveling through the subsurface or air, and it is not possible to redirect these waves into test equipment for measurement. This difficulty is exacerbated by the fact that the response of ground-coupled antennas changes continuously with varying soil properties under the antennas and varying antenna height as the antennas move along the survey line. Several methods are discussed below to surmount these problems.
Often, an antenna range is used to measure the response of antennas in air. Although there are many types of measurements that can be conducted in an antenna range, a typical experiment involves generating plane waves and measuring the antenna response to the plane waves propagating in different directions and polarizations to determine the antenna pattern. Antenna ranges are often located in a shielded anechoic room, or in a large outdoor area with minimal ambient radio frequency noise and clutter. Traditional antenna range measurements will not provide calibrations for ground-coupled GPR antennas because all measurements are made in air. For small GPR antennas, it is feasible to build an antenna range where a proxy for the soil is used to load the antennas. A large tank can be used to make near-field antenna measurements where a mixture of water, acetic acid (water and acetic acid are miscible), and salt can be used to realize a liquid with a relative permittivity in the range between 6 and 80, and a conductivity in the range between 0 and 50 mS/m (or more). Sand boxes should be avoided because it is very difficult to obtain a uniform mixture of sand, water, and other materials. Bourgeois and Smith (1996) have made antenna range measurements involving proxies for soil and using scale models of the antennas. After measuring the response of a scale model, the response of the larger system is determined through frequency scaling. Alternatively, it may be feasible to use the measured antenna response in air, and mathematically generate the response of the antenna in the presence of a lossy half-space (Meinke and Hansen, 2004; Kerns, 1981).
The experiments described above for calibrating a ground-coupled antenna can be expensive for large antennas designed for lower frequency (<250 MHz) operation. As a practical alternative, numerical models of the antenna and the soil beneath the antenna can be used to estimate the antenna response. Precise measurements of the dimensions and material properties of the components making up the antenna are needed for accurate simulations. The electrical and magnetic properties of the components used in the antenna can be measured in the laboratory using the transmission line method (Kutrubes, 1986; Canan, 1999). Once the geometry and material properties of the antenna parts have been determined, antenna response simulations can be made and verified against experimental results. For example, Oden (2006) made finite-difference time-domain (FDTD) simulations using a publicly available simulation program (Giannopoulos, 2007), and then made experimental measurements of the response of a bistatic antenna pair while it was coupled to conductive water (49 mS/m) and coupled to air. These simulations and experiments were made to calibrate a prototype GPR system (Wright et al., 2005) built by the USGS. This prototype system was designed for surveys over lossy media. This impulse radar had 50-MHz resistively loaded shielded-dipole antennas, and used real-time sampling for improved dynamic range. The antenna model used in the simulations was iteratively improved as more accurate details of the antenna construction were added to the simulated model, and as material property measurements of the antenna components became available. This iterative process continued until the model improvements did not reduce the difference between the simulated response and the measured response. With the simulations from the final antenna model, the root mean square difference between the simulated and experimentally measured waveforms was <10%.
Once an experimentally verified antenna model had been generated, a catalog of early-time received waveforms (see below) and antenna scattering matrices was created by simulating the antenna response using various soil properties and antenna heights above the soil. The antenna configuration used for these simulations was a small-offset (113 cm, which was the smallest possible with the 50-MHz prototype antennas) bistatic array. The catalog contained the simulation results for 60 different combinations of relative permittivity (
r = 4, 9, 16, and 25), conductivity ( = 0, 10, 20, 30, 50 mS/m), and standoff (height of the antennas above the ground, d = 2, 7, and 12 cm). The standoff values bracketed the antenna height settings commonly used with the USGS prototype GPR system. These parameters (permittivity, conductivity, and standoff) were arranged on a three-dimensional grid of rectangular cells. Using this gridded catalog, interpolated antenna scattering matrices (St and Sr) and the early-time waveforms can be calculated for a given set of soil properties and antenna height above the ground.
|
Estimating Shallow Soil Properties
|
|---|
Since the response of ground-coupled antennas is a function of the shallow soil properties in the reactive near-field zone of the antennas (the zone up to about a half wavelength in depth beneath the antennas), the soil properties in this region must be known to determine the antenna response and the overall calibrated response of the system. A number of methods are available for estimating the shallow soil properties using GPR. When using horizontal antenna offsets greater than a wavelength, attributes of the surface-wave (a.k.a. ground wave) arrivals such as arrival time and amplitude can be used to estimate the permittivity and conductivity of the ground (Galagedara et al., 2005a,b, 2003; Grote et al., 2003; Arcone et al., 2003; Fisher et al., 1992). A large offset allows the surface waves to be distinguished from the many other early arriving waves at the receiving antenna, and dry soils typically require larger offsets. Oftentimes, data at several offsets are needed to clearly extract the surface waves. Large offset measurements are undesirable in lossy media because larger offsets reduce the received signal level by increasing the distance waves must travel between antennas and thereby reduce the penetration depth. Amplitude vs. offset is another approach where the amplitudes of reflected waves at multiple large offsets are used to estimate soil properties (Jordan and Baker, 2002; Zeng et al., 2000). Generally, multiple offset methods are unattractive because acquiring these data sets is quite time consuming or requires expensive multichannel equipment. Soil properties determined with large and multiple offset methods may not represent the properties that affect the antenna response because of the large scale of measurement and averaging effects. Soil properties can also be estimated using standing-wave dipoles (Wakita and Yamaguchi, 1996); however, these antennas are not often used with impulse GPR systems because of their "ringy" impulse response. Finally, a number of researchers have discussed estimating shallow soil properties using air-launched antennas (Lambot et al., 2006, 2004a,b; Olhoeft and Smith, 2000; Gentili and Spagnolini, 2000). This approach has the advantage of sampling a smaller surface area than the large- and variable-offset methods. The greatest difficulty when using air-launched antennas or ground-coupled antennas with large offsets is the reduced penetration depth compared to small-offset ground-coupled antennas. A new method to estimate surface properties is discussed below that overcomes these difficulties by using small-offset (maximizing the depth of penetration) ground-coupled antennas.
The inverse model for soil properties (IMSP) algorithm presented below can be used to estimate the shallow soil properties under a ground-coupled small-offset antenna pair. The IMSP algorithm estimates the shallow soil properties using the early-time arrivals—that is, the arrivals recorded before subsurface reflections arrive. The amplitude and shape of the early-time arrivals change with soil permittivity, conductivity, and antenna standoff. The early-arriving waveforms are a complex function of several energy-coupling mechanisms between the transmitting and receiving antennas (see Fig. 2
). These mechanisms include direct waves, reflected waves, multiples, lateral or evanescent surface waves, inductively coupled energy, and reverberations within the antenna shields. It is assumed that the top layer is thick enough to prevent reflections from its base from interfering with the early-time arrivals (more discussion below). Finite-difference time-domain simulations show that early arrivals have good sensitivity to changes in soil properties, but these relationships are not straightforward. Waveform attributes such as arrival time and peak amplitude do not have simple monotonic relationships with these properties. The IMSP algorithm overcomes these difficulties using a nonlinear inversion algorithm. This algorithm is based on a catalog of FDTD simulations (see above) of early-time waveforms received by the GPR system for various model parameters (r, , and d). A forward operator is built from the catalog that describes the shape of the early-arriving waveforms as a function of the model parameters. This forward operator is inverted and then applied to the received waveform data to estimate these parameters. The algorithm is very fast because the catalog of early-time arrivals has been calculated a priori. An abbreviated description of the algorithm is given below; for more details consult Oden (2006) and Oden et al. (2006).
The IMSP algorithm estimates soil properties by comparing attributes of the recorded waveforms to attributes predicted by the forward operator. The forward operator used in the IMSP routine provides a set of 20 waveform attributes as a function of the model parameters. Having fewer attributes reduces the computation time needed for the singular value decomposition portion of the algorithm (described below). Conversely, having more attributes than model parameters increases the likelihood of finding eigenvectors that improve parameter resolution. A number of different attribute sets were tested to determine their suitability for use in the inversion:- Spectral attributes: the fast Fourier transform spectral amplitudes of the early-time waveform from 10 to 200 MHz at 10-MHz intervals
- Hilbert attributes: the Hilbert envelope (i.e., the modulus of the waveform and its Hilbert transform) of the waveform sampled in the time domain at 1.5-ns intervals
- Peak attributes: the amplitude and time of the first four waveform extrema (eight attributes)
- Temporal attributes: the arrival time of the energy peak after applying a band pass filter with a bandwidth of 10 MHz in 10-MHz intervals spanning 10 to 200 MHz
- Wave attributes: 20 waveform values taken at even intervals between 10 and 40 ns
The spectral and Hilbert sets gave better results than the other attribute sets because the objective functions based on the other attribute sets have several basins of attraction with similar minimum values to the basin containing the optimal solution (see discussion below). Using these attribute sets would require data with very low uncertainties (also discussed below). Typically, a 0- to 30-ns time window (about 1.5 cycles at 50 MHz) is used to extract early-time portions of the recorded waveform before the waveform attributes are determined, although shorter time windows can also be used. This window is implemented using a 125-MHz cosine squared taper filter to remove late-time arrivals after 30 ns.
The forward operator was constructed from the simulations catalog using a trilinear interpolation (Press et al., 1992). For the model parameters (X = [r, , d]T) at each grid cell corner in the catalog, the FDTD simulations produced a waveform from which attributes are extracted. The forward operator A returns a set of waveform attributes Y as a function of the model parameters X anywhere in model space (i.e., the space spanned by the model parameters in the catalog) by interpolating between known attributes at the grid cell corners. The Jacobian matrix needed for the inversion is constructed in similar fashion. The forward operator is mildly preconditioned by defining the model parameters such that the valid span of each parameter is roughly the same order of magnitude. The response of the attributes across model space have a widely varying character, and some typical responses are shown in Fig. 3
.
View larger version (37K):
[in this window]
[in a new window]
|
FIG. 3. Example of the inverse model for soil properties (IMSP) forward operator. Plots show the change in Hilbert envelope amplitudes of the early waveforms at 7.5, 15, 22.5, and 30 ns as a function of standoff and conductivity. The relative permittivity is 9 for all plots.
|
|
To estimate the model parameters (r, , and d), the inversion algorithm uses the forward operator and the Jacobian matrix with the Gauss–Newton method to iteratively move from an initial model to improved estimates of a solution (Zhdanov, 2002; Scales et al., 1990; Gill et al., 1986). The function A is nonlinear, but can be treated in a piece-wise linear fashion by using a local value of the Jacobian J to find the direction in which the residual (i.e., the right side of Eq. [4]) can be reduced. At each iteration, the pseudo-inverse of J (from the singular value decomposition) is used to find a value of X that further reduces the residual. Five different initial models are specified for each model parameter to uniformly span each parameter's allowed range, making 125 different initial models. The algorithm starts with each initial model X0q, and proceeds by iteratively reducing the residual of every initial model (q of 125). For a given initial model, the iterative process stops when the following relationship (the stopping criterion) is satisfied:
 | [4] |
where 
is the uncertainty between the predicted and measured waveform attributes (Scales et al., 1990). The residual for each initial model q is iteratively reduced until Eq. [4] is realized, or the number of iterations exceeds 250. An acceptable solution to the inverse problem Xq is obtained when the stopping criterion is satisfied, and an Xq value is tabulated if this criterion is met. The Xq values form the solution set to the inverse problem. This approach is similar to the multistart global nonlinear optimization technique, except that (i) minimization stops when Eq. [4] is realized rather than when the global minimum is found, and (ii) multiple solutions are sought (see rationale below).
Deterministic prior information is incorporated into the inversion algorithm by limiting the allowable range of model parameters to the range of the catalog. The permittivity range is bounded using values for dry sand (r = 4) and water-saturated sand (r = 25), and the conductivity range is from 0 to 50 mS/m. The presence of clay minerals in soil can increase relative permittivity values to >25; however, the conductivity values in this situation will usually be in excess of 50 mS/m and the GPR method will not produce usable subsurface results (Olhoeft, 1987). It may be possible to use the IMSP method to estimate surface conductivities >50 mS/m, but this has not been tested.
To illustrate the progress made by each iteration of the algorithm, the problem has been made two-dimensional by assuming that the standoff is known. In this case, the algorithm constrains the standoff from changing, and this known value is used for all initial models. This reduces the number of initial models from 125 to 25, which reduces the run time of the inversion. Figure 4
depicts the evolution of the solution set for this problem. To generate this figure, the IMSP routine was applied to an FDTD simulated waveform for soil properties of r = 9, = 20 mS/m, and d = 7 cm. The relative uncertainty ||||/||Y|| (discussed below) was 10%. As the process evolves, each initial model (triangles in Fig. 4) moves toward a local minimum until either the stopping criterion (Eq. [4]) is satisfied (squares in Fig. [4]), or the algorithm fails to reach a minimum after a large number of iterations (as on the right side of the plot). The size of the circles corresponds to the logarithm of the objective function (i.e., the right side of Eq. [4]). The models where the stopping criterion is met are plotted with a square, and these models make up the solution set. The final parameter estimates are the mean and standard deviation of the solution set (assuming an approximately normal distribution). The solution set becomes more narrowly distributed as the uncertainty between the predicted and measured waveform attributes is reduced.
View larger version (32K):
[in this window]
[in a new window]
|
FIG. 4. Inverse model for soil properties (IMSP) inversion history for a known standoff. The relative uncertainty is 10% (RDP is the relative dielectric permittivity). Triangles are initial models, squares are acceptable solutions, and circles represent the logarithm of the residual at each iteration. Initial models on the left side of the figure descend to acceptable solutions, while those on the right descend to a local minimum that does not meet the stopping criterion (Eq. [4]).
|
|
Three components of the uncertainty between the predicted and measured waveform attributes are considered. The first component, 1, is the root mean square difference between the attributes of the simulated waveforms and those determined through calibration experiments. An estimate of 1 was made using simulated and experimental results for the prototype GPR system antennas over water and in air, which resulted in 1/||Y|| = 7%. The second component of uncertainty, 2, is due to the noise in the survey data (not to be confused with noise in the calibration data), which can be estimated from the recorded trace before the first arrivals. The final component of uncertainty, 3, is due to the rather sparse sampling used to construct the interpolated forward operator. To assess 3, several FDTD simulations were made with model parameter values located inside the grid cells rather than at the corners. Although details have been omitted for brevity, 3/||Y|| 
1% except when the permittivity and conductivity parameters are less than those of the absorber used in the antennas (r = 11 and = 10 mS/m). In this case, the error can be as high as 10% (worst case) for the grid interval used in the simulations catalog described above. The uncertainty component 3 can be minimized by running more FDTD simulations and using a denser grid interval, or perhaps with the use of higher order spline functions. Therefore, the 3 uncertainty component can be made insignificant across all model space with a sufficiently dense collection of FDTD simulations. The application of Eq. [4] assumes that each component of is independent, uncorrelated, and normally distributed.
The results from the IMSP algorithm are only valid if suitable environmental conditions are met. The algorithm assumes that no subsurface reflectors exist at depths less than about half of a wavelength in the local region surrounding the antennas that could cause interference with the early-arriving waveform, and that the shallow soil is effectively homogeneous and isotropic in this region. Since the electrical properties of the soil are very sensitive to moisture content, the presence of wetting and drying fronts in the sensitive region should be avoided. The algorithm is best applied soon after a precipitation event, or after a long hiatus in precipitation. Since the catalog of FDTD simulations was created by simulating a planar soil surface, then the average height of the surface asperities should be less than about 1/16th of the predominant wavelength in air. For rough surfaces, catalogs of FDTD simulations can be created for surfaces of varying roughness. Oden (2006) discussed methods for determining when the suitable environmental conditions are not met that are based on the collected data.
A discussion of uniqueness, optimal solutions, and valid solutions to the inverse problem is in order. In an optimization problem, a single unique minimum (or maximum) is desired. For a linear inverse problem with noisy data, there is no single unique solution. All points in the neighborhood of the optimal solution that are contained in the confidence interval are valid solutions to the inverse problem. The bounds of the solution space are established by the confidence interval, which is determined from the noise distribution of the data and the sensitivity of the inverse operator. These ideas can also be applied to nonlinear problems; however, two situations occur that make solving the nonlinear inverse problem more difficult. First, the solution space bounds (i.e., confidence intervals) cannot be determined with a linear sensitivity analysis, and second, the solution space may not be contained in a single basin of attraction. The primary goal for a nonlinear inverse problem is to find a solution space contained in a single optimum basin of attraction. This is only possible if the data uncertainty is low enough to keep the solution space from "spilling out" into other basins. When the solution set contains a few members outside the optimum basin, the confidence intervals will increase appropriately (i.e., when minima on the right side of Fig. 4 become part of the solution set). The following method is used to verify that the IMSP solution space lies in a single basin of attraction. The solution set statistics are examined across a range of data uncertainty, beginning with a very small data uncertainty to ensure that the solution space lies in a single basin. As the uncertainty in the data increases, the presence of solutions from additional basins is indicated by an increase in the population and by discontinuities in the mean and standard deviation (see Fig. 5
). If the actual data uncertainty is smaller than the uncertainty value where the first population increase occurs, then the only possibility for a solution across multiple basins requires that the deepest basins have similar depths. This possibility can be tested using cluster analysis or testing for outliers. For the IMSP algorithm, this has been done visually by generating diagrams similar to Fig. 4 for small data uncertainties. In each case, the solution set resembled the quadratic bounds of a linear least-squares problem contained in a single basin.
View larger version (28K):
[in this window]
[in a new window]
|
FIG. 5. Solution set statistics vs. relative uncertainty of the data using the Hilbert attributes. Note that the population remains constant for relative uncertainties up to 12% (RDP is the relative dielectric permittivity). This indicates that the solution set is contained in a single basin of attraction.
|
|
Figure 4 shows how the inversion evolves when the optimal basin is deep enough to prevent solution set members in nonoptimal basins when the data have significant relative uncertainty (10%). The following experiment indicates that similar results are expected in other parts of model space. The IMSP algorithm was run on each of the 60 simulated waveforms from the simulations catalog using relative uncertainties of 1 and 10%. The statistics for the results of these inversions are shown in Table 1
. The median value of the solution set standard deviations is similar to the confidence intervals shown in Fig. 5 (confidence intervals taken as ±) and suggested by Fig. 4.
View this table:
[in this window]
[in a new window]
|
TABLE 1. Statistics of the solution sets calculated from all of the 60 cataloged waveform simulations. The median standard deviation ± quartile deviation of each parameter estimate are listed.
|
|
The sensitivity of the waveform attributes to the model parameters was investigated by calculating the mean sensitivity for all of the initial 125 parameter values (these values uniformly spanned each parameter range). The mean sensitivities calculated across all initial parameter values are shown in Fig. 6
. Both the spectral and Hilbert attribute sets have similar sensitivities. Note that the conductivity sensitivities are about an order of magnitude less than the standoff sensitivities, which results in a larger confidence interval for conductivity estimates. Note also that the principle components of the solution set covariance are not generally aligned with the parameter axes, and there is usually (the covariance is a function of parameter space) significant covariance between conductivity and standoff. To improve confidence intervals, the standoff could be independently measured (perhaps using an acoustic ranging sensor). Then the algorithm could more accurately estimate permittivity and conductivity, and would require less execution time. Since some of the early-time and high-frequency attributes have relatively small sensitivities, it may be possible to remove them from the algorithm to decrease execution time.
View larger version (24K):
[in this window]
[in a new window]
|
FIG. 6. The mean sensitivities calculated across the 125 initial parameter values for the spectral and Hilbert attribute sets (RDP is the relative dielectric permittivity and Cond. is conductivity).
|
|
Finally, the initial model set must be sufficiently dense to provide a solution set with a large enough population to provide quality statistics and to delineate the solution space. The density of starting models must be enough to supply all basins of attraction with a sufficient number of starting models. With a sufficient solution set population, the solution space bounds (i.e., confidence intervals) can be described statistically. The IMSP algorithm was applied to many simulated data sets spanning the parameter space. In all cases, at least half of the initial models became part of the solution set, which resulted in reliable statistics. Using 125 starting models, the IMSP algorithm requires about 3 s of computing time on a 2-GHz Pentium PC to process a GPR trace. The following modifications may decrease execution time: decreasing the number of starting models, decreasing the number of attributes, or using a global optimization approach followed by a bootstrap method for determining the confidence intervals.
|
The Mud Lake Field Experiment
|
|---|
Data from a survey conducted at the Mud Lake Research Station at the Idaho National Laboratories (43°51'1''N, 112°31'57''W) were used for testing the performance of the IMSP algorithm. The USGS prototype GPR system designed for use in conductive ground was used in the survey (this GPR is described above), with the smallest possible antenna offset (113 cm center to center). The soil at the survey site is composed of silt and clay derived primarily from Tertiary rhyolitic volcanic ash (Ralston and Chapman, 1969). The soil was dry to the touch (the site is in a high desert), and the water table is at least 50 m deep. The soil appeared to be essentially homogeneous (as determined by visual inspection) in terms of grain size, texture, and color for the shallow depths investigated by the IMSP algorithm (depths down to 1 m). Soil samples were collected (the original soil structure was not preserved) from a single boring at depths of approximately 15, 30, 45, and 60 cm. The broadband electrical properties of the soil samples were determined by laboratory tests using the transmission line method (Kutrubes, 1986; Canan, 1999), which measured the real and imaginary permittivity throughout a frequency range of 3 mHz to 3 GHz. The fact that the laboratory measurement values were all very similar supports the homogeneous soil assertion. The laboratory measurements at 50 MHz (the peak frequency of the GPR response) are shown in Fig. 7
, where the measured imaginary permittivity is represented as conductivity. The ground surface at the site had broad hummocks, which caused the antenna standoff to vary as the antenna cart rolled along.
View larger version (55K):
[in this window]
[in a new window]
|
FIG. 7. Left panels show estimates of soil properties from the inverse model for soil properties (IMSP) algorithm from the Mud Lake site. The bars indicate the uncertainty in the parameter estimate (i.e., ± one solution set standard deviation). Triangles indicate laboratory results (RDP is the relative dielectric permittivity). The results are invalid in the shaded regions due to rough surface scattering and equipment problems. Right panels show the actual ground penetrating radar (GPR) data, the waveform attributes, and the predicted attributes.
|
|
Figure 7 shows the results from the IMSP algorithm for the survey line that traversed the soil sample boring. The simulations catalog and the IMSP forward operator do not precisely reflect the actual antenna response, and the uncertainty bars shown in the figure were calculated from the uncertainties associated with the FDTD simulations and the forward operator. The laboratory results for the soil properties (triangles) are shown along with the estimated soil properties and antenna standoff. The estimates of conductivity and permittivity are similar to the laboratory values, and the standoff estimates are also similar to the expected value. The antenna heights were adjusted to 7 cm when the antenna cart was fabricated, but later damage to the steering yoke caused a small (undetermined, but estimated at <2 cm) reduction in height.
Although these test results do not rigorously validate the IMSP algorithm, they are encouraging. Note that there is significant negative correlation between the standoff and conductivity estimates, but the permittivity estimates are only weakly correlated with the other parameter estimates. A possible physical reason for this correlation is that variations in standoff and conductivity change the phase of the waveforms via changing travel time and a frequency-dependent reflection coefficient, respectively. This is not obvious when inspecting the FDTD modeled waveforms because both amplitude and phase change in the waveforms as each parameter changes. From an inversion perspective, the quality of the conductivity estimates is affected by low sensitivity and the fact that the principle components of the solution set covariance are not generally aligned with the parameter axes. An independent standoff measurement would reduce the uncertainty of the conductivity estimate that is due to cross-correlation with the standoff estimate.
Unsuitable environmental conditions appear near the hump in the middle of the line (i.e., a nonspecular surface with an asperity larger than 1/16th of the predominant wavelength in air), and also in the 20- to 22-m interval, where intermittent low-amplitude traces were recorded (faintly visible in the pseudosection) due to equipment problems. In some cases, the IMSP algorithm could not find any acceptable solutions at positions where nonideal conditions exist. In all cases, the results are considered invalid where these unsuitable conditions exist. Indeed, a major limitation to the IMSP algorithm is that there are many field situations where unsuitable environmental conditions exist.
|
The Scattering Response Function
|
|---|
Once the properties of the shallow soil under the antennas have been estimated, and the calibrated GPR system response (i.e., all scattering matrices in Eq. [3] except Sss) has been determined, a total response matrix Stotal can be defined for a specific arrangement of subsurface scatterers (which are described by Sss). This total response matrix can be used to calculate the waveform recorded by the GPR for the specific subsurface scatterers. For example, total response matrices can be defined for scatterers embedded in a homogeneous medium such as horizontal planar reflectors, horizontal cylindrical reflectors, or reflectors with any known arbitrary shape. Closely related to the total response matrix is the target response matrix Star, which is the GPR response to a single perfectly reflecting scatterer embedded in a homogeneous medium. The target response matrix is a function of the size, shape, attitude, and position of the scatterer. Thus, the target response matrix accounts for the diffraction of incident waves on an object, but not the variation in the reflection coefficient due to variations in the electrical properties of the scatterer. The target response matrix can be used to estimate the material properties of a scatterer using the method described below.
For a scatterer embedded in an effectively homogeneous medium, the subsurface scattering matrix Sss is the cascaded combination of two scattering matrices:
 | [5] |
where Sm and Ssc are the scattering matrices for the host medium and the embedded scatterer, respectively. The N-port network representing the host medium containing the scatterer will have a set of ports coupling waves traveling in all directions and with all polarizations to the antennas, and another set of ports coupling waves traveling in all directions and with all polarizations to the scatterer. The target response function Star is the GPR response to a perfect reflector in a homogeneous medium and is given by
 | [6] |
where Sm is the scattering (propagator) matrix for waves traveling through the host medium from the transmitting antenna to the scatterer, and Ssc,p is the backscattering matrix for the perfect reflector.
In general, the terms in the scattering matrix Ssc for a reflector are a function of the polarization of the incident waves, their angle of incidence, the shape of the reflecting object, and the contrast between the material properties of the scatterer and its surrounding medium. For many commonly shaped objects, the terms of Ssc can be determined from the transmission and reflection coefficients, which are treated in many texts on electromagnetic theory (e.g., Balanis, 1989). Consider the simple case where the embedded scatterer is a horizontal scattering plane in the Fraunhofer region of the antennas (many wavelengths from the antennas) so that only the reflection from the normally incident wave will produce a significant response at the receiving antenna. In this case, it is pointless to evaluate terms of the host medium's scattering matrix Sm for propagation directions other than vertical. Since we are only concerned with the reflection of the downward traveling normally incident wave with the same polarization as the transmitting and receiving antennas (transmission is ignored), the only nonzero term in Ssc will be that of the reflection coefficient at normal incidence. This coefficient is given by
 | [7] |
where 1 and 2 are the permittivity of the host medium and the scatterer, respectively; 1 and 2 are the conductivity of host medium and the scatterer, respectively; i =
(–1); and is radian frequency. Note that the permittivity and reflection coefficient 
are functions of frequency. Thus for a flat planar reflector, the total response for copolarized antennas to a horizontal plane can be defined as
 | [8] |
where the nonzero matrix entry is the geometrical optics response (i.e., Fraunhofer region; Smith, 1997; Balanis, 1989) of the waves that have traveled vertically down through the antenna ports at depth z1, down through the host medium, reflected off of the horizontal plane at depth z2, and traveled back up through the host medium. The target response function Star for a (perfectly reflecting) horizontal plane is
 | [9] |
This is the response of the GPR system to a perfect horizontal reflector at depth z2. If the response to a horizontal plane V– is measured in a survey, and the target response function Star is known, then the reflection coefficient can be determined from
 | [10] |
The target response function in Eq. [10] can be determined for a cylindrical reflector such as a pipe, a spherical object, or an arbitrarily shaped scatterer. This would allow the reflection coefficient to be determined for any scatterer. This method assumes that no transmission occurs through the scatterer that could lead to later reflections. It is also assumed that the reflection coefficient is the same for all angles of incidence, which of course is not the case. The reflection coefficient is approximately the same for waves near normal incidence, however, which in most cases provides the biggest contribution to reflected energy for small-offset antenna arrays. Since the terms of Eq. [10] are functions of frequency, this frequency-domain equation is equivalent to a time-domain convolution operation, and determining the reflection coefficient amounts to a deconvolution operation. After determining the reflection coefficient at each frequency, the material properties can be estimated by inverting Eq. [7]. Additional constraints may be placed on this inversion by using a permittivity distribution model such as the Cole–Cole model (Wtorek, 2003; Sihvola, 1999) to represent the frequency-dependent permittivity.
|
Deconvolving for Reflector Properties
|
|---|
The equations presented above (Eq. [5–10]) can be used to estimate the electrical properties of subsurface objects with known shapes. The deconvolution method discussed above is similar to the reflection coefficient method used in TDR surveys, except that it can be tailored for objects with a known shape. This process is demonstrated below using a numerical simulation of a 10-m-wide by 2-m-thick clay lens embedded in a homogeneous medium. The simulation was made using a dispersive 2.5-dimensional forward modeling code called 2D Radar (Powers, 1995). The simulated antennas were co-located Hertzian dipoles on the surface, and the transmitted wavelet has most of its energy in the band from about 75 to 325 MHz. The host medium was an ideal moist sand with a relative permittivity of 10 and zero conductivity. The clay lens had a permittivity defined by the Cole–Cole distribution (r,dc = 100, r,
= 25, 
= 3.0 ns, 
= 0.8, and = 100 mS/m; see Wtorek, 2003 for more information about the Cole–Cole distribution) as shown in the lower left panel of Fig. 8
. The simulated waveforms for a perfectly conducting plane at a depth of 2 m and for the clay lens at a depth of 8 m are shown in Fig. 9
. The signal/noise ratio for the clay reflection is 10:1. The upper waveform in Fig. 9 is the simulated time-domain GPR target response (–StarV+) for an infinite-conductivity plate (an infinite-conductivity reflector has a reflection coefficient of –1) at a depth of 2 m. In practice, the target response can be calculated using Eq. [3] (assuming that the antennas have a symmetrical radiation pattern), where the lower reference surface for the antenna networks is a plane greater than half a wavelength deep, and Sss is the identity matrix. Assuming far-field propagation, Eq. [9] can be used to calculate the response to a perfectly conducting plane at another depth, which can then be used to estimate reflector properties as described below.
View larger version (36K):
[in this window]
[in a new window]
|
FIG. 8. The upper-left panel shows the reflection coefficient at the surface of the clay lens, and the lower-left panel shows the permittivity of the clay lens (conductivity is included in the imaginary component). The upper-right panel shows the reflection coefficient calculated from low noise survey data (thick lines) and data with a signal/noise ratio of 10:1 (thin lines). The lower-left panel shows the corresponding permittivity determined from the simulated data.
|
|
View larger version (14K):
[in this window]
[in a new window]
|
FIG. 9. Simulated waveforms for a 2-m-deep perfectly conducting plane (top) and for an 8-m-deep clay lens (bottom). The signal/noise ratio for the lower waveform is 10:1.
|
|
The process of deconvolving for the properties of a planar reflector is outlined in Fig. 10
. The first step is to extract the wavelet corresponding to the reflection from the planar interface and convert it to a frequency-domain representation. Next, the target response Star corresponding to the shallow soil properties under the antennas and for a perfect horizontal planar reflector at depth z is calculated from Eq. [9]. The depth z of the reflector is systematically varied until the correlation between the system response function and the recorded reflection is minimum phase. After determining the proper depth of the reflector, Eq. [10] is used to deconvolve for the reflection coefficient. Finally, the permittivity is determined by inverting Eq. [7]. The results for the simulated clay lens are shown in Fig. 8. Note that the reflection coefficient determined by deconvolution is not identical to the actual reflection coefficient (noise aside). This may be due to Fresnel diffraction region effects, or due to difficulties in resolving both the true depth to the reflector and the phase of the reflection coefficient. For the noisy reflection, only a limited bandwidth can be used in estimating material properties. This limitation can be partially mitigated by stacking traces to reduce noise; however, the limited bandwidth of impulse radars (typically about two octaves) will always be problematic. These results indicate that it is possible to estimate the material properties of subsurface reflectors by deconvolution when using ground-coupled antennas—if a calibrated system is available.
View larger version (23K):
[in this window]
[in a new window]
|
FIG. 10. Processing outline to estimate the electrical properties of a planar reflector (TX is transmission, RX is receiving, FDTD is finite difference time domain, and Star(z, ) is the target response function).
|
|
An antenna calibration survey for the 50-MHz prototype GPR (introduced above) provided a data set that was used to examine the feasibility of the deconvolution method (in addition to verifying FDTD simulations of early-time arrivals). For this survey, the antennas were placed on floats and deployed on the surface of a lake. At the survey location, the lake bottom is comprised of sand and clay, and is nearly horizontal at a depth of about 4 m. The scattering matrices for each network in Eq. [3] were determined for this radar system a priori, with the exception of Sss. Using the process outlined above, the mean reflection coefficient determined for the lake bottom from the reflected wavelet was 0.21 ± 0.10 above the frequency range of 50 to 200 MHz; and the estimated lake bottom relative dielectric permittivity (RDP) was 34 (+18,–12). The large uncertainties are due to the large amount of surface clutter (signal/clutter ratios near 2:1) in the data, and the small amplitude reflection resulting from conductive water (49 mS/m). Clutter removal was not possible since no position data were taken during the survey. The low signal/clutter ratio and the limited bandwidth precluded making conductivity estimates with <100% uncertainty. Three lake bottom sediment samples were taken, and the calculated mean reflection coefficient from the laboratory analyses was 0.32 ± 0.03 at 100 MHz and the RDP was 18 ± 0.06. The difference between the mean material properties determined from GPR data and from laboratory test results is approximately one standard deviation. Although this example is by no means a systematic verification of the method, the results indicate that the accuracy of the system calibration is reasonable, and that it may be feasible to estimate the material properties of subsurface reflectors by deconvolution when using ground-coupled antennas. Ideally, reflection data with signal/noise ratios better than 10:1 would be used. Unfortunately, no further surveys were made with the prototype GPR, because the system was dismantled to obtain parts for other projects when funding for the prototype GPR project was exhausted.
|
Conclusions
|
|---|
The methods presented here require a calibrated GPR system. Calibrating a ground-coupled GPR system involves several techniques, and is not a trivial task because the antenna response depends on the shallow soil properties. A prototype GPR system was calibrated using laboratory measurements to calibrate the electronic parts of the system, and a combination of numerical simulations and physical experiments to calibrate the antennas. The resulting calibrations enabled fast and efficient methods for estimating soil properties such as the IMSP and deconvolution methods. These methods can be used with small-offset ground-coupled antennas, which have the advantage of increased penetration in comparison to air-launched antennas or large-offset antenna pairs. Evaluation of the IMSP algorithm using synthetic data shows that the algorithm can provide a solution from a single basin of attraction if the system calibration is sufficiently accurate. The average uncertainty in the IMSP soil property estimates at Mud Lake were ±1.9 for the relative permittivity, ±4.4 mS/m for the conductivity, and ±0.62 cm for the antenna height. There is significant cross-correlation between the standoff and conductivity estimates, and an independent measure of antenna height would reduce the conductivity uncertainty. It is important to note that the IMSP forward model depends on the GPR antenna system, and that different results should be expected for different antennas. The deconvolution example from the lake bottom, with a 2:1 signal/clutter ratio, is not accurate enough to provide more information than can be deduced by examining a pseudosection. The synthetic example with a 10:1 signal/noise ratio indicates, however, that it is possible to calculate useful material property estimates from reflected waveforms. This would be useful for planar reflectors that do not produce diffraction hyperbolae. It is unfortunate that more field surveys were not made with the USGS prototype system so that these methods could be further tested. This system was very bulky, making it difficult and expensive to deploy, and subsequently, only a small number of surveys were conducted where ground truth was obtained before the system was dismantled. Nonetheless, the preliminary results presented here are encouraging, and show that it is possible to obtain effective calibrations for ground-coupled antennas that enhance the capability of GPR. Further work is certainly needed to advance these methods to general field use.
The accuracy and uncertainty of the estimates returned by the methods presented here are affected by how well the system has been calibrated, how much the survey conditions deviate from the assumed ideal conditions, the amount of noise and clutter, the presence of anisotropic soils, and the presence of small-scale heterogeneities in the soils. The accuracy of the methods can be improved with more accurately calibrated antennas. Future work should address this issue and strive for more accurate calibration methods. Future work should also involve the full calibration of commercial ground-coupled GPR systems so that calibrated measurements of deep vadose zone regions beyond the penetration depth of air-launched or large-offset systems can be made. This, in turn, will provide more information to vadose zone studies by improving our ability to measure moisture content, measure hydraulic properties, delineate fluid flow pathways, evaluate aquitards, and monitor contaminants.
|
REFERENCES
|
|---|
- Arcone, S.A., R.P. Paige, and L. Liu. 2003. Propagation of a ground-penetrating radar (GPR) pulse in a thin-surface waveguide. Geophysics 68:1922–1933.[CrossRef]
- Balanis, C.A. 1989. Advanced engineering electromagnetics. John Wiley & Sons, New York.
- Binley, A., P. Winship, L.J. West, M. Pokar, and R. Middleton. 2002. Seasonal variation of moisture content in unsaturated sandstone inferred from borehole radar and resistivity profiles. J. Hydrol. 267:160–172.[CrossRef]
- Bourgeois, J.M., and G.S. Smith. 1996. A fully three-dimensional simulation of a ground penetrating radar: FDTD theory compared with experiment. IEEE Trans. Geosci. Remote Sens. 34:36–44.[CrossRef]
- Bowling, J.C., A.B. Rodriguez, D.L. Harry, and C. Zheng. 2005. Delineating alluvial aquifer heterogeneity using resistivity and GPR data. Ground Water 43:890–903.[Medline]
- Brewster, M.L., A.P. Annan, J.P. Greenhouse, B.H. Kueper, G.R. Olhoeft, J.D. Redman, and K.A. Sander. 1995. Observed migration of a controlled DNAPL release by geophysical methods. Ground Water 33:977–987.[CrossRef][Web of Science]
- Canan, B. 1999. Dielectric properties of mixtures of clay–water–organic compounds. Ph.D. diss. Colorado School of Mines, Golden.
- Collin, R. 1992. Foundations for microwave engineering. 2nd ed. McGraw-Hill, New York.
- Feng, W., C.P. Lin, R.J. Deschamps, and V.P. Drnevich. 1999. Theoretical model of a multisection time domain reflectometry measurement system. Water Resour. Res. 35:2321–2331.[CrossRef]
- Fisher, E., G.A. McMechan, and A.P. Annan. 1992. Acquisition and processing of wide-aperture ground penetrating radar data. Geophysics 57:495–504.[CrossRef][Web of Science]
- Friel, R., and D. Or. 1999. Frequency analysis of time-domain reflectometry (TDR) with application to dielectric spectroscopy of soil constituents. Geophysics 64:707–718.[CrossRef][Web of Science]
- Galagedara, L.W., G.W. Parkin, and J.D. Redman. 2003. An analysis of the GPR direct ground wave method for soil water content measurement. Hydrol. Proc. 17:3615–3628.[CrossRef]
- Galagedara, L.W., G.W. Parkin, J.D. Redman, P. von Bertoldi, and A.L. Endres. 2005a. Field studies of the GPR ground wave method for estimating soil water content during irrigation and drainage. J. Hydrol. 301:182–197.[CrossRef]
- Galagedara, L.W., J.D. Redman, G.W. Parkin, and A.P. Annan. 2005b. Numerical modeling of GPR to determine the direct ground wave sample depth. Vadose Zone J. 4:1096–1106.[Abstract/Free Full Text]
- Gentili, G.G., and U. Spagnolini. 2000. Electromagnetic inversion in monostatic ground penetrating radar: TEM horn calibration and application. IEEE Trans. Geosci. Remote Sens. 38:1936–1946.[CrossRef]
- Giannopoulos, A. 2007. GprMax: A ground penetrating radar simulation tool. Available at www.gprmax.org/Home/ (verified 30 July 2007). A. Giannopoulos, School of Eng. and Electronics, Univ. of Edinburgh, Edinburgh.
- Gill, P.E., W. Murry, and M.H. Wright. 1986. Practical optimization. Elsevier, Amsterdam.
- Gish, T.J., W.P. Dulaney, K.J.S. Kung, C.S.T. Daughtry, J.A. Doolittle, and P.T. Miller. 2002. Evaluating use of ground penetrating radar for identifying subsurface flow pathways. Soil Sci. Soc. Am. J. 66:1620–1629.[Abstract/Free Full Text]
- Greenhouse, J.P., M.L. Brewster, G.W. Schneider, D.J. Redman, A.P. Annan, G.R. Olhoeft, J.E. Lucius, K.A. Sander, and A.T. Mazzella. 1993. Geophysics and solvents: The Borden experiment. Leading Edge 12:261–267.[Abstract]
- Grote, K., S. Hubbard, and Y. Rubin. 2002. Monitoring spatial and temporal variations in soil water content using GPR reflection data. Leading Edge 21:482–504.[Free Full Text]
- Grote, K., S. Hubbard, and Y. Rubin. 2003. Field-scale estimation of volumetric water content using GPR groundwave techniques. Water Resour. Res. 39:1321.[CrossRef]
- Heimovaara, T.J., J.A. Huisman, J.A. Vrugt, and W. Bouten. 2004. Obtaining the spatial distribution of water content along a TDR probe using the SCEM-UA Bayesian inverse modeling scheme. Vadose Zone J. 3:1128–1145.[Abstract/Free Full Text]
- Hendrickx, J.M.H., B. Borchers, and J. Woolslayer. 2001. Spatial variability of dielectric properties in field soils. p. 398–408. In A.C. Dubey et al. (ed.) Detection and remediation technologies for mines and minelike targets VI. SPIE Proc. Vol. 4394. SPIE, Bellingham, WA.
- Huisman, J.A., S.S. Hubbard, J.D. Redman, and A.P. Annan. 2003. Measuring soil water content with ground penetrating radar: A review. Vadose Zone J. 2:476–491.[Abstract/Free Full Text]
- Johnson, R.H., and E.P. Poeter. 2005. Iterative use of the Bruggeman–Hanai–Sen mixing model to determine water saturations in sand. Geophysics 70:K33–K38.[CrossRef]
- Jordan, T.E., and G.S. Baker. 2002. A conceptual model for the detection of NAPL using amplitude and phase variation with offset (APVO) analysis of ground penetrating radar data. In Symp. Application of Geophysics to Engineering and Environmental Problems (SAGEEP) Proc., Las Vegas, NV [CD-ROM]. 10–13 Feb. 2002. Environ. Eng. Geophys. Soc., Denver, CO.
- Kerns, D.M. 1981. Plane-wave scattering-matrix theory of antennas and antenna–antenna interactions. Natl. Bur. of Standards Monogr. 162. U.S. Gov. Print. Office, Washington, DC.
- Kowalsky, M.B., S. Finsterle, J. Peterson, S. Hubbard, Y. Rubin, E. Majer, A. Ward, and G. Gee. 2005. Estimation of field-scale soil hydraulic and dielectric parameters through joint inversion of GPR and hydrological data. Water Resour. Res. 41:W11425, doi:10.1029/2005WR004237.[CrossRef]
- Kutrubes, D.L. 1986. Dielectric permittivity measurements of soils saturated with hazardous fluids. M.S. thesis. Colorado School of Mines, Golden.
- Lambot, S., J. Rhebergen, I. van den Bosch, E.C. Slob, and M. Vanclooster. 2004a. Measuring the soil water content profile of a sandy soil with an off-ground monostatic ground penetrating radar. Vadose Zone J. 3:1063–1071.[Abstract/Free Full Text]
- Lambot, S., E.C. Slob, I. van den Bosch, B. Stockbroeckx, B. Scheers, and M. Vanclooster. 2004b. Estimating soil electric properties from monostatic ground-penetrating radar signal inversion in the frequency domain. Water Resour. Res. 40(4):W04205, doi:10.1029/2003WR002095.[CrossRef]
- Lambot, S., L. Weihermuller, J. Huisman, H. Vereecken, and M. Vanclooster. 2006. Analysis of air-launched ground-penetrating radar techniques to measure the soil surface water content. Water Resour. Res. 42(11):W11403, doi:10.1029/2006WR005097.[CrossRef]
- Loeffler, O., and M. Bano. 2004. Ground penetrating radar measurements in a controlled vadose zone: Influent of the water content. Vadose Zone J. 3:1082–1092.[Abstract/Free Full Text]
- McGillem, C.D., and G.R. Cooper. 1984. Continuous and discrete signal and system analysis. HRW Ser. Elect. Comput. Eng. Holt, Rinehart, and Winston, New York.
- Meinke, P., and T.B. Hansen. 2004. Plane-wave characterization of antennas close to a planar interface. IEEE Trans. Geosci. Remote Sens. 42:1222–1232.[CrossRef]
- Oden, C.P. 2006. Calibration and data processing techniques for ground penetrating radar systems with applications in dispersive ground. Ph.D. diss. Colorado School of Mines, Golden.
- Oden, C.P., D.L. Wright, M.H. Powers, G.R. Olhoeft, J.B. Rittgers, T. Irons, and A.J. Meininger. 2006. Estimating ground properties under ground-coupled GPR antennas. In Proc. Int. Conf. on Ground Penetrating Radar (GPR 2006), 11th, Columbus, OH [CD-ROM]. 19–22 June 2006. ElectroScience Lab., Ohio State Univ., Columbus.
- Olhoeft, G.R. 1987. Electrical properties from 10–3 to 10+9 Hz: Physics and chemistry. p. 281–298. In J.R. Banavar et al. (ed.) Proc. Int. Symp. Physics and Chemistry of Porous Media, 2nd, Ridgefield, CT. October 1986. AIP Conf. Proc. 154. Am. Inst. Phys., New York.
- Olhoeft, G.R. 2000. Maximizing the information return from ground penetrating radar. J. Appl. Geophys. 43:175–187.[CrossRef]
- Olhoeft, G.R., and S.S. Smith. 2000. Automatic processing and modeling of GPR data for pavement thickness and properties. p. 188–193. In D.A. Noon et al. (ed.) Proc. Int. Conf. on Ground Penetrating Radar (GPR 2000), 8th, Gold Coast, Australia. Proc. SPIE Vol. 4084. SPIE, Bellingham, WA.
- Powers, M.H. 1995. Dispersive ground penetrating radar modeling in 2D. Ph.D. diss. Dep. of Geophysics, Colorado School of Mines, Golden.
- Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. 1992. Numerical recipes in C: The art of scientific computing. 2nd ed. Cambridge Univ. Press, New York.
- Ralston, D.R., and S.L. Chapman. 1969. Water level changes in the Mud Lake area, Idaho, 1958–1968. Water Inf. Bull. 7. Idaho Dep. of Reclamation, Boise.
- Rea, J., and R.J. Knight. 1998. Geostatistical analysis of ground-penetrating radar data: A means of describing spatial variation in the subsurface. Water Resour. Res. 34:329–339.[CrossRef]
- Robinson, D.A., S.B. Jones, J.M. Wraith, D. Or, and S.P. Friedman. 2003. A review of advances in dielectric and electrical conductivity measurements in soils using time domain reflectometry. Vadose Zone J. 2:444–475.[Abstract/Free Full Text]
- Rossiter, J.R. 1977. Interpretation of radio ionterferometry depth sounding, with emphasis on random scattering from temperate glaciers and the lunar surface. Ph.D. diss. Univ. of Toronto.
- Sander, K.A., and G.R. Olhoeft. 1994. 500-MHz Ground penetrating radar data collected during an intentional spill of tetrachloroethylene at Canadian Forces base Borden in 1991. [CD-ROM] Digital Data Ser. DDS-25. U.S. Geol. Surv., Reston, VA.
- Scales, J.A., P. Docherty, and A. Gersztenkorn. 1990. Regularization of nonlinear inverse problems: Imaging the near surface weathering layer. Inverse Probl. 6:115–131.[CrossRef]
- Schmalholz, J., M. Muller, U. Yaramanci, A. Kemna, and H. Stoffregen. 2004. Small scale determination of volumetric water content distribution in the uppermost soil. p. 489–492. In E. Slob et al. (ed.) Proc. Int. Conf. on Ground Penetrating Radar (GPR 2004), Delft, the Netherlands,10th. 21–24 June 2004. Inst. Electrical and Electronics. Eng., New York.
- Sihvola, A. 1999. Electromagnetic mixing formulas and applications. IEE Electromagnetic Waves Ser. 47. Inst. Electrical Eng., London.
- Simpson, G.R. 1981. A generalized n-port cascade connection. IEEE MTT-S Int. Microwave Symp. Dig. 81:507–509.
- Smith, G.S. 1997. An introduction to classical electromagnetic radiation. Cambridge Univ. Press, New York.
- Sneddon, K.W. 2003. Porosity and multiphase fluid saturation from time-lapse ground penetrating radar. M.S. thesis. Colorado School of Mines, Golden.
- Wakita, Y., and Y. Yamaguchi. 1996. Estimation of the soil permittivity and conductivity by a GPR antenna. p. 123–127. In Proc. Int. Conf. on Ground Penetrating Radar (GPR ‘96), 6th, Sendai, Japan. 30 Sept.–3 Oct. 1996. Dep. of Geosci. and Technol., Tohoku Univ., Sendai, Japan.
- Weber, R.J. 2001. Introduction to microwave circuits: Radio frequency and design applications. IEEE Press Ser. RF Microwave Technol. IEEE Press, New York.
- Wollschläger, U., and K. Roth. 2004. Estimating the three-dimensional hydraulic structure of soils from ground-penetrating radar measurements. In In E. Slob et al. (ed.) Proc. Int. Conf. on Ground Penetrating Radar (GPR 2004), Delft, the Netherlands,10th. 21–24 June 2004. Inst. Electrical and Electronics. Eng., New York.
- Wright, D.L., C.P. Oden, M.H. Powers, C.W. Moulton, S.R. Hutton, J.D. Kibler, G.R. Olhoeft, and W.F. Woodruff. 2005. A ground penetrating radar system for high loss environments. p. 1–10. In Symp. Application of Geophysics to Engineering and Environmental Problems (SAGEEP) Proc., Atlanta, GA. 3–7 Apr. 2005. Environ. Eng. Geophys. Soc., Denver, CO.
- Wtorek, J. 2003. EUDEM2 technology survey: Electrical and magnetic properties of soil. IST-2000-29220. Available at www.eudem.vub.ac.be/files/EM_prop_soil_vp.pdf (verified 30 July 2007). Geneva Int. Ctr. for Humanitarian Demining, Geneva.
- Zeng, X., G.A. McMechan, and T. Xu. 2000. Synthesis of amplitude-versus-offset variations in ground penetrating radar. Geophysics 65:113–125.[CrossRef][Web of Science]
- Zhdanov, M.S. 2002. Geophysical inverse theory and regularization problems. Elsevier, Amsterdam.
This article has been cited by other articles:
|
|
|
|
 
S. Lambot, A. Binley, E. Slob, and S. Hubbard
Ground Penetrating Radar in Hydrogeophysics
Vadose Zone J.,
February 25, 2008;
7(1):
137 - 139.
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
