NO843822L

NO843822L - CALCULATOR INSTRUCTION FOR EXTRAPOLATION OF HOUSE FIELD.

Info

Publication number: NO843822L
Application number: NO843822A
Authority: NO
Inventors: Olav Holberg
Original assignee: Olav Holberg
Priority date: 1984-09-24
Filing date: 1984-09-24
Publication date: 1986-07-08
Also published as: WO1986001918A1

Abstract

Den foreliggende oppfinnelse vedrører datamaskinarkitekturen til en enhet for hurtig ekstrapolasjon av bølgefelt som kan anvendes blant annet til seismisk migrasjon. Bølgefeltekstrapolasjonen foregår i rom-frekvens domenet fordi dette domenet tillater uavtiengig og parallell behandling av alle frekvenser og fordi dette domenet er det mest fleksible for ekstrapolasjon av bølgefelt. Ved å distribuere alle kontrollsignaler systolisk er behovet for global kommunikasjon fullstendig eliminert. Den foreslåtte arkitektur kan lett utvides til parallell prosessering av et hvilket som helst antall frekvenser. For den klasse problemer som kan formuleres som bandmatrise-vektor multiplikasjon eller rom-variant konvolusjon, kan den foreslåtte struktur dekomponers til en enkelt multipliser-adder celle som benyttes repetitivt med enkle grensesnitt.The present invention relates to the computer architecture of a unit for rapid extrapolation of wave fields which may be used inter alia for seismic migration. The wave field extrapolation takes place in the space-frequency domain because this domain allows for incomplete and parallel processing of all frequencies and because this domain is the most flexible for extrapolation of wave fields. By distributing all control signals systolic, the need for global communication is completely eliminated. The proposed architecture can easily be extended to parallel processing of any number of frequencies. For the class of problems that can be formulated as band matrix-vector multiplication or space-variant convolution, the proposed structure can be decomposed into a single multiplier-adder cell which is used repetitively with simple interfaces.

Description

BESKRIVELSE AV REGNEMASKINSTRUKTUR FOR EKSTRAPOLASJON AV BØLGEFELTDESCRIPTION OF CALCULATOR STRUCTURE FOR EXTRAPOLATION OF WAVE FIELD

INNLEDNINGINTRODUCTION

En bølgefeltekstrapolator simulerer effekten av bølgeforplantning.A wavefield extrapolator simulates the effects of wave propagation.

Fra (estimerte) data på ett nivå beregner ekstrapolatoren data på et annet nivå. Rekursiv ekstrapolasjon av bølgefelt danner grunnlaget for alle moderne algoritmer for seismisk migrasjon og spiller også en viktig rolle i seismisk modellering. From (estimated) data on one level, the extrapolator calculates data on another level. Recursive wavefield extrapolation forms the basis of all modern seismic migration algorithms and also plays an important role in seismic modeling.

Den viktigste inn-informas jon for seismisk migrasjon består av en hastighetsmodell for mediet. En detaljert hastighetsmodell kan normalt ikke konstrueres før migrasjon har vært utført. En iterativ tolkningsmessig prosedyre hvor hastighetsmodellen oppdateres suksessivt for å forbedre det migrerte resultat er derfor ønskelig. The most important input information for seismic migration consists of a velocity model for the medium. A detailed velocity model cannot normally be constructed until migration has been carried out. An iterative interpretation procedure where the velocity model is successively updated to improve the migrated result is therefore desirable.

Bølgefeltekstrapolasjon er en ekstremt regnekrevende oppgave. Introduksjon av vektormaskiner som Cray-1 og Cyber 205 har gjort slike beregninger en praktisk realitet. Imidlertid er regnetiden på slike maskiner for realistiske seismiske problemer altfor lang til at iterative teknikker kan brukes på rutinebasis. Wave field extrapolation is an extremely computationally demanding task. Introduction of vector machines such as the Cray-1 and Cyber 205 have made such calculations a practical reality. However, the computation time on such machines for realistic seismic problems is far too long for iterative techniques to be used on a routine basis.

VLSI (Very Large Scale Integration) teknikk har muliggjort realisering av høy-parallelle regnemaskinstrukturer. Forsøk på å frembringe høy-parallelle generelle regnemaskiner har hittil ikke vært vellykkede. Ytelsen til slike maskiner er svært algoritmeavhengig og eksisterende og foreslåtte systemer krever spesiell ikke-portabel programkode. Idag eksisterer der intet kjent alternativ til å implementere dedikerte høy-parallelle systemer for å tilfredstille ekstreme behov for regnekraft. VLSI (Very Large Scale Integration) technology has enabled the realization of highly parallel computing structures. Attempts to produce highly parallel general purpose computers have so far been unsuccessful. The performance of such machines is highly algorithm dependent and existing and proposed systems require special non-portable program code. Today, there is no known alternative to implementing dedicated highly parallel systems to satisfy extreme needs for computing power.

En hurtig bølgefeltekstrapolator vil muliggjøre iterative og interaktive migrasjonsprosedyrer. En relatert potensiell anvendelse av en hurtig bølgefeltsekstrapolator er som modelleringsdelen av et Kalman filter for estimering av parametere som karakteriserer en geologisk modell. A fast wavefield extrapolator will enable iterative and interactive migration procedures. A related potential application of a fast wavefield extrapolator is as the modeling part of a Kalman filter for estimating parameters characterizing a geological model.

SEISMISK MIGRASJON I ROM-FREKVENS DOMENETSEISMIC MIGRATION IN THE SPACE-FREQUENCY DOMAIN

I seismisk migrasjon av null-offset (stackede) data blir bølgefeltet som er målt på overflaten rekursivt ekstrapolert nedover i dyp. På hvert dybdenivå utføres en avbilding ved at man beregner amplitude for tid lik null. In seismic migration of zero-offset (stacked) data, the wavefield measured at the surface is recursively extrapolated downward at depth. At each depth level, an image is performed by calculating the amplitude for time equal to zero.

For ikke å gjøre fremstillingen unødig komplisert er dette notat begrenset til behandling av migrasjon av null-offset data med konvolusjonsoperatorer for ekstrapolasjon av akustiske bølger (Kirchoff summasjonsteknikker og eksplisitte endelig-differens teknikker). Imidlertid skal det understrekes at den foreslåtte struktur er fullstendig generell sålenge det ikke er noen direkte kopling mellom forskjellige frekvenser. In order not to make the presentation unnecessarily complicated, this note is limited to the treatment of migration of zero-offset data with convolution operators for the extrapolation of acoustic waves (Kirchoff summation techniques and explicit finite-difference techniques). However, it should be emphasized that the proposed structure is completely general as long as there is no direct connection between different frequencies.

Migrasjon av null-offset data i rom-frekvens domenet kan summarisk fremstilles som følger (Berkhout, 1982): 1. Transformer hver seismisk trase til frekvens domenet: 2. Konvolver dataene med ekstrapolatoren W for hver frekvens-komponent uj : 3. Avbilde det ekstrapolerte resultat på hvert nivå: 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: 2. Convolve the data with the extrapolator W for each frequency component uj : 3. Image it extrapolated results at each level:

Skritt 1 anvendes kun en gang. Vi skal ant at dataene er tilgjengelige i rom-frekvens domenet eller at FFT (Fast Fourier Transform) kan utføres mens dataene leses fra masselager. Step 1 is used only once. We shall assume that the data is available in the space-frequency domain or that the FFT (Fast Fourier Transform) can be performed while the data is being read from mass storage.

Skritt 2 og 3 utføres repetitivt. I praktiske anvendelser kan antall dybdeskritt forventes å variere fra noen få hundre til omtrent ett tusen. Merk at for hvert ekstrapolasjonsskritt må hele datavolumet leses og skrives. Steps 2 and 3 are performed repetitively. In practical applications, the number of depth steps can be expected to vary from a few hundred to about one thousand. Note that for each extrapolation step the entire data volume must be read and written.

Ekstrapolasjonen kan utføres uavhengig for alle w. For typiske seismiske problemer er antallet frekvenser som skal behandles størrelsesorden 300. The extrapolation can be performed independently for all w. For typical seismic problems, the number of frequencies to be treated is of the order of 300.

Laterale hastighetsvariasjoner kan handteres ved å anvende rom-variante konvolusjonsoperatorer. Operatorene kan aksesseres med tabelloppslag. Lateral velocity variations can be handled by applying space-variant convolution operators. The operators can be accessed with a table lookup.

OVERORDNET BESKRIVELSE AV REGNEMASKINSTRUKTURENGENERAL DESCRIPTION OF THE CALCULATOR STRUCTURE

Efter en triviell reformulering fremkommer den sentrale del (skritt 2 og 3) av algoritmen beskrevet ovenfor som følger: After a trivial reformulation, the central part (steps 2 and 3) of the algorithm described above appears as follows:

hvor N er antall frekvenser som skal prosesseres samtidig. where N is the number of frequencies to be processed simultaneously.

En mulig struktur for å utføre den sentrale del av algoritmen er vist i fig.1. En monokromatisk bølgefeltsekstrapolator (MBE) og en hukommelsesmodul er tilordnet til hver frekvens i det frekvensbånd som skal behandles samtidig. De private hukommelsesmodulene til hver MBE lastes med data fra en global data buss. Under ekstrapolasjon og (partiell) avbilding, ekstrapoleres det ønskede antall frekvenser samtidig ved å operere de individuelle MBE'er simultant. Vi skal anta at hver MBE virker som et systolisk ar/ray; data fra et dybdenivå pumpes inn mens data ekstrapolert til neste dybdenivå pumpes ut. Sammen med data må hver MBE fores med romlig hastighetsinformas jon. Avbilding for tid lik null utføres ved simultan summasjon av realdelen av resultatet fra hver MBE. A possible structure for executing the central part of the algorithm is shown in fig.1. A monochromatic wavefield extrapolator (MBE) and a memory module are assigned to each frequency in the frequency band to be processed simultaneously. The private memory modules of each MBE are loaded with data from a global data bus. During extrapolation and (partial) imaging, the desired number of frequencies is extrapolated simultaneously by operating the individual MBEs simultaneously. We shall assume that each MBE acts as a systolic ar/ray; data from one depth level is pumped in while data extrapolated to the next depth level is pumped out. Along with data, each MBE must be accompanied by spatial velocity information. Representation for time equal to zero is performed by simultaneous summation of the real part of the result from each MBE.

Selv om denne strukturen er konseptuelt enkel, har den enkelte alvorlige ulemper: Although this structure is conceptually simple, it has some serious drawbacks:

- Addereren må implementers som et "pipelined" adderer tre med lengde lik logaritmen til antall MBE'er. - For et stort antall MBE'er, må all global kommunikasjon implementeres med ("pipelined") tre-strukturer. Dette problem blir alvorligere jo høyere klokkefrekvensen er. - På grunn av tre-strukturene er en eventuell utvidelse av antall MBE'er ikke triviell. - The adder must be implemented as a "pipelined" adder tree with length equal to the logarithm of the number of MBEs. - For a large number of MBEs, all global communication must be implemented with ("pipelined") tree structures. This problem becomes more severe the higher the clock frequency. - Due to the three structures, a possible expansion of the number of MBEs is not trivial.

De ovenfor nevnte problemer kan unngåes ved å benytte følgende design, vist i fig.2. Her er behovet for global kommunukasjon totalt eliminert. Hver MBE er fysisk forbundet kun til sine to naboer. Alle kontrollsignaler, inkluderende systemklokke, fores fra toppen av strukturen. For hver klokkepuls propageres disse signaler en MBE nedover. Linder innlasting av data pumpes data og deres destinasjonsadresser systolisk ned gjennom hukommelsesmodulene. Hver privat hokommelsesmodul aksepterer en skrivkommando kun hvis data adressen ligger innenfor adresseområdet til den aktuelle modul. Under ekstrapolasjon, opereres strukturen som et to-nivå systolisk array. Hastighetsinformas jon og adresseinformasjon pumpes vertikalt nedover. I hver MBE pumpes bølgefeltdata og hastighetsinformas jon horisontalt. Avbilding utføres ved rekursiv summasjon. Partielle resultater pumpes vertikalt nedover og ved hver MBE adderes bidraget fra den aktuelle frekvens. Avbildede (eller delvis avbildede) resultater er tilgjengelige på bunnen av strukturen. Avbilding for andre tidspunkt enn null kan utføres ved å vekte summasjonene slik at man foretar en partiell invers DFT (Discrete Fourier Transform). Samtidig avbilding for forskjellige tidspunkt kan utføres ved å bruke flere " avbildings - kolonner". The above-mentioned problems can be avoided by using the following design, shown in fig.2. Here, the need for global communication is completely eliminated. Each MBE is physically connected only to its two neighbors. All control signals, including the system clock, are routed from the top of the structure. For each clock pulse, these signals are propagated one MBE down. When loading data, data and their destination addresses are systolically pumped down through the memory modules. Each private memory module accepts a write command only if the data address lies within the address range of the relevant module. During extrapolation, the structure is operated as a two-level systolic array. Speed information and address information are pumped vertically downwards. In each MBE, wave field data and velocity information are pumped horizontally. Mapping is performed by recursive summation. Partial results are pumped vertically downwards and at each MBE the contribution from the relevant frequency is added. Imaged (or partially imaged) results are available at the bottom of the structure. Imaging for times other than zero can be performed by weighting the summations so that a partial inverse DFT (Discrete Fourier Transform) is performed. Simultaneous imaging for different time points can be performed by using multiple "imaging columns".

Utvidelse av den foreslåtte struktur til en hvilken som helst størrelse er triviell. Antall MBE'er som kan brukes er begrenset kun av det totale antall frekvenser man vil ta hensyn til. Extension of the proposed structure to any size is trivial. The number of MBEs that can be used is limited only by the total number of frequencies that will be taken into account.

Den struktur som er skissert her har i tillegg den fordel at båndbredden på hastighetsinformasjonen kan begrenses gradvis mens informasjonen propageres nedover. Dette kan oppnås ved å la hver MBE virke som en del av et kaskadefilter. The structure outlined here also has the advantage that the bandwidth of the speed information can be gradually limited while the information is propagated downwards. This can be achieved by having each MBE act as part of a cascade filter.

På grunn av den systoliske vertikale distribusjon av klokkesignaler, vil det være en faseforskyvning proporsjonal med antall MBE'er mellom kontollenheten og bunnen av strukturen. En synkroniseringsenhet kan derfor være nødvendig for å kunne sende resultatene tilbake til kontrollenheten. Due to the systolic vertical distribution of clock signals, there will be a phase shift proportional to the number of MBEs between the control unit and the base of the structure. A synchronization unit may therefore be necessary to be able to send the results back to the control unit.

På grunn av det potensielt høye antall prosessorer må strukturen være i stand til å utføre sin oppgave selv når ikke alle MBE'er er operative. En defekt MBE må derfor fortsatt være i stand til å propagere den informasjon den mottar videre. Hver MBE bør derfor utstyres med et selv-prøvende system som kan aktivere en "bypass" funksjon hvis en feil har oppstått. Due to the potentially high number of processors, the structure must be able to perform its task even when not all MBEs are operational. A defective MBE must therefore still be able to propagate the information it receives further. Each MBE should therefore be equipped with a self-testing system that can activate a "bypass" function if a fault has occurred.

ARKITEKTURER FOR MONOKROMATISK BØLGEFELTEKSTRAPOLASJONARCHITECTURES FOR MONOCHROMATIC WAVEFIELD EXTRAPOLATION

Et systolisk array for rom-variant 1-D (en-dimensjonal) konvolusjon med fast operatorlengde (Kung, 1982) er vist i fig.3. Vekter W er på forhånd lastet inn i arrayet, et sett for hver celle. Både partielle resultater p(x,y,z+Az) og inn-data p(x,y,z) pumpes fra venstre mot høyre, men inn-dataene propageres dobbelt så hurtig som de partielle resultatene. Mer presist; hvert delresultat holdes i cellen det passerer i en periode og bruker derfor dobbelt så lang tid på å passere gjennom arrayet som inn-dataene gjør. Hvert delresultat initialiseres til null før det entrer cellen lengst til venstre og akkumulerer alle sine ledd mens det propagerer mot høyre. I hver celle adresseres vektene av hastighetsinformasjonen bundet til inn-dataene. A systolic array for space-variant 1-D (one-dimensional) convolution with fixed operator length (Kung, 1982) is shown in Fig.3. Weights W are preloaded into the array, one set for each cell. Both partial results p(x,y,z+Az) and input data p(x,y,z) are pumped from left to right, but the input data is propagated twice as fast as the partial results. More precisely; each partial result is held in the cell it passes for a period and therefore takes twice as long to pass through the array as the input data does. Each partial result is initialized to zero before entering the leftmost cell and accumulates all its terms as it propagates to the right. In each cell, the weights are addressed by the speed information bound to the input data.

For symmetriske operatorer kan omtrent halvparten av multiplikasjonene spares ved å addere datapunkter lokalisert symmetrisk om operatormidtpunktet før multiplikasjon. I det tilfellet er det fordelaktig å bruke en arkitektur som vist i fig.4. Denne arkitektur har to uavhengige datastrømmer for inn-data. I den nedre datastrøm pumpes inn-data tre ganger så hurtig som i den øvre og dobbelt så hurtig som resultatstrømmen. Inn-argumentene adderes før de multipliseres med operatorvektene. Operatorvektene adresseres av hastighetsinformasjonen som nå er bundet til resultatstrømmen. For symmetric operators, about half of the multiplications can be saved by adding data points located symmetrically about the operator midpoint before multiplication. In that case, it is advantageous to use an architecture as shown in fig.4. This architecture has two independent data streams for input data. In the lower data stream, input data is pumped three times as fast as in the upper one and twice as fast as the result stream. The input arguments are added before being multiplied by the operator weights. The operator weights are addressed by the velocity information that is now bound to the result stream.

Der er mange mulige måter å gjøre de ovenfor nevnte arkitekturer i stand til å utføre konvolusjon med operatorer som er større i lengde og dimensjon enn arrayene. En enkel løsning er vist i fig.5. I første omgang utføres konvolusjon med en del av den komplette operator. Resultatene mellomlagres i en cache-hukommelse. I neste omgang adresseres en annen del av operatoren og resultatene fra forrige omgang brukes som startverdier. Dette fortsettes til hele operatoren i en dimensjon er dekket. Så kan man starte med en ny rad (av en 2-D operator) og bruke mellomresultatene som startverdier. På denne måten akkumuleres bidraget fra alle deler av operatoren til hele operatoren har dekket dataene. I siste omgang pumpes resultatene ut. There are many possible ways to make the above mentioned architectures capable of performing convolution with operators that are larger in length and dimension than the arrays. A simple solution is shown in fig.5. In the first instance, convolution is performed with a part of the complete operator. The results are temporarily stored in a cache memory. In the next round, another part of the operator is addressed and the results from the previous round are used as starting values. This is continued until the entire operator in one dimension is covered. Then one can start with a new row (of a 2-D operator) and use the intermediate results as starting values. In this way, the contribution from all parts of the operator is accumulated until the entire operator has covered the data. In the last round, the results are pumped out.

Når en 2-D operatormatrise har rang lik en, kan en 2-D konvolusjon utføres som to påfølgende 1-D konvolusjoner i ortogonale retninger. When a 2-D operator matrix has rank equal to one, a 2-D convolution can be performed as two consecutive 1-D convolutions in orthogonal directions.

DISKUSJONDISCUSSION

En systolisk regnemaskinstruktur for hurtig ekstrapolasjon av bølgefelt er foreslått. Strukturen er dekomponert i grunnleggende elementer som brukes repetitivt med enkle grensesnitt. Dette vil betydelig forenkle utvikling, testing og produksjon. Der er er mange muliger varianter av det grunnleggende konsept. I det mest generelle tilfelle kan hver MBE implementeres som en generell regnemaskin og romlig informasjon om en rekke fysiske parametere kan brukes som inn-informas jon. A systolic calculator structure for rapid extrapolation of wavefields is proposed. The structure is decomposed into basic elements that are used repetitively with simple interfaces. This will significantly simplify development, testing and production. There are many possible variations of the basic concept. In the most general case, each MBE can be implemented as a general calculator and spatial information about a number of physical parameters can be used as input information.

Strukturen har opplagte anvendelser også utenfor områder som strengt omfattes av rekursiv bølgefeltekstrapolasjon. Eksempler er romlig datafiltrering og pre stack partiell migrasjon. I den terminologi som er brukt her, kan disse problemer utføres som en (eller noen få påfølgende) ekstrapolasjon. The structure has obvious applications also outside areas that are strictly covered by recursive wavefield extrapolation. Examples are spatial data filtering and pre stack partial migration. In the terminology used here, these problems can be performed as one (or a few consecutive) extrapolations.

Ved å bruke det foreslåtte konsept kan en stor del av den totale potensielle parallellisme utnyttes ved samtidig prosessering av et hvilket som helst antall frekvenser. En alternativ måte å introdusere parallellisme på kunne være å parallellisere i rom. Imidlertid ville en slik fremgangsmåte kreve simultan forsyning av By using the proposed concept, a large part of the total potential parallelism can be exploited by simultaneous processing of any number of frequencies. An alternative way to introduce parallelism could be to parallelize in space. However, such a method would require the simultaneous supply of

hastighetsinformas jon for flere forskjellige romkoordinater. velocity information for several different spatial coordinates.

REFERANSERREFERENCES

Berkhout, A.J., 1982, Seismic migration. Imaging of acoustic energy by wavefield extrapolation. A. Theoretical aspects. Berkhout, A.J., 1982, Seismic migration. Imaging of acoustic energy by wavefield extrapolation. A. Theoretical aspects.

Elsevier Scientific Publ. Co. Amsterdam.Elsevier Scientific Publ. Co. Amsterdam.

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

Claims

1. Multiprocessor calculator structure for wave field extrapolation in the frequency domain or in another domain that allows independent processing of individual spectral components. The structure is characterized by the fact that a memory and a processor, later called a monochromatic wave field extrapolator (MBE), are assigned to each frequency in the frequency band for which parallel processing is desired and that each MBE receives control and control signals from one and the same control unit (figure 1). Each MBE can be a general purpose calculator or a dedicated multiprocessor.

2. Calculator structure according to claim 1, further characterized in that the system clock and other control and control signals are distributed systolically; i.e. propagated from one MBE to the next for each clock pulse. The need for global communication can thus be eliminated (figure 2).

3. Calculator structure according to claim 1, further characterized in that spatial information about the medium the waves propagate in is cascaded filtered while it is pumped down through the structure of MBEs (figure 2). K.

Calculator structure according to claim 1, further characterized in that back-transformation to the time domain takes place by recursive summation in such a way that the contribution from each MBE is propagated one MBE further and added to its contribution (figure 2).

5. Computing machine structure according to claim 1, further characterized in that each MBE is operated as a systolic array where wave field data from one depth is pumped in together with spatial information about the medium the waves propagate in in such a way that the spatial information is used to address it or the operators used to extrapolate the wave field and that wave field data fully or partially extrapolated to the next depth can be pumped out (Figure 3).

6. Calculator structure according to claims 1 and 5, further characterized in that extrapolation with a space-variant symmetric folding operator takes place in the following way (figure *): The unit has two independent paths for input data and one path for (partial) results. The lower stream of input data is propagated three times as fast as the upper stream and twice as fast as the result stream. In each cell, the input arguments are added before multiplication by the operator weight. The operator weights are addressed by spatial information that flows parallel to the result stream.