## 1 Introduction  
The ground-penetrating radar (GPR) method provides high-resolution, non-destructive imaging of shallow subsurface and, among other applications, has proved to be a valuable tool for the delineation of the internal geometry of fractures, discontinuities, and voids in rock bodies [1, 2, 3, 4, 5, 6, 7]. This method has been applied by several researchers over limestone, a common base-rock, which is characterized by the generation of cavities and fractures due to inner water movement resulting in a calcium carbonate dissolution process. Even though several techniques have been used to improve the detectability and detailed imaging of the GPR method of sections, such as texture attributes [8] or signal polarization studies [9], diffraction imaging is the most common technique [2]. [PERSON] et al. and [PERSON] et al. [10, 11] agree that diffracted energy needs to be focused to arrive at laterally resolved GPR images. Diffraction focusing collapses the diffractions to their apices, enhancing, in this way, the coherency and linearity of the fracture planes. Parallel GPR study lines can provide 3D GPR data imaging, which, if not densely sampled in the off-line direction (assuming a minimum requirement of a quarter-wavelength spatial sampling), is spatially aliased and requires spatial interpolation for an efficient 3D migration process. An extensive study on this drawback was carried out by [PERSON] et al. [10]. Decimation of spatial sampling in the in-line direction to avoid aliasing, or high-cut filtering as is commonly used in seismic, also decimates the diffraction information, while on the contrary, dense parallel study lines increase the acquisition time and may lead to a problematic 3D migration scheme due to the large amount of data. These are the main reasons why researchers and engineers mainly use grids with a study line spacing of 0.25 m or, in exceptional cases, down to 0.1 m, mainly when using multi-antenna equipment. This has led to the commonly employed 2D migration schemes, free of diffraction focusing in the off-line direction, which have been proven to provide a weighted result between acquisition and processing time on the one hand and efficient 3D imaging on the other. What is common in both 2D and 3D migration processes and is unfortunately taken for granted in all of the above references is the need for an efficiently detailed velocity model of the subsurface. Velocity models for sufficient migration outputs may be extracted by multi-offset GPR data, which require time-consuming acquisition and are also characterized by the use of the conventional processing approach used in seismic reflection methods [12; 13]. As the most commonly used data acquisition method is the common-offset technique, adequate velocity models cannot be extracted, even though several approaches such as average velocity values based on ground truth [14; 15] or hyperbola fitting over diffractions [16] have been reported. An efficient migration velocity model for diffraction focusing is the multipath summation process, which involves the stacking of constant velocity migrated sections [17; 18]. Still, this approach requires subjective picking in its last stages. Multipath summation was negatively evaluated as a direct imaging method by [PERSON] and [PERSON] [17], because of the increased blurriness of the constant velocity-migrated stacked images. [PERSON] et al. [19] presented a diffraction focusing strategy to circumvent the picking subjectivity and degraded GPR sections due to both frequency-dependent attenuation effects and the decrease of the dominant frequency, consisting, in the summation, of constant-velocity migrated images. This 2D imaging strategy does not require to explicitly build EM migration velocity models and involves, except for conventional processing steps, the weighted stacking of sequentially constant velocity migrated sections over a range of EM velocity values and time-varying spectral whitening, with the aim to enhance both the lateral and vertical resolution.