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

\usepackage{amsmath}

\setfigdir{Fig}

\begin{document}
%\documentclass{segabs}
%\documentclass[paper,twocolumn,twoside]{geophysics}
%\documentclass[manuscript,revised]{geophysics}
%\newcommand{\rs}[1]{\mathstrut\mbox{\scriptsize\rm #1}}
%\newcommand{\rr}[1]{\mbox{\rm #1}}
%\usepackage{amsmath}
%\usepackage{multicol}
%\usepackage{placeins}
%\renewcommand{\figdir}{./Fig} % figure directory
%\begin{document}
%\renewcommand{\thefootnote}{\fnsymbol{footnote}} 
 
\footer{July--2009}
\section{Adjoint-state method for least-squares migration in the shot-profile domain}

\subsection{Upward-continuation modeling in the shot-profile domain}
Starting from the cross-correlation imaging condition in the shot-profile migration
\begin{eqnarray}
r({\bf x},{\bf h}) = \sum_{{\bf x}_s}\sum_{\omega} p^*({\bf x}-{\bf h},{\bf x}_s,\omega) u({\bf x}+{\bf h},{\bf x}_s,\omega), \label{eq:imgcond}
\end{eqnarray}
where $p({\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. 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} 
In a more compact notation, not explicitly writing the dependencies on $x$,  $y$ and $\bf h$, and for a given depth level, $z$, equation \ref{subeqn:adjimgcond2} can be written as:
\begin{eqnarray}
u^{z}(\omega) = {\bf{P}}^z(\omega) r^z, 
\label{eq:imgconcompact}
\end{eqnarray}
where ${\bf{P}}^z$ is a convolutional matrix. If $r=r({\bf x},h_x)$, for the case of 5 $x$ positions and $3$ subsurface offsets and subsurface-offset interval equals to the CMP interval, equation \ref{eq:imgconcompact} can be explicitly written as:
\begin{equation} 
\left[\begin{array}{l}
u_2\\ 
u_3\\ 
u_4\\ 
u_5\\ 
u_6\\ 
\end{array}\right]=
\left[\begin{array}{lllllllllllllll}
0 & p_1& 0  & p_3& 0  & 0  & 0  & 0  & 0  & 0  & 0  & 0  & 0  & 0  & 0 \\ 
0 & 0  & 0  & 0  & p_2& 0  & p_4& 0  & 0  & 0  & 0  & 0  & 0  & 0  & 0 \\ 
0 & 0  & 0  & 0  & 0  & p_1& 0  & p_3& 0  & p_5& 0  & 0  & 0  & 0  & 0 \\ 
0 & 0  & 0  & 0  & 0  & 0  & 0  & 0  & p_2& 0  & p_4& 0  & 0  & 0  & 0 \\ 
0 & 0  & 0  & 0  & 0  & 0  & 0  & 0  & 0  & 0  & 0  & p_3& 0  & p_5& 0 \\ 
\end{array}\right]
\left[\begin{array}{l}
r_{1,1}\\
r_{1,2}\\
r_{1,3}\\
r_{2,1}\\
r_{2,2}\\
r_{2,3}\\
r_{3,1}\\
r_{3,2}\\
r_{3,3}\\
r_{4,1}\\
r_{4,2}\\
r_{4,3}\\
r_{5,1}\\
r_{5,2}\\
r_{5,3}\\
\end{array}\right]
\end{equation}
where the subscripts of $u$ and $p$ represent $x$ positions. The first and second subscripts of $r$ represent $x$ positions and subsurface offset, respectively.
\par
For subsequent depth levels, $p(\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)
\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$,
which explicit dependency has not been showed.
\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}

\subsection{Weighted least-squares data fitting}

The optimum reflectivity, $\bf r$, minimizes the objective function
\begin{eqnarray}
J({\bf r})=h({\bf \underline{u}(r),r})=\frac{1}{2} \left | \left | {\bf W} \left ( {\bf d - d_{obs}}) \right ) \right | \right | ^2.
\label{eq:objfun}
\end{eqnarray}

The forward equation reads:
\begin{eqnarray}
F({\bf \underline{u}(r),r})=\left ( {\bf I} - {\bf T}_\uparrow \right ){\bf {\underline u}} - {\bf PSr} = 0.
\label{eq:forwmod}
\end{eqnarray}

\begin{itemize}
\item Build the augmented functional
\begin{eqnarray}
\begin{array}{l}
\mathcal{L}(\underline{\lambda},{\bf \underline{u}},{\bf r}) = h({\bf \underline{u},r}) - \langle \underline{\lambda}, F({\bf \underline{u},r}) \rangle \\
\\
\mathcal{L}(\underline{\lambda},{\bf \underline{u}},{\bf r}) = \frac{1}{2} \left | \left | {\bf W} \left ( {\bf d - d_{obs}}) \right ) \right | \right | ^2 - \langle \underline{\lambda},{\left ( {\bf I} - {\bf T}_\uparrow \right ) \bf {\underline u}} - {\bf PSr} \rangle.
\end{array}
\end{eqnarray}
\item Define the adjoint state equations
\begin{eqnarray}
\frac{\partial \mathcal{L}(\underline{\lambda},{\bf \underline{u}},{\bf r})}{\partial {\bf \underline{u}}} = {\bf E}^T {\bf W}^2({\bf d - d_{obs}}) - \left ( {\bf I} - {\bf T}_\uparrow \right )^* \lambda = 0.
\end{eqnarray}
The adjoint state equations define the scalar recursive downward-continuation equations 
\begin{equation}
\lambda^{z+1}(\omega)  = T^z_{_\downarrow}(\omega)\lambda^z(\omega) 
\label{eq:downlambda}
\end{equation}
with the initial conditions
\begin{equation}
\lambda^{1}(\omega)  = w(\omega) \left ( d(\omega) - d_{obs}(\omega)\right ).
\label{eq:downlambda}
\end{equation}
In matrix notation, the adjoint state equations are:
\begin{eqnarray}
{\underline \lambda} = {\bf T}_{\downarrow} {\underline \lambda}+{\bf E}^T {\bf W}^2({\bf d - d_{obs}}).
\end{eqnarray}
\item Compute the gradient with respect to $\bf r$
\begin{eqnarray}
\frac{\partial \mathcal{L}(\underline{\lambda},{\bf \underline{u}},{\bf r})}{\partial {\bf r}} = \frac{\partial \langle \underline{\lambda},{\bf PSr} \rangle}{\partial {\bf r}} = \frac{\partial \langle \left ({\bf PS}\right )^*{\underline\lambda},{\bf r} \rangle}{\partial {\bf r}} = {\bf S^* P^* }{\underline\lambda}.
\end{eqnarray}
\end{itemize}
that corresponds to the downward propagation of the residuals and cross-correlation imaging condition of the downward propagated  residuals with the source wavefield.




\end{document}
