S3D/S3Di (Propagazione acustica ed elastica in 3D e Full Waveform Inversion)

S3Di è un programma per computer sviluppato in C ++ e dedicato alla propagazione acustica ed elastica (nel dominio del tempo) delle onde nei solidi.

Ho creato il programma per scopi di ricerca. Tuttavia, a differenza di altri codici da me sviluppati, ad oggi non ne è prevista la distribuzione.

Un'interfaccia utente grafica (GUI), sviluppata utilizzando QT, viene utilizzata per interagire con due sottoprogrammi principali:

Simulazione (Modello in avanti).

La routine può simulare sia la propagazione del campo d'onde puramente acustica che elastica, in un sottosuolo generale e/o attraverso oggetti geometricamente complessi. Il simulatore di propagazione d'onda implementa il metodo "time marching" degli elementi finiti spettrali. Ispirato dal lavoro di Komatitsch, Vilotte e Trump (1998, 1999), ho ri-derivato la formulazione matematica del problema elastodinamico e ne ho creato la mia personale implementazione in C ++ (cioè il codice è una nuova implementazione, e non un porting). La precisione del risultato è stata verificata rispetto al programma "SPECFEM3D", sviluppato in fortran da Komatitsch, Vilotte e Tromp.

L'input del programma è costituito da un insieme di file ascii, mentre la geometria della mesh può essere generata sia utilizzando CUBIT, oppure uno strumento Matlab grafico appositamente sviluppato. Se si utilizza CUBIT, la geometria può essere resa molto flessibile, aprendo una vasta gamma di opportunità per l'applicazione del presente programma in campi come ingegneria, meccanica e geotecnica. Per quanto riguarda lo strumento grafico Matlab, esso comprende la visualizzazione del modello prodotto in 3D, lo studio delle proprietà di sorgenti e ricevitori e implementa uno strumento grafico per la valutazione della condizione di stabilità del metodo degli elementi spettrali, rispetto alle proprietà spettrali del segnale della sorgente (forza eccitante). Quest'ultima capacità, a mia conoscenza, non è disponibile in alcun altro codice.   

Inversione

L'inversione del campo d'onda è implementata sotto forma di ottimizzazione locale. Come nella maggior parte dei metodi locali, una "funzione costo" (o "Energia") viene minimizzata iterativamente. A tale minimizzazione corrisponde una evoluzione dei parametri che descrivono il sottosuolo/oggetto (il modello), a partire da un'ipotesi iniziale verso una valori ottimizzati.

Il gradiente dell'Energia rispetto a un cambiamento nei parametri del modello viene calcolato tramite il l' "Adjoint method", mentre la perturbazione del modello verso parametri ottimali viene calcolato tramite il metodo di discesa del gradiente coniugato.

Alcuni dettagli sulle scelte di implementazione intraprese meritano di essere citate: tipicamente, per ottenere il gradiente dell'Energia usando l' Adjoint method, vengono calcolati il campo d'onda propagato ("in avanti", o "diretto") ed il corrispondente campo "aggiunto"; il loro prodotto è quindi integrato nello spazio e nel tempo. Le sorgenti del campo aggiunto sono tipicamente situate nelle posizioni dei ricevitori del campo d'onda diretto, e la loro variazione temporale è dominata dalla forma matematica dell'energia funzionale. Il campo aggiunto viene calcolato propagando l'effetto delle sorgenti aggiunte a ritroso nel tempo. Il principale vantaggio di questo approccio è che per ottenere il gradiente dell' Energia sono sufficienti solo due utilizzi del programma di simulazione.   

Mentre i dettagli dell' Adjoint Method possono essere trovati in Plessix (2006), è importante sottolineare che il calcolo del gradiente dell Energia richiede di mantenere in memoria (o su disco), i valori di spostamento, velocità ed accelerazione di entrambi i campi d'onda per ogni nodo della mesh, ogni componente del movimento e ogni campione temporale. La quantità di risorse diventa rapidamente proibitiva per la maggior parte dei computer.

Pertanto, ho implementato la routine di inversione in modo che una prima esecuzione del simulatore venga utilizzata per registrare lo stato del campo d'onda ad intervalli di tempo specifici. Successivamente, la propagazione in avanti ed a ritroso vengono eseguite a blocchi, per porzioni di tempo finite e limitate. Il vantaggio principale di questa strategia è che la RAM necessaria per il calcolo è di circa due ordini di grandezza inferiore, al solo costo di un unico utilizzo extra del simulatore. Di conseguenza, durante ogni iterazione le "forze aggiunte" vengono calcolate dinamicamente, senza la necessità di salvarle in file di grandi dimensioni o sospendere l'algoritmo quando quesi vengono ricaricati. Inoltre, nonostante l'approccio di inversione attualmente implementato sia posto in forma di differenza "Least square" di sismogrammi (i pieno stile FWI), la forma dell'Energia da minimizzare, e quindi l'approccio di inversione stesso, pussono essere facilmente personalizzati, semplicemente sostituendo le routine di calcolo dedicate.

Sviluppi Futuri:

  • sono in corso lavori per sviluppare una versione 2D del codice. La motivazione principale è che il calcolo 3D completo non è pratico quando il sottosuolo da modellare raggiunge grandi dimensioni, ad es. per applicazioni nella geofisica di esplorazione.
  • sono in corso ricerche per applicazioni geotecniche su piccola scala, ad es. interazione suolo-struttura e valutazione degli effetti degli interventi di compattazione del suolo.  

Riferimenti:

2006: Plessix, R.E., A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophysical Journal International,  167, p. 495–503,  DOI: 10.1111/j.1365-246X.2006.02978.x

1999: Dimitri Komatitsch and Jeroen Tromp, Introduction to the spectral-element method for 3-D seismic wave propagation, Geophysical Journal International, 139, p. 806-822

1998: Dimitri Komatitsch and Jean-Pierre Vilotte, The Spectral Element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bulletin of the Seismological Society of America88, p. 368-392