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 u src for smooth background velocity v � 1 � v 2 ∂ 2 t − ∇ 2 u src ( x , t ) = δ ( x − x src ) δ ( t ) x Reflectivity: Velocity gets an oscillatory perturbation v ( x ) → v ( x )(1 + r ( x )) Scattered field: Treat this perturbation by linearization and let u scat be the perturbation of the field � 1 � u scat ( x , t ) = 2 r v 2 ∂ 2 t − ∇ 2 v 2 ∂ 2 t u src . x Inverse problem: Determine r ( x ) from the following data: u scat ( x , t ) for x = ( x 1 , x 2 , x 3 ) in a subset of x 3 = 0, t ∈ [0 , T max ]

Reverse Time Migration ◮ Numerically solve wave equation to obtain u src ( x , t ) = incoming (source) wave field u rt ( 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. � T I ( x ) = u src ( x , t ) u rt ( x , t ) dt (*) 0 or � I ( x ) = 1 u rt ( x , ω ) � u src ( x , ω ) d ω. 2 π � ◮ Analysis of Rakesh (1988) gives for (*) that I ( x ) = (ΨDO) r ( x ) ,

A simple example Medium total velocity v v r vnonlin v dv 0 0 0 1.5 1.5 100 100 100 0.08 200 200 200 0.06 1.4 1.4 300 300 300 0.04 400 400 400 0.02 1.3 = 1.3 + z 500 z 500 z 500 0 1.2 1.2 600 600 600 ! 0.02 700 700 700 ! 0.04 1.1 1.1 800 800 800 ! 0.06 900 900 900 ! 0.08 1 1 1000 1000 1000 0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000 x x x Simulated data u scat x Raypaths of waves 0 200 400 600 800 1000 0 dv x s 0 100 0.08 0.5 200 0.06 300 0.04 0.1 400 0.02 1.0 t , like expected 0 500 0 z 600 ! 0.02 from ray theory: -0.1 1.5 700 ! 0.04 800 ! 0.06 900 ! 0.08 2.0 1000 simulated data, xs=300 0 100 200 300 400 500 600 700 800 900 1000 x Problem : Invert the map F : r �→ u scat

Imaging with the adjoint Result of applying F ∗ on data from 5 sources x 0 200 400 600 800 1000 0 200 0.4 400 0.2 0 z -0.2 600 -0.4 -0.6 800 1000 imaging with the adjoint 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 � � u rt ( x , ω ) − ω − 2 v 2 ∇ � 1 u src ( x , ω ) � u src ( x , ω ) ∇ � u rt ( x , ω ) I ( x ) = d ω. u src ( x , ω ) | 2 2 π i ω | � (cf. Kiyashchenko et al., 2007) ◮ We prove a theorem that I ( x ) ≡ r ( x ) , microlocally, or, more precisely, I ( x ) = R ( x , D x ) 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 ) − 1 F ∗ d . Now F ∗ F = Ψ = pseudodifferential op. (PsDO) , � 1 R n e i � x ,ξ � S Ψ ( x , ξ ) � Ψ r = r ( ξ ) d ξ, (2 π ) n (Beylkin 1985, Ten Kroode 1998), with symbol S Ψ given by an integral involving geometrical optics rays from the scattering point x (x, , , ) ! " # x (x, , , ) ! " # s r ! " r S Ψ ( x , ξ ) = � ξ � n − 1 “ “ ” “ ”” ξ ξ B x , + B x , − , � ξ � � ξ � " # v 2 sin( θ ) n − 2 v ( x r ) v ( x s ) ZZ B ( x , ν ) = d θ d ψ GO rays 8(4 π ) 2 n − 1 cos( θ/ 2) 2 n − 1 cos( θ r ) cos( θ s ) ( x s , x r ) ∈ Acq x Requires extensive ray calculations → avoid by approximating ( F ∗ F ) − 1 using curvelets

Model problem Notation x = ( x 1 , x 2 , x 3 ), wave vector ξ = ( ξ 1 , ξ 2 , ξ 3 ) ∂ 2 ∂ t 2 2 Scattered from planar source, factor v 2 omitted ( 1 x ) u scat ( t , x ) = A δ ( t − x 3 v 2 ∂ 2 t − ∇ 2 v ) r ( x ) . � �� � u src Compute u scat , u rt in terms of Fourier transfrom � r � Av 2 2 − 2 iv � ξ � e ix · ξ − iv � ξ � t � u scat ( x , t ) = (2 π ) 3 Re r ( ξ − (0 , 0 , � ξ � )) d ξ for t s.t. u src has passed supp( r ) u rt ( x , t ) = same formula for all t !" (0,|| ! ||) in refl Wave vectors satisfy Snell’s law ! ξ outgoing wave number (0 , 0 , � ξ � ) incoming wave number ξ − (0 , 0 , � ξ � ) reflectivity wave number (0,|| ! ||)

Inversion: excitation time imaging condition Try 2 I ( x ) = v 2 A ( ∂ t + (0 , 0 , v ) · ∇ x ) u rt ( x , x 3 / v ) This gives � � � 2 1 − ξ 3 e ix · ( ξ − (0 , 0 , � ξ � )) � I ( x ) = (2 π ) 3 Re r ( ξ − (0 , 0 , � ξ � )) d ξ � ξ � R 3 Change of variables ˜ ξ = ξ − (0 , 0 , � ξ � ), !" (0,|| ! ||) � � in � � � ∂ ˜ ξ � = 1 − ξ 3 refl with � ξ � , and ∂ξ ! ξ ∈ R 2 × R < 0 domain ˜ Result is reconstruction except for ξ 3 = 0. (0,|| ! ||) � � 2 1 ξ ) e ix · ˜ ξ ) e ix · ˜ ξ d ˜ ξ d ˜ r (˜ r (˜ I ( x ) = (2 π ) 3 Re � ξ = � ξ. (2 π ) 3 R 2 × R < 0 R 2 × R � =0

Inversion: ratio imaging condition We had 2 I ( x ) = v 2 A ( ∂ t + (0 , 0 , v ) · ∇ x ) u rt ( x , x 3 / v ) Rewrite into modified ratio imaging condition ◮ Use that (0 , 0 , v ) = v 2 ∇ T ( x ), and that u src ( x , ω ) = Ae − i ω T ( x ) , � u src ( x , ω ) ≈ − i ω ∇ x T ( x ) Ae − i ω T ( x ) . ∇ x � ◮ Insert factors 2 v 2 ω 2 left out in model problem This results in � � u rt ( x , ω ) − ω − 2 v 2 ∇ � 1 u src ( x , ω ) � u src ( x , ω ) ∇ � u rt ( x , ω ) I ( x ) = d ω. u src ( x , ω ) | 2 2 π i ω | �

Tools for variable background: Microlocal analysis ◮ Study singularities of distributions with position and orientation (Hormander 1985, Duistermaat 1996, ...) L ◮ Fourier integral operators (FIO’s): operators of the form �� A ( x , y , θ ) e i φ ( x , y ,θ ) d θ u ( y ) dy . Fu ( x ) = Mapping of singularities according to the canonical relation Λ. Pseudodifferential operator ΨDO’s with symbol A ( x , ξ ) � 1 A ( x , D ) u = A ( x , ξ ) � u ( ξ ) d ξ (2 π ) n Identity canonical relation ◮ High-frequency wave packets mapped approximately the same y " ! x

Wave equation ◮ Use ray theory and Fourier integral operators ◮ Locally solution is given explicitly by �� 1 e i α ( y , t − s ,ξ ) a ( y , t − s , ξ ) � u ( y , t ) = f ( ξ, s ) ds d ξ + complex conj , (2 π ) 3 � �� � WKB solutions ◮ Global analysis: Propagation of singu- x =0 3 larities in ( x , t , ξ, ω ) space, on surface given by dispersion relation (y ,y ,t, , , ) ! ! " 1 2 1 2 � ξ � = v ( x ) − 1 | ω | t x 3 t=0 (x, ) #

Three kinds of fields Scattered field u scat Reverse time continued field u rt acquisition u inc u inc t>t S u scat u rt t<t S Continued scattered field u cont (unphysical, for purposes of the analysis only) u inc t>t S u cont t<t S

Analysis of the continued scattered field Characterization of the map r �→ u cont as a Fourier integral operator ◮ locally it is an explicit FIO �� 1 e i ϕ T ( y , t , x ,ξ ) A ( y , t , x , ξ ) r ( x ) d ξ dx . u cont ,ω< 0 ( y , t ) = (2 π ) 3 With phase function ϕ T ( y , t , x , ξ ) = α ( y , t − T ( x ) , ξ ) − ξ · x # ◮ Globally we have the FIO ! property with canonical rela- y " tion ( x , ζ ) �→ ( y , t , η, ω ) ac- t cording to Snell’s law and x propagation along rays || " || n s

Reverse time continued field: Boundary operators acquisition u inc t>t S u rt t<t S u rt is described by a reverse time wave equation with source N ( y , D y , D t )Ψ M ( y , t , D y , D t ) u scat , bdy ( y , t ) Here y 1 , y 2 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 Scattered field Reverse Time continued field acquisition u u inc inc t>t S u scat Describe the relation with u rt the continued scattered t<t S Continued scattered field field u inc t>t S u cont t<t S Theorem There are ΨDO’s Φ − , Φ + , such that u rt ( x , t ) = Φ − ( t , x , D x ) u cont , − + Φ + ( t , x , D x ) u cont , + in which u cont , ± is the part of u cont with ± ω > 0. To highest order acquisition Φ ± ( s , x , ξ ) = Ψ M ( y 1 , y 2 , t , η 1 , η 2 , ω ) if there is a reverse ray with ± ω > 0, connecting ( x , s , ξ, ω ) time cont. with ( y 1 , y 2 , t , η 1 , η 2 , ω )