An inverse problem for a waveguide Virginia Selgas Peter Monk - - PowerPoint PPT Presentation

An inverse problem for a waveguide

Peter Monk

University of Delaware USA

Virginia Selgas

University of Oviedo Spain

PICOF 2012

Goal: non-destructive testing of pipes

Problem. We want to detect and locate bounded inhomogeneous obstacles in an infinite tubular waveguide, when we are given measurements of pressure waves due to point sources inside the waveguide. We have in mind a potential application of acoustic techniques for the inspection of underground pipes such as sewers.

thanks to S.N. Chandler-Wilde for bringing this to our attention! We assume that an acoustic source and array of microphones could be lowered into a sewer and then used to collect multistatic scattering data from the mouth of the pipe.

Previous work

Engineers have suggested using acoustic means to inspect sewers (but have not considered multistatic data).

[F . Podd, M. Ali, K. Horoshenkov, et al., Water Science and Technology, 56 (2007), pp. 131–9] [M.T.Bin Ali, K.V.Horoshenkov, S.J.Tait, NOVATECH 2010]

Linear Sampling Method (LSM)

[D. Colton and A. Kirsch, Inv. Prob., 12 (1996), pp. 383–93]

Reciprocity Gap Method (RGM)

[D. Colton and H. Haddar, Inv. Prob., 21 (2005), pp. 383–98]

LSM for infinite sound hard acoustic waveguides (i.e. pipes) and sound soft obstacle.

[L. Bourgeois and E. Lunéville, Inv. Prob., 24 (2008), pp. 15–8]

Several authors have considered the use of time-reversal methods for small inclusions in waveguides.

e.g. [P . Roux and M. Fink, J. Acoust. Soc. Am., 107 (2000), pp. 2418–29]

Forward problem: geometry

Geometry. Cross-section of the waveguide: Σ ⊂ R2, bounded, smooth, simply connected later a disc Infinite tubular waveguide: R × Σ ⊂ R3 Penetrable obstacle: D ⊂ R × Σ, bounded, with connected complement R × Σ \ D ∂D denotes the boundary of D We identify each point x ∈ R3 with (x1, ˆ x) ∈ R × R2.

✓ ✒ ✏ ✑ ✄ ✂

• ✁

✲ ❄

• ✒

{ }×Σ D ν0 ν νD

Forward problem: equations

A known incident field ui is excited in the pipe (e.g. by a loudspeaker). Forward problem. Find the fields u and v such that                          △u + k2u = 0 in (R × Σ) \ D , △v + k2 n v = 0 in D , v − u = ui across ∂D ∩ (R × Σ) , ∂νDv − ∂νDu = ∂νDui across ∂D ∩ (R × Σ) , ∂νu = 0

• n (R × ∂Σ) \ ∂D ,

∂νv = 0

• n (R × ∂Σ) ∩ ∂D ,

The inhomogeneity is represented by the coefficient n ∈ C(D) with Re(n) ≥ C > 0 and Im(n) ≥ 0 in D. The total field is ut = ui + u, being u the scattered field (as above). A further radiation condition required for uniqueness!

Forward problem: radiation condition

The radiation condition requires the scattered wave to propagate (or decay) away from the inhomogeneity D. To formalize it, we need: cross-sectional modes

kn ∈ [0, +∞), where k2

n is an eigenvalue of the Neumann problem

for the negative Laplacian on Σ, sorted in such way that kn ր +∞; θn an eigenfunction associated to the eigenvalue k2

n and with

{θn}n∈N forming an orthonormal basis of L2(Σ);

waveguide modes g±

n (x1, ˆ

x) := θn(ˆ x) e±ıβn x1, where βn :=

• k2 − k2

n ∈ C and we choose Re(βn) ≥ 0 and Im(βn) ≥ 0.

Note:

• bviously ∆g±

n + k2g± n = 0 in R × Σ, and ∂νg± n = 0 on R × ∂Σ;

there are a finite number of propagating modes, and the remainder are evanescent modes.

Forward problem: radiation condition formalization

The Dirichlet–to–Neumann operator is T ±R : H1/2(Σ±R) → ˜ H−1/2(Σ±R) given by T ±Rg = ±ı

• n∈N

βn

• Σ±R

g θn dS θn , where Σ±R := {±R} × Σ. The radiation condition can be written as T ±Rus = ∂ν0us

• n Σ±R .

In addition we assume that the incident field ui is a smooth solution of

• △ui + k2ui = 0

in (−RS, RS) × Σ , ∂νui = 0

• n (−RS, RS) × ∂Σ ,

for some RS ∈ (0, R) with D ⊂ (−RS, RS) × Σ.

Forward problem: summary

Forward problem. We consider R > 0 large enough so that D ⊂ (−R, R) × Σ := BR. The forward problem is to find u ∈ H1(BR \ D) and v ∈ H1(D) such that                          △u + k2 u = in BR \ D , △v + k2 n v = in D , v − u = ui across ∂D ∩ BR , ∂νDv − ∂νDu = ∂νDui across ∂D ∩ BR , ∂νu = 0 and ∂νv =

• n ((−R, R) × ∂Σ) ∩ ∂D ,

∂ν0u = T ±Ru

• n Σ±R .

Result. Existence and uniqueness except for a discrete set of frequencies (now assumed to be avoided).

Inverse problem: geometry

We consider two surfaces ΣS := ΣaS and ΣM := ΣaM where we place sources and take measurements, respectively. We assume that −R < aS < aM < R and D ⊂ (aM, R) × Σ , We denote BS

R := (aS, R) × Σ and BM R := (aM, R) × Σ, so that

D ⊂ BM

R ⊂ BS R ⊂ BR .

✓ ✒ ✏ ✑ ✓ ✒ ✏ ✑ ✓ ✒ ✏ ✑ ✓ ✒ ✏ ✑ ✄ ✂

• ✁

Σ−R ΣS ΣM ΣR D

Inverse problem: point sources

For each source point x0 ∈ ΣS, the incident field ui

x0 is

ui

x0(x) := Φ(x0, x)

with x0 ∈ ΣS, x ∈ R × Σ , where Φ is the fundamental solution for the waveguide Φ(x, y) := Φx(y) := −

• n∈N

eıβn|x1−y1| 2ıβn θn(ˆ x) θn(ˆ y) for x = (x1, ˆ x), y = (y1, ˆ y) ∈ R × Σ with x = y. For this incident field ui

x0, the corresponding total field is

ut

x0 = us x0 + ui x0 in BR \ D

and vx0 in D .

Inverse problem

Given ut

x0(x) and ∂ν0ut x0(x) for

all source points x0 ∈ ΣS and all measurement points x ∈ ΣM, determine the boundary of D. Remark. Of course in practice we would only have a finite number of sources and a finite number of measurements.

LSM and RGM algorithms

Two possible methods considered in our paper: Linear Sampling Method (LSM) It uses only measurements of ux0, but requires the fundamental solution of the pipe and entire domain (i.e. pipe and manhole). Reciprocity Gap Method (RGM) It uses all measurements, of ut

x0(x) and ∂ν0ut x0(x), but isolates the

pipe and only requires the fundamental solution of the pipe. For the remainder of the talk we shall concentrate on the RGM.

Reciprocity Gap operator: function spaces

For B ⊆ BR, define the function space H(B) =

• u ∈ H1(B) ; △u + k2 u = 0 in B,

∂νu = 0 on ∂B ∩ ((−R, R) × ∂Σ)} , as well as its subspace HR(B) := {u ∈ H(B) ; ∂ν0u = T Ru on ΣR} if ΣR ⊂ ∂B. We also define U := {ut

x0 = us x0 + ui x0 ;

x0 ∈ ΣS} .

Reciprocity Gap operator

For sources placed on ΣS and measurements made on ΣM, we consider the RG operator R : HR(BS

R) → L2(ΣS) given by

Rv(x0) := R(ut

x0, v)

for x0 ∈ ΣS, v ∈ HR(BS

R) ,

where the Reciprocity Gap (RG) bilinear form is R(ut

x0, v) :=

• ΣM
• ut

x0 ∂ν0v − ∂ν0ut x0 v

• dS

for x0 ∈ ΣS, v ∈ HR(BS

R) .

Interior Transmission Problem

To analyze the RGM we first need to consider the following ITP . Interior Transmission Problem. Find v1 ∈ H1(D), v2 ∈ H1(D) such that ITP :                    △v1 + k2 v1 = 0 in D , △v2 + k2 n v2 = 0 in D , v1 − v2 = G across ∂D ∩ BR , ∂νDv1 − ∂νDv2 = g across ∂D ∩ BR , ∂νv1 = ∂νv2 = 0

• n ((−R, R) × ∂Σ) ∩ ∂D .

for suitable data g and G. Remark. This is a generalized interior transmission problem because of the mixed boundary conditions.

Interior Transmission Problem: analysis

Using the arguments of

[D.Colton, L.Paeivaerinta and J.Sylvester,IPI, 1 (2007) 13-28]

we can rewrite ITP as a mixed 4th order boundary valued problem for w = v1 − v2. Noticing that w satisfies an additional natural boundary condition on the boundary ((−R, R) × ∂Σ) ∩ ∂D1, we can show:

Lemma

If Re(n) > C > 1 in D, the interior transmission problem is well posed for any k except for, at most, a discrete set of k values called transmission eigenvalues. If Im(n) > 0 there are no transmission eigenvalues.

1Personal communication of F.Cakoni

First properties of R

We use this lemma and follow

[F .Cakoni, M.Cayoeren and D.Colton, Inv. Prob., 24 (2008), 065016].

Lemma

Provided k is not a transmission eigenvalue and D = ∅, the RG

• perator R : HR(BS

R) → L2(ΣS) is injective.

Single Layer operator

We need a suitable parametrization for the RG operator argument, v. This is provided by a single layer operator using densities on yet another open surface ΣSL := ΣaSL with aSL ∈ (−R, aS) . We define SΣSL: ˜ H−1

/2(Σ SL)→H±R(BR\Σ SL)

and SΣSL: ˜ H−1

/2(Σ SL)→H1 /2(Σ SL)

by (SΣSLg)(x) =

• ΣSL Φ(x, y) g(y) dSy

for a.e. x ∈ BR \ ΣSL , (SΣSLg)(x) =

• ΣSL Φ(x, y) g(y) dSy

for a.e. x ∈ ΣSL .

Single Layer operator: properties

We show that SΣSL : ˜ H−1/2(ΣSL) → H1/2(ΣSL) is linear and continuous; defines an isomorphism; provides a dense set of fields in H(D) (when one avoids mixed eigenvalues for the Laplacian).

Reciprocity Gap Method (RGM)

We propose locating and reconstructing D by means of the RGM: Reciprocity Gap Method (RGM). Given z ∈ BM

R , we decide whether

z ∈ D

• r

z ∈ BM

R \ D

by studying a regularized solution fz ∈ L2(ΣSL) of the ill–posed integral equation (R ◦ SΣSL)fz = RΦz

• n ΣS .

This is solved for each of many sampling points z ∈ BM

R , and an

indicator function given by 1/fzL2(ΣSL) is used to locate the boundary ∂D.

Justification of RGM

In order to use the RGM we would like to know that the operator R ◦ SΣSL is suitable for Tikhonov regularization; the result of this procedure can be used as an indicator function.

Justification of RGM: regularization

Lemma

Suppose k is not a transmission eigenvalue or an eigenvalue for the mixed Laplacian. Then, for any −R < aSL ≤ aS ≤ aM < R, R ◦ SΣSL : L2(ΣSL) → L2(ΣS) defines a linear continuous operator that is compact and injective, and has a dense range. Consequence. The operator R ◦ SΣSL is suitable for Tikhonov regularization.

Justification of RGM: the indicator function

Theorem

Assume that k is not a tranmission eigenvalue and not a mixed eigenvalue for D. Then, for any ε > 0 there exists f ε

z ∈ L2(ΣSL) with

||(R ◦ SΣSL)f ε

z − RΦz(x0)||L2(ΣS) < ε .

In addition (i) For any z ∈ D the corresponding sequence of functions {SΣSLf ε

z ; ε > 0} ⊂ L2(ΣSL) converge in L2(ΣSL) as ε → 0.

(ii) For each z ∈ BM

R \ D and any sequence of functions

{f ε

z ; ε > 0} ⊂ L2(ΣSL) as above, it holds that

||f ε

z ||L2(ΣSL) → ∞

for ε → 0 . Consequence. The result of the RGM, f ε

z , can be used as an indicator function.

Forward synthetic data

We take the pipe to be a cylinder with axis of rotation along the x1-axis and radius 0.075 m. Three frequencies are investigated: f = 10KHz, 6KHz and 2KHz. The sound speed in the pipe is c = 343 m/s and the density of air is ρ = 1.2 kg/m3. The equation is slightly more general than earlier: div 1 ρ(x)∇u

• +

ω2 ρ(x)c(x)2 u = 0 with piecewise constant coefficients. We use the Ultra Weak Variational Formulation of Cessenat and Després with the DtN map to solve the forward problem. f k = ω/c # prop. modes # total modes wavelength 2KHz 3.6637 3 38 1.71 6KHz 10.9910 21 38 0.57 10KHz 18.1530 53 77 0.2

Forward synthetic data

In all cases the field is measured on ΣM at x1 = aM := −4m, while the sources are at x1 = aS := −4.1m. We use point sources positioned at Mq Gauss-Jacobi points radially and Nq in angle. We assume that the Fourier coefficients of the scattered field can be measured on ΣM. wave number k Mq Nq 3.6637 & 10.9910 4 12 18.1530 6 20

−0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 y2 y3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 y2 y3

Pipe with penetrable ball

There is a penetrable ball of radius 0.03m inside the pipe. In the ball c = 600 m/s and density ρ = 20 kg/m3. An isosurface of the indicator function 1/||fz||L2(ΣSL): Exact object Reconstruction at 10KHz Reconstruction at 6KHz Reconstruction at 2KHz

Pipe with penetrable ball: 6KHz

Contour plots of 1/||fz||L2(ΣSL) on several cross-sections:

Slice at z2= 0 z1 z3 −1 −0.5 0.5 1 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.02 0.04 0.06 0.08 0.1 0.12 0.14 slice at z1= 0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.02 0.04 0.06 0.08 0.1 0.12 0.14

z1, z3 plane at z2 = 0 z2, z3 plane at z1 = 0.2

slice at z1= 0 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.02 0.04 0.06 0.08 0.1 0.12 0.14 slice at z1= −0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.02 0.04 0.06 0.08 0.1 0.12 0.14

z2, z3 plane at z1 = 0 z2, z3 plane at z1 = −0.2

Pipe with penetrable ball: 10KHz

Contour plots of 1/||fz||L2(ΣSL) on several cross-sections:

Slice at z2= 0 z1 z3 −1 −0.5 0.5 1 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.005 0.01 0.015 0.02 0.025 0.03 0.035 slice at z1= 0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.005 0.01 0.015 0.02 0.025 0.03 0.035

z1, z3 plane at z2 = 0 z2, z3 plane at z1 = 0.2

slice at z1= 0 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.005 0.01 0.015 0.02 0.025 0.03 0.035 slice at z1= −0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.005 0.01 0.015 0.02 0.025 0.03 0.035

z2, z3 plane at z1 = 0 z2, z3 plane at z1 = −0.2

Pipe with hard obstruction

Here, ∂νu = 0 on the surface of the obstruction. An isosurface of the indicator function 1/||fz||L2(ΣSL): Exact object Reconstruction at 10KHz Reconstruction at 6KHz

Pipe with hard obstruction: 6KHz

Contour plots of 1/||fz||L2(ΣSL) on several cross-sections:

Slice at z2= 0 z1 z3 −1 −0.5 0.5 1 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.02 0.04 0.06 0.08 0.1 0.12 slice at z1= 0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.02 0.04 0.06 0.08 0.1 0.12

z1, z3 plane at z2 = 0 z2, z3 plane at z1 = 0.2

slice at z1= 0 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.02 0.04 0.06 0.08 0.1 0.12 slice at z1= −0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.02 0.04 0.06 0.08 0.1 0.12

z2, z3 plane at z1 = 0 z2, z3 plane at z1 = −0.2

Pipe with hard obstruction: 10KHz

Contour plots of 1/||fz||L2(ΣSL) on several cross-sections:

Slice at z2= 0 z1 z3 −1 −0.5 0.5 1 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 slice at z1= 0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

z1, z3 plane at z2 = 0 z2, z3 plane at z1 = 0.2

slice at z1= 0 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 slice at z1= −0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

z2, z3 plane at z1 = 0 z2, z3 plane at z1 = −0.2

Pipe with penetrable obstruction

Same parameters as for the penetrable ball. An isosurface of the indicator function 1/||fz||L2(ΣSL): Exact object Reconstruction at 10KHz Reconstruction at 6KHz

Pipe with penetrable obstruction: 6KHz

Contour plots of 1/||fz||L2(ΣSL) on several cross-sections:

Slice at z2= 0 z1 z3 −1 −0.5 0.5 1 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 slice at z1= 0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11

z1, z3 plane at z2 = 0 z2, z3 plane at z1 = 0.2

slice at z1= 0 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 slice at z1= −0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11

z2, z3 plane at z1 = 0 z2, z3 plane at z1 = −0.2

Blocked pipe

We consider a semi-infinite sound hard domain (−∞, 0) × Σ. Not covered by our theory! An isosurface of the indicator function 1/||fz||L2(ΣSL): Exact object Reconstruction at 6KHz The method fails to detect the end of the pipe!

Blocked pipe: 6Khz

The failure is not caused by the wrong choice of isosurface:

Slice at z2= 0 z1 z3 −1 −0.5 0.5 1 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 slice at z1= 0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38

z1, z3 plane at z2 = 0 z2, z3 plane at z1 = 0.2

slice at z1= 0 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 slice at z1= −0.2 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38

z2, z3 plane at z1 = 0 z2, z3 plane at z1 = −0.2

Blocked pipe (near data)

To test if measuring evanescent modes might help, we move the measurement surface to the near field at x1 = −0.05m: Exact Object Reconstruction at 6KHz (near data) Now the end of the pipe is detected!

Blocked pipe: 6KHz (near data)

The end of the pipe is clearly indicated:

Slice at z2= 0 z1 z3 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.05 0.1 0.15 0.2 0.25 0.3 slice at z1= 0.192 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.05 0.1 0.15 0.2 0.25 0.3

z1, z3 plane at z2 = 0 z2, z3 plane at z1 = 0.192

slice at z1= 0 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.05 0.1 0.15 0.2 0.25 0.3 slice at z1= −0.192 z2 z3 −0.6 −0.4 −0.2 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.05 0.1 0.15 0.2 0.25 0.3

z2, z3 plane at z1 = 0 z2, z3 plane at z1 = −0.192

Comments and conclusions

The inversion is rapid and easy to implement. The method can detect some scatterers in the pipe and close to the walls. We need to understand the scatterers that the method fails to find.2 Detecting surface imperfections, joints and cracks would be interesting.

2A modal solution of the blocked pipe problem confirms our observations

concerning the blocked pipe (Personal communication of L.Bourgeois).