  We propose an approach to infer large-scale heterogeneities within a small
celestial body from measurements of its gravitational potential, provided for
instance by spacecraft radio-tracking. The non-uniqueness of the gravity
inversion is here mitigated by limiting the solutions to piecewise-constant
density distributions, thus composed of multiple regions of uniform density
(mass anomalies) dispersed in a background medium. The boundary of each anomaly
is defined implicitly as the 0-level surface of a scalar field (called the
level-set function), so that by modifying this field the shape and location of
the anomaly are varied. The gravitational potential associated with a density
distribution is here computed via a line-integral polyhedron method, yielding
the coefficients of its spherical harmonics expansion. The density distribution
is then adjusted via an iterative least-squares approach with Tikhonov
regularization, estimating at every iteration corrections to the level-set
function, the density contrast of each anomaly, and the background density, in
order to minimize the residuals between the predicted gravity coefficients and
those measured. Given the non-convexity of the problem and the lack of prior
knowledge assumed (save for the shape of the body), the estimation process is
repeated for several random initial distributions, and the resulting solutions
are clustered based on global properties independent of the input measurements.
This provides families of candidate interior models in agreement with the data,
and the spread of the local density values across each family is used to assess
the uncertainties associated with the estimated distributions.
