Seismic inverse scattering by reverse time migration Chris Stolk - - PowerPoint PPT Presentation

Seismic inverse scattering by reverse time migration

Chris Stolk

University of Amsterdam

RICAM Workshop, November 21, 2011

Seismic imaging and linearized inversion

Reflection seismic imaging

◮ Claerbout’s imaging principle (1971), reverse time migration

(RTM) (Schultz and Sherwood, 1980).

◮ Mathematical treatment as a linearized inverse problem, using

microlocal analysis (Beylkin, 1985, Rakesh, 1988, Nolan and Symes 1997, Ten Kroode Et Al., 1998)

◮ Today: RTM as linearized inverse scattering

The linearized inverse scattering problem

Source field usrc for smooth background velocity v 1 v2 ∂2

t − ∇2 x

• usrc(x, t) = δ(x − xsrc)δ(t)

Reflectivity: Velocity gets an oscillatory perturbation v(x) → v(x)(1 + r(x)) Scattered field: Treat this perturbation by linearization and let uscat be the perturbation of the field 1 v2 ∂2

t − ∇2 x

• uscat(x, t) = 2r

v2 ∂2

t usrc.

Inverse problem: Determine r(x) from the following data: uscat(x, t) for x = (x1, x2, x3) in a subset of x3 = 0, t ∈ [0, Tmax]

Reverse Time Migration

◮ Numerically solve wave equation to obtain

usrc(x, t) = incoming (source) wave field urt(x, t) = reverse time continued receiver field

◮ Claerbout’s imaging principle: Reflectors exist at points in the

ground where the first arrival of the downgoing wave is time coincident with the upgoing wave. Implemented using so called imaging condition, e.g. I(x) = T usrc(x, t)urt(x, t) dt (*)

• r

I(x) = 1 2π

• urt(x, ω)
• usrc(x, ω) dω.

◮ Analysis of Rakesh (1988) gives for (*) that

I(x) = (ΨDO) r(x),

A simple example

Medium

total velocity

z x vnonlin 100 200 300 400 500 600 700 800 900 1000 100 200 300 400 500 600 700 800 900 1000 1 1.1 1.2 1.3 1.4 1.5

= v

x z v 100 200 300 400 500 600 700 800 900 1000 100 200 300 400 500 600 700 800 900 1000 1 1.1 1.2 1.3 1.4 1.5

+ v r

z x dv 100 200 300 400 500 600 700 800 900 1000 100 200 300 400 500 600 700 800 900 1000 !0.08 !0.06 !0.04 !0.02 0.02 0.04 0.06 0.08

Simulated data uscat

0.5 1.0 1.5 2.0 t 200 400 600 800 1000 x

simulated data, xs=300

• 0.1

0.1

, like expected from ray theory:

z x dv 100 200 300 400 500 600 700 800 900 1000 100 200 300 400 500 600 700 800 900 1000 !0.08 !0.06 !0.04 !0.02 0.02 0.04 0.06 0.08

xs Raypaths of waves

Problem: Invert the map F : r → uscat

Imaging with the adjoint

Result of applying F ∗ on data from 5 sources

200 400 600 800 1000 z 200 400 600 800 1000 x

imaging with the adjoint

• 0.6
• 0.4
• 0.2

0.2 0.4

F ∗Fr has discontinuities imaged at correct position but with incorrect amplitudes.

Inverse scattering with a single source

◮ We show that RTM can be used for inverse scattering, i.e. the

singularities of r are recovered

◮ New imaging condition

I(x) = 1 2πi ω usrc(x, ω) urt(x, ω) − ω−2v2 ∇ usrc(x, ω)∇ urt(x, ω) | usrc(x, ω)|2 dω. (cf. Kiyashchenko et al., 2007)

◮ We prove a theorem that

I(x) ≡ r(x), microlocally,

• r, more precisely,

I(x) = R(x, Dx)r, Here R is a pseudodifferential operator satisfying R(x, ξ) = 1 for “visible” reflectors (cf. Beylkin, 1985).

Multisource inversion: The normal operator

Inverse given by r = (F ∗F)−1F ∗d. Now

F ∗F = Ψ = pseudodifferential op. (PsDO), Ψr = 1 (2π)n

• Rn eix,ξSΨ(x, ξ)

r(ξ) dξ, (Beylkin 1985, Ten Kroode 1998), with symbol SΨ given by an integral involving geometrical optics rays from the scattering point

SΨ(x, ξ) = ξn−1 “ B “ x,

ξ ξ

” + B “ x, −

ξ ξ

”” ,

GO rays

! " #

s

x (x, , , ) ! " #

r

! " # x "r x (x, , , )

B(x, ν) = v2 8(4π)2n−1 ZZ

(xs,xr)∈Acq

dθ dψ sin(θ)n−2v(xr)v(xs) cos(θ/2)2n−1 cos(θr) cos(θs)

Requires extensive ray calculations → avoid by approximating (F ∗F)−1 using curvelets

Model problem

Notation x = (x1, x2, x3), wave vector ξ = (ξ1, ξ2, ξ3) Scattered from planar source, factor

∂2 ∂t2 2 v2 omitted

( 1 v2 ∂2

t − ∇2 x) uscat(t, x) = A δ(t − x3

v )

• r(x).

usrc Compute uscat, urt in terms of Fourier transfrom r uscat(x, t) = 2 (2π)3 Re

• Av2

−2ivξ eix·ξ−ivξt r(ξ − (0, 0, ξ)) dξ for t s.t. usrc has passed supp(r) urt(x, t) = same formula for all t Wave vectors satisfy Snell’s law ξ

• utgoing wave number

(0, 0, ξ) incoming wave number ξ − (0, 0, ξ) reflectivity wave number

! (0,||!||) in refl !"(0,||!||)

Inversion: excitation time imaging condition

Try I(x) = 2 v2A (∂t + (0, 0, v) · ∇x) urt(x, x3/v) This gives I(x) = 2 (2π)3 Re

• R3
• 1 − ξ3

ξ

• eix·(ξ−(0,0,ξ))

r(ξ − (0, 0, ξ)) dξ Change of variables ˜ ξ = ξ − (0, 0, ξ), with

• ∂ ˜

ξ ∂ξ

• = 1 − ξ3

ξ, and

domain ˜ ξ ∈ R2 × R<0

! (0,||!||) in refl !"(0,||!||)

Result is reconstruction except for ξ3 = 0. I(x) = 2 (2π)3 Re

• R2×R<0
• r(˜

ξ)eix·˜

ξ d ˜

ξ = 1 (2π)3

• R2×R=0
• r(˜

ξ)eix·˜

ξ d ˜

ξ.

Inversion: ratio imaging condition

We had I(x) = 2 v2A (∂t + (0, 0, v) · ∇x) urt(x, x3/v) Rewrite into modified ratio imaging condition

◮ Use that (0, 0, v) = v2∇T(x), and that

• usrc(x, ω) = Ae−iωT(x),

∇x usrc(x, ω) ≈ − iω∇xT(x)Ae−iωT(x).

◮ Insert factors 2v2 ω2 left out in model problem

This results in I(x) = 1 2π i ω usrc(x, ω) urt(x, ω) − ω−2v2 ∇ usrc(x, ω)∇ urt(x, ω) | usrc(x, ω)|2 dω.

Tools for variable background: Microlocal analysis

◮ Study singularities of distributions with position and

• rientation (Hormander 1985, Duistermaat 1996, ...)

L

◮ Fourier integral operators (FIO’s): operators of the form

Fu(x) =

• A(x, y, θ)eiφ(x,y,θ) dθ u(y) dy.

Mapping of singularities according to the canonical relation Λ. Pseudodifferential operator ΨDO’s with symbol A(x, ξ) A(x, D)u = 1 (2π)n

• A(x, ξ)

u(ξ) dξ Identity canonical relation

◮ High-frequency wave packets mapped approximately the same

! x y "

Wave equation

◮ Use ray theory and Fourier integral operators ◮ Locally solution is given explicitly by

u(y, t) = 1 (2π)3

• eiα(y,t−s,ξ)a(y, t − s, ξ)
• WKB solutions
• f (ξ, s) ds dξ + complex conj,

◮ Global analysis: Propagation of singu-

larities in (x, t, ξ, ω) space, on surface given by dispersion relation ξ = v(x)−1|ω|

! ! (y ,y ,t, , , ) "

1 2 1 2

t=0 x =0

3 3

x (x, ) t #

Three kinds of fields

Scattered field uscat u uscat

inc

Reverse time continued field urt

t>tS

S

inc

u

t<t

urt

acquisition

Continued scattered field ucont (unphysical, for purposes of the analysis only)

t<tS t>tS

u ucont

inc

Analysis of the continued scattered field

Characterization of the map r → ucont as a Fourier integral

• perator

◮ locally it is an explicit FIO

ucont,ω<0(y, t) = 1 (2π)3

• eiϕT (y,t,x,ξ)A(y, t, x, ξ) r(x) dξ dx.

With phase function ϕT(y, t, x, ξ) = α(y, t − T(x), ξ) − ξ · x

◮ Globally we have the FIO

property with canonical rela- tion (x, ζ) → (y, t, η, ω) ac- cording to Snell’s law and propagation along rays

t

x

y ! " #

s

n ||"||

Reverse time continued field: Boundary operators

t>tS

S

inc

u

t<t

urt

acquisition

urt is described by a reverse time wave equation with source N(y, Dy, Dt)ΨM(y, t, Dy, Dt)uscat,bdy(y, t) Here y1, y2 are boundary coordinates, and N and ΨM are pseudodifferential operators:

◮ ΨM is a cutoff operator with three effects

◮ Smooth taper near acquisition boundary ◮ Remove direct waves ◮ Remove tangently incoming waves

◮ N is to normalize the wave field to get correct amplitudes

Reverse time continued field vs. continued scattered field

Describe the relation with the continued scattered field

Scattered field

uscat

inc

u

t<tS t>tS inc

u

Continued scattered field Reverse Time continued field

inc

u ucont

t<tS

urt

t>tS

acquisition

Theorem There are ΨDO’s Φ−, Φ+, such that urt(x, t) = Φ−(t, x, Dx)ucont,− + Φ+(t, x, Dx)ucont,+ in which ucont,± is the part of ucont with ±ω > 0. To highest order Φ±(s, x, ξ) = ΨM(y1, y2, t, η1, η2, ω) if there is a ray with ±ω > 0, connecting (x, s, ξ, ω) with (y1, y2, t, η1, η2, ω)

acquisition reverse time cont.

Inverse scattering

Theorem Suppose v is such that there is no multipathing between source and reflection points. Consider the image I(x) defined by

I(x) = 1 2π

• Ω(ω)

iω | usrc(x, ω)|2

• usrc(x, ω)

urt(x, ω) − v(x)2 ω2 ∇ usrc(x, ω)∇ urt(x, ω)

• dω

Then I(x) = ΨDO r(x) = (R−(x, Dx) + R+(x, Dx)) r(x) To highest order R±(x, ζ) = Φ±(x, ξ) if ζ and ξ are related through ζ = ξ ± ω∇T(x).

y !

inc

! " x

#

Inverse scattering, proof

To prove this we first prove the thorem with a ray theoretical approximation for usrc Lemma The result of the theorem hold for I replaced by Iexct(x) = −1 Asrc(x)∂

− n+1

2

t

• ∂t + v2∇T · ∇x
• urt(x, T(x))

(*) Secondly we use that usrc = Asrc(x)∂

n−3 2

t

δ(t − T(x)) + lower order and that Asrc = 0, to convert (*) to an expression in usrc and urt

Example: Wave packets that can be reconstructed

Example 1: Velocity perturbation, reconstructed velocity perturbation, and three traces for background medium v = 2.0 + 0.001z (z in meters and v in km/s). Error 8-10 %.

Example: Horizontal reflector with receiver side caustics

500 1000 1500 2000 200 400 600 800 1000 1200 1400 1600 1800 2000 Velocity, rays and fronts x (m) z (m) 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4

Example 2: (a) A velocity model with some rays; (b) Simulated data, with direct arrival removed; (c) Velocity perturbation; (d) Reconstructed velocity perturbation. Error 0-10%

F ∗F: amplitude effects and filters

Effects of F ∗F can be summarized as

• 1. Amplitude decay with depth by spherical spreading
• 2. Frequency dependent weight ∼ (−∆)(n−1)/2
• 3. Energy in data depends on orientation of reflector!
• 4. Further position and angle dependence because rays are

curved

• 5. Filter due to finite bandwidth

Correct for 1 by factor √t in data, and for 2 by factor ω−(n−1)/4 in data in freq. domain → redefine F, F ∗. After corrections a real, selfadjoint order 0 PsDO remains, which acts like an “orientation and position dependent multiplication”

When orientation is a function of position

(Symes, 2007) Assume r is smooth in the direction normal to a vector field ν(x) for all x, then F ∗F F ∗d = SΨ(x, ν(x)) F ∗d, i.e. regular multiplication. Then SΨ(x, ν(x))−1 ≈ F ∗d F ∗F F ∗d , r = SΨ(x, ν(x))−1F ∗d. (with in first line a regularized division).

(x) x !

Generalize this to case of intersecting discontinuities in r

Curvelets

Special basisfunctions. Notation: φµ = φµ(x), µ ∈ M: curvelets, M index set. C: curvelet transform.

k1

• 0.25
• 0.50

0.25 0.50

• 0.50

k2

• 0.25

0.25 0.50 100 200 300 400 500 100 200 300 400 500

Samples Samples

Normal operator acting on curvelets

Normal operator becomes approximately diagonal in curvelet domain

◮ Take a curvelet with central frequency size λ. Approximate

• perator with multiplication by SΨ(xµ, νµ)

◮ Error due to variation of symbol over support of curvelet

positionspace ! !

1/2

!

!1

!

!1/2

fourierspace

◮ Using this we proved

• Ψφµ
• −

SΨ(xµ, νµ)φµ

• = O(λ−1/2).
• perator

multiplication

Diagonal approximation

◮ We therefore have a diagonal approximation

Ψr ≈ C TDSΨCr, where DSΨ is an diagonal matrix whose elements should approximate SΨ(xµ, νµ).

◮ We estimate an optimal DSΨ to satisfy this relation. Notes

◮ The problem is underdetermined: Regularization with

smoothness.

◮ Require positivity, this is used to fix the regularization

parameter

◮ Use a scaled version of F ∗d as vector r for more accuracy

Some results

Results are for a combination of this method with a curvelet based ℓ1 optimization including some anisotropic smoothing. Background medium and perturbation (from EAGE salt model)

Migrated image. Derived from this is the reference vector The imaged reference vector, and the result of the diagonal approximation are very close to each other

Final result after ℓ1 inversion with continuity promotion: good reconstruction Recovery remains good with noise present.

Final remarks

Single source

◮ We modified reverse time migration to become a linearized

inverse scattering method

◮ The modified imaging condition suppresses low-frequency

artifacts

◮ Thanks to my co-authors, Tim Op’t Root, and Maarten de

Hoop Multi-source

◮ We propose a position-wavenumber dependent scaling, as

• pposed to a diagonal scaling.

◮ Published with Peyman Moghaddam and Felix Herrmann in

ACHA (2008)

◮ Further methods in this direction: Nammour and Symes

(2009) and matrix probing by Demanet Et Al. (preprint)