\title{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. Revised manuscript submitted on 4/29/2009.
}

\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 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 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 a 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 an 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 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
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}. 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:
\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. 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.
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. 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. 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, 
or more specifically, encoding function, to be specified later. 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 
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.

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:
\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, 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'}w_s({\bf x}_s')f_s^*(\omega) 
    G^*({\bf y},{\bf x}_s^{\prime},\omega) \beta^*({\bf x}_s',{\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'}w_r({\bf x}_r')
    G^*({\bf y},{\bf x}_r^{\prime},\omega) \alpha^*({\bf x}_r',{\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}. 

\par
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}

\par
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. 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:
\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. 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
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.

\par
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 (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 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, 
\ref{fig:const-hess-offd-1s2r}e and \ref{fig:const-hess-offd-1s2r}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.}

\par
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.

\par
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 when 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. 

\par
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 
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 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 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. 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,sigsb2a-hess-offd-target}a 
shows the diagonal part of the Hessian operator; 
Figure \ref{fig:sigsb2a-hess-diag-target,sigsb2a-hess-offd-target}b 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,sigsb2a-invt-target,sigsb2a-refl-target}a is the migrated image, 
and Figure \ref{fig:sigsb2a-imag-target,sigsb2a-invt-target,sigsb2a-refl-target}b 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, 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-imag-target,sigsb2a-invt-target,sigsb2a-refl-target}c shows the filtered 
true reflectivity.

\par
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 line) better predicts the true reflectivity (the dotted line) 
than the migrated result (the 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}.
The dotted line, solid line and dashed line represent traces extracted from
the true reflectivity model, the migrated image and the inverted image, respectively. }

\section{discussion}
\par
The model-space inversion approach, which requires building the Hessian 
operator ${\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. 

\par
The inverted image, 
Figure \ref{fig:sigsb2a-imag-target,sigsb2a-invt-target,sigsb2a-refl-target}b, 
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. \new{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.}

\par
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. \new{The theory also assumes single scattering 
(Born approximation), other coherent energy, such as multiple, does not 
fit into this inversion 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,sigsb2a-invt-target,sigsb2a-refl-target}a) 
becomes even more pronounced in the inverted image 
(Figure \ref{fig:sigsb2a-imag-target,sigsb2a-invt-target,sigsb2a-refl-target}b).
Hence, multiple suppression is suggested prior to inversion.}

\section{conclusions}
\par
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. 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 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}
\par
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}
\par
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}
\par
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}