\documentclass[12pt]{article}
\usepackage{wsh}
\usepackage{amsmath}

\setlength{\textwidth}{\paperwidth}
\addtolength{\textwidth}{-1.6in}
\addtolength{\oddsidemargin}{-0.8in}
\setlength{\evensidemargin}{\oddsidemargin}

\begin{document}  
\bibliographystyle{plain}

\markright{Rytov/Eikonal Wavepaths --- W.S. Harlan}

\title{Rytov/Eikonal Approximation of Wavepaths}

\author{William S. Harlan}

\date{April 1998}

\maketitle

\def \xv {{\svector x}} 
\def \xp { \xv ' }
\def \xpp { \xv '' }
\def \xppp { \xv ''' }

\def \pftwo{ \pi^2 f^2 }
\def \xz {\xv_0}
\def \grad {{\svector \nabla}}
\def \laplace {\nabla^2}

\def \xt {( \xv , t )}
\def \xf {( \xv , f )}
\def \xxf {(\xv,\xp,f)}
\def \xxpf {(\xv,\xpp,f)}
\def \xpxpf {(\xp,\xpp,f)}
\def \xxppf {(\xv,\xppp,f)}
\def \xxt {(\xv,\xp,t)}
\def \xxpt {(\xv,\xpp,t)}
\def \xx {(\xv,\xp)}
\def \xxp {(\xv,\xpp)}
\def \xpx {(\xp,\xpp)}
\def \xxz { ( \xv , \xz ) }
\def \xpxz { ( \xp , \xz ) }

\section{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.

\section{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.

\begin{eqnarray}
v(\xv )^{-2} \ddot{p} \xt  
	- \laplace p \xt &=&  w \xt ; \\
{\rm or } ~~ 4 \pftwo s(\xv)^2 p \xf + \laplace p \xf &=&
- w \xf .
\end{eqnarray}

The pressure $p$ is the non-zero diagonal of the strain tensor.  The vector
$\xv$ 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.

\begin{eqnarray}
\label{eqgreen}
4 \pftwo s(\xv)^2 G \xxf &+& \laplace G \xxf = 
	- \delta ( \xv - \xp ), \\
\label{eqgreenp}
{\rm so ~ that } ~~ p \xf &=& \iiint G \xxf w (\xp,f ) d\xp .
\end{eqnarray}

If we perturb the Green's function definition (\ref{eqgreen}) 
and discard second-order terms, we find the Born approximation

\begin{eqnarray}
 4 \pftwo s(\xv)^2 \Delta G \xxf + \laplace \Delta G \xxf &=&
-8 \pftwo s(\xv) G \xxf \Delta s(\xv),
\end{eqnarray}
with the solution (using (\ref{eqgreenp}))
\begin{eqnarray}
\Delta G \xxpf  = 
   \iiint 8 \pftwo s(\xp) G \xxf G \xpxpf \Delta s(\xp) d\xp ,
\end{eqnarray}
which defines
\begin{eqnarray}
\label{eqderg}
\frac{\partial G \xxpf}{\partial s(\xp)} =8 \pftwo s(\xp) G \xxf G \xpxpf .
\end{eqnarray}
The equivalent derivative for the wavefield is 
\begin{eqnarray}
\label{eqderp}
\frac{\partial p \xf}{\partial s(\xp)} &=& \iiint
\frac{\partial G \xxpf}{\partial s(\xp)} w(\xpp,f) d\xpp \\
&=& 
\label{eqderpp}
\iiint 8 \pftwo s(\xp) G \xxf G \xpxpf w(\xpp,f) d\xpp .
\end{eqnarray}
Now that we have a conventional Born perturbation, we can
compare the Rytov.

\section{The Rytov Approximation and Wavepaths}

To arrive at a Rytov approximation, take the wavefield derivative 
(\ref{eqderp}) and substitute the phase of the Green's function,
defined by the complex logarithm
\begin{eqnarray}
\theta \xxpf &\equiv& \log G \xxpf , { ~ \rm and } \\
\label{eqlogder}
\frac{\partial \theta \xxpf }{\partial s(\xp)} G \xxpf &=&
\frac{\partial G \xxpf}{\partial s(\xp)} ,
\end{eqnarray}
allowing us to rewrite the wavefield derivative (\ref{eqderp}) as
\begin{eqnarray}
\frac{\partial p \xf}{\partial s(\xp)} &=& \iiint
\frac{\partial \theta \xxpf }{\partial s(\xp)} G \xxpf w(\xpp,f) d\xpp, \\
\label{eqdert}
{\rm where ~~} 
\label{eqplin}
\frac{\partial \theta \xxpf }{\partial s(\xp)} &=&
8 \pftwo s(\xp) \frac{G \xxf G \xpxpf}{G \xxpf} .
\end{eqnarray}
In effect, we have only regrouped the wavefield derivative
(\ref{eqderpp}) after multiplying and dividing by the
same Green's function.  This regrouping is crucial because further
approximations will be applied to the phase derivative (\ref{eqdert})
rather than to the Green's function derivative (\ref{eqderg}).

Marta Woodward (1989, Ph.D thesis, Stanford) defines a wavepath
as the imaginary part of the phase derivative (\ref{eqdert}).
She plotted the wavepath function over different perturbed slowness
positions $\xp$ while holding frequency $f$ and endpoints 
$\xv$ and $\xpp$ 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(\xpp,f)$.

We reintroduce the source by defining the phase of the pressure field 
\begin{eqnarray}
\phi \xf &\equiv& \log p \xf
\end{eqnarray}
so that
\begin{eqnarray}
\label{eqphasegrad}
\frac{\partial \phi \xf}{\partial s(\xp)} &=& \frac
{\iiint 8 \pftwo s(\xp) G \xxf G \xpxpf w(\xpp,f) d\xpp}
{\iiint G \xxppf w(\xppp,f) d\xppp} .
\end{eqnarray}
The numerator is the wavefield derivative (\ref{eqderp}), and
denominator is the wavefield $p \xf$ from (\ref{eqgreenp})

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 (\ref{eqphasegrad}) 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 (\ref{eqphasegrad}).

\section{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
\begin{eqnarray}
G \xxt &\equiv& u \xx \delta[ t - \tau \xx ] \\
\label{eqapprox}
{\rm ~ or ~ }
G \xxf &\equiv& u \xx \exp [ -i 2 \pi f \tau \xx ] .
\end{eqnarray}
$u$ and $\tau$ are real and smooth functions.  Note that
\begin{eqnarray}
p \xt = \iiint u \xx w [ \xp , t - \tau \xx ] d\xp .
\end{eqnarray}
Substitute the approximation (\ref{eqapprox}) into the
Green's function definition (\ref{eqgreen}) for a homogeneous
equation.  The terms scaled by different
powers of frequency must each vanish, giving the Eikonal equation
\begin{eqnarray}
\label{eqeik}
\left | \grad \tau \xx \right |^2 = s (\xv)^2 ,
\end{eqnarray}
and transport equations,
\begin{eqnarray}
\label{eqtrans}
\laplace u \xx &=& - \delta ( \xv - \xp ) \\
\label{eqtranstwo}
{\rm ~ and ~ }
2 \grad \tau \xx \cdot \grad u \xx &+& u \xx \laplace \tau \xx = 0 .
\end{eqnarray}
Gradients $\grad$ and Laplacians $\laplace$ 
are taken with respect to the unprimed position $\xv$.
The Eikonal traveltime (phase delay) can be extrapolated from the Eikonal
equation alone if the corresponding term dominates.  This condition
is satisfied if
\begin{eqnarray}
\left | \grad v (\xv) \right| \ll 2 \pi \left| f \right| .
\end{eqnarray}
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 (\ref{eqeik}) and transport equations (\ref{eqtrans}) and 
(\ref{eqtranstwo}) and in whichever direction is convenient.

Applying the approximation (\ref{eqapprox}) to the phase linearization
(\ref{eqplin}), we find 
\begin{eqnarray}
\frac{\partial \theta \xxpf}{\partial s(\xp)} &=& 8 \pftwo s (\xp) 
\frac{ u \xx u \xpx}{u \xxp} \exp \{ -i 2 \pi f 
[ \tau \xx + \tau \xpx - \tau \xxp ] \} ,\\
\frac{\partial \theta \xxpt}{\partial s(\xp)} &=& -2 s (\xp)
\frac{ u \xx u \xpx}{u \xxp} 
\ddot{\delta} [ t - \tau \xx - \tau \xpx + \tau \xxp ] , {\rm and ~} \\
\frac{\partial G \xxpt}{\partial s(\xp)} &=& -2 s (\xp)
u \xx u \xpx \ddot{\delta} [t - \tau \xx - \tau \xpx + \tau \xxp ] .
\end{eqnarray}
The last equation uses the derivative (\ref{eqlogder}).

The linearization (\ref{eqderp}) for the wavefield can be rewritten as 
\begin{eqnarray}
\frac{\partial p \xt}{\partial s(\xp)} &=& \int\iiint
\frac{\partial G \xxpt}{\partial s(\xp)} w (\xpp , t - t' ) d\xpp dt' 
\nonumber \\
\label{eqwavelin}
&=& \iiint -2 s ( \xp ) u \xx u \xpx 
\ddot{w} [\xpp, t - \tau \xx - \tau \xpx + \tau \xxp ] d \xpp .
\end{eqnarray}
The time lag in this result (\ref{eqwavelin}) is 
the time difference between the
fastest path between two points and the fastest path passing through
any scattering point $\xp$ 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.

\section{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
\begin{eqnarray}
w \xf = w(f) \delta ( \xv - \xz ), {\rm ~~ then ~}
p \xf = G (\xv,\xz,f) w(f) .
\end{eqnarray}
With the Eikonal approximation, we find
\begin{eqnarray}
\frac{\partial p \xt}{\partial s(\xp)} &=& \int
\frac{\partial \theta ( \xv, \xz, t - t')}{\partial s(\xp)} p (\xv,t') dt' \\
&=& -2 s(\xp) \frac{ u \xx u \xpxz }{u \xxz } 
\ddot{p} [ \xv , t - \tau \xx - \tau \xpxz + \tau \xxz ] \\
\label{eqresult}
&=& -2 s(\xp) u \xx u \xpxz \ddot{w} [ t - \tau \xx - \tau \xpxz + \tau \xxz ] .
\end{eqnarray}
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 (\ref{eqresult})
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 \xx $) 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 (\ref{eqresult}) as 
a function of this lag.  

% \bibliography{wsh}
\end{document}