\title{\new{Target-oriented wave-equation least-squares migration/inversion with phase-encoded Hessian}}

\lefthead{Tang}
\righthead{Inversion with phase-encoded Hessian}

\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\ms{GEO-2008-0443}

\address{
\footnotemark[1]Stanford Exploration Project,\\
Geophysics Department,\\
Stanford University,\\
Stanford, CA 94305 \\
E-mail: tang@sep.stanford.edu \\
Presented at the 78th Annual Meeting, Society of Exploration Geophysicists, November 2008.
Manuscript submitted on 12/24/2008.
}

\author{Yaxun Tang\footnotemark[1]}


\maketitle

\begin{abstract}
Prestack depth migration produces blurred images due to limited acquisition apertures, complexities in the velocity model 
and bandlimited characteristics of seismic waves.
This distortion can be partially corrected using the model-space least-squares imaging/inversion approach, 
where a target-oriented wave-equation Hessian operator is first explicitly computed,
and then inverse filtering is iteratively applied
to deblur or invert for the reflectivity. 
However, one difficulty is the cost for computing the explicit Hessian operator, 
which requires storing a large number of Green's functions, making it challenging for large-scale applications.
As a remedy to this difficulty, I present a new method for computing the Hessian operator for 
the wave-equation-based least-squares migration/inversion problem. 
The proposed method modifies the original explicit 
Hessian formula, enabling efficient computation of this operator. A particular 
advantage of this method is that it eliminates on-disk storage of 
Green's functions.
The modifications, however, also introduce undesired crosstalk artifacts. 
I examine two different phase-encoding schemes, namely, 
plane-wave-phase encoding and random-phase encoding, to suppress the crosstalk. 
I apply the randomly phase-encoded Hessian operator to the Sigsbee2A synthetic data set, 
where an improved subsalt image with \new{more balanced amplitudes} is obtained.
\end{abstract}

\section{introduction}
\par
Migration is an important tool for imaging subsurface structures using reflection seismic data. The classic imaging principle 
for shot-based migration states that reflectors are 
located where the forward-propagated source wavefield correlates with the backward-propagated receiver wavefield \cite[]{claerbout1971}. 
However, this imaging principle is only the adjoint of the forward Born-modeling operator \cite[]{lailly1983}, 
which provides reliable structural information of the subsurface, but blurs the image
because of the non-unitary nature of the Born modeling operator. 
To deblur the migrated image and correct the effects of limited acquisition geometry, complex overburden and bandlimited wavefields, 
the imaging problem can be \new{posed} as an inverse problem based on the minimization of a least-squares functional.
The inverse problem can be formulated either in the data space \cite[]{lailly1983,tarantola1984,nemeth1999,kuhl2003,marie2005} or 
in the model space \cite[]{beylkin1985,chavent1999,rickett2003,sjoeberg2003,guitton2004,plessix2004,valenciano2006,yu2006,symes2008}.
The data-space approach can be solved iteratively using the gradient-based method \cite[]{nemeth1999,kuhl2003,marie2005}
without explicit construction of the Hessian, the matrix of the second derivatives of the error functional 
with respect to the model parameters. The iterative solving, however, is relatively costly and converges slowly without proper preconditioning.
\par
On the other hand, the model-space approach requires explicitly constructing the Hessian and applying its pseudo-inverse to the migrated 
image. 
The full Hessian of the least-squares functional is too big and expensive to be computed in practical applications; 
hence some authors \cite[]{chavent1999,rickett2003,plessix2004,symes2008}
approximate it by a diagonal matrix. In the case of high-frequency asymptotics, and with an infinite aperture, 
the Hessian is diagonal in most cases \cite[]{beylkin1985}.
For a finite range of frequencies and limited acquisition geometry, however, the Hessian is no longer diagonal and not even diagonally dominant 
\cite[]{pratt1998,chavent1999,plessix2004,valenciano2006}. It has been shown by \cite{albertin2004} and \cite{valenciano2008} that, 
in areas of poor illumination, e.g., subsalt regions, the Hessian's main diagonal energy
is smeared along 
its off-diagonals. Therefore, a diagonal matrix has limited effect in deblurring the migrated image, 
especially in poorly illuminated areas. That is why several authors, e.g., \cite{albertin2004} and \cite{valenciano2008}, 
suggest computing a limited number of 
the Hessian off-diagonals to compensate for poor illumination
and improve the inversion/deblurring result.

\par
Since the exact Hessian off-diagonals are expensive to compute, some attempts have been made 
to reduce the cost by computing the non-diagonal Hessian in an approximate sense. 
\cite{yu2006} introduce a lateral invariant non-diagonal Hessian by assuming a 1-D layered medium;
\cite{guitton2004} uses a bank of non-stationary filters to approximate a non-diagonal inverse of the Hessian;
with local plane-wave assumptions, \cite{lecomte1998}, \cite{gelius2002} and \cite{lecomte2008} compute 
the Hessian in the local phase domain using ray-based approach.
They demonstrate that, since ray tracing conveniently gives local propagation angles of both source and receiver rays, 
the local scattering wavenumber for image points in the subsurface can be efficiently constructed. 
However, the high-frequency asymptotic approximation and the caustics inherent in ray theory may prevent the ray-based approach accurately handle 
complex geologies \cite[]{hoffmann2001}.
To better honor the bandlimited nature of the seismic waves, 
\cite{xie2006} use one-way wave-equation-based approach to compute the phase-domain Hessian through local plane-wave decomposition.
The above theories of computing the
local phase-domain Hessian operator assume the velocity model is locally homogeneous.
Although this is a reasonable assumption for most cases, when
there are sharp velocity contrasts as in the vicinities of salt boundaries, this assumption becomes less reliable. 

\par
Another way of computing the \new{wave-equation} non-stationary Hessian operator is through crosscorrelation of the source and receiver Green's functions 
in the space domain \cite[]{plessix2004,valenciano2006}. 
This approach does not introduce any assumptions to the velocity model.
However, computing even a limited number of the Hessian off-diagonals in the space domain, by directly implementing 
the explicit Hessian formula, is very cumbersome. A huge number
of Green's functions (easily several hundred terabytes for a typical 3-D survey with a reasonable frequency band) must be 
pre-computed and stored and then read from the disk to generate the Hessian.
Such operations not only require high-volume storage, but also high-speed I/O and network communication.
Though computer speed continues to improve rapidly, computing the Hessian in such a way still presents a challenge. 

\par
To make the space-domain Hessian
more affordable, I describe a method based on the phase-encoding technique. In this method, the original explicit Hessian formula is
slightly modified to enable efficient computation of this operator. The proposed method makes the Hessian computation similar to 
the shot-profile migration, but with slightly modified imaging and boundary conditions for the wavefields. 
The new method eliminates the need to store Green's functions, but it also introduces crosstalk artifacts. 
I examine two phase-encoding schemes,
plane-wave-phase encoding \cite[]{whitmore1995,zhang2005,liu2006} and random-phase encoding \cite[]{romero2000}, to attenuate the crosstalk. 

\par
This paper is organized as follows. 
First, I briefly review the theory of formulating the inverse problem in the model space. 
Next, I discuss how the explicit Hessian can be efficiently computed 
using phase encoding. Finally, I apply the 
phase-encoded Hessian to deblur the migrated image for the Sigsbee2A model.

\section{least-squares hessian}
\par
{\new Within limits of the Born approximation of the acoustic wave equation, the seismic data $d({\bf x}_r, {\bf x}_s, \omega)$
recorded by a receiver at ${\bf x}_r=(x_r,y_r)$ due to a source at ${\bf x}_s=(x_s,y_s)$
can be modeled using a linear equation as follows \cite[]{stolt1986}}:
\begin{equation}
d({\bf x}_r, {\bf x}_s, \omega) = \omega^2 \sum_{{\bf x}} f_s(\omega) G({\bf x}, {\bf x}_s, \omega) G({\bf x}, {\bf x}_r, \omega) m({\bf x}),  \label{eq_fwd}
\end{equation}
where $\omega$ is the angular frequency; $f_s(\omega)$ is the source signature; 
and $m({\bf x})$ denotes the reflectivity (a perturbed quantity from the background velocity) at image point ${\bf x}=(x,y,z)$ in the subsurface;
$G({\bf x},{\bf x}_s,\omega)$ and $G({\bf x},{\bf x}_r,\omega)$
are the Green's functions connecting the source and receiver respectively, to the image point ${\bf x}$.
Throughout this paper, we assume the Green's functions are computed by means of 
one-way wavefield extrapolation \cite[]{claerbout1985,stoffa1990,ristow1994,biondi2002}. 
\new{Although not discussed in this paper, Green's functions obtained using other methods, 
such as by solving two-way wave equation, can also be
used under this framework.}

\par
In equation \ref{eq_fwd}, we assume ${\bf x}_s$ and ${\bf x}_r$
are infinite in extent and independent of each other. For a real survey, however, we do not have infinitely long cable and
infinitely many sources; thus we must introduce the following  acquisition mask matrix to limit the size of the modeling:
\begin{eqnarray}
w({\bf x}_r,{\bf x}_s) = \left\{ \begin{array}{ll}
                         1 & \mbox{ if } {\bf x}_r \mbox{ is within the recording} \\
                           & \mbox{ range of a shot at } {\bf x}_s; \\
                         0 & \mbox{ otherwise }. \end{array} \right. \label{eq_mask}
\end{eqnarray}
In 2-D, for a marine acquisition geometry, $w({\bf x}_r,{\bf x}_s)$ is similar to a band-limited diagonal matrix;
for an Ocean Bottom Seismometer (OBS) or a land acquisition geometry, where all shots share the same receiver array, 
$w({\bf x}_r,{\bf x}_s)$ is a rectangular matrix.
In 3-D, the structure of the matrix is slightly more complicated.

\par
To find a model that best fits the observed data, we can minimize the following objective function in the least-squares sense:
\begin{eqnarray}
F &=& \frac{1}{2}\sum_{\omega}\sum_{{\bf x}_s}\sum_{{\bf x}_r} 
|w({\bf x}_r,{\bf x}_s) \nonumber \\
& &\times [d({\bf x}_r,{\bf x}_s,\omega)-d_{\rm obs}({\bf x}_r,{\bf x}_s,\omega)]|^2, \label{eq_objfnc}
\end{eqnarray}
where $d_{\rm obs}({\bf x}_r,{\bf x}_s,\omega)$ represents the observed data. It can be written in a more compact notation as follows:
\new{
\begin{equation}
F({\bf m}) = \frac{1}{2}||{\bf W}({\bf d} - {\bf d}_{\rm obs})||_2^2 = \frac{1}{2}||{\bf W}({{\bf Lm}-{\bf d}_{\rm obs}})||_2^2, \label{eq_objfnc_compact}
\end{equation}
}
where ${\bf L}$ is the forward modeling operator defined in equation \ref{eq_fwd}, ${\bf W}$ is the acquisition mask operator described by equation \ref{eq_mask}, 
and $||\cdot||_2$ stands for the $\ell_2$ norm.
\new{
The data-space objective function $F({\bf m})$ is usually minimized with a gradient-based optimization solver, 
which iterates until an acceptable image is obtained \cite[]{nemeth1999,kuhl2003,marie2005}. 
However, the data-space inversion scheme
lacks flexibility and cannot be implemented in a target-oriented fashion.
Full-domain migration/demigration has to be carried out within each iteration;
and the optimization converges slowly without a proper preconditioner.
Therefore, the data-space inversion scheme is computationally challenging for large-scale applications.
}

\par
Alternatively, we can reformulate the inverse problem and solve it in the model space.
Because ${\bf L}$ is a linear operator, $F({\bf m})$ is a quadratic function. 
$F({\bf m})$ reaches its minimum when ${\bf m}$ satisfies \cite[]{tarantola1987}:
\begin{equation}
{\bf m} = {\bf H}^{-1}{\bf L}^{*}{\bf W}^{*}{\bf d}_{\rm obs}, \label{eq_minimum}
\end{equation}
where $^{*}$ means adjoint (complex conjugation and transpose)
and ${\bf H}={\bf L}^{*}{\bf W}^{*}{\bf W}{\bf L}$ is the weighted Hessian operator.
\new{Equation \ref{eq_minimum} has only symbolic meaning, because the Hessian is often
singular and its inverse is not easy to obtain directly.} A more practical way would be 
reconstructing the reflectivity ${\bf m}$ through iterative inverse filtering by
minimizing a model-space objective function defined as follows \cite[]{valenciano2006}:
\begin{equation}
J({\bf m}) = \frac{1}{2}|| {\bf Hm} - {\bf m}_{\rm mig} ||_2^2, \label{eq_objfnc_modl}
\end{equation}
where ${\bf m}_{\rm mig}={\bf L}^{*}{\bf W}^{*}{\bf d}_{\rm obs}$ is the migrated image.

One advantage of the model-space formulation is 
that it can be implemented in a target-oriented fashion, 
which can substantially reduce the size of the problem and hence the computational cost \cite[]{valenciano2006,valenciano2008}. 
For example, we can choose to invert only areas of particular interest, 
such as subsalt regions, where potential reservoirs are located and migration often fails to provide reliable images.
In fact, the model-space approach divides the inversion problem 
into two stages: computing the explicit Hessian (only once) and minimizing the objective function $J({\bf m})$. 
Since the cost of minimizing the objective function $J({\bf m})$ is relatively small once the explicit Hessian is obtained, 
various regularization schemes incorporating different {\it a priori} information can be easily tested at virtually little cost.

\par
Each element of the Hessian is given as follows by taking the second derivative of the objective function $F({\bf m})$ with respect to 
the model parameters \cite[]{plessix2004,valenciano2006}; 
note that $w^2({\bf x}_r,{\bf x}_s)=w^{*}({\bf x}_r,{\bf x}_s)=w({\bf x}_r,{\bf x}_s)$:
\begin{eqnarray}
& & H({\bf x},{\bf y})   \nonumber \\
& & = \Re \left\{ \right.\sum_{\omega} \omega^4 \sum_{{\bf x}_s} |f_s(\omega)|^2
       G({\bf x},{\bf x}_s,\omega)G^{*}({\bf y},{\bf x}_s,\omega) \nonumber \\
& & \times \sum_{{\bf x}_r} w({\bf x}_r,{\bf x}_s) G({\bf x},{\bf x}_r,\omega)G^{*}({\bf y},{\bf x}_r,\omega)\left. \right\}, \label{eq_hessian_exact}
\end{eqnarray}
where $\Re$ denotes taking the real part of a complex value, 
and ${\bf y}$ is a neighbor point around the image point ${\bf x}$ in the subsurface.
\new{
When ${\bf x}={\bf y}$, we obtain the diagonal elements of the Hessian; when ${\bf x}\neq{\bf y}$, we obtain the off-diagonal elements.
A target-oriented truncated Hessian is obtained by computing the Hessian for ${\bf x}$'s that are within the target zone and
a small number of ${\bf y}$'s that are close to each ${\bf x}$ \cite[]{valenciano2006,valenciano2008}.
}

\par
Hereafter, I call $H({\bf x}, {\bf y})$ in equation \ref{eq_hessian_exact} the exact Hessian, 
since it is derived strictly from the least-squares functional $F({\bf m})$ in the shot-profile domain.
\new{
However, direct implementation of equation \ref{eq_hessian_exact} requires saving Green's functions, which
may bring considerable computational issues, because the Green's functions can be huge for practical
applications, especially in 3-D.
To reduce the computational overburden, in the subsequent sections, 
I introduce an alternative method based on phase encoding to 
compute the Hessian operator. 
As I will demonstrate, by using this approach, we do not need to save any
Green's functions, and the cost is also significantly reduced.
}

\section{phase-encoded Hessian}

\subsection{Encoding of the receiver-side Green's functions}
\par
Suppose we have a new operator $\widetilde{H}({\bf x},{\bf y},{\bf p}_r)$ defined as follows:
\begin{eqnarray}
& &\widetilde{H}({\bf x}, {\bf y}, {\bf p}_r) \nonumber \\
& &= \Re \left\{ \right. \sum_{\omega} \omega^4 \sum_{{\bf x}_s}
|f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega)G^{*}({\bf y},{\bf x}_s,\omega) \nonumber\\
& & \times \sum_{{\bf x}_r}w({\bf x}_r,{\bf x}_s)
    G({\bf x},{\bf x}_r,\omega)\alpha({\bf x}_r,{\bf p}_r,\omega) \nonumber \\
& & \times \sum_{{\bf x}_r^{\prime}}w({\bf x}_r^{\prime},{\bf x}_s)
    G^{*}({\bf y},{\bf x}_r^{\prime},\omega)\alpha^{*}({\bf x}_r^{\prime},{\bf p}_r,\omega) \left. \right\}, \label{eq_hessian_approx}
\end{eqnarray}
where we introduce an extra summation $\sum_{{\bf x}_r^{\prime}}$ for the receiver-side Green's functions,
and $\alpha({\bf x}_r,{\bf p}_r,\omega)$ 
is some weighting function, \new{or more specifically, encoding function,} to be specified later. 
\new{Hereafter, $\widetilde{H}({\bf x},{\bf y},{\bf p}_r)$ is referred as the receiver-side phase-encoded Hessian.}
Though $\widetilde{H}({\bf x},{\bf y},{\bf p}_r)$ can 
be computed in a way similar to the direct implementation of the exact Hessian as defined by equation \ref{eq_hessian_exact},
it offers more flexibility and can be very efficiently computed without 
explicitly saving the Green's functions.
\par
Because of the linearity of the wave equation, 
the term 
\small
$\sum_{{\bf x}_r} w({\bf x}_r,{\bf x}_s)G({\bf x},{\bf x}_r,\omega)\alpha({\bf x}_r,{\bf p}_r,\omega)$  
\normalsize
now becomes the extrapolated wavefield at image point ${\bf x}$. 
It corresponds to the composite source $f_c(x,y,{\bf x}_s,{\bf p}_r,\omega)$ on the surface that can be defined as follows:
\begin{eqnarray}
& & f_{c}(x,y,{\bf x}_s,{\bf p}_r,{\omega}) \nonumber \\
& & = \sum_{{\bf x}_r}w({\bf x}_r,{\bf x}_s)\delta(x-x_r,y-y_r)\alpha({\bf x}_r,{\bf p}_r,\omega),
\end{eqnarray}
where $\delta(x-x_r,y-y_r)$ is the Dirac delta function centered at ${\bf x}_r=(x_r,y_r)$ and can be considered as a point source on the surface. 
Hence, the composite source is obtained by linearly combining
point sources at all receiver locations for a specific shot located at ${\bf x}_s$. 
The function $\alpha({\bf x}_r,{\bf p}_r,\omega)$ serves
as a weight when combining these point sources. 
The same thing holds for the other summation term, 
\small
$\sum_{{\bf x}_r^{\prime}}w({\bf x}_r^{\prime},{\bf x}_s)G^{*}({\bf x}_r^{\prime},{\bf y},\omega)\alpha^{*}({\bf x}_r^{\prime},\omega)$
\normalsize
, except that
it is the complex conjugate of the same extrapolated wavefield at a neighbor image point ${\bf y}$.
To make it clearer, we define a receiver wavefield $R({\bf x},{\bf x}_s,{\bf p}_r,\omega)$ corresponding to 
the receiver composite source $f_c(x,y,{\bf x}_s,{\bf p}_r,\omega)$ as
\begin{eqnarray} 
R({\bf x},{\bf x}_s,{\bf p}_r,\omega) = \sum_{{\bf x}_r}w({\bf x}_r,{\bf x}_s)G({\bf x}_r,{\bf x},\omega)\alpha({\bf x}_r,{\bf p}_r,\omega), \label{eq_r}
\end{eqnarray}
and a source wavefield corresponding to a point source located at ${\bf x}_s$ on the surface as
\begin{equation}
S({\bf x},{\bf x}_s,\omega) = f_s(\omega)G({\bf x}, {\bf x}_s, \omega). \label{eq_s}
\end{equation}
Substituting equations \ref{eq_r} and \ref{eq_s} into equation \ref{eq_hessian_approx} leads to
\begin{eqnarray}
\widetilde{H}({\bf x},{\bf y},{\bf p}_r) &=& \Re \left\{ \right. \sum_{\omega} \omega^4 \sum_{{\bf x}_s} 
S({\bf x},{\bf x}_s,\omega) S^{*}({\bf y},{\bf x}_s,\omega) \nonumber \\
&\times& R({\bf x},{\bf x}_s,{\bf p}_r,\omega) R^{*}({\bf y},{\bf x}_s,{\bf p}_r,\omega)\left. \right\}, \label{eq_hessian_approx1}
\end{eqnarray}
which means $\widetilde{H}({\bf x},{\bf y},{\bf p}_r)$ can be computed by 
crosscorrelating the source and receiver wavefields with their shifted complex conjugates 
(${\bf y}$ is a neighborhood point around ${\bf x}$).
Therefore, the algorithm for computing the receiver-side phase-encoded Hessian 
is similar to the wave-equation shot-profile migration algorithm, except that the boundary condition for the receiver wavefield 
and the imaging condition are slightly modified. In the shot-profile migration, the boundary condition for the receiver wavefield is the recorded data,
and we crosscorrelate the source wavefield and the receiver wavefield to produce the image; here, in contrast, we use the composite source as the boundary condition
for the receiver wavefield, and then invoke the imaging condition defined by equation \ref{eq_hessian_approx1} 
to produce the new operator $\widetilde{H}({\bf x},{\bf y},{\bf p}_r)$.
In other words, we do not have to save the Green's functions at all. 
Computing $\widetilde{H}({\bf x},{\bf y},{\bf p}_r)$ is also efficient because multiple Green's functions are computed at 
the same time during the wavefield extrapolation.

\new{
We can stack $\widetilde{H}({\bf x},{\bf y},{\bf p}_r)$ over ${\bf p}_r$ and after some simple algebraic manipulation, we obtain:
\begin{eqnarray}
& &\widetilde{H}({\bf x},{\bf y})=\sum_{{\bf p}_r}\widetilde{H}({\bf x},{\bf y},{\bf p}_r) \nonumber \\
& &= \Re \left\{ \right. \sum_{\omega} \omega^4
    \sum_{{\bf x}_s} |f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega) G^{*}({\bf y},{\bf x}_s,\omega) \nonumber \\
& & \times \sum_{{\bf x}_r}w({\bf x}_r,{\bf x}_s) G({\bf x},{\bf x}_r,\omega)G^{*}({\bf y},{\bf x}_r,\omega)
    \sum_{{\bf p}_r}|\alpha({\bf x}_r,{\bf p}_r,\omega)|^2 \nonumber \\
& & + \sum_{\omega} \omega^4 \sum_{{\bf x}_s}
      |f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega) G^{*}({\bf y},{\bf x}_s,\omega) \nonumber \\
& & \times \sum_{{\bf x}_r} \sum_{{\bf x}_r'\not={\bf x}_r}
           w({\bf x}_r,{\bf x}_s)G({\bf x},{\bf x}_r,\omega) w({\bf x}_r',{\bf x}_s)G^{*}({\bf y},{\bf x}_r',\omega) \nonumber \\
& & \times \sum_{{\bf p}_r}\alpha({\bf x}_r,{\bf p}_r,\omega)\alpha^{*}({\bf x}_r',{\bf p}_r,\omega) \left. \right\}.
 \label{eq_hessian_approx3}
\end{eqnarray}
}
\par
If we let the weighting function satisfy $\sum_{{\bf p}_r}|\alpha({\bf x}_r,{\bf p}_r,\omega)|^2=1$, the first term 
in equation \ref{eq_hessian_approx3} becomes the exact Hessian $H({\bf x}, {\bf y})$ (equation \ref{eq_hessian_exact}); 
however, if the second term cannot be removed, it becomes
undesired crosstalk from the crosscorrelations among different receiver-side Green's functions and we obtain:
\new{
\begin{eqnarray}
\widetilde{H}({\bf x},{\bf y}) = H({\bf x},{\bf y}) + {\rm Crosstalk}. \label{eq_hessian_approx4}
\end{eqnarray}
}
%\begin{eqnarray}
%& &\widetilde{H}({\bf x}, {\bf y}) = H({\bf x},{\bf y}) \nonumber \\
%& & + \Re \left\{ \right. \sum_{\omega} \omega^4 \sum_{{\bf x}_s}
%      |f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega) G^{*}({\bf y},{\bf x}_s,\omega) \nonumber \\
%& & \times \sum_{{\bf x}_r} \sum_{{\bf x}_r'\not={\bf x}_r}           
%    w({\bf x}_r,{\bf x}_s)G({\bf x},{\bf x}_r,\omega) w({\bf x}_r',{\bf x}_s)G^{*}({\bf y},{\bf x}_r',\omega) \nonumber \\
%& & \times \sum_{{\bf p}_r}\alpha({\bf x}_r,{\bf p}_r,\omega)\alpha^{*}({\bf x}_r',{\bf p}_r,\omega) \left. \right\}.
% \label{eq_hessian_approx4}
%\end{eqnarray}

\subsection{Simultaneous encoding of the source- and receiver-side Green's functions}
We can further encode the source-side Green's function by synthesizing composite sources from the source locations.
For simplicity, we assume OBS or land acquisition geometries, where all shots have the same receiver spread. 
Therefore, we have the following relation:
\begin{eqnarray}
w({\bf x}_r,{\bf x}_s) = w_r({\bf x}_r)w_s({\bf x}_s),
\end{eqnarray}
where $w_r({\bf x}_r)$ and $w_s({\bf x}_s)$ define the ranges of receivers and sources for a given acquisition geometry.
With the above assumption, \new{we can construct a simultaneously phase-encoded Hessian as follows}:
\begin{eqnarray}
& &\widetilde{\widetilde{H}}({\bf x},{\bf y},{\bf p}_s,{\bf p}_r)
         = \Re \left\{ \right. \sum_{\omega}\omega^4 \nonumber \\
& & \times \sum_{{\bf x}_s}        w_s({\bf x}_s)         f_s(\omega) G({\bf x},{\bf x}_s,\omega)          \beta({\bf x}_s,{\bf p}_s,\omega) \nonumber \\
& & \times \sum_{{\bf x}_s^{\prime}} w_s({\bf x}_s^{\prime})f_s^*(\omega) G^*({\bf y},{\bf x}_s^{\prime},\omega) \beta^*({\bf x}_s^{\prime},{\bf p}_s,\omega) \nonumber\\
& & \times \sum_{{\bf x}_r}          w_r({\bf x}_r)         G({\bf x},{\bf x}_r         ,\omega) \alpha({\bf x}_r,{\bf p}_r,\omega) \nonumber \\
& & \times \sum_{{\bf x}_r^{\prime}} w_r({\bf x}_r^{\prime})G^*({\bf y},{\bf x}_r^{\prime},\omega) \alpha^*({\bf x}_r^{\prime},{\bf p}_r,\omega) \left. \right\}.
 \label{eq_hessian_both}
\end{eqnarray}
where we introduce two extra summations: $\sum_{{\bf x}_s^{\prime}}$ for the source-side Green's functions and 
$\sum_{{\bf x}_r^{\prime}}$ for the receiver-side Green's functions. 
Similar to $\alpha({\bf x}_r,{\bf p}_r,\omega)$, $\beta({\bf x}_s,{\bf p}_s,\omega)$ is also a weighting function and will be specified later.
Let us once again define the composite source wavefield $S({\bf x},{\bf p}_s,\omega)$ and composite receiver wavefield $R({\bf x},{\bf p}_r,\omega)$ as follows:
\begin{equation}
S({\bf x},{\bf p}_s,\omega) = \sum_{{\bf x}_s}w_s({\bf x}_s)f_s(\omega)G({\bf x},{\bf x}_s,\omega)\beta({\bf x}_s,{\bf p}_s,\omega), \label{eq_comp_sou_obc}
\end{equation}
and
\begin{equation}
R({\bf x},{\bf p}_r,\omega) = \sum_{{\bf x}_r}w_r({\bf x}_r)G({\bf x},{\bf x}_r,\omega)\alpha({\bf x}_r,{\bf p}_r,\omega). \label{eq_comp_rec_obc}
\end{equation}
Substituting equations \ref{eq_comp_sou_obc} and \ref{eq_comp_rec_obc} into equation \ref{eq_hessian_both} leads to
\begin{eqnarray}
\widetilde{\widetilde{H}}({\bf x},{\bf y},{\bf p}_s,{\bf p}_r)
&=& \Re \left\{ \right. \sum_{\omega}\omega^4S({\bf x},{\bf p}_s,\omega)S^{*}({\bf y},{\bf p}_s,\omega) \nonumber \\
&\times& R({\bf x},{\bf p}_r,\omega)R^*({\bf y}, {\bf p}_r,\omega) \left. \right\}. \label{eq_hessian_both2}
\end{eqnarray}
Therefore, computing the simultaneously phase-encoded Hessian $\widetilde{\widetilde{H}}({\bf x},{\bf y},{\bf p}_s,{\bf p}_r)$ 
requires only two wavefield propagations: 
one for the composite source wavefield defined by equation \ref{eq_comp_sou_obc}, and the other for the composite receiver wavefield defined by
equation \ref{eq_comp_rec_obc}. 

\new{
Similar to the case of the receiver-side phase-encoded Hessian, 
we can stack $\widetilde{\widetilde{H}}({\bf x},{\bf y},{\bf p}_s,{\bf p}_r)$ over both ${\bf p}_s$ and ${\bf p}_r$. 
If the weighting functions satisfy  
$\sum_{{\bf p}_r}|\alpha({\bf x}_r,{\bf p}_r,\omega)|^2=1$ and $\sum_{{\bf p}_s}|\beta({\bf x}_s,{\bf p}_s,\omega)|^2=1$, 
the stacked result $\widetilde{\widetilde{H}}({\bf x},{\bf y})$ becomes the sum of the exact Hessian and the crosstalk: 
\begin{eqnarray}
\widetilde{\widetilde{H}}({\bf x},{\bf y}) &=& \sum_{{\bf p}_s}\sum_{{\bf p}_r}\widetilde{\widetilde{H}}({\bf x},{\bf y},{\bf p}_s,{\bf p}_r) \nonumber \\
&=& H({\bf x},{\bf y}) + {\rm Crosstalk}. \label{eq_hessian_both4}
\end{eqnarray}
}

Thus, by encoding the Green's functions for the Hessian computation, 
we face a situation similar to that encountered in phase-encoding migration \cite[]{romero2000,liu2006}, i.e., our exact Hessian is contaminated by
crosstalk artifacts, so we seek to define weighting functions $\alpha({\bf x}_r,{\bf p}_r,\omega)$ and $\beta({\bf x}_s,{\bf p}_s,\omega)$
that attenuate the crosstalk as much as possible. In the next two sections, I examine two different phase-encoding schemes
to attenuate the crosstalk, namely, plane-wave-phase encoding and random-phase encoding.

\subsection{Plane-wave-phase encoding}
\par
Suppose we choose the weighting functions to be
\begin{eqnarray}
\alpha({\bf x}_r,{\bf p}_r,\omega) &=& A_r(\omega){\rm e}^{i \omega {\bf p}_r \cdot {\bf x}_r}, \\ 
\beta ({\bf x}_s,{\bf p}_s,\omega) &=& A_s(\omega){\rm e}^{i \omega {\bf p}_s \cdot {\bf x}_s},\label{eq_wght_rec}
\end{eqnarray}
which are the well known plane-wave or delayed-shot phase-encoding functions \cite[]{whitmore1995,zhang2005,liu2006}, 
where $i = \sqrt{-1}$, $A_r(\omega)$ and $A_s(\omega)$ are real functions depending upon the angular frequency $\omega$;
and ${\bf p}_r=(p_{r_x},p_{r_y})$ and ${\bf p}_s=(p_{s_x},p_{s_y})$ are the ray parameters for the receiver and source plane waves on the surface, respectively. 
As proved by \cite{liu2006}, in continuous case, integrating over plane waves
completely attenuates the crosstalk in the plane-wave source migration, and the final migration result
is exactly equivalent to the standard shot-profile migration result. 
\new{
The same property holds here in the scenario of Hessian computation, as proved in Appendix A and B: 
integrating over ${\bf p}_r$ and ${\bf p}_s$, and 
choosing $A_r(\omega)$ and $A_s(\omega)$ to satisfy $A_r^2(\omega)=|\omega|$, $A_s^2(\omega)=|\omega|$ in 2-D 
and $A_r^2(\omega)=|\omega|^2$, $A_s^2(\omega)=|\omega|^2$ in 3-D, asymptotically, 
we have $\sum_{{\bf p}_r}\alpha({\bf x}_r,{\bf p}_r,\omega)\alpha^{*}({\bf x}_r',{\bf p}_r,\omega)$
and $\sum_{{\bf p}_s}\beta({\bf x}_s,{\bf p}_s,\omega)\beta^{*}({\bf x}_s',{\bf p}_s,\omega)$ equal $1$ 
(if ${\bf x}_r={\bf x}_r'$, ${\bf x}_s={\bf x}_s'$) and $0$ (if ${\bf x}_r\neq{\bf x}_r'$, ${\bf x}_s\neq{\bf x}_s'$).
Hence the crosstalk in both equations \ref{eq_hessian_approx4} and \ref{eq_hessian_both4} becomes zero, 
and the approximate Hessians $\widetilde{H}({\bf x},{\bf y})$ and $\widetilde {\widetilde {H}}({\bf x},{\bf y})$ 
converge to the exact Hessian $H({\bf x},{\bf y})$. 
}
In the discrete form used in practice, the number of ${\bf p}_s$ and ${\bf p}_r$
values can be chosen in a way similar to those discussed by \cite{stork2004,zhang2005,etgen2005,stork2006,zhang2006}.

\subsection{Random-phase encoding}
Instead of using the plane-wave-phase encoding function, we can use random phases to disperse unwanted crossterms \cite[]{romero2000}. 
The weighting functions can be chosen as follows:
\new{
\begin{eqnarray}
\alpha({\bf x}_r,{\bf p}_r,\omega) &=& \frac{1}{\sqrt{N_{\rm realize}}}{\rm e}^{i \gamma({\bf x}_r,{\bf p}_r,\omega)}, \\
\beta ({\bf x}_s,{\bf p}_s,\omega) &=& \frac{1}{\sqrt{N_{\rm realize}}}{\rm e}^{i \gamma({\bf x}_s,{\bf p}_s,\omega)},
\end{eqnarray}
}
where the phase functions $\gamma({\bf x}_r,{\bf p}_r,\omega)$ and $\gamma({\bf x}_s,{\bf p}_s,\omega)$ are sequences of random
numbers between $0$ and $2\pi$, they are random functions of both space and frequency; the parameters ${\bf p}_r$ and ${\bf p}_s$ 
now define the indexes of different realizations of the random sequences; $N_{\rm realize}$ is the number of realizations.
Obviously, the weighting functions satisfy $\sum_{{\bf p}_r}|\alpha({\bf x}_r,{\bf p}_r,\omega)|^2=1$ 
and $\sum_{{\bf p}_s}|\beta({\bf x}_s,{\bf p}_s,\omega)|^2=1$, therefore, 
the encoded Hessian has the same form as in equation \ref{eq_hessian_approx4} (receiver-side encoded Hessian) 
or \ref{eq_hessian_both4} (simultaneously encoded Hessian).
Because the phases of the crosstalk are randomized, when we sum over frequencies, sources and receivers to generate the final result, 
the phases of the crosstalk will not agree with each other, and consequently, the crosstalk will be attenuated by stacking. To maximize
the phase differences, a uniformly distributed random sequence could be used \cite[]{romero2000}. To further attenuate the crosstalk, 
multiple realizations of the random sequences can be used ($N_{\rm realize}>1$): the randomly phase-encoded Hessian is computed multiple times with 
different realizations of the random-phase function
and then stacked together to produce the final result. 
The cost of this process is linearly proportional to the number of realizations. 
 

\subsection{Cost comparison}
\par
In this section, I compare the cost for different methods to determine the savings generated by 
using the phase-encoding method. As discussed before, Hessian computation contains two main parts: wavefield propagation 
(Green's functions) and crosscorrelation among different Green's functions. Since the crosscorrelation parts are similar for methods 
with or without phase encoding, I will only compare the cost for the first part, i.e, wavefield propagation for Green's functions.
\new{I consider general 3-D cases; for 2-D applications, the cost for each method can also be conviniently obtained by setting
the corresponding crossline parameter(s) (with subscript $_y$) to $1$.
}
\par
Let us assume a 3-D seismic survey that has $N_s$ shots in total, covering a surface area that can be divided into $N_x \times N_y$ cells; 
and our target area for inversion is a 3-D cube with $M_x$, $M_y$ and $M_z$ samples in $x$, $y$ and $z$ directions. 
Table \ref{tbl:tb1} illustrates the cost comparison for a marine acquisition geometry, while Table \ref{tbl:tb2} compares the
cost for an OBS or a land acquisition geometry. Since simultaneous phase encoding is derived under an OBS or a land acquisition geometry, 
Table \ref{tbl:tb1} only
compares the cost of the receiver-side encoding method with that of the direct method. In both tables,
$N_{\omega}$ is the number of frequencies of the wavefield, $N_{p_{s_x}}$ and $N_{p_{s_y}}$ are the numbers of source ray parameters for the inline
and crossline directions, and $N_{p_{r_x}}$ and $N_{p_{r_y}}$ are the numbers of receiver ray parameters for the inline
and crossline directions, respectively. $N_{\rm realize}$ is the number of realizations of the random sequences, 
the same as that defined in the previous section. 

\tabl{tb1}{Cost comparison for a marine acquisition geometry.}
{
\begin{center}
\begin{tabular}{|l|l|l|}
\hline
Method & Number of wavefield propagations & Number of Green's functions to store\\
\hline
direct computation & $N_x N_y N_{\omega}$ & $M_xM_yM_zN_xN_yN_{\omega}$ \\
plane-wave, receiver-side & $(1+N_{p_{r_x}}N_{p_{r_y}})N_sN_{\omega}$ & $0$      \\
random, receiver-side & $(1+N_{\rm realize})N_sN_{\omega}$ & $0$ \\
\hline
\end{tabular}
\end{center}
}

\tabl{tb2}{Cost comparison for an OBS or a land acquisition geometry.}
{
\begin{center}
\begin{tabular}{|l|l|l|}
\hline
Method & Number of wavefield propagations & Number of Green's functions to store\\
\hline
direct computation & $N_x N_y N_{\omega}$ & $M_xM_yM_zN_xN_yN_{\omega}$ \\
plane-wave, receiver-side & $(1+N_{p_{r_x}}N_{p_{r_y}})N_sN_{\omega}$ & $0$      \\
plane-wave, simultaneous  & $N_{p_{s_x}}N_{p_{s_y}}N_{p_{r_x}}N_{p_{r_y}}N_{\omega}$ & $0$ \\
random, receiver-side & $(1+N_{\rm realize})N_sN_{\omega}$ & $0$ \\
random, simultaneous  & $2N_{\rm realize}N_{\omega}$ & $0$ \\
\hline
\end{tabular}
\end{center}
}

\par
\new{
From both tables, we can find that 
the simultaneous random-phase encoding is the most efficient method for an OBS or a land acquisition geometry, 
which requires only $2N_{\omega}$ wavefield propagations for one realization. 
If a marine acquisition geometry is used, 
it would be efficient to use the receiver-side random-phase encoding (if $N_s<N_xN_y$); 
the cost for one realization is comparable to that of a shot-profile migration. 
Plane-wave phase-encoding methods, however, may potentially need more wavefiled propagations than
the direct method, depending on how many plane waves are used for computing the Hessian. 
Also note that the phase-encoding methods do not require storage of any Green's functions, which
is a crucial benefit for large-scale applications, since the size of the Green's functions ($M_x M_y M_z N_x N_y N_{\omega}$) can easily 
reach an unaffordably large number. For example, for a target area with size $M_x=M_y=M_z=100$ and survey area with size $N_x=N_y=1000$, the
Green's functions for a single frequency ($N_{\omega}=1$) is about $8$ terabytes ($8$ comes from the fact that the Green's function is a complex value).
}

\section{numerical examples}
\par
In this section, I show several numerical examples for two different velocity models. The first is a simple constant-velocity model, which I use
to verify the proposed algorithms for computing the Hessian; 
the second is the more complicated Sigsbee2A velocity model \cite[]{paffenholz2002}. Because of the complex salt body and limited acquisition geometry, 
there are shadow zones under the salt where conventional wave-equation migration algorithms often fail to produce reliable images. 
Therefore, it is useful to test if inversion by using the phase-encoded Hessian can improve the final image.

\subsection{Verification of the algorithm: a constant velocity model}
\par
I verify the proposed phase-encoding algorithms on a constant velocity model ($v=2$ km/s) for several different acquisition geometries.
Figure \ref{fig:const-hess-diag-1s2r} shows the diagonal part of the Hessian (when ${\bf x} = {\bf y}$) for an acquisition geometry containing
one shot at $-0.6$ km and two receivers located at $0.6$ km and $1.2$ km on the surface.
A ricker wavelet with a $20$ Hz dominant frequency is used as the source function, and 
frequencies from $5$ Hz to $35$ Hz are used to generate the following results.
Figure \ref{fig:const-hess-diag-1s2r}(a) is the exact diagonal of the Hessian, uncontaminated by any crosstalk artifacts. 
Figure \ref{fig:const-hess-diag-1s2r}(b) shows the Hessian with crosstalk, which is 
obtained using equation \ref{eq_hessian_approx} with $\alpha({\bf x}_r,{\bf p}_r,\omega)=1$; 
the crosstalk can be easily identified as the vertical stripes on the right side of the image.
Figure \ref{fig:const-hess-diag-1s2r}(c) is the receiver-side plane-wave phase-encoded Hessian with $31$ plane waves stacked together.
The crosstalk has been successfully removed, and Figure \ref{fig:const-hess-diag-1s2r}(c) looks 
very similar to Figure \ref{fig:const-hess-diag-1s2r}(a).
Figure \ref{fig:const-hess-diag-1s2r}(d) shows the receiver-side randomly phase-encoded Hessian with one random realization; 
the crosstalk is dispersed throughout the image.
Figure \ref{fig:const-hess-diag-1s2r}(e) and Figure \ref{fig:const-hess-diag-1s2r}(f) show the results
for $5$ and $20$ random realizations; as expected, the crosstalk is better attenuated by using more random realizations.

Figure \ref{fig:const-hess-offd-1s2r} shows the Hessian with off-diagonals (with size $81\times81$) at the image point $x=0.5$ km, $z=1.5$ km. 
The acquisition geometry is the same as that in Figure \ref{fig:const-hess-diag-1s2r}. 
The horizontal and vertical axes in Figure \ref{fig:const-hess-offd-1s2r} denote the horizontal and 
vertical offsets away from the image point $x=0.5$ km, $z=1.5$ km. Figure \ref{fig:const-hess-offd-1s2r}(a) shows the exact Hessian.
Since the illumination in the subsurface is very poor (only one source and two receivers), we can observe two events intersecting at the origin 
\new{(the diagonal element of the Hessian)} and 
spreading out over the off-diagonals. When more sources and receivers are used, i.e., with better subsurface illumination, 
the energy in the local Hessian operator would be more focused at the origin (as illustrated in the subsequent examples).
Figure \ref{fig:const-hess-offd-1s2r}(b) is the Hessian with crosstalk, where two extra undesired events \new{present} in the operator; 
Figure \ref{fig:const-hess-offd-1s2r}(c) shows the receiver-side plane-wave phase-encoded Hessian, plane-wave-phase encoding successfully attenuates the crosstalk; 
Figures \ref{fig:const-hess-offd-1s2r}(d), (e) and (f) show the receiver-side randomly phase-encoded Hessian with different number of realizations.

\plot{const-hess-diag-1s2r}{width=\columnwidth}
{Diagonal of the Hessian for a constant-velocity model with one shot and two receivers. 
(a) The exact diagonal of the Hessian (equation \ref{eq_hessian_exact});
(b) the Hessian contaminated by crosstalk (equation \ref{eq_hessian_approx} with $\alpha({\bf x}_r,{\bf p}_r,\omega)=1$); 
(c) the receiver-side plane-wave phase-encoded Hessian;
(d), (e) and (f) are the receiver-side randomly phase-encoded Hessians obtained with $1$, $5$ and $20$ random realizations, respectively.}

\plot{const-hess-offd-1s2r}{width=\columnwidth}
{The local Hessian operator for image point $x=0.5$ km, $z=1.5$ km (a row of the Hessian matrix). 
 The acquisition geometry is the same as that in Figure \ref{fig:const-hess-diag-1s2r}. 
 The size of the Hessian operator is $81$ samples in both $x$ and $z$ directions.
 (a) The exact Hessian; 
 (b) The Hessian contaminated by crosstalk;
 (c) the receiver-side plane-wave phase-encoded Hessian;
 (d), (e) and (f) are the receiver-side randomly phase-encoded Hessian obtained with $1$, $5$ and $20$ random realizations, respectively.}

In the next example, I change the acquisition geometry used in the previous example from 
$2$ receivers to $401$ receivers. The receivers 
range from $-2$ km to $2$ km, with a spacing of $0.01$ km. The source location is changed to $-1$ km. 
Once again, Figure \ref{fig:const-hess-diag-1s401r} shows the diagonal part of the Hessian operator, 
while Figure \ref{fig:const-hess-offd-1s401r} illustrates the off-diagonals at 
the image point $x=0.5$ km, $z=1.5$ km in the subsurface.
Since the receiver spread is asymmetric with respect to the source location, the energy in the exact diagonal of Hessian 
(Figure \ref{fig:const-hess-diag-1s401r}(a)) is slightly tilted towards the right and the energy in the local Hessian operator shown 
in Figure \ref{fig:const-hess-offd-1s401r}(a) is well focused due to contributions from many receivers. 
However, when the Hessian is contaminated by crosstalk, 
the energy in the diagonal of the Hessian (Figure \ref{fig:const-hess-diag-1s401r}(b)) shows wrong orientation 
and the illumination pattern is very different from the correct one (Figure \ref{fig:const-hess-diag-1s401r}(a)); 
the local Hessian operator shown in Figure \ref{fig:const-hess-offd-1s401r}(b) is also far from focusing.
The crosstalk is successfully attenuated by the proposed phase-encoding techniques, as shown in Figure \ref{fig:const-hess-diag-1s401r}(c) 
and Figure \ref{fig:const-hess-offd-1s401r}(c) for the receiver-side plane-wave-phase encoding ($31$ receiver-side plane waves are used), and 
Figure \ref{fig:const-hess-diag-1s401r}(d) and Figure \ref{fig:const-hess-offd-1s401r}(d) 
for the receiver-side random-phase encoding with one realization. Although the crosstalk is greatly suppressed in 
Figure \ref{fig:const-hess-diag-1s401r}(d) and Figure \ref{fig:const-hess-offd-1s401r}(d), some random noise appears in the background, 
which is caused by the random nature of the phase-encoding function. The random noise can be further suppressed by using more random realizations.
Figure \ref{fig:const-hess-diag-1s401r}(e), Figure \ref{fig:const-hess-offd-1s401r}(e) 
and Figure \ref{fig:const-hess-diag-1s401r}(f), Figure \ref{fig:const-hess-offd-1s401r}(f) show the results
with $5$ and $20$ random realizations, respectively. The crosstalk is greatly attenuated.

Instead of computing more random realizations, which may substantially increase the cost, 
we can use the fact that the randomized crosstalk can be well attenuated 
when many sources are used (this is almost always true in practice, where
thousands of shots are often fired in acquiring field seismic data). 
Because the phase function is random both in space and frequency, stacking over more sources makes the phases of the crosstalk
more inconsistent with each other, hence better suppresses these artifacts. 
Figure \ref{fig:const-hess-diag-401s401r} and Figure \ref{fig:const-hess-offd-401s401r} show examples of many sources
are used for the same constant velocity model. In this case, a land acquisition geometry is used, where the receiver locations are fixed for all shots.
The receiver location is the same as in the previous example, but now we have $401$ shots from $-2$ km to $2$ km with a $0.01$ km spacing.
Figure \ref{fig:const-hess-diag-401s401r}(a) shows the exact diagonal of the Hessian, while Figure \ref{fig:const-hess-diag-401s401r}(b)
shows the diagonal of the receiver-side randomly phase-encoded Hessian with only one realization. As apposed to the previous example 
(Figure \ref{fig:const-hess-diag-1s401r}(d)), 
the crosstalk is greatly attenuated even
with one random realization. Figure \ref{fig:const-hess-offd-401s401r} compares the off-diagonals of the Hessian for the image point located at
$x=0.5$ km and $z=1.5$ km; as we can see from the comparison, even with only one realization, the encoded Hessian is almost 
identical to the exact Hessian. 

\plot{const-hess-diag-1s401r}{width=\columnwidth}
{Diagonal of the Hessian for a constant-velocity model with one shot and $401$ receivers.
(a) The exact diagonal of the Hessian;
(b) the Hessian contaminated by crosstalk;
(c) the receiver-side plane-wave phase-encoded Hessian;
(d), (e) and (f) are the receiver-side randomly phase-encoded Hessians obtained with $1$, $5$ and $20$ random realizations, respectively.}

\plot{const-hess-offd-1s401r}{width=\columnwidth}
{The local Hessian operator for image point $x=0.5$ km, $z=1.5$ km.
 The acquisition geometry is the same as that in Figure \ref{fig:const-hess-diag-1s401r}.
 (a) The exact Hessian;
 (b) The Hessian contaminated by crosstalk;
 (c) the receiver-side plane-wave phase-encoded Hessian;
 (d), (e) and (f) are the receiver-side randomly phase-encoded Hessian obtained with $1$, $5$ and $20$ random realizations, respectively.}

\plot{const-hess-diag-401s401r}{width=\columnwidth}
{Diagonal of the Hessian with $401$ shots and $401$ receivers. 
(a) The exact diagonal of the Hessian;
(b) the receiver-side randomly phase encoded Hessian with only one random realization.}

\plot{const-hess-offd-401s401r}{width=\columnwidth}
{The local Hessian operator for image point at $x=0.5$ km and $z=1.5$ km.
(a) The exact diagonal of the Hessian;
(b) the receiver-side randomly phase encoded Hessian with only one random realization.}

\subsection{Application of the phase-encoded Hessian: Sigsbee2A model}
\par
To demonstrate an application of the explicit Hessian operator, 
I apply the model-space inversion approach to the Sigsbee2A model \cite[]{paffenholz2002}.
Figure \ref{fig:sigsb2a-vmod} shows the stratigraphic velocity model, which has been used for migration and Hessian computation.
The original data is modeled using two-way wave equation with a marine type acquisition geometry, 
where the receiver spread is moving along with the source \cite[]{paffenholz2002}.
There are $500$ shots in total and the maximum offset for each shot is about $8$ km.
Though the data is modeled with two-way wave equation, in the following examples, one-way wave-equation-based 
Fourier finite-difference (FFD) \cite[]{ristow1994} propagator is used for both migration and phase-encoded Hessian computation.
The explicit Hessian operator is computed using the receiver-side
random-phase encoding method with the frequency band from $5$ Hz to $35$ Hz, equivalent to that of the migrated image ${\bf m}_{\rm mig}$. 
Only one realization has been used; because the number of shots is big, one realization is sufficient to provide reasonably good suppression of the crosstalk. 

I apply two different strategies of applying the explicit Hessian operator: 
the first is to compute only the diagonal of the Hessian operator and use it \new{as a weight} to normalize the migrated image; 
this is often known as the normalized migration \cite[]{rickett2003} or amplitude-preserving migration \cite[]{plessix2004}. 
The other strategy is to compute a limited number of Hessian off-diagonals, and then iteratively minimize the 
model-space objective function $J({\bf m})$ (equation \ref{eq_objfnc_modl})
to find an optimum reconstruction of the reflectivity ${\bf m}$. 
In this example, the number of the off-diagonal elements is chosen by trial and error and each local Hessian filter has $21\times21$ elements, 
the size of which seems to be big enough to capture most of 
the energy even for poorly illuminated areas.
A linear conjugate-gradient solver and a simple damping regularization that minimizes the energy of the model parameters
are used for the model-space iterative inversion.
 

\subsubsection{Normalized migration}
\par
Figure \ref{fig:sigsb2a-hess-diag-random} shows the diagonal of the receiver-side randomly phase-encoded Hessian. 
Note the uneven illumination below the salt caused by the complex salt body and limited acquisition geometry. For comparison, 
Figure \ref{fig:sigsb2a-hess-diag-souint} shows the source-wavefield intensity, 
computed using 
$H_{SI}({\bf x},{\bf x}) 
= \Re \left\{ \right. \sum_{\omega}\omega^4\sum_{{\bf x}_s}|f_s(\omega)G({\bf x},{\bf x}_s,\omega)|^2 \left. \right\}$.
$H_{SI}({\bf x},{\bf x})$ is a crude approximation to the exact diagonal of Hessian, because 
it assumes the constant receiver-side Green's functions and ignores the effects of the limited receiver aperture. It over-estimates 
the total energy that enters the earth and returns to be recorded by the receivers. This is why Figure \ref{fig:sigsb2a-hess-diag-souint} shows 
a stronger, but less accurate, illumination below the salt.

\plot{sigsb2a-vmod}{width=\columnwidth}
{Sigsbee2A stratigraphic velocity model.}

\plot{sigsb2a-hess-diag-random}{width=\columnwidth}
{Diagonal of the Hessian for the Sigsbee2A model. The Hessian is obtained by 
 using receiver-side random-phase encoding,
 which takes the limited receiver array into consideration. The result shows the Hessian
 for a frequency band from $5$ Hz to $35$ Hz.}

\plot{sigsb2a-hess-diag-souint}{width=\columnwidth}
{Source-wavefield intensity for the Sigsbee2A model,
 which assumes the receiver-side Green's functions to be constant;
 it ignores the effects of the limited receiver aperture and over-estimates the total energy that enters
 the earth and returns to be recorded by the receivers.}

\par
Figure \ref{fig:sigsb2a-imag-ffd} 
shows the conventional one-way wave-equation shot-profile migrated image, 
where \new{the unbalanced amplitudes caused by uneven illumination} below the salt are easily identified. 
Figure \ref{fig:sigsb2a-imag-ffd-decon-diag-random} and Figure \ref{fig:sigsb2a-imag-ffd-decon-diag-souint}
show the migrated image normalized with the diagonal of the phase-encoded Hessian (Figure \ref{fig:sigsb2a-hess-diag-random}) and
the source-wavefield intensity (Figure \ref{fig:sigsb2a-hess-diag-souint}), respectively. All three figures are plotted with the same scale.
Figure \ref{fig:sigsb2a-imag-ffd-decon-diag-random} shows more balanced amplitudes across the section than Figure \ref{fig:sigsb2a-imag-ffd-decon-diag-souint}, 
especially in areas below the salt. \new{However, also note that normalized migration does not help improve the
resolution of the image, this is because the diagonal of the Hessian only serves as a scaling factor, it has little deconvolution effects.
The resolution can be improved by using the off-diagonal elements of the Hessian matrix as shown in the subsequent examples.}

%Figure \ref{fig:sigsb2a-amplitude-comparison} compares the normalized amplitude for the horizontal reflector located at $9$ km in depth, 
%normalization with the randomly phase-encoded Hessian shows significant amplitude improvement over migration, 
%and it also provides a slightly better result than normalization with source-wavefield intensity.

\plot{sigsb2a-imag-ffd}{width=\columnwidth}
{Conventional one-way wave-equation shot-profile migration result of the Sigsbee2A model. 
 Note the uneven illumination below the salt.}

\plot{sigsb2a-imag-ffd-decon-diag-random}{width=\columnwidth}
{The migrated image (Figure \ref{fig:sigsb2a-imag-ffd}) is 
 normalized by the diagonal of the Hessian shown in Figure \ref{fig:sigsb2a-hess-diag-random}.}

\plot{sigsb2a-imag-ffd-decon-diag-souint}{width=\columnwidth}
{The migrated image (Figure \ref{fig:sigsb2a-imag-ffd}) is
 normalized by the source-wavefield intensity shown in Figure \ref{fig:sigsb2a-hess-diag-souint}.}

%\plot{sigsb2a-amplitude-comparison}{width=6in}
%{Comparison of the normalized peak amplitude of the horizontal reflector located a $9$ km in depth. 
%Dashed line represents the true amplitude; dotted line is migration; solid line is migration normalized
%with phase-encoded Hessian; dash-dotted line denotes migration normalized with source-wavefield
%intensity.}
\subsubsection{Target-oriented inversion}

\par
The model-space target-oriented inversion is performed to a selected area below the salt. 
Figure \ref{fig:sigsb2a-hess-diag-target,sigsb2a-hess-offd-target} 
shows the receiver-side randomly phase-encoded Hessian operator for the particular region of interest:
Figure \ref{fig:sigsb2a-hess-diag-target} shows the diagonal part of the Hessian operator; 
Figure \ref{fig:sigsb2a-hess-offd-target} is obtained by convolving the Hessian operator ($21\times21$ in size) with a collection of point scatterers. 
It demonstrates the varying shapes and non-stationarities of the Hessian operators across the section. Note that in well-illuminated areas, the 
Hessian operator is well focused, while in poorly illuminated areas, the Hessian operator is tilted and smeared and has a preferred dipping orientation, 
which means these image points are illuminated by only a few dip angles \cite[]{gelius2002}. 
\par 
Figure \ref{fig:sigsb2a-imag-target,sigsb2a-invt-target,sigsb2a-refl-target} 
shows the comparison between migration and inversion: Figure \ref{fig:sigsb2a-imag-target} is the migrated image, 
and Figure \ref{fig:sigsb2a-invt-target} is the inverted image obtained using the randomly phase-encoded Hessian operator. 
The result is obtained after $20$ iterations of the linear conjugate-gradient method. 
Figure \ref{fig:sigsb2a-invt-fmov-win} plots the normalized residual as a function of iteration number; the residual converges after about $12$ iterations.
In the inversion result, \new{the effects of uneven illumination are corrected, the inverted image looks slightly sharper and crisper than the migrated image;
the shadow zones that in the migrated image are now filled in with structures (circled areas), and the sediments (the area outlined by the top circle) 
and the dipping faults (indicated with arrows)
extend much closer to the salt body. The improved image helps locate sediment truncations against salt, which are important in finding
subsalt structural traps}. For comparison, Figure \ref{fig:sigsb2a-refl-target} shows the filtered true reflectivity.

\new{A closeup comparison for areas around the point diffractors is shown in Figure \ref{fig:sigsb2a-closeup}.
Figure \ref{fig:sigsb2a-spectrum} compares the amplitude spectrum of the migrated image (Figure \ref{fig:sigsb2a-closeup}(a))
and the inverted image (Figure \ref{fig:sigsb2a-closeup}(b)); the inverted image shows a broader range of spatial frequencies,
suggesting slightly improved spatial resolution. Figure \ref{fig:sigsb2a-htrace} shows a horizontal trace located 
at depth level $5.1892$ km (slicing through the diffractors)
extracted from
different images shown in Figure \ref{fig:sigsb2a-closeup}. 
The inverted result (the dashed red line) better predicts the true reflectivity (the dotted black line) 
than the migrated result (the blue solid line).}

%\plot{sigsb2a-hess-win}{width=0.65\columnwidth}
\multiplot{2}{sigsb2a-hess-diag-target,sigsb2a-hess-offd-target}{width=0.6\columnwidth}
{Receiver-side randomly phase-encoded Hessian operators for the Sigsbee2A model. 
 Panel (a) shows the diagonal part for a particular region of interest under the salt. 
 Note the uneven illumination due the complex salt body and limited acquisition geometry.
 Panel (b) shows the result obtained by convolving the Hessian operator (with a size $21\times21$) 
 with a collection of point scatterers. 
 Note the non-stationarities of the operators.}

%\plot{sigsb2a-invt-win}{width=0.65\columnwidth}
\multiplot{3}{sigsb2a-imag-target,sigsb2a-invt-target,sigsb2a-refl-target}{width=0.6\columnwidth}
{Comparison between migration and inversion. Panel (a) shows the migrated result 
(b) shows the inversion result using the receiver-side randomly phase-encoded Hessian operator and (c) is the filtered true reflectivity. 
All three panels are plotted with the same scale.}

\plot{sigsb2a-invt-fmov-win}{width=\columnwidth}
{Normalized residual vs. number of iterations for the Sigsbee2A model; the inversion converges after about $12$ iterations.}

%\plot{sigsb2a-vtrace}{width=6in}
%{A vertical trace located at horizontal distance $9.1135$ km extracted from images shown in Figure \ref{fig:sigsb2a-invt-win}. 
%Black dotted line, blue solid line and red dashed line represent traces extracted from 
%the true reflectivity model, the migrated image 
%and the inverted image, respectively.}

\plot{sigsb2a-closeup}{width=\columnwidth}{A closeup view of (a) the migrated image, (b) the inverted image and (c) the filtered true reflectivity model.}

\plot{sigsb2a-spectrum}{width=\columnwidth}{The amplitude spectrum of (a) the migrated image (Figure \ref{fig:sigsb2a-closeup}(a)) 
and (b) the inverted image (Figure \ref{fig:sigsb2a-closeup}(b)). Note the inverted image contains a broader range of spatial frequencies.}

\plot{sigsb2a-htrace}{width=\columnwidth}
{A horizontal trace at depth $5.1892$ km extracted from images shown in Figure \ref{fig:sigsb2a-closeup}.
Black dotted line, blue solid line and red dashed line represent traces extracted from
the true reflectivity model, the migrated image and the inverted image, respectively. }

\section{discussion}
The model-space inversion approach, which requires building the Hessian operator \new{${\bf L}^{*}{\bf L}$}, 
squares the condition number of the Born modeling operator ${\bf L}$. Generally speaking, 
the Hessian is often ill-conditioned, and inverting the Hessian is sensitive to roundoff errors. 
Hence, regularization should be used to stabilize the inversion process. 

The inverted image, Figure \ref{fig:sigsb2a-invt-target}, also shows increased noise, 
for example, at around horizontal distance $11$ km and depth level $4$ km, where the illumination is extremely poor.
This is partially due to the null space of the Hessian operator.
Another possible contribution might be the randomized crosstalk introduced in the randomly phase-encoded 
Hessian. This can be a concern, since the accuracy of the Hessian may affect the convergence of inversion.
However, a small amount of random noise in the operator
might also be useful if one considers it as a random regularization, hence, it may have the chance to improve the condition number
of the Hessian operator, making it easier to invert.
How the crosstalk affects the inversion is still not very clear and is beyond the scope of this paper;
it remains a research area for further investigation.
To suppress the noise, more sophisticated regularization that better predicts the inverse of the model covariance, 
such as roughening along reflection angles \cite[]{kuhl2003,marie2005} or geological dips \cite[]{prucha2001,marie2005},
can be introduced in the inversion process.

\new{
The least-squares inversion approach assumes an accurate background velocity model is known,
which, however, is not always true in practice. Therefore, future research should test
the sensitivity of the inversion approach to velocity errors. 
The theory also assumes single scattering (Born approximation), other coherent energy, such as multiple, does not fit into this framework. 
In the Sigsbee2A example, The internal multiple at around $x=11.6$ km and $z=4.2$ km in the migrated image 
(Figure \ref{fig:sigsb2a-imag-target}) becomes even more pronounced in the inverted image (Figure \ref{fig:sigsb2a-invt-target}).
Hence, multiple suppression is suggested prior to inversion.
}

\section{conclusions}
I present a method based on phase encoding that allows efficient computation of the explicit Hessian operator 
for the wave-equation least-squares migration/inversion problem. 
The proposed algorithm does not require storing Green's functions, which is often too big to be affordable in practice.
Hence, the phase-encoded Hessian has the potential to be applied to large-scale problems at a lower cost.
Although this paper demonstrates examples with one-way wave equation, the theory for computing the phase-encoded Hessian
is not limited to one-way wave equation. Two-way wave equation can also be used in a similar way to compute the phase-encoded Hessian, 
which may bring even more savings in computational cost. 
\new{
The phase-encoded Hessian via two-way wave equation may also be applied as a preconditioner for the full waveform inversion problem.
}
\par
The target-oriented inversion examples on the Sigsbee2A model show that 
wave-equation least-squares migration/inversion
can partially correct the effects of uneven illumination, 
and \new{it can produce an image with more balanced amplitudes and slightly higher spatial resolution to that produced by migration, especially
in areas with low illumination and shadow zones.} 
Thus it provides a good tool for imaging complex geologies.

\section{acknowledgments}
I thank Alejandro Valenciano and Biondo Biondi for enlightening discussions. I thank the sponsors of 
the Stanford Exploration Project for their support, and I acknowledge SMAART JV for providing the Sigsbee2A
data set. I also thank Associate Editor Isabelle Lecomte, Gerrit Toxopeus and two anonymous reviewers, whose
comments and suggestions significantly improve the quality of this paper.

\append{Proof of convergence of the receiver-side plane-wave phase-encoded Hessian}
This appendix demonstrates that the receiver-side plane-wave phase-encoded Hessian $\widetilde{H}({\bf x},{\bf y},{\bf p}_r)$ 
converges to the exact Hessian $H({\bf x},{\bf y})$ by integrating over receiver ray parameters.
Using the plane-wave phase-encoding functions, the approximate Hessian for a single ${\bf p}_r=(p_{r_x},p_{r_y})$ takes the form:
\begin{eqnarray}
& &\widetilde{H}({\bf x}, {\bf y}, {\bf p}_r) \nonumber \\
& & = \Re \left\{ \right. \sum_{\omega}\omega^4
                        \sum_{{\bf x}_s}|f_s(\omega)|^2 A_r^2(\omega) G({\bf x},{\bf x}_s,\omega) G^{*}({\bf y},{\bf x}_s,\omega) \nonumber \\
& & \times \sum_{{\bf x}_r} w({\bf x}_r,{\bf x}_s) G({\bf x},{\bf x}_r,\omega) {\rm e}^{i \omega {\bf p}_r \cdot {\bf x}_r} \nonumber \\
& & \times \sum_{{\bf x}_r^{\prime}} w({\bf x}_r^{\prime},{\bf x}_s)
                                  G^{*}({\bf y},{\bf x}_r^{\prime},\omega) {\rm e}^{-i \omega {\bf p}_r \cdot {\bf x}_r^{\prime}} \left. \right\}.
\end{eqnarray}
Integrating over ${\bf p}_r$ from $-\infty$ to $+\infty$, and changing the order of integration and summation, we have
\begin{eqnarray}
& &\widetilde{H}({\bf x},{\bf y}) = \int_{-\infty}^{+\infty}\widetilde{H}({\bf x}, {\bf y}, {\bf p}_r) d{\bf p}_r  \nonumber \\
& &= \Re \left\{ \right. \sum_{\omega}\omega^4\sum_{{\bf x}_s} |f_s(\omega)|^2 A_r^2(\omega) 
G({\bf x},{\bf x}_s,\omega) G^{*}({\bf y},{\bf x}_s,\omega) \nonumber \\
& & \times \sum_{{\bf x}_r} w({\bf x}_r,{\bf x}_s)G({\bf x},{\bf x}_r,\omega) 
\sum_{{\bf x}_r^{\prime}} w({\bf x}_r^{\prime},{\bf x}_s)G^{*}({\bf y},{\bf x}_r^{\prime},\omega) \nonumber \\
& & \times \int_{-\infty}^{+\infty}{\rm e}^{-i \omega {\bf p}_r \cdot ({\bf x}_r^{\prime}-{\bf x}_r)}d{\bf p}_r \left. \right\}. 
\end{eqnarray}
Note that
\begin{equation}
\int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_r \cdot ({\bf x}_r^{\prime}-{\bf x}_r)} d{\bf p}_r = \frac{1}{|\omega|}\delta({\bf x}_r^{\prime}-{\bf x}_r)
\end{equation}
in 2-D and
\begin{equation}
\int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_r \cdot ({\bf x}_r^{\prime}-{\bf x}_r)} d{\bf p}_r = \frac{1}{|\omega|^2}\delta({\bf x}_r^{\prime}-{\bf x}_r)
\end{equation}
in 3-D. Thus, 
if we choose the real function $A_r(\omega)$ such that it satisfies
$A_r^2(\omega)=|\omega|$ in 2-D and $A_r^2(\omega)=|\omega|^2$ in 3-D, 
then we get
\begin{eqnarray}
\widetilde{H}({\bf x},{\bf y})
= H({\bf x}, {\bf y}). \label{appendix_eq_rec_proof}
\end{eqnarray}

\append{Proof of convergence of the source- and receiver-side simultaneously plane-wave phase-encoded Hessian}
This appendix demonstrates that the simultaneously plane-wave phase-encoded Hessian $\widetilde{\widetilde{H}}({\bf x},{\bf y},{\bf p}_s,{\bf p}_r)$ 
converges to the exact Hessian $H({\bf x},{\bf y})$ by 
integrating over both source and receiver ray parameters. The simultaneously plane-wave phase-encoded Hessian takes the form
\begin{eqnarray}
& & \widetilde{\widetilde{H}}({\bf x},{\bf y},{\bf p}_s,{\bf p}_r) \nonumber \\
& & = \Re \left\{ \right. \sum_{\omega}\omega^4|f_s(\omega)|^2 A_s^2(\omega) A_r^2(\omega) \nonumber\\
& & \times \sum_{{\bf x}_s}w_s({\bf x}_s) G({\bf x},{\bf x}_s,\omega) {\rm e}^{ i \omega {\bf p}_s\cdot{\bf x}_s} \nonumber \\
& & \times \sum_{{\bf x}_s^{\prime}}w_s({\bf x}_s^{\prime})G^{*}({\bf y},{\bf x}_s^{\prime},\omega) {\rm e}^{-i \omega {\bf p}_s\cdot{\bf x}_s^{\prime}} \nonumber \\
& & \times \sum_{{\bf x}_r}w_r({\bf x}_r)G ({\bf x},{\bf x}_r,\omega) {\rm e}^{ i \omega {\bf p}_r\cdot{\bf x}_r} \nonumber \\
& & \times \sum_{{\bf x}_r^{\prime}}w_r({\bf x}_r^{\prime})G^{*}({\bf y},{\bf x}_r^{\prime},\omega) {\rm e}^{-i \omega {\bf p}_r\cdot{\bf x}_r^{\prime}} \left. \right\}.
\end{eqnarray}
Integrating over ${\bf p}_s=(p_{s_x},p_{s_y})$ and ${\bf p}_r=(p_{r_x},p_{r_y})$, and changing the order of summation and integration, we get
\begin{eqnarray}
& &\widetilde{\widetilde{H}}({\bf x},{\bf y})=
\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\widetilde{\widetilde{H}}({\bf x},{\bf y},{\bf p}_s,{\bf p}_r)d{\bf p}_s d{\bf p}_r \nonumber \\
& & = \Re \left\{ \right. \sum_{\omega}\omega^4|f_s(\omega)|^2 A_s^2(\omega) A_r^2(\omega) \nonumber\\
& & \times \sum_{{\bf x}_s} w_s({\bf x}_s)G ({\bf x},{\bf x}_s,\omega) 
    \sum_{{\bf x}_s^{\prime}}w_s({\bf x}_s^{\prime}) G^{*}({\bf y},{\bf x}_s^{\prime},\omega) \nonumber \\
& & \times \int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_s\cdot({\bf x}_s^{\prime}-{\bf x}_s)}d{\bf p}_s \nonumber \\
& & \times \sum_{{\bf x}_r} w_r({\bf x}_r) G ({\bf x},{\bf x}_r,\omega) 
    \sum_{{\bf x}_r^{\prime}} w_r({\bf x}_r^{\prime})G^{*}({\bf y},{\bf x}_r^{\prime},\omega) \nonumber \\
& & \times \int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_r\cdot({\bf x}_r^{\prime}-{\bf x}_r)}d{\bf p}_r\left.\right\}.
\end{eqnarray}
Once again note that
\begin{eqnarray}
\int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_s \cdot({\bf x}_s^{\prime}-{\bf x}_s)} d{\bf p}_s&=& \frac{1}{|\omega|}\delta({\bf x}_s^{\prime}-{\bf x}_s), \\
\int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_r \cdot({\bf x}_r^{\prime}-{\bf x}_r)} d{\bf p}_r&=& \frac{1}{|\omega|}\delta({\bf x}_r^{\prime}-{\bf x}_r),
\end{eqnarray}
in 2-D and
\begin{eqnarray}
\int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_s \cdot({\bf x}_s^{\prime}-{\bf x}_s)} d{\bf p}_s&=& \frac{1}{|\omega|^2}\delta({\bf x}_s^{\prime}-{\bf x}_s), \\
\int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_r \cdot({\bf x}_r^{\prime}-{\bf x}_r)} d{\bf p}_r&=& \frac{1}{|\omega|^2}\delta({\bf x}_r^{\prime}-{\bf x}_r),
\end{eqnarray}
in 3-D. Thus, if we choose real functions $A_s(\omega)$ and $A_r(\omega)$ satisfying $A_s^2(\omega)=|\omega|$, $A_r^2(\omega)=|\omega|$ in 2-D and
$A_s^2(\omega)=|\omega|^2$, $A_r^2(\omega)=|\omega|^2$ in 3-D, 
then we have
\begin{eqnarray}
\widetilde{\widetilde{H}}({\bf x},{\bf y})=H({\bf x},{\bf y}).
\end{eqnarray}


\bibliographystyle{seg}
\bibliography{refs}

