\chapter{MVA using image-space generalized sources}
\label{chap:chap04}
In this chapter, I extend the theory of migration-velocity analysis (MVA) by wavefield extrapolation to the image-space generalized-sources domain. In this new domain, PERM wavefields and ISPEWs are the carriers of information for defining the migration velocity. These wavefields allow faster velocity updates than when conventional wavefields are used, as, for example, in shot-profile migration. The greater computational efficiency is possible due to the small data size of image-space generalized-source wavefields, in addition to their inherent capability for being used in a target-oriented manner. Moreover, with these wavefields, we can incorporate well-established strategies used in ray-based MVA, such as horizon-based tomography, into MVA by wavefield extrapolation. This new feature gives more flexibility to MVA by wavefield extrapolation and can improve the convergence to an optimal velocity. I illustrate the use of image-space generalized wavefields in velocity optimization with 2D examples. 
\newpage
\section{introduction}
Wave-equation tomography solves for earth models that better explain observed seismograms under some norm. There are two main categories, depending on the domain in which the objective function is minimized.
\par
In one category, known as waveform inversion, the objective function is minimized in the data space; data is modeled with the current model parameters and compared with the observed data. Extensive research has been being devoted to the application of waveform inversion to seismic exploration \cite[]{lailly83, tarantola1984, tarantola1987, mora:1211, woodward1992, pratt98}. A comprehensive overview of waveform inversion can be found in \cite {virieux:WCC1}. However, in spite of its maturity as a technology, examples of waveform inversion using 3D-field data are still very limited \cite[]{vigh:VE135,plessix:WCC149}. Waveform inversion is a highly nonlinear and ill-posed problem. Nonlinearity arises because the forward-modeling operator is a function of the searched model parameters. The ill-posedness is due to the many local minima of the objective function caused by the band-limited nature of seismic data, the incomplete acquisition, the presence of noise, and the incomplete physics of the operators. Because of the non-linearity and ill-posedness, the initial model plays a crucial role in waveform inversion. Ideally, the initial model should adequately describe the lower-frequency components of the velocity model, whose higher-frequency components are to be determined in subsequent iterations.
\par
Under the Born approximation, modeling of seismic data is linearly related to the reflectivity. However, Born modeling nonlinearly depends on the background velocity. These relations allow recasting the imaging of the subsurface into two separate, but related, problems: migration, which reconstructs the reflectivity given a background velocity, and velocity analysis, which determines the background velocity used in migration. These two problems are intimately related, and the degree of accuracy of the background velocity directly influences the quality of the migrated image. This relationship is explored in the second category of wave-equation tomography, here called image-space wave-equation tomography (ISWET). In ISWET, the objective function is optimized in the image space, and the residual is represented by an image perturbation. Similar to waveform inversion, ISWET is also nonlinear and ill-posed. The nonlinearity arises because the migration operator is a function of the background velocity.
\par
Because low--spatial-frequency components (i.e., background velocity) and high--spatial-frequency components (i.e.,  reflectivity) of the velocity model can be solved separately, the ill-posedness of ISWET is less severe than that of waveform inversion, which inverts for remarkably detailed velocity models. An essential feature of ISWET is that its objective function is intimately related to the final product of seismic processing, which is an image of the subsurface. The optimal velocity is the one which gives the best image. Moreover, less expensive one-way extrapolators can be used, in contrast with the two-way extrapolator used in waveform inversion.
\par
%Despite their differences, the two wave-equation tomography categories are related \cite[]{symes2008}. This author shows that ISWET can be regarded as a strategy to solve a partially linearized waveform inversion problem. The linearization enable us to compute the velocity perturbation by back-propagating the residuals of the non-linear problem using the adjoint of the Born scattering operator. The partially linearized problem solves for the long-wavelength components of the velocity model in a similar manner as ISWET and assures global convergence. Waveform inversion computes the short-wavelength components, providing accuracy.
%\par
Two major variants of ISWET are wave-equation migration velocity analysis (WEMVA) \cite[]{biondi:1723, sava2004a, sava2004b} and differential semblance velocity analysis (DSVA) \cite[]{shen2004, shen:VE49}. Both variants seek the optimal velocity by driving an image perturbation to a minimum. However, they differ in the way the image perturbation is computed and, consequently, in the numerical optimization scheme. As \cite{biondi2008} points out, WEMVA is not easily automated. The image perturbation is computed by the linearized-residual prestack-depth migration \cite[]{sava2003}, which uses a manually picked residual-moveout parameter. Since the perturbed image computed with the linearized-residual prestack-depth migration is consistent with the application of the forward wave-equation tomographic operator, WEMVA can be solved using a two-step approach. First, in a nonlinear iteration the background image is computed with the current velocity, a residual-moveout parameter is interpreted using enhanced versions of the background image, and  the current perturbed image  for the interpreted residual moveout is computed. Then, linear iterations using conjugate-gradients search for a perturbation in velocity that better explains the current perturbed image. The corresponding velocity solution is used to compute a new background image for the next nonlinear iteration.
\par
In DSVA, the perturbed image is computed by applying the fully automated differential-semblance operator (DSO) \cite[]{symes:654} to SODCIGs or ADCIGs. When applied to SODCIGs, DSO minimizes the energy not focused at zero-offset. When applied to ADCIGs, DSO minimizes energy  of the reflectors departing from flatness. Although DSO easily automates ISWET, it produces perturbed images that do not present the depth phase-shift introduced by the forward one-way ISWET operator. Moreover, the amplitude behavior of the perturbed image computed with DSO greatly differs from that of the perturbed image computed with 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, which require the explicit computation of the gradient of the objective function.
\par
The use of horizons that represent major velocity changes and present good signal-to-noise ratio, a common practice in MVA by ray-based methods \cite[]{stork:680, kosloff:1511, kosloff:1815, billette:1889, bob00}, defines two main strategies for velocity update: grid-based tomography and horizon-based tomography. A residual moveout parameter defined along horizons is back-projected through the entire velocity model for both strategies or, alternatively, can be restricted to certain layers in the horizon-based strategy. Moreover, the horizon-based strategy enables us to apply different regularization parameters for different layers, which can improve convergence. 
\par
Like waveform inversion, ISWET is a computationally demanding process. This computational cost is commonly decreased by using generalized sources \cite[]{shen:VE49,yax2008}. Because of the smaller data size, image-space generalized wavefields can drastically decrease the cost of ISWET. Also, as discussed in the previous chapters, image-space generalized wavefields can be propagated in a limited portion of the model space. Under the framework of migration velocity analysis, this allows their use in a target-oriented manner, since the wavefield propagation can be restricted to a region where the velocity model is inaccurate. Moreover, as these wavefields are initiated at some representative horizons, a horizon-based strategy is naturally incorporated into ISWET. Hence, the velocity model can be easily computed by a layer-striping scheme using wavefields initiated at an individual reflector or, more appropriately, at group of reflectors. Solving for a group of reflectors instead of using the layer-striping scheme avoids the propagation to deeper layers of velocity errors from shallower layers. 
\par
In this chapter, I describe image-space wave-equation tomography. Then, I discuss the tomographic operator in the shot-profile domain. Next, I extend the tomographic operator to the areal-shot domain using image-space generalized wavefields. Finally, I illustrate the use of image-space phase-encoded wavefields in DVSA for optimizing the velocity of the Marmousi model.  

\section{Image-space wave-equation tomography}
Image-space wave-equation tomography is a non-linear inverse problem. It searches for an optimal background velocity that minimizes an objective function defined in the image space. The objective function is represented by the residual $\Delta I$, hereafter called perturbed image, which is derived from the background image $I$ computed with the background velocity. The perturbed image represents the residual in the data-parameter space. The minimum of the perturbed image under some norm is unlikely to be global due to the non-linearity and ill-posedness of the problem. Therefore, constraints must be added to the model-parameter space by using a regularization operator. 
\par
The perturbed image can be computed by the DSO operator \cite[]{symes:654}, in the DVSA variant of ISWET, and by  linearized-residual prestack-depth migration \cite[]{sava2003}, in the WEMVA variant of ISWET. According to \cite{biondi2008}, a general form of the perturbed image can be expressed as
\begin{eqnarray}
{\Delta I} = I - {\bf F}[I],
\end{eqnarray}
where ${\bf F}$ is a focusing operator. Its application highlights the lack of focusing of the migrated image. Here and hereafter, we use bold capital letters for operators. Square brackets indicate the application of the operator to the argument. 
\par
In DVSA \cite[]{shen2004, shen:VE49} the focusing operator assumes form
\begin{eqnarray}
{\bf F} = {\bf 1 - H}, \label{eq_focus_dso}
\end{eqnarray}
where ${\bf 1}$ is the identity operator, and ${\bf H}$ is the DSO operator either in the subsurface offset domain or in the angle domain. The subsurface-offset-domain DSO focuses the energy at zero offset, whereas the angle-domain DSO flattens the ADCIGs. Hereafter, for the sake of simplicity, instead of velocity and velocity perturbation, we use slowness and slowness perturbation in the formulation of ISWET, because of the direct relation between the slowness perturbation and the perturbed image.
\par 
In WEMVA \cite[]{sava2004a, sava2004b}, the focusing operator is the linearized-residual prestack-depth migration \cite[]{sava2003} defined as
\begin{eqnarray}
{\bf F} = {\bf 1} + {\bf K}\left(\Delta \rho \right),
\end{eqnarray}
where $\Delta \rho = 1-\rho$, and $\rho$ is the ratio between the background slowness $s_0$ and the true slowness $s$. The differential-residual-migration operator ${\bf K}\left (\Delta \rho \right)$, which phase-shifts the image for different reflection angles and geological dips, is defined as
\begin{eqnarray}
{\bf K}\left(\Delta \rho\right) = \Delta \rho \left. \frac{\partial {\bf R}\left (\rho \right)}{\partial \rho} \right|_{\rho=1},
\label{eq:diff_res_mig0}
\end{eqnarray}
where ${\bf R}\left (\rho \right)$ is the residual-prestack-depth migration. Application of the chain rule to equation \ref{eq:diff_res_mig0} gives
\begin{eqnarray}
{\bf K}\left(\Delta \rho\right) = \Delta \rho \left. \frac{d{\bf R}\left (\rho \right)}{d k_z} \frac{d k_z}{d \rho} \right|_{\rho=1},
\label{eq:diff_res_mig1}
\end{eqnarray}
where $k_z$ is the vertical wavenumber. As \cite{sava2004} shows, all the elements in the right-hand side of equation \ref{eq:diff_res_mig1} can be easily computed. Notice that $\bf K$ implicitly depends on $\Delta s$ through $\Delta \rho$.
\par
Under the $\ell_2$ norm, the ISWET objective function is
\begin{eqnarray}
J = \frac{1}{2}||\Delta I||_2 = \frac{1}{2}||I-{\bf F}[I]||_2.
\end{eqnarray}
Gradient-based optimization techniques, such as the quasi-Newton method and the conjugate-gradient method, can be used to minimize the objective function $J$. The gradient of $J$ with respect to the slowness ${\bf s}$ is
\begin{eqnarray}
{\bf \nabla} J = \left(\frac{\partial {I}}{\partial {s}} - \frac{\partial {\bf F}[I]}{\partial s} \right)^{\prime} \left( I-{\bf F}[I]\right), \label{eq_obj_grad_gen}
\end{eqnarray}
where $^{\prime}$ denotes the adjoint. 
\par
In DVSA, the DSO operator ${\bf H}$ is independent of the slowness, so we have
\begin{eqnarray}
\frac{\partial {\bf F}({\bf I})}{\partial {\bf s}} = ({\bf 1}-{\bf H})\frac{\partial {\bf I}}{\partial {\bf s}}. \label{eq_focus_grad_dso}
\end{eqnarray}
Substituting equations \ref{eq_focus_dso} and \ref{eq_focus_grad_dso} into equation \ref{eq_obj_grad_gen} and evaluating the gradient at a background slowness yields
\begin{eqnarray}
{\bf \nabla} J_{\rm DSO} = \left.\left(\frac{\partial I}{\partial s}\right|_{s=s_0}\right)^{\prime}{\bf H}^{\prime}{\bf H} I_0, \label{eq_grad_dso}
\end{eqnarray}
where $I_0$ is the background image computed using the background slowness $s_0$.
\par
In WEMVA, the focusing operator depends on the slowness $s$. To simplify the gradient computation, we apply the focusing operator to the background image $I_0$ instead of $I$, and $\Delta \rho$ is interpreted on the background image, that is 
\begin{eqnarray}
{\bf F}[I_0]=I_0+{\bf K}(\Delta \rho)I_0.
\end{eqnarray}
With these assumptions, the WEMVA gradient is
\begin{eqnarray}
\nabla J_{\rm WEMVA} = - \left.\left(\frac{\partial I}{\partial s}\right|_{s=s_0}\right)^{\prime}{\bf K}(\Delta \rho)[I_0].
\end{eqnarray}
\par
%In Appendix A, is given a detailed derivation of the gradient of the ISWET objective function according to the adjoint-state method. We use the concept of the augmented Lagrangian, following \cite[]{plessix2006}.
%\par
The linear mapping from the slowness perturbation ${\Delta s}$ to the perturbed image ${\Delta I}$ is performed by the operator $\left.\frac{\partial{\bf I}}{\partial{\bf s}}\right|_{s=s_0}$. It is derived by keeping the zero- and first-order terms of the Taylor series expansion of the image around the background slowness
\begin{eqnarray}
{\Delta I} = {\bf T} {\Delta s},
\end{eqnarray}
where ${\Delta I} = {I} - {I_0}$ and ${\Delta s} = {s} - {s_0}$ and, ${\bf T} = \left. \frac{\partial{I}}{\partial{s}}\right|_{s=s_0}$ is the wave-equation tomographic operator. The tomographic operator has been evaluated either in the source and receiver domain \cite[]{sava2004} or in the shot-profile domain \cite[]{shen2004}. Because it is natural to extend ${\bf T}$ from the shot-profile domain to the generalized-sources domain, we review next the forward and adjoint tomographic operator in the shot-profile domain. 

\subsection{Shot-profile domain wave-equation tomographic operator}
My derivation of the wave-equation tomographic operator $\bf T$ in the shot-profile domain follows \cite{yax2008}. I use the Marmousi model to illustrate the components of the shot-profile domain wave-equation tomographic operator. In the present case, the velocity perturbation occurs below the horizon in Figure \ref{fig:ismarm101}, where the true and the background velocity models are shown. Data comprise 375 shots computed with the two-way wave equation.
\plot{ismarm101}{width=0.8\textwidth}{Marmousi velocity models: a) True velocity model. b) Background-velocity model computed by smoothing and scaling down the true model below the horizon indicated by the black line.}
\par
In shot-profile migration, both source and receiver wavefields are downward continued with the following one-way wave equations \cite[]{claerbout71}:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right)D({\bf x},{\bf x}_s,\omega) = 0 \\
D(x,y,z=0,{\bf x}_s,\omega) = f_s(\omega)\delta({\bf x}-{\bf x}_s) \end{array} \right., \label{eq_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)U({\bf x},{\bf x}_s,\omega) = 0 \\
U(x,y,z=0,{\bf x}_s,\omega) = Q(x,y,z=0,{\bf x}_s,\omega) \end{array} \right., \label{eq_wveq_bU}
\end{eqnarray}
where $D({\bf x},{\bf x}_s,\omega)$ is the source wavefield for a single frequency $\omega$ at image point ${\bf x}=(x,y,z)$, with the source located at ${\bf x}_s=(x_s,y_s,0)$; $U({\bf x},{\bf x}_s,\omega)$ is the receiver wavefield for a single frequency $\omega$ at image point ${\bf x}$ for the source located at ${\bf x}_s$; ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector; $f_s(\omega)$ is the frequency-dependent source signature; $\delta({\bf x}-{\bf x}_s)$ defines the point-source function at ${\bf x}_s$, which serves as the boundary condition of equation \ref{eq_wveq_bD}; and $Q(x,y,z=0,{\bf x}_s,\omega)$ is the recorded shot gather for the shot located at ${\bf x}_s$, which serves as the boundary condition of equation \ref{eq_wveq_bU}. Snapshots of both source and receiver background wavefields for a shot position around 5000 m are shown in Figure \ref{fig:ismarm102}.
\plot{ismarm102}{width=0.6\textwidth}{Snapshots of background: a) source, and b) receiver wavefields.}
\par
The cross-correlation imaging condition produces the image $I({\bf x},{\bf h})$:
\begin{eqnarray}
I({\bf x},{\bf h}) = \sum_{{\bf x}_s}\sum_{\omega} D^*({\bf x}-{\bf h},{\bf x}_s,\omega) U({\bf x}+{\bf h},{\bf x}_s,\omega).
\label{eq_ic_bI}
\end{eqnarray}
The background image at zero-subsurface offset and some SODCIGs positioned at their approximate location are shown in Figure \ref{fig:ismarm103}.  Notice the curvature of the reflectors at the region of inaccurate velocity. 
\plot{ismarm103}{width=0.6\textwidth}{Background image and SODCIGs computed with the background wavefields of Figure \ref{fig:ismarm102}.}
\par
The perturbed image shown in Figure \ref{fig:ismarm104} is derived by a simple application of the product rule to equation \ref{eq_ic_bI}, which gives
\begin{eqnarray}
\Delta I({\bf x},{\bf h})  &=& \sum_{{\bf x}_s}\sum_{\omega} \Delta D^*({\bf x}-{\bf h},{\bf x}_s,\omega) {U_0}({\bf x}+{\bf h},{\bf x}_s,\omega) + \nonumber \\
& & {D_0}^*({\bf x}-{\bf h},{\bf x}_s,\omega) \Delta U ({\bf x}+{\bf h},{\bf x}_s,\omega), \label{eq_ic_dI}
\end{eqnarray}
where ${D_0}({\bf x}-{\bf h},{\bf x}_s,\omega)$ and ${U_0}({\bf x}+{\bf h},{\bf x}_s,\omega)$ are the background source and receiver wavefields computed with the background slowness, and $\Delta D({\bf x}-{\bf h},{\bf x}_s,\omega)$ and $\Delta U({\bf x}+{\bf h},{\bf x}_s,\omega)$ are the perturbed source wavefield and perturbed receiver wavefield, respectively. The perturbed source and receiver wavefields are the response to a slowness perturbation. These 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 D({\bf x},{\bf x}_s,\omega)= D_{\scriptscriptstyle SC}({\bf x},{\bf x}_s,\omega)
 \\
\Delta D(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right. , \label{eq_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 U({\bf x},{\bf x}_s,\omega)=U_{\scriptscriptstyle SC}({\bf x},{\bf x}_s,\omega)\\
\Delta U(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right.. \label{eq_wveq_dU}
\end{eqnarray}
Snapshots of the perturbed source and receiver wavefields are shown in Figure \ref{fig:ismarm105}. Notice that the perturbed wavefields occur only in the region where velocity perturbation is different from zero.
\plot{ismarm104}{width=0.6\textwidth}{Perturbed image computed with equation \ref{eq_ic_dI}.}
\plot{ismarm105}{width=0.6\textwidth}{Snapshots of perturbed: a) source, and b) receiver wavefields.}
\par
The wavefields in the right-hand side of equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} are the scattered source and receiver wavefields, respectively, which result from the interaction of the background wavefields with a slowness perturbation according to
\begin{eqnarray}
D_{\scriptscriptstyle SC}({\bf x},{\bf x}_s,\omega)=\frac{i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {s}_0^2({\bf x})}}}D_0({\bf x},{\bf x}_s,\omega)
\label{eq_Dsc}
\end{eqnarray}
and
\begin{eqnarray}
U_{\scriptscriptstyle SC}({\bf x},{\bf x}_s,\omega)=\frac{-i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 s_0^2({\bf x})}}}U_0({\bf x},{\bf x}_s,\omega).
\label{eq_Usc}
\end{eqnarray}
These wavefields are injected at every depth level during the recursive propagation of the perturbed  wavefields. The perturbed source and receiver wavefields are used along with the precomputed background source and receiver wavefields in equation \ref{eq_ic_dI} to generate the perturbed image. The background source and receiver wavefields are obtained by recursively solving equations \ref{eq_wveq_bD} and \ref{eq_wveq_bU} using the background slowness.
\par
To evaluate the adjoint tomographic operator ${\bf T}'$, we first apply the adjoint imaging condition to compute the perturbed source and receiver wavefields, which are represented by the convolutions
\begin{eqnarray}
\Delta D({\bf x},{\bf x}_s,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h}) U_0({\bf x}+{\bf h},{\bf x}_s,\omega) \nonumber \\
\Delta U({\bf x},{\bf x}_s,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h}) D_0({\bf x}-{\bf h},{\bf x}_s,\omega). \label{eq_aic_dI}
\end{eqnarray}
These perturbed wavefields are upward propagated using the adjoint counterparts of equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU}. At every depth level of the upward propagation, the perturbed source wavefield is cross-correlated with the scattered source wavefield, and the perturbed receiver wavefield is cross-correlated with the scattered receiver wavefield to generate the slowness perturbation according to
\begin{eqnarray}
\Delta s ({\bf x}) &=& \sum_{{\bf x}_s}\sum_{\omega} D_{\scriptscriptstyle SC}^*({\bf x},{\bf x}_s,\omega) \Delta D({\bf x},{\bf x}_s,\omega) + \nonumber \\
& & U_{\scriptscriptstyle SC}^*({\bf x},{\bf x}_s,\omega) \Delta U({\bf x},{\bf x}_s,\omega).
\label{eq_asc}
\end{eqnarray} 
The slowness perturbation for the Marmousi example is shown in Figure \ref{fig:ismarm106}. Notice that even though the slowness perturbation is restricted to the region below the black horizon, the image perturbation is also back-projected to the region of correct velocity in the model space. These back-projected residuals are gradually decreased as the velocity iteratively converges to an optimal velocity that focuses the image.  
\plot{ismarm106}{width=0.8\textwidth}{Slowness perturbation from back-projected image perturbations.}

\subsection{Image-space generalized-sources domain wave-equation tomographic operator}
We will extend the wave-equation tomographic operator from the shot-profile domain to the image-space generalized-sources domain. Although we will focus on the use of image-space generalized wavefields described in Chapters \ref{chap:chap02} and \ref{chap:chap03}, with some minor modifications the derivation presented here is also valid when using data-space generalized wavefields. We use $\widetilde D$ for the image-space generalized source wavefield and $\widetilde U$ for the image-space  generalized receiver wavefield, irrespective of whether they are PERM wavefields or ISPEWs. For the image computed with these wavefields, we use $\widetilde I$. 
\par
I also use the Marmousi model to illustrate the components of the image-space generalized-sources domain wave-equation tomographic operator. The image-space generalized source and receiver gathers were computed using 12 selected reflectors from the background image of Figure \ref{fig:ismarm102}. These selected reflectors are submitted to the rotation of the subsurface offsets according to the apparent geological dip, as discussed in Chapter \ref{chap:chap02} (Figure \ref{fig:ismarm201}). Here, we use ISPEWs initiated at a spatial-sampling period of 35 SODCIGs.
\plot{ismarm201}{width=0.6\textwidth}{Selected reflectors used to model the image-space generalized source and receiver gathers.}
\par
In the areal-shot migration of image-space generalized wavefields, both source and receiver 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){\widetilde D}({\bf x},{\bf p},\omega) = 0 \\
{\widetilde D}(x,y,z=z_{min},{\bf p},\omega) =  {\widetilde{\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){\widetilde U}({\bf x},{\bf p},\omega) = 0 \\
{\widetilde U}(x,y,z=z_{min},{\bf p},\omega) = {\widetilde{\widetilde  U}}(x,y,z=z_{min},{\bf p},\omega)
 \end{array} \right., \label{aeq_wveq_bU}
\end{eqnarray}
where ${\widetilde D}({\bf x},{\bf p},\omega)$ is the image-space generalized source wavefield for a single frequency $\omega$ at image point ${\bf x}=(x,y,z)$; ${\bf p}$ is the index of the areal shot; ${\widetilde U}({\bf x},{\bf p},\omega)$ is the image-space generalized receiver wavefield for a single frequency $\omega$ at image point ${\bf x}$; and ${\widetilde{\widetilde  D}}(x,y,z=z_{min},{\bf p},\omega)$ and ${\widetilde{\widetilde  U}}(x,y,z=z_{min},{\bf p},\omega)$ are the image-space generalized source and receiver gathers synthesized with the pre-stack exploding-reflector model, with or without phase encoding, and collected at $z=z_{min}$, which denotes the top of a target zone. These image-space generalized gathers 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) source, and b) receiver wavefields.}
\par
The cross-correlation imaging condition produces the image ${\widetilde I}({\bf x},{\bf h})$ (Figure \ref{fig:ismarm203}):
\begin{eqnarray}
{\widetilde I}({\bf x},{\bf h}) = \sum_{\bf p}\sum_{\omega} {\widetilde D}^*({\bf x}-{\bf h},{\bf p},\omega) {\widetilde 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}.}

The perturbed image is derived by applying the product rule to equation \ref{aeq_ic_bI}, which gives
\begin{eqnarray}
\Delta {\widetilde I}({\bf x},{\bf h})  &=& \sum_{\bf p}\sum_{\omega} \Delta {\widetilde D}^*({\bf x}-{\bf h},{\bf p},\omega) {\widetilde U}_0({\bf x}+{\bf h},{\bf p},\omega) + \nonumber \\
& & {\widetilde D}_0^*({\bf x}-{\bf h},{\bf p},\omega) \Delta {\widetilde U} ({\bf x}+{\bf h},{\bf p},\omega), \label{aeq_ic_dI}
\end{eqnarray}
where ${\widetilde D}_0({\bf x}-{\bf h},{\bf p},\omega)$ and ${\widetilde U}_0({\bf x}+{\bf h},{\bf p},\omega)$ are the image-space generalized background source and receiver wavefields computed with the background slowness; and $\Delta {\widetilde D}({\bf x}-{\bf h},{\bf p},\omega)$ and $\Delta {\widetilde U}({\bf x}+{\bf h},{\bf p},\omega)$ are the image-space generalized perturbed source wavefield and the image-space generalized perturbed receiver wavefield, respectively. These image-space generalized perturbed wavefields are the response to a slowness perturbation. The image-space generalized perturbed source and receiver 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 {\widetilde D}({\bf x},{\bf p},\omega)= {\widetilde D}_{\scriptscriptstyle SC}({\bf x},{\bf p},\omega)
 \\
\Delta {\widetilde 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 {\widetilde U}({\bf x},{\bf p},\omega)={\widetilde U}_{\scriptscriptstyle SC}({\bf x},{\bf p},\omega)\\
\Delta {\widetilde 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{ismarm204}{width=0.6\textwidth}{Perturbed image computed with equation \ref{aeq_ic_dI}.}
\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 source and receiver wavefields, respectively, which result from the interaction of the image-space generalized background wavefields with a slowness perturbation according to
\begin{eqnarray}
{\widetilde 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})}}}{\widetilde D}_0({\bf x},{\bf p},\omega)
\label{aeq_Dsc}
\end{eqnarray}
and 
\begin{eqnarray}
{\widetilde 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})}}}{\widetilde U}_0({\bf x},{\bf p},\omega).
\label{aeq_Usc}
\end{eqnarray}
These wavefields are injected at every depth level during the recursive 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. 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.
\par
The adjoint-tomographic operator ${\bf T_{\scriptscriptstyle IS}}'$ is obtained by applying the adjoint-imaging condition to compute the image-space generalized perturbed source and receiver wavefields given by the following convolutions:
\begin{eqnarray}
\Delta {\widetilde D}({\bf x},{\bf p},\omega) &=& \sum_{\bf h} \Delta {\widetilde I}({\bf x},{\bf h}) {\widetilde U}_0({\bf x}+{\bf h},{\bf p},\omega) \nonumber \\
\Delta {\widetilde U}({\bf x},{\bf p},\omega) &=& \sum_{\bf h} \Delta {\widetilde I}({\bf x},{\bf h}) {\widetilde 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 {\widetilde s} ({\bf x}) &=& \sum_{{\bf p}}\sum_{\omega} {\widetilde D}_{\scriptscriptstyle SC}^*({\bf x},{\bf p},\omega) \Delta {\widetilde D}({\bf x},{\bf p},\omega) + \nonumber \\
& & {\widetilde U}_{\scriptscriptstyle SC}^*({\bf x},{\bf p},\omega) \Delta {\widetilde 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}. Compare it with the slowness perturbation computed with in the shot-profile domain (Figure \ref{fig:ismarm106}). The slowness perturbation computed with ISPEWs shows the correct polarity and  a general structure similar to that of the slowness perturbation computed in the shot-profile domain. The main amplitude differences are in the left part of Figure \ref{fig:ismarm206}, where only one reflector was selected to synthesize ISPEWs.
\plot{ismarm206}{width=0.8\textwidth}{Slowness perturbation from back-projected image perturbations computed with 35 ISPEWs.}

\section{Velocity optimization using image-space generalized wavefields}
ISWET is a nonlinear optimization problem in which the optimal migration velocity is determined by driving the objective function (i.e, the perturbed image) to a minimum. The negative of the gradient of the objective function provides search directions for iterative velocity updating. In DSVA, the DSO operator applied to the current background image computed with the image-space generalized wavefields yields the perturbed image, according to
\begin{eqnarray}
\Delta {\widetilde I} ({\bf x,h}) &=& \left | {\bf h} \right | {\widetilde I} ({\bf x,h}),
\label{aeq_dso}
\end{eqnarray}
where $\bf h$ is the vector of subsurface offsets. The DSO operator penalizes the focusing of reflectors at zero-subsurface offset. Considering complete illumination and infinite frequency bandwidth, energy not focused at zero-subsurface offset indicates velocity errors. DSO easily automates ISWET. 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, such as nonlinear-conjugate gradients, L-BFGS \cite[]{nocedal}. We use a nonlinear-conjugate gradient solver, for which we need to provide the value of the objective function and its gradient for the current velocity model . 
\par
As discussed in Chapter \ref{chap:chap02}, image-space generalized wavefields produce migrated images with the wavelet squared. Because of the squaring of the wavelet, images computed with these wavefields present amplitude variations stronger than those in the image computed with the original data. This effect is also present in the DSO perturbed image. When applying $\bf T'_{\scriptscriptstyle IS}$, the DSO perturbed image is convolved with the image-space generalized background wavefields to generate the image-space generalized perturbed wavefields. Subsequently, the gradient is obtained by cross-correlating the image-space generalized perturbed wavefields with the image-space generalized scattered wavefields. These two operations introduce another squaring of the wavelet. Therefore, the gradient of the objective function computed with the image-space generalized wavefields is a fourth-power version of that computed with the original wavefields. To minimize the effects of squaring the wavelet on the gradient of the objective function, we use a signed-square-root version of the initial image as the initial conditions for synthesizing image-space generalized gathers. 
\par
In addition to the amplitude variations described above, there are also amplitude variations in the migrated image caused by uneven illumination. Considering that the gradient of the objective function is computed along the wave paths, the uneven illumination will be imprinted on the gradient. An eventual velocity update using this unbalanced amplitude gradient can originate a velocity model that violates the smoothness assumption implied by the Born approximation. 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$. Other smoothing schemes could also be applied, such as smoothing along geological dips \cite[]{bob00} computed on the migrated image using the original shot-profiles.
\par
As already described in Chapter \ref{chap:chap03}, when using image-space generalized wavefields, a target-oriented strategy can be adopted if the velocity errors are believed to be limited to a sub-region of the model. In this case, the target region, where the image-space generalized wavefields are propagated, lies below the bottom of the accurate velocity region. 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. 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 ${\bf \nabla} J_{\scriptscriptstyle IS}$ is
\begin{eqnarray}
{\bf \nabla} J_{\scriptscriptstyle IS} = {\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
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 velocity model only below the black horizon. To evaluate the influence of dispersed crosstalk on ISWET results, we use two different ISPEW datasets. One dataset is modeled with a spatial-sampling period of 35 SODCIGs, herein called 35-ISPEWs. For the other, the spatial-sampling is 11 SODCIGs, herein called 11-ISPEWs. This dataset is expected to generate more crosstalk than the 35-ISPEWs dataset. Both ISPEW datasets 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.
\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. To determine the step length, 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 initial and final background images for the 11-ISPEWs and 35-ISPEWs examples are shown in Figures \ref{fig:ismarm302} and  \ref{fig:ismarm301}, respectively. In both figures, at the top is the initial image, and at the bottom is the final image. Below the zero subsurface-offset section are shown  subsurface-offset gathers. In both figures, the final image is more focused, and the pulled-up reflectors at x = 6000 m are better positioned.
\plot{ismarm302}{width=0.6\textwidth}{a) Initial and b) final background image for the optimization with 11 ISPEWs.}
\plot{ismarm301}{width=0.6\textwidth}{a) Initial and b) final background image for the optimization with 35 ISPEWs.}
\par
The cross-plot of Figure \ref{fig:ismarm303} presents the evolution of the objective function, normalized by the highest value, for the two cases. Iteration stopped after 13 iterations in the 11-ISPEWs case and 9 iterations for the 35-ISPEWs case, because the variation of the objective function was less than the predefined value of 0.002$\%$. Overall, both cases present similar convergence. However, the final value of the 11-ISPEWs objective function is slightly greater than that of the 35-ISPEWs case. This can be explained by the different amount of dispersed crosstalk in the images of  Figures \ref{fig:ismarm302} and \ref{fig:ismarm301}. The 11-ISPEWs case has more dispersed crosstalk than the 35-ISPEWs case, which contributes to the higher value of the objective function.
\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
The final velocity models are shown in Figure \ref{fig:ismarm304}. In both cases, 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. 
\plot{ismarm304}{width=0.8\textwidth}{Optimized velocity models for: a) the 11-ISPEWs case, and b) 35-ISPEWs case.}
\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 models (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 models of Figure \ref{fig:ismarm304}: a) the 11-ISPEWs case, and b) 35-ISPEWs case.}

\plot{ismarm306}{width=0.6\textwidth}{Image computed with shot-profile migration using the true velocity model.}
\par
Although images computed with 11-ISPEWs dataset have more dispersed crosstalk than that computed with 35-ISPEWs dataset, ISWET results are insignificantly affected. ISWET using either the 11-ISPEWs set, or the 35-ISPEWs set 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 the 11-ISPEWs case would be nearly 30 times faster than ISWET in the shot-profile domain, and the 35-ISPEWs case would be 15 times faster.
\par
\section{conclusions}
This chapter presents ISWET in the image-space generalized-sources domain. It is an extension from the shot-profile domain to this new generalized-sources domain in which a dramatic gain in computational efficiency is achieved by decreasing data size and solving ISWET in a target-oriented manner. Also, by selecting key reflectors to initiate the image-space generalized source and receiver gathers, great flexibility is incorporated, since a horizon-based approach is possible, which can improve convergence.
\par
For the Marmousi velocity model, two different ISPEW datasets were used to optimize the migration-velocity model by ISWET. Using either dataset yields an equally accurate velocity model, which indicates robustness of the velocity inversion in the image-space generalized-sources domain.

\section {acknowledgments}
I would like to acknowledge Yaxun Tang for the fruitful discussions about ISWET and with whom I collaborated on extending ISWET from the shot-profile domain to the generalized-sources domain \cite[]{yax2008}.
