Wolfram Mathematica in Geophysik

Vielen Dank an den Blog-Autor Anton Ekimenko für seinen Vortrag

Einführung

Diese Notiz wurde im Anschluss an die russische Technologiekonferenz von Wolfram verfasst und enthält eine Zusammenfassung des Berichts, den ich gegeben habe. Die Veranstaltung fand im Juni in St. Petersburg statt. Angesichts der Tatsache, dass ich einen Block vom Konferenzort entfernt arbeite, konnte ich nicht anders, als an dieser Veranstaltung teilzunehmen. In den Jahren 2016 und 2017 habe ich mir die Konferenzberichte angehört und dieses Jahr eine Präsentation gehalten. Erstens ist ein interessantes (wie mir scheint) Thema aufgetaucht, das wir gemeinsam mit Kirill Belov entwickeln , und zweitens sind nach einer langen Untersuchung der Gesetzgebung der Russischen Föderation in Bezug auf die Sanktionspolitik zwei Wolfram Mathematica- Lizenzen in dem Unternehmen erschienen, in dem ich arbeite.

Bevor ich zum Thema meiner Präsentation übergehe, möchte ich auf die gute Organisation der Veranstaltung hinweisen. Die Konferenzseite verwendet das Bild der Kasaner Kathedrale. Die Kathedrale ist eine der Hauptattraktionen von St. Petersburg und von der Halle, in der die Konferenz stattfand, sehr gut sichtbar.



Am Eingang zur Staatlichen Wirtschaftsuniversität St. Petersburg wurden die Teilnehmer von studentischen Hilfskräften empfangen - sie ließen sich nicht verirren. Bei der Registrierung wurden kleine Souvenirs ausgegeben (ein Spielzeug - blinkende Adhäsionen, ein Stift, Aufkleber mit Wolfram-Symbolen). Mittag- und Kaffeepausen waren ebenfalls im Konferenzplan enthalten. Über leckeren Kaffee und Kuchen habe ich bereits an der Wand der Gruppe vermerkt - gut gemachte Köche. Mit diesem Einführungsteil möchte ich betonen, dass die Veranstaltung selbst, ihr Format und ihr Veranstaltungsort bereits positive Emotionen hervorrufen.

Der Bericht, der von mir und Kirill Belov erstellt wurde, heißt „Mit Wolfram Mathematica Probleme der angewandten Geophysik lösen. Spektralanalyse von seismischen Daten oder "wo die alten Flüsse flossen". Der Inhalt des Berichts besteht aus zwei Teilen: Erstens der Verwendung von in Wolfram Mathematica verfügbaren Algorithmen zur Analyse geophysikalischer Daten und zweitens der Platzierung geophysikalischer Daten in Wolfram Mathematica.

Seismische Erforschung

Zuerst müssen Sie einen kleinen Ausflug in die Geophysik machen. Geophysik ist eine Wissenschaft, die die physikalischen Eigenschaften von Gesteinen untersucht. Nun, da die Gesteine ​​unterschiedliche Eigenschaften haben: elektrisch, magnetisch, elastisch, gibt es entsprechende Methoden der Geophysik: elektrische Exploration, magnetische Exploration, seismische Exploration ... Im Zusammenhang mit diesem Artikel werden wir nur die seismische Exploration genauer behandeln. Seismische Exploration ist die Hauptmethode zum Auffinden von Öl und Gas. Die Methode basiert auf der Anregung elastischer Schwingungen und der anschließenden Aufzeichnung der Reaktion von Gesteinen, aus denen sich das Untersuchungsgebiet zusammensetzt. Die Anregung von Schwingungen erfolgt an Land (Dynamit oder nicht explosive Schwingungsquellen elastischer Schwingungen) oder auf See (Luftgewehre). Elastische Schwingungen breiten sich durch die Dicke der Gesteine ​​aus und werden an den Grenzen von Schichten mit unterschiedlichen Eigenschaften gebrochen und reflektiert. Reflektierte Wellen kehren an die Oberfläche zurück und werden von Geophonen an Land (typischerweise elektrodynamische Geräte, die auf der Bewegung eines in einer Spule hängenden Magneten basieren) oder Hydrophonen im Meer (basierend auf dem piezoelektrischen Effekt) aufgezeichnet. Zum Zeitpunkt des Eintreffens der Wellen kann man die Tiefen der geologischen Schichten beurteilen.

Seismische Schiffsschleppausrüstung:



Luftpistole erregt elastische Vibrationen:



Wellen gehen durch die Gesteinsmasse und werden von Hydrophonen aufgezeichnet:



Geophysikalisches Forschungsschiff „Ivan Gubkin“ am Pier der Blagoweschtschenski-Brücke in St. Petersburg:

Seismisches Signalmodell

Gesteine ​​haben unterschiedliche physikalische Eigenschaften. Für die seismische Erforschung sind die elastischen Eigenschaften am wichtigsten - die Ausbreitungsgeschwindigkeit der elastischen Schwingungen und die Dichte. Wenn zwei Schichten die gleichen oder ähnliche Eigenschaften haben, wird die Welle die Grenze zwischen ihnen nicht „bemerken“. Wenn sich die Wellengeschwindigkeiten in den Schichten unterscheiden, tritt an der Grenze der Schichten eine Reflexion auf. Je größer der Unterschied in den Eigenschaften ist, desto intensiver ist die Reflexion. Seine Intensität wird durch den Reflexionskoeffizienten (rc) bestimmt:



wobei ρ die Gesteinsdichte ist, ν die Wellengeschwindigkeit ist, 1 und 2 die obere und untere Schicht bezeichnen.

Eines der einfachsten und am häufigsten verwendeten seismischen Signalmodelle ist das Faltungsmodell, wenn eine registrierte seismische Spur als Ergebnis der Faltung einer Folge von Reflexionskoeffizienten mit einem Sondenpuls dargestellt wird:



wobei s (t) die seismische Spur ist, d.h. alles, was ein Hydrophon oder ein Geophon über eine feste Aufnahmezeit aufzeichnet, w (t) ist das Signal, das die Luftpistole erzeugt, n (t) ist zufälliges Rauschen.

Wir berechnen zum Beispiel eine synthetische seismische Vermessung. Wir werden den in der seismischen Erforschung weit verbreiteten Ricker-Impuls als Ausgangssignal verwenden.

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] 

Quelle seismischer Impuls



Wir setzen zwei Grenzen in Tiefen von 300 ms und 600 ms, und die Reflexionskoeffizienten sind Zufallszahlen.

 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]] 

Die Folge der Reflexionskoeffizienten:



Berechnen Sie die seismische Spur und zeigen Sie sie an. Da die Reflexionskoeffizienten unterschiedliche Vorzeichen haben, erhalten wir zwei abwechselnde Reflexionen auf dem seismischen Pfad.

 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 muss eine Reservierung vorgenommen werden - in der Realität wird die Tiefe der Schichten natürlich in Metern bestimmt und die Berechnung des seismischen Pfades erfolgt für den Zeitbereich. Es wäre korrekter, die Tiefen in Metern einzustellen und die Ankunftszeiten zu berechnen, wenn man die Geschwindigkeit in den Schichten kennt. In diesem Fall setze ich sofort die Ebenen auf der Zeitachse.

Wenn wir über Feldstudien sprechen, werden aufgrund solcher Beobachtungen eine große Anzahl solcher Zeitreihen (seismische Spuren) aufgezeichnet. Wenn Sie beispielsweise einen Ort erkunden, der 25 km lang und 15 km breit ist, wobei als Ergebnis der Arbeit jede Spur eine Zelle mit einer Größe von 25 x 25 Metern charakterisiert (eine solche Zelle wird als Bin bezeichnet), enthält das endgültige Datenarray 600.000 Spuren. Bei einem Zeitabtastschritt von 1 ms, einer Aufnahmezeit von 5 Sekunden beträgt die endgültige Datendatei mehr als 11 GB, und die Menge des ursprünglichen „Rohmaterials“ kann Hunderte von Gigabyte betragen.

Wie arbeite ich mit ihnen in Wolfram Mathematica ?

GeologyIO- Paket

Die Entwicklung des Pakets begann mit einer Frage an der VK-Wand der russischsprachigen Selbsthilfegruppe. Dank der Antworten der Community wurde sehr schnell eine Lösung gefunden. Infolgedessen entwickelte sich daraus eine ernsthafte Entwicklung. Der entsprechende Beitrag an der Wand der Wolfram Comunity wurde sogar von Moderatoren zur Kenntnis genommen. Derzeit unterstützt das Paket die Arbeit mit den folgenden Datentypen, die in der geologischen Industrie aktiv verwendet werden:

1. Kartendaten im ZMAP- und IRAP-Format importieren
2. Import von Messungen in LAS-Format-Wells
3. Eingabe und Ausgabe von seismischen Dateien im SEGY- Format

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 auf jedem Mathematica-Notebook aus :

 If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[ "https://wolfr.am/FiQ5oFih", FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}] ]]] 

Danach wird das Paket im Standardordner installiert, zu dem der Pfad wie folgt abgerufen werden kann:

 FileNameJoin[{$UserBasePacletsDirectory, "Repository"}] 

Als Beispiel zeigen wir die Hauptfunktionen des Pakets. Der Aufruf gilt traditionell für Pakete in Wolfram-Sprache:

 Get["GeologyIO`"] 

Das Paket wird mit Wolfram Workbench entwickelt . Auf diese Weise können Sie die Hauptfunktionalität des Pakets mit einer Dokumentation versehen, die sich im Präsentationsformat nicht von der Dokumentation von Wolfram Mathematica selbst unterscheidet, und dem Paket Testdateien für die erste Bekanntschaft bereitstellen.





Eine solche Datei ist insbesondere die Datei Marmousi.segy, ein synthetisches Modell eines vom French Petroleum Institute entwickelten geologischen Abschnitts. Mit diesem Modell testen Entwickler ihre eigenen Algorithmen zur Modellierung des Wellenfelds, zur Datenverarbeitung, zur Inversion seismischer Spuren usw. Das Marmousi-Modell selbst wird im Repository gespeichert, von wo das Paket selbst heruntergeladen wurde. Führen Sie den folgenden Code aus, um die Datei abzurufen:

 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 beinhaltet die Speicherung verschiedener Informationen über die Beobachtungen. Erstens sind dies Textkommentare. Hier werden Informationen über den Arbeitsplatz, die Namen der Unternehmen, die die Messungen durchgeführt haben, usw. eingegeben. In unserem Fall wird dieser Header durch eine Anforderung mit dem TextHeader-Schlüssel aufgerufen. Hier ist eine verkürzte Textüberschrift:

 Short[marmousi["TextHeader"]] 

"Der Marmousi-Datensatz wurde am Institut erstellt ... Mindestgeschwindigkeit von 1500 m / s und maximal 5500 m / s)"

Das tatsächliche geologische Modell kann angezeigt werden, indem die seismischen Spuren mit dem Schlüssel „Spuren“ kontaktiert werden (eines der Merkmale des Pakets ist die Groß- und Kleinschreibung der Schlüssel):

 ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"] 

Marmousi Modell:



Derzeit können Sie mit dem Paket auch Daten in Teilen aus großen Dateien herunterladen, sodass Dateien verarbeitet werden können, deren Größe mehrere zehn Gigabyte erreichen kann. Das Paket enthält auch Funktionen zum Exportieren von Daten nach .segy und zum teilweisen Umschreiben an das Ende der Datei.

Unabhängig davon ist die Funktionalität des Pakets bei der Arbeit mit der komplexen Struktur von .segy-Dateien zu beachten. Da es nicht nur den Zugriff auf Schlüssel und Indizes für einzelne Traces und Header ermöglicht, sondern auch deren Änderung beim anschließenden Schreiben in eine Datei. Viele der technischen Details der Implementierung von GeologyIO gehen über den Rahmen dieses Artikels hinaus und verdienen wahrscheinlich eine separate Beschreibung.

Die Relevanz der Spektralanalyse bei der seismischen Erforschung

Durch die Möglichkeit, seismische Materialien in Wolfram Mathematica zu importieren, können Sie die integrierte Signalverarbeitungsfunktion für experimentelle Daten verwenden. Da jede seismische Spur eine Zeitreihe ist, ist die Spektralanalyse eines der Hauptwerkzeuge für ihre Untersuchung. Zu den Voraussetzungen für die Analyse der Frequenzzusammensetzung seismischer Daten gehören beispielsweise:

1. Verschiedene Arten von Wellen zeichnen sich durch unterschiedliche Frequenzzusammensetzung aus. Dies ermöglicht es uns, nützliche Wellen auszuwählen und Interferenzwellen zu unterdrücken.
2. Gesteinseigenschaften wie Porosität und Sättigung können die Frequenzzusammensetzung beeinflussen. Auf diese Weise können Sie Gesteine ​​mit besseren Eigenschaften auswählen.
3. Schichten mit unterschiedlichen Dicken verursachen Anomalien in unterschiedlichen Frequenzbereichen.

Der dritte Absatz ist im Kontext dieses Artikels von zentraler Bedeutung. Unten finden Sie ein Code-Snippet zur Berechnung seismischer Spuren bei einer Schicht mit unterschiedlicher Dicke - ein Keilmodell. Dieses Modell wird traditionell in der seismischen Erforschung zur Analyse von Interferenzeffekten untersucht, wenn sich von vielen Schichten reflektierte Wellen überlappen.

 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 der Ausbruchsbildung:



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 die Reflexionskoeffizienten und seismischen 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 Figur gezeigte Folge von seismischen Spuren wird als seismischer Abschnitt bezeichnet. Wie Sie sehen können, kann die Interpretation auf einer intuitiven Ebene durchgeführt werden, da die Geometrie der reflektierten Wellen eindeutig dem zuvor festgelegten Modell entspricht. Wenn Sie die Spuren genauer analysieren, können Sie feststellen, dass sich die Spuren vom 1. bis zum 30. nicht unterscheiden - die Reflexion von der Oberseite der Formation und von der Sohle überlappt sich nicht. Ab Route 31 beginnen Reflexionen zu stören. Und obwohl sich im Modell die Reflexionskoeffizienten nicht horizontal ändern - seismische Spuren ändern ihre Intensität mit einer Änderung der Dicke des Reservoirs.

Betrachten Sie die Reflexionsamplitude von der oberen Grenze des Reservoirs. Ab der 60. Route beginnt die Reflexionsintensität zuzunehmen und auf der 70. Route wird sie maximal. Auf diese Weise manifestiert sich die Interferenz von Wellen vom Dach und der Sohle für Formationen, was in einigen Fällen zu signifikanten Anomalien 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, wenn das Signal eine niedrigere Frequenz hat, Interferenzen bei großen Dicken der Formation auftreten und im Fall eines Hochfrequenzsignals Interferenzen bei kleineren Dicken auftreten. Das folgende Codefragment 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 of wavelets",Black,20], LabelStyle->Directive[Black,Italic], ImageSize->Large,InterpolationOrder->2] 

Eine Reihe von Quellensignalen mit Frequenzen von 35 Hz, 55 Hz, 85 Hz



Nachdem wir die seismischen Spuren berechnet und die Amplituden der reflektierten Welle aufgetragen haben, können wir sehen, dass bei verschiedenen Frequenzen eine Anomalie bei verschiedenen Dicken der Formation 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 von der Oberkante des Keils für verschiedene Frequenzen



Die Fähigkeit, auf der Grundlage der Ergebnisse seismischer Beobachtungen Rückschlüsse auf die Dicke des Reservoirs zu ziehen, ist äußerst nützlich, da eine der Hauptaufgaben bei der Suche nach einem Ölfeld darin besteht, die vielversprechendsten Punkte für die Verlegung eines Bohrlochs zu bewerten (d. H. Die Bereiche, in denen das Reservoir eine große Dicke aufweist). Darüber hinaus kann es im geologischen Abschnitt solche Objekte geben, die aufgrund ihrer Entstehung eine starke Änderung der Dicke der Formation verursachen. Dies macht die Spektralanalyse zu einem effektiven Werkzeug, um sie zu untersuchen. Im nächsten Teil des Artikels werden wir solche gelogischen Objekte genauer betrachten.

Experimentelle Daten. Wo werden sie empfangen und worauf muss man bei ihnen achten?

Die in dem Artikel analysierten Materialien wurden auf dem Gebiet Westsibiriens erhalten. Die Region ist, wie jeder ausnahmslos weiß, die wichtigste Ölförderregion unseres Landes. Die aktive Entwicklung von U-Bahn-Geburten begann in der Region in den 60er Jahren des letzten Jahrhunderts. Die Hauptmethode zum Auffinden von Ölfeldern ist die seismische Exploration. Es ist interessant, Satellitenbilder dieses Gebiets zu betrachten. In kleinem Maßstab kann eine große Anzahl von Sümpfen und Seen festgestellt werden. Durch Vergrößern der Karte können Sie Clusterbrunnen zum Bohren von Brunnen sehen und durch Erhöhen der Karte bis zum Grenzwert können Sie auch Lichtungen von Profilen unterscheiden, mit denen seismische Beobachtungen gemacht wurden.

Satellitenbild von Yandex-Karten - das Gebiet der Stadt Noyabrsk:



Ein Netzwerk von Cluster-Standorten in einem der Felder:



Die ölhaltigen Gesteine ​​Westsibiriens kommen in einer Vielzahl von Tiefen vor - von 1 km bis 5 km. Das Hauptvolumen der ölhaltigen Gesteine ​​wird in der Jura- und Kreidezeit gebildet. Die Jurazeit ist wahrscheinlich vielen durch den gleichnamigen Film bekannt. Das Klima der Jurazeit unterschied sich erheblich vom modernen. In der Encyclopedia Britannica gibt es eine Reihe von Paläokarten, die jede gelogische Ära charakterisieren.

Die Gegenwart:



Jurazeit:



Bitte beachten Sie, dass das Gebiet Westsibiriens im Jura eine Seeküste war (Land, das von Flüssen und flachem Meer durchquert wird). Da das Klima angenehm war, kann davon ausgegangen werden, dass eine typische Landschaft dieser Zeit wie folgt aussah:

Sibirien des Jura:



In diesem Bild sind uns weniger Tiere und Vögel wichtig als das Bild des Flusses im Hintergrund. Ein Fluss ist das geologische Objekt, an dem wir früher angehalten haben. Tatsache ist, dass die Aktivität von Flüssen die Ansammlung von gut sortierten Sandsteinen ermöglicht, die dann zu einem Ölreservoir werden. Diese Stauseen können eine bizarre, komplexe Form haben (wie ein Flussbett) und sie haben eine variable Dicke - entlang der Küste ist die Dicke gering und nimmt näher an der Mitte des Kanals oder in Abschnitten der Mäander zu. Die in der Jurazeit gebildeten Flüsse befinden sich nun in einer Tiefe von etwa drei Kilometern und sind Gegenstand der Suche nach Ölreservoirs.

Experimentelle Daten. Verarbeitung und Visualisierung

Wir machen sofort einen Vorbehalt bezüglich der im Artikel gezeigten seismischen Materialien - da die für die Analyse verwendete Datenmenge erheblich ist - wird nur ein Fragment des ursprünglichen Satzes seismischer Spuren in den Text des Artikels eingefügt. Auf diese Weise kann jeder die obigen Berechnungen reproduzieren.

Bei der Arbeit mit seismischen Daten verwenden Geophysiker normalerweise spezielle Software (es gibt mehrere Branchenführer, deren Entwicklung aktiv genutzt wird, z. B. Petrel oder Paradigm), mit denen Sie verschiedene Datentypen analysieren können und die über eine praktische grafische Oberfläche verfügen. Trotz aller Bequemlichkeit haben diese Softwaretypen auch ihre Nachteile - beispielsweise nimmt die Implementierung moderner Algorithmen in stabilen Versionen viel Zeit in Anspruch, und die Möglichkeiten zur Automatisierung von Berechnungen sind normalerweise begrenzt. In einer solchen Situation wird es sehr bequem, Computermathematiksysteme und Programmiersprachen auf hoher Ebene zu verwenden, die es Ihnen ermöglichen, eine breite algorithmische Basis zu verwenden und gleichzeitig viel Routine zu übernehmen. Nach diesem Prinzip wird die Arbeit mit seismischen Daten in Wolfram Mathematica erstellt. Es ist nicht zweckmäßig, eine umfangreiche Funktion interaktiver Arbeit mit Daten zu schreiben. Es ist wichtiger, das Laden aus einem allgemein akzeptierten Format sicherzustellen, die gewünschten Algorithmen auf diese anzuwenden und sie wieder in ein externes Format hochzuladen.

Laden Sie nach dem vorgeschlagenen Schema die ursprünglichen seismischen Daten und zeigen Sie sie in Wolfram Mathematica an :

 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 auf diese Weise heruntergeladenen und importierten Daten sind die auf einer 10 x 5 Kilometer großen Website registrierten Tracks. Für den Fall, dass die Daten durch die Methode der dreidimensionalen seismischen Exploration erhalten werden (die Wellenaufzeichnung erfolgt nicht entlang separater geophysikalischer Profile, sondern gleichzeitig über das gesamte Gebiet), wird es möglich, seismische Datenwürfel zu erhalten. Dies sind dreidimensionale Objekte, deren vertikale und horizontale Abschnitte eine detaillierte Untersuchung der geologischen Umgebung ermöglichen. Im betrachteten Beispiel handelt es sich nur um dreidimensionale Daten. Wir können einige Informationen aus der Textüberschrift erhalten, zum Beispiel so

 StringPartition[seismic3DSEGY["textheader"], 80] // TableForm 

C 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 UMFRAGSNAME: ÜBERALL IN SIBIRIEN
C 7 DATEI TYP 3D SEISMISCHES VOLUMEN
C 8
C 9
C10 Z BEREICH: ERSTE 2200M LETZTE 2400M

Dieser Datensatz wird ausreichen, um die Hauptphasen der Datenanalyse zu demonstrieren. Die Spuren in der Datei werden nacheinander aufgezeichnet und jede von ihnen sieht ungefähr so ​​aus wie in der folgenden Abbildung - 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] 

Einer der seismischen Abschnitte:



Wenn Sie wissen, wie viele Spuren sich in jeder Richtung des untersuchten Bereichs befinden, können Sie ein dreidimensionales Datenarray bilden 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]] 

Dreidimensionales Bild eines Würfels seismischer Daten (vertikale Achse - Tiefe):



Für den Fall, dass geologische Objekte von Interesse intensive seismische Anomalien verursachen, können Visualisierungswerkzeuge mit Transparenz verwendet werden. "Unwichtige" Abschnitte der Aufnahme können unsichtbar gemacht werden, sodass nur Anomalien sichtbar bleiben. In Wolfram Mathematica kann dies mit Opacity [] und Raster3D [] 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] 

Bild eines seismischen Datenwürfels mit den Funktionen Opacity [] und Raster3D []



Wie im synthetischen Beispiel können auf den Scheiben des ursprünglichen Würfels einige geologische Grenzen (Schichten) mit variablem Relief unterschieden werden.

Das wichtigste Werkzeug zur Spektralanalyse ist die Fourier-Transformation. Mit ihm können Sie das Amplituden-Frequenz-Spektrum jedes Pfades oder jeder Pfadgruppe auswerten. Nach dem Übertragen von Daten in den Frequenzbereich gehen jedoch Informationen darüber verloren, zu welchen Zeiten (in welchen Tiefen gelesen) sich die Frequenz ändert. Um Signaländerungen auf der Zeitachse (tief) lokalisieren zu können, werden die Fenster-Fourier-Transformation und die Wavelet-Zerlegung verwendet. Dieser Artikel verwendet die Wavelet-Zerlegung. Die Wavelet-Analysetechnologie wurde in den 90er Jahren aktiv in der seismischen Exploration eingesetzt. Der Vorteil gegenüber der Fenster-Fourier-Transformation wird als die beste Zeitauflösung angesehen.

Mit dem folgenden Codefragment 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] 

Komponentenzerlegung



Um zu beurteilen, wie sich die Reflexionsenergie zu unterschiedlichen Ankunftszeiten der Wellen verteilt, werden Skalogramme (analog zum Spektrogramm) verwendet. In der Praxis müssen in der Regel nicht alle Komponenten analysiert werden. Wählen Sie normalerweise die Nieder-, Mittel- und Hochfrequenzkomponente.

 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. Ergebnis der WaveletScalogram [] -Funktion



Wolfram Language verwendet die Funktion ContinuousWaveletTransform [] für Wavelet-Transformationen . Die Anwendung dieser Funktion auf den gesamten Satz von Traces erfolgt mit der Funktion Table [] . Hier ist eine der Stärken von Wolfram Mathematica die Verwendung der ParallelTable [] -Parallelisierung . Im angegebenen Beispiel ist keine Parallelisierung erforderlich - das Datenvolumen ist nicht groß. Bei der Arbeit mit experimentellen Datensätzen mit Hunderttausenden von Spuren ist dies jedoch erforderlich.

 tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; 

Nach Anwendung der Funktion ContinuousWaveletTransform [] werden neue Datenfelder angezeigt, die den ausgewählten Frequenzen entsprechen. Im obigen Beispiel sind dies 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 aus Sicht eines Geologen die aussagekräftigsten aus.

Wenn Sie die Ergebnisse mit Kollegen teilen oder dem Kunden zur Verfügung stellen müssen, können Sie die SEGYExport [] -Funktion von GeologyIO 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]; 

Mit drei solchen Würfeln (Niederfrequenz-, Mittelfrequenz- und Hochfrequenzkomponenten) verwenden sie in der Regel die RGB-Mischung zur gemeinsamen Visualisierung von Daten. Jeder Komponente wird eine eigene Farbe zugewiesen - Rot, Grün, Blau. In Wolfram Mathematica kann dies mit der Funktion ColorCombine [] erfolgen .

Als Ergebnis werden Bilder erhalten, die zur Durchführung einer geologischen Interpretation verwendet werden können. Die Mäander, die an dem Abschnitt hängen, ermöglichen es Ihnen, den Paläorus zu skizzieren, der eher Reservoire ist und Ölreserven enthält. Die Suche und Analyse moderner Analoga eines solchen Flusssystems ermöglicht es uns, die vielversprechendsten Teile der Mäander zu bestimmen. Die Betten selbst zeichnen sich durch dicke Schichten aus gut sortiertem Sandstein aus und sind ein gutes Reservoir für Öl. Bereiche außerhalb der "Spitzen" -Anomalien ähneln modernen Auenablagerungen. Auensedimente werden hauptsächlich durch tonige Gesteine ​​dargestellt, und Bohrungen in diese Zonen sind unwirksam.

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 mäandrierenden Fluss verfolgen.



In einigen Fällen ermöglicht die Qualität der seismischen Daten deutlich schärfere Bilder. Dies hängt von der Methodik der Feldarbeit und der vom Geräuschreduzierungsalgorithmus verwendeten Ausrüstung ab. In solchen Fällen sind nicht nur Fragmente des Flusssystems sichtbar, sondern auch ganze ausgedehnte Paläoreken.

RGB-Mischung der drei Komponenten des seismischen Datenwürfels (horizontale Schicht). Die Tiefe beträgt ca. 2 km.



Satellitenbild der Wolga in der Region Saratow

Fazit

In Wolfram Mathematica können Sie seismische Daten analysieren und Anwendungsprobleme im Zusammenhang mit der Suche nach Mineralien lösen. Mit dem GeologyIO-Paket können Sie diesen Prozess komfortabler gestalten. Die Struktur seismischer Daten ist so aufgebaut, dass integrierte Methoden zur Beschleunigung von Berechnungen verwendet werden ( ParallelTable [] , ParallelDo [], ...) ist sehr effizient und ermöglicht die Verarbeitung großer Datenmengen. Dies wird in hohem Maße durch die Datenspeicherfunktionen des GeologyIO-Pakets erleichtert. Das Paket kann übrigens nicht nur im Bereich der angewandten seismischen Exploration eingesetzt werden. In der Georadar- und Seismologie werden fast dieselben Datentypen verwendet. Wenn Sie Vorschläge zur Verbesserung des Ergebnisses haben, welche Wolfram Mathematica-Signalanalysealgorithmen auf diese Daten anwendbar sind, oder wenn Sie kritische Kommentare haben, hinterlassen Sie Kommentare.