' &

Numerical Methods for Traveltimes and Amplitudes

William W. Symes

The Rice Inversion Project Department of Computational and Applied Mathematics Rice University Mathematical Geophysics Summer School Stanford University, August 1999
1

'
Outline
1. Hamilton-Jacobi Theory: wavefronts and rays 2. Lagrangian computation of traveltimes: Ray Tracing 2.1 Fundamental concepts and notations 2.2 Rays as solutions of initial value problems: Ray Shooting

&

2.3 Wavefront Construction 2.4 Rays as solutions of 2 point boundary value problems: Ray Bending

'
Outline contd
3. Eulerian computation of traveltimes: Eikonal Solvers 3.1 Viscosity solutions: rst arrival times 3.1 Space Marching methods: an example 3.2 Overview of algorithms 4. Reections

&

5. Amplitudes 6. Eulerian approaches to computation of multiple arrival times

'
Hamilton-Jacobi Theory: Wavefronts and Rays

&

'
Wavefronts and orthogonal trajectories Recall eikonal equation for wave phase or traveltime T (x), x R3 in velocity eld v (x): v 2 T T = 1 Wavefronts are the surfaces T = const.

&

Orthogonal trajectories are curves perpindicular at each point to wavefront, i.e. with velocity vector parallel to T
5

' &

Rays and Hamiltons Equations

Physical units natural ODE for orthogonal trajectories is x = v 2 (x)T (x) t Introduce ray slowness vector p(t) T (x(t)), so that (x(t), p(t)) forms a trajectory in phase space R6 . Claim: provided that T is of class C 2 in its domain of denition, the phase space trajectory (ray) (x(t), p(t)) satises Hamiltons Equations: dx dp = p H, = x H dt dt 2 v (x)p p in which H is the Hamiltonian H (x, p) = 1 2
6

' &

Proof of Claim
2

1 p H = v p, x H = (v 2 )p p 2 so rst equation is just ODE for orthogonal trajectories introduced above. From the denition of p, dp dx = T (x) = v 2 T (x) T (x) dt dt 1 2 1 2 = v (T (x) T (x)) = ( v )T (x) T (x) 2 2 because of the eikonal equation; the last expression is just the second Hamilton equation in view of the denition of p. q.e.d.
7

'
Observations on Hamiltons Equations Observation 1. For any solution (x(t), p(t)) of Hamiltons equations, H (x(t), p(t)) is independent of t (HW Problem 1!) Observation 2. If v 2 (x0 )p0 p0 = 1, then the solution (x(t), p(t)) of Hamiltons equations with initial data (x0 , p0 ) satises dx dx v 2 (x) dt dt (HW Problem 2!)

&

' &

The Method of Characteristics 1

Hamiltons equations are a priori independent of any solution of the eikonal equation, but may be used to construct solutions, provided that v has Lipshitz continuous partial derivatives. Ingredients: S R3 bounded smooth surface function T0 : S R = Cauchy data for eikonal equation; must satisfy v 2 S T0 S T0 < 1 (S = gradient along S ) projection x0 from R3 onto tangent space to S at x0 S
9

'
The Method of Characteristics 2 System of 2 linear and 1 quadratic equation: x0 p = S T0 (x0 ); v 2 (x0 )p p = 1 has two solutions at each x0 S ; make consistent (smooth) choice p0 : S R3 .

&

At each x0 S solve Hamiltons equations with initial data (x0 , p0 (x0 )). Easy excercise in Implicit Function Theorem: (x0 , t) x(t) is an invertible, dierentiable map for < t < , x0 S provided is small enough.
10

' &

The Method of Characteristics 3

$
dx . dt

That is, for each point x close enough to S , there is a unique choice of t, x0 so that the ray with initial data (x0 , p0 (x0 )) passes through x at time t. Dene T (x) = t. Claim: T , dened near S as above, satises the eikonal equation.
x Proof: By construction, T || d , so dt

dx dx dx d (t) = c (t) (t) 1 = T (x(t)) = T (x(t)) dt dt dt dt

where c is the ratio of the lengths of T (x(t)) and

11

'
Method of Characteristics 4 From Observation 2, this is = cv 2 (x(t)) c = v 2 (x(t)) dx dx T (x) T (x) = v (x) = v 2 (x) dt dt
4

q.e.d

&

HW Problem 3: show that the solution T (x) constructed by the method of characteristics satises p(t) = T (x(t)) as expected.
12

'
Observations on the M. of C. provides a local construction of solutions to the eikonal from Cauchy data on a smooth surface requires that v have Lipschitz derivatives (so must be modied for less regular v , eg. with discontinuities)

&

fails where rays cross - because then denition of T becomes ambiguous. Apparently T may be multivalued in general.

13

' &

Another form of the Ray Equations dx dp 2 = v (x)p, = v (x) v (x) p p dt dt

Reserve the term ray for solutions (x(t), p(t)) of this system satisfying the addtional condition p p = v 2 (x) only these solutions have expected relation with eikonal (p = T (x)) For rays, second equation is equivalent to dp = log v (x) dt

14

'
A second order ray equation Dierentiate the rst ray equation and eliminate second to get dx d2 x = 2 v (x) 2 dt dt dx v (x)v (x) dt
dp dt

$
using

&

First term: acceleration parallel to instantaneous direction of ray. Second term: acceleration perpindicular to ray

15

'
Simple consequences ray equations rays bend toward decreasing velocity If v = const (indepedent of x), then p = const along each ray, i.e. ray doesnt change direction, i.e. ray is straight line traversed at constant speed.

&

If v is independent of some space coordinate, say x, then corresponding component of p is constant along ray.

16

'
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.1 0.2

Example: Constant Velocity

10 0 2 9 0.2 1.8 8 0.4 1.6 7 0.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.8

1.2

1.4

1.6

1.8

0.3

0.4

0.5

0.6

0.7

0.8

0.9

&

Figure 1: Left: rays from xs =0.55, 101 equally spaced angles; Right: interpolation from irregularly sampled rays onto regularly sampled grid - time defualted to 0 in region of no coverage

17

' &

Yet another equivalent form...

In terms of Slowness eld u(x) = 1/v (x), arc length along ray s: dx d u(x) = u(x) ds ds Derivation from Hamiltonian form: HW Problem 4 Follows from this: if x(s) is space component of ray, parametrized by arc length, and xs = x(s0 ), xr = x(s1 ), then s1 T (xr ) T (xs ) = u(x(s)) ds
s0

(derivation: HW Problem 5)
18

' &

Regularity and Fermat

Standard Thm on local existence, uniqueness of rays: v must have Lipshitz continuous derivatives. Extension possible to less regular velocity elds (eg. piecewise smooth) via relation to Fermats Principle: (space component of) ray makes the traveltime integral u ds stationary amongst all C 0 paths connecting xs to xr . Note relaxed regularity of path in statement of Fermat: 1 ) can also relax regularity of v (eg. to piecewise CLip
19

' &

Still another form

2D case (3D similar): x = (z, x)

1 on any ray, can write Because H = 2 p = v 1 (x)(cos , sin ) and describe ray in terms of z, x, : dz = v (z, x) cos dt dx = v (z, x) sin dt d v v = sin cos dt z x This form of ray equations particularly convenient for computation.
20

'
Traveltime is:

Two Points of View

a quantity which evolves along rays a solution to the eikonal equation Corresponding families of numerical methods:

&

Lagrangian: solve the ray equations numerically (ray tracing) Eulerian: solve the eikonal equation numerically (eikonal solvers)
21

'
Lagrangian computation of traveltimes: Ray Tracing = Method of Characteristics

&

22

'
Given: Find:

Problem Statement

source point xs , receiver points xr R velocity eld v (x) traveltime(s) from xs to each xr , or equivalently a ray (x, p) and a time t for which x(0) = xs , x(t) = xr then t is (a) traveltime from xs to xr

&

variant # 1: nd the smallest such traveltime variant # 2: nd all such times

23

' &

Lagrangian Method 1: Shooting

Idea: sample a circle (sphere in 3D) of p of radius v 1 (xs ); for each sample ps , approximate a ray by solving the initial value problem dp dx 2 = v (x)p, = log v (x) dt dt x(0) = xs , p(0) = ps numerically, using either general- or special-purpose solver. Must interpolate times on rays to output points somehow.
24

' &

Lagrangian Method 1: Shooting

Example implementation: ray.m, a Matlab ray tracer Principle: apply general purpose ODE/IVP solver to ray equations Inputs: source and receiver coordinates, a function vel.m returning the velocity and its partial derivatives at any point in space (gridded velocity interpolation)

Outputs: graphical representation of rays, interpolated times at receiver points ODE solver: ode113 (Matlab variable order nonsti solver - could use others); interpolation is linear, using griddata
25

'
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.1 0.2

Example: Constant Velocity

10 0 2 9 0.2 1.8 8 0.4 1.6 7 0.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.8

1.2

1.4

1.6

1.8

0.3

0.4

0.5

0.6

0.7

0.8

0.9

&

Figure 2: Left: rays from xs =0.55, 101 equally spaced angles; Right: interpolation from irregularly sampled rays onto regularly sampled grid - time arbitrary in region of no coverage

26

'
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.1 0.2 0.3

Example: Linear Velocity

9 0 0.8

$
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.8 0.9 1

0.2 8 0.4 7 0.6 6

0.8

1.2

1.4 3 1.6 2 1.8 1 0.4 0.5 0.6 0.7 0.8 0.9 1

&

2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

Figure 3: v = 1 + 4z . Left: rays from xs =0.55, 101 equally spaced angles; Right: interpolation from irregularly sampled rays onto regularly sampled grid - time arbitrary in region of no coverage

27

' &

Special ray equation solvers

Many based on the observation: in linear velocity, rays are arcs of circles (Proof: HW Problem 6) Favorite model parametrization: v(x) piecewise linear. Fermat C 0 , piecewise C 2 generalized rays. exact (generalized) solution of ray equations in piecewise linear velocity eld: link circular arcs by Snells law at boundaries (eg. Chander, 1975) much more ecient than general purpose IVP solver for these special models Variants: approximate circular arcs by quadratics in t, lose transcendentals (eg. Langan et al. 1985)
28

' &

Eects of Ray Bending

Recall: rays bend towards slower regions example lens1: a fast region in center - rays bend away example lens2: a slow region in center - rays bend toward center, begin to cross example lens3: an even slower region - rays bent more strongly 0 - 3 rays pass over each pt interpolation very confused

example lens4: several strong lenses - complex ray pattern, interpolation is hopeless.
29

'
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.1 0.2 0.3 0.4

Example: A Fast Lens

1.2 0 1.8 1.18 0.2 1.6 1.16 0.4 1.4 1.14 0.6 1.2 1.12 0.8 1 1.1 1 0.8

$
0.6 0.4 0.2 0 0.7 0.8 0.9 1

1.08

1.2

1.06

1.4

1.04

1.6

1.02

1.8

1 0.5 0.6 0.7 0.8 0.9 1

2 0 0.1 0.2 0.3 0.4 0.5 0.6

&

Figure 4: v = 1 + d(fx sin x + fz cos z ), d = 0.2, fx = 1, fz = 0.5. Left: rays from xs =0.55, 101 equally spaced angles; Right: interpolation from irregularly sampled rays onto regularly sampled grid
30

'
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.1 0.2 0.3 0.4

Example: A Slow Lens

1 0 0.2 2 0.96 0.4

$
1.5 1 0.5 0 0.7 0.8 0.9 1

0.98

0.94

0.6

0.92

0.8

0.9

0.88

1.2

0.86

1.4

0.84

1.6

0.82

1.8

0.8 0.5 0.6 0.7 0.8 0.9 1

2 0 0.1 0.2 0.3 0.4 0.5 0.6

&

Figure 5: v = 1 + d(fx sin x + fz cos z ), d = 0.2, fx = 1, fz = 0.5. Left: rays from xs =0.55, 101 equally spaced angles; Right: interpolation from irregularly sampled rays onto regularly sampled grid
31

'
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.1 0.2 0.3

Example: A Slower Lens

1 0 0.2 0.95 2.5 0.4 0.9 0.6 0.85 2

$
1.5 1 0.5 0 0.8 0.9 1

0.8

0.8

0.75

1.2

1.4 0.7 1.6 0.65 1.8 0.6 0.4 0.5 0.6 0.7 0.8 0.9 1

2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

&

Figure 6: v = 1 + d(fx sin x + fz cos z ), d = 0.4, fx = 1, fz = 0.5. Left: rays from xs =0.55, 101 equally spaced angles; Right: interpolation from irregularly sampled rays onto regularly sampled grid
32

'
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.1 0.2 0.3

Example: A Complex Lens

1.2 0

$
2 1.5 1 0.5 0 0.8 0.9 1

0.2 1.15 0.4 1.1 0.6 1.05

0.8

0.95

1.2

1.4 0.9 1.6 0.85 1.8 0.8 0.4 0.5 0.6 0.7 0.8 0.9 1

&

2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

Figure 7: v = 1 + d(fx sin x + fz cos z ), d = 0.4, fx = 3, fz = 1. Left: rays from xs =0.55, 101 equally spaced angles; Right: interpolation from irregularly sampled rays onto regularly sampled grid

33

'
Ray Shooting Issues: Shadow Zones no rays no data for interpolation possibile solution: shoot more rays appropriate selection of rays not evident a priori, must be discovered dynamically.

&

This is a limitation of ray shooting.

34

' &

Ray Shooting Issues: Multiple Arrivals

generally, several rays pass over each point at sucient distance from source traveltime is multivalued Cannot interpolate to give a sensible single-valued output a solution (variant # 1): discard all arrivals but one another solution (variant # 2): interpolate a multivalued function (!). This a feature of the traveltime eld, not of any particular computational method.
35

' &

Wave kinematics are multivalued!

Wavefield at 350ms with 325ms traveltime contour 400

450

500

550
z (m)

600

650

700

750 200

250

300

350

400 x (m)

450

500

550

600

Figure 8: Accurate nite dierence solution of wave equation with low velocity lens, showing presence of multiple arrivals
36

' &

Lagrangian Method 2: Wavefront Construction

Vinje et al (1993): advance all rays of interest by uniform steps in time (thus implicitly advancing wavefront = constant time surface). When rays grow too far apart, interpolate new rays. When rays converge, delete rays. Thus maintain more or less uniform ray density. some use in exploration geophysics can either keep least-time ray (variant # 1) or all rays (variant # 2) complex, dynamic data structure relatively expensive
37

'
Lagrangian Method 3: Two Point Ray Tracing Idea: to nd time for wave to travel from xs to xr , solve the BVP dp dx 2 = v (x)p, = log v (x) dt dt

$
x(0) = xs ; p(0) p(0) = v 2 (xs ); x(tr ) = xr ;

&

Unknowns are trajectory x(t), p(t), and arrival time tr .

38

' &

Lagrangian Method 3: Two Point Ray Tracing

Three general approaches: A. Shooting: regard x(tr ) = Fshooting [x(0), p(0), tr ] as function of x(0), p(0), tr (with constraint that p(0) p(0) = v 2 (ps ),; solve the system Fshooting [x(0), p(0), tr ] = xr B. Bending: regard the BVP itself as a large nonlinear system, solve somehow

C. Continuation: solve BVP in simpler velocity eld (eg. v0 = const.), evolve along a 1 param. family of models v0 v . (Keller & Perozzi, 1983)
39

' &

Lagrangian Method 3B: Ray Bending

Solve discretized BVP as nonlinear system for x(t), p(t), and tr [or equivalent system]. Two general approaches: ad hoc discretization, system solver (eg. Newton with constant damping) (Julian & Gubbins 1977, Farra 1992,...) standard NA technology for 2 pt bdry problems: variable order, variable step discretization, quasiNewton method with globalization (linesearch or trust region) (Pereyra et al 1980).
40

' &

Ray Bending in PW Constant Velocity

Perhaps the most popular algorithm of all! Velocity, slowness = piecewise constant on rectangular grid Rays = straight lines inside cells, bend according to Snell at cell boundaries satisfy Fermat +: Very easy to implement - leads to ray path matrix model for linearized tomography M s = t, cf. Berryman lectures -: Leads to nonconvergent tomographic results under grid renement, see Delprat-Jannaud and Lailly 1993.
41

' &

Multiple Arrivals impede Ray Bending

0 2 0.2 1.9

0.4

1.8

0.6

1.7

0.8

1.6

1.5

1.2

1.4

1.4

1.3

1.6

1.2

1.8

1.1

2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 9: Rays in piecewise linear velocity: v(z, x) = 2 2x if x < 0.5, v(z, x) = 1 else. For any point (zr , xr ) in x > 0.5, beginning bending with the straight line from xs guarantees that the curved rays through the same point will never be found! (see Farra 1992, Snieder & Spencer 1993)
42

'
roughly,

Comparison

All 2 pt ray tracing methods approximate solution by (i) discretization, (ii) iteration (Newton,...) relatively expensive compared to ray shooting + interpolation, maybe more accurate

&

Shooting + interpolation for large number of output pts with modest accuracy, 2 pt ray tracing for small number of output pts with low high accuracy NB: many caveats - choice must be problem driven
43

'
Eulerian computation of traveltimes: Eikonal Solvers

&

44

'
Eikonal Equation 1 |T (x)| = u(x) = v (x) Requires side conditions to determine solution uniquely: give T (x) along smooth hypersurface

$
|T (x) T (xs ) u(xs ) x xs | = o( x xs )

&

point source (limit of small spheres):

45

' &

Local vs. global solutions

For smooth u(x), Method of Characteristics (= ray tracing) gives smooth solution near datum set (hypersurface, source point) Problem: like many nonlinear PDEs, global smooth solutions do not exist: method of characteristics breaks down where characteristics = rays cross Another gloss: implicit in statement of eikonal equation is single-valuedness of T (x). However we already know that traveltime is generally multivalued at sucient distance from point source or other initial value surface.
46

' &

Viscosity solutions

Since smooth solutions do not exist (globally), must generalize concept of solution. Viscosity solutions = C 0 , piecewise C 1 fcns satisfying eikonal equation whereever possible + addl conditions. Existence, uniqueness for v C 0 : Lions 1982 (defns, proofs). Related idea: entropy solutions of conservation laws Lions 1982: The viscosity solution of the point source problem is the rst arrival time eld. gridded eikonal solvers compute rst arrival time.
47

'
Other Space Marching Methods Large literature - see reference list. Very small selection: unidirectional Godunov schemes Vidales scheme (1988)

&

Sethians Fast Marching method (1996)

48

'
Unidirectional Space Marching Idea: view eikonal equation as evolution equation in a space variable, at least locally (perhaps dierent space variables in dierent places): for example, if you choose +z , get

$
T = z u2 T x
2

&

(2D version - 3D version similar).

49

'
Simple Dierence Scheme Aim: produce approximations to traveltime T (x) on a regular grid (mx, nz ):
n T (mx, nz ) Tm

&

Principle: use connection between rays, eikonal to extrapolate times from depth level n to depth level n + 1 by approximate backwards ray tracing

50

'
Geometry of Simple Scheme
m-1 x z approximate ray m x m+1 n

$
n+1

&

51

' &

Approximations

Locally, approximate slowness u as constant = un m = u(mx, nz ) time on level n as piecewise linear rays as straight lines Ray to the target point (mx, (n + 1)z ) intersects level n at ( x, nz )

n+1 = time at intersection point ( x, nz ) + time along Tm ray segment

52

' &

Time at intersection point

$
m

Three possibilities: R: ray comes from right: mx < x < (m + 1)x L: ray comes from left: (m 1)x < x < mx O: ray comes from outside [(m 1)x, (m + 1)x] Assume rst - second similar, third handled at end = T n + ( x mx)D+ T n , where Then T ( x, nz ) T
m x n n T T m m 1 n Tm = Dx , x

53

'
Time along Ray Segment Locally const slowness approximation time increment along ray segment = slowness length
2 + z 2 ( x m x ) = un m

Thus

&

+ n n n+1 Tm + un + ( x mx)Dx = Tm Tm x mx)2 + z 2 m (

54

'
Key ingredients:

Determining x

1. value of Hamiltonian: p p = u2 = v2p 2. rst ray equation x 3. basic Hamilton-Jacobi equation: T = p

&

Write p = (q, p). 2 1. q = u p2 2. slope of ray =

x z

p q

=
55

p u2 p2

'
Therefore

Determining x

$
p z p2 +1 2 2 n (um ) p

x x mx = z = z Get two conclusions from 3:

2 2 (un m) p

+ n Tm . Since (1) From 1. get p = Dx q = u2 z > 0, p = u2 x < 0 for ray as in diagram, + n Tm < 0 assumption (R) above is possible only if Dx

&
(2)

n+1 Tm

n Tm

p2 Z
2 2 (un m) p
56

un m z

'
Synthesis
+ n Tm , A little bit of algebra (HW Problem 7), plus p = Dx gives n n+1 2 + n 2 + z (un = Tm Tm m ) (Dx Tm ) + n Tm < 0. in case (R) i.e. when Dx n Tm > 0) Similarly in case (L) (Dx

$
n n+1 2 n 2 + z (un = Tm Tm m ) (Dx Tm )

&

57

' &

Choices

Note that one, both, or neither of (R) or (L) could occur. if exactly one occurs, then preceding slide gives n+1 formula for Tm
n+1 from if both occur, then select smaller value of Tm above two formulae this scheme will compute the rst arrival time = viscosity solution! + n n Tm , then Tm 0 Dx if neither occur, because Dx rays are diverging; compute time increment using vertical ray, n n+1 + zun = Tm Tm m

58

'
Properties:

Convenient form

HW Problem 8: show that the above four cases can be subsumed in a single update formula:
n n+1 + = Tm Tm + n 2 2 n 2 z (un m ) max((max(0, Dx Tm )) , (min(0, Dx Tm )) )

= Godunov scheme (eg. Osher and Sethian 1988)

&

formally rst order accurate for z -marching eikonal upwind: dierence in direction from which ray comes
59

' &
T = z

Practicalities

$
2

If this occurs, algorithm becomes unstable! Example of Courant-Friedrichs-Lewy instability of dierence schemes. Remedy: legislate such steep rays out of existence! Gray and May 1994: modify eikonal to paraxial eikonal max u2

Still dangling: case (O), in which rays are converging but from outside the interval [(m 1)x, (m + 1)x].

T x

Correctly computes rst arrivals along rays making angle < with z -axis.
60

, u2 cos2

'
Matlab 1st Order Godunov Program eik.m Input parameters: function returning velocity eld and partial derivatives as in ray.m, location of source, computational rectangle, x, aperture limit - z computed internally to satisfy CFL stability criterion for Godunov scheme, paraxial eikonal Output: traveltime eld throughout computational rectangle, displayed in color contour plot.

&

61

' &

Examples

A slower lens (lens3): note that eikonal solver returns much smoother traveltime eld estimate than shooting + interpolation. Eikonal gives rst arrival time, lower part returned by shooting is not any traveltime. Note transition locus from left-arriving rays to right-arriving rays: this is analogous to shock in uid mechanics. A piecewise linear velocity (badbend): eikonal solver returns complete eld (though incorrect for points outside aperture because of paraxial modication). No shadow zone!
62

'
0 0.2 0.4 0.6 0.8

Example: a slower lens

2.5

$
2.5 2 1.5 1 0.5 0 0.7 0.8 0.9 1

0.2
2

0.4

0.6
1.5

0.8

1
1

1.2

1.2
1.4

1.4
1.6 0.5

1.6
1.8

1.8
2 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

&

2 0 0.1 0.2 0.3 0.4 0.5 0.6

Figure 10: v = 1+ d(fx sin x + fz cos z ), d = 0.4, fx = 1, fz = 0.5. Left: eikonal solver; Right: shooting + interpolation

63

'
0 0.2 0.4 0.6 0.8

Example: piecewise linear velocity

2.5

$
2.5 2 1.5 1 0.5 0

0.2

0.4

0.6

1.5

0.8

1.2

z
0.5 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

1.2

1.4

1.4

1.6

1.6

1.8

1.8

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

&

Figure 11: piecewise linear velocity: v(z, x) = 2 2x if x < 0.5, v(z, x) = 1 else. Left: eikonal solver; Right: shooting + interpolation. Note that times are similar where only one ray passes; also absence of shadow zone in eikonal result.
64

' &

Extensions, performance

Essentially NonOscillatory order > 1 extensions of Godunov & similar schemes: Osher and Sethian 1988, Osher and Shu 1991 much more ecient than 1st order, needed for amplitude computation cost = O(N ), N = nx ny ... = number of gridpoints

limited by paraxial assumption: correctly computes times along rst-arriving rays making angle < with vertical - OK for some apps, not for others
65

'
Vidales Scheme Vidale, 1988: rst generally reliable solver, got the eld going. Features: implicit upwind scheme

&

space marching in all directions - 4 in 2D, 6 in 3D second order

66

'
Vidales scheme Concept: given times T on a square array, extrapolate to neighboring rows and columns of grid cells. [Array is initially a point.] nd minima of T on sides of box

&

rst: at minima, gradient of T should be normal to box so extrapolate T linearly: if grid step is , then Tnext = Tmin + u

67

'
2

Vidales scheme

$
2

in other grid cells with corners A, B, C, and D, solve centered dierence approximation on diagonal grid
2 2 2 u2 A + uB + uC + uD (TA TC ) +(TB TD ) = (x +z ) 4 2 2

to compute T at 4th corner given T at other three

&

start at cells next to mins where 3 values are known, work out from mins - in eect, with proper ordering system is block triangular

NB: implicit and upwind

68

'
...

Local Geometry of Vidale scheme

$
n n+1

Tmin
box box box box box

N+2

N+1

linear extrapolation

...

&

Figure 12: Extrapolation along normal at minimum time point performed rst (large circle large box), then box scheme used to update other points in order indicated
69

'
B

Box scheme
A

$
D

&

Figure 13: Given TZ , TB , TC , for example, nd TD by solving (TA TC )2 + (TB TD )2 = (x2 + z 2 )uavg
70

'
Advantages: Disadvantages:

Vidales scheme

propagates times in all directions, extends to 3D good accuracy; cost = O(N ), N = nx ny ... = number of gridpoints

&

stops (system inconsistent) when ray carrying rst arrival is tangent to side of expanding box requires nding minima on sides
71

' &

Sethians Fast Marching Method

Sethian 1996, Sethian & Popovici 1999, Sethian 1999 for many references, seismic & nonseismic applications, and demos: www.math.berkeley.edu/~sethian Based on upwind discrete eikonal:
+ n n Tm , 0))2 + Tm , Dx (max(Dx 2 ) = (un m

+ n n Tm , 0))2 Tm , Dz (max(Dz

72

' &

Fast Marching - Initialization

Presume T given on initial subgrid ( initial value surface - may be single point) - forms initial set of accepted points Assign T = at all other inactive points NB: At a gridpoint next to one or more with nite T values, only dierences with known values contribute to LHS! Solve discrete eikonal to assign values to all neighbors of initial (accepted) points: these form initial set of active points (no longer inactive)
73

' &
Thats it!

Fast Marching - Propagation

Repeat until all points are accepted: 1. Select active point with smallest value of T , move it to accepted set 2. Solve discrete eikonal to assign T values to all neighbors of accepted points, move these to active set (actually needs to be done only with neighbors of point selected in 1.)

74

'
* * * * * * * * * * * * * * * * * * * * * *

&

Figure 14: Fast Marching Preinitialization: initial data forms initial accepted point set = blue squares, all other points are inactive, assigned value T = = asterisks

75

'
* * * * * * * * * * * * * * * * *

&

Figure 15: Fast Marching Initialization: assign tentative T values to all neighbors of accepted points using discrete eikonal; these form initial active point set = red circles

76

'
* * * * * * * * * * * * * * * * *

&

Figure 16: Fast Marching Propagation, step 1A: select active point having smallest value of T = green circle

77

'
* * * * * * * * * * * * * * * * *

&

Figure 17: Fast Marching Propagation, step 1B: move active point with smallest time to accepted point set

78

'
* * * * * * * * * * * * * * * *

&

Figure 18: Fast Marching Propagation, step 2: identify all neighbors of accepted points, assign T values using discrete eikonal, move to active set

79

'
Fast Marching - Implementation sorting step 1. done once per gridpoint - with heap sort, cost is logarithmic total cost = O(N log N ), N = nx ny ... total number of gridpoints

&

higher order variants more ecient for some seismic applications - see Sethian 1999

80

'
Reections

&

81

'
Reections Weak form of wave equation + asymptotic analysis (or physics) leads to Snells law for reected waves: the time of arrival for a wave reected from an interface S is the sum of times along incident and reected rays. These rays meet at a point on the surface, making equal angles with its normal.

&

Incident ray (xd , pd ) connects source xs to x S in time td ; reected ray (xu , pu ) connects x to receiver xr in time tu ; time of arrival t = tu + td .

82

' &

Reections

Equal angles law: bisector of pd , pu parallel to surface normal S . Since pd (td ), pu (tu ) have same length (eikonal equation!), (pd (td ) + pu (tu )) S (x ) = 0 This is a system of equations! Either trace candidate incident and reected rays, or use eikonal solver to get pu = Tu etc.

In either case nd reection point x etc. by solving equal angles condition.

83

' &

Source

Receiver

incident ray s

reflected ray

Reflection point

Figure 19: Trace rays (or compute traveltimes) from source and receiver to points on the surface, and nd those where the equal angles condition holds (there may be more than one - is there here?)

84

'
Amplitudes

&

85

' &

Amplitudes 1

Recall from rst lecture: to leading order in frequency, wave amplitude A satises the transport equation 2T A + (2 T )A = 0 [NB: this A is the exponential of the A in Berrymans lecture.] Given known T , this is an advection equation for A. Just as for eikonal, two choices: method of characteristics nite dierence method on a grid
86

' &

Amplitudes 2

method of characteristics is trivial for linear scalar equations: along a ray, dA = v 2 (2 T )A dt solve by quadrature gridded solver: use equivalent form (A2 T ) = 0 plus divergence theorem (Courant & Hilbert 1962 Ch. 6, Friedlander 1958 Ch.1)
87

' &

Amplitudes 3

$
A2 v

The divergence theorem gives 0=

R

(A2 T )

for any volume R Rn , (n = 2 or 3). Take for R a ray tube (see gure); then 0= A2 v

D1

D0

Shrinking the caps D0 and D1 about points x0 and x1 on the axial ray, in limit get
88

'
source pt

$
D1 R
axial ray

D2

&

Figure 20: Ray tube construction: caps Di are perp. to ray family, so normals on caps are parallel to rays; sides consist of rays (dashed lines) , so normals on sides are perp. to rays.

89

89-1

' &

Amplitudes 4 A2 (x0 ) A2 (x1 ) = v (x1 ) v (x0 )

where is the ratio of area element on D0 to area element on D1 (transverse Jacobian) on the two caps at the axial ray. v (x1 ) A(x1 ) = A(x0 ) v (x0 )

Therefore to compute amplitudes need only compute the area element .

90

'
Amplitudes 5 As in gure, make D0 part of small circle (2D!), assume velocity constant between there and source point. Coordinate on circle is = ray takeo angle (at source). Since caps are perp. to the unit vector v 1 dx/dt, (t, ) (z, x) 1 =v =v (t, ) (z, x) = v (T ) 1
1

&

91

'
Amplitudes 6 Upshot: for point source in 2D (similar formula in 3D), v A = T 2 2 Since take-o angle is constant along rays,

$
T = 0

&

which is simpler than transport equation for A.

92

' &

Amplitudes 7

Some important observations: when rays converge (caustic), area element ratio blows up this simple version of geometric optics is limited to regions of smooth single arrivals - else must use global G. O. (Maslov Theory - see Duistermaat 1973). Numerically: must have at least 1st order convergent must use at least 3rd order scheme for T !

Also: since T is singular at a point source, in order to obtain (eective) high order convergence may need adaptive gridding.
93