ELECTROMAGNETIC MODELING OF SUBSURFACE LIGHT NONAQUEOUS PHASE-LIQUIDS SPILLS
By
Sandra Soto-Cabán

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Electrical Engineering

2011

ABSTRACT
ELECTROMAGNETIC MODELING OF SUBSURFACE LIGHT NON-AQUEOUS
PHASE-LIQUIDS SPILLS
By
Sandra Soto-Cabán
Hydrocarbon contamination from fuel leaks and spills has been one of the most common
ground water pollution problems in the United States for decades. Generally, underground leaks
and spills are difficult to detect and locate and can cause serious health problems to people living
close to the contaminated areas. The goal of this study is to simulate a contaminated region and
determine its physical limits within a tolerable level of certainty given the characteristics of the
soil in question without the environmental and financial cost of physical experimentation.
The first part of this study is focused in the determination of the dielectric properties of
soil and contaminated soil. A coaxial fixture was designed to characterize dry and contaminated
soil samples in the frequency range of 100 MHz to 1 GHz. Experimental results for contaminated
and uncontaminated soil samples are presented and analyzed. Values of permittivity obtained
using the designed coaxial waveguide were used to calculate the scattering electric fields of a
simulated contaminated underground system. The information collected from these simulations
was used in the reconstruction of the contrast profile of two lossless configurations involving
dielectric circular cylinders buried in soil. The image-reconstruction algorithm is based on the
Fourier diffraction theorem and multiple frequencies were used to improve resolution.
The problem was simulated using a two-dimensional formulation since three-dimensional
structures need not be considered in the context of the underground mapping application. This
study is not intended to achieve buried material identification of contaminants. However, a wide,
full volume scanning of a potentially contaminated area is useful for ecological impact

mitigation since it can be coupled with borehole assays to determine the specific contaminant
and extent of contamination. With such information, environmental engineers can determine the
need and the most suitable method for mitigation. Hence, the proposed radio frequency based
survey research is meant to provide a screening tool to be used in concert with other methods.

I dedicate this dissertation to my wonderful family, particularly to my understanding
and patient husband, Ferdinand, who has put up with these many years of research, and to our
children, Gabhriel, Fernando, and Frances who are the joy of our lives and my inspiration.

iv

ACKNOWLEDGEMENTS
My thanks and appreciation to Dr. Leo Kempel for persevering with me as my advisor
throughout the time it took me to complete this research and write the dissertation. Also, I would
like to thank the remaining members of my dissertation committee: Dr. Edward Rothwell for
always having the time and willingness to answer all my questions and review my work, Dr.
Dennis Nyquist for teaching me the theoretical foundations for my work and continue been part
of my dissertation committee after his retirement, and Dr. Antonello Tamburrino for his
disposition to provide help and support regardless of distance.
Very special thanks to Dr. Percy Pierre for his support throughout all these years. To Dr.
Barbara O’Kelly for her support, friendship, and having the time to edit and correct my work.
I would also like to thank Dr. Eric Law, Department of Geology at Muskingum
University, for providing unlimited access to his Soils Laboratory and equipment.
I would also like to thank my parents, Fidel and Cuca, for always supporting me and
encouraging me with their best wishes.
Finally, I would like to thank my wonderful husband, Ferdinand Avila-Medina. He was
always there cheering me up and stood by me through the good times and bad.

v

TABLE OF CONTENTS

LIST OF FIGURES .................................................................................................................... VIII
LIST OF TABLES ...................................................................................................................... XIII
KEY TO SYMBOLS AND ABBREVIATIONS ...................................................................... XIV
CHAPTER 1
INTRODUCTION ...........................................................................................................................1

CHAPTER 2
COAXIAL WAVEGUIDE FIXTURE FOR THE MEASUREMENT OF DIELECTRIC
PROPERTIES OF CONTAMINATED SOIL AT UHF FREQUENCIES .....................................5
2.1 Introduction ..................................................................................................................5
2.2 Transmission/Reflection method ..................................................................................7
2.3 Measurement procedure .............................................................................................11
2.3.1 Design of the coaxial waveguide fixture ........................................................11
2.3.2 Experimental setup .........................................................................................13
2.4 Sample preparation .....................................................................................................19
2.5 Calculation of soil electromagnetic properties ...........................................................21
2.6 Experimental results ...................................................................................................24
CHAPTER 3
ANALYTICAL FORMULATION OF ELECTROMAGNETIC SCATTERING FROM
AN UNDERGROUND HYDROCARBON SPILL ......................................................................33
3.1 Introduction ................................................................................................................33
3.2 Scattering by dielectric objects...................................................................................36
3.2.1 Normal incidence plane wave scattering by dielectric circular
cylinder: fundamental theory..........................................................................36
3.2.2 Layered dielectric cylinders: theory and analytical solutions ........................52
3.3 Simulation of electromagnetic scattering from an underground hydrocarbon
spill ........................................................................................................................62
3.3.1 Homogeneous spills illuminated by a plane wave .........................................63
3.3.2 Inhomogeneous spill illuminated by a plane wave ........................................68
3.3.3 Line source illumination .................................................................................71
CHAPTER 4
INVERSION ALGORITHM .........................................................................................................74
4.1 Introduction ................................................................................................................74
4.2 Diffraction Tomography.............................................................................................74
4.3 Fourier Diffraction Theorem ......................................................................................76
4.4 Reconstruction Algorithm ..........................................................................................85
4.5 Results ........................................................................................................................89

vi

CHAPTER 5
CONCLUSIONS AND FUTURE WORK ....................................................................................95
BIBLIOGRAPHY ..........................................................................................................................99

vii

LIST OF FIGURES
Figure 1-1. Cross-sectional view of a hypothetical LNAPL spill. Figure has been adapted
from Delin and others, 1998, USGS Fact Sheet FS-084-98. For
interpretation of the references to color in this and all other figures, the
reader is referred to the electronic version of this dissertation..................................2
Figure 1-2. A schematic of a typical leaking underground storage tank (LUST) spill site.
The leaky tank releases gasoline, or LNAPL. The gasoline descends
through the unsaturated soil zone to float on the water table (gasoline is
lighter than water) causing contamination of the aquifer [5]. ...................................2
Figure 2-1.

Coaxial sample holder. ...............................................................................................8

Figure 2-2.

Cross-sectional view of the coaxial sample holder. ...................................................8

Figure 2-3. Dimensions of the coaxial waveguide sample holder. The inner conductor is
placed inside the outer conductor for measurements as shown in the top
view. ........................................................................................................................13
Figure 2-4.

ASTM D 4935 Flanged Coaxial Fixture ..................................................................14

Figure 2-5. Measurement setup. The material sample can be seen in the interior of the
sample holder. .........................................................................................................14
Figure 2-6.

Calibration standards. ...............................................................................................15

Figure 2-7. Complex relative permittivity and permeability of air in a frequency range


from 100 to 950 MHz. a) Relative permittivity  r   r  j r  . b)



Relative permeability  r  r  j r  . .................................................................16

Figure 2-8. Complex relative permittivity and permeability of 9.50 cm long Plexiglas
sample in a frequency range from 100 to 950 MHz. a) Relative




permittivity  r   r  j r  . b) Relative permeability  r   r  j r  . ...............17
Figure 2-9. Complex relative permittivity and permeability of 4.75 cm long Plexiglas
sample in a frequency range from 100 to 900 MHz. a) Relative




permittivity  r   r  j r  . b) Relative permeability  r  r  j r  . ...............18
Figure 2-10. The red star indicates the area where the soil was collected. The map was
generated using the Web Soil Survey tool of the USDA Natural Resources
and Conservation Service - Soil Survey Area Map for Muskingum County
OH (Survey OH119) [23]. .......................................................................................19

viii

Figure 2-11. Complex relative permittivity of dry soil sample. Predicted values of
permittivity were obtained using the dielectric model for soils developed by
Peplinski et al. [26]. ................................................................................................25
Figure 2-12. Dielectric spectra of bulk/liquid water calculated using the Debye model and
measured data. Relaxation times of 9.23, 8.30, 6.48, and 5.17 ps were used
for liquid water at T = 20, 25, 35, and 45°C, respectively [36]. .............................26
Figure 2-13. Complex relative permittivity of soil with 10% moisture content. Predicted
values of permittivity were obtained using the dielectric model for soils
developed by Peplinski et al. [26]. ..........................................................................27
Figure 2-14. Complex relative permittivity of soil with 20% moisture content. Predicted
values of permittivity were obtained using the dielectric model for soils
developed by Peplinski et al. [26] ...........................................................................28
Figure 2-15. Complex relative permittivity of soil with 25% moisture content. Predicted
values of permittivity were obtained using the dielectric model for soils
developed by Peplinski et al. [26] ...........................................................................28
Figure 2-16. Real relative permittivity value for dry soil with 25% oil content. Predicted
values of permittivity were obtained using the Maxwell Garnett mixing
formula presented in [29]. .......................................................................................30
Figure 2-17. Results for a mixture of 25% volume of contaminant in soil with 10% and
25% volumetric moisture content. Measured results are compared with the
predicted values calculated using the mixing rules for moist and high-loss
materials presented in [41]. .....................................................................................31
Figure 2-18. Power attenuation in dry, wet and contaminated soil in a frequency range of
100 MHz to 1 GHz. .................................................................................................32
Figure 3-1. Schematic diagram of a possible underground spill and borehole logging
system. For simplicity, the soil is modeled as a homogeneous dielectric
space and the contaminated area is modeled as two eccentric, homogeneous
circular cylinders. ....................................................................................................34
Figure 3-2. Cross-sectional view of the spill illustrated in Figure 3-1. The spill is modeled
as an eccentric dielectric cylinder embedded in another dielectric cylinder.
The non-contaminated soil is modeled as a homogeneous dielectric space ............35
Figure 3-3. Definitions of E-parallel (TM) and E-perpendicular (TE) polarizations
relative to a cylinder, as defined in [49]. TM polarization occurs when the
electric field is parallel to the long axis of the cylinder. TE polarization
occurs when the electric field is perpendicular to the long axis of the
cylinder. ...................................................................................................................37

ix

Figure 3-4. Magnitude of the bistatic scattering field for a dielectric cylinder with radius
/2, surrounded by wet soil. Different numbers of modes were used to
calculate the scattering field and analyze the convergence of the series
expansion. The series was truncated when the percent difference between
the results of consecutive modes was less than 0.01%............................................47
Figure 3-5. Backscattering widths as a function of radius, normalized by the wavelength
(λ) of the incident field. TM backscattering widths are greater than TE
backscattering widths for most cylinders. ...............................................................49
Figure 3-6. Scattering widths for high contrast dielectric cylinders, normalized by the
wavelength (λ) of the incident field. As the radius of the cylinder becomes
small compared to wavelength, TM scattering widths become nearly
constant as a function of scattering angle. ...............................................................50
Figure 3-7. Scattering width as a function of frequency for cylinders of different radii.
The radii are varied from 0.05 to  at the highest frequency. Forward (top)
and back (bottom) scattering widths are illustrated. ................................................51
Figure 3-8. Scattering width as a function of the contrast ratio  c . Similar response is
obtained for high- and low-impedance dielectric cylinders. ...................................52
Figure 3-9. Geometry of the layered dielectric cylinders. The TM z wave traveling in a
medium with permittivity  s and permeability 0 (Region 0) is incident
on a layered dielectric circular cylinder of outer radius b with permittivity
1 (Region I), and inner radius a with permittivity  M (Region II).......................53
Figure 3-10. Scattering by a homogeneous cylinder in air. (a) Scattering width  2 D / 
as a function of scattering angle . 0  180 ,  r  2.2  j 0 , and r = . (b)
The simulated electric field magnitude. The black line denotes the cylinder
location. ...................................................................................................................64
Figure 3-11. Scattering by a homogeneous oil spill in dry soil. (a) Scattering width

 2 D /  as a function of scattering angle . 0  180 ,  r  2.2  j 0 , and
r = . (b) The simulated electric field magnitude. The black line denotes the
cylinder location. .....................................................................................................65
Figure 3-12. Scattering by a homogeneous oil spill in wet soil. (a) Scattering width

 2 D /  as a function of scattering angle . 0  180 ,  r  2.2  j 0 , and
r = . (b) The simulated electric field magnitude. The black line denotes the
cylinder location. .....................................................................................................66
Figure 3-13. Scattering by a homogeneous wet-contaminated spill in wet soil. (a)
Scattering width  2 D /  as a function of scattering angle . 0  180 ,
x

 r  12  j 0 , and r = . (b) The simulated electric field magnitude. The
black line denotes the cylinder location. .................................................................67
Figure 3-14. Simulated electric field magnitude for a layered small spill in wet soil at 500
MHz. The inner circle represents pure contamination of oil (r = 2.2) and
the outer cylinder represents wet-contaminated soil (r = 12). The cylinders
are surrounded by wet soil.......................................................................................69
Figure 3-15. Simulated electric field magnitude for an eccentric spill in wet soil at 500
MHz. The inner circle represents pure contamination of oil (r = 2.2) and
the outer cylinder represents wet-contaminated soil (r = 12). The cylinders
are surrounded by wet soil.......................................................................................70
Figure 3-16. Electric field magnitude for a spill in a vadose region. The blue dashed line
represents the position of the spill. The line source is located at 3soil from
the spill. ...................................................................................................................71
Figure 3-17. Simulated electric field magnitude for a concentric spill in wet soil
illuminated by a line source with frequency of 500 MHz. The inner circle
represents pure contamination of oil (r = 2.2) and the outer cylinder
represents wet-contaminated soil (r = 12). .............................................................72
Figure 3-18. Simulated electric field magnitude for an eccentric spill in wet soil at 500
MHz. The inner circle represents pure contamination of oil (r = 2.2) and
the outer cylinder represents wet-contaminated soil (r = 12).................................73
Figure 4-1. Illustration of the Fourier Diffraction Theorem. When 2D object is
illuminated by a plane wave, the 2D Fourier transform (FT) of the scattered
field measured along the line TT' gives the values of the 2D Fourier
transform of the contrast evaluated along a circular arc centered in -k0 and
radius k0 in the frequency domain. ..........................................................................77
Figure 4-2. Fourier domain for an incident plane wave. Under Born approximation, the
2D Fourier transform of the contrast profile is obtained in the transformed
domain over a circumference shifted on the direction opposite to the
incident plane wave. ................................................................................................81
Figure 4-3. Multiview Fourier domain. Different directions of illumination are necessary
to obtain the reconstruction of the object. The 2D Fourier transform of the
contrast profile is obtained inside a circle of radius 2k0 centered at the
origin. ......................................................................................................................82
Figure 4-4. Circular array of antennas. Different directions of illumination can be
obtained by transmitting with each one of the elements. The transmitter is
Tx and all others are the receivers Rx. .....................................................................83

xi

Figure 4-5. Contrast profile for a centered spill configuration. The contrast of the
cylinders is constant and purely real. For the large cylinder the contrast is
0.294. For the small cylinder, or core region, the contrast is 0.8167. (a)
Original image. (b) Reconstructed profile with 16 antennas. (c)
Reconstructed profile with 32 antennas. (d) Reconstructed profile with 64
antennas. ..................................................................................................................90
Figure 4-6. Contrast profile for a centered spill configuration. For the large cylinder the
contrast is 0.294. For the small cylinder, or core region, the contrast is
0.8167. (a) Original image contrast profile. (b) Reconstructed profile with
16 antennas. (c) Reconstructed profile with 32 antennas. (d) Reconstructed
profile with 64 antennas. .........................................................................................91
Figure 4-7. Contrast profile for an eccentric spill configuration. The contrast of the
cylinders is constant and purely real. For the large cylinder the contrast is
0.294. For the small cylinder, or core region, the contrast is 0.8167. (a)
Original image. (b) Reconstructed profile with 16 antennas. (c)
Reconstructed profile with 32 antennas. (d) Reconstructed profile with 64
antennas. ..................................................................................................................92
Figure 4-8. Contrast profile for a centered spill configuration. The contrast of the
cylinders is constant and purely real. For the large cylinder the contrast is
0.294. For the small cylinder, or core region, the contrast is 0.8167. (a)
Original image. (b) Reconstructed profile with 16 antennas. (c)
Reconstructed profile with 32 antennas. (d) Reconstructed profile with 64
antennas. ..................................................................................................................93

xii

LIST OF TABLES
Table 2-1. Physical Soil Properties – Muskingum County, Ohio .................................................20
Table 2-2. Measured and calculated predicted values of permittivity of soil with different
moisture content in a frequency range from 300 to 800 MHz. Predicted
values of permittivity were obtained using the dielectric model for soils
developed by Peplinski et al. [26]...............................................................................27
Table 5-1. Percent difference for permittivity values of wet- and wet-contaminated soils. .........96

xiii

KEY TO SYMBOLS AND ABBREVIATIONS
DT

Diffraction Tomography

EM

Electromagnetic

GPR

Ground penetrating radar

LNAPL

Light non-aqueous phase liquids

LUST

Leaking underground storage tank

NIST

National Institute of Standards and Technology

RCS

Radar cross section

SAR

Synthetic Aperture Radar

SNR

Signal-to-Noise Ratio

SW

Scattering width

TE

Transverse electric

TEM

Transverse electromagnetic

TM

Transverse magnetic

TR

Transmission/reflection

UHF

Ultra-high frequency

VNA

Vector network analyzer

xiv

CHAPTER 1
INTRODUCTION

Hydrocarbon contamination from fuel leaks and spills has been one of the most common
ground water pollution problems in the United States for decades [1]. In particular, leaks from
underground storage tanks and improper disposal of gasoline and oil threaten drinking water
wells, lakes, streams, and basements all over the country [2]. Hydrocarbon contamination can be
widely dispersed or confined to a local region and they can contain a variety of chemicals that
can cause serious health problems to people living close to the contaminated areas. Leaks in
underground storage tanks and underground fuel transfer lines are difficult to detect and locate.
Often the leaks are so small that routine tests are ineffective in locating the source of the
pollutant. The rate of dissipation of hydrocarbon contamination in soil is extremely slow in many
cases and even a small fuel leak can produce serious contamination that is extremely costly to
remediate [3].
Most hydrocarbons are lighter than water with a low dielectric permittivity. Because of
these properties, they are known as light non-aqueous phase liquids (LNAPL). At LNAPL
contamination sites, LNAPL can form a pool in the subsurface on top of the water table [4]. The
following figures illustrate two examples of hydrocarbon contamination. Figure 1 is a cross
sectional view of a hypothetical LNAPL spill due to improper disposal of contaminants and
Figure 2 is an example of underground water contamination due to a leaking storage tank.

1

Figure 1-1. Cross-sectional view of a hypothetical LNAPL spill. Figure has been adapted from
Delin and others, 1998, USGS Fact Sheet FS-084-98. For interpretation of the
references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation.

Figure 1-2. A schematic of a typical leaking underground storage tank (LUST) spill site. The
leaky tank releases gasoline, or LNAPL. The gasoline descends through the
unsaturated soil zone to float on the water table (gasoline is lighter than water)
causing contamination of the aquifer [5].

2

In view of the magnitude and the continuing nature of this problem, it is important that
techniques for detecting and mapping hydrocarbon contamination in soil and ground water be
continually improved. There are various methods used to detect organic contaminants in soil.
The most common are chromatography and mass spectrometry [6]. These methods are based on
chemical analysis of the soil pore-fluid and are both time consuming and expensive. As an
alternative to chemical analysis, geophysical methods, such as surface geophysics and borehole
geophysical logging, can be used to aid in characterizing changes in subsurface features [7].
These methods are based on the propagation of electromagnetic waves and their performance is
determined fundamentally by the soil electromagnetic properties and the target characteristics.
Because of that, before using any of the geophysical methods available, a better knowledge of
the dielectric properties of the media involved in the investigation is needed.
The goal of this study is to simulate a contaminated region and have an algorithm that can
be used to determine its physical limits within a tolerable level of certainty given the
characteristics of the soil in question. The first part of this study is focused in the determination
of the dielectric properties of soil and contaminated soil. A durable and cost effective coaxial
waveguide was designed to obtain permittivity values of the soils. Dielectric measurements of
dry, wet, and contaminated soil were performed in a frequency range of 100 to 1000 MHz. The
design of the coaxial waveguide and measurement results are presented in Chapter 2. Values of
permittivity obtained using the designed coaxial waveguide were used for a computer simulation
of a contaminated underground system. The soil was modeled as a homogeneous dielectric space
and the contaminated area was modeled as an inhomogeneous dielectric space. The contaminated
area was assumed to be composed of wet and contaminated soil. Analytical formulation of the
electromagnetic scattering from the contaminated system is presented in Chapter 3. Different

3

system configurations were simulated and the information collected from these simulations was
used in the reconstruction of the contrast profile of two lossless configurations involving
dielectric circular cylinders buried in soil. The image-reconstruction algorithm is based on the
Fourier diffraction theorem [8] and multiple frequencies were used to improve resolution [9].
The inversion algorithm and reconstructed contrast images are presented in Chapter 4. The
concluding remarks and the future work are issued in Chapter 5.
This study is not intended for a direct identification of contaminants. However, a wide,
full volume scanning of a potentially contaminated area is useful for ecological impact
mitigation since it can be coupled with borehole assays to determine the specific contaminant
and extent of contamination. With such information, environmental engineers can determine the
need and the most suitable method of mitigation. Hence, the proposed radio frequency based
survey research is meant to provide a screening tool to be used in concert with other methods.

4

CHAPTER 2
COAXIAL WAVEGUIDE FIXTURE FOR THE MEASUREMENT OF DIELECTRIC
PROPERTIES OF CONTAMINATED SOIL AT UHF FREQUENCIES

2.1

Introduction
Identification of organic contaminants in the subsurface can be obtained using

geophysical methods such as ground penetrating radar (GPR) [10] and borehole logging [11].
The performance of these microwave technologies is determined fundamentally by the soil
electromagnetic (EM) properties and the target characteristics at study sites as well as the sensor
equipment. From the perspective of the signal, both wave speed and power attenuation depends
on the dielectric properties of the propagation media.
Many researchers have collected data on soil dielectric properties using different
measurement techniques. Measurement of dielectric properties involves measurements of the
complex relative permittivity (r) and complex relative permeability (μr) of the materials. There
are many techniques developed for measuring the complex permittivity and permeability and
each technique is limited to specific frequencies, materials, applications, etc. by its own
constraint. Some measurement techniques include open-ended coaxial probe, free-space,
resonant, and transmission/reflection line. Baker-Jarvis et al. present a discussion of these prior
efforts in [12]. Curtis [13] used a coaxial transmission/reflection measurement system to collect
complex dielectric property data for soils at 45 MHz to 26.5 GHz. Different from his work, this
research presents a new, large-volume coaxial sample holder, designed and fabricated at
Michigan State University for the dielectric measurements of soil and contaminated soil at the

5

UHF band. This fixture was designed and used since it offered wide bandwidth and the ability to
measure a volume of material that is sufficiently large at the operating frequencies so that
homogenization of properties is assured given the dimensions of the potential solid phase of the
samples.
The transmission/reflection (TR) method is widely used for the measurement of solid
materials. In the TR method, a material sample is placed in a section of a waveguide or coaxial
line after calibration of the fixture has been performed, and an incident electromagnetic field is
applied. The technique involves measurement of the reflected (S11) and transmitted signal (S21).
Scattering equations are found from the analysis of the electric field at the sample interface. The
scattering equations relate the measured scattering parameters to the permittivity ( ) and
permeability ( ) of the material through an inversion process. In developing the scattering
equations, only the fundamental waveguide mode is assumed to propagate.
The TR inversion method was first introduced by Nicholson and Ross [14] in the early
1970’s. In their paper, a new method for obtaining the complex permittivity and permeability of
linear materials over a broad range of frequencies is presented. Before this study, such
measurements were made at fixed frequencies using slotted-line or impedance bridge
configurations. After Nicholson and Ross, Weir [15] determined the complex values of  and 
from measurements made directly in the frequency domain. This is known as the NicholsonRoss-Weir (NRW) method. This method is a commonly used technique for direct calculation of
both the permittivity and permeability from the S-parameters. However, the technique diverges
for low-loss materials at frequencies corresponding to integer multiples of one-half wavelength
in the sample. Baker-Jarvis et al. [16] improved the NRW method by allowing measurements to
be taken on samples of arbitrary length. This technique is known as NIST Iterative method.

6

NIST Iterative technique performs the calculation using a Newton-Raphson’s root finding
technique and is suitable for permittivity calculation only. It works well if a good initial guess is
available and is suitable for long samples and characterizing low loss materials. Boughriet, et al.
[17] developed a non-iterative transmission/reflection method applicable to permittivity
measurements using arbitrary sample lengths in wide-band frequencies. This method is based on
a simplified version of the NRW method and its accuracies are comparable to the iterative
technique. It does not need an initial estimation of permittivity and can perform the calculation
very fast. Later, Courtney [18] used the Nicholson and Ross method, and extended the procedure
to account for the infinite set of transmissions and reflections of the transient wave at the
interfaces of the material sample.
In this chapter, the design of a new, large-volume coaxial sample holder for the dielectric
measurements of soil and contaminated soil at the UHF band is presented. Unlike traditional
sample holders common in many electromagnetic material characterization laboratories, it can
accommodate significant volumetric samples of material. The sample holder is a symmetric
coaxial waveguide fabricated using brass. It is intended to operate in a frequency range between
100 MHz and 1 GHz. Using a Vector Network Analyzer (VNA), the values for complex relative





permittivity ( r   r  j r ) and complex relative permeability (r  r  j r ) were calculated
using the TR method.
2.2

Transmission/Reflection method
A sample of material with unknown electromagnetic parameters  r , r , and known

thickness d is placed in an air-filled coaxial sample holder. A typical cylindrical coaxial sample
holder is illustrated in Figure 2-1.

7

Figure 2-1. Coaxial sample holder.

Figure 2-2 shows a cross-sectional view of the sample holder with the material-under-test
placed at the center of the holder in the volume between the inner and outer conductors. The
calibration reference plane positions for the measurements are shown and the electric fields in
each region of the line are indicated.
A
Einc

Z0

B

Z

Z0

MUT

Eref

L1

Etrans

d

L2

Reference planes
Figure 2-2. Cross-sectional view of the coaxial sample holder.

8

A transverse electromagnetic (TEM) wave propagates on the coaxial transmission line.



Measurements of the S-parameters, indicated as S11 and S21 , are made at the reference planes.
These values can be written for a time-harmonic excitation as:
Eref
Einc

(2.1)

E

S21    tran
Einc

(2.2)


S11   

where Einc is the incident field, Eref is the total reflected field, and Etran is the total transmitted
field. Translating these parameters to the interfaces between regions yields


S11    S11   e j 20 L1


S21    S21   e

(2.3)

j  0 L1  0 L2 

(2.4)

where 0   c0 ,   2 f , f = frequency, c0 = speed of light in vacuum, and L1 and L2 are the
distances between the sample and the reference plane positions. The complex transient reflection
coefficient  can be expressed as:

12 

Z  Z0

Z  Z0

r 2
r 2

 r 2 1
r2 1





where Z0  0  0 , Z  2  2 , 2  0  r 2  j r 2  ,  2   0  r 2  j r 2  .
The transmission coefficient through the sample material can be written as:

9

(2.5)

z  e j  2 d

(2.6)

where 2    r 2  r 2 c0 .
Equations (2.3) and (2.4) can be rewritten in terms of the wave impedances as:

S11   

 Z 2  Z02  1  e j2 d 
2

 Z  Z0 

2

2
  Z  Z0  e j 2 2d

(2.7)

and
4  Z0 Z  e j  2 d
.
S21   
 Z  Z 0 2   Z  Z 0 2 e  j  2 d

(2.8)

Using transmission-line theory [14]:





1  z 2 12
EA
S11   

Einc 1   2 z 2
12



(2.9)



2
1   1    z  1  12 z
EB
S21   

2
2
Einc
1  12 z 2
1  12 z 2

(2.10)

where EA and EB are the fields at planes A and B, respectively.
From the above equations, the complex permittivity and complex permeability can be
expressed as:

 1  12   j ln 1 z  c 


 1  12  2 fd 

r  

10

(2.11)

and

 1  12   j ln(1 z )c 

.
 1  12   2 fd 

r  

(2.12)

Using these equations, the complex permeability and permittivity can be obtained from
measurements of the transmission and reflection scattering coefficients of the material. For nonmagnetic materials, r = 1. This eliminates the ill-behaved term 1  12  1  12  , which
produces undesired ripples in the extracted r and r [17]. This technique gives smooth
permittivity results, with no divergences at frequencies that are multiples of one-half wavelength,
is accurate, arbitrary length of samples can be used, and no initial guess is needed. Some
disadvantages of the technique are that is only applicable for permittivity measurements and the
challenge of solving the logarithm of a complex function.
2.3

Measurement procedure

2.3.1 Design of the coaxial waveguide fixture
The designed sample holder is a symmetrical coaxial waveguide fixture made of brass (
7

= 1.57x10 S/m,  = 4x10

-7

H/m). It was designed to operate in a frequency range between

100 MHz to 1 GHz. The length of the coaxial waveguide is slightly less than a half-wavelength,
for the dominant TEM mode, of the signal at the desired higher frequency of 1 GHz. The
distributed electromagnetic parameters of the coaxial line fixture were determined using the
following equations [19]:

11

R

Rs  b 
1
2 b  a 



L

 b
ln  
2  a 

C

(2.13)

2 
b
ln  
a

G  C tan  R 

Rs  b 
1
2 b  a 



where R is the distributed resistance, L is the distributed inductance, C is the distributed
capacitance, and G is the distributed conductance of the line. The inner diameter of the outer
conductor is b and a is the outer diameter of the inner conductor. Rs is surface resistivity,   is
the real part of the permittivity of the dielectric material, and tanθ is the loss tangent of the
dielectric filling. These EM parameters can be used to determine the propagation constant , and
the characteristic impedance Z0 of the coaxial waveguide. For the designed fixture, L = 1.704 x
-7

10

H/m and C = 6.528x10

-11

F/m. The characteristic impedance is Z0  L C  51.08  . The

propagation constant  is j2.095 rad/m at 100 MHz, j10.477 rad/m at 500 MHz, and j20.954
rad/m at 1 GHz. The transverse electric mode TE11 has the lowest cutoff frequency of all the
higher modes. For this mode, the cutoff frequency was 1.785 GHz. The intended frequency range
of operation is below the cutoff condition and the transverse electromagnetic mode (TEM) will
be the only mode of propagation within the frequency range given sufficient stand-off distance
for the measurement ports and the calibration planes (e.g. the axial boundaries of the MUT). The
dimensions of the coaxial fixture are illustrated in Figure 2-3. The total volumetric capacity of
the sample holder is approximately 552 cm3. This volume permits homogenization of the

12

material properties even with relatively large (order of millimeter) solid phase constituent
materials in the soil mixture.

13.2 cm

3.3 cm

7.6 cm

3.3 cm
15 cm

7.6 cm
13.2 cm
Top view
Outer conductor

Inner conductor

Figure 2-3. Dimensions of the coaxial waveguide sample holder. The inner conductor is placed
inside the outer conductor for measurements as shown in the top view.

2.3.2 Experimental setup
The designed coaxial waveguide fixture was connected to an HP 8753D Vector Network
Analyzer (VNA) using a set of ASTM D 4935 flanged coaxial fixtures [20] shown in Figure 2-4.
The experimental setup is shown in Figure 2-5. The VNA was used for measurement of the Sparameters of the sample. The Thru-Reflection-Line (TRL) method was used for calibration
purposes [21]. The Thru standard was realized by directly connecting the two ASTM D 4935
coaxial sections. A circular brass plate was placed at the end of the coaxial section to provide a
reflection standard. The line standard was done using the empty sample holder fixture between
the two ASTM D 4935 coaxial sections. The calibration configurations are shown in Figure 2-6.

13

After the calibration process, the network analyzer reported the S-parameters of the sample,
which then were processed to obtain the required EM properties.

Figure 2-4. ASTM D 4935 Flanged Coaxial Fixture

Sample
holder

VNA

1

2
Material
sample

Figure 2-5. Measurement setup. The material sample can be seen in the interior of the sample
holder.

14

Section
with no
sample

Short

Thru

Reflection
Line
Figure 2-6. Calibration standards.

To test the system, air was used as the first sample material. Values of permittivity and
permeability for air are well known and well approximated by the free-space values of these
parameters. The measured complex permittivity and permeability for air are shown in Figure 2-7.
Measured values are in accordance with the expected values.
The second material used for testing was Plexiglas. Two Plexiglas samples with lengths
9.5 cm and 4.75 cm were machined to fit snugly inside the coaxial fixture. The samples were
placed as shown in Figure 2-5, at the bottom of the sample holder instead of the center. This
configuration allows for a flat surface when materials with loose particles, like soil, are used. For
the two Plexiglas samples, the average measured values of permittivity and permeability were

 r  2.58  j 0.02 and r  1.01  j 0.02 in a frequency range from 100 to 950 MHz. These

15

values are in accordance with the typical values for Plexiglas [22]. Measured values for the two
Plexiglas samples are presented in Figure 2-8 and Figure 2-9.

Figure 2-7. Complex relative permittivity and permeability of air in a frequency range from 100


to 950 MHz. a) Relative permittivity  r   r  j r  . b) Relative permeability



 r  r  jr  .

16

Figure 2-8. Complex relative permittivity and permeability of 9.50 cm long Plexiglas sample in


a frequency range from 100 to 950 MHz. a) Relative permittivity  r   r  j r  .



b) Relative permeability  r   r  j r  .

17

Figure 2-9. Complex relative permittivity and permeability of 4.75 cm long Plexiglas sample in


a frequency range from 100 to 900 MHz. a) Relative permittivity  r   r  j r  .



b) Relative permeability  r  r  j r  .

18

2.4

Sample preparation
Samples of dry soil, wet soil, and contaminated soil were prepared. The soil was

collected in Muskingum County, Ohio; specifically in the Zanesville region. Map of the region
with soil physical properties is shown in Figure 2-10. Description of the terms in each region of
the map and corresponding soil properties are presented in Table 2-1.

Figure 2-10. The red star indicates the area where the soil was collected. The map was generated
using the Web Soil Survey tool of the USDA Natural Resources and Conservation
Service - Soil Survey Area Map for Muskingum County OH (Survey OH119) [23].

19

Table 2-1. Physical Soil Properties – Muskingum County, Ohio
Map symbol and
soil name

--

8-20

--

--

12-26

10-50

--

--

22-32

--

--

8-20

--

--

15-23

--

--

18-30

20-40

--

--

24-41

--

--

24-41

--

--

--

--

--

12-17

3-44

5-15

50-83

12-35

44-60

15-35

30-78

7-35

0-11

--

--

7-27

11-32

8-20

50-74

18-35

32-60

10-20

50-78

12-40

--

--

--

--

0-9

--

--

15-30

9-40

10-35

30-70

20-35

40-60

10-35

30-72

18-35

60-62

--

--

--

0-7

--

42-95

12-27

7-28

5-20

45-70

18-35

28-45

5-20

47-77

18-33

45-80

10-35

25-70

20-40

0-7

--

42-95

12-27

7-28
ZnC2: Zanesville

--

0-3

ZnB: Zanesville

22-32

60-75

WuC2:
Westmoreland

--

40-60

Ud: Udorthents

--

8-20

Ne: Newark

10-50

0-8

Me: Melvin

Clay
(%)
12-26

50-80

CtE: Coshocton

Silt
(%)
--

0-10
AfC2: Alford

Sand
(%)
--

50-80

AfB: Alford

Depth
(inches)
0-10

5-20

45-70

18-35

28-45

7-42

40-60

18-33

45-80

0-38

42-64

20-40

20

Soil physical properties were measured and results compared with the data presented in
Table 2-1. The collected soil had a composition of approximately 15% clay, 5% sand and 80%
silt. Based on these results and Table 2-1, the soil is Zanesville silt loam (ZnC2) which agreed
with the soil in the collection region.
The collected soil was dried following Standard T 265 of the American Association of
State Highway and Transportation Officials (AASHTO) [24]. The soil was heated in an oven at
105°C for approximately 24 hours. The weight of the dry soil samples was measured, and a
predetermined amount of distilled water was then added to achieve the desired volumetric
moisture content. After weighing the wet soil samples, they were sealed and allowed to settle for
a few hours. From the weight of the dry soil Wd, and the weight of the water added Ww, the
gravimetric moisture content mg, can be calculated using mg  Ww Wd . The volumetric
moisture content mv, is determined by b mg , where b is the bulk density of the soil sample
given by Wd V , where V is the volume of the dry sample [25]. Since these values are bulk
densities, the possibility of air spaces between the soil particles is ignored.
2.5

Calculation of soil electromagnetic properties
Researchers have developed semi-empirical dielectric mixing models that calculate the

electrical properties of the different components of soil. For this study, a model for predicting



both the real and imaginary parts of the relative permittivity r ( r   r  j r ) for a known
frequency band and different moisture levels was required. Peplinski et al. [26] developed a
semi-empirical dielectric model for soils, covering the frequency range of 0.3 to 1.3 GHz. The
model provides frequency-dependent expressions for the real and imaginary parts of the relative

21

dielectric constants of a soil medium in terms of the fraction of sand particles S, the fraction of
clay particles C, the volumetric water content mv, and the bulk density of the soil b, all being
available parameters. The real part of the complex relative permittivity for the bulk soil is
approximated as:
1
 b 
    m    0.68 .

 m  1.15 1 
 s  mv fw
v
 s


 

(2.14)

The imaginary part is given as:
1
 m     

 m  v fw





(2.15)



where  m   m  j m is the relative complex dielectric constant of the soil-water mixture, s =
3

2.66 g/cm is the specific density of the solid soil particles,  = 0.65 is an empirically
determined constant [27], and  ' and  " are soil-type dependent constants given by

 '  1.2748  0.519S  0.152C

(2.16)

 "  1.33797  0.603S  0.166C

(2.17)

and

The quantities  fw and   are the real and imaginary parts of the relative dielectric constant
fw
of free water given by a Debye-type dispersion equation [28] are

22

 fw   w 

  
fw

 w0   w
1   2 f  w 

2

2 f  w  w0   w   eff   s  b 

2
2 0 f  s mv
1   2 f  w 

(2.18)

(2.19)

where  0 is the permittivity of free space, f is the frequency in hertz,  w is the relaxation time
for water,  w0 is the static dielectric constant for water, and  w  4.9 is the high-frequency
limit of  fw . At room temperature (20ºC), 2  w  0.58 1010 s and  w0  80.1 [28]. The
effective conductivity  eff derived in [26] is given by

 eff  0.0467  0.2204b  0.4111S  0.6614C

(2.20)

When a contaminant is added to the soil samples, a mixing formula is used to calculate
the effective permittivity of the mixture [29]. Mixing formulas or rules are used to provide an
approximate electromagnetic description for the structural complexity of materials occurring in
nature. Mixing formulas use the concept of effective or macroscopic permittivity, which implies
that the mixture responds to electromagnetic excitation as if it were homogeneous [30]. The size
of the inclusions in the mixture and the spatial correlation length of the permittivity function
need to be small with respect to the wavelength. There are many effective medium
approximations [29]-[31]. Perhaps the most common mixing rule is the Maxwell Garnett formula
[32]. It is simple in appearance but with broad applicability. The effective permittivity of a
mixture  eff relative to the environment permittivity  e is determined by only two parameters
and obeys the following rule:

23

i   e
 eff   e  3vi e
 i  2 e  vi  i   e 

(2.21)

where  i is the inclusion permittivity and vi is the volume fraction of the inclusions. Other
important mixing rules in the remote sensing community are the Polder–van Santen formula
[33], known as Bruggerman formula [34], and the Coherent Potential formula [35]. For dilute
mixtures (vi <<1), all three mixing rules, Maxwell Garnet, Polder–van Santen, and Coherent
Potential, predict the same results [30].
2.6

Experimental results
First, the permittivity of the dry soil sample was measured. The obtained value of

permittivity was compared to the predicted value using the dielectric model for soils developed
by Peplinski et al. [26]. Figure 2-11 illustrates the measured and predicted relative complex
permittivity of dry soil. Averaged measured relative permittivity in a frequency range from 300
to 800 MHz for dry soil is 2.0954-j0.0028. In the same frequency range, the calculated relative
permittivity using the Peplinski model is 2.0863-j0.000.

24

Figure 2-11. Complex relative permittivity of dry soil sample. Predicted values of permittivity
were obtained using the dielectric model for soils developed by Peplinski et al. [26].

Once we know the permittivity of dry soil, the second step is to measure the permittivity
of wet soil. First, the dielectric spectra for bulk/liquid water were calculated using equations 2.18
and 2.19 of the Debye model. The dielectric spectra of bulk/liquid water are presented in Figure
2-12. Relaxation times of 9.23, 8.30, 6.48, and 5.17 picoseconds were used for liquid water at a
temperature T = 20, 25, 35, and 45°C, respectively [36]. The obtained permittivity values using
the Debye model and measured data are in agreement with the literature, within 0.5% of
previously reported data [37].

25

Figure 2-12. Dielectric spectra of bulk/liquid water calculated using the Debye model and
measured data. Relaxation times of 9.23, 8.30, 6.48, and 5.17 ps were used for
liquid water at T = 20, 25, 35, and 45°C, respectively [36].

Wet soil samples were prepared by adding 10%, 20%, and 25% volume of moisture to
the dry soil samples. Averaged values of permittivity for the different wet samples are tabulated
in Table 2-2. Permittivity values versus frequency are presented in Figures 2-13, 2-14, and 2-15.
Table 2-2 As expected, a reasonable increase in permittivity values is obtained as the water
content of the sample increases. As shown in the following figures, measured results are in
agreement with the model predicted values.

26

Table 2-2. Measured and calculated predicted values of permittivity of soil with different
moisture content in a frequency range from 300 to 800 MHz. Predicted values of
permittivity were obtained using the dielectric model for soils developed by Peplinski
et al. [26].

Relative permittivity
(Average)

Moisture Content
10%

20%

25%

Measured

6.9148-j0.0028

13.1200-j0.0680

17.2000-j0.0255

Calculated

6.9648-j0.0606

13.2700-j0.2038

16.8700-j0.2748

Figure 2-13. Complex relative permittivity of soil with 10% moisture content. Predicted values
of permittivity were obtained using the dielectric model for soils developed by
Peplinski et al. [26].

27

Figure 2-14. Complex relative permittivity of soil with 20% moisture content. Predicted values
of permittivity were obtained using the dielectric model for soils developed by
Peplinski et al. [26]

Figure 2-15. Complex relative permittivity of soil with 25% moisture content. Predicted values
of permittivity were obtained using the dielectric model for soils developed by
Peplinski et al. [26]

28

Next the permittivity value for contaminated soil was determined. Dry soil was polluted
with a contaminant following the same procedure as with the wet soil samples. The contaminant
used was motor oil SAE 30. This material was chosen as representative of hydrocarbons without
having onerous material handling (e.g. safety) requirements. Dielectric properties of various
petroleum oils were given by von Hippel [22]. Based on his data, the dielectric constant and loss
tangent of motor oil used in this work were estimated using the following equations:


 r  2.24  0.000727T

(2.22)

tan     0.527T  4.82  104

(2.23)

where T is the temperature of the room in C and tan() is the loss tangent. For a room
temperature of 22C, the calculated relative permittivity of motor oil is  r  2.2236  j 0.0016 .
Its permittivity value is close to the relative permittivity of dry soil. Therefore, not much change
in permittivity was expected. Samples of dry soil with 25% volumetric oil content were
measured. The results were compared with the predicted values calculated using the Maxwell
Garnett mixing rule presented in [29]. Figure 2-16 shows the permittivity values obtained. There
was insufficient difference in permittivity values to distinguish between dry and contaminated
soil.

29

Figure 2-16. Real relative permittivity value for dry soil with 25% oil content. Predicted values
of permittivity were obtained using the Maxwell Garnett mixing formula presented
in [29].

Greater contrast should be observed, though, by considering soil with a higher water
content per unit volume. This represents soil close to groundwater. General facts about
groundwater can be found in [38]. The zone of soil above the water table is known as the
unsaturated zone or the vadose zone. Soil particles in the unsaturated zone are coated with layers
of water. No matter how tightly the soil particles are packed, there will always be void spaces
between them. The void spaces between the soil particles are known as the soil pores. Below the
water table the pore spaces are filled with water. Above the water table the pore spaces are filled
with varied amounts of air and water. As contaminants move vertically through water-wet soil to
the water table, it displaces the air and some of the water from the soil pores [39]. Soil particles
remain wetted by water and the contaminant emerges as a dispersed phase. At a macro-scale, the
oil does not appear to separate from the water due to the soil. After reaching the water table, the
dissolved contaminant will move along with the groundwater in the direction of its hydraulic

30

gradient potentially causing contamination of wells [40]. To represent spills in the vadose zone,
mixtures of 25% volume of oil and soil with 10% and 25% volumetric moisture content were
simulated. Measured results were compared with the predicted values calculated using the
mixing rules for moist and high-loss materials presented in [41]. Results are presented in Figure
2-17.

Figure 2-17. Results for a mixture of 25% volume of contaminant in soil with 10% and 25%
volumetric moisture content. Measured results are compared with the predicted
values calculated using the mixing rules for moist and high-loss materials presented
in [41].

A significant difference, i.e. greater than the nominal uncertainty in the measured data,
with regards to the permittivity values, was observed between the wet soil and wet soil with oil
contamination. This is an indication that if wet soil is present, it will provide sufficient contrast
to distinguish contaminated from uncontaminated soil.
As GPR signals travel through the soil, they are attenuated at a rate determined by the
complex dielectric constant of the soil. The power attenuation in dB/m is given by

31

2 f
Power Attenuation  8.6855
c



2 
 "
1     1


2
'




 '

(2.24)

8

where c is the speed of light 2.99710 m/s.
Figure 2-18 shows how changes in frequency relate to radar wave attenuation in different
soils. As frequency increases, the attenuation of the GPR signals in both wet and contaminated
soils increases rapidly. High frequency radar is often used to enhance resolution since resolution
increases with frequency, but as shown in this figure, signal attenuation increases quite
dramatically at higher frequencies.

Figure 2-18. Power attenuation in dry, wet and contaminated soil in a frequency range of 100
MHz to 1 GHz.

The next chapter presents the analytical formulation of the algorithm used in the
simulation of contaminants in soil using ground-truth obtained in these experiments.

32

CHAPTER 3
ANALYTICAL FORMULATION OF ELECTROMAGNETIC SCATTERING FROM AN
UNDERGROUND HYDROCARBON SPILL

3.1

Introduction
Location and imaging of objects in the subsurface is an important engineering problem

facilitating the design of systems utilized for the detection of mines, underground structures,
pipelines, faults, waste pollution, etc. Among the techniques proposed for this purpose, those
based on electromagnetic wave scattering are widely used [42]. Scattering from buried dielectric
objects has been studied extensively by many researchers from the theoretical and numerical
point of view during the past 4 decades.
Large volumes of hazardous liquids and fuels are stored worldwide in surface and
underground tanks. Frequently, these tanks are found to leak resulting in loss of inventory and,
more importantly, contamination to soil and groundwater. There are two methods for detecting
leakage from tanks: monitoring the liquid level of fuel in the tank and monitoring the soil under
the tanks for leaks. Chemical detection of leaks is very expensive, because it requires that
hundreds of chemical sensors be placed around and underneath the tank to ensure reliable
detection of the chemical spill [43]. Instead, borehole logging can be used to detect leaks in
underground storage tanks and underground fuel transfer lines. Borehole logging is a
nondestructive measuring technique and logistically simple and economical compared to other
exploration methods [44]. Geophysical borehole logging involves gradually lowering a probe
down a borehole, while the probe measures a physical property of the surrounding rock or soil.

33

Probes can be designed to measure any one of a variety of physical properties [45]. Since the
measured physical property is related to the composition of the surrounding rocks and soils,
borehole logs can be used to map the subsurface.
In this study, a simulated GPR borehole logging system was used as the method for
detection. Figure 3-1 is a representation of a possible underground spill and borehole logging
system. It is essential that the soil dielectric characteristics, presented in Chapter 2, are
considered when analyzing the GPR response from dielectric targets. For the quantitative
analysis of the subsurface, modeling of the electromagnetic interaction of the waves with the
targets, is essential. Computational models provide basis for identifying target signal features and
serve as the forward model for inversion analysis and reconstruction [42].

Borehole

Non-contaminated
soil
Oil

Probe
Residual contamination
Oil volatiles

Water table

Figure 3-1. Schematic diagram of a possible underground spill and borehole logging system. For
simplicity, the soil is modeled as a homogeneous dielectric space and the
contaminated area is modeled as two eccentric, homogeneous circular cylinders.

34

In this chapter, solutions for scattering from a dielectric cylinder of infinite length (e.g. a
two-dimensional, normal incidence scattering scenario) are reviewed. Values of permittivity
obtained using the coaxial waveguide presented in Chapter 2 are used for a computer simulation
of a contaminated underground system. For simplicity, the soil was modeled as a homogeneous
dielectric space and the contaminated area was modeled as two eccentric, homogeneous circular
cylinders as shown in Figure 3-2. Cylinders with different permittivities were considered to be
sufficiently representative of a contaminated region for academic purposes.

Contaminated
soil (Oil)

Noncontaminated
soil

Residual
contamination
Oil volatiles

Figure 3-2. Cross-sectional view of the spill illustrated in Figure 3-1. The spill is modeled as an
eccentric dielectric cylinder embedded in another dielectric cylinder. The noncontaminated soil is modeled as a homogeneous dielectric space

First, scattering equations from a homogeneous dielectric cylinder of infinite length
immersed in a lossless homogeneous medium are presented. Convergence analysis and contrast
analysis was performed to assure convergence of the Bessel and Hankel functions and scattering
field intensity, respectively. Then, a system of eccentric cylinders is analyzed. The cylinders are
assumed to be infinite in length and their axes are parallel to the z-direction. The system is
35

illuminated by a cylindrical wave from a line source of infinite extent located at  0 , 0  in a
direction parallel to the cylinder axis. Hence the incidence angle is normal rather than skew
resulting in a non-depolarizing scattering scenario. Electromagnetic scattering from the cylinders
is calculated using a boundary-value mode-matching approach presented by Kishk, Parrikar, and
Elsherbeni in [46], called the KPE algorithm. The approach to the problem involves solving the
scalar Helmholtz equation in cylindrical coordinates subject to the boundary conditions imposed
by Maxwell’s equations.
3.2

Scattering by dielectric objects

3.2.1 Normal incidence plane wave scattering by dielectric circular cylinder: fundamental
theory.
Cylinders represent an important class of objects for GPR since they represent important
environmental and engineering targets and also because their scattering properties are strongly
polarization dependent [47]. Polarization in this context refers to the orientation of the electric
field vectors. The electric field of a wave traveling in a direction normal to the z-axis can be
described by two orthogonal components as given by [48] [49]:

Ex  z, t   Ex0e z cos t   z  x 



E y  z, t   E y0e z cos t   z   y



(3.1)

(3.2)

where α represents the attenuation constant, β the phase constant, ω the angular frequency,  the
phase, and Ex0 and Ey0 are the maximum amplitudes of the Ex and Ey components, respectively.

36

Polarization vectors are chosen such that one vector is oriented along the long axis of the
cylinder (E-parallel or transverse-magnetic (TM)) and the other vector oriented orthogonal to the
long axis of the cylinder (E-perpendicular or transverse electric (TE)) as shown in Figure 3-3.

y


a
x

E-parallel (TM)

y


a
x

E-perpendicular (TE)
Figure 3-3. Definitions of E-parallel (TM) and E-perpendicular (TE) polarizations relative to a
cylinder, as defined in [49]. TM polarization occurs when the electric field is parallel
to the long axis of the cylinder. TE polarization occurs when the electric field is
perpendicular to the long axis of the cylinder.

37

The scattered field is a function of the electrical properties of the cylinder and
surrounding material, the distance from the cylinder, and the scattering angle , as defined in
Figure 3-3. Since circular cylinders have a cylindrical shape, their scattering properties are
conveniently described using Hankel and Bessel functions.

TMz polarization
A TMz uniform plane wave traveling in a medium with permittivity  and permeability

 0 is normally incident on a lossless dielectric circular cylinder of radius a, permittivity  d and
permeability  0 as illustrated in Figure 3-3. The time dependence e jt is implied and


suppressed throughout the derivation. The incident field E i can be expressed as
jk x cos i  y sin i 
jk  cos  cos i sin  sin i 
ˆ i ˆ
ˆ
E i  a z E z  a z E0 e 0 
 a z E0 e 0 
jk  cos i 
ˆ
 a z E0 e 0

(3.3)

where E0 is the amplitude of the incident wave, k0     0 is the wavenumber in the
incident medium, i is the incident angle,  is the observation angle, and  is the observation
distance. Using wave transformation from Cartesian to cylindrical coordinates [48],

e jx  e j  cos  





j n J n (  )e jn

(3.4)

n 
where J n   is the Bessel function of the first kind of order n. The incident electric field can be
expressed as

38



i
Ez  E0



j n J n (k0  )e

jn i 

(3.5)

n 
The total field in the incident medium is the sum of the incident and the scattered fields, that is,

s
ˆ( i
E  Ei  E s  z Ez  Ez )

(3.6)

The scattered field Es can be represented by



E
 0
 n 
s 
Ez  



 E0
 n 





jn  i 
(2)
j n an H n (k0  )e 

,  a
(3.7)

j ncn J n (kd  )e

jn i 

,  a

(2)
where a n and cn are unknown amplitude coefficients, H n   is the Hankel function of the

second kind of integer order n, and kd    d 0 is the wavenumber inside the dielectric.
These unknown amplitude coefficients can be found by applying the following boundary
conditions at   a :

ref
i
tran
Ez  Ez  Ez

(3.8)

ref
i
tran
H  H  H

(3.9)

For the electric field, using equations (3-3) and (3-5) into equation (3-6) we obtain

E tran  E0





n 

jn  i 
(2)
j n  J n  k0a   an H n  k0a   e 





Simplifying and solving for cn,

39

(3.10)



E0



jn  i 
j ncn J n  kd a  e 
 E0

n 





n 

jn  i 
(2)
j n  J n  k0a   an H n  k0a   e 





(2)
cn J n  kd a   J n  k0a   an H n  k0a 

cn 

(2)
J n  k0a   an H n  k0a 

(3.11)

J n  kd a 

ref
i
tran
Now, for the boundary condition for the magnetic field, H  H  H , we need to find

the H-field in each region. Using Faraday’s equation

H 

1
j

1 
1 E
E 
ˆ
ˆ
 a
 a

j 
 
 

 E  

(3.12)

H  

1  E 
ˆ
 a

j   

(3.13)

H  

1 
E 
1 Ez
ˆ
 a    j 
j 


(3.14)

i k E
H  0 0
j

ref

H



jn  i 

j n J n  k0   e 



(3.15)

n 

k E
 0 0
j





jn  i 
(2)
j n an H n  k0   e 

(3.16)

n 

tran  kd E0
H
j






j nCn J n  kd   e

n 

40

jn i 

(3.17)

(2)

where J n   is the derivative of the Bessel function of the first kind of order n and H n    is

the derivative of the Hankel function of the second kind of integer order n. At the boundary of
the cylinder (   a ) the total tangential magnetic field can be written as

k0 E0
j





n 

jn  i  kd E0
(2)

j n  J n  k0a   an H n  k0a   e 



j







jn  i 

j ncn J n  kd a  e 

n 

Simplifying,

(2)


  J n  k0a   an H n  k0a    d cn J n  kd a 





 J  k a  a H (2) ' k a    d c J  k a
 n  0  n n  0 
 n n d 



(3.18)

Using equation (3-9) and defining  c   d  as the contrast ratio of the permittivity of the
dielectric and the permittivity of the medium,


(2) '

J n  k0a   an H n  k0a   d

 J k a  a H (2) k a 
 n  0  n n  0   J n  kd a 

 

J n  kd a 



 J k a J  k a  a H (2) k a J  k a 
  n  d  n n  0  n  d 
(2) '

J n  k0a   an H n  k0a    c  n 0


J n  kd a 



(2)


 c an H n  k0a  J n  kd a 
 J  k a  J n  kd a 
(2) '

an H n  k0a  
 c n 0
 J n  k0a 
J n  kd a 
J n  kd a 
 H (2) ' k a J k a   H (2) k a J   a 
 0  n  d  c n  0  n  d  
an  n


J n  kd a 




 c J n  k0a  J n  kd a   J n  k0a  J n  kd a 
J n  kd a 

41

an 



 c J n  k0a  J n  kd a   J n  k0a  J n  kd a 

(3.19)

(2) '
(2)

H n  k0a  J n  kd a    c H n  k0a  J n  kd a 

Substituting equation (3-17) into equation (3-9) we obtain
(2)


 c J n  k0a  J n  kd a   J n  k0a  J n  kd a  
H n  k0 a  


cn 

J n  kd a 
J n  kd a   H (2) '  k a  J  k a    H (2)  k a  J   k a  
0
n d
c n
0
n d 
 n
J n  k0 a 




 c J n  k0a  J n  kd a   J n  k0a  J n  kd a   (2)

H
J n  k0 a  
 k0a 
 H (2) '  k a  J  k a    H (2)  k a  J   k a   n
0
n d
c n
0
n d 
 n

J n  kd a 
  H (2) k a J k a J  k a  H (2) k a J  k a J k a 
 0  n  0  n  d  n  0  n  0  n  d 
J n  k0 a    c n
(2) '
(2)



H n  k0a  J n  kd a    c H n  k0a  J n  kd a 



J n  kd a 

(2)
(2)

J n  k0a  H n  k0a   H n  k0a  J n  k0a 
cn 
(2) '
(2)

H n  k0a  J n  kd a    c H n  k0a  J n  kd a 

(3.20)

For the TM polarization normally incident on a dielectric cylinder, the scattering field is
given by:

s
E z  E0





n 

jn



 c J n  k0a  J n  kd a   J n  k0a  J n  kd a 
(2) '
(2)

H n  k0a  J n  kd a    c H n  k0a  J n  kd a 

(3.21)

jn  i 
(2)
 H n ( k0  ) e 

To obtain the scattered field in the far-zone, the asymptotic expansion for large argument
of the Hankel function is used [49], [50].

k0  
(2)
H n  k0  


2 j n  jk0 
j e
 k0 

42

(3.22)

From equation (3-5), the reflected field is


ref
E z  E0



jn  i 
(2)
j n an H n (k0  )e 

n 


  E0



jn  i 
(2)
j n an H n (k0  )e 

(3.23)

n 

Therefore, the far-zone approximation of the scattered field is given by

s
Ez   E0

2 j  jk0 
e
 k0 





ane

jn i 

(3.24)

n 

TEz polarization
The previous analysis is repeated for a TEz plane wave normally incident on a dielectric
cylinder of radius a, permittivity  d and permeability  0 as shown in Figure 3-3. The TE
polarization case is expressed in terms of magnetic fields for convenience and one may convert
between electric and magnetic fields using Ampere's and Faraday's laws. The incident magnetic
field can be written as

jk x cos i  y sin i 
ˆ
ˆ
H i  a z H 0e 0 
 az H 0





j n J n  k0   e

jn i 

(3.25)

n 
The corresponding incident electric field can be obtained by using Ampere’s law, which reduces
to
Ei 

1
j

 Hi 

i
i
H z 
1 
1 H z
 a

 a
j 
 
 



The components of the electric field are

43

(3.26)

i
i  1 1 H z  H 0 1
E
j  
j 

i
i   1 H z   k0 H 0
E
j 
j



 nj n1Jn  k0  e jn i 

(3.27)

n 



jn  i 
'
j n J n  k0   e 



(3.28)

n 

The scattered magnetic field has a similar form to that of the scattered electric field of equation
(3-5) for the TMz polarization. It can be written as

ˆ s ˆ
H s  az H z  az H 0



(2)
 dn Hn  k0 e jn i 

(3.29)

n 

where dn represents the unknown amplitude coefficient. The scattered electric fields are

s H 1
1 1 H z
s
E 
 0
j  
j 



(2)
 njHn  k0  dne jn i 

(3.30)

n 

s
k H
1 H z
s
E 
 0 0
j 
j





jn  i 
(2)
dn H n  k0   e 

(3.31)

n 

Applying boundary conditions at   a
ref
i
tran
E  E  E

(3.32)

ref
i
tran
E  E  E

(3.33)

ref
i
tran
Hz  Hz  Hz

(3.34)

Solving for the  component of the electric field

44

tran   k0 H 0
E
j





n 

 n
jn i 
jn  i  
(2)

 d n H n  k0a  e 
 j J n  k0a  e




Simplifying,

 kd H 0
j





n 

k H
jn  i 

j nCn J n  kd a  e 
 0 0
j

(2)

  j n J n  k0a   dn H n   k0a  e jn i  

k
(2)


Cn J n  kd a   0  J n  k0a   d n H n  k0a  

d
 


kd

(2)
 

k0 d  J n  k0a   d n H n  k0a  
Cn 


kd  
J n  kd a 


 J  k a  d H (2) k a 
  n n  0 
Cn   c  n 0



J n  kd a 



(3.35)

ref
i
tran and solving
Applying the third boundary condition for the magnetic field H z  H z  H z
for dn,

tran  H
Hz
0





jnJ

n  k0a  e

jn i 



 H0

n 



H0

 Cn Jn  k0a  e

n 



jn i 
(2)
j n d n H n  k0 a  e

n 

jn i 



 H0



n 

 j n J k a  j nd H (2) k a e jn i 
n 0 
n n  0 




(2)
Cn J n  kd a   J n  k0a   dn H n  k0a 

 J  k a  d H (2) k a 
  n n  0   J k a  J k a  d H (2) k a
c  n 0
  n 0  n n  0 


 n d
J n  kd a 



45



(2)

 c J n  k0a  J n  kd a   dn  c H n  k0a  J n  kd a  
(2)


J n  kd a  J n  k0a   d n J n  k0a  H n  k0a 
(2)
(2)

dn   c H n  k0a  J n  kd a   J n  kd a  H n  k0a   






J n  k d a  J n  k0a    c J n  k0a  J n  k d a 

dn 



J n  kd a  J n  k0a    c J n  k0a  J n  kd a 
(2)
(2)

 c H n  k0a  J n  kd a   J n  kd a  H n  k0a 

(3.36)

For TEz polarization, the scattered electric field can be summarize as
ˆ
ˆ
Es   E   E

(3.37)

where

E 1
s
E  0
j 

k E
s
E   0 0
j

and  





jn  i 
(2)
njdn H n  k0  e 

(3.38)

n 




jn  i 
(2)
dn H n  k0   e 

(3.39)

n 

0
is the intrinsic impedance of the incident medium. Field equations for both, TM


and TE polarizations are expressed in terms of Bessel and Hankel functions. Before performing
any calculations, we need to analyze the convergence of the series. The magnitude of the bistatic
field for TMz and TEz polarization was calculated for various numbers of terms in the series
expansion. The series was truncated when the percent difference between the results of
consecutive modes was less than 0.01%. Results are shown in Figure 3-4.

46

TMz

TEz

Figure 3-4. Magnitude of the bistatic scattering field for a dielectric cylinder with radius /2,
surrounded by wet soil. Different numbers of modes were used to calculate the
scattering field and analyze the convergence of the series expansion. The series was
truncated when the percent difference between the results of consecutive modes was
less than 0.01%.

47

Results for both polarizations showed that the required number of terms N in the series is
approximately N  2k0r , where k0 is the wavenumber in the surrounding medium and r is the
radius of the dielectric cylinder.
Scattering width
The radar cross-section (RCS) represents a convenient way to describe the strength of
scattered fields observed in the far-field. The RCS is defined as the area intercepting the amount
of power that when scattered isotropically, produces at the receiver a density that is equal to the
density scattered by the actual target [51]. For 2D objects such as an infinite cylinder, the RCS is
represented by the scattering width (SW) or alternatively the RCS per unit length [49].

SW  lim 2 r
r 

Es
Ei

2
2

(3.40)

Most GPR surveys used to image the subsurface are conducted in common-offset mode
using closely spaced antennas [52]. This results in a bistatic angle of approximately 0°
(backscattered), depending on antenna separation and target depth. It is thus important to further
investigate backscattered scattering widths as a function of cylinder radius. TMz and TEz
scattering widths were calculated for the case of a dielectric cylinder with relative permittivity of
2.22 surrounded by a medium with permittivity of 17.0. This represents targets like hydrocarbons
surrounded by wet loamy soil. Figure 3-5 illustrates backscattering width as a function of radius
for TM and TE polarization for the intended target. The results show that, for small cylinders
immersed in the intended higher permittivity medium, TM backscattering widths are greater than
TE backscattering widths for most radii considered herein.
48

Figure 3-5. Backscattering widths as a function of radius, normalized by the wavelength (λ) of
the incident field. TM backscattering widths are greater than TE backscattering
widths for most cylinders.

Several features are observed in the scattering widths as a function of scattering angle for
dielectric cylinders. As the radius of the cylinder becomes small compared to wavelength the
scattering width amplitude oscillates less as one would expect since the contributions to the
scattered field at any given observation angle has less phase difference (e.g. in the limit, all of the
scattered field contributions emanate from the same point in space). The TM polarization
scattering widths become nearly constant as a function of scattering angle. This is illustrated in
Figure 3-6.

49

Figure 3-6. Scattering widths for high contrast dielectric cylinders, normalized by the wavelength
(λ) of the incident field. As the radius of the cylinder becomes small compared to
wavelength, TM scattering widths become nearly constant as a function of scattering
angle.

Figure 3-7 shows the scattering width as a function of frequency for cylinders of different
radii. The radii are varied from 0.05 to at the highest frequency. The scattering width
increases as the frequency increases for both, forward- and backscattering. Also, as the radius of
the cylinder increases, more oscillations appear in the scattering width amplitudes.

50

Figure 3-7. Scattering width as a function of frequency for cylinders of different radii. The radii
are varied from 0.05 to  at the highest frequency. Forward (top) and back (bottom)
scattering widths are illustrated.

51

The scattering width can also be plotted as a function of the contrast ratio  c as shown in
Figure 3-8. The response is similar for high- and low-impedance dielectric cylinders. As the size
of the cylinder increases, oscillations in the scattering width response increases. Also, if the
contrast ratio approaches unity, the scattered fields will vanish, as expected.

Figure 3-8. Scattering width as a function of the contrast ratio  c . Similar response is obtained
for high- and low-impedance dielectric cylinders.

3.2.2 Layered dielectric cylinders: theory and analytical solutions
In this section, the equations needed for the calculation of the scattering field from an
inhomogeneous spill are derived. The scattering field is calculated using a boundary-value mode-

52

matching approach presented by Kishk, Parrikar, and Elsherbeni in [46], called the KPE
algorithm. The spill is assumed to be a system composed of a circular dielectric cylinder
embedded into another and surrounded by a dielectric medium as shown in Figure 3-2. The
system is illuminated by a cylindrical wave from a line source of infinite extent located at

 0 ,0 

in a direction parallel to the cylinder axis. The TMz wave traveling in a medium with

permittivity  S and permeability 0 (Region 0) is incident on a layered dielectric circular
cylinder of outer radius b with permittivity 1 (Region I), and inner radius a with permittivity

 M (Region II) as illustrated in Figure 3-9.

Line source

Y1

YM

a

b

M

Region II

Region 0


XM
X1

Region I

S

1

Figure 3-9. Geometry of the layered dielectric cylinders. The TMz wave traveling in a medium
with permittivity  s and permeability 0 (Region 0) is incident on a layered
dielectric circular cylinder of outer radius b with permittivity 1 (Region I), and
inner radius a with permittivity  M (Region II)

53

The fields in the dielectric regions are expressed in terms of a set of cylindrical harmonic
functions with unknown coefficients. The time dependence e jt is implied and suppressed
i
throughout the derivation. The incident electric field E z is given by [49]:

 k I (2)
i
Ez   ,     0 0 H 0  k0   0 
4

(3.41)

where I is the current of the line source, k0    s 0 is the wave number in the incident
medium,  s is the permittivity of the soil, and 0  0  s is the intrinsic impedance of the

i
medium. E z can be expressed in terms of cylindrical harmonic functions as


jn  0 
(2)
 0k0 I
H n  k0 0  J n  k0   e 
,

4
n 

i
Ez   ,    


jn  0 
(2)
 0k0 I
H n  k0   J n  k0 0  e 
,

4
n 





  0
(3.42)

  0

(2)
where J n and H n are the Bessel function of the first kind and Hankel function of the second
kind, respectively. The line source excitation is given by

 k I
E0   0 0 .
4

(3.43)

The z-component of the total electric field in region 0 (soil) is expressed as

0
Ez  1, 1   E0





n 

 H (2) k  J k   B0 H (2) k   e jn1 0 
 n  0 0  n  0 1  n n  0 1 



54

(3.44)

0
where Bn are the unknown scattering coefficients. The electric field inside each region p (p = I,
II) is expressed with respect to the (XM, YM) coordinate system as

p
Ez   M , M   E0







n 







p (2)
ApJ k 
 Bn H n k p  M  e jnM
 n n p M




(3.45)

p
p
where k p    p  p and An , Bn are the unknown scattering coefficients in each region. To
translate the field to the (X1, Y1) coordinate system, the addition theorem for cylindrical
functions was used [53]. In summary, the addition theorem establishes that

 displaced wave 


 of order  

 undisplaced wave 
.
of order +m 

 Jm  kr0   
m

Therefore, the translated field equation from the (XM, YM) coordinate system to the (X1, Y1)
coordinate system is given by





p
E z  p ,  p  E0





  Ji n  k prpM 

n  i 

(3.46)

j i  i  n  pM 
p
p (2)

  An J i k p  p  Bn H i
k p p  e  p













where rpM is the distance between the origin of the (X1, Y1) coordinate system and the origin of
the (XM, YM) coordinate system and  pM is the angle between the origin of (X1, Y1) and the
origin of (XM, YM) coordinate systems.
For the magnetic field, the z component of the total magnetic field in region 0 is given by

55

E
0
H  1, 1   0
j0





n 

 H (2) k  J  k   B0 H (2) k   e jn1 0  (3.47)
 n  0 0  n  0 1  n n  0 1 



and inside each region p (p = I, II) is expressed with respect to the (XM, YM) coordinate system
as

E
H   M , M   0
j p



p



n 









p (2)
ApJ k 
 Bn H n
k p M  e jnM (3.48)
 n n p M




where 0 andp are the intrinsic impedances of region 0 and region p, respectively. Translating
the equation to the (X1, Y1) coordinate system





E
p
H  p ,  p  0
j p





  Ji n  k p rpM 

n  i 

(3.49)

j i  i  n  pM 
p
p (2)
 .
  An J i k p  p  Bn H i
k p p  e  p













Boundary conditions
The present case involves two embedded cylinders. There are two boundary surfaces
available on which the boundary conditions for the electric and magnetic fields are to be applied.
Boundary conditions at the external surface (p = I)
At the outermost surface  I  rI  b the boundary conditions are

0
I
Ez  rI , I   Ez  rI , I 

56

(3.50)

and
0
I
H  rI , I   H  rI , I  .

(3.51)

For the electric field,



E0



n 

 H (2) k r J k b  B 0 H (2) k b  e jnI 0 
 n  0 0  n  0  n n  0 




 E0



  Ji n  kI rIM 

(3.52)

n  i 

j i  i  n IM 
I
I (2)

  An J i  k I b   Bn H i  k I b   e  I




Two complex exponential functions are said to be orthogonal over [-, ] if the inner product is
equal to zero. If F ( x)  e jnx and G( x)  e jkx , the inner product is

F ( x), G ( x) 







 F ( x)  G( x) dx


 e

j (n  k ) x dx



 e j (n  k ) x 


0
 j (n  k ) 

 
Using orthogonality of complex exponential functions and equating to zero yields

57

(3.53)

 H (2) k r J k b  B0 H (2) k b  e jn0
 q  0 0  q  0  q q  0 








n 

 j q  n IM
I
I (2)
J q  n  k I rIM    An J q  k I b   Bn H q  k I b   e 
.





(3.54)

0
Solving equation (3.53) for Bq results in

0
Bq 


1
(2)
 jq0
 H q  k0r0  J q  k0b   e
J q  n  k I rIM
(2)
H q  k0b 
n 






(3.55)

 j q  n IM 
I
I (2)
  An J q  k I b   Bn H q  k I b   e 
.






0
Following the same procedure for the magnetic field and solving for Bq yields

0
Bq 


0  jq0
(2)

  H q  k0r0  J q  k0b   e
J q  n  k I rIM
(2)
I
H q  k0b 
n 
1







 j q  n IM 
I 
I (2)
  An J q  k I b   Bn H q   k I b   e 
.






0
Equating (3.54) and (3.55) to eliminate Bq ,

58

(3.56)


1
(2)
 jq0
 H q  k0r0  J q  k0b   e
J q  n  k I rIM
(2)
H q  k0b 
n 







 j q  n IM 
I
I (2)
  An J q  k I b   Bn H q  k I b   e 







(3.57)

0  jq0
(2)


  H q  k0r0  J q  k0b  
e
J q  n  k I rIM
(2)
I
H q   k0b 
n 
1







 j q  n IM 
I 
I (2)
  An J q  k I b   Bn H q   k I b   e 







or, rearranged (3.57) becomes,





n 


 jq0  j  q  n IM

 J q  n  k I rIM  e



(2)

 (2)
 

  J q  k I b   0 J q  k I b   A I   H q  k I b   0 H q  k I b   B I  

n  (2)
n 
 (2)
(2)
(2)
 H q  k0b   I H q   k0b  

 H q  k0b   I H q   k0b    
 



 


(3.58)

 J k b

J q  k0b  
q 0 
(2)
 .
 H q  k0r0  

 (2)

(2)

 H q  k0b  H q  k0b  

Boundary conditions at the core surface (p = II)
At the surface of the core (e.g. inner) cylinder, continuity of the tangential components of
the electric and magnetic fields is enforced. For the electric field at  II  rII  a
II
I
Ez  rII , II   Ez  rII , II 

Therefore,

59

(3.59)





n 

 II
 jnII
II (2)
 An J n  k II a   Bn H n  k II a   e



(3.60)






n 

 A I J k a  B I H (2) k a  e jnII .
 n n  I  n n  I 



II
The fields at the core cylinder are to be finite. Therefore, the condition Bn  0 is imposed.
Simplifying (3.60) yields,




n 

 AII J  k a   e jnII 
 n n II 







n 

 AI J k a  B I H (2) k a  e jnII
 n n  I  n n  I 



(3.61)

Applying orthogonality of exponential functions (3.53) and equating to zero,

(2)
II
I
I
J q  kII a  Aq  J q  k I a  Aq  H q  k I a  Bq  0 .

(3.62)

Following the same procedure for the magnetic field,
II
I
H  rII , II   H  rII , II 

(3.63)

II 1 J  k a AI  1 H (2) k a B I  0 .

J q  k II a  Aq 
  q
 I  q
 II
I q I
I q

(3.64)

and

1

II
Using equations (3.62) and (3.64) to eliminate Aq yields
(2)
II  J n  k I a  AI  H n  k I a  B I
An
n
n
J n  k II a 
J n  k II a 

60

(3.65)

and

(2)

II   II J n  k I a  AI   II H n   k I a  B I .
An


 I J n  k II a  n  I J n  k II a  n

(3.66)

Equating (3.65) and (3.66),
(2)
(2)
 J n  k I a   II J n  k I a   I  H n  k I a   II H n   k I a   I


 B  0.



 An 


 J n  k II a   I J n  k II a   n
 J n  k II a   I J n  k II a  





(3.67)

Equations (3.58) and (3.67) constitute a system of linear equations of the unknown

I
I
I
I
coefficients An and Bn . Once these equations are solved, values for An and Bn are used in
0
equation (3.55) or (3.56) to compute the scattering coefficients Bn . The values of the scattering
coefficients are then used to compute the fields in each region due to a line source excitation.
As established in section 3.2.1, an important parameter in scattering is the scattering
width. The scattering width is obtained by calculating the ratio of the scattered field in the far
zone to the incident field.

E s   ,  , 0 
    lim 2 z
i
 
E z   ,  , 0 

2
(3.68)

From equation (3-44), the scattered field is

s
Ez  1, 1   E0





jn  
0 (2)
Bn H n  k0 1  e  1 0 

n 

61

(3.69)

where the value of E0 depends on the excitation (line source or plane wave). Using the large
argument of the Hankel function and normalizing the resulting expression by a factor of

E0





2 j  k  e jk  , yields

s
E z  , o  





0
Bn j ne

jn 0 

.

(3.70)

.

(3.71)

n 

For the composite cylinder, the scattering cross-section is given by

   

2







0
Bn j ne

jn 0 

2

n 

Near-field distribution and scattering width for the composite cylinder are computed in
the next section for both concentric and eccentric configurations.
3.3

Simulation of electromagnetic scattering from an underground hydrocarbon spill
The expressions presented in section 3.2 were used to simulate the electromagnetic

scattering from an underground hydrocarbon spill. The spill is assumed to be a system composed
of a circular dielectric cylinder embedded into another and surrounded by a dielectric medium as
shown in Figure 3-2. The center frequency used is 500 MHz. This frequency is low enough for
good penetration into the soil, and high enough for good resolution of small objects. An object is
considered small if its radius is less or equal to the wavelength of the incident wave. Even though
a center frequency of 500 MHz is used, the system is assumed to radiate substantial energy from
100 MHz to 1 GHz. Different sizes and configurations of contaminated spills were simulated.

62

Results for concentric and eccentric configurations illuminated by a plane and a line source are
presented.
3.3.1 Homogeneous spills illuminated by a plane wave
The expressions for the electric fields in each region were verified using known
configurations. First, the case of two concentric cylinders made of the same material in free
space was considered. This was achieved by letting  2  1 and 2  1 . The incident field is a
plane wave with a frequency of 500 MHz. The incident wave has amplitude of one and is
traveling along the +y direction. As the wave propagates along +y direction, it is disturbed by a
dielectric circular cylinder with  r  2.2  j 0 at 500 MHz and radius equal to 1 wavelength of
the incident wave. A target with permittivity of 2.20 represents oil contamination. The
scattering width for TM excitation is presented in Figure 3-10(a). These results are as expected
for a dielectric cylinder in free space. The electric field magnitude scattered by the homogeneous
dielectric cylinder in air is presented in Figure 3-10(b). The incident plane wave is diffracted by
the dielectric cylinder, creating a maximum peak in the forward direction and standing waves
appear in the backward direction. The peak observed at the back side of the cylinder decreases
rapidly, leaving a shadow area in the forward direction.
Now, soil dielectric characteristics are considered to analyze the GPR response from the
dielectric target. Figures 3-11 and 3-12 show the modeled scattered field from the same dielectric
target in two different backgrounds: dry and wet soil. The scattering width for the system of
homogeneous wet soil and a spill with a radius equal to the incident wavelength in wet soil at
500 MHz is shown in Figure 3-13.

63

(a)

(b)
Figure 3-10. Scattering by a homogeneous cylinder in air. (a) Scattering width  2 D /  as a
function of scattering angle . 0  180 ,  r  2.2  j 0 , and r = . (b) The
simulated electric field magnitude. The black line denotes the cylinder location.

64

(a)

(b)
Figure 3-11. Scattering by a homogeneous oil spill in dry soil. (a) Scattering width  2 D /  as a
function of scattering angle . 0  180 ,  r  2.2  j 0 , and r = . (b) The
simulated electric field magnitude. The black line denotes the cylinder location.

65

(a)

(b)
Figure 3-12. Scattering by a homogeneous oil spill in wet soil. (a) Scattering width  2 D /  as a
function of scattering angle . 0  180 ,  r  2.2  j 0 , and r = . (b) The
simulated electric field magnitude. The black line denotes the cylinder location.

66

(a)

(b)
Figure 3-13. Scattering by a homogeneous wet-contaminated spill in wet soil. (a) Scattering
width  2 D /  as a function of scattering angle . 0  180 ,  r  12  j 0 , and r
= . (b) The simulated electric field magnitude. The black line denotes the
cylinder location.

67

Note that the field distribution in Figures 3-11(b) and 3-12(b) from the dry and wet soil
background cases, both viewed in cross section, are completely different. The dielectric constant
differences between dry and wet soil generate remarkable differences in these two cases. When
the soil is wet, the permittivity increases and more contrast between the soil and the spill is
obtained. In the vadose region, the soil and contaminated soil is mostly uniformly wet. Assuming
that the soil has 25% water content, the measured value of relative permittivity is approximately
17. The same amount of water content in the contaminated soil has a relative permittivity value
of 12. The scattering width for the system of homogeneous wet soil and a spill with a radius
equal to the incident wavelength in wet soil at 500 MHz is shown in Figure 3-13. Notice that the
magnitude of the signal is attenuated at this incident frequency. Soil moisture increases the radar
attenuation rates, limiting radar performance.
3.3.2 Inhomogeneous spill illuminated by a plane wave
A more realistic approach to the problem is to assume that the spill is not homogeneous,
but composed of separated layers with different permittivity. To simulate this, the spill is
assumed to be a layered dielectric circular cylinder of outer radius b with permittivity 1 , and
inner radius a with permittivity  2 as illustrated in Figure 3-9.
The first configuration analyzed is a concentric cylinder with small radius relative to the
wavelength in wet soil of the incident wave. The spill has a relative permittivity of 1  12 in the
outer layer, corresponding to contaminated soil with 25% water content and 25% oil content, and

 2  2.2 in the core region. The permittivity of 2.2 corresponds to fully contaminated soil or soil
with more than 50% oil. The pores of the soil in the core region are assumed to be saturated with
oil and the water had been displaced to the outer region. The total radius of the spill is 2 and the

68

core region radius is equal to  of the field wave in wet soil. The simulated electric field
magnitude for an incident plane wave with a frequency of 500 MHz is shown in Figures 3-14.
The scattering pattern in the wet soil region is similar to the scattering pattern for a homogeneous
cylinder, but a reduction in the scattered field magnitude is observed in the direction of
propagation at x = 0. There is a significant contrast between the outer cylinder and the core
cylinder. The change in wavelength is observed. This difference in permittivity produces a welldefined shadow area in the forward direction that starts in the outer cylinder. The wave fronts do
not reform after passing through the contaminated region, making this a useful characteristic for
determining the concentration of contaminant in the soil.

Figure 3-14. Simulated electric field magnitude for a layered small spill in wet soil at 500 MHz.
The inner circle represents pure contamination of oil (r = 2.2) and the outer
cylinder represents wet-contaminated soil (r = 12). The cylinders are surrounded
by wet soil.
69

Not all oil spills spread radially from the source. Depending on the soil and other
parameters, contamination could spread at different rate in different directions. A spill with a
shifted core region is illustrated in Figure 3-15. As before, the core is assumed to be fully
contaminated soil and the water particles are pushed to the surroundings. The size of the core is
1/3 the size of the total contaminated area. There is a shadow area behind the core cylinder.
When the wave transition into the wet soil area, the wave fronts start to regenerate. The direction
of propagation of these wave fronts is deflected by the position of the core.

Figure 3-15. Simulated electric field magnitude for an eccentric spill in wet soil at 500 MHz. The
inner circle represents pure contamination of oil (r = 2.2) and the outer cylinder
represents wet-contaminated soil (r = 12). The cylinders are surrounded by wet
soil.

70

3.3.3 Line source illumination
The electric field magnitude for a small spill in the vadose region illuminated by a line
source with a frequency of 500 MHz is shown in Figure 3-16 for the homogeneous case. The
wet-contaminated soil has

r = 12 and radius equal to one wavelength of the incident wave is

assumed. The wet soil has a relative permittivity of 17. The wave travels from a high permittivity
medium to a lower permittivity target.

Figure 3-16. Electric field magnitude for a spill in a vadose region. The blue dashed line
represents the position of the spill. The line source is located at 3soil from the
spill.
A concentric spill with a relative permittivity of 1  12 in the outer layer, and  2  2.2
in the core region is now presented. As before, the pores of the soil in the core region are

71

assumed to be saturated with oil and the water had been displaced to the outer region. The total
radius of the spill is 2λ and the core region radius is λ of the field wave in wet soil. The electric
field magnitude for an incident plane wave with a frequency of 500 MHz is shown in Figure 317. Standing waves appear in the wet-contaminated region of the spill due to backward scattering
from the core region. At the core, a single period of the wave can be seen, with minima and
maxima amplitude.

Figure 3-17. Simulated electric field magnitude for a concentric spill in wet soil illuminated by a
line source with frequency of 500 MHz. The inner circle represents pure
contamination of oil (r = 2.2) and the outer cylinder represents wet-contaminated
soil (r = 12).
An eccentric spill with its core region located at the bottom-right side of the contaminated
area was simulated to inspect its scattering characteristics. The size of the core is 1/3 the size of

72

the total spill. The shadow area behind the spill is disturbed by the scattering from the core that is
closer to the spill’s boundary. The waves are deflected in the direction of the shifted core. This is
shown in Figure 3-18.

Figure 3-18. Simulated electric field magnitude for an eccentric spill in wet soil at 500 MHz. The
inner circle represents pure contamination of oil (r = 2.2) and the outer cylinder
represents wet-contaminated soil (r = 12).
In this chapter, solutions for scattering from a dielectric cylinder of infinite length were
reviewed. Values of permittivity obtained using the coaxial waveguide presented in Chapter 2
were used for a computer simulation of the electric field magnitude in a contaminated
underground system. Scattered field matrices obtained will be used as input parameter in the
inversion algorithm presented in Chapter 4.

73

CHAPTER 4
INVERSION ALGORITHM

4.1

Introduction
In this chapter, an inversion algorithm for underground spill localization is presented.

Inversion techniques in the frequency domain are more common for microwave medical and
ultrasound imaging modalities, but have also been applied to stepped frequency GPR [54],[55].
The most common form of single frequency inversion is tomographic [56]. Tomography
generates an image of the entire illuminated region, making it useful for characterizing variations
of material type. Cross-well radar (CWR) tomography, otherwise known as cross-borehole
ground-penetrating radar (cross-borehole GPR) is a minimally invasive method that uses high
frequency electromagnetic waves transmitted and received by antennas in the subsurface to
image objects of contrasting dielectric properties [57]. Cross-borehole GPR measures
transmission as well as reflection signals, making it better suited for frequency domain inversion
[58],[59]. A 2D diffraction tomography inversion algorithm and results calculated using
synthetic data are presented in this chapter.
4.2

Diffraction Tomography
Diffraction tomography (DT) is based on the determination of the dielectric properties of

an object. The response of the object from both the electric and the magnetic field must be
evaluated, as well as the interaction between the sample and the probe. DT was first proposed by
Wolf in [60] and is now widely used in various forms for such applications as medical imaging,
optical imaging, geophysical tomography and radar imaging [60]-[66]. The principle of DT is

74

based on the derivation of a linear relation between the spatial Fourier transform of the contrast
function and the scattered field for weak scatterers [61]-[63],[67],[68]. A generalized DT
algorithm for multi-frequency multi-monostatic Ground Penetrating Radar (GPR) measurement
configuration was first proposed by Deming and Devaney in [67]. Most of the related works on
DT were originally focused on freespace synthetic aperture radar (SAR) imaging [65],[66] and
later on subsurface imaging [67]-[69].
The information bearing signals associated with tomographic imaging are scattering,
diffraction and significant attenuation, making many conventional imaging techniques not
applicable [69]. The success of diffraction tomography techniques relies on three assumptions:
linearity, number of transmitters and receivers, and frequencies used [70]. Linearized diffraction
tomography is based on linear approximation of the scattered field such as Born- or Rytov-type
approximations [71]. Both approximations can produce excellent reconstructions for small
objects with small refractive index changes. In many cases, the linearity assumption fails;
however, for geophysical applications and other cases, the target can be considered a weak
scatterer and therefore the linear assumption holds [72], [73]. From a resolution point of view,
the best way to implement diffraction tomography is by using the highest possible number of
receiver/transmitter units, together with the highest number of frequencies available. Different
constraints, like economic, safety, operating, geometric or physical, limit the number of
receiving and/or transmitting units, leading to insufficient or incomplete data sampling which
reduces the resolution or fidelity of the reconstructed image. Although frequency is not an
intrinsic problem of diffraction tomography, it drastically limits the resolution attainable in the
reconstruction process. As presented in [71], Born and Rytov approximations are valid only if
the inhomogeneities in the medium are relatively small compared with the wavelength of the

75

incident field. On the other hand, attenuation in a lossy medium is heavily dependent on
frequency: higher frequencies yield higher losses. In some applications, such as geophysical and
biomedical imaging, the choice of the frequency is essentially based on the highest possible
value such that scattered fields can still be detected above a certain signal-to-noise ratio (SNR)
[74].
The solution is to make some necessary approximations, and usually, a Born-type
approximation is used, which gives an approximate measurement of the scattering [71]. With this
approximation, and the use of scattered waves collected by an array of receivers, the Fourier
diffraction slice theorem can be used [75]. In this method, a plane wave propagates through the
object, and scattering occurs along with diffraction either from edges, corners, creeping waves,
etc. The field is measured along a parallel line and the Fourier Transform is applied to the
measured data. The Fourier transformed data forms an arc through the origin in the frequency
domain. The theory of this method is explained in the next section.
4.3

Fourier Diffraction Theorem
The Fourier Diffraction theorem is a fundamental theorem for the comprehension of the

object reconstruction [71]. This theorem is based on the principle that objects with large
inhomogeneities when compared to a wavelength of the incident field, cannot be analyzed using
ray theory. Instead, one must resort directly to wave propagation and diffraction phenomena
[76]. The Fourier Diffraction theorem relates the Fourier transform of the measured forward
scattered field with the Fourier transform of the contrast profile of the object under test, which is
the magnitude to be imaged. The Fourier Diffraction Theorem can be stated as:
When a 2D object is illuminated by a plane wave as shown in Figure 4.1, the 2D Fourier
transform of the forward scattered field measured along a line perpendicular to the
76

direction of propagation of the wave (line TT' in Figure 4-1) gives the values of the 2D
Fourier transform of the contrast evaluated along a circular arc centered in -k0 and radius
k0 in the frequency domain [76].

Spatial Domain

Frequency Domain

ky

y

E

inc

T
-k

0

x

k

0

kx

T’ ES measured
2D FT contrast profile

Figure 4-1. Illustration of the Fourier Diffraction Theorem. When 2D object is illuminated by a
plane wave, the 2D Fourier transform (FT) of the scattered field measured along the
line TT' gives the values of the 2D Fourier transform of the contrast evaluated along
a circular arc centered in -k0 and radius k0 in the frequency domain.

Some important remarks of this theorem can be summarized as follows [69]:


Transmission and reflection measurement geometries. When measuring scattered
fields with a planar system in transmission, only the components of the angular
spectrum propagating according to the direction of incidence can be measured,
whereas the backscattering is lost. These forward components are associated with the
right half of the circle in the spectral domain. In contrast, if another measurement line
is placed in the opposite side of the object under test, or in reflection configuration,
both forward and backscattered fields can be acquired and points over the entire circle
in the spectral domain will be recovered. Accordingly, for a circular geometry, both

77

the forward and backward scattered fields can be measured and the Fourier transform
of the contrast profile will be obtained over the entire shifted circle in the spectral
domain. However, in practice, backscattered fields are difficult to measure near the
transmitting element, and the resolution is not as good as theoretically expected.


Low pass filtering. The Fourier transform of the contrast profile will be zero for the
angular spatial frequencies greater than 2k0. For angular spatial frequencies smaller
than 2k0 this magnitude will be placed along a circle in the transformed domain.



Multiview. By changing the incidence angle, the reconstruction of the contrast profile
will be obtained in shifted circles of the same radius, all with its center placed along a
circle of radius k0 centered in the origin of coordinates.



Multifrequency. By changing the illuminating frequency, circles of different radius
centered along the same direction are obtained.



Limit. Given that the radius of the circle in spectral domain is k0 

2 f
(in air),
c

when f   the circle degenerates into a line. In practice, this approximation is
valid in X-rays, which simplifies greatly the reconstruction procedure in X-ray
imaging systems.
The importance of the theorem is that if an object is illuminated by plane waves from
many directions over 360, the resulting circular arcs in the transformed-plane will fill up the
frequency domain. The image function may be recovered by Fourier inversion [69].
In this work, a 2D geometry is assumed and hence objects under investigation and
electromagnetic fields are invariant with respect to the z-axis. If vertically polarized antennas

78

relative to the z-axis direction are used, all electric field and current vectors are z-directed, and
accordingly we can use the scalar field equations [77], [78].
Let ET  r  represent the total electric field measured at one receiving antenna in
presence of the object under test, and Ei  r  the incident field. The scattered field can be defined
as:

Es  r   ET  r   Ei  r 

(4.1)

In terms of the equivalent electric current radiating in the external medium, the scattering field
can be expressed as:

Es  r   

S j0 J  r  G  r  r dr

(4.2)

where G  r  r   is the Green’s function and
J  r    j   r    0  ET  r 

(4.3)

is the equivalent induced electric current. The 2D Fourier transform of the induced electric
current can be defined as:

 

J K  F  J  r  

S J  r  e

jkr dr

(4.4)

For weakly scattering objects, Born-type approximations can be used [79], [80]. The simplest
approximation is the first-order Born approximation, which assumes that the scattered electric
field can be written only in terms of the incident field, which is known. Assume that the cylinder
is illuminated by a plane wave incident at an angle θ0, the incident electric field is

79

ˆ
Ei  r   e jk00r

(4.5)

where k0   0 0 is the wavenumber in the surrounding medium. Note that for this work, the
surrounding medium is not assumed to be free-space and hence the wavenumber and permittivity
is not that of free-space; however, the surrounding medium is assumed to be non-magnetic. The
quantity to be imaged is the contrast profile, defined as

C r   1

 r 
  r   0

0
0

(4.6)

Inserting (4.5) into (4.4), we obtain
ˆ
J  r    j   r    0  ET  r   j 0C  r  e jk00r

(4.7)

Applying the 2D transform to (4.7),

 

J K  F  J  r   j 0

 j 0



S C  r  e

 S C  r  e 

ˆ
where C k  k00 

ˆ
j k  k00

S

ˆ
C  r  e  jK00 r e jkr dr



ˆ
j k  k00

dr

dr  j 0C  k  k0ˆ0 

(4.8)

is the 2D Fourier transform of the contrast

ˆ
profile evaluated over a circle of radius k0 and center k00 in the Fourier space, as illustrated
in Figure 4-2.

80

k0

Figure 4-2. Fourier domain for an incident plane wave. Under Born approximation, the 2D
Fourier transform of the contrast profile is obtained in the transformed domain over
a circumference shifted on the direction opposite to the incident plane wave.

The reconstruction of (4.8) gives information of the 2D Fourier transform of the contrast
profile only over a circumference of radius k0 shifted a distance of k0 in the direction opposite to
the propagation of the incident field. In order to obtain the whole reconstruction of the object,
different directions of incidence must be considered, as illustrated in Figure 4-3. In other words,
by transmitting from one of the antennas and receiving from the rest and repeating this procedure
for all transmitters, a two-dimensional sampling inside a circle of radius 2k0 in the spectral
domain will be obtained. Accordingly, the reconstructed contrast profile calculated through the

 

inversion of C k will be a low-pass filtered version of the original one.

81

-k
-2k

k

0

0

0

2k0

Figure 4-3. Multiview Fourier domain. Different directions of illumination are necessary to
obtain the reconstruction of the object. The 2D Fourier transform of the contrast
profile is obtained inside a circle of radius 2k0 centered at the origin.
Two configurations are commonly used for obtaining different directions of illumination:
planar and circular arrays. In this work, a circular array, like the one shown in Figure 4-4, is
used. For two-dimensional fields, the incident field is a cylindrical wave. The spectral
formulation for planar systems can be applied in circular geometries, if a synthetic aperture
approach is made. Specifically, to create a plane wave by superposing the existing cylindrical
waves created by each source located along the circular antenna array.
Let’s assume that the incident field in Figure 4-4 is a plane wave defined as
ˆ
ˆ
Ei r ,0  e jk00r . The electric equivalent current induced in the cylinder by the incident





field is J b . The current J b radiates the scattered field Es that is measured at the receivers.

82

y



Rx

x

Tx



Figure 4-4. Circular array of antennas. Different directions of illumination can be obtained by
transmitting with each one of the elements. The transmitter is Tx and all others are
the receivers Rx.
The 2D Fourier transform of the induced current is

  v J a  Esdva
a

ˆ
Jb k0 

(4.9)

where J a is the current in the transmitter, va is the circular path defined by the array, and
ˆ
ˆ
Jb k0  Jbe jk0 r dr

(4.10)

ˆ
Jb  r   j 0C  r  e jk00r

(4.11)

  

Using (4.7), we obtain

83

 

ˆ
Then, Jb k0 can be expressed in terms of the contrast profile

 

ˆ
ˆ
ˆ
ˆ ˆ
Jb k0  j 0 C  r  e jk00r  e jk00r dr  j 0C k0   0

 





(4.12)

For a circular antenna array of radius R centered in the xy-plane

 

ˆ ˆ
j 0C k0   0

2

  0





ˆ
J a  ,  Es  ,0 Rd

(4.13)

It can be demonstrated that the current producing a plane wave propagating in the ˆ
direction can be expressed as the summation of cylindrical modes [81]:





ˆ
I   



2
0 R





ˆ
jn  
jn
e
(2)
n  H n  k0 R 





(4.14)



ˆ
where I    is the amplitude of the current distribution of the transmitter at an angular
position  along the antenna. Equation (4.13) can be written as,

 

ˆ ˆ
j 0C k0   0

2

  0 I   ˆ  Es  ,ˆ0  Rd

(4.15)

ˆ
ˆ
When the incident field is a cylindrical wave, the plane wave Ei r ,0  e jk00r can





be expanded as the superposition of cylindrical waves Ei  r ,   generated by line sources





ˆ
I   0 located at an angular position  in the array.



2

 0 I   ˆ0  Ei  r ,  Rd

ˆ
Ei r ,0 

84

(4.16)

The corresponding scattered field can be expressed as





ˆ
where Es  ,0



2

 0 I   ˆ0  Es  ,  Rd

ˆ
Es  ,0 

(4.17)

is the scattered field measured at position  generated by an incident





ˆ
cylindrical wave produced by a line source located at  and I   0 is





ˆ
I   0 

2
0 R





ˆ
jm  0
jm
e
(2)
m  H m  k0 R 





(4.18)

Using (4.15) and (4.17), the 2D Fourier transform of the contrast profile is

 

ˆ ˆ
C k0   0

  j0
1







  Es  ,  I   ˆ  I   ˆ0  d d

(4.19)

0 0





2 2
R2



ˆ
ˆ
The term I    I   0 in equation (4.19) can be defined as

ˆ
ˆ
I ( ,  )  I    I   0
 2 


 0 R 

4.4

2











ˆ
ˆ
jn  
jm  0
jn
jm
e
e
(2)
(2)
H n  k0 R 
H m  k0 R 
n  m 

 



(4.20)

Reconstruction Algorithm
The circular array consists of N antennas in the xy-plane. As noted above, a sample of the

scattering field Es  ,   can be obtained by transmitting alternatively with each antenna and
receiving each time with the other N-1 antennas.

85

Equation (4-19) can be seen as a 2D convolution between the scattered field and the
current that generates a plane wave existing along the circular array. This convolution can be
calculated as the inverse Fourier transform of the product between the corresponding spectra.





Es  ,    I  ,    F 1 Es  ,    I  ,  

(4.21)

For an array of transmitting and receiving antennas placed along a circle of radius R, the
spectrum of the scattered field can be calculated as

 
Es  r , t  

  E s  ,   e

 jr e jt

 
 


 



 





2
k0 I

4 0
2
k0 I

4 0

2
k0 I

4 0




n 




n 





n 

2
(2)
an  H n  k0 R   e jn(  )e  jr e  jt d d





(2)
an  H n  k0 R  





2  

 e

j (n  r ) e  j (n t ) d d

 

2
(2)
an  H n  k0 R     n  r    n  t 





 k 2I
2
 0 a  H (2)  k R  
 n

0
Es  r , t    4 t 

0

0


t  r 

(4.22)

otherwise

where an is defined in (3.19). The Fourier transform of the current I ( ,  ) defined in (4.20) can
be calculated as follows. Rearranging (4.20) we obtain,

86

 2 
I ( ,  )  

 0 R 



2

j



 

n m 

(2)
(2)
n  m  H n  k0 R  H m  k0 R 

e jn e jm

(4.23)

The corresponding Fourier transform of (4.23) is

 2 
I (r , t )  

 0 R 

2

 2 


 0 R 

2





 

j

 

n m 

 e

(2)
(2)
n  m  H n  k0 R  H m  k0 R   




 

j

n m 

(2)
(2)
n  m  H n  k0 R  H m  k0 R 


 2

I (r , t )   0 R




2
r t 

j

(2)
 H r  k0 R  H t(2)  k0 R 
0

jn e jm e jr e jt d d

  n  r  m  t 

t, r 

(4.24)

otherwise

Therefore, the discrete version of equation (4.19) is

 

ˆ ˆ
C k0   0





R 2  2

j 0  N

 1
Es  ,    I  ,  
F






(4.25)

where F 1 is the inverse Fourier transform. The contrast profile can be obtained by computing
the inverse Fourier transform of (4.25). Defining

2
k0 I
A
4 0
and

87

(4.26)

 2 
B

 0 R 

2
(4.27)

The product of the spectra Es  r , t   I  r , t  can be calculated as

2
 H (2) k R 
 t  0 


f  r , t   Es  r , t   I  r , t 
 ABat
t r
(2)
(2)
H r  k0 R  H t  k0 R 
(2)

H
k R
 ABa t  0 

t (2)
f  r, t   
H t  k0 R 

0



(4.28)

t  r 
otherwise

Calculating the inverse Fourier transform of (4.28), we obtain



f  ,   



 



f  r , t  e jr e jt

r  t 

 AB

  1t at e j  t

(4.29)

t 

t r

The contrast profile is

 

ˆ ˆ
C k0   0





1
j 0

R 2 f  ,0  



jI

 k0 

2



 1t at e





ˆ ˆ
j  0 t

(4.30)

t 

For a given frequency, a sampled version of the spectrum of the contrast profile inside a
circle of radius 2k0 is obtained. In order to improve the quality of the image, a set of multiple
frequencies can be combined [82]. Previously reported results lead to the conclusion that
configurations with large contrast with respect to the surrounding medium can be reconstructed
with an acceptable resolution by using multiple-frequency reconstruction [83][84][72]. At each

88

new frequency, the final result for the previous frequency is used to estimate the unknown
permittivity profile. To avoid destructive frequency-dependent interferences, the average of the
magnitudes of the images is used. The total contrast is given by

C

f max



f  f min

 

F 1 Ck
0

(4.31)

where F 1 is the inverse Fourier transform. Using average magnitudes produces some loss of
quantitative dielectric properties information, but the size and location is preserved.
4.5

Results
The algorithm implementation was done using Matlab. The starting point of the

algorithm is the scattered field measurement obtained from simulations of the cylinders
presented in Chapter 3. This data is organized as a third-order matrix of dimensions
NTNRNf, where the number of transmitting and receiving antennas is the same (NT = NR)
and Nf is the number of frequencies acquired. In practice, the number of antennas is limited as
well as the number of frequencies. The following results are for a frequency range of 300 to 800
MHz for three separate antenna configurations: 16, 32, and 64 antennas equally-spaced in an
angular sense forming a circular array. Simulations were performed for concentric and eccentric
spill configurations as presented in Chapter 3.
For the centered configuration, the spill has a relative permittivity of 1  12 in the outer
layer, corresponding to contaminated soil with 25% water content and 25% oil content, and

 2  2.2 in the core region, corresponding to a fully contaminated soil. The surrounding medium
has a relative permittivity of 17, corresponding to wet soil. The contrast of the cylinders is
89

constant and purely real. Using (4.6) for the large cylinder, representing the outer region of
contaminated soil, the contrast is 0.294. For the small cylinder, or core region, the contrast is
0.8167. The original contrast profile and reconstructed results are shown in Figures 4-5 and 4-6.

(a)

(b)

(c)

(d)

Figure 4-5. Contrast profile for a centered spill configuration. The contrast of the cylinders is
constant and purely real. For the large cylinder the contrast is 0.294. For the small
cylinder, or core region, the contrast is 0.8167. (a) Original image. (b)
Reconstructed profile with 16 antennas. (c) Reconstructed profile with 32 antennas.
(d) Reconstructed profile with 64 antennas.

90

(a)

(b)

(c)

(d)

Figure 4-6. Contrast profile for a centered spill configuration. For the large cylinder the contrast
is 0.294. For the small cylinder, or core region, the contrast is 0.8167. (a) Original
image contrast profile. (b) Reconstructed profile with 16 antennas. (c)
Reconstructed profile with 32 antennas. (d) Reconstructed profile with 64 antennas.
Results for the centered configuration showed that the algorithm can detect the spill with
a poor resolution. The contrast profile is partially reconstructed, with expected contrast values at
the boundary of the outer cylinder. The core cylinder is detected with a high contrast profile at
the center and a contrast of almost zero at the boundary. As the number of antennas increases,
the boundaries are better defined for both cylinders. Notice in Figure 4-6 that the surrounding

91

medium shows a contrast profile of approximately 0.18 instead of zero. The size of the spill is
preserved.
For the eccentric case, the same contrast profile was used. Results are shown in Figures
4-7 and 4-8.

(a)

(b)

(c)

(d)

Figure 4-7. Contrast profile for an eccentric spill configuration. The contrast of the cylinders is
constant and purely real. For the large cylinder the contrast is 0.294. For the small
cylinder, or core region, the contrast is 0.8167. (a) Original image. (b)
Reconstructed profile with 16 antennas. (c) Reconstructed profile with 32 antennas.
(d) Reconstructed profile with 64 antennas.

92

(a)

(b)

(c)

(d)

Figure 4-8. Contrast profile for a centered spill configuration. The contrast of the cylinders is
constant and purely real. For the large cylinder the contrast is 0.294. For the small
cylinder, or core region, the contrast is 0.8167. (a) Original image. (b)
Reconstructed profile with 16 antennas. (c) Reconstructed profile with 32 antennas.
(d) Reconstructed profile with 64 antennas.

Reconstructed profile for the eccentric case showed that the image is affected by the
position of the core cylinder. The outer cylinder is blended with the core cylinder and the
contrast profile is not obtained. Hence, if the material within the resolution cell is a blend of the
two materials, the resulting image is blurred. This is to be expected with any non-parametric

93

imaging algorithm. Parametric algorithms, which are beyond the scope of this research, utilize a
priori information to sharpen the image. Such information involving the permittivity profile is
generally not available for the remote sensing scenario under investigation herein. Also, the
reconstructed image is shifted to the side opposite to the position of the core. In other words, the
reconstructed image is centered with the core cylinder. As before, the surrounding medium
appears with a contrast of 0.18.
In both cases, the reconstructed images conserve the size of the spill and the boundaries
are well defined. At the boundary of the outer cylinder, the contrast profile obtained with 64
antennas, has an average value of 0.2798. This is close to the original contrast of 0.294. The
maximum value obtained at the core cylinder for the same number of antennas is 0.849.
The inversion method presented in this chapter can be used for detection of spills based
on the contrast between contaminated soils. The results of the simulation model are, at least from
a geometrical standpoint, acceptable. The algorithm is able to determine the position of the spill
and the core cylinder. The contrast and therefore, permittivities of the soil and contaminated soil
can be determined by the color variance in the images.

94

CHAPTER 5
CONCLUSIONS AND FUTURE WORK

The primary goal of this work was to simulate a contaminated region and determine its
physical limits within a tolerable level of certainty given the characteristics of the soil in question
without the environmental and financial cost of physical experimentation. To achieve this, the
first part of this study focused in the determination of the dielectric properties of soil and
contaminated soil. A coaxial waveguide for measuring the complex permittivity and permeability
of soil and contaminated soil was designed and tested successfully. Different from existing
sample holders, the new sample holder was fabricated to accommodate significant volumetric
samples of material. This type of sample holder is attractive for large measurement programs
were larger samples are needed.
Air and Plexiglas samples were measured to verify the reliability and accuracy of the
measurement system. Results obtained for air, Plexiglas, and the dry soil samples were as
expected. The increase in relative permittivity for the soil with moisture was reasonable.
The measured contrast between dry soil and dry soil contaminated with oil was within the
nominal uncertainty of the experiment and hence it is not expected that a rise in permittivity for
dry soil contaminated with fairly high volumetric concentrations of oil will allow sufficient
contrast for either radio frequency bore-hole or GPR mapping. There was however a significant
difference between the measured permittivity for wet soil and that of wet soil contaminated with
oil as shown in Table 5-1. Hence, it is reasonable that this permittivity contrast can be used as a
discriminator for detecting contamination plumes. This is significant since contamination in an
aquifer is one of the most environmentally challenging issues facing the world.

95

Table 5-1. Percent difference for permittivity values of wet- and wet-contaminated soils.
Real relative permittivity
(Average)

Moisture Content
10%

25%

Wet soil

6.9148

17.2000

Wet–contaminated soil

5.5271

11.7855

% DIFFERENCE

22.3%

37.36%

The ground-truth obtained in these experiments was used to simulate contaminants in soil
and determine the scattering characteristics of NAPL’s spills. The problem is two-dimensional
and the effects of various geometrical and electrical parameters were examined. Results for plane
wave and line source excitations for homogeneous and inhomogeneous spills were obtained.
Scattered field matrices obtained was used as input parameter in the inversion algorithm.
The inversion algorithm is based in the Fourier Diffraction theorem and, in order to improve the
quality of the image, a set of multiple frequencies are combined. To avoid destructive frequencydependent interferences, the average of the magnitudes of the images was used. Using average
magnitudes produces some loss of quantitative dielectric properties information, but the size and
location is preserved. The results obtained for concentric and eccentric spills showed that the
reconstructed images conserved the size of the spill and the boundaries were well defined. At the
boundary of the outer cylinder, the contrast profile obtained with 64 antennas, has an average
value of 0.2798. This is close to the original contrast of 0.294. The maximum value obtained at
the core cylinder for the same number of antennas is 0.849, close to the original contrast of
0.817.

96

The inversion method presented can be used for detection of spills based on the contrast
between contaminated soils. The results of the simulation model were, at least from a
geometrical standpoint, acceptable. The algorithm was able to determine the position of the spill
and the core cylinder. The contrast and therefore, permittivities of the soil and contaminated soil
can be determined by the color variance in the images.
There are several directions in which this work may be extended. The influence of
variations in the background medium must be further investigated in order to improve the
accuracy of the system. Larger spills with respect to wavelength should be analyzed. Also, the
effects of increasing bandwidth without increasing the number of transmitters/receivers should
be analyzed. Concerning experimental measurements, it would be interesting to test the
algorithm with data form a real oil spill site instead of simulated values. Also, the design of the
borehole and antenna system is worth looking into. Low frequencies need to be used and the
system needs to be low-powered, because of the possible explosive hazard of gases situated in
the subsurface. Also, the electronic components of the system must be able to withstand the
extreme temperatures and the moisture in the boreholes.
The work presented could have useful applications in the geophysical and biomedical
areas. Future work will be directed in those areas.

97

BIBLIOGRAPHY

98

BIBLIOGRAPHY

[1]

A.M. Chrestman, G.D. Comes, S.S. Kooper, and P.G. Malone, “Rapid detection of
hydrocarbon contamination in ground water and soil,” Proceedings of the Hydraulic
Engineering, Water Forum ’92, Baltimore, Maryland, August 2–6, pp. 1165-1170, 1992.

[2]

United States Environmental Protection Agency. (2009, July 31). Drinking Water
Contaminants [Online]. Available: http://water.epa.gov/drink/contaminants

[3]

United States Environmental Protection Agency. (2005, March 23). Cost and Performance
Report for LNAPL Recovery [Online]. Available:
http://www.epa.gov/tio/download/rtdf/napl/cpbpsugarcreek.pdf

[4]

Remediation Technologies Development Forum. (2005, March 21). "The Basics"
Understanding the Behavior of Light Non-Aqueous Phase Liquids (LNAPLs) in the
Subsurface [Online]. Available: http://www.rtdf.org/public/napl/training/module1.pdf

[5]

J. Ryan. (2009). Leaking underground storage tanks [Online]. Available:
http://bcn.boulder.co.us/basin/waterworks/lust.html

[6]

F. Nadim, G. E. Hoag, S. Liu. R. J. Carley, and P. Zack, “Detection and remediation of soil
and aquifer systems contaminated with petroleum products: an overview,” Journal of
Petroleum Science and Engineering, vol. 26, pp. 168-178, 2000.

[7]

C. J. Newell, S. D. Acree, R. R. Ross, and S. G. Huling, “Light nonaqueous phase liquids,”
Ground Water Issue, EPA/540/S-95/500, U.S.EPA, R.S. Kerr Environ. Res. Lab., Ada,
OK, pp. 1-28, 1995.

[8]

A.J. Devaney, “Geophysical diffraction tomography,” IEEE Trans. Geosc.Remote Sensing,
vol. GE-22, no. 1, pp. 3-13, Jan.1984.

[9]

A.G. Tijhuis, K. Belkebir, A.C.S. Litman, and B.P. deHon, “Multiple-frequency distortedwave Born approach to 2D inverse profiling,” Inverse Problems, vol. 17, pp. 1635-1644,
2001.

[10] Standard Guide for Using the Surface Ground Penetrating Radar Method for Subsurface
Investigation, ASTM Standard D 6432 – 99 (Reapproved 2005).
[11] Timur, A. and Toksöz, M. N., “Downhole geophysical logging,” Annual Review of Earth
and Planetary Sciences, vol. 13, pp.315-344, 1985.
[12] J. Baker-Jarvis, M.D. Janezic, B.F. Riddle, R.T. Johnk, P. Kabos, C.L. Holloway, R.G.
Geyer, and C.A. Grosvenor, “Measuring the permittivity and permeability of lossy
materials: solids, liquids, metals, building materials, and negative-index materials,” NIST
Tech. Note 1536, Boulder, CO, December 2004.

99

[13] J.O. Curtis, “A durable laboratory apparatus for the measurement of soil dielectric
properties,” IEEE Trans. Instrum. Meas., vol. 50, no. 5, pp. 1364-1369, Oct. 2001.
[14] A.M. Nicholson and G.F. Ross, “Measurement of the intrinsic properties of materials by
time-domain techniques,” IEEE Trans. Instrum. Meas., vol. IM-19, pp. 377-382, Nov.
1970.
[15] W.B. Weir, “Automatic measurement of complex dielectric constant and permeability at
microwave frequencies,” Proc. IEEE, vol. 62, no.1, pp. 33-36, Jan. 1974.
[16] J. Baker-Jarvis, E.J. Vanzura, and W.A. Kissick, “Improved technique for determining
complex permittivity with the Transmission/Reflection method,” IEEE Trans. Microwave
Theory Tech., vol. 38, no. 8, pp.1096-1103, August 1990.
[17] A.H. Boughrie, C. Legrand, A. Chapoton, “Non-iterative stable transmission/reflection
method for low-loss material complex permittivity determination”, IEEE Trans.
Microwave Theory Tech., vol. 45, no.1, pp. 52-57, January 1997.
[18] C.C. Courtney, “Time-domain measurement of the electromagnetic properties of
materials,” IEEE Trans. Microwave Theory Tech., vol. 46, no. 5, pp. 517-522, May 1998.
[19] D.M. Pozar, Microwave Engineering, New York: Wiley, 2005.
[20] HVS Technologies, Inc., 2597-2 Clyde Ave., State College, PA 16801.
[21] G.F. Engen and C.A. Hoer, “Thru-Reflection-Line: An improved technique for calibrating
the dual six-port automatic network analyzer,” IEEE Trans. Microwave Theory Tech., vol.
MTT-27, no. 12, pp. 987-993, Dec. 1979.
[22] A. R. Von Hippel, Dielectric Materials and Applications, New York: Wiley, 1954.
[23] USDA Natural Resources and Conservation Service Web Soil Survey. (2009, November
11). Report OH119 Muskingum County OH [Online]. Available:
http://websoilsurvey.nrcs.usda.gov/app/
[24] American Association of State Highway and Transportation Officials Standards and
Guidelines. (2004). T-265: Laboratory Determination of Moisture of Soils [Online].
Available: http://www.dot.nd.gov/manuals/materials/testingmanual/t265.pdf
[25] J.H. Dane and G.C. Topp, Methods of Soil Analysis, Part 4 – Physical Methods, Soil
Science Society of America, Inc., Madison, WI, 2002.
[26] N.R. Peplinski, F.T. Ulaby, and M.C. Dobson, “Dielectric properties of soils in the 0.3 –
1.3 GHz range,” IEEE Trans. Geosc. Remote Sensing, vol. 33, no. 3, pp. 803-807, May
1995.

100

[27] M.C. Dobson, F.T. Ulaby, M. Hallikainen, and M. El-Rayes, “Microwave dielectric
behavior of wet soil – Part II: four-component dielectric mixing models,” IEEE Trans.
Geosc. Remote Sensing, vol. GE-23, no. 1, pp. 35-41, January 1985.
[28] F.T. Ulaby, R.K. Moore, and A.K. Fung, Microwave Remote Sensing, vol. 3, Dedham,
MA: Artech House, 1986, Appendix E.
[29] A. Sihvola, Electromagnetic Mixing Formulas and Applications, IEE Electromagnetic
Waves Series 47, Michael Faraday House, UK, 1999.
[30] A. Sihvola, “Mixing rules with complex dielectric coefficients,” Subsurface Sensing
Technologies and Applications, vol. 1, no. 4, pp. 393-415, October 2000.
[31] W.R Tinga, W.A.G Voss, D.F. Blossey, “Generalized approach to multiphase dielectric
mixture theory,” J. Appl. Phys. vol. 44 no. 9, pp. 3897, 1973.
[32] J.C. Maxwell Garnett, “Colours in metal glasses and metal films,” Trans. of the Royal
Society, vol. CCIII, pp. 385–420, 1904.
[33] D. Polder and J.H. van Santen, “The effective permeability of mixtures of solids”, Physica,
vol. 12, no. 5, p. 257–271, 1946.
[34] D.A.G. Bruggeman, “Berechnung verschiedener physikalischer Konstanten von
heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper
aus isotropen Substanzen,” Annalen der Physik, vol. 416, no. 8, pp. 665–679, 1935.
[35] L. Tsang, J.A. Kong, and R.T. Shin, Theory of microwave remote sensing, Wiley New
York, 1985.
[36] U. Kaatze, “Complex permittivity of water as a function of frequency and temperature,” J.
Chem. Eng. Data, vol. 34, no. 4, pp. 371-374, October 1989.
[37] W.J. Ellison, “Permittivity of Pure Water, at Standard Atmospheric Pressure, over the
Frequency Range 0–25 THz and the Temperature Range 0–100 °C,” J. Phys. Chem. Ref.
Data, vol. 36, pp. 1-18, 2007.
[38] W.M. Alley, T.E. Reilly, and O.L. Franke, “Sustainability of Ground-Water Resources,”
U.S. Geological Survey Circular 1186, Denver, Colorado, 1999.
[39] F.M. Francisca, and V.A. Rinaldi, “Complex dielectric permittivity of soil-organic
mixtures (20 MHz-1.3 GHz)”, Journal of Environmental Engineering, vol. 129, no. 4, pp.
347-357, April 2003.
[40] B.B.S. Singhal and R.P. Gupta, “Groundwater Contamination” in Applied Hydrogeology
of Fractured Rocks, Springer Netherlands, 2010, ch.12, pp. 221-236.

101

[41] A. Sihvola, “Model systems for materials with high dielectric losses in aquametry,” in
Electromagnetic Aquametry: Electromagnetic Wave Interaction with Water and Moist
Substances, K. Kupfer Ed., Springer-Verlag Berlin Heidelberg, 2005, ch. 5, pp. 93-111.
[42] C.M. Rappaport, “Ground penetrating radar,” in The RF and Microwave Handbook,
Volume 1, Mike Goglio Ed., CRC Press, Boca Raton, FL, 2008, ch. 17.
[43] A. Ramírez, W. Daily, A. Binley, and D. Roelant, “Detection of leaks in underground
storage tanks using electrical resistivity methods,” Journal of Environmental and
Engineering Geophysics, vol. 1, no. 3, pp. 189-203, December, 1996.
[44] D. Eisenburger and V. Gundelach, “A new direction-sensitive borehole logging tool for
the spatial reconnaissance of geological structures,” Proceedings of the 12th International
Conference on Ground Penetrating Radar, Birmingham, UK, June 16-19, 2008.
[45] S.E. Prensky, “Advances in borehole imaging technology and applications,” in Borehole
Imaging: Applications and Case Histories, Lowell, Williamson, and Harvey Eds.,
Geological Society Special Publication No. 159, UK, 1999, pp. 1-44.
[46] A.A. Kishk, R.P. Parrikar, and A.Z. Elsherbeni, “Electromagnetic Scattering from an
Eccentric Multilayered Circular Cylinder,” IEEE Trans. Antennas and Prop., vol.40, no. 3,
pp. 295-303, March 1992.
[47] S.J. Radzevicius and J.J. Daniels, “Ground penetrating radar polarization and scattering
from cylinders,” Journal of Applied Geophysics, vol. 45, no. 2, pp. 111-125, September
2000.
[48] R.F. Harrington, Time-Harmonic Electromagnetic Fields, New York: McGraw-Hill, 1961.
[49] C. Balanis, Advanced Engineering Electromagnetics, New York: Wiley, 1989.
[50] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, New York, Dover
Publications, 1972.
[51] G.T. Ruck, D.E. Barrick, W.D. Stuart, and C.K. Krichbaum, Radar Cross Section
Handbook, New York: Plenum, 1970.
[52] A.P. Annan, J.L. Davis, “Ground-penetrating radar for high-resolution mapping of soil and
rock stratigraphy,” Geophysical Prospecting, vol. 37, no. 5, pp. 531-551, 1989.
[53] J.A. Stratton, Electromagnetic Theory, New York: McGraw-Hill, 1941.
[54] S. Lambot, I. Van den Bosch, B. Stockbroeckx, P. Druyts, M. Vanclooster, E.C. Slob,
“Frequency dependence of the soil dielectric properties derived from ground penetrating
radar signal inversion,” Subsurface Sensing Technologies and Applications, vol.6, no. 1,
pp. 73-87, 2005.

102

[55] S. Lambot, E. Slob, J. Minet, K.Z. Jadoon, M. Vanclooster, and H. Vereecken, “Fullwaveform modeling and inversion of ground-penetrating radar data for non-invasive
characterization of soil hydrogeophysical properties,” in Proximal Soil Sensing, Progress
in Soil Science I, R.A.Viscarra-Rossel, A.B. McBratney, B. Minasny, Eds., Springer
Science and Business Media, 2010, ch. 25, pp. 299-311.
[56] C.M. Rappaport, “Ground Penetrating Radar,” in RF and Microwave Applications and
Systems, M. Goglio and J. Goglio Eds., CRC Press, 2nd Ed., 2008, ch. 17.
[57] M. Farid, A. N. Alshawabkeh, and C. M. Rappaport, “Electromagnetic waves in
contaminated soils,” in Electromagnetic Waves in Complex Matter, Ahmed Kishk Ed.,
InTech, 2011, ch. 5, pp. 117-154.
[58] M. Farid, A. N. Alshawabkeh, H. Zhan, and C. M. Rappaport, “Challenges and Validation
of Cross Tomography Experimentation for Inverse Scattering Problems in Soils,” in Proc.
of Geo-Frontiers 2005: GSP 138, Site Characterization and Modeling, 2005 © ASCE
Conf. Proc. doi: 10.1061/40785(164)39.
[59] H. Zhan, C. M. Rappaport, M. Farid, A. N. Alshawabkeh, and H. Raemer, “Lossy halfspace Born approximation modeling of electromagnetic wave source and scattering by soil
by cross well radar,” Proc. of Geo-Frontiers 2005: GSP 138, Site Characterization and
Modeling, 2005 © ASCE Conf. Proc. doi: 10.1061/40785(164)38.
[60] E.Wolf, “Three-dimensional structure determination of semi-transparent objects from
holography data,” Optics Communications, vol. 1, pp.153-156, 1969.
[61] T.J.Cui, and W.C. Chew, “Novel diffraction tomographic algorithm for imaging twodimensional dielectric objects buried under a lossy earth,” IEEE Trans. Geosci. Remote
Sensing, vol. 38, no. 4, pp. 2033-2041, July 2000.
[62] T.J. Cui, and W.C. Chew, “Diffraction tomographic algorithm for the detection of threedimensional objects buried in a lossy half space,” IEEE Trans. Antennas and Propagation,
vol. 50, no. 1, pp. 42-49, Jan. 2002.
[63] S.K. Lehman, “Superresolution planar diffraction tomography through evanescent fields,”
Int. J. Imaging Systems Technol., vol. 12, pp. 16-26, 2002.
[64] W.C. Chew, Waves and Fields in Inhomogeneous Media, chapter 2, IEEE Press,
Piscataway, NJ, 1997.
[65] M. Soumekh, “A system model and inversion for synthetic aperture radar imaging,” IEEE
Trans. Image Processing, vol. 1, no. 1, pp. 64-76, Jan. 1992.
[66] J.M. Lopez-Sanchez and J. Fortuny-Guasch, “3-D radar imaging using range migration
techniques,” IEEE Trans. Antennas and Propagation, vol. 48, no. 5, pp. 728-737, May
2000.

103

[67] R. Deming, and A.J. Devaney, “Diffraction tomography for multi-monostatic ground
penetrating radar imaging," Inverse Problems, vol. 13, pp. 29-45, 1997.
[68] T. Hansen, B. and P.M. Johansen, “Inversion scheme for monostatic ground penetrating
radar that takes into account the planar air-soil interface,” IEEE Trans. Geosci. Remote
Sensing, vol. 38, no. 1, pp. 496-506, Jan. 2000.
[69] A. Kak and M. Slaney, “Tomographic Imaging with Diffracting Sources,” in Principles of
Computarized Tomographic Imaging, SIAM Classics in Applied Mathematics, ch. 6, pp.
203-273, 2001.
[70] L. LoMonte, A.M. Bagci, D. Erricolo, and R. Ansari, “Spatial resolution in tomographic
imaging with diffracted fields,” Proc. URSI General Assembly 2008 [Online] Available:
http://www.ursi.org/proceedings/procGA08/papers/B07p2.pdf.
[71] M. Slaney, A. Kak, and L. Larsen, “Limitations of imaging with first-order diffraction
tomography,” IEEE Transactions on Microwave Theory and Techniques, vol. 32, pp. 860874, 1984.
[72] A.G. Tijhuis, K. Belkebir, A.C.S. Litman, and B.P. deHon, “Theoretical and computational
aspects of 2-D inverse profiling,” IEEE Trans. Geosci. Remote Sensing, vol. 39, no. 6, pp.
1316-1330, June 2001.
[73] G.A. Tsihrintziz and A.J. Devaney, “Maximum likelihood estimation of object location in
diffraction tomography, part II: strongly scattering objects,” IEEE Trans. Signal
Processing, vol.39, no. 6, pp. 1466-1470, June 1991.
[74] M. Pastorino, Microwave Imaging, New Jersey: John Wiley & Sons, 2010.
[75] K.W.Suh., S.Y. Kim, and J.W. Ra, “Electromagnetic reconstruction of 2D permittivity
profiles of inhomogeneous dielectric cylinder by the diffraction slice-projection
algorithm”, IEEE Transactions on Magnetics, vol. 29, no. 2, pp. 1799-1802, March 1993.
[76] S.X. Pan and A.C. Kak, “A computational study of reconstruction algorithms for
diffraction tomography: interpolation versus filtered backpropagation,” IEEE Trans. On
Acoustics, Speech, and Signal Proc., vol. ASSP-31, no. 5, pp. 1262-1275, Oct. 1983.
[77] J. Rius, C. Pichot, L. Jofre, J. Bolomey, N. Joachimowicz, A. Broquetas, and M. Ferrando,
“Planar and cylindrical active microwave temperature imaging: numerical simulations,”
IEEE Transactions on Medical Imaging, vol. 11, pp. 457-459, 1992.
[78] J. M. Rius, M. Ferrando, L. Jofre, E. de los Reyes, A. Elias-Fuste, and A. Broquetas,
“Microwave tomography: an algorithm for cylindrical geometries,” Electronics Letters,
vol. 23, pp. 564-565, 1987.
[79] P.M. Morse and M. Feshbach, Methods of Theoretical Physics, McGraw-Hill, New York,
1953.

104

[80] M. Born and E. Wolf, Principles of Optics, Pergamon Press, New York, 1965.
[81] A. Broquetas, “Tomografía de Microondas en geometria cilíndrica para aplicaciones
biomédicas,” Ph.D. Dissertation, Universitat Politècnica de Catalunya, Barcelona, Spain,
1989.
[82] M.G. García, “UWB Tomography for Breast Tumor Detection,” M.S. Thesis, Dept. of
Signal Theory and Communications, Universitat Politècnica deCatalunya, Barcelona,
Spain, 2009.
[83] W.C. Chew and J.H. Lin, “A frequency-hopping approach for microwave imaging of large
inhomogeneous bodies,” IEEE Microwave and Guided Wave Letters, vol. 15, pp. 439-441,
1995.
[84] Y. Chen, “Inverse scattering via Heisenberg’s uncertainty principle,” Inverse Problems,
vol. 13, pp. 253-282, 1997.

105