Kirchhoff and F-K migration to focus ground penetrating radar images

Ground penetrating radar (GPR) based land mine detection has a main challenge of having an accurate image analysis method that is capable of reducing false alarms. However this image analysis depends on having sufficient spatial resolution in the backscattered signal. This paper aims at getting better resolution by applying two migration algorithms. One is by Kirchhoff’s migration using geometrical approach and other one is F-K migration algorithms with Fourier transform. The algorithms are developed using MATLAB simulations over different scenarios for stepped frequency continuous wave (SFCW) GPR.

Use of non-uniform fast fourier transform (NUFFT) algorithm to migrate GPR data as discussed in [9] has resulted in high signal to noise ratio in the migrated data. The nonuniform nature of the wavenumber space requires linear interpolation before the regular fast Fourier transform (FFT) could be applied. However, linear interpolation usually degrades the quality of reconstructed images. The NUFFT method mitigates such errors by using high-order spatial-varying kernels. The phase-shift migration (PSM) method is first introduced and applied by Gazdag [7] to handle lateral velocity variation. This 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 performance of focusing algorithms over different parameters are compared over various migration algorithm shows F-K migration is the fastest one and the Back-Projection based Migration (BPA) is the second best while the Kirchhoff migration is the slowest [4]. From the signal to clutter ratio results, the clutter suppression ability of BPA is the best. Synthetic aperture radar (SAR) based focusing techniques have been applied to GPR data in [13] based on F-K migration. It also suggests that F-K domain SAR algorithm better focuses than hyperbolic summation method and also more efficient in computation time.

Kirchhoff migration algorithm
Data collection in GPR can be done through CMP (common mid-point) gathering and then applying the signal traces for NMO (Normal move out) correction and stacking it. The Fig. 1 depicts it. CMP gathering is one of the techniques of collecting the data showcasing the subsurface of the earth through the GPR system. The Fig. 2 shows how the data is captured using CMP gathering technique [5].
The transmitter and receiver must me placed such that for all positions of transmitter and receiver for a particular shot gather, the midpoint of the offset must remain the same. Here the offset refers to the distance between the transmitter and receiver. For each shot, we obtain the reflection which is referred to as the A-scan and the collection of such several A-scans gives the 2-D B-scan image which is used for processing. The B-Scan image turns out to be the hyperbola.
Non-zero offset data is characterized by a travel time increase with increase in offset distance from the source to the reflector. The non-zero offset to zero offset conversion is achieved through a correction called NMO (normal move out) correction. For the single constant horizontal velocity layer the trace time curve of a function of offset is a hyperbola.
The time difference between travel time at a given offset and at zero offset is called NMO. The velocity required to correct for normal move out is called the normal move where x 2 is the distance (offset) between the source2 and receiver2 position in Fig. 2. v 2 is the velocity of the medium above the reflecting interface. t 2 (0) is twice the travel time along the vertical path. Generally, the NMO correction is given by the difference between t(x) and t(0) Each CMP gather after NMO correction is summed together to yield a stacked trace. Stacking enhances the in-phase components and reduces the random noise. It yields zero offset section (in the absence of dipping layers in the subsurface) [5].
Migration is the process that moves the reflection energies from the apparent locations to the true locations. Kirchhoff migration is a non-recursive method. It uses integral form of wave equation and based on Huygens principle, according to which, the seismic reflector is viewed as if it is composed of closely placed point diffractions. Apex of the diffraction curve to be the location of the true point reflector. The purpose of Kirchhoff migration is to sum up the energy produced by every Huygens secondary source and map it into point of generation as shown in Fig. 3.
In this way, each point on the migrated section is treated independently from the other points. Each point on the output migrated section is produced by adding all data values along a diffraction that is centered at that point. The diffraction summation method and CC' respectively [5].
In Kirchhoff 's migration procedure, the aim is to find the solution to the scalar wave equation for the wave function (x, z, t). Solving the wave equation, the far-field approximation considered in Kirchhoff migration can be expressed as: Above equation clearly contains the obliquity (cos ∅) and the spherical divergence (1/r) factors among its terms. The second term in this equation, (1/r) u(r 0 , t 0 ) is usually dropped, because it is proportional to 1/r 2 [4]. Figure 5a and b shows the geometry considered for performing Kirchhoff migration. Consider scanning of surface is done beginning at point O. At point A transmitter-receiver pair gets a reflection from the dipped reflector or the object. By continuing scanning in x direction, assume that the reflection ends at point B. Now, we have the apparent position of the object as C'D' . It is necessary to apply the migration algorithm on this geometry of Fig. 5a to get true position of the object as CD [5].
The length BD' and AC' are known from reflected time. Also, from geometry in Fig. 5a, Angle BDO = 90° and Angle ACD = 90°. Consider the origin O(x o , z o ), location Thus true location is P(x p , z 0 ) and apparent location is A(x a , z 0 ) over the ground surface as in Fig. 5b. C'(x a , z a ) over depth of z a from ground is the apparent position and C(x p , z p ) is the true location over a depth of z p from ground surface.

F-K migration algorithm
Migration of dipping event in the F-K domain is shown in Fig. 6. The vertical axis is the frequency axis ω for pre-migrated event and it is the vertical wavenumber k z , associated with the depth axis for post-migrated event. We assume velocity equal to 1. F-K migration maps lines of constant frequency AB in the (ω, k x ) plane to circles AB' in the (k z , k x ) plane. Therefore, migration maps point B vertically onto point B' . Note that in this process, the horizontal wave number k x does not change as a result of mapping. When this mapping is completed, the dipping event OB is mapped along OB' after migration. Thus, the dip angle after migration (θ') is greater than before migration (θ). For comparison, these two radial lines are shown on the same plane (k z , k x ), Fig. 6.
The phase shift method can only handle vertically varying velocities. To handle lateral velocity variations, first the input wave field is extrapolated by the phase-shift method using a multiple number of laterally constant velocity functions and a series of reference wave fields are created. The imaged wave field then is computed by interpolation from the reference wave fields. This migration method is known as phase-shift-plus-interpolation. If the medium velocity is constant, then migration can be expressed as a direct mapping from temporal frequency ωω to vertical wavenumber k z.

Mathematical model
The Fourier like transform is taken for the obtained back scattered field e(x, y, z, t) where z = 0. Since the reflected EM waves are measured at the surface level i.e. zero depth which gives the following result as shown in the equation.
The inverse Fourier transform for the given signal is given as: Now, putting t = 0 in the above equation and introducing the k z term by substituting the results of the following equation, e(x, y, t)e j(k x x+k y y−ωt) dxdydt ϕ(x, y, 0, t) e x, y, t = 1 2π The above equation presents a focused image. However, the data in (k x , k y , w) domain should be transformed to (k x , k y , k z ) domain to be able to use the FFT. Hence the mapping is done in order to do this transformation. The value of k z is found out by using the following relation, Solving this equation we get: Now k z value is found out for every corresponding k x and w and substituted in the following equation; After the above procedures have been carried out, the inverse Fourier Transform is taken to obtain the migrated or focused image.
All the above equations are for C-scan data migration. If a B-scan data has to be processed then the y and the k y terms should be dropped respectively from the above equations. The resulting equation in that case will be as under: Although the mapping procedure slows down the algorithm, yet it is speedier than any other migration algorithm used for humanitarian demining. This is due to the use of FFT for fast computation of the data.
Step by step description of F-K migration algorithm 1. Obtain the time traces of wave-field amplitudes E(x, z = 0, t) 2. Discretize E(x, z = 0, t) and obtain MxN matrix 3. Take 2-D FFT on E(x, z = 0, t) 4. Interpolate E(k x , f ) to obtain desired depth E(k x , k z ) 5. Apply 2-D IFFT on E(k x , k z ) and select absolute values to get focused image.

Simulation results
The proposed algorithms are implemented in MATLAB 2014a on Intel core i3 processor at a speed of 3.5 GHz on windows 7 operating systems. B-scan data is simulated for different positions of targets. Then F-K and Kirchhoff migration algorithms are applied on e x, y, z, 0 = 1 2π 3 2 ∞ −∞ E(k x , k y , ω)e −j(k x x+k y y+k z z) dk x dk y dω these B-scan data, assuming ground to be homogeneous, without any clutter, to get true locations of targets. We have considered multiple targets at a depth of 100 to 150 cm. A stepped frequency input waveform of bandwidth 3 GHz and divided it into 256 equal steps is used. Taking more spatial points results in better accuracy but increases computational time. At each spatial point, the frequency is varied from 200 MHz to 3.06 GHz and back scattered electric field data is stored. During simulation for the B-scan generation, the following parameters were used: the relative permittivity of ground medium was chosen as 9. The backscattered electric field was collected along a straight path ranging from x = 0 to x = 6m showing 256 points on range. The data was simulated by taking origin as centre of the x-axis and also results are indicated with respect to pixel co-ordinates. Physical dimensions are calculated as (co-ordinate number*resolution). Resolution is considered as (total distance/number of sample points). The cross range (depth) is ranging from 0 to 2 m and is also divided into 256 points. The A-scan required for testing was simulated by considering the target co-ordinates and calculating the distance of SFCW GPR. The data consists of in-phase (I) and quadrature (Q) values, which corresponds to real and imaginary parts of DFT. The time varying radar signal is recaptured by taking inverse DFT of I-Q data. Then successive A-scan data is stacked by taking new distances to yield hyperbola called B-scan. The B-scan GPR image in spatial domain is obtained by taking the Inverse FFT of the back scattered signal data in the frequency domain as shown in Fig. 7a and b. The raw GPR image of the simulated object exhibits the hyperbolic defocusing behavior. As a case 1, B-scan data was generated considering three targets having locations numbered as 1, 2 and 3 and having true coordinates (50, 105) (100, 115) and (80, 110) respectively resulting in Fig. 7a. Figure 8a shows the focused image of migrated data after applying Kirchhoff migration algorithm, resulting in output data with migrated targets in locations (50, 105) (100, 115) (80, 110) respectively. This resultant image is successfully focused. The target is located at the apex of hyperbola which is converging into a single point at the location of apex of hyperbola for each of the target, after applying our Kirchhoff migration algorithm.
The Table 1 shows results of Kirchhoff migration for multiple targets for the cases discussed above. It includes the targets numbered as 1, 2 and 3 of Fig. 8a which shows the true and migrated target positional co-ordinates and also physical dimensions (lateral and buried depth) in cm. The Table 1 also includes the targets numbered as 4, 5 and 6 considered in Fig. 8b, showing the error in depth as −0.38 and −0.785 cm lateral error. Average execution time of Kirchhoff migration algorithm in Matlab is 3.45 s.
For the second set of simulation results as a case 3, we have considered same B-scan data set shown in Fig. 7a as used for Kirchhoff migration algorithm. Then, we apply the F-K migration algorithm in aiming to get a better focused image. The resultant migrated image is displayed in Fig. 9a, output data with migrated targets 1, 2 and 3 are in locations (50, 103) (100, 112) (80, 108) respectively. Where in the targets responses are seen to be more localized around their correct locations.
As a case 4, the Fig. 7b is considered as input to F-K migration algorithm. Figure 9b shows the focused image of migrated data after applying F-K migration algorithm, resulting in output data with migrated targets of 4, 5 and 6 in locations (18, 155) (53, 119) (66, 136) respectively.
The Table 2 is a result of F-K migration on multiple targets which shows that there is no error in lateral position and −2.34 cm is the maximum error in depth is for the targets numbered as 1, 2 and 3 considered in Fig. 9a. The Table 2 also includes the targets numbered as 4, 5 and 6 as shown in Fig. 9b, with maximum error of −2.34 cm in depth and maximum error of 2.3 cm in lateral positions. Average execution time of F-K migration algorithm in Matlab is 2.05 s. This shows F-K migrated data gives better results than Kirchhoff algorithm. The Kirchhoff method results in slightly more error in depth and also takes more computational time.

Conclusions
We have considered different cases for SFCW GPR with multiple targets, homogeneous ground conditions and no clutter. Both Kirchhoff and F-K migration methods showed promising results in improving resolution of GPR in focusing the images. The Kirchhoff migration is done after doing hyperbolic summation as NMO and stacking. Maximum error in lateral position is of one coordinate corresponding to 2.3 cm and maximum error in depth position is of 3 coordinates corresponding to 2.34 cm in F-K migration. This indicates that the actual and migrated depths are very close to each other. Under various scenarios, both migration algorithm provide very less error, hence these methods can be implemented on data from GPR in all cases. In F-K migration algorithm, 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. Here, F-K migration algorithm takes an average time of 2.05 s and Kirchhoff migration takes 3.45 s. It clearly says that F-K migration algorithm results in less computation time than Kirchhoff migration algorithm and also results in slightly less error in lateral and depth position.