Vielen Dank an den Blog-Autor für seinen Vortrag
Einführung
Dieser Artikel wurde im Anschluss an die Konferenz geschrieben und enthält eine Zusammenfassung des Vortrags, den ich gehalten habe. Die Veranstaltung fand im Juni in Sankt Petersburg statt. Da ich nur einen Block vom Veranstaltungsort arbeite, konnte ich es mir nicht entgehen lassen, dieses Ereignis zu besuchen. In den Jahren 2016 und 2017 habe ich Vorträge der Konferenz gehört, und in diesem Jahr trat ich selbst als Redner auf. Erstens gab es ein interessantes (wie ich finde) Thema, das wir gemeinsam mit , entwickelt haben, und zweitens hat unser Unternehmen nach einer langen Untersuchung der russischen Gesetzgebung zur Sanktionspolitik gleich zwei Lizenzen erhalten .
Bevor ich zum Thema meines Vortrags übergehe, möchte ich die gute Organisation der Veranstaltung hervorheben. Auf der Visitenkarte der Konferenz wird ein Bild der Kasaner Kathedrale verwendet. Die Kathedrale ist eine der Hauptsehenswürdigkeiten von Petersburg und ist sehr gut aus dem Saal sichtbar, in dem die Konferenz stattfand.
Am Eingang der SPbGUEU wurden die Teilnehmer von studentischen Helfern empfangen, die dafür sorgten, dass niemand sich verirrt. Bei der Registrierung gab es kleine Souvenirs (ein blinkendes Spielzeug, ein Stift, Aufkleber mit Wolfram-Symbolik). Mittagessen und Kaffeepausen waren ebenfalls im Konferenzablauf enthalten. Über den köstlichen Kaffee und die Snacks habe ich bereits in der Gruppenwand erwähnt – die Köche haben gute Arbeit geleistet. Mit diesem einleitenden Teil möchte ich betonen, dass die Veranstaltung selbst, ihr Format und der Veranstaltungsort bereits positive Emotionen hervorrufen.
Der von mir und Kirill Belov vorbereitete Vortrag trägt den Titel „Einsatz von Wolfram Mathematica zur Lösung angewandter geophysikalischer Aufgaben. Spektralanalyse seismischer Daten oder „wo flossen die alten Flüsse“. Der Inhalt des Vortrags gliedert sich in zwei Teile: erstens die Anwendung von Algorithmen, die verfügbar sind, zur Analyse geophysikalischer Daten, und zweitens, wie man geophysikalische Daten in Wolfram Mathematica einfügt.
Seismische Exploration
Zunächst ist es wichtig, einen kurzen Überblick über die Geophysik zu geben. Geophysik ist die Wissenschaft, die die physikalischen Eigenschaften von Gesteinen untersucht. Da Gesteine unterschiedliche Eigenschaften wie elektrische, magnetische und elastische besitzen, gibt es auch verschiedene geophysikalische Methoden: Elektroprospektion, Magnetprospektion, Seismik… In diesem Artikel konzentrieren wir uns jedoch näher auf die seismische Prospektion. Die seismische Prospektion ist die Hauptmethode zur Suche nach Öl und Gas. Sie basiert auf der Erzeugung elastischer Wellen und der anschließenden Registrierung der Rückantwort von den Gesteinen, die das untersuchte Gebiet bilden. Die Erzeugung dieser Wellen erfolgt entweder an Land (mit Dynamit oder nicht-explosiven vibratorischen Quellen) oder im Meer (mit pneumatischen Kanonen). Elastische Wellen breiten sich durch die Gesteinsschichten aus und werden an den Grenzen von Schichten mit verschiedenen Eigenschaften gebrochen und reflektiert. Die reflektierten Wellen kehren an die Oberfläche zurück und werden mit Geophonen an Land (in der Regel elektrodynamische Geräte, die auf der Bewegung eines in einer Spule aufgehängten Magneten basieren) oder mit Hydrofonen im Meer (die auf dem Piezoeffekt basieren) registriert. Anhand der Ankunftszeiten der Wellen lässt sich auf die Tiefen der geologischen Schichten schließen.
Das seismische Schiff zieht die Ausrüstung.
Die Druckluftpistole erzeugt elastische Wellen.
Die Wellen durchdringen die Gesteinsschichten und werden von Hydrofonen registriert.
Das Forschungsschiff für geophysikalische Erkundung „Iwan Gubkin“ liegt am Kai unter der Blagoweschtschensky-Brücke in St. Petersburg.
Modell des seismischen Signals.
Gesteine haben unterschiedliche physikalische Eigenschaften. Für die seismische Erkundung sind vor allem die elastischen Eigenschaften wichtig – die Geschwindigkeit der Ausbreitung elastischer Wellen und die Dichte. Wenn zwei Schichten die gleichen oder ähnliche Eigenschaften haben, wird die Welle die Grenze zwischen ihnen „nicht bemerken“. Wenn jedoch die Wellenfüllungen in den Schichten unterschiedlich sind, kommt es an der Grenze der Schichten zu Reflexionen. Je größer der Unterschied in den Eigenschaften ist, desto intensiver ist die Reflexion. Die Intensität wird durch den Reflexionskoeffizienten (rc) bestimmt:
wobei ρ die Dichte der Gesteine, ν die Wellen Geschwindigkeit, 1 und 2 die obere und untere Schicht bezeichnen.
Ein Modell des seismischen Signals, das besonders einfach und häufig verwendet wird, ist das Faltungssignal. Dabei wird die aufgezeichnete seismische Spur als Ergebnis der Faltung einer Reihe von Reflexionskoeffizienten mit einem sondierenden Impuls dargestellt.
wobei s(t) — die seismische Spur ist, d.h. alles, was das Hydrofon oder Geophon während eines festen Registrierungszeitraums aufgezeichnet hat, w(t) — das Signal, das durch eine Druckluftpistole erzeugt wird, n(t) — zufälliger Lärm.
Lassen Sie uns als Beispiel eine synthetische Seismogramm berechnen. Als Ausgangssignal verwenden wir den in der seismischen Exploration häufig eingesetzten Ricker-Impuls.
length=0.050; (*Signal-Länge*)
dt=0.001;(*Abtastrate des Signals*)
t=Range[-length/2,(length)/2,dt];(*Signalzeit*)
f=35;(*Zentralfrequenz*)
welle=(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)];
ListLinePlot[welle, Frame->True,PlotRange->Full,Filling->Axis,PlotStyle->Black,
PlotLabel->Style["Ursprüngliche Welle",Black,20],
LabelStyle->Directive[Black,Italic],
FillingStyle->{White,Black},ImageSize->Large,InterpolationOrder->2]Ursprünglicher seismischer Impuls
Wir legen zwei Grenzen in Tiefen von 300 ms und 600 ms fest, während die Reflexionskoeffizienten Zufallszahlen sein werden.
rcExample=ConstantArray[0,1000];
rcExample[[300]]=RandomReal[{-1,0}];
rcExample[[600]]=RandomReal[{0,1}];
ListPlot[rcExample,Filling->0,Frame->True,Axes->False,PlotStyle->Black,
PlotLabel->Style["Reflexionskoeffizienten",Black,20],
LabelStyle->Directive[Black,Italic]]Reihe der Reflexionskoeffizienten
Wir berechnen und stellen den seismischen Trace dar. Da die Reflexionskoeffizienten unterschiedliche Vorzeichen haben, erhalten wir auch auf dem seismischen Trace zwei wechselnde Reflexionen.
traceExamle=ListConvolve[wavelet[[1;;;;1]],rcExample];
ListPlot[traceExamle,
PlotStyle->Black,Filling->0,Frame->True,Axes->False,
PlotLabel->Style["Seismischer Trace",Black,20],
LabelStyle->Directive[Black,Italic]]Simulierter Trace
Für dieses Beispiel muss ich klarstellen, dass die tatsächliche Tiefe der Schichten in Metern bestimmt wird und die Berechnung des seismischen Traces im Zeitbereich erfolgt. Es wäre korrekt, die Tiefen in Metern anzugeben und die Ankunftszeiten zu berechnen, wenn die Geschwindigkeiten in den Schichten bekannt sind. In diesem Fall habe ich die Schichten direkt auf der Zeitachse festgelegt.
Wenn es um Feldforschungen geht, werden durch solche Beobachtungen eine Vielzahl ähnlicher Zeitreihen (seismischer Trassen) erfasst. Zum Beispiel, wenn man ein Gebiet von 25 km Länge und 15 km Breite untersucht, wobei jede Trasse eine Zelle von 25x25 Metern charakterisiert (solch eine Zelle nennt man Bin), enthält die finale Datensammlung 600.000 Trassen. Bei einer Zeitdiskretisierung von 1 ms und einer Aufnahmezeit von 5 Sekunden wird die endgültige Datendatei mehr als 11 GB umfassen, während das Volumen des ursprünglichen „rohen“ Materials Hunderte von Gigabyte betragen kann.
Wie man damit in ?
Das Paket
Die Entwicklung des Pakets begann mit auf der Wand der VK-Gruppe der russischsprachigen Unterstützung. Dank der Antworten der Gemeinschaft wurde eine Lösung sehr schnell gefunden, die sich in eine ernsthafte Entwicklungsarbeit verwandelte. Der entsprechende wurde sogar von Moderatoren hervorgehoben. Momentan unterstützt das Paket die Arbeit mit den folgenden Datentypen, die in der Geologiebranche aktiv verwendet werden:
- Import von kartografischen Daten im ZMAP- und IRAP-Format
- Import von Messungen in Bohrlöchern im LAS-Format
- Ein- und Ausgabe von seismischen Dateien im Format
Um das Paket zu installieren, müssen Sie den Anweisungen auf der Download-Seite des gepackten Pakets folgen, d.h. den folgenden Code in einem :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
"https://wolfr.am/FiQ5oFih",
FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]Danach wird das Paket im Standardverzeichnis installiert, dessen Pfad wie folgt abgerufen werden kann:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]Zur Veranschaulichung werden die grundlegenden Funktionen des Pakets demonstriert. Der Aufruf erfolgt traditionell für Pakete in der Wolfram-Sprache:
Get["GeologyIO`"]Das Paket wird unter Verwendung von . Dies ermöglicht es, die Hauptfunktionen des Pakets mit Dokumentation zu ergänzen, die im Format der Wolfram Mathematica-Dokumentation nicht unterscheidet, und dem Paket Testdateien für eine erste Einführung bereitzustellen.
Diese Datei, insbesondere die Datei „Marmousi.segy“, ist ein synthetisches Modell eines geologischen Schnitts, das vom französischen Institut für Erdöl entwickelt wurde. Entwickler nutzen dieses Modell, um ihre eigenen Algorithmen zur Modellierung von Wellenfeldern, Datenverarbeitung, Inversion von seismischen Spuren usw. zu testen. Das Marmousi-Modell befindet sich im Repository, von dem das Paket heruntergeladen wurde. Um die Datei zu erhalten, führen wir den folgenden Code aus:
If[Not[FileExistsQ["Marmousi.segy"]],
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]Das Ergebnis des Imports ist ein SEGYData-Objekt.
Das SEGY-Format erfordert die Speicherung verschiedener Informationen zu Beobachtungen. Zunächst gibt es Textkommentare. Hier wird Informationen über den Arbeitsort, die Namen der Unternehmen, die die Messungen durchgeführt haben, usw. aufgenommen. In unserem Fall wird dieser Header über eine Anfrage mit dem Schlüssel TextHeader aufgerufen. Hier ist ein gekürzter Textheader:
Short[marmousi["TextHeader"]]„Die Marmousi-Datenreihe wurde am Institut … minimale Geschwindigkeit von 1500 m/s und eine maximale von 5500 m/s erzeugt“
Die geologische Modellierung kann erfolgen, indem auf die seismischen Spuren mit dem Schlüssel „traces“ zugegriffen wird (eine der Funktionen des Pakets ist die Unabhängigkeit der Schlüssel von der Groß- oder Kleinschreibung):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Marmousi-Modell
Derzeit ermöglicht das Paket auch das Laden von Daten in Teilen aus großen Dateien, wodurch die Verarbeitung von Dateien, deren Größen mehrere Gigabyte erreichen können, möglich wird. Zu den Funktionen des Pakets gehören ebenfalls Exportfunktionen für Daten im .segy-Format sowie die teilweise Wiederbeschreibung am Ende der Datei.
Besonders hervorzuheben ist die Funktionalität des Pakets beim Umgang mit komplexen .segy-Dateistrukturen. Es ermöglicht nicht nur den Zugriff auf einzelne Spuren und Header über Schlüssel und Indizes, sondern auch deren Änderung mit anschließender Speicherung in der Datei. Viele technische Details zur Implementierung von GeologyIO gehen über den Rahmen dieses Artikels hinaus und verdienen wahrscheinlich eine gesonderte Beschreibung.
Die Relevanz der Spektralanalyse in der seismischen Exploration
Die Möglichkeit, seismische Materialien in Wolfram Mathematica zu importieren, ermöglicht die Nutzung der integrierten Signalverarbeitungsfunktionen für experimentelle Daten. Da jede seismische Spur eine Zeitreihe darstellt, ist die spektrale Analyse eines der grundlegenden Werkzeuge für deren Untersuchung. Zu den Voraussetzungen für die Analyse des Frequenzspektrums seismischer Daten zählen unter anderem folgende Aspekte:
- Verschiedene Wellenarten zeichnen sich durch unterschiedliche Frequenzspektren aus. Dies ermöglicht die Identifizierung nützlicher Wellen und die Unterdrückung von Störwellen.
- Eigenschaften von Gesteinen wie Porosität und Sättigung können das Frequenzspektrum beeinflussen. Dadurch können Gesteine mit besseren Eigenschaften identifiziert werden.
- Schichten unterschiedlicher Dicke verursachen Anomalien in verschiedenen Frequenzbereichen.
Der dritte Punkt ist entscheidend im Kontext dieses Artikels. Im Folgenden finden Sie einen Codeausschnitt zur Berechnung seismischer Trajektorien in einem Layer mit variabler Dicke – ein Keilmodell. Dieses Modell wird traditionell in der seismischen Exploration untersucht, um Interferenz-effekte zu analysieren, bei denen die von mehreren Schichten reflektierten Wellen übereinander liegen.
nx=200;(* Anzahl der Gitterpunkte in X-Richtung *)
ny=200;(* Anzahl der Gitterpunkte in Y-Richtung *)
T=2;(* Gesamte Propagationszeit *)
(* Geschwindigkeit und Dichte *)
modellv=Table[4000,{i,1,ny},{j,1,nx}];(* P-Wellen-Geschwindigkeit in m/s *)
rho=Table[2200,{i,1,ny},{j,1,nx}];(* Dichte in g/cm³, konstante Dichte verwendet *)
Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}];
Table[modellv[[;;70,i]]=4500;,{i,1,200}];
(* Modell plotten *)
MatrixPlot[modellv,PlotLabel->Style["Modell der Schicht",Black,20],
LabelStyle->Directive[Black,Italic]]Modell des schichtübergreifenden Layers
Die Geschwindigkeit der Wellen innerhalb des Keils beträgt 4500 m/s, außerhalb des Keils 4000 m/s, und die Dichte wurde als konstant bei 2200 g/cm³ angenommen. Für ein solches Modell berechnen wir die Reflexionskoeffizienten und seismischen Trajektorien.
rc=Table[N[(modellv[[All,i]]-PadLeft[modellv[[All,i]],201,4000][[1;;200]])/(modellv[[All,i]]+PadLeft[modellv[[All,i]],201,4500][[1;;200]])],{i,1,200}];
traces=Table[ListConvolve[wavelet[[1;;;;1]],rc[[i]]],{i,1,200}];
starttrace=10;
endtrace=200;
steptrace=10;
trasenum=Range[starttrace,endtrace,steptrace];
traserenum=Range[Length@trasenum];
tracedist=0.5;
Rotate[Show[
Reverse[Table[
ListLinePlot[traces[[trasenum[[i]]]]*50+trasenum[[i]]*tracedist,Filling->{1->{trasenum[[i]]*tracedist,{RGBColor[0.97,0.93,0.68],Black}}},PlotStyle->Directive[Gray,Thin],PlotRange->Full,InterpolationOrder->2,Axes->False,Background->RGBColor[0.97,0.93,0.68]],
{i,1,Length@trasenum}]],ListLinePlot[Transpose[{ConstantArray[45,80],Range[80]}],PlotStyle->Red],PlotRange->All,Frame->True],270Degree]Seismische Trajektorien für das Keilmodell
Die Sequenz der seismischen Linien, die in diesem Bild dargestellt ist, wird als seismischer Schnitt bezeichnet. Wie zu erkennen ist, kann ihre Interpretation auch auf intuitiver Ebene erfolgen, da die Geometrie der reflektierten Wellen eindeutig dem zuvor festgelegten Modell entspricht. Wenn man die Linien detaillierter analysiert, fällt auf, dass die Linien von der ersten bis etwa zur dreißigsten nicht voneinander abweichen — die Reflexion von der oberen und der unteren Schicht überlagert sich nicht. Ab der dreißigsten Linie beginnen die Reflexionen zu interferieren. Und obwohl im Modell die Reflexionskoeffizienten in horizontaler Richtung nicht variieren, ändern sich die seismischen Linien in ihrer Intensität mit der Dicke der Schicht.
Betrachten wir die Amplitude der Reflexion von der Oberkante der Schicht. Ab der sechzigsten Linie beginnt die Intensität der Reflexion zuzunehmen und erreicht bei der siebzigsten Linie ihren Höhepunkt. So zeigt sich die Interferenz von Wellen von der oberen und unteren Grenze der Schichten, was in einigen Fällen zu signifikanten Anomalien in der seismischen Aufzeichnung führt.
ListLinePlot[GaussianFilter[Abs[traces[[All,46]]],3][[;;;;2]],
InterpolationOrder->2,Frame->True,PlotStyle->Black,
PlotLabel->Style["Amplitude der Reflexion",Black,20],
LabelStyle->Directive[Black,Italic],
PlotRange->All]Diagramm der Amplitude der reflektierten Welle von der oberen Kante des Keils
Es ist logisch, dass bei einem niedrigfrequenten Signal die Interferenz bei größeren Schichtdicken auftritt, während sie bei hochfrequenten Signalen bei geringeren Dicken sichtbar wird. Der folgende Codeabschnitt erzeugt ein Signal mit Frequenzen von 35 Hz, 55 Hz und 85 Hz.
waveletSet=Table[(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)],
{f,{35,55,85}}];
ListLinePlot[waveletSet,PlotRange->Full,PlotStyle->Black,Frame->True,
PlotLabel->Style["Set von Wavelets",Black,20],
LabelStyle->Directive[Black,Italic],
ImageSize->Large,InterpolationOrder->2]Set von Ausgangssignalen mit Frequenzen von 35 Hz, 55 Hz, 85 Hz
Durch die Berechnung der seismischen Spuren und die Erstellung von Diagrammen der Amplitude der reflektierten Welle können wir sehen, dass für unterschiedliche Frequenzen die Anomalie bei unterschiedlichen Schichtdicken auftritt.
tracesSet=Table[ListConvolve[waveletSet[[j]][[1;;;;1]],rc[[i]]],{j,1,3},{i,1,200}];
lowFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[1]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All];
medFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[2]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All];
highFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[3]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All];
Show[lowFreq,medFreq,highFreq,PlotRange->{{0,100},All},
PlotLabel->Style["Amplituden der reflektierten Welle",Black,20],
LabelStyle->Directive[Black,Italic],
Frame->True]Diagramme der Amplituden der reflektierten Welle von der oberen Kante der Schicht für verschiedene Frequenzen
Die Möglichkeit, aus seismischen Beobachtungen Rückschlüsse auf die Schichtdicke zu ziehen, ist äußerst nützlich. Denn eine der Hauptaufgaben bei der Erschließung von Ölvorkommen besteht darin, die vielversprechendsten Punkte für die Bohrung zu bewerten (d. h. die Abschnitte, in denen die Schicht dick ist). Darüber hinaus können in der geologischen Aufschichtung Objekte auftreten, die durch ihre Genese einen plötzlichen Wechsel der Schichtdicken bedingen. Dies macht die spektrale Analyse zu einem effektiven Instrument zu ihrer Untersuchung. Im nächsten Abschnitt des Artikels werden wir solche geologischen Objekte näher betrachten.
Experimentelle Daten. Wo erhalten und wonach suchen?
Die in dem Artikel analysierten Materialien stammen aus der Region Westsibirien. Diese Region, die wohl jedem bekannt ist, ist das Hauptgebiet der Ölförderung in unserem Land. Die aktive Erschließung von Ölvorkommen begann in den 1960er Jahren. Die Hauptmethode zur Suche nach Ölfeldern ist die seismische Erkundung. Es ist interessant, Satellitenbilder dieses Gebiets zu betrachten. Im Kleinmaßstab kann man eine große Anzahl von Sümpfen und Seen erkennen, vergrößert man die Karte, sieht man die Baurampen zur Bohrung von Brunnen, und bei maximaler Vergrößerung kann man auch die Schläge erkennen, entlang derer seismische Beobachtungen durchgeführt wurden.
Satellitenbild von Yandex Maps – Gebiet der Stadt Noyabrsk
Netzwerk von Bohrstandorten in einem der Vorkommen
Die erdölhaltigen Gesteine in Westsibirien liegen in einem breiten Tiefenbereich – von 1 km bis 5 km. Der Großteil der ölhaltigen Gesteine bildete sich in der Jurazeit und der Kreidezeit. Die Jurazeit ist wohl vielen durch den gleichnamigen Film bekannt. unterscheidet sich erheblich von der heutigen. In der Britannica-Enzyklopädie gibt es eine Reihe von Paläokarten, die jede geologische Epoche charakterisieren.
Gegenwärtige Zeit
Jura
Bitte beachten Sie, dass das Gebiet Westsibirien zur Zeit des Jura eine Küstenregion darstellte (Land, das von Flüssen durchzogen und von einem flachen Meer umgeben war). Da das Klima mild war, kann man annehmen, dass die typische Landschaft dieser Zeit folgendermaßen aussah:
Sibirien im Jura
Auf diesem Bild sind für uns nicht so sehr die Tiere und Vögel wichtig, sondern die Darstellung des Flusses im Hintergrund. Der Fluss ist das geologische Objekt, auf das wir zuvor hingewiesen haben. Die Aktivität der Flüsse ermöglicht es, gut sortierte Sandsteine abzulagern, die später zu Ölreserven werden. Diese Reservoire können eine bizarre, komplexe Form haben (ähnlich wie der Flusslauf) und ihre Dicke variiert – an den Ufern ist die Dicke gering, während sie näher an der Flussmitte oder in Mäandern zunimmt. So befinden sich die im Jura entstandenen Flüsse jetzt in einer Tiefe von etwa drei Kilometern und sind Gegenstand der Ölsuche.
Experimentelle Daten. Verarbeitung und Visualisierung
Lassen Sie uns eine Anmerkung zu den in diesem Artikel gezeigten seismischen Daten machen – aufgrund des Umfangs der für die Analyse verwendeten Daten ist im Text des Artikels nur ein Fragment der ursprünglichen seismischen Datensätze enthalten. Dies ermöglicht es allen Interessierten, die genannten Berechnungen nachzuvollziehen.
Bei der Arbeit mit seismischen Daten nutzt ein Geophysiker in der Regel spezialisierte Software (einige führende Anbieter wie Petrel oder Paradigm werden aktiv eingesetzt), die es ermöglicht, verschiedene Datentypen zu analysieren und eine benutzerfreundliche grafische Schnittstelle bietet. Trotz der Benutzerfreundlichkeit haben solche Softwarelösungen auch ihre Nachteile – beispielsweise benötigt die Implementierung modernster Algorithmen in stabile Versionen viel Zeit, und die Automatisierung von Berechnungen ist in der Regel eingeschränkt. In dieser Situation erweist sich der Einsatz von Computersystemen und Hochsprachen als sehr vorteilhaft, da sie eine breite algorithmische Basis bieten und gleichzeitig viele Routinearbeiten übernehmen. So funktioniert die Arbeit mit seismischen Daten auch in Wolfram Mathematica. Es ist nicht sinnvoll, eine reichhaltige Funktionalität für die interaktive Datenverarbeitung zu entwickeln – wichtiger ist es, den Import aus gängigen Formaten zu ermöglichen, die gewünschten Algorithmen anzuwenden und die Daten dann wieder in ein externes Format zu exportieren.
Gemäß dem vorgeschlagenen Schema laden wir die originalen seismischen Daten herunter und stellen sie dar in :
Get["GeologyIO`"]
seismic3DZipPath = "seismic3D.zip";
seismic3DSEGYPath = "seismic3D.sgy";
If[FileExistsQ[seismic3DZipPath], DeleteFile[seismic3DZipPath]];
If[FileExistsQ[seismic3DSEGYPath], DeleteFile[seismic3DSEGYPath]];
URLDownload["https://wolfr.am/FiQIuZuH", seismic3DZipPath];
ExtractArchive[seismic3DZipPath];
seismic3DSEGY = SEGYImport[seismic3DSEGYPath]Die so heruntergeladenen und importierten Daten sind Strecken, die auf einem Areal von 10 mal 5 Kilometern registriert wurden. Wenn die Daten nach der Methode der dreidimensionalen Seismik erhalten wurden (die Registrierung der Wellen erfolgt nicht entlang einzelner geophysikalischer Profile, sondern über das gesamte Gebiet gleichzeitig), ist es möglich, Würfel seismischer Daten zu erhalten. Diese dreidimensionalen Objekte, deren vertikale und horizontale Schnitte eine detaillierte Untersuchung des geologischen Umfelds ermöglichen. Im betrachteten Beispiel geht es um dreidimensionale Daten. Einige Informationen können wir aus dem Textheader entnehmen, beispielsweise so
StringPartition[seismic3DSEGY["textheader"], 80] // TableFormC 1 DIES IST EINE DEMO-DATEI FÜR DEN GEOLOGYIO-PAKETTEST
C 2
C 3
C 4
C 5 DATUM BENUTZERNAME: WOLFRAM BENUTZER
C 6 UNTERSUCHUNGSNAME: IRGENDWO IN SIBIRIEN
C 7 DATEITYP 3D-SEISMISCHES VOLUMEN
C 8
C 9
C10 Z-BEREICH: ERSTER 2200M LETZTER 2400M
Dieser Datensatz reicht aus, um die grundlegenden Schritte der Datenanalyse zu demonstrieren. Die Traces in der Datei sind sequenziell aufgezeichnet und jede von ihnen sieht ungefär so aus wie in der folgenden Abbildung – eine Darstellung der Amplituden reflektierter Wellen entlang der vertikalen Achse (Tiefe).
ListLinePlot[seismic3DSEGY["traces"][[100]], InterpolationOrder -> 2,
PlotStyle -> Black, PlotLabel -> Style["Seismischer Trace", Black, 20],
LabelStyle -> Directive[Black, Italic], PlotRange -> All,
Frame -> True, ImageSize -> 1200, AspectRatio -> 1/5]Einer der Traces des seismischen Profils
Indem man weiß, wie viele Traces in jede Richtung des untersuchten Bereichs vorhanden sind, kann man ein dreidimensionales Datenarray bilden und es mit der Funktion Image3D[] visualisieren.
traces=seismic3DSEGY["traces"];
startIL=1050; EndIL=2000; stepIL=2; (*X-Koordinate des Anfangs und des Endes der Aufnahme sowie Schritt der Traces*)
startXL=1165; EndXL=1615; stepXL=2; (*Y-Koordinate des Anfangs und des Endes der Aufnahme sowie Schritt der Traces*)
numIL=(EndIL-startIL)/stepIL+1; (*Anzahl der Traces entlang der X-Achse*)
numXL=(EndXL-startXL)/stepIL+1; (*Anzahl der Traces entlang der Y-Achse*)
Image3D[ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}],ViewPoint->{-1, 0, 0},Background->RGBColor[0,0,0]]Dreidimensionale Darstellung eines Würfels seismischer Daten. (Vertikale Achse – Tiefe)
Wenn geologische Objekte von Interesse starke seismische Anomalien erzeugen, können Visualisierungswerkzeuge mit Transparenz eingesetzt werden. »Unwichtige« Teile der Aufzeichnung können unsichtbar gemacht werden, wobei nur die Anomalien sichtbar bleiben. In Wolfram Mathematica lässt sich dies mit Hilfe von und .
data = ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}];
Graphics3D[{Opacity[0.1], Raster3D[data, ColorFunction->"RainbowOpacity"]},
Boxed->False, SphericalRegion->True, ImageSize->840, Background->None]Darstellung eines Würfels seismischer Daten unter Verwendung der Funktionen Opacity[] und Raster3D[]
Wie im synthetischen Beispiel können auf den Schnitten des originalen Würfels einige geologische Grenzen (Schichten) mit variabler Topografie hervorgehoben werden.
Das Hauptwerkzeug der spektralen Analyse ist die Fourier-Transformation. Damit lässt sich das Amplituden-Frequenz-Spektrum jeder Spur oder Gruppe von Spuren bewerten. Allerdings geht beim Übertragen der Daten in den Frequenzbereich die Information darüber verloren, wann (sprich auf welchen Tiefen) sich die Frequenz ändert. Um die Veränderungen des Signals auf der zeitlichen (tiefen) Achse lokalisieren zu können, werden das Fenster-Fourier-Transformation und die Wavelet-Zerlegung eingesetzt. In diesem Artikel wird die Wavelet-Zerlegung verwendet. Die Wavelet-Analyse-Technologie fand in den 90er Jahren der letzten Jahrhunderts in der seismischen Exploration Anwendung. Ein Vorteil gegenüber der Fenster-Fourier-Transformation ist die bessere zeitliche Auflösung.
Mit dem folgenden Code-Schnipsel kann die Zerlegung einer der seismischen Spuren in einzelne Komponenten durchgeführt werden:
cwd=ContinuousWaveletTransform[seismicSection["traces"][[100]]]
Show[
ListLinePlot[Re[cwd[[1]]],PlotRange->All],
ListLinePlot[seismicSection["traces"][[100]],
PlotStyle->Black,PlotRange->All],ImageSize->{1500,500},AspectRatio->Full,
PlotLabel->Style["Wavelet-Zerlegung",Black,32],
LabelStyle->Directive[Black,Italic],
PlotRange->All,
Frame->True]Zerlegung einer Spur in Komponenten
Um zu beurteilen, wie die Energie der Reflexion über verschiedene Ankunftszeiten der Wellen verteilt ist, werden Skalenogramme (ähnlich wie Spektrogramme) verwendet. In der Praxis ist es in der Regel nicht erforderlich, alle Komponenten zu analysieren. Normalerweise wählt man die niedrig-, mittel- und hochfrequenten Anteile aus.
freq=(500/(#*contWD["Wavelet"]["FourierFactor"]))&/@(Thread[{Range[contWD["Octaves"]],1}]/.contWD["Scales"])//Round;
ticks=Transpose[{Range[Length[freq]],freq}];
WaveletScalogram[contWD,Frame->True,FrameTicks->{{ticks,Automatic},Automatic},FrameTicksStyle->Directive[Orange,12],
FrameLabel->{"Zeit","Frequenz(Hz)"},LabelStyle->Directive[Black,Bold,14],
ColorFunction->"RustTones",ImageSize->Large]Skalogramm. Ergebnis der Funktion
In der Wolfram-Sprache wird für die Wavelet-Transformation die Funktion . Die Anwendung dieser Funktion auf das gesamte Set von Traces erfolgt mithilfe der Funktion . Hier ist zu beachten, dass eine der Stärken von Wolfram Mathematica die Möglichkeit ist, Parallelisierung zu verwenden. . In dem gegebenen Beispiel ist keine Parallelisierung erforderlich – das Datenvolumen ist nicht groß, aber beim Arbeiten mit experimentellen Datensätzen, die Hunderte von Tausenden von Traces enthalten, ist dies unerlässlich.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; Nach der Anwendung der Funktion entstehen neue Datensätze, die den ausgewählten Frequenzen entsprechen. Im oben genannten Beispiel sind dies Frequenzen: 38 Hz, 33 Hz, 27 Hz. Die Auswahl der Frequenzen erfolgt meist auf der Grundlage von Tests – es werden Ergebnisse für verschiedene Frequenzkombinationen ermittelt und die informativste aus geologischer Sicht ausgewählt.
Wenn die Ergebnisse mit Kollegen geteilt oder einem Kunden zur Verfügung gestellt werden müssen, kann die Funktion SEGYExport[] aus dem GeologyIO-Paket verwendet werden.
outputdata=seismic3DSEGY;
outputdata["traces",1;;-1]=tracesCWD[[All,3]];
outputdata["textheader"]="Wavelet-Zerlegungsresultat";
outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]];
SEGYExport["D:result.segy",outputdata];Verfügt man über drei solcher Würfel (niedrigfrequente, mittel- und hochfrequente Komponenten), wird in der Regel RGB-Mischung für die gemeinsame Visualisierung der Daten verwendet. Jeder Komponente wird eine eigene Farbe zugewiesen – rot, grün, blau. In Wolfram Mathematica kann dies mit der Funktion .
Das Ergebnis sind Bilder, die eine geologische Interpretation ermöglichen. Die Mäander, die im Schnitt sichtbar sind, helfen dabei, die paläo-flüsse zu umreißen, die mit höherer Wahrscheinlichkeit Reservoire darstellen und Ölreserven enthalten. Die Suche und Analyse moderner Analogien eines solchen Flusssystems ermöglicht es, die vielversprechendsten Teile der Mäander zu bestimmen. Die eigentlichen Flüsse zeichnen sich durch mächtige Schichten gut sortierten Sandsteins aus und sind ein hervorragendes Reservoir für Öl. Bereiche außerhalb der 'Schnurdrossel'-Anomalien ähneln modernen Auenablagerungen. Auenablagerungen bestehen hauptsächlich aus tonigen Gesteinen, und das Bohren in diese Zonen wird ineffektiv sein.
RGB-Schnitt des Datenwürfels. In der Mitte (etwas links von der Mitte) kann man einen mäandernden Fluss verfolgen.
RGB-Schnitt des Datenwürfels. In der linken Hälfte lässt sich ein mäandernder Fluss nachvollziehen.
In einigen Fällen ermöglicht die Qualität der seismischen Daten klarere Abbildungen. Dies hängt von der Methodik der Feldarbeiten, der verwendeten Ausrüstung und den Algorithmen zur Geräuschunterdrückung ab. In solchen Fällen sind nicht nur Fragmente von Flusssystemen sichtbar, sondern auch ganze ausgedehnte Paläoflüsse.
RGB-Mischung der drei Komponenten des seismischen Datenwürfels (horizontale Schnittansicht). Tiefe etwa 2 km.
Satellitenbild des Flusses Wolga in der Nähe von Saratow
Fazit
Mit Wolfram Mathematica können seismische Daten analysiert und praktische Aufgaben im Zusammenhang mit der Suche nach Bodenschätzen gelöst werden, während das Paket GeologyIO diesen Prozess erleichtert. Die Struktur der seismischen Daten ist so beschaffen, dass die integrierten Methoden zur Berechnungsgeschwindigkeitssteigerung verwendet werden können (, ,…) ist äußerst effektiv und ermöglicht die Verarbeitung großer Datenmengen. Dies wird teilweise durch die Speichereigenschaften des GeologyIO-Pakets unterstützt. Übrigens kann das Paket nicht nur im Bereich der angewandten seismischen Exploration verwendet werden. Praktisch dieselben Datentypen finden sich auch in der Georadar- und Seismologieforschung. Wenn Sie Vorschläge zur Verbesserung der Ergebnisse haben, welche Signalverarbeitungsalgorithmen aus dem Arsenal von Wolfram Mathematica auf solche Daten anwendbar sind oder wenn Sie kritische Anmerkungen haben, hinterlassen Sie bitte Kommentare.
Quelle: habr.com
