%% This file is automatically generated. Do not edit!
\documentclass[12pt]{geophysics}

\usepackage{amsmath}

\setfigdir{Fig}

\begin{document}
\lefthead{Tang et al.}
\righthead{Image-space wave-equation tomography}

\title{Image-space wave-equation tomography in the generalized source domain}
\author{Yaxun Tang, Claudio Guerra, and Biondo Biondi}
\footer{SEP-136}

\maketitle

\begin{abstract}
We extend the theory of image-space wave-equation tomography to the generalized source domain, where a smaller number of synthesized shot gathers are generated either by data-space phase encoding or image-space phase encoding. We compare the gradients of the tomography objective functional obtained using both data-space and image-space encoded gathers with that obtained using the original shot gathers. We show that with those encoded shot gathers we can obtain a gradient similar to that 
computed in the original shot-profile domain, but at lower computational cost. 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. Wave-equation-based 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,mora1989,woodward1992,pratt1999} or in the image-space domain \cite[]{biondi1999,shen:VE49}. The image-space approach minimizes the residual in the image domain obtained after migration. Some measurement of the quality of the image is used to evaluate the accuracy of the velocity model.

The goal of this paper is to extend the theory of image-space wave-equation tomography from the conventional shot-profile domain \cite[]{shen2004,shen2005b} to the generalized source domain, where a smaller number of synthesized shot gathers make the tomographic velocity update substantially faster. The generalized source domain can be obtained either by data-space phase encoding \cite[]{whitmore1995,romero2000,zhang2005,duquet2006,liu2006} or image-space phase encoding  \cite[]{biondi2006,biondi2007,guerra2008}. We show that in these generalized source domains, we can obtain gradients, which are used for updating the velocity model, similar to that obtained in the original shot-profile domain, but with less computational cost. We use data modeled with one-way wave-equation on Marmousi model to illustrate our results.

\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, ${\bf \Delta I}$, defined in the image space. The residual field is derived from the background image, ${\bf I}$, which is computed with a background slowness (or the current estimate of the slowness). The general form of the residual field is \cite[]{biondi2008}
\begin{eqnarray}
{\bf \Delta 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 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}.
\par
Under $\ell_2$ norm, the tomography objective function 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}
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)^{\dag} \left( {\bf I}-{\bf F}({\bf I})\right), \label{eq_obj_grad_gen}
\end{eqnarray}
where $^{\dag}$ denotes the adjoint. 
As for the DSO method, the linear operator ${\bf O}$ is independent of the slowness, the gradient reads
\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)^{\dag}{\bf O}^{\dag}{\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}$.

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}$, 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}
{\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}, in the shot profile domain, to generalize the expression of the tomographic operator to the generalized source domain.

%--------------------------------------------------------------------
\section{tomography with the encoded wavefields}
%--------------------------------------------------------------------
To reduce the cost of wave-equation tomography in the image-space, we extend the theory to the
generalized source domain, where a smaller number of synthesized shot gathers are used for computing the gradient. 
We discuss the data-space phase-encoding method and the image-space phase-encoding method.

\subsection{Data-space encoded wavefields}
The data-space encoded shot gathers are obtained by linear combination of the original shot gathers after phase encoding \cite[]{whitmore1995,romero2000,zhang2005,duquet2006,liu2006}. 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 equations:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right){\widetilde D}({\bf x},{\bf p}_s,\omega) = 0 \\
{\widetilde D}(x,y,z=0,{\bf p}_s,\omega) = 
\sum_{{\bf x}_s}\overline{f_s(\omega)\delta({\bf x}-{\bf x}_s) \alpha({\bf x}_s,{\bf p}_s,\omega)} \end{array} \right.,
\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){\widetilde U}({\bf x},{\bf p}_s,\omega) = 0 \\
{\widetilde U}(x,y,z=0,{\bf x}_s,\omega) = 
\sum_{{\bf x}_s}Q(x,y,z=0,{\bf x}_s,\omega) \alpha({\bf x}_s,{\bf p}_s,\omega) \end{array} \right.,
\end{eqnarray}
where ${\widetilde D}({\bf x},{\bf p}_s,\omega)$ and ${\widetilde U}({\bf x},{\bf p}_s,\omega)$ are the encoded source and receiver wavefields respectively, 
 $\alpha({\bf x}_s,{\bf p}_s,\omega)$ is the phase-encoding function, and ${\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. 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}\sum_{\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 $c=\omega$ for plane-wave phase encoding and $c=1$ for random phase encoding \cite[]{yaxun2008b}. It has been shown that a smaller number of encoded data produces similar images as the original shot gathers \cite[]{romero2000,etgen2005,liu2006}. Therefore, the following relation approximately holds:
\begin{eqnarray}
I({\bf x},{\bf h}) \approx I_{\rm de}({\bf x},{\bf h}).
\end{eqnarray}
From Equation \ref{eq_ic_dI_de}, the perturbed image can be easily obtained as follows:
\begin{eqnarray}
\Delta I_{\rm de}({\bf x},{\bf h})  &=& 
\sum_{{\bf p}_s}\sum_{\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 equations defined by Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU}.

\subsection{Image-space encoded wavefields}
The image-space encoded gathers are obtained using the prestack exploding-reflector modeling method introduced by \cite{biondi2006} and \cite{biondi2007}. The general idea of this method is to model the data and the corresponding source function that are related to only one event in the subsurface, where a single unfocused SODCIG (obtained with an inaccurate velocity model) 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\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\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\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\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. 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)$.

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 one by one 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,guerra2008b}, random phase encoding could be used to attenuate the crosstalk. 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\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\right) Q_D({\bf x},{\bf p}_m,\omega) 
= \widetilde{I}_D({\bf x},{\bf h},{\bf p}_m,\omega) \\
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\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\right) Q_U({\bf x},{\bf p}_m,\omega) 
= \widetilde{I}_U({\bf x},{\bf h},{\bf p}_m,\omega) \\
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 follows:
\begin{eqnarray}
\widetilde{I}_D({\bf x},{\bf h},{\bf p}_m,\omega) &=& \sum_{x_m}\sum_{y_m}I_D({\bf x},{\bf h},x_m,y_m)\beta({\bf x},x_m,y_m,{\bf p}_m,\omega), \\
\widetilde{I}_U({\bf x},{\bf h},{\bf p}_m,\omega) &=& \sum_{x_m}\sum_{y_m}I_U({\bf x},{\bf h},x_m,y_m)\beta({\bf x},x_m,y_m,{\bf p}_m,\omega),
\end{eqnarray}
where $\beta({\bf x},x_m,y_m,{\bf p}_m,\omega)=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 $Q_D(x,y,z=0,{\bf p}_m,\omega)$ and areal receiver data $Q_U(x,y,z=0,{\bf p}_m,\omega)$, which can be collected on the surface.

The synthesized new data sets are downward continued using the same one-way wave equation defined by Equations \ref{eq_wveq_bD} and \ref{eq_wveq_bU} (with different boundary conditions) as follows:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right){\widetilde D}({\bf x},{\bf p}_m,\omega) = 0 \\
\widetilde {D}(x,y,z=0,{\bf p}_m,\omega) = 
\overline{Q_D(x,y,z=0,{\bf p}_m,\omega)} \end{array} \right.,
\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){\widetilde U}({\bf x},{\bf p}_m,\omega) = 0 \\
{\widetilde U}(x,y,z=0,{\bf x}_s,\omega) = 
Q_U(x,y,z=0,{\bf p}_m,\omega) \end{array} \right.,
\end{eqnarray}
where $\widetilde{D}({\bf x},{\bf p}_m,\omega)$ and $\widetilde{U}({\bf x},{\bf p}_m,\omega)$ are the downward continued areal source and areal receiver wavefields for realization ${\bf p}_m$. The image is produced by cross-correlating the two wavefields and summing images for all realization ${\bf p}_m$ as follows:
\begin{eqnarray}
I_{\rm me}({\bf x},{\bf h}) = \sum_{{\bf p}_m}\sum_{\omega} \widetilde{D}({\bf x},{\bf p}_m,\omega) \widetilde{U}({\bf x},{\bf p}_m,\omega). \label{eq_ic_dI_me}
\end{eqnarray}
The crosstalk artifacts can be further attenuated if the number of ${\bf p}_m$ is large; therefore, approximately,
the image obtained by migrating the image-space encoded gathers is kinematically equivalent to the image obtained in the shot-profile domain.

From Equation \ref{eq_ic_dI_me}, the perturbed image is easily obtained as follows:
\begin{eqnarray}
\Delta I_{\rm me}({\bf x},{\bf h})  &=& 
\sum_{{\bf p}_m}\sum_{\omega} 
\left( \Delta \widetilde{D}({\bf x}-{\bf h},{\bf p}_m,\omega) \widehat{\widetilde{U}}({\bf x}+{\bf h},{\bf p}_m,\omega) + \right. \nonumber \\
& & \left. \widehat{\widetilde{D}}({\bf x}-{\bf h},{\bf p}_m,\omega) \Delta \widetilde{U}({\bf x}+{\bf h},{\bf p}_m,\omega) \right),
\end{eqnarray}
where $\widehat{\widetilde{D}}({\bf x},{\bf p}_m,\omega)$ and $\widehat{\widetilde{U}}({\bf x},{\bf p}_m,\omega)$
are the image-space encoded background source and receiver wavefields;
$\Delta \widetilde{D}({\bf x},{\bf p}_m,\omega)$ and $\Delta \widetilde{U}({\bf x},{\bf p}_m,\omega)$ are the perturbed source and receiver wavefields
in the image-space phase-encoding domain, which satisfy the perturbed one-way wave equations defined by Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU}.

\section{Numerical examples}
We test the image-space wave-equation tomography in the generalized source domain on a simple model which contains only one reflector located at $z=1500$ m.
Figure \ref{fig:twin-slow} shows the correct slowness model.
The slowness model consists of a constant background slowness $1/2000$ s/m and two Gaussian anomalies 
located at $(x=-800,z=800)$ and $(x=800,z=800)$ respectively. The left anomaly has $5\%$ higher slowness, while the right one has $5\%$ lower slowness.
We modeled $401$ shots ranging from $-4000$ m to $4000$ m, with a shot interval $20$ m. 
The receiver locations also range from $-4000$ m to $4000$ m, but with a $10$ m interval. The receivers are fixed for all shots
to mimic a land acquisition geometry. 

\plot{twin-slow}{width=\textwidth,height=0.225\textheight}{The correct slowness model. The slowness model consists of a constant background slowness ($1/2000$ s/m)
and two $5\%$ Gaussian anomalies. [{\bf ER}]}

Figure \ref{fig:twin-bimg-all} shows the migrated images in different domains computed with a background slowness $\widehat s = 1/2000$ s/m. 
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.

\plot{twin-bimg-all}{width=\textwidth,height=0.8\textheight}
{Migrated image cubes with a constant background slowness ($\widehat s = 1/2000$ s/m).
Panel (a) is the result obtained in the original shot-profile domain;
Panel (b) is the result obtained by migrating $61$ plane waves, while panel (c) is obtained by migrating $40$ image-space encoded areal gathers. [{\bf CR}]}

Figure \ref{fig:twin-dimg-all} shows the image perturbations obtained by applying the forward tomographic operator ${\bf T}$ in different domains.
For this example, we assume that we know the correct slowness perturbation $\Delta {\bf s}$, which is obtained by subtracting the background slowness $\widehat {\bf s}$
from the correct slowness ${\bf s}$. Figure \ref{fig:twin-dimg-all}(a) shows the image perturbation computed with the original $401$ shot gathers; notice the relative
$90$ degree phase rotation compared to the background image shown in Figure \ref{fig:twin-bimg-all}(a). Figure \ref{fig:twin-dimg-all}(b) is the result 
obtained by using $61$ data-space plane-wave encoded gathers; the result is almost identical to Figure \ref{fig:twin-dimg-all}(a). 
Figure \ref{fig:twin-dimg-all}(c) shows the result computed with $40$ image-space encoded gathers; the kinematics are also similar to those in
Figure \ref{fig:twin-dimg-all}(a).

\plot{twin-dimg-all}{width=\textwidth,height=0.8\textheight}
{The image perturbations obtained by applying the forward tomographic operator ${\bf T}$ to the correct slowness perturbations in different domains. 
Panel (a) shows the image perturbation obtained using the original shot gathers,
while panels (b) and (c) are obtained using the data-space encoded gathers and image-space encoded gathers, respectively. [{\bf CR}]}

Figure \ref{fig:twin-dslw-all} illustrates the predicted slowness perturbations by 
applying the adjoint tomographic operator ${\bf T}^{\prime}$ to the image perturbations obtained in Figure \ref{fig:twin-dimg-all}.
For comparison, Figure \ref{fig:twin-dslw-all}(a) shows the correct slowness perturbation, i.e., $\Delta {\bf s} = {\bf s} - \widehat{\bf s}$;
Figure \ref{fig:twin-dslw-all}(b) is the predicted slowness perturbation by back-projecting Figure \ref{fig:twin-dimg-all}(a) using all $401$ shot gathers;
Figure \ref{fig:twin-dslw-all}(c) is the result by back-projecting Figure \ref{fig:twin-dimg-all}(b) 
using all $61$ data-space plane-wave encoded gathers and is almost identical to 
Figure \ref{fig:twin-dslw-all}(b); Figure \ref{fig:twin-dslw-all}(d) shows the result by back-projecting Figure \ref{fig:twin-dimg-all}(c) 
using all $40$ image-space encoded gathers. 
The result is also similar to Figure \ref{fig:twin-dslw-all}(b). However, notice that Figure \ref{fig:twin-dslw-all}(d) shows a slightly less focused result than 
Figure \ref{fig:twin-dslw-all}(b) and (c), which might be caused by the unattenuated crosstalk and the pseudo-random noise 
presented in Figure \ref{fig:twin-dimg-all}(c).

\plot{twin-dslw-all}{width=\textwidth,height=0.8\textheight}
{The slowness perturbation obtained by applying the adjoint tomographic operator ${\bf T}^{\prime}$ on the image perturbations in Figure \ref{fig:twin-dimg-all}.
Panel (a) shows the exact slowness perturbation;
Panel (b) shows the slowness perturbation estimated by back-projecting the image perturbation shown in Figure \ref{fig:twin-dimg-all}(a);
Panel (c) shows the result obtained using the data-space plane-wave encoded gathers by back-projecting Figure \ref{fig:twin-dimg-all}(b) and 
Panel (d) shows the result obtained using the image-space encoded gathers by back-projecting Figure \ref{fig:twin-dimg-all}(c). [{\bf CR}]}

The final example we show is the comparison among the gradients of the objective functional obtained in different domains.
For simplicity, here we compare only the negative DSO gradients ($-\nabla J_{\rm DSO}$) defined by Equation \ref{eq_grad_dso}
(we compare  $-\nabla J_{\rm DSO}$ instead of $\nabla J_{\rm DSO}$, 
because $-\nabla J_{\rm DSO}$ determines the search direction in a gradient-based nonlinear optimization algorithm). 
Figure \ref{fig:twin-dimg-offdso-all} shows the 
DSO image perturbations computed as follows:
\begin{eqnarray}
\Delta {I}({\bf x},{\bf h}) = |{\bf h}|^2 \widehat {I}({\bf x},{\bf h}),
\end{eqnarray}
or in matrix form:
\begin{eqnarray}
\Delta {\bf I} = {\bf O}^{\prime}{\bf O}\widehat{\bf I},
\end{eqnarray}
where ${\bf O}$ is the DSO operator. Figure \ref{fig:twin-dimg-offdso-all}(a) is the result obtained in the original shot-profile domain, 
whereas Figure \ref{fig:twin-dimg-offdso-all}(b) and (c) are obtained in the data-space phase-encoding domain and the 
image-space phase-encoding domain, respectively. The coherent energy at non-zero offests are indicators of velocity errors. 

\plot{twin-dimg-offdso-all}{width=\textwidth,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=\textwidth,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 extend the theory of image-space wave-equation tomography to the generalized source domain. 
One important advantage of this new domain 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. 
We demonstrate how these new data sets can be generated by using both the data-space
phase encoding method and the image-space phase encoding method. 
Our preliminary tests on a simple synthetic model show that
with the synthesized gathers, we are able to obtain 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.
\bibliographystyle{seg}
%\bibliography{SEP,MISC,GEOTLE,SEG,GEOPHYSICS,sep136a}
\bibliography{sep136a}
\end{document}
