\def\figdir{Fig}
\lefthead{Guerra}
\righthead{Gradient of ISWET by the adjoint-state method}
%\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 in the generalized source domain by the adjoint-state method }
\author{Claudio Guerra}

\ABS{
Optimization using gradient descent techniques requires the computation of the gradient of the objective function. The adjoint-state method is an efficient way of computing the gradient, since is not necessary to compute the expensive  Frech$\acute {\rm e}$t derivatives. Here,  I derive the gradient of the image-space wave-equation tomography using the adjoint state method and show its application with a numerical example.
}

\section{introduction}
Wave-equation tomography aims to solve for Earth models that explain the data under some norm. In general, two categories are distinguished depending on the domain 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, generally called here as image-space wave-equation tomography (ISWET) from which wave-equation migration velocity analysis (WEMVA) \cite{sava2004a, sava2004b} and differential semblance velocity analysis (DSVA) \cite{shen2004, shen:VE49} are two variants, the objective function to be minimized is defined in the image space. 
\par
Despite their differences, those two categories are linked by the concept of extended modeling \cite{symes2008}. \longcite{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.
\par
Similar to waveform inversion, ISWET is a computationally demanding process. To decrease its computational cost it is common practice to use generalized sources \cite{shen:VE49,yax2008}. Because wavefield propagation is a linear process, generalized sources can be computed by linearly combining source wavefields and receiver wavefields, using phase encoding techniques \cite{whitmore1995, romero2000}. In particular, \longcite{guerra2009a} used generalized source functions synthesized by phase encoding modeling experiments in the image space drastically decreasing the cost of DSVA. The example showed in this report also uses this kind of generalized source functions.
\par
The aforementioned ISWET methods differ on the way the objective function is computed and, consequently, the numerical optimization scheme applied. As pointed out by \longcite{biondi2008}, in WEMVA, an intermediate residual-moveout parameter picking is necessary to produce the perturbed image via linearized residual prestack depth migration \cite{sava2004a}. The picking procedure makes it difficult to automate the whole process. As the computed perturbed image is consistent with the wave-equation tomographic operator, however, conjugate gradient methods can be used to invert for the slowness perturbation. In DSVA, the perturbed image is computed by applying the fully automated differential-semblance operator (DSO) on the background image, both using the subsurface-offset gathers (ODCIG) or angle gathers (ADCIG). When applied to ODCIG, DSVA minimizes the energy not focused at zero-offset. When applied to ADCIG, DSVA minimizes energy departing from flatness of the reflectors. Notwithstanding the ability of automating ISWET, DSVA produces perturbed images not consistent with the wave-equation tomographic operator. To minimize the objective function computed with DSO, quasi-Newton algorithms can be used for which the gradient of the objective function needs to be computed.
\par
The gradient of the objective function can be obtained by computing the Frech$\acute {\rm e}$t derivatives. However, this computation even for 2D applications of ISWET can be very expensive. An efficient way of computing the gradient without the need of Frech$\acute {\rm e}$t derivatives is the adjoint-state method \cite{chavent:476,plessix2006}. \longcite{plessix2006} describes two methodologies for computing the gradient of the objective function using the adjoint-state method. The more intuitive is the one that uses the augmented Lagrangian. The augmented Lagrangian is formed by the objective function and the scalar product of the adjoint-state variables and general solutions of the forward modeling equations. The adjoint-state variables are solution of a system of adjoint-state equations defined by equating to zero the derivative 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 
Here, I show how to compute the gradient of the ISWET objective function using the augmented Lagrangian methodology. The derivation of the gradient presented here is valid whether ISWET uses areal-shot migration or shot-profile migration. \longcite{yax2008} provide a complete derivation of the forward and adjoint wave-equation tomographic operators in the generalized source domain. The goal of this paper is to provide a more detailed derivation of the gradient of the ISWET objective function using the adjoint state method than the ones available, for example, \longcite{shen:2132}. 
\par
The references on the image-space phase-encoded gathers are \longcite{biondi2006, biondi2007} and \longcite{guerra2008b}. For completeness, here I briefly describe how to compute these phase-encoded gathers. Then, I derive the ISWET gradient using the adjoint-state method and, finally show an example of DSVA using image-space phase-encoded gathers on the Marmousi model. 

\section {Image-space generalized sources}
The prestack exploding-reflector modeling, introduced by \longcite{biondi2006} and \longcite{biondi2007}, synthesizes areal data and the corresponding areal source function, having as initial condition, a wave-equation migrated prestack image obtained with an inaccurate slowness model. In general, one single reflection event from one single ODCIG is modeled by the 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 the isolated SODCIGs 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 SODCIGs according to the apparent geological dip of the reflector. This rotation maintains the velocity information needed for migration velocity analysis, especially for dipping reflectors \cite{biondi2007}. $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 SODCIG 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$ is such that it 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 be applied in a target-oriented manner.
\par
However, if one considers to model one single reflection from one single SODCIG at a time, the modeled dataset can be orders of magnitude bigger than the original data set. As discussed by \longcite{biondi2006,guerra2008a,guerra2008b}, by modeling several reflectors and several SODCIGs simultaneously using random phase encoding a much smaller dataset is obtained. 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 SODCIGs 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 chosen to be the 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. Recursively solving Equations \ref{eq_perm_bD2} and \ref{eq_perm_bU2} gives us 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)$, which can be collected at the depth level, $z_c$. 

\section{Gradient computation of ISWET by the adjoint state method}
Image-space wave-equation tomography aims to solve for the slowness model, $s=s({\bf x})$, that minimizes the linearized 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 which measures the goodness of the slowness model. Notice that the objective function depends on the slowness model through the image perturbation. $\Delta r$ is computed by applying a differential residual-focusing operator $\bf M$ \cite{biondi2008}, using either differential residual prestack-migration \cite{sava2004a,sava2004b} or differential-semblance optimization (DSO) \cite{shen:VE49} operators to the image $r = r({\bf x},{\bf h})$. Herein, operators are represented by bold capital letters.
\par
If the differential residual-focusing operator $\bf M$ is independent on 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 `$'$' 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 \widehat r}{\partial s}$ defines the mapping,  $\frac{\partial \widehat r}{\partial s} \Delta s = \Delta r$, between the slowness perturbation $\Delta s$ and the image perturbation and it is called image-space wave-equation tomographic operator.
\par
As 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 (areal shot and shot-profile migration, e.g.), 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 and ${\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 to 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:
\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 function of the slowness, $s$, and $q(\omega)$ is the source wavefield used as 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 $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}
\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, 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 representation of equations \ref{eq:pertsou}, \ref{eq:pertrec}, \ref{eq:pertimg} is
\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}
\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 a generalized source scheme. They depend on the state variables $\Delta {\underline {\bf d}}$, $\Delta {\underline {\bf u}}$ and $\Delta {\underline {\bf r}}$. Following the the augmented functional methodology for computing the adjoint-states described in \longcite{plessix2006}, 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 equal 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 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 crosscorrelated in time with the scattered wavefields.

\section{Example}
The gradient of the ISWET objective function as computed in the previous section can be used in a quasi-Newton optimization scheme. I use L-BFGS-B bound constrained optimization algorithm \cite{nocedal} to invert for slowness using the Marmousi model. The gradient was B-spline smoothed according to  smoothing of the gradient After computing the gradient375 two-way shots at every 24 m were modeled using the original Marmousi model of Figure \ref{fig:true}a. The maximum offset is 6600 m. Figure \ref{fig:true}b shows the one-way shot profile image with the correct slowness model. These data are also input to one-way shot profile migration using the background slowness of Figure \ref{fig:bckg}a to compute the background image of Figure \ref{fig:bckg}b. When comparing the zero-offset section of the two images, it is apparent the pull-up effect in the background image of migrating with a too slow velocity. 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.
\activeplot{true}{width=\columnwidth,height=0.8\textheight}{CR}{(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}]}
\activeplot{bckg}{width=\columnwidth,height=0.8\textheight}{CR}{(a) Background slowness. (b) Zero-subsurface offset section of the one-way shot-profile migration with the  background slowness model.[{\bf CR}]}
\activeplot{dslow}{width=\columnwidth,height=0.4\textheight}{ER}{ Slowness ratio.[{\bf ER}]}
\par
Following \longcite{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.
\activeplot{bwkim}{width=\columnwidth,height=0.4\textheight}{CR}{Zero-subsurface offset section showing the reflectors selected on the one-way shot-profile migration with the  background slowness model.[{\bf CR}]}
\par
Figure \ref{fig:opti}a shows the optimized slowness model after 4 iterations comprising 41 function and gradient evaluations. The optimized slowness decreased as expected, however without recovering the details of the true slowness. The smoothness of the optimized slowness can be partially credited to the B-spline filtering of the gradient and also because image-space wave-equation tomography is not able to solve for the short wavelengths of the slowness model, as pointed out by \longcite{symes2008}. Figure \ref{fig:opti}b shows the shot-profile migration of the original shots using the optimized slowness. Notice that the pull-up effect was greatly mitigated and the reflectors are more focused than in the background image.
\activeplot{opti}{width=\columnwidth,height=0.8\textheight}{CR}{(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}]}

Figure \ref{fig:fang1} shows, from top to bottom, the ADCIGs resulting from shot-profile migration of the original shots using the true slowness, the background slowness and the optimized slowness. The reflectors computed with the optimized slowness are flatter than that with the background slowness, attesting the higher accuracy of the optimized slowness.
\activeplot{fang1}{width=\columnwidth,height=0.8\textheight}{CR}{ADCIGs obtained with true slowness (top), background slowness (center), and optimized slowness  (bottom).[{\bf CR}]}

\section{conclusion}
I used the straightforward augmented Lagrangian methodology to derive the adjoint-state of the gradient of the objective function for the image-space wave-equation tomography problem. The derivation is valid whether one considers migration of shot-profiles or generalized sources. Image-space phase-encoded gathers were used to optimize the slowness model. This significantly accelerates wave-equation tomography. The quasi-Newton optimization was able to increase the accuracy of the slowness as the numerical example shows. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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}

%\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

