Rytov/Eikonal Approximation of Wavepaths

William S. Harlan

Date: April 1998

Summary

This short paper derives a simple method of constructing wavepaths (bandlimited raypaths) from traveltimes using the Eikonal approximation. The Rytov approximation linearly perturbs the phase of a wavefield with respect to model parameters such as slowness. (The Born approximation perturbs the amplitude.) When the Green's functions for point sources are replaced by Eikonal approximations, the Rytov perturbed wavefield becomes a scaled, differentiated, time-delayed version of the reference wavefield.

Preliminaries

I will illustrate these approximations with a scalar wave equation for isotropic pressure, assuming constant density. You can generalize the results for more general elasticity if you are willing to allow separation of modes.

Let us first choose our notation.

$$v({\svector x})^{-2} \ddot{p} ( {\svector x}, t ) - \nabla^2p ( {\svector x}, t ) = w ( {\svector x}, t );$$ (1)

$${\rm or } ~~ 4 \pi^2 f^2 s({\svector x})^2 p ( {\svector x}, f )+ \nabla^2p ( {\svector x}, f ) = - w ( {\svector x}, f ).$$ (2)

The pressure $p$ is the non-zero diagonal of the strain tensor. The vector ${\svector x}$ is the spatial coordinate. The reciprocal of velocity $v$ is the slowness $s$. The source term $w$ is the divergence of body forces. $f$ is the Fourier frequency, using Bracewell's conventions with $- i 2 \pi$ in the exponent for the forward transform. Transformed and untransformed functions will be distinguished by arguments.

A Green's function $G$ solves the following equation.

$$4 \pi^2 f^2 s({\svector x})^2 G ({\svector x}, {\svector x}' ,f) + \nabla^2G ({\svector x}, {\svector x}' ,f)= - \delta ( {\svector x}- {\svector x}' ),$$ (3)

$${\rm so ~ that } ~~ p ( {\svector x}, f ) = \iiint G ({\svector x}, {\svector x}' ,f)w ( {\svector x}' ,f ) d {\svector x}' .$$ (4)

If we perturb the Green's function definition (3) and discard second-order terms, we find the Born approximation

$$4 \pi^2 f^2 s({\svector x})^2 \Delta G ({\svector x}, {\svector x}' ,f)+ \nabla^2\Delta G ({\svector x}, {\svector x}' ,f) = -8 \pi^2 f^2 s({\svector x}) G ({\svector x}, {\svector x}' ,f)\Delta s({\svector x}),$$ (5)

with the solution (using (4))

$$\Delta G ({\svector x}, {\svector x}'' ,f)= \iiint 8 \pi^2 f^2 s(... ...( {\svector x}' , {\svector x}'' ,f)\Delta s( {\svector x}' ) d {\svector x}' ,$$ (6)

which defines

$$\frac{\partial G ({\svector x}, {\svector x}'' ,f)}{\partial s( {... ...x}' ) G ({\svector x}, {\svector x}' ,f)G ( {\svector x}' , {\svector x}'' ,f).$$ (7)

The equivalent derivative for the wavefield is

$$\frac{\partial p ( {\svector x}, f )}{\partial s( {\svector x}' )} = \iiint \frac{\partial G ({\svector x}, {\svector x}'' ,f)}{\partial s( {\svector x}' )} w( {\svector x}'' ,f) d {\svector x}''$$ (8)

$$= \iiint 8 \pi^2 f^2 s( {\svector x}' ) G ({\svector x}, {\svector ... ...)G ( {\svector x}' , {\svector x}'' ,f)w( {\svector x}'' ,f) d {\svector x}'' .$$ (9)

Now that we have a conventional Born perturbation, we can compare the Rytov.

The Rytov Approximation and Wavepaths

To arrive at a Rytov approximation, take the wavefield derivative (8) and substitute the phase of the Green's function, defined by the complex logarithm

$$\theta ({\svector x}, {\svector x}'' ,f) \equiv \log G ({\svector x}, {\svector x}'' ,f), { ~ \rm and }$$ (10)

$$\frac{\partial \theta ({\svector x}, {\svector x}'' ,f)}{\partial s( {\svector x}' )} G ({\svector x}, {\svector x}'' ,f) = \frac{\partial G ({\svector x}, {\svector x}'' ,f)}{\partial s( {\svector x}' )} ,$$ (11)

allowing us to rewrite the wavefield derivative (8) as

$$\frac{\partial p ( {\svector x}, f )}{\partial s( {\svector x}' )} = \iiint \frac{\partial \theta ({\svector x}, {\svector x}'' ,f)}{\... ... )} G ({\svector x}, {\svector x}'' ,f)w( {\svector x}'' ,f) d {\svector x}'' ,$$ (12)

$${\rm where ~~} \frac{\partial \theta ({\svector x}, {\svector x}'' ,f)}{\partial s( {\svector x}' )} = 8 \pi^2 f^2 s( {\svector x}' ) \frac{G ({\svector x}, {\svector x... ...)G ( {\svector x}' , {\svector x}'' ,f)}{G ({\svector x}, {\svector x}'' ,f)} .$$ (13)

In effect, we have only regrouped the wavefield derivative (9) after multiplying and dividing by the same Green's function. This regrouping is crucial because further approximations will be applied to the phase derivative (13) rather than to the Green's function derivative (7).

Marta Woodward (1989, Ph.D thesis, Stanford) defines a wavepath as the imaginary part of the phase derivative (13). She plotted the wavepath function over different perturbed slowness positions ${\svector x}'$ while holding frequency $f$ and endpoints ${\svector x}$ and ${\svector x}''$ constant. The wavepath is stationary along a Fermat raypath, but also has a measurable Fresnel width. A multi-frequency wavepath includes an arbitrary source function $w( {\svector x}'' ,f)$.

We reintroduce the source by defining the phase of the pressure field

$$\phi ( {\svector x}, f ) \equiv \log p ( {\svector x}, f )$$ (14)

so that

$$\frac{\partial \phi ( {\svector x}, f )}{\partial s( {\svector x}' )} = \frac {\iiint 8 \pi^2 f^2 s( {\svector x}' ) G ({\svector x}, {\s... ... ({\svector x}, {\svector x}''' ,f)w( {\svector x}''' ,f) d {\svector x}''' } .$$ (15)

The numerator is the wavefield derivative (8), and denominator is the wavefield $p ( {\svector x}, f )$ from (4)

In reports from the University of Utah Tomography and Migration consortium, G. Schuster and Yi Luo proposed picking and fitting unwrapped phases rather than traveltimes for seismic tomography. This gradient (15) allows a descent optimization.

The Rytov approximation is defined as a linearization of the phase of a wavefield with respect to a particular model parameter such as slowness. Here, the Rytov approximation is given by the imaginary part of the gradient (15).

The Eikonal approximation

By next adding the Eikonal approximation, we see more intuitively how Rytov perturbations alter the phase and time lag of a wavefield. We can also modify raypath-based tomography to honor the resolution of a source with limited bandwidth.

To apply the Eikonal approximation, assume that the Green's function can be described by a function of the form

$$G ({\svector x}, {\svector x}' ,t) \equiv u ({\svector x}, {\svector x}' )\delta[ t - \tau ({\svector x}, {\svector x}' )]$$ (16)

$${\rm ~ or ~ } G ({\svector x}, {\svector x}' ,f) \equiv u ({\svector x}, {\svector x}' )\exp [ -i 2 \pi f \tau ({\svector x}, {\svector x}' )] .$$ (17)

$u$ and $\tau$ are real and smooth functions. Note that

$$p ( {\svector x}, t )= \iiint u ({\svector x}, {\svector x}' )w [ {\svector x}' , t - \tau ({\svector x}, {\svector x}' )] d {\svector x}' .$$ (18)

Substitute the approximation (17) into the Green's function definition (3) for a homogeneous equation. The terms scaled by different powers of frequency must each vanish, giving the Eikonal equation

$$\left \vert {\svector \nabla}\tau ({\svector x}, {\svector x}' )\right \vert^2 = s ({\svector x})^2 ,$$ (19)

and transport equations,

$$\nabla^2u ({\svector x}, {\svector x}' ) = - \delta ( {\svector x}- {\svector x}' )$$ (20)

$${\rm ~ and ~ } 2 {\svector \nabla}\tau ({\svector x}, {\svector x}' )\cdot {\svector \nabla}u ({\svector x}, {\svector x}' ) + u ({\svector x}, {\svector x}' )\nabla^2\tau ({\svector x}, {\svector x}' )= 0 .$$ (21)

Gradients ${\svector \nabla}$ and Laplacians $\nabla^2$ are taken with respect to the unprimed position ${\svector x}$. The Eikonal traveltime (phase delay) can be extrapolated from the Eikonal equation alone if the corresponding term dominates. This condition is satisfied if

$$\left \vert {\svector \nabla}v ({\svector x}) \right\vert \ll 2 \pi \left\vert f \right\vert .$$ (22)

In simple media, the transport equations can be replaced by geometric spreading factors. Reciprocity allows us to swap the position arguments in the Green's functions; thus, we can extrapolate the Eikonal (19) and transport equations (20) and (21) and in whichever direction is convenient.

Applying the approximation (17) to the phase linearization (13), we find

$$\frac{\partial \theta ({\svector x}, {\svector x}'' ,f)}{\partial s( {\svector x}' )} = 8 \pi^2 f^2 s ( {\svector x}' ) \frac{ u ({\svector x}, {\svector... ... ( {\svector x}' , {\svector x}'' )- \tau ({\svector x}, {\svector x}'' )] \} ,$$ (23)

$$\frac{\partial \theta ({\svector x}, {\svector x}'' ,t)}{\partial s( {\svector x}' )} = -2 s ( {\svector x}' ) \frac{ u ({\svector x}, {\svector x}' )u (... ...tor x}' , {\svector x}'' )+ \tau ({\svector x}, {\svector x}'' )] , {\rm and ~}$$ (24)

$$\frac{\partial G ({\svector x}, {\svector x}'' ,t)}{\partial s( {\svector x}' )} = -2 s ( {\svector x}' ) u ({\svector x}, {\svector x}' )u ( {\svec... ...tau ( {\svector x}' , {\svector x}'' )+ \tau ({\svector x}, {\svector x}'' )] .$$ (25)

The last equation uses the derivative (11).

The linearization (8) for the wavefield can be rewritten as

$$\frac{\partial p ( {\svector x}, t )}{\partial s( {\svector x}' )} = \int\iiint \frac{\partial G ({\svector x}, {\svector x}'' ,t)}{\partial s( {\svector x}' )} w ( {\svector x}'' , t - t' ) d {\svector x}'' dt'$$

$$= \iiint -2 s ( {\svector x}' ) u ({\svector x}, {\svector x}' )u (... ...}' , {\svector x}'' )+ \tau ({\svector x}, {\svector x}'' )] d {\svector x}'' .$$ (26)

The time lag in this result (26) is the time difference between the fastest path between two points and the fastest path passing through any scattering point ${\svector x}'$ where slowness is perturbed. This kernel integrates all slowness perturbations to get a perturbation of the wavefield, as we require for tomographic inversion. Unlike the Born approximation, these scatterers cause phase delays instead of reflections.

Point Sources

Finally we can examine a point source to see how this linearization affects a single two-point raypath. A point source requires that

$$w ( {\svector x}, f )= w(f) \delta ( {\svector x}- {\svector x}_0), {\rm ~~ then ~} p ( {\svector x}, f )= G ({\svector x},{\svector x}_0,f) w(f) .$$ (27)

With the Eikonal approximation, we find

$$\frac{\partial p ( {\svector x}, t )}{\partial s( {\svector x}' )} = \int \frac{\partial \theta ( {\svector x}, {\svector x}_0, t - t')}{\partial s( {\svector x}' )} p ({\svector x},t') dt'$$ (28)

$$= -2 s( {\svector x}' ) \frac{ u ({\svector x}, {\svector x}' )u ( ... ...\tau ( {\svector x}' , {\svector x}_0) + \tau ( {\svector x}, {\svector x}_0) ]$$ (29)

$$= -2 s( {\svector x}' ) u ({\svector x}, {\svector x}' )u ( {\svect... ...au ( {\svector x}' , {\svector x}_0) + \tau ( {\svector x}, {\svector x}_0) ] .$$ (30)

We can easily compute this linearization from extrapolations of the Eikonal and transport equations. Reciprocity allows us to extrapolate the Eikonal traveltimes from each endpoint of a particular path.

The phase, width, and contours of the two-point wavepath (30) are dominated by the second-time derivative of the waveform $\ddot{w}(t)$. The stationary points on this wavepath are along the fastest raypath. Lags are simply the difference in time between a perturbed path and the fastest path. The amplitude scaling ( $u ({\svector x}, {\svector x}' )$) controls the decay of the wavefield strength away from the fastest path.

Numerically one could calculate this wavepath from explicit traveltime extrapolation methods, from either finite-differences or ray methods. Extrapolate traveltimes from both endpoints to the entire region of interest and sum the two tables. Subtract the minimum value from this total, and the wavefield values (30) as a function of this lag.