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.

1/--страниц