LIBRARY

Michigan State
UnIversIty

 

 

 

PLACE IN RETURN BOX
to remove this checkout from your record.
TO AVOID FINES return on or before date due.

 

MTE DUE MTE DUE MTE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

me amounts-mu

TIME-DOMAIN IMAGING OF RADAR TARGETS USING
ULTRA-WIDEBAND OR SHORT PULSE RADARS

by

Yingcheng Dai

A DISSERTATION
Submitted to
Michigan State University

in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering

1997

ABSTRACT

TIME-DOMAIN IMAGING OF RADAR TARGETS USING
ULTRA-WIDEBAND OR SHORT PULSE RADARS

By

Yingcheng Dai

The development of viable short-pulse radar system has renew the interest in time
domain imaging performed directly in time-domain with temporally measured signal.
Since the short-pulse response of a target provides significant information about the
positions and strengths of scattering centers, and if observations are made over a wide
range of aspect angle, one might created an image of the target using the short-pulse
response information. In this thesis, we have developed and implemented a time-domain
radar imaging technique based on a space-time magnetic field integral equation, using a
sine modulated exponential pulse, and employing the inverse Radon transform. Images
of various aircraft models were created from measured target responses over a wide band
of frequencies and over the entire range of aspect angles.

For the limited-view problem, two techniques have been proposed to process this
practical situation. One of the approach is the method of projections onto convex sets

(POCS) which has been use in image processing for a long time. We extend this

approach to radar imaging for the first time and show some useful results. Another
approach which we have demonstrated is to process the available measured projections
in order to generate an estimate of the fill set of projections, an image which is called
a sinogram. The goal .of this approach is to recover the sinogram from the available
measured data using linear prediction. Since the scattered field of a target can be written
as a superposition of distinct specular reflection arising from scattering centers on the
target, the position and strength of the scattering centers can be predicted using linear
prediction with the change of the observation angle. Thus the missing data can be
predicted before reconstructing the image.

In the imaging of complex radar target, the PO approximation is used in the
reconstruction algorithm. However, the PO approximation is inadequate for scattering
problems of a complex shaped conducting object such as aircraft. At high frequency, edge
diffractions, multiple reflections, creeping waves, and surface travelling waves may also
be important scattering mechanisms. Additionally, the spectral and angular windows for
data are usually restricted by practical constraints. Therefore the time domain image of
a aircraft may be different from their geometrical shape. We have investigated time
domain imaging of aircraft employing SMEP responses, and interpret the reconstructed
image from a new approach, based on analysis of the scattering mechanisms and the
back-projection algorithm utilized in image retrieval. The time-domain inverse scattering
identity with the incorporation of Geometrical Theory of Diffraction (GTD) is derived and

some interesting experimental results are provided.

COPyright by
Yingcheng Dai

1997

ACKNOWLEDGMENTS

I would like to express my sincerest appreciation to Dr. K. M. Chen. His
encouragement to finish my graduate studies at Michigan State University and his
teaching and guidance during my research will not be forgotten. Thank him very much
for his generous support throughout my work. I am especially grateful to Dr. E. J.
Rothwell. His sincere interest, encouragement and expert technical advice will long be
remembered. To D. P. Nyquist, thank for his outstanding teaching, deep concern and
guidance. To Dr. Byron Drachman, thank him for his participating as a member of my
committee.

I also want to thank my fellows at EM lab. It has been a great pleasure to work
and learn with them.

Finally, I would like to express my special gratitude to my wonderful wife, Qing

Li, and my parents for their constant love, financial support, and understanding.

TABLE OF CONTENTS

LIST OF TABLES ............................................. ix
LIST OF FIGURES ............................................ x
CHAPTER 1
INTRODUCTION ....................................... 1
CHAPTER 2
INTEGRAL EQUATION AND ITS SOLUTION FOR
TRANSIENT SCATTERING ............................... 5
2. 1 Introduction ....................................... 5
2.2 Time-Domain Integral Equations ......................... 6
2.2.1 Magnetic Field Integral Equation .................... 6
2.2.2 Electric Field Integral Equation .................... 10
2.3 Numerical Solution of MFIE ........................... 10
2.3.1 Rigorous Solutions .............................. 11
2.3.2 Instabilities in Marching-on-in-Time Methods ........... 15
2.3.3 Approximations ................................ 18
CHAPTER 3
TIME-DOMAIN INVERSE SCATTERING IDENTITY ........... 20
3. 1 Introduction ....................................... 20
3.2 The Radon Transform ................................ 21
3.2.1 Radon Transform in Several Dimensions ............... 21
3.2.2 Two Dimensions .............................. 22
3.3 Derivation of Space-Time Integral Equations
and Inverse Scattering Identity .......................... 25
3.3.1 The Bistatic Case ............................... 25
3.3.2 The Monostatic Case ............................ 35
3.3.3 Rotational Symmetric Targets ..................... 37
3.3.4 Flying Radar Targets ........................... 39
3.3.5 Three Dimensional Case ......................... 42

vi

CHAPTER 4

TIME-DOMAIN IMAGING OF RADAR TARGETS
USING ALGORITHMS FOR RECONSTRUCTION

FROM PROJECTIONS ................................... 46
4. 1 Introduction ....................................... 46
4.2 Backprojection ..................................... 47
4.3 Reconstruction Algorithm .............................. 47
4.4 Numerical Results and Images for a Metal Sphere ............. 52
4.5 Experimental Results and Images for Aircraft ................ 53
CHAPTER 5
IMAGING OF RADAR TARGETS USING
LIMITED-VIEW DATA .................................. 79
5. 1 Introduction ....................................... 79
5.2 The Method of Projection Onto Convex Sets (POCS) .......... 80
5.2.1 Basic Ideas .................................. 80
5.2.2 Convex Projections in Short—Pulse Radar imaging ....... 81
5.2.3 Experimental Results ........................... 86
5.3 Radar Image Reconstruction Using Sinogram Restoration
and Linear Prediction ................................. 87
5.3.1 Sinograms ................................... 87
5.3.2 Linear Prediction .............................. 98
5.3.3 Experiment Results ............................. 102
5.4 Conclusions ...................................... l 04
CHAPTER 6
IMAGING UNDERSTANDING AND MODIFICATION USING
HIGH FREQUENCY APPROXIMATION METHODS ........... l 18
6.1 Introduction ...................................... l 1 8
6.2 Scattering Field Understanding Using Ray Method ........... 1 19
6.3 High Frequency Approximation Methods .................. 124
6.3.1 Equivalent Edge Currents for Arbitrary
Aspect of Observation .......................... 125
6.3.2 Scattered Field from Equivalent Fringe
Currents of a Finite Edge ........................ 129
6.4 Time-Domain Inverse Scattering Identity with The
Incorporation of Edge Diffraction ...................... 130
6.5 Images of an Object Consisting of Conducting Plates ........ 132
6.5.1 Edge Diffracted Field .......................... 134
6.5.2 Scattering Field Under Physical Optics Approximations . . 136
6.5.3 Images of The Plates .......................... 138
6.6 lrnages of an Aircraft Model .......................... 138
6.7 Conclusions ...................................... 1 39

vii

CHAPTER 7

CONCLUSIONS ....................................... 150

7.1 Summary of The Thesis. ............................... 150

7.2 Suggestions for Further Studies .......................... 152
APPENDIX A

ANECHOIC CHAMBER MEASUREMENTS .................. 154

Al Introduction ...................................... 154

A2 Calibration Procedure ............................... 156
APPENDIX B

EXPERIMENTAL EQUIPMENT ........................... 165
BIBLIOGRAPHY ........................................... l7]

viii

Table 8.1

Table B.2

Table B.3

Table 8.4

LIST OF TABLES

Equipment used for frequency domain scattering

measurements ................................ 1 66
Low band measurement parameters ................. 166
High band measurement parameters ................ 167
Enhanced high band measurement parameters .......... 167

ix

Figure 2.1

Figure 3.1
Figure 3.2
Figure 3.3
Figure 3.4
Figure 3.5
Figure 3.6
Figure 4.1

Figure 4.2

Figure 4.3

Figure 4.4

Figure 4.5

Figure 4.6

Figure 4.7

Figure 4.8

Figure 4.9

Figure 4.10

LIST OF FIGURES

Scattering of an incident electromagnetic wave by a perfect

conducting body. ..................................... 7
Geometry of 2-D Radon transform. ....................... 23
Geometry of the scatterer and graphic View of space parameters. . . 27
Geometry of 3-D bistatic body reconstruction problem. ......... 33
Geometry of 3-D monostatic body reconstruction problem. ...... 36
Geometry of rotationally symmetric scattering problem. ......... 38
Geometry of the flying object. .......................... 4O
Geometry for obtaining the backprojection ................... 48
The filter response for the filtered backprojection algorithm. It has

been bandlimited to 1/23.. ............................. 55
The impulse response of the filter shown in Figure 4.2. ......... 56
Comparison of theoretical (a=8><10’, fc=IOGHz in (3.23)) and
synthesized (4-16GHz) SMEP. .......................... 57
Synthesized SMEP spectrum. ........................... 58
Far-zone SMEP response of a 14 inch sphere using

marching-on-in-time method. ........................... 59
Result of using stabilized formula to the scattering field shown in

Fig. 4.6. .......................................... 6O
Far-zone SMEP scattered field of a 14 inch sphere. ............ 61
Cross-sectional area function of a sphere 14 inches in diameter. . . . 62

Upper surface image of the 14 inch diameter sphere using PO

X

Figure 4.11

Figure 4.12

Figure 4.13

Figure 4.14

Figure 4.15

Figure 4.16

Figure 4.17

Figure 4.18

Figure 4.19

Figure 4.20

Figure 4.21

Figure 4.22
Figure 4.23
Figure 4.24

Figure 4.25

Figure 5.1

Figure 5.2

approximation. ..................................... 63

Upper surface image of the 14 inch diameter sphere considering the
correction terms. .................................... 64

Image of TR-l from 0°-180° data in band 4-16 GHz using
monostatic identity. .................................. 65

Image of F-l4 from 0°-180° data in band 4-16 GHz using
monostatic identity. .................................. 66

Image of B-52 from 0°-180° data in band 4-16 GHz using
monostatic identity. .................................. 67

Image of TR-l from 0°-180° data in band 4-16 GHz using bistatic

identity. .......................................... 68
Image of F-l4 from 0°-180° data in band 4-16 GHz using bistatic
identity. .......................................... 69
Image of B—52 from O°-180° data in band 4-16 GHz using bistatic
identity. .......................................... 70
Image of B-52 from 0°-90° data in band 4-16 GHz using bistatic
identity. .......................................... 71
Comparison of theoretical (a=4X109, fi=IOGHz in (3.23)) and
synthesized (7-13 GHz) SMEP ........................... 72
The spectrum of the synthesized SMEP ..................... 73
Image of 8-52 from O°-180° data in band 7-13 GHz using bistatic
identity (before windowing). ............................ 74
8-52 SMEP response (nose-on) before windowing. ............ 75
Window produced based on Figure 4.22. ................... 76
8-52 SMEP response (nose-on) after windowing. ............. 77
Image of 8-52 from O°-180° data in band 7-13 GHz using bistatic
identity (after windowing). ............................. 78
A symbolic representation of the projection off onto C. ........ 82
Discrete object, F(x,y), and its projection are shown

for an angle cl). ..................................... 85

xi

Figure 5.3

Figure 5.4

Figure 5.5

Figure 5.6

Figure 5.7

Figure 5.8

Figure 5.9

Figure 5.10
Figure 5.11

Figure 5.12

Figure 5.13

Figure 5.14

Figure 5.15

Figure 5.16

Figure 5.17

Figure 5.18

Figure 5.19

Figure 5.20

Geometry of discrete test image to be reconstructed. ........... 88

Image of F-l4 from 0°-180° data in band 4-16 GHZ using POCS

after 1 iteration. .................................... 89
Image of F-14 from O°-180° data in band 4-16 GHz using POCS

after 5 iteration. .................................... 90
Image of F-l4 from 45°-135° data in band 4-16 GHz using POCS

after 5 iteration. .................................... 91
Image of F-l4 from 45°-135° data in band 4-16 GHz using
convolution backprojection algorithm. ..................... 92
Geometry of the 2-D Radon transform. .................... 94
Domain C over which the sinogram is measurable in the

limited-view problem. ................................ 95
A 1:48 scale model TR-l aircraft. ........................ 96
The sinogram of TR-l aircraft. 35 dB range, 11 colors. ......... 97
(a) Discrete all-pole model in the frequency domain.

(b) Discrete all-pole model in the time domain. .............. 99
Measurable part of the sinogram of TR-l aircraft. ............ 106
Measurable part of the sinogram of TR-l aircraft formed

by the positions of the scattering centers. .................. 107
Sinogram of TR-l formed by the positions of

the scattering centers ................................. 108
Sinogram of TR-l formed by the predicted

positions of the scattering centers. ....................... 109
Measured and predicted positions of the tail of

the TR-l aircraft. .................................. 110
Measured and predicted positions of the engines of

the TR-l aircraft. .................................. 111
Measured and predicted amplitude of the scattered field in the

positions of the tail of the TR-l ........................ 112

Measured and predicted amplitude of the scattered field in the

xii

Figure 5.21
Figure 5.22
Figure 5.23

Figure 5.24

Figure 6.1
Figure 6.2

Figure 6.3

Figure 6.4

Figure 6.5

Figure 6.6

Figure 6.7

Figure 6.8

Figure 6.9

Figure 6.10

Figure 6.11

Figure 6.12

Figure A.l

positions of the engines of the TR-l. .....................
Measured and predicted data of the TR-l for ¢=45°. ..........
Restored sinogram of the TR-l by linear prediction method.

Image of the TR-l using restored sinogram. ................

Image of the TR-l from 45°-135° data using convolution
backprojection algorithm . ............................

Scattering by an aircraft model (composite obstacle). ..........

Wedge scattering geometry. ...........................

Geometry of a conducting plate in the laboratory coordinate system
and edge-fixed coordinate system. .......................

Diffracted SMEP response of the object with rotational angle
¢=45°. ..........................................

Diffracted SMEP response of the object with rotational angle

4>=O°. ...........................................

SMEP response of the object with rotational angle ¢=0° under PO
approximation. ....................................

SMEP response of the object with rotational angle ¢=45° under PO
approximation. ....................................

Image of the object from 0°-180° data in band 4-16GHz using the
total scattered field. .................................

Image of the object from 0°-180° data in band 4-16GHz under PO
approximation. ....................................

Image of B-52 from 0°-180° data in band 4-16GHz using
measured field. ....................................

Image of main edges of 8-52 from 0°-180° data in band 4-16GHz
Using theoretical edge diffraction field. ...................

Image of B-52 from 0°-180° data in band 4-16 GHz with the
incorporation of edge diffraction. .......................

Dual antenna free—field frequency domain transient
measurement system. ................................

xiii

113

114

115

116

117

120

127

133

141

142

143

144

145

146

147

148

155

Figure A.2

Figure A.3

Figure A.4

Figure A.5

Figure A.6

Figure A.7

Figure 8.1

Figure 8.2

Figure B.3

Spectral response of 14 inch metal calibration sphere, measured
using frequency domain system in MSU anechoic chamber. .....

Spectral response of 14 inch metal calibration sphere after time-
domain filtering. Measured using frequency domain system in MSU
anechoic chamber. ..................................

Geometry for sphere scattering problems. ..................

Theoretical response of 14 inch metal calibration sphere calculated
using Mie series for MSU anechoic chamber configuration. .....

Transfer function of frequency domain system obtained using 14
inch diameter sphere as calibration standard. ...............

Comparison of theoretical and experimental frequency responses of
a 3 inch diameter sphere. Measurements performed using frequency
domain system in MSU anechoic chamber ..................

Gain of H-8349B Microwave Amplifier measured using HP-8720B
network analyzer. Amplifier input power=-30 dBm ...........

Gain of PPL-5812 wideband amplifier measured using HP-87ZOB
network analyzer. Amplifier input power=-10 dBm ...........

Measured cable attenuation ............................

xiv

159

160

161

162

164

168

169

170

CHAPTER 1

INTRODUCTION

Microwave imaging of airborne targets has received considerable interest in recent
years. Most attention has focused on the use of inverse synthetic aperture radar (ISAR)
[50,51]. In ISAR, the target is modelled as a collection of scattering centers at sufficiently
high frequencies and the image is constructed from a 2-D inverse Fourier transform of
the Cartesian frequency spectra, or a backprojection algorithm along with a 1-D inverse
Fourier transform of the polar frequency spectra [50]. The development of viable short-
pulse radar systems has renewed the interest in time-domain imaging performed directly
in the time domain with temporally measured signals. The use of time- domain techniques
was first discussed by Kennaugh and Cosgroff in 1957 [6]. Since then, many researchers
have developed approaches to the inverse scattering problem [7-12]. They have shown
that the target impulse, step, and ramp responses are related to the target geometry based
on physical optics principles. Under the physical optics approximation, Bojarski has
established a Fourier transform relationship between the geometry of the conducting
scatterer and a form of the scattering cross section [14]. Since this approach is based on
the physical optics approximation, it is valid only in the limit of high frequency. When
the size of the scatterer is comparable to the incident pulsewidth, the physical optics
solution is inadequate for this scattering problem. Furthermore, the impulse, step, or ramp
response of the target is hard to obtain in practice.

1

The purpose of this thesis is to develop and implement an imaging technique
directly in the time domain. Since the short-pulse response of a target provides significant
information about the positions and strengths of scattering centers, and if observations are
made over a wide range of aspect angles, one might create an image of the target using
the short-pulse response information. The starting point of this thesis is the time-domain
magnetic field integral equation (MFIE). The derivation of the MFIE and its numerical
solutions are provided in chapter 2. Chapter 3 gives the time-domain inverse scattering
identity (thickness function) developed from the time-domain MFIE and the Radon
transform. It shows that the size and shape of an object can be obtained from its cross-
sectional area function, and that the problem can be reduced to the classical Radon
problem. In this process, the Sine-Modulated Exponential Pulse (SMEP) is used as the
incident field.

Since the area function can be estimated from the scattered field when the incident
field is a SMEP, the remaining problem is thus essentially that of reconstruction from
projections. The basic idea of reconstruction of an image from a series of its projections
appears to have been first discussed by Radon [16]. The techniques that exist for
reconstruction fall into two directions. The whole operation can be done in frequency
space directly, or the equivalent of these expressions can be transformed into the spatial
domain. Whether implemented in the spatial domain or in the frequency domain, the
reconstruction algorithms can be conveniently interpreted by means of a straightforward
and interesting theorem, which is the projection-slice theorem. This theorem states that
the Fourier transform of a projection is a center cross section of the Fourier transform of
the projected object. Most modern tomographic systems are based on this theorem. Thus,
by formulating the object reconstruction problem as a form of the Radon problem, the

2

many reconstruction techniques and the properties of reconstruction from projections
which are well developed in other fields for the two-dimensional case can be extended
and applied to the inverse scattering problem. Chapter 4 gives the development of the
reconstruction algorithm and provides clear images with highly-defined edges using
stepped-frequency, multi-aspect measurement data of aircraft models.

For complete reconstruction, the time-domain approach requires data for all aspect
directions. However, in practice, information usually available only for limited view-
angles; thus, it is necessary to analyze the effect of incomplete information which may
cause the inverse problem to be ill-posed in nature. The limited-view problem has
attracted considerable attention in computed tomography imaging research, and techniques
for dealing with it have been proposed. Techniques can be put into two categories:
transform techniques that incorporate no a priori information, and finite series expansion
methods that may incorporate a priori information as constraints. The transform
techniques are usually single-pass direct reconstruction [18-20] while the finite series
expansion methods are usually iterative [21-27]. Projection onto convex sets (POCS) in
particular has shown great flexibility in dealing with known geometric constraints and
with noise. We will extend POCS into radar target imaging. The details are carried out
in chapter 5. Also in chapter 5 we will demonstrate a new reconstruction algorithm for
radar imaging in the limited-angle case. The goal of this approach is to recover the
sinogram from the available measured data using linear prediction. Since the scattered
field of a target can be written as a superposition of distinct specular reflections arising
from scattering centers on the target, the trace of the scattering centers can be predicted
using linear prediction with the change of the observation angle. Thus, the missing data
may be predicted before reconstructing the image.

3

In the proposed time domain imaging technique, the physical optics approximation
is used to model the scattered field of an aircraft. However, the PO approximation is
inadequate for scattering problems of a complex shaped conducting object such as aircraft.
At high frequency, edge diffractions, multiple reflections, creeping waves, and surface
travelling waves may also be important scattering mechanisms. Additionally, the spectral
and angular windows for data are usually restricted by practical constraints. Therefore,
the time domain image of an aircraft may be different from its geometrical shape. In
chapter 6 we will investigate time domain imaging of aircraft employing SMEP responses,
and interpret the reconstructed image from a new approach based on analysis of the
scattering mechanisms and the back-projection algorithm utilized in image retrieval. The
time-domain inverse scattering identity with the incorporation of Geometrical Theory of

Diffraction (GTD) is derived and some interesting experimental results are provided.

CHAPTER 2
INTEGRAL EQUATION AND ITS SOLUTION

FOR TRANSIENT SCATTERING

2.1 Introduction

Before the advent of the high-speed digital computers, most electromagnetic
scattering problems were solved in the frequency domain. The boundary value solutions
associated with electromagnetic radiation and scattering phenomena were based upon
analytical techniques that attempted to generate closed-form solutions and were, therefore,
limited to a few classes of problems. Approximate techniques such as geometric optics
(GO), physical optics (PO), and the Rayleigh approximation were applied to other
problems to obtain estimates of the scattered response. With the advent of high-speed
computers, the frequency-domain integral equation technique has also provided many
interesting and useful results.

The integral equation and its solution could also be obtained directly in the time
domain [1,2]. Indeed, the frequency-domain integral equations have their time-domain
counterparts [3]. With the developments of high resolution radar and electromagnetic
pulse (EMP), the investigation of radiation and scattering of transient waveforms from
conducting objects has attracted increasing interest. There are two integral-equation
techniques available for solving transient electromagnetic problems [3]. One of them
involves the Fourier transform of the frequency response of the scatter to synthesize the

5

time-domain response. The alternative technique is based on the time-domain integral
equation and its subsequent solution using the marching-on-in—time method. The details
of this approach and other approximate techniques will described in Section 2.3. The
derivation of the general integral equation for the scattering problem is covered in Section

2.2.

2.2 Time-Domain Integral Equations

The integral equation method has several advantages, such as geometrical
generality, compactness, and computability, which make it an attractive approach to
solving electromagnetic boundary value problems. In this section we will provide in
detail the derivations of the time-domain integral equations directly, rather than via the

inverse Fourier transform of the frequency-domain solution [2].

2.2.] Magnetic Field Integral Equation
Consider the scattering of electromagnetic waves by good conductors, as shown
in Figure 2.1. The electric and magnetic fields of the scatted wave satisfy the time-domain

Maxwell equations in the source-free, free space region outside the conductor.

with» - mogfi‘mt)

a,“ _ a~.-

va (a) so 6:5 (at) (2.1)
V~E‘(F,r) = o

V-fi‘(r',r) = o

where so and [.10 are the permittivity and permeability of free space, respectively, and F

incident
wave

 

Figure 2.1 Scattering of an incident electromagnetic wave by a
perfect conducting body.

is the position vector.

These fields can be derived from a scaler potential (11‘ and a vector potentialzfs
such that

Err,» = yum-gr? ‘01:), firm = iVxX‘rm (2.2)
1‘0

These potentials are related to each other by the Lorentz condition
V-A (r,t)+———¢ (r,t) = O (2.3)
c2 at

where c is the free-space light speed. Applying equations (2.1) and (2.2) yield the wave

equations satisfied by (1)5 and X5

[Va—ii]¢‘(f’,t) = o, [Va-iik‘tm = o (2.4)
c2 6:2 02 6:2
The retarded solutions of the wave equations (2.4) are

lF-F/llc)ds,/
[I

 

s —o 1 dl’t-
<1> (at) = f 9"

41teo s IF-i"

(2.5)

"3 '° -'l _ ~_ ~/
A (fit) = Eli/‘10,! Ir r llc)dS/
4n 8 IF-F’l
where .7 is the current density and Q is the charge density on the surface of the body.
They are related by the continuity equation

VJ(F,t)+-§Q(F,t) = o (2-6)

Substituting (2.5) into (2.2) gives the following integral representation of the scattered

magnetic field 173011) at a point F exterior to the scatterer and at time t

“5.. 1 “.4 [1 l " 6".4 /
H r,t = — Jr,t xV — -—R —Jr,r (2-7)
() 4,,f.( ) (R1CR2 "aH )st
where

R‘ = F—F’ R = |F-F’| (2-3)

1: = r-llr-F’l (2.9)
C

The operator V’ acts on the source point F’ . Equation (2.7) expresses the magnetic

field of the scattered wave at a point in space and time as a composition of excitations
from individual points of the scatterer’s surface at specific retarded times.

The derivation of the Magnetic Field Integral Equation (MFIE) begins by applying
the appropriate boundary condition. On the surface of a perfectly conducting scatterer, the
total magnetic field is related to the surface current density at all times through the

equation

Jim = fixrfi‘rat)+fi‘(at)1 (“01

where If is the outgoing unit vector normal to S.

Applying (2.10) to the representation (2.7) yields the time-dependent MF IE [2]
_, _. . " _. T -/ "
J(F,t) = ZfixH‘(f’,t)+_§_x f [1.9. J racyflhfids’ (2.11)
2n 8 c or R R2

where the integral is taken to be a principal value (P.V.) type. As we will see later, the
special form of this integral equation plays a fimdamental role in enabling us to construct
its solution using a marching-on-in-time method.

9

2.2.2 Electric Field Integral Equation
The time-domain scattering of an electromagnetic wave by a perfect conductor can

also be formulated in terms of an Electric Field Integral Equation (EF IE). Using (2.5) and

(2.2), one can obtain the following integral representation of the scattered electric field E s(F,t)
at a point I" exterior to the scatterer and at time t, in terms of the retarded values of the
surface charge density Q and the surface current density .7 induced on the scatterer’s

surface S

_1_ 1

E (r, t-) — —{— 1 VRQ(r’1:,)+—l ‘iW ,r)VR+uO—J(r’1:)}ds’ (2.12)

_4TI 3 R 80R 80 at
The EFIE for the scattering problem involving a perfect conductor may be derived

from the following boundary condition which holds for all points on S

If x[E"(r,r) + mm] = o (2.13)

Applying condition (2.13) one arrives at the following EFIE involving Q and .7

where the integral involved is also of the RV. type. Clearly, the EFIE is a Fredholm
integral equation of the first kind while the MFIE is a Fredholm integral equation of the

second kind.

2.3 Numerical Solution of MFIE
There are two useful integral equation techniques available for solving transient

electromagnetic problems. One of them involves using the Fourier transform of the

10

frequency response of the scatter to synthesize the time-domain response [3]. The
alternative technique is based on the time-domain integral equation and its subsequent
solution using the marching-on-in-time method. Here, we will describe the marching-on-

in-time method in detail.

2.3.1 Rigorous Solutions
A commonly used technique for solving transient scattering problems is the

marching-on-in-time method. Consider the MFIE

 

X—

- _, g
1.0“ if) R is, (2.15)
R R2

"_. . “1'... II 16"../
Jr,t =2n H r,t -— ——Jr,r +
,() x (12nxf{cat ,( )

where f denotes the principal value integral, and 1: =t-R/c is the retarded time.

This equation is clearly a Fredholm integral equation of the second kind. Note that

for a given time t, the unknown .75 inside the integral has an argument ‘I.’ =t -R/c. Since

the principal value integral excludes the point R=0, r is always less than t. Thus, we can

regard (2.15) as an expression for the unknown current .750), at any given time t, in terms

of the known incident term 2rixii‘, and an integral which is also known from the past
history of the current .73. This procedure forms the basis of an iterative technique for

constructing the solution of (2.15). The solution is started at the time t=to, invoking

causality which guarantees that all past currents and their derivatives are identically zero

for t<t0. The solution is then constructed by gradually building it up for t>to, by

marching on in time using a time-stepping procedure.

11

The method of moments can be applied to the solution of equation (2.15). First,

let’s express the unknown J; in an expansion of pulse basis functions

N
2,; tr,1V,.(F6U,.(r1 (“61

EM?-

_,, { 1 for ‘r" on the surface segment center at f1
Vi(r ) _ 0 elsewhere

where
for t in the time interval centered at tj

U j“) = {o elsewhere

The unknown 11".}. are the vector weight coefficients for the space-time sampled values of

the current .75 on the surface of the structure which is assumed to be subdivided into N5
patches centered at r,. The total time is also assumed to be subdivided into N, time-steps.

Since there is (175/ a: in equation (2.15), .730) and (175/ a: should be differentiable

at the observation time t, so we use a second-order Lagrange interpolation to approximate
the current

(t- tj‘)(t 5.1) -. (t_tj-1)(t'tj.1)—.
(ti-1%)“)- 1-111) ijl+(tj—tj-1)(tj-tj+l) ij+1

 

 

m .r1= A,(t1=

I}

(t—tj_l)(t-tj) _.
_ _ ij+1

 

where A” is the sample value of J at the center point of patch As at timer .,Then

substituting equation (2.16) into (2.15), and point matching at the center point of each

patch Asi gives

12

N! NI

Js(r t,) =2fixH (r. ,t)--§— [2251—quale—(t)+qu(r)—El-3xR’Vp(F’)Uq(t)ds’
2" P=lq=1 (2.18)

Note that p=i gives the contribution of the principal value integral due to the induced
currents in the self-patch itself. However, it can be shown that this contribution is zero

if Asi is a flat surface
jsxR = riJsR - fix(jst) = fixfiJsR = o (2.19)
We can consider .75 and R to be constant within patch Asp, and thus the integral

over each patch can be approximated by multiplying the current value at point Fp with

the area Asp. Then

xkipuqm i=1,2,...Ns

 

NS NI
~ . ~ m. -1
Js(ri,t) 2an (rnt) an 2 ZASr[atqu(r)c—l—ZCRW +Am(r)R1;

 

 

ptlpai q=1
(2.20)
. . . . R.

where Rip= If} —Fp]. Pornt matching at time t, and letting q =j-im{—A'!i— +0.5 , where

C t

int[---] denotes taking the integer part of the data, equation (2.20) becomes

- . ~. .. 11 GA __(_a.,t ) 1 - ~ .
Js(r,,r,) = A” = 2n xH (w!) - z—xzj + —3qu(1:) xRip 1 =1,...,1v,
It P=1priAspchp r-rj Rip ""1

(2.21)

The sample points in time are not independent of the space sample points used.

In order to let the two adjacent current sample in space not fall in the same time interval,

13

the time sample spacing At must be related to the space sample spacing AR by

cAt 5 AR (2.22)

Finally, we can use a method of marching-on-in-time to solve the problem. First,
assuming the start time is to, and all the surface currents on S are zero for all time less
than to to view the integral equation as an initial value problem. For instance, let us
assume that at time tl the incident field has just reached the scatterer. By virtue of the

retardation and the principal value nature of the integral, a current is induced on part of

the scatterer which is equal to 2rixHi(f'l,tl). As the time progresses further to 12. the

current at each point is then given by the known incident field 2fixHi(F,t2) plus a

contribution from current at other points on the scatterer at earlier times which is also

known. Using this idea gives

for j=1, A"). = A}, = 2rixH'(r;.,rl)

l (2.23)
+ 71W(t)} XROIPASP

 

 

N
. -' - .. ~r_. it s 1 6 ~
for 1:2, Aij = A12 = 2an (ri’t2)+2rr x 2 .{ngmhg
p=1 rm CR,-p

”‘2

Thus, marching on in time will allow us to build up the current solution from
previously known values. This method has many advantages, such as case of
programming and speed. Also, the solution is usually explicit and does not involve the
inversion of large matrices. However, a major disadvantage of the method is that the
computed solution often becomes unstable as time progresses, usually taking the form of

a very regular exponentially increasing oscillation which alternates in sign at each time—

14

step [4,5]. The cause of this instability will be seen to be related to the existence of
resonant frequencies at which the corresponding frequency domain integral equation has

more than one solution.

2.3.2 Instabilities In Marching-on-in-Time Methods

The occurrence of exponentially increasing instabilities is a common feature of
time marching method for solving transient scattering problems. In this section we will
examine how the instability arises and describe a simple method to eliminate the
instability. More details can be found in Rynne [4] and Smith [5].

The system equation (2.21) can be written as

N, 1v,I

Eu : firxgx‘f’iix Z “MN-q (2-24)
p=1 q=1

NI:—

 

where i = l,2---NS, j = 1,2---N,, a]; are some constants and N q=im{ AP +0.5]
C I

Now, recast this system of coupled linear equations as follows. Introduce the auxiliary

variables

“(4) _ " _ _
AN _ ApJ-q p ‘1’m’Ns q—l,---,Nq (2-25)

These are vectors in R3. Fix 1' and j. Then equation (2.24) shows that each component of

the 3-vector AS) is a linear combination of the components of the N,N,,-length vector/Ix.)

and the component of the forcing term fiixHiJ, and, moreover,

15

“(11 _ ~10) ~(21 _ ~(11 4M) _ 4M4) 2.26
Au“ “Au-1’ Au ' iJ-l "'3 Aw 'AiJ-I ( )

In short, each component of any such vector A3) is a linear combination of the

(P)
J

components of the NJVq-length vector A; and the component riixHiJ.

Now, let YJ- denote the following column vector in R3N‘N' assembled from the

auxiliary variables Ag.) (q=1 ..... M and q=1,...,Nq)

-0(N‘)

(0) “’(1)
U, Au,...,AU

”(Nf)

“ "(01
1;. = [A .,AN J,...,AN .1 17 (2.27)

Also, let hj denote the following column vector in R3N’N‘ assembled from the forcing

ICI'II'IS

h. = [ii xi? 0, 0, ri xH , 0, 0,..., ii x ” ., 0, 0]T (2.28)
l 2 2.1 N N

if

y], = 3134.11. (2.29)

Here B is a suitable matrix in terms of the constants apiq.

The initial quiescence of the system before the transient arrives requires that Y0=O.
If Hinc represents a transient of finite duration then H0=O, and Hj=0, for all j exceeding

some fixed M. Hence hj=0 for j>M. The solution is therefore (when j >M)

Y1. = BV'm(hM+BhM-1+...+BM‘1h1) (230)
Thus, if B has an eigenvalue of modulus exceeding 1, the vector Yj will grow

16

exponentially in modulus. To demonstrate this, suppose that B has diagonal Jordan form

A = diag(1,, 1,, 1%) (2.31)

where B=U ’ AU for some non-singular matrix U. Denoting UYJ- by Zj, the system (2.29)
is equivalent to

Z. = AZ. +Uh. (2.32)

1 1-1 J

The homogeneous solution to (2.32) are of the form

2]. = [a11{,a,A’,...,aNqurNflq]T (2.33)

where a ,,a2,... are constants, and any nonzero component is unstable if its corresponding

eigenvalue 7», is of modulus exceeding 1.

There is a simple remedy to the problem of the onset of instability. Consider the

following weighted average across 3 times steps:

1
yj = 2i(32+2B+1)Yj_l + (21443)},j (2.35)

which represents the effect of doing two time marching steps followed by an averaging
process (2.34). The eigenvalues of (BZ+ZB+l)/4 are (2+1)2, where A is any eigenvalue of
B. This has the virtue of shifting the troublesome eigenvalue near -1 to around zero, and
more generally of shifting any eigenvalue which may be slightly outside the unit circle
to inside. So this process will almost shift all the eigenvalues of the unstable marching

process inside the unit circle. Sample numerical results will be provided in Chapter 4.

17

2.3.3 Approximations

For scatterers that are large compared to the wavelength, the physical optics
approximation is a well-known method in the frequency domain. Similar approximations
can also be obtained in the time domain. In this section, we will extend the idea of the
physical-optics approximation in the frequency domain into the time domain.

For a perfectly conducting body, the surface—current-density distribution derived

in the frequency domain under the physical optics approximation is

2rixfi‘(?,e) i5 ~r‘i<0 (2.36)

.7 F,tr1 = _,.
po( ) {O k'°ii>0

Where If is the unit vector normal to the surface at point F, H i(F411) is the incident

magnetic field with angular frequency to , and It.i is the incident wave vector. For a time-

function excitation, the incident fields may be decomposed into their frequency-domain

components by superposition as follows

I? ‘01:) = zAmfitEwp = 21mm...) (2.371

and

7*‘1
a}.
A
O

X 2fixfi‘~(r;um) = 21ixfi“(r,r) (2 38)

7r:
:1.
v
O

0

Under the far-zone approximations, equation (2.7) can be simplified as

18

Ii" (r,t)
p0 cR

.. _. g, -
. _1_f j(;/,,)x£. am ’11,. R 1 (2.391
47! R3 61’ CR2

1 "°_./ / l _ 1 "‘ a -.-o/ /
EfJJo ,t)xV(—R) -—2RXEJ(r @er

 

A ”-‘l
r: -Lkamdsl
47tRC 5 6‘1."

Thus I?“ can be explicitly written in terms of the incident field as

 

—o" _./
Esau) = -2 IRCRxfri’xfi—gflds’ (2.40)
1t 3

where f denotes an integral over the illuminated portion of the target.
I

The physical optics approximation allows the direct evaluation of .7 from a

knowledge of the scatterer’s geometry and orientation with respect to the incident wave.
One of the short-comings of the PO approximations is its inability to accurately to

represent the current near the shadow boundary.

19

CHAPTER 3

TIME-DOMAIN INVERSE SCATTERING IDENTITY

3.1 Introduction

Target imaging and identification using electromagnetic responses in the time and
frequency domains has attracted increasing interest, with most methods carried out in the
frequency domain. The use of time domain techniques was first discussed by Kennaugh
and Cosgroff in 1957 [6]. Since then many researchers have developed approaches to the
inverse scattering problem [7-12]. They have shown that the target impulse, step, and
ramp responses are related to the target geometry based on physical optics principles.
Under the physical optics approximation Bojarski has established a Fourier transform
relationship between the geometry of the conducting scatterer and a form of the scattering
cross section [14]. Since this approach is based on the physical optics approximation, it
is valid only in the limit of high frequency. When the size of the scatterer is comparable
to the incident pulsewidth, the physical optics solution is inadequate for this scattering
problem. Furthermore, the impulse, step, or ramp response of the target is hard to obtain
in practice.

In this chapter we start from the space-time magnetic field integral equation, and
by using a Sine-Modulated Exponential Pulse (SMEP) waveform as the incident field, an
exact two-dimensional time-domain bistatic inverse scattering identity can be obtained
based on the inverse Radon transform. The formal definition of the Radon transform is

20

given in Section 3.2. In Section 3.3, the time-domain bistatic inverse scattering identity
is derived in details. Special cases such as the monostatic case, rotationally symmetric
target, flying object and three dimensional case are treated in Section 3.3.2, 3.3.3, 3.3.4,

and 3.3.5, respectively.

3.2 The Radon Transform

Radon transform theory has become a very important mathematical operation and
its applications are well known. They include computerized tomography (CT) applications
in, for example, diagnostic medicine, radio astronomy, electron microscopy, optical
interferometry, and geophysical exploration. These well-established concepts can be

extended and applied to the radar inverse scattering problem.

3.2.1 Radon Transform in Several Dimensions
The Radon transform of a function at a given hyperplane is defined as the integral
of the function over that hyperplane [1 1-13]. For a hyperplane in n-dimensional Euclidean

space defined by

(3.1)

m)
>11
II

"a

where f is the spatial position vector, E is a unit vector orthogonal to the hyperplane,

and p is the Euclidean distance from the origin, the Radon transform F (8 , p) of a function f(5r‘)

over the hyperplane is given by

21

F(E,p) = [EM/(soars (3.2)

The above equation may be expressed more conveniently in the following form

by using the Dirac delta function 8
F(€.p) = [11121 6(p—E-F)ds (3.3)
The inversion of the Radon transform consists of expressing f0?) in terms of its

integrals F(E,p) over the hyperplanes. The inversion formula for odd dimension in n is

_ 1 (n-l) - ~, ~
fli) - mm!" (5.5 ad: (3.41

where Fpm'l) is the (n-I)th derivative of F with respect to p. For even 11, the inversion

formula is

 

‘ a F(n-1) ,“
it) = 1 f dEf dp_£__:(_’LE_) (3.5)

(2a))" '5'“ p-E-fc’
Since the two-dimensional Radon transform has been used in this work, we will
mainly discuss this case as an example which will also help understanding the Radon

transform for higher dimensional cases.

3.2.2 Two Dimensions
We will use the coordinate system defined in Fig. 3.1 to describe line integrals.
Let (x, y) designate coordinates of points in the plane, and consider some arbitrary function

f defined on some domain D of 2-dimensional Euclidean space. If L is any line in the

22

f(X,Y) D

 

Figure 3.1 Geometry of 2-D Radon transform.

23

plane, then the line integral of f along all possible lines L is the two-dimensional Radon

transform off provided the integral exists. It can be written as

F = fLflx.y)ds (3'6)

where ds is an increment of length along L.
The equation of line L in Fig. 3.1 is

xcosrb +ysin¢ = p (3-7)

Then the transform can be written as an integral over two dimensional Euclidean space

by using the delta function to select the line L fi'om the space

F(p,d>) = f_:f_:fix,y)6(xc08¢ +ysind>-p)dxdy = ffiflbw-E-fldf (33)

where F = (x,y) and E = (c0541, sin41).

If F(p,¢) is known for all p and 1]), then F(p,<[)) is the two-dimensional Radon
transform of f(x,y). A projection is formed by combining a set of line integrals. The
simplest projection is a collection of parallel integrals as is given by F (p.¢) for a fixed
angle ([1.

Since E'F = xcosd) + ysincb, the inverse Radon transform for the two-dimensional

case (12:2) can be written from (3.5) as

a A
= - 1 II N w (309)
fix,» 7,21. d¢f__———p_£.f dp

If the replacements x = recs 6, y = rsinB are made, then f in polar coordinates becomes

24

a .
—F(p’5) (3.10)

1 r «- 3p
,0 = -— d
fir ) 4n2f° d¢f-~r?-rcos(¢-6) p

 

These formulas for n=2 constitute a solution to the problem of reconstruction from

projections.

3.3 Derivation of Space-Time Integral Equations and Inverse Scattering Identity

When an electromagnetic field is incident upon an object, currents and charges are
induced on and in the object. The induced currents and charges will then maintain a
scattered electromagnetic field. Once the induced currents flowing on the conducting
scatterer surface and the scatterer geometry are given, the scattered field can be calculated
directly. The expression for the surface currents (MFIE) has been derived in chapter 2.
In this Section, we will apply this expression to derive the inverse scattering identities for
the monostatic and bistatic cases, respectively. Finally, we will also discuss the cases of

the rotationally symmetric object, flying radar targets, and the 3-D case.

3.3.1 The Bistatic Case

From chapter 2 we know that the expression for the induced surface current .7 at
the point F on the scatterer surface and at the time t can be written as

. -i .. ~ - - .
7(f‘,t) = 211x}: (r,t) +J,(r,t) If°q>0 (3.11)
Jctfit) nfi<0

where

25

_1_+1_<‘9_ ~01“)de d, (3.12)

R2 Rat

 

where Wat) is the incident magnetic field, ti is the unit vector normal to the scatterer

surface, F is the position vector to the observation point, F’ is the position vector to the
integration point, R = lF-F’ l, 0“, = (F -F’)/R, (j is the unit vector to the transmitting

antenna, t=t-R, and t denotes normalized time in meters (ct).

From equation (3.11), we can see that the first term of the right—hand side
represents the direct influence of the incident field on the current at the observation point.
When applied to the illuminated side of the scatterer, it yields the familiar physical Optics
approximation for the surface current. The second term on the right-hand side of (3.11)
represents the influence of currents at other surface points on the current at the
observation point.

Once the surface current density has been obtained, the far scattered field of the

scatterer shown in Figure 3.2 can be calculated by the expression

vxf 3:25—st ’J (3.13)

 

where R5 = IFS-F’l, and r’ = r-R

Assuming that the observation point F3 is not on S, the curl operator may be carried

inside the integral. Using the far-field approximations and standard vector identities, we

have

26

 

Z tr
c

11115335315“

 

shadow region

 

 

 

receiving
d1rectron

Figure 3.2 Geometry of the scatterer and graphic view of space parameters.

27

(3.14)

11
g—a
h
h D
X
g.
‘i
I it
Q
a:

5‘01,»

where RS = (FS—FfillFs-F’l

In the far zone, we can use the approximations Rs = rs - fs-F’ for the retarded time

r’ and RS .~. rs and Rs z is for the amplitude term. Then (3.14) becomes

_ * _./
h‘(F,,t1 = -—1 r,x——a’(""’as' (3.15)
4713's 3 61',

Now, we define the aperture function

1 fi'fi(F/)>0 (3.16)

A =
0' n) o fi'r‘i(F’)<0

which is unity on the illuminated region of scatterer surface, and zero in the shadow
region. Substituting (3.11) into (3.15), we have

6.7 F’,r’
c( )d5
61/

_. ”i _./
1250;,» = -—1— rxri’xiég—“flAagriytsh—Lf f x I (3.17)
T s

S 3
2113's 8 41trs

II

_ _ A A O—o/ I = - z _ A .-¢/ ~
I ri rs+(r‘.+rs) r , and t t RS t rs+rs r . Usrng the vector

where 1: = t-Ri-Rs

identity fsx(fi’x}i ‘) ri’(fs°li ‘1—13‘0315’), (3.17) can be written as

_, i~/ 6.7 «x, .-
h ‘(F,.t) = J—fsri'-K§”_g:il.4(r,r)¢il- lef f FALL, 112' (3.13)

8
2m; 1117's 5 61/

where I? = is — (fs'li i)li i, and h ‘(Fs,t) = If i711073;). Now, letting the same incident pulse

28

illuminate the shadow region of the scatterer, we have

 

 

 

-° ‘ "’ 6.7 F’,t ..
h23(f°s,t) = Lffil'KmA2(fi,fi)dY/-——l-—ffsx___c_z_(___/)dsl .h‘ (3.19)
21173 S at 41trs s 61/
Combining (3.18) and (3.19) yields
_, 1' _.,
TITffi/xah ((12,?)‘13/ = H ‘(F.,t) (3.20)
rt 5 s
where
6.7 F’,'(: ..
11502,!) = Lff‘s)(._c€___2dg’14,'
41trs s 31/
6.7 F’,r .- (3-21)
err.» = —‘ fax—4 “1.1
41v: s 31/

 

H ’(FSJ) = h 3(FSJ) +h25(f°,,t) +h:1(7_,,t) +hcsz(F_,,t)

Applying Gauss’ law to the left hand side of (3.20), we have

f Wary = TIE—H855,” (3.22)
V ar" (K-flcosw/Z)

where fl. +r‘s = 2cos—g-r‘, and B is the bistatic angle as shown in Fig. 3.2. Assuming the

incident magnetic field is a SMEP defined by

h "(o = sin((bct)e mum (3.23)

where U(t) = [(1, :3, we have

29

1g?) = [wccos(wct) - asin(w€t)]e'°“U(t) (3.241
and
dzdh :———(——t)= (tr2 -(r1 2)sin((o t)- -2a(o cos(w elem/(0111.190) (3-25)
I

Using the relation 1' = t-R‘n—Rs, the left hand side of (3.22) can be written as

62h‘(F’,r) 6111:" t) _ 1
fdevl= aathT —_dv/ - 5L h(F [,1)dV/ (3°26)

Then, using (3.22) and (3.26) together with (3.23) and (3.24), we can obtain the two

equations

 

sin(wct)e'“U(r)dv’ = otsH‘(F ,t)dt2 (3-27)
V

(K r)0091(l3/2)°

and

 

[Vwccos(wct)e “"U(r)dv’= [fo' H’(FS ,t)dr+af0‘fotH‘(Fs,t)dt2] (3.28)

(K f)COS(B/2)

Substituting (3.25) into (3.22) and using (3.27) and (3.28), we have

 

5 dv’- ”3 H‘” 2 '11“ d:
f V (r) - [ (r,.t1+ afo (r,t)

(Kuwait/ZN. (3.29)

+ (1.13 1. a2)fo‘fo'HS(FS,r)d12]

Now, defining the characteristic function

30

 

“ - 1 FIEV (3.30)
709 {0 ”V
we have
on " “ -' / 2firs t ..
fffY(r/)6(ts-r.r/)dv = " [HSO-gat) +2afoHS(rsrt)dt
'°° “”6 (3.31)
2 2 t t s _. 2
1. (w. + 01 )f0 [0 H 03,011: 1
where we use the relation 1: = t-Ri-RS and the properties of the Delta function
5(1) = 5(-t)
D 1 (3.32)
6(2cos—r) = —_ (r)

0
2 2eos(0/2)

r.+r —t
and t = -—'——s——.
‘ 2cosw/2)

It can be see from equation (3.31) and Fig. 3.3 that the right hand side of (3.31)

is the Radon transform of y(F’) (see Section 3.2). It denotes the projected area function

_ ~.-/
at the plane ts - r r

140.1,) = fjfy(F’)o(ts—f'F’)dv’ (3.33)

A

Here A(f,ts) is the projected area onto ts for the particular aspect direction r along ts.

Note that the cross sectional area A(f',ts) is formed with a time scale such that the

cutting plane rs = r‘-F’ used to determine A(f,t) moves with 1/(2cos(B/2)) (one half for

the monostatic case) of the velocity of the incident SMEP.

31

If the view angles are available only in the x- y plane (0’ = rt/2, 05(11‘ 321:),

then the body geometry can be related to the cross-sectional areas A(r‘,ts) through

AW.) = f f [101 ’.y ’,z%(r,—(x’cos¢‘+ y’sin¢‘))dx’dy ’dz’ (3.341

Integrating with respect to z’, (3.34) becomes

110,13) = ffI‘(x’,y’)6(ts—(x’cos¢‘+y’sin¢‘))dx’dy’ (3.35)

where I‘(x,y) = f "y(x,y,z)dz is the "thickness function" of the target in the z direction.

Taking the two-dimensional inverse Radon transform of equation (3.35) (see

Section 3.2.2 equation (3.10)), we can get the thickness function [1 1-13].

 

u .. aA(‘,t,) drsd
1‘03) = ——1—- 2 f r 4’
4n2 ° '°° 5‘. t,-p’008(¢’-<I>)
110.1,)

1 n «-
= __ d: (3.36)
219]" fr (1, — Noose/4112 ‘

 

 

1 r .. COS(BIZ)A(f,t,)
- —?fo dtdd)

'°" (t-13,-,+20’COS(13/2)<:OS(¢l>’-<l>))2

where 5’ = x’f+y’)‘1 p’cosd>’£+p’sind1’y, and pis = 91"93- Substituting (3.31) into

(3.36), we have

 

2 ,. .H‘(61.0+2a ‘H‘(i1"0dr+(wi+a21 ’ 11151111112 (337)
p(§’)=-MB_/Zfof f0 :9 fofo dtdtl)

(ii-r111: w, ‘°' [t-11,-,+20’cos(11/2)<=os(ct>’-<Ii)]2

From equation (3.37), we know ([1 depends on (11‘, (11‘ and [3. Here we only

32

A Z A transrnittin
r i direction g

section cut
by r1“ =ts

bisector
dlrrectron

 

 

 

)Y

1' receiving
“I \ 81......

 

Figure 3.3 Geometry of 3-D body reconstruction problem

33

consider the special case when both transmitting and receiving antennas are fixed and the
target rotating or the scatterer is static and the transmitting and receiving antennas are
moving around the scatterer with bistatic angle (if-(l1s = constant [3. Then (i) = ¢1‘-B/2 (see

Figure 3.2 and Figure 3.3), we have

H’(p‘_:,t) +21: [0'11 vague!) + (1113+1112)f0'f0‘1v‘(aim:2

 

 

P(6,)=_2__P___s :COS(B/2)fpn’2+w2f: dtd¢i
(K p)1t(11 [I'P,-_,+ZPICOS(l3/2)COS(¢I’¢‘+9/2)]2 (3,38)
= _ 2,50,15,27?” H5055)+2afo'H‘(6;.0dt+(wi+a21fo‘fo‘Hiwét1dfiMbi
(if-(5)1: we "' [t-01:2r>’cos(13/2)cos(<ti’-<I>‘+13/2)]2

This is the complete solution to the two-dimensional bistatic problem of recovering
a body from its scattered field; it needs the reflected fields and the correction term from
all possible directions in the x-y plane. The next step is to solve equation (3.38). We can

use an iterative method to get the 2-dimensional target geometry [9]. First, neglect the

correction terms hes,(p",t) and hcg(§’,t) in equation (3.38). This is the time domain
physical optics inverse problem. We can thus obtain an initial estimate of I‘l(p') , and the

correction terms, hcs,(p°’,t) and hc52(p°’,t), in (3.38) can be obtained by solving

76(F,t), .7c2(F,t) in (11) using the marching on in time method, and then used in (3.21).

Then the correction terms can be used in (3.38) to obtain a new estimate of P2(p‘). Then I‘l(p‘)

and 13(5) can be compared to see if the change is less than some small number, and the

procedure continued until this convergence criterion is satisfied. The numerical results for
a sphere of 14 inch radius using those procedures will be shown in Chapter 4.
Note that from (3.22) we can see that when the incident magnetic field is a ramp,

34

using the physical optics approximation by dropping the correction terms, and using
B=0,I?°f=l for the monostatic case, then the backseattered ramp response is proportional

to the cross-sectional area of the target as a function of the distance along the line of the
incident direction [15]. This is the same result which Kennaugh, Cosgriff, and Moffatt had
obtained [6], [7]. Das and Boemer [12] have shown that the size and shape of an object

can be obtained from its area functions, and the problem can be reduced to the classical

Radon problem.

3.3.2 The Monostatic Case
For the monostatic case, the same procedure could be used to derived the inverse

scattering identity as described in Section 3.3.1. We could also use our formulas (3.37)

just obtained for the monostatic identity if we let the bistatic angle B =0, then

908(13/2) = 0
r“), = r“)
(3.39)
9., = 29
1313 = [ii-(r1415? (1 =1
Applying (3.39), (3.37) becomes
H‘(5;,r) +21: fo'H‘(5;,r)dr+(w§+a2) [0' fo'H’(5:,t)dt2 (3.40)

dtd¢‘

 

_. 2p 1! -
1‘ = __
(p6 wef° f" [t-2r>+20’cos(¢1’-(I>‘)]2

Comparing Figure 3.3 and 3.4 will help understanding this identity.

35

6]

ts 'é section cut
1 \ by f0?= ts

 

 

16’) ; . >y

Figure 3.4. Geometry of 3-D monostatic body reconstruction problem.

36

3.3.3 Rotational Symmetric Targets

If the scatterer as shown in fig. 3.4 is a rotationally symmetric target about the z-
axis, then it is the simplest case to reconstruct the body shape of the scatterer, because
the contour of a rotationally symmetric object can be completely described by the radius

p which varies as a function of z, and its projected area function can be expressed as

A = 7.':pZ(z) (3.41)

where p(z) is the radius of the scatterer along 2 axis. Substituting the projected area

function into (3.31) yields the inversion equation for the rotationally symmetric case

1

p(z) = { 2" [H‘(Fs,r)+2a fo‘HS(r;,r)dr+(u§+a2)fo'fo'HS(§,r)d12]}2 (3.42)

 

—.

91,00?)

This expression needs the reflected field for only two incident directions to
reconstruct the rotational body, and we can use the iterative approach or the direct
"marching on in time" approach to solve (3.42). If the upper part and the lower part of

the scatterer is also symmetric, we can simplify (3.42) as

l

2 i s .. t r s _
9(2) = { rs . [Hls(Fs,t)+2a j; H1 (rs,t)dt+(w: +012) [0 foH1('sJ)d‘2]}2 (3.43)

 

—o

wc(K-r)

where H,‘(Fs,t) = h‘(Fs,t) +hcs,(Fs,t). If we use PO approximations and neglect the

correction terms, then (3.43) becomes
1

27' z 2 t t
= _—S h“, +2 h“, dt+ .+ 2 h“, dt2 2 3-44
9(2) {mc(E-f)[ (rs t) afo (rs t) ((0 a )fo [0 rs t) J} ( )

37

 
  
        

h

 

2}!

  

2/2

   
 

ng problem

Figure 3.5 Geometry of rotationally symmetric scatteri

38

3.3.4 Flying Radar Targets

In this special case, we assume that the upper part of the flying object is

symmetric to the lower part of the object such that only the reflected fields from the

lower part of the object are needed to recover the shape of the scatterer [52]. Also note

that the bistatic angle [3 and the distance between the scatterer and the transmitting and

receiving antennas change with the movement of the flying object. In order to simplify

the problem, we assume that the flying object has constant altitude h, and the distance

between the two antennas is a constant (1. When the target flies forward, it is equivalent

to letting the antennas move backwards at the same velocity of the target while keeping

the target static. The geometry of the scatterer is shown in figure 3.6. The following

relations can be found from the figure

 

r‘. = \Jh2+(hcot¢ +;)2

 

rs = J]? +(hcotd) - (51f

 

r = (,[112+(hcot(|>)2

l

 

sind)‘ = Lsind) =
r

 

\J 1+(cot41-2—(2)2

39

(3.45)

(3.46)

(3.47)

(3.48)

flying object

 

 

 

 

1% d 9 ——> V
I receiving

transmitting antenna

antenna

Figure 3.6 Geometry of the flying object

40

NI'OD

 

 

 

 

 

 

 

 

 

 

E = arc sin 1 -d1
2 (3.49)
\J 1+(cot<l> -—)2
. 1
are srn - (I)
l +(cottb — if
2h
arcsin 1 osrb + sind)
+ -—d- 2 + —£ 2
I (catch 2h) \J 1 (cotd) 2h)
. (3.50)
cos (1) + 81nd)

 

 

 

d
l+(°°t¢'§;)2 \J l+(cot¢ - 2%)2

d
1-—sin2
4h 4)

 

 

sindkj 1+(cotd1 - 56%)2

Since the bistatic angle is changes with the movement of the flying object, the

related terms should be inside the integral in the expression of the thickness fimction. The

thickness function can be written as

—11

fifofi

rscos(B/2)[H

3(a’,r)+2afo‘HS(F,’ ,t)dt+((.) 20.62)f ‘o‘f H‘(F find:

2
.. Jdtdd)
K 7‘ [t-ri -rs+2r ’cos( [3/2)cos((l)’-(1))]2

H ”/t
rCOS(B/2) (r dtdd)

 

Kr [—t r. -r +2r ’cos(B/2)cos(<l)’- -(]))]2

(3.51)

41

where

‘-./ : s-J ‘s-o/ 2+2 "s-o/ 2 3.52
H(r,t) H (r,,r)+2afoH (rs,t)dt+((11c a )f0f0H(rs,t)dt ( )

Then, substituting (3.45), (3.46) and (3.50) into (3.51) yields

 

_‘_I ' _ «I + _i 2
(h 4sm2(11)l'l(r,t) l (cotd) 2h) dtdd)

\ (3.53)

2

_ - + _12+ 1 _1 - 1_
(t ris)srn¢\Jl (cotd) 2h) 2r (1 4hsm2(t1)cos((l1 (l1)

 

 

34:17:

 

K-r‘

 

where

 

 

(3.54)

 

 

d 2 d 2
. = h 1 t —— h 1 t ——
r +(CO¢+2h) + +6204) 211)

V

The time-domain identity (3.54) can also be used to construct an image from short

pulse radar measurements.

3.3.5 Three Dimensional Case

Assuming the incident magnetic field is a Half Gaussian Modulated Sine (HGMS)

(1(1) : Sin(wct)€ MHZ) (1(1) = Sin(wct)e ”fill/(t) (3°55)

where a is the shape factor, a =7t/az. the second time differential of h(t) is

d2h(r)

dtz = (4a212_(,):-2a)siant-4awcrcoswctJe-atzum+wc5(t) (3.56)

substituting (3.56) into (3.22) gives

42

[111 [(40:21:2—(03-21!) Sinwcr ~411th coswCTJe -ar2U(r)

(3.57)

 

+(ncb(t)}dv’ = _. firs H ‘(Fs,t)
(K°f)COS(l3/2)

Using the characteristic function and the properties of the Delta function equation (3,57)

can be written as

 

fffY(F’)6(ts-r‘-F)dv’= g 2an ”30.107
.. (K ‘f)cos(l3/2) 9). (3.58)

_2_‘[V[(4122r2 — (a): —2a)sin(ocr -4awctcoswct]e '“T2U(t)dv/
wC

where rs = r-t/2.
It can be seen from equation (3.33) and Fig. 3.3 that the right hand side of (3.58)

is the Radon transform of y(F’). It denotes the projected area function at the plane

rs = f-F’ . In polar coordinate system, (3.33) can be written as

110,23) =f0°°f02“f_n;’2y(r’,e’,¢') 6(ts-r’[sin0’sin0‘cos(¢’-¢‘) + cos0’cosfl‘])r ’Zsine’dfl’dd1’dr
(3.59)

If the cross-section areas A(r",t3) are given for all directions (0',(|>‘) and for all ts, then

taking the 3-D inverse Radon transform of (3.59), we can obtain the body shape function

43

1112 02A(r‘, r)

sin0‘d0‘d41" (3.60)
“/2 air

 

7(r’,6’.¢9 = 0f2”_

[3:77

where rs = f-F’ r’[sin0’sin0‘cos((|>’-¢‘) +cos0’c050‘]. From (3.58) we have

27tr
A(f9 ts) _ _. s HS(F,I)
(K 'f)c08(B/2)w,

 

(3.61)

”(f—[V —-—H[(4a 7: ~01 i-Za) sinw 1: -4a(o 1: cosw re e’“2U(r)dv’

(0C

Since r = t-2r+2f-F’ = —2ts+2f-F’, and dr = dt = -2dts, the second differential of
A(f,ts) with respect to t, is

M .1" 8"“ 32mm
621.2 1(K'f)°°s“”2)“’c 82’

 

_f _8-[( )sinwcr +1;( )cosmc‘r]e MRI/(1)0.“
V (dc

 

 

 

 

 

(3.62)
+8063 +6a)f f fy(F’)6(ts-f-F,)dv’}
_[ 81v: 52 }
—H ‘(F, t) + 8(0): +6a)A(r, t )
1(K comm/2141 621 , 3.1
Substituting (3.62) into (3.60) gives
/ 21: -rs 62,13 . . . .
y’(r ,0 ,(11’)= [0 fi[ H(F, t)— —((o: +6a)A(r, t) sm0‘d0‘dd)‘
3.11:6 (K’ noose/21 612 ,

(3.63)

From (3.61),‘we have

2
A(r, r)| _ _, “'5 HS(F,:)|,=2,_2,.,,., (3.64)
" =(K'f)COS(l3/2)wc

 

Substituting (3.64) into (3.63) gives

 

/// -rs Znn/Z{:2 )HI‘SS -rrr
70.6.41) — f f —+2(u§ +6a) (F ,t)sm0d0d(l1

(K r)cos(B/2)nw 0 -1112 (3.65)

 

[if /2{__2 +2((.1 2+6(11)}(h ‘(Fs,t)+hCS,(Fs,t))sin0‘d0‘dd1‘
K f)cos(B/2)1r(.1 “/2 at
where t = 2r-2r’[sin0’sin0‘cos(¢’-(b‘) +cos0/c080‘].

This is the complete solution to the three dimensional problem of recovering a
body from its scattered field; it needs the reflected fields and the correction term from all
possible directions. If the correction term is dropped, we can get the three dimensional
time-domain PO inverse scattering identity

2n up

110’ 0’d>’)= Tiff f

wCO-n/Z

sin0‘dB‘drl)i (3-66)

r=2r—2f-F’

 

62
2( 2+5 )+— 3(7)
01 a aZ—rih r

 

 

45

CHAPTER 4

TIME-DOMAIN IMAGING OF RADAR TARGETS USING

ALGORITHMS FOR RECONSTRUCTION FROM PROJECTIONS

4.1 Introduction

The basic idea of reconstruction of an image from a series of its projections
appears to have been first discussed by Radon [16]. The techniques that exist for
reconstruction fall into two directions. The whole operation can be done in frequency
space directly, or the equivalent of these expressions can be transformed in the spatial
domain. Whether implemented in the spatial domain or in the frequency domain, the
reconstruction algorithms can be conveniently interpreted by means of a straightforward
and interesting theorem, which is the projection—slice theorem. This theorem states that
the Fourier transform of a projection is a center-cross-section of the Fourier transform of
the projected object. Most of the modern tomographic systems are based on this theorem.

The algorithm that is currently being used in almost all applications of straight ray
tomography is the filtered backprojection algorithm. It has been shown to be extremely
accurate and amenable to fast implementation. We will extend this approach to our radar
inversion problems. In this chapter we will show how to estimate the shape of a radar
target using a short pulse radars. We start with the definition of backprojection and then
develop the reconstruction algorithm based on the inverse scattering identities derived in

46

chapter 3.

4.2 Backprojection

Let (x,y) designate the coordinates of points in the x-y plane and let 1? = (x, y),
E = (cosd), sin(|1). Consideranarbitrary function g(p,E) where p = 52’ = xcos¢ +ysin(b.

The backprojection is defined by [13]

001.1) fo"g(xcos¢ +ysin41, 81441 (4-1)

If the replacements x = recs 0, y = rsin0 are made, then G in polar coordinates becomes

60. 01 = [grime-<11, 411441 (4.2)

Observe that if we consider g(p,(])) is the projection function F m (11) from f(x, y) by
the Radon transform, then for fixed (11, the incremental contribution d(G) to G at the point
(161’). or (r, 0), is given by F(p,(]))d(]1. The full contribution to G at (x, y) is found by

integrating over (i) as indicated in (4.1). The operation for fixed (11 is illustrated in fig. 4.1.

4.3 Reconstruction Algorithm

In chapter 3 we have described a method to extract the cross-sectional area
function of a radar target from SMEP responses and obtained the thickness function. Here

we will illustrate the image reconstruction algorithm.

47

9(Pr¢)

 

 

 

 

Figure 4.1 Geometry for obtaining the backprojection

48

Let’s start with the thickness fimction ms). It can be written as

7: no ”I i
Nb") = -——1—2-f0 .. _ , “7M“ ,_ .- 2 (4.3)
21: [t p.-_.+2r> cos(B/2)cos(¢ ¢+B/2)]

where

= 4p,c08(ltl2)

fifi’J) _. .
(K1091,

[H‘(p‘;,t) + 2afo'H‘(6;,r)dr + (1.13 + 62) fo‘ fo‘H 5(6;,r)d12] (4.4)

From equation (4.3), we know that the inner integral is the convolution of
f(p,-_,-29’cos(13/2)COS(¢’-<I>‘+Bl2)) and l/[pg-ZP’COS(B/2)cos(¢’-¢‘+BIZ)?- Equation

(4.3) thus can be written as

P071) = [07:71-14lam)eprw(p..—2p’cos(—‘2’-1cos(¢’—¢‘+%111dw44>" (45)

since the Fourier transform of 1/[pis-2p’cos([i/2)cos((|1’-(l1‘°+Bl2)]2 is —1t|(11|, and the

Fourier transform of f(pis-p’cos(|3/2)cos(¢’-¢‘+B/2)) is F(0)), and the convolution is the

Fourier transform of a product. We can see from equation (4.5) that the inner integral

represents a filtering operation, where the frequency response of the filter is given by

—7t |(11|. Therefore the inner integral part is called a "filtered projection". The resulting

projections for different angles are then added to form the estimate of I‘(p").

When the projections are bandlimited by the highest frequency B, the projection
data are collected at the Nyquist frequency, with a sampling interval of a = I/(ZB).

Equation (4.5) may be expressed as

49

17(0) = [0" f_:o(w1 Fm) exptiww,-2p’cos(0/21cos(¢’-¢‘+0/2111411140" (491

where

0(a)) = —l—|(r1|rect(w) (4-7)
27c

and

1 |w|<27rB (4.8)

rect((r1)
0 otherwise

Q(01), shown in Figure 4.2, represents the transfer function of a filter with which the
projections must be processed. The impulse response, q(t), of this filter is given by the
inverse Fourier transform of Q(00) and is

1~_1_

Ito |rect((.1)exp(i(r1t)d(r1

q(t) 2n 7°27: (4,9)

Zstinc(21tBt) — str'nc 2(rtBt)

This function is shown in Figure 4.3. Since the projection data are measured with a
sampling interval of 1/(2B), for digital processing the impulse response need only be

known with the same sampling interval. The samples, q(n), of q(r) are given by

 

—1—2 n=0
4a
(101) = t 0 n=even (4-10)
- “2:202 n=odd

 

This filter was first discussed by Ramachandran and Lakshminarayanan [17]. We can

replace integrals in (4.3) by sums and obtain the approximate reconstruction formula

50

given by

N-l on

P011) = 572: Zf(<1‘.t,.)q(t,—p..+2p’cos(0/21cos(¢’-¢‘+0/2» (4°11)
i=0 k=-~

where there are N angles (11i for which the scattered fields are known.

There are two different reconstruction algorithms that can be used depending on
how the filtering is done. If convolution in projection space is used to perform the
filtering, the above reconstruction algorithm is referred to as the convolution
backprojection technique. If the filtering is performed in Fourier space, however, the
algorithm is designated the filtered backprojection technique. Both techniques produce
comparable results.

For the convolution backprojection technique, the reconstruction algorithm
consists of the following steps:
Step 1. Create the cross-sectional area function f(¢1‘, tj) (also called projection) from
measured data, i=1, 2, ..., N, j=-oo, ..., -1, 0, l, +00.
Step 2. Convolve the projection with the discrete impulse response of the filter to
obtain the filtered projections.
Step 3. Take the backprojection for the data obtained in step 3.
For the filtered backprojection technique, the reconstruction algorithm consists of

the following steps:

Step 1. The same as step 1 in the convolution backprojection method.

Step 2. Take l-D Fourier transforms of the projections obtained in step 1.

Step 3. Multiply the transformed projections by the frequency response of the
filter.

51

Step 4. Take the inverse Fourier transform to obtain the filtered projections.
Step 5. Take the backprojection for the data obtained in step 4.

Since the filtered backprojection and convolution backprojection techniques
produce reconstructions of similar quality, we only use the convolution backprojection

algorithm in this work.

4.4 Numerical Results and Images for a Metal Sphere

A metal sphere is the simplest target which can be used to validate the formulas
developed in the previous sections. We use the SMEP as the incident magnetic field
pulse. Figure 4.4 shows the theoretical SMEP generated by taking or=8><109 and
fi.=10GHz in equation (3.23), and the synthesized SMEP using the frequency band 4-
16GHz. The spectrum of the synthesized SMEP is shown in Figure 4.5. Figure 4.6 shows
the SMEP response of the sphere obtained using the marching-on-in-time method
described in chapter 2. We can see that the solution is unstable. Figure 4.7 shows the
results of using the stabilized formula (2.35). The results in Figure 4.7 is obviously much
more stable than that in Figure 4.6. Figure 4.8 shows far-zone scattered field of a 14 inch
sphere computed by using the Mie series and the marching on in time method,

respectively. The cross-sectional areas obtained using equation (3.31) is shown in Figure

4.9. Note in Figure 4.9 that the cross-sectional area A(f,ts) is plotted with a time scale

such that the cutting plane rs = F" used to determine A(F,ts) moves with l/(2cos(B/2))

(one half for the monostatic case) the velocity of the incident SMEP. Figure 4.10 shows
the images of the sphere using the time domain identity (4.11) with the PO

approximation. Figure 4.11 shows the image of the same sphere when the correction

52

terms are considered. We can see that the correction term provides only a small

contribution to the reconstruction.

4.5. Experimental Results and Images for Aircraft

The time-domain identity (4.11) can be used in real time to construct an image
from short pulse radar measurements. A simulation of time-domain imaging is carried out
using data measured in the Michigan State University free-field scattering range. The field
scattered from several aircraft models was measured in the plane of the aircraft wings in
the frequency band 4-16GHz at 200 aspect angles from 0° (nose-on) to 180°, using an
HP8720B network analyzer. All measured data are bistatic with bistatic angle B=10°. The
data from the non-illuminated side was provided through symmetry. Each scattered field
response was first calibrated using a 14 inch diameter sphere as a reference target [see
appendix), multiplied by the SMEP window, and then inverse transformed into the time-
domain using the FFT to provide a SMEP response. The data we used are the derivative
of the measured SMEP response data, and thus a sharpening of the target edges is
provided.

Two different frequency truncated SMEPs haven been used in our radar target
imaging. Figures 4.12, 4.13 and 4.14 show the images of a 1:48 scale model TR-l
aircraft, a 1:32 scale model F-l4 aircraft and a 1:72 scale model B-52 aircraft,
respectively, using the monostatic inverse scattering identity (3.40) under PO
approximation and the synthesized SMEP shown in Figure 4.4. Figure 4.15, 4.16 and 4.17
shows the images of the same target models TR-l, F-14 and B-52, but using the bistatic
identity (3.3 8) under PO approximation. We can see there are improvements in the images
because the measured data are bistatic data. The edges of the aircraft are clearly visible,

53

with the fuselage and vertical stabilizers producing the largest thickness values, as
expected. There is some distortion in the aircraft shapes due to the use of the PO
approximation and an inaccurate estimate of target range. The image obtained using data
from the restricted range of aspect angles 0-90° is shown in figure 4.18. It can be seen
that the shadowed edges of the target are invisible due to the limited view-angles and use
of the physical optics approximation. Figure 4.19 shows another theoretical SMEP
generated by taking or=4><109 and fi=lOGHz in (3.23) and the synthesized SMEP for the
frequency band 7-13GHz. Figure 4.20 shows the spectrum of the synthesized SMEP.
Figure 4.21 shows the image of the same 1:72 scale model B-52 aircraft obtained by
using the second SMEP. We can see that the edges of the aircraft are not as clear as in
Figure 4.17 because the pulsewidth is bigger than that of the first SMEP. We can improve
the quality of the picture by adding a proper window to the SMEP response data in the
time domain. This window is formed by two steps. First, the biggest points in each
SMEP response data are found by comparing the values of the data; these biggest points
actually are the responses from the scattering centers of the target. Then the values in a
small range around the biggest point are set to unity and to 0.5 within the two next
biggest point ranges. Note that different data sets have different windows. Figure 4.22 is
a SMEP response data set before windowing, and the window produced based on Figure
4.22 is shown in figure 23. Figure 24 is the SMEP response after windowing. The image
of the same BSZ aircraft model found by using the windowed data is shown in figure

4.25. By using windowing, we have increased the resolution of the SMEP and obtained

clear images with highly defined edges.

54

Figure 4.2

 

The filter response for the filtered backprojection algorithm. It has been
bandlimited to 1/2a.

55

q(t)
"(4&1

 

 

 

 

 

Figure 4.3 The impulse response of the filter shown in Figure 4.2.

56

 

 

 

 

1.00 -
_.
1 —— theoretical SMEP
I ------ synthesized SMEP
0.50 3
I
z
4
(D I
8 - 1
t 0.00 ': a ‘lfi: """"""""
o 7 '1:
> 2‘ '1
-o.5o 5
I r
- I
-1000-IlrIIIIIIIIII1TII1IIIIITTIIIIIIIIIIIII1j
0.00 1.00 . 2.00 3.00 4.00
Time (ns)

Figure 4.4 Comparison of theoretical (or=8x109, fi=10GHz in (3.23)) and synthesized
(4-I6GHz) SMEP.

57

1.20

1.00

.0 .0 .0
45 CD CD
0 O 0

Voltage Magnitude (relative)
.0
DJ
0

 

b1111111111111111111Llrrrrrirrrlr11111111111111111111111111111

O

4.00 8.00 12.00 16.00 20.00
Frequency (GHZ)

0

Figure 4.5 Synthesized SMEP spectrum.

58

1.50

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.00-j
A j
7o -
Q) ..
.5 I
o J
E 0.50 :
L ..
o . 11
c :
v .1
(11 -
_O 0.00 z
.3 i i]
'61 :
0 :
2-050:

—1.00 ‘ IIIIIllIIIIIIITFTTI[TTTTIIIITIIIIIIIIIIIIIIIIIITI]
1 1.00 1 1.50 12.00 12.50 13.00 13.50

Time (as)

Figure 4.6 Far-zone SMEP response of sphere using marching-on-in-time method.

59

1.50

1.00

0.50

 

 

0.00

 

Maglitude (normalized)

I
.0
or
o

 

 

11111111111111111111111111111111111111111111411J

 

 

-1oOO IIIIIIITI'IIlIITIFIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIT]

11.00 11.50 12.00 12.50 13.00 13.50
Time (ns)

Figure 4.7 Result of using stabilized formula to the scattering field shown in Fig. 4.6.

60

1.50

 

 

 

 

 

 

 

E -— Mie Seris
1.00 7: ----- Marching on in time Method
:3 :
o :
.5 -
6 0.50 2
E :
o -
5 :
a, ..
3 0.00 _
'51 :
O _
5 :
—0.50§ ’
i
—1.00-IIIIIIIIIIIIIIIIIIIIIIWITIIIIIIIIIIIHIIIIIIIIIIII
11.00 1 1.50 12.00 12.50 13.00 13.50

Time (ns)

Figure 4.8 Far-zone SMEP scattered field of a 14 inch sphere

61

1.20

 

1.00
E
'fi 0.80
E.
.8 0.60
g
D.
E 0.40
0.20
0.00 T'TUUITFT YTIUIUIII IIIITFTII ITIFIIIIU
0.00 1.00 2.00 3.00 4.00
Time(ns)
Figure 4.9 Cross-sectional area function of a sphere 14 inches in diameter

62

0.64

O.5~
0.4~
O.3~
O.2~

I v 1,
-’lll’,’ ll

” \ \ 7,
'1, . . ‘ ‘. _
1 ,,.,.;.. ‘ ‘1 .. y .
-‘-' I ’’’’ I I
0.1N _ -,:,’;;;' "I”! ,;,I;:I/’, ’l

‘7‘-_

   

-‘-‘

 

0‘
0.2

  
     

‘ - _ a. .
‘\:: \\.§ \ -

- - ‘ _
‘ ov‘ --- ~;:.. 2'.-
‘\\ %::.s;:_‘:. -.7:. -: :.‘-_ - 7
Q - ‘ ‘. -',“-

Figure 4.10 Upper surface image of the 14 inch diameter sphere using PO
approximation.

63

O.6~

0.5w

0.4u

  
    

0.24

 
  

$ \\
I G
, H - 'I!) I, ” 0’. o . O ‘D “““\‘\\\\\
n :9, . I, I” ”O Q...‘. ‘\ \\\\\\ \
’ ’I7””1’9’0'4'0'0’03?“

     

 

          

0.1~ ;, _. I], )’Q a \ \\\‘
- '7’ ’ ' 0 ”5": %’I’I’Il’.' "'."""‘..‘G“.“\:‘ “ \\\‘\ \
I " " "'0: ’ "' " 9 '0 O O““‘ 0 0 Q“ ‘\\\““‘ ‘ ‘ ..::::‘:.~~:-
0 7 I."TI”:”:::':::::::'::'::':':'o'O'o'o’e’9‘.o‘.e“:“:::“:::::: “fiver: “55:: -‘ -
,' ’0. oo'o'o. . o 0.. O O . O. .‘ .. .‘ “‘ ““‘ .‘-\: :“::-‘--:_---
0.2 I I'o'.:§'::.:.:.:...:o:. .o.:.:::.‘:.::c:.‘,«a»; ~ - v: ......
O ..

Figure 4.11 Upper surface image of the 14 inch diameter sphere considering the
correction terms.

 

Figure 4.12 Image of TR-l from 0°-180° data in band 4-16 GHz using
monostatic identity.

65

~ 1"1‘ M
1: - 1
WW
"I:

 

Figure 4.13 Image of F-I4 from 0°-l80° data in band 4-16 6112 using
monostatic identity.

66

   
  

 

' r 1 ‘
I III -:
||."1."!"

1, w. '1

”fl",
1’

    

Figure 4.14. Image of B-52 from 0°-180° data in band 4—16 GHz using monostatic
identity.

67

.30 -

.20 -

I10“

.00 .,

.010 '1

 

'120 «

 

.030 II

 

3‘10 .030 .020 7010 .00 0‘10 02° 03° 0‘10

Figure 4.15 Image of TR-l from 0°-180° data in band 4-16 GHz using
bistatic identity.

68

”1° 1

.30 -

I20 ‘

.10-

'.10 -

'.20 -

'.30 -

 

unoJ

'.«,y‘,‘.“&“'HIl¥ 9“ ‘~ .0

‘ jun
“I fl If.:l:;'.
‘.:' “fl: 1..
IV

.. as
“fig" L,
’ u“;- ”I .
.'4.'. 11'.

‘3' ‘5".

       
    
   

 

' . 'I‘: a
’ , \‘u/i‘, . ‘x ._
. i " .-"'-’--‘ f i: “\\.
”a; A - g . 1\._
-". ' ’ ‘ 'u ‘ l: | ."‘
». h ‘
.0" ., “3"". am» - ‘
' ' ,. bu?“ ‘
s; , - ’2‘ .e
in ..1....,.‘..:. - . J '- 12v“ ‘ , 0- a..- 0
.. .—‘-»w ,' ‘ ~.q. :«u—v-o—‘A
.‘ ‘ I »
V (.ln _ 4 n”-
l‘h " II'IHy I
' 'J
¢'ff; ulslzi‘ 9:...
~ v 5 . .

1. y 45".» n33;

 

-IH°

'-80 '020 '110 .0 '10 I20 I30 ”10

Figure 4.16 Image of F-14 from O°-l80° data in band 4-16 GHz using bistatic

identity.

69

.‘10 -

  

         
 

 

 

 

.30 .
.20 .
.xo - -:
1:“ 3 \. ~. A“
.o. ' ' 2.} .~ -., (». m
_ 4 r .w-M . - A
. :3 ,. *r‘ -.
'.10 ‘ _ .‘
I:
l ' \
z. - * ;:
-.2o . ‘5
-.3o .
-.~.o .
-.so -.~10 -.so -.20 -.10 .oo .10 .20 .so .~.0 .50

Figure 4.17. Image of B-52 from 0°-180° data in band 4-16 GHz using bistatic identity

70

 

.‘10 -

 

 

 

 

13° -

02° -

01° -

0° 4

'40 ‘

'020 a

'030 - ~
, . -, ‘2‘“ \
<'.1,"'53?"”’}¢ ' ‘ ‘ - ~ . ‘~\‘t&\‘=
..‘g _V .4. 4, I“, , . ' . . . ,1 ‘ . ‘3 ."I r .

awoj “35“”: 759$». , g . g - .xt‘x‘a

'-5° ’.‘10 .030 '02° '-1° :00 :10 :20 .30 .‘00 '5'

Figure 4.18. Image of B-52 from 0°-90° data in band 4-16 GHz using bistatic identity.

71

 

 

 

 

 

 

 

 

 

1.00 :
Z —— theoretical SMEP
Z ------ synthesized SM SF
0.50 —
a 3 -
g 0.00 ‘ . “9:33;"; .. .............
0 fl ' g,
> 2‘ ::
2 i
-0.50 -
—1000-IUIIIIUUTIIU'UIITIIllierI'IIITITIIIYrI]
0.00 1.00 2.00 3.00 4.00
Time (ns)

Figure 4.19 Comparison of theoretical (a=4x109, fi=10GHz in (3.23) and synthesized
(7-13GHz) SMEP.

72

1.20

1.00

.0
an
o

.0
4;
0

Voltage Magnitude (relative)
.0
50
O

O

O)

O
lllllllllllllllllllllllllllllllllllllllllllljlllllllllllllll

 

 

0.00 WIIIjTTWIIIIIIrITIIIITITII1TIIIIIIIIIII]

6.00 8.00 10.00 12.00 14.00
Frequency (GHz)

Figure 4.20 The spectrum of the synthesized SMEP.

73

.‘10 -

03° 1

.20 -

01°!

‘.10-

-020 It

‘.30 J

 

'.‘-10 4

  
  

 
  

   
 

'u I“ I; . 'le.' . x.
' ' ' I.” [HII' "

l“. .‘ '“'°""'W"w
NIH”; 'lln' ‘0

f} , ' 'W'
' .0! n'~.‘

 

5:. fii‘-'".. .:.i

"(hag-““7“" O . ‘

   

    

 

'.SO

'.‘-00 '.30 '.20 '.10 .00 .10 .20 .30 .‘10 .50

Figure 4.21. Image of B-52 from 0°-l80° data in band 7-13 GHz using bistatic identity

(before windowing).

74

1.00

.0 9
O 01
O O

Voltage (relative)
s':
8

-1.00

LLlllllllllllllllllllIlllIllILlllllllllllllllllllll

 

—1.50
5. 7.00 9.00 11.00 13.00 15.00

Time (ns)

0

Figure 4.22. 8-52 SMEP response (nose-on) before windowing.

75

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.003 F F r- [- ”l-
3 “ii
a) 0.50:
0" 4
0 2
t -
é) :
0.00-j
“0.50 '|I||llllll|llllIlll||lll|ll ll||llllllll|ll|llllll|lllll
5.00 7.00 9.00 11.00 13.00 15.00

Time (ns)

Figure 4.23. Window produced based on Figure 4.22.

76

—I
'o
o

0.50

P
o
0

Voltage (relative)
9':
‘0’I

—1.00

LLllllIllllllllllIllllllllllllllllllllllllllllllllI

-1.50 'rn'rrn-nTrrnTrrrrrrrrnfi-ranfiTn-nTn-fiwrrrrfi-n
5.00 7.00 9.00 11.00 13.00 15.00

Time (ns)

 

Figure 4.24. B-52 SMEP response (nose-on) afier windowing.

77

.‘10 -

.30 .

.20 -

.10.

'010 I

'.20 -

'.30 -

 

 

 

'.50 5‘10 '.30 '.20 510 .00 .10 .20 .30 .‘10 .50

Figure 4.25. Image of 8-52 from 0°-l80° data in band 7-13 GHz using bistatic identity
(after windowing).

78

CHAPTER 5

IMAGING OF RADAR TARGETS USING LIMITED-VIEW DATA

5.] Introduction

In Chapter 3 and Chapter 4 we have shown that good images of radar targets may
be obtained given enough high quality data over a 180° angular range. When some data
are missing, the reconstructed image suffers and may be unsatisfactory. In practical cases,
however, it is impossible to obtain enough data over a full 1800 angular range. The
limited-view problem occurs when the data are available over an angular range less than
180°, and the sparse-angle problem occurs when only a small number of angles evenly
space over 180° are available.

The limited-view problems has attracted considerable attention in computed
tomography imaging research, and techniques for dealing with it have been proposed.
These techniques can be put into two categories: transform techniques that incorporate no
a priori information, and finite series expansion methods that may incorporate a priori
information as constraints. The transform techniques are usually single-pass, direct
reconstructions [18-20], while the finite series expansion methods are usually iterative
[21-27]. Projection onto convex sets (POCS) in particular has shown great flexibility in
dealing with known geometric constrains and with noise. In Section 5.2, we will extend
POCS into radar target imaging.

79

In addition to the POCS method, we will demonstrate a new reconstruction
algorithm for radar imaging in the limited-angle case. The goal of this approach is to
recover the sinogram from available measured data using linear prediction. Since the
scattered field of a target can be written as a superposition of distinct specular reflections
arising from scattering centers on the target, the trace of the scattering centers can be
predicted using linear prediction with the change of the observation angle. Thus, the
missing data may be predicted before reconstructing the image. This reconstruction

algorithm will be described in Section 5.3.

5.2 The Method of Projections Onto Convex Sets (POCS)

The theory of convex projections developed by Bregman [22] and Gubin [23] was
first applied to image processing by Youla and Webb [24]. Since then, various researchers
has used POCS to restore computed tomography imagery based an incomplete data. To
the author’s knowledge, this method has not been used in radar imaging. Here we will

give a brief review of POCS and apply it to radar imaging.

5.2.1 Basic Ideas
Assume a Hilbert space H consisting of the set L2(Q) of square-integrable

complex-valued functions over QER2 (R2 being the real plane). The inner product and

norm in {H are

0.22) = f f g(x,y)h ‘(xmxdy
a (5.1)

Mg“ = [09.15.01"2

The image to be restored f(x, y) is assumed to be an element of CH; Every a priori

80

known property of the unknown feil-C is formulated as a constraint that restricts f to lie in

a closed convex set C,. If m properties are known, there exist m constraint sets Ci

(i=1,2,...,m) and f = ECO =ri,.”:,c,.. Then our problem is to find a point in CO given the sets

Ci and the projection operators Pi (i=1,2,...,m) projecting onto the various C‘. A recursion
relation given by

fk,l = Pum_1...Plfk k=0,1,2,... (5.2)

with f() being arbitrary, weakly converges to a point of Co. More generally, various Pi may
be relaxed by

11,, = Tme_,...Tlf, k=0,1,2,... (5.3)

where T. = I +Ai(Pl.-l), O<Ai<2 and I is the identity operator. The A. (i=1,2,...,m) are

l
relaxation parameters and can be used to accelerate the rate of convergence. The

projection operator Pi projecting onto Ci is defined by the minimality criterion

ILf-Pfll = minlLf-S ll (5,4)

geCi

Equation (5.4) means that the projection Pf is geCi that minimized lLf—g II. A symbolic

representative of a projection is shown in Fig. 5.1.

5.2.2 Convex Projections in Short-Pulse Radar Imaging
In chapter 3 we have shown that the cross sectional area function and the thickness

function of a target form a two dimensional Radon transform pair and can be written as

81

Ilg-fll

 

Figure 5.1 A symbolic representation of the projection off onto C.

82

141404)) = f f 1‘(x.y)5(ts-(xcoscb +ysin¢))dxdy (5.5)

In a discrete case where the thickness function is defined over a cartesian grid of N,‘><Ny

pixels, (5.5) becomes

N N

X

AI‘(ts’¢) = Z

7
=1 j=l

I‘(xpyj)Lk(xpyj) (5'6)

where Lk(xi,yj) is the length of the kth line (ray) through pixel (xi,yj) as shown in Fig. 5.2.
For N views and K lines per view, the total number of discrete cross sectional area
fimctions (call it raysum) is M=NK. With a single subscript k (k=1,2,...,M), each raysum

can be written as

Al‘(tsk’¢k) = {11.03.40 (5.7)

Let Ci be the set of all functions hefl-f such that for a given (tsbok) they have given
raysum values equal to that of the true thickness value F(x,y). Denote this value by
Ar(tsk,¢k). Then

C.- = { h: A..(t .0.) = Ap(t,,,,¢,,) } k=1,2,...,M (5.8)

It has been shown that these constraint sets are closed and convex [26]. The

projection operators that project onto Ci are given by

, g(x.y) for 86C.-

+ Al‘(tsk’¢k) "A303 Ah)

1vx Ny

z 2 Li2(xm.)’,,)

m=l MI

(5.9)

 

P3000 : J 8(X,y) [q(st) fOr gECi

 

83

Other important constraints that can be used in radar image reconstruction are

listed below.
1) Let C A be the set of all real thickness function F’s in CH whose amplitudes must
lie in a prescribed closed interval [a,b]; aZO, b>0, a<b. So, the amplitude limit

constraint set is given by

 

i a g(x.y)<a
P.g = i g(x.y) asg(x.y)sb (5‘10)

. b g(x.y)>b

2) The bounded support constraint set is defined by
Cb = 1 h: 110‘.” = 0) (x,y)“ } (5°11)

This means that Cb is the set of all functions in {H that vanish outside a region A. The

projection g onto CI) is given by

8(x.y) if (x,y) 6A
P,g(x.y) = (5'12)
0 if (x,y)EA
3) The reference image constraint set
C. = l h: IIh-lelseR 1 (5-13)

In words, The reference image constraint set is used to restrict the reconstructed image

to lying within a distance 8R from f}. The projection operator onto CR is given by

f

8 lg -fR|l 56,.

PR8 = g'fk
8"53
“g’fkll

(5.14)

 

Ilg-fR">€R

 

84

I‘(tS! (b)

 

Figure 5.2 Discrete object, F(x,y), and its projection are shown for an angle 41.

85

See reference [26] for the derivation of equation (5.14). The effectiveness of this

constraint strongly depends on ER and fR. It can be used to reduce the size of the feasible

solution set C0.

5.2.3 Experimental Results

The performance of the proposed POCS method was investigated with a 1:32 scale
model F-14 aircraft shown in Figure 5.3. The aircraft image F(x,y) was calculated on a
320x 320 grid of rectangularly sampled points in the x-y plane by the cross sectional area
function A(ts,¢). The cross sectional area function was computed from measured data by
equation (3.31) for various angles.

Equation (5.2) defines the general pure projection algorithm used in POCS. The

specific POCS algorithm used in this experiment can be written as

f...l = PRPBPAPM...P1 fk (5.15)

with an initial point

f0 = 0 (5.16)

where PM, PI is the series of N projections onto N closed convex sets Ci (i=1,2,....M)
defined as in (5.9). The other projections in (5.15) have been defined in section 5.2.2.

Figure 5.4 shows the reconstruction image of aircraft model F-14 after 1 POCS
iteration without use of reference image constraints in the angular range 0-180° with step
angle A¢=O.9°. Figure 5.5 shows the reconstruction image of the same model after 5
POCS iterations using the same constraints as in Fig. 5.4. As seen in Figure 5.5, the
iteration does little to improve the reconstructions.

To demonstrate the power of using a priori constraints in the POCS approach. we

86

need to consider reconstruction in the limited-view case. Figure 5.6 shows reconstruction
images after 5 POCS iterations using Figure 5.5 as the reference image in the angular
range ¢=45°~135° with step angle A¢=O.9°. For convenience, an image obtained with the
convolution backprojection directly on the incomplete data is shown in Figure 5.7. We
can see that a significant improvement of reconstruction was achieved by the POCS

algorithm.

5.3 Radar Image Reconstruction Using Sinogram Restoration and Linear

Prediction

Image reconstruction from incomplete projections can be formulated as a
sinogram-recovery problem. The sinogram recovery problem is to find a complete
sinogram A(ts, 11>) based on the measured incomplete sinogram A’(ts, (b) and a priori
knowledge about the sinogram. Once an estimation of the complete sinogram is obtained,
image reconstruction by the ordinary convolution backprojection is possible.

In this section, we will describe a sinogram-recovery method using linear
prediction. In Section 5.3.1 we define the sinogram. In Section 5.3.2 the method of linear
prediction is briefly reviewed and a solution to the sinogram-recovery problem is
provided. In Section 5.3.3 we present experimental results that demonstrate the

performance of the proposed method.
5.3.1 Sinograms

Projection data (cross-sectional area fimctions) used for image reconstruction can

be arranged in a two-dimensional map in which one of the coordinates is the distance of

87

 

 

Figure 5.3 Geometry of discrete test image to be reconstructed.

88

  
    
    

    
 

I‘O ' . . . 0 .
. -.l d"? 'I ”l' ‘ I n '. .
.33 4.5., ~.‘-:.’-*"::; '1}. £109”? 3:: _ j .o
I -.' hl‘.-' V"‘ ‘ ‘ “..L» ‘ it £0.4i“ . 0
~ .‘L'. {a "a. ‘0 ‘ 6‘ 4" .‘l ',:. "new“ '
,3... .,- ‘i; - _ ,a.-'( “

 

.001

33

4|.

.10.

‘ ilk-{:2'

- a“?
' .

axo. f. I,

'020fi

 

 

“ .‘q r;- ‘ r.
- T! V. -"
"col ‘ .. 8‘ , 11 .
I . ‘0 - ‘
. 6;. 112. ”it. .5. - .s ,,
3’10 4 ,. I c.2313. nfl...-M0‘I'.‘Ifi‘d' .

 

'0‘0 .030 '020 ’.I.0 .0 0:0 02° 030 0‘0

Figure 5.4 Image of F-14 from O°-180° data in band 4-16 GHz using POCS afier l
iteration.

89

    
   

    
  

 

 

.. - I g ‘ Hm _ .
W‘s? 'ai'T‘as‘ a??? o
"’4: ‘V , 3 July
.00. 00“?“ -343!" , I. ‘15: ’- ’ ‘4.
n. -: qt a
v V )‘sf "1“
.20. \~ as
no.
.0. “.5 .
.... -‘ 3%.
' q "i .1“.
,1- ,. -.
. , . . .
- ao 1.1.511. '5‘ 7* ,
HM; .m’
moo]
-.-.o
-.~.'o -.50 -.20 a; 20 .30 .-;o

Figure 5.5 Image of F-l4 from 0°-I80° in band 4-16 GHz using POCS afier 5
iterations.

90

Figure 5.6

.‘10 .

‘.10.

 

3S

 

 

ado -.20 -.1'0 .'o .00 .20 .30 .-.0

Image of F-I4 from 45°-135° data in band 4—16 GHz using POCS after 5
iterations.

91

.‘00 .

.30 .

.20 .

'.10.

'.20 I

'.30-

 

35

 

 

H30 ’.20 ‘.10 .0 .10 .20 .30 .‘1

Figure 5.7 Image of F-14 from 45°-l35° data in band 4-16 6112 using convolution
backprojection algorithm.

92

the wave along which the line integral is taken from the center of the rotation of the
projection system, and the other coordinate is the angle of the wave. In this map, waves
through a fixed point in the object correspond to a sinusoidal curve, which is why a
display of this map is called a sinogram [28].

Referring to the geometry of Figure 5.8, we have defined the 2-D Radon transform
by equation (3.35) in chapter 3. In this section, we assume F(x,y) to be a real function
defined on the disk of radius T centered at the origin. A sinogram is an image of the 2-D
Radon transform, where ts and 4) form the horizontal and vertical axes, respectively, of
a cartesian coordinate system. Because of the periodicity of the 2-D Radon transform and
because of the assumed domain of l"(x,y), the sinogram can be defined over the complete

domain

0 = {(t..<l>)lt.e[-T.TJ, ¢e[0,n] } (517)

For the limited-view problem, the sinogram is measurable over a domain C, where

Q is assumed to be a subset of 11, i.e., @611. C can be written as

C = { (t3’¢)lt56[ -7277], ¢El¢L9n -¢L] } (5°18)

where 0L is a constant that represents the data’s missing angular range. In Fig. 5.9 we

illustrate the measurement domains C given by (5.18).

Figure 5.10 and Figure 5.11 show a 1:48 scale model TR-l aircraft and its
sinogram. Note that the sinogram is formed not from the cross sectional function but from
the SMEP responses of the target. We can see from the sinogram that most of the ’lines’
looks sinusoidal. Those lines actually are the traces of the scattering centers on the target.
This motivates the idea of predicting the traces of the scattering centers using linear
prediction. The following sections give a more detailed presentation of the approach.

93

 

 

projection
A (0.0)

  

 

Figure 5.8 The geometry of the 2-D Radon transform.

94

 

 

 

 

 

 

 

Figure 5.10 A 1:48 scale model TR-l aircraft.

96

(wet?

 

¢-d’

The sinogram of TR-I aircrafi. 35 dB range, 1112 colors.

Figure 5.11

97

5.3.2 Linear Prediction
Linear prediction is a particularly important topic in digital signal processing with
applications in a variety of areas, such as speech signal processing, image processing, and
noise suppression in communication systems [29, 30]. In digital signal processing, the
signal 32,, is modeled as a linear combination of its past values and some input x"
N
y" = {j bij..." (5.19)
j=l
Equation (5.19) shows that the signal y,, is predictable from linear combinations of past
outputs. We can consider x" as the discrepancy of the prediction at time step n, i.e. the
amount which must be added to the predicted value (the sum) to give the true value y".
Equation (5.19) can also be specified in the frequency domain by taking the z

transform on both sides of (5.19). If H(z) is the transfer fimction of the system, then we

have from (5.19)

H(z) = fl = ——1
X(z) N . (5.20)
1+2 biz”
i=1
where
ya) = z: ynZ-n (5.21)

is the z transform of y", X(z) is the z transform of x", and H(z) in (5.20) is the all-pole
model (autoregressive model). This model is shown in Fig. 5.12 in the time and frequency
domains. Given a particular signal y... the problem is to determine the predictor

coefficients bj in some manner. The derivations will be given using the least squares

98

 

 

 

1
H(Z)- y n

_L

"'13)
U
N

 

 

 

Yn

 

Yn .1

, linear predictor . Z-t ‘
oforderN

(b)

 

 

 

 

 

 

 

 

Figure 5.12 (a) Discrete all-pole model in the frequency domain. (b) Discrete all-pole
model in the time domain.

99

approach.

The prediction error x, can be written as
N
xn = y" +2 bjyn-j (5.22)
j=l

In the method of least squares, the parameters bj are obtained as a result of the
minimization of the mean with respect to each of the parameters. The total squared error
E is

N 2
E = 2 e: = 2 [ya +2 bjynv.) (5-23)
n j=1

The minimization of E given by setting

6E

_ = 0 o=l,2,...,N 5.24
ab. 1 ( )
1
leads to the set of linear equations
N
Z]: biz: yn-jyn-k = “Z by"-.. k=l,2,1,...,N (5.25)
J= n ,.

Letting yy(k) = Z Ynynw (5.25) can be simplified as

N
M") = ‘2 bit/yaw) k=l,2,...,N (5.26)
i=1

These are called the normal equations for the coefficients of the linear predictor. The

minimum mean-square prediction error is simply
N
E. = 7.0) +2 bjyy(-J') (5.27)
j=l
by expanding (5.23) and substituting (5.26), the normal equations may be expressed in

100

the compact form
N
2 3'00“!) = 0 J'=1.2.....N with b0=1 (5.2s)
j=0

The resulting minimum mean-square prediction error is given by (5.27). If we
augment (5.27) by the normal equations given by (5.28), we obtain the set of augmented

normal equations, which may be expressed as

N .
E ]=
_ _ = .. (5.29)
12; 1’17?“ 1) l0 j=l,2,...,N

OI'

”v,(0) v.0) Y,(N-1)"1' 05
7,0) yy(0) yy(N-2) b1 0 (5.30)

 

 

 

 

 

 

_y,(N-1) yy(Iv—2) y,(0) [b]. .0.

The matrix in (5.30) is a symmetric Toeplitz matrix, i.e. one whose elements are constant
along diagonals. There are several computationally efficient algorithms for solving the
normal equations [30]. In this work we use the Levinson-Durbin algorithm.

Afier solving for the linear prediction coefficients bj, the next issue which needs
to be considered is stability. The condition that (5.19) be stable as a linear predictor is

that II'IC characteristic polynomial in equation (5.20), INC
5 3
N- 2 : b N" = 0 . 1
Z i=1 jZ " ( )

has all N of its roots inside the unit circle, [2 ISI. There is no guarantee that the
coefficients produced by the Levinson-Durbin algorithm will satisfy this condition. When
instability is a problem, we have to modify the linear prediction coefficients. We do this

101

by the following three steps:

1) Solving (numerically) equation (5.31) for its N complex roots.

2) Moving the roots which lie outside unit circle to an appropriate position inside or
on the unit circle.

3) Reconstituting the modified linear prediction coefficient.
Another consideration is the choice of N, the number of poles. One should choose

N to be as small as practicable for you. We choose N between 5 and 15, depending on

the number of data used.

5.3.3 Experiment Results

Consider the sinogram of the TR-l aircraft model shown in Fig. 5.11. When only
part of the sinogram (domain C) as shown in Fig. 5.13 (for example 45°S¢<l35°) is
measurable, then the reconstructed image may be seriously distorted if we simply let the
missing data be zero. Our goal is to restore the complete sinogram (domain 11) using
linear prediction techniques before reconstructing the image by the convolution
backprojection algorithm. The techniques which we proposed can be summarized as the
following steps:

Step 1. Given a set of measured data P¢‘(tj) (¢LS¢i<n-¢L), i=1, 2, N¢, j=l, 2,

N‘, where N¢ is the number of views and N, is the number of samples in
each view, find the biggest point (the darkest points in Fig. 5.11) in each
SMEP responds by comparing the values of the data in each view. These
biggest points are actually the positions of the scattering centers of the

target. Those points with large amplitude form a new sinogram. Fig. 5.14

102

Step 2.

Step 3.

Step 4.

Step 5.

shows the measurable sinogram of the TR-l formed by these points, and
Fig. 5.15 gives the complete sinogram of the TR-l formed by these points.
From the new sinogram, find the traces of the scattering centers. Then one
obtains a set of functions, each function designating the movement of a
scattering center with the change of view. In Fig 5.15, for example,
number 3 and number 4 denote the movements of the tail and the head of
the aircraft model, respectively, number 1 and number 2 denote the
movements of two wings, and number 5 denotes the movement of the
engines.

Predict the movement of the scattering centers over the whole sinogram
region 7] using linear prediction techniques. Fig. 5.16 shows the predicted
traces of the scattering centers. Fig. 5.17 gives the trace of the tail with
measured and predicted data, and Fig. 5.18 shows the traces of the engines
with measured data and predicted data.

Predict the amplitude of the scattered field at the positions of the scattering
centers. We use linear prediction twice, once for the traces of the
scattering centers, and once for the amplitude of the scattered filed at the
scattering centers. Fig. 5.19 and Fig. 5.20 gives the measured and predicted
amplitudes of the scattered field at the positions of the tail and the engines.
Reconstruct the sinogram of the target over the domain 1]. For each view
over the missing region, find those predicted scattering centers. Then make
a small window (rectangular, Gaussian, or cosine taper window, we use
rectangular window in this work) around each scattering point and multiply

the window by a SMEP defined in (3.23). Then the predicted SMEP

103

responses are found for this view. Fig 5.21 shows the measured data and
the predicted data for observation angle 0:40.5°.

Step 6. After predicting all the SMEP responses over the missing data region in
domain 11, we can restore the complete sinogram over domain 11. The
restored sinogram is shown in Fig 5.22.

Step 7. Reconstruct the image by convolution backprojection using the restored
sinogram. Fig 5.23 gives the image of TR-l aircraft model using the
restored sinogram, and Fig 5.24 gives the image of TR-l using the original
incomplete sinogram shown in Fig. 5.13. We can see that there is a big
improvement in the quality of the reconstructed images by using the

sinogram restoration techniques.

5.4 Conclusions

For the limited-view problem, we have proposed two techniques to handle this
practical situation. One of the approaches is the method of projections onto convex sets
(POCS). The basic principle of POCS is that each piece of a priori knowledge must be
represented by a convex set onto which the current image estimate can be projected. It
has been shown that if the intersection of these convex sets is nonempty then the
sequence of cyclic projections will converge weakly to a point in this intersection. We
extend this approach to radar imaging for the first time and show some usefiil results.
Another approach which we have demonstrated is to process the available measured
projections in order to generate an estimate of the full set of projections, an image which
is called a sinogram. The goal of this approach is to recover the sinogram from the
_ available measured data using linear prediction. Since the scattered field of a target can

104

be written as a superposition of distinct specular reflections arising from scattering centers
on the target, the position and strength of the scattering centers can be predicted using
linear prediction with the change of the observation angle. Thus, the missing data can be
predicted before reconstructing the image. A big improvement in image reconstruction has

been achieved using this technique, and some useful results have been provided.

105

 

 

 

Figure 5.13 Measurable part of the sinogram of TR-l aircraft.

106

 

Figure 5.14 Measurable part of the sinogram of TR-l formed by the positions of the
scattering centers.

107

\ i, _ I". ”. . . J " ‘ .0
l ‘ ' 1" I . . ' ..
-:‘//i i \ x \1 . (u!)
.‘. \.\ . \ \ . . -‘ \_ "I. /\ \

H.‘ ‘\ \ '. . . I, / l

'- . \‘ \ \ \ _. \ ‘\ \ '

. _ \-. a, .. u / \ ' .

/ \ \\ k ' / k I \ \ fig

 

‘ '. ' \ ,- :' m5
// .. \,#3-

#4 / /-_’ V. . l .
l I \, I
i I /.\ , ." ‘

Figure 5.15 Sinogram of TR—l formed by the positions of the scattering centers.

108

 

I ' I '

i. l
_./ l I). '
l |

., 1
t l '1. .

I . l f
,l
I
.. \. x‘
.I x N I";

 

Figure 5.16 Sinogram of TR-l formed by the predicted positions of the scattering
centers.

109

—- measured trace of scattering center ’3
— — predicted trace of scattering center 3

2.00

1.50

1.00

0.50

Down Range (relative)

—0.00

 

111111111IILlllljilllllLlllllIlLlllllllIllllJlLllI

 

 

-O.50 iIIIIIITIIFTTTUIITTIYWITTIIIIIIIIIIIIII]

0.00 0.50 1.00 1.50 2.00
Rotation Angle (0—180 degrees)

Figure 5.17 Measured and predicted positions of the tail of the TR-l aircraft.

110

0.40

—- measured trace of scattering center ’5
- — predicted trace of scattering center 5

0.30

Down Range (relative)
o
N
o

0.10

 

 

‘~_A

0.00 IIUIIIIIFIUITjTIIIIIIIIITIIITIIIIIIIIIT]

0 0.20 0.40 0.60 0.80
Rotation Angle (72—0 degrees)

1L1!lIIILIIILJJIIIIIIIII[11111111111111]

 

 

.0
0

Figure 5.18 Measured and predicted positions of the engines of the TR-l.

lll

— measured amplitude of scattering center 3
— — predicted amplitude of scattering center 3

2.000E-002
1.500E-002
1.000E—002

5.000E—003

Down Range (relative)

\
/

8.674E—019

 

lLiiJiuJililiriiiiilLiii1111111111411111111111111]

-5.000E-003 '11me
0. 0 0.50 1.00 1.50 2.00 2.5

Rotation Angle (0—180 degrees)

C

Figure 5.19 Measured and predicted amplitude of the scattered field in the positions of
the tail of the TR-I.

112

-— measured amplitude of scattering center 5

 

 

0-01 j - - predicted amplitude of scattering center 5
0.01 -
’6? :
.2 -
*5 I
E) 2
30.00 1
C _.
o 5
CK -
c I
g -
0.00 -: ~~~~~~
I
0.00 rfiITer'I[III'IT'lIrjIIIII'VIIIIITrfIrj
0.00 0.20 0.40 0.60 0.80

Rotation Angle (72-0 degrees)

Figure 5.20 Measured and predicted amplitude of the scattered filed in the positions of
the engines of the TR-l.

113

 

 

 

 

4'OOE—OO3 : Measured data
3 — - Predicted data
2005-003 f
3 l
a) I I
'O - I” ‘ ll
3 : I'll ‘ l l‘
‘E 8.67E—019 : ‘ . 1, I, , . M .
O _ I, .,r
O - 'll I
2 j I,
I I
Z I.
—2.00E—OO3 :
I
-4000E—003 III1T'ITI]IIIIFTIIIITTUIIIUTIIIIIITITII
0.00 2.00 4.00 6.00 8.C

Time (ns)

Figure 5.21 Measured and predicted data of the TR-l for o=40.5°.

114

  
    
 

-: a} ‘ ,
fisig-d.‘ “
s“

"“E‘Wfi‘ 35
\\ f “fifth

Figure 5.22 Restored sinogram of the TR-l by linear prediction method.

115

 

 

 

'.35 -.25 '.15 00.5. 9.; .15 .25 .35

Figure 5.23 Image of the TR-l using restored sinogram.

116

.30 -

 

 

 

I0
.20 . m
.10 .
es
loo I ;f\ I!
' l‘. ”I! i
a?) 1m
-.10 .
'020 I
-.30 " ‘3' i
-.Ho -.30 -.2o -.10 .0 .10 .20 .30 .‘10

Figure 5.24 Image of the TR-l from 45°-135° data using convolution backprojection
algorithm.

117

CHAPTER 6
IMAGING UNDERSTANDING AND MODIFICATION USING

HIGH FREQUENCY APPROXIMATION METHODS

6.1 Introduction

Time-domain imaging is an imaging technique that uses the impulse response or
SMEP responses of the object to get the shapes or constitutive characteristics of the object
[31,10]. In the laboratory implementation of this imaging system. an object is seated on
a rotating pedestal and is illuminated by a plane wave. For each aspect angle a set of
pulses at different frequencies is transmitted and its echoes received. The object is then
rotated and the measurement is repeated to obtain the multiaspect stepped frequency
response of the scattering object. Each scattered field response is deconvolved using a 14
inch diameter sphere as a reference target, and inverse transformed into the time domain
using the FFT to provide a band-limited impulse response [10], or multiplied by the
SMEP window first and then inverse transformed into the time domain to get a SMEP
response [31].

In this time domain imaging technique, the physical optics approximation may be
used to model the scattering field of an aircraft [10, 31]. However, the PO approximation
is inadequate for scattering problems of a complex shaped conducting object such as an
aircraft. At high frequency, edge diffractions, multiple reflections, creeping waves, and
surface travelling waves may also be important scattering mechanisms [32], [53].

118

Additionally, the spectral and angular windows for data are usually restricted by practical
constraints. Therefore, the time domain image of an aircrafi may be different from its
geometrical shape.

In this chapter, we will investigate time domain imaging of aircraft employing
SMEP responses, and interpret the reconstructed image from a new approach, based on
analysis of the scattering mechanisms and the back-projection algorithm utilized in image
retrieval. In section 6.2 we will interpret the scattering mechanisms of aircraft using the
ray method. The Explicit expression for equivalent currents, and the scattered fields from
those currents are derived in section 6.3. The time-domain inverse scattering identity with
the incorporation of Geometrical Theory of Diffraction (GTD) is derived in section 6.4.
In section 6.5, we will compare the images of an object consisting of two plates using the

PO approximation and GTD. Images of an aircraft model will be shown in Section 6.6.

6.2. Scattering Field Understanding Using Ray Method

Consider an aircraft illuminated by a source field from an initial reference surface
A. The field from a source point p’ at an observation point p can be tracked along a ray
that passes through p and originates at point p’. It is possible that several rays pass
through p as shown in Figure 6.1. In that event, the total field at p is synthesized by the
sum of the fields reaching p along the various rays.

The field along a pencil of rays can be calculated from a knowledge of the initial

field value and the ray geometry by [33]

 

. p‘I p; . e _jksr‘ (6.1)
(93 +s ‘)(p'2+s ‘)

 

u(F) = uA

119

 

 

Figure 6.1 Scattering by an aircraft model (composite obstacle).

120

where oil and p; are the radii of curvature of the initial wavefront, si is the distance

between the initial and the observation wavefront, and 14,, is the initial field.

For the incident field in fig. 6.1 the initial values on A can be calculated from a
known source distribution. The reflected field is caused by the incident rays (ray R2 in
Fig. 6.1) reflected at the object surface. The reflected laws for the fields on a curved
surface are the same as for plane waves on an infinite plane surface tangent to the curved
surface at the point of impact of the incident ray (the canonical problem). Thus, the
incident and reflection angles 9i are equal. The initial value of the reflection field at p,
is given by the incident field at p, multiplied by the plane wave reflection coefficient
R(9,-) descriptive of the reflecting properties of the boundary surface. The amplitude
variation along ray R2' involves the surface curvature at p ,, the curvature parameters for
the incident wavefront, and the angle coordinates specifying the directions of the reflected

and incident rays. Therefore, the reflected contribution to the field at p is given by

 

 

u, = “2R(et)\J , 9,192, lei”;
(p1+R2)(p2 +R2) (6.2)

an.
“2 = “.42 kR2

 

where u,” is the initial field for ray 2 on surface A.

The incident field along ray R3 striking the conical tip of the aircraft in Fig. 6.1
excites a spherical wave front and therefore a family of rays centered at the tip. The
canonical problem for the tip diffracted field is that of a plane wave incident on an
infinite conical obstacle. That solution provides the diffraction coefficient D", by which

the incident ray field is modified upon emerging from the conical tip. Thus, the

121

f

contribution at p due to the field along ray R3' is
’1'" ’3 e '1‘”:

m,

e
k1?I
3

(6.3)

 

 

u3, = D" u3, u3 = a“
where u ,4 3 is the initial field for ray R3 on surface A.
When an incident wave along ray R4 hits an edge of an aircraft wing, a family of
diffracted rays is excited. Those rays must obey the law of edge diffraction. A diffracted
ray and the corresponding incident ray make equal angles with the edge at the point of
diffraction, provided that they are both in the same medium. They lie on opposite sides
of the plane normal to the edge at the point of diffraction. When the two rays lie in
different media, the ratio of the Sines of the angles between the incident and diffracted
rays and the normal plane is the reciprocal of the ratio of the indices of refraction of the
two media. An edge-diffracted ray from a point P' to a point P is a curve which has
stationary optical length among all curves from P’ to P with one point on the edge

(Fermat’s principal for edge diffraction). The contribution at P due to the edge-diffracted

wave along ray R4' is [34]

 

7m
“4, = DdJ—L—u u — e— (6.4)

where D" is the diffraction coefficient, a,“ is the initial field for ray R4 on surface A, and
pI is the distance from the edge to the caustic of the diffracted rays.

Tip diffraction due to incident ray R3 also excites a diffracted surface ray (also
called a creeping ray) that travels along a geodesic on the shadowed surface and sheds

energy continuously. The launching amplitude L, of a surface ray field and its amplitude

122

variation M(R5) along its geodesic path R5 are determined from the canonical problem.

Therefore, the surface wave contribution to the field at p along ray R5’ is
’jk-Rs/

as, = u,DsM(R,)e""‘"5 L‘s—E (6.5)

where u, is given in (6.3), and R5 denotes the geodesic length along ray R5 from the tip

to the shedding point p,,.

The total ray-optical field at point p in fig. 6.] is now given by

u = ul +u2/+u3/+u4/ +u5/ (6'6)
where
e ’M'
“I = “Alfi— (6’7)

1

denotes the field along the direct ray from the source. The remaining contributions, given
by (6.2)-(6.5), are due to the presence of the scatterer.

The result in (6.6) contains only the dominant contribution from each of the ray
fields. Generally, each ray field has, in addition to this leading term, a series of higher-
order terms that decay inversely with k. Validity of the leading term alone implies that
the higher-order terms, in particular the second term, are small in comparison with the
first.

By the very construction of the field in (6.6), it is evident that the ray method
decomposes a complicated composite scattering problem into a sequence of simpler
canonical problem via the following steps:

(1) Determination of the incident field over an initial surface A.

123

(2) Determination of the reflected and diffracted ray fields that contribute at an
observation point p.

(3) Identification of canonical problems that treat separately each of the ray reflection
and diffraction problems. The solutions of these problems furnish the initial
amplitudes along various species of reflected and diffracted rays.

(4) Synthesis of a composite scattering problem by interaction (along rays) between

canonical constituents.

6.3. High Frequency Approximation Methods

Diffraction problems of the type schematized in fig. 6.1 can be formulated in
various ways. The second field generated by the induced surface currents can be
represented in terms of the direct radiation from the elementary currents distributed over
the aircraft surface A’.

In the illuminated region A ’ on a perfectly conducting object away from surface
singularities (such as edges, tips, or corners) and from strongly curved regions (where the

radius of curvature is not large compared to wavelength), one may approximate the

induced surface currents by their "physical optics" values. The physical optics currentsjp

are based on the local behavior of high-frequency fields described in the previous section
and are taken at any point p on A’ to have the same value as on an infinite perfectly

conducting plane tangent to A’ at p. Thus

—~ 0 FGA/ (6 8)
J = _. - °
pm {ZfixH '(f) FEA’

where fl" is the incident vector magnetic field. The secondary field generated by the

124

 

physical optics current gives the correct geometric-optical reflected field and the diffracted
field near the shadow boundaries. It does not provide, however, good results for the
diffracted field away from the shadow boundaries.

Even on a smoothly curved convex object with large radius of curvature, the

physical optics currents becomes invalid near the shadow boundaries. One may therefore

attempt to add to .7p a correction if, named the "fringe current" by Ufimtsev [35], so that

the sum equals the exact current value. The correction term is due to all other effects not
included in P0, such as edges. Because of the validity of the physical optics currents well
inside the illuminated region, and away from edges or corners, the fringe currents are
confined to the vicinity of those portions on A ’that lie near shadow boundaries or surface

singularities.

6.3.1 Equivalent Edge Currents for Arbitrary Aspect of Observation

The method of equivalent currents is a powerful technique in the analysis of the
scattering of electromagnetic fields. In addition to predicting the scattered fields in the
direction of a caustic, the method allows us to calculate the diffracted fields from an edge
of finite length and for observation angles away from the cone of diffracted rays. In the
following section, we will introduce Michaeli’s equivalent edge currents for arbitrary
aspects of observation since his formula is more general and rigorous than those of other

researchers [3 6].

According to the method of equivalent currents, the diffracted field E d due to an

edge discontinuity C is given for the Fresnel and far zones by the radiation integral [4]

125

I?" ~ jkf [21(765x(s*x£) +M(F’)§xf]G(r"’,F) dl (6.9)

where k is the wavenumber of the incident wave, Z=‘/u/e is the impedance of the

medium, a and n are the permittivity and permeability of the medium, respectively, F and F’

are the position vectors of the point of observation and a point on C, respectively,

d1 = [Jr/l is the increment of arc length I along C; i = d7’/dl is the tangent unit vector,
.6 = $75 = (F—f”)/ IF—f” | is the direction of observation of the radiating edge element at

F’ G(F’,f) = exp(-—jks)/4rts is the 3-D Green’s function, and Hr") = I(F’)f and

,

M0”) = M(f")i are the electric and ma etic e uivalent currents res ectivel .
811 Cl P y
Consider a plane electromagnetic wave propagating in the s’ direction and incident

on a perfectly conducting infinite wedge with an exterior angle Nrr, 0<N<2. Our aim is
to obtain explicit expressions for I and M at a point 0 on the edge for an arbitrary
direction of observation at. The geometry of the problem is depicted in Fig. 6.2. The point
0 in question is chosen as the origin. The x axis is directed along the normal to the edge,

and the y axis coincides with the external normal to the face. The 2 axis is directed along

the edge, so as to form a right-hand system. The angles between the s and s" directions
and the edge are
(6.10)

B = cos"(s‘-z‘), [3’ = cos’1(§’i)

The angle between the "upper" face and the edge-fixed plane of incidence

containing the vectors 5" and z“ is (V. The corresponding angle for the edge-fixed plane

126

 

 

 

 

Y
A I
S incident
\\ direction
\\ //,,.
A II '
S * +
X
observation
direction
2
Nit

 

Figure 6.2 Wedge scattering geometry.

127

of observation is d).

The expressions for I and M are given by [36]
M = H:- 2jZ(1/N) sin¢ sinKn-ao/Nl sin(N1r-d>) sinus-«am
zoksinfisfinfl’ 581““: 0061(n-a,)/N]-COS(¢’/N) Sinaz cos{(n-a2)/N] +COS(¢"/N)

, _ E‘ 021'Y(1/N)sin(d>’/N)[ 1 + 1
ksin’fi’ [60810: -a.)/N1 -COS(¢’/N) cos{(1r -a2)/N] +COS(¢’/N) (6.11)

 

 

_Hi 2j(1/N). ulcotB’-cotflcos¢ sin[(1r -a 1M]
.0 ksinB’ Sinai costar -a,)/NJ —cos(¢’/N)

 

_ uzcotB’ -cothos(N1t -¢) sin[(1t -a2)/N]
sinaz 008K?! -a2)/N1 +COS(¢//N)

where a, = cos“u,, u, = sinBcoscblsinB’, Y=1/Z, u, = sinBcos(N1t -d>)/sin[3’, and
at2 = cos‘luz. In another paper [37], Michaeli shows that I and M can be split into PO

components and fringe components

I = 1P0+lf, M = Mm+Mf (6.12)

where the PO components are the contribution from the surfaces, while the fringe
components are the contribution from the edge. The fringe components for an edge in a

half plane are [37]

lf=Ezi0 2jY [Shah/I2) /1_ /2 +Hz‘0 2] 1
oksinzfl’ cos¢’+u [I ut/2cos(4>/ )] zaksinfi’costb’w.

 

 

-l
°cotB’cosd>’ +cothos¢+ficos(¢’/2)(ucotfl’ -cot[3cosd>)(1-u) 2 (6°13)

M, = ”:2 2stin¢ 1 il_ ficosw/z)
Icsinpsinp’cos¢’+p[ l-p.

128

where

11 = [SinB’sinBcos¢+cosl3’cosB —cos"[3’]/sin2[3’ (6-14)

6.3.2 Scattered Field from Equivalent Fringe Currents of a Finite Edge
The far-zone radiation field maintained by the equivalent sources of an edge with

finite length L in the edge-fixed coordinate system is given by [38]

I
-e

E = -jo)/IT+"—kfxf H = (fxE’)/n (6.15)
E

where 11 = (nu/k and

3107 = iii-Waxes), Fir) iflffigem (6.16)

41tr 41tr

where

ae/(B’,¢’) = f [cosB’cos¢’Ix/(F’)+cosO’sin4u’Iy/(F’)-sin6’Iz/(f°’)]ej""F'Ids
= — [Zsin6’12,(F’)ejkf'F/dr’
a,,(e’,¢6 = f [-sin¢’1,,(r’)+cos¢’ly,(?’)]ei**"“’ds’ = o
(6.17)
ye,(6’,d>’) = f[case/cos¢’Mx/(F’)+cose’sin¢’My/(F’)-sin6’Mz/(F’)]ejkf‘f/ds’
= —f_msin6’Mz/(F’)e’”'F/dr’
axe/,4») = fJ-sinthx/(F’)+cos¢’My/(F’)]ej”'F/ds’ = 0

Choosing the local phase reference to be at the center of the wedge and assuming

129

the equivalent currents are constants over the length of the wedge, (6.17) reduces to

//:_m-/-ojkr'-F’l=_ -/U2j2k5cose’d
06(6 ,¢) [_msme lz(r’)e dr I(0)sme fm e g
= -I(0)Lsin6’sinc(kLcosO’) (6,18)
Q’s/(6%) = -fu‘fsine’Mz,(75ei*"'i’dr’ = -M(O)Lsin6’sinc(kLcos6’)

where sinc(x) = sin(x)/x. From (6.15)-(6.18), the backseattered field can be express as

i(F)— —ij+j—erF=1wuoLsmc(kLcosfi’)sm6’ebI(O)6I—M—(O)<]>I]
e

'1

 

4m
(6.19)

:2)
I'm

Hm: = JkLsmc(kLcose’)smB’ e -n—:r[l(0)d>I+ III—(1)61
1'I

= I

At a specific aspect, the contribution to the total backscattered field from the
wedge diffraction in time-domain is obtained by Fourier Transforming the frequency

response.

6.4 Time-Domain Inverse Scattering Identity with The Incorporation of Edge

Diffractions

The time-domain inverse scattering identity based on the Radon transform and the
space-time magnetic field integral equation has been derived for real-time use in short
pulse radar systems in chapter 3. In that formula, even though the correction term has
been added to the PO term, the contributions from the correction term is hard to evaluate
for a complex scatterer. Here we will investigate a new approach. Since the physical
optics currents work well inside the illuminated region, and away from edges or comers,

the fringe currents are confined to the vicinity of those surface singularities and shadow

130

 

boundaries. Because there are big diffractions from the edges of an aircraft, we only add
the edge fringe currents to the PO currents in the inverse scattering identity.

For an aircraft, there may be several edges which contribute to the diffraction
field. To a first order approximation, the difi‘raction field can be considered as a
summation of the contributions from each "visible" edge. Thus, For N "visible" edges, the

total diffracted field is given by

M.’_(_0) e,

:Iif(0)¢:+ 61' (6.20)
Ti

H d(f)= E jkL. sinc(kL. case I)sin6,I:

i=1

 

 

The SMEP diffraction response is obtain by multiplying the frequency response

(6.20) by the SMEP window 0(a)), and then inverse transforming into the time-domain

using the FFT as

h d(F,,t) = 7'(Q(w)17d(FS)-li i) (6.21)

where Ii I denotes the polarization direction of the incident magnetic field. The far-zone

backscattered field can be written as

h‘(r,t)= —fn *’ rah (’ ——T—)A(r‘, fi)ds’+hd(r,t) (6.22)
21tr

at

The first term of the right hand side in (6.22) represents the contribution from the PO
currents, while the second term of the right hand side denotes the contribution from the
diffraction fringe currents. Assuming the incident magnetic field is a SMEP, the same
procedure can be applied as in chapter 3 and the "thickness function" of the target in the

z-direction can be written as

131

 

2a2+2a—

E: [h ‘(iifln -h"(rs’.t)]

 

.. 2 2n 1»
m» = -—-f 1.]

+ (6,562)]; [h‘(f5’,t) —h d(5’,z)]d:}

 

(6.23)
dzdcp"
t-ZP +29’cos(<b’-¢‘)

 

This is the complete solution to the two-dimensional problem of recovering a body
from its scattered field; it requires the scattered fields and the diffraction terms from all
possible directions in the x-y plane. The next step is to solve equation (6.23). We can first

neglect the diffraction terms in (6.23). This is the time domain physical optics inverse

problem. We can thus obtain an initial estimate of 13(6), then find all the edges of the

image and use (6.20) to evaluate the total edge diffracted field h"(F,t). Finally, h"(f',r) is

added to (6.23) to obtain the modified I‘(p°).

6.5 Images of an Object Consisting of Conducting Plates

Consider an object consisting of two conducting plates placed on a rotating
pedestal as illustrated in Figure 6.3. Points P, and P2 are two vertices of one plate and
the line P,P3 forms an edge of the plate. In the laboratory coordinate system, define the
z-axis as the direction of the rotational axis, and the x-axis to be in the observation
direction. At the starting angle, the polar coordinates of the end points P, and P2 are (r,,
6,, 6),) and (r3, 9,, 4),), respectively. As the plate is rotated with an angle 4), the
coordinates of P, and P2 become (r,, 9,, ¢,+¢) and (r2, 0,, ¢2+¢), respectively.

Next we define an edge-fixed coordinate for edge P,P2 of the object. Let the z ’-

axis be in the direction of the edge P,P2, and the x ’-axis be normal to the edge and lying

132

EA) rotaa’t‘ignal

ps

 
  
 
    

P7
oonductin
plate :6 ’

 

 

 

/

X projection

 

antenna

Figure 6.3 Geometry of a conducting plate in the laboratory coordinate system and
edge-fixed coordinate system.

133

II
the

on the plate surface. The corresponding inclination angle of the transmitter/receiver to the
edge-fixed coordinate system is 6’ (¢)=O°). As the plane is rotated through an angle 4), the
corresponding inclination angle for the edge-fixed coordinate system becomes 9’(4)). It is
noted that 6’ is not only a function of 4), but also a function of the orientation of the plate
and the edge.

Assuming that the transmitting and the receiving antennas are inside the x’-z’
plane such that 9=n-9’ (B=1r-B’ in (6.13)), for the backseattered case. In our laboratory
system, let the polarization direction of the electrical field be along the y-axis, and the
polarization direction of the magnetic field be along the z-axis. In the following part, we
calculate the edge diffracted field and the PO field of the object and compare the images

of the object with and without the incorporation of edge diffractions.

6.5.1 Edge Diffracted Field

There are a total of eight vertices of the object which form eight edges of the
object. The equivalent edge electric and magnetic currents for each edge have different
forms depending an orientation and wedge angle. For edges P,P2, P2P” P,P5 and P5P,,

the equivalent currents will have similar formulas in their own edge-fixed coordinate

systems. Since 4) = 6’ = 180° and 11;, = o, from (6.13) we have

 

I{(0) = -E,',, ”2‘5 M{(0) = o i=1,2,3,4 (6.24)
ksin26,,/l-u,
where u,- = —1 —2cot26,. The incident field at the center 0 of the edge is
E‘m) = soc-1‘06, = Eo[i,cosfi,cos4),+i,cosfi,sin4),-z‘,sin6,] (625)

Then (6.24) becomes

134

1.

. - 25
I,’(0)— — E‘, jrzf = Eosir16. ”2‘5 = J——" (6.26)

oksinzfl ,/1- )1 'Icsinze,,/2+2cotze, ‘0“

Substituting (6.26) into (6.19) yields the total diffraction fields from these four edges

                 

                

155%) = 452::
"’ (6.27)
HR?) = 1,2;
'1 i=1,-i

For edges P,P7 and P3P,,, since 6 = 6’ = 90", u = c054), and the equivalent current can

be written as

[f(0) = O Mf(0)" on Lint-1 (6.28)
k cosd)

Thus, the total diffraction fields from edges P6P, and P3P, are (in edge-fixed coordinates)

e-flve sin -1 fir.- sin .-1
¢ £fin=—£LK I
c054), ,-5 4m c054),

(6.29)

 

 

36W $1,;

For edge P,P3, ¢=41’=0 or 90". From (6.11) we have I(O)=O, M(0)=O, so there is no

 

2
contribution from edge P6P, For edge P,P,, ¢=90°, u=-2 cffizg . From (6.13) we have
sm

I «0) on_Y \/1+cos 6-sin6 M «0) = O

The diffracted field from edge P7P, is

 

 

 

‘17”: ,/1+c0826 -sin6 E 8 6.31
E:(r‘) = 5011; e sinc(kL,,cosBs)sin68 8 8 ”43(7) = 30"“ )
r8 cos268 '1

135

The total diffraction fields from this object can be considered as a summation of the

contributions from each "visible" edge. Thus, the total diffracted field is given by

e,-jkr e'l’” sin4).— —1
E6d(r‘)— - —2E °§i —n;Lll,sin6§inc(kLcos6)+ 025: LI41t—T cosd)

 

(6.32)

e "’"s . , l+c0526 -sin68
+ 501': smc(kL8c056,,)s1n68
41: r8

 

 

cos268

Equation (6.32) can be used to obtain the frequency responses over the finite
bandwidth 4-16GHz, multiplied by a SMEP window, and then inverse transformed into
the time-domain using the FFT to provided a SMEP response. Figure 6.4 shows the
diffracted SMEP response of the object at 4)=45° using the SMEP shown in Fig 4.4. The

response shows that there are two peaks located at range about :(L/Z) c056, for each edge,

which are at the differential ranges of the end points of the edge. At the rotation angle
4) such that 6i = 90°, as shown in Fig. 6.5, the SMEP response gives a single peak with
strong amplitude for each visible edge because the two end points of the edge have the
same differential ranges and all the points on the edges are equidistant to the observation

point.

6.5.2 Scattering Field Under Physical Optics Approximations

Consider the same object illuminating by a plane wave defined by

ion -.- éEoe 1'5"? fi‘m = k ere-ji’r‘ (6.33)

 

The far-zone scattered field can be written as

136

—jkr

 

53(7) = jI—k" f'x(fx&’) (6-34)
‘11:

where d(f) = fj(f°’)e‘j"f'FId.s’ is the directional weighting function. Using the PO

approximation, the induced current can be written as (in the illuminating region)

f0) = 2fi(r’)xfi’(7’) (6-35)
Rewrite the scattered field as
a E -jkr _.. -,_,
E ’(f) = éj—°—" f fx(fx[fix(k'xé)])e‘fl* my (6.36)
p 21: r s
For the backscatter case, f = -k', and
fx(r‘x[fix(l?’><é)]) = —I€’x[(IE"xé)(fi-IE°’)] = “(a-”5 (6-37)
so
_. E 7*, _.i ~,' _./ . E 'jb . .
556*) = ej—°"— ri’-k e72" "615’ = -6j—3e—fkc036e’2’“s‘“°ds
p 21: r s 21: r s
(6.38)

 

— éj kEOL5L8 e 7*,

cosBsinc(kL8sin6)
2n r

where L 5 and L, are the lengths of edge P6P, and edge P,P,.
Fig. 6.6 and Fig. 6.7 show the numerical SMEP responses of the object at two

different angles. The response shows that there are two peaks located at range about

:(Lj2)c056, for each horizontal edge which are at the differential ranges of the end points

of the edge. At the rotation angle 4) = 0°, as shown in Fig. 6.6, the SMEP response gives

a single peak with strong amplitude for the vertical plate because all the points on the

137

plate have the same differential ranges and are equidistant to the observation point.

6.5.3 Images of The Plates
The total scattered fields of the plates can be considered as the combination of the

edge diffracted field and the PO scattered field. They can be written as

6‘6) = E‘flmidm (6.39)

where 5pm denotes the PO scattered field and 5.40) denotes the total edge diffracted

field. Substituting (6.39) into (6.23) without subtracting the edge diffracted field we may
create the image of the plates. This is usually the practical case since the measured field
is usually the total field, and the edge diffracted field is very difficult to obtain for a
complicated target. Figure 6.8 gives the image of the plates shown in fig. 6.3. We can see
they match very well. Because equation (6.23) is derived under the PO approximation,
we may use the PO scattered field in (6.23). This gives the true image of the plates under
the PO approximation. The image is given in Fig 6.9. Comparing with Fig 6.8, we can
see that the PO approximation is inadequate for an object consisting of conducting plates
where the edge diffractions may be the dominant contributors to the scatted field. Thus,
the image of a metallic object might be different from its geometrical shape using the PO

approximation.

6.6 Images of an Aircraft Model
A simulation of time-domain imaging is can'ied out using data measured in the
Michigan State University free-field scattering range. We use the same measured data as

used in chapter 4. Figure 6.10 show the images of a 1:72 scale model B-52 aircraft. Since

138

 

 

we use the PO approximation in (6.23), and did not subtract the diffraction contribution,
the edges of the aircraft are clearly visible, and the surface contribution nearly neglectable
compared to the edge diffraction. This shows that edge diffractions are dominant
contributors to the scattered field when the scatterers consist of conducting plates. Next,
we can obtain the edges from Fig. 6.10 and use (6.20) to evaluated the total edge
diffracted field from all the "visible edges". Fig 6.11 shows the image of the main edges
of the aircraft model. Subtracting the edge diffraction field from the total measured
scattered field yields the PO component and the corresponding image is shown in Fig.
6.12. We can see that the wing edges are not as visible as in Fig. 6.10. Theoretically, the
wings should vanish after subtracting the edge diffracted fields, but since the scattered
field oscillates at high frequency, and the diffraction field from surface singularities such
as tips, corners, and multiple reflections, creeping waves and surface travelling waves are
not taken into account and they may also be big contributors to the scattered fields. Since

they may also distort the reconstructed image, the results are not as good as expected.

6.7 Conclusions

In this chapter, we have shown that in the high-frequency region, the scattered
field can be attributed to a combination of different canonical problems. This idea has
been used in developing a new time-domain inverse scattering identity and interpreting
the reconstructed image from a new approach, based on analysis of the scattering
mechanisms and the backprojection algorithm utilized in image retrieval. In this approach,
we added the edge diffracted field to the physical optics field so that the total is equal to
the exact scattered field value. We also showed that for those objects consisting of

conducting plates, edge diffractions are dominant contributors to the scattered field. The

139

 

’6

 

total edge diffracted field is the summation of the diffracted field from each "visible"
edge. Several numerical and experimental examples have been included to validate the

formulas developed in this chapter.

 

140

 

0.08

 

0.04

0.00

Magnhude

—0.04

11111111111111111111111111111l1111111111

 

 

—0.0B IIIIIIIII'IITTTITIIIIIITIIIIIIIIIlITjTI

4.50 6.50 8.50 10.50 1250
Time (ns)

Figure 6.4 Diffracted SMEP response of the object with rotational angle 4)=45°.

141

 

1.00

 

 

 

 

 

0.50}
q
C» i
'8 4
,J -
.E 0.00-l ‘W’A
U) I
O ..
E .4
—0.50§
3
-1.00 'IIITITFTT—IMTIIIIIIIIIITTWIIIIITIIIWTTTfl
4.50 6.50 8.50 10.50 12.50

Time (ns)

Figure 6.5 Diffracted SMEP response of the object with rotational angle 4)=O°.

142

 

1.50

0.50

 

Magnhude
s's
tn
CD

—1.50

-2.50

 

11111111111111111111111111111111111111111L1L11111111111

 

—3.50 111rrIfTrrtrt—IrlrrrthrIlIIrr]rtrIIrTIr]

4.50 6.50 8.50 10.50 12.50
Time (ns)

 

Figure 6.6 SMEP response of the object with rotational angle 4)=0°
under PO approximation.

143

V-

0.02

0.01

Magnfiude
O
O
O

-0.01

11111111111111111111111111111111114111|

 

—0.02 IIIIIITFTITWrTT—erllrrllrtrrrrtrrrrrT]

4.50 6.50 8.50 10.50 12.50
Time (ns)

 

Figure 6.7 SMEP response of the object with rotational angle 4)=45°
under PO approximations.

144

 

.‘10 -

.so- *“——————————————" i I J H I '-

 

02° '1

0101

.0101

.020 It

 

°o3° <

 

 

Figure 6.8 Image of the object from O°-180° data in band 4-16GH2
using the total scattered field. ‘

145

 

 

‘1.

.30 <

I20 1

 

.10-

 

.04

.10:

.20-

.80 -

 

...I

 

Figure 6.9. Image of the object from 0°-180° data in band 4-16 GHz under PO
approximation.

146

 

 

.‘10-

03° '1

 

.10-

.10 ‘ J.

I20 4

.30 -

 

 

.Ho-l
'.50 3‘10 '.30 '.20 '.10 .00 .10 .20 .30 .‘00 .50

 

Figure 6.10. Image of B-52 from 0°-180o data in band 4-16 GHz using measured data.

147

 

.140. .
.0

03° 1

02° 1

0101

00*

.0‘0 ‘

 

'.20 a

.030 ‘

w... 0*"... M

 

'.50 ‘.‘10 .030 .020 ‘.10 .00 0‘0 02° 03° 0‘10 .50

Figure 6.11. Image of main edges of B-52 from O°-180° data in band 4-16 6112 using
theoretical edge diffraction field

148

 

'\

 

.H0 -

03° '1

 

.20 s

u . l
I. .

0

.10. .- .

 

 

-.10.

'.20 '1

’.80-l

 

'.‘IO ‘

 

 

”.50 ‘.'-IO '.30 “.20 .‘10 .50

Figure 6.12. Image of B-52 from 0°-180° data in band 4-16 GHz with the incorporation
of edge diffraction.

149

 

CHAPTER 7

CONCLUSIONS

This chapter summarizes what is achieved in this thesis and points out what

remains to be studied further.

7.1 Summary of The Thesis

The most significant contribution of this thesis, from the view of the writer, is the
development and the implementation of a time-domain radar imaging technique. In this
technique, starting with the exact space-time magnetic field integral equation, and by
using Radon’s theory and a SMEP as the incident pulse, we obtain the complete inverse
scattering identity which considers both illuminated and shadowed range contributions.
We have developed the inverse scattering identities for both monostatic and bistatic cases,
for special cases such as rotationally symmetric targets, for flying object, and for the 3-D
case. The entire derivation is carried out directly in the time domain, and the incident
magnetic field waveform is general. We can see from equation (3.22) that when the
incident magnetic field is an impulse, step, or ramp, previous inverse scattering identities
can be obtained. The reason why we chose the SMEP is its band-limited property and its
ease of synthesis in the frequency domain.

In chapter 3 we have shown that the size and shape of an object (thickness

150

 

 

 

function) can be obtained from its cross sectional area functions, and the problem can be
reduced to the classical Radon problem. Since the cross sectional area function can be
estimated from the scattered field when the incident field is a SMEP, we can reconstruct
the thickness function of the object using backprojection techniques which have been
widely used in computed tomography, radio-astronomy, and geophysics. Chapter 4 gives
two reconstruction algorithms. One is a time-domain approach called convolution
backprojection. Another is the frequency-domain approach called filtered backprojection.
Both approaches provide similar quality reconstructions. By using synthesized SMEP
response data and the reconstruction algorithm, we have obtained very good images of
several aircraft models. A simulation based on stepped-frequency, multi-aspect
measurements of aircraft models produces clear images with highly-defined edges.
Chapter 5 deals with the limited-view problem. We have proposed two techniques
to handle this practical situation. One approach is the method of projections onto convex
sets (POCS). The basic principle of POCS is that each piece of a priori knowledge must
be represented by a convex set onto which the current image estimate can be projected.
It has been shown that if the intersection of these convex sets is nonempty then the
sequence of cyclic projections will converge weakly to a point in this intersection. We
extend this approach to radar imaging for the first time and show some useful results.
Another approach which we have demonstrated is to process the available measured
projections in order to generate an estimate of the full set of projections, an image which
is called a sinogram. The goal of this approach is to recover the sinogram from the
available measured data using linear prediction. Since the scattered field of a target can
be written as a superposition of distinct specular reflections arising from scattering centers
on the target, the position and strength of the scattering centers can be predicted using

151

 

' 1

linear prediction with the change of the observation angle. Thus, the missing data can be
predicted before reconstructing the image. Big improvements in image reconstruction have
been achieved using this technique, and some useful results have been shown in chapter
5.

In chapter 6 we point out that the PO approximation is inadequate for scattering
problems of a complex shaped conducting object such as an aircraft. At high frequency,
edge diffractions, multiple reflections, creeping waves, and surface travelling waves may
also be important scattering mechanisms. Additionally, the spectral and angular windows
for data are usually restricted by practical constraints. Therefore, the time domain image
of an aircraft may be different from its geometrical shape. In chapter 6 we have
investigated time domain imaging of aircraft employing SMEP responses, and interpret
the reconstructed image from a new approach, based on analysis of the scattering
mechanisms and the back-projection algorithm utilized in image retrieval. The time-
domain inverse scattering identity with the incorporation of Geometrical Theory of

Diffraction (GTD) is derived and some interesting experimental results are provided.

7.2 Suggestions for Further Studies
From the view of the writer, the following points are the suggestions for future

research:

1). At present the reconstruction has been limited to the 2-D case, so the next logical
step is to extend the reconstruction algorithms to a 3-D object. The time-domain
inverse scattering identity for 3-D cases has been derived in chapter 3 using a half
Gaussian pulse modulated sine as an incident field. The well-developed
reconstruction algorithms in other fields for the 3-D case can be extended and

152

 

 

2).

3).

applied to the 3-D inverse scattering problem.

In addition to the SMEP, impulse, step, and ramp, other types of incident field
waveforms may be used to obtain the cross-sectional area function.

In sinogram restoration for the limited-view problem, l-D linear prediction has
been used to fill in the missing data. It is more efficient and may be more

accurate if we use 2-D linear prediction, since the sinogram is a 2-D map.

153

 

APPENDIX A

 

APPENDIX A

ANECHOIC CHAMBER MEASUREMENTS

A.1 Introduction

There are two different schemes for performing transient measurements. The first
method uses traditional time domain reflectometry methods [39], [40], [41], [42] while
the second one utilizes a Fourier synthesis approach [43], [44], [45]. The time domain
reflectometry method has the important advantages of time-gating and rapid speed of
measurement, but suffers the disadvantage of limited dynamic range and equipment
instability which leads to timing jitter. In this dissertation, we do not exploit this method.

A vector network analyzer is capable of measuring both the magnitude and phase
of the S-parameters of two part network over a very wide frequency band. This capability
gives the vector network analyzer the ability to synthesize the transient response of a
target via the inverse discrete Fourier transform. Measurements using this configuration
take longer than with time domain systems. Note that the sampling interval constraint
limits the use of the network analyzer system when the antennas are used on a ground
plane since the reflections from the edges of the ground screen and the laboratory ceiling
cause alaising.

A typical free-field dual antenna frequency-domain measurement system is shown
in Figure A. l. A measurement of 52, produces a measure of scattering from part 1 to port

2 from both the scatterer and the environment. Amplifiers are added to the transmit

154

 

 

 

A

 

@

Transmit

Target

Q 54

l/ \l
V
Foam Pedestal

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.3- Anechoic Chamber
>
3
g
e w
_I"_ 821 VGCIOI'
Network
I Analyzer
Port 1 P0“ 2
1
4

Figure A.1 Dual antenna free-field frequency-domain transient measurement system.

155

channel to help boost the dynamic range of the measurements. With the addition of wide-
band hom antennas, the frequency range of the system has been extended from the band
1-7 GHz to the band 2-18 GHz. Also, a computer-controlled rotator allows measurements
to be made automatically at a 0.150 aspect angle increment.

This Appendix reviews the calibration and deconvolution method of the
measurement system developed by Ross [45]. Examples of measured spectral responses
of calibration targets (i.e. sphere) are presented and compared with theory. Various
methods which improve the measurements of aircraft models are described briefly and

summarized.

A.2 Calibration Procedure

Ross made a significant contribution to automatic measurements in the anechoic
chamber at the MSU EM lab. The size of chamber is 24 ft in length, 12 ft in width, and
12 ft in height. Here we will review the calibration procedure. To take full advantage of
the wide frequency sweep available with a network analyzer, the clutter and system
transfer function must be removed from the measurement. That is, all repeatable
systematic errors associated with anechoic chamber clutter, antennas, and transmission
lines should be removed from the raw measurement.

First, the background response measured when the chamber is empty is subtracted
from the raw target measurement. Figure A.2 shows the spectral response of a 14 inch
metal calibration sphere measured at 2-18 GHz. Other possible calibration targets include
a thin wire. One of the difficulties with the wire is that the wire response is strongly
aspect-dependent. It is very hard to determine the incident angle exactly. Afier the
spectral response is zero-padded and transformed to the time domain, any interaction

156

 

 

 

terms and noise that are sufficiently delayed beyond the end of the calibration target
response can be eliminated. The spectral response obtained by transforming the time-
domain response using the FFT becomes much smoother than before as shown in Figure
A.3. Next the system function can be found by dividing this response by the theoretical
value.

The theoretical scattered field response of a sphere was first solved by Mie in
1908 [46]. It has also appeared in many books [47], [48]. The notation used in the course
notes by Chen [49] is used here.

Consider a plane wave incident on a perfectly conducting sphere as illustrated in

Figure A.4. The incident electric field is given by
5"”(7) = JEEoexp(-jkz) (AJ)

The far zone scattered electric field is obtained by approximating the spherical vector

wave functions in Mie series [48]

Eg=-E

Pl
q)e___xp(—jkk)z 2n+___1_{an”»_(°°se’ .bn[%pnl(eme)]}em¢ (A.2)

,,_ , n(n +1) sin6

_e____xp(-ij) 2n+1P1(c056)+ 3 1 - (A.3)
ES, =jE0 "Z; n—(—n+l){b sin6 a"[66p" (0036)] 31nd)

where a, and b,, are found using the boundary conditions and are given by

a .
. —(Rjn(kR))|R=a
- Mk“) b = 5R (A.4)

_ 2 n
hi ’(ka) a—imhflkR» lR=a

 

For the back-scatter case, the following relations are used to the evaluation of the Mie

series in equation (A.2) and (A3)

157

 

 

Pn’(c056) I,” = -(‘21)" n(n+1) (A.5)

£13,1(e6s0) |,,:, = - “21)" n(n+ 1) (A6)

 

The expressions shown in (A.2) and (A3) are programmed with the use of widely
available generalized Bessel function subroutines. These subroutines are very accurate and
efficient.

Figure A.5 shows the theoretical waveform obtained using the above Mie series.
The system transfer function obtained using this theoretical response is shown in Figure
A.6. To provide verification of the calibration process, a second sphere 3 inches in
diameter is used. Figure A.7 shows a comparison of the measured and theoretical spectral
responses of the smaller sphere. The agreement is very good except in the high frequency
region. This is probably due to the hole in the sphere and effects of windowing of

extraneous system noise from the measured response.

158

”r

 

0.010

 

0.008
a)
‘0
3
.1; 0.006
a.
E
o
a; 0.004
.5
2
<1)
0:
0.002
0.000 I I I I r I I I I I I I I l I I I I I I l I I I I I I I I I l I 1 I I I IW
. 5.0 10 0 1 0 20.0

. 5.
Frequency (GHz)

Figure A.2 Spectral response of 14 inch metal calibration sphere, measured using
frequency domain system in MSU anechoic chamber.

159

 

0.008

 

 

 

 

:
.4
0.0065
0) " ‘
'0 I
3 -
fl: .1
EL 1
.1
Q) d
.2 ‘
+4 I
2 .
‘D I
CE0.002«~
I
1
0.000 IIWIIIIITIIIIIIIIIIIIIIIIII1IIIIIIIITIfl
0.0 5.0 10.0 15.0 20.0

Frequency (GHZ)

Figure A.3 Spectral response of 14 inch metal calibration sphere afier time-domain
filtering. Measured using frequency domain system in MSU anechoic
chamber.

160

 

 

 

 

 

 

Figure A.4 Geometry for sphere scattering problem.

161

0.036

 

 

a) -
“00.032-
3 .4
r: q
a -4
E I
0 -
Q) ..
> -I
1;; q
20.028-
(D -
0: .
0.024 lIllIll’lllIlIllllllIllllIlllIllllllI—Yfij—l
O 20.0

Frequency (GHz) '

Figure A.5 Theoretical response of 14 inch metal calibration sphere calculated using
Mie series for MSU anechoic chamber configuration.

162

0.250

 

 

 

 

0.2005
G» E
‘D ..
D :
E0150:
0- 1
E :
O _.
20.100;
2.: .
O I
a) :
0i :
0050-:
i
0.000 rIIIIIIIIIIIIIIIIIIIIIIIIIIII'IIIIIIII
0.0 5.0 10.0 15.0 20.0

Frequency (GHz)

Figure A.6 Transfer function of frequency domain system obtained using 14 inch
diameter sphere as calibration standard.

163

 

0.014

0.012

0.010

Mognhude
. .0
§

0
O
O
O}

 

can 1 '

0.002

 

-— theoretical 3 inch sphere
------ measured 3 inch sphere
i

IIIIITT]
2.0 6.0

0.000

IIIWIIIIIIrIIIIIIIIIIIIj
10.0

110 1&0
Frequency (GHZ)

Figure A.7 Comparison of theoretical and experimental frequency responses of a 3

inch diameter sphere. Measurements performed using frequency domain
system in MSU anechoic chamber.

164

APPENDIX B

APPENDIX B

EXPERIMENTAL EQUIPMENT

This appendix summarizes the apparatus used for the experimental portion of the
research. The equipment used in the frequency domain scattering measurements is listed
in Table B.1. There have been three different configurations used at EM lab in MSU so
far. The "low-band" configuration uses a the PPL-5812 amplifier to amplify the signal
from port 1 of the network analyzer. The parameters used in this configuration are shown
in Table 8.2. The "high-band" configuration uses the HP-8349B amplifier to amplify the
signal from port 1 of the network analyzer. The parameters used in this configuration are
described in Table 8.3. The enhanced "high-band" configuration uses the HP-8349B
amplifier too. The parameters used in this configuration are given in Table B.4.

To accurately characterize the measurement system, the network analyzer was used
to measure the gain of the amplifiers and attenuations of the cables over a wide range of
frequencies. The results of these measurements are given in the thesis by Ross [45] and
are reproduced here to indicate the limitation of measurement system, the improvements,
and baseline measurements to assess future equipment changes. The signal gain of the
HP-8349B amplifier is shown in Figure 8.1. Note that the gain is about 20 dB over 2
GHz. The gain of the PPL-5812 amplifier is shown in Figure 3.2. The gain decreases a

lot over 4 GHz. This attenuation

165

 

Table 8.1 Equipment used for frequency domain scattering measurements.

Hewlett-Packard HP-87ZOB vector network analyzer

Hewlett-packard HP-8349B microwave amplifier

Picosecond pulse labs PPL-58l2 10 dB broadband amplifier

American electronics laboratory AEL H-1734 wideband TEM horn antennas
23.5 foot RG-9B cable with N-type connector

22.5 foot RG-9B cable with N-type connector

B&K precision DC. power supply 1610 (for PPL amplifier)

Insulated wire inc. AEL H-1498 wideband TEM horn antennas

12 foot low loss cable with SMA connector

Table 8.2 Low band measurement parameters.

Measurement parameter: 821
Frequency sweep: 0.4 - 4.4 GHz
number of points: 401

IF bandwidth: 100 Hz
number of averages: 10

Sweep time: 30 seconds

 

limits usage of amplifier in higher frequency range. Cable losses for the system are shown

in Figure 8.3. Note that this cable is very lossy at high frequencies. Now a pair of low

166

 

loss flexible cables 12 foot in length are available to help improve range performance.
Tests of new cables show a 3.5 dB loss at 20 GHz.

Table B.3 High band measurement parameters.

Measurement parameter: 821
Frequency sweep: 1.0 - 7.0 GHz
number of points: 601

IF bandwidth: 100 Hz
number of averages: 10

Sweep time: 60 seconds

Table B.4 Enhanced high band measurement parameters.

Measurement parameter: 821

Frequency sweep: 2.0 - 18.0 GHz
number of points: 1601

IF bandwidth: 100 Hz
number of averages: 2

Sweep time: 60 seconds

167

”Om

Gain (dB)

1 0.00

 

0.00

 

-‘OOWIIIIIIIIIIIIIIIIUII
0.00 5.00

Illrrlrlrrrltulllirtfi

10.00 15.00 20.00
Frequency (GHz)

Figure 8.1 Gain of HP8349B Microwave Amplifier measured using HP-87ZOB
network analyzer. Amplifier input power = -30 dB,,n

168

 

1 5.00

 

 

H
d
I
I
I
d
10.00-
.1
d
d
q
I
5.”-
A d
m d
In ..
v ‘
q
.E 4
D d
o d
0.00-3'
q
d
d
d
d
d
-5.00-
d
.1
I
d
d
.4
-IOom IITIIIIII'IIIIIIIIIIIUTIIUUIUIIUITIIIIIIIIIUIIIIIIrTTThIIr]
0.00 1 .00 2.” 3.00 4.00 5.00 6.”

Frequency (CH2)

Figure 8.2 Gain of PPL 5812 wide band amplifier measured using HP-8720B network
analyzer. Amplifier input power = -10 dBm.

169

QI-h—I

l
3
3

I
8
8

-70.00

i:
s

0.00 2.50 5.00 7.50 10.00 12.50 15.00

lllllllllllllLlllllIlllLlllll

    

 

Cable Attenuation

Transmit Cable
------ Receive Cable

 

 

Figure B.3

  

IIIIIIIITI1TIIIIIIIrIIIIIIIIIIIIIIIIIIII

17.50 20.00
Frequency (GHz)

Measured cable attenuation.

170

BIBLIOGRAPHY

 

10.

11.

12.

BIBLIOGRAPHY

C. L. Bennett and G. F. Ross, "Time-Domain Electromagnetics and Its
Applications" Proceedings of IEEE, vol. 66, No. 3, pp. 299-318 March 1978.

L. B. F elsen, Transient Electromagnetic Fields, Springer-Verlag, New York, 1976.

R. Mittra, Computer Techniques for Electromagnetics, Pergamon Press, New
York, 1973.

B. P. Rynne, "Instability in Time Marching Method for Scattering Problems",
Electromagnetics (1986) 6, pp. 129-144.

F. D. Smith, "Instabilities in Time Marching Methods for Scattering: Cause and
Rectification", Electromagnetics (1990) 10, pp. 439-451.

E. M. Kennaugh and R. L. Cosgriff, "The use of impulse response in
electromagnetic scattering problem," REE National Convention Record, part 1. PP
72-77, 1958.

E. M. Kennaugh and D. L. Moffatt, "Transient and impulse response
approximations", Proc. IEEE, vol. 53, pp.893-901, Aug. 1965.

J. D. Young, "Radar imaging from ramp response signatures", IEEE Trans.
Antennas Propagat. vol AP~24 no. 3, pp.276-282, May 1976.

C. L. Bennett, "Time domain inverse scattering", IEEE Trans. Antennas propagat.
vol. AP-29, no. 2, pp.213-219, Mar. 1981.

E. J. Rothwell, K. M. Chen, D. P. Nyquist and J. E. Ross, "Time domain imaging
of airborne targets using ultra-wideband or short-pulse radar", IEEE Trans.
Antennas propagat. vol. AP-43, no. 3, pp.327-329, Mar. 1995.

D. Ludwig, "The Radon transform an Euclidean spaces", Comm. Pure and Applied
Math., v01. XIX, pp.49-81, 1966.

Y. Das and W. M. Boemer, "On radar target shape estimation using algorithm for

reconstruction from projections", IEEE Trans. Antennas Propagat., vol AP-26, no.
2, pp.274-279, Feb. 1978.

171

l3.

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

S. R. Deans, The Radon transform and some of its applications. John Wiley &
Sons, 1983.

N. N. Bojarski, "A survey of the physical optics inverse scattering identity", IEEE
Trans. Antennas Propagation, vol. AP-30, No. 5, September 1982.

Y. Dai, E.J. Rothwell, D.P. Nyquist and KM. Chen, "Time-domain imaging of
radar targets using ultra-wideband or short pulse radars", IEEE AP-S International
Symposium and URSI Radio Science Meeting. 1995.

J. Radon, "On the Determination of Function from their Integrals Along certain
Manifolds", Ber. Saechs. Akad. Wiss. Leipzig, Math. Physics Kl., vol. 69, pp. 262-
277, 1917.

G. N. Ramachandran and A. V. lakshminarayanan, "Three dimensional
reconstruction from radiographic and electron rnicrographic application of
convolutions instead of Fourier transforms," Proc. Nat.Acad. Sci. US, vol.68, pp.
21-24,l974.

A. K. Louis, "Picture reconstruction from projections in restricted range," Math.
Meth. Appl. Sci, vol. 2, pp.209-220, 1980.

T. Inoye, "Image reconstruction eith limited angle projection data," IEEE Trans.
Nucl. Sci, NS-26, pp.2666-2669, 1979.

J. A. Reeds and L. A. Shepp, "Limited angle reconstruction in tomography via
squashing," IEEE Trans. Med. Imaging, vol. MI-6, pp.89-97.

M. I. Sezan and H. Stark, "Tomographic Image Reconstruction from Incomplete
View Data by Convex Projections and Direct Fourier Inversion", IEEE Trans.
Medical Imaging, vol. MI-3, No. 2, June, 1984.

L. M. Bregman, "Finding the common point of convex sets by the method of
successive projections," Doki, Akad. Nauk USSR, vol. 162, no. 3, pp. 487-490,

1965.

L. G. Gubin, B. T. Polyak, and E. V. Raik, "The method of projections for finding
the common point of convex sets," USSR Comput. Math. Math. Phys, vol. 7, no.
6, pp. 1-24, 1967.

D. C. Youla and H. Webb, "Image reconstruction by the method of projections
onto convex sets-Part 1," IEEE TRANS. Med. Imaging, vol. MI-l, pp 81-94, Oct.
1982.

K. M. Hanson G. W. Wecksung, "Bayesian approach to limited-angle

172

26.

27.

28.

29.

30.

31.

32.

33.

34.

35.

36.

37.

38.

reconstruction in computed tomography," Appl. Optics. vol. 24, pp. 4028-4039,
Dec. 1980.

P. Oskoui-Fard and H. Stark, "Tomographic Image Reconstruction Using the
Theory of Convex Projections", IEEE Trans. Medical Imaging, vol. 7, No. 1,
March,1988.

J. L. Prince and A. S. Willsky, "Hierarchical Reconstruction Using Geometry and
Sinogram Restoration", IEEE Trans. Image Processing, vol. 2, No. 3, July, 1993.

G. T. Herman, Image Reconstruction from Projections: The Fundamentals of
Computerized Tomography, New York: Acdemic, 1980.

J. Makhoul, "Linear Prediction: A Tutorial Review", Proceedings of The IEEE,
vol. 63, No. 4, April, 1975.

J. G. Proakis and D. G. Manolakis, Digital Signal Processing, Principle,
Algorithm, and Applications. Macmillan, 1992.

Y. Dai, E.J. Rothwell, D.P. Nyquist and KM. Chen, "Time-domain imaging of
radar targets using algorithm for reconstruction from projections radars", IEEE
Trans. Antennas Propergat. Accepted for publication, 1997.

G. T. Ruck, D. E. barrick, W. D. Stuart and Krichbaum, Radar Cross Section
Handbook. New York: Plenum.

G L. James, Geometrical Theory of Diflraction for Electromagnetic Waves, Third
Edition, Herts SGl lHQ, United Kingdom, 1986.

J. B. Keller, "Diffraction by an Aperture", J. Appl. Physics, vol. 28. pp.426-444,
Apr. 1957.

P. Ya. Ufimtsev, "Approximate Computation of Diffraction of Plane
Electromagnetic Waves at certain metal Objects - I. Diffraction Patterns at a
wedge and Ribbon" Sov. Phys. - Tech. Phys. 27, pp. 1708-1718, 1957.

A. Michaeli, "Equivalent Edge Currents for Arbitrary Aspects of Observation",
IEEE Trans. Antennas Propagate, vol. AP-32, pp. 252-258, Mar. 1984.

A. Michaeli, "Elimination of Infmities in Equivalent Edge Currents, Part I: Fringe
Current Components", IEEE Trans. Antennas Propagate, vol. AP-34, pp. 912-918,
July- 1986.

D. P. Nyquist, "Antenna Theory", EE929A course notes, Michigan State
University.

173

39.

40.

41.

42.

43.

44.

45.

46.

47.

48.

49.

50.

51.

52.

53.

M. A. Morgan, "Time-domain scattering measurements," IEEE Trans. Antennas
Propagat. Newsletter, vol. 6, pp. 5-9, August 1984.

M. A. Morgan and B. W. McDaniel, "Transient electromagnetic scattering: Data
acquisition & signal processing," IEEE Trans. on Instrumentation and
Measurement, vol. 37, pp. 263-267, June 1988.

M. A. Morgan and N. J. Walsh, "Ultra-wide-band transient electromagnetic
scattering laboratory," IEEE Trans. Antennas Propagat., vol. 39, pp. 1230-1234,
August 1991.

M. A. Morgan, "Ultra-wideband impulse scattering measurements," IEEE Trans.
Antennas Propagat., vol. 42, pp. 840-846, June 1994.

F. T. Ulaby, T. F. Haddock, J. East and M. Whitt, "A millimeter-wave network
analyzer based scatterometer," IEEE Trans. on Geoscience and Remote Sensing,
Vol. GE-26, pp. 75-81, Jan. 1988.

B. Drachman, "Two methods to deconvolve: L'-method using simple algorithm
and Lz-method using least squares and a parameter," IEEE Trans. Antennas
Propagat., vol. 32, pp. 219-225, March 1984.

J. E. Ross, "Application of Transient Electromagnetic Fields to Radar Target
Discrimination," Ph. D. dissertation, Michigan State University, 1992.

Mie, Ann. Physik, 25, 377, 1908.

J. A. Stratton, Electromagnetic theory, New York and London: McGraw—Hill,
1941, pp. 563-573.

R. F. Harrington, Time harmonic electromagnetic fields, New York: McGraw-Hill,
1941.

K. M. Chen, EE 836-electromagnetic waves 1: course notes, Michigan State
University, 1992.

G. Dural and D. L. Moffatt, "ISAR Imaging to Identify Basic Scattering
Mechanism," IEEE Trans. Antenna Propagat., vol. 42, No. I, pp.99-110, 1994

D. L. Mensa, High Resolution Radar Cross-Section Imaging. Norwood, MA:
Artech House, 1991.

S. Reichel, "Time domain of bistatic inverse scattering identity," Department of
Electrical Engineering, Michigan State University, 1995.

H. Li, N. H. Farhat, Y. Shen, and C. Werner, "Image Understanding and

174

Interpretation in Microwave Diversity Itnaging," IEEE Trans. Antennas Propagat.,
vol. 37, No. 8, pp.lO48-1057, August 1989.

175