\chapter{MVA using image-space generalized sources}
\label{chap:chap04}
The theory of migration-velocity analysis by wavefield extrapolation has already been developed using survey-sinking migration and also using shot-profile migration. In this chapter, I extend the theory of migration-velocity analysis by wavefield extrapolation using areal-shot migration, more specifically using wavefields computed in the image-space generalized-source domain. In this new domain, PERM wavefields and ISPEW are used to compute the velocity updates. These wavefields allow fast velocity updates compared to when conventional wavefields are used, like in  shot-profile migration. The small data size of image-space generalized-source wavefields in addition to their inherent capability for being used in a target-oriented manner are the main causes for the computational efficiency. Moreover, the use of these wavefields enable us to incorporate well-established strategies in migration-velocity analysis by ray-methods, such as horizon-based tomography and possibly grid-based tomography, into migration-velocity analysis by wavefield extrapolation. This new feature gives more flexibility to migration-velocity analysis 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.

\section{introduction}

Wave-equation tomography aims at solving for earth models that explain observed seismograms under some norm. There are two main categories, depending on the domain in which the objective function is computed. In one category, known as waveform inversion, the objective function is minimized in the data space where data modeled with the current model parameters data are compared with the recorded seismograms. 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 incomplete acquisition, band limited nature of seismic data, presence of noise and incomplete physics of the operators. Because of non-linearity and ill-posedness the initial model plays a crucial role in waveform inversion. Ideally, it should well describe the lower-frequency components of the velocity model, whose higher-frequency components are to be determined. Extensive research has been being devoted to waveform inversion applied to exploration seismic \cite[]{lailly83, tarantola1984, tarantola1987, mora:1211, woodward1992, pratt98}. A comprehensive overview of waveform inversion can be found in \cite {virieux:WCC1}. In spite of its maturity as a technology, only recently with the advances of high-performance computing 3D-field data results have been published \cite[]{vigh:VE135,plessix:WCC149}.
\par
In the other category, here called image-space wave-equation tomography (ISWET), the objective function is minimized 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 image perturbation is a function of the model parameters, here represented by the migration velocity. Because under the Born approximation the migration velocity ISWET optimizes is smooth by definition, the ill-posedness of ISWET is less severe than that of waveform inversion. An essential feature of ISWET is that its objective function is intimately related to the final product of seismic processing, which is the image of the subsurface. The optimized velocity is the one which gives the best image.
\par
 Two variants of ISWET are wave-equation migration velocity analysis (WEMVA) \cite[]{sava2004a, sava2004b} and differential semblance velocity analysis (DSVA) \cite[]{shen2004, shen:VE49}. 
\par
Despite their differences, those two categories are linked by the concept of extended modeling \cite[]{symes2008}. \cite{symes2008} showed that ISWET can be regarded as a solution method for the partially linearized waveform inversion. Whereas ISWET solves for the long-wavelength components of the slowness model, assuring global convergence, waveform inversion computes the short-wavelength components, providing accuracy.
\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 wavefield propagation is a linear process, generalized sources are computed by linearly combining source wavefields and receiver wavefields, using phase-encoding techniques \cite[]{whitmore, romero}. In particular, \cite{guerra2009a} used phase-encoding modeling in the image space to synthesize generalized source functions that drastically decrease the cost of DSVA. 
\par
WEMVA and DVSA seek the optimal slowness 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 applied. As \cite{biondi2008} points out, WEMVA is not easily automated. The image perturbation is computed by the linearized residual prestack-depth migration \cite[]{sava2004a}, 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 linear conjugate-gradient methods. Linear iterations drive the the objective function to a minimum for the current residual-moveout parameter field, which gives the current velocity perturbation. The corresponding velocity solution is used to compute a new image whose focusing will be evaluated by a new run of residual-moveout parameter interpretation.
\par
In DSVA, the perturbed image is computed by applying the fully automated differential-semblance operator (DSO) 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, prevents the use of linear conjugate-gradient methods and, therefore, the objective function computed with DSO is typically minimized by nonlinear optimization methods. %, which require computation of the gradient of the objective function.
\par







\section{image-space wave-equation tomography}
Image-space wave-equation tomography is a non-linear inverse problem 
that tries to find an optimal background slowness that minimizes the residual field, ${\bf \Delta I}$, defined in the image space.
The residual field is derived from the background image, ${\bf I}$, which is computed with a background slowness (or the current estimate of the slowness). 
The residual field measures the correctness of the background slowness;  
its minimum (under some norm, e.g. $\ell_2$) is achieved when a correct background slowness has been used for migration. 
There are many choices of the residual field, such as residual moveout in the Angle-Domain Common-Image Gathers (ADCIGs), 
differential semblance in the ADCIGs, 
reflection-angle stacking power (in which case we have to maximize the residual field, or minimize the negative stacking power), etc..
Here we follow a definition similar to that in \cite{biondi2008}, and define a general form of the residual field as follows:
\begin{eqnarray}
{\bf \Delta I} = {\bf I} - {\bf F}({\bf I}),
\end{eqnarray}
where ${\bf F}$ is a focusing operator, which measures the focusing of the migrated image. For example,
in the Differential Semblance Optimization (DSO) method \cite[]{shen2004}, the focusing operator takes the following form:
\begin{eqnarray}
{\bf F}({\bf I}) = \left( {\bf 1} - {\bf O} \right) {\bf I}, \label{eq_focus_dso}
\end{eqnarray}
where ${\bf 1}$ is the identity operator and ${\bf O}$ is the DSO operator either in the subsurface offset domain or in the angle domain \cite[]{shen2004}.
The subsurface-offset-domain DSO focuses the energy at zero offset, whereas the angle-domain DSO flattens the ADCIGs.

In the wave-equation migration velocity analysis (WEMVA) method \cite[]{sava2004}, the focusing operator is the linearized residual migration operator defined as follows: 
\begin{eqnarray}
{\bf F}({\bf I}) = {\bf R}[\rho]{\bf I} \approx {\bf I} + {\bf K}[\Delta \rho]{\bf I},
\end{eqnarray}
where ${\bf \rho}$ is the ratio between the background slowness $\widehat{\bf s}$ and the true slowness ${\bf s}$, 
and $\Delta \rho = 1-\rho = 1-\frac{\widehat{\bf s}}{\bf s}$; 
${\bf R}[\rho]$ is the residual migration operator \cite[]{sava2003}, and 
${\bf K}[\Delta \rho]$ is the differential residual migration operator defined as follows \cite[]{sava2004a,sava2004b}:
\begin{eqnarray}
{\bf K}[\Delta \rho] = \Delta \rho \left. \frac{\partial {\bf R}[\rho]}{\partial \rho} \right|_{\rho=1}.
\end{eqnarray}
The linear operator ${\bf K}[\Delta \rho]$ applies different phase rotations to the image for different reflection angles and geological dips \cite[]{biondi2008}.

In general, if we choose $\ell_2$ norm, the tomography objective function to minimize can be written as follows:
\begin{eqnarray}
J = \frac{1}{2}||{\bf \Delta I}||_2 = \frac{1}{2}||{\bf I}-{\bf F}({\bf I})||_2,
\end{eqnarray}
where $||\cdot||_2$ stands for the $\ell_2$ norm.
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}$ reads as follows:
\begin{eqnarray}
{\bf \nabla} J = \Re \left( \left(\frac{\partial {\bf I}}{\partial {\bf s}} - \frac{\partial {\bf F}({\bf I})}{\partial {\bf s}} \right)^{\prime} \left( {\bf I}-{\bf F}({\bf I})\right) \right), \label{eq_obj_grad_gen}
\end{eqnarray}
where $\Re$ denotes taking the real part of a complex value and $^{\prime}$ denotes the adjoint. 
For the DSO method, the linear operator ${\bf O}$ is independent of the slowness, 
so we have
\begin{eqnarray}
\frac{\partial {\bf F}({\bf I})}{\partial {\bf s}} = ({\bf 1}-{\bf O})\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} = \Re \left( \left.\left(\frac{\partial {\bf I}}{\partial {\bf s}}\right|_{{\bf s}=\widehat{\bf s}}\right)^{\prime}{\bf O}^{\prime}{\bf O}\widehat{\bf I}\right), \label{eq_grad_dso}
\end{eqnarray}
where $\widehat {\bf I}$ is the background image computed using the background slowness $\widehat{\bf s}$.

For the WEMVA method, the gradient is slightly more complicated, because in this case, the focusing operator is also dependent 
on the slowness ${\bf s}$. However, one can simplify it by assuming that the focusing operator is applied on the background
image $\widehat{\bf I}$ instead of ${\bf I}$, and $\widehat{\Delta \rho}$ is also picked from the background image $\widehat{\bf I}$, that is 
\begin{eqnarray}
{\bf F}(\widehat{\bf I})=\widehat{\bf I}+{\bf K}[\widehat{\Delta \rho}]\widehat{\bf I}.
\end{eqnarray}
With these assumptions, we get the "classic" WEMVA gradient as follows:
\begin{eqnarray}
\nabla J_{\rm WEMVA} = \Re \left( - \left.\left(\frac{\partial {\bf I}}{\partial {\bf s}}\right|_{{\bf s}=\widehat{\bf s}}\right)^{\prime}{\bf K}[\widehat{\Delta \rho}]\widehat{\bf I} \right)
\end{eqnarray}
The complete WEMVA gradient without the above assumptions can also be derived following the method described by \cite{biondi2008}. 

No matter which gradient we choose to back-project the slowness perturbation, we have to evaluate the adjoint of the linear operator
$\left.\frac{\partial{\bf I}}{\partial{\bf s}}\right|_{{\bf s}=\widehat{\bf s}}$, which defines a linear mapping from
the slowness perturbation ${\bf \Delta s}$ to the image perturbation ${\bf \Delta I}$. This is easy to see by
expanding the image ${\bf I}$ around the background slowness ${\widehat {\bf s}}$ as follows:
\begin{eqnarray}
{\bf I} = {\widehat {\bf I}} + \left. \frac{\partial{\bf I}}{{\partial{\bf s}}}\right|_{{\bf s}=\widehat{\bf s}}({\bf s}-{\widehat {\bf s}}) + \cdots,
\end{eqnarray} 
Keeping only the zero and first order terms, we get the linear operator $\left. \frac{\partial{\bf I}}{\partial{\bf s}}\right|_{{\bf s}=\widehat{\bf s}}$ as follows:
\begin{eqnarray}
{\bf \Delta I} = \left. \frac{\partial{\bf I}}{\partial{\bf s}}\right|_{{\bf s}=\widehat{\bf s}} {\bf \Delta s} = {\bf T}{\bf \Delta s},
\end{eqnarray}
where ${\bf \Delta I} = {\bf I} - {\widehat {\bf I}}$ and ${\bf \Delta s} = {\bf s} - {\widehat {\bf s}}$. 
${\bf T} = \left. \frac{\partial{\bf I}}{\partial{\bf s}}\right|_{{\bf s}=\widehat{\bf s}}$ is the wave-equation tomographic operator.
The tomographic operator can be evaluated either in the source and receiver domain \cite[]{sava2004} or in the shot-profile domain \cite[]{shen2004}.
In next section we follow an approach similar to that discussed by \cite{shen2004} and review the forward and adjoint tomographic operator 
in the shot-profile domain. In the subsequent sections, we generalize the expression of the tomographic operator to 
the generalized source domain.

\section{the tomographic operator in the shot-profile domain}

For the conventional shot-profile migration, both source and receiver wavefields are downward continued with the 
following one-way wave equations \cite[]{claerbout1971}:
\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) = \overline{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 the overline stands for complex conjugate; 
$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$;
$s({\bf x})$ is the slowness at ${\bf x}$; ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector; $f_s(\omega)$
is the frequency dependent source signature, and $\overline{f_s(\omega)\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}. $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}.
To produce the image, the following cross-correlation imaging condition is used:
\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}
where ${\bf h}=(h_x,h_y,h_z)$ is the subsurface half offset.

The perturbed image can be derived by a simple application of the chain rule
to Equation \ref{eq_ic_bI}: 
\begin{eqnarray}
\Delta I({\bf x},{\bf h})  &=& 
\sum_{{\bf x}_s}\sum_{\omega} 
\left( \Delta D({\bf x}-{\bf h},{\bf x}_s,\omega) {\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega) + \right. \nonumber \\
& & \left. {\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega) \Delta U ({\bf x}+{\bf h},{\bf x}_s,\omega) \right), \label{eq_ic_dI}
\end{eqnarray}
where ${\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)$ and ${\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega)$ are the background source and 
receiver wavefields computed with the background slowness ${\widehat s}({\bf x})$; 
$\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, which are the results of the slowness perturbation $\Delta s({\bf x})$.
The perturbed source and receiver wavefields satisfy the following 
one-way wave equations, which are linearized with respect to slowness (see Appendix A for derivations):
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {\widehat s}^2 ({\bf x})-|{\bf k}|^2}\right) \Delta D({\bf x},{\bf x}_s,\omega)= 
\frac{- i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}^2({\bf x})}}}{\widehat D}({\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 {\widehat s}^2 ({\bf x})-|{\bf k}|^2}\right) \Delta U({\bf x},{\bf x}_s,\omega)=
\frac{- i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}^2({\bf x})}}}{\widehat U}({\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}
Recursively solving Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} gives us the perturbed source and receiver wavefields.
The perturbed source and receiver wavefields are then used in Equation \ref{eq_ic_dI} to generate the perturbed image $\Delta I({\bf x},{\bf h})$, where
the background source and receiver wavefields are precomputed by recursively solving Equations \ref{eq_wveq_bD} and \ref{eq_wveq_bU} with a
background slowness ${\widehat s}(\bf x)$.
Appendix B gives a more detailed matrix representation of how to evaluate the forward tomographic operator ${\bf T}$.

To evaluate the adjoint tomographic operator ${\bf T}'$, we first apply the adjoint of the imaging condition in Equation \ref{eq_ic_dI} 
to get the perturbed source and receiver wavefields $\Delta D({\bf x},{\bf x}_s,\omega)$ and $\Delta U({\bf x},{\bf x}_s,\omega)$ as follows:
\begin{eqnarray}
\Delta D({\bf x},{\bf x}_s,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h})\overline{{\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega)}, \\
\Delta U({\bf x},{\bf x}_s,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h})\overline{{\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)}.
\end{eqnarray}
Then we solve the adjoint equations of Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} to get the slowness perturbation $\Delta s({\bf x})$. 
Again, in order to solve the adjoint equations of Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU}, 
the background source wavefield ${\widehat D}({\bf x},{\bf x}_s,{\omega})$ and 
the background receiver wavefield ${\widehat U}({\bf x},{\bf x}_s,\omega)$ have to be computed in advance.
Appendix C gives a more detailed matrix representation of how to evaluate the adjoint tomographic operator ${\bf T}^{'}$.

\subsection{Image-space encoded wavefields}

From Equation \ref{eq_ic_dI_me}, the perturbed image is easily obtained as follows:
\begin{eqnarray}
\Delta I_{\rm me}({\bf x},{\bf h})  &=& 
\sum_{{\bf p}_m}\sum_{\omega} 
\left( \Delta \widetilde{D}({\bf x}-{\bf h},{\bf p}_m,\omega) \widehat{\widetilde{U}}({\bf x}+{\bf h},{\bf p}_m,\omega) + \right. \nonumber \\
& & \left. \widehat{\widetilde{D}}({\bf x}-{\bf h},{\bf p}_m,\omega) \Delta \widetilde{U}({\bf x}+{\bf h},{\bf p}_m,\omega) \right),
\end{eqnarray}
where $\widehat{\widetilde{D}}({\bf x},{\bf p}_m,\omega)$ and $\widehat{\widetilde{U}}({\bf x},{\bf p}_m,\omega)$
are the image-space encoded background source and receiver wavefields;
$\Delta \widetilde{D}({\bf x},{\bf p}_m,\omega)$ and $\Delta \widetilde{U}({\bf x},{\bf p}_m,\omega)$ are the perturbed source and receiver wavefields
in the image-space phase-encoding domain, which satisfy the perturbed one-way wave equations defined by Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU}.
The tomographic operator ${\bf T}$ and its adjoint ${\bf T}^{\prime}$ can be implemented in a manner similar to that discussed in Appendices B and C, by
replacing the original wavefields with the image-space phase-encoded wavefields.



\section{introduction}
Wave-equation tomography aims to solve for earth models that explain observed seismograms under some norm. There are two main categories, depending on the domain in which the objective function is computed. In one category, known as waveform inversion \cite[]{lines84,tarantola1987,woodward1992}, the objective function is defined in the data space, where the modeled data are compared with the recorded seismograms. In the other category, here called image-space wave-equation tomography (ISWET), the objective function is minimized in the image space. Two variants of ISWET are wave-equation migration velocity analysis (WEMVA) \cite[]{sava2004a, sava2004b} and differential semblance velocity analysis (DSVA) \cite[]{shen2004, shen:VE49}. 
%\par
%Despite their differences, those two categories are linked by the concept of extended modeling \cite[]{symes2008}. \cite{symes2008} showed that ISWET can be regarded as a solution method for the partially linearized waveform inversion. While ISWET solves for the long-wavelength components of the slowness model, assuring global convergence, waveform inversion computes the short-wavelength components, providing accuracy.
\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 wavefield propagation is a linear process, generalized sources are computed by linearly combining source wavefields and receiver wavefields, using phase-encoding techniques \cite[]{whitmore1995, romero2000}. In particular, \cite{guerra2009a} used phase-encoding modeling in the image space to synthesize generalized source functions that drastically decrease the cost of DSVA. 
\par
WEMVA and DVSA seek the optimal slowness 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 applied. As \cite{biondi2008} points out, WEMVA is not easily automated. The image perturbation is computed by the linearized residual prestack-depth migration \cite[]{sava2004a}, which uses a manually picked residual-moveout parameter. Because the perturbed image computed with the linearized residual prestack-depth migration is consistent with the application of the forward wave-equation tomographic operator, the forward and adjoint ISWET operators can be used in conjugate-gradient methods to invert for the slowness perturbation. In DSVA, the perturbed image is computed by applying the fully automated differential-semblance operator (DSO) to the subsurface-offset gathers (ODCIGs) or angle gathers (ADCIGs). When applied to ODCIGs, DSVA minimizes the energy not focused at zero-offset. When applied to ADCIGs, DSVA minimizes energy departing from flatness of the reflectors. Although DSVA easily automates ISWET, it produces perturbed images that do not present the depth phase-shift introduced by the forward one-way ISWET operator, neither the DVSA amplitude behavior can be modeled by this operator. Therefore, the objective function computed with DSO is tipically minimized by quasi-Newton algorithms, which require computation of the gradient of the objective function.
\par
The gradient of the objective function can be computed with the Frech$\acute {\rm e}$t derivatives. However, even for 2D applications of ISWET this computation can be very expensive. An efficient way to compute the gradient without using Frech$\acute {\rm e}$t derivatives is the adjoint-state method \cite[]{chavent:476,plessix2006}. \cite{plessix2006} describes two methodologies for computing the gradient of the objective function using the adjoint-state method. One methodology uses the perturbation theory, which states that, at first order, a perturbation of the model parameters causes a perturbation of the objective function. The other uses the augmented Lagrangian. The augmented Lagrangian is formed by the objective function and the scalar product of the adjoint-state variables with general solutions of the forward modeling equations. The adjoint-state variables are, in turn,  solutions of the adjoint-state equations. The adjoint-state equations are defined by equating to zero the derivatives of the augmented Lagrangian with respect to the state variables. For the linear case, the adjoint of the modeling operator applied to the adjoint-state variables gives the gradient of the objective function.
\par 
Although previous studies have computed the gradient of the ISWET objective function by the adjoint-state method (for example, \cite{shen:2132}), they have not provided a detailed derivation. Here, I show a detailed derivation the gradient of the ISWET objective function using the augmented Lagrangian methodology. The derivation is valid whether ISWET uses areal-shot migration or shot-profile migration. A complete description of the forward and adjoint of ISWET operators is given in \cite{yax2008}. A numerical example of slowness optimization on the Marmousi model illustrates the use of the gradient by a quasi-Newton algorithm.
\par
Because I use image-space phase-encoded gathers in the numerical example, for completeness, I first briefly describe how to compute these phase-encoded gathers. A detailed treatment of the image-space phase-encoded gathers can be found in \cite{biondi2006, biondi2007} and \cite{guerra2008b}. Then I derive the ISWET gradient using the adjoint-state method, and show the numerical example .

\section {Image-space phase-encoded gathers}
\cite{biondi2006,biondi2007} introduced the concept of the prestack exploding-reflector as a generalization of the exploding-reflector method \cite[]{loewenthal}. The prestack exploding-reflector modeling synthesizes areal data and the corresponding areal source function, having as an initial condition a prestack image computed with wave-equation migration. If the slowness is accurate and the energy is focused at zero-subsurface offset, the prestack exploding-reflector modeling reduces to the conventional exploding-reflector method. Basically, the prestack exploding-reflector method models one single reflection event from one single ODCIG by recursive upward continuation with the following one-way wave equations:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\right) d({\bf x},\omega;x_m,y_m) = r_D({\bf x},{\bf h};x_m,y_m) \\
d(x,y,z=z_{\rm max},\omega;x_m,y_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 {\hat s}^2({\bf x})-|{\bf k}| ^2}\right) u({\bf x},\omega;x_m,y_m) = r_U({\bf x},{\bf h};x_m,y_m) \\
u(x,y,z=z_{\rm max},\omega;x_m,y_m) = 0 \end{array} \right.,
\label{eq_perm_bU}
\end{eqnarray}
where $r_D({\bf x},{\bf h};x_m,y_m)$ and $r_U({\bf x},{\bf h};x_m,y_m)$ are a single ODCIG at the horizontal location $(x_m,y_m)$ with a single reflector, and are suitable initial conditions for modeling the source and receiver wavefields, respectively. They are obtained by rotating the original unfocused ODCIGs according to the apparent geological dip of the reflector. This operation maintains the velocity information needed for migration velocity analysis, especially for dipping reflectors \cite[]{biondi2007}. Note that $d(x,y,z=z_c,\omega;x_m,y_m)$ is the areal source data and $u(x,y,z=z_c,\omega;x_m,y_m)$ is the areal receiver data for a single reflector and a single ODCIG located at $(x_m,y_m)$; $z=z_c$ denotes that the wavefields can be collected at any depth level, $z_c$. This characteristic is important for accelerating ISWET, especially if $z_c$ separates regions of sufficiently accurate slowness above and inaccurate slowness below. Therefore, the synthesized gathers are naturally datumized, and the wavefield propagations during ISWET can be restricted to the region where the slowness model must be updated. This feature allows the application of ISWET to be target-oriented.
\par
As initially formulated, if one models a single reflection from one ODCIG at a time, the prestack exploding-reflector method generates a dataset that can be orders of magnitude bigger than the original dataset. As discussed by \cite{biondi2006} and \cite{guerra2008a,guerra2008b}, modeling several reflectors and several ODCIGs simultaneously and using random phase encoding generates a much smaller dataset that, when migrated, do not produce crosstalk. The randomly encoded areal source and areal receiver wavefields can be computed as follows:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\right) {\widetilde d}({\bf x},{\bf p}_m,\omega) 
= \widetilde{r}_D({\bf x},{\bf h},{\bf p}_m,\omega) \\
{\widetilde d}(x,y,z=z_{\rm max},{\bf p}_m,\omega) = 0 \end{array} \right.,
\label{eq_perm_bD2}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {\hat s}^2({\bf x})-|{\bf k}| ^2}\right) {\widetilde u}({\bf x},{\bf p}_m,\omega) 
= \widetilde{r}_U({\bf x},{\bf h},{\bf p}_m,\omega) \\
{\widetilde u}(x,y,z=z_{\rm max},{\bf p}_m,\omega) = 0 \end{array} \right.,
\label{eq_perm_bU2}
\end{eqnarray}
where $\widetilde{r}_D({\bf x},{\bf h},{\bf p}_m,\omega)$ and $\widetilde{r}_U({\bf x},{\bf h},{\bf p}_m,\omega)$ are the encoded ODCIGs after rotations. They are defined as follows:
\begin{eqnarray}
\widetilde{r}_D({\bf x},{\bf h},{\bf p}_m,\omega) &=& \sum_{x_m}\sum_{y_m}r_D({\bf x},{\bf h},x_m,y_m)\beta({\bf x},x_m,y_m,{\bf p}_m,\omega), \\
\widetilde{r}_U({\bf x},{\bf h},{\bf p}_m,\omega) &=& \sum_{x_m}\sum_{y_m}r_U({\bf x},{\bf h},x_m,y_m)\beta({\bf x},x_m,y_m,{\bf p}_m,\omega),
\end{eqnarray}
where $\beta({\bf x},x_m,y_m,{\bf p}_m,\omega)=e^{i\gamma({\bf x},x_m,y_m,{\bf p}_m,\omega)}$ is a pseudo-random phase-encoding function, with $\gamma({\bf x},x_m,y_m,{\bf p}_m,\omega)$ being a uniformly distributed random sequence in ${\bf x}$, $x_m$, $y_m$ and $\omega$; the variable ${\bf p}_m$ is the index of different realizations of the random sequence. The recursive solution of equations \ref{eq_perm_bD2} and \ref{eq_perm_bU2} gives the encoded areal source data ${\widetilde d}(x,y,z=z_c,{\bf p}_m,\omega)$ and areal receiver data ${\widetilde u}(x,y,z=z_c,{\bf p}_m,\omega)$, at the depth level, $z_c$. 

\section{Gradient of ISWET by the adjoint-state method}
The gradient of the objective function can be computed with the Frech$\acute {\rm e}$t derivatives. However, even for 2D applications of ISWET this computation can be very expensive. An efficient way to compute the gradient without using Frech$\acute {\rm e}$t derivatives is the adjoint-state method \cite[]{chavent:476,plessix2006}. \cite{plessix2006} describes two methodologies for computing the gradient of the objective function using the adjoint-state method. One methodology uses the perturbation theory, which states that, at first order, a perturbation of the model parameters causes a perturbation of the objective function. The other uses the augmented Lagrangian. The augmented Lagrangian is formed by the objective function and the scalar product of the adjoint-state variables with general solutions of the forward modeling equations. The adjoint-state variables are, in turn,  solutions of the adjoint-state equations. The adjoint-state equations are defined by equating to zero the derivatives of the augmented Lagrangian with respect to the state variables. For the linear case, the adjoint of the modeling operator applied to the adjoint-state variables gives the gradient of the objective function.
\par 
Although previous studies have computed the gradient of the ISWET objective function by the adjoint-state method (for example, \cite{shen:2132}), they have not provided a detailed derivation. Here, I show a detailed derivation the gradient of the ISWET objective function using the augmented Lagrangian methodology. The derivation is valid whether ISWET uses areal-shot migration or shot-profile migration. A complete description of the forward and adjoint of ISWET operators is given in \cite{yax2008}. A numerical example of slowness optimization on the Marmousi model illustrates the use of the gradient by a quasi-Newton algorithm.
\par



Image-space wave-equation tomography aims to iteratively solve for the slowness model, $s=s({\bf x})$, that minimizes the nonlinear objective function:
\begin{eqnarray}
J(s)=\frac{1}{2} \left | \left | \Delta r(s) \right | \right | ^2 = \frac{1}{2} \left | \left | r(s) - {\bf M} r(s) \right | \right | ^2,
\label{eq:objfun}
\end{eqnarray}
where $\Delta r = \Delta r({\bf x},{\bf h})$ is the image perturbation that measures the accuracy of the slowness model. To compute $\Delta r$, a differential residual-focusing operator $\bf M$ is applied to the image $r = r({\bf x},{\bf h})$ obtained with the current slowness \cite[]{biondi2008}, using either differential residual prestack migration \cite[]{sava2004a,sava2004b} or differential-semblance optimization (DSO) operators \cite[]{shen:VE49}. In this paper, operators are represented by bold capital letters.
\par
If the differential residual-focusing operator $\bf M$ is independent of the slowness, the gradient of this objective function evaluated at the current slowness ${\widehat s}={\widehat s}(\bf x)$ is
\begin{eqnarray}
\nabla J(s)= \left. \left ( \frac{\partial r}{\partial s} \right )'\right |_{s=\widehat s} \left ({\bf I - M}' \right )  {\Delta \widehat r},
\label{eq:gradobjfun}
\end{eqnarray}
where the prime denotes the adjoint, $\bf I$ is the identity operator, and $\Delta \widehat r=\Delta \widehat r({\bf x},{\bf h})$ is the perturbed image obtained with the current slowness model.  The linear operator $\frac{\partial r}{\partial s}$ defines the mapping $\frac{\partial r}{\partial s} \Delta s = \Delta r$ between the slowness perturbation $\Delta s$ and the image perturbation $\Delta r$; it is called the image-space wave-equation tomographic operator.
\par
Because the image-space wave-equation tomographic operator is composed of different operators, it is difficult to envision from equation \ref{eq:gradobjfun} which operations are performed to compute the gradient. Therefore, for a clear explanation of the operators involved, I use the adjoint-state method to derive the gradient of the objective function (equation \ref{eq:objfun}).
\par
In migration with generalized sources or shot-profile migration, the source and receiver wavefields are propagated independently, and the image $r_z = r_z({\bf x},{\bf h})$ at a depth level $z$, is computed by the crosscorrelation
\begin{eqnarray}
r_z({\bf x},{\bf h}) = \sum_{\omega} d_z^*({\bf x}-{\bf h},\omega) u_z({\bf x}+{\bf h},\omega), 
\label{eq:imgcond1}
\end{eqnarray}
where $d_z({\bf x},\omega)$ is the source wavefield for a single frequency $\omega$ at horizontal coordinates ${\bf x}=(x,y)$ ; $u_z({\bf x},\omega)$ is the receiver wavefield, ${\bf h}=(h_x,h_y)$ is the subsurface half-offset, and ``$^*$'' stands for the complex-conjugate. An additional summation over shots is required when migrating more than one shot. Hereafter, letters $d$ and $u$ stand for source and receiver wavefields, respectively, irrespective of the migration scheme.
\par
In a more compact notation, not explicitly writing the dependencies on $\bf x$ and $\bf h$, equation \ref{eq:imgcond1} can be re-written as follows:
\begin{eqnarray}
r_z = {\bf{S}}{\bf{D}}_z'(\omega) u_z(\omega) = {\bf{S}}{\bf{U}}_z(\omega) d_z^*(\omega), 
\label{eq:imgconcompact}
\end{eqnarray}
where ${\bf{D}}$ and ${\bf{U}}$ are convolutional operators composed of $(h_x,h_y)$-shifted versions of $d_z({\bf x},\omega)$ and $u_z({\bf x},\omega)$, respectively. Operator ${\bf{S}}$ corresponds to the summation over frequency.
\par
For subsequent depth levels, $d({\bf x},\omega)$ is computed by means of the recursive downward propagation 
\begin{eqnarray}
\left\{ \begin{array}{l}
d_{z+1}(\omega)  = {\bf T}_z^{\downarrow}(\omega,s)d_z(\omega) \\ 
d_1(\omega) = q(\omega),
\label{eq:downmod}
\end{array}
\right.
\end{eqnarray}
where ${\bf T}_z^{\downarrow}$ is the downward continuation operator, which is a function of the slowness $s$, and $q(\omega)$ is the source wavefield used as a boundary condition. In the case of conventional shot-profile migration, $q(\omega) = f_s(\omega)\delta({\bf x - x_s})$ is the source signature located at ${\bf x}_s=(x_s,y_s,0)$. If using the generalized sources in the image space, $q(\omega)$ represents  the image-space phase-encoded source wavefield of equation~\ref{eq_perm_bD2}.
\par
The downward continuation of the receiver wavefield is performed by
\begin{eqnarray}
\left\{ \begin{array}{l}
u_{z+1}(\omega)  = {\bf T}_z^{\downarrow}(\omega,s)u_z(\omega) \\ 
u_1(\omega) = w(\omega),
\label{eq:upmod01}
\end{array}
\right.
\end{eqnarray}
where $w(\omega)$ is the recorded data at the surface for shot-profile migration. If using generalized sources in the image space, $w(\omega)$ is the phase-encoded areal receiver wavefield of equation~\ref{eq_perm_bU2}. In equations \ref{eq:downmod} and \ref{eq:upmod01}, I omitted the dependencies of the wavefield with respect to ${\bf x}$. The subscript ${\it 1}$ in equations \ref{eq:downmod} and \ref{eq:upmod01} represents the surface for the shot-profile migration and the ``collection'' depth level, $z_c$, for the image-space phase-encoded wavefields.
\par
In the image-space wave-equation tomography problem, the perturbed source and receiver wavefields, and the image perturbation are used to compute the slowness perturbation that updates the current slowness model. From the perturbation theory, we have $d = \widehat d + \Delta d$, $u = \widehat u + \Delta u$, and, consequently, $r = \widehat r + \Delta r$ are physical realizations with $s = \widehat s + \Delta s$, where the $hat$ refers to fields obtained with the background slowness. To the first order (Born approximation), these perturbed fields are given by 
\begin{eqnarray}
\Delta d_{z+1}(\omega) = {\bf T}_z^{\downarrow}(\omega,\widehat s) \Delta d_z(\omega) + \widetilde {\widetilde{\bf D}}_z(\omega) \Delta s_z %\Delta {\bf T}_z^{\downarrow}(\omega,\widehat  s) {\widehat d}_z(\omega),
\label{eq:pertsou}
\end{eqnarray}
and
\begin{eqnarray}
\Delta u_{z+1}(\omega) = {\bf T}_z^{\downarrow}(\omega,\widehat s) \Delta u_z(\omega) + \widetilde {\widetilde{\bf U}}_z(\omega) \Delta s_z. %\Delta {\bf T}_z^{\downarrow}(\omega,\widehat s) {\widehat u}_z(\omega)
\label{eq:pertrec}
\end{eqnarray}
The diagonal operators $\widetilde {\widetilde{\bf D}}_z$ and $\widetilde {\widetilde{\bf U}}_z$ have in the diagonal entries the scattered source and receiver wavefields, respectively. These wavefields are given by the action of the scattering operator $\Delta {\bf T}_z^{\downarrow}$ on the background wavefields:
\begin{eqnarray}
\widetilde {\widetilde {\bf D}}_z(\omega) = \Delta {\bf T}_z^{\downarrow}(\omega,\widehat s) {\widehat d}_z(\omega) = i\frac{\omega^2 \widehat s}{\sqrt{\omega^2 \widehat s^2 - \left | {\bf k} \right | ^2}} dz~ \widehat d_z(\omega)
\label{eq:pscatt}
\end{eqnarray}
and
\begin{eqnarray}
\widetilde {\widetilde {\bf U}}_z(\omega) = \Delta {\bf T}_z^{\downarrow}(\omega,\widehat s) {\widehat u}_z(\omega) = -i\frac{\omega^2 \widehat s}{\sqrt{\omega^2 \widehat s^2 - \left | {\bf k} \right | ^2}} dz~ \widehat u_z(\omega).
\label{eq:uscatt}
\end{eqnarray}
\par
The perturbed image is given by
\begin{eqnarray}
\Delta r_{z} = {\bf S} \left ({\widehat {\bf U}}_z(\omega) \Delta d_z^*(\omega) + {\widehat {\bf D}}_z'(\omega) \Delta u_z(\omega)  \right ).
\label{eq:pertimg}
\end{eqnarray}
The matrix representations of equations \ref{eq:pertsou}, \ref{eq:pertrec}, \ref{eq:pertimg} are
\begin{eqnarray}
\Delta {\underline {\bf d}} = {\bf T}^{\downarrow} \Delta {\underline {\bf d}} + {\widetilde {\widetilde {\bf P}}} {\bf S}' \Delta {\bf s},
\label{eq:pertsou2}
\end{eqnarray}
\begin{eqnarray}
\Delta {\underline {\bf u}} = {\bf T}^{\downarrow} \Delta {\underline {\bf u}} + {\widetilde {\widetilde{\bf U}}} {\bf S}' \Delta {\bf s},
\label{eq:pertrec2}
\end{eqnarray}
and
\begin{eqnarray}
\Delta {\bf {\underline r}} = {\bf S} \left ({\widehat {\bf U}} \Delta {\underline {\bf d}}^* + {\widehat {\bf D}}' \Delta {\underline {\bf u}}  \right ),
\label{eq:pertimg2}
\end{eqnarray}
where ${\bf S}'$ is a spreading operator that replicates the slowness perturbation for every frequency.
\par
Equations \ref{eq:pertsou2}, \ref{eq:pertrec2}, and \ref{eq:pertimg2} are the forward modeling equations of the image-space wave-equation tomography problem using the generalized sources or shot-profile  schemes. They depend on the state variables $\Delta {\underline {\bf d}}$, $\Delta {\underline {\bf u}}$, and $\Delta {\underline {\bf r}}$. \cite{plessix2006} describes how to compute the adjoint states using the augmented functional methodology. By introducing the adjoint-state variables ${\underline \lambda}_d$, ${\underline \lambda}_u$, and ${\underline \lambda}_r$, the augmented Lagrangian reads
\begin{eqnarray}
\begin{array}{l}
\mathcal{L}(\Delta {\underline {\bf d}}, \Delta {\underline {\bf u}}, \Delta {\underline {\bf r}}, {\underline \lambda}_d, {\underline \lambda}_u, {\underline \lambda}_r; \Delta {\bf s}) =  \mathcal{R} \left [ \frac{1}{2} \left | \left | \Delta {\underline {\bf r}} \right | \right | ^2 - \right. \nonumber \\
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \left \langle {\underline \lambda}_d,\left ( {\bf I} - {\bf T}^\downarrow \right )    \Delta {\underline {\bf d}} - {\widetilde {\widetilde{\bf D}}} {\bf S}' \Delta {\bf s} \right \rangle - \nonumber \\
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \left \langle {\underline \lambda}_u,\left ( {\bf I} - {\bf T}^{\downarrow} \right )    \Delta {\underline {\bf u}} - {\widetilde {\widetilde{\bf U}}} {\bf S}' \Delta {\bf s} \right \rangle - \nonumber \\
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \left \langle {\underline \lambda}_r,\Delta {\underline {\bf r}} - {\bf S} \left ({\widehat {\bf U}} \Delta {\underline {\bf d}}^* + {\widehat {\bf D}}' \Delta {\underline {\bf u}} \right )  \right  \rangle. \\
\end{array}
\label{eq:augfct01}
\end{eqnarray}
The adjoint-state variables are computed by taking the derivative of $\mathcal{L}$ with respect to the state variables and equating to zero, which gives
\begin{subequations} 
\label{subeqn:list} 
\begin{align} 
\left ( {\bf I} - {\bf T}^\downarrow \right )' {\underline \lambda}_d = {\widehat {\bf U}} {\underline \lambda}_r, \label{subeqn:adjstvar1} \\
\left ( {\bf I} - {\bf T}^\downarrow \right )' {\underline \lambda}_u = {\widehat {\bf D}} {\underline \lambda}_r, \label{subeqn:adjstvar2} \\
{\underline \lambda}_r = \Delta {\bf {\underline r}}.\label{subeqn:adjstvar3} 
\end{align} 
\end{subequations} 
Notice that 
\begin{eqnarray}
\left ( {\bf I} - {\bf T}^\downarrow \right )' = \left ( {\bf I} - {\bf T}^{\downarrow '} \right ) = \left ( {\bf I} - {\bf T}^{\uparrow} \right )
\end{eqnarray}
corresponds to the recursive upward propagation operator. Therefore, equations \ref{subeqn:adjstvar1} and \ref{subeqn:adjstvar2} can be written as
\begin{subequations} 
\label{subeqn:list} 
\begin{align} 
{\underline \lambda}_p = {\bf T}^\uparrow {\underline \lambda}_d + {\widehat {\bf U}} {\underline \lambda}_r, \label{subeqn:adjstvar10} \\
{\underline \lambda}_u = {\bf T}^\uparrow {\underline \lambda}_u + {\widehat {\bf D}} {\underline \lambda}_r, \label{subeqn:adjstvar20} 
\end{align} 
\end{subequations} 
which correspond to the recursive upward propagation of the perturbed wavefields resulting from the convolution of the wavefields computed with the current slowness and the perturbed image.
\par
Finally, the gradient of $J$ is 
\begin{eqnarray}
\nabla_s J({\bf s}) = {\bf S} \left ({\widetilde {\widetilde{\bf D}}}'{\underline \lambda}_d  +  {\widetilde {\widetilde{\bf U}}}'{\underline \lambda}_u \right ).
\label{eq:adjstgrad}
\end{eqnarray}
To compute the gradient, the adjoint-state wavefields, ${\underline \lambda}_d$ and ${\underline \lambda}_u$, are upward propagated and cross-correlated in time with the scattered wavefields.

