**Bello Y. Idi**

Department of Geomatic Engineering, FKSG, Universiti Teknologi Malaysia, Johar, Malaysia

Md. N. Kamarudin

Institute of Geospatial Science and Technology, Universiti Teknologi Malaysia, Johor, Malaysia

Department of Geomatic Engineering, FKSG, Universiti Teknologi Malaysia, Johar, Malaysia

Md. N. Kamarudin

Institute of Geospatial Science and Technology, Universiti Teknologi Malaysia, Johor, Malaysia

The ultra wide bandwidth nature of ground penetrating radar antenna has made a raw data acquired with the tool prone to unwanted noise and hence low signal to noise ratio. Quantitative interpretation of the data is desirable for radar image quality enhancement. This study used multiresolution analysis to process radar image at different levels of decomposition. Daubechie wavelets family was used to decompose the image into 4-different levels of details. Level 3 diagonal details and level 4 horizontal details provide a noise-free visualization of the subsurface discontinuities. This led to the detection of buried utility and unique identification of its spatial location. The depth of the buried utility was estimated based on the image scale. The work demonstrates the effectiveness of Daubechie wavelet transform analysis as yet another technique of utility survey data acquired with GPR.

PDF Abstract XML References Citation

Bello Y. Idi and Md. N. Kamarudin, 2012. Interpretation of Ground Penetrating Radar Image Using Digital Wavelet Transform. *Asian Journal of Applied Sciences, 5: 174-182.*

**DOI:** 10.3923/ajaps.2012.174.182

**URL:** https://scialert.net/abstract/?doi=ajaps.2012.174.182

Multiresolution analysis is one of the most applicable techniques that are effectively useful in image processing. It involves decomposition of image data into several scale resolutions with the aim of analyzing localized area of significance. The analysis is practically implementable using wavelet transform. **Wavelet transform** refers to the decomposition or breaking up of signal into shifted and scaled version producing series of details and approximations at various levels. Therefore, it enables the signal to be analyzed at a wide range of resolutions leading to information extraction from signal that has widely different scale components. The approximation component of the decomposed signal represents the low frequency component and is obtained by passing the signal through a low pass filter which removes the high frequency component of the signal and therefore results into approximate variation (approximate information parameters). The detail component on the other hand, represents the high frequency component and is obtained by passing the signal through a high pass filter which suppresses the low frequency and maintains the high frequency component. The high frequency component is associated with instantaneous signal information parameters. The technique is widely useful in image analysis of non stationary signals like seismic, earthquakes etc. (Davoodi *et al*., 2009), in which detail information could be extracted which other signal analysis could miss. It however has a wider range of application as it was successfully found to be effective in various signal and **image processing** analysis including signal denoising (Jeng * et al*., 2009), object detection (Semenchuk and Lukyanov, 2009), signal compression (Brechet *et al*., 2007), discontinuity detection (Rossini, 2009), statistical analysis (Antoniadis, 2007) etc.

Wavelet analysis basically involves recursive decomposition of signal by shifting and scaling a wavelet (waveform of limited duration) over the signals. It is mathematically defined as the sum of the signal x(t) multiplied by scaled and shifted version of the wavelet ψ_{jk}. The resulting product is a set of constants C_{jk} known as wavelet coefficients which are functions of scales and positions j, k, defined as (Ni *et al*., 2010):

(1) |

The wavelet coefficient at an instance is the measure of the correlation between the wavelet and a given segment of the signal. The analysis can be performed at a continued scales and positions within the various segments of the signal. This is known as Continuous Wavelet Transforms (CWT). It could however be more efficient if the scales and positions are chosen based on a powers of two, a process known as dyadic scales and positions. Such analysis is refers to as Discrete Wavelet Transform (DWT). The discrete wavelet is a quantized version of scales and positions parameters of the CWT in which the wavelet ψ_{jk} is discretely chosen on a dyadic grid intervals given by Ni *et al*. (2010):

(2) |

The integers j and k control the scaling and the positioning. The scale is a compressing factor of the wavelet which is related to signal frequency. Lower scale implies compressed wavelet, leading to high frequency component of the signal. Higher scale on the other hand leads to wavelet stretchening, resulting to low frequency component. Low and high frequencies respectively represent the approximation and detail components of the signal.

In practice, the approximations and details are obtained by filtering process in which the signal X(t) passes through complimentary high pass and low pass filters simultaneously. The approximation is obtained as the output of the low pass filter while the detail is the output of the high pass filter. The process is depicted in Fig. 1a showing signal decomposition described by the equation:

(3) |

The process can iteratively continue until a suitably chosen number of levels are reached depending on the level of details and approximations necessary to achieve the required objective. Figure 1b is a decomposition tree for a signal X(t) at level 3 showing the approximation A_{i} and details D_{i} for the three levels.

Image processing in DWT implies two-dimensional wavelet analysis which recursively leads to decomposition of the image into series of approximations and details at different levels of resolutions. This enables the possibility of viewing the image at different bands of frequencies thereby enabling finer and better spectral characteristics of the target anomaly to be identified.

Fig. 1a: | Signal decomposition into Approximation A_{j}(t) and detail D_{j}(t) |

Fig. 1b: | 3-level decomposition tree |

Thus variation in wavelet coefficients due to spectral characteristics of anomalous object with respect to the background constituents enhances the possibility of detecting and recognizing the obscured objects from the image.

Ground Penetrating Radar (GPR) is a near surface geophysical tool that record the back scattered signal of the subsurface reflected due to contrast in the electrical properties of the earth composition. The 2-D output data of GPR is a radar image of the plot of variation in amplitudes of the reflected signal at various locations as a function of radar travelling time. The image consists of stacked signal amplitudes at points thereby giving a good resolution depending on the number of samples recorded per scan and the magnitude of the contrast between the reflector and its surrounding host soil. GPR is one of the nondestructive and noninvasive techniques applicable in various fields for the propose of soil mapping, road condition survey, fault detection, utility mapping etc. It is analogous to ultrasonic wave nondestructive testing technique in which the variation in the wave propagation parameters through a heterogeneous medium can be used for *in-situ* monitoring of the relative variation in the **physical properties** of the medium (Rhazi and Kodjo, 2010).

The raw radar image is however characterized by low signal-to-noise ratio mainly due to the ultra-wide bandwidth nature of the receiving antenna and high degree of subsurface inhomogeneity. The accuracy of information from GPR radargram therefore depends on the effectiveness of the signal processing technique adopted. Several **image processing** algorithms are applicable to GPR data (Persico *et al*., 2010). These are categorized into basic processing which mainly deals with image enhancement and signal-to-noise ratio improvement and advanced processing, which is user bias, that depends on the information required (Sensors and Software, 1999). In utility mapping, various signal and **image processing** techniques are adopted in an attempt to successfully distinguish between the signals reflected by the buried utility pipe from that reflected by the surrounding soil (Bonomo *et al*., 2011).

Multiresolution analysis was found to be very effective in processing non stationary signals and therefore adopted as a high resolution processing technique in radar, remote sensing and geophysical imagery. Analysis of detail component of decomposed subsurface image using wavelet transform was carried out by many authors for the purpose of identification and spatial location of subsurface anomaly. For instance wavelet multi scale was used by Xu *et al*. (2009) in separating regional gravity anomaly resulting from deep structure in Dagang area of Northern china. Large multi scale decomposition based on the spectral analysis was used to obtain regional anomaly from the gravity data and used to estimate the depth from the radial power spectrum. The wavelet analysis and the spectral density were used to establish the relationship between the separated anomalies and the buried geologic units. Tomiya and Ageishi (2003) used the combination of Haar and Daubechie’s wavelet analysis in image registration of remotely sensed data. They discovered that the wavelet transform method is more stable and superior than Fourier transform method. Ni *et al*. (2010) used Meyer wavelet in frequency domain to filter and enhance raw GPR radargram leading to the detection of buried pipes under different conditions.

In this study, Daubechie’s wavelet was used to filter and enhance radar image acquired over a site of known buried utilities. The objective of the work is to observe and evaluate the appearance of the spectral signature of the utilities resulting from the Daubechie wavelet transform.

**MATERIALS AND METHODS **

**Field procedure and data processing:** The study area is a football pitch of Kolej 16, in the main campus of Universiti Teknologi Malaysia, Skudai, Johor Bahru, Malaysia. A common offset single fold reflection profiling was used to obtain a GPR cross section of length 20 m at a surface of the pitch. The data was acquired using IDS DAD fast wave radar acquisiti on unit at a fixed centre frequency of 600 MHz. A profile was chosen across a buried metal pipes at approximate depth of about 0.75-0.85 m. The profile was scanned at a trace increment of 0.025 within 100 nsec^{-1} time window. A parallel broadside antenna orientation was used for the operation. Figure 2 shows the raw data obtained.

The acquired data image was processed with Reflex 2D GPR/seismic processing software. Temporal filtering (DEWOW) was applied to remove very low frequency components from the image data associated with possible instrumentation dynamic range or inductive phenomena. Static correction was applied in order to time-shift the data record to first significant arrivals. A bandpass butherworth filter was applied in time domain recursively in order to suppress noises that differ in frequency content from the signal. Temporarily consistent noises were eliminated with background removal.

Fig. 2: | Raw radargram |

Fig. 3: | Processed GPR radargram |

Four hyperbolic signatures are clearly visible in the raw radargram (Fig. 2). Hyperbolic signatures in radar image are usually associated with point or cylindrical reflectors of relatively smaller dimensions such as utility pipes. Velocity adaptation by hyperbolic curve fitting was used to estimate the magnitudes of signal velocity of the medium above the reflectors. A mean velocity value of 0.102 m n sec^{-1} was obtained and used to convert the time axis to depth. The mean velocity value was also used for f-k migration of the radargram in order to collapse the diffraction hyperbolas to their causative source. The processed and migrated image which serves as the input data for the wavelet analysis is shown in Fig. 3.

**Daubechie wavelet analysis:** The Daubechie wavelet, a revised and generalized version of the Haar transform also involves recursive decomposition of signal into average as approximation and differencing as the details. Its advantage over Haar is that it is easier to implement and invert. The transformation has four-level wavelet scaling functions or coefficients which are implemented as a succession of decomposition into two pairs of average c(n) and difference d(n) given by Addison (2002):

(4a) |

(4b) |

Where:

In this type of wavelet, if the input signal x(n) is linear, the local average c(n) will also be linear but the local difference d(n) will be zero. This provides sparse presentation for linear signal making the transform effective for noise reduction.

The process is accomplished in the case of **image processing** by applying 1-D wavelet transform to each row of the pixel values in the image resulting to a horizontal approximation and details. It is then reapplied to each column resulting to vertical approximations and details (Afshari, 2011). Events of specific scale are filtered by threshold. This involves decision making based on the signal in excess of a predefined threshold value. Decision is taken on the presence or absence of a suspected object based on whether the reflected signal (or pixel value) is above or below the threshold value. In most of non stationary signals like GPR, energy is mostly concentrated in a smaller number of wavelets associated with large contrast is reflectance characteristics compared to the remaining background and other regions of local inhomogeinity. These regions are therefore characterized by larger wavelet coefficients while the background noise energy is spread over a larger number of coefficients. The background noise can therefore be eliminated optimally by thresholding the smaller coefficients to zero. The detail coefficients d(n) which are related to the amplitude of the signals in this case are threshold to a new value d’(n) according to the following technique:

(5) |

where, the threshold θ is chosen based on observation of the noise characteristics.

The processed radar image is formatted to a portable Jpeg image with each data point of the section made up of image pixel having a value proportional to the signal intensity of the reflected EM pulse (radar reflectance energy). The image was analyzed using 2-D MATLAB wavelet tool box. Four-ordered Daubechie (db-4) wavelet was computed for both the horizontal, vertical and diagonal details. The signals were threshold at different levels according to Eq. 5, iteratively with the intensity value below the threshold set to zero. This suppressed the wavelet coefficients smaller than the desired amplitudes while preserving significant reflectance associated with stronger anomalies. The visibility of the image is therefore optimized enabling clearer observation of the location of causative bodies.

Figure 4 is the square-view mode of the decomposed radargrams. The original input image shown on the upper-right corner of the view, was decomposed into 4-level decompositions labeled horizontal (H1-H4), vertical (V1-V4) and diagonal (D1-D4) as shown. Each of the segments is a plot of the detail coefficients and therefore represents the image detail at unique resolution. Two most significant levels that yield optimal information are level 3 diagonal details (D3) and level 4 horizontal details (H4) shown in Fig. 5. The images show bright sparks of different colors within the active regions.

Fig. 4: | A window snap of the 4-level horizontal, vertical and diagonal decompositions of the radar image |

Fig. 5(a-b): | (a) Level 3 diagonal detail and (b) Level 4 horizontal detail |

Three brighter sparks occur within the depth of 1.5 m, the spots are marked A, B and C as shown in Fig. 5a and b. These sparks correspond to the three hyperbola curves that occur within the same depth range around 4.0, 9.5 and 14.0 m positions in the processed radargram (Fig. 3).

Anomaly identification from radar image using wavelet transform is based on the fact that radar signal energy is mostly concentrated in a small number of dimensions of relatively large coefficients. This is in contrast to the wavelet coefficients of the background soil whose energies are spread over larger number of coefficients. The transformed details images threshold the smaller coefficients to zero thereby optimally revealing the most important information related to strong signal reflectors.

Figure 4 shows how the background information are truncated while the information that satisfy the conditions set up in Eq. 5 are preserved. Three anomalously strong reflectors shown in Fig. 5a and b are associated with the hyperbolic signatures of Fig. 3 and are therefore likely point reflectors or dykes. The most significant of the three labeled A, coincides with the buried pipe and appears at approximate depth of 0.7 m. The point appears reddish at H4 detail (Fig. 5b) indicating relatively large contrast in wavelet coefficients. This implies correspondingly large contrast in reflective radar signal energy and thus a significant contrast in electrical property between the reflector and the surrounding soil. The remaining two sparks appear at relatively less significant color variation although still discernible. They also indicate the subsurface inhomogeinity that could be associated with variation in physical, electrical or hygrogeologic properties.

Although the depth of coverage is about 5 m, significant events are limited to a depth of about 2 m only. This implies that deeper features are masked due to signal attenuation with depth. It could be as a result of the increase in **moisture content** with depth which enhances the **electrical conductivity** of the soil and consequently attenuate the signals. The work demonstrates the effectiveness of Daubechie DWT analysis in extracting information about the spatial location of subsurface discontinuities from GPR image. The technique can thus be used to estimate the location and depth of buried utilities with appreciable level of accuracy.

In this study, we appraised the application of multiresolution analysis as processing technique in utility mapping with GPR. A 2D radargram recorded across a test side was processed, enhanced and analyzed using Daubechie DWT. The output images, which are plots of the coefficients of radargram at four different levels of decomposition, were used to identify the spatial location of discontinuities within the subsurface. Buried utility was uniquely identified and its spatial location was estimated from the image of the detail coefficients. The work therefore demonstrates the effectiveness of DWT technique in radar image interpretation.

The authors acknowledged with thanks the contribution of Jurukur Abadi of No. 06-01, Jalan Padi Emas 4/5, Pusat Bandar Tampoi, 81200, Johor Bahru for the test data acquisition.

- Afshari, M., 2011. An algorithm to analyze of two-dimensional function by using wavelet coefficients and relationship between coefficients. Asian J. Applied Sci., 4: 414-422.

CrossRef - Antoniadis, A., 2007. Wavelet methods in statistics: Some recent developments and their applications. Stat. Surv., 1: 16-55.

CrossRefDirect Link - Bonomo, N., M. de la Vega, P. Martinelli and A. Osella, 2011. Pipe-flange detection with GPR. J. Geophys. Eng., 8: 35-45.

CrossRef - Brechet, L., M.F. Lucas, C. Doncarli and D. Farina, 2007. Compression of biomedical signals with mother wavelet optimization and best-basis wavelet packet selection. IEEE Trans. Biomed. Eng., 54: 2186-2192.

CrossRefPubMedDirect Link - Davoodi, M., M.A. Sakhi and M.K. Jafari, 2009. Comparing classical and modern signal processing techniques in evaluating modal frequencies of Masjed Soleiman embankment dam during earthquakes. Asian J. Applied Sci., 2: 36-49.

CrossRefDirect Link - Jeng, Y., C.H. Lin, Y.W. Li, C.S. Chen and H.H. Huang, 2009. Application of multiresolution analysis in removing ground-penetrating radar noise. Frontiers + Innovation, CSPG CSEG CWLS Conventions, Calgary, Alberta, Canada, pp: 416-419. https://www.cseg.ca/conventions/abstracts/2009/2009abstracts/127.pdf.
- Ni, S.H., Y.H. Huang, K.F. Lo and D.C. Lin, 2010. Buried pipe detection by ground penetrating radar using the discrete wavelet transform. Comput. Geotech., 37: 440-448.

CrossRefDirect Link - Persico, R., F. Soldovieri and E. Utsi, 2010. Microwave tomography for processing of GPR data at Ballachulish. J. Geophys. Eng., 7: 164-173.

CrossRef - Rhazi, J. and S. Kodjo, 2010. Non-destructive evaluation of concrete by the quality factor. Int. J. Phys. Sci., 5: 2458-2465.

Direct Link - Rossini, M., 2009. Detecting discontinuities in two-dimensional signals sampled on a grid. J. Numer. Anal. Ind. Applied Math., 4: 203-215.

Direct Link - Semenchuk, E.V. and S.P. Lukyanov, 2009. The use of wavelet transformation for detection and recognition of the small size contrast subsurface targets. Proceedings of the International Sierian Confewrence on Control and Communications, March 27-28, 2009, Tomsk, Russia, pp: 334-349.

CrossRefDirect Link - Sensors and Software, 1999. Practical Processing of GPR Data. Sensors and Software Inc., Brevik Place, Missisauga, Canada, Pages: 18.

Direct Link - Tomiya, M. and A. Ageishi, 2003. Registration of remote sensing image data based on wavelet transform. TSPRS Arch., 34: 145-150.

Direct Link - Xu, Y., T. Hao, Z. Li, Q. Duan and L. Zhang, 2009. Regional gravity anomaly separation using wavelet transform and spectrum analysis. J. Geophys. Eng., 62: 279-287.

CrossRefDirect Link