We have explored the feasibility of estimating depths of multiple interfaces from gravity data. The strata are simulated by an aggregate of 3D rectangular prisms whose bottom depths are parameters to be estimated. In the inversion process, we have integrated geophysical constraints including the borehole information and the sharp condition described by the total variation function. The iterative residual function is also introduced to adjust the weighting of the estimated parameters so that layers of different depths have nearly equal likelihood for deviation. The inversion is processed by minimizing the Tikhonov parametric functional by the reweighted regularized conjugate gradient method. Inequality constraints are adopted to deal with the coupling of the interfaces. Synthetic tests show that such integration is conducive to restoring the multilayer depth distribution. Real data applications in Mariana confirm that the inversion method is effective in complex geologic settings in practice. We have also evaluated several issues that specifically deserve attention for obtaining satisfactory results in multilayer gravity inversion.