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

%\maketitle

\begin{abstract}
In this paper, we generalize wave-equation migration velocity analysis to composite shot domain. 
Any linear composite shot gather can be used for WEMVA. One particular case is the plane-wave shot gathers. 
I will demonstate some theoretical advantages of doing WEMVA in composite shot domain, 
\end{abstract}

\section{introduction}

\section{wave-equation tomography}
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..
For simplicity, we define a linear operator ${\bf K}$ that measures the image perturbation as follows:
\begin{eqnarray}
{\bf \Delta I} = ({\bf 1} - {\bf K}){\bf I},
\end{eqnarray}
where ${\bf 1}$ is the identity operator.
In the wave-equation migration velocity analysis (WEMVA) method, ${\bf K}$ is the differential residual migration operator \cite[]{sava2004}, 
it flattens the ADCIGs and also focuses the migrated image in the CMP direction; In the Differential Semblance Optimization (DSO) method \cite[]{shen2004},
\begin{eqnarray}
{\bf K} = {\bf 1} - {\bf R},
\end{eqnarray}
where ${\bf R}$ 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 1}-{\bf K}){\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}}\right)^{\prime}({\bf 1}-{\bf K})^{\prime}({\bf 1}-{\bf K}){\bf I}\right),
\end{eqnarray}
where $^{\prime}$ means the adjoint, $\frac{\partial{\bf I}}{\partial{\bf s}}$ is a linear operator which defines the 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 (Claerbout,???):
\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) = 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 $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 $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} \overline{D({\bf x}-{\bf h},{\bf x}_s,\omega)} U({\bf x}+{\bf h},{\bf x}_s,\omega), \label{eq_ic_bI}
\end{eqnarray}
where ${\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( \overline{\Delta D({\bf x}-{\bf h},{\bf x}_s,\omega)} {\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega) + \right. \nonumber \\
& & \left. \overline{{\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)} \Delta U ({\bf x}+{\bf h},{\bf x}_s,\omega) \right), \label{eq_ic_dI}
\end{eqnarray}
where ${\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)$ and ${\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega)$ are the background source and 
receiver wavefields computed with the background slowness ${\widehat s}({\bf x})$; 
$\Delta D({\bf x}-{\bf h},{\bf x}_s,\omega)$ and $\Delta U({\bf x}+{\bf h},{\bf x}_s,\omega)$ 
are the perturbed source wavefield and perturbed receiver wavefield, which are the results of the slowness perturbation $\Delta s({\bf x})$.
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}){\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}){\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}
%--------------------------------------------------------------------
However, the cost for computing the gradient of the objective function $J$ in the shot-profile domain 
is twice the cost of a shot-profile migration, to minimize $J$ many gradient and function evaluations have to be computed,
thus wave-equation tomography in the conventional shot-profile domain seems infeasible 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., 

\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 {biondi}, the prestack exploding reflector modeling aims to generate a smaller dataset to speed up iterations of migration velocity analysis. The method synthesizes areal source and receiver data at the surface by recursively upward propagating source and receiver wavefields, having prestack images obtained by wave-equation migration with an inaccurate velocity, ${\hat s}({\bf x})$, as initial conditions.
\par
The reduction of the dataset size to be migrated is achieved by modeling several wavefields at the same time and collecting them into super-areal data. To avoid crosstalk, due to correlation of unrelated wavefields during migration, phase encoding is applied in the modeling step. \cite{guerra} showed that random phases can be used to encode the modeling experiments such that the crosstalk is dispersed throughout the image. By performing more modeling realizations, the signal to noise ratio is improved. 
\par
For the areal shot migration of the synthesized data, both synthesized 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) = {\breve D}(x,y,z=0,{\bf p}_m,\omega) \end{array} \right.,
\label{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) = {\breve U}(x,y,z=0,{\bf p}_m,\omega) \end{array} \right.,
\label{permU}
\end{eqnarray}
where ${\breve D}$ and ${\breve U}$ defines, respectively, modeling realizations of the synthesized source and receiver wavefields at the surface. They are used as boundary conditions for the downward continuation. ${\widehat D}$ and ${\widehat U}$ are the downward propagated background source and receiver wavefields, respectively. The vector ${\bf p}_m$ represents the modeling realizations.
\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{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{permDIm}
\end{eqnarray}
where the wavefields indicated with a ~ $\widehat.$ ~ are the background wavefields and the perturbed source and receiver wavefields are $\Delta D$ and $\Delta U$, respectively. The perturbed wavefields satisfy the one-way wave equation:
\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{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 D}(x,y,z=0,{\bf p}_m,\omega) = 0 \end{array} \right.,
\label{permD}
\end{eqnarray}

\begin{eqnarray}
D({\bf x},\omega) &=& \left[D({\bf x,x_m},\omega) + I_D({\bf x_m,h}) \right ] e^{-ik_z\Delta z} \\
U({\bf x},\omega) &=& \left[U({\bf x,x_m},\omega) + I_U({\bf x_m,h}) \right ] e^{ik_z\Delta z}.
%{\tilde D}({\bf x},{\bf p}_m,\omega) &=& \left[D(x,z+\Delta z,\omega) + I_d(x-h,z+\Delta z) \right ] e^{-ik_z\Delta z} \\
%{\tilde U}({\bf x},{\bf p}_m,\omega) &=& \left[U(x,z+\Delta z,\omega) + I_u(x+h,z+\Delta z) \right ] e^{ik_z\Delta z} \\
\end{eqnarray}
Here The different initial conditions, $I_D$ and $I_U$, for the source and receiver wavefields respectively, are computed by a proper rotation of the original prestack image necessary for the case of using an inaccurate velocity and in the presence of dipping layers \cite[]{biondo}.
\par
Reduction of the dataset size to be migrated is achieved by modeling several wavefields at the same time and collecting them into a super-areal data. The wavefields resulting from the modeling realizations are computed by
\begin{eqnarray}
%{\tilde D}({\bf x},{\bf p}_m,\omega) &=& \left[D(x,z+\Delta z,\omega) + I_d(x-h,z+\Delta z) \right ] e^{-ik_z\Delta z} \\
%{\tilde U}({\bf x},{\bf p}_m,\omega) &=& \left[U(x,z+\Delta z,\omega) + I_u(x+h,z+\Delta z) \right ] e^{ik_z\Delta z} \\
{\tilde D}({\bf x},{\bf p}_m,\omega) &=& \sum_{{\bf x}_m} \left[D(x,z+\Delta z,\omega) + \alpha({\bf x}_m,{\bf p}_m,\omega) I_d(x-h,z+\Delta z) \right ] e^{-ik_z\Delta z} \\
{\tilde U}({\bf x},{\bf p}_m,\omega) &=& \sum_{{\bf x}_m} \left[U(x,z+\Delta z,\omega) + \alpha({\bf x}_m,{\bf p}_m,\omega) I_u(x+h,z+\Delta z) \right ] e^{ik_z\Delta z}
\end{eqnarray}
where $\alpha({\bf x}_m,{\bf p}_m,\omega)=e^{i\gamma({\bf x}_m,{\bf p}_m,\omega)}$ is the phase encoding function, ${\bf x}_m$ denotes model space coordinates and ${\bf p}_m$ represents the realizations. $\gamma$ is any function that minimizes the crosstalk. In the present case, we use a pseudo-random function following \cite [] {guerra}.
\par
The drawback of the data reduction is crosstalk generation when applying the imaging condition during migration. To attenuate crosstalk, source and receiver wavefields are encoded during the modeling step. \cite {guerra} showed that reasonable attenuation can be achieved by applying random phase encoding in $(x,z,\omega)$-dimensions. Furthermore, by migrating more modeling realizations, the signal to noise ratio is improved. 
\par

\section{Numerical examples}

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

