\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})$, that 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 to the image $r = r({\bf x},{\bf h})$. Here,  these operators are called indistinctively $\bf M$. Therefore, the objective function reads 
\begin{eqnarray}
J(s)=\frac{1}{2} \left | \left | {{\bf M} r} \right | \right | ^2.
\label{eq:objfun1}
\end{eqnarray}
If $\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 \widehat r}{\partial s} \right )'\right |_{s=\widehat s} {\bf M'}  {{\bf M} \widehat r}. 
\label{eq:gradobjfun}
\end{eqnarray}
where `$'$' denotes the adjoint and  $\widehat r=\widehat r({\bf x},{\bf h})$ is the image obtained with the current slowness model.  The linear operator $\frac{\partial \widehat r}{\partial s}$ defines a mapping 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
In 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$, considering just one shot, is computed by the crosscorrelation
\begin{eqnarray}
r_z({\bf x},{\bf h}) = \sum_{\omega} p_z^*({\bf x}-{\bf h},\omega) u_z({\bf x}+{\bf h},\omega), 
\label{eq:imgcond1}
\end{eqnarray}
where $p_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. 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{P}}_z'(\omega) u_z(\omega) = {\bf{S}}{\bf{U}}_z(\omega) p_z^*(\omega), 
\label{eq:imgconcompact}
\end{eqnarray}
where ${\bf{P}}$ and ${\bf{U}}$ are convolutional matrices composed of $(h_x,h_y)$-shifted versions of $p_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, $p({\bf x},\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(\omega) \\ 
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 ${\bf x}_s=(x_s,y_s,0)$. 
\par
The downward continuation of the receiver wavefield is performed by
\begin{eqnarray}
\left\{ \begin{array}{l}
u_{z+1}(\omega)  = T_z^{\downarrow}(\omega,s)u_z(\omega) \\ 
u_1(\omega) = d(\omega),
\label{eq:upmod01}
\end{array}
\right.
\end{eqnarray}
where $d(\omega)$ is the data at the surface. In equations \ref{eq:downmod} and \ref{eq:upmod01}, I omitted the dependencies of the wavefield with respect to ${\bf x}$.
\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 $p + \Delta p$, $u + \Delta u$ and, consequently, $r + \Delta r$ are physical realizations with $s + \Delta s$. To the first order, these perturbed fields are given by 
\begin{eqnarray}
\Delta p_{z+1}(\omega) = T_z^{\downarrow}(\omega,s) \Delta p_z(\omega) + \Delta T_z^{\downarrow}(\omega,s) p_z(\omega),
\label{eq:pertsou}
\end{eqnarray}
\begin{eqnarray}
\Delta u_{z+1}(\omega) = T_z^{\downarrow}(\omega,s) \Delta u_z(\omega) + \Delta T_z^{\downarrow}(\omega,s) u_z(\omega),
\label{eq:pertrec}
\end{eqnarray}

In equations \ref{eq:pertsou} and \ref{eq:pertrec}, $\Delta T_z^{\downarrow}$ is the scattering operator
\begin{eqnarray}
\Delta T_z^{\downarrow}(\omega,s) =  i\frac{\omega^2 s}{\sqrt{\omega^2 s^2 - \left | {\bf k} \right | ^2}} dz \Delta s.
\label{eq:scattering}
\end{eqnarray}
Its derivation is provided in the Apendix A.
\par
The perturbed image is
\begin{eqnarray}
\Delta r_{z} = {\bf S} \left (\Delta {\bf P}_z'(\omega) u_z(\omega)+ {\bf P}_z'(\omega) \Delta u_z(\omega) \right )= {\bf S} \left ({\bf U}_z(\omega) \Delta p_z^*(\omega) + {\bf P}_z'(\omega)\Delta U_z(\omega)  \right ).
\label{eq:pertimg}
\end{eqnarray}
Notice that if the perturbations are evaluated around the current slowness, $\widehat s$, equations \ref{eq:pertsou}, \ref{eq:pertrec} and \ref{eq:pertimg} read
\begin{eqnarray}
\Delta p_{z+1}(\omega) = T_z^{\downarrow}(\omega,\widehat s) \Delta p_z(\omega) + G_z \widehat p_z(\omega) \Delta s_z,
\label{eq:pertsou1}
\end{eqnarray}
\begin{eqnarray}
\Delta u_{z+1}(\omega) = T_z^{\downarrow}(\omega,\widehat s) \Delta u_z(\omega) + G_z \widehat u_z(\omega) \Delta s_z,
\label{eq:pertrec1}
\end{eqnarray}
\begin{eqnarray}
\Delta r_{z} = {\bf S} \left ({\widehat {\bf U}}_z(\omega) \Delta p_z^*(\omega) + {\widehat {\bf P}}_z'(\omega)\Delta u_z(\omega)  \right ).
\label{eq:pertimg1}
\end{eqnarray}
which are the non-renormalized equations. The operator $G_z$, is the scattering operator given by
\begin{eqnarray}
G_z(\omega) =  \pm i\frac{\omega^2 \widehat s}{\sqrt{\omega^2 \widehat s^2 - \left | {\bf k} \right | ^2}} dz,
\label{eq:pscatt}
\end{eqnarray}
where $\pm$ stands for the source and receiver scattering, respectively.
\par
However, if we write $p_z = {\widehat p_z}+\Delta p_z$ and $u_z = {\widehat u_z}+\Delta u_z$, equations \ref{eq:pertsou} and \ref{eq:pertrec} read:
\begin{eqnarray}
\Delta p_{z+1}(\omega) = T_z^{\downarrow}(\omega,\widehat s) \Delta p_z(\omega) + G_z(\omega,\widehat s) (\widehat p_z(\omega) + \Delta p_z(\omega)) \Delta s_z,
\label{eq:pertsou1.1}
\end{eqnarray}
and
\begin{eqnarray}
\Delta u_{z+1}(\omega) = T_z^{\downarrow}(\omega,s) \Delta u_z(\omega) + G_z(\omega,\widehat s) (\widehat u_z(\omega) + \Delta u_z(\omega)) \Delta s_z.
\label{eq:pertrec1.1}
\end{eqnarray}
The equation for the perturbed image stands the same.
\par
 The matrix representation for equations \ref{eq:pertsou1.1}, \ref{eq:pertrec1.1}, \ref{eq:pertimg1} is
\begin{eqnarray}
\Delta {\underline {\bf p}} = {\bf T}^{\downarrow} \Delta {\underline {\bf p}} + {\bf G}({\bf \widehat P}+\Delta {\bf \underline P}) {\bf S}' \Delta {\bf s},
\label{eq:pertsou2}
\end{eqnarray}
\begin{eqnarray}
\Delta {\underline {\bf u}} = {\bf T}^{\downarrow} \Delta {\underline {\bf u}} +  {\bf G}({\bf \widehat U}+\Delta {\bf \underline 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 p}}^* + {\widehat {\bf P}}' \Delta {\underline {\bf u}}  \right ). 
\label{eq:pertimg2}
\end{eqnarray}
where the diagonal operators ${\widehat {{\bf P}}}$, ${\widehat {{\bf U}}}$, $\Delta {\bf \underline P}$ and $\Delta {\bf \underline U}$ are composed of the background source wavefield, background receiver wavefield, perturbed source wavefield and perturbed receiver wavefield, respectively.
\par
Equations \ref{eq:pertsou2}, \ref{eq:pertrec2} and \ref{eq:pertimg2} are the forward equations of the image-space wave-equation tomography problem using a shot profile scheme and the renormalization strategy. They depend on the state variables $\Delta {\underline {\bf p}}$, $\Delta {\underline {\bf u}}$ and $\Delta {\underline {\bf r}}$. The augmented functional reads
\begin{eqnarray}
\begin{array}{l}
\mathcal{L}(\Delta {\underline {\bf p}}, \Delta {\underline {\bf u}}, \Delta {\underline {\bf r}}, {\underline \lambda}_p, {\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}_p,\left ( {\bf I} - {\bf T}^\downarrow \right - {\bf G})    \Delta {\underline {\bf p}} - {\bf G} {\widehat {{\bf P}}} {\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 {{\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 p}}^* + {\widehat {\bf P}}' \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}_p = {\widehat {\bf U}} {\underline \lambda}_r, \label{subeqn:adjstvar1} \\
\left ( {\bf I} - {\bf T}^\downarrow \right )' {\underline \lambda}_u = {\widehat {\bf P}} {\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}_p + {\widehat {\bf U}} {\underline \lambda}_r, \label{subeqn:adjstvar10} \\
{\underline \lambda}_u = {\bf T}^\uparrow {\underline \lambda}_u + {\widehat {\bf P}} {\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 {{\bf P}}}'{\underline \lambda}_p  +  {\widetilde {{\bf U}}}'{\underline \lambda}_u \right ).
\label{eq:adjstgrad}
\end{eqnarray}
To compute the gradient, the upward propagated perturbed wavefields, ${\underline \lambda}_p$ and ${\underline \lambda}_u$, are crosscorrelated in time with the scattered wavefields.

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

\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 matrix, $\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


