Bayesian gravity inversion by Monte Carlo methods

The inverse gravimetric problem consists in the reconstruction of the Earth mass density distribution from the observation of its gravitational field. The solution to this problem is generally ill-conditioned and non-unique, but, introducing very strong constraints or numerical regularization, a unique solution can be retrieved. This could be not representative of the actual mass distribution, because of non-physical or too restrictive constraints. Moreover, a strong regularization may cause very smooth solutions in terms of estimated density. In this case, the identification of the boundaries between different materials becomes a hard task. To overcome these limitations, a possible option is to apply the Bayesian approach that easily allows to introduce prior information on the unknown parameters. Moreover, secondary parameters can be easily introduced, e.g. to characterize the different types of material, allowing a sharp classification of the subsurface. In the present work, a Bayesian gravity inversion algorithm is developed, with the aim of estimating a mass density distribution together with a classification of the different types of material. The latter allows to identify the boundaries of the different materials and is developed by borrowing image analysis techniques. Consequently, the investigated volume is subdivided into volume units (voxels), each of them characterized by two random variables: the label defining the type of material (discrete) and the density (continuous). The a-priori geological information is translated in terms of this model, providing the mean density with the corresponding variability for each class of materials and the a-priori most probable label for each voxel with a set of neighbour rules. These rules have the aim to obtain a clustered model in terms of geometry with smooth boundary surfaces between the different materials. The final solution is retrieved by invoking the Maximum A Posteriori principle (MAP). The MAP is determined by applying Markov Chain Monte Carlo optimization algorithms. In particular, a simulated annealing performed by a Gibbs sampler is chosen to maximize the posterior distribution. The obtained result is then refined by a deterministic optimization algorithm. The proposed method has been implemented into a hybrid Matlab/C code. The software is then applied in the frameworks of oil exploration and crustal investigation. The former application is performed by a simulated scenario, avoiding border effects in the gravity data reduction and allowing a better understanding of the meaning of the parameters that define the prior probability. The latter is related to the inversion of the crustal structure in south China, with the aim of improving its knowledge just below a detector of neutrinos and geoneutrinos. The obtained Bayesian solution is also compared with a classical solution to the inverse gravimetric problem, showing a good agreement in terms of estimated geometry (differences with about 1 km of standard deviation). However, the Bayesian solution improve the smoothness of the estimated geometry and allows to infer horizontal and vertical density variations. In general, all the tests show that the method is able to retrieve an estimated model that is consistent with the given prior information and fits the gravity observations according to their accuracy. In this sense, one can conclude that the developed method is a sort of "artificial intelligence", supporting the user by automatizing the trial-and-error techniques. In fact, an initial guess on the model parameters is required, but differently from a trial-and-error solution these parameters are automatically adapted to fit the gravity observations, according to the rules introduced by the user as prior probability. If the solution is not satisfactory, this should be attributed to the weak information provided by the gravity or to the wrong or incompatible geological information supplied by the user.

Il problema gravimetrico inverso consiste nella ricostruzione della distribuzione della densità di massa della Terra a partire dall'osservazione del suo campo gravitazionale. La soluzione a questo problema è generalmente mal posta e non unica, quindi può essere ottenuta solo introducendo dei vincoli molto restrittivi oppure attraverso una forte regolarizzazione numerica. Purtroppo, una soluzione così ottenuta potrebbe non rappresentare una distribuzione realistica delle masse, dato che i vincoli introdotti potrebbero essere troppo restrittivi o incompatibili con la geologia della zona in esame. Inoltre, la soluzione regolarizzata potrebbe risultare troppo liscia nelle densità stimate, non permettendo una chiara identificazione dei bordi tra i differenti materiali sulla base della densità. Per superare questi problemi, una possibilità è l’approccio Bayesiano. Quest’ultimo permette facilmente di introdurre sia delle informazioni a priori sui parametri del modello, sia dei parametri secondari che garantiscano una classificazione dei diversi materiali che compongono il sottosuolo. Nel presente lavoro viene sviluppato un algoritmo Bayesiano di inversione gravimetrica, mutuando alcune tecniche dell’analisi di immagine per la classificazione. Questa scelta permette di suddividere il volume esaminato in unità di volume (voxel) ciascuna caratterizzata da due variabili casuali: un'etichetta che definisce il tipo di classe geologica (discreta) ed una densità (continua). Le informazioni geologiche a priori devono essere tradotte secondo questo modello, fornendo la densità media e la relativa variabilità per ciascuna delle possibili classi geologiche, l'etichetta più probabile a priori per ogni voxel ed un insieme di regole di vicinanza. Queste regole servono ad ottenere un modello che sia clusterizzato e che presenti delle superfici di confine tra i diversi materiali lisce. La soluzione finale viene stimata richiamando il principio di Massimo A Posteriori (MAP). Il MAP è determinato tramite algoritmi di ottimizzazione basati su metodi Monte Carlo e delle catene di Markov. In particolare, la massimizzazione della distribuzione di probabilità a posteriori avviene tramite un Gibbs sampler innestato in un simulated annealing. La stima così ottenuta viene successivamente affinata con un algoritmo di ottimizzazione deterministica. Il metodo proposto è stato implementato in un software scritto in linguaggio di programmazione Matlab, integrato con delle funzioni sviluppate in C. Il software viene applicato in due casi di studio relativi a due contesti differenti: l’esplorazione geofisica per la ricerca di idrocarburi e la determinazione della struttura crostale. Nel primo caso si tratta di uno scenario simulato, così da evitare effetti di bordo nella riduzione dei dati di gravità e consentendo una migliore comprensione del significato dei parametri usati per definire la probabilità a priori. Il secondo caso è legato determinazione della struttura crostale nella Cina meridionale, con lo scopo di migliorarne la conoscenza in corrispondenza di un rilevatore di neutrini e geoneutrini. La soluzione Bayesiana ottenuta è stata confrontata con una soluzione classica al problema gravimetrico inverso, mostrando un buon accordo in termini di stima delle geometrie (deviazione standard delle differenze di circa 1 km). La stima Bayesiana permette comunque di ottenere una geometria più liscia e di stimare anche delle variazioni laterali e verticali della densità. In generale, tutti i test mostrano che il metodo è in grado di stimare un modello che si adatti alle informazioni a priori e alle osservazioni di gravità, in base alla precisione loro assegnata. In questo senso si può concludere che il metodo sviluppato è una sorta di “intelligenza artificiale”, che supporta l'utente automatizzando le tecniche di trial-and-error. Infatti, si richiede che l'utente introduca una stima iniziale dei parametri del modello, ma, al contrario dei metodi trial-and-error, questi vengono poi automaticamente adattati al fine di spiegare correttamente le osservazioni di gravità secondo le regole introdotte dall'utente come probabilità a-priori. Se la soluzione ottenuta non è soddisfacente, ciò è da attribuirsi alla debolezza dalla gravità o alle informazioni geologiche a-priori introdotte in maniera errata dall’utente o incompatibili con la gravità.