# Spherical Planting Inversion of GRAIL Data
## 1 Introduction
[20] used massively parallel computing techniques to achieve rapid forward modeling of gravity and gravity gradient components. However, since such calculations do not store the kernel matrix, each inversion iteration in practical applications involves several forward calculations, resulting in relatively low efficiency. In addition, there are only a few published papers because of its relatively high entry barriers in terms of hardware and software. Second, using constraint functions to accelerate gravity inversion convergence, such as optimizing recovered density distributions by a depth-weighting function [21], ensures fast convergence of the objective function by applying a preconditioning matrix [22], and adding a density constraint to optimize the upper and lower bounds of the inversion results [23; 24]. Alternatively, inversion may be constrained by using seismic [25], magnetotelluric [26], and other geophysical investigations [27]. Various filters are also used to suppress or eliminate coherent noise [28; 29]. Finally, for reducing the dimensionality of the gravity inverse problem, numerous methods such as the wavelet transform [30], footprint inversion [31], and fast Fourier transform [32] are used to reduce the memory requirements of kernel matrices and significantly improve computational efficiency. Taking wavelet transform techniques as an example, [PERSON] and [PERSON] [33] represented the kernel coefficient matrix as a sparse matrix based on the wavelet compression technique, effectively decreasing the number of nonzero matrix elements. However, it loses signal accuracy to some extent. [PERSON] [34] proposed an open, reject and fill (Open-Reject-Fill) iterative criterion for the gravity inversion method, which can grow one element at each inversion iteration. However, this method relies heavily on geological priori information. Using a systematic search algorithm, [PERSON] et al. [35] selected optimal prisms that minimized the objective function at each iteration until the termination criteria were met. Finally, this yields a three-dimensional residual density distribution. Following [PERSON]'s method [34], [PERSON] and [PERSON] [36] proposed a \"grown\" gravity inversion with the total variation regularization function given by [PERSON] et al. [37] to avoid calculating and storing kernel matrices. The algorithm grows a prism around a specific prism (called a \"seed\") to obtain the boundary information of tectonics. Nevertheless, inversion results heavily depend on the correct location, density, and seed number. Thus, [PERSON] and [PERSON] [37] achieved a robust planting inversion based on a shape-of-anomaly misfit function, an L2-norm data-misfit function, and a horizontal weighting function. However, weighting functions are challenging for inversions with various seeds, especially for complex anomalies. This paper introduces a spherical planting inversion into a spherical coordinate system. First, we develop the principles of spherical planting inversion for GRAIL data and further improve the model-misfit function by using a modified distance function. Subsequently, a robust adjustment is applied to the inversion process using a novel seed weighting function based on the covariance matrix for the given data. Then, a spherical planting inversion is performed with a composite model from [PERSON] et al. [12] to verify its correctness and efficiency. Mare Crisium is chosen as a study area, and the residual Bouguer gravity data is derived from the lunar gravity model GRGM1200A [38]. We use [PERSON]'s algorithm [39], a radially averaged energy spectrum, and [PERSON]'s (2013) algorithm [40] to determine the locations and density contrast of the seeds. The spherical planting inversion with the residual Bouguer gravity data is performed to obtain 3-D density distributions of Mare Crisium. Finally, we discuss the inverted density anomalies of the spherical planting inversion.