\def\figdir{Fig}
%\documentclass{segabs}
%\documentclass[paper,twocolumn,twoside]{geophysics}
%\documentclass[manuscript,revised]{geophysics}
%\newcommand{\rs}[1]{\mathstrut\mbox{\scriptsize\rm #1}}
%\newcommand{\rr}[1]{\mbox{\rm #1}}
%\usepackage{amsmath}
%\usepackage{multicol}
%\usepackage{placeins}
% % figure directory
%\begin{document}
%\renewcommand{\thefootnote}{\fnsymbol{footnote}} 
 

\lefthead{Guerra et al.}
\righthead{Image-space wave-equation tomography}
%\email{\footnotemark[1]guerrac@stanford.edu \footnotemark[2]tang@stanford.edu \footnotemark[3]biondo@stanford.edu}
\footer{SEP--138}
\title{Wave-equation tomography using image-space phase-encoded data}
\author{Claudio Guerra, Yaxun Tang, and Biondo Biondi}

\maketitle

\begin{abstract}
Wave-equation tomography in the image-space is a powerful technique that promises to yield more reliable velocity models than ray-based migration velocity analysis in areas of complex overburden. Its practical use, however, has been limited because of the high computational cost. Applying a target-oriented approach and using data reduction can make wave-equation tomography in the image space of practical use. Here, we present results of applying image-space wave-equation tomography in the generalized source domain, where a small number of synthesized shot gathers are generated. Specifically, we generate synthesized shot gathers by image-space phase encoding. This technique can also be used in a target-oriented way. The comparison of the gradients of the tomography objective functional obtained using image-space encoded gathers with those obtained using the original shot gathers shows that those encoded shot gathers can be used in wave-equation tomography problems. Velocity inversion using image-space phase-encoded gathers converges to reasonable results when compared to the correct velocity model. We illustrate our method by applying it to the Marmousi model.
\end{abstract}
\section{introduction}
Wave-equation tomography has the potential to overcome the problems faced by ray-based traveltime tomography when estimating the velocity model in complex geological scenarios. This is because wave-equation tomography uses band-limited wavefields instead of infinite-frequency rays as carriers of information; thus it is robust even in the presence of strong velocity contrasts and immune to multi-pathing issues. However, despite its theoretical advantages, wave-equation tomography is still computationally challenging. 
\par
Wave-equation tomography can be performed in the data-space domain \cite[]{tarantola1987,woodward1992} or in the image-space domain \cite[]{biondi1999,shen2004}. The image-space approach minimizes the residual in the image domain obtained after migration. Regardless of the domain of application, using phase-encoded data can substantially decrease the computational cost of wave-equation tomography\cite[]{vigh:VE135,shen:VE49}. \cite{yax2008} extended the theory of image-space wave-equation tomography from the conventional shot-profile domain \cite[]{shen2004} to the generalized source domain. The generalized source domain can be obtained in two different spaces. In the data-space, shot gathers are combined and the corresponding source function is synthesized, using a convenient phase-encoding scheme, which characterizes the data-space phase encoding \cite[]{whitmore1995,romero2000}. In the image-space, source- and receiver-areal data are synthesized by upward propagating wavefields. The initial condition for the modeling is a prestack image computed with wave-equation migration, according to the prestack exploding-reflector modeling \cite[]{biondi2006,biondi2007}. The modeling experiments can be combined such that a small quantity of areal data is generated. In this case, to mitigate crosstalk during imaging, the modeling experiments and reflectors are phase-encoded, characterizing the image-space phase encoding \cite[]{guerra2008a}. To encode the reflectors, a picking step of some key reflectors is necessary.

In this paper, we show that image-space phase-encoded wavefields can be used to estimate the velocity model in image-space wave-equation tomography. We show that the gradient of the tomographic objective functional is similar to that obtained in the original shot-profile domain, but with less computational cost. Velocity inversion using image-space phase-encoded gathers converges to reasonable results when compared to the correct velocity model, provided that crosstalk has been sufficiently attenuated. We briefly discuss the theory of wave-equation tomography in the image-space domain; then we explain the prestack exploding-reflector modeling and show that the image-space phase encoding can be used to accelerate wave-equation tomography in the image domain. We use the Marmousi model to illustrate the method.

\section{image-space wave-equation tomography}
Image-space wave-equation tomography is a non-linear inverse problem that tries to find an optimal background slowness that minimizes the residual field, $\Delta {\bf I}$, defined in the image space. The residual field is derived from the background image, $\bf I$, which is computed with a background slowness. The general form of the residual field is \cite[]{biondi2008}
\begin{eqnarray}
{\Delta {\bf I}} = {\bf I} - {\bf F}({\bf I}), \label{eq_focus}
\end{eqnarray}                  
where $\bf F$ is a focusing operator, which measures the focusing of the migrated image. In particular, in the Differential Semblance Optimization (DSO) method \cite[]{shen2004}, the focusing operator takes the 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}.
\par
Under $\ell_2$ norm, the tomography objective function can be written as follows:
\begin{eqnarray}
J = \frac{1}{2}||{\Delta {\bf I}}||_2 = \frac{1}{2}||{\bf I}-{\bf F}({\bf I})||^2.
\end{eqnarray}
The gradient of $J$ with respect to the slowness ${\bf s}$ is
\begin{eqnarray}
{\bf \nabla} J = \left(\frac{\partial {\bf I}}{\partial {\bf s}} - \frac{\partial {\bf F}({\bf I})}{\partial {\bf s}} \right)^{\ast} \left( {\bf I}-{\bf F}({\bf I})\right), \label{eq_obj_grad_gen}
\end{eqnarray}
where $^{\ast}$ denotes the adjoint. 

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 ${\Delta {\bf s}}$ to the image perturbation ${\Delta {\bf I}}$, can be computed by expanding the image ${\bf I}$ around the background slowness ${\widehat {\bf s}}$. Keeping only the zeroth 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}
{\Delta {\bf I}} = \left. \frac{\partial{\bf I}}{\partial{\bf s}}\right|_{{\bf s}=\widehat{\bf s}} {\Delta {\bf s}} = {\bf T}{\Delta {\bf s}},
\end{eqnarray}
where ${\Delta {\bf I}} = {\bf I} - {\widehat {\bf I}}$, ${\widehat {\bf I}}$ is the background image computed with the background slowness ${\widehat {\bf s}}$ and ${\Delta {\bf 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 the shot-profile domain, both source and receiver wavefields are downward continued with the one-way wave equations \cite[]{claerbout1971}
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\Lambda\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\Lambda\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$; $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}. $\Lambda$ is the square-root operator
\begin{eqnarray}
\begin{array}{l}
\Lambda = \sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2},
\end{array}
\label{lambda}
\end{eqnarray}
where $s({\bf x})$ is the slowness at ${\bf x}$; ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector. The image is computed by applying the cross-correlation imaging condition:
\begin{eqnarray}
I({\bf x},{\bf h}) = \sum_{{\bf x}_s}\sum_{\omega} \overline{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 $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 and ${\bf h}=(h_x,h_y,h_z)$ is the subsurface half-offset.

%\begin{eqnarray}
%\left\{ \begin{array}{l}
%\left( \frac{\partial}{\partial z}+i\Lambda\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\Lambda\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$; $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}. 

The perturbed image can be derived by the 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( \overline{\Delta D({\bf x}-{\bf h},{\bf x}_s,\omega)} {\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega) + \right. \nonumber \\
& & \left. \overline{{\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})$. 
\par
To evaluate the adjoint of the tomographic operator, ${\bf T}^{\ast}$, we first apply the adjoint of the imaging condition to get the perturbed source and receiver wavefields, $\Delta D ({\bf x}+{\bf h},{\bf x}_s,\omega)$ and $\Delta U ({\bf x}+{\bf h},{\bf x}_s,\omega)$, as follows
\par
\begin{eqnarray}
\Delta D({\bf x},{\bf x}_s,\omega) &=& \sum_{{\bf x}_s}\sum_{\omega} \Delta I({\bf x},{\bf h}) {\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega) \nonumber \\
\Delta U({\bf x},{\bf x}_s,\omega) &=& \sum_{{\bf x}_s}\sum_{\omega} \Delta I({\bf x},{\bf h}) {\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega). \label{eq_aic_dI}
\end{eqnarray}
\par
The perturbed source and receiver wavefields satisfy the following one-way wave equations,  linearized with respect to slowness:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\Lambda\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\Lambda\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}
\par
When solving the optimization problem, the gradient of the objective function is obtained by computing the perturbed wavefields (equation \ref{eq_aic_dI}), where the image perturbation results from the application of a focusing operator (equation \ref{eq_focus}) on the background image; then the adjoint of the one-way wave equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} are solved by upward propagating the perturbed wavefields and cross-correlating  then with the scattered wavefields (right-hand side of equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} ); and, finally, summing the cross-correlation results. Figure \ref{fig:ISWET} displays the image-space wave-equation tomography flowchart. The upper gray box represents the process of obtaining the image perturbation, while the lower gray box corresponds to the application of the adjoint of the wave-equation tomography operator.    More detailed information on how to evaluate the forward and adjoint operators can be found in \cite {yax2008}.
\activeplot{ISWAT}{width=\columnwidth}{}{Image-space wave-equation tomography flowchart. The upper gray box represents the process of obtaining the image perturbation, while the lower gray box corresponds to the application of the adjoint of the wave-equation tomography operator.}

%\begin{eqnarray}
%\left\{ \begin{array}{l}
%\left( \frac{\partial}{\partial z}+i\Lambda\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\Lambda\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}
%The perturbed source and receiver wavefields are 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 need to be precomputed.

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

%\section{Data-space encoded wavefields}
%The data-space encoded shot gathers are obtained by linear combination of the original shot gathers after phase encoding with respect to the parameter ${\bf p}_s$. ${\bf p}_s$ is the ray parameter, in the case of plane-wave phase encoding, or the index of different random realizations, in the case of random phase encoding. Because of the linearity of the one-way wave equation with respect to the wavefield, the encoded source and receiver wavefields satisfy the one-way wave equation. The final image is obtained by applying the cross-correlation imaging condition and summing the images for all ${\bf p}_s$'s:
%\begin{eqnarray}
%I_{\rm de}({\bf x},{\bf h}) = \sum_{{\bf p}_s,{\omega}} |c|^2{\widetilde D}({\bf x}-{\bf h},{\bf p}_s,\omega) {\widetilde U}({\bf x}+{\bf h},{\bf p}_s,\omega), 
%\label{eq_ic_dI_de}
%\end{eqnarray}
%where $\bf h$ is the subsurface offset, and $c=\omega$ for plane-wave phase encoding and $c=1$ for random phase encoding \cite[]{yaxun2008b}. 

%From Equation \ref{eq_ic_dI_de}, the perturbed image can be obtained as follows:
%\begin{eqnarray}
%\Delta I_{\rm de}({\bf x},{\bf h})  &=& 
%\sum_{{\bf p}_s,{\omega}}|c|^2 
%\left( \Delta \widetilde{D}({\bf x}-{\bf h},{\bf p}_s,\omega) \widehat{\widetilde{U}}({\bf x}+{\bf h},{\bf p}_s,\omega) + \right. \nonumber \\
%& & \left. \widehat{\widetilde{D}}({\bf x}-{\bf h},{\bf p}_s,\omega) \Delta \widetilde{U}({\bf x}+{\bf h},{\bf p}_s,\omega) \right),
%\end{eqnarray}
%where $\widehat{\widetilde{D}}({\bf x},{\bf p}_s,\omega)$ and $\widehat{\widetilde{U}}({\bf x},{\bf p}_s,\omega)$ are the data-space encoded background source and receiver wavefields; $\Delta \widetilde{D}({\bf x},{\bf p}_s,\omega)$ and $\Delta \widetilde{U}({\bf x},{\bf p}_s,\omega)$ are the perturbed source and receiver wavefields in the data-space phase-encoding domain, which satisfy the perturbed one-way wave equation \cite[]{shen2004}. The slowness perturbation is then computed by applying the adjoint of the tomographic operator, ${\bf T} ^\ast$, to the image perturbation.

\section{prestack exploding-reflector modeling}
The general idea of prestack exploding-reflector modeling \cite[]{biondi2006} is to model the data and corresponding source function that are related to only one event in the subsurface. In this case, a single unfocused subsurface-offset-domain common-image gather (SODCIG) containing a single reflector is used as the initial condition for recursive upward continuation with the following one-way wave equations: 
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z} - i \Lambda\right) Q_D({\bf x},\omega;x_m,y_m) = I_D({\bf x},{\bf h};x_m,y_m) \\
Q_D(x,y,z=z_{\rm max},\omega;x_m,y_m) = 0 \end{array} \right.,
\label{eq_perm_bD}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z} - i \Lambda\right) Q_U({\bf x},\omega;x_m,y_m) = I_U({\bf x},{\bf h};x_m,y_m) \\
Q_U(x,y,z=z_{\rm max},\omega;x_m,y_m) = 0 \end{array} \right.,
\label{eq_perm_bU}
\end{eqnarray}
where $I_D({\bf x},{\bf h};x_m,y_m)$ and $I_U({\bf x},{\bf h};x_m,y_m)$ are the isolated SODCIGs at the horizontal location $(x_m,y_m)$ for a single reflector, and are suitable for the initial conditions for the source and receiver wavefields, respectively. As \cite{biondi2006} discusses, a rotation of the image gathers according to the apparent geological dip must be performed prior to modeling. By collecting the wavefields at the surface, we obtain the areal source data $Q_D(x,y,z=0,\omega;x_m,y_m)$ and the areal receiver data $Q_U(x,y,z=0,\omega;x_m,y_m)$ for a single reflector and a single SODCIG located at $(x_m,y_m)$. $\Lambda$ is the square-root operator defined by
\begin{eqnarray}
\begin{array}{l}
\Lambda = \sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}, \nonumber
\end{array}
\end{eqnarray}
where $s({\bf x})$ is the slowness at ${\bf x}$ and ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector.

Since the size of the migrated image volume can be very large in practice, and there are usually many reflectors in the subsurface, modeling each reflector and each SODCIG may generate a data set even larger than the original data set. One strategy to reduce the cost is to model several reflectors and several SODCIGs simultaneously \cite[]{biondi2006}; however, this process generates unwanted crosstalk. As discussed by \cite{guerra2008a}, random phase encoding can be used to attenuate the crosstalk. 
\par
One important characteristic of the prestack exploding reflector modeling is that, for velocity model building, the wavefields can be upward propagated to a certain depth level or depth horizon, provided that the velocity model above is sufficiently accurate. Therefore, a target-oriented strategy can be applied to derive the velocity model below the that depth.
%Next, we follow an approach similar to that discussed by \cite{shen2004}, in the shot profile domain, to generalize the expression of the tomographic operator to the generalized source domain using image-space phase encoding.

\section{image-space phase-encoded wavefields}
 The randomly encoded areal source and areal receiver wavefields can be computed as follows:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i \Lambda \right) {\widetilde Q}_D({\bf x},{\bf p}_m,\omega) 
= \widetilde{I}_D({\bf x},{\bf h},{\bf p}_m,\omega) \\
{\widetilde Q}_D(x,y,z=z_{\rm max},{\bf p}_m,\omega) = 0 \end{array} \right.,
\label{eq_perm_bD2}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\Lambda\right) {\widetilde Q}_U({\bf x},{\bf p}_m,\omega) 
= \widetilde{I}_U({\bf x},{\bf h},{\bf p}_m,\omega) \\
{\widetilde Q}_U(x,y,z=z_{\rm max},{\bf p}_m,\omega) = 0 \end{array} \right.,
\label{eq_perm_bU2}
\end{eqnarray}
where $\widetilde{I}_D({\bf x},{\bf h},{\bf p}_m,\omega)$ and $\widetilde{I}_U({\bf x},{\bf h},{\bf p}_m,\omega)$ 
are the encoded SODCIGs. They are defined as:
\begin{eqnarray}
\begin{array}{l}
\widetilde{I}_D({\bf x},{\bf h},{\bf p}_m,\omega) = \sum_{x_m,y_m}I_D({\bf x},{\bf h},x_m,y_m)\beta, \\
\widetilde{I}_U({\bf x},{\bf h},{\bf p}_m,\omega) = \sum_{x_m,y_m}I_U({\bf x},{\bf h},x_m,y_m)\beta,
\end{array}
\end{eqnarray}
where $\beta=e^{i\gamma({\bf x},x_m,y_m,{\bf p}_m,\omega)}$ is chosen to be the random phase-encoding function, with $\gamma({\bf x},x_m,y_m,{\bf p}_m,\omega)$ being a  uniformly distributed random sequence in ${\bf x}$, $x_m$, $y_m$ and $\omega$; the variable ${\bf p}_m$ is the index of different realizations of the random sequence. Recursively solving Equations \ref{eq_perm_bD2} and \ref{eq_perm_bU2} gives us the encoded areal source data ${\widetilde Q}_D({\bf x},{\bf p}_m,\omega)$ and encoded areal receiver data ${\widetilde Q}_U({\bf x},{\bf p}_m,\omega)$, which can be collected at any depth. 

In image-space wave-equation tomography, the image-space phase-encoded areal data sets are downward continued using the one-way wave equation. The background image is produced by cross-correlating the two wavefields and summing images for all realizations ${\bf p}_m$, as follows:
\begin{eqnarray}
\begin{array}{l}
I_{\rm me}({\bf x},{\bf h}) = \sum_{{\bf p}_m,\omega} \overline{\widetilde{D}({\bf x},{\bf p}_m,\omega)} \widetilde{U}({\bf x},{\bf p}_m,\omega). \label{eq_ic_dI_me}
\end{array}
\end{eqnarray}

The initial condition for modeling simultaneous events is set by regularly selecting SODCIGs in the prestack image. The amount of crosstalk in the image $I_{\rm me}({\bf x},{\bf h})$ can be controlled by choosing a convenient sampling interval for SODCIGs used simultaneously for the modeling.  For instance, if only one reflector is present and the correct velocity is used, no crosstalk is generated if the SODCIG interval is greater than twice the maximum subsurface offset of the prestack image. In the extreme case, when an incorrect velocity is used and the reflector's energy spreads through the whole range of subsurface offsets, crosstalk is not generated if the the SODCIG interval is greater than four times the maximum subsurface offset. In the presence of more than one reflector, crosstalk between reflectors occurs, regardless of the distance between SODCIGs input to modeling. By phase-encoding the reflectors, we can mitigate the crosstalk.
\par
To phase-encode the reflectors it is necessary to pick some significant reflectors in the prestack migrated data. This implies a horizon-based approach for the prestack exploding-reflector modeling. In velocity-model updating, the idea of using some key reflectors to extract the residual-moveout information is a established strategy \cite[]{stork01,kosloff:1511,jiao:VE183}. 
 
The perturbed image is obtained by applying the chain rule to Equation \ref{eq_ic_dI_me}. The slowness perturbation is computed by applying the adjoint of the tomographic operator, ${\bf T}^\ast$, to the image perturbation.
 
\section{Numerical examples}
We test the image-space wave-equation tomography using image-space encoded data on a smoothed version of the Marmousi model, computed by applying a 200 m 2D-median filter to the slowness model. One-way data were synthesized considering a reflectivity computed from the Marmousi stratigraphic velocity model. We modeled 376 shots, ranging from 0 to 9000 m, with 24 m spacing. We used split-spread acquisition geometry, with maximum offset of 6600 m and receiver spacing of 24 m. 
 
Figure \ref{fig:islow}(a) shows the true slowness model. The background velocity model is equal to the correct velocity model above 2400 m depth and above the anticline with apex at (x = 6000 m, z = 1850m). Therefore, the slowness perturbation is zero in this portion of the model. Below these horizons, the background model is characterized by a smoother version of the original Marmousi model, computed with a 400 m 2D-median filter and scaled down by a factor of $5\%$. Figure \ref{fig:islow}(b) shows the background slowness. By using this background slowness model, we assume that a layer striping approach has been used and that the model is accurately defined up to a certain horizon, as usually occurs in projects of velocity model building. The slowness perturbation, computed by taking the difference between the correct and background slownesses, is shown in Figure \ref{fig:dslow}(a). In the part where the slowness perturbation is different from zero, the ratio between the true and the background slowness ranges approximately from 0.8 to 0.92 (Figure \ref{fig:dslow}(b)). Notice that the minimum depth is 1500 m. Henceforth, all the figures will be displayed with a minimum depth of 1500 m.
\activeplot{islow}{width=\columnwidth,height=0.4\textheight}{ER}{a) Correct slowness; b) Background slowness.}
\activeplot{dslow}{width=\columnwidth,height=0.4\textheight}{ER}{a) Slowness perturbation; b) Ratio between the true and the background slownesses.}
 
To compute the image-space phase-encoded data, we picked $10$ reflectors in the non-zero slowness perturbation part, in the prestack image computed with the $376$ original shots using the background slowness model. Figure \ref{fig:bimg1} shows the background image (Figure \ref{fig:bimg1}a) computed with shot-profile migration. The panel on the left corresponds to the zero-subsurface-offset section, and the panel on the right is the SODCIG at CMP position 5500 m. Notice the effects of using an inaccurate background slowness. The reflector at (7000 m, 2500 m) is pulled up, as are the subjacent reflectors. In the SODCIG, the energy is not focused at the zero subsurface offset.
\par
Figure \ref{fig:bimg1}(b) shows the picked reflectors used to model the image-space phase-encoded data. This image is used as input for the rotation of the reflectors in the SODCIGs with respect to the apparent geological dip, and the results are used as the initial conditions to model the image-space phase-encoded data, as discussed by \cite{biondi2006,biondi2007}. Figure \ref{fig:rimg1} shows the initial conditions for the prestack modeling. Figure \ref{fig:rimg1}(a) shows the initial condition for modeling the receiver wavefield, and Figure \ref{fig:rimg1}(b) shows the initial condition for modeling the source wavefield.
\activeplot{bimg1}{width=\columnwidth,height=0.4\textheight}{CR}{a) Background image; b) Windowed image used to compute image-space-encoded data.}
\par
Two image-space phase-encoded data sets were synthesized using different parameters. One contains one random realization of phase-encoded areal shots initiated simultaneously with SODCIG sampling interval of 264 m and encoded according to the CMP position and reflector number, generating $11$ areal shot gathers. The other data set corresponds to two random realizations modeled with  SODCIG sampling interval of 840 m, composed of $70$ areal shot gathers. Because in the velocity inversion we consider the maximum subsurface offset to be 192 m, this data set is expected to generate less crosstalk. In some comparisons, we use just one random realization ($35$ areal shots) of the second data set. We use the two random realizations when comparing the results of the non-linear optimization of the slowness model. 
\par
\activeplot{rimg1}{width=\columnwidth,height=0.4\textheight}{CR}{a) Initial condition for modeling the receiver wavefield. b) Initial condition for modeling the source wavefield.}
In the slowness optimization problem, we compute the image perturbation by applying the DSO operator, $\bf O$, to the background image in the subsurface-offset domain. The corresponding objective functional is
\begin{eqnarray}
J = \frac{1}{2}||{\bf O \widehat{\bf I}}||^2 = \frac{1}{2}||h {\bf \widehat{\bf I}}||^2.
\end{eqnarray}
Since the DSO operator is independent of the slowness, the gradient of $J$ with respect to the slowness ${\bf s}$ is
\begin{eqnarray} 
{\bf \nabla} J = \left.\left(\frac{\partial {\bf I}}{\partial {\bf s}}\right|_{{\bf s}=\widehat{\bf s}}\right)^{\ast}{\bf O}^{\ast}{\bf O}\widehat{\bf I} = {\bf T}^{\ast}{\bf O}^{\ast}{\bf O}\widehat{\bf I}. \label{eq_grad_dso}
\end{eqnarray} 
Given that the computation of the DSO objective functional is fully automated, it can be minimized by using quasi-Newton methods. Here, we specifically use the constrained L-BFGS algorithm \cite[]{nocedal}. 
\par
To guarantee smoothness of the wave-equation tomography gradient, we use a B-spline representation with nodes located every 960 m in the $x$-direction and 16 m in the $z$-direction. 
\par
Figure \ref{fig:dimg1} shows the image perturbation computed by applying the forward tomographic operator, ${\bf T}$, and using the background slowness of Figure \ref{fig:islow}(b) and the known slowness perturbation of Figure \ref{fig:dslow}. Figure \ref{fig:dimg1}a shows the image perturbation computed in the shot-profile domain for the $376$ shots; Figure \ref{fig:dimg1}b shows the image perturbation computed in the image-space phase-encoded domain using 11 image-space phase-encoded gathers; and Figure \ref{fig:dimg1}c shows the image perturbation computed in the image-space phase-encoded domain using 35 image-space phase-encoded gathers. Notice that the dispersed crosstalk is stronger in Figure \ref{fig:dimg1}b than in Figure \ref{fig:dimg1}c.
\activeplot{dimg1}{width=\columnwidth,height=.6\textheight}{CR}{a) Image perturbation in the shot-profile domain; b) Image perturbation computed with $11$ image-space phase-encoded wavefields; c) Image perturbation computed with $35$ image-space phase-encoded wavefields.}
%Figure \ref{fig:twin-bimg-all}(a) is obtained by migrating the original $401$ shot gathers. Because of the inaccuracy of the slowness model, we can identify the mispositioning of the reflectors, especially beneath the Gaussian anomalies. Figure \ref{fig:twin-bimg-all}(b) is obtained by migrating the data-space plane-wave encoded gathers, where $61$ plane waves are migrated; the result is almost identical to that in Figure \ref{fig:twin-bimg-all}(a); Figure \ref{fig:twin-bimg-all}(c) is obtained by migrating the image-space encoded gathers. The image-space encoded areal source and receiver data are generated by simultaneously modeling $100$ randomly encoded unfocused SODCIGs, and $4$ realizations of the random sequence are used; hence we have $40$ image-space encoded areal gathers (each realization contains $10$ areal shots). The kinematics of the result look almost the same as those in Figure \ref{fig:twin-bimg-all}(a). However, notice the wavelet squeezing effect and the random noise in the background caused by the random phase encoding.
\par
Figure \ref{fig:dsadj} illustrates the normalized slowness perturbations obtained by applying the adjoint tomographic operator ${\bf T}^{\ast}$ to the image perturbations of Figure \ref{fig:dimg1}. Compare with the correct slowness perturbation of Figure \ref{fig:dslow}. Figure \ref{fig:dsadj}a is the predicted slowness perturbation found by back-projecting Figure \ref{fig:dimg1}a using all $376$ shot gathers; Figure \ref{fig:dsadj}b shows the back-projection of Figure \ref{fig:dimg1}b using $11$ image-space phase-encoded gathers; and Figure \ref{fig:dsadj}c shows the back-projection of Figure \ref{fig:dimg1}c using $35$ image-space phase-encoded gathers. Notice that we do not use a B-spline representation for the slowness perturbations. In general, the predicted slowness perturbation with image-space phase-encoded gathers shows a structure similar to that obtained with the original shot gathers. The differences can be credited, at first order, to the occurrence of residual crosstalk in the image-space phase-encoded perturbed image and to a sub-optimal number of selected reflectors for the prestack exploding-reflector modeling.
\activeplot{dsadj}{width=\columnwidth,height=.50\textheight}{CR}
{Normalized slowness perturbation obtained by applying the adjoint tomographic operator ${\bf T}^{\ast}$ on the image perturbations in Figure \ref{fig:dimg1}. a) Slowness perturbation computed from Figure \ref{fig:dimg1}a. b) Slowness perturbation computed from Figure \ref{fig:dimg1}b. c) Slowness perturbation computed from Figure \ref{fig:dimg1}c.}
\par 
Finally, we compare the optimized slowness models with the correct slowness model of Figure \ref{fig:islow}(a). After 5 non-linear iterations for both the $11$-gather image-space phase-encoded data set and $35$ areal shots (one random realization) of the $70$-gather image-space phase-encoded data set the optimization stopped because the difference between the objective functional of successive iterations was smaller than the machine precision. The number of function evaluations was 28 for the $11$-gather image-space phase-encoded data set, and 27 for the $35$-gather image-space phase-encoded data set. We also computed 2 non-linear iterations with a total of 6 function evaluations using the two random realizations of the $70$-gather image-space phase-encoded data set. To verify the accuracy of the resulting optimized slowness models, we also migrated the original shot gathers with the three optimized slownesses and also with the correct slowness.
\par
Figure \ref{fig:plot} displays the evolution of the objective functional with the non-linear iterations for the $11$-gather image-space phase-encoded data set (Figure \ref{fig:plot}(a)) and $35$-gather image-space phase-encoded data set (Figure \ref{fig:plot}(b)). For comparison, the value of the objective functional for the true velocity is also shown as dashed lines. The values are normalized according to the value of the objective functional of the first iteration for each case. The objective functional was reduced in 22$\%$ and 36$\%$ for the $11$-gather image-space phase-encoded data set and for the $35$-gather image-space phase-encoded data set, respectively. Notice that those values are 23$\%$ and 47$\%$, respectively, when using the true slowness model. The smaller difference between the final optimized value of the objective functional and the objective functional computed with the true slowness model for the optimization with the $11$-gather image-space phase-encoded data set, can be credited to the more severe crosstalk generated by this data set than the $35$-gather image-space phase-encoded data set. Even if the correct slowness model is used, residual crosstalk is amplified when applying the DSO operator. 
\activeplot{plot}{width=\columnwidth,height=.3\textheight}{NR}
{Evolution of the objective functional with the non-linear iterations. The dashed lines represent the value of the objective functional for the true slowness model. a) Normalized objective functional for the $11$-gather image-space phase-encoded data set. b) Normalized objective functional for the $35$-gather image-space phase-encoded data set. }
\par 
Figure \ref{fig:sfperm} shows the optimized slownesses and, for comparison purposes, the true slowness. Figure \ref{fig:sfperm}(a) displays the slowness model; Figure \ref{fig:sfperm}(b) is the slowness model obtained with the $11$-gather image-space phase encoded data; Figure \ref{fig:sfperm}(c) is the slowness model obtained with the $35$-gather image-space phase encoded data; and Figure \ref{fig:sfperm}(d) is the slowness model obtained with the $70$-gather image-space phase encoded data. In general, the predicted slownesses are reasonable. The predicted slowness using the $11$ image-space phase-encoded gathers shows slightly lower values than the other two predicted slownesses. Notice that the detailed slowness variation present in the true slowness is not recovered in the optimized slownesses, due to the B-spline representation of the gradient of the wave-equation tomography objective functional. In addition, as we are solving for the deeper portion of the model with dipping reflectors, it is likely that deficient illumination prevents us to obtain a more accurate slowness model. However, the slowness model obtained with the $70$ image-space phase encoded gathers recovers the low slowness values on the left of model better than the other two predicted slownesses.
\activeplot{sfperm}{width=\columnwidth,height=.75\textheight}{CR}
{True and optimized slownesses. a) True slowness model; b) Slowness model obtained with the $11$-gather image-space phase-encoded data. c) Slowness model obtained with the $35$gather -image-space phase-encoded data. d) Slowness model obtained with the $70$-gather image-space phase-encoded data.}
\par
Figure \ref{fig:hsfperm} shows the histograms of the ratio between the true and background interval slowness (continuous line) and between the true and predicted interval slownesses obtained with the $11$ (fine dash), $35$ (fine dot), and $70$ image-space phase-encoded gathers  (large dash) below the depth of 2400 m. 
\activeplot{hsfperm}{width=4in}{CR}
{Histograms of the slowness ratios between the true and background interval slowness (continuous line) and between the true and predicted interval slownesses obtained with the $11$ (fine dash), $35$ (fine dot), and $70$ image-space phase-encoded gathers (large dash).}
\par
The mean and standard deviation of the corresponding distributions are summarized in Table \ref{tab:medev}. In general, the predicted slownesses vary between $94\%$ to $100\%$ of the true slowness. The background slowness varies between $110\%$ to $116\%$ of the true slowness.
\begin{table}[h]
 \begin{center}
  \begin{tabular}{|c|c|c|}
  \hline
  Slowness ratio & mean & $\sigma$ \\ \hline
  background&0.884 & 0.025     \\ \hline
  $11$ gathers &1.030 & 0.033 \\ \hline
  $35$ gathers &1.027 & 0.032 \\ \hline
  $70$ gathers &1.013& 0.027 \\ \hline
  \end{tabular}
  \caption{Mean and standard deviation of the interval slowness ratio.}
  \label{tab:medev}
 \end{center}
\end{table} 

%\plot{dsfperm}{width=\columnwidth,height=.50\textheight}
%{Ratio between true slowness and the optimized slownesses of Figure \ref{fig:sfperm}. a) Slowness ratio obtained with the $11$-gather image-space phase-encoded data. b) Slowness ratio obtained with the $35$-gather image-space phase-encoded data. c) Slowness ratio obtained with the $70$-gather image-space phase-encoded data.}
\par
Figure \ref{fig:fimg1} displays the zero-subsurface-offset section after migration of the $376$ original shot gathers using the true slowness model (Figure \ref{fig:fimg1}(b)), the predicted slowness using $11$ image-space phase encoded gathers (Figure \ref{fig:fimg1}(c)), the predicted slowness using $35$ image-space phase encoded gathers (Figure \ref{fig:fimg1}(d)), and the predicted slowness using $70$ image-space phase encoded gathers (Figure \ref{fig:fimg1}(e)). For comparison, we also display in Figure \ref{fig:fimg1}(a) the zero-subsurface-offset section after migration with the background slowness of Figure \ref{fig:islow}(b). Notice that reflectors in the central portion of Figure \ref{fig:fimg1}(a) are pulled up when comparing to  Figure \ref{fig:fimg1}(b). The image obtained with the optimized slowness model computed with the $11$ image-space phase encoded gathers (Figure \ref{fig:fimg1}(c))presents pushed down reflectors around (4000 m, 2500 m) as a consequence of the lower values of the optimized slowness. In addition, in this image the undulating character of the reflector at (7000 m, 2600 m) reflects some velocity inaccuracy, when compared to Figures \ref{fig:fimg1}(b) and (e). 
\par
From top to bottom, Figure \ref{fig:fimg11} displays SODCIGs at 1500 m, 3500 m, 5500 m and 7500 m after migration of the $376$ original shot gathers, using the background slowness of Figure \ref{fig:islow}(b) in the first row, using the true slowness model in the second row, using the predicted slowness with $11$ image-space phase encoded gathers in the third row, using the predicted slowness with $35$ image-space phase encoded gathers in the fourth row, and using the predicted slowness with $70$ image-space phase encoded gathers in the fifth row. The subsurface-offset ranges from -192 m  to 192 m. The analysis of the SODCIGs in Figures \ref{fig:fimg11}(c) to (e) shows that better focusing is achieved when more image-space phase-encoded gathers are used in the wave-equation tomography.  
\activeplot{fimg1}{width=0.75\columnwidth,height=.750\textheight}{CR}
{Zero-subsurface-offset section after migration of the $376$ original shot gathers using: a) the true slowness model; b) the predicted slowness model of Figure \ref{fig:sfperm}(a); c)  the predicted slowness model of Figure \ref{fig:sfperm}(b); and d) the predicted slowness model of Figure \ref{fig:sfperm}(c).} 
\activeplot{fimg11}{width=0.75\columnwidth,height=.750\textheight}{CR}
{Subsurface-offset gathers after migration of the $376$ original shot gathers. From top to bottom: in the first row, using the background slowness model; the second row, using the true slowness model; in the third row, using the predicted slowness model of Figure \ref{fig:sfperm}(b); in the fourth row, using the predicted slowness model of Figure \ref{fig:sfperm}(c); and in the fifth row, using the predicted slowness model of Figure \ref{fig:sfperm}(d).}
\par
Figure \ref{fig:fang1} displays the angle-domain common-image gathers (ADCIGs) taken at the same CMP position as SODCIGs of Figure \ref{fig:fimg11}. From top to bottom, Figure \ref{fig:fang1} displays ADCIGs  after migration using the background slowness in the first row, using the true slowness model in the second row, using the predicted slowness with $11$ image-space phase encoded gathers in the third row, using the predicted slowness with $35$ image-space phase encoded gathers in the fourth row, and using the predicted slowness with $70$ image-space phase encoded gathers in the fifth row. Notice that migration with the predicted slowness using $70$ image-space phase encoded gathers shows virtually no residual moveout. For the case of predicted slowness using $11$ and $35$ image-space phase encoded gathers some residual moveout occurs for CMP position 5500 m.
\activeplot{fang1}{width=0.75\columnwidth,height=.750\textheight}{CR}
{ADCIGs after angle transformation of the SODCIGs of Figure \ref{fig:fimg11}. From top to bottom: in the first row, using the background slowness model; the second row, using the true slowness model; in the third row, using the predicted slowness model of Figure \ref{fig:sfperm}(b); in the fourth row, using the predicted slowness model of Figure \ref{fig:sfperm}(c); and in the fifth row, using the predicted slowness model of Figure \ref{fig:sfperm}(d).}

The accuracy of the optimized slowness improves when using more phase-encoded gathers in the wave-equation tomography, or, in other words, when the crosstalk in the perturbed image is less severe, as Figures \ref{fig:dimg1} and \ref{fig:dsadj} suggest. Figure \ref{fig:dso} shows the perturbed image computed by applying the DSO operator on the image migrated with the background slowness of Figure \ref{fig:islow}(b). The panel on the left corresponds to the subsurface-offset -144 m and the panel on the right is the SODCIG taken at 5500 m. Figure \ref{fig:dso}(a) shows the perturbed image using $11$ image-space phase-encoded gathers; Figure \ref{fig:dso}(b) shows the perturbed image using $35$ image-space phase-encoded gathers; and Figure \ref{fig:dso}(c) shows the perturbed image using $70$ image-space phase-encoded gathers.  Notice how the signal-to-noise ratio improves as more phase-encoded gathers are used.  The SODCIG of the perturbed image of Figure \ref{fig:dso}(a) presents coherent events, related to unattenuated crosstalk, curving upward at z = 2700 m; these events are not present in Figures \ref{fig:dso}(b) and (c). If these events are sufficiently incoherent along the subsurface-offset sections, a two-dimensional filter could be applied to attenuate them. In that case, a new objective function should be defined. This deserves future investigation.

\activeplot{dso}{width=\columnwidth,height=.6\textheight}{CR}
{Perturbed images computed with the DSO operator. a) Perturbed image using $11$ image-space phase-encoded gathers. b) Perturbed image using $35$ image-space phase-encoded gathers. c) Perturbed image using $70$ image-space phase-encoded gathers.}

%\plot{twin-dimg-offdso-all}{width=\columnwidth,height=0.8\textheight} {The DSO image perturbations. The coherent energy at non-zero offsets indicates velocity errors. Panel (a) is obtained using the original shot gathers; Panels (b) and (c) are obtained using the data-space encoded gathers and the image-space encoded gathers, respectively. [{\bf CR}]}

%Figure \ref{fig:twin-dslw-offdso-all} shows the negative gradients of the DSO objective functional ($-\nabla J_{\rm DSO}$) obtained by back-projecting the DSO image perturbations shown in Figure \ref{fig:twin-dimg-offdso-all}. For comparison, Figure \ref{fig:twin-dslw-offdso-all}(a) shows the exact slowness perturbation, which is the same as Figure \ref{fig:twin-dslw-all}(a); Figure \ref{fig:twin-dslw-offdso-all}(b) shows the result obtained in the original shot-profile domain; Figure \ref{fig:twin-dslw-offdso-all}(c) shows the result obtained in the data-space phase-encoding domain, which is almost identical to Figure \ref{fig:twin-dslw-offdso-all}(b); Figure \ref{fig:twin-dslw-offdso-all}(d) shows the result obtained in the image-space phase-encoding domain. The result is also similar to Figure \ref{fig:twin-dslw-offdso-all}(b), though the unattenuated crosstalk and the random noise make the gradient less well behaved than those in Figure \ref{fig:twin-dslw-offdso-all}(b) and (c). Most important, the gradient in Figure \ref{fig:twin-dslw-offdso-all}(d) is pointing towards the correct direction, which is crucial for a gradient-based optimization algorithm to converge to the correct solution.

%\plot{twin-dslw-offdso-all}{width=\columnwidth,height=0.8\textheight}
%{The negative DSO gradients obtained using different methods. Panel (a) shows the exact slowness perurbation; Panel (b) shows the result obtained using the original shot gathers; Panels (c) and (d) show the results obtained using the data-space phase encoded gathers and the image-space phase encoded gathers, respectively. [{\bf CR}]} 

\section{Conclusions}
We present a cost-effective method to perform image-space wave-equation tomography using image-space phase-encoded shot gathers. One important advantage is that we are able to synthesize a much smaller data set while still keeping necessary velocity information for migration velocity analysis; hence the computational cost of performing image-space wave-equation tomography can be significantly reduced. Our results show that using the image-space phase-encoded wavefields provides a gradient of the tomography objective functional similar to that computed using the original shot gathers, but at significantly lower cost. We also show that the accuracy of optimized slowness depends on the amount of residual crosstalk in the perturbed image.

\bibliographystyle{seg}
%\bibliography{SEP,MISC,GEOTLE,SEG,GEOPHYSICS,sep136a}
\bibliography{sep136a}
%\end{document}
