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 Cjk known as wavelet coefficients which are functions
of scales and positions j, k, defined as (Ni et al.,
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):
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:
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 Ai and
details Di 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.
|| Signal decomposition into Approximation Aj(t)
and detail Dj(t)
|| 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.,
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 Daubechies 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, Daubechies 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.
|| Raw radargram
||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):
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:
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.
||A window snap of the 4-level horizontal, vertical and diagonal
decompositions of the radar image
||(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.
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.