## 1.

## Introduction

Ground-penetrating radar (GPR) is a nondestructive method using electromagnetic radiation to locate shallow geological subsurface features and underground utilities buried in the ground. It has become a valuable tool in several applications,^{1}^{,}^{2} such as archaeological explorations, environmental engineering, and geological problems. The effective imaging of buried objects is a key part of GPR, and the efficiency and resolution of the imaging results are the measure of the imaging algorithm.^{3} The theories of the present imaging methods are based on diffraction tomography (DT), reverse time migration (RTM), range migration (RM), and back projection (BP). The principle of the DT algorithm is based on the first-order Born approximation which assumes that the buried object of interest is a weak scatterer.^{4} A few additional assumptions are also invoked during the process of DT derivation to simplify and linearize the nonlinear electric field integral equation. These assumptions incur a trade-off to the reconstruction of the buried objects especially for the practical usage when noise is present in the collected field data. Taking advantage of the multiple reflections in the propagation medium, the RTM algorithm allows high-resolution focusing.^{5} However, the number of transmitting and receiving antennas must be more than the number of scatterers in the medium. The RM algorithm can work well only when the imaging scene can be modeled as a single homogeneous medium.^{6} When the GPR antennas and buried objects are distributed in different media, the imaging result of the RM algorithm will be blurred or possibly not focused at all. The standard two-dimensional (2-D) depth migrations^{7} can recover the location and shape of buried objects with arbitrary precision, depending on the accuracy of the velocity model used. The BP algorithm is one of the most practical imaging methods because of its convenience and robustness, particularly when the imaging scene can be modeled as layered media.^{8} Based on the aforementioned theories, some of the improvements and optimization imaging methods have been advanced to distinguish the shape of buried objects in GPR imaging. A modified split-step method^{9} was applied to extract structural information from a complex synthetic data set as accurately as possible, based on the standard 2-D depth migrations. Furthermore, a synthetic aperture radar technique^{10} was implemented for GPR image reconstruction, which can recover the shape of buried objects. However, the present imaging methods depend too much on the application environment or prior knowledge of the medium being imaged, which limits the popularization and application of GPR technology.

The BP theory can accurately compensate the distortion of the wave path caused by GPR pulse signal when GPR pulse signal passes an interface of two media. It has become one of the potential GPR imaging algorithms. In order to recover the shape of buried objects in GPR imaging, the self-correlation back-projection (SBP) algorithm is proposed in this paper. The following improvements are made in the proposed algorithm. First, the valid echo information sequence of the imaging points is adaptively chosen by setting up a correlation threshold. Then, the imaging result is postprocessed by the depth energy compensation algorithm. Finally, the performance of the proposed algorithm is verified through serial contrast experiment.

## 2.

## Self-Correlation Back-Projection Algorithm

For the convenience of discussion, the 2-D imaging configuration of the GPR system is shown in Fig. 1. The scene is divided into two regions by $z=0$. The upper region is air with relative permittivity ${\u03f5}_{1}$ and conductivity ${\sigma}_{1}$. The lower region is a homogeneous medium with relative permittivity ${\u03f5}_{2}$ and conductivity ${\sigma}_{2}$. The GPR system works in a monostatic way. The antennas transmit and receive signals in each of the $M$ positions on a survey line $l$.

## 2.1.

### Two-Way Traveltime

In the 2-D imaging configuration shown in Fig. 1, the current concerned antenna position is represented by a black cube with sequence number $k$, the coordinates of which are $({x}_{k},-h)$. For a given point $A$ with coordinates $({x}_{0},{z}_{0})$ in the imaging area, the transmitting signal travels from $({x}_{k},-h)$ to $({x}_{0},{z}_{0})$, with a turning at the inflection point $({x}_{r},0)$, and returns along the same path in reverse direction. According to Snell’s law, the geometry between the incidence angle ${\theta}_{1}$ and the refraction angle ${\theta}_{2}$ satisfies

## (1)

$$\frac{\mathrm{sin}\text{\hspace{0.17em}}{\theta}_{1}}{\mathrm{sin}\text{\hspace{0.17em}}{\theta}_{2}}=\frac{c/\sqrt{{\u03f5}_{1}}}{c/\sqrt{{\u03f5}_{2}}}=\sqrt{\frac{{\u03f5}_{2}}{{\u03f5}_{1}}},$$## (2)

$$\frac{{({x}_{r}-{x}_{k})}^{2}}{{h}^{2}+{({x}_{r}-{x}_{k})}^{2}}/\frac{{({x}_{0}-{x}_{r})}^{2}}{{z}_{0}^{2}+{({x}_{0}-{x}_{r})}^{2}}=\sqrt{\frac{{\u03f5}_{2}}{{\u03f5}_{1}}}.$$Then, ${x}_{r}$ can be obtained by solving Eq. (2). If ${x}_{r}$ is known, the two-way traveltime ${\tau}_{A,k}$ from the imaging point $A$ to the $k$’th antenna position can be given by

## 2.2.

### Self-Correlation Back-Projection Algorithm Developing Steps

For a given point $A$ with coordinates $({x}_{0},{z}_{0})$ in the imaging area shown in Fig. 1, its projection point on the survey line $l$ is defined as $B$ with coordinates $({x}_{0},-h)$. The $i$’th antenna position with coordinates $({x}_{i},-h)$ is the nearest to $B$ on the survey line $l$. At the $i$’th antenna position, the A-scan of GPR data can be given by

## (4)

$$r(p,q)=\frac{\sum _{i=1}^{n}({p}_{i}-\overline{p})({q}_{i}-\overline{q})}{\sqrt{\sum _{i=1}^{n}{({p}_{i}-\overline{p})}^{2}}\times \sqrt{\sum _{i=1}^{n}{({q}_{i}-\overline{q})}^{2}}},$$Step 1: The two-way traveltime ${\tau}_{A,i}$ from the imaging point $A$ to the $i$’th antenna position is computed by Eq. (3). Both the valid echo information ${u}_{A,i}$ and the target related sequence ${C}_{i}$ are chosen from ${S}_{i}$, where the length of ${C}_{i}$ is $N$.

In vector ${S}_{i}$, ${T}_{1}$ is defined as a sampling point nearest to ${\tau}_{A,i}$. The valid echo information corresponding to ${\tau}_{A,i}$ is given by ${u}_{A,i}=|{s}_{i,{T}_{1}}|$, where $|\xb7|$ is the absolute value operator. Extracting vector ${S}_{i}$ sequentially from ${T}_{1}$, we obtain $N$ sampling points, namely the target-related sequence ${C}_{i}$. Each A-scan signal is the stack of multiple reflection echo signals. For a given imaging point, the echo information extracted by the BP algorithm from each A-scan signal may be the stack of the echo information of many imaging points. For this reason, when $N$ exceeds the permitted values, ${C}_{i}$ may contain additional echo information. Otherwise, when $N$ is too small, ${C}_{i}$ is unstable and cannot characterize the features of the reflection echo signals from the same imaging point. A Gaussian pulse, wideband Ricker wavelet, differentiated Gaussian pulse, and Blackman–Harris window function are usually used as the pulse source signal of GPR. If the pulse duration is $T$, the pulse source signal is symmetrical about its midpoint $T/2$. On the basis of the aforementioned considerations, $N$ is defined as $N=\lceil T/(2\xb7\mathrm{\Delta}t)\rfloor $, where $\lceil \xb7\rfloor $ can round the element $\xb7$ to the nearest integer, and $\mathrm{\Delta}t$ is the time interval of GPR sampling.

Step 2: For a given imaging point $A$, its valid echo information is chosen adaptively to the left starting from the $i-1$th antenna position, which is expressed as

$${U}_{A,\text{left}}=[\begin{array}{cccc}{u}_{A,i-L}& {u}_{A,i-L+1}& \cdots & {u}_{A,i}\end{array}]\mathrm{.}$$By using step 1, ${\tau}_{A,i-1}$, ${u}_{A,i-1}$, and ${C}_{i-1}$ are computed. The correlation coefficient $r({C}_{i},{C}_{i-1})$ between ${C}_{i}$ and ${C}_{i-1}$ is computed by Eq. (4). If $r({C}_{i},{C}_{i-1})>\psi $, ${u}_{A,i-1}$ is a valid echo information of the imaging point $A$, and then ${u}_{A,i-2}$ will be judged. Otherwise, the search to the left is stopped. $\psi $ is the correlation threshold related to the specific detection environment. The stronger the background noise, the smaller the $\psi $. The more appropriate $\psi $ is the more accurate the chosen valid echo information sequence of the object and better resolution imaging results will be obtained. $\psi $ is usually adjusted by known objects in the specific detection environment before imaging processing of GPR signals. It usually ranges from 0.85 to 0.95 per Ref. 2. Its value is 0.9 in this paper.

Step 3: For a given imaging point $A$, its valid echo information is chosen adaptively to the right starting from the $i+1$th antenna position, which is expressed as ${U}_{A,\text{right}}=[\begin{array}{cccc}{u}_{A,i}& {u}_{A,i+1}& \cdots & {u}_{A,i+R}\end{array}]$. It has the same process of problem solving as step 2.

Step 4: Calculating the temporary value ${E}_{A}^{*}$ of $A$ for the imaging result.

The valid echo information sequence ${U}_{A}$ of $A$ can be obtained by merging ${U}_{A,\text{left}}$ and ${U}_{A,\text{right}}$, which is formulated as

$${U}_{A}={U}_{A,\text{left}}\cup {U}_{A,\text{right}}=[\begin{array}{cccc}{u}_{A,i-L}& {u}_{A,i-L+1}& \cdots & \begin{array}{cccc}{u}_{A,i-1}& {u}_{A,i}& {u}_{A,i+1}& \begin{array}{cc}\cdots & {u}_{A,i+R}\end{array}\end{array}\end{array}].$$To suppress the artifacts, ${E}_{A}^{*}$ is calculated by adding an additional cross-correlation procedure,

^{8}which is given byStep 5: When the high-frequency electromagnetic signal propagates in the medium, the signal declines in an approximately exponential function. To improve the imaging resolution of deeper underground objects, the imaging result in step 4 is amended by using a depth energy compensation algorithm.

^{11}The amended imaging result of point $A$ can be formulated aswhere $\eta (z)=\mathrm{exp}(\alpha \xb7z)$ is the depth energy compensation coefficient, and $\alpha $ is determined by the medium dielectric property.## (6)

$${E}_{A}=\eta (z)\xb7{E}_{A}^{*}{|}_{z={z}_{0}}=\mathrm{exp}(\alpha \xb7{z}_{0})\xb7{E}_{A}^{*},$$

If the received GPR data have been amplified, the depth energy compensation algorithm will not be necessary. Otherwise, step 5 should be performed.

The aforementioned steps will be repeated until all the points in the imaging scene are considered.

## 3.

## Experimental Results and Analysis

To verify the performance of the proposed SBP algorithm, we design a synthetic and a practical field experiments for the objects imaging of complex shape. In addition, a serial contrast experiments are performed with the classic BP, the fast back projection (FBP)^{8}, the SBP, and the 2-D depth migration algorithm with straight ray assumption.^{12} To solve Eq. (2), the binary search algorithm is used by both the classic BP and SBP algorithms; the approximate algorithm proposed in Ref. 8 is used by the FBP algorithm. We evaluate the imaging error by the absolute cumulative error (ACE), which is defined as

## 3.1.

### Tests with Synthetic Data

According to Fig. 1, the geoelectric model of synthetic GPR data is formulated. The antennas are located 1 m above the surface of the surveyed structure. The scene is divided into two regions by $z=0$. The upper region is air with relative permittivity 1 and conductivity $0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{S}/\mathrm{m}$. The lower region is homogeneous medium with relative permittivity 16 and conductivity $0.001\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{S}/\mathrm{m}$. The buried objects are shown in Fig. 2 by the blue curve with relative permittivity 9. The GPR wave fields of the geoelectric model are simulated by using a finite-difference time-domain^{13} solution of the Maxwell’s equations. In the synthetic experiment, a Ricker wavelet is taken as a pulse source with a center frequency of 100 MHz. The B-scans of the synthetic GPR data are shown in Fig. 3. By using the mean removal algorithm, we can preprocess the synthetic GPR data to suppress the interference partially from the terrain echo and the coupled wave between receiving and transmitting antennas. The normalized imaging results of the classic BP, the FBP, and the SBP are shown in Figs. 4, 5, and 6, respectively. The aforementioned imaging results are shown in Table 1. The velocity model of the 2-D depth migration algorithm is shown in Fig. 7, where the electromagnetic wave travels at $3\times {10}^{8}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}/\mathrm{s}$ in air and $7.5\times {10}^{7}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}/\mathrm{s}$ in earth, respectively. The normalized imaging results of the 2-D depth migration algorithm are shown in Fig. 8. The theoretical distributions of buried objects are labeled with the white circle curve, which is convenient to analyze the imaging results.

## Table 1

Performance comparisons table about different imaging methods.

Imaging methods | y-shape object | z-shape object | ||
---|---|---|---|---|

Absolute cumulative errors (ACEs) | Running times (s) | ACEs | Running times (s) | |

Classic back-projection | 616.4 | 52.92 | 1017.4 | 52.95 |

FBP | 250.3 | 3.71 | 471.1 | 3.72 |

SBP | 28.9 | 2.53 | 51.6 | 2.54 |

The experimental results show that the objects shape cannot be distinguished from Fig. 3. It may be recovered by the GPR signal-processing techniques. Due to the fact that buried objects have complex shapes in this experiment, each A-scan signal is the stack of multiple reflection echo signals. For a given imaging point, the echo information extracted by the BP algorithm from each A-scan signal may be the stack of the echo information of many imaging points. It reduces the imaging accuracy of the classic BP and FBP algorithms. For the same B-scan in GPR imaging, at the same condition of imaging resolution, the larger the number of the sensing locations, the longer the running time of both the BP and FBP algorithms. The SBP algorithm can adaptively choose the valid echo information sequence, with its running time related to the specific B-scan. The larger the number of the valid echo information sequence, the longer the running time of the SBP algorithm. However, the SBP algorithm usually has faster calculation speed because the valid echo information sequence is generally located in the multiple nearest neighbors of the imaging point. The 2-D depth migration algorithm can distinguish the objects shape effectively. However, the locations of the objects are moved upwards compared with their theoretical distribution. The GPR antennas and objects are distributed in different media. The precision of velocity model is reduced by the air between the GPR antennas and the earth. It is the main reason for the aforementioned imaging error. The experimental results illustrate that the proposed SBP algorithm is superior to the existing BP algorithms in terms of computing speed and imaging accuracy. Compared with the 2-D depth migration algorithm, the proposed SBP algorithm has a significant advantage in providing a rough outline of buried objects without prior knowledge of the velocity distribution.

## 3.2.

### Field Test Case

The experimental data are collected at the Georgia Institute of Technology^{14} and is publicly available in Ref. 15 in MATLAB format files, and the field map in the detecting location is shown in Fig. 9. The data consist of different burial and no-object scenarios, and are taken with a multistatic stepped-frequency continuous-wave GPR. The GPR is scanned over a $1.8\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}\times 1.8\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ region at a constant height above the surface of the ground. The scan region is discretized into a grid of 91 points. At each scan position, GPR takes 401 equally spaced frequency measurements from 60 MHz to 8.06 GHz with 20-MHz increments. The pulse source is a differentiated Gaussian pulse with a center frequency of 2.5 GHz. The GPR B-scan is shown in Fig. 10, which is selected from the aforementioned experimental data. There are five buried objects at $-45$, $-20$, 0, 20, and 45 in cross-range ($x$) dimension and $y=-40$. The distance between the transmitting and receiving antennas is 12 cm. The height of antennas is 27.8 cm above the surface of the surveyed structure. Each A-scan is acquired every 2 cm. The electromagnetic wave travels at $2.998\times {10}^{8}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}/\mathrm{s}$ in air and $1.5\times {10}^{8}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}/\mathrm{s}$ in sand.^{16} The field GPR data are preprocessed by the removing mean value method to suppress the interference partially from the terrain echo and the coupled wave between receiving and transmitting antennas. The normalized imaging result of the SBP algorithm is shown in Fig. 11. The theoretical distributions of buried objects are labeled with the black curve, which is convenient to analyze the imaging results.

The experimental result shows that the proposed SBP algorithm can effectively recover the multiple objects’ position information. The calculated ACE for Fig. 11 is 426.5, and the running time of the SBP algorithm is 3.5 s. It further validates the effectiveness of the proposed SBP algorithm.

## 4.

## Conclusion

Based on the existing BP algorithms, the SBP algorithm is proposed in this paper, which can reconstruct the shape of buried objects in GPR imaging. By setting up the correlation threshold, the SBP algorithm can adaptively choose the valid echo information sequence of the imaging points, which improves the image resolution. Because the valid echo information sequence is generally located in the multiple nearest neighbors of the imaging point, the SBP algorithm usually has a faster calculation speed. In addition, the imaging result is postprocessed by the depth energy compensation algorithm to improve the imaging resolution. The experimental results show that the SBP algorithm is superior to the classic BP and FBP algorithms in terms of computing speed and imaging accuracy. It has a significant advantage in providing a rough outline of buried objects without prior knowledge of the velocity distribution. The application of the proposed SBP algorithm is valuable in the imaging of underground objects with fast speed and high quality.

## Acknowledgments

This work was funded by the National Natural Science Foundation of China (Nos. 61102115, 61362020, and 61371186), the Natural Science Foundation of Guangxi Province (Nos. 2012GXNSFBA053177, 2013GXNSFAA019327, and 2013GXNSFFA019004), and the Guangxi Experiment Center of Information Science, Guilin University of Electronic Technology (No. 20130112).

## References

## Biography

**Hairu Zhang** received his BEng degree from Wuhan University of Science and Technology, Wuhan, China, in 2009 and obtained his MEng degree from Guilin University of Electronic Technology, Guilin, China, in 2012. Currently, he is a PhD student at Xidian University, Xi’an, China. His fields of interest include signal processing for ground-penetrating radar, computational intelligence, and adaptive signal processing.

**Shan Ouyang** received his BEng degree from Guilin University of Electronic Technology, Guilin, China, in 1986 and obtained his MEng and PhD degrees in electronic engineering from Xidian University, Xi’an, China, in 1992 and 2000, respectively. He received the National Excellent Doctoral Dissertation of China in 2002. Currently, he is a professor and PhD supervisor at Guilin University of Electronic Technology. His fields of interest include signal processing for communications and radar.

**Guofu Wang** received his BEng degree from Wuhan University of Technology, Wuhan, China, in 2002 and obtained his MEng and PhD degrees from Xi’an Institute Optics and Precision Mechanics of CAS, Xi’an, China, in 2005 and 2007, respectively. Currently, he is a professor at Guilin University of Electronic Technology. His fields of interest include signal processing for ground-penetrating radar and adaptive signal processing.

**Jingjing Li** received his BEng degree from Henan Polytechnic University, Jiaozuo, China, in 2010 and obtained his MEng degree from Guilin University of Electronic Technology, Guilin, China, in 2014. Currently, she is a PhD student at Guilin University of Electronic Technology, Guilin, China. Her fields of interest include signal processing for radar.