%% This file is automatically generated. Do not edit!
\documentclass[manuscript]{geophysics}

\usepackage{color,amsmath}

\setfigdir{Figures}

\begin{document}
\title{Fast 3D migration-velocity analysis by wavefield extrapolation using the prestack exploding-reflector model}

\lefthead{Guerra and Biondi}
\righthead{MVA with the prestack exploding-reflector}

\renewcommand{\thefootnote}{\fnsymbol{footnote}}
%\ms{GEO-2008-0443}

\address{
\footnotemark[1]Formerly Stanford Exploration Project, Geophysics Department, 
Stanford University, Stanford, CA 94305. Presently Petrobras, Rio de Janeiro, Brazil.\\
E-mail: claudio.guerra@petrobras.com.br \\
\footnotemark[2]Stanford Exploration Project, Geophysics Department, 
Stanford University, Stanford, CA 94305. E-mail: biondo@sep.stanford.edu \\
%Manuscript submitted on 12/24/2008. Revised manuscript submitted on 4/29/2009.
}

\author{Claudio Guerra\footnotemark[1] and Biondo Biondi\footnotemark[2]}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
In areas of complex geology, migration-velocity estimation should use methods that describe the complexity of wavefield propagation, such as focusing and defocusing, multipathing, and frequency-dependent velocity sensitivity. Migration-velocity analysis by wavefield extrapolation has the ability to address these issues because, in contrast to ray-based methods, it uses wavefields as carriers of information. However, its high cost and lack of flexibility with respect to model parametrization and to target-oriented analysis have  prevented its routine industrial use. We overcome those limitations by using new types of wavefield as carriers of information: the image-space generalized wavefields. These wavefields are synthesized from a prestack image computed with wavefield-extrapolation methods, using the pre-stack exploding-reflector model. Cost of migration-velocity analysis by wavefield extrapolation is decreased because only a small number of image-space generalized wavefields are necessary to accurately describe the kinematics of velocity errors and because these wavefields can be easily used in a target-oriented way. Flexibility is naturally incorporated because modeling these wavefields has as the initial conditions  selected reflectors, which allows use of a horizon-based parametrization of the model space. In a 3D example of the North Sea, we show that using wavefields synthesized by the prestack exploding-reflector model greatly improves efficiency of migration-velocity analysis by wavefield extrapolation, while yielding a final migration-velocity model that is accurate as evidenced by well focused and structurally reasonable reflectors.
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{introduction}
Migration-velocity analysis (MVA) solves for earth models that explain, under some norm, a residual-moveout parameter measured in the post-migrated domain. Ray-based MVA is the most widely used method because of its low cost and flexibility in parametrizing the model space, which can improve convergence to a geologically reasonable velocity model. However, in the presence of high lateral velocity variation and irregular interfaces, ray-based MVA methods are prone to fail because they do not describe the entire complexity of wavefield propagation. Hence, in those situations it is desirable to use wavefields to define the velocity model.
\par
The past decade has seen the development of MVA by wavefield extrapolation \cite[]{biondi:1723,sava2004,shen2004,shen:VE49}, which promises to overcome the limitations imposed by the high-frequency approximation present in ray methods. MVA by wavefield extrapolation is solved in the image space, using wavefields as carriers of information. Despite its theoretical superiority to ray methods, this relatively new technology has been rarely used in 3D projects \cite[]{fei:2283} because of its higher cost and because it is less flexible than its ray-based counterpart in parameterizing the velocity model and in solving velocity problems in a target-oriented fashion. We overcome these limitations by using image-space generalized ~ wavefields \cite[]{yax2008, guerra2009a, guerra_thesis}. These wavefields are initiated in the image space, having as the initial condition key reflectors selected from a prestack image obtained with a wavefield-extrapolation migration. 
\par
The image-space generalized wavefields are computed with the pre-stack exploding-reflector model (PERM) \cite[]{biondi2006}, which is an extension of the exploding-reflector model \cite[]{loewenthal01} in the sense that migration of PERM data can produce a prestack image. The fundamental idea of PERM is to model downgoing and upgoing wavefields necessary to provide a prestack image of a single point on a reflector with kinematic information consistent with that of the original prestack image. Synthesizing data sets from a limited portion of the subsurface has been previously explored by \cite{Rietveld02} and further developed into the concept of Common-Focus Point (CFP) \cite[]{berkhout:1895, thorbecke:S273}. 
\par
To image the entire reflector or group of reflectors, many independent PERM experiments would have to be performed, which could potentially generate a volume of data considerably larger than the original. Because of the linearity of wavefield propagation, however, the modeling experiments can be combined so as to yield a dramatic data reduction, especially in 3D. When linearly combining wavefields after encoding with stochastic phase functions, the ${\it image-space ~ phase-encoded~ wavefields}$ (ISPEW) are generated. Since image-space generalized wavefields are initiated at selected reflectors, a horizon-based strategy to parametrize the velocity model can be naturally incorporated into MVA by wavefield extrapolation. Moreover, these wavefields can be generated for a small subset of the image cube and collected at any depth during their propagation, also making a target-oriented approach easily incorporated. The goal of this paper is to show that using PERM wavefields lends flexibility and efficiency to MVA by wavefield extrapolation.
\par
This paper is organized as follows. First, we describe PERM and how to synthesize ISPEW. Next, we introduce MVA by wavefield extrapolation in the image-space generalized-sources domain, illustrating with the Marmousi model. Finally, we run MVA by wavefield extrapolation, using only 30 pairs of ISPEW, to estimate the 3D-migration-velocity model for an area with complex salt body in the North Sea. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{the prestack exploding-reflector model}
A prestack image computed by wavefield-extrapolation methods with an inaccurate velocity model presents energy that undesirably departs from the zero subsurface offset. The zero subsurface offset represents the zero lag of the spatial crosscorrelation of the downgoing and upgoing wavefields when applying the multi-offset imaging condition \cite[]{rickett:883}. In conventional Born modeling, this entire prestack image is spatially convolved with the source wavefield, generating the receiver wavefield, whose back-propagation and collection at the surface creates data. Source and receiver wavefields must be  propagated with the same inaccurate velocity model used in the original migration to preserve the correct kinematics. Now, let us consider the modeling of data necessary to correctly image a reflector in a subsurface-offset-domain common-image gather (SOCIG). Since we do not know beforehand which shots contribute to forming the image at a point in the subsurface, we would have to perform several Born modeling experiments to guarantee illumination similar to that of the original SOCIG. Ideally, however, we would like to synthesize a small amount of data with the condition that migration produces kinematics consistent to that of the initial prestack reflector. This is the fundamental idea of the prestack exploding-reflector model (PERM) \cite[]{biondi2006}, whose primary objective is to accelerate migration-velocity analysis.
\par
PERM operates in a way similar to the exploding-reflector model. Areal data are generated having as the initial conditions prestack images computed with wavefield-extrapolation methods. The initial conditions differ from that of the exploding-reflector model, which considers only the zero subsurface offset. The prestack information contained in the initial conditions of PERM enables migration of PERM data to produce a prestack image. In this sense, PERM can be considered a generalization of the exploding-reflector model.% CFP technology aims at synthesizing data gathers that contain all the events that would have been recorded at the surface if we had placed a point source in the subsurface. A CFP gather is composed of traces resulting from the time convolution between the shot gather and a focusing operator plus stacking of the stationary phases. The focusing operator is computed by modeling the seismic response of one point at the subsurface. In contrast with CFP, PERM data result from the wave propagation of the prestack image at one point at the subsurface, without operating on the original shot gathers. In this sense, PERM data is more closely related to the CFP operator itself than to the CFP data.
\par
Heuristically, we can think of synthesizing PERM data as generating wavefields with a new acquisition configuration. Conventional acquisition is characterized by sources and receivers located at, or close to, the surface. Born modeling and its adjoint (migration) explain how these data are generated and how a migrated image is produced by crosscorrelating source and receiver wavefields, respectively. The PERM acquisition configuration is characterized by source (downgoing) and receiver (upgoing) wavefields that are initiated on reflectors of the conventional migrated image and collected at a certain depth, in the form of areal shots. However, Born modeling does not explain the generation of PERM wavefields. Differently from Born modeling, in which a prestack image is convolved with the downgoing wavefield to generate the upgoing wavefield, in PERM both wavefields are the result of upward propagating the prestack image: time goes backwards for the downgoing wavefield and forwards for the upgoing wavefield.
\par
A fundamental assumption of PERM is that the prestack energy is present in the SOCIG, irrespective of velocity error, so that it can be used as the initial conditions. This assumption is easily fulfilled in the case of shot-receiver migration, in which the subsurface-offset range equals the data offset. However, in shot-profile migration the subsurface-offset range must be large enough so that the unfocused prestack energy is adequately imaged.
\par
The modeling of PERM  wavefields can be carried out by any wavefield-extrapolation 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({\bf x},\omega;{\bf x_m}) = I_{\scriptscriptstyle D}({\bf x_m},{\bf h}) \\
D(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({\bf x},\omega;{\bf x_m}) = I_{\scriptscriptstyle U}({\bf x_m},{\bf h}) \\
U(x,y,z=z_{\rm max},\omega;{\bf x_m}) = 0 \end{array} \right.,
\label{eq:perm_bU}
\end{eqnarray}
where $\omega$ is the radial frequency, $s_0$ is the background slowness used to produce the initial prestack image, $D$ and $U$ represent the downgoing and the upgoing PERM wavefields, respectively, $I_{\scriptscriptstyle D}$ and $I_{\scriptscriptstyle U}$ are subsets of the prestack image centered at ${\bf x_m}$ for a single reflector, suitable as the initial conditions for modeling downgoing and upgoing PERM wavefields, respectively. The information along the subsurface-offset $\bf h$ dimension of the initial conditions is projected onto the $\bf x$ dimension for the modeling. The subsurface offset is parameterized as ${\bf h}=(h_x,h_y)$, where $h_x$ and $h_y$ are the inline- and the crossline-subsurface offsets, respectively. Notice that if the prestack image has energy focused at zero subsurface offset, the initial conditions can be parameterized only by the spatial coordinates, and PERM is equivalent to exploding-reflector modeling.  
\par
Shot-profile and areal-shot migrations by wavefield extrapolation compute pre-stack images by means of the multi-offset imaging condition (Rickett and Sava, 2002), in which source and receiver wavefields are laterally shifted prior to time correlation. Ideally, wavefields should be shifted along the geological dip direction, so that SOCIGs do not suffer from image-point dispersal in the presence of dip and inaccuracies in the migration velocity  \cite[]{biondi:1283}. The image-point dispersal causes events with different reflection angles from the same reflection point in the subsurface to be imaged at different locations. Unfortunately, obtaining SOCIGs along the geological dip is computationally demanding since wavefields must be stored at various depths. Furthermore, accurate dip information is difficult to obtain, especially when events cross because of velocity inaccuracy. The image-point dispersal in 2D is illustrated in Figure \ref{fig:pdisp} for migration velocity that is slower (Figure \ref{fig:pdisp}a) and faster (Figure \ref{fig:pdisp}b) than the true velocity. For simplicity, let us consider constant velocity in the vicinity of the image point, so source and receiver rays are straight. When the migration velocity is too low, the reflector is imaged at a shallower depth. The image point $\mathsf {I_{hx}}$ computed with horizontal shifts of the wavefields is shifted down-dip with respect to the image point $\mathsf {I_{hg}}$ computed with shifts along the apparent geological dip. The geological dip is called apparent because of the migration velocity error. The point $\mathsf I$ is where source and receiver rays cross at an angle that is  twice the apparent reflection angle $\gamma$. Source and receiver rays cross deeper than the image points, causing events to curve downward in the SOCIG. When the migration velocity is too high, the reflector is imaged at a greater depth. The image point $\mathsf {I_{hx}}$ is shifted up-dip with respect to $\mathsf {I_{hg}}$. Source and receiver rays cross at a depth shallower than the image points, causing events to curve upward in the SOCIG. 
\par
\plot{pdisp}{width=0.9\textwidth}{Geometry for the computation of SOCIGs. Source, receiver and image points are labeled with $\mathsf S$, $\mathsf R$ and $\mathsf I$, respectively. The subscript $\mathsf {hx}$ corresponds to subsurface offsets computed with horizontal shift. The subscript $\mathsf {hg}$ corresponds to subsurface offsets computed by shifting along the apparent geological dip $\alpha$. a) Underestimated velocity, and b) overestimated velocity. Modified from \cite{biondi:1283}.}
\par
Differently from Born modeling, in which the adjoint multi-offset imaging-condition operator is applied to initiate the upgoing wavefield, in PERM  this operator is not present. For this reason, to account for the image-point dispersal the initial conditions need to be preprocessed by rotating the original unfocused SOCIGs according to the apparent geological dip of the reflector.  The derivation of the rotation according to the apparent geological dip, which yields dip-independent initial conditions,  is given in the Appendix A. Figure \ref{fig:app01} illustrates the effect that rotating the initial conditions has on the kinematics of the PERM migrated image. Figure \ref{fig:app01}a is the initial image computed with a velocity 10 $\%$ slower than the true velocity of 1000 m/s. The left panel is the SOCIG at  $x = 0$ m, and the right panel is the zero subsurface-offset section.  Figure \ref{fig:app01}b is the migrated image using PERM wavefields synthesized from the non-rotated initial conditions of Figure \ref{fig:app01}a with the correct velocity. Figure \ref{fig:app01}c is the rotated image used as the initial conditions for the modeling of PERM wavefields.  Figure \ref{fig:app01}d is the migrated image using PERM wavefields synthesized from the rotated initial conditions of Figure \ref{fig:app01}c with the correct velocity.  Note how well focused is the dipping reflector in the SOCIG of Figure \ref{fig:app01}d compared to that of Figure \ref{fig:app01}b. Since the horizontal reflector is not affected by the image-point dispersal, its focusing is similar in both situations. 
\plot{app01}{width=0.8\textwidth}{Rotation of the initial conditions.  a) Shot-profile migration of 801 split-spread shots 10 m apart with velocity $10\%$ slower than the true velocity; b) areal-shot migration of PERM wavefields initiated using the non-rotated initial conditions of Figure \ref{fig:app01}a; c) rotated initial conditions for the modeling of upgoing PERM wavefields; and d) areal-shot migration of PERM wavefields initiated using the rotated initial conditions of Figure \ref{fig:app01}c.}
\par 
PERM as described by equations \ref{eq:perm_bD} and \ref{eq:perm_bU} generates sufficient data to image a SOCIG. However, to fully image a reflector, we would have to model much more than the original amount of data. Using the linearity of wavefield propagation, the modeling experiments can be combined so as to decrease PERM data size. Hence, a group of equally spaced SOCIGs is simultaneously injected into the modeling according to
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 s_0^2({\bf x})-|{\bf k}| ^2}\right) \widehat {D}({\bf x},\omega) = \widehat {I_D}({\bf x-h}) \\
\widehat{D}(x,y,z=z_{\rm max},\omega) = 0 \end{array} \right.,
\label{eq:comb_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) \widehat{U}({\bf x},\omega) = \widehat {I_U}({\bf x+h}) \\
\widehat{U}(x,y,z=z_{\rm max},\omega) = 0 \end{array} \right.,
\label{eq:comb_bU}
\end{eqnarray}
where $\widehat{I_D}$ and $\widehat{I_U}$ are the combination of SOCIGs for a single reflector to be used as the initial conditions for the modeling of combined downgoing  and upgoing wavefields, $\widehat{D}$ and $\widehat{U}$, respectively. 
\par
The selection of SOCIGs can be thought of as the multiplication of the pre-stack image by a 2D--comb function, which is shifted laterally to select new sets of SOCIGs to initiate the modeling of other pairs of combined wavefields. All the points on the reflector are used in the modeling after shifting along one period of the sampling function in the $x$ and $y$ directions. Consequently, the number of modeling experiments equals the number of lateral shifts of the sampling function.
\par
Depending on the spatial periods of the sampling function, crosstalk between wavefields modeled from different SOCIGs can occur during migration. If the periods of the 2D--comb function are chosen to be the decorrelation distances of twice the $x$ and $y$ subsurface-offset range,  no crosstalk is generated. Otherwise, the image will be contaminated by crosstalk, as described in Appendix B. If more than one reflector is used, each reflector must be separately injected into the modeling so that wavefields from different reflectors do not cross-correlate during imaging, avoiding reflector crosstalk. However, these procedures to avoid crosstalk limit data size reduction with PERM. To overcome this limitation, two different strategies can be used.
\par
One strategy mitigates the reflector crosstalk, benefiting from the fact that the image of a reflector is formed at the zero lag of the cross-correlation between downgoing and upgoing wavefields, which, when using PERM data, corresponds to the zero time of wavefield propagation when the migration velocity is accurate. When migration velocity is inaccurate, focusing of the image departs from the zero time \cite[]{sava:S209,yang:3715}. This can be used to update the velocity model \cite[]{wang:WCA159}. Depending on the magnitude of the velocity errors and the distance between reflectors, selecting events with propagation time close to zero can avoid cross-correlation of unrelated events. This is the principle of the time-windowed imaging condition, which for a single pair of PERM wavefields reads
\begin{eqnarray}
I_{\scriptscriptstyle {W}}({\bf x},{\bf h}) = \sum_{\frac{-t_{\scriptscriptstyle w}}{2} \leq t \leq \frac{t_{\scriptscriptstyle w}}{2}} {\mathcal F}^{-1} \left [{\widehat D}^{*}({\bf x}-{\bf h}) \right ] {\mathcal F}^{-1}\left [{\widehat U}({\bf x}+{\bf h}) \right ],
\label{eq:wic_bI}
\end{eqnarray}
where $t_{\scriptscriptstyle w}$ is the length of the time window. When using one-way propagators, the wavefields are inverse Fourier transformed to time by ${\mathcal F}^{-1}$.  
\par
Another strategy is to phase-encode the modeling experiments. Phase-encoding is a well stablished technique for decreasing the cost of seismic imaging by linearly combining the shot records \cite[]{Schultz,whitmore,romero,Sun01,liu01,duquet2006}. In general, phase-encoding is performed in the data space. PERM wavefields  phase-encoded in the image space are called the ${\it image-space~ phase-encoded~ wavefields}$ (ISPEW). We achieve this by phase-encoding the initial conditions  ${\widetilde I}$, according to
\begin{eqnarray}
{\widetilde I}(\widehat{\bf x}{\bf ,h},{\bf q},\omega) = \sum_m \sum_j \delta({\widehat{\bf x}-m\Delta{\bf x}}) \beta(\widehat{\bf x},j,{\bf q},\omega) W_j(\widehat{\bf x}{\bf ,h}) I(\widehat{\bf x}{\bf ,h}),
\label{eq:ispew_Itil}
\end{eqnarray}
where $W_j$ selects reflector $j$ by identifying and windowing it in the pre-stack image $I$, $\Delta{\bf x}$ is the period of the 2D sampling function, $\beta(\widehat{\bf x},j,{\bf q},\omega)$ is a pseudo-random phase-encoding function defined as
\begin{eqnarray}
\beta(\widehat{\bf x},j,{\bf q},\omega) = e^{i\epsilon(\widehat{\bf x},j,{\bf q},\omega)},
\label{eq:ispew_rnd}
\end{eqnarray}
with $\epsilon(\widehat{\bf x},j,{\bf q},\omega)$ usually being a uniformly distributed pseudo-random sequence with zero mean, and $\bf q$ is the vector of random realizations .Other sequences, Gold codes for instance, used in third- and fourth-generation cellphones, can also be used to phase-encode the modeling  \cite[]{guerra2008b}. Phase-encoding the modeling experiments mitigates crosstalk between wavefields modeled from different reflectors as well as crosstalk between wavefields modeled from different SOCIGs. Downgoing and upgoing ISPEW initiated at the same SOCIG and from the same reflector are encoded identically, whereas downgoing and upgoing ISPEW initiated at different SOCIGs and from the different reflectors have different codes assigned to them. In accordance with the law of large numbers, computing more random realizations of ISPEW increases the attenuation of  crosstalk. 
\par
Image-space generalized wavefields can be collected at any depth, such as the top of a region of especially inaccurate velocity. Unfocused reflectors due to the velocity inaccuracy within this region are selected and used as the initial conditions for the modeling. By doing so, velocity update can be restricted to just the region of inaccurate velocity. The use of  key reflectors characterizes a target-oriented strategy for performing MVA by wavefield extrapolation.  These two features (using the top of a target region to collect the wavefields and initiate them from key reflectors) along with the decrease of data size significantly decreases the cost of MVA by wavefield extrapolation. 
\par
We use the Marmousi model to illustrate the generation of PERM data and ISPEW. A prestack image computed with the inaccurate velocity model of Figure \ref{fig:ismarm01}b is used as the initial conditions after rotation according to the apparent geological dip. The data consist of 375 split-spread shots, with 6000-m maximum offset, computed by one-way Born modeling using the smooth Marmousi velocity model of Figure \ref{fig:ismarm01}a, along with the reflectivity derived from the true Marmousi model. The inaccurate velocity model is identical to the correct velocity model down to the black horizon. From this horizon down, velocity is strongly smoothed and multiplied by 0.9 (Figure \ref{fig:ismarm01}b). The curvature of the reflectors in the SOCIG and the pull-up of reflectors in the center of Figure \ref{fig:ismarm02}a reflect the velocity inaccuracy. We use the four reflectors shown in Figure \ref{fig:ismarm02}b to model 35 pairs of PERM data (Figure \ref{fig:ismarm03}) and 11 pairs of ISPEW (Figure \ref{fig:ismarm06a}). Since 17 subsurface offsets are computed in the migration, which spans less than half of the spatial sampling period, no SOCIG crosstalk is expected to occur when using 35 pairs of PERM data. Migration of the 35 pairs of PERM data using the conventional imaging condition is shown in Figure \ref{fig:ismarm04}a and using the time-windowed-imaging condition is shown in Figure \ref{fig:ismarm04}b. Note that this modified imaging condition almost completely avoids reflector crosstalk. Migration using the conventional imaging condition of the 11 pairs of ISPEW is shown in Figure \ref{fig:ismarm07}. Reflectors are satisfactorily imaged and crosstalk is dispersed throughout the image.  
\par
\plot{ismarm01}{width=0.8\textwidth}{Edited velocity models for the Marmousi example: a) Smooth velocity model used to model the one-way Born data. b) Background velocity model used to migrate the Born data, and to model and migrate PERM data.}
\plot{ismarm02}{width=0.8\textwidth}{a) Pre-stack image computed with the background velocity model. b) Selected reflectors from the background image to perform modeling of wavefields. The panels on the left are SOCIGs taken at distance = 5 km. The panels on the right are the zero subsurface-offset sections.}  
\plot{ismarm03}{width=0.8\textwidth}{PERM wavefields for the Marmousi example: a) Downgoing wavefield. b) upgoing wavefield.}  
\plot{ismarm06a}{width=0.8\textwidth}{ISPEW for the Marmousi example: a) Downgoing wavefield. b) upgoing wavefield.}  
\plot{ismarm04}{width=0.8\textwidth}{Pre-stack image computed with 35 pairs of PERM wavefields and background velocity model using: a) the conventional imaging condition, and b) the time-windowed imaging condition. Reflector crosstalk is avoided when reflectors are sufficiently separated. However, some residual crosstalk is still present, labeled as RC in the panel on the left.}
\plot{ismarm07}{width=0.8\textwidth}{Pre-stack  images computed with 11 ISPEWs.}
To evaluate the correctness of the kinematic information, we compute the residual parameter $\rho$, using equation D-7 in \cite{biondi:1283}. Velocity is slower than the correct one for $\rho < 1.0$ and faster for $\rho > 1.0$. Residual-moveout panels from images computed with the original shot records (Figure \ref{fig:ismarm08}a), with 35 pairs of PERM data (Figure \ref{fig:ismarm08}b), and with 11 ISPEW (Figure \ref{fig:ismarm08}c), are used to analyze the quality of the moveout information. Note that in the angle common-image gathers (ACIGs) corresponding to images computed with 35 pairs of PERM data and 11 pairs of ISPEW, fewer events are present than in the shot-profile migration since only four reflectors were selected to model these wavefields. The corresponding reflectors in the original shot-profile image are highlighted by green boxes in Figure \ref{fig:ismarm08}a. The residual moveout information in the three panels is similar; however, some crosstalk not entirely rejected by the time-windowed imaging condition causes the residual-moveout information to lose vertical resolution in Figure \ref{fig:ismarm08}b when compared to Figure \ref{fig:ismarm08}c.  
\plot{ismarm08}{width=0.8\textwidth}{ACIGs (top) and $\rho$-panels (bottom) corresponding to images computed by: a) Shot-profile migration of 375 shot gathers, b) areal-shot migration of 35 PERM wavefields using the time-windowed imaging condition, c) areal-shot migration of 11 ISPEWs. The moveout information is basically the same.} 
\par
The cost in computing the images in Figure \ref{fig:ismarm08} varies widely. For instance, migrating 11 ISPEW is approximately 30 times faster than migrating all the 375 original shots. This difference in performance takes an especial significance in MVA by wavefield extrapolation, which typically involves  several iterations of migration and gradient computation. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{MVA by wavefield extrapolation}
Optimization of the migration velocity is a nonlinear inverse problem. When using wavefield-extrapolation methods, we search for an optimal background velocity that minimizes an objective function defined in the image space. The residual used to compute the gradient of the objective function with respect to velocity is represented by the perturbed image $\Delta {\widehat I}$, which, in turn, is derived from the background image ${\widehat I}_0$ computed with the background slowness $s_0$. 
\par
The perturbed image can be computed by linearized residual prestack-depth migration \cite[]{sava2003} and used in wave-equation migration-velocity analysis \cite[]{biondi:1723, sava2004}. Alternatively, the perturbed image can be computed by the differential-semblance-optimization (DSO) operator \cite[]{symes:654} and used in the differential-semblance velocity optimization (DSVA) \cite[]{shen2004, shen:VE49}. 
\par
Under the $\ell_2$ norm, the DSVA objective function $J_{\rm DSO}$ is
\begin{eqnarray}
J_{\rm DSO} = \frac{1}{2}||\Delta {\widehat I}||_2 = \frac{1}{2}||{\bf H}[{\widehat I}_0]||_2,
\label{jdso}
\end{eqnarray}
where $\bf H$ is the  DSO operator either in the subsurface-offset domain or in the angle domain. In the subsurface-offset domain,  DSO penalizes energy departing from zero offset by weighting the background image with the absolute value of the subsurface offset. In the angle domain, DSO penalizes the lack of flattening of reflectors by taking the derivative of the background image along the aperture angle.
\par
DSVA is solved by gradient-based optimization techniques, such as the nonlinear conjugate-gradient method, in which the gradient must be explicitly computed. The gradient ${\bf \nabla}J_{\rm DSO}$ of the objective function (equation \ref{jdso}) with respect to the slowness ${s}$ is
\begin{eqnarray}
{\bf \nabla} J_{\rm DSO} = {\bf T}^{\prime}{\bf H}^{\prime}{\bf H} I_0, \label{eq_grad_dso}
\end{eqnarray}
where ${\bf T} = \left.\frac{\partial I}{\partial s}\right|_{s=s_0}$ is the wave-equation tomographic operator (WETOM) evaluated at $s=s_0$. It is composed of several operators and linearly maps the slowness perturbation ${\Delta s}$ into the perturbed image.
\par
 WETOM has been evaluated in the source-receiver domain \cite[]{sava2004}, in the shot-profile domain \cite[]{shen2004}, and in the generalized-sources domain \cite[]{tang2008}. We apply WETOM in the image-space generalized-sources domain. To illustrate the process we use 35 ISPEW modeled from 12 reflectors  selected from the background image computed with the velocity model of Figure \ref{fig:ismarm101}b. This velocity model is obtained by smoothing the Marmousi velocity model (Figure \ref{fig:ismarm101}a) and multiplying it by 0.9 only below the black horizon. In Figure \ref{fig:ismarm201}, is represented the initial conditions for the modeling of upgoing PERM wavefields. 
\plot{ismarm201}{width=0.8\textwidth}{Selected reflectors after rotation according to the apparent geological dip used to model 35 pairs of ISPEWs. SOCIGs below are shown beneath their approximate x-position. Notice that dipping reflectors from the top panel are not symmetric in the SOCIGs below as a consequence of the rotation, and that horizontal reflectors are not affected by the rotation.}
\plot{ismarm101}{width=0.8\textwidth}{Velocity models used in the explanation of MVA by wavefield extrapolation with the image-space generalized wavefields: a) True Marmousi velocity model. b) Background velocity model used in the initial migration and in the modeling of ISPEW.}
\par
In the areal-shot migration of image-space generalized wavefields, they are downward continued with the following one-way wave equations:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right){\widehat D_0}({\bf x},{\bf p},\omega) = 0 \\
{\widehat D_0}(x,y,z=z_{min},{\bf p},\omega) =  {\widetilde  D}(x,y,z=z_{min},{\bf p},\omega)
\end{array} \right., \label{aeq_wveq_bD}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right){\widehat U_0}({\bf x},{\bf p},\omega) = 0 \\
{\widehat U_0}(x,y,z=z_{min},{\bf p},\omega) = {\widetilde  U}(x,y,z=z_{min},{\bf p},\omega)
 \end{array} \right., \label{aeq_wveq_bU}
 \end{eqnarray}
 where ${\widehat D_0}({\bf x},{\bf p},\omega)$ is the image-space generalized background downgoing wavefield for a single frequency, ${\widehat U_0}({\bf x},{\bf p},\omega)$ is the image-space generalized background upgoing wavefield for a single frequency, ${\bf p}$ is the index of the areal shot, and ${\widetilde  D}(x,y,z=z_{min},{\bf p},\omega)$ and ${\widetilde  U}(x,y,z=z_{min},{\bf p},\omega)$ are the data synthesized with PERM and collected at $z=z_{min}$, which denotes the top of a target zone. These data serve as the boundary conditions of equations \ref{aeq_wveq_bD} and \ref{aeq_wveq_bU}, respectively. Snapshots of the background ISPEWs taken at -1.1 s for the downgoing wavefield and 1.1 s for the upgoing wavefield are shown in Figure \ref{fig:ismarm202}.
 \plot{ismarm202}{width=0.8\textwidth}{Snapshots of background ISPEWs: a) downgoing wavefield at -1.1 s, and b) upgoing wavefield at 1.1 s.}
 \par
 The crosscorrelation imaging condition produces the background image ${\widehat I}({\bf x},{\bf h})$ (Figure \ref{fig:ismarm203})
 \begin{eqnarray}
 {\widehat I}({\bf x},{\bf h}) = \sum_{\bf p}\sum_{\omega} {\widehat D}^*({\bf x}-{\bf h},{\bf p},\omega) {\widehat U}({\bf x}+{\bf h},{\bf p},\omega).
 \label{aeq_ic_bI}
 \end{eqnarray}
 \plot{ismarm203}{width=0.8\textwidth}{Background image computed with the image-space generalized background wavefields of Figure \ref{fig:ismarm202}.}
\par
 The perturbed image is derived by applying the differentiation-product rule to equation \ref{aeq_ic_bI}, which gives
 \begin{eqnarray}
 \Delta {\widehat I}({\bf x},{\bf h})  &=& \sum_{\bf p}\sum_{\omega} \Delta {\widehat D}^*({\bf x}-{\bf h},{\bf p},\omega) {\widehat U}_0({\bf x}+{\bf h},{\bf p},\omega) + \nonumber \\
 & & {\widehat D}_0^*({\bf x}-{\bf h},{\bf p},\omega) \Delta {\widehat U} ({\bf x}+{\bf h},{\bf p},\omega), \label{aeq_ic_dI}
 \end{eqnarray}
 where $\Delta {\widehat D}({\bf x}-{\bf h},{\bf p},\omega)$ and $\Delta {\widehat U}({\bf x}+{\bf h},{\bf p},\omega)$ are the image-space generalized perturbed downgoing wavefield and the image-space generalized perturbed upgoing wavefield, respectively. These perturbed wavefields are the response to a slowness perturbation. They satisfy the following one-way wave equations linearized with respect to the slowness:
 \begin{eqnarray}
 \left\{ \begin{array}{l}
 \left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {s}_0^2 ({\bf x})-|{\bf k}|^2}\right) \Delta {\widehat D}({\bf x},{\bf p},\omega)= {\widehat D}_{\scriptscriptstyle SC}({\bf x},{\bf p},\omega)
  \\
  \Delta {\widehat D}(x,y,z=z_{min},{\bf p},\omega) = 0 \end{array} \right. , \label{aeq_wveq_dD}
  \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) \Delta {\widehat U}({\bf x},{\bf p},\omega)={\widehat U}_{\scriptscriptstyle SC}({\bf x},{\bf p},\omega)\\
  \Delta {\widehat U}(x,y,z=z_{min},{\bf p},\omega) = 0 \end{array} \right.. \label{aeq_wveq_dU}
  \end{eqnarray}
  Snapshots of the image-space generalized perturbed wavefields taken at -1.1 s for the downgoing wavefield and for 1.1 s for the upgoing wavefield are shown in Figure \ref{fig:ismarm205}.
  \plot{ismarm205}{width=0.8\textwidth}{Snapshots of image-space generalized perturbed wavefields: a) downgoing wavefield at -1.1 s, and b) upgoing wavefield at 1.1 s.}
 No scattering occurs above the black horizon of Figure \ref{fig:ismarm101} because the initial conditions above it are null. The wavefields in the right-hand side of equations \ref{aeq_wveq_dD} and \ref{aeq_wveq_dU} are the image-space generalized scattered downgoing and upgoing wavefields, respectively, which result from the interaction of the image-space generalized background wavefields with a slowness perturbation given by
  \begin{eqnarray}
  {\widehat D}_{\scriptscriptstyle SC}({\bf x},{\bf p},\omega)=\frac{i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {s}_0^2({\bf x})}}}{\widehat D}_0({\bf x},{\bf p},\omega)
  \label{aeq_Dsc}
  \end{eqnarray}
  and 
  \begin{eqnarray}
  {\widehat U}_{\scriptscriptstyle SC}({\bf x},{\bf p},\omega)=\frac{- i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 s_0^2({\bf x})}}}{\widehat U}_0({\bf x},{\bf p},\omega).
  \label{aeq_Usc}
  \end{eqnarray}
  These wavefields are injected at every depth level during the recursive downward propagation of the perturbed  wavefields. The inverse square-root operator in the computation of the image-space scattered wavefields can be implemented either in the wavenumber domain \cite[]{dehoop} or in the space domain \cite[]{huang}.   We implement the scattering operator in the space domain by using the Taylor series  expansion of the denominator, which enables use of a laterally variable velocity.  The image-space generalized perturbed downgoing and upgoing  wavefields are used along with the precomputed image-space generalized background downgoing and upgoing wavefields in equation \ref{aeq_ic_dI} to generate the perturbed image (Figure \ref{fig:ismarm204}). 
  \plot{ismarm204}{width=0.8\textwidth}{Perturbed image computed with equation \ref{aeq_ic_dI}.}
  \par
  The adjoint WETOM operator ${\bf T}'$ is obtained by the following operations. First, the adjoint-imaging condition is applied to compute the image-space generalized perturbed downgoing and upgoing wavefields according to the following convolutions:
  \begin{eqnarray}
  \Delta {\widehat D}({\bf x},{\bf p},\omega) &=& \sum_{\bf h} \Delta {\widehat I}({\bf x},{\bf h}) {\widehat U}_0({\bf x}+{\bf h},{\bf p},\omega) \nonumber \\
  \Delta {\widehat U}({\bf x},{\bf p},\omega) &=& \sum_{\bf h} \Delta {\widehat I}({\bf x},{\bf h}) {\widehat D}_0({\bf x}-{\bf h},{\bf p},\omega). \label{aeq_aic_dI}
  \end{eqnarray}
  The image-space generalized perturbed wavefields are then upward propagated using the adjoint counterparts of equations \ref{aeq_wveq_dD} and \ref{aeq_wveq_dU}. At every depth of their upward propagation, the image-space generalized perturbed downgoing wavefield is cross-correlated with the image-space generalized scattered downgoing wavefield, and the image-space generalized perturbed upgoing wavefield is cross-correlated with the image-space generalized scattered upgoing wavefield to generate the slowness perturbation,
  \begin{eqnarray}
  \Delta {\widehat s} ({\bf x}) &=& \sum_{{\bf p}}\sum_{\omega} {\widehat D}_{\scriptscriptstyle SC}^*({\bf x},{\bf p},\omega) \Delta {\widehat D}({\bf x},{\bf p},\omega) + \nonumber \\
  & & {\widehat U}_{\scriptscriptstyle SC}^*({\bf x},{\bf p},\omega) \Delta {\widehat U}({\bf x},{\bf p},\omega).
  \label{aeq_asc}
  \end{eqnarray}
  The slowness perturbation for the Marmousi example computed in the image-space generalized-sources domain using ISPEWs is shown in Figure \ref{fig:ismarm206}. Notice that the slowness perturbation has values different from zero above the black line of Figure \ref{fig:ismarm206}, where no slowness perturbation is expected to occur. This is because the wavefields ${\widehat D}_{\scriptscriptstyle SC}$ and ${\widehat U}_{\scriptscriptstyle SC}$ in equation \ref{aeq_asc} are actually a phase-shifted and amplitude-modified version of the background wavefields. To improve convergence when updating the velocity model, energy above the black line in Figure \ref{fig:ismarm206} should be muted out. 
  \plot{ismarm206}{width=0.8\textwidth}{Slowness perturbation from the back-projected image perturbations computed with 11 ISPEWs.}
\par
When performing velocity optimization using DSVA, the perturbed image in equation \ref{aeq_aic_dI} is computed with DSO at the begining of every nonlinear iteration.
DSO easily automates migration-velocity optimization, but neither the phase nor the amplitudes of the DSO perturbed image is consistent with those of the perturbed image computed by the forward one-way WETOM \cite[]{vyas:4077,fei:4071}. These differences prevent the use of linear conjugate-gradient methods, so the objective function computed with DSO is typically minimized by nonlinear optimization methods.
  \par
  The gradient of the objective function in MVA by wavefield extrapolation is sensitive to amplitude variations of the prestack image caused by uneven illumination. Since these amplitude variations are not related to velocity inaccuracy, we ideally could attenuate them using some sort of illumination compensation scheme \cite[]{valenciano:WCA109,tang:WCA95}. Instead of compensating for illumination variations, we reduce the amplitude unbalancing by applying a B-spline smoothing to the gradient, which consists of representing the gradient as B-spline basis functions, using the adjoint operator $\bf B'$, and transforming it back to the Cartesian space, using the forward operator $\bf B$. 
  \par
 When using image-space generalized wavefields, a target-oriented strategy can be adopted if the velocity model is sufficiently accurate for shallower layers. A mask operator $\bf M$ is applied to the gradient, zeroing out amplitudes in the accurate velocity region, preventing the velocity model from being updated there. 
  \par
  Since the gradient is not properly scaled, we normalize it by the inverse of its maximum absolute value multiplied by the minimum slowness. This scaling is represented by the diagonal operator $\bf F$. To improve convergence, we would like to limit the size of the velocity update in each iteration relative to its previous values. This is intended to limit the new velocity to vary within a range defined by a percentage of the velocity from the previous iteration. This can be implemented by applying to the gradient either a nonlinear (since it depends on the velocity) diagonal operator, or a diagonal operator $\bf W$ linearized around the initial velocity. Therefore, in contrast with the gradient of equation \ref{eq_grad_dso}, the one  we use to update the velocity model in the 3D example that follows is
  \begin{eqnarray}
  {\bf \nabla} J_{\scriptscriptstyle DSO} = {\bf WFMBB^{\prime}} \left . {\bf T}^{\prime}\right |_{s=s_0} {\bf H}^{\prime}{\bf H} \left .{\widetilde I}\right |_{s=s_0} \label{eq_grad_is}.
  \end{eqnarray}
  \par

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{3D example }
We estimate the migration-velocity model with 3D ISPEW synthesized from prestack images of North Sea field data. The challenges for defining the velocity model for this dataset arise from a possibly irregular salt body, intense faulting, amplitude variations caused by irregular acquisition, short source-receiver offsets, and limited source-receiver azimuths. Because of the narrow azimuthal configuration, the 3D dataset was transformed to common-azimuth data by azimuth-moveout (AMO), and common-azimuth-migration (CAM)  images are used as the initial conditions for the modeling of relatively few 3D ISPEW. 
Thirty nodes of Dual Nehalem 5520 with 24Gb RAM were used, amounting to 240 CPUs. To generate 30 pairs of 3D ISPEWs, using 196 frequencies, took 10 minutes. On average, each iteration of the velocity optimization, consisting of one function evaluation, the gradient computation, and two additional function evaluations for the line search, took approximately two hours.   The velocity estimation presented here could readily have been performed interactively if a larger -- but commonly available in industrial setting -- number of computer nodes had been used to perform the computations. 
\par
The 3D North Sea dataset spans an area of approximately 55 km$^2$, with 13.5-km inlines and 4-km crosslines. It was acquired using dual sources at intervals of 25 m in the inline direction and 50 m in the crossline direction, with three cables 100-m apart and a maximum offset 3600 m. The limited crossline offsets resulting from this acquisition configuration impose limitations in the azimuthal distribution. The fold of coverage is spatially quite irregular, as can be seen in the fold of coverage map (Figure \ref{fig:3dex1025}a). Despite the regularization provided by AMO, the fold variation imposes an amplitude unbalancing, as evidenced by the acquisition footprint in a low-fold region occurring for inlines around 3500 m in Figure \ref{fig:3dex1025}b. This figure shows a  time slice at  2.8 s through the trace envelope for the stacked volume of the regularized data. The fold irregularity and amplitude imbalance in the data-space are evident in the amplitudes of the gradient of the objective function at approximately the same spatial positions. 

\plot{3dex1025}{width=0.8\textwidth}{a) Fold of coverage map, and b) time slice through the trace-envelope cube of stacked AMO data.} 
\par
The velocity model provided along with the data (Figure \ref{fig:3dex104}), which we call the ${\it original ~velocity~ model}$, has a general layered structure with an overhanging salt dome connected to a deeper layer with the same velocity as that of the salt dome. The shallower layer with a domed structure is a chalk layer, herein called ${\it chalk ~layer}$. As we shall see, the final velocity model we derive shows remarkable differences from the original velocity model.
\plot{3dex104}{width=0.8\textwidth}{Slices through the original velocity model. The black line represents the top of chalk and the blue line represents the base of chalk.} 
\par
To derive the initial velocity model, we refined the sediment velocity above the chalk layer using residual prestack depth-migration scans. Below the top of chalk, velocity was heavily smoothed using a 5000-m wide 2D median smoother. In addition, to increase the inaccuracy of the initial velocity model, velocity was scaled down by a factor of 0.9 (Figure \ref{fig:3dex105}). 
\par
A layer-stripping approach was used to define the velocity for the chalk layer, considering the velocity to be sufficiently accurate for the sediments above. Then, the top salt was interpreted, and a salt flooding procedure enabled the interpretation of the base of salt. Since no previous geological  information was available, the interpretation of the base of salt can be considered the main source of uncertainty of this example. Finally, a group of reflectors below the salt is used to define the velocity structure for the deeper part. 
\par
\plot{3dex105}{width=0.8\textwidth}{Slices through the  initial velocity model used in the velocity optimization. The black line corresponds to the top of chalk.} 
The initial velocity model produces the common-azimuth migration (CAM) image of Figure \ref{fig:3dex1066}, which shows sections of the volume for the zero subsurface offset, and Figure \ref{fig:3dex405}, which shows the zero subsurface offset on the top, and ACIGs at the bottom for the inline at 3120 m. The effects of migrating with a too low velocity are evidenced by poorly collapsed diffractions close to the salt flank, poorly imaged faults, and reflectors curving up in the ACIGs.
\plot{3dex1066}{width=0.8\textwidth}{Slices through the CAM image computed with the initial velocity model of Figure \ref{fig:3dex105}, showing poorly collapsed diffractions close to the salt flank and poorly imaged faults.} 
\plot{3dex405}{width=0.8\textwidth}{Inline 3120 m of the CAM image computed with the initial velocity model. On the top is the zero-subsurface  offset section, and ACIGs at the bottom.} 
\par
A window around the 3D pre-stack interpretation of the base of chalk selects amplitudes to input for the modeling of 30 ISPEW. We then interpolated the inline and crossline intervals of the CAM image from 20 m to 30 m, which are the inline and crossline intervals used for optimizing the migration velocity. 
\par
Since we only computed inline offsets, the CAM initial conditions can be continuously sampled in the crossline direction. Consequently, in this direction the modeling can be approximated by the exploding-reflector modeling. The continuous sampling in the crossline direction can reduce by at least an order of magnitude the number of image-space generalized wavefields to be synthesized \cite[]{guerra_thesis}. For the base of chalk, only 30 3D ISPEWs were modeled and  collected at 600-m depth. A pair of 3D ISPEW is shown in Figure \ref{fig:3dex109}. Since the number of subsurface offsets that will be used during velocity optimization is 25 and a spatial sampling in the inline direction of 30 SOCIGs is used to model the 3D ISPEWs, crosstalk is expected to be sufficiently attenuated. 
\plot{3dex109}{width=0.8\textwidth}{A pair of 3D downgoing  (a) and upgoing (b) ISPEWs computed for the base of chalk.} 
\par
In the velocity optimization for the chalk layer, velocity update is constrained to a maximum of 10$\%$ variation between iterations. All the wavefield propagation is performed between the depth at which the 30 3D ISPEWs were collected (600 m) and a maximum depth of 3300 m. Velocity is updated up to the top of the chalk layer.
\par
 To illustrate the problem of amplitude variations in the raw gradient, we compute the slowness perturbation without applying smoothing (Figure \ref{fig:3dex2102}a). The acquisition footprint is evident, specially around the inline coordinate 3500 m, which is a region of low-fold of coverage (Figure \ref{fig:3dex1025}). The gradient of the objective function must be smooth to yield velocity updates consistent with the Born approximation. We apply a B-spline smoothing with node intervals of 420 $\times$ 420 $\times$ 160 m in the inline, crossline, and depth directions, respectively. Applying B-spline smoothing on the DSO slowness perturbation of Figure \ref{fig:3dex2102}a mitigates the amplitude-variation problems and yields consistent slowness perturbations (Figure \ref{fig:3dex2102}b). 
\plot{3dex2102}{width=0.8\textwidth}{a) Slowness perturbation without smoothing, and b) Slowness perturbation after B-spline smoothing.} 
\par
Two runs of velocity optimization were necessary to define the velocity model for the chalk layer. The first one used wavefields modeled from the entire extension of base of chalk. Figure \ref{fig:3dex213} shows the evolution of the objective function for the first run. Since CAM with the resulting velocity revealed velocity inaccuracies close to the salt flanks, we modeled new wavefields using initial conditions limited to an area approximately 2-km wide around the salt edge. By doing this, we explore the localized nature of these wavefields, and the second run can be targeted for updating velocity close to salt flanks. In both runs, two function evaluations were used in the line search. 
\plot{3dex213}{width=0.8\textwidth}{Evolution of the DSVA objective function for the first run of velocity optimization for the base of chalk.}
\par
The optimized velocity models for the first and second runs of velocity optimization are shown in Figure \ref{fig:3dex216a}. Detailed views of CAM images computed with these velocity models are shown in Figure \ref{fig:3dexzoom}. Only slices through the zero subsurface offset are shown. Overall, the image computed with the velocity of the second run of velocity optimization is somewhat better focused. 
\par
\plot{3dex216a}{width=0.8\textwidth}{Slices through the optimized velocity for the chalk layer from: a) the first run of velocity optimization, and b) the second run of velocity optimization. The black line approximately represents the base of chalk.} 
\plot{3dexzoom}{width=0.8\textwidth}{Detailed view of slices through the zero subsurface offset of the CAM migrated images computed with velocity from the a) first run of velocity optimization, and b) second run of velocity optimization. On the left is an inline and on the right, a crossline. The orange line represents the top of chalk and the yellow line is the base of chalk. }
Once a sufficiently accurate velocity for the chalk layer was defined, we used salt flooding to  delineate the salt body. The top salt was interpreted and, below it, velocity was replaced by a constant value of 4500 m/s. CAM image computed with the salt-flooded velocity supported the interpretation of the base of the salt. Insufficient illumination resulting from the limited azimuthal coverage and irregular shape of the imaged salt body caused the base of salt to be discontinuous. In this situation, prior  geological information would be helpful to constrain the interpretation. The lack thereof is a source of uncertainty for defining the velocity model below the salt.
\par
After the base of salt was interpreted, salt velocity was delineated and the initial velocity was inserted below the salt and the chalk layer (Figure \ref{fig:3dex313}). This new velocity model produced the CAM  image in Figure \ref{fig:3dex400} from which seven deeper reflectors are used to model new ISPEWs to optimize sub-salt and sub-chalk velocity. These wavefields were collected at a depth of 1650 m, which is the minimum depth used in the velocity optimization. The maximum depth is 4800 m. 

\plot{3dex313}{width=0.8\textwidth}{Slices through the velocity volume after interpretation of the base of salt.} 
\plot{3dex400}{width=0.8\textwidth}{Slices through the CAM image computed with the migration velocity of Figure \ref{fig:3dex313}, showing in yellow the selected reflectors for the modeling of 30 3D ISPEWs to be used in the sub-salt velocity optimization: a) the zero subsurface offset, and b) inline 3120 m.}
\par
In the sub-salt velocity optimization, the intervals between B-spline nodes are 1050 m in the inline and crossline directions and 150 m in the depth direction. A maximum of 5$\%$ local velocity variation is allowed between iterations, and two function evaluations are performed in the line search. The final objective function dropped only 10$\%$ with respect to its initial value. This decrease is small compared to the 40$\%$ decrease for the velocity optimization for the chalk layer. Apart from the poorer illumination in the sub-salt and sub-chalk velocity optimization, which prevents the complete focusing at the zero subsurface offset, this difference can be explained by the different levels of crosstalk. For the chalk layer, crosstalk only  from different SOCIGs occurs, whereas for the sub-salt case, crosstalk from both different SOCIGs and different reflectors was generated. To decrease the influence of crosstalk on the objective function for the sub-salt-velocity optimization, more random realizations and an increase in the SOCIG sampling interval could have been used.  
%\sideplot{3dex402}{width=\textwidth}{Evolution of the DVSA objective function of the sub-salt velocity  optimization.}
\par
 The evolution of the velocity model through iterations is shown in Figure \ref{fig:3dex403}. Velocity gradually increased during the iterations, as expected from the moveouts observed in Figure \ref{fig:3dex400}b. CAM with the optimized velocity model (Figure \ref{fig:3dex403}d) can be seen in Figure \ref{fig:3dex404}. For comparison, CAM with the  original velocity model (Figure \ref{fig:3dex104}) and CAM with the inital velocity model (Figure \ref{fig:3dex105}) are shown in Figures \ref{fig:3dex406} and \ref{fig:3dex405}, respectively. When compared with the results using the initial velocity model, the improvements obtained with the optimized velocity model are clear: flatter angle gathers, better focusing of the reflectors, and better imaging of the faults. The improvements compared to the original (unmodified) velocity model are also clear: better focusing of reflectors and slightly flatter angle gathers below the salt as well as close to its flanks.
\plot{3dex403}{width=0.8\textwidth}{Slices through velocity models for: a) first iteration, b) third iteration, c) fifth iteration, and d) seventh iteration.} 
\plot{3dex404}{width=0.8\textwidth}{Inline 3120 m of the CAM image with the final velocity model after optimization for the chalk layer, salt flooding, and sub-salt velocity optimization. On the top is the zero subsurface offset section, and at the bottom ACIGs.} 
\plot{3dex406}{width=0.8\textwidth}{Inline 3120 m of the CAM image with the original velocity model. On the top is the zero subsurface offset section, and at the bottom ACIGs.} 
\par
Figure \ref{fig:3dex502} shows additional comparisons between the CAM images obtained with the initial, original, and final velocity models, for different inlines and crosslines. The images are displayed as vertical slices through the migrated cube and show the zero subsurface offset. The initial image is shown on the top, the original image in the middle, and the final image at the bottom.  The left panel of each figure shows the inline, and the right panel shows the crossline. Ovals highlight the main differences. Overall, the final image exhibts better quality than do the initial and original images, expressed by better focusing and continuity of reflectors in the chalk layer as well as in the sub-salt region.  The optimized velocity model allowed imaging of a complex fault system on the  right of the salt body, collapsing diffractions from the salt flank. In addition, overmigration effects at the base of chalk in the original image (easily detected in the right panel) are absent in the final image.
\plot{3dex502}{width=0.8\textwidth}{CAM images obtained with: a) the initial velocity model, b)  the original velocity model, and  c) the final velocity model. The final image shows better focusing and continuity for the sub-salt reflectors, better definition of subtle faults in the chalk layer.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{conclusions}
We introduced the image-space generalized wavefields, which are synthesized with the prestack exploding-reflector model. When phase encoded, these wavefields originate the image-space phase-encoded wavefields. 
Using 3D ISPEW has proven to greatly accelerate 3D-MVA by wavefield extrapolation thanks to the small number of wavefields needed to satisfactorily describe the kinematics of the prestack image and the fact that they can be computed in a target-oriented manner. Because of these computational efficiencies, we could solve a 3D migration-velocity problem employing only limited computational resources. Considering the larger computational resources available in the industry, using 3D ISPEW can turn MVA by wavefield extrapolation into an interactive process, which offers the promise of yielding more accurate and geologically reasonable solutions.
\par
In a 3D example from the North Sea, with only 30 ISPEWs we were able to obtain an accurate migration-velocity model using different strategies such as layer stripping and based on horizons. We also used a target-oriented strategy, which  allowed updating the velocity model within a both limited depth range and limited lateral extent. This is achieved by selecting as the initial conditions a specific portion of a reflector or group of reflectors that exhibit residual moveout.
\par
In addition to the computational gain, 3D ISPEWs here were able to provide reliable velocity updates. As expected, using the final velocity model produces a CAM image with quality superior to that obtained with the initial velocity model. Moreover, the image computed with the final velocity model is more accurate than that computed with the ray-based-derived original velocity model, with better focusing and continuity  of the sub-salt reflectors and better fault imaging. The computational efficiency, flexibility, and velocity accuracy obtained with data computed by the prestack exploding-reflector model enables the use of MVA by wavefield extrapolation as a routine procedure to define migration velocity in areas of complex geology. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Acknowledgements}
We thank the sponsors of the Stanford Exploration Project for the financial support. The authors thank TotalFinaElf for the field data and the Stanford Center for Computational Earth $\&$ Environmental Science in whose computational infrastructure all the computations were carried out. Bin Wang, Ken Larner, and two anonymous reviewers greatly contributed to improve the quality of this paper. Claudio Guerra is indebted to Petrobras for all kinds of support during his studies at Stanford University. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\append{Dip-independent initial conditions}
To accurately model PERM data it is crucial that the initial conditions are free of image-point dispersal. This is accomplished by rotating SOCIGs according to the geological dip. In Figure \ref{fig:pdisp},  $\gamma$ is the reflection angle,  $\alpha$ is the apparent geological dip,  $\gamma + \alpha$ is the source-ray angle, and $\gamma - \alpha$ is the receiver-ray angle. From the angle relationships, we can write 
\begin{eqnarray}
\tan \alpha = -\frac{dz_m}{dx_m}
\label{eq:talpha}
\end{eqnarray}
and
\begin{eqnarray}
\tan \gamma = -\frac{dz}{dh_x},
\label{eq:tgamma}
\end{eqnarray}
where the subscript $m$ in equation \ref{eq:talpha} refers to the local nature of the relationship. The solutions of differential equations \ref{eq:talpha} and \ref{eq:tgamma} define slant-stack paths, which allow us to transform the 2D pre-stack image $I(x,z,h_x)$ into $I(x,z,\alpha,\gamma)$ by angle decomposition according to the following integrals:
\begin{eqnarray}
\left . I(x,z,\alpha,\gamma) = \int_{x_{m_{i}}}^{x_{m_{f}}} \int_{-h_x}^{h_x} W(x_m-x) \frac{dI(x,z,h_x)}{dz} dx_{m} dh_{x} \right |_{\substack{z=z_h+h_x\tan\gamma\\~z=z_m+x_m\tan\alpha}},
\label{eq:sltstk1}
\end{eqnarray}
where the derivative with respect to $z$ is performed to recover the correct phase. A  local window $W(x_m-x)$ is used in the local slant-stack integral on $x_m$.
\par
Again, using simple trigonometry, we have 
\begin{eqnarray}
\tan (\gamma + \alpha) = \frac{\tan \gamma + \tan \alpha}{1 - \tan \gamma \tan \alpha},
\label{eq:theta0a}
\end{eqnarray}
\begin{eqnarray}
\tan (\gamma - \alpha) = \frac{\tan \gamma - \tan \alpha}{1 + \tan \gamma \tan \alpha}.
\label{eq:theta0b}
\end{eqnarray}
\par
To align the initial conditions with the geological dip, we need to change the dip along the subsurface-offset axis in accordance with the apparent geological dip, yielding the new subsurface offset $\widetilde{h_{x_s}}$ and $\widetilde{h_{x_r}}$ for the initial conditions of the modeling of source and receiver wavefield, respectively. This is accomplished by solving the following differential equations:
\begin{eqnarray}
\tan (\gamma + \alpha) = -\frac{dz}{d\widetilde {h_{x_s}}}, 
\label{eq:theta1}
\end{eqnarray}
\begin{eqnarray}
\tan (\gamma - \alpha) = \frac{dz}{d\widetilde {h_{x_r}}}.
\label{eq:theta2}
\end{eqnarray}
Equations \ref{eq:theta1} and \ref{eq:theta2} define new slant-stack operations which, in combination with equations \ref{eq:theta0a} and \ref{eq:theta0b}, reduce the dimensionality of the decomposed pre-stack image (equation \ref{eq:sltstk1}) by transforming $I(x,z,\alpha,\gamma)$ into $I_D(x,z,\widetilde{h_{x_s}})$ and $I_U(x,z,\widetilde{h_{x_r}})$.
\par
In 3D, the crossline offsets also must be rotated based on the apparent geological dip in the crossline direction in addition to the inline rotation. By assuming that source and receiver rays are coplanar such that they cross, the 3D transformation to the reflection-angle domain \cite[]{tisserant}  is given by 
\begin{eqnarray}
k_{h_x}=-k_{z} \sec\alpha_{y'}\tan\gamma,
\label{eq:3Dkhx0}
\end{eqnarray} 
and
\begin{eqnarray}
k_{h_y}=-k_{y'} \tan\gamma \tan \alpha_{x'},
\label{eq:3Dkhy0}
\end{eqnarray}
where $\alpha_{x'}$ and $\alpha_{y'}$ are the apparent geological dips in the inline and crossline  directions, respectively. Notice that, in contrast to the 2D case, the 3D transformation is dependent on the apparent geological dip. In the spatial domain, equations \ref{eq:3Dkhx0} and \ref{eq:3Dkhy0} define slant-stack transformations along the paths
\begin{eqnarray}
z=z_{h_{x}}+\frac{h_x}{\cos\alpha_{y'}}\tan\gamma,
\label{eq:3Dkhx1}
\end{eqnarray}
and
\begin{eqnarray}
z=z_{h_{y}}+h_y \tan\alpha_{x'}\tan\alpha_{y'}\tan\gamma,
\label{eq:3Dkhy1}
\end{eqnarray}
respectively. The term $\cos\alpha_{y'}$ in equation \ref{eq:3Dkhx1} stretches the inline-subsurface-offset axis, while the combination of terms $\tan\alpha_{x'}$ and  $\tan\alpha_{y'}$ in equation \ref{eq:3Dkhy1} can stretch or shrink the crossline-subsurface-offset axis.  Under the common-azimuth approximation, the 3D rotation is similar to the 2D rotation (equation \ref{eq:sltstk1}) except for the subsurface-offset stretching factor (equation \ref{eq:3Dkhx1}).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\append{Crosstalk generation}
 To decrease the number of modeling experiments, we can apply the concept of generalized sources and use the linearity of wavefield propagation to combine isolated SOCIGs and inject them simultaneously into one single model experiment as in equations \ref{eq:comb_bD} and \ref{eq:comb_bU}. For one particular frequency and a plane reflector so that the initial conditions are the same for modeling downgoing $\widehat D$  and upgoing  $\widehat U$ wavefields, the modeling of a pair of PERM wavefields starting from combined SOCIGs can be described by
\begin{eqnarray}
\widehat{D}(\boldsymbol\xi;\Delta {\bf x}) = \sum_{\bf x}\sum_{\bf h}\sum_m G_0(\boldsymbol\xi,{\bf x-h}) \delta({\widehat{\bf x}-m\Delta{\bf x}})I(\widehat{\bf x}{\bf ,h}),
\label{eq:perm_opDdx}
\end{eqnarray}
and
\begin{eqnarray}
\widehat{U}(\boldsymbol\xi;\Delta {\bf x}) = \sum_{\bf x}\sum_{\bf h}\sum_m G_0(\boldsymbol\xi,{\bf x+h}) \delta({\widehat{\bf x}-m\Delta{\bf x}})I(\widehat{\bf x}{\bf ,h}),
\label{eq:perm_opUdx}
\end{eqnarray}
where  $\displaystyle{\sum_m\delta({\widehat{\bf x}-m\Delta{\bf x}})}$ is the 2D--sampling function. The Green's function $G_0$ upward propagates the wavefields from the subsurface (represented by the the coordinates $\bf x$) up to the depth where the wavefields are collected (represented by the coordinates $\boldsymbol\xi$), using the background slowness $s_0$.
\par
The wavefields are recursively downward propagated in depth according to 
\begin{eqnarray}
\widehat{D}({\bf x};\Delta {\bf x}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x}) \widehat{D}({\boldsymbol\xi}),
\label{eq:perm_opDzdx} 
\end{eqnarray}
and
\begin{eqnarray}
\widehat{U}({\bf x};\Delta {\bf x}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x}) \widehat{U}({\boldsymbol\xi}),
\label{eq:perm_opUzdx}
\end{eqnarray}
where $G_1$ denotes the Green's function computed with a different slowness $s_1$. 
\par
The lateral shifts of the wavefields for the multi-offset imaging condition are represented by
\begin{eqnarray}
\begin{array}{lll}
\widehat{D}({\bf x-h};\Delta {\bf x})&=&\displaystyle{\sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x-h})} \\
&\times& \displaystyle{\sum_{\bf x'}\sum_{\bf h'}\sum_m G_0(\boldsymbol\xi,{\bf x'-h'}) \delta({\widehat{\bf x}-m\Delta{\bf x}})I(\widehat{\bf x}{\bf ,h'})},
\end{array}
\label{eq:perm_lshDdx} 
\end{eqnarray}
and
\begin{eqnarray}
\begin{array}{lll}
\widehat{U}({\bf x+h};\Delta {\bf x})&=&\displaystyle{\sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x+h})} \\
&\times& \displaystyle{\sum_{\bf x'}\sum_{\bf h'}\sum_m G_0(\boldsymbol\xi,{\bf x'+h'}) \delta({\widehat{\bf x}-m\Delta{\bf x}})I(\widehat{\bf x}{\bf ,h'})}.
\end{array}
\label{eq:perm_lshUdx}
\end{eqnarray}
\par
Applying the cross-correlation imaging condition to the wavefields of equations \ref{eq:perm_lshDdx} and \ref{eq:perm_lshUdx} gives
\begin{eqnarray}
\begin{array}{lll}
\widehat{I}({\bf x},{\bf h};\Delta {\bf x}) & = & \displaystyle{\sum_{\boldsymbol\xi'}\sum_{\bf x'}\sum_{\bf h'}\sum_m}\displaystyle{\sum_{\boldsymbol\xi''}\sum_{\bf x''}\sum_{\bf h''}\sum_n}G_0(\boldsymbol\xi',{\bf x'-h'}) G^*_1({\boldsymbol\xi'},{\bf x-h}) \\
&\times& G^*_1({\boldsymbol\xi''},{\bf x+h}) G_0(\boldsymbol\xi'',{\bf x''+h''}) \delta({\widehat{\bf x}}-m\Delta{\bf x}) \delta({\widehat{\bf x}}-n\Delta{\bf x}) \\
&\times& I(\widehat{\bf x}{\bf ,h'})I(\widehat{\bf x}{\bf ,h''}),
\end{array}
\label{eq:I_perm2dx}
\end{eqnarray}  
which can be recast as
\begin{eqnarray}
\begin{array}{lll}
\widehat{I}({\bf x},{\bf h};\Delta {\bf x}) & = & I({\bf x},{\bf h};\Delta {\bf x}) \\
&+& \displaystyle{\sum_{\boldsymbol\xi'}\sum_{\bf x'}\sum_{\bf h'}}\displaystyle{\sum_{\boldsymbol\xi''}\sum_{\bf x''}\sum_{\bf h''}\sum_{n \ne m}}G_0(\boldsymbol\xi',{\bf x'-h'}) G^*_1({\boldsymbol\xi'},{\bf x-h}) \\
&\times& G^*_1({\boldsymbol\xi''},{\bf x+h}) G_0(\boldsymbol\xi'',{\bf x''+h''}) \delta({\widehat{\bf x}}-m\Delta{\bf x}) \delta({\widehat{\bf x}}-n\Delta{\bf x}) \\
&\times& I(\widehat{\bf x}{\bf ,h'})I(\widehat{\bf x}{\bf ,h''}).
\end{array}
\label{eq:I_perm2dx1}
\end{eqnarray} 
The first term in the right-hand side of equation \ref{eq:I_perm2dx1} is the desired image we would obtain by independently modeling and migrating PERM wavefields. The second term represents crosstalk. To obtain a crosstalk-free image, the sampling period must be large enough that wavefields initiated at different SOCIGs do not correlate. This is achieved by setting $\Delta{\bf x} >$ twice the subsurface-offset range. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliographystyle{seg}
\bibliography{allrefs}

\end{document}
