%% This file is automatically generated. Do not edit!
\documentclass[manuscript,revised]{geophysics}

\usepackage{color,amsmath}

\setfigdir{fig_matlab_new}

\begin{document}
\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.

I examine two different phase-encoding schemes to attenuate the crosstalk, 
namely, plane-wave-phase encoding and random-phase encoding.

\bibliographystyle{seg}
\bibliography{refs}

\end{document}
