The image-space phase encoded gathers are obtained by the prestack exploding reflector modeling \cite[]{biondi2006,biondi2007}, 
which 
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$.

