3D gravity inversion using a selection of constraints including minimum distance, smoothness and compactness

In gravity interpretation, inversion algorithms have been used widely over the years, but as the potential follows the Gauss theorem, there are many equivalent source distributions that can produce the same known field. So to obtain a unique solution, suitable constraints should be introduced to the algorithm. During the last decades many authors have used several approaches to introduce a priori information into the inversion. Green (1975) found the model closest to the initially fixed model, Last and Kubik (1983) minimized the volume of the causative body, Guillen and Menichetti (1984) concentrated the solution about a geometric element, such as an axis. Li and Oldenburg (1996, 1998) used a constraint called ‘smoothness’ to find a model with minimum spatial variation of the physical property. Also they counteracted the decreasing sensitivity of the cells with depth by weighting it with an inverse function of depth.
In this paper we have presented a method to interpret gravity data using a selection of constraints including minimum distance, smoothness and compactness that can be combined using a Lagrangian formulation. In this approach the earth is divided into a large number of rectangular prismatic blocks of fixed size where each block side is equal to the distance between two observation points and the problem has been solved by calculating the model parameters linearly (i.e. the densities of each block). Since the number of parameters can be many thousands, the linear system of equations is inverted using a conjugate gradient approach. The given weights to each block depend on depth, a priori information on density and the density ranges allowed for the region under investigation.
A MATLAB-based inversion code for the presented method was prepared. The program uses a primary density model in the input file and calculates densities of blocks at each iteration. The program was tested on two different synthetic models. The first model includes two vertical dikes with different densities and the second model has encircled multiple bodies with different geometries and densities. The results on the synthetic models seem to be acceptable with a suitable convergence. The calculated density contrasts are according to the model contrasts and the horizontal boundaries are fairly reconstructed by the algorithm. Finally the inversion procedure has been applied on the real gravity data from the Golmandareh dam site (the north-eastd Khorasan, Iran). The computations show severe karsting of the area that makes the regional stabilization uneconomical and impossible.