UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Deconvolution and wavelet estimation in exploration seismology Lines, Laurence Richard 1976

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1976_A1 L55.pdf [ 7.09MB ]
Metadata
JSON: 831-1.0052513.json
JSON-LD: 831-1.0052513-ld.json
RDF/XML (Pretty): 831-1.0052513-rdf.xml
RDF/JSON: 831-1.0052513-rdf.json
Turtle: 831-1.0052513-turtle.txt
N-Triples: 831-1.0052513-rdf-ntriples.txt
Original Record: 831-1.0052513-source.json
Full Text
831-1.0052513-fulltext.txt
Citation
831-1.0052513.ris

Full Text

DECONVOLUTION AND WAVELET E5 Tift A I I J N IM  EXPLORATION SEISMOLOGY  by  Laurence Richard L i n e s B.Sc,  University  o f A l b e r t a , 1371  M.Sc. , U n i v e r s i t y  o f AlL>erta, 197J  THESIS SUBMITTED IN PARTIAL FULFILMENT THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY  in  t h e Department of  Geophysics  We a c c e p t t h i s  and A s t r o n o m y  t h e s i s as conforming t o  reguired  The  University  standard  Of B r i t i s h  Columoia  A u g u s t , 1976  (c)  Laurence Richard Lines, 1976  In p r e s e n t i n g t h i s  thesis  an advanced degree at  further  fulfilment  of  the  requirements  the U n i v e r s i t y of B r i t i s h Columbia, I agree  the L i b r a r y s h a l l make it I  in p a r t i a l  freely  available  for  this  thesis  f o r s c h o l a r l y purposes may be granted by the Head of my Department  of  this thesis for  It  financial  or  i s understood that c o p y i n g o r p u b l i c a t i o n gain s h a l l not  written permission.  Depa rtment The U n i v e r s i t y o f B r i t i s h Columbia  2075 Wesbrook Place Vancouver, Canada V6T 1W5  that  r e f e r e n c e and study.  agree t h a t p e r m i s s i o n for e x t e n s i v e c o p y i n g o f  by h i s r e p r e s e n t a t i v e s .  for  be allowed without my  ABSTRACT Seismic d e c o n v o l u t i o n  i s investigated f o r signals  by e x p l o s i o n s and mechanical  vibrators.  the s e i s m i c t r a c e i s modelled impulse  response  In t h i s  investigation,  as the c o n v o l u t i o n of the e a r t h ' s  with a n o n s t a t i o n a r y wavelet. - By use of  s e v e r a l d e c o n v o l u t i o n techniques, the g e o l o g i c a l l y impulse  response  generated  interesting  i s estimated.  Wavelet e s t i m a t i o n i s o f t e n an e s s e n t i a l part of deconvolution. of  T h i s study g i v e s the f i r s t e x t e n s i v e  the Wiener-Levinson,  deconvolution  wavelet  comparison  Wold-Kolmogorov, and homomorphie e s t i m a t o r s i n e x p l o r a t i o n seismology.  In  t h i s d i s c u s s i o n , t h e Wold-Kolmogorov method i s shown to be e q u i v a l e n t t o H i l b e r t transform wavelet  estimation.  wavelet  comparisons on n o i s e l e s s s y n t h e t i c s i n d i c a t e t h a t homomorphic deconvolution  shows c o n s i d e r a b l e promise as a wavelet  However, homomorphic f i l t e r i n g a d d i t i v e noise i s present.  encounters  difficulty  estimator. when  Hence, c e p s t r a l s t a c k i n g i s used to  reduce the noise problem. T h i s t h e s i s a l s o i n v e s t i g a t e s multichannel Wiener f i l t e r s which e x h i b i t i d e a l performance when the wavelets are known. Despite such i d e a l performance on s y n t h e t i c data, the advantage of  m u l t i c h a n n e l Wiener deconvolution over s i n g l e channel  becomes marginal or nonexistant when wavelet A case h i s t o r y of - deconvolving r e f l e c t i o n data i s shown.  methods  estimates are used.  explosion-generated  Comparisons of deconvolutions  with  w e l l l o g s y n t h e t i c s are used i n order t o provide-an i n t e r p r e t a t i o n of the e a r t h s impulse 1  response  i n the r e g i o n of  ii interest.  Observed  and decon v o l u t i o n s  differences  are i n t e r p r e t e d  c h a r a c t e r i s t i c s of s e i s m i c A new created  between s y n t h e t i c i n terras of txie  s i g n a l s and v e l o c i t y  spectrum i s removed. response's F o u r i e r prediction.  seismograms  Using c e p s t r a l f i l t e r i n g  s p e c t r a l d i v i s i o n , the wavelet p o r t i o n  increase  logs.  approach i s presented f o r deconvolving  by v i b r a t o r sources.  seismoqrams  of the t r a c e ' s  The undetermined p o r t i o n s  transform are f i l l e d  i n by  Frequency domain p r e d i c t i o n can  the time domain r e s o l u t i o n of s e i s m i c  and amplitude  of the impulse autoregressive  sUDstautially arrivals.  This  d e c o n v o l u t i o n method i s p a r t i c u l a r l y well s u i t e a t o v i b r o s e i s data because of the phase c h a r a c t e r i s t i c s of the wavelet and the known b a n d l i m i t e d spectrum of the signal.  vibroseis vibroseis  Such an approach o b v i a t e s r e s t r i c t i v e assumptions  are i n h e r e n t  which  i n c o n v e n t i o n a l approaches to v i b r o s e i s  deconvolution.  The  u s e f u l n e s s of t h i s new  i s demonstrated f o r s y n t h e t i c  d e c o n v o l u t i o n method  and r e a l d a t a .  iii TABLE OF CONTENTS Page ABSTRACT  i  TABLE OF CONTENTS  i i i  LIST OF FIGURES  iv  ACKNOWLEDGEMENTS  vii  PREFACE  viii  CHAPTER 1. INTRODUCTION TO DECONVOLUTION 1.1 The Impulse Response of a Layered E l a s t i c Earth 1.2 Wavelet Model of a Seismic Trace 1.3 A t t a n u a t i o n and D i s p e r s i o n of a Seismic Pulse 1.4 Approaches to Deconvolution  1 9 13 14  CHAPTER 2.1 2.2 2.3 2.4  2. SEISMIC WAVELET ESTIMATION The Wiener-Levinson Double Inverse Method The Wold-Kolmogorov F a c t o r i z a t i o n Method Homomorphic Deconvolution Comparisons of Wavelet E s t i m a t o r s  19 23 27 32  CHAPTER 3.1 3.2 3.3  3. THE DESIGN OF INVERSE FILTERS FOR S i n g l e Channel Wiener Deconvolution M u l t i c h a n n e l Wiener Deconvolution Time Adaptive Deconvolution  52 61 82  WAVELETS  CHAPTER 4.. A CASE HISTORY OF DECONVOLVING REAL SEISMIC DATA 4.1 Deconvolutions Performed 4.2 Impulse Response Models Obtained from S y n t h e t i c Seismograms 4.3 R e l a t i n g Deconvolutions and S y n t h e t i c Seismograms CHAPTER .5. 1 5.2 5.3 5.4 5,5'  5. A NEW APPROACH TO VIBROSEIS DECONVOLUTION I n t r o d u c t i o n ' to the V i b r o s e i s Method A V i b r o s e i s Trace Model The V i b r o s e i s Deconvolution Problem The Deconvolution Method A p p l i c a t i o n s to S y n t h e t i c and A c t u a l Data  CHAPTER 6. REFERENCES  SUMMARY  103 111 119 130 131 136 137 151 170 173  iv L I S T OF FIGURES  Figure  Page  1  P wave propogating a t v e r t i c a l i n c i d e n c e through a G o u p i l l a u d earth model  2a)  P o s s i b l e displacement waveform and recorded wavelet f o r a s e i s m i c d i s t u r b a n c e .  12  2b)  V a r i a t i o n i n the shape of the s e i s m i c wavelet as a f u n c t i o n of propogation d i s t a n c e .  12  3  Flow c h a r t o u t l i n i n g approaches to s e i s m i c d e c o n v o l u t i o n and wavelet e s t i m a t i o n .  18  4  A comparison o f source wavelet estimates f o r a synthetic trace.  35  5  An example demonstrating noise on phase curves.  38  6  A s e t of f i v e t r a c e s from a given shot used i n comparison of wavelet e s t i m a t o r s .  40  7  Use of c e p s t r a l averaging on r e a l data o f Figure 3.  43  8  Use of l o g s p e c t r a l and raw phase averaging t o determine wavelet.  47  9  Comparison of wavelet obtained by l o g s p e c t r a l averaging with minimum phase wavelets.  48  10  Comparisons of deconvolutions of r e a l using d i f f e r e n t wavelet e s t i m a t e s .  50  11  The e f f e c t o f prewhitening i n the design of s i n g l e channel channel wiener deconvolution filters.  60  12  A flow diagram o f the Wiggins-Bobinson s o l u t i o n f o r the multichannel Wiener f i l t e r .  66  13a)  A s y n t h e t i c example e x h i b i t i n g the i d e a l performance o f the multichannel wiener deconvolution f i l t e r .  69  13b)  Sum of s i n g l e channel Wiener deconvolutions using the a c t u a l wavelets.  70  14a)  Multichannel deconvolution using estimates.  74  e f f e c t s of a d d i t i v e  data  3  wavelet  V I  14b)  Stacks of s i n g l e channel deconvolutions designed using wavelet estimates  76  Real data used i n a comparison between s i n g l e channel and multichannel d e c o n v o l u t i o n .  79  16  Comparison o f s i n g l e and multichannel d e c o n v o l u t i o n s on r e a l data.  81  17  Comparison of a u t o c o r r e l a t i o n s using the MEM and zero e x t e n s i o n methods.  91  18a)  Noisy s y n t h e t i c t r a c e used i n comparison of Wiener and HEM d e c o n v o l u t i o n .  93  18  b),c) A comparison of .Wiener and MEM prediction error f i l t e r i n g .  95  19  An a p p l i c a t i o n o f R i l e y - B u r g time deconvolution.  99  20  A bandpassed time a d a p t i v e MEM deconvolution.  102  21a)  Bandpass f i l t e r e d Gulf Coast.  106  21b)  Data d i s p l a y i n g e f f e c t of s p h e r i c a l correction.  divergence  106  22a)  Wiener minimum phase d e c o n v o l u t i o n of the data i n F i g u r e 21a) using no prewhitening.  108  22b)  Wiener minimum phase d e c o n v o l u t i o n o f the same data with 10% prewhitening.  108  23a)  CDP s t a c k s o f data i n F i g u r e 21 b ) .  HO  23b)  Wiener deconvolution of CDP s t a c k s with no prewhitening.  110  23c)  Wiener d e c o n v o l u t i o n s of CDP s t a c k s with prewhitening.  110  24  Check shot and CVL v e l o c i t y models with s y n t h e t i c seismograms d e r i v e d from the v e l o c i t y curves.  25  Comparison of the a u t o c o r r e l a t i o n s o f the CVL s y n t h e t i c before and a f t e r bandpass f i l t e r i n g .  118  26  I t e r a t i v e d e c o n v o l u t i o n o f common o f f s e t s t a c k s f o r near v e r t i c a l i n c i d e n c e .  121  15  adaptive  data from an area near  the  '\Q%  114  vi 27  Comparison o f a u t o c o r r e l a t i o n a t v a r i o u s stages i n the i t e r a t i v e deconvolution process.  123  28  D e c o n v o l u t i o n r e s u l t s o b t a i n e d by Amoco R e s e a r c h f o r t h e d a t a s e t examined h e r e .  125  29  Summary o f b e s t r e f l e c t i o n p i c k s from t h e d e c o n v o l u t i o n s and s y n t h e t i c seismograms o b t a i n e d from t h e G u l f C o a s t d a t a .  30  Example o f an i d e a l i z e d trace.  31 a) Example o f a K l a u d e r spectrum. 31  model f o r a v i b r o s e i s  wavelet, and i t s a m p l i t u d e  b ) c ) D e c o n v o l u t i o n s and s p e c t r a r e s u l t i n g from a u t o r e g r e s s i v e p r e d i c t i o n o f t h e F o u r i e r - . transform o f the Klauder wavelet. #  32 a) I n p u t  s p i k e sequence  and i t s s p e c t r u m .  32 b) M o d i f i e d t r a c e p r o d u c e d spectrum.  by b a n d l i m i t i n g t h e  32 c) D e c o n v o l u t i o n which r e s u l t s the F o u r i e r t r a n s f o r m .  from  predicting  33 a) The a m p l i t u d e s p e c t r u m o f a v i b r o s e i s w a v e l e t shown w i t h t h e t r a c e model and t h e t r a c e spectrum. 33 b) C o m p a r i s o n between t h e t r a c e spectrum and a c e p s t r a l e s t i m a t e of t h e wavelet spectrum. ' 33 c) D e c o n v o l u t i o n o f t r a c e i n F i g u r e 33 a) and r e s u l t i n g spectrum. 33 d) D e c o n v o l u t i o n o f a t r a c e s y n t h e s i z e d by u s i n g a minimum phase " e a r t h f i l t e r " . 34  Minimum phase d e c o n v o l u t i o n t r a c e model.  35  a) ,b) ,c) , I n p u t v i b r o s e i s d a t a shown subsequent stages of d e c o n v o l u t i o n .  1 2  8  1 3  4  143 145  147 148  149  1 5  4  1  5  5  1  5  7  158  of the vibroseis w i t h two  164  35 d) C o m p a r i s o n o f a m p l i t u d e s p e c t r a a t v a r i o u s stages of deconvolution.  166  36  169  A processed v i b r o s e i s s e c t i o n deconvolutions.  with  subsequent  I wish to express my Professor and  Tad  s i n c e r e thanks to my  U l r y c h , whose constant  encouragement, guidance,  generous a s s i s t a n c e mads t h i s p r o j e c t I thank Rob  applying  Clayton  Geophysical  possible.  for his invaluable assistance in  homomorphic deconvolution  prediction.  supervisor,  I a l s o thank Dr.  and  autoregressive  Ralph Wiggins of Western  f o r p r o v i d i n g s e v e r a l ideas on time sequence  analysis. The  a s s i s t a n c e of other  acknowledged.  I am  colleagues  g r a t e f u l to Dr.  s y n t h e t i c seismogram program and Gumming, John B.  Davies, and  Mat  Ron  at UBC  i s also  Clowes f o r p r o v i d i n g  to Dr. Paul  Soiaerville,  a  Bill  Y e d l i n f o r t h e i r programming  suggestions. I am indebted d i s c u s s i o n s on  to Dr.  wavelet e s t i m a t i o n  a l s o l i k e to thank Dr. Ken  Peacock, and  s e i s m i c and  Bill  and  deconvolution.  I would  L a r r y Wood, Dr. Paul Gutowski, Bob Thompson of Amoco f o r supplying  well l o g data and  V i b r o s e i s data and supplied  Sven T r e i t e l of Amoco f o r many h e l p f u l  data processing  me  with  suggestions.  s p e c i a l p l o t t i n g f a c i l i t i e s sere  by George Brinkworth, B i l l  Lucas,  kindly  D a v i t t and Dr. Stan Jones of  Chevron Standard. F i n a l l y , I would l i k e to thank the  National  Research  C o u n c i l of Canada f o r t h e i r generous f i n a n c i a l support.  T h i s t h e s i s i s e s s e n t i a l l y a d e s c r i p t i o n of approaches to deconvolution. extensively  discussed  various  Some of these approaches have been  i n the l i t e r a t u r e ,  a p p l i c a t i o n s of c e p s t r a l a n a l y s i s and  but  many are  autoregressive  prediction  to e x p l o r a t i o n seismology t h a t are presented for the f i r s t In the l a t t e r case, the research had my  a large influence.  To  of my  s u p e r v i s o r , Tad  In Chapter 1, I provide  A l s o , a general  techniques i s given  d i s c u s s i o n s i n Chapters 2 and  and  show t h a t Robinson's formulation  as a prelude to  of Wold-Kolmogorov  s p e c t r a l and and  Wiggins  raw  In e s t i m a t i n g data,  wavelets from  phase averaging which was  (1976) f o r the e s t i m a t i o n  inverse  algorithm  log  noisy  introduced  of t e l e s e i s m i c  by  Clayton  wavelets.  multichannel  Wiener  showed t h a t c e r t a i n p a t h o l o g i c a l cases of  f i l t e r i n g r e q u i r e prewhitening.  c r i t e r i o n f o r a multichannel  filter  examples of i d e a l multichannel time adaptive  of the  I a l s o a p p l i e d a technique of l o g  In Chapter 3, a c l o s e i n s p e c t i o n of the filter  the  In that d i s c u s s i o n , I  f a c t o r i z a t i o n i n v o l v e s t a k i n g a H i l b e r t transform  exploration seismic  o u t l i n e of  minimum phase wavelet  are compared i n Chapter 2.  amplitude spectrum.  preface.  3.  Homomorphic deconvolution estimators  i n the  text,  models of the s e i s m i c t r a c e which  are necessary i n l a t e r d i s c u s s i o n s . processing  Olrych,  preserve c o n t i n u i t y w i t h i n the  o r i g i n a l c o n t r i b u t i o n s are emphasized only  deconvolution  time.  deconvolutions  A l s o , Sven T r e i t e l ' s  length allowed me  deconvolution.  to  find  In d i s c u s s i n g  i n Chapter 3, I show t h a t  ix differences filtering  between wiener and  e s s e n t i a l l y depend on the  Following the  discussions  d e c o n v o l u t i o n s and region  maximum entropy p r e d i c t i o n s i z e of data gates.  of i n v e r s e  well l o g synthetics  are compared.  I t was  filtering,  from a complex  encouraging to note that  p l a u s i b l e i n t e r p r e t a t i o n could  error  be e x t r a c t e d  from such  geological a difficult  data. In Chapter 5, a new developed by myself and  v i b r o s e i s deconvolution method Hob  exploits certain properties cepstral analysis  and  C l a y t o n i s presented. of the  high r e s o l u t i o n  properties  approach  v i b r o s e i s t r a c e which allow  a u t o r e g r e s s i v e p r e d i c t i o n to be  Deconvolutions of both s y n t h e t i c the  The  and  of our  used,  r e a l v i b r o s e i s data technique.  exhibit  1 CHAPTEB 1  1*1  I I I I O D O C T I O N TO DECOHVOLOTIO-N-  i l J B i i i H Besfionse of a Layered  E l a s t i c E-arth  In r e c e n t y e a r s , d e c o n v o l u t i o n has r e c e i v e d much a t t e n t i o n as a means of enhancing i n f o r m a t i o n contained seismic recordings.  in reflection  In e x p l o r a t i o n seismology,  has been used to r e s o l v e g e o l o g i c a l boundaries As an i n t r o d u c t i o n to d e c o n v o l u t i o n ,  deconvolution w i t h i n the e a r t h .  the r e f l e c t i o n of a  plane compressional  wave by a s e r i e s of h o r i z o n t a l rock l a y e r s  i s considered.  case of v e r t i c a l i n c i d e n c e i s used so t h a t  shear  The  waves need not  wave r e f l e c t i o n has T r e i t e l and Claerbout  i n t o account.  T h i s model of P  been d e s c r i b e d by s e v e r a l authors i n c l u d i n g  Hobinson  (1976).  be taken  (1966),  Kanasewich  (1973),  and  A diagram of the r e f l e c t i o n s e i s m i c method  f o r a G o u p i l l a u d e a r t h model i s given i n Figure 1. modelconsists mathematically  The  of a s e r i e s of h o r i z o n t a l l a y e r s and i s convenient  equals a constant, T/2.  s i n c e the t r a v e l time f o r each l a y e r As shown i n F i g u r e 1, r e f l e c t i o n s  t r a n s m i s s i o n s of the P wave occur at each i n t e r f a c e . m u l t i l a y e r e d medium, a downward propagating seguence of pulse a r r i v a l s t o be recorded Arrivals resulting  In  and  this  pulse causes a  at the s u r f a c e .  from s i n g l e r e f l e c t i o n s of the s e i s m i c pulse  are termed p r i m a r i e s , while a r r i v a l s r e s u l t i n g from two r e f l e c t i o n s are termed m u l t i p l e s .  or more  Figure  1.  P wave propagating  through a G o u p i l l a u d non-vertical  earth model.  at v e r t i c a l  incidence  The rays are drawn at  i n c i d e n c e t o simulate t h e change of the  wave's v e r t i c a l  p o s i t i o n with time.  0"^ ,  denote  upgoing and downgoing wave amplitudes a t the top of each layer.  3  0  .  T  2T  The  r e f l e c t i o n and t r a n s m i s s i o n of a s e i s m i c wave a t an  i n t e r f a c e i s d e s c r i b e d by two boundary c o n d i t i o n s . c o n d i t i o n s s t a t e t h a t the wave's displacement s t r e s s a s s o c i a t e d with t h i s displacement boundary.  These  and the normal  are continuous  at a  From these c o n d i t i o n s , r e l a t i o n s h i p s between the  i n c i d e n t , r e f l e c t e d and t r a n s m i t t e d displacements  can be  deduced. The  r e f l e c t i o n c o e f f i c i e n t f o r a wave t r a v e l l i n g  from l a y e r  k t o l a y e r k + 1 i s given by:  where  , A  r  are the amplitudes of the i n c i d e n t and r e f l e c t e d  waves, v i s the P wave v e l o c i t y , and  i s the d e n s i t y of the  layer. The  transmission  coefficient  A  1  ftV«t  f o r t h e wave i s given by:  f> , V«„ kt  where A-,: i s the amplitude of the t r a n s m i t t e d wave. Claerbout  (1976) p o i n t s out t h a t the r e f l e c t e d wave  amplitudes can be deduced from r e f l e c t i o n c o e f f i c i e n t s , and conversely,  the r e f l e c t i o n c o e f f i c i e n t s can be deduced from wave  amplitudes.  T o : e l a b o r a t e on t h i s , the d i s c u s s i o n of s e i s m i c  wave propagation  given by-Robinson  (1967a) i s f o l l o w e d .  Using the model of F i g u r e 1, the z transform i s employed to relate  wave amplitudes  a t the top and bottom of a l a y e r .  I f the  upgoing wave a t the top of l a y e r k i s represented by the time series  ( u^(0) , u ( 1 ) , . . . , u^ {ti)) , i t s z transform i s given by K  N  u  H  I f z=e  LO) + U d) K  Z -f-  (1.2)  , the z transform equals the d i s c r e t e F o u r i e r transform.  The z transform i s a very convenient  way o f e x p r e s s i n g a time  s e r i e s s i n c e m u l t i p l i c a t i o n by z g i v e s a u n i t time d e l a y and d i v i s i o n by z s i g n i f i e s a time  advance.  I f the z transform of the downward propagating  wave at the  top of l a y e r k ' i s denoted by D (z) , a d e s c r i p t i o n of the K  r e f l e c t i o n s and t r a n s m i s s i o n s o c c u r i n g a t the kth i n t e r f a c e i s given by the f o l l o w i n g r e l a t i o n s i n z transform n o t a t i o n (Kanasewich  (1973)) .  Repeated a p p l i c a t i o n o f the above formulae  f o r n+1 l a y e r s  r e s u l t s i n the d e r i v a t i o n o f a r e l a t i o n s h i p between the i n i t i a l p u l s e sent i n t o the ground, D ( z ) , and the upgoing and downgoing 0  6 s a v e s i n any  layer.  (1.4)  is  "communication  matrix"  Assuming t h a t the  initial  that and  U 0  o  h + )  pulse  (z)=0  (z) c a n  l a y e r n+1 sent  and be  for the  D  Q  i t h layer.  i s an  i n t o the  infinite  half-space  earth i s a u n i t spike  {z) = 1-rj (z) . 0  determined  the  Hence, t h e  from the  values  values  of the  and  that  implies of  D^,  (z)  reflection  coefficients. Consideration input  and  output  reflection (Claerbout to the  flux  of  the  Kunetz e q u a t i o n ,  energy f l u x e s of t h i s  coefficients  t o be  (1976)). . Equating escaping  which  model, a l l o w s  determined  from  the  the  flux  of  through the  r e l a t e s the  energy  lower h a l f space  the  waves the  top  yields  layer  7  Y, [ (I - U (^)(l-U.i''i)  U . W U o C V ' j - Y„ D^Wnj*-;(1.5)  -  0  where Y r e p r e s e n t s t h e a d m i t t a n c e Equation  constant  f o r the layers,  (1.5) i m p l i e s t h a t t h e a u t o c o r r e l a t i o n o f t h e e s c a p i n g  wave, D ,  (z) D ,  reflected  to the surface,  nt  The  (z~*),  n+  escaping  reflection  wave,  i s directly  D  related  t o t h e wave  U (z). c  n +  coefficients.  ,j-  (z) , c a n be used  Robinson  to determine  (1967a) shows t h i s  the  by use o f  the f o l l o w i n g polynomial.  /VH)=  The  K  polynomial  K U )  Claerbout is  + , A (k) z *  1 + A (l)z  M  i s related  to the escaping  ° z"  o  (1976)  ?  f  ±,  ( 1  coefficient  A (z) i s shown t o be t h e z t r a n s f o r m K  the  circle  coefficients reflected  i n t h e complex  o f a minimum d e l a y  z plane.  of the non-neqative  p u l s e s t o be r e l a t e d  of z  in A  f o r the kth layer.  which i m p l i e s t h a t a l l t h e z e r o s  unit  wave b y :  shows t h a t t h e c o e f f i c i e n t  egual to the r e f l e c t i o n  series,  n.6)  H  K  .  7 )  (z)  Also,  time  of A {z) are outside k  Equating the  powers o f z i n (1.5) a l l o w s the  t o the r e f l e c t i o n  coefficients  by  8 the f o l l o w i n g matrix  agnation,  Since the matrix c o n t a i n i n g the r e f l e c t e d pulse values i s in  T o e p l i t z form, the a l g o r i t h m formulated  can be used to s o l v e the values o f c solution. (1.8)  n  (1.8).  by Levinson  I n a p p l y i n g the Levinson  (1947) algorithm,  are determined at each step i n the r e c u r s i v e  As w i l l be seen l a t e r , the equations c o n t a i n e d i n  a r e i d e n t i c a l t o those d e f i n i n g a p r e d i c t i o n e r r o r f i l t e r ,  where u (t)  values a r e r e p l a c e d by a u t o c o r r e l a t i o n v a l u e s , r ( t ) .  Q  In f a c t , s i n c e the values i n the T o e p l i t z matrix of (1.8) equal the a u t o c o r r e l a t i o n of the escaping  wave m u l t i p l i e d by a  c o n s t a n t , A (z) r e p r e s e n t s the p r e d i c t i o n e r r o r f i l t e r n  escaping wave. error f i l t e r  f o r the  Hence, the l a s t c o e f f i c i e n t i n a p r e d i c t i o n  i s o f t e n termed a r e f l e c t i o n c o e f f i c i e n t .  Since  the a u t o c o r r e l a t i o n of a wave recorded a t the s u r f a c e i s used i n the design of the p r e d i c t i o n e r r o r f i l t e r ,  this f i l t e r  could be  r e l a t e d to the r e f l e c t i o n c o e f f i e n t s only under the c o n d i t i o n  that the a u t o c o r r e l a t i o n of the escaping wave i s equal to the a u t o c o r r e l a t i o n of the r e f l e c t e d wave.  T h i s i s an u n r e a l i s t i c  c o n d i t i o n of s t a t i o n a r i t y f o r s e i s m i c waves which would not be s a t i s f i e d f o r the r e a l e a r t h .  Hence, i t would not be f e a s i b l e  to determine  the r e f l e c t i o n c o e f f i c i e n t s by use of p r e d i c t i o n  error f i l t e r  coefficients.  Although i t i s d i f f i c u l t t o d i r e c t l y the r e f l e c t i o n interest. termed  determine e s t i m a t e s of  c o e f f i c i e n t s , the u ( t ) time sequence i s o f r e a l 0  T h i s sequence,  which w i l l be denoted as i ( t ) , i s  the r e f l e c t i v i t y sequence or t h e impulse response of a  layered earth. determined, be o b t a i n e d .  I f the primary a r r i v a l s i n t h i s sequence can be  then a good i n t e r p r e t a t i o n of subsurface geology can F o r t u n a t e l y , the task of determining primary  a r r i v a l s i n i ( t ) can be performed  by suppressing m u l t i p l e s  through use of p r e d i c t i v e d e c o n v o l u t i o n (Peacock and Treitel  (1969) ) .  _1.2 l a v e l e t Model of a Seismic Trace In the p r e v i o u s d e s c r i p t i o n of the earth's impulse response, the i d e a l i z e d case was used i n which the s e i s m i c s i g n a l was a sequence of d e l t a f u n c t i o n s .  Because of the nature  o f s e i s m i c s o u r c e s , the s e i s m i c r e c o r d i n g system, and the a b s o r p t i o n p r o p e r t i e s of the e a r t h , a r e a l i s t i c model r e p l a c e s  10 the d e l t a f u n c t i o n s fay source generated wavelets. r e c o r d e d by seismometers  Savelets  are p r o p o r t i o n a l to the displacement  v e l o c i t y of s e i s m i c waves and r e p r e s e n t the manner i n which the e a r t h has f i l t e r e d  a source generated s i g n a l .  F i g u r e 2a)  d i s p l a y s a p o s s i b l e displacement waveform generated by an e x p l o s i o n as well as the recorded wavelet which i s i t s time derivative.  A s i m i l a r displacement waveform has been presented  by Grant and West (1965) to d e s c r i b e the r a d i a l displacement generated by a p r e s s u r e pulse w i t h i n a s p h e r i c a l  cavity.  Many of the present concepts of s e i s m i c wavelets have evolved from the work of R i c k e r (1953).  In B i c k e r ' s paper,  s e i s m i c wavelets are considered to be the fundamental the seismogram. experiments Ricker  Using s e i s m i c wave theory as well as  performed  over the P i e r r e s h a l e i n C o l o r a d o ,  (1953) formulated fundamental  propagation.  u n i t s of  laws governing  With i n c r e a s e i n propagation time, the  wavelet wavelet's  amplitude decays while i t s width i n c r e a s e s .  A s i m u l a t i o n of  t h i s i s given i n F i g u r e 2 (b).  d e f i n i t i o n of the  wavelet i s given by Robinson  A statistical  (1967b)., The wavelet i s d e s c r i b e d  as a t r a n s i e n t time f u n c t i o n of f i n i t e energy which r e p r e s e n t s the p r e d i c t a b l e part of the s e i s m i c t r a c e .  Figure  2a).  P o s s i b l e d i s p l a c e m e n t waveform and  recorded wavelet f o r a seismic Figure  2b)  wave.  V a r i a t i o n i n t h e shape o f the s e i s m i c  w a v e l e t as a f u n c t i o n o f p r o p a g a t i o n d i s t a n c e  when  travelling  geophone.  b e t w e e n t h e s h o t and t h e r e c o r d i n g  Numbers 1 t o 4 r e l a t e t o f o u r d i f f e r e n t  wavelet shapes  which are f u n c t i o n s of propogation time a f t e r the shot explosion.  V.  REFLECTING HORIZON  13 The s e i s m i c t r a c e may  be c o n s i d e r e d as a sum  delayed wavelets which have been r e f l e c t e d from interfaces.  The  of the impulse  wavelet amplitudes  response seguence.  the s e i s m i c t r a c e , x ( t ) , may a wavelet, w ( t ) ,  are weighted  of  time  subsurface by the v a l u e s  Hence, a s t a t i o n a r y model of  be expressed as the c o n v o l u t i o n of  with the e a r t h ' s impulse response, i ( t ) .  (1.9)  1*3  A t t e n u a t i o n and Dj^£ersion of a Seismic P u l s e  The time v a r y i n g nature of the wavelet  i s a r e s u l t of the  a t t e n u a t i o n and d i s p e r s i o n of s e i s m i c waves propagating the e a r t h .  Assuming a l i n e a r wave theory and a l i n e a r  of a t t e n u a t i o n with frequency, Futterman  through variation  (1962) showed that  the  d i s p e r s i o n of s e i s m i c waves i s a conseguence of a t t e n u a t i o n . the complex wavenumbar, k(oj-), i s a f u n c t i o n of frequency, the displacement  of a plane wave can be expressed  as  If  14  ix(Rd)  =  J u(o,*o e ' "  where k(w-) = K(ur)  + i <^{U)") .  components,  U.CO.^'/ e  velocities.  The f u n c t i o n  [or di-  In t h i s case,  (  1  -  1  0  )  the F o u r i e r  have d x f f e r e n t phase K ( w ) determines the phase o f each  F o u r i e r component and hence governs the wave's d i s p e r s i o n , while C*{6J-)  determines the a t t e n u a t i o n .  D i s p e r s i o n and a t t e n u a t i o n  are r e l a t e d by the Kramers - Kronig r e l a t i o n s h i p which e s s e n t i a l l y expresses The  K ( u f ) as the H i l b e r t transform  of <?C(ur) .  r e l a t i o n s h i p i s d e r i v e d by using the p r i n c i p l e of c a u s a l i t y .  A s u f f i c i e n t c o n d i t i o n f o r u(R,t) to be c a u s a l i s t h a t has no zeros i n the upper h a l f of the o& plane.  n(0,o^)  This condition  on u a l s o i m p l i e s t h a t 'Futterman - s a t t e n u a t i o n laws a r e minimum phase.  Moreover, the a t t e n u a t i o n and d i s p e r s i o n laws imply  that  the s e i s m i c wavelet i s n o n s t a t i o n a r y .  1"H lJE££°l£iiSS t o Since  wavelets contained  Deconvolution  i n t h e seismic t r a c e a r e time  v a r y i n g and s i n c e a d d i t i v e noise i s present,  a more r e a l i s t i c  but l e s s t r a c t a b l e model of the s e i s m i c t r a c e i s given by:  15  X(t)  c(ir) w(t~t  =  }  t)  (1.11)  +  Here n (t) r e p r e s e n t s the a d d i t i v e noise term,-and w(2r,t) i s convolved  with i ( t ) .  The form w(-7j,t) i s used t o i n d i c a t e  that t h e form o f - t h e wavelet f u n c t i o n v a r i e s with time  (Clarke  (1968)).  Excluding the noise term,  (1.11) r e p r e s e n t s a d i s c r e t e form of the Booton Although  propagation expression  equation.  the g e n e r a l s o l u t i o n of (1.11) i s n o t determined i n  s e i s m i c p r o c e s s i n g , steps can be taken t o reduce t h i s to  the simpler c o n v o l u t i o n form given by (1.9).  T h i s can be  done by using the s e i s m i c data p r o c e s s i n g techniques outlined  by Wood and T r e i t e l  (1975).  redundancy i n common depth p o i n t  equation  which a r e  Because of the data  (CDP) r e c o r d i n g s , normal  moveout c o r r e c t i o n s (NMO) f o l l o w e d by s t a c k i n g can decrease the n o i s e t o s i g n a l r a t i o i n the s e i s m i c t r a c e .  Also , t h e  assumption o f a s t a t i o n a r y wavelet can be j u s t i f i e d i f restricted Given  time'gates  of a t r a c e are processed.  x (t) as d e s c r i b e d by (1,9), t h e d e c o n v o l u t i o n  problem  i n v o l v e s e s t i m a t i n g the g e o l o g i c a l l y : i n t e r e s t i n g . i m p u l s e , response,  i { t ) . To do t h i s , an operator must be designed t o  remove, or deconvolve, two  fold  the wavelet from t h e t r a c e .  This i s a  problem which i n v o l v e s e s t i m a t i n g s e i s m i c wavelets as  "well as designing d i g i t a l f i l t e r s which w i l l shape t h e wavelet to a s p i k e . • • Various approaches f o r deconvolving  explosion  generated  16 s e i s m i c data are o u t l i n e d i n F i g u r e 3, taken from L i n e s and x  Ulrych  (1976).  These d e c o n v o l u t i o n techniques are examined and  compared i n the f o l l o w i n g c h a p t e r s on wavelet e s t i m a t i o n and inverse  filtering.  Figure  3.  deconvolution NMO  T h i s flow c h a r t o u t l i n e s approaches t o and  wavelet e s t i m a t i o n  = normal moveout c o r r e c t i o n and  method. applied.  Input data should  t h a t are examined. MEM  = maximum entropy  have s t a t i c s c o r r e c t i o n s  SOME  APPROACHES  Input:  Optiona 1 Preprocessing: bandpass  CDP  filter, N M O ,  I reflection!  TO  DECONVOLUTION  MEM time adaptive d econvolution  s t a c k , gain appi ication  data Single  Division into  of  data  gates  ion  Option to apply N M O and sIac k  1 Single channel Wiener spiking of wavelets  1  Estimation of wavelets using: Homomorphic decon v o lu t i o n or Wold - Kolmogorov factorization  DECONVOLVED Multichannel Wiener or K a l m a n deeo nvoln t ion W i e n e r • LevinsonJ double  WAVELET  channel  Wiener or M E M deconvplut  ESTIMATES  inverse  TRACES  19  CHAPTER 2 Two One  SEISMIC HAVEIET ESTIMATION  g e n e r a l approaches f o r wavelet e s t i m a t i o n a r e examined.  approach uses the assumptions  a random impulse response. use these assumptions O'Brien  (1975).  assumptions  of a minimum phase wavelet and  V a r i o u s e s t i m a t i o n t e c h n i g u e s which  have been examined by White and  A second approach, which does not r e q u i r e the  of the minimum phase approach, i s homomorphic  d e c o n v o l u t i o n , which has been s u c c e s s f u l l y applxed by (1971) and U l r y c h wavelets.  et a l (1972) to the e s t i m a t i o n of  Ulrych teleseismic  However, apart from the work of Buhl e t a l  O t i s and Smith  (1974), and Buttkus  (1974),  (1975), t h e a p p l i c a t i o n s of  t h i s technique to r e f l e c t i o n seismology have not been e x p l o r e d . The Miener-Levinson and wold-Kolmogorov minimum phase e s t i m a t i o n t e c h n i q u e s are i n v e s t i g a t e d here and are compared with the homomorphic d e c o n v o l u t i o n method f o r e x p l o r a t i o n , seismology  problems.  2•I The S i e n e r - L e V i n s o n Double Wiener f i l t e r i n g  Inverse Method  i s a very common method f o r performing  d e c o n v o l u t i o n , and i s o c c a s i o n a l l y used f o r wavelet estimation,. The theory concerning Wiener f i l t e r s  and t h e i r design has been  thoroughly d e s c r i b e d by Robinson and T r e i t e l filter,  f (t) , shapes a time sequence,  approximates In  the f i l t e r  (1967).  The  x ( t ) , i n t o an output  Wiener which  the d e s i r e d output, z ( t ) , i n a l e a s t squares sense. design, the mean squared e r r o r given by  20 (2.1)  -t  i s minimized by choice o f f i l t e r  2B  for  ^  values.  That i s .  0  (2.2)  s=0,1,2,...,m. This yields  t  u  (2.3)  t  Assuming s t a t i o n a r i t y , we can r e p l a c e t - s by t t o get  q(s)  -  £  +  In the above, r (s) = JJ  (ol]  rls-u)  x (t+s) x (t)  a u t o c o r r e l a t i o n of the i n p u t , cj (s) =  (2.4)  i s the z (t+s) x(t)  i s the  •t  c r o s s c o r r a l a t i o n between the i n p u t and the d e s i r e d output, and N i s the l e n g t h of the input sequence. presented  The normal  i n (2,4) a r e u s u a l l y w r i t t e n i n matrix  equations form as  21  r  r  (O)  r(-D  (I)  rCO)  ...  rC~m) : (2.5)  Cm)  r(0)  For our purposes,  we design a Wiener i n v e r s e f i l t e r  which  t r i e s t o shape the wavelet t o a s p i k e so that c o n v o l u t i o n of the filter  with t h e s e i s m i c t r a c e w i l l provide an estimate of the  impulse  response.  I f two assumptions are made, a Wiener d e c o n v o l u t i o n  filter  can be designed without e x p l i c i t l y e v a l u a t i n g the source wavelet.  I f the impulse  sequence, the wavelet trace autocorrelation.  response  i s assumed to be a white  noise  a u t o c o r r e l a t i o n can be r e p l a c e d by the T h i s can be seen  by w r i t i n g down the  t r a c e a u t o c o r r e l a t i o n i n z transform n o t a t i o n and r e c o g n i z i n g t h a t the z transform of the c o n v o l u t i o n o f two f u n c t i o n s equals the product o f t h e i r z transforms. Chapter 1,  Using the n o t a t i o n of  22  W(z)  Since  I(z)I(z~')=1  (2.6)  (z")  f o r i ( t ) being a white noise seguence, the  t r a c e a u t o c o r r e l a t i o n can be if  I  WW)  the d e s i r e d output  used i n the f i l t e r  design.  i s s e t to be a spike at zero d e l a y ,  the f i r s t  member of the c r o s s c o r r e l a t i o n vector i n the  equations  w i l l be nonzero.  be optimum only i n cases T r e i t e l and Sobinsoh  The  output  only  normal  zero delay s p i k i n g p o s i t i o n w i l l  where the wavelet i s minimum  (1966) show that f i l t e r  delay.  performance  improve i f d e l a y s i n the spike p o s i t i o n of the d e s i r e d are allowed.  Also,  will  output  However, to use a delay of n u n i t s i n the d e s i r e d  r e q u i r e s that we have a good estimate  of the f i r s t  n  p o i n t s of the s e i s m i c wavelet to determine (9 (0) i 9 O) / • - • #9 (a) , 0,. . . ,0)  i n the normal  equations.  To o b t a i n an a c t u a l wavelet e s t i m a t e , another Wiener i s a p p l i e d to i n v e r t the o r i g i n a l minimum delay i n v e r s e Consequently, t h i s wavelet e s t i m a t i o n technique Wiener-Levinson double i n v e r s e method.  filter  filter.  i s termed the  In using t h i s  approach,  t h e r e e x i s t s the problem of using the r e s t r i c t i v e minimum phase and  randomness assumptions as well as the problem of  a s u i t a b l e l e n g t h f o r the Wiener The  Wiener i n v e r s e f i l t e r  filter.  i s designed  a u t o c o r r e l a t i o n of the t r a c e i n the f i r s t e s t i m a t i o n process.  determining  d i r e c t l y from step of the  the  wavelet  Therefore, the Wiener-Levinson double  i n v e r s e method i s r a r e l y used i n p r a c t i c e , s i n c e i t i s the inverse f i l t e r ,  not the wavelet, which i s u s u a l l y d e s i r e d .  23 2.2 l o l d - K o l r a o g o r o v F a c t o r i z a t i o n The  Wold-Kolmogorov f a c t o r i z a t i o n  described  in detail  phase s p e c t r u m amplitude  by R o b i n s o n  amplitude spectrum.  and i t s s p e c t r u m  approach. wavelet  T h i s approach  i s equivalent  correlation First  into  t h e minimum  G i v e n an n  which  will  have  this  one o f t h e s e w a v e l e t s i s minimum  may be d e d u c e d  by t h e f o l l o w i n g  o f d e t e r m i n i n g t h e minimum  to f a c t o r i n g  minimum  spectrum.  phase s p e c t r a  Only  determines  o f l e n g t h n+1, t h e r e a r e 2  f o r a wavelet  wavelets with d i f f e r e n t  phase,  (1967b),  f o r a given amplitude  spectrum  method, M h i c h i s  phase  t h e z transform of the auto-  phase b i n o m i a l f a c t o r s .  of a l l the l o g of the wavelet's F o u r i e r transform,  a ( o/) , i s g i v e n by  (2.7)  where  |W(u>-) | i s t h e a m p l i t u d e  spectrum. Fourier  For a causal  wavelet,  spectrum  and  O^iuf)  i s t h e phase  l o g (W ( ur) ) i s g i v e n by a  expansion:  (2.8)  As  p o i n t e d o u t by R o b i n s o n ,  condition that  (2.8) i s d e r i v e d  l o g ( w (z)) be a n a l y t i c  inside  under the  the u n i t  circle i n  24 the complex transform,  plane.  T h i s r e q u i r e s that the wavelet's z  I (z), have a l l i t s r o o t s o u t s i d e the unit  circle  which i s the c o n d i t i o n of minimum phase. The F o u r i e r s e r i e s f o r the even f u n c t i o n , loglW(o-) I, i s given by  oO =  /talWMl  oC(0)  +  2  COS (St  g  (2.9)  where  oc(±)  =  J  TT  Tr  /ojl  WCA)|  COS  Jl-t  Comparison of ( 2 . 8 ) and ( 2 . 9 ) shows that t>1,  #{0)=c<(0).  and  M  (2.10)  # (t) = 2 ^ (t) f o r  A l s o , comparison of ( 2 . 7 ) and ( 2 . 8 ) then  shows that oO  &(<*) =  ~  2  Z  ^  (2.11)  M  Hence, the minimum phase spectrum,  &(kr),  i s given i n  terms of the amplitude spectrum, J W ( ur) | , by the f o l l o w i n g relation.  25 IT  &(ur)  =  ~2  J T t=>  1 7  I t can r e a d i l y  O'(u-)  est  sin  v/  0 (&)  transform,  h j lw((*r/}  with  cosM  ** *  fo IW(JOldfl  2  S  12  0  be shown t h a t the above e x p r e s s i o n f o r  d e r i v e d by Robinson  the l o g amplitude  f  spectrum.  (1967b) i s the H i l b e r t transform of First  of a l l , s i n c e the H i l b e r t  , i s essentially  the c o n v o l u t i o n of  , t a k i n g the F o u r i e r transform of  — (f)(ur)  y i e l d s the f o l l o w i n g ,  F(<J>) =  - fe"  Since and  TT ,  '  M  f^  I'}  M (  i s an even f u n c t i o n d e f i n e d between  2.i3,  —IT  (2.13) can be r e w r i t t e n as: TT  F(0)  =  2i  j  C O S A *  loj\VJLJl)\dJl  > S r>Lt) 3  (2.14)  Inverse F o u r i e r t r a n s f o r m a t i o n g i v e s the o r i g i n a l value of the H i l b e r t t r a n s f o r m .  The transform i s d i s c r e t e because we are  d e a l i n g with d i g i t a l values i n the time domain. given by:  Hence, (f)(<*) i s  26 1*9  (2.15)  Since sgn (t)  i s an odd f u n c t i o n of t , the above may  w r i t t e n i n the same form as  (2.12).  be  Hence, the Wold -  Kolmogorov f a c t o r i z a t i o n method d e s c r i b e d by Robinson  (1967b) i s  o f t e n termed the H i l b e r t transform method. The wavelets o b t a i n e d by use of the H i l b e r t  transform are  o f t e n q u i t e s i m i l a r to those o b t a i n e d by using the Wiener Levinson double i n v e r s e method s i n c e both methods are determined by using an a u t o c o r r e l a t i o n estimate and the minimum assumption.  Galbraith  (1971) p o i n t s out that the methods would  be e q u i v a l e n t i f an i n f i n i t e number of f i l t e r used i n the Wiener  points could  be  filter.  The a n a l y s i s given-here i m p l i c i t l y wavelet as a t r a n s i e n t time sequence f o l l o w i n g Robinson  phase  (1967b) .  uses the d e f i n i t i o n of a  with f i n i t e  energy,  In a p p l y i n q Wiener-Levinson wavelet  e s t i m a t i o n t h i s d e f i n i t i o n r e q u i r e s the c h o i c e of a p p r o p r i a t e filter  l e n g t h s , whereas the wavelet obtained by using the  H i l b e r t transform must be •truncated a t some s u i t a b l e  point.  The minimum phase wavelet e s t i m a t o r s examined here have used a s t a t i s t i c a l d e f i n i t i o n of the impulse response as being a white noise sequence.  Recently, Robinson  (1975b) has used a  dynamic  modal f o r the impulse response fay a p p l y i n g the wave  theory presented i n Chapter 1.  2.3  S2!°!9J.EhiLc: Deconvolution  Homomorphic d e c o n v o l u t i o n , which was i n t r o d u c e d by Oppenheim imposed  (1965), does not r e q u i r e t h e r e s t r i c t i v e  conditions  by the minimum phase methods of wavelet e s t i m a t i o n .  In usinq homomorphic deconvolution, a c h a r a c t e r i s t i c t r a n s f o r m a t i o n i s a p p l i e d to a time sequence x ( t ) , i n order to y i e l d the complex cepstrum.  The complex  cepstrum of as time  sequence i s d e f i n e d as the i n v e r s e F o u r i e r transform o f the l o g a r i t h m of the time sequence's F o u r i e r transform.  For a  d i g i t a l l y sampled seguence, x ( t ) , t h e complex cepstrum, x(t)  is  given by:  The terms "quefrency" and cepstrum ,,  Bogert e t a l (1963).  ,,  were i n t r o d u c e d by  The term "quefrency" d e s c r i b e s the  " f r e q u e n c i e s " or r e p e t i t i o n r a t e s of a frequency domain and the complex  sequence  cepstrum i s e s s e n t i a l l y a measure of the  quefrency content o f the l o g c f a sequence's F o u r i e r C a l c u l a t i o n of the complex  transform.  cepstrum transforms sequences  28 which are convolved a d d i t i v e i n a new The  i n the time domain i n t o sequences which are  domain termed the quefrency  complex cepstrum of the c o n v o l u t i o n a l model f o r the  s e i s m i c t r a c e w i l l be c o n s i d e r e d . V  x(t) =  domain.  w(s)  i(t-s)  yields  The  Fourier transform  X(^)  ~  U/(w)  of  J.  Taking the l o g a r i t h m g i v e s  laj  X((JS)  -  /oj  W(w)  +  The i n v e r s e F o u r i e r transform of  x(-t)  =  where w (t) and i (t) impulse complex  response  ~h  1(UT)  (2.17) then  (2.17)  produces  tit)  (2.18)  are the c e p s t r a of the wavelet and  as d e f i n e d by  (2.16).  the  The computation of the  cepstrum i n v o l v e s c e r t a i n e s s e n t i a l steps which are  o u t l i n e d by Schafer Stoffa  wit)  /oj  et a l  (1969),  Ulrych  (1971),  and  (1974),  F i r s t , the e x p r e s s i o n f o r the l o g of the F o u r i e r transform of  x (t)  i s considered.  29 (2.19)  To make the phase spectrum,  0*(ux) , an a n a l y t i c  continuous  f u n c t i o n of freguency r e q u i r e s t h a t the raw phase values computed by the Arctan f u n c t i o n be unwrapped by adding m u l t i p l e s of  2TT. The l i n e a r component should a l s o be removed from the  phase s i n c e i t may obscure the i n t e r e s t i n g i n f o r m a t i o n i n the cepstrum.  Removal o f the l i n e a r phase component, bur , from the  phase, units.  , simply s h i f t s the o r i g i n time of the t r a c e by b A l s o , the mean o f l o g X(or)  i s u s u a l l y removed.  This  o p e r a t i o n s e t s x(0) to zero and i s e q u i v a l e n t to m u l t i p l y i n g the o r i g i n a l time sequence of  l o g X(cj-).  by a constant exp(-a)  where a i s the mean  Hence, the t=0 value i n the complex  cepstrum i s  only a f f e c t e d by t h e s c a l i n g i n x (t) . As seen by (2.18), sequences  which a r e convolved i n the  time domain are a d d i t i v e i n the guefrency domain and an estimate of  the wavelet cepstrum can be obtained by l i f t e r i n g  (i.e. linear f i l t e r i n g cepstrum  i n the guefrency domain).  i s o f t e n estimated by low c u t l i f t e r i n g .  wavelet cepstrum i s s e t egual to complex  The wavelet That i s , the  cepstrum v a l u e s f o r  J t | l e s s than some c u t o f f value. The s e p a r a t i o n o f the wavelet cepstrum  and the impulse  response cepstrum i s p a r t i c u l a r l y convenient when the impulseresponse i s minimum phase. form of i (t) .  T h i s can be seen by c o n s i d e r i n g the  30  The l o g of the F o u r i e r t r a n s f o r m of i ( t ) i s given by  Lit) Lit) ee  If  \2L  J  e  l o g ( I ( u r ) ) can be used.  < 1, a Laurent expansion f o r Thus,  "t) e — ; = r  /o,  (2.20)  j  i«»  ( 2  e  <2.2i,  A  The value of i ( t ) i s given by r  2lT  fc+/  OO  o"  l  f  r < » ^ T  (2. 22)  or  ( 2 . 23)  -IMS  Hence, i f 1 ^ ^ s ) e  | i s much l e s s than one, the l a r g e s t  s p i k e s i n the impulse response w i l l correspond t o those i n the o r i g i n a l sequence,  i (t) .  T h i s c o n d i t i o n on I ( c j ) ensures  minimum phase and can be imposed by e x p o n e n t i a l weighting.  31 E x p o n e n t i a l weighting i s performed fay m u l t i p l y i n g terms i n t sequence x{t) by °< r e p l a c i n g z=e effectively  , where 0 < ° < < 1 . 0 .  This i s equivalent to  x.n the z t r a n s f o r m , X (z),  moves the zeros of X (z)  by z=e  J  , and  outward from the o r i g i n .  E x p o n e n t i a l weighting makes the impulse response appear t o be minimum phase and x { t ) = x ( t ) .  T h e r e f o r e , i f the f i r s t  and second  a r r i v a l s c o n t a i n e d i n a minimum phase impulse response are a t t=0 and t=s the  r e s p e c t i v e l y , the f i r s t  s p i k e s i n the cepstrum of  impulse response are a t t=0 and t=s  also.  An i d e a l  d e c o n v o l u t i o n can then be performed i f the wavelet cepstrum i s localized of  w i t h i n t=s, s i n c e i n such cases a complete  separation  the wavelet cepstrum and the cepstrum of the impulse response  can be made. As emphasized  by S t o f f a et a l (1974), e x p o n e n t i a l weighting  can a l s o be used t o e l i m i n a t e c e p s t r a l a l i a s i n g caused by the n o n l i n e a r o p e r a t i o n s of the l o g a r i t h m , a b s o l u t e value and the arctangent. A f t e r the cepstrum has bean c a l c u l a t e d and l i f t e r e d an estimate of the wavelet cepstrum, the i n v e r s e  to give  characteristic  t r a n s f o r m a t i o n i s a p p l i e d to y i e l d the source wavelet e s t i m a t e . T h i s i n v e r s e c h a r a c t e r i s t i c t r a n s f o r m a t i o n i s performed by t a k i n g the F o u r i e r t r a n s f o r m , the e x p o n e n t i a l and the i n v e r s e Fourier transform. the  I f we denote the wavelet cepstrum by w(t),  wavelet estimate i s given by  32  ^  { iurt +  2  w(t)  Ws)}  e  S=-oa  Although the cepstrum  (2.24)  of the impulse response c o u l d  transformed d i r e c t l y - to o b t a i n a deconvolved difficult because may  to o b t a i n the a c t u a l cepstrum  trace, i t i s often  of the impulse  the impulse response s p i k e s c o n t a i n e d i n the  be obscured by high guefrency n o i s e .  cepstrum  low c u t l i f t e r  cepstrum  cepstrum  wavelet o b t a i n e d from t h i s cepstrum  2.4  i n the low  guefrency  of the t r a c e . .  The  can then be used i n the  t o deconvolve the t r a c e .  Comparisons of Wavelet Estimators The  performance  of the p r e v i o u s l y mentioned  wavelet  e s t i m a t i o n t e c h n i q u e s i s compared by using r e a l and data.  cepstrum  and can be obtained by a p p l y i n g a  to the complex  design of an i n v e r s e f i l t e r  response  On the other hand, the  of the wavelet i s o f t e n l o c a l i z e d  p a r t of the complex  be  synthetic  To the author's knowledge, no other comparison  wavelet e s t i m a t o r s has appeared  of these  i n e x p l o r a t i o n seismology  literature. F i g u r e 4 compares wavelet e s t i m a t i o n methods by synthetic trace.  The t r a c e i s formed  sequence with a wavelet. the author by G.L,  by c o n v o l v i n g a spike  T h i s wavelet, which  Cumming, was  using a  was suggested  s y n t h e s i z e d by using a z  to  transform of the form average,  (-1,10 + z ) * (1.75+ z)'" and removing the 4  The r e s u l t i n g time sequence s i m u l a t e s a p o s s i b l e  e x p l o s i o n generated wavelet, s i m i l a r to those found by Ricker  (1953).  Because of removal o f the average, the wavelet  i s not q u i t e minimum phase. not  A l s o , s i n c e the impulse response i s  a white noise sequence, the a u t o c o r r e l a t i o n used i n the  Wold-Kolmogorov and Wiener-Levinson wavelet estimates c o n s i s t s of  the t r a c e a u t o c o r r e l a t i o n m u l t i p l i e d by a Parzen window of  l e n g t h 20.  Although the windowed t r a c e a u t o c o r r e l a t i o n i s a  good estimate of the source wavelet a u t o c o r r e l a t i o n , the minimum  phase assumption  g i v e s both wavelet estimates a d i f f e r e n t  c h a r a c t e r from that o f the a c t u a l wavelet. The homomorphic d e c o n v o l u t i o n estimate o f the wavelet was obtained by low c u t l i f t e r i n g complex  the cepstrum.  In t h i s case the  cepstrum of the t r a c e was c a l c u l a t e d a f t e r using an  e x p o n e n t i a l weighting of 0.98 on the t r a c e .  As seen from the  c e p s t r a of F i g u r e 4, the a c t u a l source wavelet cepstrum be completely separated from the complex by a low cut l i f t e r .  cepstrum of the t r a c e  Hence, the t h i r d lobe o f the homomorphic  estimate i s not q u i t e the same as t h a t of the a c t u a l although the f i r s t duplicated.  wavelet  two l o b e s of the wavelet are v i r t u a l l y  Such t e s t s on n o i s e l e s s s y n t h e t i c s provide  encouragement f o r the use of c e p s t r a l methods i n wavelet estimation.  cannot  34  Figure found  4.  A comparison  of source wavelet  by t h e Wold-Kolmogorov f a c t o r i z a t i o n  (WKFACT),  the Wiener-Levinson  double  estimates  method  i n v e r s e method,  t h e method o f homomorphic d e c o n v o l u t i o n .  and  35  C O M P A R I S O N OF WAVELET E S T I M A T O R S SOURCE WAVELET  IMPULSE  —i  i  0.0  SOURCE AUTOCORRELATION  SOURCE CEPSTRUM  0.2  RESPONSE  1  1  1  1  0.4  0.6  O.B  1.0  TIME  TRACE . AUTOCORRELATION  TRACE CEPSTRUM  WKFflCT ESTIMATE  WIENER-LEVINSON DOUBLE INVERSE  HOMOMORPHIC ESTIMATE  For  the case of zero a d d i t i v e n o i s e , the performance of  homomorphic d e c o n v o l u t i o n as a wavelet estimator depends on w e l l the complex the  complex  cepstrum of the wavelet can be separated from  cepstrum of the t r a c e by low cut l i t t e r i n g .  S i m i l a r l y , the performance of the minimum phase depends cn how by windowing  how  approaches  w e l l the wavelet a u t o c o r r e l a t i o n can be estimated  the t r a c e a u t o c o r r e l a t i o n .  A d d i t i v e noise can cause major problems when u s i n g homomorphic d e c o n v o l u t i o n because i t i s d i f f i c u l t to determine i t s e f f e c t on the complex  cepstrum.  T h i s i s because  c a l c u l a t i o n s involve nonlinear transformations.  cepstral  Furthermore,  a d d i t i v e noise i n the s e i s m i c t r a c e can cause severe i n s t a b i l i t i e s when unwrapping demonstrated i n Figure 5.  the phase curve.  The n o i s e l e s s t r a c e shown i n the  f i g u r e has a maximum amplitude of 1.0. of  This problem i s  Different  "white n o i s e " with a standard d e v i a t i o n o f 0.3  this trace.  realizations were added to  A f t e r the usual procedures o f phase unwrapping  removal of the l i n e a r phase component, the phase c u r v e s show marked  differences.  and  Figure  5. These examples demonstrate t h a t d i f f e r e n t  noise r e a l i z a t i o n s added to a t r a c e can produce significant  d i f f e r e n c e s i n the unwrapped phase curve,  "white n o i s e " r e a l i z a t i o n s , 0.3, 1.0.  with standard  d e v i a t i o n s of  were added to a t r a c e with a maximum amplitude o f  38  TRACE  i.e NOISY TRACE WT=1.00  NOISY TRACE WT=1.00  NOISY TRACE WT=1.00  PHASE  Clayton and and  Wiggins  (1976) have i n v e s t i g a t e d t h i s problem  have shown t h a t d e v i a t i o n s i n the phase caused by  additive  n o i s e become severe when s i g n a l amplitudes are s m a l l . demonstrate t h a t 2 Tf  e r r o r s i n the phase unwrapping can  occur when the amplitude T h e r e f o r e , from  spectrum  They a l s o easily  values approach z e r o .  the work of C l a y t o n and Wiggins (1976) and  author's e x p e r i e n c e , i t i s suggested  the  that data redundancy be  used to suppress the e f f e c t s of n o i s e on the phase c u r v e s when using homomorphic d e c o n v o l u t i o n . Recently, O t i s and Wiggins  Smith  (1974) and Clayton  and  (1976) have devised methods of l o g s p e c t r a l  averaging  which decrease the e f f e c t s of n o i s e on the phase c u r v e .  These  methods of c e p s t r a l averaging h o p e f u l l y w i l l enhance the complex  cepstrum  of the wavelet.  V a r i a t i o n s of c e p s t r a l s t a c k i n g  methods are a p p l i e d to the windowed p o r t i o n of the f i v e t r a c e s shown i n F i g u r e 6.  Traces from the same shot are used i n order  to maintain redundancy i n wavelet  shape.  T h i s example i s  p a r t i c u l a r l y i n t e r e s t i n g because the f i f t h t r a c e c o n t a i n s c o n s i d e r a b l y more a d d i t i v e n o i s e than the other t r a c e s .  40  Figure 6 .  A s e t of f i v e t r a c e s which are used i n the  comparison of wavelet e s t i m a t o r s . t r a c e was  The windowed p o r t i o n of the  used i n the c a l c u l a t i o n of wavelets.  41 For c e p s t r a l computations, the  windowed p o r t i o n s of  t r a c e s of Figure 6 were e x p o n e n t i a l l y weighted using f a c t o r of 0.96.  F i g u r e 7 shows the  phase curves of  a  the weighting  the  e x p o n e n t i a l l y weighted t r a c e s , the  wavelets obtained  homomorphic d e c o n v o l u t i o n ,  of the complex c e p s t r a ,  the  wavelet obtained  Due  to the e f f e c t s of n o i s e , the phase curve of the f i f t h  i s considerably  by low  the sum  by  cut l i f t e r i n g of the c e p s t r a l sum  d i f f e r e n t from those of the other  Averaging the c e p s t r a corresponds to averaging and  phase s p e c t r a before  T h i s can  c e p s t r a f o r N channels which i s given  2.TT where X ( u r - ) K  -  .  trace  traces.  the l o g amplitude  t a k i n g the i n v e r s e F o u r i e r  be seen from the expression  ir  and  transform.  f o r the average of  the  by  W M  TT  denotes the F o u r i e r transform  of x  (t) .  I  2  -  2  5  )  Figure 7.  The use of c e p s t r a l averaging i n wavelet  e s t i m a t i o n i s shown using the t r a c e s of Figure 6.  This  f i g u r e shows the phase curves of these t r a c e s , the wavelets obtained by l i f t e r i n g  the c e p s t r a o f the t r a c e s ,  the c e p s t r a l average, and the wavelet obtained by a p p l y i n g a low c u t l i f t e r o f width 60 msec t o the c e p s t r a l average.  Scale f o r c e p s t r a l average and average  wavelet i s i n seconds.  43  PHASE  CURVES  HOMOMORPHIC DECONVOLUTION WAVELETS  fa  fafa  V— * r  H  -125 hz  125 hz FREQUENCY  C E P S T R A L  A V E R A G E  -0.4 S O U R C E  0.2  r  0.0  W A V E L E T  0.4  200 msec  44  Rewriting  (2.25)  i n the f o l l o w i n g form shows t h a t averaging  the  l o g amplitude spectra, i s e q u i v a l e n t t o f i n d i n g t h e l o g of  the  geometric mean of the amplitude  spectra.  (2.26)  When f i n d i n g the c e p s t r a l average, i t i s o b v i o u s l y l e s s expensive to average the l o g amplitude and phase s p e c t r a before t a k i n g the i n v e r s e F o u r i e r transform than to average the c e p s t r a themselves.  In t a k i n g a c e p s t r a l average, we are e s s e n t i a l l y  using an average of the unwrapped phase. Wiggins  However, C l a y t o n and  (1976) have proposed a method o f averaging t h e raw phase  values before unwrapping. greater s t a b i l i t y unwrapped  T h i s produces a phase curve with  than that obtained through averaging the  phase.  Figure 8(a) d i s p l a y s the l o g amplitude spectrum mean removed as w e l l as an unwrapping of the averaged values.  The corresponding complex  liftered  cepstrum  cepstrum  used i n wavelet computation  with i t s raw phase  and the low cut are a l s o shown.  F i g u r e 8(b) shows that s e t t i n g the l o g amplitude and phase s p e c t r a to zero when the s i g n a l ' s amplitude spectrum  becomes  l e s s than a c e r t a i n value has a n e g l i g i b l e e f f e c t on the cepstrum. raw  In F i g u r e 9, the wavelet obtained by averaging the  phase and l o g amplitude s p e c t r a i s compared with wavelets  o b t a i n e d by use of the H i l b e r t t r a n s f o r m .  From F i g u r e s 7 and 9,  45 it  can be seen that wavelets estimated by c e p s t r a l averaging  methods are d e f i n i t e l y mixed phase, as compared to the minimumphase wavelets obtained by use of t h e H i l b e r t t r a n s f o r m .  It  should be noted that the minimum phase wavelets shown i n F i g u r e 9 are delayed one p o i n t s i n c e the f i r s t not zero.  sample o f the wavelet i s  T h i s property, which i s common to a l l minimum phase  wavelets, i s not p h y s i c a l l y r e a l i s t i c .  F i g u r e 10 d i s p l a y s  d e c o n v o l u t i o n s which r e s u l t from the a p p l i c a t i o n of Wiener filters  designed t o s p i k e the estimated wavelets.  Although the  obvious r e f l e c t i o n s are shown on a l l three s e t s of deconvolved t r a c e s , d e c o n v o l u t i o n s using the wavelet estimated from of the raw phase and l o g amplitude s p e c t r a g i v e s l i g h t l y  averages better  r e s o l u t i o n of o v e r l a p p i n g r e f l e c t e d p u l s e s , e s p e c i a l l y between times of 1.5 and 1.7 seconds on the near v e r t i c a l  reflections.  F i g u r e 8.  a) E s t i m a t i o n of the wavelet using l o g  s p e c t r a l averaging f o r the t r a c e s of Figure 6.  The  average of the l o g amplitude s p e c t r a i s shown as well as the  phase curve obtained by unwrapping  the  raw phase curve.  was used.  An e x p o n e n t i a l weighting o f 0.98  The phase and l o g amplitude averages were used  to compute the cepstrum. cepstrum  the averages o f  The l i t t e r e d  v e r s i o n o f the  (also shown) i s used t o produce the wavelet  estimate.  F i g u r e 8b) shows that s e t t i n g the phase and  amplitude s p e c t r a t o zero a f t e r a c e r t a i n freguency produces a very s i m i l a r c e p s t r a t o t h a t obtained by using the  e n t i r e frequency range.  shown i s a l i t t e r e d  In b ) , the second cepstrum  v e r s i o n of the f i r s t  cepstrum.  D i f f e r e n c e s i n the wavelets o b t a i n e d by using the s p e c t r a i n a) and b) c o u l d not be determined by v i s u a l inspection.  The wavelet obtained i s shown i n Figure 9.  a)  LOG S P E C T R A L " RVERRGING  LOG RMP o. •r.  SPEC  PHRSE CURVE  ft  UJ  a  o _lo 0 0  6  ^ \ j  COMPLEX  1  125 hz  "  COMPLEX  CEPSTRUM  0.4  LOG S P E C T R A L LOG RMP  CEPSTRUM  0.4  TIME  b>  5  SPEC  RVERRGING  PHASE CURVE  UJ  Vi  <x X  a.  i  125hz  COMPLEX  CEPSTRUM  125 hz  COMPLEX  -0.4  CEPSTRUM  0.4  HILBERT TRANSFORM WAVELETS  i  — j 200  F i g u r e 9.  Comparison 'of the wavelet obtained  s p e c t r a l averaging transform  msec  with  wavelets obtained  by  log  by use of the  method to the t r a c e s shown i n F i g u r e  6.  Hilbe  Figure  10,  Comparison of  deconvolutions  t r a c e s o f F i g u r e 6 u s i n g a Wiener f i l t e r msec t o s u p p l y a)  the  an i n v e r s e f i l t e r  wavelets  of  the  of l e n g t h  80  for;  o b t a i n e d fay Wold-Kolmogorov  factorization b)  the  wavelet  obtained  by  cepstral  averaging  c)  the  wavelet  obtained  by  log spectral  averaging,  \  50  DECON WITH WKFRCT WAVELETS  Tt  eo r  »~i  7*  1.4  TIME.  i.t  —r  2.0  2.j  1  1  1  1  a.<  s.t  J.i  »•»  k)  DECON WITH CEP RV WAVELETS  —4'  .5  H  E  K  K  £  n  DECON WITH LOG SPEC RV WRVELET  "  "  "  "  "  51 In d i s c u s s i n g the r e l a t i v e merits of the e s t i m a t i o n schemes, i t should be mentioned that despite  restrictive  assumptions used i n i t s design, the H i l b e r t transform r e q u i r e s only an estimate of the wavelet produce a reasonable  wavelet.  method  l e n g t h i n order to  On the other hand, the  homomorphic d e c o n v o l u t i o n method may r e q u i r e more "parameter j u g g l i n g " s i n c e one must choose an a p p r o p r i a t e e x p o n e n t i a l weighting  and a s u i t a b l e low cut l i f t e r .  deconvolution  A l s o , homomorphic  r e s u l t s are more dependent on the c h o i c e of data  gates than the H i l b e r t transform presents a problem. homomorphic wavelet  method, and t h e . a d d i t i v e noise  These f a c t o r s g e n e r a l l y mean t h a t estimates are more expensive  to o b t a i n but  can be more s e n s i t i v e to the data when the noise l e v e l i s low and  the source generated  wavelet  i s not minimum phase.  Having d i s c u s s e d d i f f e r e n t approaches f o r e s t i m a t i n g wavelets,  a t t e n t i o n i s now turned to the second h a l f of the  d e c o n v o l u t i o n problem, the design of i n v e r s e f i l t e r s these  wavelets  to a spike.  which shape  52 CHAPTER 3  THE  3.1,  DESIGN OP  S i n q l e Channel Wiener  Wiener f i l t e r i n g  represents  approaches to d e c o n v c l u t i o n equations d e f i n i n g the outlined previously. filters  INVERSE FILTERS FOR  i s now  one  Deconvolutioi].  of the most common  i n exploration  design  and  Although shaping f i l t e r s  T r e i t e l and  Robinson  Wiener f i l t e r s spike  may  vere  discussed. choosing a spike as a  output allows the Wiener deconvolution  designed.  The  a p p l i c a t i o n of these  Osing a wavelet estimate as i n p u t and desired  geophysics.  s i n g l e channel Wiener f i l t e r The  M&VELETS  (1966) and  which shape the  filter  are not discussed  DeVoogd  here,  wavelet to some form other than a  filters  a l s o p r e s e n t s the problem of d e c i d i n g  desired  output.  of t h e i r s i m p l i c i t y , but  be  (1974) have shown t h a t  a l s o prove u s e f u l f o r deconvolution.  Spiking  to  filters  Use on  are considered  of  shaping  a suitable not  only  because  a l s o because of t h e i r i n t e r e s t i n g  r e l a t i o n s h i p to a u t o r e g r e s s i v e  models and  prediction error  filters. In assuming minimum phase, the obtained  d i r e c t l y from estimates of the  Peacock and inverse filtering  Wiener s p i k i n g f i l t e r can  Treitel  wavelet a u t o c o r r e l a t i o n .  (1969) r e l a t e the minimum pnase Wiener  f i l t e r to a p r e d i c t i o n e r r o r f i l t e r . i s examined here s i n c e i t provides  l a t e r d i s c u s s i o n of the  Burq  algorithm.  be  Prediction error the  basis for a  53  If by t h e  the u n i t d i s t a n c e p r e d i c t i o n seguence { ^ ( l ) ,  the seguence, x ( t ) ,  filter  o ( ( 2 ) , . » . , o£ (m) ) , m  m  i s d e n o t e d by  of l e n g t h  m is  given  a predicted value  of  :  -t-l  If  the p r e d i c t i o n  sequence, x ( t ) , +  o^  m  (in) x (t-m)  e r r o r s are regarded as a white  w r i t t e n a s x (t)-= + a  / > 1  (t)  (1) x ( t - 1 ) + < ^ { 2 ) x ( t - 2 )  The p r e d i c t i o n  with t h i s  filter  i s given  error  m  ...  seguence a s s o c i a t e d  by:  = X U ) ~ £ <*«(•<:-*) ><M  a (i)  +  n>  , c a n be d e s c r i b e d a s a n a u t o r e g r e s s i v e  p r o c e s s o f o r d e r m. prediction  noise  (3.2)  or  (3.3)  where the  )f {t) m  is  following.  the  m point  prediction  error  filter  defined  by  54 (3.4)  By  -°*  =  Yjt)  minimizing the  Treitel  c-t)  M  o  p r e d i c t i o n e r r o r power, Peacock  (1969) have shown that t h i s p r e d i c t i o n e r r o r  s a t i s f i e s the eguations d e f i n i n g filter  t>  )  to within  and  filter  a minimum phase Wiener  a scale factor.  inverse  That i s ,  (3.5)  where P^  i s the  of l e n g t h  m.  presented i n Levinson  The  algorithm  (1.8),  (1947).  e x p l o i t s the  p r e d i c t i o n e r r o r power a s s o c i a t e d  (2.5)  and  (3.5)  T h i s algorithm  was  filter  i s chosen such  of the f i l t e r that:  s o l u t i o n which  autocorrelation  by combining the p r e d i c t i o n e r r o r and m through use  given by  i s a recursive  prediction error f i l t e r  length  filter  f o r s o l v i n g the normal eguations  T o e p l i t z form of the  t h i s scheme, the  with a  matrix.  of l e n g t h  hindsight  coefficient  m+1  In  i s formed  error  filters  o^ (m) .  The  m  of  and  where  Therefore,  °^_(m) i s chosen so that  56 (3.10)  The given  recursive  s o l u t i o n f o r the  prediction f i l t e r  i s then  by  (3.11)  As seen from the reflected  previous d i s c u s s i o n  waves, oC (m)  can  coefficient calculated  recursion  a l s o has  out  filters on  prediction distance out  by Peacock and  digressing  i t should  be  filters  than  one  (1969), such  water bottom  multiples  the s o l u t i o n of the s i n g l e channel Wiener  the q u e s t i o n of determining a s u i t a b l e f i l t e r  considered.  inverse  Before  records.  Having discussed  If  last  Levinson  i s greater  Treitel  be very u s e f u l i n e l i m i n a t i n g  marine seismic  filter,  , be  that an important c l a s s of p r e d i c t i o n e r r o r  As pointed can  by each step i n the  of p r e d i c t i o n e r r o r f i l t e r s ,  e x i s t s f o r which the sample.  That i s , the  a physical interpretation.  from the d i s c u s s i o n pointed  vertically  , under c e r t a i n c o n d i t i o n s  c o n s i d e r e d as a r e f l e c t i o n c o e f f i c i e n t . filter  of  length  is  f  a wavelet i s minimum phase, the z transform of i t s exact filter  i s given  by:  57 (3,12) •t=o  In other words, a c a u s a l f i l t e r  of i n f i n i t e  length i s  r e g u i r e d i n the s i n g l e channel case t o perform an i d e a l deconvolution. filter  Despite the f a c t t h a t performance of the i n v e r s e  should improve as i t s l e n g t h i n c r e a s e s , expense and  roundoff e r r o r s a r e computational length.  Although  f a c t o r s which r e s t r i c t  filter  the choice of a s u i t a b l e l e n g t h f o r the s i n g l e  channel  Wiener f i l t e r  monitor  the normalized  i s s u b j e c t i v e , i t i s i n s t r u c t i v e to mean square  e r r o r , E(m).  This e r r o r  f u n c t i o n equals t h e previous e x p r e s s i o n given i n (2.1) d i v i d e d by a n o r m a l i z i n g f a c t o r , r ^ ( 0 ) , which i s the zero l a g value of 2  the a u t o c o r r e l a t i o n of the d e s i r e d output.  For a f i l t e r of  length m, E(m) reduces to  (3. 13)  The  use o f t h i s e r r o r f u n c t i o n i n designing Wiener  i s d e s c r i b e d by T r e i t e l and Robinson  filters  (1966) and Wiggins  (1968).  Often an a p p r o p r i a t e length f o r the Wiener i n v e r s e f x l t e r somewhere between the wavelet  lies  l e n g t h and twice the wavelet  length. In the design of i n v e r s e f i l t e r s , prewhitening.  Prewhitening  i t i s u s e f u l t o apply  i n v o l v e s adding a constant t o the  58 spectrum  of the wavelet.  spectrum  so t h a t the spectrum's r e c i p r o c a l  inverse f i l t e r ) frequency ranges  T h i s e l i m i n a t e s zeros i n the (the spectrum  i s w e l l d e f i n e d at a l l f r e q u e n c i e s . where the wavelet's  wavelet of the  In  s p e c t r a l values are s m a l l ,  i n v e r s e f i l t e r s enhance n o i s y p o r t i o n s of the t r a c e ' s spectrum. In  Wiener f i l t e r  design, prewhitening can a l s o be  accomplished  by m u l t i p l y i n g the wavelet a u t o c o r r e l a t i o n ' s zero l a g value by some number g r e a t e r than  1.0.  white noise to the wavelet. prewhitening of  This i s equivalent to A demonstration  i s shown i n F i g u r e 1 1 .  wavelet e s t i m a t e s from CDP  Wiener i n v e r s e f i l t e r s wavelets to a s p i k e .  s t a c k s of data, as well as the  with no prewhitening  using 10% prewhitening.  The  filtered  g i v e outputs more h i g h l y r e s o l v e d than  d e c o n v o l u t i o n s obtained from  a time gate.  are  output the  wavelets o b t a i n e d using no those  However, these h i g h l y r e s o l v e d  outputs a l s o c o n t a i n high frequency n o i s e .  prewhitening  The  the c o n v o l u t i o n of the f i l t e r s with  obtained using prewhitening. filter  This figure d i s p l a y s a set  F i l t e r s designed  wavelets i s a l s o shown. prewhitening  of the e f f e c t s of  which have been designed to shape these  compared to those designed which r e s u l t s from  adding  filters  designed  Furthermore,  with no  can d i s p l a y u n d e s i r a b l e edge e f f e c t s at the end of Hence, i n choosing a l e v e l of prewhitening, a  t r a d e o f f between r e s o l u t i o n and  the amount o f high frequency  n o i s e must be c o n s i d e r e d . When a wavelet's  spectrum  f r e q u e n c i e s , bandpass f i l t e r i n g  i s limited  to the  lower  i s o b v i o u s l y an a l t e r n a t i v e  method of s u p p r e s s i n g high frequency  noise.  Figure 11.  The a f f e c t  of Wiener f i l t e r s . from  a c t u a l data.  different  a) The wavelets Inverse f i l t e r s  l e v e l s of prewhitening  prewhitening  i s used  Outputs r e s u l t i n g filters  of prewhitening i n the d e s i g n shown were estimated  designed  with  are d i s p l a y e d .  i n fa) and "\0% prewhitening  from  No i n c) .  the c o n v o l u t i o n s o f the i n v e r s e  with the wavelets are shown.  a)  SOURCE  WRVFLET  b) WIENER  FILTER  1*-+—  c)  FILTERS  OUTPUT  ^ju>-  ESTIMATES  ^V~l~A  j^A-  ^ - i ^ -  WIENER  FILTERS  FILTER  OUTPUT  61 3. 2  M u l t i c h a n n e l Wiener  Thus f a r , only s i n g l e channel  Deconvolution  deconvolution  of s e i s m i c  t r a c e s has been c o n s i d e r e d .  Since modern s e i s m i c data have a  great d e a l of redundancy due  to m u l t i f o l d coverage,  filtering CDP  multichannel  can be used to e x p l o i t the i n f o r m a t i o n contained i n  recordings.  The  i s t h a t the output  b a s i c assumption i n multichannel  s i g n a l s from the f i l t e r i n g  to i n f o r m a t i o n contained i n more than one  filtering  system are  input s i g n a l .  related Hence,  the i t h channel of the f i l t e r e d output i s given by  k (3.14) J  Consequently,  multichannel f i l t e r i n g  involves operations  with matrices i n s t e a d of the s c a l a r o p e r a t i o n s used i n s i n g l e channel f i l t e r i n g .  The  f o l l o w i n g d i s c u s s i o n , which i s very  s i m i l a r to t h a t given by L i n e s of  Wiener m u l t i c h a n n e l  (1974), o u t l i n e s the a p p l i c a t i o n  deconvolution.  The extension of the Levinson algorithm to i t s m u l t i c h a n n e l form was  originally  qiven by S i q g i n s and Robinson  (1965).  S e v e r a l F o r t r a n programs f o r multichannel data a n a l y s i s are given by Robinson  (1967a),  and  an e x c e l l e n t review of  mathematics of multichannel f i l t e r i n g Treitel  (1970).  In t h a t review,  d e f i n i n g the m u l t i c h a n n e l f i l t e r . satisfy  i s presented  Treitel  by  d e r i v e s the  The f i l t e r  the  equations  matrices, f j t )  the f o l l o w i n g system of normal equations  for a f i l t e r  of  62 length  (f(0)  nu  1  f(/)... ffrOj  j  r(0)  r(0  ...  /g(o) 3O)...  aft*);  (3.15)  The  matrix f (t) c o n t a i n s elements f,j  (t) which operate  the j t h input channel to give a c o n t r i b u t i o n channel.  The  t o the i t h output  matrices r ( t ) and g (t) represent the  autocorrelation  and c r o s s c o r r e l a t i o n matrices r e s p e c t i v e l y .  Using the assumption  of e r g o d i c i t y , the values of r ( t ) and g (t)  are given by the f o l l o w i n g r e l a t i o n s .  r(t)  on  1  =  J(t)  5  •  Here (x , (t) , x (t) ,..., 2  sequences on k, channels  and  y3.i7)  S  5"  I  the d e s i r e d output  5  s  s  :  *•  (t) ) r e p r e s e n t s the i n p u t  time  (Y (t) , Y^ (t) ,..., Y^ (t) ) r e p r e s e n t s (  sequences on each of the J[ output  channels.  In the case of multichannel d e c o n v o l u t i o n , k i s u s u a l l y the number of t r a c e s i n a CDP d e s i r e a s i n g l e impulse  andj^  response  i s s e t equal to one s i n c e f o r the redundant d a t a .  The f a c t that r ( t ) i s i n the form of a T o e p l i t z allowed Wiggins and Robinson  we  matrix  (1965) to o b t a i n a r e c u r s i v e  s o l u t i o n to the equations d e f i n i n g a multichannel f i l t e r . flow diagram of t h i s s o l u t i o n i s given i n Figure 12.  A  T h i s flow  diagram c o r r e c t s m i s p r i n t s i n the previous diagram g i v e n i n W i g g i n s and Robinson by Robinson  (1967a).  (1965) by f o l l o w i n g the d e s c r i p t i o n One  of the m i s p r i n t s was  d e s c r i p t i o n of the a l g o r i t h m i n i t i a l i z a t i o n . failed  to mention t h a t the matrix _r (0) was  given  i n the This d e s c r i p t i o n  inverted.  By  using  the m u l t i c h a n n e l s o l u t i o n o u t l i n e d i n Figure 12, an i n v e r s e filter  may  be designed  to convert the source  wavelet  to a s p i k e .  M u l t i c h a n n e l c o n v o l u t i o n of f.(t) with the t r a c e s from a given CDP  then g i v e s an estimate of the impulse  response.  However,  64 the source wavelet  estimates which are the m u l t i c h a n n e l i n p u t s  to t h i s program must not be i d e n t i c a l .  As seen i n F i g u r e 1 2 ,  i n i t i a l i z a t i o n of the r e c u r s i v e process f o r the f i l t e r r e g u i r e s t h a t f(0)=g(0)r_  ( 0 ) , where r (0)  design  i s the i n v e r s e of the  a u t o c o r r e l a t i o n matrix at zero l a g .  I f a l l inputs to the f i l t e r  are i d e n t i c a l then r_(0) i s a matrix  with egual elements and i t s  i n v e r s e does not e x i s t .  T h i s problem of the zero l a g a u t o c o r r -  e l a t i o n matrix being s i n g u l a r i m p l i e s t h a t prewhitening necessary i n the m u l t i c h a n n e l case. accomplished  may be  Prewhitening can e a s i l y be  by m u l t i p l y i n g the d i a g o n a l s of r (0) by a number  s l i g h t l y g r e a t e r than one.  Figure 12,  T h i s flow diagram  the m u l t i c h a n n e l Wiener f i l t e r by Wiggins and Robinson  f o l l o w s the diagram  given  (1965) as w e l l as a l a t e r  d e s c r i p t i o n given by Robinson a^Jz)  of the r e c u r s i v e t o  (1967a).  The polynomials  and b.(z) r e p r e s e n t z transforms of the p r e d i c t i o n  and h i n d s i g h t e r r o r f i l t e r s . coefficients  The s i z e s of the matrix  of these z transforms are determined  number of i n p u t and output channels. used here f o r the f i l t e r  by the  In the convention  c o e f f i c i e n t s , f ( k ) , k denotes m  the time index and the s u b s c r i p t m r e f e r s to the i t e r a t i o n number i n the r e c u r s i v e process.  THE  WIGGINS -  ROBINSON  RECURSIVE  SOLUTION  FOR THE MULTICHANNEL  r (0),r (I), . . . . r(N)  WIENER  FILTER  g ( 0 ) , g (I),... g (N)  INPUTS INITIALIZATION Q  =b  (0)  0  (0)  0  = I  I  m=  0  r (0)  f (0) 0  = r"' ( 0 ) g  (0)  IN  o (Z)  =  m  b ( Z ) = bm m  (0)  + Q  (0)  +b  m  m  (1)  Z +.. . + Q'm (m) Z  (1)  Z + .. . + b ( m ) Z  f  m  (Z) = f  m  a  s  U bm = U  K  a m +1  =  (0)  m  r (m +  +... + o  (m) r (I)  m  (0)  +f  m  (I) + . . . + f_ (m) Z  m  -A  m  m  %m*-t Uam  m  °<m*t  = o t  m K +  o m t  :  f  m ( 0 ) r (m+11 + ^ ( 1 ) r (m)+..>f (m)r(l) m  ,Ubm  T a  (?m + i <?m bnrn-i om  m  :  ^ bm -r i ~ U b m ^ m '  1  u  K  •  -  z  +K  fm +1 = (g (m +i) - 2f m + i) /9 m 11  3. a + i (Z) = Qm (Z) + K  a  +  I  bm+i (Z) = b  b m +  I  m  a  m  +  , (Z), b  m  m  (Z) + K  m  Z">*' b Z  m  +  'a  m  m  (Z (Z  )  fm  +  .(Z)= f ( Z ) + K m  )  + |(Z)  ' m * i (Z) SET  m  EQUAL  TO'm+i,.'  STOP  WHEN  m=N  f  m  t  l  Z<""b  m + l  (Z-)  67 An i n t e r e s t i n g its  f e a t u r e of the multichannel Wiener f i l t e r i s  p o s s i b i l i t y f o r i d e a l performance as pointed out  Treitel  (1974).  By examining  by  the number of eguations  and  unknowns f o r a system of k input channels and a s i n g l e channel, T r e i t e l Wiener f i l t e r  (1974) demonstrates that the m u l t i c h a n n e l  has an exact s o l u t i o n whenever the f i l t e r  a l e n g t h egual to i n p u t sequence.  (N-k)/(k-1)|  A  an exact s o l u t i o n t o be  no  performance e x i s t s for the  channel case s i n c e the f i l t e r  similar single  must become i n f i n i t e i n length f o r  found.  the s y n t h e t i c example shown i n F i g u r e 13(a), the use of  Wiener multichannel d e c o n v o l u t i o n was conditions.  s t u d i e d under i d e a l  In order to make r(0) n o n s i n g u l a r , d i f f e r e n t  r e a l i z a t i o n s of white noise were added t o the wavelet 4 t o produce the wavelets of F i g u r e 13. with an impulse  response  l e n g t h of the f i l t e r  was  the source wavelet.  to produce the s y n t h e t i c t r a c e s shown.  chosen t o be one  l e n g t h , the  l e s s than the l e n g t h  When a c t u a l source wavelets are used i n  the design of the i n v e r s e f i l t e r ,  the c o n v o l u t i o n of the  with the s y n t h e t i c t r a c e y i e l d s a deconvolved d u p l i c a t e s the impulse  of Figure  These were convolved  As r e q u i r e d by the preceeding c r i t e r i o n f o r f i l t e r  of  reaches  , where N i s the l e n g t h of the  As seen from the case where k=1,  c o n d i t i o n of i d e a l f i l t e r  In  output  response.  filter  t r a c e which almost  As shown i n Figure  13(b),  s i n g l e channel d e c o n v o l u t i o n f o l l o w e d by t r a c e s t a c k i n g does not produce i d e a l r e s u l t s even when an i d e a l s p i k i n g chosen and the a c t u a l wavelets are  used.  position i s  Figure  13a).  T h i s s y n t h e t i c example  performance o f the m u l t i c h a n n e l Wiener filter  under i d e a l  synthetic gives  response.  trace  deconvolution  When c o n v o l v e d  traces, a multichannel f i l t e r  a deconvolved  impulse  conditions.  exhibits the  which a l m o s t  with  of length  the  14  d u p l i c a t e s the  SOURCE WAVELETS  IMPULSE  RESPONSE  Y -  i—  1  0.0  SYNTHETIC  i 0.0  1— 0.2  0.2  I 0.4  TIME  1— 0.6  i 0.8  TRACE  I  1 0.6  1 0.2  1 0.4  — S r  —i  1—  0.6  0.4  TIME  TRACES  DECONVOLVED  0.0  A  TIME  1 0.6  1 1.0  1 1.0  1  1.2  1 C.3  70  Figure of l e n g t h  13b).  A sura of s i n g l e deconvolutions  14 designed  spiking position.  using a  using the a c t u a l wavelets and an  filter  optimum  71 Although the r e s u l t s of these i d e a l i z e d  t e s t s of the m u l t i -  channel Wiener f i l t e r are encouraging, the t e s t was under the u n r e a l i s t i c assumption known.  conducted  that the source wavelets were  As shown i n Figure 14(a), the deconvolution d e t e r i o r a t e s  from the i d e a l when wavelet e s t i m a t e s are used xn d e s i g n i n g a multichannel f i l t e r . 4 and  13(a)  T h i s example uses two t r a c e s from F i g u r e s  to demonstrate  this.  The average  wavelet e s t i m a t e ,  obtained by homomorphic d e c o n v o l u t i o n , was used t o produce deconvolved  traces.  Deconvolutions were produced  p o s i t i o n s at the f i r s t  the  using spiking  peak and f i r s t trough of the wavelet with  the l e v e l of prewhitening s e t at 5%.  The obvious g u e s t i c n which  a r i s e s a t t h i s stage i s whether or not multichannel d e c o n v o l u t i o n o f f e r s a s i g n i f i c a n t improvement over of  s i n g l e channel d e c o n v o l u t i o n and s t a c k i n g .  flercado  (1968) pursued  combinations  Davies  t h i s q u e s t i o n f o r r e a l data and  and found  K  that the improvements r e s u l t i n g from  the use of m u l t i c h a n n e l  d e c o n v o l u t i o n were i n s i g n i f i c a n t . When the impulse response i s known, a q u a n t i t a t i v e measure of  the success of a d e c o n v o l u t i o n i s given by the energy  normalized c r o s s c o r r e l a t i o n , c ( t ) .  I f i ( t ) i s the a c t u a l  impulse response and x(t) i s the deconvolved t r a c e , c (t)  is  given by  (3.18) The sums of s i n g l e channel d e c o n v o l u t i o n s f o r the t r a c e s of F i g u r e 14(a) are shown i n F i g u r e 14(b)  for a f i l t e r  of the same  72 l e n g t h and with the same l e v e l of prewhitening a s used f o r the m u l t i c h a n n e l d e c o n v o l u t i o n of F i g u r e 14(a).  The maximum c ( t )  values f o r t h e m u l t i c h a n n e l d e c o n v o l u t i o n s were 0.666 and 0.700 f o r the s p i k i n g p o s i t i o n s of 3 and 6 r e s p e c t i v e l y .  The maximum  c ( t ) values f o r s i n g l e channel deconvolutions with t h e same spiking  p o s i t i o n s were 0.653 and 0.691 r e s p e c t i v e l y .  optimum s p i k i n g p o s i t i o n i s chosen f o r t h e s i n g l e  When an  channel  d e c o n v o l u t i o n of a t r a c e s t a c k , the maximum c (t) value was 0.694.  Figure 14a).  M u l t i c h a n n e l deconvolution  from the use of source wavelet  resulting  e s t i m a t i o n i s shown f o r  the cases where the d e s i r e d outputs are s p i k e s a t d e l a y s of two u n i t s and f i v e  units.  74  DELAY=5  Figure 14 b). deconvolutions  s t a c k s of s i n g l e  designed  channel  using the same wavelet  estimates,  the same s p i k i n g p o s i t i o n s and f i l t e r l e n g t h s and the same prewhitening  l e v e l as used i n the multichannel, case  of F i g u r e 14a). The t h i r d t r a c e shows the r e s u l t i n g deconvolution chosen,  when the optimum s p i k i n g p o s i t i o n was  Delay = 2  Delay = 5  .0  Optimum Spike  T h e r e f o r e , as can  be seen from F i g u r e s 13 and  values of the energy normalized d i f f e r e n c e s between s i n g l e and less significant actual  when wavelet  crosscorrelation,  14 and  the  the  multichannel d e c o n v o l u t i o n become estimates are used i n s t e a d of the  wavelets.  Let  us c o n s i d e r an example which compares the r e s u l t s of  s i n g l e channel  and  multichannel deconvolution on noisy data.  F i g u r e 15 shows 15 s e i s m i c t r a c e s of near v e r t i c a l i n c i d e n c e r e f l e c t i o n s from 5 CDPs of a Gulf Coast Multichannel deconvolutions  seismic l i n e .  and s t a c k s of s i n g l e  channel  deconvolutions are computed f o r these bandpassad data,  although  s t a c k i n g f o l l o w e d by deconvolution i s the l e a s t expensive " m u l t i c h a n n e l " methods, problems can develop wavelets Levin  from stacked data.  when e s t i m a t i n g  as pointed out by Dunkin  and  (1973), normal moveout c o r r e c t i o n s f o l l o w e d by s t a c k i n g  w i l l d i s t o r t the s e i s m i c pulse by enhancing content.  the low  freguency  F i g u r e 16 g i v e s a comparison of the s i n g l e and m u l t i -  channel d e c o n v o l u t i o n approaches. have the same l e n g t h and  The  deconvolution  Figure 16(b)  d i s p l a y s CDP  d e c o n v o l u t i o n s , while Figure 16(c) deconvolutions. prewhitening, position.  filters  s p i k i n g p o s i t i o n s , and are designed  the s e t of t r u n c a t e d minimum phase wavelets 16(a).  of the  on  shown i n F i g u r e  s t a c k s of s i n g l e  channel  d i s p l a y s multichannel  Both approaches used a 1% l e v e l of  the same f i l t e r l e n g t h s , and the same s p i k i n g  Filter  l e n g t h s were s a t at the optimum performance  l e n g t h f o r a system of two  input channels, which i s the s m a l l e s t  number of channels  in this set.  per CDP  The normalized mean  78 square  e r r o r i s always monitored  difficulties too l o n g .  i n order to avoid the  numerical  which can a r i s e i f the multichannel f i l t e r  The  r e s u l t s of the deconvolutions  very s i m i l a r , and deconvolutions  becomes  i n F i g u r e 16 are  the p i c k i n g of r e f l e c t i o n s on the  multichannel  i s no e a s i e r than on the stacks of s i n g l e  channel  deconvolutions. Although  the performance of the multichannel approach i s  b e t t e r than the s i n g l e channel  approach f o r soma s y n t h e t i c  examples, the improvement found  by using multichannel  d e c o n v o l u t i o n on a c t u a l data i s o f t e n e i t h e r nonexistent.  T h i s multichannel Wiener f i l t e r  redundancy by using wavelet channel  marginal  Wiener or  e x p l o i t s data  c r o s s c o r r e l a t i o n s , whereas s i n g l e  methods employ s t a c k i n g .  T h i s accounts f o r the g r e a t e r  expense r e q u i r e d i n using m u l t i c h a n n e l deconvolution and higher b e n e f i t to c o s t r a t i o which i s o f t e n obtained of  s i n g l e channel  filtering.  the  through  use  79  flGC ON TRACES SI  7~.  oi  I.I  T.  i.o  T. I.J  1  r.t  ~i  COP 2 --^^\f\J\Af\^  r •••  *•*  —\A/^  -^y\/-^/\f\J\jJ\f\A^^^  -^^v« f\j\jVv^^  CDP 3  s  COP  —vA^V^\/\/\/W\/Wvy\AAAA^—\/w-xA/ CDP 5  F i g u r e 15. channel  Data used i n a comparison between s i n g l e  and m u l t i c h a n n e l d e c o n v o l u t i o n .  of near v e r t i c a l i n c i d e n c e from 5 CDPs. l e s s than  12 msec.  attenuation e f f e c t s .  Traces show r e f l e c t i o n s NMO c o r r e c t i o n s are  h gain f u n c t i o n has been a p p l i e d t o minimize The s e p a r a t i o n o f the CDPs i s 120 f e e t .  Figure 16 a) T r u n c a t i o n s of minimum  phase  wavelets  estimated from two egual gates of the t r a c e s of Figure channel  15.  F i g u r e 16b) shows CDP s t a c k s of s i n g l e  Wiener deconvolutions using f i l t e r s of l e n g t h 72  msec which were designed to spike the wavelets shown i n a) , while 16c) shows multichannel deconvolutions using the same f i l t e r b) .  Prewhitening  designed  l e n g t h and s p i k i n g p o s i t i o n as i n  was 1% i n both  cases.  81  b )  STACKED DECON  ESTIMATED WAVELETS 0.2  i  1  0.3  0.2  TIME  V-  Y y Y  ^  V V-  C )  I  1 1.2  1 1.4  1 1.6  1  l.l  1 2.0  i 2.2  1.3  1.2  TIME  1  1.6  1  l.l  1  2 3  ,  *-2  w\|yyvvA^^  Y  0.0  1 1.0  MULTICHANNEL I 1 1 "WIENER 1 DECON 1 0.1 1.4  0.6  Y  I  0.1  TIME  /j^  0.2  1  i O.S  1 _  0.2  TINE  I  0.1  1  0.1  1  1.0  1  1.2  1  1.4  TIME  1  l.l  1  l.l  1  2.0  1  2.2  82 The Crump use  Kalman d e c o n v o l u t i o n  (1974),  also  a p p l i e s a type  o f t h e Kalman f i l t e r  results  (1975)).  wavelet  estimates  data  The Kalman  filter  shown i n F i g u r e  The problem time  and t h e a c c u r a c y time  adaptive  faced  the  that  we r e c a l l  Time A d a p t i v e  statistical  l e a d s us t o t h e n e x t  that o f designing a  Deconvolutign  of nonstationarity.  B i c k e r ' s 1953 t r e a t i s e  speaking,  a l l statistical  stationarity  any o f t h e  estimates.  of the wavelet  shape o f a s e i s m i c p u l s e changed  invariant.  the property  we a r e i n e v i t a b l y To a p p r e c i a t e  I f the process  which d e s c r i b e d how  during  propagation.  of s t a t i o n a r i t y  examined  i s assumed t o be  r e q u i r e s that only the f i r s t moments ba t i m e  invariant.  and s e c o n d  This  requires  weak o r wide s e n s e  ( J e n k i n s and W a t t s  Gaussian, order  statistical stationarity,  r e q u i r e s t h a t t h e mean, v a r i a n c e , and a u t o c o v a r i a n c e o f time  this  moments o f a t i m e s e q u e n c e be t i m e  p r o p e r t y , w h i c h i s termed  independent  on t h e  depends on t h e g u a l i t y  designing f i l t e r s f o r seismic data,  Strictly  deconvolution  filter.  with t h e problem  concept,  i t s success  The  produces  designed  16, and l i k e  of inverse f i l t e r s ,  3.3 In  here  Wiener  was a l s o  o f the wavelet  v a r y i n g nature  i n the design  analysed  to multichannel  methods o f i n v e r s e f i l t e r i n g , of  a s d e s c r i b e d by  of multichannel f i l t e r .  on t h e d a t a  which a r e s i m i l a r  (Yedlin  approach,  (1969)).  be  For a zero  mean  process such as a s e i s m i c t r a c e , t e s t s for. weak s t a t i o n a r i t y r e q u i r e examination  of t h e ' a u t o c o r r e l a t i o n f u n c t i o n .  Hence, the  a u t o c o r r e l a t i o n f u n c t i o n i s not only necessary i n the design of a Wiener f i l t e r  but i s a l s o a means of monitoring  nonstationarity.  Since s t a t i o n a r i t y  c o n s t r u c t i o n of a Wiener f i l t e r ,  i s assumed i n the  such f i l t e r s  are designed on  c e r t a i n time qates of the data i n order t o accommodate assumption. in  Although  this  the choice cf qate lengths i s o f t e n done  a s u b j e c t i v e f a s h i o n , s o p h i s t i c a t e d techniques which u t i l i z e  a u t o c o r r e l a t i o n e s t i m a t e s may be used to determine gate l e n g t h  (Wang  a suitable  (1969)).  S t a t i o n a r i t y t e s t s , however, are not the only c r i t e r i o n f o r choosing gate l e n g t h s . refrain traces.  from s p l i t t i n g  In choosing a gate, we should  also  prominent r e f l e c t i o n events i n the  furthermore, although gate l e n g t h s must be s u f f i c i e n t l y  s h o r t t o i n s u r e some semblance of s t a t i o n a r i t y w i t h i n the gate, the gates must a l s o be s u f f i c i e n t l y of  a random impulse In  response  long so that the  i s not badly  assumption  violated.  choosing gate l e n g t h s , i t i s u s e f u l t o monitor the  autocorrelation function.  C o n v e n t i o n a l estimates of the  a u t o c o r r e l a t i o n , r (t) , are given by '  84  r(t)  = 7777  N  ^  x(t-i-z)xc^)  .  (3  19)  0  Ts  (x (0) , x (1) ,.,. ,x (N) ) r e p r e s e n t s a t i m e g a t e o f t h e  where  o r i g i n a l seismic trace.  T h i s v a l u e o f r (t) r e p r e s e n t s t h e  a u t o c o r r e l a t i o n o f a sequence which gate.  i s zero o u t s i d e t h e time  E x t e n s i o n o f a g i v e n t i m e g a t e by t h e a d d i t i o n o f z e r o s  i m p l i e s t h a t t h e a c t u a l a u t o c o r r e l a t i o n o f a time seguence i s multiplied r(t)  by a B a r t l e t t  window.  An a l t e r n a t i v e  estimate of  t o t h a t g i v e n by (3.19) r e p l a c e s t h e n o r m a l i z i n g f a c t o r ,  N+1, b y N+1-t.  However, t h i s e s t i m a t e c a n h a v e t h e u n f o r t u n a t e  property t h a t r(0) i s not always property  which  incorrectly  s e g u e n c e c a n be n e g a t i v e .  the largest  value of r ( t ) ,  a  i m p l i e s t h a t t h e energy  of a time  The i n c o r r e c t a s s u m p t i o n  o f zero  e x t e n s i o n f o r t h e t i m e g a t e may be r e m e d i e d by u s i n g t h e maximum entropy  method  (MEM) p r o p o s e d  by B u r g  (1967).  The t e c h n i q u e h a s  been d e s c r i b e d i n d e t a i l  by Burg  Bishop  a l g o r i t h m h a s been o u t l i n e d  (1975).  (1974).  The B u r g  (1975) and uMrych a n d by A n d e r s e n  B u r g ' s method f o r c a l c u l a t i n g t h e a u t o c o r r e l a t i o n  only the available assumptions  uses  d a t a i n t h e t i m e g a t e and m i n i m i z e s t h e  made a b o u t t h e d a t a b e h a v i o u r  outside the time  gate.  T h i s i s e q u i v a l e n t t o s a y i n g t h a t t h e e s t i m a t e d i s p l a y s maximum entropy.  The MEM c a l c u l a t i o n  estimates the autocorrelation i n a  r e c u r s i v e f a s h i o n by u s i n g t h e p r e d i c t i o n e r r o r f i l t e r . Burg  a l g o r i t h m c a n be u n d e r s t o o d  on p r e d i c t i o n f i l t e r i n g  by r e f e r r i n g  The  to the d i s c u s s i o n  and a u t o r e g r e s s i v e p r o c e s s e s .  For Burg's a l g o r i t h m o n l y v a l u e s o f the p r e d i c t i o n  error  85 seguence between t=m and t=N are considered s i n c e i n t h i s range a ( t ) i s d e f i n e d only i n terms of a v a i l a b l e values of x ( t ) , m  That i s ,  (3.20)  i s considered  f o r values of t=m,m+1,...„ ,'N.  As seen from the above, a ( t ) i s a p o r t i o n of t h e m  deconvolved t r a c e s i n c e filter.  Y {^) i s  a  m  minimum phase  deconvolution  The h i n d s i g h t e r r o r seguence, b ^ ( t ) , i s determined by  c o n v o l v i n g x(t) with the h i n d s i g h t e r r o r f i l t e r , hindsight error f i l t e r  where the  i s d e f i n e d as the time reverse of the  prediction error f i l t e r .  Again,  only values of b ( t ) m  been c a l c u l a t e d from a v a i l a b l e data have been used.  which have That i s ,  the f o l l o w i n g e x p r e s s i o n  is  used f o r t=0,1,2,...,N-m. C o n s i d e r i n g values of a ^ f t ) and b ^ f t ) i n a r e s t r i c t e d range  of  t does not r e q u i r e that the data be appended with  Consequently the p r i n c i p l e of maximum entropy  zeros,  i s not v i o l a t e d .  86 The of  the  Burg a l g o r i t h m  Levinson algorithm  constraint that f i l t e r sum  for calculating  of the  to s o l v e  c o e f f i c i e n t s be  p r e d i c t i o n and  hindsight  c o n s t r a i n t i s necessary since concerning the the  As length  m+1  given  error f i l t e r s  by  In using b (t) m  can  the  chosen to m i n i m i z e  use  e r r o r powers.  The  {Andersen  the  the  extra  a p r i o r i assumptions are  made  (1974)),  autocorrelation i s calculated  L e v i n s o n scheme, t h e  i s f o r m e d by  coefficient are  the  requires  in  In  a  fashion.  i n the  hindsight  no  (t)  normal e q u a t i o n s with  autocorrelation estimates  Burg a l g o r i t h m ,  recursive  the  %  be  °^ {®) • m  the  the  The  prediction error  combining the  filter  prediction error  o f l e n g t h ra t h r o u g h use  of t h e  of  and filter  updated p r e d i c t i o n f i l t e r c o e f f i c i e n t s  recursive  r e l a t i o n s h i p of  above, r e c u r s i o n  (3.11).  r e l a t i o n s between a ^ ( t )  and  deduced.  (3.22)  t=m,  m + 1,. . , ,N.  87 (3.23)  t=0,1 ,.  In the Burg a l g o r i t h m , sum  N-m.  the p r e d i c t i o n f i l t e r minimizes the  of the p r e d i c t i o n and h i n d s i g h t  Although the a u t o c o r r e l a t i o n  e r r o r powers.  This  does not have t o be  explicitly  e v a l u a t e d to determine the Burg p r e d i c t i o n e r r o r f i l t e r , prediction f i l t e r  the  normal e q u a t i o n s given i n (3.6) allows the a u t o c o r r e l a t i o n  Burg  by  (3.19) and  using  the  function  defined  requires  to be r e c u r s i v e l y determined.  (1975) that r (t)  relations.  (3.22) as well as  I t i s demonstrated  can be c a l c u l a t e d by use of the  by  following  r ( o ) = -i_ A l l  ^  i  (3.25)  *  and  =  rim)  r^m-t)  As shown by U l r y c h power s p e c t r a  (3.26)  j  (1972), Burg's method f o r c a l c u l a t i n g  (or i t s time domain e q u i v a l e n t ,  the  a u t o c o r r e l a t i o n ) i s e s p e c i a l l y u s e f u l f o r short data  lengths.  T h i s makes the Burg method i d e a l f o r c a l c u l a t i n g p r e d i c t i o n e r r o r f i l t e r s f o r s h o r t time i n t e r v a l s , which i s a property of a time adaptive  deconvolution  necessary  scheme.  For longer data l e n g t h s , d i f f e r e n c e s between Burg's method and  c o n v e n t i o n a l estimates  become l e s s apparent.  T h i s i s shown  f o r the a u t o c o r r e l a t i o n s computed f o r three gates of the t r a c e shown i n F i g u r e  17.  Moreover, i f only the f i r s t  the t r a c e a u t o c o r r e l a t i o n are r e l i a b l e estimates  few  lobes of  of the  wavelet  a u t o c o r r e l a t i o n because of l a c k of randomness i n the reflectivity  seguence, i t can be argued that the  does not o f f e r s i g n i f i c a n t improvement over  MEM  approach  conventional  approaches f o r the gate lengths used i n e x p l o r a t i o n data processing. F i g u r e 18 compares the use of MEM deconvolution  and  on a noisy s y n t h e t i c t r a c e .  minimum phase Wiener Figure 18(a)  displays  89 a  s y n t h e t i c s e i s m i c t r a c e formed by a d d i n g  convolution  o f a minimum  delay  wavelet  18(b) shows t h a t t h e Wiener a n d MEM  differ  significantly  differences prediction chosen  which  described  by A k a i k e  s e e n from  the values  better  estimate  deconvolution. while is  used i n b o t h  error  o f c ( 0 ) , t h e MEM  response  The c (0)  f o r MEM  t h e c ( 0 ) f o r Wiener d e c o n v o l u t i o n  accurate.  5, has been (1975).  deconvolution  of the impulse  than  18(b) i s  (FPE) c r i t e r i o n .  (1969) and U l r y c h and B i s h o p  n o t u n e x p e c t e d b e c a u s e t h e MEM  more  The l e n g t h o f t h e  parts of Figure  i s d i s c u s s e d i n Chapter  value  will  l e n g t h s , due t o t h e  by use o f t h e f i n a l p r e d i c t i o n  This criterion,  seguence.  deconvolutions  i n the a u t o c o r r e l a t i o n estimates. error f i l t e r  noise to the  with a spike  Figure  f o r s h o r t data  white  As  gives a  t h e Wiener  deconvolution i s 0.268.  i s 0.348  This  result  autocorrelation estimate i s  Figure 17. and  Comparison of a u t o c o r r e l a t i o n s using  zero e x t e n s i o n methods.  Nonstationary s y n t h e t i c  t r a c e i s formed by j o i n i n g c o n v o l u t i o n s o f wavelets the impulse  responses shown.  MEM  Autocorrelations are  compared f o r the t h r e e gates shown.  with  91  WAVELETS  IMPULSE  RESPONSE  Figure  18.  a) Noisy s y n t h e t i c t r a c e used i n the  comparison of Wiener and The  Burg p r e d i c t i o n e r r o r  t r a c e i s formed by convolving  with a spike seguence and  filtering.  a minimum delay  then adding 1 0 1  noise.  wavelet  a)  WAVELET  F i g u r e 18 b ) , c ) .  A comparison between HS.H and  Wiener d e c o n v o l u t i o n s i s shown.  In b) the f i r s t  38  samples of the t r a c e are used and i n c) the e n t i r e t r a c e i s used.  The minimum delay d e c o n v o l u t i o n f i l t e r s were  designed from the a u t o c o r r e l a t i o n estimates shown. Filter  l e n g t h s were chosen by use of the F.PE c r i t e r i o n ,  and the deconvolutions were bandpassed freguency  noise.  to reduce high  96  For the  MEM  Riley  calculating autocorrelations approach i s e s s e n t i a l .  and  Burg  (1972) i n t h e  f o r short  Such an  design of  data  approach  was  design.  the  filter  prediction error  tion error f i l t e r  filter  seguences c a l c u l a t e d approach a p p l i e s hindsight so  that  throughout the from the  i s set e g u a l to the  i s deduced from the  predicted  from w i t h i n the  exponential  e s s e n t i a l l y the a p p l i c a t i o n of filter  Riley-Burg  sometimes s h i f t e d  r e a s o n , the  envelope of the  compared t o t h e arise so  with  small  down and of  the  the  that the  the trace  closest  sliding The error  Riley -  Burg  The  Figure  19  to  the  v a l u e of in  cf the  shows  the  an adaptive The  a p p e a r a n c e but  wrong s i g n .  For  d e c o n v o l u t i o n s i s a l s o shown impulse response.  a p p r o a c h when t h e  assumption of  the this and  Problems  gate l e n g t h  can  becomes  a random i m p u l s e r e s p o n s e b r e a k s  autocorrelation  wavelet a u t o c o r r e l a t i o n .  the  is  t i m e g a t e o f 40 msec.  and  and values  been c a l c u l a t e d  a very spiked  envelope of the Riley-Burg  The  a p p r o a c h w h i c h u s e s an  c a l c u l a t e d over a short  are  a  c a l c u l a t i o n of the  deconvolved trace. the  predic-  hindsight  prediction error f i l t e r ,  r e s u l t a n t d e c o n v o l u t i o n has spikes  window.  of  weighting to prediction e r r o r  p r e d i c t i o n e r r o r s e g u e n c e , w h i c h has  by  updates  f i l t e r length.  w e i g h t e d most h e a v i l y .  continuous updating of the  use  The  data within  prediction errors for points  value are  trace.  p r e d i c t i o n e r r o r and  e r r o r seguences i n the  the  the  T h i s approach c o n t i n u o u s l y  i s calculated  window whose l e n g t h  applied  a time adaptive  d e c o n v o l u t i o n scheme w h i c h t o seme e x t e n t o b v i a t e s time gates i n f i l t e r  lengths,  does not  In that  give  a good  case, the  estimate  prediction-  97 error f i l t e r trace.  i s not  s u i t a b l e f o r removing the  A l s o , to o b t a i n  a deconvolution that  impulse response shown i n Figure filtering  be a p p l i e d  to the  19 r e q u i r e s  approach, the  t e l e s e i s m i c data,  (Davies  w i l l resemble that  wavelets by  method i s a p p r o p r i a t e  short MEM  Although the  (1976) ), s i n c e " g l i t c h e s " are  prediction error  filters.  the  for "deglitchinq"  i s o l a t e d p r e d i c t i o n e r r o r s i n the data which can w e l l by  the  bandpass  p r e d i c t i o n e r r o r seguence.  problems e x i s t i n t r y i n g to remove s e i s m i c Riley-Burg  wavelet from  be p i c k e d  large out  Figure  19.  An a p p l i c a t i o n of Riley-Burg  adaptive d e c o n v o l u t i o n . filter  A continuously  time  time  of l e n g t h 20 msec (or 10 samples) was  adapting a p p l i e d to a  bandpassed i n p u t t r a c e formed by the c o n v o l u t i o n  of a  wavelet with an impulse response deduced from a well log synthetic.  The  envelope of the d e c o n v o l u t i o n ,  p r e d i c t i o n e r r o r seguence, and impulse response are  compared.  or  the envelope of the  SOURCE WAVELET IMPULSE  i  »•»  1  INPUT  BURG  RESPONSE  TRACE  DECON  it i . i.J.A Jill j! , J i J .ll 1  ,  1—  .0  0.2  p^'^j^N", i M 1' 0.4  0.6  A  1  0  i  ll-  'T  i,Lli. 11 i. Ii. .L  I.i L.l 1  •  L  2.2  TIME  "]  IMPULSE RESPONSE ENVELOPE i  0.0  DECON  0.0  r  0.2  T 0.4  0.6  0.8  T 1.0  0.6  0.8  1.0  i  r  — r 1.2  1.4  1.6  1.8  2.0  2.2  1.2  1.4  1.6  l.B  2.0  2.2  TIME  ENVELOPE  0.2  0.4  TIME  VO VO  100 Gutowski  (personal communication) has suggested  an HEM  approach which provides a good compromise between the time a d a p t i v e deconvolution suggested  by R i l e y and Burg  (1972) and  the c o n v e n t i o n a l deconvolution of l a r g e r data gates.  This  modified MEM approach a l s o uses a s l i d i n g r e c t a n g u l a r time in  which the p r e d i c t i o n e r r o r f i l t e r  gate  i s c o n s t a n t l y updated.  U n l i k e the R i l e y - B u r g approach, the f i l t e r i s not the same s i z e as the window, but i s u s u a l l y s m a l l e r .  T h i s modified  time  a d a p t i v e MEM approach was a p p l i e d to a n o n s t a t i o n a r y t r a c e model shown i n F i g u r e 20,  The CVL s y n t h e t i c from 0.4 t o 2.2 sec was  used, as well as a wavelet  which changed shape every  600 msec.  Noise was added t o the c o n v o l u t i o n o f t h e time v a r y i n g with the impulse in  response  .  A window s i z e o f 120 msec was used  the design of a 50 msec p r e d i c t i o n e r r o r o p e r a t o r .  filtering  wavelet  was a p p l i e d to reduce high frequency  Bandpass  noise i n the  deconvolution. It  i s o f t e n u s e f u l t o c a l c u l a t e the envelope  a d a p t i v e MEM d e c o n v o l u t i o n s Bracewell  (1965)).  the envelope  by using the H i l b e r t transform (  Clayton e t a l (1975) have demonstrated t h a t  can be very u s e f u l i n determining t e l e s e i s m i c  a r r i v a l times. deconvolution  of time  F i g u r e 20 compares the envelope with the envelope  of t h e  of the impulse response t o  demonstrate the u s e f u l n e s s o f t h i s c a l c u l a t i o n f o r the d e t e r m i n a t i o n of r e f l e c t e d s e i s m i c a r r i v a l s . envelope  Although the  c a l c u l a t i o n s a c r i f i c e s phase i n f o r m a t i o n , i t provides  complementary i n f o r m a t i o n to the time adaptive d e c o n v o l u t i o n s of reflection  data.  F i g u r e 20. deconvolution.  An a p p l i c a t i o n of time adaptive MEM wavelets of d i f f e r e n t l e n g t h s were  convolved with s e c t i o n s of an impulse response from a CVL.  deduced  The p o s i t i o n of the wavelets i n d i c a t e s  whic  p o r t i o n of the impulse response they were convolved with A time adaptive f i l t e r of l e n g t h 50 msec was a p p l i e d t o the t r a c e and bandpassed  t o produce the d e c o n v o l u t i o n .  The envelopes of the deconvolution and the bandpassed impulse response are a l s o compared.  102  SOURCE  WAVELETS  IMPULSE  INPUT  TRACE  i 0.0  RESPONSE  r 0.2  0.4  lullJUJilk II I! |r|lI P i f  llhkM  i l l f Mil/ flip Iff  1  (My * | vv  DECONVOLUTION  • 0.0  1 0.2  ENVELOPE  r  J  f  3.4  BANDPflSSED  1  IMPULSE  RESPONSE  2.0  DECON  l  0.0  M  2.2  ENVELOPE  l  0.2  r—  0.4  •  1  0.6  1  0.0  —i  "—l  1.0  TIME  1.2  1  1  1.4  1.6  :  1  l.B  r  2 0  2 2  103 CHAPTEE 4  A CASE HISTORY QF DECONVOLVING REAL SEISHIC DATA  H'l Deconvolutions  Performed  S e v e r a l d e c o n v o l u t i o n techniques were a p p l i e d t o s e i s m i c data from the Gulf Coast area.  The d e c o n v o l u t i o n r e s u l t s were  compared to s y n t h e t i c s obtained from i n t e r p r e t q e o l o g i c boundaries.  well loqs i n order t o  The data s e t analysed i s from  f i v e common depth p o i n t s of dynamite generated l a n d data.  Parts  of t h i s data s e t were shown e a r l i e r i n F i g u r e s 6 and 15. In order to f o l l o w the assumptions used of Chapter  i n the t r a c e model  1, most of the p r o c e s s i n g was done on t r a c e s of  s e i s m i c waves r e f l e c t e d at near v e r t i c a l i n c i d e n c e .  F o r these  t r a c e s normal moveout c o r r e c t i o n s are s m a l l and the s t a t i c s c o r r e c t i o n s are 4 msec or l e s s .  In most e x p l o r a t i o n work, l a r g e  amounts of data are used and the assumption  of v e r t i c a l  i n c i d e n c e i s s a c r i f i c e d so that the noise l e v e l i n the t r a c e s can be decreased  by s t a c k i n g .  The purpose o f t h i s study  was to  compare the e f f e c t s of s e v e r a l deconvolution techniques on a l i m i t e d amount of data r a t h e r than t o process a l a r g e amount of data with any g i v e n method. Due t o the n o i s e i n the s i g n a l , the data were bandpass filtered  using an Ormsby f i l t e r  10 and 55 hz.  which passed  f r e g u e n c i e s between  As shown i n F i g u r e 21, the a t t e n u a t i o n l o s s e s i n  the data a l s o had to be compensated f o r .  T h i s was  accomplished  by f i n d i n g the l e a s t sguares f i t c f an e x p o n e n t i a l curve to the  104  absolute amplitude by  o f t h e t r a c e and then  multiplying  the r e c i p r o c a l of the exponential curve.  gain c o n t r o l or s p h e r i c a l divergence  the trace  T h i s i s a form o f  correction.  I t s e f f e c t can  be s e e n by c o m p a r i n g t h e t r a c e s o f F i g u r e 2 1 . Some o f t h e d e c o n v o l u t i o n s f o r t h e s e t o f f i v e common p o i n t s were shown e a r l i e r .  Deconvolutions  performed i n the comparison o f wavelet comparison o f s i n g l e  channel  o f these  data  deconvolution  and m u l t i c h a n n e l  deconvolutions.  F i g u r e 2 2 ( a ) shows a Wiener  w i t h no p r e w h i t e n i n g  W i e n e r minimum p h a s e d e c o n v o l u t i o n  resolution. identified advisable.  data i s  minimum  phase  w h i l e F i g u r e 22(b) g i v e s a w i t h 10% p r e w h i t e n i n g .  f i g u r e s demonstrate that prewhitening deconvolution  were  e s t i m a t o r s and i n t h e  A c o m p a r i s o n o f W i e n e r d e c o n v o l u t i o n s on u n s t a c k e d given i n Figure 22.  depth  with l e s s high freguency  These  allows for a stable noise while  sacrificing  A l s o , s i n c e few r e f l e c t i o n s c a n be e a s i l y on t h e s e c t i o n i n F i g u r e 2 2 , s t a c k i n g o f t h e d a t a i s Before  proceeding  t o compare t h e d e c o n v o l u t i o n s  with  s y n t h e t i c s , one more e x a m p l e o f t h e e f f e c t s o f p r e w h i t e n i n g i s F i g u r e 23 d i s p l a y s d e c o n v o l u t i o n s f o r t h e CDP  shown.  s t a c k s of d a t a from F i g u r e 15. trace of the deconvolutions occur  data  As c a n be s e e n f r o m t h e f i r s t  i n Figure 23(b),  when u s i n g no p r e w h i t e n i n g .  edge e f f e c t s c a n  On t h e o t h e r h a n d , t h e u s e o f  1 0 1 p r e w h i t e n i n g , a s shown i n F i g u r e 2 3 ( c ) , e l i m i n a t e s t h i s problem.  L e t us now e x a m i n e s y n t h e t i c s e i s m o g r a m s w t i i c h  o b t a i n e d from v e l o c i t y  l o g s near t h e s e i s m i c  line.  were  F i g u r e 21a).  Bandpass f i l t e r e d  near the G u l f Coast. F i g u r e 21b)  data from an area  Distance between CDPs i s 120  feet.  shows the data a f t e r s p h e r i c a l divergence  c o r r e c t i o n s were a p p l i e d .  Such c o r r e c t i o n s compensate  f o r the energy l o s s e s i n s p h e r i c a l wavefront  amplitudes,  which are p r o p o r t i o n a l to the r e c i p r o c a l of the d i s t a n c e squared.  These data were used i n subsequent  deconvolution  tests.  106  BANDPASS  a)  F I L T E R E D  TRACES  CDP 1  COP 2 —  ^™*J\f^S\l^^^ COP  3  ^\/\A^AAA/^ V V v ^ v w v ~ n ^ v '  - w - ^ ~ ^ .  COP 4  "  b>AGC  -  n  ON  i  i  O.t  0.9  °/?  -  -  J.  TRACES 1 l  COP 1  1  1  —r I  Jf t.O  1.2  1.6  1.4  I  i  1.0  2.0  1 1.1  A / V -  ^^^/\f\j\^  X  =  X /  -^AJ\I\!\/YA!\^  • ^ -wvyyvwW\/\/\AMA/v\^ •^^\J\J\y-\j\ApJ^^  ^-N/—/v-  — • — • * s ^ J \ J \ J \ f \ J \ [ \ / ^ ^ — v / w w \ , COP :  s^\/W^r\f^^  5  ^  VVvv\/\/\/W^ 0-6  0.8  1.8  1.2  I  TI&  1 1,6  I  '•'  !  1 "  F i g u r e 22 a ) .  Wiener minimum phase d e c o n v o l u t i o n of  the data i n F i g u r e 21a) using no prewhitening filter  l e n g t h o f 80 msec.  and a  F i g u r e 22b) shows a minimum  phase d e c o n v o l u t i o n of the data i n F i g u r e 21a) using a filter  of t h e same l e n g t h with 10% prewhitening.  108  a>  WIENER  1  PHASE  DECON  1  I  1  1  :  1  1  1  0.8  0.8  1.0  1.2  1.4  1.8  l.l  2.0  1 0.9  1 1.0  1 1.2  1 1.4  1 1.6  1 1.0  1 2.0  i  O.I  to  MIN  WIENER  i 0.0  MIN  TIME  PHASE  1  1  1  1  1  1  1.0  1.2  1.4  1.6  1.6  i  1  1  0.8  1.0  1 1.2  1  2.2  DECON  0.8  0.1  2.2  1 -  2.0  1 2.2  1  1  1  1  1  1.4  1.6  1.9  2.0  2.2  TIME  F i g u r e 23 a ) . 21b).  CDP s t a c k s c f the data i n Figure  Wiener d e c o n v o l u t i o n s of s t a c k s with no  prewhitening  are shown i n b ) , while  similar  d e c o n v o l u t i o n s of s t a c k s with 10% prewhitening a r e shown in c ) .  a,  C D P STACKS  i —  0.6  fo  ~1 1.0  1.2  1 1 .*  TIME  -1— 1.6  i'.o  i.e  WIENER DECON ON SOURCE EST  AAAJ  c) WIENER DECON ON SOURCE EST  2.?  111 4.2  Impulse Response Models Obtained From  Synthetic  Seisaojrams Estimates interest  of the earth's  sere obtained  information.  impulse response i n the region  by using two  s e t s of seismic v e l o c i t y  These i n c l u d e d a continuous v e l o c i t y l o g  a v e l o c i t y - d e p t h i n t e r p r e t a t i o n obtained S y n t h e t i c seismograms of the e a r t h ' s were obtained  response to a s p i k e  by using the approach of Berryman et a l  (1969).  s u r f a c e , was  r e f l e c t i o n s from the  and  data. impulse a  (1960).  presented  A v e r s i o n of t h i s program, r e v i s e d by  to i n c l u d e a source at depth and  (CVL)  from check shot  g e n e r a l program which a p p l i e d t h i s approach was Clowes  of  by  C. Chapman free  used to o b t a i n s y n t h e t i c s f o r both sets of v e l o c i t y  information.  Although d e n s i t y data  ware a l s o a v a i l a b l e , the  method of Berryman et a l (1960) uses information v e l o c i t y v a r i a t i o n s i n the l a y e r s only. from t h i s area  are o f t e n u n r e l i a b l e .  p e r t a i n i n g to  A l s o , the d e n s i t y  Furthermore, i n t h i s  logs case,  s y n t h e t i c s which used v e l o c i t y l o g s were very s i m i l a r to those obtained  by using v e l o c i t y and  communication)).  density  This s i m i l a r i t y  (Lucas  occurs  (personal  because the  percentage  change i n v e l o c i t y i s u s u a l l y l a r g e r than the percentage change i n d e n s i t y so t h a t the  value of the r e f l e c t i o n c o e f f i c i e n t i s  l a r g e l y determined by v e l o c i t y v a r i a t i o n s . v a r i a t i o n s i n t h i s case o f t e n c o i n c i d e d with changes.  The  Also, density large v e l o c i t y  l a r g e s t d e n s i t y v a r i a t i o n f o r the well  logs  occured at a depth of 2200 f e e t where the d e n s i t y decreased from 2.08  g i / c c t o 1.48  d e n s i t y change was  gm/cc i n a d i s t a n c e of 20 f e e t .  This large  not c l e a r l y shown i n the seismic s e c t i o n s  but  112 d i d c o i n c i d e with a v e l o c i t y change found i n the check v e l o c i t y data.  I t should  be noted t h a t using only  v a r i a t i o n s to c o n s t r u c t s y n t h e t i c s i s not of gas  presence  d e p o s i t s because the l a r g e d e n s i t y v a r i a t i o n f o r these  A l s o , i t should constructed  divergence  reflections.  be pointed out that the s y n t h e t i c s were  such t h a t a t t e n u a t i o n was  value of Q equal to 1000  was  used.  negligible.  For t h i s , a  In a d d i t i o n , s p h e r i c a l  c o r r e c t i o n s were a p p l i e d to the s e i s m i c t r a c e s i n  21.  The  two  models of i n t e r v a l v e l o c i t y versus depth are shown  i n F i g u r e 24 along The  velocity  v a l i d i n the  " b r i g h t spots" g i v e s r i s e to l a r g e amplitude  Figure  shot  with s y n t h e t i c s deduced from these  f i r s t v e l o c i t y curve  using a borehole  was  obtained by r e c o r d i n g s u r f a c e  geophone placed at v a r i o u s depths.  v e l o c i t y data are termed check shot data s i n c e the curve obtained  models.  by t h i s technique  can  r e s u l t s of i n t e g r a t i n g a continuous  These velocity  be used to check the velocity log.  shots  Figure 24.  CVL  and  i n t e r p r e t e d check shot  velocity  models with s y n t h e t i c seismograms derived from these v e l o c i t y curves.  The f i r s t  s y n t h e t i c of the check shot  v e l o c i t i e s p l a c e s the shot a t the s u r f a c e and c o n t a i n s n ghosts.  While the second check shot s y n t h e t i c c o n t a i n s  m u l t i p l e s f o r the case where the shot i s b u r i e d a t a depth of 100  feet.  The CVL  synthetic contains multiples  f o r a shot a t the same depth.  INTERPRETED VELOCITY LU  0.0  1.0  _J  '.-CVL  2.0  _j  3.0 _ l  DEPTH IN KFT 4.0 1  5.0  L _  VELOCITY  UJ  1.0  2.0  _J  3.0 _J  INTER VEL SYNTHETIC  DEPTH IN KFT 4.0 l  PRIMARIES ONLY  o.o  o.o  o.o  5.0  i  115 The  second v e l o c i t y - d e p t h  model was  continuous v e l o c i t y l o g which was  sampled at one  In order t o o b t a i n a model which was and  g e o l o g i c a l l y r e a l i s t i c , the  obtained  by  foot i n t e r v a l s .  computationally  v e l o c i t y values  Since the logging  s u p p l i e s r e l i a b l e v e l o c i t y values  only below the  initial  part of the CVL  v e l o c i t y values  model was  from check shot  feasible  of the CVL  averaged over 25 f o o t i n t e r v a l s .  the  using a  were  sonde  well  found by using  casing, the  data.  Using the i n t e r p r e t e d v e l o c i t y from check shot  data,  s y n t h e t i c s were c a l c u l a t e d t o give only primary r e f l e c t i o n s i n one  case and  the f i r s t  primaries  case the shot  second case the shot The included  plus m u l t i p l e s i n the second c a s e .  CVL  was  curve was  primaries  was  and  placed  at the s u r f a c e and  buried at 100  The  a flat  s y n t h e t i c obtained  spectrum between 0 and using  a CVL  250  line.  In the  minimum phase deconvolution  examined p r e v i o u s l y , the impulse response was random or white noise sequence. should  Therefore,  be a d e l t a f u n c t i o n at zero delay.  shown i n Figure the s p i k e at t=0  was hz.  s u p p l i e s a model f o r  impulse response of a l a y e r e d e a r t h i n the v i c i n i t y o f seismic  an  region.  d e l t a f u n c t i o n source used i n the s y n t h e t i c  determined by using The  which  l a t t e r synthetic gives  i n d i c a t i o n of the g e o l o g i c complexity i n t h i s The  i n the  feet.  used to c a l c u l a t e a s y n t h e t i c multiples.  In  the  the  approaches  assumed to be a  i t s autocorrelation The  autocorrelation  25 g i v e s some support f o r t h i s assumption dominates the a u t o c o r r e l a t i o n .  e f f e c t of a ghost r e f l e c t o r i s a l s o shown in the  since  However, the  116 autocorrelation. second  A ghost r e f l e c t o r  i s a l s o evident i n the  s y n t h e t i c i n Figure 24 which i n c l u d e s the m u l t i p l e  r e f l e c t i o n s f o r the check shot v e l o c i t y model. r e f l e c t i o n r e s u l t s when energy  from  A ghost  the shot propagates  s u r f a c e and i s r e f l e c t e d down again.  to the  Ghost r e f l e c t i o n s have the  o p p o s i t e sign from the primary r e f l e c t i o n s and are separated i n time frem the primary  r e f l e c t i o n s by an amount equal to twice  the shot t o s u r f a c e t r a v e l time. was  s e t at 100  In t h i s case the shot  f e e t i n the s y n t h e t i c .  As w i l l be seen  depth later,  the qhost r e f l e c t i o n i s not e v i d e n t i n the a u t o c o r r e l a t i o n of the s e i s m i c t r a c e or the a u t o c o r r e l a t i o n s of subsequent deconvolutions.  The  ghost r e f l e c t o r i s present i n the  s y n t h e t i c s because of the o v e r s i m p l i f i e d had  to be used  f o r the f i r s t  1630  v e l o c i t y model which  f e e t of sediment.  T h i s should  be kept i n mind when c o r r e l a t i n g the s y n t h e t i c s with s e i s m i c data or d e c o n v o l u t i o n s . S i n c e the data ware bandpass f i l t e r e d frequency  to suppress  surface, wave i n f o r m a t i o n and high frequency  best t h a t can be expected  response.  " c o n v e n t i o n a l " i s used here s i n c e Chapter  The  term  5 c o n s i d e r s an  unconventional method of extending the spectrum The  n o i s e , the  from c o n v e n t i o n a l d e c o n v o l u t i o n i s a  bandpassed v e r s i o n of the impulse  response.  low  of the  impulse  bandpassed s y n t h e t i c and i t s a u t o c o r r e l a t i o n are  a l s o shown i n F i g u r e 25.  As expected, bandpass f i l t e r i n g  has  the u n d e s i r a b l e e f f e c t of decreasing the randomness of the impulse  response  while having the d e s i r a b l e e f f e c t of decreasing  the noise contained i n the s e i s m i c t r a c e .  117  Figure  25,  Comparison of the a u t o c o r r e l a t i o n s of  the CVL s y n t h e t i c before Bandpass f i l t e r i n g Ormsby f i l t e r of t h i s f i l t e r 45, and 55 hz.  and a f t e r bandpass  was accomplished by using  as used on the s e i s m i c data.  filtering. the same The  corners  i n the frequency domain were at 10, 15, Note t h e . e f f e c t of the "ghost "  r e f l e c t i o n on the a u t o c o r r e l a t i o n .  A U T O C O R R E L A T I O N  0.2*  '  0.4  0.6  0.8  1.0  . . . . "ghost e f f e c t "  BANDPASS  i — 0.0  0.2  1.4  1.6  n 1.6  1.6 "  n—' 1.8  :  2.0  1 2.2  2.0  ~1 2.2  —  i  TIME:  AUTOCORRELATION  0.2V  1.2  D.4  AFTER  0.6  BANDPASS  o.  i  ' ° TIME  1  1.4 ,  4  J  6  OF SYNTHETIC.  7$  i.2  Co  119  j4.3  R e l a t i n g Be c o n v o l u t i o n s and S y n t h e t i c  Sei.siograms  An i n t e r p r e t a t i o n of the t r a v e l times f o r s e i s m i c r e f l e c t i o n s from v a r i o u s h o r i z o n s i s obtained by using the s y n t h e t i c seismograms and d e c o n v o l u t i o n s .  First  s y n t h e t i c obtained from the CVL i s c o r r e l a t e d coherent  s e t of deconvolved  traces.  Secondly,  of a l l , the  with t h e most t h e best  r e f l e c t o r s from a l l the deconvolutions and s y n t h e t i c s a r e picked and  compared.  The most coherent  s e t of s e i s m i c t r a c e s was  obtained by s t a c k i n g t r a c e s with the same shot t o r e c e i v e r offset distance.  These common o f f s e t t r a c e gathers a r e  mentioned by Wood and T r e i t e l  (1S75).  The s t a c k i n g o f t r a c e s  with a common o f f s e t d i s t a n c e i s v a l i d over short d i s t a n c e s i n which the l a y e r s can be c o n s i d e r e d to be h o r i z o n t a l . s t a c k s were deconvolved designed  on wavelets  by using Wiener f i l t e r s which were  obtained by the H i l b e r t transform  Successive deconvolutions i n c r e a s e the s i m i l a r i t y impulse  response.  estimated  The o f f s e t  method.  were performed on the data i n order to  between the deconvolution and the  With each d e c o n v o l u t i o n , a new wavelet i s  and removed.  This approach, which i s s i m i l a r t o the  i t e r a t i v e d e c o n v o l u t i o n d e s c r i b e d by Robinson a method of compressing  (1975),  represents  s e i s m i c pulses i n an i t e r a t i v e f a s h i o n .  F i g u r e s 26(a) and 26(b) d i s p l a y the f i r s t and second deconvolutions the wavelet  of the o f f s e t s t a c k s .  The f i g u r e s a l s o show t h a t  estimate i s compressed with each d e c o n v o l u t i o n .  As  shown by the a u t o c o r r e l a t i o n s of F i g u r e 27, the randomness and high frequency  content of the t r a c e s i n c r e a s e s with each step i n  the i t e r a t i v e deconvolution  process.  120  Figure 26.  I t e r a t i v e deconvolution  s t a c k s f o r near v e r t i c a l i n c i d e n c e .  The  compressed with each i t e r a t i v e s t e p , and frequency  content  of the  offset  wavelet s i z e i s the  of the t r a c e s i s i n c r e a s e d .  high  121  i 0.5  1  1 1.0  0.8  OFFSET STACKS  O  1 1.1  ft  1  1 1.6  1.4  Tlie  .  «  •  .  !  1 1.8  1 2.0  1 2.2  •  WIENER DECONVOLUTIONS  ^UTQCQRR V  TIME  rir°  0.2  0.0  0.2  TIME  rtr  i.«  BANDPRSS OF IMPULSE RESPONSE  1.0  -1—  1.2  — i  1.4  TIME  '  r  I. 6  1.8  —i 2.0  1 2.2  0.2  TIME  WAVELETS  0.2 1  Jo  o'.S  1  0,V  O^tt  I—  1  * TIME ttJTQCORR  WIENER DECONVOLUTIONS  0.6  0  WAVELETS  0.2  •!  F i g u r e 27.  Comparison of a u t o c o r r e l a t i o n s at  v a r i o u s stages i n the i t e r a t i v e deconvolution  process.  I t can be seen t h a t the d e c o n v o l u t i o n s i n c r e a s e the high frequency content and  randomness.  No "ghost"  reflector  i s e v i d e n t i n the a u t o c o r r e l a t i o n s of the d e c o n v o l u t i o n s .  TRACE 3 ty  AUTOCORRELATIONS 0.2  t  0.4  O.E  0.6  f  1.0  1.2  1.4  TIME  1 W °VH ^ l||  V  W  v  W!  YT  a  '  O  2  0.2  W  - 1 ' rt A  °-  rt  V  ^  r-  0.4  ^ ~  0-6  4  Oi  U -  °V*  - |  ~ „ ~ . r - .  1.0  ">  r  * - - ^ ^ r i ^  1.2  1.4  TiME  0.6  0.8  TIME  JECQN  <-»  l.o  1 2  1.4  l.O  1.2  1.4  AUTOCORRELATIONS °-  O'-B  2  0.6  TIME  ^•_.'~»L .' v yV,r»,—p -  T  >-P  0.2  0.4  0.6  0.2  0.4  0.6  0.8  1.0  0.8  l.o  TIME  TIME  ^ECOND DECON  k  1.2  1.4  1.2  1.4  AUTOCORRELATIONS ?  -  1  r-  1  TI°M E 8  r|7yV\ft\/vH^AV vA- • " i ^ / ^ ^ A ^ ^ v " ^ /Vl  v  H>  r  0.2  0.4  0.6  0.2  0.4  0.6  IMPULSE RESPONSE  1  w  i""-"' — ^ r -  OB  1.0  0.6  1.0  TIME  TIME  •  AUTOCORRELATION  T r  .  •  1  2  1 2  ~i— 1.4  124 As c a n ba seen seismic energy  from  F i g u r e 26, t h e agreement between t h e  d a t a a n d t h a s y n t h e t i c s i s n o t good. normalized  be s h i f t e d results.  crosscorrelation  t o the l e f t Even  by 20 msec t o a g r e e  s o , t h e maximum v a l u e  crosscorrelation Although  suggests  for a shift  deconvolutions  note  that similar  were i n d e p e n d e n t l y  were poor i n t h i s  Amoco  obtained  t h e same a r r i v a l  group e x p e r i e n c e d  synthetics  with  with  seismic  the s e i s m i c  of the energy  0.15.  c a s e , i t was e n c o u r a g i n g  by Amoco B e s e a c h .  times  normalized  between t h e s y n t h e t i c s and  seismic deconvclution results  which a r e shown i n F i g u r e 28 d i s p l a y at almost  t h a t the s y n t h e t i c  o f 20 msec was o n l y  the c r o s s c o r r e l a t i o n s  the  In f a c t , t h e  f o r these T h e Amoco  the prominent  as t h o s e  deconvolutions.  data results  reflections  i n F i g u r e 26.  t h e same d i f f i c u l t y  to  The  in correlating  125  F i g u r e 28.  Deconvolution  results  o b t a i n e d by Amoco  Research f o r the data set analysed i n t h i s  chapter.  126 There are s e v e r a l reasons why impulse response and  deconvolutions  may not be compatible.  First  synthetics  of a l l , the  w e l l ' s p o s i t i o n i s not at the same p o s i t i o n as the s e i s m i c In the case examined here, the s e i s m i c information  line.  data and the w e l l l o g  were gathered at a d i s t a n c e o f about 800 f e e t  apart.  Secondly the nature of a CVL i s such that no v e l o c i t y information well casing.  i s a v a i l a b l e f o r sediments above the bottom of the Hence, r e f l e c t i o n s and m u l t i p l e s r e s u l t i n g from  near surface v e l o c i t y d i s c o n t i n u i t i e s w i l l not be i n c l u d e d i n the s y n t h e t i c , although such a r r i v a l s can be very prominent on an a c t u a l t r a c e .  A l s o , w e l l l o g measurements are taken i n the  g e o l o g i c a l l y d i s t u r b e d environment of a borehole, and the problem of mud seepage can cause the v e l o c i t y i n f o r m a t i o n inaccurate.  Furthermore, the a c t u a l seismic  to be  t r a c e w i l l not  c o n s i s t s t r i c t l y of a P wave r e f l e c t e d a t v e r t i c a l i n c i d e n c e by horizontal layers.  The t r a c e w i l l i n c l u d e s u r f a c e wave energy,  S wave a r r i v a l s , as w e l l as r e f r a c t e d , d i f f r a c t e d and converted arrivals. At t h i s p o i n t , an i n t e r p r e t a t i o n i s obtained deconvolution  and s y n t h e t i c s e i s i o g r a m  by merging the  information.  Figure 29  summarizes the most probable r e f l e c t o r s which have been from the previous  figures.  picked  Figure 29.  Summary of the best r e f l e c t i o n picks  from the d e c o n v o l u t i o n s and s y n t h e t i c seismograms presented. the  picks.  The d i f f e r e n t  arrows denote r e l i a b i l i t y  of  R E F L E C T I O N  P I C K S  e .  F R O M  o8  D E C O N V O L U T I O N S  1.2  1.O  A N D  IA5  1.<*  S Y N T H E T I C  1.3  S E I S M O G R A M S  2O  22 TIME  SOURCE OF  INFORMATION  DECON  OF  (ofrer  DECON DECON  OF OF  DECON  O F FIG. 2 6  DECON  O F FIG. 2 8  CHECK SHOT SYNTHETIC ( F I G . 2<*) SYNTHETIC (FIG. 2 6 )  • 1 i  t 1  t  t  t t tt  I"  1 1  t  .1  t  ( 3 0 5 0 Ft.)  ( 4 I S O Ft.)  ( 5 3 7 0 Ft.)  t  1  t  • 1 1  l  i  jI  4 i I  t t t • 1  •  t J  t  • 1 •  FIG. 2 2  O F FIG 2 3  t  1  FIG. 16  DECON  CVL  i1  1  FIG. IO  N M O )  t t  *  1 1  1  1  1  11 :  t! i  t ( 6 0 0 0 Ff.)  Q U E S T I O N A B L E  I  ( 8 3 0 0 Ft.)  ho 00  129 From the well  log synthetics  and t h e seismic  deconvolutions, c e r t a i n r e f l e c t o r s are consistently interval  between 0.600 and 2.200 seconds.  shown i n the  These appear a t times  of 0.92, 1.20, 1.48, 1.63, and 2.12 seconds.  For the time  used, the depth range of r e f l e c t o r s was approximately 1900 and 8700 f e e t .  gate  between  Osing the average v e l o c i t i e s deduced from  check shot data, a l l o w s us t o d e t e r i n e the depth of t h e reflecting interfaces. of approximately  These i n t e r f a c e s  were s i t u a t e d a t depths  3050, 4190, 5370, 6000 and 8300 f e e t .  wall l o g s y n t h e t i c s  Both  a l s o i n d i c a t e an i n t e r f a c e a t a depth of  6700 f e e t but t h i s i s not evident i n the seismic d a t a . Although  a p l a u s i b l e i n t e r p r e t a t i o n can be obtained  from  the i n f o r m a t i o n a v a i l a b l e , i t i s c l e a r that no s i n g l e deconvolution noisy  method w i l l provide miraculous  r e s u l t s on such  data. Having c o n s i d e r e d deconvolutions on dynamite  data, a t t e n t i o n sources.  i s now focused on data generated  generated by v i b r a t o r  130 5  CH1FTJB  A JEW  5*1  The Clayton  APPROACH TO VIBROSEIS  Isiioduction  following  DECONVOLUTION  to the Vibroseis  discussion,  Method  t a k e n from L i n e s and  (1976), o u t l i n e s a new a p p r o a c h f o r i m p r o v i n g t h e  resolution  of vibroseis  records.  Certain properties  v i b r o s e i s s i g n a l make i t p a r t i c u l a r l y  of the  suitable for using a  d e c o n v o l u t i o n t e c h n i q u e w h i c h c o m b i n e s homomorphic d e c o n v o l u t i o n and  autoregressive prediction.  been a p p l i e d 'Wiggins  In  The t e c h n i q u e h a s a l s o  t o o t h e r d e c o n v o l u t i o n p r o b l e m s by C l a y t o n and  (1976) a n d C l a y t o n  thevibroseis  and U l r y c h  (1976) .  approach, a freguency  generated using a mechanical v i b r a t o r . by  recently  dynamite e x p l o s i o n s ,  modulated s i q n a l i s  Unlike s i q n a l s  produced  t h e frequency range o f t h e v i b r o s e i s  s i g n a l c a n be c o n t r o l l e d .  Since t h e amplitude o f a v i b r o s e i s  sweep i s s u b s t a n t i a l l y l e s s t h a n t h e a m p l i t u d e o f a d y n a m i t e generated pulse,  t h e time d u r a t i o n  order that t h e enerqy associated compatible with explosions.  vibroseis  with  t h e energy o f t r a c e s  thevibroseis  t r a c e be  generated by dynamite  T h i s c o n c e p t was i m p l e m e n t e d by K l a u d e r e t a l .  (1960) i n t h e r e s o l u t i o n the  o f thesiqnal i s increased i n  of chirp radar signals.  signal requires  s i g n a l be c r o s s c o r r e l a t e d  with  that  t h einput  Analysis of  frequency  modulated  the recorded s i g n a l i n order t o  "compress" t h e energy i n t h e v i b r o s e i s t r a c e appearance i s s i m i l a r t o a t r a c e  so that  g e n e r a t e d by a h i g h  its energy  131 dynamite  source.  5.2  A V i b r o s e i s Trace  A noiseless uncorrelated considered input  as the c o n v o l u t i o n  Sodel  vibroseis trace, y (i),  c a n be  of the frequency modulated  sweep p ( t ) , w i t h t h e i m p u l s e r e s p o n s e o f a l a y e r e d  (FM) earth  i(t) .  y(t)  p(t)  =  *  c(t)  The s y m b o l * i s used t h r o u g h o u t  (5.D  to denote  convolution.  The f r e q u e n c y m o d u l a t e d sweep p ( t ) , c a n be r e p r e s e n t e d  p(t)  =  The i n s t a n t a n e o u s JTT -J— , where  KJ  A  sin  2n-F ( t t  x s t h e phase.  |t*)  0  frequency i s qiven  by f ( 1 + k t ) , 0  ( 5 - 2 )  which i s  I f T xs t h e d u r a t x o n o f t h e  sweep, t h e f r e q u e n c y o f t h e i n p u t FM sweep r a n g e s b e t w e e n f f  0  Q  and  (1+kT) . C r o s s c o r r e l a t i o n o f y (t) w i t h t h e i n p u t  i d e a l i z e d model o f a v i b r o s e i s t r a c e , x ( t ) . may  as  be r e p r e s e n t e d  FM sweep y i e l d s  an  The f u n c t i o n x ( t )  i n the f o l l o w i n g c o n v o l u t i o n a l form.  132  x(t)  =  X(i)  =  p(-t)  *  p U )  xid)  (5.3)  or  In and  Wit) $ (Ct)  (5.4)  (5.4), w (t) i s the a u t o c o r r e l a t i o n o f the i n p u t FM sweep  represents  a zero  phase sequence termed the Klauder wavelet,  A model of such an i d e a l i z e d t r a c e i s shown i n Figure  30.  This  f i g u r e shows a v i b r a t o r sweep, p (t) , i n which f =14 hz, T=1.0 0  sec, and k=3. one  Although the sweep i s g e n e r a l l y much l o n g e r  than  second, t h i s example i l l u s t r a t e s the p r o p e r t i e s o f the t r a c e  model.  Figure 30 shows the Klauder wavelet or a u t o c o r r e l a t i o n  of p ( t ) , and an impulse response of random s p i k e s .  The  recorded  v i b r o s e i s t r a c e , y (t) , and the c r o s s c o r r e l a t e d t r a c e , x ( t ) , are a l s o given.  The power spectrum o f x(t) mainly shows the  i n f l u e n c e o f the impulse response, i ( t ) .  Figure 30. vibroseis trace. the  The f i g u r e shows t h e input FM sweep,  sweep's power spectrum, the Klauder wavelet or  autocorrelation the  Example of an i d e a l i z e d model f o r a  of the sweep, an impulse response model,  recorded v i b r o s e i s t r a c e , the c r o s s c o r r e l a t i o n of the  sweep with the recorded t r a c e , and t h e power spectrum of the c r o s s c o r r e l a t e d  trace.  134  POWER  NPUT FM  SPECTRUM  j.O 3.0  AUTOCORRELATION i  1  0.6  1  -0.4  IMPULSE  r  -0.2  0.0  TIME  1  0.2  1  0.4  -1  62.5  FREQUENCY  1—  123.0  1  0.6  RESPONSE  1 1  VIBROSEIS  1.0  TIME TRACE A  0.0  0.2  POWER  SPECTRUM  CROSSCORRELATED,TRACE J.D  0.2  TIM  o.6 V V  b.a  r*-  i.o  0.0  "i 62.5  FREQUENCY  1 125.0  135 As shown by the  f i g u r e , the power spectrum of the  sweep i s r a t h e r f l a t w i t h i n the sweep's freguency 56 h z ) .  The  limits  power spectrum of the i n p u t sweep i s the  spectrum of the Klauder  wavelet.  Therefore,  input (14-  amplitude  deconvolution  t h i s i d e a l i z e d v i b r o s e i s t r a c e w i t h i n the freguency  of  l i m i t s of  the sweep i s f u t i l e , s i n c e a wavelet with a f l a t spectrum  will  produce a t r a c e which has the same amplitude spectrum as that of the bandpassed impulse response. Since any d i s p e r s e d and  s e i s m i c s i g n a l t r a n s m i t t e d through the e a r t h i s attenuated,  a realistic  the input FM s i g n a l to be f i l t e r e d d e s c r i b e s the a t t e n u a t i o n and earth.  by a f u n c t i o n e{t)  which  d i s p e r s i o n p r o p e r t i e s of  the  The f a c t t h a t the input v i b r a t o r s i g n a l i s a l s o not  i d e a l may  be i n c o r p o r a t e d i n t o a ( t ) ,  represents tend  model of y (t) c o n s i d e r s  the f i l t e r i n g  E s s e n t i a l l y e(t)  e f f e c t s of the earth on p(t) and  to decrease the high freguency  content  i n y (t).  In  model, the s i g n a l r e t u r n i n g from a p a r t i c u l a r r e f l e c t o r represented  by p ( t ) * e (t)  the t r a c e i s given  i n s t e a d of p (t)  and  the new  will our  will  be  model f o r  by  (5.5)  Hence, the c r o s s c o r r e l a t e d v i b r o s e i s t r a c e , x ( t ) , i s then given  by  136  x L t )  There In  this  phase  = w(t) *  *  ec-t)  i s evidence to suggest that  respect,  (5.6)  i(-t)  e (t)  i s minimum  t h e p r e v i o u s d i s c u s s i o n s of Futterman's  attenuation  laws  5.3  The V i b r o s e i s J j e c o n v g l u t i g n  of geological interest,  that  we wish t o remove.  (1975), t h i s  wavelet  nor z e r o p h a s e .  Theoretically,  t h e seguence  As p o i n t e d  will  This  Probl-em  i s excited  attenuation  o u t by S i s t o w and J u r c z y k  , i ngeneral,  t h e spectrum  because  be n e i t h e r  limiting  trace.  Various approaches  to  resolve  reduce  approaches within  radar s i g n a l s .  of the impulse response i s  effects  a certain  of the earth.  These  have been p r o p o s e d  to s o l v e  (1960) w e i g h t e d t h e s p e c t r u m Gurbuz  only  f r e g u e n c y band  (1972) a p p l i e d wavelets.  the vibroseis i s dealt  vibroseis this i n order  tapering to  While trace  with.  effects  lobes,  of the crosscorrelated  lobes of vibroseis  are e f f e c t i v e ,  a certain  phase  o f t h e b a n d l i m i t e d s o u r c e and t h e  the r e s o l u t i o n  Klauder et a l .  the s i d e  minimum  w a v e l e t t o have d i s t u r b i n g s i d e  thereby  problem.  i (t) ,  makes t h e d e c o n v o l u t i o n d i f f i c u l t .  and d i s p e r s i o n  cause the v i b r o s e i s  seguence,  w ( t ) * e ( t ) i s t h e wavelet  c o n t i n u o u s up t o t h e N y g u i s t f r e g u e n c y , b u t o n l y bandwidth  minimum  are o f i n t e r e s t .  Since t h e impulse response or r e f l e c t i v i t y is  phase.  these spectrum  The p r o p o s e d  137  deconvolution trace the  method  deals  not only  with  w i t h i n t h e b a n d p a s s on t h e i n p u t  The d e c o n v o l u t i o n  o f the  signal, but also  optimum a s s u m p t i o n a b o u t t h e s p e c t r u m  band.  the spectrum  method d e s c r i b e d  outside  makes  this  freguency  h e r e has t h e f o l l o w i n g  objectives, 1,  To remove t h e e f f e c t s the  v i b r o s e i s spectrum  of t h e impulse the  2.  of a t t e n u a t i o n  FM  i n order  To i n c r e a s e  r e s p o n s e ' s spectrum  the F o u r i e r transform  5.4  employed.  Using  the f i r s t  of the impulse  Method-  objective, cepstral analysis i s  s o , we assume t h a t t h e s p e c t r u m  t o t h e low g u e f r e n c y  the previous  vibroseis  w i t h i n t h e band o f  o f t h e FM sweep.  The D e c p n v p l u t i o n  with  In doing  contributes  estimate  t h e r e s o l u t i o n o f t h e v i b r o s e i s t r a c e by  r e s p o n s e beyond t h e l i m i t s  dealing  t o o b t a i n an  sweep.  extending  In  and d i s p e r s i o n on  o f e (t)  p a r t o f the spectrum  of x ( t ) .  model f o r x (t) i n ( 5 . 6 ) , and d e n o t i n g  w a v e l e t a s v ( t ) = w ( t ) * e ( t ) , we c a n w r i t e  the  log(X(ur))  138  ioa X(ur) = loa  where V(ar) i(t)  V(u)  +  loo  Ifo)  (5.7)  anil I (to-) denote the F o u r i e r transforms of v(f) and  respectively. Hence, as shown b e f o r e , the cepstrum  of the t r a c e , x (t) ,  can be r e w r i t t e n as:  (5.8)  where v ft) i s the cepstrum the cepstrum  of the v i b r o s e i s wavelet  of the impulse  In using the low wavelet, l o g (X (<*r) )  and i (t) i s  response.  g u e f r e n c i e s of x(t) to estimate  the  i s e s s e n t i a l l y c o n s i d e r e d as a s i g n a l  and  the low "freguency" content of t h i s s i g n a l i s used t o produce an estimate of l o g ( V ( u r ) ) .  The  low guefrency p o r t i o n of l o g (X ( o r ) )  can be e x t r a c t e d by using x(t) values f o r | t|<T , where T 0  the " c u t o f f  is  Q  quefrency".  As shown by C l a y t o n and  Wiggins  (1976), computation  of the  phase curve presents a major problem when computing the complex cepstrum.  Since the Klauder  wavelet  i s zero phase and s i n c e  phase unwrapping of noisy data i s u n s t a b l e , the wavelet  i s assumed to be zero phase.  method to t h i s phase assumption  vibroseis  The s e n s i t i v i t y of the  i s tested later.  I t should  noted t h a t Wood (1976) has a l s o used the zero phase  be  assumption  i n a p p l y i n g Wiener deconvolution f i l t e r s to symmetrical  pulse  shapes. In the f i r s t amplitude  step of the deconvolution process,  spectrum i s f l a t t e n e d  the  by c e p s t r a l f i l t e r i n g  of the  zero phase cepstrum while the phase i s l e f t untouched.  The  zero  phase cepstrum of the t r a c e i s given by 2TI  X  j  (t) =  0  ^  /©j  T  I X Cur) I e'  dur  (5>9)  T h i s cepstrum can a l s o be w r i t t e n as  Xo  Ct)  =  ~±  cfur  The f i r s t term i n (5.10) r e p r e s e n t s the l o g  amplitude  spectrum, or zero phase cepstrum, of the wavelet. cepstrum of the wavelet guefrency  be estimated  Since  g e n e r a l l y c o n t r i b u t e s to the  content of the cepstrum, the f i r s t by a p p l y i n g a low cut l i f t e r  JI(u>)| can be estimated  that  the  low  term i n  (5.10) can  to the cepstrum.  by a p p l y i n g the F o u r i e r transform  and e x p o n e n t i a l o p e r a t i o n s to the high cut p o r t i o n of cepstrum.  (5.10)  A second approach f o r e s t i m a t i n g l l f o r H  the  requires  |V (or) | be estimated by a p p l y i n g the F o u r i e r transform  e x p o n e n t i a l o p e r a t i o n s to the low cut p o r t i o n of the Subsequent d i v i s i o n of Although  these two  \X(u)\  by  | V (w) | then y i e l d s  approaches are mathematically  and  cepstrum. j I (<j) ].  e q u i v a l e n t and  140  g i v e s i m i l a r r e s u l t s f o r s y n t h e t i c s , the l a t t e r approach  was  used when a n a l y z i n g r e a l data s i n c e the problem of h i g h guefrency noise i s avoided.  A l s o , one d e a l s d i r e c t l y  with the  amplitude spectrum of the v i b r o s e i s wavelet, which i s e a s i e r t o d e a l with than the amplitude spectrum o f the impulse response. Assuming the phase of the t r a c e , &iUf),  • - fo be the same  as the phase of the impulse response allows the impulse response to be computed by the f o l l o w i n g i n v e r s e F o u r i e r t r a n s f o r m .  Ut)  The  -  Tw  JKw)/  e  e  dur  (  5  -  1  1  )  impulse response which r e s u l t s from t h i s t y p e o f  c e p s t r a l f i l t e r i n g g e n e r a l l y has a spectrum  which i s a f l a t t e n e d  v e r s i o n of the spectrum o f the o r i g i n a l t r a c e .  Nevertheless,  the f a c t that the spectrum i s b a n d l i m i t e d causes pulses i n i ( t )  t o have d i s t u r b i n g s i d e l o b e s .  second o b j e c t i v e of the proposed  reflected  T h i s l e a d s t o the  d e c o n v o l u t i o n method, which i s  t h a t o f extending the F o u r i e r t r a n s f o r m o f the impulse response. C o n v e n t i o n a l data p r o c e s s i n g would l e a v e zeros i n the spectrum sweep.  f o r f r e g u e n c i e s o u t s i d e t h e range of tne i n p u t FM However, such a s p e c t r a l model i s not i n accordance  with  the continuous nature of the impulse response spectrum, nor i s it  i n agreement with the p r i n c i p l e of maximum entropy.  Furthermore, the Paley-fliener c o n d i t i o n r e g u i r e s that values of s t o c h a s t i c processes be nonzero throughout  spectral the e n t i r e  141 frequency range (Smylie e t a l . (1973)). Let us b r i e f l y  review the maximum entropy p r i n c i p l e .  This  p r i n c i p l e r e g u i r e s t h a t the entropy or i n f o r m a t i o n c o n t e n t of a process be maximized when computing a power spectrum c o n s t r a i n t t h a t the spectrum The  under the  be c o n s i s t e n t with a v a i l a b l e data.  maximum entropy approach can be r e l a t e d to a u t o r e g r e s s i v e  m o d e l l i n g , and  i t r e p r e s e n t s the optimum approach f o r  c a l c u l a t i n g power s p e c t r a of f i n i t e data segments. In order to p l a c e nonzero values i n the impulse spectrum used  between zero and  the Nyquist frequency,  response  an approach i s  which i s c o n s i s t e n t with the p r i n c i p l e of maximum entropy.  T h i s approach,  which has been s u c c e s s f u l l y a p p l i e d to  t e l e s e i s m i c data by Clayton and Wiggins of the impulse  response's  (1976), uses estimates  F o u r i e r transform w i t h i n the  frequency  range of the FM sweep i n order t o p r e d i c t F o u r i e r transform values o u t s i d e the range of the sweep. i s determined  The p r e d i c t i o n  operator  by use of the complex Burg a l g o r i t h m , as d e s c r i b e d  by Smylie et a l .  (1973) and Burg  (1975).  Use of maximum  entropy p r e d i c t i o n provides a F o u r i e r transform which i s c o n s i s t e n t with known transform values w i t h i n the range of the i n p u t sweep while.remaining with r e s p e c t t o unknown v a l u e s .  maximally  Ulrych and Bishop  frequency noncommital (1975)  e s t a b l i s h t h i s property of the maximum entropy method a l s o p o i n t out that use of MEM  (MEM)  i n designing a p r e d i c t i o n  and  error  operator i s e q u i v a l e n t to using an a u t o r e g r e s s i v e model f o r the p r e d i c t e d sequence. to f i l l  Use of a u t o r e g r e s s i v e p r e d i c t i o n a l l o w s us  i n unknown p a r t s of the spectrum  by use of known  142 i n f o r m a t i o n i n s t e a d of using the i n c o r r e c t assumption of l e a v i n g zeros i n the spectrum.  The  use of a u t o r e g r e s s i v e p r e d i c t i o n to  r e s o l v e a s i n g l e a r r i v a l i s shown i n F i g u r e In d e c o n v o l u t i o n , left  one  with a d e l t a f u n c t i o n .  wishes to remove t h i s wavelet and  the Nyguist  frequency.  and  deconvolutions  F o u r i e r transform FM  spectrum between 0 hz  F i g u r e s 31 (b) and  31 (c) show s p e c t r a  r e s u l t i n g from d i f f e r e n t p r e d i c t i o n s of the of i ( t ) .  P r e d i c t i o n using the l i m i t s of the  p r e d i c t i o n using the transform  The  flat  sweep g i v e s s l i g h t l y b e t t e r r e s o l u t i o n of the a r r i v a l  yields a result  be  In other words, i d e a l s p e c t r a l  p r e d i c t i o n would g i v e a completely and  31.  values between 20 and  50  while hz  which i s e s s e n t i a l l y the desired d e l t a f u n c t i o n .  l e n g t h of the a u t o r e g r e s s i v e  p r e d i c t i o n f i l t e r had  a length  equal t o h a l f the l e n g t h of the known transform. In our a u t o r e q r a s s i v e modal, the impulse response of e a r t h i s represented  as a d i s c r e t e sum  the  of d e l t a f u n c t i o n s .  (5. 12)  where  ^  i s the sampling i n t e r v a l i n the time domain.  143  n 125.0  62.5  i 187.5  i 250.0  0.6  0.B  FREQUENCY  0.2  0.0  Figure  31a)  Example o f a K l a u d e r  spectrum.  The  wavelet  generating  freguencies  wavelet  i s the a u t o c o r r e l a t i o n i n the  14-56  hz  range.  and  1.0  i t s amplitude  o f an FM  sweep  144  Figure  31b).  autoregressive  The spectrum a f t e r a p p l y i n g an  p r e d i c t o r designed from the F o u r i e r  t r a n s f o r m i n the 14-56 hz range.  F i g u r e 31c) shows the  spectrum and wavelet r e s u l t i n g from the use of the p r e d i c t o r designed from within the 20-50 hz range.  to  SPECTRUM  OF DECON  JINIII 0.0  62.5  125~D  FREQUENCY  187.5  2S0.0  DECONVOLUTION  0.0  o  0.0  0.2  0.4  T  SPECTRUM  1  62.5  1  125.0  FREQUENCY  DECONVOLUTION  0.6  0.8  OF DECON  -T 187.5  1 250.0  146 The  F o u r i e r t r a n s f o r m , I (<«>•) , i s then a sum  of complex  sinusoids  (5.13) K  In p r e d i c t i n g values of I ( (*r ) beyond the frequency  limits  of the v i b r o s e i s sweep, the assumption i s made t h a t a sum s i n u s o i d s can be represented  as an a u t o r e g r e s s i v e p r o c e s s .  U l r y c h et a l (1973) and U l r y c h and C l a y t o n which agree with t h i s assumption.  (1976) g i v e examples  A d i s c u s s i o n on the use of  a u t o r e g r e s s i v e models to represent harmonic processes a d d i t i v e noise has  of  with  been presented by U l r y c h and C l a y t o n  (1976).  An example of a s e r i e s of t h r e e d e l t a f u n c t i o n s with a d d i t i v e noise i s given i n F i g u r e 32(a) corresponding  amplitude  spectrum.  1%  along with i t s  F i g u r e 32 (b) shows that  b a n d l i m i t i n g the output of the d e l t a f u n c t i o n s produces s i d e l o b e s on the a r r i v a l s .  As shown i n F i g u r e 3 2 ( c ) , p r e d i c t i o n of  the b a n d l i m i t e d F o u r i e r transform produces a F o u r i e r transform very s i m i l a r to the o r i g i n a l transform of Figure 32(a).  This  e x t e n s i o n of the t r a c e ' s F o u r i e r transform provides an e x c e l l e n t deconvolution.  147  w  AMPLITUDE SPECTRUM  a)  1  f— 0.0  1  1  33.25  62.5  93.75  ~1  FREQUENCY  125.0  INPUT SEQUENCE  F i g u r e 32a). shown.  The i n p u t s p i k e seguence and i t s spectrum are  The s p i k e sequence i s a s e r i e s of three d e l t a f u n c t i o n s  with 1% a d d i t i v e n o i s e . sinusoids.  I t s F o u r i e r transform  i s a sum of noisy  b)  MODIFIED  T 62.5  31.25  0.0  TRACE  SPECTRUM  93.75  "1 125.0  FREQUENCY  MODIFIED  1 F i g u r e 32b).  I  2  TU M E  Spectrum  range produces the modified  TRACE  0.4  which i s l i m i t e d to the t r a c e shown.  14-56  149  ci  MEM P R E D  ~i  o.o  SPEC  1  31.25  62.5 FREGUENCY  ~1  93.75  1  125.0  DECONVOLUTION  T 7 ~  0.4  TIME  F i g u r e 32c).  A p r e d i c t i o n o p e r a t o r o f l e n g t h 10 designed  on the b a n d l i m i t e d F o u r i e r transform and  the r e s u l t i n g  deconvolution.  produces the spectrum shown  150 Let us c o n s i d e r the design of the p r e d i c t i o n o p e r a t o r itself.  If  $ UT  i s the frequency  domain sampling i n t e r v a l ,  an  rath order a u t o r e g r e s s i v e r e p r e s e n t a t i o n of the d i s c r e t e F o u r i e r transform  of the impulse response would have the f o l l o w i n g form.  where I ( n S u r ) transform  r e p r e s e n t s v a l u e s of the d i s c r e t e F o u r i e r  of i (t) , ot(r\$u>)  o p e r a t o r , and  N( £w)  1  S  a  n  i s the complex Burg p r e d i c t i o n  white noise seguence r e p r e s e n t i n g  the  p r e d i c t i o n e r r o r or i n n o v a t i o n of the model. In a p p l y i n g the Burg a l g o r i t h m , (1, - c < { $OJ) ,-c^{2 Sur)  oC (m $(*>)) , and the complex  of i t s time r e v e r s e are convolved yield  the p r e d i c t i o n e r r o r and  respectively. the sum  Values  of  (1976), t h i s  operator,  implying  conjugate  h i n d s i g h t e r r o r sequences  h i n d s i g h t e r r o r powers  algorithm.  Claerbout  filter  with known values of I ( n £ « ) t o  o^(n £«•>-) are determined by  of the p r e d i c t i o n and  a p p l y i n g the Levinson  the p r e d i c t i o n e r r o r  As pointed out  minimizing and  by  produces a minimum phase p r e d i c t i o n e r r o r  t h a t p r e d i c t e d energy w i l l  not i n c r e a s e  o u t s i d e the known passband. I t i s important designed  to note t h a t the p r e d i c t i o n f i l t e r  over t h a t p o r t i o n of the spectrum f o r which the  to n o i s e r a t i o i s l a r g e . w i l l be  As shown by Figure 31,  i n s i d e the l i m i t s of the FM  v i b r o s e i s technique  sweep.  this  is signal  portion  Fortunately,  the  allows us to know g e n e r a l l y which p o r t i o n of  151 the  spectrum c o n t a i n s For  filter  of the  be  (FPE)  by  Akaike  criterion the  sum  process  of the  of  the  reliable  The  predicted  ( Ulrych  and may  placed  hence c o n t r o l s t h e  5.5 us  of  this  illustrates  a s e r i e s of spikes  |V(co^)'J.  different  G a u s s i a n s and ratio  as  (1975).  phase Klauder  This  the  error  not  the  in a  the  by  length  subjective amount  of  spectrum,  and  deconvolution.  deconvolution  model f o r  the  may  d e t e r m i n e s the  s y n t h e t i c examples  application  Jurczyk  length  chosen  i n t o unknown p a r t s o f  degree of  consider  order  (1976)) and  have t o be  has  coefficients.  criterion  Jk££iications t o S y n t h e t i c  Let  spectrum  Clayton  the  prediction  prediction f i l t e r FPE  the  process  autoregressive  however, t h e  prediction filter  energy  Akaike's f i n a l  the  of  A means o f  p r e d i c t i o n e r r o r power and  p r e d i c t i o n operator  manner.  be.  autoregressive  determines the  harmonic p r o c e s s e s ,  a l w a y s be  should  f o r the  (1969).  power i n v o l v e d i n e s t i m a t i n g For  t o know what l e n g t h  u s e d , o r e g u i v a l e n t l y , what  a s u i t a b l e order  been p r e s e n t e d  minimizing  should  autoregressive  determining  error  signal.  p r e d i c t i o n , i t i s important  prediction order  the  and  Figure  whose s p e c t r u m  modal i s f o r m e d by gives a similar  A perfectly flat  by  Data  which i l l u s t r a t e  method.  example g i v e n  Actual  33(a)  i s m u l t i p l i e d by  the  h a l f s of  t r a c e spectrum Ristow  the  t o sweep  and  s p e c t r u m i s used f o r t h e  wavelet i n t h i s case.  The  two  effect  zero  of m u l t i p l y i n g  a  152 the t r a c e spectrum o f Figure 32 (a) by the wavelet spectrum, |V(ur)|,  and  of F i g u r e 33(a) i s shown by the modified t r a c e spectrum  the modified t r a c e i n Figure 33 (a).  deconvolve the modified spike s e r i e s , good estimate  The problem i s t o  t r a c e i n order t o recover the i n p u t  s i n c e {V (or) \ i s a low guefrency of j V ( w ) | can be obtained  cepstral f i l t e r .  phenomenon, a  by using a low cut  Low c u t c e p s t r a l f i l t e r i n g  estimates the  wavelet spectrum from the amplitude spectrum of the t r a c e as shown i n F i g u r e 33 (b).  T h i s wavelet spectrum c o n t a i n s the low  g u e f r e n c i e s of the o r i g i n a l amplitude spectrum. F i g u r e 3 3 ( c ) , the estimate  o f \V(ur)\  As shown by  i s very s i m i l a r to the  a c t u a l value used i n c o n s t r u c t i o n of the s y n t h e t i c shown i n Figure 33(a).  In t h i s case, a low c u t f i l t e r of l e n g t h 23,  which was symmetrical about t=0, was a p p l i e d to a cepstrum o f l e n g t h 256.  The t r a c e cepstrum i s then d i v i d e d by t h e estimate of \I(ur)\,  of IV(itr)J t o produce an estimate  phase of the t r a c e g i v e s an estimate F o u r i e r transform, hz.  I (or),  regressive prediction f i l t e r  are f i l l e d  of l e n g t h 10.  r e s u l t i n g from c e p s t r a l f i l t e r i n g i s a l s o shown i n F i g u r e 33 (c). r e s u l t i n g deconvolution of s p i k e s .  the o r i g i n a l  of the impulse response's  i n the frequency  The unknown values of I (or)  Using  band from  18 to 50  i n by using an autoThe spectrum  and a u t o r e g r e s s i v e  prediction  As shown by F i g u r e 33 (c), the '  i s very s i m i l a r t o the d e s i r e d seguence  F i g u r e 33a).  The amplitude spectrum  of a v i b r o s e i s  wavelet i s shown a l o n g with the modified spectrum trace. the  and  The amplitude spectrum of t h i s wavelet i n c l u d e s  e f f e c t of " e a r t h  filtering".  a)  WAVELET SPECTRUM  125.0  FREQUENCY  MODIFIED TRRCE SPECTRUM 1  0.0  33.25  r  62.5  i  93.75  FREQUENCY  MODIFIED TRRCE n  I r~4  1  125,0  155  0.0  15.625  31.22  43.275  62.5  FREQUENCY  Figure 33b). the modified obtained  A comparison i s given between the spectrum of  t r a c e and an estimate  by c e p s t r a l  filtering.  of the wavelet spectrum  156  Figure 33c).  The spectrum r e s u l t i n g  f l a t t e n i n g and subseguent a u t o r e g r e s s i v e shown as well as the r e s u l t i n g amplitude spectrum estimated  from  prediction i s  deconvolution.  by. c e p s t r a l  cepstral  The  filtering i s  very s i m i l a r t o the o r i g i n a l wavelet spectrum shown i n F i g u r e 33a).  157'  RMP  C)  SPEC  33.25  REMOVED  E2.5  93.75  FREQUENCY  MEM,PRED  0.0  SPEC  r  1  T 31.25  P2.5  93.75  FREQUENCY  DECONVOLUTION  0.2  TIME  125.0  •  0.4  9  1  0.6  i  125.0  158  d)  . MODIFIED  TRACE  DECONVOLUTION  aw^t  0.0  . .  0.3  .  0.4  0.6  TIME  Figure 33d),  The a s s o c i a t i o n of a minimum phase spectrum  with the wavelet amplitude shown.  spectrum g i v e s the modified  trace  A p p l i c a t i o n of zero phase c e p s t r a l f i l t e r i n g and  a u t o r e g r e s s i v e p r e d i c t i o n g i v e s the r e s u l t i n g  deconvolution.  159 The zero  deconvolution  test  phase a s s u m p t i o n .  deconvolution  Ristow  wavelet  minimum  phase f u n c t i o n .  i s zero  with  with e ( t ) .  the effects  phase s h i f t e d  g i v e s a minimum  as a  phase The  phase s h i f t  Hilbert associated  t r a c e shown i n  o f t h e same  deconvolution  (a) s u p p l i e s a d e c o n v o l u t i o n  the r e s u l t  the  which i s  o f F i g u r e 3 3 ( c ) b u t which  defines  desired spikes. Figure  minimum and  from  a minimum  of "earth f i l t e r i n g " .  The a p p l i c a t i o n  a p p r o a c h as i n F i g u r e 33  the  of having  T h i s produces the modified  33(d).  Although  phase, e (t) i s o f t e n c o n s i d e r e d  of log(|E(w)|)  transform  using the  (1975) a p p l i e d w i e n e r  phase a s s u m p t i o n s .  us e x a m i n e t h e e f f e c t  associated  Figure  and J u r c z y k  and made minimum  Klauder  Let  i n F i g u r e 33 was p e r f o r m e d  phase w a v e l e t  then  spike.  34 shows t h e r e s u l t  using  i t i s bandlimited.  reasons  restricted  why c o n v e n t i o n a l  i n their  therefore  i t i s not minimum  seguence.  autocorrelation  from  a filtered  phase.  o f t h e e a r t h may assumption  approaches  33(d), to a  trace i s deconvolution T h e r e a r e two  minimum phase d e c o n v o l u t i o n s a r e  i s essentially  The l a t t e r  deconvolution  i n F i g u r e 33 ( d ) .  performance c a p a b i l i t i e s .  wavelet  response  shapes t h e wavelet  This results i n a  to that given  vibroseis  impulse  which  trace of Figure  t h e spectrum o f t h e deconvolved  which i s i n f e r i o r other  the modified  a Siener f i l t e r  Although  flattened,  from  of estimating a truncated  First  Klauder  Secondly,  of a l l , the wavelet,  the d e s i r e d  not approximate a white i s necessary  when e s t i m a t i n g  the  and  i n minimum wavelet  the trace a u t o c o r r e l a t i o n .  noise phase  Figure 3 4 .  Minimum phase deconvolution of the  v i b r o s e i s t r a c e model: a) Input t r a c e i d e n t i c a l to t h a t used i n F i g u r e 33d). b) T r u n c a t i o n of minimum phase wavelet  estimated  u s i n g the H i l b e r t transform method. c) Spectrum of the d e c o n v o l u t i o n . d) Deconvolution r e s u l t i n g from the a p p l i c a t i o n a Hiener i n v e r s e f i l t e r  designed  to spike the  wavelet.  of  SPECTRUM OF DECONVOLUTION  A  0.0  62.5  J25.0  FREQUENCY  DECONVOLUTION WWRVELET  ESTIMATE  i  4  0.0  •A>«"»''VW| ,  I  0 r  TJIH  1  0.4  J67.S  250.0  162 Having the  obtained  method o f c e p s t r a l  with a s e t of  a  two  filtering  data  from  and  statics  F i g u r e 35(a) initially  35(b),  cepstral  Figure  35(c)  o f the  noise.  filtering  domain e f f e c t s l e n g t h 27  filter  15  to f i n d  compromise  the  Figure  cepstral  the  high  and  43 h z .  a  and  freguency the  filter  prediction using  different  stage i n  determining s e t s of  autoregressive orders are  deconvolution  the  I n s p e c t i o n of  t r a c e s a t each  Often,  Each  frequency  msec)was d e s i g n e d  p r o v i d e s a good means o f  between r e s o l u t i o n  filtering  A cepstral  s p e c t r u m , and  deconvolved  l e n g t h s and  of  spectra of  approach.  the  (or 60  shown.  Fourier transform.  amplitude  v a l u e s between 16  a suitable  In  p r o c e s s i n g t o show t h e  p r o c e s s i n g parameters.  cepstral filter order  of  samples  process  The  i n c r e a s e s the r e s o l u t i o n  to f l a t t e n  s p e c t r a and  deconvolution  suitable  of  compares t h e  used  transform  the amplitude the  was  Standard.  s t a g e s i s shown.  of the d e c o n v o l u t i o n  of length  Fourier  Chevron  and  which were o b t a i n e d  expense of i n c r e a s i n g  at various stages  shows  divergence  deconvolution after  deconvolution  applied  F i g u r e 35  provides the deconvolution  d i s p l a y s the  F i g u r e 35(d)  traces  of  by  different  v i b r o s e i s t r a c e at the  spherical  used  'bright  d e c o n v o l u t i o n was  shows t h e i n p u t d a t a  processed  one  n o r m a l moveout c o r r e c t i o n s ,  autoregressive prediction  stage  s p e c t r a l e x t e n s i o n was  The  after  removal,  d e c o n v o l u t i o n a t two  and  s y n t h e t i c data,  p r o c e s s i n g seguence.  d e c o n v o l u t i o n of the data  stacking.  and  i s located.  stages i n the  corrections,  r e s u l t s from  v i b r o s e i s t r a c e s i n which a t l e a s t  s p o t ' gas r e f l e c t o r at  encouraging  used  in  which p r o v i d e s a good  stability.  Figure 35. applied  a f t e r NMO  a) Input stacked v i b r o s e i s and  statics  data  corrections.  (b) F i r s t stage of d e c o n v o l u t i o n a f t e r c e p s t r a l filtering. (c) Second stage of d e c o n v o l u t i o n a f t e r c e p s t r a l f i l t e r i n g and a u t o r e g r e s s i v e p r e d i c t i o n transform.  of the  Fourier  165  Figure 35 d).  Comparison of amplitude s p e c t r a ;  before d e c o n v o l u t i o n , (iii)  (ii) after cepstral  after autoregressive  prediction.  flattening,  (i) and  166 l)  S P E C T R A  BEFORE  DECONVOLUTION  i i )  AFTER  C E P S T R A L  F I L T E R I N G  i i i )  AFTER  AUTOREGRESSIVE  PREDICTION  rREOUtXCT  A  A. A  4i  A A A  A^ A_ .A. A A J  A..  A.  Jfr\  A jvli  ,»||.| ~\v  167 In Figure 36, the two step deconvolution the  process  decreases  width of a r e f l e c t i o n event which was caused by t h e presence  of n a t u r a l gas.  The l o c a t i o n of the r e f l e c t i o n from the n a t u r a l  gas d e p o s i t i s shown by the arrows on the s e i s m i c s e c t i o n s o f F i g u r e 36.  F i g u r e 36(a) shows a s e c t i o n i n which the ' b r i g h t  s p o t ' r e f l e c t i o n from a gas d e p o s i t was enhanced by Wiener deconvolution 36(b) in  and amplitude s c a l i n g .  The s e c t i o n i n F i g u r e  d i s p l a y s an a d d i t i o n a l • d e c o n v o l u t i o n of the s e c t i o n shown  F i g u r e 36 (a).  This deconvolution  was a p p l i e d t o provide a  b e t t e r r e s o l u t i o n o f the " b r i g h t spot" r e f l e c t i o n . deconvolution  was obtained  27 and a p r e d i c t i o n f i l t e r  The  by using a c e p s t r a l f i l t e r o f length 25*  of length  Deconvolution  compresses the r e f l e c t i o n frcm the gas d e p o s i t t o a s i n g l e w e l l d e f i n e d event on the s e c t i o n to give a sharper d e f i n i t i o n of the desired  horizon.  F i g u r e 36.  a) F i n a l processed  section  i n d i c a t e d " b r i g h t s p o t " emphasized by Wiener and amplitude  adjustments  with filtering  of data i n F i g u r e 35.  b) Subseguent d e c o n v o l u t i o n compresses pulse from " b r i g h t s p o t " r e f l e c t i o n shown i n F i g u r e 36a).  As i n  F i g u r e 35, the v a r i a b l e area - wiggle p l o t has i t s p o l a r i t y r e v e r s e d t o emphasize the negative r e f l e c t i o n c o e f f i c i e n t r e s u l t i n g from the gas d e p o s i t .  170 CHAPTER 6  SOgflaSI  T h i s t h e s i s i n v e s t i g a t e s s e v e r a l methods of deconvolving r e f l e c t i o n r e c o r d i n g s i n e x p l o r a t i o n seismology and makes recommendations as t o the usage of these methods. Chapter 1 presents models f o r t h e impulse response of a l a y e r e d earth and wavelets generated by e x p l o s i v e s o u r c e s .  The  s e i s m i c t r a c e i s c o n s i d e r e d to be,a c o n v o l u t i o n of t h e e a r t h ' s impulse response with a source generated wavelet.  The purpose  of d e c o n v o l u t i o n i s t o remove the wavelet and y i e l d the impulseresponse of a l a y e r e d The  earth.  problem of deconvolving e x p l o s i o n generated s e i s m i c  data i s i n v e s t i g a t e d i n Chapters 2 t o 4. e s t i m a t i o n methods are c o n s i d e r e d .  Three  wavelet  Wiener-Levinson  e s t i m a t i o n i s a minimum phase time domain approach r e g u i r e s an e s t i m a t e of the wavelet l e n g t h .  wavelet which  The Wold-Kolmogorov  e s t i m a t i o n method i s a frequency domain approach  which uses the  same assumptions as the Wiener-Levinson method.  Wold-Kolmogorov  f a c t o r i z a t i o n e s s e n t i a l l y uses a H i l b e r t transform of the l o g amplitude spectrum  to y i e l d a minimum phase spectrum.  The  Wiener-Levinson and Wcld-Kolmogorov methods provide s i m i l a r wavelet e s t i m a t e s f o r t h e s y n t h e t i c example used.  The  homomorphic d e c o n v o l u t i o n method provides a wavelet e s t i m a t e by s e p a r a t i n g the complex cepstrum trace.  of the wavelet from t h a t of the  For n o i s e l e s s s y n t h e t i c s , homomorphic d e c o n v o l u t i o n i s  extremely promising when compared to the minimum phase approaches.  Problems  with phase curve unwrapping a r i s e with  t h i s method when a d d i t i v e noise i s present i n the s e i s m i c t r a c e .  171 Thus some form of c e p s t r a l s t a c k i n g i s u s e f u l when e s t i m a t i n g wavelets from noisy data.  On a c t u a l e x p l o s i o n data, the  wavelets o b t a i n e d by homomorphic d e c o n v o l u t i o n a r e c l e a r l y n o t minimum phase, c a u s i n g the minimum phase assumption  t o remain  suspect. A f t e r e s t i m a t i o n o f a s e i s m i c wavelet,  the problem of  d e s i g n i n g s i n g l e channel and multichannel Wiener d e c o n v o l u t i o n filters  i s e s s e n t i a l l y one o f choosing a s u i t a b l e f i l t e r  and a l e v e l of prewhitening. deconvolution f i l t e r when the wavelet  Although  length  the m u l t i c h a n n e l Wiener  i s s u p e r i o r t o the s i n g l e channel  filter  i s known, the d i f f e r e n c e s i n the methods  becomes much l e s s s i g n i f i c a n t when estimated  wavelets  and n o i s y  data are used. D i f f e r e n c e s i n the Hiener and MEM deconvolution r e s u l t from the d i f f e r e n t a u t o c o r r e l a t i o n s being design.  filters  used i n t h e i r  While the d i f f e r e n c e i n a u t o c o r r e l a t i o n e s t i m a t e s  becomes i n s i g n i f i c a n t f o r t y p i c a l gate l e n g t h s used i n s e i s m i c p r o c e s s i n g , the MEM approach f o r c a l c u l a t i n g a time a d a p t i v e filter  i s p r e f e r a b l e when using s h o r t data segments.  Deconvolution in  Chapter  4.  of n o i s y dynamite generated data i s presented  Although the d e c o n v o l u t i o n r e s u l t s were very  s i m i l a r to those obtained by another problems are encountered impulse response  independent  when attempting  r e s e a r c h group,  to c o r r e l a t e the  s y n t h e t i c s and s e i s m i c data.  The reasons f o r  these problems l i e both i n the data q u a l i t y and the i n h e r e n t d i f f e r e n c e s between what the well l o g s and the s e i s m i c r e f l e c t i o n s are measuring.  172  In Chapter 5 ,  the p r o p e r t i e s of v i b r o s e i s r e c o r d i n g s  e x p l o i t e d to develop a new cepstral f i l t e r i n g  approach to deconvolution  and a u t o r e g r e s s i v e p r e d i c t i o n .  uses the assumption of a low  guefrency,  zero  are  which uses  T h i s approach  phase wavelet.  The  zero phase assumption avoids the problem of phase unwrapping which i s o f t e n encountered when using homomorphic on r e a l data, and i s j u s t i f i e d  by the f a c t t h a t the v i b r o s e i s  t r a c e i s composed of f i l t e r e d Klauder earth f i l t e r ,  e (t),  is successful.  deconvolution  wavelets.  Even when the  i s minimum phase, the d e c o n v o l u t i o n  Cepstral f i l t e r i n g  removes the low  guefrency  c o n t r i b u t i o n of the wavelet to the spectrum w i t h i n a passband.  An a u t o r e g r e s s i v e  extend the F o u r i e r transform l i m i t s of the FM sweep.  restricted  p r e d i c t i o n f i l t e r i s then used to of the impulse response beyond  T h i s s p e c t r a l extension  on the i d e a that the F o u r i e r transform f u n c t i o n s can  method  of a sum  the  scheme i s based of d e l t a  be p r e d i c t e d by using an a u t o r e g r e s s i v e  model.  S i n c e the bandwidth of the v i b r o s e i s s i g n a l i s w e l l known from a knowledge of sweep f r e q u e n c i e s , i t i s not d i f f i c u l t t o  decide  which p o r t i o n of the spectrum i s to be used i n d e s i g n i n q prediction f i l t e r .  As shown by the s y n t h e t i c and  analysed, the deconvolution vibroseis traces.  provides  real  the  data  qood r e s o l u t i o n of  Since the method does not r e g u i r e  the  assumption of a random impulse response, i t can be used on s h o r t segments of the s e i s m i c t r a c e .  A l s o , because of the  r e s o l u t i o n p o s s i b i l i t i e s of the method, i t can be  high  particularly  u s e f u l i n r e s o l v i n g v i b r o s e i s r e f l e c t i o n s from t h i n beds.  173 REFERENCES Akaike,H., 1969, prediction: 247.  F i t t i n g a u t o r e g r e s s i v e models f o r Ann. I n s t . Statist. Math. , v . 2 1 ,  p.243-  A n d e r s e n , N. , 1974, On t h e c a l c u l a t i o n o f f i l t e r c o e f f i c i e n t s f o r maximum e n t r o p y s p e c t r a l a n a l y s i s : G e o p h y s i c s , v.39, p.69-72. B e r r y m a n , L.H., G o u p i l l a u d , P.L., and W a t e r s , K.H., 1 958, R e f l e c t i o n s from m u l t i p l e t r a n s i t i o n l a y e r s : Geophysics, v.23, p.223-243. Bogert,  R.P., H e a l y , M.J., and l u k e y , J . W . , 1963, The g u e f r e n c y a l a n y s i s o f time s e r i e s f o r echoes: Proc. Symp. On Time S e r i e s A n a l y s i s , M. R o s e n b l a t t , e d . , New Y o r k , S i l e y , p.209-243.  Bracewell,R., 1965, applications:  The F o u r i e r t r a n s f o r m and i t s New Y o r k , McGraw - H i l l .  B u h l , P . , S t o f f a , P . L . . And Bryan,G.M., 1974, The a p p l i c a t i o n o f homomorphic d e c o n v o l u t i o n t o s h a l l o w - w a t e r m a r i n e s e i s m o l o g y - p a r t 2: r e a l d a t a : G e o p h y s i c s , v . 3 9 , p. 417426. B u r g , J . P . , 1967, Maximum e n t r o p y s p e c t r a l a n a l y s i s : p a p e r p r e s e n t e d a t t h e 1967 SEG m e e t i n g i n Oklahoma C i t y , Okla. B u r g , J . P . , 1975, Maximum e n t r o p y s p e c t r a l a n a l y s i s : Stanford Univ., Palo A l t o , C a l i f . Buttkus,  B., Homomorphic f i l t e r i n g : p. 723-748.  Geophys.  thesis,  Prosp., v.23,  C l a e r b o u t , J . , 1976, F u n d a m e n t a l s o f g e o p h y s i c a l d a t a New Y o r k , M c G r a w - H i l l . Clarke,G.K.C. 1968, Time v a r y i n g d e c o n v o l u t i o n G e o p h y s i c s , v . 3 3 , p.936-944. Clayton,  PhD  processing:  filters:  R.W., and U l r y c h , T . J . , 1976, R e s t o r a t i o n o f i m p u l s i v e f u n c t i o n s : IEEE Trans. Info. Theory, submitted f o r publication.  Clayton,R.W., M c C l a r y , B . , and l i g g i n s , R . A . , 1975, Comments on •Phase d i s t o r t i o n and H i l b e r t t r a n s f o r m a t i o n i n m u l t i p l y r e f l e c t e d and r e f r a c t e d body waves' by G.L. Choy and P.G. R i c h a r d s : B u l l . Seism. S o c . Am., i n p r e s s .  174  Clayton,R.W., and Wiggins, R.A., 1976, Source shape e s t i m a t i o n and deconvolution of t e l e s e i s m i c bodywaves: Geophys. J.R. Astr. S o c , i n press. Clowes, R.M., 1 9 6 9 , Seismic r e f l e c t i o n i n v e s t i g a t i o n s of c r u s t a l s t r u c t u r e i n Southern A l b e r t a : PhD t h e s i s , Univ. Of A l b e r t a , Edmonton, A i t a . Crump,N., 1 9 7 4 , seismic  A Kalman f i l t e r approach to the d e c o n v o l u t i o n of s i g n a l s : Geophysics, v . 3 9 , p . 1 - 1 3 .  Davies,E,B., and Mercado, E.J., 1 9 6 8 , Multichannel d e c o n v o l u t i o n f i l t e r i n g of f i e l d recorded s e i s m i c data: Geophysics, v.33, p.711-722. Davies, J.C., 1 9 7 6 , Maximum entropy s p e c t r a l a n a l y s i s of f r e e o s c i l l a t i o n s of the e a r t h : M.Sc. T h e s i s , U n i v e r s i t y of B r i t i s h Columbia. DeVoogd, N., 1 9 7 4 , Wavelet shaping and Prosp,, v . 2 2 , p . 3 5 4 - 3 6 9 .  noise  reduction:  Geophys.  F r y e r , G . J . , Odegard,M.E., and Sutton,G.H., 1 9 7 5 , Deconvolution and s p e c t r a l e s t i m a t i o n using f i n a l p r e d i c t i o n e r r o r : Geophysics, v . 4 0 , p . 4 1 1 - 4 2 5 . Grant,F., and West,G., 1 9 6 5 , I n t e r p r e t a t i o n theory i n geophysics: Toronto, McGraw - H i l l . Futterman, W.I., 1962, Dispersive v . 6 7 , p. 5279-5291.  applied  body waves: J.Geoph.Res.,  G a l b r a i t h , J . N . , 1 9 7 1 , P r e d i c t i o n e r r o r as a c r i t e r i o n f o r operator l e n g t h : Geophysics, v . 3 5 , p . 2 5 1 - 2 6 5 . Gurbuz, B.M., 1 9 7 2 , S i g n a l enhancement of v i b r a t o r y s o u r c e data i n the presence of a t t e n u a t i o n : Geophys. Prosp., v . 2 0 , p.421-438. Jenkins,G.M. And Watts, D.G., 1969, a p p l i c a t i o n s : San F r a n c i s c o ,  S p e c t r a l A n a l y s i s and i t s Holden-Day.  Kanasewich, E.R., 1 9 7 3 , Time sequence a n a l y s i s i n geophysics: U n i v e r s i t y of A l b e r t a Press, Edmonton, A l b e r t a .  175 Klauder,J,R., P r i c e , A.C., D a r l i n g t o n , S., and A l b e r s h e i m , W.J. 1960, The t h e o r y and d e s i g n o f c h i r p r a d a r s ; B a l l System Tech. J o u r . , v.39, p. 145-807. f  L e v i n s o n , N., 1947, The Wiener RMS ( r o o t mean square) e r r o r c r i t e r i o n i n f i l t e r d e s i g n and p r e d i c t i o n : J o u r . Of Math and P h y s i c s , v.25, p.261-278. Lines,  L.R., 1974, A n o t e on t h e a p p l i c a t i o n multichannel deconvolution: J . Can. G e o p h y s i c i s t s , v.10, p.65-70.  of Wiener Soc. Expl.  Lines,  L.R. And U l r y c h , T . J . , 1976, T e c h n i g u e s f o r s e i s m i c d e c o n v o l u t i o n and w a v e l e t e s t i m a t i o n ; p a p e r p r e s e n t e d a t t h e 1975 SEG m e e t i n g i n D e n v e r , C o l o r a d o and a c c e p t e d f o r p u b l i c a t i o n by Geophys. Prosp. Pending r e v i s i o n .  Lines,  L.R. And C l a y t o n , R.W., 1976, A new a p p r o a c h f o r v i b r o s e i s d e c o n v o l u t i o n : paper presented at the 1976 CSEG c o n v e n t i o n i n C a l g a r y , A l b e r t a , and s u b m i t t e d f o r p u b l i c a t i o n t o Geophys. Prosp.  Oppenheim, A.V., 1965, S u p e r p o s i t i o n i n a c l a s s o f n o n - l i n e a r s y s t e m s ; R e s e a r c h Lab o f E l e c t r o n i c s , MIT, T e c h . Rep. 432. . Oppenheim, A.V., S c h a f e r , R.W., and Stockham, T.G., 1968, N o n l i n e a r f i l t e r i n g o f m u l t i p l i e d and c o n v o l v e d s i g n a l s : P r o c . I E E E , v.65, p.1264-1291. Otis,  S.M. And S m i t h , R.B., 1974, Homomorphic d e c o n v o l u t i o n by l o g s p e c t r a l a v e r a g i n g : p a p e r p r e s e n t e d a t t h e 1974 SEG meeting i n D a l l a s , Texas.  P e a c o c k , K.L., theory Ricker, Riley,  and and  N., 1953, wavelets:  T r e i t e l , S., 1969, P r e d i c t i v e d e c o n v o l u t i o n : p r a c t i c e : G e o p h y s i c s , v.34, p.155-169. The f o r m and l a w s o f p r o p o g a t i o n G e o p h y s i c s , v.18, p.10-40.  of  seismic  D. And B u r g , J.P., 1972, Time and s p a c e a d a p t i v e d e c o n v o l u t i o n f i l t e r s : paper p r e s e n t e d a t t h e 42nd m e e t i n g o f t h e SEG, Anaheim, C a l i f . , Nov. 28, 1972.  R i s t o w , D., and J u r c z y k , D., 1974, V i b r o s e i s Geophys. P r o s p . , v.23, p.363-379. R o b i n s o n , E,A., digital  deconvolution:  1967a, M u l t i c h a n n e l time s e r i e s a n a l y s i s w i t h computer p r o g r a m s : San F r a n c i s c o , H o l d e n - D a y .  176 Robinson, E.A., 1967b, P r e d i c t i v e decomposition of time s e r i e s with a p p l i c a t i o n to s e i s m i c e x p l o r a t i o n : Geophysics, v.32, p.418-484. Robinson, E.A., and T r e i t e l , S., 1967, P r i n c i p l e s of d i g i t a l Wiener f i l t e r i n g : Geophys. Prosp., v.15, p.311-333. Robinson, E.A., 1975a, Deconvolution methods and s e i s m i c models: paper presented a t the 1S75 SEG meeting i n Denver, Colorado. Robinson, E.A., 1975b, Dynamic p r e d i c t i v e d e c o n v o l u t i o n : Geophys. Prosp., v.23, p.780-798. S c h a f e r , R.W., 1969, Echo removal by d i s c r e t e g e n e r a l i z e d l i n e a r f i l t e r i n g : Research Lab o f E l e c t r o n i c s , MIT, Tech Rep, 466. Smylie, D,E., C l a r k e , G.K.C., and U l r y c h , T.J., 1973, A n a l y s i s of i r r e g u l a r i t i e s i n the earth's r o t a t i o n : Methods of Comp. P h y s i c s , v.13, Acad. Press Inc., New York. S t o f f a , P., Buhl, P., and Bryan, G., 1974, The a p p l i c a t i o n o f homomorphic d e c o n v o l u t i o n to shallow water marine seismology - p a r t 1: theory: Geophysics, v.39, p.417436. Treitel,  S., and Robinson, E.A., 1966, The design of high r e s o l u t i o n d i g i t a l f i l t e r s : IEEE Trans. On G e o s c i . E l e c t r o n i c s , v.GE-4, p25-38.  Treitel,  S., 1970, P r i n c i p l e s of d i g i t a l multichannel f i l t e r i n g : Geophysics, v.35, p.785-811.  Treitel,  S., 1975, An i d e a l performance property of the multichannel Wiener f i l t e r : J . Can. Soc. Of E x p l . G a o p h y s i c i s t s , v.10, p.71-74.  U l r y c h , T . J . , 1971, A p p l i c a t i o n of homomorphic d e c o n v o l u t i o n to seismology: Geophysics, v,36, p.650-660. U l r y c h , T . J . , 1972, Maximum entropy power spectrum of t r u n c a t e d s i n u s o i d s : J . Geophys. Res., v.77, p.1396-1400. U l r y c h , T . J . , Jensen, O.J., E l l i s , R.M., and S o m e r v i l l e , P. G. , 1972, Homomorphic d e c o n v c l u t i o n of some t e l e s e i s m i c events: B u l l . Seism. Soc. Am., v.62, p.1269-1281.  177 U l r y c h , T . J . , Smylie, B,E., Jensen, O.G., and C l a r k e , G. K.C,, 1973, P r e d i c t i v e f i l t e r i n g and smoothing of s h o r t r e c o r d s by using maximum entropy: J.Geophys.fies., v.77, p.1396-1400. U l r y c h , T.J. And Bishop, T., 1975, Maximum entropy s p e c t r a l a n a l y s i s and a u t o r e g r e s s i v e decomposition: Reviews o f Geophysics, v.13, p.183-200. U l r y c h , T . J , , and C l a y t o n , R.W., 1976, Time s e r i e s modelling and maximum entropy: Phys. E a r t h P l a n e t . I n t . , i n press. Wang, R., 1969, The determination o f optimum gate l e n g t h s f o r time v a r y i n g Wiener f i l t e r i n g : Geophysics, v.34, p.683695, White, R.E., and O'Brien, P.U.S., 1974, E s t i m a t i o n of the primary s e i s m i c p u l s e : Geophys. Prosp., v.22, p.627651. Wiggins,  R.A., and Robinson, E.A., 1965, Recursive s o l u t i o n t o the multichannel f i l t e r i n g problem: J . Geoph Res., v,70, p.1885-1891.  Wiggins,  R.A,, 1968, The use o f expected e r r o r i n the design of l e a s t squares optimum f i l t e r s , Geophys. Prosp., v.15, p.288-296.  Wood, L.C. And T r e i t e l , S,, 1975, Seismic s i q n a l p r o c e s s i n g : Proc. IEEE, v.€3, p.649-661. Wood, L . C , 1976, High r e s o l u t i o n deconvolution of symmetrical pulses using Wiener f i l t e r s : paper presented a t the 1976 CSEG convention i n C a l g a r y , A l b e r t a . Y e d l i n , M., 1975, Kalman f i l t e r i n g and t h e d e c o n v o l u t i o n .problem: unpublished manuscript.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052513/manifest

Comment

Related Items