Non-stationary filter design and application for seismic signals

Abstract

A method of designing and applying filters for geophysical time series data comprising obtaining a plurality of spatially related geophysical time series and transforming the time series using a wavelet transformation. The wavelet transformation coefficients may be organized into a plurality of subband traces. The method includes modifying one or more transform coefficients within a plurality of the subband traces or within all but one of the subbands of traces and then inverting the subband traces using an inverse transform to produce a filtered version of the transformed portion of the geophysical time series. The method allows for design and analysis of non-stationary filters and filter parameters in untransformed data space. Non-stationary signals may be filtered in all or in windowed portions of the geophysical time series data.

Description

FIELD OF THE INVENTION

This invention relates to the field of geophysical prospecting and, more particularly to methods for acquiring and processing seismic data.

BACKGROUND OF THE INVENTION

In the oil and gas industry, geophysical seismic data analysis and processing techniques are commonly used to aid in the search for and evaluation of subterranean hydrocarbon deposits. Generally, a seismic energy source is used to generate a seismic signal that propagates into the earth and is at least partially reflected by subsurface seismic reflectors (i.e., interfaces between underground formations having different acoustic impedances). The reflections are recorded in a geophysical time series by seismic detectors located at or near the surface of the earth, in a body of water, or at known depths in boreholes, and the resulting seismic data may be processed to yield information relating to the location of the subsurface reflectors and the physical properties of the subsurface formations.

Seismic or acoustic energy may be from an explosive, implosive, swept-frequency (chirp) or random source. A geophysical time series recording of the acoustic reflection and refraction wavefronts that travel from the source to a receiver is used to produce a seismic field record. Variations in the travel times of reflection and refraction events in these field records indicate the position of reflection surfaces within the earth. The analysis and correlation of events in one or more field records in seismic data processing produces an acoustic image that demonstrates subsurface structure. The acoustic images may be used to aid the search for and exploitation of valuable mineral deposits.

These seismic data are processed to be useful for delineating features in the earth's subsurface. Geophysical time series data acquired in the field often contain unwanted energy as well as desired energy. Some of the unwanted energy may interfere with, overwhelm or otherwise mask desired seismic energy useful for subsurface investigations. There are many known methods for filtering or separating unwanted energy from seismic signals.

Prior art processing methods for eliminating unwanted energy include frequency filtering, windowed frequency filtering, and f-k filtering. Filtering may be applied to each trace individually or to multiple traces concurrently. Elimination of unwanted energy traditionally results from the application of filtering to entire traces regardless of whether the unwanted energy is localized in only part of the data trace.

Traditional filtering is usually implemented in frequency-time space using orthogonal basis functions having perfect localization in frequency but infinite extent in time. Use of the Fourier transform is well known in the art and the most common type of filtering for geophysical data. The underlying assumption of this type of filtering is that the seismic signals are stationary (e.g., do not have time dependent statistical characteristics). However, it is well known in the art that many types of seismic signals recorded in the field are non-stationary. It would be beneficial to filter geophysical time series data containing unwanted non-stationary energy without adversely impacting desired seismic signal.

An alternative to the Fourier transform is the continuous wavelet transform (CWT). In the continuous wavelet transform, a function ψ is used to create a family of wavelets ψ (at+b) where a and b are real numbers, a dilating the function ψ and b translating the function. The CWT turns a signal f(t) into a function with two variables (scale and time), which may be called c(a,b) such that c(a,b)=∫f(t)ψ(at+b)dt. The CWT allows for joint time-frequency localization ideal for non-stationary filtering. The discrete wavelet transform (DWT) is a special case of the CWT.

Wavelet methods may be designed in transform space rather than data space. However, as is evident from comparing FIG. 2A , data prior to transformation, with FIG. 2B , data from the same seismic record after transformation with wavelet methods, the wavelet transform alters coordinate relationships in the original geophysical time series data ensemble. This altered coordinate relationship is due to an inescapable limitation on joint time-frequency resolution imposed by the transform technique. Though the transformed data contains the required non-stationary time-frequency signal characteristics, signal components are difficult to recognize from their original space-time geometry. To address this dilemma, some practitioners have developed graphing techniques to approximately restore original data-space coordinate relationships to the transformed signals. The details of the technique depend on the particular choice of transform. Deighan and Watts (1997) give an example for a wavelet transform in FIG. 6 of their paper. (GEOPHYSICS, VOL. 62, NO. 6, NOVEMBER-DECEMBER 1997; P. 1896-1903) U.S. Pat. No. 5,781,502 to Becquey proposes a method to discriminate elliptical waves propagating in a geologic formation by a combined process that includes measuring and comparing ratios of components along several wave axes. The method calls for the use of multi-component seismic data. The method is directed to removing elliptical waves from seismic data after wavelet analysis is used to determine the surface or tube waves in the data. The detected signals are subtracted from the field records.

U.S Pat. No. 5,850,622 to Vassiliou et al. disclose a method to filter seismic data using very short-time Fourier transforms. The method utilizes a very great number of short overlapping Fourier transform windows, together with a Gaussian weight or taper function, to produce a plurality of near single-frequency “sub-band” traces for each seismic trace so analyzed. The transformed data may be used for selective removal of coherent noise events, analysis for seismic attributes related to subsurface features of interest, seismic trace creation by interpolation, and automatic identification and removal of noisy seismic traces.

SUMMARY OF THE INVENTION

A method of designing and applying filters for geophysical time series data is described which comprises obtaining a plurality of spatially related geophysical time series and transforming the time series using a wavelet transformation. The wavelet transformation coefficients may be organized into a plurality of subband traces. The method includes modifying one or more transform coefficients within a plurality of the subband traces or within all but one of the subbands of traces and then inverting the subband traces using an inverse transform to produce a filtered version of the transformed portion of the geophysical time series. The method allows for design and analysis of non-stationary filters and filter parameters in untransformed data space. Non-stationary signals may be filtered in all or in windowed portions of the geophysical time series data.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which:

FIGS. 1A and 1B illustrate a 3D common shot ensemble of geophysical time series seismic signals before and after non-stationary filtering by wavelet transform;

FIG. 2A illustrates a shot record containing non-stationary seismic energy.

FIG. 2B illustrates the data of FIG. 2A after a wavelet transform has been applied.

FIGS. 3 a through 3 j illustrates seismic data with varying frequency filters applied.

FIG. 4 illustrates filter design methodology in non-transformed data space.

FIG. 5 illustrates an embodiment of the wavelet filtering process.

FIG. 6 illustrates the application of filter parameterization in transform-data space.

FIG. 7 illustrates a flow chart of the present invention.

While the invention will be described in connection with its preferred embodiments, it will be understood that the invention is not limited thereto. On the contrary, it is intended to cover all alternatives, modifications, and equivalents that may be included within the spirit and scope of the invention, as defined by the appended claims.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The present invention is a method for designing and applying non-stationary filters using wavelet transforms that overcome the above-described deficiencies of current seismic processing practices. The filtering technique of the present invention provides for very efficient filter design and parameterization as well as full or partial suppression of non-stationary seismic energy in field records. The filters may be applied in a time-localized manner to further reduce impacting desired signal at other record times. Other advantages of the invention will be readily apparent to persons skilled in the art based on the following detailed description. To the extent that the following detailed description is specific to a particular embodiment or a particular use of the invention, this is intended to be illustrative and is not to be construed as limiting the scope of the invention.

This invention provides effective non-stationary filters for attenuating non-stationary and stationary noise components in geophysical time series data signals. Noise components often interfere with desired signal within time and frequency data windows that vary according to source and receiver (spatial) location. This interference makes the composite signal non-stationary with respect to time and frequency. Groups of such signals with variable spatial location thus form an ensemble that is also non-stationary with respect to space. Ground roll, airwaves, swell noise, and cable noises are some examples of the noise components typically present in seismic data. The parameterization of the present invention prescribes filters adapted to the non-stationary characteristics of such noise.

Many practitioners have applied non-stationary filters over the time coordinate of seismic signals. However, these filters lack an adequate process for determining time support of the filter within an ensemble of data. For example, the onset and termination times of the strong and unwanted airwave energy vary significantly over a seismic data record. Considering that typical seismic surveys contain hundreds of such ensembles, each one with different unwanted energy time windows, this is a serious shortcoming of previous filter parameterizations. A method that would parameterize analytic trajectories prescribing nominal onset and termination times to eliminate unwanted energy in a time or space localized fashion is highly desirable. Under the new parameterization presented here, analytic trajectories for nominal onset and ending times of noise components (and thus filters) for each subband and spatial coordinate, may be designed and applied for every ensemble of geophysical time series data.

Many mathematical transforms exist for implementing time-frequency data decompositions necessary for non-stationary filters. Common transformations include wavelet transforms, wave packet transforms, sliding-window Fourier transforms, and matching pursuit decompositions. Each of these has different strengths and weaknesses. They provide different combinations of time-frequency resolution, accuracy and stability, and computational cost.

The wavelet transform may be applied as a discrete wavelet transform (DWT). The DWT is a convenient way of using the wavelet transform with digital geophysical time series data. In a discrete wavelet transform, a wavelet is translated and dilated only by discrete values. For example, the dilation may be a power of 2 and be of the form ψ (2 j t+k), with j and k whole numbers. Orthogonal wavelets are a special case of discrete wavelets that give a representation without redundancy and lend themselves to fast algorithms.

Some prior art filtering algorithms, Fourier transform methods for example, mix transform-data coefficients along the space coordinate of geophysical time series ensembles to detect and attenuate noisy coefficients. These approaches assume the desired signal has greater spatial coherence than the noise component. Other approaches zero or mute, to a predetermined level, all coefficients within the time subband support of the filter. This practice can easily harm or delete desirable signal within that same time subband window.

FIG. 1A shows an example of non-stationary seismic energy, for instance airwaves 101 or ground roll, in a 3-D common-shot ensemble of geophysical time series data. The time-space distribution of the non-stationary energy component varies considerably, but predictably, throughout the ensemble. This noise concentrates, moreover, within particular subbands of the frequency spectrum. The present invention limits filter action within specified subband(s) to a computed time-space neighborhood of the undesirable signal present in geophysical time series data ensembles. Two analytic trajectories per subband define this neighborhood throughout the ensemble, wherever the noise component resides. The resulting filter passes data outside the prescribed time-space-subband window, unchanged and unharmed.

Desirable signal may also be present within that time-space-subband window. This parameterization employs a two-stage prescription for the filter action within the window of data to be filtered. This method of design biases filtering towards those data most likely to be noise-saturated.

This new filter design methodology presented here applies to invertible non-stationary transforms. Given an ensemble of signals, it conveys non-stationary characteristics in data space with the efficacy and cost afforded by that particular transform. One may therefore choose and evaluate filter design parameters directly in data space, incorporating the strengths and limitations imposed by one's choice of transform.

The new filter application parameterization of the present invention establishes a data-adaptive threshold with the assumption that only noisy data can exceed it. A complementary parameter defines a second data-dependent level. This two-stage filtering mechanism reduces to the second, lower level only those coefficients that exceed the first threshold, established for noise. This new filter parameterization applies to invertible non-stationary transforms. It prescribes a non-stationary filter consistent with the efficacy and cost afforded by that particular transform.

FIG. 1A compared to FIG. 1B gives a particular example of the application of a wavelet transform filter designed with the present invention. FIG. 1A contains non-stationary energy 101 , for example airwave or ground roll energy. After filtering particular subbands of the data using the wavelet transform, the strong non-stationary energy 101 of FIG. 1A has been removed 201 in FIG. 1B demonstrating non-stationary filtering as prescribed by this new parameterization and implemented via wavelet transform.

This invention processes ensembles of seismic signals to produce output that is optimal for non-stationary filter design. The processed output allows analysis of non-stationary time-frequency characteristics in the non-transformed data space, thus permitting filter specification and testing directly in non-transformed data space. Non-transformed data space naturally preserves space-time coordinate relationships required for effective filter design. Moreover, practitioners ultimately judge effectiveness of the filter from its impact in non-transformed data space.

Undesirable noise signal components often combine and interfere with the desired signal at variable times and with various frequencies. This interference results in non-stationary composite seismic signals. Non-stationary filters are therefore required to attenuate these undesirable components. Ground roll, airwaves, swell noise, and cable noises are examples of recorded seismic record components usually considered undesirable within the context of data processing objectives. FIG. 2A illustrates a shot record containing non-stationary seismic energy, ground roll and airwave components 101 .

This invention specifies the support of a non-stationary filter by a nominal beginning time t b and ending time t e as, for example, illustrated in FIG. 4 . Parameter functions t b (j) and t c (j) give a pair of times at every spatial coordinate x in the data ensemble, and apply for every subband j of the filter. These functions encompass families of two-parameter trajectories, including the linear and hyperbolic examples depicted in FIG. 4 . In FIG. 4 the beginning filter time t b is an example of a hyperbolic trajectory 401 and the ending filter time t e trajectory is an example of a linear function 403 . In these two examples, scalars τ b (j) and τ e (j) take on the interpretation of beginning and ending times, respectively, at spatial coordinate x=0. For these linear and hyperbolic forms, scalars p b (j) and p e (j) assume the role of reciprocal velocity.

The horizontal axis of FIG. 4 represents the space coordinate 415 . The vertical axis of FIG. 4 represents the time coordinate 413 . Any arbitrary space coordinate position x 417 will have an associated time t b 405 and a t e 407 representing the beginning and end of the filter support region 409 for that space coordinate. The data areas 411 outside of the filter support region 409 are not filtered.

FIG. 6 is a transform-data space representation of the filter application parameters. This figure is a generalization of how the parameters are used to apply the filtering. The vertical axis in FIG. 6 is the transform time with an example 16 transform data coefficients. ( 1 ) through ( 16 ). The FIG. 6 horizontal axis represents the transformed data coefficient values with positive values to the right of transform coefficient {circumflex over (d)} j =0 and negative values to the left. In FIG. 6 , the example filter parameterization acts on only three of these 16 transform coefficients ( 1 ) through ( 16 ) along the FIG. 6 {circumflex over (t)}transform time vertical axis. This filter parameterization treats the remaining thirteen as desirable signal, preserving them without change. ρ j prescribes the filter attenuation level as a fraction of the {overscore (s)} j data statistic. These two parameters act in concert with the filter time-support parameters to detect and attenuate noise components in a data-dependent fashion.

FIG. 6 graphs the roles of transform-data scalar parameters σ, and ρ j in this invention σ j specifies a noise detection threshold as a percentage of {overscore (s)} j , a data statistic computed from the transform coefficients in subband j of the ensemble. Reasonable choices for {overscore (s)} j include the mean, the median and variance of the data transform coefficients, and many others are possible. When transform data coefficients {circumflex over (d)} j in the filter region 619 exceed the threshold determined by −σ j |{overscore (s)}| 601 or +σ j |{overscore (s)} j | 607 , the muting function may reduce the data coefficient {circumflex over (d)} j to the values −ρ j |{overscore (s)} j | 603 or +ρ j |{overscore (s)} j | 605 . For example, data coefficient {circumflex over (d)} j ( 2 ) has a magnitude 621 greater than +σ j |{overscore (s)} j |. However, data coefficient {circumflex over (d)} j ( 2 ) is not within window filter region 619 , so no filtering of the data occur. Data coefficient {circumflex over (d)} j ( 7 ) also has a magnitude 613 which is greater than +σ j |{overscore (s)} j | 611 . Data coefficient {circumflex over (d)} j ( 7 ) is between the transform time coordinates t b at the beginning of the window filter region 619 window and the end of the window t e . The filtering operation will attenuate the magnitude of the data coefficient {circumflex over (d)} j ( 7 ) to +ρ j |{overscore (s)} j | 609 along the line represented by +ρ j |{overscore (s)} j | 605 . Similarly, data coefficient {circumflex over (d)} j ( 8 ) with transform coefficient 617 with a greater magnitude than −σ j |{overscore (s)} j | 601 would be attenuated to −ρ j |{overscore (s)} j | 603 . Data coefficient {circumflex over (d)} j ( 10 ) exceeds the threshold magnitude so the data coefficient will be attenuated to +ρ j |{overscore (s)} j | 605 . The other data coefficients within the filter region 619 do not exceed the threshold magnitudes and would not be changed.

Importantly, one can readily determine this suite of parameters for the seismic noise components enumerated above. One gleans this information from analysis of a few representative data ensembles, selected from the recorded survey. To attenuate the airwaves or other non-stationary seismic energy in FIG. 1A or FIG. 2A we may choose, for example, linear forms for the filter trajectories. That is, t b (j) (x)=τ b (j) +p b (j) x and t e (j) (x)=τ e (j) +p e (j) x, where x specifies the scalar distance from the source position to every receiver position in this ensemble. The remaining time-support parameters are determined as:

τ b (j) =0, τ e (j) =0.25, p b (j) =1/300, p e (j) =1/300. The airwave speed is about 300 m/s. The time support of the airwave is about 0.25 s. Further analysis shows this airwave to be concentrated in subbands j=5 and j=6 of the wavelet transform. For this example the data statistic, s 1 , may be chosen to be the mean of the transform data. The data coefficient threshold and attenuation parameters may, again for this example, be chosen respectively as σ 5 =σ 6 =10, ρ 5 =ρ 6 =1/0. This parameter suite produces the filtered ensemble in FIG. 1B where the non-stationary energy present in FIG. 1A has been greatly attenuated or eliminated.

The present invention enables non-stationary filter design in data space, the natural setting to observe space-time relationships. The present method restores space-time coordinate relationships to the data ensemble via the inverse transform, rather than through graphing in transform space. The problems inherent in designing filters in transform space may be seen by comparing FIG. 2B with FIG. 2 A. FIG. 2B contains the data in FIG. 2A after applying a wavelet transform. There is a limitation on joint time-frequency resolution that makes the signal components difficult to recognize from their original space-time geometry. Importantly, the present invention preserves time-frequency signal characteristics required for non-stationary filter specification, which facilitates convenient filter specification and design. The method utilizes a sequence of muting operators in transform space, each one followed by the inverse transform. FIGS. 3 a - 3 j are used to illustrate an example.

FIG. 3 a through FIG. 3 j is a common shot ensemble of filter panels from the seismic data of FIG. 2 A. FIG. 3 a is the same record as FIG. 2 A. FIG. 3 b is the wavelet-transformed data of FIG. 3 a (i.e. the same record as FIG. 2 B). FIGS. 3 c through 3 j are a progression of individual frequency bands inverse transformed after muting has been applied to wavelet transformed data of FIG. 2 A. Each figure panel is produced by selectively muting the wavelet-transformed data followed by inverting the transformed data.

In the FIG. 3 a - 3 j examples, the restored, linear geometry of the undesirable signal components primarily identifies them with frequency bands represented by FIG. 3 e and FIG. 3 f . The practitioner would therefore specify a non-stationary filter that attenuates subbands in FIG. 3 e and FIG. 3 f of the wavelet transform. Moreover, filter attenuation may be limited to a neighborhood of linear space-time coordinates occupied by the ground roll and airwave, such that no other data are filtered.

To use this invention for design of non-stationary filters for ensembles of seismic signals, operator T −1 M j T describes the entire filtering procedure in matrix notation. One applies non-stationary transform T to each input data vector in the ensemble, followed by filter operator M j and inverse transform T −1 . Operator M j modifies or otherwise filters all but the jth subband of the transform. FIGS. 3 c - 3 j demonstrate the operator for various subbands of a wavelet transform. For this transform the (i,l)th matrix element of M j is: m il (j) =δ ip δ pl where p=2 j +k and k indexes the time samples within subband j.

To use this invention for applying the designed parameterization, operator T −1 {tilde over (M)} j T likewise describes the entire filtering procedure in matrix notation. One applies non-stationary transform T to each input data vector in the ensemble, followed by filter operator {tilde over (M)} j and inverse transform T −1 . Operator {tilde over (M)} j modifies or otherwise filters one or more (jth) subband of the transform. The details of this matrix operator depend on the details or characteristics of the signal and noise in T For a wavelet transform, the (i,l)th element of {tilde over (M)} j is: {tilde over (m)} il (j) =δ il +δ ip ({tilde over (m)} pq −δ pq )δ ql . Scalars {tilde over (m)} pq hold the data-dependent attenuation weights, where p=q=2 j +k and k indexes the time samples within subband j.

A description of the elements of the non-stationary parameterization follows:

T −1 {tilde over (M)} j T: Filter operator, where {tilde over (M)} j is the filter, T is a non-stationary transform and j indexes the subband.

{overscore (s)} j : A statistic of the transform data (like mean, median, etc.; units are units of the transform data).

σ j : Unitless scalar parameter ≧0.

σ j |{overscore (s)} l |: Threshold for triggering filter.

ρ j : Unitless scalar parameter ≧0.

ρ j |{overscore (s)} l |: Replacement magnitude of filtered data.

μ j : Unitless muting coefficient that affects the magnitudes of data coefficient replacement values. $$ μ j = 1 ⁢   ⁢ for ⁢   ⁢ &LeftBracketingBar; d ^ j &RightBracketingBar; < σ j ⁢ &LeftBracketingBar; s _ j &RightBracketingBar; ⁢   ⁢ and ⁢   ⁢ μ j = ρ j ⁢ &LeftBracketingBar; s _ j d ^ j &RightBracketingBar; ⁢   ⁢ for ⁢   ⁢ &LeftBracketingBar; d ^ j &RightBracketingBar; ≥ σ j ⁢ &LeftBracketingBar; s _ j &RightBracketingBar; .

{circumflex over (d)} j : Transform data in subband j. $$ α j = ρ j ⁢ &LeftBracketingBar; s _ j d ^ j &RightBracketingBar;

We may express μ j as: μ j =H(σ j |{overscore (s)} j ||{circumflex over (d)} j |)+α j H ({circumflex over (d)} j |−σ j |{overscore (s)} j |) where H is the Heaviside function defined as H (t−t′)=0 for t<t′ and H(t−t′)=1 for t>t′.

Rewriting μ j :

μ j =1 −H ({circumflex over (d)} j |−σ j |{overscore (s)} j |)+α j H ({circumflex over (d)} j |−σ j |{overscore (s)} j |)=1−(1−α j )H ({circumflex over (d)} j |−σ j |{overscore (s)} j |).

Note that μ j =μ j ({circumflex over (t)}) that is, μ j depends on the transform time coordinate {circumflex over (t)} k .

μ j ({circumflex over (t)} k )=1=(1−α j ({circumflex over (t)} k ) H ({circumflex over (d)} j ({circumflex over (t)} k )|−σ j |{overscore (s)} j |) shows that muting coefficient μ j depends on {circumflex over (t)} k thru {circumflex over (d)} j (and α j ). Note: k indexes time coordinate {circumflex over (t)}within subband j.

Travel time trajectories are defined to parameterize regions of time-space coordinates over which muting coefficient μ j ({circumflex over (t)} k ) will be applied.

t b : Nominal beginning time of filter.

t e : Nominal ending time of filter.

t b.j (X; τ b.i , p b.j ): Scale dependent travel time trajectory for beginning times as function of space-coordinate x.

x: Space coordinate.

τ b.j : Beginning time at x=0.

p b.j : Parameterizes trajectory shape.

Line trajectory: t b.j =τ b.i +p b.j X.

Parabola trajectory: t b.j =τ b.j +p b.j X 2 .

Hyperbola trajectory: t b.j 2 =τ b.j 2 +p b.j 2 X 2 .

t e.j (x; τ e.j , p e.j ): Scale dependent travel time trajectory for ending times.

With nominal beginning and ending times given, respectively, as t b.j (x; τ b.j , p b.j ) and t e.j (x; τ e.j , p e.j ), we may write that {tilde over (m)} j ({circumflex over (t)} k )=1−(1−μ j ( {circumflex over (t)} k )H({circumflex over (t)} k −t b )H(t e −{circumflex over (t)} k ) =1−(1−μ j ({circumflex over (t)} k ))(1−H({circumflex over (t)} k −t e ))H({circumflex over (t)} k −t b ) where {tilde over (m)} j ({circumflex over (t)} k ) is an element of the operator {tilde over (M)} j that applies transform-time dependent muting to the transform data {circumflex over (d)} j ({circumflex over (t)} k ).

FIG. 5 is an example of a muting operator for a Discrete Wavelet Transform (DWT). Take N=2 j =16, the number of samples in {circumflex over (d)}, the discrete wavelet transform data vector 503 where J=4. For j=2, the subband to be filtered, FIG. 5 is {tilde over (M)} 2 , an example of {tilde over (M)} j . To generalize the {tilde over (M)} j DWT example for N=2 j data samples, vector {tilde over (m)} j ({circumflex over (t)} k ) may be mapped into 501 matrix {tilde over (m)} il (j) , the (i,l)th element of {tilde over (M)} j (i,l=1,2, . . . ,2 j ).

Define {tilde over (m)} j ({circumflex over (t)} k ) Δ{tilde over (m)} pq , with p=q=2 j +k, where j=1,2, . . . , J−1 lists eligible subbands (j=0 is not a wavelet subband but is necessary for reconstruction), and k=1,2, . . . ,2 j lists eligible time indices within subband j. Therefore, p and q range over all integers from 2 j +1 so (k=1) to 2 1+l where(k=2 j ).

These definitions give {tilde over (m)} il (j) =δ il +δ ip ( {tilde over (m)} pq −δ pq )δ pl for i,l=1,2, . . . ,2 j . p=q=2 j +k; p,q=2 j +1,2 j +2, . . . ,2 j+l . j=1,2, . . . , J−1; and k=1,2, . . . ,2 j .

A flow chart of the present invention is illustrated in FIG. 7 . Data are obtained 701 , for example a plurality of spatially related geophysical time series distributed over the earth, said geophysical time series containing data, the data being characterized by a travel time, a position, and an amplitude. These data may contain non-stationary energy difficult or impossible to remove using prior art methods. Portions of geophysical time series data may be more likely to contain undesirable energy components. At least a portion of these data transformed 703 using a wavelet transform that can produce transform coefficients. The transform coefficients are then organized 705 into a plurality of sub-band traces. The transform coefficients can then be modified 707 within at least one of the sub-band traces. For example, the modification may take place for a time window, a frequency window or a spatial window, to name a few choices. As disclosed by the explanations of the methods above, the chosen filter action may be targeted or limited within specified subband(s) to a computed time-space neighborhood of the undesirable signal present in geophysical time series data ensembles. The data are then inverted 709 using an inverse transform to produce filtered versions of the originally transformed data. As previously disclosed above, operator T −1 M j T describes how the entire filtering procedure in matrix notation may be applied.

Persons skilled in the art will understand that the methodology and processes for filter design and application described herein may be practiced in various forms including but not limited to the particular embodiments described herein. Further, it should be understood that the invention is not to be unduly limited to the foregoing which has been set forth for illustrative purposes. Various modifications and alternatives will be apparent to those skilled in the art without departing from the true scope of the invention, as defined in the following claims.

Claims

1. A method of filtering geophysical time series, the method comprising:(a) obtaining a plurality of spatially related geophysical time series distributed over the earth, said geophysical time series containing data, said data being characterized by a travel time, a position, and an amplitude; (b) transforming at least a portion of said spatially related geophysical time series using a wavelet transformation, said wavelet transformation producing transform coefficients; (c) organizing said transform coefficients into a plurality of sub-band traces; (d) modifying at least one of said transform coefficients within at least one of said sub-band traces; and, (e) inverting said sub-band traces using a wavelet transform inverse, thereby producing a filtered version of said transformed portion of said spatially related geophysical time series.

2. A method according to claim 1 wherein modifying at least one of said transform coefficients includes the step of eliminating at least one sub-band of said wavelet transformation.

3. A method according to claim 1 wherein modifying at least one of said transform coefficients includes eliminating all but one sub-band of said wavelet transformation.

4. A method according to claim 1 wherein modifying at least one of said transform coefficients includes increasing a value of at least one transform coefficient.

5. A method according to claim 1 wherein modifying at least one of said transform coefficients includes decreasing a value of at least one transform coefficient.

6. A method according to claim 1 wherein said wavelet transform is a discrete wavelet transform and said transform coefficients are wavelet transform coefficients.

7. A method according to claim 1 wherein transforming at least a portion of said spatially related geophysical time series further comprises the steps of:(a) selecting at least one geophysical time series from said plurality of spatially related geophysical time series; and, (b) applying a wavelet transformation to each geophysical time series so selected, said wavelet transformation i) being applied to a windowed portion of said geophysical time series, said windowed portion having beginning and ending times, ii) being characterized by at least one wavelet subband, and iii) producing a plurality of transform coefficients associated with said at least one wavelet subband when so applied.

8. A method according to claim 7 wherein for each selected spatially related geophysical time series:(a) said plurality of transform coefficients are attenuated if the magnitudes of said coefficients exceed a pre-determined threshold value, and (b) said coefficients that exceed a predetermined threshold value are attenuated to a pre-determined data value level.

9. A method according to claim 8 wherein said predetermined threshold value is determined from at least one of the group consisting of i) mean of the transform coefficient values, ii) median of transform coefficient values, iii) variance of transform coefficient values, and iv) a ratio of transform coefficient value statistics.

10. A method according to claim 8 wherein determination of attenuation value level of said coefficients is chosen with an arbitrary multiplicative parameter multiplied by a data dependent parameter, said data dependent parameter chosen from at least one of:(1) mean of transform coefficient values in subband, (2) median of transform coefficient values in subband, (3) variance of transform coefficient values in subband, and (4) ratios of transform coefficient values.

11. A method of determining filter parameters for geophysical time series, the method comprising:(a) obtaining a plurality of spatially related geophysical time series distributed over the earth, said geophysical time series containing data, said data being characterized by a travel time, a position, and an amplitude; (b) transforming at least a portion of said spatially related geophysical time series using a wavelet transformation, said wavelet transformation producing transform coefficients; (c) organizing said transform coefficients into a plurality of subband traces; (d) modifying at least one of said transform coefficients within at least one of said subband traces; (e) inverting said sub-band traces using a wavelet transform inverse, thereby producing a filtered version of said spatially related geophysical time series, and (f) determining which said subband traces contain non-stationary signals.

12. A method according to claim 11 wherein modifying at least one of said transform coefficients includes eliminating at least one sub-band of said transform.

13. A method according to claim 11 wherein modifying at least one of said transform coefficients includes eliminating all but one sub-band of said transform.

14. A method according to claim 11 wherein modifying at least one of said transform coefficients includes increasing at least one transform coefficient value.

15. A method according to claim 11 wherein modifying at least one of said transform coefficients includes decreasing at least one transform coefficient value.

16. A method according to claim 11 wherein said wavelet transform is a discrete wavelet transform and said transform coefficients are wavelet transform coefficients.

17. A method according to claim 11 wherein inverting said sub-band traces further comprises:(a) selecting at least one geophysical time series from said plurality of spatially related geophysical time series; and, (b) determining a subband of a wavelet transformation applied to each geophysical time series so selected, said subband containing non-stationary seismic signal.

18. A method according to claim 17 wherein for each selected spatially related geophysical time series:(a) said non-stationary signal is present within a windowed region of said selected geophysical time series, and (b) determining beginning and ending time parameters for said windowed region containing non-stationary signal of said selected geophysical time series.

19. A method according to claim 11 further comprising determining at least one statistic of transform data coefficients.

20. A method according to claim 19 wherein the at least one statistic of transform data coefficients is selected from the group consisting of i) mean of the transform coefficient values, ii) median of transform coefficient values, iii) variance of transform coefficient values, and iv) a ratio of transform coefficient value statistics.

21. A method according to claim 18 further comprising determining at least one statistic of transform data coefficients within said window.

22. A method according to claim 21 wherein the at least one statistic of transform data coefficients within said window is selected from the group consisting of i) mean of the transform coefficient values, ii) median of transform coefficient values, iii) variance of transform coefficient values, and iv) a ratio of transform coefficient value statistics.

23. A method according to claim 11 further comprising determining filter parameters to suppress non-stationary signals in said subband traces.

24. A method according to claim 11 further comprising determining filter parameters to modify non-stationary signals in said subband traces.

25. A method of determining filter parameters for geophysical time series, the method comprising:(a) obtaining a plurality of spatially related geophysical time series distributed over the earth, said geophysical time series containing data, said data being characterized by a travel time, a position, and an amplitude; (b) transforming at least a portion of said spatially related geophysical time series using a wavelet transformation, said wavelet transformation producing transform coefficients; (c) organizing said transform coefficients into a plurality of subband traces; (d) inverting said sub-band traces using a wavelet transform inverse, and (e) determining which said subband traces contain non-stationary signals.

26. A method according to claim 25 wherein said wavelet transform is a discrete wavelet transform and said transform coefficients are wavelet transform coefficients.

27. A method according to claim 25 wherein inverting said sub-band traces further comprises:(a) selecting at least one geophysical time series from said plurality of spatially related geophysical time series; and, (b) determining a subband of a wavelet transformation applied to each geophysical time series so selected, said subband containing non-stationary seismic signal.

28. A method according to claim 27 wherein for each selected spatially related geophysical time series:(a) said non-stationary signal is present within a windowed region of said selected geophysical time series, and (b) determining beginning and ending time parameters for said windowed region containing non-stationary signal of said selected geophysical time series.

29. A method according to claim 25 further comprising determining at least one statistic of transform data coefficients.

30. A method according to claim 29 wherein the at least one statistic of transform data coefficients is selected from the group consisting of i) mean of the transform coefficient values, ii) median of transform coefficient values, iii) variance of transform coefficient values, and iv) a ratio of transform coefficient value statistics.

31. A method according to claim 28 further comprising determining at least one statistic of transform data coefficients within said window.

32. A method according to claim 31 wherein the at least one statistic of transform data coefficients within said window is selected from the group consisting of i) mean of the transform coefficient values, ii) median of transform coefficient values, iii) variance of transform coefficient values, and iv) a ratio of transform coefficient value statistics.

33. A method according to claim 25 further comprising determining filter parameters to suppress non-stationary signals in said subband traces.

34. A method according to claim 25 further comprising determining filter parameters to modify non-stationary signals in said subband traces.

35. A method of filtering geophysical time series, the method comprising:(a) obtaining a plurality of spatially related geophysical time series distributed over the earth, said geophysical time series containing data said data being characterized by a travel time, a position, and an amplitude; (b) transforming at least a portion of said spatially related geophysical time series using a wavelet transformation, said wavelet transformation producing sub-band traces with transform coefficients; (c) modifying at least one of said transform coefficients within at least one of said sub-band traces; and, (d) inverting said sub-band traces using a wavelet transform inverse, thereby producing a filtered version of said spatially related geophysical time series.

36. A method according to claim 35 wherein modifying at least one of said transform coefficients includes the step of eliminating at least one sub-band of said wavelet transformation.

37. A method according to claim 35 wherein modifying at least one of said transform coefficients includes eliminating all but one sub-band of said wavelet transformation.

38. A method according to claim 35 wherein modifying at least one of said transform coefficients includes increasing a value of at least one transform coefficient.

39. A method according to claim 35 wherein modifying at least one of said transform coefficients includes decreasing a value of at least one transform coefficient.

40. A method according to claim 35 wherein said wavelet transform is a discrete wavelet transform and said transform coefficients are wavelet transform coefficients.

41. A method according to claim 35 wherein transforming at least a portion of said spatially related geophysical time series comprises the steps of:(a) selecting at least one geophysical time series from said plurality of spatially related geophysical time series; and, (b) applying a wavelet transformation to each geophysical time series so selected, said wavelet transformation i) being applied to a windowed portion of said geophysical time series, said windowed portion having beginning and ending times, ii) being characterized by at least one wavelet subband, and iii) producing a plurality of transform coefficients associated with said at least one wavelet subband when so applied.

42. A method according to claim 41 wherein for each selected spatially related geophysical time series:(a) said plurality of transform coefficients are attenuated if the magnitudes of said coefficients exceed a pre-determined threshold value, and (b) said coefficients that exceed a predetermined threshold value are attenuated to a pre-determined data value level.

43. A method according to claim 42 wherein said predetermined threshold value is determined from at least one of the group consisting of i) mean of the transform coefficient values, ii) median of transform coefficient values, iii) variance of transform coefficient values, and iv) a ratio of transform coefficient value statistics.

44. A method according to claim 42 wherein determination of attenuation value level of said coefficients is chosen with an arbitrary multiplicative parameter multiplied by a data dependent parameter, said data dependent parameter chosen from at least one of:(1) mean of transform coefficient values in subband, (2) median of transform coefficient values in subband, (3) variance of transform coefficient values in subband, and (4) ratios of transform coefficient values.

45. A method according to claim 42 further comprising obtaining the predetermined level based upon a data dependent parameter, said data dependent parameter chosen from at least one of:(1) mean of transform coefficient values in subband, (2) median of transform coefficient values in subband, (3) variance of transform coefficient values in subband, and (4) ratios of transform coefficient values.