\documentclass{segabs}
\renewcommand{\figdir}{./Fig} % figure directory

\begin{document}
\title{Time-shift imaging condition in wave-equation tomography}

\renewcommand{\thefootnote}{\fnsymbol{footnote}} 
\author{Claudio Guerra and Biondo Biondi}

\footer{Wave-equation tomography with time-shift}

\righthead{Wave-equation tomography}

\maketitle
\lefthead{Guerra and Biondi}
\righthead{  }
\author{Claudio Guerra and Biondo Biondi}

\maketitle

\begin{abstract}
We introduce a new  objective function for wave-equation migration velocity analysis. This objective function borrows the main idea of the differential semblance optimization, but instead of applying this operator in the subsurface-offset or aperture angle gathers, we apply it in the time-shift gathers computed by the time-shift imaging condition. Using a simple model of two Gaussian velocity anomalies, we show the image and slowness perturbations by applying the forward and adjoint tomographic operators and image perturbation and the gradient of the objective function after differential semblance operator on time-shift gathers. The results are promising.
\end{abstract}

\section{Introduction}
Wave-equation migration velocity analysis (WEMVA)\cite[]{savaphd,sava04} promises to overcome the limitations of migration velocity analysis due to the high-frequency approximation, which consider rays as carriers of information. Instead, WEMVA uses band-limited source and receiver wavefields to compute the gradient of the objective function that back-propagates the slowness perturbation through the model space. The objective function is represented by some norm of the image perturbation and reflects the correctness of the slowness model. In WEMVA, the image perturbation is computed by applying an operator 
%Instead, WEMVA uses band-limited perturbed source and receiver wavefields and scattered source and receiver wavefields to compute the gradient of the objective function that back-propagates the slowness perturbation through the model space. The scattered wavefields are computed by applying the scattering operator \cite[]{huang:1524} on the background wavefields, while  perturbed wavefields are computed by convolving the background receiver and source wavefields, respectively, with a perturbed image. The objective function is represented by some norm of the perturbed image and reflects the correctness of the slowness model.
\par

% In the differential semblance velocity analysis (DVSA) \cite[]{shen:VE49}, the perturbed image is obtained via differential semblance optimization (DSO). DSO can be applied on both subsurface-offset gathers or aperture-angle gathers. Since there is no residual moveout  interpretation, the objective function is fully automated and the velocity inversion can be carried out by Quasi-Newton methods, like L-BFGS \cite[]{shen:VE49}.
\par
more still to come .... \cite {sava:S209} introduce the time-shift imaging condition, in which source and receiver wavefields are phase shifted prior to correlation. The kinematic properties of these image gathers have been explored to perform velocity update, after transformation to reflection-angle gathers \cite[]{liu:2797}.
\par
In this paper we investigate the use of time-shift image gathers in DVSA. We start by briefly describing the forward and adjoint WEMVA operators; then we explain the time-shift imaging condition; proceed to develop the new objective function resultant of applying DSO-like operator on those image gathers; and finalize by showing some numerical examples.
\section{Forward and adjoint WEMVA operators}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Wave-equation migration velocity analysis is a non-linear inverse problem that solves for an optimum slowness by minimizing an objective function which measures the correctness of the slowness. The objective function is represented by the perturbed image, which can be computed in different ways such as, the difference between the background image and the residual prestack migrated image using an improved background slowness \cite[]{savaphd,sava04}; and the background image subjected to the differential semblance operator, in the subsurface-offset or reflection-angle domain, that characterizes DVSA \cite[]{shen:VE49}.  
\par
Considering sufficiently good illumination, the energy that departs from the zero-subsurface offset is an indicative of inaccurate slowness. DSO in the subsurface-offset domain applies a $\left | h \right |$-weight, the absolute value of the subsurface-offset, such that  down   method \cite[]{shen2004}, the focusing operator takes the following form:
\begin{eqnarray}
{\bf F}({\bf I}) = \left( {\bf 1} - {\bf O} \right) {\bf I}, \label{eq_focus_dso}
\end{eqnarray}
where ${\bf 1}$ is the identity operator and ${\bf O}$ is the DSO operator either in the subsurface offset domain or in the angle domain \cite[]{shen2004}.
The subsurface-offset-domain DSO focuses the energy at zero offset, whereas the angle-domain DSO flattens the ADCIGs.

In the wave-equation migration velocity analysis (WEMVA) method \cite[]{sava2004}, the focusing operator is the linearized residual migration operator defined as follows: 
\begin{eqnarray}
{\bf F}({\bf I}) = {\bf R}[\rho]{\bf I} \approx {\bf I} + {\bf K}[\Delta \rho]{\bf I},
\end{eqnarray}
where ${\bf \rho}$ is the ratio between the background slowness $\widehat{\bf s}$ and the true slowness ${\bf s}$, 
and $\Delta \rho = 1-\rho = 1-\frac{\widehat{\bf s}}{\bf s}$; 
${\bf R}[\rho]$ is the residual migration operator \cite[]{sava2003}, and 
${\bf K}[\Delta \rho]$ is the differential residual migration operator defined as follows \cite[]{sava2004a,sava2004b}:
\begin{eqnarray}
{\bf K}[\Delta \rho] = \Delta \rho \left. \frac{\partial {\bf R}[\rho]}{\partial \rho} \right|_{\rho=1}.
\end{eqnarray}
The linear operator ${\bf K}[\Delta \rho]$ applies different phase rotations to the image for different reflection angles and geological dips \cite[]{biondi2008}.

In general, if we choose $\ell_2$ norm, the tomography objective function to minimize can be written as follows:
\begin{eqnarray}
J = \frac{1}{2}||{\bf \Delta I}||_2 = \frac{1}{2}||{\bf I}-{\bf F}({\bf I})||_2,
\end{eqnarray}
where $||\cdot||_2$ stands for the $\ell_2$ norm.
Gradient-based optimization techniques such as the quasi-Newton method and the conjugate gradient method can be used to minimize the objective function $J$.
The gradient of $J$ with respect to the slowness ${\bf s}$ reads as follows:
\begin{eqnarray}
{\bf \nabla} J = \Re \left( \left(\frac{\partial {\bf I}}{\partial {\bf s}} - \frac{\partial {\bf F}({\bf I})}{\partial {\bf s}} \right)^{\prime} \left( {\bf I}-{\bf F}({\bf I})\right) \right), \label{eq_obj_grad_gen}
\end{eqnarray}
where $\Re$ denotes taking the real part of a complex value and $^{\prime}$ denotes the adjoint. 
For the DSO method, the linear operator ${\bf O}$ is independent of the slowness, 
so we have
\begin{eqnarray}
\frac{\partial {\bf F}({\bf I})}{\partial {\bf s}} = ({\bf 1}-{\bf O})\frac{\partial {\bf I}}{\partial {\bf s}}. \label{eq_focus_grad_dso}
\end{eqnarray}
Substituting Equations \ref{eq_focus_dso} and \ref{eq_focus_grad_dso} into Equation \ref{eq_obj_grad_gen} and evaluating the gradient at a background slowness yields
\begin{eqnarray}
{\bf \nabla} J_{\rm DSO} = \Re \left( \left.\left(\frac{\partial {\bf I}}{\partial {\bf s}}\right|_{{\bf s}=\widehat{\bf s}}\right)^{\prime}{\bf O}^{\prime}{\bf O}\widehat{\bf I}\right), \label{eq_grad_dso}
\end{eqnarray}
where $\widehat {\bf I}$ is the background image computed using the background slowness $\widehat{\bf s}$.

For the WEMVA method, the gradient is slightly more complicated, because in this case, the focusing operator is also dependent 
on the slowness ${\bf s}$. However, one can simplify it by assuming that the focusing operator is applied on the background
image $\widehat{\bf I}$ instead of ${\bf I}$, and $\widehat{\Delta \rho}$ is also picked from the background image $\widehat{\bf I}$, that is 
\begin{eqnarray}
{\bf F}(\widehat{\bf I})=\widehat{\bf I}+{\bf K}[\widehat{\Delta \rho}]\widehat{\bf I}.
\end{eqnarray}
With these assumptions, we get the "classic" WEMVA gradient as follows:
\begin{eqnarray}
\nabla J_{\rm WEMVA} = \Re \left( - \left.\left(\frac{\partial {\bf I}}{\partial {\bf s}}\right|_{{\bf s}=\widehat{\bf s}}\right)^{\prime}{\bf K}[\widehat{\Delta \rho}]\widehat{\bf I} \right)
\end{eqnarray}
The complete WEMVA gradient without the above assumptions can also be derived following the method described by \cite{biondi2008}. 

No matter which gradient we choose to back-project the slowness perturbation, we have to evaluate the adjoint of the linear operator
$\left.\frac{\partial{\bf I}}{\partial{\bf s}}\right|_{{\bf s}=\widehat{\bf s}}$, which defines a linear mapping from
the slowness perturbation ${\bf \Delta s}$ to the image perturbation ${\bf \Delta I}$. This is easy to see by
expanding the image ${\bf I}$ around the background slowness ${\widehat {\bf s}}$ as follows:
\begin{eqnarray}
{\bf I} = {\widehat {\bf I}} + \left. \frac{\partial{\bf I}}{{\partial{\bf s}}}\right|_{{\bf s}=\widehat{\bf s}}({\bf s}-{\widehat {\bf s}}) + \cdots,
\end{eqnarray} 
Keeping only the zero and first order terms, we get the linear operator $\left. \frac{\partial{\bf I}}{\partial{\bf s}}\right|_{{\bf s}=\widehat{\bf s}}$ as follows:
\begin{eqnarray}
{\bf \Delta I} = \left. \frac{\partial{\bf I}}{\partial{\bf s}}\right|_{{\bf s}=\widehat{\bf s}} {\bf \Delta s} = {\bf T}{\bf \Delta s},
\end{eqnarray}
where ${\bf \Delta I} = {\bf I} - {\widehat {\bf I}}$ and ${\bf \Delta s} = {\bf s} - {\widehat {\bf s}}$. 
${\bf T} = \left. \frac{\partial{\bf I}}{\partial{\bf s}}\right|_{{\bf s}=\widehat{\bf s}}$ is the wave-equation tomographic operator.
The tomographic operator can be evaluated either in the source and receiver domain \cite[]{sava2004} or in the shot-profile domain \cite[]{shen2004}.
In next section we follow an approach similar to that discussed by \cite{shen2004} and review the forward and adjoint tomographic operator 
in the shot-profile domain. In the subsequent sections, we generalize the expression of the tomographic operator to 
the generalized source domain.

\section{the tomographic operator in the shot-profile domain}

For the conventional shot-profile migration, both source and receiver wavefields are downward continued with the 
following one-way wave equations \cite[]{claerbout1971}:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right)D({\bf x},{\bf x}_s,\omega) = 0 \\
D(x,y,z=0,{\bf x}_s,\omega) = \overline{f_s(\omega)\delta({\bf x}-{\bf x}_s)} \end{array} \right., \label{eq_wveq_bD}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right)U({\bf x},{\bf x}_s,\omega) = 0 \\
U(x,y,z=0,{\bf x}_s,\omega) = Q(x,y,z=0,{\bf x}_s,\omega) \end{array} \right., \label{eq_wveq_bU}
\end{eqnarray}
where the overline stands for complex conjugate; 
$D({\bf x},{\bf x}_s,\omega)$ is the source wavefield for a single frequency $\omega$ at image point ${\bf x}=(x,y,z)$ 
with the source located at ${\bf x}_s=(x_s,y_s,0)$; $U({\bf x},{\bf x}_s,\omega)$ is the 
receiver wavefield for a single frequency $\omega$ at image point ${\bf x}$ for the source located at ${\bf x}_s$;
$s({\bf x})$ is the slowness at ${\bf x}$; ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector; $f_s(\omega)$
is the frequency dependent source signature, and $\overline{f_s(\omega)\delta({\bf x}-{\bf x}_s)}$ defines the point source function at ${\bf x}_s$, 
which serves as the boundary condition of Equation \ref{eq_wveq_bD}. $Q(x,y,z=0,{\bf x}_s,\omega)$ is the recorded shot gather for the shot located at ${\bf x}_s$, 
which serves as the boundary condition of Equation \ref{eq_wveq_bU}.
To produce the image, the following cross-correlation imaging condition is used:
\begin{eqnarray}
I({\bf x},{\bf h}) = \sum_{{\bf x}_s}\sum_{\omega} D({\bf x}-{\bf h},{\bf x}_s,\omega) U({\bf x}+{\bf h},{\bf x}_s,\omega), \label{eq_ic_bI}
\end{eqnarray}
where ${\bf h}=(h_x,h_y,h_z)$ is the subsurface half offset.

The perturbed image can be derived by a simple application of the chain rule
to Equation \ref{eq_ic_bI}: 
\begin{eqnarray}
\Delta I({\bf x},{\bf h})  &=& 
\sum_{{\bf x}_s}\sum_{\omega} 
\left( \Delta D({\bf x}-{\bf h},{\bf x}_s,\omega) {\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega) + \right. \nonumber \\
& & \left. {\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega) \Delta U ({\bf x}+{\bf h},{\bf x}_s,\omega) \right), \label{eq_ic_dI}
\end{eqnarray}
where ${\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)$ and ${\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega)$ are the background source and 
receiver wavefields computed with the background slowness ${\widehat s}({\bf x})$; 
$\Delta D({\bf x}-{\bf h},{\bf x}_s,\omega)$ and $\Delta U({\bf x}+{\bf h},{\bf x}_s,\omega)$ 
are the perturbed source wavefield and perturbed receiver wavefield, which are the results of the slowness perturbation $\Delta s({\bf x})$.
The perturbed source and receiver wavefields satisfy the following 
one-way wave equations, which are linearized with respect to slowness (see Appendix A for derivations):
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {\widehat s}^2 ({\bf x})-|{\bf k}|^2}\right) \Delta D({\bf x},{\bf x}_s,\omega)= 
\frac{- i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}^2({\bf x})}}}{\widehat D}({\bf x},{\bf x}_s,\omega) \\
\Delta D(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right. , \label{eq_wveq_dD}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {\widehat s}^2 ({\bf x})-|{\bf k}|^2}\right) \Delta U({\bf x},{\bf x}_s,\omega)=
\frac{- i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}^2({\bf x})}}}{\widehat U}({\bf x},{\bf x}_s,\omega) \\
\Delta U(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right.. \label{eq_wveq_dU}
\end{eqnarray}
Recursively solving Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} gives us the perturbed source and receiver wavefields.
The perturbed source and receiver wavefields are then used in Equation \ref{eq_ic_dI} to generate the perturbed image $\Delta I({\bf x},{\bf h})$, where
the background source and receiver wavefields are precomputed by recursively solving Equations \ref{eq_wveq_bD} and \ref{eq_wveq_bU} with a
background slowness ${\widehat s}(\bf x)$.
Appendix B gives a more detailed matrix representation of how to evaluate the forward tomographic operator ${\bf T}$.

To evaluate the adjoint tomographic operator ${\bf T}'$, we first apply the adjoint of the imaging condition in Equation \ref{eq_ic_dI} 
to get the perturbed source and receiver wavefields $\Delta D({\bf x},{\bf x}_s,\omega)$ and $\Delta U({\bf x},{\bf x}_s,\omega)$ as follows:
\begin{eqnarray}
\Delta D({\bf x},{\bf x}_s,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h})\overline{{\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega)}, \\
\Delta U({\bf x},{\bf x}_s,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h})\overline{{\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)}.
\end{eqnarray}
Then we solve the adjoint equations of Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} to get the slowness perturbation $\Delta s({\bf x})$. 
Again, in order to solve the adjoint equations of Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU}, 
the background source wavefield ${\widehat D}({\bf x},{\bf x}_s,{\omega})$ and 
the background receiver wavefield ${\widehat U}({\bf x},{\bf x}_s,\omega)$ have to be computed in advance.
Appendix C gives a more detailed matrix representation of how to evaluate the adjoint tomographic operator ${\bf T}^{'}$.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Time-shift imaging condition}
Based on the theory of reflector mapping \cite[]{claerbout71}, the image is formed at points where source and receiver wavefields coincide both in time of wavefield propagation and spatially. In the context of migration by downward-continuation, this translates into point-wise multiplying the receiver wavefield, $R$, with the complex conjugate of the source wavefield, $S^{\ast}$, at every image location, $\bf{x}$, and summing over all frequencies, $\omega$, which represents to extract the image, $I$,at zero time. For a single shot we have
\begin{equation}
I(\bf{x}) = \sum_{\omega} \it{S}^{\ast}(\bf{x},\omega) \it{R}(\bf{x},\omega).
\label{eq:e1}
\end{equation}
\par
This imaging condition, however, does not yield an image suited for estimating velocity errors. To overcome the limitation of the above imaging condition, different generalizations have been proposed \cite[]{rickett,biosymes,sava:2052}. They consist of spatially shifting the wavefields with respect to each other prior to the point-wise multiplication. The spatial shifts are interpreted as subsurface offsets, $\bf{h}$, which can be computed in any direction. For the case of complete acquisition aperture and correct velocity, negligible amplitudes are computed for subsurface offsets different from zero. Departures from the zero offset are  indicative of illumination problems or velocity errors. In general, the image is computed by
\begin{equation}
I(\bf{x},\bf{h}) = \sum_{\omega} \it{S}^{\ast}(\bf{x}-\bf{h},\omega) \it{R}(\bf{x}+\bf{h},\omega).
\label{eq:e2}
\end{equation}
\par
\cite {sava:S209} introduced the space- and time-shift imaging condition, from which the time-shift imaging condition is a special case. It is computed by applying a phase-shift to the wavefields prior to summation over frequencies
\begin{equation}
I(\bf{x},\tau) = \sum_{\omega} \it{S}^{\ast}(\bf{x},\omega) \it{R}(\bf{x},\omega)e^{i2\omega\tau}.
\label{eq:e3}
\end{equation}
\par
Time-shift image gathers present linear events with slope equal to the local velocity at the reflector. If migration is performed with a faster velocity, energy concentrates at negative time-shifts and shallower than the correct depth, indicating that source and receiver rays cross above the imaging point (where $\tau = 0$). If migration is performed with a slower velocity, energy concentrates at positive time-shifts and deeper than the correct depth, indicating that rays cross below the imaging point. For the case of correct velocity, energy concentrates at  $\tau = 0$.

\section{Using time-shift imaging condition for velocity update}
\cite {sava:S209} show that time-shift image gathers can be transformed into angle-domain common image gathers by slant-stacking. In 2D, the transformation can be carried out by a Radon transform
\begin{equation}
\it{I}(z,x,\gamma) = \int \it{I}(z,x,\tau) d\tau \arrowvert_{\zeta(\gamma,\alpha,\tau) },
\label{eq:e4}
\end{equation}
where $\gamma$ is apparent the reflection angle, $\alpha$ is the apparent geological dip, and
\begin{equation}
\zeta(\gamma,\alpha,\tau) = z + v\tau \sec \alpha \sec \gamma
\label{eq:e5}
\end{equation}
is the stacking path.
\par
This transformation yields images suited for velocity update. Using angle gathers, in wave-equation migration velocity analysis 
\cite[]{sava05}, a residual moveout parameter is picked and used to compute the image perturbation. \cite{shen:VE49} use angle gathers, computed from subsurface offset gathers, to obtain the perturbed image by applying the differential semblance operator (DSO). They show that the results of velocity inversion using DSO in the angle domain ($\bf{D}_{\gamma}$) to compute the image perturbation are of lower quality than that obtained using DSO in the offset domain ($\bf{D}_{h}$). They claim that the different performances are due to the fact that the Hessian is better conditioned in the case of the offset domain DSO.
\par
To apply a DSO-like strategy when using time-shift image gathers, transformation to the angle domain requires information about the geological dip. If this information is not available or not reliable, the slant-stack transformation should include the apparent geological dip dimension, turning the computation of the perturbed image more expensive. So, it is desirable to derive an offset-domain DSO-like operator to be used in the time-shift gathers.
\par
The equivalence between the time-shift DSO ($\bf{D}_{\tau}$) with $\bf{D}_{\gamma}$, if exists, can be expressed by
\begin{equation}
\bf{D}_{\tau} \it{I}(\bf{x},\tau) = \bf{R}_{\gamma}^{-1}  \bf{D}_{\gamma} \bf{R}_{\gamma} \it{I}(\bf{x},\tau),
\label{eq:e6}
\end{equation}
where $\bf{R}_{\gamma}$ is the time-shift to angle transformation (equation \ref{eq:e4}). The perturbed image, after application of forward and adjoint operators on the background image, is
\begin{equation}
\bf{D}_{\tau}^{\ast} \bf{D}_{\tau} \it{I}(\bf{x},\tau) = \bf{R}_{\gamma}^{\ast} \bf{D}_{\gamma}^{\ast} \bf{D}_{\gamma} \bf{R}_{\gamma} \it{I}(\bf{x},\tau).
\label{eq:e7}
\end{equation}
\par
Recalling that $\bf{D}_{\gamma} = \frac{\partial}{\partial q}$, where $q = \sec \gamma $, and that equation \ref{eq:e5} can be written as
\begin{eqnarray}
\bf{R}_{\gamma} \it{I}(z,x,\tau)& = & \int \it{I}(z + \it{qv}\tau\sec\alpha,x,\tau) d\tau \nonumber \\
 & = & \int \int \it{I}(\varsigma,x,\tau) \delta(\varsigma - z - \it{qv}\tau\sec\alpha) d\varsigma d\tau \nonumber \\
 & = & \int \int \int I(\varsigma,x,\tau) e^{ik_{\varsigma}(\varsigma - z - \it{qv}\tau\sec\alpha)} d\varsigma d\tau dk_{\varsigma},
\label{eq:e8}
\end{eqnarray}
we have
\begin{eqnarray}
\bf{D}_{\tau}^{\ast} \bf{D}_{\tau} \it{I}(z,x,\tau) & = & \bf{R}_{\gamma}^{\ast} \int \int \int{-\frac{\partial^2}{\partial q^2} 
\it{I}(\varsigma,x,\tau) e^{ik_{\varsigma}(\varsigma - z - \it{qv}\tau\sec\alpha)} d\varsigma d\tau dk_{\varsigma}} \nonumber \\
 & = & \bf{R}_{\gamma}^{\ast} \int \int \int{(\it{v}\tau \sec \alpha \it{k}_{\varsigma})^2 \it{I}(\varsigma,x,\tau) e^{ik_{\varsigma}(\varsigma - z - \it{qv}\tau\sec\alpha)} d\varsigma d\tau dk_{\varsigma}}.
\label{eq:e9}
\end{eqnarray}
\par
By moving the adjoint of the Radon transform inside the integral, and writing $\xi = \it{v} \sec \alpha$, we obtain
\begin{eqnarray}
\bf{D}_{\tau}^{\ast} \bf{D}_{\tau} \it{I}(z,x,\tau) & = & \int \int \left [ \int \int \int (\xi\tau \it{k}_{\varsigma})^2 \it{I}(\varsigma',x,\tau') e^{ik'_{\varsigma}(\varsigma' - z - \xi\it{q}\tau')} d\varsigma' d\tau' dk'_{\varsigma} \right ] e^{-ik_{\varsigma}(\varsigma - z - \xi\it{q}\tau)} dq dk_{\varsigma} \nonumber \\
 & = & \int \int \int \int (\xi\tau \it{k}_{\varsigma})^2 \it{I}(\varsigma',x,\tau') e^{ik_{\varsigma}(\varsigma' - \varsigma)} e^{ik_{\varsigma} \it{q}\xi(\tau' - \tau)} dq d\tau' d\varsigma' dk_{\varsigma} \nonumber \\
 & = & \int \int \int (\xi\tau\it{k}_{\varsigma})^2 \it{I}(\varsigma',x,\tau') e^{ik_{\varsigma}(\varsigma' - \varsigma)} \left [ \int e^{ik_{\varsigma}\xi q(\tau' - \tau)} dq \right ] d\tau' d\varsigma' dk_{\varsigma}.
\label{eq:e10}
\end{eqnarray}
\par
In the bracketed term, one can recognize the Fourier transform of $\frac{\delta(\tau'- \tau)}{\left [k_{\varsigma}\xi \right ]}$. So, equation \ref{eq:e10} reads
\begin{eqnarray}
\bf{D}_{\tau}^{\ast} \bf{D}_{\tau} \it{I}(z,x,\tau) = \int \int \tau^2 \left|\xi\it{k}_{\varsigma} \right | \it{I}(\varsigma',x,\tau) e^{ik_{\varsigma}(\varsigma' - \varsigma)} d\varsigma' dk_{\varsigma}.
\label{eq:e11}
\end{eqnarray}
\par
Finally, using the relation 
\begin{equation}
\int \left| k \right | e^{ikx} dk = \frac{-2}{x^2} \nonumber
\end{equation}
we get
\begin{eqnarray}
\bf{D}_{\tau}^{\ast} \bf{D}_{\tau} \it{I}(z,x,\tau) & = & -2 \tau^2 \it{v} \left|\sec \alpha \right | \int {\frac{\it{I}(\varsigma',x,\tau)}{(\varsigma' - \varsigma)^2} d\varsigma'}.
\label{eq:e12}
\end{eqnarray}
\par
Equation \ref{eq:e12} shows that the operator $\bf{D}_{\tau}$ is similar to $\bf{D}_{h}$, except for the scaling by the absolute value of the secant of the apparent geological dip and local velocity. The perturbed image can be computed by applying the $\left | \tau \right |$ operator. However, the spread of energy along the line with slope equal to the local velocity poses a problem to $\bf{D}_{\tau}$. Even if it is used the correct velocity, the perturbed image will present strong amplitudes away from the zero time.
\par
To circumvent this problem, we use a simple thresholding strategy prior to $\bf{D}_{\tau}$, such that the spread energy is zeroed out.
more still to come ...

\section {Numerical examples}
To test the potential of using time-shift image gathers in wave-equation migration velocity analysis, we use a simple model of two Gaussian anomalies located at (-800,800) and (800,800) embedded in a background velocity of 2000m/s. A reflector is positioned at a depth of 1500m. The Gaussian anomaly at negative x position has 10$\%$ slower velocity than the background and the Gaussian anomaly located at positive x position has 10$\%$ faster velocity than the background. We used one-way Born modeling to generate 81 shots with 100m spacing, maximum offset of 6000m, and receiver interval equals to 10m. Figure \ref{fig:f0} shows the image computed with a background velocity of 2000m/s and time-shift imaging condition.
\plot{f0}{width=\textwidth,height=0.225\textheight}{Background image computed time-shift imaging condition.$\bf{[CR]}$}
\par
Using a background velocity of 2000m/s, we computed the perturbed image by applying the forward tomographic operator in the shot-profile domain. In this case, the slowness perturbation (Figure \ref{fig:f1}) is known beforehand. Figure \ref{fig:f2} shows the perturbed image.
\plot{f1}{width=\textwidth,height=0.225\textheight}{Slowness perturbation used to compute the image perturbations of Figures \ref{fig:f2} and \ref{fig:f4}.$\bf{[ER]}$}
\plot{f2}{width=\textwidth,height=0.225\textheight}{Image perturbation computed with the slowness perturbation of Figure \ref{fig:f1} and time-shift imaging condition.$\bf{[CR]}$}
\par
The perturbed image of Figure \ref{fig:f2} was used to compute the slowness perturbation by applying the adjoint of the tomographic operator. The predicted slowness perturbation is shown in Figure \ref{fig:f3}. Strong artifacts occur due to the poor shot sampling. However, it is noticeable that the correct signal of the slowness perturbation is recovered.
\plot{f3}{width=\textwidth,height=0.225\textheight}{Predicted slowness perturbation computed with the image perturbation of Figure \ref{fig:f2} and time-shift imaging condition.$\bf{[CR]}$}
\par
For comparison, we show the perturbed image (Figure \ref{fig:f4}) and the the slowness perturbation (Figure \ref{fig:f5}) computed for the case of space-shift imaging condition. The clip of the amplitudes of the space-shift slowness perturbation is different from the time-shift slowness perturbation to highlight its general structure. Overall, both slowness perturbations show a similar pattern. However, the effect of poor shot sampling is more pronounced in the time-shift slowness perturbation than in the space-shift slowness perturbation.
\plot{f4}{width=\textwidth,height=0.225\textheight}{Image perturbation computed with the slowness perturbation of Figure \ref{fig:f1} and using space-shift imaging condition.$\bf{[CR]}$}
\plot{f5}{width=\textwidth,height=0.225\textheight}{Predicted slowness perturbation computed with the space-shift image perturbation of Figure \ref{fig:f4} and space-shift imaging condition.$\bf{[CR]}$}

\section{Conclusions}

\bibliographystyle{seg}
\bibliography{SEP,ref137a}
