ADDIS ABABA UNIVERSITY

ADDIS INSTITUTE OF TECHNOLOGY

SCHOOL OF ELECTRICAL AND COMPUTER ENGINEERING

SUB ARRAY PROCESSING TECHNIQUES FOR DETECTION OF

SUBSURFACE BURIED UTILITIES

A Thesis Submitted to the School of Graduate Studies of Addis Ababa

University in Partial Fulfillment of the Requirements for the Degree of
Masters of Science in Communication Engineering

Advisor: Ephrem Teshale(PhD)

By: Naol Worku

Addis Ababa, Ethiopia

July, 2019
 Addis Ababa University
Addis Institute of Technology
School of Electrical and Computer Engineering

A Thesis Submitted to the School of Graduate Studies of Addis Ababa University in

Partial Fulfillment of the Requirements for the Degree of Masters of Science in
Communication Engineering

Sub Array Processing Techniques for Detection of Subsurface Buried

Utilities

By:Naol Worku

Approval by Board of Examiners

Chairman, School Graduate Signature Date

Committee

Ephrem Teshale(PhD) Signature Date

Advisor

Internal Examiner Signature Date

External Examiner Signature Date

i
 Abstract
Currently population growth and industrial expansion resulted in increased infras-
tructure spending particularly in the Addis Ababa city. The city underground has
becoming a spider’s web of utility lines, electricity, fiber optics, drainage and sanitary
sewers and water mains. From daily observation and reports new construction in addis
Ababa and a growing number of rehabilitation and replacement projects undertaken
to maintain and improve the aging infrastructure have often resulted in increased
instances of damages to underground utilities, and undesirable consequences to con-
tractors, project owners and citizens. Hence identification of buried cables, pipes,
conduits and other cylindrical utilities are very important task. Performing routine
identification of buried materials are time consuming and costly endeavor with destruc-
tive method. Before the new construction and/or replacement of aging infrastructures
performing underground utility detection survey will minimize the damage made to
the subsurface utilities. Today in Addas Ababa city, preconstruction utility detec-
tion is unknown and due to this the broken water lines and damaged telecom cables
during excavation are repeatedly reported. All this detailed information is considered
priceless for construction agencies as it saves them from unnecessary excavations and
damages. Locating the buried utility services not only avoids disruptions during the
construction but also save additional cost of getting those services installed again if
they remain undetected. Thus it is advisable to get the utility detection and locating
survey done before embark on the construction project to make it the most profitable
venture of life.
This study proposed to identify the imaging technologies that have potential for being
applied in locating underground utilities, and to analyze the conditions under which
the use of these technologies is most appropriate because not all technologies can lo-
cate all types of utilities, or be used in all types of soil or at all depths.
On this regard, a nondestructive, accurate, cost-effective, and speedy method which
will benefit the construction agencies were studied. Furthermore, feasibility and accu-
racy of using Ground Penetrating Radar (GPR) by Sub Array Processing (SAP) which
is used to estimate the location of the buried object is studied. The whole receiving
array is partitioned in small sub-arrays and for each sub-array the dominant DOA is
estimated. In particular, an efficient algorithms for the estimation of the Directions
of Arrival (DOA) of the electromagnetic waves scattered by the targets, MUSIC algo-
rithm used to estimate the dominant arrival direction of the signal. Angle of arrival at
-20◦ ,10◦ and 60◦ is used for simulation using random generated signal and MUlti SIg-
nal Classification algorithm shows sharp peak pseudo spectrum around this angle of
arrival which used to estimate the location of source signal.
The fundamentals of GPR responses from arbitrarily complex targets modelled using
software tool called GprMax which helps for deeper understanding of the operation
and detection mechanism of GPR for both 2D and 3D models. Homogeneous concrete
simulation medium with single and multiple metallic cylinder at the same and differ-
ent depth embedded in it is considered scenario. The antennas received signals with
less difference in field strength shows inverted U shape around buried cylinder with
response time difference according to the depth of the target.
Keywords: Buried object Detection, Sub Array Processing; Direction of Arrival es-
timation, MUSIC Algorithm.

ii
 Declaration
This thesis is a presentation of my original research work. Wherever contributions of
others are involved, every effort is made to indicate this clearly, with due reference
to the literature, and acknowledgement of collaborative research and discussions. The
work was done under the guidance of Dr.Ephrem Teshale at the Addis Ababa Institute
of Technology, Ethiopia.

Naol Worku
Name Signature

Place: Addis Ababa Institute of Technology, Addis Ababa,Ethiopia

Date of Submission:

This thesis has been submitted for examination with my approval as a university ad-
visor.

Ephrem Teshale(PhD)
Advisor’s Name Signature

iii
 Acknowledgement
First and foremost I always thank The Most High Holy Heavenly God for giving me all
packages of life and words can only inadequately to express my deep and sincere grat-
itude to my research Advisor Dr.Ephreme Teshale, for his meticulous care, kindness,
invaluable suggestions and generosity. His fruitful comments and insightful sugges-
tions have been a crucial formative influence on the present study. He has supported
me in every possible way since the beginning of this research. His critical and careful
reading of my writing has saved me from a lot of errors. Without his guidance and
encouragement, this research would have never come out in the present form.
Secondly, I want to thank my families who have always been the sources of encour-
agement and strength throughout my life. I am also thankful for my friends in AAiT,
for all the unforgettable moments we shared together.
I would also like to acknowledge Dire Dawa Institute of Technology, for granting
and sponsoring me for this scholarship program and finally, I must also acknowledge
Ethiopian Construction Project Management Institute (ECP M I) for raising this re-
search idea and partially funding this thesis.
Lastly I would love to sincerely thank all those who have contributed in one way or
another to this study.

iv
 Contents

1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Statement of the Problem . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.1 General . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.2 Specific . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.5 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.6 Scopes and Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.6.1 Scopes of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . 8
1.6.2 Limitations of the Thesis . . . . . . . . . . . . . . . . . . . . . . 8
1.7 Contributions of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Ground Penetrating Radar 9

2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1 History of Ground Penetrating Radar . . . . . . . . . . . . . . . 10
2.2 GPR Main Components . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.1 Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.2 Constitutive equations . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.3 Material properties . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Visualizing GPR Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.1 A-scan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.2 B-scan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3.3 C-scan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3 Sub Array Processing Technique 21

3.1 Array Signal Processing . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2 Sub Array Processing in DOA Estimation . . . . . . . . . . . . . . . . 22
3.2.1 Direction of Arrival(DOA) . . . . . . . . . . . . . . . . . . . . . 23
3.2.2 DOA Estimation Algorithms . . . . . . . . . . . . . . . . . . . . 24
3.2.3 MUSIC Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3 Factors Influencing the Efficiency of MUSIC Algorithm . . . . . . . . . 25

4 GPR Modeling and Utility Detection 27

4.1 GPR modelling Basic concepts . . . . . . . . . . . . . . . . . . . . . . . 27
4.1.1 gprMax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.1.2 Coordinate System and Conventions . . . . . . . . . . . . . . . 30
4.1.3 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . 31
4.1.4 Absorbing Boundary Conditions . . . . . . . . . . . . . . . . . . 32
4.2 Modelling of single buried metallic object . . . . . . . . . . . . . . . . 34

v
 4.3 Modeling of of multiple buried metallic objecs at the same depth . . . . 35
4.3.1 Multiple buried metallic objects at different depth and size . . . 35
4.3.2 3D Model of very complex environment . . . . . . . . . . . . . 36

5 Result and Discussion 41

5.1 Detection of buried metallic object . . . . . . . . . . . . . . . . . . . . 41
5.1.1 Result I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.1.2 Result II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.1.3 Result III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.1.4 Horizontal Location Estimation . . . . . . . . . . . . . . . . . . 42
5.1.5 Result IV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

6 Conclusion and Recommendations 48

6.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.2 Recommendation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.3 Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

vi
 List of Figures

2.1 The main components associated with a ground penetrating radar system 14
2.2 Configuration and representation of an A-scan . . . . . . . . . . . . . . 18
2.3 (a) Multiple A-scans forming a B-scan (b) Representation of a B-scan
on a grey-scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4 Configuration and representation of C-scan . . . . . . . . . . . . . . . . 20

3.1 Example of geometrical configuration SAP . . . . . . . . . . . . . . . . 22

3.2 Angle of arrival and the sharp pseudo spectrum around angle of arrival 25

4.1 Single FDTD Yee cell showing electric (red) and magnetic (green) field
components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.2 Single FDTD Yee cell showing electric (red), magnetic (green), and
zeroed out (grey) field components for 2D transverse magnetic (TM)
z-direction mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.3 gprMax coordinate system and conventions. . . . . . . . . . . . . . . . 31
4.4 GPR forward problem showing computational domain bounded by Ab-
sorbing Boundary Conditions (ABCs) . . . . . . . . . . . . . . . . . . . 33
4.5 gprMax2D Model of first scenario . . . . . . . . . . . . . . . . . . . . . 35
4.6 Three cylinder buried in concrete at the same depth . . . . . . . . . . . 36
4.7 Different size cylinder at different depth . . . . . . . . . . . . . . . . . 36
4.8 3D First Model of heterogeneous environment . . . . . . . . . . . . . . 38
4.9 3D Model of very complex environment . . . . . . . . . . . . . . . . . . 40

5.1 B-scan of model of a single metallic cylinder buried in concrete . . . . 42

5.2 B-scan of model of a three metallic cylinder buried in concrete at the
same depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.3 B-scan of model of a three metallic cylinder buried in concrete at dif-
ferent depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.4 B-scan for model of metal cylinder buried in heterogeneous medium . . 46

vii
 List of Tables

2.1 Table shows antenna frequency, approximate depth penetration and ap-
propriate application . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

4.1 gprMax2D Model list of used materials . . . . . . . . . . . . . . . . . . 34

4.2 3D Heterogeneous model list of used materials . . . . . . . . . . . . . . 37
4.3 3D Very complex medium model list of used materials . . . . . . . . . . 39

6.1 2D/3D GPR Modelling commands of gprMax . . . . . . . . . . . . . . 56

viii
 List of Acronyms

ABC Absorbing Boundary Conditions

AET Acoustic Emission Testing
ACFM Alternating Current Field Measurement
AoA Angle of Arrival
CFL Courant, Freidrichs and Lewy
CFS Complex frequency shifted
CUDA Compute Unified Device Architecture
DOA Direction of Arrival
EM Electromagnetic
ET Electromagnetic Testing
FDOA Frequency Difference of Arrival
FDTD Finite-Difference Time-Domain
GPR Ground Penetrating Radar
GPU Graphics processing unit
MUSIC MUlti SIgnal Classification
NDT Non Destructive Method
PML Perfectly Matched Layer
SAP Sub Array Processing
SNR Signal-to-Noise Ratio
TDOA Time Difference of Arrival
TEM Transverse Electromagnetic
TM Transverse Magnetic

ix
 Chapter 1

Introduction

1.1 Background
Currently population growth and industrial expansion resulted in increased infras-
tructure spending particularly in the Ethiopia. The urban underground has becoming
a spider’s web of utility lines, telecom cables, electricity,fiber optics, drainage and
sanitary sewers and water mains. New construction in urban areas and a growing
number of rehabilitation and replacement projects undertaken to maintain and im-
prove the aging infrastructure have often resulted in increased instances of damages
to underground utilities, and undesirable consequences to contractors, project owners
and citizens. These consequences include construction delays, design changes, claims,
property damages, service breakdowns, disruption of neighboring business.
This thesis study to identify, through literature review and case studies, the imaging
technologies that have potential for being applied in locating underground utilities,
and to analyze the conditions under which the use of these technologies is most ap-
propriate because not all technologies can locate all types of utilities, or be used in all
types of soil or at all depths. Hence identification of buried cables, pipes, conduits and
other cylindrical utilities are a very important task. Performing routine identification
of buried materials are time consuming and costly endeavor with destructive method.
Therefore, a nondestructive, accurate, cost-effective, and speedy method can benefit
the construction agencies. In this thesis, feasibility and accuracy of using Ground
Penetrating Radar (GPR) by Sub Array Processing (SAP) are studied through a sys-
tematic research methodology.
Ground Penetrating Radar (GPR) illuminates the ground through a set of antennas
and the collected echo is analyzed, in order to extract information about the scenario
and to localize the sought objects. The use of GPR for NDT is just one of the many
different areas, where radar is being applied as an important tool assisting the engi-
neers in their efforts to determine the presence or absence as well as the kind of key
underground features. The main advantage s of GPR are: its fast data acquisition ca-
pability, its high resolving ability and the fact that it responds equally well to metallic
and non-metallic targets.
In applications involving smart antennas, and in the presence of several transmitters
operating simultaneously, it is important for a receiving array to be able to estimate
the angles of arrival, in order to decode how many emitters are present and what are
their possible locations.In the scenario considered in this paper, the sources are the
currents induced on buried objects and a set of electric-field values, measured by a
linear array of receiving antennas parallel to the air–soil interface, is used to estimate
the directions of arrival of the electric field back-scattered by the targets.

1
Directions of arrival (DOA) algorithms assume that the sources are in the far-field
region of the receiving array, so that the received wave front can be considered as
planar and the main angular direction of the field can be estimated. However, in elec-
tromagnetic sensing of subsurface utilities, the scatterers are usually quite near to the
antennas. In order to be able to work in near-field conditions, a sub-array processing
(SAP) approach is adopted in this paper: the whole receiving array is partitioned in
small sub-arrays and for each sub-array the dominant direction of arrivals is estimated.
In this way, it is possible to successfully locate also objects that are in the far-field of
the sub-array, although being in the near-field of the whole array.
Successful modelling attempts of ground penetrating radar have been reported by
many authors. Most of the proposed approaches are based on the finite-difference
time-domain (FDTD) method. The main reasons for such widespread use of the FDTD
method are: its ease of implementation in a computer program at least at a simple
introductory level and its good scalability when compared with other popular electro-
magnetic modelling methods such the finite- element and integral techniques. This
thesis work is comprised of six chapters. The first chapter deals with an introduction
of the thesis. It introduces problems of the statement, thesis objectives, methodolo-
gies used, scopes, and limitations of the study, second chapter is about overview of
GPR and historical background GPR, third chapter is about Sub Array Processing
Techniques, Direction of Arrival Estimation and MUSIC Algorithm, Fourth chapter
is all about GPR modelling and Utility detection basic concepts with different mod-
els. Chapter five is about Result and Discussion and the finally sixth chapter will be
Conclusion and Recommendation.

2
1.2 Statement of the Problem
Identification of buried cables, pipes, conduits are some of an important consideration
for construction quality assurance of new different constructions and structural capac-
ity estimation of existing constructions. Currently, Ethiopian construction agencies
performs identification of buried objects using traditional methods which is labor in-
tensive, cost ineffective and destructive.
No subsurface utility maps for preconstruction process and no label indicators spe-
cially for water line which result in high number of broken water line and flow of
water, sewage water linkage due to heavy load vehicles, destruction of telecom cable
during preconstruction process and lost revenue and high maintenance cost of damage
made are some of the problems and including the following problems. Utility dam-
ages during construction are very significant and on the rise, resulting in construction
delays, design changes, claims,property damages, service breakdowns, disruption of
neighboring businesses and even injuries and lost lives.There is no governmental and
non governmental organization which facilitate the subsurface utility related issues,no
permanent publicly-available site records (drawings) are furnished and utility depths
are often not provided. Given the drawbacks of mentioned above, it is important to
find a nondestructive and cost effective alternative system to maintain safety, service-
ability, and durability of subsurface utilities.
GPR provides continuous information in a rapid and cost effective manner, and most
importantly can save time, money and provide the department with the possibility of
continuous quality control on subsurface utilities.
Estimation of direction of arrival of scattered wave from targets, finding number of
sources and signal waveforms are most difficult scenario. DOAs algorithm [5] is promis-
ing algorithm with SAP techniques for signal processing.

1.3 Objectives
1.3.1 General
The main objective of this thesis is to apply GPR system to investigate subsurface
buried objects using SAP technique.

1.3.2 Specific
The specific objectives centered on data processing flow and the interpretation of the
images. This study is carried out specifically to accomplish the followings:

• Study DOA algorithms and DOA estimation.

• Estimate location of burials.

• Detection of single buried objects like cable, metals pipe.

3
1.4 Methodology
In this study, related secondary sources of data on GPR, SAP techniques, DOA Al-
gorithms with the help of standardization documents, previous researcher’s studies,
IEEE articles, development in subsurface utility technology, general concept and ba-
sics of GPR are used. The methodology of this study have been used to explore the
overall importance of SAP techniques with DOA algorithm to get high resolution of
GPR image based on the inputs retrieved from review of related problem in mind, the
methodology of this study is described as follows:

GPR Modelling:

• Identification of input parameter: (input file commands)

– Essential - required to run any model, such as the domain size and spatial
discretization.
– General - provide further control over the model.
– Material - used to introduce different materials into the model.
– Object construction - used to build geometric shapes with different consti-
tutive parameters.
– Source and output - used to place source and output points in the model.
– PML - provide advanced customisation and optimisation of the absorbing
boundary conditions.
– Identification of visualization tool kits: (Paraview tool )

• Identification of image processing techniques with DOA estimation algorithms:SAP

techniques, MUSIC algorithm

• Decision of Detection /Scan type:A-scan, B-scan, C-scan

Simulations:

• Identification of most appropriate simulation tools:(gprMax, MATLAB )

• Simulate using all the above input parameters with single/multiple cylindrical
metal embedded in it.

4
1.5 Literature Review
Different studies have been made in the area of Subsurface Utility Detection, among
others the following are the ones reviewed to frame the current study in line with the
problem on hand.
Acoustic Emission Testing is considered quite unique among the non-destructive test-
ing methods. Acoustic Emission (AE) refers to the generation of transient elastic waves
produced by a sudden redistribution of stress in a material [1]. When a structure is
subjected to an external stimulus (change in pressure, load, or temperature), localized
sources trigger the release of energy, in the form of stress waves, which propagate to
the surface and are recorded by sensors. With the right equipment and setup, motions
on the order of pico meters (10 -12 pm) can be identified. Sources of AE vary from
natural events like earthquakes and rock bursts to the initiation and growth of cracks,
slip and dislocation movements, melting, twinning, and phase transformations in met-
als. In composites, matrix cracking and fiber breakage and debonding contribute to
acoustic emissions. AE’s have also been measured and recorded in polymers, wood,
and concrete, among other materials. Detection and analysis of AE signals can sup-
ply valuable information regarding the origin and importance of a discontinuity in a
material. Because of the versatility of Acoustic Emission Testing (AET), it has many
industrial applications (e.g. assessing structural integrity, detecting flaws, testing for
leaks, or monitoring weld quality) and is used extensively as a research tool. Acoustic
Emission is unlike most other nondestructive testing (NDT) techniques in two regards.
The first difference pertains to the origin of the signal. Instead of supplying energy
to the object under examination, AET simply listens for the energy released by the
object. AE tests are often performed on structures while in operation, as this provides
adequate loading for propagating defects and triggering acoustic emissions. The sec-
ond difference is that AET deals with dynamic processes, or changes, in a material.
This is particularly meaningful because only active features (e.g. crack growth) are
highlighted. The ability to discern between developing and stagnant defects is signifi-
cant. However, it is possible for flaws to go undetected altogether if the loading is not
high enough to cause an acoustic event. Furthermore, AE testing usually provides an
immediate indication relating to the strength or risk of failure of a component. Other
advantages of AET include fast and complete volumetric inspection using multiple
sensors, permanent sensor mounting for process control, and no need to disassemble
and clean a specimen. Unfortunately, AE systems can only qualitatively gauge how
much damage is contained in a structure. In order to obtain quantitative results about
size, depth, and overall acceptability of a part, other NDT methods (often ultrasonic
testing) are necessary. Another drawback of AE stems from loud service environments
which contribute extraneous noise to the signals. For successful applications, signal
discrimination and noise reduction are crucial. A reliable analysis of of acoustic emis-
sion signals and the interpretation of the data in material testing are usually only
possible in case where the signals have been localized successfully.
Electromagnetic testing is a general test category that includes Eddy Current testing,
Alternating Current Field Measurement (ACFM), and Remote Field testing magnetic
particle testing is also an electromagnetic test, due to its widespread use it is consid-
ered a stand-alone test method rather as than an electromagnetic testing technique.
All of these techniques use the induction of an electric current or magnetic field into
a conductive part, then the resulting effects are recorded and evaluated [2]. Electro-
5
magnetic Testing (ET), as a form of nondestructive testing, is the process of inducing
electric currents and magnetic fields inside a test object and observing the electro-
magnetic response.A defect inside a test object creates a measurable response that
differs from background noise and allows us to detect and characterize surface and
sub-surface flaws in conductive materials.
According to [3] using a thermal model and measured IR images to detect the presence
of buried objects and characterize them in terms of thermal and geometrical proper-
ties. The inverse problem is mathematically stated as an optimization one using the
well-known least-square approach. The main difficulty in solving this problem comes
from the fact that it is severely ill posed due to lack of information in measured data.
And NDT methods such as Laser Testing Methods,Thermal/Infrared Testing,Ultrasonic
Testing and others have been developed in [4].
There are still few published works dealing with the automatic detection of patterns
associated with buried objects. The classical Hough transform is used in order to
identify linear segments in the image, representing transitions between layers of differ-
ent electrical impedances. The authors in [5] proposed also a method for extracting
hyperbolic signatures of buried there forestages: 1) preprocessing step to reduce noise
and undesired system effects; 2) image segmentation with an artificial neural network
classifier to identify areas potentially containing object reflections; and 3) Hough trans-
form to detect hyperbolic patterns.
In [6] some preprocessing steps aiming at enhancing the signature of buried targets
are implemented. Then, automatic image interpretation is carried out by a detector
based on artificial neural networks.
The authors in [7] applied a fuzzy clustering approach to identify hyperbolas from
GPR images beforehand de-noised. In this paper, they proposed a novel system to
identify and classify buried objects from GPR imagery. The entire process is subdi-
vided into four steps. In the first one, a preprocessing procedure is implemented for: 1)
reducing noise; 2) eliminating the undesired presence of the ground surface echo; and
3) compensating propagation losses. Noise reduction is performed with a median filter,
while the elimination of the ground surface echo is done by a simple average operation.
A time-gain filter is used to compensate signal amplitude for losses due to spreading
and attenuation. In the second step, the resulting preprocessed image is thresholded
to discriminate between objects and background. This binarization operation, which
allows to put under light the parts of the imagecontaining potential targets, is based
on the fact that buried objects are generally associated with large amplitude echoes.
It is implemented by means of the Kapur’s thresholding technique, which relies on the
entropy maximization principle [8]. The third step of the system consists to identify
the targets in the obtained binary image in a completely unsupervised way. This is
done by means of a search of linear and hyperbolic patterns representing potential
targets.
Several techniques have been proposed for the localization of targets from the elec-
tromagnetic field measured above the ground, ranging from synthetic aperture tech-
niques, inverse-scattering-based algorithms and synthetic aperture focusing technique
[9, 10].In synthetic aperture approaches, an image representing the sub-surface region
is achieved by refocusing the backscattered data with a coherent summation of the
measured scattered field after proper phase compensation.
On the contrary, inverse-scattering-based methods usually aim at directly retrieving
the full distribution of the dielectric properties (e.g., the relative dielectric permittivity

6
and the electric conductivity) of a predefined investigation area starting from the mea-
sured electromagnetic field. Such task is pursued by solving an inverse problem, which
turns out to be nonlinear and ill-posed. Moreover, in order to tackle the nonlinearity
of the involved equations, linearization techniques have also been proposed [11, 12] .
In order to tackle the nonlinearity of the involved equations, linearization techniques
have also been proposed (e.g., those using the Born approximation) [13] Innovative
techniques must still be developed in order to mitigate the drawbacks of existing ap-
proaches, especially when real-time operations are needed. In this paper a statistical
method for the localization of objects will be employed. It is based on the use of smart
antenna sub-array processing (SAP) techniques. The presence of several antennas op-
erating simultaneously is important for a receiving array to be able to estimate the
direction of arrival, in order to locate the area to scan and what are possible locations
of targets.The method of this thesis work is to merge the concept of sub array process-
ing to estimate where the target located with the help direction of arrival estimating
algorithms.Once the direction of the incoming signal from the target is estimated, pos-
sible to know the area of target to scan using GPR system which minimize time of
scanning unwanted area.

7
1.6 Scopes and Limitations
1.6.1 Scopes of the Thesis
The subsurface utility detection in this study is done by using only simulation tools.Thus,
the actual on sites conditions, field tests and related issues are not taken into consid-
eration. Hence the outcomes of the study is basically based simulation on results.

1.6.2 Limitations of the Thesis

The practical study of thesis might have so many limitations and most of all are the
followings

• Data acquisition may be slow over very difficult terrain.

• Depth of penetration is limited in materials with high electrical conductivities

like clays.

• Energy may be reflected and recorded from above-ground features,like walls,

canopies,etc., unless antennae are well shielded.

• Artifacts in the near surface (reinforcing bars, boulders, components of made

ground etc.) may scatter the transmitted energy and complicate the received
signal and/or reduce depth of penetration.

• The assumption made is homogeneous and heterogeneous environments.

1.7 Contributions of the Thesis

The main contribution of this thesis work is to study sub surface utility, how to iden-
tify, estimate location so that the detection of subsurface utility can be proposed for
implementation.
Since the thesis work is initiated as a result of interests form a sponsor Institute to
implement efficient subsurface utility detection, the results of this work can be used
as an input to the intended institute (Ethiopian Construction Project Management
Institute).It will help them through their way in achieving their goal of implementing
fast and reliable subsurface utility detection system.

8
 Chapter 2

Ground Penetrating Radar

2.1 Overview
Ground penetrating radar (GPR) is a relatively new geophysical technique. The last
decade has seen major advances and there is an overall sense of the technology reaching
a level of maturity. The history of GPR is intertwined with the diverse applications
of the technique. GPR has the most extensive set of applications of any geophysical
technique. As a result, the spatial scales of applications and the diversity of instru-
ment configurations are extensive. Both the value and the limitations of the method
are better understood in the global user community.
This chapter is to provide a brief history of the method, a discussion of current trends
and give a sense of future developments. Before delving into the history, GPR needs
definition. GPR uses electromagnetic fields to probe lossy dielectric materials to detect
structures and changes in material properties within the materials.Most applications
to date have been in natural geologic materials, but widespread use in man-made
composites such as concrete, asphalt and other construction materials also occurs. In
such lossy dielectric materials, electromagnetic fields will penetrate to some depth be-
fore being absorbed. With GPR, the electromagnetic fields propagate as essentially
non-dispersive waves. The signal emitted travels through the material, is scattered
and/or reflected by changes in impedance giving rise to events similar to the emitted
signal. In other words, signal recognition is simple because the return signal looks
like the emitted signal. One has to contrast GPR measurements with electromagnetic
induction sounding methods where the fields are diffusive and dispersive in character.
GPR field behavior occurs over a finite frequency range generally referred to as the
GPR plateau where velocity and attenuation are frequency independent. The GPR
plateau usually occurs in the 1 MHz to 1000 MHz frequency range. At lower frequencies
the fields become diffusive in character and pulses are dispersed. At higher frequencies
several factors increase signal absorption such that penetration is extremely limited.
The foundations of GPR lie in electromagnetic (EM) theory. The history of this field
spans more than two centuries and is the subject of numerous texts.
This overview outlines the basic building blocks needed to work quantitatively with
GPR. Maxwell’s equations mathematically describe the physics of EM fields, while
constitutive relationships quantify material properties. Combining the two provides
the foundations for quantitatively describing GPR signals. GPR uses electromagnetic
wave propagation to image and identify changes in electrical and magnetic properties
in the ground. GPR systems are used to locate underground utility lines, reinforcing
and post tensioning in concrete, monitor airplane runways for structural integrity, con-
ducting groundwater studies, detect unexploded land mines as well as forensic research

9
and surveying land for Construction Purposes.
The tools and methods of locating buried utilities are quite diverse. The most com-
mon approach is energizing metal pipes and cables with electric currents and using a
magnetic field sensor to detect the current. Provided the target object can be exposed
for connection or current can be induced, sufficient current remains on the object and
the magnetic field at the detector is strong enough, then this technique works well
and is very cost effective. When access is difficult, the electrical current does not flow
(i.e. non-metallic element or broken connection) or external noise makes detection
impossible, GPR provides an alternative. GPR provides its own source of energy,
locates both metallic and non-metallic objects, detects disturbed soil conditions and
other buried structures. Other direct approaches are to trench, hand dig, or vacuum
excavate to expose features. A priori knowledge and accurate as-built drawings are
needed to be effective with these techniques. Generally, these are not available or suf-
ficiently accurate. GPR can also be used to Locate Underground Utilities. Traditional
electromagnetic induction utility locating equipment require utilities to be conductive.
They are ineffective for locating plastic conduits or concrete storm and sanitary sew-
ers. Because GPR can detect variations in dielectric properties in the subsurface, it is
highly effective for locating non-conductive utilities.
Resolution and penetration: The depth range of GPR is limited by the electrical
conductivity of the ground, the transmitted center frequency and the radiated power.
As conductivity increases, the penetration depth decreases. Lower frequencies can
reach depths up to tens of meters (e.g. 100 MHz can travel up to 20 meters) with
a resolution of tens of centimeters, while higher frequencies can give a resolution of
centimeters but up to depths of several meters

2.1.1 History of Ground Penetrating Radar

The following is necessarily brief and intended to give high lights. References lead to
other perspectives for those interested in a more extensive understanding of GPR[14].
It is interesting to note that accounts of some activities are published many years later
and sometimes not all. During this time a great deal of research on radio wave prop-
agation above and along the surface of the earth occurred. Although several hints at
the possibility of using radio waves to probe the subsurface are mentioned, there are
no reports of successfully making this type of measurement. Vast number of papers
appeared on the subject of communications, direction finding and radar. In this time
frame, the first reported attempt at measuring subsurface features with radio wave
signals was reported. El Said (1956) attempted to use the interference between direct
air transmitted signals and signals reflected from the water table to image the water
table depth. The next reported observation of radio frequency sounding of geological
materials came about when the USAF reported altimeter errors when attempting to
land aircraft on the Greenland ice sheet (Waite and Schmidt (1961)). This was the
first time that repeatable indications of penetration into the subsurface through a nat-
urally occurring material were reported. This spawned the era of researchers focused
on developing radio echo sounding in ice. The majority of activity during this interval
involved the radio echo sounding in ice. Groups, such as the Scott Polar Research
Institute at Cambridge, Bailey et al (1964) and the Geophysical and Polar Research
Center at the University of Wisconsin, Bentley (1964), Walford (1964) were active in
polar regions and also on glaciers. During this time the ice radio echo sounding activ-

10
ity continued. In addition, applications in other favorable geologic materials started
to be explored. Cook (1973) explored the use in coal mines since coal can be a low
loss dielectric material in some instances.
Similarly, Holser et al (1972), Unterberger (1978) and Thierbach (1973) initiated eval-
uations in underground salt deposits for similar reasons During this same period Morey
and others formed Geophysical Survey Systems Inc. which has been manufacturing
and selling ground penetrating radar since that time (Morey (1974)). In addition a
better understanding of electrical properties of geologic materials at radio frequencies
started to become available. Work such as that by Olhoeft (1975) led to a much bet-
ter understanding of the electrical character of natural occurring geological materials
and the relationship between electrical conductivity and dielectric polarization of these
materials, applications started to grow because of the availability of technology and a
better understanding of geology. The Geological Survey of Canada explored a number
of applications, the primary one being a better understanding of permafrost terrain
in the Canadian Arctic. Proposals for pipelines out of the Arctic to carry oil and gas
to southern markets drove a great deal of interest in engineering in frozen soil and
environments. GPR was a tool which offered great promise and some of the initial
results are reported by Annan and Davis (1976).
Experiments with GPR were reported by the Stanford Research Institute where mea-
surements were made by Dolphin et al (1978) for archeological applications. Other
work carried out in this period which paralleled the Geological Survey of Canada
permafrost efforts was lead by Olhoeft at the United States Geological Survey who
worked on the Alaska pipeline routes. Extensive work was carried out in potash mines
in western Canada. This led to a whole series of ever improving GPR measurements
and work in this geological setting by the Geological Survey of Canada. These results
were reported by Annan et al (1988). Further coal mine developments were reported
by Coon et al (1981).
In addition, the potential for use of borehole radar to investigate rock quality in poten-
tial hard rock nuclear waste disposal sites became a topic of interest. The Geological
Survey of Canada and Atomic Energy of Canada supported this work (Davis and An-
nan (1986)). Commercial units were used for most of this work and the number of
activities spawned commercial interest. Geophysical Survey Systems Inc. remained
the dominant supplier at this time but Ensco/Xadar was spawned in an attempt to
create an alternate commercial product. One of the major issues that was noted by the
Geological Survey of Canada work was the great difficulty in using existing equipment
in remote areas. Equipment was heavy, bulky and consumed too much power. In
addition there was a need to get data into a digital form to exploit the digital seismic
processing advances which were rapidly evolving in the petroleum seismic field at the
time. interest in GPR waned to a degree. The initial optimism for the technology
gave way to the reality that many environments weren’t favorable for GPR. There was
considerable confusion over whether failures were equipment related or due to natural
material responses. In addition, very little money for technology development was
available. During this time OYO Corporation of Japan developed a radar product
called “Georadar” spawned by association with Xadar developments. This instrument
met some initial commercial success in Europe.
A-Cubed Inc. was formed in 1981 and started development of ground penetrating
radars. The low frequency digital GPR developments were reported by Davis et al
(1985). This technology development led to the pulseEKKO series of GPR’s. The nu-

11
clear waste disposal problem was continually studied and a number of countries funded
the Swedish Geological Survey in the development of borehole radar. This work is re-
ported by Olsson et al (1987). Other applications for GPR such as road investigations
and utility mapping met with mixed success. In general, the technology was quite
new and not optimized for these applications. Work by Ulriksen (1982) provided a
good foundation for some of these applications. Many non-commercial developments
occurred with prototypes that embodied the ideas for portability, digital recording and
the use of fiber optics cables. Other little reported work was conducted by Southwest
Research and the U.S. Army on borehole GPR to detect tunneling in sensitive military
areas (Owen (1981)). Interest in GPR waned to a degree. The initial optimism for
the technology gave way to the reality that many environments were not favorable for
GPR. There was considerable confusion over whether failures were equipment related
or due to natural material responses. In addition, very little money for technology
development was available. During this time OYO Corporation of Japan developed
a radar product called “Georadar” spawned by association with Xadar developments.
This instrument met some initial commercial success in Europe.
A-Cubed Inc. was formed in 1981 and started development of ground penetrating
radars. The low frequency digital GPR developments were reported by Davis et al
(1985). This technology development led to the pulseEKKO series of GPR’s. The
nuclear waste disposal problem was continually studied and a number of countries
funded the Swedish Geological Survey in the development of borehole radar. This
work is reported by Olsson et al (1987).
Other applications for GPR such as road investigations and utility mapping met with
mixed success. In general, the technology was quite new and not optimized for these
applications. Work by Ulriksen (1982) provided a good foundation for some of these
applications. Many non-commercial developments occurred with prototypes that em-
bodied the ideas for portability, digital recording and the use of fiber optics cables.
Other little reported work was conducted by Southwest Research and the U.S. Army
on borehole GPR to detect tunneling in sensitive military areas (Owen (1981)).
The real explosion in the advancement of GPR occurred during this period. Many
groups worldwide became interested in the technology. On the commercial side, Geo-
physical Survey Systems Inc. exhibited strong commercial success and was bought
by OYO Corporation. During this period, Mala Geosciences was spawned from the
Swedish Geological Survey roots. ERA in the UK also became more active using its
research into unexploded ordinance and landmine detection to create commercial prod-
ucts. Sensors and Software Inc. grew rapidly broadening its pulse EKKO product line.
On the research side, much attention started to be paid by both the geophysical and
electrical engineering community. Developments such as multi-fold data acquisition
(Fisher et al (1992)), digital data processing (Maijala (1992), Gerlitz et al (1993)),
and 2D numerical simulation (Zeng et al (1995), Cai and McMechan (1995)) occurred.
Advances in applications in archeology (Goodman (1994)), environmental (Brewster
and Annan (1994)), geological stratigraphy using radar facies (Jol (1996)) and many
other areas expanded. Environmental borehole GPR development was reported by
Redman et al (1996).
Ground penetrating radar user meetings became more formalized and were held every
2 years at various locations around the world. This meeting provided a forum for
the leading players in this field to meet, present results and discuss problems. These
meetings led to series of proceeding publications which are listed as references. These

12
proceedings provide a great deal of information for new users to the GPR field. In
early 20th century’s the evolution of the computers drove all of GPR advances. Nu-
merical modelling of full 3D problems became possible albeit still with large computers
(Holliger and Bergmann (2000), Lampe and Holliger (2000)). The ability to manage
the large volumes of information in digital form and manipulate them quickly became
routine. As a result, acquisition of data on grids to make maps and grids and 3D visu-
alization became practical (Grasmueck (1996), Annan et al (1997)). The commercial
market and demand resulted in a variety of different and simpler systems such as the
Noggin from Sensors and Software Inc.
Very recent history; In January 2018, US RadarIn the leading innovator of ground
penetrating radar (GPR) systems,introduced the GP Rover, which combines US Radar’s
advanced triple-bandwidth ground-penetrating radar technology with precision global
GPS connectivity to create subsurface maps in real time, with no base station or con-
trol points required and accurate within two in. Maps created using the GP Rover
provide a permanent and updateable record that can be shared by all contributors to a
project,The GP Rover overcomes the limitations of impermanent ’mark-outs’ that of-
ten lead to costly errors in subsequent excavations. The newly patented tilt-correction
technology on the GP Rover automatically compensates for slope and cross-slope of
the GPR and GPS antennas to increase overall accuracy on uneven terrain. Field
data collected by the GP Rover may be saved and shared in any of the five most
common file formats – SHP, DXF, KMZ, CSV and RTF are used in popular mapping
and engineering applications such as AutoCad, Google Earth, GIS and point files.
The information and maps the GP Rover generates will be valuable to users across
a variety of industries.With its combination of accuracy and ease of use for locating
and mapping, the GP Rover is poised to become the new measure of excellence in
ground-penetrating radar systems. Among the potential users are contractors, engi-
neers, utility-locating service providers, public works organizations, utility companies,
municipalities, planners, architects, surveyors, property owners, oil fields, ports and
inspectors. The GP Rover combines high accuracy and differential GPR locating ca-
pabilities with the most accurate and responsive GPS location-positioning software in
an easy-to-learn, easy-to-use operator interface. The GP Rover does not require com-
plicated setup and network configuration. It automatically connects to the worldwide
GPS network and configures as soon as the radar is turned on. With the GP Rover, any
operator can create an accurate point file, preserve location points without destructive
surface markings, display location points on the integrated computer screen, export
points to mapping applications and view in real-time locations on as-built maps.

2.2 GPR Main Components

A GPR system is made up of three main components: Control unit,Antenna, Power
Supply. The control unit contains the electronics which trigger the pulse of radar en-
ergy that the antenna sends into the ground. It also has a built-in computer and hard
disk/solid state memory to store data for examination after fieldwork. Some systems,
such as the GSSI SIR 30 are controlled by an attached Windows laptop computer with
pre-loaded control software.
This system allows data processing and interpretation without having to download
radar files into another computer. The antenna receives the electrical pulse produced

13
Figure 2.1: The main components associated with a ground penetrating radar system
[15]

by the control unit, amplifies it and transmits it into the ground or other medium at a
particular frequency. Antenna frequency is one major factor in depth penetration.The
higher the frequency of the antenna, the shallower into the ground it will penetrate.
A higher frequency antenna will also ‘see’ smaller targets. Antenna choice is one of
the most important factors in survey design. The following table which shows antenna
frequency, approximate depth penetration and appropriate application
The GPR antenna, houses the transmitter and receiver; and a profiling recorder, which
processes the received signal and produces a graphic display of the data. The trans-
mitter radiates repetitive short-duration EM signals into the earth from an antenna
moving across the ground surface. Electromagnetic waves are reflected back to the
receiver by interfaces between materials with differing dielectric constants. The inten-
sity of the reflected signal is a function of the contrast in the dielectric constant at
the interface, the conductivity of the material, which the wave is traveling through,
and the frequency of the signal. Subsurface features which may cause such reflections
are: 1) natural geologic conditions such as changes in sediment composition, bedding
and cementation horizons, voids, and water content; or 2) man-introduced materials
or changes to the subsurface such as soil backfill, buried debris, tanks, pipelines, and
utilities. The profiling recorder receives the signal from the antennae and produces
a continuous cross section of the subsurface interface reflections, referred to as reflec-
tors. Depth of investigation of the GPR signal is highly site specific, and is limited
by signal attenuation (absorption) of the subsurface materials. Signal attenuation is
dependent upon the electrical conductivity of the subsurface materials. Signal attenu-
ation is greatest in materials with relatively high electrical conductivity such as clays
and brackish groundwater, and lowest in relatively low conductivity materials such
as unsaturated sand or rock. Maximum depth of investigation is also dependent on
antenna frequency and generally increases with decreasing frequency; however, the
14
 Table 2.1: Table shows antenna frequency, approximate depth penetration and appro-
priate application
Appropriate
Application Pr.Antenna Choice Sec.Antenna Appr.Depth Range
Structural Concrete,
Roadways, Bridge Decks 2600 MHz 1600 MHz 0-0.3 m (0-1.0 ft)
Structural Concrete,
Roadways, Bridge Decks 1600 MHz 1000 MHz 0-0.45 m (0-1.5 ft)
Structural Concrete,
Roadways, Bridge Decks 1000 MHz 900 MHz 0-0.6 m (0-2.0 ft)
Concrete, Shallow Soils,
Archaeology 900 MHz 400 MHz 0-1 m (0-3 ft)
Shallow Geology, Utilities,
USTs, Archaeology 400 MHz 270 MHz 0-4 m (0-12 ft)
Geology, Environmental,
Utility, Archaeology 270 MHz 200 MHz 0-5.5 m (0-18 ft)
Geology, Environmental,
Utility, Archaeology 200 MHz 100 MHz 0-9 m (0-30 ft)
Geologic Profiling 100 MHz MLF (16-80 MHz) 0-30 m (0-90 ft)
Geologic Profiling MLF (16-80 MHz) None >30 m (90 ft)

ability to identify smaller features is diminished as frequency decreases. The various

GPR antennas used are internally shielded from above ground interference sources.
Accordingly, the GPR signal is minimally affected by nearby above ground conductive
objects such as metal fences, overhead power lines, and vehicles.
A GPR survey is performed by towing an antenna across the ground along predeter-
mined transect lines. The antenna is either pulled by a person or towed behind a
vehicle. Preliminary GPR transects are performed over random areas of the site to
calibrate the GPR equipment and characterize overall site conditions. The optimum
time range settings are selected to provide the best combination of depth of investiga-
tion and data resolution for the subsurface conditions at the site. Ideally, the survey is
performed along a preselected system of perpendicular or parallel transect lines. The
configuration of the transect lines is designed based on the geometry and size of the
target and the dimensions of the site. The beginning and ending points of the transect
lines and grid intersection points, or nodes, are marked on the ground with spray paint
or survey flags.
A grid system is used to increase the probability of crossing the short axis of a target
providing a more definitive signature in the data. The location of the antenna along
a transect line is electronically marked on the cross section at each grid intersection
point to allow correlation of the data to actual ground locations. The location of the
targets can be marked on the ground surface using spray paint or survey flags.

15
2.2.1 Maxwell’s equations
In mathematical terms, EM fields and relationships are expressed as follows:

∂B
∇×E = (2.1)
∂t
∂D
∇×H =J + (2.2)
∂t
∇·E =q (2.3)
∇·B =0 (2.4)
where E is the electric field strength vector (V /m); q is the electric charge density
(C/m3 ); B is the magnetic flux density vector (T); J is the electric current density
vector (A/m2 ); D is the electric displacement vector (C/m2 ); t is time (s); and H
is the magnetic field intensity (A/m). Maxwell succinctly summarized the work of
numerous researchers in this compact form. From these relationships, all classic EM s
(induction, radio waves, resistivity, circuit theory, etc.) can be derived when combined
with formalism to characterize material electrical properties [16, 17, 18].

2.2.2 Constitutive equations

Constitutive relationships are the means of describing a material’s response to EM
fields. For GPR, the electrical and magnetic properties are of importance. Constitutive
equations ((2.4 − 2.6)) provide a macroscopic (or average behavior) description of how
electrons, atoms, and molecules respond enmasse to the application of an EM field.

J = σE (2.5)

D = εE (2.6)
B = µH (2.7)
Electrical conductivity σ characterizes free charge movement (creating electric cur-
rent) when an electric field is present. Resistance to charge flow leads to energy
dissipation.Dielectric permittivity ε characterizes displacement of charge constrained
in a material structure to the presence of an electric field. Charge displacement results
in energy storage in the material. Magnetic permeability µ describes how intrinsic
atomic and molecular magnetic moments respond to a magnetic field. For simple ma-
terials, distorting intrinsic magnetic moments store energy in the material.
σ, ε and µ are tensor quantities and can also be nonlinear (i.e. σ = σ(E)). For virtu-
ally all practical GPR issues, these quantities are treated as field-independent scalar
qualities. (In other words, the response is in the same direction as the exciting field
and is independent of field strength.) Although these assumptions are seldom fully
valid, to date, investigators working on practical applications have seldom been able
to discern such complexity. Material properties can also depend on the history of the
incident field. Time history dependence manifests itself when the electrical charges in
a structure have a finite response time, making them appear as fixed for slow rates of

16
field change and free to move for faster rates of field change. To be fully correct, Equa-
tions (2.4), (2.5) and (2.6) should be expressed in the following form (only Equation
(2.4) is written for brevity):
Z ∞
J(t) = σ(β) · E(t − β)dβ (2.8)
0

This more complex form of the constitutive equations must be used when physical
properties are dispersive. For most GPR applications, assuming the scalar constant
form for ε, µ, σ suffices with and being the most important. For GPR, the dielectric
permittivity is an important quantity. Most often, the terms relative permittivity or
“dielectric constant” are used and defined as follows:
ε
κ= (2.9)
ε0

2.2.3 Material properties

The subject of electrical properties (ε, µ, σ) of materials is a wide-ranging topic. Back-
ground can be found in Olhoeft (1981, 1987) and Santamarina et al. (2001). Discussion
here is limited to the common basic issues. In most GPR applications, variations in ε
and σ are most important while variations in µ are seldom of concern.
Ground penetrating radar is most useful in low-electrical-loss materials. If σ = 0, GPR
would see very broad use since signals would penetrate to great depth. In practice,
low-electrical-loss conditions are not prevalent. Clay-rich environments or areas of
saline groundwater can create conditions where GPR signal penetration is very lim-
ited.
Earth materials are invariably composites of many other materials or components.
Water and ice represent the few cases where a single component is primarily present.
A simple beach sand is a mixture of soil grains, air, water, and ions dissolved in water.
Soil grains will typically occupy 60-80% of the available volume. Understanding the
physical properties of mixtures is thus a key factor in the interpretation of a GPR
response
Mixtures of materials seldom exhibit properties directly in proportion to the volume
fraction of the constituent components. In many respects, this complexity can make
quantitative analysis of GPR data impossible without ancillary information.
Although the subject of mixtures is complex, the big picture of GPR perspective is
simpler. In the 10-1000 MHz frequency range, the presence or absence of water in the
material dominates behavior with the general picture being as follows:

• Bulk minerals and aggregates in mixtures generally are good dielectric insulators.
They typically have a permittivity in the range of 3-8 (depending on mineralogy
and compaction) and are usually insulating with virtually zero conductivity.

• Soils, rocks, and construction materials have empty space between the grains
(pore space) available to be filled with air, water, or other material.

• Water is by far the most polarizable, naturally occurring material (in other words,
it has a high permittivity with κ = 80).

17
 Figure 2.2: Configuration and representation of an A-scan

• Water in the pore space normally contains ions, and the water electrical conduc-
tivity associated with ion mobility is often the dominant factor in determining
bulk material electrical conductivity. Resulting soil and rock conductivities are
typically in the 1–1000 mS/m range.

• Since water is invariably present in the pore space of natural (geologic) materials,
except in such unique situations where vacuum drying or some other mechanism
assures the total absence of water, it has a dominant effect on electrical proper-
ties.

Empirically derived forms such as the Topp relationship and variations of Archie’s law
have long demonstrated the relationship between permittivity, electrical conductivity,
and volumetric water content for soils. More advanced relationships, such as the BHS
model use effective media theory models to derive a composite material property from
constituents. Referring to the reference materials and other chapters of this text will
provide a more substantive view of this subject.

2.3 Visualizing GPR Data

The objective of GPR data presentation is to provide a display of the processed data
that is closely approximates an image of the subsurface, with the anomalies that are
associated with the objects of interest located in their proper spatial positions. Data
display is central to data interpretation. In fact, producing a good display is an integral
part of interpretation. There are three types of displays of surface data, these are: a
one-dimensional trace ( A-scan), a two dimensional cross section ( B-scan), and a three
dimensional display (C-scan) [19, 20] described briefly each of these types as follows.

2.3.1 A-scan
An A-scan (or one dimensional data presentation ) is obtained by a stationary measure-
ment, emission and collection of a signal after placing the antenna above the position
of interest. The collected signal is presented as signal strength vs.time delay. A single
waveform b(xi , yj , t) recorded by a GPR, with the antennas at a given fixed position
(xi , yj ) is referred to as an A-scan in the above Figure 2.2. The only variable is the
time, which is related to the depth by the propagation velocity of the EM waves in
the medium.

18
2.3.2 B-scan
B-scan (or two dimensional data presentation ) signal is obtained as the horizontal
collection from the ensemble of A-scans. The horizontal axis of the two dimension
image is surface position, and the vertical axis is the round-trip marvel time of the
electromagnetic wave. When moving the GPR antennas on a line along the x-axis, one
can gather a set of A-scans, which form a two dimensional data set b(x, yj, t) , called a
B-scan Figure 2.3 (a)). dimensional dataset, denominated by the acoustic terminology
A- B- and C-scans. When the amplitude of the received signal is represented by a
colour scale (or grey-scale) a 2D image as shown in Figure 2.3 (b) is obtained. The
2D image represents a vertical slice in the ground. The time axis or the related depth
axis is usually pointed downwards. Reflections on a point scatterer located below the
surface appear, due to the beamwidth of the transmitting and the receiving antenna,
as hyperbolic structures in a B-scan.

2.3.3 C-scan
C-scan (or three dimensional data presentation ) signal is obtained from the ensemble of
B-scans, measured by repeated line scans along the plane. Three dimensional displays
are fundamentally block views of GPR traces that are recorded at different positions on
the surface. Obtaining good three-dimensional images are very useful for interpreting
specific targets. Targets of interest are generally easier to identify and isolate on three
dimensional data sets than on conventional two dimensional profile lines. Finally,
when collecting multiple parallel B-scans or in other words, when moving the antenna
over a (regular) grid in the xy-plane, a three dimensional data set b(x, y, t) can be
recorded, called a C-scan (Figure 2-4). Usually a C-scan is represented as a two
dimensional image by plotting the amplitudes of the recorded data at a given time ti .
The image b(x, y, ti) represents then a horizontal slice at a certain depth, parallel to the
recording plane Figure 2.3. Nowadays, many user-software packages have integrated
functions to plot directly three-dimensional representations of the recorded C-scans.
GPR technologies have the following advantageous features.

• high-resolution electromagnetic technique designed primarily to investigate the

shallow subsurface of the earth, buried materials .

• Fast data acquisition high resolving and non destructive electromagnetic method

• Rapid ground coverage afforded by towing the antennae either by hand or from
a vehicle.

• High lateral resolution of targets, even for larger surveys.

• The instant graphic display offered by most GPR systems allows on-site inter-
pretation

19
Figure 2.3: (a) Multiple A-scans forming a B-scan (b) Representation of a B-scan on
a grey-scale

Figure 2.4: Configuration and representation of C-scan[21]

20
 Chapter 3

Sub Array Processing Technique

3.1 Array Signal Processing

Sensor array have been in use for several decades in many practical signal processing
applications. Such an array consists of a set of sensors that are spatially distributed
at known locations with reference to a common reference point.These sensors collect
signal from source in their field of view.Depending on the sensor characteristics and
the path of propagation,the source wave forms undergo deterministic and/or random
modifications. The sensor outputs are composed of these source components and ad-
ditive noise such as measurement and thermal noise.
In active sensing situation such as radar and sonar,a known waveform of finite dura-
tion is generated which in turn propagates through a medium and is reflected by some
target back to the origin.The transmitted signal is usually modified both in amplitude
and phase by the target characteristics,which by themselves might be changing with
time and its position in space.These disturbances give rise to a random return signal.In
the passive context the signal received at the array is self generated by the target,such
as propeller or engine noise from submarines in the case of sonar.Once again the signals
are random in nature.In addition to these target-generated direct signals,there may also
be spurious returns such as cluster in the case of radar.Moreover,signals from a target
can undergo reflection,creating multiple returns that are delayed,amplitude-weighted
replicas of the direct signal to the array. These,as well as international jamming sig-
nals generate coherent interference.For example ,in radar, multipath returns give rise
to secondary signals that are completely coherent with the original signal. In general,
the signals may be uncorrelated, partially correlated or completely coherent with each
other.Similarly,the additive noise contained at the sensor outputs may also be uncorre-
lated or correlated with each other and may have equal/unequal noise variances.Since
the natural causes and responsible for signals and noise are often unrelated,it is cus-
tomary to assume that the signals and noise are uncorrelated with each other.
The practical problem of interest in array signal processing are extracting the desired
parameter such as direction of arrival, power levels and cross correlations of the signals
present in scene from the available information including the measured data [22].Often
one may also be specifically interested in the actual signal of one of these sources,and
in the case it is necessary to estimate the actual waveform associated with the desired
signal by improving the overall reception in an environment having several sources.To
achieve this,ideally it should be possible to suppress the undesired signal and enhance
the desired signal

21
3.2 Sub Array Processing in DOA Estimation
In applications involving smart antennas, and in the presence of several transmitters
operating simultaneously, it is important for a receiving array to be able to estimate
the angles of arrival, in order to decipher how many emitters are present and what
are their possible locations. In the scenario considered in this paper, the sources are
the currents induced on buried objects and a set of electric-field values, measured by a
linear array of receiving antennas parallel to the air– soil interface, is used to estimate
the directions of arrival (DOAs) [23] of the electric field back-scattered by the targets.
DOA algorithms assume that the sources are in the far-field region of the receiving
array, so that the received wave front can be considered as planar and the main an-
gular direction of the field can be estimated. However, in electromagnetic sensing of
subsurface utilities, the scatterers are usually quite near to the antennas. In order to
be able to work in near-field conditions, a sub-array processing (SAP) approach [24]
is adopted in this paper: the whole receiving array is partitioned in small sub-arrays
and for each sub-array the dominant DOA is estimated. In this way, it is possible
to successfully locate also objects that are in the far-field of the sub-array, although
being in the near-field of the whole array. The geometrical configuration of the which

Figure 3.1: Example of geometrical configuration SAP

depicted in Figure 3.1 [25] A set of M circular cylinders with radii am and centers in
(zm , xm ), m=1,. . . ,M, (where zm is the ground distance of the mth cylinder axis from
the center of the array and xm the cylinder burial depth) is buried in a half-space
characterized by relative dielectric permittivity r . The targets are illuminated by
means of a plane wave and the scattered electric field is collected by means of an array
of Narray equally-spaced antennas, which is located on the air–ground interface; it is
portioned in sub-arrays composed by N elements. We consider the electromagnetic
field scattered by the M buried objects and impinging on the receiving antennas as
composed by M narrow-band plane waves (i.e., it is assumed that the objects are lo-
cated far enough from the sub-arrays to allow the use of the far-field approximation),
with center angular frequency ω. The voltages at the output of the N elements of a
sub-array can be expressed as
     
y1 s1 n1
 ..   ..   .. 
y =  .  = A  .  +  .  = As+n (3.1)
yn sM nN
22
where sm are the signals impinging on the first element of the sub-array assumed as
reference, n is the vector of the output noise at the sub-array elements, and the steering
matrix A in equation (3.1) is defined as
h i
A = a(θ1 )a(θ1 ) · · · a(θM ) (3.2)

where the steering vectors are given by a(θ1 ) = [1, e−jkdsinθ , · · · e−jk(N −1)dsinθ]t , being d
1
the array spacing, (. )t the transposition operator, and k = ω(r 0 µ0 ) 2 the wave number
in the ground. The vectory is pre-processed in order to compute the correlation matrix

R = E{yyH } = ARSS AH + Rnn (3.3)

where RSS = E{SSH } and Rnn = E{nnH } are the source and noise correlation
matrices,E{. } indicates the expected value, and (. )H symbolizes Hermitian transpo-
sition. The DOAs found by the sub-arrays can be triangulated, obtaining a set of
crossings with intersections condensed around object locations. This approach had
already been proposed for the localization of a perfectly-conducting [26] or dielectric
[27] cylindrical object buried in a dielectric half-space; different DOA algorithms had
been implemented and compared. The localization procedure showed a good capability
in estimating the obstacle position and the well-known MUltiple SIgnal Classification
(MUSIC) algorithm[25] [28] showed to be one of the most robust approaches for DOA
estimation.

3.2.1 Direction of Arrival(DOA)

The direction-of-arrival (DOA) source localization problem is directly applicable to
many fields including radar, sonar, communications, and the geological and biomed-
ical sciences. The goal of the DOA problem in radar applications is to accurately
locate all sources with an array. Various spectral estimation methods have been ap-
plied to solving this problem and may be classified as nonparametric, parametric, and
semiparametric [29]. Nonparametric methods include the well-known delay- and-sum
(DAS) beamformer and Capon beamformer. The DAS beamformer is a computation-
ally efficient and robust method, but suffers in performance due to high sidelobes, low
resolution, and poor accuracy. The data adaptive Capon beamformer improves on
the DAS method, but it requires noncoherent sources. Parametric methods, such as
multiple signal classification (MUSIC), can provide higher resolution than the non-
parametric methods. Unfortunately, these methods are typically very sensitive to
modeling errors and can even fail if coherent sources are present [30]. The sources
in the underwater acoustic environment can be coherent due to multipath reflections
from the ocean floor or surface. Recently, a new semiparametric algorithm, referred to
as the sparse iterative covariance based estimation (SPICE), has been proposed [31].
SPICE offers superior resolution and low sidelobe levels, simultaneously retaining ro-
bustness against correlated sources. In signal processing literature,[32] direction of
arrival (DOA) denotes the direction from which usually a propagating wave arrives at
a point, where usually a set of sensors are located.Various techniques for calculating
the direction of arrival, such as Angle of Arrival (AoA), Time Difference of Arrival
(TDOA), Frequency Difference of Arrival (FDOA) and/or other similar associated
techniques.
Angle of arrival (AoA) measurement is a method for determining the direction of
23
propagation of a radio-frequency wave incident on an antenna array or determined
from maximum signal strength during antenna rotation. AoA determines the direc-
tion by measuring the Time Difference of Arrival (TDOA) at individual elements of
the array from these delays the AoA can be calculated.

3.2.2 DOA Estimation Algorithms

The DOA estimation in array signal processing is one of the important and emerging
research area. The effectiveness of this direction of arrival estimation greatly deter-
mines the performance of smart antennas. It works on the digitized output from each
sensor array. The estimation results for coherent signals, broadband signals and multi-
ple signals are of greater consideration. Various information of the parameters relating
to a particular wave can be obtained by analyzing the incoming wave that is received
by a sensor or an array of sensors. The accuracy in the estimation of direction of ar-
rival is very crucial in array signal processing [33]. In the DOA estimation algorithm,
the spectral technique is a method of estimating the DOA as θ̂ that has the maximum
value on the spectrum with the incident direction θ as a variable. The spectral meth-
ods are divided into beamforming technique, subspace-based method and Artificial
Neural Networks depending on the estimation method of the spectrum. This paper
presents an overview of direction of arrival estimation using the MUSIC algorithm.

3.2.3 MUSIC Algorithm

The name MUSIC is an acronym for MUltiple SIgnal Classification. The MUSIC
algorithm estimates the pseudospectrum from a signal or a correlation matrix using
Schmidt’s eigenspace analysis method [25, 34]. The algorithm performs eigenspace
analysis of the signal’s correlation matrix in order to estimate the signal’s frequency
content. This algorithm is particularly suitable for signals that are the sum of sinu-
soids with additive white Gaussian noise. The eigenvalues and eigenvectors of the
signal’s correlation matrix are estimated if not supply the correlation matrix. MUSIC
is a well-known high-resolution technique able to estimate the number of signals and
the DOAs as maxima of the pseudospectrum function, and also the strengths of the
waveforms. MUSIC makes the assumption that the noise in each element of the array
is uncorrelated, so the noise correlation matrix is diagonal. The incident signals, in-
stead, may be somewhat correlated creating a nondiagonal signal correlation matrix;
however, under high signal correlation the traditional MUSIC algorithm breaks down
and other methods must be implemented to correct this weakness.
In the presence of M signals, the number of signal eigenvalues and eigenvectors is M
and the number of noise eigenvalues and eigenvectors is N − M . MUSIC exploits the
noise eigenvector subspace. The N − M smallest eigenvalues of R are associated with
the noise and they are all equal to the noise variance. Therefore, it is possible to con-
struct a N ×(N −M ) subspace spanned by the noise eigenvectors E , orthogonal to the
array steering vectors at the angles of arrival θ1 , θ2 · · · ,hbecause of this orthogonality
i
condition, it can be shown that the Euclidean distance aH (θ)EEH aH (θ) is equal to
zero for each DOA. Finally, the MUSIC pseudo spectrum is given by
h i−1
P (θ) = aH (θ)EEH aH (θ) (3.4)

24
 Figure 3.2: Angle of arrival and the sharp pseudo spectrum around angle of arrival

with sharp peaks at the angles of arrival.As indicated in the Figure 3.2 the algorithm
make sharp peak around the three angle of arrival which is (-20◦ , 10◦ , 60◦ ) based on
random signal generated for simulation purpose.

3.3 Factors Influencing the Efficiency of MUSIC

Algorithm
The application environment and the incoming signals are not the only parameters to
influence the operation of MUSIC algorithm.The following parameters can vary the
estimation efficiency.

• Spacing between the array elements.

• Angular separation between the incident signals.

• Number of samples taken for the incident signals.

• Signal-to-Noise Ratio.

• Number of array elements.

• Mutual Coupling between the sensor arrays.

MUSIC algorithms has advantages over other estimation algorithms because of the
sharp needle spectrum peaks which can efficiently estimate the independent source
25
signals with high precisions unlike the other estimation processes which are limited
with low precisions. It has many practical applications as it provides unbiased estima-
tion results. The MUSIC algorithm to estimate the direction has even proved to have
better performance in a multiple signal environment. MUSIC algorithm has better
resolution, higher precision and accuracy with multiple signals. But this algorithm
achieves high resolution in DOA estimation [35] only when the signals being incident
on the sensor array are non-coherent. It losses efficiency when the signals are coher-
ent. For coherent signals the conventional MUSIC algorithm fails to obtain narrow
and sharp peaks as indicated in Figure 3.2.

26
 Chapter 4

GPR Modeling and Utility

Detection

4.1 GPR modelling Basic concepts

All electromagnetic phenomena, on a macroscopic scale, are described by the well-
known Maxwell’s equations. These are first order partial differential equations which
express the relations between the fundamental electromagnetic field quantities and
their dependence on their sources.
∂B
∇×E =− (4.1)
∂t
∂D
∇×H =− + Js (4.2)
∂t
∇·B =0 (4.3)
∇ · D = qv (4.4)
where t is time (seconds) and qv is the volume electric charge density (coulombs/cubic
metre), Js is current density. In Maxwell’s equations, the field vectors are assumed to
be single-valued, bounded, continuous functions of position and time. In order to sim-
ulate the GPR response from a particular target or set of targets the above equations
have to be solved subject to the geometry of the problem and the initial conditions[36].
The nature of the GPR forward problem classifies it as an initial value – open bound-
ary problem. This means that in order to obtain a solution you have to define an
initial condition (i.e. excitation of the GPR transmitting antenna) and allow for the
resulting fields to propagate through space reaching a zero value at infinity since, there
is no specific boundary which limits the geometry of the problem and where the elec-
tromagnetic fields can take a predetermined value. Although the first part is easy
to accommodate (i.e. specification of the source), the second part cannot be easily
tackled using a finite computational space.
The FDTD approach to the numerical solution of Maxwell’s equations is to discretize
both the space and time continua. Thus the discretization spatial ∆x, ∆y and ∆z
and temporal ∆t steps play a very significant role – since the smaller they are the
closer the FDTD model is to a real representation of the problem. However, the val-
ues of the discretization steps always have to be finite, since computers have a limited
amount of storage and finite processing speed. Hence, the FDTD model represents a
discretized version of the real problem and is of limited size. Building block of this
discretized FDTD grid is the Yee cell [37] named after Kane Yee who pioneered the
27
Figure 4.1: Single FDTD Yee cell showing electric (red) and magnetic (green) field
components

FDTD method. This is illustrated for the 3D case in Figure 4.1

By assigning appropriate constitutive parameters to the locations of the electromag-
netic field components complex shaped targets can be included easily in the models.
However, objects with curved boundaries are represented using a staircase approxima-
tion

4.1.1 gprMax
gprMax is open source software that simulates electromagnetic wave propagation. It
solves Maxwell’s equations in 3D using the Finite-Difference Time-Domain (FDTD)
method. gprMax was designed for modelling Ground Penetrating Radar (GPR) but
can also be used to model electromagnetic wave propagation for many other applica-
tions. gprMax is currently released under the GNU General Public License v3 or higher
[36]. GprMax2D and GprMax3D are two computer programmes that implement the
FDTD method for GPR modelling in 2D and 3D, respectively. Some of their key fea-
tures are: an easy to use command interface, the ability to model dispersive materials,
the modelling of complex shaped targets as well as the simulation of unbounded space
using powerful absorbing boundary conditions. GprMax3D allows the simulation of
GPR antennas and even the introduction of their feeding transmission lines into the
model. GprMax2D is mainly used for GPR “signature” simulation whereas GprMax3D
is used for more detail and realistic simulations especially when comparisons with real
GPR data are important. Both GprMax2D and 3D programmes use a simple ASCII
(text) file to define the models parameters. In this file special commands are used
which instruct the software to perform specific functions that are required by the type
of the model the user wants to create. Some of the commands of GprMax2D are shown

28
in Table 6.1. The software package includes a comprehensive Users. Manual in which
details of the functionality of the programmes can be found. GprMax is principally
written in Python 3 with performance-critical parts written in Cython. It includes a
CPU-based solver parallelised using OpenMP, and a GPU-based solver written using
the NVIDIA CUDA programming model.
gprMax is fundamentally based on solving Maxwell’s equations in 3D using the FDTD
method - transverse electromagnetic (TEM) mode. However, it can also be used to
carry out simulations in 2D using the transverse magnetic (TM) mode. This is achieved
through specifying a single cell slice of the domain, i.e. one dimension of the domain
must be equal to the spatial discretization in that direction. When this occurs the
electric and magnetic field components on the two faces of single cell slice in the in-
variant direction are set to zero. This is illustrated for the 2D TM case in Figure 4.2.

Figure 4.2: Single FDTD Yee cell showing electric (red), magnetic (green), and zeroed
out (grey) field components for 2D transverse magnetic (TM) z-direction
mode

Using this approach means that Maxwell’s equations in 3D, shown in equation (4.5).
as six coupled partial differential equations, reduce to the corresponding 2D form - in
this case 2D TMz, shown in equation (4.6).
!
∂Ex 1 ∂Hz ∂Hy
= − − (JS )x − σEx
∂t  ∂y ∂z
!
∂Ey 1 ∂Hz ∂Hz
= − − (JS )y − σEy
∂t  ∂z ∂x
!
∂Ez 1 ∂Hy ∂Hx
= − − (JS )z − σEz
∂t  ∂z ∂x
29
 !
∂Hx 1 ∂Ey ∂Ez
= − − (MS )y − σ ∗ Hx
∂t µ ∂z ∂y
!
∂Hy 1 ∂Ez ∂Ex
= − − (MS )y − σ ∗ Hy (4.5)
∂t µ ∂x ∂z
!
∂Hy 1 ∂Ex ∂Ey
= − − (MS )z − σ ∗ Hz
∂t µ ∂x ∂z

!
∂Ez 1 ∂Hy ∂Hx
= − − (JS )z − σEz
∂t  ∂z ∂x
!
∂Hx 1 ∂Ey ∂Ez
= − − (MS )y − σ ∗ Hx (4.6)
∂t µ ∂z ∂y
!
∂Hy 1 ∂Ex ∂Ey
= − − (MS )z − σ ∗ Hz
∂t µ ∂x ∂z
These equations are discretized in both space and time and applied in each FDTD
cell. The numerical solution is obtained directly in the time domain in an iterative
fashion. In each iteration the electromagnetic fields advance (propagate) in the FDTD
grid and each iteration corresponds to an elapsed simulated time of one ∆t. Hence by
specifying the number of iterations you can instruct the FDTD solver to simulate the
fields for a given time window.
The price you have to pay for obtaining a solution directly in the time domain using
the FDTD method is that the values of ∆x, ∆y, ∆z and ∆t can not be assigned inde-
pendently. FDTD is a conditionally stable numerical process. The stability condition
is known as the CFL condition after the initials of Courant, Freidrichs and Lewy and
is given by,
1
∆t ≤ q
1 1 1
c (∆x)2 + (∆y)2 + (∆z)2

where c is the speed of light. Hence ∆t is bounded by the values of ∆x, ∆y and ∆z.
The stability condition for the 2D case is easily obtained by letting ∆z −→ ∞.

4.1.2 Coordinate System and Conventions

A right-handed Cartesian coordinate system is used with the origin of space coordinates
in the lower left corner at (0,0,0). Figure 4.3 illustrates the coordinate system of
gprMax. Only one row of cells in the x direction is depicted. The space coordinates
range from the left edge of the first cell to the right edge of the last one. Assuming
that ∆x = 1 metre, if you wanted to allocate a rectangle with its x dimension equal
to 3 metres and its lower x coordinate at 1 then the x range would be [1..4] but the
3D cells allocated by gprMax would be [1..3]. In the 3D FDTD cell there are no field
components located at the centre of the cell. Electric field components are tangential
to, and magnetic field components normal to the interfaces between cells.The field
components depicted in Figure 4.3 correspond to space coordinate 1. Source and
output points defined in space coordinates are directly converted to cell coordinates
30
 Figure 4.3: gprMax coordinate system and conventions.

and the corresponding field components.

The actual positions of field components for a given set of space coordinates (x, y, z)
are:
∆x

Ex x + , y, z
2
∆y

Ey x, y + ,z
2
∆z

Ez x, y, z +
2
∆y ∆z

Hx x, y + ,z +
2 2
∆x ∆z

Hy x, + ,z +
y 2
∆x ∆y

Hz x, + + ,z
y 2
Hertzian dipole sources as well as other electric field excitations (i.e. voltage sources,
transmission lines) are located at the corresponding electric field components.

4.1.3 Spatial Discretization

There is no specific guideline for choosing the right spatial discretization for a given
problem. In general, it depends on the required accuracy, the frequency content of the
source pulse and the size of the targets. Obviously, all targets present in a model must
be adequately resolved. This means, for example, that a cylinder with radius equal to
one or two spatial steps does not really look like a cylinder.
An other important factor which influences the spatial discretization is the errors
associated with numerically induced dispersion. This means that contrary to the real
world where electromagnetic waves propagate with the same velocity irrespectively of
their direction and frequency (assuming no dispersive media and far-field conditions) in
the discrete one this is not the case. This error (details can be found in [38, 39, 40] can
31
be kept in a minimum if the following rule-of-thumb is satisfied: The discretization step
should be at least ten times smaller than the smallest wavelength of the propagating
electromagnetic fields. The discretization step should be at least ten times smaller
than the smallest wavelength of the propagating electromagnetic fields.
λ
∆l =
10
Note that in general low-loss media wavelengths are much smaller compared to free
space.

4.1.4 Absorbing Boundary Conditions

One of the most challenging issues in modelling open boundary problems, such as
GPR, is the truncation of the computational domain at a finite distance from sources
and targets where the values of the electromagnetic fields can not be calculated directly
by the numerical method applied inside the model. Hence, an approximate condition
known as Absorbing boundary condition (ABC) is applied at a sufficient distance from
the source to truncate and therefore limit the computational space. The role of this
ABC is to absorb any waves impinging on it, hence simulating an unbounded space.
The computational space (i.e the model) limited by the ABCs should contain all im-
portant features of the model such as sources and output points and targets. Figure
4.4 illustrates this basic difference between the problem to be modelled and the actual
FDTD modelled space. It is assumed that the half-space which contains the target(s)
is of infinite extent.
Therefore, the only reflected waves will be the ones originating from the target. In
cases where the host medium is not of infinite extent (e.g. a finite concrete slab) the
assumption of infinite extent can be made as far as the actual reflections from the
slab termination are not of interest or its actual size is large enough that any reflected
waves which will originate at its termination will not affect the solution for the re-
quired time window. In general, any objects that span the size of the computational
domain (i.e. model) are assumed to extend to infinity. The only reflections which will
originate from their termination at the truncation boundaries of the model are due
to imperfections of the ABCs and in general are of a very small amplitude compared
with the reflections from target(s) inside the model. The ABCs employed in gprMax
will, in general, perform well (i.e. without introducing significant artificial reflections)
if all sources and targets are kept at least 15 cells away from them.
GprMax uses Perfectly Matched Layer (PML) ABCs based on a recursive integration
approach to the complex frequency shifted (CFS) PML [38]. A general formulation
of this RIPML, which can be used to develop any order of PML, has been used to
implement first and second order CFS stretching functions. One of the attractions of
the RIPML is that it is easily applied as a correction to the field quantities after the
complete FDTD grid has been updated using the standard FDTD update equations.
The cells of the RIPML, which have a user adjustable thickness, very efficiently absorb
most waves that propagate in them. Although, source and output points can be spec-
ified inside these cells it is wrong to do so from the point of view of correct modelling.
The fields inside these cells are not of interest to GPR modelling.Placing sources in-
side these cells could have effects that have not been studied and will certainly provide
erroneous results from the perspective of GPR modelling. The requirement to keep
32
Figure 4.4: GPR forward problem showing computational domain bounded by Ab-
sorbing Boundary Conditions (ABCs)

sources and targets at least 15 cells away for the PML has to be taken into account
when deciding the size of the model domain. Additionally, free space (i.e. air) should
be always included above a source for at least 15-20 cells in GPR models. Obviously,
the more cells there are between observation points, sources, targets and the absorbing
boundaries, the better the results will be. gprMax now offers the ability (for advanced
users) to customize the parameters of the PML which allows its performance to be
better optimized for specific applications. All other boundary conditions which apply
at interfaces between different media in the FDTD model are automatically enforced
in gprMax. In appendix some of the commands in Table 6.1

33
4.2 Modelling of single buried metallic object
The flexibility of GprMax2D allows the modelling of complex what-if scenarios. In
the following a simple example of modelling considered scenario was a homogeneous
dielectric medium concrete hosting circular section cylindrical targets. The geometry
of the model consists of a 0.3m wide concrete slab, where the metallic cylinder of
0.005m diameter is located under the receiving array at the coordinate of (x0 , y0 , z0 ) =
(0.15, 0, 0.0625)m and (x1 , y1 , z1 ) = (0.15, 0.002, 0.0625)m and at a depth h = 0.0625m
from the slab’s surface which also located at the midst of slab,Waveform gaussian of
type gaussiandot with maximum amplitude scaling 1, frequency 1.5GHz. Gaussiandot
is first derivative of a Gaussian waveform given by
2
W (t) = 2ζ(t − χ)e−ζ(t−χ) (4.7)

where ζ = 2π 2 f 2 , χ = f1 and f is the frequency. Hertzian dipole with polarity z

at 0.144m, 0.074m, 0.134m, using waveform gaussian.Hertzian dipole allows you to
specify a current density term at an electric field location the simplest excitation,
often referred to as an additive or soft source.
I∆l
Js = (4.8)
∆y∆x∆z
, where Js is the current density, I is the current, ∆l is the length of the infinitesimal
electric dipole, and ∆x, ∆y, and ∆z are the spatial resolution of the grid. To determine
the smallest wavelength, the highest frequency and lowest velocity present in the model
are required. In this model by examining the spectrum of a waveform it is evident
much higher frequencies are present,the highest significant frequency present in the
model is likely to be around 4GHz. The wavelength at 1.5GHz and r = 7 in the
half-space (which has the lowest velocity) would be:
c 299792458
λ= √ = √ ≈ 28.33mm
f r 4 × 109 7
This would give a minimum spatial resolution of approximately 3 mm. An illustration
of the model’s geometry is presented in Figure. 4.5 and Table 4.1 describes modelling
materials.

Table 4.1: gprMax2D Model list of used materials

ID Materials Type r σ[S/m] µr σr∗ [Ohm/m] Dielectric
Name smoothable
0 Pec Builtin 1 inf 1 0 False
1 free space Builtin 1 0 1 0 True
2 concrete 7 0.05 3 0 True
3 free dielectric 4 0.025 2 0 True
space+concrete smoothed

Receiver array 0.164m, 0.072m, 0.134m, to 0.164m, 0.076m, 0.134m with steps 0.002m,
34
 Figure 4.5: gprMax2D Model of first scenario

0.002m, 0.002m, Material concrete with r = 7, σ = 0.05 S/m; µr = 3, σ∗ = 0Ω/m

created. Geometry view from 0m, 0m, 0m, to 0.999m, 0.15m, 0.2m, discretisation
0.002m, 0.002m, 0.002m, with file name base Cylinder in concrete created. Box from
0m, 0m, 0m, to 0.3m, 0.15m, 0.124m of material(s) concrete created, dielectric smooth-
ing is on. Cylinder with face centres 0.15m, 0m, 0.062m and 0.15m, 0.15m, 0.062m,
with radius 0.005m, of material(s) pec created.

4.3 Modeling of of multiple buried metallic objecs at

the same depth
Second model is created using the same materials used in first model with different
domain size and spatial discretization. Spatial discretisation:0.003 × 0.003 × 0.003m,
Waveform gaussian of type gaussiandot with maximum amplitude scaling 1, frequency
1.5GHz created. Hertzian dipole with polarity z at 0.039m, 0m, 0.321m, using wave-
form gaussiandot created. Receiver array 0.051m, 0m, 0.321m, to 0.057m, 0m, 0.321m
with steps 0.003m, 0.003m, 0.003m, Material concrete with r = 7, σ = 0.05 S/m;
µr =3, σ ∗ = 0 Ohm/m created. In box from 0m, 0m, 0m, to 0.999m, 0.003m, 0.3m of
material(s) concrete three cylinder;as depicted in figure 4.6 cylinder 1, cylinder 2 and
cylinder 3 with the same radius of 0.01m, located horizontally varying and at the same
depth at (0.15 0 0.0625),(0.5 0 0.0625 )and (0.85 0 0.0625)m respectively created.

4.3.1 Multiple buried metallic objects at different depth and size

Waveform gaussian of type gaussiandot with maximum amplitude scaling 1, frequency
1.5GHz created. Hertzian dipole with polarity y at 0.04m, 0m, 0.32m, using waveform
gaussian created. Receiver array 0.05m, 0m, 0.32m, to 0.054m, 0m, 0.32m with steps
0.002m, 0.002m, 0.002m, Material concrete with r = 7, σ = 0.05 S/m; µr =3, σ ∗ = 0
Ohm/m created, as depicted on figure 4.7 Geometry view from 0m, 0m, 0m, to 1m,
35
 Figure 4.6: Three cylinder buried in concrete at the same depth

Figure 4.7: Different size cylinder at different depth

0.002m, 0.4m, discretisation 0.002m, 0.002m, 0.002m, with Box from 0m, 0m, 0m, to
1m, 0.002m, 0.3m of material(s) concrete created, Cylinder with face centres 0.15m,
0m, 0.15m and 0.15m, 0.002m, 0.15m, with radius 0.03m, of material(s) pec created,
Cylinder with face centres 0.5m, 0m, 0.25m and 0.5m, 0.002m, 0.25m, with radius
0.01m, of material(s) pec created, Cylinder with face centres 0.85m, 0m, 0.2m and
0.85m, 0.002m, 0.2m, with radius 0.02m, of material(s) pec created.

4.3.2 3D Model of very complex environment

This model demonstrates how to build a more realistic soil model using a stochastic
distribution of dielectric properties with s metallic cylinder buried in heterogeneous
medium with zero magnetic loss (σ ∗ =0 Ohm/m) as shown in Figure 4.8 and and Table
4.2. A mixing model for soils proposed by Peplinski [41] is used to define a series of
36
dispersive material properties for the soil.

Table 4.2: 3D Heterogeneous model list of used materials

ID Name Type r σ ∆r tau[s] µr Dielectric
smoothable
0 Pec builtin 1 inf 1 False
1 free builtin 1 0 1 True
Space
2 0.0148 debye 4.215 0.0212865 0.156189 9.231e-12 1 False
3 0.0425 debye 5.29302 0.0347966 0.731532 9.231e-12 1 False
4 0.0702 debye 6.2254 0.0439739 1.52628 9.231e-12 1 False
5 0.0978 debye 7.0781 0.0513561 2.48534 9.231e-12 1 False
6 0.1255 debye 7.87761 0.0576883 3.58129 9.231e-12 1 False
7 0.1532 debye 8.6382 0.0633114 4.79683 9.231e-12 1 False
8 0.1808 debye 9.36871 0.0684149 6.11979 9.231e-12 1 False
9 0.2085 debye 10.0751 0.0731167 7.54102 9.231e-12 1 False
10 0.2362 debye 10.7615 0.0774962 9.05329 9.231e-12 1 False
11 0.2638 debye 11.4311 0.0816098 10.6507 9.231e-12 1 False

Waveform ricker of type ricker with maximum amplitude scaling 1, frequency 1.5GHz
created. Waveform ricker of type ricker with maximum amplitude scaling 1, frequency
1.5e+09Hz created. A Ricker (or Mexican Hat) waveform which is the negative, nor-
malized second derivative of a Gaussian waveform given by
2
W (t) = −(2ζ(t − χ)2 − 1)e−ζ(t−χ) (4.9)
√
where ζ = π 2 f 2 , χ = f2 and f is the frequency. Hertzian dipole with polarity y at
0.024m, 0.024m, 0.08m, using waveform ricker created. The inclusion of improved
models of soils is important for many GPR simulations. gprMax can now be used to
create soils with more realistic dielectric and geometrical properties. A semi-empirical
model, initially suggested by [42], is used to describe the dielectric properties of the soil.
The model relates relative permittivity of the soil to bulk density, sand particle den-
sity, sand fraction, clay fraction and water volumetric fraction. Using this approach,
a more realistic soil with a stochastic distribution of the aforementioned parameters
can be modelled. The real and imaginary parts of this semi-empirical model can be
approximated using a multipole Debye function plus a conductive term. This can now
be achieved in gprMax using the new dispersive material functionality. gprMax has
always included the ability to represent dispersive materials using a single-pole Debye
model. Many materials can be adequately represented using this approach for the
typical frequency ranges associated with GPR. However, multi-pole Debye, Drude and
Lorentz functions are often used to simulate the electric susceptibility of materials such
as: water[43] and soils[44]. Electric susceptibility relates the polarization density to
the electric field, and includes both the real and imaginary parts of the complex electric
permittivity variation. In the new version of gprMax a recursive convolution based
37
 Figure 4.8: 3D First Model of heterogeneous environment

method is used to express dispersive properties as apparent current density sources.A

major advantage of this implementation is that it creates an inclusive susceptibility
function that holds, as special cases, Debye, Drude and Lorentz materials. gprMax
Allows you to introduce a material into the model described by a set of constitutive pa-
rameters based on a multiple pole Debye formulation. For example, the susceptability
function for a single pole Debye material is given by:
∆rp −t
χp (t) = e τp , (4.10)
τp
where ∆rp = rsp − r∞ , rsp is the zero-frequency relative permittivity for the pole,
r∞ is the relative permittivity at infinite frequency, and τp is the pole relaxation time.
Mixing model (Peplinski) used to create soil properties with sand fraction 0.5, clay
fraction 0.5, bulk density 2g/cm3, sand particle density 2.66g/cm3, and water volu-
metric fraction 0.001 to 0.25 created. Fractal box soil from 0m, 0m, 0m, to 0.4m, 0.05m,
0.074m with soil properties, fractal dimension 1.5, fractal weightings 1, 1, 1, fractal
seeding None, with 10 material(s) created. Fractal surface from 0m, 0m, 0.074m, to
0.4m, 0.05m, 0.074m with fractal dimension 1, fractal weightings 1, 1, fractal seeding
None, and range 0.05m to 0.0525m, added to soil. Cylinder with face centres 0.2m,
0m, 0.038m and 0.2m, 0.05m, 0.038m, with radius 0.015m, of material pec created.
The Last model is metallic cylinder buried in a very complex environment with zero
magnetic loss (σ ∗ =0 Ohm/m) as shown in Figure 4.9 and and table 4.3. Waveform
ricker of type ricker with maximum amplitude scaling 1, frequency 1.5GHz created.
Hertzian dipole with polarity y at 0.1m, 0.15m, 0.152m, using waveform my ricker
created. Receiver array 0.11m, 0.15m, 0.152m, to 0.112m, 0.15m, 0.152m with steps
0.002m, 0.002m, 0.002m. Mixing model (Peplinski) used to create soil properties
with sand fraction 0.5, clay fraction 0.5, bulk density 2g/cm3, sand particle density
2.66g/cm3, and water volumetric fraction 0.001 to 0.25 created,Fractal box soil from
0m, 0m, 0m, to 0.2m, 0.2m, 0.15m with soil properties, fractal dimension 1.5, fractal
weightings 1, 1, 1, fractal seeding None, with 10 material(s) created. Fractal surface
from 0m, 0m, 0.15m, to 0.2m, 0.2m, 0.15m with fractal dimension 1, fractal weightings
1, 1, fractal seeding None, and range 0.1m to 0.15m, added to soil. Water on surface
from 0m, 0m, 0.15m, to 0.2m, 0.2m, 0.15m with depth 0.124m, added to soil. Water is
modelled using a sigle pole Debye formulation equation 4.10 with properties rsp =80.1,
38
r∞ =4.9, and a relaxation time of τ =9.231 × 10− 12 seconds. 5000 blades of grass on
surface from 0m, 0m, 0.15m, to 0.2m, 0.2m, 0.15m with fractal dimension 1, fractal
seeding None, and range 0.151m to 0.1515m, added to soil. The grass is modelled
using a sigle pole Debye formulation with properties rsp =18.5087, r∞ =12.7174, and
a relaxation time of τ =1.0793 × 10− 11 seconds. The geometry of the blades of grass

Table 4.3: 3D Very complex medium model list of used materials

ID Name Type r σ ∆r tau[s] µr Dielectric
smoothable
0 Pec builtin 1 inf 1 False
1 free builtin 1 0 1 True
Space
2 0.0148 debye 4.215 0.0212865 0.156189 9.231e-12 1 False
3 0.0425 debye 5.29302 0.0347966 0.731532 9.231e-12 1 False
4 0.0702 debye 6.2254 0.0439739 1.52628 9.231e-12 1 False
5 0.0978 debye 7.0781 0.0513561 2.48534 9.231e-12 1 False
6 0.1255 debye 7.87761 0.0576883 3.58129 9.231e-12 1 False
7 0.1532 debye 8.6382 0.0633114 4.79683 9.231e-12 1 False
8 0.1808 debye 9.36871 0.0684149 6.11979 9.231e-12 1 False
9 0.2085 debye 10.0751 0.0731167 7.54102 9.231e-12 1 False
10 0.2362 debye 10.7615 0.0774962 9.05329 9.231e-12 1 False
11 0.2638 debye 11.4311 0.0816098 10.6507 9.231e-12 1 False
12 water builtin, 4.9 0 75.2 9.231e-12 1 False
debye
13 grass builtin, 12.7174 0 5.7913 1.0793e-11 1 False
debye

are defined by the parametric equations:

!2 !2
t t
x = xc + s x , y = y c + sy
bx by

and z = t where sx and sy can be -1 or 1 which are randomly chosen, and where the
constants bx and by are random numbers based on a Gaussian distribution. Cylinder
with face centres 0.1m, 0m, 0.074m and 0.1m, 0.2m, 0.074m, with radius 0.01m, of
material(s) pec created.

39
Figure 4.9: 3D Model of very complex environment

40
 Chapter 5

Result and Discussion

5.1 Detection of buried metallic object

5.1.1 Result I
Under this Chapter we are going to discuss the response to what’s modelled in Chapter
Four. The first simulation scenario consists of a single metallic cylinder with radius
0.005m buried in a concrete with relative dielectric permittivity r = 7. From the geom-
etry shown in Figure 4.5 a trace number 450 or B-scan is created.A B-scan is composed
of multiple traces (A-scans) recorded as the source and receiver are moved over the
target, in this case the metal cylinder.Figure 4.5 shows array antenna with three ele-
ments used to scan the buried cylindrical target embedded in homogeneous concrete
environment. Three Receiver antennas located at 0.164m, 0.072m, 0.134m,0.164m,
0.074m, 0.134m and at 0.164m, 0.076m, 0.134m with output component(s) Ex , Ey , Ez
Hx , Hy , Hz and Ex , Ez , Hx , Hy and Hz = 0 responses are created. Again from Fig-
ure 5.1, the initial part of the signal (∼ 0.5 − 1ns) represents the direct wave from
transmitter and ground surface reflected received. Then (∼ 3 − 3.85ns) indicates the
refelected wave from the metal cylinder which creates the hyperbolic shape between
approximately 175 and 275 trace number.

5.1.2 Result II
The response of second model depicted in Figure 4.6 three cylinder buried in concrete
at the same depth reflects received signal at (0.051m, 0m, 0.321), (0.054m, 0m, 0.321)
and (0.057m, 0m, 0.321) with (Ex , Ey , Ez ), (Hx , Hy , Hz) with (Ex , Ez , Hx , Hy and
Hz = 0) responses are created. Figure 5.2 shows the received signal at the range of
time which clearly indicates that the array outputs from target of interest are inverted
U shape. the initial part of the signal (∼ 0.5 − 1ns) represents the direct wave from
transmitter and ground surface reflected received. Then (∼ 2 − 2.75ns) indicates the
refelected wave from the metal cylinder which creates the hyperbolic shape between 0
to 100, 200 to 280 and 380 to 460 trace number. From the geometry shown in Figure
4.1 a B-scan is created. A B-scan is composed of multiple traces (A-scans) recorded
as the source and receiver are moved over the target, in this case the metal cylinder.
Figure 5.2 shows array antenna with three elements used to scan the buried cylindrical
target embedded in homogeneous concrete environment and the B-scan (of the Ey field
component). Again, the initial part of the signal (∼ 0.5 − 1ns) represents the direct
wave from transmitter to receiver. Then comes the refelected wave (∼ 2 − 2.75ns)
from the metal cylinder which creates the hyperbolic shape.

41
 (a) First receiver antenna output

(b) Second Antenna output (c) Third antenna output

Figure 5.1: B-scan of model of a single metallic cylinder buried in concrete

5.1.3 Result III

The third scenario was modelled to analyze the response at different depth as depicted
in Figure 4.3 and as we understand from the modelled geometry three cylindrical
metallic having radius 0.03m, 0.01 and 0.02m at depth of 0.15m, 0.25m and 0.2m and
located at different horizontal location. From the result depicted in Figure 5.3 it is
clearly visible that the response of the three metallic cylinder found at different depths
are having different time response which is helpful to estimate the depth of utility.The
first target of interest found between 0 to 100 trace number approximately received
at 4.25-4.75ns. The second target approximately found at 190 to 280 trace of number
and the received time approximately 1.85-2.35ns.and the last target of interest found
at 360 to 450 trace of number and the received time approximately 3-3.5ns.

5.1.4 Horizontal Location Estimation

From Result I Figure 5.1 the horizontal location can be estimated using;Antenna po-
sition,stepping size and trace number. The first 150 trace or B- scan have no signal
received by the array antenna. Antenna array started receiving signal approximately
around 175 trace of number and the electric field strength increases as scanner move
forward up to approximately 275 trace of number. After 275 trace of number no signal
42
 (a) First antenna output

(b) Second Antenna output

(c) Third antenna output

Figure 5.2: B-scan of model of a three metallic cylinder buried in concrete at the same
depth

43
 (a) First antenna output

(b) Second Antenna output

(c) Third antenna output

Figure 5.3: B-scan of model of a three metallic cylinder buried in concrete at different
depth
44
received at all and hence we can see that the target found clearly in between 175-275
trace of number from 450 total trace of number. So we can approximate the horizontal
location of the target by averaging the location of received signal from initial to final
point. Horizontal location of cylinder is (0.5 0 0.0625)- (0.5 0.002 0.0625)m from input
file of simulation with step size 0.002m.
Therefore:

Initial location of each antenna + Initial number of trace × step size = initial
location of received signal

From the simulation result: First Locaion: location of antenna1 + 175 trace number
× 0.002 = 0.03 + 0.05 = 0.38 this is the initial location to receive the signal
Second location: 0.03 + 200 trace number ×0.002 = 0.03 + 0.4 = 0.43
Third location: 0.03 + 250 trace number ×0.002 = 0.53 is the exact horizontal location
of the target as per simulation result.
Fourth Location: 0.03 + 275 trace number×0.002 = 0.58 which is very visible that the
received signal strength found between 200-250 trace number from the above signal.
So we can approximate the horizontal location of the target by averaging the location
of received signal from initial to final point.The maximum field strength found when
exactly the target of interest is right under the moving scanner.The width of the hy-
perbola also indicates the size of the target.The wider the hyperbola, the larger the
size of the target.
To estimate depth, the travel time (two-way travel time) information provided by a
radars used in conjunction with ground wave velocity, which depends on the dielectric
constant of materials, where it is usually assumed to be constant for the area under
investigation. This procedure provides satisfactory results in most cases. However,
wrong depth estimates can result in damage to public utilities, rupturing pipes, cut-
ting lines and so on. These cases occur mainly in areas that have a marked variation
of water content and/or soil lithology, thus greater care is required to determine the
depth of the targets.
The accuracy of depth estimates using GPR depends on precise knowledge of the ve-
locity of electromagnetic wave in the medium. Once the velocity is deter-mined, the
double time is converted to depth, and consequently, the depth of the object is deter-
mined. This velocity information usually comes from the adjustment of the hyperbola
equation in radar, where the identified diffraction hyperbolas can be associated with
buried targets. Depth estimates can also be obtained by means of direct information
from previously known targets, together with assumption that the velocity is constant
for the entire GPR profile. In addition to improving the depth estimates of buried
underground utilities, the interval velocity technique holds great potential for under-
standing wave velocity in the soil matrix, leading to more precise information about
soil compaction and concentration of water, thus expanding the field of applications
beyond the context of urban planning.

45
 (a) Ex component of antenna output

(b) Ey Component of antenna output

(c) Ez Component of antenna output

Figure 5.4: B-scan for model of metal cylinder buried in heterogeneous medium

46
5.1.5 Result IV
This result is found from the modelled scenario of Figure 4.8,3D First Model of het-
erogeneous environment Receiver array 0.034m, 0.022m, 0.08m, to 0.034m, 0.026m,
0.08m with steps 0.002m, 0.002m, 0.002m Receiver at (0.034m, 0.022m, 0.08m) to
(0.034m, 0.026m, 0.08m) with output component(s) Ex, Ey, Ez, created as depicted in
Figure 5.4. Due to the nature of the environment, the response to the target can not
be easily identified.Since the model is approached to realist medium it requires careful
reprocessing,ground feature extraction and very shielded antenna.

47
 Chapter 6

Conclusion and Recommendations

6.1 Conclusion
Ground penetrating radar and homogeneous medium with metallic objects were sim-
ulated using software,gprMax shows significant potential to provide genuine solutions
to subsurface sensing problems in subsurface utility detection. Direction of arrival
estimation is done using MUSIC algorithm. The very oversimplified 2D concrete slab
with metallic object embedded it modelled and the response to the target identified.
Multiple targets at the same and different depth have been modelled. The responses
to the target of interest at the same depth have the same time responses and which
at located at different depth have different time response as discussed in results. The
simulation results show that the detection of metallic buried objects is possible through
the collected energy at the receiver. The application and importance of the realistic
ground model in the 3-D FDTD simulations of GPR scenarios are presented.The earth
is modeled as a Debye medium with single poles and static conductivity and an ef-
ficient FDTD scheme is used to simulate the wave propagation inside stratified and
inhomogeneous dispersive soils. Direction of arrival Estimation with the help of Sub
Array processing techniques . The scope of applications has not yet been fully explored
or properly transitioned into the industry here in Ethiopia. Research and development
into ground penetrating radar technology needs to continue, and specific applications
should be identified where the technology can be best applied. Finally GPR method
with SAP techniques, with help of DOA estimation algorithm, MUSIC algorithms are
very useful method in subsurface utility detection and provide high image resolution.

6.2 Recommendation
Ground penetrating radar shows significant potential to provide a new class of sensing
method to facilitate many important aspects of sub surface probing. To date, how-
ever, ground penetrating radar has not made a significant impact on the Construction
agencies here in Addis Ababa, Ethiopia. Some of the key reasons include: Lack of
researches on subsurface utility detection where ground penetrating radar could in
fact be of use.Lack of clarity regarding how ground penetrating radar is deployed.The
complex nature of the received radar signals, which often require both knowledge of
electromagnetic propagation characteristics as well as geological expertise in order to
faithfully interpret the data in all but the simplest imaging configurations. The rel-
atively high capital cost of the radar equipment and antennas, compounded by its
relative fragility. The fact that there are no commonly available ground penetrating

48
radar units that are approved for use in the the Ethiopian Construction Agencies. A
number of measures are recommended in an attempt to address the above: 1) In-
creased researches through education and demonstrations have to directly given 2)
More research and development should be made into fundamental evaluations of the
technology to further characterize capability with respect to the particular ground
conditions that are typical of the subsurface scanning.

6.3 Future Works

The result of 3D Model of very complex environment depicted in Figure 4.8: re-
quired very shielded antenna and prepossessing system to extract some background
features,and to model the very realistic model with shielding antenna it also requires
13.8GB RAM.So It’s suggested to study the prepossessing methods in future work.
In addition to identifying buried targets, GPR can provide a good estimate of tar-
get depth, which may be critical to avoid accidental rupture of buried utilities, which
can cause serious damages.The accuracy of depth estimates using GPR depends on
precise knowledge of the velocity of electromagnetic wave in the medium.It’s also sug-
gested to study different approaches used to estimate depth of buried utilities from
the hyperbolic response of ground penetrating radar.

49
 Bibliography
[1] Christian U. Grosse and Masayasu Ohtsu Masayasu Ohtsu, eds. Acoustic Emis-
sion Testing. Springer-Verlag Berlin Heidelberg, 2008.
[2] Charlie Chong/ Fion Zhang. Nondestructive Testing Handbook Electromagnetic
Testing. Ed. by Judy Corbett 2017. Vol. 5. 2004.
[3] N. T. Thanh, H. Sahli, and D. N. Hao. “Infrared Thermography for Buried Land-
mine Detection: Inverse Problem Setting”. In: IEEE Transactions on Geoscience
and Remote Sensing 46.12 (Dec. 2008), pp. 3987–4004. issn: 0196-2892.
[4] J. R. Wait. “The Electromagnetic Basis for Nondestructive Testing of Cylindrical
Conductors”. In: IEEE Transactions on Instrumentation and Measurement 27.3
(Sept. 1978), pp. 235–238. issn: 0018-9456.
[5] E. Pasolli et al. “Automatic Detection and Classification of Buried Objects
in GPR Images Using Genetic Algorithms and Support Vector Machines”. In:
IGARSS 2008 - 2008 IEEE International Geoscience and Remote Sensing Sym-
posium. Vol. 2. July 2008, pp. II-525-II-528. doi: 10 . 1109 / IGARSS . 2008 .
4779044.
[6] P. Gamba and S. Lossani. “Neural detection of pipe signatures in ground pene-
trating radar images”. In: IEEE Transactions on Geoscience and Remote Sensing
38 (Mar. 2000), pp. 790–797. doi: 10.1109/36.842008.
[7] S. Delbo, P. Gamba, and D. Roccato. “A fuzzy shell clustering approach to
recognize hyperbolic signatures in subsurface radar images”. In: IEEE Transac-
tions on Geoscience and Remote Sensing 38.3 (May 2000), pp. 1447–1451. issn:
0196-2892. doi: 10.1109/36.843039.
[8] J.N. Kapur, P.K. Sahoo, and A.K.C. Wong. “A new method for gray-level pic-
ture thresholding using the entropy of the histogram”. In: Computer Vision,
Graphics, and Image Processing 29.3 (1985), pp. 273–285. issn: 0734-189X. doi:
https : / / doi . org / 10 . 1016 / 0734 - 189X(85 ) 90125 - 2. url: http : / / www .
sciencedirect.com/science/article/pii/0734189X85901252.
[9] G L Wang et al. “3D near-to-surface conductivity reconstruction by inversion of
VETEM data using the distorted Born iterative method”. In: Inverse Problems
20.6 (2004), S195.
[10] Q. Dong and C. M. Rappaport. “Microwave Subsurface Imaging Using Direct
Finite-Difference Frequency-Domain-Based Inversion”. In: IEEE Transactions
on Geoscience and Remote Sensing 47.11 (Nov. 2009), pp. 3664–3670. issn:
0196-2892. doi: 10.1109/TGRS.2009.2028740.
[11] Arthur B Weglein et al. “Inverse scattering series and seismic exploration”. In:
Inverse Problems 19.6 (2003), R27.

50
[12] C. Gilmore, I. Jeffrey, and J. LoVetri. “Derivation and comparison of SAR and
frequency-wavenumber migration within a common inverse scalar wave problem
formulation”. In: IEEE Transactions on Geoscience and Remote Sensing 44.6
(June 2006), pp. 1454–1461. issn: 0196-2892. doi: 10.1109/TGRS.2006.870402.
[13] Chao Ma, Lin Zhang, and Hu Zheng. “Extended born approximation technique
applied to linear inversion method”. In: 2008 World Automation Congress. Sept.
2008, pp. 1–3.
[14] A. P. Annan. “GPR—History, Trends, and Future Developments”. In: Springer
Nature,Subsurface Sensing Technologies and Applications: An International Jour-
nal 12 (Jan. 2002), pp. 2–13.
[15] Jonathon Ralston and Chad Hargrave. Forward Looking Radar for Early Sensing
of Ground Conditions. Tech. rep. 2006.
[16] D. J. Daniels, D. J. Gunton, and H. F. Scott. “Introduction to subsurface radar”.
In: IEEE Proceedings 135.54 (1988), pp. 278–320.
[17] David J. Daniels. Ground Penetrating Radar. 2nd Ed. 2004.
[18] David J. Daniels. Ground Penetrating Radar Theory and Application. Ed. by
Harry M.Jol. Elvester B.V, 2009.
[19] Peters. Jr, J. J. Daniels, and J. D. Young. “Ground penetrating radar as a sub-
surface environmental sensing tool”. In: Proceedings of the IEEE 82.12 (1994),
pp. 1802–1822.
[20] Jeffrey J. Daniels. “Ground Penetrating Radar Fundamentals”. In: 2000.
[21] Ahmet Serdar Turk et al. “High-resolution signal processing techniques for through-
the-wall imaging radar systems”. In: International Journal of Microwave and
Wireless Technologies 8.6 (2016), pp. 855–863. doi: 10.1017/S1759078716000593.
[22] S. U. Pillai and B. H. Kwon. “Forward/backward spatial smoothing techniques
for coherent signal identification”. In: IEEE Transactions on Acoustics, Speech,
and Signal Processing 37.1 (Jan. 1989), pp. 8–15. issn: 0096-3518. doi: 10.1109/
29.17496.
[23] Sathish Chandran. Advances in Direction-of-Arrival Estimation. Ed. by Sathish
Chandran. Artech House, December 31, 2005.
[24] Michel Bouvet, Bienvenu, and G. Georges. High-resolution methods in underwa-
ter acoustics. Berlin, New York, Springer-Verlag, 1991.
[25] S. Meschino and L. Pajewski. “Application of a SAP-DoA method to GPR data
for the location of reinforcing elements in concrete”. In: 2015 IEEE 15th Mediter-
ranean Microwave Symposium (MMS). Nov. 2015, pp. 1–4. doi: 10.1109/MMS.
2015.7375408.
[26] Giuseppe Schettini Simone Meschino Lara Pajewski. “Use of a sub-array statis-
tical approach for the detection of a buried object”. In: Near Surface Geophysics
8.5 (2010), pp. 365–375.
[27] S. Meschino et al. “Sub-surface localization of buried cylinders by means of sub-
array processing”. In: 2012 6th European Conference on Antennas and Propaga-
tion (EUCAP). Mar. 2012, pp. 1373–1377. doi: 10.1109/EuCAP.2012.6206193.

51
[28] R. Kumaresan and D. W. Tufts. “Estimating the Angles of Arrival of Multiple
Plane Waves”. In: IEEE Transactions on Aerospace and Electronic Systems AES-
19.1 (Jan. 1983), pp. 134–139. issn: 0018-9251. doi: 10 . 1109 / TAES . 1983 .
309427.
[29] P. Stoica and R. L. Moses. SPECTRAL ANALYSIS OF SIGNALS. Prentice-Hall
Upper Saddle River NJ, 2005.
[30] H. L. Van Trees. Optimum Array Processing Part IV of Detection and Estimation
and Modulation Theory. Wiley New York, 2002.
[31] P. Stoica, P. Babu, and J. Li. “New Method of Sparse Parameter Estimation
in Separable Models and Its Use for Spectral Analysis of Irregularly Sampled
Data”. In: IEEE Transactions on Signal Processing 59.1 (Jan. 2011), pp. 35–47.
issn: 1053-587X. doi: 10.1109/TSP.2010.2086452.
[32] Q. Zhang et al. “Fast implementation of sparse iterative covariance-based es-
timation for array processing”. In: 2011 Conference Record of the Forty Fifth
Asilomar Conference on Signals, Systems and Computers (ASILOMAR). Nov.
2011, pp. 2031–2035. doi: 10.1109/ACSSC.2011.6190383.
[33] P. Gupta and S. P. Kar. “MUSIC and improved MUSIC algorithm to estimate
direction of arrival”. In: 2015 International Conference on Communications and
Signal Processing (ICCSP). Apr. 2015, pp. 0757–0761. doi: 10.1109/ICCSP.
2015.7322593.
[34] R. Schmidt. “Multiple emitter location and signal parameter estimation”. In:
IEEE Transactions on Antennas and Propagation 34.3 (Mar. 1986), pp. 276–
280. issn: 0018-926X. doi: 10.1109/TAP.1986.1143830.
[35] F. Wen et al. “Improved MUSIC Algorithm for Multiple Noncoherent Subar-
rays”. In: IEEE Signal Processing Letters 21.5 (May 2014), pp. 527–530. issn:
1070-9908. doi: 10.1109/LSP.2014.2308271.
[36] A. Giannopoulos. “Modelling ground penetrating radar by GprMax”. In: Con-
struction and Building Materials 19.10 (2005), pp. 755–762.
[37] K. Yee. “Numerical solution of inital boundary value problems involving maxwell’s
equations in isotropic media”. In: IEEE Transactions on Antennas and Propa-
gation 14 (May 1966), pp. 302–307. doi: 10.1109/TAP.1966.1138693.
[38] A. Giannopoulos. “Unsplit Implementation of Higher Order PMLs”. In: IEEE
Transactions on Antennas and Propagation 60.3 (Mar. 2012), pp. 1479–1485.
issn: 0018-926X. doi: 10.1109/TAP.2011.2180344.
[39] R. J. Luebbers, K. S. Kunz, and K. A. Chamberlin. “An interactive demonstra-
tion of electromagnetic wave propagation using time-domain finite differences”.
In: IEEE Transactions on Education 33.1 (Feb. 1990), pp. 60–68. issn: 0018-
9359. doi: 10.1109/13.53628.
[40] A. Giannopoulos, C. Warren, and I. Giannakis. “An advanced GPR modelling
framework: The next generation of gprMax”. In: 2015 8th International Work-
shop on Advanced Ground Penetrating Radar (IWAGPR). July 2015, pp. 1–4.
doi: 10.1109/IWAGPR.2015.7292621.

52
[41] N. R. Peplinski, F. T. Ulaby, and M. C. Dobson. “Dielectric properties of soils
in the 0.3-1.3-GHz range”. In: IEEE Transactions on Geoscience and Remote
Sensing 33.3 (May 1995), pp. 803–807. issn: 0196-2892. doi: 10 . 1109 / 36 .
387598.
[42] M. C. Dobson et al. “Microwave Dielectric Behavior of Wet Soil-Part II: Dielec-
tric Mixing Models”. In: IEEE Transactions on Geoscience and Remote Sensing
GE-23.1 (Jan. 1985), pp. 35–46. issn: 0196-2892. doi: 10.1109/TGRS.1985.
289498.
[43] M. Pieraccini et al. “Propagation of Large Bandwidth Microwave Signals in
Water”. In: IEEE Transactions on Antennas and Propagation 57.11 (Nov. 2009),
pp. 3612–3618. issn: 0018-926X. doi: 10.1109/TAP.2009.2025674.
[44] I. Giannakis, A. Giannopoulos, and N. Davidson. “Incorporating dispersive elec-
trical properties in FDTD GPR models using a general Cole-Cole dispersion
function”. In: 2012 14th International Conference on Ground Penetrating Radar
(GPR). June 2012, pp. 232–236. doi: 10.1109/ICGPR.2012.6254866.

53
 Appendices
First Scenario(2D Model with single cylinder )
# title: Model of Cylinder in concrete
# material: 7 0.05 3 0 concrete
# domain: 0.3 0.002 0.2
# dx dy dz: 0.002 0.002 0.002
# time window: 8.0e− 9
# box: 0 0 0 0.3 0.002 0.125 concrete
# cylinder: 0.15 0 0.0625 0.15 0.002 0.0625 0.010 pec
# waveform: gaussiandot 1 1.5e9 my gaussian
# hertzian dipole: y 0.022 0 0.135 my gaussian
# rx array: 0.03 0 0.135 0.034 0 0.135 0 0 0
# src steps: 0.002 0 0
# rx steps: 0.002 0 0
# geometry view: 0 0 0 0.3 0.002 0.2 0.002 0.002 0.002 Bscan Cylinder in concrete n

Second Scenario(2D Model with 3 cylinder)

#title: Model of Cylinder in concrete
#material: 7 0.05 3 0 concrete
#domain: 1 0.002 0.4
#dx dy dz: 0.003 0.003 0.003
#time window: 8.0e − 9
#box: 0 0 0 1 0.002 0.3 concrete
#cylinder: 0.15 0 0.2 0.15 0.002 0.2 0.010 pec
#cylinder: 0.5 0 0.2 0.5 0.002 0.2 0.010 pec
#cylinder: 0.85 0 0.2 0.85 0.002 0.2 0.010 pec
#waveform: gaussiandot 1 1.5e9 my gaussian
#hertzian dipole: y 0.04 0 0.32 my gaussian
#rx array: 0.05 0 0.32 0.058 0 0.32 0 0 0
#src steps: 0.002 0 0
#rx steps: 0.002 0 0
#geometry view: 0 0 0 1 0.002 0.4 0.003 0.003 0.003 Bscan3 Cylinder in concrete n
3D Model with single cylinder in heterogeneous environment
#domain: 0.4 0.05 0.1
#dx dy dz: 0.002 0.002 0.002
#time window: 10e-9
#waveform: ricker 1 1.5e9 my ricker
#hertzian dipole: y 0.025 0.025 0.08 my ricker
#rx array: 0.035 0.022 0.08 0.035 0.026 0.08 0 0 0
#rx steps: 0.002 0 0
#soil peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 soil properties
#fractal box: 0 0 0 0.4 0.05 0.075 1.5 1 1 1 10 soil properties soil #add surface roughness:
0 0 0.075 0.4 0.05 0.075 1 1 1 0.05 0.0525 soil
#cylinder: 0.2 0 0.0375 0.2 0.05 0.0375 0.0150 pec
#geometry view: 0 0 0 0.4 0.05 0.1 0.002 0.002 0.002 CylinderInCompEnv n

3D Model with single cylinder in very complex environment

54
− − − − −− Define Model Domain
#domain: 0.2 0.2 0.2
− − − − −−Spatial discretisation Mostly recommended to use < 0.003 to minimize
error or increase efficiency
#dx dy dz: 0.002 0.002 0.002
− − − − −−Time step required for iteration
#time window: 10e− 9
−−−−−−Simple source with ricker wavefore i.e Mexican Hat waveform ,it is possible
to use one of gaussian,gaussiandot, gaussiandotnorm · · · etc built in wave form
#waveform: ricker 1 1.5e9 my ricker
− − − − −− Z oriented transmitter and it always must be located inside domain
#hertzian dipole: y 0.10 0.15 0.153 my ricker
− − − − −− Array Antenna with enter
element space 0.001,it must be located inside domain
#rx array: 0.11 0.15 0.153 0.114 0.15 0.153 0 0 0
− − − − −− Allows you to use a mixing model for soils proposed by Peplinski and
rough surface,water and grass added
#soil peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 soil properties
#fractal box: 0 0 0 0.2 0.2 0.15 1.5 1 1 1 10 soil properties soil
#add surface roughness: 0 0 0.15 0.2 0.2 0.15 1 1 1 0.1 0.15 soil
#add surface water: 0 0 0.15 0.2 0.2 0.15 0.125 soil
− − − − −− Buried Cylinder in complex environment
#cylinder: 0.1 0 0.075 0.1 0.20 0.075 0.010 pec
#add grass: 0 0 0.15 0.2 0.2 0.15 1 0.151 0.1515 5000 soil
− − − − −− python code to invoke antennas from built in librery
− − − − −− #python:
− − − − −−from user libs.antennas import antenna like GSSI 1500
− − − − −−antenna like GSSI 1500(0.5, 0.104, 0.6)
− − − − −−#end python:
− − − − −− provide the geometry output which can be build on visualization tool
kit(e.g.Paraview)
#geometry view: 0 0 0 0.2 0.2 0.2 0.002 0.002 0.002 CylinderInCompEnv n

55
 Table 6.1: 2D/3D GPR Modelling commands of gprMax
Command Function
#domain: Determine size of the model in the x, y, and z directions
#dx dy dz: Allows you to specify the discretization of space in the
x, y and z directions respectively
#time window: Allows you to specify the total required simulated time
#material: Allows you to introduce a material into the model de-
scribed by a set of constitutive parameters
#box: Allows you to introduce an orthogonal parallelepiped
with specific properties into the model
#add dispersion debye: Allows you to add dispersive properties to an already
defined material based on a multiple pole Debye formu-
lation
#add dispersion lorentz: Allows you to add dispersive properties to an already
defined material based on a multiple pole Lorentz for-
mulation
#add dispersion drude: Allows you to add dispersive properties to an already
defined material based on a multiple pole Drude formu-
lation
#soil peplinski: Allows you to use a mixing model for soils proposed by
Peplinski
#fractal box: Allows you to introduce an orthogonal parallelepiped
with fractal distributed properties which are related to
a mixing model or normal material into the model
#add surface roughness: Allows you to add rough surfaces to a fractal box in the
model.
#cylinder: Allows you to introduce a circular cylinder into the
model
#add surface water: Allows you to add surface water to a fractal box in the
model that has had a rough surface applied
#add grass: Allows you to add grass with roots to a fractal box in
the model
#hertzian dipole: Allows you to specify a current density term at an elec-
tric field location - the simplest excitation, often referred
to as an additive or soft source.
#rx: Allows you to introduce output points into the model.
These are locations where the values of the electric and
magnetic field components over the number of iterations
of the model will be saved to file
#rx array: Provides a simple method of defining multiple output
points in the model
#src steps and Provides a simple method to allow you to move the lo-
#rx steps: cation of all simple sources and receivers

56