A Multi-Scale Weighted Back Projection Imaging Technique for Ground Penetrating Radar Applications

In this paper, we propose a new ground penetrating radar (GPR) imaging technique based on multi-scale weighted back projection (BP) processing. Firstly, the whole imaging region is discretized by large scale and low-resolution imaging result is obtained by using traditional BP imaging technique. Secondly, the potential targets regions (PTR) are delineated from low-resolution imaging result by using intensity detection method. In the PTR, small scale discretization is implemented and higher resolution imaging result is obtained by using weighted BP imaging technique. A weight factor is designed by analyzing the statistical characteristics of scattering data on the time-delay curve. The above “discretization-imaging-PTR delineation” processing continues until the imaging resolution reaches the specified requirement. In the multi-scale imaging result, the resolution in other regions is not as high as that in PTR. This algorithm can get higher resolution imaging results with much lower computation compared with traditional BP imaging algorithm. The simulation of this algorithm is processed and experimental results validate the feasibility of this method.


Introduction
Wide band radar has been widely used in subsurface detection applications [1,2].The application of GPR in the detection and discrimination of underground objects, such as buried landmines and UXO (Unexploded Ordnance), is currently receiving much attention [2][3][4] By using synthetic aperture technique, GPR can obtain high resolution imaging results and it's very useful for geophysical prospecting, landmines reclamation, civil engineering applications, and near-surface probing [1][2][3][4].
In recent years, different GPR imaging techniques are introduced based on different scattering models.Based on electromagnetic scattering model, subsurface targets are modelled as Δε( → r ) and the imaging algorithm can be established from electromagnetic scattering integral equation [5][6][7][8][9][10].
By posing the imaging problem in the wavelet transform domain, Eric et al. and Chiappinelli et al. presented multiresolution imaging techniques and showed some inversion results [5,6].Based on forward scattering iterations, the inversion results can be obtained by solving a optimization problem, accompanied with large computation.Caorsi et al. introduced an approach for imaging inhomogeneous cylinders of arbitrary cross sections based on iterative multiresolution processing [7].By including in an efficient way all the available a priori information, the author established a multiscale inversion framework.The problems such as the occurrence of local minima and improper selection of stopping criterion may make the inversion algorithm invalid.Catapano et al. presented joint inversion of unknown embedded targets' permittivity and subsurface roughness profile by solving a local minimization problem [8].The approach proposed in [8] provides a multi-resolution reconstruction of the region under test by using wavelet Haar as basis functions to represent the unknown electric contrast.By doing so, one can accommodate the sought coefficients in a non-uniform way into investigated domain and this chance is useful for at least two different reasons.First, one can focus on the regions wherein the scatterers are actually located or a more variable permittivity profile is expected.Second, in some cases such as the half-space one, a variable degree of resolution is expected with different depths, so that it comes out to be useful to look for a large number of coefficients in the shallow region while using a lower resolution for the deeper part of the region under test.Moreover, in [9] a two steps procedure is proposed, in which the scatterers support is reconstructed by means of a qualitative imaging approach firstly and then such an information is exploited to reduce the investigated region wherein the quantitative imaging is performed.The imaging methods presented in references [7][8][9] all belong to the class of quantitative imaging strategies, that is to say, they aim at reconstructing the dielectric permittivity and electric conductivity of the targets and provide information on their geometrical features.Soldovieri et al. faced the problem of embedded rebar localization and transformed it to a linear inverse problem by exploiting the Born scattering model [10].This inversion method falls within the framework of the sparse minimization and accounts for the sparse nature of the scatterers in the investigated domain by exploiting a distributional representation of the unknown function.Based on Born approximation, the relationship between the scattering data and the unknown parameters is established under the assumption of neglection of mutual scattering components.Therefore, this method does not provide a quantitative reconstruction of the investigated scenario but it aims at localizing unknown targets and estimating their size and shape.
Another kind of imaging techniques are based on point scattering model [11][12][13][14].Fauzia et al. presented a new idea for sparse target imaging based on compressed sensing technique and showed good imaging results [11].The main advantage of this method lies in not demanding requirements of data collection.As shown in [11], it's possible to reconstruct the corresponding range profile from a reduced set of randomly selected frequencies.Based on compressed sensing framework, this imaging algorithm needs the construction of dictionary matrix and effective solving technique for optimazation problem.
As a representative imaging method corresponding to point scattering model, traditional GPR BP algorithm is suitable for subsurface imaging of targets located in lossy layer media.Lei et al. analyzed GPR wave's radiation characteristic in subsurface layered medium and provided the corresponding BP imaging algorithm based on basic electromagnetic radiation theory [12].Zhou et al. proposed a novel artifact suppression technique and provided a modified BP imaging algorithm [13].Lei et al. analyzed the statistical feature of scattering data at time-delay curve and proposed a new windowed and weighted BP imaging algorithm [14].Based on "delay and sum" operation, all these algorithms are single-scale and thus take heavy computation burden.
To improve the efficiency and resolution of BP imaging, a multi-scale weighted BP imaging algorithm is studied in this paper.This paper is organized as follows.In Section 2, the procedure of traditional BP imaging algorithm is described and a multi-scale weighted BP imaging algorithm is presented.In Section 3, the simulation experiments and on-site experiment of multi-scale weighted BP imaging algorithm are processed.The imaging results of multi-targets are presented and contrasted with traditional BP imaging results at the aspects of resolution and computation burden.

Traditional BP Imaging Algorithm
The scene of a mono-static GPR measurement setup is shown in Figure 1.An impulse signal is transmitted through GPR antenna to underground region.The antenna is positioned at point (x a , −h a ) and synthesizes an aperture on a line parallel to the x axis at a distance h a .The main lobe width of the antenna is 2β and the whole scanning length of targets located at depth of h b is x b − x a .As the T/R (transmitter/receiver) antenna pair move along the aperture line z = −h a with an interval of Δx, multi scattering signals e 1 (x, t) at each focal aperture point x p = x a + (p − 1)Δx, p=1, …, P can be recorded.The total scan number is P=1+(x b − x a )/Δx.The purpose of GPR imaging algorithm is to obtain e 2 (x, z) from the original radargram e 1 (x, t), that is to say, by process e 1 (x, t) profile, the scattering density profile of the scanning region can be obtained.
Traditional BP imaging algorithm based on "delay and sum" operation is a well-known method [13,14].The basic idea of this imaging algorithm is to sum all data coherently at a time-delay curve of one focal point (x n , z m ) and repeat for all points in the ROI (range of interest).Before imaging, the ROI needs to be discretized.For example, a 2-dimensional ROI can be divided into M × N pixels.The purpose of BP imaging is to obtain the scattering density of each grid, shown as e 2 (x n , z m ).At each focal point (x n , z m ), all time-shifted responses are coherently summed up.The relevant equation is shown in Equation (1).
where e 1 ′(x, t) is the radargram after pre-processing procedure, and τ m,n,p is the travel time from the transmitter to the focal point and back to the receiver: where c = 3 × 10 8 m/s is the speed of light, v represents the propagation velocity of electromagnetic wave in the soil, R a and R b are the distances from the refraction points to the T/R pair point (x p ,z = −h a ) and focal point (x n , z m ), respectively.Based on "delay and sum" operation, the BP imaging algorithm is easy to realize.But it needs M × N times of "delay and sum" processing and the computation burden is very heavy.

Multi-Scale Processing
In GPR detection scenario, the potential targets are often located at small region of the whole ROI, that is to say, the targets have spatial sparse feature.This means that only a small region of ROI needs high-resolution imaging, the other region needn't high-resolution imaging.In traditional BP imaging procedure, the whole ROI is divided by single-scale and all grids are considered as the same.This division method brings much computation and large parts of the computation are useless.
By using multi-scale division, the whole ROI, denoted as Ω 0 , is divided into large-scale pixels and low-resolution imaging result is obtained firstly.Suppose the size of large-scale pixel is Δz 1 × Δx 1 , the whole ROI is divided into M 1 × N 1 pixels.The imaging result of M 1 × N 1 pixels, denoted as e I (x, z)| Ω0 , can be obtained by using Equation (1).
Then by using intensity detection technique, the potential target regions are extracted and delineated from low-resolution imaging result.The whole ROI can be divided into two parts, one is PTR, denoted as Ω I-obj , the other is rest, denoted as Ω I-non-obj .In Ω I-obj , smaller scale discretization is implemented.Suppose the size of discrete grid is Δz 2 × Δx 2 , higher resolution imaging result can be obtained by using Equation (1), denoted as e II (x, z)| ΩI-obj .The imaging result of whole ROI Ω 0 can be expressed as the combination of e II (x, z)| ΩI-obj and e I (x, z)| ΩI-non-obj .After several times iterations, the size of imaging grid in PTR becomes smaller.When it reaches the resolution requirement, the iteration terminated.The imaging result of whole ROI Ω 0 may has the combination form like {e S (x, z)| ΩI-non-obj ; …; e II (x, z)| ΩI-non-obj ; e I (x, z)| ΩI-non-obj }.
In the multi-scale BP imaging procedure, the whole region is divided into discrete grids with different sizes hierarchal.Only regions where the potential targets located at are divided into pixels with small size and other regions are divided into pixels with large size.By using this multi-scale technique, the total computation is decreased greatly while the resolution of PTR is guaranteed.
The selection criteria of key parameters in multi-scale BP imaging algorithm are discussed below.
(1) Initial sampling grid Based on above analysis, the min-max range of initial sampling grid is listed below: where T represents the sampling points of each A-scan.When we let M 1 = T, N 1 = P, the imaging algorithm is based on single-scale processing actually, that is to say, it degenerated into traditional BP imaging algorithm.When we let M 1 = 2, N 1 = 2, the imaging result e I (x, z)| Ω0 only contains 4 pixels.The selection of initial sampling grid depends on specific GPR detection scenario and radargram.On occasions of subsurface sparse targets detection such as pipe detection,tunnel detection, UXO (unexploded ordance) detection, underground grave detection and etc., the selection of initial sampling grid follows the following general guideline: M 1 = N 1 = round(P/α), where α = 5~10, the operator round(a) returns a value to the nearest integers of a.
On occasions of dense targets detection such as reinforced bar mesh detection, M 1 = N 1 = P is preferred.
(2) PTR delination A intensity threshold detection method is used for PTR delineation.Suppose the current imaging reuslt is described as e s (x, z)| Ω0 , we select the following region as PTR: where k = 0.4~0.6.
(3) Hierarchical refinement coefficient If the grids number of current PTR is M s × M s and the imaging resolution in PTR haven't meet the required resolution, thus the imaging pixel in a new PTR after PTR extraction operation needs to be refined.The hierarchical refinement coefficient follows the following general guideline:

) Iteration stop condition
With the iteration procedure running, the size of imaging pixels in PTR becomes smaller.Suppose the size of current discrete grid is Δz i × Δx i , the iteration stop condition is shown below:

Weighted Factor
In addition, the statistic features of scattering data at each time curve can be used to design a weighted factor to improve the imaging resolution.By designing a new additional weighting factor w(x n , z m ), the imaging quality can be improved significantly.The weighting factor w can be interpreted as a mask to the traditional BP imaging results.The selection of w is based on the fact that at such points (x n , z m ) where there are real targets, the 1-dimensional scattering data has good correlation.It can be calculated in three steps.Firstly, for each focal point, the scattering data e 1 ′(x p , t = τ m,n,p ) at time-delay curve need to be extracted from the original radargram.Next, the mean value and standard deviation is calculated as shown below.
  , In the third step, the weighting factor w(x n , z m ) is calculated as By multiplying the above weighted factor at each pixel, higher resolution imaging results can be obtained as the following formula:

Experiments and Imaging Results
In this section, we apply multi-scale weighted BP imaging algorithm to some simulation data and experimental data.Simulation data are obtained by using FDTD (Finite Difference Time Domain) simulation technique [15] while real GPR data are obtained from on-site experimental collection.
GPR simulation procedure is briefly described below.Firstly, a 2-dimensional scanning scenarios is established virtually.A cube, with infinite length in y direction and a defined rectangular in oxz plane as its cross section, is designated as background medium.Subsurface embedded targets, also with infinite length in y direction, are buried at known locations in this medium.The electrical properties of both background medium and embedded targets are designated.Secondly, GPR parameters, such as exciting signal, antenna form, scanning procedure and antenna height are also designated.Suppose GPR antenna scans the subsurface region along the x direction.Thirdly, the whole 2-D region is discretized into small grids and time variable is also discretized, all according to FDTD discretization criterion.And then we use 2-D FDTD program to implement GPR forward simulation and obtain GPR echo profile.The purpose of this process is to simulate real-world GPR scanning procedure and obtain ground-penetrating radargram.

Example 1
Figure 2a shows the original GPR detection scenario of single target.One metal bar with infinite length is located at point (1.0, 0.7).The relative permittivity of the background medium is 16.The transmitting antenna is located at point (1.0, 0) and is excited by a Ricker wavelet with center frequency of 400 MHz.The receiving antenna is moved from point (0.2, 0) to point (1.8, 0), with the interval space of 0.02 m.Thus, the total aperture number is 81.At each aperture point, the scattering signal is received by the receiving antenna.Each A-scan contains 4240 sampling points in the total 40 ns time window.After finishing the scan procedure along the aperture line, a B-scan profile of the scanning region can be obtained.Before the imaging processing, a preprocessing procedure such as direct wave subtracting or inverse filtering is needed [1]. Figure 2b shows the B-scan profile after preprocessing.The iteration imaging results and PTR delineation results in the first 3 rounds of iteration are shown in Figure 2c-f.The iteration parameters in the imaging processing is shown in Table 1.

GPR Scanning
After 2nd iteration, Δx 3 < Δx, the iteration processing stopped.The imaging results of multi-scale weighted BP and traditional BP imaging algorithm are shown in Figure 2g,h, respectively.
The processing of hierarchal division and multi-scale imaging is shown clearly.By contrast with traditional BP imaging result, it can be shown that the imaging result in Figure 2g shows lower side lobe level and higher focusing quality than that of imaging result in Figure 2h.In order to quantitatively assess the imaging results, a focusing parameter is used here [16].The focusing parameter of a 2D imaging result O(x, z) can be calculated as below: Here we only calculate focusing parameters of image pixels corresponding to the same PTR shown in Figure 2f.The focusing parameters of the above two imaging results defined in the same PTR are 0.131 and 0.0167 respectively.The focusing quality of the imaging result improved above 680%.Both imaging results are obtained via a PC with CPU Pentium(R) Dual-Core 2.00 GHz and RAM 2.0 GB.The time requirements are 7.99 s and 504.7 s respectively.By using multi-scale iteration technique, the computation burden is decreased greatly and the spending time reduces to 1.58% of traditional BP imaging algorithm.

Example 2
As shown in Figure 3a, three metal bars with infinite length are located at points (0.4, 0.5), (1.0, 0.7) and (1.6, 0.7) respectively.The scanning parameters are the same as that in Example 1. White Gauss noise are added to GPR original signal of each A-scan with SNR = 0 dB.The radargram after pre-processing is shown in Figure 3b.The 1st iteration imaging results and PTR delineation results are shown in Figures 3c,d.After 2nd iteration, the resolution reaches the specified requirement.The 2nd imaging result is shown in Figure 3e, while imaging result obtained by using traditional BP imaging algorithm is shown in Figure 3f.The iteration parameters in the imaging processing is shown in Table 2.The focusing parameters of the above two imaging results corresponding to the same PTR shown in Figure 3d are 0.0766 and 0.0071 respectively.The focusing quality of the imaging result improved above 970%.Both algorithms are implemented with the same PC as that in above example and time consuming are 20.9 s and 248.5 s respectively.

Example 3
A GPR survey was conducted along a provincial road in Hunan Province for road quality assessment.An on-site contained several cement plates newly overhauled was chosen.Several reinforced rods are embedded in these cement plates.We use a impulse GPR system equipped with a 2 GHz antenna pair for subsurface targets detection.The B-scan line is selected along the road with a length of 3.45 m.Each A-scan contains 512 sampling points in the total 26 ns time window.Original GPR echo profile and B-scan after preprocessing are shown in Figure 4a,b respectively.The 1st and 2nd iteration imaging results and PTR delineation results are shown in Figure 4c-f.After 2nd iteration,Δx 3 < Δx, the iteration processing stopped.The imaging results of multi-scale weighted BP and traditional BP imaging algorithm are shown in Figure 4g,h, respectively.The iteration parameters in the imaging processing is shown in Table 3.

Conclusions
A multi-scale weighted back projection imaging algorithm is presented in this paper.By analyzing the imaging mechanism of traditional back projection imaging algorithm, a multi-scale back projection imaging algorithm is presented base on a hierarchal discretization technique and intensity detection technique.By reducing imaging resolution of regions where there haven't potential targets, the multi-scale back projection imaging algorithm can reduce computation burden greatly.By introducing a weighted factor, the imaging resolution in potential target region is also improved.Simulation results and on-site experiments demonstrated its validity in ground penetrating radar fast and high resolution imaging occasions.The quantitative analysis of the performance of presented imaging technique depends on specific ground penetrating radar scanning scenarios and radar-gram.In example 1, the focusing quality of the imaging result improved above 680% and the spending time reduces to 1.58% compared with traditional back projection imaging algorithm.In example 2, the focusing quality improved above 970% and the spending time reduces to 0.36%.In real data experiment, the focusing quality improved above 3000% and the spending time reduces to 16% respectively.From the above examples, we can see that the imaging accuracy improved greatly with much lower computation burden.Key parameters such as initial sampling grid and hierarchical refinement coefficient have much influence in imaging accuracy and iteration speed of imaging procedure.Future work will focus on the selection criteria of initial sampling grid and hierarchical refinement coefficient.