i: L.
{5%

 

 

 

 

 

L)
«a. g

0.23.14

 

a: I 7. El:

542. sifiéxgmwflfigfi

[x

myo

LIBRARY
Michigan State
University

This is to certify that the
dissertation entitled

INVERSE MEDIUM SCATTERING FOR
ELECTROMAGNETIC WAVE PROPAGATION

presented by

Peijun Li

has been accepted towards fulfillment
of the requirements for the

PhD. degree in Applied Mathematics

 

 

 

Major Professor’? Signature

07/23 /0}‘—

Date

 

MSU is an Affirmative Action/Equal Opportunity Institution

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

03370360793

 

-AA

Y} I j

 

AUG 3 02005

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 c:/CI'RC/oateoue_.ind¢'p'.1s‘

 

 

INVERSE MEDIUM SCATTERING FOR
ELECTROMAGNETIC WAVE PROPAGATION

By

Peijun Li

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mathematics

2005

ABSTRACT

INVERSE MEDIUM SCATTERING FOR ELECTROMAGNETIC
WAVE PROPAGATION

By

Peij un Li

This thesis focuses on the study of continuation methods for solving inverse
medium scattering problems from electromagnetic wave propagation.

The first Considers a time—harmonic electromagnetic plane wave incident on a
medium enclosed by a bounded domain in R3. A continuation method for the inverse
medium scattering problem, which reconstructs the scatterer of an inhomogeneous
medium from boundary measurements of the scattered wave, is developed. The algo-
rithm requires multi-frequency scattering data. Using an initial guess from the Born
approximation, each update is obtained via recursive linearization on the wavenumber
k by solving one forward problem and one adjoint problem of Maxwell’s equations.

In part two, we consider the inverse medium scattering problem for Helmholtz’s
equation at fixed frequency. A new continuation method for the inverse medium
scattering is developed. The algorithm requires only singlefrequency scattering data.
Using an initial guess from the Born approximation, each update is obtained via
recursive linearization on the spatial frequency of a one-parameter family of plane

waves by solving one forward and one adjoint problem of the Helmholtz equation.

To my son, Ray.

iii

ACKNOWLEDGMENTS

With great pleasure, I would like to express my sincere gratitude to all those who
assisted and helped me while I was studying and doing research for this dissertation.

First and foremost I am deeply grateful to Prof. Gang Bao, my thesis advisor, for
his guidance and assistance, for his encouragement and support, and for his kindness
and hospitality during my graduate study at Michigan State University. I would like
to thank him for suggesting problems and giving directions which sheds light not only
on my current research but also on my future career.

I wish to express special thanks to Prof. Tien—Yien Li, who shows me some
mathematical insights and gives me his support without reservation. I would like
to thank Prof. Zhengfang Zhou for his inspiring discussions and suggestions. Also
I would like to thank Prof. Chichia Chiu, Prof. Guowei Wei, and Prof. Keith
Promislow for their time and concern.

I am forever indebted to my wife, Ying Zhang, and my family for their enduring

patience, constant support, and encouragement.

iv

TABLE OF CONTENTS

LIST OF FIGURES
LIST OF TABLES
Introduction

1 Inverse Medium Scattering Problems for Electromagnetic Waves
1. 1 Introduction ................................

1.2 Analysis of the Variational Problem for Forward Scattering ......
1.3 Inverse Medium Scattering . . .......................
1.3.1 Low-Frequency Modes of the Scatterer .............

1.3.2 Recursive Linearization ......................

1.4 Numerical Experiments ..........................

1.5 Concluding Remarks ...........................

2 Inverse Medium Scattering at Fixed Frequency

2. 1 Introduction ................................
2.2 Analysis of the Scattering Map .....................
2.3 Inverse Medium Scattering ........................

2.3.1 Born Approximation .......................

2.3.2 Recursive Linearization ......................
2.4 Numerical Experiments ..........................
2.5 Concluding Remarks ...........................

BIBLIOGRAPHY

vi

ix

A

14
15
18
25
35

42
42
48
57
58
60
66
87

88

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

LIST OF FIGURES

Parsity pattern of an FE matrix with 1820 unknowns: (a) original
ordering; (b) RCM ordering ........................ 27
Example 1: Integrals at different wavenumbers for the fixed incident
angle 6 = 7r/ 3 and ti = 7r/ 3. Solid curve: the exact integral value of
the left-hand side of (1.3.2), +: the computed integral value of the first
term of the right—hand side of (1.3.2), at: the computed integral value of
the second term of right-hand side of (1.3.2), o: the computed integral
value of the right-hand side of (1.3.2) ................... 29
Example 1: Integrals with different 6 for the fixed wave number k = 2.0
and (b = 7r / 3. Solid curve: the exact integral value of the left-hand side
of (1.3.2), o: the computed integral value of the right hand-side of (1.3.2). 30
Example 1: Integrals with different ¢ for the fixed wave number k = 2.0
and 9 = 7r / 3. Solid curve: the exact integral value of the left-hand side
of (1.3.2), o: the computed integral value of the right-hand side of (1.3.2). 31
Example 1: (a) the slice :1: = 0 of the true scatterer; (b) the slice :1: = 0
of the reconstruction ............................ 32
Example 1: (a) the slice 3] = 0 of the true scatterer; (b) the slice y = 0
of the reconstruction ............................ 33
Example 1: (a) the slice 2 = 0 of the true scatterer; (b) the slice 2 == 0
of the reconstruction ............................ 34
Example 2: Integrals at different wavenumbers for the fixed incident
angle 6 = 7r/ 3 and (25 = 7r / 3. Solid curve: the exact integral value of
the left-hand side of (1.3.2), +: the computed integral value of the first
term of the right-hand side of (1.3.2), *: the computed integral value of
the second term of right-hand side of (1.3.2), o: the computed integral
value of the right-hand side of (1.3.2) ................... 36
Example 2: Integrals with different 6 for the fixed wave number k = 3.0
and 45 = 7r / 3. Solid curve: the exact integral value of the left-hand side
of (1.3.2), o: the computed integral value of the right hand-side of (1.3.2). 37

vi

1.10

1.11

1.12

1.13

2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
2.11

2.12

2.13
2.14
2.15
2.16
2.17

Example 2: Integrals with different (76 for the fixed wave number k = 2.0
and 9 = 7r / 3. Solid curve: the exact integral value of the left-hand side
of (1.3.2), o: the computed integral value of the right-hand side of (1.3.2). 38

Example 2: (a) the slice :r = 0 of the true scatterer; (b) the slice :7: = 0

of the reconstruction ............................ 39
Example 2: (a) the slice y = -0.5 of the true scatterer; (b) the slice

y = —0.5 of the reconstruction. ..................... 40
Example 2: (a) the slice 2 = 0 of the true scatterer; (b) the slice 2 = 0

of the reconstruction ............................ 41
Evanescent plane wave at k0 = 4.0 with 77 = 4.7 ............. 44
Evanescent plane wave at k0 = 4.0 with 77 = 8.0 ............. 45
Example 1: true scatterer ql ........................ 68
Example 1: Born approximation. .................... 69
Example 1: reconstruction of ql at 77 = 10.2 ............... 70
Example 1: reconstruction of ql at 77 = 8.4 ................ 71
Example 1: reconstruction of q] at 77 = 6.6 ................ 72
Example 1: reconstruction of ql at 77 = 4.8 ................ 73
Example 1: reconstruction of ql at 77 = 3.0 ................ 74
Example 1: reconstruction of ql at 77 = 0.0 ................ 75
Example 1. Relative error of the reconstruction ql at different

wavenumber kg. 0: reconstruction at [to = 10.0; *: reconstruction
at k0 = 8.0; Cl: reconstruction at .170 = 6.0; +: reconstruction at
k0 = 4.0 ................................... 76

Example 1. Relative error of the reconstruction ql at different step size
677. o: reconstruction at (577 = 0.6; *: reconstruction at 677 = 1.2; D:
reconstruction at 677 = 2.0 ......................... 77
Example 2: true scatterer of q2. ..................... 79
Example 2: reconstruction of qg ...................... 80
Example 2: contour view of the true scatterer q2. ........... 81
Example 2: contour view of the reconstructed scatterer q2. ...... 81
Example 2. Evolution of slice for the reconstruction q2. Solid line:

true scatterer; Circle: reconstruction. Top row from left to right: re—
construction at 77 = 14.45; reconstruction at 77 = 13.60; reconstruc-
tion at 77 = 12.75; middle row from left to right: reconstruction at
77 = 10.20; reconstruction at 77 = 8.50; reconstruction at 77 = 6.80; bot-
tom row from left to right: reconstruction at 77 = 5.10; reconstruction
at 77 = 2.55; reconstruction at 77 = 0.0. ................. 82

vii

2.18 Example 3:
2.19 Example 3:
2.20 Example 3:
2.21 Example 3:

true scatterer of q3. ..................... 83

reconstructed scatterer of q3 ................. 84
image view of the true scatterer q3. ............ 85
image view of the reconstructed scatterer q3. ....... 86

viii

1.1

1.2
1.3

2.1

LIST OF TABLES

Recursive linearization reconstruction algorithm for inverse medium

scattering. ................................. 24
Relative error at different wavenumbers of Example 1 .......... 28
Relative error at different wavenumbers of Example 2 .......... 35
Recursive linearization reconstruction algorithm ............. 65

ix

Introduction

Consider an electromagnetic plane wave propagating in a homogeneous medium. In
the absence of any inhomogeneities, the wave will continue to propagate and nothing
of physical interest will happen. However, if there are inhomogeneities present, the
wave will be scattered and we can express the total field as the sum of original incident
wave and the scattered wave. The behavior of the scattered wave will depend on both
the incident wave and the nature of the inhomogeneities in the medium. The direct
problem, given this information, is to find the scattered wave. The inverse problem
takes this answer to the direct scattering problem as its starting point and ask, what
is the nature of the inhomogeneities that gave rise to such scattered field behavior?

This thesis focuses on the study of continuation methods for solving inverse scat-
tering problems. These inverse scattering problems arise naturally from diverse ap—
plications such as medical imaging, nondestructive testing, and geophysical explo-
ration [10]. Two major difficulties for solving these inverse problems by optimization
methods are the ill—posedness and the presence of many local minima. Based on
multi-experimental data with physical parameters, we have developed regularized
continuation methods to solve the inverse medium scattering for three-dimensional
time-harmonic Maxwell’s equations and the inverse medium scattering for Helmholtz
equation at fixed frequency.

In Chapter 1, we consider a time—harmonic electromagnetic plane wave incident

on a medium enclosed by a bounded domain in R3. Existence and uniqueness of the

variational problem for direct scattering are established. An energy estimate for the
scattered field with a uniform bound with respect to the wavenumber is obtained in
the case of low frequency on which the Born approximation is based.

The inverse medium scattering problem is to determine the scatterer from the
measurements of near field currents densities on the boundary, given the incident
field. Although this is a classical problem in inverse scattering theory, little is known
on reconstruction methods, especially in the three dimensional case, due to the non-
linearity, ill-posedness, and the large scale computation associated with the inverse
scattering problem. Our goal is to present a recursive linearization method that solves
the inverse medium scattering problem of Maxwell’s equations in three dimensions.
The algorithm requires multi-frequency scattering data, and the recursive lineariza-
tion is obtained by a continuation method on the wavenumber It. It first solves a linear
equation (Born approximation) at the lowest k, which may be done by using the Fast
Fourier Transform. Updates are subsequently obtained by using higher and higher
wavenumber k. Using the idea of Kaczmarz method, we use partial data to perform
the nonlinear Landweber iteration at each stage of the wavenumber It. For each iter-
ation, one forward and one adjoint state of the Maxwell’s equations are solved. This
may be implemented by using the symmetric second order edge (Nédélec) elements.

Chapter 2 considers a time-harmonic electromagnetic plane wave incident on a
medium enclosed by a bounded domain in R2. Existence and uniqueness of the
variational problem for the direct scattering are established. An energy estimate for
the scattered field is obtained on which the Born approximation is based. Fréchet
differentiability of the scattering map is examined.

The main purpose of this work is to present a single-frequency inversion method,
and to demonstrate the efficiency of a new continuation method for the inverse
medium scattering. The illuminating fields, including the high spatial frequency

evanescent plane waves, form a one-parameter family of plane waves. When the ob-

ject is probed with the high spatial frequency of the evanescent plane waves, only a
thin layer of the object is penetrated. Corresponding to this exponentially decaying
incident field, the scattered field which is measured on the boundary contains infor-
mation of the object in that thin layer. Such a measurement is entirely inadequate
to determine the whole object. However, the measurement may be used to obtain an
approximation. The evanescent plane waves with less spatial frequency are needed
to illuminate the object. While the probing energy penetrates a thicker layer of the
object, the relation between the measurement and the scatterer to be recovered in
the thicker layer becomes more nonlinear. These nonlinear equations can be consid-
ered as perturbations to the already solved equations at the previous thicker layers,
and therefore can be continually and recursively linearized with standard perturba-
tional techniques. Thus, the recursive linearization is a continuation method on the

transverse direction of the incident waves, which controls the depth of its penetration.

CHAPTER 1

Inverse Medium Scattering
Problems for Electromagnetic

Waves

1. 1 Introduction

Consider the systems of time harmonic Maxwell’s equations in three dimensions

v x Et = iwu*Ht, (1.1.1)

V x Ht = —iw5*Et, (1.1.2)

where Et and H t are the total electric field and magnetic field, respectively; w > 0 is
the frequency; and 5* and 71* are the electric permittivity and the magnetic perme-
ability, respectively. Denote by 60 > 0, 71.0 > 0 the permittivity and permeability of
the vacuum. The fields are further assumed to be nonmagnetic; i.e. 71* = 710. Rewrit-
ing 5* == 805, e = 1 + q(:r) is the relative permittivity, where q(a:) is the scatterer,
which is assumed to have a compact support, and 93(q(2:)) > —1.

Taking the curl of (1.1.1) and eliminating the magnetic field H t, we obtain the

uncoupled equation for the electric field Et:
V x (V x Et) — kQEEt = 0, (1.1.3)

where k = ca 50770 is called the wavenumber, satisfying 0 < kmin S k S kmax < 00.

The total electric field Et consists of the incident field Bi and the scattered field E:
Et = Ei + E.
Assume that the incident field is a plane wave of the normalized form [10]
Ei = ikfieikl‘ ‘ ’53, (1.1.4)

where 717 E 82 is the propagation direction and 75' E 82 is the polarization satisfying

75' - 7'7' 2 0. Evidently, such an incident wave satisfies the homogeneous equation
V x (V x E?) — k213i = 0. (1.1.5)
It follows from the equations (1.1.3) and (1.1.5) that the scattered field satisfies
v x (v x E) — WEE = k2q(z)Ei. (1.1.6)

In addition, the scattered field is required to satisfy the following Silver-Muller radi-

ation condition:

In 7“[VxEx-:—:—ikE]=0,

li
r -> 00
where 7‘ = |:i:|. In practice, it is convenient to reduce the problem to a bounded domain

by introducing an artificial surface. Let Q be the compact support of the scatterer

q(:r). Assume that R > 0 is a constant, such that the support of the scatterer, Q,

is included in the ball B = (a: 6 1R3 : |:1:| < R}. Let S be the sphere of the ball,
i.e. S = {as 6 R3 : [ml = R}. Denote V the outward unit normal to S. A suitable
boundary condition then has to be imposed on S. For simplicity, we employ the first

order absorbing boundary condition (impedance boundary condition) [22] as
ux(VxE)+ikux(1/><E)=0 onS. (1.1.7)

Given the incident field E', the forward problem is to determine the scattered field
E for the known scatterer q(:c), which is assumed further to be in L°°(B). Based on
the Helmholtz decomposition and a compact imbedding result, the forward problem
is shown to have a unique solution for all but possibly a discrete set of wavenumbers.
Furthermore, an energy estimate for the scattered field, with a uniform bound with
respect to the wavenumber, is given in the low frequency case. The estimate provides
a theoretical basis for our linearization algorithm. For numerical solution of the for-
ward scattering problem in an open domain, the reader is referred to [25, 26, 27, 32]
and references therein. The inverse medium scattering problem is to determine the
scatterer q(:r) from the measurements of near field current densities, the tangential
trace of the scattered field 1/ x E I 3, given the incident field. Although this is a clas-
sical problem in inverse scattering theory, little is known on reconstruction methods,
especially in the three dimensional case, due to the nonlinearity, ill-posedness, and
large scale computation associated with the inverse scattering problem. We refer the
reader to [1, 13, 19, 20, 34] for related results on the inverse medium problem. See
[10] for an account of recent progress on the general inverse scattering problem.

The goal of this work is to present a recursive linearization method that solves
the inverse medium scattering problem of Maxwell’s equations in three dimensions.
The reader is referred to [3, 8] for recursive linearization approaches for solving the

inverse medium scattering problems in two dimensions. Our algorithm requires multi—

frequency scattering data, and the recursive linearization is obtained by a continuation
method on the wavenumber. It first solves a linear equation (Born approximation)
at the lowest wavenumber, which may be done by using the fast Fourier transform
(FFT). Updates are subsequently obtained by using higher and higher wavenumbers.
Following the idea of the Kaczmarz method [29, 30, 13], we use partial data to per-
form the nonlinear Landweber iteration at each wavenumber. For each iteration,
one forward and one adjoint state of Maxwell’s equations are solved, which may be
implemented by using the symmetric second order edge (Nédélec) elements.

The plan of this paper is as follows. Analysis of the variational problem for
forward scattering is presented in section 1.2. Based on the Helmholtz decomposition,
a compact imbedding result, and the Lax—Milgram lemma, the well-posedness of the
forward scattering is proved. An important energy estimate is given. Section 1.3 is
devoted to the numerical study of inverse medium scattering. Using the initial guess
of the reconstruction derived from the Born approximation, a regularized iterative
linearization algorithm is proposed. Numerical examples are presented in section 1.4.

The paper is concluded with some remarks and future directions in section 1.5.

1.2 Analysis of the Variational Problem for For-
ward Scattering

In this section, the variational formulation for the forward scattering problem is dis-
cussed. The analysis provides a criterion for weak scattering, which plays an impor-
tant role in the inversion algorithm.

To state our boundary value problem, following [28], we first introduce the stan-

dard Sobolev spaces:

L?(S) = {u E (L2(S))3: l/ - u = 0 onS},
H7](B) = {u E H1(B): u = 0 onS},
H(curl, B) = {u 6 (L2(B))3 : V x u E (L2(B))3},

Himp(cur1, B) = {u E H(curl, B) : 1/ x u E L?(S)},

where H,mp(curl, B) is an appropriate subspace of H (curl, B) for solving problems
involving the impedance boundary condition. Correspondingly, these spaces are

equipped with the norms

H u ”143(5):“ u |](L2(S))3

+ H W H2

2 _ , 2 2 .
H “ ”H(c11r1,B)—-ll ” ”(L2(B))3 + “ V X “' ”(112(3))3’
, 2 _ , 2 2
H u llHimp(curl, 8)-” u llH(curl, B) + H V X 11 IILflS) '

For convenience, denote the (L2(B))3 and (L2(S))3 inner products by

(u, v) =/u-Ed:c and (11,21) 2/71-6ds,
B S

respectively, where the overline denotes the complex conjugate. Introduce the bilinear

form a : Himp(curl, B) x Himp(curl, B) ——> (C;

a(E,q5) = (v x E, v x a) — k2(sE,d>) + my x E,u x a),

and the linear functional on Himp(curl, B);
has) = can} <1).

Then we have the weak form of the boundary value problem (1.1.6) and (1.1.7): find
E E Himp(cur1, B) such that

a<E, a = bus) w e H...p(cur1. B). (12.1)

Throughout the paper, C stands for a positive generic constant whose value may
change step by step but should always be clear from the context.
Before presenting the main result for the variational problem, we state several

useful lemmas. The reader is referred to [28] for detailed discussions and proofs.
Lemma 1.2.1 (Helmholtz decomposition). The spaces X and Y are closed sub-
spaces of Himp(curl, B). which is the direct sum of the spaces X and Y, i.e.,
Himp(curl, B) = X EB Y.
Here
X = {u E Himp(curl, B) : div(eu) = 0, in B}

and

Y = {vg: g e H3(B)}.

Lemma 1.2.2 (compact imbedding). The space X is compactly imbedded into the
space (L2(B))3.

Lemma 1.2.3 (Friedrichs inequality). There exists a positive constant C, inde-

pendent of the wavenumber, such that for all u E X

“ “ "<L2<B>)35 C (“ V x ”(122(8))“ + “ ” X ’“ “(Inert)“

Next we prove the well-posedness of the variational problem (1.2.1) and obtain
an energy estimate for the scattered field with a uniform bound with respect to the

wavenumber in the case of low frequency.

Theorem 1.2.1. If the wavenumber is sufiiciently small, the variational problem
(1.2.1) admits a unique weak solution in Himp(curl,B) given by E = u + Vp, while

u E X, p E H6(B). Furthermore, we have the estimate

H E “Himmcuflflfi Cthll/Q H q ”30(3), (1202)

where the constant C is independent of k and Q is the compact support of the scatterer.

PROOF. Using the Helmholtz decomposition, we take E = u + Vp and ct = u + V5,
for any u E X,f E H3(B). Observe that a(u,V§) = 0, for any 6 E H6(B), by the
definition of X. Therefore, we decompose the variational equation (1.2.1) into the

form
a(u,v) + a.(Vp,v) + a(Vp, V5) 2 b('v) + b(V{) Vt) E X, 6 E H3(B). (1.2.3)
First, we determine p E H6(B) by the solution of
am), Vs) = blvé) V6 6 HEB),
which gives explicitly

—<sz.vo = (qEévo V6 6 Haw).

10

The existence and uniqueness of the solution p in H6(B) may be proved by a direct

application of the Lax—Milgram lemma with the estimate
H W ”mamas Cum“? ll q “Loom, (1.2.4)
Rewrite (1.2.3) as
a(u, u) = b(u) — a(Vp, u) Vu E X (1.2.5)
and decompose the bilinear form a into a = a1 + k2a2, where

a1(u,u) = (V X u,V x v) + ik(1/ x u,V x u),

Using the inequality of arithmetic and geometric means, we conclude from Lemma

1.2.3 that al is coercive

Ia1(u,u)] Z Ck (H V X u ”(L2(B))3 + H V X 71. IIEL2(S))3) 2 Ck H u ”2Himp(curla B)

VuEX.

The continuity of the bilinear form a1 follows from the Cauchy—Schwarz inequality.

Next we prove the compactness of a2. Define an operator A : (L2(B))3 -—> X by
a1(.Au,v) = a2(u,u) Vu E X,
which gives
(V x Au,V x v) +ik(7/ x Au,1/ x u) = —(eu,u) Vu E X.

11

Using the Lax—Milgram lemma again, it follows that

H A“ ”Himp(cur1,B)S % II U ”(L2(B))3’ (1.2.6)

where the constant C is independence of 19. Thus A is bounded from (L2(B))3 to X,
and X is compactly imbedded into (L2(B))3. Hence A : (L2(B))3 —> (L2(B))3 is a
compact operator.

Define a function w E (L2(B))3 by requiring w E X and satisfying
a1(w, u) = (2(1)) — a(Vp, v) V?) E X.
More specifically, we have by using the Stokes formula that
a1(w, u) = k2(qu,u) +lk2(eVp,u) V7) 6 X.
It follows from the Lax-Milgram Lemma that
H w 11H,,,p(.,,.1, B): C (72191”? H «1 “130(3) +k H W ||( 12(3))3) .
An application of (1.2.4) yields
H w 11H,,,p(c,,,1, 3,: other” H q ”1400(3), (127)

Using the operator A, we can see that the problem ( 1.2.5) is equivalent to finding
u E (L2(B))3 such that
(I + k2A)u = to. (1.2.8)

When the wavenumber k: is small enough, the operator T + k2A has a uniformly

12

bounded inverse. We then have the estimate

l] u ”(L2(B))3S C I] w ”(L2(B))3’ (129)

where the constant C is independent of k. However, rearranging (1.2.8), we have

u = w — k2Au, so u E X and, by the estimate (1.2.6) for the operator A, we have

N U ||H,,,,p(cur1, B)S|| w “Himp(curl, B) +Ck “ “ ”(L2(B))3.

Combining the estimates (1.2.9) and (1.2.7) leads to

H U ”Himp(curl’ B)S CkQIQII/2 H q ”LOO(B) (1'2'10)

Finally, it. follows from the definition of the norm in H,,,,,,(curl, B) that

I] E ”Himp(curL B)SH U ”Himp(curl, B) + ]] VP ]|(L2(B))3

The proof is complete by noting the estimates (1.2.10) and (1.2.4) for sufficiently

small wavenumbers. Cl

Remark 1.2.1. The energy estimate of the scattered field (1.2.2) provides a criterion
for weak scattering. From this estimate, it is easily seen that firing any two of the
three quantities, i.e. the wavenumber, the compact support of the scatterer Q, and
the L°°(B) norm of the scatterer, the scattering is weak when the third one is small.
Especially for the given scatterer q(;r), i.e. the norm and the compact support are

fired, the scattering is weak when the wavenumber is small.

Remark 1.2.2. For a general wavenumber, from (1.2.8), the uniqueness and exis-
tence follow from the Fredholm alternative. If the scatterer q(;I:) is more regular, say

of 03 (B) [15], unique continuation may be used to prove the uniqueness and thus the

13

existence of the forward scattering problem (1.1.6), (1.1.7) for all k > 0. Otherwise,
if k2 is not the eigenvalue for Maxwell’s equations in the domain B, then the operator
I + k2A has a bounded inverse. However, the bound depends on the wavenumber.

Therefore, the constant C in the estimate (1.2.2) depends on the wavenumber.

From the above discussion, we have the following theorem on the well-posedness

of the variational problem (1.2.1).

Theorem 1.2.2. Given the scatterer q E L°°(B), for all but possibly a discrete set
of wavenumbers, the variational problem (1.2.1) admits a unique weak solution in

Himp(curl, B), given by E = u + Vp, while u E X, p E H6(B).

1.3 Inverse Medium Scattering

In this section, a regularized recursive linearization method for solving the inverse
medium scattering problem of Maxwell’s equations in three dimensions is proposed.
The algorithm, obtained by a continuation method on the wavenumber, requires mul-
tifrequency scattering data. At each wavenumber, the algorithm determines a forward
model which produces the prescribed scattering data. At a low wavenumber, the scat-
tered field is weak. Consequently, the nonlinear equation becomes essentially linear,
known as the Born approximation. The algorithm first solves this nearly linear equa-
tion at the lowest wavenumber to obtain low-frequency modes of the true scatterer.
The approximation is then used to linearize the nonlinear equation at the next higher
wavenumber to produce a better approximation which contains more modes of the
true scatterer. This process is continued until a sufficiently high wavenumber, where

the dominant modes of the scatterer are essentially recovered.

14

1.3.1 Low-Frequency Modes of the Scatterer

Rewrite (1.1.6) as
v x (v x E) — k2E = k2q(x)(Ei + E), (1.3.1)

where the incident wave is taken as Ei = ikp‘leikx ' 77.]. Consider a test function
F = ikp'geikx ' n2, where 732, 772 E S2 satisfy [7'2 - 7'52 = 0. Hence F satisfies (1.1.5).

Multiplying (1.3.1) by F and integrating over B on both sides, we have

/F-[V x (v x E)]dr—k2/ F-Edrz k2] q(;c)F.Eidx+k2/ q(:1:)F-Ed:z:.
B B B 8

Integration by parts yields

fE-[Vx(VXF)]dr+/[EX(VXF)-Fx(VXE)]-uds-—k2/F~Ed.r
13 s B

= 132/ q(;r)F-Eid:c + k2/ q(:1:)F-Ed2‘.
B B

W'e have, by noting (1.1.5),

q(.r)F - E‘ida: + k2] q(:1:)F - Edzr.

/[Ex(V><F)—FX(V><E)]~Vds=k2f
S B

8

Using the boundary condition (1.1.7) of the scattered field and the special form of

the incident wave E’ and F, we get

— /(V x E) - (fig x 772)eik$ ' 77st + / [1/ x (1/ x E)] ~figeik$ ' 732ds
3 S

=/Bq($)F.Eidx+/Bq(x)F-Edr.

15

A simple calculation yields

/ «new ' (a + Wake = —..——1:——, [(12 x E) - (it x a + u x 5219”“ ' fi2ds
B (P1'P2)k s
i _. 'ka: - it
+—_,——_,——— q 117 p ' Eel 2611C.
(P1'P2)k/E ( ) 2

(1.3.2)

From Theorem 1.2.1 and Remark 1.2.1, for a small wavenumber, the scattered
field is weak and the inverse scattering problem becomes essentially linear. Dropping

the nonlinear (second) term of (1.3.2), we obtain the linearized integral equation

- ., .-* _. 1 . . _.
/ (10(1))81kI - (711 + ”flit : —.,—_:——§ /(VXE)-(fi2X]Tg+VXfi2)elkx ”2618, (1.3.3)
B (P1'I92)k s

which is the Born approximation. The function (10(2) will be used as the starting
point for our recursive linearization algorithm.

Since the scatterer q0(:r) has a compact support, we use the notation
40(5) ___/ q0(17)eil"x'(n1+ n2)d:1:,
B
where (jobs) is the Fourier transform of q0(:1:) with E = k(fi1 + rig). Choose
7i,- 2 (sindj cos (75,-,sin 6,- sinqfij,cosdj), j = 1,2,

where 6], gb,- are the latitudinal and longitudinal angles, respectively. It is obvious that
the domain [0, 7r] x [0, 2n] of (6,, dj),j = 1, 2, corresponds to the ball {5 E 1R3 : If] 3
2k}. Thus, the Fourier modes of (70(5) in the ball {5 : [6| S 2k} can be determined.
The scattering data with the higher wavenumber must be used in order to recover

more modes of the true scatterer.

16

Define the data

—:——1:.——, /(u x E) - (112 x 52 + u x 172).?” ' 772113 for ICI 3 2k,
a(g) = (171192)]? 5

0, otherwise,

where C = C (k, 61, (1)1, 62, d2) E R3. The linear integral equation (1.3.3) can then be

formally reformulated as

reim'cat: . ..
([2,,qu d 0(4) (134)

Taking the inverse Fourier transform of (1.3.4) leads to

1 —i$'C - iy-C, : 1 _ig;.C
(27Tl3 /R3 6 [A23 q0(y)e dy]dC (2703 /Rs e C(CldC-

By the Fubini theorem, we have

 

 

I , ei(y — (p) -C = I e_i$ . C

Using the inverse Fourier transform of the Dirac delta function

 

 

1
(2703

 

go—x)< = _$
[1,, dc 6o ).

we deduce
1

—:r = e—ias-§
[Rheum )dy (2).)?» [11,, made

 

which gives

1 .
1' :: e—IIII'C . . .
(10(l (2703/1113 G<<>d< (135)

 

In practice, the integral equation (1.3.5) is implemented by using the FFT.

17

1.3.2 Recursive Linearization

As discussed in the previous section, when the wavenumber is small, the Born ap-
proximation allows a reconstruction of those Fourier modes less than or equal to 2k
for the function q(a:). We now describe a procedure that recursively determines qk
at k = k,- for j = 1,2, with the increasing wavenumbers. Suppose now that the
scatterer q]; has been recovered at some wavenumber k, and that the wavenumber k
is slightly larger than k. We wish to determine qk, or equivalently, to determine the

perturbation

5Q = (11; - (1,1.-

For the reconstructed scatterer q,;, we solve at the wavenumber k the forward

scattering problem

v x (v x E) — 1:2(1+ q&)E = quZEl, :1: e B, (1.3.6)

I/ x (V X E) +ik1/ x (1/ x E) = 0 011 S. (1.3.7)
For the scatterer qk, we have

v x (v x E) — 1:2(1+ qk)E =1~.2qu"', :1: e B, (1.3.8)

ux(VxE)+ik1/x(uxE)=0 onS. (1.3.9)

Subtracting (1.3.6), (1.3.7) from (1.3.8), (1.3.9), and omitting the second order small-

ness in 6q and in 6E = E — E, we obtain

v x (v x 6E) — 18(1 + q§)6E = k25q(Ei + E), x e B, (1.3.10)

ux(Vx6E)+ika(z/X6E)=0 onS. (1.3.11)

18

For the scatterer qk and the incident wave El, we define the map S (qk, E2) by
S (qk, E2) = E,

where E is the scattered field at the wavenumber k. Let '7 be the trace operator to

the boundary S of the ball B. Define the scattering map
M(qk, El) = 75(qk, El)-

It is easily seen that the scattering map 1W(qk, Ei) is linear with respect to El but
is nonlinear with respect to qk. For simplicity, denote M (qk, Ei) by M (qk). By the

definition of the trace operator, we have
ll’l(qk) = I/ X E's.

We refer to [1] for the Fréchet differentiability of the scattering map. Let DM(ql} ) be

the Fréchet derivative of hf(qk), and denote the residual operator
R(qk) = l/ X 6Els.

It follows from [1] that

DA/[(qk)6q = R(q];) (1.3.12)
The regularized least-squares solution of (1.3.12) is

52 = [01 + DM*(Q,;)DM(q,;.)]"DM*((I,~€)R(qk),

where DM*(q]~€) is the adjoint operator of DhI(ql~c), I is the identity operator, and

a is some suitable positive number. In practice, the main difficulty is the enormous

19

computational cost of solving linear systems with huge full matrix. Here, we consider
an alternative way of solving (1.3.12) which is much less computationally demanding.

To state the approach, we first examine the boundary data 1/ x E(:1:; 0, d); k). Here,
the variable a: is the observation point, which has two degrees of freedom since it is
on the sphere S. The terms 0, a are latitudinal and longitudinal angles respectively
of the incident wave Ei. At each frequency, we have four degrees of freedom, and
thus data redundance, which may be addressed by fixing one of the incident angles,
say 6. Define d),- = (j — 1) >1: 27r/m,j = 1, ..., m, and the residual operator

31(9):.) = V X E($;9.¢j;k)|s — V X E(I§91¢ji fills)

where m is the total number of the incident waves or sweeps, and E (.r; 9, (p); k.) is the
solution of (1.3.6), (1.3.7) with the incident wave of longitudinal angel ()6, and the
scatterer ql} . Instead of solving (1.3.12) for all incident waves simultaneously, we may
solve it for one incident wave at a time while updating the residual operator after

each determination of the incremental correction (Sq. Thus, for each incident wave

with incident angle c5}, we consider the equation
llIJ-(qk) = 1/ x E(;I:;6,c'oj; k)]5, (1.3.13)

where A1j(qk) is the scattering map corresponding to the incident wave with longitu-

dinal angle (75,-. It follows from [1] that
DMACIIQCSQJ = Rflqé)» (1-3-14)

where Dll/Ij( ) is the Fréchet derivative of the scattering map ll! 1((1kl' The nonlinear

‘11:.

20

[Landweber iteration for (1.3.13) yields
5% = filelIflql‘clRAQg), (1.3.15)

where DM;(ql~€) is the adjoint operator of DAIJ-(qg), and 5;, is some relaxation pa-

rameter [14].

Remark 1.3.1. For a fixed wavenumber, the stopping index of nonlinear Landweber
iteration (1.3.15) could be determined from the discrepancy principle. However, in
practice, it is not necessary to do many iterations. Our numerical results indicate
that the iterative process for different incident angles oj,j = 1, ..., m, is sufficient to

obtain reasonable accuracy.

Next, we discuss the role of the relaxation parameter [3,. in the iteration (1.3.15),
which may be understood more clearly by considering the iteration from a different
point of view.

Consider the optimization problem of (1.3.13),

nan I] ll'IJ-(qk) —— 1/ x E(x;6. at]; k) “$11209”? (1.3.16)

The first order optimality condition for the problem (1.3.16) is given by
17113017,) (MAGIC) — V X ECU; 9. 6.51; k)) IS = 0- 03-17)

To solve the optimality equation (1.3.17), the time marching scheme proposed in [33]

consists of finding the steady state of the following parabolic equation:

dq. *
7:- = DM, ((1);) (u x Ema, a; k) — M.<q,,)) Is.

21

The numerical solution could be computed from the explicit method
502‘ = TDMflqklRAqfc),

where r is the discretized time step. Thus, the relaxation parameter Bk is essentially
the step size of time marching, whose length is restricted by the stability of the explicit
method.

In order to compute the correction 6q,, we need some efficient way to compute

Dh/[;(q,~€)RJ(q,}), which is given by the following theorem.
Theorem 1.3.1. Given the residual R,(q,~c), there exits a function F,- satisfying the

adjoint equations

v x (v x 1“,) —1:2(1+q)F, = 0, x e B, (1.3.18)

V x F, — iku x F,- = R,(q§) onS, (1.3.19)

such that the adjoint Fréchet derivative DilI,‘.'(q,~i‘) satisfies

[DM;-'<q,~,,)R.-<q,~,)l (a) = Mia) + in» - 8(1), (1.3.20)

where E;- is the incident wave with the longitudinal angle a),- and E,- is the solution of
(1.3.6), (1.3.7) with the incident wave E}.
PROOF. Let. E,- be the solution of (1.3.6), (1.3.7) with the incident wave E}. Consider

the equations as follows:

v x (v x 6E) — k2(1+ qk)6E =1:25q(E§?+ E,), x e B, (1.3.21)

1/ x (V x (5E) + iku x (1/ x (SE) = 0 on S, (1.3.22)

22

and the adjoint equations (1.3.18) and (1.3.19), which take the variational form:

(v x F,,V x (p) — k2((1+q;)F,-,¢) —ik(u x F,,1/ x a)

= (Rj((1j€)aV X 49), V45 6 Himp(curl, B).

The existence and uniqueness of the weak solution for the adjoint equations may be
proved in the same way as for the scattered field. The proof is omitted.
Multiplying equation (1.3.21) with the complex conjugate of F,- integrating over

B on both sides, we obtain

/F,-- [Vx (VX6E)]dx—k2/
B

8(1-1- q,;)F, - dde = k72/136q(E,i+ E,) -F,-dx.

Integration by parts yields
[[6E x (V X F,) —F,- x (V X 6E)] ~1/ds = 1172/ 6q(E,2-’ + E,) - F,dx.
s B

Using the boundary condition (1.3.22), we deduce

 

/(1/ x 6E) - (V x F,- +ik1/ x F,)ds = 11:2/6q(E,f+ E) 'Fjdx.
s B
It follows from (1.3.14) and the boundary condition (1.3.19) that

/S[Dlll,- (q,~€)6q]-R —_k_)ds= k2L( 6q(E F,- dx.

We know from the adjoint operator Dll/I;(q,;,) that

 

f3 6qDM"(q ,C)R,(q,;)dz=k2 f3 6q(Ej+E,).E,-dx.

23

 

Initialization:
k = kmin smallest kmin
qo Born approximation
Reconstruction loop:
FOR k = kmin : kmax march along wavenumbers
FOR j = 1 : m perform m sweeps over incident angles
solve (1.3.6)-(1.3.7) for E,- one forward problem
solve (1.3.18)-(1.3.19) for F,- one adjoint problem

64;, = new; +5313
qi := q]: + (ngc
END
qr. == 2,7,”
END

q :2 qkmax final reconstruction

 

Table 1.1. Recursive linearization reconstruction algorithm for inverse medium scat-
tering.

Since this holds for any dq, we have

 

DM;(CI,})RJ'((17€)= k2(E; + E1) ' F,»

Taking the complex conjugate of the above equation yields the result. El

Using this theorem, we can rewrite (1.3.15) as
6q, = uremia) + E_,-(x)) - 13(3). (1.3.23)

Thus, for each incident wave with a longitudinal angle 45,-, we solve one forward
problem (1.3.6), (1.3.7) and one adjoint problem (1.3.18), (1.3.19). Since the adjoint
problem has a variational form similar as the forward problem, we need to compute
essentially two forward problems at each sweep. Once (Sq,- is determined, q,~c is updated

by qk + 6q,. After completing the mth sweep, we get the reconstructed scatterer Q},

24

at the wavenumber k.
The recursive linearization for inverse medium scattering of Maxwell’s equations

can be summarized in Table 1.1.

1.4 Numerical Experiments

In this section, we discuss the numerical solution of the forward scattering problem
and the computational issues of the recursive linearization algorithm.

As for the forward solver, we adopt the edge elements which were developed orig-
inally for the finite element solution of Maxwell’s equations [31, 22] in the beginning
of the 19808. From the mathematical point of view, these are natural approxima-
tion spaces for the Hilbert space H (curl, B), which is the adequate functional space
for the variational formulation of Maxwell’s equations. Vector fields in such finite
element (FE) spaces have continuous tangential traces, which is consistent with the
physics. Therefore, the natural degrees of freedom for these elements are related to
tangential traces along edges or faces. Here, we take the symmetric second order
tetrahedral edge elements [23]. When the unknowns are ordered according to the re-
verse Cuthill-McKee (RCM) ordering [16], the profile of FE matrix is highly banded
which improves the condition number of the FE coefficient matrix. The sparse large-
scale linear system can be most efficiently solved if the zero elements of the coefficient
matrix are not stored. We use the commonly used compressed row storage (CRS)
format, which makes no assumptions about the sparsity structure of the matrix and
does not store any unnecessary elements. In fact, from the variational formula of our
direct problem (1.2.1), the coefficient matrix is complex symmetric. Hence, only the
lower triangular portion of the matrix needs to be stored. Figure 1.1 shows a typical
sparsity pattern of an FE matrix with 1820 unknowns from the symmetric second

order edge element. Regarding the linear solver, either biconjugate gradient (BiCG)

25

or quasi-minimal residual (QMR) algorithms with diagonal preconditioning may be
employed to solve the sparse, symmetric, and complex system of the equations. It
appears for our experiments that the QMR is more efficient.

In the following, we present two numerical examples where the number of the inci-
dent wave 777. = 20, the incident latitudinal angle 6 = 0, and the incident longitudinal
angle (7’), = (j — 1) =1: 27r/m, j = 1, ...,m. The relaxation parameter 6;, is taken to be
0.1 / k for the tested examples. For stability analysis, some relative random noise is

added to the data, i.e. the tangential trace of the electric field takes the form
1/ x E|5 := (1+ orand) - (1/ x EIS).

Here, rand gives uniformly distributed random numbers in [—1,1], and a is a noise
level parameter taken to be 0.02 in our numerical experiments. Define the relative

error by

_ 1/2

_ (214,1: lqijk - qijkIQ)

e? _ 2 1/2
(22',ch Iqijkl )

 

where (j is the reconstructed scatter and q is the true scatterer.

Example 1. Reconstruct a scatterer defined by

 

$2 yz Z2 x2 yz 22
1" —+—+— f0r-—+——.—+—<1,
q(zv. y, 2) = 12 0.82 0.52 12 0,82 .52 —
0: otherwise.

The compact support of this scatterer is an ellipsoid contained in the unit ball. For
simplicity, we take r'il 2 fig and p} = 772 to test the forward solver. The numerical
results are shown in Figure 1.2-1.4. In Figure 1.2, for the fixed incident latitudinal
angle 6 = 7r / 3 and the longitudinal angle (75 = 7r / 3, the forward problem is solved at
different wavenumbers. In Figure 1.3 and 1.4, for the fixed wavenumber k = 2, the

numerical results are shown with different latitudinal angles 6 E [0, 7r] (fix (75 2 rr / 3)

26

 

 

 

ii-
200.. .6 J
”'Iuj.
it"‘} “h
400i:.,a '."'J;I ]
. I'
1 ., -.I
iif-I .' 5: are.
600‘”. I '-. ‘ ‘0' g .]
f" 4 ‘ 3'
hr. vulrfi.s:tur, "'Sw'r.“
800 "o r“ O... ." ... . .f‘n {C q
5’... .. I u b
".' ' ’2. " T.)"‘.".’
B
1 - ‘ , -
°°° ., . 1:
".1";
IF- .3 t: 'i i.e.li:
1200» . ' I p .5”; i. .
3 S ' lo -
. ' 5 SJ- "Q,
' - u . " ' . ' .' 3'2“. 0‘s
1400» - I: . f; ,fjtih .
t ’ .‘ri' I _'.'}o.’
"' ‘ . l" .. . in." o ..
.- - .- .. .. ‘rarfie- ..-.-..-r..-;-
1600’ ' n 1' I' f...- “b . '?I"r -
J I ' .‘o .5 ' ' ‘ .* .. .. ‘0 I '0‘“, i-
. ' - . . . ‘3'“00'0' b ." -.o.r... ,.4 U
u - F ‘ ' .3 ..I . U; o.:.. I. ..r
1800» _ , .. l - . I... 1. .b- ,.

0 200 400 600 800 ~ 1000 1200 1400 1600 1800

(a)

 

400

600

800

1200

1400

1600

 

 

 

1 800 P I l L L l I ’
0 200 400 600 800 1000 1200 1400 1600 1800

(b)

Figure 1.1. Parsity pattern of an FE matrix with 1820 unknowns: (a) original order-
ing; (b) RCM ordering.
27

 

k 1 2 3 4 5 6
62 0.5494 0.4876 0.3197 0.1856 0.1534 0.0895

 

 

Table 1.2. Relative error at different wavenumbers of Example 1.

and with different longitudinal angles qb E [0, 2n] (fix 6 = 7r/3), respectively. It is
easily seen from Figure 1.2 that the first term of the right-hand side of the integral
equation (1.3.2) is dominant compared with the second (nonlinear) term when the
wavenumber is small, which validates the Born approximation. Figure 1.5-1.7 (a)
show the slices of the true scatterer and Figure 1.5-1.7 (b) give the reconstruction
at the wavenumber k = 6. The relative errors are shown in Table 1.2 at different
wavenumbers.

Example 2. Reconstruct a scatterer defined by

511157;) - Sin(($2 + (y + 05)? + 2%) for 2:2 + (v + 0.5)2 + z"’ s 0.42,
4 .
q(x, y, 2) = sin(2—75r) — sin((x2 + (y —— 0.5)2 + zz)7r) for x2 + (y - 0.5)2 + Z2 _<_ 0.42,
0, otherwise.

The compact support of this scatterer is two isolated balls with the same radius of
0.4 and the centers at (0, -—0.5,0) and (0,0.5,0). For simplicity, we take ii] = 772
and 771 = {2‘2 in the test of the forward solver. The numerical results are given in
Figure 1.8-2.0. In Figure 1.8, for the fixed incident latitudinal angle 6 2 rr / 3 and the
longitudinal angle gt 2 7r / 3, the forward problem is solved at different wavenumbers.
In Figure 1.9 and 2.0, for the fixed wavenumber k = 3, the numerical results are shown
with different latitudinal angles 6 E [0, 7r] (fix (b = 7r / 3) and with different longitudinal
angles (,0 E [0, 27r] (fix 6 = 7r / 3), respectively. It is easily seen from Figure 1.8 that the
first term of the right-hand side of the integral equation (1.3.2) is dominant compared

with the second (nonlinear) term when the wavenumber k is small, which once again

28

 

0.45 1 r I l r

0.4

0.35

0.3

0.25

0.2

0.15

0.1

0.05

 

 

 

—

 

_005 l l l l
0

wave number k

Figure 1.2. Example 1: Integrals at different wavenumbers for the fixed incident
angle 6 = 7r/ 3 and (b = rr/ 3. Solid curve: the exact integral value of the left-hand
side of (1.3.2), +: the computed integral value of the first term of the right-hand side
of (1.3.2), *2 the computed integral value of the second term of right-hand side of
(1.3.2), o: the computed integral value of the right-hand side of (1.3.2).

29

 

0.32 r r 1 r r r

0.3 - o g 1
0.28 - d
0.26 - 0 E
0.24 - -
0.22 - g . -

0.2 ~— 2

 

 

1 l I

1.5V 2 2.5 3 3.5
incident latitudinal angle (radian)

 

0.18 1
0 0.5

_Lr—

Figure 1.3. Example 1: Integrals with different 6 for the fixed wave number k = 2.0
and (75 = 7r / 3. Solid curve: the exact integral value of the left-hand side of (1.3.2), o:
the computed integral value of the right hand-side of (1.3.2).

30

 

0.23 1 77 l I 1 I

l
l

0.22

0.21 - 1

0.2 r _

0.19

I
l

0.18 ’ s ' ' _

 

 

0.166 L 1 10 m 1 1 O

1 2 3 4 5 6 7
incident longitudinal angle (radian)

 

Figure 1.4. Example 1: Integrals with different a) for the fixed wave number k = 2.0
and 6 = 7r / 3. Solid curve: the exact integral value of the left-hand side of (1.3.2), o:
the computed integral value of the right-hand side of (1.3.2).

31

 

Figure 1.5. Example 1: (a) the slice x = 0 of the true scatterer; (b) the slice x = 0 of
the reconstruction.
32

 

Figure 1.6. Example 1: (a) the slice y = 0 of the true scatterer; (b) the slice y = 0 of
the reconstruction.

33

 

Figure 1.7. Example 1: (a) the slice 2 = 0 of the true scatterer; (b) the slice 2 = 0 of
the reconstruction.

34

 

k 1 2 3 4 5 6 7
e2 0.6963 0.6479 0.5891 0.4951 0.3376 0.2568 0.2221

 

 

Table 1.3. Relative error at different wavenumbers of Example 2.

validates the Born approximation. Figure 1.11-1.13 (a) show the slices of the true
scatterer, and Figure 1.11-1.13 (b) give the reconstruction at the wavenumber k = 7.

The relative errors are shown in Table 1.3 at different wavenumbers.

1.5 Concluding Remarks

The proposed recursive linearization algorithm is stable and efficient for solving the
inverse medium scattering problem with multiple frequency scattering data in three
dimensions. Theoretically, scattering data with even higher—wavenumbers could be
used to recover more complicated scatterers which contain higher frequency features,
i.e. more Fourier modes. However, the difficulty lies in the fact that the forward
model becomes difficult to solve due to the highly oscillatory nature of the solution.
For a larger k, the mesh size has to be smaller, which makes numerical solution more
expensive. Finally, we point out two important future directions of this research. The
first concerns with the convergence analysis of the recursive linearization algorithm,
which is currently in progress and will be reported elsewhere. Another challenging
project is to develop an efficient algorithm for the inverse medium scattering with

fixed frequency scattering data.

 

0.1

0.08

0.06

0.04

0.02

—0.02

-0.04

 

 

.—

J l

1 2 3 4 5 6
wave number k

_

 

 

-0.06 l
0

Figure 1.8. Example 2: Integrals at different wavenumbers for the fixed incident
angle 6 2 rr/ 3 and d) = 7r/ 3. Solid curve: the exact integral value of the left—hand
side of (1.3.2), +: the computed integral value of the first term of the right-hand side
of (1.3.2), *: the computed integral value of the second term of right-hand side of
(1.3.2), o: the computed integral value of the right-hand side of (1.3.2).

36

 

 

 

 

0.08 1 1 1 r I I
c C
0.06 " _,
c "
0.04 *- "
0.02 r -.
o “ 0
-0.02 _ .1
e ’l
—0.04 ~ 7
s 4
o
-0.06 1 I 1 I 1 1
0 0.5 1 1.5 2 2.5 3 3.5

incident latitudinal angle (radian)

Figure 1.9. Example 2: Integrals with different 6 for the fixed wave number k = 3.0
and (,1) = 7r / 3. Solid curve: the exact integral value of the left—hand side of (1.3.2), o:
the computed integral value of the right hand—side of (1.3.2).

37

 

0.08 ] fir I I 17 fl

0.06 ~ _

0.04 - j

I
l

0.02

0
It '
0 - a

I
l

-0.02

I
i

-0.04

O a S ’l

 

 

 

_0.% l 1 l l l
0 1 2 3 4 5 6 7

incident longitudinal angle (radian)

Figure 1.10. Example 2: Integrals with different 4) for the fixed wave number k = 2.0
and 6 = 7r / 3. Solid curve: the exact integral value of the left-hand side of (1.3.2), o:
the computed integral value of the right-hand side of (1.3.2).

38

 

Figure 1.1]. Example 2: (a) the slice x = 0 of the true scatterer; (b) the slice x = 0
of the reconstruction.

39

 

Figure 1.12. Example 2: (a) the slice y = —0.5 of the true scatterer; (b) the slice
y = —0.5 of the reconstruction.

40

 

Figure 1.13. Example 2: (a) the slice 2 = 0 of the true scatterer; (b) the slice 2 = 0
of the reconstruction.

41

CHAPTER 2

Inverse Medium Scattering at

Fixed Frequency

2. 1 Introduction

Consider the Helmholtz equation in two dimensions
A¢+k§(1+q(x))¢= 0, (2.1.1)

where (75 is the total field; k0 is the wavenumber, and q(x) > —1, which has a compact
support and a lower bound, is the scatterer.

Assume that the scatterer lies in the upper half plane R1 = {(x1,x2) E R2 : x2 >
0}. Denote the wave vector k = (r7, k(r7)), where 77 is the transverse part of the wave

vector and

k8 — 772 for kg 2 |77|,

k(n) =
i1/772 — 1:3 for kg < [77].

The number [77] is known as the spatial frequency.

42

The scatterer is illuminated by a one—parameter family of plane waves
.%=em*fi min

which gives explicitly

' / .2 _ 2
81(77371'1' k0 7] $2) for ’00 Z [Tl],

¢0($la$2) = . F7
emf1 _ 77 _ k0“ for k0 < [77].

The modes for which |r7| S k0 correspond to propagating plane waves while the
modes with In] > k0 correspond to evanescent plane waves. Therefore, the illuminat-
ing field could consist of high spatial frequency evanescent plane waves. They may be
generated at the interface of two media by total internal reflection [12, 18], which has
been in practical use for decades and primarily been used in near-field optics [6, 7]. A
recent review 011 the near-field microscopy and near-field optics may be found in [11].
These waves are oscillatory parallel to the x1 axis and decay exponentially along the
x2 axis in the upper half plane R3. The higher the spatial frequency of the evanes-
cent plane waves used to probe the scatterer is, the more rapidly the field decays as a
function of depth into the scatterer. See Figure 2.1 and 2.2 for examples. Evidently,

such incident waves satisfy the homogeneous equation
na+fia=a 913
The total electric field 90 consists of the incident field <20 and the scattered field w:

¢=¢0+1lk

43

 

 

_1 0 X2

Figure 2.1. Evanescent plane wave at k0 = 4.0 with 77 = 4.7.

44

 

-10 x2

Figure 2.2. Evanescent plane wave at k0 = 4.0 with 77 = 8.0.

45

It follows from the equations (2.1.1) and (2.1.3) that the scattered field satisfies
All + k3(1 + (1)10 = -kgq¢o- (2-1-4)

Remark 2.1.1. In this paper, we adopt the non—global approach, i.e. the scattered
field resulting from the interaction of the incident field with the scatterer is analyzed
in the absence of other medium or tip. The scattering problem may be formulated in
the free space. The global approach which takes into account the entire system is the

subject of our ongoing research.

In the free space, the scattered field is required to satisfy the following Sommerfeld
radiation condition

r 13,1th (5; —1k0w) — 0, r — Ix],

uniformly along all directions x / Ix]. In practice, it is convenient to reduce the problem
to a bounded domain by introducing an artificial surface. Let D 2 [—L1, L1] x [0, L2]
be a square, which contains the compact support of the scatterer, 9. Let 8D be
the boundary of D. Denote n the unit outward normal to 8D. A suitable boundary
condition needs to be imposed on 8D. For the sake of simplicity, we employ the first

order absorbing boundary condition [22] as

('9 .
8—1: —1k0w = 0, on 3D. (2.1.5)

Given the incident field (to, the direct problem is to determine the scattered field
1b for the known scatterer q(x). Using the Lax-Milgram lemma and the Ffedholm
alternative, the direct problem is shown in this paper to have a unique solution for all
100 > 0. An energy estimate for the scattered field is given, which provides a criterion

for the weak scattering. Furthermore, properties on the continuity and the Fréchet

46

differentiability of the nonlinear scattering map are examined. For the regularity
analysis of the scattering map in an open domain, the reader is referred to [2, 24, 10].
The inverse medium scattering problem is to determine the scatterer q(x) from the
measurements of near-field currents densities, wlap, given the incident field 050. The
inverse medium scattering problems arise naturally in diverse applications such as
radar, sonar, geophysical exploration, medical imaging, and nondestructive testing
[10]. However, there are two major difficulties associated with these inverse problems:
the ill-posedness and the presence of many local minima. In' [8, 4], stable and efficient
continuation methods with respect to the wavenumber were proposed to solve the
two-dimensional Helmholtz equation and the three-dimensional Maxwell’s equations,
respectively, in the case of full aperture data. A homotopy continuation method with
limited aperture data may be found in [3] These approaches require multi-frequency
scattering data and are based on recursivellinearization along wavenumbers.

The main purpose of this paper is to study the inverse medium problem for
Helmholtz’s equations at a single-frequency. We present a new continuation method
for the inverse medium scattering problem. In the case of radially symmetric scatter—
ers, Chen in [9] developed a recursive linearization algorithm with single-frequency
data, where spherical incident waves were used. In this paper, we attempt to remove
the radially symmetric assumption on the medium. Our approach is motivated by the
recent studies of near-field optics. As a special feature, the illuminating fields used
in this paper including the high spatial frequency evanescent plane waves are a one-
parameter family of plane waves. When a medium is probed with an evanescent plane
wave at a high spatial frequency, only a thin layer of the medium is penetrated. Cor-
responding to this exponentially decaying incident field, the scattered field measured
on the boundary contains information of the medium in that thin layer. Although
such a measurement is entirely inadequate to determine the whole medium, it does

give rise to an approximation. To accurately determine the medium, information

47

at lower spatial frequencies of the evanescent plane waves is needed to illuminate
the medium. While the probing field penetrates a thicker layer of the medium, the
relation between the measurement and the scatterer to be recovered in the thicker
layer becomes more nonlinear. These nonlinear equations can be considered as per-
turbations to the already solved equations in the previous thicker layers. Therefore,
they can be continually and recursively linearized by a standard perturbation tech-
nique. Thus, the recursive linearization is a continuation method along the transverse
direction of the incident wave, which controls the depth of its penetration.

The plan of this paper is as follows. The analysis of the variational problem
for direct scattering is presented is Section 2.2. In particular, the well—posedness of
the direct scattering is proved. The Fréchet differentiability of the scattering map
is also given. In Section 2.3, an initial guess of the reconstruction from the Born
approximation is derived in the case of weak scattering. Section 2.4 is devoted to
numerical study of a regularized iterative linearization algorithm. Numerical examples
are presented. The paper is concluded with some general remarks and directions for

future research in Section 2.5.

2.2 Analysis of the Scattering Map

In this section, the direct scattering problem is studied to provide some criterion for
the weak scattering, which plays an important role in the inversion method. The
Fréchet differentiability of the scattering map for the problem (2.1.4), (2.1.5) is ex-

amined.

Remark 2.2.1. Some analysis of the scattering map was given previously by Keys
and Weglein [24] based on the integral equation approach and contraction mapping
theorem. The assumption of small perturbation of the potential is necessary for their

approach. Our approach is diflerent. Based on the Fredholm alternative and a unique-

48

ness result, we develop a variational approach to prove the existence of the scattered
field for all k0 > 0, given q E L°°(D), the continuity of the scattering map, the bound-
edness of the formal linearized map, and the Fréchet difierentiability of the scattering
map. The assumption of small perturbation is not needed in our analysis. More im-
portantly, we give an explicit energy estimate for the scattered field, which provides
a criterion for weak scattering hence plays a central role in the development of the
inversion algorithm of Section 2.3. An analysis of the Fréchet difierentiability on the
scattering map for the equation (2.1.4) along with the Sommerfeld radiation condition

may also be found in [2] using the integral equation approach.

To state our boundary value problem, we introduce the bilinear form a : H 1( D) x
H1(D) —> (C
a(u, v) = (Vu, Vv) -— kg((1 + q)u, v) — ik0(u,v),

and the linear functional 011 H1(D)

bf”) = (“601490, "Ul-

Here, we have used the standard inner products

(u,v)=/u-deand(u,v)=/ u-fids,
D a!)

where the overline denotes the complex conjugate.
Then, we have the weak form of the boundary value problem (2.1.4) and (2.1.5):
find ll) E H1(D) such that

a(t/J,€) = 0(5), vg 6 111(0). (2.2.1)

Throughout the paper, the constant C stands for a positive generic constant whose

49

value may change step by step, but should always be clear from the contexts.
For a given scatterer q and an incident field (to, we define the map S (q, (750) by
1,1) = S (q, (750), where 112 is the solution of the problem (2.1.4), (2.1.5) or the variational
problem (2.2.1). It is easily seen that the map S (q, 030) is linear with respect to (750
but is nonlinear with respect to q. Hence, we may denote S (q, (to) by S (q)q§o.
Concerning the map S (q), we have the following regularity results. Lemma 2.2.3
gives the boundedness of S (q), while a continuity result for the map S (q) is presented

in Lemma 2.2.4.

Lemma 2.2.1. Given the scatterer q E L°°(D), the direct scattering problem (2.1.4),

(2.1.5) has at most one solution.

PROOF. It suffices to show that w = 0 in D if 920 = 0 (no source term). From the
Green’s formula

_‘_ _ G_ fa. I . ' ‘
0 = / (wA'lp — wAvir)dx = / (lit-i) — w—u—jhls = —21k0/ |w|2ds,
D at) 671 a" an

we get w = 0 on 8D. The absorbing boundary condition on (9D yields further that

6
77): = 0 on 8D. By the Holmgren uniqueness theorem, 11.1 = 0 in R2 \ D. A unique
"n

continuation result [21] concludes that w = 0 in D. C]

Lemma 2.2.2. If the wavenumber k0 is sufiiciently small, the variational problem
(2.2.1) admits a unique weak solution in H1(D) and S (q) is a bounded linear map

from L2(D) to H1(D). Furthermore, there is a constant C dependent of D, such that

)1 saw. 11 H1( D): or. H q 130(0)) to limp, . (2.2.2)

50

PROOF. Decompose the bilinear form a into a = a1 + kgag, where

611(4)!) = (Vi/1370 41200146).
LIN/1,6): -((1+ (IMO-

We conclude that al is coercive from

2 0130 II V ”in1”.

where the last inequality may be obtained by applying standard elliptic estimates
[17]. Next, we prove the compactness of (1.2.. V Define an operator A : L2(D) -—+ H1(D)
by

amino = art/1,5). V6 6 111(1)),

which gives
(Wu), Vt) — MAM) = —((l + «005). V6 6 H1(D).
Using the Lax—Milgram Lemma, it follows that
“Aw” k911i“ (223)
H1(D)— k0 L2(D)’ ' '

where the constant C is independent of k0. Thus A is bounded from L2(D) to H1(D)
and H1(D) is compactly imbedded into L2(D). Hence A : L2(D) —> L2(D) is a

compact operator.

51

Define a function u E L2(D) by requiring u E H1(D) and satisfying
a1(u,§) = b(€), V{ E H1(D).

It follows from the Lax—Milgram Lemma again that

H U “H1(D)S Cko H q ||L°°(D)H (J50 ”L2(D) - (224)

Using the operator A, we can see that the problem (2.2.1) is equivalent to find it E
L2(D) such that
(I + kgA)t/) = 11. (2.2.5)

When the wavenumber k0 is small enough, the operator I + kgA has a uniformly

bounded inverse. We then have the estimate

'1 ill ”L2(D)S C H u “L2(D)’ (2-2-6)

where the constant C is independent of k0. Rearranging (2.2.5), we have 17’) = u —

kgAt/J, so it E H1(D) and, by the estimate (2.2.3) for the operator A, we have

I] 76 ”H1(D)S“ u “H1(D) +0130 II '16 “L2(D)'

The proof is complete by combining the estimates (2.2.6) and (2.2.4) and observing
that 712 = S(q)<750. C]

For a general wavenumber k0 > 0, from the equation (2.2.5), the existence follows
from the Fredholm alternative and the uniqueness result. However, the constant C

in the estimate (2.2.2) depends on the wavenumber.

Lemma 2.2.3. Given the scattererq E L°°(D), the variational problem (2.2.1) admits

52

a unique weak solution in H1(D) for all k0 > 0 and S (q) is a bounded linear map

from L2(D) to H1(D). Furthermore, it holds the estimate

1] S(‘lMO “H1(D)S C H q “L°°(D)” (DO ”L2(D)’ (2-2-7)

where the constant C depends on kg and D.

Remark 2.2.2. It follows from the explicit form of the incident field (2.1.2) and the

estimate (2.2.7) that

 

 

 

is» s CIQI‘” n q ”Loom,

111(0)

where Q is the compact support of the scatterer q and the constant C depends on

k0, D. Moreover, we have for [77] > k0 that

H 11> ll H1( 0,: 0(7)? — kill-V4 II 0 “Loom (228)

where the constant C depends on kg and D.

Remark 2.2.3. The estimate of the scattered field in Remark 2.2.1 provides a cri-
terion for the weak scattering. For a fixed wavenumber kg and a scatterer q, the

scattered field is weak if the spatial frequency of the incident wave, |77|, is large.

Lemma 2.2.4. Assume that q1,q2 E L°°(D). Then

ll S(qlwo - sumo u H1( D): C H q. — q. 11mm” it 152(1),, (229)

where the constant C depends on kg, D, and [I q2 ]]L00(D).

PROOF. Let 1111 = S(q1)cbo and 1132 = S(q2)qbo. It follows that for j = 1,2

53

By setting w = tbl — tbg, we have

Aw + 190(1 + (11)“) = —k3(91 — Q2)(¢0 1' ‘62).

The function w also satisfies the boundary condition (2.1.5).

We repeat the procedure in the proof of Lemma 2.2.3 to obtain

1] w ”H1(D)S C I] (11 — (I2 HL°°(D)H <00 +113’211L2(D) .

Using Lemma 2.2.3 again for 1112 yields

1] "(’2 HH’(D)S 01102 llL°°(D)H €50 ”L2(D)’

which gives

”5010650 — 50121650 “H1(D)S C H (11 — Q2 l|L°°(D)H €50 ”L2(D)’

where the constant C depends 011 D, kg, and I] qg [I Loo ( D)‘ D
Let y be the restriction (trace) operator to the boundary 8D. By the trace
theorem, 7 is a bounded linear operator from H1(D) onto H1/2(6D). We can now
define the scattering map It! (q) = 75 (q)
Next, consider the Pféchet differentiability of the scattering map. Recall the map
S (q) is nonlinear with respect to q. Formally, by using the first order perturbation
theory, we obtain the linearized scattering problem of (2.1.4), (2.1.5) with respect to

a reference scatterer q,

Av + 1:3(1 + q)v = —t:g(5q(¢0 + a), (2.2.10)
8v _
57—; —1k0v = 0, (2.2.11)

54

where 16 = S(q)¢0.

Define the formal linearzation T (q) of the map S (q) by v = T(q)(6q, (750), where
v is the solution of the problem (2.2.10), (2.2.11). The following is a boundedness
result for the map T (q) A proof may be given by following step by step the proofs

of Lemma 2.2.2 and Lemma 2.2.3. Hence we omit it here.

Lemma 2.2.5. Assume that q,6q E L°°(D) and do is the incident field. Then v =
T(q)((5q, do) E H1(D) with the estimate

H T(q)(5q,¢0) IIHI(D)S C 1] (Sq [ll/30(1))“ <50 “L2(D)’ (2-2-12)

where the constant C depends on k0, D, and [I q ||Loo(D).
The next lemma is concerned with the continuity property of the map.

Lemma 2.2.6. For any q1,q2 E L°°(D) and an incident field 020, the following esti-

mate holds

(2.2.13)

where the constant C depends on k0, D, and [I q2 ]]LOO(D).

PROOF. Let v,- = T(q,)(6q, 450), for i = 1, 2. It is easy to see that

AWI — ”2) + A712)“ + €11)(V1 - ’02) =

_ k65qwi — 1V2) — 193011 — Q2)V2.

where ”(bi = S ((17100-

Similar to the proof of Lemma 2.2.3, we get

|| 0102— 1) 1,1,0): 0 (II 5a “Loom“ it — 112|IH1(D,+Hq1— (12 “30(0)” ||H1(D))-

55

From Lemma 2.2.2 and Lemma 2.2.3, we obtain

ll vi - v2IIH1(D)S C H (11 - (12 “L°°(D)“ 5Q ”L°°(D)” (I50 “L2(D)’

which completes the proof. [:1

The following result concerns the differentiability property of S (q)

Lemma 2.2.7. Assume that q, (Sq E L°°(D). Then there is a constant C dependent

of k0, D, and [I q ]]L00(D), for which the following estimate holds

H 5(9 + 5(000 — 506950 — T(CIWSCI) <90) “H1(D)S 01169112 00(1))“ $0 ”L2(D)'
(2.2.14)

Proof. By setting w] = S(q)cj')0, v.12 = S(q + (Smog, and v = T(q)(6q, (110), we have

A991 + (33(1 + Q)'l+‘/l1 = —li’(2)(1¢o,
A1112 +l'3(1+ q + 5mm = 4301 + 50m).

A?) + k3<1+ Q)’U '2 —k(2,6q1f)1 — ligtsquo.

In addition, 1/11, 1112, and v satisfy the boundary condition (2.1.5).

Denote U = tbg — 7,01 — v. Then
AU '1' 196(1 + (DU = —k735Q(’ll-’2 — 4’1)-

Similar arguments as in the proof of Lemma 2.2.3 give

”UHH1(D)S C I] 6‘1 HL°°(D)” 762 ‘ ill] “H1(D)'

56

From Lemma 2.2.3, we obtain further that

IIUIIH1(D)S C II 5C] [Ii-100(1))” $0 ”122(0) .

Finally, by combining the above lemmas, we arrive at

Theorem 2.2.1. The scattering map M (q) is Fréchet differentiable with respect to q
and its Fréchet derivative is

D.M(q) = 7T(q). (2.2.15)

2.3 Inverse Medium Scattering

In this section, a regularized recursive linearization method for solving the inverse
medium scattering problem of the Helmholtz equation in two dimensions is proposed.
The algorithm, obtained by a continuation method on the spatial frequency of a 011e-
parameter family of incident plane waves, requires only single-frequency scattering
data. At each transverse part of the incident wave, the algorithm determines a forward
model which produces the prescribed scattering data. Since the incident wave at a
high spatial frequency can only penetrate a thin layer of the scatterer, the scattered
field is weak. Consequently, the nonlinear equation becomes essentially linear, known
as the Born approximation. The algorithm first solves this nearly linear equation at
the largest I77] to obtain an approximation of the scatterer. This approximation is then
used to linearize the nonlinear equation at the next smaller spatial frequency of the
incident wave, which can penetrate a thicker layer of the scatterer, to produce a better
approximation. When the spatial frequency, [77], is smaller than the fixed wavenumber
k0, the incident wave becomes usual propagating plane wave, and the whole scatterer

is illuminated. This process is continued until the spatial frequency is zero, where the

57

approximation of the scatterer is considered as the final reconstruction.

2.3.1 Born Approximation

Rewrite (2.1.4) as
All + k 01/)— _ —k3q(¢o + a). (2.3.1)

Consider a test function 7110 = eikO‘T ' d, d: (cos 6, sin 6), 6 E [0, 2n]. Hence we satisfies
(2.1.3).
Multiplying the equation (2.3.1) by 1,190, and integrating over D on both sides, we

have

/ 'tfioAu’idx + 113/ wowdx = —k(2, / (1(920 +u’2)v’10dx.
D D D

Integration by parts yields

, (9d) 811’”, I ' '
/ @I’All’odiv + / (ii/’OT— — Til—ON £18 + k3] 11.1017de = —kg/ q(cfm + 11))u’0dx.
D OD 0” an) D D

We have by noting (2.1.3) and the boundary condition (2.1.5) that

1 0
/ (1(00 ‘1’ w>w0d$= kj/D if) (Eli—O —lli‘01,l’0) d8.
D

Using the special form of the incident wave and the test function, we then get

/ q($)ei(n + kO COS 9)$1ei(k(77) + [6081116)3326133 : i— 1601' J_ 1)elk01: ' (ids
D

60
“ / (1164506117-
D

(2.3.2)

From Lemma 2.2.3 and Remark 2.2.2, using an evanescent incident plane wave at

a high spatial frequency, the scattered field is weak and the inverse scattering problem

58

becomes essentially linear. Dropping the nonlinear (second) term of (2.3.2), we obtain

the linearized integral equation

. _ /. 2_ 2 - - . _, . ~
[q(x)el(77+k0C086)$16( ’l k0+lk0Sln6)$2dx i w(n°d—I)€lk0$.dd3,
D

= k0 80
(2.3.3)

which is the Born approximation.

Since the scatterer q(x) has a compact support, (2.3.3) can be rewritten as

/L2(j(€312)e(fl‘/722 _ kg + iko sin 0>$2dx2 ___ i 11/)(71 . J— 1)eik$ . dds
0 0 BD
where f = r) + ko cos6 and a(g, x2) is the Fourier transform of q(x) with respect to
.131. When the spatial frequency I77I is large, the incident wave penetrates a thin layer
of the scatterer. Tints, the Born approximation allows a reconstruction containing
information of the true scatterer in that thin layer. In [8, 4], the inversion involves data
related to the scatterer through the Fourier transform in the case of weak scattering.
Here, due to the presence of the evanescent wave, the inversion involves data related
to the scatterer through a Fourier (with respect to x1)-Laplace (with respect to $2)
transform in the case of the weak scattering. Since the inversion of the Laplace
transform is ill-posed, we consider the Landweber iteration to implement the linear
integral equation (2.3.3) in order to reduce the computation cost and instability [29].
Define the data
i

10). 6) = k0
0 for I711 < "max,

/ Wn - d— 1)eik$ ' dds for In] 2 Timax,
oD

where ”max is some large positive number.

59

The integral equation (2.3.3) can be written as the operator form

AM, 9;:I:)q($) = f (17, 9)- (2-3-4)

Following the idea of the Kaczmarz method, we use partial measurement data in-
stead of using all them simultaneously for each sweep. Let 77,,i = 1, ...,I, be the

discretization of 77, where I is the number of sweeps. Then we can rewrite (2.3.4) as

A(77,-,6;x)q(x) = f(77,-,6), i = 1, ...,I,

or in short

Aiq=f,', i=1,...,1.

For each sweep i, the Landweber iteration takes the form

0;” = q;'-1>+ aAgu. — Amy—"n, 1e 8.

o
1

where a is a relaxation parameter. Since we just need an initial guess for the iteration
in the recursive linearization, we only take one step Landweber iteration for each

sweep, which yields
q,- = q..,-_1+ aA;(f,- — A(q,:_1)), i = 1, ..., 1, (2.3.5)
where q] is used as the starting point of the following recursive linearization algorithm.

2.3.2 Recursive Linearization

As discussed in the previous section, when the spatial frequency I77| is large, the Born

approximation allows a reconstruction of the thin layer for the true scatterer. In this

60

section, a regularized recursive linearization method for solving the two-dimensional
Helmholtz equation at fixed frequency is proposed.

Choose a large positive number ”max and divide the interval [0, Umaxl into N
subdivisions with the endpoints {770, 771, ..., 7710}, where 770 = 0, 771v =-- ”max, and
77,- _ 1 < 77,; for 1 S i S N. We intend to obtain qn recursively at 77 = 771v, 7710-1, ..., 770.

Suppose now that the scatterer qfi has been recovered at some 77 = 77,- + 1 and that
77 = 77,- is slightly less than 7'). We wish to determine qn, or equivalently, to determine

the perturbation

50 = (177 — (1;,-

For the reconstructed scatterer q,~,, we solve at the spatial frequency 77 the forward

scattering problem

 

A620”) + 1:3(1 + (1.,*,)'l;’(j’ l) = —k8q,~,c‘)((,]‘2), (2.3.6)
3,7,,(111') , - - -
an — 11601.0()" 2) = 0, (2.3.7)
where the incident wave 0653’ Z) = eilli‘l’l + i“7(1):”, I j] Z i.

For the scatterer q,,, we have

 

A1120”) + k‘3(1 + (177111104) = —kgq,3§,]’ 2), (2.3.8)
Jan —1k0w(]’7‘) = 0. (2.3.9)

Subtracting (2.3.6), (2.3.7) from (2.3.8), (2.3.9) and omitting the second-order small-

ness in 6g and in 61/10) = 71/2012) — do” 7’), we obtain

 

A330) + k3(1+ (1,932,001 —_- —1:33q(¢§j’ 7") + 11303 '0), (2.3.10)
(1') .
36:” 41.052110) = 0. (2.3.11)

61

For the scatterer qn and the incident wave (153] ’ i), we define the map 5,-(6177, <65], 3))
by

SJ(CIn,¢c()j’i)) = TAU”),

where $01 i) is the scattering data corresponding to the incident wave ((9th ’ 2). Let 7

be the trace operator to the boundary 8D. Define the scattering map

114110.05”): 751(0)),29’20

(1» 73)

For simplicity, denote M,(q77, Q50 ) by ll’I,(q77). By the definition of the trace oper-
ator, we have

nan) = W ”in.

Let D11! ,(qfl) be the Fréchet derivative of lib-(rm) and denote the residual operator
by

Rj(Qf]) :¢1(J17')IOD _ 7,3(J Z)I0D'

It follows from Theorem 2.2.1 that

Similarly, in order to reduce the computation cost and instability, we consider the

Landweber iteration of (2.3.12), which has the form

(Sq = fiDlW,’-'(q,~,)R,-(q,~,), for all IjI 2 i, (2.3.13)

where B is a relaxation parameter and DM,‘-‘(q.;,) is the adjoint operator of DM,(q,~,).
In order to compute the correction 6g, we need some efficient way to compute

DM,*‘(qfi)R,-(q,~,), which is given by the following theorem.

62

Theorem 2.3.1. Given residual R,(q,7), there exits a function (1501i) such that the

adjoint Fréchet derivative DM,i(q,~,) satisfies

 

 

[DM,:(q,-,)R,(q,,)] (x) = 163(3),?" 0 (x) + 3(1,i)(x))¢(ja 21(3), (2.314)

where ¢(()],z) is the incident wave and 62012) is the solution of (2.3.6), (2.3.7) with
(332')

the incident wave 0‘10 .

PROOF. Let 173(1) 2) be the solution of (2.3.6), (2.3.7) with the incident wave 63]”).

Consider the following problem

 

Add“) + 1:3(1 + q,~,)6¢’~(j) = _k36q(<b,(,]’ "l + 160’ 21), (2.3.15)
06 (j) , , °
5” — 11031.9(] ) = 0. (2.3.16)

and the adjoint problem

 

A3031) + 13(1 + q,~,)o(j’ 7") = 0, (2.3.17)
005(1) i) , , - -
on + 11.6001 2) = R,(q,~,). (2.3.18)

Since the existence and uniqueness of the weak solution for the adjoint problem may
be established by following the same proof of Lemma 2.2.2, we omit the proof here.
Multiplying the equation (2.3.15) with the complex conjugate of (19013) and inte-

grating over D on both sides, we obtain

 

 

 

/ 3(1»11A63(J’)dx+1:3/ (l+q,~,)6¢)(j)q§(j~.ildx= 413/ 6q(¢(()j’i)+1l‘(j"i))¢(jiildx.
D D D

63

Integration by parts yields

 

 

[0190,50, 0663:) _ 510(1)3¢;:1))d8 2 _kg/qunga) + d3(j,i))¢(j,i)dx.

 

Using the boundary condition (2.3.16), we deduce

 

 

 

/ 6,0)(3350 _,,,0,,(j,1))d,= kg / W331)1,36,25,60,...
at) n D

It follows from (2.3.12) and the boundary condition (2.3.18) that

 

[a lDll’lj((I.,7)5(1]33(0)?)618 = 63/ 6q(c’9((,J’l) + .,,T,(j.i))¢(j,i)d$.
'D D

We know from the adjoint operator D11’I,"(q,r,) that

 

 

j(quM;(q,,)Rj(q,~,)dx=kg/ 6q(@3j,'i)+,;,(1.i))¢(j. 0,133,
o D

Since it holds for any dq, we have

 

 

DM;<q.~,)R.-<q~) = 142,653”) + 00* (1)601 2').

Taking the complex conjugate of the above equation yields the result.

Using this theorem, we can rewrite (2.3.13) as

 

 

(Sq = kg,6(¢,(,j’ 2') + 1,003 i))d>(j’ 6. (2.3.19)

So for each incident wave with a transverse part 77,-, we have to solve one forward

problem (2.3.6), (2.3.7) along with one adjoint problem (2.3.17), (2.3.18). Since the

adjoint problem has a similar variational form as the forward problem. Essentially, we

need to compute two forward problems at each sweep. Once (Sq is determined, qfi is

64

 

Initialization

771v = nmax largest nmax
q'rmax Born approximation
Reconstruction loop:
FOR i = N : 0 (77, = ”max : 770) march along spatial frequency
FOR 7 = N : i(|77,-I = ”max : 772-) perform refinement

solve (2.3.6), (2.3.7) for ($013) one forward problem
solve (2.3.17), (2.3.18) for 3(3) 2) one adjoint problem

66' = 333037 7') + new 2)
3’: == at + 5,7
END .
Z

(17' 3: (11'
END

q :2 qO final reconstruction

 

 

Table 2.1. Recursive linearization reconstruction algorithm.

updated by q,~, + dq. After completing sweeps with [77,-I Z 77, we get the reconstructed

scatterer (177 at the spatial frequency 77.

Remark 2.3.1. For given 77,, iterations for [77,-I Z 77,- could be repeated to improve
the accuracy of the approximation for qui' However, in practice, this refinement is
usually unnecessary because of the slow convergence of the Landweber iteration at the
same stage [14], i.e. without using essentially difierent data. Numerical results show
that the iterative process described as the reconstruction loop in Table 2.1 is suflicient

to obtain reasonable accuracy.

The recursive linearization for inverse medium scattering at fixed frequency is

summarized in Table 2.1.

65

2.4 Numerical Experiments

In this section, we discuss the numerical solution of the forward scattering problem
and the computational issues of the recursive linearization algorithm.

The scattering data are obtained by numerical solution of the forward scattering
problem. As for the forward solver, we adopt the finite element method (FEM), which
leads to a sparse matrix. The sparse large-scale linear system can be most efficiently
solved if the zero elements of coefficient matrix are not stored. We used the commonly
used compressed row storage (CRS) format which makes no assumptions about the
sparsity structure of the matrix and does not store any unnecessary elements. In fact,
from the variational formula of our direct problem (2.2.1), the coefficient matrix is
complex symmetric. Hence, only the lower triangular portion of the matrix needs
be stored. Regarding the linear solver, either biconjugate gradient (BiCG) or quasi-
minimal residual (QMR) algorithms with diagonal preconditioning may be used to
solve the sparse, symmetric, and complex system of the equations. For our examples,
it appears that the QMR is more efficient.

In the following, to illustrate the performance of the algorithm, three numerical
examples are presented for reconstructing the scatterer of the Helmholtz equation in
two dimensions. For stability analysis, some relative random noise is added to the

date, i.e. the electric field takes the form

wIaD :2 (1+ orand)w|30.

Here, rand gives uniformly distributed random numbers in [—1, 1] and a is a noise level

parameter taken to be 0.02 in our numerical experiments. The relaxation parameter

66

fl is taken to be 0.01. Define the relative error by

(Em lqz'j - Quiz)”2

(22', j lqrijlzll/2

 

’

where (j is the reconstructed scatter and q is the true scatterer.

Example 1. Let

v 2 2 2 2
q(r1,:1:2) = 0.3(1— 3:1)26‘351—(552 +1) _ (fl __ $3 _ mike—($1 + 272)

5 l
—ie—($1 + ”2 — 333,
30
reconstruct a scatterer defined by
(1101.122) = (1(3921, 3(272 — 1)) (2.4.1)

inside the domain D = [—1, 1] x [0, 2]. See Figure 2.3 for the surface plot of the
scatterer function. Figure 2.10 is the final reconstruction using the wavenumber
k0 = 10.0 and the step size of the spatial frequency 617 = 0.6. Figures 2.4-2.9 show
the evolution of reconstructions at different spatial frequencies. Figure 2.11 presents
the effect of the wavenumber he on the result of reconstruction, which illustrates
clearly that the inversion using a larger wavenumber k0 is better than that using a
smaller one. This result may be explained by Heisenberg’s uncertainty principle [8, 9].
Figure 2.12 shows the relative error by using different step size of the spatial frequency,
which suggests that we may use a large step size in order to save computation cost

since the final reconstruction is not really sensitive to the step size.

67

 

 

Figure 2.3. Example 1: true scatterer ql.

68

, t a 9' .
. $““‘ 0..
’ .m-‘with‘li“i“l\‘3‘$’h’i’
t-1‘l\‘QXo‘h‘b‘:\“:\‘:“‘9‘99?!
k-\‘\&\\\“\‘:ts‘\§ \‘g‘h if

\\ \\\\\\‘ ‘6

\ ‘ \
\\§{\‘R‘:“ ‘
\

\ ‘ \\
i . I
\\ \\\\\ V§§QT
‘ k k i" ii “ . .

 

 

Figure 2.4. Example 1: Born approximation.

69

 

 

Figure 2.5. Example 1: reconstruction of ql at 17 = 10.2.

70

.::.

‘33“? ’9’,

‘3\\\“

 

 

Figure 2.6. Example 1: reconstruction of ql at 77 = 8.4.

71

“e.
‘\
I 20.06.:3‘

’ was s
i ’r’o‘q‘fix

8% a

a ‘ ‘\;\\ \I‘Q-Kfs ”5], I.
‘ a if; ,“ ;

’;{\
Jxfifig ~
.u

 

 

Figure 2.7. Example 1: reconstruction of ql at 17 = 6.6.

72

“Eh“ \\\\\\\:

Ilo’u:,7y\\\\\\\us

 

 

Figure 2.8. Example 1: reconstruction of q1 at 77 = 4.8.

73

 

 

Figure 2.9. Example 1: reconstruction of ql at 77 = 3.0.

74

   

       

 

0.8\

”l " i it uh

0.4\ [I'M .““‘:‘“§§§ v“.
'“

 

Figure 2.10. Example 1: reconstruction of ql at 77 = 0.0.

75

 

Relative error

 

 

 

l l l l L

. 0
12 1O 8 6
Spatial frequency

Figure 2.11. Example 1. Relative error of the reconstruction ql at different wavenum—
ber kg. 0: reconstruction at k0 = 10.0; *: reconstruction at k0 = 8.0; E]: reconstruc-

tion at k0 = 6.0; +: reconstruction at k0 = 4.0.

76

 

Relative error

 

 

 

l l l l l

0
1 2 1 0 8 6 4 2 0
Spatial frequency

Figure 2.12. Example 1. Relative error of the reconstruction ql at different step size
677. o: reconstruction at (577 = 0.6; *: reconstruction at 677 = 1.2; D: reconstruction
at 577 = 2.0.

77

Example 2. Reconstruct a scatterer defined in D by

(11(231/0-8, 132/08) for ac? + (x2 — 1)2 5 0.7472,
q2(:c1, x2) = -0.3, for 0.7472 < as? + (1:2 — 1)2 3 0.8532,

0, for 113% + (3:2 — 1)2 > 0.8532.

See Figures 2.13 and (2.15) for the surface and contour plots of the function. It is
easily seen that this scatterer is difficult to reconstruct because of the discontinuity
across two circles. The example could be regarded as a model problem for ultrasound
tomography of a human head, where the skull is represented by the thin layer of
denser material in the narrow annulus region. Figures 2.14 and 2.16 Show the surface
and contour plots of the reconstructed scatterer using the wavenumber k0 = 15.0
and the step size 677 = 0.85. Figure 2.17 gives the evolution of reconstruction hor—
izontally across the diameter. An examination of the plots shows that the error of
the reconstructions occurs largely around the discontinuities, while the smooth part
is recovered more accurately. As expected, the Gibbs phenomenon appears in the
reconstructed scatterer near the discontinuity.

Example 3. Reconstruct a scatterer defined in D by

cos(2.57rr1) for r1 3 0.2,
(13(1131, 1‘2) = eos(2.57rr2) for r2 3 0.2,

0 otherwise,

 

 

where r1 = ,/(:c1+ 0.25)2 + (272 —1.0)2 and r2 = \/(;1:1 - 0.25)2 + (.732 —1.0)2. The
compact support of this scatterer is two isolated disks with the same radius of 0.2 and
the centers at (—0.25, 1.0) and (0.25,1.0). See Figures 2.18 and 2.20 for the surface
plot and image of the function. Figures 2.19 and 2.21 are the final reconstruction

using the wavenumber k0 = Br and the step size of the spatial frequency 677 = 0.6.

78

.. ’7 ‘ \R.
.l/ l. I b “‘“M
”W 1111““:
"lllll‘llxtt ll:
mile“

/ I
ally! '
”lilo?

{’9‘

 

 

Figure 2.13. Example 2: true scatterer of q2.

79

\\\l\\\\\{
l 1‘ \\\\\ll \ r.
li\\\\\\\ll\li“

(17:11 ”‘1.
”4.1.1 \

 

 

Figure 2.14. Example 2: reconstruction of (12.

80

 

 

 

 

 

 

2

1.5-

 

 

0.5 -

 

 

 

 

Figure 2.16. Example 2: contour view of the reconstructed scatterer (72.

81

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 2.17. Example 2. Evolution of slice for the reconstruction q2. Solid line:
true scatterer; Circle: reconstruction. Top row from left to right: reconstruction at
77 = 14.45; reconstruction at 77 = 13.60; reconstruction at 77 = 12.75; middle row from
left to right: reconstruction at 77 = 10.20; reconstruction at 77 = 8.50; reconstruction
at 77 = 6.80; bottom row from left to right: reconstruction at r7 = 5.10; reconstruction
at 77 = 2.55; reconstruction at 77 = 0.0.

82

12R

 

 

Figure 2.18. Example 3: true scatterer of q3.

This example is used to examine the resolution of the reconstructed image. In this
numerical experiment, the wavelength of the incident plane waves is 27r/k0 = 0.6.
The distance of the centers for the compact support is 0.5, which is less than one
wavelength. From the well separated bumps, the resolution of the image is clearly in
the scale of subwavelength. The subwavelength resolution is expected since evanescent

waves are used for illumination.

83

1.2\

1\

0.8 x,

0.6\

b“ ,1

  

0.4 a

AIIWWLL.
I «I I‘ " '

 

A ..);_~___

. ““:" v.1-

. ‘1‘; lFtI'Il'l “‘4‘.‘l 4 l‘ ‘ ‘
n . ‘ ff

- n..
'i.3l}_' 1‘ ‘HI... 1
“unwindin-
W' '

 

Figure 2.19. Example 3: reconstructed scatterer of q3.

84

 

 

Figure 2.20. Example 3: image view of the true scatterer q3.

85

 

 

Figure 2.21. Example 3: image view of the reconstructed scatterer q3.

86

2.5 Concluding Remarks

We have presented a new continuation method with respect to the spatial frequency of
a one-parameter family of plane waves. The recursive linearization algorithm is robust
and efficient for solving the inverse medium scattering at fixed frequency. Finally, we
point out some future directions along the line of this work. The first is concerned
with the convergence analysis. Although our numerical experiments demonstrate the
convergence and stability of the inversion algorithm, no rigorous mathematical result
is available at present. Another direction is to investigate inverse medium problems
for Maxwell’s equations at fixed frequency. We are currently attempting to extend the
approach in this paper to the more complicated 3D model problems and will report

the progress elsewhere.

87

BIBLIOGRAPHY

[1] H. Ammari and G. Bao, Analysis of the scattering map of a linearized inverse

medium problem for electromagnetic waves, Inverse Problems, 17(2001), 219-234.

[2] G. Bao, Y. Chen, and F. Ma, Regularity and stability for the scattering map of a
linearized inverse medium problem, J. Math. Anal. Appl., 247(2000), 255—271.

[3] G. Bao and J. Liu, Numerical solution of inverse scattering problems with multi-
experimental limited aperture data, SIAM J. Sci. Comput., 25(2003), 1102-1117.

[4] G. Bao and P. Li, Inverse medium scattering problems for electromagnetic waves,

SIAM J. Appl. Math, to appear.

[5] G. Bao and P. Li, Inverse medium scattering for Helmholtz equation at fired

frequency, Inverse Problems, submitted.

[6] P. Carney and J. Schotland, Three-dimensional total internal reflection mi-
croscopy, Opt. Lett., 26(2001), 1072-1074.

[7] P. Carney and J. Schotland, Near-field tomography, MSRI Series in lV’Iathematics
and Its Applications, 47(2003), 131-166.

[8] Y. Chen, Inverse scattering via Heisenberg uncertainty principle, Inverse Prob-
lems, 13(1997), 253-282.

[9] Y. Chen, Inverse scattering via skin effect, Inverse Problems, 13(1997), 649-667.

[10] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory,
2nd ed., Appl. Math. Sci. 93, Springer-Verlag, Berlin, 1998.

[11] D. Courjon, Near-field Microscopy and Near-field Optics, Imperial College Press,
2003.

[12] D. Courjon and C. Bainier, Near field microscopy and near field optics, Rep.
Prog. PhyS., 57(1994), 989-1028.

88

[13] O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, A nonlin-
ear inversion method for 3D electromagnetic imaging using adjoint fields, Inverse
Problems, 15(1999), 1523-1558.

[14] H. Eng], M. Hanke, and A. Neubauer, Regularization of Inverse Problems, Dor-
drecht: Kluwer, 1996.

[15] M. Eller, V. Isakov, G. Nakamura, and D. Tataru, Uniqueness and stability
in the Cauchy problem for Maxwell’s and elasticity system, Nonlinear Partial

Differential Equation 14, College de France Seminar, Eds. D. Cioranescu, J. L.
Lions, Studies in Mathematics and its Applications, 31(2002), 329—350.

[16] N .E. Gibbs, W.G. Poole, Jr., and PK. Stockmeyer, An algorithm for reducing
bandwidth and profile of a sparse matrix, SIAM J. Numer. Anal, 13(1976), 236-
250.

[17] D. Gilbarg and N. Trudinger, Elliptic Partial Differential Equations of Second
Order, Springer-Verlag, New York, 1983.

[18] C. Girard and A. Dereux, Near-field optics theories, Rep. Prog. Phys, 59(1996),
657-699.

[19] T. Hohage, 0n the numerical solution of a three-dimensional inverse medium
scattering problem, Inverse Problems, 17(2001), 1743-1763

[20] H. Haddar and P. Monk, The linear sampling method for solving the electromag-
netic inverse medium problem, Inverse Problems, 18(2002), 891-906.

[21] D. Jerison and C. Kenig, Unique continuation and absence of positive eigenvalues
for Schrodinger operators, Ann. Math., 121(1985), 463-488.

[22] J. Jin, The Finite Element Methods in Electromagnetics, John Wiley & Sons,
2002.

[23] A. Kameari, Symmetric second order edge elements for triangle and tetrahedra,
IEEE "Hans. Magn., 35(1999) 1394-1397.

[24] R. G. Keys and A. B. Weglein, Generalized linear inversion and the first Born
theory for acoustic media, J. Math. Phys, 24(1983), 1444-1449.

[25] A. Kirsch and P. Monk, A finite element/spectral method for approximating the
time-harmonic Maxwell ’3 system in R3, SIAM J. Appl. Math, 55(1995), 1324-
1344.

89

 

[26] A. Kirsch and P. Monk, Corrigendum: A finite element/spectral method for
approximating the time-harmonic Maxwell’s system in R3, SIAM J. Appl. Math.,
58(1998), 2024-2028.

[27] A. Kirsch and P. Monk, A finite element method for approximating electro-
magnetic scattering from a conducting object, Numer. Math., 92(2002), 501-534.

[28] P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford University
Press, Oxford, 2003.

[29] F. Natterer, The Mathematics of Computerized Tomography, Stuttgart: Teubner,
1986.

[30] F. Natterer and F. Wiibbeling, A propagation-backpropagation method for ultra-
sound tomography, Inverse Problems, 11(1995), 1225-1232.

[31] J. C. Nédélec, Mixed finite elements in R3, Numer. l\=‘lath., 35(1980), 315-341.

[32] J. C. Nédélec, Acoustic and Electromagnetic Equations: Integral Representations
for Harmonic Problems, Springer, New York, 2000.

[33] L. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal
algorithms, Physica D, 60(1992), 259-268.

[34] M. V6geler, Reconstruction of the three-dimensional refractive index in elec-
tromagnetic scattering by using a propagation-backpropagation method, Inverse
Problems, 19(2003), 739-753.

90

 

 

MlCHlG

l

um
7
l
]

ER

lll[llll[ll/lllllllll

3 1293 02735