\lefthead{Guerra}
\righthead{Research Proposal}
\email{guerrac@stanford.edu}
\footer{PhD Proposal}
\title{3D Prestack Exploding-Reflector Modeling for Migration Velocity Analysis}
\author{Claudio Guerra}

\section{ABSTRACT}
 
\subsection{Problem} 
In 
Velocity analysis for depth migration is an iterative process with basically three steps: 1) prestack depth migration with a presumed velocity model; 2) measurement of misfocusing of image gathers; and 3) velocity model update based on the misfocusing property from step 2. Depending on the migration algorithm used, the first step can be, by far, the most time-consuming of all. For velocity analysis in 3D depth migration using wave-equation algorithms to be cost effective, it is therefore necessary to reduce the turnaround time of the step 1. I propose a procedure to expand to 3D the method recently proposed by \longcite{Biondi.sep.124.biondo1,Biondi.sep.131.biondo1}, which aims to decrease the amount of data to be migrated, but maintain the kinematic information necessary to perform velocity update. 

\subsection{Importance}
In the past decade, oil exploration has been being conducted in increasingly more geologically complex areas. The imaging challenges along with the continuous growth of computer power has augmented the relative importance of depth migration, in particular wave-equation migration. A variety of wave-equation migration algorithms have been developed to balance cost and final image quality. Velocity model definition for depth migration, however, remains strongly dependent on Kirchhoff algorithms because of their flexibility and low cost when compared with wave-equation migration velocity analysis. In areas where the high-frequency approximation is violated and the wavefields become too complex, this dependence might result in inaccurate velocity models that yield poor quality or unreliable images even when a state-of-the-art wave-equation algorithm is used to produce the final image. In these areas it is therefore desirable to use wave-equation migration to generate image gathers in order to evaluate the accuracy of the velocity model.

\subsection {Proposal}
I propose to extend to three dimensions the originally proposed 2D prestack modeling for migration velocity analysis \cite{Biondi.sep.124.biondo1,Biondi.sep.131.biondo1}. 
\par
In these papers, Biondi showed that many common image gathers (CIGs) can be combined together prior to modeling to simulate one single modeling experiment. Unfortunately, the gain in effectiveness is followed by degradation of the image obtained after migration of the modeled data. The artifacts are caused by cross-talk of events when CIGs are combined together. I will propose a way to mitigate this problem of cross-talk. Random phase-encoding of shot gathers \cite{Romero} is a viable strategy to use during the migration step. I will also investigate an optimal way to combine CIGs before the modeling to minimize the cross-talk. 
\par
To verify the validity of the modeling approach, after migration of the synthesized dataset I will apply the complete procedure to perform velocity update by estimating residual moveout and 3D-reflection tomography. This should enable a complete cost effective procedure that does velocity model estimation for 3D prestack depth migration using one-way shot profile migration.
\par
After establishing the procedure, I will apply it to a real 3D dataset.

\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 to build 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 giving rise to complex wavefields, inaccurate updates of the velocity model and even divergence of the solution can occur. To circumvent this problem, geological constraints are frequently applied in the inversion of the misfocusing parameters into velocity variations or even directly on the velocity model. A way to overcome the complex overburden problem is to use migration methods that better handle velocity model complexity. Full one- or two-way wave equation migrations are expected to overcome this problem, but at expense of a higher cost. Therefore, to make wavefield-continuation methods practically applicable it is necessary to reduce the computational effort. 
\par
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
\longcite{Biondi.sep.124.biondo1,Biondi.sep.131.biondo1} proposed a method to model prestack wavefields for use in velocity updating. The method can be considered as a prestack extension of the exploding reflector concept \cite{loewenthal}. Initially, consider common image gathers (CIG) computed with an assumed velocity model. Ignoring illumination problems and assuming that the migration was performed with an inaccurate velocity model, some misfocusing will occur, showing up as reflector amplitudes positioned away from the zero subsurface offset. After a preprocessing step that includes the alignment of the prestack CIGs along the normal to the reflector instead of along the vertical; that is, by interpreting the offset as offset along the geological-dip instead of along the horizontal, we then inject the image as the initial condition to synthesize source and receiver wavefields at the surface, where they are recorded by an areal array that can span the whole survey. Many CIGs can be combined together to speed up the modeling process, but at the expense of resulting cross-talk of events. The 3D extension of this method has the potential to dramatically decrease the cost of 3D velocity model building using shot-profile migration because of the areal shot configuration used to record the synthesized source and receiver wavefields.



\section{PROPOSED SCHEME}
I plan to use one-way wave equation migration scheme to generate CIGs and, again, the one-way wave equation to propagate to the surface the source and receiver wavefields. These algorithms, coded by Brad Artman, are available at SEP and have been already used in two PhD thesis to migrate 3D data \cite{Artman,Alvarez130} as well as numerous times to migrate 2D data. This algorithm is therefore considered quite stable and reliable. The migration code allows use of a generalized source function, so the source wavefield resulting from the modeling step can be easily incorporated in the migration. A more concise version of this algorithm is also available but it is still being evaluated. I will code the 3D preprocessing step required to align the CIGs along the normal to the reflector. 
\par
As mentioned earlier, many CIGs can be combined together and used in a single modeling step to make the whole modeling procedure more cost effective. Depending on how close are the CIGs, however, a severe cross-talk between reflectors can occur, degrading the prestack image after migration. To mitigate the cross-talk I will apply random phase-encoding during the migration step and investigate the effect of randomly select CIGs to be used as initial condition in the modeling step.
\par
The 3D extension of the prestack modeling for migration velocity analysis relies on the co-planarity condition \cite{biondi2}, which states that the source and receiver rays must lie on the same plane, with all the possible propagation planes sharing the line that connects the source and receiver location at each depth level. The derivation of the 3D formulas is work still to be done.

\section{WORK COMPLETED}
To the present, I have been working on implementing the preprocessing step in 2D and investigating the synthesis and migration of areal shots.

\subsection {Rotation of Common Image Gathers}
 In the preprocessing step a rotation is applied to the CIGs that changes 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, one ensures that the apparent source and receiver ray local directions are respected when injecting the rotated CIGs as initial conditions for the modeling of source and receiver wavefields. The rotation can be implemented in the wavenumber domain ($k_z-k_x-k_h$). 
\par
As seen in Figure \ref{fig:sodcig}, 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:sodcig} 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 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}

\plot{sodcig}{width=5.00in,height=7.0833in}{Geometric representation of the rotation of an image event from horizontal subsurface offset to geologic-dip subsurface offset (from \longcite{Biondi.sep.131.biondo1}).}

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. Therefore, after rotation, two new set of 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. These initial conditions are expressed by the following equations:}
\begin{equation}
P_s(t=0,x=x{_{_\xi}}-\hat{h}{_{_\xi}}^s,z=z{_{_\xi}})=\hat{I}_s(z{_{_\xi}},x{_{_\xi}},\hat{h}{_{_\xi}}^s)
\label{eq:ic_sou}
\end{equation}
\begin{equation}
P_g(t=0,x=x{_{_\xi}}+\hat{h}{_{_\xi}}^g,z=z{_{_\xi}})=\hat{I}_g(z{_{_\xi}},x{_{_\xi}},\hat{h}{_{_\xi}}^g)
\label{eq:ic_rec}
\end{equation}
Note that because of the negative sign for subsurface offset in equation \ref{eq:ic_sou}, the new image gather for the source case is flipped along the offset direction prior to injection in the modeling. Conventional zero-offset exploding-reflector modeling will take place if reflectors are perfectly focused at zero subsurface offset.

\par
Here I show results of the rotation scheme on two different data. I computed them analytically to easily evaluate the action and correctness of the rotation. One emulates a horizontal geology illuminated by $45^{\circ}$-reflection angle; the other can be described as a zero-reflection angle experiment illuminating a $45^{\circ}$-dip reflector. For the horizontal case illuminated by $45^{\circ}$-aperture angle, after rotation the reflector is expected to have the same orientation after rotation in both new image gathers because, as mentioned before, the image gather for the source wavefield has the offset direction reversed. Therefore, the initial conditions show opposite direction. The same reasoning applied to the $45^{\circ}$-dip reflector illuminated by zero reflection angle shows that, after rotation, the new image gathers must exhibit opposite dip directions. 
\par
Figures \ref{fig:geo0} and \ref{fig:geo45} pertain to the horizontal and $45^{\circ}$-dipping reflector, respectively. The orientation of the reflectors are in accordance to what is expected. As can be seen in Figure \ref{fig:geo45t}, however, for the $45^{\circ}$-dipping reflector the images show an unexpected behavior, becoming asymmetric related to the zero offset as we shift to different CMP positions and being even wrapped-around. I included in both Figures \ref{fig:geo45} and \ref{fig:geo45t} the padded traces to illustrate the `wrap-around'-like problem. I don't yet understand this behavior; it still must be corrected.

\subsection {Synthesis and Migration of Areal Shots}
As mentioned earlier, 3D depth migration by means of wavefield extrapolation, specially using the shot-profile scheme, is still considered too expensive to be applied routinely in complex areas. There exist ways to decrease the cost of 3D shot-profile migration by sampling and/or combining shot records in such a way that the number of sampled/combined shot records are considerably smaller than the original number of shot profiles. Plane-wave migration (Liu***) and phase-encoding migration (Romero***) are the most widely used methods to achieve this goal. \par
Here, I present results of plane-wave migration of areal shots synthesized according to the concept of controlled illumination introduced by (Rietveld***). 
examples of 



\section{WORK TO BE DONE}
The first task still to be solved is to get correct rotation results for various simple models. After that, I will increase the model complexity to evaluate the validity of the method. 
\par
The extension to the 3D situation will be done afterwards, along with the investigation of a method to select CIGs to be optimally combined together prior to modeling so as to minimize the cross-talk and implementing the phase-encoding routine to the migration.
\par
I will use 3D-reflection tomography to evaluate if the kinematic information is sufficiently reliable to perform velocity updates. Finally, I will apply the methodology to a 3D marine real data. 

\section{TIME-LINE}
\begin{enumerate}
\item \underline{Winter 2007-2008}\\
Correct the 2D rotation results\\
Test the rotation for complex models\\
Course requirements
\item \underline{Spring 2007-2008}\\
Investigate an optimal selection of CIGs\\
Implement phase-encoding routine\\
Course requirements
\item \underline{Summer 2007-2008}\\
Extend to 3D rotation and prestack modeling\\
Test with SEG-EAGE salt model\\
Course requirements
\item \underline{Autumn 2008-2009}\\
3D Migration of real data\\
Application of 3D-reflection tomography\\
Start writing thesis
\item \underline{Winter 2008-2009}\\
Finalize deliverables for graduation\\
\end{enumerate}

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

\plot{geo0}{width=6.00in,height=4.00in}{Image gathers after rotation to be used as the initial condition to model the source wavefield (a) and the receiver wavefield (b) for the case of horizontal geology illuminated by $45^{\circ}$-aperture angle.}

\plot{geo45}{width=6.00in,height=4.00in}{Image gathers after rotation to be used as initial condition to model the source wavefield (a) and the receiver wavefield (b) for the case of $45^{\circ}$-dipping geology illuminated by $0^{\circ}$ aperture angle.}

\plot{geo45t}{width=6.00in,height=4.00in}{Image gathers after rotation to be used as initial condition to model the source wavefield (a) and the receiver wavefield (b) for the case of $45^{\circ}$-dipping geology illuminated by $0^{\circ}$ aperture angle. The reflector was supposed to be centered around zero-offset and 900m depth.}
