High-Resolution Coherency Functionals for Improving the Velocity Analysis of Ground-Penetrating Radar Data

Abstract

We aim at verifying whether the use of high-resolution coherency functionals could improve the signal-to-noise ratio (S/N) of Ground-Penetrating Radar data by introducing a variable and precisely picked velocity field in the migration process. After carrying out tests on synthetic data to schematically simulate the problem, assessing the types of functionals most suitable for GPR data analysis, we estimated a varying velocity field relative to a real dataset. This dataset was acquired in an archaeological area where an excavation after a GPR survey made it possible to define the position, type, and composition of the detected targets. Two functionals, the Complex Matched Coherency Measure and the Complex Matched Analysis, turned out to be effective in computing coherency maps characterized by high-resolution and strong noise rejection, where velocity picking can be done with high precision. By using the 2D velocity field thus obtained, migration algorithms performed better than in the case of constant or 1D velocity field, with satisfactory collapsing of the diffracted events and moving of the reflected energy in the correct position. The varying velocity field was estimated on different lines and used to migrate all the GPR profiles composing the survey covering the entire archaeological area. The time slices built with the migrated profiles resulted in a higher S/N than those obtained from non-migrated or migrated at constant velocity GPR profiles. The improvements are inherent to the resolution, continuity, and energy content of linear reflective areas. On the basis of our experience, we can state that the use of high-resolution coherency functionals leads to migrated GPR profiles with a high-grade of hyperbolas focusing. These profiles favor better imaging of the targets of interest, thereby allowing for a more reliable interpretation.

1. Introduction

The Ground-Penetrating Radar (GPR) method is based on electromagnetic wave (EM) propagation and its response to changes in the EM properties of the subsurface. An EM impulse (GPR wave) generated by a transmitting antenna (transmitter, TX) propagates in the subsurface and is reflected by an interface or “scatters” on a discrete object. The energy back-reflected/scattered to the surface is recorded by a receiving antenna (receiver, RX) [1,2]. The estimation of the propagation velocity field of GPR waves (hereinafter velocity field) is an important issue to obtain a reliable and focused data migration, a crucial processing step to achieve a correct final interpretation [2,3,4,5]. This is particularly true in complex subsurface contexts, where the existence of numerous and different targets (e.g., walls, voids, pavements, fractures, metal rebars, burials) interspersed in heterogeneous materials may generate several reflections and diffraction patterns. In such cases, adopting a constant velocity for migration may produce unsatisfactory results, potentially reducing the data interpretation. The subsurface contexts whose imaging may benefit from an accurate velocity field estimation are thereby numerous, from archeology to civil engineering and structural geology settings. Velocity analysis on GPR data is usually carried out using an approximate approach in which a synthetic hyperbola of known velocity is imposed on hyperbolic-shaped real data, i.e., a real diffraction evident in the GPR profile. By a visual fitting of the curvature of the synthetic hyperbola with the observed real diffraction, the GPR wave velocity is determined. Alternatively, the semblance functional [6] on a data subset extracted from the radar profile is used. These operations are performed on a point or in a selected area of the GPR profile, where the signal-to-noise ratio (S/N) of the diffractions is high. Recently, automatic tools of hyperbola detection have been introduced assisting the operators in data analysis and processing [7,8]. However, avoiding false identifications may require some manpower effort and expertise. The resulting velocity field is constant or at most varying only with depth (1D). The estimated velocity is then used for data migration, thus collapsing the diffractions and moving the observed events in the correct position on the migrated section. When a basic processing sequence is adopted the estimated velocity is frequently used just for a vertical time-to-depth conversion. This way of proceeding disregards many factors that can actually result in an inaccurate velocity field estimation, thereby not favoring the production of an optimal migrated or depth converted section. The most important of these factors is that the propagation velocity of the electromagnetic radiation can vary both vertically and laterally along the survey profile. Although in many case studies the GPR investigated medium can be approximately considered homogeneous, there are a lot of settings where this assumption is not applicable a priori, in particular where a survey covers extensive areas and the subsurface may vary for instance in terms of sediment grain-size, material density/compaction, water content, target composition, voids occurrence among the other. Another relevant issue worth mentioning is the form of the recorded event. It is commonly assumed to be an ideal point-source hyperbolic diffraction, but frequently the scatterer is not point-like but has a flat or slightly inclined top part. Thereby the recorded event can appear like a misshapen hyperbola, composed by a planar/slightly inclined reflection and one or two lateral arms. This leads to an incorrect velocity estimation caused by difficulties in fitting a synthetic hyperbola. In this case, and also when a sharp lateral target termination is present in the subsurface (e.g., the edge of wall, interruption of stratigraphic layers caused by burials) [9], it should be more appropriate to use only one arm (left or right) of the fitting hyperbola to avoid introducing significant errors in the estimated velocity.

In this study, we propose the application to GPR data of the methodology based on high-coherency functionals computation. The computation of signal coherency (the similarity of the signal between adjacent traces) according to these functionals has long been proven to be effective in seismic exploration to estimate reflection event parameters in a common midpoint gather (CMP), even in cases of low S/N [10,11,12,13,14,15,16,17,18]. Moreover, this methodology offers a higher resolution of the estimated parameters with respect to the semblance functional [6] that is frequently used in various applications, including GPR velocity analysis [3]. The technique consists of the computation of high-resolution coherency maps by means of functionals that have been adapted from a pre-stack (source and receiver location some offset apart) to a post-stack (source and receiver coincident) velocity analysis. Applying this technique to GPR data may offer an additional processing tool capable of improving the correct location and focusing of subsurface objects in GPR profiles thanks to a more accurate velocity field estimation. We apply this technique first to a synthetic radar profile to test the effectiveness of the proposed method and evidence some possible issues potentially manifesting when real data are processed. Then we estimate the velocity field on a real data set composed of 41 parallel GPR profiles acquired in an archaeological site (Badia Pozzeveri, Lucca, Italy) and migrate the data in time with a Kichhoff algorithm. The archeological site was used because there is the possibility of ground-truthing the results of GPR data processing thanks to the evidence of excavations in the same area previously covered by the survey [9]. The migrated GPR lines were used to build some time-slices, depicting spatial coherencies of reflected energy consistent with archaeological targets [19,20,21,22,23,24,25].

2. Materials and Methods

2.1. Site Description and GPR Survey

The real GPR data of this study come from the archaeological site of Badia Pozzeveri (43°49′20″N, 18°38′39″E), located about 10 km south-east of the city of Lucca (north-west Italy). Here, an 11th-century church was believed to be larger than it is today, and buried structures were expected [26]. A GPR survey was conducted with a radar system of the IDS Georadar Company [27] equipped with a monostatic transmitter and receiver that operate at 600 MHz (nominal peak frequency). A common offset acquisition of 41 parallel survey lines 50 cm-spaced was performed covering an area of about 660 m². The in-line trace spacing of 1.2 cm was controlled by an odometer wheel. The configuration for acquisition provided 1024 samples in a time window of 100 ns (Figure 1). Previous GPR processing revealed the existence of walls buried at variable depths in the modern churchyard together with other reflective areas potentially related to anthropogenic features (e.g., tombs, burials, channels) [9]. These results found a good match with the excavation results undertaken after the GPR survey. Indeed, the excavation operation made it possible to discover the foundations and parts of the structures of the church, as well as several medieval and modern tomb features. Finally, in the area of the modern churchyard, the excavation revealed a channel for the drainage of rainwater, contained in a funnel-shaped section and originating near the entrance of the church.

2.2. High Coherency Functionals Rationale

The travel time equation of electromagnetic waves recorded at the surface from a point diffractor in a homogeneous medium is given by:

$$t=\sqrt{t_0{}^2+4{\left(x-x_0\right)}^2/V^2}$$

where t is travel time, x is the horizontal coordinate of the transmitting and receiving antennae, (x₀,z₀) the coordinates of the diffractor with the vertical axis z pointing downward, t₀ = 2 z₀/V the two-way travel time for x = x₀ and V the velocity of the electromagnetic wave in the medium.

A common way to estimate the velocity V is to use the semblance functional (Cs) that measures the ratio between the signal energy and the total energy in a given time window that follows a predefined hyperbolic trajectory within a user-defined aperture or number of summed traces:

$$C_s=\frac{1}{M}\frac{{{\displaystyle \sum }}_{t=t_0-T/2}^{t=t_0+T/2}{\left({{\displaystyle \sum }}_{i=1}^M{d}_i\left(t_0,V_s\right)\right)}^2}{{{\displaystyle \sum }}_{t=t_0-T/2}^{t=t_0+T/2}{{\displaystyle \sum }}_{i=1}^M{d}_i^2\left(t_0,V_s\right)}$$

where (t₀,Vs) are the parameters of the hyperbola, T is the window length in ns, M is the number of traces considered, d_i are the data of the i_th trace. The semblance is known to be a robust coherency measure against noise [6]. However, the semblance coherency functional is only able to compute low-resolution panels, limiting the precision and accuracy of the estimated subsurface velocity.

In signal processing, and in particular in seismic exploration, a wide variety of tools has been developed to detect a signal that is buried in the noise [13,15,16,17]. Such tools exploit the signal coherency among different traces gathered in common midpoint gathers (CMP). Even in the CMP case, the travel time equation is hyperbolic, and thereby the techniques that have been developed can easily be adapted to the GPR velocity estimation. Due to the high number and variety of existing tools found in literature, it is not possible to comment on them all in detail here. Therefore, we only give a brief description of those, in addition to the semblance, we have found more promising in the tests we performed on GPR data in terms of resolution, noise rejection and computational time.

The first one we consider is the Complex Matched Analysis (Ccm) functional proposed by Spagnolini et al. [28]. It makes use of an approximate estimation of the wavelet to obtain a more efficient rejection of the random noise and is able to better discriminate between interfering events. The Ccm is given by:

$$C_{cm}=\frac{1}{M}\frac{{{\displaystyle \sum }}_{t=t_0-T/2}^{t=t_0+T/2}{\left|{{\displaystyle \sum }}_{i=1}^M{D}_i\left(t_0,V_s\right)\right|}^2}{{{\displaystyle \sum }}_{t=t_0-T/2}^{t=t_0+T/2}{{\displaystyle \sum }}_{i=1}^M{\left|D_i\left(t_0,V_s\right)\right|}^2}$$

where (t₀,Vs), T, and M have the same meaning as before and D_i is the i_th trace data filtered in the Hilbert domain by an estimated or known wavelet. Its efficacy has been proved on seismic data also by Tognarelli et al. [17,29]; however, we expect a resolution similar to the estimated wavelet that is employed and, ultimately, to the semblance.

The second functional we consider is given by the Complex Matched Coherency Measure (Ccmcm) algorithm, originally described by Key and Smithson [11]. It exploits the high resolution that can be achieved using the eigenvalue analysis of the data covariance matrix computed in a short time window that encompasses the signal. Considering the first eigenvalue as due to the signal and the others as due to uncorrelated random noise, the S/N can be computed as:

$$S/N=\frac{{\lambda }_1-{{\displaystyle \sum }}_{i=2}^M\frac{{\lambda }_i}{M-1}}{{{\displaystyle \sum }}_{i=2}^M\frac{{\lambda }_i}{M-1}}$$

where λ_i are the eigenvalues of the data covariance matrix in decreasing order of magnitude and M is the number of traces used. The computed S/N re-scaled in the [0 1] interval gives the Ccmcm measure. The resolution that is achieved with this functional is very high. However, the combined use of a complex matched filter and the eigenstructure analysis of the data covariance matrix allows for better discrimination of interfering or crossing events [17].

For this reason, a third functional has been introduced, given by a combination of the Ccmcm and Ccm algorithms, trying to benefit from the positive properties of both. In other words, the Ccmcm has been integrated with the Ccm by Grandi et al. [14,15], thus combining the robustness inherited by the Ccm with the resolution capabilities of Ccmcm. Therefore, the new functional is given by:

$$C_{cmKS}=C_{cmcm}{C}_{cm}$$

Overall, we tested more than 15 functionals, built on the basis of the original description given in the cited references and also obtained through their combination. Our choice has represented the best trade-off between the required characteristics of a signal coherency estimator (i.e., high resolution and noise rejection capability) and allowed us to pick the velocity spectra on GPR data with good reliability.

3. Results

3.1. Velocity Analysis of the Synthetic Data

We tested the proposed methodology on a synthetic data set built by using the 2D forward modelling software GPRSIM 3.0 by Geophysical Archaeometry Laboratory [30] based on ray-tracing techniques [31]. Via a discrete grid, we built a ground model (Figure 2) composed by two regions of different relative permittivity (ε_r), resulting in a velocity of approximately 12.2 cm/ns and 6.7 cm/ns on the left and right sides, respectively. Inside each region, two iron objects of different sizes are placed at different depths. We run a GPR simulation according to the same parameters (antenna frequency, temporal window, sampling frequency in the spatial and temporal domains) adopted in the GPR survey composing the real dataset (see later in the text). Furthermore, and consistently with the constructive parameter of the employed GPR, we simulated a transmitter radiating GPR waves over a 170° angle (symmetrical half aperture of 85° with respect to the vertical). The impulse shape was set as Ricker wavelet, resulting in the computed radar profile shown in Figure 3. Note that in the asymptotic hyperbolic arms far away from the top of the events, some arrivals are not computed. This is due to the discretization of the grid or of the starting angles of the rays, and some arrivals are missed even if smaller discretizations than the ones used are set. However, we can consider this problem as a more severe test that mimics the absence of the signal in the velocity analysis because of missing traces. As it will be evident in the computation of the coherency spectra, the proposed methodology is solid enough to overcome this issue if it is not too severe.

The velocity analysis using the previously-described functionals has been carried out on the events of the GPR profile in Figure 3 with the parameters illustrated in

Table 1. Parameters used in the velocity analysis on synthetic data by the coherency functionals.

n	64	Number of Time Samples of the Velocity Analysis Window
win_ap	101	Aperture of the velocity analysis window (number of traces)
v	4–20	Velocity interval scanned by velocity analysis (cm/ns)
dv	0.1	Increment in the velocity scanned (cm/ns)
dt0	0.098	Time increment in the functional computation (ns)
noise	0.1	Standard deviation of the added white noise to avoid functional singularity
wavelet	600 MHz	Maximum of the Ricker wavelet frequency spectrum

. As examples, coherency maps computed for the positions A and D are displayed in Figure 4a,b, along with the picked events marked by a red circle. Each of these (t₀,Vs) pairs corresponds to the hyperbolic trajectory imposed on the GPR profile at the analysis location (hyperbola with magenta and green arms, Figure 3).

As it can be observed in Figure 4a,b, the resolution of Ccmcm and CcmKS functionals is much higher than the one given by the semblance and by Ccm, and consequently, the velocity picking is much easier and accurate. It is well worth noting that the velocity analysis described has been performed using only the portion of the data to the left of the apex (the data pertaining to the magenta hyperbola arm). The reason for this choice is that the full hyperbola does fit an event in the radargram if it originates from an object that can actually be considered a point diffractor. The object in A is sufficiently small to fulfill this requirement, as testified by the fitting of the green right arm also to the data. Instead, the dimensions of the object in D do not satisfy this point-like requirement and generate an event that differs from a classic hyperbola. Also, including in the velocity analysis the data in the correspondence of the green arm, we obtain a less accurate coherency map, and the velocity picking becomes inevitably more difficult or even impossible. For the event in A, the manual picking allows us to estimate a velocity of 11.8 cm/ns at 9.2 ns, and for the event in D a velocity of 6.6 cm/ns at 29.5 ns using the Ccmcm functional. These values are in good agreement with the background velocities given above (12.2 cm/ns and 6.7 cm/ns respectively).

The picked (t₀,Vs) pairs are used to time migrate the radar profile with a Kirchhoff algorithm [32]. Figure 5 displays in frames (a) and (b) the outcomes obtained using a constant velocity field built with the lower (6.7 cm/ns) and higher (12.2 cm/ns) velocity, respectively, and in frame (c) the outcomes when the input to migration is the varying velocity model built with the picked values (in A, B, C, and D). It is clear that only in the latter case do we have a satisfactory result from migration, and that in (a) and (b) only a portion of the radar profile is correctly migrated.

From this application of high-coherency functionals to synthetic data, it is evident that the migrated GPR profile can benefit from an accurate velocity analysis carried out both across time in a determined location and along the whole profile where hyperbolic diffractions are present.

3.2. Velocity Analysis on the Real Data

We performed the velocity analysis using the previously described functionals on data that have undergone the following processing steps to increment the S/N ratio and to reduce the reverberations: (1) background removal; (2) a band-pass filter in the range of 80-100-600-700 MHz to select only the frequencies related to the transmitted electromagnetic energy; (3) an automatic gain control on a time window of 10 ns to recover the low amplitudes of the signal at a later time; (4) a predictive deconvolution (prediction lag 2 ns, filter length 4 ns) to reduce the reverberations observed on the radargram, thus facilitating the velocity analysis; (5) a more selective band-pass filter in the range of 80-120-450-550 MHz.

Coherency maps of GPR data can only be significant in correspondence with the diffracted energy that has a hyperbolic trajectory. Due to the time required to select the locations and to compute the coherency maps, we performed the velocity analysis only on some selected GPR profiles (L-lines) of the data set where the diffraction events are more evident. Figure 6a,b show a portion of the L74 and L92 lines with imposed the estimated trajectories after the velocity picking on the coherency panels. The different colors of the left and right arms of the hyperbolas allow for easier identification of the data used in the analysis. Indeed, after some attempts, we chose to employ a 64-sample time window around the right arm for the L74 line, and an identical window around the left arm for the L92 line. In Figure 7a–d we show the coherency panels used to pick the (t₀,Vs) pairs and the values actually picked with red asterisks in positions E and F of L74, and G and H of L92. The velocity analysis in point E highlights the difficulties that may be found to match real diffraction patterns on GPR data. We selected the right arm (the use of the left one resulted in a too high velocity and an over-migrated section) and performed the analysis every 6 cm, eventually choosing the best velocity function on the basis of hyperbola fitting and migrated results. This made it possible to make a good compromise regarding the issue of possible uncertainties in velocity estimation due to the finite dimensions of the subsurface targets [8]. We made use of both the Ccmcm and the semblance functionals since the contribution of each of them was significant in the velocity analysis phase. Semblance was helpful for giving a rough indication of the regions in the (t₀,Vs) map where the signal concentrates, while Ccmcm was useful for its higher resolution, which allows us to select more precise values and discard artefacts introduced by the semblance algorithm.

In Figure 7 we also show the CcmKS functional that is generally less noisy than the Ccmcm due to the Ccm weighting, but at least in case E was less effective. A possible explanation can be given by looking at the Semblance. The sample-by-sample weights introduced by the Ccm in the CcmKS, which are similar to the semblance, filter out the coherency observed around 22 ns at 10 cm/ns that was detected by the Ccmcm functional. Note that having more than a single coherency measure available (the ones displayed plus the Ccm which is very similar to the Semblance and therefore not shown) helps to better discriminate the (t₀,Vs) pairs related to signal, leading to an improved velocity analysis. The parameters used for the real data case are indicated in

Table 2. Parameters used in the velocity analysis on real data by the different functionals.

n	64	Number of Time Samples of the Velocity Analysis Window
win_ap	81	Aperture of the velocity analysis window (number of traces)
v	4–20	Velocity interval scanned by velocity analysis (cm/ns)
dv	0.1	Increment in the velocity scanned (cm/ns)
dt0	0.098	Time increment in the functional computation (ns)
wavelet	400 MHz	Maximum of the Ricker wavelet frequency spectrum

. The velocity values obtained are used to build a velocity field and to migrate with a Kirchhoff algorithm the GPR profiles (L-lines) without imposing limits on the dip angle and on the available frequencies. This permits on the one hand to collapse all the hyperbola arms and on the other hand to preserve the maximum possible resolution allowed by the data. The migrated GPR profiles in Figure 6a,b are displayed in Figure 8a,b respectively, for an analysis of the results. As can be observed, along both profiles the diffracted energy is fairly collapsed in correspondence with the apex region (only point E appears slightly under-migrated), allowing us to better define the discontinuities in the subsurface.

This is more evident in the close-up of Figure 9a,b pertaining to the L92 line before and after migration, where the arrows indicate some interesting locations. It can be observed that after migration, the dimension of the objects illuminated by the electromagnetic energy is laterally well-delimited, thus giving a more realistic description of the size and shape of the buried targets. Moreover, thanks to the multi-channel (multi-trace) filtering process operated by the migration operator, random noise is reduced after migration.

3.3. The Effect of Varying Velocity Field on Time-slices

The migrated GPR lines were used to build time-slices (amplitude maps) parallel to the surface and cut at different times. It is well known that this kind of representation is largely used to correlate spatially coherent reflected energy (i.e., reflections) with targets at different depths. Testing whether the high-coherency functionals increase the readability of reflective areas in a time-slice representation is thereby crucial. To do this, we gridded and sliced three GPR datasets: (i) unmigrated data, (ii) data migrated with a constant velocity (10 cm/ns), and (iii) data migrated according to the varying velocity field reconstructed via the high-coherency functionals described above. The value of 10 cm/ns, used for the constant velocity migration, comes from the average of the picked velocity values on the coherency maps relative to several GPR profiles. In building the time-slices, we took care to keep slicing and gridding parameter constant (slice thickness: 4 ns; the square amplitude of amplitude values; cell dimension: 0.05 m; search area: circular; search radius: 1 m), as well as the interpolation statistical method (i.e., inverse weighted distance).

All these operations were performed with the GPRSlice software via Geophysical Archaeometry Laboratory [30]. For the sake of comparison, we illustrate only three time-slices in the temporal interval of 26–36 ns, centered at 28, 31, 34 ns and 4 ns thick (Figure 10). Instead of the same scale of reflection strength (amplitude), we used a relative normalization for each time-slice because this allows us to better imagine the subsurface targets and make the comparison more effective.

It is evident that the spatially-coherent reflective areas are clearer and of higher resolution in the time-slices built using the GPR lines migrated with a variable velocity field reconstructed via the high-resolution functionals (Figure 10). Moreover, the reflection intensity appears to have increased, and spot reflections, associated with local subsurface conditions (e.g., localized increases in density due to the presence of boulders/blocks concentrations or fragments of man-made remnants), are reduced. The time-slices built with GPR profiles migrated with a constant velocity (10 cm/ns) show content in amplitude higher than those built with un-migrated GPR profiles. However, in some sectors, it appears that the migration process has not favored the continuity of the reflections even if compared with those depicted in the un-migrated time-slices. This is an effect of the non-optimal velocity used in the migration that hampers the focusing of the reflected energy at the correct position.

As can be seen from the comparison with the excavation results (Figure 1), the linear-shaped reflections are consistent with the archaeological findings unearthed in the investigation area. There is a clear correspondence both with masonry structures composed of stones/bricks and cuts consisting in channels for rainwater drainage. Relative to the channel EW crossing the central part of the time-slices (especially visible in the interval 29–33 ns), the backscattered energy from the filling material (coarse gravel and boulders) has returned in the GPR profiles hyperbolic diffractive figures, although with an irregular form [9]. Note that we selected only three time-slices ranging between 26–36 ns, and some of the archeological structures may be not evident in Figure 10 because they are shallower or deeper. Indeed, the evident NS-oriented wall developed between 22 and 30 m of horizontal distance is at about 5–10 ns.

The test of high-resolution functionals on such a kind of data was quite demanding. This is because the structure of the buried walls is composed laterally by interlocked decimetric boulders (interconnected boulders and open voids) and internally by boulders dispersed in an abundant mortar. In GPR profiles, the tops of these structures correspond to a planar reflection, while the boulders placed on the sides and corners of the wall to arms of hyperbolas (i.e., half hyperbolas) (see the black arrow pointing to 22 ns and 1700 cm of distance in Figure 9a). Moreover, the boulders in the mortar-dominant core generated small hyperbolic diffractions with interfering arms [9]. While the latter is difficult to model and manage even with a variable velocity field, the approach based on high-resolution functionals has clearly made possible a more than satisfactory focusing of the principal diffraction hyperbolas, and therefore a better definition of the wall structures in the time-slices.

4. Conclusions

The aim of this paper was to verify whether the use of high-resolution coherency functionals could improve the quality in terms of S/N ratio and interpretability of GPR data by introducing a variable velocity field in the migration process.

After carrying out tests on synthetic data to schematically simulate the problem, assessing the types of functionals most suitable for GPR data analysis, a varying velocity field relative to a real dataset was calculated. This dataset was acquired in an archaeological area for which the position of the targets was already available because these have been unearthed after the GPR survey.

The best functionals turned out to be the Ccmcm and the CcmKS, because they are able to compute coherency maps characterized by high resolution and strong noise rejection that make it possible to better define the time-velocity pairs of the observed diffracted events. In the velocity analysis, particular attention was paid to checking the fitting between the reconstructed diffraction hyperbola arms and those generated by buried structures, which do not consist of ideal scattering points but targets with a certain lateral extension as in most real cases. In the selected data, we did not include the top planar part of the event, limiting the application of the method to the hyperbola arms, reducing the possibility of wrong velocity estimation. We observed that with the precisely-picked 2D velocity field migration algorithms perform better than in the case of constant or 1D velocity field, with satisfactory collapsing of the diffracted events and moving of the reflected energy to the correct position.

The varying velocity field was applied to all the GPR profiles composing the survey and covering the entire archaeological area. The time-slices constructed are of a higher quality than those obtained from non-migrated or migrated GPR profiles at a constant velocity. The improvements are inherent to resolution, continuity, and energy amplitude of linear reflective areas.

On the basis of our experience, we can affirm that the use of high-resolution functionals allows us to obtain migrated GPR profiles with a high grade of hyperbola focusing. On these profiles, the targets of interest are better-imaged and the interpretation can be carried out more reliably. A varying velocity field, even sparsely sampled like the one we have built, can allow for more accurate positioning in time of the objects of interest with respect to the constant or 1D velocity model routinely used, or when no data migration is applied. This could be relatively important when studying an archaeological site or a geological setting, cases in which a slightly incorrect positioning of targets/layers is acceptable. However, accurate reconstruction of the target geometry and positioning becomes crucial when it comes to studying other topics, e.g., those related to civil engineering and micro-structural settings, where even modest velocity variations can lead to reconstruction errors which may be quite relevant.

The main drawbacks we experienced were the time required by the functional computation (approximately 5 min on an i7 laptop for each location) and the time required to choose the locations for performing the velocity analysis. Concerning this last issue, one suggestion could be to perform the velocity analysis on a whole line (or the more promising portions of a line) and then inspect the velocity volume (time, velocity, inline distance) to select the best coherency maps on which to perform the velocity analysis. This can be done at a very high computational cost but can save some manpower time. An alternative option could be exploiting automatic tools of hyperbola detection to generate a dataset of locations on which to perform the high-resolution velocity analysis. By applying the functionals only to this limited number of points, selected along the GPR profiles, allows us to limit computational costs.

Notwithstanding these drawbacks, we believe that in cases where an accurate and detailed image of the subsurface is desired, this methodology is worth the effort.