Method for Seismic Trace Interpolation Using a Forward and Backward Application of Wave Equation Datuming

1. Field of the Invention

This invention relates to a method for interpolating between measured waveforms obtained from an acoustic, elastic or other experiment which obeys the wave equation. Interpolation between measured waveforms is accomplished by successive forward and backward applications of wave equation datuming. This invention further relates to seismic exploration and to a method for seismic trace interpolation by applying wave equation datuming techniques.

2. Description of the Prior Art

Physical experiments are often conducted in which recorded measurements which pertain to the outcome of the physical experiment obey the wave equation. Most often, the measurements recorded in such an experiment are related to the time history of the particle motion, displacement, acceleration or pressure at the location of a transducer. For example, in the field of underwater acoustics, a source may be used to generate sound waves into the ocean. An array of hydrophones towed behind a boat will be used to record the sound waves impinging upon the array. The recorded sound waves, which may be used to provide information regarding the location of objects beneath the ocean surface, will obey the wave equation. In the field of exploration seismology, an example of a physical experiment which obeys the wave equation would be the generation of seismic energy by a source and the use of an array of geophones spread out on the surface of the earth for recording elastic waves which are reflected from the subsurface and which impinge upon the array.

In seismic exploration, it is common practice to deploy a large array of geophones on the surface of the earth and to record the vibrations of the earth at each geophone location to obtain a collection of seismic traces. The traces are sampled and recorded for further processing. When the vibrations so recorded are caused by a seismic source activated at a known time and location, the recorded data can be processed by a computer in known ways to produce an image of the subsurface. The image thus produced is commonly interpreted by geophysicists to detect the possible presence of valuable hydrocarbons.

Seismograms are commonly recorded as digital samples which represent the amplitude of a received seismic signal as a function of time. Since seismograms are usually obtained along a line of exploration on the surface of the earth, the digital samples can be formed into x-t arrays with each sample in the array representing the amplitude of the seismic signal as a function of horizontal distance and time. When such arrays are visually reproduced, by plotting or the like, seismic sections are produced. A seismic section depicts the subsurface layering of a section of the earth. It is the principal tool which the geophysicist studies to determine the nature of the earth's subsurface formations. Before an array of seismic samples or traces can be converted into a seismic section for interpretation by geophysicists, the array must be extensively processed to remove noise and to make reflection events discernable.

It frequently occurs, however, that insufficient data is recorded data for proper interpretation of the characteristics of the subsurface formation. The most common data insufficiency problem is that the acquired data includes either seismic traces with no recorded data or seismic traces which clearly contain severe noise contamination such as spikes. These and other insufficiencies in the acquired data may be traced to any one of numerous causes. For example, a data insufficiency may result when seismic traces are recorded in the wrong location due to intrabed reflections of the generated seismic energy. Another typical data insufficiency arises when it is determined after seismic exploration has been completed that the collected traces need to be interpolated to a finer spatial sampling for proper interpretation. Here, the selection of too large a spacing between geophones prior to commencing exploration is typically the cause. A seismic trace without data may be produced when one or more geophones fail during exploration. Seismic traces severely contaminated with noise may be inadvertently produced due to multiple reflections or ground roll.

Too often, cost or time restraints will prevent the reshooting of the unsatisfactory data. Under such circumstances, the most acceptable option available to the geophysicist who must complete an evaluation of the formation based upon the insufficient data, is to estimate the missing or bad traces from the set of collected traces.

There are numerous methods for determining estimates of data points from measured physical data which may be classified as interpolation methods. For example, many methods exist for the interpolation of functions which are dependent on only one parameter such as time. These include both the well known linear interpolation methods and the equally well known finite difference interpolation methods. Linear interpolation is based upon the proposition that for an argument "x", the value of a dependent variable "y" will vary linearly over the interval under investigation. Finite difference interpolation methods are based upon the proposition that given "y=f(x)", inverse interpolation may be used in a process for finding "x" at a value of "y" intermediate to two tabulated values. For finite difference interpolation methods, any one of several alternate procedures such as Lagrange's formula, successive approximation, or reversion of series may be used to obtain the desired result.

Standard practice among geophysicists faced with seismic traces which contain either no recorded data or which are severely contaminated has been to exclude such traces, commonly referred to as "null" traces, from the otherwise satisfactory data set. The collected seismic data would be processed using conventional methods without the excluded data. When the missing trace was necessary for proper processing of the seismic data, prior attempts to restore the missing trace and create a trace with events consistent with nearby coherent events focused upon combining traces near the missing trace in the x-t domain to create the missing trace. The use of such prior methods restricted conventional multi-dimensional trace interpolation upon the apparent characteristics of the measured traces and did not permit trace interpolation to use the characteristics of the underlying physical experiment being performed.

The above mentioned trace interpolation methods are inherently 1-dimensional, for example, f(x)=b. These methods may not be applied to data sets which are 2 or 3-dimensional, such as f(x,t)=A, or f(x,y,t)=B. For example, a seismic trace measures amplitude as a function of time. Thus, each seismic trace in the gather is used to analyze subsurface formations which are separated by a distance dependent on the location of the geophones along the surface. This spatial variation may be only in one, two or three directions, depending on the location of the geophones utilized in the gather.

Alternate methods to the interpolation methods described above were developed which address the problems of interpolating functions having spatial and temporal dependencies. Too often, however, such alternate methods are based upon the observed statistics of the trace gather rather than upon the physics from which the experiment came. One such conventional solution which relies upon the statistics of the trace gather and is popular in the field of seismic exploration utilizes t frequency-wavenumber transform approach.

In the processing of seismograms, x-t arrays have been transformed into arrays representing amplitude as a function of frequency and wave number. This is commonly referred to as a "frequency-wave number" or "f-k" transformation. The f-k transformation has been used as a tool to study seismic filtering. F-k transforms are routinely used to represent data collected by a large array of sensors. Usually, the f-k representations are computed by Fast Fourier Transforms (hereafter referred to as FFTs) which transform the data from the time and spatial domain to the frequency and wavenumber domain by utilizing the corresponding multi-dimensional Fourier transform. The resulting data representations are parameterized by frequencies, wave numbers (spatial frequencies), amplitudes, and phases. More specifically, for each frequency there is a collection of wave numbers, and for each frequency-wave number pair there is a complex number representative of an amplitude and a phase. The data is then inverse Fast Fourier transformed back to the time and spatial domain. After this inverse transform, the traces have a different spatial sampling and possible a different temporal sampling.

In U.S. Pat. No. 4,594,693 issued to Pann et al, seismic trace interpolation is carried out by inserting zero amplitude traces between the seismic traces in a section where spatial aliasing is a problem. The traces are then transformed into an f-k array. The f-k array is filtered with a filter which rejects samples in a region of frequency and wave number which exhibits aliasing. The filtered f-k array is then transformed back into a seismic section representing amplitude as a function of time and distance. One limitation on the use of such methods, however, it that seismic interpolation using f-k transforms requires that the measured traces be uniformly sampled in space. This requirement of a regular trace spacing limits the usefulness of f-k transforms for seismic interpolation applications.

Finally, it is known to use a technique generally referred to as "velocity replacement" to remove the distortion of information gathered during seismic exploration which is caused by changes in the velocity of seismic energy between irregular, i.e. non-planar, layers in the subsurface formation being explored. Due to the varied characteristics of the layers, seismic energy travels at different velocities in each layer. To remove distortion from seismic data recorded at the surface which has traveled through such an irregular subsurface layer, the seismic data is propagated at a first velocity to a first real datum, corresponding to the bottom of the irregular layer, and then propagated at a second velocity to a second real datum, corresponding to the top of the layer. Further, as each layer in the surface has different characteristics which affect the velocity of seismic energy passing through that layer, the actual velocity of the irregular layer is used as the first velocity, and the velocity of the medium below the layer is used as the second velocity. For this reason, velocity replacement methods have consistently used different velocities for each propagation.

It is an object of this invention to provide a method of performing trace interpolation using the characteristics of the underlying physical experiment conducted.

It is another object of this invention to provide a method of successively using wave equation datuming in the forward and backward directions for the interpolation of data.

Yet another object of this invention is to provide a method of propagating a data set located at a first datum at a first velocity to a fictitious datum and to repropagate the data at the same velocity to the original datum to interpolate data missing from the original data set.

In accordance with the present invention, there is provided herein a method for interpolating traces between existing traces obtained by an experiment which is governed by the acoustic or elastic wave equation. Such a method is particularly applicable to the interpolation of seismic traces between existing traces obtained during the exploration of a subsurface formation. Traces are recorded or measured upon some datum surface in accordance with standard data acquisition methods for the particular experiment. Forward and backward wave equation datuming of the data set, which typically includes both satisfactory and unsatisfactory data points, utilizes the existing traces to interpolate, i.e. estimate new measurements for the unsatisfactory data originally acquired.

More specifically, the method of the present invention is accomplished by the propagation of existing measured traces to a fictitious datum plane using wave theory. The original traces are thus datumed to the fictitious plane to create a new set of traces. The new traces are propagated, again by the application of wave theory, back to the original datum plane to create a third set of traces. This process is hereby defined as forward and backward wave equation datuming. The third set of traces produced by forward and backward wave equation datuming will include new traces to be used in place of the originally acquired unsatisfactory data.

In an alternate embodiment of the invention, the fictitious datum plane may be located above or below the measurement datum plane. In yet another embodiment of the invention, the contribution to the interpolated trace at the desired location from each original trace may be computed directly without having to specifically compute the traces on the fictitious datum. In still yet another embodiment of the invention, the method of the present invention is applied iteratively. On subsequent iterations, however, the estimate of the desired traces previously determined is included with the original traces when they are datumed a second time to and from the fictitious datum.

The above and other object, advantages and features of the invention will be more readily understood from the following detailed description of the invention, which is provided in connection with the accompanying drawings.

FIG. 1 illustrates, in flow chart form, the method of interpolating traces using a forward and backward application of wave equation datuming of the present invention;

FIG. 2 illustrates a typical seismic gather which includes at least one null trace;

FIG. 3 illustrates the seismic gather of FIG. 2 where the null trace or traces have been replaced by interpolated traces produced according to the method of FIG. 1; and

FIG. 4 illustrates, in flow chart form, an alternative embodiment of the method of the present invention.

For purposes of this application, the terms "trace" and "seismic trace" refer to the time history of the measurements at a transducer. The terms "trace gather" and "seismic section" refer to the collection of traces from all of the transducers. The term "datum" generally refers to the series of numerical quantities along the "z" axis which serves as a base reference for the seismic traces which take up a seismic section. The term "original datum" refers to the collection of seismic traces acquired during seismic exploration which are placed according to the curvature of the earth's surface. The term "wave equation datuming" refers to as an upward or downward continuation of U(x, z=z.sub.1, t) to produce U(x, z=z.sub.2, t). The term "forward and backward wave equation datuming" refers to the propagation of data from a first datum to a second datum followed by the propagation back to the first datum by successive applications of wave equation datuming techniques using the same velocity.

The purpose of wave equation datuming is to redefine the reference surface on which the sources and receivers appear to be located. Wave equation datuming in two dimensions is accomplished according to the equation: ##EQU1## where: U.sub.out(t) is one of a plurality of traces included as part of the output datum;

i is the location of the acquired data traces on the input datum;

.DELTA.x.sub.i is the separation between adjacent data trace locations for the input datum;

t.sub.i is the traveltime between an input trace on the input datum and an output trace on the output datum;

r.sub.i is the distance between the input trace location and the output trace location;

.theta..sub.i is the angle between r.sub.i and the normal to the input datum; and

Q(t-t.sub.i) is a representation of the input trace at location i delayed by the traveltime t.sub.i and convolved with a time-domain shaping operator. This shaping operator, for the 2-dimensional case, effects a half derivative of the original input trace, which is equivalent to multiplying each frequency component of the trace by .sqroot.iw , where i=.sqroot.-1, and w is frequency.

Wave equation datuming in three dimensions is accomplished according to the equation: ##EQU2## where: U.sub.out(t) is one of a plurality of traces included as part of the output datum;

i is the x location of the acquired data traces on the input datum surface;

k is the y location of the acquired data traces on the input datum surface;

.DELTA.x.sub.i is the x separation between adjacent data trace locations for the input datum surface;

.DELTA.y.sub.k is the y separation between adjacent data trace locations for the input datum surface;

t.sub.i is the traveltime between an input trace on the input datum and an output trace on the output datum;

r.sub.ik is the distance between the input trace location and the output trace location;

.phi.i.sub.k is the angle between r.sub.ik and the normal to the input datum; and

P(t-t.sub.ik) is a representation of the input trace at x location i and y location k delayed by traveltime t.sub.ik and convolved with a time-domain shaping operator. This shaping operator, for the three dimensional case, effects a derivative of the original input trace.

While the method of trace interpolation which follows is described with respect to the interpolation of seismic traces, it should be clearly understood that the application to seismic traces is exemplary only and that the methods of the present invention are equally applicable for the interpolation of missing traces for other collections of data which are related to the acoustic or elastic wave equation.

Turning first to FIG. 1, the method of seismic trace interpolation using a forward and backward application of wave equation datuming is hereby described. The method of the present invention commences at step 100 with the acquisition of seismic data using any one of numerous well known seismic exploration techniques. For example, an artificial disturbance may be generated along the earth by the use of dynamite or the like. The resulting acoustic waves travel downwardly in the earth and are reflected upward from subsurface reflecting interfaces. The reflected waves are received at geophones or other detectors located along the surface and recorded in reproducible form as seismograms. Seismic traces which would be acquired by a typical seismic survey may be seen by reference to FIG. 2. FIG. 2 shows the amplitude of seismic reflections as a function of time and distance along a line of exploration in the x direction of the earth's surface. These traces have been gathered into an x-t array commonly referred to as a "seismic section" or "seismic record". During a typical seismic exploration, a seismic trace would be produced for each geophone, the collection of seismic traces making up a seismic gather. For any number of the reasons previously set forth, the seismic gather may include a number of insufficient or otherwise unacceptable seismic traces. For example, FIG. 2 shows at details "a", "b" and "c", unacceptable seismic traces. For traces "a", "b" and "c", a blank or null trace has been substituted in the seismic gather for the original unacceptable seismic trace.

Proceeding to step 110, the insufficient or bad seismic traces are removed from the seismic gather and replaced with null traces. At step 120, the seismic gather it related to the original datum. Preferably, the original datum selected will be the earth surface. The earth surface datum, which should show the variations in the altitude of the earth surface between the geophones along the line of exploration, may be determined using conventional means. Each trace is related to the datum according to the location of the geophone which recorded the corresponding seismic trace.

Proceeding to step 130, a second datum, also referred to as the "fictitious" datum, is selected for later propagation of the seismic traces of the original datum. The selected next datum may be located above or below the original datum and may correspond to any arbitrary line equation. At step 140, the trace locations for the next datum corresponding to the trace locations on the earth datum are selected. Successive traces on the next datum shall be separated by a distance .DELTA.x.sub.i corresponding to the distance between geophone locations on the earth surface. Proceeding to step 150, .DELTA.x.sub.i, r.sub.i t.sub.i and .theta..sub.i are determined for each trace location on the next datum.

At step 160, propagation of the seismic traces from the original earth datum to the next, or fictitious, datum is commenced using the appropriate calculations as set forth in Equation (1) where the original or earth datum is the input datum and the next or fictitious datum is the output datum. Each seismic trace is forward datumed to the next datum where the fictitious datum may take on any arbitrary shape. In propagating the seismic traces located at the original earth datum to the fictitious next datum, each original trace and each data point included as part of an original trace will provide a component of each next trace produced by propagation of the original data to the next datum. The null traces are not propagated to the fictitious datum and will therefore not contribute to the series of next traces produced by propagation. Propagation of each trace to the fictitious datum continues until all original data traces have been propagated to the next datum. In such a manner, the set of new traces produced for the fictitious datum will include new traces which correspond with the null traces of the original seismic gather.

Proceeding to step 170, the set of seismic traces which have been produced by the propagation of the acquired data from the original datum to the fictitious datum are next propagated back to the original datum. Propagation of the produced seismic traces back to the original datum, a step which is also generally referred to as "backward datuming" of the seismic traces, is accomplished using the same calculations set forth in Equation (1) which were required for forward datuming of the original traces to the fictitious datum with the exception that the fictitious datum is now the input datum and the original datum is now the output datum. As before, each input trace and each data point included as part of the input trace will provide a component of each output trace produced by the propagation of produced data to the original datum. Propagation, or backward datuming, of each produced trace back to the original datum continues until all produced data traces have been backward datumed. In such a manner, a second set of new traces which corresponds to the original datum will include traces in place of the original null traces. The traces in the second set of new traces which correspond to the original null traces produced by backward datuming may be seen by reference to details "a", "b" and "c" of FIG. 3. As previously stated, the second set of new traces include traces which had been null prior to the forward and backward datuming steps.

At step 180, the second set of new seismic traces which correspond to the null or otherwise bad traces of the original seismic gather are substituted for the null traces of the seismic gather of FIG. 2. Proceeding to step 190, the seismic gather which include the newly interpolated seismic traces in place of the null traces are processed by any one of the numerous data processing techniques to produce useful information regarding the characteristics of the subsurface formation explored.

Turning next to FIG. 4, another embodiment of the present invention of a method for seismic trace interpolation using a forward and backward application of wave equation datuming is hereby described. Commencing at step 200, seismic data is acquired using any one of the previously mentioned well known seismic exploration techniques. Again, it is expected that the seismic gather acquired during exploration will include a number of seismic traces which are either missing or otherwise deemed unacceptable. Proceeding to step 210, the unacceptable traces are removed from the seismic gather and replaced with null traces. Null traces are also inserted into the data set where the missing traces should have been located. At step 220, the seismic gather is placed on an original datum, i.e. the earth's surface.

Proceeding to step 230, a second datum is selected for later propagation of the seismic traces of the original datum. At step 240, the trace locations for the next, fictitious datum corresponding to the trace locations on the original datum are selected. Proceeding to step 250, .DELTA.x.sub.i, r.sub.i t.sub.i and .theta..sub.i are determined for each trace location on the fictitious datum. At step 260, propagation of the seismic traces from the original earth datum to the fictitious datum is commenced using the wave equation datuming techniques where the original datum is the input datum and the fictitious datum is the output datum. The null traces are not, however, propagated to the fictitious datum. Propagation of each trace to the fictitious datum continues until all original data traces have been propagated to the fictitious datum.

Proceeding to step 270, the set of seismic traces which have been produced by the propagation of the acquired data from the original datum to the fictitious datum are next propagated back to the original datum by the aforementioned wave equation datuming techniques where the fictitious datum is the input datum and the original datum is the output datum. Propagation, or backward datuming of each produced trace back to the original datum continues until all produced data traces have been backward datumed. In such a manner, a second set of new traces which corresponds to the original datum and which will include traces in place of the original null traces is produced.

At step 280, the new seismic traces which correspond to the null or otherwise bad traces of the original seismic gather are substituted for the null traces of the seismic gather of FIG. 2. Proceeding to step 290, a determination is made whether the produced data traces are satisfactory for seismic interpretation. If the interpolated traces which replaced the original null traces are satisfactory, the seismic gather which include the newly interpolated seismic traces in place of the null traces are processed at step 300 by any one of the numerous data processing techniques to produce useful information regarding the characteristics of the subsurface formation explored. If the interpolated traces are not satisfactory, the method of the invention returns to step 260 for another propagation of the newly interpolated seismic gather to the fictitious datum followed by the return propagation of the next set of produced traces back to the original datum.

While numerous methods of determining whether a next iteration of the forward and backward datuming steps should be performed, one such method would be to determine whether the interpolated traces produced are coherent with the adjoining original seismic traces. One method of determining coherency would be to measure the amplitude of seismic events on the interpolated trace and the adjoining traces which would suggest to one skilled in the interpretation of seismic gathers to be suggestive of a coherent event and comparing the difference to a preselected reference value. If the difference is less than the preselected reference value, then the interpolated trace may be considered to be coherent and the interpolated trace satisfactory for seismic analysis. The preselected value which would give an accurate determination of coherency would be dependent on the characteristics of the particular seismic gather and may vary depending on factors such as the strength of the coherent event being examined or the size of spatial sampling between traces.

The methods of the present invention which were previously described by reference to the flow charts of FIGS. 1 and 4, while satisfactory, typically require extensive data processing time for the forward and backward datuming of the entire seismic gathers. In yet another embodiment of the present invention, such processing time is reduced by limiting the amount of the seismic sections for which propagation is required. We envision that processing time would be reduced by direct calculation of the contribution to the interpolated trace from each original, non-null trace. Direct calculation of the contribution to the interpolated trace from each original, nonnull trace, a method hereby defined as "redatuming" is determined according to the following equation: ##EQU3## where: .DELTA.x.sub.i is the separation between adjacent data trace locations for the input datum;

.DELTA.x.sub.j is the separation between adjacent data trace locations on the fictitious datum;

j is the location of the trace on the fictitious datum

R.sub.i (t-((r.sub.i -r.sub.j)/v)) is a representation of the original trace at the i.sup.th location to which a time shift of (r.sub.i -r.sub.j)/v has been applied and a derivative with respect to time has been taken;

f.sub.o is the center frequency of the input data traces;

X.sub.i (t) is the component of the interpolated trace which is produced by the ith one of the plurality of traces included as part of the original seismic section;

t is time;

r.sub.i is the distance between the original trace location and the jth location on the fictitious datum subject to contribution from the ith original trace;

r.sub.j is the distance between the jth location on the fictitious datum and the location of the interpolated trace;

.theta..sub.i is the angle between r.sub.i and the normal to the original datum;

.theta.j is the angle between r.sub.j and the normal to the fictitious datum; and

v is the velocity in the medium between the original and fictitious datums.

Equation (3) will produce the portion of the interpolated trace which is attributable to the information contained in a single seismic trace of the seismic gather. To produce the entire interpolated trace, the contribution of each seismic trace of the gather to the interpolated trace should be determined and the interpolation data attributable to all of the seismic traces combined.

Thus, there has been described and illustrated herein methods for interpolating missing seismic traces using a forward and backward application of wave equation datuming techniques. However those skilled in the art will recognize that many modifications and variations besides those specifically set forth may be made in the techniques described herein without departing substantially from the concept of the present invention. Accordingly, it should be clearly understood that the form of the invention described herein is exemplary only, and is not intended as a limitation on the scope of the claims.