\def\figdir{./Figs}
One of the major difficulties with
depth tomography  is that our reflectors may be significantly mispositioned.
To illustrate this fact I constructed a synthetic model
(Figure~\ref{fig:synth-model}).  The model is an anticline with six reflectors,
one above the anticline, four within the anticline, and one flat reflector
representing
basement rock. For added difficulty there is a low velocity layer
between the second and third reflector.
The model was used to do acoustic wave modeling, with the
resulting dataset having 32 meter  common midpoint (CMP) 
spacing and   80 offsets spaced
64 meter apart.
If we use as our initial estimation of the
slowness  $s(z)$ function from outside the anticline,
the migrated reflector positions are pulled up
due to using too low  a velocity within
the anticline (Figure~\ref{fig:mig0}).
Our reflector mispositioning causes significant errors
in our linearized tomography problem.  For example, if we compare a set
of true raypaths  (Figure~\ref{fig:rays})
versus our estimated raypaths, we can see that we are
back projecting our moveout errors to the wrong positions in the slowness model.
Another way to think about it is that our linear approximation, equation
(\ref{eq:linear-approx}) is poor, causing us to have to take smaller
step sizes per
non-linear iteration, and increasing the chance that we will get stuck
in local minima or maxima.

\activeplot{synth-model}{height=3.3in,width=6.0in}{ER}{Left panel is  the 
synthetic velocity model with six reflectors spanning the anticline.
The right panel shows the data generated from this  model.}

\activeplot{mig0}{height=3.0in, width=6.0in}{CR}{Left panel is our initial guess
at the velocity function, the right panel shows the zero offset
ray parameter reflector
position using this migration velocity. The correct reflector
positions are shown as `*'.  Note that reflectors are
significantly mispositioned.}

\activeplot{rays}{height=3.0in, width=6.0in}{CR}
{The solid lines show raypath through the correct velocity and the correct
reflector position.  The dashed lines are raypaths through
our initial model and our initial reflector positions.  Note that
our estimated raypaths have significant error, therefore our tomography
operator will have significant error.}

\par
When doing time domain velocity analysis, reflector mispositioning
is not as much of a problem.
Reflector positioning in depth is done in
two stages  \cite{GPR25.04.07380745,GEO46.05.07340750}.
The first step is to find the velocity, $v_{\rm rms}$,
that best focuses the data.  
This can be relatively easier
in time because reflector position is much less effected by velocity.
After $v_{\rm rms}$ has been estimated, the reflectors are mapped
in depth using a velocity function, sometimes different than $v_{\rm rms}$.
The disadvantage of
time domain methods is that they cannot accurately define wave behavior
in complex media. This is clearly indicated in Figure~\ref{fig:time}.
A single hyperbola can not accurately describe the wave behavior even
in this relatively simple model.
\activesideplot{time}{width=3.0in,height=3.0in}{ER}{The solid curve represents
the correct travel times from $x=8$  along the bottom reflector to the
surface.  The
dashed lines show the arrival times that time migration would use for three
different
$v_{\rm rms}s$.
Note that both curves differ significantly
from the correct curve.}
\par
The tau domain offers a good alternative that combines
the best of both depth and time velocity analysis. The tau, or
vertical travel time, domain  is defined by
the transform
\beqa
\tau(z,x) &=& \int_0^z  2 s(z',x)  dz' \label{eq:transform1} \\
x'(z,x)  &=& x \label{eq:transform2} .
\eeqa
In this domain, flat reflectors are independent of the slowness
field and 
dipping reflector movement is much more limited
(Figure~\ref{fig:tau-initial}).
Additionally the coordinate transform  has not diminished our ability
to handle complex wave behavior.

\activeplot{tau-initial}{height=3.0in,width=6.0in}{ER}{Left panel is the same model as in Figure~\ref{fig:synth-model} in tau, the right is our initial guess in tau.
Note how our reflectors positions are much more accurate in the tau
case compared to the depth case (the right panel of Figure~\ref{fig:mig0}).}

%Before deriving the a ray based tau tomography it is useful to derive
%tau domain Eikonal equation
%in this domain.
%First, write the Eikonal equation in terms of both a mapping $\vmap$
%and a focusing $\sfoc$ slowness
%\beq
%{ 1 \over \smap(z,x)^2 }\left[{\partial t (z,x) \over \partial z}\right]^2 +
%%{ 1 \over \sfoc(z,x)^2} \left[{\partial t (z,x) \over \partial x}\right]^2  = 1 \label{eq:eik-depth} .
%\eeq

%We can substitute $\sigma_m$ into 
%(\ref{eq:eik-depth})  and
%end up with a new Eikonal equation in tau,
%\beq
%4\left[ { \part(\tau,x') \over  \partau } \right] +
%{1 \over \sfoc(\tau,x')^2 }\left[ {\part(\tau,x') \over  \parxp} | \sigma_m(\tau,x') 
%{\part (\tau,x') \over \partau} \right]^2  =1 .
%\eeq
%The advantage of the 
%new Eikonal equation is that while it depends directly on our focusing velocity,
%but  indirectly on the mapping velocity through $\sigma$.

\subsection{Tau operator}
In order to perform tomography in the tau domain, we must
find an operator equivalent to our depth domain operator 
(\ref{eq:tomo-base}) in the
tau domain.  This amounts to applying some simple coordinate transforms
to our already derived depth equations.
We will start from our relationship for 
a single ray segment, equation (\ref{eq:segment}).
We will begin  by rewriting it in terms of focusing and mapping slowness
($s_f$ and $s_m$):
\begin{equation}
\widetilde{dt} = 
\sqrt{\left(\widetilde{dz}\widetilde{S_m}\right)^2+
\left(\widetilde{dx}\widetilde{S_f}\right)^2} 
\label{eq:tau-segment}\end{equation}
where $\widetilde{}$ represents a ray segment quantity.

By applying the chain rule to (\ref{eq:transform1}) and 
(\ref{eq:transform2}) we find a relationship
between our derivatives:
 
\beqa
{\part \over \parz} &=& {\part \over \partau}{\partau \over \parz} +
{\part \over \parxp}{\parxp \over \parz} 
 = {\part \over \partau} {2 \smap(z,x)} \\
{\part \over \parx} &=& {\part \over \parxp}{\parxp \over \parx} +
{\part \over \partau}{\partau \over \parx} 
= {\part \over \parxp} + {\part \over \partau} \int_0^z {\partial \over \parx}
\left[{ 2 \smap(z',x) } \right] dz' .
\eeqa
Rather than carrying around the integral we can
define a differential mapping factor $\sigma$, 
\begin{equation}
\sigma(z,x)= \int_0^z {\partial \over \parx} \left[{ 2 \smap(z',x) } \right] dz' .\end{equation}
The differential mapping factor has non-zero values only in
areas at and below where the medium deviates from $s(z)$.   For the synthetic
model it is only inside the anticline (Figure~\ref{fig:sigma}).
\activesideplot[t]{sigma}{width=3.0in,height=3.0in}{ER}{Sigma values for the synthetic model in Figure~\ref{fig:tau-initial}}

We can use equations 
(\ref{eq:transform1}) and (\ref{eq:transform2}) to write a relationship
between change in $(z,x)$ and ($\tau,x')$ space,
\beqa
\widetilde{dz}&=&{\widetilde{d\tau} \over 2 \widetilde{S_m}}  - 
{\widetilde{\sigma} \widetilde{dx'} \over 2 \widetilde{S_m}} \\
\widetilde{dx}&=&\widetilde{dx} .
\eeqa
We  apply these definitions to equation (\ref{eq:tau-segment}) and come up
with a new definition for $dt$ in tau space,
\begin{equation}
\widetilde{dt} = 
\sqrt{\left(  
\left(\widetilde{S_f} \widetilde{dx'}\right)^2 +
{ \widetilde{d\tau} -  \widetilde{\sigma}
\widetilde{dx'} \over 2} \right)^2}  .\end{equation}
We  can then take the derivative with respect to $s_f$,
\begin{equation}
{\widetilde{dt} \over ds_f} ={ \widetilde{S_f} \widetilde{dx'}^2  \over 
\widetilde{dt}}
{d\widetilde{S_f} \over ds_f} - {(\widetilde{d\tau} - \widetilde{dx'} 
\widetilde{\sigma} ) 
\widetilde{dx'}
\over 4 \widetilde{dt}} {d\widetilde{\sigma} \over ds_f} .\end{equation}
Unfortunately, our operator isn't quite as simple as its depth counterpart.
In this case we do not have a linear relationship between 
$d\widetilde{S_f} \over ds_f$ 
and ${\widetilde{dt} \over ds_f}$ because of an
intermediate ${d\widetilde{\sigma}\over ds_f}$ term. 
We can  find a linear relationship between ${d\widetilde{\sigma}\over ds_f}$
and $d\widetilde{S_f} \over ds_f$ by first writing
      
\beqa
\frac {\partial \tau }{ \partial x'}&=&
\frac{\partial \tau}{\partial x}
\frac{\partial x}{\partial x'}
+ \frac{\partial \tau}{\partial z}
\frac{\partial z}{\partial x} \\ 
\frac{\partial \tau}{\partial x} &=&
\frac {\partial \tau }{ \partial x'} -
\frac {\partial}{\partial z}
\int_0^z 2 s(z',x) dz'  \nonumber \\
\sigma(z,x) &=& - 2 s_f(\tau, x')
\frac{\partial z}{\partial x}  \nonumber \\
\sigma(z,x) &=& -  s_f(\tau, x')
\int_0^\tau  \frac{\partial}{\partial x'} \frac{1}{ s_f(\tau',x')}  d \tau'  \nonumber .
\eeqa

We then take the partial derivative  with the respect to $s_f$ and evaluate
along the ray segment
\begin{equation}
\frac{\partial  \widetilde{\sigma}}{\partial s_f}  
= {\widetilde{\sigma}\over \widetilde{s_f}} \frac{\partial \widetilde{S_f}} { \partial s_f} -
s_f(\tau,x') \int_0^\tau \frac{\partial^2 s_f(\tau',x') }{\partial x'
\partial s_f} 
\partial \tau' .  \label{eq:dtauop} \end{equation}
With this relationship
(\ref{eq:dtauop}) we can then create a new complex tomography operator  in
tau space:
\begin{equation}
\dt \approx \tautomo \ds \label{eq:tomo-base-tau} .\end{equation}
