\title{Fast 3D wave-equation migration-velocity analysis using the prestack exploding-reflector model}

\lefthead{Guerra and Biondi}
\righthead{WEMVA 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}
An accurate depth-velocity model is the key for obtaining good quality and reliable depth images in areas of complex geology. In such 
areas, velocity-model definition should use methods that describe the complexity of wavefield propagation, such as focusing and defocus
ing, multiple arrivals, and frequency-dependent velocity sensitivity. Wave-equation tomography in the image space has the ability to ha
ndle these issues because it uses wavefields as carriers of information. However, its high cost and low flexibility for parametrizing t
he model space has prevented its routine industrial use.
This thesis aims at overcoming those limitations by using new wavefields as carriers of information: the image-space generalized wavef
ields. These wavefields are synthesized by using a pre-stack generalization of the exploding-reflector model. Cost of wave-equation tom
ography in the image space is decreased because only a small number of image-space generalized wavefields are necessary to accurately d
escribe the kinematics of velocity errors and because these wavefields can be easily used in a target-oriented way. Flexibility is natu
rally incorporated into wave-equation tomography in the image space by using these wavefields because their modeling have as the initia
l conditions some key  selected reflectors, allowing a layer-based parametrization of the model space.
To use the image-space generalized wavefields in wave-equation tomography in the image space, the method is extended from the shot-pro
file domain to the image-space generalized-sources domain. In this new domain, the velocity updates are very fast. Migration with the o
ptimized velocity model provides good quality and reliable depth images, as can be seen in a 3D-field data example.
\end{abstract}

\section{introduction}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{prestack exploding-reflector model} 
A prestack image computed with an inaccurate velocity model by wavefield-extrapolation methods presents energy departing from the zero subsurface offset. In conventional Born modeling, the entire prestack image is spatially convolved with the source wavefield, generating the receiver wavefield, whose back-propagation and collection at the surface originate data. In order to guarantee illumination similar to the original, several shot records need to be modeled. Notice that both 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 single subsurface-offset common-image gather (SODCIG) of a reflector. Since we do not know beforehand which shots contribute to forming the image at a point in the subsurface, again we would have to perform several modeling experiments to synthesize these shot records. However, ideally, we would like to synthesize a small amount of data with the condition that migration produces the same kinematics 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
Areal data is generated by PERM, propagating energy of a prestack image.  PERM operates in a way similar to the exploding-reflector model (ERM). However, since migration of PERM data can produce a prestack image, PERM can be considered a generalization of ERM. The modeling of PERM  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({\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 same background velocity used to produce the initial prestack image, $D$ and $U$ represent the downward and the upward 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 for the initial conditions for modeling downward and upward 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 its spatial coordinates, and PERM is equivalent to ERM.  
\par
The initial conditions are obtained by rotating the original unfocused SODCIGs according to the apparent geological dip of the reflector.  This operation corrects for the image-point dispersal due to velocity inaccuracy, causing events with different reflection angles from the same reflection point to be imaged at different locations. This is a consequence of not shifting the wavefields along the apparent geological-dip direction \cite[]{biondi:1283} when applying the multi-offset image condition \cite[]{james2002}. The 2D geometry that explains the rotation is depicted in Figure \ref{fig:pdisp} for the cases of migration velocity which is slower (Figure \ref{fig:pdisp}a) and faster (Figure \ref{fig:pdisp}b) than the correct one.
\par
\plot{pdisp}{width=0.9\textwidth}{Geometry for the computation of SODCIGs. 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}.}
The angles $\gamma + \alpha$ and $\gamma - \alpha$ are the 
source and receiver ray angles, respectively. They are the propagation directions of the wavefields locally at the image point. In 2D, $\alpha$ and $\gamma$ are
 related to slopes in the pre-stack image according to
 \begin{eqnarray}
 \tan \alpha = -\frac{dz_m}{dx_m}
 \label{eq:talpha}
 \end{eqnarray}
 \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 the 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)$.
 To align the initial conditions with the geological dip, the dip along the subsurface-offset axis is changed, 
 yielding the new subsurface offset $\widetilde{h_{x_s}}$ and $\widetilde{h_{x_r}}$ for the initial conditions of the modeling of downward and upward PERM wavefields,
  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}
 The solutions of equations \ref{eq:theta1} and \ref{eq:theta2} define new slant-stack operations, which 
  reduce the dimensionality of the decomposed pre-stack image 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}})$.
In 3D, the cross-line subsurface offsets also must be rotated.  
A complete description of the 3D rotation is given by \cite{guerra_thesis}.
\par
PERM as described by equations \ref{eq:perm_bD} and \ref{eq:perm_bU} generates sufficient data to image a single SODCIG. However, to fully image a reflector, we would have to model data that could potentially be even greater than the original data. Using the linearity of wavefield propagation, the modeling experiments can be combined in a way that PERM data size is decreased. Hence, a group of SODCIGs equally spaced by certain $x$ and $y$ decorrelation distances 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 SODCIGs for a single reflector to be used as the initial conditions for the modeling of combined downward and upward wavefields, $\widehat{D}$ and $\widehat{U}$, respectively. The selection of SODCIGs can be thought of as the multiplication of the pre-stack image by a 2D--comb function, which is 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. If the periods of the 2D--comb function are chosen to be the decorrelation distances, no crosstalk between wavefields from different SODCIGs occurs during migration. The decorrelation distances correspond to twice the $x$ and $y$ subsurface-offset range. If more than one reflector is used, each reflector is separately injected into the modeling, so that wavefields from different reflectors do not cross-correlate during imaging, avoiding reflector crosstalk. These procedures to avoid crosstalk limit data size reduction with PERM. To overcome this limitation, two different strategies can be used, mitigating the effects of both kinds of crosstalk. 
\par
One strategy mitigates the reflector crosstalk, benefiting from the fact that the image of a reflector is formed at the zero time of wavefield propagation when the migration velocity is accurate. When the migration velocity is inaccurate, the focusing of the image departs from the zero time, and this can be potentially used to update the velocity model \cite[]{sava:S209,yang:3715}. Therefore, depending on the magnitude of the velocity errors and the distance between reflectors, selecting events with propagation time close to zero can avoid cross-correlating 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}. Phase-encoding is usually performed in the data space. However, wavefields can also be phase-encoded in the image space, using stochastic phase functions parametrized according to the model-space coordinates and selected reflector. Phase-encoding the modeling experiments mitigates not only reflector crosstalk but also crosstalk from different SODCIGs. This is achieved 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 the reflector $j$ by identifying and windowing it in the pre-stack image $I$, and $\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. Sequences used in third- and fourth-generation cellphones can also be used to phase-encode the modeling \cite[]{guerra2008b}. The pseudo-random phase-encoding function causes the frequency components of the initial conditions to be randomly injected into the modeling. Downward and upward image-space phase-encoded wavefields (ISPEW) initiated at the same SODCIG and from the same reflector are equally encoded, whereas downward and upward ISPEW initiated at different SODCIGs and from the different reflectors have different codes assigned to them. According to the law of large numbers, crosstalk attenuation is more effective for more random realizations. 
\par
PERM wavefields can be collected at any depth, which can be the top of a region of inaccurate velocity, having key reflectors within this region as the initial conditions. By doing so, velocity update will be restricted to the inaccurate velocity region. The use of  key reflectors characterizes a target-oriented strategy to perform MVA by wavefield extrapolation.  These two features along with the decrease of data size yielded by PERM drastically decreases the cost of MVA by wavefield extrapolation. 
\par
To illustrate the computation 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 the apparent geological dip. Data comprise 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 exactly the same as the correct velocity model up to the black horizon. From this horizon down, velocity is strongly smoothed and multiplied by 0.9. The curvature of the reflectors in the SODCIG 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}). 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. Notice that crosstalk is almost completely avoided. Migration 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}{Velocity models for the Marmousi example: a) Smooth velocity model used to model the 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.}  
\plot{ismarm03}{width=0.8\textwidth}{PERM wavefields for the Marmousi example: a) Upward wavefield. b) Downward wavefield.}  
\plot{ismarm06a}{width=0.8\textwidth}{ISPEW for the Marmousi example: a) Upward wavefield. b) Downward wavefield.}  
\plot{ismarm04}{width=0.8\textwidth}{Pre-stack image computed with PERM wavefields and background velocity model using: a) the conventional imaging condition (equation \ref{eq:I_perm}), and b) the time-windowed imaging condition (equation \ref{eq:wic_bI}). Reflector crosstalk is avoided when reflectors are sufficiently separated. However, some residual crosstalk is still present (RC). Notice the phase difference of the PERM image due to the squaring of the wavelet when compared to the windowed reflectors of Figure \ref{fig:ismarm02}b.}
\plot{ismarm07}{width=0.8\textwidth}{Pre-stack  images computed with: a) Four random realizations of ISPEWs, and b) a single random realization.}
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. Notice that in the ADCIGs corresponding to images computed with PERM wavefields and ISPEWs, less events are present since we selected four reflectors to model these wavefields. The corresponding reflectors in the original shot-profile image are highlighted by the green boxes in Figure \ref{fig:ismarm08}a. The residual moveout information in the four panels is very similar. However, some crosstalk not entirely rejected by the time-windowed imaging condition causes the residual-moveout information lose resolution in Figure \ref{fig:ismarm08}b when compared to Figures \ref{fig:ismarm08}c.  
\plot{ismarm08}{width=0.8\textwidth}{ADCIGs (top) and $\rho$-panels (bottom) corresponding to images computed by: a) Shot-profile migration of 360 shot gathers, b) areal-shot migration of 35 PERM wavefields using the time-windowed imaging condition, c) areal-shot migration of 44 ISPEWs corresponding to four random realizations, and d) areal-shot migration of 11 ISPEWs corresponding to a single random realization. The moveout information is basically the same.} 
\par
The cost for obtaining the images in Figure \ref{fig:ismarm08} widely varies. For instance, migrating 11 ISPEW is approximately 30 times faster than migrating all the 375 original shots. Using the same wavefields, computing the gradient of the objective function of MVA by wavefield-extrapolation can be 60 times faster.  This difference in performance is even more dramatic when we consider that several iterations of migrations and gradient computations are performed during migration-velocity  optimization by wavefield-extrapolation.

\section{MVA by wavefield-extrapolation}
Optimization of the migration velocity is a non-linear inverse problem. When using wavefield-extrapolation methods, it searches for an optimal background velocity which minimizes an objective function defined in the image space. The residual that characterizes the objective function 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 the linearized-residual prestack-depth migration \cite[]{sava2003}, and it is used in the w
ave-equation migration-velocity analysis \cite[]{biondi:1723, sava2004}. Also, the perturbed image can be computed by the differential-semblance-optimization (DSO) operator \cite[]{symes:654}, and it is 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,
\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 at zero offset, by weighting the background image with the absolute value of the subsurface offset. In the angle domain, DSO penalizes flattening of events, by taking the derivative of the background image along the aperture angle.
\par
Gradient-based optimization techniques, such as the nonlinear conjugate-gradient method, can be used in the optimization, to which the gradient must be explicitly computed. The gradient ${\bf \nabla}J_{\rm DSO}$ of the objective function  with respect to the slowness ${\bf 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), which is composed of several operators and linearly maps the slowness perturbation ${\Delta s}$ to the perturbed image.
\par
 WETOM has been evaluated in the source and receiver domain \cite[]{sava2004}, in the shot-profile domain \cite[]{shen2004}, and in the generalized-sources domain both in the data and in the image space \cite[]{tang2008}. We discuss WETOM in the image-space generalized-sources domain, using 35 ISPEW modeled from the 12 reflectors of Figure \ref{fig:ismarm201}, selected from the background image of Figure \ref{fig:ismarm102}a. 
\plot{ismarm201}{width=0.6\textwidth}{Selected reflectors used to model 35 pairs of ISPEW.}
\par
In the areal-shot migration of image-space generalized wavefields, wavefields 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}({\bf x},{\bf p},\omega) = 0 \\
{\widehat D}(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}({\bf x},{\bf p},\omega) = 0 \\
{\widehat U}(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}({\bf x},{\bf p},\omega)$ is the image-space generalized downward wavefield for a single frequency, ${\widehat U}({\bf x},{\bf p},\omega)$ is the image-space generalized upward 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 using phase encoding or not, and collected at $z=z_{min}$, which denotes the top of a target zone. The PERM data serve as the boundary conditions of equations \ref{aeq_wveq_bD} and \ref{aeq_wveq_bU}, respectively. Snapshots of image-space generalized background wavefields are shown in Figure \ref{fig:ismarm202}.
 \plot{ismarm202}{width=0.6\textwidth}{Snapshots of image-space generalized background wavefields: a) downward, and b) upward wavefields.}
 \par
 The cross-correlation imaging condition produces the 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.6\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 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 ${\widehat D}_0({\bf x}-{\bf h},{\bf p},\omega)$ and ${\widehat U}_0({\bf x}+{\bf h},{\bf p},\omega)$ are the image-space generalized background downward and upward wavefields computed with the background slowness; and $\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 downward wavefield and the image-space generalized perturbed upward wavefield, respectively. These image-space generalized perturbed wavefields are the response to a slowness perturbation. The image-space generalized perturbed wavefields 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 are shown in Figure \ref{fig:ismarm205}.
  \plot{ismarm205}{width=0.6\textwidth}{Snapshots of image-space generalized perturbed wavefields: a) source, and b) receiver.}
  Since the boundary conditions at $z_{min}$ are null, no scattering occurs at depth levels shallower than $z_{min}$. The wavefields in the right-hand side of equations \ref{aeq_wveq_dD} and \ref{aeq_wveq_dU} are the image-space generalized scattered downward and upward wavefields, respectively, which result from the interaction of the image-space generalized background wavefields with a slowness perturbation according to
  \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 image-space generalized perturbed source and receiver wavefields are used along with the precomputed image-space generalized background source and receiver wavefields in equation \ref{aeq_ic_dI} to generate the perturbed image (Figure \ref{fig:ismarm204}). The image-space generalized background source and receiver wavefields are obtained by recursively solving equations \ref{aeq_wveq_bD} and \ref{aeq_wveq_bU} using the background slowness.
  \plot{ismarm204}{width=0.6\textwidth}{Perturbed image computed with equation \ref{aeq_ic_dI}.}
  \par
  The adjoint WETOM ${\bf T}'$ is obtained by the folowing operations. First, the adjoint-imaging condition is applied to compute the image-space generalized perturbed source and receiver 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 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 source wavefield is cross-correlated with the image-space generalized scattered source wavefield, and the image-space generalized perturbed receiver wavefield is cross-correlated with the image-space generalized scattered receiver wavefield to generate the slowness perturbation according to
  \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}.
  \plot{ismarm206}{width=0.8\textwidth}{Slowness perturbation from back-projected image perturbations computed with 35 ISPEWs.}
\par
When performing velocity optimization using DSVA, for instance, 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. However, neither the phase nor the amplitudes of the DSO perturbed image are consistent with those of the perturbed image computed by the forward one-way ISWET operator. These differences prevent the use of linear conjugate-gradient methods, and therefore the objective function computed with DSO is typically minimized by nonlinear optimization methods.
  \par
  The gradient of the wave-equation tomography objective function is sensitive to amplitude variations of the prestack image caused by uneven illumination. Since these amplitude variations are not related to velocity inaccuracy, we should ideally attenuate them using some sort of illumination compensation scheme \cite[]{valenciano:WCA109,tang:WCA95}. Instead, to prevent these amplitude variations we apply 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. 
  \par
  Since the gradient is not properly scaled, we normalize it with the diagonal operator $\bf F$, which is the smallest value of the initial slowness divided by the maximum absolute value of the gradient. To improve and sometimes guarantee convergence, we would like to limit the velocity update from one iteration with respect to its previous values. In other words, we would like 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, the final gradient 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
  In the next section, we use the operators described above to smooth and scale the gradient in a velocity-optimization problem for a 3D field-data example.

  \section{3D field-data example }
%  We illustrate the use of image-space generalized wavefields in ISWET for of the Marmousi velocity model. The initial velocity model (Figure \ref{fig:ismarm101}b) differs from the true one (Figure \ref{fig:ismarm101}a) only below the black horizon. We modeled 35 ISPEW, which were collected at a depth of 1500 m. Therefore, the wavefield propagation in ISWET is performed between this depth and the maximum depth of 3000 m, characterizing a target-oriented strategy. We show some of the optimization results having 1500 m as the initial depth.
%  \plot{ismarm101}{width=0.8\textwidth}{Optimized velocity model.}
%  \par
%  A nonlinear conjugate-gradient is used in the optimization. The maximum allowed velocity change between iterations is 10 $\%$. The nodes of the B-spline smoothing of the gradient are separated by 480 m and 160 m in x and z, respectively. Two function evaluations are performed in each iteration, and if the objective function does not decrease, a 50$\%$-smaller step length is used.
%  \par
%  The final velocity model is shown in Figure \ref{fig:ismarm304}. Velocity has increased by nearly a 12 $\%$. Since the bottom of the velocity model is poorly constrained by reflectors, in this part of the model the velocity update is almost zero.  
%  The initial and final background images are shown in Figure  \ref{fig:ismarm301}. The final image is more focused, and the pulled-up reflectors at x = 6000 m are better positioned.
%  \plot{ismarm304}{width=0.8\textwidth}{Optimized velocity model.}
%  \plot{ismarm301}{width=0.6\textwidth}{a) Initial and b) final background image.}
%  \par
%  The cross-plot of Figure \ref{fig:ismarm303} presents the evolution of the objective function, normalized by the highest value. Iteration stopped after nine iterations, because the variation of the objective function was less than the predefined value of 0.002$\%$.
%  \sideplot{ismarm303}{width=0.8\textwidth}{Evolution of the objective function for the 11-ISPEWs case (blue diamonds) and 35-ISPEWs case (red squares).}
%  \par
%  As previously mentioned, ISWET solves for the long-wavelength component of the velocity model, which is consistent with the Born approximation. The long-wavelength component of the velocity model is responsible for the kinematics of the wavefield propagation. To evaluate how accurate the results are, we smooth the slowness derived from the original velocity model of Figure \ref{fig:ismarm101}a in a way similar to the gradient of the objective function. This smoothed version of the original velocity model is compared to the initial velocity model and the optimized velocity models within the box of Figure \ref{fig:ismarm307}a. Histograms of the velocity ratio between the smoothed true velocity model and the initial velocity model, the smoothed true and the 11-ISPEWs optimized velocity model, and the smoothed true and the 35-ISPEWs optimized velocity model are shown in Figure \ref{fig:ismarm307}b-d. The concentration around one in Figures \ref{fig:ismarm307}c-d indicates that the long-wavelength components of the velocity model were appropriately recovered by ISWET.
%  \plot{ismarm307}{width=0.8\textwidth}{a) Smoothed version of the true velocity model, b) histogram of the velocity ratio between the smoothed true velocity model and the initial velocity model, c) histogram of the velocity ratio between the smoothed true velocity model and the 11-ISPEWs optimized velocity model, and d) histogram of the velocity ratio between the smoothed true velocity model and the 35-ISPEWs optimized velocity model.}
%  \par
%  The images computed with shot-profile migration using the optimized velocity model (Figure \ref{fig:ismarm305}) show focused reflectors, and the pull-up has been corrected. Compare with  the image computed with the true velocity model in Figure \ref{fig:ismarm306}.
%  \plot{ismarm305}{width=0.6\textwidth}{Images computed with shot-profile migration using the optimized velocity model of Figure \ref{fig:ismarm304}.}

 % \plot{ismarm306}{width=0.6\textwidth}{Image computed with shot-profile migration using the true velocity model.}
 % \par
 %MVA by wavefield-extrapolation using image-space generalized sources is extremely inexpensive when compared to that with the original shot profiles. Considering shot profiles datumized to a depth of 1500 m, we estimate that, in the present case, velocity optimization would be nearly 15 times faster than in the shot-profile domain.


\bibliographystyle{seg}
\bibliography{allrefs}

