\footer{July--2009}
\section{Image-space wave-equation tomography in the shot-profile domain according to the adjoint-state method}

\subsection{Image-space wave-equation tomography}
Image-space wave-equation tomography aims to solve for the slowness model, $s=s({\bf x})$, which minimizes the linearized objective function
\begin{eqnarray}
J({\bf s})=\frac{1}{2} \left | \left | \Delta r \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. $\Delta r$ is computed either by applying wave-equation migration-velocity analysis (WEMVA) \cite[]{sava2004a,sava2004b} or differential-semblance optimization (DSO) \cite[]{shen:VE49} operators. Here,  these operators are called indistinctively $\bf M$. Operator $\bf M$ is independent of the slowness and operates on the image obtained with the current slowness model, $\widehat r=\widehat r({\bf x},{\bf h})$. Therefore, the objective function reads 
\begin{eqnarray}
J(s)=\frac{1}{2} \left | \left | {{\bf M} \widehat r} \right | \right | ^2.
\label{eq:objfun1}
\end{eqnarray}
The gradient of this objective function with respect to the slowness evaluated at the current slowness, ${\widehat s}={\widehat s}(\bf x)$, is
\begin{eqnarray}
\nabla J(s)= \left. \left ( \frac{\partial \widehat r}{\partial s} \right )'\right |_{s=\widehat s} {\bf M'}  {{\bf M} \widehat r}. 
\label{eq:gradobjfun}
\end{eqnarray}
where $'$ denotes the adjoint. The linear operator $\frac{\partial \widehat r}{\partial s}$ defines a linear map between the slowness perturbation $\Delta s$  to the image perturbation $\Delta r$, and it is called image-space wave-equation tomographic operator.
\par
The image-space wave-equation tomographic operator is composed of different operators. However, this is not clear from the representation of equation \ref{eq:gradobjfun}. Therefore, for a clear explanation of the operators involved, I use the adjoint-state method to derive the gradient of the objective function in equation \ref{eq:objfun1}.
\par
Considering the shot profile migration, 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_{{\bf x}_s}\sum_{\omega} p_z^*({\bf x}-{\bf h},{\bf x}_s,\omega) u_z({\bf x}+{\bf h},{\bf x}_s,\omega), 
\label{eq:imgcond1}
\end{eqnarray}
where $p_z({\bf x},{\bf x}_s,\omega)$ is the source wavefield for a single frequency $\omega$ at horizontal coordinates ${\bf x}=(x,y)$ with the source located at ${\bf x}_s=(x_s,y_s,0)$; $u_z({\bf x},{\bf x}_s,\omega)$ is the receiver wavefield and ${\bf h}=(h_x,h_y)$ is the subsurface half-offset, and `$^*$' stands for the complex-conjugate. 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{P}}_z^*(\omega) u_z(\omega), 
\label{eq:imgconcompact}
\end{eqnarray}
where ${\bf{P}}$ is a convolutional matrix composed of $(h_x,h_y)$-shifted versions of $p_z({\bf x},{\bf x}_s,\omega)$.
\par
For subsequent depth levels, $p({\bf x},{\bf x}_s,\omega)$ is computed by means of the recursive downward propagation 
\begin{eqnarray}
\left\{ \begin{array}{l}
p_{z+1}(\omega)  = T_z^{\downarrow}(\omega,s)p_z \\ 
p_1(\omega) = f_s(\omega)\delta({\bf x - x_s}),
\label{eq:downmod}
\end{array}
\right.
\end{eqnarray}
where $T_z^{\downarrow}$ is the downward continuation operator, which is function of the slowness, $s$, and $f_s(\omega)$ is the source signature located at $x_s$. In equation \ref{eq:downmod}, I omitted the dependencies on ${\bf x}$ and ${\bf x}_s$.
\par
By taking the adjoint of equation \ref{eq:imgconcompact}, the receiver wavefield can be computed by convolving the source wavefield with the reflectivity
\begin{eqnarray}
u_z(\omega) = {\bf{P}}_z^(\omega) r_z, 
\label{eq:rec01}
\end{eqnarray}
where I omitted again the dependencies on $\bf x$ and, for ${\bf{P}}_z^(\omega)$ and  $r_z$, on $\bf h$.
\par
The upward propagation of the receiver wavefield and recording at the surface are performed by 
\begin{eqnarray}
\left\{ \begin{array}{l}
u_{z-1}(\omega)  = T_z^{\uparrow}(\omega,s)u_z + {\bf{P}}_z^(\omega) r_z \\ 
d(\omega) = u_1(\omega),
\label{eq:upmod01}
\end{array}
\right.
\end{eqnarray}
where $d(\omega)$ is the data collected at the surface.
\par
The use of equation \ref{eq:upmod01} for deriving the adjoint state equations and computing the adjoint state variables is cumbersome. When using one-way wave-equation, \cite{biondi:1340} shows the equivalence between shot-profile and source-receiver migration provided that the source function, $f_s(\omega)$, is an impulse at zero time and a delta function in space and crosscorrelation imaging condition is used in shot-profile migration. Without proving (yet), I consider that this equivalence is also valid for modeling and rewrite equation \ref{eq:upmod01} as
\begin{eqnarray}
\left\{ \begin{array}{l}
u_{z-1}(\omega)  = T_z^{\uparrow}(\omega,s)u_z + r_z \\ 
d(\omega) = u_1(\omega),
\label{eq:upmod02}
\end{array}
\right.
\end{eqnarray}
which is the modeling equation of the receiver wavefield in the source-receiver strategy.

\begin{eqnarray}
u^{z}(\omega) = {\bf{P}}^z(\omega) r^z, 
\label{eq:imgconcompact}
\end{eqnarray}
where ${\bf{P}}^z$ is a convolutional matrix. 

\par
The recursive upward propagation of $u^z(\omega)$ and recording at the surface are given by the following equations
\begin{eqnarray}
\left\{ \begin{array}{l}
u^{z-1}(\omega)  = P^z(\omega) r^z + T^z_{_\uparrow}(\omega,s)u^z \\ 
d(\omega) = u^1(\omega)
\label{eq:upwmod}
\end{array}
\right.
\end{eqnarray}
where $T^z_{_\uparrow}$ is the upward continuation operator, which is function of the slowness, $s$.
\par
In a matrix form, the terms in the above equations can be writen as
\begin{equation}
{\bf u}(\omega)=\left [ \begin{array}{c}
u^1(\omega)\\
u^2(\omega)\\
\vdots\\
u^{n_z-1}(\omega)\\
u^{n_z}(\omega)
\label{eq:uvec1}
\end{array}
\right ]; {\bf T}_{\uparrow}(\omega)=\left [ \begin{array}{cccccc}
0      & T^2_{\uparrow}(\omega)  & \hdots  & 0                            & 0                           \\
\vdots & \vdots                  & \ddots  & \vdots                       & \vdots                      \\
0      & 0                       & \hdots  & T^{n_z-1}_{\uparrow}(\omega) & 0                           \\
0      & 0                       & \hdots  & 0                            &T^{n_z-1}_{\uparrow}(\omega) \\
0      & 0                       & \hdots  & 0                            & 0                           
\end{array}
\right ]; \nonumber
\end{equation}
\par
\begin{equation}
{\bf r}=\left [ \begin{array}{c}
r^1\\
r^2\\
\vdots\\
r^{n_z-1}\\
r^{n_z}
\label{eq:rvec1}
\end{array}
\right ]; {\bf P}(\omega)=\left [ \begin{array}{ccccccc}
P^1(\omega)    & 0             & \hdots    & 0         & 0        \\
0      & P^2(\omega)           & \hdots    & 0         & 0        \\
\vdots & \vdots        & \ddots    & \vdots    & \vdots   \\
0      & 0             & \hdots    & P^{n_z-1}(\omega) & 0        \\
0      & 0             & \hdots    & 0         & P^{n_z}(\omega)                           
\end{array}
\right ]. \nonumber
\end{equation}
\par
In a furthermore compact notation, we can write
\begin{equation}
{\bf {\underline u}}=\left [ \begin{array}{c}
{\bf u}(\omega_1)\\
{\bf u}(\omega_2)\\
\vdots\\
{\bf u}(\omega_{n_{\omega}-1})\\
{\bf u}(\omega_{n_{\omega}})\\
\label{eq:uvec2}
\end{array}
\right ]; {\bf T}_{\uparrow}=\left [ \begin{array}{cccccc}
{\bf T}_{\uparrow}(\omega_1)     & 0                                & \hdots                       & 0                                         & 0                                         \\
 0                               & {\bf T}_{\uparrow}(\omega_2)     & \hdots                       & 0                                         & 0                                         \\
\vdots                           & \vdots                           & \ddots                       & \vdots                                    & \vdots                                    \\
0                                & 0                                & \hdots                       & {\bf T}_{\uparrow}(\omega_{n_{\omega}-1}) & 0                                         \\
0                                & 0                                & \hdots                       & 0                                         & {\bf T}_{\uparrow}(\omega_{n_{\omega}}) 
\end{array}
\right ]; \nonumber
\end{equation}
\par
\begin{equation}
{\bf P}=\left [ \begin{array}{ccccccc}
{\bf P}(\omega_1)    & 0             & \hdots    & 0         & 0        \\
0      & {\bf P}(\omega_2)           & \hdots    & 0         & 0        \\
\vdots & \vdots        & \ddots    & \vdots    & \vdots   \\
0      & 0             & \hdots    & {\bf P}(\omega_{n_{\omega}-1}) & 0        \\
0      & 0             & \hdots    & 0         & {\bf P}(\omega_{n_{\omega}})                           
\end{array}
\right ]; \nonumber
\end{equation}
\par
The reflectivity is replicated for every frequency by the spreading matrix, $\bf S$,
\begin{equation}
{\bf S}=\left [ \begin{array}{c}
{\bf I}(n_x n_z)\\
{\bf I}(n_x n_z)\\
\vdots\\
{\bf I}(n_x n_z)\\
{\bf I}(n_x n_z)
\label{eq:sprd1}
\end{array}
\right ],
\end{equation}
composed of $n_{\omega}$ block-identity matrices of $n_x n_z$-dimension. Finally the restriction matrix , $\bf E$,
\begin{equation}
{\bf E}=\left [ \begin{array}{ccccccc}
{\bf e}    & 0             & \hdots    & 0     & 0        \\
0      & {\bf e}           & \hdots    & 0     & 0        \\
\vdots & \vdots        & \ddots    & \vdots    & \vdots   \\
0      & 0             & \hdots    & {\bf e}   & 0        \\
0      & 0             & \hdots    & 0         & {\bf e}                           
\end{array}
\right ]; \nonumber
\end{equation}
where 
\begin{equation}
{\bf e}=\left [ \begin{array}{ccccc}
1    & 0      & \hdots    & 0     & 0 \end{array} \right ], \nonumber
\end{equation}
selects the wavefield ${\bf {\underline u}}$ at $z=0$. The resultant matrix representation for equation \ref{eq:upwmod} is
\begin{eqnarray}
\begin{array}{l}
{\underline{\bf u}} = {\bf P} {\bf Sr} + {\bf T}_\uparrow {\underline{\bf u}}\\ 
{\bf d} = {\bf E}{\underline{\bf u}}
\label{eq:upwmodmatrix}
\end{array}
\end{eqnarray}



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



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{sep136a}



