1.

Introduction

A popular method for nondestructive detection of underground targets is the remote sensing techniques. One of the most successful modalities for remote sensing of underground targets is the ground-penetrating radar (GPR), and it has been extensively used in many fields, such as military,^1,2 civil engineering,^3,4 and geological survey.⁵ GPR emits electromagnetic wave toward the ground, due to the different relative permittivity of the underground medium, the electromagnetic wave will be reflected back. When GPR antenna moves across the buried target, the time of arrival for responses from the target constitutes a hyperbola in the B-SCAN image. Therefore, the buried target detection can be viewed as mapping from the object space to an image space. Then, according to the characteristic curve of the target, the geometric parameters such as buried depth, width and height can be calculated.

In recent years, the main methodologies used to detect hyperbola in GPR B-SCAN images include curve fitting,⁶ machine learning,⁷ and deep learning.⁸ Hough transform technique⁹ is used successfully for deformed objects fitting, but the discretization of a large number of parameters would lead to greater computational complexity. Hough transform and least-squares (LS)-based hyperbolic asymptotes fitting algorithms were proposed to improve the computational efficiency.^10,11 However, the interaction between multiple hyperbolas in B-SCAN profiles will greatly reduce the accuracy of LS curve fitting results, and it is limited in practical application. Machine learning methodologies,^12,13 which use machine vision technology to extract image attributes and classify them, have been used in underground target detection. Two methods were most commonly applied in buried threat detection, i.e., Viola–Jones algorithm¹⁴ and image feature descriptor algorithm.^15,16 For these works, feature descriptors need to be designed by manual operation with experience, and the detection accuracy is not satisfactory.^17,18 Recently, deep convolutional neural network (DCNN) algorithm has achieved impressive performance for image segmentation and recognition tasks on optical and medical images,^19–21 and the use of the DCNN for the identification of hyperbolic curves in B-SCAN profiles has been increasingly frequent.^22–24 Due to the complexity of the geological environment, the structural features of B-SCAN images are usually complex, which also makes the application of CNN in GPR face some challenges. For example, there is a lack of data sets for training network models and automatically extracting the characteristic curve of the target accurately. To address the scarcity of GPR training data, Lei et al.²⁵ enhanced the training data set through rotation, gray-scale transformation, and scale transformation. Zhang et al.²⁶ and Wang et al.²⁷ used generative adversarial networks (GAN) model-based strategies to generate more realistic GPR images. In addition, to improve the accuracy of hyperbolic detection, the researchers studied the recombination of classic feature extraction and classification networks. Zhang et al.²⁸ proposed DCNN and incremental random sampling method, in which Resnet50 network was used to extract features, and YOLO V2 network was used for feature recognition, to detect moisture damages in asphalt pavements. Ozkaya et al.²⁹ adopt residual CNN + Bi-LSTM model to improve classification accuracy of GPR type, scanning frequency, and soil type. However, these methods are data-driven and do not use the physical mechanism of GPR work, and the accuracy of the results largely depends on the size of the sample data set used for training the network model. In this work, we have proposed a method based on the combination of data- and model-driven. On the one hand, the physical mechanism of GPR is used as a constraint to design a feature extraction network to obtain the salient feature hyperbola of the target in the B-SCAN image; on the other hand, the CycleGAN³⁰ is used to enhance data for the training of the CNN model. Further, on the basis of extracting the salient characteristic hyperbola of the target, it is proposed to detect the weak useful information (the lower boundary characteristic curve under the upper boundary characteristic curve of the target) in the local space.

In general, this article adopts a method based on the combination of data and GPR physical model to detect and identify the characteristic information of B-SCAN images, and discusses in detail the following three aspects: (1) The characteristics of the target scattering electric field in the GPR echo data; (2) the scarcity of GPR training data; (3) effective extraction of target features, especially the extraction of weak feature information. The main contributions could be summarized as follows.

The remaining part of this paper is organized as follows. Section 2 presents the GPR target detection model. In Sec. 3, we propose CNN-based schemes to detect the features of B-SCAN images, and GPR data augmentation method. In Sec. 4, the experimental results are shown and discussed in Sec. 5. Finally, conclusions are drawn in Sec. 6.

2.

GPR Target Detection Model

GPR detects subsurface structure by emitting electromagnetic wave pulse signals into the ground, there are reflection, refraction, and other electromagnetic physical phenomena at the interfaces of different mediums. The reflected wave is received by the GPR receiving antenna, and the refracted wave enters the next layer of media. According to the Snell’s law, the reflection coefficient $\gamma$ and transmission coefficient $T$ at a dielectric boundary are given by

Fig. 1

Sample GPR data from GPRMAX simulation software.

In Fig. 1, $C_U$ and $C_L$ are the upper and lower boundary characteristic curve of the target, respectively. $D$ is the buried depth, $W_U$ and $W_L$ are the width of the target, $H$ is the higher of the target. $D$ can be calculated according to the vertex position of $C_U$, $H$ can be calculated according to the vertex position of $C_L$ and $C_U$, $W_U$ and $W_L$ can also be calculated according to the effective widths of curves $C_U$ and $C_L$, respectively, which are shown as follows:

The intensity of $C_L$ is weak compared with the characteristic curve of $C_U$. To improve the accuracy of the detection results of $C_L$, the lower boundary characteristic curve of the target is extracted based on the amplitude parameter information of the scattering field combined with the phase parameters in this paper. When the electromagnetic wave is reflected, the phase information changes in the location of the medium boundary in the scattering field. The spherical electromagnetic wave propagating is considered in subsurface medium, which can be expressed as

According to the above analysis of the physical mechanism of GPR, the proposed CNN feature extraction network based on GPR physical model constraints is well suited for the GPR subsurface targets detection for the following reasons:

3.

Deep Learning Network Architecture for GPR B-SCAN Image Processing

The acquired GPR echo data usually contain various clutters. As the transmitting and receiving antennas of GPR work synchronously, the amplitude of the generated direct coupled waves far outweigh that of the target. Moreover, there are attenuation and diffraction phenomena combined with the complex geological environment factors when electromagnetic wave penetrates a lossy medium. These issues pose a challenge to detect the underground targets. To address these issues, we detect underground targets on the basis of building a deep learning network model through the data-driven and GPR physic-based methods.

The flowchart of the data processing is shown in Fig. 2, which could be divided into three steps. The DCNN is first designed to extract the high-level semantic feature information through convolution and pooling operations to process the direct wave interference with high intensity in GPR echo data, on this basis, the regional location of the direct wave in B-SCAN image is determined and removed. Then, the B-SCAN image is input into the multilayer feature fusion DCNN to extract the feature map that contains more complete target information, hereafter outputting the feature map to the output layer. Finally, training the classification network based on the obtained feature maps; thus, the geometric parameters of the target can be further calculated using the extracted feature curves of the subsurface targets.

Fig. 2

GPR data processing flowchart based on deep learning.

Additionally, in the process of training the classifier network, the CycleGAN is used for generating realistic B-SCAN images to achieve data augmentation in processing GPR field data due to the lack of real data for training. In this paper, we use the classification network based on Faster R-CNN,³² the network structure and its application in GPR B-SCAN hyperbolic image detection can be referred to Ref. 27. Next, we will introduce the B-SCAN image feature extraction network and data augmentation network, respectively.

3.1.

GPR B-SCAN Image Feature Extraction Network

The process of extracting B-SCAN image features through convolution operation based on CNN can be calculated as

Eq. (7)

$$\{\begin{array}{l}{\mathcal{F}}_{\text{target}}=f_{\Sigma }[{\mathcal{F}}_s^{(L)}(s=1,2,3,\dots ,S)|\mathcal{I}\stackrel{L}{\to }{\mathcal{F}}_S^{(L)}]\\ {\mathcal{F}}^{(l)}={\mathcal{I}}_{w\times h\times M}*{\mathcal{C}}_{k,N,D}\\ {\mathcal{F}}^{(l+1)}=\text{pooling}({\mathcal{F}}^{(l)})\end{array},$$

where ${\mathcal{F}}_{\text{target}}$ is the feature map extracted by CNN, $f_{\Sigma }[·]$ means feature fusion, $\stackrel{L}{\to }$ represents the convolution and pooling feature mapping through $L$ times, ${\mathcal{F}}_S^{(L)}$ represents $S$ feature maps obtained after $L$ times of mapping. ${\mathcal{F}}^{(l)}$ is the feature map output by the $l$’th layer of CNN, “*” is the convolutional operation, $\mathcal{I}$ is the input image, and $w$, $h$, and $M$ are the width, height, and number of channels of the input image, respectively. As the B-SCAN image corresponding to GPR echo data is single channel, the channel number $M=1$. $\mathcal{C}$ is the convolution kernel in the convolution layer, $k$, $N$, and $D$ are the size, number, and channel of the convolution kernel, respectively, where the number of channels in the convolution kernel equal to the number of channels in the input image, i.e., $D=M=1$, $\text{pooling}(·)$ is the max pooling function. The extracted feature image needs to further identify the target and noise in ${\mathcal{F}}_{\mathrm{B}\text{-}\mathrm{SCAN}}$, in reflectance measurement mode, the target in GPR echo data usually has specific structural characteristics, and the noise is generally random. Therefore, the sample set is constructed based on ${\mathcal{F}}_{\mathrm{B}\text{-}\mathrm{SCAN}}$ and manual labeling to train the classification network, which is used to identify the noise and target with supervised learning method.

Designed a cascade structure of CNN, and according to the characteristics of direct wave interference signal strength, the first-level DCNN is used for the extraction of deep semantic features of the B-SCAN image to determine the location of the direct waves, and the corresponding network structure is shown in the upper half region (a) in Fig. 3. Then for the GPR B-SCAN image with the direct wave noise removed, the DCNN with multiscale feature fusion is adopted to extract more complete target feature information, the corresponding network structure is shown in the lower half region (b) in Fig. 3.

Fig. 3

The brief architecture of proposed CNN for GPR B-SCAN image feature extraction.

The feature map extraction network is based on CNN, where the convolution kernel size of the convolution layer is $3\times 3$, the activation function is a rectified linear unit function,³³ and the filter size of the max pooling layer is $2\times 2$. The network structures shown in Figs. 3(a) and 3(b) are designed by the characteristics of the network input images, respectively. In Fig. 3(a), the input signal is the raw GPR echo data, in which the intensity of the direct wave is far outweigh the target, and high-level semantic features are extracted through multiple convolution and pooling operations, i.e., through 10-layer convolution and pooling processing, 512 feature maps with a size of $53\times 15$ are obtained. In the output layer, considering that the characteristics of the direct wave signal are usually horizontal and linear, and the texture structure is relatively simple, the pooling result S5 is up-sampled through cubic interpolation to make the output feature map with the same spatial dimensions as the input B-SCAN image. As shown in Fig. 3(b), through 14-layer convolution and pooling processing, 512 feature maps with a size of $126\times 30$ are obtained. The input of the network is a B-SCAN image containing clutters with the direct wave removed (zero setting), where the clutters are formed by scattering and diffraction when electromagnetic waves propagate in complex geological scenes. A multilayer feature fusion strategy is adopted in network design to extract feature maps containing more target information. Accordingly, the M1 layer is the result of feature fusion of the C8 layer and the S7 layer. The size of the feature maps of the C8 layer and the S7 layer are $1008\times 240$ and $504\times 120$, respectively. For the feature map of the C8 layer, a pooling core with a size of $2\times 2$ is used for average pooling to make it consistent with the size of the feature map of the S7 layer. The M2 layer is the result of feature fusion of the C10 layer and the S8 layer, it is similar to the process of M1 layer.

The underground target feature detection models were trained based on 4600 labeled images. Figure 4 shows the process of proposed model fitting.

Fig. 4

The train loss plot of the proposed model fitting.

The main parameter settings in the network model training process are as follows: The initial learning rate is set to 0.0001, the max epoch is equal to 25, a total of 115,000 iterations, and output the current epoch loss after every 100 iterations. The loss function takes the cross-entropy loss, and the optimization process uses the stochastic gradient descent momentum method. As shown in Fig. 4, when epoch is approximately equal to 13, the model tends to converge after 50,000 iterations.

Remark 1:

The CNN of the cascade structure is designed with GPR physical model as a priori condition. Direct wave and target signal are processed based on the amplitude and structural features in B-SCAN image, respectively:

1. Direct wave signal has horizontal linear structure characteristic and large amplitude. The focus of processing is to locate the area where it is located in the B-SCAN image to remove this interference signal.

2. In the B-SCAN image after removing the direct wave, the target exhibits the characteristics of a hyperbolic structure with a downward opening and a large amplitude. The semantic feature detection network is mainly used to extract the hyperbola with complete structural features.

3.2.

GPR B-SCAN Data Augmentation Network

To transform two groups of GPR B-SCAN images with different styles to generate new images to achieve data augmentation, we assume that the real GPR background image data set, which is denoted as $\mathcal{A}$, which mainly contains reference scenario information of various complex texture structure features. GPRMAX toolbox is used to simulate objects (underground cavity and aquifer) in different spatial positions in uniform media, and the obtained B-SCAN image data set, which is denoted as $\mathcal{B}$, which mainly contains the hyperbolic geometric structure characteristic of the target. The data set composed of B-SCAN images with complex background texture structure features of $m$ pieces is denoted as $\mathcal{A}=\{T^1,T^2,\dots ,T^m\}$, where $T^i(i=\mathrm{1,2},\dots ,m)$ is the $i$’th image complex background. The data set constituted by $n$ B-SCAN images with hyperbolic geometric structure features is denoted a $\mathcal{B}=\{H^1,H^2,\dots ,H^n\}$, where $H^j(j=\mathrm{1,2},\dots ,n)$ is the $j$’th B-SCAN image containing the target. The B-SCAN image containing only complex background texture structure features can be represented as $T_A=\{T|T_C\}$, where $T_C$ is the complex texture features in image $T$. The B-SCAN image with hyperbolic geometric features and a simple background texture can be represented as $H_B=\{H|H_C;T_S\}$, where $H_C$ and $T_S$ are the hyperbolic geometric structure feature and simple background texture structure in image $H$, respectively. Furthermore, the B-SCAN images in data set $\mathcal{A}$ and $\mathcal{B}$ are mapped and transformed, where the images in data set $\mathcal{A}$ have complex background texture features, the images in data set $\mathcal{B}$ have both hyperbolic geometric structure features and simple background texture features. Thus the generated B-SCAN image has both complex background texture features and hyperbolic geometric structure features. This mapping transformation can be described by

Eq. (8)

$$f(T_A\to H_B)=\{\widehat{H}|H_C;T_S\leftrightarrow T_C\},$$

where → is mapping transformation, $\widehat{H}$ is the generated image, $\leftrightarrow$ is transfer. Equation (8) indicates that the integrity of the hyperbolic structure should be maintained in the generated B-SCAN image, and simple background texture and complex background texture are transferred.

Furthermore, CycleGAN can be used for image style transfer with different texture features in data set $\mathcal{A}$ and $\mathcal{B}$ (unpaired), i.e., the mapping transformation can be established by learning the features of images between data set $\mathcal{A}$ and $\mathcal{B}$, and the hyperbolic structure features are retained in the generated B-SCAN image, but the background texture features are transferred. The structure of CycleGAN is shown in Fig. 5.

Fig. 5

Architecture of CycleGAN.

The upper and lower parts of a CycleGAN are dual, and form a circular structure. For the mapping function $G_{A\to B}$ and its discriminator $D_B$, we can introduce an adversarial loss function $L_{\mathrm{GAN}}(G_{A\to B},D_B,A,B)$. Similarly, the mapping function $G_{B\to A}$ and its discriminator $D_A$ corresponding adversarial loss function is $L_{\mathrm{GAN}}(G_{B\to A},D_A,B,A)$, which is calculated as

Eq. (9)

$$L_{\mathrm{GAN}}(G_{A\to B},D_B,A,B)=E_{H\sim p_{\text{data}}(H)}[\log D_B(H)]+E_{T\sim p_{\text{data}}(T)}[\log (1D_B(G_{A\to B}(T)))],$$

Eq. (10)

$$L_{\mathrm{GAN}}(G_{B\to A},D_A,B,A)=E_{T\sim p_{\mathrm{data}}(T)}[\log D_A(T)]+E_{H\sim p_{\mathrm{data}}(H)}[\log (1-D_A(G_{B\to A}(H)))],$$

where $E$ is the mean function. To ensure that the learned mapping functions are cycle-consistent, we can use a cycle consistency loss

Eq. (11)

$$L_{\text{cycle}}(G_{A\to B},G_{B\to A})=E_{T\sim p_{\text{data}}(T)}[{\Vert G_{B\to A}(G_{A\to B}(T))-T\Vert }_1]+E_{H\sim p_{\text{data}}(H)}[{\Vert G_{A\to B}(G_{B\to A}(H))-H\Vert }_1],$$

where ${\Vert •\Vert }_1$ is L1 norm, the full objective is

Eq. (12)

$$L(G_{A\to B},G_{B\to A},D_A,D_B)=L_{\mathrm{GAN}}(G_{A\to B},D_B,A,B)+L_{\mathrm{GAN}}(G_{B\to A},D_A,B,A)+\lambda L_{\text{cycle}}(G_{A\to B},G_{B\to A}),$$

where $\lambda$ is usually an empirical value, which controls the relative importance of adversarial loss and cycle consistency loss. We aim to solve the optimization objective function $L(G_{A\to B},G_{B\to A},D_A,D_B)$ through maximizing $D_A$ and $D_B$ to improve the discrimination ability of the network, and minimizing $G_{A\to B}$ and $G_{B\to A}$ to generate realistic images, i.e.,

Eq. (13)

$$(G_{A\to B}^{*},G_{B\to A}^{*})=\arg \underset{G_{A\to B},G_{B\to A}}{\min }\underset{D_A,D_B}{\max }L(G_{A\to B},G_{B\to A},D_A,D_B).$$

The Adam optimization algorithm³⁴ can be further used to optimize Eq. (13).

Remark 2:

The data augment method based on the GANs can generate realistic GPR data to address the scarcity of GPR data, which is different from traditional data augment methods, such as image translation, scale, cropping, rotation, etc. Especially in dealing with the challenging problem of lack of field data, GAN-based data sets can be better used for the training of classification networks. On this context, the proposed method can be better applied to target detection in complex scenes.

4.

GPR Target Detection Experiment

4.1.

Data Description

A total of 5050 GPR B-SCAN images are used for experiments, where 4600 images are used for training the classification network, and 450 images are used for testing. The data are generated through GPRMAX software simulation, CycleGAN generation, and GPR field measurement. Specifically, 4000 images are obtained using GPRMAX toolbox, where 1000 images are simulated for cavity and aquifer with 1 and 2 targets, respectively. Due to the lack of field data, the data generation method based on CycleGAN is adopted to generate 1000 realistic GPR B-SCAN images, which uses the following two data sets, i.e., 1000 background B-SCAN images are collected in real geological scenes, and 1000 target foreground B-SCAN images are simulated in scenes with uniform distribution of background medium. 50 images for experimental verification obtained by noise addition, sandpit experiment, and field measurement. GPRMAX simulation are carried out under the framework shown in Fig. 6.

Fig. 6

Geometry of the simplified 2-D half-space.

Where $\mathrm{\Omega }c$ is the concrete slab layer, $\mathrm{\Omega }s$ is soil layer, and $S_1$ is the contact surface of the two layers. The upper half-space is filled with air, and in the lower half-space we consider a two-layered background medium where the upper one is concrete slab layer $\mathrm{\Omega }c$, the lower one is soil layer $\mathrm{\Omega }s$. The contact surface of the two layers $S_1$ is rough. The targets, including cavities and aquifers, are located in layer $\mathrm{\Omega }s$. A coordinate system shown in Fig. 6 is set up to facilitate description. The transmitting and receiving antennas of GPR move along the $x$ direction closely to the concrete pavement. The scenario is set up as follows.

• 1. The width along the $x$ direction is 5 m, the depth along the $y$ direction is $d1+d2=2\mathrm{m}$, the average thickness of concrete is $d1=0.25\mathrm{m}$.

• 2. There are four to six concave and convex blocks with a maximum width and height 0.3 and 0.05 m, respectively, which randomly generated at the lower boundary of concrete to simulate the rough interface ($S_1$) between concrete and soil.

• 3. The dielectric constant of soil along the $y$ direction is generated in a random way, considering the influence of water content in soil at different depths, the range of is set as 9 to 25.

• 4. The number of targets is set as 1 or 2, and the shape is set to rectangle. (The size of the cavity is randomly generated between 0.4 to 0.5 m in width and 0.3 to 0.4 m in height; the width and height of the aquifer are randomly generated between 0.5 to 0.7 m and 0.25 to 0.35 m, respectively.) A single target is randomly distributed in $\mathrm{\Omega }s$ region, and two targets are randomly distributed in the left and right half of $\mathrm{\Omega }s$, respectively.

4.2.

GPR Data Augmentation

It is difficult to acquire a large number of B-SCAN images containing cavities and aquifers by field measurement. Data set $B$ contains 1000 B-SCAN images simulated by GPRMAX, where the target is set in a uniform media scenario. Data set $A$ containing 1000 B-SCAN images of the underground background (without targets and removed the direct-wave) is obtained by field measurement. CycleGAN is used to convert images with different styles in data set $A$ and $B$, to generate realistic GPR B-SCAN images. Four B-SCAN images are selected from data set $A$ and $B$ for analysis.

Figure 7(a) shows GPR field images with complex background texture features. Figure 7(b) shows simulated B-SCAN images with hyperbolic geometric structure features of the target and texture features of the simple background. CycleGAN is used to transform simple background texture and complex background texture, and the generated B-SCAN images have both hyperbolic geometric structure features and complex texture background features, which is as shown in Fig. 7(c). Thus, the generated realistic B-SCAN images, which can be used to train the classification network to improve the accuracy of target detection results in GPR field data.

Fig. 7

Generated realistic GPR B-SCAN based on CycleGAN. (a) Field background B-SCAN images; (b) simulated cavity B-SCAN images; and (c) generated B-SCAN images.

4.3.

Direct Wave Detection

In the raw GPR echo data, the direct wave signal has high amplitude. According to this salient feature, the designed DCNN is used to detect and remove the direct wave in a GPR B-SCAN image.

As shown in Fig. 8(a), the intensity of the direct wave signal far outweighs the target signal; it is difficult to judge whether there is a target in the B-SCAN image from subjective vision. Figure 8(b) shows the extracted direct wave feature map. It can be seen that the method of feature extraction using DCNN can accurately locate the region of the direct wave in the B-SCAN image. Figure 8(c) shows that the B-CAN image after the direct wave is removed in Fig. 8(a), where the texture feature information can be seen. It should be noted that the “B-SCAN images” in the following refers to B-SCAN images all preprocessed with the direct wave removed.

Fig. 8

Example of removing direct wave from B-SCAN image. (a) The raw GPR B-SCAN image; (b) the direct wave feature map; and (c) GPR B-SCAN image with direct wave removed.

4.4.

Target Feature Curve Extraction

The real underground environment is complex and variable, and the collected GPR B-SCAN images contained various clutters. Figures 9(a) and 9(b) are B-SCAN images obtained by simulation, which correspond to one cavity and two aquifer targets, respectively. Taking them as input to the DCNN with structure shown in Fig. 3(b), and the extracted feature maps are shown in Figs. 9(c) and 9(d).

Fig. 9

Feature map is extracted from the B-SCAN image. (a) B-SCAN image with one target; (b) B-SCAN image with two targets; (c) feature map of Fig. (a); and (d) feature map of Fig. (b).

From the extracted feature map, in Fig. 9(c), the characteristics of the interface between concrete and soil medium, as well as the upper boundary characteristics of the cavity are obtained, but the lower boundary feature of the cavity have not been detected. In Fig. 9(d), the characteristics of the interface between concrete and soil media, as well as the upper and lower boundary characteristics of the aquifer, are extracted, and clutters are also included, which is labeled in Fig. 9(d). Therefore, the extracted feature map also needs to be classified, meanwhile, it is necessary to further detect the weak feature information of the target in the B-SCAN image, i.e., the lower boundary characteristic curve of the target.

Regarding the problem of feature recognition, it should be noted that, as shown in Figs. 9(c) and 9(d), semantic feature is embodied in the extracted feature map. To improve the accuracy of feature image recognition by the classification network, in this paper, the feature map is first mapped to the B-SCAN image to obtain more details, and then manual labeling is adopted to train the classification network.

Detection of weak characteristic of the target in B-SCAN image, for the convenience of presentation, the B-SCAN image with a cavity, as shown in Fig. 9(a), is used to analyze extracting the weak feature curve of the target in local space. Mapping transformation was conducted between the upper boundary feature map of the cavity, which is shown in Fig. 9(c), and the B-SCAN image, which is shown in Fig. 9(a), to determine the local region. The mapping transformation is denoted as $I_{\text{sub}}$ and is shown in the dotted line area in Fig. 10(a). In local space, the characteristic curve corresponding to the maximum amplitude is first detected. Then, the region corresponding to the negative peak characteristic curve of the scattering field at the medium interface is detected. Finally, the two adjacent regions are merged, and LS algorithm is used for curve fitting, and the fitting result is used as the final target lower boundary characteristic curve. In the real scene, the lower boundary characteristic curve of the target with weaker amplitude is susceptible to various interference factors, the accuracy of the detected lower boundary curve can be improved by region merging.

Fig. 10

Detection of weak feature curve of target in local space. (a) Local space of a B-SCAN image; (b) the detected positive peak curve; and (c) the detected negative peak curve.

The intensity of the lower boundary feature information of the target in the B-SCAN image is weak, and the characteristics of the geometrical structure of the hyperbola with downward opening is not obvious and relatively flat. The convolution kernel ${\mathit{C}}_{\mathrm{hori}}$ shown in Eq. (14) is designed, which can detect the features of horizontal linear structure well, and the convolution operation of ${\mathit{I}}_{\text{sub}}$ and ${\mathit{C}}_{\mathrm{hori}}$ is carried out to effectively segment the lower boundary feature curves with weak intensity and relatively flat structure feature

Eq. (14)

$${\mathit{C}}_{\mathrm{hori}}=\left[\begin{array}{ccccc}0& 0& -1& 0& 0\\ 0& 0& -2& 0& 0\\ 0& 0& +6& 0& 0\\ 0& 0& -2& 0& 0\\ 0& 0& -1& 0& 0\end{array}\right].$$

The characteristic figure ${\mathit{I}}_{\mathrm{seg}}$ is the result of the convolution operation between ${\mathit{I}}_{\text{sub}}$ and ${\mathit{C}}_{\mathrm{hori}}$, which can be calculated as

Eq. (15)

$${\mathit{I}}_{\mathrm{seg}}={\mathit{I}}_{\mathrm{sub}}*{\mathit{C}}_{\mathrm{hori}},$$

where ${\mathit{C}}_{\mathrm{hori}}$ convolved with the B-SCAN image in local space, and the signals without horizontal structure characteristics will be suppressed, while the flat characteristic curve of the target can be well detected. Then, the largest connected domain in ${\mathit{I}}_{\text{sub}}$ is taken as the characteristic curve of the target, the detection result is shown in Fig. 10(b). In local space, clutter is suppressed by ${\mathit{C}}_{\mathrm{hori}}$ and the segmentation of the target’s lower boundary characteristic curve, shown in Fig. 11.

Fig. 11

Segment the lower boundary feature curve of the target. (a) B-SCAN image with two cavities; (b) local space of the left cavity; and (c) the result of feature curve segmentation by convolution.

There are two cavities in the B-SCAN image, shown in Fig. 11(a). When detecting the lower boundary characteristic curve of the target in local space, it is affected by the upper feature hyperbola of the adjacent target. As it shown in Fig. 11(b), we can see that target and clutter are indistinguishable. Figure 11(c) is the result of ${\mathit{C}}_{\mathrm{hori}}$ convolved with the B-SCAN image in local space, which shows that the clutter is suppressed while the target features are enhanced.

Further, to detect the negative peak characteristic curve of medium interface in local space, ${\mathit{I}}_{\text{sub}}$ is first inverted and denoted as ${\mathit{I}}_{\mathrm{sub}}^{\prime }$ (${\mathit{I}}_{\mathrm{sub}}^{\prime }=-{\mathit{I}}_{\mathrm{sub}}$). Then, the convolution operation between ${\mathit{I}}_{\mathrm{sub}}^{\prime }$ and ${\mathit{C}}_{\mathrm{hori}}$ is carried out in the same process to detect the lower boundary characteristic curve of the target, i.e.,

Eq. (16)

$${\mathit{I}}_{\mathrm{seg}}^{\prime }={\mathit{I}}_{\mathrm{sub}}^{\prime }*{\mathit{C}}_{\mathrm{hori}}.$$

The maximum connected domain in ${\mathit{I}}_{\mathrm{seg}}^{\prime }$ is regarded as the characteristic curve of the negative peaks, and the detection result is shown in Fig. 10(c). It can be seen from Figs. 10(b) and 10(c), the convolution kernel ${\mathit{C}}_{\mathrm{hori}}$ can be effectively used to segment the positive and negative peak characteristic curves of the scattering field in local space.

5.

Analysis of Experimental Results

The test sample set has a total of 450 B-SCAN images, i.e., a single aquifer (S1), a single cavity (S2), two aquifers (S3), and two cavities (S4) are simulated with 100 images each, respectively. The other 50 B-SCAN images are obtained by setting spherical stones as clutters in the simulation scenario, sand pit experiments, and field measurements. Four simulation B-SCAN images, including target S1–S4, are selected for analysis, which is shown in Fig. 12(a).

Fig. 12

Target feature curve detection. (a) Sample images of experiment used for target feature extraction; (b) feature maps corresponding to the sample images; (c) the detected positive peak feature curve in local space; (d) the detected negative peak feature curve in local space; (e) combine the region of feature curves detected in Figs. (c) and (d); and (f) fitting results of target upper and lower boundary feature curves.

It can be seen from the extracted feature map shown in Fig. 12(b), the reflections of significant magnitude, and the feature map corresponding to S1 and S3 contains part of the lower boundary characteristic curve $C_L$ of the aquifer. The result shown in Fig. 12(c) is the extraction results of the characteristic curve (denoted by $C_P$) of the largest amplitude region in the local space. Compared with aquifer’s feature maps shown in Fig. 12(b), $C_L$ is consistent in spatial position, which also verifies the effectiveness of using the convolution kernel of Eq. (14) to extract feature curves in the local space. The result shown in Fig. 12(d) is the negative peak characteristic curve, denoted by $C_N$, of the scattered field at the medium interface in local space, compared with the result in Fig. 12(c), they are similar in structure and contiguous. Figure 12(e) is the result of the merged regions corresponding to $C_P$ and $C_N$, denoted as ${\mathrm{\Omega }}_{P\_N}$. Further, curve fitting based on LS algorithm is carried out for the upper boundary characteristic curve region (${\mathrm{\Omega }}_U$) and lower boundary characteristic curve region (${\mathrm{\Omega }}_{P\_N}$) of targets in Fig. 12(e), and the fitted results are shown in Fig. 12(f). It should be note that when the widths W1 and W2 of the upper and lower boundaries of the target are calculated according to Eq. (4), the accuracy is mainly determined by their fitting results. So, we need to determine the principal value interval (denoted as ${\mathrm{\Omega }}_{\text{up}}$ and ${\mathrm{\Omega }}_{\text{down}}$) of the target in ${\mathrm{\Omega }}_U$ and ${\mathrm{\Omega }}_{P\_N}$. Specifically, by setting thresholds $T_{\text{up}}$ and $T_{\text{down}}$, only regions with amplitude values greater than $T_{\text{up}}$ and $T_{\text{down}}$ are retained in ${\mathrm{\Omega }}_U$ and ${\mathrm{\Omega }}_{P\_N}$, respectively, which can be calculated as

Our approach focuses on effectively extracting the feature curve of the target in the B-SCAN image to further determine its buried depth and size. Compared with the upper boundary characteristic curve of the target, the amplitude of the lower boundary characteristic curve is usually small, and the detection is more difficult. Combining the physical mechanism of GPR, we proposed a method to detect the weak characteristic curve of the target in a local space, and obtained satisfactory detection results. The main reason for these results is that, on the one hand, according to the possible changes in phase information at the interface between the target and the background medium [the phase change due to the reflection coefficient is shown in Eq. (6)], the positive and negative peak curve regions detected in the local space are fused to fit the lower boundary characteristic curve of the target with a more complete structure [Fig. 12(f)]. On the other hand, the corresponding convolution kernel, as shown in Eq. (14), is designed with the structural characteristics of the target characteristic curve as a prior condition, and the target characteristic curve is segmented by convolution operation, which can largely suppress random interference signals with large amplitude (Fig. 11). Putting all these together, we can effectively detect the characteristic curve of the lower boundary of the target with weak intensity in B-SCAN image.

To further verify the effectiveness of the DCNN network designed in this paper for extracting the hyperbolic feature curve of the target, the B-SCAN images collected by adding clutters, where spherical stones with a radius of 0.1 to 0.25 m are set around the underground target as clutters, sand pit experiment, and field measurement are tested, respectively. The proposed method is compared with traditional histogram of oriented gradients (HOG) algorithm and Faster R-CNN algorithm. The result is shown in Fig. 13.

Fig. 13

GPR target detection experiment in complex environment. (a) Experimental pictures for testing; (b) the results of HOG algorithm; (c) the results of Faster R-CNN algorithm; and (d) the results of the proposed algorithm.

Figures 13(a1)–13(a3) are the spherical stone interference source set at different positions around the aquifer and the cavity target, Fig. 13(a4) is two iron balls set in the sand pit, Fig. 13(a5) is the measured steel bar in the concrete, and Fig. 13(a6) is the measured spherical limestone buried under the soil layer.³⁵ Figures 13(b)–13(d) are the detection results of HOG, Faster R-CNN, and the proposed algorithm, respectively. In Fig. 13(a1), the interference source is in the middle of two aquifers, there are false detection by HOG algorithm, and the Faster R-CNN algorithm misses the upper boundary characteristic curves of the target on the right, while the proposed algorithm effectively detects the upper boundary characteristic curves of the two targets. In Fig. 13(a2), the interference source is located directly above the aquifer, which has a serious influence on the structural features of the hyperbola of the target. Faster R-CNN algorithm failed, while HOG algorithm and the proposed algorithm detect part of the target feature curves. In Fig. 13(a3), the interference source is far to the upper right of the cavity, and all the three methods can effectively detect the upper boundary characteristic curve of the target. Figures 13(a4)–13(a6) are B-SCAN images obtained by sand pit experiments and field measurements, both HOG and Faster R-CNN algorithms failed to detect the target, with a high omission factor, and the DCNN network designed in this paper can well extract hyperbola features.

The above experimental results on GPR field data and simulation data with interference sources show that our approach is more robust. The main reason is that the proposed method is based on the analysis of the GPR physical model and obtained the prior information that the amplitude of the target scattering electric field in the echo data corresponds to the salient characteristic hyperbola in the B-SCAN image. On this context, CNN is designed to extract the high-level semantic features of B-SCAN images and used for the training of the classification network. Since the semantic features contain less interference information, the network model is easy to train and can obtain a higher recognition rate. Especially for the field B-SCAN images with complex background structure features shown in Figs. 13(a4)–13(a6), the influence of background interference can be greatly reduced by extracting semantic characteristics. Meanwhile, the GAN-based data enhancement strategy used in this paper also improves the accuracy of the network’s detection of field data.

To quantitatively analyze the accuracy of target detection results, for the fitted characteristic hyperbola, the comprehensive evaluation index $F$ is used for measurement, which can be calculated as

Fig. 14

Mean $F$ of detection result of target boundary feature curves. (a) The detection results of three algorithms for different types of targets and (b) mean $F$ of the three algorithms on the entire data set.

To objectively compare the detection results of the three algorithms, the detection of the target’s lower boundary characteristic curves are all carried out in local space. It can be seen from the statistical results in Fig. 14(a), in general, the mean $F$ of the three algorithms is that single target is higher than two targets, and aquifer is higher than cavity. This is mainly due to the fact that the completeness of the target characteristic curve detection result is a major factor, especially when the two targets are close to each other, the detection results are affected. Meanwhile, compared with air medium, the relative permittivity of the water medium is larger, and the speed of electromagnetic waves propagating in it is smaller, which corresponds to the longer two-way travel time interval in the B-SCAB image, and it is beneficial to reduce the mutual interference between the hyperbola of the upper and lower boundary of the target. As shown in Fig. 14(b), compared with the deep learning algorithm Faster R-CNN and the traditional machine learning algorithm HOG, the experiments on the test set show that our method has a high detection accuracy rate for both simulated data and field data, and the mean $F$ of the proposed algorithm is 93.23% and is 5.54% higher than Faster R-CNN algorithm. This also verifies the effectiveness of the proposed method, which is used to extract feature curve of the target in the GPR B-SCAN image.

On the basis of obtaining the upper and lower boundary characteristic curves of the target, the embedded depth $D$ and height $H$ of the target can be further calculated according to Eq. (4), and the result is measured by the accuracy rate, which is defined as

Fig. 15

Mean accuracy of the buried depth and height of the target.

According to Fig. 15, the proposed algorithm has the highest accuracy for $D$ and $H$ detection results, with an average accuracy of 94.57% and 93.72%, respectively. According to Eq. (4), the factors affecting the accuracy of $D$ and $H$ mainly include the two-way travel time and the relative permittivity of the medium. The accuracy of the two-way travel time is mainly determined by the accuracy of the vertex detection results of the upper and lower boundary characteristic curves of the target. For the relative permittivity of the medium, the estimation of ${\epsilon }_r$ is always rough, and the approximate of the inhomogeneous background medium will also bring errors. Some other factors will also affect the calculation results of $D$ and $H$, such as the frequency of the radar, the conductivity of the medium, the size and depth of the target, etc. The method in this paper processes the original echo data detected by GPR, and the cascade structure CNN designed to extract the salient features of B-SCAN images can effectively detect the target feature hyperbola under certain attenuation conditions. Note that it is constrained by the range resolution of GPR, and its calculation formula can be expressed as $D_r=c/4f_c\sqrt{{\epsilon }_r}$. It can be seen that under a certain $f_c$, the medium with slower detection wave speed can usually obtain a higher $D_r$. For a cavity target with a small height $H$, it will be difficult to effectively detect the lower boundary feature information due to the influence of the upper boundary characteristic curve with larger intensity information and the range resolution. This is also the main reason why aquifers have higher average accuracy than cavities at the target average height $H$.

Overall, the experimental results show that both HOG network and Faster R-CNN can obtain satisfactory detection results for hyperbolic features with complete structure. To effectively detect the features of the target in a complex background, CNN semantic feature extraction network is designed in this paper and used for the training of classification network, achieving better recognition results. Meanwhile, when the distance between the upper and lower boundary characteristic curves of the target is large, better detection results of the weak feature of the lower target can be obtained in the local space. Finally, for the limited amount of field GPR data, the GAN-based data enhancement method can better train the network model and obtain higher detection accuracy.

6.

Conclusion

Buried target detection from B-SCAN image is achievable using the CNN model based on GPR physical mechanism. Our approach consisting of three main components: automatic detection of salient feature information in B-SCAN image, GAN-based generation of realistic B-SCAN images, and detection of weak feature curve of target in local space. The proposed method has been evaluated in both simulated and field data, and has been compared with other CNN methods, i.e., traditional HOG algorithm and Faster R-CNN algorithm. The experimental results show that, on the comprehensive evaluation index $F$, the proposed method achieves the optimal result of 93.23%, which is 5.54% higher than the suboptimal result. In short, our method shows promising potential to use deep learning-based feature extraction to detect buried target from GPR field data. Notable, the CNN model designed in this paper focuses on the extraction of feature information from B-SCAN images, and combining the characteristics of GPR target feature information to design a more efficient classification and recognition network to improve the accuracy of target detection results is a work that needs further research in the future.