\def\figdir{./Fig}
\lefthead{Guerra}
\righthead{Ph.D Proposal}
\nokeywords

\begin{center}
\hrule
\vspace{0.3in}
\title{RESEARCH PROPOSAL}
\hrule
\vspace{1.0in}
\title{Wave equation velocity inversion \\
using exploding reflectors modeling}
\vspace{0.5in}
\nokeywords
\email{guerrac@stanford.edu}
\author{Claudio Guerra}
\par
10:00am, Monday June 16, 2008
\par
\vspace{1.0in}
\includegraphics[width=0.15\textwidth]{Fig/stanfordseal.ps}
\end{center}

\par

\hyphenation{Cla-er-bout}

\newpage

\section{ABSTRACT}
 
\subsection{Problem} 
To obtain a seismic image of the Earth's subsurface, estimates of the seismic velocity are required at different steps of seismic processing. The accuracy of the results depends on how precisely the velocity estimates represent the subsurface, based on the physical assumptions of the wave propagation.
\par
Velocity model estimation involves some transformation of the data and interpretation of velocity parameters. For instance, the normal-moveout (NMO) velocity is obtained by fitting hyperbolas to data sorted to common-midpoint (CMP).  Interval velocity for depth migration is determined iteratively by: 1) transforming data into common-image gathers (CIG) by means of prestack-depth migration with a presumed velocity model; 2) measuring the misfocusing on CIG; and 3) updating the velocity model based on the misfocusing parameter of step 2.
\par
Ray methods are widely used in depth migration velocity analysis because of its flexibility and low cost compared to wave-equation methods. However, in areas where the high-frequency approximation is violated and the wavefields become too complex, the dependence of velocity-model definition on Kirchhoff migration might result in inaccurate velocity models. This may lead to unreliable results even if wave-equation migration is used to produce the final image. In this case, it is therefore desirable to use wave-equation methods during depth migration velocity analysis.

%Most of the velocity estimation methods rely on the kinematics of primary reflections. However, data do not comprise solely primary reflections. Multiple reflections are, for instance, events that can obscure primaries making it difficult velocity model definition. Additionally, data is acquired subjected to some areal aperture limitation which decreases the velocity resolution. Besides non-linear, these factors make velocity estimation an ill-posed problem. Therefore, using geological a priori knowledge and making  assumptions on the Earth model help to constrain the velocity model solutions. 
%\par
% It is an expensive process because every iteration involves depth migration and back-projection of the misfocusing property into the slowness (velocity) model.

\subsection{Importance}
Oil exploration is being conducted in geologically complex areas. In the last ten years, huge oil reserves have been continuously discovered below salt bodies in Gulf of Mexico. Recently, it was announced 5 billion barrels of light oil reserve in sub-salt offshore Santos Basin, Brazil. These sub-salt oil discoveries can be partially credited to depth imaging.
\par
As a reliable image in depth is strongly dependent on a accurate velocity model, velocity estimation is one of the most important problems in the exploration community. Velocity model estimation by wave-equation methods is a very active area of research. To be applied in a production environment, however, its cost is still a limiting factor. Therefore, one of the challenges exploration geophysics faces is to decrease the cost of wave-equation methods, specially wave-equation migration, for velocity analysis.

\subsection {Proposal}
Data reduction has a great impact on the final cost of migration. For wave-equation migration velocity analysis, one can take advantage of the linearity of the wave propagation to linearly combine shots using phase encoding functions \cite{romero,liu}. Decomposition and combination of data in plane-waves is capable of decreasing the amount of data by one order of magnitude in 2D problems \cite{liu}.
\par
Aiming at a significantly higher data reduction, \longcite{Biondi.sep.124.biondo1,biondo2} presented an extension of the exploding reflector model to generate prestack data suitable for wave-equation migration velocity analysis. This method uses the exploding reflector concept \cite{loewenthal} to synthesize areal sources and areal receivers along the entire survey at the surface. The initial conditions for the modeling are prestack images obtained with wave-equation migration.
\par
Here, I propose a hybrid approach to wave-equation migration velocity analysis under the framework of the prestack exploding reflector modeling with one-way extrapolator. By using one-way operators, phase encoding methods can be implemented to adequately reduce input data to migration. I will also investigate the application of modified imaging conditions to attenuate crosstalk due to data reduction. I show that by selecting a few reliable reflectors in conjunction with crosstalk attenuation strategies makes prestack exploding reflector modeling applicable to field data. A question to be addressed is how wave-equation velocity inversion performs under this horizon-based prestack exploding reflector model.
\par
This research will contribute to make wave-equation migration velocity estimation of practical use.
%\par
%\longcite{Biondi.sep.124.biondo1,biondo2} uses the two-way wave equation to propagate the wavefields. I propose to use the one-way wave equation to achieve additional computational saving. Furthermore, under the one-way framework, I propose a phase encoding-like scheme during the modeling to reduce cross-talk between events, resulting from the combination of individual modeling experiments, while applying the imaging condition during the areal-shot migration. I will also investigate the use of different imaging conditions in the areal shot migration to evaluate their impact on crosstalk attenuation. 
%\par
%The goals of this research are to extend the method, initially derived in 2D, to 3D and to define a methodology for velocity update under the framework of the prestack-exploding-reflector model. This should enable a complete cost effective procedure that performs velocity model estimation for 3D prestack depth migration.

\section {INTRODUCTION}
In depth-imaging projects, specially in geologically complex areas, the definition of an accurate velocity model is decisive for achieving a reliable final image. Much effort is expended on building a geologically reasonable velocity model. As velocity model building for depth-imaging projects is strongly dependent on Kirchhoff algorithms, however, in areas where ray theory fails, inaccurate updates of the velocity model and even divergence of the solution can occur. 
\par
Wave-equation methods have the potential to overcome the limitations ray-based methods face when the velocity model is complex with sharp and irregular boundaries. In wave-equation methods, velocity model updates are obtained by minimizing a measure of image misfocussing. The differential semblance optimization (DSO) \cite{symes:654,shen:2132,shen:2249} and wave-equation migration velocity analysis (WEMVA) are the most promising wave-equation methods. In particular, WEMVA \cite{savaphd,sava03,sava04} measures the misfocusing in a non-automated way by picking a velocity parameter that best flattens angle-Domain common image gathers (ADCIGs) after residual prestack depth migration or, alternatively, residual depth moveout. In the following I give a brief review of WEMVA.
\par

\subsection{Wave-equation migration velocity analysis}
\longcite{savaphd} and \longcite{sava03,sava04} introduced WEMVA which iteratively updates the velocity model using an image perturbation, $\Delta I$, as a function of the slowness parameter, $\rho=\frac{s}{s_0}$. This parameter is obtained by interactively picking the values that best flatten reflectors on residual migrated ADCIGs. $s$ and $s_0$ are the trial slowness and the reference slowness, respectively. The slowness perturbation, $\Delta S = s - s_0$, is determined by minimizing the objective function: 
\begin{equation}
J(\Delta S) =  \left | \right | \Delta I - {\bf L} \Delta S  \left| \right | ^2.
\label{eq:wemva1}
\end{equation}
In equation \ref{eq:wemva1}, the image perturbation is obtained by applying a linearized version of the residual prestack depth migration onto the reference image, $I_0$, using the interpreted slowness parameter, $\rho$. Doing so, one ensures that the Born approximation is not violated when the velocity changes are too large. The operator ${\bf L}$ is the chaining of an imaging operator and the Born scattering operator, ${\bf S}$, represented by:
\begin{equation}
{\bf S}(U_0,\Delta S) \approx -i \frac{dk_z}{ds} \Big{|}_{s=s_0} \Delta S \Delta z U_0.
\label{eq:wemva2}
\end{equation}
Note that the operator ${\bf S}$ depends on the reference wavefield ,$U_0$, and on the slowness perturbation, $\Delta S$. The term $\frac{dk_z}{ds} \Big{|}_{s=s_0}$ is given by:
\begin{equation}
\frac{dk_z}{ds} \Big{|}_{s=s_0} = \frac{\omega ^2 s_0}{\sqrt{\omega ^2 s_0^2 - \left | {\bf k} \right | ^2}},
\label{eq:wemva3}
\end{equation}
where $\bf k$ is the horizontal wavenumber and $\omega$ is the angular frequency.
\par
The inversion of the local information given by the image perturbation into global slowness updates is performed by an iterative algorithm, such as conjugate gradients. The image perturbation is back-propagated into slowness perturbation along wave-paths. 
\par
In 3D, WEMVA has considerable storage and computational costs \cite{savaphd}. Every iteration involves roughly the cost of two wave-equation migrations. The practical implementation of wave-equation migration velocity analysis, therefore, is still a challenge. There are opportunities for improvements like automatic picking of the velocity parameter \cite{biondo3} and reduction of the storage and computational cost by decreasing the amount of input data to migrate.
%Downward-continuation migration is performed in one of two to main ways: shot-profile and source-receiver migration. Provided that some requirements are fulfilled, these two schemes are equivalent \cite{biondi1}. The efficiency of one method compared to the other heavily depends on the acquisition geometry. Normally, 3D shot-profile migration is more efficient than 3D source-receiver migration when the number of shots is relatively small compared to the number of receiver positions, as in land and ocean-bottom-cable acquisition geometries. The use of shot-profile migration in 3D streamer marine surveys, however, can be cost prohibitive. 
%\par
%For the case of migration velocity analysis, the prestack image used as initial condition in the modeling is supposed to have been computed with an inaccurate velocity model. The modeled areal data are to be used in subsequent migrations through the velocity updating process. These areal data must be a smaller dataset than the one used in the initial migration. This is achieved by combining individual modeling experiments into super-areal data. Moreover, the areal data must maintain the kinematic information necessary to perform migration velocity analysis.

\subsection{Previous data reduction schemes}
After the seminal work of \longcite{Schultz}, in which they introduced the concept of  synthesis of plane waves, different strategies have been formulated to reduce the computational cost of wave-equation migration by combining shot records based on the linearity of wavefield propagation. The combination of shots drastically reduces the cost of wave-equation migration but generates crosstalk between uncorrelated shots during imaging. The crosstalk can be attenuated by using sufficient linear combinations of shots. In the case of plane-wave migration, if the set of synthesized plane-waves adequately samples the range of propagation angles present in the original data the crosstalk is largely attenuated \cite{liu}.
\par
\longcite{Rietveld01} used the idea of predefined wavefields of arbitrary shape at the subsurface from \longcite{Schultz} to synthesize areal sources at the surface. These areal sources convolved with shots synthesize areal shots to be migrated using an areal shot migration scheme. \longcite{Rietveld02} further developed the method to prestack migration using a target-oriented controlled illumination strategy. The phase encoding migration of \longcite{romero} linearly combine shots using phase functions different from the plane-wave. They use random phase encoding which disperses the crosstalk energy throughout the image and linear phase encoding which shifts the crosstalk outside the migration domain. They show that the migration results deteriorate as the number of encoded shots increases. 
\par
The combination of more than one phase encoding technique achieves additional reduction of the computational effort. \longcite{Sun01} and \longcite{zhang:1031} show good results on Marmousi data using mixed phase encoding and sign-opposite phase encoding on previously plane-wave encoded data. \longcite{liu} 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. Additionally, the use of plane-wave migration enables the use of different coordinate systems suited to image steep dips \cite{shan:969,shan:2190}.
\par
These previous works synthesize a smaller dataset by combining the recorded shot records. Under a different perspective, \longcite{Biondi.sep.124.biondo1,biondo2} introduces the prestack exploding reflector modeling. This method synthesizes source and receiver wavefields along the entire survey at the surface, in the form of areal data, starting from a prestack migrated image cube represented by subsurface-offset domain common-image gathers (SODCIGs). For the case of migration velocity analysis, the aim is to generate a small dataset by combining several individual modeling experiments, while maintaining the kinematics necessary to perform migration velocity analysis. As in the other data reduction schemes, crosstalk is also generated during imaging. I will present strategies to attenuate it.
%\par
%If the migration velocity is correct and assuming sufficiently good illumination, the reflectors focus around the zero subsurface offset and the method reduces to the conventional exploding reflector modeling. However, when the migration velocity is inaccurate reflectors do not focus close to the zero subsurface-offset anymore and show up as curved events in the SODCIG. In this case, subsurface offsets different from zero must be considered during the modeling. Moreover, in the presence of geologic dip, a pre-processing step is needed before modeling. \longcite{biondo2} describes how to generate dip-independent image gathers to be used as initial condition for the upward propagation of source and receiver wavefields. This is achieved by rotating the SODCIGs according to the geological dip in the wavenumber domain.
%\par
% Conceptually, in the prestack exploding reflector modeling, the synthesized areal data is computed by upward propagating source and receiver wavefields using as initial condition one SODCIG. However, I take advantage of the linearity of the wave propagation to combine several experiments into a set of composite records, therefore decreasing the amount of data to migrate.
%\par
%Combining several experiments, though, gives rise to crosstalk during imaging. \longcite{Biondi.sep.124.biondo1} uses a decorrelation distance between SODCIGs sufficiently large to prevent crosstalk. Using this criterion, however, limits the maximum data reduction to the ratio between the number of SODCIGs and the decorrelation distance. Furthermore, he mentions that an additional saving can be achieved by using concepts similar to phase encoding \cite{romero} during the modeling step. This allows us to combine more shots into one set. From now on, upward propagated areal shots combined into one single areal shot will be called a super-areal shot.
%\par
%\longcite{Biondi.sep.124.biondo1,biondo2} uses the two-way wave equation to propagate the wavefields. Here, I use the one-way wave equation to achieve additional computational saving. Furthermore, under the one-way framework, I introduce a phase encoding-like scheme during the modeling to reduce cross-talk when applying the imaging condition.

\subsection{Summary}
In the next section, I will first present the prestack exploding reflector model and discuss the crosstalk problem and attenuation methods. Particularly, I show that in the presence of several reflectors the combination of random phase encoding and a depth-window around selected reflectors virtually eliminates the crosstalk. Finally, I will point out the future work as well as the time-line.

\section{WORK COMPLETED}
\subsection{Prestack exploding reflector modeling}
The basic idea behind the prestack exploding reflector model \cite{Biondi.sep.124.biondo1,biondo2} is that, under the single scattering assumption, the synthesis of 2D source and receiver wavefields at the surface can be described by:
\begin{equation}
\begin{array}{cc}
S(x,\omega) = G(z_{\xi},x_{\xi}-h_{\xi};x,z=0,\omega)*I_s(z_{\xi},x_{\xi},h_{\xi}), \\
R(x,\omega) = G(z_{\xi},x_{\xi}+h_{\xi};x,z=0,\omega)*I_r(z_{\xi},x_{\xi},h_{\xi}), 
\end{array}
\label{eq:eq1}
\end{equation}
where $S(x,\omega)$ is the source wavefield at $z=0$, $R(x,\omega)$ is the receiver wavefield at $z=0$, $I_s(z_{\xi},x_{\xi},h_{\xi})$ is the prestack image used as initial condition to model the source wavefield, $I_r(z_{\xi},x_{\xi},h_{\xi})$ is the prestack image used as initial condition to model the receiver wavefield, and $G(z_{\xi},x_{\xi} \pm h_{\xi};x,z=0,\omega)$ represents an operator which extrapolates the wavefields from the subsurface to the surface. Here, I use one-way extrapolators. Additionally, $h_{\xi}$ is the subsurface offset, $z_{\xi}$ is depth, $\omega$ is the temporal frequency, $x_{\xi}$ is the spatial position of the selected prestack image, and $x$ is the spatial coordinate in the data space coinciding with $x_{\xi}$.
\par
The prestack images used as initial condition for the source and receiver wavefield extrapolation are supposed to be dip-independent gathers computed by changing the dip along the offset direction according to the apparent geological dip \cite{biondo2}. In Appendix A, this transformation is described in detail.
\par
To illustrate how the prestack exploding reflector modeling works, I consider a simple 0.75km deep flat reflector in a 2km/s constant velocity model. Figure \ref{fig:nodipshtpr} shows the result of shot-profile migration using a too slow velocity of 1.8km/s. The front panel represents the zero-subsurface offset section and the panel on the right is a SODCIS taken at 2km. Notice, on the SODCIG panel, that the energy does not focus at zero subsurface offset and the reflector is mapped at a shallower depth.
\activeplot{nodipshtpr}{width=3in,height=2in}{CR}{Shot-profile migration using a slower velocity.}
\par
In Figure \ref{fig:migar1}a is shown the areal source and Figure \ref{fig:migar1}b the areal receiver resulting from the modeling of a single SODCIG as initial condition located at 1.75km. Since the source wavefield is backward propagated, it occurs at negative times. The modeling was performed with the same slower velocity used in migration.
\activeplot{migar1}{height=2.5in}{CR}{Areal data modeled from a single SODCIG as initial condition: a) areal source; and b) areal receivers.}
\par
The modeled areal data is to be used in a subsequent areal shot migration using an updated velocity model. However, in this example I show areal shot migration results with the same slower velocity, just for clarification. Figure \ref{fig:migscig1}a presents the SODCIG taken at 1.75km after the areal shot migration. For comparison, in Figure \ref{fig:migscig1}b is shown the SODCIG resulting from the shot profile migration.
\par
\activeplot{migscig1}{width=4in}{CR}{SODCIG from areal shot migration of the modeled data (a) and from shot profile migration of the original data(b).}
Figures \ref{fig:migacig1}a and \ref{fig:migacig1}b show the angle-domain common-image gather (ADCIG) computed from data in Figures \ref{fig:migscig1}a and \ref{fig:migscig1}b, respectively. Notice that the kinematic information relative to migration with too low velocity was maintained.
\activeplot{migacig1}{width=4in}{CR}{ADCIG from areal shot migration of the modeled data (a) and from shot profile migration of the original data(b). Notice that the kinematic information was maintained.}
\par
To decrease the amount of input data to migration, areal shots need to be combined  into super areal shots after the upward propagation. This can be achieved by regularly selecting individual experiments and adding them up into their respective super areal shots, after being upward propagated. Pairs of super areal source and super areal receiver are to be used as the source and receiver wavefields, respectively, in areal shot migration.
\par
To illustrate the potential data reduction the prestack exploding reflector model can achieve in 3D, I use the geometry of the 3D field data \longcite{valencia133} imaged with shot-profile migration:
\begin{itemize}
\item number of shots: 10,000 (200 along inline  and 50 along crossline)
\item migration aperture: 140 ilines and 330 crosslines
\item model space: 50 ilines and 150 crosslines
\item number of subsurface offsets: 31
\end{itemize}
\noindent{Consider generating 3D areal data at every 10th SODCIGs in the inline and crossline directions. So, the number of super areal shots would be 100. As will be shown later, more realizations improve the quality of the areal shot migration when using random phase encoding during the modeling. Consider 4 realizations, so the final number of super areal shots would be 400.}
\par
 The cost of modeling one areal shot is compared to the cost of the areal migration without the application of imaging condition. Therefore, the total cost to obtain an areal shot migrated image (including the modeling of the areal data) should roughly correspond to the migration of 800 super areal shots. As the aperture of the prestack exploding reflector modeling and areal shot migration spans the whole migrated output space, the relative cost comparison with shot-profile migration is difficult. The cost comparison is easier to be done with plane-wave migration as the migrated domain is the same. For instance, \longcite{shan132} used 2,700 plane-waves in a 3D plane-wave migration. 
\par
In a migration velocity estimation workflow, the areal data is likely to be modeled once. Therefore, in this example, the cost of subsequent migrations should correspond to 1/7 of plane-wave migration.
\par
\longcite{savaphd} pointed out that WEMVA has a huge storage cost because the wavefields need to be saved on disk. As in the prestack exploding reflector model the number of areal shots is about to be 1 or 2 orders of magnitude less than that of the original dataset, the storage costs are expected to be dramatically decreased.  
\par
As I show next, due to the data reduction, crosstalk arises from correlation of unrelated events during imaging. Crosstalk degrades the image, obscuring the kinematic information we are interested in to perform migration velocity analysis.

%\par
%In principle, in the velocity update workflow under the framework of the prestack exploding reflector model, the areal data is modeled once. Nevertheless, as the velocity model is refined and the energy becomes more focused around the zero-subsurface offset, a new modeling can be done parameterized in such a way that more individual experiments are combined together in the same super-areal shot, therefore decreasing the migration cost for later iterations.

\subsection{Crosstalk attenuation}
Wavefield propagation is a linear process. This allows to combine wavefields to decrease the size of input data to migration. However, during imaging, unrelated source and receiver wavefields crosscorrelate, generating crosstalk. Conventionally, to attenuate crosstalk some type of phase encoding technique is applied \cite{romero, Sun01, zhang:1031, liu}. 
\par
In the prestack exploding reflector modeling, let us consider the migration of a super areal data with a single pair of individual experiments. The prestack image is formed by crosscorrelating the source, $\widetilde{S_m}$, and receiver, $\widetilde{R_m}$, wavefields according to:
\begin{equation}
\widetilde{I_m}(z_{\xi},x_{\xi},h_{\xi}) = \sum_{\omega} \widetilde{S_m}^{\star}(z_{\xi},x_{\xi}-h_{\xi},\omega) 
\widetilde{R_m}(z_{\xi},x_{\xi}+h_{\xi},\omega),
\label{eq:eq3}
\end{equation}
where $^{\star}$ represents complex conjugation. If $\widetilde{S_m}(x,\omega)$ and $\widetilde{R_m}(x,\omega)$ are comprised by two summed areal shots, the image $\widetilde{I_m}(z_{\xi},x_{\xi},h_{\xi}) $ will be given by:
\begin{eqnarray}
\lefteqn{\widetilde{I_m}(z_{\xi},x_{\xi},h_{\xi}) = I_1(z_{\xi},x_{\xi},h_{\xi}) + I_2(z_{\xi},x_{\xi},h_{\xi}) + } \nonumber \\
 &  &  \sum_{\omega}S_1^{\star}(z_{\xi},x_{\xi}-h_{\xi},\omega)R_2(z_{\xi},x_{\xi}+h_{\xi},\omega) + \nonumber \\
 &  &  \sum_{\omega}S_2^{\star}(z_{\xi},x_{\xi}-h_{\xi},\omega)R_1(z_{\xi},x_{\xi}+h_{\xi},\omega).
\label{eq:eq4}
\end{eqnarray}
In equation \ref{eq:eq4}, the last two summation terms represent the crosstalk. 
\par
To illustrate the crosstalk problem in the prestack exploding reflector strategy let us make use of a simple 2km/s constant velocity model with two intersecting reflectors, one horizontal and the other dipping 15$^{\circ}$. The dataset is comprised of 200 split-spread shots spaced every 20 m with a maximum offset of 2000 m. Figure~\ref{fig:mig08} shows the shot-profile migration results using the correct velocity. The trace spacing in the migrated result is 20 m and the number of subsurface offsets in the prestack image is 41, ranging from -400m to 400m. The front panel corresponds to the zero-subsurface offset section, the side panel is a SODCIG selected at x=1.5 km. This prestack image is used to model the areal data.
\activesideplot{mig08}{width=3in,height=2in}{CR}{Prestack cube obtained with shot profile migration with the correct velocity.}
\par
Figure \ref{fig:ardata} shows super areal sources (top row) and super areal receivers (bottom row) computed with different distances between SODCIGs. Since the areal source wavefield is backward propagated it exists for negative times. The super areal source labeled a) and super areal receiver labeled d) are modeled from SODCIGs separated every 51 traces (1000 m); the super areal source labeled b) and super areal receiver labeled e) are modeled from SODCIGs separated every 41 traces (800 m); and the super areal source labeled c) and super areal receiver labeled f) are modeled from SODCIGs separated every 11 traces (200 m). These distances in number of traces defines the number of super areal shots to be migrated and, consequently, the migration efficiency. As I show next, this  also influences the intensity of the crosstalk.
\activeplot{ardata}{width=6in}{CR}{Super areal data computed with different distances between SODCIGs. Top row: source wavefields. Bottom row: receiver wavefields. Left column: areal data computed with a SODCIG distance of 51 traces (a and d). Central column: areal data computed with a SODCIG distance of 41 traces (b and e). Right column: areal data computed with a SODCIG distance of 11 traces (c and f).}
\par
Figure \ref{fig:armigW} shows the areal shot migration using cross-correlation imaging condition of the super areal data of Figure \ref{fig:ardata}. The number of subsurface offsets is 41. Figure \ref{fig:armigW}a shows the migration of 51 super areal data from Figure \ref{fig:ardata}a and \ref{fig:ardata}d), Figure \ref{fig:armigW}b shows the migration of 41 super areal shots (Figure \ref{fig:ardata}b and \ref{fig:ardata}e) and Figure \ref{fig:armigW}c shows the migration of 11 super areal shots (Figure \ref{fig:ardata}c and \ref{fig:ardata}f). 
\par
\longcite{Biondi.sep.124.biondo1} shows that using a decorrelation distance greater than twice the maximum subsurface offset, the crosstalks can be eliminated. This explains why migration of 51 super areal shots shows no crosstalk in the SODCIG. The migration of the data with 41 super areal shots shows an acceptable crosstalk in the SODCIG. However, the crosstalk is strong in the SODCIG of the data with 11 super areal shots. Notice that for all of the three results a strong crosstalk occurs at zero-subsurface offset (x,z-domain). 

\activeplot{armigW}{width=6in}{CR}{Areal shot migration using conventional cross-correlation imaging condition of 51 super areal shots (a), 41 super areal shots (b) and 11 super areal shots (c).}
\par
The crosstalk described in the previous paragraph have two distinct origins. Crosstalk in the zero-subsurface offset section results from the cross-correlation of events in the source wavefield with that in the receiver wavefields modeled by the same SODCIG, which herein I call type-1. Crosstalk in non-zero subsurface offsets is related to cross-correlation of events in the source wavefields with that in the receiver wavefields modeled by different SODCIGs, which herein I call type-2.
\par
%\par
%Because we performed migration with the correct velocity and the energy is concentrated at zero-offset, the pre-processing step of rotating the SODCIGs according to the geological dip is not necessary before modeling. Therefore, to model the areal source and receiver, SODCIGs are upward propagated using the one-way wave equation (equation~\ref{eq:eq1}) without any pre-processing step. The result of one modeling experiment is collected at the surface as an areal shot and assigned to its specific super areal shot. 

%The cross-correlation occurs at times different from zero propagation time. Therefore, one possible strategy is to use a time-windowed imaging condition about zero propagation time.
\longcite{biondo2} and \longcite{guerra} show that type-1 crosstalk is strongly attenuated by windowing the wavefields around zero propagation time. This method yields good results if the migration velocity is accurate. In this case, reflectors `explode' at zero time \cite{iei}. However, the window size must be chosen with care if the migration velocity is inaccurate to include the lags where the wavefields correlate.
\par
Here, I describe the use of random phase encoding during the modeling as a way to attenuate crosstalk. This strategy can be combined with the selection of specific reflectors to be modeled for the case of more complex geology. This procedure will be explained at the end of this section.

%. \longcite{romero} attenuate these migration artifacts by multiplying the source and the receiver wavefields by a function of frequency and space. By doing so, the crosstalks are dispersed (random phase encoding) throughout the image after stacking over frequency, during imaging, or shifted outside of the image space (linear phase encoding). Here, we use a similar approach while modeling the areal shots.
%\subsection{Crosstalk attenuation with time-windowed imaging condition}
%The fact that the areal source wavefield also contains reflections originates crosstalk. According to the exploding reflector model \cite{loewenthal}, reflectors explode at time zero. Therefore, the reflections in the areal source wavefield should be focused at zero time of wavefield propagation. We can use this information to restrict the cross-correlation in a time-window about zero time of wavefield propagation. As \longcite{biondo2} shows, the cross-correlation of areal wavefields within a small window around zero time of wavefield propagation can diminish the type-1 crosstalk at zero-subsurface offset. Equation \ref{eq:eq5} shows that, using the time-windowed imaging condition, crosstalk at the zero-subsurface offset only occurs between events having similar traveltimes.
%\begin{eqnarray}
%\lefteqn{\widetilde{I}(z_{\xi},x_{\xi},h_{\xi}) = I_1(z_{\xi},x_{\xi},h_{\xi}) + I_2(z_{\xi},x_{\xi},h_{\xi}) + }   \nonumber \\
% &  &  \sum_{-\Delta t \leq t=0 \leq \Delta t}S_1^{\star}(z_{\xi},x_{\xi},t)R_2(z_{\xi},x_{\xi},t) + \nonumber \\
% &  &  \sum_{-\Delta t \leq t=0 \leq \Delta t}S_2^{\star}(z_{\xi},x_{\xi},t)R_1(z_{\xi},x_{\xi},t).
%\label{eq:eq5}
%\end{eqnarray}
%\par
%Applying this imaging condition using one-way propagators in frequency domain, requires storing both wavefields in memory for all the frequencies, which increases the memory requirements. On the other hand, applying the imaging condition for a small number of time samples saves much computational effort. Figure \ref{fig:armigT} shows the areal shot migration using the time-windowed cross-correlation imaging condition from equation \ref{eq:eq4}. Again, the number of subsurface offsets is 41. Figure \ref{fig:armigT}a is the migration of 51 super areal shots, Figure \ref{fig:armigT}b is the migration of 41 super areal shots and Figure \ref{fig:armigT}c is the migration of 11 super areal shots. Notice how the crosstalk in the zero-subsurface offset sections have been largely attenuated in all the migration results. 
%\plot{armigT}{width=6in}{Areal shot migration using time-windowed imaging condition of 51 super areal shots (a), 41 super areal shots (b) and 11 super areal shots (c). Compare with Figure \ref{fig:armigW}.}
%\par
%For the migration in Figure \ref{fig:armigT}c, the time-windowed imaging condition, however, does not attenuate the severe type-2 crosstalk in the SODCIG. Moreover, wavefields cross-correlate at zero lags in space and time \cite{iei,sava:S209} when the illumination is sufficiently good and is used the correct migration velocity. However, when the migration velocity is not accurate, the peak of the cross-correlation between wavefields does not occur at zero lags. Therefore, care must be taken when choosing the time interval to apply the imaging condition when migration velocity is inaccurate. By choosing a sufficiently wide window, events migrated with wrong velocity can be adequately imaged. However, the wider the time-window, the less effective is the crosstalk attenuation.

\subsection{Crosstalk attenuation with random phase encoding}
Ideally, the areal shots pertaining to the same super areal shot should be uncorrelated. A possible way to decrease the correlation between areal shots is by encoding them with a random phase function in such a way that the crosstalk is dispersed throughout the image. \longcite{romero} linearly combine shots, using as encoding function $a(x_{\xi},\omega) = e^{i f(x_{\xi},\omega)}$. $f(x_{\xi},\omega)$ is a random sequence varying between $(-\pi,\pi)$.
\par 
The modeling with $f(x_{\xi},\omega)$-phase encoding synthesizes data according to
\begin{equation}
\begin{array}{cc}
S_a(x,\omega) = a(x_{\xi},\omega) G(z_{\xi},x_{\xi}-h_{\xi};x,z=0,\omega)*I_s(z_{\xi},x_{\xi},h_{\xi}), \\
R_a(x,\omega) = a(x_{\xi},\omega) G(z_{\xi},x_{\xi}+h_{\xi};x,z=0,\omega)*I_r(z_{\xi},x_{\xi},h_{\xi}), 
\end{array}
\label{eq:eq6}
\end{equation}
and the super areal shots of random phase encoded areal sources and receivers, $\widetilde{S_n^a}(x,\omega)$ and $\widetilde{R_n^a}(x,\omega)$, respectively, are given by:
\begin{eqnarray}
\widetilde{S_n^a}(x,\omega)& = &\sum_{n=1}^k  \sum_{i=n,k,N} S_i^a(x,\omega), \nonumber \\
\widetilde{R_n^a}(x,\omega)& = &\sum_{n=1}^k  \sum_{i=n,k,N} R_i^a(x,\omega), 
\label{eq:eq7}
\end{eqnarray}
\par
Figures \ref{fig:arrphw}a and \ref{fig:arrphw}b show $(x_{\xi},\omega)$-random phase encoded super areal source and receiver wavefields, respectively. 
\par
\par
\activeplot{arrphw}{width=4in}{CR}{Areal source (a) and receiver (b) wavefields computed with random phase encoding, $a(x_{\xi},\omega)$.}
\par
The image, $\widetilde{I_m}(z_{\xi},x_{\xi},h_{\xi})$, after areal shot migration of super areal data comprised by two randomly phase encoded areal sources and receivers is given by:
\newpage
\begin{eqnarray}
\lefteqn{\widetilde{I_m}(z_{\xi},x_{\xi},h_{\xi}) = I_1(z_{\xi},x_{\xi},h_{\xi}) + I_2(z_{\xi},x_{\xi},h_{\xi}) + } \nonumber \\
 &  &  \sum_{\omega}a_1^{\star}a_2 S_m^{a_1^{^\star}}(z_{\xi},x_{\xi}-h_{\xi},\omega)R_m^{a_2}(z_{\xi},x_{\xi}+h_{\xi},\omega) + \nonumber \\
 &  &  \sum_{\omega}a_2^{\star}a_1 S_m^{a_2^{^\star}}(z_{\xi},x_{\xi}-h_{\xi},\omega)R_m^{a_1}(z_{\xi},x_{\xi}+h_{\xi},\omega).
\label{eq:eq8}
\end{eqnarray}
\par
The cross-correlation between the phase encoding functions $a_1(x_{\xi},\omega)$ and $a_2(x_{\xi},\omega)$ in the summation terms must provide random values that disperse the crosstalk throughout the image.
\par
Figure \ref{fig:arm11phw} shows the areal shot migration of 11 super areal shots of  $(x_{\xi},\omega)$-randomly encoded areal data. The reduction of the  type-2 crosstalk in the SODCIG is significant. However, the type-1 crosstalk in the zero-subsurface offset is much less affected. This crosstalk originates from uncorrelated reflections modeled from the same SODCIG as described before. A possible way to attenuate this crosstalk is to decorrelate the exploding reflectors by using a different random function $a(z_{\xi},x_{\xi},\omega) = e^{i f(z_{\xi},x_{\xi},\omega)}$. 
\par
The modeling with $a(z_{\xi},x_{\xi},\omega)$-random phase encoding synthesizes data according to
\begin{equation}
\begin{array}{cc}
S_a(x,\omega) = \varepsilon a(z_{\xi},x_{\xi},\omega) G(z_{\xi},x_{\xi}-h_{\xi};x,z=0,\omega)*I_s(z_{\xi},x_{\xi},h_{\xi}), \\
R_a(x,\omega) = \varepsilon a(z_{\xi},x_{\xi},\omega) G(z_{\xi},x_{\xi}+h_{\xi};x,z=0,\omega)*I_r(z_{\xi},x_{\xi},h_{\xi}).
\end{array}
\label{eq:eq9}
\end{equation}
In equation \ref{eq:eq9}, the $\varepsilon$ term is used to control the randomness of the wavefields. As \longcite{romero} pointed out, the encoding results deteriorate as the amount of encoded data increases. By letting $\varepsilon=1$, the migration results become too noisy. A typical value for $\varepsilon$ is 0.1. %Figure \ref{fig:arrphz}a and \ref{fig:arrphz}b show a super areal shot of random phase encoded areal source and receiver wavefields, respectively.

%\plot{arrphz}{width=5in}{Areal source (a) and receiver (b) wavefields computed with random phase encoding, $a(z_{\xi},x_{\xi},\omega)$.}

Figure \ref{fig:arm11phz} show the areal shot migration results of 11 super areal shots of randomly encoded areal data using $(z_{\xi},x_{\xi},\omega)$-randomly encoded areal data. The comparison with Figure \ref{fig:arm11phw} and Figure \ref{fig:armigW}c shows that encoding in $(z_{\xi},x_{\xi},\omega)$ yields a significant reduction of the type-1 crosstalk in the zero-subsurface offset.
\par
The dispersed noise is not negligible. Using more realizations of $(z_{\xi},x_{\xi},\omega)$-randomly phase encoded data, the noise amplitudes can be largely attenuated. Figure \ref{fig:arm11phz4} shows the superior result of migrating 4 different realizations of 11 $(z_{\xi},x_{\xi},\omega)$-randomly encoded data. Besides virtually eliminating the crosstalk, the speckled noise has much lower amplitudes than in Figure \ref{fig:arm11phz}.

\activesideplot{arm11phw}{width=3in,height=2in}{CR}{Areal shot migration of randomly encoded data. 11 super areal shots of $(x_{\xi},\omega)$-randomly encoded areal data were migrated. Compare with Figure \ref{fig:armigW}c.}
\activesideplot{arm11phz}{width=3in,height=2in}{CR}{Areal shot migration of randomly encoded data. 11 super areal shots of $(z_{\xi},x_{\xi},\omega)$--randomly encoded areal data were migrated. Compare with Figure \ref{fig:arm11phw}.}
\activesideplot{arm11phz4}{width=3in,height=2in}{CR}{Areal shot migration of randomly encoded data. 4 groups of 11 super areal shots of $(z_{\xi},x_{\xi},\omega)$--randomly encoded areal data were migrated. Compare with Figures \ref{fig:arm11phw} and \ref{fig:arm11phz}.}

\par
In the previous examples, I used simple models consisting of two reflectors to illustrate crosstalk attenuation methods. The crosstalk originated by the correlation of reflections modeled from the same SODCIG are fairly well attenuated by the random phase encoding in  $(z_{\xi},x_{\xi},\omega)$. However, in the presence of several reflectors this type of crosstalk is severe and using only random phase encoding techniques do not perform successfully. So, I introduce the horizon-based prestack exploding reflector model which is the basis for the hybrid approach of wave-equation velocity inversion. I illustrate the horizon-based prestack exploding reflector model using the sediment part of Sigsbee2a \cite{paffenholz02} dataset.

\subsection{Prestack exploding reflector modeling - {S}igsbee2a}
I use part of the Sigsbee2a dataset \cite{paffenholz02} to show how the crosstalk attenuation methods described previously perform when more reflectors are present. I selected the sediment part to facilitate the analysis. Again, the correct migration velocity was used in all wavefield extrapolations. The areal shots were generated using every $15^{th}$ SODCIG as initial condition for the modeling. Therefore, 15 super areal shots were input to migration. Here, I show the results of areal shot migration using the conventional imaging condition in frequency and migration of $(z_{\xi},x_{\xi},\omega)$-randomly encoded areal data. 
\par
Figure \ref{fig:sigshp} shows the shot profile migration of the sediment portion of Sigsbee2a data. The main feature in the zero-subsurface offset  section is represented by a faulted syncline with reflectors dipping in different directions and a fault on the left side. Figure \ref{fig:signicW} shows the areal shot migration using the frequency-domain imaging condition of the non-encoded data. This image is dominated by type-1 and type-2  crosstalk represented by the extra `reflectors' when comparing with Figure \ref{fig:sigshp}, making the true reflectors unrecognizable. 
\activesideplot{sigshp}{width=3in,height=2in}{CR}{Shot profile migration of the sediment portion of Sigsbee2a dataset.}
\activesideplot{signicW}{width=3in,height=2in}{CR}{Areal shot migration of the non-encoded data in the sediment portion of Sigsbee2a using frequency-domain imaging condition. Compare with Figure \ref{fig:sigshp}.}
%\par
%Figure \ref{fig:signicT} shows the areal shot migration using the time-windowed imaging condition of the non-encoded data. Using the time-windowed imaging condition has the ability to attenuate the type-1 crosstalk in the zero-subsurface offset. However, as can be seen in the subsurface-offset gather, the type-2 crosstalk is still severe.
%\sideplot{signicT}{width=3.5in}{Areal shot migration of the modeled data in the sediment portion of Sigsbee2a using time-windowed imaging condition of the non-encoded data. Compare with Figures \ref{fig:sigshp} and \ref{fig:signicW}.}
\par
I show the areal shot migration of 4 realizations of both $(x_{\xi},\omega)$-random encoding and $(z_{\xi},x_{\xi},\omega)$-random encoding in  Figures \ref{fig:signphw4} and \ref{fig:signphz4}, respectively. In general, type-2 crosstalk is reasonably attenuated in the subsurface-offset gathers. Type-1 crosstalk in the zero-subsurface-offset section, however, is still severe even in the migration of $(z_{\xi},x_{\xi},\omega)$-randomly encoded data. This is somewhat related to the degradation of the image quality when the number of encoded shots increases, as observed by \longcite{romero}. 
\par
\activesideplot{signphw4}{width=3in,height=2in}{CR}{Areal shot migration of 4 realizations of $(x_{\xi},\omega)$-randomly encoded data in the sediment portion of Sigsbee2a.}
\activesideplot{signphz4}{width=3in,height=2in}{CR}{Areal shot migration of 4 realizations of $(z_{\xi},x_{\xi},\omega)$-randomly encoded data in the sediment portion of Sigsbee2a.}
\par
To circumvent the crosstalk problem in Sigsbee2a example, I use a sparse representation of the initial conditions by selecting only a few reliable, high amplitude reflectors in the original SODCIGs to model. To justify this sparse representation of the migrated image, remind that for velocity model updating we can use some key reflectors, generally high signal-to-noise ratio (S/R), to extract the residual-moveout information \cite{stork01,kosloff:1511,marfurt}. Herein I call this procedure horizon-based prestack exploding reflector model. 
%\par
%Ideally, these horizons must limit layers with significant differences in the acoustic properties and they are generally chosen along reflectors with good signal-to-noise ratio (S/R). They are also used to conveniently  constrain the null space of reflection-tomography problems by parameterizing the model space with layers \cite{kosloff:1511}.  Additionally, this representation agrees with the formulation of the velocity determination problem in a target oriented way, considering that a reasonable velocity model is available for shallower depths. Finally, dip information from horizons can be used in further investigations on the application of different imaging conditions \cite{sava02}. Consequently, for velocity update purposes, just the migrated image along these key horizons is sufficient.
\par
Figure \ref{fig:sigpick2} shows two horizons I choose to constrain the modeling. As the imaged reflectors are characterized by a wavelet stretched to depth, I use a small window around the horizon to select the reflector amplitude to input to modeling.
\activesideplot{sigpick2}{width=3.4in,height=2in}{NR}{Two picked horizons used to constrain the prestack exploding reflector modeling.}
\par
 Figure \ref{fig:sigphz4} shows the areal shot migration of 4 realizations of $(z_{\xi},x_{\xi},\omega)$-randomly encoded areal data constrained by the horizons. The crosstalk is almost completely eliminated. 
\activesideplot{sigphz4}{width=3in,height=2in}{CR}{Areal shot migration of 4 realizations of $(z_{\xi},x_{\xi},\omega)$-randomly encoded areal data using horizons.}
\par
Before proceeding to the next section, where I point out the work to be done, I give a rationale for using the horizon-based prestack exploding reflector model in wave-equation velocity inversion.
\par
The main application of WEMVA is velocity model update under complex overburden, for instance, irregular salt bodies. To obtain accurate updates, it is crucial that the velocity model for the overburden is reasonably defined, such that observed defocusing effects are due to velocity errors in the underlying sediments. 
\par
The velocity contrast between salt bodies and sediments is commonly strong. But in general, seismic velocity in sedimentary layers is horizontally smooth. This is because sediment deposition is mainly controlled by the gravitational force \cite{sediment}, so the main velocity variations occur in the vertical direction. Strong lateral velocity variations can occur across big geological faults which put rocks of different ages in lateral contact. Even so, velocity varies smoothly within the faulted blocks.
\par
Apart from the geological reasons, velocity uncertainties are generally big because of complex overburden causing uneven illumination of the subsurface, limited acquisition geometry and presence of noise. These factors make the velocity inversion problem ill-posed and the use of regularization, penalizing roughness of the velocity model, is mandatory \cite{bob106}. Additionally, velocity resolution is limited by the Fresnel zone size \cite{thore:593}. For these reasons, velocity interpreters generally prefer smooth velocity models in the sediment portion \cite{guo:1210}. This characteristic is frequently observed in most of the depth migration velocity models from Gulf of Mexico.
\par
Theoretically, by using a few reflectors, the slowness perturbations are supposed to be smoother than that resulting from the inversion of the complete image perturbation. However, given that velocity models are generally smooth below the complex overburden, I expect this possible additional smoothness will not be a problem. In the other hand, these few reliable reflectors can possibly give rise to more reliable velocity updates. These are questions I expect to address in this research.

\section{WORK TO BE DONE}
In the section "Work Completed", I explained the prestack exploding reflector model and discussed the crosstalk problem. Finally, I introduced the horizon-based approach to prestack exploding reflector model to make the method applicable to field data. The results on a realistic synthetic seismic data confirms that the method is promising. 
\par
To further illustrate the good results of the horizon-based prestack exploding reflector model, Figure \ref{fig:sigphz4ang} shows the data from Figure \ref{fig:sigphz4} transformed to the angle domain. For comparison, the ADCIG of the original shot profile migration is shown in Figure \ref{fig:sigshtang}. Notice that the kinematic information is maintained. 
\par
\activesideplot{sigphz4ang}{width=3in,height=2in}{CR}{Angle-domain data computed from the areal shot migration of 4 realizations of $(z_{\xi},x_{\xi},\omega)$-randomly encoded areal data constrained by the horizons.}
\activesideplot{sigshtang}{width=3in,height=2in}{CR}{Angle-domain data computed from the original shot-profile migration.}
\par
As the main goal of this research is to explore WEMVA using the prestack exploding reflector concept, in the following I outline the general structure of WEMVA under the framework of the horizon-based prestack exploding reflector model by exploring the topics below.
\par
Some important questions have to be addressed, like:
\begin{itemize}
\item WEMVA in the presence of noise,
\item WEMVA in the presence of illumination problems.
\end{itemize}
Other important issue is how crosstalk attenuation performs when using:
\begin{itemize}
\item different phase encoding schemes,
\item alternative imaging conditions.
\end{itemize}
These questions can be properly addressed using 2D algorithms. The 3D extension of the algorithms will enable to tackle more interesting problems from the exploration point of view. So, other developments to be made are to extend to 3D:
\begin{itemize}
\item the rotation of SODCIGs,
\item the prestack exploding reflector modeling algorithm,
\item the areal shot migration.
\end{itemize}
{As I showed, the horizon based approach to the prestack exploding reflector model plays an important role in crosstalk attenuation. So, to incorporate this approach in the 3D workflow it will be necessary to:
\begin{itemize}
\item investigate 3D residual prestack migration of the picked horizons,
\item investigate 3D and 5D automatic horizon picking strategies.
\end{itemize}
\par

\subsection{WEMVA in the presence of noise}
Similarly to DSO, one potential problem WEMVA may face is the presence of coherent noise. \longcite{mulder:1184} showed that velocity inversion with DSO is very sensitive to the presence of multiples and converted waves. To obtain reasonable results, a 2D-wavenumber filter was incorporated into the DSO objective function. 
\par
Another issue is how WEMVA will perform in the presence of the randomness of the wavefields when using random phase encoding. Or, how the scattered wavefield computed by the interaction of the reference wavefields with velocity perturbations is affected by the randomness of the wavefields? Will the image perturbations, resulting from the migration of the scattered wavefield, be stable enough to provide reliable velocity updates?
\par
To deal with the noise problem, a possibility I want to explore in this research is to apply WEMVA on a modeled dataset with only on a few select reflectors with higher S/R.

\subsection{WEMVA in the presence of illumination problems}
Reflectors do not focus around zero-subsurface offset under uneven illumination. When the illumination problem is too severe, shadow zones occur and migration artifacts dominate the migration result \cite{valencia133}. In the presence of illumination problems, reflectors in ADCIGs can show departures from flatness not related to velocity errors. 
\par
One solution to uneven illumination is to pose imaging as an inverse problem \cite{Mclapp122,valencia133}. Because of ill-posedness, this inverse problem needs regularization. However, regularization in the prestack inversion enforces reflectors to be focused at zero-subsurface offset or, alternatively, flatness in the ADCIG. When the velocity is not accurate this regularization can destroy the kinematic information necessary to velocity updates.
\par
\longcite{savaphd} shows that the angular coverage in the subsurface influences the slowness perturbation. The example of WEMVA applied on Sigsbee clearly shows that the quality of the velocity update deteriorates below the salt body. \\
\par
\vspace{0.2in}
I will investigate the potential of WEMVA constrained to some specific horizons. This characterizes a hybrid approach to velocity analysis using wave-equation methods. Ideally, these horizons must limit layers with significant differences in the acoustic properties and correspond to   reflectors with good S/R. 
\par
The horizon-based technique is used to conveniently  constrain the null space of reflection-tomography problems by parameterizing the model space with layers \cite{stork01,kosloff:1511,marfurt}. Additionally, this representation agrees with the formulation of the velocity determination problem in a target oriented way, considering that a reasonable velocity model is available for shallower depths. Consequently, for velocity update purposes, just the migrated image along these key horizons is sufficient.

\subsection{Using different phase encoding schemes}
One drawback of applying random phase encoding is that it makes the reflectors randomly ``explode'' at times different from the zero time of the wavefield propagation. Consequently, the time-windowed imaging condition can not be applied.
\par
The linear phase encoding \cite{romero} aims to shift the crosstalk out of the imaging domain. The crosstalk starts occurring at propagation times greater than the time-shift applied to the wavefields, enabling us to use the time-windowed imaging condition. Another phase encoding schemes to be tested are the mixed-phase encoding and the sign-opposite phase encoding \cite{Sun01,zhang:1031}.

\subsection{Using alternative imaging conditions}
Type-2 crosstalk is similar to the one described by \longcite{sava02}. \longcite{sava02} uses an imaging condition that cross-correlates decomposed source and receiver wavefields as a function of an a priori known local reflector slope at every position and time. If the reflector slope is not available, it is necessary to loop over a range of possible dip angles in a vicinity of the expected reflector dip. This procedure yields good results, but at the expense of a high  computational cost. The horizon-based approach can yield the dip angle information to implement this imaging condition in a more efficient way. \\
\par
\vspace{0.2in}

The research on phase encoding schemes and alternative imaging conditions is extremely interconnected and I expect to determine the most appropriate combination between them.

\subsection{3D extension of rotation, modeling and migration algorithms}
The 3D extension of the prestack exploding reflector modeling and areal shot migration algorithms along with code parallelization is, in principle, trivial. However, more understanding is needed to formulate the 3D SODCIG rotation in the case of inaccurate migration velocity. That is because, except for the case of common-azimuth migration \cite{biondi:1822}, when the migration velocity is inaccurate source and receiver rays do not necessarily cross and the co-planarity condition is not respected.

\subsection{Picking strategies}
As already mentioned, using a horizon-based approach to prestack exploding reflector model virtually eliminates crosstalk. Its efficiency, however, may depend on the S/R of the reflectors and on the vertical distance between them. There seems to be, therefore, a vertical decorrelation distance between reflectors. Its determination is a work still to be done.
\par
For the 2D case of a prestack image obtained with an accurate velocity, just two-dimensional information is needed. If the 2D prestack images are computed with an inaccurate velocity model, however, it is necessary to expand the dimensionality of the picks, using the subsurface-offset dimension to correctly capture the kinematic deviations from the zero-subsurface offset. In 3D, the picks are five-dimensional.  To implement the horizon-based approach in 3D, automatic picking strategies will be needed, at least to pick the original depth migrated data.
\par
Horizon picking after velocity updates are also needed. I intend to use 3D residual prestack depth migration \cite{sava:634} to transform the original 5D picks into updated picks after velocity updates.
\par
3D residual prestack depth migration accounts for remapping, in the wavenumber domain, the vertical wavenumber axis, $k_z{_0}$, of the migrated image into the vertical wavenumber axis, $k_z$, of the residual migrated image according to:
\begin{eqnarray}
\lefteqn{k_z = \sqrt{\frac{1}{\rho^2}\frac{\left (k^2_{z_0}+ \left |{\bf k_m} \right |^2 \right )\left (k^2_{z_0}+ \left |{\bf k_h} \right |^2 \right )}{4k^2_{z_0}}-\frac{1}{4}\left ({\bf k_m}+{\bf k_h} \right )\left ({\bf k_m}+{\bf k_h} \right )}} \nonumber \\
 & & \sqrt{\frac{1}{\rho^2}\frac{\left (k^2_{z_0}+ \left |{\bf k_m} \right |^2 \right )\left (k^2_{z_0}+ \left |{\bf k_h} \right |^2 \right )}{4k^2_{z_0}}-\frac{1}{4}\left ({\bf k_m}-{\bf k_h} \right )\left ({\bf k_m}-{\bf k_h} \right )}
\end{eqnarray}
where ${\bf k_m}$ and ${\bf k_h}$ are the midpoint and subsurface-offset wavenumber vectors, respectively, and $\rho = \frac{v_0}{v}$ is the ratio between the initial velocity, $v_0$, and the updated velocity, $v$.

\section{Tools and data}

\subsection{Software}
In the proposed method an initial migration must be available. To generate the initial images I plan to use the one-way wave equation shot-profile migration, coded by Brad Artman. This algorithm has been already used in three PhD thesis to migrate 3D data \cite{Artman,Alvarez130,valencia133}. This algorithm is therefore considered quite stable and reliable. To perform WEMVA I will run the codes Paul Sava used in his thesis.
\par  
2D prestack exploding reflector modeling and 2D areal shot migration algorithms are already tested. Their 3D extension is still a work to be done along with code parallelization.
\par
SEPlib and SEP3D programs complete the suite of programs needed for my research.

\subsection{Hardware}
Even with the data reduction achievable by using the prestack exploding reflector model, 3D modeling and migration algorithms are computationally intensive. Moreover, 3D WEMVA is a big challenge because of the storage and computational requirements.
\par
The computational resources for this research will be provided by SEP's computer power. Currently, SEP has 1 shared memory machine with 16 nodes and 32Gb RAM (pompei), 3 shared memory machines with 8 nodes and 8Gb RAM (mad, sad and glad), 1 Linux PC cluster with 32 nodes and 2GB RAM (sep200's) and 1 Linux PC cluster with 32 nodes and 4/8GB RAM (sep400's). Two additional Linux PC clusters are expected to be configured soon.
\par
I expect to count on this computational resources to develop my research. 

\subsection{Seismic Data}
To prove the concepts and test the programs, besides of modeling specific datasets, SEP data library has interesting data to be used. 
\par
I will use Sigsbee2a to test 2D algorithms and 3D SEG-EAGE salt model for the 3D extensions.
\par
Petrobras has elected a challenging offshore 3D data acquired over an area of the Santos Basin, Brazil, to benchmark depth-imaging algorithms and velocity-model estimation strategies. The data was depth migrated with a Kirchhoff algorithm and the velocity was estimated using reflection tomography. In spite of the effort the contractor expended to produce an accurate velocity model, Petrobras evaluates that there are opportunities for improvements. Especially, imaging of steep salt-flanks and characterization of a possible allochtonous salt body in the middle of the section along with sediments below.
\par
Along with the seismic data itself, the depth migration velocity model used to migrate data from Figure \ref{fig:fig01}, observer's report and 9 time-migrated horizons are available. This data could to be used in this research.
\activeplot{fig01}{width=5in}{NR}{Section of the depth migrated data from Santos Basin, Brazil.}
\par
If this data proves not be adequate for the purposes of this research, because of lack of sediments below the salt, Elf dataset in the North sea also presents challenges to velocity definition. It is also a well known data set and it has been extensively  used by SEP students.

\newpage

\section{TIME-LINE}
\begin{enumerate}
\item \underline{Summer 2007-2008}\\
Further investigation on crosstalk attenuation and imaging conditions\\
2D tests of WEMVA using horizon-based prestack exploding reflector model \\
Investigation on picking strategies \\
3D migration of SEG-EAGE salt model\\
Course requirements
\item \underline{Autumn 2008-2009}\\
Further investigation on crosstalk attenuation and imaging conditions\\
Investigation on picking strategies \\
3D extension of the areal shot migration\\
3D extension of the prestack exploding reflector model\\
Test 3D algorithms SEG-EAGE 3D salt model\\
Course requirements
\item \underline{Winter 2008-2009}\\
3D extension of the areal shot migration\\
3D extension of the prestack exploding reflector model\\
Investigation on picking strategies \\
Test 3D algorithms SEG-EAGE 3D salt model\\
Course requirements
\item \underline{Spring 2008-2009}\\
3D migration of field data\\
Picking of the field data migration \\
Course requirements
\item \underline{Summer 2008-2009}\\
3D prestack exploding reflector modeling of field data\\
3D velocity model update of field data\\
3D areal shot migration of field data\\
Start writing thesis
\item \underline{Autumn 2009-2010}\\
3D prestack exploding reflector modeling of field data\\
3D velocity model update of field data\\
3D areal shot migration of field data\\
Writing thesis
\item \underline{Winter 2009-2010}\\
Finalize deliverables for graduation\\
\end{enumerate}

\newpage

\APPENDIX{A}
\subsection {Rotation of Common Image Gathers}
In shot-profile migration, SODCIGs are computed by applying the multi-offset imaging condition \cite{rickett}, in which the wavefields are shifted horizontally prior to correlation. If the migration velocity is inaccurate, the image point does not coincide with the point where source and receiver rays intersect \cite{biosymes}. Consequently, in the presence of dip, source- and receiver-ray end points do not lie on the imaged reflector. To correct for the effects these imperfections have on the prestack exploding reflector modeling results \longcite{biondo2} proposes to apply a preprocessing step before modeling.
\par
In the preprocessing step, a rotation is applied to the CIGs, changing the dips in the $z-h$ plane (depth-offset) based on the apparent dips in the $z-x$ plane (depth$-$x coordinate). By rotating the CIGs, the upward propagation of source and receiver wavefields starts along the apparent source and receiver ray local directions. The rotation can be implemented in the wavenumber domain ($k_z-k_x-k_h$). 
\par
As Figure \ref{fig:cig-rot} shows, local directions of the source and receiver rays are given by $\theta _s = \alpha - \gamma$ and $\theta _g = \alpha + \gamma$, respectively, where  $\alpha$ is the apparent dip and $\gamma$ is the apparent reflection angle. %I call them apparent because the situation represented in Figure \ref{fig:cig-rot} corresponds to data migrated with a velocity slower than the correct, so the angles are not the true ones. 
The aperture angle is related to the dips along the offset direction by \cite{Sava.sep.103.paul2}:
\begin{equation}
\tan \gamma = -\frac{k_h{_{_\xi}}}{k_z{_{_\xi}}},
\label{eq:gamma}
\end{equation}
\noindent
{where $k_z{_{_\xi}}$ is the vertical wavenumber and $k_h{_{_\xi}}$ is the offset wavenumber. The apparent geological dip is computed as} 
\begin{equation}
\tan \alpha = -\frac{k_x{_{_\xi}}}{k_z{_{_\xi}}},
\label{eq:alpha}
\end{equation}
\noindent
{where $k_x{_{_\xi}}$ is the midpoint wavenumber. Using equations \ref{eq:gamma} and \ref{eq:alpha} and trigonometric expressions, $\tan \theta _s$ and $\tan \theta _g$ are computed as} 
\begin{equation}
\tan \theta _s (k_z{_{_\xi}},k_x{_{_\xi}},k_h{_{_\xi}})= k_z{_{_\xi}} \frac{k_h{_{_\xi}}-k_x{_{_\xi}}}{k_z{_{_\xi}}^2+k_x{_{_\xi}} k_h{_{_\xi}}},
\label{eq:thetas}
\end{equation}
\begin{equation}
\tan \theta _g (k_z{_{_\xi}},k_x{_{_\xi}},k_h{_{_\xi}})= -k_z{_{_\xi}} \frac{k_h{_{_\xi}}+k_x{_{_\xi}}}{k_z{_{_\xi}}^2-k_x{_{_\xi}} k_h{_{_\xi}}},
\label{eq:thetag}
\end{equation}
\activesideplot[]{cig-rot}{width=3.0in}{NR}{Geometric representation of the rotation of an image event from horizontal subsurface offset to geologic-dip subsurface offset (from \longcite{biondo2}).}

The rotation of the reflectors in the CIGs is realized via a simple remapping of subsurface-offset wavenumber according to
\begin{equation}
\hat{k}_h{_{_\xi}}^s = -k_z{_{_\xi}} \tan \theta _s
\label{eq:mapkhs}
\end{equation}
\begin{equation}
\hat{k}_h{_{_\xi}}^g = k_z{_{_\xi}} \tan \theta _g
\label{eq:mapkhg}
\end{equation}
\noindent
{where $\hat{k}_h{_{_\xi}}^s$ and $\hat{k}_h{_{_\xi}}^g$ are the wavenumbers corresponding to the new subsurface-offset axes $\hat{h}{_{_\xi}}^s$ and $\hat{h}{_{_\xi}}^g$, respectively. 
%After rotation, two new images gathers are generated. One refers to the initial condition for the modeling of the source wavefield and the other to the initial condition for the modeling of the receiver wavefield. By using these rotated gathers as initial conditions, the modeled areal data is dip-independent.
\par
To illustrate the importance of the SODCIG rotation, I use a simple model of a $25^{\circ}$-dipping reflector. Data were migrated with a 10 percent slower velocity. The SODCIG at 2km was used as non-rotated (Figure \ref{fig:dipshtpr1}) initial condition and after rotation  (Figure \ref{fig:datarot}) to generate areal data. For the non-rotated case, the initial condition is the same for the source and receiver wavefield modeling.  
\activeplot{dipshtpr1}{width=4in,height=1.75in}{CR}{Initial condition to model the source wavefield and the receiver wavefield for the non-rotated case.}
\activeplot{datarot}{width=4in,height=3.5in}{CR}{SODCIG after rotation according to the geological dip used as initial condition to model the source wavefield (a) and the receiver wavefield (b).}
\par
The modeled areal data were migrated using the correct velocity (2km/s). After migrating with the correct velocity, the energy is expected to focus in the up-dip direction, since it was initially under-migrated \cite{thore:840}. Figures \ref{fig:dipnorot} and \ref{fig:diprot} show the areal shot migration using areal data modeled from non-rotated SODCIG and areal data modeled from rotated SODCIG, respectively. Areal migration using data modeled from rotated SODCIG correctly positions the event and focus the reflection around the zero-subsurface offset as can be seen in the right panel. Areal migration using data modeled from non-rotated SODCIG does not laterally shift the reflector and the reflection energy shows some curvature in the subsurface-offset gather. 
\activesideplot{dipnorot}{width=4in,height=2in}{CR}{Areal shot migration using areal data modeled from non-rotated SODCIG of Figure \ref{fig:dipshtpr1}.}
\activesideplot{diprot}{width=4in,height=2in}{CR}{Areal shot migration using areal data modeled from rotated SODCIG of Figure \ref{fig:datarot}.}

\newpage

\bibliographystyle{sep}
\bibliography{SEP,ref134a}
