GPR background removal using a directional total variation minimisation approach

Abstract

Ground penetrating radar (GPR) is a leading geophysical subsurface imaging tool for various purposes. This efficiency, however, is compromised by the interference of different types of noise. Background noise (clutter) is one of the nagging types of noise that undermines the high-resolution imaging capabilities of GPR. This study presents the experience of applying a directional total variation minimisation (DTVM) filter to attenuate clutter in GPR data. The application of DTVM to both synthetic and field GPR data proves its great capability to attenuate clutter without affecting the features of interest of a GPR section. The results also show that the proposed DTVM method affords superior image quality than both the most commonly used and the most recently published background removal techniques.

ground penetrating radar (GPR), clutter, image denoise, total variation (TV)

1. Introduction

Among the geophysical techniques for shallow subsurface imaging, ground penetrating radar (GPR) is the most common approach with a wide variety of applications. GPR has several geophysical and engineering applications, including groundwater, soils and mineral investigation. The main reason behind this broad range of applications and the growing use of GPR is its unique competence in acquiring high-resolution subsurface images in a relatively small time frame at reasonable cost. This proficiency, however, is weakened by the noise interference that degrades the high-resolution images provided by GPR and complicates the target identification process. The well-known background noise (clutter) is a typical artifact in almost all GPR images. This category of noise is presented as horizontal periodical straps that are sometimes so strong that they mask the genuine features of interest in the GPR section. Major source of clutter is the neighboring surface objects, especially those with metallic nature (Dojack 2012). Another provenance of clutter is antenna ringing. Antenna ringing can be caused by either the residual currents reverberating within the antenna or by the impedance inconsistency between the antenna and the ground surface (Daniels et al2008). Subsurface metals also appear as recurrent features that look like ringing but these are easily distinguishable because they are bounded by the object region (Radzevicius 2000).

Background removal is a routine procedure in GPR data processing. The most common techniques to diminish the background are the average trace subtraction techniques (Nobes 1999) and modified versions such as moving average trace subtraction. Despite their simplicity, these traditional background removal techniques are not always successful. Consequently, several different background removal techniques have appeared in the past two decades. Some of these techniques are based on mathematical algorithms such as domain filtering based on local differences between signal and noise (Young and Sun 1999). Deterministic deconvolution is used to improve the temporal resolution of GPR sections (Xia et al2003). A comparison study show that Radon transform (t-p) and discrete wavelets transform methods are superior to classical space-domain methods (Nuzzo and Quarta 2004). A fuzzy weighted background calculation through a sliding window is proposed in (Sezgin 2011). Very recently, alternative seismic stacking techniques (Rashed 2013) and background matrix subtraction (Rashed and Harbi 2014) were presented. Others are based on image processing techniques such as eigenimage processing (Cagnoli and Ulrich 2001). Eigenimage splitting using singular value decomposition is presented in Kim et al (2007). Jeng et al (2009) uses a multiresolution wavelet analysis to enhance the signal-to-noise ration GPR data.

The approach of total variation (TV) is a common signal processing technique that was introduced for the first time more than two decades ago (Rudin et al1992). TV is based on the principle of reducing the global variation of a signal, while preserving a close match to its original form. It has been used in several image processing applications such as image recovery and restoration (Chambolle and Lions 1997, Chan et al2000), pattern recognition (Vese and Osher 2003), image segmentation (Chan et al2006) and tomographic image reconstruction (Sidky et al2006). Although TV is a powerful approach in image denoising, it is hardly being used in GPR data processing. Very recently, a TV-based split Bregman method was proposed for GPR image restoration (Wang et al2013). Their results showed that a TV model of split Bregman can suppress the staircase artifacts while retaining the edges effectively. However, this method is mainly used to defeat low-frequency artifacts and might not be of the same effectiveness in cases where the clutter reaching amplitude values are very close to the original data, which is a truly common case in real GPR data.

In this study, a directional total variation minimisation (DTVM) algorithm is proposed. The DTVM algorithm is used to suppress the GPR clutter with minimal effect on the targeted features. Experimental studies using both synthetic and real GPR data are implemented to evaluate the proposed approach. The results are presented in comparison to those of the conventional average trace subtraction technique and to those of the most recently published background subtraction technique (Rashed and Harbi 2014).

2. Theory

2.1. Total variation and image denoise

The concept of TV minimisation is a powerful approach for image restoration (Caselles et al2011). A major advantage over conventional least squares methods in image denoising applications is the ability to preserve edges. For a two-dimensional image x, the TV-norm is defined as:

$${\displaystyle \parallel x{\parallel }_{\mathrm{T}\mathrm{V}}={\int }_{\mathrm{\Omega }}\mid \mathrm{\nabla }x\mid \mathrm{d}u\mathrm{d}v,}$$

1

where $\mathrm{\nabla }$ is the gradient, u and v are the orthogonal Cartesian coordinates and $\mathrm{\Omega }$ is the image domain. The so-called isotropic TV-norm can be calculated in the discrete form as:

$${\displaystyle \parallel x{\parallel }_{\mathrm{T}\mathrm{V}}\approx {\displaystyle \sum _u}{\displaystyle \sum _v}\sqrt{(x_{u+1,v}-x_{u,v}{)}^2+(x_{u,v+1}-x_{u,v}{)}^2+ϵ},}$$

2

where ε is a small positive number to ensure differentiability at points where $\mid \mathrm{\nabla }x\mid =0$⁠. In image denoising applications, the common general model of TV is formulated using the following cost function:

$${\displaystyle \underset{x}{min}f\left(x\right)=\frac{1}{2}\parallel x-y{\parallel }^2+\alpha \parallel x{\parallel }_{\mathrm{T}\mathrm{V}}}$$

3

where y is the measured noisy image and α is the regularisation parameter that balances the strength of the TV regularisation term. The cost function in equation (3) aims to find the solution image x^* that is close to the original image y in terms of the ${\mathrm{\ell }}_2$ norm distance and at the same time contains minimised TV value. The value of the parameter α can favour one term over another based on the target image and imaging application.

2.2. Problem formulation

Considering the GPR ringing noise, the following data model is used:

$${\displaystyle \mathbf{y}=\mathbf{x}+\mathbf{r}}$$

4

where $\mathbf{x}$ is the $m\times n$ matrix representing the original noise-free GPR data to be estimated, $\mathbf{y}$ is the acquired noisy data and $\mathbf{r}$ is the unwanted signal generated from different signal acquisition procedures. Model the incorrect data measurements $\mathbf{r}={\mathbf{r}}_n+{\mathbf{r}}_c$⁠, where ${\mathbf{r}}_n$ is random noise and ${\mathbf{r}}_c$ is clutter. To estimate $\mathbf{x}$ accurately, an interesting feature of the clutter ${\mathbf{r}}_c$ is considered here. The clutter is presented in almost consistent pixel values within the regions located within the analogical depth value (i.e. values located within the same horizontal line in GPR sections). However, the clutter has a nature of periodical value variation within the same trace (i.e. values located within the same vertical line in GPR sections). The objective is to extract the clutter ${\mathbf{r}}_c$ from the acquired data $\mathbf{y}$ by forcing the TV regularisation term over the horizontal direction. Then, a primary noise-free signal can be obtained by direct subtraction using equation (4). For simplicity, we ignore the effect of random noise (⁠${\mathbf{r}}_n=0$⁠) as the goal of this work is to mitigate clutter ${\mathbf{r}}_c$⁠. Now, consider dividing the matrix $\mathbf{x}$ into s sub-blocks such that $\mathbf{x}=[x_{1,n},x_{2,n},\dots ,x_{s,n}]$⁠, where each block x_{i, n} is a $q\times n$ matrix and $m=qs$⁠. The size of the image block can be varied based on the nature of the clutter. The corresponding matrices $\mathbf{y}$ and ${\mathbf{r}}_c$ are processed in the same way. Then, the cost function in equation (3) is solved for each matrix sub-block independently. This is achieved by solving the following optimisation problem:

$${\displaystyle \underset{r}{min}g(r)={\displaystyle \sum _{i=1}^s}g(r_i)}$$

5

$${\displaystyle g\left(r_i\right)=\frac{1}{2}\parallel y_i-r_i{\parallel }^2+{\alpha }_i\parallel r_i{\parallel }_{\mathrm{T}\mathrm{V}}}$$

6

where y_i and r_i are the i^th horizontal block of the matrices $\mathbf{y}$ and ${\mathbf{r}}_c$⁠, respectively. In other words, the main problem is divided into a set of s TV minimisation problems assigned to each image block. Each block contains the same number of columns and a smaller number of rows. The proposed algorithm is called directional TV minimisation (DTVM) to describe the approach of considering the TV minimisation through a single direction (horizontal direction). To solve the problem in equation (5), the lagged diffusivity fixed-point iteration (FP) method presented in Vogel and Oman (1996) is used. The FP iteration for DTVM algorithm is expressed as:

1 for$i=1:s$do

2  k  =  0

3  $r_i^{(0)}$  =  initial guess;

4  repeat

$d^{(k)}=-A(r_i^{(k)}{)}^{-1}g(r_i^{(k)})$⁠;

$r_i^{(k+1)}=r_i^{(k)}+d^{(k)}$

k  =  k  +  1;

untilstopping criterion

end

while the gradient g(r_i) is defined as:

$${\displaystyle g(r_i)=r_i-\alpha \mathrm{\nabla }\cdot \left(\frac{\mathrm{\nabla }{r}_i}{\sqrt{\mid \mathrm{\nabla }{r}_i{\mid }^2+ϵ}}\right)-y_i,}$$

7

and operator A(t) is given by:

$${\displaystyle A(t)t=(1+\alpha L_{ϵ}(t))t,L_{ϵ}(a)b=-\mathrm{\nabla }\cdot \left(\frac{\mathrm{\nabla }b}{\sqrt{\mid \mathrm{\nabla }a{\mid }^2+ϵ}}\right).}$$

8

This method is based on a quasi-Newton form, which tends to be less sensitive to round off error than the fixed-point form. The FP method is robust and proved to be globally convergent. The detailed computational algorithm can be found in Chan et al (2011).

3. Application

Evaluating the efficiency of a new processing technique requires applying this algorithm to both synthetic and field datasets and comparing its outcome to that of state-of-the-art techniques. Since there is a large number of existing background removal techniques, comparing the results of the proposed DTVM technique is limited to those of the most frequently used and the most recently published techniques. The most commonly used technique is the conventional average trace subtraction (ATS) technique, while the most recently published technique is the background matrix subtraction (BMS) technique (Rashed and Harbi 2014). One synthetic and two field GPR datasets are used for evaluation. These are the same datasets used by Rashed and Harbi (2014).

3.1. Synthetic data

The synthetic dataset is generated using 2.5D algorithm that uses finite-difference expression to solve Maxwell’s equations. The tested model is 25 m long and 15 m deep and was gridded each 2 cm in both directions. The model has three layers split by two interfaces dipping approaching the right-hand side of the section. Two ball-shaped objects are added in the middle and lower layers of the model, as shown in figure 1(a). The noise-free synthetic dataset is composed of 100 scans. The data are acquired using a source of Ricker wavelet at 200 MHz infinitesimal dipole in the 2D plane with a spacing of 0.25 m (figure 1(b)). The clutter presented as horizontal strips, of the same average amplitude, are generated using the identical approach employed to generate the synthetic clean data (figure 1(c)) and then merged with the clean data (figure 1(d)).

Figure 1.

Experiment data setup for synthetic data. (a) synthetic data velocity model, (b) clean data, (c) clutter and (d) data with clutter (summation of (b) and (c)). The colour scale in (d) and following figures represents the signal amplitude. Reprinted from Rashed and Harbi (2014), Copyright 2014, with permission from Elsevier.

The ATS, BMS and DTVM techniques were applied to the synthetic dataset and the results are illustrated in figure 2. The ATS technique is applied using a window width of 55 scans, while the BMS technique is applied using an identical window width, a trimming percentage of 25% and a weighting factor of 1.2. These are the same parameters used in Rashed and Harbi (2014). The DTVM technique was employed using 200 iterations with ${\alpha }_i=1(i=1,\dots ,s)$⁠, q  =  1 and $ϵ={10}^{-6}$⁠. Selecting the block size of a single row is the most extreme case, where the 2D image is processed as a set of m vectors.

Figure 2.

Results obtained from three background remove techniques. (a), (c) and (e) presenting a signal matrix computed using the ATS, BMS and DTVM techniques, respectively. The corresponding residual matrices are shown in (b), (d) and (f).

A careful inspection of figure 2 clearly shows the difference between the outputs of the three background removal techniques. Both the BMS and the DTVM techniques yield far better results than those of the ATS technique. The GPR section treated using the ATS technique has some streaks, especially on both sides of the peak of the superficial hyperbola. Streaks are also recognised in the proximity of the left end of the shallow dipping reflector as shown in figure 2(a). Moreover, the residual matrix resulted from subtracting the extracted clutter matrix using ATS from the synthetic background matrix shows many residues as presented in figure 2(b). These residues appear as strong horizontal streaks close to regions of high-amplitude features on the synthetic data. Moreover, residues of the dipping reflectors can be clearly seen on the residual matrix, indicating the extraction of a good portion of signal while removing clutter.

On the other hand, the two GPR sections treated using the BMS and the DTVM techniques look very similar, as shown in figures 2(c) and (e). The residual matrices, however, show significant differences. While the residual matrix of the BMS technique still has some distortion near time samples 80, 170, 230 and 280, the residual matrix of the DTVM contains mainly zeros (green on the colour scale) except near sample 280, which corresponds to the peak of the large hyperbola shown in figures 2(d) and (f). The DTVM residual matrix contains mostly zeros, indicating that the DTVM technique has succeeded in extracting the clutter. The extracted clutter is closely matched to the synthetic clutter with no or minimal effect on the signal we desire to preserve, as demonstrated in figure 2(f).

In addition to visual inspection, a measurement test is applied to check the similarity between the clean data, and the data resulted from applying the three background removal techniques. The correlation coefficient computed between the original clean matrix and the matrices resulted from applying the ATS, BMS and DTVM techniques are 0.9907, 0.9990 and 0.9997, respectively. Both the BMS and the DTVM matrices have a good correlation coefficient that is higher than that of the ATS technique. However, the data treated using the proposed DTVM technique are a closer match to the noise-free data than those treated using the BMS technique. Observing the results obtained from the two dipping reflectors using the DTVM method indicates that the proposed method works better with a steeply dipping reflector. This behaviour is expected since the proposed method discriminates against horizontal events.

It is interesting to evaluate the performance of the DTVM method with different imaging scenarios. Another synthetic data experiment is performed to test the validity of the proposed method with different dip angles of the reflector and different signal-to-noise ratios. Another experiment is conducted to demonstrate the effects of DTVM parameters such as α and the number of iterations k. To clearly show how the DTVM method works in different experiment conditions, a simple model is used. The model consists of two layers split by a single reflector with a dip angle ranging from 10° to 80° with an angle increment of 5° in each example, as shown in figure 3 (only the results of the 10 degrees increment are displayed for space adequacy). Horizontal clutter of the same average amplitude is generated and added to each model. The DTVM method is implemented with 1000 iterations and the other parameters are set to the same values as in the above experiment. It is clear from the output shown in figure 3 that the performance of the proposed method increases as the original signal is presented in the vertical direction (i.e. a dip angle close to 90°). When the original signal is presented in horizontal direction (i.e. a dip angle close to 0°), it is almost mixed with the clutter and extracting original signal is challenging. Quantitatively, the clutter reduction is measured using the relative root mean square error (RRME) defined as:

$${\displaystyle \mathrm{R}\mathrm{R}\mathrm{M}\mathrm{E}=\sqrt{\frac{{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}(x_{i,j}-x_{i,j}^{*}{)}^2}{{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}(x_{i,j}^{*}{)}^2}},}$$

9

where $\mathbf{x}$ denotes the output matrix computed using the DTVM method and ${\mathbf{x}}^{\mathbf{*}}$ is the corresponding original clean matrix. The RRME values computed from different dip angles are plotted in figure 5(a).

Figure 3.

Performance of the DTVM method with a different dipping angle. The left-hand column is the clean data with the dipping angle ranged from 10° to 80° (the dipping angle is shown in the top right-hand corner of each image). The middle column is the clean data mixed with clutter. The right-hand column is the corresponding output of the DTVM method.

Figure 4.

Assessment of the DTVM method with different values of the iteration number (left-hand column), α (middle column) and clutter amplitude corresponding to the clean signal (right-hand column). The parameter values are shown in the top right-hand corner of each image.

Figure 5.

Measured RRME with different values of (a) dip angle, (b) number of iterations, (c) regularisation parameter α and (d) clutter amplitude.

To trace the convergence behaviour of the DTVM method, the dip angle is fixed to ${45}^{°}$ and measures the RRME value with different iteration numbers. The output images are shown in figure 4 and the RRME plot is shown in figure 5(b). It can be observed that the output matrix becomes closer to the clean matrix as the iteration proceeds; however, it reaches stability after about 1000 iterations. Additional iterations would gain limited quality improvement. The effect of regularisation parameter α is also tested by fixing the number of iterations to 1000 and testing different values of α corresponding to ${10}^t,t=-3,-2,\dots ,3$⁠. The GPR sections obtained from different values of α are shown in figure 4 and the RRME plot is shown in figure 5(c). When α is relatively small, the effect of the TV is rather weak and consequently the separation between the original signal and the clutter does not work effectively. By increasing the value of α up to the value of 1.0 a significant signal separation is observed. Increasing α to higher values would not have an additional notable improvement. Finally, the value of $\alpha =1.0$ is used and the DTVM method is tested with different fractions of relative signal-to-noise ratios. Different fractions of average clutter amplitude are used corresponding to $1\times ,2\times ,3\times ,5\times ,8\times ,10\times \mathrm{a}\mathrm{n}\mathrm{d}16\times$ of the average clean signal amplitude. It is clearly shown from figure 4 that almost the same image is achieved even if the clutter average amplitude is 16 times the clean signal average amplitude. The RRME values are shown in figure 5(d) and they indicate that the change in clutter amplitude has a very limited effect on the performance of the DTVM method.

3.2. Field data

Synthetic data are important for evaluating a novel algorithm and for explaining how it works. However, a new processing algorithm is evaluated by the amount of enhancement it brings to the field data. The proposed DTVM filtering technique is tested using two field datasets. A SIR-3000 GSSI system with a 400 MHz antenna and transmit rate of 25 kHz is used to generate the first GPR dataset. The data measurements are acquired using 512 samples per scan and 50 scans per metre with low/high-cut filters of 100/800 MHz, respectively (figure 6(a)). The second field dataset is also acquired using the same equipment but with a 100 MHz antenna, a transmit rate equal to 100 kHz and a time range of 500 ns. The dataset is collected using a sampling rate of 1024 samples per scan and 50 scans per metre with a low-cut filter of 20 MHz and a high-cut filter of 400 MHz (figure 6(b)). A set of 400 scans of each dataset is used for the evaluation study, and the same processing routine is used for both datasets. More details on the modeling parameters of the synthetic data and field data specification can be found in Rashed and Harbi (2014).

Figure 6.

The first and second field GPR original data in (a) and (b), respectively.

Experiments conducted on the field GPR data confirm the results obtained earlier through the synthetic data experiment and further show the capability of the DTVM technique to suppress clutter with minimal effect on the targeted signal. The two field datasets shown in figure 6 are treated using the same three background removal techniques mentioned above and the results are displayed in figures 7 and 8.

Figure 7.

The first field GPR data computed using different background removing techniques. (a), (c) and (e) present the signal matrix computed using the ATS, BMS and DTVM techniques, respectively. The corresponding residual matrices are shown in (b), (d) and (f). The arrows (1)–(5) indicate identical regions in residual matrices computed using different techniques for comparison.

Figure 8.

The second field GPR data computed using different background removing techniques. (a), (c) and (e) present the signal matrix computed using the ATS, BMS and DTVM techniques, respectively. The corresponding residual matrices are shown in (b), (d) and (f). The arrows (1)–(3) indicate identical regions in residual matrices computed using different techniques for comparison.

The results of applying the ATS, BMS and DTVM background removal methods to the first field GPR data are illustrated in figure 7. Again, a quick look at the GPR sections treated with the three techniques shows that both the BMS and DTVM techniques yielded very similar results that are much better than those of the ATS technique, as shown in figures 7(a), (c) and (e). However, a careful look at the large hyperbola located in the middle region of the GPR section shows some difference between the results of the BMS and DTVM techniques. In the section treated using the DTVM technique, this hyperbola has a larger amplitude, more continuity and sharper edges than those of the BMS technique. These characteristics can be seen to a lesser extent on the smaller hyperbola located close to the right-hand edge of the section.

The extracted clutter shows more noticeable differences as demonstrated in figures 7(b), (d) and (f). While the clutter removed by implementing the ATS method is overwhelmed by feature signatures, the clutter matrix obtained using the BMS technique is less influenced by the high-amplitude features of the GPR section. On the other hand, the clutter matrix computed using the DTVM technique is originated purely from clutter with almost no effect on the high-amplitude features. The five marks (1)–(5) point to the same spots on the clutter matrices extracted using the three techniques. Mark (1) points to a spot that corresponds to the ringing shallow feature close to the left-hand border of the GPR section. Obviously, the clutter matrix computed using the ATS method is influenced by this high-amplitude feature, where the extracted ringing has both a curvature and amplitude change, which indicates that a good portion of the signal is extracted along with the clutter. This impact, however, does not exist in the two clutter matrices computed using the BMS and the DTVM techniques. The same description can be seen on mark (4), which corresponds to the small shallow hyperbola and mark (5), which points to the strong reflection at time sample 300 close to the right-hand border of the section.

Mark (2) points to a spot corresponding to the large hyperbola located in the middle of the GPR section. Clearly, the clutter computed using the ATS technique contains a notable portion of this hyperbola, while a smaller portion is included in the clutter computed using the BMS technique, as shown in figures 7(b) and (d). The clutter obtained using the proposed DTVM technique has no indication of influence of this high-amplitude feature, and the spot is mainly consisting of horizontal strips. The same can be said to the spot of mark (3) corresponding to the chain-like feature, which strongly influences the ATS results, and to a lesser extent the BMS results, while it has no effect at all on the DTVM results. It can be concluded that the implementation of the three background removal techniques to the first field GPR data agree with the findings of the synthetic data experiment and demonstrates the superiority of the DTVM technique to both the ATS and BMS techniques.

The implementation of the three background removal techniques to the second field dataset further confirms the above-mentioned results (figure 8). Although it is hard to see significant differences between the GPR sections extracted using different techniques, the extracted clutter matrices have clear differences. Again, high-amplitude features have a great impact on the clutter matrix computed using the ATS technique, a less but noticeable influence on the BMS clutter matrix and no influence at all on the DTVM matrix. Marks (1)–(3) corresponding to the fault location, the dipping basement and the high-frequency noise, respectively, clearly show the difference of these features on the computed clutter artifacts, as shown in figures 8(b), (d) and (f). Again, the superiority of the DTVM technique over both the BMS and the ATS techniques is clear in this field GPR dataset.

4. Summary and conclusions

This study presents an original approach to attenuate clutter in GPR data using directional total variation minimisation (DTVM). The main advantage of the new background removal technique is its capability to attenuate clutter without affecting features of interest. Synthetic and field GPR data are used to evaluate the proposed technique. The results indicate that the proposed method is effective in reducing clutter even with a relatively high amplitude. The original features presented in the vertical direction are more preserved than the features presented in the horizontal direction. The experimental results are discussed in contrast to those of the ATS and BMS techniques. The ATS technique is the typical background removal technique, while the BMS technique is the most recently published one. The results indicate that the proposed DTVM technique achieves superior results than either of these techniques.

A major reason behind the poorer results of the ATS technique is that it is based on a simple moving averaging process. Consequently, anomalously high and low amplitude values ruin the averaging process and cause the ATS technique to remove a good portion of the signal along with the clutter. To avoid this problem, the BMS technique uses a series of trimming, exclusion and weighting of these values, but it is basically based on the same sliding window mechanism. That is why the BMS technique reduces the effect of the anomalous values but does not completely terminate it. The DTVM technique, however, uses a totally different approach to deal with the clutter problem. It deals smartly with the natural characteristics of the ringing noise, which is consistent within the same depth. The DTVM algorithm enforces minimising the TV of the measured data to suppress primary signals while keeping the clutter. Then a noise-free signal can be extracted through a direct subtraction. This explains the strong capability of the DTVM filtering technique to attenuate clutter artifacts in a satisfactory manner without affecting the features of interest. In conclusion, the proposed DTVM background filtering technique is capable of removing most clutter with minimal effect on the hyperbolic and dipping features of interest, beating both the most commonly used ATS technique and the most recently published BMS technique.

Acknowledgments

The author would like to thank Prof M Rashed and Dr H Harbi of King Abdulaziz University for providing the GPR data and for their kind support during the course of this work.