CA2244714C - Spectral decomposition for seismic interpretation - Google Patents

Abstract

The present invention is directed generally toward a method of processing seismic data so as to provide improved quantification and visualization of subtle seismic thin bed tuning effects and other sorts of lateral rock discontinuties A reflection from a thin bed has a characteristic expression in the frequency domain that is indicative of the thickness of the bed: the reflection has a periodic sequence of notches in its amplitude spectrum, said notches being spaced a distance apart that is inversely proportional to the temporal thickness of the thin bed. Further, this characteristic expression may be used to track thin bed reflections through a 3-D volume and estimate their thicknesses and lateral extent. The usefulness of this invention is enhanced by a novel method of frequency domain whitening that emphasizes the geologic information present within the spectrum. Although the present invention is preferentially applied to a 3-D seismic volume, it is alternatively applied to any collection of spatially related seismic traces.

Description

CA 02244714 2004-12-14 ..
' . . f . . _ , .,WO gg/25161 ~ ' PCT/US97I21195 SPECTRAL DECOMPOSITION FOR SEISMIC INTERPRETATION
FIELD OF THE INVENTION
The present invention is directed generally toward a method of quantifying and visualizing subtle seismic thin bed tuning effects. This method disclosed herein is implemented by decomposin_ the seismic reflection signal into its frequency components through the use of a discrete Fourier (or other orthonormal) transform of length dependent upon the thickness of the bed which is to be resolved. :after decomposition by said discrete transform, the coefficients obtained thereby are c~r2anized and displayed in a fashion which reveals and enhances the characteristic .frequency domain expression of thin bed reelection events.
making variations in subsurface layer thici;ness visible that otherwise might not have been observed. The presecit invention allows the ,s~ismir interpreter to analyze and map subsurface geologic arid > ' stratieraphic features as a t~unction oi~ spatial position. travel time, and frequency to an extent that has heretofore not been possible.
BACKGROUND
Bv most standards exploration Leophysics is a relatively young science. with some of the earlicst~work in th~° suh,ject area dating back to the 19?U"s and the rcncvvned C1~I1' approach datin~~ from only tlt~ l9:iU's. In the years since its genesis. however. it has become the oil Btu industry's preemincm approach to tindin~= subsurface petroleum deposits.
Although exploration geophyics Lenerally encompasses the three broad subject, areas of gravity, magnetics. and seismic. today it is the seismic niethod that dominatca almost to the pint of wcludin~_ the other disciplines. In ' fact, a simple count of the number of seismic crews in the field has become one accepted measure of the health of the entire oil industry A seismic survey ~ represents an attempt to map the subsurface of the earth by sending sound energy down into the ~_round and recording the "echoes'' that return from the rock layers below. The source of the down-going sound energy might come from e~cplosions or seismic vibrators on land. and air guns in marine environments. Durin_ a seismic survey. the energy source is moved across the land above the geologic structure of interest. Each time the source is CA 02244714 2004-12-14 ' detonated, it is recorded at a great many locations on the surface of the earth. Multiple explosionirecording combinations are then combined to create a near continuous profile of the subsurface that can extend for many miles. In a two-dimensional (?-D) seismic survey, the recording locations are generally laid out along a straight Iine. whereas in a three-dimensional s (3-D) survey the recording locations are distributed across the surface in a grid pattern. In simplest terms. a 2-D seismic line can be thought of as giving a cross sectional picture of the earth (avers as they exist directly beneath 'the recording locations. A 3-D
survey produces a data "cube" or volume that is. at least conceptually, a 3-D picture of the subsurface that lies beneath the survey area. Mote that it is possible to ewract individual 2-D line surveys from «rithin a 3-D
no data volume.
n seismic sutvew is composed of a vcn~ large number of individual seismic recordings or traces. In a mpical ?-D survey. there will usually be several tens of thousands of traces. whereas in a 3-D survey the number of individual traces may run into the multiple millions of traces. A
seismic trace is a digital recordinL of the sound energy reflecting back from inhomogeneities in the subsurface, a partial reflection occurring each time there is a chan_e in.
the acoustic impedance of the subsurface materials. The di~_ital samples are usually acquired at Q.t704 second t-1 millisecond) internals, although ? millisecond and l millisecond intervals are also eo.mmon.
Thus, each sample in a seismic trace is associated with a travel time. a two-way travel time in the case of reflected energy. Further. the surface position of every trace in a seismic survey is =o carefully recorded and is generally made a part of the trace itself tas pan of the trace header informationl. '(his allows the seismic information contained v~ithin the traces to be later correlated v<-ith specific subsurface locations. thereby providing a means for posting and contouring seismic data, and attributes extracted therefrom. on a map (i.e..
"mapping~~). The signs! that is sent down into the earth is called a seismic waveform or seismic wavelet. Different ~5 ~ seismic waveforms are generated depending on whether the source is an air gun, dynamite or vibrator. The .term "source signature" or "source pulse" is generally used to describe the recorded seismic character of a particular seismic waveform.

. ~~ CA 02244714 2004-12-14 PCT/US97/2119s . , , wo 9srislsl .;_ A seismic source generated at the surface of the earth immediately begins to-travel outward and downward from its point of origin. thereafter encountering and passing through rock units in the subsurface. At each interface between two different rock units, there is the .
potential for a seismic reflection to occur. The amount of seismic energy that is reflected at an interface is dependent upon the acoustic impedance contrast between the units and the reflection coefficiem is one conventional measure of that contrast. The reflection coeff cient can be thought of as the ratio of the amplitude of the reflected wave compared with the amplitude of the incident wave. In terms of rock properties:
to reflection coefficient acoustic impedance below - acoustic impedance above-acoustic impedance below r acoustic impedance ak~.ove-=
P~~~~ - P~V~
PTV= t P~V~
where. the acoustic impedance of a rock unit is defined to be the mathematical product of the rock densiy (pl and p~ being the densities of the upper .ioWer rock units.
respcctivelyl multiplied ~timcs the velocity in the same rock unit. V ~ and V~ corresponding to the upper and lower rock unit velocities. t Strictly speakin~'. this equation is exactly correct only when the wavelet strikes the rock interface at vertical incidence. However.' it is ~!enerally accepted in the .o industy that the requirement of verticality is Satisfied if the wavelet strikes the interface within about ?0° of the vertical.) Reflected energy that is recorded at the surface can be represented conceptually as the convolution of the seismic wavelet with a subsurface reflectivity. function:
the so-called "convolutional model". In brief. the convolutional model attempts to explain the seismic signal recorded at the surface as the mathematical.. convolution of the down going source wavelet with a reflectivity function that represents the reflection coefficients at the interfaces between different rock layers in the subsurface. In terms of equations.
x(t) = w(t) * e(t) + n(t) .

WU 98/ZS161 ' ', ~ ~ PG"T/U897I21195 . . , -~-where x(t) is the recorded seismogram, w(t) is the seismic source wavelet.
e(t) is the earth's reflectivity .function. n(t) is random ambient noise. and "*~'' represents mathematical convolution. Additionally. the convolutional model requires, in part. that ( 1 ) the source wavelet remains im~ariant as it travels through the subsurface (i.e., that it is stationary and unchangingj, s and {?) the seismic trace recorded at the surface can be represented as the arithmetic sum of the separate convolutions of the source wavelet with each interface in the subsurface (the principle of "superposition'', i.e.. that wavelet reflectivity and propagation is a linear system.) Although few truly believe that the convoiutional model fully describes the mechanics of wave propagation. the model is sufficiently accurate for, most purposes: accurate enough to make the 1o model very useful in practice. The convolutional model is discussed in some detail in Chapter 3.3 ~f .Seismic Ua~a I'rocesvifyf by Ozdogan Yilmaz. Society of Exploration Geophysicists.
1987. ; ' Seismic data that have been properly acquired and processed can provide. a wealth of information to the expIorationist, one of the individuals within an oil company whose job it is to vs Iocate potential drilling= sites: For example. a seismic profile ~.:ives the explorationist a broad view of the subsurface structure of the rock (avers and .otien reveals important features associated vi~ith the entrapment and storage of hydrocarbons such as faults.
folds. anticlines.
unconform'iti.es. and sub-surface salt domes and reefs. among many others.
During the computrr processing= of seismic data. estimates of subsurface velocity are routinely generated and near ~t> surface inhomogcneities are detected and displayed. In some cases. seismic data can be used to directly estimate rock porosity. water saturation. and hydrocarbon content.
Less obviously, seismic waveform attributes such as phase, peak amplitude, peak-to-trough ratio. and a host of others. can often be empirically correlated with known hydrocarbon occurrences and that correlation applied to seismic data collected over new exploration targets, In brief. seismic data t5 pro~~ides some of the best subsurface structural and stratigraphic information that is available.
short of drilling a well.
r Seismic data are limited. through, in one crucial regard: rock units that are relatively CA 02244714 2004-12-14 ' . _~ l wo 9snsi6i ~ ' rcrivs9~mi9s _;_ "thin" are often not clearly resolved. In more particular. whereas seismic reflection data can provide a near "geologic. cross section" representation of the subsurface when the lithologic layers are relatively "thick". the seismic image that results when the layers are ''thin'' is much less clear. This phenomenon is known to those skilled in the art as the seismic resolution s problem.
Seismic resolution in the present context will be taken to refer to vertical resolution within a single seismic trace. and is loosely defined to be to the minimum separation between two seismic reflectors in the subsurface that can be recognized as separate interfaces - rather than as a single composite reflection - on the seismic record. By way of explanation. a ~ta subsurface unit is ideally recognized on a seismic section as a combination of two things: a distinct ret7ection orieinatin~.ae the for of the unit and a second distinct reelection. possibly -of opposite polarity. originatin!s tcom its base. In the ideal case. both the top and the bottom of the unit appear on the recorded seismo~=ram as distinct and isolated reflectors that can be individually "time picked' (i.e.. marked and identified) on the seismic section. the seismic data i ~ within the inten~al between the two time picks containin~~ information about the intervening rock unit. On the other hand. where the seismic unit is not sufficiently thick, the returning reflections from the top and the bottom of the unit overlap. thereby causing interference between the vyo veilectiotoevents and blurring the image of the subsurface. This blurred image of the subsurface is one evamplc of,the phenomena kllowtl t0 those skilled in the art as the ''thin bed" problem.
.o rigure 1 illustrates in a very Lencral way how the thin bed probiern arises under, the axioms of the convolutional model. Consider first the "thick" bed reelection depicted in Figure la. On the left side of this figure is represented a source wavelet which has been generated on the surface of the earth. The source wavelet travels downward unchanged through the earth alone path P1 until it encounters the rock unit interface labeled "A.'' (Note that the wave paths in this figure are actually vertical. but have been illustrated as angled for purposes of clarity.
This is in keeping with the general practice in the industry.) Iri Figure la.
when the down-going seismic vyaveform encounters Interface "A'' a portion of its energy is reflected back toward the W'O 98125161 PCT/US97/Z1195 surface along path P2 and is recorded on the surface as the reflected event Rl . Note that wavelet R 1 is reversed in polarity as compared with the source wavelet, thereby indicating a negative reflection coefficient at the "A" interface. This polarity reversal is offered by way of example only and those skilled in the an recognize that reflection coefficients of either polarity are possible.
The remainder of the down-Going energy (after the partial reflection at the interface ''A") continues through the thick bed until it strikes Interface "B'' at the base of the thick lithographic unit. Upon reaching the "B" interface, pan of the wavelet energy continues deeper into the earth along path P~. while the remainder of its enere~~ i,s reelected back to the surface along path to P4 where it is recorded as reelection R?. dote that the reflection from interface "B'~ occurs at'a later point in time than the reflection from interface ":~". The exact time separation between the two events depends on the thickness and velocity of the intervening liver between the two interfaces. with thick layers and / or slow velocities creatin~~ a ~~reater time separation between the top and base reelections. The temporal thickness of this layer is the time that is required for the seismic wavelbrm to traverse it.
Uc~ the surface. the composite thick hed reelection - the expression actually recorded -is the arithmetic .S'll)11 superposition) of the two returning reelections.
takin~_ into account the time separation between the events. Because the two returnine wavelets do not overlap in time.
the recorded seismic record clearly displays both events as indicative of two discrete horizons.
(Note that the time separation between the two reflected events has llUt been accurately portrayed in this figure. Those skilled in the art know that the time separation should actually be twice the temporal thickness of the layer.) Turning now to Figure I b. wherein a thin bed reflection is illustrated. once again a source wavelet is generated on the surface of the earth which then travels along path P6 until it encounters the rock unit interface labeled "C.'' (As before. the wave paths in the figure are actually vertical.) .~s is illustrated in Figure lb, when the down-going seismic wavefotm encounters Interface "C" a portion of its energy is reflected back toward the surface along path ' , , , . , WO 98!25161 ~ I PCT/US97/21195 _7_ P'7. where it i~ recorded as reflection R3. The remainder of the down-going energy continues through the thin bed until it strikes Interface "D". Upon reaching the "D"
interface. pan of the wavelet energy continues deeper into the earth along path P 10, while the remainder of its energy is reflected back to the surface along path P9 where it is recorded as reflection R4.
Note once again, that the reelection from interface "D" occurs at a later time than the reflection from interface "C ". however the temporal separation between the two reflections in the case of a thin bed is less because the distance the waveform must travel before being reflected from interface "D" is less. In fact. the time separation between the two reflections is so small that the returning (upward Going) wavelets overlap. Since the composite thin bed lU reflection is once again the arithmetic sum of the two returning reflections, the actual-record~eti sisnal is an event that is not clearly representative oC either the reflection from the top or the b~ a ., , .
of the unit and its interpretation is correspondinLly difficult. This indefinite composite reflected event exempli ,fees chc typical thin bed problem.
Needless to say. the thickness of a subsuri'ace exploration target is of considerable economic imponancc to the oil conipany explorationist because. other things bein~~ equal. the thicker the lithographic unit the greater the volume of hydrocarbons -it might potentially contain.
Given the importance of~ accurately determining layer thickness. it should come as no surprise thai a variem of approaches have heen employed in an effort to ameliorate the thin bed problem.
:1 first technique that is almost universally applied is shortenin~_ the length of the seismic wavclet, longer «~avclets gcncraliv offering worse resolution than shorter ones. During the data processing phase the recorder! seismic waveform may often be shortened dramatically by the application of well known signal processing techniques. By way of example. it is well known to those skilled in the art that conventional predictive deconvolution can be used to whiten the ?5 spectrum of the vayelet. thereby decreasing its effective length.
Similarly, general wavelet processing techniques. including source signature deconvolution and any number of other approaches. might altematiyely be appiied to attempt to reach a similar end result: a more compact waveform. Although any of these processes could result in dramatic changes to the 'gyp 9g/25161 ~ PCTlU597/21195 _g_ character of the seismic section and may shorten the length of the wavelet significantly, it is often the case that further steps' must be taken.
Even the best signal processing efforts ultimately only postpone the inevitable: no matter how compact the wavelet is. there will be rock layers of economic interest that are too s thin for that wavelet to resolve properly. Thus. other broad approaches have been utilized that are aimed more toward analysis of the character of the composite reflection.
These approaches are based generally on the observation that. even when there is only a single composite reflection and the thickness of the layer cannot be directly observed. there is still information to be found within the recorded seismic data that may indirectly be used to estimate the actual thickness of io the lithographic unit.
By way of example. Figure da illustrates a familiar "pinch out" seismic model.
wherein the stratigraphic unit of interest (here with its thickness measured in travel time rather than in length) gradually decreases in thickness until it disappears ( i.e.. "pinches out" ) at the left most end of the display. Figure 4b is a collection of mathematically generated synthetic seismograms ) s computed from this model that illustrate the noise tree convolution of a seismic vyavelet with the interfaces that bound. this layer. Notice that at the right most edge of rigure -1b. the composite sicnal recorded on the first trace shows that the reflector is cleary delimited by a negative reflection at the top of the unit and a positive reflection at its base.
Moyin<_ now to the left within Figure fib. the individual rcf7ections at the top and base begin to merge into a single composite rrilection and eventually disappear as the thickness of the interval goes to zero. Note.
however, that t'~ie composite reflection still continues to chance in character even after the event has degenerated into, a single reflection. Thus. even though there is little direct visual evidence that the reflection arises from two interfaces. the changes the reflections exhibit as the thickness decreases suggest that there is information contained in these reflection that might, in the proper 3s circumstances, provide some information related to the thickness of the thin bed.
The pioneering work oi~ Widess in 1973 (Widess. Hov thin i.r a tlzir~ bed?.
Geophysics.
Vol. 38, p. 1176-1180) has given birth to one popular approach to thin bed analysis wherein ,wo 9sns><61 PC'Tl~TS97n1195 calibration curves are developed that 'rely on the peak-to-trough amplitude of the composite reflected thin bed event, together with the peak-to-trough time separation. to provide an estimate of the approximate thickness of the ''thin'' layer. (See also. Neidell and Poggiagliomi, Strarigraphic Modeling and Interpretation - Geophysical principles aJZd Techniques, in SEISMIC STRATIGRAPHY APPLICATIONS TO HYDROCARBON EXPLORATION, A.A.P.G. Memoir 26.
1977). A necessary step in the calibration process is to establish a "tuning"
amplitude for the thin bed event in question. said tunin~,_ amplitude occurrin~~ at the layer thickness at which maximum constructive interference occurs between the reflections from the top and base of the unit. In theory at least. the tuning thickness depends only on the dotriinant wavelength of the wavelet, i., and is equal to i. '_' where the'reflection coefficients on the top and bottom of the unit are the same sib=n. and i../~I where the ret'lection coefficients are opposite in sign. ;~
f3ocausc of the flexibifitv of calibration-ype approaches. ti-iev have been used with~some success in rather diverse exploration settings. However. these amplitude and time based calibration methods are heavily dependent on careful seismic processing to establish the correct wavelet phas,c and to control the relative trace-to-trace seismic trace amplitudes. l~hose skilled in the sciYmic processin;,_ arts know. however. how difficult it can be to produce a seismic section that truly maintains relative amplitudes thrauQhout. Further, the calibration based method disclosed above is not well suited for examining thin bed responses over a tare 3-D
survey: the method works best when it can be applied to an isolated reflector on a single seismic ~o line. It is a difficult enou~=h task to develop the calibration curve for a singly line: it is much more difficult io find a calibration curve that will work reliahlv throughout an entire 3-D grid of seismic data-.
Heretofore. as is well known in the seismic processing and seismic interpretation arts.
there has been a need for a method extracting useful thin bed information from conventionally . acquired seismic data which does suffer from the above problems. Further.
the method should also preferably provide attributes For subsequent seismic stratigraphic and structural analysis.
Accordingly, it should now be recognized. as was recognized btv the present inventors. that there exists. and has existed for some time.~a very real need for a method of seismic data processing WO 98/251b1 . PCT/US97/21195 .gyp.
that would address and solve the above-described problems.
Before proceeding to a description of the present invention. however. it should be noted and remembered that the description of the invention which follows, together with the accompanying drawings, should not be construed as limiting the invention to the examples (or preferred embodiments) shown and described. This is so because those skilled in the art to which the invention pertains will be able to devise other forms of this invention within the ambit of the appended claims. Finally. although the invention disclosed herein may be illustrated by reference .to various aspects of the convolutional model. the methods taught below do not rely on any particular model of the recorded seismic trace and work equally well under hroad deviations n from the standard convolutional model.. . .

The present inventors have discovered a novel means of utilizing the discrete Fourier transform tc~ imas~e and map the event of thin beds and other lateral rock discontinuities in .
conventioi;al ?-U and 3-I7 seismic data. In more particular. the invention disclosed herein is motivated by the observation that the reflection ttom a thin bed has a characteristic expression in .o the frequency domain that is indicative of the thickness of the bed: . a homo'eneous thin bed introduces a periodic sequence of notches into the amplitude spectrum of the composite rellecticm. said notches bcin~= spaced a distance apart that is inversely .
prt~ponional to the temporal thickness of the thin 'bed. rurther, if the Fourier transform coefficients are properly displayed This characteristic expression may he exploited by the interpreter to track thin bed reflections through a 3-D volume and estimate their thicknesses and extent to a degree not heretofore possible. More generally, the method disclosed herein may be applied to detect and identify vertical and lateral discontinuities in the local rock mass. Also.
the usefulness of the present invention is enhanced by a novel method of frequency domain whitening that emphasizes the geologic information present within the spectrum. Finally, the instant invention 30 is also directed toward uncovering seismic attributes that can be correlated with subsurface w CA 02244714 2004-12-14 PCT/US9T/2i195 wo 9snslsi . -i 1-structural and.stratieraphic features of interest. thereby providing quantitative values that can be mapped by the explorationist and used to predict subsurface hydrocarbon or other mineral accumulations.
By way of general background. the present invention preferably utilizes a relatively short discrete Fourier transform to determine the frequency components of a seismic trace. As is known to those skilled in the art. calculation of a Fouriec transform of a time series, even one that is exclusively real valued. results in a collection of Fourier transform coeff cients that are complex.data values of the form ":~ + Bi". where "i" represents the "imaginary" number or the square root of a ne~ativ~e one. Further. it is yell know that the expression A
= Bi may t U equivalently written as:
~1+Bi=rci~
where.
r =~ A + fii ~= A-+ B-and 8=tan-'(~~.
The quarnimE3 is known as the phase. angle Ior just the "phase") of the complcv quantiy A + Bi, the quantity "r" its ma~~nitude, and the expression i.-~ - Bi~ is the mathematical notation for the magnitude of a comply valued quantiy: also called its absolute value. A
frequency spectrum is obtained from the l~ourier transform coefficients by calculating= the complex magnitude of each ~u transform coecient. Further, the numerical size of each coefficient in the frequency spectrum is proportional tc~ the stren~~th of that fz'equencv in the original data.
i~inally, after application of a Fourier transform to some .particular time series. the resulting series of complex coefficients are said to be in the frequency domain. whereas the untransformed data are referred to as being in the time domain.
~S
Returning now to a discussion of the instant invention. the invention disclosed herein relies on the general observation that a frequency spectrum calculated from a whole trace Fourier ~ .. ~
~V1>O 98125161 - ~ PG"T/US97/21195 transform tends to resemble the spectrum of the source wavelet. whereas shorter window spectra tend to reflect more of the underlying geological information. This is because long analysis windows encompass a great deal of geological variation. said variations tending over the long term to exhibit a "white" (or random and uncorrelated) reflectivity function.
which has a ''flat'' amplitude spectrum. Thus. the shape of a frequency spectrum calculated from an entire seismic trace is largely dependent on the frequency content of the source wavelet.
(See. for e;cample.
Chapter ?.?.1 of Seismic Data Prncess~ing by Ozdogan Yilmaz. Society of Exploration Geophysicists; 1987. the disclosure of which is incorporated herein by reference.) On the other hand. where the analysis windotv is so short that the earth reflectivity function is non-white. the resulting Fourier spectrum contains influences from both the wavelet and local geology. Over such small windows. ~eoloey acts as a tilter, attenuating the spectrum of the source wavelet and thereby creatin~~ non-stationary short-window spectra. ' The foregoing ideas are illustrated generally in l-figure 2. wherein a typical seismic trace and some frequency spectra calculated therefrom have been plotted. :1t the top of that fi~urc is i s the Fourier transform frequency spectntm of the entire seismic trace. The appearance of that spectrum ~is ~=enerallv that of a typical field wavelet. However. the spectra calculated over shorter windows. pictured at the bottom of Figure 2. are nonstationarv.
tending to reflect the underlying geology which may potentially chan'e dramatically over fairly short intervals.
The importance of this observation for the present invention is illustrated in Fi~_ure 3.
wherein two representative spectra are ~~enerically pictured. The upper frequency spectrum represents that of a typical broad bandwidth source wavelet. The, lower spectrum: however.
represents in a general way the frequency domain expression of a composite reflection from a thin bed. In this later case. the geology of the thin bed has tended to act as a frequency domain 25 filter and has produced its own mark on the frequency content of the reflected wavelet. As is generally illustrated in this figure. the present inventors have discovered that a homogeneous thin bed affects the amplitude spectrum of the reflection event by introducing ''notches". or narrow bands of attenuated frequencies. into it. thereby producing~a characteristic appearance. A

' CA 02244714 2004-12-14 ~. -' Wt? 98115161, ~ ~ PCTlUS97I21195 -!3-homogeneous-bed is one that has a constant velocity and density throughout its extent..~.Further.
the distance between the notches so intraduced is equal to the inverse of the "temporal thickness" of the thin bed layer, temporal thickness being the length of time that it takes for a wavelet to traverse the layer in one direction (the thickness of the layer divided by the its velocity). Thus, attenuated frequencies in the amplitude spectrum may be used to identify a thin bed reflection and to gauge its thickness.
Turning now to FiLUre :~. the results suggested in the previous paragraph are extended to the analysis of a simplistic '_-D geological model. wherein the frequency domain e~cpression of a thin bed is investigated: In Fi_ure Via. a typical "pinch out" reflectivity function geological n model) is presented. Fi~:ure Jc contains a grey-scale display of the Fourier transform frequetT~:y spectrum amplitudes calculated from this model. This display was produced by ereati:n!-~aa discrete time series at fifth' equally- spaced positions across the model.
each of which has oniv two non-zero values: one corresponding to the reflection caefticient at the top of the iaver. and the other to the reflectian coefficient at its base. Each of the non-zero values is placed within the otherwise zero-filled trace in positions that reflect the location in time of the top and bottom of the reflector. resyctivcly. :1 standard discrete Fourier transform was then calculated for the time series followed by caiculation oi'thc complex ma~_nitude of each coefficient.
In Figure -te, the lighter portions correspond to LarLer values of the amplitude spectra.
whereas the darker portions represent small a'alues. Thus. "notches" in the amplitude spectra =a are represented bsr the darker values in the plot. 'l~his ligurc displays.
in a most lateral sense. the Fourier transfa- w.t~ the geology and. more particularly, the characteristic siLnaturc that is ci :-impressed on the mavelet by this event. V4'hat is most significant about this plot relative to the instant invention is that. as the thickness of the model decreases, the spacing between the notches increases. Further, for a given model thickness the notches are periodic. with a period 35 equal to the reciprocal of the temporal thickness of the layer. Thus. if this signature - periodic frequency domain notches - can be located within a seismic survey. it is strong evidence that a thin bed is present.

' CA 02244714 2004-12-14 . . ~ wo 9sns~lsa PCT/US97/2119s .
x-.1 ~t.
According to one aspect of the present invention there has been pro~~ided a system for .
interpreting seismic data containing thin bed events. wherein the data are decomposed into a series of Fourier transform 2-D lines or 3-D volumes. thereby providing enhanced imaging of the extent of said thin bed layers. The instant embodiment utilizes a single Fourier transform s window which is separately applied to the portion of each seismic trace that intersects a zone of interesi.
Those skilled in the art will realize that the same method could also be practiced to advantage on a 2-D collection of seismic traces to enhance the visibility of thin bed reflections contained therein.
As a first. step a collection of spatialU related seismic traces are assembled. These tra~~s might be. for purposes of illustration only. one or more shot records. a constant offset gathei. a CMP gadter: a VSP survey. a w~o-dimensional seismic line. a m~o-dimensional stacked seismic line extracted from a 3-D seismic survey or. preferably.. a 3-D portion of a 3-D seismic survey.
Furthermore. the .present invention mieht also be applied to a~ ?-D or 3-D
sutwev wherein the n a data have been transposed. i.e.. where an "offset' or spatial axis ("X" or "~"~ axis for 3-D data) has been oriented so as to replace the vertical or "tirrtc~ axis. ~~ore ~_enerally, am~ ~-I~ volume of digital data may be processed by the methods disclosed herein. That bein~
said. for purposes of clarity the .vertical axis will be referred to as the "time" axis hereinatier. although those skilled in the art understand that the digital samples might not be separated by units of time.
~u Whatever thr choice. the invention disclosed herein is most effective when applied to a croup of seismic traco~ ~ftat have an underlyin~_ spatial retationship with respect to some subsurface geological feature. :\gain for purposes or illustration only. the discussion that follows will be couched in terms of traces contained within a stacked 3-D survey. although any assembled group of spatially related seismic traces could conceivably be. used.
A zone of interest is next selected within a particular 3-D volume. The zone of interest might be, by way of example, the undulating region bounded by two picked reflectors . In that case. the reflector is preferentially flattened or datumized (i.e.. made flat by shifting individual traces up or down in WO 98123161 ~ ~ PCT/US97/21195 ' -IS-time) before analysis. and possibly also palinspastically reconstructed. More conventionally, a specific bounded time inten~al (for example, from 2200 ms to 2404 ms) might be specified.
thereby defining a "cube" or, more accurately. a ''box" of seismic data within the 3-D volume:
a sub-volume. Additionally. the lateral extent of the zone of interest may be limited by speciying bounding "in-line" arid "cross-line" trace limits. Other methods of specifying the zone of interest are certainly possible and have been contemplated by the inventors.
The selection and extraction of the data corresponding to the zone of interest is known as "subsettinn'~ the datai. Ond criterion that Guides the selection of the zone of interest is the desire to keep the zone as short (in time) as possible. This is in 4ceepinn with the general t o philosophy espoused above red=arding the tendency of long-window Fourier transform spectra~~ ro resemble the wavelet and short-window Fourier transform spectra to contain snore geologically related information. '~i'ote that there is a "hidden" window enlargement that is olicn applied automaticaliv and unthinkingly to Fourier.transform windows: extension of the window size to a length that is a power of two. This lenethenine of the window is done for purposes of is computational efficiency. as mindow Ieneths that are pw~ers of two are candidates for the Fast FourierTransfcarm tFFTI al~_orithm. 1-lowever. the presrnt inventors specifically counsel against this industry-wide practice and prefer to use a more general. if less computationally efficient.
discrete Fourir:r transform algorithm. thereby keeping the length of the analysis window to its minimum possible value. Given the computational power of computers today.
there is little _e reason not to tr~tsform only the data within the zone of interest.

WO 98IZ5161 P.CT/US97I21195 vt Once the spectra have been calculated and stored. they are ready to be used in the geophysical exploration for thin beds. Note that it is critical that. when the data are subsequently displayed. each spectrum be organized and viewed in the same spatial relationship with the other spectra as the traces from which they were calculated. That is. spatial relationships that are present in the untransformed data must be preser<~ed in the arrangements of the transform coefficients. The presently preferred method of viewing the transform coefficients is to begin by .17_ forming them into a 3-D "volume" (tuning cube), provided. of course. that the input data were originally taken from a 3-D volume. Note. however, that the vertical ("z") axis is no longer "time'' as it was before the transformation. but rather now represents. by convention. units of frequency. as Fourier transform coefficients are stored therein.
The tuning, cube may now be viewed in any manner that could be used to view a conventional 3-D volume of seismic data.
That being said.
the present inventors have discovered that viewing successive horizontal slices through the volume of coefficients is a preferr:d way to locate and visualize thin bed effects. Note that in the tuning cube arrangement. a horizontal slice represents all of the coefficients corresponding to IU a single Fourier frequency and. thus is a constant frequency cross section.
Further, as an aid._u .
the analysis o.f the data contained within this volume, the inventors preferably animate a series ~~f horizontal views throw=h the volume. In the event that the zone ol' interest is a portion of .an individual seismic line rather thin a volume. the resulting display. beine a collection of spatially related seismic trace Fourier transform spectra displayed in their orieinal spatial relationship.
I; mill still be referred to as a tunin~~ cube herein. even though.
technically. it may not be a "cube"
of data. . .

WO 98125161 ~ PG"T/US97/Z1195 _18.
According to a second aspect of the present invention. there has been provided a system for processing seismic data to enhance the appearance on the seismic record of thin bed events, wherein the data are decomposed into a series of Fourier transform 2-D lines or 3-D volumes by using a series of overlapping short-window Fourier transforms, thereby providing enhanced imaging of thin bed layers.
As disclosed previously, the first step in the present embodiment involves the interpreter mapping the temporal bounds of the seismic zone of interest. As was described previously. the mapping: may result in a seismic data cube-or rectangular piece of an individual seismic line.
In the present embodiment. rather than applyin~~ a sin=l~~-window' Fouricr transform to each trace. instead a series of ovcrlappin~~ short window Fourier transforms arc utilized. The len~~th of the window and the amount of overlap varies with the particular application. but once U again the window length need not be equal to a power ot'''?" but rather should be chosen so as to host inr'zs:c the underlvin~ ~=eolow. Note that a wei!_ht might optionally be applied to the data within each short window before transformation and. as before. a Gaussian weight is the preferred choice.
As ~cach short-w-ind~w Fourier transform is calculated, the .o coefticients resultin~_ therefrom are separately stored within an individual tuning cube that remains asso~$a~ed vrith the short-window that gave rise to it. TvTote that in the instant case there will be as many tuning= cubes produced as there were overlapping ~rindows in the analysis.
Scaling. if it is applied. is applied separately to each frequency plane in each tunin~T cube.

2; Each short-window tuning cube produced by a sliding window may now be individually examined in exactly the same fashion as that proposed previously for first embodiment. Once again. each cube is preferably viewed in horizontal slices or constant frequency images. thereby providing a means for visualizing geological changes with frequency. Further.
since there is ' CA 02244714 2004-12-14 WO 98125161 ~ PCT/US97/Z1195 now a collection of tuning cubes calculated at different time points in the trace. in effect a collection of tuning cubes have been produced that span a range of depths in the subsurface Finally, according to a third aspect of Ehe present invention there has been provided a system for processing seismic data to enhance the appearance on the seismic record of thin bed events, wherein the data are decomposed into a series of Fourier transform 2-D
lines or 3-D
volumes by using a short-window Fourier transform and are then reorganized into single frequency tuning cubes. thereby prow idin' enhanced imaging of thin bed layers..
The first steps in this present embodiment mirror those of the previous two embodiments: the data are first interpreted then subsetted. Thereafter.
to a series of overlapping shoe windo4v. Fourier transforms are calculated from the seismic data wlithin the zone of interest. optionally preceded, by the.application of a weight or taper wi~tin each window betor~ eolculatin<r the transform. ;1s in the previous embodiment.
the coefficients .
Irom each short window transform are accumulated. In the instant case however.
rather than viewing the calculated Fourier transform coefficients as tuning cubes. the data are reorganized 15 into single frequency ener<_y cubes which can thereafter be examined either in a horizontal or vertical plane for evidence of thin bed effects.
In more particular. the reorganization contemplated by the present im~entors conceptually involves extractin~_ tom all of the tunine cubes ever<~ horizontal slice that corresponds to a particular frequ~ncv. 'hhen. thesr individual same-frequency slices are "stacked" together. the .o topmost slice containinc: coefficients calculated from the topmost sliding windcaw. the next slice containinu coefficients calculated from the first sliding,window below the top, etc. Mote nhat.
after rror~~anization. the volume: of coefficients is or~:anizcd into units of "x-y" and time:. This is because the vertical axis. is ordered by the "time~~ of the sliding window that gave rise to a particular coefficient.
To utilize the information with the single frequency tuning cubes constructed by the previous step. a seismic interpreter would select a frequency and the seismic ~~olume corresponding thereto (e.g., he or she might select the coefficient voltune corresponding to l,Ohz.

WO 98125161 ~ ' '~ PCT/US97/21195 -?0-andlor the volume for l lhz. etc.). Each constant-frequency cube may be viewed in plan or horizontal view, or in any other manner, thereby providing a means for visualizing geological charities with lateral extent for a particular frequency.
It is important to note that. in all of the above described embodiments. the fact that the original untransformed traces were spatially related provides additional utility to the invention disclosed herein. In more particular, it is well known that short-window Fourier uansform coefficients are inherently quite noisy and have poorer frequency resolution in comparison with a longer window transform. One approach that the present inventors have used to improve the reliability. of the transformed values is to appl a Gaussian ~ weight function to the .pre-transformed data values. However.. another equally important step taken by the present inventors is to display the coefficients within a volume in the same spatial relationship as that'of the input data. Since tile traces so.displayed contain spatiafU correlzted information. displaying=
them new to each other allows the observe to visually "smooth out" the noise and perceive the underlvins_: coherent si~rial information.
Finally, although the present invention is discussed herein in terms of the discrete Fourier transform. in reality the Fourier transform is just one of any number of discrete time data transformations that could used in exactly the same fashl0(t. The general sups of ( 11 computing a short vt~indow transformation (? ) associating the resulting transform coefficients into a volume.
and (s) examining the volume for evidence of thin bed effects. could be accomplished with a 'o wide varicm of discrete data transformations other than the Fourier. If' the transformation is other than a Fourier, tuning volumes would be formed by 'roupinL to~~ether ctiefficienrs corresponding to the same basis function. Thus, when ~~ single frequency tuning cube ~~ is aced hereinafter. the inventors intend that phrase to apply. not .just to a tuning cube. composed of traditional Fourier transform.coefficients. but also to' any cube formed to same-basis function ?5 coefficients. : , 'I'hose skilled in the art will understand that a discrete Fourier transform is just one of many discrete linear unitary transformations that satisfy the following properties: (1) they are linear operators that are (2) exactly invertible. and (3) their basis functions form an orthonormal WO 98/25161 , PCT/US971Z1195 set. In terms of equations. if ~c(k). k = 1. L, represents a time series, and X(n) its"nth"
transformed value. n = 1. L, then the forward transform of the time series may be generally written for this class of transformations as L_t X(n) _ ~ x(k)A(k:n).
A=~I
where A(k:n) represents the forward transform kernel or collection of basis functions. Further, there is an inverse transform which maps the transformed values back into the ori'inal data values:
1 _, x(k)=~X(n)B(k:n).
lu where R(k:n1 is inverse transform ,kernel. The requirement of orthonormalim implies that the inner products bem~een nvo different basis functions must be equal to zero, and the magnitude of each basis function must br equal to unity. This requirement may be succinctly summarized by the followin~~ cduations:
~A~,~;n)A'Ik:n)=cal j-kJ
I L~
n:a i .:
~A(k:n)A'(k:m.)=oln-m?
where ~(Il- 111) - (~,~~", I I.n-m ' ~u and A*(k:n1 represents the complex conjugate of A(k;nJ. For the discrete Fourier transform. the basis functions corresponding to a forward transform of length L are conventionally chosen to be the set of cornplea eaponentials:
A(lc;n) - ~e_,~'~"~L~k = O,L-11, There are thus L basis functions (or basis vectors in this case), one basis function for each value CA 02244714 2004-12-14 . - , of "n"
n=-~; ..,0,...~-1.
To summarize: each transform coefficient. X(n), calculated from a data window corresponds to a particular basis function. and a tunin~T wiume is formed by collecting all of the transform coefficients corresponding to a particular zone of interest and storing those coefficients in an auxiliay storaLe area in the same spatial arrangement as the traces from which each window was computed.
By way of another specific example. those skilled in the art understand that a discrete Vv'alsh tr;mst'orm could br used in place of , the Fourier transform and the ~t'alsh coefficients 1o similar!v ~~rouped_ displayed. anti anayzeu. fn the manner disclosed above.
a Welsh transform ma~~ b~~ computed within an overlappin~~ series of sliding= windows and the coefficients resulting therefrom or~,anized and stored into tunin~_ cubes. Rather than the calculated transform coefficients repr~sentin~~ ti-equency. of course. these coefficients instead represent a similar quantim caflcd "sequence" by those skilled in the art. Thus. ,"single sequencv'v tunin~~ cubes I ~ may he formed from the Waish transform coefticients in a manner exactly anaio~ous to that used in the construction of Fouricr tuning cubes. These "same Crequencv (or more '_enerally, same-hasis function) tunin~_ cubes" will be referred to as a single orthonortnai basis function tuning=
cube. hereinafter.
Finally. althou~_h the discrete hourier transform iv a transform that is characterized by a set of onhonormal basis functions. application of a non-trivial -wei~_ht function to the basis functions prior to computation of a transformation destroys their orthonormality. Under conventional theory. a weight function that is applied within a windovr is viewed as being applied to the basis functions rather than the data. thereby preserving the integrity of the underlying data. However. basis functions that were orthogonal before application of the weight function will generally no longer be so thereafter. Thar being,said. in point of fact whether the weight function is applied to the data or to the basis functions. the end computational result after . wo 9srislsl transformation is exactly the same.
One means of avoiding the minor theoretical dilemma that arises when a weight function is used with a discrete orthonormal transform is to select ari orthonormal transform l weight combination which is not so affected. By way of example. the local cosine (and local sine) s transform is a discrete orthonormal transform wherein the weight function of choice is a smooth, specially designed taper that preserves the orthonormality of the basis functions at the expense of some loss in frequency resolution. Further. the underlying rationale of the local cosine i sine transform provides a natural theoretical bride to the field of general wavelet transforms.
io BRIEF DESCRITTION OF THE DRAWINGS-Figure 1 is a sch~a~~atic dia_ram that illustrates generally the nature of the thin bed . .
problem.
Fieure 2 displays a mpical seismic trace and compares long and short window spectra computed theret~om.
Figure 3 illustrates generally how the response of a seismic wavelet to a thin bed is expressed in the frequency domain.
Fi~:ure .t contains a simple seismic pinch out model. the convoiutional response thereto.
and the frequenc~~ domain representation ot'said convolutional response.
~u Fieurc ~ is a schematic diaLram that illustrates thu general approach of a presently preferred endiment.
Figure 6 contains a schematic illustration o!' how a presently preferred embodiment of the present invention is used in an exploration setting.
Figure 7 is a schematic illustration of another presently preferred embodiment.
Figure 8 is a slow chart that illustrates a presently preferred embodiment.
Figure 9 is a schematic illustration that describes the appearance of a thin bed during animation of constant frequency slices.
Figure 10 illustrates the general approach utilized to scale the constant frequency slices . WO 98125161 . pCTlUS97121195 _y_ so as to enhance the geologic content of the transformed data.
Figure 11 is a schematic illustration of another presently preferred embodiment.
DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
The instant invention provides a method of processing seismic data using a discrete Fourier transform. whereby its utility as a detector of thin beds is enhanced.
According to a first aspect of the present invention. there has been provided a method of enhancing and viewing thin bed effects using a discrete Fourier transform wherein a single Fourier transform is calculated for a window .spanning the zone of interest and the coefficients to obtained thereti'om are thereafter displayed in a novel manner. As is illustrated generally in Figure ~. let x(k.j.nt) represent a 3-D seismic data volume. where k = 1. l~.
and j = 1. J, represent indices that identify aspecific trace within a liven 3-I) volume. By way ot'~xample only. these indices might be in-line and cross-line position numbers, although other location schemes are also possible. The wriablc "nt" will be used to represent the time for dcpth~
position of within t5 each seismic trace. nt = 0. NTOT-1. the total number of samples in each individual trace. The time separation between successive values of x(k.j.nt) (i.e.. the sample rate) will be denoted by ,fit. where ,fit is customarily measured in milliseconds. Each trace in the 3-D volume. therefore.
contains a recordation of (NTOT)* ~t milliseconds of data. the first sample conventionally taken to occur at "zero" time. That bein~_ said, those skilled in the art understand that some seismic o data that arc eminently sui.tablc for analysis by the invention disclosed herein arc not ordered in terms of "time".~ By way of example only. seismic data samples that have been processed by a depth migration program are stored within a seismic trace in order of increasing depth. ,~z.
However. the instant invention can and has been applied in exactly the same fashion to this sort of data. Thus. in the text that follows 0t land "time") will be used in the broader sense of ?5 referring to the separation between successive digital samples, whatever measurement form that separation might take.
As a first step.'the explorationist or seismic interpreter selects a zone of interest within --,) CA 02244714 2004-12-14 . W0.98n51i1 ' PG"F/US97J21195 the 3-D volume. This might be done, by way of example only, by digitizing time picks ("picking") seismic events either on a digitizing table or. more commonly, at a~: seismic workstation. When an event is picked. the explorationist attempts to pinpoint the same reflector feature (e.g., peak, trough. zero crossing. etc.) on every seismic trace in which it appears. the ultimate coal being the production of a computer tile that contains time and surface location information that tracks the event across a 2-D section or through a 3-D
volume. As illustrated in Figure 11. given this information a computer program can be designed to read the picks and find the zone.of interest for any trace within the data volume. and / or perform the method of the -present invention. Said program might be transported into the computer. for example. by magnetic disk. tape. optical disk or CD-ROM.
Alternatively. the interpreter might simply specifyr constant starting and endir~ tires - , which bound the event of interest throughout the entire volume. thereby creating a "cube-w of interest. where "cube' is used in the Generic sense to represent a 3-D sub-volume of the original 3-D survey volume. For purposes of illustration only. the diseussion that follows will assume that a ~-D sub-cube has been extracted, although those skilled in the art will recognize that the same techniques discussed below can easily be adapted to a window that is not constant in time.
Again. just for purposes of illustrating the disclosed techniques. the temporal zone of interest.
after extraction. will be assumed to extend from the first sample of the 3-D
sub-volume. to last samplc,~sample number "''" hereinafter. Similarly. it wits also be assumed hereinafter that the ~u zone of intecl~.<i~ present on.evey trace in the sub-volume, althouLh those skilled in the art will recognize th'~:i~ often the case that the zone of interest extends to only a ponion of the 3-D
volume.
Given a zone of interest. the next step is to select a Fourier transform window length, 25 "L"hereinafter. Generally speaking, the length of the transform should be no longer than is absolutely necessary to encompass the zone of interest. Conventionally. the length of the Fourier transform is selected for purposes of computational efficiency and is usually restricted to an integer power of ? (e.~T.. 3?. 64. 128. -etc.). thereby allowing the highly efficient FFT calculation .~ CA 02244714 2-004-12-14 _ . wo 9snslsl ~ .

-?6-algorithm to be utilized. rather than a somewhat less efficient mixed radix Fourier transform or a much less efficient general discrete Fourier transform. However, within the context of the present invention the inventors specifically recommend against extending the chosen window length. as is conventionally done, to an integer power of two: a general discrete Fourier transform should be used instead. That being said. in the discussion that follows. it is understood by those skilled in the art that whenever a discrete Fourier transform is called for, an FFT will be calculated if appropriate. Otherwise. a general discrete Fourier transform, or some mixed radix variant, will be selected if the chosen window length is not an integer power, of 2.
Before beginning the Fourier transformations, an auxiliary storage volume must be I o established in which to store the calculated Fourier coefficients.
Auxiliary storage least as laxge as L computer words in event must be provided for each trace in which to save the calculated transform coefficients. with even more storage being required if the seismic data values or~the transformed results are to be kept as double (or higher)-precision. By way of explanation, a Fourier transform of a real time series of length L requires storage for Lf' complex data values.
each of which normally requires two computer words o.f storage. (There are actually only [(L/?) -I] uniqtu: complex data values. rather than L, because for a real time series the Hourier transform coeCfieients corresponding to positive and negative frequencies are directly related:
they arc complex conjugate pairs. In addition, there are two real. values: the coefficient at zero ("dc'-) hertz and the coefficient at the Nyquist frequency. both of which could be stored in a ''o single complex data value. Finally. if L is an odd inteter. the number of unique data values is (L+I1:'?). If th~tc: arc a total of J times h seismic traces within the zone (cubef of interest, the total amount of auxiliary stora~_e required will be equal to. at minimum. the product of L, J. and K measured in computer words. Let, the array A(k.j,nt) represent an auxiliary storage area for the present embodiment.
?s As a .first computational step. and as illustrated in Figure 8, the data values within the zone of interest are extracted from an input trace x(j,k.nt) taken from the sub-volume:
y(nl) = x(j,k.nl), nl = 0, L - 1 wo ~nsi6i rc"r~s9~nm9s and the weight function is optionally applied: w y(nl) = y(nl) w(nl), nl = 0. L -I , where the array y(nl) is a temporary storage area. (Note that in this present embodiment. the length of the analysis window is equal to the length of the zone of interest.) The weight function w(t). or data window as it is referred to by some. could take any number of forms' Some of the more popular data windows are the Hamming, Harming. Parzen. Bartlett. and Blackman windows. )each window function has certain advantages and disadvantages. The present inventors. however. have discovered that the use of a Gaussian window is in many ways optimal for this application. The Ciaussian weight function is defined by the following expression:
IU _y -1 n I - 1..' 21' i tr, w(nl)=a.~e , nl=O.L-1 where. " ' _L ~ , 1 6, - . ~, - a,, a~ -6 ~ 2na , In general. though. the weir=ht function should be a real function and non-zero over its ran~ee.
I ~ After tht: weight function has been applied. the discrete Fourier transform is then calculated according to~ the lbllowing standard e;cpression: ' X(n) - ~ v(k)e_,-'~~":~. ti = 1= . 0 .. L -1 , . . , , . ,~ , ,." _ _ where :~I n 1 represents the complex Fourier transform cocfficienn at the frequency, ftl, said ~u frequency being dependent on the length of the window L. In general. it is well known that the Fourier transform produces coefficients that provide estimates of the spectral amplitude at the following Fourier frequencies:
- ~(at i l ooo) ' n - - L .. . 0.... z - I.
It should be noted in passing that the nominal sample rate of a seismic trace.
0t, may not be the CA 02244714 2004-12-14 _ WO 98n5161 ) ~ ~ PCT/US97/21195 _28_ same sample rate at which the data were .acquired in the field. For example.
it is common practice to resample a seismic trace to a larger sample rate to save storage when there is little useful information at the highest recorded frequencies. On the other hand, on occasion a seismic trace may be resampled to a smaller sampling rate when. for example. it is to be combined with other -- higher sample rate -- lines. In either case, the nominal sample rate of the data may not accurately reflect its true spectral bandwidth. A simple modification of the previous equation will accomriuodate that contingency:
f n - n Fm ax ~ n= L, ..., ~, .. L-1, ., 2 . .
. __ where Fmax is the highest frequency contained in the data.
Since a seismic trace is a "real" function (i.e.. rion-ima~_inary), its Fourier transform is symmetric and the Fourier coefficients corresponding to the positive and negative frequencies are related as follows:
RE[X(f )J = RE[X(f ~)J.
o and _.
IM[X(f )] =-IM[X(f.~~)J.
where RL:[z] is a function that' extracts the real portion of thz complex value z. and 1M[z]
extracts the imaginary portion. As a consequence of this relationship. only L,'2 - 1 unique values arc produced within each Fourier transform window. Thus, for purposes of specificiy, 'o. only the positi~c frcquencics will be considered in the discussion that follows. although those skilled in the art understand that the same results could have been obtained by utilizin_ only the negative frequencies.
A nest step in the process involves placing the calculated complex frequency values into the auxiliary storage array. These traces are filled with the calculated complex Fourier coefficients as indicated below:
A(j,k,i) _ x(i), i = o. L/2.
wherein. "j" and "k" match the indices corresponding to the original data trace. In practice. the ' ~ . ' vwo 9srzsisi rc~rws9~nm9s _?9.
array A(j,k,i) .may not ever actually be kept entirely in R.AM (random access memory) at one time. but may be located, in whole or in part. on tape, disk. optical disk. or other storage means.
Additionally, because the presently preferred thin bed display requires the use of the frequency spectrum rather than the complex values. it would be convenient at the same time to calculate s the complex magnitude as each coefficient is placed into the auxiliary storage array:
A(j,k,i) _ ~X(i)~. i = 0, L/2.
Homever. there are many circumstances in which the complex coefficients would be needed and useful. so, as indicated in Figure 8, the complex coefficients are preferentially stored in the auxiliary storage area.
a.
The procedure described above is repeated for every trace in the deiined sub-voluti~e.
thereby filling the auxiliary storaLe array with transform coefficients in preparation forT'vicwing by the explorationist. Before viewing the results. however. the data are preferentially scaled in a novel fashion. whereby the '_=eological inFormation within the transform coefficients is emphasized relative to the contribution of the wavelet. This general. method involved in this frequencwdomain scaling is illustrated in Figure 1U. The scalin~~ method disclosed herein is designed to equalize the avers«e spectral amplitude in each ti-equencv slice.
thereby tending to produce a whitened wavelet spectrum. As ,illustrated in some detail in Figure 8. Let Tfj,k.ii represent a temporary storage array into which an entire tuning cube will_be stored. For a given qtr frequency slice. t. calculate the averate spectral amplitude therein:
YAV~i= -~ ~ IT(j.k.i,i~ .
~~ i.~
25 The spectral magnitude has been calculated because the T(j,k,i) are potentially complex valued.
As a next step. the values in this particular frequency slice are adjusted so that their average is equal to some user specified constant value, represented by the variable AVG:

_.

T(~'k'') TAVGTtJ,k,t)~ .l=t,J, k=1.K, where the primed notation has been used to indicate that the T(j,k,i) array has been modified. In practice. AVG will be set to some specific numerical value. 100. for example.
This scaling operation is repeated separately for every frequency slice (i = 0. L/2) in the tuning cube volume.
At the conclusion of this operation. every slice has the same average amplitude and a kind of spectral ~baIancine has been performed. Note that this form of single-frequency scaling.is gust one scaling algorithm that could be applied to the tuning cube data and the instant inventors have contemplated that other methods might also be used to advanta=e. By way of example. rather to than computing the arithmetic average of the items in the a slice, another measure of central tendency or any other statistic could have been equalized instead (e.'., median. mode. ~=eome~ic mean. variance. etc.). .rls another example. rather setting the averase value in each frequency slice equal to the same constant, each slice could be set equal to a different constant average value. thereby enhancing some frequencies in the spectrum arid suppressing others.
if the scaled ~tunin~; cube data are now inverted back into the time domain using a standard Pourier transform inverse. a spectrally balanced version of the original input seismic traces are thereby obtained. Let X(k) represent a scaled collection of transform coefficients obtained by the previously disclosed process and taken from location (j,k) within the scaled tuning cube array. Then. a spectrally whitened version of the input data may be obtained by means of the fallowin~, equation:
1 1 L-' +2~ik(nl) / L
x' ( j; k, nl) _ - ~ ~ .~ X(k)e ,~nl = O;~L -1 L, w(nl) k-0 ?s where x'(j,k;nl) represents the now modified (spectrally balanced) version of the input data x(j,k.nl). The divisor w(nl) is there to remove the effects of the' weight function that was applied prior to transformation. That term may be omitted if no weight w-as applied in the fot~~ard ' . ' WO 98J25161 . ~ PCT/US97/Z1195 .;1.
transform direction.
However, rather than inverting the scaled tuning cube, the presently preferred use for it is as an exploration tool for detecting thin beds. After all of the traces have been processed and placed in auxiliay storage, horizontal (constant frequency) amplitude slices.
Si(j,k), corresponding to the ''ith" frequency may be extracted from A(j,k,i) for viewing and / or animation:
siU.k) = IAG.k,i)I.
When these slices arc animated (i.e., viewed rapidly in succession) thin beds will be ~ ._.
recognizable as those events that successively alternate between high and low amplitude values.
t a Further: for many sons of thin beds. there will be a characteristic pattern of moving notches that clearly signal that an event is generated by a thin bed. Notc that it is preferable for the method disclosed herein that the slices be ordered in terms of frequency (either strictly increasing~or decreasing) when then are animated and viewed.
t: Figure 9~ illustrates the source of this diagnostic moving pattern. Figure 9a contains a lens-mpg ~eoloLic thin bed model and Figure 9b a sylized Fouricr transform of said model.
wherein only die notches have been drawn. As discussed previousy. the notches are periodic with period equal to the inverse of the temporal thickness of the model at that point. Now.
consider the model in Figure 9a as representinL a ?-D cross section of a 3-D
{disk-type) radially ~o symmetric m~e~l. and rigure 9b as a similarly radially symmetric collection of one dimensional Fouricr transf~s of said 3-D model. If the constant frequency plane labeled Plane 1 is passed through thewolumc as indicated. the plan view display of said plane twill reveal a low amplitude circular region corresponding to the first notch. Plane 2 passes through two notches and exhibits two fow amplitude circular regions. ~ Finally. Plane 3 contains three low amplitude circular ?5 regions. corresponding to the three notches, that it intersects. Now, when these slices are viewed in rapid succession in order of increasing frequency, there is a visual impression of a growing ''bulls eye" pattern wherein the rings move outward from the center. This pattern of moving notches is diagnostic for thin beds.

A~ CA 02244714 2004-12-14 __ When the thin bed is not circular, a related pattern is observed. Rather than concentric circles though, there will appear a series of moving notches that progress from away from the thicker areas and toward the thinner ones. For example. consider the model of Figure 9 as a cross section of a lens-shaped stream channel. When viewed in successive plan view frequency slices. a pattern of outward moving notches - moving from the center of the channel toward its periphery - will be observed all along its length.
It should be noted that if the thin bed is ,not homogeneous. for example if it contains a gradational velocity increase or decrease, it may not exhibit the characteristic "notch" pattern of the homogeneous thin bed. but rather have some different frequency domain expression. In 1 o these cases. the preferred method of identifying the characteristic response is to create a model of the event and calculate its Fourier transform. as was illustrated previously in Figure db.
Armed with this information. an explorationist may then examine an animated tuning cube~for instances of the predicted response.
Not only is the pattern of notches a qualitative indication of a homogeneous thin bed. but is ' it is also yields a quantitative measure of the e~ctent of the thin bed.
Returning to fiLUre 9. note that notches arc limited in lateral ewent~by the outer most edges of the model. . Thus, by panning through a stack of frequency slices and noting the outermost limits of movement by the notches.
a quantitative estimate of the event of the bed may be obtained.
'the foregoing is a striking visual effect that can bc: readily obsen~ed in actual seismic data volumes. Sine the typical non-thin bed event will have a somewhat consistent and slowly changing amplitude spectrum. the thin bed response is distinctive and easily identified. l~otc that in the present embodiment where a single window is calculated .for the entire zone of interest. the actual time position (i.e:. depth) of the thin bed within the zone of interest is not 35 ~ particularly important. If the thin bed is located anywhere within the temporal zone of interest.
the spectrum for that window will exhibit the characteristic moving notch pattern. Those skilled in the art will understand that moving the location of an e~rent in time does not change its amplitude spectrum. Rather. it only introduces a chance in the phase which will not be apparent wo 9snsi6i rcr~IS9znm~s if the amplitude spectrum is calculated and. viewed. , ....
As an alternative to displaying the amplitude spectra in animated plan view, the present embodiment may also be used with any number of other attributes calculated from the complex values stored in the tuning cube. By way of example. the phase of the complex transform Wcoefficients provides another means of identifying thin bed events and. more generally, lateral discontinuities in the rock mass. The phase tuning cube is calculated as follows:
,( IM(A(j,k.i))l P( j, k, i ) = tan ' RE(A( j, k. i))J .
where. P(j.k.i t contains the phase portion of the complex Fourier transform coefficients for evey point in the original tuning cube. Phase sections have long been used by those skilled in the, art to assist in picking indistinct reflectors, a phase section tending to emphasize continui~es in the seismic data. In the present embodiment however, it is lateral discontinuities in the spectral phase response that are indicative of lateral variability in the local rock mass. of which truncation of thin beds is a prime example. When viewed in animated plan view.
the phase values in the vicinim of a lateral edge will tend to be relatively"unstable":
tending to have an ill-behaved first derivative. Thus. the edges of thin beds and, more ~~eneralt, lateral discontinuities in the rock mass le.g.. faults. tiactures. non-conformities.
unconformities. etc.) will tend to have a phase that contrasts with surrounding phase values and ~wiil be. therefore.
relatively easy to identity. This behavior may used either by itself to identify lateral boundaries or in tandem with the amplitude spectrum tuning cube as a confirmation of the presence of local rock mass variabiliy.
Finally. it is anticipated by the instant inventors that the tuning cube technology disclosed herein might yield additional insights into seismic reflection data.
The tuning cube 35 (either containing phase or amplitude data) might be displayed and examined for empirical correlations with subsurface rock contents. rock properties. subsurface structure or layer stratigraphy. Alternatively. the Fourier transform values stored in the tuning cube may be further manipulated to generate new seismic attributes that can be useful in exploration settinES.

By way of example only, attributes that could be calculated from the tuning cube values include the average spectral magnitude or phase, and any number of other attributes.
The importance of this aspect of the present invention is best described as follows. It is.well known in the seismic interpretation arts that spatial variations in a seismic reflector's character may often be empirically correlated with changes in reservoir lithology or fluid content.
Since the precise physical mechanism which gives rise to this variation in reflection character may not be well understood. it is common practice for interpreters to calculate a variey of seismic attributes and then plot ar map them. looking for an attribute that has some predictive value. The attributes produced from the tuning cube calculations represent localized analyses of reflector properties n (beine calculated, as they are, from a shoe windowl and. as such. are potentially of consideral~fr importance to the advancement of the interpretation arts.
.Iii Figure S. the ''COVIPtJTL'' step. as applied to the present embodiment.
consists of at least one operation: calculating a discrete Fouri,er transform over the zone of interest. The resulting . .coefficients: the spectral decomposin tort of the zone of interest. are then stored as pan of the output spectral decomposition volume f "tuning cube") for subsequent viewing. Note that there will be one trace (i.c.: collection of Fourier transform coefficients) in the output tuning cube volume for, each seismic . trace processed as part of they input. Also note that in this.
presently .preferred output arrangement, horizontal plane slices, through the volume contain coefficients curiesponding to~a single common Fourier frequency..
.,o Optionally. -the "COMPUTE'' step may contain additional operations ~avhich . have the potential ~to enhance the quality of.the,~output volume and subsequent analysis., First. a weight function inay be applied to the seismic data within the zone of interest prior to calculating the transform. The 'purpose of the weight function is to taper or smooth the data within the Fouriei analysis window. thereby lessening the frequency-domain distortions that can arise with a "box-car" type analysis window. The use of a weight function prior to transformation ~is well known _5 , .
to those skilled in the art. 'The preferred wei~Thin~ function for the invention disclosed herein is Gaussian in shape and is in inanj~ ways optimal for this application. M'hat being said: note that rnam~ other weight functions could also potentially be used.

WO 98/25161 . PCT/US97~1195 Additionally. .since it is usually the amplitude spectrum that is of greatest interest to fee ..
explorationist. the amplitude spectrum may he calculated from the transform coefficients as then are moved into an auxiliary storage area. Alternativel. a phase spectrum. .or some other derived attribute. can be calculated from the transform coefficients before .storaLe and. indeed. these sorts of calculations have, been made by 'the present inventors.
Fiaallv. u, part of the computation step. :in indi~~idual frequency scalin~u may he applied to each plane (i.e.. frequency> in the output volume. :~s illustrated generalft~ iwFigure 10, the inventors lye found.it preferable to separately scale each frequency slice in the output volume to have the same average value before yievine it. This scalin~_~ is.just one-of many that might be applied. .but tl~e inventors prefer this method bocause it tends to cnhancc~
the Lrola~~ical content of the stored freqttcncv, spectra at~ihe cxpcrisc of the carrtmon wavelet information.
Animatin~~ successive horizontal slices through the spectral volume is the preferred method of vicwinL . .and analyzing ~thc_ transform coefficients. said..
animation preferably taking place un the computer moniu~r of a high speed workstation. As is well known to those skilled in, the art. aiiimati~n in the form of interactive panning through the volume is a fast and. efficient way to view large amounts of data. The data volume miglit be a~ievved ip horizontal. vertical. or, oblique 'slices. 'each of~vvhich provides a unique.view of the data: More importantly. however. in the context of the, present invention rapidly viewing successive horizontal, slices in succession provides . a; diagnostic means to survey a large volume of data and identify the ~ thin bed reflections therein, the details of which will be discussed' below°.
Note that it is preferable for the method disclosed hereiw that the slices be ordered in. terms of frequency (either strictly increasing or decreasing) when they are animated and i~iewed.
According to a second aspect of the present invention. there has been prow ided a method of enhancin~~thin bed effects using a discrete Fouriertransform wherein a series of sliding short-window ~Fouricr. transforms arc calculated aver ~a window spanning the. zone of interest .and thereafter displayed in a novel manner.This method is illustrated 'enerallj~
in Figure 6 and in more detail-~in 1"iLUre 8. Conceptuallr~. the nrcsent embodiment may, be thought of as producing . a series of tuning cubes of the sort disclosed previously. one tuning cube for-each Fourier transform w~itidorv position specified by the user. . ~ ..

WO 98125161 PCTNS97f21195 . _ . -36-(>ncc again .a(k j:n) represents a 3-I) scisyie data volume and . "L"- the Ieneth of the choscwslidin~:-window houriei transform. In thispresent cmbodiinent '.'l."
will generally be substantially shorter than the length bf the zone 'of interest. ~. As before..
the length of the Fourier transform window is to be selected. not on the basis of computational efficiency. but rather with the intent. of imaging particular classes of thin bed events in the subsurface. Bv way of example. a reasonable starting point for the transform length is one that 'is just long enough to span the "thickest" thin bed within the zone of interest. Note that it may be necessary to wincrease this minimum length in circumstances where, for instance. the waveforrn ,is, not particularly compact. In this ,later case. the minitimin window length might be increased by as .
several windows iriight possibly be applied to each individual trace. Let AM(j,k,i) represent the volume of collected of Fourier.transform coefficients taken.from. all~traces in the~zone of interest . for the ''M"th calculated window gositioti. Note that the amount bf storage that must be allocated to this array .has increased markedly. Now. the total amount of storage depends on the number of. sliding windows calculated for each trace. say NW. and must be at least as may words .
of storage.as .the product of NW L. J, and K: _ , Storage = (NW)(L)(J)(K).
As was mentioned previously, it is entirely possible that AM.(j,k:i) may never. be kept completely in RAM. but instead kept partially in RAM and the remainder on disk.
Using the array notation introduced above and again assuming that the Fourier transfoiniz of the v~eighted data is stored in X(i). the transform coefficients for the Mth window of trace (~ j) are stored in array location: . ~ . ~ .
AM(j.k.i) = X(i), : i = 0, L/2. . .
Once again. the individual frequency slices within the numerous tunin<,~ cubes .stored in ~a AM(j:k.il are 'preferably scaled by the procedure disclosed in Figure 8 prior to their examination 25 ~. for thin bed artifacts. In each case. the horizontal frequency slices arc individually scaled so that theiraveraee value is set to some particular constant. thereby~whitening the spectra.-wo 9sns~6i rcrrvs~~nm9s ' -37-After processing the seisrnic traces within the zone of interest. each tuning cube may be individually examined for evidence of thin bed effects. As . before. . thin bed effects may be identified in the amplitude. spectra by viewing a scrics of horizontal slices corresponding to different _ frequencies. Furihermorc. this may 'be 'separately done for the tuning cube S corresponding to each window position. thereby obtaining some general indication as to the Temporal and spatial extent of a particular thin~bed event.
The embodiment illustrated generally in FIG.6 as applied to 3-D seismic data, but those skilled in the art will realized that the same method could also be practiced to advantage on a 2-D seismic line to enhance the visibility of thin bed reflections contained therein.
~ o. ~ ~ According to a third .aspect of the present invention., there has been provided a method of enhancing thin bed effects using a discrete Fourier transform in. the manner described above for the second embodiment, but containing the additional step of organizing the Fourier transform w: coefficients into single-frequency volumes prior to display and analysis.
This method ~is illustrated.generally in Figure 7. As disclosed supra in connection with the second embodiment, the auxiliar~r storage array AM(j,k,i) will be filled with Fourier transform coefficients and will be preferably scaled.
Let F(j,k.m) represent a single-frequency volume extracted from A~(~,k,i).
Theie will be L/2 + I different volumes ((L+I)/2 values if L is odd), one for each Fourier frequency t produced by a transform of length ''L".. A. volume corresponding to the "i"th Fourier frequency 'o , is extracted' from AM(j,k.i)~ as follows:
F(j,k,m~ = Am(j,k,i), m = I. NW. j = 1, J, k = 1, K.
In effect.. the array F(j.k.m) may be viewed conceptually as being constructed by taking horizontal slices from each of nhe sliding window volumes and stacking them in order, of . .
increasing sltori Window counter. M.
?5 ' -3 8-The advantage_c~'this present data organization for purposes of.thin bed recognition is that it provides a means by which the location of the thin bed e~~ent in time and space may be determined. By way of explanation. as was indicated previously the temporal location of the thin bed within the zone of interest does not affect its response: thin beds near the top of the zone of interest and thin beds near the bottom produce the same characteristic amplitude spectra, This is advantageous from the standpoint of identifying thin beds, but it is a disadvantage; in terms of determinin_=their potential for hydrocarbon accumulation-higher bed elevations being generally preferable.
I-iowev~r. in the present embodiment the volume of same frcquc:ncy slices has "time' as its vertical ~axis~~- the variable M being n counter that roughly corresponds to distance. dawn the seismic tract. This organisation provides additional utility in that an approximate time duration of a thin bed event can be established.
For purposes of illustration, assume that a given thin bed event has a frequency domain notch as 10 hertz. Then, every short window Fourier transform that includes that bed will exhibit the same notch. If a 10 hertz volume of slices is examined. there will be a range of slices that contain the notch. Thus. by viewing successive slices in the. constant frequency volume: it is possible to localize in time the reflector of interest. More importantly, if it i5 known that a particular notch occurs 'at, say. 10 hertz. the 10 hertz tuning cube can be animated and viewed as an aid in determining the lateral extent of the thin bed: the limits of the notch as observed in this, frequency tuning cube defining the, terminus of the bed. _ ,. . . -39-In the previous discussion: the language has been expressed in terms of operations performed on conventional seismic data. But, it is understood by those skilled in the art that the invention lierein. described, could be applied advantageously in other subject matter areas; and used to locate other subsurface minerals besides hydrocarbons. By way of example only. the same approach described herein could be used to process andlor analyze mufti-component seismic data. shear wave data, triaeneto-telluric data. cross well survey data. full wavefortn sotFic logs. or model-based digital simulations of any of the foreeoin~. In short.
the process disclosed l herein can potentially be applied to am' single geophysical time series, hut it is preferably applied to a collection of spatially related time series containing thin bed events. Thus. in the text that follows those skilled in the arc will understand that "seismic trace" is used herein in a generic sense to apply to geophysical time series in general.
While the inventive device has been described and illustrated herein by reference to certain preferred embodiments in relation to .the drawings attached hereto.
~~arious changes and further modifications. apart from those shown or suggested hereiri. may be made therein by those skilled in the art. without departin. from, the spirit of the inventive concept, the scope of which is to be determi~ted by the followin~~ claims.

Claims (18)

40 WHAT IS CLAIMED IS:

1. A method for the exploration of hydrocarbons, comprising the steps of:
(a) accessing a set of spatially related seismic traces, said spatially related seismic traces containing digital samples, said digital samples being characterized by at least a time, a position, and an amplitude;
(b) selecting a part of said set of spatially related seismic traces to define a zone of interest;
(c) transforming at least a portion of said seismic traces within said zone of interest using a Fourier transformation, that is characterized by a plurality of orthonormal basis functions, and that is applied to a window containing said digital samples to produce a plurality of transform coefficients associated with said orthonormal basis functions;
(d) organizing said transform coefficients into a tuning cube;
(e) multiplying said transform coefficients by a scaling value to form a scaled tuning cube, said scaling value being determined by:
(i) selecting at least two transform coefficients corresponding to a same said basis function, (ii) calculating a complex magnitude of all transform coefficients so selected, (iii) calculating an average value from all transform coefficient magnitudes so calculated, and, (iv) calculating a scaling value from said average value; and, (f) displaying said scaled tuning cube.

2. A method for the exploration of hydrocarbons, comprising the steps of:
(a) obtaining a representation of a set of spatially related seismic traces distributed over a pre-determined volume of the earth, said seismic traces containing digital samples, said digital samples being characterized by at least a time, a position, and an amplitude;

(b) selecting a part of said volume and the spatially related seismic traces contained therein to define a zone of interest within said volume;
(c) defining a window within said zone of interest, said window having a starting sample number and encompassing digital samples;
(d) transforming at least a portion of said spatially related seismic traces within said zone of interest using a discrete orthonormal transformation that is characterized by a plurality of orthonormal basis functions, and that is applied to a window containing said digital samples to produce a plurality of transform coefficients associated with said orthonormal basis functions;
(e) organizing said transform coefficients into a tuning cube, said tuning cube and the transform coefficients therein being associated with said starting sample number;
(f) repeating steps (c) and (d) for a plurality of window definitions, thereby producing a plurality of tuning cubes; and, (g) displaying one or more of said plurality of tuning cubes.

3. A method according to claim 2, wherein step (g) comprises the steps of:
(i) selecting an orthonormal basis function;
(ii) selecting a tuning cube from said plurality of tuning cubes;
(iii) extracting from said selected tuning cube a plurality of the transform coefficients associated with said selected orthonormal basis function;
(iv) repeating steps (ii) and (iii) for at least one other selected tuning cube;
(v) organizing said extracted transform coefficients into a single frequency tuning cube; and, (vi) displaying said single frequency tuning cube.

4. A method according to claim 3, wherein said organization in step (v) includes ordering said extracted transform coefficients by said starting sample number associated therewith.

5. A method according to claim 2, wherein said discrete orthonormal transform is a Fourier transform.

6. A method according to claim 2, wherein a weight function is applied within said window containing digital samples prior to transformation by said discrete orthonormal transform.

7. A method according to claim 6, wherein said weight function is a Gaussian weight function.

8. A method according to claim 2, wherein the step (g) includes the further step of recording visually perceptible images representative of one or more of said tuning cubes on a generally flat medium.

9. A method according to claim 8, further including the step of using said visually perceptible images to identify subsurface structural and sedimentological features commonly associated with the entrapment and storage of hydrocarbons.

10. In the exploration for hydrocarbons, a seismic attribute map prepared by a process, said process comprising the steps of:
(a) accessing by means of a computer, a data set comprising seismic traces distributedover a pre-determined volume of the earth, said seismic traces containing digital samples, said digital samples being characterized by at least a time, a position, and an amplitude;
(b) selecting a plurality of spatially related seismic traces from said seismic trace data set;
(c) selecting a zone of interest within said selected plurality of spatially related seismic traces;
(d) transforming at least a portion of said spatially related seismic traces within said zone of interest using a discrete orthonormal transformation, said discrete orthonormal transformation producing transform coefficients from said spatially related seismic traces so transformed; and, (e) organizing said transform coefficients into a tuning cube;
(f) calculating a plurality of seismic attribute values from said transform coefficients organized into said tuning cube; and, (g) displaying said seismic attribute values at locations representative of said positions.

11. In the exploration for hydrocarbons, wherein seismic data comprising reflected seismic energy are recorded as a function of time over a pre-determined volume of the earth to produce a plurality of spatially related seismic traces, said spatially related seismic traces containing samples, said samples being characterized by at least a time, a position, and an amplitude, a map for the exploration of oil and gas produced by the process of claim 10, comprising:
(a) a generally flat medium for recording visually perceptible images thereon;
and, (b) at least one visually perceptible image on said generally flat medium, said visually perceptible image representative of said calculated seismic attribute values.

12. A method for the generation of seismic attributes for use in the exploration of hydrocarbons, comprising the steps of:
(a) obtaining a representation of a set of seismic traces distributed over a pre-determined volume of the earth, said seismic traces containing samples, said samples being characterized by at least a time, a position, and an amplitude;
(b) selecting a part of said volume and the seismic traces contained therein to define a zone of interest within said volume;
(c) transforming at least a portion of said seismic traces within said zone of interest using a discrete orthonormal transformation, said very short time discrete orthonormal transformation producing transform coefficients;
(d) organizing said transform coefficients into a tuning cube; and, (e) calculating, from said tuning cube, a plurality of seismic attributes.

13. A computer based method of filtering geophysical time series, comprising the steps of:
(a) obtaining a representation of a set of spatially related seismic traces distributed over a pre-determined volume of the earth, said spatially related seismic traces containing digital samples, said digital samples being characterized by at least a time, a position, and an amplitude;
(b) selecting a part of said volume and the spatially related seismic traces contained therein to define a zone of interest within said volume;
(c) transforming at least a portion of said spatially related seismic traces within said zone of interest using a discrete orthonormal transformation, said discrete orthonormal transformation producing transform coefficients from said spatially related seismic traces so transformed; and, (d) organizing said transform coefficients into a tuning cube;
(e) scaling said transform coefficients within said tuning cube; and, (f) inverting said tuning cube using a discrete orthonormal transformation inverse, thereby producing a filtered version of said transformed portion of said spatially related seismic traces.

14. A method according to claim 13, wherein the step of scaling said transform coefficients within said tuning cube comprises the steps of:
(i) selecting at least two transform coefficients corresponding to a same basis function;
(ii) calculating a complex magnitude of all transform coefficients so selected;
(iii) calculating a statistical value from all transform coefficient magnitudes so calculated;
(iv) calculating a scaling value from said statistical value; and, (v) applying said scaling value to a plurality of transform coefficients corresponding to said same basis function.

15. A method according to claim 14, wherein said statistical value is an arithmetic average of all transform coefficient magnitudes so calculated.

16. In a digital computer wherein seismic traces obtained over a pre-determined volume of the earth are read into memory, wherein a plurality of spatially related seismic traces are selected from said seismic traces, and wherein a zone of interest within said spatially related seismic traces has been defined, a digital computer programmed to perform a process comprising the steps of:
(a) transforming at least a portion of said spatially related seismic traces within said zone of interest using a discrete orthonormal transformation, said discrete orthonormal transformation producing transform coefficients from said spatially related seismic traces so transformed;
(b) organizing said transform coefficients into a tuning cube;
(c) scaling said transform coefficients within said tuning cube; and, (d) inverting said tuning cube using a discrete orthonormal transformation inverse, thereby producing a filtered version of said transformed portion of said spatially related seismic traces.

17. A device adapted for use by a digital computer wherein a plurality of computer instructions defining the process of claim 16 are encoded, said device being readable by said digital computer, and said computer instructions programming said computer to perform said process.

18. The device of claim 17, wherein said device is selected from the group consisting of a magnetic tape, a magnetic disk, an optical disk and a CD-ROM.