A Review on Migration Methods in B-Scan Ground Penetrating Radar Imaging

Abstract

Even though ground penetrating radar has been well studied and applied by many researchers for the last couple of decades, the focusing problem in the measured GPR images is still a challenging task. Although there are many methods offered by different scientists, there is not any complete migration/focusing method that works perfectly for all scenarios. This paper reviews the popular migration methods of the B-scan GPR imaging that have been widely accepted and applied by various researchers. The brief formulation and the algorithm steps for the hyperbolic summation, the Kirchhoff migration, the back-projection focusing, the phase-shift migration, and the ω -k migration are presented. The main aim of the paper is to evaluate and compare the migration algorithms over different focusing methods such that the reader can decide which algorithm to use for a particular application of GPR. Both the simulated and the measured examples that are used for the performance comparison of the presented algorithms are provided. Other emerging migration methods are also pointed out.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Academic Editor:Fatih Yaman

Department of Electrical and Electronics Engineering, Mersin University, Ciftlikkoy, 33343 Mersin, Turkey

Received 14 March 2014; Accepted 5 May 2014; 10 June 2014

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Ground penetrating radar (GPR) is remote sensing technique that can nondestructively sense and detect objects inside the visually opaque environment or underneath ground surface. GPR has gained its fame thanks to its high resolution capability and applicability in numerous fields such as detecting mines and unexploded ordinances, finding water leakages, investigating archeological substances, spotting asphalt/concrete cracks in highways, searching buried victims after an earthquake or an avalanche, and imaging behind the wall for security applications [1-4]. Depending on the application, different scanning schemes, namely, A-scan, B-scan, and C-scan, are being employed [1]. In the B-scan measurement situation, a downward looking GPR antenna is moved along a straight path on the top of the surface while the GPR sensor is collecting and recording the scattered field at different spatial positions. This static measured data collected at single point is called an A-scan [1].

The data collection in a typical GPR operation can be either employed in the time domain by recording the scattered response of a time-domain pulse or in the frequency domain by recording the frequency response of the scattered field. For the former case, the two-dimensional (2D) space-time GPR image I(x,t) is attained by conveying a time-domain pulse towards the surface at consecutive, distinct synthetic aperture points. For the latter case, an inverse Fourier transform (IFT) operation should be accommodated to carry the collected data from the spatial-frequency domain to space-time GPR image. In any case, the depth resolution is achieved by the transmitted signal's frequency diversity. The resolution along the scanning direction is attained by the synthetic processing of the received data collected at different spatial points of the B-scan. While a fine resolution in the depth axis is usually easy to get by utilizing a wide-band transmitted signal, the resolution along the scanning direction is much harder to realize and requires special treatment. On the other hand, some specific GPR applications like preservation tasks in ancient constructions or stone masonry structures and detection of antipersonnel landmines for humanitarian demining require high resolution images of the investigated structures or regions especially in the cross-range dimensions.

In a typical space-time B-scan GPR image, any scatterer within the image region shows up as a hyperbola because of the different trip times of the EM wave while the antenna is moving along the scanning direction. For this reason, the resolution along the synthetic aperture direction depicts undesired low resolution features owing to the long tails of the hyperbola. Therefore, one of the most applied problem for the B-scan GPR image is to transform (or to migrate) the unfocused space-time GPR image to a focused one showing the object's true location and size with corresponding EM reflectivity. The common name for this task is called migration or focusing [5-18]. In fact, migration methods were primarily developed for processing seismic images [19] and they were also applied to the GPR thanks to the likenesses between the acoustic and the electromagnetic wave equations [7-11]. Kirchhoff's wave-equation [5] and frequency-wave number (ω -k ) based [6, 7] migration algorithms are widely recognized and employed. The wave number domain focusing techniques, for instance, were first formulated for seismic imaging applications [7] and then adapted to the contemporary synthetic aperture radar (SAR) imaging practices [12-15]. These algorithms are also named as seismic migration [12, 13] and frequency-wave number (or ω -k ) migration [16-18] by different researchers.

In this work, we present a brief review of the B-scan GPR migration/focusing methods that are commonly used by GPR research community. Although we have done some preliminary studies on the performance comparison of only two GPR focusing algorithms earlier [20] by utilizing very limited performance parameters, in this paper, we extend our studies to compare a total of five popular algorithms by assessing different aspects of these algorithms. In Section 2, we address the problem of migration together with the well-known exploding source concept . In fact, most migration algorithms make use of this concept as we shall explain in this work. In Section 3, most popular and functional migration methods for the B-scan GPR imaging applications are reviewed. For this purpose, (i) hyperbolic summation, (ii) Kirchhoff's migration, (iii) phase-shift migration, (iv) frequency-wavenumber (or the ω-k ) migration, and (v) back-projection method are presented together with the detailed algorithm steps. In Section 4, the success and the performance of these algorithms were compared to each other by testing them over the simulated and the measured B-scan GPR data. Our main objective in this work is to present a qualitative comparison of these popular and widely used algorithms such that the reader can benefit from the performance results over different parameters ranging from resolution to computation time. Section 5 is dedicated to discussions and the conclusion.

2. The Problem of Migration/Focusing

The common goal in a typical GPR image is to display the information of the spatial location and the reflectivity of an underground object. As illustrated in Figure 1, a single point scatterer appears as a hyperbola in the space-time GPR image for the monostatic operation. Since this is the expected image pattern for the B-scan operation, such information can be thought as sufficient if the main goal of the GPR application is just to sense a pipe or comparable objects. However, the information about the scatterer including its depth, its size, and its EM reflectivity can be very important in most GPR applications. Therefore, the hyperbola or dispersion in the space-time B-scan GPR image ought to be converted to a focused one that demonstrates the object's factual location and size together with its scattering amplitude. A focused or migrated image is gathered after the elimination of this hyperbolic type of diffraction or any other kind of dispersion [5-18].

Figure 1: (a) Typical GPR measurement setup (B-scan) and (b) resultant space-time image that contains the well-known hyperbolic aliasing.

[figure omitted; refer to PDF]

2.1. Exploding Source Model

Most of the migration methods are based on the concept called exploding source model (ESM). In 1985, Claerbout [21] came with the clever idea of thinking of the scattered field at the radar receiver as if it is originated from the source at the target location. Instead of assuming a two-way trip convention of the EM wave, therefore, it is imagined that a fictitious source "explode" at a reference time of t=0 around the target location and send EM wave to the receiver as illustrated in Figure 2. The real data collection scheme is shown in Figure 2(a) where, in fact, the two-way propagation between the radar and the object exists. When the ESM is utilized, however, the collected data is assumed to be originally radiated from the source on the object. Therefore, one-way propagation is assumed in the ESM depicted in Figure 2(b). Since the trip time of the EM wave would be half of the original problem the compensation should be made for the velocity of the EM wave by just dividing it by two to have _vm =v/2 where v is the speed of the EM wave within the medium of propagation.

Geometry for (a) B-scan GPR data collection scheme and (b) utilizing "exploding source model."

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

The migration using the ESM is essentially carried out by applying these two actions:

(i) the received signal is extrapolated back to exploding source points;

(ii) the migrated image is realized by forming the back extrapolated EM wave at the time of t=0 .

3. Migration/Focusing Methods

In this section, we will briefly review the basic steps of the mostly applied GPR algorithms, namely, the hyperbolic summation, the Kirchhoff migration, the phase-shift migration, the ω-k (Stolt) migration, and the back-projection focusing. Other migration methods will also be mentioned at the end of this section.

3.1. Hyperbolic (Diffraction) Summation

In a typical GPR application, the radar antenna collects the scattered or back-scattered EM wave from the air-to-ground interface and subsurface objects together with many cluttering effects mainly attributed to inhomogeneities within the ground. For the idealistic case, the phase of the scattered signal is directly proportional to the trip time (or distance) that the EM wave possesses if the propagation medium is homogenous. The monostatic backscattered signal from a single point-like scatterer experiences different round-trip distances while the antenna is moving over the surface for the B-scan operation. For each static measurement along the B-scan axis, one-dimensional (1D) range profile (or depth profile) of the subsurface scene is obtained by taking the inverse Fourier transform (IFT) of the frequency-diverse back-scattered signal. After putting all the depth profiles aside, the 2D space-time (or space-depth) B-scan GPR image is obtained. As the GPR antenna has a finite beamwidth, any subsurface object is illuminated for a finite length along the B-scan axis as demonstrated in Figure 1. Therefore, this object shows up as a parabolic hyperbola in the space-time GPR image due to different trip distances of the EM wave between the radar and the illuminated scattering object. The true location of the object is, in fact, at the apex of this hyperbola.

Let us assume a perfect point scatterer situated at (_xo ,_zo ) in the 2D plane where x -axis corresponds to the scanning direction and the z -axis is the depth as illustrated in Figure 1. If the propagating medium is homogeneous, the parabolic hyperbola in the GPR image can be characterized by the following equation when the radar is moving on a straight path along X -axis. Consider [figure omitted; refer to PDF] In the above equation, X is the synthetic aperture vector along the B-scan and R represents the path length vector from the antenna to the scatterer. Assuming that the resultant B-scan GPR image can be regarded as the contribution of finite number of hyperbolas that correspond to different points on the object(s) below the surface, the following methodology can be applied to migrate the defocused image structures, that is, hyperbolas to the focused versions [22].

(i) For each pixel point, (_xi ,_zi ) on the 2D original B-scan space-depth GPR image matrix, find the associated hyperbolic template using (1) and trace all the pixels under this template.

(ii) Record the image data for the traced pixels under that template. At this stage of the algorithm, we have 1D field data ^Es , whose length N is the same as the number of sampling points along the synthetic aperture X .

(iii): Then, calculate the root-mean-square (rms) of the entire energy contained inside this 1D complex field data as follows: [figure omitted; refer to PDF]

(iv) The calculated rms figure is recorded on the new image matrix at the point (_xi ,_zi ) . This process is reiterated until all pixels on the original GPR image are passed through the algorithm.

3.2. Kirchhoff's Migration

Kirchhoff's migration (KM), also known as reverse-time wave equation migration , is mathematically equivalent to hyperbolic summation method with some correcting factors included in the solution [22]. In Kirchhoff's migration procedure, the aim is to find the solution to the following scalar wave equation for the wave function [straight phi](x,z,t) within the propagating medium: [figure omitted; refer to PDF] Here, _vm is the velocity of the EM wave within the propagating medium and taken as v/2 by utilizing the "exploding source" concept. The solution to this differential equation is available for the far-field approximation from the Kirchhoff integral theorem [23] as [figure omitted; refer to PDF] Here, θ is the angle of the incident wave to the depth axis (z) and r=^{((x-_xm )2 +z2 )1/2} is the path from the target point at (x,z) to the observation point at (_xm ,0) .

In KM method, the following correction factors that are not considered in the hyperbolic summation are accounted.

(i) Compensation for the spherical spreading is taken into account. For this purpose, a correction factor of 1/_vm r is used for the 2D propagation of the EM wave, whereas this factor becomes 1/_(vm r) for the 3D propagation.

(ii) The directivity factor cos......θ is also considered. This factor corrects the diffraction amplitudes.

(iii): The phases and amplitudes of the wave are also corrected. The phase is corrected by π/2 and π/4 for 2D and 3D propagation cases, respectively. The amplitude of the EM wave is corrected with a factor proportional to the square of the frequency for the 2D propagation case and to the frequency for the 3D propagation case [23-25].

3.3. Phase-Shift Migration

The phase-shift migration (PSM) method is first introduced and applied by Gazdag [6]. Similar to Kirchhoff's migration, this method also utilizes the ESM concept [26]. In brief, the algorithm iteratively puts a phase-shift to migrate the wave field to the exploding time of t=0 such that all the scattered waves are drawn back to object site to have a focused image.

The main aim of the PSM algorithm is to calculate the wave field at t=0 by extrapolating the downward (z -directed) EM wave with the phase factor of exp......(j_kz z) .

The PSM algorithm can be briefly summarized via the following steps.

(i) First, 2D measured raw dataset in frequency wave-number or (ω,_kx ) domain is multiplied by a phase-shift factor, K , along the depth axis (or z -axis) as given below [figure omitted; refer to PDF] This factor can also be written in terms of wave-number along the data collecting axis _(kx ) as [figure omitted; refer to PDF] Putting k=ω/_vm (where ω is the angular frequency and _vm is the speed of wave inside the propagation medium) into (5), one can get [figure omitted; refer to PDF] In the above equation, the incremental depth parameter Δz is proportional to the time sampling interval Δt of the input data via Δz=_vm ·Δt . After this modification, the final form of the phase-shift factor becomes [figure omitted; refer to PDF] 2D measured raw dataset ^Es (ω,_kx ) is multiplied by K along the depth axis for the time steps of Δt .

(ii) By utilizing the ESM concept, the imaging task is accomplished by taking the inverse Fourier transform (IFT) of the ^Es (ω,_kx ) after selecting the time variable as Δt=0 . Therefore, only single FT operation is required at one point (when Δt=0 ) for the focused image.

After updating the factor K for every value of Δt and ω , we have the data in 3D (ω,_kx ,_kz ) domain as [figure omitted; refer to PDF] The new dataset ^{Es [variant prime]} (ω,_{kx ,kz} ) is summed up along the frequency axis and indexed for different values of Δt as [figure omitted; refer to PDF]

(iii) Setting Δt=0 and taking the 2D IFT with respect to _kx and _kz , the focused the image in the (x,z) domain can be obtained as given below: [figure omitted; refer to PDF]

3.4. Frequency-Wavenumber (Stolt) Migration

Frequency-wavenumber (ω -k ) migration, also known as Stolt migration or the f-k migration, utilizes the ESM idea and the scalar wave equation [7]. The algorithm behind the frequency-wavenumber method works faster than the previously presented migration methods. The ω -k migration method has proven to be working well for the constant-velocity propagation mediums [6, 7]. The solution of the ω -k migration can be rewritten to be the same as the solution of the Kirchhoff migration [18]. Below is the brief explanation of the algorithm.

The algorithm begins with the 3D scalar wave equation for the wave function [straight phi](x,y,z,t) within the constant-velocity propagation medium [figure omitted; refer to PDF] In the Fourier space, spatial wave-numbers and the frequency of operation are related with the following equation: [figure omitted; refer to PDF] Stratton [27] demonstrated that any given wave function can be written as the summation of infinite number of plane wave functions, say E(_kx ,_ky ,ω) , as [figure omitted; refer to PDF] For the GPR operation, the scattered field is taken to be measured on the z=0 plane above the surface. When carefully treated, the above equation offers a Fourier transform pair where e(x,y,t) can be regarded as the time-domain measured field on the z=0 plane. Consider [figure omitted; refer to PDF] It is very important to notice that this equation designates a 3D forward FT relationship between e(x,y,t) and E(_kx ,_ky ,ω) for the negative values of time variable, t . Then, the inverse FT can be dually defined in the following way: [figure omitted; refer to PDF] Afterwards, we now use the ESM to focus the image by setting t=0 in (15) and using [figure omitted; refer to PDF] Therefore, one can get the time-domain measured field via [figure omitted; refer to PDF] The above equation presents a focused image. However, the data in (_kx ,_ky ,ω) domain should be transformed to (_kx ,_ky ,_kz ) domain to be able to use the FFT. Therefore, a mapping procedure from ω domain to _kz domain is required for fast processing. The relationship between the ω - and -_kz and dω - and -_dkz can be easily obtained from (13) as [figure omitted; refer to PDF] Substituting these equations into (18), one can obtain the following: [figure omitted; refer to PDF] Here, ^Em (_kx ,_ky ,_kz ) is the mapped version of the original data E(_kx ,_ky ,ω) . After this mapping, the new data set does not lie on the uniform grid due to nonlinear feature of the transformation. Therefore, an interpolation procedure should also be applied to be able to use the FFT for fast processing of the collected dataset. Equation (20) suggests a well-focused image of the subsurface region that may contain a finite numbers of scatterers. Since the FFT routine is utilized, the GPR image in the ω -k migration is obtained quite fast. It is also important to express that the mapped dataset is scaled by the factor of _{^vm2 kz} /ω . This scaling is sometimes called the "Jacobian transformation from ω to _kz . "

The above focusing equation is valid for the 3D GPR geometry or the C-scan problem. In most subsurface imaging problems, however, the raw data is collected in 2D space, that is, space-time (space-depth) or space-frequency. Therefore, the focusing equation in (20) can be easily reduced to 2D B-scan GPR problem in the space-depth domain via the following equation: [figure omitted; refer to PDF]

3.5. Back-Projection Based Migration

The last category of the GPR focusing algorithms is based on the tomographic principles used in medical imaging and collectively termed as the back-projection (BP) algorithms. Since its first formulation for the 2D monostatic SAR processing [28], the BP algorithm has gathered an increasing interest in radar community thanks to its distinct features that are proved to be very useful for various SAR imaging applications. For example, the algorithm does not require a straight and uniformly sampled synthetic scan aperture due to its serial processing nature. More clearly, each 1D range profiles are serially processed and spread or back-projected over the entire 2D image independently. This sequential processing means "real-time" operation capability and hence the imaging process can begin before acquiring the entire synthetic aperture data. Furthermore, the specific subsections of the region to be imaged can be easily selected to investigate these subsections more closely. In applications where the approximate location of the target is known a priori, the detailed image of the region around this target can be easily formed by the algorithm.

Before briefly explaining the algorithm, let us first consider the B-scan, monostatic GPR geometry shown in Figure 3. A radar antenna, situated at a height of h from the ground, transmits a frequency-diverse waveform at each distinct point m along the synthetic axis. The relative permittivity of the subsurface medium is denoted by _{[straight epsilon]r} and the reflectivity function of the scatterers is represented by e(x,z) . The instantaneous position of the antenna (_xm ,_zm ) is defined by a unit vector _um pointing from the scene center towards this location. The corresponding view-angle _θm is defined as the angle between the unit vector _um and the depth axis z . Assuming a stepped-frequency continuous waveform (SFCW) transmission, the back-scattered signal at a specific observation angle can be written as [figure omitted; refer to PDF] where _kr is the wavenumber defined for the two-way propagation as _kr =4πf/v , f is the frequency, v is the speed of the propagation, and _rm is the range from the instantaneous antenna location (_xm ,_zm ) to the points (_x0 ,_z0 ) within the imaging scene. The range profile of the illuminated scene can be obtained by employing 1D IFT to (22) and can be mathematically conveyed as [figure omitted; refer to PDF] which is nothing but the Radon transform of the scene.

Figure 3: The B-scan geometry of the 2D monostatic GPR application.

[figure omitted; refer to PDF]

The derivation of the BPA [29] begins with the implementation of IFT of the scattering function e(x,z) written in Cartesian coordinates as [figure omitted; refer to PDF] where E(_kx ,_kz ) is the 2D FT of e(x,z) . Equation (24) can be reformed to be rewritten in the polar coordinates (_kr ,_θm ) as follows: [figure omitted; refer to PDF] Now, the projection-slice theorem [30] can be used to relate the target's FT E(_kx ,_kz ) to the collected measured data _Sθ (_kr ) . For the 2D problem, the theorem essentially states that 1D FT of the projection at the angle θ represents the slice of the 2D FT of the projected (original) scene at the same angle; that is, _Sθ (_kr )=E(_kr ,θ) . Therefore, the sampled representation of E(_kx ,_kz ) can be obtained from the FT of the projections _Sθ (_kr ) measured at several observation aspects. By the help of this principle, (25) becomes [figure omitted; refer to PDF] The bracketed integral term in (26) can be regarded as the 1D IFT of a function _Qθm (_kr )=_Sθm (_kr )_kr calculated at _rm . Defining _qθm (r) as the IFT of this function, (26) can be represented as [figure omitted; refer to PDF] Equation (27) is the final focused image of the 2D filtered back-projection algorithm. For the SFCW system, the execution of the algorithm can be summarized as follows.

(1) Preallocate an image matrix of zeros e(x,z) to hold the values of the scene reflectivity.

(2) Multiply the acquired spatial-frequency data _Sθm (_kr ) with _kr .

(3) Take 1D IFT of the result to obtain _qθm (r) which represents the filtered version of the range profile _sθm (r) .

(4) For each pixel position in the image data, evaluate the corresponding range value _rm and acquire its _qθm (_rm ) value by the help of an appropriate interpolation method.

(5) Iteratively add interpolated data values to e(x,z) .

(6) Repeat the above steps from 2 to 5 to cover all the observation angles _θm .

3.6. Some Other Methods

In addition to the above mentioned common migration methods, there are many other techniques that have been introduced and studied by different GPR researchers. Fisher et al. [31] applied reverse-time migration procedure for GPR profiles. Capineri et al. [32] applied a technique based on Hough transformation to the B-scan GPR data to attain better resolved images of pipe structures. Leuschen and Plumb [9] and Morrow and van Genderen [33] realized back-propagation procedures that rely on finite difference time-domain (FDTD) reverse-time focusing methods to resemble the focusing matter in GPR images.

4. Performance Comparison of the Algorithms

Performance of the algorithms is evaluated by the help of following various parameters that are resolution, integrated sidelobe ratio, signal to clutter ratio, and computation speed.

Resolution . GPR measurement can provide fine resolutions both in depth and azimuth resolutions. The term depth (range) corresponds to the line of sight distance from the radar to the target to be imaged. The term azimuth (transverse range, cross range) is used for the dimension that is perpendicular to range or parallel to the radar's along-track axis [34]. The depth resolution is defined as an ability of the radar equipment to separate two reflecting objects on a same line of sight, but at different ranges from the antenna. And the azimuth resolution can be defined as a capability of the separation of two reflecting objects on the same ranges, but different aspects from the radar. In GPR operation, the high resolution in depth is obtained by utilizing a transmitted signal of wideband. Fine azimuth resolution is achieved by coherently processing the target's electromagnetic scattering measured at different aspects while the radar is moving along a straight path. However, the measured resolutions after the postprocessing are dependent on the focusing ability of the migration algorithm. That is why the resolution performance of the algorithms is tested by measuring the -4 dB contours of a point scatterer (whose coordinate is x=0 m and z=1.5 m) in the resultant images of simulations. Normally, algorithms with better focusing abilities are expected to produce lower resolved images.

Integrated Sidelobe Ratio (ISLR) . Since ISLR is defined as the ratio of the total energy within the sidelobes to the peak energy of the main lobe [35], it can be regarded as an important parameter when assessing the focusing abilities of migration algorithms [36]. Based on Sanchez's definition of ILSR, it is the ratio of the -3 dB width of the main lobe to the rest of the energy within all lobes. Therefore, the ISLR of a 2D image function I(x,y) can be formulated as follows: [figure omitted; refer to PDF] where ^{∫-3 dB} Id^{x[variant prime]} d^{y[variant prime]} is the energy within the -3 dB width of the main lobe and ^∫-∞+∞ Id^{x[variant prime]} d^{y[variant prime]} is the total energy.

Signal to Clutter Ratio (SCR) . Signal to clutter ratio (SCR) is also a crucial parameter when evaluating the quality of the focusing ability. In a well-focused image, the clutter/noise level should be much lower than the signal level to have a well-contrasted image. The SCR can be defined as the ratio of the received target signal power to the received clutter power as shown below: [figure omitted; refer to PDF] where _Ptar is the received target signal power and _Ptot is the total power within the image.

Computation Speed . The processing time of the computation of the algorithms can be crucial if the GPR application requires processing of vast data such as scanning and cleaning a minefield.

4.1. Comparison over Simulation

The migration algorithms were first tested through an imaging simulation performed in MATLAB. Considering the classical B-scan geometry shown in Figure 4, the subsurface was assumed to have a completely homogeneous soil medium structure (i.e., constant-velocity medium) with a dielectric constant of 2.4. The medium is also assumed to be nonmagnetic; that is, _μr =1 . The monostatic operation was considered and the -3 dB beamwidth of the antenna was assumed to be large enough such that all targets to be imaged fall inside the beam for the whole aperture with constant antenna gain. The electric field data of the isotropic point scatterers with locations depicted in Figure 4 were then collected along a straight path ranging from x=-2.5 m to x=2.5 m and with 2 cm steps. At each spatial point, the frequency response of the subsurface was acquired for the 1.25 GHz-3.75 GHz frequency range that was uniformly sampled at 168 points. Therefore, 251 × 168 spatial-frequency B-scan data E(x,f) were generated for the investigated scene. The frequency data were subsequently preprocessed by applying a Hanning window for sidelobe control and 4 times zero padding for interpolation purposes.

Figure 4: Simulation geometry of the B-scan imaging example.

[figure omitted; refer to PDF]

Figure 5(a) shows the raw, unfocused B-scan image simply formed by applying a 1D inverse Fourier transform to the collected data, E(x,f) . As expected, the target reflections are observed as hyperbolas with different intensities and curvature slopes. To collapse these hyperbolic signatures at the apex of the hyperbolas, the migration algorithms were applied and the corresponding obtained results were given in Figures 5(b)-5(f). As a quick interpretation, it can be noticed from all images that the scattering mechanisms are almost concentrated around the exact locations of the targets.

Simulation results for the -40 dB dynamic range.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

(f) [figure omitted; refer to PDF]

To compare the performances of the algorithms, Table 1 lists the success merits of the algorithms including the depth resolution, the azimuth resolution, the ISLR, and the computational time. By looking at the resolution outcomes in both depth and azimuth direction, KM and ω -K algorithms seem to provide better resolved images. On the other hand, KMA is the worst in terms of computation time while ω-k A is the fastest among all others. In terms of the ISLR, BPA is superior to any other algorithm by providing a sharp ISLR that is at least 6 dB better than the others. When comparing all parameters, one can conclude that the BPA and ω-k A seem to have a little bit better in terms of the migration ability.

Table 1: Performances of focusing algorithms over different parameters for the simulated experiments.

Algorithm	Computation time (s)	Depth resolution (cm)	Azimuth resolution (cm)	ISLR (dB)
HSA	46.32	6.78	12	-13.15
KMA	9,372.70	5.43	4	-12.17
PSA	2.13	8.70	10	-13.02
ω -k A	0.41	6.30	4	-12.75
BPA	4.29	7.70	5.8	-6.13

When the resultant focused images of these algorithms are visually compared in Figure 5, the KM, the BPA, and the ω -k A look more successful than others as the results of the HSA and PSA expose some image degradations. If those well-localized images of Figures 5(c), 5(e), and 5(f) are compared, some minor differences can even be discerned between each other. For example, the target reflections in the BP migrated image have more regular pattern and also have stronger intensities. Also, the images for the KM and the ω -k A have some relatively higher sidelobe levels around the top and bottom targets. Additionally, if the image for the ω -k A is carefully checked (see Figure 5(e)), the depth of the bottom target is also shown to be mapped to z=4 m which indicates a little deviation from its true value of z=4.5 m.

By evaluating the GPR images in Figure 5, we can compare the focusing algorithms against each other as follows. While the ω -k A may experience some problems in imaging deep targets with low reflectivities, successful imaging of shallow targets can be problematic for the KM. Secondly, the images for the HSA and the PSA look considerably poorer than those of the other algorithms. This can be seen from the significantly defocused image that corresponds to PSA (see Figure 5(d)). This image also has some extra undesired scattering features such as blurring and spreading. On the other hand, the image for the HSA exhibits a major degradation only around the uppermost target. The other targets are observed to be almost well focused.

To comprehend the sources of abovementioned defocusing effects for these two algorithms, another simulation was performed by utilizing different simulation parameters. For this simulation, the total number of frequency samples was increased to 302 points. The bandwidth was then correspondingly changed to B=4.5 GHz to satisfy the same unambiguous range of the previous simulation, that is, 6.5 m. Corresponding results for the simulation are demonstrated in Figure 6. It is shown from the Figure 6(a) that the HSA is not affected from the variation of the sampling points along the depth direction. The migrated image is similar to the previous result given in Figure 5(b), except for the definite resolution enhancement thanks to wider bandwidth. In any case, HSA is shown to have difficulty only in providing a reasonable representation of the uppermost target along the spatial direction. Although this aliasing is also somewhat observed for the deeper targets, they can also be regarded as well focused within the dynamic range of -40 dB . As stated before, this spatial aliasing for the shallower target is also slightly observed in the result of the KA as given in Figure 5(c). Thus, noting the similarities between these algorithms, this fact can be assumed to stem from the characteristics of these methods. Nevertheless, this spatial aliasing can be avoided by utilizing different implementations, such as the ones that incorporate more points along a hyperbolic trajectory, rather than using a single point during integration.

Results of the hyperbolic summation and Kirchhoff's migration algorithms for a decreased frequency sampling interval (i.e., frequency bandwidth B=4.5 GHz with 302 sampling points).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

The PSA result seen in Figure 6(b) conversely shows drastic improvement when compared to the previous one shown in Figure 5(d). All of the four scatterers are shown to be well localized for the increased value of the frequency sampling points. This notable change can be attributed to the algorithms' iterative processing nature along the depth direction. Since the PSA first records the data at the surface and repeatedly multiply this data with a phase factor for each downward range points, the short sampling intervals would result in better quality images. Hence, this synthetic example illustrates that the PSA is highly sensitive to the employed frequency sampling interval. Lastly, due to the algorithmic similarities between the PSA and the ω -k A, the corresponding images given in Figures 6(b) and 5(e) are shown to have similar distortion patterns especially seen at the bottom of these images.

Finally, the computational efficiency of the algorithms is tested by looking at the simulation runtime and the required memory during implementation. Based on our evaluation, the computational efficiency of the algorithms is listed from the best to the worst as follows: the ω -k A, the PSA, the BPA, the HSA, and the KM.

4.2. Comparison over Measurement

The performance of the focusing algorithms is also tested by the help of measured data sets. For this purpose, two different experiments were conducted to better understand the differences between the focusing capabilities of the algorithms for real datasets. These experiments were conducted by the help of the Agilent E5071B ENA Vector Network Analyzer (VNA) that can generate SFCW signals within the frequency range of 0.8-8.5 GHz. In both experiments, a C-band double-ridged pyramidal rectangular horn antenna was used as a monostatic transceiver.

4.2.1. Experiment Number 1

In the experiment, we have constructed a test bed by building a big wooden pool with a size of 190 cm × 100 cm × 80 cm and filled it with homogeneous and dry sand material. The dielectric constant of the sand was measured to be almost constant around 2.4 for the frequency range between 1 GHz and 8 GHz. It is also assumed that the sand material has the unit magnetic permeability for the frequencies of operation. Two cylindrical metal rods with different sizes were selected as targets in the scene. The thin rod with 4.5 cm in diameter and 45 cm in length was buried at (x=-12 cm, z=65 cm) and a thick rod with 6 cm in diameter and 32 cm in length was buried at (x=18 cm, z=75 cm). A B-scan measurement was accomplished by moving the antenna along the perpendicular direction of the target axes and by spanning a synthetic aperture length of L=1.34 m sampled at 63 spatial points. The stepped frequency response of the medium was then obtained for a bandwidth of 0.8-4.5 GHz that was uniformly sampled at 751 points.

After collecting the data, the classical space-depth B-scan GPR image is obtained as shown in Figure 7(a) from which the unfocused target signatures can be clearly seen. The governing scattering mechanism from air-ground boundary is also easily detected around at z=40 cm and seen throughout the entire synthetic aperture. Then, we apply the above presented migration algorithms in aiming at getting a better focused image. The resultant migrated images are displayed in Figures 7(b) to 7(f) wherein the target responses are seen to be more localized around their correct locations. To better compare the algorithms against each other, all resultant focused images were displayed within the same dynamic range of -60 dB. The focusing performance of all algorithms seems to be very similar to each other as the tails (or the sidelobes) of the targets' images are comparable dispersion features. In terms of the image contrast, the ω -k A and the PSA outperforms the BPA, the HSA, and the KM as it can be clearly deducted from the images. For this experiment, the numerical noise generated by the iterative implementation of the BPA seems to be very high when compared to other algorithms. In fact, the real data contains infinite number of scatterers (with high or low reflectivities) under the ground; it seems that the BPA is more sensitive to reflectivity amplitude than the others.

Results for the real soil Experiment 1. Note the different dynamic range (i.e., -60 dB) for the Kirchhoff migration result.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

(f) [figure omitted; refer to PDF]

4.2.2. Experiment Number 2

In the second experiment, the algorithms were tested under a more heterogeneous soil in another test bed of the indoor environment. The targets were selected as a water-void target, a metal pipe, and a metal plate whose sizes and burial locations were as depicted in Figure 8. The synthetic aperture length was set to 134 cm for a total of 68 discrete points and the frequency was changed from 0.8 GHz to 5 GHz for 501 discrete points. A picture during the experiment is presented in Figure 9.

Figure 8: Geometric representation of the targets used in Experiment 2.

[figure omitted; refer to PDF]

Figure 9: A scene from real soil measurements.

[figure omitted; refer to PDF]

After collecting the data, the raw image obtained by taking the IFT of the measured back-scattering data is shown in Figure 10(a). As expected, the image is unfocused around the buried objects. After applying the migration algorithms, the focused images of the three buried objects are acquired as seen in Figures 10(b) to 10(f). For this experiment, very similar features were observed as in the case of the first experiment. Again, focusing performances of the algorithms are alike as the sidelobe contours of target's images are comparable. The contrast performances of the images also yield similar finding as we found in the first experiment. Again, the ω -k A and the PSA produce more quality images when compared to other three algorithms. As the final comments, the BPA seems to have problems again when dealing with real-world data as it produces significant numerical noise as it can be seen from Figure 10(f).

Results for the real soil Experiment 2. Note the different dynamic ranges (i.e., -50 and -60 dB) for the hyperbolic summation and Kirchhoff's migration results.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

(f) [figure omitted; refer to PDF]

Table 2 summarizes the performance of the presented five different algorithms in terms of processing time and the SCR over the conducted measurements. In assessing the speed of the algorithms, we observe similar performances as in the case of the simulation results that are listed in Table 1. Again, ω-k A is the fastest one and the BPA is the second best while the KMA is the worst. If we look at the SCR results, the clutter suppression ability of BPA is the best while the others demonstrate moderate performances when compared to BPA.

Table 2: Performances of focusing algorithms over different parameters for the measured experiments.

	Experiment number 1 results	Experiment number 2 results
	SCR (dB)	Time (s)	SCR (dB)	Time (s)
HSA	-6.27	11.12	-12.48	7.61
KMA	-4.55	5,723.80	-9.47	26,76.30
PSA	-5.41	13.99	-10.40	3.74
ω -k A	-5.66	0.57	-9.92	0.34
BPA	-4.48	1.75	-9.01	2.19

5. Conclusion and Discussions

In this paper, we have reviewed and compared fundamental algorithms that are used to focus the B-scan GPR data by presenting the algorithm steps. Namely, the hyperbolic summation, the Kirchhoff migration, the phase-shift migration, the frequency-wavenumber (ω -k ) migration, and the back-projection based migration algorithms have been explained and applied to both the simulation and the measurement experiments. Based on the results obtained from the simulated and measured images, we have the following remarks for the investigated migration algorithms.

(i) The HSA is conceptually very simple and therefore it is very easy to implement. On the other hand, numerous hyperbolic template vectors should be applied for every pixel in the image which, in turn, is an expensive process in terms of the computation time. If the image size is relatively big, then the whole process of calculating the energy under the hyperbolic templates of every pixel may take a very long time. Although the focusing ability of the algorithm is not as good as the ω -k A and the BPA, the measurement results have showed that focusing success of the HSA also has some merit as it can be seen from the constructed GPR images. Also, the HSA produce fair contrast images when compared to others.

(ii) Although KM is conceptually similar to HSA, the algorithm is based on the scalar wave equation and it corrects the amplitude and the phases of some parameters such as spherical spreading and directivity that are not taken in to account in HSA. Therefore, the results obtained by KM are a little bit better than the ones with HSA that can also be seen from the resulted images. This is, of course, at the price of a longer simulation time. The improvement obtained by this amplitude correction can be clearly observed in Figures 7(c) and 10(c) to be compared to the images in Figures 7(b) and 10(b). The targets with smaller reflectivity are more visible in the images obtained with KM. By looking at both the simulation and the measurement results, KMA seems to provide the best resolution figures when compared to other four algorithms.

(iii): The PSA utilizes the ESM concept and iteratively tries to add phase-shifts to migrate the wave field to the exploding time of t=0 . The algorithm uses the advantage of FFT; therefore, it is much faster than the abovementioned algorithms. The PSA's performance in terms of computation time, focusing, and image quality is modest when compared to others.

(iv) Stolt migration algorithm (or the ω -k A) also makes use of the ESM concept and the scalar wave equation for the scattered field. It can be shown mathematically that the ω -k A provides the same solution set as in the case of KM. The ω -k A tries to constitute a 3D Fourier transform relationship between the image at the object space and the collected scattered field. Before applying the FFT routine, a mapping procedure from frequency-wavenumber domain to wavenumber-wavenumber domain is necessary. Although this mapping procedure may slow down the execution time of the algorithm, it is still fast thanks to the FFT step. This study suggests that the focusing performance and the image quality feature of the ω -k A are the best among all five algorithms. The algorithm is also fast when compared to HSA, KM, and BPA.

(v) As the BPA is not conceptually as simple as HSA or KM, the implementation of the algorithm is more complex. On the other hand, the BPA is much faster than these two algorithms since it takes and acts on the data as a block. Therefore, it requires less computation resources when compared to HSA or KM. This study has showed that the BPA is quite fast and its focusing ability is also very good. This can be viewed from both the simulated and the measured images and from the listed tables where BPA is the best for almost all the focusing parameters. One drawback of the BPA is due to its vulnerability to numerical noise. The algorithm produces comparably higher noise floor when compared to other four algorithms. This situation may cause resultant images of the low reflectivity targets to lie under the exaggerated noise floor level. Therefore, such targets may not be imaged within selection of dynamic range for the image display.

Acknowledgment

This work is supported by the Scientific and Research Council of Turkey (TUBITAK) under Grant no. EEEAG-104E085.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] D. J. Daniels Surface-Penetrating Radar , IEEE Press, 1996.

[2] L. Peters Jr., J. J. Daniels, J. D. Young, "Ground penetrating radar as a subsurface environmental sensing tool," Proceedings of the IEEE , vol. 82, no. 12, pp. 1802-1822, 1994.

[3] S. Vitebskiy, L. Carin, M. A. Ressler, F. H. Le, "Ultra-wideband, short-pulse ground-penetrating radar: simulation and measurement," IEEE Transactions on Geoscience and Remote Sensing , vol. 35, no. 3, pp. 762-772, 1997.

[4] L. Carin, N. Geng, M. McClure, J. Sichina, L. Nguyen, "Ultra-wide-band synthetic-aperture radar for mine-field detection," IEEE Antennas and Propagation Magazine , vol. 41, no. 1, pp. 18-33, 1999.

[5] W. A. Schneider, "Integral formulation for migration in two and three dimensions," Geophysics , vol. 43, no. 1, pp. 49-76, 1978.

[6] J. Gazdag, "Wave equation migration with the phase-shift method," Geophysics , vol. 43, no. 7, pp. 1342-1351, 1978.

[7] R. H. Stolt, "Migration by Fourier transform," Geophysics , vol. 43, no. 1, pp. 23-48, 1978.

[8] E. Baysal, D. D. Kosloff, J. W. C. Sherwood, "Reverse time migration," Geophysics , vol. 48, no. 11, pp. 1514-1524, 1983.

[9] C. J. Leuschen, R. G. Plumb, "A matched-filter-based reverse-time migration algorithm for ground-penetrating radar data," IEEE Transactions on Geoscience and Remote Sensing , vol. 39, no. 5, pp. 929-936, 2001.

[10] K. Gu, G. Wang, J. Li, "Migration based SAR imaging for ground penetrating radar systems," IEE Proceedings: Radar, Sonar and Navigation , vol. 151, no. 5, pp. 317-325, 2004.

[11] J. Song, Q. H. Liu, P. Torrione, L. Collins, "Two-dimensional and three-dimensional NUFFT migration method for landmine detection using ground-penetrating radar," IEEE Transactions on Geoscience and Remote Sensing , vol. 44, no. 6, pp. 1462-1469, 2006.

[12] C. Cafforio, C. Prati, F. Rocca, "Full resolution focusing of Seasat SAR images in the frequency-wave number domain," International Journal of Remote Sensing , vol. 12, no. 3, pp. 491-510, 1991.

[13] C. Cafforio, C. Prati, F. Rocca, "SAR data focusing using seismic migration techniques," IEEE Transactions on Aerospace and Electronic Systems , vol. 27, no. 2, pp. 194-207, 1991.

[14] A. S. Milman, "SAR imaging by ω -k migration," International Journal of Remote Sensing , vol. 14, no. 10, pp. 1965-1979, 1993.

[15] H. J. Callow, M. P. Hayes, P. T. Gough, "Wavenumber domain reconstruction of SAR/SAS imagery using single transmitter and multiple/receiver geometry," Electronics Letters , vol. 38, no. 7, pp. 336-338, 2002.

[16] A. Gunawardena, D. Longstaff, "Wave equation formulation of synthetic aperture radar (SAR) algorithms in the time-space domain," IEEE Transactions on Geoscience and Remote Sensing , vol. 36, no. 6, pp. 1995-1999, 1998.

[17] Z. Anxue, J. Yansheng, W. Wenbing, W. Cheng, "Experimental studies on GPR velocity estimation and imaging method using migration in frequency-wavenumber domain," in Proceedings of the 5th International Symposium on Antennas, Propagation and EM Theory (ISAPE '00), pp. 468-473, Beijing, China, 2000.

[18] C. Gilmore, I. Jeffrey, J. LoVetri, "Derivation and comparison of SAR and frequency-wavenumber migration within a common inverse scalar wave problem formulation," IEEE Transactions on Geoscience and Remote Sensing , vol. 44, no. 6, pp. 1454-1461, 2006.

[19] J. Gazdag, P. Sguazzero, "Migration of seismic data," Proceedings of the IEEE , vol. 72, no. 10, pp. 1302-1315, 1984.

[20] C. Ozdemir, S. Demirci, E. Yigit, "A review on the migration methods in B-scan ground penetrating radar imaging," in Proceedings of the Progress in Electromagnetics Research Symposium (PIERS '12), pp. 789-793, Kuala Lumpur, Malaysia, March 2012.

[21] J. F. Claerbout Imaging the Earth's Interior , Blackwell Scientific Publications, 1985.

[22] C. Ozdemir, S. Demirci, E. Yigit, A. Kavak, "A hyperbolic summation method to focus B-scan ground penetrating radar images: an experimental study with a stepped frequency system," Microwave and Optical Technology Letters , vol. 49, no. 3, pp. 671-676, 2007.

[23] O. Yilmaz Seismic Data Processing , Society of Exploration Geophysicists, 1987.

[24] D. S. Jones The Theory of Electromagnetism , Pergamon Press, New York, NY, USA, 1964.

[25] P. M. Morse, H. Feshbach Methods of Theoretical Physics , McGraw-Hill, New York, NY, USA, 1953.

[26] I. Lecomte, S.-E. Hamran, L.-J. Gelius, "Improving Kirchhoff migration with repeated local plane-wave imaging? A SAR-inspired signal-processing approach in prestack depth imaging," Geophysical Prospecting , vol. 53, no. 6, pp. 767-785, 2005.

[27] J. A. Stratton Electromagnetic Theory , McGraw-Hill, New York, NY, USA, 1941.

[28] D. C. Munson Jr., J. D. O'Brien, W. K. Jenkins, "A tomographic formulation of spotlight-mode synthetic aperture radar," Proceedings of the IEEE , vol. 17, no. 8, pp. 917-925, 1983.

[29] J. L. Bauck, W. K. Jenkins, "Tomographic processing of spotlight-mode synthetic aperture radar signals with compensation for wavefront curvature," in Proceedings of the International Conference on Acoustics, Speech, and Signal Processing (ICASSP '88), vol. 2, pp. 1192-1195, New York, NY, USA, 1988.

[30] R. Mersereau, A. Oppenheim, "Digital reconstruction of multidimensional signals from their projections," Proceedings of the IEEE , vol. 62, no. 10, pp. 1319-1338, 1974.

[31] E. Fisher, G. A. McMechan, A. P. Annan, S. W. Cosway, "Examples of reverse-time migration of single-channel, ground- penetrating radar profiles," Geophysics , vol. 57, no. 4, pp. 577-586, 1992.

[32] L. Capineri, P. Grande, J. A. G. Temple, "Advanced image-processing technique for real-time interpretation of ground-penetrating radar images," International Journal of Imaging Systems and Technology , vol. 9, no. 1, pp. 51-59, 1998.

[33] I. L. Morrow, P. A. van Genderen, "2D polarimetric backpropagation algorithm for ground-penetrating radar applications," Microwave and Optical Technology Letters , vol. 28, pp. 1-4, 2001.

[34] C. Ozdemir Inverse Synthetic Aperture Radar Imaging with Matlab Algorithms , John Wiley & Sons, Hoboken, NJ, USA, 2012.

[35] A. Martinez, J. L. Marchand, "SAR image quality assessment," Revista de Teledetección , no. 2, pp. 1-7, 1993.

[36] E. Yigit, S. Demirci, C. Ozdemir, M. Tekbas, "Short-range ground-based synthetic aperture radar imaging: performance comparison between frequency-wavenumber migration and back-projection algorithms," Journal of Applied Remote Sensing , vol. 7, 2013.