Vielen Dank an den Autor des Blogs für seinen Bericht
Einführung
Diese Notiz wurde im Anschluss an die Konferenz verfasst und enthält eine Zusammenfassung des Berichts, den ich gegeben habe. Die Veranstaltung fand im Juni in St. Petersburg statt. Wenn man bedenkt, dass ich einen Block vom Konferenzort entfernt arbeite, konnte ich nicht anders, als an dieser Veranstaltung teilzunehmen. 2016 und 2017 habe ich mir Konferenzberichte angehört und dieses Jahr einen Vortrag gehalten. Erstens ist ein (meiner Meinung nach) interessantes Thema aufgetaucht, mit dem wir uns weiterentwickeln , und zweitens erschienen nach einem langen Studium der Gesetzgebung der Russischen Föderation zur Sanktionspolitik in dem Unternehmen, in dem ich arbeite, gleich zwei Lizenzen .
Bevor ich zum Thema meiner Rede übergehe, möchte ich die gute Organisation der Veranstaltung hervorheben. Auf der Besucherseite der Konferenz wird ein Bild der Kasaner Kathedrale verwendet. Die Kathedrale ist eine der Hauptattraktionen von St. Petersburg und vom Saal, in dem die Konferenz stattfand, sehr gut sichtbar.
Am Eingang der Staatlichen Wirtschaftsuniversität St. Petersburg wurden die Teilnehmer von Assistenten aus der Mitte der Studenten empfangen – sie ließen nicht zu, dass sie sich verirrten. Bei der Registrierung wurden kleine Souvenirs verteilt (ein Spielzeug – ein blinkender Dorn, ein Stift, Aufkleber mit Wolframsymbolen). Auch Mittag- und Kaffeepausen waren im Tagungsprogramm vorgesehen. Ich habe bereits an der Wand der Gruppe von leckerem Kaffee und Kuchen gemerkt – die Köche sind großartig. Mit diesem Einführungsteil möchte ich betonen, dass die Veranstaltung selbst, ihr Format und ihr Ort bereits positive Emotionen hervorrufen.
Der von mir und Kirill Belov erstellte Bericht trägt den Titel „Verwendung von Wolfram Mathematica zur Lösung von Problemen in der angewandten Geophysik“. Spektralanalyse seismischer Daten oder „wo alte Flüsse verliefen“. Der Inhalt des Berichts umfasst zwei Teile: Erstens die Verwendung der verfügbaren Algorithmen zum Analysieren geophysikalischer Daten, und zweitens, wie man geophysikalische Daten in Wolfram Mathematica einfügt.
Seismische Erkundung
Zunächst müssen Sie einen kurzen Ausflug in die Geophysik machen. Geophysik ist die Wissenschaft, die die physikalischen Eigenschaften von Gesteinen untersucht. Nun, da Gesteine unterschiedliche Eigenschaften haben: elektrisch, magnetisch, elastisch, gibt es entsprechende Methoden der Geophysik: elektrische Prospektion, magnetische Prospektion, seismische Prospektion ... Im Rahmen dieses Artikels werden wir nur auf die seismische Prospektion näher eingehen. Die seismische Exploration ist die wichtigste Methode zur Suche nach Öl und Gas. Die Methode basiert auf der Anregung elastischer Schwingungen und der anschließenden Aufzeichnung der Reaktion der Gesteine, aus denen das Untersuchungsgebiet besteht. Die Schwingungsanregung erfolgt an Land (mit Dynamit oder nicht-explosiven Schwingungsquellen elastischer Schwingungen) oder auf See (mit Luftgewehren). Elastische Schwingungen breiten sich im Gestein aus und werden an den Grenzen von Schichten mit unterschiedlichen Eigenschaften gebrochen und reflektiert. Reflektierte Wellen kehren zur Oberfläche zurück und werden von Geophonen an Land (normalerweise elektrodynamische Geräte, die auf der Bewegung eines in einer Spule schwebenden Magneten basieren) oder Hydrophonen im Meer (basierend auf dem piezoelektrischen Effekt) aufgezeichnet. Anhand des Zeitpunkts des Eintreffens der Wellen kann man die Tiefe der geologischen Schichten beurteilen.
Schleppausrüstung für seismische Schiffe
Das Luftgewehr regt elastische Schwingungen an
Die Wellen durchdringen das Gestein und werden von Hydrophonen aufgezeichnet
Geophysikalisches Forschungsschiff „Ivan Gubkin“ am Pier in der Nähe der Blagoweschtschenski-Brücke in St. Petersburg
Seismisches Signalmodell
Gesteine haben unterschiedliche physikalische Eigenschaften. Für die seismische Erkundung sind vor allem elastische Eigenschaften wichtig – die Ausbreitungsgeschwindigkeit elastischer Schwingungen und die Dichte. Wenn zwei Schichten die gleichen oder ähnliche Eigenschaften haben, wird die Welle die Grenze zwischen ihnen „nicht bemerken“. Unterscheiden sich die Wellengeschwindigkeiten in den Schichten, kommt es an der Schichtgrenze zu Reflexionen. Je größer der Unterschied in den Eigenschaften ist, desto intensiver ist die Reflexion. Seine Intensität wird durch den Reflexionskoeffizienten (rc) bestimmt:
Dabei ist ρ die Gesteinsdichte, ν die Wellengeschwindigkeit, 1 und 2 geben die obere und untere Schicht an.
Eines der einfachsten und am häufigsten verwendeten seismischen Signalmodelle ist das Faltungsmodell, bei dem die aufgezeichnete seismische Spur als Ergebnis der Faltung einer Folge von Reflexionskoeffizienten mit einem Sondierungsimpuls dargestellt wird:
wo s(t) — seismische Spur, d. h. alles, was während einer festgelegten Aufnahmezeit von einem Hydrophon oder Geophon aufgezeichnet wurde, w(t) - das vom Luftgewehr erzeugte Signal, n(t) - zufälliges Geräusch.
Lassen Sie uns als Beispiel eine synthetische seismische Spur berechnen. Als Ausgangssignal werden wir den Ricker-Impuls verwenden, der in der seismischen Erkundung weit verbreitet ist.
length=0.050; (*Signal lenght*)
dt=0.001;(*Sample rate of signal*)
t=Range[-length/2,(length)/2,dt];(*Signal time*)
f=35;(*Central frequency*)
wavelet=(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)];
ListLinePlot[wavelet, Frame->True,PlotRange->Full,Filling->Axis,PlotStyle->Black,
PlotLabel->Style["Initial wavelet",Black,20],
LabelStyle->Directive[Black,Italic],
FillingStyle->{White,Black},ImageSize->Large,InterpolationOrder->2]Erster seismischer Impuls
Wir werden zwei Grenzen in Tiefen von 300 ms und 600 ms festlegen und die Reflexionskoeffizienten werden Zufallszahlen sein
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["Reflection Coefficients",Black,20],
LabelStyle->Directive[Black,Italic]]Reihenfolge der Reflexionskoeffizienten
Lassen Sie uns die seismische Spur berechnen und anzeigen. Da die Reflexionskoeffizienten unterschiedliche Vorzeichen haben, erhalten wir zwei abwechselnde Reflexionen auf der seismischen Spur.
traceExamle=ListConvolve[wavelet[[1;;;;1]],rcExample];
ListPlot[traceExamle,
PlotStyle->Black,Filling->0,Frame->True,Axes->False,
PlotLabel->Style["Seismic trace",Black,20],
LabelStyle->Directive[Black,Italic]]Simulierte Strecke
Für dieses Beispiel ist eine Reservierung notwendig – in Wirklichkeit wird die Tiefe der Schichten natürlich in Metern bestimmt und die Berechnung der seismischen Spur erfolgt für den Zeitbereich. Es wäre richtiger, die Tiefen in Metern anzugeben und die Ankunftszeiten unter Kenntnis der Geschwindigkeiten in den Schichten zu berechnen. In diesem Fall habe ich die Ebenen sofort auf der Zeitachse festgelegt.
Wenn wir von Feldforschung sprechen, dann werden als Ergebnis solcher Beobachtungen eine Vielzahl ähnlicher Zeitreihen (seismische Spuren) aufgezeichnet. Wenn Sie beispielsweise einen Standort mit einer Länge von 25 km und einer Breite von 15 km untersuchen, bei dem jede Spur als Ergebnis der Arbeit eine Zelle mit den Maßen 25 x 25 Meter charakterisiert (eine solche Zelle wird als Behälter bezeichnet), enthält das endgültige Datenfeld 600000 Spuren. Bei einer Abtastzeit von 1 ms und einer Aufnahmezeit von 5 Sekunden wird die endgültige Datendatei mehr als 11 GB groß sein, und das Volumen des ursprünglichen „Rohmaterials“ kann Hunderte von Gigabyte betragen.
Wie man mit ihnen arbeitet ?
Package
Die Entwicklung des Pakets begann an der VK-Pinnwand der russischsprachigen Selbsthilfegruppe. Dank der Antworten der Community konnte sehr schnell eine Lösung gefunden werden. Und in der Folge entwickelte sich daraus eine ernsthafte Entwicklung. Dazugehörigen Es wurde sogar von Moderatoren markiert. Derzeit unterstützt das Paket die Arbeit mit den folgenden Datentypen, die in der geologischen Industrie aktiv verwendet werden:
- Import von Kartendaten in den Formaten ZMAP und IRAP
- Import von Messungen in Wells im LAS-Format
- Eingabe und Ausgabe des seismischen Dateiformats
Um das Paket zu installieren, müssen Sie den Anweisungen auf der Download-Seite des zusammengestellten Pakets folgen, d. h. Führen Sie den folgenden Code in beliebiger Reihenfolge aus :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
"https://wolfr.am/FiQ5oFih",
FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]Anschließend wird das Paket im Standardordner installiert, dessen Pfad wie folgt ermittelt werden kann:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]Als Beispiel demonstrieren wir die Hauptfunktionen des Pakets. Der Aufruf erfolgt traditionell für Pakete in der Wolfram Language:
Get["GeologyIO`"]Das Paket wird mit entwickelt . Dies ermöglicht es Ihnen, die Hauptfunktionalität des Pakets mit einer Dokumentation zu begleiten, die sich im Präsentationsformat nicht von der Dokumentation von Wolfram Mathematica selbst unterscheidet, und das Paket mit Testdateien für die erste Bekanntschaft bereitzustellen.
Eine solche Datei ist insbesondere die Datei „Marmousi.segy“ – dabei handelt es sich um ein synthetisches Modell eines geologischen Abschnitts, das vom französischen Erdölinstitut entwickelt wurde. Mithilfe dieses Modells testen Entwickler ihre eigenen Algorithmen zur Wellenfeldmodellierung, Datenverarbeitung, seismischen Spureninversion usw. Das Marmousi-Modell selbst wird im Repository gespeichert, von dem das Paket selbst heruntergeladen wurde. Um die Datei abzurufen, führen Sie den folgenden Code aus:
If[Not[FileExistsQ["Marmousi.segy"]],
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]Ergebnis importieren – SEGYData-Objekt
Das SEGY-Format beinhaltet die Speicherung verschiedener Informationen über Beobachtungen. Zum einen handelt es sich hierbei um Textkommentare. Dazu gehören Informationen über den Ort der Arbeiten, die Namen der Firmen, die die Messungen durchgeführt haben, etc. In unserem Fall wird dieser Header durch eine Anfrage mit dem TextHeader-Schlüssel aufgerufen. Hier ist eine verkürzte Textüberschrift:
Short[marmousi["TextHeader"]]„Der Marmousi-Datensatz wurde am Institut erstellt ... Minimalgeschwindigkeit von 1500 m/s und Maximalgeschwindigkeit von 5500 m/s)“
Sie können das tatsächliche geologische Modell anzeigen, indem Sie mit der Taste „Spuren“ auf die seismischen Spuren zugreifen (eine der Funktionen des Pakets besteht darin, dass bei den Schlüsseln die Groß-/Kleinschreibung nicht beachtet wird):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Modell Marmousi
Derzeit ermöglicht das Paket auch das Laden von Daten in Teilen aus großen Dateien, wodurch Dateien verarbeitet werden können, deren Größe mehrere zehn Gigabyte erreichen kann. Zu den Funktionen des Pakets gehören auch Funktionen zum Exportieren von Daten nach .segy und zum teilweisen Anhängen an das Ende der Datei.
Unabhängig davon ist die Funktionalität des Pakets bei der Arbeit mit der komplexen Struktur von .segy-Dateien erwähnenswert. Denn es ermöglicht nicht nur den Zugriff auf einzelne Traces und Header über Schlüssel und Indizes, sondern auch deren Änderung und anschließendes Schreiben in eine Datei. Viele der technischen Details der Implementierung von GeologyIO würden den Rahmen dieses Artikels sprengen und verdienen wahrscheinlich eine gesonderte Beschreibung.
Relevanz der Spektralanalyse bei der seismischen Erkundung
Durch die Möglichkeit, seismische Daten in Wolfram Mathematica zu importieren, können Sie integrierte Signalverarbeitungsfunktionen für experimentelle Daten nutzen. Da jede seismische Spur eine Zeitreihe darstellt, ist die Spektralanalyse eines der wichtigsten Werkzeuge zu ihrer Untersuchung. Als Voraussetzungen für die Analyse der Frequenzzusammensetzung seismischer Daten können wir beispielsweise Folgendes nennen:
- Verschiedene Wellentypen zeichnen sich durch unterschiedliche Frequenzzusammensetzung aus. Dadurch können Sie Nutzwellen hervorheben und Störwellen unterdrücken.
- Gesteinseigenschaften wie Porosität und Sättigung können die Frequenzzusammensetzung beeinflussen. Dadurch ist es möglich, Gesteine mit den besten Eigenschaften zu identifizieren.
- Unterschiedlich dicke Schichten verursachen Anomalien in unterschiedlichen Frequenzbereichen.
Der dritte Punkt ist der wichtigste im Kontext dieses Artikels. Nachfolgend finden Sie ein Codefragment zur Berechnung seismischer Spuren bei einer Schicht mit unterschiedlicher Dicke – ein Keilmodell. Dieses Modell wird traditionell in der seismischen Erkundung untersucht, um Interferenzeffekte zu analysieren, wenn von vielen Schichten reflektierte Wellen einander überlagern.
nx=200;(* Number of grid points in X direction*)
ny=200;(* Number of grid points in Y direction*)
T=2;(*Total propagation time*)
(*Velocity and density*)
modellv=Table[4000,{i,1,ny},{j,1,nx}];(* P-wave velocity in m/s*)
rho=Table[2200,{i,1,ny},{j,1,nx}];(* Density in g/cm^3, used constant density*)
Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}];
Table[modellv[[;;70,i]]=4500;,{i,1,200}];
(*Plotting model*)
MatrixPlot[modellv,PlotLabel->Style["Model of layer",Black,20],
LabelStyle->Directive[Black,Italic]]Modell einer Pinch-Out-Formation
Die Wellengeschwindigkeit innerhalb des Keils beträgt 4500 m/s, außerhalb des Keils 4000 m/s und die Dichte wird als konstant 2200 g/cm³ angenommen. Für ein solches Modell berechnen wir Reflexionskoeffizienten und seismische Spuren.
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 Spuren für das Keilmodell
Die in dieser Abbildung dargestellte Abfolge seismischer Spuren wird als seismischer Abschnitt bezeichnet. Wie Sie sehen, ist die Interpretation auch auf einer intuitiven Ebene möglich, da die Geometrie der reflektierten Wellen eindeutig dem zuvor angegebenen Modell entspricht. Wenn Sie die Spuren genauer analysieren, werden Sie feststellen, dass sich die Spuren von 1 bis etwa 30 nicht unterscheiden – die Reflexion vom Dach der Formation und vom Boden überlappen sich nicht. Ab der 31. Spur beginnen die Reflexionen zu interferieren. Und obwohl sich im Modell die Reflexionskoeffizienten horizontal nicht ändern, ändern die seismischen Spuren ihre Intensität, wenn sich die Dicke der Formation ändert.
Betrachten wir die Reflexionsamplitude von der oberen Grenze der Formation. Ab der 60. Route beginnt die Intensität der Reflexion zuzunehmen und erreicht bei der 70. Route ihr Maximum. Auf diese Weise manifestiert sich die Interferenz der Wellen vom Dach und Boden der Schichten, was in einigen Fällen zu erheblichen Anomalien in der seismischen Aufzeichnung führt.
ListLinePlot[GaussianFilter[Abs[traces[[All,46]]],3][[;;;;2]],
InterpolationOrder->2,Frame->True,PlotStyle->Black,
PlotLabel->Style["Amplitude of reflection",Black,20],
LabelStyle->Directive[Black,Italic],
PlotRange->All]Diagramm der Amplitude der reflektierten Welle von der Oberkante des Keils
Es ist logisch, dass bei einem Signal mit niedrigerer Frequenz Interferenzen bei großen Formationsdicken auftreten und bei einem Hochfrequenzsignal Interferenzen bei geringeren Dicken auftreten. Der folgende Codeausschnitt erstellt 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 of wavelets",Black,20],
LabelStyle->Directive[Black,Italic],
ImageSize->Large,InterpolationOrder->2]Eine Reihe von Quellsignalen mit Frequenzen von 35 Hz, 55 Hz, 85 Hz
Durch die Berechnung seismischer Spuren und das Zeichnen von Diagrammen der reflektierten Wellenamplituden können wir erkennen, dass für verschiedene Frequenzen eine Anomalie bei unterschiedlichen Formationsdicken beobachtet wird.
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["Amplitudes of reflection",Black,20],
LabelStyle->Directive[Black,Italic],
Frame->True]Diagramme der Amplituden der reflektierten Welle vom oberen Rand des Keils für verschiedene Frequenzen
Die Fähigkeit, aus den Ergebnissen seismischer Beobachtungen Rückschlüsse auf die Mächtigkeit der Formation zu ziehen, ist äußerst nützlich, da eine der Hauptaufgaben bei der Ölexploration darin besteht, die vielversprechendsten Punkte für die Verlegung einer Bohrung (d. h. die Bereiche, in denen sich die Formation befindet) zu ermitteln dicker). Darüber hinaus kann es im geologischen Abschnitt Objekte geben, deren Entstehung eine starke Änderung der Mächtigkeit der Formation verursacht. Dies macht die Spektralanalyse zu einem effektiven Werkzeug für ihre Untersuchung. Im nächsten Teil des Artikels werden wir solche geologischen Objekte genauer betrachten.
Versuchsdaten. Woher hast du sie und worauf solltest du bei ihnen achten?
Die im Artikel analysierten Materialien wurden in Westsibirien gewonnen. Die Region ist, wie wahrscheinlich ausnahmslos jeder weiß, die wichtigste Ölförderregion unseres Landes. In den 60er Jahren des letzten Jahrhunderts begann in der Region die aktive Erschließung von Lagerstätten. Die wichtigste Methode zur Suche nach Ölvorkommen ist die seismische Erkundung. Es ist interessant, sich Satellitenbilder dieses Gebiets anzusehen. In kleinem Maßstab können Sie eine große Anzahl von Sümpfen und Seen erkennen. Wenn Sie die Karte vergrößern, können Sie Cluster-Bohrlochbohrungen sehen, und wenn Sie die Karte bis zur Grenze vergrößern, können Sie auch die Lichtungen der Profile erkennen, entlang derer seismische Messungen durchgeführt wurden Es wurden Beobachtungen durchgeführt.
Satellitenbild von Yandex-Karten – Stadtgebiet Nojabrsk
Ein Netzwerk von Bohrlöchern auf einem der Felder
Ölhaltiges Gestein Westsibiriens kommt in einem breiten Tiefenbereich vor – von 1 km bis 5 km. Der Großteil der ölhaltigen Gesteine entstand in der Jura- und Kreidezeit. Die Jurazeit dürfte vielen aus dem gleichnamigen Film bekannt sein. unterschied sich deutlich vom modernen. Die Encyclopedia Britannica enthält eine Reihe von Paläokarten, die jede helogische Ära charakterisieren.
Derzeit
Jurazeit
Bitte beachten Sie, dass das Gebiet Westsibiriens in der Jurazeit eine Meeresküste war (von Flüssen und einem flachen Meer durchzogenes Land). Da das Klima angenehm war, können wir davon ausgehen, dass eine typische Landschaft dieser Zeit so aussah:
Jura-Sibirien
Bei diesem Bild sind uns nicht so sehr die Tiere und Vögel wichtig, sondern das Bild des Flusses im Hintergrund. Der Fluss ist dasselbe geologische Objekt, an dem wir zuvor Halt gemacht haben. Tatsache ist, dass sich durch die Aktivität der Flüsse gut sortierte Sandsteine ansammeln, die dann zu einem Ölreservoir werden. Diese Stauseen können eine bizarre, komplexe Form haben (wie ein Flussbett) und ihre Dicke variiert – in Ufernähe ist die Dicke gering, aber näher an der Mitte des Kanals oder in Mäanderbereichen nimmt sie zu. So liegen die im Jura entstandenen Flüsse mittlerweile in einer Tiefe von etwa drei Kilometern und sind Gegenstand einer Suche nach Ölvorkommen.
Versuchsdaten. Verarbeitung und Visualisierung
Machen wir sofort einen Vorbehalt bezüglich der im Artikel gezeigten seismischen Materialien – da die für die Analyse verwendete Datenmenge erheblich ist, ist im Text des Artikels nur ein Fragment des ursprünglichen Satzes seismischer Spuren enthalten. Dies ermöglicht es jedem, die obigen Berechnungen zu reproduzieren.
Bei der Arbeit mit seismischen Daten verwendet ein Geophysiker normalerweise spezielle Software (es gibt mehrere Branchenführer, deren Entwicklungen aktiv genutzt werden, zum Beispiel Petrel oder Paradigm), die die Analyse verschiedener Datentypen ermöglicht und über eine praktische grafische Oberfläche verfügt. Bei aller Bequemlichkeit haben solche Softwaretypen aber auch ihre Nachteile – zum Beispiel nimmt die Implementierung moderner Algorithmen in stabile Versionen viel Zeit in Anspruch und die Möglichkeiten zur Automatisierung von Berechnungen sind meist begrenzt. In einer solchen Situation ist es sehr praktisch, computermathematische Systeme und höhere Programmiersprachen zu verwenden, die den Einsatz einer breiten algorithmischen Basis ermöglichen und gleichzeitig viel Routine erfordern. Dies ist das Prinzip, nach dem in Wolfram Mathematica mit seismischen Daten gearbeitet wird. Es ist unangemessen, umfangreiche Funktionen für die interaktive Arbeit mit Daten zu schreiben – es ist wichtiger, das Laden aus einem allgemein akzeptierten Format sicherzustellen, die gewünschten Algorithmen darauf anzuwenden und sie wieder in ein externes Format hochzuladen.
Nach dem vorgeschlagenen Schema werden wir die ursprünglichen seismischen Daten laden und in anzeigen :
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]Bei den so heruntergeladenen und importierten Daten handelt es sich um die aufgezeichneten Routen auf einer Fläche von 10 mal 5 Kilometern. Wenn die Daten mit einer dreidimensionalen seismischen Vermessungsmethode gewonnen werden (Wellen werden nicht entlang einzelner geophysikalischer Profile, sondern gleichzeitig über das gesamte Gebiet aufgezeichnet), ist es möglich, seismische Datenwürfel zu erhalten. Hierbei handelt es sich um dreidimensionale Objekte, deren vertikale und horizontale Schnitte eine detaillierte Untersuchung der geologischen Umgebung ermöglichen. Im betrachteten Beispiel haben wir es mit dreidimensionalen Daten zu tun. Wir können einige Informationen aus der Textüberschrift erhalten, etwa 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 NAME DER UMFRAGE: Irgendwo in Sibirien
C 7 DATEITYP 3D SEISMISCHES VOLUMEN
C 8
C 9
C10 Z-REIHE: ERSTE 2200 M, LETZTE 2400 M
Dieser Datensatz reicht aus, um die Hauptphasen der Datenanalyse zu demonstrieren. Die Spuren in der Datei werden nacheinander aufgezeichnet und jede sieht in etwa wie in der folgenden Abbildung aus – dies ist die Verteilung der Amplituden der reflektierten Wellen entlang der vertikalen Achse (Tiefenachse).
ListLinePlot[seismic3DSEGY["traces"][[100]], InterpolationOrder -> 2,
PlotStyle -> Black, PlotLabel -> Style["Seismic trace", Black, 20],
LabelStyle -> Directive[Black, Italic], PlotRange -> All,
Frame -> True, ImageSize -> 1200, AspectRatio -> 1/5]Eine der seismischen Abschnittsspuren
Wenn Sie wissen, wie viele Spuren sich in jeder Richtung des untersuchten Bereichs befinden, können Sie ein dreidimensionales Datenarray erstellen und es mit der Funktion Image3D[] anzeigen
traces=seismic3DSEGY["traces"];
startIL=1050;EndIL=2000;stepIL=2; (*координата Х начала и конца съёмки и шаг трасс*)
startXL=1165;EndXL=1615;stepXL=2; (*координата Y начала и конца съёмки и шаг трасс*)
numIL=(EndIL-startIL)/stepIL+1; (*количество трасс по оис Х*)
numXL=(EndXL-startXL)/stepIL+1; (*количество трасс по оис Y*)
Image3D[ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}],ViewPoint->{-1, 0, 0},Background->RGBColor[0,0,0]]XNUMXD-Bild eines seismischen Datenwürfels. (Vertikale Achse – Tiefe)
Wenn die interessierenden geologischen Merkmale starke seismische Anomalien verursachen, können Visualisierungswerkzeuge mit Transparenz verwendet werden. „Unwichtige“ Bereiche der Aufnahme können unsichtbar gemacht werden, sodass nur noch Anomalien sichtbar bleiben. In Wolfram Mathematica kann dies mit erfolgen и .
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]Seismisches Datenwürfelbild mit den Funktionen Opacity[] und Raster3D[].
Wie im synthetischen Beispiel kann man auf Abschnitten des ursprünglichen Würfels einige geologische Grenzen (Schichten) mit variablem Relief identifizieren.
Das wichtigste Werkzeug der Spektralanalyse ist die Fourier-Transformation. Mit seiner Hilfe können Sie das Amplituden-Frequenz-Spektrum jeder Spur oder Gruppe von Spuren auswerten. Nach der Übertragung der Daten in den Frequenzbereich gehen jedoch Informationen darüber verloren, zu welchen Zeiten (in welchen Tiefen) sich die Frequenz ändert. Um Signaländerungen auf der Zeitachse (Tiefenachse) lokalisieren zu können, werden die gefensterte Fourier-Transformation und die Wavelet-Zerlegung verwendet. Dieser Artikel verwendet die Wavelet-Zerlegung. Die Wavelet-Analysetechnologie wurde in den 90er Jahren erstmals aktiv in der seismischen Erkundung eingesetzt. Der Vorteil gegenüber der gefensterten Fourier-Transformation wird in einer besseren Zeitauflösung gesehen.
Mithilfe des folgenden Codefragments können Sie eine der seismischen Spuren in einzelne Komponenten zerlegen:
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 decomposition",Black,32],
LabelStyle->Directive[Black,Italic],
PlotRange->All,
Frame->True]Zerlegung einer Spur in Komponenten
Um zu beurteilen, wie sich die Reflexionsenergie bei verschiedenen Wellenankunftszeiten verteilt, werden Skalogramme (analog zu einem Spektrogramm) verwendet. In der Praxis ist es in der Regel nicht erforderlich, alle Komponenten zu analysieren. Typischerweise werden tiefe, mittlere und hohe Frequenzkomponenten ausgewählt.
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->{"Time","Frequency(Hz)"},LabelStyle->Directive[Black,Bold,14],
ColorFunction->"RustTones",ImageSize->Large]Skalogramm. Funktionsergebnis
Die Wolfram Language verwendet die Funktion zur Wavelet-Transformation . Und die Anwendung dieser Funktion auf den gesamten Spurensatz erfolgt über die Funktion . Hier ist eine der Stärken von Wolfram Mathematica hervorzuheben – die Fähigkeit, Parallelisierung zu verwenden . Im obigen Beispiel ist keine Parallelisierung erforderlich – das Datenvolumen ist nicht groß, aber wenn mit experimentellen Datensätzen gearbeitet wird, die Hunderttausende von Spuren enthalten, ist dies eine Notwendigkeit.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; Nach dem Anwenden der Funktion Entsprechend den ausgewählten Frequenzen erscheinen neue Datensätze. Im obigen Beispiel sind diese Frequenzen: 38 Hz, 33 Hz, 27 Hz. Die Auswahl der Frequenzen erfolgt meist auf der Grundlage von Tests – sie erhalten effektive Karten für verschiedene Frequenzkombinationen und wählen die aus Sicht eines Geologen aussagekräftigste aus.
Wenn Sie die Ergebnisse mit Kollegen teilen oder dem Kunden zur Verfügung stellen müssen, können Sie die Funktion SEGYExport[] des GeologyIO-Pakets verwenden
outputdata=seismic3DSEGY;
outputdata["traces",1;;-1]=tracesCWD[[All,3]];
outputdata["textheader"]="Wavelet Decomposition Result";
outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]];
SEGYExport["D:result.segy",outputdata];Bei drei dieser Würfel (Niederfrequenz-, Mittelfrequenz- und Hochfrequenzkomponenten) wird normalerweise RGB-Mischung verwendet, um die Daten gemeinsam zu visualisieren. Jeder Komponente ist eine eigene Farbe zugeordnet – Rot, Grün, Blau. In Wolfram Mathematica kann dies über die Funktion erfolgen .
Das Ergebnis sind Bilder, die eine geologische Interpretation ermöglichen. Die auf dem Abschnitt aufgezeichneten Mäander ermöglichen die Abgrenzung von Paläokanälen, bei denen es sich eher um Reservoirs handelt und Ölreserven enthalten. Die Suche und Analyse moderner Analoga eines solchen Flusssystems ermöglicht es uns, die vielversprechendsten Teile der Mäander zu bestimmen. Die Kanäle selbst zeichnen sich durch dicke Schichten aus gut sortiertem Sandstein aus und sind ein gutes Ölreservoir. Gebiete außerhalb der „Spitzen“-Anomalien ähneln modernen Auenablagerungen. Ablagerungen in Überschwemmungsgebieten bestehen hauptsächlich aus tonigen Gesteinen und Bohrungen in diese Zonen werden wirkungslos sein.
RGB-Slice des Datenwürfels. In der Mitte (etwas links von der Mitte) können Sie den mäandrierenden Fluss verfolgen.
RGB-Slice des Datenwürfels. Auf der linken Seite können Sie den sich schlängelnden Fluss verfolgen.
In einigen Fällen ermöglicht die Qualität der seismischen Daten deutlich klarere Bilder. Dies hängt von der Feldarbeitsmethode und der vom Rauschunterdrückungsalgorithmus verwendeten Ausrüstung ab. In solchen Fällen sind nicht nur Fragmente von Flusssystemen sichtbar, sondern auch ganze ausgedehnte Paläoflüsse.
RGB-Mischung von drei Komponenten eines seismischen Datenwürfels (horizontale Scheibe). Tiefe ca. 2 km.
Satellitenbild der Wolga bei Saratow
Abschluss
Mit Wolfram Mathematica können Sie seismische Daten analysieren und angewandte Probleme im Zusammenhang mit der Mineralexploration lösen. Das GeologyIO-Paket macht diesen Prozess komfortabler. Die Struktur seismischer Daten ist so, dass integrierte Methoden zur Beschleunigung von Berechnungen verwendet werden können (, ,…) ist sehr effizient und ermöglicht die Verarbeitung großer Datenmengen. Dies wird zu einem großen Teil durch die Datenspeicherfunktionen des GeologyIO-Pakets erleichtert. Das Paket kann übrigens nicht nur im Bereich der angewandten seismischen Erkundung eingesetzt werden. Nahezu die gleichen Arten von Daten werden im Bodenradar und in der Seismologie verwendet. Wenn Sie Vorschläge zur Verbesserung des Ergebnisses haben, welche Signalanalysealgorithmen aus dem Wolfram Mathematica-Arsenal auf solche Daten anwendbar sind oder wenn Sie kritische Kommentare haben, wenden Sie sich bitte an uns Hinterlasse einen Kommentar.
Source: habr.com