\chapter{Introduction}\label{chapter:Intro}
\def\w{\omega}
\def\xiv{\mbox{\boldmath$\xi$}}
\INT
Passive seismic imaging is an example of wavefield interferometric
imaging.  In this case, the goal is the production of subsurface
structural images by recording the ambient noisefield of the Earth
with surface arrays of seismometers or geophones.  The images produced
with this technique are directly analogous to those produced with
conventional reflection seismic data.  Interferometric imaging is not
new science.  This statement is true for two reasons.  
\par
First, the idea of imaging the subsurface by correlating the 
recordings from subsurface sources was first  
introduced by \longcite{GEO33-02-02640269}.  That work 
provided a one-dimensional proof that the auto-correlation of time
series collected by a receiver array can produce the
equivalent to a zero-offset time section.  Subsequently, 
\longcite{Zhang.sep.60.85.art}, used planewave decomposition 
to prove the result in 3D over a homogeneous medium.  \longcite{fink}
developed the Green's function of a heterogeneous medium
with acoustic waves via correlation and validated the observational 
theory with an ultrasonic experiment.  \longcite{kees04}, through 
one-way reciprocity, 
proved that cross-correlating traces of the transmission response of
an arbitrary medium synthesizes the complete reflection response,
in the form of shot-gathers, collected in a conventional active source
experiment.  In fact, in recent years, there has been a sudden
increase in interest around interferometry that has generated a wealth
of literature \cite{wapenaar:SI1}.  
\par
Second, interferometry is simply a generalization of the seismic 
problem formulated by earthquake and active-source 
exploration seismology.  On the of the continuum of purposefulness,
active experiments occupy an endpoint where humans 
initiate an experiment over a specific patch of Earth that is 
of interest.  In this case the source function is (hopefully) zero
phase, impulsive, and time zero is precisely known.  Toward the middle of the 
continuum, humans install receivers over an area of interest and 
capitalize on naturally occurring far-field (teleseismic) 
earthquake energy to probe the subsurface.  Here, the source energy is 
areal, and time zero is defined by windowing data beginning at the 
onset of the earthquake arrival.  At the opposite end of 
this spectrum is this work, which asks if the complicated,
low amplitude ambient noisefield in the Earth, about which one has no 
knowledge, can provide sufficient energy to image the subsurface
\cite{artman:SI177}.  However, the physics, and mathematics used to 
manipulate these various data sets are the same.  The more purposeful 
the experiment, many of the details of seismic interferometry can be 
simplified or are implicitly included in the form of the 
acquisition.  The less one knows about the source energy probing the 
subsurface, the more precise the processing needs to be in terms of 
the rigorous formulation of seismic interferometry.
\section{Theory}
\subsection{Natural phase encoding}
A simple explanation of how correlation of subsurface
sources returns the conventional active seismic data volume is in
terms of a phase encoding experiment.  Phase encoding shot records for
simultaneous processing was discussed in
\longcite{GEO65-02-04260436}.  The the motivation for phase encoding
is to manipulate the independent shot gathers in a predictable manner
such that the entire data volume can be processed (in this case
migrated) simultaneously.  To that end, each shot gather is convolved
with a unique function that will minimize cross-talk
interference between the shots during a summation over the shot axes
in a 3D seismic survey.  Random time series are quite effective
encoding functions.  Planewave migration strategies can also be viewed
as a phase encoding migration operation.  Consider a set of arbitrary encoding
functions $a({\bf x};{\bf x}_s,\w)$ for each shot ${\bf x}_s$ 
over the map coordinates ${\bf x}$ with
which shot records are convolved.  Normally, source functions for
active seismic data are bandlimited impulses, $I({\bf x}={\bf
  x}_s,\w)$, at the source location, and zero elsewhere.   After 
convolving the impulses with encoding functions the source 
wavefield, $\hat{S}$, required for simultaneous migration is
\begin{equation}\label{source}
\hat{S}_{z=0}({\bf x},\w)=\sum_{{\bf x}_s} 
a({\bf x};{\bf x}_s,\w)I_{z=0}({\bf x}={\bf x}_s,\w)\;.
\end{equation}
The five dimensional active seismic data, $R({\bf x},{\bf x}_s,\w)$, 
are similarly convolved with encoding functions and summed to reduce its
dimensionality
\begin{equation}\label{encode}
\hat{R}_{z=0}({\bf x},\w)=\sum_{{\bf x}_s} 
a({\bf x};{\bf x}_s,\w)R_{z=0}({\bf x};{\bf x}_s,\w)\;,
\end{equation}
where $R$ is the complete active-source data volume.  If the encoding
functions are perfectly uncorrelable, migrating all shot gathers
simultaneously,$\hat{R}$, provides a substantial cost savings.  Unfortunately,
perfect encoding functions do not exist, and a final image is produced
by summing the results of several encoded migrations with different
encoding functions.  With this technique, the cross-talk incurred 
by each migration will be out of phase and stacking the images will
reduce their effects.
\par
For seismic interferometry, we assume unknown and
arbitrarily complicated encoding functions arrive from the
subsurface.  Therefore a transmission record $T({\bf x},\w)$ recorded
by an array at the surface is $\hat{R}_{z=0}$.  If it is possible to
limit or window a passive recording around the arrival of a single
source, we can parameterize the transmission wavefield
\begin{equation}\label{encode2}
      T_{z=0} ({\bf x};\xiv,\w)=
\hat{R}_{z=0} ({\bf x};{\bf x}_s,\w)=
a({\bf x};\xiv,\w)R_{z=0}({\bf x};{\bf x}_s,\w)\;,
\end{equation}
where subsurface source $\xiv$ and active source ${\bf x}_s$ play
analogous roles in illuminating the subsurface.
\par
To recover $R_{z=0}({\bf x};{\bf x}_s,\w)$, we need 
to invert equation~\ref{encode} or~\ref{encode2}.  Since the adjoint
of convolution is correlation, correlating every trace with every
other trace would return a scaled version, depending on the amplitude
of $a$, of the active data $R$ except for the presence
of the sum in equation~\ref{encode}.  Unfortunately, the correlation
of a sum is not equal to the sum of correlations.  Therefore, the
summation over ${\bf x}_s$ will introduce cross-talk artifacts that
are convolutions of the encoding functions and the shot records one is
trying to recover.  The simple way to avoid this problem is to
parameterize the transmission wavefields as a function of individual
subsurface (encoding) functions, $\xiv$, and move the summation
outside of the correlation.  Thus, synthesizing conventional shot
gathers from transmission wavefields is predicated on
equation~\ref{encode2}, such that
\begin{equation}\label{uncode}
\tilde{R}_{z=0}({\bf x};{\bf x}_s,\w)=\sum_\xiv T({\bf x};\xiv,\w)T^*({\bf
  x}_s,\xiv,\w)=\sum_\xiv|a({\bf x},\xiv,\w)|^2R_{z=0}({\bf x};{\bf
  x}_s,\w)\;. 
\end{equation}
Above, $\tilde{R}$ includes a scalar amplitude term that
can be removed by normalizing by the amplitude spectrum, though is
kinematically equivalent to the active source data assuming that the
sum over $\xiv$ is dense.  If too few sources are captured,
$\tilde{R}$ may not accurately reproduce the kinematics of 
$R({\bf x}_s)$.
\par
Active seismic experiments produce data that is intrinsically 
real.  However, the above explanation includes causal and acausal 
parts.  If perfect illumination is achieved by capturing a full suite 
of subsurface sources ${\bf x}_s$, the two parts will be completely
symmetric and further processing with only the causal part of the result 
is appropriate.  However, at the limit of only a single subsurface 
source, the correlation is not symmetric as the kinematics
of the result is very close to the original wavefield, with a static
time shift applied.  Therefore, the degree to which the causal and
acausal parts of the result are asymmetric is indicative of 
the illumination problems associated with the distribution of 
subsurface sources captured in the experiment. 
\par
Similarly, the source wavefield, $\tilde{S}$, required to migrate $\tilde{R}$
will be sum of the correlations of all of the encoding functions $a$.  To 
numerically calculate $\tilde{S}$, it would be necessary to separate
the initial upcoming energy, $a({\bf x},\w)$, from the recorded data,
$T$, for correlation.  However, these
autocorrelations will produce the bandlimited impulse functions, 
$I_{z=0}({\bf x}={\bf x}_s,\w)$ in equation~\ref{source}.  Therefore, this 
procedure need not be performed, but instead a zero-phase wavelet
can be modeled that shares the bandwidth of the recorded data.  This 
means that
if sufficient source energy is captured and the resulting transmission
wavefields are correlated and stacked, the appropriate source functions
for migration are impulses at every recording station.  This result 
recapitulates that produced by \longcite{gerry} where it is shown 
that the Kirchoff migration kernel for migrating correlated gathers 
is the same as that used to image prestack active data when one assumes 
that impulsive virtual sources are located at all receiver locations.
\par
For convenience, the unknown source terms, $a$, can be thought
of as planewave constituents from all dips and azimuths or a set of point
sources completely surrounding the volume of interest.  If one is 
specifically interested in body waves probing subsurface structure, 
neither nearly vertical planewaves nor point sources near the surface are 
required to recover the information desired.  This recognition 
provides the possibility that synthesized gathers from subsurface
sources could be less prone to problematic wavefield kinematics
in particular applications.  For imaging subsurface structure with
body waves, synthesizing gathers with only small-valued ray-parameter 
source constituents would avoid direct arrivals and ground 
roll.  Alternatively, \longcite{louie}, described using surface waves
generated by cultural noise for near-surface refraction 
velocity analysis.
\par
Viewing passive seismic processing as an inverse phase-encoding 
problem reinforces the intuitive conception that a transmission 
wavefield is some convolution of the incoming energy from the 
subsurface and the resultant 
reflections from the subsurface after the energy has been mirrored 
at the free surface.  Correlation finds similar energy packets at
all times that correspond to the initial upcoming energy, and  
redefines time zero at each receiver location to correspond to the 
two-way traveltime to subsurface reflectors:  Initial 
waiting time for the arrival of energy at any location is removed
as the correlation calculates time differences between the initial
arrival and subsequent reflections.  Viewed in this light, passive
imaging is an identical operation, with a possibly more complicated
wavelet, to imaging primaries with multiples ADD REF. TO LATEST DELPHI WORK
\cite{guitton:1296}.
\par
Figure~\ref{fig:sketchP}, simplified for clarity,
shows the basic kinematics exploited in processing transmission
data.  The figure includes two recording stations 
capturing an approximately planar wavefront emerging from a two-layer
subsurface.  Panel ({\tt a}) shows the raypaths associated with the direct
arrival and one reflected both at the free surface and a subsurface
interface.  The second travel path (labeled reflection ray) has the 
familiar kinematics of the reflection seismic experiment if a source
were excited at receiver one.  The transmission
wavefield is depicted in panel ({\tt b}).  Wavelet polarity is
appropriate for direct arrivals and reflection.  The three main
features of transmission wavefields can be appreciated here.  First, the
exact timing of the arrival is unknown.  Second, the phase 
and duration of transmitted energy are unknown and likely
complicated.  Third, if the incident wavetrain coda is long,
arrivals in the transmission record can interfere.
\par
%%\sideplot{figure1}{width=3.3in, height=2.4in}{ ({\tt a}) Approximately
\sideplot{sketchP}{width=3.3in,height=2.4in}{({\tt a}) Approximately
  planar arrival with rays showing important propagation paths for
  passive imaging.  ({\tt b}) Idealized traces from a transmission
  wavefield. ({\tt c}) Shot-gather (reflection wavefield) synthesized
  using trace \textit{r1} as the source.  Many details, explained in
  the text, are neglected for simplicity.} 
Choosing trace \textit{r1} as the comparison trace, panel ({\tt c})
depicts the correlation spikes associated with the arrivals in the
data panel ({\tt b}), where $\otimes$ is correlation.  A solid line
with linear move-out is super-imposed across the correlated traces
corresponding to the direct arrival recorded at each receiver
location.  The dashed line on panel ({\tt c}) has hyperbolic moveout.
However, no correlation peak exists on the \textit{r1}$\otimes$
\textit{r1} trace under the hyperbola.  Not drawn, the second arrival
on \textit{r2} will have a counterpart on \textit{r1} from a ray
reaching the free surface farther to the left of the model.  
\par
In fact, the correlations produced from a single planewave will 
produce another planewave.  However, each planar reflection is moved 
to the lag-time associated with
a two-way trip from the surface to the reflector.  Correlation 
removes the wait time for the initial arrival and maintains the 
time differences between the direct arrival and reflections.  Summing the 
correlations from a full suite of planewaves builds hyperbolic events 
through constructive and destructive interference.  Analyzing seismic 
data in terms of planewave constituents is a commonly invoked tool in 
seismic processing.  Summing the correlations from incident planewaves 
is a planewave superposition process that builds hyperbolas.
\par
Passive seismic imaging is predicated on raypaths bouncing every
which-way from every direction.  Cartoons depicting the transmission 
experiment always leave something out that
causes an inconsistency that needs more raypaths and receivers to
explain.  Unfortunately the trend continues nearly forever.  Such a
complication arises with the inclusion of a second 
reflector.  The two reflection rays correlate with each other with a
positive coefficient.  The two travel paths share the time through the
shallow layer, so they correlate at a lag (or time difference), equal
to the two-way travel time through the deep layer.  This correlation
is not a problem however.  
\par
An upcoming interbed multiple event arrives after the direct arrival
delayed by the two-way travel-time through the deeper layer.  This event
has the opposite polarity compared to the direct arrival after once
changing its propagation direction from $\uparrow$ to $\downarrow$.  The 
delay of its arrival at receiver \textit{r2} compared to the direct
arrival at receiver \textit{r1} is the same lag as the correlation
between the reflectors and has opposite sign.  Therefore, the internal
multiple cancels potential artifacts of the correlation.  This shows
the importance of modeling transmission data with a two-way 
extrapolator.  Without all possible multiples, correlation artifacts
will quickly overwhelm the Earth structure within the correlated output. 
\par
Cross-correlation of each trace with every other trace handles the
three main difficulties of transmission data: timing, waveform, and
interference.  First, the output of the correlation is in lag units,
that when multiplied by the time sampling interval, provide the
time delays between like events on different traces.  The zero lag of the
correlation is the zero time for the synthesized 
shot-gathers.  Second, each trace records the character and duration
of the incident energy as it is reflected at the surface.  This
becomes the source wavelet analogous to a recorded vibrator 
sweep.  Third, overlapping wavelets are separated by correlation. 
\subsection{Generalization of interferometry}
To calculate the Fourier transform of the reflection response of a
loss-less subsurface, $R({\bf x}_r,{\bf x}_s, \omega)$,
\longcite{kees04}, proved
\begin{equation}\label{rigor}
2\Re [R({\bf x}_r,{\bf x}_s,\omega)]=
\delta({\bf x}_s-{\bf x}_r)-\int_{\partial\! D_m} 
T  ({\bf x}_r,\xiv,\omega)
T^*({\bf x}_s,\xiv,\omega)\;
\mbox{d}^2\xiv \;,
\end{equation}
where $^*$ represents conjugation and subscripts $r$ and
$s$ indicate different station locations from a transmission
wavefield.  After correlation $r$ and $s$ acquire the meaning of receiver
and source locations, respectively, associated with an active survey.  The 
RHS represents summing correlations of windows of passive data around
the  occurrence of individual sources from locations $\xiv$.  The
symbol $\partial \! D_m$ represents the domain boundary that surrounds
the subsurface region of interest on which the sources are located.  The 
transmission wavefields $T({\bf x}_r,\xiv,\omega)$ contain the arrival and
reverberations due to only one subsurface source.  To synthesize the
reflection experiment exactly, impulsive sources should completely
surround the volume of the subsurface one is trying to image.  Alternatively, 
many impulses can be substituted for a full suite of planewaves
emerging from all angles and azimuths as in the kinematic explanation
above.
\par
\longcite{wapenaar:2378} derives the correlation result above assuming a
loss-less medium, source-receiver reciprocity, and uses Gauss' theorem
which stipulates that the net flux through a domain boundary must be equal to
the net energy production within the volume of the domain.  The detailed 
mathematics of the derivation are inconsequential, but Figure~\ref{fig:domain} 
cartoons the starting point of the development and provides important 
understanding of the interferometric principle.  Of the many
derivations of seismic interferometry produced in the literature, the
lack of constraints on the geometry of the measurements makes this
result very general.   
\par
The first panel of Figure~\ref{fig:domain} shows a volumetric domain with
two sources surrounded by a surface containing receivers.  Three down-going
raypaths are drawn out of the nearly infinite set.  The two raypaths that 
hit the diffractor share the travelpath between the sources that is drawn 
as a solid line.  The second sketch switches the roles of sources and 
receivers via source-receiver reciprocity and the three raypaths drawn
travel inward.  The final sketch shows the special case when the
domain is truncated by a free surface.  The raypath recorded at the
left receiver from the top source is now recovered by a reflection
from the free surface due to the lowest source.  Similarly, the right-most 
source provides the raypath via reflection from the free surface that 
is lost by truncating the domain.  The total traveltime of these two 
new rays is different than in the second sketch.  However, the path
between the receivers is the same.  Therefore, correlation of the traces
recorded by the two receivers will return the same two-way traveltime 
to the diffractor in both cases.

\sideplot{domain}{width=3in}{Conservation of energy (Gauss' theorem)
foundation of the correlational proof \cite{wapenaar:2378}.  Net flux through
a domain boundary must equal the energy generated within.  Second figure
depicts source-receiver reciprocity.  Bottom figure shows the special
case when a free surface is included in the derivation.}
\par
For any receiver location in the cartoon, the physical
experiment shown would simply quantify the size of the energy sink
needed to remove all the acoustic energy within the domain due to the
sources surrounding it.  However, the energy conservation relations
can be rewritten in terms of Green's functions instead of particle
velocities.  The final result, equation~\ref{rigor}, is arrived at
by considering the difference between the recordings at any two 
receiver locations.  This is intuitively satisfying as energy that 
passes from one edge of the domain, through a single receiver, and out the 
opposite side does not contribute to correlations between two 
receivers.  In fact, a source from the opposite side of the boundary
will cancel such energy during summation over source location.  The 
stationary phase energy that is maintained is only the energy that 
passes through both receiver locations.  In the simple cartoon above, 
only two such raypaths exist.  However, in an arbitrarily complicated 
medium, scattering and reflection events are all maintained to synthesize
the Greens' function between the two receiver locations.
\par 
This generalized development allows several important conclusions.  First,
the receivers can be anywhere within the Earth volume.  Second, the 
correlation is intrinsically a two-station operator such that arrays 
need not be coplanar.  Last, the perfectly reflecting free-surface is not 
a requirement for the experiment, but rather a fortuitous mirror that 
allows complete synthesization of $R$ when sources above the Earth's 
surface are not physically practical.
\section{Overview of thesis}
While a wealth of literature is currently available describing many
different starting points to justify the correlation processing of
transmission wavefields to synthesize active seismic shot gathers,
very few real trials have been fielded.  This thesis presents several
experiments with synthetic data to explore potential problems and 
limitations of imaging the subsurface with passive data, and presents
results obtained from processing some 45 hours of passive data
recorded over the Valhall oil reservoir in the North Sea provided by
BP.  Because the method is intrinsically rooted in a poor
signal-to-noise regime, this data set provides a unique opportunity to
investigate the efficacy of extracting subsurface information from the
ambient noisefield given its size.  The array consists of 2300
four-component geophones distributed in 13 ocean-bottom cable (OBC)
lines over an area 6x12~km.  
\par
The basic correlation processing produces a 5D data volume that should
provide the kinematics of the active seismic experiment.  This product
is then roughly equivalent, and undoubtedly of lower quality, to data
that has been the intense focus of scientist and engineers for a very
long time.  Therefore, there is seemingly no end to the options and
strategies available for continued processing of this complicated data
volume.  
\par
I will focus on two main approaches in processing this field data.  First
is simple correlation, filtering, NMO, and stack.  The second approach
is based on the shot-profile migration method.  Because the
synthesized shot gathers likely suffer from low signal-to-noise
ratios, I apply migration as a method to combine the weak, redundant
signal distributed throughout the 5D volume of gathers to appropriate
subsurface reflection locations.  Migration not only collects the 
distributed subsurface information, but also drastically reduces the volume 
and complexity of data to analyze for subsurface information.
\par
\subsection{Chapter description}
Chapter~\ref{chapter:TimeSynth} will analyze the effectiveness of
correlation processing using a two-layer, 2D synthetic model.  Issues
such as source character, distribution, and the inclusion of multiple
source functions within a transmission wavefiel will be
illustrated with examples.  
\par
Chapter~\ref{chapter:TimeShot} introduces the Valhall field data and
describes and interprets the shot gathers generated by correlation and
spectral whitening.  Chapter~\ref{chapter:TimeCMP} extends the analysis
of the data by further processing the gathers with NMO and stacking.  
\par
Chapter~\ref{chapter:DepthSyn} Presents migrations of the correlated
gathers shown in Chapter~\ref{chapter:TimeSynth} to explore the
effects of the artifacts within synthesized shot gathers in the image
space.  I then introduce the concept of passive seismic processing as an
imaging problem within the course of shot profile migration.  The 2D
synthetic data is used to illustrate the procedure with the suite of
modeled data used previously.
\par
Chapter~\ref{chapter:DepthVal} presents migration results from the
Valhall passive data set produced by both correlation/imaging and
direct migration strategies.  An image provided by BP, produced with
data from an active seismic survey, is used as ground truth to
interpret the subsurface content of the final results produced with
the passive data. 
\par
Chapter~\ref{chapter:DepthMult} extends the concept of data processing
at the imaging condition of shot profile migration presented in
Chapter~\ref{chapter:DepthSyn} where the transmission wavefield was
used as both the upcoming and downgoing wavefields.  In this case,
shot gathers are used for both cases, and the resulting image is a 
prediction of the surface-related multiples in the depth domain.  The
prediction volume is equivalent to migrating the multiple prediction
produced by surface-related multiple prediction (SRMP) in the data
domain. 
%%\section{Conclusion}
%%There are many ways to arrive at the result that correlation of traces 
%%from a transmission wavefield can produce the kinematics of shot gathers
%%from an actively acquired seismic experiment.  Viewed in the framework
%%of Nature's own phase encoding, the result is immediately 
%%derived.  Beginning with first principles of energy conservation and 
%%source-receiver reciprocity provides comforting 
%%corroboration of the result as well as the understanding that 
%%acquisition geometry and the presence of the free surface are 
%%immaterial.  
%%\par
%%All derivations corroborate that the Earth volume can be arbitrarily
%%heterogeneous, but that the distribution of subsurface sources 
%%needs to be sufficiently distributed to recover the desired data.  If 
%%body waves investigating the subsurface structure are desired, 
%%low-angle planewave energy sources are mandatory (or a distribution
%%of subsurface point sources across the area of interest).  All 
%%derivations assume that the volume of investigation is loss-less which
%%may in practice prove a limiting constraint.  In effect, energy originating
%%at depth must make roughly three trips to an array on the surface of 
%%the Earth.  First energy must arrive at the array, then it needs to 
%%reflect down to make a two-way travelpath to discontinuous 
%%structure.  Given the absorption characteristics of the local Earth, 
%%and the energy budget of the ambient noisefield, this requirement
%%will constrain the usefulness of the technique to image subsurface
%%structure.
%%\bibliographystyle{sep}                                                         
%%\bibliography{SEP,GEOTLE,MISC,EAEG,SEG,ref}
