Kirchhoff migration with one-dimensional wavelets

William S. Harlan

November 1996

Introduction

“Kirchhoff migration of seismic data compressed by matching pursuit decomposition,” by Bin Wang and Keh Pann [2], 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 [1]. (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 broader 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.

Legal wavelets

I will closely imitate the continuous transform originally provided by Goupillaud, Grossman, and Morlet [1]. 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)$:

$$\tilde w (s) = \int e^{-i 2 \pi s t } w(t) dt {\rm ~and}$$ (1)

$$w(t) = \int e^{i 2 \pi s t } \tilde w (s) ds .$$ (2)

(All integrals with unlabeled extremes are assumed to be over the entire real line.)

The following integral must exist for a valid wavelet.

  $$C_w = \int { 1 \over \vert s\vert } \tilde w (s) ds .$$ (3)

Thus, the spectrum of a valid wavelet must approach a zero value at zero frequency: $\tilde w (s) / \vert s\vert^\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.

  $$\int w(ut) du = C_w \cdot \delta (t) .$$ (4)

In effect, this integral says that many wavelets, stretched uniformly over all scales, will sum destructively everywhere but at $t=0$.

  $${\bf Proof:~} \int w(ut) du = \int \left \{ \int {1 \over \vert u\vert} \tilde w ( { s \over u } ) e^{i 2 \pi s t } ds \right \} du \ldots$$ (5)

Change variables from $(s,u)$ to $(s'=s,u'=s/u)$. The Jacobian is
$\partial (s,u) / \partial (s',u') = \vert s'/u'^2\vert$. Thus,

$$\int w(ut) du = \int \int \left \vert {u' \over s'} \right \vert \tilde w ( u' ) e^{i 2 \pi s' t } \left \vert { s' \over u'^2 } \right \vert ds' du'$$ (6)

$$= \int \int {1 \over \vert u'\vert} \tilde w ( u' ) e^{i 2 \pi s' t } ds' du'$$ (7)

$$= \left \{ \int {1 \over \vert u'\vert} \tilde w ( u' ) du' \right \} \cdot \left \{ \int e^{i 2 \pi s' t } ds' \right \}$$ (8)

$$= C_w \cdot \delta (t) .$$ (9)

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:

  $$w(t) = w(-t) .$$ (10)

Forward and inverse transforms

Define a wavelet transform $F(u,a)$ of an $L^2$ function $f(t)$ by the following

  $$F(u,a) = \int w[u(a - t)] f(t) dt .$$ (11)

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 [1] 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$:

  $$f(t) = C_w^{-1} \int F(u, a=t) du .$$ (12)

$${\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$$ (13)

$$= C_w^{-1} \int \left \{ \int w[u(t - t' )] du \right \} f(t') dt'$$ (14)

$$= C_w^{-1} \int C_w \delta(t-t') f(t') dt'$$ (15)

$$= f(t).$$ (16)

Here I make use of equation 4.

Summation imaging

The previous sections were implicit in Goupillaud et al [1], although I personally found it nontrivial to fill in the missing derivation of the results. I was also happy to stumble on equation 4, 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:

  $$g(\tau) = f[t(\tau)].$$ (17)

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.

$$t(\tau) \approx t(\tau_0 ) + (\tau - \tau_0 ) \cdot {dt \over d\tau} ( \tau_0 ) ~~~~{\rm and }$$ (18)

$$\tau(t) \approx \tau_0 + [ t - t (\tau_0 ) ] \cdot \left [ {dt \over d\tau} ( \tau_0 ) \right ]^{-1} .$$ (19)

Define forward and inverse wavelet transforms of $g(\tau)$ with the same wavelets as in definition 11 and inverse 12.

  $$G(v, b ) = \int w [ v ( b - \tau ) ] g (\tau) d\tau ~~~~ {\rm and}$$ (20)

  $$g(\tau) = C_w^{-1} \int G(v, b = \tau ) dv.$$ (21)

We would like to be able to express the transform $G(v, b)$ as a simple function of $F(u,a)$. The approximation 18 allows us to write

  $$G(v,b) \approx \left \vert {dt \over d\tau} ( b ) \right \vert^{... ...v \cdot \left [ {dt \over d\tau} ( b ) \right ]^{-1} , ~ a = t(b) \right \} .$$ (22)

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.

$${\bf Proof: ~} G ( v, b) = \int w[v(b-\tau)] \cdot f[t(\tau)] d \tau$$ (23)

$$\approx \int w [ v(b-\tau)] \cdot f \left [ t(b) + (\tau-b) \cdot {dt \over d\tau} ( b ) \right ] d\tau \ldots$$ (24)

Substitute the variable of integration:

$$t' \equiv t(b) + (\tau-b) \cdot {dt \over d\tau} ( b ) ~~~~ {\rm and}$$ (25)

$$\tau = b + [ t' - t(b) ] \cdot \left [ {dt \over d\tau} ( b ) \right ]^{-1} .$$ (26)

Thus,

$$G ( v, b) = \int w \left \{ v \cdot [t(b) - t'] \cdot \left [ {dt \over d\tau... ...\} \cdot f(t') \cdot \left \vert {dt \over d\tau} ( b ) \right \vert^{-1} ~ dt'$$ (27)

$$= \left \vert {dt \over d\tau} ( b ) \right \vert^{-1} \int w \left... ...t \over d\tau} ( b ) \right ]^{-1} \cdot [t(b) - t'] \right \} \cdot f(t') d t'$$ (28)

$$= \left \vert {dt \over d\tau} ( b ) \right \vert^{-1} ~ F \left \{ u=v \cdot \left [ {dt \over d\tau} ( b ) \right ]^{-1} , ~ a = t(b) \right \} .$$ (29)

The symmetry of $w(t)$ was not necessary here.

Equation 22 is the result which allows us to perform the stretch directly on the wavelet-transformed trace.

References

1

P. Goupillaud, A. Grossman, and J. Morlet.
Cycle–octave and related transforms in seismic signal analysis.
Geoexploration, 23:85–102, 1984.

2

Bin Wang and Keh Pann.
Kirchhoff migration of seismic data compressed by matching pursuit decomposition.
In 66th Annual International Meeting, SEG, Expanded Abstracts, pages 1642–1645. Soc. Expl. Geophys., 1996.