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

\usepackage{amsmath}

\setfigdir{Fig}

\begin{document}
%\documentclass[manuscript,revised]{geophysics}
%\begin{document}

\lefthead{Guerra}
\righthead{ISWET Gradient}
%\email{\footnotemark[1]guerrac@stanford.edu \footnotemark[2]tang@stanford.edu \footnotemark[3]biondo@stanford.edu}
\footer{SEP--139}
\title{Gradient of image-space wave-equation tomography by the adjoint-state method }
\author{Claudio Guerra}

\maketitle

\begin{abstract}
Optimization with gradient-descent techniques requires computing the gradient of the objective function. The gradient can be determined by using the Frech$\acute {\rm e}$t derivatives, which, for practical problems, can be very expensive. The gradient can be efficiently computed by using the adjoint-state method, which does not require the use of the Frech$\acute {\rm e}$t derivatives. Here, I derive the gradient of the image-space wave-equation tomography using the adjoint-state method. I also show its application with a numerical example using image-space phase-encoded gathers.
\end{abstract}

\section{introduction}
Wave-equation tomography aims to solve for Earth models that explain the observed seismograms under some norm. Two main categories exist, 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, called here 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. The numerical example shown in this report also uses this kind of generalized sources.
\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. Notwithstanding the ability of automating ISWET, DSVA perturbed images do not contain the depth phase-shift introduced by the forward one-way ISWET operator, neither the DVSA amplitude behavior is modeled by this operator. Therefore, to minimize the objective function computed with DSO, quasi-Newton algorithms are used, which requires computation of the gradient of the objective function.
\par
The gradient of the objective function can be obtained by computing 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 formulation uses the perturbation theory and 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, in turn, are 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 
While in other papers the gradient of the ISWET objective function is computed by the adjoint-state method (for example, \cite{shen:2132}), they do not provide a detailed derivation though. Here, I show the 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}. 
\par
Because I use image-space phase-encoded gathers in the numerical example, for completeness, I 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 a numerical example of DSVA on the Marmousi model.

\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} that synthesizes areal data and the corresponding areal source function, having as initial condition a prestack image computed with wave-equation migration. If the slowness is accurate and the energy is focused at zero-subsurface offset, it 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)$ for a single reflector, and are suitable for the initial conditions for 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 to accelerate ISWET, especially if $z_c$ separates regions of sufficiently accurate slowness above and inaccurate slowness below. As a consequence, 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 ISWET to be applied in a target-oriented manner.
\par
As initially formulated, however, 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}
Image-space wave-equation tomography aims to iteratively solve for the slowness model, $s=s({\bf x})$, that minimizes the 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. Notice that the objective function depends on the slowness model through the image perturbation. 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}. Herein, 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 image-space wave-equation tomographic operator.
\par
Because the image-space wave-equation tomographic operator is composed of different operators, it is difficult from equation \ref{eq:gradobjfun} to envision 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 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 image perturbations are used to compute the slowness perturbation to update the current slowness model. From the perturbation theory, we have that $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
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.

\section{Numerical example}
The gradient of the ISWET objective function computed in the previous section can be used in a quasi-Newton optimization scheme. I use the L-BFGS-B bound constrained optimization algorithm \cite[]{nocedal} to invert for slowness using the Marmousi model. A B-spline smoothing, with nodes spaced at 240 $m$ in $x$ and 16 $m$ in $z$, is applied to the gradient to prevent problems caused by poor  illumination and noise. I model 375 two-way shots with 24 $m$ spacing using the original Marmousi model of Figure \ref{fig:true}a. The maximum offset is 6600 $m$.
\par
The one-way shot profile image with the correct slowness model can be seen in Figure \ref{fig:true}b. The modeled data were also migrated using the background slowness of Figure \ref{fig:bckg}a to compute the background image of Figure \ref{fig:bckg}b. The background slowness differs from the true slowness in a limited region, as can be seen in Figure \ref{fig:dslow}, which shows the ratio between true and background slowness. When comparing the two images, there is a clear pull-up effect caused by migrating with a too slow velocity at the center-bottom of the background image.
\plot{true}{width=\columnwidth,height=0.8\textheight}{(a) True slowness of the Marmousi model. (b) Zero-subsurface-offset section of the one-way shot-profile migration with the true slowness model.[{\bf CR}]}
\plot{bckg}{width=\columnwidth,height=0.8\textheight}{(a) Background slowness. (b) Zero-subsurface-offset section of the one-way shot-profile migration with the  background slowness model.[{\bf CR}]}
\plot{dslow}{width=\columnwidth,height=0.4\textheight}{ Ratio between the true and the background slowness.[{\bf ER}]}
\par
Following \cite{guerra2009a}, representative reflectors were selected from the background image (Figure \ref{fig:bwkim}) to synthesize eleven image-space phase-encoded source and receiver wavefields using the background slowness. Prior to modeling, the selected reflectors were subjected to rotation according to the apparent geological dip as mentioned before. The use of only eleven areal shots decreases the computational cost of the slowness optimization in $95 \%$ if it was carried out using the original 376 original shot-profiles.
\plot{bwkim}{width=\columnwidth,height=0.4\textheight}{Zero-subsurface-offset section showing the reflectors selected on the one-way shot-profile migration with the  background slowness model.[{\bf CR}]}
\par
The optimization stopped after 4 iterations comprising 41 function and gradient evaluations. The optimized slowness model can be seen in Figure \ref{fig:opti}a. The optimized slowness decreased as expected, but did not recover the details of the true slowness. The smoothness of the optimized slowness is partially because of the B-spline filtering of the gradient and also because image-space wave-equation tomography cannot solve for the short wavelengths of the slowness model, as pointed out by \cite{symes2008}.
\par
The accuracy of the optimized slowness model can be evaluated by comparing the images computed with true, background and optimized slownesses. The shot-profile migration of the original shots using the optimized slowness is shown in Figure \ref{fig:opti}b. Notice that the pull-up effect, which occurs in the background image, is greatly mitigated. In addition, reflectors are more focused in the image computed with the optimized slowness than in the image computed with the background slowness. Another slowness accuracy indicator is flatness of the reflectors in the ADCIGs. The ADCIGs resulting from shot-profile migration of the original shots using, from top to bottom, the true slowness, the background slowness and the optimized slowness are shown in Figure \ref{fig:fang1}. The reflectors computed with the optimized slowness are flatter than that with the background slowness, demonstrating the higher accuracy of the optimized slowness.

\plot{opti}{width=\columnwidth,height=0.8\textheight}{(a) Optimized slowness after 4 iterations with 41 function and gradient evaluations. (b) Zero-subsurface offset section of the one-way shot-profile migration with the  optimized slowness model.[{\bf CR}]}

\plot{fang1}{width=\columnwidth,height=0.8\textheight}{ADCIGs obtained with true slowness (top), background slowness (center), and optimized slowness  (bottom).[{\bf CR}]}

\section{conclusion}
I showed the detailed derivation of the gradient of the objective function for the image-space wave-equation tomography by means of the adjoint-state method. I used the augmented Lagrangian methodology, which simplifies the understanding of the different operators image-space wave-equation tomography is composed of. The derivation is valid whether one considers migration of shot-profiles or generalized sources. Image-space phase-encoded gathers significantly accelerate image-space wave-equation tomography. In addition to decreasing the computational cost, the numerical example shows that image-space phase-encoded gathers can be used to generate sufficiently accurate optimized slowness.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The image obtained with the correct slowness model, ${\bf r} = {\bf r}({\bf x},{\bf h})$, can be expanded in Taylor series with respect to slowness around the background slowness, $\widehat{\bf s}=\widehat{\bf s}({\bf x})$, as
%\begin{eqnarray}
%{\bf r} = {\widehat{\bf r}} + \left. \frac{\partial {\bf r}}{\partial {\bf s}} \right |_{{\bf s}={\bf \widehat{s}}} \Delta {\bf s} + ..., 
%\label{eq:taylor1}
%\end{eqnarray}
%where ${\widehat{\bf r}}$ is the background image computed with the known background slowness and $\Delta {\bf s} = \Delta {\bf s} ({\bf x})$ is the slowness perturbation.
%\par
%By keeping up to the first order terms, we define the linear operator, $\bf T(\widehat{\bf s})$, that maps perturbations in slowness to perturbations in the image
%\begin{eqnarray}
%\Delta {\bf r} = {\bf T} \Delta {\bf s}.  
%\label{eq:tomo_op}
%\end{eqnarray}
%Starting from the cross-correlation imaging condition in the shot-profile migration
%\begin{eqnarray}
%{\bf r}({\bf x},{\bf h}) = \sum_{{\bf x}_s}\sum_{\omega} {\bf d}^*({\bf x}-{\bf h},{\bf x}_s,\omega) {\bf u}({\bf x}+{\bf h},{\bf x}_s,\omega), \label{eq:imgcond}
%\end{eqnarray}
%where ${\bf 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 and ${\bf h}=(h_x,h_y,h_z)$ is the subsurface half-offset, and `$^*$' stands for the complex-conjugate. Considering only one shot, by applying the chain rule to equation \ref{eq:imgcond}, we obtain, after dropping the dependencies on $({\bf x},{\bf h})$,
%\begin{eqnarray}
%\Delta {\bf r} = \sum_{\omega} \Delta {\bf d}^*(\omega) {\bf u}(\omega) + {\bf d}^*(\omega) \Delta {\bf u}(\omega),
%\label{eq:pert_imgcond}
%\end{eqnarray}
%where $\Delta {\bf d}(\omega)$ is the perturbed source wavefield and $\Delta u(\omega)$ is the perturbed receiver wavefield. 
%
%We can compute the adjoint of equation \ref{eq:imgcond} as:
%\begin{subequations} 
%\label{subeqn:list} 
%\begin{align} 
%%\begin{eqnarray}
%p({\bf x},{\bf x}_s,\omega) = \sum_{\bf h} u({\bf x}+{\bf h},{\bf x}_s,\omega) r({\bf x},{\bf h}) \label{subeqn:adjimgcond1} \\
%u({\bf x},{\bf x}_s,\omega) = \sum_{\bf h} p({\bf x}-{\bf h},{\bf x}_s,\omega) r({\bf x},{\bf h}) . \label{subeqn:adjimgcond2}
%%\end{eqnarray}
%\end{align} 
%\end{subequations} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\bibliographystyle{seg}
\bibliography{sep139a}
%\end{document}

%\append{}
%The perturbed wavefields satisfy the one-way wave equation linearized with respect to the slowness.
%%\begin{eqnarray}
%%\left\{ \begin{array}{l}
%%\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {\widehat s}^2 -|{\bf k}|^2}\right) \Delta d(\omega)= 
%%\frac{- i\omega \Delta s}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}^2}}}{\widehat d}(\omega) \\
%%\Delta d(z=0,\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 k}|^2}\right) \Delta u(\omega)=
%%\frac{- i\omega \Delta s}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}^2}}}{\widehat u}(\omega) \\
%%\Delta u(z=0,\omega) = 0 \end{array} \right., \label{eq_wveq_dU}
%%\end{eqnarray}
%%where ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector.
%%\par 
%Take, for instance, the recursive solution of the one-way wave equation for the source wavefield
%\begin{eqnarray}
%{\bf d}^{z+\Delta z} = e^{ik_z\Delta z}{\bf d}^{z} \label{eq:1way_sol0}
%\end{eqnarray}
%to be linearized with respect to the slowness evaluated around the background slowness. The vertical wavenumber, $k_z$, can be defined as the sum of the background wavenumber, $\widehat{k}_z$, and the wavenumber perturbation, $\Delta k_z$. By substituting this definition for $k_z$ into equation \ref{eq:1way_sol0} we get
%\begin{eqnarray}
%{\bf d}^{z+\Delta z} &=& e^{i{\widehat k}_z\Delta z}e^{i\Delta k_z\Delta z}{\bf d}^z \nonumber \\
%      & = & e^{i\Delta k_z\Delta z} {\widehat{\bf d}}^{z+\Delta z}, \label{eq:1way_sol1}
%\end{eqnarray}
%where
%\begin{eqnarray}
%\Delta k_z = \left. \frac{dk_z}{d{\bf s}} \right|_{{\bf s}={\widehat{\bf s}}} \Delta {\bf s}, \label{eq:1way_sol2}
%\end{eqnarray}
%and
%\begin{eqnarray}
%\left. \frac{dk_z}{d{\bf s}} \right|_{{\bf s}={\widehat{\bf s}}} = \frac{\omega^2 {\widehat s}}{\sqrt{\omega^2 {\widehat s}^2-\left |{\bf k} \right |^2}}, \nonumber
%\end{eqnarray}
%where ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector.
%\par
%The exponential in equation \ref{eq:1way_sol1} can be linearized using $e^{i\theta}\approx1+\theta$. With this linearization and equation \ref{eq:1way_sol2}, the linearized version of equation \ref{eq:1way_sol1} reads 
%\begin{eqnarray}
%{\bf d}^{z+\Delta z} &=& (1+i\frac{dk_z}{d{\bf s}}\Delta {\bf s} \Delta z){\widehat{\bf d}}^{z+\Delta z}. \label{eq:1way_sol3}
%\end{eqnarray}
%By rearranging the terms we finally get
%\begin{eqnarray}
%\Delta {\bf d}^{z+\Delta z} &=& i\frac{dk_z}{d{\bf s}}\Delta {\bf s} \Delta z {\widehat{\bf d}}^{z+\Delta z} = e^{i{\widehat k}_z\Delta z} \{i\frac{dk_z}{d{\bf s}}\Delta {\bf s} \Delta z {\widehat{\bf d}}^{z} \}. \label{eq:1way_sol4}
%\end{eqnarray}
%The perturbed source wavefield, $\Delta {\bf d}^{z+\Delta z}$, in equation \ref{eq:1way_sol4} is computed by applying the scattering operator, $i\frac{dk_z}{d{\bf s}}\Delta {\bf s} \Delta z$, to the background source wavefield at the previous depth level generating the scattered source wavefield. Then, the scattered source wavefield is propagated to the next depth level using the background slowness. The same reasoning can be applied to compute the perturbed receiver wavefield, which gives
%\begin{eqnarray}
%\Delta {\bf u}^{z+\Delta z} &=& e^{-i{\widehat k}_z\Delta z} \{-i\frac{dk_z}{d{\bf s}}\Delta {\bf s} \Delta z {\widehat{\bf u}}^{z} \}. \label{eq:1way_sol5}
%\end{eqnarray}
%\par
%


%\append{}
%In matrix form, for every depth level, the terms in the equations \ref{eq:pertsou1}, \ref{eq:pertrec1}, \ref{eq:pertimg1}, \ref%{eq:pscatt} and \ref{eq:uscatt} can be written as
%\begin{equation}
%\Delta{\bf p}(\omega)=\left [ \begin{array}{c}
%0\\
%\Delta p_2(\omega)\\
%\vdots\\
%\Delta p_{n_z-1}(\omega)\\
%\Delta p_{n_z}(\omega)
%\end{array}
%\right ]; 
%\Delta{\bf u}(\omega)=\left [ \begin{array}{c}
%0\\
%\Delta u_2(\omega)\\
%\vdots\\
%\Delta u_{n_z-1}(\omega)\\
%\Delta u_{n_z}(\omega)
%\end{array}
%\right ]; 
%\Delta{\bf r}=\left [ \begin{array}{c}
%0\\
%\Delta r_2\\
%\vdots\\
%\Delta r_{n_z-1}\\
%\Delta r_{n_z}
%\end{array}
%\right ]; 
%\nonumber
%\end{equation}
%\begin{equation}
%{\bf \widetilde p}(\omega)=\left [ \begin{array}{c}
%0\\
%\widetilde p_2(\omega)\\
%\vdots\\
%\widetilde p_{n_z-1}(\omega)\\
%\widetilde p_{n_z}(\omega)
%\end{array}
%\right ]; 
%{\bf \widetilde u}(\omega)=\left [ \begin{array}{c}
%0\\
%\widetilde u_2(\omega)\\
%\vdots\\
%\widetilde u_{n_z-1}(\omega)\\
%\widetilde u_{n_z}(\omega)
%\end{array}
%\right ]; 
%\nonumber
%\end{equation}
%\begin{equation}
%{\bf \widehat P}(\omega)=\left [ \begin{array}{ccccccc}
%\widehat P_1(\omega)    & 0             & \hdots    & 0         & 0        \\
%0      & \widehat P_2(\omega)           & \hdots    & 0         & 0        \\
%\vdots & \vdots        & \ddots    & \vdots    & \vdots   \\
%0      & 0             & \hdots    & \widehat P_{n_z-1}(\omega) & 0        \\
%0      & 0             & \hdots    & 0         & \widehat P_{n_z}(\omega)                           
%\end{array}
%\right ]; \nonumber
%\end{equation}
%
%\begin{equation}
%\Delta {\bf P}(\omega)=\left [ \begin{array}{ccccccc}
%0      & \Delta P_2(\omega)           & \hdots    & 0         & 0        \\
%\vdots & \vdots        & \ddots    & \vdots    & \vdots                  \\
%0      & 0             & \hdots    & \Delta P_{n_z-1}(\omega) & 0        \\
%0      & 0             & \hdots    & 0         & \Delta P_{n_z}(\omega)  \\
%0    & 0             & \hdots    & 0         & 0        
%\end{array}
%\right ]; \nonumber 
%\end{equation}
%
%\begin{equation}
%\Delta {\bf U}(\omega)=\left [ \begin{array}{ccccccc}
%0      & \Delta U_2(\omega)           & \hdots    & 0         & 0        \\
%\vdots & \vdots        & \ddots    & \vdots    & \vdots                  \\
%0      & 0             & \hdots    & \Delta U_{n_z-1}(\omega) & 0        \\
%0      & 0             & \hdots    & 0         & \Delta U_{n_z}(\omega)  \\
%0    & 0             & \hdots    & 0         & 0        
%\end{array}
%\right ]; 
%\nonumber
%\end{equation}
%
%\begin{equation}
%{\bf T}^{\downarrow}(\omega)=\left [ \begin{array}{cccccc}
%0      & T_2^{\downarrow}(\omega)  & \hdots  & 0                              & 0                             \\
%\vdots & \vdots                    & \ddots  & \vdots                         & \vdots                        \\
%0      & 0                         & \hdots  & T_{n_z-1}^{\downarrow}(\omega) & 0                             \\
%0      & 0                         & \hdots  & 0                              &T_{n_z}^{\downarrow}(\omega) \\
%0      & 0                         & \hdots  & 0                              & 0                           
%\end{array}
%\right ];  \nonumber
%\end{equation}
%
%\begin{equation}
%\Delta {\bf T}^{\downarrow}(\omega)=\left [ \begin{array}{cccccc}
%0      & \Delta T_2^{\downarrow}(\omega)  & \hdots  & 0                                     & 0                                   \\
%\vdots & \vdots                           & \ddots  & \vdots                                & \vdots                              \\
%0      & 0                                & \hdots  & \Delta T_{n_z-1}^{\downarrow}(\omega) & 0                                   \\
%0      & 0                                & \hdots  & 0                                     &\Delta T_{n_z}^{\downarrow}(\omega)\\
%0      & 0                                & \hdots  & 0                                     & 0                           
%\end{array}
%\right ]; \nonumber
%\end{equation}
%\par
%In a furthermore compact notation, now for every frequency, we can write
%\begin{equation}
%\Delta{\bf \underline p}=\left [ \begin{array}{c}
%\Delta \bf p(\omega_1)\\
%\Delta \bf p(\omega_2)\\
%\vdots\\
%\Delta \bf p(\omega_{n_{\omega}-1})\\
%\Delta \bf p(\omega_{n_{\omega}}))
%\end{array}
%\right ]; 
%\Delta{\bf \underline u}=\left [ \begin{array}{c}
%\Delta \bf u(\omega_1)\\
%\Delta \bf u(\omega_2)\\
%\vdots\\
%\Delta \bf u(\omega_{n_{\omega}-1})\\
%\Delta \bf u(\omega_{n_{\omega}}))
%\end{array}
%\right ]; 
%\nonumber
%\end{equation}
%\begin{equation}
%{\widetilde {\underline {\bf p}}}=\left [ \begin{array}{c}
%{\widetilde {\bf p}}(\omega_1)\\
%{\widetilde {\bf p}}(\omega_2)\\
%\vdots\\
%{\widetilde {\bf p}}(\omega_{n_{\omega}-1})\\
%{\widetilde {\bf p}}(\omega_{n_{\omega}})
%\end{array}
%\right ]; 
%{\widetilde {\underline {\bf u}}}=\left [ \begin{array}{c}
%{\widetilde {\bf u}}(\omega_1)\\
%{\widetilde {\bf u}}(\omega_2)\\
%\vdots\\
%{\widetilde {\bf u}}(\omega_{n_{\omega}-1})\\
%{\widetilde {\bf u}}(\omega_{n_{\omega}})
%\end{array}
%\right ]; 
%\nonumber
%\end{equation}
%\begin{equation}
%{\widehat {{\bf P}}}=\left [ \begin{array}{ccccccc}
%{\widehat {\bf P}}(\omega_1)    & 0             & \hdots    & 0         & 0        \\
%0      & {\widehat {\bf P}}(\omega_2)           & \hdots    & 0         & 0        \\
%\vdots & \vdots        & \ddots    & \vdots    & \vdots   \\
%0      & 0             & \hdots    & {\widehat {\bf P}}(\omega_{n_{\omega}-1}) & 0        \\
%0      & 0             & \hdots    & 0         & {\widehat {\bf P}}(\omega_{n_{\omega}})                           
%\end{array}
%\right ]; \nonumber
%\end{equation}
%
%\begin{equation}
%\Delta {\bf P}=\left [ \begin{array}{ccccccc}
%\Delta {\bf P}(\omega_1)&0           & \hdots    & 0         & 0        \\
%0      & \Delta {\bf P}(\omega_2)           & \hdots    & 0         & 0        \\
%\vdots & \vdots        & \ddots    & \vdots    & \vdots                  \\
%0      & 0             & \hdots    & \Delta {\bf P}(\omega_{n_{\omega}-1}) & 0        \\
%0      & 0             & \hdots    & 0         & \Delta {\bf P}(\omega_{n_{\omega}})  
%\end{array}
%\right ]; \nonumber 
%\end{equation}
%
%\begin{equation}
%\Delta {\bf U}=\left [ \begin{array}{ccccccc}
%\Delta {\bf U}(\omega_1)&0           & \hdots    & 0         & 0        \\
%0      & \Delta {\bf U}(\omega_2)           & \hdots    & 0         & 0        \\
%\vdots & \vdots        & \ddots    & \vdots    & \vdots                  \\
%0      & 0             & \hdots    & \Delta {\bf U}(\omega_{n_{\omega}-1}) & 0        \\
%0      & 0             & \hdots    & 0         & \Delta {\bf U}(\omega_{n_{\omega}})  
%\end{array}
%\right ]; 
%\nonumber
%\end{equation}
%
%\begin{equation}
%{\bf T}^{\downarrow}=\left [ \begin{array}{cccccc}
%{\bf T}^{\downarrow}(\omega_1)  & 0      & \hdots  & 0                              & 0                             \\
%0      & {\bf T}^{\downarrow}(\omega_2)  & \hdots  & 0                              & 0                             \\
%\vdots & \vdots                    & \ddots  & \vdots                         & \vdots                        \\
%0      & 0                         & \hdots  & {\bf T}^{\downarrow}(\omega_{n_{\omega}-1}) & 0                      \\
%0      & 0                         & \hdots  & 0                              &{\bf T}^{\downarrow}(\omega_{n_{\omega}}) 
%\end{array}
%\right ];  \nonumber
%\end{equation}
%
%\begin{equation}
%\Delta {\bf T}^{\downarrow}=\left [ \begin{array}{cccccc}
%\Delta {\bf T}^{\downarrow}(\omega_1)  & 0      & \hdots  & 0                              & 0                             \\
%0      & \Delta {\bf T}^{\downarrow}(\omega_2)  & \hdots  & 0                              & 0                             \\
%\vdots & \vdots                    & \ddots  & \vdots                         & \vdots                        \\
%0      & 0                         & \hdots  & \Delta {\bf T}^{\downarrow}(\omega_{n_{\omega}-1}) & 0                      \\
%0      & 0                         & \hdots  & 0                              &\Delta {\bf T}^{\downarrow}(\omega_{n_{\omega}}) 
%\end{array}
%\right ]; \nonumber
%\end{equation}
%
%The contribution of every frequency to the image perturbation is achieved by the summation operator, $\bf S$,
%\begin{equation}
%{\bf S}=\left [ \begin{array}{ccccc}
%{\bf I}(n_{\bf x} n_z) & {\bf I}(n_{\bf x} n_z) & \hdots & {\bf I}(n_{\bf x} n_z) & {\bf I}(n_{\bf x} n_z)
%\label{eq:sprd1}
%\end{array}
%\right ],
%\end{equation}
%composed of $n_{\omega}$ block-identity matrices of $n_{\bf x} n_z$-dimension. 
%\par


\end{document}
