The Leading Edge; October 2009; v. 28; no. 10;
p. 1234-1239; DOI: 10.1190/1.3249780
© 2009 Society of Exploration Geophysicists
SPECIAL SECTION: Hydrogeophysics |
Identifying dispersive GPR signals and inverting for surface wave-guide properties
Jan van der Kruk
formerly ETH Zurich
H. Vereecken
Forschungszentrum Juelich
Robert W. Jacob
Bucknell University
Corresponding author: j.van.der.kruk{at}fz-juelich.de
The vadose zone is a dynamic environment in which water is retained or transferred into the saturated zone or atmosphere. Knowledge of water content in the near-surface soil layers is important for improving our understanding of groundwater recharge, evaporation, and uptake by crops or natural vegetation (Vereecken et al., 2008). Ground-penetrating radar (GPR) measurements are capable of estimating the subsurface radar velocity, which can be converted to water content using Topp's equation. The direct ground wave is often used to provide radar velocities for the very shallow part of the vadose zone (Huisman 2003). However, this method can only be applied when the subsurface can be approximated by a homogeneous half-space (e.g., a subsurface with a uniform water content profile). Due to the large permittivity of water (
r = 80), the permittivity of soils can change dramatically when precipitation or irrigation is occurring at a site.
At locations where a thin surface layer overlies a substrate medium that has a lower permittivity, or a much larger permittivity/conductivity than the surface layer (Figure 1), pronounced dispersion of GPR waves can be observed and the surface layer acts as a wave guide. For the first case, total reflection occurs beyond the critical angle on the upper and lower interfaces and is called a low-velocity wave guide. For the second case, total reflection only occurs at the upper interface. Although, the lower interface is a strong reflector, some energy is still transmitted across the interface; hence it is a leaky wave guide. In both cases, the electromagnetic waves are trapped within the wave guide and the radar energy is internally reflected, resulting in a series of interfering multiples that manifest themselves as a package of dispersed waves which can propagate over large distances.
These dispersed EM waves contain important information about the nature of the wave guide. However, because the different phases cannot be clearly identified, standard traveltime techniques for estimating the velocity and thickness of the surface wave guide from a common-midpoint (CMP) measurement cannot be applied. While dispersion is a typical feature of seismic data at all scales, so much so that dispersive Rayleigh waves are commonly used to infer subsurface material properties, most GPR users do not always realize that similar dispersion phenomena may be present in GPR data. Recently, several data sets were identified as containing strong dispersion due to low velocity and leaky wave guides (Arcone et al., 2003; van der Kruk, 2006; Strobbia and Cassiani, 2007; van der Kruk et al., 2007).
Here, a three-dimensional finite-difference time domain (FDTD) visualization of the wave propagation within a wave guide is presented showing the trapped electromagnetic waves and the modal behavior. In addition, we show how to correctly identify dispersion present in GPR data, and how to invert for the medium properties. Hopefully, this paper will enable GPR users to correctly identify and invert their dispersive GPR data.
3D modeling of dispersive EM wave propagation in a surface wave guide
To show the underlying wave phenomena, we performed three-dimensional FDTD modeling of electromagnetic waves propagating through a low-velocity surface wave guide with 1 = 23 and h = 0.35 m overlying a half-space with 2 = 11. A dipole current source (Tx) is present at (x,y,z) = (0,0,0) and emits a Gaussian pulse. Figure 2 shows for (a) t = 50 and (b) t = 100 ns the electric field (Ex) distribution in the plane x = 0 (enhanced mpg file available in the online version of TLE). The body wave traveling in air (including the direct air wave) has already traveled outside of the plotting range. The body wave traveling in the ground (lower half-space) is clearly visible (Figure 2a). The head wave (critical refraction from the lower half-space) is also indicated within the wave guide.
Within the surface wave guide, the GPR signal does not consist of one single pulse but of a more elongated wave train due to multiple internal reflections within the wave guide. Figure 3 shows the modal electric field distributions, TE0, TE1, TE2, and TE3 that can propagate through a dielectric wave guide due to the constructive interference caused by the multiple internal reflections. The propagating modes depend on the TE reflection coefficients at either interface and the height and permittivity of the wave guide. The TE0 and TE1 modes can be clearly identified in Figure 2 within the yellow circles. The fundamental TE0 mode has a positive (red) or negative (blue) amplitude over the whole vertical range of the wave guide, whereas the higher-order modes contain alternating positive and negative (or vice versa) amplitudes. Note that the overall shape of the electric field distributions (pattern and amplitude in the y-direction within the wave guide) in Figures 2a and 2b differ because the different phases travel with different phase velocities.
It is obvious that CMP measurements collected at the surface (for z = 0) can show complicated electric field patterns that look completely different from an ordinary CMP data set and therefore require specially adapted processing and inversion tools. Note that similar modal wave propagation occurs for the endfire configuration, where the source and receiver antennas are oriented in the y-direction. In this case, TM modes that travel through the wave guide have different behavior because the dispersion depends on the TM reflection coefficients.
Identifying wave-guide dispersion in GPR data
Figure 4 shows a typical dispersive CMP data set collected with a 100-MHz pulseEKKO 100 system across a terrace of braided river sediments deposited on and adjacent to the Alpine Fault Zone in New Zealand. Nearby trenching revealed a thin layer of organic-rich sandy silt overlying gravel units with varying amounts of finer material. The data, normalized to the maximum of each trace, show that the dispersive waves enclosed in the dotted lines contain most of the energy in the CMP data. This is due to the trapped waves having a geometrical spreading of 
1
r, whereas typical body waves have a geometrical spreading of 1/r. Moreover, the dispersive signal is identified by the shingled reflections present in the data, appearing with a pattern similar to tiles on a roof. The dashed arrows indicate the apparent phase velocities of these shingling events (v = 0.087 m/ns), whereas the solid arrow in Figure 4 is the group velocity and indicates with which velocity the energy is traveling through the wave guide (v = 0.061 m/ns). Note that there is a clear difference between the phase and group velocities, which indicates a frequency-dependent phase velocity. These properties are characteristic for the presence of dispersive GPR signals in CMP soundings.
View larger version (73K):
[in this window]
[in a new window]
|
Figure 4. Measured CMP data set showing clear dispersion characteristics enclosed between the dotted lines. The data are normalized to the maximum of each trace. The solid arrow is the group velocity and indicates with which velocity the energy is traveling through the wave guide (v=0.061 m/ns). Observed phase velocities are indicated by the dashed arrows (v=0.087 m/ns).
|
|
The dispersion observed in Figure 4 is investigated in more detail in Figure 5 which shows the application of different band-pass filters. The yellow dashed lines shown in Figure 5a–e indicate a decrease of phase velocity with increasing frequency v = 
y/t. Figure 5f shows a sudden increase in phase velocity and for increasing frequency again a decreasing phase velocity can be observed. Note that the apparent phase velocity observed in Figure 4 (v = 0.087 m/ns) was dominated by the center frequency which is approximately 60 MHz.
View larger version (94K):
[in this window]
[in a new window]
|
Figure 5. Frequency-band analysis of data set shown in Figure 4 for frequency ranges (a) 24<f<44, (b) 44<f<63, (c) 63<f<83, (d) 83<f<103, (e) 103<f<122, (f) 122<f<142, (g) 142<f<161 and (h) 161<f<181 MHz. The data are normalized to the maximum of each trace. The yellow dashed lines indicate the dominant phase velocity for the specified frequency range. The phase velocity is clearly reducing with increasing frequency in (a)-(e), whereas in (f) a sudden increase of phase velocity is observed, This indicates that a higher-order mode is starting to propagate for these frequencies that again shows a reducing phase velocity with increasing frequency (f)–(h).
|
|
In a process that is similar to a normal moveout stack (NMO), the traveltimes for all traces at each distance y in Figure 5 can be corrected for a linear moveout across a range of phase velocities according to t = y/v and then the traces can be added. A maximum will occur for each panel in Figure 5 due to constructive interference at v = (a) 0.103, (b) 0.089, (c) 0.082, (d) 0.077, (e) 0.075, (f) 0.115, (g) 0.105, and (h) 0.1 m/ns; for other phase velocities, destructive interference occurs within the corrected stack and lower amplitudes will be obtained. When this procedure is carried out in the frequency domain for each frequency (Park et al., 1998), we obtain the phase-velocity spectrum D(v,f)
 | (1) |
where Ê (y,f) are the CMP data E(y,t) shown in Figure 4, transformed to the frequency domain. Figure 6 shows the phase-velocity spectrum D(v,f); red indicates high amplitudes due to constructive interference and blue low amplitudes due to destructive interference. A dispersion curve that shows the phase velocity as a function of frequency is obtained by picking these maxima for each TE mode, which is defined as
 | (2) |
View larger version (73K):
[in this window]
[in a new window]
|
Figure 6. Phase-velocity spectrum for the TE-mode data shown in Figure 4. The thin yellow curves are picked dispersion curves. The yellow dots indicate the phase velocities obtained from analyzing the different frequency ranges in Figure 5. The thin black curves are the dispersion curves for the TE0 inverted model; and the thin blue curves are the dispersion curves for the TE0-TE1 inverted model.
|
|
The TE0 and TE1 dispersion curves are indicated by yellow lines in Figure 6. For low frequencies, the first mode that can propagate above a certain cutoff frequency is the fundamental mode TE0, which is present for frequencies between 20 and 120 MHz. At frequencies above 120 MHz, a sudden increase in phase velocity occurs, similar to that observed in Figure 5e and 5f. For this frequency the first higher-order mode, TE1, is starting to propagate and has a higher amplitude than the fundamental mode. Higher-order TE modes start to propagate at frequencies above a specific cutoff that depends on the permittivity and thickness of the effective surface wave guide and the permittivity of the material below it (van der Kruk, 2007). Note that the phase velocities obtained from Figure 5 are plotted in Figure 6 as yellow dots.
Inversion of the dispersion curves
These dispersion curves form the basis for the processing tools used to interpret dispersive GPR data. The curves picked for the TE0 and TE1 modes (Figure 6) are used to invert for the medium properties of the wave guide and the underlying half-space, using similar techniques as developed in the seismic community, and involve the minimization of the cost function
 | (3) |
where vTEm (fi, 1, 2, h) are the calculated theoretical TEm dispersion curves for a range of models. The example data have the fundamental and first higher-order TE modes, so a combined TE0-TE1 can be performed by minimizing the following cost function
 | (4) |
We employ a combined global and local optimization approach to solve the inverse problem because of the strongly nonlinear nature of the forward problem. For a variety of starting models, a local minimization algorithm based on the simplex search method is initiated. The result is a local minimum for each starting model. The local minimum that produces the closest match between the model and observed dispersion curve is regarded as the global minimum.
Inversion of experimental data
Separate TE0 and TE0-TE1 inversions were carried out on the picked dispersion curves shown in Figure 4 by minimizing the cost functions in Equations 3 and 4. For each suite of inversions, we used 27 different starting models incorporating all possible combinations of 1 = 10, 15, 20; 3 = 3, 6, 9; and h = 0.2, 0.5, and 0.8 m. For each starting model, a local minimum was obtained. The local minimum with the smallest cost function is assumed to be the global minimum. The global minimum obtained for the TE0 and TE0-TE1 inversions are 1 = 23.88, 1 = 7.03, h = 0.32 m and 1 = 23.75, 2= 6.35, and h = 0.33 m, respectively (Tables 1 and 2). Most other local minima were very close to these global maxima. Out of the 27 iterations for the TE0 and TE0-TE1 inversions, 18 and 21 local minima, respectively, were within 10% error from the global minimum, indicating that these inversions were well constrained. The mean value and standard deviations of these local minima are given in Tables 1 and 2. The higher number of local minima close to the global minimum for the TE0-TE1 inversion indicates that including higher-order modes in the inversion results in a better constrained subsurface model (van der Kruk, 2006).
Conclusions
Three-dimensional modeling shows the amplitude distribution of the different modes traveling through a surface wave guide. Due to the constructive and destructive interference of the many multiples that are trapped within the wave guide, the dominant waves traveling through the wave guide have a frequency-dependent phase velocity.
These dispersive GPR signals can be identified using three key characteristics:
- Normalizing the data on the maximum amplitude for each trace shows that most energy is contained within the dispersive waves.
- Shingling reflections are present in the data and indicate different phase and group velocities.
- The phase-velocity spectrum clearly indicates the presence of a frequency-dependent phase velocity.
A frequency-dependent phase velocity will be represented on a phase-velocity spectrum as a general decreasing phase velocity with increasing frequency. A sudden increase in the phase velocity (Figure 6) indicates that higher-order modes are also propagating in the wave guide, visible due to their higher amplitudes in comparison to lower-order modes. Traditional techniques to interpret shallow subsurface properties, such as normal moveout velocity scans or semblance analysis, do not work for dispersive GPR data. Similar processing and inversion techniques, commonly used in seismics (Socco and Strobbia, 2004, Xia et al., 2003) and adapted for the GPR case (van der Kruk et al., 2006) should be used to invert these data. These techniques involve picking dispersion curves from the phase-velocity spectrum and inverting them for the subsurface material properties by using a combined local- and global-minimization algorithm. Additionally, using higher-order modes (van der Kruk, 2006) or performing a combined TE-TM inversion (van der Kruk et al., 2006) results in better constrained inversions.
These wave-guide phenomena should be expected in all situations where horizontal layering occurs and relatively large velocity contrasts are present. For example, low-velocity wave guides can be expected when sandy silt is overlying gravel, saturated soil is overlying dry gravel, a thawed soil layer is overlying a frozen soil layer, or precipitation and/or irrigation resulted in a wet shallow soil layer overlying a dry soil layer. Leaky wave guides can be expected when frozen saturated sand is overlying moist sand, or ice is overlying salt or fresh water.
Suggested reading
"Propagation of a ground-penetrating radar (GPR) pulse in a thin-surface wave guide" by Arcone et al. (GEOPHYSICS, 2003). "Measuring soil water content with ground penetrating radar: A review" by Huisman et al. (Vadose Zone, 2003). "Imaging dispersion curves of surface waves on multichannel record" by Park et al. (SEG 1998 Expanded Abstracts). "Surface-wave method for near-surface characterization: A tutorial" by Socco and Strobbia (Near Surface Geophysics, 2004). "Multilayer ground-penetrating radar guided waves in shallow soil layers for estimating soil water content" by Strobbia and Cassiani (GEOPHYSICS, 2007). "Properties of surface wave guides derived from separate and joint inversion of dispersive TE and TM GPR data" by van der Kruk et al. (GEOPHYSICS, 2006). "Properties of surface wave guides derived from inversion of fundamental and higher mode dispersive GPR data" by van der Kruk (IEEE TGRS, 2006). "Fundamental and higher mode inversion of dispersed GPR waves propagating in an ice layer" by van der Kruk et al. (IEEE TGRS, 2007). "On the value of soil moisture measurements in vadose zone hydrology: A review" by Vereecken et al. (Water Resources Research, 2008). "Inversion of high-frequency surface waves with fundamental and higher modes" by Xia et al. (Journal of Applied Geophysics, 2003).
Acknowledgments:
This research was supported by grants from the Swiss National Science Foundation and ETH Zurich. We thank R. Streich for acquiring the data.
Copyright © 2009 by Society of Exploration Geophysicists
1, 
