We develop a novel gravity inversion algorithm in spherical coordinates based on adaptive inversion mesh refinement and multiphysical parameter constraints. The inversion mesh is discretized into tesseroids (spherical prisms) to take the curvature of the Earth into account. To reduce the number of unknowns and computational cost, the inversion mesh is adaptively refined according to the spatial variation of parameters at each iteration. Wavelet compression is used to further reduce the computational requirement. Besides, to alleviate the nonuniqueness of inversion, the algorithm is capable of incorporating a priori models from other geophysical methods via cross‐gradient coupling and/or direct parameter coupling. For cross‐gradient coupling, the vector product of spatial gradients of two model parameters is included in the objective function. For direct parameter coupling, an a priori model is converted to a density model using the relation between two physical properties, which is then used as the initial and reference model in the inversion. A synthetic example is presented to demonstrate the effectiveness of our inversion method. Finally, we invert the gravity data over the Tibetan Plateau to obtain the density variations in the upper mantle, with a global S wave tomographic model used as a constraint. The inversion model shows that the mantle lithosphere beneath the Himalayan collision zone varies from west to east. Low‐density anomalies are observed beneath the southern and the northern Tibetan Plateau, which are consistent with published low‐velocity anomalies and magmatism distributions, possibly indicating two asthenospheric upwellings.