A Controlled-Site Comparison of Microwave Tomography and Time-Reversal Imaging Techniques for GPR Surveys

Abstract

This paper provides a comparative study between microwave tomography and synthetic time-reversal imaging techniques as applied to ground penetrating radar (GPR) surveys. The comparison is carried out by processing experimental data collected at a controlled test site, with different types of buried targets at given subsurface depths and representative soil conditions. It is shown that the two techniques allow us to obtain complementary information about position, depth and size of the targets from a single GPR survey.

1. Introduction

Ground penetrating radar (GPR) finds a large number of applications related to detection and imaging of subsurface targets and anomalies, such as underground utilities, pipes, chemical spills, groundwater levels, etc. Historically, GPR data has been typically analyzed and interpreted based on a “visual” analysis of the radargram [1,2]; however, this analysis is able to provide reliable interpretation only in simple scenarios. In order to improve the interpretability of GPR data, inverse scattering and migration algorithms [3,4,5,6,7,8], including microwave tomography (MT) techniques [9,10,11], have also been used in different scenarios related to GPR surveys. MT has been successfully used, for example, to investigate water leakage from a large metallic pipe in [12], where it was able to give more reliable information than a direct radargram analysis. In addition, several GPR surveys at different time intervals were done for a controlled oil spill and analysed through MT in [13].

In all these scenarios, environmental clutter can affect the measured data and deteriorate imaging results, making it more difficult to estimate the target geometry [14]. MT was compared to conventional migration techniques for imaging metallic and plastic pipes in [9], under controlled experimental conditions, where it was shown that MT is capable of providing improved accuracy and resolution. A similar study was carried out in [15] in a forensic application scenario. A theoretical analysis of the reconstruction performance of the MT and migration was carried out in [16], based on the Singular Value Decomposition of the relevant operator. MT also proved capable of providing good results for multi-frequency systems [17]. The effect of background electrical conductivity values on tomographic images was addressed [18]. This latter study showed that the more accurate the estimated background conductivity value, the higher the contrast function image sharpness associated with the anomaly of interest, suggesting that this effect can be alternatively exploited to estimate the background medium conductivity as well. Structural degradation assessment, using a combination of MT and seismic tomography, was discussed by [19], for cultural heritage monitoring.

Finally, the flexibility enabled by MT has permitted its use for processing data acquired from airborne surveys [20], where it was shown that very good results can be obtained, especially in areas with no vegetation covers, such as glaciers.

Time-reversal-based (TR) techniques were first developed for acoustics [21,22], and later successfully used in several applications including medicine [23,24], non-destructive testing and evaluation [25], atmospherics studies [26], microwave remote sensing [27], as well as near-subsurface geophysics [28,29,30,31,32,33]. As the name suggests, TR-based techniques exploit the invariance of the wave equation under time reversal. The received data can be synthetically time-reversed (in a first-in last-out fashion) and, after physical or synthetic re-transmission to the region of interest, used to create wavefields that automatically focus on reflective targets and/or anomalies. Under certain conditions, TR techniques can be used for detection and localization of obscured targets in cluttered and rich-scattering environments. Variants on the basic TR algorithm exist, which allow for selective focusing (on secondary or weaker targets) [34,35] or tracking of obscured moving targets as well [36].

The objective of this work is to compare MT and TR by analyzing reconstruction results obtained by experimental data from GPR surveys on a controlled site with known subsurface targets. This comparison has the objective of examining the relative weaknesses and strengths of each method under conditions pertaining to realistic GPR field acquisitions. The consideration of a controlled site with known targets allows for a better assessment of the capabilities and drawbacks of the two approaches. It is worth noting that the test site is at the scale of the realistic conditions and is more challenging than usual laboratory conditions. The targets under investigation are representative examples of objects normally found in archaeological sites (disturbed soil and ceramic vases), geotechnical evaluations (concrete tubes) and environmental surveys (storage tanks). The comparison between MT and TR is presented here for the first time in literature and suggests that the two approaches could be potentially used in a complementary fashion.

The paper is organized as follows: Section 2 and Section 3 provide a brief overview of the MT and TR methodologies, respectively. Section 4 describes the controlled site characteristics and the set of buried targets considered. Section 5 compares the results from experimental GPR data. Section 6 presents a discussion of the results and, finally, some concluding remarks are provided in Section 7.

2. Microwave Tomography

Since both MT and TR have been widely studied in the past, we will only discuss them very briefly here.

MT formulates GPR data processing as an inverse scattering problem [37,38]. Consider, for simplicity, the 2D geometry depicted in Figure 1. Each triangle indicates the position of the transmitter-receiver pair. In a common-offset configuration, each trace of the radargram is acquired at a different point, $\overrightarrow{r}\left(x,z\right)$, at the air-soil interface. The profile is coincident with this interface, as shown in Figure 1. The electric field irradiated by the transmitter antenna in the absence of targets corresponds to the incident field, ${\overrightarrow{E}}_i$. The targets are assumed to be buried inside the region of interest, $I$, which is discretized by a regular grid of points, ${\overrightarrow{r}}^{\prime }\left(x^{\prime },z^{\prime }\right)$. The interaction of the incident field with a buried target within the incident field generates the scattered electric field, ${\overrightarrow{E}}_s$. The summation of scattered and incident fields results in the total electric field, $\overrightarrow{E}$. The scattered field, at the measurement domain, $O$, conveys information about the buried targets and represents the input data to the processing scheme. According to [39] the scattered field can be written as:

$${\overrightarrow{E}}_s\left(\overrightarrow{r},\omega \right)=k_b^2{\int }_I\overline{\overline{G}}\left(\overrightarrow{r}|{\overrightarrow{r}}^{\prime }\right)·\overrightarrow{E}\left({\overrightarrow{r}}^{\prime }\right)\chi \left({\overrightarrow{r}}^{\prime }\right)d{\overrightarrow{r}}^{\prime };\hspace{1em}\overrightarrow{r}\in O$$

(1)

where $k_b=\sqrt{{\omega }^2{\mu }_0{\epsilon }_b-i\omega {\mu }_0{\sigma }_b}$ is the complex wavenumber of the background medium, $\omega$ is the angular frequency, ${\mu }_0$ is the free-space magnetic permeability, ${\epsilon }_b$ is the dielectric permittivity, ${\sigma }_b$ is the electric conductivity and $\overline{\overline{G}}$ is the background Green’s function [39]. The contrast function, $\chi$, is the unknown of the tomographic imaging inverse problem and is defined as:

$$\chi \left({\overrightarrow{r}}^{\prime }\right)=\frac{\epsilon \left({\overrightarrow{r}}^{\prime }\right)}{{\epsilon }_b}-1$$

(2)

where $\epsilon \left({\overrightarrow{r}}^{\prime }\right)={\epsilon }_0{\epsilon }_r\left({\overrightarrow{r}}^{\prime }\right)-i\left[\sigma \left({\overline{r}}^{\prime }\right)/\omega \right]$ is the complex dielectric permittivity inside the region of interest, $I$, with ${\epsilon }_r\left({\overrightarrow{r}}^{\prime }\right)$ corresponding to its dielectric constant and $\sigma \left({\overline{r}}^{\prime }\right)$ corresponding to its electrical conductivity, and ${\epsilon }_b={\epsilon }_0{\epsilon }_{rb}-i{\sigma }_b/\omega$ is the complex dielectric permittivity of the homogenous background medium. The problem in (1) is nonlinear with respect to the electrical field. The most common linearization strategy for this problem is to employ the Born approximation [39]. Under its assumptions, the total field is approximated by the incident field in (1), so that the expression for the scattered electric field becomes:

$${\overrightarrow{E}}_s\left(\overrightarrow{r},\omega \right)\approx k_b^2{\int }_I\overline{\overline{G}}\left(\overrightarrow{r}|{\overrightarrow{r}}^{\prime }\right)·{\overrightarrow{E}}_i\left({\overrightarrow{r}}^{\prime }\right)\chi \left({\overrightarrow{r}}^{\prime }\right)d{\overrightarrow{r}}^{\prime };\hspace{1em}\overrightarrow{r}\in O$$

(3)

The problem stated by (3) of estimating the contrast function, $\chi$, from the scattered field, E_S, is generally ill-posed, in particular because the stability of the solution is affected by the noise present in the data. In order to achieve a stable solution, regularization schemes can use used. Here, we adopt the Truncated Singular Value Decomposition (TSVD) [40]. The Equation (3) can be written as ${\overrightarrow{E}}_s=ℒ\chi$, where $ℒ$ represents the linear operator connecting the contrast function to the scattered field data. In the presented approach, $ℒ$ is discretized using the method of moments [41]. By applying an SVD on this linear operator, we can write:

$${\overrightarrow{E}}_s={{\displaystyle \sum }}_{n=0}^{\infty }{s}_n\langle \chi ,u_n\rangle v_n$$

(4)

from which the solution for $\chi$ can be obtained by inverting (4), i.e.,

$$\chi ={\displaystyle \sum }_{n=0}^{N_T}\frac{\langle {\overrightarrow{E}}_s, v_n\rangle }{s_n}{u}_n$$

(5)

where ${\left\{v_n\right\}}_{n=0}^{\infty }$ is the set of singular vectors (orthonormal basis) in the data space, ${\left\{u_n\right\}}_{n=0}^{\infty }$ is the set of singular function vectors (orthonormal basis) in the space of unknowns. ${\left\{s_n\right\}}_{n=0}^{\infty }$ is the set of singular values, ordered in a decreasing order, and $N_T$ is the truncation index. $N_T$ is a problem-dependent parameter, selected to achieve a good balance between resolution and stability against noise [37,38]. The truncation indexes, $N_T$, used here, are discussed in Section 5.

3. Time-Reversal-Based Technique

TR-based techniques were first introduced for ultrasonic waves [21,22] and later extended to electromagnetic waves. They explore the invariance of the wave equation under time reversal [42]. This invariance is only exact in reciprocal and lossless media but, under certain conditions, the techniques can be also applied to lossy media as well. In rich scattering scenarios, TR can achieve super-resolution and for wideband signals, it can provide statistical stability for imaging in random media [43,44,45,46]. The basic TR process is given by the following steps [42]:

(i)

A short pulse is transmitted from one “or more” transceivers to the region of interest, where it is scattered by one or more targets;

(ii)

The scattered signal is registered by the transceivers;

(iii)

The received signal waveform is time-reversed (first-in, last-out) and retransmitted (either physically or synthetically) to the region of interest;

(iv)

Due to time invariance of the wave equation, the retransmitted waveform will tend to focus around the original target location(s).

In the TR process, the focusing of the retransmitted waveform suffers from limitations in the resolution, because the receivers typically comprise a limited-aspect aperture in practice (i.e., they do not capture the scattered field in all directions). In addition, the evanescent spectrum (present only in the very near-field) is not captured, nor is the presence of losses, if any, in the region of interest.

Nevertheless, under some conditions [42], the resolution enabled by TR can go beyond the one dictated by the conventional diffraction limit. The TR process can also be understood as a matched filter operation in both time and space [42,43]. Assuming that a transmitter at a location, ${\overrightarrow{r}}_0$, sends a pulse, s(t), the signal measured at a receiver located at ${\overrightarrow{r}}_i$ can be expressed in terms of the convolution for $i=1, \dots , N$ being N the number of receivers (for simplicity, we neglect the antenna response in this discussion, which can be separately compensated):

$$f_i\left(t\right)=s\left(t\right)\times h_{{\overrightarrow{r}}_0{\overrightarrow{r}}_i}\left(t\right)$$

(6)

where $h_{{\overrightarrow{r}}_0{\overrightarrow{r}}_i}\left(t\right)$ is the impulse response (time-domain Green’s function) between ${\overrightarrow{r}}_0$ and ${\overrightarrow{r}}_i$. From the reciprocity theorem, we can write $h_{{\overrightarrow{r}}_0{\overrightarrow{r}}_i}\left(t\right)=h_{{\overrightarrow{r}}_i{\overrightarrow{r}}_0}\left(t\right)$, and the time-reversed retransmitted signal at ${\overrightarrow{r}}_0$ due to a source at ${\overrightarrow{r}}_i$ is given by:

$$p_i\left(r_0,t\right)=s\left(-t\right)\times h_{{\overrightarrow{r}}_0{\overrightarrow{r}}_i}\left(-t\right)\times h_{{\overrightarrow{r}}_i{\overrightarrow{r}}_0}\left(t\right)$$

(7)

where $s\left(-t\right)\times h_{{\overrightarrow{r}}_0{\overrightarrow{r}}_i}\left(-t\right)$ is $f_i\left(-t\right)$ Generally, for a time reversal array (TRA) with N transceivers, the received signal is given by the following equation:

$$p(\overrightarrow{r_0},t)={\displaystyle \sum }_i^N{p}_i(\overrightarrow{r_0},t)={\displaystyle \sum }_{i=1}^{N}s\left(-t\right)\times h_{{\overrightarrow{r}}_0{\overrightarrow{r}}_i}\left(-t\right)\times h_{{\overrightarrow{r}}_i{\overrightarrow{r}}_0}\left(t\right)$$

(8)

A more extensive discussion on TR-based techniques and their variants can be found in [42].

4. Test Site

The comparative study comprised eight different targets, representative of objects found in archaeological, geotechnical, and environmental studies. The experimental data was collected at a controlled test site (at the scale of the realistic situations) situated in the Institute of Astronomy, Geophysics and Atmospheric Science at the University of Sao Paulo (IAG/USP), Brazil, see Figure 2, during dry weather. The test site is situated at the border of a sedimentary basin in the southwest of Brazil, characterized by clay soil and clay-sand sediments overlapped to a granite-gneissic basement. The test site is outdoors and is affected by climatic events.

A commercial 200 MHz GPR system, manufactured by GSSI (Geophysical Survey Systems, Inc., Nashua, NH, USA), was used to collect the data. We used 512 samples in each A-scan, 100 ns for the time-window and a total of 50 A-scans/m of sampling. The GPR system employs shielded bow-tie antennas and has a nominal frequency range from about 50 MHz to about 325 MHz.

Table 1. Targets description.

	Targets	Description
1	Disturbed soil	Disturbed soil, with 1 m3 volume
2	Ceramic vase	Empty ceramic vase, with 0.55 m diameter and 1.0 m depth
3	Concrete tube	Horizontal concrete tube (with iron structure), with 0.7 m diam. and 1.0 m depth
4	Concrete tube	Vertical concrete tube (with iron structure), with 0.7 m diam. and 1.0 m depth
5	Concrete tube	Horizontal concrete tube, with 0.26 m diameter and 0.5 m depth
6	Metallic tank	Horizontal metallic tank with 0.59 m diameter and 0.5 m depth
7	Metallic tank	Double horizontal metallic tanks with 0.59 m diameter and 1.0 m depth
8	Metallic tank	Vertical metallic tank with 0.86 m height and 1.0 m depth

provides the list of the targets present underneath each of the GPR tracks considered. Figure 3 shows a schematic view of the target distribution in the subsoil. Before data processing with either MT or TR, the acquired radargrams were pre-processed using ReflexW^TM software [47], using a conventional sequence based on header gain removal, zero-time correction, background removal, gain function, and frequency filtering. Since the targets are buried at different locations in a wide area, the background medium (soil) may exhibit some variation in its permittivity, due to soil content variations in shallow geologic material. Because of this, the (mean) permittivity was first retrieved using the relationship, ${\mathit{\epsilon }}_{\mathit{r}}={\left(\mathit{c}/\mathit{v}\right)}^{\mathbf{2}}$, where $\mathit{c}$ is the speed of light and $\mathit{v}$ is the phase velocity in the subsoil. The depths of targets are already known, and using the two-way GPR time travel, we can determine the velocity, v. The (mean) conductivity was also retrieved a posteriori, by analyzing the images, after performing the inversion with different conductivity values and selecting the value that gave the better focused tomographic image. This procedure is similar to that described by [18]. The values adopted for the electrical properties in subsoil and the target size, in terms of the wavelength propagation in the corresponding medium, in order to show the expected resolution, are summarized in

Table 2. Electromagnetic properties of the background media surrounding each target, and size of the targets compared to the propagating wavelengths, (λ₁₈ = 0.35 m, λ₁₁ = 0.45 m). φ, h and s correspond to diameter, height and separation, respectively, of the objects.

Target	${\mathit{\epsilon }}_{\mathit{r}\mathit{b}}$	${\mathit{\sigma }}_{\mathit{b}}$ [S/m]	Size (m)
disturbed soil	18.0	0.007	2.85λ18
ceramic vase	18.0	0.007	1.57λ18 (φ)
horizontal concrete tube	11.1	0.001	1.55λ11 (φ)
vertical concrete tube	11.1	0.001	1.55λ11 (φ) 2.22λ11 (h)
horizontal concrete tube	11.1	0.001	0.57λ11 (φ)
horizontal metallic storage tank	18.0	0.007	1.68λ18 (φ)
pair of horizontal metallic storage tanks	18.0	0.007	1.68λ18 (φ) 2.86λ18 (s)
vertical metallic storage tank	18.0	0.007	1.68λ18 (φ) 2.46λ18 (h)

.

5. Comparative Results

In order to apply MT, we considered 19 frequencies, equally distributed in the 50 to 325 MHz nominal bandwidth of the antenna. The investigation domain ($I$) was discretized in 0.025 m × 0.025 m pixels for all of the considered cases. The TSVD regularization parameter was chosen, case by case, according to the best reconstruction; the values of the regularization parameters are listed in

Table 3. Threshold values adopted for the Truncated Singular Value Decomposition (TSVD) regularization.

Target	TSVD Threshold Value [dB]	TSVD Threshold Index
disturbed soil	−27.5	728
ceramic vase	−30.7	782
horizontal concrete tube	−25.6	630
vertical concrete tube	−25.1	626
horizontal concrete tube	−16.3	169
horizontal metallic storage tank	−31.5	1010
pair of horizontal metallic storage tanks	−42.3	1150
vertical metallic storage tank	−24.8	643

.

TR images can be obtained according to the procedure described in [32,33], where the standard deviation of the back-propagated TR wavefield amplitudes, sampled in either time or space (i.e., depth or co-range), is used to provide the dispersion of each data set. In this manner, larger variations in the TR amplitudes caused by focusing effects are emphasized without the need for a precise a priori estimate of the focusing instant (as it turns out to be necessary in other TR approaches). The standard deviation is computed for the whole data set, comprising a three-dimensional TR matrix. After data acquisition, the TR retransmission step (backpropagation) is carried out numerically. The FDTD (finite-difference time-domain) algorithm [48] is employed here for this purpose.

In [32,33] two types of data processing are applied for computation of the standard deviation, denoted as Mode 1 and Mode 2. Mode 1 aims at emphasizing the cross-range resolution or the target position along the GPR track, whereas Mode 2 tends to emphasize the co-range resolution or the target depth. Although these modes provide data that can, in principle, be applied and interpreted separately, the single mode results can be affected by spurious artifacts, such as ringing effects, that can confound interpretation. In order to decrease the image artifacts and improve the results, we exploit the best features of both modes in the present study, by cross-correlating the two resulting images (ccTR). The basic sequence of steps is shown in the diagram in Figure 4. In all of the presented results, the first image (a) shows the conventional radargram, the second image (b) shows the image obtained with MT, and the third image (c) shows the ccTR image. White dashed lines indicate the geometric shape and actual location of each target. In all results, each technique yields “artifacts” in the reconstructed images, in addition to the target response. These anomalies (secondary reflections) are prevalent in most field data, due to the presence of geologic stratification and other inhomogeneities in the subsurface. The indication of these artifacts in all three techniques is provided by red arrows.

Target 1: Disturbed soil. Results from the processing of the GPR data, acquired over the region comprising disturbed soil, are shown in Figure 5. Figure 5a shows the pre-processed radargram, where one can appreciate a reduction in the backscattered field amplitude for the region highlighted by the dashed white rectangle indicating the region of disturbance. The reflection amplitudes are reduced because any significant reflection to geological stratification was removed by the excavation and refilling process. In particular, a strong reflector, related to a shallow geologic interface, is present at a depth of about 0.5 m, beyond the 2.5 m position along the GPR line (indicated by red arrows). This reflector appears at lower depths at the near end of the radargram and is interrupted by the soil disturbance. Figure 5b shows the MT results. The strong responses in the region between the surface and a depth of about 1 m are again associated with reflections caused by geologic stratification on the test site. Although the discontinuity in the response does not reproduce the exact disturbance width, the evidence of soil disturbance is also clear in this image, with a sharper delineation, compared to the radargram result. Figure 5c shows the image obtained based on cross-correlation of the two TR modes. Again, the vertical interface is easy to identify, and the excavation boundaries are well-defined at 1.5 and 2.5 m. Notice that the bottom and the top of the disturbance region cannot be precisely determined in any of the cases.

Target 2: Ceramic vase. Figure 6 shows the results from GPR data for an investigation domain containing a ceramic vase buried at a depth of 1.0 m. A shallow geologic interface at a depth of about 0.5 m can be again discerned from the radargram shown in Figure 6a (indicated by red arrows). There is also a clear reduction in the reflection amplitude above the target location, which results from the soil excavation from burying the target. The anomaly related to the target appears as a hyperbolic feature at a depth of about 1.0 m. However, this hyperbolic feature is merged to another strong horizontal reflection at a depth of around 1.2 m, which is related to geological features. The MT results in Figure 6b show a clearer separation between the target and the geologic background. The discontinuity of the geologic reflectors due to excavation is also evident, as indicated by the yellow dashed lines. The delocalization in depth for the MT reconstruction can be explained as follows: The ceramic vase is a low-loss dielectric object; therefore, the electric field is able to penetrate into the object and this entails the possibility that the tomographic approach is able to reconstruct both the upper and lower edges of the target (see, for example, the reconstruction in [49]). These two reconstructed edges are well detected and visible in the tomographic image if the resolution along the depth is adequate. When the resolution in depth is not adequate (as in the case at hand), the two spots accounting for the edge combine and give the reconstructed spot at the center of the target.

The image obtained based on the ccTR is shown in Figure 6c. The excavation marks are defined and the presence of geologic layers between the ground and a depth of 1.5 m is visible as high correlation values. In the target location, we can see an anomaly with a better lateral delineation than in the conventional radargram result. In particular, the TR results can provide good estimates for depth, position along track, and lateral extent of the target, in this case.

Target 3: Horizontal concrete tube with 0.7 m diameter. Figure 7 shows the results from the GPR survey over a horizontal concrete tube with a diameter of 0.7 m. The radargram presented in Figure 7a exhibits a clear hyperbolic anomaly related to the target at a depth of about 0.8 m. A strong geologic reflector is also visible at ~0.6 m depth, and has been interrupted by the excavation to bury the pipe. The image retrieved from the MT is shown in Figure 7b, where a strong anomaly related to the target is visible close to the top of the true target position, as indicated by the dashed line circle, with a slight vertical deviation with respect to the target contour. A discontinuity is seen in the anomalies related to the geologic reflector, as indicated by the vertical dashed lines. The corresponding ccTR image is shown in Figure 7c, where the anomaly exhibits a similar pattern to the previous cases but with a slight deviation in depth, and the vertical smearing (right red arrow) of the geologic layer cannot provide a good contrast from the target. This may have been caused by the presence of the small anomalies seen above the target. Again, the lateral boundaries of the excavation region can be well discerned.

Target 4: Vertical concrete tube with 0.7 m diameter. Results corresponding to the vertical concrete tube are presented in Figure 8. The response associated with this target is strong, but mostly confined to the top boundary of the target (shallow end), due to the strong reflection and the fact that the electric field is not able to penetrate into the metallic structure. The target response in the MT image, presented in Figure 8b, arises from the top boundary (upper side) as well. The TR processing result is shown in Figure 8c, which yields a good match for both the target position, as well as vertical and horizontal dimensions. The top of the anomaly is rightly located, and the bottom is clearly defined at 3.0 m approximately. However, together with the target vertical extent recovery, secondary artifacts arise in other portions of the image, associated with deeper geological features.

Target 5: Horizontal concrete tube with 0.26 m diameter. Figure 9 shows the results for the GPR survey over a horizontal concrete tube with a diameter of 0.26 m. A clear hyperbolic anomaly is present at the depth of 0.5 m in the radargram presented in Figure 9a. The corresponding tomographic image in Figure 9b shows a strong anomaly, coincident to the true target location. Weaker anomalies related to the geologic background are indicted by the red arrows and exhibit a spatial arrangement roughly delineating the lateral and lower edges of the excavation done for target installation in the test site (dashed lines). Despite a visible anomaly associated with the concrete tube, the ccTR result in Figure 9c does not show a good resolution for this particular target, because it is not distinguishable from the strong anomaly positioned at ~1.5 m. There is an anomaly in the same position of the target, but, despite its barely-seen circular shape, it is weaker than the anomaly associated with the geologic reflector. In this case, as one can see in all three results, there is a strong geological reflector at a location near 1.5 m along the track at about 0.8 m in depth, which can be confused with another target.

Target 6: Single horizontal metallic storage tank. Figure 10 depicts the results from the GPR survey over a horizontal metallic storage tank. A clear, well-defined hyperbolic anomaly is observed in the radargram, shown in Figure 10a, at about 0.5 m in depth and 2.0 m along the track. A horizontal reflection, related to a geologic interface, is observed as well near 1.0 m deep (indicated by red arrows). The image retrieved from MT in Figure 10b shows a single strong anomaly, coincident with the upper edge of the target position. Anomalies caused by the geologic reflector are visible at 1.0 m in depth, but are considerably weaker compared to the target. This is because of the high contrast between the target and the background medium in this case. Figure 10c shows the TR result, where the response from the geologic layers is clearly visible again. Although the estimate target location along the track coincides well with the actual one, this target generates an image anomaly larger in size than the true target dimensions, and it is not possible to accurately determine the actual depth of the storage tank.

Target 7: Pair of horizontal metallic storage tanks. Figure 11 shows the images related to the pair of metallic tanks. The radargram in Figure 11a shows the targets at 1.0 m in depth, appearing as well-marked hyperbolic anomalies. The MT result in Figure 11b shows two anomalies related to the true targets’ locations, plus a secondary one in the midpoint, slightly below them. The centers of the target anomalies are somewhat vertically displaced from the actual centers of the targets. Similar to the prior horizontal target case, the ccTR signal processing shown in Figure 11c does not yield a good result for these horizontal targets. The strong amplitude anomaly shows the position of the storage tanks along the track, but it is not possible to recover the depth or size, because the anomaly extends well above the targets. Shallow anomalies, as indicated by red arrows, are also seen in these results, due to the geological features on the subsoil.

Target 8: Single vertical metallic storage tank. The results from the survey over the vertical metallic storage tank are shown in Figure 12. The radargram, seen in Figure 12a, shows an anomaly located at the upper portion of the target and mostly confined within the true lateral edges of the target. The MT image, seen in Figure 12b, shows a strong anomaly coinciding with the top of the target; however, there are also similar anomalies above the target position, at around 0.4 m in depth. Lower amplitude anomalies in the contrast function suggest some sort of continuity between these shallow anomalies and the ones related to the top boundary of the target, at a depth of 1.0 m. There are also low-amplitude anomalies at the target position, between 1.0 m and 1.9 m in depth. This whole set of anomalies may induce a misinterpretation based on the MT results, as they appear to be related to a single target located between 0.4 m and 1.9 m in depth. Figure 12c shows the image based on the ccTR signal, which provides good estimates for the position and lateral size. The depth is slightly underestimated in this case, with anomalies above the target not precluding the target identification.

6. Further Discussion

All buried targets considered in this study, except for Target 1, are either metallic targets—and thus the target bottom reflections are negligible relative to the target top reflections—or are hollow targets filled with air—thus the top and bottom reflections tend to become more overlapped, due to the faster phase velocity of the EM waves in air. In this latter case, the reflections may appear as a blurred, single hyperbolic feature and are reconstructed as only the shallowest diffraction. It should also be noted also that, in regard to MT, the relationship between the target and the signal wavelength is an important aspect to be considered: Since the problem is linearized based on the Born approximation, this implies that, for better accuracy of reconstruction, the targets should not be too large compared to the signal wavelength in the background soil (unless the dielectric contrast is very weak, which is not the case here). As alluded to before, MT is based on a backpropagation problem, in which the observed data are diffraction hyperbolas. Thus, the size of the reconstructed target also depends on the wavelength of the diffracted signal created by each target. In this sense, diffraction hyperbolas composed of a signal of wavelength, comparable to the target size, will converge to point images about the same size of this target. On the other hand, if the target is bigger, but the signal wavelength is short, the hyperbolic features will converge to spot features smaller than the actual target, the so-called “scattering centers” [50,51]. If two or more scattering centers corresponding to the same target can be identified in the radargram, they may be resolved individually depending on their separation. A different approach is employed in TR, where an iterative backpropagation process may allow different interfaces to be resolved, even when a reflection is barely seen in the raw radargram. This iterative process is responsible for modeling the signature from the bottom of the target.

MT is not a quantitative technique, since it does not always provide a precise estimate of the target permittivity or conductivity. The contrast function only gives the order of magnitude for the difference between the electric properties of the target and the surrounding medium. In general, MT is capable of better reconstructing low-contrast (“weak”) scatterers than TR, as in the case of targets 2 and 3. These are not strictly “weak scatterers” but have lower contrast compared to target 6, for example.

In this work, the artifacts shown in the images are from geological materials presented in the medium and generally, in a real field study, we cannot predict these “artifacts”. The background velocity was estimated through the two-way travel time observed in the radargram and known depth of each target. It is known that this velocity does not correspond exactly to the phase velocity in the medium, but it gives an acceptable velocity for imaging purposes. MT uses the information from the diffraction hyperbola as if it was generated in a homogeneous background, thus, the position of the target in the image corresponds to the two-way travel time registered in the radargram. On the other hand, TR considers different models and solves them through a synthetic FDTD model, simulating the environment that takes into account both the physical properties of the medium and the target. In this sense, TR may obtain better resolved images along the range (depth) than MT for some targets.

7. Conclusions

This paper provided a comparison between microwave tomography (MT) and a time-reversal-based technique (TR) applied to process GPR data obtained from a controlled site. The MT results were based on a first-order Born approximation with TVSD regularization. TR results were based on computing the cross-correlation between Modes 1 and 2 [32,33] of back-propagated TR signals. For the examples considered, MT gave better results for horizontal targets that have a length comparable to the signal wavelength. On the other hand, TR gave the best results for imaging vertical targets which have a length greater than the signal wavelength. Anomalies present in the MT images for targets are associated mostly to the top boundary of the targets with the above-mentioned characteristics, with the exception of targets 1 and 2. MT allowed better reconstruction of shallow or closely-spaced targets. In those cases, hyperbolic-like anomalies are still seen in TR imaging. Because of the somewhat complementary nature of their performance, the combination of MT and TR techniques can provide valuable information that either one, if used separately, might not able to provide, and hence, to potentially improve, GPR data interpretation. For the cases considered, both techniques yield artifacts in the images, which are caused by the presence of geological stratification and other inhomogeneities in the subsoil.

Finally, it is important to stress that MT and TR-based techniques are instantiated here using particular implementation choices. These choices do not fully exhaust the range of options for applying MT or TR techniques to GPR imaging. In particular, we reiterate that the TR-based results presented here were based on the use of the standard deviation of back-propagated TR signals. Different TR-based processing techniques are available that can also be applied to GPR problems, including TR-DORT [35], TR-MUSIC [42,52] and TR hybrid methods [53]. Likewise, several augmentations can be made to the MT implementation considered here, including the use of high-order Born approximations and iterative MT reconstructions [54,55,56]. The use of controlled-site data to compare all such options is beyond the scope of the present work and will be the subject of future studies.