Seismic Inversion Chaiwoot Boonyasiriwat October 21, 2020 - PowerPoint PPT Presentation

Seismic Inversion Chaiwoot Boonyasiriwat October 21, 2020

Petroleum Exploration ▪ Petroleum was generated in organic-rich sedimentary rocks and then migrated upward and was trapped in various structures such as anticline and fault. ▪ Images of subsurface structures can be used to interpret where petroleum reservoirs are located in the earth. Image from http://media-3.web.britannica.com/eb-media/69/669-004-7B2A0E51.jpg 2

Seismic Exploration ▪ Mechanical waves generated by several sources propagate into the earth and reflect back to the surface. ▪ The wavefields recorded by receivers at the surface are the observed data used to estimate the earth structures. ▪ Seismic exploration can be performed on both land and marine environments. Subsurface Structure Seismic Data 3

Marine Seismic Acquisition Air gun source https://www.rigzone.com/training/insight.asp?insight_id=303&c_id= https://9ca7a3fb47-custmedia.vresp.com/ 303b3a0f74/Bolt_e_source_airgun.jpg Hydrophone http://filemaker2-server.cbl.umces.edu/ sensorimages/9663-TC4032.gif 4

Land Seismic Acquisition Vibroseis Geophone http://img.diytrade.com/cdimg/1296480/15832570/0/1286613918/LAND http://static.panoramio.com/photos/large/27525469.jpg _GEOPHONE_STRING.jpg 5 http://www.oilinuganda.org/wp-content/media/2013/06/Seismic-Illustration1.bmp

Traveltime Inversion A seismic velocity model v can be estimated by first- arrival travel time T by solving the constrained minimization problem subject to where T o is the observed data, x r and x s are the receiver and source locations, respectively, and s = 1/ v is slowness. 6

Traveltime Inversion Using the method of Lagrange multiplier, we first form the Lagrangian as where When , 7

Traveltime Inversion We have Setting yields the state problem Setting yields the adjoint problem 8

Traveltime Inversion Since , The gradient of the objective functional is the product between the slowness s and the Lagrange multiplier  . The state and adjoint problems can be solved using ▪ Ray tracing ▪ Fast marching method (Sethian and Popovici, 1999) ▪ Fast sweeping method (Zhao, 2004) ▪ Fast iterative method (Jeong and Whitaker, 2008) 9

Ray Tracing Let t be unit vector along ray path. We have where is slowness So or Rawlinson et al. (2008, p. 207) 10

Ray Equation Taking gradient of the last equation and using the definition of tangent vector t yields ray equation The ray equation is the Euler-Lagrange equation of the Lagrangian which is called the Fermat integral . Fermat’s principle of stationary time: ray paths correspond to extremal curves of Fermat integral. 11

Hamiltonian Formalism ▪ Ray paths are characteristic curves of the Hamiltonian which can be written in various form. ▪ In the Hamiltonian formalism of ray tracing, the momentum is corresponding to the slowness vector defined as . ▪ The eikonal equation can then be written as with a Hamiltonian 12

Characteristic System ▪ The 6-vector is called the canonical coordinate vector in a 6D phase space. ▪ is a hypersurface in the phase space. ▪ On the hypersurface, we have ▪ The Einstein summation convention is always used. ▪ Along a characteristic, we have Cerveny (2001, p. 104) 13

Characteristic System We then obtain the first 6 ODEs of the system The last equation for T can be obtained by which can be written in vector form as Cerveny (2001, p. 104) 14

Characteristic System ▪ This system of 7 first-order ODEs are the Hamilton equations and can be solved for the characteristic curve r , the slowness vector p , and the traveltime T . ▪ Note that the first 6 equations are independent of the last equation. ▪ The parameter u depends on the form of H . Cerveny (2001, p. 104) 15

General Ray Tracing Systems ▪ Starting from a general Hamiltonian ▪ The corresponding characteristic system is Cerveny (2001, p. 105) 16

Examples of Ray Tracing Systems ▪ When , ▪ When , ▪ When , Cerveny (2001, p. 105) 17

Examples of Ray Tracing Systems Because , the 3 systems become ▪ When , ▪ When , ▪ When , Cerveny (2001, p. 106) 18

Shooting Method: Initial-Value Problem Let’s use the system with Initial conditions Initial point is source point In 2D, the system becomes Initial momentum is slowness vector where is the shooting angle 19

Hamilton-Jacobi Equation General form of Hamilton-Jacobi equation (HJE) where is the generalized coordinates is the Hamilton’s principal function is the Hamiltonian is the generalized momentum Hamilton-Jacobi equation is first-order nonlinear PDE while Euler-Lagrange equations are second-order PDE. It can be shown that and, therefore, is the classical action. https://en.wikipedia.org/wiki/Hamilton-Jacobi_equation 20

Static Hamilton-Jacobi Equation Static Hamilton-Jacobi equation (HJE) can be written as So the eikonal equation is a special case of HJE when Gremaud and Kuster (2006) 21

Numerical Hamiltonian On a 2D Cartesian grid, the Hamilton-Jacobi equation can be discretized as Let Godunov-flux Hamiltonian (first-order upwind): Lax-Friedrichs-flux Hamiltonian Gremaud and Kuster (2006, p. 2) 22

2D Problem We consider the 2D boundary-value problem “The slowness field can be unbounded, corresponding to the presence of obstacles.” (Gremaud and Kuster, 2006) 23

Upwind FD Scheme In 2D, the upwind scheme can be written as where In practice, we use a square grid with Sethian and Popovici (1999) 24

Godunov Scheme Godunov upwind scheme can be written as where Zhao (2004) 25

Update Cases Let Traveltime at a node can be updated in one of these cases: Case 1: Case 2: Case 3: Case 4: 26

Fast Sweeping Method ▪ Apply the Godunov upwind scheme to all grid points except the boundary points and replace the old value if the new value is smaller. ▪ 4 sweeping schemes are used: Zhao (2004) 27

Frequency-domain FWI ▪ Inversion of first-arrival traveltime data can only provide a smoothed velocity model. ▪ Inversion of full-waveform data can provide a higher- resolution velocity model. ▪ Suppose the wave propagation is governed by ▪ Inserting a planewave into the wave equation yields the Helmholtz equation where s is slowness. 28

2D Helmholtz Problem In 2D, wavefield is governed by The direct scattering problem can be solved using the finite-difference frequency-domain (FDFD) method:

2D Helmholtz Problem Absorbing boundary conditions are used at the boundary For example, the absorbing boundary condition at the bottom is x z Using the finite difference method, we obtain

Example Result Ajo-Franklin (2005, p. 13)

Frequency-domain FWI ▪ Consider the constrained minimization problem subject to ▪ The corresponding Lagrangian is where 32

Frequency-domain FWI The gradient of objective functional with respect to slowness is when and in this case 33

Frequency-domain FWI Setting yields the state equation Since setting yields the adjoint equation 34

Example Result: Wave Path Ajo-Franklin (2005, p. 14)

References ▪ Ajo-Franklin, J. B., 2005, Frequency-domain modeling techniques for the scalar wave equation: An introduction, Earth Resources Laboratory Technical Report, MIT. ▪ Cerveny, V., 2001, Seismic Ray Theory, Cambridge University Press. ▪ Gremaud, P. A., and C. M. Kuster, 2006, Computational study of fast methods for the eikonal equation: SIAM Journal on Scientific Computing, 27, 6, 1803-1816. ▪ Jeong, W., and R. T. Whitaker, 2008, A fast iterative method for eikonal equations: SIAM Journal of Scientific Computing, 30(5), 2512 – 2534. ▪ Rawlinson, N., J. Hauser, and M. Sambridge, 2008, Seismic ray tracing and wavefront tracking in laterally heterogeneous media: Advances in Geophysics, 49, 203-273. ▪ Sethian, J. A., and A. M. Popovici, 1999, 3-D traveltime computation using the fast marching method: Geophysics 64(2). ▪ Zhao, H., 2004, A fast sweeping method for eikonal equations: Mathematics of Computation, 74, 603-627. 36