Wolfram Mathematica in der Geophysik

Danke dem Blogautor Anton Ekimenko fĂŒr seinen Vortrag

EinfĂŒhrung

Dieser Beitrag wurde im Anschluss an die Konferenz Wolfram Russian Technology Conference verfasst und enthÀlt eine Zusammenfassung des Vortrags, den ich gehalten habe. Die Veranstaltung fand im Juni in St. Petersburg statt. Da ich nur einen Block vom Konferenzort arbeite, konnte ich es mir nicht leisten, dieses Ereignis zu verpassen. In den Jahren 2016 und 2017 habe ich den VortrÀgen der Konferenz gelauscht, und in diesem Jahr habe ich selbst einen Vortrag gehalten. Zuerst gab es ein interessantes (wie ich finde) Thema, das wir mit Kirill Beloventwickeln, und zweitens, nach einer langen Analyse der Gesetzgebung der Russischen Föderation im Bereich der Sanktionspolitik, hat das Unternehmen, in dem ich arbeite, tatsÀchlich zwei Lizenzen erhalten Wolfram Mathematica.

Bevor ich zum Thema meines Vortrags komme, 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 Hauptattraktionen Petersburgs und ist vom Saal, in dem die Konferenz stattfand, sehr gut sichtbar.

Wolfram Mathematica in der Geophysik

Am Eingang der SPbGEU wurden die Teilnehmer von Assistenten aus den Reihen der Studenten empfangen – sie sorgten dafĂŒr, dass niemand sich verlief. Bei der Registrierung wurden kleine Souvenirs (Spielzeug – blinkender Spike, Stift, Aufkleber mit Wolfram-Symbolik) verteilt. Mittagessen und Kaffeepause waren ebenfalls im Konferenzzeitplan enthalten. Zum leckeren Kaffee und den Snacks, die ich bereits in der Gruppenwand erwĂ€hnt habe – die Köche waren großartig. Mit diesem einleitenden Teil möchte ich betonen, dass die Veranstaltung selbst, ihr Format und der Veranstaltungsort bereits positive Emotionen hervorrufen.

Der Vortrag, den ich zusammen mit Kirill Belov vorbereitet habe, trĂ€gt den Titel „Verwendung von Wolfram Mathematica zur Lösung von Aufgaben der angewandten Geophysik. Spektralanalyse seismischer Daten oder „wo flossen die alten FlĂŒsse“. Der Inhalt des Vortrags umfasst zwei Teile: erstens die Verwendung von Algorithmen, die in Wolfram Mathematica zur Analyse geophysikalischer Daten zur VerfĂŒgung stehen, und zweitens, wie geophysikalische Daten in Wolfram Mathematica integriert werden können.

Seismik

ZunĂ€chst ist es notwendig, einen kleinen Exkurs in die Geophysik zu machen. Geophysik ist die Wissenschaft, die sich mit den physikalischen Eigenschaften von Gesteinen befasst. Da Gesteine unterschiedliche Eigenschaften besitzen – elektrische, magnetische, elastische – existieren entsprechende Methoden der Geophysik: Elektroexploration, Magnetexploration, Seismik... In diesem Artikel werden wir uns jedoch nĂ€her mit der seismischen Exploration befassen. Seismik ist die Hauptmethode zur Suche nach Öl und Gas. Das Verfahren basiert auf der Erzeugung elastischer Schwingungen und der anschließenden Registrierung der RĂŒckmeldung von den Gesteinen, die das untersuchte Gebiet bilden. Die Erzeugung der Schwingungen erfolgt an Land (mit Dynamit oder nicht detonierenden vibrierenden Quellen elastischer Schwingungen) oder auf dem Meer (mit Pneumpattern). Elastische Schwingungen breiten sich durch die Schichten der Gesteine aus, refragiert und reflektiert an den Grenzen von Schichten mit unterschiedlichen Eigenschaften. Die reflektierten Wellen kehren zur OberflĂ€che zurĂŒck und werden mit Geophonen an Land (ĂŒblicherweise elektrodynamische GerĂ€te, die auf der Bewegung eines Magneten beruhen, der in einer Spule hĂ€ngt) oder Hydrofonen auf dem Meer (die auf dem Fotoeffekt basieren) registriert. Anhand der Ankunftszeiten der Wellen kann auf die Tiefen der geologischen Schichten geschlossen werden.

Das seismische Schiff schleppt AusrĂŒstung
Wolfram Mathematica in der Geophysik

Die Pneumpattern erzeugt elastische Schwingungen
Wolfram Mathematica in der Geophysik

Die Wellen durchdringen die Gesteinsschichten und werden von Hydrofonen registriert
Wolfram Mathematica in der Geophysik

Das Forschungsschiff fĂŒr geophysikalische Erkundung ‚Iwan Gubkin‘ am Kai unter der Blagoweschtschenski-BrĂŒcke in St. Petersburg
Wolfram Mathematica in der Geophysik

Modell des seismischen Signals

Gesteine weisen verschiedene physikalische Eigenschaften auf. FĂŒr die seismische Exploration sind vor allem die elastischen Eigenschaften wichtig – die Geschwindigkeit der elastischen Schwingungen und die Dichte. Wenn zwei Schichten identische oder Ă€hnliche Eigenschaften haben, wird die Welle die Grenze zwischen ihnen ‚nicht bemerken‘. Wenn sich jedoch die Wellen Geschwindigkeiten in den Schichten unterscheiden, kommt es an der GrenzflĂ€che der Schichten zu einer Reflexion. Je grĂ¶ĂŸer der Unterschied in den Eigenschaften, desto intensiver wird die Reflexion sein. Ihre IntensitĂ€t wird durch den Reflexionskoeffizienten (rc) bestimmt:

Wolfram Mathematica in der Geophysik

wobei ρ die Dichte der Gesteine und Μ die Wellen geschwindigkeit ist, 1 und 2 die obere und untere Schicht bezeichnen.

Eines der einfachsten und am hÀufigsten verwendeten Modelle seismischer Signale ist das Faltungsmodell, bei dem die aufgezeichnete seismische Spur als Ergebnis der Faltung einer Folge von Reflexionskoeffizienten mit einem anregenden Impuls dargestellt wird:

Wolfram Mathematica in der Geophysik

wo s(t) — seismische Spur, d.h. alles, was das Hydrofon oder Geophon wĂ€hrend der festen Registrierungszeit aufgezeichnet hat, w(t) — das Signal, das die Pneumatik erzeugt, n(t) — zufĂ€lliger LĂ€rm.

Berechnen wir zum Beispiel eine synthetische Seismospur. Als Ausgangssignal verwenden wir den in der seismischen Exploration weit verbreiteten Ricker-Impuls.

length=0.050; (*Signal-LĂ€nge*)
dt=0.001;(*Abtastrate des Signals*)
t=Range[-length/2,(length)/2,dt];(*Signal-Zeit*)
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["Initial wavelet",Black,20],
LabelStyle->Directive[Black,Italic],
FillingStyle->{White,Black},ImageSize->Large,InterpolationOrder->2]

UrsprĂŒnglicher seismischer Impuls
Wolfram Mathematica in der Geophysik

Wir setzen zwei Grenzen in Tiefen von 300 ms und 600 ms, und die Reflexionskoeffizienten werden zufÀllige Zahlen 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["Reflexionskoeffizienten",Black,20],
LabelStyle->Directive[Black,Italic]]

Folge der Reflexionskoeffizienten
Wolfram Mathematica in der Geophysik

Wir berechnen und visualisieren die seismische Spur. Da die Reflexionskoeffizienten unterschiedliche Vorzeichen haben, erhalten wir auch in der seismischen Spur zwei wechselnde Reflexionen.

traceExamle=ListConvolve[welle[[1;;;;1]],rcExample];
ListPlot[traceExamle,
PlotStyle->Black,Filling->0,Frame->True,Axes->False,
PlotLabel->Style["Seismische Spur",Black,20],
LabelStyle->Directive[Black,Italic]]

Simulierte Spur
Wolfram Mathematica in der Geophysik

FĂŒr dieses Beispiel muss ich darauf hinweisen, dass in der RealitĂ€t die Tiefe der Schichten natĂŒrlich in Metern bestimmt wird, wĂ€hrend die Berechnung der seismischen Spur im Zeitbereich erfolgt. Es wĂ€re richtiger, die Tiefen in Metern anzugeben und die Ankunftszeiten zu berechnen, wenn man die Geschwindigkeiten in den Schichten kennt. In diesem Fall habe ich die Schichten sofort auf der Zeitachse festgelegt.

Wenn es um Feldforschung geht, werden durch solche Beobachtungen eine enorme Menge Ă€hnlicher Zeitreihen (seismische Trassen) registriert. Zum Beispiel, bei der Untersuchung eines Gebiets von 25 km LĂ€nge und 15 km Breite, wo jede Trasse aufgrund der Arbeiten eine Zelle von 25x25 Metern charakterisiert (eine solche Zelle wird als Bin bezeichnet), wird die endgĂŒltige Datenmenge 600.000 Trassen enthalten. Bei einer Zeitdiskretisierung von 1 ms und einer Aufzeichnungszeit von 5 Sekunden wird die endgĂŒltige Datendatei ĂŒber 11 GB groß sein, wĂ€hrend das Volumen des ursprĂŒnglichen „rohen“ Materials Hunderte von Gigabyte betragen kann.

Wie man damit in Wolfram Mathematica?

Paket GeologyIO

Die Entwicklung des Pakets begann auf Stack Overflow dazu vor einem Monat gestellt, aber anscheinend ist die Zielgruppe des Projekts nicht besonders geeignet fĂŒr solche Fragen, daher musste ich selbst die Antwort aufgrund der Analyse geben. in der VK-Gruppe der russischsprachigen UnterstĂŒtzung. Dank der Antworten der Community wurde die Lösung sehr schnell gefunden. Dies hat sich zu einer ernsthaften Entwicklung entwickelt. Der entsprechende Beitrag in der Wolfram Community wurde sogar von Moderatoren hervorgehoben. Momentan unterstĂŒtzt das Paket die Arbeit mit folgenden Datentypen, die in der Geologiebranche aktiv verwendet werden:

  1. Import von kartografischen Daten im ZMAP- und IRAP-Format
  2. Import von Messungen in Bohrlöchern im LAS-Format
  3. Ein- und Ausgabe von seismischen Dateien im SEGY

Um das Paket zu installieren, befolgen Sie die Anweisungen auf der Seite zum Herunterladen des gebauten Pakets, d.h. fĂŒhren Sie den folgenden Code in einem Mathematica-Notizbuch:

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

Nach dem erfolgt die Installation des Pakets in den Standardordner, dessen Pfad folgendermaßen abgerufen werden kann:

FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]

Um die Hauptfunktionen des Pakets zu demonstrieren, wird der Aufruf traditionell fĂŒr Pakete in der Wolfram-Sprache durchgefĂŒhrt:

Get["GeologyIO`"]

Das Paket wird mit dem Einsatz von Wolfram Workbenchentwickelt. Dies ermöglicht es, die Hauptfunktionen des Pakets mit Dokumentationen zu begleiten, die in der PrĂ€sentationsform nicht von der Dokumentation selbst von Wolfram Mathematica abweichen und dem Paket Testdateien fĂŒr eine erste Bekanntschaft zur VerfĂŒgung zu stellen.

Wolfram Mathematica in der Geophysik

Wolfram Mathematica in der Geophysik

Eine solche Datei ist unter anderem die Datei „Marmousi.segy“ – ein synthetisches Modell eines geologischen Querschnitts, das vom französischen Institut fĂŒr Erdöl entwickelt wurde. Die Entwickler testen mithilfe dieses Modells ihre eigenen Algorithmen zur Modellierung des Wellenfelds, zur Datenverarbeitung, zur Inversion seismischer Spuren usw. Das Modell Marmousi selbst wird in einem Repository gespeichert, von wo aus 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
Wolfram Mathematica in der Geophysik

Das SEGY-Format sieht die Speicherung unterschiedlicher Informationen ĂŒber die Beobachtungen vor. ZunĂ€chst handelt es sich um Textkommentare. Hier werden Informationen ĂŒber den Arbeitsort, die Namen der Firmen, die die Messungen durchgefĂŒhrt haben, usw. eingetragen. In unserem Fall wird dieser Header durch eine Abfrage mit dem SchlĂŒssel TextHeader aufgerufen. Hier ist der verkĂŒrzte Textheader dargestellt:

Short[marmousi["TextHeader"]]

„Das Marmousi-Dataset wurde am Institut 
 Mindestgeschwindigkeit von 1500 m/s und eine Höchstgeschwindigkeit von 5500 m/s erzeugt“

Die geologische Modell kann angezeigt werden, indem man auf die seismischen Spuren mit dem SchlĂŒssel „traces“ zugreift (eine der Funktionen des Pakets ist die Groß-/Kleinschreibung der SchlĂŒsselunabhĂ€ngigkeit):

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

Modell Marmousi
Wolfram Mathematica in der Geophysik

Derzeit ermöglicht das Paket auch das schrittweise Laden von Daten aus großen Dateien, wodurch die Verarbeitung von Dateien mit einer GrĂ¶ĂŸe von bis zu mehreren zehn Gigabyte möglich wird. Zu den Funktionen des Pakets gehören auch Funktionen zum Export von Daten in .segy und zum teilweisen AnhĂ€ngen an das Ende der Datei.

Besonders hervorzuheben ist die FunktionalitĂ€t des Pakets im Umgang mit komplexen .segy-Dateistrukturen. Es ermöglicht nicht nur den Zugriff auf einzelne Spuren, Header und deren Änderungen ĂŒber SchlĂŒssel und Indizes, sondern auch deren folgende Speicherung in der Datei. Viele technische Details der GeologyIO-Implementierung 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 wichtigsten Werkzeuge zu ihrer Untersuchung. Zu den Faktoren, die die Analyse des Frequenzinhalts seismischer Daten beeinflussen können, gehören beispielsweise folgende:

  1. Verschiedene Wellentypen weisen unterschiedliche Frequenzzusammensetzungen auf. Dies ermöglicht es, nĂŒtzliche Wellen zu isolieren und Störwellen zu dĂ€mpfen.
  2. Eigenschaften von Gesteinen wie PorositÀt und SÀttigung können den Frequenzinhalt beeinflussen. Dies ermöglicht die Unterscheidung von Gesteinen mit besseren Eigenschaften.
  3. Schichten unterschiedlicher Dicke verursachen Anomalien in verschiedenen Frequenzbereichen.

Der dritte Punkt ist im Kontext dieses Artikels der entscheidende. Im Folgenden ist ein Codeausschnitt zur Berechnung seismischer Spuren im Fall einer Schicht mit variierender Dicke – das Keilmodell – aufgefĂŒhrt. Dieses Modell wird traditionell in der Seismik zur Analyse von Interferenzeffekten untersucht, wenn Wellen, die von vielen Schichten reflektiert werden, ĂŒbereinanderliegen.

nx=200;(* Anzahl der Gitterpunkte in X-Richtung *)
ny=200;(* Anzahl der Gitterpunkte in Y-Richtung *)
T=2;(* Gesamte Ausbreitungszeit *)
(* 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}];
(* Modellzeichnung *)
MatrixPlot[modellv,PlotLabel->Style["Modell der Schicht",Schwarz,20],
LabelStyle->Directive[Schwarz,Italic]]

Modell eines keilförmigen Schicht
Wolfram Mathematica in der Geophysik

Die Geschwindigkeit der Wellen innerhalb des Keils betrĂ€gt 4500 m/s, außerhalb des Keils 4000 m/s, und die Dichte wird konstant mit 2200 g/cmÂł angenommen. FĂŒr dieses 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],Schwarz}}},PlotStyle->Directive[Grau,DĂŒnn],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->Rot],PlotRange->All,Frame->True],270Degree]

Seismische Spuren fĂŒr das Keilmodell
Wolfram Mathematica in der Geophysik

Die Sequenz der seismischen Linien, die auf 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. Bei einer detaillierteren Analyse der Linien kann festgestellt werden, dass die Linien von 1 bis etwa 30 sich nicht unterscheiden — die Reflexion von der Oberkante und der Unterkante der Schicht ĂŒberlagern sich nicht. Ab der 31. Linie beginnen die Reflexionen zu interferieren. Und obwohl in dem Modell die Reflexionskoeffizienten horizontal unverĂ€ndert bleiben — Ă€ndern sich die seismischen Linien in ihrer IntensitĂ€t mit der Änderung der Schichtdicke.

Betrachten wir die Amplitude der Reflexion von der oberen Grenze der Schicht. Ab der 60. Linie beginnt die IntensitĂ€t der Reflexion zu steigen und erreicht bei der 70. Linie ihren Höchstwert. So zeigt sich die Interferenz der Wellen von der Ober- und Unterkante fĂŒr die Schichten, was in einigen FĂ€llen zu erheblichen Anomalien im seismischen Datensatz 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 oberen Kante des Keils
Wolfram Mathematica in der Geophysik

Es ist logisch, dass bei einem tieferfrequenten Signal die Interferenz bei grĂ¶ĂŸeren Schichtdicken auftritt, wĂ€hrend bei einem hochfrequenten Signal die Interferenz bei geringeren Dicken auftritt. 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 of wavelets",Black,20],
LabelStyle->Directive[Black,Italic],
ImageSize->Large,InterpolationOrder->2]

Set von Eingangssignalen mit Frequenzen von 35 Hz, 55 Hz, 85 Hz
Wolfram Mathematica in der Geophysik

Nachdem wir die seismischen Linien berechnet und die Diagramme der Amplituden der reflektierten Welle erstellt haben, können wir sehen, dass die Anomalie bei unterschiedlichen Frequenzen 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 Reflexion",Black,20],
LabelStyle->Directive[Black,Italic],
Frame->True]

Diagramme der Amplituden der reflektierten Welle von der oberen Kante der SchrĂ€ge fĂŒr verschiedene Frequenzen
Wolfram Mathematica in der Geophysik

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 Erkundung von Ölfeldern ist die Bewertung der vielversprechendsten Stellen fĂŒr den Bohrlochbau (d.h. die Abschnitte, in denen die Schicht eine große Dicke hat). DarĂŒber hinaus können im geologischen Profil Objekte vorkommen, die durch ihre Genese einen abrupten Wechsel der Schichtdicken verursachen. Dies macht die spektrale Analyse zu einem effektiven Werkzeug zu ihrer Untersuchung. Im nĂ€chsten Teil des Artikels werden wir solche geologischen Objekte genauer betrachten.

Experimentelle Daten. Wo wurden sie gewonnen und was ist darin zu suchen?

Die in dem Artikel analysierten Materialien stammen aus dem Gebiet der Westsibirischen Region. Diese Region ist, wie wohl jeder weiß, das Hauptölproduktionsgebiet unseres Landes. Die aktive Erschließung von LagerstĂ€tten begann in den 60er Jahren des letzten Jahrhunderts. Die Hauptmethode zur Suche nach Ölvorkommen ist die seismische Exploration. Es ist interessant, Satellitenbilder dieses Gebiets zu betrachten. Im kleinen Maßstab kann man eine riesige Anzahl von SĂŒmpfen und Seen bemerken, vergrĂ¶ĂŸert man die Karte, so sind Bohrstellen in Form von Kuppeln zu sehen, und wenn man die Karte bis ins kleinste Detail vergrĂ¶ĂŸert, kann man auch die DurchgĂ€nge der Profile erkennen, entlang derer seismische Beobachtungen durchgefĂŒhrt wurden.

Satellitenbild von Yandex Maps – Gebiet der Stadt Noyabrsk
Wolfram Mathematica in der Geophysik

Netzwerk von Kuppelbohrstellen in einem der LagerstÀtten
Wolfram Mathematica in der Geophysik

ÖlfĂŒhrende Gesteine in der Westsibirischen Region liegen in einem breiten Tiefenbereich – von 1 km bis 5 km. Der Hauptteil der ölhaltigen Gesteine wurde im Jura- und Kreidezeitalter gebildet. Viele wissen wahrscheinlich, dass die Jurazeit durch den gleichnamigen Film bekannt ist. Klima der Jurazeit war wesentlich anders als die moderne. In der Encyclopedia Britannica gibt es eine Reihe von PalĂ€okarten, die jede geologische Epoche charakterisieren.

GegenwÀrtige Zeit
Wolfram Mathematica in der Geophysik
Jurazeit
Wolfram Mathematica in der Geophysik

Beachten Sie, dass das Gebiet Westsibirien zur Jurazeit eine KĂŒstenlinie darstellte (LandflĂ€chen durch FlĂŒsse und ein flaches Meer durchzogen). Da das Klima angenehm war, kann man annehmen, dass die typische Landschaft dieser Zeit folgendermaßen aussah:

Sibirien in der Jurazeit
Wolfram Mathematica in der Geophysik

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 dem wir zuvor bereits verweilten. Der Grund dafĂŒr ist, dass die AktivitĂ€t von FlĂŒssen gut sortierte Sandsteine anhĂ€uft, die spĂ€ter zu einem Reservoir fĂŒr Erdöl werden. Diese Reservoirs können eine kurvenreiche, komplexe Form haben (wie auch das Flussbett) und sie haben eine wechselnde Dicke – an den Ufern ist die Dicke gering, wĂ€hrend sie nĂ€her zur Mitte des Flussbetts oder in MĂ€anderabschnitten zunimmt. Die in der Jurazeit gebildeten FlĂŒsse befinden sich nun auf einer Tiefe von etwa drei Kilometern und sind Ziel von Erdölsuchen.

Experimentelle Daten. Verarbeitung und Visualisierung

Wir möchten sofort anmerken, dass bezĂŒglich der seismischen Materialien, die in diesem Artikel gezeigt werden – aufgrund des erheblichen Datenvolumens, das fĂŒr die Analyse verwendet wurde – nur ein Ausschnitt des Originalsatzes seismischer Strecken in den Text des Artikels aufgenommen wurde. Dies ermöglicht allen Interessierten, die dargestellten Berechnungen nachzuvollziehen.

Wenn Geophysiker mit seismischen Daten arbeiten, nutzen sie in der Regel spezialisierte Software (einige MarktfĂŒhrer, deren Entwicklungen aktiv genutzt werden, sind zum Beispiel Petrel oder Paradigm), die es ermöglicht, verschiedene Datentypen zu analysieren und eine benutzerfreundliche grafische OberflĂ€che bietet. Trotz aller Bequemlichkeit haben solche Softwarelösungen auch ihre Nachteile – so dauert die Implementierung moderner Algorithmen in stabile Versionen oft viel Zeit, und die Möglichkeiten zur Automatisierung von Berechnungen sind in der Regel begrenzt. In einer solchen Situation ist die Verwendung von Systemen der Computer-Mathematik und Hochsprache-Programmiersprachen sehr praktisch, die eine breite algorithmische Basis nutzen und zudem viel Routinearbeit ĂŒbernehmen. Nach diesem Prinzip funktioniert auch die Arbeit mit seismischen Daten in Wolfram Mathematica. Es ist nicht sinnvoll, eine umfangreiche FunktionalitĂ€t fĂŒr die interaktive Arbeit mit Daten zu schreiben – wichtiger ist es, den Import aus einem gebrĂ€uchlichen Format zu gewĂ€hrleisten, die gewĂŒnschten Algorithmen anzuwenden und sie anschließend in ein externes Format zurĂŒckzuladen.

GemĂ€ĂŸ dem vorgeschlagenen Schema laden wir die ursprĂŒnglichen seismischen Daten und stellen sie dar in Wolfram Mathematica:

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 Trassen, die auf einem Gebiet von 10 mal 5 Kilometern registriert wurden. Wenn die Daten nach der Methode der dreidimensionalen Seismik gewonnen 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 erzeugen. Dies sind dreidimensionale Objekte, deren vertikale und horizontale Schnitte eine detaillierte Untersuchung des geologischen Umfelds ermöglichen. Im betrachteten Beispiel haben wir es mit genau dreidimensionalen Daten zu tun. Einige Informationen können wir aus dem Textheader entnehmen, zum Beispiel so

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

C 1 DIES IST DIE DEMO-DATEI FÜR DEN TEST DES GEOLOGYIO-PAKETS
C 2
C 3
C 4
C 5 DATUM BENUTZERNAME: WOLFRAM BENUTZER
C 6 BEFRAGUNGSNAME: IRGENDWO IN SIBIRIEN
C 7 DATEITYP 3D SEISMISCHER 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 jeder von ihnen sieht ungefÀhr so aus wie im folgenden Bild - es handelt sich um 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 Traces des seismischen Profils
Wolfram Mathematica in der Geophysik

Wenn man weiß, wie viele Traces in jeder Richtung des untersuchten Gebiets vorhanden sind, kann man ein dreidimensionales Datenarray erstellen und es mit der Funktion Image3D[] darstellen.

traces=seismic3DSEGY["traces"];
startIL=1050;EndIL=2000;stepIL=2; (*X-Koordinate des Anfangs und Endes der Aufnahme und Schritt der Traces*)
startXL=1165;EndXL=1615;stepXL=2; (*Y-Koordinate des Anfangs und Endes der Aufnahme und 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)
Wolfram Mathematica in der Geophysik

Falls geologische Objekte von Interesse intensive seismische Anomalien erzeugen, können Visualisierungstools mit Transparenz verwendet werden. Unwichtige Teile der Aufzeichnung können unsichtbar gemacht werden, wÀhrend nur die Anomalien sichtbar bleiben. In Wolfram Mathematica kann dies erreicht werden mit Opacity[] und Raster3D[].

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 des WĂŒrfels seismischer Daten mit den Funktionen Opacity[] und Raster3D[]
Wolfram Mathematica in der Geophysik

Wie im synthetischen Beispiel können an den Schnitten des ursprĂŒnglichen WĂŒrfels einige geologische Grenzen (Schichten) mit variierendem Relief hervorgehoben werden.

Das Hauptwerkzeug der Spektralanalyse ist die Fourier-Transformation. Damit lĂ€sst sich das Amplitudenfrequenzspektrum jeder Spur oder Gruppe von Spuren bewerten. Allerdings geht beim Übergang der Daten in den Frequenzbereich die Information darĂŒber verloren, zu welchen Zeiten (d.h. in welchen Tiefen) sich die Frequenz Ă€ndert. Um die Änderungen des Signals auf der Zeit- (Tiefen-)Achse lokalisieren zu können, werden die Fenster-Fourier-Transformation und die Wavelet-Zerlegung verwendet. In diesem Artikel wird das Wavelet-Zerlegung verwendet. Die Technologie der Wavelet-Analyse begann in den 90er Jahren aktiv in der seismischen Erkundung Anwendung zu finden. Ein Vorteil gegenĂŒber der Fenster-Fourier-Transformation ist die bessere zeitliche Auflösung.

Mit folgendem Codefragment kann die Zerlegung in einzelne Komponenten einer der seismischen Spuren 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 der Spur in Komponenten
Wolfram Mathematica in der Geophysik

Um zu bewerten, wie die Energie der Reflexionen zu verschiedenen Ankunftszeiten der Wellen verteilt ist, werden Skalenogramme (Àhnlich wie Spektrogramme) verwendet. In der Praxis ist es in der Regel nicht notwendig, alle Komponenten zu analysieren. Normalerweise werden die niedrigen, mittleren und hohen 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->{"Zeit","Frequenz(Hz)"},LabelStyle->Directive[Black,Bold,14],
ColorFunction->"RustTones",ImageSize->Large]

Skalogramm. Ergebnis der Funktion WaveletScalogram[]
Wolfram Mathematica in der Geophysik

In der Wolfram-Sprache wird fĂŒr die Wavelet-Transformation die Funktion ContinuousWaveletTransform[]. Die Anwendung dieser Funktion auf den gesamten Satz von Spuren erfolgt mit Hilfe der Funktion Table[]. Hier sollte man einen der großen Vorteile von Wolfram Mathematica, die Möglichkeit der Parallelisierung, erwĂ€hnen. ParallelTable[]. Im gegebenen Beispiel ist keine Parallelisierung erforderlich – das Datenvolumen ist nicht groß, aber beim Arbeiten mit experimentellen DatensĂ€tzen, die Hunderte von Tausenden 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 der Anwendung der Funktion ContinuousWaveletTransform[] Es erscheinen neue DatensĂ€tze, die den gewĂ€hlten Frequenzen entsprechen. Im obigen Beispiel sind dies die Frequenzen: 38 Hz, 33 Hz, 27 Hz. Die Auswahl der Frequenzen erfolgt hĂ€ufig auf Basis von Tests, wobei effektive Karten fĂŒr verschiedene Frequenzkombinationen erstellt und die informativste aus geologischer Sicht ausgewĂ€hlt wird.

Wenn die Ergebnisse mit Kollegen geteilt oder einem Kunden zur VerfĂŒgung gestellt werden mĂŒssen, kann die Funktion SEGYExport[] des GeologyIO-Pakets verwendet werden.

outputdata=seismic3DSEGY;
outputdata["traces",1;;-1]=tracesCWD[[All,3]];
outputdata["textheader"]="Wavelet-Dekompositions-Ergebnis";
outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]];
SEGYExport["D:result.segy",outputdata];

Mit drei solchen Kuben (niedrigfrequente, mittel- und hochfrequente Komponenten) wird normalerweise eine RGB-Mischung zur gemeinsamen Visualisierung der Daten verwendet. Jeder Komponente wird eine Farbe zugewiesen – rot, grĂŒn, blau. In Wolfram Mathematica kann dies mit der Funktion gemacht werden. ColorCombine[].

Das Ergebnis sind Bilder, anhand derer geologische Interpretationen vorgenommen werden können. Die MĂ€ander, die im Schnitt sichtbar werden, ermöglichen es, alte Flussbetten zu konturieren, die mit höherer Wahrscheinlichkeit Reservoirs enthalten und Öl fĂŒhren können. Die Suche und Analyse moderner Analogien solcher Flusssysteme ermöglicht es, die vielversprechendsten Teile der MĂ€ander zu bestimmen. Eigentliche Flussbetten sind durch mĂ€chtige Schichten gut sortierten Sandsteins charakterisiert und sind gute Reservoirs fĂŒr Öl. Bereiche außerhalb der "Schnur-Anomalien" sind Ă€hnlich den modernen Auenablagerungen. Auenablagerungen bestehen hauptsĂ€chlich aus tonigen Gesteinen, und Bohrungen in diese Zonen werden ineffektiv sein.

RGB-Schnitt des Datenkubus. In der Mitte (etwas links von der Mitte) kann man den mÀandrierenden Fluss verfolgen.
Wolfram Mathematica in der Geophysik
RGB-Schnitt des Datenkubus. In der linken BildhÀlfte kann man den mÀandrierenden Fluss verfolgen.
Wolfram Mathematica in der Geophysik

In einigen FĂ€llen ermöglicht die QualitĂ€t der seismischen Daten wesentlich klarere Bilder. Dies hĂ€ngt von der Methode der Feldarbeiten, der eingesetzten AusrĂŒstung und den angewandten Algorithmen zur GerĂ€uschunterdrĂŒckung ab. In solchen FĂ€llen sind nicht nur Fragmente von Flusssystemen sichtbar, sondern auch ganze ausgedehnte alte FlĂŒsse.

RGB-Mischung der drei Komponenten des seismischen Datenkubus (horizontaler Schnitt). Die Tiefe betrÀgt etwa 2 km.
Wolfram Mathematica in der Geophysik
Satellitenbild des Flusses Wolga in der NĂ€he von Saratow.
Wolfram Mathematica in der Geophysik

Fazit

In Wolfram Mathematica können seismische Daten analysiert und praktische Aufgaben im Zusammenhang mit der Suche nach BodenschĂ€tzen gelöst werden. Das Paket GeologyIO macht diesen Prozess einfacher. Die Struktur der seismischen Daten ist so, dass die Verwendung integrierter Methoden zur Beschleunigung von Berechnungen (ParallelTable[], ParallelDo[],
) sehr effizient ist und große Datenmengen verarbeitet werden können. Dies wird erheblich durch die Speichereigenschaften des GeologyIO-Pakets unterstĂŒtzt. Übrigens kann das Paket nicht nur im Bereich der angewandten seismischen Untersuchung eingesetzt werden. Praktisch dieselben Datentypen werden in der Georadaruntersuchung und Seismologie verwendet. Wenn Sie VorschlĂ€ge zur Verbesserung der Ergebnisse haben, welche Signalverarbeitungsalgorithmen aus dem Sortiment von Wolfram Mathematica auf solche Daten anwendbar sind oder wenn Sie kritische Anmerkungen haben – hinterlassen Sie bitte Ihre Kommentare.

Quelle: habr.com

60GB SSD 8Gb DDR4