\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 image-space phase-encoded wavefields (ISPEWs). ISPEWs are computed using PERM along with phase-encoding techniques to further improve data reduction achieved with PERM. In this sense, ISPEWs can also be defined as image-space generalized source wavefields. The phase encoding is performed during the modeling, in which the source function of every event is coded using a particular coding sequence. Phase encoding during the modeling allows injection of more closely spaced SODCIGs than the technique showed in Chapter \ref{chap:chap02} and multiple reflectors, while diminishing the prejudicial effect of crosstalk during imaging. Modeling of ISPEWs 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 use of phase encoding combined with the target-oriented modeling  dramatically decreases the cost of migration velocity analysis, especially in 3D projects.
%\end{abstract}
%
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
Chapter \ref{chap:chap02} introduced the pre-stack exploding-reflector model (PERM). Using the concept of exploding reflectors, PERM synthesizes new wavefields that enables us to compute pre-stack images by wavefield extrapolation methods at a low cost. In 3D, PERM data size can be up to two orders of magnitude smaller than plane-wave data.
\par
The strategy used in Chapter \ref{chap:chap02} to ensure that crosstalk of wavefields does not occur during migration is, using of the linearity of wavefield propagation, to simultaneously inject SODCIGs separated by a decorrelation distance equal to twice the subsurface-offset range. This decorrelation distance in PERM modeling plays a role similar  to the grating interval in an optical fiber, which influences optical crosstalk in optical transmission systems and networks \cite[]{yamada}. In this optical fiber communication, crosstalk is further reduced (and bandwidth increased) by phase encoding the optical information \cite[]{teh}.
\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}. 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 for accessing communication channels, allow several users to share the same communication channel without crosstalk. This is achieved by encoding the information during transmission using sequences with unique correlation properties \cite[]{gold,dinan}. Information is recovered with minimal distortion after decoding with the corresponding sequence.
\par
Phase encoding has been long used in seismic exploration to enable simultaneous shooting for acquisition with seismic vibrators \cite[]{ward:938,martin:943,bagaini:70}. To decrease seismic imaging cost, using the concept of generalized sources, wavefields are usually phase-encoded using phase functions like the plane-wave phase function \cite[]{Schultz,whitmore,liu01,duquet2006} and random phase functions \cite[]{romero,Sun01}. Recently, phase-encoded wavefields have also 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, ideally, orthogonal phase functions, 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, yields strong amplitudes corresponding to the real reflectors. The cross-correlation of unrelated wavefields, encoded with different phase functions, yields 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}. A different pseudo-random sequence is assigned to each SODCIG used in the modeling of a pair of PERM wavefields, so that during migration correlation of unrelated wavefields is attenuated. This characterizes the image-space phase-encoded wavefields (ISPEWs), 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 when using generalized sources, crosstalk must be incoherent. This has been observed by \cite{krebs:WCC177} in the context of waveform inversion using phase-encoded wavefields, for which using different codes for different iterations yields more accurate models and improves convergence. Similarly, in the context of migration velocity analysis by wavefield extrapolation using ISPEWs, if necessary, several random realizations of ISPEWs 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, especially 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 has 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}. Hence, 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. Therefore, limiting the velocity update to the inaccurate velocity region is easily achieved using these wavefields. After being upward propagated with the inaccurate velocity of deeper levels, PERM wavefields are collected at the top of the inaccurate velocity region. Thus, 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 from combined SODCIGs separated by an interval that is shorter than the decorrelation distance and with more than one reflector in the initial conditions. I will show that phase-encoding techniques can be applied to the modeling of ISPEWs, yielding further data size reduction while attenuating the deleterious effects of crosstalk. I will use an example with the Marmousi model to illustrate the usefulness of ISPEWs to migration-velocity analysis.
 %
%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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{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 along midpoint 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 to subsurface offsets different from zero. As we will see in this section, if these restrictions are not observed, two different kinds of crosstalk are generated: one related to the correlation of wavefields from different reflectors, and the other related to the correlation of wavefields from unrelated SODCIGs. The first originates during the cross-correlation in time, whereas the second originates during the cross-correlation in space.
\par
Let us first consider the crosstalk from different reflectors. Since reflectors are simultaneously injected in the initial condition, several events are present both in the source wavefield and in the receiver wavefield, as in Figure \ref{fig:dip06} where we can see two events coming from the dipping and the deeper flat reflector. During migration, in addition to the expected cross-correlations between the source and receiver wavefields synthesized from the same reflector, cross-correlation of wavefields from different reflectors will occur, generating crosstalk between reflectors. Areal-shot migration of a pair of PERM wavefields synthesized from a set of SODCIGs with sampling period of 163 is shown in Figure \ref{fig:comb02}. The intermediate events between the dipping and the deeper reflectors are due to the cross-correlation of wavefields initiated at different reflectors. The panel on the left is the SODCIG taken at x = 0 m and the panel on the right is the zero-subsurface offset section. 
\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, labeled with `C', 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. 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 local dips, distance and propagation velocity between reflectors $\overline {v}$. For the simple case of parallel reflectors, the wavefield propagation time at which reflector crosstalk $t_{rc}$ occurs is 
\begin{eqnarray}
t_{rc} = \frac{0.5\Delta z_r}{\overline {v}},
\label{eq:t_rc}
\end{eqnarray}
where $\Delta z_r$ is the distance between reflectors.
\par
\plot{comb03}{width=0.80\textwidth}{Snapshots of propagation of wavefields used to compute the image of Figure \ref{fig:comb02}. Wavefields are labeled `D' (Downgoing) for the source wavefield and `U' (Upgoing) for the receiver wavefield. The panels on the left are selected at horizontal positions 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.}
Now let us consider the crosstalk formed by cross-correlating wavefields from unrelated SODCIGs of the same areal shot. From equation \ref{eq:I_perm2dx1}, repeated here,
\begin{eqnarray}
\begin{array}{lll}
\widehat{I_{\scriptscriptstyle P}}({\bf x},{\bf h};\Delta {\bf x}) & = & I_{\scriptscriptstyle P}({\bf x},{\bf h};\Delta {\bf x}) + \displaystyle{\sum_{\boldsymbol\xi'}\sum_{\bf x'}\sum_{\bf h'}}\displaystyle{\sum_{\boldsymbol\xi''}\sum_{\bf x''}\sum_{\bf h''}\sum_{n \ne m}}G_0(\boldsymbol\xi',{\bf x'-h'}) G^*_1({\boldsymbol\xi'},{\bf x-h}) \nonumber \\
&\times& G^*_1({\boldsymbol\xi''},{\bf x+h}) G_0(\boldsymbol\xi'',{\bf x''+h''}) \delta({\widehat{\bf x'}}-m\Delta{\bf x}) \delta({\widehat{\bf x''}}-n\Delta{\bf x}) \\
&\times& I(\widehat{\bf x'}{\bf ,h'})I(\widehat{\bf x''}{\bf ,h''}).
\end{array}
\label{eq:I_perm2dxa}
\end{eqnarray} 
we see that the crosstalk has structure similar to that of the desired image. In equation \ref{eq:I_perm2dxa} $G$ are the Green's functions from a point in the subsurface $\bf x$ to a point $\boldsymbol\xi$ at the datum where the wavefields are collected. The Green's functions can be propagated with different velocities, which is indicated by the subscripts. The sampling functions $\delta({\bf x}-n\Delta{\bf x})$ select the image $I$ at periods of $\Delta {\bf x}$. Crosstalk is not formed for $\Delta {\bf x} > \boldsymbol \eta$, where $\boldsymbol \eta$ is the decorrelation distance equal to twice the subsurface-offset range. If this criterion is not observed, crosstalk will occur in the SODCIGs according to a period of $\Delta {\bf x}$.
\par
To see how crosstalk from unrelated SODCIGs is formed, let us use the same two-reflectors model as in Chapter \ref{chap:chap02}. PERM data were modeled starting from the rotated images of Figure \ref{fig:dip05} using SODCIGs combined into sets with sampling period of 41 and 81 SODCIGs. Equation \ref{eq:I_perm1} shows that no crosstalk is generated if the sampling period is chosen to be the decorrelation distance of twice the subsurface-offset range, which in the present case must be greater than the distance spanned by 162 SODCIGs. Recall that the number of subsurface offsets in the original image is 81. As can be seen in Figures \ref{fig:ispew01}a and \ref{fig:ispew01}b, crosstalk occurs according periods of one-fourth and one-half of the sampling period, respectively. The corresponding ADCIGS at $x = 0$ m and the ADCIG computed from the image with no crosstalk are shown in the top panels of Figure \ref{fig:ispew02}. In the bottom panels we can see the corresponding $\rho$ scans, computed using equation D-7 in \cite{biondi:1283}. Notice that manual picking can identify the correct $\rho = 0.9$ in Figures \ref{fig:ispew02}a-b. Therefore, ray-based methods for velocity update can back-project the correct moveout information. However, when wavefield-extrapolation methods are used for velocity update, perturbed images computed from Figures \ref{fig:ispew01}a-b or \ref{fig:ispew02}a-b will potentially provide incorrect gradients.
\plot{ispew01}{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.}
\plot{ispew02}{width=0.85\textwidth}{ADCIGs (top) and $\rho$-panels (bottom) corresponding to images computed by wavefields modeled with sampling period of: a) 41, b) 81, and c) 163. Velocity information has been destroyed by the crosstalk in a) and b).} 




%, we consider a plane reflector such that the initial conditions are the same for modeling the source and the receiver wavefields. The modeling of the source and receiver PERM wavefields, ${\widehat D_{\scriptscriptstyle P}}$ and ${\widehat U_{\scriptscriptstyle P}}$, respectively, initiated simultaneously at SODCIGs separated by the 2D-sampling period $\Delta {\bf x}=(\Delta x,\Delta y)$ is represented by 


%\begin{eqnarray}
%{\widehat D_{\scriptscriptstyle P}}(\boldsymbol\xi) = \sum_n{\mathcal L}_0(\boldsymbol\xi,{\bf x-h-}n\Delta{\bf x})I({\bf x-h-}n\Delta{\bf x,h}),
%\label{eq:ispew_opD}
%\end{eqnarray}
%and
%\begin{eqnarray}
%{\widehat U_{\scriptscriptstyle P}}(\boldsymbol\xi) = \sum_n{\mathcal L}_0(\boldsymbol\xi,{\bf x+h-}n\Delta{\bf x})I({\bf x+h-}n\Delta{\bf x,h}),
%\label{eq:ispew_opU}
%\end{eqnarray}
%where, for conciseness, here $\bf x$ is the horizontal coordinates $(x,y)$, and the depth coordinate is not represented. In Chapter \ref{chap:chap02}, to avoid crosstalk the 2D-sampling period was chosen to be larger than the subsurface-offset range. Here, we use shorter sampling periods. 
%\par
%Following the same reasoning used to obtain equation \ref{eq:I_perm1}, the pre-stack image $\widehat{I_{\scriptscriptstyle P}}$ computed by the areal-shot migration of a pair of combined PERM wavefields is given by 
%\small
%\begin{eqnarray}
%\begin{array}{lll}
%\widehat{I_{\scriptscriptstyle P}}({\bf x},{\bf h};\Delta{\bf x}) &=& \displaystyle\sum_n {\mathcal L}'_0({\boldsymbol\xi},{\bf x-h-h'-}n\Delta{\bf x}){\mathcal L}_1(\boldsymbol\xi,{\bf x-h})I({\bf x-h-h'-}n\Delta{\bf x,h'}) \\
%&\times& \displaystyle\sum_m {\mathcal L}'_1({\boldsymbol\xi},{\bf x+h}){\mathcal L}_0(\boldsymbol\xi,{\bf x+h+h'-}m\Delta{\bf x})I({\bf x+h+h'-}m\Delta{\bf x,h'}).
%\end{array}
%\label{eq:ispew01}
%\end{eqnarray} 
%\normalsize
%\par
%For simplicity, let us assume that the wavefield propagations are performed with the same velocity. Therefore, the chain of the ${\mathcal L}$ operator and its adjoint approximates the identity operator, and equation \ref{eq:ispew01} can be recast as
%\small
%\begin{eqnarray}
%\widehat{I_{\scriptscriptstyle P}}({\bf x},{\bf h};\Delta{\bf x}) = I^2({\bf x,h};\Delta{\bf x}) + \displaystyle \sum_{n \neq m} I({\bf x-h-h'-}n\Delta{\bf x}) I({\bf x+h+h'-}m\Delta{\bf x,h'}),
%\label{eq:ispew02}
%\end{eqnarray} 
%\normalsize
%where the first term in the left-hand side is the PERM image that would have been computed by individually migrating wavefields generated at the selected SODCIGS, and the second term corresponds to the crosstalk.
%\par
%Equation \ref{eq:ispew02} can be further written as
%\small
%\begin{eqnarray}
%\begin{array}{lll}
%\widehat{I_{\scriptscriptstyle P}}({\bf x},{\bf h};\Delta{\bf x}) &= &I^2({\bf x,h};\Delta{\bf x}) + \displaystyle \sum_{n \neq m} \mathrm{\widetilde{rect}}(\frac{{\bf x}-n\Delta{\bf x}}{h_r}) \mathrm{\widetilde{rect}}(\frac{{\bf x}-m\Delta{\bf x}}{h_r}) I^2({\bf x,h};\Delta{\bf x}),
%\end{array}
%\label{eq:ispew03}
%\end{eqnarray} 
%\normalsize
%where  $h_r$ is the subsurface-offset range, and $\mathrm{\widetilde{rect}}$ is a scaled version of the $h_r$---wide 2D---rectangular function with period $\Delta {\bf x}$. From signal analysis, we know that the convolution of two periodic rectangular functions results in a periodic triangular function if for $\Delta {\bf x} \ge h_r$. In this case, the crosstalk term is zero.  
%\par
%From equation \ref{eq:ispew01} we can see that the image $I_{\scriptscriptstyle P}$ is formed by contributions from wavefields generated within a neighborhood around ${\bf x}$ of $(max\left |2h_x \right |,max\left |2h_y \right |)$ for symmetric SODCIGs with respect to subsurface offset or, more generally, twice the subsurface offset range. Therefore, to ensure that the areal-shot migrated image has kinematics at all subsurface offsets similar to those in the original isolated SODCIG, we need to model PERM data from a set of SODCIGs within that neighborhood around the central SODCIG.
%\par

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Crosstalk attenuation}
Let us now consider two different strategies to attenuate the two different types of crosstalk: the reflector crosstalk and the crosstalk from unrelated SODCIGs. In the first strategy, we will take advantage of the imaging principle to perform time cross-correlations within a window centered at zero time of wavefield propagation to avoid reflector crosstalk. In the second strategy, we will use random-phase encoding to combine the modeling experiments to attenuate both types of crosstalk. These are mutually excluding strategies, because when wavefields are phase encoded, the reference for the zero time of wavefield propagation is lost, since the frequency components of the wavefields are randomly injected in time into the modeling. We illustrate the crosstalk attenuation using a smoothed version of the Marmousi model.
\subsection{Time-windowed imaging condition}
From 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. This observation can be used to avoid crosstalk by applying a modified imaging condition. 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 time zero with a length that excludes the times at which crosstalk is formed, reflector crosstalk can be avoided \cite[]{biondi07}. The time-windowed imaging condition for a single pair of 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
The efficacy of the time-windowed imaging condition for attenuating the reflector crosstalk is independent of the number of modeled events; however, it does depend on the distance between reflectors. This is illustrated using a smoothed Marmousi model example (Figure \ref{fig:ismarm01}a). This model was input into  one-way Born modeling, along with a 2D reflectivity derived from the true Marmousi velocity model, to generate 360 split-spread shot gathers with 6000 m maximum offset. 
\par
Shot-profile migration of the Born data was performed with a migration velocity approximately $5\%$ slower than the modeling velocity for depth levels below the horizon shown in Figure \ref{fig:ismarm01}b. Above this horizon, migration velocity is equal to the modeling velocity. All the left panels shown in this example are SODCIGs selected at $x = 5000$ m. The right panel is the zero-subsurface-offset section. The background image contains 17 subsurface offsets 24 m apart, starting at -192 m. Four reflectors were selected to initiate the modeling of PERM wavefields. They were rotated according to the transformation described in Chapter \ref{chap:chap02}. The background image and the selected reflectors are shown in Figures \ref{fig:ismarm02}a-b, respectively, only for depths between 1500 and 3000 m. The rotated initial conditions for the modeling of source and receiver wavefields are shown in Figures \ref{fig:ismarm05}a-b, respectively.
\par
The horizontal distance between SODCIGs injected simultaneously into the modeling of one areal shot is 840 m, which means that 35 pairs of PERM wavefields are generated. This number corresponds to $\left . 1 \right / 10$ of the original shot gathers. The spatial sampling interval is greater than the decorrelation distance, which in this case is 792 m. Hence, we do not expect crosstalk from unrelated SODCIGs. 
\plot{ismarm01}{width=0.8\textwidth}{Velocity models for the Marmousi example: a) Smooth velocity model used to model the Born data. b) Background velocity model used to migrate the Born data, and to model and migrate PERM data.}
\plot{ismarm02}{width=0.8\textwidth}{a) Pre-stack image computed with the background velocity model. b) Selected reflectors from the background image to perform modeling of wavefields.}  
\plot{ismarm05}{width=0.8\textwidth}{Rotated initial conditions for modeling: a) source wavefields, and b) receiver wavefields.}  
\par
PERM wavefields were modeled and collected at a depth of 1500 m. A pair of PERM wavefields are shown in Figure \ref{fig:ismarm03}. These wavefields were migrated using the conventional (equation \ref{eq:I_perm}) and the time-windowed  (equation \ref{eq:wic_bI}) imaging conditions. The areal-shot migration using the conventional imaging condition is strongly contaminated with reflector crosstalk (Figure \ref{fig:ismarm04}a). The areal-shot migration of PERM wavefields using the time-windowed imaging condition with a time-window length of 0.016 s successfully attenuated the reflector crosstalk in regions where reflectors are sufficiently separated (Figure \ref{fig:ismarm04}b). However, as can be seen in the SODCIG in the left panel, some residual crosstalk from closely spaced reflectors still persists. This is because the time of wavefield propagation at which reflector crosstalk occurs is within the time window used for cross-correlation. This is easily seen by inserting the values for the local background velocity of 2540 m$\left/ \right .$s and the distance of 160 m between the two reflectors responsible for the reflector crosstalk into equation \ref{eq:t_rc}. The result, 0.125 s, is slightly smaller than the time-window length. Using a shorter time window could completely avoid crosstalk, but also could cause phase distortion.
\plot{ismarm03}{width=0.8\textwidth}{PERM wavefields for the Marmousi example: a) Receiver wavefield. b) Source wavefield.}  
\plot{ismarm04}{width=0.8\textwidth}{Pre-stack image computed with PERM wavefields and background velocity model using: a) the conventional imaging condition (equation \ref{eq:I_perm}), and b) the time-windowed imaging condition (equation \ref{eq:wic_bI}). Reflector crosstalk is avoided when reflectors are sufficiently separated. However, some residual crosstalk is still present (RC). Notice the phase difference of the PERM image due to the squaring of the wavelet when compared to the windowed reflectors of Figure \ref{fig:ismarm02}b.}
\par
The imaging principle states that the reflector image is formed at the zero time of wavefield propagation when the migration velocity is accurate. When the migration velocity is inaccurate, the focusing of the image departs from the zero time, and this can be potentially used to update the velocity model \cite[]{sava:S209,yang:3715}. Depending on the magnitude of the velocity errors and the distance between reflectors, the use of the time-windowed imaging condition can corrupt the velocity information.  
\par 
Next, we will see how phase-encoding the modeling experiments can attenuate crosstalk without the risk of affecting the velocity information.

\subsection{Phase encoding the modeling experiments}
Phase encoding is a well-established technique for decreasing the cost of seismic imaging by linearly combining the shot records, while maintaining the image quality \cite[]{Schultz,whitmore,romero,Sun01,liu01,duquet2006}. Wavefields are usually phase-encoded in the data space simply because conventional wavefields are initiated and recorded at the boundaries of the reflection seismic problem. Therefore, the phase functions are parametrized according to the source index or source coordinates, which are data space parameters.
\par
In a similar way, PERM wavefields initiated on the reflectors can also be phase-encoded. In this case, the phase functions are parametrized according to the model space coordinates and selected reflector, which are model space parameters. The parametrization of these phase functions characterizes the image-space phase encoding. The resulting wavefields are the image-space phase-encoded wavefields (ISPEWs). We will see that randomly phase-encoding the modeling experiments enables us to have more than one reflector in the initial conditions and to use a spatial sampling period  smaller than the decorrelation distance. This is not possible when using the time-windowed imaging condition strategy.
\par
To describe how ISPEWs are generated, let us slightly change the previous formulation of PERM to include the simultaneous modeling of more than one reflector. Horizontal reflectors are used for simplicity. For all the frequencies, using different random realizations ${\bf q}$ of the initial conditions $\widetilde I$ for modeling  ISPEWs can be described by
\begin{eqnarray}
\widetilde{D_{\scriptscriptstyle I}}(\boldsymbol\xi,{\bf q},\omega) = \sum_{\bf x}\sum_{\bf h} G_0(\boldsymbol\xi,{\bf x-h},\omega) {\widetilde I}(\widehat{\bf x}{\bf ,h},{\bf q},\omega),
\label{eq:ispew_opDdx}
\end{eqnarray}
and
\begin{eqnarray}
\widetilde{U_{\scriptscriptstyle I}}(\boldsymbol\xi,{\bf q},\omega) = \sum_{\bf x}\sum_{\bf h} G_0(\boldsymbol\xi,{\bf x+h},\omega) {\widetilde I}(\widehat{\bf x}{\bf ,h},{\bf q},\omega).
\label{eq:ispew_opUdx}
\end{eqnarray}
The subscript of the Green's function $G$ denotes propagation with a background velocity $s_0({\bf x})$. The initial conditions ${\widetilde I}$ are computed as
\begin{eqnarray}
{\widetilde I}(\widehat{\bf x}{\bf ,h},{\bf q},\omega) = \sum_m \sum_j \delta({\widehat{\bf x}-m\Delta{\bf x}}) \beta(\widehat{\bf x},j,{\bf q},\omega) W_j(\widehat{\bf x}{\bf ,h}) I(\widehat{\bf x}{\bf ,h}),
\label{eq:ispew_Itil}
\end{eqnarray}
where $W_j$ selects the reflector $j$ by identifying and windowing it in the pre-stack image $I$, and $\beta(\widehat{\bf x},j,{\bf q},\omega)$ is a pseudo-random phase-encoding function defined as
\begin{eqnarray}
\beta(\widehat{\bf x},j,{\bf q},\omega) = e^{i\epsilon(\widehat{\bf x},j,{\bf q},\omega)},
\label{eq:ispew_rnd}
\end{eqnarray}
with $\epsilon(\widehat{\bf x},j,{\bf q},\omega)$ usually being a uniformly distributed pseudo-random sequence with zero mean. \cite{guerra2008b} also use Gold codes \cite[]{gold}, sequences widely used in third- and fourth-generation cellphones, to phase-encode the modeling. The pseudo-random phase-encoding function causes the frequency components of the initial conditions to be randomly injected into the modeling. Source and receiver wavefields initiated at the same SODCIG and from the same reflector are equally encoded, whereas source and receiver wavefields initiated at different SODCIGs and from the different reflectors have different codes assigned to them.
\par
The recursive downward propagation with a different velocity $s_1({\bf x})$ is performed according to
\begin{eqnarray}
\widetilde{D_{\scriptscriptstyle I}}({\bf x},\omega,{\bf q}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x},\omega) \widetilde{D_{\scriptscriptstyle I}}({\boldsymbol\xi}\omega,{\bf q}),
\label{eq:ispew_opDzdx} 
\end{eqnarray}
and
\begin{eqnarray}
\widetilde{U_{\scriptscriptstyle I}}({\bf x},\omega,{\bf q}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x},\omega) \widetilde{U_{\scriptscriptstyle I}}({\boldsymbol\xi}\omega,{\bf q}).
\label{eq:ispew_opUzdx}
\end{eqnarray}
\par
The lateral shifts of the wavefields for the multi-offset imaging condition are represented by
\begin{eqnarray}
\begin{array}{lll}
\widetilde{D_{\scriptscriptstyle I}}({\bf x-h},\omega,{\bf q})&=&\displaystyle{\sum_{\boldsymbol\xi}\sum_{\bf x'}\sum_{\bf h'} \sum_m \sum_j G^*_1({\boldsymbol\xi},{\bf x-h},\omega)G_0(\boldsymbol\xi,{\bf x'-h'},\omega)} \\
&\times&\displaystyle{ \delta({\widehat{\bf x}-m\Delta{\bf x}}) \beta(\widehat{\bf x},j,{\bf q},\omega) W_j(\widehat{\bf x}{\bf ,h'}) I(\widehat{\bf x}{\bf ,h'})}, 
\end{array}
\label{eq:ispew_lshDdx} 
\end{eqnarray}
and
\begin{eqnarray}
\begin{array}{lll}
\widetilde{U_{\scriptscriptstyle I}}({\bf x+h},\omega,{\bf q})&=&\displaystyle{\sum_{\boldsymbol\xi}\sum_{\bf x'}\sum_{\bf h'} \sum_m \sum_j  G^*_1({\boldsymbol\xi},{\bf x+h},\omega) G_0(\boldsymbol\xi,{\bf x'+h'},\omega)} \\
&\times&\displaystyle{\delta({\widehat{\bf x}-m\Delta{\bf x}}) \beta(\widehat{\bf x},j,{\bf q},\omega) W_j(\widehat{\bf x}{\bf ,h'}) I(\widehat{\bf x}{\bf ,h'})}, 
\end{array}
\label{eq:ispew_lshUdx}
\end{eqnarray}
\par
Applying the cross-correlation imaging condition to the wavefields of equations \ref{eq:ispew_lshDdx} and \ref{eq:ispew_lshUdx} and summing over frequency and over realizations gives
\begin{eqnarray}
\begin{array}{lll}
\widetilde{I_{\scriptscriptstyle I}}({\bf x},{\bf h})&=\displaystyle{\sum_{\omega}\sum_{\boldsymbol\xi'}\sum_{\bf x'}\sum_{\bf h'} \sum_{\boldsymbol\xi''}\sum_{\bf x''}\sum_{\bf h''}\sum_{\bf q} \sum_{m} \sum_{n} \sum_{j} \sum_{l}}&\\
&\times G^*_0({\boldsymbol\xi'},{\bf x'-h'},\omega) G_1({\boldsymbol\xi'},{\bf x-h},\omega)G^*_1({\boldsymbol\xi''},{\bf x+h},\omega) G_0(\boldsymbol\xi'',{\bf x''+h''},\omega)&\\
&\times \delta({\widehat{\bf x}-m\Delta{\bf x}}) \delta({\widehat{\bf x}-n\Delta{\bf x}}) \beta(\widehat{\bf x},j,{\bf q},\omega) \beta(\widehat{\bf x},l,{\bf q},\omega) &\\
&\times W_j(\widehat{\bf x}{\bf ,h'}) W_j(\widehat{\bf x}{\bf ,h''}) I(\widehat{\bf x}{\bf ,h'})I(\widehat{\bf x}{\bf ,h''}), &\\
\end{array}
\label{eq:I_ispew2dx}
\end{eqnarray} 
which can be recast as
\begin{eqnarray}
\begin{array}{lll}
\widetilde{I_{\scriptscriptstyle I}}({\bf x},{\bf h})&=I_{\scriptscriptstyle {\Delta x},j}({\bf x},{\bf h}) &\\
&+ \displaystyle{\sum_{\omega}\sum_{\boldsymbol\xi'}\sum_{\bf x'}\sum_{\bf h'} \sum_{\boldsymbol\xi''}\sum_{\bf x''}\sum_{\bf h''}\sum_{\bf q} \sum_{m \ne n} \sum_{j \ne l}}&\\
&\times G^*_0({\boldsymbol\xi'},{\bf x'-h'},\omega) G_1({\boldsymbol\xi'},{\bf x-h},\omega)G^*_1({\boldsymbol\xi''},{\bf x+h},\omega) G_0(\boldsymbol\xi'',{\bf x''+h''},\omega)&\\
&\times \delta({\widehat{\bf x}-m\Delta{\bf x}}) \delta({\widehat{\bf x}-n\Delta{\bf x}}) \beta(\widehat{\bf x},j,{\bf q},\omega) \beta(\widehat{\bf x},l,{\bf q},\omega) &\\
&\times W_j(\widehat{\bf x}{\bf ,h'}) W_l(\widehat{\bf x}{\bf ,h''}) I(\widehat{\bf x}{\bf ,h'})I(\widehat{\bf x}{\bf ,h''}). &\\
\end{array}
\label{eq:I_ispew2dx1}
\end{eqnarray} 
The first term in the right-hand side of equation \ref{eq:I_ispew2dx1} is the desired image we would obtain by  modeling and migrating PERM wavefields initiated at isolated reflectors of isolated SODCIGs and is the result of the cross-correlation of encoded wavefields with $m=n$ and $j=l$. The second term represents the attenuated crosstalk, which includes the attenuated reflector crosstalk and the attenuated crosstalk from unrelated SODCIGs. Crosstalk attenuation is achieved by the cross-correlation of quasi-orthogonal sequences, whose values are small if compared with the zero lag of the auto-correlation. Ideally, the attenuated crosstalk is unstructured and occurs as speckled noise throughout the image.  
\par
Theoretically, according to the law of large numbers, if the number of random realizations is sufficiently large, the crosstalk term is negligible. When summing migration results of different random realizations, reflectors constructively interfere, whereas crosstalk destructively interferes. In practice, a satisfactory crosstalk attenuation is achieved using a small number of realizations. 
\par
I use the Marmousi data to illustrate the generation and migration of ISPEWs. The rotated SODCIGs of Figure \ref{fig:ismarm05} are phase-encoded to generate the initial conditions. We phase-encode the SODCIGs as well as the reflectors. The sampling period is 264 m, which roughly corresponds to $\left . 1 \right /3$ of the decorrelation  distance. We compute four random realizations, corresponding to 44 pairs of ISPEWs. Two pairs of ISPEWs from different random realizations initiated at the same SODCIGs are shown in Figure \ref{fig:ismarm06}.
\plot{ismarm06}{width=0.8\textwidth}{ISPEWs from different random realizations initiated at the same SODCIGs for the Marmousi example: a,c) Receiver wavefields. b,d) Source wavefields.}  
\par
The areal-shot migration of ISPEWs is shown in Figure \ref{fig:ismarm07}. At the top, we see the result of migrating 44 wavefields from the four random realizations. At the bottom, we see the result of migrating 11 wavefields from one single random realization. Although the two images have a comparable quality, more realizations give a cleaner image.
\par 
\plot{ismarm07}{width=0.8\textwidth}{Pre-stack  images computed with: a) Four random realizations of ISPEWs, and b) a single random realization.}
To analyze the quality of the kinematic information, we compute residual-moveout panels from images computed with the original data (Figure \ref{fig:ismarm08}a), with PERM wavefields and the time-windowed imaging condition (Figure \ref{fig:ismarm08}b), with all four random realizations (Figure \ref{fig:ismarm08}c), and with a single realization (Figure \ref{fig:ismarm08}d). Notice that in the ADCIGs corresponding to images computed with PERM wavefields and ISPEWs, less events are present since we selected some reflectors to model these wavefields. The corresponding reflectors in the original shot-profile image are highlighted by the green boxes in Figure \ref{fig:ismarm08}a.
\par
The residual moveout information in the four panels is very similar. However, some crosstalk not entirely rejected by the time-windowed imaging condition causes the residual-moveout information lose resolution in Figure \ref{fig:ismarm08}b when compared to Figures \ref{fig:ismarm08}c-d.  
\par
The cost for obtaining the images in Figure \ref{fig:ismarm08} widely varies. For instance, migrating one random realization of ISPEWs is approximately 30 times faster than migrating all the 360 original shots. Using the same wavefields, computing the gradient of the objective function of migration-velocity analysis by wavefield-extrapolation can be 60 times faster.  This difference in performance is even more dramatic when we consider that several iterations of migrations and gradient computations are performed during migration velocity  optimization by wavefield-extrapolation. In Chapter \ref{chap:chap04}, we will see the usefulness of ISPEWs in migration velocity  optimization by wavefield-extrapolation. Chapter \ref{chap:chap05} will use these wavefields to optimize the migration velocity for a 3D survey from the North Sea.
\plot{ismarm08}{width=0.8\textwidth}{ADCIGs (top) and $\rho$-panels (bottom) corresponding to images computed by: a) Shot-profile migration of 360 shot gathers, b) areal-shot migration of 35 PERM wavefields using the time-windowed imaging condition, c) areal-shot migration of 44 ISPEWs corresponding to four random realizations, and d) areal-shot migration of 11 ISPEWs corresponding to a single random realization. The moveout information is basically the same.} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%\section{3D-{ISPEW} from common-azimuth migrated images}
%In Chapter \ref{chap:chap02}, we saw that using initial conditions computed under the common-azimuth assumption enables us to continuously sample the SODCIGs in the cross-line direction. Consequently, the number of modeling experiments is determined by the sampling period of SODCIGs in the in-line direction. In the 3D example of Chapter \ref{chap:chap02} this sampling period was chosen to be 48 times the cross-line spacing to not generate crosstalk, which is greater than the decorrelation distance given that the number of subsurface offsets is 21. Here, we use the same pre-stack image to compute ISPEW according to a sampling period smaller than the decorrelation distance.
%\par
%For the same reasons we can continuously sample in the cross-line direction the SODCIGs from images computed with common-azimuth migration (CAM), the phase-encoding function used to compute ISPEW is only a function of the x--coordinate and the reflector index. No encoding is necessary in the cross-line direction, since no spatial cross-correlation is performed in this direction during the areal-shot migration of ISPEW computed from CAM images.
%\par
%To model 3D-ISPEW, the sampling period in the in-line direction was chosen to be 16 times the cross-line spacing. With this sampling period we model 16 pairs of ISPEW, which is $\left . 1 \right / 3$ the number of pairs of PERM wavefields used to generate Figure \ref{fig:cam05}. A 3D-receiver ISPEW is shown in Figure \ref{fig:ispew3d00}.
%\par 
%Areal-shot migration of 3D-ISPEW is shown in Figure \ref{fig:ispew3d01}. Compare it with Figure \ref{fig:cam05}. 
%\plot{cam03}{width=0.8\textwidth}{3D-receiver ISPEW. The left panel is the in-line at $y = 1200$ m, the right panel is the cross-line at $x = 1400$ m, and the top panel is the time-slice at $t = 0.5$ s.}
%The areal-shot migration of 3D-ISPEW is shown in Figure \ref{fig:ispew3d01}.   
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
\section{Conclusions}
In this chapter, PERM was used to model wavefields without observing the decorrelation distance between SODCIGs and using more than one reflector in the initial conditions. Doing this generates two types of crosstalk are generated: the reflector crosstalk and the crosstalk from unrelated SODCIGs. The reflector crosstalk is generated when performing the cross-correlation  of the wavefields in time, whereas the crosstalk from unrelated SODCIGs is generated when performing the cross-correlation of the wavefields in space. 
\par
Since reflectors are imaged at time zero of wavefield propagation, an effective strategy to almost completely attenuate the reflector crosstalk is to cross-correlate wavefields only within a small time window around the zero time of wavefield propagation. To apply the time-windowed imaging condition, the frequency slices of the wavefields must be stored, which can be impractical for 3D applications. This strategy does not avoid the crosstalk from unrelated SODCIGs and, therefore, the decorrelation distance has to be used. In addition, in the presence of velocity inaccuracy, this strategy can corrupt the velocity information. 
\par
A more general method for attenuating crosstalk, independent of its origin, is to phase-encode the modeling experiments. Since the phase-encoding sequences are defined as a function of image parameters, namely spatial coordinates and reflector index, the wavefields are called image-space phase-encoded wavefields -- ISPEWs. We exemplified the usefulness of ISPEWs with an example using the Marmousi data. The migrated images using these wavefields have kinematic information similar to that obtained with shot-profile migration for the selected reflectors, as shown in the residual-moveout panels, but at a much lower cost. In the next chapter, we will see that they also provide a similar gradient for the migration-velocity analysis by wavefield-extrapolation objective function. 

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

