Seismic migration method

Abstract

A seismic data processing method in which seismic traces (each trace containing a signal resulting from reflection of a seismic signal at a location within the earth, and each trace being associated with at least one point in a two-dimensional spatial grid (x,y)) are subjected to Fourier transformations, the coefficients of the Fourier-Transformed traces are subjected to a recursive KF migration operation, and the migrated traces are thereafter inverse-Fourier-transformed. When displayed, the processed seismic data accurately represents the position within the earth of whatever caused the reflection. The method may be employed to process stacked seismic traces, each associated with a single point (x,y) in the grid, or may be employed to process unstacked seismic traces, each associated with both a seismic source location (x.sub.s,y.sub.5) and a different seismic receiver location (x.sub.r,y.sub.r) in the grid. In performing the method, the earth is modeled as a stack of M horizontal layers each characterized by a seismic wave velocity. The recursive KF migration step is iterated M-1 times for each trace, where part of the output of each iteration is stored and part discarded.

Description

FIELD OF THE INVENTION

This invention is a method for processing seismic data. More particularly, the invention is an improved seismic migration technique of the type in which seismic data are Fourier-Transformed (or subjected to cosine or sine transformations), the coefficients of the transformed data are thereafter processed ("migrated"), and the processed coefficients are thereafter inverse-Fourier-transformed (or subjected to inverse-cosine or inverse-sine transformations).

BACKGROUND OF THE INVENTION

There are several types of conventional techniques for migrating seismic data. The objective of all migration techniques is to generate a display of processed seismic data in which the position of each seismic reflection indicated by the display accurately represents the position within the earth of whatever caused the reflection.

In one conventional class of seismic migration techniques, known as "KF-migration" (or "FK-migration") techniques, the seismic data to be processed are first Fourier-transformed, then the Fourier-coefficients are migrated, and thereafter, the migrated coefficients are inverse-Fourier-transformed. Examples of such KF-migration methods are described in the following papers: R. H. Stolt, "Migration by Fourier Transform," Geophysics, V. 43, No. 1 (February, 1978), pp. 23-48; J. Gazdag, "Wave Equation Migration with the Phase-Shift Method," Geophysics, V. 43, No. 7 (December, 1978), pp. 1342-1351; and J. Gazdag, "Wave Equation Migration with the Accurate Space Derivative Method," Geophysical Prospecting, V. 28, (1980), pp. 60-70.

A drawback of conventional KF-migration techniques of the type described in the Stolt paper is that they require that the relationship between depth of a reflector in the earth and the seismic wave traveltime to and from the reflector must be represented as a single fixed ratio. In other words, those methods require that the velocity of seismic waves in the earth be represented as a constant.

The Gazdag (1978) technique has the disadvantage that it is time-consuming and expensive to implement, for a number of reasons, including the reason that the number of operations required to migrate each seismic trace in accord with the Gazdag (1978) technique is proportional to N.sup.2, where N is the number of data points (or samples) comprising the trace. The Gazdag (1980) technique, too, has the disadvantage of being time-consuming and expensive to implement.

The present invention represents an improvement over conventional KF-migration techniques because it facilitates efficient migration of seismic data without requiring that the velocity of the associated seismic waves be represented as a constant. Rather, the invention requires only that the earth be represented as a sequence of uniform horizontal layers, each layer having a characteristic velocity, or ratio of depth increment to traveltime. The number of operations required to migrate a single trace in accord with the present invention is proportional to N log N, where N is the number of data points (or samples) comprising the trace.

SUMMARY OF THE INVENTION

The preferred embodiment of the invention includes the following nine steps:

1. Determining the reference time-depth relation;

2. Conforming each seismic trace of the input data to the reference time-depth relation;

3. Performing a Discrete Fourier Transform over the x-coordinate of each conformed trace;

4. Performing a Discrete Fourier Transform over the y-coordinate of each conformed trace;

5. Performing a Discrete Fourier Transform on the time coordinate of each conformed trace;

6. Performing the process known as "Recursive KF Migration" on the data after they have been processed in accordance with steps 1-5;

7. Performing a Discrete Inverse Fourier Transform over the y-coordinate of each migrated data trace;

8. Performing a Discrete Inverse Fourier Transform over the x-coordinate of each migrated data trace; and

9. After processing the data in steps 1-8, restoring the processed data's individual time-depth relations.

The earth is modeled as a stack of M horizontal layers, each having, in general, a different seismic wave velocity. The sixth step of the preferred embodiment ("recursive KF migration") is a sequence of operations performed M-1 times in succession (as an index j assumes values from 2 through M) for each transformed seismic trace. Recursive KF migration reassigns new values to the data points comprising the traces to approximate changes in the wave field as the associated source-receiver pair is "moved" downward into the earth. No interpolation is required in performing the recursive KF migration step of the invention. Rather, during each of the M-1 iterations for each trace, each data point of the trace is reassigned a single new value (that may be stored in a single computer storage location), and some of these new values are stored while the rest of these new values are discarded. The stored values remaining after all M-1 iterations are the result of applying the recursive KF migration operation to the relevant trace.

In alternative embodiments, steps 2 and 9 are omitted. The output traces produced in accord with any embodiment of the invention may be displayed in form meaningful to an interpreter using any of several suitable conventional methods and means.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram showing an unprocessed set of seismic traces, intermediate states of representative ones of these traces as they are processed in accord with the inventive method, and a set of traces in processed form after the invention has been applied thereto.

FIG. 2 is a diagram showing the iterations performed in performing the step of the invention known as "recursive KF migration".

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

The steps of the preferred embodiment will first be described in detail with reference to FIG. 1. Variations on this preferred embodiment will than be described. Thereafter, the sixth step of the preferred embodiment (known as "recursive KF migration") will be described in greater detail with reference to FIG. 2.

The seismic data to be processed in accordance with the preferred embodiment will be described with reference to Cartesian coordinate system (x,y,z), where the x-axis and y-axis are horizontal, and the z-axis is vertical. The seismic data are a set of seismic traces g(x,y,t.sub.i). There is one trace for each horizontal position (x,y). An example of a set of traces to be processed in accord with the invention is set A of traces shown in FIG. 1. In FIG. 1, set A includes twenty traces. It should be appreciated that a more typical set of traces to be processed in accord with the invention will be a 256 trace by 512 trace array of 131,072 traces, or an even larger array of traces. The parameters t.sub.i for each trace comprise a set of N time instants separated by sampling interal .DELTA.t. Each set t.sub.i may be described using the notation: {t.sub.i =i.DELTA.t:1.ltoreq.i.ltoreq.N}. It is within the scope of the invention to process data associated with only a single value of x (or y), or with many values of x and y (i.e., the method may be performed to process 2-D or 3-D data). For simplicity, it will be assumed for purposes of describing the preferred embodiment, that each trace represents the stacked electrical output of a number of acoustic receivers. For each instant t.sub.i, the outputs of the receivers are stacked together so that each stacked data point is associated with a single position (x,y) and a particular sample time t.sub.i, and trace g(x,y,t.sub.i) is a stacked seismic trace associated with single position (x,y). For a given trace, the sum of the receiver outputs at instant t.sub.i is represented by the numerical value g.sub.i =g(x,y,t.sub.i).

Despite the actual positions of the seismic source and receiver (or receivers) associated with each trace g(x,y,t.sub.i), the trace may be regarded as representing the output of an acoustic receiver at position (x,y) over an interval of time (0.ltoreq.t.ltoreq.N.DELTA.t) resulting from actuation of a seismic source at the same position (x,y) at t=0. It is presumed that each trace is indicative of at least one signal resulting from reflection by a subterranean reflector of the seismic wave proceeding outward from the source. It is further presumed that there exists a known or assumed relationship between the arrival time (t.sub.i) of each reflection at mesh point (x,y) and the depth (z.sub.i) of any real or hypothetical reflecting object directly below mesh point (x,y). The relationship between time and depth specified for each such regularly spaced time instant may be described using the notation {t.sub.i, z.sub.i :1.ltoreq.i.ltoreq.N}.

The first step of the preferred embodiment of the invention is performed as follows. For each time sample (t.sub.i), add together the depth values specified for each (x,y) mesh point and divide the sum by the number of addends. In other words, for each sample t.sub.i, generate ##EQU1## where the set of all mesh points (x,y) is {x.sub.mn, y.sub.mn :1.ltoreq.m.ltoreq.N.sub.x, 1.ltoreq.n.ltoreq.N.sub.y }, and z.sub.imn is the member of set z.sub.i associated with mesh point x.sub.mn, y.sub.mn. This produces the average time-depth relation, a set of depth values (z.sub.i.sup.av) corresponding one to each time sample of the seismic traces. From the resulting set of time-depth pairs {t.sub.i =t.sub.i.sup.av, z.sub.i.sup.av :1.ltoreq.i.ltoreq.N} select the reduced set {t.sub.j.sup.ref, z.sub.j.sup.ref :1.ltoreq.j.ltoreq.M, M<N}, constituting the reference time-depth relation, to satisfy the rules t.sub.1 =t.sub.1.sup.ref =t.sub.1.sup.ref =z.sub.1.sup.av, .vertline.t.sub.i.sup.av -[t.sub.j-i.sup.ref +(z.sub.i.sup.av -z.sub.j-1.sup.ref)(t.sub.j.sup.ref -t.sub.j-1.sup.ref)/(z.sub.j.sup.ref -z.sub.j-1.sup.ref)].vertline.<E, where t.sub.j-1.sup.ref <t.sub.i.sup.av .ltoreq.t.sub.j.sup.ref, and E is a selected acceptable tolerance. The feasibility of the invention proceeds from the observation that M may be of order 10 while N is of order 1000, if E is of order 10.DELTA.t.

Each seismic trace of the input data is associated with its individual time-depth relationship {t.sub.i, z.sub.i :1.ltoreq.i.ltoreq.N} defining the time and depth for each sample of the trace. The second step of the preferred embodiment is conforming the input seismic traces to the reference time-depth relationship. This second step includes the sub-step of assigning each sample a modified time t'.sub.i according to the scheme

t'.sub.i =t.sub.j-1.sup.ref +(t.sub.j.sup.ref -t.sub.j-1.sup.ref)(z.sub.i -z.sub.j-1.sup.ref)/(z.sub.j.sup.ref -z.sub.j-1.sup.ref),

where z.sub.j-1.sup.ref .ltoreq.z.sub.i .ltoreq.z.sub.j.sup.ref. The resulting seismic trace {t'.sub.i, g.sub.i } is not sampled at regular time intervals. The second sub-step of the second inventive step is to construct from each trace {t'.sub.i, g.sub.i } the regularly sampled conformed trace {(i-1).DELTA.t, f.sub.i } according to the scheme f.sub.i =g.sub.j-1 +(g.sub.j -g.sub.j-1) ((i-1).DELTA.t-t'.sub.j-1)/(t'.sub.j -t'.sub.j-1), where t'.sub.j-1 .ltoreq.(i-1).DELTA.t.ltoreq.t'.sub.j.

The set of conformed input traces {(i-1).DELTA.t, f.sub.i } may be represented as f[(m-1).DELTA.x, (n-1).DELTA.y, (i-1).DELTA.t]; where 1.ltoreq.m.ltoreq.N.sub.x, 1.ltoreq.n.ltoreq.N.sub.y, 1.ltoreq.i.ltoreq.N. The third step of the preferred embodiment is to subject this set of N.sub.x N.sub.y N numbers to a Discrete Fourier Transform on index m to produce set of values e[(m-1).DELTA.k.sub.x, (n-1).DELTA.y, (i-1).DELTA.t], where .DELTA.k.sub.x =1/(N.sub.x .DELTA.x) and 1.ltoreq.m.ltoreq.N.sub.x.

In the fourth step of the preferred embodiment, the set of signals e is subjected to a Discrete Fourier Transform on index n to become set of signals d[m-1).DELTA.k.sub.x, (n-1).DELTA.k.sub.y, (i-1).DELTA.t], where .DELTA.k.sub.y =1/(N.sub.y .DELTA.y) and 1.ltoreq.n.ltoreq.N.sub.y.

In the fifth step of the preferred embodiment, the set of signals d is subjected to a Discrete Fourier Transform on index i to become set of signals c[(m-1).DELTA.k.sub.x, (n-1).DELTA.k.sub.y, (i-1).DELTA.w] where .DELTA.w=1/(N .DELTA.t) and 1.ltoreq.i.ltoreq.N. The Discrete Fourier Transform in any of steps 3, 4, and 5 may be performed using any of several suitable conventional implementations thereof. See, for example, G. D. Bergland, "A Guided Tour of the Fast Fourier Transform," IEEE Spectrum, 6, pp 41-52 (July, 1969) for a discussion of the Discrete Fast Fourier Transform.

The steps two, three and four are represented in FIG. 1 by generation of all the traces d[(m-1).DELTA.k.sub.x, (n-1).DELTA.k.sub.y, (i-1).DELTA.t] represented by single trace B from set A. Step five is represented in FIG. 1 by transformation of all traces represented by trace B into set of transformed traces c[(m-1).DELTA.k.sub.x, (n-1).DELTA.k.sub.y, (i-1).DELTA.w] represented by single signal C.

The set of signals represented by signal C is next transformed into the set of numerical signals a[(m-1).DELTA.k.sub.x, (n-1).DELTA.k.sub.y, (i-1).DELTA.t] according to the "Recursive KF Migration" method described below. This recursive KF migration method (step six of the preferred embodiment) is a key element of all embodiments of the invention, and is represented in FIG. 1 by transformation of single signal C into signal D. In the FIG. 1 example, should it be understood that during step six, the recursive KF migration operation will be repeated to process each of the signals (all represented by signal C) resulting from performance of step five.

In the seventh step of the preferred embodiment, the set of signals a is subjected to a Discrete Inverse Fourier Transform on index n to become set of signals h[(m-1).DELTA.k.sub.x, (n-1).DELTA.y, (i-1).DELTA.t]. The Discrete Inverse Fourier Transform of this step seven (and of step eight below) may be performed using any of several suitable conventional implementations thereof.

In the eighth step of the preferred embodiment, the set of numerical signals h is subjected to a Discrete Inverse Fourier Transform on index m to become set of signals q[(m-1).DELTA.x, (n-1).DELTA.y, (i-1).DELTA.t]. The resulting set q of numerical values is a set of migrated conformed output traces. To select a pair of particular values (m,n) is to specify a particular migrated conformed output trace {(i-1).DELTA.t, q.sub.i :1.ltoreq.i.ltoreq.N} associated with the time-depth relationship at mesh point ((m-1).DELTA.x, (n-1).DELTA.y).

In the ninth step of the preferred embodiment, to compensate for the effect of step two described above, each sample is assigned a modified time t.sub.i according to the scheme t.sub.i =t.sub.j-1 +(t.sub.j -t.sub.j-1)(z'.sub.i -z.sub.j-1)/(z.sub.j -z.sub.j-1) where z.sub.j-1 <z'.sub.i .ltoreq.z.sub.j ; and {t.sub.j, z.sub.j :1.ltoreq.j.ltoreq.N} is the time-depth relationship at the mesh point; and z'.sub.i =z.sub.j-1.sup.ref +(z.sub.j.sup.ref -z.sub.j-1.sup.ref)((i-1).DELTA.t-t.sub.j-1.sup.ref)/(t.sub.j.sup.ref -t.sub.j-1.sup.ref), with t.sub.j-1.sup.ref <(i-1).DELTA.t.ltoreq.t.sub.j.sup.ref. The resulting migrated trace (t.sub.i, q.sub.i) is not sampled at regular time intervals. From it is constructed the regularly sampled migrated output trace {(i-1).DELTA.t, p.sub.i } according to the scheme p.sub.i =q.sub.j-1 +(q.sub.j -q.sub.j-1)((i-1).DELTA.t-t.sub.j-1)/(t.sub.j -t.sub.j-1). Steps seven, eight, and nine are represented in FIG. 1 by transformation of the set of signals represented by single signal D into set E, and will be carried out to process all the signals produced as a result of step six.

Performance of all nine steps of the preferred embodiment of the invention results in generation of a migrated output trace {(i-1).DELTA.t,p.sub.i :1.ltoreq.i.ltoreq.N} for each of the mesh points (x,y) associated with an input trace. An example of a set of migrated output traces generated in accord with the invention is set E of traces shown in FIG. 1. The output traces produced in accord with the invention may be displayed in form meaningful to an interpreter using any of several suitable and well-known display methods and means.

The data processing steps described above should be translated into a series of specific computer instructions in a manner that will be apparent to those ordinarily skilled in the art of computer programming for seismic data processing. The method steps preferably should be performed using any of many commercially available, appropriately programmed digital computers.

It will be apparent to those of ordinary skill in the field of seismic data processing that numerous variations on the above-described preferred embodiment may be implemented. The inventors consider that it is still within the scope of the invention to vary the above-described method as follows:

(a) Step 5 precedes set 4 or step 3;

(b) Steps 3 and 4 are interchanged;

(c) Steps 7 and 8 are interchanged;

(d) Steps 2 and 9 are omitted;

(e) Nonlinear interpolation is performed where, for simplicity, linear interpolation has been indicated; or

(f) Discrete Cosine Transforms or Discrete Sine Transforms are employed where Discrete Fourier Transforms have been indicated (any of several suitable conventional implementations of the Discrete Cosine Transforms or Discrete Sine Transform may be employed).

Other variations may be within the scope of the invention. Throughout this application (including in the claims), the phrase "Discrete Fourier Transform" (or "Discrete Inverse Fourier Transform") should be understood to denote any discrete transform selected from the group including Discrete Fourier Transforms, Discrete Cosine Transforms, and Discrete Sine Transforms (and the corresponding discrete inverse transformation).

The sixth step of the preferred embodiment ("Recursive KF Migration") will next be described. In this description, the symbol W.sub.N.sup.ij will represent the complex number cos Q+.sqroot.-1 sin Q, where Q=2.pi.ij/N. The reference time-depth relation is defined by the set of numbers {t.sub.j.sup.ref, z.sub.j.sup.ref :1.ltoreq.j.ltoreq.M, M.gtoreq.2}. The recursive KF migration step is a sequence of operations carried out M-1 times in succession, as index j assumes values from 2 through M. This process is carried out separately for every possible combination of m and n. Each performance of recursive KF migration on a trace (identified by index pair m,n) reassigns a new value to each sample comprising the trace, to approximate changes in the seismic wave field as the source-receiver pair (or pairs) associated with index pair n,n are "moved" downward into the earth. The earth is modeled as a stack of M horizontal layers, each having, in general, a different seismic wave velocity, v.sub.j (to be defined precisely below). FIG. 2 shows the iteration cycles performed during application of recursive KF migration to a set of traces.

The first operation is known as downward continuation. In the first execution, when j=2, this sub-step is bypassed. In subsequent executions, each number c[m-1).DELTA.k.sub.x, (n-1).DELTA.k.sub.y, (i-1).DELTA.w] as it currently exists is multiplied by the corresponding factor W.sub.N.sup.pq where q=Int[(t.sub.j-1.sup.ref -t.sub.j-2.sup.ref)/.DELTA.t], understanding that the symbol Int[ ] indicates that the expression within the brackets is to be rounded to the nearest integer, and p=[(-1).sup.2 -k.sup.2 V.sub.j-1 ].sup.1/2, V.sub.j-1 =(z.sub.j-1.sup.ref -z.sub.j-2.sup.ref)/(t.sub.j-1.sup.ref -t.sub.j-2.sup.ref), k.sup.2 =(m-1).sup.2 (.DELTA.k.sub.x /.DELTA.w).sup.2 +(n-1).sup.2 (.DELTA.k.sub.y /.DELTA.w).sup.2, provided that i.gtoreq.1+kV.sub.j-1 ; otherwise, c[(m- 1).DELTA.k.sub.x, (n-1).DELTA.k.sub.y, (i-1).DELTA.w]=0.

The next operation is known as relinearization. For each sample index i, a quantity i' is computed according to i'Int[1+{(i-1).sup.2 -k.sup.2 V.sub.j.sup.2 .sup.1/2 ], 1+kV.sub.j.ltoreq. i.ltoreq.N/2+1, or i'=Int[N+1-{(N+1-i).sup.2 -k.sup.2 V.sub.j.sup.2 }.sup.1/2 ], N/2+1<i.ltoreq.N-kV.sub.j. The objective is to define c'[(m-1).DELTA.k.sub.x, (n-1).DELTA.k.sub.y, (i'-1).DELTA.w]=c[(m-1).DELTA.k.sub.x,(n-1).DELTA.k.sub.y, (i-1).DELTA.w]. If some particular integer value i' is not among those defined in this manner, the corresponding value of c' is taken to be zero. Various interpolation methods are permissible to distribute any sample of c among several samples of c' . The restriction that must be observed is that ##EQU2##

In the next operation, the set of numbers c' is subjected to a Discrete Inverse Fourier Transform on index i' to become b[(m-1).DELTA.k.sub.x,(n-1).DELTA.k.sub.y,(i'-1).DELTA.t], 1.ltoreq.i'.ltoreq.N.

In the next operation, one stores the partial result: a[(m-1).DELTA.k.sub.x,(n-1).DELTA.k.sub.y,(i-1).DELTA.t]=b[(m-1).DELTA.k.s ub.x, (n-1).DELTA.k.sub.y,(i'-1).DELTA.t], 1+Int[t.sub.j-1.sup.ref /.DELTA.t].ltoreq.i<1+Int[t.sub.j.sup.ref /.DELTA.t]. Thus, for each iteration associated with a particular value of j, some but not all of the values b[(m-1).DELTA.k.sub.x,(n-1).DELTA.k.sub.y, (i-1).DELTA.t]] are stored.

With reference to FIG. 2, the order in which different values of m and n are employed in performing recursive KF migration is immaterial. It is necessary that each possible combination of m and n be employed once and only once.

In the above description of the recursive KF migration step, let V.sub.max be defined as the maximum of the ratio (z.sub.j.sup.ref -z.sub.j-1.sup.ref)/(t.sub.j.sup.ref -t.sub.j-1.sup.ref), considering all permissible values of j. Before carrying out the recursive KF migration operation, it may be useful to set c[(m-1).DELTA.k.sub.x, (n-1).DELTA.k.sub.y,(i-1).DELTA.w]=0 for 1.ltoreq.i.ltoreq.kV.sub.max /sin d.sub.mig, where k.sup.2 =(m-1).sup.2 (.DELTA.k.sub.x /.DELTA.w).sup.2, d.sub.mig is arbitrary but restricted by the requirement 0< sin d.sub.mig .ltoreq.1. It may furthermore be useful, before thus modifying the set of numbers c, to generate and set aside the set of numerical signals

c.sub.save (i)=c[(m-1).DELTA.k.sub.x,(n-1).DELTA.k.sub.y,(i-1).DELTA.w]

for (kV.sub.max /tan d.sub.noise).ltoreq.i-1.ltoreq.(kV.sub.max /tan d.sub.mig) where d.sub.noise is arbitrary but restricted by the requirement tan d.sub.mig .ltoreq.tan d.sub.noise. Then, during recursive KF migration, each time c' is computed, c.sub.save is added to c', sample by sample, over the range of samples for which c.sub.save is defined.

Numerous other variations on the inventive embodiment described above are with the scope of the invention. For example, the seismic data to be processed in accord with the invention may be stacked or unstacked. The preferred embodiment described above is suitable for processing zero-offset (i.e., stacked) data. With unstacked data, each input trace g(x.sub.s, y.sub.s, x.sub.R, y.sub.R, t.sub.i) will be associated with both a source position (x.sub.s, y.sub.s) and a receiver position (x.sub.R, y.sub.R); not with a single position (x,y) representing both source and receiver positions. Accordingly, where unstacked data are to be processed in accordance with the invention, modifications to the above-described must be made in a manner that will be apparent to one of ordinary skill in the seismic data processing art. In particular, the relations employed for downward continuation and relinearization will change to reflect what is known to those of ordinary skill in the art as the "Double Square Root" equation.

The above description is merely illustrative of the invention. It is contemplated that various changes in the details of the structures and methods described may be within the scope of the invention as defined by the appended claims.

Claims

1. A method for processing a set of seismic traces, each trace corresponding to at least one location in a two-dimensional spatial grid (x,y), and each trace including a set of N time samples corresponding to time instants {t.sub.i =i.DELTA.t:1.ltoreq.i.ltoreq.N}, given a time-depth relationship {t.sub.i,z.sub.i:1.ltoreq.i.ltoreq.N} for each trace, said method including the steps of:

(a) determining a reduced set of time-depth pairs {t.sub.j.sup.ref,z.sub.j.sup.ref:1.ltoreq.j.ltoreq.M, M<N}, where t.sub.1 =t.sub.1.sup.av =t.sub.1.sup.ref, z.sub.1.sup.ref =z.sub.1.sup.av,.vertline.t.sub.i.sup.av -[t.sub.j-1.sup.ref +(z.sub.i.sup.av -z.sub.j-1.sup.ref)(t.sub.j.sup.ref -t.sub.j-1.sup.ref)/(z.sub.j.sup.ref -z.sub.j-1.sup.ref)].vertline.<E, where t.sub.j-1.sup.ref <t.sub.i.sup.av.ltoreq.t.sub.j.sup.ref, z.sub.i.sup.av is determined by adding together the depth values z.sub.i for each trace and dividing by the number of traces, and E is a selected acceptable tolerance on the order of about 10.DELTA.t;

(b) subjecting each trace to Discrete Fourier Transforms on the spatial coordinates corresponding to the trace, and on time index i;

(c) performing recursive KF migration on the traces after they have been transformed in step (b); and

(d) subjecting each migrated trace generated in step (c) to Discrete Inverse Fourier Transforms on the spatial coordinates corresponding thereto.

2. The method of claim 1, where each trace is a stacked trace g(x,y,t.sub.i), corresponding to a single point (x.sub.mn, y.sub.mn) in grid (x,y), and (x.sub.mn, y.sub.mn) is determined by the position coordinates of the seismic sources and receivers associated with the individual traces comprising the stacked trace.

3. The method of claim 1, where each trace is an unstacked trace g(x.sub.s, y.sub.s, x.sub.R, y.sub.R, t.sub.i) corresponding to a seismic source position (x.sub.s, y.sub.s) and to an acoustic receiver position (x.sub.R, y.sub.R).