Abstract

Effective diagnostics with ground penetrating radar (GPR) is strongly dependent on the amount and quality of available data as well as on the efficiency of the adopted imaging procedure. In this frame, the aim of the present work is to investigate the capability of a typical GPR system placed at a ground interface to derive three-dimensional (3D) information on the features of buried dielectric targets (location, dimension, and shape). The scatterers can have size comparable to the resolution limits and can be placed in the shallow subsurface in the antenna near field. Referring to canonical multimonostatic configurations, the forward scattering problem is analyzed first, obtaining a variety of synthetic GPR traces and radargrams by means of a customized implementation of an electromagnetic CAD tool. By employing these numerical data, a full 3D frequency-domain microwave tomographic approach, specifically designed for the inversion problem at hand, is applied to tackle the imaging process. The method is tested here by considering various scatterers, with different shapes and dielectric contrasts. The selected tomographic results illustrate the aptitude of the proposed approach to recover the fundamental features of the targets even with critical GPR settings.

1. Introduction

In a large variety of civil, forensic, geophysical, and planetary applications, ground penetrating radar (GPR) is used in standard configurations close to a ground interface with the aim of locating and imaging shallow-buried targets [14]. The reliability of the reconstructed images to provide information on the geometrical and physical features of the targets (i.e., location, size, shape, and contrast) depends both on the amount and quality of the data available with GPR systems [19] and on the adopted imaging procedures, mostly based on the solution of inverse scattering problems [14, 1021].

In this context, the study presented here is focused on the assessment of the GPR technique to provide valid imaging performance, specifically referring to the critical problem of reconstructing targets buried in shallow subsurface. These objects can have different shape and dielectric contrast; their typical dimensions can be compared to the wavelengths of the propagating GPR signal, and their location is not necessarily in the far-field region with respect to the antenna.

Full three-dimensional (3D) target reconstructions will be attempted by gathering scattering “synthetic” numerical data from proper sets of GPR “B-scan” analysis (“radargrams”). The investigated scenarios will take into account several operational factors that affect the quality of imaging, such as the limited signal spectrum, the realistic antenna system, the target size comparable to the resolution limits of the instrument, and the near-field illumination. The detailed description of the numerical setup for the analysis of the scattering problem and the relevant discussion of the synthetic data achievable from radargrams will be given in Section 2.

These GPR scattering data will then be processed by means of a tailored microwave tomographic approach that faces the imaging as a suitable solution of a linear inverse scattering problem. The distinctive features of the inversion method will be outlined in Section 3. The relevant results of tomography from 3D imaging, presented and discussed in Section 4, will illustrate the quality of the information recoverable on the geometrical and physical features of the targets, for canonical sets of parameters (shape, size, location, and contrast). Conclusive remarks will be outlined in Section 5, in the frame of issues and perspectives of this class of problems.

2. Numerical Setup and Synthetic Radargram Data

In order to have a wide set of data to be efficiently processed for imaging purposes, the “forward” scattering problem is analyzed here with a numerical implementation able to simulate a typical GPR environment of interest. A recently developed setup (whose main details can be found in [8]), based on the electromagnetic time-domain computer-aided-design (CAD) tool CST Microwave Studio [22], enables us to explore different realistic complex GPR scenarios with affordable computational efforts and excellent accuracy.

The investigated 3D domain here consists of two half-space media separated by a flat interface, representing an air/soil environment, in which various types of scatterers can be buried. The GPR system can be simulated with various types of antennas, both in bistatic and in monostatic configurations. In order to keep the GPR system as simple as possible in the analysis of the imaging performance, in the present study a monostatic configuration is chosen, with the antenna placed at the interface between air and ground. A specific wideband printed monopole antenna is designed (reference details on the procedure followed can also be found in [8]), whose operational frequencies are quite similar to commercially available GPR systems [8, 9]. The numerical GPR environment under analysis is shown in the 3D view of Figure 1(a). The antenna scanning domain at the interface is sketched in the top view of Figure 1(b), with the numbering of the scanning lines, the grid points for gathering synthetic data, and the proper coordinate system.

In accordance with the specifics of a GPR instrument such as PulseEkko Pro by Sensors and Software [23] (already used in a number of experimental studies [8, 9, 21]), the simulated system is fed with an input Gaussian-like signal having a spectrum between 0.5 and 1.5 GHz, generated as the modulation of a pulse through a 1 GHz sinusoid; this waveform is shown by the time-domain trace of Figure 1(c), with the relevant frequency spectrum of Figure 1(d). The matching characteristics of the antenna in the considered operational range are summarized through the return loss (RL), that is, the magnitude of the antenna scattering reflection coefficient in dB, versus frequency . It is seen from Figure 1(e) that the RL curve for our antenna in realistic operating conditions (i.e., located at the interface of the air/ground external environment) lies in the range of interest well below the typical matching threshold of −10 dB, as desired.

The far-field radiative features of the antenna, again operating in the air/ground environment, are described by means of the radiation patterns. The pair of examples of Figure 1(f) show the polar-form plots versus the elevation angle , for the two principal planes, or xz plane (left) and or plane (right), at the central frequency of the spectrum ( GHz). In both cases it is noticed that, as expected, radiation is focused more in the ground dielectric region and less in the air, with a certain forward squint of the main lobe (around 45° in the plot) due to the asymmetry of the radiating printed element with respect to the central plane .

Since in general the GPR system will not operate in far-field conditions, it is also important to have precise information on the actual radiation in the near field. The CAD tool can compute the distribution of the fields in the volume of interest close to the antenna at the interface. A section view in the plane ( fixed) of the magnitude of the electric field for our antenna in the two half-spaces (air/ground) environment is illustrated by the color plot of Figure 1(g), for a fixed frequency (again at the central frequency  GHz). The spatial focusing of the near-field radiation is clearly emphasized, which is placed at an angle, as already observed in the far-field radiation patterns as well. This piece of information on the field distribution is particularly significant since in our cases the shallow scatterers will just be located in such near-field region.

Our environment is then completed by the introduction of suitable targets that are buried in the ground medium. In the forward problem, shape, size, depth, and electromagnetic parameters of the scatterer are fixed. In each simulation of the CAD tool, after the signal is launched by the antenna, the scattering effect of such a target is evaluated by means of the time-domain trace of the received signal, which is due to the scattering from the buried object. (It is noted that our choice for an interfacial antenna allows us to eliminate most of the signal reflection due to the media discontinuity; in addition to this, absorbing boundary conditions have been implemented in order to avoid any further reflection effect from the surfaces enclosing the 3D investigation domain.) The numerical evaluations of the received signal traces are repeated in a multimonostatic configuration, that is, for a number of positions of the antenna with respect to the fixed target inside a 2D rectangular domain (area on the interface plane at ), scanning along the direction for different parallel lines by changing , in a way similar to the actual measurements performed by a GPR instrument. In the numerical results presented here the scanning step along is fixed at 3 cm; the number of points along is , for an overall investigated length of 81 cm. The distance between parallel lines along is also fixed at 3 cm; the number of lines along is , for an overall width of 21 cm (refer to the sketch of Figure 1(b) for the investigated grid, the relevant numbering of the lines, and the antenna position). Each time-domain trace of the received signal is therefore available in the points of a grid with unit cell  cm (the relevant distance between grid points compared to the free-space wavelength at the central frequency being equal to 0.1).

As is typical in GPR instruments, the results of the time-domain traces giving the amplitude of the scattered wave received by the antenna can be plotted in the form of B-scan radargrams (time delay of the received signal versus scan position, in grey-scale form for its intensity) [1, 2, 8, 9]. The environment in the simulations presented here is chosen as air (vacuum) for the upper medium and as a “dry sandy soil” (, , ) for the lower medium [8, 9, 21].

In order to emphasize the effects of different radar cross sections, the buried scatterers are chosen of cubic and spherical shapes and can have different electromagnetic contrasts: specifically, computation is presented here for (a basalt block) and for (air cavity).

The effects of the dimensions can be also analyzed, for example, with possible different cubic sides or sphere diameters. Finally, also the influence of the scatterer location can be considered, choosing variable depths under the surface ( being the distance of the top point of the target from the air/sand interface). We remind that in our conditions the characteristic dimensions of the scatterers are generally comparable to the wavelengths of the signal (of the order of 10 cm) and are illuminated in near-field conditions.

First examples of the GPR synthetic radargrams are shown in Figure 2 for a cubic basalt scatterer of 10 cm side and buried 10 cm deep in the sand: specifically, Figures 2(a), 2(b), and 2(c) refer to B-scans versus for three chosen lines along (lines numbers 1, 5, and 8). The scattering events, emphasized by the received-signal amplitudes in grey scale, present expected hyperbolic shapes along . It is seen that, considering also the squint effect of the monopole antenna and its relevant orientation in the scanning process (see Figure 1(b)), the reflection effect correctly appears to be in general lower for the upper positions (i.e., line 1), where radiation is not properly focused on the target, and sensitively increases as the antenna scans towards the center (i.e., line 5), where the scattering effect is maximal; at lower positions (i.e., line 8) the effect is higher with respect to the upper position (line 1) due to the asymmetric squinted illumination (in accordance with the field distributions reported in Figures 1(f) and 1(g)). Moreover, it is observed that the position in time of the main upper and lower hyperbolic events is fully consistent with the traveling speed of the signal scattered by the upper and the lower sides of the penetrable object. Also, the echoes from the lateral sides of the block are present, even if weaker than those from its top: this is related to the fact that the scattering from the top of the cube is mainly due to strong reflection by a flat interface and edge points, while the scattering from the side walls is related to diffracted fields spread over wider angles.

In order to illustrate the specific differences that can be observed in the synthetic radargrams depending on the shape of the scatterer, Figure 3 presents the results for a basalt sphere of radius 5 cm, all the other parameters being fixed as in Figure 2. Again, Figures 3(a), 3(b), and 3(c) show B-scans for parallel lines with different (lines numbers 1, 5, and 8). By using the same grey scale as in Figure 2 for the received signal intensity, it is clearly seen that the reflection events from the sphere are generally less strong than the cube, as expected from the reduced scattering effect of the round shape with respect to the flat and wedged shapes. The different radar cross section of the two geometries affects in part also the time delays in which the maxima of the scattered waves are found [8, 9, 21].

The effect of different electromagnetic contrast has also been evaluated through the radargram forms. The results of Figure 4 refer to the case of a cube as in Figure 2 (same dimensions and location) but having permittivity of a vacuum (representing an air cavity in the sandy soil): again, Figures 4(a), 4(b), and 4(c) show the B-scans for three different lines (numbers 1, 5, and 8), with the relevant amplitudes in grey scale. Even if in this case the scattering effect is related to a target which is less dense than the outer host environment (sand), the radargrams of Figure 4 appear quite similar to the ones of Figure 2. Actually, from a careful analysis of the single-trace numerical data (not reported explicitly here for brevity), it is seen that the main reflection hyperbolic event, due to the upper side of the cube, is in general slightly stronger for the vacuum cube than for the basalt one. This is related to the fact that the reflection amounts in the hosting medium (sand) are different for the two scattering materials, being greater for vacuum than for basalt: in fact, for the sand/basalt case the magnitude of the reflection coefficient for a normally-incident plane wave is about 0.35, while for the sand/air case the magnitude is about 0.5. Furthermore, the thickness of the scattered hyperbola appears to be different for the two objects (greater for the basalt than for the vacuum cube), due to the propagation velocity of the wave transmitted through the penetrable scatterers ( being lower in the basalt than in the vacuum), which generates a dissimilar composition with the wave scattered from the bottom of the target.

The results of Figure 5 refer to the case of a sphere having permittivity of a vacuum, buried in the sand: dimensions and location are the same as in Figure 3. Figures 5(a), 5(b), and 5(c) show the B-scans for the usual different lines (numbers 1, 5, and 8), with the same grey scale of Figure 4. It is seen that the distinctive features already presented for the vacuum cube can be applied to the vacuum sphere as well. Even if the scattering effects are confirmed to be less significant for the sphere if compared to the equivalent cube, it is again observed that the main reflection hyperbolic event is in general slightly stronger for the vacuum sphere than for the equivalent basalt one.

From all the shown synthetic radargrams (Figures 2 to 5) it is seen that certain amounts of spurious signal occur (“flat” events found around 2 ns), due to slight reflections at the input port of the antenna, whose matching features in the operational bandwidth are affected also by the presence of nearby dielectric discontinuities. These spurious contributions are particularly visible when the basalt scatterer is in the closest position with respect to the antenna, whilst the undesired signals become almost invisible particularly for the vacuum targets located far from the antenna (see, e.g., Figure 5, specifically (a) and (c) cases).

From this analysis it is seen that all such synthetic data concerning the forward scattering problem appear to be pretty consistent and regular, so that they are particularly suitable to be processed by an efficient inversion algorithm, as discussed next.

3. Full 3D Microwave Tomographic Approach

The simulated GPR data have been processed by means of a “microwave tomographic approach” capable of providing full 3D images of buried objects from data gathered on a planar surface, that is, along several parallel linear traces. Such an approach basically exploits the Born approximation to model the underlining scattering phenomenon [16]. Therefore, the imaging is faced as the solution of a linear inverse scattering problem, wherein the data/unknown relationships account for the dyadic nature of the interaction between the radiation and the probed materials. In particular, the surveyed environment is assumed as made of two homogeneous half-spaces separated by a planar interface (, being the main horizontal scanning axis, as already introduced in the previous analysis). The upper half-space is air, while the lower half-space represents the medium hosting the objects, and it is assumed to be homogeneous with known relative permittivity and conductivity . In the following, these parameters are referred to the already described ground sandy material. The antenna system works in a multimonostatic measurement configuration, and it is moved with a uniform step along evenly spaced linear traces, which are parallel to the air/soil interface and span a planar domain Γ at height above the ground. The probing signal is modelled as the time-harmonic electric field radiated by an elementary electric dipole oriented along the -axis and operating in the frequency range [, ]. It is worth noting that this configuration idealizes the position of the GPR dipole antennas having their major length just oriented along the -axis (similarly to the CAD implementation). The processed data are the component of the electric field backscattered by the unknown targets as evaluated at measurement points on Γ and discrete frequencies uniformly spaced in the range [] of the chosen signal (0.5–1.5 GHz here). The targets to be reconstructed are enclosed in a parallelepiped domain below the air/ground interface.

According to the previous assumptions, the data/unknown relationship is expressed theoretically through the integral equation [24]: where is the component of the field scattered by the targets, is the wavenumber in air, G() is the dyadic Green function, that is, the electric field radiated on Γ by an arbitrary oriented electric dipole located in , and () is the incident field in due to a -oriented electric dipole. Moreover, , , and specify the position of the source, of the receiver, and of a generic point in , respectively. The quantity is the unknown “contrast function,” relating the permittivity of the target to that of the hosting medium. It is worth pointing out that, since the antenna system is very close to the air/ground interface, the Green function as well as the incident field is herein approximated as those referred to a homogeneous medium having the electromagnetic properties of the ground hosting the targets. Their analytical expression can be found in [24].

Equation (1) states a linear integral relationship between the meaningful signal and the unknown function , which is expressed through the compact operator , where the subscript is to underline that both the transmitting and receiving antennas are here placed along the -axis.

The imaging problem being formulated as the inversion of the relationship in (1), a linear inverse problem must be carefully handled. This task is herein faced by using the truncated singular value decomposition (TSVD) as a regularized inversion scheme. Therefore, a stable regularized solution is obtained as [11] where is the -dimensional data vector, is the singular system of the matrix arising from the discretization of the compact operator in (1), being equal to and denoting the number of voxels used to discretize , denotes the scalar product in the data space, and is the truncation threshold. The choice of the index is performed in order to ensure a trade-off between the contrasting needs for accuracy and resolution from one hand (which should push to increase such an index) and for the stability of the solution from the other hand (which should push to limit the increase of the index) [21]. The “regularized reconstruction” is then “edged” to limit spurious effects due to noise and smoothing introduced by regularization.

4. Results of 3D Imaging from GPR Synthetic Data

The results of the microwave tomographic approach outlined in Section 3 are derived here by processing the set of forward synthetic data discussed in Section 2. The following images represent the amplitude of the reconstructed contrast function as provided by (2) and normalized with respect to its maximum value inside the investigation domain . Such a domain has size 0.8 m × 0.4 m × 0.5 m, and it is discretized by means of = 13090 cubic cells, whose side is 2.24 cm long. The TSVD truncation threshold is fixed in such a way to filter out all the singular values which are 20 dB lower than the maximum one.

The first processed data sets are referred to the buried basalt cube and basalt sphere already analyzed in Figures 2 and 3, respectively. The depth slices of the reconstructed contrast functions in color-plot form (related to its intensity according to a proper scale from red to blue) are shown in Figure 6 for the cube and in Figure 7 for the sphere. These figures present 20 planar () color plots corresponding to values which regularly increase moving from the interface towards the bottom of the ground medium, with a spatial step of 2.24 cm, for an overall investigated height of about 50 cm.

The analysis of Figures 6 and 7 shows that the two dielectric objects are approximately located at the same depth, which is properly retrieved, and allows us to obtain a satisfactory estimation of the upper cross-section size of both the targets, while their size along the depth direction is overestimated. This agrees with the fact that the objects are characterized by a larger relative permittivity than the hosting medium and thus by an actual wave propagation velocity that is lower than the modeled one (this latter being that corresponding to the background [18]). Moreover, even if the actual shape is not accurately reconstructed and no information on their relative permittivity is available, from Figures 6 and 7 it is possible to infer that the objects have different shape and are made by the same material. As a matter of fact, while the reconstructed cross sections for the cube do not significantly vary with ranging from about 12 cm to 16 cm, this correctly does not happen for the sphere (in the latter case the radius of the circular scattering section changes as varies). In addition, both the reconstructions exhibit similar features as far as their behavior along the axis is concerned.

To corroborate the previous considerations, Figures 8 and 9 show the depth slices of contrast functions retrieved by processing the data sets corresponding to vacuum cubic and spherical objects (with data analyzed in Figures 4 and 5), respectively, which have the same size and location of the already-considered basalt ones. These figures show that the targets depth is again properly retrieved, as well as the size of their upper cross section. Moreover, it is evident that, as in the case of basaltic objects, the reconstructed cross sections of the cubic object do not significantly change for ranging from about 12 cm to 16 cm, while this is not true for the sphere. Furthermore, the reconstructions in Figures 8 and 9 have a similar behavior along the axis, which is different from that shown in Figures 6 and 7. In this respect, it is worth nothing that now the targets size along the depth direction is underestimated in agreement with the lower relative permittivity of the objects than the one of the hosting medium.

5. Conclusion

Information retrieval by operating with canonical GPR systems has been investigated for cases of remote sensing which find application in several practical diagnostic problems. A full 3D imaging has been attempted for dielectric targets, which can have dimensions comparable to the typical wavelengths of the signals and are buried in the shallow subsurface in near-field conditions. The forward scattering data have been obtained with a tested numerical approach considering a variety of realistic scenarios simulating typical multimonostatic GPR operation with proper antenna, signal waveforms, media environment, and targets.

An efficient tomographic algorithm for the inverse problem, suitably tailored for our GPR setup, has successfully been applied by processing of the various synthetic data of radargrams. A consistent analysis of the imaging results for canonical buried dielectric scatterers has thus been possible.

It has been noted that intrinsic limits related to the adopted simulation configurations cannot allow for the identification of very fine geometrical details. This is related both to the procedure of gathering data, which derive from a two-dimensional and finite scanning domain of the monostatic GPR instrument, and to the fixed range of the signal frequency spectrum. Moreover, the angular spread of the radiated field by the GPR antenna, particularly significant when operating in the near-field regions, is an additional parameter which adversely affects the resolution features of the procedure.

Despite these aspects, it has been quantitatively tested that, starting from an accurate set of forward data and using a stable reconstruction algorithm for the inversion, in various typical conditions this approach possesses a good capability to prove the correct spatial position of the reflecting parts of the targets, in conjunction with a quite satisfactory prediction of the shape and size (cross section).

Acknowledgment

The authors gratefully acknowledge financial support from the Italian Space Agency through Contract ASI I/060/10/0 EXOMARS SCIENCE phase C2/D.