\documentclass[12pt]{article}
\usepackage{wsh}

\begin{document}  
\bibliographystyle{plain}

\markright{Wavelet Kirchhoff migration --- W.S. Harlan} 
\title{Kirchhoff migration with one-dimensional wavelets}

\author{William S. Harlan}

\date{November 1996}

\maketitle

\section {Introduction}

``Kirchhoff migration of seismic data compressed by matching pursuit 
decomposition,''  
by Bin Wang and Keh Pann \cite{SEG-1996-wavelet}, 
shows how to perform Kirchhoff summation migration directly in
the wavelet domain.
I am happy finally to see this simple idea in print.  The same
idea came to me years ago with a previous employer.  Now all
details are public knowledge, and I can talk about it again.
This paper also suggests many useful ideas that came from
practical implementation.  (For example, inverse wavelets can
include necessary rho filters.)

First I will provide my own favorite definition of a continuous
wavelet transform, which I wrote down in 1986 to rederive the results in 
Goupillaud et al
\cite{GoupillaudOne}.  (At the time in Zhuo Xian, China, 
I had an office down the hall from Pierre Goupillaud.)
This continuous transform seems much better suited 
to geophysical applications than the now more popular discrete
wavelet transform using Daubechies wavelets.  The latter have
a broad spectrum than necessary, more suitable for compressing
photographic images with sharp edges.   The continuous transform
can use simple Gaussian tapered monochromatic waves (one or multi-dimensional)
which maximizes the locality in both time/space and frequency/wavenumber.  

This new paper points out that the innermost loop of a Kirchhoff
migration can be described as a simple shift and sum, which is
a trivial operation in the wavelet domain.  (Although trivial,
the analytic details are given below, since I already have them
worked out.)  More precisely, input seismic traces are stretched, resampled, 
scaled, and then finally summed. If I perform
an equivalent operation directly on a wavelet-transformed trace,
then the cost could also be reduced manyfold. 

First, I will give the definition of what constitutes a
valid ``wavelet'' for the transform.   Second, the inverse
of a wavelet transform will be derived to demonstrate
why the transform works.  Last, I will show how the wavelet-transform
of a stretched trace can be written as a simple function
of the wavelet-transform of the original unstretched trace.

\section {Legal wavelets}

I will closely imitate the continuous transform originally
provided by Goupillaud, Grossman, and Morlet 
\cite{GoupillaudOne}.  They provide forward
and inverse transforms for related transforms, but
without derivation.

Use the following Fourier convention to link a 
wavelet $w(t)$ and its Fourier transform $\tilde w (s)$:
\begin{eqnarray}
\tilde w (s) &=& \int e^{-i 2 \pi s t }  w(t) dt {\rm ~and} \\
w(t) &=& \int e^{i 2 \pi s t }  \tilde w (s) ds . 
\end{eqnarray}
(All integrals with unlabeled extremes are assumed to be
over the entire real line.)

The following integral must exist for a valid wavelet.
\begin{equation}
C_w = \int { 1 \over |s| } \tilde w (s) ds .
\end{equation}
Thus, the spectrum of a valid wavelet must approach a zero
value at zero frequency: $\tilde w (s) / |s|^\epsilon $ $\rightarrow 0 $
as $s \rightarrow 0$ for some $\epsilon > 0 $.

This wavelet need not correspond to a seismic waveform.
Rather, when convolved with a trace, the wavelet should
suppress all but a band of frequencies.  Because of the
uncertainty principle, I must balance the narrowness
of the bandwidth with the narrowness in time.

I recommend a Gaussian-tapered sinusoid, such as
$w(t) = {\rm cos}(2\pi t) (e^{- \pi t^2 / c^2} )$. 
Such a wavelet is reasonably compact in both
the time and frequency domain. (The constant $c$
adjusts the relative compactness, within
the limits of the uncertainty principle.)

A valid wavelet can be used in an unfamiliar construction
of a delta function.
\begin{equation}
\label{eq:delta}
\int w(ut) du = C_w \cdot \delta (t) .
\end{equation}
In effect, this integral says that many wavelets, stretched
uniformly over all scales, will sum destructively everywhere but
at $t=0$.
\begin{equation}
{\bf Proof:~}  \int w(ut) du 
= \int \left \{ \int  {1 \over |u|} \tilde w ( { s \over u } ) 
e^{i 2 \pi s t } ds \right \} du \ldots
\end{equation}
Change variables from $(s,u)$ to $(s'=s,u'=s/u)$.
The Jacobian is \\
$\partial (s,u) / \partial (s',u') = |s'/u'^2|$.
Thus,
\begin{eqnarray} 
\int w(ut) du 
&=& \int \int  \left | {u' \over s'} \right | \tilde w ( u' ) e^{i 2 \pi s' t } 
\left | { s' \over u'^2 } \right | ds' du' \\
&=& \int \int  {1 \over |u'|} \tilde w ( u' ) e^{i 2 \pi s' t } ds' du' \\
&=& \left \{ \int {1 \over |u'|} \tilde w ( u' ) du'  \right \}  \cdot 
\left \{ \int e^{i 2 \pi s' t } ds' \right \} \\
&=& C_w \cdot \delta (t) . 
\end{eqnarray}
Here, I use the familiar Fourier definition of a delta function,
with the same restrictions.

The constraint of symmetry will also simplify the definition
of a wavelet transform in the following section, but
this assumption is optional:
\begin{equation}
w(t) = w(-t) .
\end{equation}

\section {Forward and inverse transforms}

Define a wavelet transform $F(u,a)$ of an $L^2$ function $f(t)$ by 
the following 
\begin{equation}
\label{eq:defn}
F(u,a) = \int w[u(a - t)] f(t) dt .
\end{equation}
This transform decomposes the function as a function of
position $a$ and local frequency $u$.  When calculated
discretely, the sampling need not be uniform over $a$ and $u$,
but the weights should reflect the above integration.

Goupillaud et all 
\cite{GoupillaudOne} prefer $2^u$ or $u^{-1}$ instead of $u$ 
as the stretch factor.
If you prefer an asymmetric wavelet, then you must use
two transforms---the above, and the integral with the time-reversed
wavelet.)

I find the following inverse, which simply scales and sums the transform
with uniform weight over local frequency $u$:
\begin{equation}
\label{eq:inverse}
f(t) = C_w^{-1} \int F(u, a=t) du .
\end{equation}
\begin{eqnarray}
{\bf Proof:~} 
C_w^{-1} \int F(u,a=t) du 
&=& C_w^{-1} \int \left \{ \int w[u(t - t' )] f(t') dt' \right \} du  \\
&=& C_w^{-1} \int \left \{ \int w[u(t - t' )] du \right \} f(t') dt' \\
&=& C_w^{-1} \int C_w \delta(t-t') f(t') dt' \\
&=& f(t). 
\end{eqnarray}
Here I make use of equation~\ref{eq:delta}.

\section {Summation imaging}

The previous sections were implicit in Goupillaud et al 
\cite{GoupillaudOne},
although I personally found it nontrivial to fill in the 
missing derivation of the results.  I was also happy 
to stumble on equation \ref{eq:delta}, which continues
to fascinate me.

This last section contains the only new idea particular
to Kirchhoff migration, and this idea is now public.
The idea is so obvious, really, that I'm sure everyone
who has considered applying Kirchhoff migration
to wavelets has thought of it.  I write down the
details just as a way of making the notation conform to
that of the preceding continuous wavelet transform.

The inner loop of a summation imaging procedure (Kirchhoff migration)
can be written as a simple stretching function:
\begin{equation}
g(\tau) = f[t(\tau)].
\end{equation}
The stretching function $t(\tau)$ gives the input data time $t$ 
as a monotonically increasing
function of the output image depth $\tau$.  
Moreover, I will assume that the stretching
function can be well approximated locally by a straight line.
\begin{eqnarray}
\label{eq:approx}
t(\tau) &\approx& t(\tau_0 ) + (\tau - \tau_0 ) 
\cdot {dt \over d\tau} ( \tau_0 ) ~~~~{\rm and }\\
\tau(t) &\approx& \tau_0 + [ t - t (\tau_0 ) ] 
\cdot \left [ {dt \over d\tau} ( \tau_0 ) \right ]^{-1} . 
\end{eqnarray}
Define forward and inverse wavelet transforms 
of $g(\tau)$ with the same wavelets as in definition~\ref{eq:defn}
and inverse~\ref{eq:inverse}.
\begin{equation}
G(v, b ) = \int w [ v ( b  - \tau ) ]  g (\tau) d\tau  ~~~~ {\rm and}  
\end{equation}
\begin{equation}
g(\tau) = C_w^{-1} \int G(v, b  = \tau ) dv.
\end{equation}
We would like to be able to express the transform $G(v, b)$
as a simple function of $F(u,a)$.  The approximation~\ref{eq:approx} 
allows us to write
\begin{equation}
\label{eq:stretch}
G(v,b) \approx 
\left | {dt \over d\tau} ( b ) \right |^{-1} 
 ~ F \left \{ u=v \cdot  
\left [ {dt \over d\tau} ( b ) \right ]^{-1} , ~ a = t(b) \right \} .
\end{equation}
We can perform an equivalent stretch on the transformed
input by stretching the position $a$, and by scaling
the local frequency $u$.  Both operations require a simple 
two-dimensional mapping of the transformed function.
The amplitude must be scaled as well.
\begin{eqnarray} 
{\bf Proof: ~} 
G ( v, b) &=& \int w[v(b-\tau)] \cdot f[t(\tau)] d \tau \\
&\approx& \int w [ v(b-\tau)] \cdot f \left [ t(b) + (\tau-b) 
\cdot {dt \over d\tau} ( b ) \right ]  d\tau \ldots 
\end{eqnarray}
Substitute the variable of integration:
\begin{eqnarray}
t' & \equiv& t(b) + (\tau-b) \cdot  {dt \over d\tau} ( b ) 
 ~~~~ {\rm and} \\
\tau &=& b + [ t' - t(b) ] 
\cdot \left [ {dt \over d\tau} ( b ) \right ]^{-1} .
\end{eqnarray}
Thus,
\begin{eqnarray}
G ( v, b) &=& \int w \left \{ v \cdot [t(b) - t'] \cdot
\left [ {dt \over d\tau} ( b ) \right ]^{-1} \right \} \cdot f(t') 
\cdot \left | {dt \over d\tau} ( b ) \right |^{-1} ~ dt'  \\
&=& \left | {dt \over d\tau} ( b ) \right |^{-1} 
\int w \left \{ v 
\left [ {dt \over d\tau} ( b ) \right ]^{-1} 
\cdot [t(b) - t'] 
\right \}  \cdot f(t') d t'  \\
&=& \left | {dt \over d\tau} ( b ) \right |^{-1} 
 ~ F \left \{ u=v \cdot  
 \left [ {dt \over d\tau} ( b ) \right ]^{-1} , ~ a = t(b) \right \} . 
\end{eqnarray}
The symmetry of $w(t)$ was not necessary here.
     
Equation~\ref{eq:stretch} is the result which allows us to perform
the stretch directly on the wavelet-transformed trace.

\bibliography{wsh}
\end{document}