Accounting for the Different Propagation Velocities for the Focusing and Time–Depth Conversion in a Layered Medium

Abstract

Featured Application

Imaging in stratified media and imaging of buried cavities.

Abstract

This contribution deals with the focusing and time–depth conversion of GPR data in a two-layer scenario. In a layered scenario, the different propagation velocities of the electromagnetic waves in different layers can induce distortion in the reconstruction. In particular, it is not possible, with common data processing, to achieve the focusing of both the targets embedded in the upper and deeper layers at the same time. We will present a strategy for mitigating this problem. The strength points of this strategy are its simplicity and the fact that it only requires a commercial code for the GPR data processing plus some steps easily implementable through a homemade code. This work extends a contribution of the same authors at the 20th International GPR conference held in Chanchung, China, in June 2024, and the conclusions are updated with three-dimensional numerical results, including depth slices determined directly in the spatial domain. Both numerical and experimental data will be shown.

1. Introduction

Layered scenarios constitute a topic of large interest within Ground-Penetrating Radar (GPR) prospecting. Common situations are asphalt or concrete pavements, within and under which metallic rods, pipes, voids, and several kinds of discontinuities could be displaced [1,2,3,4]. In other contexts, the internal coverage of the road tunnels also generates an intrinsically layered medium [5,6]. Another relevant case of large interest is the imaging of buried cavities [7,8], which is of interest for cultural and/or safety purposes. In particular, a cavity can be seen locally as a particular three-layered medium, according to the sequence “soil–air–soil”. Further cases of interest are archaeological heritage buried beyond the first layer of the soil or through-the-wall imaging, where the first layer of the investigated scenario is constituted by masonry.

In many cases, the inhomogeneity of the background scenario is neglected. However, a good choice of the background scenario provides better results. Two relevant challenges to this are the nonlinearity of the underlying inverse scattering problem and the computational burden required to afford it in a rigorous way. More specifically, while slow variations of the properties of the soil can be accommodated by linear inverse scattering algorithms [9,10] or reverse migrations [11,12,13], abrupt variations are more challenging. This can be the case of cavities and layered media with strong discontinuity of the electromagnetic characteristics between the layers. In particular, we consider here the case of a two-layered medium with a non-flat interface between the layers. Unlike the case of layers with a flat interface parallel to the air–soil interface [14,15,16,17], in this case, there is no analytical solution for the inverse scattering problem.

Two issues consequently arise, namely, the focusing of the buried targets and the time–depth conversion to apply to the achieved image. In homogeneous soil, these two issues are strictly linked to each other because the “true” value of the propagation velocity is the value driving both the best focusing and the most correct time–depth conversion [18,19,20]. However, things are trickier in stratified media, especially when a curved interface separates the layers. In particular, modeling is complicated and, in general, the commercial codes available for data processing are not able to accommodate such situations [21,22]. Moreover, we have to consider the fact that in GPR prospecting, the extension of the inversion domain can be electrically very large (even thousands of central wavelengths or more), so there is an obvious need for fast and efficient algorithms. This excludes, in most cases, the concrete possibility of exploiting nonlinear full wave inversions [23,24], which are incidentally also prone to the problem of the local minima [25].

For all these reasons, we have introduced two simple algorithms, aimed at the focusing and time–depth conversion in layered media, which we have labeled “joined migration” and “combined time–depth conversion” [4,26]. These algorithms correspond to a rough model of the propagation of the waves and of the scattering, but on the other hand, are computationally non-demanding, easy to implement, and do not require more than a commercial code for the post-processing plus a home-made code that is quite easy to implement. In spite of the simplicity of the model, in several situations of interest, the imaging of the background scenario can be meaningfully improved and the spatial relationship between the buried targets can be better reproduced with respect to the results achievable from more customary algorithms. With regard to this aspect, let us outline the fact that the customary non-invasive methods (diffraction hyperbolas, common mid-point) for retrieving the propagation velocity of the waves work correctly with regard to the shallower layer of a stratified soil, but they do not work satisfyingly any longer with regard to the deeper layers. Last but not least, it is impossible to find analytical methods for the diffraction curves in the layers of soil beyond the shallowest one, unless the interfaces between the layers are flat and parallel to the air–soil interface. In particular, in this contribution, we will implicitly rely on different kinds of measurements of the propagation velocity of the waves in the deeper layer, coming, e.g., from local ground truthing [27], Time Domain Reflectometry data [28], or (more theoretically) a cooperative target [29,30].

This work follows other papers of these authors on the same topic [4,26]. The novel aspect proposed here is the slicing in a two-layered medium achieved after joined migration and combined time-domain conversion. This will show that the proposed algorithm allows us to determine horizontal slices directly in the spatial domain instead of (as commonly performed) in the time domain.

In the next section, we will provide a brief description of the joined migration, namely the focusing algorithm adopted for targets embedded in a two-layered medium. In Section 3, we will describe the combined time-depth conversion algorithm, namely the time–depth conversion accounting for both the propagation velocities in the upper and lower mediums. In Section 4, the numerical results will be shown, based on 3D GPRMax (version 3.15) data [31]. In Section 5, some experimental results will be shown. Finally, our conclusions and future developments will be proposed.

2. The Joined Migration

As is well known, in homogeneous soil, the correct value of the propagation velocity of the waves allows the best focalization of the buried targets, achievable from migration or from an inverse scattering algorithm [18,19,20]. The same value also allows the most correct time–depth conversion available. In a multilayered medium, this happens only with regard to the shallowest layer (some exceptions are possible if the target in the first layer is somehow “shadowed” by the interface between the two first layers, but these cases are not frequent and, in most cases, recognizable from the data), and in particular, the diffraction hyperbolas generated by the targets in the shallowest layer provide the correct propagation velocity value according to the well-known procedure exploited in any customary GPR data processing [18]. However, this is no longer true for targets embedded in the deeper layers because the refraction of the electromagnetic layers at the buried interface deforms the diffraction curves (in general, in the deeper layers, they are no longer hyperbolas). In particular, it is possible to have a model in closed form of the diffraction curve only if the interface between the layers is horizontal [15,16,17] but no longer in the general case of a curved interface. In such a case, one should apply a numerical full-wave inversion in order to correctly invert the targets in both the upper and lower layers. However, a full-wave model leads to a nonlinear ill-posed problem that is mathematically difficult [32,33] and computationally challenging [34], and so here we propose a simplified solution for which the user only needs any commercial code for GPR data processing and some minimal further home-made software that is easy to implement and computationally non-demanding.

In particular, we will use trials to determine a value of propagation velocity that optimizes the focusing of the targets in the second layer, making use of a common migration algorithm in the time domain [20]. This optimal value is of course different from the value most suitable for focusing in the upper layer. Moreover, the optimal value to focus the targets in the second layer does not coincide with the actual value of the propagation layer in the second medium because the refraction of the waves at the buried interface is not accounted for. So, if we want to reconstruct the target in the second layer with good focusing and good depth positioning, then we need two propagation velocity values for this layer: the first one for achieving good focusing and the second “actual” one for achieving the correct time–depth conversion. However, while the optimal value for focusing is easily retrieved by trials, retrieving the actual propagation velocity in the lower layer from the GPR data is a more complex problem. At this stage, as said, we implicitly assume there is a possibility of conducting localized ground truthing or a TDR measurement (possibly locally removing the upper layer). That said, the joined migration consists of performing and retaining two migration results, one (let us label it M₁) optimal for focusing the targets in the upper layer and another one (let us label it M₂) optimal for focusing the anomalies buried in the lower layer. Afterward, we will create a matrix whose values are that of M₁ above the layer interface as detected in the time domain and those of M₂ in the complementary part, i.e., below the interface. More details on this joining can be found in refs. [4,26].

3. The Combined Time–Depth Conversion

After achieving joined migration in the time domain, we have to convert it into a graph represented by the abscissa and the spatial depth. The schematic of Figure 1 displays how this can be achieved. In particular, Figure 1A schematizes a result in abscissa “x” and time “t”, and the (not-flat) buried interface is represented in red: the representation is piecewise constant because the pixels are put into evidence. Three targets are embedded in the soil, represented by three small circles. After performing the joined migration, or let us say more generally, any elaboration on the data that provides a focused result in the (x–t) plane, we have to represent the buried scenario in the abscissa and spatial depth, i.e., in the (x–z) plane. The essential point is that the vertical size of the pixels in the x–t plane is the same in the two layers, but the corresponding vertical size in the x–z domain is not the same because it is given by the well-known round-trip formula, which should, however, be applied separately to the upper and lower layers. Therefore, we have:

(1)$\left\{\begin{array}{c}\Delta z_1={\displaystyle \frac{c_1\Delta t}{2}}\\ \Delta z_2={\displaystyle \frac{c_2\Delta t}{2}}={\displaystyle \frac{c_2}{c_1}}\Delta z_1\end{array}\right.$

where $c_1$ is the actual propagation velocity of the electromagnetic waves in the upper layer and $c_2$ is the actual propagation velocity in the lower layer. $\Delta t$ is the time step of the data, $\Delta z_1$ is the ideal vertical thickness of the pixel in the upper layer, and $\Delta z_2$ is the homologous quantity in the lower layer. So, the converted pixels should appear as in Figure 1B. More precisely, Figure 1B represents the case $c_1>c_2$. In the opposite case, the pixels of the lower layer should be vertically dilatated with respect to those of the upper layer. As can be seen from Figure 1A, in the time domain, the central and left-hand targets appear at the same depth, whereas the left-hand target appears deeper. Instead, in Figure 1B, i.e., “in reality”, the left-hand and central targets are at the same depth, whereas the right-hand one is shallower. Whatever the case, the final result has to be represented by a matrix, where the values are necessarily associated with a regular rectangular grid, and this is not the case in Figure 1B. In order to achieve pixels of the same size, interpolation is needed so that the cells of the upper medium and those of the lower medium have the same size.

With some easy algebra, this means that in the case of $c_1>c_2,$ for each column of the “matrix”, starting from the air–soil interface and ending at the interface between the two layers, the initial number of pixels $n\left(j\right)$ (indexed after the column index j) should be interpolated into $N\left(j\right)={\displaystyle \frac{c_1}{c_2}}n\left(j\right)$ pixels. In the dual case $c_2>c_1,$ for each column of the “matrix”, starting from the interface between the two layers and ending at the (local) spatial bottom scale, the initial number of pixels $n\left(j\right)$ should be interpolated into $N\left(j\right)={\displaystyle \frac{c_2}{c_1}}n\left(j\right)$ pixels. In general, the ratio ${\displaystyle \frac{c_1}{c_2}}$ (${\displaystyle \frac{c_2}{c_1}})$ is not an integer number, and so the integer number closest to the quantity ${\displaystyle \frac{c_1}{c_2}}n\left(j\right) \left({\displaystyle \frac{c_2}{c_1}}n\left(j\right)\right)$ should be chosen. The resulting residual distortion is customarily negligible. After interpolation, the pixels will appear as in Figure 1C. This set of pixels can be completed with zero padding as shown in Figure 1D. In this way, we compose a complete matrix that depicts the scenario in the x–z domain. The necessity of zero padding corresponds to the fact that the time bottom scale of the data is constant vs. the abscissa but the corresponding spatial bottom scale is not, and this is natural because the average velocity along any vertical path is not the same vs. the abscissa. Figure 1 also allows us to appreciate how different velocities can distort the spatial relationships among the targets. It is also easily understood that if the three targets were not point-like but had some vertical thickness, the ratios among their thicknesses would appear distorted without a combined time–depth conversion.

The flux diagram of the algorithm is provided in Figure 2.

4. Numerical Results

We have considered a three-dimensional two-layer scenario, simulated by means of the open access GPRMax code [31], based on the method of the finite difference in time domain FDTD [35]. The simulated scenario is shown in Figure 3.

In particular, the observation surface is extended to 4.2 × 3 m, and three pipes are buried and labeled P₁, P₂, and P₃, respectively. Pipe P₁ is the upper one: it is parallel to the long side of the box and is buried at a depth of 45 cm. Pipes P₂ and P₃ are both at the depth level of 90 cm. Pipe P₂ follows a direction orthogonal to the upper pipe P₁, whereas P₃ is slightly tilted (7.91 degrees) with respect to P₂. The background medium is composed of two layers separated by a flat but slanted interface that starts from the air–soil interface on the left-hand side and arrives at a depth of 1 m on the right-hand side. The dielectric permittivity is equal to 4 in the upper layer and 16 in the lower one, which means that the propagation velocity of the electromagnetic waves is equal to 15 cm/ns in the upper layer and 7.5 cm/ns in the lower layer. No magnetic properties are present, and a small electrical conductivity equal to 0.005 S/m has been introduced in both layers. The data have been taken along an orthogonal grid with interline spacing of 30 cm and in-line spacing of 2 cm. The time bottom scale is 80 ns and the time step is 0.0192 ns. The pixels used for implementing the FDTD were sized 1 cm in any direction. The transmitting and receiving antennas are two dipoles located 1 cm from the air–soil interface and displaced at negligible distances from each other. The radiated signal is a Ricker wavelet whose 3 dB band ranges from 137.5 MHz to 450 MHz. The data have been processed with the Reflexw code [36] through zero timing, background removal up to 5 ns (in order to avoid distortion of the shallowest pipe), linear and exponential gain vs. depth, Butterworth filtering, and Kirchoff migration (with two different velocities as mentioned). We have put into the inversion algorithm c₁ = 16 cm/ns for the upper layer, as retrieved from the diffraction hyperbola of pipe P₃ (the exact value is 15 cm/ns instead). With regard to the second layer, we have exploited a value of 12 cm/ns for the focusing, retrieved by trials, and then we have exploited a value of 8 cm/ns for the time–depth conversion, voluntarily introducing a small error with respect to the exact value of 7.5 cm/ns. In this way, some uncertainty in the measurement of the propagation velocity of both layers has been implicitly accounted for.

In Figure 4, a Bscan processed according to customary data elaboration is shown. In particular, the 6th (out of 11) Bscan taken along the long side of the grid is shown. The length of the Bscan is 4.2 m.

In Figure 4, we have also converted the return time into the spatial depth according to the estimated propagation velocity in the upper medium, namely 16 cm/ns. From Figure 4, we can see that we achieved good focusing of the pipe in the upper medium, namely P₃. However, the focusing of pipe P₂, embedded entirely in the lower layer, is not satisfying at all. Moreover, while the depth of P₃ is reliably retrieved (indeed, its depth appears equal to 100 cm instead of the actual value of 90 cm because of the overestimation of the propagation velocity), the apparent depth of P₂ is 150 cm, which is meaningfully larger than the actual value. This happens because the propagation velocity of the waves in the deeper layer is lower than that in the upper layer. We also appreciate the signature of pipe P₁, crossing the two layers. Since this pipe is parallel to the direction of the Bscan, it appears as a straight line with regard to the part of the pipe contained in the upper layer. However, when the pipe enters the second layer on the left-hand side, it appears to progressively dip, coherently with the progressive decrease in the average velocity of the waves along their path from the transmitting dipole up to the pipe and back to the receiving dipole. Finally, the interface between the two buried layers is clearly visible as a slanted line. In Figure 5, the same Bscan of Figure 4 is shown, but in this case, focusing was achieved according to the propagation velocity allowing for better focusing of the pipe in the lower layer (P₂). This velocity, found by trials, is 12 cm/ns. In Figure 5, we can see that good focusing of pipe P₂ is achieved, but obviously at the cost of defocusing pipe P₃.

In Figure 5, the same time–depth conversion as in Figure 4 is applied. The apparent bending of pipe P₁ also remains, as can be seen. It is also to be said that the problem is more complex than what Figure 5 represents. For example, if we had more pipes in the second layer, the propagation velocity allowing the best focusing of one of them would not necessarily guarantee the best focusing of the other ones. We are exploiting a heuristic model, and so in more complex cases, we should choose a compromised solution optimizing the focusing of the presumably most important targets while minimizing the defocusing of the other ones as much as possible. In particular, targets displaced in the second layer but very close to the interface between the two layers would be indeed better focused making use of the propagation velocity of the upper layer. In fact, in this case, the propagation of the waves occurs mostly in the upper layer.

In Figure 6, we propose the application of a combined time–depth conversion to the result of Figure 4. As can be seen, the combined time–depth conversion allows us to restore the correct depth of the pipes because now pipes P₂ and P₃ appear to be at about the same depth. Moreover, the apparent bending of pipe P₁ is much more mitigated. However, pipe P₂ appears strongly distorted because we have not applied joined migration.

In particular, if a slicing procedure were applied based on all the Bscans processed in this way, we would achieve severe distortion of the image of pipe P₂.

In Figure 7, we propose the result of joined migration without combined time–depth conversion. As can be seen, this time pipes P₂ and P₃ are both well-focused. However, P₂ is imaged deeper than its actual position, and P₁ is again bent on the left-hand side. The conclusion is that only the application of both joined migration (or more generally, an algorithm devoted to the focusing of the target in both layers in the time domain) joined to combined time–depth conversion can provide good imaging in the case at hand. In Figure 8, both joined migration and combined time–depth conversion are applied. As said, the propagation velocity exploited for the focusing of pipe P₂ is 12 cm/ns, whereas that exploited for the combined time–depth conversion is 8 cm/ns. As can be seen from Figure 8, now both pipes P₂ and P₃ are well focused and appear to be displaced at the same depth level as they are. Moreover, the bending effect of pipe P₁ is meaningfully mitigated. It can be seen that a small distortion is introduced with regard to pipe P₂, whose image is no longer “roundish” and is influenced by the shape of the interface between the two layers. It can be also seen that P₂ now appears smaller than pipe P₃: this, however, is natural because it is related to the smaller wavelength in the second layer, which drives a better resolution. The same processing applied to the Bscan in Figure 8 was then applied to all of the Bscans, identifying the interface in the time domain between the two media each time. Of course, the shape of this interface changes for the Bscans directed along the short side of the observation grid with respect to those directed along the long side and generally changes from Bscan to Bscan. This processing produced 26 Bscans (11 along the long side and 15 directed along the short side) all properly focused and time–depth-converted, and from them, we retrieved depth slices. Let us specify that since the Bscans have undergone a combined time–depth conversion, they are vertically calibrated directly in meters, and so the slices achieved from them are determined directly in the spatial domain, unlike the customary slices that are determined in the time domain.

Notwithstanding, the slicing in the spatial domain does not require the implementation of an algorithm. In particular, the slices can be produced with any commercial code accepting ASCII files as input: in this way, we only have to input the spatial bottom scale in meters to the code. Then, the code will “understand” the bottom scale as a number expressed in nanoseconds, but this is not relevant. A sequence of depth slices achieved from Bscans processed with a joined migration and combined time–depth-converted Bscan is shown in Figure 9. From the slices at subsequently deeper levels, we note a sort of belt directed along the short side of the slice, shifting from the left-hand side to the right-hand side progressively vs. the depth. This echo is the reflection provided by the interface between the two media at the considered depth level. Due to the shape of the interface, descending from the left- to the right-hand side of each (long) Bscan, this belt shifts toward the right-hand side at progressively deeper levels.

Beyond this reflection, pipe P₁ is well identified between the depth levels of 40 and 50 cm, whereas pipes P₂ and P₃ are imaged between the depth levels of 90 and 100 cm, and correctly, they are at the same depth level. Pipe P₃ provides a stronger echo because it is displaced in the upper layer and there is no further interface beyond the air–soil one between the antennas and the pipe. On the other hand, it appears broken in the middle because of the overflying pipe P₁. This breaking effect is instead not evident with regard to pipe P₂ because the further interface between the upper and lower layer introduces more complicated scattering phenomena.

For comparison, in Figure 10, we also propose the imaging achieved from Bscans focused through joined migration but not processed with a combined time–depth conversion. The depth levels are evaluated according to the propagation velocity in the upper layer. From Figure 10, we can see that while pipe P₃ is still correctly imaged, pipe P₂ is imaged more intensely between 140 and 150 cm. Moreover, the non-compensated apparent bending of pipe P₁ on the left-hand side generates an artifact resembling a target different from the pipe itself and extending along a depth of up to 60 cm or more, which makes it tricky to determine the correct interpretation of the data.

5. Experimental Results

An experimental setup was organized at the Laboratory of Hydrogeophysics of the Institute of Methodologies for Environmental Analyses IMAA-CNR, in Marsico Nuovo (Potenza, Italy). In particular, a small tank of 70 × 50 × 45 cm was filled up, constituting a two-layered medium with a curved interface between the layers. The lower layer constituted wet sand, whereas the upper one constituted dry gravel. In Figure 11, a photograph of the tank during its filling is shown. Four metallic targets were inserted according to the scheme shown in Figure 12. In particular, the targets “of interest” for our purposes are A and B, whereas targets C and D are exploited for the direct evaluation of the propagation velocity. In fact, they are displaced in a flat zone of the interface between the two layers at a known depth. Consequently, the return time from their echoes will be exploited in order to evaluate the propagation velocity of the waves in the two layers.

The GPR system exploited was a Ris Hi-mode equipped with an antenna at a nominal central frequency of 2 GHz. A single Bscan crossing all the metallic targets (which are essentially rods parallel to each other) was gathered. The measurement line was 63 cm long, the time bottom scale was 16 ns, the spatial step was 2 mm, and the time step was 0.0156 ns. In Figure 13, the preprocessed (not focused) data are shown. Pre-processing constituted zero timing, subtracting the average on 51 traces (i.e., shifting background removal performed on 51 traces), gain vs. depth, and 1D Butterworth filtering, plus a final time cut reducing the time bottom scale to 8 ns. In particular, in Figure 13, we indicate the vertexes of the diffraction curves of the four metallic targets with the same symbols used in Figure 12. In Figure 13, the interface between the two layers is clearly appreciated. In particular, the subtracted average has been adopted in order to avoid the production of spurious horizontal belts due to the background removal [37] on all the traces. Offline from this preprocessing, in order to ensure a more reliable evaluation of the propagation velocity in the second medium, in Figure 14, we show the result of another preprocessing performed on the same data, where the subtracted average has been replaced by background removal on all the traces. This has generated the spurious horizontal belt clearly visible in the figure. The imaging provided by Figure 14 is certainly worse than that provided by Figure 13; however, with regard to targets C and D, the horizontal belts make their time depth levels clearer. In Figure 14, we have also indicated target B, whereas target A disappears because it is covered by the strength of the horizontal belt generated by the diffraction curve of target D.

As shown in Figure 12, the depth of target D is equal to 27 cm (bottom of the first layer) and the depth of target C is equal to 45 cm, which represents 18 cm in the deeper layer. Their time depth levels are instead (from Figure 14) 2.8 ns and 6.2 ns, respectively (which means 3.4 ns of time depth in the second layer). From these data, we retrieve a propagation velocity of 19.2 cm/ns for the upper layer (made up of dry gravel) and 10.6 cm/ns for the second layer (made up of wet sand). The diffraction hyperbola of target B in Figure 12 provides a propagation velocity of 20 cm/ns for the upper layer, which is in good agreement with the values retrieved from target D in Figure 14. Instead, from Figure 13, the diffraction curve of target A matches the propagation velocity of 13 cm/ns, with a meaningful discrepancy with respect to the evaluation provided by target C, as expected.

In Figure 15, the result of joined migration is shown. For joined migration, the propagation velocities of 20 cm/ns and 13 cm/ns have been exploited, namely, those providing the best focusing of targets B and A, respectively. Target A is focused correctly but appears to be surrounded by further two targets on the left-hand and right-hand sides: this is an effect of the subtracted average, which mitigates (namely, reduces its extension) but does not completely erase the production of the horizontal belt. In particular, a similar effect is also visible with regard to target B.

Instead, targets C and D are “lost” because they are very close to the interface between the layer and the bottom of the tank, respectively.

In Figure 16, the result of a combined time–depth conversion applied to the data of Figure 15 is applied. For the time–depth conversion, we instead exploited the values 20 cm/ns (provided by the diffraction hyperbolas) and 10.6 cm/ns, somehow retrieved in another way. As can be seen, the depths of both targets A and B are now correctly retrieved.

The seam line visible in Figure 15 and Figure 16 is an artifact due to the separation line between the zone migrated with the “shallower” propagation velocity and the “deeper” propagation velocity. It is well-recognized because it is a discontinuity that is too spatially sharp to be ascribable to the real scenario. However, due to the improvement in the focusing and the displacement of the targets of interest achieved in both layers, this is a tolerable artifact.

6. Conclusions and Future Developments

In this contribution, we have shown some potentialities of the application of joined migration and a combined time–depth conversion to GPR data, which has allowed us to improve the imaging in a two-layered medium. The algorithm improved both the vertical Bscans and the slices because the different propagation velocities in the two layers have been accounted for. In particular, we have proposed a simplified algorithm that neglects the vectorial refraction of the waves at the interface but can account for the relatively complex shapes of the interface between the two layers, which is a case for which no analytic solution is currently available.

Theoretically, the exposed considerations also easily extend to the case of N-layered media, but it is naturally progressively more difficult (and less reliable) to retrieve the electromagnetic characteristics of the subsequently deeper layers. Whether this is practically possible depends on the particular application, the losses, the thickness of the layers, and the difference between the electromagnetic characteristics of the adjacent layers, namely the amplitude of the first kind of discontinuity jump shown by the propagation velocity of the waves.

Several possible future developments are being thought of. First of all, the mathematical description of the interface between two different layers (or the top and bottom of a buried cavity) can be made faster by means of some computer graphic algorithm enabling the user to rapidly pick up these curves in the time domain by clicking points on their computer screen with the mouse and/or by moving the mouse along the interface on the screen. This is especially important if the interface between two buried layers does not have a simple shape.

A further development related to the efficiency of the algorithm is the automatization of the exportation of the Bscans from a commercial code (we exploited the Reflexw, version 9.1.3) to an external code for the interpolations and the zero paddings (we dialed a home-made code in the MATLAB (version 9.3.0) environment) and then the re-importation into commercial code for the slicing. In particular, we exported the Bscans in ASCII format from Reflexw to MATLAB in order to apply joined migration and the combined time–depth conversion and then reimported the combined files again in ASCII format from MATLAB to Reflexw in order to perform the slicing. This implied the loss of the headers of the files each time (essentially the length of the Bscans and their bottom scale) that have been manually rewritten each time.

Another aspect of interest is the retrieval of the permittivity of the second layer from the same GPR data. This aspect has been studied in the case of a sequence of flat layers [15,16,17], where analytical calculations are possible. In the case of a curved interface, no analytical solution is available, but notwithstanding, suitable forward modeling could help in retrieving the propagation velocity of the waves in the second layer. Forward modeling can be implemented, making use of the above-quoted open-access code GPRMax [31], using trials to determine the value of permittivity of the second medium and allowing the best graphical matching of the diffraction curves visible in the second layer (and presumably attributable to small targets). This analysis might also provide an evaluation of the available precision on the basis of the sensibility of the diffraction curve to the trial value of dielectric permittivity. This point is relevant because we might not have the possibility of performing local ground truthing and we might not have a TDR probe at our disposal or, even with it, the depth of the first layer could result in a difficult local excavation up to the second layer.

Future work is also necessary with regard to a specific comparison of this method with more analytic algorithms, like time reverse migration or regularized linear inverse scattering methods. It is expected that more analytic methods perform better for slight discontinuities between the layers (because analytical models are more mathematically rigorous than the method exploited here), but strong discontinuities can be problematic (because analytical models are approximated). In particular, it is certainly rare to find analytic focusing methods applied to a layered medium with strong discontinuities of the permittivity between two adjacent layers. Moreover, we have not found a combined time–depth conversion or something similar applied for achieving the final image. This is an important point because, regardless of which focusing algorithm one exploits, if the layers show a different propagation velocity of waves, the spatial depth formally considered is not the same vs. the abscissa, and this should be accounted for in order to achieve a more correct time–depth conversion.

These considerations are also valid for the case of soil in which the characteristics change more smoothly (e.g., because of some moisture gradient), producing a continuous variation in the permittivity either along the horizontal or vertical directions. In these cases, the focusing of targets will likely be better accommodated with time reverse migration, but the issue of correct time–depth conversion remains, and a suitable combined time–depth conversion is the tool that can address this problem properly.

Author Contributions

Conceptualization, R.P.; methodology, I.C., G.E., L.C. and G.D.M.; software, D.Y.; validation, D.Y. and R.P.; formal analysis, R.P., I.C., L.C., G.D.M., G.M. and D.Y.; investigation, I.C., L.C., G.D.M. and G.M.; resources, G.M.; data curation, D.Y.; writing—original draft preparation, R.P.; writing—review and editing, I.C., G.E. and D.Y. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Conflicts of Interest

Author Gianfranco Morelli was employed by the company Geostudi Astier s.r.l. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Representation of a scenario with two layers. (A) Pictorial of the matrix corresponding to the data in abscissa and time domain. (B) Corresponding pixels in abscissa and depth. (C) Interpolated equally sized pixels. (D) Zero-padded final complete matrix.

Figure 2. Flux diagram of the algorithm.

Figure 3. 3D Simulation scenario. (A) Geometric details of the perspective view scene. (B) Geometric details of the frontal view scene. (C) Geometric details of the top-down view scene.

Figure 4. Sixth processed Bscan on the long side. Focusing achieved after the propagation velocity of the upper layer. Common time–depth conversion is applied.

Figure 5. Sixth processed Bscan on the long side. Focusing achieved after the propagation velocity allowing the best focusing of the pipe in the lower layer. Common time–depth conversion is applied.

Figure 6. Sixth processed Bscan on the long side. Combined time–depth conversion is applied to the image in Figure 4.

Figure 7. Sixth processed Bscan on the long side. Joined migration is performed and applied but not combined time–depth conversion.

Figure 8. Sixth processed Bscan on the long side. Both joined migration and combined time–depth conversion are applied.

Figure 9. Depth slices achieved from all the Bscans after joined migration and combined time–depth conversion.

Figure 10. Depth slices achieved from all the Bscans after joined migration but without combined time–depth conversion.

Figure 11. Tank filled up with sand gravel and metallic targets.

Figure 12. Schematic of the ground truth.

Figure 13. Preprocessed data with subtracted average on 51 traces.

Figure 14. Preprocessed data with background removal.

Figure 15. Result of joined migration.

Figure 16. Result of the combined time–depth conversion.

References

1. Elseicy, A.; Alonso-Diaz, A.; Solla, M.; Rasol, M.; Santos-Assunçao, S. Combined Use of GPR and Other NDTs for Road Pavement Assessment: An Overview. Remote Sens.; 2022; 14, 4336. [DOI: https://dx.doi.org/10.3390/rs14174336]

2. Rodes, J.P.; Reguero, A.M.; Perez-Gracia, V. GPR Spectra for Monitoring Asphalt Pavements. Remote Sens.; 2020; 12, 1749. [DOI: https://dx.doi.org/10.3390/rs12111749]

3. Joshaghani, A.; Shokrabadi, M. Ground penetrating radar (GPR) applications in concrete pavements. Int. J. Pavement Eng.; 2021; [DOI: https://dx.doi.org/10.1080/10298436.2021.1954182]

4. Persico, R.; Morelli, G. Combined Migrations and Time-Depth Conversions in GPR Prospecting: Application to Reinforced Concrete. Remote Sens.; 2020; 12, 2778. [DOI: https://dx.doi.org/10.3390/rs12172778]

5. Kang, X.; Xie, X.; Zou, J.; Lei, S.; Zhou, B.; Dong, Q. A new stage control device for open TBM tunnel lining support in water-rich weakly cemented soft rock. Int. J. Rock Mech. Min. Sci.; 2024; 179, 105754. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2024.105754]

6. Prego, F.J.; Solla, M.; Núñez-Nieto, X.; Arias, P. Assessing the Applicability of Ground-Penetrating Radar to Quality Control in Tunneling Construction. J. Constr. Eng. Manag.; 2015; 142, 06015006. [DOI: https://dx.doi.org/10.1061/(ASCE)CO.1943-7862.0001095]

7. Utsi, E.; Colls, K. The GPR investigation of the Shakespeare family graves. Archaeol. Prospect.; 2017; 24, pp. 335-352. [DOI: https://dx.doi.org/10.1002/arp.1573]

8. Persico, R.; D’Amico, S.; Matera, L.; Colica, E.; De Giorgio, C.; Alescio, A.; Sammut, C.; Galea, P. GPR Investigations at St John’s Co-Cathedral in Valletta. Near Surf. Geophys.; 2019; 17, pp. 213-229. [DOI: https://dx.doi.org/10.1002/nsg.12046]

9. Pierri, R.; Leone, G.; Soldovieri, F.; Persico, R. Electromagnetic inversion for subsurface applications under the distorted Born approximation. Nuovo C.; 2001; 24, pp. 245-261.

10. Salucci, M.; Oliveri, G.; Massa, A. GPR Prospecting Through an Inverse-Scattering Frequency-Hopping Multifocusing Approach. IEEE Trans. Geosci. Remote Sens.; 2015; 53, pp. 6573-6592. [DOI: https://dx.doi.org/10.1109/TGRS.2015.2444391]

11. Alashloo, S.M.; Ghosh, D.P.; Bashir, Y.; Yusoff, W.W. Depth Imaging Enhancement Using Reverse Time Migration. IOP Conf. Ser. Earth Environ. Sci.; 2017; 88, 012017. [DOI: https://dx.doi.org/10.1088/1755-1315/88/1/012017]

12. Zhang, X.; Li, R.; Cui, L.; Liu, D.; Liu, G.; Zhang, Z. Least squares reverse time migration imaging with illumination preconditioned based on improved PRP conjugate gradients. Sci. Rep.; 2023; 13, 13623. [DOI: https://dx.doi.org/10.1038/s41598-023-40578-8]

13. Patsia, O.; Giannopoulos, A.; Giannakis, I. Background Removal, Velocity Estimation, and Reverse-Time Migration: A Complete GPR Processing Pipeline Based on Machine Learning. IEEE Trans. Geosci. Remote Sens.; 2023; 61, 2003311. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3300276]

14. Chew, W.C. Waves and Fields in Inhomogeneous Media; Institute of Electrical and Electronics Engineers: Piscataway, NJ, USA, 1995.

15. Lambot, S.; Slob, E.C.; van den Bosch, I.; Stockbroeckx, B.; Vanclooster, M. Modeling of Ground-Penetrating Radar for Accurate Characterization of Subsurface Electric Properties. IEEE Trans. Geosci. Remote Sens.; 2004; 42, pp. 2555-2568. [DOI: https://dx.doi.org/10.1109/TGRS.2004.834800]

16. Dix, C.H. Seismic velocities from surface measurements. Geophysics; 1955; 20, pp. 68-86. [DOI: https://dx.doi.org/10.1190/1.1438126]

17. Persico, R.; Leucci, G.; Matera, L.; De Giorgi, L.; Soldovieri, F.; Cataldo, A.; Cannazza, G.; De Benedetto, E. Effect of the height of the observation line on the diffraction curve in GPR prospecting. Near Surf. Geophys.; 2014; 13, pp. 243-252. [DOI: https://dx.doi.org/10.3997/1873-0604.2014042]

18. Lesselier, D.; Duchene, B. Wavefield inversion of objects in stratified environments: From back-propagation schemes to full solutions. Review of Radio Science 1993–1996; Stone, R. Oxford University Press: Oxford, UK, 1996.

19. Stolt, R.H. Migration by Fourier Transform. Geophysics; 1978; 43, pp. 23-48. [DOI: https://dx.doi.org/10.1190/1.1440826]

20. Schneider, W.A. Integral formulation for migration in two and three dimensions. Geophysics; 1978; 43, pp. 49-76. [DOI: https://dx.doi.org/10.1190/1.1440828]

21. Sala, J.; Linford, N. Processing stepped frequency continuous wave GPR system so to obtain maximum value from archaeological data sets. Proceedings of the XIII International Conference on Ground Penetrating; Lecce, Italy, 21–25 June 2010; pp. 74-79.

22. Dong, Z.; Feng, X.; Zhou, H.; Zou, L.; Sato, M. 3D Migration Depth Focus Velocity Analysis of Hand-Held Ground Penetrating Radar. Geosciences; 2022; 12, 178. [DOI: https://dx.doi.org/10.3390/geosciences12040178]

23. Catapano, I.; Crocco, L.; Persico, R.; Pieraccini, M.; Soldovieri, F. Linear and Nonlinear Microwave Tomography Approaches for Subsurface Prospecting: Validation on Real Data. IEEE Trans. Antennas Wirel. Propag. Lett.; 2006; 5, pp. 49-53. [DOI: https://dx.doi.org/10.1109/LAWP.2006.870363]

24. Wapenaar, K.; Broggini, F.; Slob, E.; Snieder, R. Three-Dimensional Single-Sided Marchenko Inverse Scattering, Data-Driven Focusing, Green’s Function Retrieval, and their Mutual Relations. Phys. Rev. Lett.; 2013; 10, 084301. [DOI: https://dx.doi.org/10.1103/PhysRevLett.110.084301]

25. Norton, S.J. The inverse-scattering problem and global convergence. J. Acoust. Soc. Am.; 2005; 118, pp. 1534-1539. [DOI: https://dx.doi.org/10.1121/1.1979407]

26. Persico, R.; Morelli, G.; Esposito, G.; Catapano, I.; Capozzoli, L.; De Martino, G.; Yang, D. An innovative time-depth conversion for the management of buried scenarios with strong discontinuities. J. Appl. Geophys.; 2024; 227, 105435. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2024.105435]

27. Conyers, L. Interpreting Ground-Penetrating Radar for Archaeology; 1st ed. Routledge Taylor and Francis Group: New York, NY, USA, 2012; ISBN 978-1-61-132217-0

28. Persico, R.; Pieraccini, M. Measurement of dielectric and magnetic properties of Materials by means of a TDR probe. Near Surf. Geophys.; 2018; 16, pp. 118-126. [DOI: https://dx.doi.org/10.3997/1873-0604.2017046]

29. Allen, C.T.; Shi, K.; Plumb, R.G. The use of ground-penetrating radar with a cooperative target. IEEE Trans. Geosci. Remote Sens.; 1998; 36, pp. 1821-1825. [DOI: https://dx.doi.org/10.1109/36.718649]

30. Soldovieri, F.; Prisco, G.; Persico, R. A strategy for the determination of the dielectric permittivity of a lossy soil exploiting GPR surface measurements and a cooperative target. J. Appl. Geophys.; 2009; 67, pp. 288-295. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2008.09.007]

31. Warren, C.; Giannopoulos, A.; Giannakis, I. gprMax: Open source software to simulate electromagnetic wave propagation for Ground Penetrating Radar. Comput. Phys. Commun.; 2016; 209, pp. 163-170. [DOI: https://dx.doi.org/10.1016/j.cpc.2016.08.020]

32. Colton, D.; Kress, R. Inverse Acoustic and Electromagnetic Scattering Theory; Springer: Berlin/Heidelberg, Germany, 1992.

33. Bertero, M.; Boccacci, P. Introduction to Inverse Problems in Imaging; Institute of Physics Publishing: Bristol, UK, Philadelphia, PA, USA, 1998.

34. Meles, G.; Greenhalgh, S.; van der Kruk, J.; Green, A.; Maurer, H. Taming the non-linearity problem in GPR full-waveform inversion for high contrast media. J. Appl. Geophys.; 2012; 78, pp. 31-43. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2011.12.001]

35. Sarkar, D. FDTD Analysis of Guided Electromagnetic Wave Interaction with Time-Modulated Dielectric Medium; Springer: Berlin/Heidelberg, Germany, 2022.

36. Available online: https://www.sandmeier-geo.de/reflexw.html (accessed on 8 November 2024).

37. Rashed, M.; Harbi, H. Background matrix subtraction (BMS): A novel background removal algorithm for GPR data. J. Appl. Geophys.; 2014; 106, pp. 154-163. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2014.04.022]