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

\begin{document}  
\bibliographystyle{plain}
\markright{A convenient anisotropy approximation --- W.S. Harlan}

\title{
A convenient approximation of transverse
isotropy for higher-order moveout,
prestack time migration, 
and depth calibration}

\author{William S. Harlan}

\date{August 1998}
\maketitle

\section{Introduction}
Processors use different RMS velocity models for 
three steps of time imaging: NMO, DMO, and poststack time migration.
To perform prestack time migration in a single step, we must
use a single velocity model.  A single step avoids an extra
stationary-phase approximation and should produce more accurate
results.  Nevertheless, results are usually worse with a single 
velocity model, unless different velocities are used for flat and dipping 
reflections \cite{SEG-1991-1151}.  Those velocities which best
fit prestack normal moveouts over offset (flat reflections) do not best 
focus the tails of diffractions over midpoint (dipping reflections).  
Conventional processing hides this difference with inconsistent 
velocity models for prestack moveout analysis and poststack migration.

Occasionally, processors want to use a higher-order 
normal-moveout equation to flatten prestack reflections with long offsets.
Conventional moveout analysis does a good job of fitting the
difference in traveltime between near and far offsets,
but a higher-order moveout can better flatten any residual
bulge in the middle.  Our parameterization of this normal moveout 
should be consistent with the model of anisotropy 
used in full prestack time imaging.

It is recognized that the kinematics of surface reflection seismic data
are insensitive to component of transverse isotropy that
is essential for an accurate conversion of time to depth.
We can isolate the parameters needed to fit surface 
reflection times and allow depths to be calibrated independently.

A single convenient analytic approximation of transverse isotropy
will allow us to generalize prestack moveout, time imaging, and
depth conversion.  No more parameters will be introduced than
necessary.  Parameters are decoupled so that each can be
estimated in turn, if the additional degrees of freedom are
required to fit the data.  

Much of this material originally appeared in an appendix to
Harlan \cite{ha:95a}.

\section{Parameters for approximate transverse isotropy}
\label{app:aniso}
Assume that anisotropic velocities have a vertical axis of symmetry,
like the transversely isotropic (TI) media described in 
Thomsen \cite{GEO51-10-19541966}.  Although that paper is
titled ``Weak elastic anisotropy,'' the same parameterizations
can be applied to very strong anisotropy \cite{GEO59-08-12901304}.

Three of Thomsen's parameters, $V_z$, $\delta$, and $\epsilon$,
are defined by the elastic constants of a general TI medium.
These constants can be used to specify three different
effective velocities at a single point in the model.
$V_z$ is the velocity of a wave traveling vertically along
the axis of symmetry.  The velocity in any horizontal direction
$V_x$ is defined by
\begin{eqnarray}
\label{eq:epsilondef}
\epsilon &=& V_x^2 (V_z^{-2} - V_x^{-2})/2 , \\
V_x^2 &=& (1 + 2 \epsilon) V_z^2, {\rm ~ and } \\
V_x &\approx& (1+\epsilon) V_z . 
\end{eqnarray}
A ``normal moveout velocity'' (NMO) velocity $V_n$ is defined by
\begin{eqnarray}
\delta &\equiv& V_n^2 (V_z^{-2} - V_n^{-2})/2 , \\
V_n^2 &=& (1 + 2 \delta) V_z^2, {\rm ~ and } \\
V_n &\approx& (1+\delta) V_z .
\label{eq:deltadef}
\end{eqnarray}
If these TI properties represent the equivalent medium of
many isotropic layers \cite{backus-equivalent-media,GEO54-05-05810589}, 
then we can expect $\epsilon > \delta$ \cite{Berryman79}.
Using Backus averaging, Phil Anno of Conoco has also shown that
we can expect $\epsilon > 0$ and $\delta < 0$, if
the $V_s/V_p$ ratio and 
$V_s$ have a positive correlation.  These inequalities imply
that $V_n \leq  V_z \leq V_x$.  Well calibrations have shown
that shales can produce $\delta > 0$.  Such shales
possess ``intrinsic'' anisotropy that cannot be
modeled as the macroscopic equivalent of isotropic layers.

Researchers at the Colorado School of Mines
\cite{GEO59-08-12901304,GEO60-05-15501566} have also defined a constant
\begin{eqnarray}
\label{eq:etadefa}
\eta  &\equiv& (\epsilon - \delta)/(1 + 2 \delta) \\
\label{eq:etadefaa}
 &=& V_x^2 (V_n^{-2} - V_x^{-2})/2 ,\\
\label{eq:etadefb}
V_x^2 &=& (1 + 2 \eta) V_n^2, {\rm ~ and } \\
V_x &\approx& ( 1 + \eta ) V_n.
\label{eq:etadefc}
\end{eqnarray}
For an equivalent medium of isotropic layers $\eta > 0$.

Many combinations of three of these parameters 
$V_z, V_x, V_n, \epsilon, \delta, \eta$ can be used to describe
a TI medium for compressional P waves with a known axis of symmetry. 
Such an approximation has already dropped a fourth constant (shear
wave velocity) to which compressional waves are insensitive.

The parameter $\eta$ is an excellent parameter choice for maximum
sensitivity to the kinematics of surface measurements only.
If TI velocities are parameterized by $V_x$, $\eta$, and $\delta$,
we find that surface measurements are very insensitive to $\delta$.
We could not make the same claim if $\epsilon$ were used instead of
$\eta$.

I prefer the three velocities $V_z, V_x,$ and $V_n$ because
they share the same units.  Surface measurements are very insensitive
to $V_z$, given values for $V_x$ and $V_n$.
This choice is not the most convenient for processing, however.

\section {Phase and group velocities}
The exact equations for TI phase velocity as a function of
angle are rather clumsy, and {\sl no explicit form} is available
for group velocity.  Explicit approximate equations can 
fit the same family of curves almost as well
as the original correct equations \cite{GPR41-04-03810412}.  
I use an approximate equation for group velocity
which appears to emulate closely the exact curves for
large ranges of positive $\epsilon$ and negative $\delta$.
Estimated curves usually have larger statistical errors from noisy data
than introduced by these approximations.  

Kirchhoff migrations can calculate traveltimes by integrating 
the group velocity along a stationary Fermat raypath.
Fourier-domain implementations can use only
the equation for phase velocities.  Some raytracing methods
use both, because phases must be matched across
discontinuous boundaries.

I choose approximate curves with the three velocities $V_z, V_x,$ and $V_n$.
Let $\phi$ be the group angle of a raypath from the vertical.
Then the group velocity $V(\phi)$ can be expressed as 
\begin{eqnarray}
\label{eq:group}
V(\phi)^{-2} &=& V_z^{-2}  \cos^2 \phi
+ ( V_n^{-2} - V_x^{-2} ) \cos^2 \phi \sin^2 \phi
+  V_x^{-2}  \sin^2 \phi ,\\
V(\phi)^{-1} &=& V_x^{-1} 
\sqrt{ 1+2 \eta \cos^2 \phi  \sin^2 \phi +2 \epsilon \cos^2 \phi } \\
\label{eq:groupapprox}
&\approx& V_x^{-1} 
( 1+ \eta \cos^2 \phi \sin^2 \phi + \epsilon \cos^2 \phi ) \\
\label{eq:groupapproxc}
&\approx& V_x^{-1} 
[ 1+ \eta \cos^2 \phi  (1 + \sin^2 \phi) + \delta \cos^2 \phi ] .
\end{eqnarray}
Compare Byun et al \cite{GEO54-12-15641574}, who use
the same approximation with different parameters.

Greg Lazear of Conoco found that a symmetric equation
approximates the phase velocity $v(\theta)$
as a function of the phase angle $\theta$,
but with reciprocals of the same velocity parameters:
\begin{eqnarray}
v(\theta)^{2} = V_z^{2} \cos^2 \theta 
+ ( V_n^{2} - V_x^{2} ) \cos^2 \theta  \sin^2 \theta 
+ V_x^{2}  \sin^2 \theta .
\label{eq:phasevel}
\end{eqnarray}
Compare this phase equation closely to the group equation (\ref{eq:group}).
I know of no other approximation that allows such symmetry.

\section{Extended stacking moveouts}
The normal-moveout (NMO) velocity $V_n$ has a physical
interpretation to justify its name.  
Imagine an experiment on a homogeneous
and anisotropic medium, or imagine a small scale experiment
on a smooth model.  Measure the traveltime $t_0$ between two
points placed on a vertical line, separated by a vertical
distance $z = V_z t_0$.  Now displace the
upper point a distance $h$ along a horizontal line (a normal-moveout)
and measure the new traveltime $t_h$.   

Then according to equation (\ref{eq:group}) the traveltime $t_h$
as a function of offset $h$ is exactly
\begin{eqnarray}
t_h^2 = t_0^2 + 
\left [ V_n^{-2} + (V_x^{-2} - V_n^{-2}) { h^2 \over {h^2 + V_z^2 t_0^2}}
\right ] h^2 .
\label{eq:nmoeq}
\end{eqnarray}
For small offsets $h \ll V_z t_0$, the value of $t_h$ in this
``moveout equation'' is 
dominated by the NMO velocity $V_n$ rather than $V_x$.  For large
offsets $h \gg V_z t_0$, the raypath is almost horizontal
and $V_x$ dominates.

I find it convenient to define a stacking velocity $V_h(h)$
as a function of the offset $h$ for a fixed vertical distance
$z = V_z t_0$:
\begin{eqnarray}
V_h(h)^{-2} &\equiv&  (t_h^2 - t_0^2)/h^2 \\
\label{eq:stackdefv}
&=& V_n^{-2} + (V_x^{-2} - V_n^{-2}) { h^2 \over {h^2 + V_z^2 t_0^2}} \\
\label{eq:stackdef}
&=& V_n^{-2} \left( 1 - 
\frac{2 \eta}{1 + 2 \eta} \cdot \frac{h^2}{h^2 + V_z^2 t_0^2}\right) .
\end{eqnarray}
I use the term stacking velocity because I want to suggest
the best-fitting curve over a finite range of offsets, as
you would prefer for a stacking or semblance analysis.

Simplify the moveout equation (\ref{eq:nmoeq}) to fit a pseudo-hyperbola:
\begin{eqnarray}
t_h^2 = t_0^2 + h^2 / V_h(h)^2 .
\label{eq:stackeq}
\end{eqnarray}
The stacking velocity covers the range 
$V_n \leq V_h(h) \leq V_x$ for a Backus equivalent medium with 
negative $\delta$, increasing in value as $h$ increases.
When $\eta = 0$, the curve is exactly hyperbolic, and $V_n = V_x$.
Notice that this stacking velocity can measure a local property
as well as an average to the surface.
To use two-way reflection times in (\ref{eq:stackeq})
we need only replace the half offset
$h$ by the full offset.

Three measurements of traveltimes at three
different offsets $h$ should uniquely determine the three
velocity constants {$V_z$, $V_x$, $V_n$}.  The
traveltimes are much more sensitive to $V_n$, which determines
moveouts at small offsets, and to $V_x$, which determines
moveout at larger offsets.  The vertical velocity $V_z$ affects
only the rate at which the stacking velocity (\ref{eq:stackdefv}) 
changes from one limit to the other.  As long as $V_z$ has roughly the
correct magnitude, then we can fit all measured traveltimes
very well.  

For imaging surface data in time, we acknowledge our insensitivity to
$V_z$ and can approximate it with another value.  
We can approximate $V_z \approx V_n$ and simplify
stacking velocity (\ref{eq:stackdef}) even further, as in 
\begin{eqnarray}
\label{eq:vhapprox}
V_h(h)^{-2} \approx V_n^{-2} \left( 1 - 2 \eta \frac{h^2}{h^2 + V_n^2 t_0^2} \right) .
\end{eqnarray}
This new equation depends only on two parameters, $V_x$ and $V_n$.
A better approximation might be $V_z \approx V_n^2/V_x$, or
equivalently $\epsilon \approx 2 \delta$, which is probably
closer to commonly observed values.  
NMO equation (\ref{eq:nmoeq}) now is equivalent to
\begin{eqnarray}
\label{eq:grechka}
t_h^2 \approx t_0^2 + 
\left [ 
 V_n^{-2} + (V_x^{-2} - V_n^{-2}) { h^2 \over {h^2 + V_n^4 V_x^{-2} t_0^2}}
\right ] h^2 .
\end{eqnarray}
which is equivalent to equation (5) in
Alkhalifah \cite{alk-97} and
equation (7) in Grechka and Tsvankin 
\cite{grechka-tsvankin}.  Both these publications derive from 
Tsvankin and Thomsen \cite{GEO59-08-12901304}, which
uses an asymptotic correction of a Taylor expansion to arrive
at this approximation.

The equivalent stacking velocity is then
\begin{eqnarray}
\label{eq:vhapproxtwo}
V_h(h)^{-2} \approx V_n^{-2} \left( 
1 - 2 \eta \frac{h^2}{h^2 + V_n^4 V_x^{-2} t_0^2} 
\right) .
\end{eqnarray}
The differences between these two approximate 
stacking velocities (\ref{eq:vhapproxtwo})
and (\ref{eq:vhapprox}) are negligible for numerical work.  
I will use the simpler version (\ref{eq:vhapprox}).

\def \Vi {V_{\rm iso}}
\def \hm {h_{\rm max}}
Moveout analyses determine the stacking velocity for a specific aperture of 
offsets.  Define our best {\em isotropic} approximation of the velocity
to be the stacking velocity $\Vi$ at the maximum offset $\hm$
of the aperture:
\begin{eqnarray}
\label{eq:viso}
\Vi &\equiv& V_h(\hm) \\
&\approx& V_n \left( 1 + \eta \frac{\hm^2}{\hm^2 + V_n^2 t_0^2} \right) .
\end{eqnarray}
Or we can assume that we know the aperture angle $\alpha$ from the vertical,
so that
\begin{eqnarray}
\label{eq:vivn}
\Vi &\approx& V_n ( 1 + \eta \sin^2 \alpha ) ,
{\rm ~ where ~}  \tan \alpha = h / V_z t_0 .
\end{eqnarray}
The best isotropic velocity $\Vi$ is a simple function of 
infinitesimal-offset NMO velocity $V_n$ and the anisotropic parameter 
$\eta$.  
Similarly we can use the definition of $\eta$ in (\ref{eq:etadefaa})
to rewrite the above equation (\ref{eq:vivn}) as 
\begin{eqnarray}
\Vi \approx V_x ( 1 - \eta \cos^2 \alpha ) .
\label{eq:vxvi}
\end{eqnarray}
This form will prove to be very useful when rewriting 
our group velocity equations 
(\ref{eq:groupapprox}) and (\ref{eq:groupapproxc}).

A conventional velocity analysis produces
densely picked values for $\Vi$.  
The anisotropic parameter $\eta$ adjusts the
moveout between near and far offsets.
For anisotropic moveout analysis, we could use the
following offset-dependent stacking velocity to scan for $\eta$,
holding $\Vi$ constant:
\begin{eqnarray}
V_h(h)^{-2}  
\approx \Vi^{-2} 
\left[ 1 + 2 \eta 
\left(
\frac{\hm^2}{\hm^2 + \Vi^2 t_0^2} 
- \frac{h^2}{h^2 + \Vi^2 t_0^2} 
\right) 
\right] .
\end{eqnarray}
Conventional moveout analysis is not very sensitive
to the anisotropic parameter $\eta$ except
for unusually wide-aperture data, with offsets greater than depth.  
Much more anisotropic information is available by performing a full prestack 
time migration.

For prestack time migration, we can expect that conventional
analysis for $\Vi$ will best describe the moveouts of flat
reflections.  Dipping reflections are difficult to pick in
prestack semblance analysis because they are sparser and
move across midpoints, from gather to gather, as migration
velocity changes.

Holding $\Vi$ constant, we can perturb $\eta$ at a low
spatial resolution until the imaging of steep reflections improves.
As $\eta$ changes, there will be a small adjustment of the
bulge in flat reflections over offset, with little effect
on a fully time-migrated stack.  The moveouts of dipping reflections, 
on the other hand, will change drastically as $\eta$ changes, with swings
from positive to negative residual moveouts, and with lateral
movement over midpoint.  As $\eta$ improves, you should see fault-plane 
reflections sharpen and focus in targeted prestack time-migrated images
%(See thesis work by Tang).
Flat reflections should change little, and $\Vi$ should require little 
revision after updating $\eta$.  By contrast, an optimization of 
$V_n$ and $\eta$ requires both to be adjusted simultaneously,
with equal resolution.

\section{Adjustment for depth ties}
The group velocity equation (\ref{eq:groupapprox}) is useful
for Kirchhoff depth imaging.
We can focus images very well with good values for $V_x$ and $V_n$,
(or $\Vi$ and $\eta$),
then adjust imaged depths to tie wells with $V_z$ 
(or $\epsilon$, or $\delta$) while
holding the other two parameters constant.

If we have used tomographic methods to estimate 
isotropic {\em interval} velocities, then
we should attempt also to estimate the effective aperture angle $\alpha$
at each depth in the model.

To recapituate our algebra, we combine the definitions of 
$\eta$ (\ref{eq:etadefb}) and of stacking velocity $V_h (h)$
(\ref{eq:stackdef}) to solve for the horizontal velocity:
\begin{eqnarray}
V_x^{-1} 
&\approx& V_h(\hm)^{-1} 
\left ( 1 - \eta { {V_z^2 t_0^2} \over {\hm^2 + V_z^2 t_0^2}} \right ) \\
&\approx& \Vi^{-1} ( 1 - \eta \cos^2 \alpha ) .
\end{eqnarray}
At maximum offset, we recover our previous 
relationship (\ref{eq:vxvi}) between $\Vi$ and $V_x$.
Substitute this horizontal velocity into the group velocity 
(\ref{eq:groupapprox}).
We can then adjust the isotropic velocities $\Vi$ with $\eta$ and
$\epsilon$ according to
\begin{eqnarray}
\label{eq:stackb}
V( \phi )^{-1} \approx \Vi^{-1} 
(1 - \eta \cos^2 \alpha ) 
(1 + \eta \cos^2 \phi  \sin^2 \phi + \epsilon \cos^2 \phi ) . 
\end{eqnarray}
The isotropic velocity $\Vi$ best explains the moveouts and traveltimes 
of relatively flat reflections over the finite aperture.  
The parameter $\eta$ modifies these velocities at high dips, to image steep 
reflections better without degrading the imaging of low-dip reflections.
The third parameter $\epsilon$ has little effect on measured
surface traveltimes at any dip (holding $\Vi$ and $\eta$ constant), 
but can be adjusted as necessary to tie wells.  
We can also use layered medium theory to 
predict this $\epsilon$ from estimated
$\eta$ and $V_x$.  Or if shale dominates, with strong
intrinsic anisotropy and $\delta > 0$, then correlations
can be calibrated for a given area.
At worst, we know $0 < \epsilon$, so we can assume a default
value of $\eta/2 < \epsilon  < 2 \eta$, as appropriate
for a given area.  Such a default value is
still better than a default value of $0$.   

These three parameters $\Vi$, $\eta$, and $\epsilon$
do not completely decouple the steps of anisotropic 
velocity analysis, but they should minimize the number
of iterations necessary for revisions.

If you prefer to use $\delta$ instead of $\epsilon$
as the third degree of freedom, then simply use
the definition of $\eta$ in (\ref{eq:etadefa}) for

\begin{eqnarray}
\label{eq:stackc}
V( \phi )^{-1} \approx \Vi^{-1} 
(1 - \eta \cos^2 \alpha ) 
[1 + \eta \cos^2 \phi (1+ \sin^2 \phi) + \delta \cos^2 \phi ] . 
\end{eqnarray}

\section{Adjustment of narrow-aperture velocities}

Although I greatly prefer the approach in the preceding
section, many prefer to treat their estimated isotropic
velocities as equivalent to NMO velocities $V_n$.
Such an assumption is not a bad one if angles are limited
during interval velocity estimation.  Dix inversion
of stacking velocities may be closer to $V_n$ if stacking
velocities were consciously optimized for inner offsets.
The Common Reflecting Surface tomography of Karlsr\"uhe University
inverts only the curvature of reflection traveltimes
around zero-offsets.

For such approaches, one might prefer a different triplet
of velocity parameters: $V_n$, $\eta$, and $\delta$.
Group velocity can be described by replacing $V_x$ in 
the approximation (\ref{eq:groupapproxc}) with 
$V_n$ and $\eta$ as in equation (\ref{eq:etadefc}).
Additionally we can replace $\epsilon$ with $\delta$,
using the definition of $\eta$ in (\ref{eq:etadefa}):

\begin{eqnarray}
V(\phi)^{-1} 
&\approx& V_n^{-1} 
( 1 - \eta \sin^4 \phi + \delta \cos^2 \phi ) \\
&\approx& V_n^{-1} 
[ 1 - \eta (1 - \cos^2 \phi \sin^2 \phi ) + \epsilon \cos^2 \phi ] .
\end{eqnarray}

Again $V_n$ and $\eta$ should be sufficient to model
all surface reflection traveltimes, for all dips.
Holding these two constant, we can adjust either $\delta$
or $\epsilon$ to tie known well depths.

\section{Acknowledgments}
Many thanks to Greg Lazear and Phil Anno, who spent 
months with me at Conoco trying to reconcile and simplify 
different parameterizations of transverse isotropy.
Thanks to Andreas R\"uger for a very careful review of this paper.
\bibliography{wsh}

\end{document}

