\documentclass{segabs}
\newcommand{\rs}[1]{\mathstrut\mbox{\scriptsize\rm #1}}
\newcommand{\rr}[1]{\mbox{\rm #1}}  
\renewcommand{\figdir}{./Fig} % figure directory
\usepackage{multicol} 
\usepackage{natbib} 
\usepackage{alltt}
%\input{preamble} 

\begin{document}
\title{Fast 3D velocity updates using the pre-stack exploding reflector model} 
%:
\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\author{Claudio Guerra and Biondo Biondi \footnotemark[1], Stanford University.}
%: 
\footer{Fast 3D velocity updates using PERM wavefields}
\lefthead{Guerra \& Biondi}   
\righthead{Fast 3D velocity updates using PERM wavefields}
\maketitle

\begin{abstract}
We show 3D real data results of migration velocity analysis by wavefield extrapolation using data synthesized by the pre-stack exploding-reflector model (PERM). PERM is a generalization of the exploding reflector model (ERM) in a sense that migration of PERM data can generate a pre-stack image, which is not achievable with ERM. PERM data allow fast migration velocity analysis, since its size can be orders of magnitude smaller than, for instance, plane-wave data. Also, because the modeling can be limited to the region where the velocity is to be updated, additional savings are possible by solving for the velocity in a target-oriented way. 
\end{abstract} 
%
\section{Introduction}
%
In areas of complex geology, pre-stack depth migration is required not only for imaging purposes but also for velocity estimation. In such areas, migration by wavefield extrapolation has been widely used to produce the final image because it properly handles complex distortions of the wavefields. Migration velocity analysis by wavefield extrapolation \cite[]{sava2004a,shen:VE49} promises to produce more reliable results in those areas. However, due to the high computational cost, wavefield extrapolation methods are rarely used to estimate the migration velocity model in 3D projects \cite[]{fei:2283}, where ray-based methods are the industry standards. In addition to the lower cost, ray-based methods are very flexible with respect to strategies for defining the velocity model \cite[]{stork:680,kosloff:1815,billette:1889}. But despite their advantages, ray methods do not satisfactorily describe complex wave propagation in the presence of large lateral velocity contrasts. In this case, a more complete description of the wavefield complexities is needed, and therefore we face the challenge of decreasing the cost of migration velocity analysis by wavefield extrapolation while maintaining its robustness. 
\par
A typical way of decreasing the cost of wavefield extrapolation is to reduce the amount of input data by linearly combining the wavefields using plane-wave phase encoding \cite[]{whitmore,liu01} and random phase encoding \cite[]{romero,Sun01}. Combination of wavefields is exploited by the prestack-exploding reflector model (PERM) \cite[]{biondi:3056} to significantly reduce the data size while keeping the migrated image crosstalk-free. This is achieved by selecting subsurface-offset common-image gathers (SODCIGs) separated by a decorrelation distance such that the wavefields from different SODCIGs in the same modeling experiment do not correlate during migration. \cite{guerra:3964} show that further reduction can be achieved by randomly phase encoding the modeling experiments, significantly decreasing the cost of migration velocity analysis iterations.  
\par   
An interesting feature of PERM is that, because the wavefields are upward propagated, data can be collected at any depth level. In migration velocity analysis, PERM wavefields can be propagated only in the region with velocity inaccuracies, such that velocity update can be performed in a target-oriented way, which contributes to an additional cost reduction of migration velocity analysis. 
\par
In 3D, reduction of the data size can be considerable if the initial image used to model PERM wavefields has only in-line subsurface offsets, as in the common-azimuth approximation. We show that 3D-PERM data size can be potentially up to two orders of magnitude smaller than 3D-plane wave data. The usefulness of PERM data to migration-velocity analysis with wavefield extrapolation is illustrated using a North Sea 3D dataset. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{pre-stack-exploding-reflector model} 
The fundamental idea of PERM is to model data that describes the correct kinematics of an isolated SODCIG. When using Born modeling, since we do not know beforehand which shots contribute to forming the image at a point in the model space, we have to model several shots. Ideally, instead of performing many modeling experiments, we would like to synthesize a small amount of data with the condition that migration has the same kinematics as the initial SODCIG. This can be achieved by extrapolating source and receiver wavefields stating from a pre-stack image computed with wave-extrapolation methods.
\par
The modeling of PERM source $D_{\scriptscriptstyle P}$ and receiver $U_{\scriptscriptstyle P}$ wavefields can be carried out by any wavefield-continuation scheme. Here, we use the following one-way wave equations:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 s_0^2({\bf x})-|{\bf k}| ^2}\right) D_{\scriptscriptstyle P}({\bf x},\omega;{\bf x_m}) = I_D({\bf x_m},{\bf h}) \\
D_{\scriptscriptstyle P}(x,y,z=z_{\rm max},\omega;{\bf x_m}) = 0 \end{array} \right.,
\label{eq:perm_bD}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s_0^2({\bf x})-|{\bf k}| ^2}\right) U_{\scriptscriptstyle P}({\bf x},\omega;{\bf x_m}) = I_U({\bf x_m},{\bf h}) \\
U_{\scriptscriptstyle P}(x,y,z=z_{\rm max},\omega;{\bf x_m}) = 0 \end{array} \right.,
\label{eq:perm_bU}
\end{eqnarray}
where $I_D({\bf x_m},{\bf h})$ and $I_U({\bf x_m},{\bf h})$ is the isolated SODCIG at the horizontal location ${\bf x_m}$ for a single reflector, suitable for the initial conditions for the source and receiver wavefields, respectively. The subsurface-offset $\bf h$ is be parameterized as ${\bf h}=(h_x,h_y)$. The initial conditions are obtained by rotating the original unfocused SODCIGs according to the apparent geological dip of the reflector \cite[]{biondi07}. By doing so, image-point dispersal is corrected such that the velocity information needed for migration velocity analysis is maintained. 
\par
To decrease the number of modeling experiments, linearity of wave propagation can be used to combine isolated SODCIGs and inject them simultaneously into one single model experiment, using the same modeling equations as above with the initial conditions replaced by the combined SODCIGs. The selection of SODCIGs can be thought of as the multiplication of the pre-stack image by 2D $Comb$ functions, which are shifted laterally to select new set of SODCIGs to initiate the modeling of another pair of combined wavefields. After shifting along one period of the sampling function in the $x$ and $y$ directions, all the points on the reflector are used in the modeling. Consequently, the number of modeling experiments equals the number of lateral shifts of the sampling function.
\par
The choice of the sampling period determines the amount of crosstalk in the migrated image. To obtain a crosstalk-free image, the sampling period must be large enough that wavefields initiated at different SODCIGs do not correlate. PERM wavefields generated from SODCIGs at an interval equals to twice the subsurface-offset range still contribute to the image at the central SODCIG. For the same reason, no crosstalk is generated during migration if the period of the sampling function is larger than that interval. Since the focusing of energy in the SODCIG is velocity-error dependent, so is the the period of the sampling function. Therefore, for small velocity errors a small sampling period can be used and, consequently, a smaller number of combined modeling experiments is needed.
\par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection*{PERM using a common-azimuth migrated image}
In the way PERM is formulated, there is no restriction on the dimensionality of the pre-stack image used as the initial condition for the modeling, which means that if the original data have sufficient cross-line offsets as in the acquisition geometries with wide range of azimuths \cite[]{regone:SM231,kapoor:926}, the initial conditions are a five-dimensional hypercube on $\bf x$, $h_x$ and $h_y$. To synthesize PERM data starting with the five-dimensional initial conditions such that no crosstalk is generated during migration, the minimum number of modeling experiments is $2n_{h_{x}}n_{h_{y}}$, where $n_{h_{x}}$ and $n_{h_{y}}$ are the number of subsurface offsets in the $x$ and $y$ directions. Considering usual parameters, the number of modeling experiments may be as low as several hundreds. This data reduction is very substantial if we compare, for instance, with data reduction achieved by 3D-plane-wave decomposition. Using plane waves, to obtain artifact-free SODCIGs due to the lack of illumination from some propagation directions we need to migrate few thousands of plane waves. This means that 3D-PERM data size can be one order of magnitude smaller than 3D-plane wave data.
\par
Despite the recent good migration results obtained in geologically complex areas using wide-azimuth data, narrow-azimuth acquisition is still the industry standard. Narrow-azimuth data can be efficiently imaged by common-azimuth wave-equation migration (CAM) \cite[]{biondi:1822}. By assuming zero cross-line offset in contrast with the full-azimuth migration, instead of a five-dimensional hypercube, CAM images are four-dimensional hypercubes in $\bf x$ and $h_x$. Because of the lower dimensionality of CAM images, when using them as the initial conditions to synthesize PERM data, the SODCIGs in the cross-line direction can be sampled continuously, as depicted in Figure \ref{fig:cam01}b. Contrast this case with the five-dimensional initial conditions for the full azimuth case of Figure \ref{fig:cam01}a. The continuous sampling of SODCIGs in the cross-line direction yields one more order of magnitude of data reduction. Therefore, under the common-azimuth assumption, 3D-PERM data size can be two orders of magnitude smaller than 3D-plane wave data.
\plot{cam01}{width=1.1\columnwidth}{The initial conditions for synthesizing PERM data from CAM images can be specified as in b) because no pre-stack information exists in the cross-line direction, in contrast with the full azimuth situation in a).}
\par
To illustrate the validity of the above assumptions, a split-spread data with maximum offset of 1587.5 m was computed using 3D-Born modeling on a 30$^{\circ}$ dipping reflector with 45$^{\circ}$ azimuth with respect to the acquisition direction. There are 96 in-lines and cross-lines spaced 25 m apart. The offset interval is 25 m.  The velocity used in the modeling is the 1D function $v(z) = (1500 + 0.5z)$ m/s. The Born data was input to CAM with a 5$\%$ slower velocity. Migration results can be seen in Figures \ref{fig:cam02}a and \ref{fig:cam02}b for SODCIGs positioned at $(x = 750$ m, $y = 600$ m$)$ and  $(x = 750$ m, $y = 1000$ m$)$, respectively. The panel on the left is the SODCIG, which contains 21 subsurface offsets ranging from $-250$ m to $250$ m. The panel in the middle is the in-line at zero subsurface offset, with $y = 600$ m (Figure \ref{fig:cam02}a) and $y = 1000$ m (Figure \ref{fig:cam02}b). The panel on the right is the cross-line  at zero subsurface offset, with $x = 750$ m.
\par
In the common-azimuth regime, the computation of the dip-independent initial conditions is performed by simply rotating the SODCIGs in the in-line direction, since no cross-line offset is computed in migration. %The underlying assumption that source and receiver rays cross, which makes the rotation valid in 3D, implies that they are coplanar. This assumption is in agreement with the co-planarity condition present in the formulation of CAM.
\par
Source and receiver wavefields were modeled using continuous sampling along the cross-line direction and sampling period of 48 in the in-line direction. This period is sufficient to avoid crosstalk during the areal-shot migration, given that the number of subsurface-offsets of the pre-stack image is 21. The synthesized 3D receiver wavefield is shown in Figure \ref{fig:cam03}. The left panel is the in-line at $y = 1200$ m, the right panel is the cross-line at $x = 1400$ m, and the top panel is the time-slice at $t = 0.5$ s.
\par
The 3D migration of the 48 areal shots with the velocity underestimated by 5$\%$ is shown in Figures \ref{fig:cam04}a and and \ref{fig:cam04}b for SODCIGs positioned at $(x = 750$ m, $y = 600$ m$)$ and  $(x = 750$ m, $y = 1000$ m$)$, respectively. The kinematics of the SODCIGs computed with PERM wavefields matches those of the SODCIGs computed with CAM, especially between $-200$ m and $200$ m. This enables the use of 3D PERM wavefields computed from CAM images in optimization of migration velocity, as will be shown next for a 3D survey from the North-Sea.
\plot{cam02}{width=.8\columnwidth}{CAM of the 3D-Born data. Middle: in-line at zero-subsurface offset, and $y = 600$ m (Figure \ref{fig:cam02}a) and $y = 1000$ m (Figure \ref{fig:cam02}b). Right: cross-line at zero-subsurface offset, and $x = 750$ m. Left: SODCIGs.}
\plot{cam03}{width=.8\columnwidth}{3D-PERM receiver wavefield. The left panel is the in-line at $y = 1200$ m, the right panel is the cross-line at $x = 1400$ m, and the top panel is the time-slice at $t = 0.5$ s.}
\plot{cam04}{width=.8\columnwidth}{3D-areal-shot migration of PERM data. Middle: in-line at zero-subsurface offset, and $y = 600$ m (Figure \ref{fig:cam04}a) and $y = 1000$ m (Figure \ref{fig:cam04}b). Right: cross-line at zero-subsurface offset, and $x = 750$ m. Left: SODCIGs.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
\section*{Example -- 3D velocity update}
Migration velocity analysis by wavefield extrapolation using 3D PERM data was performed on a 3D real dataset from the North Sea with 3600 m maximum offset. This dataset was submitted to azimuth-moveout and imaged with CAM with an initial velocity sufficiently accurate up to the top of the chalk layer. From this reflector to deeper levels velocity is underestimated (Figure \ref{fig:nsea03}) as revealed by reflectors curving up in angle-domain common-image gathers at the bottom of Figure \ref{fig:nsea01}. 
\par
We computed velocity updates for the chalk layer in part of the 3D volume, which has 12 km in the in-line direction and 4 km in the cross-line direction. The velocity model definition for the entire 3D cube is ongoing work. We start by interpreting and windowing the base of chalk, indicated by $\bf BC$ in Figure \ref{fig:nsea01}, in the pre-stack image. The top of chalk, indicated by $\bf TC$ in Figure \ref{fig:nsea01}, is used as the upper limit for the velocity updates. The windowed reflector is then submitted to the rotation according to the apparent geological dip to correct for the image-point dispersal and it is used as the initial conditions to model 45 pairs of PERM source and receiver wavefields. Figure \ref{fig:nsea00} shows the PERM receiver wavefield of one pair of wavefields. The wavefields are initiated by SODCIGs 1350 m apart. Since the offset range computed in the velocity optimization is 600 m, no crosstalk is generated.  The wavefields are collected at 600 m depth, which is close to the top of the target-zone. Therefore, during the migration velocity analysis the wavefields are propagated between this depth level and the deeper level, which is 2800 m.
\par 
To optimize the velocity, we use a nonlinear conjugate gradient algorithm. The objective function we minimize computed via differential-semblance optimization (DSO) \cite[]{symes:654,shen:VE49}, which corresponds to weighting the pre-stack image computed using the current velocity with the absolute value of the subsurface offset. The updated velocity is constrained to be within bounds 50$\%$-lower and 50$\%$-higher the initial velocity. Because of unbalanced amplitudes in the original data, is necessary to strongly smooth the gradient. It is B-spline smoothed, with node $xy$-spacing of 1000 m. The optimization stopped after 5 iterations. The value of the objective function dropped 50$\%$ with respect to the initial value. 
\par
The optimized velocity model is shown in Figure \ref{fig:nsea04} for the in-line at 3760 m and cross-line at 3500 m. Notice that, as expected, the optimized velocity is higher than the initial velocity. CAM using the optimized velocity confirms the correctness of the velocity update, by showing more focused reflectors and flatter ADCIGs than that obtained with the initial velocity (Figure \ref{fig:nsea02}). Compare with Figure \ref{fig:nsea01}. Some residual moveout is still present as can be seen in the rightmost ADCIGs of Figure \ref{fig:nsea02}. In this region, the amplitudes in the initial pre-stack image are weaker, so are the velocity updates. We expect to obtain more even velocity updates by pre-conditioning the amplitudes in the original data. This is still under investigation.
\plot{nsea03}{width=1\columnwidth}{Initial velocity model. Left: cross-line at 3500 m. Right: in-line at 3760 m.}
\plot{nsea01}{width=1\columnwidth}{Angle-gathers after CAM with the initial velocity. Top: zero-angle section. Bottom: ACDCIGs selected roughly at the same x position in the section above.}
\plot{nsea00}{width=1\columnwidth}{PERM receiver wavefield of one areal shot.}
\plot{nsea04}{width=1\columnwidth}{Optimized velocity model. Left: cross-line at 3500 m. Right: in-line at 3760 m.}
\plot{nsea02}{width=1\columnwidth}{Angle-gathers after CAM with the optimized velocity. Top: zero-angle section. Bottom: ACDCIGs selected roughly at the same x position in the section above.}
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Conclusions}
We show that pre-stack exploding-reflector model synthesizes wavefields that provides migrated images with correct kinematics while decreasing data size. Data reduction is achieved by combining the modeling experiments and is controlled by the number of subsurface offsets that will be computed during areal-shot migration of PERM data. Implicit to PERM is that reflectors must be identified. 3D Pre-stack interpretation can be cumbersome, but it allows, for instance, avoiding the use of reflectors with low signal-to-noise ratios in migration velocity estimation. In 3D, PERM data size can be one order of magnitude smaller if cross-line subsurface-offsets are to be computed. Further data size reduction by one order of magnitude is achieved if the initial conditions are computed with common-azimuth migration.
\par
In addition to the reduced data size, the ability PERM wavefields naturally have to be used in a target-oriented way greatly accelerates the definition of the velocity model. The 3D example with real data confirms the accuracy of the velocity solution achieved by migration velocity analysis by wavefield extrapolation using PERM wavefields.
%

\section*{Acknowledgments}
The authors thank TotalFinaElf for providing the data used in this paper and the sponsors of the Stanford Exploration Project for their financial support.
\par
\cleardoublepage
\newpage
\onecolumn
\bibliography{allrefs}
\bibliographystyle{seg} 
\end{document}
