3D data-space gravity inversion using compactness constraint

In this paper the 3D inversion of  gravity data is considered. The goal is to reconstruct models of subsurface density distribution using a set of known gravity observations measured on the earth surface. The subsurface under the survey area is divided into large number of rectangular blocks of known sizes and positions. The unknown density contrasts within each prism define the parameters to be estimated. This kind of parameterization is flexible for the reconstruction of the subsurface model, but requires more unknown model parameters than observations (here N << M, where N is the number of data and M is the number of model parameters). The final density distribution will be obtained by minimizing a global objective function consists of data misfit and a regularization term. The inverse problem is solved in data space, which needs inverse of matrix with N×N dimension, as compared with M×M  dimension system in model space inversion. This methodology was used by Pilkington (2009) in 3D inversion of magnetic data. To solve the resulting set of linear equation, the conjugate gradient method is used. Combination of data-space method with conjugate gradient leads to keep the storage and computational time to a minimum. The iteratively-defined regularization matrix, which is used in objective function, is a combination of three diagonal matrix; namely depth weighting, compactness and hard constraint matrices. The compactness constraint was introduced in Last and Kubik (1983) and developed in Portniaguine and Zhdanov (1999), who used term "minimum support stabilizer", is considered here to produce models with non-smooth features. It is a suitable and well-known constraint for identifying geologic structures which have material properties that vary over relatively short distances. The depth weighting matrix, introduced in Li and Oldenburg (1998),  is used in regularization term to counteract the natural decay of the kernel with depth. The hard constraint allows us to incorporated priori geological and geophysical information into inversion process. While depth weighting and hard constraint matrices both are independent of the iteration index, the compactness depends on iterations. In order to recover a feasible image of the subsurface, realistic lower and upper density bounds are imposed during the inversion process. The computer program is written in MATLAB and tested on synthetic data produced by a model consists of two cubes. The cubes have same dimension and density, but located at different depths. The results indicate that the algorithm is efficient to handle large-scale gravity inverse problems. For the shallow cube the geometry and density of the reconstructed model are close to those of the original model, but for the deeper body the resolution decrease and a smooth image of  subsurface obtained. The gravity data acquired over the Safo mining camp in the north-west of Iran, which is well-known for manganese ores, are used as a real modeling case. The results show a density distribution in the subsurface from about 5 to 35-40 m in depth and about 35 m extent in the x direction, which are close to those obtained by bore-hole drilling on the site.