Published online 25 February 2008
Published in Vadose Zone J 7:160-170 (2008)
DOI: 10.2136/vzj2006.0147
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: GROUND PENETRATING RADAR IN HYDROGEOPHYSICS
Application of Microwave Tomography in Hydrogeophysics: Some Examples
Francesco Soldovieria,*,
Giancarlo Priscoa and
Raffaele Persicob
a Istituto per il Rilevamento Elettromagnetico dell'Ambiente (IREA-CNR), Via Diocleziano 328, I-80124 Napoli, Italy
b Istituto per i Beni Archeologici e Monumentali (IBAM-CNR), Via per Monteroni, Campus Universitario, 73100 Lecce, Italy
* Corresponding author (soldovieri.f{at}irea.cnr.it).
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 30 September 2006.
 |
ABSTRACT
|
|---|
We have developed a new strategy for indirect measurement of the dielectric permittivity of the soil using ground penetrating radar (GPR) data. In particular, the GPR data are provided as field scattered by a pipe buried in the soil; the scattered field data is collected at the air–soil interface and under a multimonostatic configuration (B-scan). The strategy is based on the use of a microwave tomographic approach that achieves more robust results with respect to the simple calculation of the round-trip time of the signal. Preliminary results are presented with synthetic data, achieved by means of a finite-difference time-domain code, for both the simple case of a pipe buried in a homogeneous soil and the more complicated case of multiple pipes in a layered medium.
Abbreviations: EM, electromagnetic • FDTD, finite-difference time-domain • GPR, ground penetrating radar • TDR, time domain reflectometry • TSVD, truncated singular value decomposition • TWT, two-way time
|
INTRODUCTION
|
|---|
The possibilities offered by a GPR system in hydrology and underground stratigraphy have been gaining interest among the scientific community in the last few years (Nakashima et al., 2001; Daniels et al., 1995). In particular, there is a strong applicative interest in the fields of intensive agriculture (Hubbard et al., 2002; Huisman et al., 2003) and monitoring of subsoil pollution.
Ground penetrating radar is a fast and simple tool for the determination of the water content of the soil (Huisman et al., 2003; Lambot et al., 2004a). In particular, it allows the creation of spatial maps of water content at an intermediate scale (and with an intermediate resolution) between satellite systems and time domain reflectometry (TDR) systems (Friel and Or, 1999; Mattei et al., 2006).
A first class of radar techniques uses borehole data, which can provide information about the permittivity of the soil up to a depth of 10 m and more (Binley et al., 2001). Borehole measurements, however, are more difficult than GPR prospecting in reflection mode.
A second class of GPR measurements is based on data collected at the air–soil interface. The GPR data are then processed in two steps to provide information about the water content of the soil. The first step consists of determination of the velocity of the electromagnetic wave in the soil (and therefore the dielectric permittivity of the soil), obtained from the two-way time (TWT) of the signal (Huisman et al., 2003; Loeffler and Bano, 2004; Greaves et al., 1996; Tillard and Dubois, 1995). The second step consists of determination of the water content of the soil starting from the estimated dielectric permittivity (Topp et al.,1980; Herkelrath et al., 1991). In this study, we focused on the first step of the procedure, that is, the determination of the dielectric permittivity of the soil starting from GPR data at the air–soil interface.
Two kinds of configuration are customarily used for this procedure (Huisman et al., 2003). The first configuration uses a single measurement point and the electromagnetic (EM) propagation velocity is derived from knowledge of the depth of a reflecting target buried on purpose. This configuration is mainly based on a one-dimensional model of the EM scattering, and thus problems can arise if there is some spurious buried object in the neighborhood of the known object. Moreover, under this configuration, the resolution available in the determination of the object's depth is quite poor, and this affects the precision of the determination of the TWT. Still, in many cases it is not correct to consider the depth of the reflector as a perfectly known quantity, in particular if the same reflector is used for months or years to monitor the time behavior of soil characteristics. In that case, weathering and vegetation or agricultural activity can modify the upper layers of the soil, thus making the target slightly shallower or slightly deeper, or even slightly oblique with respect to the interface. In such a situation, it seems better advised to assume that the depth of the reflector is unknown, even if its presence is known and we have an approximate idea of its depth.
The second configuration uses several observation points. Within this framework, further choices are offered by common offset, common midpoint, or wide-angle radar reflection (WARR) configurations (Daniels, 2004; Galagedara et al., 2003). The characteristics of the soil, in this case, can be worked out from the comprehensive image of the buried object obtained by several GPR traces. In particular, several procedures have been proposed to determine the EM propagation velocity in the soil from a common offset configuration or B-scan "radargram" (Huisman et al., 2003; Tillard and Dubois, 1995). The B-scan configuration is faster and simpler to implement than common midpoint or WARR, and probably is the most commonly used GPR configuration.
For each of these configurations, the result is often provided by direct interpretation of the collected data. The interpretation of the radar signal is a relatively simple task in the case of a single target buried in a homogeneous half-space, However, it can become substantially more difficult if several targets are present (including clutter and different soil layers).
To mitigate these difficulties, in this study we have developed some preliminary examples where synthetic GPR data were processed by means of a microwave tomographic algorithm (Leone and Soldovieri, 2003; Persico et al., 2005). In particular, we made use of a multimonostatic configuration where the source and observation points coincide, i.e., the configuration used was common offset with an offset of zero.
The tomographic approach provides a simple strategy, similar to semblance analysis (Yilmaz, 1987), to determine the dielectric permittivity of the soil (the velocity of the EM wave) through a focusing procedure. This strategy determines the effective dielectric permittivity of the soil as the value that, when inserted into the microwave tomographic algorithm, provides the most focused image of a given buried pipe.
With respect to other methods, the strategy proposed here has some advantages, which are briefly addressed below. First of all, it is suitable for depths up to a few meters, which is an intermediate scale between the one obtainable with TDR and transfer electromagnetic horns (Gentili and Spagnolini, 2000; Serbin and Or, 2004; Lambot et al., 2004b) and that achievable with borehole prospecting (Binley et al., 2001; Rucker and Ferré, 2005). In particular, at the depth of one or a few meters it may be an oversimplification to adopt the one-dimensional geometry customary in TDR because the geometrical spreading of the radiated field is not negligible. For the same reason, the "spatial diversity" of the data can play a relevant role in the proper reconstruction of the buried object and this, in turn, can allow a more reliable determination of the dielectric permittivity of the soil. Also, it is important to stress the fact that an inverse scattering method can take into account the characteristics of the source and of the propagation phenomenon in a fashion that is theoretically more refined than a model that describes the propagation with rays (Soldovieri et al., 2005). On the scale of a few meters, in fact, we may have near-field effects not fully accounted for in ray tracing or one-dimensional models (Rucker and Ferré, 2005).
Another point worth mentioning is that the method proposed here theoretically doesn't need to assume a priori any dispersion law, and it may even help with determining the variation in the permittivity with the frequency, as will be seen. In other methods, some a priori information about the dispersion law is often used, in most cases by assuming a Deybe or a Q-constant dispersion model (Gentili and Spagnolini, 2000; Lambot et al., 2004b). Theoretically, an inverse scattering approach could also provide an estimation of the electrical conductivity of the soil, but this probably would be more difficult to do reliably with respect to the permittivity determination (Lambot et al., 2004b), and thus we do not currently pursue this goal.
The proposed approach requires data that can be gathered with a "simple" GPR, without the need for more complex setups or of particularly refined human expertise. And, finally, the proposed method can partially mitigate the bad effect of an incorrect estimation of the zero time, which is also a relevant goal (Galagedara et al., 2003).
To use a tomographic approach, however, an inverse scattering problem must be dealt with, which entails some mathematical difficulty due to its inherent nonlinearity and ill-posedness (Colton and Kress, 1992). In particular, the ill-posed nature of the problem requires adoption of a regularization technique to make the solution more robust (even if somehow less accurate) against the effects of noise and the parametric uncertainties affecting the model. In addition, nonlinearity of the data–unknown relationship requires attention to the problem of false solutions that can make the inversion results quite different from the ground truth.
These difficulties can be dealt with by making use of a simplified model of the EM scattering. In this study, in particular, we adopted the linearized model founded on the Born approximation (Cui and Chew, 2002; Meincke, 2001; Leone and Soldovieri, 2003; Persico et al., 2005). Under the Born approximation, a stable solution of the problem can be easily achieved and the problem of false solutions is automatically avoided because the adopted model is linear. Moreover, the computational burden is considerably reduced with respect to that customarily needed within different (nonlinear) models, and this makes it possible to investigate large areas.
Last, but not least, the linear Born model allows analysis of the reconstruction capabilities of the solution algorithm in terms of the retrievable spatial variation of unknown objects and thus, ultimately, in terms of the achievable resolution limits, which is a substantially more difficult task in the framework of a nonlinear model.
|
The Microwave Tomographic Approach
|
|---|
The microwave tomographic approach can be broken down into three subdivisions. First, the electromagnetic scattering inverse problem must be formulated in a two-dimensional geometry and in the frequency domain. Second, the problem must be implemented numerically and then inverted. Finally, steps are needed to move from the total field in the time domain (i.e., the usual datum provided by a commercial GPR) to the scattered field in the frequency domain (i.e., the datum suitable for the inversion algorithm).
Formulation of the Electromagnetic Inverse Scattering Problem
To formulate the inverse scattering problem, we have to define the background scenario. This entails the definition of the geometry of the problem that, for this case, consists of two half-spaces (air and soil) separated by a planar interface at depth z = 0 (see Fig. 1
). The upper half-space is made up of free space (dielectric permittivity and magnetic permeability 
0 and µ0, respectively). The lower half-space is assumed homogeneous. In particular, we assume that it shows a relative dielectric permittivity b, and electrical conductivity 
b. The incident field source is a time-harmonic [time dependence exp(j2
ft), where j is the imaginary unit 
(–1), f is frequency, and t is time] filamentary y-directed electric current (transverse magnetic polarization) of infinite extent and invariant along the y axis (two-dimensional source current). The radiation is multifrequency in the band [fmin, fmax]. A multimonostatic observation configuration is assumed, where the locations of the transmitting and receiving antennas coincide. The targets are invariant along the y axis and their cross-section is assumed to be included in a rectangular investigation domain D = [–a,a] x [zmin,zmin + 2b], where a is the semi-extent along the x axis of the investigation domain, zmin is the minimum depth of the investigation domain, and b is the semi-extent of the investigation domain along the depth. The unknowns of the problem are the relative dielectric permittivity profile r(x,z) and the conductivity profile (x,z) inside D. In less technical terms, the buried objects are modeled as inhomogeneities or anomalies of the dielectric and conductivity properties of the soil.
The scattering phenomenon is then recast in terms of the so-called contrast function, defined as
 | [1] |
where
 | [2] |
and
 | [3] |
which are the equivalent "comprehensive" dielectric permittivity to be reconstructed and the equivalent dielectric permittivity of the soil, respectively. Since the buried object is included in the investigation domain, the contrast is different from zero only within the investigation domain D.
Under the Born approximation, the relationship between the unknown contrast function and the scattered field data is provided by the following integral equation (Cui and Chew, 2002; Meincke, 2001; Leone and Soldovieri, 2003; Persico et al., 2005):
 | [4] |
where Es denotes the scattered electric field, measured along the air–soil interface, at the abscissa xs (datum of the problem), 
= 2f; ks is the wave number in the lower half-space (i.e., in the soil), Ge is the external Green's function, r' is the generic point within the investigation domain, and Einc is the incident electric field in the lower half-space. The scattered field is collected over a rectilinear observation domain at the air–soil interface for xs ranging from –xM to xM, where xM is the semi-extent of the measurement (observation) domain (see Fig. 1).
The scattered field is given as the difference between the total field and the incident field. The total field is the field reflected by the soil when the objects are present (Fig. 2b
). The incident field is the field reflected by the soil when the objects are absent and thus it accounts for the reflection at the air–soil interface (Fig. 2a).
View larger version (10K):
[in this window]
[in a new window]
|
FIG. 2. The meaning of the (a) incident field and (b) the total field as the sum of the incident field and the field scattered by object (pipes) for a half-space geometry (Einc is the incident field; Escat is the field scattered by the buried objects; TX is the transmitting antenna; and RX are the receiving antennas).
|
|
The Born approximation in Eq. [4] does not account for the mutual interactions between buried objects and thus it considers each object as if it were alone in the inhomogeneous space. Accordingly, the scattered field due to different objects is assumed equal to the sum of the scattered fields due to each of them (see Fig. 3
).
View larger version (13K):
[in this window]
[in a new window]
|
FIG. 3. Born approximation assuming that no mutual interaction arises between the objects (Einc is the incident field; Escat is the field scattered by the buried objects; TX is the transmitting antenna; and RX are the receiving antennas).
|
|
The external Green's function, Ge, accounts for the electric field radiated in the upper half-space (air) due to a line source electric current with unit strength and located in the lower half-space (soil). By resorting to a spectral relationship (Persico et al., 2005), the external Green's function is written as
 | [5] |
where
 | [6] |
and
 | [7] |
where k0 is the wave number in the upper half-space (i.e., in the air) and u is the generic point in the spectral domain (the conjugate variable with respect to xs). The imaginary part of the square roots in Eq. [6] and [7] is to be chosen non-positive.
The incident electric field in the lower half-space (i.e., the electric field in the soil when the scattering objects are absent), Einc, is given by (Persico et al., 2005)
 | [8] |
where I0 is the strength of the source current and v is the generic point in the spectral domain. The variables w0 and ws are functions of the spectral variable v in Eq. [8], while in Eq. [5–7] they are functions of the variable u.
The problem of reconstructing the targets inside D amounts to solving the integral relationship Eq. [4], where the unknown is the contrast function and the datum is the scattered field. After substituting Eq. [5] and [8] into Eq. [4], it can be rewritten as
 | [9] |
Note that Eq. [9] takes into account the fact that the contrast is a null function outside the investigation domain D, and therefore we have extended all the integration limits for the spatial variables x' and z' from –
to .
Actually, the contrast function as defined in Eq. [1] changes with the frequency. We will neglect this dependence in the inversion algorithm, but the simulated data will account for it. Finally, let us state explicitly that the dielectric and conductive characteristics of the soil influence Eq. [9] because they are contained (through the wavenumber ks of the soil) in the expressions of both the incident field (Eq. [8]) and Green's function (Eq. [5]).
Numerical Implementation and Solution of the Inverse Problem
A discretized version of the integral relationship in Eq. [4] (or equivalently Eq. [9]) is achieved via a method of moments (Leone and Soldovieri, 2003); the contrast function is expanded by means of a Fourier series along the horizontal variable x and by means of rectangular pulses of equal length 
z along the depth z. Moreover, the actual data used for the inversion are discrete and are points (equally spaced) in the frequency and observation domains.
Since, as said, the problem is ill posed, the inversion of the discretized version of the integral relationship Eq. [9], that is the solution of an algebraic system of equations, can be unstable (Colton and Kress, 1992). This requires a regularized inversion procedure. In this study, we regularized the solution by making use of the truncated singular value decomposition (TSVD) expansion of it:
 | [10] |
where the brackets denote the scalar product in the functional space of the data, and un, vn, and n are the left singular functions, the right singular functions, and the singular values of the linear integral relationship Eq. [9], respectively (Bertero and Boccacci, 1998). This method provides a stable solution of the problem. The choice of the index N in Eq. [10] is a tradeoff between the accuracy of the solution (which potentially increases with N) and its stability with respect to noise and uncertainties (which potentially decreases with N). In this study, we have chosen to retain in the expansion of Eq. [10] the terms corresponding to the singular values >0.1 times the first singular value (i.e., we have applied a threshold of –20 dB to the singular values). More details about this inversion algorithm are provided in Leone and Soldovieri (2003).
Passing from Time Domain Total Field to Frequency Domain Scattered Field
The adopted processing poses two relevant problems when time domain GPR systems are used to collect the data. The first one is that a GPR system provides data proportional to the total field present in the observation points, whereas the inversion algorithm processes data proportional the scattered field in the observation points. The second problem is that a GPR system provides data in the time domain, whereas the adopted inversion algorithm is devised in the frequency domain.
Therefore, before performing the inversion, preprocessing operations are needed to pass from raw data (total-field data in the time domain) to scattered-field data in the frequency domain. The first step of the preprocessing is to make the first part of all the radar traces equal to zero, which corresponds to erasing the direct and surface wave contributions at the interface. Alternatively, this first step can be replaced by background removal. Of course, the muting of the first part erases the shallower objects, whereas background removal creates in general some spurious deformation in the image. Whichever the choice, this step essentially provides an estimation of the scattered field from the data. From a theoretical point of view, this method of retrieving the scattered field is not rigorous, especially if the soil is not homogeneous. In many practical cases, however, it is virtually impossible to follow a rigorous procedure, i.e., to calculate the scattered field as the difference between the field measured in the presence of the buried objects and the field measured in their absence, unless one works in a well-controlled situation (Persico and Soldovieri, 2006).
The second preprocessing step is a Fourier transform of the time-domain gated data, so as to achieve scattered-field data in the frequency domain. The result is the data set that will subsequently be "inverted." The synthetic, raw time-domain data used in this study were simulated by means of the GprMax 2D code (Giannopoulos, 2003) based on the finite-difference, time-domain (FDTD) method. Therefore, this preprocessing procedure was followed. In particular, it is relevant to state that the choice of the frequency band used in the inversion algorithm was driven by the spectrum of the probing time signature of the source.
|
Strategy to Retrieve the Permittivity from a Tomographic Approach
|
|---|
We have previously investigated, in one-dimensional geometry, some effects of inaccuracy in the knowledge of the dielectric permittivity of the medium when a microwave tomographic approach is used (Persico and Soldovieri, 2004). In that case, the main effect is a spurious shift in the depth of the targets. Here, we investigate this effect for two-dimensional objects when a multimonostatic, multifrequency measurement configuration is used.
In this section, we will assume some simplifying hypotheses. In particular, we assume that the reference scenario is a homogeneous lossless medium with a relative dielectric permittivity b, and we consider the case of a point-like scattering object located at (0,z0). The scattered field is collected along a line at z = 0 ranging across the interval (–xM,xM).
Under the multimonostatic configuration, the field scattered by a point-like object is given (apart from factors weakly dependent on the frequency) by [H0(2)(βR0)]2 (Harrington, 1961), where H0(2)(x) is the Hankel function of second kind and zero order, β = (0bµ0) is the wavenumber of the medium, and R0 = (x2 + z02) is the distance between the point-like target and the current observation point. Multifrequency data create a variable wavenumber βmin 
β βmax. Under the hypothesis that the scattering target is in the far field zone of the source [βR0 >> 1 
x 
(–xM,xM)] with respect to the observation domain, and by using the large argument approximation of the Hankel function (Harrington, 1961), we can determine that the scattered field is approximately proportional to exp(–2jβR0). In this case, it is known that the reconstruction by means of a tomographic approach under the Born approximation essentially reduces to a backpropagation algorithm, where the scattering point is determined by (Devaney, 1984; Witten and Molyneux, 1988; Meincke, 2001)
 | [11] |
where R = [(x' – x)2 + z'2] is the distance between the observation point and a generic point (x',z') in the investigation domain D.
At this point, we can analyze the effect of an inaccurate knowledge of the dielectric permittivity of the medium by considering in Eq. [11] an estimated value of the wavenumber β1 = β(es/b), where es represents the estimated value of the relative dielectric permittivity of the medium, different from the actual value b. In this case, Eq. [11] becomes
 | [12] |
This equation is derived from optics (Mansuripur, 2002). It has been verified that the inaccuracy in the knowledge of the dielectric permittivity of the medium has two main effects. The first one concerns the delocalization along the depth of the target. This is easily understood by observing that the wavenumber assumed in the focusing procedure [i.e., β1 = β(es/b)] is different from the actual wave number β = (0bµ0). This requires that the velocity of the EM wave assumed in the model is different from the actual one, which translates into a localization of the point-like target at the wrong depth zes=
ztrue, where zes is the estimated depth and ztrue is the actual depth.
In addition, a second effect arises, consisting of some "defocusing" of the object, i.e., a "smearing" of the spot that represents the reconstructed object. This defocusing is more pronounced as the error in the estimated dielectric permittivity increases. The effect of the defocusing can be roughly inferred by considering the evaluation of the integral Eq. [12]. To this end, we can rewrite the exponential term in the integral of Eq. [12] as exp{–2jβR0[1 – (R/R0)(es/b)]}. Therefore, with βR0 >> 1, the exponential term in Eq. [12] is a rapidly oscillating function (and so it gives a negligible contribution to the overall integral) except in the region of the points, where it results in 1 – (R/R0)(es/b) 
0. Accordingly, when (es/b) = 1, the only region that will be imaged will be that about the point (0,z0), because only for (x',z') 
(0,z0), we have 1 – (R/R0)(es/b) = 0 identically on the observation domain (i.e., for every observation point x). It follows that for (es/b) 
1, we can't have the previous condition identically satisfied, i.e., there is no value of the couple (x',z') in the investigation domain D for which 1 – (R/R0)(es/b) = 0 for all the points of the observation domain, and this causes a defocusing effect in the reconstructed image.
Also, it is worth noting that a defocusing effect arises when the frequency behavior and the radiation pattern of the source are not correctly accounted for in the model. Anyway, in the cases presented below, we have not considered these effects.
The considerations above are made clear by the following numerical result, which takes into consideration a medium with an actual relative dielectric permittivity b = 4. A point-like scattering target is located at the position (x,z) = (0,2) m and the scattered field has been measured across the line (–2,2) m in a work frequency band 100 to 1500 MHz. An investigation domain of extent 2 by 4.5 m has been considered.
Figures 4 and 5
depict the modulus of the reconstructed contrast function through Eq. [12] for eight values of the estimated relative dielectric permittivity es, ranging from 1 to 8. In particular, Fig. 4 depicts the reconstructed contrast function for es = 1, 2, 3, and 4, and we can observe the two effects mentioned above. The first one is an incorrect determination of the depth of the reconstructed reflector, which appears deeper than the true position. The incorrect depth location of the object is well described by the "rule" zes=ztrue . Moreover, a defocusing effect is also observed, and it is more pronounced with increasing error in the estimated dielectric permittivity. Figure 5
depicts the reconstructed contrast function for es = 5, 6, 7, and 8. In this case, the reconstruction of the reflector is shallower than its real position. Again, it can be verified that the incorrect depth location obeys to the rule zes=ztrue, and the defocusing effect is again observed. To appreciate the defocusing effect in a quantitative way, in Fig. 6
the maximum value of the modulus of the reconstruction of the point scattering object, Vmax, vs. the estimated dielectric permittivity of the medium is depicted. We can observe than Vmax attains its maximum value vs. es for the exact guess es = b. Finally, to better appreciate the sensitivity of the focusing procedure for values around the true value b = 4, the same analysis has been performed for estimated relative dielectric permittivities of the model, es, ranging from 3.5 to 4.5 with a step of 0.1. Figure 7
depicts the behavior of Vmax vs. the estimated dielectric permittivity of the medium in this case.
View larger version (24K):
[in this window]
[in a new window]
|
FIG. 4. Reconstruction of a point-like target by backpropagation in a homogeneous medium. Actual permittivity = 4. Model permittivities are (a) 1, (b) 2, (c) 3, and (d) 4. The gray scale is normalized to the maximum reconstructed level of (d).
|
|
View larger version (19K):
[in this window]
[in a new window]
|
FIG. 5. Reconstruction of a point-like target by backpropagation in a homogeneous medium. Actual permittivity = 4. Model permittivities are (a) 5, (b) 6, (c) 7, and (d) 8. The gray scale is normalized to the maximum reconstructed level of (d).
|
|
View larger version (16K):
[in this window]
[in a new window]
|
FIG. 6. Normalized graph of the maximum level of contrast vs. the model permittivity for the cases of Fig. 4 and 5. The actual relative dielectric permittivity is 4.
|
|
View larger version (18K):
[in this window]
[in a new window]
|
FIG. 7. Normalized graph of the maximum level of contrast vs. the model permittivity; relative dielectric permittivity ranges from 3.5 to 4.5 with a step of 0.1. The actual relative dielectric permittivity is 4.
|
|
These considerations suggest a possible strategy to determine the correct value of the dielectric permittivity of the soil in absence of a priori information about the location of the reflector. The strategy consists in performing different reconstructions of the buried target while varying the value of the dielectric permittivity of the soil assumed in the Born-based inversion model. The permittivity of the medium (soil) is determined to be the one achieving the most focused image of the target, characterized by the maximum level of the modulus of the reconstructed contrast.
It is also worth noting that, in principle, the best focusing is obtained for es = b independently of the extent of the used work frequency band. In fact, the criterion adopted for determining the permittivity of the soil is based on the achievement of the best focalization and is not related to the intrinsic "quality" of the best focusing.
This explains why, in principle, the proposed method does not require strong a priori information about the frequency behavior of the relative dielectric permittivity of the soil being tested. Actually, the methods could even be applied on suitable sub-bands of the probing signal so as to provide some insight about the frequency behavior of the dielectric permittivity of the soil. It has to be said, however, that as the work frequency band decreases, the size of the focused spot increases, and this makes a proper localization of the object troublesome and the sensitivity of the focusing procedure becomes critical. In other words, the sub-bands cannot be too small.
Finally, it can now be explained why the procedure presented here is not dramatically sensitive to the choice of the zero time of the traces. In fact, a small shift in the hyperbola corresponds, at a first approximation, with a coherent change in both the distances R and R0. This does not affect the fact that the object is best focused only for es = b, even if different choices of the zero time involves different guessed depths of the buried objects.
|
Numerical Results
|
|---|
Here we present more realistic cases. In particular, the simplifying hypothesis of a homogeneous lossless medium, which underlies the results of Fig. 4 to 7
, is removed. In this case we must account for the frequency behavior of the source because, in the simulation of the data, the source-time signature is given by a Ricker wavelet. We have accounted for this fact in a rigorous way by appropriately weighting the subparts of the matrix of the system related to each frequency.
For all the results here, the synthetic, raw time-domain data have been simulated by means of the GprMax 2D code (Giannopoulos, 2003) based on the FDTD method. This makes the inversion code totally independent from the code that creates the data.
A Single Scattering Object Buried in a Homogeneous Soil
In the first example, we consider a multimonostatic measurement configuration with the measurements collected along a line located at the interface and 3 m wide. The measurements are collected with a spatial step of 0.0312 m, and the work frequency band ranges from 450 to 1350 MHz with a frequency step of 25 MHz.
The actual relative dielectric permittivity of the soil is 10 and its conductivity is 0.005 S/m. The scattering object used for the dielectric permittivity estimation is a pipe that is a perfect electrical conductor; its radius is 0.05 m and (with reference to Fig. 8
) its center is located at the point (0,0.5) m.
For this inversion, the investigation domain has been divided into three adjacent subdomains, each of which is 1 m across and extends in depth from 0.1 to 1.6 m. For each investigation subdomain, we have used measurements collected along an observation line 1 m wide. The contrast function has been expanded along a functional basis made up of 31 horizontal Fourier harmonics and 51 vertical steps, and a numerical TSVD truncated at –20 dB has been used in the inversion.
In Fig. 9
, the results obtained for three different estimated relative dielectric permittivity values, es, assumed in the inversion model are shown. The 3-m-wide reconstruction result has been obtained by joining the three tomographic reconstructions for each of the single, 1-m-wide investigation domains. As can be seen, even without the simplifying approximations underlying Eq. [12], the object is most focused when the permittivity assumed in the model matches the actual permittivity. In addition, when the model assumes the true dielectric permittivity b = 10, the reconstructed contrast function attains the maximum level; this can be seen from the gray scale in Fig. 9, which is the same for each of the three images.
View larger version (17K):
[in this window]
[in a new window]
|
FIG. 9. Reconstruction of a single target in a half-space. Actual permittivity = 10. Model permittivity equals (a) 4, (b) 10, and (c) 20. The gray scale is normalized to the maximum reconstructed level of (b).
|
|
To provide a comparison with more classical methods to determine the permittivity of the soil, making use of the same radar traces used to obtain the result in Fig. 9, we have determined the permittivity from the shape of the hyperbolas. In particular, following Daniels (2004), we have performed the following procedure. Preliminarily, we added noise to the raw time domain data (by simulating a white Gaussian noise that gives an average signal/noise ratio of 40 dB on the GPR traces). Then, we made the first part of the traces equal to zero, corresponding to direct coupling through the air–soil interface. Finally, we have filtered the time-domain traces by a band pass filter at 450 to1350 MHz (i.e., the same frequency band used for the tomographic inversion). At this point, the relative permittivity of the soil has been determined from the formula (Daniels, 2004)
 | [13] |
where xm is the mth observation abscissa, tm is the TWT gathered at x = xm, tmin is the minimum return time from the target, gathered at x = xmin, c0 is the propagation velocity of the light in free space, and M is the number of traces considered in Eq. [13]. The minimum point x = xmin is avoided in Eq. [13] because it would provide an indeterminate 0/0 value.
In the framework of this method, the first thing to determine is the choice of the number of points M to be used in Eq. [13]. In fact, as long as M increases, the solution may become more and more unstable; this happens because the scattered signal coming from the target becomes weaker when the trace is gathered farther from the object and thus it is more sensitive to the noise. In particular, in these cases, we have considered M = 31 traces about the central point x = xmin. In Fig. 10
, the radar traces are shown. In particular, Fig. 10a shows the noiseless data, Fig. 10b depicts the noisy data, and finally Fig. 10c shows the noisy data after the filtering procedure. In this case, by using Eq. [13], the value determined for the relative dielectric permittivity of the soil is 8.6 in the case of noiseless data, 73 in the case of noisy data, and about 12 in the case of noisy and filtered data. It can be noted that even with noiseless data, this procedure does not achieve the correct value of the permittivity, which is 10. This arises from the intrinsic uncertainties introduced by the choice of the zero time and by the finite frequency band of the data.
View larger version (63K):
[in this window]
[in a new window]
|
FIG. 10. Radar traces of a single target: (a) noiseless data, (b) noisy data, (c) noisy data filtered between 450 and 1350 MHz.
|
|
In addition, we processed the data of Fig. 10 by an approach based on the generalized Hough transform method (Windsor et al., 2005a,b). This method attempts to estimate four parameters, such as the (x,z) location of the center of the buried pipe, its radius, and the wave propagation velocity in the soil. This is performed by using nonlinear equations relating the four unknown parameters to four points [in the (x,t) plane] falling on the hyperbola diffracted by the pipe. The results of the generalized Hough transform analysis are dependent on the choice of the points on the hyperbola (Windsor et al., 2005a). Thus, to mitigate the effect of the choice of the four points, following a guideline similar to that given in Windsor et al. (2005b), we randomly chose 1000 different quadruples of points and, for every quadruple, we applied the Hough transform method. Finally, we calculated the average and the standard deviation of the simulation results. The results are reported in Tables 1
and 2
for the case of noiseless (Fig. 10a) and noisy (Fig. 10b) data. As can be seen, even in the noiseless case, the estimation of the soil dielectric permittivity is not very accurate, whereas it completely fails when noise is added on the data. The results in the case of the noisy and filtered data (Fig. 10c) are very similar to those in the noiseless case, and therefore are not shown.
Now, let us show the results of the tomographic approach obtained from the same data of the previous example. In particular, we will make use of the noisy traces of Fig. 10b. Figure 11
depicts a close-up of the reconstruction obtained by means of a tomographic approach for two different guessed values for the dielectric permittivity. In particular, Fig. 11a depicts the reconstruction obtained with a model relative permittivity of 10 (the correct value), whereas in Fig. 11b the reconstruction with a wrong model permittivity of 12 is shown (which is the value obtained from Eq. [13]). As can be seen, the reconstruction of the pipe achieved with an exact value of the permittivity is more correctly positioned and, above all, is more focused and exhibits a greater maximum level. A similar result (not shown) was obtained comparing the focusing obtained using a permittivity value of 8.55, obtained from the generalized Hough transform procedure.
View larger version (41K):
[in this window]
[in a new window]
|
FIG. 11. Tomographic reconstruction from the data of Fig. 10b. Actual permittivity = 10. Model permittivity equals (a) 10 and (b) 12. The gray scale is normalized to the maximum reconstructed level of (a).
|
|
Three Scattering Objects Buried in a Two-Layer Soil
In this case, the simulated scenario is represented in Fig. 12
. We inverted the data for this scenario by assuming the simple model above based on the half-space geometry. This example is presented to show the potentialities of the method in less common cases than a single object in a homogeneous medium.
View larger version (10K):
[in this window]
[in a new window]
|
FIG. 12. A stratified scenario with three target in a two layered soil with flat buried interface. The permittivity is 9 in the shallower layer and 20 in the deeper one (p.e.c. = perfectly electric conducting).
|
|
In particular, in Fig. 12, the shallower layer has a relative dielectric permittivity of 9 and a conductivity of 0.005 S/m, whereas the second layer (which starts at the depth of 80 cm) has a relative permittivity of 20 and a conductivity of 0.01 S/m. The parameters of the measurement configuration are the same as the previous example. In the inversion model, we assumed a relative dielectric permittivity of 10 and a conductivity of 0.005 S/m. Figure 13
depicts the radar time-domain traces (after the gating procedure), while the tomographic reconstruction is depicted in Fig. 14
. To better show the reconstruction of the buried interface, Fig. 14 was created by joining the three reconstructed images, each of which was normalized to its maximum. We can see that the three targets are satisfyingly focused. In particular, the target placed in the first layer is the best focused because the relative permittivity assumed in the model (i.e., 10) is very near to the true relative permittivity of the first layer (i.e., 9). The deeper targets, in contrast, are less focused and are imaged deeper than their actual positions. This is consistent with the fact that the effective permittivity of the second layer is larger than the one assumed in the Born model. In fact, the relative dielectric permittivity is 10 in the inversion while the actual relative dielectric permittivity of the second layer is 20. This entails that, for the second layer, the EM propagation velocity assumed in the model is v = 30/(10) cm/ns 9.5 cm/ns, whereas the actual EM velocity in the second layer is v = 30/(20) cm/ns 6.7 cm/ns. Thus the actual wave is slower than we have foreseen, and therefore the inversion algorithm images the objects as deeper than they are.
View larger version (26K):
[in this window]
[in a new window]
|
FIG. 13. Radar traces in the stratified scenario with three target in a two layered soil with flat buried interface.
|
|
View larger version (52K):
[in this window]
[in a new window]
|
FIG. 14. Tomographic reconstruction obtained from the traces of Fig. 13. Model permittivity = 10. The gray scale is normalized to the maximum reconstructed level.
|
|
When we assume that the distance between the buried interface (i.e., the interface between the two layers of soil) and the deeper objects is known, however, we can still determine the permittivity of the second layer by inverting the background permittivity from the formula zes = (b/es)ztrue, where zes and ztrue denote the estimated and true distance in depth of the objects from the buried interface between the two layers. The permittivity assumed in the model used for the inversion procedure is es, while b is the "actual" permittivity of the second layer. Consequently, therefore, in the formula zes = (b/es)ztrue, the quantity to be determined is b. In this case, this method has provided a relative permittivity of the second layer of 21 (compared with an actual value of 20, so that the result seems quite satisfying). The rationale of this procedure is the fact that the buried interface between the two layers of the soil is considered less variable with time than the air–soil interface.
Three Scattering Objects Buried in a Two-Layered Soil with an Irregular Buried Interface
Finally, in Fig. 15
, the case of a two-layered soil with a curved buried interface is shown. The measurement configuration and the inversion parameters of this case are the same as the previous example. In Fig. 16
, the time-domain traces are shown, whereas in Fig. 17
, the tomographic inversion (with model permittivity of 10) is depicted. Figure 17 was created by joining the three reconstructed images, each of which was normalized to its maximum.
View larger version (29K):
[in this window]
[in a new window]
|
FIG. 15. A stratified scenario with three targets in a two layer soil with an irregular buried interface. The permittivity is 9 in the shallower layer and 20 in the deeper one (p.e.c. = perfectly electric conducting).
|
|
View larger version (33K):
[in this window]
[in a new window]
|
FIG. 16. Radar traces in the stratified scenario with three targets in a two-layer soil with an irregular buried interface.
|
|
View larger version (28K):
[in this window]
[in a new window]
|
FIG. 17. Tomographic reconstruction obtained from the traces of Fig. 16. Model permittivity = 10. The gray scale is normalized to the maximum reconstructed level.
|
|
As can be seen from Fig. 16, in this case the raw time-domain traces provided a confused image of the underground scenario. In particular, unlike for the simpler case of a layered soil with flat interfaces (see Fig. 13), this time the buried interface seems indistinguishable from a set of isolated scattering objects. Instead, the tomographic reconstruction provides a more detailed and clearer image of the curved buried interface, especially in the zone where the objects are buried in the second layer. In the left-hand zone, instead, we have some masking effect of the shallower object with respect to the buried interface.
|
Conclusions
|
|---|
We have presented a new strategy to determine the dielectric permittivity of soil by making use of GPR measurements collected at the air–soil interface for a pipe buried on purpose. The strategy is based on the use of a microwave tomographic approach; the unknown dielectric permittivity of the soil is determined to be the one that creates the best "focusing" of the buried pipe.
There are several advantages to the proposed strategy. First, a fully automatic procedure has been developed where we do not need to make a correct choice among the data points (as occurred in the method based on the hyperbolas and also in the method based on the generalized Hough transform [Windsor et al., 2005a]). Second, our strategy mitigates the effect of an incorrect choice of time zero. The approach benefits from an inversion scheme that is able to achieve a stable solution to the problem with respect to uncertainties in the data. Finally, the strategy, along with determining the value of the dielectric permittivity, also creates a clearer image of the investigated domain than that of the raw data. This also allows determination of the time behavior of the investigated buried area.
We have presented preliminary results here, and a more accurate analysis of the strategy will be addressed in future research. In particular, it will be necessary to address the question of how the parameters of the measurement configuration (the extent of the measurement domain and the work frequency band) affect the accuracy and the sensitivity of the results. We also need to address an experimental validation of the approach.
|
REFERENCES
|
|---|
- Bertero, M., and P. Boccacci. 1998. Introduction to inverse problems in imaging. Inst. of Physics Publ., Philadelphia, PA.
- Binley, A., P. Winship, R. Middleton, M. Pokar, and J. West. 2001. High-resolution characterization of vadose zone dynamics using cross-borehole radar. Water Resour. Res. 37:2639–2652.[CrossRef]
- Cui, T.J., and W.C. Chew. 2002. Diffraction tomographic algorithm for the detection of three-dimensional objects buried in a lossy half-space. IEEE Trans. Antennas Propag. 50:42–49.[CrossRef]
- Colton, D., and R. Kress. 1992. Inverse acoustic and electromagnetic scattering theory. Springer-Verlag, Berlin.
- Daniels, D. 2004. Ground penetrating radar. 2nd ed. IEE Press, London.
- Daniels, J.J., R. Roberts, and M. Vendl. 1995. Ground penetrating radar for the detection of liquid contaminants. J. Appl. Geophys. 33:195–207.[CrossRef]
- Devaney, A.J. 1984. Geophysical diffraction tomography. IEEE Trans. Geosci. Remote Sens. 22:3–13.[CrossRef]
- 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 ground-penetrating radar direct ground wave method for soil water content measurement. Hydrol. Processes 17:3615–3628.[CrossRef]
- 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. 2003. GprMax2D, version 1.5. Available at www.gprmax.org (verified 10 Sept. 2007). GprMax Simulator, Edinburgh, UK.
- Greaves, R.J., D.P. Lesmes, J. Lee, and M.N. Toksöz. 1996. Velocity variations and water content estimated from multi-offset, ground-penetrating radar. Geophysics 61:683–695.[CrossRef][Web of Science]
- Harrington, R.F. 1961. Time-harmonic electromagnetic fields. McGraw-Hill, New York.
- Herkelrath, W.N., S.P. Hamburg, and F. Murphy. 1991. Automatic, real time monitoring of soil moisture in a remote field area with time domain reflectometry. Water Resour. Res. 22:857–864.
- Hubbard, S., K. Grote, and Y. Rubin. 2002. Mapping the volumetric soil water content of a California vineyard using high-frequency GPR ground wave data. Leading Edge 21:552–559.[Free Full Text]
- Huisman, J., S. Hubbard, J. Redman, and A. Annan. 2003. Measuring soil water content with ground penetrating radar: A review. Vadose Zone J. 2:476–491.[Abstract/Free Full Text]
- 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, and M. Vanclooster. 2004b. Modeling of ground-penetrating radar for accurate characterization of subsurface electric properties. IEEE Trans. Geosci. Remote Sens. 42:2555–2568.[CrossRef]
- Leone, G., and F. Soldovieri. 2003. Analysis of the distorted Born approximation for subsurface reconstruction: Truncation and uncertainties effects. IEEE Trans. Geosci. Remote Sens. 41:66–74.[CrossRef]
- Loeffler, O., and M. Bano. 2004. Ground penetrating radar measurements in a controlled vadose zone: Influence of the water content. Vadose Zone J. 3:1082–1092.[Abstract/Free Full Text]
- Mansuripur, M. 2002. Classical optics and its applications. Cambridge Univ. Press, Cambridge, UK.
- Mattei, E., A. Di Matteo, A. De Santis, E. Pettinelli, and G. Vannaroni. 2006. Role of dispersive effects in determining probe and electromagnetic parameters by time domain reflectometry. Water Resour. Res. 42:W08408, doi:10.1029/2005WR004728.[CrossRef]
- Meincke, P. 2001. Linear GPR inversion for lossy soil and a planar air–soil interface. IEEE Trans. Geosci. Remote Sens. 39:2713–2721.[CrossRef]
- Nakashima, Y., H. Zhou, and M. Sato. 2001. Estimation of groundwater level by GPR in an area with multiple ambiguous reflections. J. Appl. Geophys. 47:241–249.[CrossRef]
- Persico, R., R. Bernini, and F. Soldovieri. 2005. The role of the measurement configuration in inverse scattering from buried objects under the Born approximation. IEEE Trans. Antennas Propag. 53:1875–1887.[CrossRef]
- Persico, R., and F. Soldovieri. 2004. Effects of uncertainty on background permittivity in one dimensional linear inverse scattering. J. Opt. Soc. Am. A: Opt. Image Sci. 21:2334–2343.[CrossRef][Web of Science]
- Persico, R., and F. Soldovieri. 2006. A microwave tomography approach for a differential configuration in GPR prospecting. IEEE Trans. Antennas Propag. 54:3541–3548.[CrossRef]
- Rucker, D.F., and T.P.A. Ferré. 2005. Automated water content reconstruction of zero-offset borehole ground penetrating radar data using simulated annealing. J. Hydrol. 309:1–16.[CrossRef]
- Serbin, G., and D. Or. 2004. Ground-penetrating radar measurement of soil water content dynamics using a suspended horn antenna. IEEE Trans. Geosci. Remote Sens. 42:1695–1705.[CrossRef]
- Soldovieri, F., R. Persico, and G. Leone. 2005. Effect of source and receiver radiation characteristics in subsurface prospecting within the DBA. Radio Sci. 40:RS3006, doi:10.1029/2004RS003034.[CrossRef]
- Tillard, S., and J.C. Dubois. 1995. Analysis of GPR data: Wave propagation velocity determination. J. Appl. Geophys. 33:77–91.[CrossRef]
- Topp, G., J. Davis, and A.P. Annan. 1980. Electromagnetic determination of soil water content: Measurements in coaxial transmission lines. Water Resour. Res. 16:574–582.[CrossRef]
- Windsor, C., L. Capineri, and P. Falorni. 2005a. The estimation of buried pipe diameters by generalized Hough transform of radar data. p. 345–349. In Progr. in Electromagn. Res. Symp. (PIERS) 2005, Hangzhou, China. 22–26 Aug. 2005. Electromagnetics Acad., Cambridge, MA.
- Windsor, C., L. Capineri, P. Falorni, S. Matucci, and G. Borgioli. 2005b. The estimation of buried pipe diameters using ground penetrating radar. Insight 47:394–399.[CrossRef][Web of Science]
- Witten, A., and J.E. Molyneux. 1988. Geophysical imaging with arbitrary source illumination. IEEE Trans. Geosci. Remote Sens. 26:409–419.[CrossRef]
- Yilmaz, O. 1987. Seismic data processing. SEG Investig. Geophys. 2. Soc. Explor. Geophys., Tulsa, OK.
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
