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

\begin{document}
\title{Wave-equation tomography using image-space phase encoded data}

\renewcommand{\thefootnote}{\fnsymbol{footnote}} 
\author{Claudio Guerra *, Yaxun Tang, and Biondo Biondi, Stanford University}

\footer{Generalized sources and wave-equation tomography}

\lefthead{Guerra et al.}
\righthead{Image-space wave-equation tomography}


\maketitle

\begin{abstract}
We present the image-space wave-equation tomography to the generalized source domain, where a smaller number of synthesized shot gathers are generated. Specifically, we generate synthesized shot gathers by image-space phase encoding. The comparison of the gradients of the tomography objective functional obtained using image-space encoded gathers with that obtained using the original shot gathers shows that those encoded shot gathers can be used in wave-equation tomography problems. The advantages of using image-space encoded data is the decreased computational when compared to that of the original shot gathers and the ability of performing velocity analysis in a target-oriented way. We illustrate our examples on Marmousi model.
\end{abstract}
\section{introduction}
Wave-equation tomography has the potential to accurately estimate the velocity model in complex geological scenarios where ray-based traveltime tomography is prone to fail. That 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 from multi-pathing issues. However, despite its theoretical advantages, wave-equation tomography is still computationally challenging.

Wave-equation tomography can be performed on 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. Independently on the domain of application of wave-equation tomography, the use of phase encoded data can substantially decrease the computational cost \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 either by data-space phase encoding \cite[]{whitmore1995,romero2000} or image-space phase encoding \cite[]{guerra2008a}. Image-space phase encoding is used to combine data synthesized by the prestack exploding-reflector model \cite[]{biondi2006,biondi2007}. The prestack exploding-reflector model synthesizes data and the corresponding source function by upward propagating wavefields, having  as initial condition a prestack image computed with wave-equation migration.

In this paper, we show that image-space encoded data can provide the gradient of the tomographic objective functional similar to that obtained in the original shot-profile domain, but with less computational cost. We briefly discuss the theory of wave-equation tomography in the image-space domain, then we explain the prestack exploding-reflector modeling and how the image-space phase encoding can be used to speed up wave-equation tomography in the image domain. We use Marmousi model to illustrate the theory.

\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}),
\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 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 ${\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 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}
{\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 equation \cite[]{claerbout1971} and the image is computed by applying the cross-correlation imaging condition:
\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 $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; ${\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( \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 perturbed one-way wave equation \cite[]{shen2004}.%, which are 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}
%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.

To evaluate the adjoint tomographic operator ${\bf T}^\ast$, 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 perturbed one-way wave equation to get the slowness perturbation $\Delta s({\bf x})$ \cite[]{shen2004}. More detailed information on how to evaluate the forward and adjoint operators can be found in \cite {yax2008}.

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




%\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 the prestack exploding-reflector modeling \cite[]{biondi2006} is to model the data and the 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 the 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}$; ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector.

Since the size of the migrated image volume can be very big in practice and there are usually many reflectors in the subsurface, modeling each reflector and each SODCIG may generate a data set even bigger 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 could be used to attenuate the crosstalk. 
%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(x,y,z=0,{\bf p}_m,\omega)$ and encoded areal receiver data ${\widetilde Q}_U(x,y,z=0,{\bf p}_m,\omega)$, which can be collected at any depth level. 

In the image-space wave-equation tomography, the image-space encoded areal data sets are downward continued using the one-way wave equation. The 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} \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. Due to the localized nature of the initial conditions used to model prestack exploding reflector data, the amount of crosstalk in the image $I_{\rm me}({\bf x},{\bf h})$ can be controlled by conveniently choosing the sampling interval of SODCIGs used simultaneously for the modeling.  For the case of one reflector, if the SODCIG interval is greater than twice the maximum subsurface offset of the prestack image, no crosstalk is generated. 

The perturbed image is obtained by applying chain rule on 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 200m 2D-median filter on 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 9000m, at every 24m. We used split-spread acquisition geometry, with maximum offset of 6600m and receiver spacing of 24m. 

Figure \ref{fig:islow}(a) shows the slowness model. The background velocity model is equal to the correct velocity model above 2400m depth and the anticlinal with apex at $(x=6000m,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 400m 2D-median filter and scaled down by a factor of $5\%$. Figure \ref{fig:islow}(b) shows the slowness perturbation computed by taking the difference between the correct and the background slowness.
\plot{islow}{width=\columnwidth,height=0.165\textheight}{a) Correct slowness; b) Slowness perturbation.}
To compute the image-space encoded data, we picked $8$ reflectors in the non-zero slowness perturbation part, on 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 and the windowed data (Figure \ref{fig:bimg1}b) used as initial condition for the image-space encoded data. Two image-space encoded data were synthesized with different parameters. For both, just one random realization per areal shot was generated. One contains areal shots initiated simultaneously at every $11$ CMP positions and encoded according to the CMP position and reflector number, originating $11$ areal shot gathers. The other dataset was modeled having SODCIG sampling interval of $31$, or 744m. Given that the SODCIG sampling interval is greater than twice the maximum subsurface offset, 720m, the second data is expected to generate less crosstalk.
\plot{bimg1}{width=\columnwidth,height=0.165\textheight}{a) Background image; b) Windowed image used to compute image-space encoded data.}
Figure \ref{fig:dimg1} shows the image perturbation computed by applying the forward tomographic operator, ${\bf T}$, and using the background slowness and the known slowness perturbation of Figure \ref{fig:islow}b. The panel on the left corresponds to the zero-subsurface offset section and the panel on the right is the SODCIG at CMP position 6500m. 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 encoded domain using 11 image-space encoded gathers; and Figure \ref{fig:dimg1}c shows the image perturbation computed in the image-space encoded domain using 31 image-space encoded gathers. Notice that the dispersed crosstalk is stronger in Figure \ref{fig:dimg1}b than in Figure \ref{fig:dimg1}c.
\plot{dimg1}{width=\columnwidth,height=0.33\textheight}{a) Image perturbation in the shot-profile domain; b) Image perturbation computed with $11$ image-space encoded wavefields; c) Image perturbation computed with $31$ image-space 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.

Figure \ref{fig:dsadj} illustrates the predicted slowness perturbations 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:islow}(b). Figure \ref{fig:dsadj}a is the predicted slowness perturbation by back-projecting Figure \ref{fig:dimg1}a using all $376$ shot gathers; Figure \ref{fig:dsadj}b refers to the by back-projection of Figure \ref{fig:dimg1}b using $11$ image-space encoded gathers; and Figure \ref{fig:dsadj}c shows the back-projection of Figure \ref{fig:dimg1}c using $31$ image-space encoded gathers. The predicted slowness perturbation with image-space encoded gathers, in general, shows a similar structure to that obtained with the original shot gathers. The differences can be credited, at first order, to a sub-optimal number of selected reflectors for the prestack exploding reflector modeling. It is possible that, by selecting more reflectors, we can decrease this difference.  In addition, in spite of the stronger crosstalk residual in the image perturbation for the $11$ image-space encoded gathers, the predicted slowness perturbation is very similar for the image-space encoded data
\plot{dsadj}{width=\columnwidth,height=0.33\textheight}
{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.}
Finally, in Figure \ref{fig:dsdso}, we show the gradients of the DSO objective function computed by back-projecting the image perturbation computed with DSO. Figure \ref{fig:dsdso}a is the DSO gradient obtained in the original shot-profile domain, using $376$ shots; Figure \ref{fig:dsdso}b shows the DSO gradient computed with $11$ image-space encoded wavefields, and Figure \ref{fig:dsdso}c shows the DSO gradient computed with $31$ image-space encoded wavefields. Even though the DSO gradient computed with image-space phase encoded wavefields is less well behaved than that of the shot-profile domain, it is important to notice that, the gradient is pointing towards the correct direction, which is crucial for a gradient-based optimization algorithm to converge to the correct solution.
\plot{dsdso}{width=\columnwidth,height=0.33\textheight}
{Gradients of the DSO objective function. a) using $376$ original gathers; b) using $11$ image-space encoded wavefields; c) using $31$ image-space encoded wavefields.}

%\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 presented 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, provide a gradient of the tomography objective functional similar to that computed using the original shot gathers, but at significantly lower cost. The correct gradient is thus important for the  gradient-based optimization algorithm to converge to the correct solution.

\section{acknowledgments}
The authors would like to thank the sponsors of the Stanford Exploration Project for their financial support.
\bibliographystyle{seg}
%\bibliography{SEP,MISC,GEOTLE,SEG,GEOPHYSICS,sep136a}
\bibliography{sep136a}
\end{document}
