a quick introduction to seismic imaging
play

A Quick Introduction to Seismic Imaging Alison Malcolm Memorial - PowerPoint PPT Presentation

A Quick Introduction to Seismic Imaging Alison Malcolm Memorial University of Newfoundland September 12, 2017 A Sense of Scale Terrabytes of data! Images from:


  1. A Quick Introduction to Seismic Imaging Alison Malcolm Memorial University of Newfoundland September 12, 2017

  2. A Sense of Scale Terrabytes of data! Images from: http://www.offshoreenergytoday.com/pgs-seismic-vessel-transfers-data-to-shore-via-12-mbits-link/ and https://www.pgs.com/marine-acquisition/tools-and-techniques/the-fleet/flexible-capacity/sanco-sword/

  3. Data Acquisition Images courtesy Prof. Jeremy Hall

  4. Schematic View of Wave Propagation (measurement movie) From: https://giphy.com/gifs/refraction-pbGzTWOkabGFy/download

  5. Why is this hard? material c(x) sedimentary rocks: 2-4 km/s igneous rocks: 2-7 km/s metamorphic rocks: 1-4 km/s at depth: 8-10 km/s Stolk & Symes (2004) travel distance: tens of wavelengths

  6. Mathematical Model e.g. Achenbach (73), Landau & Lifshitz (86), Aki & Richards (02) Conservation of momentum (F = ma): ρ Dv j Dt = ρ f j + ∂ i σ ij Da Dt = ∂ a ∂ t + v · ∇ a Hooke’s Law (linearly elastic, isotropic material): σ ij = λǫ kk δ ij + 2 µǫ ij σ ij stress tensor � � ∂ x j + ∂ u j ǫ ij = 1 ∂ u i strain tensor 2 ∂ x i

  7. Mathematical Model General Assumptions: • long wavelength compared to amplitude • smooth displacement For Today: • linear elasticity • constant density • isotropy

  8. Mathematical Model Elastic Wave Equation: ρ∂ 2 u j ∂ t 2 = ( λ + µ ) ∂ j ∂ k u k + µ ∇ 2 u j

  9. Mathematical Model Elastic Wave Equation: ρ∂ 2 u j ∂ t 2 = ( λ + µ ) ∂ j ∂ k u k + µ ∇ 2 u j Helmholtz decomposition: u = ∇ φ + ∇ × ψ

  10. Mathematical Model Elastic Wave Equation: ρ∂ 2 u j ∂ t 2 = ( λ + µ ) ∂ j ∂ k u k + µ ∇ 2 u j Helmholtz decomposition: u = ∇ φ + ∇ × ψ ∂ 2 t φ = c 2 p ∇ 2 φ ∂ 2 t ψ = c 2 s ∇ 2 ψ � c p = ( λ + 2 µ ) /ρ � See Aki & Richards 2002 book c s = µ/ρ

  11. Acoustic Simplification Acoustic (really P-wave only) assumption ∇ 2 φ − 1 c 2 ∂ 2 t φ = f φ = 0 t < 0 ∂ z φ | z=0 = 0 Reasons: • sources and receivers often in the ocean • computational cost

  12. Contrast Formulation Acoustic (really P-wave only) assumption L φ := ∇ 2 φ − 1 c 2 ∂ 2 t φ = f Linearize: c(x) = c 0 (x) + δ c(x) L φ = f L 0 φ 0 = f note that L 0 and φ 0 use c 0 (x)

  13. Contrast Formulation Acoustic (really P-wave only) assumption L φ := ∇ 2 φ − 1 c 2 ∂ 2 t φ = f Linearize: c(x) = c 0 (x) + δ c(x) L φ = f L 0 φ 0 = f note that L 0 and φ 0 use c 0 (x) subtract L o δφ = δ L φ Symes 09 and Stolk 00 give estimates on linearization error

  14. Contrast formulation Born approximation L 0 δφ = δ L φ 0 2 ∇ 2 δφ − 1 t δφ = 2 δ c ∂ 2 ∂ 2 t φ 0 c 3 c 0 0 δφ is called the scattered field

  15. Separation of Scales ∇ 2 δφ − 1 t δφ = 2 δ c c 02 ∂ 2 c 03 ∂ 2 t φ 0 We assume that on the scale of the wavelength: • δ c is oscillatory s r • c 0 is smooth G 0 (x , t , s) G 0 (r , t , x) V(x)

  16. Data Model ∇ 2 δφ − 1 t δφ = 2 δ c c 02 ∂ 2 c 03 ∂ 2 t φ 0 Assume a δ -source � G 0 (r , ω, x) 2 δ c(x) ω 2 G 0 (x , ω, s)dx δφ (s , r , ω ) = − c 0 (x) 3 X s r � �� � V(x) G 0 (x , t , s) G 0 (r , t , x) V(x)

  17. ‘Typical’ Processing Steps Filtering/Signal Processing/Geometry/Statics 1 (we will ignore these steps) Velocity analysis, i.e. find c 0 , 2 usually via iterative imaging Imaging (a.k.a. Migration), 3 i.e. find δ c From: Fehler & Larner TLE, 2008, http://tle.geoscienceworld.org/content/27/8/1006 and Spectrum http://www.spectrumgeo.com/imaging-services/land-environment/depth-processing/pre-stack-depth-migration

  18. Imaging Methods Increasing complexity ⇒ Stacking (aka averaging) 1 Kirchhoff Migration 2 One-way methods 3 Reverse-time methods 4

  19. Velocity Analysis Methods Increasing complexity ⇒ Normal Moveout Analysis/Semblance 1 Iterative Kirchhoff Migration 2 Iterative One-way methods 3 Iterative Reverse-time methods 4 –Full-Waveform Inversion (FWI)

  20. Ray Tracing (very briefly) Assume solution form: A k (x , t) φ 0 (x , t) = e i ωψ (x , t) � (i ω ) k k • A k , and ψ SMOOTH • e i ωψ (x , t) oscillatory • remove frequency dependence Developed by Wentzel, Kramers, Brillouin, independently in 1926 and by Jeffreys in 1923; see Appendix E of Bleistein, Cohen and Stockwell for more details.

  21. Ray Tracing (very briefly) Apply Helmholtz to A k (x , t) φ 0 (x , t) = e i ωψ (x , t) � (i ω ) k k Eikonal equation: 1 ( ∇ ψ ) 2 = c(x) 2 Transport equations: 2 ∇ ψ · A k + A k ∇ 2 ψ = 0

  22. Ray Tracing (very briefly) Eikonal equation: 1 ( ∇ ψ ) 2 = c(x) 2 Transport equations: 2 ∇ ψ · A k + A k ∇ 2 ψ = 0 Nonlinear! 1 Method of characteristics (ray-tracing) 2 Usually use Runge-Kutta (requires smooth c) 3 Note that for constant velocity rays are straight lines 4

  23. Simplest Approach Stacking and Normal-Moveout Analysis � � x � 2 t(x) = 2 h 2 1 + c 2 Correct for the x-dependence �� x � 2 � t(0) � 2 t corr (x) = t meas (x) − 2 + c 2 2 This is called the Normal Moveout Correction (NMO).

  24. NMO Velocity Analysis Minimizing ∂ x t(x) = 0 gives a measure of velocity.

  25. NMO Velocity Analysis

  26. NMO Velocity Analysis

  27. A more refined method Differential Semblance Minimize: � J h = 1 h 2 I 2 (x , z , h)dxdzdh 2 where I(x , z , h) is the NMO-corrected data before stacking (averaging). DSO (Differential Semblance Optimization) has been extensively studied by Symes and co-authors. Of particular importance are Santosa & Symes (1986) and Shen & Symes (2008)

  28. Kirchhoff Migration WKBJ Modeling � � G 0 (r , t − t 0 , x)2 δ c(x) c 0 (x) 2 ∂ 2 δφ (s , r , t)= t G 0 (x , t 0 , s)dxdt 0 X T � A(x , s , ω )e i ω (t − T(x , s)) d ω G 0 (x , t 0 , s) ≈ B(x , r , s ,ω ) � �� � � � A(x , s , ω )A(r , x , ω )2 δ c(x) ω 2 c 0 (x) 2 e i ω (t − T(s , x) − T(x , r)) dxd ω δφ (s , r , t)= X R

  29. Kirchhoff Migration WKBJ Modeling Formula � � ω 2 B(x , r , s , ω )e i ω (t − T(s , x) − T(x , r)) dxd ω δφ (s , r , t) = X R Assume B independent of ω � B(x , r , s) δ ′′ (t − T(s , x) − T(x , r))dx δφ (s , r , t) = X This is a Generalized Radon Transform Note that a Radon transform is often called a τ -p transform in exploration seismology.

  30. Kirchhoff Migration WKBJ Modeling Formula � B(x , r , s) δ ′′ (t − T(s , x) − T(x , r))dx δφ (s , r , t) = X distance (m) receiver position (m) 0 500 1000 0 500 1000 0 0 depth (m) 500 time (s) 0.5 1000

  31. Kirchhoff Migration Goal: Locate the singularities of δ c from δφ Requires F − 1 Recall: data are redundant Least Squares: F − 1 LS = (F ∗ F) − 1 F ∗ � � � F ∗ [ δφ ](x)= ω 2 B(x , r , s , θ )e − i ω (t − T(s , x) − T(x , r)) d θ dsdr R S R 2n − 1 (Beylkin (85), Rakesh (88), Symes (95))

  32. Velocity Analysis: Kirchhoff Methods Data are redundant, exploit the redundancy to find the velocity model c(x) �→ c(x , h), we find ( ∂ h F ∗ [c](d(s , r , t))) argmin c h can be • offset (almost NMO) • scattering angle • subsurface offset • time • . . . Symes’ 2009 review paper has an overview of this Symes 1999, 2001 justifies the use of local optimization for layered partially linearized inversions

  33. Imaging Methods – Summary • Kirchhoff ◮ Integral technique, usually uses ray theory ◮ Linearized with Kirchhoff approximation ◮ Related to X-ray CT imaging ◮ Generalized Radon Transform • For velocity analysis, iterate over ‘flatness’

  34. One-Way Methods Physical Motivation • downward continuation s r • imaging condition Claerbout 71, 85

  35. One-Way Methods Approximating the Wave Equation Idea (1D, c constant=1): ( ∂ 2 x − ∂ 2 t ) φ = ( ∂ x − ∂ t )( ∂ x + ∂ t ) φ c not constant: (c(x) 2 ∂ 2 x − ∂ 2 t ) φ = (c(x) ∂ x − ∂ t )(c(x) ∂ x + ∂ t ) φ − c(x)( ∂ x c(x)) ∂ x φ c(x) smooth ⇒ better approximation Taylor (81), Stolk & de Hoop (05) give more detail and more dimensions

  36. Imaging Methods – Summary • Kirchhoff ◮ Integral technique, usually uses ray theory ◮ Linearized with Kirchhoff approximation ◮ Related to X-ray CT imaging ◮ Generalized Radon Transform • One-way ◮ Based on a paraxial approximation ◮ Usually computed with finite differences • For velocity analysis, iterate over ‘flatness’

  37. Reverse-Time Migration Forming an Image Procedure: Whitmore (83), Loewenthal & Mufti (83), Baysal et al (83) • back propagate in time • imaging condition s r G 0 (x , t , s) G 0 (r , t , x) V(x)

  38. Reverse-Time Migration an Adjoint State Method Lailly (83,84), Tarantola (84,86,87) Symes (09) For a fixed source, s, � (c − 2 ∂ 2 t − ∇ 2 )q(x , t; s) = δφ (r , t; s) δ (x − r)dr R s q( · , t , · ) = 0 for t > T receivers act as sources, reversed in time � � 2 q(x , t; s) ∂ 2 Im(x) = t G 0 (x , t , s)dtds c 2 (x)

  39. Reverse-Time Migration Example Liu et al (07)

  40. Reverse-Time Migration Example Liu et al (07)

  41. Reverse-Time Migration Example Liu et al (07)

Recommend


More recommend