Transcript Document

Gabor Deconvolution
Gary Margrave and Michael Lamoureux
=
Outline
•
•
•
•
•
•
The Gabor idea
The continuous Gabor transform
The discrete Gabor transform
A nonstationary trace model
Gabor deconvolution
Examples
=
Gabor Papers in the Research Report
1. Gabor Deconvolution by Margrave and Lamoureux
2. Constant-Q wavelet estimation via a nonstationary
Gabor spectral model by Grossman, Margrave,
Lamoureux, and Aggarwalla
3. Gabor deconvolution applied to a Blackfoot dataset by
Iliescu and Margrave
4. A ProMAX implementation of
nonstationary deconvolution
by Henley and Margrave
=
The Gabor Idea
• Dennis Gabor proposed (1946) proposed the
expansion of a wave in terms of “Gaussian
wave packets”.
• The mathematics for the continuous Gabor
transform emerged quickly.
• The discrete Gabor transform
emerged in the 1980s due
to Bastiaans.
=
The Gabor Idea
A seismic signal
Multiply
A shifted Gaussian
A Gaussian slice or
wave packet.
The Gabor Idea
A seismic signal
A suite of
Gaussian slices
Remarkably, the suite of Gaussian slices can be designed such that
they sum to recreate the original signal with high fidelity.
The Gabor Idea
Window center time
Window center time
The Gabor transform
Fourier
transform
time
Suite of Gabor slices
frequency
Gabor spectrum or
Gabor transform
The Gabor Idea
Window center time
Window center time
The inverse Gabor transform done two ways
frequency
time
Inverse
Fourier
transform
Sum
Sum
Inverse
Fourier
transform
The Continuous Gabor Transform
Forward transform
Vg s  , f  

Gaussian
“analysis” window
 s t  g t    e

The Gabor transform
The signal
2 if t
dt
Fourier
exponential and
integration over t
The Gabor transform is a 2D time-frequency decomposition of
a 1D signal.
The Continuous Gabor Transform
Inverse transform
s t  
 
Gaussian
“synthesis” window
  Vg s  , f    t    e
 
The signal
The Gabor transform
2 if t
dfd
Fourier exponential
and integration over f
and 
The inverse Gabor transform is a 2D integration over the timefrequency plane.
A discrete Gabor transform
In the discrete case, Gabor transforms have proven
difficult to implement without loss. The issue was
finally “solved” in the 1980’s using “frame theory”.
We have implemented a simpler, approximate
transform that is based on a convenient summation
property of Gaussians.
A discrete Gabor transform
The Gaussian summation property.
Summation curve
Gaussian half-width .1 s
window increment .05
A discrete Gabor transform
The Gaussian summation property.
We show the the sum of a discretely spaced set on Gaussians
Summation curve
on an infinite interval is given by

  1half-width
 g  t  k 
Gaussian
.1 s
k
2cos  2 t /   e
window increment .05
 T /  
First-order error term
T is the Gaussian half-width
 is the Gaussian spacing
2

A discrete Gabor transform
The Gaussian summation property.
T

is the Gaussian half-width
is the Gaussian spacing
Summation curve
T / 
First-order error term
Gaussian half-width .1 s
.5 window -21
decibels.05
increment
1.0
-85 decibels
1.5
-150 decibels
2.0
-340 decibels
The inter-Gaussian spacing should well less than the
Gaussian half-width.
Gaussian summation test
Gaussian half-width .1 seconds.
T / 
2
1.5
1
0.5
0.25
0.1
Gabor transform test
Unnormalized transform
Difference
After forward and
inverse Gabor transform
Seismic signal
Gabor transform test
Normalized transform
Difference
After forward and
inverse Gabor transform
Seismic signal
Normalization is a simple process that reduces end effects.
Gabor transform example
Minimum phase
wavelet
Reflectivity
Nonstationary
synthetic Q = 25
Gabor transform example
Reflectivity
Q decay curve
Wavelet
No clipping
Gabor transform example
Q decay curve
Clipping
Stationary Convolution Model
s t  

 w t    r   d

Stationary Convolution Model
in the Fourier domain
sˆ  f   wˆ  f  rˆ  f 
Nonstationary Convolution Model
sˆ  f   wˆ  f 

 Q  , f  r   e

 Q  , f   e
 f  / Q
nonstationary
stationary
2 if 
d
Nonstationary Convolution Model
We have proven that
Vg s  , f   wˆ  f Q  , f Vg r  , f 
Gabor transform
of seismic signal
is approximately given by
Fourier transform
of wavelet
Q-filter
Gabor transform
of reflectivity
Proof
2 i  u
2e
 i u 



ˆw
VgVs s , , w
ˆ  


u
,

r


u
g

u





du  du 
  ,  r u  g u    e
g
1 

1   fT 2
1   fT 2   T /  2  2 it / 


2 it /  
F   f 
e


f

e

e
e

e




 
 2i t '
   f   2 i   f
 





I   sg tt e wˆ  f     u, fdt r
ee2 if u du  e2gifˆ tdf2 ift 
u

 Vgs  ,    s  t  g  t    e
dt

 

wˆ 



ˆ
Vg s V ,s , 
w



f

2 i u





ˆ

w


u
,

r
u
g
u


e
du 








 

g

2 if t



sˆ  f   w
 tf,t /fQr  t  e
dt
ˆf ft /Q
iH
 t, f   e

therefore
Vg s  ,   wˆ    , Vg r  ,  
Nonstationary Convolution Model
example
Gabor transform of seismic signal
Nonstationary Convolution Model
example
Gabor transform of reflectivity
Nonstationary Convolution Model
example
Q filter transfer function
Nonstationary Convolution Model
example
Q filter times Gabor transform of reflectivity
Nonstationary Convolution Model
example
wavelet surface
Nonstationary Convolution Model
example
wavelet times Q filter times Gabor transform of reflectivity
Nonstationary Convolution Model
example
Gabor transform of seismic signal
Gabor Deconvolution
Gabor transform of
the seismic trace
Given this
Gabor transform of
the reflectivity
Estimate this
Recall the nonstationary trace model:
Vg s  , f   wˆ  f Q  , f Vg r  , f 
Gabor Deconvolution
Vg r  , f  

Vg s  , f 
wˆ  f  Q  , f 
Gabor Deconvolution
So, given this:
Vg s  , f 
We must estimate this:
wˆ  f  Q  , f 
Gabor Deconvolution
Estimation of
wˆ  f  Q  , f 
“the propagating wavelet”
Can be done by:
• Smoothing
Vg s  , f 
• Time-variant Burg + smoothing
• Least squares
modeling of
Vg s  , f 
• Lotsa other ways …
Gabor Deconvolution
Phase calculation
The phase can be locally minimum, implemented
by a Hilbert transform over frequency, or locally
zero.
Gabor deconvolution example
nonstationaryseismic
seismicsignal
signal
nonstationary
reflectivity
reflectivity
Initial Gabor spectrum
Initial Gabor/Burg spectrum
Smoothed Gabor spectrum
Smoothed Gabor/Burg spectrum
Deconvolved Gabor spectrum
Fourier algorithm
Deconvolved Gabor spectrum
Burg algorithm
Deconvolved Gabor spectrum
Actual reflectivity
Gabor deconvolution result
Fourier algorithm
Gabor deconvolution result
reflectivity
Gabor deconvolution result
Burg algorithm
Gabor deconvolution result
reflectivity
Least squares Q and wavelet estimates
Constant-Q wavelet estimation via a nonstationary Gabor
spectral model by Grossman et al.
We fit a constant-Q wavelet model to the Gabor transform
of the seismic signal to obtain expressions for a Q
estimate and the wavelet amplitude spectrum.
Least squares Q and wavelet estimates
The nonstationary convolution model predicts
Vg s  , f   wˆ  f    , f  Vg r  , f 
This, together with the nearly random nature of
suggests that we impose the model
Vg r  , f 
Vg s  , f   wˆ  f    , f 
where the equality is taken to hold in the least-squares sense.
Least squares Q and wavelet estimates
f


The Q and wavelet estimates
emerge as integrations of the
Gabor transform of the signal over
a domain in the time-frequency
plane.
The integration domain is bounded by hyperbolae
that are level-lines of the constant Q operator:
e
 f  / Q
Least squares Q and wavelet estimates
This method is still in development but we are hopeful that is
will prove superior to simple smoothing as the method of
spectral separation in Gabor deconvolution.
ProMAX Implementation
A ProMAX implementation of nonstationary deconvolution
by Henley and Margrave
This module is in the current CREWES software release. The
FORTRAN computation module is self-contained and can be
easily ported to other systems.
Both Fourier and Burg methods are included with spectral
separation by smoothing. The Q estimator is not there yet.
A word of caution: unwise choice of parameters can mean very
long run times.
Blackfoot testing
Gabor deconvolution applied to a Blackfoot dataset
by Iliescu and Margrave
Gabor/Burg deconvolution
Wiener spiking deconvolution
Blackfoot testing
Gabor
Blackfoot testing
Wiener
Blackfoot testing
Typical trace spectra
0 dB
-10 dB
-20 dB
Amplitude
Amplitude
0 dB
-10 dB
-20 dB
-30 dB
-30 dB
-40 dB
-40 dB
-50 dB
0 Hz
20 Hz
40 Hz
60 Hz
80 Hz
After Wiener
100 Hz
120 Hz
0 Hz
20 Hz
40 Hz
60 Hz
80 Hz
After Gabor
100 Hz
120 Hz
Blackfoot testing
Gabor/Burg stack
Blackfoot testing
Wiener stack
Blackfoot testing
Gabor/Burg (pre and post stack) and final migration
Blackfoot testing
Wiener (pre stack), TVSW and final migration
Conclusions
• The discrete Gabor transform can be developed
by simple Gaussian slicing.
• The nonstationary extension of the convolution
model predicts that the Gabor transform of a
seismic signal decomposes into the product of
the Fourier transform of the wavelet, the symbol
of the Q filter, and the Gabor transform of the
reflectivity.
Conclusions
• Gabor deconvolution is implemented by a spectral
factorization applied to the Gabor transform of a
trace.
• Gabor deconvolution generalizes Wiener
deconvolution to the nonstationary setting.
• Within the Wiener design gate Gabor and Wiener
are similar; elsewhere, Gabor seems better.
• A new Gabor deconvolution code is released and
shows promising results.
Acknowledgements
All of the following provided support
CREWES: Consortium for Research in Elastic Wave
Exploration seismology
NSERC: Natural Sciences and Engineering Research
Council of Canada
MITACS: Mathematics of Information Technology and
Complex Systems
PIMS: Pacific Institute of the Mathematical Sciences