close

Вход

Забыли?

вход по аккаунту

?

Broadband microwave inverse scattering: Theory and experiment

код для вставкиСкачать
INFORMATION TO USERS
This manuscript has been reproduced from the microfilm master. UMI
films the text directly from the original or copy submitted. Thus, some
thesis and dissertation copies are in typewriter face, while others may
be from any type of computer printer.
The quality of this reproduction is dependent upon the quality of the
copy submitted. Broken or indistinct print, colored or poor quality
illustrations and photographs, print bleedthrough, substandard margins,
and improper alignment can adversely affect reproduction.
In the unlikely event that the author did not send UMI a complete
manuscript and there are missing pages, these will be noted. Also, if
unauthorized copyright material had to be removed, a note will indicate
the deletion.
Oversize materials (e.g., maps, drawings, charts) are reproduced by
sectioning the original, beginning at the upper left-hand comer and
continuing from left to right in equal sections with small overlaps. Each
original is also photographed in one exposure and is included in
reduced form at the back of the book.
Photographs included in the original manuscript have been reproduced
xerographically in this copy. Higher quality 6” x 9" black and white
photographic prints are available for any photographs or illustrations
appearing in this copy for an additional charge. Contact UMI directly
to order.
UMI
University Microfilms International
A Bell & Howell Information Company
300 North Z eeb Road. Ann Arbor. Ml 48106-1346 USA
313/761-4700 800/521-0600
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
O rder N u m b er 9503847
B roadband m icrow ave inverse scattering: T h eory and
experim ent
Weedon, William Harold, Ph.D.
University of Illinois at Urbana-Champaign, 1994
UMI
300 N. Zeeb Rd.
Ann Arbor, MI 4S106
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
B R O A D B A N D MICROWAVE IN V ER SE SCATTERING :
THEORY A N D EX PE R IM E N T
BY
W IL L IA M H A R O L D W E E D O N
B .S ., N o rth ea stern U niversity, 1989
M .S ., N o rth ea stern U n iversity, 1990
T H E SIS
S u b m itte d in partial fulfillm ent o f th e req uirem en ts
for th e d egree o f D o cto r o f P h ilosop h y in E lectrical E n gin eerin g
in th e G rad u ate C ollege o f th e
U n iv ersity o f Illinois a t U rbana-C ham paign, 1994
U rbana, Illinois
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN
t
THE GRADUATE COLLEGE
M A Y 1994
W E HEREBY RECOMMEND TH AT THE THESIS BY
W ILLIA M H A R O LD W E E D O N
fmtttt
r n B R O A D B A N D M ICROW AVE IN V E R S E S C A T T E R IN G :
T H E O R Y A N D E X P E R IM E N T
BE ACCEPTED IN PARTIAL FULFILLM ENT OF TH E REQUIREMENTS FOR
D O C T O R OF P H IL O S O P H Y
TH E DEGREE O F.
Vi . c
Director df Thesis Research
Cl
Head oFDepartment
oflDi
Committee on Final Examination!
. C . c A t^ J
Chairperson
rfgCJLuLanj.
— ,__________
t Required for doctor’s degree but not for master’s.
0-517
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
© Copyright by W illiam Harold W eedon, 1994
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
B R O A D B A N D MICROWAVE IN V E R SE SCATTERING :
TH EO RY A N D E X PE R IM E N T
W illiam H arold W eedon, P h .D .
D ep a rtm en t o f E lectrical an d C om p uter E n gineering
U n iv ersity o f Illinois a t U rban a-C h am p aign, 1994
W en g C ho C h ew , A dvisor
This thesis presents both theoretical formulations and experimented methods for
performing broadband time-domain inverse scattering. The inverse scattering prob­
lem is very important for a number of application areas including nondestructive
evaluation, geophysical probing, medical imaging and military target identification.
Emphasis is placed on the use of microwaves to probe the unknown object, although
much of the theory presented applies to other types of waves such as acoustic and
elastic waves.
Distorted-Born iterative method (DBIM) inverse scattering algorithms are pre­
sented for solving 2-D TM and TE problems. The TM algorithm is shown to be
capable of accurately inverting objects with contrast as high as 10:1, but the TE
algorithm breaks down when the object contrast exceeds 2:1 due to the buildup of
polarization charges inside the object.
A local-shape-function (LSF) inverse scattering algorithm is presented for imag­
ing very strong scattering objects such as metallic scatterers. The LSF algorithm has
a higher resolution capability than the DBIM algorithm for reconstructing closely
spaced metallic scatterers. The LSF algorithm also converges faster in the metallic
scatterer case.
Broadband time-domain data are preferable to continuous-wave (CW) data at
just a few discrete frequencies due to the higher information content inherent in a
broadband pulse and the ability to use time gating to eliminate unwanted earlytime and late-time arrival signals. Broadband time-domain data may be collected
in a practical microwave measurement system using either a step-frequency radar
approach or an impulse radar. There are advantages and disadvantages to using
iii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
both data collection methods, and the choice of one technology over the other
depends primarily on the application.
A prototype step-frequency radar system has been developed to demonstrate
the capability of our inverse scattering algorithms with real experimental data.
Reconstructions of both metallic and dielectric objects including metallic cylinders
and plastic PVC pipes in air from experimental data are shown. A commercial
monostatic impulse radar system is described and plans are discussed for building
a rudimentary bistatic impulse radar.
A method is presented for implementing an efficient finite-difference time-domain
(FDTD) electromagnetic scattering algorithm on a massively parallel supercom­
puter. The main challenge in designing an efficient algorithm is in the implemen­
tation of an absorbing boundary condition at the edge of the FDTD grid. Since
the inverse scattering methods that we present here rely on the solution of for­
ward scattering problems at each step in an iterative algorithm, the efficient FDTD
algorithm allows us to solve very large inverse scattering problems quickly on a
massively parallel supercomputer.
iv
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
D E D IC A T IO N
To my grandmother, Mildred L. Weedon, for her love
and support, and most of all for giving me a home.
v
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ACKNOW LEDGEM ENTS
First and foremost, I would like to thank my advisor, Professor Weng Cho
Chew, for the support, encouragement and guidance he has provided me over the
past three and a half years. I have learned many things during my stay at Illinois,
and it has been a great privilege to work with this man who I regard as one of the
brightest in my profession.
I would also like to thank the members of my doctoral committee, Professors
Shun-Lien Chuang, Paul E. Mayes, William D. O’Brien and Jose E. Schutt-Aine,
for their helpful discussions and suggestions.
A “h at’s off’ goes to my colleagues, past and present, in our research group, FuChiamg Chen, Yong-Hua Chen, Olivier Franza, Hong Gan, Levent Giirel, Jiun-Hwa
Lin, Caicheng Lu, Mahta Moghaddam, Muhammad Nasir, Gregory Otto, Jiming
Song, Robert Wagner, Yiming Wang and Xuguang Yang. I have shared an office
with many of these fine scholars at one time or another, and have benefited from the
knowledge that I have gained from working with each of them. I am particularly
indebted to my predecessors, Yiming Wang, Mahta Moghaddam, and Gregory Otto,
for providing the foundation on which much of the work in this thesis is based.
I would especially like to thank the undergraduate student Chad Ruwe for
spending two semesters of independent study with me and helping to build the
microwave imaging array presented in Chapter 4 of this thesis. I would also like
to thank the undergraduate student Bryan Hoon for the assistance that he has
provided near the end of this project.
For providing the prototype Vivaldi antenna design used in all of the microwave
measurements presented in this thesis, I thank Professor Paul E. Mayes and Dr.
David Tammen. I am indebted to Professor Jose Schutt-Aine for teaching me the
basics of microwave data collection and calibration techniques for helping me to get
started in the laboratory. I thank Dr. Frangois Colomb for showing me how to
collect data in the antenna test range and how to make good microwave connectors.
I would also like to thank Professor Carey M. Rappaport at Northeastern University,
vi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
my undergraduate Alma Mater, for many useful technical discussions on absorbing
material boundary conditions.
Thanks to Shirley Dipert, our Electromagnetics Laboratory secretary, for her
administrative assistance at various times during my stay at the Electromagnetics
Laboratory. Thanks also to Yan-Hua Pan for her assistance in preparing this thesis.
A thanks goes to the “gang” for providing me with relief on the weekends,
especially Jeff, Scott, Art, “G,” Joe, Jim, Brian, Dave, Jerry, Lisa, Marty and
John. Thanks also to the Lotito family for giving me a place of refuge on many
weekends and holidays.
Finally, I wish to thank my family for their support and love, especially my
parents and my grandmother to whom this thesis is dedicated. I am particularly
grateful to my father for the motivation that he provided me early on to excel in
academics and his financial support. I also thank my mother and stepfather for
their financial support and encouragement. In the words of Meng Chio (751-814),
“How could a blade of grass repay the warmth from the spring sun?”
This work was supported by the Office of Naval Research under grant N00014-89-J1286, the Army Research Office under contract DAAL03-91-G-0339, the
National Science Foundation under grant NSF-ECS-92-24466, and NASA under
grant NASA-NAG-2-871. Computer time was provided by the National Center for
Supercomputer Applications at the University of Illinois, Urbana-Champaign.
vii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
TABLE OF C O N TEN TS
CHAPTER
1
2
PA G E
IN T R O D U C T IO N
........................................................................................ 1
1.1 B ack grou n d on Inverse S catterin g T h e o r i e s ................................ 2
1.2
T im e-D o m a in Inverse S catterin g
..............................................5
1.3
1 .4
F D T D C o m p u ta tio n s on a M assively P arallel
S u p ercom p u ter
.......................................................................................6
T h esis O r g a n iz a t io n ...................................................................................7
1.5
R eferen ces
.................................................................................................. 9
D IS T O R T E D -B O R N IT E R A T IV E M E T H O D (D B IM )
.
.
14
2.1
C on ju gate G radient O p tim ization P r o c e d u r e ...............................15
2.2
2-D T M ( E z) P ola riza tion C a s e ................................................
2.3
2-D T E ( H z ) P o la riza tion C a s e .......................................................25
2 .4
B o rn Itera tiv e M eth o d (B IM )
2.5
2.6
V alid ity T est
...........................................................................................31
C o m p u ter S im u lation R e s u l t s ............................................................32
18
........................................... 30
2.6.1 TM polarization r e s u l t s ..................................................................... 32
2.6.2 TE polarization results
3
..................................................................... 38
2 .7
C o n c l u s i o n s ................................................................................................ 42
2.8
R eferen ces
................................................................................................43
LO C A L S H A P E F U N C T IO N (L SF) M E T H O D ............................... 45
3.1
A T -m a trix F o r m u l a t io n ...................................................................... 46
3.2
F rechet D erivative O perator
3.3
F rechet T ransposed O p e r a t o r ............................................................ 51
3 .4
Inverse S olu tio n U sin g C onjugate G r a d i e n t ...............................54
3.5
C om p u ter Sim u lation R e s u l t s ............................................................ 55
3.6
C o n c l u s i o n s ................................................................................................ 63
3 .7
R eferen ces
............................................................ 49
................................................................................................ 63
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4
B R O A D B A N D M ICR O W AVE IN V E R S E S C A T T E R IN G
M E A S U R E M E N T S .................................................................................... 65
4.1 S tep -F req u en cy R adar (S F R ) M icrow ave
M easu rem en t S y s t e m .................................................................68
4.1.1 Description of the SFR s y s t e m ........................................68
4.1.2 Broadband Vivaldi antenna e l e m e n t ..............................72
4.1.3 Early prototype antenna a r r a y ....................................
73
4.1.4 New switched antenna array
.................................................
4 .2 O verview o f M easu rem en t D a ta P r o c e s s i n g ....................82
4 .3 R eco n stru ctio n s U sin g E xp erim en tal S F R D a ta . . .
4.3.1 Metallic object reconstructions using early array
. . . .
4.3.2 Metallic object reconstructions using new a r r a y .......... 84
4.3.3 Dielectric object reconstructions using new array
. . . .
4 .4 C om m ercial Im p u lse R adar D a ta C ollection S y stem
. .
4 .5 R u d im en ta ry B ista tic Im p u lse R adar S y stem . . . .
4 .6 C o n c l u s i o n s ...................................................................................
97
4 .7 R eferen ces
.........................................................................................
5
6
E F F IC IE N T F D T D C O M P U T A T IO N O N A
M A S S IV E L Y PA R A LL EL S U P E R C O M P U T E R
. . . .
79
82
82
87
87
93
99
101
5.1
M od ified M axw ell E q u a t i o n s ................................................. 102
5.2
S in gle Interface P r o b l e m ............................................................104
5.3
P erfectly M atch ed Interface
5 .4
M od ified M axw ell E quations in th e T im e D om ain
5.5
C o m p u ter S im ulation R e s u l t s ..................................................109
5.6
C o n c l u s i o n s ......................................................................................113
5 .7
R eferen ces
..........................................................106
.
.
107
..............................................................................................114
C O N C L U S IO N S A N D D IR E C T IO N S F O R F U T U R E
W O R K ............................................................................................................ 115
V IT A
............................................................................................................ 120
ix
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
L IST OF SY M B O L S
E ( r , t)
Electric field vector (V/m).
H { r ,t)
Magnetic field vector (H/m).
9 ( r ,r ',t)
Scalar wave equation Green’s function.
G (r, r ', t )
Dyadic Green’s function.
e(r)
Inhomogeneous permittivity profile.
o-(r)
Inhomogeneous conductivity profile.
7 (r)
Inhomogeneous shape function profile.
M
Continuous model space.
M
Discretized model space.
V
Time-domain data space.
V
Frequency-domain data space.
F
Frechet derivative operator mapping elements of the continuous
model space M. to the frequency-domain data space V .
F
Frechet derivative operator mapping elements of the continuous
model space M to the time-domain data space V.
F
Frechet derivative operator mapping elements of the discretized
model space M to the frequency-domain data space V .
F
Frechet derivative operator mapping elements of the discretized
model space M to the time-domain data space V.
F1
Frechet transposed operator mapping elements of the frequencydomain data space T) to the continuous model space A i.
F1
Frechet transposed operator mapping elements of the time-domain
data space V to the continuous model space M .
F*
Frechet transposed operator mapping elements of the frequencydomain data space T> to the discretized model space M .
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
pt
Frechet transposed operator mapping elements of the time-domain
data space V to the discretized model space M .
K ( r ,r ', r n,t)
Kernel of Frechet derivative operator T .
K ( r ,r ', r„, t)
Kernel of Frechet derivative operator T .
F
Matrix representation of Frechet derivative
xi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 1
IN T R O D U C T IO N
In inverse scattering, one attempts to determine the internal profile of an inho­
mogeneous object by probing the unknown object with waves or fields. The inverse
scattering problem is then to generate a quantitative image of the unknown object
based on measurement data collected away from the scatterer.
Many different types of waves and fields are commonly used to image objects
including microwaves, acoustic waves, magnetic fields, and X-rays. These differ­
ent types of waves result in different imaging modalities and effectively add new
dimensions of artificial sight to our sensing capabilities. In addition to obtaining
the image of an object, a quantitative description of the scatterer such as its per­
mittivity, velocity, or conductivity profile is also obtainable from inverse scattering
methods and can contribute valuable diagnostic information. Different types of
waves carry different information about an object’s internal structure. In fact, the
different imaging technologies are complementary, and the more information the
investigator has about a structure or body, the better decisions that can be made
about its internal integrity or the location of an inhomogeneity.
Some example application areas of inverse scattering include nondestructive
evaluation, geophysical probing, medical imaging and military target identification.
In nondestructive evaluation, inverse scattering theory may be used to locate and
image possible cracks or defects in civil structures such as bridges, buildings, roads,
aircraft runways and t u n n els. Inverse scattering may also be used to generate images
of geophysical formations for locating minerals or buried objects such as hazardous
wastes. In medical imaging, X-rays as well as ultrasonic CAT scans and MRI are
commonly used to generate images of the human body. Military applications include
locating and identifying targets from radar data and finding unexploded ordnance.
The choice of a particular sensing technology for an application depends heav­
ily on how well the waves penetrate the background medium to be investigated.
1
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Acoustic waves axe good at penetrating metallic structures such as steel reactor
containment vessels, whereas microwaves do not propagate in solid metallic struc­
tures. Microwaves are good at penetrating many civil structures made of concrete
or brick, whereas it is difficult to couple acoustic waves into these materials. Each
technology has its advantages and disadvantages including resolution of the result­
ing image, cost, complexity of the measurement apparatus, portability and hazards.
In this thesis, we shall be concerned primarily with microwave inverse scatter­
ing. This technology is well-suited for many nondestructive evaluation applications,
particularly when investigating civil structures. Microwave inverse scattering is also
important for ground penetrating radar applications such as locating underground
tunnels, utilities, unexploded ordnance, and toxic wastes. The inverse scattering
theories developed here, however, are not limited to microwave inverse scattering
and may be applied in many other sensing applications.
1.1
B ack grou n d o n Inverse S catterin g T heories
The inverse scattering problem is often quite difficult, especially when wave
interactions are present. The reason is that, unlike visible light, many waves such
as microwaves and acoustic waves do not always travel in straight lines. This is due
to the phenomenon of diffraction, in which waves bend around objects to satisfy
boundary conditions and the underlying partial differential equations describing
the wave propagation. The inverse scattering problem is usually nonunique because
high spatial frequency portions of the object give rise to evanescent waves th at
cannot be measured. Hence, high spatial frequency information of the object is
often lost. Another phenomenon known as multiple scattering causes the scattered
field to be nonlinearly related to the scatterer. In addition to being nonunique and
nonlinear, inverse scattering problems are often ill-posed as well, due to the limited
measurement data as enforced by the problem geometry [1].
In developing new inverse scattering methods, we must keep in mind that the
goal of this research is to develop new methods and ideas th at will have a significant
impact on technology. Therefore, we propose th at new inverse scattering methods
2
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
should be judged against the following criteria:
(1) resolution of the reconstructed object;
(2) maximum object contrast that may be inverted accurately;
(3) size (in wavelengths) of the object that may be inverted;
(4) computational complexity (speed) of the algorithm;
(5) validity for objects with multidimensional spatial variation, arbi­
trary shape and arbitrary number of scatterers;
(6) robust, computationally stable implementation th at can withstand
the test of noisy or imperfect experimental data;
(7) requirement of a minimal amount of a priori information about
the object;
(8) validity for practical scattering measurement geometries where
measurement data are available only at discrete transmitter and
receiver locations that may not surround the object completely;
(9) sufficiently general method that may be applied to many wave
equations and scattering phenomena.
Unfortunately, there is presently no inverse scattering theory that meets all of the
above criteria perfectly. Below we summarize the work th at has been done to date
in the area of inverse scattering.
In general, the relationship between the scattered field and the scattering object
is a nonlinear one. This nonlinearity comes from multiple scattering effects within
the object [1]. However, a linear relationship can be found sometimes for certain
limited cases. Hence, inverse scattering theories can be categorized into linear
inverse scattering theories and nonlinear inverse scattering theories.
3
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
In the linear inverse scattering theories, approximations are made such th at a
linear relationship is found between the measured data and the scattering object.
The inverse scattering problem here amounts to determining the object function by
solving a set of linear equations. In X-ray tomography or computed tomography
(CT), the attenuation of an X-ray is linearly related to an integral summation of
the attenuations experienced by the X ray as it traverses a straight line path [2].
In nuclear magnetic resonance imaging (MRI), the received radio frequency (RF)
signal is proportional to an integral summation of the resonating spin densities along
a straight line [3-6]. The summation of an object function along a straight line
is also known as a projection. A backprojection algorithm using the projectionslice theorem can be used to reconstruct the object function efficiently from its
projections [1].
Another linear theory known as diffraction tomography takes into account diffrac­
tion effects when the scatterer is weak, or of low contrast [7-12]. A linear relation­
ship exists between the scattered field and the object because only single scattering
is important. Either the Bom or Rytov approximation [13] can be used to calcu­
late the scattered field and consequently a Fourier transform relationship between
the scattered field and the object function. The Fourier spectrum of the object
function may be recovered using a filtered backpropagation algorithm as in [7] and
the object therefore determined. Diffraction tomography is often used in ultrasonic
imaging [14] to achieve good reconstructions when multiple scattering effects are
small. However, when multiple scattering is important as is the case for objects
with strong contrasts, the technique breaks down [15].
When multiple scattering effects are important, different nonlinear theories have
been proposed to unravel the multiple scattering effects [16-25]. Many of these
n o n lin ea r theories have been proposed to solve the inverse scattering problem “ex­
actly” for 1-D objects, but have not been verified as practical methods for multidi­
mensional inverse scattering. Of particular interest is the work of Newton [19, 20]
who generalized the Gelfand-Levitan-Marchenko integral equation [1, 26] to multi­
dimensions. However, this method requires an impractical measurement apparatus
and has never been numerically implemented successfully. Yagle claims to have
4
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
a working 2-D layer-stripping algorithm [27], but this method applies only to the
Schrodinger equation. In [28], a ray tracing method is used for the forward problem
solver, but such a method is targeted for high-frequency applications.
The only general nonlinear inverse scattering theory valid for higher dimensions
that has been verified to date involves numerical methods [29-49]. These numerical
methods are usually iterative, and a forward scattering solver is needed in each
iteration.
1.2
T im e-D o m a in Inverse S catterin g
Broadband data, and in particular, time-domain data axe important for inverse
scattering for a number of reasons. The first is obviously th at the information
content available from a broadband transient pulse is usually much greater than
that available from CW measurements at a few discrete frequencies. W ith the
added information, the problem is usually better conditioned, or more “well-posed.”
A second reason is perhaps less widely known. One of the biggest problems with
conventional diffra c tio n tomography, especially in the area of medical ultrasound
imaging [12], is the “phase wrapping” problem. When only a few frequencies are
used to probe an object, and the object space occupies many wavelengths, phase
changes from one period to the next become wrapped and are hence ambiguous.
This problem does not exist when transient data are used along with a time-domain
inverse scattering algorithm such as the distorted-Bom iterative method (DBIM).
Another advantage to using time-domain data is that time gating may be used to
eliminate unwanted early-time and late-time arrival signals.
From a measurement perspective, it is usually more advantageous to collect data
in the frequency domain. One reason is that a transient or pulsed data collection
system is more subject to random electrical noise because of the high bandwidth.
Other problems with time-domain data collection include timing jitter due to noise
and uncertainties in the trigger and time base circuits as well as zero-level drift
and zero ambiguity [50, 51]. Narrowband microwave devices and electronics have
been in existence since the 1950s and are highly developed primarily due to the
5
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
radar industry [52]. Picosecond pulse measurement systems really did not reach
a mature development stage until the early 1970s [51-54]. Because of the wider
availability of narrowband microwave components, it is therefore often preferable
to use a frequency-domain measurement system.
To achieve the advantages of broadband time-domain measurement data using
a frequency-domain data collection system, a step-frequency radar (SFR) measure­
ment system was developed and used in this thesis. This system will be described
in Chapter 4 of this thesis, in which we will also describe and compare it with a
commercial impulse radar.
1.3
F D T D C o m p u ta tio n s o n a M assively P arallel S u percom puter
Although there are many advantages to using nonlinear iterative inverse scat­
tering algorithms, perhaps the biggest disadvantage is th at they require vast com­
putational resources. The reason is that the iterative algorithms rely on forward
scattering solvers to compute the electromagnetic fields incident to the scattering
object and compute the Frechet derivative and transposed operators. Many differ­
ent types of forward scattering solvers are available and may be used for the inverse
scattering problem. For example there are methods based on the finite element
method (FEM) [56, 57], the method of moments (MOM) [58, 59], and fast wave
scattering solvers such as the recursive aggregate T-matrix algorithm (RATMA) [59,
60]. We rely exclusively on the finite-difference time-domain (FDTD) technique [1,
61, 62] in this thesis since this method is very efficient for time-domain scattering
problems.
Full-wave forward scattering solvers such as those mentioned above are CPU
intensive due to the size of the problems that must be solved for many inverse
scattering problems. All of these methods require th at the scattering object re­
gion be discretized into cells on the order of 0.1 Amin where Amin is the minimum
wavelength of the excitation source in the scattering region. Hence, a moderate-size
scatterer that is 10 wavelengths in diameter would require that fields be computed
at (100)2 = 104 space locations for a 2-D scattering problem and (100)3 = 106 space
6
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
locations for a 3-D problem. This is one of the main reasons that the majority of
nonlinear inverse inverse scattering work to date has dealt with 2-D problems only.
Fortunately, the FDTD method is very amenable to implementation on a mas­
sively parallel supercomputer. The reason is that the computational molecule of the
FDTD algorithm involves only nearest-neighbor communications. The requirement
of only nearest-neighbor communications is very important for massively parallel
machines such as the Thinking Machines Corporation Connection Machine CM-5
supercomputer and allows the algorithm to be implemented very efficiently.
It turns out the most challenging aspect of programming an efficient FDTD code
on the Connection Machine is not in the core FDTD operations, but rather in the
absorbing boundary condition. Most conventional absorbing boundary conditions
[63, 64] rely on interpolation schemes for deriving the grid edge elements from
elements normal to or surrounding the grid boundary. This interpolation ru in s the
nearest-neighbor communications requirement of the algorithm and means that the
algorithm must use other, more costly types of communication to derive the grid
boundary values.
This problem with conventional absorbing boundary conditions suggests an al­
ternative approach whereby a synthetic absorbing material is placed at the edge of
the computation domain to absorb the waves. The advantage to using an absorbing
material boundary condition is that this approach does not require any additional
communications operations other than those of the core FDTD operations, and
may be implemented by modifying the material values at the nodes surrounding
the computation domain. A recent breakthrough in absorbing material boundary
conditions [65, 66] allows the specification of highly absorbing borders with border
widths of only a few grid points. This method has been verified in a 3-D FDTD
forward solver and implemented efficiently on the Connection Machine CM-5 [66].
1.4
T h esis O rganization
In this thesis, we present two nonlinear time-domain inverse scattering scatter­
ing algorithms that may be used to reconstruct objects directly from time-domain
7
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
data. These algorithms are the the distorted-Bom iterative method (DBIM) and the
local shape function (LSF) method. Both algorithms are considered time-domain
algorithms and use the full waveform of the probing wave to reconstruct an object
profile. Reconstructions using computer simulated data axe shown and demonstrate
th at the algorithms axe capable of achieving superresolution imaging.
The finite-difference time-domain (FDTD) algorithm is used as a workhorse in
all of the time-domain inverse scattering algorithms presented here. We will not
explicitly present the FDTD method since this basic numerical technique has been
presented many times in the past. Actually, a detailed knowledge of the FDTD
method is not necessary for understanding the development presented here. The
implementation of the algorithms discussed here would require such a knowledge,
however. The interested reader may refer to References [1, 61, 62] for excellent
presentations of the FDTD method.
Chapter 2 presents the theory of the DBIM for 2-D geometries for both trans­
verse magnetic (TM) and transverse electric (TE) polarizations. In Chapter 3, the
LSF algorithm is derived by first using a frequency-domain T-matrix method, and
then interpreting the results in the time domain to give a time-domain implemen­
tation of the LSF method. Inversions using computer generated scattering data
axe shown for both the DBIM and LSF algorithms. Chapter 4 discusses broadband
microwave scattering measurement techniques. A prototype step frequency radar
(SFR) im a g in g system developed at the University of Illinois is described. A com­
mercial impulse radar system for performing time-domain microwave measurements
is also described and compared against the prototype SFR system. Image recon­
structions are demonstrated using real experimental data for both dielectric and
metallic objects. Chapter 5 discusses efficient FDTD computations on a massively
parallel computer. Efficient forward scattering solutions are particularly impor­
tant for nonlinear inverse scattering because the iterative algorithms solve forward
scattering problems at each iteration. Finally, Chapter 6 concludes and suggests
directions for future work.
8
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1.5
References
[1] W. C . Chew, Waves and Fields in Inhomogeneous Media. New York: Van
Nostrand, 1990.
[2] A. C. Kak, “Computerized tomography with x-ray, emission and ultrasound
sources,” Proc. IEEE, vol. 67, no. 9, pp. 1245-1272, 1979.
[3] P. C. Lauterbur, “Image formation by induced local interactions: Examples
employing nuclear magnetic resonance,” Nature, vol. 242, pp. 190-191, 1973.
[4] P. C. Lauterbur and C. M. Lai, “Zeugmatography by reconstruction from pro­
jections,” IEEE Trans. Nucl. Sci., vol. NS-27, pp. 1227-1231,1980.
[5] W. V. House, “Introduction to the principles of NMR,” IEEE Trans. Nucl.
Sci., vol. NS-27, pp. 1220-1226,1980.
[6] I. L. Pykett, “NMR imaging in medicine,” Sci. Am., vol. 246, pp. 78-88, May,
1982.
[7] A. J. Devaney, “A filtered backpropagation algorithm for diffraction tomogra­
phy,” Ultrason. Imaging, vol. 4, pp. 336-360, 1982.
[8] A. J. Devaney, “A computer simulation study of diffraction tomography,” IEEE
Trans. Biomed. Eng., vol. BME-30, pp. 377-386, 1983.
[9] A. J. Devaney, “Geophysical diffraction tomography,” IEEE Trans. Geosci.
Remote Sensing, vol. GE-22, pp. 3-13, Jan. 1984.
[10] K. T. Ladas and A. J. Devaney, “Generalized ART algorithm for diffraction
tomography,” Inverse Probl., vol. 7, pp. 109-128,1991.
[11] K. T. Ladas and A. J. Devaney, “Application of an ART algorithm in an
experimental study of ultrasonic diffraction tomography,” Ultrason. Imaging,
vol. 15, pp. 48-58, 1993.
[12] T. J. Cavicchi and W. D. O’Brien, “Numerical study of high-order diffraction
tomography via the sine basis moment method,” Ultrason. Imaging, vol. 11,
pp. 42-74, 1989.
[13] M. Born and E. Wolf, Principles of Optics, 6th Ed. New York: Pergamon Press,
1980.
[14] M. Slaney, A. C. Kak, and L. E. Larsen, “Limitations of imaging with first-order
diffraction tomography,” IEEE Thins. Microwave Theory Tech., vol. MTT-32,
pp. 860-874, 1984.
[15] H. Lee and G. Wade, Modem Acoustic Imaging. New York: IEEE Press, 1986.
[16] F. C. Lin and M. A. Fiddy, “Image estimation from scattered field data,” Int.
J. Imaging Syst. Technoi, vol. 2, pp. 76-95, 1990.
[17] E. Baysal, D. D. Kosloff, and J. W. C. Sherwood, “Reverse time migration,”
Geophysics, vol. 48, no. 11, pp. 1514-1524,1983.
9
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
[18] J. G. Berryman and R. R. Greene, “Discrete inverse methods for elastic waves
in layered media,” Geophysics, vol. 45, no. 2, pp. 213-233, 1980.
[19] R. G. Newton, “Inverse scattering II: Three dimensions,” J. Math. Phys.,
vol. 21, no. 7, pp. 1698-1715,1980.
[20] R. G. Newton, “Inverse scattering HI: Three dimensions, continued,” J. Math.
Phys., vol. 22, no. 10, pp. 2191-2200, 1981.
[21] J. H. Rose, M. Cheney, and B. DeFacio, “The connection between timeand frequency-domain three-dimensional inverse scattering methods,” J. Math.
Phys., vol. 25, pp. 2995-3000, 1984.
[22] A. E. Yagle and B. C. Levy, “Layer-stripping solutions of multidimensional
inverse scattering problems,” J. Math. Phys., vol. 27, no. 6, pp. 1701-1701,
1986.
[23] K. P. Bube and R. Burridge, “The one-dimensional inverse problem of reflection
seismology,” SIAM Review, vol. 25, no. 4, 1983.
[24] V. H. Weston, “Invariant imbedding for the wave equation in three dimensions
and the applications to the direct and inverse problems,” Inverse Probl., vol. 6,
pp. 1075-1105, 1990.
[25] S. Coen, M. Cheney, and A. Weglein, “Velocity and density of a twodimensional acoustic medium from point source surface data,” J. Math. Phys.,
vol. 25, no. 6, pp. 1857-1861, 1984.
[26] T. M. Habashy and R. Mittra, “On some inverse methods in electromagnetics,”
J. Electromag. Waves A ppl, vol. 1, pp. 25-28,1987.
[27] A. E. Yagle and P. Raadhakrishnan, “Numerical performance of layer stripping
algorithms for two-dimensional inverse scattering problems,” Inverse Probl.,
vol. 8, pp. 645-665, 1992.
[28] J. G. Berryman, “Weighted least-square criteria for seismic traveltime tomog­
raphy,” IEEE Trans. Geosci. Remote Sensing, vol. GE-27, no. 3, pp. 302-309,
1989.
[29] Y.-M. Wang and W. C. Chew, “An iterative solution of two-dimensional elec­
tromagnetic inverse scattering problem,” Int. J. Imaging Syst. Technol., vol. 1,
pp. 100-108, 1989.
[30] W. C. Chew and Y.-M. Wang, “Reconstruction of two-dimensional permittiv­
ity using the distorted Bom iterative method,” IEEE Trans. Med. Imaging,
vol. MI-9, no. 2, pp. 218-225, 1990.
[31] M. Moghaddam, W. C. Chew, and M. Oristaglio, “Comparison of the Bom
iterative method and Tarantola’s method for an electromagnetic time-domain
inverse problem,” Int. J. Imaging Syst. Technol., vol. 3, pp. 318-333, 1991.
[32] M. Moghaddam and W. C. Chew, “Nonlinear two-dimensional velocity pro­
file inversion using time domain data,” IEEE Trans. Geosci. Remote Sensing,
vol. 30, Jan. 1992.
10
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[33] A. Tarantola, “The seismic reflection inverse problem,” in Inverse Problems of
Acoustic and Elastic Waves, F. Santosa, Y. H. Pao, W. Symes, and C. Holland,
Eds., Philadelphia: SIAM, 1984.
[34] E. Crase, A. Pica, M. Noble, J. McDonald, and A. Tarantola, “Robust nonlin­
ear waveform inversion: Application to real data,” Geophysics, vol. 55, no. 5,
pp. 527-538, 1984.
[35] A. Tarantola, Inverse Problem Theory. New York: Elsevier, 1987.
[36] R. E. Kleinman and P. M. van den Berg, “Nonlinearized approach to profile
inversion,” Int. J. Imaging Syst. Technol., vol. 2, pp. 119-126,1990.
[37] N. Joachimowicz, C. Pichot, and J.-P. Hugonin, “Inverse scattering: an itera­
tive numerical method for electromagnetic imaging,” IEEE Trans. Antennas
Propagat., vol. AP-39, no. 12, pp. 1742-1752, 1991.
[38] W. C. Chew and G. P. Otto, “Microwave imaging of multiple conducting cylin­
ders using local shape functions,” IEEE Microwave Guided Wave Lett.., vol. 2,
pp. 284-286, July 1992.
[39] G. P. O tto and W. C. Chew, “Microwave inverse scattering-local shape function
imaging for improved resolution of strong scatterers,” IEEE Trans. Microwave
Theory Tech., vol. 42, no. 1, pp. 137-141, 1994.
[40] G. P. Otto and W. C. Chew, “Inverse scattering of H . waves using local shape
function imaging: A T-matrix formulation,” Int. J. Imaging Syst. Technol.,
1994. Accepted for publication.
[41] S. A. Johnson and M. L. Tracy, “Inverse scattering solutions by a sine basis,
multiple source, moment method—part I: Theory,” Ultrason. Imaging, vol. 5,
pp. 361-375, 1983.
[42] S. A. Johnson and M. L. Tracy, “Inverse scattering solutions by a sine basis,
multiple source, moment method—part II: Numerical evaluations,” Ultrason.
Imaging, vol. 5, pp. 376-392, 1983.
[43] P. Mora, “Nonlinear two-dimensional elastic inversion of multioffest seismic
data,” Geophysics, vol. 52, no. 9, pp. 1211-1228, 1987.
[44] T. M. Habashy, W. C. Chew, and E. Y. Chow, “Simultaneous reconstruction of
permittivity and conductivity profiles in a radially inhomogeneous slab,” Radio
Sci., vol. 21, no. 4, pp. 635-645, 1986.
[45] T. M. Habashy, E. Y. Chow, and D. G. Dudley, “Profile inversion using the
renormalized source-type integral equation approach,” IEEE Trans. Antennas
Propagat, vol. AP-38, no. 5, pp. 668-682, 1990.
[46] A. Roger, “Newton-Kantorovitch algorithm applied to an electromagnetic in­
verse problem,” IEEE Trans. Antennas Propagat., vol. AP-29, no. 2, pp. 232238, 1981.
[47] D. T. Borup, S. A. Johnson, W. W. Kim, and M. J. Berggren, “Nonperturbative
diffraction tomography via Gauss-Newton iteration applied to the scattering
integral equation,” Ultrason. Imaging, vol. 14, pp. 69-85, 1992.
11
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[48] W. H. Weedon and W. C. Chew, “Time-domain inverse scattering using the
local shape function (LSF) method,” Inverse Probl., vol. 9, pp. 551-564,1993.
[49] W. C. Chew, G. P. Otto, W. H. Weedon, J. Lin, C. Lu, Y. Wang, and
M. Moghaddam, “Nonlinear diffraction tomography-the use of inverse scat­
tering for imaging,” Int. J. Imaging Syst. Technol., 1994. Submitted for pub­
lication.
[50] J. R. Andrews, “Comparison of ultra-fast rise sampling oscilloscopes.” Applic.
Note AN-2a, Picosecond Pulse Labs, 1989.
[51] G. H. Bryant, Principles of Microwave Measurements. London: Peter Peregrinus Ltd., 1988.
[52] M. I. Skolnik, Introduction to Radar Systems. New York: McGraw-Hill, 1980.
[53] J. R. Andrews, “Pulse measurements in the picosecond domain.” Applic. Note
AN-3a, Picosecond Pulse Labs, 1988.
[54] R. Lawton, S. Riad, and J. Andrews, “Pulse & time-domain measurements,”
Proc. IEEE, vol. 74, pp. 77-81, 1986.
[55] E. K. Miller, Ed., Time Domain Measurements in Electromagnetics. New York:
Van Nostrand Reinhold, 1986.
[56] J. Jin, The Finite Element Method in Electromagnetics. New York: Wiley
Interscience, 1993.
[57] H. Gan, “A model based inverse scattering method for nonlinear diffraction
tomography,” Ph.D. dissertation, Worcester Polytechnic Institute, 1993.
[58] R. F. Harrington, Field Computation by Moment Methods. Malabar, Florida:
Krieger, 1968.
[59] Y.-M. Wang, “Theory and applications of scattering and inverse scattering
problems,” Ph.D. dissertation, University of Illinois at Urbana-Champaign,
1991.
[60] W. C. Chew, L. Giirel, Y.-M. Wang, G. P. Otto, R. L. Wagner, and Q. Liu,
“A generalized recursive algorithm for wave-scattering solutions in two dimen­
sions,” IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 716-723, Apr.
1992.
[61] K. S. Yee, “Numerical solution of initial boundary value problems involving
Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat.,
vol. AP-14, pp. 302-307, 1966.
[62] A. Taflove, “Review of the formulation and applications of the finite-difference
time-domain method for numerical modeling of electromagnetic wave interac­
tions with arbitrary structures,” Wave Motion, vol. 10, pp. 547-582,1988.
[63] Z. P. Liao, H. L. Wong, B. P. Yang, and Y. F. Yuan, “A transmitting bound­
ary for transient wave analysis,” Scientia Sinica. (Series A), vol. 27, no. 10,
pp. 1063-1076, 1984.
12
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[64] G. Mur, “Absorbing boundary conditions for the finite-difference approxima­
tion of the time-domain electromagnetic field equations,” IEEE Trans. Elec­
tromag. Compat., vol. EMC-23, pp. 377-382, 1981.
[65] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic
waves,” J. Comput. Phys., Mar. 1994.
[66] W. C. Chew and W. H. Weedon, “A 3-D perfectly matched medium from modi­
fied Maxwell’s equations with stretched coordinates,” Microwave Opt. Technol.
Lett., 1994. Submitted for publication.
13
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CH A PTER 2
DISTORTED-BORN ITERATIVE METHOD (DBIM)
We present here the theoretical formulation for the time-domain distorted-Bom
iterative method (DBIM) for solving nonlinear 2-D inverse scattering problems. The
DBIM is an iterative optimization scheme whereby the distorted-Bom approxima­
tion is applied at each iteration. The distorted-Bom approximation linearizes the
nonlinear integral equation of scattering and is similar to the Bom approximation
with the exception that the distorted-Bom approximation assumes an inhomogeneous background medium, whereas the Bom approximation assumes that the back­
ground medium is homogeneous [1]. Another algorithm known as the Bom iterative
method (BIM) that employs the standard Bom approximation will be discussed as
a special case of the DBIM.
The most general method of solving nonlinear inverse scattering problems that
is valid for high contrasts, arbitrary number and shapes of scatterers and measure­
ment configurations is to use an optimization approach. Using an optimization
technique such as the conjugate gradient method, the Frechet derivative of the
scattered field measured a t the receivers is required for computing the gradient and
step size for the parameter update. Here, an operator formulation of the conjugate
gradient method is employed whereby the action of the Frechet derivative may be
computed as a single call to a forward scattering solver. We shall see th at the
Frechet transposed operator is required in addition to the Frechet derivative op­
erator. The Frechet transposed operator may be computed as a backpropagation
followed by a correlation.
We consider the DBIM algorithm in the context of both the transverse magnetic
(TM) or E z polarization case and the transverse electric (TE) or Hz case. The
TM case assumes a cylindrical z-directed source of electric current incident on a
cylindrical scatterer that is infinite in the z-direction. The TE case assumes a
cylindrical z-directed source of magnetic current incident on a cylindrical scatterer
14
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
infinite in the ^-direction. Although these problems are artificial in the sense that
it is impossible to produce current sources and scatterers that are infinitely long,
they do serve as a fair approximation for many problems. For example, the field
produced by a long wire or an antenna that produces a linearly polarized electric
field with a broad bandwidth in the E-plane would approximate the incident field
for the TM case fairly well. The TE case is also important because the variable
density acoustic wave equation is of the same form [2]. The advantage to solving a
2-D inverse is that the computational cost is much reduced as compared to solving
a 3-D inverse problem.
Both BIM and DBIM have been implemented previously for both CW and
transient excitations [3-7]. For the CW case, the recursive aggregate T-matrix
algorithm (RATMA) [8] has been used for the fast forward scattering solver. For
the transient case, a finite-difference time-domain (FDTD) algorithm is usually used
as the forward solver [9-14]. In this thesis, we shall concentrate exclusively on
time-domain inverse scattering using the FDTD algorithm as a forward scattering
solver.
2.1
C on ju gate G radient O p tim ization P roced u re
Before discussing the details of the DBIM inverse scattering algorithm, we shall
first briefly review an operator formulation of the conjugate gradient optimization
(CG) algorithm [15-17].
To use the CG algorithm, we first have to define a
functional or cost function to minimize. The parameter to optimize such as c(r),
We let the vector
( t ( t ) or e-1 (r) is represented by the vector * where [*]i =
0 meas represent a measured scalar field, or the exact value of the scalar field, at the
receivers. The vector <f>represents the computed field at the same location under
the current estimate for the model parameter x. Note that the vectors <f>meas and
<j>are indexed as a function of transmitter location, receiver location and time.
We define our cost function as
S(m) = \ ||0 - -T 'llg + | ||z> 15
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.1)
where the norms axe defined as
||4> - *■— 1|5 = {<j>-
■C ^ 1 ■(4> - < T '“ )
(2.2)
and
il®fc *fc—1||/^ = (*fc
*fc—1) ’ C m *(*fc
*fc—1) •
(2-3)
The parameters X k and X k - i represent parameter values at the current and previous
iteration steps. The matrices C D and C M are a priori data and model information
matrices. Usually, we set C M = e l and let C D be a function of time only
because we do not know the object location in the inverse problem. The second
term in Equation ( 2 .1 ) is known as the regularization term [ 1 8 , 1 9 ] and is added to
improve convergence of the algorithm. Hence, the information matrix C D is used
to implement a time-varying gain (TVG) in order to boost late time arrivals, while
e in
is used as a tuning parameter to improve convergence of the algorithm.
The iteration process is started by assuming some initial estimate for the pa­
rameter *o- Usually, this starting value is taken to be th at of the homogeneous
background medium. The parameter x is updated along the conjugate direction
hk,
x k = x k- i - OLk h k
( 2 .4 )
where a k is the step size for the update at the kth iteration. Initially, the conjugate
direction is set to the gradient g,
ho = go,
(2 -5 )
but is later updated as
hk —9h “I- 0k hk—x,
k > 1.
(2.6)
The step size for computing the conjugate direction may be computed as
A, =
( 2 .7 )
9 k - 1 *9 k -1
while the step size for the parameter update may be computed as
Qk =
p g)
h \ ’H k •h k
16
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
K
j
where H is the Hessian matrix.
From our definition of the cost function in Equations (2.1) - (2.3), the gradient
may be computed as
9 k = ^jh T
+ e ( * * - * * _ i)
»=»«.
(2-9)
where 8<f>= <f>— <j>meas. The matrix F is the matrix representation of the Frechet
derivative and is given explicitly as
' d<t>i
9X\
<f>2
F = 99x2
d<t>i
dX2
94>2
9x2
(2.10)
Similarly, the Hessian is given as
9fS (x)
_
‘
d x d x ‘ m=xk
= F •C
n
• F + s.
(2.11)
Substituting Equation (2.11) into Equation (2.8), we have
9k
h k
___________
(2.12)
(F • h ky • C DX• ( F . hk) + e h tk hk
Note in Equations (2.9) and (2.12) that the Frechet transposed F l is used
in computing the gradient, which defines the search direction, while the Frechet
derivative F is used only in computing the step size for the parameter update. In
practice, operator forms of the Frechet derivative and transposed operators are used
in computing the gradient and step size for update. The explicit operator forms of
the Frechet derivative and transposed operators will be presented in later sections.
However, we note th at using the Frechet derivative operator T , the gradient in
Equation (2.9) is computed as
gk = F { C ^ 1 ■5<f>} + e (xk - x k- i) .
(2.13)
Similarly, the vector F • h k in Equation (2.12) may be computed as the action of
the Frechet transposed operator J* on the vector h k.
17
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2.2
2-D TM (E x) Polarization Case
Consider the two-dimensional (2-D) scattering problem illustrated in Figure 2.1.
A line source of current J Zin(r,i) supported by the volume Q produces the elec­
tric field E Ztn ( r , t ) that is scattered by a 2-D cylindrical scatterer. We will assume
throughout this chapter that r e l 2. We shall use the subscript n to parameter­
ize the transmitter number because generally in an inverse scattering measurement
there will be multiple transmitter locations. The scatterer is characterized by the
permittivity and conductivity profile e(r) + <Se(r), <x(r) + 8cr(r) and exists in an
inhomogeneous background medium e(r), <r(r). That is, the scatterer consists of
a perturbation 8e(r), 8a(r) in the inhomogeneous background. In the formula­
tion that follows we shall assume that 8e(r) and 8a(r) are nonzero only within the
support volumes M e and M ° of the scatterers. Hence, the permittivity and con­
ductivity everywhere may be written as e(r) 4- 8e(r), cr(r) + 8cr(r). This is known
as the 2-D transverse-magnetic (TM) or ^--polarization scattering problem in an
inhomogeneous background medium.
Since both the line source and scatterer in our model have infinite extent in the
2-direction, the electric field will have only a ^-component. The vertical component
of electrical field E z,n(r,t) produced by the line source Jz>n(r,t) is given as the
solution to the scalar wave equation
V2 - w e M g j j - w r ( r ) ^ j £ „ „ ( r , t) = « , - • / „ „ ( r, t )
d2
d
+ fi0d e ( r ) ^ E z,n(r,t) + {iQ8a(r)— E z,n(r,t)
where we have assumed an isotropic, dispersionless medium. We define the inho­
mogeneous medium Green function g(r, r ', t) as the solution to
, ,d 2
_2
,
g (r,r ',t) = - 8 ( r - r') 8(t).
(2.15)
Then, we may employ the distorted Bom approximation to write the solution to
Equation (2.14) as
E z , n ( r , t )
» E z,n>o(r,t) + 8E ^n(r,t) + 8E°tTl(r,t).
(2.16)
18
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
z(r), a(r)
Cylindrical electric
current source
0.
Scatterer
E z, n ( r , i )
Electric field
+
e(r) 8e(r)
c(r) + da(r)
F ig u re 2.1 Two-dimensional TM scattering problem where the 2-D
scatterer Se(r), 5cr(r) consists of a perturbation of the background
inhomogeneous medium e(r), cr(r). The scatterer is excited by the
z-directed cylindrical source of electric current J z, n ( r , t ) .
19
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
In the above,
where Q represents the support of Jz,n{T',t) and E Zjnjo(r,t) is the incident field
in the presence of the background inhomogeneous medium e(r), cr(r). The terms
SEl n(r,t) and 5E^n (r,t) approximate the scattered fields induced by the permit­
tivity perturbation Se(r) and conductivity perturbation 8a(r). They are given as
and
where M € and M ° represent the model spaces where the parameters <5e(r') and
8<r(r') are defined. We assume in Equations (2.17) - (2.19) that the source Jz<n(r \ t')
is turned on at time t = 0 and generates a pulse that is for all practical purposes
zero for any time t > T at any space point r of interest. Equations (2.17) - (2.19)
are strictly valid in the limit T —¥ oo.
The integral equation given by (2.16) is only approximate because the distortedBom approximation [1, 9,10] has been used in writing Equations (2.18) and (2.19).
The approximation amounts to the fact th at the incident field E z,n$ inside integrals
in Equations (2.18) and (2.19) has been substituted in place of the total field E 2,n.
This approximation is equivalent to assuming that the scattered fields 8 E l n(r, t)
and 8Ez n (r,t) are weak compared to the incident field E Ztnto(r,t). The distorted
Bom approximation also linearizes the integral equation.
Comparing Equations (2.18) and (2.19) with Equation (2.17), we see that Equa­
tions (2.18) and (2.19) represent forward propagations of the induced sources
(2 .20 )
and
(2 .21 )
20
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
from the object spaces A i € and M ° to the receivers. Hence, they may be computed
losing a FDTD forward scattering solver. We shall expound on the significance of
these operators below.
Equations (2.18) and (2.19) can be thought of as operator forms of the Frechet
derivatives that map perturbations 5e(r) and 5a(r) into the field variations SE ^n(r, t )
and 6EZn (r, t). Wedefine the time-domain Frechet derivativeoperators
and J-a
as continuous linear operators mapping elements of the model spaces6e(r') € M e
and 5a(r') e M a to the data space T>:
: M 6 (->• T>
(2.22)
: M a ^ V.
(2.23)
and
The Frechet derivative operators
and
in Equations (2.18) and (2.19) may be
expressed as
*EJfB( r , t ) « 5 * { fc ( r ')}
r
~
(2.24)
= j dr1K e(r, r', r n, t) Se(rf)
Jm'
j
and
5 E Z J r ,t ) = r { 6 * ( r ' ) }
t
,
(2-25)
= I dr K a( r ,r , r n,t)8(r(r)
JM ”
where K r- (r,r' , r n,t) and K a{T ,r',rn,t) are the kernels of the Frechet derivative
integral operators. The kernels may be obtained from Equations (2.18) and (2.19)
as
K €( r , r ' , r n,t) = ~ /i0 Jq dt1g ( r , r ' , t - t') J ^ J 5 * >„f0(r/,t')
(2.26)
and
K a{r, r', r n, t) = -fj,0 jf dt' g(r, r', t - 1')
*')-
(2.27)
The Frechet transposed operators T*'1 and 7:a't corresponding to the Frechet
derivative operators map the field perturbations 5El n(r, t) and SE ^n(r,t) back
into the permittivity and conductivity spaces. That is,
:V ^ M €
21
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.28)
and
:
V
i—^ M
a.
( 2 .2 9 )
The Frechet transposed operators are defined to satisfy the inner products
t), F {& (r')} ) 5 = { Se(r%
{«%„(>•,«)} ) M,
(2.30)
{ < ^ , ( r , ( ) , r { J I,( r ') } ) 5 = < M r ' ) , r ' ' { ^ ( r , ( ) } ) Jl<.
(2.31)
(
and
for arbitrary 5E* n(r, t), 5e(r'), 6E ^n(r,t), Scr(r'). We define the various inner
products as
Nt Nr rT
{ f { r ,t ) ,9 { T ,t ) ) v = E E I d.tf{r,t)g{r,t),
( 2 .3 2 )
n=l m=l
{x{r,) ,y ( r ,) ) M<= f dr'x(r') y(r'),
JA/i<
( 2 .3 3 )
and
f d r ' x ( r ' ) y ( r ’).
(2.34)
JM.a
Since the Frechet transposed operators are linear operators, we may express them
in integral form as
«£( r ' ) = ^ - , {«£;,n(r>()}
Nt Nr
T _
= E E I & K €'t{ r \ r ~ t,Tr„ t ) 8 E tZtn(Tit)
( 2 .3 5 )
n=l m=l ®
and
fo (r ') = j * * {« ? ,„ (.•,* )}
Nt Nr
t _
= E E / # ^ ' V . » - m , r n,* )$££n(r,t)
( 2 .3 6 )
n=l m=l ■
/°
where K €,t{r' ,Tm, v n,t) and
, r m, r n,t) axe the kernels of the transposed
operators. We can easily show that the Frechet derivative and transposed operators
have the same kernel. That is,
K ‘( r , r ' , r n,t) = K ^ ( r f, r , r nyt)
( 2 .3 7 )
22
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
and
K a( r , r ', r n, t ) = K a,t(r' , r , r n,t).
(2.38)
Then substituting Equations (2.26) and (2.27) into Equations (2.35) and (2.36),
respectively, we have [9,10]
Nt
rT
Se(r) = -tto £
J
Q2
dt g j j £ 2,„,o(r',t)
" =1
(2-39>
,T
52
X
f
m
n=l
g(r',Tm,l! - t ) 6 E l fn(Tm ,1?)
0
and
Nt
8a{r)
=
-hq
fT
Y J
Q
dt —E z^ 0{r',t)
n=1 ° Nv r
X 52
(2.40)
-x
I
dt* g(r',Tm,t —t’) 5Egn{rmit >).
m = lJ°
Performing the change of variables r = T —t and r ' = T —1\ we then have
Nt ,T
Q2
5e(r) = - f i 0 Y JQ dT — .E z,„)0(r/, T - r )
N
r
*Y J
.T
dr' g(r', r m, r —r ;) <SE'jn(rm, T —r')
771=1
and
x
" =1
^
«,
Nr -X
x Jr_i
II %I/o
(2.42)
dT*g{r' ,Tm,T - t ’) 8EZn [rm, T - r ' ) .
771=1
Equations (2.41) and (2.42) are the desired representation of the Frechet trans­
posed operators. They may be implemented as a “backpropagation,” or timereversed propagation, followed by a correlation. In Equation (2.41), the induced
sources
f d tS F ^ r^ T -t),
JO
m = l,...,N a
(2.43)
23
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
axe backpropagated from the measurement space T> to the model space M.€. Al­
though the summation on receiver locations occurs outside the inner time integral
in Equation (2.41), linearity allows us to backpropagate all of the induced sources
'^z’ind(r rn5t ') , tti = 1 ,... ,N r simultaneously. This will cut down on computa­
tion time since it would otherwise require that the FDTD solver be nm separately
for each receiver location. The second step in computing the Frechet transposed
operator requires that the backpropagated field be correlated with the second time
derivative of the incident field inside the object, or multiplied by QT*Ez,nfl{r 'i T t )
and integrated over time. Finally, the process is repeated for each transmitter lo­
cation and the result summed.
Similarly, Equation (2.42) requires a backpropagation of the induced sources
<fnd(»*m ,r')= f Q d t8 E ^ n(rm, T —t),
m = l,...,iV *
(2.44)
from the measurement space f> to the model space M ? . This backpropagated field
is then correlated with the first time derivative of the object field. The process is
again repeated and summed over the N t transmitter locations.
The distorted Bom approximation is used frequently in diffraction tomography
to perform inverse scattering on objects with weak contrast compared to a known
background. But instead of applying the distorted Bom approximation only once,
this approximation may be applied repeatedly if the background medium is updated
at each step. When the distorted Bom approximation is used in an iterative fashion,
the resulting algorithm is known as the distorted Bom iterative method (DBIM).
The solution th at is obtained from the DBIM solves the nonlinear inverse problem
and hence is valid for much larger contrasts than if the distorted Bom approximation
were to be applied only once.
Both the Frechet derivative and transposed operators are required in an iterative
optimization scheme such as the conjugate gradient method. The Frechet derivative
operator is used in computing the conjugate gradient step size for update along a
given search direction. The Frechet transposed operator is used to compute the
gradient and hence the search direction. In the DBIM, we may replace e(r>) and
24
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
<r(r) with €k(r ) and crfc(r), the permittivity and conductivity at the Aith iteration.
The field quantity Ez,n,k is the incident field at the fcth iteration in the presence of
the background medium efc(r), <Tk(r). The object model parameters efc(r), crk(r)
are then updated at each iteration.
2 .3
2-D T E ( H z) P o larization C ase
For the TE polarization case, we assume the same geometry as the TM case
with the exception th at the cylindrical source of electric current JZyTl is replaced
by a cylindrical source of magnetic current M Zj„ supported by volume Q. The
electromagnetic field produced by the source M 2j„ is incident on a 2-D dielectric
scatterer e-1 (r) + Je-1 (r) embedded in the inhomogeneous background medium
e~x(r) as shown in Figure 2.2. We assume that the perturbation <5e-1 (r) is nonzero
only over the region of support M e. We shall assume here that the medium is
lossless, or a(r) = 0 since we cannot write down directly a second-order wave
equation for H z alone when <r(r) ^ 0 without neglecting the displacement term.
The wave equation for Hz (r,t) assuming cr(r) = 0, n(r) = (iq and an isotropic,
dispersionless medium may be written as
(2.45)
We define the inhomogeneous Green function g(r, r ', t) to satisfy
(2.46)
Using the principle of linear superposition and our definition of the Green function
above, we may write the solution to Equation (2.45) as
H s,n(»*51) « H s,nfi{r, t) + 5Hz,n{r, t)
where
and
25
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Cylindrical magnetic
current source
e-J(r)
z,n(r >0
Scatterer
t
Hz>n(r,
Magnetic field
e~J(r) + 8 e~J(r)
F ig u re 2.2 Two-dimensional TE scattering problem where the 2-D
scatterer <5e- 1 (r) consists of a perturbation of the background inho­
mogeneous medium e-1 (r). The scatterer is excited by the z-directed
cylindrical source of magnetic current M s,n{r,t).
26
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The approximation in Equation (2.46) above is due to the fact we have used the
distorted-Bom approximation to replace the total field H z,n (r,t) with the incident
field HZinto(r,t) in Equation (2.48) above.
Comparing Equation (2.48) with (2.47), it is clear that Equation (2.48) is a
forward propagtion of the induced source
M z,ind(r,, 0 = - /* dt V' • ( f c r V ) V'iyS)n,o(r',i)
J0
(2.49)
from the model space M € to the data space V. However, from Ampere’s law, the
gradient of H z is given as
V H z (r,t) = - x e ( r ) ^ E y{r,t) + y e ( r ) ^ E x {r,t).
(2.50)
The components E x and E y are available if the 2-D TE Yee FDTD algorithm is
used as a forward scattering solver, rather than a FDTD algorithm that solves the
scalar wave equation for H~ only. Substituting Equation (2.50) into (2.48), we have
(2.51)
Hence, the induced source may be written in the equivalent form
= £ !S t - 1(T')e(T')Ey(T',t’) - ^ 5 £- 1(r')£ (r')E I (r',t').
(2.52)
The derivatives ^ and ^ in Equations (2.52) may be computed numerically using
standard forward differencing.
Equation (2.48) is an expression for the Frechet derivative
T :fc -V )
6Hz,n(r,t).
(2.53)
This expression is convenient because the it allows the Frechet derivative to be com­
puted using a forward scattering solver with the induced source V' • <fe-1 V H z^n,o*
27
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
However, Equation (2.48) is not convenient for deriving the Frechet transposed op­
erator because the perturbation 5c"1( r 1) is embedded inside the V'- operator and
cannot be separated from V H ZfTl$ ( r ' ,t' ) in the present form.
To circumvent this problem, we use integration by parts along with Gauss’
divergence theorem to rewrite Equation (2.48) as
= f T dt' f
0 ,T
- 1
d r ' [VV(r,r',t-(')]-*-1(>-')V'H'2l„,o(r',i')
;
* ' JeM ■ d r ' s ( r y , t - o & - V )
8
<2-54>
The second term in Equation (2.54) above is zero because the integration volume
A4e and surface d /A € may be deformed slightly such that the support of <Je-1 (r) is
included wholly inside of the volume M . €. The value of 8e ~ 1( r /) is then zero on the
surface d M e. Equation (2.54) becomes
5H z n ( r , t ) = Jf0 d t' J.M*
f d r , \ ^ ,g ( r y i t - t ,) ) - S e - l { r ,) V H Ztn,0( r , i tf).
(2.55)
The Frechet derivative implied by Equation (2.55) may be written as
S H Ztn( r , t ) = ^‘{5€- 1(r,)|
= jr d r ' K ( r , r ' , r n , t ) 8e i (r ')
(2-56)
K ( r , r ' , r n , t ) = f d t ' V 9( t , t ' , t - t ' ) - V H ^ r ' , t ' ) .
(2.57)
where
jo
The Frechet transposed operator J* is defined to satisfy the inner product
relation
{ f c - V ') } ) 5 = ( f c - ' M , ? » , ( - ■ , ( ) } ) „ ,
(2.58)
for arbitrary 8e ~ l ( r r) and 8H z ,n ( r , t ) . Since the Frechet transposed operator is a
linear operator, we may express it in the form
«5e-1(r/) = ^ t {5H2(n(r,t)}
Nt
N
r
.t
= £ J L I d . t K t ( r ' , r m , r n , t ) 8H z ,n ( r m , t )
n=l m =l
Nt
N
(2.59)
r
= £ £ / d t K ( r m , r \ r n , t ) 8H Zjn( r m , t )
n = l m = l J°
28
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where we have used the fact that
K ( r , r ' , r n,t) = K ^ r 1, r , r n,f)
(2.60)
Substituting Equation (2.57) into Equation (2.59), we have
Nt rT
S t - \ r ' ) = £ J d f V J W r ' ,* ' )
"=l ° Nr fT
Y l f dt V 'g(rm, r', t - tf) 8Hz>n(rm, t).
(2.61)
i J0
Performing the change of variables r = T — t and t ‘ = T —t ’, using reciprocity to
replace g ( r , r ',t —t') with g ( r ', r ,t —t’) and moving the V' operator outside the
inner time integral, we have
Nt fT
&-V) =S /
"-1 °
Nr fT
• v ' 53 f d r g i r ' t V r n y - T ) 8 H 2tn(rm, T - r ) .
(2.62)
m =l
Equation (2.62) is the desired expression for the Frechet transposed operator.
For numerical computations, however, we use Equation (2.50) to replace V if, with
Ex and E y field components. First, let us examine the gradient of the backpropa­
gated field
V ’H Z,bp = V'
Nr rT
£ / d T g (r \T m, r ' - r) 8Hz,n(rm, T - r).
m=l J°
(2.63)
= x e ( r ' ) ^ E ytbp(r', t ’) + y e{r’) - jp E x,bp{r', /)■
(2.64)
From Equation (2.50),
V 'tf2)bp(r', t ')
The field components E x^ p(r',T') and E J,)bp((^’^ r ,) would be produced if the in­
duced source
M2>ind(rm, r) = f dt 8Hz,n(rm, T - t )
(2.65)
Jo
was propagated from the data space V to the model space M €. If the induced
source
^ M 2)ind(rm, r) = 8HZtn(rm, T - t )
(2.66)
29
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
was used instead, the fields ^ p E x^p{T', t >) and -^pEy^ v{r',r') would be produced
in the object space. Hence, the x- and y-components of V'fi^bpC^', r') may be com­
puted by propagating the induced source
,ind(**mj t ) given by Equation (2.66)
from the data space V to the model space M €, taking the E x and E y components,
respectively, of the backpropagated field, and multiplying by the permittivity e(r').
As in the TM polarization case, the field data at all of the receiver locations may
be backpropagated simultaneously to reduce the computational cost.
The internal incident field
o(r ', T') may also be computed from the E x
and Ey components with the aid of Equation (2.50). The inner product of
V'iJz.bpC**', t ') and V fl’j )„lo (r',T —r*) is computed and integrated over time, which
in effect correlates the various components of the backpropagated field with the
internal incident field. A final summation is taken over the transmitter locations.
2 .4
B orn Itera tiv e M eth o d (B£M )
Another method of solving nonlinear inverse scattering problems in the time
domain is to use the Bom iterative method (BIM) [10, 12]. The major difference
between the BIM and the DBIM is th at in the BIM, the background Green function
is not updated at each iteration. Rather, the background Green function is chosen
to be that of a homogeneous medium, which may be computed in closed form.
For both BIM and DBIM, however, the field internal to the object is computed
numerically and updated at each iteration.
A BIM algorithm for transient data was presented previously [10] by converting
the time-domain field data to the frequency domain and solving the inverse problem
in the frequency domain using a known analytical expression for the Green function.
We will not present this method here. Rather, we consider the BIM algorithm to
be a special case of DBIM when the inhomogeneous Green function gk( r ,r ',t) is
replaced by a homogeneous medium Green function.
The main advantage to using the BIM algorithm is th at there is evidence that
the BIM performs better than the DBIM when the scattering data are noisy [12].
However, the DBIM can be shown to have second-order convergence and is hence
30
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
valid for larger contrasts than the BIM.
2.5
V alid ity T est
In a computer implementation of the algorithms presented thus far, it would be
very easy to make a programming error in computing the Frechet derivative and
transposed operators and miss, for example, a factor of A t or A x. The inverse
scattering algorithm would then either not converge or give an erroneous result.
This section presents an easy way to check for this type of programming error. In
the computer programs used in this thesis, the Frechet derivative and transposed
operators were computed as separate subroutine calls. This modular programming
approach allows for easy error checking.
First, the Frechet derivative operators given by Equations (2.2.18), (2.2.19)
and (2.3.51) axe tested. We note that the Frechet derivative operators are valid
in the limit 8e —> 0 and 8a
0 due to the distorted-Bom approximation. The
Frechet derivative operators give the perturbation in the scattered field measured
at the receivers due to a perturbation in either Se or 8a. Therefore, an easy way to
test if the Frechet derivative operator computation is correct is to specify a small
perturbation <Se(r) or 8a(r) and use an FDTD forward solver to actually compute
the scattered field due to the perturbation. The result of Frechet derivative operator
operating on the perturbation is then compared against the forward solver result.
To test the Frechet transposed operators, we note that the Frechet trans­
posed operators axe defined by the inner products (2.2.30), (2.2.31) and (2.3.58).
Moreover, the inner products are valid for arbitrary 8Ez n (r,t), 8e(r') in Equa­
tion (2.2.30), 8Ez n {r,t), 8a(r') in Equation (2.2.31) and arbitrary 8e~1(r/) and
8Hx<n(r, t) in Equation (2.3.58). Thus, it is allowable to let 5Et n {r, t) = ^ { ^ ( r ') } ,
8 E °n(r,t) = Jr{8a{r')} and 8Hz,n{v,t) = T{8€~1(r')}. The Frechet transposed
operator test then reduces to specifying arbitrary values of <5e(r'), 8a{r') and
<Je- 1 (r') and testing the equality of the three inner products
L
<2-67)
( T * {*<’ (’• ' ) } ' ? * {fo(r')} ) . I ( S ^ ( T ' ) , ^ ‘ { ^ { S ^ r ' ) } } ) M, ,
31
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
(2.68)
and
(* -{ fc -V )}
2.6
)- 2 ( f c - V ) .* 4 {? {fc -‘M }} ) M. •
(2-69)
C o m p u te r S im u latio n R e su lts
2.6.1 TM polarization results
For all of the computer simulation results presented in this chapter, the fullangle scattering geometry shown in Figure 2.3 was used. A 2-D Yee FDTD forward
scattering solver [20, 21] was used with a iV* x Nj FDTD grid. Internal to the FDTD
grid was a subgrid region of dimension Oi x Oj grid points where the object function
was located and solved for in the inverse scattering problem. A total of 4 transmitter
locations and 8 receiver locations were placed around the object region at a distance
dtr from the center (see Figure 2.3). At the edge of the FDTD grid, a perfectly
matched absorbing material boundary condition [22, 23] occupying a border of
width 8 grid points was used to absorb outgoing waves with minimal reflection
from the FDTD grid boundary. Both the forward and inverse scattering problems
were solved on a Connection Machine CM-5 massively parallel supercomputer.
For all of the computer simulations presented in this chapter, the grid parame­
ters used were A x = A y = 2.5 m m and A t = 5.5 ps. The electric current density
source pulse used for the TM cases was of the form
J ‘ = A x A y [4 (t/T )3 " (t/T)“] e~’/r (A/m2)'
(2' 70)
The magnetic current density source pulse used for the TE cases was also of the
form
=
t4(i/T)3 “ (</T)4]
(V/m2)’
(2’71)
In Equations (2.70) and (2.71), r = l/47r/o, where the center frequency of the pulse
was chosen to be fo = 2.0. The background medium was assumed to be that of free
space in all cases.
Figure 2.4 shows a DBIM reconstruction of a crescent-shaped (in cross section)
smooth object with m a x im u m permittivity €r<max = 2 .0 using a TM-polarized inci­
dent wave. The object had a radius of 6 grid points and the reconstructions were
32
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
dtr
dtr
-------A
Absorbing
border
— E3-------------B ------------ 0 - ~
__s
s
0 __
Optimization
subgrid
I
I
□
: Receiver location
H
: Transmitter/receiver location
F ig u re 2.3 Full angle scattering geometry used for computer simula­
tion results. FDTD grid contains IV* x Nj grid points which contain
an optimization subgrid that contains 0* x O j grid points, transmit­
ter and receivers located around the optimization subgrid, and an
absorbing border region.
33
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
computed using a 64 x 64 FDTD grid with a 21 x 21 subgrid. The transmitters and
receivers were located at a distance of dtr = 16 grid points (see Figure 2.3). The
original object is shown in (a), while the reconstructions after 10 iterations and 20
iterations, are shown in (b) and (c), respectively. The relative residual error (RRE),
defined as
n p p _
~ tin e a s )* *
Axt
.
~ <ftmeas)
n ~ X. A\____
(2.72)
is shown for the two reconstruction cases. After only 10 iterations, the maximum
contrast of the object is reconstructed, but the object appears “washed out,” and
many of the fine details are not resolved. After 20 iterations, the resolution is much
better. The object could be resolved to an even greater detail with more iterations,
but the program was stopped after 20 iterations.
In Figure 2.5, the reconstruction of an object similar to the one in Figure 2.4, but
with maximum contrast erimax = 10.0 is shown. All of the other object parameters
are the same as those used in Figure 2.4. Note th at it takes many more iterations
for the algorithm to converge for this high-contrast case, since the distorted-Bom
approximation th at is applied at each iteration is stretched to a region of question­
able validity. After 60 iterations, the shape and maximum contrast of the object are
reconstructed, but the object reconstruction appears to be noisy. Apparently, this
case is near the maximum contrast th at the DBIM algorithm is capable of resolving.
Figure 2.6 shows a reconstruction of a much larger object than the ones in
Figures 2.4 and 2.5. The object is the same shape as the ones in the two previous
cases, but here the object radius is 30 grid points, and is computed using a 128 x 128
FDTD grid with an 81 x 81 optimization subgrid. The maximum permittivity for
this case is er)max = 2.0. After only 20 iterations, the object is reconstructed well.
This case is important because it demonstrates th at large objects that are several
wavelengths in dimension can be reconstructed. The reconstruction of large objects
has been a problem in the past for diffraction tomography algorithms, particularly
in ultrasound tomography [24], due to a phenomenon known as “phase wrapping.”
For objects with contrasts around er — 2.0, the TM DBIM algorithm performs
34
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
(b) 10 Iter., RMS Err.= 0.113
Relative Permittivity
(a) Original Object
1.5-
Y
0 0
X
Relative Permittivity
(c) 20 Iter., RMS Err.= 0.0574
Y
0 0
F ig u re 2.4 Original object and DBIM TM-Polarization permittivity
optimization reconstructions of crescent-shape smooth object with
er,max = 2.0. The object had a radius of 6 grid points and was com­
puted using a 64 x 64 FDTD grid with a 21 x 21 optimization subgrid.
Full-angle scattering geometry shown in Figure 2.3 containing 4 trans­
mitters and 8 receivers is used with dtr = 16 grid points.
35
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(b) 10 Iter., RMS Err.= 0.761
Relative Permittivity
(a) Original Object
20
0 0
(d) 30 Iter., RMS Err.= 0.266
Relative Permittivity
(c) 20 Iter., RMS Err.= 0.511
Relative Permittivity
(e) 40 Iter., RMS Err.= 0.200
(f) 60 Iter., RMS Err.= 0.165
lOi
10 -
Y
0 0
X
F ig u re 2.5 Original object and DBIM TM-Polarization permittivity
optimization reconstructions of crescent-shaped smooth object with
£r,max = 10-0- The object had a radius of 6 grid points and was com­
puted using a 64 x 64 FDTD grid with a 21 x 21 optimization subgrid.
Full-angle scattering geometry shown in Figure 2.3 containing 4 trans­
mitters and 8 receivers is used with dtr = 16 grid points.
36
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(b) 10 Iter., RMS Err.= 0.258
Relative Permittivity
(a) Original Object
Relative Permittivity
(c) 20 Iter., RMS Err.= 0.0434
F ig u re 2.6 Original object and DBIM TM-Polarization permittivity
optimization reconstructions of crescent-shaped smooth object with
er ,m ax = 2.0. The object had a radius of 30 grid points and was
computed using a 128 x 128 FDTD grid with a 81 x 81 optimization
subgrid. Full-angle scattering geometry shown in Figure 2.3 contain­
ing 4 transmitters and 8 receivers is used with dtr = 50 grid points.
37
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
extremely well. Large objects of several wavelengths in diameter do not pose a
problem for reconstruction, with the exception that they take more computer time.
For high contrast objects, the TM DBIM algorithm can resolve objects with per­
mittivities of er = 10.0, but more iterations are required since the validity of the
distorted-Bom approximation applied at each iteration is stretched.
2.6.2 TE polarization results
The same cases th at were considered above for TM-polarization are repeated
here for TE-polarization. In Figure 2.7, the crescent-shaped smoothed object with
radius 6 grid points and maximum permittivity er,max = 2.0 was reconstructed.
After 20 iterations, the object is reconstructed fairly well with the exception that
the top of the object appears noisy or “spiky.” After 30 iterations, the top of the
object is still not reconstructed properly.
In Figure 2.8, the object contrast is increased to er,max = 10.0. After 10 itera­
tions, the object is poorly resolved and has a single spike up to er = 10.0. After 20
iterations, the object is still poorly resolved, and now the spike extends nearly up
to a permittivity value of er = 20.0.
If Figure
a large object with radius 3 0 grid points and maximum permittiv­
ity €r ,m ax = 2 . 0 is reconstructed using T E polarization. The object is well resolved
after 20 iterations, demonstrating that large, low-contrast objects may be recon­
structed well using the T E DBIM algorithm.
2 .9 ,
For objects with maximum contrasts around er = 2.0, the TE DBIM algorithm
performs well even for large objects of several wavelengths in diameter. Objects
with contrasts greater than er = 2.0, are not reconstructed properly using this
implementation of the 2-D TE DBIM algorithm. The reason is th a t the TE problem
is really a vector inverse scattering problem in disguise, due to the polarization
charges th at are built up inside the object. This is evident from the fact th at the
gradient of the H z field, or equivalently the E x and E y fields inside the object were
required in computing the Frechet derivative and transposed operators.
38
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(b) 10 Iter., RMS Err.= 0.135
Relative Permittivity
(a) Original Object
£ 1.5-
Y
0 o
Y
X
Relative Permittivity
(c) 20 Iter., RMS Err.= 0.0649
0 0
X
(d) 30 Iter., RMS Err - 0.0424
>
1
Y
0 0
Y
X
0 0
X
F ig u re 2.7 Original object and DBIM TE-Polarization permittivity
optimization reconstructions of crescent-shaped smooth object with
€r,max = 2.0. The object had a radius of 6 grid points and was com­
puted using a 64 x 64 FDTD grid with a 21 x 21 optimization subgrid.
Full-angle scattering geometry shown in Figure 2.3 containing 4 trans­
mitters and 8 receivers is used with dtr = 16 grid points.
39
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
5 Iter., RMS Err= 0.739
Relative Permittivity
Original Object
10
,
» 1-5
Y
0 0
Y
20 Iter., RMS Err.= 0.503
10 Iter., RMS Err= 0.550
Relative Permittivity
0 0
10-
•3 15-
10
20
20
Y
0 0
Y
0 0
F ig u re 2.8 Original object and DBIM TE-Polarization permittivity
optimization reconstructions of crescent-shaped smooth object with
er,max = 10-0. The object had a radius of 6 grid points and was com­
puted using a 64 x 64 FDTD grid with a 21 x 21 optimization subgrid.
Full-angle scattering geometry shown in Figure 2.3 containing 4 trans­
mitters and 8 receivers is used with dtr = 16 grid points.
40
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(b) 10 Iter., RMS Err.= 0.258
Relative Permittivity
(a) Original Object
Relative Permittivity
(c) 20 Iter., RMS Err.= 0.0434
F ig u re 2.9 O rigin al object and DBIM TE-Polarization permittivity
optimization reconstructions of crescent-shaped smooth object with
€r, m ax = 2 .0 . The object had a radius of 30 grid points and was
computed using a 128 x 128 FDTD grid with a 81 x 81 optimization
subgrid. Full-angle scattering geometry shown in Figure 2.3 contain­
ing 4 transmitters and 8 receivers is used with dtr = 50 grid points.
41
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2.7
Conclusions
This chapter presented the distorted-Bom iterative (DBIM) a lg o rith m s for 2 D TM and TE polarization electromagnetic inverse scattering. The iterative al­
gorithms were developed from an operator formulation of the conjugate gradient
optimization algorithm. The distorted-Bom approximation was used to linearize
the integral equation of scattering whereby the total field internal to the object was
replaced with the incident field at the kth. iteration. The linearized integral equa­
tion of scattering allowed the computation of the Frechet derivative as a single call
to a forward scattering solver. The Frechet transposed operator was derived from
the Frechet derivative operator and could also be computed as a backpropagation
followed by a correlation. This allows the IV-dimensional parameter search direc­
tion and update step size to be computed with only 3 calls to a forward scattering
solver.
Computer simulation results were shown to demonstrate the effectiveness of the
DBIM algorithms. The TM DBIM algorithm was shown to be capable of inverting
objects with contrasts a large as er = 10.0. The larger contrast object example took
longer to converge than the lower contrast object, as expected. Large objects of
several wavelengths in diameter pose no reconstruction problems.
The TE DBIM algorithm was seen to be capable of inverting objects with con­
trast as large as er = 2 .0 , but did not give a satisfactory reconstruction when the
contrast was increased to er = 10.0. The reason is that the 2-D TE problem is really
a vector inverse scattering problem due to the buildup of polarization charges inside
the object. Large objects of several wavelengths in diameter pose no reconstruction
problems for the TE polarization, as long as the contrast is below er = 2 .0 .
The results presented here are limited since the DBIM and the BIM algorithms
have been studied quite extensively already [3-7,10-13]. Other cases of interest and
importance such as the limited angle inverse problem, the maximum resolution of
the algorithms, and conductivity optimization have been considered in the references
listed above.
42
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2.8
R eferences
[1] M. Bom and E. Wolf, Principles of Optics, 6th Ed. New York: Pergamon Press,
1980.
[2] W. C. Chew, Waves and Fields in Inhomogeneous Media. New York: Van
Nostrand, 1990.
[3] Y.-M. Wang and W. C. Chew, “An iterative solution of two-dimensional elec­
tromagnetic inverse scattering problem,” Int. J. Imaging Syst. Technol., vol. 1,
pp. 100-108, 1989.
[4] W. C. Chew and Y.-M. Wang, “Reconstruction of two-dimensional permittiv­
ity using the distorted Bom iterative method,” IEEE Trans. Med. Imaging,
vol. MI-9, no. 2, pp. 218-225, 1990.
[5] Y.-M. Wang and W. C. Chew, “Limited angle inverse scattering problems and
their applications for geophysical explorations,” Int. J. Imaging Syst. Technol.,
vol. 2, no. 2, pp. 96-111, 1990.
[6] Y.-M. Wang and W. C. Chew, “Accelerating the iterative inverse scattering
algorithms by using the fast recursive aggregate T-matrix algorithm,” Radio
Sci., vol. 27, no. 2 , pp. 109-116, 1992.
[7] Y.-M. Wang, “Theory and applications of scattering and inverse scattering
problems,” Ph.D. dissertation, University of Illinois at Urbana-Champaign,
1991.
[8] W. C. Chew, L. Giirel, Y.-M. Wang, G. P. Otto, R. L. Wagner, and Q. Liu,
“A generalized recursive algorithm for wave-scattering solutions in two dimen­
sions,” IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 716-723, Apr.
1992.
[9] W. H. Weedon and W. C. Chew, “Time-domain inverse scattering using the
local shape function (LSF) method,” Inverse Probl., vol. 9, pp. 551-564,1993.
[10] M. Moghaddam, W. C. Chew, and M. Oristaglio, “Comparison of the Bom
iterative method and Tarantola’s method for an electromagnetic time-domain
inverse problem,” Int. J. Imaging Syst. Technol., vol. 3, pp. 318-333,1991.
[11] M. Moghaddam and W. C. Chew, “Nonlinear two-dimensional velocity pro­
file inversion using time domain data,” IEEE Trans. Geosci. Remote Sensing,
vol. 30, Jan. 1992.
[12] M. Moghaddam and W. C. Chew, “Study of some practical issues in inver­
sion with the Bom iterative method using time-domain data,” IE EE Trans.
Antennas Propagat., vol. 41, no. 2, pp. 177-184, 1993.
[13] M. Moghaddam, “Forward and inverse scattering problems in the time do­
main,” Ph.D. dissertation, University of Illinois at Urbana-Champaign, 1991.
[14] A. Tarantola, “The seismic reflection inverse problem,” in Inverse Problems of
Acoustic and Elastic Waves, F. Santosa, Y. H. Pao, W. Symes, and C. Holland,
Eds., Philadelphia: SIAM, 1984.
43
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[15] J. E. Dennis and R. B. Schnabel, Numerical Methods for Unconstrained Opti­
mization and Non-linear Equations. Englewood Cliffs: Prentice-Hall, 1983.
[16] E. Polak, Computational Methods in Optimization. New York: Academic Press,
1971.
[17] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical
Recipes in FORTRAN, 2nd Ed. New York: Cambridge University Press, 1992.
[18] A. N. Tikhonov, “On the problems with approximately specified information,”
in Ill-Posed Problems in the Natural Sciences, A. N. Tikhonov and A. V. Goncharsky, Eds., Moscow: MIR Publishers, 1987.
[19] A. V. Goncharsky, “Ill-posed problems and their solution methods,” in Ill-Posed
Problems in the Natural Sciences, A. N. Tikhonov and A. V. Goncharsky, Eds.,
Moscow: MIR Publishers, 1987.
[20] K. S. Yee, “Numerical solution of initial boundary value problems involving
Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat,
vol. AP-14, pp. 302-307, 1966.
[21] A. Taflove, “Review of the formulation and applications of the finite-difference
time-domain method for numerical modeling of electromagnetic wave interac­
tions with arbitrary structures,” Wave Motion, vol. 10, pp. 547-582, 1988.
[22] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic
waves,” J. Comput. Phys., Mar. 1994.
[23] W. C. Chew and W. H. Weedon, “A 3-D perfectly matched medium from modi­
fied Maxwell’s equations with stretched coordinates,” Microwave Opt. Technol.
Lett., 1994. Submitted for publication.
[24] T. J. Cavicchi and W. D. O’Brien, “Numerical study of high-order diffraction
tomography via the sine basis moment method,” Ultrason. Imaging, vol. 11,
pp. 42-74, 1989.
44
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 3
LOCAL S H A P E F U N C T IO N (LSF) M E T H O D
In this chapter, we present another nonlinear inverse scattering algorithm that
uses a local shape function (LSF) approximation to parameterize very strong scatterers in the presence of a transient excitation source. The local shape function
(LSF) method was presented previously for CW excitation and derived using a Tmatrix formulation [1-4]. Derivation in the time domain is not so straight-forward
because the wave equation describing the scattering process becomes undefined as
either the permittivity or conductivity approaches infinity. Consequently, it is only
appropriate to describe the scattering from perfect conductors in terms of boundary
conditions on each of the scatterers. A time-domain LSF algorithm is developed
here using the T-matrix formulation and interpreting the results in the time domain
[5].
The time-domain LSF formulation allows the inverse scattering problem to be
solved with a finite-difference time-domain (FDTD) forward solver [6 , 7]. The lo­
cal shape function may be implemented as a volumetric boundary condition in the
FDTD algorithm. The inverse scattering problem is then cast as a nonlinear opti­
mization problem similar to the DBIM where the
dimensional Frechet derivative
of the scattered field is computed as a single forward propagation of the FDTD
forward solver. The Frechet transposed operator is similarly computed as a backpropagation and correlation using the FDTD forward solver.
The LSF formulation that we present here is sufficiently general th at inverse
scattering may be performed when many metallic scatterers are involved and when
the scatterers are of arbitrary shape, size and location. Although much work has
been done for imaging dielectric scatterers, there has been relatively little work
done for imaging metallic scatterers. Severed previous methods [8-13] parameterize
the shape of the surface of scatterers. However, these methods generally require
knowledge of the location of the scatterer and cannot handle multiple scatterers.
45
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.1
A T -m atrix Formulation
We first consider the scattering due to a single CW excitation source in a region
O possibly containing metallic scatterers. The scattering region O is discretized into
N subvolumes Vi, i = 1, 2 , . . . ,N and a binary local shape function 7 * is assigned
to each subregion V. such that
-y. =
(31)
7l
10,
V in ^ = 0
1
'
where S represents the total volume occupied by the metallic scatterers. Now,
assuming that the subregions V. are sufficiently small, we may approximate the
scattering solution as the field emanating from N filamental metallic scatterers, one
located at the center of each of the regions V{. The total field may be written in a
discretized form as
N
<j>(r) = tf(r a) ea + Y , I)ir i) ai
(3*2)
i= 1
where the first term corresponds to the incident field (the field generated with
no scatterers present), and the second term represents the scattered field. The
monopole wave functions used to expand the incident and scattered fields are given
as V’f c ) = H ^ ik o lr —r^j) where r t is the location of the ith scatterer, r is the
observation point and r , = r —r[. Similarly, r s = r —r's where r '3 is the location
of the illuminating source field. We assume th at the incident field is generated by a
filamental line source oriented along the 2-axis so that the problem is entirely twodimensional (see Figure 3.1). Note that the shape of the scatterer is not necessary
in writing (3.2). A formulation is presented below such th at a* = 0 when j i = 0.
This formulation is valid for both TM and TE polarized incident waves.
The boundary condition at each scatterer Vi may now be written as
4>sca.(r i) = 7i^i(l)^inc(r i)i
* = lj 2 ,..., JV
(3.3)
where Tj(i) is the single-scatterer transition coefficient for Vi. Imposing the multiple
scattering boundary condition of (3.3) on each of the scatterers, we have [2, 3]
(
N
O-ia
\
“i" y & ij Qj
j- 1
3&
(3.4)
)
46
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
o
F ig u re 3.1 Discretization of a metallic scatterer S on a finite difference
grid. The shaded region indicates 7 * = 1. Inversion is performed in
object region O.
47
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where c e ij represents a coordinate transformation from the j th to the ith scatterer
and Tt(i) is a single-scatter transition (or reflection) coefficient. In deriving (3.4)
from (3.3), we have essentially redefined the incident field as the total field impinging
on the ith scatterer and the scattered field as that emanating from the ith scatterer
only.
We now convert (3.4) to a matrix-vector notation to simplify the problem. Defin­
ing
[o], = a*,
[A
] ij
=
i = 1,2, ... , 1V
o c ij (1 - S i j )
[ T ] ij = l i T n \ ) & i j ,
i, j =
i, j
1,2 ,...,
(3.5)
N
(3.6)
= 1 ,2,...,N
(3.7)
[®s]i =
0^'®)
Equation (3.4) becomes
a = y • (ea + A -a)
which reduces to
(3.9)
_
_
a = ( I - t • A )~ x - 7 • es
=
7
• ( I - A ■7 ) - 1 • es
(3.10)
= T 9-
Equation (3.10) relates the incident field vector es to the scattering amplitude vector
a . The ith component of the vector g appearing in Equation (3.10) represents the
incident field on the ith scatterer before it is multiplied by the shape function 7 *.
We will refer to g as the “ghost field” because the field inside the scatterer that
may be located at V i may or may not exist depending on the value of 7 ,.
The time-domain equivalent of Equation (3.10) may be conveniently imple­
mented using an FDTD forward solver. At each time step after the ghost field g
is computed, the product 7 • g is computed by applying the boundary condition
specified by Equation (3.3) inside the scatterer.
48
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.2
F rech et D eriv a tiv e O perator
The above formulation expresses the scattered field 0Sca(**) as a function of the
shape function 7 . To solve the inverse problem using a gradient search approach,
we have to know how the scattered field changes as a function of a change in 7 . In
other words, we need the Frechet derivative.
To obtain the Frechet derivative of the scattered field, we let 7 = j 0 + 87 and
a = cio + 8a in Equation (3.10). We also relax the assumption that 7 is a binary
variable, i.e., we let 7 be a continuous real variable over the interval [0,1]. The value
7 0 can be thought of as a background shape function and <$7 a small perturbation in
the background value, with the corresponding background and perturbed scattered
field amplitudes a 0 and 8a.
Using a binomial approximation to Equation (3.10), we can show that
5o« ( J -
70
-A ) _1 • £ 7 • ( I —.4. • 7o ) - 1 *
(3.H)
The scattered field may then be written as
<tyscaM = '»l>t ’8a
= il>t - ( I - 7 0 ’ A r 1 ' 8i ’ ( I - A - j 0) - 1 . e a
= ^ * . (I where
),
70
(3.12)
• A ) - 1 ■8 7 -g
i = 1, 2 , . . . , N .
We now define the frequency-domain Frechet derivative F as a semi-discrete
linear operator mapping a set of discrete shape functions in a model space M into
the continuous frequency-domain measurement space V ,
(3.13)
The Frechet derivative F maps the perturbation 87 € M intothe corresponding
change in the scattered field 8<f>€ V ,
to c a (r) = F{<$7}.
(3.14)
49
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The operation of the Frechet derivative F acting on Sy may be computed from
Equation (3.12). In Equation (3.12), the ghost field g is first multiplied by the
shape function perturbation £ 7 and then operated on by tj}1•(I —y 0 •A ) -1 . Clearly,
this operator maps a field internal to the object out to the receivers, but its exact
physical interpretation is not obvious at first.
We again consider the scattering problem from N volumetric scatterers V*, but
instead of a single line source we assume that there are N filamental line sources in
the same positions as the scatterers Vi. The total field is then (similar to Equation
(3.2))
N
<P(r ) =
N
13^(r«) *i +t13
^(ri)
=l
i= l
(3-15)
where s, and a* represent the source and scattered field amplitudes, respectively, at
the ith scatterer. We can show that the scattered field may now be written as
<f>(r)
(7 —7 0 -A ) _1
s
(3.16)
where [s]i = Sf.
Comparing Equation (3.16) with Equation (3.12), we see that the Frechet deriva­
tive is simply a forward propagation of the induced source <£7 • g back out to the
receivers. This operation may be conveniently implemented in the time domain us­
ing an FDTD forward solver with a volumetric boundary condition. First the ghost
field is computed inside the object space by running the FDTD code and storing
the total (ghost) field at each time step before the volumetric boundary condition
is applied. Then the ghost field is multiplied by 6 7 , treated as an induced source
inside the object, and forward propagated back out to the receivers. Care must be
taken when including this induced source because, assuming a background shape
function 7 0, this induced source must be added after the boundary condition is
applied at each time step so that the induced source is not modified. The value
of 5<j>sea thus obtained is identical to the scattered field that one would obtain by
solving the forward scattering problem in the presence of scatterer £ 7 in the limit
£7
—» 0 .
In the continuous limit where N —»■00 and ^ -)• 0, Equation (3.12) may be
50
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
written in integral form. Noting that V’i =
H ^ik o lr'i — r '|) ea, Equation (3.12) becomes
—*•*[) and [e3]i = a tsea =
fysc*(r) = J d r'i h 2(r,r'i) 6y(r<) h i(r',r's) ea(k0).
(3.17)
where the frequency dependence of the source ea has been included. In the above, hi
is a propagator mapping the source ea into the ghost field internal to the scatterers,
while h-i represents the propagator mapping the induced source S'f • g out to the
receivers. From Equation (3.12) it is also apparent that hi and hi are reciprocal
operators, or that
h 2 (r,p ') = h i(r ',r ).
(3.18)
Equation (3.17) is an expression for the continuous-space frequency-domain Frechet
derivative T . We define this operator formally as the operator mapping functions
in the continuous model space M. into the continuous frequency-domain data space
V,
F - .M v ^ V .
(3.19)
5 0 sc a k )= ^ { 5 7 k )} -
(3.20)
Hence, we have that
Comparison of with Equation (17) with the corresponding Frechet derivative op­
erator in the DBIM method shows that both algorithms allow the Frechet derivative
to be computed as a forward propagation of a forward scattering solver. In Equa­
tion (17) the ghost field, analogous to the incident field in the DBIM, must be
computed first and then product of the ghost field and Sy(r') is propagated out to
the receivers. Hence, new LSF algorithm and the DBIM algorithms derived in the
previous chapter have sim ilar structure, but differences lie in their implementation.
3.3
Frechet T ransposed O perator
In the previous section, we derived the Frechet derivative operator T that maps
a perturbation Sy(r) into the corresponding change in scattered field 6<f>sca(r). Also
needed to solve the inverse problem is the Frechet transposed operator F 1 that maps
$<t>sc&(T) back into the shape function space,
(3.21)
51
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The frequency-domain transposed operator is defined by the relation
< f A » ( r ) , ^ { # y (r)} ) v = { f 7 (r), ^ { W ) } ) M
(3-22)
where the first inner product is over the data space and the second is over the model
space.
Although the Frechet operator was derived using a frequency-domain T-matrix
formulation in the previous section, it may be written in the time domain by Fourier
transforming Equation (3.17). From Equation (3.17), we have
5<j>Sca.{rm , r n , t ) = J dr' J d t'h 2(rm, r ', t - t') S y ( r ')
x [ d t"h \{r',T n,t' - t")sn{rn,t")
(3.23)
J
= J d r' J dt' h 2(rm, r \ t - t')5'y{r')<j>g (r', r„ , f')
where sn (t) is the time waveform transmitted from the n th source and (f>g( r ',r n,t')
is the ghost field inside the object. Note that in Equation (3.23), 5<j)sc& is still a
continuous function of space, but is assumed to be evaluated at receiver location
r m. The transm itter location r n is a parameter.
The kernel hi appearing in Equation (3.23) corresponds to the computation of
the time-domain ghost field, while the kernel h 2 corresponds to the propagation, in
the presence of background shape function 70 , of the induced source
back out
to the receivers. Intuitively, one can see th at the kernel h 2 is the reciprocal kernel
of hi because while hi is the field inside the object before the boundary condition
is applied, h 2 isthe field at the receivers due to a source inside the object that is
applied after the boundary condition is applied. Hence, as in Equation (3.18),
h 2(r,r',< ) = h i(r;,r ,f ) .
(3-24)
We now define the time-domain Frechet operator T that maps shape functions
defined on the continuous model space M. into the continuous time-domain data
space 2>,
T :M
V.
52
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.25)
Equation (3.23) may be expressed as
^sc a(r m) r n,t) = F { 8'y{r')}
=j
dr' K{Tm,Tny ■ ,t) 8'r{r>)
where now T is the time-domain Frechet derivative. In Equation (3.25),
K (rm ,T n,r ',t) is the kernel of the integral operator and may be obtained from
Equation (3.23) as
K { r m, r ny , t ) = J d t 'h 2 (rmy , t - l f ) 4>g(T',Tn,t?).
(3.26)
The Frechet transposed operator maps a perturbation in the scattered field back
into the object space,
J* : T > ^ M .
(3.27)
This operator may be expressed as
* 7 ( 0 = ^{^sca(T*m ,rn,t)}.
(3.28)
The time-domain equivalent of Equation (3.22) is then
( 8<f>sea.(vrni Tn, t), ^ { * 7 ( 0 } )g = ( * 7 (O i ^*{*0sca(**m, **»»<)} ) M .
(3.29)
Since T 1 is also an integral operator, we may express Equation (3.28) in the general
form
Nt
Nr
*7(0 = £ £
,
_
I d t K t { r ',r m, r n,t) 8<i>sc*{1*77111*711t)
(3.30)
n=l m =l
where now K l { r \ r m, r n, t) is the kernel of the transposed operator and N? and N r
are the numbers of transmitters and receivers. However, from Equation (3.29) it
follows directly that JF and J* have the same kernel, a general property of transposed
operators [14], or that
-K'(r,, r m, 7*n,i) = Jift ( r ', r m, r n,t).
(3.31)
Substituting Equation (3.26) into Equation (3.30), we have
NT
f
Nr
f
_
* 7 ( 0 = £ J dt* <j>g( r ',r n,t') Y l I d th 2(rm, r ', t - t') 8<j>sca(rm, r n,t).
n=l
(3.32)
m =1
53
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Equation (3.32) is not very useful as written because h 2(rm, r ', t ) maps a function
in r ' to r m, but 6<f>sc&(rm, r n,t) is a function of r m. Also, we have to perform a
change of variables in time because the last integral to the right of Equation (3.32)
is a time-reversed convolution. Hence, we let r = T —t, r ' = T —t' and use Equation
(3.24) to replace h% with its reciprocal kernel to obtain
Nt
fT
S7(r ') - Z ) [
n = 1 J°
dT><t>g(r\ rn, T N
r
,T
x V /
m=1 J°
r')
drfei(r/,rTO,r/ -
t ) 5<j>scak{rm, r n, T - r ) .
Equation (3.33) is the desired equation for computing the Frechet transposed
operation. The integral to the right of Equation (3.33) is a backpropagation of
the time-reversed field S<f>sca.{rm, r n,T — r) from the receivers back into the object
space. Since the kernel hi is specified in the integral, the backpropagated field at
each time step in the FDTD algorithm is stored prior to the application of the
volumetric boundary condition, forming, in effect, a “backpropagated ghost field.”
This backpropagated field is computed and summed for all of the receiver locations
and correlated with a time-reversed version of the original ghost field inside the
object with a final sum over the transmitter locations. Again, the structure of
Equation (3.33) is similar to those derived in the previous chapter for computing
the DBIM Frechet derivatives, but the implementation is different.
3 .4
Inverse S olu tion U sin g C onjugate G radient
The inverse problem, that of reconstructing the shape function y from measure­
ments of the scattered field <f>meas, may be posed as the problem of minimizing a
certain cost function 5(7) with respect to 7. For the LSF method, we define the
cost function at the kth iteration as
S f c ( 7 ) = ^ | | 0 - < £m eas||5 +
= \{<t> ~ tineas)* '
\\\lk ~
7 fc-llU
‘ (<t>~ 4>meas) + |CTk ~ I k - 1)* ' (ik ~ Tfc-l)
(3.34)
54
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where || • ||^ and || • ||.m indicate norms in the data and model spaces. The diagonal
matrix C jj is used to weight the data corresponding to various receiver locations
and may be used as a “boosting procedure” in the limited angle problem. The
parameter e is a “tuning parameter” used to regularize the inverse problem. Note
th at a minimal amount of a priori information is used here.
The functional defined by Equation (3.34) may be minimized numerically using
a conjugate gradient algorithm similar to the one presented in the previous chapter.
Note again that the Frechet derivative T is used only in computing the parameter
update step size, whereas the Frechet transposed operator 7* is used to compute
the gradient which defines the search direction.
3.5
C o m p u ter Sim u lation R esu lts
The LSF algorithm was verified by reconstructing time-domain data obtained
from the inverse Fourier transform of data obtained from a two-cylinder frequencydomain 2-D TM-polarization T-matrix code [15]. Two geometries were used for the
inverse scattering simulation. In geometry (1) shown in Figure 3.2, the scattering
object O was considered to be contained inside a 214.4 mm x 214.4 mm (27 x 27
grid points) region. The object O and all transmitters were contained inside the
FDTD grid. The background medium was assumed to be th at of free space and an
absorbing (radiation) boundary condition was used at the edges of the grid. The
transmitted waveform had a —3 dB cutoff of 1.22 GHz, corresponding to a minimum
wavelength Amjn = 245.0 mm.
The scattering data of two metallic scatterers, each of radius a = 7.942 mm
(0.0645 Am in ) and separated by d = 112.3 mm (0.458 Amin ) were computed using the
analytic formulation. The reconstruction obtained from the LSF method is shown
in Figure 3.3 and that obtained using the DBIM with a conductivity optimization
is shown in Figure 3.4. No convergent solution was obtainable using DBIM with
permittivity optimization. In Figure 3.3, the two metallic scatterers are easily
distinguishable, while those in Figure 3.4 are not resolved.
55
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
FDTD grid boundary
F ig u re 3.2 Geometry (1) for bistatic inverse scattering simulation.
Object region is 214.4 mm x 214.4 mm (27 x 27 grid points), dt =
127.1 mm (0.517 Am in ) , d n = 135.0 mm (0.549 Am jn ) . Entire geometry
including transmitters and receivers is contained within the FDTD
grid.
56
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(a) Original Object
(b) Shape Function Reconstruction
F ig u re 3.3 Local shape function reconstruction of two metallic cylin­
ders, each of radius a = 7 .9 4 2 mm ( 0 .0 6 4 5 Am jn ) , separation d =
1 1 2 . 3 m m ( 0 .4 5 8 Am i„ ) using geometry ( 1 ) shown in Figure 3 . 2 . The
LSF algorithm was applied for 2 0 iterations and no regularization was
used.
(a) Original Object
(b) Conductivity Reconstruction
%
F ig u re 3.4 Conductivity reconstruction using the distorted Bom it­
erative method of the same object as used in Figure 3.3. The DBIM
was applied for 20 iterations and no regularization was used.
57
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The relative residual error (RRE), defined as
1/2
(0
$ zn e a s)
'
‘ (0
0meas)
(3.35)
' V'meas
1
is shown in Figure 3.5 for each iteration step. Only twenty iteration steps were
computed and it is clear from Figure 3.5 that while both algorithms converged, the
LSF algorithm has a much faster convergence than the DBIM with conductivity
optimization. The DBIM solution leveled off with a RRE of about (0.5), whereas
the LSF solution leveled off at about (0.15).
RRE =
Although the inverse solution formulation allows for regularization of the so­
lution, no regularization was used in any of the reconstructions presented here.
Regularization is not needed here because the time-domain inverse solution is suf­
ficiently well-posed th at the inverse mapping contains a small null space. Also, the
inversion was performed for only a small number of iterations (20 ) compared to the
number of unknowns (625) and the small eigenvalues do not have a chance to cause
instability. In our experience, it is always best to use the least amount of regular­
ization necessary to obtain a convergent solution because additional regularization
only causes slower convergence.
Figure 3.6 shows a LSF reconstruction of two metallic scatterers that are of the
same radius as those used in Figure 3.3, but the separation has been reduced to
d = 84.2 mm (0.344 Amjn). The two cylinders are still clearly resolved, but it is
apparent th at this separation is at about the resolution limit of this algorithm.
In geometry (1) shown in Figure 3.2, the closest receiver was placed at 0.55lAmin.
To demonstrate that the reconstruction shown in Figure 3.3 can also be obtained
in a regime where evanescent waves are negligible, geometry (2) shown in Figure
3.6 was used whereby the transmitters and receivers were moved outward such that
the closest receiver was located at approximately 3 Amin. To make the problem
computationally feasible, the transmitters and receivers were moved outside the
FDTD grid and propagated into the grid using an analytic formulation. Figure 3.8
shows the resultant reconstruction. While there is some resolution loss in the size
of the scatterers, the two scatterers are still well-resolved.
58
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
s0 5
0.9
'— '
0.8
o
t °-7
W 0.6
0.5
C
O
0)
0.4 •
Q)
0.3-
>
0.2
0.1
Iteration Number
F ig u re 3.5 Comparison of relative residual error (RUE) for LSF (solid
curve) and DBIM (dashed curve) reconstructions shown in Figures 3.3
and 3.4.
59
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(a) Original Object
(b) Shape Function Reconstruction
F ig u re 3.6 Local shape function reconstruction of two metallic cylin­
ders more closely spaced than those shown in Figure 3.3. Cylinders
are of radius a = 7.942 m m (0.0645 Amin ), separation d = 84.2 mm
(0.344 Amin) using geometry (1) shown in Figure 3.2. The LSF algo­
rithm was applied for 20 iterations and no regularization was used.
60
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
R
R
Y
E
E
E
E E
FD TD g rid b o u n d ary
E
E
E
F ig u re 3.7 Geometry (2) for bistatic inverse scattering simulation.
Object region is 214.4 m m x 214.4 mm (27 x 27 grid points), d r —
762.4 m m (3.10 Amin), dR = 810.1 m m (3.29 Amin). Transmitters and
receivers are placed outside the FDTD grid and propagated in using
an analytic formulation.
61
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(a) Original Object
(b) Shape Function Reconstruction
*
■jf
F ig u re 3.8 Local shape function reconstruction of two metallic cylin­
ders, each of radius a = 7.942 mm (.0645 Amjn), separation d = 112.3
mm (.458 Amin) using geometry (2) shown in Figure 3.7. The LSF al­
gorithm was applied for 20 iterations and no regularization was used.
62
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.6
Conclusions
We presented a new inverse scattering method for time-domain inverse scatter­
ing that uses the local shape function method to parameterize very strong scatterers.
This new method was derived from a similar frequency-domain algorithm by inter­
preting the results in the time domain and implementing the local shape function
as a volumetric boundary condition in an FDTD forward solver.
The new time-domain LSF algorithm is similar in structure to the DBIM al­
gorithm presented previously [16, 17], but the new algorithm uses the local shape
function approximation rather than the distorted Bom approximation. In terms
of implementation, the two algorithms differ in the computation of the forward
scattering solution and Frechet derivative and transposed operators.
The computer simulation results show that while a convergent solution may be
obtained using a distorted Bom iterative algorithm with conductivity optimization,
the convergence using DBIM is extremely slow compared to LSF. A higher resolution
solution was obtained using the new LSF method algorithm for the same number
of iterations.
3 .7
R eferen ces
[1] W. C. Chew and G. P. Otto, “Microwave imaging of multiple conducting cylin­
ders using local shape functions,” IEEE Microwave Guided Wave Lett., vol. 2,
pp. 284-286, July 1992.
[2] G. P. Otto and W. C. Chew, “Microwave inverse scattering-local shape function
imaging for improved resolution of strong scatterers,” IEEE Trans. Microwave
Theory Tech., vol. 42, no. 1, pp. 137-141, 1994.
[3] G. P. Otto and W. C. Chew, “Inverse scattering of H z waves using local shape
function imaging: A T-matrix formulation,” Int. J. Imaging Syst. Technol.,
1994. Accepted for publication.
[4] G. P. Otto, “Electromagnetic inverse scattering problems for strong scatterers,”
Ph.D. dissertation, University of Illinois at Urbana-Champaign, 1993.
[5] W. H. Weedon and W. C. Chew, “Time-domain inverse scattering using the
local shape function (LSF) method,” Inverse Probl., vol. 9, pp. 551-564, 1993.
63
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[6 ] K. S. Yee, “Numerical solution of initial boundary value problems involving
Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat.,
vol. AP-14, pp. 302-307, 1966.
[7] A. Taflove, “Review of the formulation and applications of the finite-difference
time-domain method for numerical modeling of electromagnetic wave interac­
tions with arbitrary structures,” Wave Motion, vol. 10, pp. 547-582, 1988.
[8] D. Colton, “The inverse electromagnetic scattering problem for a perfectly con­
ducting cylinder,” IEEE Trans. Antennas Propagat., vol. AP-29, no. 2, pp. 364368, 1981.
[9] D. Colton and P. Monk, “A novel method for solving the inverse scattering
problem for time-harmonic acoustic waves in the resonance region II,” SIAM
J. Appl. Math., vol. 46, no. 3, pp. 506-523, 1986.
[10] A. Roger, “Newton-Kantorovitch algorithm applied to an electromagnetic in­
verse problem,” IEEE Trans. Antennas Propagat., vol. AP-29, no. 2 , pp. 232238, 1981.
[11] C. C. Chiu and Y. W. Kiang, “Inverse scattering of a buried conducting cylin­
der,” Inverse Probl., vol. 7, pp. 187-202, 1991.
[12] C. C. Chiu and Y. W. Kiang, “Microwave imaging of multiple conducting
cylinders,” IEEE Trans. Antennas Propagat., vol. 40, no. 8 , pp. 933-941,1992.
[13] P. Maponi, L. Misici, and F. Zirilli, “An inverse problem for the three dimen­
sional vector helmholtz equation for a perfectly conducting obstacle,” Comput­
ers Math. Applic., vol. 22, no. 4, pp. 137-146, 1991.
[14] A. Tarantola, “The seismic reflection inverse problem,” in Inverse Problems of
Acoustic and Elastic Waves, F. Santosa, Y. H. Pao, W. Symes, and C. Holland,
Eds., Philadelphia: SIAM, 1984.
[15] W. C. Chew, Waves and Fields in Inhomogeneous Media. New York: Van
Nostrand, 1990.
[16] W. C. Chew and Y.-M. Wang, “Reconstruction of tw o -d im en sio n a l permittiv­
ity using the distorted Bom iterative method,” IEEE Trans. Med. Imaging,
vol. MI-9, no. 2 , pp. 218-225, 1990.
[17] M. Moghaddam, W. C. Chew, and M. Oristaglio, “Comparison of the Bom
iterative method and Tarantola’s method for an electromagnetic time-domain
inverse problem,” Int. J. Imaging Syst. Technol., vol. 3, pp. 318-333, 1991.
64
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 4
B R O A D B A N D MICROWAVE IN V E R SE
SCATTERING M EA SU R EM EN TS
Up to this point in this thesis, we have been concerned primarily with inverse
scattering imaging algorithms. Another important part of inverse scattering is
the measurement techniques. Most theoretical works on inverse scattering assume
that measurements can be made with infinite signal-to-noise ratio and ignore many
sources of measurement error and modeling errors that would be present in a prac­
tical situation. This chapter is devoted to broadband microwave inverse scattering
measurement technology. We will discuss several sources of measurement and mod­
eling errors present in practical measurement situations and ways to reduce the
measurement errors. Step-frequency radar (SFR) and time-domain impulse radar
measurement systems are discussed and compared, and reconstructions are shown
using real experimental data collected with a prototype SFR system.
High-quality inverse scattering measurements are actually quite difficult to ob­
tain primarily because scattered fields resulting from unknown objects are usually
weak compared to the exciting incident field. Broadband microwave scattering mea­
surements are even more difficult because all of the system components including
signal sources, receiver electronics, signal separators, microwave switches, anten­
nas and transmission lines must operate over a broad bandwidth. In addition, the
frequency-dependent transfer function of the transmission lines and antennas in the
measurement system must be accounted for.
There are two basic approaches to obtain broadband scattering measurements.
One approach is to collect data in the frequency domain using a vector (magnitude
and phase) network analyzer or step-frequency radar system. The data collected
from a SFR system could then be converted to the time domain by means of a dis­
crete Fourier transform. The other basic approach is to collect data directly in the
time domain using an impulse radar consisting of a pulse generator and a data col65
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
lection system employing a sampling head such as an digital sampling oscilloscope.
Both types of data collection systems have their advantages and disadvantages, and
the choice of one over the other depends on the application.
Some of the advantages of frequency-domain SFR data collection over the cor­
responding impulse radar system include:
(1) Higher signal-to-noise ratio (SNR) due to the reduced bandwidth.
Time-domain data collection must operate at full bandwidth in
order to collect the scattered pulse with little distortion. Stepfrequency radar systems employing superheterodyne receiver de­
signs use narrowband intermediate-frequency amplifiers and elec­
tronics and therefore have much higher SNR [1, 2].
(2) Very stable signal sources are available for CW measurements such
as synthesized sources that drift by as little as 10_ 9/H z/day and
have frequency resolutions between 1 and 4 Hz. Nonsynthesized
sources have stabilities and frequency resolution on the order of
±10 MHz [3, Chap. 4]. Time-domain measurements, on the other
hand, are subject to timing jitter due to inaccuracies in the trigger
and time base circuitry. This problem is troublesome because the
jitter varies with time, and therefore signed averaging is of limited
effectiveness for reducing jitter. Time-domain data collection sys­
tems also suffer from zero level drift and zero ambiguity [3, Chap.
11].
(3) Various waveform shapes may be simulated and spectral shaping
may be used by filtering the frequency-domain data collected from
an SFR system. Frequency-domain filtering may be performed
with data collected from a pulsed system, but is generally more
difficult [4-6].
(4) Many systematic errors due to reflections in transmission lines
and connectors may be removed easily using network analyzer
66
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
calibration techniques [3, 7].
(5) A large percentage of the total source power is available at each
measurement frequency in a SFR system, whereas the source
power is spread out over the pulse bandwidth in a time-domain
data collection system.
Impulse radar systems do offer some advantages over SFR systems:
( 1) Measurement time is generally reduced.
(2) It may be easier to make a portable, cost efficient impulse radar
system using time domain data collection techniques. This argu­
ment, however, is open for some debate.
(3) Time domain aliasing that is present in SFR systems due to the
frequency spacing is not present in impulse radar systems. The
maximum time window is determined by the pulse repetition rate
rather than the frequency spacing.
For laboratory measurements at least, it would seem that an SFR system is
preferable to an impulse radar system for broadband microwave measurements due
to the increased accuracy and flexibility. Most commercial ground-penetrating radar
systems, however, use impulse radar. This is most likely due to the ability to
make compact, portable, cost-efficient impulse radar systems. W ith the increasing
demands for accuracy in areas such as nondestructive evaluation, SFR systems are
likely to play a larger role in the future.
This chapter discusses three different broadband microwave inverse scattering
measurement systems. First, a prototype SFR system developed at the University
of Illinois is discussed, and image reconstructions are shown for both metallic and
dielectric scattering objects from real experimental data. A commercial monostatic
impulse radar system purchased from Geophysical Survey Systems, Inc. is also dis­
cussed. Finally, we discuss a rudimentary bistatic impulse radar system constructed
67
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
from a picosecond step generator, broadband antennas and a sampling oscilloscope.
4.1
S tep-F requency R ad ar (S F R ) M icrow ave M easu rem en t S y stem
4.1.1
Description of the SFR system
prototype broadband microwave inverse scattering measurement system has
been designed and built to demonstrate the practicality and effectiveness of stepfrequency radar inverse scattering for nondestructive evaluation (NDE). A block
diagram of basic components of the prototype SFR measurement apparatus is shown
in Figure 4.1. The system consists of a broadband switched antenna array, an
HP 8510B automated network analyzer, microwave switches and controller, and
an optional broadband amplifier. The entire measurement system is automated
and controlled by a computer workstation. Custom software was written in the
C programming language to control the measurement apparatus via an IEEE-488
(GPIB) interface.
A
The H P 8 5 1 0 B automated network analyzer serves as both the transmitter and
receiver and allows us to collect both amplitude and phase information by stepping
through various frequencies. A breakdown of the various system components in the
H P 8 5 1 0 B is shown in Figure 4.2 . The network analyzer contains the H P 8 5 1 0 B
m a.infra.mft unit, an H P 8 5 1 4 A S-parameter test set, and an 8 3 6 2 1 A synthesized
sweeper. A photograph of the H P 8 5 1 0 B is shown in Figure 4 .3 . The synthesized
sweeper is not visible in the photograph.
The test set is a two-port device that separates incoming and outgoing waves
and allows one to measure any of the various parameters £ u , £ 21, £ l 25 £ 22- The pa­
rameters £ n and £22 are used if monostatic measurements are desired and measure
the reflected signal in ports 1 and 2 , respectively, due to a wave transmitted from
the same port. The parameters £21 and £12 measure the signal received in port
2 due to a signal transmitted from port 1, and vice versa, and are used if bistatic
measurements are desired.
The custom control software allows the HP 83621A synthesized frequency sweeper
68
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
HP 8510B
Network
Analyzer
(.5-18 GHz)
Data
IN
Computer
Workstation
Control
OUT
Port 2 (RX)
Broadband
Amplifier
Port 1 (TX)
RX Microwave
Switch (DC-18 GHz)
DC Power
Supply
TX Microwave
Switch (DC-18 GHz)
Switch
Controller
Broadband
Antenna Array
(2-12 GHz)
Scattering
Objects
F ig u re 4.1 Block diagram of prototype step-frequency radar (SFR)
broadband inverse scattering measurement system.
69
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
o
HP 8510B
Mainframe
Cfl
3
HP8514A
Test Set
CO
00
00
■'tIfm
w
w
Port 2
Port 1
•
•
HP 83621A
Synthesized
Sweeper
F ig u re 4.2 Block diagram of various system components of the HP
8510B vector network analyzer.
70
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F ig u re 4.3
Photograph of HP 8510B automated network analyzer.
71
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
to be operated in one of three modes. In the first mode, a ramp sweep mode is used
th at employs the “lock and roll” technique [3, Chap. 4] in which phase-locking
occurs only for the first frequency and the other frequency points are collected in
rapid succession. This is the fastest data collection mode available, but the least
accurate. The second mode, known as the step mode, performs phase-locking at
each individual frequency. This mode is more accurate than the ramp mode, but
takes more time for the data collection. For both the ramp and step modes, the
frequency sweep/step is performed automatically by the HP 8510B. However, these
modes axe both limited to a fixed number of frequency points: 51, 101, 201, 401
or 801. A third mode has been added in the custom software that performs the
frequency stepping via software and allows an arbitrary number of frequency points.
This mode has the same accuracy as the step mode, since phase-locking is performed
at each frequency, but takes longer than even the automatic step mode to collect
data. For most inverse scattering measurements, the ramp sweep mode employing
the “lock and roll” technique provides sufficient accuracy. This allows for accurate,
rapid inverse scattering measurements.
Time averaging may be performed within the HP 8510B to reduce random
noise. We have found that the averaging does not significantly increase the data
collection time since the averaging is performed in the HP 8410B firmware. For
bistatic measurements, an optional broadband low-noise preamplifier may be placed
in series with the receiver port to improve the dynamic range of the system.
4.1.2
Broadband Vivaldi antenna element
To collect broadband scattering data, an antenna array with broadband radiat­
ing elements is required. Several types of broadband antennas are available for this
purpose including bowtie antennas [8, 9], spiral antennas [10], broadband horns
and tapered slotline antennas [11-14].
The bowtie antenna is used frequently
in ground-penetrating radar systems, but has an operating bandwidth of approxi­
mately only 2:1. Furthermore, bowtie antennas have a low directional gain, similar
to a dipole antennas. Usually, the bowtie antennas axe arranged in an array or have
a resonant cavity to boost the gain. Spiral antennas provide a circularly polarized
72
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
wave, as opposed to the linear polarized wave that we desire here. They are difficult
to fabricate. Broadband, horn antennas provide linearly polarized waves, have high
gain, and broad bandwidth, but are usually very expensive.
In the work presented here, a tapered slotline antenna known as the Vivaldi
antenna [12] was used. The Vivaldi antenna has an exponentially tapered slot,
provides a linearly polarized wave, has higher directional gain than a bowtie antenna
and a bandwidth of roughly 6:1. In addition, the Vivaldi antennas are inexpensive
and can be fabricated easily using a metalized substrate material [11].
A pair of prototype Vivaldi antennas were provided to us by Professor Paul
Mayes and Dr. David Tammen at the University of Illinois, Urbana-Champaign.
One of these prototypes is shown in Figure 4.4. This particular design employs a
wedge feed for performing an impedance transformation from the coaxial feed to
the slotline. The design used by Mayes and Tammen was duplicated and refined in
order to create an array of Vivaldi elements. The photos in Figure 4.5 show one of
the new Vivaldi elements in the final stages of the fabrication process.
4.1.3 Early prototype antenna array
For microwave inverse scattering measurements, independent measurements are
required at transmitter and receiver locations that are placed at many points around
the scattering object. A particular arrangement of transmitter and receiver loca­
tions is known as a “look angle.” As a general rule, the more look angles th at are
used and the broader the range of look angles, the better the reconstruction that
may be obtained from the inverse scattering algorithm.
For inverse scattering measurements, it is preferable to have several antenna
elements arranged in a fixed antenna array. Unlike a conventional antenna array
in which many antenna elements are excited simultaneously, the inverse scattering
measurement array operates in a switched mode in which only one transmitting
element is excited and one receiving element is active. By employing broadband
microwave switches in the array to multiplex the transmitting and receiving ele­
ments, a full suite of look angles can be measured automatically without the need
73
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F ig u re 4.4 Photograph of broadband Vivaldi antenna designed and
build by Prof. Paul Mayes and Dr. David Tammen at the University
of Illinois, Urbana-Champaign.
74
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(a)
F ig u re 4.5 Photographs of new Vivaldi elements based on the design
of Mayes and Tammen in final stages of fabrication process.
75
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
for repositioning of the antenna or disconnecting and reconnecting elements.
Some early experimental measurements were made using a primitive antenna
array constructed for this purpose. The early array consisted of a wooden frame,
shown in Figure 4.6, with 9 slots cut in the frame in which the antenna element
locations are arranged in a linear fashion. Two 2-12 GHz broadband Vivaldi an­
tennas were used as the transmitting and receiving elements and were moved by
hand between the various slots to obtain the various measurement configurations.
Microwave switches were not employed in this early prototype since the array was
used for proof of principle only.
The primitive array is capable of operating in either a bistatic mode, in which
separate antennas were used as the transmitter and receiver, or a monostatic mode
in which the same antenna element operates as both a transmitter and receiver.
For bistatic measurements, 4 slots were used as transmitter locations and 5 slots
were used as receiver locations totaling 20 different measurement configurations.
The transmitter and receiver locations were interlaced in the bistatic mode and
separated by 2.0 in, or 5.08 cm, as shown in Figure 4.7. In the monostatic mode
of operation, all 9 transmitter/receiver locations were used to collect measurement
data.
The antenna elements in the array were aimed at a fixed range of 20.0 cm,
where the center of the object was expected to be located. The reason for aiming
the antenna elements at a fixed position, rather than aiming them straight ahead,
was simply to maximize the signal-to-noise ratio at all measurement frequencies. A
general principle of directed radiators is that the higher the frequency of operation,
for the same radiator, the narrower the beamwidth of the main radiation lobe. We
found by experience that by originally aiming all of the antenna elements straight
ahead, the elements in the center of the array would be aimed directly at the fixed
range of 20.0 cm, but the elements at the edge of the array would be directed away
from the fixed range. As a result, the high-frequency components of the received
signal from the end elements were smaller than those of the center elements. The
end elements experienced a broadening, or smoothing of the received pulse obtained
76
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F ig u re 4.6 Photograph of early antenna array consisting of a slotted
wood frame and two broadband Vivaldi antennas that were moved
between the various slots by hand.
77
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
40.6 cm
5.08 cm
19.8 cm
8.89 cm
8.89 cm
F ig u re 4.7 Arrangement of tr a n sm itte r s, receivers and object grid us­
ing early antenna array.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
from Fourier transforming the SFR data. This pulse broadening was due solely to
the frequency dependence of the directional gain of the antenna elements.
In principle, the pulse broadening effects can be removed by inserting a lownoise preamplifier in tandem with the receiving element, measuring the frequency
dependence of the directional gain of the antennas, and calibrating out the effect.
Aiming the antennas at a fixed direction is a cheap and simple way to boost the SNR
and ensure that all of the antenna elements have the same frequency-dependent gain.
Of course, this solution is valid only for scatterers located in the vicinity of the range
location th at the antennas are aimed. But it is better to aim the antenna elements
at the center of the object space, or the location where the object is expected to
occur, rather than aiming the elements straight ahead.
4.1.4
New switched antenna array
Recently a new 11-element broadband switched antenna array, shown in Figure
4.8, was designed and built. This new antenna array contains 11 identical 2-12 GHz
Vivaldi antennas arranged in a linear array and two DC-18 GHz SP6T microwave
switches th at are computer controlled. One switch is connected to 5 array elements
while 6 elements are connected to the other. Hence, the array may be configured via
computer control to operate as either an 11-element monostatic array or a multi­
bistatic array consisting of 30 different measurements. The switches automatically
terminate the antenna elements at 50 fi when they are switched off, reducing the
coupling among the elements. The antenna elements and microwave switches are
enclosed in a polystyrene housing.
The arrangement of transmitters and receivers for new switched antenna array
is shown in Figure 4.9. The separation of the antenna elements in the new antenna
array has been increased from 5.08 cm to 8.0 cm, giving a total baseline of 80.0
cm. This baseline length was chosen with the goal in mind of resolving objects at
a range of R = 40.0 cm. As in the early prototype array, the antennas were aimed
at a fixed range, here given as R = 40.0 cm.
79
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F ig u re 4.8 Photograph of new 11-element broadband switched an­
tenna array containing 11 identical broadband Vivaldi antennas and
two microwave switches enclosed in a polystyrene housing.
80
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
80.0 cm
8.0 cm
8.5 cm
8.5 cm
F ig u re 4.9 Arrangement of transmitters, receivers and object grid for
new switched antenna array. Note: dimensions are correct in the
figure, but drawing not to scale.
81
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4.2
O verview of M e asu rem en t D a ta P rocessing
The various steps in the data processing of the measurement data using both
the DBIM and LSF algorithms are summarized in the flow diagram of Figure 4.10.
The switch on the left indicates that either measured data or computer generated
scattering data (synthetic data) may be used in the inverse scattering algorithm.
The algorithm begins with specification of the initial parameters, which are set to
zero because we wish to use a minimal amount of a priori information. Using the
current computer model, forward scattering data are generated and subtracted from
the measured data. This difference is then used to compute a measure of the resid­
ual field error. If the difference is below a specified tolerance, the current model
parameters are displayed on a graphics workstation. If the field error is not below a
specified tolerance, the field error is sent to a conjugate gradient optimization pro­
cedure which returns an update to the model parameters. The process is repeated
until a convergent solution is attained.
4.3
4.3.1
R e c o n stru c tio n s U sing E x p e rim e n ta l S F R D a ta
Metallic object reconstructions using early array
Early experimental results were obtained using the primitive wooden frame
antenna array shown in Figure 4.6. The arrangement of transmitters, receivers and
the object grid for the early experiment is as shown in Figure 4.7. The measurement
consists of 4 transmitter locations indicated by a “T” in Figure 4.7, and 5 receiver
locations indicated by an “R.” For the object functions and reconstructions shown
below, the object space consists of the 35 x 35 grid point (g.p.) region indicated by
a dashed box in Figure 4.7. A grid space step of size A x = 2.54 mm and time step
size A t -- 5.5 ps were used in the FDTD code. The source pulse used to drive the
FDTD forward solver was obtained from a Fourier transform of the product of the
gains of the two Vivaldi antennas at boresight, measured in the antenna test range.
We note th at in all of the multibistatic measurements using the early system,
data were collected using only the first two transmitter locations, but all receiver
locations. These data were then reflected about the center of the array and copied
82
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Initial
Parameters
Measured Data
(Frequency-Domain)
Inverse
FFT
Computer
Model
Compute
Forward Scattering
Solution
Se(r), 8o(r), 8y(r)
Parameter
Update
Conjugate
Gradient
Optimization
,
Computer
Generated
Data
(Time-Domain)
Compute
Residual
Field Error
,
| Frechet
i Transposed
| Operator
!-
Frechet
Derivative
Operator
Error
Below
Tolerance?
Display
Model Parameters
s(r), a(r), y(r)
F ig u re 4.10 Block diagram of processing of measured and computer
simulated scattering data.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
to the data space for the two other transmitter locations. This was done because
the measurement process was extremely tedious since the transmitting and receiver
antennas had to be disconnected, moved to different slots in the array, and then
reconnected for each measurement. An entire measurement took approximately 40
m in with this symmetrization. This also limited us to measuring objects that were
se lf-sy m m e tr ic about the center of the array.
Four different cases of metallic scattering objects were considered, all consisting
of two metallic cylinders. First, we examined two metallic cylinders aligned hori­
zontally, or parallel to the array. The cylinders were separated by 3.2 cm, each with
a diameter of 0.4 cm. The reconstruction is shown in Figure 4.11 and it is clear th at
we can resolve this case very well. The second case that we considered consisted of
the same cylinders with 3.2 cm separation, but in this case aligned vertically. The
reconstruction is shown in Figure 4.12, and both cylinders axe still well-resolved.
Next, we used the same cylinders, but reduced the separation to 1.5 cm. First,
we looked at the horizontal arrangement, shown in Figure 4.13 and then the vertical
arrangement, shown in Figure 4.14. Note that in both cases, an image is obtained,
but we cannot distinguish the two scatterers from each other. This particular
config uration of scatterers is just beyond the resolution of our experiment at present.
4.3.2
Metallic object reconstructions using new array
The experimental results performed using the early antenna array were repeated
using the new switched antenna array for comparison purposes. The measurement
geometry used in the new array is shown in Figure 4.9. The object space again
consisted of a 35 x 35 subgrid, but with A x = 2.5 mm and A t = 5.5 ps. Again, the
source pulse used to drive the FDTD forward solver was obtained from a Fourier
transform of the product of the gains of the two Vivaldi antennas a t boresight,
measured in the antenna test range.
Figures 4.15 and 4.16 show the results for the horizontal and vertical alignment
with 3.2 cm spacing while Figures 4.17 and 4.18 show the results with 1.5 cm
spacing. The new array seems to give a better resolution, but the 1.5 cm case is
84
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Shape Function Reconstruction
Original Object
F ig u re 4.11 Original object and shape function reconstruction of two
metallic cylinders of diameter 0.4 cm aligned horizontally with sepa­
ration 3.2 cm. D ata were collected using early antenna array.
Shape Function Reconstruction
Original Object
F ig u re 4.12 Original object and shape function reconstruction of two
metallic cylinders of diameter 0.4 cm aligned vertically with separation
3.2 cm. D ata were collected using early antenna array.
85
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Original Object
Shape Function Reconstruction
F ig u re 4.13 Original object and shape function reconstruction of two
metallic cylinders of diameter 0.4 cm aligned horizontally with sepa­
ration 1.5 cm. Data were collected using early antenna array.
Original Object
Shape Function Reconstruction
F ig u re 4.14 Original object and shape function reconstruction of two
metallic cylinders of diameter 0.4 cm aligned vertically with separation
1.5 cm. Data were collected using early antenna array.
86
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
still not resolved.
No symmetrization of the data were done with any of the data collected from
the new array. Apparently, the symmetrization used to produce the images from
the early antenna array did not add an unfair advantage to the reconstructions.
4.3.3
Dielectric object reconstructions using new array
Metallic objects are more difficult for inverse scattering algorithms to image
than dielectric objects because the inverse scattering problem is more nonlinear for
metallic objects. However, dielectric objects are more difficult to measure because
the scattered field produced by a dielectric object is much weaker than th at of a
metallic object. We present reconstructions of dielectric objects in this section to
demonstrate that accurate scattering data can be collected from dielectric objects
and th at high-quality images may be generated.
Figures 4.19 and 4.20 show reconstructions of plastic PVC pipes of diameters 2.7
cm and 4.8 cm, respectively. Both of these pipes were located in an air background.
The DBIM permittivity optimization algorithm was used for both cases and the
algorithm was terminated at 60 iterations. For both pipes, high quality images
were produced. The bottoms of the pipes are not reconstructed as well as the tops
because scattering data was collected from the top only.
Figure 4.21 shows a DBIM permittivity reconstruction of an empty glass grad­
uated cylinder of diameter 5.25 cm. An good reconstruction was obtained after 60
iterations. The graduated cylinder was located in air.
4 .4
C om m ercial Im p u lse R adar D a ta C ollection S ystem
Another approach to collect broadband microwave scattering is to use an im­
pulse radar data collection system. One commercially available impulse radar sys­
tem is produced by Geophysical Survey Systems, Inc. (GSSI) in Hudson, NH. Figure
4.22 show a picture of the mainframe unit of the GSSI System 8 impulse radar sys­
tem. This system has been used at the University of Illinois, Urbana-Champaign
87
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Range, cm
Original Object
Shape Function Reconstruction
16
16 ■
18
18
20
a
£002 0
?
•
•
^
Ife;
0
':
i
-x
c
CO
OS
22
22
3#
24
24
- 4 - 2
0
- 4 - 2
2
■%
.......
0
2
Transverse Axis, cm
Transverse Axis, cm
F ig u re 4.15 Original object and shape function reconstruction of two
metallic cylinders of diameter 0.4 cm aligned horizontally with sepa­
ration 3.2 cm. D ata were collected using new switched antenna array.
Shape Function Reconstruction
Original Object
16
16
I18
18
Range, cm
M
•
ao
20
&20
00
c
es
04
•
22
'
22
JA
*'
24
24
- 4 - 2
0
- 4 - 2
2
0
2
Transverse Axis, cm
Transverse Axis, cm
F ig u re 4.16 Original object and shape function reconstruction of two
metallic cylinders of diameter 0.4 cm aligned vertically with separation
3.2 cm. D ata were collected using new switched antenna array.
88
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
*
Shape Function Reconstruction
Original Object
Range, cm
16
-4
- 4 - 2
0
2
Transverse Axis, cm
-2
0
2
Transverse Axis, cm
4
F ig u re 4.17 Original object and shape function reconstruction of two
metallic cylinders of diameter 0.4 cm aligned horizontally with sepa­
ration 1.5 cm. D ata were collected using new switched antenna array.
Shape Function Reconstruction
Original Object
Range, cm
16
- 4 - 2
0
2
Transverse Axis, cm
- 4 - 2
0
2
Transverse Axis, cm
F ig u re 4.18 Original object and shape function reconstruction of two
metallic cylinders of diameter 0.4 cm aligned vertically with separation
1.5 cm. Data were collected using new switched antenna array.
89
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
DBIM Permittivity Reconstruction
16
18
&20
s(3
22
24
-4
-2
0
2
Transverse Axis, cm
4
F ig u re 4.19 DBIM permittivity reconstruction after 60 iterations of
a hollow plastic PVC pipe of diameter 2.7 cm located in air.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
DBIM Permittivity Reconstruction
.k
0
-* m » -
5
Transverse Axis, cm
F ig u re 4.20 DBIM permittivity reconstruction after 60 iterations of
a hollow plastic PVC pipe of diameter 4.8 cm located in air.
91
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
DBIM Permittivity Reconstruction
15
S
o
&20
a03
QC
25
-5
0
5
Transverse Axis, cm
F ig u re 4.21 DBIM permittivity reconstruction after 60 iterations of
an empty glass graduated cylinder of diameter 5.25 cm located in air.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
for the nondestructive evaluation of civil structures [15,16]. Figure (4.23) shows the
GSSI System 8 being used in the laboratory to examine a concrete test specimen.
Figure 4.24 shows a closeup of the 900 MHz transducer provided with the impulse
radar system and Figure 4.25 shows a concrete specimen with a styrofoam filled
void th at was examined. We do not present any results from the GSSI system since
extensive data from this system has been processed using diffraction tomography
techniques [15].
Attempts to process the data from the GSSI system using nonlinear inverse
scattering methods [16] were of limited success because the measurement system
was too crude and contained many modeling errors that could not be removed. For
example, there was no way to measure the source pulse or the antenna gain in order
to incorporate their effect into the computer model. Nonlinear inverse scattering
methods have the potential to achieve much higher resolution and invert larger
contrasts than diffra c tio n tomography methods, but the measurement data must
be of sufficient quality for this advantage to be realized.
4 .5
R u d im en ta ry B ista tic Im p u lse R adar S ystem
To learn more about time-domain data collection, a PSPL 4050B step generator
was purchased from Picosecond Pulse Labs, Inc. in Boulder, CO. The pulse gener­
ator shown in Figure 4.26, along with the remote pulse head shown in Figure 4.27,
has a 10 Volt amplitude and a 45 ps risetime with a timing jitter of only 1.5 ps RMS.
This pulse generator is well-suited for impulse radar because of its high-amplitude
pulse and extremely stable timing circuitry.
A rudimentary bistatic impulse radar system was assembled using the PSPL
4050B step generator, the broadband Vivaldi antennas, and an old Tektronics sam­
pling oscilloscope that was available in the laboratory. The Tektronix model 7603
oscilloscope, shown in Figure 4.28, contained a 7S-12 sampling unit th at included
an S-6 sampling head and an S-53 trigger recognizer. The PSPL step generator and
remote sampling head were connected to a broadband Vivaldi antenna to generate a
high-amplitude pulse. A second Vivaldi antenna was connected to the S-6 sampling
93
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F ig u re 4.23 GSSI System 8 system being used at the University of
Illinois, Urbana-Champaign to test concrete specimens for defects.
94
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F ig u re 4.24
MHz.
GSSI System 8 transducer with center frequency of 900
95
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F ig u re 4.26
Picosecond Pulse Labs model 4050B step generator.
F ig u re 4.27 Remote pulse head for PSPL 4050B step generator, shown
connected to broadband Vivaldi antenna.
96
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
head with the loop-thru connection terminated at 50 ft.
Attempts to measure a transmitted pulse from the bistatic impulse radar ar­
rangement were not very successful due to the noise and timing jitter of the sam­
pling oscilloscope. The noise and jitter effects of the scope were so bad th at it
was not possible to measure a stable waveform without using the “high-resolution”
smoothing feature of the sampling scope. The waveform was fairly stable with the
high-resolution turned on, but the waveform measured under this condition was
extremely distorted. Apparently, this feature reduces the sampling efficiency below
its normal 100% value by lowering the forward gain in the sampling postamplifier,
and it is recommended that it not be used [17]. Another problem with using the
Tektronix sampling oscilloscope for an impulse radar is th at the horizontal offset
had a very limited range, and many times it was difficult or impossible to shift the
pulse so that it was visible on the CRT display. The 7603 also did not have an
external computer interface from which to extract data.
Due to the problems listed above with the 7603 sampling oscilloscope, this
project had to be terminated until funds could be obtained to purchase a new
digital sampling oscilloscope. Newer digital sampling scopes such as the HewlettPackard HP 54120 Series scopes and Tektronics 11801A scope with a Tektronics
SD-22 or SD-26 sampling head are much more stable and contain digital averag­
ing capabilities. Furthermore, the newer scopes have built-in IEEE-488 computer
interfaces for extracting data.
4 .6
C onclusions
We have discussed two basic approaches for broadband inverse scattering data
collection: step-frequency radar and impulse radar. A prototype SFR data collec­
tion system developed at the University of Illinois was presented and reconstruc­
tions of both metallic and dielectric objects were shown from data collected with
the SFR system. Two impulse radar systems, a commercial ground-penetrating im­
pulse radar system as well as a rudimentary system constructed from a picosecond
step generator and sampling oscilloscope, were discussed. Advantages and disad97
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F ig u re 4.28 Tektronix model 7603 oscilloscope with 7S-12 sampling
unit containing S-6 sampling head and S-53 trigger recognizer.
98
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
vantages of the three systems were discussed and particular attention was paid to
the accuracy of the measurements.
Two prototype broadband antenna arrays were shown along with the SFR sys­
tem. Both the early prototype array and the newer switched array presented were
linear arrays. The array shape chosen was strictly for convenience for NDE mea­
surements of objects beneath flat surfaces. These designs could easily be modified
to semicircular arrays for measuring objects such as cylindrical concrete support
beams. These linear arrays could also be used for 3-D measurements by making
measurements at several locations in the vertical direction.
To have a significant impact on technology, inverse scattering measurement
technology should be developed at a rate commensurate with the development of
inverse scattering imaging algorithms. Furthermore, in developing inverse scattering
im a g in g a lg o r ith m s , we should keep in mind the many sources of error th at would be
present in a practical measurement situation so that the effects of the unavoidable
measurement errors can be minimized.
4 .7
R eferen ces
[1] M. I. Skolnik, Introduction to Radar Systems. New York: McGraw-Hill, 1980.
[2] S. Haykin, Communication Systems, 2nd Ed. New York: Wiley, 1983.
[3] G. H. Bryant, Principles of Microwave Measurements. London: Peter Peregrinus Ltd., 1988.
[4] S. Riad, “The deconvolution problem: An overview,” Proc. IEEE, vol. 74, no. 1,
p p . 82-85, 1986.
[5] N. S. Nahman, “Software correction of pulse measurement data,” in High Speed
Electrical and Optical Pulse Measurements, J. E. Thompson, Ed., pp. 351-417,
Boston: Martinus Wijhoff, 1986.
[6] J. A. Landt, “Typical time domain measurement configurations,” in Time Do­
main Measurements in Electromagnetics, E. K. Miller, Ed., New York: Van
Nostrand Reinhold, 1986.
[7] W. Kruppa and K. F. Sodomsky IEEE Trans. Microwave Theory Tech., Jan.
1971.
[8] G. H. Brown, “Experimentally determined radiation characteristics of conical
and triangular antennas,” RCA Rev., Dec. 1952.
99
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[9] H. J. Jasik, Antenna Engineering Handbook. New York: McGraw-Hill, 1961.
[10] J. D. Dyson, “The equiangular-spiral antenna,” IR E Trans. Antennas Propag.,
vol. AP-7, pp. 181-187, 1959.
[11] K. M. Frantz, “An investigation of the Vivaldi flared radiator.” M.S. thesis,
University of Illinois at Urbana-Champaign, 1992.
[12] P. J. Gibson, “The Vivaldi aerial,” in Ninth European Microwave Conference,
(Brighton, U.K.), 1979.
[13] K. L. Kollberg, T. Thungren, and K. S. Yngvesson, “Vivaldi antenna/Snline circuit for SIS mixers,” in Sixth International Conference on IR and MM
Waves, (Miami Beach, FL), 1981.
[14] K. S. Yngvesson, D. H. Schaubert, T. L. Korzeniowski, E. L. Kollberg, T. Thun­
gren, and J. F. Johansson, “Endfire tapered slot antennas on dielectric sub­
strates,” IEEE Trans. Antennas Propagat., vol. AP-33, no. 12, pp. 1392-1400,
1985.
[15] J. E. Mast, “Microwave pulse-echo radar imaging for the nondestructive eval­
uation of civil structures.” Ph.D. dissertation, University of Illinois at UrbanaChampaign, 1993.
[16] W. H. Weedon, J. E. Mast, W. C. Chew, H. Lee, and J. P. Murtha, “Inver­
sion of real transient radar data using the distorted-Bom iterative algorithm,”
in IEEE Antennas and Propagation Society International Symposium Digest,
(Chicago, IL), 1992.
[17] J. R. Andrews, “Comparison of ultra-fast rise sampling oscilloscopes.” Applic.
Note AN-2a, Picosecond Pulse Labs, 1989.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 5
E F F IC IE N T F D T D C O M P U T A T IO N O N A
M A S S IV E L Y P A R A L L E L S U P E R C O M P U T E R
The finite-difference time-domain method [1, 2] is widely regarded as one of the
most popular computational electromagnetics algorithms. Although FDTD is con­
ceptually very simple and relatively easy to program, the method is actually quite
efficient since it involves 0 ( N 1-5) computational complexity in 2-D and 0 ( N 1,33)
computational complexity in 3-D [3]. In fact, FDTD can be considered an optimal
algorithm since 0 ( N a) numbers are produced in 0 ( N a) operations.
FDTD is also ideally suited for implementation on a single-instruction multipledata (SIMD) massively parallel computer because the stencil operations that must
be computed at each node of the space grid involve only nearest-neighbor interac­
tions and may be implemented at a m in im u m communication cost [4]. A major
challenge, however, is in implementing absorbing boundary conditions (ABCs) at
the edges of the FDTD grid. On scalar and vector computers, these boundary con­
ditions are typically computed using methods such as the Engquist-Majda [5], Mur
[6], Liao [7] or Higdon [8] ABC. However, these methods are not ideal for parallel
supercomputers since they all involve communication with many elements normal to
the grid boundary. Such communication can easily surpass the time spent comput­
ing core FDTD operations in the grid interior, especially for higher-order boundary
conditions, and hence cam become a bottleneck in the FDTD code. Also, they do
not allow for SIMD operation on a parallel machine without the use of masking.
An alternate method of implementing an ABC is to use a conventional absorbing
material boundary [4, 9-14] . For SIMD parallel computation, these methods have
the advantage th at the ABC may be implemented with the same FDTD stencil
operation as the interior nodes by modifying the conductivity material parameter
at the edge of the FDTD grid. The disadvantage is that the reflection coefficient
at the absorbing border is zero only at normal incidence and is both angle and
101
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
frequency dependent. Consequently, the absorbing material border region must
be made quite large—typically 20-100 grid points along each edge to minimize
reflections.
Recently, Berenger [15] suggested a more general method of implementing an ab­
sorbing material boundary condition. Berenger proposed a procedure for 2-D wave
propagation whereby the Maxwell equations are generalized and added degrees of
freedom axe introduced. The added degrees of freedom allow the specification of
absorbing interfaces with zero reflection coefficient at all angles of incidence and
all frequencies. Moreover, the generalized Maxwell equations reduce to the famil­
iar Maxwell equations as a special case and hence the same generalized equations
can be used to propagate fields in both the interior region and absorbing region.
Although the interface between the interior region and the absorbing boundary is
reflectionless, there is still a reflection from the edge of the grid. The advantage of
using Berenger’s procedure is that much larger conductivity values may be specified
in the absorbing region, leading to a drastic reduction in the number of grid points
required for the absorbing boundary.
A formulation similar to the Berenger idea is derived for 3-D wave propagation
from first principles using a coordinate stretching approach. The advantage of the
new method for SIMD parallel computation is stressed. The method is validated
with 3-D FDTD numerical computations on a Thinking Machines Corporation Con­
nection Machine CM-5.
5.1
M od ified M a x w ell E quations
For a general medium, we define the modified Maxwell equations in the fre­
quency domain, assuming e~tut time dependence, as
V e x £ = iujjiH
(5.1)
Vh x H = —iueE
(5.2)
Va • eE = p
(5.3)
Ve • fiH = 0
(5.4)
102
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where
_
. 1 9
. 1 9
.i a
.
v ’ = xTxTx +vTvTy + z7xTz'
( 5
.
' 5 )
and
h
Xhx d x
y hy dy
hz dz
^
^
In the above, ez,ft*, i = x ,y ,z are coordinate-stretching variables that stretch the
x , y , z coordinatesfor Ve and V^. It shall be shown later th at when e, and hi are
complex numbers, the medium can be lossy. Note that (5.3) and (5.4) are derivable
from (5.1) and (5.2). A general plane wave solution to Equations (5.1) - (5.4) has
the form
E = E 0 eikr
(5.7)
and
H = H 0 eikv
(5.8)
where k = x k x + yky + zkz. Substituting Equations (5.7) and (5.8) into Equations
(5.1) and (5.2) above gives
k e x E = uifiH
(5-9)
kh x H = —ueE ,
where fce = x ^ -1- y ^ + z ^ and kh = x j£ + y
have
—lo ueH = k e x kh x H
(5.10)
+ zj£ . Combining the above, we
5.11
= k h(ke - H ) - ( k e - k h)H .
But from Equation (5.9), k e • H = 0 for a homogeneous medium. This gives the
dispersion relation
a}2/j,€ = k e -kh
(5-12)
or
k2
=
~ \~ k l
exhx
+ —r - f c y +
eyhy y
~ \ ~ ks
ezh~
(5 -1 3 )
where «2 = a
Equation (5.13) is the equation of an ellipsoid in 3-D and is satisfied by
kx —
exhx sin 0 cos <f>,
(5.14)
103
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(5 .1 5 )
and
(5.16)
Note th at when a , hi, i = x, y, z are complex, the wave in the x, y, and z directions
are attenuative and can be independently controlled. Under the matching condition,
ex = hx, ey = hy, and e2 = h~, we have |fce]2 = \kh\2 = k2. The wave impedance is
then given by
(5.17)
irrespective of the values for et-, i = x ,y ,z and the direction of propagation.
5 .2
S in gle Interface P rob lem
Assume that a plane wave is obliquely incident on the interface z = 0 in Figure
5.1. Furthermore, we may assume th at the plane wave is of arbitrary polarization.
The incident field may be decomposed into a sum of two components, one with
electric field transverse to z (TE2) and the other with magnetic field transverse to
z (TMZ). We will examine these two components individually.
In the (TE2) case, we let the incident field in region 1 be given as
E i = E 0 eik i'r
(5.18)
In the above, khi • E q = 0, and E q is in the xy plane. Similarly, we define the
reflected field in region 1 as
E r = R ^ E o r eikr'r
(5.19)
and the transmitted field in region 2 as
E t = T t e E o t eik t'r
(5.20)
Phase matching requires that kiX = krx = ktx and k,y — kry = kty. Hence, we can
define E qt = Eot = E q since they ail point in the same direction. Applying the
104
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
incident
wave
reflected
wave
region 1
region 2
z-0
transmitted
wave
P^6
F igu re 5. 1 Plane wave with arbitrary polarization incident on the
plane z = 0 .
105
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
boundary condition that the tangential electric field must be continuous across the
plane z = 0, we have
l + i2TE = TTE.
(5.21)
This boundary condition follows from the modified Maxwell Equation (5.1).
The magnetic field may be determined using Equation (5.9) for regions 1 and 2
as
H l = fe»« X E 0 eik , r + ^ T E fcre X E 0 ^
ujfii
and
H 2 = r TE kte X E ° eikt r
0Jfi2
(5.23)
where k;e — x ^c*+ y C
^y+ z ez
^ and similarly for k re and k te. We also define k iz = ki~,
k2z = ktz and note that krz = —k\z. Then equating the tangential components of
Equations (5.22) and (5.23), we have
k \ze2z[x2 [l - -ftTE] = r TE k2ze\z(i\.
(5.24)
Combining Equations (5.21) and (5.24), we have
# te
k ize2zfi2 - k2ze i z fii
_
^
^
k ise2zfjt2 4- k2ze\sn\
and
/yiTE_
Tte =
____________
.
k ize2zV>
2 + k2ze\~H\
(5 .2 6 )
Applying a similar procedure to the TMZ component, we have
jjT M
_
k i z h 2z e2 — k 2z h i ~ e i
k\zh2ze2 + k2zh\ze\
and
T ™ = ____ 2klzh2z*2--------.
k izh2zt2 + k2zh\~€\
5.3
(5.28)
P erfectly M atch ed Interface
The phase matching condition requires that k \x = k2x and kiy = k2y, or
Kiy/eixhix sm0icos<f>i = k 2\Je2xh2x sind2 cos 4>2
(5.29)
106
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
and
K \^e \y h iy sin 61 sin fa — k 2yje 2y/i2y sin 02 sin <7^2
(5.30)
where «i = Uy/fiiex and « 2 = For a perfectly matched medium, we
choose
Ci = e2,H\ = fJ-2, ex = ft* and ey = hy. Equations (5.29) and (5.30) become
eia, sin 0i cos 4>x =
* sin 62 cos 0 2
(5.31)
eiysin 0 isin 0 i = e2y sin 02 sin 02 .
(5.32)
62
and
If we now choose e\x = e2x and eiy = e2y, then 9\ = 02, 0 i = 02 and we can show
th at both R t e = 0 and i2™ = 0 for all angles of incidence and all frequencies.
If region 1 is a vacuum, then n = n 0, e = eo> and
(eia,,eiy,e i2,
hl2) =
(l,l,l,l,l,l).
(5.33)
In order to have a lossy region 2 with no reflections at the region 1/region 2 interface,
we choose
(^2*5 C2y>C2X, h 2x,/l 2y, h2xJ = (1,1, S2, 1,1, £2)
(5.34)
where s2 is a complex number. In this case,
kix = k 2x = «o sin 0 cos 0
(5.35)
kiy = fc2y = Ko sin 6 sin 0
(5.36)
k \2 = Ko cos 0
(5.37)
k 2z = KoS2COS0
(5.38)
where kq = Uy/fioeQ. If s2 = s '2 + is2» the wave will attenuate in the 2 direction.
This kind of interface is useful for building material ABCs in an FDTD simulation.
5.4
M o d ified M a x w ell E quations in th e T im e D om ain
For the general case of a matched medium, we let ex = hx = sx, ey = hy = sy
and ez = hz = sx. Then, Vc = V*. =
+ y ± j^ +
In Equation (5.1),
we write the curl as
1
a
^Q
Ve x £ = — — x x E H — y x E
sx ox
Syay
1
^
— z x E.
s~ oz
107
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(5.39)
Then, defining H„x, H 3y, and H 3z in terms of the components of Equation (5.39),
we let
iu>nH3z = — — x x E
(5.40)
vX
. „
i a .
(5.41)
^ H '- “ 7t a ,? x E
„
and
,,
1 o ^
tujfiHaz — —— z x E
sz dz
where H = H Sx + H Sy + H 3z. Similarly, we can write Equation (5.2) as
. „
1 d „
__
—zueE. = — — x x H
sx ox
r,
and
1
„
-zu>eE3z =
^ - xH
,,
(5.42)
(5.43)
(5.44)
o
1
.
—— 2
x IT,
(5.45)
sr
where E = E Sx + E Sy + E 8z. Note that H Si, E Si, i = x ,y ,z are two-component
vectors.
We now let sx = 1 + iax/we, sy = 1 + iay/u e and sz = 1 4- iaz/ue. Writing
Equations (5.40) - (5.42) and (5.43) - (5.45) in the time domain, we have
(5.46)
(5.47)
dJEf.
<r2u
5
1 a t ' + « * ‘. = - a z ‘ * E
(5.48)
(5.49)
(5.50)
dEa
d
£ at + ' ■ * - - a . * x j r
(5.51)
108
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Equations (5.46) - (5.51) describe 3-D wave propagation in a perfectly matched
medium. The wave propagation phenomenon described by these equations is very
similar to that described by the Maxwell equations with the exception that atten­
uation may be controlled through the ax, <ry and az variables. The FDTD imple­
mentation of these equations on a Yee FDTD grid is straightforward. Absorbing
boundaries at the edges of the simulation region may be created by choosing ap­
propriate values of crx , ay and az. Equations (5.46) - (5.51) may be seen to include
Berenger’s equations [15] as a subset for the 2-D TE or TM case.
The above equations involve 12 components of electromagnetic fields. For a
free-space/lossy medium interface, a scheme may be devised using only 10 field
components for the 3-D case, and only 3 components for the 2-D case. However,
this is achieved at the loss of SIMD operation on a parallel machine.
5.5
C om p u ter Sim ulation R esu lts
To demonstrate the new method, a 3-D orthogonal grid FDTD algorithm was de­
veloped based on Equations (5.46) - (5.51). The FDTD algorithm was implemented
as a SIMD code on the Thinking Machines Corporation Connection Machine CM5. The algorithm operates very efficiently on the CM-5 because the FDTD stencil
operations that have to be computed at each node involve only nearest-neighbor
interactions. The communication operations resulting from the nearest-neighbor
interactions are at a minimum cost since the neighboring processors are for the
most part at the bottom of the fat-tree communication network, where communi­
cation bandwidth is maximum.
To validate our 3-D FDTD algorithm, we solved a simple problem of computing
the field radiated from an infinitesimal electric dipole in free space. An analytic
solution was also computed in the frequency domain for many excitation frequen­
cies. The frequency-domain solution was then multiplied by the spectrum of FDTD
source pulse and inverse Fom.c. transformed to yield a time-domain analytic solu­
tion for comparison with the FDTD solution.
The FDTD solution was solved in a cubic region of dimension (Nx,N y,N .) =
109
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(128,128,32) grid points. The grid parameters chosen were Ax = Ay = Az =
2.5 mm, A t = 4.5 ps and N t — 512 time steps were computed.
The infinitesimal electric dipole was simulated by exciting the E y field in a single
grid cell with the source pulse
(5.52)
A xA yA z
where r = 1/47t/ o and a value of fo = 1.0 GHz was chosen. The dipole source
was located at grid location (nx,n y,n z) = (91,64,16). The E x and E y fields were
obtained by sampling the fields at grid location (nx,n y,n z) = (37,91,16).
The absorbing boundaries used for the FDTD simulation consisted of planar
layers of 8 grid points thickness on all surfaces. Along the borders parallel to the
x-axis, the value of <rx was specified, while ay and <rz were specified on the borders
parallel to the y- and 2-axes, respectively. The conductivity values were chosen
with a parabolic taper decreasing from the maximum value towards the center of
the grid such that the reflection coefficient at normal incidence was R q = 0.0001.
The E x field computed using both the analytic formulation and the FDTD
algorithm are overlaid in Figure 5.2. The curves due to the analytic and numeri­
cal solutions are barely distinguishable, indicating excellent agreement. Sim ilarly,
the E y fields due to the analytic and numerical solutions are overlaid in Figure
5.3. Again, we see excellent agreement. Any difference between the analytic and
numerical solutions in Figures 5.2 and 5.3 may be attributed to modeling errors
such as the finite size of the dipole source and the discrete approximation of the
Maxwell equations in addition to reflections due to imperfections in the absorbing
boundaries.
The CM-5 machine used to solve the FDTD problem is located at the National
Center for Supercomputing Applications (NCSA) at the University of Illinois. The
program was written in CM Fortran and compiled using CMF version 2.1. The
CM-5 at the NCSA has 512 nodes with vector units. The CPU times were deter­
mined by running the problem on 32, 64, 128 and 256 node partitions. For this
problem, a total of 1.5 million unknown field quantities (128 x 128 x 32 grid x 3
110
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4000
A m p litu d e
of Ex
fie ld
2000
-2000
-4000
-6000
-8000 -
-10000
100
200
300
400
500
T im e s te p
F ig u re 5.2 Analytic and numerical FDTD solutions overlaid for the
E x field resulting from an infinitesimal electric dipole.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
600
x 10
A m p litu d e
of Ey f i e l d
0.5
-0.5
-1.5.
100
200
300
400
500
T im e ste p
F ig u re 5.3 Analytic and numerical FDTD solutions overlaid for the
E y field resulting from an infinitesimal electric dipole.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
600
field components) were determined for 512 time steps. The CPU times axe shown
in Table 5.1.
T able 5.1: CPU times for FDTD Problem on CM-5
Nodes
32
64
128
256
5.6
CPU sec (Run 1, Run 2, Run 3, Avg.)
50.5,50.2,50.6; 50.4
29.9,30.0,30.0; 30.0
17.9,18.4,18.4; 18.2
12.4,13.2,12.7; 12.8
C onclusions
A modified set of Maxwell equations has been introduced using complex coordi­
nate stretching factors along the three cartesian coordinate axes. This modification
introduces additional degrees of freedom in the Maxwell equations such that absorb­
ing boundaries may be specified with zero reflection coefficient at all frequencies and
all angles of incidence. The formulation was shown to be related to the perfectly
matched layer that was recently derived by Berenger for 2-D wave propagation.
A 3-D FDTD algorithm was developed from the modified Maxwell equations th at
uses the reflectionless absorbing interface property to implement radiation bound­
ary conditions at the edges of the FDTD grid. The accuracy of the algorithm was
validated by computing the field radiated from an infinitesimal electric dipole and
comparing against a known analytical expression. The FDTD algorithm was imple­
mented on the Connection Machine CM-5 and timing results were presented. This
breakthrough in absorbing material boundary conditions allows EM scattering to be
computed very efficiently on SIMD parallel computers. This fast 3-D forward solver
will allow us to solve 3-D inverse scattering problems very quickly on a massively
parallel supercomputer.
113
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
5.7
References
[1] K. S. Yee, “Numerical solution of initial boundary value problems involving
Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat.,
vol. AP-14, pp. 302-307, 1966.
[2] A. Taflove, “Review of the formulation and applications of the finite-difference
time-domain method for numerical modeling of electromagnetic wave interac­
tions with arbitrary structures,” Wave Motion, vol. 10, pp. 547-582, 1988.
[3] W. C. Chew, Waves and Fields in Inhomogeneous Media. New York: Van
Nostrand, 1990.
[4] W. H. Weedon, W. C. Chew, and C. M. Rappaport, “Computationally efficient
FDTD simulation of open-region scattering problems on the connection ma­
chine CM-5,” in IEEE Antennas and Propagation Society International Sym­
posium Digest, (Seattle, WA), June 19-24, 1994.
[5] B. Engquist and A. Majda, “Absorbing boundary conditions for the numerical
simulation of waves,” Math. Comput., vol. 31, pp. 629-651, 1977.
[6] G. Mur, “Absorbing boundary conditions for the finite-difference approxima­
tion of the time-domain electromagnetic field equations,” IEEE Trans. Electromagn. Compat., vol. EMC-23, pp. 377-382,1981.
[7] Z. P. Liao, H. L. Wong, B. P. Yang, and Y. F. Yuan, “A transmitting bound­
ary for transient wave analysis,” Scientia Sinica. (Series A), vol. 27, no. 10,
pp. 1063-1076, 1984.
[8] R. L. Higdon, “Numerical absorbing boundary conditions for the wave equa­
tion,” Math. Comput., vol. 49, pp. 65-90,1987.
[9] I. Katz, D. Parks, A. Wilson, M. Rotenberg, and J. Harren, “Non-reflective free
space boundary conditions for SGEMP codes,” Systems, Science and Software,
vol. SSS-R-76-2934, 1976.
[10] R. Holland and J. W. Williams, “Total-field versus scattered-field finitedifference codes: A comparative assessment,” IEEE Trans. Nuclear Sci.,
vol. NS-30, pp. 4583-4588, 1983.
[11] J.-P. Berenger in Actes du Collogue CEM, (Tregastel, France), 1983.
[12] C. Cerjan, D. Kosloff, R. Kosloff, and M. Reshef, “A nonreflecting boundary
condition for discrete acoustic and elastic wave equations,” Geophysics, vol. 50,
pp. 705-708, 1985.
[13] R. Kosloff and D. Kosloff, “Absorbing boundaries for wave propagation prob­
lems,” J. Comput. Phys., vol. 63, pp. 363-376,1986.
[14] C. M. Rappaport and L. Bahrmasel, “An absorbing boundary condition based
on anechoic absorber for EM scattering computation,” J. Electromag. Waves
Appl., vol. 6, no. 12, pp. 1621-1634, 1992.
[15] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic
waves,” J. Comput. Phys., Mar. 1994.
114
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 6
C O N C L U S IO N S A N D D IR E C T IO N S F O R F U T U R E W O R K
In our introduction, we proposed a list of nine criteria against which inverse
scattering algorithms should be judged. It was argued and shown in this thesis that
nonlinear iterative inverse scattering algorithms are capable of achieving very high
resolution images and inverting larger contrast objects than can any other general
method th at we know of. Hence, nonlinear inverse scattering algorithms satisfy
criteria (1) and (2) very well.
The nonlinear optimization technique th at we use in our algorithms is very gen­
eral and therefore may be applied to many different types of wave equations such as
the acoustic and elastic wave equations, in addition to electromagnetic wave equa­
tions. The approach is also valid for multidimensional wave equations and is not
limited to 1-D or 2-D problems. Arbitrary numbers and shapes of scatterers do not
pose any added difficulties in our formulation and a minimal amount of a priori in­
formation about object locations and shapes is assumed. The algorithms presented
here have been tested and verified using real experimental data with practical mea­
surement geometries. High-quality reconstructions were obtained in the presence
of both random (time-varying) systematic (not time-varying) experimental errors.
Criteria (5) - (9) are therefore satisfied as well.
The most difficult criteria for our nonlinear inverse scattering algorithms to
satisfy are (3) and (4), the maximum-size object that may be inverted and the
speed of the algorithms. Since a forward scattering problem must be solved at each
iteration of the inversion algorithm, our inverse scattering algorithms are limited
by the computational complexity of the forward scattering solution. This, in effect,
places a limit, on the maximum size object th at may be inverted on a given computer
in an acceptable amount of time. Fortunately, the computational complexity of
the FDTD algorithm, 0 ( N 1'5) for 2-D problems and 0(JV133) for 3-D problems,
is very reasonable. One disadvantage of using the FDTD algorithm is that the
115
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
forward scattering problem must be solved separately for each, transm itter location.
The implementation of the FDTD algorithm on the Connection Machine CM-5 has
allowed us to solve larger problems and has made it practical to solve 3-D nonlinear
inverse scattering problems. It is anticipated that future advances in computer
technology and parallel processing will allow us to solve even larger problems.
The distorted-Bom iterative method was presented for the 2-D TM and TE
polarization cases. In the TM polarization case, it was shown th at objects with
permittivity contrast as high as 10:1 could be reconstructed faithfully. Large ob­
jects several wavelengths in diameter also posed no problem for the reconstruction
algorithm.
The 2-D TE DBIM algorithm could invert low-contrast objects with permittivity
contrasts of 2:1, but had problems reconstructing the high-contrast case. In fact, the
algorithm did not converge to the correct object when an object with contrast 10:1
was considered. The reason is that the 2-D TE inverse scattering problem is really
a vector problem due to the buildup of polarization charges inside the scattering
object.
A local-shape-function inverse scattering algorithm was presented for imaging
very strong scatterers such as metallic scatterers. This algorithm parameterizes
the scattering object using a shape function variable that varies between 0 and 1,
and makes the inverse scattering problem more linear. The LSF algorithm was
shown to be capable of resolving more closely spaced metallic objects, and hence,
have a higher resolution capability, than the DBIM algorithm with conductivity
optimization. The LSF algorithm was also found to converge faster than the DBIM
algorithm for the metallic scatterer case.
Broadband time-domain data are generally preferred for inverse scattering be­
cause a much greater information content is contained in a broadband pulse than is
available from frequency-domain data at a few discrete frequencies. Time-domain
data also allow time gating to be used to eliminate unwanted early-time and late­
time arrival signals. The “phase wrapping” problem present in diffraction tomo116
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
graphy, whereby the phase of the wave becomes ambiguous for objects of several
wavelengths in diameter, is not present when transient data are used along with an
iterative algorithm such as DBIM or the LSF method.
While it is preferable to have time-domain data available for the inverse scat­
tering algorithms, it is usually easier to collect data in the frequency domain.
Frequency-domain data collection systems, such as step-frequency radar, offer high
signal-to-noise ratio due to the inherent narrowband electronics. Very stable signal
sources such as synthesized sweepers axe available for CW measurement systems,
whereas time-domain data collection systems are subject to timing jitter and zerolevel drift. Systematic errors inherent in microwave measurements may be removed
effectively in frequency-domain systems using various calibration techniques. Cali­
bration is much more difficult with time-domain systems. Finally, a large percentage
of the total source power is available at each measurement frequency in a CW sys­
tem. The source power in a time-domain system is distributed across the various
frequency components of the source pulse in time-domain systems.
A prototype step-frequency radar system targeted for nondestructive evalua­
tion applications was presented. Reconstructions of both metallic and dielectric
objects including metallic rods, plastic PVC pipes and a hollow glass cylinder were
shown from experimental data collected with the prototype system. A commercially
available monostatic impulse radar system used for ground-penetrating radar appli­
cations was also discussed. The commercial impulse radar system was found to be
not suitable for nonlinear inverse scattering because the transmitted source pulse
and the gain of the antenna were not known and therefore could not be modeled
properly. A rudimentary bistatic impulse radar was constructed using a picosecond
step generator, a pair of broadband Vivaldi antennas and a sampling oscilloscope.
Unfortunately, however, this project had to be abandoned because the sampling os­
cilloscope available in our laboratory could not resolve our source pulse, and funds
were not available to purchase a newer digital sampling oscilloscope.
Finally, a method was presented for implementing a very efficient FDTD code
on a massively parallel supercomputer such as the Connection Machine CM-5. The
117
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
main challenge in designing an efficient code was actually in the implementation
of the absorbing boundary condition at the edge of the FDTD grid. Conventional
absorbing boundary conditions break the nearest-neighbor structure of the FDTD
algorithm and therefore require a much greater amount of communication opera­
tions than the core FDTD operations. A recent breakthrough in absorbing material
boundary conditions using a perfectly matched medium allows for an extremely ac­
curate absorbing material boundary using a border consisting of a few grid points
around the computation domain. The new approach allows FDTD computations
to be performed fast and accurately on a massively parallel supercomputer. This
new FDTD code also allows us to solve much larger inverse scattering problems in
a shorter amount of time.
In this thesis, we have touched on many areas of inverse scattering including
inverse scattering theory and the development of new imaging algorithms, data
collection methods, and the efficient implementation of a forward scattering solver.
There are obviously many directions in which to extend this work in each of these
areas. In the development of new imaging algorithms, the present 2-D imaging
algorithms should be extended to the 3-D vector case. Most of the theoretical
work on inverse scattering algorithms up until now have been limited to the 2-D
case because we previously did not have the computational resources to solve 3-D
inverse problems. There are several salient features of the vector problem, such as
the buildup of polarization charge inside the scattering object, that were not present
in the 2-D scalar inverse scattering problems.
To generate the best possible image reconstructions using experimental data,
improved calibrations techniques have to be developed to remove many of the sys­
tematic errors that are present in inverse scattering measurements. One of the
major sources of systematic error is the frequency-dependent, position-dependent
response of the array. In this thesis, very simple calibration techniques were used,
such as measuring the gain of the antennas at boresight only. We should be able
to achieve higher resolution images from the experimental apparatus by taking the
frequency-dependent, position-dependent response of the array into account in a
calibration procedure.
118
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Up until now, we have considered very simple scattering objects to test our
inverse scattering algorithms and measurement systems. Even in the computer
simulation results, simple “test objects” were used to demonstrate the resolving
power of our algorithms for various sizes and contrasts of objects. If our a lg o r ith m s
and measurement apparatus are going to have a technological impact on the many
application areas, we must consider more realistic scatterers similar to the ones that
would be found in a practical situation. This would require focusing efforts towards
particular applications and determining the specific needs and requirements of those
applications.
119
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
V IT A
William Harold Weedon was bom on August 29, 1966 in Port Jefferson, NY,
and raised in the town of Warwick, RI, where he attended Toll Gate High School. In
the Fall of 1984, Mr. Weedon enrolled in the College of Engineering at Northeastern
University, Boston, MA. It was at Northeastern that he acquired a deep interest in
electromagnetics and wave propagation.
As a cooperative education student, Mr. Weedon worked for four years (19851989) at the now defunct Submarine Signal Division of Raytheon Company, in
Portsmouth, RI. At Raytheon, he had an opportunity to work on many interesting
research and development problems such as the MK-50 torpedo proposal and the
AN/SQQ-89I towed array sonar. He also participated in field performance evalua­
tion and at-sea trials of the AN/SQQ-32 minehunting sonar abord a Navy ship in
the Gulf of Mexico.
After graduating third in his class of the College of Engineering, with a Bachelor
of Science in Electrical Engineering, Mr. Weedon was persuaded by several senior
faculty members and the assistance of a very attractive fellowship from the Cen­
ter for Electromagnetics Research at Northeastern to remain for graduate study.
He remained at Northeastern for another year and obtained his Master of Science
degree in Electrical Engineering under the supervision of Professor Anthony J. Devaney. However, he later decided to pursue his Ph.D. degree at another university
in the interests of academic diversity. He spent the summer of 1990 at The MITRE
Corporation, Bedford, MA, where he worked on the phased array antenna tracking
problem.
In January of 1991, Mr. Weedon joined Professor Weng Cho Chew’s research
group at the University of Illinois at Urbana-Champaign, where he remained for his
Ph.D. While at Illinois, he had the rare opportunity to work on both theoretical
and experimental aspects of electromagnetics. He also had an opportunity to gain
some teaching experience by helping Professor Chew write lecture notes and teach
120
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
several lectures of a new undergraduate course in computational electromagnetics
in the Fall of 1993.
A list of Mr. Weedon’s technical publications is given below:
R efereed P apers:
1. W.H. Weedon, S.W. McKnight and A.J. Devaney, “Selection of optimal angles
for inversion of multiple-angle ellipsometry and reflectometry equations,” J.
Opt. Soc. Am. A., Vol. 8, p. 1881-91, Dec. 1991.
2. W. H. Weedon, W. C. Chew, J. H. Lin, A. Sezginer and V. Druskin, “A 2.5D
scalar Helmholtz wave solution employing the spectral Lanczos decomposition
method (SLDM),” Microwave and Optical Tech. Lett., Vol 6, pp. 587-92, Aug.
1993.
3. W. H. Weedon and W. C. Chew, “Time-domain inverse scattering using the
local shape function (LSF) method,” Inverse Probl., Vol. 9, pp. 551-64, Oct.
1993.
4. W. C. Chew, G. P. Otto, W. H. Weedon, J. H. Lin, C. C. Lu, Y. M. Wang and M.
Moghaddam, “Nonlinear diffraction tomography-the use of inverse scattering
for imaging,” Int. J. Imaging Sys. Technol., 1994. Submitted for publication.
5. W. C. Chew and W. H. Weedon, “A 3-D perfectly matched medium from mod­
ified Maxwell’s equations with stretched coordinates,” Microwave and Optical
Tech. Lett., 1994. Submitted for publication.
C on feren ce P a p ers and T echnical R eports:
1. C.A. Dimarzio, A.J. Devaney, R.E. Grojean, W. Crosby, W.H. Weedon, “Dis­
crimination between coherent and incoherent light,” Interim Report Prepared
for Army Natick R.D. and E. Center, (Natick, MA), June 1989.
2. W.H. Weedon, S.W. McKnight and A.J. Devaney, “Selection of optimal an­
gles for inversion of multiple-angle ellipsometry and reflectometry equations,”
Annual Meeting of the Optical Society of America, (Boston, MA), Nov. 4-9,
1990.
3. W.H. Weedon, “Phased array sensor allocation for highly maneuverable targets
using simplified utility analysis,” working paper, The MITRE Corp., (Bedford,
MA), June 1991.
121
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4. W.H. Weedon, J.E. Mast, W.C. Chew, H. Lee and J.P. Murtha, “Inversion of
real transient radar eata using the distorted-Bom iterative algorithm,” IEEE
Antennas and Propagation Society International Symposium Digest, (Chicago,
IL), July 18-25, 1992.
5. W. H. Weedon, W. C. Chew, J. H. Lin, A. Sezginer and V. Druskin, “A 2.5D
scalar Helmholtz wave solution employing the spectral Lanczos decomposition
method (SLDM),” IEEE Antennas and Propagation Society International Sym­
posium Digest, (Ann Arbor, ME), June 28-July 2, 1993.
6. W. H. Weedon and W. C. Chew, “A local shape function (LSF) method for
time-domain inverse scattering,” IEEE Antennas and Propagation Society In­
ternational Symposium Digest, (Ann Arbor, MI), June 28-July 2, 1993.
7. W. H. Weedon and W. C. Chew, “Broadband microwave inverse scattering for
nondestructive evaluation (NDE),” Proceedings of the 20th Annual Review of
Progress in Quantitative NDE, (Brunswick, ME), Aug. 25-Sept. 2, 1993.
8. W. C. Chew, G. P. Otto, J. H. Lin, W. H. Weedon and C. C. Lu, “Nonlinear
inverse scattering—past, present, and future” X XIVth General Assembly of the
International Union of Radio Science, (Kyoto, Japan), Aug. 1993.
9. W. C. Chew, G. P. Otto, W. H. Weedon, J. H. Lin, C. C. Lu, Y. M. Wang and
M. Moghaddam, “Nonlinear diffraction tomography-the use of inverse scatter­
ing for imaging,” Twenty-Seventh Asilomar Conference on Signals, Systems, &
Computers, (Pacific Grove, CA), Nov. 1-3, 1993.
10. W. H. Weedon, W. C. Chew and C. M. Rappaport, “Computationally effi­
cient FDTD simulation of open-region scattering problems on the Connection
Machine CM-5,” IEEE Antennas and Propagation Society International Sym­
posium Digest, (Seattle, WA), June 19-24, 1994. Accepted for presentation.
11. W. H. Weedon, W. C. Chew and C. Ruwe, “Step-frequency radar imaging for
NDE and GPR applications,” SPIE Advanced Microwave and Millimeter Wave
Detectors, (San Diego, CA), July 24-29,1994. Accepted for Presentation.
M iscella n eo u s C ontributions:
1. I Aksun, J. H. Lin, C. C. Lu, G. Otto, Y. M. Wang, R. Wagner and W. H.
Weedon, Solution Manual for Waves and Fields in Inhomogeneous Media, W.
C. Chew, Ed., 1993.
122
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2. W. C. Chew and W. H. Weedon, Introduction to Computational Methods for
Electric Fields, unpublished lecture notes for course ECE 371 in the Department
of Electrical and Computer Engineering at the University of Illinois, Urbana,
IL, Fall 1993.
Mr. Weedon in the Electromagnetics Laboratory at the University of
Illinois, Urbana-Champaign.
123
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Документ
Категория
Без категории
Просмотров
0
Размер файла
5 443 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа