1. Introduction
The proper recognition of foundation soil is essential for designing construction investments, both cubature and infrastructure. It increases safety and reduces additional costs at the realization stage and during the exploitation of engineering structures. One of the tools supporting ground recognition is a geophysical survey, which, among other things, recognizes the tectonic structure of the layers during sulfur extraction [
1], landslide monitoring in the infrastructure areas [
2] and tectonic movements [
3]. Several research projects have been conducted on the geometrization of changes to the ground-based on 2D geophysical profiles [
4,
5,
6,
7]. The possibility of contouring the boundaries of soil layers allows for a 3D analysis (i.e., determining the range and volume of the identified structures. This is particularly important for ground reinforcement design, stability assessment and soil replacement.
In the field of construction investment, however, the principles of subsoil recognition using geophysical methods are not standardized. In 2015, the Polish General Directorate for National Roads and Motorways (GDDKiA), together with the National Center for Research and Development (NCBiR), under the joint undertaking “Development of Road Innovation”, undertook the task of structuring and adapting the regulations to the current technological state. The outcome of the grant completed in 2018 is a project of “Guidelines for the performance of subsoil research for road construction” [
8,
9]. The scope of the geophysical methods was included very broadly, from the issues, through the selection of research methods and parameters, to their development and documentation. The recommended research methods included ground-penetrating radar (GPR), although priority was given to electrical resistivity tomography (ERT) due to the possible penetration range.
Geophysical surveys are a form of remote sensing, the accuracy and resolution of which are determined by the influence of the propagation environment on the electromagnetic signal. In the case of GPR, dependence on the value of the recorded signal amplitude on the wave propagation medium is simplified to the electromagnetic wave reflection coefficient Γ, according to the equation [
10]:
where
εr1,
εr2 are the relative electric permeability of the adjacent media.
Distinguishing the boundaries of the medium in GPR is possible when the power of the reflected signal
PΓ is described by the relation (Annan, 2000, 2001):
As Formula (2) shows, the greater the differentiation of the relative dielectric permeabilities of the adjacent layers, the greater the chance of unambiguously mapping the boundary between them. The relative dielectric permeability of air equals 1, and values for other media are at least three times higher. Such a parameter ratio provides good visualization of voids in natural or anthropogenic forms on radargrams, especially for media with low electrical conductivity (i.e., with a low electromagnetic wave propagation attenuation coefficient).
Due to clutter [
11], which interferes with the useful signal [
12], or wave attenuation, it is difficult to develop interpretation keys for GPR data even if based on calibration measurements.
Many research examples confirm the relatively high effectiveness of the GPR method in the detection of underground utility networks [
13,
14] and the localization and diagnosis of shallow crypts and tunnels [
15,
16,
17]. However, some differences between anthropogenic and natural forms should be noted. The examples reported in the mentioned publications presented results obtained on anthropogenic objects with regular shapes. Moreover, the calibration of the results of geophysical surveys carried out on anthropogenic objects can be performed based on information such as the typical geometry of the structure, maps and projects (underground utility networks, road layers, stratified walls, reinforcement elements, flood embankments). The deposition depth of these objects is not large, varying between 2 and 4 m. Similar values were found for typical crypts, cellars and tunnels in both modern and historical buildings.
Among the studies of GPR application to the recognition of natural phenomena, many have referred to the detection of karst formations [
18,
19,
20]. These occur in carbonate rocks for which the relative electrical permeability is in the range of 6–11 [
21], which affects the reflection coefficient at a level of 0.42–0.54, and the power of the reflected signal in the range of 0.18–0.29. Therefore, there are favorable conditions for applying GPR at a media boundary like air–rock (especially of high compactness), mainly to detect shallow natural and anthropogenic voids [
22,
23]. However, for anthropogenic objects such as irregular or unidentified excavations and post-mining voids [
24,
25], and natural objects, like lithological layers [
26] or karst phenomena in the subsoil [
26,
27], the necessary calibration material is usually information from boreholes. Occasionally, calibration was performed on exposed geological structures in the study area, but this approach also had significant limitations because of edge effects [
21].
In areas with varied geological structures, the confrontation of GPR and ERT results is also a subject of research. As in the aforementioned “Guidelines”, some authors underlined the significance and superiority of ERT in the hierarchy of geophysical methods for ground recognition [
26,
28]. However, the advantage of GPR is that it gives a result with better resolution and efficiently determines the boundaries of materials having large variations in the relative dielectric permittivity of low conductivity, which was emphasized in the publication [
29]. In addition, research [
30] has noted that GPR was more successful in identifying previously recognized karst phenomena (
Figure 1).
According to Ortyl [
30], a wide discussion on the application and effectiveness of geophysical methods was carried out. In particular, he investigated the GPR method for recognizing and mapping changes in the structure of the subsoil. Furthermore, an algorithm for converting GPR data into a point cloud file based on geopositioned radargrams was presented. This concept is described as underground radar scanning (UgRS) (
Figure 2). Such an approach to GPR data processing eliminates the need for direct access to voids or caves to measure their geometrization and quantification by laser scanning [
26,
28]. The paper also offered an extensive evaluation of the effectiveness of selected filtering procedures and wave speed calibrations based on terrestrial laser scanning (TLS). Analyses were performed based on synthetic data generated by the Finite Differential Time Domain (FDTD) method and real data.
The following paper is the result of a practical implementation of GPR measurements and data processing to recognize karst formations. During the construction of a national road investment with highway parameters, karst phenomena were identified that had not been indicated in the geological documentation. GPR surveys were carried out on a 7 km section of the newly built road, giving separated areas with different geological structures.
Section 2 presents the methodology of data acquisition and postprocessing, in particular the filtering algorithm and the integration of GPR with geodetic data. As a representation of the digital terrain model, the point cloud was acquired using low-range photogrammetry and an unmanned aerial vehicle (UAV). The strong GPR signal responses may have indicated an anomaly caused by a void, rock outcrops, or variations in soil moisture; therefore, post-processing with calibration was required to interpret the GPR data correctly. This paper will then analyze and validate the GPR results based on control trenches with very wide cross-sections. This method allows the exact calibration of the GPR signal response visualization with the actual layout of the ground layers.
2. Materials and Methods
2.1. Description of Research Area
The entire research area was located in the southwestern part of the Lublin Upland, within the mesoregion of Urzędów Heights, close to the border with Western Roztocze and the Biłgoraj Plain. It is a plateau with levelled hilltops cut by river valleys and gullies with steep slopes, which are characteristic of loess areas [
31]. The only exposure of Cretaceous and Upper Jurassic formations is located within the Rachowska and Gościeradowa anticlines on a relatively small area. Upper Cretaceous rock predominates, favoring the maintenance of the Tertiary surface plan. The presence of tertiary rocks determines the landscape of the southern part of the region. Especially in the southwestern zone, it is revealed by the stepwise arrangement of monadnock morphological elements [
32] (
Figure 3).
The research area covered the national expressway using the construction parameters of the motorway, which was under construction. GPR and geodetic surveys were carried out on five separate polygons, having a total length of the study route axis of about 2 km distributed over a 7 km length (
Figure 4).
2.2. Methodology of Data Acquisition
Field surveys for polygons 1–2 and 4–5 were conducted using Mala RAMAC CU II GPR with a shielded antenna of 250 MHz, which provided an assumed depth of effective penetration up to 5 m. The choice of the antenna frequency was a compromise between the range and resolution of measurements. It enabled ground penetration to capture possible boundaries of native rock, changes in the type of substrate and the presence of karst phenomena. The time window was assumed to be 200 ns, the lowest theoretical speed of electromagnetic wave propagation in limestone equal to 0.09 m/ns, and reached the potential range of penetration up to 8 m. For the highest theoretical speed in limestone (0.12 m/ns), the range increased to 12 m, and echogram traces were recorded at 0.1 m intervals. The signal sampling frequency was 2526 MHz, and the number of signal stacks was assumed to be 32. The Leica TCRA 1002plus robotic tachymeter was used for real-time geopositioning. Measurements were carried out with reference to Poland’s national coordinate system 2000 (CS2000) zone 7, and the Kronstadt’86 height reference frame.
Preliminary analysis of the results from Mala GPR confirmed suitable conditions for conducting GPR measurements to identify the karst phenomena. Therefore, in polygon 3, the Leica DS2000 GPR was used. This dual-channel unit allowed the simultaneous acquisition of radargrams from two antennas at a frequency of 700 and 250 MHz. The use of the higher frequency, 700 MHz, ensured a higher resolution image of the karst phenomena. A similar approach using Mala GPR to acquire data was performed, especially regarding the assumed values of the measurement parameters. Thus, the time window was equal to 200 ns, the signal transmission trigger was based on a distance interval of 0.1 m, and a signal sampling frequency of 2564 MHz. The maximum number of signal stacks for this GPR unit was assumed to be equal to 16. In this solution, the positioning of the GPR antenna was conducted using GNSS satellite measurements in the RTK–RTN mode (referenced to the HxGN SmartNet base stations) in the global system WGS84. As part of the data processing and for further integration, all results were transformed to CS2000 zone 7 and the Kronstadt’86 height reference system.
Georeferencing results of GPR data allowed for their integration with the Digital Terrain Model in the form of a point cloud obtained from low-altitude photogrammetric measurements conducted using a UAV. A typical grid of profiling of the research polygons is presented below (
Figure 5). The line spacing between adjacent profiles was 3 m, which resulted in an average number of parallel profiles equal to 11. The transverse and diagonal profiles can also be seen in the figure below. Their implementation aimed to assess the geometric compatibility and content repeatability of GPR results at the intersection of profiles with different driving-up directions.
2.3. Methodology of Data Post-Processing
In the case of the GPR method, the extraction of information necessary for correct data analysis required advanced processing algorithms. In the first step, the GPR data was filtered using ReflexW software by Sandmeier & Software. The basic sequence of procedures included:
Moving the start time to zero as the point of first contact between the electromagnetic wave and the ground;
Shifting the average signal value to zero—subtract-DC-shift;
The removal of the low-frequency signal—subtract-mean (dewow);
The removal of horizontal reflections by the subtraction of average trace—background removal;
The removal of the bottom part of the radargram without significant reflections—time cut.
The selection of procedures for the second part of the filtering was not accidental but based on years of experience and others’ research [
30]. Finally, the following filters were combined in the following processing step:
Subtracting the average filter.
Energy Decay as gaining filter.
Bandpass butterworth filter for band range 50–400 MHz.
Average xy-filter 3 × 3 for smoothing.
The filtering procedure resulted in better visibility in the radargrams of regions with variable response characteristics to the electromagnetic field in the medium. However, the analysis of a single vertical B-scan image allowed identifying anomalies or zones of variable parameters but only within the range of a single measurement profile. The full picture of the magnitude of the phenomenon and the variability of the ground nature was observable due to spatial visualization and the integration of GPR and geodetic data.
Consequently, using the accurate georeferencing of results, GPR radargrams were also converted into the form of an underground radar scanning (UgRS) point cloud. The EM velocity value equals 0.1 m/ns has been used to convert the TWT time into real depth. In contrast to laser scanning measurements, UgRS is based on a recording of the emitted electromagnetic wave amplitude, which carries information about the intensity of the signal reflected from the structures in the ground. The conversion of the time-to-depth value performed based on the propagation wave velocity, allowed us to calculate the full spatial position of the reflection point. The accuracy of the point position was adequate for the measurement method and its geopositioning. GPR data were converted to the UgRS point cloud file based on an original script prepared in MATLAB and ReflexW software. The algorithm consists of two modules that perform the following tasks:
The conversion of field data into a file format of dedicated software and possible simultaneous coordinates correction [
33].
Processing and saving the filtered GPR radargram in a point cloud file format.
The first procedure concerns Mala’s GPR device and is performed based on data acquired during field measurement, conducted with GPR-GNSS or GPR-robotic total station set. This module is responsible for converting coordinate files, which regards the GPR antenna position as having been recorded using GroundVision software (*.lcg or *.cor file format), to ReflexW software format (*.utm file format). Furthermore, the application allowed the smoothing of the z-coordinate values. This option was particularly important in the case of a positioning solution based on satellite measurements when the precision of the adjacent results is variable. For that purpose, the algorithm used a subtracting average with the possibility of selecting a filter window size. There was also the possibility of a coordinate geometric correction implementation related to the slope of the GPR profile.
The second procedure was based on processed radargrams generated in ReflexW. It allowed the filtered radargram with georeferencing to be written to an ASCII file. The first line of the result file found the total number of points, followed by a stream of xyz coordinates, intensity values (laser beam reflection intensity) and colors in RGB for all points. One advantage of this procedure was the ability to choose whether the traces would be presented as vertical or tilted, according to the surface normal. The slope value was calculated based on the z-values from coordinates in *.utm file.
The final result of the process was the ability to handle GPR data as a point cloud and enable their integration with laser-scanning or vector and raster data. The figure below (
Figure 6) presents a visualization of the UgRS point cloud with a digital terrain model as a point cloud from one of the research polygons.
The next crucial step in data processing, which contributed to the increased reliability of the ground variability assessment, was to conduct GPR data calibration. It consisted in comparing selected radargram fragments with corresponding fragments from visible soil structures made possible by, e.g., boreholes and trenches. For this purpose, on the basis of the spatial overview of GPR and DTM data, the accurate locations for a few control trenches were determined on each surface (
Figure 7). They were selected so that the confrontation of the real situation with the radargram view would correctly identify potential karst phenomena in the entire research area.
4. Discussion
The purpose of this study was to present the methodology of conducting and analyzing GPR surveys on selected parts of a newly built expressway to determine the interpretation keys for identifying the location of karst forms such as sinkholes or karst funnels.
The work was conducted under specific conditions while the construction of the project was in progress. The measurement area mostly covered the road section levelled to the road alignment. The distribution and intensity of the carbonate rocks and overburdened materials were evident in the exposed soil regions. Carbonate outcrops dominated; hence, in accordance with the mentioned “Guidelines”, the ERT method was excluded in favor of GPR. Furthermore, as presented in
Figure 1, GPR also allowed the identification of subsurface voids, which would have been challenging for the ERT method.
In the examples, despite the use of an antenna frequency of 250 MHz, which was a compromise between the measurement resolution and the expected range of effective penetration, the layer boundaries with large differences in relative electric permeability were very clearly mapped on radargrams. The analyzed fragments showed that the soil variability, even including the variability of the carbonate rocks, themselves (compact/fine/loose/aggregate), correlated with the distribution of signal variability in the GPR radargrams. In places where strong reflections appeared, the rock was compact. The GPR response lost strength in the areas of loosened or weathered rock. Completely weak signals indicated a correlation with areas of overburden with clay or clay–sandy material, including karst funnels refilled with clay material.
However, it should be highlighted that an accurate determination of the nature of the observed soil variability was not possible because only GPR data were available. The remote sensing nature of the method is naturally associated with a significant uncertainty of identification and forces the development of interpretation keys through direct measurements in trenches. Strong reflections on radargrams may have indicated a void as well as a rock with crumbly nature. Hence, it was necessary to compare the GPR material with the exposed ground fragments. To categorize the type of soil over large areas even roughly, it was necessary to assume the homogeneity of the identified anomalies throughout the study area. Additionally, the physical principles of the GPR method made some phenomena impossible to identify. This was confirmed by the example of the ineffectiveness of identifying a few fractures in trenches up to a depth of 3 m. Despite using an antenna with a high frequency of 700 MHz, the cracks were too small for their boundaries to be determined on the radargram among the signal responses from crumbly heterogeneous rock.
An important aspect supporting ground recognition in GPR data analysis is the novel approach to their visualization.
For a long time, the integration of GPR and geodetic data focused on the method and geolocation of radargrams and their assessment [
34,
35,
36]. The ability to precisely determine the position of GPR data allowed for the implementation of further steps, i.e., combining the GPR data with the results of other sensors and derived products, especially rasters, such as orthophoto or DTM. In many cases, GPR data, through interpolation algorithms, were transformed to cube form, which allowed for the generation of horizontal cross-sections–C scan. This approach is effective for anthropogenic objects that exhibit a certain regularity of form in a spatial system [
37,
38]. However, it limits the analysis to a simultaneous view of a single cross-section at a predetermined depth or inclined plane, which significantly affects the possibility of full spatial interpretation of the results. The approach presented in our work is based on GPR data visualization as an UgRS point cloud. It avoids the problem of limited possibilities of correct interpolation of irregular shapes of natural objects. The data is still presented spatially, but without the redundancy of artificial interpolation data.
Furthermore, due to the accurate geolocation and radargram conversions using the author’s algorithm, it was possible to spatially compile the GPR and geodetic data in software specific to 3D object analysis. The spatial combination of the results of these remote sensing methods in one visualization made it possible to identify the correlation between the ground surface, visible on the DTM data in the point cloud form, and the recorded GPR data. In particular, a strong relationship was identified between areas of clay–sandy soil or visible limestone outcrops and the distribution of GPR signal amplitude on radargrams, which showed strong attenuation.
The essence of the joint data presentation in the form of a point cloud and supporting the analysis of the subsurface situation through the premises observed on the surface data strengthens the quality and interpretational credibility, which is emphasized in the publication [
39]. An additional value of the spatial analysis of UgRS data was the possibility of assessing the consistency of radargram content distribution on intersecting profiles, as in the case of Polygon 1. High coherence in the arrangement of signal-response amplitudes proved both the high relative accuracy of the GPR data location as well as their non-random local character. Thus, it increased the reliability of the method as stable and repeatable. The integration of DTM and UgRS point clouds also enabled the efficient and unambiguous selection of optimal sites for control trenches. The analysis of the material focused on those fragments of the GPR data that represented areas of an anomaly with the compact and strong amplitude of the back signal. However, while the nature of the GPR signal response in some fragments of the study areas did not clearly indicate the presence of large voids, the amplitude analysis in the surface approach indicated a small void in the crack networks of the rock material, which significantly affected ground load-bearing capacity. This was a serious indication for the execution of a control trench in this place.
Unfortunately, no soil samples were collected during the work, and thus it was not possible to determine the variability of ground structure and composition in the road laboratory. Ground analysis performed in this way could lead to enhanced calibration of radargrams and the determination of precise interpretation keys in a numerical approach. This would allow for the automation or semi-automation of the classification of the nature of soil changes, including even a quantitative assessment of anomaly distribution in the surveyed area. This assessment for Polygon 3 is as follows:
Characteristics revealed in the trenches 3A, 3B covered 20% of the surveyed area
In trench 3C it was about 9%
In trench 3D, about 47%
In trenches 3G and 3H, 10%.
An embankment covers about 7% of the area.
Minor local areas of strong attenuation cover about 6% of the area.
It was also concluded that with a large base of accurate interpretation keys, automation of ground classification would reduce the necessity of assuming a homogeneous nature of change within a whole study area.
The necessity of conducting complementary research resulting from the discrepancy between the existing geological materials and the actual state in the field only highlights the importance of correct soil recognition, as early as the design stage, both for infrastructural and cubature investment. As our research showed, it resulted not only in a safety issue and reduced additional costs, but also a much lower possibility of effective identification when construction is at an advanced level. Where GPR measurement was carried out on a road surface under which slag material was located, only a strong attenuation of the signal was visible on the radargram. Thus, it was concluded that the implementation of non-destructive geophysical surveys should occur at the latest before the formation of such a layer. The demonstration of this relation by comparison with the respective radargrams (
Figure 32) had a broader meaning. If the survey aim were to conduct subsoil recognition for an existing road modernization, and there were indications that the supporting substructure layer included slag, limitations can be expected from the very beginning.
5. Conclusions
The “Guidelines” systematize and introduce a clear methodology for conducting geophysical surveys and classifying the cases of their application.
GPR research should be performed before conducting invasive tests so that they can be targeted at areas of identified anomalies or variations in signal behavior. It was possible to conduct GPR measurements during investment construction to obtain additional information after removing the layers to the level of the road alignment but before laying the structure layers with slag elements. Carrying out a survey on a completed road structure with the presence of such a layer strongly reduced the effectiveness of the GPR method. This should be carefully considered when planning the work. In some cases, measurements should be designed outside the roadbed if possible due to site accessibility.
Requiring the geopositioning of GPR measurements had positive consequences in processing and interpreting the GPR data. First of all, it enabled the implementation of the proposed methodology by converting radargrams into a UgRS point cloud. This was an original approach, which as presented in the publication in combination with DTM data, gave an improved interpretation and understanding of GPR–substrate relations.
This practical example, which implemented GPR, demonstrated that we can all conduct a full and efficient integration of remote sensing data from various sources. The digitization of the construction process seems to be an irreversible trend. In addition, the example confirmed the high effectiveness and usefulness of using the GPR method to identify the variability of the substrate, especially in areas of carbonate rocks where karst phenomena and forms may occur.