Multiscale Analysis of Wave Propagation and Imaging in Random Media - - PowerPoint PPT Presentation

Multiscale Analysis of Wave Propagation and Imaging in Random Media Josselin Garnier (Ecole Polytechnique, France)

• In random media, correlation-based imaging aims at exploiting the information

carried by incoherent wave fluctuations when the coherent wave is vanishing.

• Fourth-order moment analysis makes it possible to understand when

correlation-based imaging is effective.

HKUST May 2019

Wave propagation in random media

• Wave equation:

1 c2( x) ∂2u ∂t2 (t, x) − ∆

xu(t,

x) = F(t, x),

• x = (x, z) ∈ R2 × R
• Time-harmonic source in the plane z = 0: F(t,

x) = δ(z)f(x)e−iωt. ֒ → u(t, x) = ˆ u( x)e−iωt and ˆ u satisfies ω2 c2( x) ˆ u( x) + ∆

xˆ

u( x) = −δ(z)f(x),

• x = (x, z) ∈ R2 × R
• Random medium model:

1 c2( x) = 1 c2

• 1 + µ(

x)

• co is a reference speed,

µ( x) is a zero-mean random process.

• The statistical properties of the wave field ˆ

u can be studied by multiscale analysis [1,2].

[1] G. Papanicolaou, SIAM J. Appl. Math. 21 (1971) 13. [2] J.-P. Fouque et al., Springer, 2007.

Wave propagation in the random paraxial regime

• Consider the time-harmonic form of the scalar wave equation (with

x = (x, z)): (∂2

z + ∆⊥)ˆ

u + ω2 c2

• 1 + µ(x, z)
• ˆ

u = −δ(z)f(x).

• The function ˆ

φ (slowly-varying envelope of a plane wave) defined by ˆ u(x, z) = ico 2ω ei ωz

co ˆ

φ

• x, z
• satisfies

∂2

z ˆ

φ +

• 2i ω

co ∂z ˆ φ + ∆⊥ ˆ φ + ω2 c2

• µ(x, z)ˆ

φ

• = 2i ω

co δ(z)f(x).

• If µ is mixing and smooth, if

ω → ω ε4 , µ(x, z) → ε3µ x ε2 , z ε2

• ,

f(x) → f x ε2

• ,

then ˆ φε converges in distribution in C0(R+, L2(R2)) (or C0(R+, Hk(R2))) as ε → 0 to the solution of: dˆ φ = ico 2ω ∆⊥ ˆ φdz + iω 2co ˆ φ ◦ dB(x, z) with B(x, z) Brownian field E[B(x, z)B(x′, z′)] = γ(x − x′) min(z, z′), γ(x) = ∞

−∞ E[µ(0, 0)µ(x, z)]dz, and ˆ

φ(z = 0, x) = f(x) [1].

[1] J. Garnier et al., Ann. Appl. Probab. 19 (2009) 318.

First-order moments in the paraxial regime Consider dˆ φ = ico 2ω ∆⊥ ˆ φdz + iω 2co ˆ φ ◦ dB(x, z) starting from ˆ φ(x, z = 0) = f(x).

• By Itˆ
• ’s formula,

d dz E[ˆ φ] = ico 2ω ∆⊥E[ˆ φ] − ω2γ(0) 8c2

• E[ˆ

φ] and therefore E ˆ φ(x, z)

• = ˆ

φhom(x, z) exp

• − γ(0)ω2z

8c2

• ,

where γ(x) = ∞

−∞ E[µ(0, 0)µ(x, z)]dz and ˆ

φhom is the solution in the homogeneous medium.

• Strong damping of the coherent wave.

= ⇒ Identification of the scattering mean free path Zsca =

8c2

• γ(0)ω2 [1].

= ⇒ Coherent imaging methods (such as Kirchhoff migration or Reverse-Time migration) fail.

[1] A. Ishimaru, Academic Press, 1978.

Second-order moments in the paraxial regime

• The mean Wigner transform defined by

W(r, ξ, z) =

• R2 exp
• − iξ · q
• E
• ˆ

φ

• r + q

2 , z ˆ φ

• r − q

2 , z

• dq,

is the angularly-resolved mean wave energy density. By Itˆ

• ’s formula, it solves a radiative transport equation

∂W ∂z + co ω ξ · ∇rW = ω2 4(2π)2c2

• R2 ˆ

γ(κ)

• W(ξ − κ) − W(ξ)
• dκ,

starting from W(r, ξ, z = 0) = W0(r, ξ), the Wigner transform of f. = ⇒ Identification of the scattering cross section

ω2 4c2

• ˆ

γ(κ) [1].

• The fields at nearby points are correlated and their correlations contain information

about the medium. = ⇒ One should use cross correlations for imaging in random media (CINT [2], ...).

[1] A. Ishimaru, Academic Press, 1978. [2] L. Borcea et al., Inverse Problems 21 (2005) 1419.

Remarks on fourth-order moments

• Physical conjecture: the wave field has Gaussian statistics; therefore we know

everything when the first two moments are characterized. → the conjecture may be wrong.

• The statistical second-order moments and the mean Wigner transform are not
• bserved directly, only empirical quantities are observed.
• Calculations of fourth-order moments are useful to:
• test the Gaussian conjecture.
• quantify the statistical stability of empirical second-order moments, Wigner

transforms, and correlation-based imaging methods.

• implement intensity-correlation-based imaging methods (when only intensities can

be measured, as in optics).

HKUST May 2019

Fourthd-order moments in the random paraxial regime

• Consider

dˆ φ = ico 2ω ∆⊥ ˆ φdz + iω 2co ˆ φ ◦ dB(x, z) starting from ˆ φ(x, z = 0) = f(x).

• Let us consider the fourth-order moment:

M4(r1, r2, q1, q2, z) = E

• ˆ

φ r1 + r2 + q1 + q2 2 , z ˆ φ r1 − r2 + q1 − q2 2 , z

• ׈

φ r1 + r2 − q1 − q2 2 , z ˆ φ r1 − r2 − q1 + q2 2 , z

• By Itˆ
• ’s formula,

∂M4 ∂z = ico ω

• ∇r1 · ∇q1 + ∇r2 · ∇q2
• M4 + ω2

4c2

• U4(q1, q2, r1, r2)M4,

with the generalized potential U4(q1, q2, r1, r2) = γ(q2 + q1) + γ(q2 − q1) + γ(r2 + q1) + γ(r2 − q1) −γ(q2 + r2) − γ(q2 − r2) − 2γ(0). = ⇒ One can get a general characterization of the fourth-order moment [1].

[1] J. Garnier et al., ARMA 220 (2016) 37.

Scintillation Assume that f(x) = exp

• − |x|2

2r2

• .
• The scintillation index:

S(r, z) := E

• ˆ

φ(r, z)

• 4

− E

• ˆ

φ(r, z)

• 22

E

• ˆ

φ(r, z)

• 22

satisfies (in the paraxial regime): S(r, z) = 1 − 1

• 1

4π

• R2 exp
• ω2

4c2

• z

0 γ

• u coz′

ωro

• dz′ − |u|2

4

+ iu · r

ro + |r|2 r2

• du
• 2 .

The physical conjecture is that S ≃ 1 when the propagation distance is larger than the scattering mean free path, as it should be for a (complex) Gaussian process.

HKUST May 2019

Scintillation

1 2 3 4 5 6 0.2 0.4 0.6 0.8 1 z / Zsca scintillation index Zc / Zsca=0.125 Zc / Zsca=0.25 Zc / Zsca=0.5 Zc / Zsca=1

Scintillation index at the beam center S(0, z) as a function of the propagation distance for different values of Zsca =

8c2

• ω2γ(0) and Zc = ωroℓc

co . Here

γ(x) = γ(0) exp(−|x|2/ℓ2

c).

→ The physical conjecture that S ≃ 1 when z ≫ Zsca is true in the random paraxial regime.

HKUST May 2019

Stability of the Wigner transform of the field

• The Wigner transform

W(r, ξ, z) :=

• R2 exp
• − iξ · q

ˆ φ

• r + q

2 , z ˆ φ

• r − q

2 , z

• dq

is not statistically stable (i.e. standard deviation > mean).

• Let us consider the smoothed Wigner transform (for rs, ξs > 0):

Ws(r, ξ, z) = 1 (2π)2r2

s ξ2 s

• R2×R2 W(r − r′, ξ − ξ′, z) exp
• − |r′|2

2r2

s

− |ξ′|2 2ξ2

s

• dr′dξ′.

Its coefficient of variation: Cs(r, ξ, z) :=

• E[Ws(r, ξ, z)2] − E[Ws(r, ξ, z)]2

E[Ws(r, ξ, z)] determines its statistical stability. ֒ → Analysis of high-order moments of ˆ φ [1].

[1] J. Garnier et al., ARMA 220 (2016) 37.

Stability of the Wigner transform of the field

ξs

0.5 1 1.5 2

rs

0.5 1 1.5 2

0.33 0.5 . 5 . 7 5 0.75 0.75 1 1 1 1.25 1.25 1.5 1.5 2 2 4 4

Contour levels of the coefficient of variation of the smoothed Wigner transform. Here rs = rs/ρ, ξs = ξsρ, and ρ = ρ(z; ω, ro, ℓc, Zsca). → This result makes it possible to achieve optimal trade-off between stability and resolution for correlation-based imaging [1,2]. Example: ultrasonic non-destructive testing in concrete.

[1] L. Borcea et al., Inverse Problems 27 (2011) 085004. [2] J. Garnier et al., ARMA 220 (2016) 37.

Wave propagation in random waveguides z x d/2 −d/2 Wave propagation in two-dimensional waveguides:

• (∂2

x + ∂2 z) + k2n2(x, z)

• ˆ

u(x, z) = δ(z)fs(x)

• Ideal waveguide:

n(0)(x) =    n if x ∈ (−d/2, d/2), 1

• therwise,

where n > 1 is the relative index of the core and d > 0 is its diameter.

• Type I perturbation: the index of refraction within the core is randomly perturbed:

n(ε)(x, z) =    n + εν(x, z) if x ∈ (−d/2, d/2) and z ∈ (0, L/ε2), 1

• therwise.
• Type II perturbation: the boundaries of the core are randomly perturbed:

n(ε)(x, z) =    n if x ∈

• − d/2 + εdν−(z), d/2 + εdν+(z)
• and z ∈ (0, L/ε2),

1

• therwise.

HKUST May 2019

Ideal waveguide

• The Helmholtz operator ∂2

x + k2n(0)(x) has spectrum (−∞, k2) ∪ {β2 N−1, . . . , β2 0}.

• the N modal wavenumbers βj are positive and k2 < β2

N−1 < · · · < β2 0 < n2k2.

• φj, j = 0, . . . , N − 1, are the eigenfunctions associated to β2

j .

• φ(t)

γ , t ∈ {e, o} are the generalized eigenfunctions associated to γ ∈ (−∞, k2).

• ˆ

u can be expanded on the complete set of orthonormal eigenmodes: ˆ u(x, z) =

N−1

• j=0

uj(z)φj(x) +

• t∈{e,o}

k2

−∞

u(t)

γ (z)φ(t) γ (x)dγ

The mode amplitudes satisfy for z = 0: ∂2

zuj + β2 j uj = 0,

j = 0, . . . , N − 1, ∂2

zu(t) γ + γu(t) γ

= 0, γ ∈ (−∞, k2) Therefore, for z > 0: ˆ u(x, z) =

N−1

• j=0

aj,s

• βj

eiβjzφj(x)

• guided mode

+

• t∈{e,o}

k2 a(t)

γ,s

γ1/4 ei√γzφ(t)

γ (x)

• radiating mode

dγ +

−∞

a(t)

γ,s

|γ|1/4 e−√

|γ|zφ(t) γ (x)

• evanescent mode

dγ The coefficients a·,s are constant and determined by the source.

HKUST May 2019

Random waveguide

• ˆ

u can be expanded on the complete set of eigenmodes of the ideal waveguide: ˆ u(x, z) =

N−1

• j=0

uj(z)φj(x) +

• t∈{e,o}

k2

−∞

u(t)

γ (z)φ(t) γ (x)dγ

→ uj, u(t)

γ

satisfy coupled equations.

• For large propagation distances:

ˆ u(x, z ε2 ) ≃

N−1

• j=0

cp(z)eiβj z

ε2 φj(x)

and the moments of cp satisfy closed-form equations. This result is obtained in the limit ε → 0, in distribution, in C0(R+, L2

weak(R)). HKUST May 2019

First- and second-order moments in random waveguides

• E[cp] decays exponentially with z (with a mode-dependent scattering mean free

path ℓp),

• (E[|cp|2])N

p=1 satisfies the coupled system

d dz E[|cp|2] =

• q=p

Γpq(E[|cq|2] − E[|cp|2]) − ΛpE[|cp|2] where Γjl = k4 2βjβl ∞ Rjl(z) cos

• (βl − βj)z
• dz,

Rjl(z) =   

• R
• R φjφl(x)E[ν(x, 0)ν(x′, z)]φjφl(x′)dxdx′

type I (n2 − 1)2d2 φ2

jφ2 l

• − d

2

• E[ν−(0)ν−(z)] + φ2

jφ2 l

d

2

• E[ν+(0)ν+(z)]
• type II

Λj = k2 k4 2√γβj

• t∈{e,o}

∞ R(t)

j,γ(z) cos

• (√γ − βj)z
• dzdγ

→ exchange of energy between the guided modes, → loss of energy towards the radiating modes.

HKUST May 2019

Fourth-order moments in random waveguides

• The moments

Rjl(z) = E

• |cj|2(z)|cl|2(z)
• ,

j, l = 0, . . . , N − 1, satisfy the closed equations ∂zRjj = − 2ΛjRjj +

• n=j

Γjn(4Rjn − 2Rjj), ∂zRjl = − (2Γjl + Λj + Λl)Rjl +

• n=l

Γln(Rjn − Rjl) +

• n=j

Γjn(Rnl − Rjl), j = l.

• Analysis of the linear system: Var(|cp|2)/E[|cp|2]2 grows exponentially with z !
• Remark: When N → ∞, the exponential growth rate goes to zero.
• This may happen in underwater acoustics (Pekeris waveguide).
• Correlation-based imaging is challenging in such conditions.

HKUST May 2019

Remarks on fourth-order moments

• Calculations of fourth-order moments are useful to:
• test the Gaussian conjecture.
• quantify the statistical stability of empirical second-order moments, Wigner

transforms, and correlation-based imaging methods.

• implement intensity-correlation-based imaging methods (when only intensities can

be measured, as in optics).

HKUST May 2019

Speckle intensity correlation imaging through a scattering medium − →

Speckle pattern

Experimental set-up [1]

• The light source is a time-harmonic plane wave.
• The object to be imaged is a mask that can be shifted transversally.
• For each position of the object the spatial intensity of the transmitted field (speckle

pattern) can be recorded by the camera.

[1] J. A. Newmann et al., Phys. Rev. Lett. 113 (2014) 263903.

Speckle intensity correlation imaging through a scattering medium

• The field just after the mask (in the plane z = 0) is (for a transverse shift r):

fr(x) = f(x − r), where f is the indicator function of the mask.

• The field in the plane of the camera (in the plane z = L) is denoted by ˆ

φr(x).

• The measured intensity correlation is

Cr,r′ = 1 |A0|

• A0

|ˆ φr(x)|2|ˆ φr′(x)|2dx − 1 |A0|

• A0

|ˆ φr(x)|2dx 1 |A0|

• A0

|ˆ φr′(x)|2dx

• ,

where A0 is the spatial support of the camera.

HKUST May 2019

Speckle intensity correlation imaging through a scattering medium

• Result (in the paraxial regime):

E

• Cr,r′

=

• A0

dX

• dY
• 1

(2π)2 f

• x + r′ − r

2

• f
• x − r′ − r

2

• exp
• − iζ · x
• dx
• × exp
• iζ ·
• X − r + r′

2

• exp

ω2 4c2

• L

γ coζ ω z − Y

• − γ(0)dz
• dζ
• 2

−

• 1

(2π)2

• A0

dX f

• x + r′ − r

2

• f
• x − r′ − r

2

• exp
• − iζ · x
• dx
• × exp
• iζ ·
• X − r + r′

2

• exp
• − ω2

4c2

• γ(0)L
• dζ
• 2

, with γ(x) = ∞

−∞ E[µ(0, 0)µ(x, z)]dz. HKUST May 2019

Speckle intensity correlation imaging through a scattering medium

• Result: When L ≫ Zsca :=

8c2

• γ(0)ω2 and

|A0|(∼ diam(camera)2) ≫ ρ2

L := Zscaℓ2 c

L (ρL = speckle radius), we have Cr,r′ ≃ E

• Cr,r′

≈

• | ˆ

f(κ)|2 exp

• iκ · (r′ − r)
• dκ
• 2

, up to a multiplicative constant, where ˆ f(κ) =

• f(x) exp
• − iκ · x
• dx.

֒ → It is possible to reconstruct the mask indicator function f.

HKUST May 2019

Speckle intensity correlation imaging through a scattering medium

• We have

Cr,r′ ≃ E

• Cr,r′

≈

• | ˆ

f(κ)|2 exp

• iκ · (r′ − r)
• dκ
• 2

֒ → It is possible to reconstruct the incident field f by a two-step phase retrieval algorithm (Gerchberg-Saxon-type). 1) Given Cr,r′, we know the modulus of the (I)FT of | ˆ f(κ)|2, and we know the phase

• f | ˆ

f(κ)|2 (zero) → we can extract | ˆ f(κ)|2. 2) Given | ˆ f(κ)|2, we know the modulus of the FT of f(x), and we know the phase of f(x) (zero) → we can extract f(x). → not very stable.

HKUST May 2019

Speckle intensity correlation imaging through a scattering medium (II) Experimental set-up

• A laser beam with incident angle θ is shined on the scattering medium.
• The object to be imaged is a mask.
• The total intensity of the light that goes through the mask is collected by a bucket

detector. → For each incident angle θ the total transmitted intensity Eθ is measured.

HKUST May 2019

Speckle intensity correlation imaging through a scattering medium (II) Consider: C(∆θ) = 1 Θ

• Θ

EθEθ+∆θdθ − 1 Θ

• Θ

Eθdθ 2

• Result (in the paraxial regime):

E[C(∆θ)] = 1 (2π)2

• exp

ω2 2c2

• L

γ

• x + ∆θ(z + Lo)
• dz
• e−ix·κ|

f 2(κ)|2dκdx × exp

• − ω2γ(0)L

2c2

• − |

f 2(0)|2 exp

• − ω2γo(0)L

2c2

• ,

with γ(x) = ∞

−∞ E[µ(0, 0)µ(x, z)]dz. HKUST May 2019

Speckle intensity correlation imaging through a scattering medium (II)

• Result: When L ≫ Zsca :=

8c2

• γ(0)ω2 , ρ2

L = Zscaℓ2

c

L

is small enough and Lo is large enough, then C(∆θ) ≃ E[C(∆θ)] ≈

• |

f 2(κ)|2 exp

• − iκ∆θLo
• dκ

֒ → One can extract | f 2(κ)|2 (by Fourier transform or FFT) and then f 2(x) by a (one-step) phase retrieval algorithm.

HKUST May 2019

Speckle intensity correlation imaging through a scattering medium (III)

• Noise source (laser light passed through a rotating glass diffuser).
• without object in path 1; a high-resolution detector measures the spatially-resolved

intensity I1(t, x).

• with object (mask) in path 2; a single-pixel detector measures the

spatially-integrated intensity I2(t). Experiment: the correlation of I1(·, x) and I2(·) is an image of the mask [1,2].

[1] A. Valencia et al., PRL 94 (2005) 063601; [2] J. H. Shapiro et al., Quantum Inf. Process 1 (2012) 949.

Speckle intensity correlation imaging through a scattering medium (III)

• Wave equation in paths 1 and 2:

1 c2

j(

x) ∂2uj ∂t2 − ∆

xuj = e−iωotn(t, x)δ(z) + c.c.,

• x = (x, z) ∈ R2 × R,

j = 1, 2

• Noise source (with mean zero):
• n(t, x)n(t, x′)
• = F(t − t′) exp
• − |x|2

r2

• δ(x − x′)
• Wave fields: uj(t,

x) = vj(t, x)e−iωot + c.c., j = 1, 2

• Intensity measurements:

I1(t, x) = |v1(t, (x, L))|2 in the plane of the high-resolution detector I2(t) =

• R2 |v2(t, (x′, L + L0))|2dx′ in the plane of the bucket detector
• Correlation:

CT (x) = 1 T T I1(t, x)I2(t)dt − 1 T T I1(t, x)dt 1 T T I2(t)dt

• HKUST

May 2019

Speckle intensity correlation imaging through a scattering medium (III)

• If the propagation distance is larger than the scattering mean free path, then

CT (x)

T →+∞

− →

• R2 H(x − y)f(y)4dy,

where f(x) is the mask indicator function and H(x) is a convolution kernel [1].

• If the medium is homogeneous:

H(x) = r4

• 28π2L4 exp
• − |x|2

4ρ2

gi0

• ,

ρ2

gi0 = c2

• L2

2ω2

• r2
• .
• If the medium in path 1 and 2 is random (independent realizations):

H(x) = r4

• ρ2

gi0

28π2L4ρ2

gi1

exp

• − |x|2

4ρ2

gi1

• ,

ρ2

gi1 = ρ2 gi0 +

4c2

• L3

3ω2

• Zscaℓ2

c

. ֒ → Scattering (only slightly) reduces the resolution.

• If the medium in path 1 and 2 is random (same realization):

H(x) = r4

• 28π2L4 exp
• − |x|2

4ρ2

gi2

• ,

1 ρ2

gi2

= 1 ρ2

gi0

+ 16L Zscaℓ2

c

. ֒ → Scattering enhances the resolution !

[1] J. Garnier, Inverse Problems and Imaging 10 (2016) 409.

Speckle intensity correlation imaging through a scattering medium (IV)

• The medium in path 2 is randomly heterogeneous.
• There is no other measurement than I2(t).
• The realization of the source is known (use of a Spatial Light Modulator) and the

medium is taken to be homogeneous in the “virtual path 1” → one can compute the field (and therefore its intensity I1(t, x)) in the “virtual” output plane of path 1. ֒ → a one-pixel camera can give a high-resolution image of the mask!

HKUST May 2019

Final remark: On the role of the random medium Is random medium good or bad for imaging ? Random medium in region 0 is good. Random medium in regions 1 and 2 is bad. Random medium in region 3 plays no role.

HKUST May 2019

Conclusion

• Correlation-based imaging allows for imaging in randomly scattering media.
• One needs to process well-chosen cross correlations of the data.
• Fourth-order moment of the wave field is useful.
• First application: Scintillation index and stability of Wigner transform.
• Second application: Intensity correlation-based imaging, ghost imaging.

HKUST May 2019