NEAR-FIELD IMAGING VIA INVERSE SCATTERING
By
Junshan Lin

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Applied Mathematics
2011

ABSTRACT
NEAR-FIELD IMAGING VIA INVERSE SCATTERING
By
Junshan Lin
Near-ﬁeld optics is an emerging research topic in the past few decades, mostly motivated by applications in the near-ﬁeld microscopy in an eﬀort to break the diﬀraction
limit. In the far-ﬁeld imaging, only the propagating wave components with spatial
frequency below the wavenumber are available, and it is well-known that the resolution of the image is approximately λ/2 (diﬀraction limit) [17, 51, 53]. In the near
ﬁeld, however, the bandwidth of the spatial frequency may be expanded by taking account of the evanescent (exponentially decayed) waves. Nowadays there exist various
conﬁgurations for the near-ﬁeld microscopies, see for example [51]. However, it is recognized that the images obtained from the near-ﬁeld microscopies are problematic by
visualizing the object in an analogical way [18, 47]. Therefore, the inverse scattering
theory is applied to understand how the structure of the scattering object is encoded
in the measured scattered ﬁeld. When single scattering (or Born approximation)
is assumed, the studies are complete for the near-ﬁeld scanning optical microscopy
(NSOM) and the total internal reﬂection microscopy (TIRT) within the framework
of the inverse scattering theory [19, 20, 57].
In this thesis, we focus on one speciﬁc problem where the imaging target is a
ground plane with some local disturbance. The data is collected in the near-ﬁeld
regime with a distance above the surface displacement that is smaller than the wavelength. In the recent paper [29], a linearized model has been introduced for the
nonlinear inverse scattering problem by the single scattering assumption. The authors also proposed a broadband imaging strategy for denoising and improving the
resolution of the image. In the thesis, we investigate the more general case by considering the full scattering model, for which the linearized model in [29] is no longer

valid. By the analysis of the scattered ﬁeld, it is conﬁrmed that the evanescent wave
modes which are not accessible in the far-ﬁeld regime become signiﬁcant in the near
ﬁeld. Evanescent wave modes make it possible to break the diﬀraction limit. It is
shown that such exponentially decayed modes of the scattered wave contain the high
spatial frequency information (ﬁne features) of the proﬁle. We formulate explicitly the
connection between the evanescent wave modes and the high frequency components
of the surface displacement, and present a new numerical scheme to reconstruct the
surface displacement from the boundary measurements. By extracting the information carried by the evanescent modes eﬀectively, it is shown that the resolution of the
reconstructed image is signiﬁcantly improved in the near ﬁeld. Numerical examples
show that images with a resolution of λ/10 are obtained.
To overcome the ill-posedness and the presence of local minima associated with
this nonlinear imaging problem, we propose to use multiple frequency data to image
the proﬁle of the surface displacement in the second part of the thesis. The main
idea is to march from the lowest wavenumber to the highest wavenumber. At the
ﬁxed wavenumber, by an analysis of the domain derivative for the forward scattering
map, a vector ﬁeld is chosen such that the deﬁned cost functional decreases. The
reconstructed proﬁle evolves with the chosen vector ﬁeld at the ﬁxed wavenumber
and the evolution process continues until it reaches the highest wavenumber. The
proposed reconstructed scheme is able to capture the main feature of the proﬁle at
low frequency and recover the ﬁne details at higher frequency. In particular, for a
multiple scale proﬁle, it resolves the ﬁne scale with suﬃciently high frequency.

ACKNOWLEDGMENT
I am deeply grateful to my thesis advisor Professor Gang Bao, whose continuous support and encouragements were invaluable for my four-year PhD studies at Michigan
Sate University. He introduced to me the project on near-ﬁeld imaging, and generously shared his wide expertise with me. His view on applied mathematics also helps
me to shape my own understanding in the ﬁeld, and will continue to inﬂuence me in
the long run. I am also much obliged to him for helping me to improve my English
writing skill, and for providing opportunities of attending conferences where I met
the leading ﬁgures in the ﬁeld and did presentations.
I would like to thank Professors Zhengfang Zhou, Andrew Christlieb, Chichia
Chiu, Di Liu and Jianliang Qian for serving in my thesis committee. I am so indebted
to them and Professor Habib Ammari for providing letters of reference. I am also
grateful to Professor Baisheng Yan and Barbara Miller for taking care of everything,
in particular for their help during my ﬁrst year study as a fresh PhD student.
Many thanks also go to the research group here at Michgan State University. We
have various interesting discussions and a lot of fun together. Thanks all the group
members for sharing your expertise and happiness with me: Justin Droba, Guanghui
Hu, Jun Lai, Peijun Li, Songting Luo, Yuanchang Sun, Yuliang Wang, Zhengfu Xu,
KiHyun Yun, and Hai Zhang. Thanks in particular Yuanchang, Yuliang and Maureen
for taking care of me when I ﬁrst came to US in 2007. Thanks all my friends here,
in particular Lianzhang Bao, Langhua Hu, Jie Niu, Xun Wang, Tianshuang Wu and
Qiong Zheng.
Finally, I would like to thank my parents and parents in law for their strong
support. I owe my deepest gratitude to my wife Yu Chen for all the sacriﬁce she
made. I can always feel her strong support behind me!

iv

TABLE OF CONTENTS

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1 Introduction
1.1 Near-ﬁeld
1.2 Near-ﬁeld
1.3 Near-ﬁeld

. . . . . . . . . . . . . .
Optics . . . . . . . . . . . . .
Imaging Microscopy . . . . . .
Imaging via Inverse Scattering

.
.
.
.

. .
. . .
. . .
. . .

.
.
.
.

2 Mathematical Theory of Inverse Problem . . .
2.1 Concepts of Ill-Posed Problem . . . . . . . . . . .
2.2 Singular Value Decomposition for Linear Compact
2.3 Regularization Methods for Linear Problems . .
2.4 Regularization Methods for Nonlinear Problems .

.
.
.
.

. .
. . .
. . .
. . .

.
.
.
.

.
.
.
.

. . . . .
. . . . . .
Operators
. . . . . .
. . . . . .

vii

. .
. . .
. . .
. . .

.
.
.
.

.
.
.
.

1
1
3
4

. .
. .
. .
. .
. .

.
.
.
.
.

.
.
.
.
.

7
7
8
10
13

.
.
.
.

3 Near-ﬁeld Imaging of the Surface Displacement on an Inﬁnite Ground
Plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1 Forward and Inverse Scattering Problems . . . . . . . . . . . . . . . . 15
3.2 Analysis of the Scattered Wave . . . . . . . . . . . . . . . . . . . . . 18
3.2.1 Layer Potential and Boundary Integral Equations . . . . . . . 18
3.2.2 Scattered Wave in Near-ﬁeld Regime . . . . . . . . . . . . . . 23
3.3 Near-ﬁeld Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1 Inversion Method . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.2 Fr´chet Diﬀerentiability of the Nonlinear Operator . . . . . . . 32
e
3.4 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.5 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4 Imaging with Multiple Frequency Data . . . . . . .
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Problem Formulations . . . . . . . . . . . . . . .
4.1.2 Why Multiple Frequency Data . . . . . . . . . . .
4.2 Analysis of the Forward Scattering Problem . . . . . . .
4.2.1 Notations . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Dirichlet-to-Neumann Map . . . . . . . . . . . . .
4.2.3 Well-posedness of the Forward Scattering Problem
Domain . . . . . . . . . . . . . . . . . . . . . . .
4.3 Domain Derivative of the Forward Scattering Map . . . .
4.4 Imaging with Multiple Frequency Data . . . . . . . . . .
4.4.1 Descent Direction for the Cost Functional . . . .
4.4.2 Reconstruction Scheme . . . . . . . . . . . . . . .
4.5 Numerical Examples . . . . . . . . . . . . . . . . . . . .
v

. . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
in Bounded
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .

45
45
45
47
47
47
49
51
53
59
59
61
64

4.5.1
4.5.2

Numerical approximation of the forward scattering problem .
Imaging of the local surface displacement . . . . . . . . . . . .

64
65

5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

Bibliography

. . . . . . . . . . . . . . . . . . . . . . . . . .

vi

78

LIST OF FIGURES

Images in this dissertation are presented in color

1.1

Evanescent wave. For interpretation of the references to color in this
and all other ﬁgures, the reader is referred to the electronic version of
this thesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

3.1

Setup of the problem. . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.2

eik2 (κ)d for d ∈ (0, 1.5λ) when the spatial frequency κ = 1.2k and
1.5k respectively . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.3

kc /k versus the distance varying from λ/10 to λ. . . . . . . . . . . . .

29

3.4

Real part (left) and the imaginary part (right) of the scattered ﬁeld. .

40

3.5

Near-ﬁeld (top) and far-ﬁeld (bottom) images. The solid line represents
the real image, and the dotted line is the reconstruction. . . . . . . .

41

3.6

Relative error with respect to Newton iterations. . . . . . . . . . . . .

41

3.7

Near-ﬁeld (left) and far-ﬁeld (right) images. The solid line represents
the real image, and the dotted line is the reconstruction. . . . . . . .

42

3.8

Real proﬁle, the distance between neighboring bumps is λ/10. . . . .

42

3.9

Comparison of the reconstruction (dotted line) with the real proﬁle
(solid line) for the near-ﬁeld case. . . . . . . . . . . . . . . . . . . . .

43

Problem geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

4.1

vii

4.2

Top: real proﬁle. Bottom: reconstruction (dash line) at k = 16 compared with the real curve (solid line). . . . . . . . . . . . . . . . . .

67

4.3

Relative error with respect to the wavenumber. . . . . . . . . . . . .

68

4.4

Evolution of the reconstruction at k = 1, 5, 8, 10, 12, 16 (dash line, left
to right, top to bottom). . . . . . . . . . . . . . . . . . . . . . . . . .

69

Top: relative residual with respect to the iteration number at k = kM
when only single frequency data is used. Bottom: reconstruction (dash
line) with single frequency only. . . . . . . . . . . . . . . . . . . . .

70

Reconstruction (dash line) at k = 30 compared with the real curve
(solid line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

4.7

Relative error with respect to the wavenumber. . . . . . . . . . . . .

72

4.8

Reconstruction (dash line) at k = 16 compared with the real curve
(solid line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

Reconstruction at k = 10, 20, 33, 40 (dash line, top to bottom) compared with the real curve (solid line). . . . . . . . . . . . . . . . . . .

73

4.5

4.6

4.9

viii

Chapter 1
Introduction

1.1

Near-ﬁeld Optics

Near-ﬁeld optics is the study of the evanescent wave ﬁelds and their interactions with
matter on a sub-wavelength scale. Usually the evanescent ﬁelds are localized to the
optical source region or the surface of the scattering object, and the study of the nearﬁeld optics is mostly concerned with the localized region within one wavelength. The
modern interests on the topic was mainly motivated by its applications in near-ﬁeld
imaging microscopies. We refer the reader to [51, 62] and the references therein for
detailed discussions.
The study of the near-ﬁeld optics has its origin in an eﬀort to break the diﬀraction
limit imposed by the far-ﬁeld imaging. In far-ﬁeld optics, the cut-oﬀ of the spatial
spectrum is very strict: only the propagating wave components with spatial frequency
below the wavenumber can be used. The loss of higher spatial frequencies leads to
the diﬀraction limit, which is also known as the Rayleigh resolution limit. At the
end of the nineteenth century, Abbe and Rayleigh [1, 53] derived a criterion for this
limit. The minimum distance ∆x between two point sources at which they can still
1

Figure 1.1: Evanescent wave. For interpretation of the references to color in this and
all other ﬁgures, the reader is referred to the electronic version of this thesis.

be unambiguously distinguished as two separate sources reads as

∆x = 0.61λ/N A,

where N A is the numerical aperture, and the best possible N A for optical glasses
is roughly 1.5. Therefore, the spatial resolution for a far-ﬁeld optical microscopy is
approximately λ . In near-ﬁeld optics, however, the bandwidth of spatial frequency
2
is expanded by taking account of the evanescent waves.
The evanescent waves can be described by plane waves of the form ei(k·x−ωt) ,
where at least one component of the wave vector k describing the direction of propagation is imaginary. They play a central role in near-ﬁeld optics. In the spatial
direction deﬁned by the imaginary component of k, the wave no longer propagates,
but decays exponentially. Figure 1.1 is the plot of an evanescent wave that propagates
in the xy plane, but decays exponentially in the z direction.
In fact, the total internal reﬂection at the surface of a dielectric medium generates
2

such an evanescent wave. Let us consider a plane wave impinging on a ﬂat surface
between two media characterized by the refraction index n1 and n2 respectively and
n1 > n2 . By the boundary conditions on the interface, the wave vectors k1 and k2
in the two media take the form of

k1 = (kx , ky , kz1 ),
k2 = (kx , ky , kz2 ),

n ω
|k1 | = 1 ,
c
n2 ω
.
|k2 | =
c

Here c is the speed the light in the vacuum. If the incident angle θi is the larger than
n
the critical angle θc = arcsin n2 , then kz2 becomes imaginary, and the exponential
1
n1 2 2
decay constant is |k2 | ( ) sin θi − 1.
n2
The evanescent waves are always restricted to the surfaces of imaging objects,
thus the study of the optical interactions on a subwavelength scale becomes very
signiﬁcant. This usually complicates the analysis as well as the computation — the
price we have to pay for the inclusion of evanescent waves to obtain images with
better spatial resolution.

1.2

Near-ﬁeld Imaging Microscopy

The central idea of near-ﬁeld microscopies is to retain the spatial frequencies associated with evanescent waves, thereby increasing the bandwidth of spatial frequencies.
In principle, arbitrary resolution can be achieved provided the bandwidth is inﬁnite.
However, this is at the expense of strong coupling between the source and the imaging
object. Here we brieﬂy summarize several existing near-ﬁeld imaging modalities, and
refer to the monograph [51] and [30, 25, 26, 54] for references. The list is of course
by no means to be complete.
Near-ﬁeld scanning optical microscopy (NSOM) [51, 25, 30] applies the excitation
3

beam emanating from a tiny aperture (e.g., pointed optical ﬁber) or a tiny metal
tip. The electric ﬁeld distribution changes signiﬁcantly as the aperture size becomes
smaller. In particular, when the aperture size is below the wavelength, the z component of the wave vector becomes imaginary. The strong localization of such evanescent
waves yields the subwavelength resolution that can be achieved in near-ﬁeld scanning
optical microscopy.
The scanning tunneling optical microscope (STOM) [26], also called the phonon
scanning tunneling microscope (PSTM) [54] uses far-ﬁeld illumination and near-ﬁeld
detection. To illuminate the sample, a laser beam undergoes total internal reﬂection
at the surface of the sample support. A bare tapered glass ﬁber is dipped into this
evanescent ﬁeld to locally couple some of the light into the probe where it is converted
into propagating modes that are guided towards a detector.
There are other imaging modalities such as the near-ﬁeld illumination and nearﬁeld detection conﬁguration which makes use of near-ﬁeld interactions for both excitation and detection, and the energy-transfer microscopy where the near-ﬁeld interaction between probe and sample is achieved through dipole-dipole coupling. There are
new methods being developed continuously. Apparently large diversity of methods
are usually categorized according to their speciﬁc illumination and detection conditions, and in practice, it is usually desirable to have diﬀerent specialized modalities
combined together to give more eﬃcient imaging.

1.3

Near-ﬁeld Imaging via Inverse Scattering

In an optical microscopy, the object is usually visualized in an analogical way with
little or no numerical treatment. Therefore, the images of the scattering object is usually problematic [18, 47]. There is an alternative type of near-ﬁeld imaging technique
which relies entirely on a numerical inversion procedure to reconstruct the sample
4

from the data of the scattered ﬁeld, i.e., to solve the inverse scattering problem.
In [15], M. Bertero, P. Boccacci and M. Piana discussed the concept of resolution
for the inverse diﬀraction problems that arise in the acoustic holography. They showed
that the super-resolution is possible due to the information conveyed by evanescent
waves in near-ﬁeld.
P. Carney and J. Schotland also studied the inverse scattering problem for the
near-ﬁeld microscopy and three-dimensional total internal reﬂection microscopy (TIRM).
They also showed subwavelength spatial resolution [19, 20]. Their framework on the
near-ﬁeld imaging is all based on the weak-scattering approximation (or Born approximation).
More recently, A. Sentenac et al carried out the whole nonlinear inversion process
for the total internal reﬂection microscopy [5, 6]. Unfortunately, in optics the highest
refractive index available for the prism function in TIRM is close to two, which limits
the spatial frequency of illuminating ﬁeld. In order to increase the spatial frequencies
of the illuminating ﬁeld beyond that reachable with a prism, they proposed to deposit
the sample on an optimized grating [21]. This is the so called grating-assisted optical
diﬀraction tomography.
The starting point of this thesis is also based on the idea of imaging by applying the inverse scattering theory. We study the inverse diﬀraction problem for
an unbounded obstacle which is a ground plane with some local disturbance. New
numerical schemes are developed to reconstruct the surface displacement from the
boundary measurements.
In Chapter 2, we brieﬂy summarize the mathematical theories of inverse problems,
which serves as the building blocks for the study of the inverse scattering problem.
Chapter 3 is devoted to the analysis the scattered ﬁeld in near-ﬁeld regime, and
the design a numerical method that makes use of the evanescent modes eﬀectively to
improve the resolution of the near-ﬁeld image. Most of Chapter 3 comes from [10].
5

In Chapter 4, we use multiple frequency data to image the proﬁle of the surface
displacement. A stable and accurate reconstruction method is presented and investigated with numerical simulations. In particular, for a multiple scale proﬁle, the
proposed method resolves the ﬁne scales with suﬃciently high frequency information.
Chapter 4 is mostly extracted from [11].

6

Chapter 2
Mathematical Theory of Inverse
Problem
In this chapter, we brieﬂy summarize the mathematical theories for the inverse problems, with focus on the linear compact operators in Hilbert spaces. It is intended as
an introduction to the basic ideas on the ill-posed problems and regularization methods for their solutions. We refer to the monographs [31, 43, 24, 37] and references
therein for detailed discussions.

2.1

Concepts of Ill-Posed Problem

In his lecture published in [38], Hadamard claims that a mathematical model for
a physical problem has to be properly posed or well-posed in the sense that it has
following properties:
• There exists a solution of the problem (existence).
• There exists at most one solution of the problem (uniqueness).
• The solution depends continuously on the data (stability).
7

If any of the above three criteria is violated, the problem is called ill-posed. In practice,
the instability is one of the primary interest in the study of ill-posed problems. The
violations of stability always creates serious numerical problems: an inﬁnitesimal noise
in the measurement will give rise to large errors in the solution. No mathematical trick
can make an inherently unstable problem stable. A remedy is to use the regularization
methods. But all that a regularization method can do is to recover partial information
of the solution as stable as possible.
Consider the operator equation

Kx = y,

(2.1)

where K is an operator between Hilbert spaces X and Y . The typical example of
an ill-posed problem is when K is a compact operator, since the inverse K −1 is
unbounded. We will base our exposition on (2.1) in the following sections.

2.2

Singular Value Decomposition for Linear Compact Operators

A better way to understand the structure of linear compact operator K (X → Y )
is from the spectral theory. For an self-adjoint operator K, all of its eigenvalues are
real. Moreover, K has at least one but at most a countable number of eigenvalues
with 0 as the only possible accumulation point.
Assume that the sequences {λn }∞ of the nonzero eigenvalues is ordered such
n=1
that
|λ1 | ≥ |λ2 | ≥ |λ3 | ≥ · · ·
8

Let xn be a sequence of corresponding orthonormal eigenvectors. Then for x ∈ X,
∞
λn (x, xn )xn .

Kx =
n=1

For a proof of this spectral decomposition for self adjoint operators, see for example
[45].
The spectral theorem for compact self-adjoint operators has an extension to nonselfadjoint operators. If K : X → Y is a linear compact operator, then its adjoint operator K ∗ : Y → X is also compact. We call the nonnegative square roots of the
eigenvalues for the self-adjoint compact operator K ∗ K the singular values.
Theorem 2.2.1. (Singular Value Decomposition) Let K : X → Y be a linear compact
operator, K ∗ : Y → X be its adjoint operator, and µ1 ≥ µ2 ≥ µ3 · · · > 0 be ordered
sequences of the positive singular values of K. Then there exist orthonormal systems
{xn }∞ ⊂ X and {yn }∞ ⊂ Y with the following properties:
n=1
n=1
Kxn = µn yn

and K ∗ yn = µn xn

for all n ∈ N.

The system (µn , xn , yn ) is called a singular system for K. For each x ∈ X, there
exists a singular value decompostion
∞
x = x0 +

(x, xn )xn
n=1

for some x0 ∈ Ker(K) and
∞
Kx =

µn (x, xn )yn .
n=1

The following theorem expresses the solution to the equation Kx = y in terms of
the singular system.
9

Theorem 2.2.2. (Picard) Let K : X → Y be a linear compact operator with the
singular system (µn , xn , yn ). The solution to
Kx = y
is solvable if and only if f ∈ (Ker(K ∗ ))⊥ and
∞

1
|(y, yn )|2 < ∞.
2
n=1 µn
In this case

∞
x=
n=1

1
(y, yn )xn .
µn

Picard’s theorem demonstrate the ill-posed nature of the equation Kx = y. If we
perturb the right hand side by δy, the perturbation of the solution x can be made
arbitrarily large due to the fact that the singular values tend to zero.

2.3

Regularization Methods for Linear Problems

Since the inverse of the operator K is not bounded, one remedy is to use the regularization methods. A regularization strategy is a family of linear and bounded
operators
Rα : Y → X,

α>0

such that
lim Rα Kx = x for all x ∈ X.
α→0
Let y δ be the measured data (with error) with y − y δ ≤ δ, then
xα,δ := Rα y δ
10

is an approximation of the solution of (2.1). The error x − xδ can be split into two
parts:
δ Rα

and

Rα Kx − x ,

where the ﬁrst term described the error in the data multiplied by Rα and the
second term denotes the approximation error (Rα − K −1 )y . The art of choose
the regularization parameter α will always be to ﬁnd the right compromise between
accuracy and stability.
A convenient way to construct admissible regularization strategies is given by
ﬁltering singular value systems. The idea of using ﬁlters has a long history and is
very convenient for theoretical purposes [36, 58]. Let (µn , xn , yn ) be a singular system
for K. The regularization strategies is constructed by damping the factors 1/µn .
Theorem 2.3.1. (Regularization by damping [43]) Let K : X → Y be a linear
compact operator with singular system (µn , xn , yn ).
q : (0, ∞) × (0, K ) → R

is a function with the following properties:
• |q(α, µ)| ≤ 1 for all α > 0 and 0 < µ < K .
• For each α > 0 there exists c(α) such that

|q(α, µ)| ≤ c(α)µ

and

0<µ≤ K .

• limα→0 q(α, µ) = 1 for each 0 < µ < K .
Then the operator Rα : Y → X, deﬁned by
∞
x = Rα y =
n=1

q(α, µn )
(y, yn )xn
µn

11

is a regularization strategy with Rα ≤ c(α). The function q is called a regularizing
ﬁlter for K.

There are various choices of the function q that satisﬁes the above properties.
Typically, we list three ﬁlter functions as follows:

(1) q is deﬁned by



 1, µ2 ≥ α,
q(α, µ) =
 0, µ2 < α.


This choice of q is also known as the spectral cutoﬀ.
(2) q(α, µ) = µ2 /(α + µ2 ). This choice of the ﬁlter function q is equivalent to the
Tikhonov regularization strategy [59, 60]. The Tikhonov regularization tries to
minimize the functional

Jα (x) := Kx − y + α x ,
and the minimum xα of the Tikhonov functional is the unique solution of the
normal equation
αxα + K ∗ Kxα = K ∗ y.
It is clear that the solution xα can be written in the form xα = Rα y with
Rα := (αI + K ∗ K)−1 K ∗ =

∞

µn

2
n=1 α + µn

(·, yn )xn .

(3) q(α, µ) = 1 − (1 − aµ2 )1/α . This choice of the ﬁlter function q leads to the
so-called Landweber iteration method [16, 35, 48]. The iteration scheme is the
steepest descent algorithm for the quadratic functional x → Kx − y 2 , where
12

xm is computed recursively by
xm = xm−1 − aK ∗ (Kxm−1 − y),

m = 1, 2, · · ·

The resulting regularized operator has the form
∞
1 − (1 − aµ2 )m
∗ K)j K ∗ =
Rm := a
(I − aK
(·, yn )xn .
µn
n=1
j=0
There are many other regularization strategies for the solution of the linear illposed problems, for example, the conjugate gradient method [27], the regularization
by discretization [14, 45], etc.
One important issue for the regularization strategy is the choice of the regularization parameter α such that the accuracy and stability are balanced. In practice,
usually the posteriori parameter choice rules are applied. One of the most widely used
strategies is the so-called discrepancy principle due to Morozov [49]. We also refer
to [31, 43] for the discussions of other posteriori rules to choose the regularization
parameter.

2.4

Regularization Methods for Nonlinear Problems

Let F be a nonlinear operator between Hilbert spaces X and Y , and we want to solve

F (x) = y.

(2.2)

The regularization theory for the nonlinear ill-posed problem (2.2) is by far not so
well developed as in the linear case.
13

As in the linear case, the Tikhonov regularization minimizes the functional

F (x) − y + α x .

(2.3)

The minimization problem (2.3) admits a solution, though the solution may not be
unique in general. Moreover, the Tikhonov regularization strategy is stable in the
sense of continuous dependence of the solution on the data y [31, 55]. If y in (2.3)
is replaced by the perturbed data y δ , then the corresponding solution xα,δ of (2.3)
converges to the minimum norm solution of (2.2) for general nonlinear problems when
δ → 0[55]. The rates of convergence are presented in [32, 50].
The nonlinear Landweber iteration to the nonlinear ill-posed problem updates the
solution along the steepest descent direction of the functional F (x) − y :
xm = xm−1 − aF (xm−1 )∗ (F (xm−1 ) − y),

m = 1, 2, · · ·

For nonlinear problems, the iteration will not have a global convergence property.
However, under suitable condition on the nonlinear operator F with an appropriate
stopping rule, the convergence is obtained in [39].
There are also Newton’s type methods for solving (2.2). The main idea is to
repeatedly linearize the nonlinear operator equation by solving
F (xm−1 )(x − xm−1 ) = y − F (xm−1 ).
Usually F (xm−1 ) is compact, and regularization strategy is applied when solving
the linearized problem. We refer to [31] for detailed discussions.

14

Chapter 3
Near-ﬁeld Imaging of the Surface
Displacement on an Inﬁnite
Ground Plane
3.1

Forward and Inverse Scattering Problems

Consider the scattering of the time harmonic electromagnetic wave that impinges on
an unbounded obstacle. The boundary of the obstacle is assumed to be a perturbation

Figure 3.1: Setup of the problem.
15

of the x1 x3 plane. It is also assumed that the surface displacement is invariant along
the x3 direction, and is local with respect to x1 . We further restrict the study to
the TM polarization or E-parallel case, i.e., the electric ﬁeld E = (0, 0, u(x1 , x2 )).
Consequently, the Maxwell equations are reduced to the two dimensional Helmholtz
equation.
Before introducing the forward scattering model, we describe the geometry of the
obstacle shown in Figure 3.1. Let γ ⊂ R be bounded, open, and ∂γ be the boundary
of γ. Denote the closure of γ as γ. The ground plane Γ0 := R\γ, and the local surface
displacement is represented by Γ := {x = (x1 , x2 ) | x1 ∈ γ, x2 = f (x1 )}, where the
function f is deﬁned on γ:

f (x1 ) > 0 for x1 ∈ γ,

f (x1 ) = 0 for x1 ∈ ∂γ.

By requiring f > 0 on γ, the surface displacement is directed upward. Clearly, ∂D :=
Γ∪Γ0 is the boundary of the whole unbounded obstacle on which the electromagnetic
wave impinges. The domain above ∂D is denoted as D (⊂ R2 ).
The incident wave ﬁeld ui = eikq·x is a plane wave that propagates along the
direction q = (sin θ, − cos θ)T , where θ is the incident angle and k = ω is the
c
wavenumber. Here ω is the angular frequency, and c is the speed of the wave propagating in the vacuum. Let λ = 2π denote the wavelength. If the obstacle is a ﬂat
k
perfect conductor, then the reﬂected ﬁeld ur = −eikq ·x produced by the ﬂat surface
is a plane wave propagating along the direction q = (sin θ, cos θ)T . In general, the
total ﬁeld ut from the scattering by Γ ∪ Γ0 consists of three parts: the incident wave
ui , the reﬂected wave ur , and the scattered ﬁeld. The scattered ﬁeld u satisﬁes the
Helmholtz equation
∆u + k 2 u = 0 in D.

(3.1)

Assuming that the obstacle is a perfect conductor, the total ﬁeld vanishes on the
16

boundary. Hence
u = −(ui + ur ) on ∂D.

(3.2)

It is easily seen that u = 0 on Γ0 . Moreover, the scattered ﬁeld satisﬁes the Sommerﬁeld radiation condition:

lim
r→∞

√

r

∂u
− iku
∂r

= 0,

r = |x| .

(3.3)

The forward scattering problem (4.1)-(4.3) admits a unique solution u ∈ C 2 (D) ∩
¯
C(D) if ∂D is C 2 and boundary data u|∂D is continuous [61]. Clearly here the plane
wave ui , ur ∈ C ∞ (R2 ). It is worth mentioning that there are many results on related
scattering problems in the literature. The well-posedness of the scattering problem
for an obstacle with locally downward surface displacement (f (x1 ) < 0 for x1 ∈ γ)
was studied in [3, 4]. There are also general studies on the scattering by a non-local
perturbed half plane. See for example [64] and the references therein.
In our framework for the inverse scattering problem, data is collected on the
line x2 = d above the surface displacement with a distance that is smaller than
the wavelength λ (near-ﬁeld regime). To be more precise, it is required that 0 <
d − maxx ∈¯ f (x1 ) < λ. The inverse problem is to reconstruct f from the scattered
1 γ
ﬁeld u(·, d) collected on the line x2 = d. Our work is originally motivated by the
recent paper [29], in which a linearized model has been introduced for the nonlinear
inverse scattering problem by the single scattering assumption. The authors also
proposed a broadband imaging strategy for denoising and improving the resolution
f
1 and the modulus
of the image. However, this linearized model is valid only if
λ
of its derivative f
1 simultaneously. Here we investigate the more general case
by considering the full scattering model, for which the linearized model in [29] is no
longer valid. For the reconstruction of the star-like local disturbance from the far
ﬁeld pattern of the scattered ﬁeld, we refer to [46].
17

This imaging problem shares many of the well-known diﬃculties with other inverse boundary value problems, particularly nonlinearity and ill-posedness. However,
by collecting data in the near-ﬁeld regime, the evanescent wave modes which are not
accessible in the far-ﬁeld regime (d − max f

λ) become signiﬁcant. This crucial

fact may be conﬁrmed by the analysis of the scattered ﬁeld in Section 3.2. Evanescent wave modes make it possible to break the diﬀraction limit. It is shown that
such exponentially decayed modes of the scattered wave contain exactly the high spatial frequency information (ﬁne features) of the proﬁle f . Our study is to analyze
the scattered ﬁeld carefully, and design a numerical method that makes use of the
evanescent modes eﬀectively, thus to improve the resolution of the image. Numerical
examples conﬁrm that a resolution of λ/10 is obtained in the near ﬁeld.

3.2

Analysis of the Scattered Wave

3.2.1

Layer Potential and Boundary Integral Equations

Introduce Green’s function

G(x, y) := Φ(x, y) − Φ(xr , y),
i (1)
H (k|x − y|) is the fundamental solution for the Helmholtz
4 0
equation in R2 and xr is the reﬂection of x by the x1 axis, i.e., xr = (x1 , −x2 ).
where Φ(x, y) =

Denote Γ = Γ ∪ ∂γ, Γr := { (x1 , −x2 ) | x ∈ Γ }, Γ ∪ Γr = Γ ∪ Γr ∪ ∂γ, and
Dr := { (x1 , −x2 ) | x ∈ D }. For a function ψ ∈ C(Γ), we deﬁne the single layer
potential:
u(x) =

Γ

G(x, y)ψ(y)dsy ,

x ∈ R2 \Γ ∪ Γr .

The following lemma is concerned with the limit of the normal derivative of the single
18

layer potential, when it is extended from above and below the boundary Γ. The limit
for the case when Γ is the smooth boundary of a bounded obstacle is well known
[23, 24]. Here Γ is not a closed curve. For completeness, the proof of the lemma is
provided.

Lemma 3.2.1. Assume that Γ is C 2 . For the single layer potential with continuous
density ψ, the following holds:
∂G(x, y)
∂u
(x) =
ψ(y)dsy
∂ν ±
Γ ∂νx
where ν is the unit normal directed into D,
hν(x)).

1
ψ(x),
2

x ∈ Γ,

∂u
(x) := limh→0+ ν(x) ·
∂ν ±

u(x ±

Proof. For x ∈ Γ, > 0, denote Γx, := { y ∈ Γ, |y − x| < }. We ﬁrst show that
lim lim
→0 h→0+ Γx,

∂Φ(x + hνx , y)
1
dsy = − .
∂νx
2

Let ∂Ω (the boundary of some bounded connected domain Ω) be a C 2 closed curve
such that Γ ⊂ ∂Ω. Moreover, for x ∈ Γ, the unit outward normal of x on ∂Ω coincides
with νx . Let ψ ≡ 1 on ∂Ω , from the classical results in [23, 24],
∂u
∂Φ(x, y)
1
(x) =
dsy − ,
∂ν +
2
∂Ω ∂νx

x ∈ Γ.

(3.4)

On the other hand, for any small ﬁxed number > 0,

lim
h→0+ ∂Ω

∂Φ(x + hνx , y)
dsy =
lim
∂νx
h→0+
=

+

∂Φ(x + hνx , y)
dsy
∂νx

Γx,
∂Ω\Γx,
∂Φ(x, y)
∂Φ(x, y)
dsy + lim
dsy .
∂νx
∂νx
∂Ω\Γx,
h→0+ Γx,
19

Now letting → 0, we get

lim
h→0+ ∂Ω

∂Φ(x, y)
∂Φ(x + hνx , y)
∂Φ(x, y)
dsy =
dsy + lim lim
dsy .
∂νx
∂νx
→0 h→0+ Γx,
∂Ω ∂νx
(3.5)

From (3.4), (3.5), we have

lim lim
→0 h→0+ Γx,

∂Φ(x + hνx , y)
1
dsy = − .
∂νx
2

(3.6)

. For the integral on Γ, for any ﬁxed small number > 0,

lim
h→0+ Γ

∂Φ(x + hνx , y)
ψ(y)dsy
∂νx

∂Φ(x + hνx , y)
∂Φ(x + hνx , y)
ψ(y)dsy + lim
dsy ψ(x)
∂νx
∂νx
h→0+ Γx,
∂Φ(x + hνx , y)
(ψ(y) − ψ(x))dsy
+ lim
∂νx
h→0+ Γx,
∂Φ(x, y)
∂Φ(x + hνx , y)
=
ψ(y)dsy + lim
dsy ψ(x)
∂νx
∂νx
Γ\Γx,
h→0+ Γx,
=

lim
h→0+ Γ\Γx,

∂Φ(x + hνx , y)
(ψ(y) − ψ(x))dsy .
∂νx

+ lim
h→0+ Γx,

On the other hand, there exists some constant M that
∂Φ(x + hνx , y)
dsy ≤ M
∂νx
Γx,
and
|ψ(y) − ψ(x)| → 0,

as → 0 for y ∈ Γx, .

Therefore, by letting → 0 and noting (3.6), we obtain

lim
h→0+ Γ

∂Φ(x + hνx , y)
∂Φ(x, y)
1
ψ(y)dsy =
ψ(y) − ψ(x).
∂νx
2
Γ ∂νx

20

For x ∈ Γ,

lim
h→0+ Γ

∂Φ(xr + hνx , y)
∂Φ(xr , y)
ψ(y)dsy =
ψ(y).
∂νx
∂νx
Γ

By combining the above, we have
∂u
∂G(x, y)
1
(x) =
ψ(y)dsy − ψ(x),
∂ν +
2
Γ ∂νx
The proof for

x ∈ Γ.

∂u
can be carried out in the same fashion.
∂ν −

The following representation result may serve as a starting point for the analysis
of the scattered ﬁeld.
Lemma 3.2.2. If ∂D is C 3 , then there exists ψ ∈ C(Γ) such that the solution to
(4.1)-(4.3) can be expressed as the single layer potential

u(x) =

Γ

G(x, y)ψ(y)dsy

x ∈ D.

(3.7)

Proof. Let u be the solution of (4.1)-(4.3). Note u = 0 on Γ0 , by Green’s theorem
and the radiation condition, the scattered ﬁeld takes the following form

u(x) =

∂G(x, y)
∂u(y)
u(y) − G(x, y)
dsy ,
∂νy
Γ ∂νy

x ∈ D.

(3.8)

˜
Denote the domain bounded by Γ ∪ Γr as D. Let u0 := ui + ur . Then u0 ∈ C ∞ (R2 )
˜
satisﬁes the Helmholtz equation in D. By Green’s theorem, we have

Γ

+

Γr

∂u (y)
∂Φ(x, y)
u0 (y) − Φ(x, y) 0
∂νy
∂νy
21

dsy = 0,

x ∈ D,

(3.9)

where νy is the unit normal directed to D for y ∈ Γ, νy is the unit normal directed
to Dr for y ∈ Γr .
By noting that u0 (y) = −u0 (yr ), and

∂u (yr )
∂u0 (y)
=− 0
for y ∈ Γ, we get
∂νy
∂νyr

∂Φ(x, y)
u0 (y)dsy =
Γr ∂νy
∂u (y)
Φ(x, y) 0 dsy =
∂νy
Γr

∂Φ(xr , y)
(−u0 (y))dsy ,
∂νy
Γ
∂u (y)
Φ(xr , y)(− 0 )dsy .
∂νy
Γ

(3.10)

Substituting (3.10) into (3.9) yields
∂u (y)
∂G(x, y)
u0 (y) − G(x, y) 0 dsy = 0,
∂νy
Γ ∂νy

x ∈ D.

(3.11)

A combination of (3.8) and (3.11) leads to
∂u
∂u
∂G(x, y)
(u + u0 )(y) − G(x, y)
+ 0
∂νy
∂νy
Γ ∂νy
∂u
∂u
= −
G(x, y)
+ 0 (y)dsy . x ∈ D.
∂νy
∂νy
Γ

u(x) =

(y)dsy ,

∂u
∂u
+ 0 . ψ ∈ C(Γ) follows by the standard regularization theory
∂νy
∂νy
for second-order elliptic equations and the Sobolev imbedding theorems [34]. The
Let ψ = −

proof is now complete.

Introduce the integral operator K : C(Γ) → C(Γ),

(Kψ)(x) :=

Γ

G(x, y)ψ(y)dsy ,

x ∈ Γ.

(3.12)

Denote g = −(ui + ur )| . Then the existence of the solution to the integral equation
Γ
Kψ = g follows from Lemma 3.2.2. Note that since the single layer potential can
be continuously extended to Γ, the density function ψ deﬁned in Lemma 3.2.2 is a
22

solution of the integral equation Kψ = g. Regarding the uniqueness, we have the
following result.
Proposition 3.2.3. If ∂D is C 3 , then Kψ = g is uniquely solvable if k 2 is not the
˜
˜
eigenvalue of −∆ in D for the Dirichlet problem. Here D is the domain bounded by
Γ ∪ Γr .
Proof. If Kψ = 0 for some ψ ∈ C(Γ), then the single layer potential

u(x) =

Γ

G(x, y)ψ(y)dsy

solves the exterior problem

 ∆u + k 2 u = 0


1
1


u =0
 1


√ ∂u1
 lim

− iku1
r→∞ r
∂r

in D,
on ∂D,

(3.13)

=0

and the interior problem


 ∆u + k 2 u = 0
˜
in D,
2
2
 u =0
 2
on Γ ∪ Γr .

(3.14)

˜
If k 2 is not the eigenvalue of −∆ in D for the Dirichlet problem, then (3.13) and
∂u1
∂u
(3.14) attains a unique solution respectively. Hence
− 2 = 0. By the jump
∂v
∂v
∂u1 ∂u2
condition in Lemma 3.2.1 we have −ψ =
−
= 0 on Γ. Note that ψ ∈ C(Γ),
∂v
∂v
thus ψ ≡ 0 on Γ.

3.2.2

Scattered Wave in Near-ﬁeld Regime

Now we are ready to examine the scattered ﬁeld in the near-ﬁeld regime. We point
out that the exponentially decayed (evanescent) wave modes which are localized to
23

the surface the local disturbance are signiﬁcant in the near ﬁeld, and formulate explicitly the connection between the evanescent wave modes and the high frequency
components for the proﬁle of the local disturbance.
For convenience, the Green’s function G(x, y) and the fundamental solution Φ(x, y)
are written with an explicit dependence on the ﬁrst variable:
˜
˜
˜
G(x, y) = G(x1 − y1 ; x2 , y2 ) := Φ(x1 − y1 , x2 − y2 ) − Φ(x1 − y1 , x2 + y2 ) ,
i (1)
˜
where Φ(x1 , x2 ) = H0 (k x2 + x2 ).
1
2
4
˜
For the free space fundamental solution Φ(x1 − y1 , x2 − y2 ), we have the following
plane wave decomposition:
1 i(x −y )·κ ik (κ)|x −y |
i
˜
2 2 dκ,
e 1 1 e 2
Φ(x1 − y1 , x2 − y2 ) =
4π R k2 (κ)
where




k2 (κ) =

 i


k 2 − |κ|2

|κ| < k (propagating modes),

|κ|2 − k 2 |κ| > k (evanescent modes).

(3.15)

(3.16)

The decomposition may be viewed as the sum of the plane waves that consist of
propagating and evanescent modes. It is clear that all the wave modes propagate along
the x1 direction. Along the x2 direction, when the magnitude of the spatial frequency
|κ| is below k, the wave mode also propagates; otherwise, it decays exponentially along
the x2 direction and is denoted as the evanescent mode.
˜
A similar plane wave decomposition holds for G(x1 − y1 ; x2 , y2 ) evaluated at
x2 = d:
1 i(x −y )·κ ik (κ)|d−y |
i
˜
2 − eik2 (κ)|d+y2 | dκ.
e 1 1
e 2
G(x1 − y1 ; d, y2 ) =
4π R k2 (κ)
(3.17)
Thus for the scattered ﬁeld at x = (x1 , d), by noting the single layer potential (3.7),
24

some simple calculations yield
i
1 i(x −y )·κ ik (κ)|d−f (y )|
1 − eik2 (κ)|d+f (y1 )|
e 1 1
e 2
4π R γ k2 (κ)
ψ(y1 , f (y1 )) J dy1 dκ,

u(x1 , d) =

where J =

1 + f 2 . This implies that the measured scattered ﬁeld u on the line

x2 = d can also be viewed as the superposition of the propagating and evanescent
wave modes. It is important to note that the evanescent modes with spatial frequency beyond the wavenumber k decay exponentially along the x2 direction, and
are localized to the surface of the obstacle within one wavelength. Therefore, in the
far-ﬁeld regime, such evanescent modes carried by the scattered ﬁeld are lost. However, in the near-ﬁeld regime, the evanescent modes are signiﬁcant, and the measured
scattered ﬁeld carries more information for the proﬁle of the local disturbance to be
reconstructed.
˜
On the other hand, following from the Taylor expansion of G(x1 − y1 ; d, y2 ) at
y2 = 0, the scattered ﬁeld can also be expanded as
u(x1 , d) =

˜
˜
∂ 2 G(x1 − y1 ; d, 0)
∂ G(x1 − y1 ; d, 0)
f (y1 ) +
(f (y1 ))2
2
∂y2
∂y2
γ
+O(f 3 )] ψ(y1 , f (y1 )) J dy1 .
˜
[G(x1 − y1 ; d, 0) +

where
˜
∂ 2 G(x1 − y1 ; d, 0)
˜
˜
= 0 by the symmetry property of G;
(a) G(x1 − y1 ; d, 0) =
2
∂y2
(b) for the high order term, by a direct calculation, asymptotically
f
O(f 3 ) = O ( )2
λ
We assume that

f 2
λ

∂G(x1 − y1 ; d, 0)
f.
∂y2

1, and denote ϕ = ψ
25

1 + f 2 . Then the scattered

ﬁeld is simpliﬁed as

u(x1 , d) ≈
=

˜
∂ G(x1 − y1 ; d, 0)
f (y1 )ϕ(y1 )dy1
∂y2
γ
˜
∂ G(x1 − y1 ; d, 0)
f (y1 )ϕ(y1 )dy1 .
∂y2
R

(3.18)

The equality holds since the surface displacement is supported on γ. Here we implicitly extend the deﬁnition of f and ϕ to R by setting f as 0 for x1 ∈ R \ γ. On the
other hand, from (3.17) a simple calculation yields
˜
∂ G(x1 − y1 ; d, 0)
1
ei(x1 −y1 )·κ eik2 (κ)d dκ.
=
∂y2
2π R
Therefore, by taking the Fourier transform of (3.18), we arrive at
u(κ, d) ≈ eik2 (κ)d (f ϕ)(κ),
ˆ

κ∈R,

(3.19)

where · denotes the Fourier transform.
Remark 2.1. The expression (3.19) formulates explicitly the connection between
the evanescent wave modes and the high frequency components of the proﬁle f . It
indicates that a high spatial frequency mode of the scattered ﬁeld u(κ, d) carries high
ˆ
spatial frequency information (ﬁne features) of f to be reconstructed.
f
1 and the modulus of its derivative f
λ
reduced automatically to the linear model discussed in [29].

Remark 2.2. If

1, then (3.19) is

Now we distinguish the far-ﬁeld and the near-ﬁeld cases based on (3.19). When
the spatial frequency |κ| > k, eik2 (κ)d decreases exponentially with respect to the
distance d and its value vanishes when d exceeds one wavelength (see Figure 2). Thus
in the far-ﬁeld regime (d

λ), u(κ, d) ≈ 0 for |κ| > k, i.e., the high spatial frequency
ˆ

information of f is lost in the far-ﬁeld measurement. In the context of imaging, this
26

1
κ=1.2 k
κ=1.5 k

0.8
0.6
0.4
0.2
0

0.5

d/λ

1

1.5

Figure 3.2: eik2 (κ)d for d ∈ (0, 1.5λ) when the spatial frequency κ = 1.2k and 1.5k
respectively .

implies that it is impossible to recover f with very high resolution when any noise is
present. However, in the near-ﬁeld regime with d < λ, eik2 (κ)d is not close to 0 and
the exponentially decayed modes are still signiﬁcant in the scattered ﬁeld. Therefore,
the higher spatial frequency components of f can still be retrieved by inverting the
evanescent modes of the scattered ﬁeld.

3.3
3.3.1

Near-ﬁeld Imaging
Inversion Method

Assume that the measurement u(·, d) is polluted with some additive noise n(x1 ),
which takes the following form

n(x1 ) = σ · rand(x1 ) · u(x1 , d).
27

σ is the noise level, and rand(x1 ) is a uniformly distributed random variable in [−1, 1]
for each x1 ∈ R. Moreover, rand(x1 ) is mutually independent for diﬀerent values of
x1 .
Based on the previous analysis, we present a reconstruction method. From (3.19),
we introduce the pseudo-inverse operator Id as follows:

 −ik2 (κ)d
 e
|κ| ≤ kc ,
Id (κ) =
 0

|κ| > kc ,

(3.20)

where Id (κ) is a cut-oﬀ regularized operator, and kc is a regularization parameter.
In the far-ﬁeld case, as we discuss in the previous section, only the propagating
modes can be used for imaging if noise is present, thus the cutoﬀ frequency kc = k.
In the near-ﬁeld regime, the bandwidth of the spatial frequency is expanded beyond
the wavenumber k by taking account of the evanescent waves. Note that e−ik2 (κ)d
is an exponentially increasing function with respect to |κ| when |κ| > k. Hence, the
noise may be exponentially ampliﬁed for large |κ|. For ﬁxed distance d, the cutoﬀ
frequency kc depends on the noise level (or signal-to-noise ratio). Here, following [29],
we choose kc in such a way that
2
− kc −k 2 d
eik2 (kc )d = e
= σ.

(3.21)

That is, the spatial frequency with the transfer function eik2 (κ)d below the noise
level σ is cut oﬀ. More explicitly,

k c = k 2 +

1
log σ
d


2 1/2

.

(3.22)

In view of (3.21) or (3.22), the pseudo-inverse (3.20) oﬀers a regularization strategy
28

5

kc/k

4
3
2
1

0.2

0.4

d/λ

0.6

0.8

1

Figure 3.3: kc /k versus the distance varying from λ/10 to λ.

for the inverse problem. We plot the function kc /k for various distance d in Figure 3
at 5% noise level. It is easily seen that at the ﬁxed noise level, the cutoﬀ frequency
kc

k when d < λ, i.e., the bandwidth of the spatial frequency in the near ﬁeld

is much larger than in the far ﬁeld. This guarantees better resolution for the ﬁnal
reconstruction in the near-ﬁeld regime, since the higher spatial frequency components
of f are recovered.
Denote
ˆ
h(κ) = Id (κ)ˆ(κ, d).
u

(3.23)

To compute h, the FFT may be applied to compute the inverse Fourier transform,
where h is an approximation of f ϕ. To reconstruct f from h, we need to take into
account of the boundary data on Γ, which turns out to be a (well-posed) nonlinear
problem.
We next introduce some notations for representing the surface displacement f .
Denote by C 0,1 (γ) the set of Lipschitz continuous functions on γ. Introduce the
0,1
Banach space C0 (γ) := {f | f ∈ C 0,1 (γ), f = 0 on ∂γ} with the usual norm
29

f 0,1 = f ∞ + supx ,y ∈γ,x =y
1 1
1 1

|f (x1 ) − f (y1 )|
.
|x1 − y1 |

For a ﬁxed small number δ, we deﬁne
γδ := {x1 ∈ γ | dist(x1 , ∂γ) < δ}.
For x1 ∈ γδ , let xb ∈ ∂γ such that x1 − xb = dist(x1 , ∂γ). Assume that
{ x0 , x1 , x2 · · · xN } is a set of grid points deﬁned on γ. We represent f by a piecewise
1 1 1
1
linear function, where
j−1 j
f is linear on [x1 , x1 ] for j = 1, 2, · · · , N , and is continuous on γ globally.
Moreover, it is strictly greater than 0 for the interior grid points x1 , x2 · · · xN −1 ,
1 1
1
and is 0 on the boundary x0 , xN . We denote the set of all such functions by P1 . It
1 1
˜ 0,1
is clear that P1 is a subset of C0 (γ), which is deﬁned as follows:
0,1
˜ 0,1
C0 (γ) := {f ∈ C0 (γ) | f (x1 ) > 0 for x1 ∈ γ; ∃ , δ > 0, s.t. f (x1 ) ≥
x1 − xb for x1 ∈ γδ }.
For ﬁxed grid points { x0 , x1 , x2 · · · xN }, δ = min{ x1 − x0 , xN − xN −1 },
1
1
1
1 1 1
1
1
N −1 )
1)
f (x1
f (x1
˜ 0,1
= min{
,
}. On the other hand, C0 (γ) is an open subset of
1 −x0
N −xN −1
x1 1
x1
1
0,1
C0 (γ).
˜
We rewrite the integral operator (3.12) on γ by introducing the operator Kf
deﬁned as
˜
(Kf ϕ)(x1 ) :=

γ

G(x1 − y1 ; f (x1 ), f (y1 ))ϕ(y1 )dy1 ,

x1 ∈ γ,

(3.24)

˜
1 + f 2 . The kernel of the integral operator Kf has explicit depen˜
dence on f . Here we adopt the subscript for Kf to emphasize its dependence on f .
Similarly, deﬁne gf (x1 ) := g(x1 , f (x1 )), where g = −(ui + ur )| is the boundary
˜
Γ
data.
where ϕ = ψ

30

A natural way of separating f from h computed in (3.23) is to minimize the
functional
min
f ∈P1

˜
Kf ϕ − gf
˜

L2 (γ)

+ fϕ − h

L2 (γ)

.

(3.25)

However, this minimization problem is diﬃcult to solve in practice, since f and ϕ are
both unknowns. One alternative is to solve

min f ϕ − h 2
L (γ)
f ∈P1

˜
subject to Kf ϕ = gf .
˜

It should be pointed out that in general ϕ is not Fr´chet diﬀerentiable with respect
e
˜
to f since the operator Kf is compact. Therefore, existing numerical methods such
as Newton’s method can not be applied directly.
We introduce a new function h0 such that h0 is Lipschitz continuous on γ , more¯
over



 h(x ) x ∈ γ\γ , for suﬃciently small δ;
1
1
δ
h0 (x1 ) =
 0

x1 ∈ ∂γ.

(3.26)

By choosing a small δ (usually the length of two neighboring grid points), h0 − h 2
L (γ)
is small. In practice, h0 may be chosen as the following piecewise linear function:


 h(xj ) j = 1, 2, · · · , N − 1;
j
1
h0 (x1 ) =
 0

j = 0, N.
Now we match the data on the boundary by solving the minimization problem with
a constraint:

˜
min Kf ϕ − gf 2
˜
L (γ)
f ∈P1

where ϕ ∈ C(γ) satisﬁes f ϕ = h0 .

(3.27)

Remark 3.1. The minimization problem (3.27) is a special case of (3.25). By letting
f ϕ = h0 , (3.25) becomes
31

˜
min Kf ϕ − gf 2
˜
+ h0 − h 2 ,
L (γ)
L (γ)
and h0 − h 2
is small by the deﬁnition of h0 .
L (γ)
˜ 0,1
Remark 3.2. For f ∈ P1 ⊂ C0 (γ), the function ϕ ∈ C(γ) that satisﬁes f ϕ = h0
is well and uniquely deﬁned. Moreover
1
ϕ ∞≤(
+ 1/ ) h0 0,1 .
minx ∈γ\γ f (x1 )
1
δ

f
The problem (3.27) can be solved by Newton’s method. Since ( )2
1, the
λ
iteration is expected to converge fast to the real solution, which is conﬁrmed by
our numerical examples. To linearize the problem, we require the mapping F (f ) :=
˜
Kf ϕ − gf be Fr´chet diﬀerentiable with respect to f ∈ P1 .
˜
e

3.3.2

Fr´chet Diﬀerentiability of the Nonlinear Operator
e

˜
Here and thereafter, M and M stand for some generic positive constants, whose values
may vary from step to step but should be clear from the contexts. Let L(C(γ), C(γ))
be the set of all bounded linear operators that map the functional space C(γ) to itself.
0,1
˜ 0,1
˜
Lemma 3.3.1. If f ∈ C0 (γ) ⊂ C0 (γ), then the mapping f → Kf is Fr´chet
e
0,1
diﬀerentiable from C0 (γ) to L(C(γ), C(γ)). Moreover, the Fr´chet derivative is the
e
0,1
˜
˜
linear mapping δf → (Kf ) (δf ) for δf ∈ C0 (γ), where (Kf ) (δf ) ∈ L(C(γ), C(γ))
is deﬁned as
˜
[(Kf ) (δf )]ϕ(x1 ) =

˜
∂ Φ(x1 − y1 , f (x1 ) − f (y1 ))
(δf (x1 ) − δf (y1 ))
∂x2
γ
˜
∂ Φ(x1 − y1 , f (x1 ) + f (y1 ))
(δf (x1 ) + δf (y1 )) ϕ(y1 )dy1
−
∂x2

for ϕ ∈ C(γ).
32

˜
Proof. It is clear that the mapping δf → (Kf ) (δf ) is linear.
˜
Step1: the mapping δf → (Kf ) (δf ) is bounded.
Denote

r1 =

(x1 − y1 )2 + (f (x1 ) − f (y1 ))2 , r2 =

(x1 − y1 )2 + (f (x1 ) + f (y1 ))2 .

We ﬁrst estimate the ﬁrst part:
˜
∂ Φ(x1 − y1 , f (x1 ) − f (y1 ))
∂x2

ik
4

=
≤

k
4

(1)
H0 (kr1 )
(1)
H0 (kr1 )

f (x1 ) − f (y1 )
r1
.

For a small ﬁxed constant τ0 , if |y1 − x1 | ≥ τ0 (away from the singularity), then
(1)
H0 (kr1 )

˜
≤M .

It follows that
˜
∂ Φ(x1 − y1 , f (x1 ) − f (y1 ))
(δf (x1 ) − δf (y1 ))ϕ(y1 )dy1
∂x2
γ\{|y1 −x1 |<τ0 }
≤ M δf
ϕ ∞.
0,1

On the other hand, if |y1 − x1 | < τ0 ,

(1)
H0 (kr1 )

small. We have

≤

∼ O(

1
) for τ0 suﬃciently
r1

˜
∂ Φ(x1 − y1 , f (x1 ) − f (y1 ))
(δf (x1 ) − δf (y1 ))ϕ(y1 )dy1
∂x2
{|y1 −x1 |<τ0 }
M
ϕ ∞.
ϕ ∞ ≤ M δf
|x1 − y1 | dy1 δf
0,1
0,1
{|y1 −x1 |<τ0 } r1
33

For the second part,
˜
∂ Φ(x1 − y1 , f (x1 ) + f (y1 ))
(δf (x1 ) + δf (y1 ))ϕ(y1 )dy1
∂x2
γ
˜
∂ Φ(x1 − y1 , f (x1 ) + f (y1 ))
=
[δf (y1 ) − δf (x1 )]ϕ(y1 )dy1
∂x2
γ
˜
∂ Φ(x1 − y1 , f (x1 ) + f (y1 ))
+2
ϕ(y1 )dy1 δf (x1 )
∂x2
γ
=: A1 + A2
From the estimate of the ﬁrst term, the following inequality also holds:

|A1 | ≤ M δf
ϕ ∞.
0,1
For x1 ∈ γ\γ δ ,
2
|A2 | ≤ Mδ δf
ϕ ∞,
0,1

where Mδ =

k
min
2 x ∈γ\γ
1
δ
2

(1)
H0 (kf (x1 ))

.

For x1 ∈ γ δ ,
2
|A2 | ≤

k
2

γ\[x1 −δ/2,x1 +δ/2]

+

(1)
H0 (kr2 )

[x1 −δ/2,x1 +δ/2]

dy1

δf (x1 ) ϕ ∞
(1) kδ
H0 (
2

˜
≤ M
˜
+M

˜
≤ M

x1 +δ/2
x1 −δ/2
(1) kδ
H0 ( )
2

ϕ ∞
δf
0,1
1
(y1 − x1 )2 + f (x1 )2
+ ln 2 − ln(1 −

34

dy1 δf (x1 ) ϕ ∞

δ
δ 2 + ( /2)2 x1 − xb 2

) x1 − xb 




˜
≤ M

(1) kδ
H0 ( )
2

+ ln 2 − ln(1 −

δ
δ 2 + ( /2)2 x1 − xb 2

) x1 − xb 

δf
ϕ ∞.
0,1
The last inequality follows from the fact that there exists an

such that f (x1 ) ≥

x1 − xb for x1 ∈ γδ , and δf (x1 ) ≤ δf
x − xb .
0,1 1
Therefore
|A2 | ≤ M ( , δ) δf
ϕ ∞ ,
0,1
where
˜
M ( , δ) = M
˜
−M

(1) kδ
(H0 ( )
2


˜
+ M ln 2

δ

ln(1 −
inf
x1 ∈γδ/2

δ 2 + ( /2)2 x1 − xb 2

) x1 − xb 

is a positive constant.
Therefore
˜
[(Kf ) (δf )]ϕ
≤ M ( , δ, f ) δf
ϕ ∞
∞
0,1
for any ϕ ∈ C(γ), i.e.,
˜
(Kf ) (δf )
≤ M ( , δ, f ) δf
.
0,1
L(C(γ),C(γ))

35

˜
˜
˜
Step2: An estimate of the remainder term Kf +δf − Kf − (Kf ) (δf ).

For any ϕ ∈ C(γ), by Taylor’s formula,

=

˜
˜
˜
[Kf +δf − Kf − (Kf ) (δf )]ϕ(x1 )
˜
1 ∂ 2 Φ(x1 − y1 , f (x1 ) − f (y1 ) + t(δf (x1 ) − δf (y1 )) )
∂x2
2

γ 0

dt

[δf (x1 ) − δf (y1 )]2 ϕ(y1 )dy1
˜
1 ∂ 2 Φ(x1 − y1 , f (x1 ) + f (y1 ) + t(δf (x1 ) + δf (y1 )) )
−
dt
∂x2
γ 0
2
2 ϕ(y )dy
[δf (x1 ) + δf (y1 )]
1 1
= B1 + B2 .
δf
0,1
is suﬃciently small. For a ﬁxed small constant τ0 , if |y1 − x1 | ≥ τ0 (away from
By a similar argument as in Step 1, we can estimate B1 . Assume that

singularity), then
˜
1 ∂ 2 Φ(x1 − y1 , f (x1 ) − f (y1 ) + t(δf (x1 ) − δf (y1 )) )
dt
∂x2
γ\{|y1 −x1 |<τ0 } 0
2
2
[δf (x1 ) − δf (y1 )]2 |ϕ(y1 )| dy1 ≤ M δf
ϕ ∞.
0,1

If |y1 − x1 | < τ0 ,

˜
∂ 2Φ
1
∼ O( ) for τ0 suﬃciently small. Thus
2
∂x2
r1
2

˜
1 ∂ 2 Φ(x1 − y1 , f (x1 ) − f (y1 ) + t(δf (x1 ) − δf (y1 )) )
dt
∂x2
{|y1 −x1 |<τ0 } 0
2
2
M
ϕ ∞
|x1 − y1 |2 dy1 δf
[δf (x1 ) − δf (y1 )]2 |ϕ(y1 )| dy1 ≤
2
0,1
{|y1 −x1 |<τ0 } r1
2
ϕ ∞.
≤ M δf
0,1
36

1
2
Next, the term B2 can also be split into two parts B2 and B2 :
˜
1 ∂ 2 Φ(x1 − y1 , f (x1 ) + f (y1 ) + t(δf (x1 ) + δf (y1 )) )
dt
∂x2
γ 0
2
2 ϕ(y )dy
[δf (x1 ) − δf (y1 )]
1 1

1
B2 =

and
˜
1 ∂ 2 Φ(x1 − y1 , f (x1 ) + f (y1 ) + t(δf (x1 ) + δf (y1 )) )
dt
∂x2
γ 0
2
[δf (x1 )δf (y1 )]ϕ(y1 )dy1 .

2
B2 = 4

2
It suﬃces to estimate B2 .

For x1 ∈ γ\γ δ ,
2

2
2
B2 ≤ Mδ δf
ϕ ∞,
0,1

where

Mδ =

min k 2
x1 ∈γ\γ δ
2

(1) 2k
)
H0 (
f (x1 )

+

2k
f (x1 )

(1) 2k
)
H0 (
f (x1 )

is suﬃciently small.
if δf
0,1

For x1 ∈ γ δ ,
2
2
B2

≤

γ\[x1 −δ/2,x1 +δ/2]
δf (x1 )

+

[x1 −δ/2,x1 +δ/2]

ϕ ∞
δf
0,1
37

1 ∂ 2Φ
˜
dt dy1
0 ∂x2
2

x1 +δ/2
2
1
˜
˜
≤ Mδ δf
ϕ ∞+M
dy1
0,1
x1 −δ/2 (y1 − x1 )2 + ( /2)2 x1 − xb 2
δf (x1 )

δf
ϕ ∞
0,1

˜
2
M
˜
δ (x )
≤ Mδ δf
ϕ ∞+
0,1
/2 x1 − xb f 1
2
≤ M ( , δ) δf
ϕ ∞.
0,1

δf
ϕ ∞.
0,1

Therefore
2
˜
˜
˜
Kf +δf − Kf − (Kf ) (δf )
= O( δf
),
L(C(γ),C(γ))
0,1
for suﬃciently small δf
, which completes the proof.
0,1
0,1
˜ 0,1
Let h0 ∈ C0 (¯ ) be deﬁned as in (3.26). For f ∈ P1 ⊂ C0 (γ), let ϕ ∈ C(γ)
γ
satisfy f ϕ = h0 . The following lemma concerns the Fr´chet derivative of the mapping
e
0,1
f → ϕ ( C0 (γ) → C(γ) ).
0,1
˜ 0,1
Lemma 3.3.2. If f ∈ C0 (γ) ⊂ C0 (γ), then the mapping f → ϕ deﬁned as the
0,1
above is Fr´chet diﬀerentiable from C0 (γ) to C(γ). Moreover, its Fr´chet derivative
e
e
0,1
is the linear mapping δf → ϕ for δf ∈ C0 (γ), where ϕ ∈ C(γ) and
ϕ(x1 )
ϕ (x1 ) = −
δf (x1 ) for x1 ∈ γ.
f (x1 )
Proof. It is easy to show that
1
ϕ (x1 ) ≤ (Mδ + ) ϕ ∞ δf
0,1

for x1 ∈ γ ,

1
where Mδ = minx ∈γ\γ
is a constant. Therefore
1
δ f (x1 )
ϕ

1
≤ (Mδ + ) ϕ ∞ δf
0,1
∞

0,1
and the mapping δf → ϕ is bounded from C0 (γ) → C(γ).
38

0,1
For a perturbation of f with δf ∈ C0 (γ), a perturbation of ϕ satisﬁes
(ϕ + δϕ )(f + δf ) = h0 .
If δf
is suﬃciently small, for x1 ∈ γ, the following estimate for the high order
0,1
term holds:

(ϕ + δϕ )(x1 ) − ϕ(x1 ) − ϕ (x1 )
1
=
|ϕ(x1 )| δf (x1 )
(f (x1 ) + δf (x1 ))f (x1 )
2
1
1
2
≤
(Mδ + ) ϕ ∞ δf
.
2
2
0,1
Thus
(ϕ + δϕ ) − ϕ − ϕ

2
= O( δf
),
∞
0,1

for suﬃciently small δf
.
0,1
By Taylor’s expansion, it is easily seen that the mapping f → gf (:= g(x1 , f (x1 )) )
˜
0,1
e
is also Fr´chet diﬀerentiable from C0 (γ) to C(γ). We denote its Fr´chet derivative
e
as the mapping δf → gf . Combining Lemma 3.3.1 and Lemma 3.3.2 and using the
˜
product rule, we have the following theorem:
0,1
˜ 0,1
˜
Theorem 3.3.3. If f ∈ C0 (γ) ⊂ C0 (γ), F (f ) := Kf ϕ − gf is Fr´chet dif˜
e
0,1
ferentiable from C0 (γ) to C(γ). Moreover, the Fr´chet derivative maps δf to
e
˜
˜
˜
DF (δf ) = [(K)f (δf )]ϕ + Kf (ϕ ) − gf .

3.4

Numerical Examples

First, let us consider the solution of the forward scattering problem. By Proposition
˜
3.2.3, if k 2 is not an eigenvalue of −∆ in D, to solve the forward scattering problem
(4.1)-(4.3) eﬃciently in our numerical simulation, we can ﬁrstly solve the integral
39

Figure 3.4: Real part (left) and the imaginary part (right) of the scattered ﬁeld.

equation Kψ = g and substitute ψ into (3.7). If k 2 is an eigenvalue, then the forward
scattering problem can be solved by introducing the artiﬁcial boundary (e.g. half
circle) and solving the problem in a bounded domain. Since our focus is on inverse
problem, without loss of generality, we assume that k 2 is not an eigenvalue in our
numerical examples.
In the following examples, an incident wave ui = eikq·x /1002 with normal incident direction impinges on the obstacle. The wavenumber k = 100, λ ≈ 6.28 cm,
q = (0, −1)T . In all the ﬁgures, the plots are rescaled with respect to the wavelength
λ.

Example 4.1. The real surface displacement is represented by two bumps, each
one with the size of order λ. The two bumps are close to each other, and separated
with distance λ/10. The scattered ﬁeld in the region [−3λ, 3λ] × [0, 3λ] is plotted
in Figure 3.4. Data are collected above the obstacle with distance d = λ/5 (nearﬁeld). We also assume that 5% noise is added to the simulation data. It follows from
40

y/λ

0.1
0.05
0

−1

−0.6

−0.2

0

0.2

0.6

1

0.2

0.6

1

x/λ

y/λ

0.1
0.05
0

−1

−0.6

−0.2

0

x/λ

Relative error

Figure 3.5: Near-ﬁeld (top) and far-ﬁeld (bottom) images. The solid line represents
the real image, and the dotted line is the reconstruction.

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

10

Iterations

20

Figure 3.6: Relative error with respect to Newton iterations.

41

30

y/λ

0.1
0.05
0

−1

−0.8 −0.6 −0.4 −0.2

0

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

x/λ

y/λ

0.1
0.05
0

−1

−0.8 −0.6 −0.4 −0.2

0

x/λ
Figure 3.7: Near-ﬁeld (left) and far-ﬁeld (right) images. The solid line represents the
real image, and the dotted line is the reconstruction.

Y/λ

0.2
0.1
0

−10

−8

−6

−4

−2

0

2

4

6

8

X/λ
Figure 3.8: Real proﬁle, the distance between neighboring bumps is λ/10.

42

10

Y/λ

0.2
0.1
0

−10

−8

−6

−4

−2

0

2

4

6

8

10

X/λ
Figure 3.9: Comparison of the reconstruction (dotted line) with the real proﬁle (solid
line) for the near-ﬁeld case.

(3.22) that kc ≈ 2.6k. The near-ﬁeld image fn and the real image f are plotted in
Figure 3.5 (left). Though two bumps are close to each other (λ/10), they are clear
distinguishable. Therefore, super-resolution is achieved with near-ﬁeld measurements.
To conﬁrm the convergence of Newton’s method for solving the minimization problem
(3.27), the relative error is shown with respect to the iteration number in Figure 3.6.
Here the relative error is deﬁned as
1/2
N f (xj ) − f (xj ) 2
n 1
j=0
1
N f (xj )2 1/2
j=0
1

.

The reconstruction converges fast, which leads to the real surface displacement after
the ﬁrst 20 iterations.
To compare with the far-ﬁeld image, we collect the data again with d = 5λ and 5%
noise. It is obvious that kc = k in the far-ﬁeld case. The image fn is shown in Figure
3.5 (right). It is clear that the two bumps can not be distinguished, which is due to
the fact that the high spatial frequency information of the surface displacement is lost.

Example 4.2 We consider a non-smooth proﬁle in this example. Two bumps are
separated with distance λ/10. The measurement is polluted with 5% noise. Figure
3.7 is the reconstructed near-ﬁeld and far-ﬁeld images when d = λ/5 and 5λ respec43

tively. Super-resolution is also achieved via near-ﬁeld imaging. On the other hand, it
is also observed that the accuracy of the reconstruction is deteriorated when data is
collected at far ﬁeld.

Example 4.3 The real surface displacement is a long periodic structure (see Figure
3.8). Each period is a bump with size of order λ, and two neighboring bumps are
separated with distance λ/10. The measurement distance d again is λ/5, where 5%
noise is added to the simulation data. We compare the near-ﬁeld image and the real
image in Figure 3.9. The periodic structure is also accurately reconstructed with
super-resolution.

3.5

Discussions

The proposed numerical method in the near-ﬁeld regime, which recovers high spatial
frequency modes of the surface displacement by taking account of the evanescent
modes, yields super-resolution for the reconstructed image.
Numerically, the Newton iteration for the minimization problem (3.27) converges
fast to the real solution. However, no rigorous theoretical convergence analysis is
presently available. Another issue is denoising in the near-ﬁeld regime. Noise will
be exponentially ampliﬁed in near-ﬁeld regime. The cutoﬀ wavenumber kc , which is
the critical parameter for the resolution of the reconstructed image, strongly depends
on the noise level. A denosing technique based on the broadband signal has recently
been proposed in [29] for the linearized model. For the nonlinear imaging, the problem
becomes much more challenging and is completely open.
f
To reconstruct more general local surface displacement (without assumption ( )2
λ
1), we propose to use multiple frequency data. This is discussed fully in the next
Chapter.

44

Chapter 4
Imaging with Multiple Frequency
Data
4.1
4.1.1

Introduction
Problem Formulations

We consider the very general case without speciﬁc restriction on the surface displacement in this Chapter. Figure 4.1 is the schematic setup of the scattering problem and
data collection. As in Chapter 3, we let x = (x1 , x2 ), and f (x1 ) be a function deﬁned
on the real line R with compact support γ. The curve ∂D := { x | x1 ∈ R, x2 =
f (x1 ) } represents the boundary of the whole obstacle on which the electromagnetic
wave impinges, and the domain above ∂D is denoted as D.
Let R be some positive constant such that γ ⊂⊂ (−R, R) and Γ := { x | x1 ∈
(−R, R), x2 = f (x1 ) } represent the local disturbance. Let BR be a disk with radius
+
R centered at origin. ∂BR := ∂BR ∩ D is the half circle above the ground plane,
where ∂BR is the boundary of the disk BR . Without loss of generality, we assume
+
that Γ lies below ∂BR ( see Figure 4.1). We also denote the domain bounded by
45

Figure 4.1: Problem geometry.

+
∂BR ∪ Γ as DR .
The mathematical model for the forward scattering problem is the same as in
Chapter 3. The total ﬁeld u satisﬁes the Helmholtz equation
∆u + k 2 u = 0 in D.

(4.1)

u = 0 on ∂D.

(4.2)

and vanishes on the boundary:

In addition, at inﬁnity the scattered ﬁeld us satisﬁes the Sommerfeld radiation condition:
lim
r→∞

√

r

∂us
− ikus
∂r

= 0,

r = |x| .

(4.3)

Now we are ready to present the inverse scattering problem. The total ﬁeld u is
+
collected on the half circle ∂BR , and the measurements are assumed to be available
for a set of wavenumbers { kj | j = 1, 2, · · · , M ; kj < kj+1 }. The inverse problem
is to reconstruct the local disturbance Γ from the measured multiple frequency data
+
on ∂BR .
46

4.1.2

Why Multiple Frequency Data

In the presence of noise, the resolution of the reconstructed image is usually limited
by the operating frequency of the electromagnetic wave. When the height of the
displacement is small compared to the wavelength, it is already observed in Chapter
3 that a high spatial frequency mode of the scattered wave contains exactly the
high spatial frequency information (ﬁne features) of the proﬁle. Theoretical studies
also reveal that higher frequency information may yield increased stability (less illposedness).
In fact, in our recent paper on a multiple frequency inverse source problem for
the Helmholtz equation [12], it is proved that a H¨lder type stability estimate may
o
be obtained with suﬃciently high wavenumber. In [40], Isakov also proved increased
stability for reconstructing the Schr¨dinger potential from the Dirichlet-to-Neumann
o
map with higher wavenumber. On the other hand, when only single frequency data
is available, classical regularized iterative optimization methods such as Newton’s
method ( see for example [31] ) applied to the inverse scattering problem usually fail
to compute the global minimizer. Multiple frequency data overcomes the diﬃculty of
reaching some local minimum. In [7, 8, 9, 22], a stable recursive linearization method
is developed for the inverse medium scattering problem with multiple frequency data.
The method applies the Born approximation at the lowest frequency to obtain the
initial guess of the medium, and sequentially updates at higher frequency until the
dominant modes of the medium are essentially recovered. The convergence of the
recursive linearization method has been analyzed recently in [13].

4.2
4.2.1

Analysis of the Forward Scattering Problem
Notations

47

We begin with some standard notations that will be used throughout this Chapter.
Let
H1 (DR ) := { u | u ∈ L2 (DR ), ∂j u ∈ L2 (DR ) }
be the standard Sobolev space equipped with the norm

u 1,D = u 0,D + |u|1,D ,
R
R
R
where
2
u 2
0,DR = D |u(x)| dx,
R

|u|2
1,DR =

2
j=1 DR

∂j u

2

dx.

Let (r, θ) represent the polar coordinates. Deﬁne the Sobolev space H1/2 (∂BR ) as
H1/2 (∂BR ) := { ϕ ∈ L2 (∂BR ) |

∞

1 + n2

ϕs 2 + ϕc 2
n
n

< +∞ },

n=0

where ϕs and ϕc are the Fourier coeﬃcients of the function ϕ on the circle ∂BR :
n
n
ϕs =
n
ϕc =
0

1 π
ϕ(θ) sin(nθ)dθ
π −π

π
1
ϕ(θ)dθ,
2π −π

ϕc =
n

(n ≥ 0),

1 π
ϕ(θ) cos(nθ)dθ
π −π

(n ≥ 1).

+
For a function ϕ deﬁned on the half circle ∂BR , it is extended to ∂BR by odd
reﬂection:



 ϕ(θ)
θ ∈ (0, π),
ϕ(θ) =
˜
 −ϕ(−θ) θ ∈ (−π, 0).


Deﬁne
+
+ ˜
˜
H1/2 (∂BR ) := { ϕ ∈ H1/2 (∂BR ) | ϕ ∈
H1/2 (∂BR ), where ϕ is the odd extension of ϕ }.
˜
48

Let
ϕn =

2 π
ϕ(θ) sin(nθ)dθ.
π 0

Clearly ϕs = ϕn and ϕc = 0. Hence
˜n
˜n
+
+
˜
H1/2 (∂BR ) = { ϕ ∈ H1/2 (∂BR ) |

∞

1 + n2 |ϕn |2 < +∞ },

n=1

with the norm

ϕ

4.2.2


+ =
˜
H1/2 (∂BR )

∞

1/2
1 + n2 |ϕn |2 

.

n=1

Dirichlet-to-Neumann Map

We reformulate the forward scattering model (4.1) - (4.3) in the bounded domain
+
DR by introducing the Dirichlet-to-Neumann map on ∂BR . A similar Dirichlet-toNeumann map is also introduced in [63] for the scattering from an overﬁlled cavity.
In fact, the idea of using pseudo-diﬀerential operators to reduce the inﬁnite domain
to a bounded domain has been applied to various wave simulation problems. See for
example [3, 4, 33].
Let ϕ(θ) = us | + . From (4.1) - (4.3), it follows that the scattered ﬁled us
∂BR
satisﬁes

 ∆us + k 2 us = 0

in D\BR ,




 s
+

 u =ϕ
on ∂BR ,
 us = 0
on Γ0 ,




s


 limr→∞ √r ∂u − ikus = 0, r = |x| .

∂r
The solution of the scattering problem takes the form
us (r, θ) =

∞

π
2
(1)
ϕ(θ) sin(nθ)dθ Hn (kr) sin(nθ),
(1)
n=1 πHn (kR) 0
49

(1)
where Hn is the ﬁrst kind Hankel function of order n. Let ν be the unit normal on
+
+
∂BR directed into the inﬁnite domain D. Then the normal derivative on ∂BR can
be written as
∂us
2k
(R, θ) =
∂r
π

(1)
(Hn ) (kR) π
ϕ(θ) sin(nθ)dθ sin(nθ).
(1)
Hn (kR) 0
n=1
∞

Deﬁne the Dirichlet-to-Neumann map T ( ϕ →

2k
(T ϕ)(θ) =
π

∂us
) as
∂r

(1)
(Hn ) (kR) π
ϕ(θ) sin(nθ)dθ sin(nθ).
(1)
Hn (kR) 0
n=1
∞

(4.4)

+
˜
Lemma 4.2.1. The Dirichlet-to-Neumann map T is bounded from H1/2 (∂BR ) to
+
˜
its dual space H1/2 (∂BR ) .
(1)
˜ 1/2 (∂B + ), let cn = (Hn ) (kR) . Then
Proof. For any ϕ, ψ ∈ H
R
(1)
Hn (kR)
∞
π ∞
¯ dθ = kπR
¯
¯ ds = kR
cn ϕn sin(nθ)ψ
cn ϕn ψn .
< T ϕ, ψ >=
+Tϕψ
2
0 n=1
∂BR
n=1
Here,
ϕn =

2 π
2 π
ϕ(θ) sin(nθ)dθ and ψn =
ψ(θ) sin(nθ)dθ.
π 0
π 0

From the recurrence relation of Hankel functions [2],
(1)
(1)
(1)
(1)
−Hn+1 (kR) + n Hn (kR)
−Hn+1 (kR)
(Hn ) (kR)
n
kR
cn =
=
=
+
,
(1)
(1)
(1)
kR
Hn (kR)
Hn (kR)
Hn (kR)
Thus |cn | ≤ 1 +

n
≤c
kR

1 + n2 for some positive constant c.
50

Now

|< T ϕ, ψ >| ≤

ckπR 
2

≤

∞

ckπR
2

∞
n=1

¯
1 + n2 ϕn ψn
1/2 

1 + n2 |ϕn |2 

∞



1/2
1 + n2 |ψn |2 

,

n=1

n=1

i.e.,
|< T ϕ, ψ >| ≤

ckπR
ϕ 1/2 ψ 1/2 .
˜
˜
2
H
H

The proof is now complete.

4.2.3

Well-posedness of the Forward Scattering Problem in
Bounded Domain

The normal derivative of the total ﬁeld is given by
∂u
∂us ∂ui ∂ur
∂ui ∂ur
=
+
+
= T (u) − T (ui + ur ) +
+
.
∂ν
∂ν
∂ν
∂ν
∂ν
∂ν
∂ui
∂ur
+
− T (ui + ur ), then the total ﬁeld satisﬁes the following
∂ν
∂ν
boundary value problem:
Denote g =


 ∆u + k 2 u = 0 in D ,


R


u=0
on Γ,

 ∂u


+

= T (u) + g on ∂BR .
∂ν

(4.5)

Next, we address the well-posedness of the scattering problem (4.5), which also
provides the basis for the analysis of the domain derivative of the forward scattering
map in the next section.
Deﬁne a subspace of H1 (DR )
51

+
˜1
˜
H0 (DR ) := { u ∈ H1 (DR ) | u = 0 on Γ, u| + ∈ H1/2 (∂BR ) }.
∂BR
By introducing the bilinear form

a(u, w) =

DR

u·

w − k 2 uw dx− < T u, w >
¯
¯

˜1
˜1
for u, w ∈ H0 (DR ), then u ∈ H0 (DR ) is a weak solution of the boundary value
problem (4.5) if
˜1
for all w ∈ H0 (DR ).

a(u, w) = < g, w >

(4.6)

Theorem 4.2.2. The variational problem (4.6) attains a unique solution. Moreover,
u 1,D ≤ C g
R

+
˜
H1/2 (∂BR )

for some positive constant C.

Proof. First we prove a G˚
arding type inequality. It is easy to show that
kπR
< T u, u >=
2
where

∞

cn |un |2 ,

(4.7)

n=1

(1)
(Hn ) (kR)
2 π
u(θ) sin(nθ)dθ.
cn =
, and un =
(1)
π 0
Hn (kR)

(1)
Note that Hn = Jn + i Yn , where Jn and Yn are the ﬁrst and second kind
Bessel functions respectively, and the modulus of the Hankel function is a decreasing
function [2], then

cn =
Here

2
2
1 (Jn ) (kR) + (Yn ) (kR)
Jn (kR)Jn (kR) + Yn (kR)Yn (kR)
=
< 0.
2
2
2
2
2 Jn (kR) + Yn (kR)
Jn (kR) + Yn (kR)

denotes the real part of a function.

From (4.7) and (4.8) it follows that

(< T u, u >) < 0. Consequently

2
a(u, u) ≥ u 2
˜
1,DR − c u 0,DR
52

(4.8)

for some positive constant c.
˜
To prove the existence of a weak solution, from the Fredholm alternative we only
˜1
need to prove the uniqueness. If a(u, w) = 0 for any w ∈ H0 (DR ), then the imaginary
part of a(u, u)
kπR
a(u, u) = − < T u, u >= −
2

∞

cn |un |2 = 0.

n=1

By the Wronskian formula [2],

cn =

Jn (kR)Yn (kR) − Yn (kR)Jn (kR)
2
1
=
> 0.
2 (kR) + Y 2 (kR)
2 (kR) + Y 2 (kR)
kπR Jn
Jn
n
n

+
˜1
Hence un = 0 for all n. Note that u ∈ H0 (DR ), therefore u = 0 on ∂BR , and
∂u
+
= T u = 0 on ∂BR .
T u = 0 by deﬁnition. On the other hand, g = 0 implies that
∂ν
Now u ≡ 0 in DR follows from the unique continuation result [41].
The continuous dependence of the solution on g follows from the Lax-Milgram
lemma and the Fredholm alternative.

4.3

Domain Derivative of the Forward Scattering
Map

The theory of shape sensitivity analysis has been studied extensively for various shape
optimization problems. We refer the reader to [52, 56] for detailed discussions and
references. For the inverse scattering problem, Kirsch rigorously derived the domain
derivative of the far-ﬁeld operator for a C 2 bounded obstacle in [44]. In this section,
we derive the domain derivative of the operator that maps Γ to the measurement
u| + (Theorem 4.3.1).
∂BR
53

Deﬁne the forward mapping M : Γ → u| + . It maps the local surface displace∂BR
+
ment to the measurement on the half circle ∂BR . The vector ﬁeld V (x) is deﬁned on
2
Γ. It is assumed that V (x) ∈ C0 (Γ; R2 ), i.e., twice continuously diﬀerentiable with
compact support suppV ⊂⊂ Γ. For a given vector ﬁeld V (x), denote
2
Γδ = {x + δ · V (x) | x ∈ Γ, V (x) ∈ C0 (Γ; R2 ) }
as the perturbation of Γ with V (x). Now the domain derivative of the forward mapping M with respect to the direction V (x) is deﬁned as

M (Γ) := lim
δ→0

M (Γδ ) − M (Γ)
.
δ

(4.9)

2
Theorem 4.3.1. Let u be the solution to (4.1)-(4.3). If Γ is C 2 , V ∈ C0 (Γ; R2 ),
then the domain derivative M (Γ) exists. Moreover, M (Γ) = u | + , where u
∂BR
solves

 ∆u + k 2 u = 0 in D ,


R


∂u
u = −(V · ν)
on Γ,
(4.10)

∂ν

 ∂u

+

= T (u )
on ∂BR .
∂ν
Here ν is the unit normal on Γ directed into the inﬁnite domain D.
Proof. Let δ be a small real number. The scattered ﬁeld uδ produced by the
scattering from the perturbed proﬁle Γδ satisﬁes the boundary value problem


 ∆uδ + k 2 uδ = 0 in Dδ ,


R


uδ = 0
on Γδ ,


 ∂uδ


+

= T (uδ ) + g on ∂BR ,
∂ν
+
δ
where DR is the domain bounded by Γδ and ∂BR . The weak solution of this bound54

˜1 δ
ary value problem uδ ∈ H0 (DR ) satisﬁes
aδ (uδ , wδ ) = < g, wδ >

˜1 δ
for all wδ ∈ H0 (DR ),

(4.11)

wδ − k 2 uδ wδ dx− < T uδ , wδ > .

(4.12)

where
aδ (uδ , wδ ) =

δ
DR

uδ ·

2
We extend the deﬁnition of V ∈ C0 (Γ; R2 ) to the closure of the bounded domain
DR , which is denoted by DR , such that V (x) is C 2 for x ∈ DR and V (x) = [0, 0]T
if |x| > R − α for some small positive constant α. Deﬁne the mapping Ψ(y) by letting

x = Ψ(y) = y + δV (y) for y ∈ DR .
δ
Clearly, Ψ is a C 2 mapping from DR → DR . Denote the inverse map of Ψ(y) as
δ
Φ(x), which maps DR → DR .
∂uδ
∂ uδ ∂Φm
˜
Let uδ (y) = uδ (Ψ(y)), wδ (y) = wδ (Ψ(y)). It follows that
˜
˜
= 2
m=1 ∂ym ∂x ,
∂xi
i
where Φm is the mth component of Φ. Therefore, the bilinear form (4.12) can be
written as

aδ (uδ , wδ ) =

Here, bmn =

2


DR m,n=1


δ ∂ wδ
∂u
˜
˜
bmn
− k 2 uδ wδ  J dy− < T uδ , wδ > .
˜ ˜
˜ ˜
∂ym ∂yn

¯
2 ∂Φm ∂ Φn , the Jacobian J = det DΨ.
i=1 ∂x ∂x
i
i

Deﬁne a new bilinear form aδ (˜δ , w) by letting
˜ u

aδ (˜δ , w) =
˜ u

2


DR m,n=1


δ ∂w
∂u
˜
¯
bmn
− k 2 uδ w J dy− < T uδ , w >
˜ ¯
˜
∂ym ∂yn
55

˜1
for uδ , w ∈ H0 (DR ). Then the variational problem (4.11) is equivalent to ﬁnding
˜
˜1
uδ ∈ H0 (DR ) such that
˜
aδ (˜δ , w) = < g, w >
˜ u

˜1
for all w ∈ H0 (DR ).

(4.13)

From (4.6) and (4.13), it is easily seen that uδ − u satisﬁes
˜
a(˜δ − u, w) = −(˜δ (˜δ , w) − a(˜δ , w)).
u
a u
u

(4.14)

On the right hand side,

aδ (˜δ , w) − a(˜δ , w) =
˜ u
u

2


DR m,n=1
−

uδ ·
˜

DR


δ ∂w
¯
∂u
˜
− k 2 uδ w J dy
˜ ¯
bmn
∂ym ∂yn
w − k 2 uδ w dy.
¯
˜ ¯

(4.15)

The coeﬃcient matrix (bmn ) and the Jacobian J can be further written as
J = 1+δ

· V + O(δ 2 ),

(bm,n )J = I − δ(˜mn ) + O(δ 2 ),
b

(4.16)

where I is the 2 × 2 identity matrix and the matrix
(˜mn ) =
b

V + ( V )T − (

· V )I.

(4.17)

Therefore, from (4.14) - (4.17) and the continuous dependence of uδ − u on the right
˜
hand side, it follows that
uδ − u
˜

1,DR

→ 0 as δ → 0.

56

Now

a(

uδ − u
˜
1
, w) = − (˜δ (˜δ , w) − a(˜δ , w))
a u
u
δ
δ
2
˜δ ¯
˜mn ∂ u ∂ w + k 2 (
=
b
∂ym ∂yn
DR m,n=1

Since uδ − u
˜

1,DR

· V )˜δ w dy + O(δ).
u ¯

→ 0 as δ → 0,

2
uδ − u
˜
¯
˜mn ∂u ∂ w + k 2 (
a(
, w) →
b
δ
∂ym ∂yn
DR m,n=1

· V )uw dy
¯

as δ → 0.

uδ − u
˜
exists, which we denote as u0 .
Hence the limit limδ→0
δ

Clearly, u0 is the solution to the following variational problem:
2
a(u0 , w) =

DR m,n=1

¯
˜mn ∂u ∂ w + k 2 (
b
∂ym ∂yn

· V )uw dy
¯

(4.18)

˜1
for any w ∈ H0 (DR ). By the assumption that Γ is C 2 , the scattered ﬁeld u ∈
¯
˜1
C 2 (D) ∩ C(D) [61]. Apply the formula (4.17), then for any w ∈ H0 (DR ) ∩ H2 (DR ),
2

¯
˜mn ∂u ∂ w
b
∂ym ∂yn

m,n=1
=
uT ( V + ( V )T − (
= [ (V ·

w) ·
¯

−( ( u)V ) ·
=

(V ·

w) ·
¯

· V )I) w
¯

u − ( ( w)V ) ·
¯
w] −
¯
u+

uT (

(V ·
57

u] + [ (V ·

u) ·

w
¯

· V )I w
¯

u) ·

w−
¯

· [( u ·

w)V ].
¯

+
By the Green’s formula and noting that V = 0 on ∂BR , (4.18) can be reduced to
a(u0 , w) =

DR
−

Γ

−(V ·
(V ·

w)∆u +
¯
w)
¯

(V ·

∂u
−( u·
∂ν

u) ·

w + k2(
¯

· V )uw dy
¯

w)(V · ν)dsy .
¯

On the other hand, u = w = 0 on Γ, we have (V ·

w)
¯

∂u
−( u·
∂ν

w)(V · ν) = 0 on
¯

Γ. Therefore,

a(u0 , w) =
=

where D k 2
R
a(u0 , w) =

DR

DR
DR

k 2 u(V ·
(V ·

w) +
¯
u) ·

(V ·

u) ·

w − k 2 (V ·
¯

w + k2(
¯
u)w dy +
¯

· V )uw dy
¯

DR

k2

· (uwV ) dy,
¯

· (uwV ) dy = 0 by the Divergence theorem. We ﬁnally obtain
¯
(V · u) · w − k 2 (V · u)w dy
¯
¯

˜1
for any w ∈ H0 (DR ) ∩ H2 (DR ).
(4.19)

It is easy to check that u0 deﬁned above is the weak solution of the following
boundary value problem

 ∆u + k 2 u = (∆ + k 2 )(V ·


0
0


u =0
 0
 ∂u


0 = T (u )

0
∂ν
Let u = u0 − V ·

u) in DR ,
on Γ,
+
on ∂BR ,

+
u, then u solves (4.10). Further, note that V = 0 on ∂BR ,

therefore
u = u0 = lim
δ→0

uδ − u
˜
+
= M (Γ) on ∂BR .
δ

The proof is now complete.
58

4.4

Imaging with Multiple Frequency Data

4.4.1

Descent Direction for the Cost Functional

+
For a ﬁxed wavenumber k, denote um as the measurement collected on ∂BR . For a
given curve Γ, deﬁne the cost functional

F (Γ) :=

1
M(Γ) − um 2 2
+ .
2
L (∂BR )

A descent direction on Γ is a vector ﬁeld V such that the cost functional decreases,
i.e., F (Γδ ) < F (Γ) if δ ∈ (0, t0 ] for some small positive constant t0 . We also call
this descent direction V descent velocity. The perturbation of Γ may be viewed as an
evolution process of assigning each point x on Γ a descent velocity V (x) and moving
Γ with the velocity V (x) in a small artiﬁcial time interval [0, t0 ]. The characterization
of a descent velocity for the cost functional F (Γt ) is established in the next theorem.
Theorem 4.4.1. Let Γ be C 2 , ν be the unit normal on Γ directed into the inﬁnite
domain D. u is the solution to the forward scattering problem (4.1)-(4.3), and u∗ is
the solution of the following boundary value problem:

 ∆u∗ + k 2 u∗ = 0

in DR ,



u∗ = 0
on Γ,


 ∂u∗

+

= T ∗ (u∗ ) + u − um on ∂BR .
∂ν

(4.20)

2
where T ∗ is the adjoint operator of T . If V ∈ C0 (Γ; R2 ), and satisﬁes
−

then

dF (Γt )
< 0 . Here
dt
t=0

Γ

(V · ν) ·

∂u ∂u∗
ds < 0,
∂ν ∂ν

denotes the real part of a function.
59

(4.21)

Proof. Clearly, u also solves the boundary value problem (4.5). Let u and u∗ be
the solution of (4.10) and (4.20), respectively. By the Green’s formula,
∂u∗
∂u
∂u
∂u∗
u − u∗
ds =
u − u∗
ds.
+ ∂ν
∂ν
∂ν
∂BR
Γ ∂ν

(4.22)

Using the boundary conditions in (4.10) and (4.20), the identity (4.22) is reduced to
∂u∗ ∂u
T ∗ (u∗ )u ds+
(u−um )u ds−
u∗ T (u )ds = − (V ·ν)
ds.
+
+
+
∂ν ∂ν
∂BR
∂BR
∂BR
Γ
Since T ∗ is the adjoint operator of T , we have
∂u∗ ∂u
(u − um )u ds = − (V · ν)
ds.
+
∂ν ∂ν
∂BR
Γ
Therefore by Theorem 4.3.1,
dF (Γt )
=
dt
t=0

m
+ (u − u )u ds
∂BR

=−

Γ

(V · ν) ·

∂u ∂u∗
ds.
∂ν ∂ν

The assertion of the theorem holds.

From (4.21), it is easily seen that in practice there are many possible choices for
the descent direction V deﬁned on Γ. In our case, since the local surface displacement
2
is represented by the function f , we let V = [0, v]T , where v ∈ C0 (Γ) is a compactly
supported function on Γ satisfying the inequality

−

Γ

(v · ν2 ) ·

∂u ∂u∗
ds < 0.
∂ν ∂ν

(4.23)

Here ν2 is the second component of the unit normal ν. In particular, if v · ν2 =
∂u ∂u∗
, then V is the steepest descent direction for the cost functional F (Γt ).
∂ν ∂ν
60

∂u ∂u∗
may not be a smooth and compactly supported function on
∂ν ∂ν
Γ. Therefore, we need to modify its deﬁnition to make it smooth and compactly
However,

supported. This can be accomplished by the cubic spline interpolation [28]. We
denote the new smooth, compactly supported function as vs . The choice of the
smooth version vs of v is a regularization procedure in the regularization theory for
ill-posed problems.

4.4.2

Reconstruction Scheme

Suppose that multiple frequency measurements { um | j = 1, 2, · · · , M } on the half
kj
+ are available for a set of wavenumbers { k | j = 1, 2, · · · , M }, where
circle ∂BR
j
kj > ki if j > i. Starting from the ﬂat surface as our initial guess, the proposed
reconstruction method marches from k1 to kM . At the ﬁxed wavenumber k = kj ,
by the cubic spline interpolation, a smooth version of the descent vector ﬁeld V is
deﬁned, where V = [0, v]T and v satisﬁes (4.23). The reconstructed proﬁle evolves
with the chosen descent velocity at the ﬁxed wavenumber. This evolution process
continues until k = kM .
Let Γrec be the reconstruction, Drec is the domain above the curve Γrec ∪ Γ0 .
+
rec
DR := Drec ∩ BR is the domain bounded by Γrec and ∂BR . The reconstruction
method with multiple frequency data is proposed as follows:
(1) (Initialization)
Let k = k1 , initially set Γrec = { (x1 , x2 ) | x1 ∈ (−R, R), x2 = 0 } (ﬂat
surface).
(2) (Update the reconstruction by marching along the wavenumbers)
FOR j = 1, 2, 3, · · · , M
Let k = kj ,
61

FOR n = 1, 2, 3, · · · , N (Evolution process at ﬁxed wavenumber k ):
2
Choose the descent direction Vn = [0, vn ] s.t. vn ∈ C0 (Γ) and
−

Γ

(vn · ν2 ) ·


 ∆u + k 2 u = 0 in Drec ,



R

u=0
on Γrec ,

 ∂u


+

= T (u) + g on ∂BR .
∂ν

∂u ∂u∗
ds < 0,
∂ν ∂ν

where


 ∆u∗ + k 2 u∗ = 0
rec

in DR ,



u∗ = 0
on Γrec ,


 ∂u∗

+

= T ∗ (u∗ ) + u − um on ∂BR .
k
∂ν

Update: Γrec ←− Γrec + αn Vn , αn is the step length.
END
END
Remark 5.1 The lowest wavenumber k1 has to be small to guarantee that the main
features of the proﬁle are captured. Based on the numerical investigation, one basic rule to follow is 1/k1 ≥ R, i.e., the wavelength λ1 corresponding to the lowest
wavenumber k1 is at least comparable to R. The restriction may be due to the reason
that the reconstruction with the measurement um usually captures the feature of order
k1
λ1 , which needs to be suﬃciently large in order that the main feature of the proﬁle
(of order R) is recovered in our approach. Without loss of generality, it is assumed
that R = 1 and k1 = 1. On the other hand, 1/kM has to be smaller than the scale of
the ﬁne features of the proﬁle such that the small details of the target are accurately
reconstructed as well.

Remark 5.2 Theorem 4.4.1 does not require that Γ be parameterized by some function f . In fact the assertion of the theorem is valid for any C 2 local disturbance Γ.
Therefore, our proposed method can be modiﬁed to deal with more general cases. The
corresponding descent direction may be chosen to be parallel to ν such that V = v · ν
62

∂u ∂u∗
. The evolution process of Γ with the descent velocity V can be
∂ν ∂ν
simulated by the level set method.

and v =

At low wavenumbers, the reconstruction captures the main features of the local
disturbance (low frequency modes). It also serves as the starting point for the reconstruction at the next higher wavenumber, where the evolution process continues to
recover the ﬁne features of the local disturbance (high frequency modes). Ideally, the
smaller the increment kj+1 −kj is, the better reconstruction image would be. In fact,
the increment parameter kj+1 − kj depends on the scale feature of the real proﬁle.
If kj+1 − kj is too large, the reconstruction procedure may fail to recover the Fourier
modes of the real proﬁle between kj and kj+1 . From numerical simulations discussed
in the next section, kj+1 − kj = 2 is suﬃcient to obtain accurate reconstruction.

It is also worth pointing out that the reconstruction with only the single highest
wavenumber k = kM will not yield a satisfactory image, as the convergence to the
global minimizer of the cost functional strongly relies on the initial guess. In practice,
a good initial guess is hard to choose without a priori information of the imaging
target. This may be conﬁrmed by numerical simulations presented in the next section.
However, the reconstruction with multiple frequency data overcomes the diﬃculty
of reaching some local minimum. Staring from the ﬂat surface as an initial guess,
the reconstructed image captures the large scale features of the proﬁle at a lower
wavenumber kj , which serves as the initial guess for the reconstruction at a higher
wavenumber kj+1 .
63

4.5

Numerical Examples

4.5.1

Numerical approximation of the forward scattering problem

∂u
∂u∗
and
on Γrec to deﬁne a descent direction at a ﬁxed
∂ν
∂ν
wavenumber, the solution u and u∗ in the interior domain need not be calculated.
We need to calculate only

Hence the boundary integral equation method may be applied to solve the forward
problem, which is fast and reduces the complexity of the computation by solving a
much smaller linear system compared to the ﬁnite element method. In particular, the
boundary integral equation method can handle the high wavenumber problem.
i (1)
H (k |x − y|) be the fundamental solution for the Helmholtz
4 0
equation in R2 . The solution u to (4.5) satisﬁes the following system of integral
Let Φ(x, y) =

equations:
1
∂u(y)
u(x)+
Φ(x, y)dsy −
+ (T u)(y)Φ(x, y)dsy
2
Γrec ∂νy
∂BR
∂Φ(x, y)
+
+
+ ∂νy u(y)dsy = ∂B + g(y)Φ(x, y)dsy x ∈ ∂BR ,
∂BR
R

∂u(y)
∂Φ(x, y)
Φ(x, y)dsy −
+ (T u)(y)Φ(x, y)dsy + ∂B + ∂νy u(y)dsy
Γrec ∂νy
∂BR
R
=

+
∂BR

g(y)Φ(x, y)dsy
64

x ∈ Γrec .

Similarly, we calculate

∂u∗
+
on ∂BR by solving the system
∂ν

1 ∗
∂u∗ (y)
∗ ∗
u (x)+
Φ(x, y)dsy −
+ (T u )(y)Φ(x, y)dsy
2
Γrec ∂νy
∂BR
∂Φ(x, y) ∗
m
+
+ ∂νy u (y)dsy = ∂B + (u(y) − u (y))Φ(x, y)dsy
∂BR
R

+
x ∈ ∂BR ,

∂u∗ (y)
∂Φ(x, y) ∗
Φ(x, y)dsy −
(T ∗ u∗ )(y)Φ(x, y)dsy +
u (y)dsy
+
+
Γrec ∂νy
∂BR
∂BR ∂νy
m
=
+ (u(y) − u (y))Φ(x, y)dsy x ∈ Γrec .
∂BR

4.5.2

Imaging of the local surface displacement

Several numerical examples are presented to illustrate the eﬃciency of the proposed
method. Example 1 considers a smooth upward proﬁle (f ≥ 0), the convergence
of the reconstruction is highlighted. We also compare the reconstruction with the
image when only the single highest wavenumber is used. In Example 2, we remove
the restriction that f ≥ 0. The surface displacement is represented by a general
smooth function f , where Γ+ := { (x1 , f ) | x1 ∈ (−R, R), f (x1 ) > 0 } and
Γ− := { (x1 , f ) | x1 ∈ (−R, R), f (x1 ) < 0 } are both nonempty sets. Example
3 discusses the reconstruction for the piecewise linear (nonsmooth) local surface displacement. In the last example, the reconstruction of a multiscale proﬁle is presented.
In all simulations, we use noisy data by adding 5% additive noise to the measurements.

Example 1. Smooth upward proﬁle
65



 0.2 + 0.2 cos(4πx + π) x ∈ [−0.5, 0),

1
1


f (x1 ) =
0.1 + 0.1 cos(4πx1 + π) x1 ∈ [0, 0.5],




 0
elsewhere.
The boundary of the whole obstacle is C 1 (Figure 4.2 left). An incident wave ui =
eikd·x with normal incident direction impinges on the obstacle. The wavenumbers
k1 = 1 and kM = 20. The reconstruction fn when k = 16 and the real proﬁle f are
plotted in Figure 4.2 (right). It is observed that the reconstruction is accurate even
though noise is present in measurements. To test the convergence of the proposed
method, the relative error with respect to the wavenumber is shown in Figure 4.3.
Here the relative error is deﬁned as
1/2
R
2
−R |f (x1 ) − fn (x1 )| dx1
.
R |f (x )|2 dx 1/2
1
1
−R
It is clear that the relative error decreases until the main Fourier modes of f are
recovered. Figure 4.4 illustrates the reconstructions for various wavenumbers. At low
wavenumbers, the reconstruction captures the main features of the real proﬁle, while
the ﬁne features of the proﬁle is recovered as the evolution process continues.
Next is the case when only the single highest frequency measurement is used.
M(Γrec ) − um
kM 2
+
L (∂BR )
Figure 4.5 (left) presents the relative residual
for the
m
u
kM L2 (∂B + )
R
evolution process at k = kM when only the single frequency is applied. It is clear that
the reconstruction reaches some local minimum of the cost functional. The residual
decreases and stagnates around 0.2 after 40 iterations. The corresponding image is
shown in Figure 5 (right, dash line), which deviates greatly from the real proﬁle. As
pointed out previously, in this case the convergence to the global minimizer of the
66

0.6
0.4
0.2
0
−1

−0.5

0

0.5

1

−0.5

0

0.5

1

0.6
0.4
0.2
0
−1

Figure 4.2: Top: real proﬁle. Bottom: reconstruction (dash line) at k = 16 compared
with the real curve (solid line).

cost functional strongly relies on the initial guess, which is hard to choose without
a priori information of the imaging target. Multiple frequency data overcomes the
diﬃculty, since the reconstruction at each lower frequency serves as the initial guess
for the reconstruction at the higher frequency.

Example 2. General smooth proﬁle
In this example, the surface displacement is represented by a general smooth function
f = f1 + f2 (solid line in Figure 4.6), where


 0.2(cos( 5π x ) − 1) + 0.05(1 + cos( 5π x + π)) x ∈ [− 4 , 0),

1

2 1
2 1
5

4
5π x ) − 1)
f1 (x1 ) =
x1 ∈ [0, 5 ],
0.2(cos( 2 1




 0
elsewhere
67

1
0.9

relative error

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20

wavenumber
Figure 4.3: Relative error with respect to the wavenumber.

and



 exp

f2 (x1 ) =


 0

1
5
( 4 x1 )2 − 1

4
|x1 | < ,
5
elsewhere.

Here Γ+ := { (x1 , f ) | x1 ∈ (−R, R), f (x1 ) > 0 } and Γ− := { (x1 , f ) | x1 ∈
(−R, R), f (x1 ) < 0 } are both nonempty sets. We set kM = 30, and the ﬁnal
reconstruction with 5% noise in measurements is plotted in Figure 4.6 (dash line).
Figure 4.7 shows the relative error for various wavenumbers.
The proposed reconstruction method with multiple frequency data again gives a
stable and accurate reconstruction of the general proﬁle. However, compared with
the proﬁle considered in Example 1, the convergence rate decreases when the nonnegativity or non-positivity assumption is removed from the imaging target. (see
Figure 4.3 and Figure 4.7).

68

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−1

−0.5

0

0.5

1

−1

0.6

1

−0.5

0

0.5

1

−0.5

0

0.5

1

0.2

0

0.5

0.4

0.2

0

0.6

0.4

−0.5

0

−1

−0.5

0

0.5

1

−1

0.6

0.6

0.4

0.4

0.2

0.2

0
−1

−0.5

0

0.5

1

0
−1

Figure 4.4: Evolution of the reconstruction at k = 1, 5, 8, 10, 12, 16 (dash line, left to
right, top to bottom).

69

0.8

Residual

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

10

20

30

40

50

60

70

Iterations

0.4
0.2
0
−0.2
−0.4
−1

−0.5

0

0.5

1

Figure 4.5: Top: relative residual with respect to the iteration number at k = kM
when only single frequency data is used. Bottom: reconstruction (dash line) with
single frequency only.

70

0.4
0.2
0
−0.2
−0.4
−1

−0.5

0

0.5

1

Figure 4.6: Reconstruction (dash line) at k = 30 compared with the real curve (solid
line).

Example 3. Piecewise linear proﬁle
The previous examples assume some regularity on the local surface displacement.
Next we consider a piecewise linear proﬁle (the solid line in Figure 4.8). By using
our method, the reconstruction fn captures the position and height of bumps (where
f > 0) accurately when kM = 16 (the dash line in Figure 4.8). The corner of the proﬁle is also approximated with reasonable accuracy by the C 2 smooth reconstructed
function. It will be interesting to develop a regularization strategy to approximate
the corner with some non-smooth function, which will be investigated in our future
work.

Example 4. Multiscale proﬁle
The multiscale proﬁle is represented by



f (x1 ) =

0.13 + 0.1 cos(


 0

8π
3 1
x1 + π) + 0.03 cos(16πx1 + π) x1 ∈ [− , ],
5
4 2
elsewhere.

It consists of two scales. The macroscale feature of the proﬁle is represented by the
71

1.2

relative error

1
0.8
0.6
0.4
0.2
0
0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30

wavenumber
Figure 4.7: Relative error with respect to the wavenumber.

0.6
0.4
0.2
0
−1

−0.5

0

0.5

1

Figure 4.8: Reconstruction (dash line) at k = 16 compared with the real curve (solid
line).

72

0.4
0.2
0
−1

−0.5

0

0.5

1

−0.5

0

0.5

1

−0.5

0

0.5

1

−0.5

0

0.5

1

0.4
0.2
0
−1

0.4
0.2
0
−1

0.4
0.2
0
−1

73

Figure 4.9: Reconstruction at k = 10, 20, 33, 40 (dash line, top to bottom) compared
with the real curve (solid line).

8π
x + π) ; while the microscale is represented by function
5 1
1
0.03 cos(16πx1 + π) with period . The reconstruction captures the macroscale fea8
tures when k = 10 (Figure 4.9, top left). To recover the ﬁne details of the proﬁle with
1
the period , k needs to be suﬃciently high. Here, microscale features are captured
8
when k = 40 (Figure 4.9, bottom right). The whole local disturbance is accurately
function 0.13 + 0.1 cos(

reconstructed with noisy data. On the other hand, it is observed that the resolution
of the reconstruction does not improve much from k = 10 to k = 33. This is due
to the fact that other than the macroscale feature, no scale length of the proﬁle is
comparable with the corresponding wavelength for k ∈ (10, 33).

74

Chapter 5
Conclusion
This thesis is the initial endeavor in the whole picture of the near-ﬁeld imaging. We
focus on one speciﬁc problem, where the local surface displacement on an inﬁnite
ground plane is the imaging target. The main contribution is the explicit formulation
of the connection between the evanescent wave modes and the high frequency components of the surface displacement, and a new numerical scheme to reconstruct the
surface displacement that extracts the information carried by the evanescent modes
eﬀectively. Numerical examples show that images with a resolution of λ/10 are obtained. For the general local surface displacement, a reconstruction scheme with
multiple frequency data is proposed that captures the main (large scale) features at
low frequencies and recovers the ﬁne details at high frequencies.
Numerical evidences conﬁrm that the resolution in the near-ﬁeld regime can be
signiﬁcantly improved. However, theoretical study on the uniqueness and the stability estimates remain open. In particular, the stability estimate which is able to
incorporates explicitly the dependence on the distance d would help to understand
the ill-posedness of the inverse scattering problem better in the near ﬁeld.
There are also many related imaging problems that arise in biological and nanosciences. For example, the imaging of the human cancer cell in the near-ﬁeld regime,
75

which turns out to be an inverse medium scattering problem. Though the evanescent
wave modes carried by the scattered ﬁeld is available in the near-ﬁeld measurement,
the relation between such evanescent wave modes and the high frequency components
of the medium is completely open. Such relation is the key to design eﬃcient numerical
algorithms that make full use of the evanescent waves at hand.
On the other hand, when multiple frequency data is available, one theoretically
challenging problem is to investigate the stability for the inverse scattering problem.
It is believed that a H¨lder type stability estimate can be obtained with suﬃciently
o
large band of frequencies.

76

BIBLIOGRAPHY

77

BIBLIOGRAPHY

[1] E. Abbe, Beitr¨ge zur Theorie des Mikroskops und der mikroskopischena
Wahrnehmung, Archivf. Miroskop. Anat., 9, (1873), 413.
[2] M. Abramowitz and I. Stegun (Eds.), Handbook of Mathematical Functions with
Formulas, Graphs, and Mathematical Tables, NBS Applied Mathematics Series
55, National Bureau of Standards, Washington, DC (1964).
[3] H. Ammari, G. Bao, and A.W. Wood, An integral equation method for the electromagnetic scattering from cavities, Math. Meth. Appl. Sci. 23 (2000), 10571072.
[4] H. Ammari, G. Bao, and A.W. Wood, Analysis of the electromagnetic scattering
from a cavity, Japan J. Indust. Appl. Math. 19 (2002), 301-310.
[5] K. Belkebir, P Chaumet, and A. Sentenac, Superresolution in total internal
reﬂection tomography, J. Opt. Soc. Am. A, 22 (2005), 1889-1897.
[6] K. Belkebir and A. Sentenac, High-resolution optical diﬀraction microscopy, J.
Opt. Soc. Am. A, 20 (2003), 1223-1229.
[7] G. Bao and P. Li, Inverse medium scattering for three-dimensional time harmonic Maxwell equations, Inverse Problems, 20 (2004), L1-L7.
[8] G. Bao and P. Li, Inverse medium scattering problems for electromagnetic waves,
SIAM J. Appl. Math., 65 (2005), 2049-2066.
[9] G. Bao and P. Li, Inverse medium scattering problems in near-ﬁeld optics, J.
Comput. Math., 25 (2007), 252-265.
[10] G. Bao and J. Lin, Near-ﬁeld imaging of the surface displacement on an inﬁnite
ground plane, submitted.
[11] G. Bao and J. Lin, Imaging of local surface displacement on an inﬁnite ground
plane: the multiple frequency case, submitted.
78

[12] G. Bao, J. Lin, and F. Triki, A multi-frequency inverse source problem, J. Differential Equations, 249 (2010), 3443-3465.
[13] G. Bao and F. Triki, Error estimates for the recursive linearization for solving
inverse medium problems, J. Comput. Math., 28 (2010), 725-744.
[14] J. Baumeister, Stable Soltions of Inverse Problems, Vieweg-Verlag, Braunschweig, (1987).
[15] M. Bertero, P. Boccacci and M. Piana, Resolution and super-resolution in inverse diﬀraction, Lecture Notes in Physics 486 (1997), 1-17.
[16] H. Bialy, Iterative Behandlung linearer Funktionalgleichungen, Arch. Rat. Mech.
Analy. 4 (1959), 166.
[17] M. Born and E. Wolf, Principles of Optics (6th ed.), Cambridge Univ. Press,
1980.
[18] R. Carminati and J. Greﬀet, Inﬂuence of dielectric contrast and topography on
the near ﬁeld scattered by an inhomogeneous surface, J. Opt. Soc. Am. A, 12
(1995), 2716-2725.
[19] P. Carney and J. Schotland, Inverse scattering for near-ﬁeld microscopy, Appl.
Phys. Lett. 77 (2000), 2798-800.
[20] P. Carney and J. Schotland, Three-dimensional total internal reﬂection microscopy, Opt. Lett. 26 (2001), 1072-1074.
[21] P Chaumet, K. Belkebir, and A. Sentenac, Numerical study of grating-assisted
optical diﬀraction tomography, Phy. Rev. A, 76 (2007), 013814.
[22] Y. Chen, Inverse scattering via Heisenberg’s uncertainty principle, Inverse Problems, 13 (1997), 1-13.
[23] D. Colton and R. Kress, Integral Equation Methods in Scattering Theory, Pure
and Applied Mathematics, Wiley, New York (1983).
[24] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Applied Mathematical Sciences, Vol. 93, Springer-Verlag, Berlin (1998).
[25] D. Courjon and C. Bainier, Near ﬁeld microscopy and near ﬁeld optics, Rep.
Prog. Phys. 57 (1994), 989-1028.
[26] D. Courjon, K. Sarayeddine, and M. Spajer, Scanning tunneling optical microsocpy, Opt. Commun. 71 (1989), 23-28.
[27] J. W. Daniel, The conjugate gradient emthod for linear and nonlinear operator
equaitons, SIAM J. Numer. Anal. 4 (1967), 10-26.
79

[28] C. de Boor, A Practical Guide to Splines, Applied Mathematical Sciences, vol.
27, Springer-Verlag, New York (2001).
[29] G. Derveaux, G. Papanicolaou and C. Tsogka, Resolution and denoising in
near-ﬁeld imaging, Inverse Problems 22 (2006), 1437-1456.
[30] R. Dunn, Near-Field Scanning Optical Microscopy, Chem. Rev. 99 (1999), 28912927.
[31] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, Mathematics and Its Application, Kluwer Academic Pubishers, New York
(1996).
[32] H. W. Engl, K. Kunisch, and A. Neubauer,, Convergence rates for Tikhonov
Regularization of nonliear ill-posed problems, Inverse Problems, 5 (1989) 541547.
[33] B. Engquist and A. Majda, Absorbing boundary conditions for the numerical
simulation of waves, Math. Comput., 31 (1977), 629-651.
[34] L. C. Evans, Partial Diﬀerential Equations, Graduate Studies in Mathematics,
Vol. 19, American Mathematical Society (1997).
[35] V. Fridman, A method of successive approximations for Fredholm integral equations of the ﬁrst kind, Uspki Mat. Nauk., 11 (1956), 233-234.
[36] J. Graves and P. M. Prenter, Numerical iterative ﬁlters applied to ﬁrst Freholm
integral equations, Numer. Math., 30 (1978), 281-299.
[37] C. W. Groetsch, The Theory of Tikhonov Regularization for Fredholm Equations of the First Kind. Pitman, Boston (1984).
[38] J. Hadamard, Lectures on Cauchy’s problem in Linear Partial Diﬀerential Equations, Yale University Press, New Haven (1923).
[39] M. Hanke, A. Neubauer, and O. Scherzer A convergence analysis of the Landweber iteration for nonlinear ill-posed problems, Numer. Math., 72 (1995), 21-37.
[40] V. Isakov, Increasing stability for the Schr¨dinger potential from the Dirichleto
to-Neumann map, Discrete Contin. Dynam. Systems-S, to appear.
[41] D. Jerison and C. Kenig, Unique continuation and absence of positive eigenvalues for Schr¨dinger operators. Ann. Math., 121(1985), 463-488.
o
[42] J. Jin, The Finite Element Method in Electromagnetics (2nd ed.), Jonh Wiley
& Sons (2002).
[43] A. Kirsch, An Introduction to the Mathematical Theory of Inverse Problems,
Applied Mathematical Sciences, Vol. 120, Springer-Verlag, New York, 1996.
80

[44] A. Kirsch, The domain derivative and two applications in inverse scattering
theory, Inverse Problems, 16 (1993), 81-96.
[45] R. Kress, Linear Integral Equations, Applied Mathematical Sciences, Vol. 82,
Sringer-verlag, New York (1999).
[46] R. Kress and T. Tran, Inverse scattering for a locally perturbed half-plane, Inverse Problems 16 (2000), 1541-1559.
[47] D. Labeke and D. Barchiesi, Scanning-tunneling optical microscopy: a theoretical macroscopic approach, J. Opt. Soc. Am. A 9 (1992), 732-739.
[48] L. Landweber, An iteration formula for Freholm integral equations of the ﬁrst
kind, Am. J. Math. 73 (1951), 615-624.
[49] V. A. Morozov, On the soltuion of functional equations by the method of regularization, Soviet Math. Dokl 7 (1966), 414-417.
[50] A. Neubauer, Tikhonov Regularization for nonliear ill-posed problems: optimal
convergence and ﬁnite-dimensional approximation, Inverse Problems, 5 (1989)
541-547.
[51] L. Novotny and B. Hecht, Principles of Nano-Optics, Cambridge University
Press (2006).
[52] O. Pironneau, Optimal Shape Design for Elliptic Systems, Springer-Verlag, New
York (1984).
[53] L. Rayleigh, On the theory of optical images with special reference to the optical
microscope, Phil. Mag. 5 (1896), 167-195.
[54] R. C. Reddick, R. J Warmack, D. W. Chilcott, S. L. Sharp , and T. L. Ferrel,
Photon scanning tunneling microscopy, Rev. Sci. Instrum. 61 (1990) 3669-3677.
[55] T. I. Seidman and C. R. Vogel, Well-posedness and convergence of some regularization methods for nonlinear ill-posed problems, Inverse Problems, 5 (1989)
227-238.
[56] J. Sokolowski and J. Zolesio, Introduction to Shape Optimization: Shape Sensitivity Analysis, Springer-Verlag, Berlin (1992).
[57] J. Sun, P. Carney, and J. Schotland, Near-ﬁeld scanning optical tomography: a
nondestructive method for three-dimensional nanoscale imaging, IEEE J. Sel.
Top. Quant. 12 (2006), 1072-1082.
[58] S. Twomey, The application of numerical ﬁltering to the solution of integral
equations encountered in indirect sensing measurements, J. Franklin Inst. 279
(1965), 95-109.
81

[59] A. V. Tikhonov, On the solution of incorrectly formulated problems and the
regularization method, Soviet Math. Doklady, 4 (1963), 1035-1038.
[60] A. V. Tikhonov, Regularization of incorectly posed problems, Soviet Math. Doklady, 4 (1963), 1624-1627.
[61] A. Willers, The Helmholtz equation in disturbed half-spaces, Math. Meth. Appl.
Sci. 9 (1987), 312-323.
[62] E. Wolf (Ed.) Progress in Optics 50, Elsevier, Amsterdam, The Netherlands
(2007).
[63] A. Wood, Analysis of electromagnetic scattering from an overﬁlled cavity in the
ground plane, J. Comput. Phys., 215 (2006), 630-641.
[64] B. Zhang and S. N. Chandler-Wilde, Integral equation methods for scattering by
inﬁnite rough surfaces, Math. Meth. Appl. Sci. 26 (2003), 463-488.

82