Minimum phase bandpass filtering of seismic signals

Abstract

A process for better defining first break seismic energy travel times for correlated seismic return signals that consists of constructing a minimum phase bandpass filter having bandwidth narrower than the vibrational seismic return signals and convolving the filter with the return signals to convert the seismic return signals to a causal wave form.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention

The invention relates generally to a method of filtering correlated seismic signals in order to define the first break with greater accuracy thereby to improve derivation of static correction data and, more particularly, but not by way of limitation, it relates to an improved form of digital filtering for defining travel times of near surface data.

2. Description of the Prior Art

Applicant knows of no prior art that is directed to a method for accurately determining the near surface or refraction signal arrival times on correlated vibrational type seismic data. Current methods involve choosing of an arbitrary position or amplitude on the first positioned vibrational wavelet, i.e. correlated wavelet, or in using only relative time changes from trace to trace rather than attempting to pick the true travel times.

SUMMARY OF THE INVENTION

The present invention relates to an improved method for computer implemented digital filtering for use with correlated vibrational seismic signals wherein filtering with a selected minimum phase bandpass filter effectively transforms the two-sided symmetric first break wavelet to a one-sided causal waveform having the attribute that the exact onset time of first arrival is located at the first point of departure of the return signal trace from the base line. Thus, the process functions to construct a bandpass filter in the frequency domain with zero phase; then to inverse Fourier transform to the time domain and to truncate the zero phase, symmetric filter pulse operator to a specified length; thereafter, a Hilbert transform is performed to obtain minimum phase spectrum for final formation of a filter operator in the time domain having minimum phase for convolution with the seismic traces under consideration.

Therefore, it is an object of the present invention to provide a filtering method for correlated seismic data that more accurately defines the first arrival time of seismic return signals, particularly refracted signals.

It is also an object of the present invention to provide a relatively simple digital filtering method that is readily operative with correlated seismic signals to accurately process refracted return data.

It is still further an object of the present invention to provide a method for establishing accurate first break times of correlated seismic signals despite the presence of correlation side lobes and other arrivals.

Finally, it is an object of the invention to design a digital filter which can change the phase of the vibrational seismic signal wavelet, i.e. as correlated, and alter its amplitude spectrum in such a way that it appears to become a causal wavelet, i.e. uncorrelated.

Other objects and advantages of the invention will be evident from the following detailed description when read in conjunction with the accompanying drawings which illustrate the invention.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1A is a graphic depiction in amplitude versus time of a first arrival portion of seismic signal return from a causal source;

FIG. 1B is a graphic depiction in amplitude versus time of an idealized symmetric seismic signal, as results from the crosscorrelation process;

FIG. 1C is a graphic depiction in amplitude versus time of a portion of near surface vibrational seismic signal after correlation;

FIG. 2 is a flow diagram illustrating first steps of the process of the present invention;

FIG. 3 is a flow diagram illustrating final steps of the process, including convolution of a seismic trace set;

FIG. 4A is a computer printout of a typical correlated seismic wavelet;

FIG. 4B is an amplitude spectrum for the wavelet illustrated in FIG. 4A;

FIG. 5A is a computer printout of a causal wavelet resulting from application of the method of the invention to the correlated wavelet of FIG. 4A;

FIG. 5B is an adjusted amplitude spectrum making up the minimum phase filter that produced the causal waveform of FIG. 5A.

DETAILED DESCRIPTION

The present invention functions to produce a minimum phase bandpass filter that is effective to alter the time/phase relationship of a two-sided or symmetric type of correlation wave to produce a one-sided type of wave that is causal in appearance. That is, a causal wave similar to the characteristic noncorrelated seismic return wave forms. Thus, in using specific program implemented computer, a filter is designed which can change the phase of the correlated signal wavelet and alter its amplitude spectrum in such a way that it appears to become a causal wavelet. Very simply, this enables more accurate determination of return signal first break as it may be picked at the point of first trace departure from the trace base line. Thus, refracted seismic data travel time can be more accurately determined for static corrections.

As is well known, the vibrational method of seismic exploration involves the transmitting of a sinusoidal sweep of selected frequency range and duration into the earth, and thereafter recording the return signals from specified transducer locations. The sweep signal data is then further defined by crosscorrelating recorded return seismic signals with a replica of the input sweep or control signal. This crosscorrelation process produces a zero phase, two-sided symmetric waveform which has its peak centered at the location of the event. That is, the center of the primary or largest amplitude pulse indicates actual travel time or event location while the entire wavelet proceeds on both sides in equally diminishing side lobes. Thus, if many reflection events are present at the same or near same time, the wavelets can become superimposed in interference with one another, tending to alter or mask over any peak wave form which might be examined for exact time or event location. Thus, as shown in FIG. 1A, a causal waveform as processed without use of crosscorrelation, e.g. return signal from a dynamite source, is illustrated at waveform 10 and it may be readily seen that first break signal travel time may be picked exactly at the point of departure from the base line 12, i.e. point 14 as illustrated by dash-line 16.

FIG. 1B illustrates a correlated wavelet 18 as might result from crosscorrelation of vibrational seismic data and in this case the event location or exact wave time is at peak center 20 as illustrated by dash-line 22. Thus, the typical correlation wavelet such as wavelet 18 has the primary center lobe plus a plurality of symmetrically disposed side lobes on each side. In the usual case of seismic signal return, there is not a clear cut signal pulse return but a grouping of plural return pulses which tend to partially overlay and interfere each with the other as far as concise wave form definition is concerned. FIG. 1C illustrates a typical type of first break interference on a seismic trace has been crosscorrelated with control signal to produce the final trace. The seismic trace 24 shows signal on baseline 12 at point 26; however, the actual arrival time for static correction purposes would lie at the center of the peak of the first correlation wavelet, and trace 24 indicates a slurred and indistinct peak 28 showing presence of several first break pulses in close time relationship. Thus, the actual time zero or first break time will fall some place between the dash-lines indicated by bracket 30.

In order to more accurately process the profile of seismic traces through later specific enhancement procedures, it is first desirable to derive an exact time relationship for the near surface, from first arrival signals such as might be refracted from the weathering layer and other shallow interfaces. In some cases, the reflectors or refractors lie as a layered medium, and in still other cases the shallow refracting medium may appear in a gradational velocity relationship as where sucessive compaction of earth particles has taken place such that shallow surface wave travel is generally arcuate in propagation and return. Whether the weathering velocity is gradational or exhibits fairly sharp layering, travel times for energy withing this low velocity sector become very critical when processing to develop the static correction times that may then be applied to the overall seismic section data to eliminate the effect of variations in elevation, weathering thickness or velocity.

The present invention processes correlated seismic signals, such as trace 24 of FIG. 1C, to alter the signal in the frequency domain for reconstitution in time domain as a causal-type wavelet, such as trace 10 of FIG. 1A, wherein the first break or first departure from baseline 12 will give a true first arrival time indication. The digital filter which converts the two-sided, symmetric vibrational wavelet to a causal-type one-sided wavelet is a minimum phase bandpass filter in which the bandpass is narrower than the bandwidth of the vibrational seismic wavelets. This digital filtering is effected by the following steps:

(1) Constructing a bandpass filter in the frequency domain having zero phase;

(2) inverse Fourier transforming to the time domain and truncating the zero phase, symmetric filter pulse to a specified operator length, typically 101 points;

(3) performing a Hilbert transform on the logarithm of the amplitude spectrum of the filter pulse to obtain the minimum phase spectrum;

(4) forming a new filter pulse with the same amplitude spectrum but minimum phase; and

(5) convolving the minimum phase bandpass wavelet with seismic traces.

This process is illustrated in the flow diagrams of FIGS. 2 and 3 as taken in series. The program list out for the process is enumerated hereinafter, and the numbers to the right of the flow blocks in FIGS. 2 and 3 illustrate the program number for the particular flow step procedure. A first flow stage 32 carries out inputting of corner frequencies and filter length after initial header information is cleared through the program. Program steps 46 through 83 carry out the input functions in accordance with selected filter attributes as determined from initial considerations of the seismic survey. Thus, initial frequency and final frequency limits will be input in accordance with the original vibrator control signal band width including the on and off tapers. It should be considered, too, that return seismic signal bandwidth may be somewhat narrower than the original control signal bandwidth and this should be taken into consideration when selecting the minimum phase filter parameters. Next flow stage 34 then functions to compute the bandpass filter in accordance with the prior set parameters.

Flow stage 36 performs inverse Fourier transform on the filter data to convert the data to the time domain, i.e. time is the independent variable upon which operation may be made. The time domain data is then processed in flow stage 38 to truncate the symmetric filter to a desired operator length. Thus, the time domain length of the impulse response in the nature of a convolution operator may be specified by a certain number of points, typically 101 points, which has specific relationship to the seismic date sampling rate and number of sampling intervals. Flow stage 40 then functions to Fourier transform the filter operator back to the frequency domain.

Flow stage 42 then functions to process the frequency domain data to compute the power spectrum as the magnitude squared at each frequency. The power spectrum is the square of the amplitude-frequency response, i.e. the Fourier cosine transform of the autocorrelation function. Thereafter, stage 44 computes one half times the natural logarithm of the power spectrum and flow stage 46 again inverse Fourier transforms the data to the time domain. Flow stage 48 then multiplies the time function by a step function selected to be of amplitude 2 and onset at t=0. This data is again Fourier transformed in stage 50 to the frequency domain, where frequency is the independent variable.

Referring to FIG. 3, the data is operated on in stage 52 to form a complex exponential of the frequency spectrum and this data is inverse Fourier transformed to the time domain in stage 54 and truncated in stage 56 to establish a digital filter pulse or filter operator of desired length for the particular seismic section undergoing processing. Thereafter, flow stage 58 inputs all traces of the section for processing with the minimum phase bandpass filter operator to convert the crosscorrelated seismic trace data to data having causal attributes wherein event time may be accurately picked. Given this output as to seismic energy travel time, the data can be further processed to establish static corrections for the trace data through various subsequent data enhancement processes.

The basic computer program performing the process of the present invention may be written in Fortran 5 language for operation on such as a Control Data Corporation computer commercially available as CYBER model 760. It should also be understood that the program process of the invention is fully compatible with hybrid computer applications utilizing relatively simple microprocessor circuitry. Such selection of hardware and requisite variances in program are well within the skill of the artisan.

A full listout for the minimum phase bandpass filter is as follows: ##SPC1## ##SPC2##

The result of applying a minimum phase bandpass filter to a zero phase wavelet may be illustrated by means of computer printouts drawn in accordance with the present invention.

FIG. 4A illustrates graphically a zero phase wavelet having an amplitude spectrum such as shown in FIG. 4B. The amplitude spectrum tapers on from three to seven cycles, maintains a unit value from seven to fifty cycles, and tapers off from fifty to fifty-four cycles. A printout of a zero phase wavelet dictated by such amplitude spectrum is indicated by wavelet 60 of FIG. 4A. Wavelet 60 is printed out on a baseline 62 representative of time and the event time for wavelet 60 is indicated as 0.5 seconds, the center most position of symmetrical wave peak 64.

As shown in FIGS. 5A and 5B, a similar type of computer printout illustrates that a casual-type wave 66 results from application of a minimum phase filter operator constructed from an amplitude spectrum such as that shown in FIG. 5B. Thus, the amplitude spectrum outlines a narrower bandwidth of taper on from eight cycles to twelve cycles, unit amplitude from twelve to forty-five cycles, and taper off from forty-five to forty-nine cycles. Computer printout of the filtered casual wave 66 indicates that the first breakpoint from the baseline 68 is now at 0.50 seconds which is the same as the event indication for the symmetrical wavelet of FIG. 4A. Thus, location of the event has been shifted to the onset of the wavelet in casual form.

The foregoing discloses a novel method utilizing a digital filter and selected minimum phase operator for converting vibrational seismic signal refraction arrivals into a casual form so that the first break or onset of the refraction arrival may be accurately determined. Such accurate time determinations are valuable for use in subsequent processing to compute the static corrections for the particular set of reflection data along whatever the selected survey site. The invention relies on the discovery that the minimum phase bandpass filter has the capability of effecting the necessary changes in the wave form thereby to produce a casual appearance, and it should be understood that the process may be variously effected by equivalent hardware and program implementation.

Changes may be made in combination and arrangement of elements as heretofore set forth in the specification and shown in the drawing; it being understood that changes may be made in the embodiments disclosed without departing from the spirit and scope of the invention as defined in the following claims.

Claims

1. A method of filtering crosscorrelated seismic data comprising:

constructing a selected bandpass filter operator in the frequency domain having zero phase;

transforming said bandpass filter operator and obtaining a minimum phase spectrum, and forming a frequency domain filter operator with the same amplitude spectrum but minimum phase;

again transforming said frequency domain filter operator to the time domain to form said minimum phase bandpass filter operator; and

convolving said minimum phase bandpass filter operator with said crosscorrelated seismic data to define the data first arrival time.

2. A method as set forth in claim 1 wherein said step of transforming comprises:

inverse Fourier transforming said bandpass filter operator to the time domain and truncating to obtain a zero phase, symmetric filter operator of selected length; and

transforming the logarithm of the amplitude spectrum of the symmetric filter operator to obtain the minimum phase frequency domain filter operator.

3. A method as set forth in claim 1 wherein said step of again transforming includes:

an inverse Fourier transform of said frequency domain filter operator.

4. A method as set forth in claim 2 wherein said step of again transforming includes:

an inverse Fourier transform of said frequency domain filter operator.

5. A method as set forth in claim 1 wherein said step of transforming comprises:

inverse Fourier transforming said bandpass filter operator to the time domain and truncating to produce a zero phase filter operator; and

computing one-half times the natural logarithm of the power spectrum of the zero phase filter operator and performing a Hilbert transform to produce said minimum phase frequency domain filter operator.

6. A method as set forth in claim 5 wherein said step of again transforming comprises:

forming the complex exponential of the power spectrum; and

inverse Fourier transforming to the time domain and truncating to selected operator length to produce said minimum phase bandpass filter operator for convolution with said crosscorrelated seismic data.

7. A method as set forth in claim 1 wherein said step of transforming comprises:

inverse Fourier transforming said bandpass filter operator to the time domain and truncating to obtain a zero phase, symmetric filter operator of selected length;

Fourier transforming said symmetric filter operator to the frequency domain and computing the power spectrum as the magnitude squared at each frequency;

computing one-half times the natural logarithm of the power spectrum;

inverse Fourier transforming said power spectrum logarithm to the time domain and multiplying by a step function of selected amplitude two and onset at time equals zero; and

Fourier transforming said multiplied time function.

8. A method as set forth in claim 2 wherein said step of again transforming comprises:

forming the complex exponential of the logarithm of the power spectrum of the minimum phase frequency domain filter operator; and

inverse Fourier transforming to the time domain and truncating to selected operator length to produce said minimum phase bandpass filter operator for convolution with said crosscorrelated seismic data.

9. A method as set forth in claim 7 wherein said step of again tansforming comprises:

forming the complex exponential of the logarithm of the power spectrum of the frequency domain filter operator; and

inverse Fourier transforming to the time domain and truncating to selected operator length to produce said minimum phase bandpass filter operator for convolution with said crosscorrelated seismic data.