\chapter{Image-space phase-encoded wavefields}
\label{chap:chap03}
\def\w{\omega} 
\def\bX{\boldsymbol{x}}
\def\bxi{\boldsymbol{\xi}}

%\begin{abstract}
This chapter introduces the image-space phase-encoded wavefields (ISPEW). ISPEW are computed using PERM along with phase-encoding techniques to further improve data reduction achieved with PERM. The phase encoding is performed during the modeling, in which every SODCIG in the initial conditions has its particular code. An important feature of ISPEW is that its modeling can be confined to a region of the subsurface where the velocity model is inaccurate, allowing migration velocity analysis using wavefield extrapolation to be easily solved in a target-oriented manner. The combination of these two characteristics dramatically decreases the cost of migration velocity analysis, especially in 3D projects.
%\end{abstract}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
In Chapter \ref{chap:chap02}, I introduced the pre-stack exploding-reflector model (PERM). Using the concept of exploding reflectors, PERM synthesizes new wavefields which enable computing pre-stack images by wavefield extrapolation methods at a low cost. This is achieved by carefully combining the modeling experiments, such that unrelated wavefields do not correlate during migration. In 3D, PERM data size can be up to two orders of magnitude smaller than plane-wave data. I will show in this chapter that phase-encoding techniques can be applied to the modeling of the wavefields, yielding further data size reduction while attenuating the deleterious effects of crosstalk.
\par
Phase encoding is a well established technique in radar \cite[]{levanon2004radar}, medical imaging \cite[]{bernstein:1452, weishaupt}, cryptography \cite[]{ruikang,javidi:2506} and wireless communication \cite[]{castoldi2002multiuser,zigangirov2004theory} communities. Phase encoding enables faster data acquisition, larger bandwidth and more reliable signal recovery. In wireless communication, for instance, systems using Code Division Multiple Access (CDMA), a method very popular to access the communication channels, allow several users to share the same communication channel without unwanted interferences (crosstalk). This is possible by encoding the information during transmission using sequences with unique correlation properties \cite[]{gold,dinan}. Because of the correlation properties, the information can be recovered with minimal distortion after decoding with the same sequence.
\par
The same idea has been long used in exploration seismic community to enable simultaneous shooting for acquisition with seismic vibrators \cite[]{ward:938,martin:943,bagaini:70}. In seismic processing, based on the linearity of wavefield propagation, phase-encoded wavefields are frequently used in the form of plane-wave phase encoding \cite[]{Schultz,whitmore,liu01,duquet2006} and random phase encoding \cite[]{romero,Sun01} to decrease the imaging cost. Only recently, phase-encoded wavefields have been applied to velocity estimation by waveform inversion \cite[]{vigh:VE135,ben-hadj-ali:2288,krebs:WCC177} and migration-velocity analysis using wavefield extrapolation \cite[]{shen:VE49,guerra:3964}. Phase-encoded wavefields can also be used to decrease the cost of computing the Hessian operator in least-squares migration \cite[]{tang:WCA95}.
\par
Usually, phase encoding is applied in the data space. Shot gathers are weighted with appropriate phase functions, ideally orthogonal, and combined into areal shots. The corresponding point sources are also encoded with the same phase functions and combined into areal source functions. During migration, the cross-correlation of related wavefields, encoded with the same phase function, originates strong amplitudes corresponding to the real reflectors. In the other hand, the cross-correlation of unrelated wavefields, encoded with different phase functions, originates attenuated amplitudes corresponding to the attenuated crosstalk.
\par
A similar strategy can be applied in the model space to encode PERM experiments \cite[]{guerra2008a,guerra2008b}. To each SODCIG used in the modeling of a pair of PERM wavefields is assigned a different pseudo-random sequence, such that, during migration, correlation of unrelated wavefields is attenuated. This characterizes the image-space phase-encoded wavefields (ISPEW), and it allows decreasing the sampling interval of the combined SODCIGs and using more than one reflector in the initial conditions.
\par
To provide relevant information for velocity updates, crosstalk must be incoherent. This has been observed by \cite{krebs:WCC177} in the context of waveform inversion using phase-encoded wavefields, for which velocity model updates improve if different codes are used along iterations. Similarly, in the context of migration velocity analysis by wavefield extrapolation using ISPEW, if necessary, several random realizations of ISPEW can be computed to be used in different velocity iterations.
\par
In migration velocity analysis using ray methods, it is common practice to limit the velocity update to a certain portion of the model space. This is especially true for deeper levels when the velocity model is accurate enough at shallower depths. However, when performing pre-stack depth migration during velocity iterations, the entire model space have to be imaged. To avoid the worthless computation of the migrated image in regions where the velocity is sufficiently accurate, the wavefields can be downward extrapolated up to a datum at the bottom of this region \cite[]{berryhill:1329, bevc:577, wang:E1}. Therefore, if the wavefields are datumized migration can be restricted to the region where the velocity accuracy needs to be improved.
\par
Since PERM wavefields are upward propagated, they can be collected at any depth level. In a situation similar to that described in the preceding paragraph, after being upward propagated with the inaccurate velocity of deeper levels, the wavefields could be collected at a the top of the inaccurate velocity region. Therefore, PERM naturally datumizes the wavefields. ISPEW takes advantage of this characteristic to decrease the cost of both the modeling and migration velocity analysis iterations.
\par
In this chapter, I describe how crosstalk is generated when migrating PERM data computed without respecting the minimum decorrelation distance and with more than one reflector in the initial conditions. Then, I introduce ISPEW and illustrate its effectiveness in reducing crosstalk with 2D and 3D examples.
 %
%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.

%% To decrease the prejudicial effect crosstalk has on the image quality, different strategies have been formulated to combine shot records based on the linearity of wavefield extrapolation. After the introduction of synthesis of plane waves \cite{Schultz}, several authors used this technique to perform imaging and velocity analysis using asymptotic methods \cite[]{akbar:1068} and wavefield extrapolation methods \cite[]{whitmore,mosher1998,zhang:E21}. \cite{liu01} provide a general framework to evaluate plane-wave decomposition in prestack source plane-wave migration and show the equivalence of shot profile and plane-wave migrations. Further generalizations of Schultz and Claerbout's method use wavefields of arbitrary shape in the subsurface to be upward propagated to the surface, where they define synthesis operators to be convolved with the shot profiles to generate areal shots \cite[]{Rietveld01,Rietveld02}. Migration of randomly phase-encoded shot-records \cite[]{romero} is a technique that enables migrating combined shots with minimal  cross-talk. Recently, random-phase encoding has been used to decrease the computational cost of waveform inversion \cite[]{ben-hadj-ali:2288}. The combination of phase-encoded and plane-wave areal shots can be used to drastically decrease the computational effort \cite[]{Sun01}. 
%\par

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Crosstalk generation}
In the combination of modeling experiments of Chapter \ref{chap:chap02}, two basic restrictions were applied to the initial conditions. First,  the initial conditions contained only one reflector. Second, the sampling period was sufficiently large to avoid crosstalk during migration (Figures \ref{fig:comb01b} and \ref{fig:cam04}). Considering that we probably need more than one reflector to define the velocity structure, the first restriction would increase the number of the modeling experiments by a factor equal to the number of selected reflectors. Moreover, because the choice of the sampling period depends on how many subsurface offsets are needed for performing velocity updates data size reduction is partly conditioned by the velocity inaccuracy, which causes energy to spread through subsurface offsets greater than zero. As we will see in this section, if those restrictions are not observed, two different kinds of crosstalk are formed: one related to the correlation of wavefields from different reflectors, other related to the correlation of wavefields from unrelated SODCIGs.
\par
The areal-shot migration of PERM data using wavefields computed with combined SODCIGs (equations \ref{eq:comb_bD and comb_bU}) produces a pre-stack image by applying the multi-offset imaging condition (equation \ref{eq:ic_bI}) (repeated here for the PERM case)
\begin{eqnarray}
\widehat{I_{\scriptscriptstyle P}}({\bf x},{\bf h}) = \sum_{\omega} \sum_{{\bf x}_s} \widehat{D_{\scriptscriptstyle P}}^{*}({\bf x}-{\bf h},{\bf x}_s,\omega) \widehat{U_{\scriptscriptstyle P}} ({\bf x}+{\bf h},{\bf x}_s,\omega).
\label{eq:ic_bIp}
\end{eqnarray}
\par
For conciseness, I dropped the dependency of the wavefields on the sampling period $\Delta {\bf x}=(\Delta x, \Delta y)$.
\par
Let us analyze the image formation for the case of crosstalk generated by correlating wavefields from different SODCIGs. In this case, the sampling period is shorter than the decorrelation distances $2(max\left | h_x \right |,max\left | h_y \right |)$ equations \ref{eq:comb_bD0} and \ref{eq:comb_bU0}
where $x_s$ are the pairs of PERM data.


To understand how crosstalk is formed using PERM wavefields, I use the same two reflectors model of Chapter \ref{chap:chap02}. Let us first consider the crosstalk formed by cross-correlating wavefields I modeled PERM data starting with the rotated images of Figure \ref{fig:dip05} and SODCIGs combined into sets with sampling period of 41 and 81 SODCIGs. As can be seen in Figures \ref{fig:comb01a}a and \ref{fig:comb01a}b, crosstalk occurs in the SODCIG according to a period which is half of the sampling period. Compare with Figure \ref{fig:comb01b} in which no crosstalk occurs. 
\plot{comb01a}{width=0.85\textwidth}{Areal-shot migration of PERM data synthesized from sets of SODCIGs selected with sampling period of: a) 41 and, b) 81.}

\par
We would like to have the ability of combining more modeling experiments without corrupting the velocity information with crosstalk. 

\par
To further decrease the number of modeling experiments, we can combine the reflectors in the initial condition. By doing so, two events are present both in the source wavefield and in the receiver wavefield, as in Figure \ref{fig:dip06}. During migration, in addition to the expected cross-correlations, cross-correlation of wavefields from different reflectors will occur, generating crosstalk between reflectors. Areal-shot migration of PERM data synthesized from a set of SODCIGs with sampling period of 163 is shown in Figure \ref{fig:comb02}. Both reflectors are simultaneously injected in the initial condition. Notice the presence of reflector crosstalk labeled with `C'.
\plot{comb02}{width=0.85\textwidth}{Areal-shot migration of PERM data synthesized from a set of SODCIGs selected with sampling period of 163. The two reflectors are simultaneously injected to the model. Notice the reflector crosstalk resulting from the cross-correlation of the wavefields from the horizontal reflector with  that from the dipping reflector.}
\par
To illustrate the reflector crosstalk generation, snapshots of the wavefield propagation are shown in (Figure \ref{fig:comb03}). Both wavefields are represented in the same panel and indicated by `D' (Downgoing) for the source wavefield and `U' (Upgoing) for the receiver wavefield. The panels on the left show the wavefields in the time-depth domain selected at horizontal positions where reflector crosstalk occurs in Figure \ref{fig:comb02}. The panels on the right are taken at the propagation time when the wavefields cross on the left panels. The crossing times vary from 0.14 s to 1.24 s, for the SODCIGs from right to left, respectively. The crossing times are a function of the distance between reflectors and the propagation velocity. In Figure \ref{fig:comb03}a, notice that for propagation times less than 0.14 s minus the period of the wavelet in time, no crosstalk will occur.
%For the simple case of constant velocity, the time at which reflector crosstalk occurs is  
\plot{comb03}{width=0.80\textwidth}{Snapshots of the wavefield propagation of Figure \ref{fig:comb02}. The panels on the left are selected at a horizontal position where the crosstalk occurs in Figure \ref{fig:comb02}. The panels on the right are taken at the propagation time when the wavefields cross on the left panel.}
\par
Crosstalk can be avoided if we take advantage of the imaging principle, which states that the reflector image is formed at the zero time of wavefield propagation. This is especially true when using migration velocity is accurate. When the migration velocity is inaccurate, the focusing of the image departs from the zero time, and this can potentially be used to update the velocity model \cite[]{sava:S209,yang:3715}. As Figure \ref{fig:comb03} shows, crosstalk is formed at times different from zero. Therefore, if the wavefields are cross-correlated within a time window centered at zero, reflector crosstalk might not occur depending on the distance in depth between the reflectors \cite[]{biondi07}. The time-windowed imaging condition for a single areal shot reads
\begin{eqnarray}
I{\scriptscriptstyle {P_w}}({\bf x},{\bf h}) = \sum_{\frac{-t_{\scriptscriptstyle w}}{2} \leq t \leq \frac{t_{\scriptscriptstyle w}}{2}} {\mathcal F}^{-1} \left [ D^{*}({\bf x}-{\bf h}) \right ] {\mathcal F}^{-1}\left [ U({\bf x}+{\bf h}) \right ],
\label{eq:wic_bI}
\end{eqnarray}
where $t_{\scriptscriptstyle w}$ is the length of the time window. When using one-way propagators, the wavefields are inverse Fourier transformed to time by ${\mathcal F}^{-1}$. Migration using the time-windowed imaging condition of  equation \ref{eq:wic_bI}  is shown in Figure \ref{fig:comb04}. The length of the time window is 0.06 s, which avoids the cross-correlation of events from different reflectors. The image is completely free of reflector crosstalk. 
\plot{comb04}{width=0.8\textwidth}{Areal-shot migration using the time-windowed imaging condition (equation \ref{eq:wic_bI}). The reflector crosstalk is completely eliminated.}
\par


In the example using the time-windowed imaging condition, the number of independent experiments is one fifth as many as in the original dataset. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Phase encoding the modeling experiments}
The image-space encoded gathers are obtained using PERM. 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.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\par

\cleardoublepage
%\bibliographystyle{seg}
%\bibliography{SEP,SEG,MISC,GEOTLE,allrefs}

