\title{Image-space wave-equation tomography in the generalized source domain}
\author{Yaxun Tang, Claudio Cardoso 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 either generated 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 these new domains.
We compare the gradients of the image-space misfit 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 similar gradient as 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 the well known
triplication model with a strong Gaussian anormaly in the subsurface.
\end{abstract}

\section{introduction}
Wave-equation tomography has the potential to accurately estimate the velocity model in complex geological senarios 
where ray-based traveltime tomography often fails. Wave-equation based tomography uses band-limited wavefields instead of wide-band rays as carriers
of information, thus it is robust even in the presence of strong velocity contrast and immune from multi-pathing issues. 
Generally speaking, wave-equation tomography can be classified into two different categaries based on the domian where it minimizes the residual. 
The domain can either be the data space or the image space. For the data-space based approach, it directly compares the modeled waveform with the recorded
waveform, and is widely known as the 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 based approach, also known as image-space
wave-equation tomography, minimizes the residual 
in the image domain obtained after migration. The image domain is often much simpler than the data domain, 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
tompgraphy has the potential to mitigate some of the difficulties that we encounter in the data-space based approach. 
Another advantage of the image-space based approach is that the more efficient one-way wave-equation extrapolator can be utilized, whereas
in waveform inversion, the one-way propagator is difficult (if not imposible) to find its usage because of its inability to model the multiple arrivales, though
some tweaks can be employed so that the one-way propagater can be applied in the case of turning wave tomography \cite[]{shragge2007}.

However, despite of its theoretical advantages, 
image-space wave-equation tomography suffers from the computational cost.
Each iteration of tomographic velocity updating is computationally
expensive and often converges slowly. Practial 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 original shot-profile
domain to the generalized source domain, where a smaller number of synthesized shot gathers make the tomographic velocity
update substentially faster.
The generalized source domain can be obtained either by data-space phase encoding
or image-space phase encoding. 
By data-space phase encoding, we mean that the synthesized shot gathers are obtained by linear combination of
the original shot gather 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.
By image-space phase encoding, we mean that the synthesized gathers are obtained by prestack exploding reflector modeling \cite[]{biondi2006,cardoso2008}.
where several subsurface-offset-domain common-image gathers (SODCIGs) are simultaneously demigrated to generate area source and area receiver gathers. 
To attenuate the cross-talk, the SODCIGs
have to be encoded, e.g., by random phase encoding. This encoding process is done in the image space, so 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 brieflly review the theory of image-space wave-equation tomography, 
we review the formulation of the objective function and its gradient; then we discuss how to evaluate the forward
tomographic operator and its adjoint in the orignal
shot-profile domain, the latter is an important component in computing the gradient; then we 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, 
it is based on an optimization 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 semblence in the ADCIGs, 
reflection angle stacking power (in this case we have to maximize the residual field, or minimize the minus stacking power), and etc..
Here we follow a similar definition as 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 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; 
${\bf K}[\Delta \rho]$ is the differential residual migration operator defined as follws:
\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 angle and geological dips \cite[]{biondi2008}.
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 O}$ is the DSO operator either in the subsurface offset domain or in the angle domain. 
Subsurface offset domain DSO focuses the energy at zero offset while angle domain DSO flattens the ADCIGs. 

In general, if we choose $\ell_2$ norm, the 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 Quasi-Newton method and Conjugate Gradient method can be used to minimize the objective function $J$, 
which is non-linear because the image ${\bf I}$ is a non-linear function of the slowness ${\bf s}$. 
To use the gradient based methods, we have to provide the
gradient of $J$ with respect to the slowness ${\bf s}$.
The gradient reads
\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 $^{\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} yields
\begin{eqnarray}
{\bf \nabla} J_{\rm DSO} = \Re \left( \left(\frac{\partial {\bf I}}{\partial {\bf s}}\right)^{\prime}{\bf O}^{\prime}{\bf O}{\bf I}\right)
\end{eqnarray}

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 performed on the background
image $\widehat{\bf I}$ instead of ${\bf I}$, i.e., 
\begin{eqnarray}
{\bf F}(\widehat{\bf I})=\widehat{\bf I}+{\bf K}[\widehat{\Delta \rho}]\widehat{\bf I},
\end{eqnarray}
With this assumption, we get the "classic" WEMVA gradient as follows:
\begin{eqnarray}
\nabla J_{\rm WEMVA} = \Re \left( - \left(\frac{\partial {\bf I}}{\partial {\bf s}}\right)^{\prime}{\bf K}[\widehat{\Delta \rho}]\widehat{\bf I} \right)
\end{eqnarray}
The complete WEMVA gradient can also be derived following the method described by \cite{biondi2008}. 

No matter which  gradient we choose to back-project the slowness perturbeation, we have to evaluate the adjoint of the linear operator
$\frac{\partial{\bf I}}{\partial{\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 a background slowness ${\widehat {\bf s}}$ as follows:
\begin{eqnarray}
{\bf I} = {\widehat {\bf I}} + \frac{\partial{\bf I}}{{\partial{\bf s}}}({\bf s}-{\widehat {\bf s}}) + \cdots,
\end{eqnarray} 
where ${\widehat{\bf I}}$ is the background image computed with the background slowness ${\widehat {\bf s}}$. 
Keeping only the zero and first order terms, we get the linearized operator $\frac{\partial{\bf I}}{\partial{\bf s}}$ defined as follows:
\begin{eqnarray}
{\bf \Delta I} = \frac{\partial{\bf I}}{\partial{\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} = \frac{\partial{\bf I}}{\partial{\bf s}}$ is the wave-equation tomographic operator.
The tomographic operator can either be defined in the source and receiver domain \cite[]{sava2004} or in the shot-profile domain \cite[]{shen2004}.
In next section we follow a similar approach disscussed 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 
generalized source domains.

\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 freqnency 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 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 at some image point 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 
linearized (with respect to slowness) one-way wave equations, respectively(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}({\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}({\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}

Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} are linear with respect to slowness perturbation, 
thus they define linear mappings from the slowness perturbation $\Delta s({\bf x})$ to
the perturbed source wavefield $\Delta D({\bf x},{\bf x}_s,\omega)$ and the perturbed receiver wavefield $\Delta U({\bf x},{\bf x}_s,\omega)$, respectively.
Therefore, recursively solving equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} give 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 $I({\bf x},{\bf h})$, where
the background source and receiver wavefields are precomputed by solving the recursive 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 compute 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 \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 \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 B 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 as that of a shot-profile migration, because to compute the perturbed wavefields, the background
wavefields are required. As 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 nowadays computer resources.
To save the cost and to make this powerful method more practical, we extend the wave-eqaution tomography theory to
generalized source domains, or phase encoded source domain. We discuss two different strategies generate 
the generalized shot gathers, i.e., the data-space phase encoding and the image-space phase encoding.

\subsection{Data-space encoded wavefields}
We encode the source wavefield and the corresponding receiver wavefield
\begin{eqnarray}
{\widetilde D}({\bf x},{\bf p}_s,\omega) &=& \sum_{{\bf x}_s} D({\bf x},{\bf x}_s,\omega) \alpha({\bf x}_s,{\bf p}_s,\omega) \\
{\widetilde U}({\bf x},{\bf p}_s,\omega) &=& \sum_{{\bf x}_s} U({\bf x},{\bf x}_s,\omega) \alpha({\bf x}_s,{\bf p}_s,\omega)
\end{eqnarray}
where $\alpha({\bf x}_s,{\bf p}_s,\omega)$ is some 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.
Because of the linearity of the wave equation, the en

Imaging in the generallized source domain
\begin{eqnarray}
I({\bf x},{\bf h}) = \sum_{{\bf p}_s}\sum_{\omega} {\widetilde D}({\bf x}-{\bf h},{\bf p}_s,\omega) {\widetilde U}({\bf x}+{\bf h},{\bf p}_s,\omega)
\end{eqnarray}



\subsection{Image-space encoded wavefields}
As introduced by \cite {biondi2006,biondi2007}, the prestack exploding reflector modeling aims to generate a smaller dataset to speed up iterations of migration velocity analysis. The method synthesizes areal source ($D$) and areal receiver ($U$) data at the surface by recursively upward propagating source and receiver wavefields, having prestack images, $I_D$ and $I_U$, obtained by wave-equation migration with an inaccurate velocity, ${\hat s}({\bf x})$, as source terms
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\right) D({\bf x},\omega) = f_s(\omega)I_D({\bf x},{\bf h}) \\
D(x,y,z=z_{max},\omega) = 0 \end{array} \right.,
\label{eq:perm1}
\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) U({\bf x},\omega) = f_s(\omega)I_U({\bf x},{\bf h}) \\
U(x,y,z=z_{max},\omega) = 0 \end{array} \right.,
\label{eq:perm2}
\end{eqnarray}
where $f_s(\omega)$ is a frequency dependent source wavelet, such that the source terms are not strictly impulsive. The rationale of using different source terms for the source and the receiver propagation is due to the need of rotating the image when the background velocity is inaccurate and in the presence of an apparent geological dip, as \cite{biondo} shows.
\par
The reduction of the dataset size to be migrated is achieved by combining the modeled wavefields into super areal data. To avoid crosstalk, due to correlation of unrelated wavefields during migration, phase encoding is applied in the modeling step, according to
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\right){\check D}({\bf x},{\bf p_m};\omega) = f_s(\omega)\alpha({\bf x},{\bf p_m};\omega)I_D({\bf x},{\bf h}) \\
{\check D}(x,y,z=z_{max},{\bf p_m};\omega) = 0 \end{array} \right.,
\label{eq:perm3}
\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){\check U}({\bf x},{\bf p_m};\omega) = f_s(\omega)\alpha({\bf x},{\bf p_m};\omega)I_U({\bf x},{\bf h}) \\
{\check U}(x,y,z=z_{max},{\bf p_m};\omega) = 0 \end{array} \right.,
\label{eq:perm4}
\end{eqnarray}
where $\alpha({\bf x},{\bf p_m};\omega) = e^{i\gamma({\bf x},{\bf p_m};\omega)}$. \cite{guerra2008a} showed that random phases can be used to encode the modeling experiments such that the crosstalk is dispersed throughout the image. In this case, $\gamma$ is a random function varying in the interval $(-\pi,\pi)$. In a companion paper, \cite{guerra2008b} use random phase functions with particular correlation properties to encode the modeling experiments. ${\check D}$ and ${\check U}$ defines, respectively, the synthesized source and receiver wavefields at the surface represented by super areal data. The vector ${\bf p}_m$ represents the modeling realizations. 
\par
For the areal shot migration of the synthesized data, both source and receiver wavefields are downward continued 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){\widehat D}({\bf x},{\bf p}_m,\omega) = 0 \\
{\widehat D}(x,y,z=0,{\bf p}_m,\omega) = {\check D}(x,y,z=0,{\bf p}_m,\omega) \end{array} \right.,
\label{eq:permD}
\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){\widehat U}({\bf x},{\bf p}_m,\omega) = 0 \\
{\widehat U}(x,y,z=0,{\bf p}_m,\omega) = {\check U}(x,y,z=0,{\bf p}_m,\omega) \end{array} \right.,
\label{eq:permU}
\end{eqnarray}
where ${\widehat D}$ and ${\widehat U}$ are the downward propagated background source and receiver wavefields, respectively. 
\par
The background image is computed by applying the cross-correlation imaging condition according to:
\begin{eqnarray}
{\widehat I}({\bf x},{\bf h}) = \sum_{{\bf p}_m}\sum_{\omega} \overline{{\widehat D}({\bf x}-{\bf h},{\bf p}_m,\omega)} {\widehat U}({\bf x}+{\bf h},{\bf p}_m,\omega).
\label{eq:permIm}
\end{eqnarray}
Ideally, if the phase encoding scheme is effective, the image ${\widehat I({\bf x},{\bf h})}$ shows no crosstalk.
\par
Similarly, to the shot-profile case, a perturbed image, $\Delta I$, can be obtained by applying the chain rule to equation~\ref{eq:permIm}:
\begin{eqnarray}
\Delta I({\bf x},{\bf h})  &=& 
\sum_{{\bf p}_m}\sum_{\omega} 
\left( \overline{\Delta D({\bf x}-{\bf h},{\bf p}_m,\omega)} {\widehat U}({\bf x}+{\bf h},{\bf p}_m,\omega) + \right. \nonumber \\
& & \left. \overline{{\widehat D}({\bf x}-{\bf h},{\bf p}_m,\omega)} \Delta U ({\bf x}+{\bf h},{\bf p}_m,\omega) \right),
\label{eq:permDIm}
\end{eqnarray}
where $\Delta D$ and $\Delta U$ are the perturbed source and receiver wavefields, respectively. The perturbed wavefields satisfy the one-way wave equation linearized with respect to slowness:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\right){\Delta D}({\bf x},{\bf p}_m,\omega) = \frac{+i\omega \Delta {\hat s}(\bf x)}{\sqrt{1-\frac{\left |{\bf k} \right |^2}{\omega ^2 {\hat s}^2({\bf x})}}}{\widehat D}({\bf x},{\bf p}_m,\omega) \\
{\Delta D}(x,y,z=0,{\bf p}_m,\omega) = 0 \end{array} \right.,
\label{eq:permDD}
\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){\Delta U}({\bf x},{\bf p}_m,\omega) = \frac{-i\omega \Delta {\hat s}(\bf x)}{\sqrt{1-\frac{\left |{\bf k} \right |^2}{\omega ^2 {\hat s}^2({\bf x})}}}{\widehat U}({\bf x},{\bf p}_m,\omega) \\
{\Delta U}(x,y,z=0,{\bf p}_m,\omega) = 0 \end{array} \right..
\label{eq:permDU}
\end{eqnarray}
The slowness perturbation is obtained by applying the adjoint of  the tomographic operator to the image perturbation. To compute the adjoint of the tomographic operator, we firstly obtain the perturbed wavefields by applying the adjoint of the imaging operator
\begin{eqnarray}
\Delta D({\bf x},{\bf p}_m,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h}){\widehat U}({\bf x}+{\bf h},{\bf p}_m,\omega) \\
\Delta U({\bf x},{\bf p}_m,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h}){\widehat D}({\bf x}-{\bf h},{\bf p}_m,\omega).
\end{eqnarray}
The perturbed image, $\Delta I$, is computed according to equation \ref{eq:permDIm}. Then the perturbed wavefields are upward propagated by solving the adjoint equations of \ref{eq:permDD} and \ref{eq:permDU}.
Notice that, as stated in equations~\ref{eq:permDD} and \ref{eq:permDU}, scattered wavefields are injected at every depth step prior to the upward propagation. The scattered wavefields results from the interaction of the scattering operator  with the background wavefields, $\widehat D$ and $\widehat U$. 


\section{Numerical examples}


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

\append{}
This appendix derives the linearized one-way wave equation

\append{}
This appendix demonstrates a matrix representation of the forward tomographic operator ${\bf T}$
\begin{eqnarray}
{\bf D}_{z+\Delta z} = {\bf E}_z({\bf s}_z){\bf \Delta D}_z, \label{eq_app_down}
\end{eqnarray}
where
\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){\bf D}_z, \label{eq_app_mat_dD}
\end{eqnarray}
where ${\bf \Delta E}_z$ represents the perturbed extrapolator, it 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
\end{eqnarray}
where ${\widehat k}_z=k_z({\widehat {\bf s}}_z)$. So the perturbed extrapolator reads
\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}
Substitute equation \ref{eq_app_dE} into \ref{eq_app_mat_dD}, we get
\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) {\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({\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) {\bf D}_z 
= \frac{i\omega \Delta z}{\sqrt{1-\frac{|{\bf k}|^2}{\omega ^2 {\widehat {\bf s}}_z^2}}} {\bf D}_z
\end{eqnarray}
Then the perturbed source wavefield for depth level $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 depth levels 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 receiver wavefield satisfy the following recursive relation:
\begin{eqnarray}
{\bf U}_{z+\Delta z} = \overline{{\bf E}_z({\widehat {\bf s}}_z)}{\bf \Delta U}_z, \label{eq_app_down_U}
\end{eqnarray}
where the over line stands for the complex conjugate, thus
\begin{eqnarray}
\overline{{\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}\Delta z}.
\end{eqnarray}
Apply the chain rule to euqation \ref{eq_app_down_U}, we obtain the recursive equation for the scattered receiver wavefield:
\begin{eqnarray}
{\bf \Delta U}_{z+\Delta z} = \overline{{\bf E}_z({\widehat {\bf s}}_z)}{\bf \Delta U}_z - 
\overline{{\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 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 further 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}  \\
\overline{{\bf E}_0} & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                  & {\bf 0}  \\
{\bf 0}              & \overline{{\bf E}_1} & {\bf 0}  & {\cdots} & {\bf 0}                  & {\bf 0}  \\
{\vdots}             & {\vdots}             & {\vdots} & {\vdots} & {\vdots}                 & {\vdots} \\
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & \overline{{\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}\\
\overline{{\bf E}_0} & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                   & {\bf 0}\\
{\bf 0}              & \overline{{\bf E}_1} & {\bf 0}  & {\cdots} & {\bf 0}                   & {\bf 0}\\
{\vdots}             & {\vdots}             & {\vdots} & {\vdots} & {\vdots}                  & {\bf 0}\\
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & \overline{{\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} = \overline{{\bf E}({\widehat {\bf s}})}{\bf \Delta U} - 
\overline{{\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} - \overline{{\bf E}({\widehat {\bf s}})} \right)^{-1} \overline{{\bf E}({\widehat {\bf s}})}
{\bf G}({\widehat {\bf U}},{\widehat {\bf s}}){\bf \Delta s}. \label{eq_app_solution_dU}
\end{eqnarray}

Therefore, 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}
\overline{{\bf \Delta D}_0}\\
\overline{{\bf \Delta D}_1}\\
\overline{{\bf \Delta D}_2}\\
{\vdots}\\
\overline{{\bf \Delta D}_n}
\end{array} \right)
+
\left( \begin{array}{ccccc}
{\overline{\widehat{\bf D}}_0} & {\bf 0}                        & {\bf 0}                         & {\cdots} & {\bf 0}  \\
{\bf 0}                        & \overline{{\widehat{\bf D}}_1} & {\bf 0}                         & {\cdots} & {\bf 0}  \\
{\bf 0}                        & {\bf 0}                        & \overline{{\widehat{\bf D}}_2}  & {\cdots} & {\bf 0}  \\
{\vdots}                       & {\vdots}                       & {\vdots}                        & {\vdots} & {\vdots} \\
{\bf 0}                        & {\bf 0}                        & {\bf 0}                         & {\cdots} & \overline{{\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)\overline{{\bf \Delta D}} + 
\overline{{\rm diag}\left( {\widehat{\bf D}}\right)}{\bf \Delta U}. \label{eq_app_compact_dI}
\end{eqnarray}
After substituting equations \ref{eq_app_solution_dD} and \ref{eq_app_solution_dU} into equation \ref{eq_app_compact_dI} 
\begin{eqnarray}
{\bf \Delta I} &=& \left( {\rm diag}\left( {\widehat {\bf U}}\right)
                    \overline{ \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.  \overline{{\rm diag}\left({\widehat {\bf D}} \right)}
                    \left( {\bf 1} - \overline{{\bf E}({\widehat {\bf s}})} \right)^{-1} 
                    \overline{{\bf E}({\widehat {\bf s}})}{\bf G}({\widehat {\bf U}},{\widehat {\bf s}})
             \right){\bf \Delta s}, \label{eq_app_compact_dI2}
\end{eqnarray}
we finally get the following matrix representation of the forward tomographic operator ${\bf T}$: 
\begin{eqnarray}
{\bf T} &=& {\rm diag}\left( {\widehat {\bf U}}\right)
                    \overline{ \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} 
                    {\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf D}},{\widehat {\bf s}})} + \nonumber \\
           & &     \overline{{\rm diag}\left({\widehat {\bf D}} \right)}
                    \left( {\bf 1} - \overline{{\bf E}({\widehat {\bf s}})} \right)^{-1} 
                    \overline{{\bf E}({\widehat {\bf s}})}{\bf G}({\widehat {\bf U}},{\widehat {\bf s}})
\end{eqnarray}


\append{}
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 need to 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} = {\rm diag}\left( {\widehat {\bf U}} \right) {\bf \Delta I} \label{eq_app_adjoint_dD}
\end{eqnarray}
and
\begin{eqnarray}
{\bf \Delta U} = {\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}
\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 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}
Substitute equation \ref{eq_app_adjoint_dD} into equation \ref{eq_app_dPD}, we get:
\begin{eqnarray}
{\bf \Delta P}_{D} = {\bf E}^{\prime}({\widehat {\bf s}}) {\bf \Delta P}_{D} + 
{\bf E}^{\prime}({\widehat {\bf s}}) {\rm diag}\left( {\widehat {\bf U}} \right) {\bf \Delta I}
\end{eqnarray}
Therefore, ${\bf \Delta P}_{D}$ can be obtained by recursively upward continuation, where ${\bf \Delta D} = {\rm diag}\left( {\widehat {\bf U}} \right) {\bf \Delta I}$ 
serves as the boundary condition. 
The back-projected slowness perturbation from the perturbed source wavefield is then obtained by applying the ajoint of the scattering
operator ${\bf G}({\widehat {\bf D}},{\widehat {\bf s}})$ to the wavefield ${\bf \Delta P}_{D}$
\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
\begin{eqnarray}
{\bf \Delta s}_{U} = -{\bf G}^{\prime}({\widehat {\bf U}},{\widehat {\bf s}})\left(\overline{{\bf E}({\widehat {\bf s}})}\right)^{\prime}
                     \left( {\bf 1} - \left(\overline{{\bf E}({\widehat {\bf s}})}\right)^{\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} = \left(\overline{{\bf E}({\widehat {\bf s}})}\right)^{\prime} \left( {\bf 1} - 
\left(\overline{{\bf E}({\widehat {\bf s}})}\right)^{\prime}\right)^{-1}{\bf \Delta U}
\end{eqnarray}
After rewriting it, we get the following recursive form
\begin{eqnarray}
{\bf \Delta P}_{U} &=& \left(\overline{{\bf E}({\widehat {\bf s}})}\right)^{\prime} {\bf \Delta P}_{U}+
  \left(\overline{{\bf E}({\widehat {\bf s}})}\right)^{\prime} {\bf \Delta U} \nonumber \\
&=& \left(\overline{{\bf E}({\widehat {\bf s}})}\right)^{\prime} {\bf \Delta P}_{U}+
  \left(\overline{{\bf E}({\widehat {\bf s}})}\right)^{\prime} {\rm diag}\left( {\widehat {\bf D}} \right) {\bf \Delta I}.
\end{eqnarray}
And 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}$
\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}
