Ringraziamo l'autore del blog per la sua presentazione
Introduzione
Questa nota è scritta a seguito della conferenza e contiene un riassunto della presentazione che ho tenuto. L'evento si è svolto a giugno nella città di San Pietroburgo. Considerando che lavoro a un quartiere dal luogo della conferenza, non potevo non partecipare a questo evento. Nel 2016 e 2017 ho ascoltato le presentazioni della conferenza, e quest'anno ho presentato io stesso. Prima di tutto, c'era un argomento interessante (che ritengo) che stiamo sviluppando con , e in secondo luogo, dopo un lungo studio della legislazione russa riguardo la politica sanzionatoria, nell'azienda in cui lavoro sono state ottenute ben due licenze .
Prima di passare all'argomento della mia presentazione, vorrei sottolineare l'ottima organizzazione dell'evento. Nella pagina di benvenuto 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 erano accolti da assistenti tra gli studenti – non permettevano di perdere la strada. Durante la registrazione venivano distribuiti piccoli gadget (giocattolo – un lampo di Spikey, penna, adesivi con la simbologia di Wolfram). Il pranzo e il coffee break erano anche inclusi nel programma della conferenza. Ho già accennato al buon caffè e ai dolcetti nel gruppo – i cuochi sono bravi. Con questa introduzione, vorrei sottolineare che l'evento stesso, il suo formato e il luogo di svolgimento portano già emozioni positive.
La presentazione che ho preparato insieme a Kirill Belov si chiama «Uso di Wolfram Mathematica per risolvere problemi di geofisica applicata. Analisi spettrale dei dati sismici o 'dove scorrevano i fiumi antichi'. Il contenuto della presentazione comprende due parti: prima di tutto, l'uso degli algoritmi disponibili in per l'analisi dei dati geofisici, e in secondo luogo, come integrare i dati geofisici in Wolfram Mathematica.
Sismica
Per cominciare, è necessario fare una breve introduzione alla geofisica. La geofisica è la scienza che studia le proprietà fisiche delle rocce. Poiché le rocce hanno proprietà diverse: elettriche, magnetiche, elastiche, esistono metodi geofisici corrispondenti: prospezione elettrica, prospezione magnetica, prospezione sismica… Nel contesto di questo articolo parleremo più nel dettaglio solo della prospezione sismica. La prospezione sismica è il principale metodo per la ricerca di petrolio e gas. Questo metodo si basa sull'innesco di oscillazioni elastiche e sulla successiva registrazione della risposta delle rocce che compongono il territorio in esame. Le oscillazioni sono generate sulla terraferma (con dinamite o fonti di vibrazione non esplosive) o in mare (con cannoni pneumatici). Le onde elastiche si propagano attraverso gli strati di rocce 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). Dall'orario di arrivo delle onde è possibile dedurre la profondità degli strati geologici.
La nave sismica traina l'attrezzatura
Il cannone pneumatico genera oscillazioni elastiche
Le onde attraversano gli strati di rocce e vengono registrate dagli idrofoni
La nave di ricerca geofisica "Ivan Gubkin" ormeggiata presso il ponte di Blagoveščensk a Pietroburgo
Modello del segnale sismico
Le rocce hanno diverse proprietà fisiche. Per la prospezione sismica, le proprietà elastiche sono soprattutto importanti: la velocità di propagazione delle oscillazioni elastiche e la densità. Se due strati hanno proprietà uguali o simili, l'onda "non riconoscerà" il confine tra di esse. Se le velocità delle onde negli strati sono diverse, allora ci sarà una riflessione al confine degli strati. Maggiore è la differenza nelle proprietà, più intensa sarà la riflessione. La sua intensità sarà determinata dal coefficiente di riflessione (rc):
dove ρ è la densità delle rocce, ν è la velocità delle onde, 1 e 2 indicano lo strato superiore e quello inferiore.
Uno dei modelli di segnale sismico più semplici e comunemente utilizzati è il modello di convoluzione, in cui la traccia sismica registrata è rappresentata come il risultato della convoluzione di una sequenza di coefficienti di riflessione con un impulso di sondaggio:
dove s(t) è la traccia sismica, ovvero tutto ciò che ha registrato un idrofono o un geofono durante un periodo fisso di registrazione, w(t) è il segnale generato da un cannone pneumatico, n(t) è il rumore casuale.
Calcoliamo, per esempio, una traccia sismica sintetica. Come segnale di partenza utilizzeremo l'impulso di Ricker, ampiamente utilizzato 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["Wavelet iniziale",Black,20],
LabelStyle->Directive[Black,Italic],
FillingStyle->{White,Black},ImageSize->Large,InterpolationOrder->2]Impulso sismico iniziale
Stabilirà 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 di coefficienti di riflessione
Calcoliamo e visualizziamo la traccia sismica. Poiché i coefficienti di riflessione hanno segni diversi, anche sulla traccia sismica otteniamo due riflessioni alternative.
traceExamle=ListConvolve[wavelet[[1;;;;1]],rcExample];
ListPlot[traceExamle,
PlotStyle->Black,Filling->0,Frame->True,Axes->False,
PlotLabel->Style["Traccia sismica",Black,20],
LabelStyle->Directive[Black,Italic]]Traccia simulata
Per questo esempio bisogna fare una precisazione: in realtà la profondità dei sedimenti è definita in metri, mentre il calcolo della traccia sismica avviene nell'area temporale. Sarebbe più corretto definire le profondità in metri e calcolare i tempi di arrivo conoscendo le velocità nei sedimenti. In questo caso, ho subito definito i livelli sull'asse temporale.
Quando si parla di ricerche di campo, si registra un'enorme quantità di serie temporali simili (tracce sismiche) come risultato di tali osservazioni. Ad esempio, nell'analisi di un'area lunga 25 km e larga 15 km, dove ogni traccia caratterizza una cella delle dimensioni di 25x25 metri (questa cella è chiamata bin), il file finale dei dati conterrà 600.000 tracce. Con un passo di campionamento nel tempo pari a 1 ms e un tempo di registrazione di 5 secondi, il file finale dei dati avrà una dimensione superiore a 11 GB, mentre il volume del materiale
Come lavorare con loro in ?
Pacchetto
L'inizio dello sviluppo del pacchetto è stato sulla bacheca del gruppo VK di supporto in lingua russa. Grazie alle risposte della comunità, la soluzione è stata trovata molto rapidamente. E questo si è evoluto in un serio sviluppo. Il è stato persino evidenziato dai moderatori. Attualmente, il pacchetto supporta il lavoro con i seguenti tipi di dati, ampiamente utilizzati nel settore geologico:
- importazione di dati cartografici nei formati ZMAP e IRAP
- importazione di misurazioni nei pozzi nel formato LAS
- input e output di file sismici nel formato
Per installare il pacchetto, è necessario 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 principali funzionalità del pacchetto. La chiamata viene effettuata in modo tradizionale per i pacchetti in Wolfram Language:
Get["GeologyIO`"]Il pacchetto è sviluppato utilizzando . Questo consente di accompagnare la funzionalità principale del pacchetto con una documentazione che non differisce nel formato da quella di Wolfram Mathematica e di fornire al pacchetto file di test per un primo approccio.
Un file di questo tipo è il file «Marmousi.segy» - un modello sintetico di un profilo geologico sviluppato dall'Institut Français du Pétrole. Utilizzando questo modello, gli sviluppatori testano i propri algoritmi di modellazione del campo d'onda, elaborazione dei dati, inversone delle tracce sismiche, ecc. Il modello Marmousi è conservato in un repository, da cui è stato scaricato il pacchetto stesso. 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. Innanzitutto, ci sono commenti testuali. Qui vengono riportate informazioni sul luogo delle operazioni, i nomi delle aziende che hanno effettuato le misurazioni, ecc. Nel nostro caso, questo intestazione è richiamata da una richiesta con la chiave TextHeader. Qui è riportato un'intestazione testuale abbreviata:
Short[marmousi["TextHeader"]]«Il set di dati Marmousi è stato generato presso l'Istituto … la velocità minima di 1500 m/s e una massima di 5500 m/s)»
È possibile visualizzare il modello geologico accedendo alle tracce sismiche con la chiave «traces» (una delle caratteristiche del pacchetto è l'indipendenza delle chiavi dalle maiuscole):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Il modello Marmousi
Attualmente, il pacchetto consente anche di caricare dati in parti da file grandi, il che rende possibile l'elaborazione di file le cui dimensioni possono raggiungere decine di gigabyte. Inoltre, tra le funzioni del pacchetto ci sono quelle per esportare dati in .segy e la scrittura parziale alla fine del file.
Vale la pena notare la funzionalità del pacchetto nel lavorare con strutture complesse di file .segy. Infatti, consente non solo di accedere a singole tracce, intestazioni tramite chiavi e indici, ma anche di modificarle con successiva scrittura nel file. Molti dettagli tecnici sull'implementazione di GeologyIO vanno oltre l'ambito di questo articolo e meritano probabilmente una descrizione separata.
La rilevanza dell'analisi spettrale nella sismica
La possibilità di importare materiali sismici in Wolfram Mathematica consente di utilizzare la funzionalità integrata per l'elaborazione dei segnali sui dati sperimentali. Poiché ogni traccia sismica rappresenta una serie temporale, uno degli strumenti principali per il loro studio è l'analisi spettrale. Tra le premesse per l'analisi della composizione in frequenza dei dati sismici si possono citare, ad esempio, le seguenti:
- I diversi tipi di onde sono caratterizzati da diverse composizioni frequenziali. Questo consente di separare onde utili e di sopprimere onde di disturbo.
- Proprietà delle rocce, come la porosità e la saturazione, possono influenzare la composizione frequenziale. Questo permette di isolare rocce con migliori caratteristiche.
- Strati di diversa altezza causano anomalie in diversi intervalli di frequenza.
Il terzo punto è fondamentale nel contesto di questo articolo. Di seguito è riportato un frammento di codice per calcolare le tracce sismiche nel caso di uno strato con spessore variabile — un modello a cuneo. Questo modello è tradizionalmente studiato nella sismica per analizzare gli effetti di interferenza, quando le onde riflesse da molti strati si sovrappongono.
nx=200;(* Numero di punti di griglia nella direzione X*)
ny=200;(* Numero di punti di 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³, usata come densità costante *)
Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}];
Table[modellv[[;;70,i]]=4500;,{i,1,200}];
(* Plotting 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, all'esterno del cuneo è di 4000 m/s, mentre la densità è considerata costante a 2200 g/cm³. Per questo modello calcoliamo 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]Tracce sismiche per il modello a cuneo
La sequenza di traiettorie sismiche rappresentata in questa immagine è chiamata sezione sismica. Come può essere notato, la sua interpretazione può avvenire anche a livello intuitivo, poiché la geometria delle onde riflesse corrisponde chiaramente al modello precedentemente definito. Analizzando più nel dettaglio le traiettorie, si può notare che le traiettorie dalla 1° alla 30° non presentano differenze: il riflesso dal tetto del strato e dal fondo non si sovrappongono. A partire dalla 31° traiettoria, i riflessi iniziano a interferire. E, sebbene nel modello i coefficienti di riflessione non cambino orizzontalmente, le traiettorie sismiche modificano la loro intensità al variare dello spessore dello strato.
Consideriamo l'ampiezza del riflesso dal confine superiore dello strato. A partire dalla 60° traiettoria, l'intensità del riflesso inizia ad aumentare e alla 70° traiettoria raggiunge il massimo. Questo è come si manifesta l'interferenza delle onde dal tetto e dal fondo degli strati, portando in alcuni casi a significative anomalie nella registrazione sismica.
ListLinePlot[GaussianFilter[Abs[traces[[All,46]]],3][[;;;;2]],
InterpolationOrder->2,Frame->True,PlotStyle->Black,
PlotLabel->Style["Ampiezza del riflesso",Black,20],
LabelStyle->Directive[Black,Italic],
PlotRange->All]Grafico dell'ampiezza dell'onda riflessa dal bordo superiore del cuneo
È logico che quando il segnale è a bassa frequenza, l'interferenza inizia a manifestarsi a spessori maggiori dello 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["Set di wavelet",Black,20],
LabelStyle->Directive[Black,Italic],
ImageSize->Large,InterpolationOrder->2]Insieme di segnali originali con frequenze di 35 Hz, 55 Hz, 85 Hz
Dopo aver calcolato le traiettorie sismiche e costruito i grafici dell'ampiezza dell'onda riflessa, possiamo vedere che per diverse frequenze l'anomalia si osserva a vari spessori dello 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 dell'onda riflessa dal bordo superiore del cuneo per diverse frequenze
La possibilità di trarre conclusioni sulla profondità del giacimento dai risultati delle osservazioni sismiche è estremamente utile, poiché uno dei compiti principali nella ricerca di giacimenti di petrolio è valutare i punti più promettenti per la perforazione (cioè le aree in cui il giacimento ha uno spessore maggiore). Inoltre, nel profilo geologico possono trovarsi oggetti che, per la loro genesi, comportano un cambiamento brusco dello spessore del giacimento. Questo rende l'analisi spettrale uno strumento efficace per il loro studio. Nella parte successiva dell'articolo esamineremo in dettaglio questi oggetti geologici.
Dati sperimentali. Dove sono stati ottenuti e cosa cercare in essi?
I materiali analizzati nell'articolo sono stati ottenuti nel territorio della Siberia Occidentale. La regione, come probabilmente sanno tutti senza eccezioni, è il principale territorio produttore di petrolio del nostro paese. Lo sviluppo attivo dei giacimenti è iniziato nella regione negli anni '60 del secolo scorso. Il principale metodo di ricerca dei giacimenti petroliferi è la sismica. È interessante osservare le immagini satellitari di questo territorio. Su scala ridotta si può notare un'enorme quantità di paludi e laghi; ingrandendo la mappa si possono vedere le piattaforme di perforazione e, ingrandendo ulteriormente, si riescono a distinguere anche i sentieri dei profili lungo i quali sono state effettuate le osservazioni sismiche.
Immagine satellitare delle mappe Yandex — area della città di Noyabrsk
Rete di piattaforme su uno dei giacimenti
Le rocce contenenti petrolio della Siberia Occidentale giacciono a profondità che variano da 1 km a 5 km. La maggior parte delle rocce contenenti petrolio si è formata durante il Giurassico e il Cretaceo. Il periodo giurassico è probabilmente noto a molti per il film omonimo. era significativamente diverso da quello attuale. Nell'Enciclopedia Britannica esiste una serie di paleocarte che caratterizzano ogni epoca geologica.
Tempo presente
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ò supporre che il paesaggio tipico di quel tempo apparisse nel seguente modo:
Siberia del periodo Giurassico
In questa immagine, per noi non contano tanto gli animali e gli uccelli, quanto l'immagine del fiume sullo sfondo. Il fiume è proprio quell'oggetto geologico su cui ci siamo soffermati in precedenza. Infatti, l'attività dei fiumi consente l'accumulo di sabbie ben classificate, che in seguito diventeranno serbatoi per il petrolio. Questi serbatoi possono avere forme bizzarre e complesse (così come il letto del fiume) e hanno uno spessore variabile: vicino alle sponde lo spessore è ridotto, mentre più verso il centro del letto o in tratti meandrati aumenta. Quindi, i fiumi formati durante il Giurassico si trovano ora a una profondità di circa tre chilometri e sono oggetto di ricerca di serbatoi di petrolio.
Dati sperimentali. Elaborazione e visualizzazione
Facciamo subito una premessa riguardo ai materiali sismici mostrati nell'articolo: data la grande quantità di dati utilizzati per l'analisi, nel testo dell'articolo è stato inserito solo un frammento del set originale di tracciati sismici. Questo permetterà a chiunque di riprodurre i calcoli presentati.
Lavorando con dati sismici, i geofisici utilizzano generalmente software specializzati (ci sono alcuni leader di settore, le cui soluzioni sono ampiamente utilizzate, come Petrel o Paradigm), che permette di analizzare diversi tipi di dati e presenta un'interfaccia grafica intuitiva. Nonostante la loro comodità, tali software hanno anche i loro svantaggi: l'implementazione di algoritmi moderni nelle versioni stabili richiede molto tempo e le possibilità di automazione dei calcoli sono generalmente limitate. In questa situazione, l'utilizzo di sistemi di matematica computazionale e linguaggi di programmazione ad alto livello si rivela molto utile, poiché consentono di utilizzare una vasta base di algoritmi e, allo stesso tempo, si occupano di molte attività ripetitive. Su questo principio si basa anche il lavoro con i dati sismici in Wolfram Mathematica. Non è sensato scrivere una funzionalità ricca per un’interazione interattiva con i dati — è più importante garantire il caricamento da un formato comune, applicare gli algoritmi desiderati e scaricare nuovamente 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 importati e caricati in questo modo sono tracciati registrati su un'area di 10 per 5 chilometri. Nel caso in cui i dati siano stati ottenuti utilizzando la metodologia di sismica 3D (in cui le registrazioni delle onde vengono effettuate non lungo singoli profili geofisici, ma su tutta l'area contemporaneamente), è possibile ottenere cubi di dati sismici. Questi sono oggetti tridimensionali, i cui tagli verticali e orizzontali consentono di studiare dettagliatamente l'ambiente geologico. Nell'esempio analizzato, stiamo trattando proprio dati tridimensionali. Alcune informazioni possono essere ricavate 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: UTENTE WOLFRAM
C 6 NOME DELLA RICERCA: DA qualche parte in SIBERIA
C 7 TIPO DI FILE VOLUME SISMICO 3D
C 8
C 9
C 10 INTERVALLO Z: PRIMO 2200M ULTIMO 2400M
Questo insieme 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 approssimativamente come nella seguente illustrazione: si tratta della distribuzione delle ampiezze delle onde riflesse lungo l'asse verticale (asse delle profondità).
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]Una delle tracce del profilo sismico
Conoscendo il numero di tracce situate in ciascuna direzione dell'area studiata, è possibile formare un array tridimensionale di dati e visualizzarlo utilizzando la funzione Image3D[]
traces=seismic3DSEGY["traces"];
startIL=1050;EndIL=2000;stepIL=2; (*coordinata X dell'inizio e della fine della rilevazione e passo delle tracce*)
startXL=1165;EndXL=1615;stepXL=2; (*coordinata Y dell'inizio e della fine della rilevazione e passo delle tracce*)
numIL=(EndIL-startIL)/stepIL+1; (*numero di tracce lungo l'asse X*)
numXL=(EndXL-startXL)/stepIL+1; (*numero di tracce lungo l'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 del cubo di dati sismici. (L'asse verticale è la profondità)
Nel caso in cui gli oggetti geologici di interesse generino anomalie sismiche intense, è possibile utilizzare strumenti di visualizzazione con trasparenza. Le sezioni 'non importanti' della registrazione possono essere rese invisibili, lasciando visibili solo le anomalie. In Wolfram Mathematica ciò 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 del cubo di dati sismici utilizzando le funzioni Opacity[] e Raster3D[]
Come nell'esempio sintetico, nei tagli del cubo originale è possibile evidenziare alcuni confini geologici (strati) con topografia variabile.
Lo strumento principale per l'analisi spettrale è la trasformata di Fourier. Con essa è possibile valutare lo spettro ampiezza-frequenza di ogni tracciato o gruppo di tracciati. Tuttavia, dopo la conversione dei dati nell'area delle frequenze, si perde l'informazione su quali tempi (leggi, a quali profondità) varia la frequenza. Per poter localizzare le variazioni del segnale sull'asse temporale (di profondità), si utilizzano la trasformata di Fourier a finestra e la scomposizione wavelet. In questo articolo si utilizza la scomposizione wavelet. La tecnologia di analisi wavelet è stata attivamente applicata nella sismica negli anni '90. Un vantaggio rispetto alla trasformata di Fourier a finestra è una migliore risoluzione temporale.
Con il seguente frammento di codice è possibile eseguire la scomposizione in singole componenti di uno dei tracciati sismici:
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["Scomposizione wavelet",Black,32],
LabelStyle->Directive[Black,Italic],
PlotRange->All,
Frame->True]Scomposizione del tracciato in componenti
Per valutare come è distribuita l'energia di riflessione nei diversi tempi di arrivo delle onde si utilizzano le scalogrammi (l'analogo della spettrogramma). In pratica, di solito non è necessario analizzare tutte le componenti. Si scelgono generalmente 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 trasformata wavelet si utilizza la funzione . L'applicazione di questa funzione all'intero set di tracciati si realizza utilizzando la funzione . Qui è opportuno notare uno dei punti di forza di Wolfram Mathematica: la possibilità di utilizzare il parallelismo . Nell'esempio fornito non è necessaria la parallelizzazione — il volume dei dati non è elevato, ma quando si lavora con set di dati sperimentali contenenti centinaia di migliaia di tracciati, 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 Appaiono nuovi set di dati corrispondenti alle frequenze selezionate. Nell'esempio sopra, queste frequenze sono: 38Hz, 33Hz, 27Hz. La selezione delle frequenze avviene più spesso sulla base di test: si ottengono mappe efficaci per diverse combinazioni di frequenze e si sceglie la più informativa dal punto di vista del geologo.
Se è necessario condividere i risultati con i colleghi o fornirli al cliente, è possibile utilizzare la funzione SEGYExport[] del pacchetto GeologyIO.
outputdata=seismic3DSEGY;
outputdata["traces",1;;-1]=tracesCWD[[All,3]];
outputdata["textheader"]="Risultato della decomposizione dell'onda";
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), in genere si utilizza la mescolanza RGB per la visualizzazione congiunta dei dati. A ciascuna delle componenti viene assegnato un colore: rosso, verde, blu. In Wolfram Mathematica questo può essere fatto utilizzando la funzione .
Ne risultano immagini su cui è possibile eseguire un'interpretazione geologica. I meandri, che vengono visualizzati nella sezione, consentono di delineare i paleocorsi d'acqua che con maggiore probabilità possono essere serbatoi e contenere riserve di petrolio. La ricerca e l'analisi di analoghi moderni di tale sistema fluviale permette di determinare le parti più promettenti dei meandri. I corsi d'acqua propriamente detti sono caratterizzati da spesse stratificazioni di arenaria ben selezionata e rappresentano un buon serbatoio di petrolio. Le aree al di fuori delle anomalie 'cordonate' sono simili agli attuali depositi alluvionali. I depositi alluvionali sono principalmente costituiti da materiali argillosi e la perforazione in queste zone sarà inefficace.
Taglio RGB del cubo di dati. Al centro (un po' a sinistra del centro) si può seguire il fiume meandrante.
Taglio RGB del cubo di dati. Nella parte sinistra si può seguire il fiume meandrante.
In alcuni casi, la qualità dei dati sismici consente di ottenere immagini notevolmente più chiare. Ciò dipende dalla metodologia dei lavori di campo, dall'attrezzatura e dagli algoritmi di riduzione del rumore utilizzati. In tali casi, sono visibili non solo frammenti di sistemi fluviali, ma anche interi paleoflussi prolungati.
Mescolanza RGB delle tre componenti del cubo di dati sismici (sezione orizzontale). Profondità di circa 2 km.
Immagine satellitare del fiume Volga nella zona di Saratov.
Conclusione
In Wolfram Mathematica è possibile analizzare i dati sismici e risolvere problemi pratici legati alla ricerca di risorse minerarie, e il pacchetto GeologyIO rende questo processo più comodo. 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. Non poco contribuiscono a questo le caratteristiche di archiviazione dei dati del pacchetto GeologyIO. A dire il vero, il pacchetto può essere utilizzato non solo nel campo della prospezione sismica applicata. Tipi di dati praticamente identici sono utilizzati nella geo-radar e nella sismologia. Se avete suggerimenti su come migliorare i risultati, quali algoritmi di analisi del segnale dell'arsenale di Wolfram Mathematica sono applicabili a questi dati, o se avete commenti critici - lasciate un commento.
Fonte: habr.com
