\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 demonstrate how to evaluate the wave-equation forward tomographic operator 
and its adjoint in this new domain.
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. 
The saving in cost is  
important for putting this theory into practical applications. We illustrate our examples on 
a simple model with Gaussian anomalies in the subsurface.
\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. 
Generally speaking, wave-equation tomography can be classified into two different categories based on the domain where it minimizes the residual. 
The domain can be either the data space or the image space. The data-space approach directly compares the modeled waveform with the recorded
waveform, and is widely known as waveform inversion, or data-space wave-equation tomography \cite[]{tarantola1987,mora1989,woodward1992,pratt1999}. 
The main disadvantage of the data-space approach is that in complex areas, the recorded waveforms can be 
very complicated and are usually of low signal-to-noise ratio (S/N), 
so matching the full waveform might be extremely difficult. On the other hand, the image-space approach, also known as image-space
wave-equation tomography, minimizes the residual 
in the image domain obtained after migration. The migrated image is often much simpler than the original data, because even with
a relatively inaccurate velocity, migration is able to (partially) collapse diffractions and enhance the S/N; thus the image-space wave-equation
tomography has the potential to mitigate some of the difficulties that we encounter in the data-space approach. 
Another advantage of the image-space approach is that the more efficient one-way wave-equation extrapolator can be used. 
In waveform inversion, however, the one-way propagator is difficult (if not impossible) to use because of its inability to model the multiple arrivals, although
some tweaks can be employed so that the one-way propagator can be applied to turning-wave tomography \cite[]{shragge2007}.

However, despite its theoretical advantages, 
image-space wave-equation tomography is still computationally challenging.
Each iteration of tomographic velocity updating is computationally
expensive and often converges slowly. Practical applications are still rare and small in scale \cite[]{biondi1999,shen2005b,albertin2006}.
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
or image-space phase encoding. 
For the data-space phase encoding, the synthesized shot gathers are obtained by linear combination of
the original shot gathers after some kind of phase encoding; in particular, here we mainly consider plane-wave phase encoding 
\cite[]{whitmore1995,zhang2005,duquet2006,liu2006} and random phase encoding \cite[]{romero2000}.
As the encoding process is done in the data space, we call it data-space phase encoding.
For the image-space phase encoding, the synthesized gathers are obtained by prestack exploding-reflector modeling \cite[]{biondi2006,biondi2007,guerra2008a},
where several subsurface-offset-domain common-image gathers (SODCIGs) and several reflectors are simultaneously demigrated to 
generate areal source and areal receiver gathers. 
To attenuate the cross-talk, the SODCIGs and the reflectors
have to be encoded, e.g., by random phase encoding. Because the encoding process is done in the image space, we call it image-space phase encoding.
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.

This paper is organized as follows: We first briefly review the theory of image-space wave-equation tomography. 
Then we discuss how to evaluate the forward
tomographic operator and its adjoint in the original
shot-profile domain. The latter is an important component in computing the gradient of the tomography objective functional. 
We then extend the theory
to the generalized source domain.
Finally, we show examples on a simple synthetic model.


\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 residual field measures the correctness of the background slowness;  
its minimum (under some norm, e.g. $\ell_2$) is achieved when a correct background slowness has been used for migration. 
There are many choices of the residual field, such as residual moveout in the Angle-Domain Common-Image Gathers (ADCIGs), 
differential semblance in the ADCIGs, 
reflection-angle stacking power (in which case we have to maximize the residual field, or minimize the negative stacking power), etc..
Here we follow a definition similar to that in \cite{biondi2008}, and define a general form of the residual field as follows:
\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. For example,
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}.
The subsurface-offset-domain DSO focuses the energy at zero offset, whereas the angle-domain DSO flattens the ADCIGs.

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

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

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

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

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

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

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

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

%--------------------------------------------------------------------
\section{tomography with the encoded wavefields}
%--------------------------------------------------------------------
It is clear from previous sections that the cost for computing the gradient of the objective function $J$ in the original shot-profile domain 
is at least twice the cost of a shot-profile migration, because to compute the perturbed wavefields, the background
wavefields are required. Because minimizing the objective function $J$ requires a considerable number of gradient and function evaluations,
image-space wave-equation tomography in the conventional shot-profile domain seems to be infeasible for large-scale 3-D applications, 
even with modern computer resources.
To reduce the cost and make this powerful method more practical, we extend the theory of image-space wave-equation tomography to the
generalized source domain, where a smaller number of synthesized shot gathers are used for computing the gradient. 
We discuss two different strategies to generate 
the generalized shot gathers, i.e., the data-space phase-encoding method and the image-space phase-encoding method, 
both of which can achieve considerable data reduction while still keeping the necessary kinematic information for
velocity analysis.

\subsection{Data-space encoded wavefields}
The data-space encoded shot gathers are obtained by linear combination of the
original shot gathers after phase encoding.  
For simplicity, we mainly consider
plane-wave phase-encoding \cite[]{whitmore1995,zhang2005,duquet2006,liu2006} and random phase-encoding \cite[]{romero2000}. 
Because of the linearity of the one-way wave equation with respect to the wavefield, 
the encoded source and receiver wavefields also satisfy the same one-way wave equations defined by Equations \ref{eq_wveq_bD} and \ref{eq_wveq_bU}, 
but with different boundary conditions:
\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, 
and $\alpha({\bf x}_s,{\bf p}_s,\omega)$ is the phase-encoding function. In the case of plane-wave phase encoding, 
$\alpha({\bf x}_s,{\bf p}_s,\omega)$ is defined as
\begin{eqnarray}
\alpha({\bf x}_s,{\bf p}_s,\omega) = e^{i\omega{\bf p}_s{\bf x}_s},
\end{eqnarray}
where ${\bf p}_s$ is the ray parameter for the source plane waves on the surface. In the case of random phase encoding, 
the phase function is
\begin{eqnarray}
\alpha({\bf x}_s,{\bf p}_s,\omega) = e^{i\gamma({\bf x}_s,{\bf p}_s,\omega)},
\end{eqnarray}
where $\gamma({\bf x}_s,{\bf p}_s,\omega)$ is a random sequence in ${\bf x}_s$ and $\omega$. The parameter ${\bf p}_s$ defines the index
of different realizations of the random sequence \cite[]{yaxun2008b}. 
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 by \cite{etgen2005} and \cite{liu2006} that plane-wave phase-encoding migration, 
by stacking a considerable number of ${\bf p}_s$, produces a migrated image almost identical to the shot-profile migrated image.
If the original shots are well sampled, the number of plane waves required for migration is generally much smaller than the number of
the original shot gathers \cite[]{etgen2005}. Therefore plane-wave source migration is widely used in practice.
Random-phase encoding migration is also an efficient tool, but the random phase function is not very effective in attenuating the crosstalk, 
especially when many sources are simultaneously encoded \cite[]{romero2000,yaxun2008b}. Nevertheless, if many realizations
of the random sequences are used, the final stacked image would also be approximately the same as the shot-profile migrated image.
Therefore, the following relation approximately holds:
\begin{eqnarray}
I({\bf x},{\bf h}) \approx I_{\rm de}({\bf x},{\bf h}).
\end{eqnarray}
That is, with the data-space encoded gathers, we obtain an image similar to that
computed by the more expensive shot-profile migration.
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}.
The tomographic operator ${\bf T}$ and its adjoint ${\bf T}^{\prime}$ can be implemented in a manner similar to that discussed in Appendices B and C by
replacing the original wavefields with the data-space phase encoded wavefields.

\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. 
They are obtained by rotating the original unfocused SODCIGs 
according to the apparent geological dip of the reflector. 
This rotation maintains the velocity information
needed for migration velocity analysis, especially for dipping reflectors \cite[]{biondi2007}.
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 after rotations. 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}.
The tomographic operator ${\bf T}$ and its adjoint ${\bf T}^{\prime}$ can be implemented in a manner similar to that discussed in Appendices B and C, by
replacing the original wavefields with the image-space phase-encoded wavefields.

\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}



%\append{}
\renewcommand{\theequation}{A-\arabic{equation}}
\setcounter{equation}{0}
\section*{APPENDIX A}

This appendix derives the perturbed one-way wave equation with respect to the slowness perturbation. 
Let us start with the one-way wave equation for the source wavefield 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)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_app_wveq_bD}
\end{eqnarray}
We can rewrite the slowness and the source wavefield as follows:
\begin{eqnarray}
s({\bf x}) &=& \widehat{s}({\bf x}) + \Delta s({\bf x}) \label{eq_app_perurb_s}\\
D({\bf x},{\bf x}_s,\omega) &=& \widehat{D}({\bf x},{\bf x}_s,\omega) + \Delta D({\bf x},{\bf x}_s,\omega), \label{eq_app_perturb_D}
\end{eqnarray}
where $\widehat{s}({\bf x})$ and $\widehat{D}({\bf x},{\bf x}_s,\omega)$ are the background slowness and background wavefield, and
$\Delta s({\bf x})$ and $\Delta D({\bf x},{\bf x}_s,\omega)$ are small perturbations in slowness and source wavefield, respectively.
If $\Delta s({\bf x})$ is small, then the square root in the first equation of \ref{eq_app_wveq_bD} can be approximated using Taylor expansion as follows:
\begin{eqnarray}
\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2} \approx \sqrt{\omega ^2 \widehat{s}^2({\bf x})-|{\bf k}| ^2} + 
\frac{\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat{\bf s}}^2({\bf x})}}}. \label{eq_app_perturb_sqrt}
\end{eqnarray}
Substituting Equations \ref{eq_app_perurb_s}, \ref{eq_app_perturb_D} and \ref{eq_app_perturb_sqrt} into Equation \ref{eq_app_wveq_bD} and 
ignoring the second-order terms yield the following linearized one-way wave equation for the perturbed source wavefield: 
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {\widehat s}^2 ({\bf x})-|{\bf k}|^2}\right) \Delta D({\bf x},{\bf x}_s,\omega)=
\frac{- i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}^2({\bf x})}}}{\widehat D}({\bf x},{\bf x}_s,\omega) \\
\Delta D(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right. .
\end{eqnarray}
Similarly, we can also obtain the linearized one-way wave equation for the perturbed receiver wavefield as follows:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {\widehat s}^2 ({\bf x})-|{\bf k}|^2}\right) \Delta U({\bf x},{\bf x}_s,\omega)=
\frac{- i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}^2({\bf x})}}}{\widehat U}({\bf x},{\bf x}_s,\omega) \\
\Delta U(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right..
\end{eqnarray}


%\append{}
\renewcommand{\theequation}{B-\arabic{equation}}
\setcounter{equation}{0}
\section*{APPENDIX B}

This appendix demonstrates a matrix representation of the forward tomographic operator ${\bf T}$.
Let us start with the source wavefield, where the source wavefield ${\bf D}_z$ at depth $z$
is downward continued to depth $z+\Delta z$ by the one-way extrapolator ${\bf E}_z({\bf s}_z)$ as follows:
\begin{eqnarray}
{\bf D}_{z+\Delta z} = {\bf E}_z({\bf s}_z){\bf D}_z, \label{eq_app_down}
\end{eqnarray}
where the one-way extrapolator is defined as follows:
\begin{eqnarray}
{\bf E}_{z}({\bf s}_z) = e^{-ik_z({\bf s}_z)\Delta z} = e^{-i\sqrt{\omega^2{\bf s}_z^2-|{\bf k}|^2}} \label{eq_app_exp}
\end{eqnarray}

The perturbed source wavefield at some depth level can be derived from the background wavefield by a simple
application of the chain rule to equation \ref{eq_app_down}:
\begin{eqnarray}
{\bf \Delta D}_{z+\Delta z} = {\bf E}_z({\widehat{\bf s}}_z){\bf \Delta D}_z + {\bf \Delta E}_z({\widehat{\bf s}}_z)\widehat{{\bf D}}_z, \label{eq_app_mat_dD}
\end{eqnarray}
where $\widehat{{\bf D}}_z$ is the background source wavefield and ${\bf \Delta E}_z$ represents the perturbed extrapolator, 
which can be obtained by a formal linearization 
with respect to slowness of the extrapolator defined in Equation \ref{eq_app_exp}:
\begin{eqnarray}
{\bf E}_z({\bf s}_z) = e^{-ik_z({\bf s}_z)\Delta z}  
&\approx& e^{-i\Delta z {\widehat k_z}} + 
e^{-i\Delta z {\widehat k_z}}\left( -i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z} \right) {\bf \Delta s}_z \nonumber \\
&=&{\bf E}_z(\widehat{{\bf s}}_z) + {\bf E}_z(\widehat{{\bf s}}_z) 
\left( - i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z}\right)  {\bf \Delta s}_z, \label{eq_app_exp_perturb}
\end{eqnarray}
where ${\widehat k}_z=k_z({\widehat {\bf s}}_z)$ and $\widehat{{\bf s}}_z$ is the background slowness at depth $z$. 
From Equation \ref{eq_app_exp_perturb}, the perturbed extrapolator reads as follows:
\begin{eqnarray}
{\bf \Delta E}_z({\widehat {\bf s}}_z) = {\bf E}_z(\widehat{{\bf s}}_z) 
\left( - i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z}\right)  {\bf \Delta s}_z. \label{eq_app_dE}
\end{eqnarray}
Substituting Equation \ref{eq_app_dE} into \ref{eq_app_mat_dD} yields
\begin{eqnarray}
{\bf \Delta D}_{z+\Delta z} = {\bf E}_z({\widehat{\bf s}}_z){\bf \Delta D}_z + 
{\bf E}_z(\widehat{{\bf s}}_z) 
\left( - i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z}\right) \widehat{{\bf D}}_z {\bf \Delta s}_z. \label{eq_app_mat_dD2}
\end{eqnarray}
Let us define a scattering operator ${\bf G}_z$ that interacts with the background wavefield as follows:
\begin{eqnarray}
{\bf G}_z(\widehat{{\bf D}}_z,{\widehat {\bf s}}_z) = \left( - i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z} \right) \widehat{{\bf D}}_z 
= \frac{- i\omega \Delta z}{\sqrt{1-\frac{|{\bf k}|^2}{\omega ^2 {\widehat {\bf s}}_z^2}}} \widehat{{\bf D}}_z.
\end{eqnarray}
Then the perturbed source wavefield for depth $z+\Delta z$ can be rewritten as follows:
\begin{eqnarray}
{\bf \Delta D}_{z+\Delta z} = {\bf E}_z({\widehat {\bf s}}_z){\bf \Delta D}_z + 
{\bf E}_z({\widehat {\bf s}}_z) {\bf G}_z({\widehat {\bf D}}_z,{\widehat {\bf s}}_z) {\bf \Delta s}_z. \label{eq_app_mat_dD3}
\end{eqnarray}
We can further write out the recursive Equation \ref{eq_app_mat_dD3} for all depths in the following matrix form:
\tiny
\begin{eqnarray}
\left( \begin{array}{c}
{\bf \Delta D}_0 \\
{\bf \Delta D}_1 \\
{\bf \Delta D}_2 \\
{\vdots}\\
{\bf \Delta D}_n
\end{array} \right)
&=&
\left( \begin{array}{cccccc}
{\bf 0}   & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf 0}       & {\bf 0}  \\
{\bf E}_0 & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf 0}       & {\bf 0}  \\
{\bf 0}   & {\bf E}_1 & {\bf 0}  & {\cdots} & {\bf 0}       & {\bf 0}  \\
{\vdots}  & {\vdots}  & {\vdots} & {\vdots} & {\vdots}      & {\vdots} \\
{\bf 0}   & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf E}_{n-1} & {\bf 0}
\end{array} \right)
\left( \begin{array}{c}
{\bf \Delta D}_0 \\
{\bf \Delta D}_1 \\
{\bf \Delta D}_2 \\
{\vdots}\\
{\bf \Delta D}_n
\end{array} \right)
+ \nonumber \\
& & \left( \begin{array}{cccccc}
{\bf 0}   & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf 0}        & {\bf 0}\\{\bf E}_0 & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf 0}        & {\bf 0}\\
{\bf 0}   & {\bf E}_1 & {\bf 0}  & {\cdots} & {\bf 0}        & {\bf 0}\\
{\vdots}  & {\vdots}  & {\vdots} & {\vdots} & {\vdots}       & {\bf 0}\\
{\bf 0}   & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf E}_{n-1}  & {\bf 0}
\end{array} \right)
\left( \begin{array}{ccccc}
{\bf G}_0 & {\bf 0}   & {\bf 0}    & {\cdots} & {\bf 0}  \\
{\bf 0}   & {\bf G}_1 & {\bf 0}    & {\cdots} & {\bf 0}  \\
{\bf 0}   & {\bf 0}   & {\bf G}_2  & {\cdots} & {\bf 0}  \\
{\vdots}  & {\vdots}  & {\vdots}   & {\vdots} & {\vdots} \\
{\bf 0}   & {\bf 0}   & {\bf 0}    & {\cdots} & {\bf G}_n
\end{array} \right)
\left( \begin{array}{c}
{\bf \Delta s}_0 \\
{\bf \Delta s}_1 \\
{\bf \Delta s}_2 \\
{\vdots}\\
{\bf \Delta s}_n
\end{array} \right), \nonumber
\end{eqnarray}
\normalsize
or in a more compact notation,
\begin{eqnarray}
{\bf \Delta D} = {\bf E}({\widehat {\bf s}}){\bf \Delta D} + 
{\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf D}},{\widehat {\bf s}}){\bf \Delta s}. \label{eq_app_compact_dD}
\end{eqnarray}
The solution of Equation \ref{eq_app_compact_dD} can be formally written as follows:
\begin{eqnarray}
{\bf \Delta D} = \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} 
{\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf D}},{\widehat {\bf s}}){\bf \Delta s}. \label{eq_app_solution_dD}
\end{eqnarray}

Similarly, the perturbed receiver wavefield satisfies the following recursive relation:
\begin{eqnarray}
{\bf \Delta U}_{z+\Delta z} = {\bf E}_z({\widehat {\bf s}}_z){\bf \Delta U}_z + 
{\bf E}_z({\widehat {\bf s}}_z) {{\bf G}_z({\widehat {\bf U}}_z,{\widehat {\bf s}}_z)} {\bf \Delta s}_z,
\end{eqnarray}
where ${\bf G}_z({\widehat {\bf U}}_z,{\widehat {\bf s}}_z)$ is the scattering operator, which interacts with the background receiver wavefield as follows:
\begin{eqnarray}
{\bf G}_z({\widehat {\bf U}}_z,{\widehat {\bf s}}_z) = 
\left( - i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z} \right) {\widehat {\bf U}}_z 
= \frac{- i\omega \Delta z}{\sqrt{1-\frac{|{\bf k}|^2}{\omega ^2 {\widehat {\bf s}}_z^2}}} {\widehat {\bf U}}_z. \label{eq_app_dU}
\end{eqnarray}
We can also write out the recursive Equation \ref{eq_app_dU} for all depth levels in the following matrix form:
\tiny
\begin{eqnarray}
\left( \begin{array}{c}
{\bf \Delta U}_0 \\
{\bf \Delta U}_1 \\
{\bf \Delta U}_2 \\
{\vdots}\\
{\bf \Delta U}_n
\end{array} \right)
&=&
\left( \begin{array}{cccccc}
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                  & {\bf 0}  \\
{\bf E}_0            & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                  & {\bf 0}  \\
{\bf 0}              & {\bf E}_1            & {\bf 0}  & {\cdots} & {\bf 0}                  & {\bf 0}  \\
{\vdots}             & {\vdots}             & {\vdots} & {\vdots} & {\vdots}                 & {\vdots} \\
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf E}_{n-1}            & {\bf 0}
\end{array} \right)
\left( \begin{array}{c}
{\bf \Delta U}_0 \\
{\bf \Delta U}_1 \\
{\bf \Delta U}_2 \\
{\vdots}\\
{\bf \Delta U}_n
\end{array} \right)
+ \nonumber \\
& & \left( \begin{array}{cccccc}
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                   & {\bf 0}\\
{\bf E}_0            & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                   & {\bf 0}\\
{\bf 0}              & {\bf E}_1            & {\bf 0}  & {\cdots} & {\bf 0}                   & {\bf 0}\\
{\vdots}             & {\vdots}             & {\vdots} & {\vdots} & {\vdots}                  & {\bf 0}\\
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf E}_{n-1}             & {\bf 0}
\end{array} \right)
\left( \begin{array}{ccccc}
{\bf G}_0 & {\bf 0}   & {\bf 0}    & {\cdots} & {\bf 0}  \\
{\bf 0}   & {\bf G}_1 & {\bf 0}    & {\cdots} & {\bf 0}  \\
{\bf 0}   & {\bf 0}   & {\bf G}_2  & {\cdots} & {\bf 0}  \\
{\vdots}  & {\vdots}  & {\vdots}   & {\vdots} & {\vdots} \\
{\bf 0}   & {\bf 0}   & {\bf 0}    & {\cdots} & {\bf G}_n
\end{array} \right)
\left( \begin{array}{c}
{\bf \Delta s}_0 \\
{\bf \Delta s}_1 \\
{\bf \Delta s}_2 \\
{\vdots}\\
{\bf \Delta s}_n
\end{array} \right), \nonumber
\end{eqnarray}
\normalsize
or in a more compact notation,
\begin{eqnarray}
{\bf \Delta U} = {\bf E}({\widehat {\bf s}}){\bf \Delta U} + 
{\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf U}},{\widehat {\bf s}}){\bf \Delta s}. \label{eq_app_dU2}
\end{eqnarray}
The solution of Equation \ref{eq_app_dU2} can be formally written as follows:
\begin{eqnarray}
{\bf \Delta U} = \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} {\bf E}({\widehat {\bf s}})
{\bf G}({\widehat {\bf U}},{\widehat {\bf s}}){\bf \Delta s}. \label{eq_app_solution_dU}
\end{eqnarray}

With the background wavefields and the perturbed wavefields, the perturbed image can be obtained as follows: 
\tiny
\begin{eqnarray}
\left( \begin{array}{c}
{\bf \Delta I}_0 \\
{\bf \Delta I}_1 \\
{\bf \Delta I}_2 \\
{\vdots}\\
{\bf \Delta I}_n
\end{array} \right)
=
\left( \begin{array}{ccccc}
{\widehat {\bf U}}_0 & {\bf 0}              & {\bf 0}               & {\cdots} & {\bf 0}  \\
{\bf 0}              & {\widehat {\bf U}}_1 & {\bf 0}               & {\cdots} & {\bf 0}  \\
{\bf 0}              & {\bf 0}              & {\widehat {\bf U}}_2  & {\cdots} & {\bf 0}  \\
{\vdots}             & {\vdots}             & {\vdots}              & {\vdots} & {\vdots} \\
{\bf 0}              & {\bf 0}              & {\bf 0}               & {\cdots} & {\widehat {\bf U}}_n
\end{array} \right)
\left( \begin{array}{c}
{\bf \Delta D}_0\\
{\bf \Delta D}_1\\
{\bf \Delta D}_2\\
{\vdots}\\
{\bf \Delta D}_n
\end{array} \right)
+
\left( \begin{array}{ccccc}
{\widehat{\bf D}}_0            & {\bf 0}                        & {\bf 0}                         & {\cdots} & {\bf 0}  \\
{\bf 0}                        & {\widehat{\bf D}}_1            & {\bf 0}                         & {\cdots} & {\bf 0}  \\
{\bf 0}                        & {\bf 0}                        & {\widehat{\bf D}}_2             & {\cdots} & {\bf 0}  \\
{\vdots}                       & {\vdots}                       & {\vdots}                        & {\vdots} & {\vdots} \\
{\bf 0}                        & {\bf 0}                        & {\bf 0}                         & {\cdots} & {\widehat{\bf D}}_n
\end{array} \right)
\left( \begin{array}{c}
{{\bf \Delta U}_0}\\
{{\bf \Delta U}_1}\\
{{\bf \Delta U}_2}\\
{\vdots}\\
{{\bf \Delta U}_n}
\end{array} \right), \nonumber
\end{eqnarray}
\normalsize
or in a more compact notation,
\begin{eqnarray}
{\bf \Delta I} = {\rm diag}\left({\widehat {\bf U}}\right){\bf \Delta D} + 
{\rm diag}\left( {\widehat{\bf D}}\right){\bf \Delta U}. \label{eq_app_compact_dI}
\end{eqnarray}
Substituting Equations \ref{eq_app_solution_dD} and \ref{eq_app_solution_dU} into Equation \ref{eq_app_compact_dI} yields
\begin{eqnarray}
{\bf \Delta I} &=& \left( {\rm diag}\left( {\widehat {\bf U}}\right)
                    \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} 
                    {\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf D}},{\widehat {\bf s}}) \right. + \nonumber \\
           & &    \left.  {\rm diag}\left({\widehat {\bf D}} \right)
                    \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} 
                    {\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf U}},{\widehat {\bf s}})
             \right){\bf \Delta s}, \label{eq_app_compact_dI2}
\end{eqnarray}
from which we can read the forward tomographic operator ${\bf T}$ as follows: 
\begin{eqnarray}
{\bf T} &=& {\rm diag}\left( {\widehat {\bf U}}\right)
                     \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} 
                    {\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf D}},{\widehat {\bf s}}) + \nonumber \\
           & &     {\rm diag}\left({\widehat {\bf D}} \right)
                    \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} 
                    {\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf U}},{\widehat {\bf s}}).
\end{eqnarray}


%\append{}
\renewcommand{\theequation}{C-\arabic{equation}}
\setcounter{equation}{0}
\section*{APPENDIX C}

This appendix demonstrates a matrix representation of the adjoint tomographic operator ${\bf T}^{'}$.
Since the slowness perturbation ${\bf \Delta s}$ is linearly related to the perturbed wavefields, ${\bf \Delta D}$ and ${\bf \Delta U}$,
to obtain the back-projected slowness perturbation, we first must get the back-projected perturbed wavefields from the perturbed image ${\bf \Delta I}$.
From Equation \ref{eq_app_compact_dI}, the back-projected perturbed source and receiver wavefields are obtained as follows:
\begin{eqnarray}
{\bf \Delta D} = \overline{{\rm diag}\left( {\widehat {\bf U}} \right)} {\bf \Delta I} \label{eq_app_adjoint_dD}
\end{eqnarray}
and
\begin{eqnarray}
{\bf \Delta U} = \overline{{\rm diag}\left( {\widehat {\bf D}} \right)} {\bf \Delta I}. \label{eq_app_adjoint_dU}
\end{eqnarray}
Then the adjoint equations of Equations \ref{eq_app_solution_dD} and \ref{eq_app_solution_dU} are 
used to get the back-projected slowness perturbation ${\bf \Delta s}$.
Let us first look at the adjoint equation of Equation \ref{eq_app_solution_dD}, which can be written as follows:
\begin{eqnarray}
{\bf \Delta s}_{D} = {\bf G}^{\prime}({\widehat {\bf D}},{\widehat {\bf s}}){\bf E}^{\prime}({\widehat {\bf s}})
                     \left( {\bf 1} - {\bf E}^{\prime}({\widehat {\bf s}})\right)^{-1}{\bf \Delta D}.
\end{eqnarray}
We can define a temporary wavefield ${\bf \Delta P}_{D}$ that satisfies the following equation:
\begin{eqnarray}
{\bf \Delta P}_{D} = {\bf E}^{\prime}({\widehat {\bf s}}) \left( {\bf 1} - {\bf E}^{\prime}({\widehat {\bf s}})\right)^{-1}{\bf \Delta D}.
\end{eqnarray}
After some simple algebra, the above equation can be rewritten as follows:
\begin{eqnarray}
{\bf \Delta P}_{D} = {\bf E}^{\prime}({\widehat {\bf s}}) {\bf \Delta P}_{D} + {\bf E}^{\prime}({\widehat {\bf s}}) {\bf \Delta D}. \label{eq_app_dPD}
\end{eqnarray}
Substituting Equation \ref{eq_app_adjoint_dD} into equation \ref{eq_app_dPD} yields
\begin{eqnarray}
{\bf \Delta P}_{D} = {\bf E}^{\prime}({\widehat {\bf s}}) {\bf \Delta P}_{D} + 
{\bf E}^{\prime}({\widehat {\bf s}}) \overline{{\rm diag}\left( {\widehat {\bf U}} \right)} {\bf \Delta I}.
\end{eqnarray}
Therefore, ${\bf \Delta P}_{D}$ can be obtained by recursive upward continuation, 
where ${\bf \Delta D} = \overline{{\rm diag}\left( {\widehat {\bf U}} \right)} {\bf \Delta I}$ 
serves as the initial condition. 
The back-projected slowness perturbation from the perturbed source wavefield is then obtained by applying the adjoint of the scattering
operator ${\bf G}({\widehat {\bf D}},{\widehat {\bf s}})$ to the wavefield ${\bf \Delta P}_{D}$ as follows:
\begin{eqnarray}
{\bf \Delta s}_{D} = {\bf G}^{\prime}({\widehat {\bf D}},{\widehat {\bf s}}) {\bf \Delta P}_{D}.
\end{eqnarray}

Similarly, the adjoint equation of Equation \ref{eq_app_solution_dU} reads as follows:
\begin{eqnarray}
{\bf \Delta s}_{U} = {\bf G}^{\prime}({\widehat {\bf U}},{\widehat {\bf s}}){\bf E}({\widehat {\bf s}})^{\prime}
                     \left( {\bf 1} - {\bf E}({\widehat {\bf s}})^{\prime}\right)^{-1}{\bf \Delta U}.
\end{eqnarray}
We can also define a temporary wavefield ${\bf \Delta P}_{U}$ that satisfies the following equation:
\begin{eqnarray}
{\bf \Delta P}_{U} = {\bf E}({\widehat {\bf s}})^{\prime} \left( {\bf 1} - 
{\bf E}({\widehat {\bf s}})^{\prime}\right)^{-1}{\bf \Delta U}.
\end{eqnarray}
After rewriting it, we get the following recursive form:
\begin{eqnarray}
{\bf \Delta P}_{U} &=& {\bf E}({\widehat {\bf s}})^{\prime} {\bf \Delta P}_{U}+
  {\bf E}({\widehat {\bf s}})^{\prime} {\bf \Delta U} \nonumber \\
&=& {\bf E}({\widehat {\bf s}})^{\prime} {\bf \Delta P}_{U}+
  {\bf E}({\widehat {\bf s}})^{\prime} \overline{{\rm diag}\left( {\widehat {\bf D}} \right)} {\bf \Delta I}.
\end{eqnarray}
The back-projected slowness perturbation from the perturbed receiver wavefield is then obtained by applying the adjoint of
the scattering operator ${\bf G}({\widehat{\bf U}},{\widehat{\bf s}})$  to the wavefield ${\bf \Delta P}_{U}$ as follows:
\begin{eqnarray}
{\bf \Delta s}_{U} = {\bf G}^{\prime}({\widehat {\bf U}},{\widehat {\bf s}}) {\bf \Delta P}_{U}.
\end{eqnarray}
The total back-projected slowness perturbation is obtained by adding ${\bf \Delta s}_D$ and ${\bf \Delta s}_U$ together:
\begin{eqnarray}
{\bf \Delta s} = {\bf \Delta s}_D + {\bf \Delta s}_U.
\end{eqnarray}
