GNGTS 2019 - Atti del 38° Convegno Nazionale

GNGTS 2019 S essione 1.4 241 inside the voxel (e.g. Upper Crust), and a mass density ρ i . The prisms are disposed on a regular grid in Cartesian coordinates, and the forward modelling is performed in planar approximation. This geometry allows to easily introduce neighbourhood relationships. The inversion algorithm is based on the Bayes theorem: P ( x | y ) ∝ L  ( y | x ) P ( y ) (1) where y is the vector of observables, i.e. the gravity signal, x is the vector of parameters L , ρ for all voxels, P ( x | y ) and P ( x ) are the posterior and the prior distribution, respectively, and L  ( y | x ) is the likelihood. Gravity is observed, therefore the likelihood represents the degree of fit between the observed and modelled signal. The prior distribution is defined by considering the available geological and geophysical information on the study region, integrated with some regularization conditions. This information is supplied to the algorithm in the following way: • a range of variation of each boundary surface between two layers with different labels; • neighbourhood rules between the possible couple of labels; • the density of each material, i.e. of each label, in terms of the most probable value and its range of variation; • constraints on lateral and vertical variation of the density (e.g. maximum admissible lateral and vertical variation or increasing/decreasing density with depth). The shape of the prior distribution is chosen to highlight the dependency of each density ρ i on the label L i : (2) The density of each voxel ρ i is assumed to be normally distributed once the label L i is given. Its mean μ ρ ( L i ) and its variance σ  2 ρ ( L i ) are given as a-priori information. Furthermore, the distribution is truncated at μ ρ ( L i ) ± 3 σ  2 ρ ( L i ), or even with a stricter range when required to satisfy the constraints. On the other hand, the labels L are modelled as a Markov Random Field. Therefore, their probability distribution assumes the shape of a Gibbs distribution where the energy is the sum of the clique potential (Azencott, 1988), depending from two penalty functions s 2 i ( L i ) and q 2 ( L i , L j ), where L j is the label of a neighbour voxel. These functions are used to define the limits of the boundary surfaces and the neighbour rules between different materials, respectively. Moreover, q 2 ( L i , L j ) inherently imposes the regularity of the boundary surfaces between two layers. Then, invoking the MaximumA Posteriori (MAP) principle, the most probable set of labels and densities is chosen as the solution. This corresponds to find the minimum of the following target function: (3) conditioned to the constraints defined by prior information (e.g. maximum lateral density variation). Where y o is the vector of observed gravity, C –1 νν its noise covariance matrix and η , γ , and λ weights balancing the contribution of the prior terms among them and with respect to the likelihood. The minimum of Eq. 3 is retrieved by using a stochastic optimization method, i.e. a Simulated Annealing aided by a Gibbs Sampler (Robert and Casella, 2004). Gravity data and prior modelling. The model to be estimated was chosen with a horizontal resolution of 50 km × 50 km and a vertical one of 100 m. The observations to be inverted are the gravity anomalies synthesized from the GOCE-only space-wise model up to degree and order 330 (Gatti and Reguzzoni, 2017). We choose this model since a satellite-only solution can avoid the introduction of useless high frequencies that cannot be interpreted by a voxel model with the given geometrical resolution.