HRS A New Look at Seismic Attributes [PDF]

Citation preview

A new look at seismic attributes Brian Russell



Introduction  Seismic attributes have evolved into a set of seemingly unrelated methods:     



Instantaneous attributes, Windowed frequency attributes, Coherency and semblance attributes, Curvature attributes, Phase congruency, etc.



 In this talk, I want to make sense of all these attributes by showing how they are inter-related.  In particular, I will “close the loop” between instantaneous attributes and trace-by-trace correlation attributes.  To illustrate the various methods, I will use the Boonsville dataset, a karst limestone example from Texas. November, 2012



2



The Boonsville dataset



November, 2012



Hardage et al., 1996



 The Boonsville gas field in north Texas (Hardage et al., 1996) is controlled by karst collapse features, making it ideal for studying attributes.  A dataset over the field was made available by the Texas Bureau of Economic Geology.  The study area is shown on the left. 3



Boonsville Geology  In the Boonsville gas field, production is from the Bend conglomerate, a middle Pennsylvanian clastic deposited in a fluvio-deltaic environment.  The Bend formation is underlain by Paleozoic carbonates, the deepest being the Ellenburger Group of Ordovician age.  The Ellenburger contains numerous karst collapse features which extend up to 760 m from basement through the Bend conglomerate, shown here in schematic:



November, 2012



Lucia (1995)



4



A seismic volume  The 3D seismic volume is shown here in grey level variable density format.  It consists of 97 inlines and 133 crosslines, each with 200 samples (800– 1200 ms).  The karst features are illustrated by the red ellipses. November, 2012



Crossline or y direction



Inline or x direction Time or t direction



5



Slices through the volume Here are vertical and horizontal slices through the seismic volume:



Inline 146 1000 ms



Inline seismic section 146 November, 2012



Time slice at 1000 ms



Note the karst feature at the east end of the line.



6



A brief history of attributes  The first attributes were single trace instantaneous attributes, introduced by Taner et al. (1979).  Barnes (1996) extended these attributes to 2D and 3D.  From 1995-1999, Amoco researchers introduced two new types of spatial attributes:  Spectral decomposition, based on the Fourier transform.  Coherency, based on trace-to-trace correlation.



 Roberts (2001) introduced curvature attributes.  Marfurt et al. (2006) then showed the relationship between instantaneous attributes and dip and curvature.  We shall also discuss the phase congruency attribute (Kovesi, 1996) which was developed in image analysis.  A good place to start is with the Fourier transform. November, 2012



7



The Fourier transform  The Fourier transform (FT) takes a real time signal s(t) and transforms it to a complex frequency signal S(w):



s(t )  S (w )  S (w ) e (w ) , where



w  2f , f  frequency, S (w )  the Fourier transform of s(t ), S (w )  amplitude spectrum, and  (w )  phase spectrum, of s(t ).  We can also perform the inverse transform from frequency back to time.  In using the FT for attribute analysis, note that it operates on the complete seismic trace, so cannot be localized at a time sample to give an “instantaneous” attribute value.  Partyka et al. (1999) used the FT over small seismic windows to compute the spectral decomposition attribute. November, 2012



8



Spectral decomposition



y



y



x time



x Fourier Transform



frequency



As shown here, spectral decomposition involves taking a small window around a zone of interest (here: 1000 ms +/- 50 ms) and transforming each trace window to the frequency domain. November, 2012



9



Spectral slices



Data slice



5 Hz



10 Hz



25 Hz



30 Hz



35 Hz



15 Hz



40 Hz



20 Hz



45 Hz



Spectral slices at various frequencies, compared with seismic amplitude. November, 2012



10



Spectral decomposition limitations  As seen in the previous slide, each frequency slice in a spectral decomposition reveals a different underlying channel structure in the original seismic amplitude slice.  However, this can be seen as a limitation since we really don’t know which slice to use.  Also, as the size of the window decreases, the frequency resolution also decreases.  This resolution was improved by Sinha et al. (2005) using the continuous wavelet transform.  However, the fact that each frequency slice gives a different result still remained.  This leads us to the concept of instantaneous attributes, including instantaneous frequency. November, 2012



11



The complex time signal  Instantaneous attributes are based on the concept of the complex time signal, invented by Gabor in 1947 in a classic paper called “Theory of communication”.  Since the Fourier transform of a real signal has a symmetric shape on both the positive and negative frequency side, Gabor proposed suppressing the amplitudes belonging to negative frequencies, and multiplying the amplitudes of the positive frequencies by two.  An inverse transformation to time gives the following complex signal: z (t )  s(t )  ih (t ) where :



z (t )  the complex signal, i   1, and h(t )  the Hilbert transform of s(t ). November, 2012



12



The Hilbert transform  The Hilbert transform is a filter which applies a 90o phase shift to every sinusoidal component of a signal.  Notice that the complex trace can be transformed from rectangular to polar coordinates, as shown below, to give the instantaneous amplitude A(t) and instantaneous phase F(t): z (t )  s(t )  ih (t )  A(t )eiF ( t )  A(t ) cos F (t )  iA(t ) sin F (t ), where :



h(t) A(t)



f(t) s(t)



time November, 2012



 h(t )  A(t )  s(t )2  h(t )2 , F (t )  tan 1  .   s(t ) 



 The complex trace was introduced into geophysics by Taner et al. (1979), who also discussed its digital implementation. 13



Instantaneous frequency  Taner et al. (1979) also introduced the instantaneous frequency of the seismic trace, which was initially derived by J. Ville in a 1948 paper entitled: “Théorie et applications de la notion de signal analytique”.  The instantaneous frequency is the time derivative of the instantaneous phase: dh(t ) ds(t ) s(t )  h(t ) dF(t ) d tan1 (h(t ) / s(t ) ) dt dt w (t )    A(t )2 dt dt  Note that to compute w(t) we need to differentiate both the seismic trace and its Hilbert transform.  Like the Hilbert transform, the derivative applies a 90o phase shift, but it also applies a high frequency ramp. November, 2012



14



Instantaneous attributes +1.0



+1.0



This figure shows the instantaneous attributes associated with the seismic amplitude slice at 1000 ms.



Data slice



-1.0 +180o



November, 2012



Inst. Phase



-180o



Inst. Amp. Inst. Freq.



0.0 120 Hz



0 Hz



15



New attributes  Attribute analysis until the mid-90s was thus based on the three basic attributes introduced by Gabor and Ville in the 1940s: instantaneous amplitude, phase, and frequency.  Then, in the space of two years, papers on two new approaches to attribute analysis appeared:  The coherency method (Bahorich and Farmer, 1995)  2-D (and 3-D) complex trace analysis (Barnes, 1996)  The coherency method was a new approach which utilized cross-correlations between traces.  The paper by Barnes actually showed how instantaneous attributes and aspects of coherency were related.  However, let us first discuss the coherency method. November, 2012



16



The coherency method  The first coherency method involved finding the maximum correlation coefficients between adjacent traces in the x and y directions, and taking their harmonic average.  Marfurt et al. (1998) extended this by computing the semblance of all combinations of J traces in a window.  This involves searching over all x dips p and y dips q, over a 2M + 1 sample window: x y q dip November, 2012



Marfurt et al. (1998)



2M+1 samples



p dip 17



Coherency mathematics  If the covariance matrix between all locations i and j is:



 c11  c1J  Mt C ( p, q)      , cij   si (t  pxi  qyi ) s j (t  px j  qy j ),   t   Mt cJ 1  cJJ  then the two coherency measures are as follows: 1/ 2   a T C ( p, q ) a  c12 c13     and coh2  max  coh1  max  ,  1/ 2 1/ 2   (c11c22 ) (c11c33 )    TrC ( p, q)  where a T  1, , 1 and TrC ( p, q)  sum of main diagonal of C ( p, q).



 A third measure, called eigen-coherence (Gursztenkorn and Marfurt, 1999), is:   1 coh3  max  , 1  first eigenvalue of C ( p, q).  November, 2012  TrC ( p, q)



18



Coherency slice 1.1



0.7 +1.0



The x, y and z slices through a coherency volume with the 1000 ms data and coherency slices on the right. Note the low amplitude discontinuities and the highlighted event. November, 2012



Data slice



-1.0



19



Dip and azimuth attributes  At first glance, there appears to be no relationship between instantaneous attributes and coherency.  However, notice that in the coherency method it is important to find the x and y dips, called p and q.  The paper by Barnes (1996) introduced two new spatial instantaneous attributes and showed how these attributes could be used to find both the dips and the azimuth.  These concepts are even more important when we look at the relationship between instantaneous attributes and curvature.  Barnes (1996) also discussed instantaneous bandwidth and phase versus group velocity based on the work of Scheuer and Oldenburg (1988).  However, in the following slides I will only discuss his new spatial attributes.



November, 2012



20



New instantaneous attributes  Re-visiting our earlier 3D dataset, notice that we only computed the frequency attribute in the time direction  Since the Hilbert transform is actually a function of three coordinates (i.e. F(t,x,y)) we can also compute spatial frequencies in the inline (x) and crossline (y) directions. November, 2012



Crossline or y direction



Inline or x direction



Time or t direction



21



Instantaneous wavenumber  Analogous to instantaneous frequency, Barnes (1996) therefore defined the instantaneous wavenumbers kx and ky (using the Marfurt (2006) notation):



h s s h F(t , x, y ) kx   x 2 x , and x A h s s h F (t , x, y ) y y ky   . 2 y A  As with instantaneous frequency, the derivative operation can be done either in the frequency domain or using finite differencing up to a given order. November, 2012



22



Dip and azimuth attributes  The instantaneous time dips in the x and y direction, p and q, are given as:



p



kx



w







T



x



and q 



ky



w







T



y



,



where : T  period and   wavelengt h.  Using the dips p and q, the azimuth f and true time dip q are then given by:



f  tan1 ( p / q)  tan1 (k x / k y ), and q  November, 2012



p 2  q2 . 23



The dip attribute in 2D



 = 36 traces T = 50 ms dip = 1.4 ms/trace



To understand this concept, I created a seismic volume that consisted of a dipping cosine wave. The period, wavelength and time dip are shown on the plot. November, 2012



24



Dip and azimuth attributes  Next, let’s look at the dipping cosine wave in 3D.  In this display, we have sliced it along the x, y and time axes.  For this dipping event, it is clear how the dips and azimuths are related to the instantaneous frequency and wavenumbers. November, 2012



ky



kx q



w



p



ky



q



kx



w f



ky kx



25



Examples  Next, we will illustrate these instantaneous spatial attributes using the Boonsville dataset.  There are many possible displays based on the “building blocks” for these attributes.  These include the derivatives of seismic amplitude, Hilbert transform and phase in the t, x, and y directions, the three frequencies, and p, q, true dip and azimuth.  Note that the many divisions involved in the computations make this method very sensitive to noise when compared to correlation based methods.  In the following displays, I show the true dip and azimuth. November, 2012



26



True dip volume



Seismic with 3 x 3 alpha-trim mean November, 2012



True dip 27



Azimuth volume



Seismic with 3 x 3 alpha-trim mean November, 2012



Azimuth 28



Curvature attributes  Roberts (2001) shows that curvature can be estimated from a time structure map by fitting the local quadratic surface given by:



t ( x, y )  ax 2  by 2  cxy  dx  ey  f  This is a combination of an ellipsoid and a dipping plane.



November, 2012



29



Curvature attributes  Roberts (2001) computes the curvature attributes by first picking a 3D surface on the seismic data and then finding the coefficients a through f from the map grid shown below: x t1



t2



t3



t4



t5



t6



t7



t8



t9



y



 As an example, the coefficient d, which is the linear dip in the x direction, is computed by:



t (t3  t6  t9 )  (t1  t4  t7 ) d  x 2x



 Klein et al. (2008) show how to generalize this to each point on the seismic volume by finding the optimum time shift between pairs of traces using cross-correlation. November, 2012



30



Curvature attributes  The coefficients d and e are identical to dips p and q defined earlier, so when a = b = c = 0, we have a dipping plane and can also define the true dip and azimuth as before.  For a curved surface, Roberts (2001) defines the following curvature attributes (Kmin and Kmax are shown on the surface): 2 K max  K mean  K mean  K gauss , 2 K min  K mean  K mean  K gauss,



4ab  c 2 where : K gauss  , and 2 2 1/ 2 (1  d  e ) K mean November, 2012



a (1  e2 )  b(1  d 2 )  cde  . 2 2 3/ 2 (1  d  e ) 31



Curvature attributes  Note that the inverse relationships are: K  K max Kmean  min and K gauss  K min K max . 2  Also, can compute most positive and negative curvature K+ and K-, which are equal to Kmin and Kmax with d = e = f = 0 (or, the eigenvalues of the quadratic involving a, b and c): K  (a  b)  (a  b)2  c2 and K  (a  b)  (a  b)2  c2 .



 Finally, we have Keuler, where d is an azimuth angle: Keuler (d )  K min cos 2 d  K max sin 2 d  Kmean  Keuler (45o ), K gauss  Keuler (0o ) Keuler (90o )



 The next few slides show some examples from our dataset. November, 2012



32



Minimum Curvature



Seismic with 3 x 3 alpha-trim mean November, 2012



Minimum curvature 33



Maximum Curvature



Seismic with 3 x 3 alpha-trim mean November, 2012



Maximum curvature 34



Azimuth from Curvature



Seismic with 3 x 3 alpha-trim mean November, 2012



Azimuth from curvature by correlation 35



Azimuth comparison



Instantaneous Azimuth November, 2012



Azimuth from curvature by correlation 36



Curvature and instantaneous attributes  Differentiating the Roberts quadratic, we find that at x = y = 0 we get the following relationships for the coefficients:



1  p  1  q  1  p q  a   , b   , c    , d  p and e  q. 2  x  2  y  2  y x   This leads to the following quadratic relationship:



1  p  2 1  q  2 1  p q  t ( x, y )    x    y     xy  px  qy  f 2  x  2  y  2  y x   Thus, all of the curvature attributes can be derived from the instantaneous dip attributes described earlier, using a second differentiation.  The next figure shows a comparison between maximum curvature derived the two different ways.



November, 2012



37



Max Curvature Comparison High



Low Data slice



Instantaneous maximum curvature



Correlation maximum curvature



Notice that the curvature features are similar, but that instantaneous curvature shows higher frequency events. November, 2012



38



Possible instability  However, there is a lot of complexity hidden in the second differentiation.  To see this, let’s expand the p term:



s   h s   h   s  h   s  h   x   t t  p  (k x / w )  x  , where :   x x x s  seismic volume and h  Hilbert transform volume.  Performing this differentiation will produce a large number of seismic and Hilbert transform derivatives in both the numerator and denominator, which can cause instability in some datasets. November, 2012



39



Amplitude curvature  To avoid this instability, Chopra (2012) recommends a new approach to curvature called amplitude curvature (as opposed to the previous structural curvature), which involves first and second derivatives of only the seismic data:



s s  2 s  2 s , , 2, 2  x y x y s s s 1  2s 2s  a   cos  sin and emean   2  2 .  x y 2  x y   In the above options, a represents the amplitude curvature at an azimuth angle  and emean is the mean amplitude curvature.  Examples are shown in the next slides. November, 2012



40



X derivative



Seismic with 3 x 3 alpha-trim mean November, 2012



Amplitude derivative in x-direction 41



Y derivative



Seismic with 3 x 3 alpha-trim mean November, 2012



Amplitude derivative in y-direction 42



o



45 azimuth amp curvature



Seismic with 3 x 3 alpha-trim mean November, 2012



Amplitude curvature at 45o azimuth 43



The phase congruency method  This concludes our discussion of the most common attributes available to interpreters.  But I want to finish with yet another spatial attribute, which is not as commonly used.  This attribute is called the phase congruency algorithm (Kovesi, 1996), and was developed to detect corners and edges on 2D images.  The algorithm is applied on a slice-by-slice basis to the seismic volume and produces results that are similar to coherency.  The following few slides give a brief overview of the theory and then an application to our dataset. November, 2012



44



The phase congruency method  The concept behind the phase congruency algorithm in the 1D situation is shown below: Imaginary An



E(x)



f(x)



f (x) Real



Note that the Fourier components in the above figure are all in phase for a step. November, 2012



By lining these components up as vectors, we can measure the phase congruency. Figures from Kovesi, 1996



45



Phase congruency flow chart Input data slice in x-y space



Create N radial filters in kx-ky space



Create M angular filters in kx-ky space



2D FFT



Data slice in kx-ky space



Multiply to create N*M filters



Apply filters Inverse 2D FFT Normalize and sum radial terms for each angle



Minimum moment = corners



Apply moment analysis



Maximum moment = edges



Phase Congruency on 3D Seismic Seismic Volume



Time slices



FFT



Frequency partitioning Phase spectrum at Freq Partitions



Map the Phase



The phase congruency algorithm is applied to seismic slices, as shown above. November, 2012 47



Movie of seismic slices



Animation of seismic amplitude Animation of seismic amplitude for each vertical slice in the for each time slice in the Boonsville dataset. Boonsville dataset. 48



Coherency and phase congruency



Animation of seismic phase congruency for each time slice in the Boonsville dataset.



Animation of seismic coherency for each time slice in the Boonsville dataset. 49



Seismic vs Phase Congruency



A vertical slice showing karst features superimposed on a horizontal slice at 1080 ms, roughly halfway through the karst collapse. Seismic on left and the phase congruency on right. Note the good definition of the collapse on phase congruency. November, 2012



50



Seismic vs Coherency



A vertical slice showing karst features superimposed on a horizontal slice at 1080 ms, roughly halfway through the karst collapse. Seismic on left and coherency on right. Note the definition of the collapse on the coherency volume. November, 2012



51



Conclusions  In this document, I have given an overview of most of the commonly used attributes and how they are related.  The most basic attributes are derived from the Fourier transform, but cannot be localized.  Instantaneous attributes involve combinations of the derivatives of the seismic amplitude volume and its Hilbert transform, and were initially just done in the time direction.  Correlation attributes involve computing the cross-correlation between pairs of traces in the inline and crossline directions.  Coherency is based on the correlation amplitude and curvature on the correlation time-shift.  Instantaneous attributes in all three seismic directions can be combined to give dip, azimuth and curvature. November, 2012



52



Conclusions (continued)  Finally, we discussed the phase congruency algorithm which involved both the 2D Fourier transform and radial filters, and looks at where the Fourier phase components line up.  I illustrated these methods using the Boonsville dataset, which was over a gas field trapped in karst topography.  It is important to note that the seismic volume itself is probably our best seismic attribute!  However, each of the attributes that I discussed has its own advantages and disadvantages when used to extend the interpretation process.  By helping to show the relationships among the various attributes, I hope that I have shed some light on when these attributes can be used to best advantage by interpreters. November, 2012



53



Acknowledgements  I first met instantaneous attributes while working at Chevron Geophysical in Houston in 1979 and am indebted to Tury Taner, Fulton Koehler and Bob Sheriff for introducing them.  I became fascinated by the coherency and spectral decomp attributes introduced by Mike Bahorich, Kurt Marfurt, Greg Partyka and their colleagues at Amoco in the 1990s, but did not at first see their relationship to the earlier methods.  I also failed to see the importance of Art Barnes’ work on 2D and 3D instantaneous attributes when it first appeared, but now appreciate how important his work has been.  The ongoing work by Professor Kurt Marfurt and his students also provided the inspiration for much of this talk.  Finally, my thanks go to Satinder Chopra for both his work on attributes and his book, co-authored with Kurt Marfurt. November, 2012



54



References al-Dossary, S., and K. J. Marfurt, 2006, 3D volumetric multispectral estimates of reflector curvature and rotation: Geophysics, 71, 41–51. Bahorich, M. S., and S. L. Farmer, 1995, 3-D seismic discontinuity for faults and stratigraphic features, The coherence cube: The Leading Edge, 16, 1053–1058. Barnes, A. E., 1996, Theory of two-dimensional complex seismic trace analysis: Geophysics, 61, 264–272. Cohen, L., 1995, Time-Frequency Analysis: Prentice-Hall PTR.



Gabor, D., 1946, Theory of communication, part I: J. Int. Elect. Eng., v. 93, part III, p. 429-441. Gersztenkorn, A., and K. J. Marfurt, 1999, Eigenstructure based coherence computations as an aid to 3D structural and stratigraphic mapping: Geophysics, 64, 1468–1479. Klein, P., L. Richard and H. James, 2008, 3D curvature attributes: a new approach for seismic interpretation: First Break, 26, 105–111. November, 2012



55



References Marfurt, K. J., 2006, Robust estimates of 3D reflector dip: Geophysics, 71, 29-40 Marfurt, K. J., and R. L. Kirlin, 2000, 3D broadband estimates of reflector dip and amplitude: Geophysics, 65, 304–320. Marfurt, K. J., R. L. Kirlin, S. H. Farmer, and M. S. Bahorich, 1998, 3D seismic attributes using a running window semblance-based algorithm: Geophysics, 63, 1150–1165. Marfurt, K. J., V. Sudhakar, A. Gersztenkorn, K. D. Crawford, and S. E. Nissen, 1999, Coherency calculations in the presence of structural dip: Geophysics 64, 104–111. Roberts, A., 2001, Curvature attributes and their application to 3D interpreted horizons: First Break, 19, 85–100. Taner, M. T., F. Koehler, and R. E. Sheriff, 1979,Complex seismic trace analysis: Geophysics, 44, 1041–1063.



November, 2012



56