Inverse wave scattering problems: fast algorithms, resonance and - - PDF document

Inverse wave scattering problems: fast algorithms, resonance and applications

Wagner B. Muniz Department of Mathematics Federal University of Santa Catarina (UFSC) w.b.muniz@ufsc.br III Col´

• quio de Matem´

atica da Regi˜ ao Sul 2014

Inverse scattering (acoustics, EM)

i

u

s

u

D

ui(x) = known incident wave us(x) = measured scattered wave

incident ui + scattered us = total field u Time-harmonic assumption: ω = frequency acoustics: p(x, t) = ℜe

{

u(x)e−iωt} , EM: (E, H)(x, t) = ℜe

{

(E, H)(x)e−iωt}

1

Inverse scattering (acoustics, EM)

i

u

s

u

D

ui(x) = known incident wave us(x) = measured scattered wave

Direct problem: Given D (and its physical properties) describe the scattered field us Inverse ill-posed problem : Determine the support (shape) of D from the knowledge of us far away from the scatterer (far field region)

2

Outline

• 1. Approaches for inverse scattering:

− Traditional methods − Qualitative sampling methods

• 2. Forward scattering

− Radiating (outgoing) solutions − Rellich’s lemma

• 3. Elements of inverse scattering theory

− Far field operator − Herglotz wave function

• 4. Sampling formulation

− Fundamental solution − Linear sampling method − Factorization method

• 5. Resonant frequencies

− Modified Jones/Ursell far-field operator − Object classification algorithm

• 6. Applications

− Real experimental data − Buried obstacles detection

3

• 1. Approaches for inverse scattering

Qualitative/sampling schemes Goal: try to

• recover shape as opposed to physical properties
• recover shape and possibly some extra info

Fixed frequency of incidence ω:

i

u

s

u

D

Sampling: Collect the far field data u∞ (or the near

field data us) and solve an ill-posed linear integral equation for each sample point z

4

Inverse Scattering Methods

Nonlinear optimization methods

Kleinmann, Angell, Kress, Rundell, Hettlich, Dorn, Weng Chew, Ho- hage, Lesselier ...

• need some a priori information

− parametrization, # scatterers, etc

• flexibility w.r.t. data
• need forward solver (major concern)
• full wave model
• inverse crimes not uncommon!

Asymptotic approximations (Born, iterated- Born, geometrical optics, time-reversal/mi- gration, ...) Bret Borden, Cheney, Papanicolaou, ...

• need a priori information so linearizations

be applicable (not for resonance region)

• (mostly) linear inversion schemes
• radar imaging with incorrect model?

Qualitative methods (sampling, Factoriza- tion, Point-source, Ikehata’s, MUSIC?...)

Colton, Monk, Kirsch, Hanke-Bourgeois, Cakoni, Pot- thast, Devaney, Hanke, Ikehata, Ammari, Haddar, ...

• no forward solver
• no a priori info on the scatterer
• no linearization/asymptotic approx.:

– full nonlinear multiple scattering model

• need more data
• do not determine EM properties (σ, ϵr)

5

• 2. Forward wave propagation 101

Wave equation (pressure p = p(x, t), velocity c) ∂2 ∂t2 p − c2△p = 0 Time-harmonic dependency: ω = frequency p(x, t) = ℜe

{

u(x)e−iωt} Helmholtz (reduced wave) equation: (−i ω)2u − c2 △u = 0 ⇒ −△u − k2u = 0 where k = ω/c is the wavenumber. Plane wave incidence ’Plane wave’ in the direction d, |d| = 1, p(x, t) = cos {k(x · d − c0t)} = ℜe

{

eikx·de−iωt} Plane wave ui(x) = eikx·d satisfies −△ui − k2ui = 0 em R3, where k = ω/c0

6

Forward scattering

Incident field (say plane wave or point source) −△ui − k2ui = f in R3, where k = ω/c0 Helmholtz equation for the total field −△u − k2u = 0 in R3 \ D, Bu = 0 on ∂D, Total field u = ui + us, us perturbation due to D Boundary condition (impenetrable) Bu := ∂νu + iλu impedance (Neumann λ = 0) = u Dirichlet/PEC Analogous to Maxwell with ∇ × ∇ × E − k2E = F in R3 \ D

7

Sommerfeld/Silver-M¨ uller conditions

Exterior boundary value problem for us Uniqueness: us travels away from the obstacle −△us − k2us = 0 in R3 \ D, Bus = f := −Bui on ∂D, lim

R→∞

∫

r:=|x|=R

• ∂

∂rus − ikus

• 2

ds(x) = 0 (Sommerfeld radiation condition) Here x = |x|ˆ x = rˆ x, ˆ x ∈ Ω Notation: Ω unit sphere Sommerfeld: ”... energy does not propagate from infinity into the domain ...”

8

Radiating solutions II

Sommerfeld radiation condition on us −△us − k2us = 0 in R3 \ D, Bus = f := −Bui on ∂D, lim

R→∞

∫

r:=|x|=R

• ∂

∂rus − ikus

• 2

ds(x) = 0 Asymptotic behavior of radiating solutions

• Def. us is radiating if it satisifies

– Helmholtz outside some ball and – Sommerfeld radiation condition Theor. If us is radiating then us(x) = eik|x| |x| u∞(ˆ x) + O

(

1 |x|2

)

0.5 1 1.5 30 210 60 240 90 270 120 300 150 330 180

9

Rellich’s lemma [1943]

Key tool in scattering theory: Identical far field patterns ⇓ Identical scattered fields (in the domain of definition) Rellich’s lemma (fixed wave number k > 0) If v1

∞(ˆ

x) = v2

∞(ˆ

x) for infinitely many ˆ x ∈ Ω then vs

1(x) = vs 2(x), x ∈ R3 \ D.

That is, if v1

∞(ˆ

x) = 0 for ˆ x ∈ Ω then vs

1(x) = 0, x ∈ R3 \ D.

Remark: R >> 1,

∫

|x|=R |vs(x)|2ds(x) ≈

∫

Ω |v∞(ˆ

x)|2ds(ˆ x)

10

• 3. Inverse Scattering Theory

Inverse problem: ill-posed and nonlinear Given several incident plane waves with dir. d ui(x, d) = eikx·d, measure the corresponding far-field pattern u∞(ˆ x, d), ˆ x ∈ Ω and determine the support of D

Re 100 200 300 50 100 150 200 250 300 350 Im 100 200 300 50 100 150 200 250 300 350

11

Far field operator (data operator): F : L2(Ω) → L2(Ω) (Fg)(ˆ x) :=

∫

Ω u∞(ˆ

x, d)g(d)ds(d) Remark 1: F is compact (smooth kernel u∞) Remark 2: F is injective and has dense range whenever k2 ̸= interior eigenvalue Proof: Fg = 0 implies (Rellich)

∫

Ω us(x, d)g(d)ds(d) = 0, x ∈ R3 \ D

−B

∫

Ω ui(x, d)g(d)ds(d) = 0, x ∈ ∂D

that is, − Bvg(x) = 0, x ∈ ∂D where Herglotz wave function: vg(x) :=

∫

Ω eikx·dg(d)ds(d), kernel g ∈ L2(Ω)

so that vg satisfies the interior e-value problem −△vg − k2vg = 0 in D, Bvg = 0 on ∂D and vg = 0, g = 0, if k2 ̸= eigenvalue

• 12

Far field operator (data operator): ( ↗ ) F : L2(Ω) → L2(Ω) (Fg)(ˆ x) :=

∫

Ω u∞(ˆ

x, d)g(d)ds(d) Obs.: F normal in the Dirichlet, Neumann and non-absorbing medium cases

13

Herglotz wave function

Superposition with kernel g

∫

Ω

eikx·dg(d)ds(d) ❀

∫

Ω

us(x, d)g(d)ds(d) ❀

∫

Ω

u∞(ˆ x, d)g(d)ds(d) ∥ ∥ ∥ vg(x) ❀ vs(x) ❀ (Fg)(ˆ x)

By superposition the incident Herglotz func- tion vg(x) induces the far field pattern (Fg)(ˆ x) The fundamental solution (R3): Φ(x, z) := eik|x−z| 4π|x − z|, x ̸= z, is radiating in R3 \ {z}. Fixing the source z ∈ R3 as a parameter, then Φ(·, z) has far field pattern Φ(x, z) := eik|x| |x| Φ∞(ˆ x, z) + O

(

1 |x|2

)

, withΦ∞(ˆ x, z) = 1 4πe−ikˆ

x·z

14

• 4. Linear Sampling Method (LSM)

Far field equation Let z ∈ R3. Consider Fgz(ˆ x) = Φ∞(ˆ x, z) It is solvable only in special cases, if z = z0 and D is a ball centered at z0. In general a solution doesn’t exist.

• Ex. 2D Neumann obstacle: (k = 3.4, k = 4)

k =3.4 −2 2 −3 −2 −1 1 2 3 10 20 30 40 50 60 k =4 −2 2 −3 −2 −1 1 2 3 10 20 30 40 50 60

z inside D, ||gz|| remains bounded z outside D, ||gz|| becomes unbounded

Nevertheless the regularized algorithm is nu- merically robust and the following approxima- tion theorem holds

15

LSM theorem

( ↗ ) Theorem If −k2 ̸= Dirichlet eigenvalue for the Laplacian in D then

(1) For any ϵ > 0 and z ∈ D, there exists a gz ∈ L2(Ω) such that

• ∥Fgz − Φ∞(·, z)∥L2(Ω) < ϵ,

and

• limz→∂D ∥gz∥L2(Ω) = ∞,

limz→∂D ∥vgz∥H1(D) = ∞. (2) For any ϵ > 0, δ > 0 and z ∈ R3 \ D, there exists a gz ∈ L2(Ω) such that

• ∥Fgz − Φ∞(·, z)∥L2(Ω) < ϵ + δ

and

• limδ→0 ∥gz∥L2(Ω) = ∞, limδ→0 ∥vgz∥H1(D) = ∞

where vgz is the Herglotz function with kernel gz.

16

LSM motivation (Dirichlet)

• Assume u∞(ˆ

x, d) known for ˆ x, d ∈ Ω corresponding to ui(x, d) = eikx·d

• Let z ∈ D and g = gz ∈ L2(Ω) solve Fg = Φ∞(·, z):

∫

Ω u∞(ˆ

x, d)g(d)ds(d) = Φ∞(ˆ x, z)

• Rellich’s lemma:

∫

Ω us(x, d)g(d)ds(d) = Φ(x, z),

x ∈ R3 \ D

• Boundary condition us(x, d) = −eikx·d on ∂D implies:

−

∫

Ω eikx·dg(d)ds(d) = Φ(x, z),

x ∈ ∂D, z ∈ D. If z ∈ D and z → x ∈ ∂D then ||g||L2(Ω) → ∞ since |Φ(x, z)| → ∞ Same analogy: Neumann, impedance, inho- mogeneous medium

17

Factorization method (Dirichlet)

Generalized scattering problem: f ∈ H1/2(∂D) ∆v + k2v = 0 in R3 \ D, v = f on ∂D, v radiating Data to far-field operator: takes f into v∞ G : H1/2(∂D) → L2(Ω), f ❀ Gf := v∞ Theorem z ∈ D iff Φ∞(·, z) ∈ Range(G) Proof: Rellich + singularity of Φ(·, z) at z.

18

Factorization: characterizes range of G (and therefore D by the previous theorem) in terms

• f the data operator F, i.e.

in terms of the singular system of F Theorem Let k2 ̸=Dirichlet e-value of −∆ in

• D. Let {σj, ψj, ϕj} be the singular system of F.

Then z ∈ D iff

∞

∑

1

|(Φ∞(·, z), ψj)|2 σj < ∞

( ↗ )

Factorization method (Dirichlet) II

Factorization of the far field operator: F = −GS∗G∗ where S is the adjoint of the single layer po- tential

• Obs. This corresponds to solving in L2(Ω)

(F ∗F)1/4g = Φ∞(·, z) i.e. Range(G) = Range(F ∗F)1/4

• 5. And resonant frequencies?

2 Dirichlet eigenvalues (peanut) Lack of injectivity of F

k =1.6805 k =2.6 k =3.0418

• Is it a true failure?
• Can we get some extra info about the

scatterer at eigenfrequencies?

• First an algorithm that works for all k.

19

Modified far field operator

( ↗ )

Back to Jones, Ursell (1960s), Kleinman & Roach and Colton & Monk (1988, 1993)

Find a ball BR(0) of radius R > 0, BR ⊂ D. Define amn, n = 0, 1, ..., |m| ≤ n, such that (1) |1+2amn| > 1 for all n = 0, 1, . . . , , |m| ≤ n (2)

∞

∑

n=0 n

∑

m=−n

( 2n

keR

)2n

|amn| < ∞,

R

O

D

Define a series of far field patterns u0

∞(ˆ

x, d) := 4π ik

∞

∑

n=0 n

∑

m=−n

amnY m

n (ˆ

x)Y m

n (d),

where Y m

n

= spherical harmonics

20

Modified far field operator ( ↗ ) (F0 g)(ˆ x) :=

∫

Ω

(

u∞(ˆ x, d) − u0

∞(ˆ

x, d)

)

g(d)ds(d) Each term of the series of far field patterns 4π ik amnY m

n (d)Y m n (ˆ

x) corresponds to radiating Helmhotz solutions of the form us,0

mn(x) = 4πinamnY m n (d) h(1) n

(k|x|)Y m

n (ˆ

x)

21

Modified LSM valid for all k > 0

Theor. F0 : L2(Ω) → L2(Ω) is injective with dense range. Theor.

(as before with F0, without restriction on k)

Jones/Ursell modification F0:

k =1.6805 k =2.6 k =2.8971

Before:

k =1.6805 k =2.6 k =3.0418

22

Object classification at e-frequencies

Claim: at eigenfrequencies, imaging ||gz|| in- dicates the zeros of the corresponding eigen- functions (easy to see in the 2D/3D spherical case) Corollary: Given the far field data for k ∈ [k0, k1] (containing e-freq.) then one can classify a scatterer as either a PEC (Dirichlet) or not. Dirichlet

k =4.3934 k =5 k =5.3551

Neumann

k =2.7096 k =3 k =3.3694

23

• 6. Applications

Landmine detection: near field inversions Real far-field 2D data inversions

24

Landmine detection

Carl Baum: ”... we detect everything, we identify nothing! ” Metal detectors : high rate of false alarms (non landmine artifacts)

?

sand air

• high cost (due to false alarms) :

USD 3 to buy, USD 200–1000 to clear

• requires high level of detection accuracy (deminers

safety) as opposed to military demining ≈ 100 million landmines world-wide ≈ 2000 victims per month

25

Humanitarian Demining Project

(HuMin/MD: http://www.humin-md.de) Our goal: Decrease the number of false alarms through fast new imaging algorithms. 1. Local 3D imaging: Karlsruhe, Mainz, Cologne, G¨

• ttingen, & des Saarlandes
• 2. Signal analysis
• 3. Hardware and soil

Our frequency domain approach:

• Factorization Method

(Kirsch, Grinberg, Hanke-Bourgeois)

• Linear Sampling Method

(Colton, Kirsch, Monk, Cakoni) (Multi-static/array data setting)

26

3D EM inversions: synthetic data

Multi-static measurement on 12 x 12 grid (40 x 40 cm) Frequency 1 kHz, k− = k+ ≈ 2.1 · 10−5, PEC objects

Reconstruction in perspective Zoomed reconstruction

27

2D inversions: synthetic data

Two-layered background. Frequency 10 kHz. Soil EM properties: σ− = 10−3 S/m, ϵ−

r = 10

k− ≈ 0.0063(1 + i) (δ = O(100m)) k+ ≈ 2.1 · 10−4 30 meas./source points along Γ = [−0.4, 0.4] × {0.05},

Two penetrable obstacles

−0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1

σD = 105 (high), ϵD

r = 8

U-shape metal

Linear sampling −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 Factorization −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1

σD = 106 (high) ϵr = 2.

28

Plastic only mine.

Linear sampling −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 Factorization −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1

σout = σin = 10−1 (weakly conductive) ϵin

r = 3, ϵout r

= 3 (plastic/TNT) Metal trigger.

Linear sampling −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 Factorization −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1

Further multiple PEC scatterers

−0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1 −0.4 −0.2 0.2 0.4 −0.3 −0.2 −0.1

29

Experimental 2D far-field data

Free-space parameters Frequency 10 GHz, λ = 3 cm, L = 15 cm Ipswich data (US Air Force Research Lab) Multi-static setting: 32 incident and measurement dir. Aluminum triangle Plexiglas triangle

FM −15 −10 −5 5 10 15 −15 −10 −5 5 10 15 FM −15 −10 −5 5 10 15 −15 −10 −5 5 10 15

Cavity

FM −15 −10 −5 5 10 15 −15 −10 −5 5 10 15

30

Remark

Superposition of the array data via

∫

Ω

u∞(ˆ x, d)g(d)ds(d) allows us to devise a criterion to determine whether a sampling point z belongs to the scatterer.

• This is done by testing the data against the back-

ground Green’s function (or dyadic in 3D) Φ(x, z) through a linear equation for each point z.

• Scattering data from an obstacle D is compatible

with the field due to a point source when z is in- side D and not compatible when z is outside D (ranges...)

References:

The factorization method for inverse problems (2008), Kirsch and Grinberg, Springer Qualitative methods in inverse scattering the-

• ry (2007), Cakoni and Colton , Springer

Inverse acoustic and EM scattering theory (2013), 3rd ed., Colton and Kress, Springer Stream of papers in Inverse problems journal

31

Recapping

Sampling methods

• No forward solver
• No a priori info on the scatterer
• No asymptotic approximation (full EM)
• Potentially fast
• Eigenfrequencies exploitable
• Robust within various settings

Drawbacks

• Too much data – multi-static setup
• Cannot easily incorporate extra info
• Does’t determine scatterer properties
• Needs background Green’s function

− Approximately − Greens tensor in 3D − Hankel transforms in the layered case

32