Confronto di algoritmi Markov Chain Monte Carlo con applicazioni a distribuzioni di probabilità analitiche e a tomografia sismica a trasmissione

Uno degli aspetti cruciali di ogni problema inverso è la quantificazione dell’incertezza connessa ai parametri del modello stimati. A questo fine il problema inverso viene spesso formulato in un’ottica Bayesiana, in cui le informazioni a priori sui parametri del modello, combinandosi con le informazioni fornite dal dato osservato su tali parametri (la così detta likelihood function), danno vita alla distribuzione a posteriori; tale distribuzione sarà la soluzione finale del problema inverso Bayesiano. Le caratteristiche statistiche della distribuzione a posteriori (es. media e varianza) possono essere analiticamente calcolate solo nel caso di distribuzioni Gaussiane dei parametri del modello e dei dati e solo per problemi inversi lineari. In tutti gli altri casi per una stima accurata della distribuzione a posteriori è necessario utilizzare metodi numerici. I più utilizzati tra essi sono gli algoritmi Markov Chain Monte Carlo (MCMC), che trasformano un problema inverso in un problema di campionamento nello spazio dei modelli, in cui la densità di campionamento è pari alla distribuzione a posteriori desiderata. Tali metodi richiedono di solito un elevato costo computazionale, dovuto ai numerosi modelli (e quindi valutazioni del forward modelling) che è necessario campionare per giungere a stime stabili della distribuzione a posteriori. Tale costo computazionale aumenta con il numero di parametri del modello (problema del curse of dimensionality) e con la complessità (es. multi-modalità) della distribuzione a posteriori desiderata. Per diminuire tale sforzo computazionale numerosi metodi MCMC e diverse strategie di campionamento sono state proposte nel corso degli anni. In questa Tesi si testeranno 5 metodi MCMC, con l’intento di analizzare le loro capacità di fornire stime affidabili della distribuzione a posteriori e la loro velocità di convergenza. I metodi testati sono: Random Walk Metropolis (RWM), Adaptive Metropolis (AM), Metropolis algorithm with adaptation of the scaling factor (AM_sd), Differential Evolution Markov Chain (DE-MC) e Differential Evolution Adaptive Metropolis (DREAM). Tali metodi saranno testati sia su distribuzioni di probabilità analitiche, sia su un problema geofisico altamente non-lineare come la tomografia a trasmissione 2D cross-hole, che sarà effettuata su dati sismici sintetici ricavati su un modello di sottosuolo semplificato. Tutti gli algoritmi sono stati implementati in Matlab e le loro performances di convergenza saranno valutate tramite il calcolo del potential scale reduction factor, che ci offre una stima del numero di modelli che è necessario campionare per giungere ad una stima accurata della distribuzione a posteriori. Per migliorare l’esplorazione dello spazio dei modelli ciascun algoritmo farà uso di più catene MCMC, i cui modelli, campionati dopo la fase iniziale di campionamento, detta burn-in, saranno utilizzati per il calcolo della distribuzione a posteriori (es. media, varianza, distribuzioni di probabilità marginali). Una differenza sostanziale tra i vari metodi testati sta proprio sul loro uso delle diverse catene: gli algoritmi RWM, AM e AM_sd le catene evolvono in maniera totalmente indipendente senza nessuna comunicazione tra di esse (quindi senza nessuna combinazione “mixing” dei modelli campionati dalle varie catene). Diversamente DE-MC e DREAM combinano i modelli estratti da ciascuna catena per migliorare il campionamento dello spazio dei modelli e quindi per velocizzare la convergenza verso la distribuzione a posteriori. Per quanto riguarda le distribuzioni analitiche, sono state considerate due distribuzioni Gaussian-mixture sia in una che due dimensioni, sia con mode vicine sia con mode separate da una zona di bassa probabilità. Questi test analitici hanno un basso costo computazionale (non vi è in sostanza nessun operatore di forward modelling) e quindi permettono una valutazione rapida ed accurata sia delle performances dei vari metodi analizzati che delle loro potenzialità e criticità. La tomografia cross-hole è invece computazionalmente molto più onerosa perché richiede la soluzione numerica dell’equazione eikonale per mezzi acustici per ogni modello campionato. Le maggiori difficoltà del problema risiedono sia nel forte mal-condizionamento del problema inverso (cioè più modelli saranno caratterizzati da valori simili della funzione di likelihood), sia nel suo elevato numero di incognite (cioè il numero di celle che formano il modello di velocità compressionale). In tale contesto la distribuzione a posteriori cercata sarà altamente localizzata in uno spazio ad elevate dimensioni e quindi sarà difficilmente campionabile. Nel primo test tomografico si è cercato di attenuare tali problemi andando ad utilizzare due griglie per la parametrizzazione del sottosuolo: una più lasca, in cui ogni nodo costituisce una incognita del modello, ed una più fitta, ottenuta interpolando quella più lasca, che sarà utilizzata per la fase di forward modelling. In questo primo test si trascura inoltre ogni tipo di correlazione spaziale a priori tra le velocità del modello; in sostanza si considera che la velocità P di ogni cella sia totalmente indipendente da quella delle celle ad essa adiacenti. Anche in questa configurazione il problema risulta ancora fortemente mal-condizionato e quindi i vari algoritmi mostrano seri problemi di convergenza. Pertanto, si è adottata anche un’altra strategia di inversione in cui si utilizza un'unica griglia di parametrizzazione, sia per definire le incognite, sia per la fase di forward modelling, ma in questo caso si va ad inserire nel modello a priori l’informazione circa la correlazione spaziale delle velocità delle singole celle. In questo modo si va a ridurre drasticamente il mal-condizionamento del problema e si facilita pertanto la fase di campionamento perché gli algoritmi andranno ad esplorare solo quei modelli di sottosuolo in accordo con il modello di correlazione spaziale inserito. In tutti i test tomografici si considera una distribuzione del noise e dei parametri del modello Gaussiana. In aggiunta oltre al potential scale reduction factor, e al confronto tra il modello vero e le caratteristiche statistiche delle distribuzioni a posteriori fornite dei vari metodi, si è valutata l’attendibilità dei risultati andando a confrontare i risultati prodotti per diverse configurazioni di sorgenti e ricevitori e verificando che la deviazione standard sui parametri del modello fosse coerente con i principi di illuminazione sismica del modello di sottosuolo. Nei test analitici 1D o 2D tutti i metodi campionano accuratamente la distribuzione di probabilità desiderata quando essa è una Gaussian-mixture caratterizzata da mode della distribuzione molto vicine, l’unica differenza risiede nella velocità di convergenza: i migliori risultati sono forniti dai metodi DREAM e DE-MC, seguiti poi dagli algoritmi AM e AM_sd e in ultimo dal metodo RWM. Le differenze nei risultati forniti dai vari metodi aumentano, però, all’aumentare della complessità della distribuzione desiderata come confermano i test eseguiti sulle distribuzioni Gaussian-mixture caratterizzate da mode separati da una zona di bassa probabilità. Nel caso del problema tomografico per ridurre il costo computazionale dovuto al confronto di ben 5 metodi, si è deciso di sfruttare i risultati ottenuti nei test analitici in cui le performance di AM_sd e AM, da un lato, e di DE-MC e DREAM dall’altro, erano molto simili- Pertanto si è scelto in questo caso di confrontare solo DE-MC, AM e RWM. Questi test confermano le migliori performance di convergenza di DE-MC rispetto ad AM e soprattutto rispetto a RWM, sebbene tutti i metodi forniscano stime comparabili delle caratteristiche statistiche della distribuzione a posteriori desiderata (es. media e varianza). Tali caratteristiche, a loro volta, sono in accordo con l’illuminazione sismica del sottosuolo determinata dalla geometria di acquisizione (es. si ha maggiore varianza a posteriori dove l’illuminazione è minore). Infine, l’inserimento del modello di correlazione spaziale nella distribuzione a priori si è dimostrata una strategia vincente per diminuire drasticamente il mal-condizionamento del problema e per migliorare, dunque, le performance di tutti gli algoritmi