Modeling of GPR Clutter Caused by Soil Heterogeneity

In small-scale measurements, ground-penetrating radar (GPR) often uses a higher frequency to detect a small object or structural changes in the ground. GPR becomes more sensitive to the natural heterogeneity of the soil when a higher frequency is used. Soil heterogeneity scatters electromagnetic waves, and the scattered waves are in part observed as unwanted reflections that are often referred to as clutter. Data containing a great amount of clutter are difficult to analyze and interpret because clutter disturbs reflections from objects of interest. Therefore, modeling GPR clutter is useful to assess the effectiveness of GPR measurements. In this paper, the development of such a technique is discussed. This modeling technique requires the permittivity distribution of soil (or its geostatistical properties) and gives a nominal value of clutter power. The paper demonstrates the technique with the comparison to the data from a GPR time-lapse measurement. The proposed technique is discussed in regard to its applicability and limitations based on the results.


Introduction
During the last years, ground-penetrating radar (GPR) has been used more and more commonly for small-scale measurements in civil engineering and hydrology.Small changes in a medium need to be captured with a high sensitivity in such applications.Therefore, high frequencies (typically around 1 GHz or higher) are often employed.With increasing the frequency, GPR becomes more sensitive to heterogeneity of the medium surrounding objects, which results in unwanted scattering (commonly referred to as clutter) in the data.Clutter degrades the quality of GPR data and makes their analysis and interpretation difficult.Assessing the effectiveness of GPR by using rapid measurements and analysis in a site prior to the actual survey may help saving time and costs.
Studies on GPR and electromagnetic wave scattering in relation to random media have already been carried out in the past.These studies may be categorized into two types: numerical or analytical simulations of scattering for modeling the GPR response of heterogeneous media (e.g., [1][2][3][4]), and characterization of heterogeneous subsurface media by analyzing GPR data (e.g., [5][6][7][8]).The first type of studies requires a model of the heterogeneous medium.The simulated GPR response may be exact; however, exactly the same response cannot be obtained from different media.The latter type of studies has demonstrated the usefulness of the GPR data to characterize subsurface media in terms of dielectric and electric properties.A challenging point of studies on GPR with heterogeneous media is the fact that each heterogeneous distribution of soil is unique and exactly the same distribution cannot be found.Furthermore, to simulate full wave GPR response, three-dimensional distribution of electric, dielectric, and magnetic properties of soil need to be known which requires tremendous time and destructive measurements.In case of the detection of buried objects by GPR, the response of heterogeneous soil is considered geological noise (i.e., clutter) and it disturbs the target signature.The peak clutter power can be used to assess the detection performance in this case and the full wave information of clutter is not necessary.Therefore, a model that provides clutter power from soil heterogeneity has become of our interest, and simulating full wave GPR signature of heterogeneous soil is not the scope of the work.Random media in this regard may also include rough ground surface, however the work is focused only on the heterogeneous distribution of dielectric property of soil.
The authors developed a simple modeling method of GPR clutter caused by heterogeneous soils and demonstrated the effectiveness of the method by comparing its results with data acquired during an infiltration experiment [9,10].In these previous works, a model was constructed using parameters obtained by geostatistical analysis of GPR data, which were assumed to reflect soil heterogeneity.The model consists of a dielectric sphere, and the radar cross-section (RCS) of the sphere was calculated.The RCS is proportional to the backscattering power of the dielectric sphere and it is considered to be modeled clutter power.While the model calculation with the Rayleigh approximation showed significant difference from clutter actually observed in GPR data, the modeling with the Mie solution well agreed with the experiment.However, the agreement was achieved in a certain range of soil heterogeneity, and it was not in entire range.The results might indicate some limitations of the technique and the necessity of modifications for wider ranges of applicability and more accurate modeling.
The modeling method requires the parameters of soil heterogeneity, namely, correlation length and variability of the permittivity that can be determined by a geostatistical analysis of spatial distribution obtained by field measurements.In this study, the technique is demonstrated by other data sets that were repeatedly acquired during a few months on an outdoor test site.Time domain reflectometry (TDR) data acquired at the same time as the GPR measurements are used for modeling GPR clutter.Based on these modeling results, the applicability of the method and the characteristics of GPR clutter caused by natural soil heterogeneity are discussed in this paper.The experimental setup is almost the same as in the previously demonstrated studies [9,10].However, in the previous studies, the technique was demonstrated in the comparison with GPR data collected under irrigation of the ground which was an artificially produced change in water content and it can be seen unnatural.In this study, the changes in water content of the soils were only due to the seasonal changes of the hydrologic balance, that is, no artificial precipitations were applied.

GPR Clutter Modeling
2.1.Characterization of Soil Heterogeneity.The modeling technique requires input parameters that characterize soil heterogeneity and these parameters first need to be obtained.The easiest way to measure spatial variations of dielectric permittivity of soils in the field might be by carrying out TDR measurements at multiple locations on the ground surface.TDR gives the permittivity (or associated water content) of soil around the probes, and thus the measurement scale depends on the probe configuration.If relatively small probes (e.g., 10 cm in length and 2 cm in separation of the rods) are employed, the resolution may be similar to that of high-frequency GPR.Further, TDR typically uses frequencies ranging from 500 MHz to 1 GHz [11] that corresponds to relatively high frequency GPR.Therefore, TDR measurements can be used to measure and characterize soil heterogeneity in the similar spatial scale and frequency range as relatively high frequency GPR.
A semivariogram is a geostatistical analysis tool and it can be used to quantify heterogeneity.For soil permittivity data measured on a 1D profile, the semivariogram γ(h) is computed as [12] follows: where h is the lag distance, or separation between two data points, z(x i + h) and z(x i ), and N is the number of data pairs with a constant lag distance h from all data points.Often, the semivariance γ(h) increases with the lag distance h up to a certain value and then it becomes constant.The lag distance h and the semivariance γ(h) where γ(h) becomes constant are called range a and sill C, respectively.The range indicates the mean distance at which a data pair does not correlate anymore and thus it is equivalent to correlation length and characteristic length [13].The sill corresponds to the maximum variance within a data set and thus it is the indication of the variability.The exponential semivariogram model with no nugget effect given as [12] follows: is fitted to the obtained experimental semivariograms in order to determine range a and sill C. Since the exponential model asymptotically reaches its sill, the practical range a is defined as the distance where the model value reaches 95% of the sill [14,15].An example of the exponential semivariogram model is illustrated in Figure 1. this study is conceptually illustrated in Figure 2. It consists of a dielectric sphere embedded in a homogeneous space.The homogeneous background is defined having a dielectric permittivity equal to the mean permittivity, that is,

Dielectric
The circumference of the dielectric sphere is chosen to be the correlation length and thus the diameter d is defined as d = a/π, where a is the correlation length of the soil permittivity distribution determined by the geostatistical analysis as written in the previous section.The permittivity of the sphere ε 2 is set so that the difference to the ambient medium is equal to the square root of variability as follows: With this model, the radar cross-section (RCS) that is proportional to the backscattering power of the dielectric sphere is theoretically calculated, assuming a monostatic configuration and plane wave incidence.There are two ways to calculate RCS for the model: the Mie solution and Rayleigh approximation.The Mie solution is the exact solution of Maxwell's equations and describes the scattering by an arbitrarily sized particle [16].Therefore, it is valid in all frequency and sphere size ranges including not only the Mie scattering region but also the Rayleigh scattering and optical regions.By the Mie solution, the RCS σ s of a dielectric sphere is given as [17] follows: where x = kr, k is the wavenumber in the ambient medium, and r is the radius of the sphere.In case there is no change in the magnetic permeability between the dielectric sphere and the ambient medium which is assumed in this study, the coefficients a n and b n are given by the following: where m denotes the refractive index.The j n (x) and h (1)  n (x) are the spherical Bessel function of first kind of order n and the spherical Hankel function of order n, respectively.Primes mean derivatives with respect to the argument.The summation in ( 4) is theoretically from n = 1 to ∞, but the infinite series can be truncated after n max = x + 4x 1/3 + 2 [17].
The Rayleigh approximation describes scattering by a small particle compared to a wavelength [17], and it is generally valid only in Rayleigh scattering region.In such a case (x 1), the Mie solution ( 4) and ( 5) can be simplified by the Rayleigh approximation as [17] as follows: The equation exhibits that the scattering power due to the particle is inversely proportional to the fourth power of the wavelength (or proportional to the fourth power of the sphere size) in Rayleigh scattering region.In this study, both Mie solution and Rayleigh approximation were examined.

Repetitive GPR Measurements
GPR and TDR measurements were carried out on an outdoor test site at the Leibniz Institute for Applied Geophysics (LIAG) in Hannover, Germany.The soil at the test site is a natural soil material that was developed in postglacial sedimented aeolian sand.These deposits cover the quaternary sediments and form the uppermost layer in wide areas of North German lowland plain.The texture of the soil is medium sand with 1.0, 6.7, and 92.3% of mineral soil of clay, silt, and sand contents, respectively.The soil also has high humus content that is 6.3%.At the time of data acquisition, the ground surface was covered by grass, which was considered to increase the heterogeneity of soil moisture distribution and the associated dielectric permittivity distribution due to the water consumption by the irregular root system.A GPR system (GSSI SIR-3000) with 1.5 GHz antennas was employed.GPR data were collected on the ground along a 5 m long profile.At the same time as the GPR measurements, the dielectric permittivity of the topsoil was measured by a TDR (FOM/mts, Institute of Agrophysics of the Polish Academy of Sciences) on the same profile.The TDR uses two 10 cm-long stainless probes with the separation of 1.6 cm.GPR scanned the full 5 m long profile and sampled data every 1 cm, while TDR measurements were .Travel time was converted to depth using mean permittivity measured by TDR (Figure 4).Amplitudes of the profiles were normalized throughout the measurements.Only a bandpass filter was applied.
carried out every 2 cm along only the first 1 m of the profile.The measurements were repeated seven times from April to June 2010 with the same equipment, under the same set up and on the same profile.
Figure 3 shows all seven GPR profiles.In order to keep the original characteristics of raw GPR data, no advanced signal processing but only a bandpass filter for the removal of DC and high frequency noise components was applied to the data.The travel time was converted to depth using the mean permittivity measured by the TDR in each case (shown in Figure 4) and the amplitudes were normalized to the maximum throughout all the measurements.Almost the same patterns can be observed in all seven radar profiles.Only small changes in pattern and intensity differences are observed.Figure 4 shows the spatial change of relative permittivity measured by TDR at the same time with the GPR measurements.The associated volumetric water content was given by Topp's equation [18].It can be observed that the water content and permittivity changed due to evapotranspiration and some precipitation events.It has to be noted that the site was left completely natural during the measurements (i.e., without any water artificially applied).

GPR Data
4.1.Geostatistical Analysis of the TDR Data.Experimental semivariograms were calculated according to (1) for the TDR permittivity data as shown in Figure 5(a).By semiautomatically fitting the exponential model with no nugget effect given by ( 2) to the experimental semivariograms in the least square manner, the correlation length and variability of dielectric permittivity distributions were obtained as shown in Figure 5(b).The result showed that the correlation length and variability changed in the ranges of 10-25 cm and 0.5-1.5, respectively, and thus the heterogeneity of soil permittivity had changed naturally during the experiments.

Clutter Modeling.
The determined correlation length and variability as well as the mean permittivity were set into the model shown in Figure 1 and clutter power was calculated with this model as described.In this study, both Mie solution and Rayleigh solution were examined and the calculations are shown in Figure 6.Modeled clutter power was normalized by the mean for comparison.In order to confirm the modeling, clutter observed in the GPR data was extracted for the comparison with modeled clutter.The 20 highest amplitudes in a depth section deeper than 7.5 cm that was just below the ground surface reflection were picked from the GPR data shown in Figure 3.The extracted clutter power and their mean values are plotted with crosses in Figure 7. Clutter modeled with the Mie solution is similar to the mean of the picked clutter power.They do not perfectly fit, but the tendencies (ups and downs) agree.On the other hand, the modeled clutter with the Rayleigh approximation does not really fit to the experiments.Mean squared errors of modeled clutter with the Mie solution and Rayleigh approximation are 0.15 and 0.87, respectively.The result indicates that the power of clutter signals caused by soil heterogeneity can roughly be calculated by the Mie solution for a rather simple model, shown in Figure 2, with the support of measurements of soil permittivity spatial distribution (i.e., TDR measurements) and a geostatistical analysis.

Discussion and Conclusions
In this study, GPR clutter was modeled from the soil heterogeneity that was characterized by the geostatistical analysis of TDR measurements.The modeling results agreed with clutter power observed in the GPR data.
The method could reasonably model the GPR clutter variation caused by the change in soil heterogeneity over a few months.However, modeling did not fit perfectly to the experiments.Especially, the errors for the first (on April 15th), fourth (on May 5th), and sixth (on June 3rd) measurements are relatively large as shown in Figure 7.The possible reason can be found in relative correlation length (Figure 8) and variability (red circles in Figure 5(b)).The results on these dates showed longer correlation lengths and larger variability.The second measurements on April 21st showed a longer correlation length, but variability was not large and thus the modeling error is not large.To illustrate the modeling error in relation to the geostatistical properties of soil heterogeneity, the errors are plotted with respect to correlation lengths relative to wavelength multiplied by variability (i.e., Ca/λ) in Figure 9.It is clearly seen that the modeling error becomes larger as the quantity increases.The observation also supports the results of our previous work; the modeling did not fit well to the experiments immediately after irrigation where long correlation lengths and large variability were observed [9,10].A possible explanation is that a longer correlation length may indicate gradual changes of permittivity in space that contradict a sharp boundary assumed in the dielectric sphere model.Moreover, a larger variability may cause larger scattering and thus more complex wave interactions in the soil and they are not taken into account in the simple modeling.
Comparing the modeled (Figure 6) or measured clutter (Figure 7) to the determined correlation length and variability of dielectric permittivity of the heterogeneous soils (Figure 5(b)), a simple relationship between these parameters and clutter power cannot be found.This is because the scattering by heterogeneous soils is in the Mie scattering region as the relative correlation lengths are larger than 1 (Figure 8) and suggested by [19].The backscattering power oscillates as correlation length varies in this region and thus a simple relationship cannot be found.If the scattering was dominated by Rayleigh scattering, clutter power and these parameters characterizing the soil heterogeneity should have a direct correlation since scattering in the Rayleigh scattering region can be expressed by a relatively simple form as given in (6).In fact, the modeled clutter power with the Rayleigh approximation is much larger than that with the Mie scattering if they are not normalized as shown in Figures 6 and 7. Since the Mie solution is also valid in the Rayleigh scattering region, both calculations should give similar values if the scattering problem was in the Rayleigh scattering region.Therefore, it is evident that the scattering caused by soil heterogeneity, which is observed as clutter and represented by the dielectric sphere model in this study is dominated by Mie scattering.This may be the case that a relatively higher frequency (e.g., higher than 1 GHz) is employed.It is different to the situation with a lower frequency, which is commonly used for large-scale measurements.In the case of a lower frequency GPR, the scattering due to heterogeneous media may be governed by Rayleigh scattering [11].

Figure 1 :
Figure 1: An example of the exponential semivariogram model with the practical range a = 10 cm and sill C = 0.7.

InternationalFigure 3 :
Figure 3: GPR profiles on a 5 m long profile acquired on (a) April 15th, (b) April 21st, (c) April 28th, (d) May 5th, (e) May 26th, (f) June 3rd, and (g) June 17th.Travel time was converted to depth using mean permittivity measured by TDR (Figure4).Amplitudes of the profiles were normalized throughout the measurements.Only a bandpass filter was applied.

Figure 4 :Figure 5 :
Figure4: (a) Spatial variation of the relative permittivity of topsoil up to 10 cm depth measured by a time domain reflectometry (TDR) on the first 1 m of the profile with 2 cm spacing and (b) the mean permittivity of the TDR measurements.Volumetric water content is calculated from permittivity by Topp's equation[18].

Figure 6 :Figure 7 :
Figure 6: Modeled clutter power using the Mie solution (blue dots) and Rayleigh approximation (red circles).Values are normalized by the mean.

Figure 8 :Figure 9 :
Figure 8: Correlation length relative to wavelength determined from the semivariograms.