\title{Wave-equation tomography in the generalized source domain}
\author{Yaxun Tang, Claudio Cardoso and Biondo Biondi}
\footer{SEP-136}

\maketitle

\begin{abstract}
In this paper, we generalize wave-equation migration velocity analysis to composite shot domain. 
Any linear composite shot gather can be used for WEMVA. One particular case is the plane-wave shot gathers. 
I will demonstate some theoretical advantages of doing WEMVA in composite shot domain, 
\end{abstract}

\section{introduction}

\section{wave-equation tomography}
Wave-equation tomography is a non-linear inverse problem, 
it is based on an optimization that tries to find an optimal background slowness that minimizes the residual field, ${\bf \Delta I}$, defined in the image space.
The residual field is derived from the background image, ${\bf I}$, which is computed with a background slowness (or the current estimate of the slowness). 
The residual field measures the correctness of the background slowness,  
its minimum (under some norm, e.g. $\ell_2$) is achieved when a correct background slowness has been used for migration. 
There are many choices of the residual field, such as residual moveout in the Angle-Domain Common-Image Gathers (ADCIGs), 
differential semblence in the ADCIGs, 
reflection angle stacking power (in this case we have to maximize the residual field, or minimize the minus stacking power), and etc..
For simplicity, we define a linear operator ${\bf K}$ that measures the image perturbation as follows:
\begin{eqnarray}
{\bf \Delta I} = ({\bf 1} - {\bf K}){\bf I},
\end{eqnarray}
where ${\bf 1}$ is the identity operator.
In the wave-equation migration velocity analysis (WEMVA) method, ${\bf K}$ is the differential residual migration operator (Sava and Biondi???), 
it flattens the ADCIGs and also focuses the migrated image in the CMP direction; In the Differential Semblance Optimization (DSO) method,
\begin{eqnarray}
{\bf K} = {\bf 1} - {\bf R},
\end{eqnarray}
where ${\bf R}$ is the DSO operator either in the subsurface offset domain or in the angle domain. 
Subsurface offset domain DSO focuses the energy at zero offset while angle domain DSO flattens the ADCIGs. 

In general, if we choose $\ell_2$ norm, the objective function to minimize can be written as follows:
\begin{eqnarray}
J = \frac{1}{2}||{\bf \Delta I}||_2 = \frac{1}{2}||({\bf 1}-{\bf K}){\bf I}||_2,
\end{eqnarray}
where $||\cdot||_2$ stands for the $\ell_2$ norm.
Gradient based optimization techniques such as Quasi-Newton method and Conjugate Gradient method can be used to minimize the objective function $J$, 
which is non-linear because the image ${\bf I}$ is a non-linear function of the slowness ${\bf s}$. 
To use the gradient based methods, we have to provide the
gradient of $J$ with respect to the slowness ${\bf s}$.
The gradient reads
\begin{eqnarray}
{\bf \nabla J} = \Re \left( \left( \frac{\partial {\bf I}}{\partial {\bf s}}\right)^{'}({\bf 1}-{\bf K})^{'}({\bf 1}-{\bf K}){\bf I}\right),
\end{eqnarray}
where $^{'}$ means the adjoint, $\frac{\partial{\bf I}}{\partial{\bf s}}$ is a linear operator which defines the mapping from
the slowness perturbation ${\bf \Delta s}$ to the image perturbation ${\bf \Delta I}$. This is easy to see by
expanding the image ${\bf I}$ around a background slowness ${\widehat {\bf s}}$ as follows:
\begin{eqnarray}
{\bf I} = {\widehat {\bf I}} + \frac{\partial{\bf I}}{{\partial{\bf s}}}({\bf s}-{\widehat {\bf s}}) + \cdots,
\end{eqnarray} 
where ${\widehat{\bf I}}$ is the background image computed with the background slowness ${\widehat {\bf s}}$. 
Keeping only the zero and first order terms, we get the linearized operator $\frac{\partial{\bf I}}{\partial{\bf s}}$ defined as follows:
\begin{eqnarray}
{\bf \Delta I} = \frac{\partial{\bf I}}{\partial{\bf s}} {\bf \Delta s} = {\bf T}{\bf \Delta s},
\end{eqnarray}
where ${\bf \Delta I} = {\bf I} - {\widehat {\bf I}}$ and ${\bf \Delta s} = {\bf s} - {\widehat {\bf s}}$. 
${\bf T} = \frac{\partial{\bf I}}{\partial{\bf s}}$ is the wave-equation tomographic operator.
The tomographic operator can either be defined in the source and receiver domain (Sava,????) or in the shot-profile domain (Shen,???).
In next section we follow a similar approach disscussed by Shen,?? and review the forward and adjoint tomographic operator 
in the shot profile domain, in the subsequent sections, we generalize the expression of the tomographic operator to 
generalized source domains.

\section{the tomographic operator in the shot-profile domain}

For the conventional shot-profile migration, both source and receiver wavefields are downward continued with the 
following one-way wave equations (Claerbout,???):
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right)D({\bf x},{\bf x}_s,\omega) = 0 \\
D(x,y,z=0,{\bf x}_s,\omega) = f_s(\omega)\delta({\bf x}-{\bf x}_s) \end{array} \right., \label{eq_wveq_bD}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right)U({\bf x},{\bf x}_s,\omega) = 0 \\
U(x,y,z=0,{\bf x}_s,\omega) = Q(x,y,z=0,{\bf x}_s,\omega) \end{array} \right., \label{eq_wveq_bU}
\end{eqnarray}
where $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 for a single frequency $\omega$ at image point ${\bf x}$ for the source located at ${\bf x}_s$.
$s({\bf x})$ is the slowness at ${\bf x}$, ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector, $f_s(\omega)$
is the freqnency dependent source signature and $f_s(\omega)\delta({\bf x}-{\bf x}_s)$ defines the point source function at ${\bf x}_s$, 
which serves as the boundary condition of equation \ref{eq_wveq_bD}. $Q(x,y,z=0,{\bf x}_s,\omega)$ is the shot gather for the shot located at ${\bf x}_s$, 
which serves as the boundary condition of equation \ref{eq_wveq_bU}.
To produce the image, the following cross-correlation imaging condition is used:
\begin{eqnarray}
I({\bf x},{\bf h}) = \sum_{{\bf x}_s}\sum_{\omega} \overline{D({\bf x}-{\bf h},{\bf x}_s,\omega)} U({\bf x}+{\bf h},{\bf x}_s,\omega), \label{eq_ic_bI}
\end{eqnarray}
where ${\bf h}=(h_x,h_y,h_z)$ is the subsurface half offset.

The perturbed image at some image point can be derived by a simple application of the chain rule
to equation \ref{eq_ic_bI}: 
\begin{eqnarray}
\Delta I({\bf x},{\bf h})  &=& 
\sum_{{\bf x}_s}\sum_{\omega} 
\left( \overline{\Delta D({\bf x}-{\bf h},{\bf x}_s,\omega)} {\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega) + \right. \nonumber \\
& & \left. \overline{{\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)} \Delta U ({\bf x}+{\bf h},{\bf x}_s,\omega) \right), \label{eq_ic_dI}
\end{eqnarray}
where ${\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)$ and ${\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega)$ are the background source and 
receiver wavefields computed with the background slowness ${\widehat s}({\bf x})$; 
$\Delta D({\bf x}-{\bf h},{\bf x}_s,\omega)$ and $\Delta U({\bf x}+{\bf h},{\bf x}_s,\omega)$ 
are the perturbed source wavefield and perturbed receiver wavefield, which are the results of the slowness perturbation $\Delta s({\bf x})$.
The perturbed source and receiver wavefields satisfy the following 
linearized (with respect to slowness) one-way wave equations, respectively(see Appendix A for derivations)
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-\sqrt{\omega ^2 {\widehat s}^2 ({\bf x})-|{\bf k}|^2}\right) \Delta D({\bf x},{\bf x}_s,\omega)= 
\frac{+ i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}({\bf x})}}}{\widehat D}({\bf x},{\bf x}_s,\omega) \\
\Delta D(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right. , \label{eq_wveq_dD}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+\sqrt{\omega ^2 {\widehat s}^2 ({\bf x})-|{\bf k}|^2}\right) \Delta U({\bf x},{\bf x}_s,\omega)=
\frac{- i\omega \Delta s({\bf x})}{\sqrt{1-\frac{|{\bf k}|^2}{\omega^2 {\widehat s}({\bf x})}}}{\widehat U}({\bf x},{\bf x}_s,\omega) \\
\Delta U(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right.. \label{eq_wveq_dU}
\end{eqnarray}

Equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} are linear with respect to slowness perturbation, 
thus they define linear mappings from the slowness perturbation $\Delta s({\bf x})$ to
the perturbed source wavefield $\Delta D({\bf x},{\bf x}_s,\omega)$ and the perturbed receiver wavefield $\Delta U({\bf x},{\bf x}_s,\omega)$, respectively.
Therefore, recursively solving equations \ref{eq_wveq_dD} and \ref{eq_wveq_dU} give us the perturbed source and receiver wavefields.
The perturbed source and receiver wavefields are then used in equation \ref{eq_ic_dI} to generate the perturbed image $I({\bf x},{\bf h})$, where
the background source and receiver wavefields are precomputed by solving the recursive equations \ref{eq_wveq_bD} and \ref{eq_wveq_bU} with a
background slowness ${\widehat s}(\bf x)$.
Appendix B gives a more detailed matrix representation of how to evaluate the forward tomographic operator ${\bf T}$.

To compute the adjoint tomographic operator ${\bf T}'$, we first apply the adjoint of the imaging condition in equation \ref{eq_ic_dI} 
to get the perturbed source and receiver wavefields $\Delta D({\bf x},{\bf x}_s,\omega)$ and $\Delta U({\bf x},{\bf x}_s,\omega)$ as follows
\begin{eqnarray}
\Delta D({\bf x},{\bf x}_s,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h}){\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega) \\
\Delta U({\bf x},{\bf x}_s,\omega) &=& \sum_{\bf h} \Delta I({\bf x},{\bf h}){\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)
\end{eqnarray}
Then we solve the adjoint equations of \ref{eq_wveq_dD} and \ref{eq_wveq_dU} to get the slowness perturbation $\Delta s({\bf x})$. 
Again, in order to solve the adjoint equations of \ref{eq_wveq_dD} and \ref{eq_wveq_dU}, 
the background source wavefield ${\widehat D}({\bf x},{\bf x}_s,{\omega})$ and 
the background receiver wavefield ${\widehat U}({\bf x},{\bf x}_s,\omega)$ have to be computed in advance.
Appendix B gives a more detailed matrix representation of how to evaluate the adjoint tomographic operator: ${\bf T}^{'}$.

%--------------------------------------------------------------------
\section{tomography with the encoded wavefields}
%--------------------------------------------------------------------
However, the cost for computing the gradient of the objective function $J$ in the shot-profile domain 
is twice the cost of a shot-profile migration, to minimize $J$ many gradient and function evaluations have to be computed,
thus wave-equation tomography in the conventional shot-profile domain seems infeasible with nowadays computer resources.
To save the cost and to make this powerful method more practical, we extend the wave-eqaution tomography theory to
generalized source domains, or phase encoded source domain. We discuss two different strategies generate 
the generalized shot gathers, i.e., 

\subsection{Data-space encoded wavefields}
We encode the source wavefield and the corresponding receiver wavefield
\begin{eqnarray}
{\widetilde D}({\bf x},{\bf p}_s,\omega) &=& \sum_{{\bf x}_s} D({\bf x},{\bf x}_s,\omega) \alpha({\bf x}_s,{\bf p}_s,\omega) \\
{\widetilde U}({\bf x},{\bf p}_s,\omega) &=& \sum_{{\bf x}_s} U({\bf x},{\bf x}_s,\omega) \alpha({\bf x}_s,{\bf p}_s,\omega)
\end{eqnarray}
where $\alpha({\bf x}_s,{\bf p}_s,\omega)$ is some phase encoding function. In the case of plane-wave phase encoding, 
$\alpha({\bf x}_s,{\bf p}_s,\omega)$ is defined as
\begin{eqnarray}
\alpha({\bf x}_s,{\bf p}_s,\omega) = e^{i\omega{\bf p}_s{\bf x}_s},
\end{eqnarray}
where ${\bf p}_s$ is the ray parameter for the source plane waves on the surface.
Because of the linearity of the wave equation, the en

Imaging in the generallized source domain
\begin{eqnarray}
I({\bf x},{\bf h}) = \sum_{{\bf p}_s}\sum_{\omega} {\widetilde D}({\bf x}-{\bf h},{\bf p}_s,\omega) {\widetilde U}({\bf x}+{\bf h},{\bf p}_s,\omega)
\end{eqnarray}

\subsection{Image-space encoded wavefields}

\section{Numerical examples}

\append{}
This appendix derives the linearized one-way wave equation

\append{}
This appendix demonstrates a matrix representation of the forward tomographic operator ${\bf T}$
\begin{eqnarray}
{\bf D}_{z+\Delta z} = {\bf E}_z({\bf s}_z){\bf \Delta D}_z, \label{eq_app_down}
\end{eqnarray}
where
\begin{eqnarray}
{\bf E}_{z}({\bf s}_z) = e^{ik_z({\bf s}_z)\Delta z} = e^{i\sqrt{\omega^2{\bf s}_z^2-|{\bf k}|^2}} \label{eq_app_exp}
\end{eqnarray}

The perturbed source wavefield at some depth level can be derived from the background wavefield by a simple
application of the chain rule to equation \ref{eq_app_down}:
\begin{eqnarray}
{\bf \Delta D}_{z+\Delta z} = {\bf E}_z({\widehat{\bf s}}_z){\bf \Delta D}_z + {\bf \Delta E}_z({\widehat{\bf s}}_z){\bf D}_z, \label{eq_app_mat_dD}
\end{eqnarray}
where ${\bf \Delta E}_z$ represents the perturbed extrapolator, it can be obtained by a formal linearization 
with respect to slowness of the extrapolator defined in equation \ref{eq_app_exp}
\begin{eqnarray}
{\bf E}_z({\bf s}_z) = e^{ik_z({\bf s}_z)\Delta z}  
&\approx& e^{i\Delta z {\widehat k_z}} + 
e^{i\Delta z {\widehat k_z}}\left( i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z} \right) {\bf \Delta s}_z \nonumber \\
&=&{\bf E}_z(\widehat{{\bf s}}_z) + {\bf E}_z(\widehat{{\bf s}}_z) 
\left( i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z}\right)  {\bf \Delta s}_z
\end{eqnarray}
where ${\widehat k}_z=k_z({\widehat {\bf s}}_z)$. So the perturbed extrapolator reads
\begin{eqnarray}
{\bf \Delta E}_z({\widehat {\bf s}}_z) = {\bf E}_z(\widehat{{\bf s}}_z) 
\left( i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z}\right)  {\bf \Delta s}_z \label{eq_app_dE}
\end{eqnarray}
Substitute equation \ref{eq_app_dE} into \ref{eq_app_mat_dD}, we get
\begin{eqnarray}
{\bf \Delta D}_{z+\Delta z} = {\bf E}_z({\widehat{\bf s}}_z){\bf \Delta D}_z + 
{\bf E}_z(\widehat{{\bf s}}_z) 
\left( i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z}\right) {\bf D}_z {\bf \Delta s}_z, \label{eq_app_mat_dD2}
\end{eqnarray}
Let us define a scattering operator ${\bf G}_z$ that interacts with the background wavefield as follows:
\begin{eqnarray}
{\bf G}_z({\bf D}_z,{\widehat {\bf s}}_z) = \left( i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z} \right) {\bf D}_z 
= \frac{i\omega \Delta z}{\sqrt{1-\frac{|{\bf k}|^2}{\omega ^2 {\widehat {\bf s}}_z^2}}} {\bf D}_z
\end{eqnarray}
Then the perturbed source wavefield for depth level $z+\Delta z$ can be rewritten as follows:
\begin{eqnarray}
{\bf \Delta D}_{z+\Delta z} = {\bf E}_z({\widehat {\bf s}}_z){\bf \Delta D}_z + 
{\bf E}_z({\widehat {\bf s}}_z) {\bf G}_z({\widehat {\bf D}}_z,{\widehat {\bf s}}_z) {\bf \Delta s}_z. \label{eq_app_mat_dD3}
\end{eqnarray}
We can further write out the recursive equation \ref{eq_app_mat_dD3} for all depth levels in the following matrix form:
\tiny
\begin{eqnarray}
\left( \begin{array}{c}
{\bf \Delta D}_0 \\
{\bf \Delta D}_1 \\
{\bf \Delta D}_2 \\
{\vdots}\\
{\bf \Delta D}_n
\end{array} \right)
&=&
\left( \begin{array}{cccccc}
{\bf 0}   & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf 0}       & {\bf 0}  \\
{\bf E}_0 & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf 0}       & {\bf 0}  \\
{\bf 0}   & {\bf E}_1 & {\bf 0}  & {\cdots} & {\bf 0}       & {\bf 0}  \\
{\vdots}  & {\vdots}  & {\vdots} & {\vdots} & {\vdots}      & {\vdots} \\
{\bf 0}   & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf E}_{n-1} & {\bf 0}
\end{array} \right)
\left( \begin{array}{c}
{\bf \Delta D}_0 \\
{\bf \Delta D}_1 \\
{\bf \Delta D}_2 \\
{\vdots}\\
{\bf \Delta D}_n
\end{array} \right)
+ \nonumber \\
& & \left( \begin{array}{cccccc}
{\bf 0}   & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf 0}        & {\bf 0}\\{\bf E}_0 & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf 0}        & {\bf 0}\\
{\bf 0}   & {\bf E}_1 & {\bf 0}  & {\cdots} & {\bf 0}        & {\bf 0}\\
{\vdots}  & {\vdots}  & {\vdots} & {\vdots} & {\vdots}       & {\bf 0}\\
{\bf 0}   & {\bf 0}   & {\bf 0}  & {\cdots} & {\bf E}_{n-1}  & {\bf 0}
\end{array} \right)
\left( \begin{array}{ccccc}
{\bf G}_0 & {\bf 0}   & {\bf 0}    & {\cdots} & {\bf 0}  \\
{\bf 0}   & {\bf G}_1 & {\bf 0}    & {\cdots} & {\bf 0}  \\
{\bf 0}   & {\bf 0}   & {\bf G}_2  & {\cdots} & {\bf 0}  \\
{\vdots}  & {\vdots}  & {\vdots}   & {\vdots} & {\vdots} \\
{\bf 0}   & {\bf 0}   & {\bf 0}    & {\cdots} & {\bf G}_n
\end{array} \right)
\left( \begin{array}{c}
{\bf \Delta s}_0 \\
{\bf \Delta s}_1 \\
{\bf \Delta s}_2 \\
{\vdots}\\
{\bf \Delta s}_n
\end{array} \right) \nonumber
\end{eqnarray}
\normalsize
Or in a more compact notation
\begin{eqnarray}
{\bf \Delta D} = {\bf E}({\widehat {\bf s}}){\bf \Delta D} + 
{\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf D}},{\widehat {\bf s}}){\bf \Delta s}. \label{eq_app_compact_dD}
\end{eqnarray}
The solution of equation \ref{eq_app_compact_dD} can be formally written as follows:
\begin{eqnarray}
{\bf \Delta D} = \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} 
{\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf D}},{\widehat {\bf s}}){\bf \Delta s}. \label{eq_app_solution_dD}
\end{eqnarray}

Similarly, the receiver wavefield satisfy the following recursive relation:
\begin{eqnarray}
{\bf U}_{z+\Delta z} = \overline{{\bf E}_z({\widehat {\bf s}}_z)}{\bf \Delta U}_z, \label{eq_app_down_U}
\end{eqnarray}
where the over line stands for the complex conjugate, thus
\begin{eqnarray}
\overline{{\bf E}_{z}({\bf s}_z)} = e^{-ik_z({\bf s}_z)\Delta z} = e^{-i\sqrt{\omega^2{\bf s}_z^2-|{\bf k}|^2}\Delta z}.
\end{eqnarray}
Apply the chain rule to euqation \ref{eq_app_down_U}, we obtain the recursive equation for the scattered receiver wavefield:
\begin{eqnarray}
{\bf \Delta U}_{z+\Delta z} = \overline{{\bf E}_z({\widehat {\bf s}}_z)}{\bf \Delta U}_z - 
\overline{{\bf E}_z({\widehat {\bf s}}_z)} {{\bf G}_z({\widehat {\bf U}}_z,{\widehat {\bf s}}_z)} {\bf \Delta s}_z,
\end{eqnarray}
where ${\bf G}_z({\widehat {\bf U}}_z,{\widehat {\bf s}}_z)$ is the scattering operator interacts with the background receiver wavefield as follows:
\begin{eqnarray}
{\bf G}_z({\widehat {\bf U}}_z,{\widehat {\bf s}}_z) = 
\left( i\Delta z \left. \frac{dk_z}{d{\bf s}_z} \right|_{{\bf s}_z={\widehat {\bf s}}_z} \right) {\widehat {\bf U}}_z 
= \frac{i\omega \Delta z}{\sqrt{1-\frac{|{\bf k}|^2}{\omega ^2 {\widehat {\bf s}}_z^2}}} {\widehat {\bf U}}_z. \label{eq_app_dU}
\end{eqnarray}
We can further write out the recursive equation \ref{eq_app_dU} for all depth levels in the following matrix form:
\tiny
\begin{eqnarray}
\left( \begin{array}{c}
{\bf \Delta U}_0 \\
{\bf \Delta U}_1 \\
{\bf \Delta U}_2 \\
{\vdots}\\
{\bf \Delta U}_n
\end{array} \right)
&=&
\left( \begin{array}{cccccc}
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                  & {\bf 0}  \\
\overline{{\bf E}_0} & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                  & {\bf 0}  \\
{\bf 0}              & \overline{{\bf E}_1} & {\bf 0}  & {\cdots} & {\bf 0}                  & {\bf 0}  \\
{\vdots}             & {\vdots}             & {\vdots} & {\vdots} & {\vdots}                 & {\vdots} \\
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & \overline{{\bf E}_{n-1}} & {\bf 0}
\end{array} \right)
\left( \begin{array}{c}
{\bf \Delta U}_0 \\
{\bf \Delta U}_1 \\
{\bf \Delta U}_2 \\
{\vdots}\\
{\bf \Delta U}_n
\end{array} \right)
- \nonumber \\
& & \left( \begin{array}{cccccc}
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                   & {\bf 0}\\
\overline{{\bf E}_0} & {\bf 0}              & {\bf 0}  & {\cdots} & {\bf 0}                   & {\bf 0}\\
{\bf 0}              & \overline{{\bf E}_1} & {\bf 0}  & {\cdots} & {\bf 0}                   & {\bf 0}\\
{\vdots}             & {\vdots}             & {\vdots} & {\vdots} & {\vdots}                  & {\bf 0}\\
{\bf 0}              & {\bf 0}              & {\bf 0}  & {\cdots} & \overline{{\bf E}_{n-1}}  & {\bf 0}
\end{array} \right)
\left( \begin{array}{ccccc}
{\bf G}_0 & {\bf 0}   & {\bf 0}    & {\cdots} & {\bf 0}  \\
{\bf 0}   & {\bf G}_1 & {\bf 0}    & {\cdots} & {\bf 0}  \\
{\bf 0}   & {\bf 0}   & {\bf G}_2  & {\cdots} & {\bf 0}  \\
{\vdots}  & {\vdots}  & {\vdots}   & {\vdots} & {\vdots} \\
{\bf 0}   & {\bf 0}   & {\bf 0}    & {\cdots} & {\bf G}_n
\end{array} \right)
\left( \begin{array}{c}
{\bf \Delta s}_0 \\
{\bf \Delta s}_1 \\
{\bf \Delta s}_2 \\
{\vdots}\\
{\bf \Delta s}_n
\end{array} \right) \nonumber
\end{eqnarray}
\normalsize
Or in a more compact notation
\begin{eqnarray}
{\bf \Delta U} = \overline{{\bf E}({\widehat {\bf s}})}{\bf \Delta U} - 
\overline{{\bf E}({\widehat {\bf s}})}{\bf G}({\widehat {\bf U}},{\widehat {\bf s}}){\bf \Delta s}. \label{eq_app_dU2}
\end{eqnarray}
The solution of equation \ref{eq_app_dU2} can be formally written as follows:
\begin{eqnarray}
{\bf \Delta U} = \left( {\bf 1} - \overline{{\bf E}({\widehat {\bf s}})} \right)^{-1} \overline{{\bf E}({\widehat {\bf s}})}
{\bf G}({\widehat {\bf U}},{\widehat {\bf s}}){\bf \Delta s}. \label{eq_app_solution_dU}
\end{eqnarray}

Therefore, the perturbed image can be obtained as follows: 
\tiny
\begin{eqnarray}
\left( \begin{array}{c}
{\bf \Delta I}_0 \\
{\bf \Delta I}_1 \\
{\bf \Delta I}_2 \\
{\vdots}\\
{\bf \Delta I}_n
\end{array} \right)
=
\left( \begin{array}{ccccc}
{\widehat {\bf U}}_0 & {\bf 0}              & {\bf 0}               & {\cdots} & {\bf 0}  \\
{\bf 0}              & {\widehat {\bf U}}_1 & {\bf 0}               & {\cdots} & {\bf 0}  \\
{\bf 0}              & {\bf 0}              & {\widehat {\bf U}}_2  & {\cdots} & {\bf 0}  \\
{\vdots}             & {\vdots}             & {\vdots}              & {\vdots} & {\vdots} \\
{\bf 0}              & {\bf 0}              & {\bf 0}               & {\cdots} & {\widehat {\bf U}}_n
\end{array} \right)
\left( \begin{array}{c}
\overline{{\bf \Delta D}_0}\\
\overline{{\bf \Delta D}_1}\\
\overline{{\bf \Delta D}_2}\\
{\vdots}\\
\overline{{\bf \Delta D}_n}
\end{array} \right)
+
\left( \begin{array}{ccccc}
{\overline{\widehat{\bf D}}_0} & {\bf 0}                        & {\bf 0}                         & {\cdots} & {\bf 0}  \\
{\bf 0}                        & \overline{{\widehat{\bf D}}_1} & {\bf 0}                         & {\cdots} & {\bf 0}  \\
{\bf 0}                        & {\bf 0}                        & \overline{{\widehat{\bf D}}_2}  & {\cdots} & {\bf 0}  \\
{\vdots}                       & {\vdots}                       & {\vdots}                        & {\vdots} & {\vdots} \\
{\bf 0}                        & {\bf 0}                        & {\bf 0}                         & {\cdots} & \overline{{\widehat{\bf D}}_n}
\end{array} \right)
\left( \begin{array}{c}
{{\bf \Delta U}_0}\\
{{\bf \Delta U}_1}\\
{{\bf \Delta U}_2}\\
{\vdots}\\
{{\bf \Delta U}_n}
\end{array} \right) \nonumber
\end{eqnarray}
\normalsize
Or in a more compact notation
\begin{eqnarray}
{\bf \Delta I} = {\rm diag}\left({\widehat {\bf U}}\right)\overline{{\bf \Delta D}} + 
\overline{{\rm diag}\left( {\widehat{\bf D}}\right)}{\bf \Delta U}. \label{eq_app_compact_dI}
\end{eqnarray}
After substituting equations \ref{eq_app_solution_dD} and \ref{eq_app_solution_dU} into equation \ref{eq_app_compact_dI} 
\begin{eqnarray}
{\bf \Delta I} &=& \left( {\rm diag}\left( {\widehat {\bf U}}\right)
                    \overline{ \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} 
                    {\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf D}},{\widehat {\bf s}})} \right. + \nonumber \\
           & &    \left.  \overline{{\rm diag}\left({\widehat {\bf D}} \right)}
                    \left( {\bf 1} - \overline{{\bf E}({\widehat {\bf s}})} \right)^{-1} 
                    \overline{{\bf E}({\widehat {\bf s}})}{\bf G}({\widehat {\bf U}},{\widehat {\bf s}})
             \right){\bf \Delta s}, \label{eq_app_compact_dI2}
\end{eqnarray}
we finally get the following matrix representation of the forward tomographic operator ${\bf T}$: 
\begin{eqnarray}
{\bf T} &=& {\rm diag}\left( {\widehat {\bf U}}\right)
                    \overline{ \left( {\bf 1} - {\bf E}({\widehat {\bf s}}) \right)^{-1} 
                    {\bf E}({\widehat {\bf s}}){\bf G}({\widehat {\bf D}},{\widehat {\bf s}})} + \nonumber \\
           & &     \overline{{\rm diag}\left({\widehat {\bf D}} \right)}
                    \left( {\bf 1} - \overline{{\bf E}({\widehat {\bf s}})} \right)^{-1} 
                    \overline{{\bf E}({\widehat {\bf s}})}{\bf G}({\widehat {\bf U}},{\widehat {\bf s}})
\end{eqnarray}

\append{}
This appendix demonstrates a matrix representation of the adjoint tomographic operator ${\bf T}^{'}$
It is easy to get the adjoint equation from equation \ref{eq_app_compact_dI2}, 
