Spherical Planting Inversion of GRAIL Data

Abstract

In large-scale potential field data inversion, constructing the kernel matrix is a time-consuming problem with large memory requirements. Therefore, a spherical planting inversion of Gravity Recovery and Interior Laboratory (GRAIL) data is proposed using the L1-norm in conjunction with tesseroids. Spherical planting inversion, however, is strongly dependent on the correct seeds’ density contrast, location, and number; otherwise, it can cause mutual intrusion of anomalous sources produced by different seeds. Hence, a weighting function was introduced to limit the influence area of the seeds for yielding robust solutions; moreover, it is challenging to set customized parameters for each seed, especially for the large number of seeds used or complex gravity anomalies data. Hence, we employed the “shape-of-anomaly” data-misfit function in conjunction with a new seed weighting function to improve the spherical planting inversion. The proposed seed weighting function is constructed based on the covariance matrix for given gravity data and can avoid manually setting customized parameters for each seed. The results of synthetic tests and field data show that spherical planting inversion requires less computer memory than traditional inversion. Furthermore, the proposed seed weighting function can effectively limit the seed influence area. The result of spherical planting inversion indicates that the crustal thickness of Mare Crisium is about 0 km because the Crisium impact may have removed all crust from parts of the basin.

1. Introduction

Kernel matrices are usually generated using partial differential and integral equation methods in gravity field inversion. Partial differential equations, such as the finite-difference [1,2], finite-element [3], and finite-volume methods [4,5], are used to discretize Poisson’s equation to obtain potential fields. However, these methods are difficult to apply to gravity inversion because higher-order derivatives involve higher-order differences or the fast Fourier transform. Comparatively, integral equations utilize the superposition principle to calculate the gravity contribution of complex geological bodies.

The observation surface of gravity data in small areas can be approximated as a plane; in the Cartesian coordinate system, the potential field calculation in the spatial domain is usually performed using prisms or general polyhedral cells [6,7,8,9]. In large regional and even global gravity problems, the observation surface is generally expressed and processed in spherical coordinates because of the need to consider the observation surface’s curvature [10]. As a result, interpretation models are usually divided into tesseroids in a spherical coordinate system using the parameters radius, latitude, and longitude. The gravitational effect generated by each unit is calculated and then summed to obtain global gravity/tensor gravity data.

As a result of the CHAMP (Challenging Mini-Satellite Payload), GRACE (Gravity Recovery and Climate Experiment), GOCE (Gravity Field and Steady-State Ocean Circulation Explorer), and GRACE Follow-On satellite launches [11], it is now possible to obtain large-scale high-precision gravity data economically.

Inversion of the gravity field is a powerful tool for understanding the density distribution inside planets. Many gravity inversions have been proposed in recent years. Liang et al. applied tesseroid units to the lunar GRAIL gravity inversion and obtained the 3D density structure of the Moon [12]. Uieda et al. [13] analyzed the Moho surface depth of the South American continent using tesseroid grids for a fast nonlinear inversion of GOCO5S model data in a spherical coordinate system. Zhang et al. [14] developed an efficient algorithm for gravity gradient full-tensor data, which was used to study the density distribution of the Mare Smythii mascon on the Moon. Zhao et al. [15] performed a 3D spherical gravity inversion with tesseroid units for the GL1500E model data to obtain the 3D density distribution of the Imbrium and Serenitatis mascon basins.

ETOPO1 is a one-arc-minute global relief model of the Earth’s surface that integrates land topography and ocean bathymetry provided by the National Oceanic and Atmospheric Administration [16]. Take a spherical gravity inversion with the ETOPO1 model as an example. When a 1- by 1- minute latitude-longitude grid is used to cover the globe, the number of observation points is 180 × 60 × 360 × 60 = 2.3328 × 10⁸. The number of kernel matrix elements for potential inversion is the product of the number of tesseroid units and observation points. The memory requirements of the kernel matrix, whether storing double-precision or single-precision data, are far more than the physical memory size of existing mainstream workstations.

Constructing kernel matrices is a time-consuming problem with high memory requirements for potential inversion, especially for large-potential data [17]. Fast and accurate three-dimensional gravity inversion of large-scale, multi-scale and massive data is currently performed by improving the computational performance of equipment, reducing data dimensionality, and speeding up the convergence of the inversion iterations.

First, high-performance computation relies on multi-core CPUs, clustered computer systems, and GPUs. Then, massively parallel acceleration studies are performed to reduce time-consuming forward modeling and inversion [18]. For example, Moorkamp et al. [19] and Chen et al. [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, Li and Oldenburg [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. René [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, Camacho 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 René’s method [34], Uieda and Barbosa [36] proposed a “grown” gravity inversion with the total variation regularization function given by Silva Dias 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, Uieda and Valeria [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 Liang 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 Parker–Oldenburg’s algorithm [39], a radially averaged energy spectrum, and Wieczorek’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.

2. Materials and Methods

2.1. Spherical Planting Inversion

Following the ideas of Uieda and Barbosa [36,37], to avoid solving an extensive system of linear equations constructed by density imaging, the spherical planting inversion of GRAIL data is proposed using an L1-norm-based systematic search algorithm to reduce memory requirements and facilitate fast convergence by replacing iterative inversion with cumulative summation analysis. For this reason, a residual vector $r$ between the observation data $d$ and the prediction data $g$ is constructed by the L1-norm based on Equation (1).

(1)$r=d-g$

The data-misfit function is defined (based on the L2-norm) by least squares estimation, which can be expressed as

(2)${\varphi }_d={\Vert r\Vert }_2={\left({\displaystyle \sum _{i=1}^{N_d}{\left(g_i-d_i\right)}^2}\right)}^{1/2}$

where the subscript i corresponds to the i-th observation point. $N_d=n_x\times n_y$. $n_x$ $n_y$, and $n_r$ represent the section number of the interpretation model along the x-, y- and radial axes, respectively.

Compared to the L2-norm, the L1-norm is not sensitive to data with significant errors so it can provide a robust estimation strategy.

(3)${\varphi }_d={\Vert r\Vert }_1={\displaystyle \sum _{i=1}^{N_d}\left|g_i-d_i\right|}$

Conventional inversion applies regularization terms to reduce the instability and non-uniqueness of the solutions. Many inverse strategies have been proposed for recovering compact anomalous sources, such as maximizing their compactness [41], concentrating them along one or more axes [42,43], and minimizing their area/volume [44]. For example, René [34] proposed a systematic search inversion of gravity data using the “shape-of-anomaly” function to reconstruct 2D compact bodies (i.e., with no hollows inside). However, this inversion algorithm only uses a single sign/density contrast to recover density distributions. Compared with René’s (1986) method, Uieda and Barbosa interpret multiple geologic sources with density contrasts of different signs [36].

Consider that recovering a density distribution from GRAIL data is an ill-conditioned problem, it can be converted into a well-conditioned problem by adopting the following constraints:

(1) A compact scheme ensures that there are no voids in the recovered density distributions, which are obtained from spherical gravity inversion with a seed [36,41,42,43].

(2) An optimal prism is selected one by one from a neighboring tesseroids list, which is derived from a seed with density contrast ${\rho }_s$, and shares faces with the seed or its adjacent tesseroids. Then, the optimal prism’s density value ${\rho }_j$ is set to the same as the seed.

(3) For other uninvolved tesseroids, their density value equals zero. Based on the above constraint strategies, the objective function of the spherical planting inversion of GRAIL data is given by

(4)$P^a\left(m\right)={\varphi }_d+\mu {\varphi }_m\to \min$

where, the regularization parameter $\mu$ balances the trade-off between the data-misfit function ${\varphi }_d$ and the model-misfit function ${\varphi }_m$, and is defined in the model space that imposes geologic constraints on the recovered density distribution.

(5)${\varphi }_m={\displaystyle \sum _{k=1}^{N_s}\left({\displaystyle \sum _{j=1}^{N_m}\frac{{\rho }_j}{{\rho }_j+\epsilon }\varsigma \left(j,k\right)l_{j,k}^{\beta }}\right)}$

where, ${\rho }_j$ is the j-th element of the density vector m, whose size is N_m, and N_s is the number of seeds. $\epsilon$ is a small positive count and will avoid a singularity error when ${\rho }_j=0$. The sign function $\varsigma \left(j,k\right)$ returns 1 when the j-th tesseroid is derived from the k-th seed, as shown in Figure 1; otherwise, it returns 0. When N_s = 1, Formula (5) is similar to Formula (10) in Uieda and Barbosa’s paper [36]. $l_{j,k}^{\beta }$ is the distance from the j-th tesseroid’s centroid to the k-th seed’s centroid.

Regardless of how a tesseroid’s length and width vary with its position, the modified version $l_{j,k}^{\beta }$ is given by

(6)$l_{j,k}^{\beta }=\sqrt{{\left(j_x-k_x\right)}^2+{\left(j_y-k_y\right)}^2+{\left(j_r-k_r\right)}^2}$

where $\left(j_x,j_y,j_r\right)$ and $\left(k_x,k_y,k_r\right)$ stand for the order of the j-th tesseroid and the k-th seed along the longitude, latitude, and radius, respectively. Here, $j=j_x+(j_y-1)\times n_x+(y_r-1)\times n_x\times n_y$.

The solution of the spherical planting inversion of GRAIL data is recovered through a growth procedure divided into the following four steps, as shown in Figure 1.

1. Initially, this algorithm requires a set of N_s seeds, each of which is a tesseroid of the geophysical model. Seed locations can be obtained by Euler deconvolution or by manual selection, and the user can specify the density contrast of seeds. Meanwhile, the density contrast of all other tesseroids equals zero.

2. Secondly, for each seed, a list is constructed of neighboring tesseroids that share a face with the seed or neighboring ones. In constraint 3, only seeds will have a specific density contrast, and all other tesseroids’ density contrasts are set to zero. Then the initial total objective function $P^{\left(0\right)}$ is given by

(7)$P^{\left(0\right)}={\varphi }_d^{\left(0\right)}+\mu {\varphi }_m^{\left(0\right)}$

where, the initial total data-misfit function is given by

(8)${\varphi }_d^{\left(0\right)}=d-{\displaystyle {\sum }_{i=1}^{N_s}{\rho }_s}{a}_s$

Here, ${\rho }_s$ is the density-contrast value of the s-th seed, and the vector $a_s$ consists of kernel functions which are gravitational effects at all observation points due to the s-th seed. Due to there being no neighboring tesseroids in the initial step, ${\varphi }_m^{\left(0\right)}=0$.

3. An iteration involves Ns accretions in a spherical planting inversion. This algorithm finds a tesseroid from its neighboring tesseroid list for every accretion by minimizing the objective function and guaranteeing convergence of the following two conditions:

(9)$\left|{\varphi }_d^{\left(new\right)}-{\varphi }_d^{\left(old\right)}\right|/{\varphi }_d^{\left(old\right)}\ge \delta$

where the total data-misfit function ${\varphi }_d^{\left(new\right)}$ is calculated for the chosen neighboring tesseroids in the current accretion and all tesseroids whose density contrast is not equal to zero, ${\varphi }_d^{\left(old\right)}$ is the total data-misfit function at the previous accretion. The small positive scalar $\delta$ ranges from 10⁻² to 10⁻⁶.

(10)$\left|P_{}^{\left(new\right)}-P_{}^{\left(old\right)}\right|/P_{}^{\left(old\right)}\ge \tau$

where $P^{\left(new\right)}$ and $P^{\left(old\right)}$ are total global functions of current accretion and precious accretion, respectively.

Ensuring the chosen tesseroid satisfies Equations (9) and (10), and treating it as the optimal tesseroid in the current iteration, then, the algorithm changes its density contrast from zero to the seed’s and removes it from the seed’s neighboring tesseroids list to guarantee constraint 3; otherwise, the algorithm will traverse over other seeds.

4. In any accretion at the previous iteration, if the convergence conditions (Equations (9) and (10)) are met, the growth process continues. Rather, the spherical planting inversion should be terminated, and a 3D density distribution should be obtained.

To avoid the recovery density distribution of spherical planting inversion being strongly influenced by the parameters of seeds, Uieda et al. [37] proposed a “shape-of-anomaly” data misfit that measures the difference in shape between the observed and predicted data and is insensitive to differences in amplitude. By introducing “shape-of-anomaly” data misfit, the data-misfit function is defined as:

(11)${\varphi }_d=\sqrt{{\displaystyle {\sum }_{i=1}^{N_d}{\left(a{g}_i-d_i\right)}^2}}$

where $\alpha$ is the scaling factor; for a given predicted data $g$, $a$ can be calculated by transforming the above equation:

(12)$a={{\displaystyle {\sum }_{i=1}^{N_d}\left(g_i{d}_i\right)}/{\displaystyle {\sum }_{i=1}^{N_d}\left(d_i\right)}}^2$

In the spherical planting inversion, the reconstructed density distributions are determined by the seeds’ parameters. As the reconstruction model’s size increases, the number of neighboring tesseroids gradually increases, and the distance between seeds and neighboring tesseroids increases. This causes neighboring tesseroids in a particular direction always to be determined as optimal tesseroids, eventually leading to mutual intrusion between adjacent anomalous sources. This issue is addressed by the weighting function, which limits the seeds’ influence area and makes them ignore far abnormal sources.

(13)$w_h=\exp {\left(\frac{{\left(x_s-x\right)}^2+{\left(y_s-y\right)}^2}{{\left(r_s\right)}^2}\right)}^{-{\beta }_h}$

where (x_s, y_s) and (x, y) are the coordinates of the s-th seed and observation point, respectively, and r_s is the influence radius of the s-th source. For the inversion of gravity data, the horizontal decay factor ${\beta }_h$ ranges from 1 to 3. Then, the data-misfit function (11) can be rewritten as:

(14)${\varphi }_d=\sqrt{{\left(a{w}_{h}g-d\right)}^{\mathrm{T}}\left(ag-d\right)}$

where T is the transpose operator.

However, it is challenging to choose the right values for Equation (13), particularly when the quantity of seeds is large or the gravity anomaly is complicated. To overcome this difficulty, by referring to the data covalence matrix $W_d$ in the density image, we can define a novel seed weighting function:

(15)${\overline{w}}_h=diag\left(W_d\right)=diag\left(\left\{1/{\sigma }_1,\cdots ,1/{\sigma }_n\right\}\right)$

where ${\sigma }_i$ is the error standard deviation associated with the i-th observation data.

The process of the spherical planting inversion of GRAIL data is detailed in a flow chart, shown in Figure 2.

2.2. Model Studies

2.2.1. Verification of the Validity of the Novel Seed Weighting Function ${\overline{w}}_h$

First, a simple synthetic model with an isolated anomalous source was proposed by Zhang et al. [14] and is used as an example to compare the reconstruction effect of the two weighting functions mentioned in Equations (13) and (15), respectively. The isolated anomalous source is in the blue frame as in Figure 3 and with 0.5 g/cm³ density contrast. Then, a rectangular survey grid, with 40 points, with a 0.5° grid interval along the longitude and latitude directions, is employed to obtain gravity data. A seed with locations 35°, 35°, and 1673 km along latitude, longitude, and radius directions, respectively, and with a density contrast of 0.5 g/cm³, is employed for performing the spherical gravity inversion for gravity data, which is contaminated by 3% Gaussian noise of the gravity data.

Using a seed with locations at 35°, 35°, and 1673 Km along latitude, longitude, and radius directions, respectively, and with a density contrast of 0.5 g/cm³, we perform spherical gravity inversion for gravity data, which is contaminated by 3% Gaussian noise in the forward result.

Regarding relatively simple anomalous sources, it is feasible to set a horizontal weighting function for multiple seeds, which can effectively limit the seed’s influence area and avoid the intrusion of the estimation grid of adjacent seeds into each other. Despite this, as far as the characteristics of the horizontal weighting function $w_h$ are concerned, its influence range is a circle. Consequently, it is too expensive to artificially set the horizontal weighting function for multiple seeds, especially for complex gravity data. The recovered density distributions in Figure 3a are far beyond the actual model in the depth direction, which is caused by the gravitational contribution of the tesseroid becoming smaller with the increase in depth. Meanwhile, the corresponding data misfit function becomes smaller, so the optimal tesseroid will gradually increase along the depth direction. Eventually, the reconstructed density model will be far beyond the actual model in the depth direction. For this reason, drawing on the role of the data weighting function W_d in density imaging, we construct a novel seed weighting function. Compared to the recovered density distribution derived from $w_h$ as shown in Figure 3a, the spherical planting inversion with the horizontal weighting function yields a more compact result in Figure 3b.

2.2.2. Verification of the Validity of Spherical Planting Inversion

As shown in Figure 4, we verify the proposed algorithm using the composite model following Zhang et al. [14] with two tesseroids at the same depth for easy comparison. One model occupied the volume from 32° to 34° longitude, 34° to 36° north latitude, and 1678 to 1718 km in depth, with a density contrast of 0.2 g/cm³. The other model had a density contrast of 0.7 g/cm³, and its geological setting ranged from 36° to 38° longitude, 34° to 36° latitude, and 1658 to 1698 km in depth. A total of 62,500 data points is calculated 10 km above the lunar surface with an interval of 0.8° × 0.8° along latitude and longitude. The field source space is divided into 20 × 20 × 10 = 4000 tesseroids, and each tesseroid is 0.8° × 0.8° × 10 km along latitude and longitude, and radius direction, respectively. Figure 5 shows the data contaminated with pseudo-random Gaussian noise with zero means and five mGal standard deviations. Then, we invert the contaminated data with two seeds; one seed is located at [33°, 35°, 1698 km] with 0.2 g/cm³ density contrast, and the coordinates of the other seed were 35°, 37° and 1678 km along latitude, longitude and radius direction, respectively, with a density contrast 0.7 g/cm³.

It is difficult for the Occam inversion to recover depth resolution from surface gravity data because there is no depth resolution in the data [21,45,46]. The interfaces of the two models in the spherical Occam inversion results are ambiguous, even though the depth weighting matrix is used. In contrast, the density slices of the spherical planting inversion are consistent with those of the forward model in Figure 6, which verifies the correctness of the algorithm presented in this paper.

3. Results

The discovery of lunar mascons was one of the most significant results of early gravity field studies in lunar exploration [47,48]. Moreover, Mare Crisium is one of the largest mascons on the Moon [49] and is an impact-formed poly-cyclic basin [50,51] that formed during the Nectarian Period [52,53], located between 10° and 30° N and 50° and 70° E, in the northeastern part of the lunar front (centered near 17.0° N and 58.8° E), with a basin diameter of 555 km and covering an area of about 17,600 km². The main ring of the crater rim is about 740 km in diameter and is one of the deepest basins on the lunar surface (4.57 km). Therefore, we selected Mare Crisium as an example to elucidate the mechanism of spherical planting inversion.

3.1. Data

The gravity model GRGM1200A, applied in Bouguer correction by NASA, expanded to degree and order 1200 with sensitivity down to <5 km resolution from the GRAIL mission, is used in this paper [54]. Deliberating the very high correlation between gravity and topography at L = 7~700 [55], following Deng et al. [56], the spatial resolution of the gravity field is related to the spherical harmonic coefficient and defined by the half wavelength in spherical inversion [57]. When the grid size of the interpretation model is 0.8° × 0.8° on a spherical earth, according to the relation between half wavelength lambda and the suitable degree N of the gravity model in Liang et al. [12], the suitable degree N = 450 [58]. Therefore, in Figure 7, the residual Bouguer gravity is expanded between spherical degrees 7 and 450 [57] to suppress the influence of noise from the high degrees and filter out signals from the deep interior (such as the crust–mantle boundary or Moho) from the low degrees for highlighting lunar mascons [59,60,61]. The topography of the Crisium mascon is derived from the spherical harmonic topography model MoonTopo2600p, as shown in Figure 8. Then, 20 × 20 = 400 observation data points are calculated at 10 km height relative to the reference radius of 1738 km using the gravity field model mentioned above, as shown in Figure 9.

In Figure 8, Mare Crisium is surrounded by an upland massif [53], with much of its interior and a small part of its marginal upland filled with basaltic rocks and a ring of folded ridges on its outer boundary [62,63]. A folded ridge is a low, linear, or arcuate broad ridge, usually tens of kilometers long. It is placed symmetrically/asymmetrically on an arch close to the circular Maria edge [64,65].

3.2. Methodology

A new Bouguer gravity data covering 4°–34° N and 50°–70° E is gridded at an interval of 2 km with a Mercator projection centered at (180 E, 0 N) to determine the depth and density contrast of the Mare Crisium mascons. This area is designed to eliminate the effect of distorted shapes due to the Mercator projection. Assuming that the density contrast is 0.65 g/cm³, the convergence criterion is 0.025 km, the mean crustal reference depth is 30 km [40], the truncation window data length is 2.5%, and the smaller and greater cut-off frequencies are 1/100 and 1/83.3, Parker–Oldenburg’s algorithm [39] was employed to invert the residual Bouguer gravity anomaly to reconstruct the crust thickness in this study area. Figure 10 shows that the crust thickness ranges from 0 to 45 km, and the topography variation outlines the impact crater basins.

The power spectrum method performs regional and residual separations based on a fast Fourier transform for outlining subsurface structures. As an extension of the power spectrum method, the radially averaged energy spectrum yields a typical energy spectrum plot containing three parts corresponding to deep source, superficial source, and noise contribution. Hence, the power spectrum method can also be employed to determine the depths of shallow and deeper structures [66,67,68,69].

By following [66,70,71], the radially averaged power spectrum according to frequencies was calculated and is displayed with blue dots in Figure 11. Then, the slopes of the portions can be employed to estimate the depths of anomalous sources. The depth (z) for each source corresponding to each segment was calculated by presenting the slope of this part:

$Z=\frac{-D}{4\pi }$

where the slope of the dashed line D is the ratio between the logarithmic power spectrum and the wavenumber.

The radially averaged power spectrum curve can be interpreted as two linear slope segments and has been applied to estimate the depth of the gravity interface at 0.4 km, representing the average depth to the surface of the mascon (Figure 11).

In order to further confirm the depth and density of the mascon in this study area, the spherical harmonic density model density_no_mare_n3000_f3050_719.sh is proposed by Wieczorek [72], and the spherical harmonic topography model MoonTopo2600p [73] and the spherical harmonic gravity model GRGM900C [55] are employed to calculate the crustal density, topography, and gravity anomaly, respectively. The crust is assumed to have an average thickness of 35 km, a porosity of 0.12, and a density of 3.22 g/cm³. By following [40] and [74], we used homogeneous density to minimize the correlation between surface topography and Bouguer gravity and determine the crustal thickness. From the Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12, the crustal thickness in this study area is described as 5 to 40 km. When combined with the topographic map, we can find that the crustal thickness in the middle of this study area is 0.

3.3. Spherical Planting Inversion of Mare Crisium

To achieve the spherical planting inversion of residual Bouguer gravity at the Crisium mascon, the geophysical interpretation model covers an area of 9°–29°, 50°–70° and 1638–1738 km in latitude and longitude as well as radial directions. The field source space is divided into 20 × 20 × 10 = 4000 cells, and each cell is 0.8° × 0.8° × 10 km along latitude and longitude, and radius directions, respectively.

Because of the low gravity at each corner and the similar linear background field in Figure 9, following Uieda and Barbosa (2012) [36], only positive gravity anomalies of interest associated with the mascon are targeted for exploration. Therefore, in this paper, ten seeds (asterisk marker in Figure 13) were employed to invert only target sources (black contours in Figure 14) of the Mare Crisium mascon.

In spherical planting inversion, seed parameters are generally determined empirically or based on drilling information. Considering that previous studies [40,75,76] used 2.580 ± 0.170 g/cm³ as the mean crustal density and 3.22 g/cm³ as the mean mantle density, and 2.50 g/cm³ was employed to Bouguer correction for GRGM1200A, the seed density is now set to 0.45 g/cm³. Because the Mare Crisium mascon is a mantle uplift product, according to the crustal thickness of Mare Crisium in Figure 12, the depth of the seeds is set at 30 km.

For larger or more complex anomalous sources, it is challenging to set individual horizontal weighting functions for each seed, as in Figure 10 and Figure 11 in Uieda and Barbosa (2012) [36], where the recovered density distribution in the right branch is well beyond the limit areas of the horizontal weighting function. Further considering the gravity anomaly characteristics (Figure 9), ten seeds with yellow asterisk markers are now arranged in the positive and negative transition regions of the Bouguer gravity anomaly (Figure 13). Then, the spherical planting inversion is performed using a single seed horizontal weighting function with $\delta$ = 0.05 and $\tau$ = 5.0 × 10⁻⁴.

Compared with Figure 14, the difference between the predicted anomalies of the spherical planting inversion and the exploration target anomalies is tiny, and its norm is 1.2, indicating the present algorithm’s stability. Figure 13 and Figure 14 show that the novel seed weighting function fits the predicted data range, indicating that it can effectively constrain the spherical planting inversion [40,77].

4. Discussion and Conclusions

Since the tesseroid’s length and width vary with its position, we revise the distance calculation formula to avoid preferentially selecting tesseroids closer to their corresponding seeds in the spherical planting inversion.

Because the spherical planting inversion depends heavily on the user’s seed parameters, it can incorporate geological priori knowledge to constrain reconstructed density distributions [36,37]. It is like the initial model in iterative inversion. The initial model is not a controlling factor in traditional iterative inversion, and many inversion reports use random initial models. However, the location of seeds in the spherical planting inversion largely determines the spatial spread of the reconstructed density distribution.

Since the gravity field originates from the volume effect of various subsurface anomalous sources [78], incorrect seed densities lead to the recovered density model deviating from the actual one. If the seed density is smaller than the theoretical value, the reconstructed density distribution becomes inflated compared to the actual model. Conversely, the recovered density distribution becomes more compact than the actual one. Because the gravity response rapidly decays with tesseroid depth [79], we found in repeated experiments that, when the threshold is low, the spherical planting inversion tends to select tesseroids near the seeds. As a result, the predicted gravity anomaly obtained from the inversion is similar to the observed one. However, the total mass of the recovered density distribution is significantly larger than the actual model. For this reason, in this paper, we analyze the density/geology studies in this study area and give a more reasonable residual density value for setting seeds.

When a single seed or a few seeds (N_s < 4) are used, it was found that the recovered density distribution grew in certain directions and eventually formed a trunk-like/tentacle-like structure, as shown in Figure 10 and Figure 11 in Uieda et al. [36]. In contrast, a compact geological structure is formed in the middle of the recovered density distribution when numerous seeds are used, as shown in Figure 15 and Figure 16, with its top and bottom appearing as scattered tesseroids because there is no tight constraint between seeds according to Equation (4).

Although tentacles can help optimize seed settings in the subsequent inversion iteration, choosing individual parameters for each seed is difficult, especially when processing large or complex gravity anomalies. For this reason, we propose a novel horizontal weighting function based on W_d, as illustrated in Figure 13. Figure 15 and Figure 16 show that the novel horizontal weighting function effectively constructs the subsurface geological structure. It also effectively avoids tentacles in the recovered density distribution.

The algorithm involves calculating the kernel matrix at most once so that the spherical planting inversion of gravity data is faster and has lower memory requirements/consumption than traditional density imaging. The inversion algorithm in this paper is a constrained system search algorithm that achieves fast and efficient spherical planting inversion with few memory requirements. There may be a significant increase in searches in one iteration if the number of seeds is too large.

Correct seed densities are crucial for a spherical planting inversion. Conversely, the recovered density distribution is easily inflated or shrunken compared to the actual geological structure. When the seed’s density is smaller than the theoretical value, it is easy to cause the recovered density distribution range to be much more extensive than the actual one. This results in dramatically increasing inversion iterations and reducing inversion efficiency. Therefore, when determining the seed’s density, it is necessary to scrutinize the density research history in the study area.

Since the smoothness effect originates from the L2 regulation term [80], the recovered density distributions are smooth everywhere, making it difficult to distinguish the boundaries of geological structures. Also, their density value is far from the theoretical value due to the volume effect of the gravity–magnetic field. There is no depth resolution in surface gravity data. It is challenging to use Occam-like inversion to recover depth information [21,45,46]. Nevertheless, the spherical planting inversion employed in this paper can preserve depth resolution.

Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 illustrate that the depth of the top interface of the Mare Crisium basin is about 4–6 km. Figure 15 and Figure 16 show that the mascon has spread from the surface, showing that the crustal thickness of Mare Crisium is about 0 km because the Crisium impact may have removed all crust from parts of the basin [40,77].

Author Contributions

G.L.: Conceptualization, Supervision, Writing—review and editing; D.Z.: Software, Validation, Writing—original draft, review and editing; S.C.: Formal analysis, Methodology, Software, Writing—review and editing; Y.D.: Software, Validation; G.X.: Software, Investigation; Y.L.: Software, Investigation; Z.Z.: Resources; P.C.: Validation, Visualization. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is unavailable.

Acknowledgments

We appreciate the GRAIL mission team distributing the gravity models GL0900C and GRGM1200A. The program in SHTools, which is freely accessible at https://shtools.github.io/SHTOOLS/ (accessed on 4 March 2023), is used to compute gravity anomalies. We appreciate M. Wieczorek for sharing the spherical harmonic model of the shape of Earth’s Moon MoonTopo2600p and the software CTplanet for estimating crustal thickness, which is accessible for free at https://github.com/MarkWieczorek/ctplanet (accessed on 4 March 2023). The authors would like to thank D. Gomez-Ortiz and BNP Agarwal for making Parker–Oldenburg’s algorithm available. Valuable comments by two anonymous reviewers are gratefully acknowledged. Moreover, we would like to thank Tiffany Liu and the Applied Sciences editorial team for their persistent efforts.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Four stages of the spherical planting algorithm: (a) Initialization of seeds (in cyan) at the first stage; (b) Initialization of neighboring tesseroids (in green) of seeds at the second stage; (c) Finding the optimal tesseroid (in red) at first accretion, and (d) updating neighboring tesseroids of the first seed at the third stage; (e) Finding the optimal tesseroid at second accretion, and (f) updating neighboring tesseroids of the second seed at the fourth stage.

Figure 2. Flow chart for the spherical planting inversion.

Figure 3. Recovered density distribution of spherical planting inversion of gz data. (a) With the horizontal weighting function [Forumla omitted. See PDF.]. (b) With the novel seed weighting function [Forumla omitted. See PDF.]. The weight distributions as shown in the z = 120 km plane. [Forumla omitted. See PDF.] = 50, [Forumla omitted. See PDF.] = 1.2, [Forumla omitted. See PDF.] = 0.05, and [Forumla omitted. See PDF.] = 5.0 × 10−4.

Figure 4. Composite model: (a) A horizontal slice of the composite model with a radial length of 1698 km. (b) A vertical profile of the composite model with a latitude of 35°. [Forumla omitted. See PDF.] = 2.0 × 10−5 and [Forumla omitted. See PDF.] = 1.0 × 10−6.

Figure 5. Gravity maps: (a) Forward modeling result. (b) Prediction data obtained by spherical planting inversion. (c) Residual data between forward modeling result (a) and prediction data (b). [Forumla omitted. See PDF.] = 0.05, and [Forumla omitted. See PDF.] = 5.0 × 10−4.

Figure 6. Density-contrast distribution maps. (a) Horizontal slice of the reconstruction model along with a radial length of 1698 km, and (b) vertical profile of the composite model along a latitude of 35°.

Figure 7. A filtered Bouguer map (l = 7−450) highlights the lunar mascons.

Figure 8. The topography of the Crisium mascon.

Figure 9. The residual Bouguer gravity data.

Figure 10. The crustal thickness was obtained from the residual Bouguer gravity (Figure 9) using Parker–Oldenburg’s method.

Figure 11. Radially averaged power spectrum.

Figure 12. The crustal thickness of Mare Crisium.

Figure 13. New horizontal weighting function applied to field data.

Figure 14. Gravity maps. (a) Observed data (color−scale map) and predicted data (black contour lines). The observation data is derived from GRGM1200 and is also shown in Figure 9. Due to the low gravity at each corner and the similar linear background field in Figure 9, following Uieda and Barbosa (2012) [36], only positive gravity anomalies of interest associated with the mascon are targeted for exploration. Predicted data is obtained using a spherical planting inversion with ten seeds (Figure 13). (b) Targeted data (color-scale map) and predicted Bouguer gravity anomalies (black contour lines).

Figure 15. Vertical profiles of the recovered density distribution along (a) longitude 56°, (b) longitude 64°, (c) latitude 15°, and (d) latitude 25° at Mare Crisium are obtained by spherical planting inversion.

Figure 16. Relative density distribution for different depths, (a–d) are 10 km, 25 km, 40 km, and 50 km, respectively, at Mare Crisium area obtained by spherical planting inversion. Figure 10 shows that the depth of the Mare Crisium basin is about 4–6 km. Figure 15 and Figure 16 show that the mascon has spread from the surface, indicating that the crustal thickness of Mare Crisium is about 0 km because the Crisium impact may have removed all crust from parts of the basin [40,77].

References

1. Farquharson, C.G.; Mosher, C.R.W. Three-dimensional modelling of gravity data using finite differences. J. Appl. Geophys.; 2009; 68, pp. 417-422. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2009.03.007]

2. Jahandari, H.; Bihlo, A.; Donzelli, F. Forward modelling of gravity data on unstructured grids using an adaptive mimetic finite-difference method. J. Appl. Geophys.; 2021; 190, 104340. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2021.104340]

3. Martyshko, P.; Ladovskii, I.; Byzov, D.; Tsidaev, A. Gravity Data Inversion with Method of Local Corrections for Finite Elements Models. Geosciences; 2018; 8, 373. [DOI: https://dx.doi.org/10.3390/geosciences8100373]

4. Lelièvre, P.G.; Farquharson, C.G.; Hurich, C.A. Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration. Geophysics; 2012; 77, pp. K1-K15. [DOI: https://dx.doi.org/10.1190/geo2011-0154.1]

5. Jahandari, H.; Farquharson, C.G. Forward modeling of gravity data using finite-volume and finite-element methods on unstructured grids. Geophysics; 2013; 78, pp. G69-G80. [DOI: https://dx.doi.org/10.1190/geo2012-0246.1]

6. Ren, Z.; Chen, C.; Pan, K.; Kalscheuer, T.; Maurer, H.; Tang, J. Gravity Anomalies of Arbitrary 3D Polyhedral Bodies with Horizontal and Vertical Mass Contrasts. Surv. Geophys.; 2017; 38, pp. 479-502. [DOI: https://dx.doi.org/10.1007/s10712-016-9395-x]

7. Zhang, Y.; Chen, C. Forward calculation of gravity and its gradient using polyhedral representation of density interfaces: An application of spherical or ellipsoidal topographic gravity effect. J. Geod.; 2018; 92, pp. 205-218. [DOI: https://dx.doi.org/10.1007/s00190-017-1057-3]

8. Ren, Z.; Chen, C.; Zhong, Y.; Chen, H.; Kalscheuer, T.; Maurer, H.; Tang, J.; Hu, X. Recursive Analytical Formulae of Gravitational Fields and Gradient Tensors for Polyhedral Bodies with Polynomial Density Contrasts of Arbitrary Non-negative Integer Orders. Surv. Geophys.; 2020; 41, pp. 695-722. [DOI: https://dx.doi.org/10.1007/s10712-020-09587-4]

9. Zhang, J.; Jiang, L. Analytical expressions for the gravitational vector field of a 3-D rectangular prism with density varying as an arbitrary-order polynomial function. Geophys. J. Int.; 2017; 210, pp. 1176-1190. [DOI: https://dx.doi.org/10.1093/gji/ggx230]

10. Ren, Z.; Tang, J.; Kalscheuer, T.; Maurer, H. Fast 3-D large-scale gravity and magnetic modeling using unstructured grids and an adaptive multilevel fast multipole method. J. Geophys. Res. Solid Earth; 2017; 122, pp. 79-109. [DOI: https://dx.doi.org/10.1002/2016JB012987]

11. Zheng, W.; Xu, H. Progress in satellite gravity recovery from implemented CHAMP, GRACE and GOCE and future GRACE Follow-On missions. Geod. Geodyn.; 2015; 6, pp. 241-247. [DOI: https://dx.doi.org/10.1016/j.geog.2015.05.005]

12. Liang, Q.; Chen, C.; Li, Y. 3-D inversion of gravity data in spherical coordinates with application to the GRAIL data. J. Geophys. Res. Planets; 2014; 119, pp. 1359-1373. [DOI: https://dx.doi.org/10.1002/2014JE004626]

13. Uieda, L.; Barbosa, V.C.F. Fast nonlinear gravity inversion in spherical coordinates with application to the South American Moho. Geophys. J. Int.; 2016; 208, pp. 162-176. [DOI: https://dx.doi.org/10.1093/gji/ggw390]

14. Zhang, Y.; Wu, Y.; Yan, J.; Wang, H.; Rodriguez, J.A.P.; Qiu, Y. 3D inversion of full gravity gradient tensor data in spherical coordinate system using local north-oriented frame. Earth Planets Space; 2018; 70, pp. 1-23. [DOI: https://dx.doi.org/10.1186/s40623-018-0825-5]

15. Zhao, G.; Chen, B.; Uieda, L.; Liu, J.; Kaban, M.K.; Chen, L.; Guo, R. Efficient 3-D Large-Scale Forward Modeling and Inversion of Gravitational Fields in Spherical Coordinates With Application to Lunar Mascons. J. Geophys. Res. Solid Earth; 2019; 124, pp. 4157-4173. [DOI: https://dx.doi.org/10.1029/2019JB017691]

16. Amante, C.; Eakins, B.W. ETOPO1 1 arc-minute global relief model: Procedures, data sources and analysis. NOAA technical memorandum NESDIS NGDC-24. Natl. Geophys. Data Cent. NOAA; 2009; 10, V5C8276M.

17. Čuma, M.; Wilson, G.A.; Zhdanov, M.S. Large-scale 3D inversion of potential field data. Geophys. Prospect.; 2012; 60, pp. 1186-1199. [DOI: https://dx.doi.org/10.1111/j.1365-2478.2011.01052.x]

18. Čuma, M.; Zhdanov, M.S. Massively parallel regularized 3D inversion of potential fields on CPUs and GPUs. Comput. Geosci.; 2014; 62, pp. 80-87. [DOI: https://dx.doi.org/10.1016/j.cageo.2013.10.004]

19. Moorkamp, M.; Jegen, M.; Roberts, A.; Hobbs, R. Massively parallel forward modeling of scalar and tensor gravimetry data. Comput. Geosci.; 2010; 36, pp. 680-686. [DOI: https://dx.doi.org/10.1016/j.cageo.2009.09.018]

20. Zhao-Xi, C.; Xiao-Hong, M.; Liang-Hui, G.; Guo-Feng, L. Three-dimensional fast forward modeling and the inversion strategyfor large scale gravity and gravimetry data based on GPU. Chin. J. Geophys.; 2012; 55, pp. 4069-4077. [DOI: https://dx.doi.org/10.6038/j.issn.0001-5733.2012.12.019]

21. Li, Y.; Oldenburg, D. 3-D inversion of gravity data. Geophysics; 1998; 63, pp. 109-119. [DOI: https://dx.doi.org/10.1190/1.1444302]

22. Pilkington, M. 3-D magnetic imaging using conjugate gradients. Geophysics; 1997; 62, pp. 1132-1142. [DOI: https://dx.doi.org/10.1190/1.1444214]

23. Commer, M.; Newman, G.A. New advances in three-dimensional controlled-source electromagnetic inversion. Geophys. J. Int.; 2008; 172, pp. 513-535. [DOI: https://dx.doi.org/10.1111/j.1365-246X.2007.03663.x]

24. Abubakar, A.; Habashy, T.M.; Li, M.; Liu, J. Inversion algorithms for large-scale geophysical electromagnetic measurements. Inverse Probl.; 2009; 25, pp. 1541-1548. [DOI: https://dx.doi.org/10.1088/0266-5611/25/12/123012]

25. Sampietro, D.; Capponi, M. Seismic Constrained Gravity Inversion: A Reliable Tool to Improve Geophysical Models Away from Seismic Information. Geosciences; 2021; 11, 467. [DOI: https://dx.doi.org/10.3390/geosciences11110467]

26. Zhang, R.; Li, T.; Zhou, S.; Deng, X. Joint MT and Gravity Inversion Using Structural Constraints: A Case Study from the Linjiang Copper Mining Area, Jilin, China. Minerals; 2019; 9, 407. [DOI: https://dx.doi.org/10.3390/min9070407]

27. Evrard, M.; Dumont, G.; Hermans, T.; Chouteau, M.; Francis, O.; Pirard, E.; Nguyen, F. Geophysical Investigation of the Pb–Zn Deposit of Lontzen–Poppelsberg, Belgium. Minerals; 2018; 8, 233. [DOI: https://dx.doi.org/10.3390/min8060233]

28. Pilkington, M.; Shamsipour, P. Noise reduction procedures for gravity-gradiometer data. Geophysics; 2014; 79, pp. G69-G78. [DOI: https://dx.doi.org/10.1190/geo2014-0084.1]

29. Liao, J.; Liu, H.; Li, W.; Guo, Z.; Wang, L.; Wang, H.; Peng, S.; Hursthouse, A. 3-D Butterworth Filtering for 3-D High-density Onshore Seismic Field Data. J. Environ. Eng. Geophys.; 2018; 23, pp. 223-233. [DOI: https://dx.doi.org/10.2113/JEEG23.2.223]

30. Davis, K.; Li, Y. Efficient 3D inversion of magnetic data via octree-mesh discretization, space-filling curves, and wavelets. Geophysics; 2013; 78, pp. J61-J73. [DOI: https://dx.doi.org/10.1190/geo2012-0192.1]

31. Cox, L.H.; Wilson, G.A.; Zhdanov, M.S. 3D inversion of airborne electromagnetic data using a moving footprint. Explor. Geophys.; 2018; 41, pp. 250-259. [DOI: https://dx.doi.org/10.1071/EG10003]

32. Shin, Y.H.; Choi, K.S.; Xu, H. Three-dimensional forward and inverse models for gravity fields based on the Fast Fourier Transform. Comput. Geosci.; 2006; 32, pp. 727-738. [DOI: https://dx.doi.org/10.1016/j.cageo.2005.10.002]

33. Li, Y.; Oldenburg, D. Fast inversion of large scale magnetic data using wavelets. SEG Tech. Program Expand. Abstr.; 1997; pp. 490-493. [DOI: https://dx.doi.org/10.1190/1.1885941]

34. Rene, R.M. Gravity inversion using open, reject, and “shape-of-anomaly” fill criteria. Geophysics; 1986; 51, pp. 988-994. [DOI: https://dx.doi.org/10.1190/1.1442157]

35. Camacho, A.G.; Montesinos, F.G.; Vieira, R. Gravity inversion by means of growing bodies. Geophysics; 2000; 65, pp. 95-101. [DOI: https://dx.doi.org/10.1190/1.1444729]

36. Uieda, L.; Barbosa, V. Robust 3D gravity gradient inversion by planting anomalous densities. Geophysics; 2012; 77, pp. G55-G66. [DOI: https://dx.doi.org/10.1190/geo2011-0388.1]

37. Uieda, L.; Barbosa, V.C.F. Use of the “shape-of-anomaly” data misfit in 3D inversion by planting anomalous densities. Proceedings of the SEG Technical Program Expanded Abstracts 2012; Las Vegas, NV, USA, 4–9 November 2012; [DOI: https://dx.doi.org/10.1190/segam2012-0383.1]

38. Goossens, S.; Lemoine, F.G.; Sabaka, T.J.; Nicholas, J.B.; Mazarico, E.; Rowlands, D.D.; Loomis, B.D.; Chinn, D.S.; Neumann, G.A.; Smith, D.E. A Global Degree and Order 1200 Model of the Lunar Gravity Field Using GRAIL Mission Data; In Proceedings of the 47th Lunar and Planetary Science Conference, The Woodlands, TX, USA, 21-25, March, 2016; p. 1484.

39. Gómez-Ortiz, D.; Agarwal, B.N.P. 3DINVER.M: A MATLAB program to invert the gravity anomaly over a 3D horizontal density interface by Parker–Oldenburg’s algorithm. Comput. Geosci.; 2005; 31, pp. 513-520. [DOI: https://dx.doi.org/10.1016/j.cageo.2004.11.004]

40. Wieczorek, M.A.; Neumann, G.A.; Nimmo, F.; Kiefer, W.S.; Taylor, G.J.; Melosh, H.J.; Phillips, R.J.; Solomon, S.C.; Andrews-Hanna, J.C.; Asmar, S.W. et al. The Crust of the Moon as Seen by GRAIL. Science; 2013; 339, pp. 671-675. [DOI: https://dx.doi.org/10.1126/science.1231530]

41. Last, B.; Kubik, K. Compact gravity inversion. Geophysics; 1983; 48, pp. 713-721. [DOI: https://dx.doi.org/10.1190/1.1441501]

42. Barbosa, V.; Silva, J. Generalized compact gravity inversion. Geophysics; 1994; 59, pp. 57-68. [DOI: https://dx.doi.org/10.1190/1.1443534]

43. Ghalehnoee, M.H.; Ansari, A.; Ghorbani, A. Improving compact gravity inversion based on new weighting functions. Geophys. J. Int.; 2016; 208, pp. 546-560. [DOI: https://dx.doi.org/10.1093/gji/ggw413]

44. Stocco, S.; Godio, A.; Sambuelli, L. Modelling and compact inversion of magnetic data: A Matlab code. Comput. Geosci.; 2009; 35, pp. 2111-2118. [DOI: https://dx.doi.org/10.1016/j.cageo.2009.04.002]

45. LaBrecque, D.J.; Owen, E.; Dailey, W.; Ramirez, A.L. Noise and occam’s inversion of resistivity tomography data. SEG Technical Program Expanded Abstracts 1992; Society of Exploration Geophysicists: New Orleans, LA, USA, 1992; pp. 397-400. [DOI: https://dx.doi.org/10.1190/1.1822100]

46. Siripunvaraporn, W.; Sarakorn, W. An efficient data space conjugate gradient Occam’s method for three-dimensional magnetotelluric inversion. Geophys. J. Int.; 2011; 186, pp. 567-579. [DOI: https://dx.doi.org/10.1111/j.1365-246X.2011.05079.x]

47. Muller, P.M.; Sjogren, W.L. Mascons: Lunar mass concentrations. Science; 1968; 161, pp. 680-684. [DOI: https://dx.doi.org/10.1126/science.161.3842.680]

48. Klokočník, J.; Kostelecký, J.; Cílek, V.; Kletetschka, G.; Bezděk, A. Gravity aspects from recent gravity field model GRGM1200A of the Moon and analysis of magnetic data. Icarus; 2022; 384, 115086. [DOI: https://dx.doi.org/10.1016/j.icarus.2022.115086]

49. Yuan, F.; Xiao-hong, M.; Jun, W.; Zhao-xi, C. Study of the characteristics of lunar gravity field based on the model GRGM660PRIM. Prog. Geophys.; 2018; 33, pp. 2211-2218. [DOI: https://dx.doi.org/10.6038/pg2018CC0118]

50. Spudis, P.D.; Sliz, M.U. Impact melt of the lunar Crisium multiring basin. Geophys. Res. Lett.; 2017; 44, pp. 1260-1265. [DOI: https://dx.doi.org/10.1002/2016GL071429]

51. Rodionova, Z.F.; Zharkova, A.Y.; Grishakina, E.A.; Shevchenko, V.V. Topographic Features of the Lunar Maria and Basins. Sol. Syst. Res.; 2021; 55, pp. 183-199. [DOI: https://dx.doi.org/10.1134/S0038094621030059]

52. Wilhelms, D.E.; McCauley, J.F.; Trask, N.J. The Geologic History of the Moon; Professional Paper 1348; United States Government Printing Office: Washington, DC, USA, 1987.

53. Lu, X.; Cao, H.; Ling, Z.; Fu, X.; Qiao, L.; Chen, J. Geomorphology, Mineralogy, and Geochronology of Mare Basalts and Non-Mare Materials around the Lunar Crisium Basin. Remote Sens.; 2021; 13, 4828. [DOI: https://dx.doi.org/10.3390/rs13234828]

54. Smith, D.E.; Zuber, M.T.; Neumann, G.A.; Lemoine, F.G.; Mazarico, E.; Torrence, M.H.; McGarry, J.F.; Rowlands, D.D.; Head, J.W.; Duxbury, T.H. et al. Initial observations from the Lunar Orbiter Laser Altimeter (LOLA). Geophys. Res. Lett.; 2010; 37, pp. 1-6. [DOI: https://dx.doi.org/10.1029/2010GL043751]

55. Lemoine, F.G.; Goossens, S.; Sabaka, T.J.; Nicholas, J.B.; Mazarico, E.; Rowlands, D.D.; Loomis, B.D.; Chinn, D.S.; Neumann, G.A.; Smith, D.E. et al. GRGM900C: A degree 900 lunar gravity model from GRAIL primary and extended mission data. Geophys. Res. Lett.; 2014; 41, pp. 3382-3389. [DOI: https://dx.doi.org/10.1002/2014GL060027]

56. Deng, Q.; Li, F.; Yan, J.; Ye, M.; Barriot, J.P. Degree-Depth Relation for Planetary Gravity Field Model Based on Wavelength. Earth Space Sci.; 2022; 9, e2021EA002143. [DOI: https://dx.doi.org/10.1029/2021EA002143]

57. Featherstone, W.E.; Hirt, C.; Kuhn, M. Band-limited Bouguer gravity identifies new basins on the Moon. J. Geophys. Res. Planets; 2013; 118, pp. 1397-1413. [DOI: https://dx.doi.org/10.1002/jgre.20101]

58. Goossens, S.; Sabaka, T.J.; Wieczorek, M.A.; Neumann, G.A.; Mazarico, E.; Lemoine, F.G.; Nicholas, J.B.; Smith, D.E.; Zuber, M.T. High-Resolution Gravity Field Models from GRAIL Data and Implications for Models of the Density Structure of the Moon’s Crust. J. Geophys. Res. Planets; 2020; 125, e2019JE006086. [DOI: https://dx.doi.org/10.1029/2019JE006086]

59. Jansen, J.C.; Andrews-Hanna, J.C.; Li, Y.; Lucey, P.G.; Taylor, G.J.; Goossens, S.; Lemoine, F.G.; Mazarico, E.; Head, J.W.; Milbury, C. et al. Small-scale density variations in the lunar crust revealed by GRAIL. Icarus; 2017; 291, pp. 107-123. [DOI: https://dx.doi.org/10.1016/j.icarus.2017.03.017]

60. Jozwiak, L.M.; Head Iii, J.W.; Neumann, G.A.; Wilson, L. Observational constraints on the identification of shallow lunar magmatism: Insights from floor-fractured craters. Icarus; 2017; 283, pp. 224-231. [DOI: https://dx.doi.org/10.1016/j.icarus.2016.04.020]

61. Andrews-Hanna, J.C.; Head, J.W.; Johnson, B.C.; Keane, J.T.; Kiefer, W.S.; McGovern, P.J.; Neumann, G.A.; Wieczorek, M.A.; Zuber, M.T. Ring faults and ring dikes around the Orientale basin on the Moon. Icarus; 2018; 310, pp. 1-20. [DOI: https://dx.doi.org/10.1016/j.icarus.2017.12.012]

62. Hiesinger, H. New Views of Lunar Geoscience: An Introduction and Overview. Rev. Mineral. Geochem.; 2006; 60, pp. 1-81. [DOI: https://dx.doi.org/10.2138/rmg.2006.60.1]

63. Stoffler, D. Cratering History and Lunar Chronology. Rev. Mineral. Geochem.; 2006; 60, pp. 519-596. [DOI: https://dx.doi.org/10.2138/rmg.2006.60.05]

64. Golombek, M.P.; Anderson, F.S.; Zuber, M.T. Martian wrinkle ridge topography: Evidence for subsurface faults from MOLA. J. Geophys. Res. Planets; 2001; 106, pp. 23811-23821. [DOI: https://dx.doi.org/10.1029/2000JE001308]

65. Taylor, S.R.; McLennan, S. Planetary Crusts: Their Composition, Origin and Evolution; Cambridge University Press: NY, UK, 2009; Volume 10.

66. Spector, A.; Grant, F.S. Statistical models for interpreting aeromagnetic data. Geophysics; 1970; 35, pp. 293-302. [DOI: https://dx.doi.org/10.1190/1.1440092]

67. Mousa, A.; Mickus, K.; Al-Rahim, A. The thickness of cover sequences in the Western Desert of Iraq from a power spectrum analysis of gravity and magnetic data. J. Asian Earth Sci.; 2017; 138, pp. 230-245. [DOI: https://dx.doi.org/10.1016/j.jseaes.2017.02.022]

68. Ganguli, S.S.; Pal, S.K.; Singh, S.L.; Rao, J.; Balakrishna, B. Insights into crustal architecture and tectonics across Palghat Cauvery Shear Zone, India from combined analysis of gravity and magnetic data. Geol. J.; 2020; 56, pp. 2041-2059. [DOI: https://dx.doi.org/10.1002/gj.4041]

69. Ganguli, S.S.; Pal, S.K.; Sundaralingam, K.; Kumar, P. Insights into the crustal architecture from the analysis of gravity and magnetic data across Salem-Attur Shear Zone (SASZ), South-ern Granulite Terrane (SGT), India: An evidence of accretional tectonics. Episodes; 2021; 44, pp. 419-442. [DOI: https://dx.doi.org/10.18814/epiiugs/2020/020095]

70. de Souza, J.; Pires Rostirolla, S. A fast MATLAB program to estimate the multifractal spectrum of multidimensional data: Application to fractures. Comput. Geosci.; 2011; 37, pp. 241-249. [DOI: https://dx.doi.org/10.1016/j.cageo.2010.09.001]

71. Essa, K.S.; Munschy, M.; Youssef, M.A.S.; Khalaf, E.E.D.A. Aeromagnetic and Radiometric Data Interpretation to Delineate the Structural Elements and Probable Precambrian Mineralization Zones: A Case Study, Egypt. Min. Metall. Explor.; 2022; 39, pp. 2461-2475. [DOI: https://dx.doi.org/10.1007/s42461-022-00675-0]

72. Wieczorek, M.A.; Jolliff, B.L.; Khan, A.; Pritchard, M.E.; Weiss, B.P.; Williams, J.D.; Hood, L.L.; Righter, K.; Neal, C.R.; Shearer, C.K. New views of the moon. Rev. Min. Geochem; 2006; 60, pp. 221-364. [DOI: https://dx.doi.org/10.2138/rmg.2006.60.3]

73. Wieczorek, M.A. Gravity and topography of the terrestrial planets. Treatise Geophys.; 2015; 10, pp. 153-193. [DOI: https://dx.doi.org/10.1016/B978-0-444-53802-4.00169-X]

74. Besserer, J.; Nimmo, F.; Wieczorek, M.A.; Weber, R.C.; Kiefer, W.S.; McGovern, P.J.; Andrews-Hanna, J.C.; Smith, D.E.; Zuber, M.T. GRAIL gravity constraints on the vertical and lateral density structure of the lunar crust. Geophys. Res. Lett.; 2014; 41, pp. 5771-5777. [DOI: https://dx.doi.org/10.1002/2014GL060240]

75. Wieczorek, M.A.; Beuthe, M.; Rivoldini, A.; Van Hoolst, T. Hydrostatic Interfaces in Bodies With Nonhydrostatic Lithospheres. J. Geophys. Res. Planets; 2019; 124, pp. 1410-1432. [DOI: https://dx.doi.org/10.1029/2018JE005909]

76. Kiefer, W.S.; Macke, R.J.; Britt, D.T.; Irving, A.J.; Consolmagno, G.J. The density and porosity of lunar rocks. Geophys. Res. Lett.; 2012; 39, L07201. [DOI: https://dx.doi.org/10.1029/2012GL051319]

77. Corley, L.M.; McGovern, P.J.; Kramer, G.Y.; Lemelin, M.; Trang, D.; Gillis-Davis, J.J.; Taylor, G.J.; Powell, K.E.; Kiefer, W.S.; Wieczorek, M. Olivine-bearing lithologies on the Moon: Constraints on origins and transport mechanisms from M3 spectroscopy, radiative transfer modeling, and GRAIL crustal thickness. Icarus; 2018; 300, pp. 287-304. [DOI: https://dx.doi.org/10.1016/j.icarus.2017.09.012]

78. May, D.A.; Knepley, M.G. Optimal, scalable forward models for computing gravity anomalies. Geophys. J. Int.; 2011; 187, pp. 161-177. [DOI: https://dx.doi.org/10.1111/j.1365-246X.2011.05167.x]

79. Tian, Y.; Wang, Y. Inversion of the density structure of the lithosphere in the North China Craton from GOCE satellite gravity gradient data. Earth Planets Space; 2018; 70, pp. 1-22. [DOI: https://dx.doi.org/10.1186/s40623-018-0942-1]

80. Bühlmann, P.; Yu, B. Boosting with the L 2 loss: Regression and classification. J. Am. Stat. Assoc.; 2003; 98, pp. 324-339. [DOI: https://dx.doi.org/10.1198/016214503000125]