\title{Wave-equation inversion using 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}

\author{Yaxun Tang\footnotemark[1]}

\maketitle

\begin{abstract}
I present a new method for computing the Hessian operator for the wave-equation least-squares imaging problem, 
which is also known as the resolution function or the point-spread function, under the 
Born approximation. The proposed new method modifies the original explicit 
Hessian formula, enabling efficient computation of this operator. A particular 
advantage of this method is that it reduces or eliminates on-disk storage of 
Green's functions, a major obstacle for computing the Hessian
operator using the original formula.
The modifications, however, also introduce undesired cross-talk artifacts. 
I examine two different phase-encoding schemes, namely, 
plane-wave phase encoding and random phase encoding, to suppress the cross-talk. 
I apply the randomly phase-encoded Hessian operator to the Sigsbee2A synthetic data set, 
where an improved subsalt image with higher resolution 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 distorts
the amplitude of the reflectors because of the non-unitary nature of the Born modeling operator. 
To improve relative amplitude behavior and correct for the effects of limited acquisition geometry, complex overburden and band limited wavefields, 
the imaging problem can be treated 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,plessix2004,valenciano2006,yu2006}.
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.  

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 \cite{chavent1999} and \cite{plessix2004}
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 dominated 
\cite[]{pratt1998,chavent1999,plessix2004,valenciano2006}. It has been shown by \cite{valenciano2006} 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, the migrated image pre-multiplied by a diagonal matrix cannot perfectly recover the amplitude information, 
especially in poorly illuminated areas. That is why \cite{valenciano2006} and \cite{valenciano2008} suggest computing a limited number of 
the Hessian off-diagonals to compensate for poor illumination
and improve the inversion result.
However, computing even a limited number of the Hessian off-diagonals by directly implementing the explicit Hessian formula is very challenging. 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 explicit 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 cross-talk artifacts. 
I examine two phase-encoding schemes,
plane-wave phase encoding and random phase encoding, to attenuate the cross-talk. 
\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 the Sigsbee2A model.

\section{least-squares hessian}
\par
In general, primaries can be modeled by the Born modeling equation as follows:
\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 $d({\bf x}_r, {\bf x}_s, \omega)$ is the modeled data for a single frequency $\omega$ with the source and receiver located
at ${\bf x}_s=(x_s,y_s,0)$ and ${\bf x}_r=(x_r,y_r,0)$ on the surface; $f_s(\omega)$ is the source signature; 
and $m({\bf x})$ denotes the reflectivity 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}$.
In the examples throughout this paper, Green's functions are computed by means of one-way wavefield extrapolation. 
\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 range of a shot at } {\bf x}_s; \\
                         0 & \mbox{ otherwise }. \end{array} \right. \label{eq_mask}
\end{eqnarray}
For the marine acquisition geometry, $w({\bf x}_r,{\bf x}_s)$ is similar to a band-limited diagonal matrix;
for Ocean Bottom Seismometer (OBS) or land acquisition geometry, where all shots share the same receiver array, $w({\bf x}_r,{\bf x}_s)$ is a rectangular matrix.
Figure \ref{fig:acquisition-mask} illustrates the acquisition mask matrices for these two typical geometries in 2-D cases.

\plot{acquisition-mask}{width=6in}{Acquisition mask matrices for different geometries in 2-D cases. Grays denote ones while whites denote zeros.
The left panel shows the matrix for a typical marine acquisition geometry; the right panel shows the matrix for a typical OBS or land acquisition geometry.}

\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)(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}) = ||{\bf W}({\bf d} - {\bf d}_{\rm obs})||_2^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.
Equation \ref{eq_objfnc_compact} can be minimized iteratively using the conjugate-gradient based methods \cite[]{nemeth1999,kuhl2003,marie2005}.
The main obstacle to its wide use is its cost;
each iteration costs about two migrations, which makes it impractical for large-scale problems. Another drawback is its
relatively slow convergence rate, since
only the gradient information of the objective function $F$ is used in the optimization process.

\par
Instead of solving equation \ref{eq_objfnc_compact} iteratively in the data space, we can reformulate it and solve it in the model space.
Because ${\bf L}$ is a linear operator, $F$ is a quadratic function. Its minimum is obtained when ${\bf m}$ satisfies:
\begin{equation}
{\bf m} = {\bf H}^{-1}{\bf L}^{\dagger}{\bf W}^{\dagger}{\bf d}_{\rm obs}, \label{eq_minimum}
\end{equation}
where $^{\dagger}$ means adjoint and ${\bf H}={\bf L}^{\dagger}{\bf W}^{\dagger}{\bf W}{\bf L}$ is the Hessian operator.
By recognizing that ${\bf L}^{\dagger}{\bf W}^{\dagger}{\bf d}_{\rm obs}$ is the migrated image ${\bf m}_{\rm mig}$, 
we can rewrite equation \ref{eq_minimum} as follows 
\begin{equation}
{\bf Hm} = {\bf m}_{\rm mig}.
\end{equation}
Now we can define our new objective function in the model space as follows:
\begin{equation}
J({\bf m}) = || {\bf Hm} - {\bf m}_{\rm mig} ||_2. \label{eq_objfnc_modl}
\end{equation}
Equation \ref{eq_objfnc_modl} implies that once we obtain the explicit Hessian operator $\bf H$, 
we can use any direct or iterative solvers to minimize the new objective function $J$.

\par
The main advantage of minimizing the objective function $J$ in the model space is 
that we can invert the model parameters in a target-oriented fashion \cite[]{valenciano2006,valenciano2008}, 
which substantially reduces the size of the problem. 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 and minimizing the objective function $J$. 
Since the cost of minimizing the objective function $J$ is relatively small once the explicit Hessian is obtained, 
various regularization schemes can be easily tried to improve the inversion result without any additional cost. 
This is another important advantage of formulating inverse problem in the model space.

\par
We obtain the Hessian as follows by taking the second derivative of the objective function $F$ with respect to 
the model parameters \cite[]{plessix2004,valenciano2006}; 
note that $w^2({\bf x}_r,{\bf x}_s)=w^{\prime}({\bf x}_r,{\bf x}_s)=w({\bf x}_r,{\bf x}_s)$:
\begin{eqnarray}
H({\bf x},{\bf y})  &=& \Re \left( \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) \right. \times \nonumber \\
                    & & \left. \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)\right), \label{eq_hessian_exact}
\end{eqnarray}
where $\Re$ denotes taking the real part of a complex value and $'$ means complex conjugation; 
${\bf y}$ is a neighbor point around the image point ${\bf x}$ in the subsurface.

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$ in the shot-profile domain.
It is cumbersome and very expensive to compute the exact Hessian by directly implementing the above equation,
because it requires saving a large number of Green's functions.
In the subsequent sections, I introduce an alternative method based on phase encoding to 
compute the Hessian operator, which reduces the computational requirement. 
As I will demonstrate, by using this approach, we do not need to save any
Green's functions, and the cost for computing the Hessian is also significantly reduced.


\section{phase-encoded Hessian}

\subsection{Encoding of the receiver-side Green's functions}
\par
Suppose we have the following operator:
\begin{eqnarray}
\widetilde{H}({\bf x}, {\bf y}) &=& \Re \left( \sum_{\omega} |\omega|^4 \sum_{{\bf x}_s} 
                       |f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega)G^{\prime}({\bf y},{\bf x}_s,\omega) \times \right. \nonumber\\
                         & & \left( \sum_{{\bf x}_r}         w({\bf x}_r,{\bf x}_s)
                                                      G ({\bf x},{\bf x}_r         ,\omega)\alpha ({\bf x}_r         ,\omega) \right) \times \nonumber \\
                         & & \left. \left( \sum_{{\bf x}_r^{\prime}}w({\bf x}_r^{\prime},{\bf x}_s)
                                       G^{\prime}({\bf y},{\bf x}_r^{\prime},\omega)\alpha^{\prime}({\bf x}_r^{\prime},\omega) \right) \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; 
$\alpha({\bf x}_r,\omega)$ and $\alpha({\bf x}_r^{\prime},\omega)$ 
are some weighting functions to be specified later.
Though equation \ref{eq_hessian_approx} can be implemented in a way similar to the direct implementation of equation \ref{eq_hessian_exact},
it offers more flexibility and can be very efficiently implemented without 
explicitly saving the Green's functions.
\par
With the extra summation and because of the linearity of the wave equation, 
the term $\sum_{{\bf x}_r} w({\bf x}_r,{\bf x}_s)G({\bf x},{\bf x}_r,\omega)\alpha({\bf x}_r,\omega)$  
can now be seen as the extrapolated wavefield at the image point ${\bf x}$, with the following composite source as the source function on the surface: 
\begin{eqnarray}
f_{c}({\bf x},{\bf x}_s,{\omega}) = \sum_{{\bf x}_r}w({\bf x}_r,{\bf x}_s)\delta({\bf x}-{\bf x}_r)\alpha({\bf x}_r,\omega)
\end{eqnarray}
The same thing holds for the other summation term, 
$\sum_{{\bf x}_r^{\prime}}w({\bf x}_r^{\prime},{\bf x}_s)G^{\prime}({\bf x}_r^{\prime},{\bf y},\omega)\alpha^{\prime}({\bf x}_r^{\prime},\omega)$, 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,\omega)$ corresponding to 
the receiver composite source $f_c({\bf x},{\bf x}_s,\omega)$ as
\begin{eqnarray} 
R({\bf x},{\bf x}_s,\omega) = \sum_{{\bf x}_r}w({\bf x}_r,{\bf s}_s)G({\bf x}_r,{\bf x},\omega)\alpha({\bf x}_r,\omega), \label{eq_r}
\end{eqnarray}
and a source wavefield corresponding to a point source as the source function 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}) = \Re \left( \sum_{\omega} |\omega|^4 \sum_{{\bf x}_s} S({\bf x},{\bf x}_s,\omega) S^{\prime}({\bf y},{\bf x}_s,\omega)
                           R({\bf x},{\bf x}_s,\omega) R^{\prime}({\bf y},{\bf x}_s,\omega)\right), \label{eq_hessian_approx1}
\end{eqnarray}
which means $\widetilde{H}({\bf x},{\bf y})$ can be computed by cross-correlating the source and receiver wavefields with their shifted complex conjugates 
(${\bf y}$ is the neighborhood point around ${\bf x}$).
Therefore, the algorithm for computing the receiver-side encoded Hessian 
is similar to the wave-equation shot-profile migration algorithm, except 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 correlate 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 invoke the imaging condition defined by equation \ref{eq_hessian_approx1} to produce the Hessian.
In other words, we do not have to save the Green's functions at all. 
Computing $\widetilde{H}({\bf x},{\bf y})$ is also efficient because multiple Green's functions are computed at 
the same time during the wavefield extrapolation.
After some simple algebraic manipulation, equation \ref{eq_hessian_approx} can be rewritten as follows:
\begin{eqnarray}
& & \widetilde{H}({\bf x}, {\bf y}) = \Re \left( \sum_{\omega} |\omega|^4 
                             \sum_{{\bf x}_s} |f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega) G^{\prime}({\bf y},{\bf x}_s,\omega) \times \right. \nonumber \\
                         & & \sum_{{\bf x}_r}w({\bf x}_r,{\bf x}_s) G({\bf x},{\bf x}_r,\omega)G^{\prime}({\bf y},{\bf x}_r,\omega)|\alpha({\bf x}_r,\omega)|^2 +\nonumber \\ 
                         & & \sum_{\omega} |\omega|^4 \sum_{{\bf x}_s} |f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega) G^{\prime}({\bf y},{\bf x}_s,\omega) \times \nonumber \\
                         & & \left. \sum_{{\bf x}_r} \sum_{{\bf x}_r^{\prime}({\bf x}_r\not={\bf x}_r^{\prime})}w({\bf x}_r,{\bf x}_s)w({\bf x}_r^{\prime},{\bf x}_s)
                             G({\bf x},{\bf x}_r,\omega)G^{\prime}({\bf y},{\bf x}_r^{\prime},
                             \omega)\alpha({\bf x}_r,\omega)\alpha^{\prime}({\bf x}_r^{\prime},\omega)\right) .
 \label{eq_hessian_approx3}
\end{eqnarray}
If we let the weighting function $\alpha({\bf x}_r,\omega)$ satisfy $|\alpha({\bf x}_r,\omega)|=1$, the first term 
in equation \ref{eq_hessian_approx3} becomes the exact Hessian $H({\bf x}, {\bf y})$; however, the second term can not be removed and it becomes
undesired cross-talk from the cross-correlations among different receiver-side Green's functions: 
\begin{eqnarray}
& &\widetilde{H}({\bf x}, {\bf y}) = H({\bf x},{\bf y}) +\nonumber \\
                         & & \Re \left( \sum_{\omega} |\omega|^4 \sum_{{\bf x}_s} 
                              |f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega) G^{\prime}({\bf y},{\bf x}_s,\omega) \times \right. \nonumber \\
                         & & \left. \sum_{{\bf x}_r} \sum_{{\bf x}_r^{\prime}({\bf x}_r\not={\bf x}_r^{\prime})}w({\bf x}_r,{\bf x}_s)w({\bf x}_r^{\prime},{\bf x}_s)                              G({\bf x},{\bf x}_r,\omega)G^{\prime}({\bf y},{\bf x}_r^{\prime},
                             \omega)\alpha({\bf x}_r,\omega)\alpha^{\prime}({\bf x}_r^{\prime},\omega) \right) .
 \label{eq_hessian_approx4}
\end{eqnarray}

\subsection{Simultaneous encoding of the source-side 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 the shots share the same receiver array. 
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, the following operator can be constructed:
\begin{eqnarray}
& &\widetilde{\widetilde{H}}({\bf x},{\bf y}) 
         = \Re \left( \sum_{\omega}|\omega|^4 \right. \times \nonumber \\
& & \left(\sum_{{\bf x}_s}          w_s({\bf x}_s)         f_s(\omega) G({\bf x},{\bf x}_s,\omega)          \beta({\bf x}_s,\omega) \right)
  \left(\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},\omega) \right)^{\prime} \times \nonumber\\
& & \left.  \left(\sum_{{\bf x}_r}          w_r({\bf x}_r)         G({\bf x},{\bf x}_r         ,\omega) \alpha({\bf x}_r,\omega) \right)
            \left(\sum_{{\bf x}_r^{\prime}} w_r({\bf x}_r^{\prime})G({\bf y},{\bf x}_r^{\prime},\omega) \alpha({\bf x}_r^{\prime},\omega) \right)^{\prime} \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. 
Let us once again define the composite source wavefield $S({\bf x},\omega)$ and composite receiver wavefield $R({\bf x},\omega)$ as follows:
\begin{equation}
S({\bf x},\omega) = \sum_{{\bf x}_s}f_s(\omega)G({\bf x},{\bf x}_s,\omega)\beta({\bf x}_s,\omega), \label{eq_comp_sou_obc}
\end{equation}
and
\begin{equation}
R({\bf x},\omega) = \sum_{{\bf x}_r}G({\bf x},{\bf x}_r,\omega)\alpha({\bf x}_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{equation}
\widetilde{\widetilde{H}}({\bf x}, {\bf y}) = \Re \left( \sum_{\omega}|\omega|^4S({\bf x},\omega)S'({\bf y},\omega)
                             R({\bf x},\omega)R'({\bf y}, \omega) \right) . \label{eq_hessian_both2}
\end{equation}
Equation \ref{eq_hessian_both2} tells us that to compute the simultaneously encoded Hessian $\widetilde{\widetilde{H}}({\bf x},{\bf y})$ 
for OBS or land acquisition geometry, only two wavefield propagations are required: 
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}. 
Once again, if we can choose proper weighting functions $\alpha({\bf x}_r,\omega)$ and $\beta({\bf x}_s,\omega)$ such 
that $|\alpha({\bf x}_r,\omega)|=1$ and $|\beta({\bf x}_s,\omega)|=1$, 
equation \ref{eq_hessian_both2} becomes the sum of the exact Hessian and the cross-talk as follows: 
\begin {eqnarray}
& &\widetilde{\widetilde{H}}({\bf x}, {\bf y}) = H({\bf x},{\bf y}) + \nonumber \\
     & & \Re \left( \sum_{\omega}|\omega|^4 \sum_{{\bf x}_s}\sum_{{\bf x}_s^{\prime}({\bf x}_s\not={\bf x}_s^{\prime})}|f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega)
     G^{\prime}({\bf y},{\bf x}_s^{\prime},\omega)\beta({\bf x}_s,\omega)\beta^{\prime}({\bf x}_s^{\prime},\omega)  \times \right. \nonumber \\
                         & & \left. \sum_{{\bf x}_r}\sum_{{\bf x}_r^{\prime}({\bf x}_r\not={\bf x}_r^{\prime})}G({\bf x},{\bf x}_r,\omega)
                G^{\prime}({\bf y},{\bf x}_r^{\prime},\omega)\alpha({\bf x}_r,\omega)\alpha^{\prime}({\bf x}_r^{\prime},\omega) \right) . \label{eq_hessian_both4}
\end{eqnarray}

From equations \ref{eq_hessian_approx4} and \ref{eq_hessian_both4}, 
we face a situation similar to that encountered in phase-encoding migration \cite[]{romero2000}; i.e., our exact Hessian is contaminated by
cross-talk artifacts, so we seek to define weighting functions $\alpha({\bf x}_r,\omega)$ and $\beta({\bf x}_s,\omega)$
that attenuate the cross-talk as much as possible. In the next two sections, I examine two different phase-encoding schemes
to attenuate the cross-talk, 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,\omega) &=& \alpha({\bf x}_r,{\bf p}_r,\omega) = A_r(\omega){\rm e}^{i \omega {\bf p}_r {\bf x}_r} \\ 
\beta ({\bf x}_s,\omega) &=& \beta ({\bf x}_s,{\bf p}_s,\omega) = A_s(\omega){\rm e}^{i \omega {\bf p}_s {\bf x}_s},\label{eq_wght_rec}
\end{eqnarray}
which are the well known plane-wave phase-encoding functions \cite[]{whitmore1995,zhang2005,liu2006}, 
where $A_r(\omega)$ and $A_s(\omega)$ are real functions depending upon the angular frequency $\omega$,
$i = \sqrt{-1}$; 
and ${\bf p}_r=(p_{r_x},p_{r_y})$ and ${\bf p}_s=(p_{s_x},p_{s_y})$ are the surface ray-parameters for the receiver and source plane waves, respectively. 
As proved by \cite{liu2006}, summing an infinite number of densely sampled plane waves
completely attenuates the cross-talk 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: stacking 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,
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})$.

\subsection{Random phase encoding}
Instead of using the plane-wave encoding phase function, we can use random phases to disperse unwanted cross-terms \cite[]{romero2000}. 
The weighting functions can be chosen as follows:
\begin{eqnarray}
\alpha({\bf x}_r,\omega) &=& {\rm e}^{i \gamma({\bf x}_r,\omega)} \\
\beta ({\bf x}_s,\omega) &=& {\rm e}^{i \gamma({\bf x}_s,\omega)},
\end{eqnarray}
where the phase functions $\gamma({\bf x}_r,\omega)$ and $\gamma({\bf x}_s,\omega)$ are sequences of random
numbers between $0$ and $2\pi$.
When we sum over $\omega$ and aslo over 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 from each frequency, a uniformly distributed random sequence could be used \cite[]{romero2000}. To further attenuate the cross-talk, 
multiple realizations of the random sequences can be used: 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 cross-correlation among different Green's functions. Since the cross-correlation 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.
\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. 
Table \ref{tb_cost} illustrates the number of wavefield propagations and the size of Green's functions 
for an image cube with the size $N_x'\times N_y' \times N_z'$ for different methods. 
$N_{\omega}$ is the number of frequencies of the wavefield, $N_{p_{s_x}}$ and $N_{p_{s_y}}$ are the number of source ray parameters for the in-line 
and cross-line directions, and $N_{p_{r_x}}$ and $N_{p_{r_y}}$ are the number of receiver ray parameters for the in-line
and cross-line directions respectively. $N_{\rm realize}$ is the number of realizations for random phase encoding. 
For simultaneous encoding methods, OBS or land recording geometry has been assumed. 
\begin{table}[h]
\caption{Comparison of the costs for different methods}
\centering
\begin{tabular}{|c|c|c|}
\hline
Method & Wavefield propagations & Size of Green's functions \\
\hline
direct computation & $N_x N_y N_{\omega}$ & $N_xN_yN_x'N_y'N_z'N_{\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 & $2N_sN_{\omega}N_{\rm realize}$ & $0$ \\
random, simultaneous  & $2N_{\omega}N_{\rm realize}$ & $0$ \\
\hline
\end{tabular}
\label{tb_cost}
\end{table}
\par
Considering that $N_{p_{s_x}}, N_{p_{s_y}}, N_{p_{r_x}}$ and $N_{p_{r_y}}$ are usually much smaller than $N_x$ and $N_y$ in practice, 
the phase-encoding methods generally need far fewer wavefield propagations. The most efficient method would be 
the simultaneous random encoding method, which requires only $2N_{\omega}$ wavefield propagations for one realization. 
If the acquisition geometry is not OBS or land geometry, 
it would be efficient to use the receiver-side random phase-encoding method; its cost for one realization is comparable to that of a shot-profile migration. 
Besides needing fewer wavefield propagations, 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 ($N_xN_yN_x'N_y'N_z'N_{\omega}$) can easily 
reach a unaffordably large number.

\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 Hessian computation; 
the second is the more complicated Sigsbee2A velocity model. 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 show how 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=2000$m/s) for two different acquisition geometries.
Figure \ref{fig:const-hess-comparison} shows the diagonal of the Hessian (when ${\bf x} = {\bf y}$).
The diagonal of the Hessian is computed using different methods, with an acquisition geometry 
containing only one shot located at $-600$m and two receivers located at $600$m and $1200$m on the surface.
The frequency band of the seismic experiment is from $5$ Hz to $35$ Hz, and all the frequencies are used to generate the following results.
Figure \ref{fig:const-hess-comparison}(a) is the exact diagonal of the Hessian, $H({\bf x},{\bf x})$, uncontaminated by any cross-talk artifacts. 
Figure \ref{fig:const-hess-comparison}(b) is obtained using equation \ref{eq_hessian_approx} with $\alpha({\bf x}_r,\omega)=1$, 
so the result is contaminated by strong cross-talk, which can be easily identified as the vertical stripes on the right side of the image.
Figure \ref{fig:const-hess-comparison}(c) is the receiver-side plane-wave phase-encoded Hessian; $61$ receiver-side plane waves have been
stacked to generate this result. 
The cross-talk has been successfully removed, and Figure \ref{fig:const-hess-comparison}(c) looks 
very similar to Figure \ref{fig:const-hess-comparison}(a).
Figure \ref{fig:const-hess-comparison}(d) shows receiver-side randomly phase-encoded Hessian; only one realization has been computed. 
The cross-talk is greatly reduced by 
the random phase functions, but is not completely attenuated.
The result can be further improved by stacking more receivers or shots or using multiple realizations of the random phase functions.

Figure \ref{fig:const-hess-offd-comparison} shows the Hessian with off-diagonals (with size $21\times21$) at the image point $x=680$m, $z=800$m. 
The acquisition geometry is the same as that in Figure \ref{fig:const-hess-comparison}. 
The horizontal and vertical axes in Figure \ref{fig:const-hess-offd-comparison} show the horizontal and 
vertical offsets away from the image point $x=680$m, $z=800$m. Figure \ref{fig:const-hess-offd-comparison}(a) is the exact Hessian; 
Figure \ref{fig:const-hess-offd-comparison}(b) shows the Hessian with strong cross-talk; 
Figures \ref{fig:const-hess-offd-comparison}(c) and (d) show the results obtained with receiver-side plane-wave phase encoding and 
receiver-side random phase encoding methods, respectively.
Both phase-encoding schemes successfully remove the cross-talk shown in Figure \ref{fig:const-hess-offd-comparison}(b).

%\plot{const-hess-comparison}{width=6in,height=3.25in}
\plot{const-hess-comparison}{width=6in}
{Diagonal of the Hessian for a constant-velocity model with only one shot and
                                        two receivers. 
                                        (a) The exact diagonal of the Hessian (equation \ref{eq_hessian_exact});
                                        (b) the Hessian contaminated by cross-talk (equation \ref{eq_hessian_approx} with $\alpha({\bf x}_r,\omega)=1$); 
                                        (c) the receiver-side plane-wave phase-encoded Hessian;
                                        (d) the randomly phase-encoded Hessian.}


%\plot{const-hess-offd-comparison}{width=6in,height=1.5in}
\plot{const-hess-offd-comparison}{width=6in}
{The local Hessian operator for image point $x=680$m, $z=800$m.
 The acquisition geometry is the same as that in Figure \ref{fig:const-hess-comparison}. 
 The size of the Hessian operator is $21$ samples in both $x$ and $z$ directions.
 (a) The exact Hessian operator (equation \ref{eq_hessian_exact}); 
 (b) The Hessian contaminated by cross-talk (equation \ref{eq_hessian_approx} with $\alpha({\bf x}_r,\omega)=1$);
 (c) the receiver-side plane-wave phase-encoded Hessian;
 (d) the randomly phase-encoded Hessian.}

In the next example, I slightly change the acquisition geometry used in the previous example from 
two receivers to $401$ receivers. The receivers 
range from $-2000$m to $2000$m, with a spacing of $10$m. 
Once again, Figure \ref{fig:const-hess-more-comparison} shows the diagonal of the Hessian operator, 
while Figure \ref{fig:const-hess-more-offd-comparison} illustrates the off-diagonals at 
the image point $x=680$m, $z=800$m.
The results shown in Figure \ref{fig:const-hess-more-comparison}(c) and Figure \ref{fig:const-hess-more-offd-comparison}(c) 
are obtained by using receiver-side plane-wave phase encoding, which perfectly removes the cross-talk.
Figure \ref{fig:const-hess-more-comparison}(d) and Figure \ref{fig:const-hess-more-offd-comparison}(d) show the results of 
receiver-side random phase encoding. The cross-talk in the bottom left and right corners is also well removed. Also note that
some random noise has been introduced in the background (Figure \ref{fig:const-hess-more-comparison}(d)).

%\plot{const-hess-more-comparison}{width=6in,height=3.25in}
\plot{const-hess-more-comparison}{width=6in}
{Diagonal of the Hessian for a constant velocity model with only one shot but $401$ receivers. 
 (a) The exact diagonal of the Hessian (equation \ref{eq_hessian_exact});
 (b) the Hessian contaminated by cross-talk (equation \ref{eq_hessian_approx} with $\alpha({\bf x}_r,\omega)=1$);
 (c) the receiver-side plane-wave phase-encoded Hessian;
 (d) the randomly phase-encoded Hessian.}

%\plot{const-hess-more-offd-comparison}{width=6in,height=1.5in}
\plot{const-hess-more-offd-comparison}{width=6in}
{The local Hessian operator for image point at $x=680$m and $z=800$m. 
The acquisition geometry is the same as that in Figure \ref{fig:const-hess-more-comparison}.
The size of the Hessian operator is $21$ samples in both $x$ and $z$ directions.
(a) The exact Hessian operator (equation \ref{eq_hessian_exact});
(b) The Hessian contaminated by cross-talk (equation \ref{eq_hessian_approx} with $\alpha({\bf x}_r,\omega)=1$);
(c) the receiver-side plane-wave phase-encoded Hessian;
(d) the randomly phase-encoded Hessian.}


\subsection{Inversion with the phase-encoded Hessian: Sigsbee2A}
\par
To demonstrate an application of the explicit Hessian operator, I apply the model-space inversion approach to the Sigsbee2A model.
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. 
I apply two different strategies: the first is to compute only the diagonal of the Hessian operator and use it to normalize the migrated image; 
this is often known as the 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$
to find an optimum reconstruction of the reflectivity ${\bf m}$. 
A linear conjugate-gradient solver and a simple damping regularization term has been used for the model-space iterative inversion.
 
\par
Figure \ref{fig:sigsb2a-hess-diag-random} shows the diagonal of the receiver-side randomly 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 the following equation:
\begin{eqnarray}
H({\bf x},{\bf x})\approx H_{SI}({\bf x},{\bf x}) &=& \Re \left( \sum_{\omega}|\omega|^4
                              \sum_{{\bf x}_s}|f_s(\omega)|^2 G({\bf x},{\bf x}_s,\omega) G^{\prime}({\bf x},{\bf x}_s,\omega) \right) \nonumber \\
                       &=& \Re \left( \sum_{\omega}|\omega|^4\sum_{{\bf x}_s}|f_s(\omega)|^2 |G({\bf x},{\bf x}_s,\omega)|^2 \right). \label{eq_souint}
\end{eqnarray}
Comparing equation \ref{eq_souint} with equation \ref{eq_hessian_exact}, 
$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 arrays. In other words, 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 inaccurate, illumination below the salt.

%\plot[t]{sigsb2a-hess-diag-random}{width=6in,height=3.5in}
\plot{sigsb2a-hess-diag-random}{width=6in,height=4in}
{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 the entire frequency band ($5$ Hz - $35$ Hz).}

%\plot[t]{sigsb2a-hess-diag-souint}{width=6in,height=3.5in}
\plot{sigsb2a-hess-diag-souint}{width=6in,height=4in}
{Source-wavefield intensity for the Sigsbee2A model obtained by using equation \ref{eq_souint},
 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 shadow zones below the salt are easily identified; Figure \ref{fig:sigsb2a-imag-ffd-decon-diag-random}
shows the result of normalizing the migrated image
with the diagonal of the Hessian shown in Figure \ref{fig:sigsb2a-hess-diag-random}.
For comparison, Figure \ref{fig:sigsb2a-imag-ffd-decon-diag-souint} shows the result of normalizing the migrated image with the source-wavefield intensity shown
in Figure \ref{fig:sigsb2a-hess-diag-souint}.
Figure \ref{fig:sigsb2a-imag-ffd-decon-diag-random} shows more balanced amplitude across the section than Figure \ref{fig:sigsb2a-imag-ffd-decon-diag-souint}, 
especially in areas below the salt. This is because Figure \ref{fig:sigsb2a-hess-diag-random} takes the limited receiver arrays into consideration, and hence 
better predicts the illumination pattern in the subsurface.

%\plot[t]{sigsb2a-imag-ffd}{width=6in,height=3.5in}
\plot{sigsb2a-imag-ffd}{width=6in,height=4in}
{Conventional one-way wave-equation shot-profile migration result of the Sigsbee2A model. 
 Note the shadow zones under the salt.}

%\plot[t]{sigsb2a-imag-ffd-decon-diag-random}{width=6in,height=3.5in}
\plot{sigsb2a-imag-ffd-decon-diag-random}{width=6in,height=4in}
{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[t]{sigsb2a-imag-ffd-decon-diag-souint}{width=6in,height=3.5in}
\plot{sigsb2a-imag-ffd-decon-diag-souint}{width=6in,height=4in}
{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}.}
\par
Figure \ref{fig:sigsb2a-hess-win} shows the receiver-side randomly encoded Hessian operator for a particular region of interest under the salt:
(a) shows the diagonal part of the Hessian operator; 
(b) is obtained by convolving the Hessian operator (with a size $21\times21$) with a collection of point scatterers in the model 
space. It demonstrates the varying shapes and non-stationarities of the Hessian operators across the model space. 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. 
\par 
Figure \ref{fig:sigsb2a-invt-win} shows the comparison between migration and inversion: (a) is the migrated image, 
and (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 residuals as a function of iteration number; the residual converges after about $12$ iterations.
In the inversion result, the vertical resolution is greatly enhanced; 
the shadow zones that in the migrated image are now filled in with structures, and the sediments and the large dipping fault
extend much closer to the salt body. However, the inversion result also shows increased noise, which might be caused by either the null space 
or the random noise introduced in the randomly phase-encoded Hessian, or both. 
How the cross-talk affects the inversion is still not very clear and is beyond the scope of this paper, 
it remains an active research area for further investigation.
To suppress the noise, a more sophisticated regularization term other than damping can be introduced in the inversion process.

\plot{sigsb2a-hess-win}{width=6in,height=6.5in}
{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=6in,height=6.5in}
{Comparison between migration and inversion. Panel (a) shows the conventional 
                                   one-way wave-equation shot-profile migration result; 
                                   (b) shows the inversion result using the receiver-side randomly phase-encoded Hessian operator.}

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

\section{conclusions}
\par
I introduce a method based on phase encoding that allows efficient computation of the explicit Hessian operator.
The proposed algorithm closely resembles shot-profile migration, 
except that slightly different imaging and boundary conditions for the source and receiver wavefields are used. 
With these modifications, no Green's functions need to be stored
on disk; hence, the phase-encoded Hessian has the potential to be applied to large-scale problems with a much lower cost. However, this method
also generates undesired cross-talk. I examine plane-wave phase encoding and random phase encoding to attenuate the cross-talk. 
Numerical examples demonstrate that these two phase-encoding schemes work well in suppressing the cross-talk. 
The inversion examples of using randomly phase-encoded Hessian on the Sigsbee2A model show that least-squares inversion
can produce an image superior to that produced by migration in areas with low illumination and shadow zones. 
Thus it provides a good tool for imaging complex geologies.

\section{acknowledgment}
I thank Alejandro Valenciano and Biondo Biondi for enlightening discussions. I thank the sponsors of 
the Stanford Exploration Project for their financial support, and I acknowledge SMAART JV for providing the synthetic
data set.

\append{Proof of convergence of the receiver-side plane-wave encoded Hessian}
This appendix demonstrates that the receiver-side plane-wave-encoded Hessian $\widetilde{H}({\bf x},{\bf y})$ 
converges to the exact Hessian $H({\bf x},{\bf y})$ by integrating over receiver ray parameters.
Using the weighting functions defined in equation \ref{eq_wght_rec}, 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) &=& \Re \left( \sum_{\omega}|\omega|^4
                                  \sum_{{\bf x}_s}|f_s(\omega)|^2 A_r^2(\omega) G({\bf x},{\bf x}_s,\omega) G^{\prime}({\bf y},{\bf x}_s,\omega) \times \right. \nonumber \\
                  & & \left( \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} \right) \times \nonumber \\
                  & & \left. \left( \sum_{{\bf x}_r^{\prime}} w({\bf x}_r^{\prime},{\bf x}_s)
                                  G^{\prime}({\bf y},{\bf x}_r^{\prime},\omega) {\rm e}^{-i \omega {\bf p}_r \cdot {\bf x}_r^{\prime}} \right) \right).
\end{eqnarray}
Integrating over ${\bf p}_r$ from $-\infty$ to $+\infty$, and changing the order of integration and summation, we have
\begin{eqnarray}
& &\int_{-\infty}^{+\infty}\widetilde{H}({\bf x}, {\bf y}, {\bf p}_r) d{\bf p}_r 
= \Re \left( \sum_{\omega}|\omega|^4\sum_{{\bf x}_s} |f_s(\omega)|^2 A_r^2(\omega) 
G({\bf x},{\bf x}_s,\omega) G^{\prime}({\bf y},{\bf x}_s,\omega) \times \right. \nonumber \\
& & \left. \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^{\prime}({\bf y},{\bf x}_r^{\prime},\omega) 
    \int_{-\infty}^{+\infty}{\rm e}^{-i \omega {\bf p}_r ({\bf x}_r^{\prime}-{\bf x}_r)}d{\bf p}_r \right). 
\end{eqnarray}
Note that
\begin{equation}
\int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_r ({\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 ({\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}
\int_{-\infty}^{+\infty}\widetilde{H}({\bf x}, {\bf y}, {\bf p}_r) d{\bf p}_r
= 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 encoded Hessian}
This appendix demonstrates that the simultaneously plane-wave-encoded Hessian $\widetilde{\widetilde{H}}({\bf x},{\bf y})$ 
converges to the exact Hessian $H({\bf x},{\bf y})$ by 
stacking over both source and receiver ray parameters. The simultaneously plane-wave-encoded Hessian takes the form
\begin{eqnarray}
\widetilde{\widetilde{H}}({\bf x},{\bf y},{\bf p}_s,{\bf p}_r) &=& \Re \left(\sum_{\omega}|\omega|^4|f_s(\omega)|^2 A_s^2(\omega) A_r^2(\omega) \times \right. \nonumber\\
& & \sum_{{\bf x}_s}w_s({\bf x}_s) G({\bf x},{\bf x}_s,\omega) {\rm e}^{ i \omega {\bf p}_s{\bf x}_s}
    \sum_{{\bf x}_s^{\prime}}w_s({\bf x}_s^{\prime})G^{\prime}({\bf y},{\bf x}_s^{\prime},\omega) {\rm e}^{-i \omega {\bf p}_s{\bf x}_s^{\prime}}\times \nonumber \\
& & \left. \sum_{{\bf x}_r}w_r({\bf x}_r)G ({\bf x},{\bf x}_r,\omega) {\rm e}^{ i \omega {\bf p}_r{\bf x}_r}
    \sum_{{\bf x}_r^{\prime}}w_r({\bf x}_r^{\prime})G^{\prime}({\bf y},{\bf x}_r^{\prime},\omega) {\rm e}^{-i \omega {\bf p}_r{\bf x}_r^{\prime}} \right).\nonumber \\
\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}
& &\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 
= \Re \left( \sum_{\omega}|\omega|^4|f_s(\omega)|^2 A_s^2(\omega) A_r^2(\omega) \times \right. \nonumber\\
& & \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^{\prime}({\bf y},{\bf x}_s^{\prime},\omega) \int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_s({\bf x}_s^{\prime}-{\bf x}_s)}d{\bf p}_s\times \nonumber \\
& & \left. \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^{\prime}({\bf y},{\bf x}_r^{\prime},\omega) \int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_r({\bf x}_r^{\prime}-{\bf x}_r)}d{\bf p}_r\right). \nonumber
\end{eqnarray}
Once again note that
\begin{eqnarray}
\int_{-\infty}^{+\infty} {\rm e}^{-i \omega {\bf p}_s ({\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 ({\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 ({\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 ({\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}
\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 = H({\bf x},{\bf y}).
\end{eqnarray}


\bibliographystyle{seg}
\bibliography{refs}

