Back-projection algorithm based on self-correlation for ground-penetrating radar imaging

Abstract. In ground-penetrating radar imaging, the classic back-projection (BP) algorithm has an excellent reputation for imaging in layered media with convenience and robustness. However, it is time-consuming and generates many artifacts, which have adverse effects on detection and recognition. A self-correlation back-projection (SBP) algorithm is proposed, which is fast imaging and can distinguish the object’s shape. It improves the existing BP algorithms in the following aspects. First, the reflection echo signals of a specific imaging point obtained from its nearest exploration point have high correlation with the one from its multiple nearest neighbors. By setting up a correlation threshold, the valid echo information sequence of the imaging points can be adaptively chosen, which enables the SBP algorithm to have faster calculation speed and better resolution. Then, the imaging result is amended by using a depth energy compensation algorithm. It can improve the imaging resolution of the deep underground objects. The experimental results show that the proposed SBP algorithm is superior to the existing BP algorithms in terms of computing speed and imaging accuracy, which can effectively recover objects with complex shapes. It has a significant advantage in providing a rough outline of buried objects without prior knowledge of the velocity distribution.


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. 3The 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. 4A 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. 5However, 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. 6When 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. 8Based 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.

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 ε 1 and conductivity σ 1 .The lower region is a homogeneous medium with relative permittivity ε 2 and conductivity σ 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.

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 θ 1 and the refraction angle θ 2 satisfies Fig. 1 Imaging configuration.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 7 3 5 where c is the velocity of the electromagnetic wave in free space.Equation (1) can be turned into Eq.( 2) based on the coordinates in Fig. 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 6 7 1 Then, x r can be obtained by solving Eq. ( 2).If x r is known, the two-way traveltime τ A;k from the imaging point A to the k'th antenna position can be given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 5 9 6

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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 2 ; 1 1 6 ; 4 7 2 where L is the number of sampling points; t j is the sampling instant of s i;j .The correlation coefficient rðp; qÞ between the vector p and the vector q is defined as ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 4 1 6 rðp; qÞ where n is the length of p and q, and p and q are the mean values of p and q, respectively.In the following, the SBP algorithm will be described for GPR imaging.
Step 1: The two-way traveltime τ 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 τ A;i .The valid echo information corresponding to τ A;i is given by u A;i ¼ js i;T 1 j, where j • j 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 ¼ dT∕ð2 • ΔtÞc, where d•c can round the element • to the nearest integer, and Δ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 − 1th antenna position, which is expressed as ; t e m p : i n t r a l i n k -; s e c 2 . 2 ; 1 1 6 ; 1 1 7 By using step 1, τ A;i−1 , u A;i−1 , and 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.ψ is the correlation threshold related to the specific detection environment.The stronger the background noise, the smaller the ψ.The more appropriate ψ is the more accurate the chosen valid echo information sequence of the object and better resolution imaging results will be obtained.ψ 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 þ 1th antenna position, which is expressed as 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;left and U A;right , which is formulated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 2 ; 1 1 6 ; 5 7 2 To suppress the artifacts, E Ã A is calculated by adding an additional cross-correlation procedure, 8 which is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 4 9 6 Step 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. 11The amended imaging result of point A can be formulated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 3 8 9 where ηðzÞ ¼ expðα • zÞ is the depth energy compensation coefficient, and α is determined by the medium dielectric property.
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.

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. 12To 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 1 5 2 where B and B Ã are the normalized matrices of the objects position matrix and imaging results matrix, and b i;j and b Ã i;j are the elements in the matrices B and B Ã , respectively.In this paper, all the programs are run by MATLAB 7.1 on the same hardware.

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 S∕m.The lower region is homogeneous medium with relative permittivity 16 and conductivity 0.001 S∕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 × 10 8 m∕s in air and 7.5 × 10 7 m∕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.
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      Fig. 9 Field maps in the detecting location. 16nd 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.

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 m × 1.8 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 × 10 8 m∕s in air and 1.5 × 10 8 m∕s in sand. 16he 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.

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.

Fig. 3 B
Fig. 3 B-scans of the synthetic GPR data: (a) y -shape object and (b) z-shape object.

Fig. 4
Fig. 4 Imaging result of the classic back-projection algorithm: (a) y -shape object and (b) z-shape object.

Fig. 5
Fig. 5 Imaging result of the FBP algorithm: (a) y -shape object and (b) z-shape object.

Fig. 6
Fig. 6 Imaging result of the self-correlation back-projection (SBP) algorithm: (a) y -shape object and (b) z-shape object.
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

Fig. 8
Fig. 8 Imaging result of the 2-D depth migration algorithm: (a) y -shape object and (b) z-shape object.

Fig. 10 B
Fig. 10 B-scan of the field GPR data.

Fig. 11
Fig. 11 Imaging result of the SBP algorithm.

Table 1
Performance comparisons table about different imaging methods.