\lefthead{Valenciano}
\righthead{Ph.D Proposal}
\nokeywords

\begin{center}
\hrule
\vspace{0.3in}
\title{RESEARCH PROPOSAL}
\hrule
\vspace{1.0in}
\title{Amplitudes and resolution in multidimensional \\ 
shot-profile imaging}
\vspace{0.5in}
\nokeywords
\email{valencia@sep.stanford.edu}
\author{Alejandro A. Valenciano }
\par
10:00am, Monday June 9, 2003
\par
\vspace{1.0in}
\includegraphics[width=0.15\textwidth]{Fig/stanfordseal.ps}
\end{center}

\newpage

\hyphenation{Cla-er-bout}

\section{Abstract}

\subsection{Problem}
The seismic data requirements for structural definition and for characterization of hydrocarbon reservoirs are different.  
Usually, reservoir characterization requires better resolution. It also needs seismic amplitudes reflecting the rock properties and the type and pressure of the fluids contained in the rock pores. 
Thus, for a reservoir geophysicist a good resolution image whose amplitudes are proportional to the reflection coefficient in the subsurface is the ideal data. This ideal image is almost never obtained by seismic data imaging techniques.

Shot-profile migration is one of the most accurate seismic imaging techniques used today in structurally complex areas. Its ability to preserve the correct amplitudes during the wavefield's extrapolation has only recently become a subject of study. Even though it is known that a deconvolution imaging condition for shot-profile migration is the way to correctly handle the amplitudes and increase the resolution, it hasn't been applied in practice.
 As the understanding of the data dimensionality and computer power increases, the use of deconvolution imaging conditions has become possible. Thus, there is an opportunity to make migrated images more suitable for reservoir characterization, in terms of resolution and amplitudes.   

\subsection{Importance}

As exploration targets move to areas with structurally complex geology (i.e., sediments under salt bodies), wave equation techniques like shot-profile migration become necessary to obtain a good image. After the reservoir is discovered, true-amplitude seismic data is integrated with other high resolution measurements (well logs) to obtain the reservoir properties. Thus,  any effort aimed at increasing the resolution and getting better amplitudes from seismic data  has a big impact in the reservoir exploitation. 
 
\subsection{Proposal}

In order to get suitable images for reservoir characterization, I propose to pose the imaging condition as a multidimensional deconvolution  in the shot-position/time domain. Posing the imaging condition in this new domain allows the use of the curvature and the time of the events to improve the image quality. 
  
I intend to study different aspects of this imaging condition. Mainly, how it affects the image amplitudes, image resolution. 
%I also going to study if the proposed imaging condition can reduce the effect on the amplitudes of the uneven illumination of the reflectors. 
I will also address image processing issues related to the implementation of the multidimensional imaging condition (like the domain where is going to be implemented, Fourier vs. space, and different schemes for the stabilization of the deconvolution). Finally, I will address  computational issues like cost and the parallel implementation, in order to make this new imaging condition  applicable to 3-D seismic data.

\newpage

\section{Introduction}
After a hydrocarbon reservoir is discovered, seismic attributes are employed to describe its characteristics. But many factors complicate the attributes estimation \cite{GEO49-10-16371648}. The lack of resolution in a finely layered medium is among them, because tuning makes attribute interpretation difficult \cite{GEO64-06-19391948, GEO64-02-04940503, GEO66-06-17141720}.   
There are also seismic processing issues that make seismic attribute analysis unreliable, especially in complex subsurface environments. 

For a long time, Kirchhoff migration methods have been used to image the earth's interior. Their ability to correctly handle amplitudes have been shown by many authors \cite{beylkin,GEO52-07-09310942}; but it also has been proven that Kirchhoff methods do not correctly handle wavefield multipathing \cite{SEG-1998-1365}, which makes them inappropriate for imaging under complex structures.

One alternative to Kirchhoff migration is shot-profile migration \cite{GEO36.03.04670481}. But it is only recently that some theoretical work has been done in order to correctly handle the amplitudes.  
\longcite{zhangyseg2002}  provide the theory for properly managing amplitudes and phases in wavefield extrapolation  during migration. But shot-profile migration consists of two parts: the first is the extrapolation of the source  and the receiver wavefields for each shot position, and the second the application of the imaging condition. Zhang's work focused on the wavefield extrapolation, assuming that the imaging condition (deconvolution of the receiver by the source wavefield) is a solved problem. 

\subsection{Previous imaging conditions}
\longcite{GEO36.03.04670481} states the imaging principle supporting shot-profile migration imaging condition and poses the imaging condition as a deconvolution of the receiver wavefield  by the shot wavefield in the time domain. Considering that there are important  factors affecting the amplitudes (in the extrapolation of the wavefields) and the instability of deconvolution, 
he suggests replacing deconvolution by the crosscorrelation of the shot and the receiver wavefields as the imaging condition. In his scheme, for each shot position, a partial image is obtained by crosscorrelating the source and the receiver wavefields along the time dimension. Then, the image is formed by stacking the partial images at each subsurface location.

 An attempt to incorporate more than one dimension in the imaging condition was done by \longcite{GEO55-09-12231234}.  He uses a 2-D deconvolution imaging condition in the $x$ position/time domain $(x,t)$ for conventional shot-profile migration. This imaging condition implies a spatially-invariant reflectivity (with $x$), which is the case of a horizontal reflector in a constant velocity medium.   In case of dipping reflectors de Bruin's method needs an independent dip estimate for a proper placement of reflectors \cite{SEG-1990-1503}. All of these makes it difficult to implement in a realistic model. 

%A different approach to migration is proposed by \longcite{GEO43-04-06910714}. It consists of synthesizing a plane wave source from  multiple shots 
% as a pre-migration processing.  \longcite{Ji.sepphd.90} uses a 1-D deconvolution imaging condition to form partial images corresponding to each plane wave. 
% To form the final image he stacks the contribution of each plane wave (form by the synthesis of multiple shots) rather than stacking the contribution of the shots as is the case of the conventional 1-D and de Bruin's imaging conditions. But the cost of plane wave synthesis is huge \cite{Audebert.sep.80.191}, since the synthesis of one plane wave requires the browsing of all the input prestack data.       
 
\subsection{My proposed imaging condition}
I propose a different imaging condition that also satisfies Claerbout's imaging principle and combines multiple shots without stacking. I first assume that the receiver wavefield is equal to the multidimensional convolution of the source wavefield with the image in the shot-position/time  domain. Then I can remove the "degradation" of the image due to the source wavefield by deconvolving the receiver wavefield by the source wavefield in the imaging condition. As an early result, I show that this multidimensional deconvolution imaging condition has the advantage of improving the image resolution while keeping the resulting image stable.   
\par
But not only the resolution can be improved. Although deconvolution as the imaging condition  adds more complexity and potential instability in the computation of the image, it also better approximates the definition of the reflection coefficient (ratio between incoming and reflected wave amplitude). Thus more "physical" amplitudes can be obtained.  

\subsection{Summary}

In the next sections, I will present the theory and practical applications of the multidimensional deconvolution imaging condition. 
%First I introduce topics like: the data and wavefield dimensionality and Claerbout's imaging principle. This makes it easier to understand the advantage of the proposed methodology.
 First, I review the work previously done: Claerbout's imaging principle and 1-D imaging conditions. This is  followed by the proposed multidimensional deconvolution imaging condition. I also show an example with synthetic data demonstrating the advantages of the  multidimensional deconvolution imaging condition over the 1-D imaging conditions. After that, I describe the future work, give a time-line, and talk about the resources I count on to finish the project. 


\section{Previous work: 1-D imaging conditions}
\subsection{Claerbout's imaging principle}

According to Claerbout's (1971) imaging principle, a reflector exists where the source and the receiver wavefields coincide in time and space.
Claerbout expresses the imaging condition as follows:
\begin{equation}
{\bf r}(x,z)=\frac{ {\bf u}(x,z,t_{d})}{ {\bf d}(x,z,t_{d})},
\label{jonip}
\end{equation}
where $x$ is the horizontal coordinate, $z$ is the depth, and $t_{d}$ is the time at which the source wavefield ${\bf d}(x,z,t_{d})$ and the receiver wavefield ${\bf u}(x,z,t_{d})$ coincide in time and space. This principle states that the reflectivity strength ${\bf r}(x,z)$  depends only on the source and the receiver wavefields at time $t_{d}$.
The time $t_{d}$ is not known {\it a priori}, therefore we need a practical way to locate the reflector position in the $(x,z)$ plane and compute its strength.

Claerbout discusses the fact that the  ${\bf r}(x,z)$ given by equation (\ref{jonip}) should be proportional to the reflection coefficient. He also points out that the definition of  equation (\ref{jonip}) is invariant to changes in source position. He dismisses this ambiguity considering that there are other  problems affecting wavefield amplitudes (attenuation and scattering in weathering layers).

%Deconvolution and crosscorrelation are operations kinematically equivalent. But  posing the imaging condition as a deconvolution has the advantage  that it remove the "degradation" of the image due to the source wavefield from the receiver wavefield. The following two subsections show 

\subsection{1-D Crosscorrelation}
A practical way to compute the reflectivity strength in equation (\ref{jonip}) is discussed in  \longcite{GEO36.03.04670481}. 
 He computes the reflector strength and position as the zero lag of the  crosscorrelation of the source and the receiver wavefields  in the time dimension (Figure~\ref{fig:match1}). 

\par

\sideplot{match1}{width=3.5in}{Source and receiver wavefields to be matched in the time dimension.}

This is expressed in the formula :
\begin{equation}
{\bf r}(x,z)= \sum_{x_s}({\bf u}\stackrel{t}{\star}{\bf d})(x,z,\tau=0,x_s),
\label{mde1}
\end{equation}
where ${\bf r}(x,z)$ is the sum over shot position $(x_s)$ of the zero lag coefficient $(\tau=0)$ of the one dimensional crosscorrelation of the receiver wavefield ${\bf u}(x,z,t,x_s)$ and the source wavefield ${\bf d}(x,z,t,x_s)$.

\par
\subsection{1-D Deconvolution}
The division in equation (\ref{jonip}) is better approximated by implementing a 1-D deconvolution imaging condition in the time dimension. It adds more complexity and potential instability in the computation of the image, but better approximates the definition of the reflection coefficient (ratio between incoming and reflected wave amplitude).

\par
This is expressed in the formula :
\begin{equation}
{\bf r}(x,z)= \sum_{x_s} \left( \frac{{\bf u}}{{\bf d}}\right)(x,z,\tau=0,x_s).
\label{mde}
\end{equation}
where ${\bf r}(x,z)$ is the sum over shot position $(x_s)$ of the zero lag coefficient $(\tau=0)$ of the deconvolution of the receiver wavefield ${\bf u}(x,z,t,x_s)$ by the source wavefield ${\bf d}(x,z,t,x_s)$.

A 1-D imaging condition implicitly makes the assumption that each shot contributes to the image with the same weight [equations (\ref{mde1}) and (\ref{mde})]. This assumption is far from reality, since  even when the subsurface geometry is not complex, reflectors are illuminated in a different way according to the source position. 

\section{Proposal: Multidimensional imaging conditions in the shot-position/time space}
I propose an imaging condition that does not assume a uniform contribution of the shots to the image. This imaging condition uses the curvature and time of the events in the  shot-position/time domain $(x_{s},t)$ to find the position and reflectivity value of the reflectors.
This can be done by posing the imaging condition as either a 2-D crosscorrelation or a 2-D deconvolution in the shot-position/time domain $(x_{s},t)$. 

Using the curvature information improves the ``discerning'' power when matching the source and receiver wavefields. Figure~\ref{fig:Wavefields_t-xs} shows the source (Figure~\ref{fig:Wavefields_t-xs}a) and the receiver (Figure~\ref{fig:Wavefields_t-xs}b) wavefields for a synthetic model of five dipping reflectors.  Notice the coincidence in position and curvature of the hyperbolas corresponding to the same reflector.   
\plot{Wavefields_t-xs}{width=6in} {Source and receiver wavefields, at a fixed $(z,x)$ in the image, to be matched in the shot-position/time plane.}

\par

\subsection{Multidimensional crosscorrelation}
The simplest way to match the source and the receiver wavefields in the shot-position/time dimensions is a multidimensional crosscorrelation.   
The zero lag coefficient of the 2-D crosscorrelation is computed as:
\begin{equation}
{\bf r}(x,z)= ({\bf u}\stackrel{t,x_s}{\star}{\bf d})(x,z,\tau=0,\xi=0).
\label{mde2}
\end{equation}
Note that whereas we sum over the shot positions ($x_s$) in equation (\ref{mde1}), we now take the zero lag coefficient $(\xi=0)$  in equation (\ref{mde2}). 

But crosscorrelation doesn't fully exploit the advantages the new imaging condition domain offers. It is easy to show that the the zero lag of the 2-D crosscorrelation in the shot-position/time dimensions is equivalent to the zero lag of the 1-D crosscorrelation in the time domain with a stack across the shot position axis.   
\subsection{Multidimensional deconvolution}
In the case of deconvolution, the relation between 2-D and 1-D (plus stacking), if it exists, is not  straightforward to illustrate. This is because deconvolution in the $(x_{s},t)$ plane should be implemented by recursive filtering. However, we know that 2-D deconvolution in the $(x_{s},t)$ plane compresses the information in both the shot and the time dimensions. This is a better way of compressing information than 1-D deconvolution followed by a stack in the shot dimension, as I will show later.

\par
Similar to equation (\ref{mde2}), the zero lag coefficient of the 2-D deconvolution is computed as:
\begin{equation}
{\bf r}(x,z)= \left( \frac{{\bf u}}{{\bf d}} \right) (x,z,\tau=0,\xi=0).
\label{mde3}
\end{equation}
Note that whereas we sum over the shot positions ($x_s$) in equation (\ref{mde}), we now take the zero lag coefficient $(\xi=0)$  in equation (\ref{mde3}). 

\subsection{How the imaging condition affects the image resolution}
I argue that, the higher the resolution in the shot-position/time domain (where the imaging condition is performed), the higher the resolution of the image. Deconvolution and crosscorrelation are operations that are kinematically equivalent. But it is advantageous to  pose the imaging condition as a deconvolution, removing the "degradation" of the image due to the source wavefield from the receiver wavefield. Figure~\ref{fig:resolution}  shows  a cartoon illustrating this idea. 

\plot{resolution}{width=4.5in}{Resolution in the imaging condition space implies resolution in the image space.}

Consider a point scatterer A in the image space (Figure~\ref{fig:resolution}a). The result of computing the 2-D crosscorrelation of the source and the receiver wavefields, at the scatterer position (A), is shown in figure~\ref{fig:resolution}b. Notice that the crosscorrelation result is centered at the zero lag ($\tau=0$ and $\xi=0$). The result of computing the 2-D deconvolution of the receiver by the source wavefield, at the scatterer position (A), is shown in figure~\ref{fig:resolution}c. Also, notice that the crosscorrelation result is centered at the zero lag. 

Now, consider a point away from the scatterer (B). The result of computing the 2-D crosscorrelation of the source and receiver wavefields, away from the scatterer position (B), is shown in figure~\ref{fig:resolution}d. Notice that the crosscorrelation center is shifted from the zero lag ($\tau=0$ and $\xi=0$), but there is some energy at the zero lag. This means that even though at this location there is no scatterer, there is going to be some energy in the image at point B. The result of computing the 2-D deconvolution of the receiver by the source wavefield, away from the scatterer position (B), is shown in figure~\ref{fig:resolution}e. Notice that no energy lies at zero lag thus no energy is going to be present in the image at point B.    

\par
%\plot{mig_tune2}{width=6.in}{Comparison of 3 different imaging conditions, right panels are close-ups of the left panels. (a) Image obtained with
% the crosscorrelation imaging condition. (b) Image obtained with the 1-D deconvolution imaging condition in the time dimension. (c) Image obtained with the 2-D deconvolution imaging condition in the $(x_{s},t)$ dimensions.}

\section{Work completed}
In this section I show the advantages of the 2-D deconvolution over crosscorrelation and 1-D deconvolution imaging conditions.
To do this, I use a synthetic  model with five dipping reflectors. Using this model, I demonstrate the advantages of 2-D deconvolution not only for the image resolution but also for estimating the correct moveout of reflectors in the angle-domain.  

\subsection{Crosscorrelation and deconvolution in the Fourier domain}

In practice, the crosscorrelation / deconvolution can be implemented as a polynomial multiplication / division in the Fourier domain. Appendix A explains this idea in detail.
The following results, are obtained  with the formulas presented in Appendix A. 
 
\subsection{Zero offset image resolution}
To compare the result of the different imaging conditions, I built a constant velocity model of five dipping layers pinching-out to the right of the model.
The deepest reflector has the steepest dip angle $(\approx 63^{\circ})$ and the shallowest has zero dip.  Figure~\ref{fig:mig_tune2} shows a comparison of the images obtained with three different imaging conditions. Figure~\ref{fig:mig_tune2}a uses the crosscorrelation imaging condition,  Figure~\ref{fig:mig_tune2}b the 1-D deconvolution imaging condition along the time dimension and Figure~\ref{fig:mig_tune2}c the 2-D deconvolution imaging condition in the $(x_{s},t)$ plane. 
 Notice the better resolution of the 2-D deconvolution image. 
\plot{mig_tune2}{width=6.in}{Comparison of 3 different imaging conditions, right panels are close-ups of the left panels. (a) Image obtained with
 the crosscorrelation imaging condition. (b) Image obtained with the 1-D deconvolution imaging condition in the time dimension. (c) Image obtained with the 2-D deconvolution imaging condition in the $(x_{s},t)$ dimensions.}
\par
The better stability of the 2-D deconvolution image  is demonstrated by comparing Figure~\ref{fig:Comp_1d_epsst} and Figure~\ref{fig:Comp_epsst}. They display the result of the 1-D deconvolution and the 2-D deconvolution for different dimensionless regularization parameter ($\lambda $) values [equations (\ref{landa1}) and (\ref{landa2})]  at a fixed x position (1.96 km in Figure~\ref{fig:mig_tune2}).
\par
Notice that when $\lambda $ decreases, the  1-D deconvolution result presents low frequency noise. The stacking across the shot position [equation (\ref{mde})] reduces the high frequency noise but might increase some low spatial frequency noise that is periodic with the shot positions. In the case of  2-D deconvolution, when $\lambda $ decreases, the noise contaminates the image in all the spatial frequency bandwidth. However, the signal to noise ratio for this imaging condition is much higher than for the 1-D deconvolution imaging condition.    

\sideplot[t]{Comp_1d_epsst}{width=3.in}{Effect of lambda on 1-D deconvolution. From bottom to top: crosscorrelation, $\lambda=0.05$, $\lambda=0.005$ and $\lambda=0.0005$} 

\par
%From these results, I conclude that the 2-D deconvolution imaging condition in the $(x_{s},t)$ plane gives a better resolution than the other imaging conditions. The 2-D deconvolution final image is less sensitive to the choice of $\lambda$ and is less affected by the low frequency noise visible in the 1-D deconvolution result. The less sensitivity to $\lambda$ choice is important, since deconvolution major handicap is the selection of the regularization parameter.

\par

\sideplot{Comp_epsst}{width=3.in}{Effect of lambda on 2-D deconvolution. From bottom to top: crosscorrelation, $\lambda=0.05$, $\lambda=0.005$ and $\lambda=0.0005$}

\subsection{CIGs and velocity analysis resolution}
Not only the zero-offset image resolution is important. Having a better image resolution for a range of offsets and angles can lead to better estimation of the reservoir properties and better velocity analysis. The angle-domain Common Image Gather (CIG) is the domain where AVA analysis is performed, because it has information about the reflectivity variation with angle \cite{shuey}.   
Angle-domain CIGs have also information about how well events are focused at depth, thus providing a natural domain for migration-focusing velocity analysis \cite{Biondi.sep.100.biondo2}.

\plot[t]{test-art-mig}{width=6in}{Comparison between the 1-D crosscorrelation and the 2-D deconvolution imaging conditions for the correct and the low velocity. (a) Angle-domain and ART CIGs computed with the 1-D crosscorrelation and the correct velocity. (b) Angle-domain and ART CIGs computed with the 2-D deconvolution and the correct velocity. (c) Angle-domain and ART CIGs computed with the  1-D crosscorrelation and the low velocity. (d) Angle-domain and ART CIGs computed with the  2-D deconvolution and the low velocity.}

I computed offset-domain CIGs with the 2-D deconvolution imaging condition using an approach analogous  to \longcite{jamespaul2000}.
After that, I transformed offset-domain CIGs  into angle-domain CIGs following the method presented by \longcite{Sava.sep.103.paul2}.  
Finally, I use the angle-domain Radon transform (ART) methodology described in  \longcite{Sava.sep.113.paul1} to estimate the curvature of the reflectors in the angle-domain CIGs. Since  the curvature of the reflectors in the ART domain is related to velocity errors, a better curvature estimation can lead to a better velocity estimation.  A flat reflector in the angle-domain CIG maps at zero q (curvature index) in the ART CIGs. Conversely, a reflector with moveout in the angle-domain CIG should map away from zero q in the ART CIG.  

\par

In Figure~\ref{fig:test-art-mig} I compare the ART CIGs with the correct velocity (Figure~\ref{fig:test-art-mig}a and Figure~\ref{fig:test-art-mig}b) and with a lower velocity (Figure~\ref{fig:test-art-mig}c and Figure~\ref{fig:test-art-mig}d). In both cases I calculated the angle-domain CIGs with 1-D crosscorrelation (Figure~\ref{fig:test-art-mig}a and Figure~\ref{fig:test-art-mig}c) and 2-D deconvolution (Figure~\ref{fig:test-art-mig}b and Figure~\ref{fig:test-art-mig}d).   
Notice that the 2-D deconvolution imaging condition gives ART CIGs  that are easier to interpret. Also notice that due to the lack of resolution, the ART CIGs computed with 1-D crosscorrelation  and the correct velocity Figure~(\ref{fig:test-art-mig}a) show events that can be misinterpreted as reflectors with the wrong velocity. In addition, the ART CIGs computed with 1-D crosscorrelation and the wrong velocity Figure~(\ref{fig:test-art-mig}c) show events that can be misinterpreted as reflectors with the correct velocity.  This problem is reduced with the 2-D deconvolution imaging condition (Figure~\ref{fig:test-art-mig}b and Figure~\ref{fig:test-art-mig}d) . 
\par

\subsection{Conclusions from work completed}
By posing the  shot-profile migration imaging condition as a 2-D deconvolution in the $(x_{s},t)$  plane the final image resolution is enhanced.
The 2-D deconvolution gives better resolution than the crosscorrelation and the 1-D deconvolution imaging conditions.
\par
The 2-D deconvolution imaging condition also gives a  more robust result than the 1-D deconvolution imaging condition because it is less dependent on the regularization parameter and it is less affected by low frequency noise. This is an important property of the 2-D deconvolution,
 since one major handicap of deconvolution application is the selection of the regularization parameter.
\par
In addition, the angle-domain CIGs computed with the 2-D deconvolution imaging condition are less affected by low resolution effects. Consequently,  we can estimate better velocities with migration-focusing velocity analysis.

\section{Work to be done}

The "work completed" shows some advantages of doing multidimensional deconvolution imaging condition, especially for image resolution. But, there are still  some questions to answer in the way this imaging condition handles:
\begin{itemize}
\item true-amplitudes,
%\item different types of noise (coherent or incoherent), 
%\item aliasing
\item and uneven illumination of the reflectors (effects of the overburden),
%\item multiple arrivals due to multipathing trough a complex geology,
%\item and missing shots (irregular data),
\end{itemize}
 among others.

There are also other problems related with the implementation of the multidimensional deconvolution that deserve some thought, for example:
\begin{itemize}
\item Fourier domain implementation 

      Considering alternative ways to avoid division by zero in the implementation of the deconvolution as a polynomial division. 

\item Space domain implementation 

      Setting the deconvolution as a least-squares inverse filtering problem makes possible to handle irregular data and also to include {\it a priori} information to regularize the inversion in the case of uneven illumination. 
\end{itemize}

The multidimensional deconvolution in the shot-position/time domain implies a change in the logistic of the shot-profile migration algorithms. That is why my research also needs to address the cost, and the implementation with 3-D data of the new imaging condition.
 
In the next subsections I present how I intend to explore the previous topics.

\subsection{True-amplitude angle dependent reflectivity}
Since the goal of doing  the new multidimensional imaging condition is to obtain images suitable for reservoir characterization, a main concern should be the adequacy of the amplitudes resulting from this process for attribute analysis. 

Although deconvolution as the imaging condition  adds more complexity and potential instability in the computation of the image, it also better approximates the definition of the reflection coefficient (ratio between incoming and reflected wave amplitude). Thus more "physical" amplitudes can be obtained.

%Other authors had obtained angle dependent reflectivity for one-way wave equation migration.  \longcite{Ji.sepphd.90} shows that using "conventional" 1-D deconvolution imaging condition in plane wave synthesis migration gives angle dependent reflectivity. \longcite{GEO55-09-12231234} use a different approach with the same objective. They perform a 2-D deconvolution in the $x$ position/time $(x,t)$ plane as the imaging condition for conventional shot-profile migration. de Bruins method has the disadvantage that in case of dipping reflectors needs an independent dip estimate for a proper placement of reflectors \cite{SEG-1990-1503}.

 I plan to test the accuracy of the 2-D deconvolution imaging condition in the $(x_{s},t)$  plane  following \longcite{Ji.sepphd.90} strategy. Ji compares four different imaging conditions (analytical, 1-D deconvolution, de Bruins, and plane wave synthesis) over a synthetic model with two  flat layers with no-lateral velocity variations. He arrives at the conclusion that plane wave synthesis has the AVA behavior  closest to the analytical result. 

By testing the 2-D deconvolution imaging condition in the $(x_{s},t)$ domain, for that simple example will allow me to isolate the image AVA behavior from other effects on amplitudes, like uneven illumination and multipathing.   
 But I also plan to introduce more complexity in the the seismic data to test how amplitudes are preserved. The next subsection address this subject.
  
\subsection{Uneven illumination of the reflectors} 

The uneven illumination of the reflectors  can lead to misinterpretation of theirs amplitudes. In the worst case, shadow zones can be seen in the seismic sections due to this effect \cite{Muerdter,Prucha.sep.97.marie1}. 
Some attempts to solve this problem have been done using migration as an imaging operator in an iterative least-squares inversion scheme \cite{nemeth99,duquet,ronen00}. 
In areas with poor illumination,
the inversion problem is ill-conditioned, therefore the inversion needs to be regularized. Since the inversion can be expensive, other authors propose to approximate the inverse of the Hessian to reduce the cost \cite{GEO66-03-09390952,SEG-2001-20812084,Guitton.sep.113.antoine3}.   

\longcite{Prucha.sep.111.marie1} showed that the use of regularized inversion with smoothing along the angle axis results in better images in shadow zones in two-way wave equation inversion. Multidimensional deconvolution can be used as an alternative to the previous approach in shot-profile migration.  The inversion problem could be solve at each image point without having to propagate the wavefields more than one. 
As I show in the next subsections the deconvolution imaging condition can be set as a least-squares inverse filtering problem [equation (\ref{spa})]. 

To regularize the inversion in the angle-domain, the inversion problem has to be set to estimate the angle dependent reflectivity. 
Analogous to \longcite{jamespaul2000}, the 2-D deconvolution imaging condition can be extended to compute a range of offsets. This is done by shifting the source and the receiver wavefields in the $x$ dimension (Appendix B). Following, an offset to offset-ray parameter slant-stack (${\bf S}$) transformation can be applied to build the angle-domain CIG \cite{Prucha.sep.100.marie1}.  

The fitting goal could be stated as:  
\begin{eqnarray}
\bf 0 &\approx & {\bf d}(h) \stackrel{t,x_s}{*} {\bf S^{'}} {\bf r}(p_h)-{\bf u}(h) \\ \nonumber
\bf 0 &\approx & \varepsilon {\bf A r}(p_h)
\end{eqnarray}
where $\stackrel{t,x_s}{*}$ states for multidimensional dimensional convolution in the shot-position dimension, $\bf A$ is the regularization operator in the angle direction, $\bf S^{'}$ is the offset-ray parameter to offset transformation operator, ${{\bf r}(p_h)}$ is the offset-ray parameter  dependent reflectivity, $h$ is the subsurface half offset, ${\bf d}(x+h,z,t,x_s)$ and ${\bf u}(x+h,z,t,x_s)$ are the shifted source and receiver wavefields. 

I intend to prove the benefit of the previous approach for uneven illumination of the reflectors using a synthetic dataset provided to SEP by Bill Symes at The Rice Inversion Project (TRIP). It consist of a flat reflector bellow a Gaussian low  velocity anomaly (Figure~\ref{fig:rays}). This velocity anomaly produces focusing defocussing effects affecting the amplitudes of the  reflector. Figure~\ref{fig:rays} also shows the rays emerging from the reflector at $x=800$ m. The reflector receive more illumination for certain angles than for others,  since the velocity anomaly produces focusing defocussing effects and the maximum recording offset is 2000 m. 

\sideplot{rays}{width=3in} {Velocity model for the TRIP synthetic. The anomaly is low velocity.} 

\subsection{Smooth polynomial division in the Fourier domain}
As we discussed before, the deconvolution can be computed as a polynomial division in the Fourier domain. My strategy to avoid the division by zero was the addition of a small regularization parameter (\ref{mde1a}). The inconvenient of this approach is that there is not a criterion for the selection of this parameter. 

Another possible approach to avoid polynomial division instability could be the smoothing the denominator and the numerator of the division. 
The zero lag coefficient of the 2-D deconvolution computed with this approach is:
\begin{equation}
{\bf r}(x,z)= \sum_{k_{x_s}}\sum_{\omega}\frac{{\bf \bar{T}}\stackrel{k_{x_s},w}{*}{\bf \bar{U}}(x,z,\omega,k_{x_s}) {\bf \bar{D}}^{*}(x,z,\omega,k_{x_s})}{{\bf \bar{T}}\stackrel{k_{x_s},w}{*}{\bf \bar{D}}(x,z,\omega,k_{x_s}) {\bf \bar{D}}^{*}(x,z,\omega,k_{x_s})},
\label{mde11}
\end{equation}
where $\stackrel{k_{x_s},w}{*}$ means multidimensional convolution in the $(\omega,k_{x_s})$ domain, ${\bf \bar{D}}^{*}(x,z,\omega,k_{x_s})$ and ${\bf \bar{U}}(x,z,\omega,k_{x_s})$ are the two dimensional Fourier transforms of the source and the receiver wavefields respectively,  and $\bf \bar{T}$ is a multidimensional smoothing filter. 

This approach presents the problem of not having a   criterion for the selection of the size of the smoothing filter. But, has the advantage of widening the wavefields in the $(\omega,k_{x_s})$ domain, what implies an increase in resolution of the deconvolution result (in the space domain). Its practical application is going to depend on the sensitivity of the deconvolution result to the filter size.  

\subsection{1-D deconvolution by the total illumination}
Another possibility for the imaging condition is to pose it as the division in Fourier domain (deconvolution in time) of two new functions. In the numerator the sum over shot-position of all the frequency components of the crosscorrelation of the source and the receiver wavefields, and in the denominator the sum  over shot-position of all the frequency components of the autocorrelation   of the source wavefield. Then the zero lag coefficient is:   
\begin{equation}
{\bf r}(x,z)= \sum_{\omega}\frac{\sum_{x_s}{\bf U}(x,z,\omega,x_s) {\bf D}^{*}(x,z,\omega,x_s)}{{\sum_{x_s}\bf D}(x,z,\omega,x_s) {\bf D}^{*}(x,z,\omega,x_s)+\varepsilon^2}.
\label{mdeato}
\end{equation}
where, ${\bf U}(x,z,\omega)$ and ${\bf D}(x,z,\omega)$ are the one-dimensional Fourier transforms of the receiver and the source wavefields respectively and $\varepsilon$ is a small number added to avoid division by zero.

This imaging condition goes beyond the normalization by the total  illumination presented by \longcite{SEG-2001-20812084}. In Rickett work the division is performed after taking the zero lag (sum over frequency). This is equivalent to do the migration (crosscorrelation) and after normalize by the illumination.     

Figure~\ref{fig:mig_tune3}  presents the result of applying this new imaging condition to our five dipping layers synthetic model (compare with Figure~\ref{fig:mig_tune2}). Notice that the results are very similar in terms of resolution to the results given by the multidimensional deconvolution (Figure~\ref{fig:mig_tune2}c).  

\plot{mig_tune3}{width=6.in}{ 1-D deconvolution by the total illumination imaging condition. The right panel is a close-up of the left panel.}

This imaging condition has the implementation advantage that it only needs to storage on disk the cumulative sums $\sum_{x_s}{\bf U}(x,z,\omega,x_s) {\bf D}^{*}(x,z,\omega,x_s)$ and $\sum_{x_s}{\bf D}(x,z,\omega,x_s) {\bf D}^{*}(x,z,\omega,x_s)$ for each frequency. It makes its practical aplication atractive.   

\subsection{Least-squares inversion imaging condition in the space domain}

In the "work completed" the deconvolution was performed as a polynomial division in the Fourier domain. This has the disadvantage of requiring the data to be defined over a regular grid in the space domain. But this condition is not likely to happen especially in a 3-D case. Furthermore, one of the advantages of shot-profile migration is that it can use sparse shots to obtain a good quality image without the need of a data space regularization. Another inconvenience of regularization is that it increases the amount of data.
     
I plan to implement the deconvolution in the data space domain \cite{Claerbout.blackwell.76,GEO29-03-03950404}. It can be stated as a least-squares inverse filtering problem \cite{GEO27-01-00040018,Claerbout.blackwell.76}:
\begin{equation}
 0 \approx {\bf d}\stackrel{t,x_s}{*}{\bf r-u} ,
\end{equation}
where $\stackrel{t,x_s}{*}$  symbol means multidimensional convolution. This approach can handle deconvolution with non-minimum phase source wavefield, but this subject has to be studied since errors can increase as the source wavefield gets far the minimum phase \cite{GEO29-01-01180120}.    

To address the  multidimensional convolution we can make use of the helix concept \cite{claerbout98}. If we put the source and the receiver wavefields in helical coordinates,  we will be able to treat the multidimensional convolution as a 1-D convolution.

Then equation (\ref{mde1a}) that is the Fourier implementation of equation (\ref{mde3})  is equivalent to the fitting goal
\begin{eqnarray}
\label{spa}
\bf 0 &\approx & \bf Dr-u \\ \nonumber
\bf 0 &\approx & \bf \varepsilon I r,
\end{eqnarray}
where $\bf I$ is the identity matrix that is used here as the regularization operator, and $ \bf D$ is a convolution matrix. 

\subsection{3-D extension} 
The addition of another dimension to the imaging condition space (from  $(x_{s},t)$ to $(x_{s},y_{s},t)$) should make the imaging condition more robust, since instead of matching hyperbolas in  $(x_{s},t)$ space, hyperboloids in  $(x_{s},y_{s},t)$ space have to be matched.   

The extension to 3-D of the multidimensional deconvolution might look trivial if we only consider the formula to handle more dimensions. But, there are important data acquisition geometries and cost related issues that need to be addressed. These are:
\begin{itemize}
\item Are all 3-D acquisition geometries suitable for this imaging condition?
 
Up to this point I have only considered the case of orthogonal acquisition geometry. In this case the deconvolution could be implemented as an extension to more dimensions of the polynomial division in the Fourier domain already presented. But the least-squares inversion implementation of the deconvolution in the data space could allow for other less regular acquisition geometries.  

\item How can we include the aperture concept in this scheme? 

Not all the shots in the survey contribute to a point in the image with the same weight. Then there is no need to store every survey shot for every image point. Thus, like in Kirchhoff migration methods, the aperture concept plays an important role in the migration cost.
\end{itemize}

\subsection{Cost reduction and implementation strategies}
In the proposed imaging condition, to build the image at each position $(x,y,z)$ we need to access the source and the receiver wavefield dimensions $(x_s,y_s,\omega)$ (in Fourier). In the usual shot-profile migration implementation, since shots are migrated one after the other, the wavefield dimensions are in this order $(x,y,\omega,z,x_s,y_s)$. Thus, to work in the shot-position/frequency dimensions, the  source and the receiver wavefield need to be saved on disk.
This is not desirable, since  the source and the receiver wavefields can be huge in a 2-D case, and impossible to handle in a 3-D case.    
But, shot-profile migration can be parallelized in different ways. Assigning the propagation of a different shot to each node of a parallel machine would allow the propagation of all the shots at the same time. Then each node could send the corresponding wavefields for a constant $(x,y,z)$ to the main node where the imaging condition is performed. This will allow access to each shot contributing to the image at a certain point without having to  store the wavefields as a whole on disk. 
  
An advantage of the multidimensional deconvolution in the shot-position/time domain is that it can be performed in a  target oriented fashion.  During a conventional shot-profile migration the wavefields can be saved in the image interval of interest (the reservoir), with a considerable reduction of the amount of data to save on disk.
 
More understanding of the requirements of communication between nodes, Input/Output, CPU nodes speed,  and  other issues is very important to implement the proposed imaging condition.
 
\section{Tools and data}

\subsection{Software}
Even though there is a shot-profile migration parallel code in the SEP, the requirements of the proposed imaging condition are going to require substantial modification of this code.
   
There is already an out-of core optimization library that can handle big volumes of data in SEP environment. It is designed to deal with 
 complex fitting goals (regularization, preconditioning etc).
 I will use this library to implement the space domain multidimensional deconvolution, allowing me to tackle the uneven illumination and the non regular shot coverage problems. 

\subsection{Hardware}
My plan involves the application of the proposed imaging condition in shot-profile migration of 3-D data. It is, therefore, very demanding as far as
computer power needs.

SEP already has two LINUX clusters with the 
speed and capacity to handle the 2D data (SEP000 and SEP100). Moreover, we have two new 
LINUX clusters that can be used to run 3-D data processes (SEP200 with 40 nodes and SEP300 with 8 nodes but a big RAM memory). Furthermore, SEP continuously increases it's computer power, thus by the time I will start to work with 3-D  data I expect it to  
 be sufficient.  

\subsection{Field Data}
There is a data library at SEP that will allow me to prove some of the concepts. We have in the library a 3-D dataset acquired by Elf in the North sea (figure~\ref{fig:L7d-ComAz-mig-all}) that presents some important features I want to study. For example sediments pinching out against a salt body. The structure of the salt dome is highly 3-D, with serious illumination problems under its edges. The horizontal section at depth 900 m highlights complex patterns that had been interpreted as turbidite channels. Those turbidite channels are ideal to test the improvements of resolution of the multidimensional imaging condition. 

This is also a well known data set, with a well  known velocity model that has been extensively  used by SEP students.
The velocity model was created by Elf \& IFP (Project ``Helios'') using the so-called ``SMART'' method (Sequential Migration Aided Reflection Tomography). 
 \plot[t]{L7d-ComAz-mig-all}{width=5.in}{Elf North sea 3-D dataset. Vertical section shows the sediments pinching out against the salt and the horizontal section shows the turbidite channels. }

To test all the possible benefits of the proposed imaging condition I need at least one well with enough log measurements  to tie the migration results with. SEP is expecting to receive the MC311 Gulf
of Mexico dataset from ExxonMobil. MC311 dataset includes: pre-stack surface seismic data, velocity model, well log data, check-shot data and
offset vsp data. I expect these dataset to be perfect to test my results. 
 
\section{Time-Line}
The time line for this project is:

\begin{itemize}

\item July/2003 - September/2003: Changes in shot-profile migration strategy 

\item October/2003 - December/2003: True-amplitude angle dependent reflectivity

\item January/2004 - June/2004 : Least-squares inversion imaging condition implementation  

\item July/2004 - December/2004: Uneven illumination of the reflectors 

\item January/2005 - December/2005: 3D extension and 3-D real data.

\item January/2006 - June/2006: Write thesis.

\item July/2006: Graduate

\end{itemize}

\bibliographystyle{sep}
\bibliography{SEP,SEG,MISC,GEOTLE,bei}

\APPENDIX{A}
\subsection{Crosscorrelation and deconvolution in the Fourier domain}

The zero lag of the 1-D crosscorrelation can be expressed as:
\begin{equation}
{\bf r}(x,z)= \sum_{x_s}\sum_{\omega}{\bf U}(x,z,\omega,x_s) {\bf D}^{*}(x,z,\omega,x_s),
\label{md4a}
\end{equation}
where ${\bf r}(x,z)$ is the zero lag coefficient of the crosscorrelation, which is computed by summation over the frequencies and ${\bf U}(x,z,\omega)$ and ${\bf D}(x,z,\omega)$ are the one-dimensional Fourier Transforms of the receiver and the source wavefields, respectively. The contribution of each shot (located at $x_s$) is added to form the final image.


The 1-D deconvolution can be computed in the Fourier domain as a polynomial division. The zero lag coefficient is computed as the sum over frequencies: 
\begin{equation}
{\bf r}(x,z)= \sum_{x_s}\sum_{\omega}\frac{{\bf U}(x,z,\omega,x_s) {\bf D}^{*}(x,z,\omega,x_s)}{{\bf D}(x,z,\omega,x_s) {\bf D}^{*}(x,z,\omega,x_s)+\varepsilon^2(x,z,x_s)}.
\label{mdea}
\end{equation}
where, ${\bf U}(x,z,\omega)$ and ${\bf D}(x,z,\omega)$ are the one-dimensional Fourier transforms of the receiver and the source wavefields respectively and $\varepsilon$ is a small number added to avoid division by zero. The contribution of each shot (located at $x_s$) is added to form the final image.

\par

\longcite{Jacobs.sepphd.34} discuss the difficulties of choosing a spatially variable $\varepsilon$. In practice, defining $\varepsilon$ as a function of a dimensionless parameter $\lambda$ makes its selection easier \cite{Claerbout.blackwell.92}. 
I define this dependence as  
\begin{equation}
\varepsilon^2(x,z,x_s)=\lambda < {\bf D}(x,z,\omega,x_s) {\bf D}^{*}(x,z,\omega,x_s) >, 
\label{landa1}
\end{equation}
where $<>$ indicates the mean in the frequency dimension. $\lambda$ can be set constant for the whole image because it is independent of the data scale. 

The 2-D deconvolution can be computed as a polynomial division in the Fourier domain. But in this case, the wavefields that are divided, ${\bf \bar{U}}(x,z,\omega,k_{x_s})$ and ${\bf \bar{D}}(x,z,\omega,k_{x_s})$, are the two-dimensional Fourier transforms of the receiver and the source wavefields, respectively. 
\par
The zero lag coefficient of the 2-D deconvolution is computed as:
\begin{equation}
{\bf r}(x,z)= \sum_{k_{x_s}}\sum_{\omega}\frac{{\bf \bar{U}}(x,z,\omega,k_{x_s}) {\bf \bar{D}}^{*}(x,z,\omega,k_{x_s})}{{\bf \bar{D}}(x,z,\omega,k_{x_s}) {\bf \bar{D}}^{*}(x,z,\omega,k_{x_s})+\varepsilon^2(x,z)}.
\label{mde1a}
\end{equation}
Note that whereas we sum over the shot positions ($x_s$) in equation (\ref{mdea}), we now sum over the shot-position spatial frequency $(k_{x_s})$ in equation (\ref{mde1a}). 

\par

Also notice that the regularization parameter $\varepsilon(x,z)$ is spatially variable but constant in the $(\omega,k_{x_s})$ plane. It is defined as
\begin{equation}
\varepsilon^2(x,z)=\lambda < {\bf \bar{D}}(x,z,\omega,k_{x_s}) {\bf \bar{D}}^{*}(x,z,\omega,k_{x_s}) >
\label{landa2}
\end{equation}
where $<>$ is the mean in the $(\omega,k_{x_s})$ plane.

%\APPENDIX{B}



\APPENDIX{B}
The nonzero-offset reflectivity can be computed as follows:
 \begin{equation}
{\bf r}(x,z,h)= \sum_{k_{x_s}}\sum_{\omega}\frac{{\bf \bar{U}}(x-h,z,\omega,k_{x_s}) {\bf \bar{D}}^{*}(x+h,z,\omega,k_{x_s})}{{\bf \bar{D}}(x+h,z,\omega,k_{x_s}) {\bf \bar{D}}^{*}(x+h,z,\omega,k_{x_s})+\varepsilon^2(x+h,z)}.
\label{mde1h}
\end{equation}
where $h$ is the subsurface offset.
% 
%\section{Data and wavefield dimensionality}
%To implement a better imaging condition that is also feasible in practice, it is important to understand 3-D prestack data and wavefield dimensionality. 
%3-D prestack seismic data are defined in a 5-D continuum $(t,x_s,y_s,x_g,y_g)$ \cite{Biondi.sep.98.biondo}, where $t$ is time, 
%$x_s,y_s$ and $x_g,y_g$  are the surface coordinates of the sources and the receivers, respectively.
%
%\par 
%
%After applying the first step of shot-profile migration (source and receiver wavefield construction) the dimensionality of the data increases. Thus, for each shot position $(x_s,y_s)$ the source and the receiver wavefields ${\bf u}(x,y,z,t)$ and  ${\bf d}(x,y,z,t)$ are defined in four dimensions (Figure~\ref{fig:match3}), where $(x,y,z)$ are the image space dimensions.
%
%\par
%
%\sideplot{match3}{width=3.5in} {Data dimensionality. The source and receiver wavefields are 6-D datasets. }
%
%In the following analysis, I explain how to combine the source and the receiver wavefields to obtain an image. For simplicity, I restrict my analysis to 2-D prestack data. In this case, the source and the receiver wavefields are 4-D datasets defined in $(x,z,t)$ for each shot position $x_s$.
%      
 
