Ground Penetrating Radar (GPR) is a well assessed diagnostic instrumentation that finds application in a very large range of sectors where the aim is of achieving information about the presence, the position and the morphological properties of targets embedded in opaque media is c crucial task. Therefore, it is of interest in many fields which range from civil engineering diagnostics (Hugenschmidt & Kalogeropulos, 2009; Soldovieri et al. 2006), to archaeological prospecting (Conyers & Goodman, 1997; Piro et al., 2003; Orlando & Soldovieri, 2008), to cultural heritage management, to geophysics, only to quote a few examples.
Besides these classical fields of application, GPR methodology/technology is arising interest in new fields related to the security framework such as: through-wall imaging (Baranoski, 2008, Ahmad et al. 2003; Solimene et al., 2009), tunnel and underground facility detection (Lo Monte et al., 2009), borders surveillance.
Another field of recent interest is concerned with the exploration of planets; in particular, satellite platform radars are just investigating the Marsis subsurface (Watters et al., 2006) for the water search and other similar missions are under development for Lunar exploration.
GPR is based on the same operating principles of classical radars (Daniels, 1996). In fact, it works by emitting an electromagnetic signal (generally modulated pulse or Continuous Wave) into the ground; the electromagnetic wave propagates through the opaque medium and when it impinges on a non-homogeneity in the electromagnetic properties (representative of the buried target) a backscattered electromagnetic field arises. Such a backscattered field is then collected by the receiving antenna located at the air/opaque medium interface and undergoes a subsequent processing and visualization, usually as a 2D image.
Anyway, significant differences arise compared to the usual radar systems.
The first one is that while usually radar systems act in a free-space scenario, GPR works in more complicated scenarios with media that can have losses and exhibit frequency-dependent electromagnetic properties thus introducing dispersive effects.
Secondly, while the range of radars may be also of hundred of kilometres, in the case of GPR we have ranges, also in the most favourable cases, of some meters due to the limited radiated power and the attenuation of the signal.
Moreover, GPR resolution limits, which ranges from some centimetres to some metres (in dependence of the working frequency), usually are smaller as compared to the ones of the usual radar systems.
According to the working principles described above, the architecture of the radar system can be schematised very briefly by the following blocks:
Electronic unit that has the task to drive and command the transmission and the reception of the antennas and to send the collected data to a monitor or to a processing unit;
Two antennas that have the task of radiating the field impinging on the target and collecting the field backscattered by the target, respectively;
A monitor/processing unit that has the task of the subsequent processing and visualization of the measurements collected by the receiving antenna.
From the schematization above, it emerges that a first advantage of the GPR instrumentation resides in the moderate cost and easiness of employ; in fact no particular expertise is required to collect the data. Secondly, the instrumentation is easily portable (unless very low frequencies are exploited thus increasing the physical size of the antennas) and allows to survey regions also of hundreds of square metres in reasonable time. Finally, the flexibility of the GPR system is ensured by the adoption of antennas (mostly portable) working at different frequencies and that can be straightforwardly changed on site.
It is worth noting that the necessities of probing lossy medium and of achieving resolution of the order of centimetres poses a very challenging task for the antennas deployed in the survey. On the other hand, as mentioned above, the portability of the system has to be ensured so that no complicated and large antenna systems can be employed. In particular, the trade-off between the necessity to have a large investigation range (that pushes to keep low the operating frequency) and the aim of achieving good spatial resolution makes the overall working frequency band exploited by GPR systems ranging from some tens of MHz to some GHz.
Therefore, the antennas designed for GPR applications have to comply with the following requirements: small dimension, ultra-wide-band behaviour and low levels of the direct coupling between the transmitting and receiving antennas.
The most widespread antennas employed for GPR surveys fall within the class of electrically small antennas such as dipole, monopole or bow-tie antennas.
Since these antennas are designed electrically small, the radiation pattern is broad and so low gain is achieved; keeping broad the radiation pattern allows to exploit synthetic aperture antenna which are particularly suitable to obtain focussing effect in the reconstructions. In addition, these simple antennas radiate a linear polarisation but also exhibit the drawback that relatively small frequency bandwidth is achieved. The working frequency band can be enlarged by loading techniques and the most-widely used technique is based on the application of resistive loading (Wu & King, 1965). The main disadvantage of resistive loading is the considerable reduction of the radiation efficiency due to the ohmic losses. In the last years, some approaches exploiting reactive elements to achieve the loading effect have been presented. A very recent development is concerned with the loading approaches based on the combination of resistive and capacitive loads; this offers the best performances in terms of trade off between small late-time ringing (increase of the frequency bandwidth) and high radiation efficiency (Lestari et al., 2004).
In addition, very recent advances are in course for adaptive (Lestari et al., 2005) and reconfigurable antennas, already exploited in communication field. In particular, the development and the implementation of RF/microwave switches have permitted the birth of a new concept of the antenna as a device that can dynamically adapt its behaviour to different measurement situations and operational contexts thanks to the change of its main parameters such as: input impedance, frequency band, radiation pattern. A design of an antenna for GPR purposes, with a geometry completely reconfigurable within the frequency band 0.3–1 GHz, and its performances compared to a reference bow-tie antenna are reported in (Romano et al., 2009). The structure has been “synthesized” by means of the total geometry morphing approach that represents the most structurally complicated but also versatile design method exploited to achieve the reconfigurability of an antenna.
Another class of antennas very widespread is the frequency independent one made up of spiral antennas (equiangular and conical). This kind of antennas are characterised by a geometry that repeats itself and in most cases the geometry (as in the case of spiral antennas) is completely described by an angle. Besides the frequency independent behaviour, spiral antennas are able to radiate a circular polarization; this is useful when the aim of the survey is to achieve information on objects that are elongated along a preferential direction. Another kind of frequency independent antennas is the Vivaldi antenna where the travelling wave energy is strongly bounded between conductors when the separation is very small, while becomes progressively weaker and more coupled with the radiation field when the separation increases. As main features, the Vivaldi antenna has an end-fire radiation mechanism, and the beamwidth is about the same in both planes E and H; the gain is proportional to the whole length and to the rate at which the energy is leaked away.
2. The measurement configuration and the radargram
Different configurations are adopted in GPR surveys and their choice is dictated by different applicative motivations such as: the type of investigation to be made; the kind of targets to be detected; the extent of the investigated region.
In general, all the reasons above push towards the adoption of the simplest acquisition mode for which the GPR system works in a monostatic or bistatic configuration. In the former case, the locations of the transmitting and the receiving antennas are coincident (or very close in terms of the radiated wavelength), whereas in the second case the transmitting and the receiving antennas are spaced by a fixed offset that is constant while they move along the survey line.
By moving the antenna system along a selected profile (line) above the ground surface a two-dimensional reflection profile (radargram) is obtained in which, for each location of the antenna system, a trace is achieved where the amplitude and the delay time of the recorded echoes (that can be related to the depth of the underground reflectors) is drawn (Daniels, 1996).
Such a radargram provides a rough information on the presence and the location of the targets but the actual shape of the target is blurred due to the effect of the propagation/scattering of the electromagnetic wave in the soil.
Such a statement is made clearer by considering the scattering by a point-object (an object small in terms of the probing wavelength) that is imaged as a hyperbola in the radargram. In fact, by denoting with x the position of the TX/RX system along the measurement line and with (0,d) the position of the point-object, the two-way travel-time is given by
where is the velocity of the electromagnetic wave in the soil. Therefore, the recorded data is represented as a hyperbola with a vertical axis and the apex at (0, 2d/v). Figure 1 is a pictorial image of the consideration above.
As a result, the shape of the hyperbola depends on the electromagnetic properties of the investigated medium (electromagnetic velocity), the configuration (monostatic or bistatic) and on the depth of the scatterer.
Besides the configuration said above, others acquisition modes can be deployed in the survey. Among them, we recall the common midpoint (CMP) where the transmitting and receiving antennas are further spaced during the survey but leaving the same the midpoint between them. The adoption of this sounding configuration is particularly useful when one aims at determining the velocity of the electromagnetic wave in the medium (Huisman et al., 2003).
Further configurations are also exploited, but they have the drawback that are not particularly simple, as the multi-fold ones where one transmitting antenna is located at a fixed position while the backscattered field is collected by the receiving antenna moving at different locations. The same measurement configuration can be further complicated by repeating the measurement procedure for different location of the transmitting antennas.
All the configurations presented above are concerned with a reflection mode where both the transmitting and receiving antennas are from the same side with respect to the investigation domain (as for example in the case of subsurface prospecting).
A different mode can be considered for the case of the masonry structure as the trasillumination or transmission mode where the transmitting antenna and the receiving one are at the opposite sides of the structure under investigation.
3. GPR subsurface imaging
The aim of subsurface imaging is to obtain images of a scene which is buried in the soil in order to reconstruct non-homogeneities perturbing the background medium.
This can be achieved at different degree of difficulty depending on what one may want to retrieve. For example, one can be interested in detecting the buried scatterers “only” or even in determining their geometrical sizes or, finally, the nature of scatterers can be of concern.
From a mathematical point of view, subsurface imaging is a special case of inverse scattering problem where, at microwave frequencies, the scatterers are reconstructed in terms of a spatial map of their dielectric and/or conductivity features.
Even though much work has been done in this field non-linear imaging schemes still suffer from reliability problems due to the false solutions (i.e. local minima of the cost functional to be minimized which can trap the optimization procedure) and are computationally demanding both as far as time and memory resources are concerned. As a matter of fact, these drawbacks hinder their use when imaging of a large (in terms of wavelength) spatial region is required and when the time to achieve imaging is a constraint.
Fortunately, in many practical contexts detecting and localizing the buried scatterers is the primary concern. This can be successfully achieved by adopting imaging algorithms based on simplified scattering models which neglect high order effects due to mutual interactions and which linearize thus the scattering equations. Indeed, the majority of the subsurface imaging algorithms are implicitly or explicitly based in the linear assumption. This is because, in the framework of scatterers localization and shape reconstruction linear imaging algorithms work well far beyond the limits of validity of the scattering linear models upon which they are founded (Pierri et al., 2006).
Here, we will consider only imaging algorithms based on linear models that are the most widely employed in practical scenarios. Such algorithms can be substantially categorized in two groups: migration algorithms and inverse filtering algorithms. Therefore, in the sequel we will describe and compare such algorithms. In particular, starting from the linearized scattering equations the way they are linked and the expected performance are derived under an unified framework.
Before attempting to pursue such a task it must be said that the imaging algorithms are only a part of the signal processing required in subsurface imaging. Most of the effort in signal processing concerns primarily with the reduction of the clutter. In fact, while system noise in GPR can be reduced by averaging, generally GPR data are heavily contaminated by clutter which affects the capability of detecting the buried scatterers. The strongest clutter mainly arises from the air/soil interface and different methods exist to mitigate such a contribution (Daniels, 2004). Connected with this question is the problem of subsurface dielectric permittivity and conductivity estimation. This problem is by its own of significant interest in several applicative fields, ranging from agriculture to monitoring of soil pollution and so on, and for this reason many estimation approaches have been developed. For example the dielectric permittivity can be exploited as an indicator of the soil surface water content (Lambot et al., 2008; Huisman et al., 2003).
However, for imaging purposes it is clear that if the electromagnetic properties of the soil were known one could estimate the field reflected from the air/soil interface and hence can completely cancel the corresponding clutter contribution. What is more, any imaging schemes requires the knowledge of the soil properties to effectively achieve “focusing”; otherwise blurred images where the scatterers appear smeared and delocalized.
Migration procedures essentially aim at reconstructing buried scattering objects from measurements collected above or just at the air/soil interface (from the air side). They were first conceived as a graphical method (Hagendoorn, 1954) based on high frequency assumption. After this approach has found a mathematical background based on the wave equation.
We illustrate the basic idea by assuming that the transmitting and the receiving antennas are just located at the air/soil interface where data are taken accordingly to a multimonostatic configuration (i.e., the receiver is at the same position as the transmitter while the latter moves in order to synthesize the measurement aperture). Moreover, for the sake of simplicity and according to the content of many contributions in the field of migration algorithms, we develop the discussion neglecting the air/soil interface. This means that we consider a homogeneous background scenario with the electromagnetic features of the subsurface. Moreover, the arguments are developed within a two-dimensional and scalar geometry so that the phenomena are described in terms of scalar fields obeying the two-dimensional and scalar wave equation.
Consider a point-like scatterer located in the object space at and denote as the observation variable, that is the positions where the scattered field is recorded. In particular, we assume and , where is the synthesized measurement aperture (see Fig. 2).is the transmitted signal (ideally a delta impulse), the corresponding backscattered field (for the case at hand a B-scan measurement) is given by
where is the soil propagation speed assumed constant and the amplitude factor due to the propagation spreading has been neglected. Accordingly to the considerations said above (see eq. (1)), the backscattered signal in the data space appears as a diffraction hyperbola whose apex is in which can be translated in the image space by exploiting that and (see Fig. 1).
The hyperbolic nature of the B-scan is due to the finite directivity of the antennas: migration aims at compensating for such a spreading by re-focalizing each segment of the hyperbola to its apex. To do this it is observed that the travel-time cannot be directly translated into depth because equal travel-times imply equal distances but the direction of arrival is not specified.
Hence, for each trace (i.e., A scan), the scatterer position should lie on a semicircle centred on the source-receiver position and whose radius is equal to the distance obtained by multiplying the travel-time by half the wave-speed in the soil. Accordingly, each data point is equally distributed over a semicircle in the image space (Hogendoorn, 1954; Bleinstein & Gray, 2001) so that all the semicircle intersect at .
In order to estimate the travel-time, such a method requires that the hyperbolic curve be well discernable from other data features; this entails that a good Signal to Noise Ratio (SNR) is needed and that the scattering scenario is not much complex. In any case, an extended scatterer is considered as an ensemble of point scatterers which makes it clear that linearity of the scattering model is implicitly assumed.
A counterpart of Wave Interference Migration is the so-called Diffraction Summation (Gazdag & Sguazzero, 1984) also knew as pixel-driven approach (Marklein et al., 2002).
In this method the object space is divided in pixels and for each of them a diffraction hyperbola is constructed in the data (image) space. The reconstruction at each pixel is then obtained by summing all the traces (A-scans) that the synthetic hyperbola intersects in the data space.
This procedure can be implemented automatically and requires the evaluation of the following summation integral for each pixel
where and are the measurement aperture and the interval of time during data are collected, respectively, and is the corresponding migrated data.
If we denote as the Fourier transform of (by considering the temporal Fourier kernel ), then previous equation can be recast as
with being the adopted frequency bandwidth. Now, by putting , eq. (4) can be rewritten in terms of integration in the wave-number domain as
where now denotes the frequency band in the domain and an unessential scalar factor has been omitted. Eq. (5) allows us to point out the equivalence between the Diffraction Summation and the Range Migration Technique presented in (Lopez-Sanchez & Fortuny-Guasch, 2000) and also to the Synthetic Aperture Focusing Technique (SAFT) (Markelein et al., 2002).
As can be seen, these methods are essentially based on a convolution in and an integration in . The convolution can be conveniently computed in the spatial Fourier domain as
where is the selected frequency band in the spatial spectral , is an amplitude factor where has been neglected, . is the Fourier transform with respect to of (the spatial Fourier kernel exp(jk x ) is considered) whereas the Fourier transform of the exponential term has been evaluated by means of the stationary phase method. Now, by substituting the frequency wavenumber k by the integration in k z , eq. (6) can be recast as a double Fourier transform as
Eq. (7) has the computational advantage that it can be computed by Fast Fourier Transform (FFT) algorithms directly for each point in the object space but also requires data to be interpolated and re-sampled according to a rectangular grid in the k x -k z spatial spectral domain (Stolt, 1978).
When the amplitude factor is ignored, Eq. (7) becomes identical to the Synthetic Aperture Radar (SAR) imaging algorithm presented in (Soumekh, 1999) and also very similar to the F-K Migration as outlined in (Gilmore et al., 2006).
Interestingly, all the migration algorithms (a part the F-K Migration) have not been obtained by invoking some mathematical theory of the scattering phenomenon. Rather, they have been developed by simply adopting intuitive physically based arguments.
A more rigorous way to derive migration methods can be obtained by invoking on the wave equation (scalar for the case at hand) and by adopting the so-called “exploding source model”. According to this model, the scattered field is thought as being radiated by a source at time t=0 embedded within a medium characterized by a wave velocity being half the actual one (this corresponds to consider the wave number of this equivalent medium twice the one of the subsurface medium). Therefore, outside the source region the scattered field is a solution of the following wave-equation (in the frequency domain)
with and where we still adopted to denote the scattered field.
By Fourier transforming with respect to x one obtains that the field spatial spectrum at quota is related to the one at quota , with , as
where k z is the same as previously defined and only up-travelling waves (that is along the negative z direction) have been considered. Accordingly, the field in the objects space can be determined by forcing the boundary condition over the measurement line at z=0 and by a back-propagation procedure. Finally, a superimposition along the frequency ω returns
This is the so-called F-K Migration which has been shown to be a generalization of the Doppler compression technique typical in SAR imaging which, in turn, can be considered a small angle (of view) approximation version (Cafforio et al., 1991). Moreover, it is not difficult to understand why this migration algorithm is also addressed in the literature as Phase Shift Migration (Gazdag & Sguazzero, 1984).
As can been seen, the Phase Shift Migration is very similar to the Summation Diffraction (in the frequency domain) and can be readily recast as a double Fourier integral in the spatial/spectral domain. However, it has been obtained by requiring less approximations on the underlying model. In particular, the amplitude spreading terms have not been ignored and the stationary phase method has not been employed. Indeed, the Phase Shift Migration derivation can be considered as the mathematical rationale supporting the intuitive discussion under which previous migration schemes have been developed.
Unfortunately, this mathematical model contains an implicit contradiction: while the field is back-propagated as a solution of a homogeneous wave equation the exploding source assumes it being radiated by a localized source. From a mathematical point of view, this contradiction manifests in the assumption of considering the scattered field as being composed of only up-going (towards the air/soil interface) waves. Therefore, it is expected that the field “extrapolation” in-depth in the object space works as long as the scatterers are reached. This becomes particularly apparent when the scatterers are buried within a dispersive and dissipative medium. In this case, while back-propagating the aperture field, if evanescent waves (which reflect the presence of the scatterers) are not filtered out, the migration rapidly degrades for points in the object plane beyond the scatterer (Oden et al., 2007).
By inverting the Fourier transformation with respect to , eq. (10) can be rewritten as
where is the complex conjugated Green’s function. Migration scheme reported in eq. (11) is known as the Rayleigh-Sommerfeld holography which is a particular case (when data are collected over an infinite line) of the so-called Generalized Holography ( Langenberg, 1987) which in turn is founded on the Porter-Bojarski integral equation.
The connection with the Generalized Holography is important because it establishes, in rigorous way, the relationship between the migrated field and the secondary sources, that is the ones arising from the interaction incident field-scatterers, and hence between the migrated field and the scatterers. This is especially true under scattering linear models (Langenberg, 1987). Note that in all the previous migration schemes the relationship between the migrated field and the scatterers has been implicitly assumed as an ansatz.
This question will be further addressed in the next section.
Finally, we observe that the time domain version of eq. (11) is nothing else that the well known Kirchhoff Migration (Berkhout, 1986). In fact, by Fourier transforming eq. (11) with respect to we roughly obtain
where takes into account the angle between the unit normal vector at the measurement line and the vector .
5. Scattering equations and the Born approximation
In the previous section we described different migration algorithms and we showed that they are all very similar since one can pass from the different migration implementations by Fourier (spatial and/or temporal) transform operators.
However, the link between the migrated field and the scatterers to be reconstructed has not been clearly shown and remained only supported by intuitive arguments.
To cope with this question the equations describing accurately the scattering phenomenon are needed. Therefore, in this section, we just introduce the equations governing the scattering phenomena. Furthermore, we also shown as they simplify under the Born approximation for the case of penetrable scatterers (Chew, 1995).
Accordingly, the subsurface imaging problem is cast as an inverse scattering problem where one attempts to infer the electromagnetic properties of the scatterer from the scattered field measured outside the scatterer.
The statement of the problem is then the following: given an incident field, E inc , that is the electromagnetic field existing in the whole space (the background medium) in absence of the scattering object and generated by a known sources, by the interaction with the object the scattered field E s arises; from the knowledge of the scattered field some properties about the scatterer, either geometrical and/or structural, have to be retrieved.
Hence, it is mandatory to introduce the mathematical equations subtending the scattering phenomena to solve the above stated problem. To this end, we refer to a two-dimensional and scalar geometry as done in the previous section.
We consider a cylindrical dielectric object (that is invariant along the axis out-coming from the sheet) enclosed within the domain D illuminated by an incident field linearly polarized along the axis of invariance. The scattered field is observed over the domain (not rectilinear necessarily). Moreover, we denote by e by the equivalent dielectric permittivity function of the unknown object and that of the background medium, respectively. In particular, the latter must not be necessarily constant (i.e., non-homogeneous background medium is allowed) but has to be known. The magnetic permeability is assumed equal the one of the free space everywhere. The geometry of the problem is detailed in Fig. 2.
The problem, thus, amounts to retrieving the dielectric permittivity profile of the unknown object by the knowledge of the scattered field E s .
It can be proven that the frequency domain mathematical relationship between the data and the unknown is furnished by the following scalar Helmoltz equation (Chew, 1995)
with is the total field, is the subsurface (background) wave-number and is the dimensionless contrast function. By adopting the Green’s function method, eq. (13) leads to the pair of scalar integral equations (Colton & Kress, 1992)
where is the pertinent Green’s function, is the observation point and is the source’s position.
In accordance to the volumetric equivalence theorem (Harrington, 1961), the above integral formulation permits to interpret the scattered field as being radiated by secondary sources (the ”polarization currents”) which are just located within the scatterers spatial region.
The reconstruction problem thus consists of inverting the “system of equations (14)” for the contrast function. However, since (from the first equation) the field inside the scatterers depends upon the unknown contrast function, the relationship between the contrast function and the scattered field is nonlinear. The problem can be cast as a linear one if the first line equation is arrested at the first term of its Neumann expansion (Krasnov et al., 1976). In this way, is assumed within the scatterer region and the so-called Born linear model is obtained (Chew, 1995). Accordingly, the scattering model now becomes
We just remark that, within the linear approximation, the internal field does not depend on the dielectric profile which is the same as to say that mutual interactions between different parts of the object and between different objects are neglected. In other words, this means to consider each part of the object as elementary scatterer independent on the presence of the other scatterers. It is worth noting that this is the same assumption as we made a priori while discussing about the migration algorithms.
Now, by considering a homogeneous background medium and a monostatic data collection ( ), as done in the previous section, and by resorting to asymptotic arguments as far as the Green’s function is concerned, eq. (15) can be rewritten as
where a two-dimensional filamentary current with current has assumed as source of the incident field and a scalar factor has been omitted. If is meant as in the previous section
Then by Fourier transforming the scattered field with respect to and by exploiting the plane-wave spectrum of the Green’s function we obtain
with is defined as in the previous section. Accordingly, when the spatial Fourier transformed data are rearranged in the spectral plane, the contrast function can be obtained as
It is important to note that eq. (19) coincides just to the F-K Migration when integration in is replaced by the integration in . Therefore, within the framework of the linear Born approximation we have established a clear connection between the migrated field and scatterers in terms of their contrast functions.
6. Inverse filtering imaging algorithms
As stated in the previous section, under the Born approximation and for a two-dimensional and scalar geometry, the scattering phenomenon may be modelled through a linear scalar operator
where and are the contrast function (the unknown of the problem) and the scattered field (the data of the problem), respectively.
Moreover, and represent the functional sets within we search for the contrast function and the one we assume the scattered field data belong to, respectively. In particular, we assume them as Hilbert spaces of square integrable functions, the first one of complex valued functions defined on the investigation domain , whereas the second one of such functions supported over . However, in general the data space depends on the adopted configuration, that is on the choice of the measurement domain and the strategy of illumination and observation.
It is worth remarking that the choice of and as Hilbert spaces of square integrable functions accommodates the fact that no a priori information are available on the unknown except that one on the finiteness of its ”energy” as dictated by physical consideration. On the other side, it assures that is ”broad” enough to include the effect of uncertainties and noise on data.
Thus, the problem amounts to inverting equation (20) to determine the contrast function. Since the kernel of the operator in (20) is a continuous function on , then the linear operator is compact (Taylor and Lay, 1980). As stated above this means that the inverse problem is an-illposed problem (Bertero, 1989). For compact and non-symmetric operator (as the one at hand) the singular value decomposition is a powerful tool to analyze and to solve the problem.
Therefore, we denote as the singular system of operator . In particular, is the sequence of the singular values ordered in non-increasing way, and are orthonormal set of functions solution of the following shifted eigenvalue problems
where is the adjoint operator of [Taylor & Lay, 1980], which span the orthogonal complement of the kernel of , , and the closure of the range of , , respectively.
where denotes the scalar product in the data space .
By virtue of the compactness of , is not a closed set (Krasnov et al., 1987). This implies that the Picard’s conditions is not fulfilled for any data functions (i.e., the ones having non null component orthogonal to ) ; thus the solution may not exist and does not depend continuously on data (Bertero, 1989). This is just a mathematical re-statement of ill-posedness. From another point of view, we have to take into account that the actual data of the problem are corrupted by uncertainties and noise , hence
Now, because of the compactness, the singular values tend to zero as their index increases.
This implies that, the second series term (noise-related) in (23) in general does not converge. This leads to an unstable solution since even small error on data are amplified by the singular values close to zero.
The lack of existence and stability of solution can be remedied by regularizing the addressed ill-posed inverse problem (De Mol, 1992). For example, this can be achieved by discarding, in the inversion procedure, the ”projections” of data on the singular functions corresponding to the ”less significant” singular values, that is by filtering out the singular functions corresponding to the singular values which are below to a prescribed noise dependent threshold. This regularizating scheme is known as Numerical Filtering or Truncated Singular Value Decomposition (TSVD) and is the simplest one within the large class of windowing based regularizating algorithms (Bertero, 1989).
Accordingly, the finite-dimensional approximate stable solution is then given by
More in general, the basic idea of regularization theory is to replace an ill-posed problem by a parameter dependent family of well-posed neighbouring problems
with being the so-called regularization parameter (in the TSVD this corresponds to the truncation index in eq. (24)) and the noise level, so that to establish a compromise between accuracy and stability (Groetsch, 1993). In particular, as also and the regularized reconstruction must tends to the generalized inverse whose outcome is just shown in eq. (22).
The Tikhonov regularization scheme is a widespread regularization scheme which takes advantages from exploiting a priori information about the unknown (Tikhonov, 1977). In this case the inversion problem is cast as a constrained optimization problem
Here the constraint arises from the available a priori information and in general can be different from the energy constraint reported in the equation above.
The Landweber regularization scheme recasts the first kind integral equation to be inverted as an integral equation of second kind so that a well-posed problem is obtained (Groetsch, 1993). Accordingly, eq. (26) is recast as
and a solution is then obtained by an iterative procedure. In this case the regularization parameter is the number of iterations exploited in the iterative minimization .
These regularization scheme can be compared if all are analyzed in terms of the operator properties. This can be done by expressing the different regularized reconstruction in terms of the singular system. By doing so, one obtains
for the Tikhonov regularization, and
for the Landweber method.
As can be seen, all this regularization methods result (but in a different) filtering of the unknown spectral expansion.
It is clear that, a part the computational convenience which can dictate the regularization algorithm to adopt, the key question is the choice of the regularization parameter. As remarked, this choice has to be done by accounting for the noise level, the mathematical features of the operator to be inverted and the available a priori information about the unknown. Different methods exist to select the regularization parameter. Such methods can explicitly exploit the knowledge of the noise level, (such as the Morozov discrepancy principle) or not (such as the generalized cross validation) (Hansen at al., 2008).
In inverse scattering problem, the scattered field has finite number of degrees of freedom (Solimene & Pierri, 2007) which in turn depends on the parameter of the configuration. Therefore, such an information can be exploited to select the regularization parameter.
The singular system formalism can be also employed to compare migration and inverse filtering. By looking at the scattering eq. (17), it is seen that the Diffraction Summation migrated field reconstruction, substantially, corresponds to achieve inversion by means of the adjoint operator, that is
whose expression in terms of the singular system is
It is readily noted that migration allows to obtain a stable reconstruction because the singular values now appear at the numerator. However, it is not a regularization scheme in the sense of Tikhonov (Tikhonov & Arsenin, 1977) because even in absence of noise the actual unknown is not retrieved. From a practical point of view, this entails an intrinsic limit on the resolution achievable in the reconstructions regardless of the noise level. To explain this we observe that and its reconstructed version are linked by the following integral relationship
This entails that, for each of the inverse filtering procedure described above, the model resolution kernel tends to a Dirac delta, that is, , when noise is absent and (Dudley et al., 2000). This is not true for migration algorithms. Therefore, this means that, depending on the noise level, inverse filtering can achieve better resolution in the reconstructions.
This is shown in Fig. 3 where TSVD and F-K migration performances are compared.
However, as long as lossless medium is considered, as migration algorithms can be easily recast in the spatial Fourier domain they generally have a lower computational cost than inverse filtering. Indeed, in lossy scenarios, losses should be considered in the inversion procedure and this impairs the possibility of adopting only FFT based imaging algorithms (Meincke, 2001).
7. Reconstruction results
In order to give an idea of how subsurface reconstructions look like, in this Section we report a brief survey of reconstruction performances against synthetic data. The scatterd field data have been computed thanks to the Finite Difference Time Domain (FDTD) GPRMAX code (Giannopoulos, 2005) in time domain; these time-domain data are then Fourier transformed to pass to the frequency domain.
The TSVD inverse filtering approach based on the Born model is adopted to achieve the reconstructions; in particular, to emphasize the usefulness of linear reconstruction schemes test cases concerning applicative cases where the objects have electromagnetic properties very different from the ones of the soil (test cases outside of the Born regime) are considered.
The inversion algorithm assumes an investigation domain D whose extent ranges from 0.1m to 1.5m in depth and with the semi-horizontal extent a=0.75m. The soil has a relative dielectric permittivity equal to 9 and an electrical conductivity equal to 0.01 S/m. The measurement domain is located at the air/soil interface and has extent 1.5m. Over such a domain, 51 data are taken at a spatial step of 3cm. The working frequency Band Ω ranges from 300 to 900 MHz.
Different cases of a single object have been considered; in particular, we have considered the same geometry of the object, i.e., a circular cylindrical object of radius of 0.1 m and whose center is located at the depth of 0.6m. Three different values of the relative dielectric permittivity filling the object have been considered, i.e., 9.1 (low Born model error), 4 and 16.
The reconstruction results are given in terms of the modulus of the retrieved contrast function normalised with respect to its maximum in the investigation domain. The regions where the normalised modulus of the contrast function is significantly different from zero give information about the support of the targets.
Figure 4 depicts the reconstruction for the case of the target with relative dielectric permittivity equal to 9.1; the reconstruction permits to localize and reconstruct well the shape of the upper and lower sides of the circle. The features of the reconstruction agree with the theoretical expectations reported in (Persico et al. 2005) where the filtering effect of the regularized inversion algorithm consists in a low-pass filtering along the antenna movement direction (x-axis) and a band-pass filtering along the depth (z-axis). In this case, when we exploit the TSVD scheme for the inversion we retained in the summation (24) the singular values larger than 0.1 the maximum one, that is the TSVD threshold is equal to 0.1.
The second case is concerned with the circular target having a relative dielectric permittivity equal to 4 so that the Born model error is significant. Figures 5 and 6 depict the corresponding reconstructions when the TSVD threshold is equal to 0.1 and 0.01, respectively.
As can be seen, both the results permit to point out in a correct way the location and the shape of the upper side of the target. Moreover, by adopting a smaller value of the TSVD threshold (more singular values in the TSVD summation (24)), the achievable resolution
improves and both the upper and the lower sides of the scatterers are discernable in the reconstruction. In particular, for the lower side of the target, we observe that the reconstruction approach is able to detect it but does not provide to its actual position. This can be explained by noting that the velocity of the electromagnetic wave in the target is equal to v=15 cm/ns and different from the one assumed in the model equal to v=10 cm/ns. Therefore, the inversion model assumes a velocity smaller than the actual one inside the target and this leads to a reconstruction of the lower side at a depth smaller than the actual one.
Finally, we consider the case with the circular target having a relative dielectric permittivity equal to 16 (high Born model error). Figures 7 and 8 depict the tomographic reconstruction results when the TSVD threshold is equal to 0.1 and 0.01, respectively.
We can observe that both the results permit to point out in a correct way the location and the shape of the upper side of the target; also, by adopting a smaller value of the TSVD threshold (more singular values in the TSVD summation (24)), the reconstruction is slightly improved.
For the lower side of the target, we observe that the reconstruction approach is able to detect it but the reconstruction does not correspond to the actual target’s shape. Similarly to the explanation given above, we note that the velocity of the electromagnetic wave in the target is now equal to v=7.5 cm/ns and different from the one assumed in the model equal to v=10 cm/ns. Therefore, the inversion model assumes a velocity larger than the actual one inside the target and this arises in a reconstruction more deeply located. Also, it is worth noting that for the smaller TSVD threshold, the reconstruction of the lower side has a shape more similar to the one of the target.
Finally, we present the case (see fig. 9) of a wrong choice of the TSVD threshold equal to 0.0001, in this case the necessity to improve the reconstruction (by lowering the TSVD threshold) is made it not possible by the effect of the uncertainties on data (in this particular case related only to numerical errors) that produce a reconstruction very difficult to be interpreted.
This chapter has dealt with the imaging of the buried targets thanks to Ground Penetrating Radar. After a brief description of the operating principles of GPR we have focussed the attention on the data processing that is one of the main questions in order to achieve interpretable images of the subsurface.
In particular, we have described the main versions of the migration technique, usually exploited in data processing, and a relation between them has been pointed out. Then, we have presented a more recent data processing approach based on microwave tomography. Such a tomographic approach faces the inverse scattering problem and the adoption of a physics-based model allows us not only to obtain more interpretable reconstruction but also to analyse the performances of the reconstruction approach.
In particular, the performances of a microwave tomography based approach exploiting a simplified model of the electromagnetic scattering was investigated in detailed way. Such a linear inverse approach has revealed suitable in realistic applications where the aim is the geometry of the targets. Finally, the feasibility and limitations of this linear approach have been outlined by an analysis against synthetic data.