Ringraziamo l'autore del blog per la sua presentazione
Introduzione
Questo articolo è scritto a seguito della conferenza e contiene un riassunto della relazione che ho presentato. L'evento si è tenuto a giugno nella città di San Pietroburgo. Considerando che lavoro a pochi passi dal luogo della conferenza, non ho potuto mancare a questo evento. Nel 2016 e 2017 ho seguito le presentazioni della conferenza, e quest'anno ho avuto l'onore di presentare. In primo luogo, è emersa un'idea interessante (a mio avviso) che stiamo sviluppando con , e in secondo luogo, dopo un lungo studio della legislazione russa in materia di politiche sanzionatorie, presso l'azienda dove lavoro sono state ottenute ben due licenze .
Prima di passare al tema del mio intervento, vorrei sottolineare la buona organizzazione dell'evento. Sulla pagina di presentazione della conferenza è utilizzata un'immagine della Cattedrale di Kazan. La cattedrale è una delle principali attrazioni di San Pietroburgo ed è ben visibile dalla sala in cui si è svolta la conferenza.
All'ingresso del SPbGEU i partecipanti sono stati accolti da assistenti tra gli studenti, che hanno impedito di perdersi. Durante la registrazione sono stati distribuiti piccoli souvenir (un giocattolo – un mitico spiky lampeggiante, penne, adesivi con il simbolo Wolfram). Il pranzo e il coffee break erano anche inclusi nell'agenda della conferenza. Ho già menzionato il delizioso caffè e i pasticcini sulla bacheca del gruppo – i cuochi sono stati fantastici. Con questa introduzione vorrei sottolineare che l'evento stesso, il suo formato e il luogo di svolgimento già portano emozioni positive.
La relazione che ho preparato insieme a Kirill Belov si intitola «Utilizzare Wolfram Mathematica per risolvere problemi di geofisica applicata. Analisi spettrale dei dati sismici o «dove correvano i fiumi antichi». Il contenuto della relazione copre due parti: in primo luogo, l'uso di algoritmi presenti in per l'analisi dei dati geofisici, e in secondo luogo, come inserire i dati geofisici in Wolfram Mathematica.
Sismica
Per iniziare, è necessario fare una breve introduzione alla geofisica. La geofisica è la scienza che studia le proprietà fisiche delle rocce. Poiché le rocce presentano diverse caratteristiche: elettriche, magnetiche, elastiche, esistono metodi corrispondenti della geofisica: prospezione elettrica, prospezione magnetica, prospezione sismica… In questo articolo, ci concentreremo solo sulla prospezione sismica. La prospezione sismica è il metodo principale per la ricerca di petrolio e gas. Il metodo si basa sull'eccitazione di oscillazioni elastiche e sulla successiva registrazione della risposta delle rocce che compongono l'area studiata. L'eccitazione delle oscillazioni avviene sulla terra (utilizzando dinamite o fonti vibranti non esplosive) o in mare (con cannoni pneumatici). Le oscillazioni elastiche si propagano attraverso gli strati di roccia, rifrangendosi e riflettendosi ai confini degli strati con proprietà diverse. Le onde riflesse tornano in superficie e vengono registrate da geofoni sulla terra (di solito sono strumenti elettrodinamici basati sul movimento di un magnete sospeso in una bobina) o da idrofoni in mare (basati sull'effetto piezoelettrico). Da quando arrivano le onde, è possibile determinare le profondità degli strati geologici.
La nave sismica trascina l'attrezzatura
La pistola pneumatica genera oscillazioni elastiche
Le onde attraversano lo spessore delle rocce e vengono registrate dai idrofoni
La nave di ricerca geofisica "Ivan Gubkin" al molo presso il ponte di Blagoveshchensky a San Pietroburgo
Modello del segnale sismico
Le rocce hanno diverse proprietà fisiche. Per la sismica, le proprietà elastiche sono particolarmente importanti: la velocità di propagazione delle onde elastiche e la densità. Se due strati hanno proprietà identiche o simili, l'onda non percepirà il confine tra di essi. Tuttavia, se le velocità delle onde negli strati sono diverse, si verificherà un riflesso al confine degli strati. Maggiore è la differenza nelle proprietà, più intenso sarà il riflesso. La sua intensità sarà determinata dal coefficiente di riflessione (rc):
dove ρ è la densità delle rocce, ν è la velocità delle onde, 1 e 2 indicano i livelli superiore e inferiore.
Uno dei modelli di segnale sismico più semplici e comunemente utilizzati è il modello convoluzionale, dove il tracciato sismico registrato viene rappresentato come il risultato della convoluzione di una sequenza di coefficienti di riflessione con un impulso di sondaggio:
dove s(t) è il tracciato sismico, cioè tutto ciò che è stato registrato dal idrofoni o geofoni durante un tempo di registrazione fisso, w(t) è il segnale generato da una pistola pneumatica, n(t) è il rumore casuale.
Calcoliamo, per esempio, un tracciato sismico sintetico. Come segnale iniziale utilizzeremo l'impulso di Ricker, ampiamente usato nella sismica.
length=0.050; (*Lunghezza del segnale*)
dt=0.001;(*Frequenza di campionamento del segnale*)
t=Range[-length/2,(length)/2,dt];(*Tempo del segnale*)
f=35;(*Frequenza centrale*)
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["Iniziale wavelet",Black,20],
LabelStyle->Directive[Black,Italic],
FillingStyle->{White,Black},ImageSize->Large,InterpolationOrder->2]Impulso sismico iniziale
Definiamo due confini a profondità di 300 ms e 600 ms, e i coefficienti di riflessione saranno numeri casuali.
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["Coefficienti di riflessione",Black,20],
LabelStyle->Directive[Black,Italic]]Sequenza dei coefficienti di riflessione
Calcoleremo e mostreremo il tracciato sismico. Poiché i coefficienti di riflessione hanno segni diversi, nel tracciato sismico otteniamo quindi due riflessioni alternate.
traceExamle=ListConvolve[wavelet[[1;;;;1]],rcExample];
ListPlot[traceExamle,
PlotStyle->Black,Filling->0,Frame->True,Axes->False,
PlotLabel->Style["Tracciato sismico",Black,20],
LabelStyle->Directive[Black,Italic]]Tracciato simulato
Per questo esempio è necessario fare una precisazione: in realtà la profondità degli strati è definita in metri, mentre il calcolo del tracciato sismico avviene in un'area temporale. Sarebbe più corretto impostare le profondità in metri e calcolare i tempi di arrivo conoscendo le velocità negli strati. In questo caso, ho subito impostato gli strati sull'asse temporale.
Если говорить о полевых исследованиях, то в результате таких наблюдений регистрируется огромное количество подобных временных рядов (сейсмических трасс). Например, при исследовании участка длиной 25 км и шириной 15 км, где в результате работ каждая трасса характеризует ячейку размером 25х25 метров (такая ячейка называется бин), финальный массив данных будет содержать 600000 трасс. При шаге дискретизации по времени равном 1 мс, времени записи 5 секунд окончательный файл данных составит более 11 Гб, а объём исходного «сырого» материала может составить сотни гигабайт.
Как работать с ними в ?
Il pacchetto
Началом разработки пакета стал на стене VK группы русскоязычной поддержки. Благодаря ответам сообщества решение было найдено очень быстро. И в результате переросло в серьёзную разработку. Соответствующий даже был отмечен модераторами. На текущий момент пакет поддерживает работу со следующими типами данных, которые активно используются в геологической отрасли:
- импорт картографических данных формата ZMAP и IRAP
- импорт измерений в скважинах формата LAS
- input e output di file sismici in formato
Per installare il pacchetto, seguire le istruzioni sulla pagina di download del pacchetto assemblato, ovvero eseguire il seguente codice in qualsiasi :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
"https://wolfr.am/FiQ5oFih",
FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]Dopo di che, il pacchetto verrà installato nella cartella predefinita, il cui percorso può essere ottenuto nel seguente modo:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]Per esempio, dimostriamo le funzionalità principali del pacchetto. La chiamata è eseguita tradizionalmente per i pacchetti in Wolfram Language:
Get["GeologyIO`"]Il pacchetto è sviluppato utilizzando . Ciò consente di accompagnare la funzionalità principale del pacchetto con documentazione che non differisce nel formato di presentazione dalla documentazione stessa di Wolfram Mathematica e fornire al pacchetto file di test per un primo approccio.
Il file "Marmousi.segy" è un modello sintetico di sezione geologica sviluppato dall'Institut Français du Pétrole. Utilizzando questo modello, gli sviluppatori testano i propri algoritmi di simulazione del campo ondoso, elaborazione dei dati, inversione delle tracce sismiche e altro. Il modello Marmousi stesso è conservato in un repository, dal quale è stato scaricato il pacchetto. Per ottenere il file, eseguiamo il seguente codice:
If[Not[FileExistsQ["Marmousi.segy"]],
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]Il risultato dell'importazione è un oggetto SEGYData.
Il formato SEGY prevede la conservazione di diverse informazioni sulle osservazioni. Prima di tutto, ci sono commenti testuali. Qui viene annotata l'informazione riguardante il luogo di lavoro, i nomi delle aziende che hanno effettuato le misurazioni e così via. Nel nostro caso, questo intestazione è richiamata tramite una query con la chiave TextHeader. Qui è riportato un'intestazione testuale abbreviata:
Short[marmousi["TextHeader"]]«Il set di dati Marmousi è stato generato presso l'Institut ... velocità minima di 1500 m/s e una massima di 5500 m/s)»
È possibile visualizzare il modello geologico richiedendo le tracce sismiche utilizzando la chiave «traces» (una delle funzionalità del pacchetto è l'indipendenza delle chiavi dalla maiuscole/minuscole):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Modello Marmousi
Attualmente, il pacchetto consente anche il caricamento dei dati in parti da file di grandi dimensioni, il che rende possibile l'elaborazione di file che possono raggiungere decine di gigabyte. Inoltre, le funzionalità del pacchetto includono funzioni per l'esportazione dei dati in formato .segy e per l'aggiunta parziale alla fine del file.
Vale la pena notare la funzionalità del pacchetto nell'elaborazione di strutture complesse di file .segy. Infatti, consente non solo di accedere alle singole tracce e intestazioni tramite chiavi e indici, ma anche di modificarle e registrarle nel file. Molti dettagli tecnici sull'implementazione di GeologyIO vanno oltre il contenuto di questo articolo e meritano probabilmente una descrizione separata.
Rilevanza dell'analisi spettrale nella sismica
La possibilità di importare materiali sismici in Wolfram Mathematica consente di utilizzare le funzionalità integrate per l'elaborazione dei segnali su dati sperimentali. Poiché ogni traccia sismica è una serie temporale, uno degli strumenti principali per il loro studio è l'analisi spettrale. Tra i presupposti per l'analisi del contenuto in frequenza dei dati sismici si possono considerare, ad esempio, i seguenti:
- I diversi tipi di onde presentano composizioni in frequenza differenti. Questo consente di identificare onde utili e di sopprimere onde di disturbo.
- Proprietà delle rocce come la porosità e la saturazione possono influenzare il contenuto in frequenza. Questo permette di distinguere le rocce con migliori caratteristiche.
- Strati di spessore diverso causano anomalie in diverse gamme di frequenza.
Il terzo punto è fondamentale nel contesto di questo articolo. Di seguito un frammento di codice per calcolare le tracce sismiche nel caso di uno strato con spessore variabile — il modello a cuneo. Questo modello è tradizionalmente studiato nella sismica per analizzare gli effetti di interferenza, quando le onde riflesse da più strati si sovrappongono.
nx=200;(* Numero di punti della griglia nella direzione X*)
ny=200;(* Numero di punti della griglia nella direzione Y*)
T=2;(*Tempo totale di propagazione*)
(*Velocità e densità*)
modellv=Table[4000,{i,1,ny},{j,1,nx}];(* Velocità delle onde P in m/s*)
rho=Table[2200,{i,1,ny},{j,1,nx}];(* Densità in g/cm³, densità costante utilizzata*)
Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}];
Table[modellv[[;;70,i]]=4500;,{i,1,200}];
(*Tracciamento del modello*)
MatrixPlot[modellv,PlotLabel->Style["Modello di strato",Black,20],
LabelStyle->Directive[Black,Italic]]Modello di strato a cuneo
La velocità delle onde all'interno del cuneo è di 4500 m/s, al di fuori del cuneo 4000 m/s, e la densità è assunta costante a 2200 g/cm³. Per tale modello calcoleremo i coefficienti di riflessione e le tracce sismiche.
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]Traçes sismiche per il modello a cuneo
La sequenza delle tracce sismiche mostrata in questa figura è chiamata sezione sismica. Come si può notare, la sua interpretazione può avvenire anche a un livello intuitivo, poiché la geometria delle onde riflesse corrisponde chiaramente al modello che è stato definito in precedenza. Se si analizzano le tracce in modo più dettagliato, si può osservare che le tracce dalla 1ª alla 30ª circa non differiscono — le riflessioni dalla partesuperiore e dalla base non si sovrappongono. A partire dalla 31ª traccia, le riflessioni iniziano a interferire. E, anche se nel modello i coefficienti di riflessione non variano orizzontalmente, le tracce sismiche cambiano la loro intensità al variare dello spessore del strato.
Consideriamo l'ampiezza della riflessione dalla superficie superiore del strato. A partire dalla 60ª traccia, l'intensità della riflessione inizia ad aumentare e alla 70ª traccia raggiunge il massimo. Questo manifesta l'interferenza delle onde dalla copertura e dalla base per i strati, portando in alcuni casi a notevoli anomalie nella registrazione sismica.
ListLinePlot[GaussianFilter[Abs[traces[[All,46]]],3][[;;;;2]],
InterpolationOrder->2,Frame->True,PlotStyle->Black,
PlotLabel->Style["Ampiezza della riflessione", Black, 20],
LabelStyle->Directive[Black, Italic],
PlotRange->All]Grafico dell'ampiezza dell'onda riflessa dall'orlo superiore del cuneo
È logico che quando il segnale ha una frequenza più bassa, l'interferenza inizia a manifestarsi a spessori maggiori del strato, mentre nel caso di un segnale ad alta frequenza l'interferenza si verifica a spessori minori. Il seguente frammento di codice crea un segnale con frequenze di 35 Hz, 55 Hz e 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["Serie di wavelet", Black, 20],
LabelStyle->Directive[Black, Italic],
ImageSize->Large,InterpolationOrder->2]Set di segnali originali con frequenze di 35 Hz, 55 Hz, 85 Hz
Eseguendo il calcolo dei tracciati sismici e costruendo grafici delle ampiezze dell'onda riflessa, possiamo vedere che per diverse frequenze l'anomalia si osserva a spessori diversi del strato.
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]Grafici delle ampiezze delle onde riflesse dal bordo superiore del cunetto per diverse frequenze
La possibilità di trarre conclusioni sulla profondità del giacimento dai risultati delle osservazioni sismiche è estremamente utile, poiché uno degli obiettivi principali nella ricerca di giacimenti petroliferi è valutare i punti più promettenti per la perforazione del pozzo (cioè le aree dove il giacimento ha una maggiore profondità). Inoltre, nel profilo geologico potrebbero trovarsi oggetti che, attraverso la loro genesi, determinano un cambiamento brusco della profondità del giacimento. Questo rende l'analisi spettrale uno strumento efficace per il loro studio. Nella parte successiva dell'articolo esamineremo più da vicino tali oggetti geologici.
Dati sperimentali. Dove sono stati ottenuti e cosa cercare in essi?
I materiali analizzati nell'articolo provengono dalla regione della Siberia occidentale. Questa area, come probabilmente tutti sanno, è il principale regione estrattiva di petrolio del nostro paese. Lo sviluppo attivo dei giacimenti è iniziato negli anni '60 del secolo scorso. Il metodo principale per la ricerca dei giacimenti petroliferi è la sismica. È interessante osservare le immagini satellitari di questa area. A una bassa scala, si può notare un enorme numero di paludi e laghi; aumentando la mappa, si possono vedere le piattaforme di trivellazione e ingrandendo ulteriormente, si possono distinguere anche i profili dei tracciati lungo i quali sono state effettuate le osservazioni sismiche.
Immagine satellitare delle mappe Yandex - zona della città di Novy Urengoy
Rete di piattaforme in uno dei giacimenti
Le rocce ricche di petrolio della Siberia occidentale si trovano a una vasta gamma di profondità, da 1 km a 5 km. La maggior parte delle rocce contenenti petrolio si è formata nel periodo giurassico e cretaceo. Il periodo giurassico è probabilmente ben noto a molti grazie al film omonimo. era significativamente diverso da quello attuale. Nell'Enciclopedia Britannica c'è una serie di paleocarte che caratterizzano ogni epoca geologica.
Tempi moderni
Periodo Giurassico
Si noti che durante il periodo giurassico il territorio della Siberia occidentale rappresentava una costa marina (terra attraversata da fiumi e un mare poco profondo). Poiché il clima era mite, si può ipotizzare che il paesaggio tipico di quell'epoca si presentasse come segue:
Siberia del periodo giurassico
In questa immagine, ciò che ci interessa non sono tanto gli animali e gli uccelli, quanto l'immagine del fiume sullo sfondo. Il fiume è quell'oggetto geologico su cui ci siamo soffermati in precedenza. Infatti, l'attività dei fiumi permette l'accumulo di arenarie ben ordinate, che diventeranno poi serbatoi per il petrolio. Questi serbatoi possono avere forme bizzarre e complesse (proprio come il letto del fiume) e presentano uno spessore variabile: vicino alle rive lo spessore è ridotto, mentre al centro del letto o in tratti di meandro aumenta. Pertanto, i fiumi formati nell'era giurassica ora si trovano a una profondità di circa tre chilometri e sono oggetto di ricerca per i serbatoi di petrolio.
Dati sperimentali. Elaborazione e visualizzazione
Facciamo subito una precisazione riguardo ai materiali sismici mostrati nell'articolo: poiché il volume di dati utilizzati per l'analisi è significativo, nel testo dell'articolo è stato inserito solo un frammento del set originale di tracciati sismici. Questo consentirà a chiunque desideri riprodurre i calcoli forniti.
Lavorando con i dati sismici, i geofisici utilizzano generalmente software specializzato (ci sono diversi leader di settore, come Petrel o Paradigm, le cui soluzioni sono ampiamente adottate) che consente di analizzare diversi tipi di dati e offre un'interfaccia grafica intuitiva. Nonostante la praticità, tali software presentano anche dei limiti: per esempio, l'integrazione di algoritmi moderni nelle versioni stabili richiede molto tempo e le possibilità di automazione dei calcoli sono generalmente limitate. In questa situazione, risulta molto utile l'uso di sistemi di matematica computazionale e linguaggi di programmazione di alto livello, che offrono una vasta base algoritmica e si occupano di molte operazioni ripetitive. Questa è la filosofia su cui si basa l'elaborazione dei dati sismici in Wolfram Mathematica. Non è produttivo scrivere funzionalità complesse per l'interazione con i dati; è più importante garantire il caricamento da formati standard, applicare gli algoritmi desiderati e quindi esportare di nuovo nel formato esterno.
Seguendo lo schema proposto, caricheremo i dati sismici originali e li visualizzeremo in :
Get["GeologyIO`"]
seismic3DZipPath = "seismic3D.zip";
seismic3DSEGYPath = "seismic3D.sgy";
If[FileExistsQ[seismic3DZipPath], DeleteFile[seismic3DZipPath]];
If[FileExistsQ[seismic3DSEGYPath], DeleteFile[seismic3DSEGYPath]];
URLDownload["https://wolfr.am/FiQIuZuH", seismic3DZipPath];
ExtractArchive[seismic3DZipPath];
seismic3DSEGY = SEGYImport[seismic3DSEGYPath]I dati caricati e importati in questo modo sono tracciati registrati su un'area di 10 per 5 chilometri. Nel caso in cui i dati siano stati ottenuti tramite la metodologia della sismica 3D (la registrazione delle onde avviene non lungo singoli profili geofisici, ma sull'intera area contemporaneamente), è possibile ottenere cubi di dati sismici. Questi sono oggetti tridimensionali, le sezioni verticali e orizzontali delle quali consentono di studiare in dettaglio l'ambiente geologico. Nel caso considerato, trattiamo proprio dati tridimensionali. Alcune informazioni possiamo ottenerle dall'intestazione testuale, ad esempio in questo modo
StringPartition[seismic3DSEGY["textheader"], 80] // TableFormC 1 QUESTO È UN FILE DIMOSTRATIVO PER IL TEST DEL PACCHETTO GEOLOGYIO
C 2
C 3
C 4
C 5 DATA NOME UTENTE: WOLFRAM USER
C 6 NOME DEL SAGGIO: DA QUALCHE PARTE IN SIBERIA
C 7 TIPO DI FILE VOLUME SISMICO 3D
C 8
C 9
C10 INTERVALLO Z: PRIMO 2200M ULTIMO 2400M
Questo set di dati sarà sufficiente per dimostrare le fasi principali dell'analisi dei dati. Le tracce nel file sono registrate in sequenza e ognuna di esse appare circa come nell'immagine seguente: questa rappresenta la distribuzione delle ampiezze delle onde riflettute lungo l'asse verticale (asse delle profondità).
ListLinePlot[seismic3DSEGY["traces"][[100]], InterpolationOrder -> 2,
PlotStyle -> Black, PlotLabel -> Style["Traccia sismica", Black, 20],
LabelStyle -> Directive[Black, Italic], PlotRange -> All,
Frame -> True, ImageSize -> 1200, AspectRatio -> 1/5]Una delle tracce del profilo sismico
Conoscendo il numero di tracce situate in ogni direzione del sito studiato, è possibile formare un array tridimensionale di dati e visualizzarlo utilizzando la funzione Image3D[]
traces=seismic3DSEGY["traces"];
startIL=1050; EndIL=2000; stepIL=2; (*coordinate X di inizio e fine rilevamento e passo delle tracce*)
startXL=1165; EndXL=1615; stepXL=2; (*coordinate Y di inizio e fine rilevamento e passo delle tracce*)
numIL=(EndIL-startIL)/stepIL+1; (*numero di tracce sull'asse X*)
numXL=(EndXL-startXL)/stepIL+1; (*numero di tracce sull'asse Y*)
Image3D[ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}], ViewPoint->{-1, 0, 0}, Background->RGBColor[0,0,0]]Immagine tridimensionale di un cubo di dati sismici. (L'asse verticale rappresenta la profondità)
Nel caso in cui gli oggetti geologici di interesse generino anomalie sismiche intense, è possibile utilizzare strumenti di visualizzazione con trasparenza. Le aree di registrazione considerate "non significative" possono essere rese invisibili, lasciando visibili solo le anomalie. In Wolfram Mathematica, questo può essere fatto con e .
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]Immagine di un cubo di dati sismici utilizzando le funzioni Opacity[] e Raster3D[]
Come nel caso di esempio sintetico, in alcune sezioni del cubo originale è possibile evidenziare alcuni confini geologici (strati) con un rilievo variabile.
Lo strumento principale per l'analisi spettrale è la trasformata di Fourier. Con essa è possibile valutare lo spettro ampiezza-frequenza di ciascuna traccia o di un gruppo di tracce. Tuttavia, una volta traslati i dati nell'area delle frequenze, si perde l'informazione riguardo a quali tempi (cioè a quali profondità) varia la frequenza. Per localizzare le variazioni del segnale sull'asse temporale (o di profondità), si utilizza la trasformata di Fourier a finestre e la scomposizione wavelet. Questo articolo utilizza la scomposizione wavelet. La tecnologia dell'analisi wavelet ha iniziato a essere applicata attivamente nella sismica negli anni '90. Rispetto alla trasformata di Fourier a finestre, offre una risoluzione temporale migliore.
Con il seguente frammento di codice è possibile eseguire la scomposizione in componenti di una delle tracce sismiche:
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["Decomposizione wavelet",Black,32],
LabelStyle->Directive[Black,Italic],
PlotRange->All,
Frame->True]Scomposizione della traccia in componenti
Per valutare come è distribuita l'energia di riflessione nei diversi tempi di arrivo delle onde, si utilizzano le scalogramme (simili alle spettrogrammi). Di solito, nella pratica, non è necessario analizzare tutti i componenti. Solitamente si selezionano le componenti a bassa, media e alta frequenza.
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->{"Tempo","Frequenza(Hz)"},LabelStyle->Directive[Black,Bold,14],
ColorFunction->"RustTones",ImageSize->Large]Scalogramma. Risultato della funzione
In Wolfram Language, per la trasformazione wavelet si utilizza la funzione . L'applicazione di questa funzione a tutto il set di tracce avviene utilizzando la funzione . Qui è importante notare uno dei punti di forza di Wolfram Mathematica: la possibilità di utilizzare il parallelismo . Nell'esempio riportato, il parallelismo non è necessario — l'ammontare dei dati non è elevato, ma quando si lavora con set di dati sperimentali che contengono centinaia di migliaia di tracce, questo diventa una necessità.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; Dopo aver applicato la funzione si generano nuovi insiemi di dati corrispondenti alle frequenze selezionate. Nell'esempio sopra, queste sono le frequenze: 38Hz, 33Hz, 27Hz. La scelta delle frequenze avviene principalmente sulla base di test, ottenendo mappe efficaci per diverse combinazioni di frequenze e selezionando la più informativa dal punto di vista geologico.
Se è necessario condividere i risultati con i colleghi o fornire loro ai clienti, è possibile utilizzare la funzione SEGYExport[] del pacchetto GeologyIO.
outputdata=seismic3DSEGY;
outputdata["traces",1;;-1]=tracesCWD[[All,3]];
outputdata["textheader"]="Risultato della Decomposizione dell'Ondabola";
outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]];
SEGYExport["D:result.segy",outputdata];Avendo a disposizione tre di questi cubi (componenti a bassa, media e alta frequenza), si utilizza generalmente la mescolanza RGB per una visualizzazione congiunta dei dati. A ciascuna delle componenti viene assegnato un colore specifico: rosso, verde, blu. In Wolfram Mathematica, questo può essere fatto usando la funzione .
Di conseguenza, si ottengono immagini che consentono un'interpretazione geologica. I meandri, che si evidenziano nel taglio, permettono di delineare i paleo-alveo che con maggiore probabilità potrebbero essere serbatoi contenenti riserve di petrolio. La ricerca e l'analisi di analoghi moderni di tale sistema fluviale consentono di identificare le parti più promettenti dei meandri. Gli alvei propriamente detti sono caratterizzati da spessi strati di arenaria ben selezionata e costituiscono un buon serbatoio per il petrolio. Le aree al di fuori delle anomalie «a stringa» sono simili alle moderne depositi alluvionali. I depositi alluvionali sono principalmente rappresentati da rocce argillose e la perforazione in queste zone risulterà inefficace.
Sezione RGB del cubo dati. Al centro (un po' a sinistra del centro) si può vedere un fiume meandrante.
Sezione RGB del cubo dati. Nella parte sinistra si può vedere un fiume meandrante.
In alcuni casi, la qualità dei dati sismici consente di ottenere immagini notevolmente più chiare. Ciò dipende dalla metodologia delle operazioni di campo, dagli strumenti e dagli algoritmi di soppressione del rumore utilizzati. In tali casi, non sono visibili solo frammenti dei sistemi fluviali, ma anche interi e lunghi paleoreti.
Miscelazione RGB delle tre componenti del cubo di dati sismici (sezione orizzontale). Profondità di circa 2 km.
Immagine satellitare del fiume Volga nella regione di Saratov.
Conclusione
In Wolfram Mathematica è possibile analizzare i dati sismici e risolvere problemi pratici legati alla ricerca di risorse minerarie, mentre il pacchetto GeologyIO rende questo processo più agevole. La struttura dei dati sismici è tale che l'uso dei metodi integrati per accelerare i calcoli (, ,…) è molto efficace e consente di elaborare grandi volumi di dati. Questo è in gran parte dovuto alle caratteristiche di memorizzazione dei dati del pacchetto GeologyIO. Vale la pena notare che il pacchetto può essere utilizzato non solo nel campo della geofisica applicata. Tipi di dati praticamente identici vengono utilizzati nella georadar e nella sismologia. Se hai suggerimenti su come migliorare il risultato, quali algoritmi di analisi del segnale dall'arsenale di Wolfram Mathematica siano applicabili a tali dati, o hai osservazioni critiche - lascia i tuoi commenti.
Fonte: habr.com
