GPR based soil electromagnetic parameters determination for subsurface imaging

The problem of estimating the dielectric permittivity and the electric conductivity of the soil starting from GPR measurements is addressed. A new estimation procedure is proposed and checked against synthetic data generated by a FDTD forward solver. A two-dimensional geometry and a two-layered background medium are considered.


Introduction
Subsurface imaging consists in the identification of certain physical properties of objects which are embedded in a layered background medium.As a typical framework, we recall the problem of localizing scatterers which are buried beneath the air-soil interface.
Many physical principles have been exploited to device reconstruction schemes that, starting from measurements collected in a non invasive way, aimed at determining some parameters of interest.In this realm, those procedures based on wave phenomena (e.g., acoustic, electromagnetic waves) are particularly appealing as they allow to set up high resolution imaging procedures from a remote not necessarily in contact with the air/soil interface measurements.
The ground penetrating radar (GPR) is a radar system suitably developed for subsurface imaging which exploits the electromagnetic field to probe the buried spatial region of interest (Daniels, 2004).Basically, it consists of a transmitting Correspondence to: R. Solimene (raffaele.solimene@unina2.it)antenna which radiates an electromagnetic field towards the interface separating the air from the soil.The backscattered field is then collected over a set of positions by means of an array or a moving receiving antenna.
The simplest way to achieve an image of the subsurface scenario consists in displaying raw data as function of time for each measurement position.This way the so-called radargrams are achieved where typically buried objects are reproduced as hyperbolas from which scatterers can be detected and localized.
However, better results in terms of the achievable resolution and, above all, in terms of an increased interpretability of the scene (which becomes very difficult to achieve through radagrams when the scene encompasses different scatterers) are obtained by employing the so-called focussing approaches.
Many focussing methods are widespread in the literature (see for example, Stolt, 1978;Bleinstein et al., 2000;Gilmore et al., 2006;Lopez-Sanchez and Fortuny-Guasch, 2000).All such methods can be considered as a sort of migration where the principle of time-reversing (Fink, 1992) or equivalently of backpropagating (Devaney, 1984) the received field is exploited to collapse the hyperbolas over the scatterers they arose from.Focusing can be also achieved by adopting a tomographic approach where the subsurface imaging problem is addressed as an electromagnetic inverse scattering one (Pierri et al., 2004).Electromagnetic inverse scattering problems are particularly difficult to address because they are non-linear and ill-posed (Colton and Kress, 1992).However, since the primary task in subsurface imaging is to detect, localize and possibly to determine the geometrical features of the buried scatterers, linear approximations for the scattering phenomenon (e.g., Born, Kirchhoff approximations) can be exploited (Pierri et al., 2006).Accordingly, under a linear formulation the problem is cast as the inversion of a linear operator which can be effectively achieved through its Truncated Singular Values Decomposition (TSVD); this way regularized solutions that are stable against uncertainties and noise on data are obtained (Bertero and Boccacci, 1998).
A more detailed discussion and comparison about the various focussing techniques is not addressed herein as it is beyond the plan of the paper.
It is worth noting that even though subsurface imaging shares most of the scientific/technical aspects of free-space imaging, it is definitively more difficult from an applicative point of view due to the presence of the soil which, in general, is not homogeneous and entails dealing with the roughness of the air/soil interface.
Despite of that, many subsurface imaging methods have been developed and proven to work (also against experimental data) under the assumptions of a homogeneous half-space soil layer separated from the air by a planar air/soil interface (Hansen and Johansen, 2000;Piscitelli et al., 2007).However, even under such a simplified scenario, a crucial question still has to be addressed.It is concerned with the a priori knowledge of the soil's electromagnetic parameters (i.e., the dielectric permittivity as well as the electric conductivity of the soil) which are key factors for imaging procedures.
Many techniques have been devised to cope with the problem of estimating the electromagnetic properties of the soil (see Fortuny-Guasch, 2002;Zeng et al., 2000;Soldovieri et al., 2008;Lambot et al., 2004, and references therein).
Here, we present and test a novel and simple technique to estimate the soil parameters from GPR measurements.Such a technique is based on an integral equation which relates the field reflected from the air/soil interface to the Fresnel reflection coefficient.Accordingly, the reflection coefficient is retrieved by inverting such a relationship; then the parameters of the soil are estimated via an algebraic relation from the reflection coefficient.
The plan of the paper is the following.In Sect.2, we describe the adopted scattering configuration.In Sect.3, we introduce the estimation procedure.Section 4 is devoted to showing the validation of the estimation procedure against the synthetic data generated by a FDTD forward solver.Finally, conclusions follow.

Scattering configuration
We consider the scattering configuration depicted in Fig. 1.The background is a two-layered medium.The upper layer is assumed to be the free-space where 0 and µ 0 are its dielectric permittivity and magnetic permeability, respectively.The lower half-space layer is representative of the soil which is assumed as nonmagnetic (i.e., its magnetic permeability is equal to the one of the free-space) homogeneous with dielectric permittivity and electric conductivity denoted as s and σ s , respectively.The two half-spaces are separated by a planar interface (the air/soil interface) at z=0 (see Fig. 1).
In principle, the lower half-space, i.e., the soil, is allowed to be dispersive and containing ohmic losses.Accordingly, in the frequency domain the dielectric permittivity is actually a complex function which depends on the angular frequency ω, that is s (ω)= sr (ω)−j si (ω) and the equivalent permittivity is then given by eq (ω)= sr (ω)−j [ si (ω)+σ s /ω].
A multistatic measurement configuration is considered.This means that for each position of the transmitting antenna the scattered field is collected over a set of different positions occupied by the receiving antenna.We assume that the receiving positions are all characterized by the same distance from the air/soil interface so that the measurements are collected over an aperture synthesized along the x-axis parallel to and at distance d from the air/soil interface (see Fig. 1).It is worth noting that the analysis we are showing implicitly encompasses both the cases of a single source position as well as the case of several source positions.In the first case we refer to the measurement configuration as a single view/multistatic configuration; in the second case we term the measurement configuration as multi-view/multistatic.In any case, both the receiving and the transmitting antennas are assumed to be far enough from the air/soil interface so that the interaction between the soil and the antennas can be considered negligible.This means that the antennas work as if they were within the free-space and their plane wave spectra can be retrieved by a preliminary characterization stage which can be achieved experimentally or by numerical simulations.Herein, however, we consider the canonical case of filamentary electric currents infinitely long along the y-axis.
Scatterers, if any, are buried in the lower half-space and are assumed to be invariant along the y-axis.Accordingly, we deal with the subsurface imaging problem for a twodimensional and scalar geometry.

Estimating the soil parameters
In this section we tackle the problem of estimating the soil's electromagnetic parameters starting from GPR measurements.In particular, we present an estimation procedure which does not require taking a sample of the soil structure in order to proceed to an off-line soil characterization nor identifying a spatial region free from buried scatterers.Moreover, it does not exploit auto-focusing procedures, which generally are time consuming, and does not need plane wave assumption for the incident field.Indeed, it can accommodate a generic incident field with its own plane wave spectrum.
Hence, the method we propose is particularly appealing as it allows an in-situ soil characterization which is achieved by means of the same measurements to be exploited for the (subsequent) imaging problem.
The subtended idea is very simple and consists in establishing a relationship between the scattered field and the Fresnel reflection coefficient at the air/soil interface.By inverting such a relationship one first recovers the reflection coefficient.Then the electromagnetic properties of the soil are inferred from the estimated reflection coefficient.In particular, the reflection coefficient is retrieved for each plane wave contributing to achieve the incident field.Such values are then exploited to estimate the soil parameters (which are always the same regardless of the particular plane wave at hand) through an averaging procedure which also helps in curtailing the effect of the field scattered by the buried objects which can be actually present on the measurements.
As mentioned in the previous section, we consider filamentary electric currents infinitely long along the y-axis as sources of the incident field.Accordingly, the field impinging on the air/soil interface, in the frequency domain, is given by the free-space Green's function of the two-dimensional Helmholtz equation whose spectral representation is In Eq. (1), r≡(x, z) is the generic field point, r s ≡(x s , z s ) is the source position and Moreover, the contribution due to the evanescent waves (i.e., the ones corresponding to |k 0x |>k 0 ) is assumed to be negligible.
When no scatterers are buried in the soil the scattered field is only due to the reflection occurring at the air/soil interface.Accordingly, the scattered field can be expressed as where the incident field (1) has been used, r O ≡(x O , z O ) is the observation point at which the scattered field is collected and is the Fresnel reflection coefficient at the air/soil interface when no scatterers are buried in the lower half-space.If we particularize Eq. (2) to the configuration described in the previous section, we obtain the relationship to be inverted for soil parameters estimation.In detail, Eq. ( 2) can be rewritten as 1 whose kernel is equal to where we recall that d is the quota over the air/soil interface of the transmitting antenna and of the receiving points, whereas the interface is located at z=0.Equation ( 3) is a linear integral operator from which, starting from the scattered field observations (datum of the problem), we can retrieve the reflection coefficient .
Note that, since both E S and depend on the frequency through k 0 , such an equation should be solved separately for each adopted frequency.However, as our final aim is imaging buried scatterers, we are not interested in characterizing the soil parameters as a function of frequency.Rather, we need to know the soil parameters only for a discrete set of frequencies as dictated by the optimum frequency sampling step discussed in (Persico et al., 2006).Of course, such an optimum number of frequencies can be determined when an upper bound on the real part of the dielectric permittivity is known.
Let us assume that the frequency step has been determined and let k 01 , k 02 , • • •k 0h , • • •k 0N k be the corresponding wavenumbers at N k frequencies uniformly taken within the frequency band [k 0 min , k 0 max ].Equation (3) can be rewritten parameterized with respect to k 0 at the h-th frequency as with Hence, the estimation procedure amounts to solving Eq. ( 5) for h for all the N k adopted frequencies.Accordingly, we now focus on the solution procedure.

R. Solimene et al.: GPR based soil em parameters determination
While solving Eq. ( 5) for , two important issues have to be addressed.First, a suitable discretization scheme has to be adopted in order to devise a numerical procedure for solving such equations.Second, the extent of the synthesized measurement aperture has to be chosen so that the filtering introduced by the regularization algorithm, mandatory to obtain a stable reconstruction, is not too severe.
As to the first point, we note, from Eqs. ( 3) and ( 4), that the scattered field is a bandlimited function of bandwidth k 0h with respect to the variable x O −x s .This entails acquiring the scattered field in order to sample the observation variable x O −x s at a step at least equal to π/k 0h .
As to the extent of the measurement aperture, it should be large enough so that the data it captures represent a significant part of the reflected field amplitude dynamic.This extent can be roughly estimated by considering | |=1 in Eq. ( 3) so that a Hankel function is obtained.From the numerical analysis, we have observed that an aperture extent which allows for retaining a Hankel function amplitude of about −13 dB under the maximum is sufficient.
Consistently with previous discussion, both the observation spatial step and the measurement aperture extent are dependent on the working frequency.In order to simplify the procedure and fulfill the criteria mentioned above, in the following examples we will consider fixed the spatial step and the aperture extent.In particular, the spatial step is chosen as π/k 0 max , corresponding to the highest adopted frequency, whereas the aperture is chosen according to the lowest frequency.
At this point, the reflection coefficient at the h-th frequency is recovered by inverting Eq. ( 5) with a TSVD scheme (Bertero and Boccacci, 1998).In detail, said u hn , v nh and σ hn the n-th left singular function, the n-th right singular function and the n-th singular value, respectively, of the operator in Eq. ( 5), the reflection coefficient is estimated as where ˜ h is the reconstructed version of the reflection coefficient and N T is a suitable truncation index which controls the propagation of uncertainties from data to the reconstruction (Bertero and Boccacci, 1998).
The TSVD reconstruction scheme in Eq. ( 6) is then repeated for all the adopted frequencies.
The last step to complete the estimation procedure is to retrieve eq from the estimated reflection coefficient.At the h-th frequency this can be done through the following algebraic relationship where ω h is the h-th adopted angular frequency.
The above relation means that in the ideal case the eq is a function that does not depend on k 0x .However, due to the necessity of adopting the TSVD scheme in order to mitigate the model error and the noise, different values of eq generally can be achieved for the different k 0x .Therefore, as we recover the reflection coefficient for all the k 0x ∈ [−k 0h , k 0h ], we can achieve an estimation of eq for each k 0x and then take an average value.Accordingly, we can estimate eq , at the h-th frequency and for all the N k frequencies by the following average value Most importantly, the averaging procedure has the additional advantage of mitigating the effect of the field scatterers from buried scatterers which entails a model error that can impair the goodness of the soil permittivity estimation given by Eq. ( 6) and founded on Eq. ( 2).Indeed, it is expected that such an averaging procedure will be successful in removing the effect of the field scattered by buried objects.This is because for large scatterers (in terms of the wavelength) the scattered field is expected to be mainly focussed over a limited spectral part.On the contrary, small objects scatter almost isotropically but the level of their scattered field is generally lower than the one resulting from the air/soil interface.
We finally observe that the averaging procedure in Eq. ( 8) is also useful in curtailing the effect of noise still retained in the TSVD reconstructions.

Soil estimation numerical examples
In this section we shown numerical examples to check the effectiveness of the proposed estimation procedure.
Although the proposed estimation procedure can be applied for generic soil parameters, here, we consider only the case of a non-dispersive half-space with ohmic losses.
A frequency band of [0.3, 1] GHz sampled at a step of 50 MHz is employed.The reflected field is taken at a spatial step of 10 cm over a line having a quota (from the air/soil interface) d=0.5 m.The observation variable x O −x s synthesizes an aperture of [−10, 10] m.This is actually achieved by considering a measurement aperture so that -5 m≤x O ≤5 m and two different positions for the transmitting antenna at −5 m and 5 m, respectively.
We start by dealing with situations where no buried scatterers are present.Two reconstruction examples are reported in Figs. 2 and 3 where the real and the imaginary part of eq are reported as functions of the frequency, respectively.
As can be seen, the proposed estimation technique works perfectly in retrieving the soil parameters when non noise corrupts data (see blue lines reported in such figures).Therefore, in order to test the estimation procedure's robustness against noise, we consider a further example where data are corrupted by an additive zero mean white Gaussian noise in order to have a signal to noise ratio (SNR) of 10 dB.The corresponding reconstructions are reported in Figs. 2 and 3 as red lines.For the noisy case the reconstructions are no longer perfect.However, they are still very close to the actual solution for most of the adopted frequencies and the frequency behavior is well mimed.We now turn to tackle the soil's parameters estimation problem when a scatterer is buried in the lower half-space.In particular, we consider a perfect conducting cylinder whose axis is parallel to y-axis and of circular cross section of radius equal to 0.3 m and centered at (1, 0.5) m.The corresponding reconstructions are reported in Figs. 4 and 5.As can be seen, also in this case the estimation of the soil's parameters is very good.More in detail, the proposed procedure returns almost www.adv-geosci.net/19/39/2008/Adv.Geosci., 19,[39][40][41][42][43][44]2008 the true values of sr and σ s for all the frequencies when no noise corrupts the data, whereas when the data are noisy the achieved performance are similar to the ones shown in Figs. 2  and 3. Hence, we conclude that the proposed estimation procedure works well also when scatterers are buried in the soil.

Conclusions
In this paper we have considered the problem of estimating the soil's electromagnetic parameters for subsurface imaging purposes.In particular, we have considered a two-layered background medium where the upper layer was representative of the air whereas the lower one schematized the soil.
An estimation scheme based on a frequency domain integral relationship linking the reflected field and the reflection coefficient at the air/soil interface has been proposed.In particular, such a relationship has been inverted by means of a TSVD procedure.After recovering the reflection coefficient the dielectric permittivity and the electric conductivity of the soil have been retrieved.
We outline that, such an approach can be applied by exploiting the same measurement survey for the subsurface imaging.In fact we have shown that it works also for situations where buried non-cooperative scatterers are present.
Moreover, it can account for the plane wave spectrum of the source and does not require neither a model for the soil dispersive law nor a rough a priori knowledge about the soil parameters.
Synthetic examples, obtained by exploiting the data generated by a FDTD forward solver, have shown the effectiveness of the method.
Edited by: L. V. Eppelbaum Reviewed by: two anonymous referees

Fig. 2 .Fig. 3 .
Fig. 2.Reconstructed real part of the soil equivalent permittivity normalized to 0 reported as function of frequency.Noiseless data have been used for the blue line, whereas the case of noisy data with SNR=10 dB have been considered for achieving the reconstruction reported by the red line.

Fig. 4 .Fig. 5 .
Fig. 4. Reconstructed real part of the soil equivalent permittivity normalized to 0 reported as function of frequency for the case of a buried scatterer.Blue line: noiseless data reconstruction; red line: noisy data reconstruction with SNR=10 dB.