DE69623975T2

DE69623975T2 - Verfahren zum verarbeiten von seismischen daten in der frequenzdomäne auf einem massiv-parallelen rechner

Info

Publication number: DE69623975T2
Application number: DE69623975T
Authority: DE
Inventors: J. Maher; E. Willen
Original assignee: ExxonMobil Upstream Research Co
Current assignee: ExxonMobil Upstream Research Co
Priority date: 1995-06-16
Filing date: 1996-06-13
Publication date: 2003-05-22
Anticipated expiration: 2016-06-14
Also published as: WO1997000485A1; NO975883D0; JP2001519932A; DE69623975D1; CA2219478C; CA2219478A1; NO975883L; EP0846294B1; EP0846294A4; MY132149A; EP0846294A1; AU707757B2; AU6278796A

Description

GEBIET DER ERFINDUNG

Diese Erfindung betrifft das Gebiet eines geophysikalischen Schürfbetriebs. Insbesondere enthält die Erfindung ein Frequenzdomänenverfahren zum Bearbeiten seismischer Daten unter Verwendung paralleler Prozessoren.

HINTERGRUND DER ERFINDUNG

Die Suche nach unterirdischen Kohlenwasserstoffablagerungen enthält typischerweise eine Sequenz aus einer Datenerfassung, einer Analyse und von Interpretationsprozeduren. Die Datenerfassungsphase enthält eine Verwendung einer Energiequelle zum Erzeugen von, Signalen, die sich in die Erde ausbreiten und sich von verschiedenen unterirdischen geologischen Strukturen reflektieren. Die reflektierten Signale werden durch eine Mehrzahl von Empfängern auf oder nahe der Oberfläche der Erde oder in einem darüber liegenden Körper von Wasser aufgezeichnet. Die empfangenen Signale, die oft seismische Spuren genannt werden, bestehen aus Amplituden akustischer Energie, die sich als Funktion der Zeit ändern, einer Empfängerposition und einer Quellenposition, und die sich, was am wichtigsten ist, als Funktion der physikalischen Eigenschaften der Strukturen ändern, von welchen die Signale reflektieren. Der Datenanalytiker verwendet diese Spuren zusammen mit einem geophysikalischen Modell zum Entwickeln eines Bildes der unterirdischen geologischen Strukturen.
Die Analysephase enthält Prozeduren, die sich in Abhängigkeit von der Art der geologischen Struktur ändern, die untersucht wird, und von den Charakteristiken bzw. Eigenschaften des Datensatzes selbst. Allgemein besteht jedoch der Zweck einer typischen Anstrengung bezüglich einer Bearbeitung seismischer Daten im Erzeugen eines Bildes der geologischen Struktur aus den aufgezeichneten Daten. Das Bild wird unter Verwendung theoretischer und empirischer Modelle von der Art entwickelt, auf welche die Signale in die Erde gesendet werden, durch die unterirdischen Schichten bzw. Formationen gedämpft werden und von den geologischen Strukturen reflektiert werden. Die Qualität des Endprodukts der Datenverarbeitungssequenz ist stark abhängig von der Genauigkeit dieser Analyseprozeduren.
Die Endphase ist die Interpretation der Analysenergebnisse. Insbesondere besteht die Aufgabe eines Interpretierers darin, das Ausmaß festzulegen, bis zu welchem unterirdische Kohlenwasserstoffablagerungen vorhanden sind, um dadurch bei solchen Entscheidungen zu helfen, ob ein zusätzliches Forschungsbohren erwünscht ist, oder was ein Szenario einer optimalen Kohlenwasserstoffwiedergewinnung sein kann. Bei dieser Festlegung enthält die Interpretation des Bildes eine Vielfalt von unterschiedlichen Anstrengungen. Beispielsweise studiert der Interpretierer oft die abgebildeten Ergebnisse, um die regionale unterirdische Geologie zu verstehen. Dies kann ein Markieren von Hauptstrukturmerkmalen enthalten, wie beispielsweise von Fehlern, Synklinalen und Antiklinalen. Danach kann eine vorläufige Konturierung von Zonen durchgeführt werden. Ein darauf folgender Schritt zum kontinuierlichen Verfolgen von Zonen über die verschiedenen vertikalen Schnitte mit Korrelationen der interpretierten Fehler kann auch erfolgen. Wie es im Stand der Technik klar verstanden wird, haben die Qualität und die Genauigkeit der Ergebnisse des Datenanalyseschritts der seismischen Sequenz einen signifikanten Einfluss auf die Genauigkeit und Nützlichkeit der Ergebnisse dieser Interpretationsphase.
Prinzipiell kann das seismische Bild unter Verwendung eines dreidimensionalen geophysikalischen Modells einer seismischen Wellenausbreitung entwickelt werden, um dadurch eine genaue Tiefe und eine azimutale Skalierung aller Reflexionen bei den Daten zu ermöglichen. Genau spezifizierte Reflexionen vereinfachen eine Dateninterpretation stark, da der Interpretationsfokus auf die Art der geologischen Struktur gerichtet werden kann, die beteiligt ist, und nicht auf die Genauigkeit des Bildes. Unglücklicherweise erfordern dreidimensionale geophysikalische Modelle häufig nicht tolerierbar lange Rechenzeiten, und seismische Analytiker werden gezwungen, die Datenverarbeitungsanstrengung weitest möglich zu vereinfachen, um die Belastungen von sowohl einer Analysezeit als auch von Kosten zu reduzieren.
Zusätzlich zu der 3-D-Berechnungsherausforderung steht der Analytiker einer Verarbeitungsvolumenherausforderung gegenüber. Beispielsweise kann eine Ausübung einer typischen Datenerfassung mehrere hundert bis tausende von Quellenstellen enthalten, wobei jede Quellenstelle hunderte von Empfängerstellen hat. Weil jedes Quellen-Empfänger-Paar einen wertvollen Beitrag zum erwünschten Ausgangsbild bilden kann, kann die Datenbehandlungslast (d. h. die Eingabe/Ausgabe-Datentransferanforderung) in sich selbst eine Belastung sein, und zwar unabhängig von der Berechnungsbelastung.
Analytiker von seismischen Daten haben historisch mehrere unterschiedliche Ansätze zum Managen dieser Belastungen verwenden, und zwar direkt oder indirekt. Diese Ansätze beziehen sich hauptsächlich auf entweder die Art, auf welche die Ausübung einer Datenerfassung entwickelt ist und ausgeführt wird, oder auf die während der Datenanalyseanstrengung gemachten Annahmen. In beiden Fällen kann die Qualität der Ausgabe der Dateninterpretationsprozedur direkt beeinflusst werden. Diese Ansätze werden am einfachsten in Zusammenhang mit Fig. 1 diskutiert, die eine perspektivische Ansicht eines Bereichs 20 der Erde zeigt, für welchen ein geophysikalisches Bild erwünscht ist. Auf der Oberfläche 18 der Erde ist eine Anzahl von Aufnahmelinien 2 gezeigt, entlang welchen die seismischen Daten erfasst werden. Wie es in Fig. 1A gezeigt ist, bestehen die Aufnahmelinien 2 aus einer Sequenz von Positionen, bei welchen eine seismische Quelle 3 angeordnet ist und von welchen seismische Signale 5 in die Erde gesendet werden. Entlang jeder Linie angeordnete Empfänger 4 empfangen die Signale von jeder Quellenposition nach einer Reflexion von verschiedenen unterirdischen Reflektoren 6.
Ein erstes Verfahren zum Managen der oben diskutierten Belastungen bezüglich seismischer Daten enthält eine sorgfältige Definition des Bereichs, über welchen die Daten erfasst werden. Insbesondere kann eine Verwendung irgendeiner verfügbaren vorläufigen geologischen und geophysikalischen Information die Minimierung des Oberflächengebiets erleichtern, über welchem es nötig sein kann, dass seismische Daten erfasst werden. Eine solche Minimierung wird die Menge an Daten direkt reduzieren, die letztlich erfasst werden. Weiterhin wird ein gleichermaßen sorgfältiges Planen der Beabstandung zwischen Aufnahmelinien, die Analyseanstrengung durch Reduzieren eines Datenvolumens optimieren. Schließlich wird auch eine Optimierung der Anzahl von Quellen und Empfängern, die verwendet werden, und der Beabstandung zwischen benachbarten Quellen- und Empfängerpositionen auch einen Vorteil für den Datenanalytiker haben.
Keine dieser Anstrengungen kann ohne einen Gegenaufwand bzw. Kompromiss erreicht werden. Beispielsweise reduziert eine relativ breite Beabstandung zwischen Aufnahmelinien oder zwischen Quellen und Empfängern die Auflösung des berechneten seismischen Bildes, was eine Interpretation schwieriger macht. Zusätzlich kann es sein, dass komplexe geologische Eigenschaften bzw. Merkmale ohne eine relativ nahe Beabstandung nicht auflösbar sind. Schließlich lassen bestimmte Ausübungen einer Datenerfassung, wie beispielsweise in relativ unerforschten Gebieten, keine Optimierung des Oberflächengebiets zu, über welchem Daten erfasst werden sollen. Als Ergebnis kann die Belastung einer Datenhandhabung nicht vollständig durch eine Planung der Datenerfassung eliminiert werden.
Verfahren zum Minimieren der Berechnungsbelastung sind oft während einer Datenanalyse implementiert. Eine gemeinhin verwendete Technik enthält eine Verwendung eines zweidimensionalen geophysikalischen Modells. Beispielsweise sind in Fig. 1A die Signale für jede Quelle als in der Ebene wandernd gezeigt, die direkt unterhalb der Aufnahmelinie ist, auf welcher die Quelle liegt. Somit wird angenommen, dass sich das Signal unabhängig von geologischen Strukturen außerhalb der Ebene ausbreitet. Diese vereinfachende Annahme lässt eine Verwendung von zweidimensionalen geophysikalischen Modellen beim Bilderzeugungsprozess zu, und, wie es wohlbekannt ist, können zweidimensionale Analyseprozeduren sehr viel recheneffizienter als dreidimensionale Analyseprozeduren sein.
Es existieren Beschränkungen bezüglich der 2-D-Analyseannahme. Geologische Strukturen sind selten, wenn überhaupt, zweidimensional; diese Annahme kann daher zu ungenau spezifizierten Bildern führen. Weil allgemein wenig über die geologische Struktur bekannt ist, die untersucht wird, kennt der Analytiker normalerweise nicht das Ausmaß, bis zu welchem das Bild fehlerhaft ist. Zusätzlich muss der Interpretierer deshalb, weil jede Ebene unabhängig analysiert wird, die Bilder für jede Ebene an ein jeweiliges der anderen durch eine Interpolation oder andere ähnliche Interpretationsverfahren anhängen, wenn ein kontinuierliches Bild über den gesamten kubischen Bereich erwünscht ist. Schließlich können einige komplexe Strukturen, wie beispielsweise fehlerhafte Bereiche und Salzeigenschaften, lediglich durch eine Verwendung von zweidimensionalen Verfahren nicht genau analysiert werden.
Aufgrund dieser und anderer Beschränkungen, die Analytiker für seismische Daten lange eingeschränkt haben, ist die Mineralölindustrie typischerweise ein früher Anwender einer neu entwickelten Hochgeschwindigkeits-Computerhardware gewesen. Da jede neue Generation einer Anlage verfügbar geworden ist, sind Analyseroutinen, die vollständige dreidimensionale Analysefähigkeiten implementieren, allgemeiner verwendet worden. Nichtsdestoweniger ist es nicht ungewöhnlich, dass signifikante Computerzeiten an komplexen Analysen beteiligt sind, was oft Wochen oder Monate einer tatsächlichen Verarbeitungszeit enthält.
Die letzte Verfügbarkeit von massiv parallelen Prozessoren bietet eine signifikante Gelegenheit für Analytiker von seismischen Daten. Massiv parallele Prozessoren (MPPs) haben mehrere Zentralverarbeitungseinheiten (CPUs), die gleichzeitige Berechnungen durchführen können. Durch effiziente Verwendung bzw. Ausnutzung dieser CPUs können die für komplexe Analysen bislang erforderlichen Wochen oder Monate auf einige Tage reduziert werden, oder vielleicht ein paar Stunden. Jedoch kann dieser signifikante Vorteil nur realisiert werden, wenn effiziente Berechnungsalgorithmen in der MPP-Software codiert sind. Somit schafft die Gelegenheit, die die MPPs Analytikern von seismischen Daten bieten, auch eine Herausforderung für die Entwicklung von geeigneten Berechnungsalgorithmen, die einen Vorteil aus den Mehrfach- CPUs ziehen.
Diese Herausforderung kann einfach dadurch diskutiert werden, dass man die Art betrachtet, auf welche Berechnungsalgorithmen für existierende seismische Analyseroutinen am weitesten verbreitet geschrieben worden sind. Bis zuletzt beruhten Computer auf einem Operationsmode, der sequentielle Berechnung genannt wird. Eine sequentielle Berechnung enthält eine Verwendung von analytischen Routinen, die nur eine einzige Prozedur durchführen, oder sich vielleicht auf eine einzelne Untergruppe der Daten oder eines Bildes richten, und Zwar zu einer gegebenen Zeit. Dies ist ein direktes Ergebnis eines Computers mit nur einer CPU. Aus diesem Grund sind die einzigen Optimierungsprozeduren, die auf Computern mit einer einzelnen CPU verwendet werden können, diejenigen, die die Effizienz der Verarbeitung bezüglich der Prozedur oder der Untergruppe erhöhen. Weil jedoch alle Berechnungen letztlich durch diese einzelne CPU durchgeführt werden müssen, sind die Optionen zum Erhalten einer hohen Leistungsfähigkeit von Natur aus begrenzt.
Andererseits bietet die Mehrfach-CPU-Kapazität von MPPs einen offensichtlichen Vorteil einer gleichzeitigen Berechnung. Dieser Vorteil besteht darin, dass die zum Lösen eines Rechenproblems erforderlich Gesamtzeit durch Unterteilen der durchzuführenden Arbeit auf die verschiedenen CPUs reduziert werden kann, vorausgesetzt, dass die Unterteilung zulässt, dass eine jeweilige CPU eine nützliche Arbeit durchführt, während die anderen CPUs auch eine Arbeit durchführen. Unglücklicherweise besteht der Nachteil von Mehrfach-CPU- Hardware darin, dass die sequentiellen Verarbeitungsverfahren, die in der Softwareentwicklung lange Zeit verwendet worden sind, durch geeignetere parallelisierte Berechnungsverfahren ersetzt werden müssen. Einfach gesagt erfordern MPPs, dass Verarbeitungsverfahren entwickelt werden, die die Mehrfach-CPU-Hardware effizient verwenden. Idealerweise sollen diese Verfahren die Verteilung von Arbeit relativ gleichmäßig auf die Prozessoren organisieren, und sicherstellen, dass alle Prozessoren nötige Berechnungen eher alle gleichzeitig durchführen, als auf Zwischenergebnisse von anderen Prozessoren zu warten.
Die Herausforderung, parallelisierte Verarbeitungsverfahren zu definieren und diese einmal definierten parallelisierte Verfahren zu optimieren, ist auf dem Gebiet der Verarbeitung bzw. Bearbeitung von seismischen Daten besonders akut. Seismische Daten bestehen allgemein aus einer großen Anzahl von individuellen Spuren, die, jeweils etwas unabhängig von den anderen Spuren aufgezeichnet sind. Es ist somit logisch, dass sequentielle Berechnungsverfahren, die erfordern, dass der Analysefokus auf eine einzige Berechnung zu einer Zeit gerichtet, wird, gut an einer Analyse dieser unabhängigen Spuren angepasst wird. Dies gilt sogar dann, wenn Berechnungsengpässe existieren können. Beispielsweise können Teile der Analysesequenz relativ mehr Rechenzeit als andere Teile erfordern, müssen fertig gestellt sein, bevor andere Berechnungen fortfahren können, oder können auf gleichen Eingangsdaten wie andere Spuren beruhen wie beispielsweise Wanderzeiten. Da keine gleichzeitigen Berechnungen bei einer sequentiellen Verarbeitung auftreten, führt keiner dieser Engpässe zu einer Reduzierung bezüglich einer Berechnungseffizienz bei einer einzelnen CPU, außer bezüglich der gesamten Verarbeitungszeit, die erforderlich ist. Mit Ausnahme dieser Anforderung bezüglich der Gesamtzeit erlegt die Existenz solcher Berechnungsengpässe sonst keine Probleme für den Analytiker auf. Zum Ziehen des vollen Vorteils aus den Berechnungsfähigkeiten von MPP jedoch, wo das Ziel im Durchführen einer gleichzeitigen Verarbeitung in allen CPUs besteht, müssen Verfahren zum Optimieren der seismischen Analysephase durch Eliminieren von solchen Engpässen entwickelt werden.
Dieser Vorteil eines MPP wird dadurch deutlich, dass man die Begrenzung betrachtet, die eine Rechenzeit einer Bildbereichsgröße in Einzel-CPU-Computern auferlegt. Ein Erhöhen der Größe des Bildes, z. B. durch Ausdehnen der Größe eines Kubus 20 in Fig. 1, oder ein Erhöhen der zu verarbeitenden Menge an Daten, z. B. durch Hinzufügen zusätzlicher Quellen 3 und Empfänger 4 zu Aufnahmelinien 2, erhöht die Gesamtberechnung. Dieser direkte Einfluss auf eine Berechnungszeit legt einem seismischen Analytiker eine schwere Belastung auf, um eine Bildgröße zu optimieren, und zwar insbesondere deshalb, weil auch kleine Bildbereiche Wochen einer Berechnungszeit auf selbst den Computern mit einer sequentiellen Verarbeitung mit höchster Geschwindigkeit erfordern können. Gegensätzlich dazu sollte eine effiziente Verarbeitung auf MPPs, die so viele wie 256 einzelne CPUs oder noch mehr haben können, nur minimal lange Berechnungszeiten enthalten, da jede CPU nur einen Bruchteil, wie beispielsweise 1/256 der zusätzlichen Arbeit annehmen würde, die durch den größeren Berech erforderlich ist. Dieses Potential zur Skalierbarkeit des Bildbereichs und die Arbeitslast, die bei einer Bilderzeugung erforderlich ist, ist ein Hauptvorteil von MPPs, was ein Vorteil ist, der nur dann realisiert werden kann, wenn parallelisierte seismische Bearbeitungsverfahren, die eine solche Arbeitslastskalierung zulassen, entwickelt werden.
Grundsätzliche Betrachtungen zum Bestimmen effizienter parallelisierter seismischer Bearbeitungsverfahren werden durch erneutes Betrachten des obigen Überblick über den seismischen Analyseprozess offensichtlich. Wie es angemerkt ist, besteht der Zweck einer seismischen Analyse im Analysieren gemessener seismischer Daten unter Verwendung geophysikalischer Modelle zum Entwickeln von Bildern des Unterirdischen. Daher kann jede der drei Hauptverarbeitungskomponenten -- Daten, Modell und Bild -- derart betrachtet werden, dass sie ein Kandidat zum Verteilen von Berechnungsarbeit auf die verschiedenen Prozessoren in einem MPP ist. Eine Option zum Verteilen von Arbeit auf die Prozessoren wäre ein Zuordnen unterschiedlicher Gruppen der eingegebenen seismischen Spurdaten bzw. Verfolgungsdaten zu unterschiedlichen Prozessoren. Beispielsweise können Spuren durch Quellenstellen gruppiert sein, wobei unterschiedliche Prozessoren unterschiedlichen Gruppen zugeordnet sind. Gleichermaßen könnte das ausgegebene Bild unterteilt und unterschiedlichen Prozessoren zugeordnet werden. Schließlich kann es auch möglich sein, dass zum Erzeugen des ausgegebenen Bildes verwendete geophysikalische Modelt in Gruppierungen zu unterteilen, die den Verschiedenen Prozessoren zugeordnet werden können. (Dieses Modell wird allgemein derart angesehen, dass es in den arithmetischen Operationen verkörpert ist, die durch das mathematische Modell erforderlich sind, das der Gegenstand der Verarbeitungsanstrengung ist. Beispielsweise basiert das mathematische Modell bei einer seismischen Analyse oft auf der Wellengleichung.) Beispielsweise können die Daten in den Frequenzbereich bzw. die Frequenzdomäne transformiert werden, und zwar mit einzelnen Frequenzen einzelnen Prozessoren zugeordnet. Es kann auch möglich sein, Kombinationen dieser Ansätze zu entwickeln. Beispielsweise können Gruppen von Prozessoren einer kollektiven Verantwortlichkeit für spezifische Frequenzen im Modell und alle Tiefen im Bild zugeordnet werden, während sie eine individuelle Verantwortlichkeit für spezifische horizontale Stellen im Bild haben. Die Herausforderung für den Analytiker von seismischen Daten besteht im Bestimmen von Verfahren zum Unterteilen der seismischen Daten, des Modells und des Bilds in Komponenten, die einzelnen Prozessoren im MPP zugeordnet werden können, um dadurch zuzulassen, dass Berechnungen in jedem Prozessor unabhängig von anderen Prozessoren durchgeführt werden. Diese Unterteilung einer Analyse von seismischen Daten in individuelle Komponenten wird gemeinhin seismische Zerlegung genannt.
Ein Typ von MPP hat von tausenden bis zu zehntausenden relativ unausgeklügelten Verarbeitungselementen. Bei dieser Art von Maschine führen die Verarbeitungselemente typischerweise dieselbe Operation auf mehreren Datenströmen durch, und sie ist eine Einzelbefehls-Mehrfachdaten-Strom- (SIMD-)Maschine. Ein Beispiel ist CM2, ein Produkt der Thinking Machines Corporation. Diese Arten von Maschinen hatten typischerweise keinen gemeinsam genutzten Speicher, d. h. jeder Prozessor hat seine eigene separate Speichereinheit, und auf die Information im Speicher kann von anderen Prozessoren nicht direkt zugegriffen werden. Die einzelnen Prozessoren haben typischerweise eine begrenzte Rechenkapazität und einen begrenzten Speicher. Aufgrund der großen Anzahl von Verarbeitungselementen und einem Mangel an gemeinsam genutzten. Speicher ist ein Datentransfer zwischen, den Verarbeitungselementen ein größerer Engpass bei einer effizienten Ausnutzung der Kapazität bzw. Fähigkeit der Maschinen. Selbst mit ausgeklügelten Verbindungstechniken, wie beispielsweise einem Hyperkubusaufbau, ist ein Transfer von Daten zwischen Prozessoren ein Hauptfaktor bei der Laufzeit von Programmen.
Andere Computer haben viel leistungsfähigere Elemente in Feldern bzw. Matrizen von zehnfachen oder hundertfachen. Der T3D, ein Produkt von Cray Research Corporation, ist ein Beispiel für diese Art von Maschine. Daneben, dass er einzelne Verarbeitungselemente hat, die sehr viel mächtiger als diejenigen im CM2 sind, hat der T3D weniger der Elemente und einen physikalisch verteilten, logisch gemeinsam genutzten Speicher. Diese Mehrfachbefehls-Mehrfachdatenstrom- (MIMD-)Maschine hat verschiedene Elemente, die verschiedene Operationen an unterschiedlichen Teilen der Daten gleichzeitig durchführen. Die reduzierte Anzahl von Verarbeitungselementen bedeutet, dass Daten nicht zu so vielen Elementen wie in einer SIMD-Maschine transferiert werden müssen. Aufgrund der erhöhten Ausgeklügeltheit und der Kosten der einzelnen Elemente und aufgrund ihrer geringeren Anzahl erfordert eine effiziente Ausnutzung, dass die Last an den Verarbeitungselementen ausgeglichen wird. Ein zusätzlicher Faktor besteht darin, dass jedes Verarbeitungselement eine größere Untergruppe des gesamten Datenvolumens unterbringen muss; Berechnungen, die ein Sortieren der Daten enthalten, könnten komplizierter werden.
Das US-Patent Nr. 5,404,296, das Moorhead erteilt ist, offenbart und beansprucht ein Verfahren zur Wanderung auf einem MPP, wobei es eine große Anzahl von Verarbeitungselementen gibt, die in einem vorausgewählten, regelmäßigen Muster angeordnet sind. Die. Daten werden anfangs aufnahmesortiert und in Blöcke von Aufnahmen aufgeteilt, so dass das Produkt der Anzahl von Empfängerpositionen in der Datengruppe und der Anzahl von Aufnahme in einem Block gleich der Gesamtanzahl von verfügbaren Verarbeitungselementen ist. Weil seismische 3-D-Datenvolumen typischerweise hunderte von Aufnahme- und Empfängerpositionen enthalten, sind die Offenbarung und die Ansprüche auf SIMD-Maschinen beschränkt.
Ein Ansatz, der die Menge der 3-D-Verarbeitung reduziert, besteht im Abbilden einer teilverarbeiteten Datengruppe. Dies ist insbesondere dann nützlich, wenn die Teilverarbeitungsoperation berechnungsmäßig einfacher als ein vollständiger Abbildungsprozess ist. Beispielsweise könnten die Daten vor einer Wanderung unter Verwendung einer herkömmlichen Geschwindigkeitsanalyse für ein normales Bewegen nach außen (NMO = Normal Moveout) gestapelt werden. Wie es Fachleuten auf dem Gebiet vertraut sein wird, sind Algorithmen, die eine Migration von Null-Offsetdaten durchführen, berechnungsmäßig weniger intensiv als diejenigen, die eine Migration von nicht gestapelten Daten durchführt. Die Stapeloperation verbessert allgemein das Signal-zu-Rausch-Verhältnis der Daten im Vergleich mit demjenigen der einzelnen Spuren. Zu dem Ausmaß, dass die gestapelten Daten genau eine hypothetische Null-Offset-Spur darstellen, wird die Berechnungsbelastung reduziert.
Im US-Patent Nr. 5,349,527, offenbaren Pieprzak und Highnam ein Verfahren zur Migration von einer solchen Datengruppe auf einem SIMD-MPP. Das Verfahren beginnt mit einem 3-D- Datenvolumen, das aus einer zeitlichen Serie (vertikale Achse) von seismischen Reflexionsdaten auf einem regelmäßigen rechteckigen Gitter (horizontalen Achsen) auf der Oberfläche der Erde [(t, x, y) Datenbereich bzw. Datendomäne] besteht. Eine Fouriertransformation wird an den seismischen Reflexionsdaten durchgeführt, um ein 3-D-Volumen zu ergeben, das aus Amplituden als Funktion einer Frequenz (vertikale Achse) auf dem rechteckförmigen Oberflächengitter [(ω, x, y) Datendomäne, wobei ω die temporale Frequenz ist] besteht. Die Daten werden in Frequenz-Unterbänder aufgeteilt, und innerhalb jedes Frequenz-Unterbandes werden Daten bei ausgewählten Frequenzen in Untergruppen des rechteckförmigen Oberflächengitters aufgeteilt und individuellen bzw. einzelnen Prozessoren zugeordnet. Dies ist aufgrund der begrenzten Speicherkapazität der einzelnen Prozessoren nötig. Die Beschränkungen der einzelnen Prozessoren benötigt auch eine Hyperkubusanordnung der Prozessoren, um die Zeit zu reduzieren, die bei einem Datentransfer zwischen den Prozessoren beteiligt ist.
Bei der Erfindung von Pieprzak und Highnam wird die darauf folgende Migration durch einen iterativen zweistufigen Prozess erreicht. Der erste Schritt ist die Fortführung einer einzelnen Frequenzkomponente über ein Tiefenintervall in Richtung nach unten bzw. die Abwärtsfortsetzung. Die Beiträge aller Frequenzen werden dann summiert, um ein migriertes Bild am unteren Ende des Tiefenintervalls zu erzeugen. Die zwei Schritte werden für aufeinander folgende Tiefenintervalle wiederholt, bis die maximale Tiefe erreicht ist.
Die Fortführung in Richtung nach unten wird unter Verwendung des McClellen Transformationsverfahrens durchgeführt, das in Hale, D., 3-D Depth Migration via McClellan Transformations, 56 Geophysics 1778 (1991), angegeben. Die aktuelle Implementierung wird mittels eines Faltungsfilters durchgeführt. Wie es Fachleuten auf dem Gebiet bekannt sein wird, hat das Hale-Verfahren eine Schwierigkeit beim Behandeln von steilen Gefällen, und die Rechenbelastung wird bei Gefällen sehr groß, die größer als 60º sind. Zum Reduzieren der Rechenbelastung lehren. Pieprzak und Highnam die Verwendung eines rekursiven Tschebyscheff-Filters für den Faltungsschritt. Die Koeffizienten des Filters werden aus dem Geschwindigkeitsmodell des Unterirdischen im Voraus berechnet und redundant im lokalen Speicher der Verarbeitungselement gespeichert. Ein solcher Ansatz ist bei SIMD-Computern geeignet, bei welchen die Verarbeitungselemente nur ein kleines Untergitter von (x, y) Punkten verfolgen müssen und nicht die Fähigkeit zum Berechnen der Filterkoeffizienten "im Vorübergehen" haben.
Fachleute auf dem Gebiet erkennen, dass für große Probleme Frequenzdomänenverfahren signifikante Berechnungskosteneinsparungen bieten können. Weiterhin sind die Migrationsverfahren für ein zweifaches Durchlassen nach dem Stand der Technik rechenmäßig schnell, weil sie die Migration in der x- und y-Richtung separat durchführen. Pieprzak und Highnam verwenden einen Einzeldurchlauf- Frequenzraumdomänenmigrationsprozess, der berechnungsmäßig genauer als der im Stand der Technik verwendete Zweifachdurchlass-Migrationsprozess ist. Jedoch erfolgt die erhöhte Genauigkeit des Einzeldurchlassverfahrens auf die Kosten einer sehr erhöhten Berechnungslast. Eine Option zum Überwinden dieser erhöhten Berechnungslast besteht im Verwenden einer Frequenzwellenzahlmigration, die signifikant schneller als eine Frequenzraummigration sein kann. Jedoch wäre ein auf der Offenbarung von Pieprzak und Highnam basierendes Frequenzwellenzahlschema aus dem Grund nicht praktikabel, dass die Transformation vom räumlichen Bereich bzw. der räumlichen Domäne zur Wellenzahldomäne teuer wird, wenn jeder Prozessor einen Zugriff auf nur einen kleinen Teil des x-y-Gitters hat.
Die 1988 Berryhill, Gonzalez und Kim erteilte US-gesetzmäßige Erfindungsregistrierung H 482 offenbart ein Verfahren zur rekursiven Zeitmigration im Frequenzwellenzahlbereich. Dieses Verfahren nimmt an, dass das Unterirdische der Erde durch Schichten modelliert werden kann, in welchen die Geschwindigkeit konstant ist. Jedoch ist das Verfahren ineffizient zum Implementieren auf einem MPP, weil es aufeinander folgende Fouriertransformationen in unterschiedlichen Richtungen über ein Datenvolumen erfordert, das gänzlich innerhalb eines Prozessorspeichers ist, und die geophysikalischen Datengruppen von Interesse für die Industrie nicht vollständig innerhalb eines Prozessorspeichers gespeichert werden können.
Es wäre daher wünschenswert, ein Verfahren zu haben, das eine Einzeldurchlassmigration im Frequenzwellenzahlbereich auf einem MPP auf eine kosteneffektive Weise durchführen kann. Die vorliegende Erfindung erfüllt diese Notwendigkeit.

ZUSAMMENFASSUNG DER ERFINDUNG

Die vorliegende Erfindung, wie sie beansprucht ist, ist ein Verfahren zum Bearbeiten von seismischen Daten auf parallelen Prozessoren, und zwar vorzugsweise auf einem massiv parallelen Prozessor (MPP). Das Endprodukt der Erfindung ist ein migrierter seismischer Schnitt in der Form eines Bildes des Unterirdischen der Erde. Die Eingabe zum Prozess können entweder gestapelt seismische Daten oder durch ein Herausbewegen aus einer Steigung (DMO = Dip Moveout) korrigierte allgemeinen bzw. gemeinsame seismische Offset- Daten sein. Der Prozess migriert seismische Daten in der Frequenzwellenzahldomäne, um dadurch die Rechengeschwindigkeit über Frequenzraumverfahren signifikant zu erhöhen.
Ein Geschwindigkeitsmodell wird durch die einzelnen Prozessoren des MPP verwendet, um Parameter für die Migrationsoperation zu erhalten. Die einzelnen Prozessoren führen Migrationsberechnungen unabhängig von den anderen Prozessoren durch. Die Erfindung teilt die Daten auf die einzelnen Prozessoren auf, um die darauf folgende Verarbeitung effizient zu machen. Der Prozess verwendet eine Reihe von Matrizentranspositionen, um die Fouriertransformationen effizient zu implementieren, die für die Migration erforderlich sind, und zwar auf eine effiziente Weise auf einem MPP.
Bei einem Ausführungsbeispiel werden die Daten migriert, um ein Bild bezüglich der Tiefe von einem 3-D-Datenvolumen zu erzeugen. Dieses Ausführungsbeispiel verwendet ein Berechnungsschema bei welchem die Geschwindigkeit einer Ausbreitung von seismischen Wellen im Unterirdischen sich lateral und vertikal ändern kann. Die Daten werden in Richtung nach unten in einem Frequenzraumbereich von einer Tiefe zur nächsten unter Verwendung des Geschwindigkeitsmodells fortgeführt; eine Summierung der verschiedenen Frequenzkomponenten erzeugt ein migriertes Tiefenbild. Die Fortführung in Richtung nach unten bzw. die Abwärtsfortsetzung wird in zwei Schritten durchgeführt: Zuerst wird eine paaraxiale Annäherung an die Wellengleichung unter Verwendung eines impliziten finiten Differenzverfahrens mit einem Aufteilen der Gleichung in der gleichgerichteten x- und der quer gerichteten y-Richtung gelöst. Dies ist analog zum Stand der Technik, der einen Zweifachdurchlass- Migrationsprozess verwendete. Das Aufteilen reduziert die Rechenbelastung um mehrere Größenordnungen, führt aber einen Fehler in der Lösung ein. Als Zweites wird eine Fehlerkorrektur durchgeführt, die den durch den Prozess eines Aufteilens der Wellengleichung in x- und y-Komponenten eingeführten Fehler kompensiert. Die Fehlerkorrektur ist berechnungsmäßig schnell, weil sie im Frequenzwellenzahlbereich durchgeführt wird, und sie nicht bei jedem Tiefenintervall ausgeführt werden muss. Das Ergebnis ist eine genaue Lösung der Wellengleichung, die für steile Gefälle gültig ist.
Die Lösung der Wellengleichung für den Schritt für ein nach unten gerichtetes Fortführen erfordert die Lösung eines großen tri-diagonalen Systems von Gleichungen. Die vorliegende Erfindung enthält ein neues "Brennen-an-beiden- Enden"-Verfahren für eine Lösung eines solchen Systems von Gleichungen auf einem MPP. Durch effiziente Verwendung der Prozessoren führt das Verfahren die Fortführung in Richtung nach unten in etwa der Hälfte der Berechnungszeit durch, die bei Verfahren nach dem Stand der Technik zum Lösen eines solchen Gleichungssystems nötig wäre.
Bei einem zweiten Ausführungsbeispiel ist das Endprodukt des Prozesses ein migriertes Zeitpunkt. Die seismischen Daten (gestapelte Spuren oder DMO-korrigierte Spuren) für ein 3-D- Volumen werden nach einer Vorverarbeitung auf einen MPP- Computer in Richtung nach unten fortgeführt. Das Verfahren beginnt an der Oberfläche und geht Zeitstufe für Zeitstufe unter Verwendung eines Phasenverschiebungsverfahrens weiter. Eine Summierung der Beiträge der einzelnen Frequenzen erzeugt ein migriertes Bild. Das zweite Ausführungsbeispiel ist im Frequenzwellenzahlbereich implementiert, wodurch eine signifikante Erhöhung bzw. der Rechengeschwindigkeit gegenüber dem Stand der Technik geschaffen wird.
Bei einem dritten Ausführungsbeispiel ist das Endprodukt des Prozesses auch ein migrierter Zeitschnitt. Das Geschwindigkeitsmodell ist eine stückweise kontinuierliche Funktion der Tiefe. Beginnend an der Oberfläche mit Daten in einem (t, x, y)-Kubus wird die so genannte Stolt-Migration über jedes Intervall angewendet. Dies wird rekursiv durchgeführt, und zwar unter einem Wechsel von der Stolt- Migration und einer Fortführung in Richtung nach unten des Wellenfeldes in jeder Schicht. Dieses Ausführungsbeispiel ist auch im Frequenzwellenzahlbereich implementiert. Für unterirdische Modelle, die relativ dicke Schichten konstanter Geschwindigkeit haben, bietet dieses Ausführungsbeispiel eine signifikante Verbesserung gegenüber dem zweiten Ausführungsbeispiel.

KURZE BESCHREIBUNG DER ZEICHNUNGEN

Die vorliegende Erfindung und ihre Vorteile werden durch Bezugnahme auf die folgende detaillierte Beschreibung und die beigefügten Zeichnungen und Tabellen einfacher verstanden werden, wobei:
Fig. 1 und 1A die bei einer seismischen 2-D-Erforschung typischerweise verwendete Datenerfassungskonfiguration schematisch darstellen.
Fig. 2 die Vorverarbeitung der für ein Tiefenmigrations-Ausführungsbeispiel erforderlichen Daten zeigt.
Fig. 3A und 3B den gesamten Datenfluss bei der Erfindung und bei einem ersten Tiefenmigrations- Ausführungsbeispiel beteiligte konzeptmäßige Schritte zeigen.
Fig. 4 die bei der Fortführung in Richtung nach unten und Abbildungsberechnungen für das in den Fig. 3A und 3B gezeigte Ausführungsbeispiel beteiligten Schritte darstellt.
Fig. 5A, 5B und 5C das Rundsenden von Geschwindigkeitsscheiben zwischen Prozessorgruppen, die Aufteilung einer Frequenzscheibe auf die PEs in einer zu einer anderen Prozessorgruppe und die Fortführung in Richtung nach unten einer einzelnen Frequenzscheibe von der Oberfläche, die um ein Tiefenintervall darunter ist, piktoriell zeigen.
Fig. 6A und 6B die Prozedur zum Wiedergewinnen zuvor berechneter Teilbilder und das Verfahren zum Kombinieren von Teilbildern in einem Tiefenmigrations-Logarithmus-Fan-in zeigen.
Fig. 7A, 7B, 7C, 7D und 7E den gesamten Datenfluss bei Zeitmigrations-Ausführungsbeispielen der vorliegenden Erfindung zeigen.
Fig. 8 die Berechnungszeitverbesserungen für ein Tiefenmigrations-Ausführungsbeispiel dieser Erfindung im Vergleich mit einer Technik nach dem Stand der Technik zeigt.

DETAILLIERTE BESCHREIBUNG DER ERFINDUNG

Die vorliegende Erfindung ist ein Verfahren zum Bearbeiten von seismischen Daten auf einem parallelen Prozessor. Das Verfahren kann auf irgendeinem Typ eines parallelen Verarbeitungscomputers verwendet werden, der mehr als einen Prozessor zum Durchführen einer Analyse und/oder für Eingabe/Ausgabe-Aufgaben verfügbar hat. Beispielsweise kann ein Cray-Y-MP, der bis zu acht Prozessoren hat, verwendet werden. Vernetzte Personalcomputer oder Workstations können auch verwendet werden. Bei einem bevorzugten Ausführungsbeispiel wird das Verfahren auf einem massiv parallelen Prozessor (MPP) mit 64 oder mehr Verarbeitungselementen verwendet werden. Ein Beispiel für einen solchen MPP ist T3D, der von Cray Research Corporation hergestellt ist. Der Annehmlichkeit halber und nicht als beschränkend gedacht, wird hierin die Abkürzung MPP zur Bezugnahme auf irgendeinem parallelen Verarbeitungscomputer verwendet werden, der für das gegenwärtige Verfahren geeignet ist. Jedoch ist die Verwendung von MPP auf diese Weise keineswegs als beschränkende Anwendung des vorliegenden Verfahrens auf nur diejenigen Computer gedacht, die durch den kommerziellen Marktplatz als massiv parallel bestimmt sind.
Während die Implementierung in Bezug auf eine MIMD-Maschine diskutiert wird, ist die Erfindung nicht auf sie beschränkt. Abweichend von Geschwindigkeitsdaten enthält das vorliegende Verfahren eine Verwendung einer relativ kleinen Anzahl von Eingangsverarbeitungsparametern und kann daher auch nützlich auf SIMD MPPs implementiert werden. Bei sowohl MIMD- als auch SIMD-Implementierungen werden die Eingangsparameter in einem Prozessor gelesen und zu den anderen Prozessoren rundgesendet.
Wie es wohlbekannt ist, beginnen seismische Datenverarbeitungsverfahren allgemein mit Daten von einem Teil der Erdoberfläche. Fig. 1 zeigt eine perspektivische Ansicht eines Bereichs 20 der Erde, für welchen ein geophysikalisches Bild erwünscht ist. Auf einer Oberfläche 18 der Erde ist eine Anzahl von Aufnahmelinien 2 gezeigt, entlang welchen die seismischen Daten erfasst werden. Wie es in Fig. 1A gezeigt ist, bestehen die Aufnahmelinien 2 aus einer Sequenz von Positionen, bei welchen eine seismische Quelle 3 angeordnet ist und von welchem seismische Signale 5 in die Erde gesendet werden. Entlang jeder Linie angeordnete Empfänger 4 empfangen die Signale von jeder Quellenposition nach einer Reflexion von verschiedenen unterirdischen Reflektoren 6. Wie es Fachleuten auf dem Gebiet bekannt sein wird, ist es nun weit verbreitet, dass die Erfassung mit mehreren Zeilen bzw. Reihen von Quellen und Empfängern ausgeführt wird. Die vorliegende Offenbarung hat die Bedeutung, dass diese Situationen auch umfasst sind und die vorläufigen Verarbeitungsschritte, die am "Binärmachen" einer solchen Datengruppe beteiligt sind, sind Teil des Standes der Technik.
Diese Reflexionssignale werden typischerweise vorläufig unter Verwendung von einer von zwei bekannten Techniken verarbeitet. Zuerst kann sich die vorläufige Verarbeitung auf ein Erzeugen eines gestapelten Datenvolumens richten. Der Prozess zum Stapeln besteht aus einem Nehmen von Daten von einer einzelnen Linie 2. Durch einen Prozess einer Geschwindigkeitsanalyse und einer Korrektur für ein normales Herausbewegen (NMO) wird ein ungefährer seismischer Null- Offset-Schnitt erhalten. In diesem 2-D-Schnitt ist die horizontale (x) Koordinate die Position entlang der Linie 2, ist die vertikale Koordinate die Zeit und stellen die Daten Reflexionsamplituden dar. Durch eine gleiche Verarbeitung einer Anzahl von parallelen Linien wird ein 3-D-Datenvolumen erhalten.
Eine zweite vorläufige Verarbeitungsalternative kann verwendet werden, bei welcher die Daten entlang den seismischen Linien einer Korrektur bezüglich eines Herausbewegens aus einem Gefälle (DMO) unterzogen werden können, um eine Reihe von gemeinsamen Offset-Panelen mit demselben Aufnahme-Empfänger-Abstand zu erzeugen. Das Endergebnis von sowohl einer Stapel- als auch einer DMO- Korrektur ist, wie es Fachleuten auf dem Gebiet wohlbekannt ist, ein 3-D-Datenvolumen von Reflexionsamplituden im (x, y, t)-Raum.
Wie es Erfahrenen auf dem Gebiet vertraut sein würde, muss die Verarbeitungsrichtung nicht dieselbe wie die Richtung einer Datenerfassung sein. Weiterhin ist es möglich, seismische Daten von einer einzelnen Quelle mit einer Vielzahl von Linien von Empfängern zu erfassen. Eine solche Erfassungsgeometrie macht es möglich, ein 3-D-Bild von der Erde zu erhalten und das Stapeln und eine DMO-Korrektur in beliebigen Richtungen durchzuführen. Die vorliegende Offenbarung soll solche Datensätze enthalten.
Der Annehmlichkeit halber wird die detaillierte Beschreibung von mehreren Ausführungsbeispielen der vorliegenden Erfindung wie folgt organisiert werden: Zuerst diskutiert ein Einführungsabschnitt Vorverarbeitungsbetrachtungen. Als zweites wird eine. Übersicht über die Hauptverarbeitungsschritte für das Ausführungsbeispiel präsentiert. Als zweites wird eine detaillierte Beschreibung der Datenmanipulationsprozeduren präsentiert, die für eine effiziente Implementierung des Ausführungsbeispiels auf einem MPP nötig sind. Schließlich wird das beim Ausführungsbeispiel verwendete Abbildungsschema diskutiert. Die Organisation ist jedoch nicht als beschränkend gedacht, sondern wird lediglich dazu verwendet, eine Diskussion der verschiedenen Aspekte der vorliegenden Erfindung zu vereinfachen.

TIEFENMIGRATION

Tiefenmigrations-Ausführungsbeispiele der vorliegenden Erfindung lassen sowohl eine effiziente Implementierung des Tiefenmigrations-Abbildungsschemas im Frequenzraumbereich als auch die optionale Implementierung eines Fehlerkompensationsschemas zum Korrigieren von durch die Migration eingeführten Fehlern zu. Die Fehlerkompensation wird im Frequenzwellenzahlbereich durchgeführt, um dadurch zuzulassen, dass genaue Migrationsergebnisse in reduzierten MPP-Berechnungszeiten erhalten werden. Die Ausführungsbeispiele enthalten ein neues Lösungsschema für den innersten Teil der Migrationsberechnung, was Berechnungszeiten verglichen mit dem Stand der Technik reduziert.

Vorverarbeitungsbetrachtungen

Wie es oben angegeben ist, wird die vorläufige Verarbeitung der in Fig. 1 gezeigten seismischen Daten allgemein entweder ein Stapeln oder eine DMO-Analyse enthalten. Das, Ergebnis dieser vorläufigen Verarbeitung ist ein 3-D-Datenvolumen, das aus einer Zeitserie von Reflexionsamplituden besteht, die auf einem Oberflächengitter einer gleichgerichteten x- und einer quer gerichteten y-Richtung aufgezeichnet sind, um dadurch das in Fig. 2 gezeigte (t, x, y)-Datenvolumen 30 zu definieren. Das Datenvolumen 30 besteht aus einer Zeitreihe. (einer vertikalen Achse), die bei verschiedenen Punkten auf einem horizontalen Oberflächengitter 32 aufgezeichnet sind. Tiefenmigrations-Ausführungsbeispiele der Erfindung werden allgemein mit einem zusätzlichen Daten- Vorverarbeitungsschritt fortfahren, der auf einem Computer, wie beispielsweise dem Cray Y-MP durchgeführt wird, um die Daten in ein 3-D-Volumen in (x, y, ω) zu transformieren. Dieser Schritt enthält eine Fouriertransformation in ein Volumen 34, wobei die vertikale Achse eine Frequenz ist und die horizontalen Achsen unverändert sind. Dieser Schritt wird von Fachleuten auf dem Gebiet verstanden werden.
Als nächstes wird das Volumen 34 senkrecht zu der ω-Achse in Scheiben geteilt, um Frequenzscheiben 36 zu ergeben, um bei der Migration verwendet zu werden. Die Anzahl von Scheiben, die erzeugt werden, und die Frequenzen für jede Scheibe werden eine Funktion der durchzuführenden Analyse sein, wobei die Art der Datengruppe bzw. des Datensatzes analysiert wird und andere Faktoren, die wohlbekannt sind. Die Scheiben können optional komprimiert werden, um Speicherbelastungen und Daten-Eingabe/Ausgabe-Kommunikationsbelastungen zu minimieren. Die Datenkompression könnte ohne Beschränkung aus irgendeiner Anzahl von wohlbekannten Prozeduren bestehen, wie beispielsweise einem Darstellen der Daten mit weniger Bits.
Eine Tiefenmigration-Vorverarbeitung kann auch ein vorberechnetes Tiefenmodell für die Migration in Scheiben teilen. Das Tiefenmodell wird aus Geschwindigkeiten als Funktion einer Tiefe auf demselben (x, y)-Gitter wie das Datenvolumen bestehen, und die Tiefen und die Anzahl von Geschwindigkeitsscheiben werden auch durch den Analytiker bestimmt werden. Die Geschwindigkeitsscheibe. 38 in Fig. 2 können auch komprimiert werden. Die Frequenzscheiben und die Geschwindigkeitsscheiben werden auf eine Diskette oder eine andere Hochgeschwindigkeitsvorrichtung für eine nachfolgende Verwendung durch, den MPP geschrieben. Hierin nachfolgend soll der Ausdruck "Diskette" andere Hochgeschwindigkeitsvorrichtungen enthalten, auf welchen Daten geschrieben und von welchen sie wiedergewonnen werden können. Jeder der obigen Schritte wäre Erfahrenen auf dem Gebiet vertraut.
Eine anfängliche Zuteilung der Frequenzscheiben zu den Prozessoren im MPP wird basierend auf der Anzahl von x- und y-Punkten bezüglich der Daten, der Anzahl von bei der Migration zu verwendenden Frequenzen und des auf jedem Prozessor im MPP verfügbaren Speichers durchgeführt. Der Annehmlichkeit halber wird diese Zuteilung als Teil der Hochfahrroutine durchgeführt, die einem PE innerhalb des MPP zugeordnet ist, welches als Steuerprozessor für den MPP beim Bestimmen der Zuteilung wird. Die Zuteilung wird typischerweise ein Bestimmen der minimalen Anzahl von Prozessoren enthalten, die pro Frequenzscheibe erforderlich sind, wobei jeder Prozessor eine gleiche Anzahl von x- und y- Datenpunkten hat, bei welchem Berechnungen durchgeführt werden. Jedoch ist das vorliegende Verfahren nicht auf die Verwendung von nur dieser minimalen Anzahl oder auf ein Beibehalten dieser gleichen Anzahl pro Prozessor beschränkt, welche hierin der Annehmlichkeit halber ausgewählt sind. Auf ein Beendigen dieser Zuteilung hin kehrt der Steuerprozessor zu einem Funktionieren als Analyseprozessor zurück, und zwar auf dieselbe Weise wie alle anderen Prozessoren im MPP. Der Zweck der Zuteilung besteht im Bestimmen der Anzahl von Prozessoren, die zum Durchführen von Berechnungen für jede Frequenzscheibe erforderlich sein werden. Das Ergebnis der Zuteilung besteht darin, dass bei den Berechnungsprozeduren und Datenmanipulationen, die weiter unten diskutiert werden, jeder Prozessor einer oder mehreren Frequenzscheiben zugeordnet wird, oder alternativ dazu eine Gruppe von Prozessoren einer einzelnen Frequenzscheibe zugeordnet wird. Wenn ein einzelner Prozessor mehr als einer Frequenzscheibe zugeordnet werden kann, werden die Frequenzen derart ausgewählt, dass sie die Rechenlast auf den Prozessoren im Gleichgewicht halten. Normalerweise wird dies bedeuten, dass eine Niederfrequenz- und eine Hochfrequenz-Scheibe in Paaren dem Prozessor zugeordnet werden. Wie es Erfahrenen auf dem Gebiet vertraut sein würde, ist eine Migration von niedrigen Frequenzen weniger eine Rechenbelastung als eine Migration von hohen Frequenzen; daher hilft eine solche Paarbildung ein Ausgleichen der Rechenlast unter den Prozessoren.

Verarbeitungsübersicht

Datenfluss - Allgemein

Gemäß Fig. 3A werden bei "START" die Daten aus einer Diskette zu den PEs gelesen. Wenn die Anzahl von Frequenzscheiben pro Prozessor 42 kleiner als Eins ist, folgt die Verarbeitung den Schritten in Fig. 3B. Wenn diese Anzahl Eins oder darüber ist, dann werden die Frequenzscheiben in Gruppen bzw. Sätzen 44 gesammelt, werden Verarbeitungselementen (PEs) zugeordnet und werden in den Speicher jedes PE 46 gelesen. Wenn die PEs eine Anfangsgruppe von Frequenzen bei einer ersten Durchlaufanalyse 48 verarbeiten, wird der Prozess einer Fortführung in Richtung nach unten und einer Bildgebung (die nachfolgend in Verbindung mit Fig. 4 zu beschreiben ist) angewendet, 50. Das Ergebnis dieser Fortführung nach unten und dieser Bildgebung besteht im. Erzeugen einer Teilbildebene für jede Tiefe.
Wenn zusätzliche Frequenzen zu verarbeiten bleiben, wie es normalerweise der Fall ist, werden zusätzliche Frequenzscheiben zugeordnet und von der Diskette in jedes PE gelesen, 46. Für die Verarbeitung der zweiten und nachfolgenden Gruppe von Frequenzen, 48, aktualisiert der Prozess für ein Fortführen nach unten und eine Bildgebung 54, die zuvor erzeugte Bildebene für jede Tiefe. Diese Sequenz wiederholt sich solange, wie zusätzliche Frequenzen zu verarbeiten sind, 52. Die Migration ist beendet, wenn alle Frequenzen verarbeitet worden sind, 55.
Wenn mehr als ein Prozessor zum Halten der Daten für eine einzige Frequenzscheibe erforderlich ist, wird den Schritten in Fig. 3B gefolgt. Fig. 3B ist analog der Fig. 3A, außer dass Gruppen von Prozessoren zum Verarbeiten jeder Frequenz erforderlich sind. Jede Frequenzscheibe wird einer Gruppe von PEs zugeordnet und in der Gruppe in jedes PE gelesen, 56. Für die erste Gruppe von Frequenzscheiben, 58, wird der Prozess für eine Fortführung in Richtung nach unten und eine Bildgebung (die später in Zusammenhang mit Fig. 4 beschrieben werden, angewendet, 60. Wie oben besteht das Ergebnis dieser Fortführung in Richtung nach unten und dieser Bildgebung im Erzeugen einer Bildebene für jede Tiefe. Wenn zusätzliche Frequenzen zu verarbeiten sind, wie es normalerweise der Fall ist, werden diese zusätzlichen Frequenzen zugeordnet und in jede Gruppe von PEs gelesen, 56. Ebenso wie oben aktualisiert der Prozess für eine Fortführung nach unten und eine Bildgebung, 64, für die Verarbeitung der zweiten und nachfolgenden Gruppe von Frequenzscheiben, 58, die zuvor erzeugte Bildebene für jede Tiefe. Wenn alle Frequenzen verarbeitet worden sind, 62, wird dieser Prozess beendet, 65.
Eines der bemerkenswerten Merkmale dieser Erfindung besteht darin, dass Tiefenmigrations-Ausführungsbeispiele die Fähigkeit zum Speichern von Zwischenergebnissen haben, nachdem eine Gruppe von Frequenzen verarbeitet worden ist, 55 in Fig. 3A und 65 in Fig. 3B. Die Migration kann daher fortgeführt oder zu einer späteren Zeit neu gestartet werden, was eine Korrektur oder eine Beschränkung von Eingangsparametern und eine Analyse von Zwischenergebnissen zulässt. Die Erwünschtheit eines solchen Neustartmerkmals wird Erfahrenen auf dem gebiet vertraut, sein obwohl die kosteneffektive Implementierung eines solchen Merkmals zuvor nicht erreicht worden ist.

Datenfluss - Abwärtsfortsetzung und Bildgebung

Der Schritt für eine Abwärtsfortsetzung bzw. Fortführung nach unten und eine Bildgebung (50 und 54 in Fig. 3A, 60 und 64 in Fig. 3B) wird gemäß Fig. 4 detaillierter beschrieben. Dieser Prozess wird durch jeden Prozessor oder jede Gruppe von Prozessoren durchgeführt, d. h. 46 in Fig. 3A und 56 in Fig. 3B. Der Prozess beginnt, 70, bei einer Anfangstiefe im Unterirdischen, und zwar normalerweise bei der Oberfläche der Erde. Als nächstes werden die Geschwindigkeitsdaten für das darunter liegende Tiefenintervall von einer Diskette in eine der Prozessorgruppen gelesen, 72. (Es ist von hier an zu verstehen, dass dann, wenn auf Prozessorgruppen Bezug genommen wird, eine Gruppe aus genau einem Prozessor bestehen kann, wenn ein Prozessor eine oder mehrere Frequenzen hält.)
Weil alle Prozessorgruppen im MPP die Fortführung nach unten durchführen werden, erfordert eine effiziente Ausnutzung des MPP, dass diese Geschwindigkeitsinformation allen Prozessorgruppen effizient verfügbar gemacht wird. Beim bevorzugten Ausführungsbeispiel wird das Rundsenden von Geschwindigkeiten, 74, durch einen Prozess eines logarithmischen Fanouts durchgeführt. Dieser Schritt eines logarithmischen Fanouts ist in Fig. 5A dargestellt. Ein logarithmischer Fanout ermöglicht, dass die für die Fortführung nach unten nötige Geschwindigkeitsinformation zu allen Prozessoren schnell verbreitet wird. Eine komprimierte Geschwindigkeitsscheibe, 38, von der Diskette wird in eine Gruppe von Prozessoren geladen, die mit #1 bezeichnet ist, und dekomprimiert. Darauf folgend werden die dekomprimierten Geschwindigkeitsdaten von #1 zu einer zweiten Prozessorgruppe, #2, transferiert. Bei der nächsten Stufe sendet #1 die Geschwindigkeitsscheibe zu #3, während gleichzeitig #2 sie zu #4 sendet. Bei der nächsten Stufe sendet #1 Geschwindigkeitsinformation zu #5, sendet #2 sie zu #6, #3 zu #7 und #4 zu #8. Auf diese Weise wächst die Anzahl von Prozessoren, die die Geschwindigkeitsinformation haben, geometrisch. Dies ist als logarithmischer Fanout bekannt.
Die Bezeichnung der Prozessorgruppen in Fig. 5A dient nur illustrativen Zwecken. In Wirklichkeit kann der Transfer von Geschwindigkeitsdaten zwischen Prozessoren in derselben physikalischen Nachbarschaft innerhalb des MPP oder dem, was sonst durch die internen Kommunikationspfade des MPP eng angeschlossen sein kann, oder basierend auf einer anderen Kommunikationsbeziehung, die auf das MPP-System anwendbar ist, das verwendet wird, erfolgen. Es sollte auch verstanden werden, dass andere Schemen für eine Verbreitung der Geschwindigkeitsdaten verwendet werden könnten und die Erfindung nicht auf eine logarithmischen Fanout beschränkt ist. Dies wäre Erfahrenen auf dem Gebiet vertraut.
Bei einigen MPPs haben die Verarbeitungselemente eine gepufferte Lese- und Schreibkapazität bzw. -fähigkeit. Dieses Ausführungsbeispiel der vorliegenden Erfindung zieht einen Vorteil aus dieser Fähigkeit durch asynchrones Lesen der Geschwindigkeitsdaten aus der Diskette. Während die Geschwindigkeitsdaten für eine Tiefe zu allen Prozessoren rundgesendet werden, werden die Geschwindigkeitsdaten für das nächste Tiefenintervall von der Diskette in den ersten Prozessor gelesen. Diese. Prozedur wäre Erfahrenen auf dem Gebiet vertraut.
Bei der Beendigung des Schritts 74 der Fig. 4 wird jedes PE Frequenzscheibendaten in seinem Speicher haben, und zwar über die Zuordnungs- und Leseprozeduren des Schritts 45 in Fig. 3A und des Schritts 56 in Fig. 3B. Jedes PE wird auch die Geschwindigkeitsdaten vom Schritt 74 der Fig. 4 haben. Wenn die Anzahl von Frequenzen pro RE größer als oder gleich Eins ist, dann wird jedes PE die Frequenz- und Geschwindigkeitsscheibendaten für alle x- und y-Positionen im Datenvolumen in seinem Speicher haben. Wenn jedoch mehr als ein Prozessor für eine Frequenzscheibe nötig ist, wird jedes PE zugeordnet werden, und es wird in seinem Speicher Frequenzscheiben- und Geschwindigkeitsdaten für den vollen Bereich des Datenvolumens entlang der x-Achse zurückhalten, aber für nur einen beschränkten Bereich von y-Werten. Beispielsweise ist eine einzelne Frequenzscheibe 36 von Fig. 2 in Fig. 5B gezeigt. Die Aufteilung der Scheibe 36 ist bei 36' gezeigt, wo jeder Teil sich über die gesamte Länge der x- Achse erstreckt, aber über weniger als die gesamte Länge der y-Achse.
Der Schritt 76 in Fig. 4 für eine Fortführung nach unten besteht daraus, dass die PEs für jede Frequenzscheibe 36 die Fortführungsberechnung nach unten für ihre Frequenz ωn durchführen. Der Prozess 76 für eine Fortführung nach unten nimmt diese Daten bei der Anfangstiefe, kombiniert sie mit der Geschwindigkeitsinformation im darunter liegenden Tiefenintervall und erzeugt Frequenzdaten für den unteren Teil des Tiefenintervalls. Fig. 5C zeigt die Fortführung nach unten von der Oberfläche um ein Tiefenintervall nach unten.
Eine Frequenzscheibe 36' entsprechend der Frequenz ωn wird von der Tiefe z = 0 (92) zur Tiefe z = Δz (94) nach unten fortgeführt. Jedes PE wird die Fortführung nach unten für eine gesamte Frequenzscheibe durchführen, oder für den gesamten Teil, dem es zugeordnet worden ist, wenn mehr, als ein PE für eine jeweilige Scheibe erforderlich ist.
Die Prozedur für eine Fortführung nach unten im Schritt 76 kann irgendeine Anzahl von bekannten Prozeduren enthalten, wie es weiter unten diskutiert ist. Bei einem Ausführungsbeispiel wird eine Annäherung an die Wellengleichung verwendet was die Berechnungskomplexität des Prozesses für eine Fortführung nach unten reduziert. Die Annäherung beruht auf einer paaraxialen Form der Wellengleichung und einem Teilen in der gleichgerichteten x- und der quer gerichteten y-Richtung. Diesem folgt ein Fehlerkorrekturschema 78 in Fig. 4. Das Fehlerkorrekturschema kompensiert durch die paaraxiale Annäherung eingeführte Fehler durch eine Gitterdispersion und durch das Aufteilen der Wellengleichung in x- und y-Richtungen. Die Fehlerkorrektur wird allgemein nicht bei jeder Tiefe durchgeführt werden müssen, für welche eine Fortführung nach unten durchgeführt wird, und zwar in Abhängigkeit von den Eigenschaften bzw. Charakteristiken des Geschwindigkeitsmodells. Die Prozesse einer Fortführung nach unten und einer Fehlerkorrektur werden weiter unten diskutier.
Wie es von Fachleuten auf dem Gebiet verstanden wird, erzeugt die Summierung 80 der Beiträge von den verschiedenen Frequenzscheiben ein Tiefenbild am unteren Ende des Tiefenintervalls. Wo frühere Bilder von anderen Frequenzen existieren (als Ergebnis einer früheren Analyse bei 50 und 45 in Fig. 3 oder bei 60 und 64 in Fig. 3B), werden diese früheren Bilder mit dem bei diesem Schritt 80 erzeugten Bild kombiniert. Die Schritte 76, 78 und 80 werden von Fachleuten auf dem Gebiet manchmal als der Migrationsprozess bezeichnet.
Der Prozess einer Summierung der Frequenzkomponenten und der früheren Bilder ist in den Fig. 6A und 6B dargestellt. Figur CA zeigt eine Reihe von Teilbildern 100 im Datenvolumen, das in früheren Schritten für Tiefen von z&sub1; - zn berechnet worden ist. Das Teilbild entsprechend der Tiefe, bei welcher die letzten Berechnungen beendet worden sind, wie es in Fig. 6A als 100' gezeigt ist, wird von der. Diskette wiedergewonnen und zu einem PE transferiert, das in Fig. 6A als #1 bezeichnet ist. Wie es oben angegeben ist, kann PE #1 auch eine Gruppe von PEs sein. Fig. 6B zeigt PE #1 und andere PEs, die mit #2 bis #8 bezeichnet sind und die zuletzt eine Fortführung nach unten für eine Gruppe von Frequenzen ω&sub1; bis ω&sub8; durchgeführt haben. Das zuvor existierende Teilbild bei der Tiefe von Interesse, nämlich 100', wird zum berechneten Bild 110 im PE #1 addiert, um ein aktualisiertes Teilbild 110" zu erzeugen. Gleichermaßen werden die Bilder einer Fortführung nach unten für die Frequenzen ω&sub2; bis cos auf Prozessorgruppen #2 bis #8 (von welchen alle der Annehmlichkeit halber in Fig. 6B als 110 gezeigt sind) mit dem Teilbild 110" in einem Schema kombiniert, das das umgekehrte des Schemas des logarithmischen Fanout der Fig. 5A ist, das beim Rundsenden von Geschwindigkeiten verwendet wird, mit dem Ergebnis, dass eine neue Bildabschätzung 100''' erhalten wird. Die Bildabschätzung 100''' in Fig. 6B ist eine aktualisierte Abschätzung des Bildes 100' von Fig. 6A bei der Tiefe entsprechend der Tiefe des Bildes 100'. Wo mehr als eine Frequenzscheibe auf einen einzelnen Prozessor passt, werden die Bilder von diesen Frequenzen innerhalb des Prozessors kombiniert, bevor sie mit den Bildern auf anderen Prozessoren kombiniert werden.
Das aktualisierte Bild wird zur Diskette oder zu, einer anderen Hardwarevorrichtung 82 in Fig. 4 geschrieben. Es sollte hier verstanden werden, dass die Anzahl und die Bezeichnung der Prozessoren und der Frequenzen illustrativen Zwecken dienen. Das vorliegende Verfahren kann mit Bildaktualisierungen implementiert sein, die unter irgendeiner Anzahl von Prozessoren auftreten, und das umgekehrte Fanout-Schema kann entsprechend angepasst sein Andere Schemen können auch zum Erzeugen von Bildaktualisierungen verwendet werden. Schließlich wird eine Prüfung durchgeführt, um zu sehen, ob die maximale erwünschte Tiefe abgebildet worden ist, 84. Wenn es nicht so ist, wird der Prozess für nachfolgende Tiefenintervalle 86, 72 wiederholt.
Als Teil der Ausgleichung der Last auf die Prozessoren führt dieselbe Prozessorgruppe die folgenden Aufgaben nicht durch: (i) ein Einlesen und Rundsenden einer Geschwindigkeitsscheibe, wie es in Fig. 5A gezeigt ist, (ii) Einlesen einer zuvor berechneten Bildscheibe, wie es in Fig. 6A gezeigt ist, und (iii) Ausschreiben einer Bildscheibe zur Diskette, wie es in Fig. 6B gezeigt ist. Statt dessen werden diese Aufgaben durch unterschiedliche Prozessorgruppen durchgeführt, d. h. die Prozessorgruppen, die in den Fig. 5A, 6A und 6B als #1 bezeichnet sind, können physikalisch unterschiedlich sein. Ein Grund für diese unterschiedlichen Zuordnungen besteht im Erhöhen der Effizienz bei der Verwendung eines Prozessorspeichers. Insbesondere ist es beim vorliegenden Verfahren effizienter, wenn jeder Prozessor nur einen Teil seines Speichers derart zugeteilt hat, dass er ein Daten-Eingabe/Ausgabe-Puffer ist. Weil auf diesen Puffer nur für eine Aufgabe gleichzeitig zugegriffen werden kann, ist es effizienter, wenn jeder Prozessor entweder einer Eingabe- oder einer Ausgaberolle gleichzeitig zugeordnet ist, und daher enthält das Lastenausgleichsschema des vorliegenden Verfahrens ein Sicherstellen, dass jeder Prozessor nur einer I/O-Aufgabe gleichzeitig zugeordnet ist.

Datenmanipulation

Eine wohlverstandene Herausforderung für den Analytiker von seismischen Daten besteht im Entwickeln von Datenmanipulations- und Speicherungsschemen, die die Verarbeitungsprozeduren von Interesse auf dem Computersystem effizient und kosteneffizient implementieren, das dem Analytiker verfügbar, ist. Die vorliegende Erfindung enthält eine eindeutige Reihe von Datenmanipulationen, die einen Vorteil aus den Fähigkeiten des MPP ziehen, und die die die Kosten und die Zeit reduzieren, die zum Vervollständigen einer Analyse erforderlich sind, im Vergleich mit Techniken nach dem Stand der Technik. Diese Datenmanipulationen werden am einfachsten in Zusammenhang mit der Tabelle I diskutiert, die die Ablauffolge von Schritten zeigt, in welchen Transformations- und Transpositionsoperationen durch die vorliegende Erfindung bezüglich der Anordnung des Datenvolumens im PE-Speicher während der Tiefenmigration durchgeführt werden. Die Tabelle I zeigt auch die Ablauffolge, in welcher die Daten am Ende jeder Operation im Prozessorspeicher gespeichert werden. Diese Operationen entsprechen allgemein 76, 78 und 80 in Fig. 4, wie es oben diskutiert ist.
Die Frequenzscheibendaten werden im PE-Speicher in einem (x, y, ω)-Matrizenformat gespeichert, welches in der Tabelle I der Annehmlichkeit halber allgemein als die 1-, 2- und 3- Achsen bezeichnet ist. Weil jede Frequenzscheibe eine feste Frequenz ω enthalten wird, enthält dieses Ausführungsbeispiel der vorliegenden Erfindung eine Sequenz von Transformations- und Transpositions-Operationen in Bezug auf nur die 1- und 2-Achsen.
Das anfängliche Sortieren von Daten, das im Schritt 1 der Tabelle I gezeigt ist, resultiert aus der Fouriertransformation des Datenvolumens und der Spezifikation von Frequenzscheiben, die oben in Zusammenhang mit Fig. 2 diskutiert sind. Im Schritt 2 wird die Prozedur für eine Fortführung nach unten ausgeführt, und zwar auch wie es oben diskutiert ist. Diese Prozedur tritt für alle x- und y- Gitterpunkte auf, dies in jedem PE gespeichert sind. Wenn die Frequenzscheibe über mehr als ein PE aufgeteilt ist, wird die Prozedur für alle x- und y-Gitterpunkte im Teil in dem PE wiederholt. Der Schritt 2 in der Tabelle I entspricht 76 in Fig. 4.
Die nächsten auszuführenden Analyseprozeduren, nämlich die Prozeduren einer Fehlerkorrektur und einer Umlauffilterung, entsprechen 78 in Fig. 8. Diese Prozeduren, nämlich die Schritte 3 bis 9 der Tabelle I, sind bis zu dem Ausmaß optional, dass sie nicht bei jeder Tiefe durchgeführt werden, wie es unten diskutiert wird. Diese Prozeduren werden am effizientesten im Wellenzahlbereich bzw. in der Wellenzahldomäne ausgeführt und ihnen geht daher eine Reihe von Matrizentranspositionen und Fouriertransformationen voraus, wie es in den Schritten 3, 4 und 5 der Tabelle I gezeigt ist. Die Transpositionen sind erforderlich, weil es in einer effizienteren Verwendung jedes Prozesses im MPP resultiert, wenn alle Fouriertransformationen entlang nur der Achse - 1, 2 oder 3 in der Tabelle I - durchgeführt werden, bei welcher das PE den gesamten Bereich von Daten hat. Beim vorliegenden Verfahren treten alle solchen Transformationen entlang der 1-Achse auf, wie es folgt. Wie es oben in Zusammenhang mit der Fig. 5B diskutiert ist, hat jeder Prozessor in einer Gruppe den gesamten Bereich von x-Daten in einer Frequenzscheibe, hat aber nur einen Teil der y-Daten (d. h. alle Daten entlang der 1-Achse der Tabelle I, aber nur einen Teil der Daten entlang der 2-Achse). Jeder Prozessor kann daher die Fouriertransformation des Schritts 3 der x- Daten unabhängig von den anderen Prozessoren durchführen. Die Ausgabe dieses Schritts besteht in (kx, y, ω)-Formatdaten. Als nächstes wird eine Matrixtransposition der 1- und 2- Achsen durchgeführt, die auf Matrix-Transpositionsroutinen der MPP-Bibliothek beruht. Diese Transposition tritt unter allen Prozessoren in der Gruppe auf, mit dem Ergebnis, dass jeder Prozessor einen gesamten Bereich von y-Daten in seinem Speicher hat, aber für nur einen Teil der kx-Achse, und dadurch die Anwendung der FFT auf y-Daten in jedem Prozessor im Schritt S zulässt. Das Ergebnis vom Schritt S besteht darin, dass sowohl die x-Daten als auch die y-Daten in der Wellenzahldomäne sind. Der Schritt 6 ist die Anwendung der Fehlerkorrektur und einer optionalen Umlauffilterung, die keine Transformation der (kx, ky, ω)-Daten enthält.
Eine Verwendung der Aufteilung führt Fehler in die Wellenfeldextrapolation ein, wie es in Brown, D. L., Applications of operator separation in exploration seismology, 48 Geophysics 288 (1983) beschrieben ist und wäre Erfahrenen auf dem Gebiet vertraut. Wie es von Li, Zhiming, Compensating finite difference errors in 3-D modeling and migration, 56 Geophysics 1650 (1991) gezeigt ist, kann die Fehlerkorrektur der Wellengleichungaufteilung bis zur ersten Ordnung durch einen zusätzlichen Phasenverschiebungsausdruck bei der Wellenfeldextrapolation behandelt werden. Dieser Fehlerkorrekturschritt (Fig. 4, Schritt 78) erfordert, dass die Daten von der (x, y, ω)-Domäne zu der Frequenzwellenzahldomäne transformiert werden, weshalb die Schritte 3, 4 und 5 der Tabelle I erforderlich sind. Dieser Fehlerkorrekturschritt muss normalerweise nur einmal bei allen Tiefenextrapolationsschritten angewendet werden, wie es von Fachleuten auf dem Gebiet verstanden werden wird. Wo es große vertikale Geschwindigkeitsgradienten gibt, kann die Fehlerkorrektur häufiger angewendet werden müssen. Es soll verstanden werden, dass Korrekturen höherer Ordnung, obwohl sie bei diesem Ausführungsbeispiel nicht verwendet werden, dort implementiert werden könnten, wo es große Geschwindigkeitsgradienten in der horizontalen Richtung gibt.
Nachdem die Fehlerkorrektur angewendet ist, kann auch eine optionale Umlauffilterung angewendet werden. Diese Umlauffilterung wäre Versierten auf dem Gebiet vertraut.
Nach der Fehlerkorrektur und der Umlauffilterung werden inverse Prozeduren durchgeführt, um die (x, y, ω)-Domäne erneut zu erhalten. Die Inverse FFT im Schritt 7 verlegt die Daten in den (y, kx, ω)-Raum. Der Algorithmus für die Inverse FFT wird auch entlang der 1-Achse angewendet. Danach ermöglicht die Transposition des Schritts 8 die Inverse FFT des Schritts 9, welche die Daten in den (x, y, ω)-Raum bringt. Die Summierung über ω im Schritt 10 erzeugt das tiefenmigrierte Bild im (x, y, z)-Raum für eine einzelne Frequenz, wie es oben in Zusammenhang mit dem Schritt 80 der Fig. 4 diskutiert ist.
Ein wichtiger Vorteil der vorliegenden Erfindung besteht darin, dass das obige Datenmanipulationsschema zulässt, dass jedes PE eine direkte Lösung "im Vorbeigehen" der Wellengleichung im Schritt für eine Fortführung nach unten berechnet, was entgegengesetzt zu Verfahren nach dem Stand der Technik ist, wobei die PEs eine indirekte Nachschautabellenlösung durch Beruhen auf einer im Voraus berechneten Filterkoeffiziententabelle berechneten. Die Verfahren nach dem Stand der Technik wurden gezwungen, diese indirekte Lösung durchzuführen, und zwar aufgrund der hohen Rechenbelastung, die zuvor erforderlich war, um die Wellengleichung direkt zu lösen. Das vorliegende Schema eliminiert diese hohe Belastung und vermeidet auch diese im Voraus berechneten Tabellen, wodurch genauere Verarbeitungsergebnisse zur Verfügung gestellt werden, als sie zuvor bei seismischen Datenverarbeitungsverfahren erreicht worden sind. Diese Genauigkeit ist unabhängig von dem Ausmaß eines Gefälles, wodurch eine Hauptsorge der Verfahren nach dem Stand der Technik vermieden wird. Zusätzlich lässt das vorliegende Schema durch Vermeiden der Vorausberechnung und des Speicherns von Koeffiziententabellen auch zu, dass ein höherer Prozentsatz an Prozessorspeicher Rechenaufgaben zugeteilt wird, d. h. durch Zulassen, dass größere und/oder mehr Teile in jedem Prozessor gespeichert werden, wodurch ein erhöhter Berechnungsdurchsatz im Vergleich mit den Verfahren nach dem Stand der Technik ermöglicht wird. Tabelle I: Parallelisierung von Tiefenmigrations- Datenmanipulationsaufgaben Prozessoren sind in der 1-2-Ebene

Bildgebungsberechnungen

Die Daten werden durch Verwenden eines impliziten finiten Differenzschemas nach unten fortgeführt (Fig. 4, Schritt 76), um die paaraxiale Wellengleichung zu lösen, wie es beispielsweise in Claerbout, J. F., Imaging the earth's interior, Scientific Publications, Inc. (1985) diskutiert ist. Wie es Vertrauten mit dem Gebiet bekannt sein würde, erfordert eine direkte implizite Lösung der finiten Differenz-Wellengleichung die Lösung einer Gruppe von (nx·ny · nx·ny), wobei nx und ny die Anzahl von Punkten in den x- und y-Richtungen des Datenvolumens sind. Gegensätzlich dazu erfordert ein explizites finites Differenzschema weniger Berechnungen, ist aber numerisch instabil. Die vorliegende Erfindung erhält eine Näherungslösung für die Wellengleichung unter Verwendung eines impliziten Schemas und durch Aufteilen der Wellengleichung in den x- und y-Richtungen. Dieses Schema erfordert die sukzessive Lösung einer viel kleineren Gruppe von Matrixgleichungen der Größenordnung nx · nx und ny · ny. Die Differenz zwischen impliziten und expliziten finiten Differenzlösungsschemen wird von Fachleuten auf dem Gebiet verstanden werden.
Beide Gruppen von Matrixgleichungen, die an dem Prozess für eine Fortführung nach unten beteiligt sind, erfordern die Lösung von Gleichungen der Form:
Ax = b (1)
wobei die Matrix A tri-diagonal ist. Der Vektor x bezieht sich hier auf allgemein Unbekannte in der Gleichung und nicht nur auf die räumliche Koordinate x im Datenvolumen. Wie es wohlbekannt ist, geht die herkömmliche Lösung dieser Gleichung damit weiter, dass mit der ersten Gleichung begonnen wird und dass darauf folgend eine Reihe von Eliminationen durchgeführt wird, um zur letzten Gleichung zu kommen, und dass dann die an der Stelle erhaltene Lösung in einer umgekehrten Reihenfolge rücksubstituiert wird, um alle Gleichungen zu lösen. Diese herkömmliche Lösung arbeitet zufrieden stellend auf einem MPP, wo alle Daten für eine einzelne Frequenzscheibe auf einem einzigen Prozessor sind, weil jeder der Prozessoren unabhängig von dem anderen Prozessoren arbeiten kann. Dies ist auch der Fall für die Lösung der x-Gleichungen, weil, wie es oben beschrieben ist, alle Daten in einer Reihenfolge innerhalb des Speichers jedes Prozessors verfügbar sind. Jedoch dort, wo mehr als ein Prozessor für eine einzige Frequenzscheibe nötig ist, d. h. ein individueller Prozessor für nur einen Teil der y-Achse verantwortlich ist, wie in 36' der Fig. 5B, ist eine direkte Anwendung des herkömmlichen Verfahrens auf einem MPP ineffizient. Insbesondere wird ein Berechnungs-Engpass erzeugt, weil nur ein PE die y-Gleichungen für jede ω- Scheibe lösend sein kann, obwohl jedes PE alle x-Gleichungen lösen kann. Als Ergebnis können dann, wann immer eine Frequenzscheibe einer Gruppe von mehr als einem PE zugeordnet ist, die PEs in der. Gruppe nicht effizient genutzt werden.
Die vorliegende Erfindung ist eine Verbesserung gegenüber dem Stand der Technik diesbezüglich, dass die y-Gleichungen gelöst werden können, so dass zwei PEs nahezu immer an einer jeweiligen ω-Scheibe aktiv sind. Insbesondere während das den kleinsten y-Werten zugeordnete PE ein Eliminieren des Wellenfeldes bei diesen Werten beginnt, beginnt auch das den größten y-Werten zugeordnete PE ein Eliminieren des Wellenfeldes beginnend bei dem größten y-Wert und fortschreitend in Richtung zu kleineren y-Werten. Wie es von Fachleuten auf dem Gebiet verstanden werden wird, enthält dieser Prozess die sequentielle Elimination von unbekannten Wellenfeldparametern. Da jedes PE die Grenze seines Bereichs von zugeordneten y-Werten erreicht, lässt es Information, die die modifizierten Gleichungen beschreibt, zu dem PE durch, das für den benachbarten Bereich von y-Werten verantwortlich ist. Dieser Prozess dauert an, bis das Wellenfeld bei einem inneren Wert von y vollständig bestimmt werden kann. Wie es von Fachleuten auf dem Gebiet auch verstanden werden wird, wird die Gruppe von Gleichungen, die das Wellenfeld charakterisieren, bei diesem inneren Wert von einer unbestimmten Gruppe zu einer bestimmten Gruppe reduziert worden sein, was zulässt, dass das Wellenfeld an dieser Stelle bestimmt wird. Danach wird das Wellenfeld bei diesem Wert von y bei einer umgekehrten Elimination verwendet, um das Wellenfeld bei sowohl größeren als auch kleineren Werten von y zu bestimmen, was wiederum zulässt, dass zwei PEs gleichzeitig aktiv sind. Ein spezieller Fall entsteht dann, wenn eine ungerade Zahl von PEs den Bereich von y-Werten überspannen. In diesem Fall ist das für den zentralen y-Teil verantwortliche PE durch sich selbst für einen Teil der Berechnung aktiv. Diese "Brennen-an-beiden-Enden" Duallösungstechnik, wie sie gerade beschrieben ist, vermeidet die Transposition der ω-Scheibe, die sonst zum Platzieren des gesamten y-Bereichs innerhalb jedes PE sowie der entsprechenden inversen Transposition nötig wäre. Die "Brennen-an-beiden-Enden"-Technik könnte auch bei feinkörnigen MPPs in sowohl der x- als auch der y-Richtung zum Reduzieren von Rechenzeit verwendet werden. Der Computercode für dieses Verfahren ist im Anhang angegeben.
Zusammenfassend ist das Endprodukt des Tiefenmigrationsprozesses ein genaues migriertes Bild mit beachtlich niedrigeren Kosten als sie durch eine einfache Frequenzraummigration erhalten werden können. Dies ist ein Ergebnis einer Verwendens einer Näherungslösung für die Wellengleichung unter Verwendung der x-y-Aufteilung unter Verwendung des "Brennen-an-beiden-Enden"-Verfahrens zum Lösen einer Gruppe von tri-diagonalen Gleichungen und unter Anwendung einer Fehlerkorrektur in der (ky, kx, ω)-Domäne, wovon alles effizient auf einem MPP implementiert ist.

ZEITMIGRATION

Ein zweites und ein drittes Ausführungsbeispiel der vorliegenden Erfindung führen eine Zeitmigration auf einem MPP durch. Das zweite Ausführungsbeispiel führt eine Zeitmigration in der Frequenzwellenzahldomäne durch, was eine Verbesserung gegenüber Zeitmigrationen nach dem Stand der Technik ist, die auf die Frequenzraumdomäne begrenzt sind. Ein drittes Ausführungsbeispiel bietet eine signifikante Verbesserung bezüglich der Rechengeschwindigkeit gegenüber dem zweiten Ausführungsbeispiel, wenn die Geschwindigkeit in der Erde über relativ dicke Intervalle konstant ist. Der Unterschied zwischen diesen zwei Ausführungsbeispielen liegt in den Abbildungsberechnungen, die durchgeführt werden.

Vorverarbeitungsbetrachtungen

Auf eine Weise, die gleich dem oben diskutierten Ausführungsbeispiel einer Tiefenmigration ist, wird eine Anfangsbestimmung bezüglich der Anzahl von Prozessoren pro Gruppe durchgeführt. Es gibt jedoch Unterschiede bezüglich der Art, auf welche die Datenscheiben bestimmt werden, und bezüglich des Datenflusses. Darüber hinaus kann bei unterschiedlichen Stufen der Verarbeitung die Anzahl von Prozessoren in einer Gruppe unterschiedlich sein. Diese Ausführungsbeispiel werden weiter unten unter Verwendung einer Dreiteiligkeit analog zu der oben verwendeten diskutiert. Bei der folgenden Diskussion werden jedoch zuerst die Datenmanipulationsprozeduren diskutiert werden, die bei einer Zeitmigration beteiligt sind, gefolgt durch die Diskussion der Prozeduren des Datenflusses und der Abbildungsberechnung bzw. Bildgebungsberechnung. Die Ablauffolge vereinfacht die Diskussion der komplexeren Serie von Transpositionen und Transformationen, die bei den Zeitmigrations-Ausführungsbeispielen beteiligt sind, im Vergleich mit dem oben diskutierten Tiefenmigrations- Ausführungsbeispiel.

Datenmanipulation

Ein wichtiger Teil der Zeitmigrations-Ausführungsbeispiele ist die Manipulation der Daten und die Weise, auf welche sie im Speicher des MPP sitzen. Beide Ausführungsbeispiele beginnen mit Daten in einem (t, x, y)-Datenvolumen, wie es beispielsweise in Fig. 2 gezeigt ist. Die MPP- Zeitmigrationssequenz enthält eine Serie von Transpositionen und Transformationen analog zu denjenigen, die beim Tiefenmigrations-Ausführungsbeispiel durchgeführt werden. Jedoch dehnt sich die Transpositions- und Transformationssequenz über das hinausgehend aus, was für das Tiefenmigrations-Ausführungsbeispiel erforderlich war, weil Fouriertransformationen in sowohl zeitlichen als auch zwei räumlichen Domänen durchgeführt werden müssen, im Vergleich mit den Tiefenmigrations-Transformationen, die nur in den zwei räumlichen Domänen erfolgten. Wie es oben diskutiert ist, besteht der Zweck der Transpositionen im Sicherstellen, dass alle Transformationen entlang der 1-Achse durchgeführt werden, wenn diese Achse wie oben in Zusammenhang mit der Tabelle I definiert wurde und wie sie nachfolgend in Verbindung mit der Tabelle II verwendet wird. Der Wunsch zum Durchführen aller Transformationen entlang derselben Achse leitet sich von der Präferenz ab, dass nur ein PE an der Transformation beteiligt ist, wie es auch oben bei dem Tiefenmigrations-Ausführungsbeispiel angemerkt wurde. Die Tabelle II zeigt die Effekte der Transformations- und Transpositionsoperationen bezüglich des Aufbaus des Datenvolumens im Speicher während der verschiedenen Schritte einer Zeitmigration. Es ist zu beachten, dass die in der Tabelle II gezeigte Ablauffolge auch zulässt, dass Daten optional zur Diskette geschrieben werden, wenn es erwünscht ist, um die Analyse zu stoppen und später neu zu starten.
Die Datenmanipulation bei der Erfindung ist entwickelt, um einen Vorteil aus der Struktur des MPP zu ziehen. Anfangs sind die Daten in einem (t, x, y)-Volumen, Schritt 1. Gegensätzlich zu den in Fig. 2 für eine Tiefenmigration gezeigten allgemein horizontalen Frequenzscheiben sind diese Scheiben allgemein vertikal mit einem konstanten Wert entlang der y-Achse. Das Datenvolumen wird in y-Ebenen in Scheiben geschnitten, so dass eine Gruppe von Prozessoren eine t-x- Ebene hat. Die Daten werden im MPP-Speicher transponiert, Schritt 2, und fouriertransformiert, um (kx, t, y)-Volumen zu ergeben, Schritt 3.
Als nächstes lässt eine Sequenz von drei Transpositionen, nämlich die Schritte 4, 5 und 6 eine Fouriertransformation von der räumlichen y-Domäne zu der ky-Wellenzahldomäne zu, Schritt 7. Eine zusätzliche Transposition, nämlich der Schritt 8, ermöglicht eine Fouriertransformation von einem zeitlichen in einen Frequenzraum (ω, ky, ky), Schritt 9. Die Daten sind nun in einer Position, um zeitmigriert zu werden, Schritt 10. Wie es oben angemerkt wurde, zieht die Reihenfolge, in welcher die Transpositionen und die Fouriertransformationen durchgeführt werden, einen Vorteil aus der Struktur des MPP und der Speicherzuteilung in ihm.
Nach dem Zeitmigrationsprozess des Schritts 10 ermöglichen die Transpositionen und Transformationen die Rückkehr bzw. den Rücksprung zur räumlichen Domäne [(τ, x, y)-Raum]. Diese Schritte kehren die Ablauffolge von Operationen um, die in den Schritten 2-8 durchgeführt sind.
Wie es oben in Verbindung mit dem Tiefenmigrations- Ausführungsbeispiel der vorliegenden Erfindung angezeigt ist, lässt das obige Zeitmigrations-Datenmanipulationsschema zu, dass jedes PE eine direkte Lösung der Wellengleichung im Migrationsschritt berechnet und das Beruhen des Standes der Technik auf indirekten Nachschautabellenlösungen vermeidet, wobei beispielsweise separate Operatoren für jede Frequenz erforderlich waren. Die Zeitmigrations-Ausführungsbeispiele der vorliegenden Erfindung erhalten dadurch genauere, kosteneffektivere Ergebnisse, als sie zuvor möglich gewesen sind. Tabelle II: Parallelisierung von Zeitmigrations- Datenmanipulationsaufgaben Prozessoren sind in der 1-2-Ebene
*Optionales Schreiben zur Diskette, gefolgt durch ** optionales Lesen von der Diskette

Verarbeitungsübersicht

Datenfluss - Allgemein

Der Gesamtfluss von Daten ist in den Fig. 7A-7E angegeben. Diese Figuren richten sich allgemein auf eine Transformation der x-Koordinate zu der kx-Domäne (Fig. 7A, Schritte 156- 164), eine Transformation der y-Koordinate zu der ky-Domäne (Fig. 7B, Schritte 172-182), und eine Transformation der t- Koordinate zu der Frequenzdomäne vor einer Bildgebung, zusammen mit einer optionalen schließlichen Vorabbildungsfilterung (Schritt 7D, Schritt 182 und Fig. 7D Schritte 220-224). Die Bildgebungssequenzen sind auf dem Übrigen der Fig. 7D und in Fig. 7E umrissen, und die Nachabbildungsprozeduren, die zum Transformieren der Daten zurück zur migrierten Zeit erforderlich sind, die auf dem Gebiet allgemein τ-, x- und y-Domänen genannt werden, sind in der Fig. 7B (Schritte 184-190) und in der Fig. 7C gezeigt.
Die Daten sind anfangs in einem (t, x, y)-Kubus 150. Der erste Schritt besteht im Teilen der Daten in y-Scheiben, z. B. in Ebenen in der t-x-Domäne, 152. Gegensätzlich zu den allgemein horizontalen Frequenzscheiben, die in Fig. 2 für eine Tiefenmigration gezeigt sind, sind diese Scheiben allgemein vertikal mit einem konstanten Wert entlang der y- Achse. Eine anfängliche Zuteilung dieser Scheiben zu Prozessorgruppen wird durchgeführt, 154. Wie oben, soll verstanden werden, dass Gruppen aus einem oder mehreren Prozessoren bestehen können. Als nächstes werden die Daten geprüft, um zu sehen, ob irgendwelche Scheiben unverarbeitet bleiben, 156. Anfangs würde dies natürlich der Fall sein, und die nächste Scheibe von t-x-Daten wird in die nächste Gruppe von Prozessoren gelesen, 158.
Eine FFT wird in der x-Richtung angewendet, um y-Scheiben von t-kx-Daten zu ergeben, 160. Wie es oben diskutiert ist, ist diese FFT am effizientesten auf einem MPP implementiert, wenn ein jeweiliger Prozessor, der die Berechnung durchführt, alle Daten entlang der x-Richtung für eine spezifische y-Scheibe enthält. Daher erfordert die FFT eine Transposition der t-x- Daten in x-t-Daten, eine Anwendung der FET in der x-Richtung und eine Transposition zurück, um t-kx-Daten zu ergeben. In Abhängigkeit davon, ob das PE die Daten innerhalb seines Speichers zurückhalten kann und eine Verarbeitung fortführen kann, 162, springt eine Verarbeitung entweder zurück zu 156 oder schreibt die Daten aus zur Diskette, 164, bevor sie zurück zu 156 geht.
Nachdem alle y-Scheiben von Daten in die kx-Domäne verarbeitet worden sind, wird eine Transformation der y- Koordinate zu der ky-Domäne auf eine analoges Weise durchgeführt (Fig. 7B, Schritte 172-182). Jedoch erfordert dieser Teil der Verarbeitung, dass die (t, kx, y)-Daten in allgemein vertikale kx-Scheiben unterteilt werden, 172. Wie bei der anfänglichen, y-Scheibenbildung wird die Anzahl von kx-Scheiben und die Zuteilung von diesen Scheiben zu Prozessorgruppen, 172, von der Größe des Datensatzes (d. h. der Anzahl von Abtastungen in jeder der drei Richtungen t, kx und y), der Anzahl von verfügbaren Prozessoren und der Menge an Speicher, der pro Prozessor verfügbar ist, abhängen, wie es von Fachleuten auf dem Gebiet verstanden werden wird.
Wenn die Daten im Speicher sind, 174, dann werden sie in kx- Scheiben transponiert, 176. Ein Schritt 178 prüft, um zu sehen, ob alle kx-Scheiben verarbeitet worden sind, und wenn es so ist, wird die Endverarbeitungssequenz in Fig. 7C initiiert. Wenn irgendeine der kx-Scheiben nicht verarbeitet ist, dann wird die nächste Gruppe von kx-Scheiben in den Speicher von der Diskette gelesen, wenn es nötig ist, 180. Während des Lesens von der Diskette tritt die Transposition zu kx-Scheiben, welche zuvor im Schritt 176 für Daten im Speicher auftragen, auf. Bei dieser Stufe des Prozesses sind die Daten in der (t, y, kx)-Form. Ein Schritt 182 enthält die schließlichen Transpositionen und Transformationen zu (t, ky, kx) und dann durch eine FFT zu dem (ω, ky, kx)-Raum.
Die Daten werden dann dem Schema in den Fig. 7D und 7E folgend verarbeitet, um ein migriertes Bild im (τ, ky, kx)- Raum zu erzeugen, wobei τ die migrierte Zeit ist. Dieser Migrationsprozess ist unten bei der Diskussion der Fig. 7D und 7E angegeben. Die Prozedur wiederholt sich über die zusätzliche Scheibenverarbeitungssequenz der Schritte 188 und 178, wobei dann, wenn es nötig ist, zur Diskette geschriebene Daten y und kx transponieren, 190. Es ist zu beachten, dass die Schritte 186 und 190 die Transformation und die Transpositionen in den Schritten 11, 12, 13 und 14 in Tabelle II durchführen.
Fig. 7C umreißt die Schritte, die bei einer Transformation vom (τ, y, kx)-Raum zum (τ, x, y)-Raum beteiligt sind, und ist gleich zu und umgekehrt zu den Schritten, die bei den Transpositionen und Transformationen beteiligt sind, die bis zudem Migrationsschritt führen. Noch einmal wird das Datenvolumen in y-Scheiben aufgeteilt, 200, und die Daten werden auf Prozessorgruppen aufgeteilt, 202. Wenn die Daten im Speicher sind, 204, werden sie transponiert, um y-Scheiben zu ergeben, 206. Die Scheiben werden geprüft, um zu sehen, ob irgendwelche weiteren y-Scheiben zurückbleiben, 208; wenn es nicht so ist, wird der Prozess beendet. Wenn es jedoch zusätzliche zu verarbeitende Scheiben gibt, werden sie in die Prozessorgruppen von der Diskette als y-Scheiben gelesen, wenn es nötig ist, 210. Die Daten werden dann in den (τ, x, y)-Raum transformiert, 212. Dies erfordert eine Transposition, eine inverse FFT und eine weitere Transposition zur effizienten Implementierung auf einem MPP. Die migrierten Zeitscheiben werden zur Diskette ausgeschrieben, 214, und zwar als das Endprodukt des Migrationsprozesses.

Datenfluss - Migration

Die Fig. 7D und 7E zeigen den Migrationsprozess. Das (ω, ky, kx)-Datenvolumen von 182, Fig. 7B, wird optional in der ω- Domäne gefiltert, und das Geschwindigkeitsmodell in Schichten konstanter Geschwindigkeit geteilt, 222. Es ist zu beachten, dass diese Schichten konstanter Geschwindigkeit weniger kompliziert als das Geschwindigkeitsmodell des Tiefenmigrations-Ausführungsbeispiels sind, das oben diskutiert ist, und daher wird die Kommunikation der Geschwindigkeiten zu den Prozessoren allgemein nicht die Rundsendetechnik erfordern, die bei einer Tiefenmigration verwendet wird, obwohl diese Technik verwendet werden kann. Die erwünschte migrierte Zeitdomäne τ wird auch in Schichten aufgeteilt, 224. Wenn der Phasenverschiebungsalgorithmus beim zweiten Ausführungsbeispiel nicht zu verwenden ist, 226, sind Verarbeitungsschritte in der Fig. 7E gezeigt. Wenn der Phasenverschiebungsalgorithmus zu verwenden ist, werden die Daten geprüft, um zu sehen, ob es irgendwelche übrigen τ- Schichten gibt, die zu verarbeiten sind, 228. Wenn es nicht so ist, ist die Frequenzdomänenmigration fertig und werden die Daten zurück zu 184 in Fig. 7B gesendet. Wenn zusätzliche τ-Schichten zu verarbeiten bleiben, wird die Geschwindigkeit für diese Schicht erhalten; wenn sie unterschiedlich von der Geschwindigkeit für die vorherige τ-Schicht ist, werden die frequenzabhängigen Phasenverschiebungen berechnet, 230, und angewendet, 232, um von der migrierten Zeit τ zu einer Zeit τ + Δτ nach unten fortzufahren. Die Daten werden über ω summiert, um ein Bild bei τ + Δτ zu bilden, 234. Dies wird von 228 an für zusätzliche Δτ-Intervalle wiederholt.
Das dritte Ausführungsbeispiel führt eine rekursive Stolt- Migration durch, die in Fig. 7E gezeigt ist. Eine vollständige Diskussion des herkömmlichen Stolt-Algorithmus ist in Stolt, R. H., Migration by Fourier Transforms, 43 Geophysics 23 (1978) angegeben. Andere Verfahren sind im Stand der Technik bekannt, und die Offenbarung hierin hat nicht die Bedeutung, dass sie auf den rekursiven Stolt- Algorithmus beschränkt ist. Das dritte Ausführungsbeispiel nimmt hier an, dass die Geschwindigkeitsvariation in der vertikalen Richtung durch Schichten angenähert werden kann, innerhalb welcher die Geschwindigkeit konstant ist. Das rekursive Stolt-Verfahren ist auch in der gesetzmäßigen Erfindungsregistrierung H 482 von Berryhill, Gonzalez und Kim diskutiert, erteilt im Jahre 1988. Jedoch ist das dritte Ausführungsbeispiel von der von Berryhill et. al. offenbarten Erfindung diesbezüglich unterschiedlich, dass das Verfahren für eine Zeitmigration verwendet wird. Während der Effekt von FFT-Schritten bei der vorliegenden Erfindung gleich demjenigen ist, der von Berryhill et al. gelehrt wird, hat ihre Erfindung nicht die Matrix-Transpositionen zwischen aufeinander, folgenden Fouriertransformationen, die beim vorliegenden Verfahren beteiligt sind.
Eine Prüfung wird für irgendwelche übrigen Geschwindigkeitsschichten durchgeführt, 236. Wenn keine Geschwindigkeitsschichten zu verarbeiten bleiben, ist die Frequenzdomänenmigration fertig und eine Verarbeitung geht zurück zu 184. Wenn zusätzliche Geschwindigkeitsschichten zu verarbeiten bleiben, wird die Geschwindigkeit für diese V- Schicht erhalten, 240. Eine Stolt-Migration wird für alle τ in dieser Geschwindigkeitsschicht und darunter durchgeführt, 242. Es gibt an dieser Stelle keine explizite Fortführung nach unten. Die τ in dieser Geschwindigkeitsschicht werden gesichert, 244, was tiefere migrierte Zeiten wegwirft bzw. entfernt. Als nächstes werden die nicht migrierten Daten von der obersten Stelle der Schicht zur untersten Stelle der Schicht unter Verwendung der Schichtengeschwindigkeit nach unten fortgeführt. Dies erfordert eine Berechnung einer Phasenverschiebung, 246, gefolgt durch eine Fortführung nach unten, 248. Dies wird für zusätzliche τ-Schichten wiederholt, wobei am Ende der Vervollständigung der Frequenzwellenzahlmigration zu 184 in Fig. 7B zurückgegangen wird.

Bildgebungsberechnungen

Das Phasenverschiebungsverfahren von Gazdag, J., Wave- Equation Migration with the Phase Shift Method, 43 Geophysics 1342, (1978) ist die Basis für den Zeitmigrationsprozess beim zweiten Ausführungsbeispiel. Der seismische Null-Offset- Schnitt, p(x, t, τ = 0) für den 2-D-Fall kann als ein an der Oberfläche der Erde gemessenes Wellenfeld angesehen werden. Die variablen x, t und τ sind die horizontale Position, die Zwei-Wege-Laufzeit bzw. die migrierte Zwei-Wege-Laufzeit. Die Migration wird durch Verwenden des Ausdrucks für eine Fortführung nach unten in der Frequenzdomäne durchgeführt:
P(kx, ω, τ + Δτ) = P(kx, ω, τ)exp(-i kz v Δτ/2) (2)
wobei
kz = (4ω²/v² - k )1/2 (3),
wobei kx die horizontale Wellenzahl und v die Geschwindigkeit im Intervall Δτ ist.
Diese Erfindung verwendet die Gleichung (2) in einer Schicht durch eine Schichtenextrapolation, was k in der Gleichung (3) durch (k + k ) für das 3-D-Äquivalent der Gazdag- Ableitung ersetzt. Das zeitmigrierte Bild wird als p(x, t = 0, τ) erhalten, was von Fachleuten auf dem Gebiet der Bildgebungszustand genannt wird. In der Frequenzdomäne wird dies durch eine Summierung der verschiedenen Frequenzkomponenten in der Lösung der Gleichung (2) bei 234 der Fig. 7D erhalten. Wie es Vertrauten mit dem Gebiet bekannt sein würde, könnten alternative Lösungen für die Wellengleichung eher als die Gleichung (2) verwendet werden. Die vorliegende Implementierung der Gleichung (2) hat nicht die Bedeutung, dass sie eine Beschränkung in Bezug auf andere Mittel für eine Wellenfeldextrapolation ist.
Das dritte Ausführungsbeispiel verwendet einen rekursiven Frequenzwellenzahl-Migrationsprozess basierend auf der so genannten Stolt-Migration. Eine Stolt-Migration verwendet das Ergebnis, dass der nach unten geführte P (kx, τ, ω) in Bezug auf das nicht migrierte Wellenfeld P (kx, 0, ω) und ein Extrapolationsoperator E(τ) gegeben ist durch:
P(kx, τ, ω) = E(τ)P(kx, 0, ω) (4).
Der Extrapolationsoperator ist gegeben durch
E(τ) = exp [-i(ω² - W²)1/2 τ]h(ω²-W²) (5),
wobei h(...) die Einheitsschrittfunktion ist und
W = kx v/2 (6),
wobei v die Geschwindigkeit in der Schicht ist.
Das zeitmigrierte Bild wird erhalten als die inverse Fouriertransformation von
Ω/(Ω² + W²)1/2 P(kx, τ, Ω), wobei
Ω = ω(1 - W²/ω²)1/2 (7).
Das rekursive Schema führt zuerst die Migration mit einer Geschwindigkeit durch, die für die erste Schicht geeignet ist, was migrierte Zeiten unter dieser Schicht wegwirft. Dem folgt eine Fortführung nach unten des nicht migrierten Wellenfeldes zum unteren Ende der ersten Schicht. Die Fortführung nach unten wird durch das beim zweiten Ausführungsbeispiel angegebene Phasenverschiebungsverfahren durchgeführt. Dieses erzeugt Daten, die am unteren Ende der ersten Schicht aufgezeichnet worden wären. Die Stolt- Migration wird dann mit der Geschwindigkeit für die nächste Schicht durchgeführt, und der Prozess einer Fortführung nach unten und einer Migration wird wiederholt, bis die erwünschte migrierte Zeit erreicht worden ist.
Das rekursive Migrationsverfahren ist signifikant schneller als frühere Verfahren einer Migration, weil die Migration durch Verwenden einer FFT durchgeführt wird, was entgegengesetzt zu einer Faltung ist und die Migration und die Fortführung nach unten werden über relativ große Intervalle durchgeführt, in welchen die Geschwindigkeit konstant ist. Das zweite Ausführungsbeispiel erfordert eine Summierung des extrapolierten Wellenfeldes über alle Frequenzen, um ein migriertes Bild im Frequenzbereich bzw. in der Frequenzdomäne zu ergeben. Das dritte Ausführungsbeispiel ergibt das zeitmigrierte Bild für alle τ innerhalb der Geschwindigkeitsschicht durch eine inverse Fouriertransformation. Die bei der Fortführung nach unten verwendete Phasenverschiebung ist dieselbe beim zweiten und beim dritten Ausführungsbeispiel. Wenn die Geschwindigkeitsschichtdicke gleich der τ-Schichtdicke ist, sind die zwei Verfahren theoretisch äquivalent, obwohl eine Berechnungseffizienz in Abhängigkeit von Implementierungsbetrachtungen variieren wird, wie beispielsweise der FFT-Berechnungsroutine, die verwendet wird.

Tiefenmigrations-Beispiel

Die Leistungsfähigkeit der vorliegenden Erfindung im Vergleich mit dem Stand der Technik ist in Fig. 8 demonstriert, die die Reduzierung bezüglich der Taktzeit für eine 3-D-Nachstapel-Tiefenmigration unter Verwendung des vorliegenden Verfahrens darstellt, das auf einem T3D MPP implementiert ist, verglichen mit demjenigen für ein Verfahren nach dem Stand der Technik, das auf einem Y-MP implementiert ist. Die Fig. 8 zeigt die Taktzeit für eine relativ kleine 3-D-Übersicht, für welche eine Nachstapel- Tiefenmigration auf einem Cray-Y-MP mit 4 Prozessoren implementiert wurde. Selbst für das relativ kleine Datenvolumen von 1 · 10¹¹ waren mehr als 60 Stunden zum Durchführe der Migration erforderlich. Im Vergleich dazu kann ein Datenvolumen, das nahezu fünfmal so groß ist, auf einem T3D mit 128 Prozessoren in weniger als 20 Stunden verarbeitet werden. Ein Datenvolumen von 4 · 10¹¹ entspricht beispielsweise 100 Frequenzen einer Übersichtsgröße von 20 km · 20 km mit einer Abtastung für alle 10 Meter und 1000 Tiefen alle 5 Meter. Das Ergebnis besteht in drastisch niedrigeren Kosten für eine seismische Datenbearbeitung.
Es sollte verständen werden, dass die Erfindung nicht ungebührlich auf das Vorangehende beschränkt sein soll, das zu illustrativen Zwecken aufgezeigt worden ist. Verschiedene Modifikationen und Alternativen werden Fachleuten auf dem Gebiet klar werden, ohne vom Schutzumfang der Erfindung abzuweichen, wie sie beansprucht ist. ANHANG

Claims

1. Frequenzdomänenverfahren zum Bearbeiten seismischer Daten, wobei die seismischen Daten einer unterirdischen Region der Erde entsprechen, wobei das Bearbeiten auf einem Computersystem durchgeführt wird, welches mehrere Bearbeitungselemente hat, wobei das Verfahren folgende Schritte umfasst:

a) Bestimmen von Scheiben in den Daten;

b) Zuordnen von zumindest einem Teilbereich auf zumindest einer der Scheiben zu zumindest jedem einer Vielzahl der Bearbeitungselemente, wobei die Teilbereichszuordnung zu jedem der Vielzahl der Bearbeitungselemente unabhängig ist von der Teilbereichszuordnung zu jedem anderen von der Vielzahl der Bearbeitungselemente;

c) Vorberechnen eines Geschwindigkeitsmodells entsprechend der unterirdischen Region;

d) Durchführen einer direkten Migration des Teilbereichs in das Bearbeitungselement in jedem der Vielzahl der Bearbeitungselemente unter Verwendung des Geschwindigkeitsmodels, wobei die Migration in der Frequenzdomäne durchgeführt wird, wobei jede der Vielzahl der Bearbeitungselemente die Migration unabhängig von jedem anderen der Bearbeitungselemente durchführt, wobei die Migration ermöglicht wird durch Fourier- Transformation jedes der zwei räumlichen Parameter, welche die seismischen Daten charakterisieren, in die Wellenzahldomäne, wobei jede der Transformationen von jedem der Vielzahl von Bearbeitungselementen unabhängig von jedem anderen der Vielzahl von Bearbeitungselementen durchgeführt wird; und

e) Erzeugen einer Karte von Reflektoren in den seismischen Daten unter Verwendung der migrierten Teilbereiche.

2. Verfahren nach Anspruch 1, wobei

a) die seismischen Daten von der Zeitdomäne in die Frequenzdomäne Fourier-transformiert werden und die Scheiben Frequenzscheiben in den transformierten Daten sind;

b) sich das Geschwindigkeitsmodell bis zu einer maximalen Tiefe von Interesse in der unterirdischen Region erstreckt und unterteilt ist in eine Familie von Geschwindigkeitsscheiben;

c) für jedes der Vielzahl von Bearbeitungselementen ein erster der zugeordneten Teilbereiche zu dem Bearbeitungselement übertragen wird;

d) die Migration einschließt

(i) Übertragen einer ersten der Familie von vorberechneten Geschwindigkeitsscheiben an jedes der Vielzahl von Bearbeitungselementen;

(ii) Abwärtsfortsetzen jedes der Teilbereiche in jedem der Bearbeitungselemente unter Verwendung der übertragenen Geschwindigkeitsscheibe, wobei die Abwärtsfortsetzung eine Getrennte - Wellengleichung - Lösungstechnik einschließt, in der ein erster der zwei räumlichen Parameter für jeden der Teilbereiche durch das zugeordnete Bearbeitungselement unabhängig sowohl von jedem anderen der Teilbereiche und von jedem anderen der Bearbeitungselemente gelöst wird; und

(iii) Wiederholen der Prozeduren des Geschwindigkeitsscheiben-Übertragens und des Abwärtsfortsetzens bis die maximale Tiefe von Interesse erreicht ist;

e) Wiederholen der Teilbereichsübertragung und der Migration, bis alle Teilbereiche auf allen Scheiben migriert sind; und

f) die Karte eine Reflektortiefenkarte ist.

3. Verfahren nach Anspruch 2, wobei für einen zweiten der zwei räumlichen Parameter die Abwärtsfortsetzung für jede der Frequenzscheiben eine Von-Beiden-Seiten- Duallösungs-Technik einschließt, welche eine simultane Kalkulation der Lösung der getrennten Wellengleichung für den zweiten räumlichen Parameter für die Scheibe durch zumindest zwei Bearbeitungselemente erlaubt.

4. Verfahren nach Anspruch 2, wobei die Abwärtsfortsetzung auch eine Fehlerberichtigungskalkulation in der Wellenzahldomäne einschließt.

5. Verfahren nach Anspruch 4, wobei für jede Frequenzscheibe die unabhängigen Transformationen durch zumindest einen Austausch von Daten für einen der räumlichen Parameter zwischen den Bearbeitungselementen, zu welchen die Teilbereiche auf der Scheibe zugeordnet sind, ermöglicht werden.

6. Verfahren nach Anspruch 2, wobei die Abwärtsfortsetzung auch eine Umlauf-Filter-Kalkulation in der Wellenzahldomäne einschließt.

7. Verfahren nach Anspruch 6, wobei für jede Frequenzscheibe die unabhängigen Transformationen durch zumindest einen Austausch von Daten für einen der räumlichen Parameter zwischen den Bearbeitungselementen, zu welchen die Teilbereiche auf der Scheibe zugeordnet sind, ermöglicht werden.

8. Verfahren nach Anspruch 1, wobei

a) die bestimmten Scheiben Zeitscheiben in den seismischen Daten sind;

b) das Geschwindigkeitsmodell eine Familie von Schichten konstanter Geschwindigkeit einschließt, wobei sich das Modell zu einer maximalen Tiefe von Interesse in der unterirdischen Region erstreckt, wobei die Familie von Schichten zu jedem der Vielzahl der Bearbeitungselemente kommuniziert wird;

d) die Migration in jedem der Vielzahl der Bearbeitungselemente einschließt:

(i) Fourier-Transformieren der übertragenen Teilbereiche von der Zeitdomäne in die Frequenzdomäne und

(ii) Zeit-Migrieren des transformierten Teilbereichs unter Verwendung der Familie von Geschwindigkeitsschichten, bis die maximale Tiefe von Interesse erreicht ist;

e) die Migration wiederholt wird, bis alle Teilbereiche auf allen Scheiben migriert sind; und

f) die Karte eine Reflektorzeitkarte ist.

9. Verfahren nach Anspruch 8, wobei die unabhängigen Transformationen durch zumindest einen Austausch von Daten zwischen den Bearbeitungselementen, zu welchen die Teilbereiche auf der Scheibe zugeordnet sind, ermöglicht werden.

10. Verfahren nach Anspruch 8, wobei die direkte Migration eine Phasenverschiebungs-Abwärtsfortsetzung verwendet.

11. Verfahren nach Anspruch 8, wobei die direkte Migration eine rekursive Frequenz-Wellenzahl-Abwärtsfortsetzung verwendet.