BACKGROUND AND OBJECTS OF THE INVENTION

A wavefield extrapolator simulates the effect of wave propagation. From [estimated) recordings at one depth level the extrapolator computes estimates of recordings which would have been recorded at another depth level. Recursive wavefield extrapolation forms the basis of any modern seismic migration technique and also plays an important role in seismic forward modelling.

The most important input information for migration techniques consists of a velocity model of the subsurface. A detailed velocity model normally cannot be constructed before migration has been performed. An iterative interpretive procedure where the velocity model is successively updated in order to improve the migrated result is therefore highly desirable.

Wavefield extrapolation is an. extremely demanding computational task. The introduction of high performance vector processors like the Cray-1 and the Cyber 205 has at least made such computations a practical reality. However, the run time on such computers for realistic seismic problems is far too long to make iterative techniques feasible on routine basis.

The emergence of Very Large Scale Integration technology and automated design tools have made highly parallel computing structures feasible. Attempts at producing highly parallel general purpose computers have not yet been successfull. The performance of such computers is highly algorithm dependent and existing and proposed systems require special non-portable code. Today, there is no known way to meet extreme computational requirements other than to implement dedicated highly parallel systems.

A fast wavefield extrapolator would make iterative and interactive migration schemes feasible. A related potential application for a fast wavefield extrapolator is as the modelling part of a Kalman filter for estimation of the parameters in a geological depth model.

SEISMIC MIGRATION IN THE SPACErFREQUENCY DOMAIN

In seismic migration of zero-offset data the upcoming wavefield recorded at the surface is recursively downward inverse extrapolated. At each depth level imaging is performed by computing the pressure for zero time.

For simplicity, this note is limited to the treatment of migration of zero-offset data using convolutional operators (Kirchoff summation techniques and explicit finite difference techniques). However, as long as there is no coupling between different temporal frequencies, the proposed overall structure is completely general.

The total procedure of migration of zero-offset data in the space-frequency domain can be summarized as follows (Berkhout, 1982):

1. Transform each seismic trace to the frequency domain:

^' p(x .,y ,z=0,t) ^■ * p.x_.,y. ,z=0,u.) j=1,

J ^K J ^K ι_--1

2. ^' Convolve the data with the inverse extrapolator W for each frequency component ω : plx.y.z.iu.) ^■ * p(x,y,z+Δz,u). ) i=1,.... p.x.y.z+Δz.ω..) = p(x,y,z,ω.) * W(c(x,y,z))

3. Image the inverse extrapolated result at each depth level: p(x ,y ,z+Δz,t=0) = 1/ττ E p(x ., ,z+Δz.ω)

Step 1 is applied only once. The problem of performing fast Fouriertransformation will not be addressed here. We shall assume that the data are availible in the space-frequency domain or that the FFT can be performed as the data are read from mass storage.

Step 2 and 3 are executed repeatedly. In practical applications the total number of depth steps can be expected to range from a few hundred to about one thousand. Note that for each extrapolation step the complete data volume must be retrieved and stored.

The extrapolation can be performed independently for all ω. For typical seismic problems the number of frequencies that should be

Lateral velocity variations can be coped with by using space variant convolutions. The operators can be precomputed and accessed by table » look-up.

OVERALL DESCRIPTION OF THE COMPUTATIONAL STRUCTURE

After a trivial reformulation the central part (step 2 and 3) of the algorithm described above appears as follows:

where N is the number of temporal frequencies that is to be processed concurrently.

A possible structure for performing the central part of the algorithm is shown in fig.1. One monochromatic wavefield extrapolator (MWE) and one memory module is assigned to each temporal frequency in the frequency band that is to be processed concurrently. The private memories of each MWE are loaded with data from a global data bus. During extrapolation and (partial) imaging the desired number of frequencies are extrapolated simultaneously by operating the individual MWE's in lock-step. We shall assume that each MWE works like a systolic array; data from one depth level flows into the MWE while data extrapolated to the next depth level flows out. Together with the data the MWE's must be fed with velocity information. Imaging for t=0 is performed by simultaneous summation of part of the output from each MWE.

Although conceptually simple, this structure has some serious disadvantages:

The adder must be implemented as a pipelined adder tree with length equal to the logarithm of the number of MWE's.

_» For a high number of MWE's all global communication must be implemented as (pipelined) tree structures. This problem becomes increasingly important with decreasing cycle time.

Because of the tree structures, an eventual increase in the number of MWE's is a nontrivial task.

The problems mentioned above are all circumvented in the following proposed design, displayed in fig.2. Here the need for for global communication is totally eliminated. Each MWE is physically connected only to its two neighbours. All control signals, including system clock, are fed from the top of the structure. For each clock cycle these signals propagate one MWE downwards. During loading, the input data and their destination addresses are pumped systolically down through all the memories. Each private memory accepts a write only if the data address ^' lies within the address field of that particular memory. During extrapolation the structure is operated as a two level systolic array. Seismic velocity and data address information flow vertically downward. In each MWE velocity information and wavefield data flow horizontally. Imaging is performed by recursive summation. Partially imaged results are pumped vertically downwards and at each MWE the contribution from that particular frequency is added. Imaged (or partially imaged) results are availible at the bottom. Imaging for other values of t than zero can be performed by weighting the summations to yield partial inverse DFT's. Concurrent imaging for different values of t might be performed by adding more "imaging columns" .

Extension of the proposed structure to any desired size is trivial. The number of MWE's that can be added is limited only by the total number of temporal frequencies that is to be taken into account.

The approach sketched here has the additional advantage that the bandwidth of the velocity information can be progressively shrinked as the velocity information propagates down the structure. This can be obtained by letting each MWE act as part of a cascade filter.

Because of the systolic vertical distribution of clock signals, there

will be a clock skew proportional to the number of MWE's between the control processor and the bottom of the array. A synchronization unit might therefore be required to transmit the output results to their destinations.

Because of the potentially high number of processors involved, the computational structure must be able to perform its task even when not all the MWE's are operational. A defect MWE must therefore still be able to propagate the information it receives further downward. Each MWE should thus be equipped with a self testing system that is able to activate a bypass function if an error has occured.

ARCHITECTURES FOR MONOCHROMATIC WAVEFIELD EXTRAP0LAT0RS

A systolic array for space variant 1-D convolution with fixed operator length (Kung, 19Θ2) is displayed in fig.3. Weights W are preloaded into the array, one set for each cell. Both partial result ^' s p(x,y,z+Δz) and inputs p(x,y,z) flow from left to right, but the inputs move twice as fast as the partial results. More precisely, each partial result stays inside every cell it passes for one cycle, thus it takes twice as long to march through the array as an inputs does. Each partial result is initialized to zero before entering the leftmost cell, accumulating all its terms while marching to the right. In each cell the weight is addressed by the velocity information tied to the input.

For symmetric operators roughly half the multipliers can be saved by adding data points located symmetrically about the operator midpoint before rather than after multiplication. In this case it is advantageous to .use an architecture as displayed in fig.4. This design has two independent input streams. The lower input stream moves three times as fast as the upper input stream and twice as fast as the output stream. The input arguments are added prior to multiplication with the operator weight. The operator weights are addressed by the velocity information that is now tied to the output stream.

There are many ways of making the above structures capable of performing convolution with operators that are larger in size and dimension than the array. A simple approach is shown in fig.5. In the first pass, a convolution is performed with one part of the complete

operator. The results are stored in a cache memory. In the next pass, another part of the operator is addressed and the results of the previous pass, rather than zeroes, are used as starting values. This continues till the whole operator length in one dimension is covered. Then we start with a new row of the 2-D data and use our partial results as starting values. In this way the contributions from all parts of the operator are accumulated till the complete operator has covered the data. In'the last pass the extrapolated data are output.

In the case that the 2-D operator matrix has rank one, or approximately so, the 2-D convolution can be carried out as two subsequent 1-D convolutions in orthogonal directions.

CONCLUSIONS AND DISCUSSIONS

A systolic computer architecture for ultra fast seismic migration by recursive wavefield extrapolation is proposed. The structure is decomposed into a basic building block which is used repetitively with simple interfaces. This is expexted to significantly simplify development, testing and production. There are a number of possible variations on the proposed concept, involving a variety of tradeoffs. In the most general case each MWE could be implemented as a general computer and spatial information about a number of physical properties could be used as input.

The concept also has obvious applications to areas that strictly speaking are not embraced by recursive wavefield extrapolation, for instance spatial data filtering and pre stack partial migration. In terms of wavefield extrapolation, as the word is used here, these operations can be performed as one or a few wavefield extrapolations.

Using the concept proposed here, a large amount of the total potential parallelism can be utilized by concurrent processing of any number of temporal frequencies. An alternative, or additional, way of introducing concurrency would be to parallelize in space. However, any such approach would inherently suffer from the serious drawback that velocity information for a number of different spatial coordinates would have to be supplied in parallel.

REFERENCES

Berkhout, A.J., 1982, Seismic migration. Imaging of acoustic energy by wavefield extrapolation. A. Theoretical aspects. Elsevier Scientific Publ. Co. Amsterdam.

Kung, H.T., 1982, Why systolic architectures? IEEE Computer, Jan 82,

FIGURE CAPTIONS

Figure 1. A possible computing structure for the migration of one frequency band.

Figure 2. Proposed computing structure for the migration of one frequency band. By distributing all control signals systolically, the need for global communication is totally eliminated.

Figure 3. A systolic array for band matrix - vector multiplication.

Figure 4. A systolic array for space variant convolution with symmetric operators.

Figure 5. ^Extension to make the linear systolic arrays capable of performing convolution with operators that are larger in size and dimension than the arrays.