Underground Object Classification for Urban Roads Using Instantaneous Phase Analysis of Ground-Penetrating Radar (GPR) Data

Abstract

Ground-penetrating radar (GPR) has been widely used to detect subsurface objects, such as hidden cavities, buried pipes, and manholes, owing to its noncontact sensing, rapid scanning, and deeply penetrating remote-sensing capabilities. Currently, GPR data interpretation depends heavily on the experience of well-trained experts because different types of underground objects often generate similar GPR reflection features. Moreover, reflection visualizations that were obtained from field GPR data for urban roads are often weak and noisy. This study proposes a novel instantaneous phase analysis technique to address these issues. The proposed technique aims to enhance the visibility of underground objects and provide objective criteria for GPR data interpretation so that the objects can be automatically classified without expert intervention. The feasibility of the proposed technique is validated both numerically and experimentally. The field test utilizes rarely available GPR data for urban roads in Seoul, South Korea and demonstrates that the technique allows for successful visualization and classification of three different types of underground objects.

1. Introduction

Sinkholes are one of the most severe threats to urban roads. Sudden ground collapse can lead to road closure, economic loss, and human fatalities and injuries. Several major sinkhole incidents, including in urban roads, have been reported in recent years, e.g., in Hangzhou, China (April 2016) [1], Fukuoka, Japan (November 2016) [2], and Cheltenham township, PA, USA (January 2017) [3]. Extraction of groundwater, changes in water drainage patterns, and water main breakages are considered to be the main causes of the underground cavities that can potentially develop into urban sinkholes [4,5]. As sinkholes often appear without any forewarning, the demand for their early detection in urban areas is increasing.

Nondestructive sensing via ground-penetrating radar (GPR) is gaining popularity for underground visualization, because it offers noncontact testing, fast scanning speeds, and deep penetration [6]. GPR emits high-frequency electromagnetic waves into the subsurface and analyzes the reflected electromagnetic waves. As reflection only occurs if there is an abrupt change in the electromagnetic characteristics of the medium, underground structures and their boundaries can be identified by analyzing the received signal. The main applications of GPR have been in the archeological and geological fields, where it is used for mapping subsurface features [7,8,9]. It has been also applied for other fields [10], including demining [11], identifying veins of the mineral [12], monitoring of mines and tunnels [13], and forensic applications [14]. Its use has recently been extended to detecting underground cavities and potential sinkholes [15,16,17].

Currently, the identification and classification of underground objects with GPR data heavily rely on the time-consuming and labor-intensive work of well-trained experts. In urban areas, an expert may be able to analyze no more than a few kilometers of three-dimensional GPR (3D GPR) data in a week. Moreover, the interpretation of these data depends on the subjective judgment of the expert, meaning that the results may be unreliable at times. GPR signals reflected from the pavement layer of urban roads typically dominate the reflection profile because of the large permittivity difference between the pavement interface and air. Underground objects cause relatively weak and noisy reflected signals, making data interpretation more challenging. Although the dominant signals can be partially eliminated by subtracting a reference signal from the GPR signal (known as the subtraction or background removal method [18,19]), this method is highly susceptible to noise and small signal disturbances in the time domain, e.g., variations in pavement thickness or unexpected arbitrary boundary conditions in the pavement. Time-varying gains have been used to enhance weak reflection signals [6,20], but this requires the careful selection of the gain values on the part of the user. A further technical challenge in GPR data interpretation is that similar GPR features are often generated by different underground objects. Even if an underground object is identified, it is challenging to classify it. In particular, the reflected signal from an underground cavity is often indistinguishable from the signals from other underground objects, such as pipes, manholes, or gravel. Without cumbersome underground boreholing, it can be difficult to confirm the presence of a cavity.

This study addresses these issues by proposing an instantaneous phase analysis technique for underground object classification. This technique aims to enhance the visibility of underground objects and provide objective criteria for GPR data interpretation, so that the objects can be automatically classified without expert intervention. Background filtering is first applied to 3D GPR data to enhance the visibility of objects beneath urban roads and hence reduce the need for human intervention in data analysis. The background signals, including dominant reflections from the pavements, are filtered out using a basis pursuit approach, and the underground objects are visualized by reconstructing the filtered 3D GPR data. An instantaneous phase analysis technique is then applied to classify the visualized underground objects. Reflections from cavities are in-phase with the direct GPR waves, whereas reflections from other underground objects are out-of-phase, meaning that underground objects can be classified while using phase analysis. Although the proposed technique uses a simple electromagnetic wave propagation model, numerical modeling and experimental validation with rare field-test data from urban roads in Seoul, South Korea confirm that it is effective for underground object identification.

This paper is organized, as follows. Section 2 details the working principles of GPR and explains the proposed instantaneous phase analysis technique. In Section 3, the technique is numerically validated using models of various types of underground objects, such as an underground cavity and a buried pipe. Section 4 then examines the feasibility of the proposed technique with experimental field data obtained from urban roads in Seoul, South Korea. The paper concludes with a brief summary and discussion of its application in Section 5.

2. Underground Object Classification Using GPR Data

This section sets out the proposed instantaneous phase analysis method for underground object classification. The working principles of GPR and the technical limitations of conventional GPR data analysis techniques are first briefly introduced. A novel underground object classification technique to automatically classify the underground objects is then proposed. This involves enhancing the visibility of underground GPR signals through the application of a basis pursuit filter and then using instantaneous phase analysis to classify the visualized object.

2.1. Working Principles of GPR

Figure 1 schematically illustrates a typical GPR setup for underground object visualization. The GPR system consists of an antenna and a data acquisition system (upper part of Figure 1). In the antenna system, a transmitter (T) emits high frequency electromagnetic wave pulses and a receiver (R) detects the reflected electromagnetic waves. Electromagnetic wave pulse excitation is commonly achieved while using stepped-frequency technology by sweeping the desired frequency band [21]. Any change in the electromagnetic characteristics of the medium, such as its permittivity and permeability, causes the electromagnetic waves to be reflected [22]. The incident electromagnetic field of the waves at the reflection location can be represented by

$${\mathbf{E}}_{\mathrm{i}}\left(z,t\right)=E_{\mathrm{i}}{e}^{\mathrm{j}\left(k_{1}z-\omega t\right)}\widehat{\mathbf{x}},$$

(1)

$${\mathbf{B}}_{\mathrm{i}}\left(z,t\right)=\frac{E_{\mathrm{i}}}{v_1}{e}^{\mathrm{j}\left(k_{1}z-\omega t\right)}\widehat{\mathbf{y}},$$

(2)

where ${\mathbf{E}}_{\mathrm{i}}\left(z,t\right)$ and ${\mathbf{B}}_{\mathrm{i}}\left(z,t\right)$ indicate the incident electric and magnetic fields at a specific depth $z$ and point in time $\mathrm{t}$, and $E_{\mathrm{i}}$ is the amplitude of the incident electric field. $k_1$, $\omega$, and $v_1$ refer to the wavenumber, frequency, and velocity of the electromagnetic waves in the corresponding medium, respectively. The reflected electromagnetic field can be expressed as

$${\mathbf{E}}_{\mathrm{r}}\left(z,t\right)=E_{\mathrm{r}}{e}^{\mathrm{j}\left(-k_{1}z-\omega t\right)}\widehat{\mathbf{x}},$$

(3)

$${\mathbf{B}}_{\mathrm{r}}\left(z,t\right)=-\frac{E_{\mathrm{r}}}{v_1}{e}^{\mathrm{j}\left(-k_{1}z-\omega t\right)}\widehat{\mathbf{y}}.$$

(4)

here, $E_{\mathrm{r}}$ represents the amplitude of the reflected electric field:

$$E_{\mathrm{r}}=\left(\frac{n_1-n_2}{n_1+n_2}\right)E_{\mathrm{i}}=\frac{\sqrt{{\epsilon }_1{\mu }_1}-\sqrt{{\epsilon }_2{\mu }_2}}{\sqrt{{\epsilon }_1{\mu }_1}+\sqrt{{\epsilon }_2{\mu }_2}}{E}_{\mathrm{i}},$$

(5)

where $n_1$ and $n_2$ are the refractive indices of each medium. ${\epsilon }_1$ and ${\mu }_1$, and ${\epsilon }_2$ and ${\mu }_2$ refer to the relative permittivity and relative permeability of the initial medium (1) and the new medium (2), respectively. These reflected electromagnetic waves are then collected by the data acquisition system. The collected signal contains multiple pulses corresponding to each reflection. The amplitude of the reflected electromagnetic field is larger where there is a greater difference in electromagnetic properties. Note that the above explained electromagnetic wave propagation is a simplified model and it does not perfectly match with the actual wave propagation on the soil. However, even though the proposed technique uses a simple electromagnetic wave propagation model, its superiority for underground object classification will be verified while using numerical and experimental methods in the following sections.

The reflection data are next analyzed for underground object detection. When GPR scans along the region of interest, three different data representations, A-, B-, and C-scans, are typically obtained. These are illustrated in the lower part of Figure 1. The A-scan provides punctual time domain information (the z-axis in Figure 1) at a certain spatial point. The B-scan image, representing the x–z plane, is constructed by stacking multiple A-scan data along the scanning direction (the x-axis in Figure 1). If there is an abrupt change in permittivity due to an underground object in a specific area, a parabola typically appears in the B-scan image. This phenomenon occurs because the transmitter and receiver are not located at the exact same spatial location. This parabola feature is often used to detect underground objects. The C-scan image is obtained in the x–y domain by stacking B-scan images that were measured by multiple antennae. C-scan images can be obtained along the z-axis to physically represent a certain depth. C-scan images revealing cross-sectional images of objects are also useful for detection purposes.

2.2. GPR Signal Analysis for Underground Object Classification

Figure 2 provides an overview of the proposed technique for underground object classification, which proceeds according to four steps: (1) GPR data collection, (2) background filtering to enhance underground object visualization, (3) instantaneous phase analysis for underground object classification, and (4) decision making.

2.2.1. GPR Data Collection

Time domain GPR signals are collected from each spatial scanning location along the region of interest. A two-dimensional GPR (2D GPR) image (B-scan image) is obtained by representing the collected data in the time and spatial domains. Electromagnetic wave attenuation along the depth direction is compensated by applying gain factors. Electromagnetic waves are typically attenuated as they propagate through a medium [23]. This phenomenon can be represented as [24]

$$E_{\mathrm{r},\mathrm{measured}}=E_{\mathrm{r}}{e}^{-\alpha l},$$

(6)

where $E_{\mathrm{r},\mathrm{measured}}$ is the measured amplitude of the reflected electric field, $l$ is the distance propagated by the electromagnetic waves, and $\alpha$ is the attenuation constant given by

$$\alpha =\omega \sqrt{\frac{\mu \epsilon }{2}\left[\sqrt{1+{\left(\frac{\sigma }{\omega \epsilon }\right)}^2}-1\right]}.$$

(7)

here, $\sigma$ refers to the charge density of the medium. The attenuation is compensated by applying the following gain factor to the measured signal ${\mathbf{E}}_0$:

$$E\left(t\right)=e^{\alpha v_{1}t}{E}_0\left(t\right)=e^{\alpha l}{E}_0\left(\frac{l}{v_1}\right).$$

(8)

2.2.2. Background Filtering to Enhance Underground Object Visualization

To extract the reflected signals from underground objects from the GPR data, it is necessary to remove the undesired background signals, including (1) pavement surface reflections, (2) reflections from underground layers, and (3) background noise. Interactions between the electromagnetic waves and underground objects can be identified more easily by representing the measured time domain signal $\mathbf{s} \left(T\times 1\mathrm{vector}\right)$ as a weighted linear combination of bases via the following transformation [25,26]:

$$\begin{gathered} \mathbf{s}=\mathbf{D}\alpha , \\ \mathbf{D}=\left\{{\mathbf{d}}_1, {\mathbf{d}}_2, \dots , {\mathbf{d}}_L\right\}, \end{gathered}$$

(9)

where $\mathsf{\alpha }$ is a representation of $\mathbf{s}$ in the transformed domain with the dimension $L\times 1$. D is a $T\times L$ dictionary matrix and consists of ${\mathbf{d}}_i$ bases (i = 1, …, L). Each basis ${\mathbf{d}}_i$ is a time domain signal with the same dimensions as $\mathbf{s}$.

We define a reference point as a pristine scanning location without underground objects in the inspection area. The reference signal $E_{\mathrm{ref}}$ measured from the reference point is expressed, as follows:

$$E_{\mathrm{ref}}\left(t\right)=E{e}^{-\mathrm{j}\omega t},$$

(10)

where $\omega$ refers to the electromagnetic wave frequency. In the dictionary used in this study, each basis represents a reference signal with a specific time delay. For example, ${\mathbf{d}}_i$ denotes the reference signal with a time delay of $\left(L/2-i\right)\mathsf{\Delta }t$:

$${\mathbf{d}}_i=E\left(t+\left(L/2-i\right)\mathsf{\Delta }t\right)=E\left(t+\left(L^{\prime }-i\right)\mathsf{\Delta }t\right)=E{e}^{-\mathrm{j}\omega \left(t+\left(L^{\prime }-i\right)\mathsf{\Delta }t\right)}.$$

(11)

in other words, the measured GPR data $\mathbf{s}$ is represented as a sum of reference signals with different time delay. Any change in $\mathbf{s}$ from the reference signal $E_{\mathrm{ref}}$, such as underground object reflections, is compensated by adding weighted ${\mathbf{d}}_i$ with proper time delay. This is as a time-domain function is represented by a sum of multiple sinusoidal signals with different amplitude and frequency in Fourier transform.

A large L value is preferable for achieving a high resolution in the transformed domain, or a shorter time delay difference $\mathsf{\Delta }t$ in this case. Then, Equation (9) typically represents an underdetermined system of equations (L > T), creating nonunique solutions for α. Additional information is required to select single α from the nonunique multiple possibilities. A unique solution for α can be obtained via the assumption that the actual solution has the sparsest representation, minimizing the number of nonzero entities in α [27]. On the basis of this assumption, Chen and Donoho proposed the basis pursuit approach to solve the underdetermined system in Equation (9) [28]. The sparsest representation is obtained by solving the following with a given dictionary D:

$$\min {||\mathsf{\alpha }||}_1 \mathrm{s}.\mathrm{t}. \mathbf{s}=\mathbf{D}\alpha ,$$

(12)

where ${||\mathsf{\alpha }||}_1$ denotes the ${\ell }_1$ norm of α. Candes provides an open-source implementation solving Equation (12) coded with MATLAB [29], and the proposed technique has been realized based on the code. Then, time domain GPR signals can be sparsely represented while using the dictionary.

Figure 3 presents how to filter background signals using the transformation in Equation (9). $\mathbf{s}$ is the measured signal, $\mathbf{D}$ is defined in Equation (10), and α is solved using Equation (12) so that there is no unknown left. If the current scanning location does not include an underground object (as is the case in the center of Figure 3, the measured GPR signal, $\mathbf{s}$, can be represented by only a single basis corresponding to the reference signal (left of Figure 3). This is the only nonzero element in $\mathsf{\alpha }$, and a sparse representation of the time domain signal is possible in the dictionary domain. $\mathbf{s}$ can then be represented as

$$s=E\left(t\right)=d_{L^{\prime }}=d_{L^{\prime }}{\alpha }_{L^{\prime }}.$$

(13)

The GPR signal measured from a scanning location where there is an underground object in the dictionary domain (right of Figure 3) is represented by Equation (14). The underground object reflection now exists in the GPR signal. The measured signal $\mathbf{s}$ can be represented as the superposition of the reference signal and the reflected signal from the underground object:

$$\mathbf{s}=E\left(t\right)+E\left(t+m\mathsf{\Delta }t\right)={\mathbf{d}}_{L^{\prime }}{\alpha }_{L^{\prime }}+{\mathbf{d}}_{L^{\prime }-m}{\alpha }_{L^{\prime }-m},$$

(14)

where $m\mathsf{\Delta }t$ indicates the depth of the underground object.

As shown in Equations (13) and (14), a new basis ${\mathbf{d}}_{L^{\prime }-m}$ appears only if an underground object is present. Thus, the background signals can be filtered out by removing the basis corresponding to the reference signal ${\mathbf{d}}_{L^{\prime }}$. The resultant dictionary domain data are then inversely transformed into the time domain. Finally, an image of the subsurface region is constructed from the filtered time domain data. Only the reflection signals from underground objects are highlighted in this image, and the other background signals are eliminated.

2.2.3. Instantaneous Phase Analysis for Underground Object Classification

Although the filtered images contain only reflections from underground objects, the nature of these objects is still unknown. Phase information provides a useful means of classifying them, because the phase depends on the permittivity of the corresponding underground object.

Because the permeability of most of the underground objects is similar, Equation (5), which expresses the reflected electromagnetic waves from an underground object, can be rewritten as

$$E_{\mathrm{r}}=\left(\frac{n_1-n_2}{n_1+n_2}\right)E_{\mathrm{i}}=\frac{\sqrt{{\epsilon }_1}-\sqrt{{\epsilon }_2}}{\sqrt{{\epsilon }_1}+\sqrt{{\epsilon }_2}}{E}_{\mathrm{i}}.$$

(15)

If the permittivity of the underground object ${\epsilon }_2$ is lower than that of the surroundings ${\epsilon }_1$, $\frac{\sqrt{{\epsilon }_1}-\sqrt{{\epsilon }_2}}{\sqrt{{\epsilon }_1}+\sqrt{{\epsilon }_2}}>0$ and the reflected electromagnetic waves are in phase with the incident waves. However, a higher ${\epsilon }_2$ will lead to out-of-phase reflected electromagnetic waves. As cavities have low permittivity relative to other underground objects, they can be detected on the basis of the relationship of the phase of the underground reflections to that of the direct waves (Figure 4).

To this end, it is necessary to analyze the phase information at the location of the underground object, but this location is usually unknown. Here, the location of the underground object is estimated as the location of the maximum amplitude of the A-scan signal that passes through the corresponding parabola. A specific number of points are then selected above and below the estimated location of the object. This number of points is determined on the basis of the data collection resolution of the GPR device in the depth direction. The corresponding GPR signal $I\left(x,z\right)$ is then converted to the complex domain while using the Hilbert transform,

$$H\left(x,z\right)=\frac{1}{\pi }P{{\displaystyle \int }}_{-\infty }^{\infty }\frac{I\left(x,z\right)}{z-\tau }d\tau ,$$

(16)

where P denotes the Cauchy principal value and x and z are the spatial coordinates in the scanning direction and the depth direction, respectively. The instantaneous phase value at each spatial point can then be calculated, as follows:

$$\theta \left(x,z\right)={\tan }^{-1}\left(\frac{\mathrm{Im}\left[H\left(x,z\right)\right]}{\mathrm{Re}\left[H\left(x,z\right)\right]}\right),$$

(17)

where Re and Im represent the real and imaginary components of a complex value, respectively.

2.2.4. Decision Making

In the final step, the phase information of the underground reflection is compared with the phase of the direct GPR wave on the basis of the ratio of instantaneous phase change, $\mathsf{\Delta }\theta /\mathsf{\Delta }z$. This ratio is calculated along the region where Equation (16) is applied. The underground object is considered to be in-phase with the direct GPR wave if both of the phase change ratios have the same sign (Figure 4a). Conversely, if their signs are different, they are considered to be out-of-phase, and it may be a high-permittivity object, e.g., a pipe (Figure 4b).

3. Simulation Result: Numerical Validation Using GPR Simulation Data

3.1. Simulation Setup

The proposed underground object classification technique was numerically validated using simulated 2D GPR data for various underground models. The proposed technique was applied to a pristine model, a model with an underground cavity, and a model with an underground pipe. The visibility of the underground object was enhanced by removing background signals, and the cavity and pipe were distinguished from each other on the basis of phase information.

Open-source software, gprMax [30], was used to simulate electromagnetic wave propagation while using the finite-difference time domain method. Figure 5 depicts the composition of the simulated underground model. The model has a depth of 0.8 m and a width of 2.0 m, with a 0.2-m-thick upper air layer. The spatial resolution of the model is 2 mm. To simulate the process of GPR data collection, a virtual GPR device scans the model surface from left to right at intervals of 20 mm, resulting in 100 GPR signals. The virtual GPR device consists of a single channel transmitter and a single channel receiver, 40 mm apart. The electromagnetic wave is a normalized first derivative of a Gaussian curve with a center frequency of 1.6 GHz.

The white circle at 0.35 m depth in Figure 5 represents an underground object. For a pristine model, this region is filled with soil (${\epsilon }_{\mathrm{s}}=5)$. For models with an underground cavity and pipe, the region is filled with air (${\epsilon }_{\mathrm{a}}=1)$ and a perfect conductor (${\epsilon }_{\mathrm{p}}=\infty )$, respectively, to simulate the corresponding underground object.

3.2. Classification Results for the Pristine Model

The proposed technique was first applied to the pristine model. The raw GPR image in Figure 6a represents a cross-sectional view of the model. As expected, no reflection from an underground object is observed, but dominant surface reflections are seen near the ground surface. The background filter was applied to this region to identify any reflections from underground objects while using the leftmost scanning location of the model as the reference point and the corresponding A-scan response as the reference signal. No distinguishable features are observed in Figure 6b after filtering, indicating that the proposed technique is valid for the pristine model and does not give a false positive.

3.3. Classification Results for Models with Underground Objects

Next, the feasibility of the proposed technique for classifying underground objects was examined. First, visualization enhancement via the proposed technique was studied while using models with a cavity (Figure 7a) and a pipe (Figure 7b). Figure 7c,d shows that the background filtering effectively removes the dominant surface reflections, allowing the underground objects to be clearly visualized. In addition, the location of the first reflection wave packet perfectly matches the actual location of the underground object, which is represented by a red circle.

However, it is not possible to classify the type of underground object on the basis of the filtered images, as both objects are represented by parabolas. The instantaneous phase information for each case was therefore examined in the filtered images and compared with the direct wave phase, corresponding to the blue circles in Figure 7c,d. For the cavity (Figure 7e), the phase of the direct wave (blue box) and that at the location of the underground object (red box) match each other. This in-phase reflection implies that the cavity has lower permittivity relative to the surrounding soil. This can be distinguished more easily after the binarization of the raw phase signals (Figure 7g). Binarization transforms the phase value to −1, 0, or 1, simplifying the comparison with the direct wave phase. The phase change ratio is positive for both direct waves and the cavity reflection, indicating an in-phase relationship. However, at the underground pipe (Figure 7f), the reflection is out-of-phase with the direct wave because of its higher permittivity. The phase change ratio of the pipe reflection is negative, as shown in Figure 7h. This numerical simulation demonstrates that cavities and pipes can be clearly distinguished while using instantaneous phase information, as proposed in the previous section.

3.4. Robustness of the Proposed Technique to Complicated Underground Features

The proposed technique was applied to various simulation models with complicated underground features. First, the sensitivity of the proposed technique was studied with different diameters (2 cm and 25 cm) of the underground cavity. As shown in Figure 8, the predominant surface reflections were effectively removed after background filtering and the smallest cavity with a 2 cm diameter was clearly visualized. In addition, the location of the first reflection wave packet matched very well with the actual cavity location (red circles).

Second, the effect of the soil permittivity to the proposed technique was studied. Here, the soil permittivity ${\epsilon }_{\mathrm{s}}$ in Figure 5 was changed from 5 to 2 and 10. For the case ${\epsilon }_{\mathrm{s}}=2$, smaller permittivity difference between the soil and the cavity weakened cavity reflections (Figure 9a). But, the proposed technique still highlights them (Figure 9b). Note that the GPR signals are dispersed to the depth axis due to the slower electromagnetic wave speed in low permittivity media. On the other hand, for the case ${\epsilon }_{\mathrm{s}}=10$, stronger cavity reflections are observed, as predicted in Equation (5). There is no dramatic change in pipe reflections as the pipe permittivity was already set very high (${\epsilon }_{\mathrm{p}}=\infty$) and the permittivity difference did not largely change after soil permittivity was changed. The dispersed GPR signals are still observed in Figure 9e,g.

Next, multiple underground objects existed in a single model. In Figure 10a, a cavity and a pipe were placed in parallel with 10 cm distance between them. The interactions of each electromagnetic wave reflections are clearly visualized in raw GPR image (Figure 10b), and still exists after removing the background signals (Figure 10c). However, the binarization of the phase signals clearly indicates the in-phase cavity reflection (Figure 10d) and out-of-phase pipe reflection (Figure 10e). This result implies that the proposed technique can identify multiple underground objects in the inspection region.

Finally, the effect of complex underground compositions was also studied. An additional sand layer (${\epsilon }_d=3$) is located below (Figure 11a) and above (Figure 11b) the cavity, respectively. Though this layer creates another dominant reflection signals (Figure 11c,d), the proposed technique effectively removes them and demonstrates the robustness of the proposed technique to the underground disturbances (Figure 11e,f). The phase of the direct wave (blue box) and that at the location of the underground object (red box) match each other in both cases, implying their in-phase reflections. Note that additional phase packets.

4. Experimental Result: Field Validation Tests Using 3D GPR Data

4.1. Experimental Setup

The proposed underground classification technique was experimentally validated by generating visualizations of underground cavities, buried pipes, and underground gravels from 3D GPR data collected from urban roads in Seoul, South Korea. The performance of the proposed method in analyzing a pristine region was also investigated to validate its robustness.

The inspected region was a 0.7 km stretch of an urban road in a residential area near Yaksu subway station, located in Jung-gu, Seoul, as shown in Figure 12a A DXG1820 GPR antenna (3D-RADAR), which includes 20 T and R channels corresponding to a 1.5 m scanning width, was used for this experiment. The measured GPR signals were collected with a GEOSCOPE MK IV data acquisition system that was provided by 3D-RADAR, which has a time resolution of 0.35 ns and a maximum sampling rate of 13,000 Hz. The input signal has a frequency range of 200–3000 MHz, with a step-frequency input waveform. Figure 12b shows the equipment that was used for field data collection. Automatic data collection was carried out by integrating the GPR device into a vehicle. The average scanning speed was approximately 5 km/h.

The region inspected was composed of four layers: a surface layer, a basecourse, a sub-base, and a road bed. The uppermost surface layer is a 10-cm-thick asphalt pavement. This is underlain by a 20-cm-thick basecourse (usually consisting of construction aggregate). The third layer is a 30-cm-thick sub-base of unbound granular or cement-bound materials. The road bed layer is typically much thicker than the other layers, depending on the ground conditions. Four representative regions of the 0.7 km inspection area were selected for the field validation of the proposed technique, as shown in Figure 13. The pristine region does not contain any underground objects. In Region 1, there is an underground cavity (gray circle with a solid outline) in the subbase layer. A pipe (brown circle with a dashed outline) is buried in Region 2 and it is perpendicular to the scanning direction. This region is inspected as an example of common artificial underground objects in urban road. There is a punctual gravel (black circle with a dotted outline) in the basecourse layer of Region 3. Gravels widely exist in urban roads but usually misclassified as cavities with their similar size and shape. The classification robustness of the proposed technique is examined with this region.

4.2. Underground Object Classification in the Pristine Region

The proposed technique was first applied to the pristine region. This region comprises the four underground layers, but there are no underground objects. The raw GPR image that is shown in Figure 14a is a representative B-scan side view at a width of 1.275 m. There are no visible reflections from underground objects, but there are clear dominant surface pavement reflections near the ground surface. Background filtering was applied to the image to locate any reflections from underground objects while using the leftmost scanning location of the region as the reference point and the corresponding response as the reference signal. This results in Figure 14b, in which no distinguishable features are observed. Moreover, no objects are identified in the three C-scan images on the right side of Figure 14b, which correspond to the red lines in the B-scan image. These results indicate that the proposed technique is suitable for analyzing pristine regions without giving false positives.

The instantaneous phase information of the direct waves corresponding to the blue circle in Figure 14b was analyzed and is represented by a blue box in Figure 14c. The phase change ratio is positive for the direct waves (Figure 14d).

4.3. Underground Object Classification in the Regions with Underground Objects

The proposed technique was applied to the three regions with different underground objects. For all three regions, the proposed technique effectively visualized and classified the underground objects.

In Region 1, a cavity is located in the sub-base layers. However, there is no clear reflected signal from a cavity in the side view of the raw GPR image at 0.375 m width (Figure 15a), owing to the dominance of surface pavement reflections. After background filtering, the cavity reflections become much more apparent (Figure 15b). These reflections are also consistent with the actual cavity location, which is indicated by a red circle. A circular object can be seen in three C-scan images around the location of the observed object (red lines in Figure 15b). The instantaneous phase is shown in Figure 15c, with the phase corresponding to the cavity reflection being indicated by a red box. It has a positive phase change ratio (Figure 15d).

The next example is a pipe; these are very common under urban roads. The pipe in Region 2 is not visible in the raw GPR image (Figure 16a), but its location is clear after background filtering (Figure 16b). The C-scan images show that the pipe runs perpendicular to the scanning direction. The C-scan images show a pipe shape that is distinct from the cavity case (Figure 15), thus making it possible to distinguish a pipe from a cavity in 3D GPR images. In addition, the reflection from the pipe (Figure 16c,d) is out of phase with the direct wave in Figure 14 and it has a negative phase change ratio.

The last example is a gravel. Gravels are widely observed in urban underground GPR images, but it is challenging to distinguish them from a cavity on the basis of 3D GPR images only. While other underground structures, including pipes, can be easily differentiated from their shape in the C-scan images, a gravel is often misclassified as a cavity with their similar size and shape. Although an object is visualized in both the B-scan and C-scan (Figure 17b) images after background filtering, there are no grounds for classifying it as a cavity or a different type of underground object. However, the reflection from the identified object (Figure 17c,d) is out-of-phase with the direct wave in Figure 14, while the cavity reflections show in-phase reflections in Figure 15.

5. Discussion

5.1. Underground Object Classification Compared with the Conventional Subtraction Method

The underground object classification performance of the proposed technique is compared with the conventional subtraction method [12,13]. The subtraction method, which is being widely used in GPR data and image processing, eliminates the dominant surface reflections by subtracting a reference signal from the raw GPR data. However, this method is highly susceptible to noise and small signal disturbances in the time domain, e.g., variations in pavement thickness or unexpected arbitrary boundary conditions in the pavement.

Figure 18 clearly represents the limitation of the subtraction method. For each region, the dominant surface reflections are not fully eliminated, even after applying the subtraction method (Figure 18a,c,e). As the subtraction method assumes identical underground background reflections from every spatial point, even a small measurement disturbance from field condition affects its elimination performance. Another limitation is that the underground layer boundaries are also highlighted after applying gain enhancement. While the proposed technique extracts the object reflections information while using the basis pursuit approach, the subtraction method cannot distinguish the object reflection and the underground layer reflection and highlights both of them if they are placed nearby. The red boxes in Figure 18a,c,e indicate the uneliminated soil layer boundaries even after applying the subtraction method, while they are successfully eliminated using the proposed technique in Figure 18b,d,f. The basis pursuit approach can compensate the disturbance between the reference signal and the actual measurement, leading the effective elimination of the background reflections. As already presented in Section 3.4, the proposed technique has an effective underground layer reflection elimination performance and highlights only the underground object.

5.2. Underground Object Classification Using a Phase Change Ratio

Another limitation in the conventional GPR analysis, including the subtraction method, is that they do not tell about the nature of the identified underground objects. The proposed phase analysis technique is thus applied to classify the object. The phase information and the phase change ratio are analyzed for Regions 1 to 3 in Figure 15, Figure 16 and Figure 17, respectively, and the results are summarized in

Table 1. Comparisons of the phase change ratio in each region with the reference region.

Region	Underground Object	Phase Change Ratio	Comparison with the Reference
Reference	-	+25.4 m−1	-
1	Cavity	+50.9 m−1	In-phase
2	Pipe	−43.7 m−1	Out-of-phase
3	Gravels	−25.4 m−1	Out-of-phase

.

For the cavity reflection in Region 1, it has a positive phase change ratio (Figure 15d) and it matches the direct wave phase in the reference region (Figure 14d). It is thus an in-phase reflection, implying that the object has low permittivity in comparison with the surrounding soil and can therefore be classified as a cavity.

On the other hand, the reflection from the pipe (Figure 16c,d) is out of phase with the direct wave in Figure 14 and has a negative phase change ratio. Though the actual pipe material information is not provided to the authors, the pipe material is estimated to have higher permittivity than soil, possibly metal, from the observed phase inversion. Please note that the reflection from the pipe can be in phase if it has lower permittivity than the soil (e.g., polyethylene or Teflon) and the only use of phase information can lead to false alarms. However, the proposed technique still identifies the pipe while using both its background filtered C-scan images and the phase information, making it possible to be complementary to each other in underground object identification process. The negative phase change ration in Region 3 (Figure 17d) also indicates that the object is not a cavity and must be an underground object with a high permittivity, possibly a gravel.

Therefore, as the phase of the underground object reflection is related with the object permittivity, it is possible to have an additional information for underground object classification via the proposed phase analysis.

5.3. Limitations of the Proposed Technique

However, certain technical limitations of the proposed technique should be addressed before it can be applied to more realistic conditions. Careful selection of the reference signal is required, as the performance of the background filtering technique strongly depends on reference signal selection and dictionary construction. For example, if the reference signal includes any undesired underground object, the background filter will not work properly. Follow-up studies are underway to automate the reference signal selection process by using the machine-learning techniques.

Development of a fully automated identification system through integration with a GPR scanning system and further validation with multiple in situ datasets are also required.

6. Conclusions

This study proposes a robust technique for classifying underground objects while using (1) background filtering of 3D GPR data to remove environmental noise and field deviations and (2) phase analysis to distinguish underground cavities and possible sinkholes from other underground objects. The proposed technique enhances the visibility of underground objects by highlighting the electromagnetic waves reflected only from underground objects. The permittivity of the identified objects is then compared with that of the surrounding soils using the corresponding instantaneous phase information. The object is classified as a cavity if its permittivity is lower than that of its surroundings. The performance of the proposed technique was validated while using numerical simulations and field data obtained from urban roads in Seoul, South Korea. The use of such field-test data is particularly useful because of its rarity in previously published literature.

The authors are now extending the proposed technique to utilize a machine-learning-based automatic cavity classification system. The proposed technique provides additional information with which to train such a system, such as enhanced images and phase change ratios, promoting faster and more robust classification results than is possible through the time-consuming efforts of human experts.