Fast Magnetization Vector Inversion Method with Undulating Observation Surface in Spherical Coordinate for Revealing Lunar Weak Magnetic Anomaly Feature

Abstract

The three-dimensional magnetic vector structure (magnetization intensity and direction) of the planet can be effectively used to analyze the characteristics of its formation and operation. However, the quick acquisition of a large region of the magnetic vector structure of the planet with bigger observation surfaces undulation is hard and indispensable. We firstly proposed a fast magnetization vector inversion method for the inversion of a magnetic anomaly with the undulating observation surfaces in the spherical coordinate system, which first transforms the data to a plane when the data are distributed on a surface. Then, it uses a block-Toeplitz-Toeplitz-block (BTTB)-FFT to achieve fast inversion with the constraint that the magnetization intensities of the grids between the transformed observation surfaces and the terrain are zero. In addition, Gramian constraint term is used to reduce the ambiguity of the magnetic vector inversion. The theoretical model tests show that the proposed method can effectively improve the computational efficiency by 23 times in the 60 × 60 × 10 grid division compared to the conventional inversion method, and the accuracy of the two computation methods is comparable. The root-mean-square error of the magnetization intensity is only 0.017, and the angle error is within 1°. The magnetization vector structure shows that the largest crater diameter does not exceed 340 km in the Mare Australe region, the amplitude of the magnetic anomaly is much higher than the current meteorite impact simulation results, and the depth of the magnetic source is less than 10 km, which cannot be explained by the impact simulation experiments. In addition, the magnetization directions of adjacent sources differ by 122° (or 238°), and the high-frequency dynamics of the Moon as well as the short-lived dynamics may be responsible for this phenomenon. The magnetization directions of the three adjacent sources in the Mare Crisium region are close to each other and differ in depth with different cooling times, making it difficult to record the transient fields produced by meteorite impacts. In addition to the above characteristics, the magnetization direction of the magnetic sources in both regions is uniformly distributed without reflecting the dispersion of the magnetization direction of the meteorite impact magnetic field. Therefore, it can be inferred that the magnetic anomalies in these two regions are related to the generator hypothesis.

1. Introduction

Compared with the conventional magnetic anomalies inversion method, which can only obtain the magnitude of magnetization, the three-dimensional magnetic vector inversion can obtain the three-dimensional magnetic vector structure, which contains the magnitude and the direction of magnetization [1,2,3,4,5]. The three-dimensional magnetic vector structure can obtain more information of the subsurface magnetic structure, which is more helpful to explain the formation and evolution of the subsurface structure. The three-dimensional magnetic vector inversion method was first proposed by Lelièvre and Oldenburg [1], but it requires inversion of three physical quantities, which triples the ambiguity compared with the conventional magnetic anomaly inversion method. Jorgensen and Zhdanov [2] proposed a magnetization vector inversion method based on Gramian constraints to reduce the non-uniqueness of the inversion due to the increase in the inversion parameters. However, the current magnetic vector inversion method only considers the case of Cartesian coordinate system discretization, but the influence of planet curvature cannot be neglected when inverting satellite data in large regions [6,7,8,9]. Thus, a three-dimensional magnetic vector inversion system in a spherical coordinate system is needed. When the inversion data amount is large, in order to improve the practicality of the inversion algorithm, computer hardware (GPU and CPU) can be utilized for parallel computation to increase the inversion speed [10,11,12,13]. Due to the large amount of large region satellite data and the slow inversion speed, Zhao et al. [14,15] proposed a highly efficient gravity inversion method based on the block-Toeplitz Toeplitz-block (BTTB) structure and the fast Fourier transform (FFT) method for spherical coordinate systems that reduces the inversion time and improves the practicality of satellite data inversion. In this paper, we mainly use the BTTB-FFT algorithm to improve the efficiency; of course, parallel computation can also be applied on this basis to further improve the efficiency. In addition, when the large regional data inversion is carried out, the large observation surface undulation may affect the inversion, so the inversion strategy that considers the observation surface undulation is more consistent with the real situation [16,17,18].

The origin of the lunar weak magnetic anomaly has always been a scientific concern. In the 1960s and 1970s, the measurements by the Apollo missions revealed that the Moon does not have a global basic magnetic field [19], which made scientists have a strong interest in the origin of the weak lunar magnetic anomaly [20,21,22]. Daily and Dyal considered two hypotheses, impact-related shock and liquid nuclear dynamics, as the most likely causes [23].

At first, the lunar magnetic field was considered to have characteristics similar to the geomagnetic field. Runcorn et al. [20] and Strangway et al. [21] proposed that magnetized fields are global in nature and originally formed when a dynamo was active in the core of the Moon, and the core dynamo was powered by heat convection in the Moon’s liquid iron core. Moreover, most paleomagnetic studies support the idea that the Moon once had a global magnetic field, and impact magnetization is secondary [24,25,26,27,28,29]. However, some authors had a different view and suggested that lunar materials have acquired magnetization only via transient magnetic fields generated by hypervelocity impacts [22,30,31,32,33,34,35,36]. Hood et al. revised some of the earlier proposed dominance of magnetization from meteorite impacts and concluded that numerous Nectarian basins contain thermal remanent magnetization from a dynamo field rather than magnetization from transient impact fields [37].

In recent years, Evans et al. [38] performed numerical simulations of possible lunar dynamos with thermochemical convection, and the simulations showed that lunar crustal magnetization was obtained via the ancient dynamo field. Crawford [39] concluded from numerical simulations that magnetic anomalies in the lunar Crisium, Nectaris, Serenitatis, Humboldtianum, and Mendel-Rydberg basins can be explained by magnetization generated by charge separation associated with a 20 km radius impactor. Oran et al. [40] performed magnetohydrodynamic and impact simulations, and the results showed that the magnetic fields generated by impacts are too weak to be effectively explained by the lunar crustal magnetic anomalies, so the core dynamo is considered to be the only reasonable origin of most magnetization on the Moon. Tarduno et al. [41] performed impact simulations and found that the field produced by the impact itself was consistent with the ancient intensity value of the 64,455 Apollo sample, suggesting that the Apollo lunar samples may have recorded an external field produced by the impact rather than an ancient core dynamo. The debate over the above two hypotheses is still ongoing.

Satellite magnetic data are an effective way to reveal the internal structure of the planet. The three-dimensional magnetic vector structure can explain the causes of lunar magnetic anomalies by analyzing the depth of the magnetic source and the information of the direction of the magnetization to reveal the mechanism of the magnetic field evolution [42].

We firstly proposed a fast magnetization vector inversion method for the data with undulating observation surfaces in a spherical coordinate system. The influence of the undulating observation surfaces is considered, the BTTB-FFT algorithm is used to improve the efficiency of the inversion to enhance the practicality, and the Gramian constraint is introduced to reduce the ambiguity of the inversion. The three-dimensional magnetic vector structures of the Mare Australe and the Mare Crisium regions of the Moon are obtained by this method, and the results show that the dynamo hypothesis may be a more reasonable explanation for the origin of the lunar magnetic anomaly.

2. A Fast Magnetization Vector Inversion Method with an Undulating Observation Surface in a Spherical Coordinate System

In real situations, due to the fact that spherical coordinate system inversion is generally used under large regional conditions, there may be significant terrain undulation in the research region, which can cause undulation of the observation surface when measuring on the surface or along the terrain. As shown in Figure 1a, when the observation surface is undulating, the conventional method is to directly perform a regular grid division from the undulating observation surface and then obtain its kernel function matrix by accumulating observation points one by one before performing inversion calculations. When conducting a precise inversion of the research region, it is necessary to refine the grid, and the generation of kernel functions takes a long time. The above method does not facilitate the application of fast calculation methods because the grid is not uniformly distributed in space, resulting in a correspondence of the observation points and the grid that is not one-to-one.

We developed a fast magnetization vector inversion method for the inversion of a magnetic anomaly with the undulating observation surface in the spherical coordinate system. Firstly, the data are transformed from the undulating surface to the plane by the equivalent source method [43], as shown in Figure 1b. Then, block-Toeplitz-Toeplitz-block (BTTB)-FFT is used with the constraint that the magnetization intensities of the grids between the real observation surface and the terrain are zero to achieve fast inversion. The objective function for inversion in this method is the following:

$$\begin{array}{ll}\phi & ={\phi }_d+\alpha {\phi }_m+\beta {\phi }_G\\ & ={‖KM-B‖}_2^2+\alpha {‖W(M-M_{\mathrm{ref}})‖}_2^2+\beta {\phi }_G\\ & ={‖\left[\begin{array}{ccc}{K}_{{\theta }_{M_{\theta }}}& K_{{\theta }_{M_{\varphi }}}& K_{{\theta }_{M_r}}\\ K_{{\varphi }_{M_{\theta }}}& K_{{\varphi }_{M_{\varphi }}}& K_{{\varphi }_{M_r}}\\ K_{r_{M_{\theta }}}& K_{r_{M_{\varphi }}}& K_{r_{M_r}}\end{array}\right]\left[\begin{array}{c}{M}_{\theta }\\ M_{\varphi }\\ M_r\end{array}\right]-\left[\begin{array}{c}{B}_{\theta }\\ B_{\varphi }\\ B_r\end{array}\right]‖}_2^2+\\ & \alpha {‖diag(W_v{W}_d)(\left[\begin{array}{c}{M}_{\theta }\\ M_{\varphi }\\ M_r\end{array}\right]-M_{\mathrm{ref}})‖}_2^2+\beta \left[\begin{array}{c}{\phi }_{G_{\theta }}\\ {\phi }_{G_{\varphi }}\\ {\phi }_{G_r}\end{array}\right].\end{array}$$

(1)

${\phi }_d$ is the data objective function, and ${\phi }_m$ is the model objective function. α represents the regularization parameters. The influence of planetary curvature cannot be ignored when inverting satellite data over large regions, so the tesseroid is used as the grid cell in the spherical coordinate system inversion, as shown in Figure 1c [44,45]. Since the Moon has a very small oblateness of 1/963.7256, and its equatorial and polar radii are very close to each other, so it is closer to a sphere. Considering the planetary curvature, the kernel matrix for magnetic vector inversion in the spherical coordinate system is as follows [46]:

$$\begin{array}{l}{K}_{i_{M_j}}=\frac{{\mu }_0}{4\pi }{\displaystyle \underset{{\theta }_1}{\overset{{\theta }_2}{\int }}{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{M}_j(\frac{3{\Delta }_i{\Delta }_j}{l^5}-\frac{{\delta }_{ij}}{l^3})\tau d{r}^{\prime }d{\varphi }^{\prime }d}}}{\theta }^{\prime }\\ {\Delta }_{\theta }=r^{\prime }[\cos \varphi \sin {\varphi }^{\prime }-\sin \varphi \cos {\varphi }^{\prime }\cos ({\theta }^{\prime }-\theta )]\\ {\Delta }_{\varphi }=r^{\prime }\cos {\varphi }^{\prime }\sin ({\theta }^{\prime }-\theta )\\ {\Delta }_{\lambda }=r^{\prime }\cos \psi -r\\ \tau ={r^{\prime }}^2\cos {\varphi }^{\prime }\\ l=\sqrt{{r^{\prime }}^2+r^2-2r^{\prime }r\cos \psi }\\ \cos \psi =\sin \varphi \sin {\varphi }^{\prime }+\cos \varphi \cos {\varphi }^{\prime }\cos ({\theta }^{\prime }-\theta )\end{array}$$

(2)

$K_{i_{M_j}}$ (i, j = θ, ϕ, r, and θ collectively represents longitudinal direction, ϕ represents latitudinal direction, r represents radial direction) is the kernel matrix and represents the response of the tesseroid center (θ, ϕ, r) to the observation (θ′, ϕ′, r′), and ${\delta }_{ij}$ is the Kronecker function (when i = j, ${\delta }_{ij}$ = 1; when i ≠ j, ${\delta }_{ij}$ = 0). In the inversion process, the large condition number of $K_{i_{M_j}}$ makes it an ill-posed problem, and the slightest perturbation of the data can make the inversion results very different. Tikhonov and Arsenin proposed the regularized inversion technique, which improves the ill-posed problem of the inversion and improves the stability of the inversion. At the same time, it introduces constraints through the regularization term and then improves the multisolvability problem of the inversion, which to some extent solves the instability and the multisolvability of the inversion problem at the same time [47]. M_j (j = θ, ϕ, r) and B_j (j = θ, ϕ, r) are the magnetization intensity and the magnetic anomaly in different directions, respectively. M_ref is the reference model.

Regarding the determination of the regularization parameter α, currently, there are mainly generalized cross validation (GCV) and L-curve methods. The optimal parameter is obtained when the function obtains the minimum value in the GCV method, and the L-curve method is generally considered to be optimal at the inflection point of the curve (at the maximum curvature). Although GCV is a good solution to many problems, in some cases, it is difficult to find a good regularization parameter for this method. The variation in the GCV function is sometimes too flat when it is difficult to locate its minimum value. Another problem is that GCV sometimes mistakes correlated noise for a signal, and if the errors are highly correlated, the method is likely to give unsatisfactory results [48]. In contrast, the calculation of the inflection points in the L-curve method is a well-defined numerical problem, and the method is rarely “fooled” by correlation errors and has good stability. Therefore, this method is adopted as the regularization parameter selection [49,50].

As can be seen in Figure 1b, the volume of the tesseroid is the same in the latitudinal direction, and the volume changes only with longitude and depth. To effectively overcome the inconsistency of volume and the difference of the depth, which leads to low resolution, we use the depth and the volume as the weight functions for inversion, and the formulas are as follows:

$$W_v={(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$V$}\right.)}^m,W_d={(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$D$}\right.)}^n.$$

(3)

W_v is the volume weight, V and m are the volume of tesseroids and the factor, respectively, and m is usually chosen to be one. In addition, we introduce a depth weighting matrix, W_d, in order to balance the effect of depth variation in inversion, where D and n are the depth of tesseroids and the factor, respectively, and n is usually taken to be 0.5.

${\phi }_{G_i}(i=\theta ,\varphi ,r)$ collectively represent the Gramian constraints ${\phi }_{G_i}(i=\theta ,\varphi ,r)$, which are used to improve the accuracy of the inversion results. β is the weighting coefficient for the Gramian constraint term, which is usually chosen as 1000 based on experience.

Next, a regular grid division is performed based on the transformed observation surface. Under this grid division, as shown in Formula (4), we introduce the block-Toeplitz-Toeplitz-block structure and the fast Fourier transform (BTTB-FFT) fast algorithm [14,15], which converts matrix multiplication operations into convolution operations to achieve acceleration effects by reducing the computational complexity. In Figure 1d, we can see that the kernel function matrix satisfies the properties of the Toeplitz matrix, and the same color elements have equivalence. However, this method requires a one-to-one correspondence between the grids and the observation points [14,15], so the conventional method shown in Figure 1a is not suitable for this fast algorithm. In addition, when performing the fast algorithm, we use the positional relationship between the terrain and the center point of the grid as a constraint strategy. Grids with center points higher than the terrain are used as air layers, and grids lower than the terrain are involved in the computation. By setting the physical property parameter, M_j, of the air layer between the observation surface and the undulating terrain to zero, we can eliminate the invalid grids, which is more in line with the real situation.

$$B_i^{ext}={\displaystyle \sum _j{K}_{i_{M_j}}^{ext}{M}_i^{ext}}=\sum f^{-1}[f(k_i^{ext})f(M_i^{ext})],(i,j=\theta ,\varphi ,r)$$

(4)

$$K_{i_{M_j}}^{ext}=\left[\begin{array}{ccccc}{k}_0& k_{n-1}& \cdots & k_2& k_1\\ k_1& k_0& k_{n-1}& \ddots & k_2\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ k_{n-2}& \ddots & k_1& k_0& k_{n-1}\\ k_{n-1}& k_{n-2}& \cdots & k_1& c_0\end{array}\right]$$

(5)

where $B_i^{ext}={\left[\begin{array}{cc}{B}_i& 0_{n\times 1}\end{array}\right]}^T$ and $M_i^{ext}={\left[\begin{array}{cc}{M}_i& 0_{n\times 1}\end{array}\right]}^T$ are the extension vectors of B_i and M_i, respectively. $K_{i_{M_j}}^{ext}$ is the cyclic matrix, which is a special type of Toeplitz matrix, where the elements can be written in the form of a vector $k^{ext}=(k_0,k_1,\dots ,k_{n-1},0,k_{1-n},\dots k_{-1})$.

3. Theoretical Model Tests

To verify the accuracy of the inversion method, a single model is firstly set up. The model positions are at 39°N–41°N, 39°E–41°E, and depths of 20–50 km. The magnetization intensity of the model is 0.5 A/m, and the magnetization declination and inclination are 45° and −45°, respectively. Satellite data measurements are made in the spherical coordinate system, which expresses coordinate positions in terms of radial distance (the distance from a spherical coordinate point to the center of the sphere, denoted by r), polar angle (the angle between the z-axis and r, generally denoted by ϕ), and the azimuth angle (the angle on the equatorial plane (the plane defined by the x-axis and the y-axis) that begins at the x-axis and is measured in an anticlockwise direction, usually denoted by θ). Therefore, data in the radial direction, in the direction of the polar angle, and in the direction of the azimuth angle are expressed in terms of r, θ, and φ, respectively, in the measurement data. The B_θ, B_ϕ, and B_r components and the magnetic anomaly of the model are shown in Figure 2a–d.

From the total inversion objective function, it can be seen that the regularization factor, α, serves to weigh the relative degree of data fitting and model fitting. When α is small, the data fitting part has a larger weight, and the model fitting part has a smaller weight, so the inversion model fits the observed data more, and the inversion result tends to the surface; when α is larger, the inversion process fits the model objective function more, and the inversion result tends to the depth. Therefore, the selection of the regularization parameter α is crucial for the inversion results. Compared with the GVC method, the calculation of the inflection point of the L-curve with respect to the selection of the regularization parameter α is a well-defined numerical problem, which is seldom fooled by the correlation error and has better stability [48,49].

The L-curve obtained in the inversion process is shown in Figure 3d. Its maximum curvature corresponds to the best regularization factor α = 10, and the corresponding inversion result is also considered to be the best inversion result.

In this model, the condition number of the kernel matrix is as high as 3.9 × 10¹³, which is a seriously ill-posed problem, and it is difficult to get a reliable solution using the general solution method. Furthermore, it is easy to cause the relative error of the solution to be very large and deviate from the real solution, so it is necessary to add the regularization method to solve the problem. After adding regularization, the matrix condition number is reduced to 65.5, and a relatively reliable solution can be obtained by the inverse algorithm [47].

We performed the fast inversion after the transformation of the observation surface and the conventional inversion, respectively, and the results are shown in Figure 3a,b. It is difficult to distinguish the difference between the two results with the naked eye. The residual plot between the two is shown in Figure 3c, with a maximum absolute value of 0.02 A/m and a root-mean-square error of 0.017. Overall, the difference is relatively small. The declinations of the two results are 43.11° and 44.03°, and the inclinations are −47.34° and −46.47°, respectively. The angles of the two results are also close to each other, within 1°, and both are close to the real model, with an error from the real inclination and declination of less than 5°. From the above two aspects of magnetization intensity and magnetization direction, the inversion accuracies of the two methods are equivalent.

However, as shown in Figure 3e, compared to the conventional method, the efficiency of the fast magnetization vector inversion method is higher, the improvement of the inversion efficiency is increased with the increasing number of grid subdivisions, and the efficiency is improved by 23 times with the grids of 60 × 60 × 10. Therefore, the method in this paper can improve the efficiency without decreasing the accuracy, and it has better practicality.

In the real data inversion process, we cannot guarantee that the data do not contain noise, and the existence of noise has a great impact on the inversion results. Therefore, we need to test the anti-noise aspect of the method. To simulate the noise, we add 5% Gaussian noise to the magnetic anomaly component data in Figure 2a–c. Figure 4a–d displays the noise-corrupted component of the magnetic anomalies and the magnetic anomalies B, and the noise has an obvious impact on the anomaly shape.

From the inversion slice at 40°N in Figure 5e, it can be seen that the convergence of both the magnetization intensity results and the magnetization direction results is reduced compared to Figure 3a due to the influence of noise, and the assignment of the magnetization intensity is only recovered to about 30%. However, the overall shape and the position of the magnetization intensity results are still acceptable without much deviation, and the average magnetization declination and inclination within the model range are 36.29° and −51.16°, respectively, which are within 10° of the true inclination and declination. Furthermore, the error is within the acceptable range compared to the variation range of inclination of 180° and declination of 360°.

In addition to this, we also compared the results for different noise conditions. Shown in Figure 5 is a box plot of the results of inclination and declination under different noise conditions. Box plots are commonly used statistical charts that show the dispersion of data, which can be used to visually compare the distribution characteristics of multiple sets of data. Figure 5 shows a complete box plot style, which has five characteristic positions as upper and lower limits, upper and lower quartiles, and median. The median is the number in the middle of the dataset in ascending order; when the data set is even, the middle is two numbers, and then the average of the two numbers is taken as the median. The median position is indicated by the red line in Figure 5. The positions of the upper and the lower quartiles are the positions of 3(n + 1)/4 and (n + 1)/4 when the data are sorted in ascending order, where n is the total number of numbers in the dataset. These are shown at the top and bottom of the blue box in Figure 5. The upper and lower quartiles contain about half of the data in the entire dataset, reflecting the main distributional characteristics of the data. The difference between the upper and lower quartiles is the interquartile range (IQR). The upper and lower limits are the positions of the upper and lower quartiles after expanding the IQR distance 1.5 times to the ends of the data arranged in ascending order, which are marked with black lines in Figure 5. It can be seen that the smaller the IQR is in the box plot, the more concentrated the data are at the median position. Due to the fact that inversion is the process of solving large underdetermined equations, we see that the inversion results are in the neighborhood of the true values, as shown by the green line in Figure 5. However, as the noise increases, the ranges of upper and lower limits and IQR become larger, and the values become more scattered with instability. From Figure 5a, it can be seen that the dispersion of the declination changes significantly when the noise reaches more than 10%. In addition, it can be seen from Figure 5b that the dispersion of the inclination changes significantly when the noise reaches 15% or more. It indicates that, when the noise reaches a certain level, the inversion results are more sensitive to changes in the noise and are more affected by the noise.

Considering the complexity of the real situation, we designed a model that contains three anomalous bodies to test the applicability of the method. The model positions are at 39°N–41°N, 33°E–36°E, 39°N–41°N, 39°E–42°E, and 44°E–47°E. Models one and three are at depths of 30–60 km with magnetization intensities of 0.75 A/m, model two is at a depth of 20–50 km with a magnetization intensity of 1 A/m, and the inclinations and declinations of all models are 45° and −45°. The component of the magnetic anomalies and the magnetic anomalies B of the model are shown in Figure 6a–d.

The magnetic vector structure computed by the fast magnetization vector method is shown in Figure 6e. The declination results are 44.37°, 44.36°, and 44.31°, and the inclination results are −47.36°, −47.37°, and −47.36°, respectively. The errors of inclination and declination are both less than 5°, which are all close to the true model. The degree of recovery of the magnetization intensity can reach about 40% for models one and three and about 60% for model two. The overall position of the three models is close to that of the real model, and the maximum value of the magnetization intensity is located inside the model.

From the above results, it can be seen that the method in this paper can obtain good results, even under multiple model conditions with different burial depths and different magnetization intensities, which further validates the applicability of the method.

4. The Origin of the Lunar Weak Magnetic Anomalies

The latest IGRF model shows that the strength of the earth’s magnetic field can reach 65,000 nT [51], while the magnetic field on the surface of the Moon can reach up to about 450 nT [52], which is very weak in comparison. In order to study the origin of lunar weak magnetic anomalies, we selected two different regions, Mare Australia and Mare Crisis. Therefore, we use the fast magnetization vector inversion method in the spherical coordinate system to obtain the three-dimensional magnetic vector structure.

4.1. Mare Australe

The topographic map drawn from the Lunar Orbiter Laser Altimeter topographic data (the data were downloaded from https://planetarymaps.usgs.gov/mosaic/Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.tif (accessed on 21 October 2023)) is shown in Figure 7a, and it can be seen that there are many craters in the region, and the highest point of the topography is 1349 m. The geological map is shown in Figure 7b, which shows the presence of material from three geologic periods, Aitkenian, Imbrian, Nectarian, and Ti-basalts in the region [53]. We used the satellite magnetic anomaly data for the lunar surface solved by Ravat et al. [52] and transformed the data to magnetic anomalies at an elevation of 1349 km, as shown in Figure 7c–f, using the continuation of the data from the curved surface to the plane technique. This revealed that the magnetic anomalies are concentrated near the large craters (the black dashed line).

In addition, combining Figure 7b,c, it can be found that the major part of the anomaly concentration region (within the red line) exists during the Aitkenian geological period, and the material is mainly ejecta, while the blue region is mainly composed of basalt with a relatively weak magnetic field [54]. In addition, the whole region is located within the lunar Mare Australe tectonics.

The underground space is divided into tesseroids with a latitude and longitude of 0.3° and a depth of 1 km, and the number of grids is 101 × 41 × 10, which improves the calculation efficiency by 16 times compared with the conventional method. The magnetic vector structure is shown in Figure 8 and Figure 9a–c, which shows that the highest magnetization magnitude is 0.5 A/m, which, compared to the maximum magnetization intensity of the Earth’s crust, is 12 times smaller [55]. Nicholas et al. [56] studied the magnetic field in the Reiner Gamma region and found that a magnetization intensity of 1 A/m is required if the ejecta is 1 km thick, and Carley et al. [57] calculated a maximum magnetization intensity of 0.2 A/m on the Moon but constructed a model with a minimum layer thickness of 10 km. In our results, the magnetic sources are buried at depths between 4–9 km and a thickness of 6 km, so the magnetization intensity of 0.5 A/m is possible. Dyal et al. [19] found that surface magnetometers show variations in the intensity and the direction of the lunar magnetic field over distances of only a few kilometers, which suggests that they are associated with local magnetic sources close to the surface, indicating that the depths of burial of the magnetic sources we obtain are likely to exist.

In Figure 8 and Figure 9a–c, we can see that the magnetization directions within the magnetic sources are relatively consistent, the magnetization directions between different sources are inconsistent, and the magnetic sources are distributed near the craters. The declinations of magnetic sources Ⅰ, Ⅱ, and Ⅲ in Figure 8 are 342°, 331°, and 93°, respectively, where the closest distances of sources Ⅱ and Ⅲ are only 0.3° in latitude and longitude, respectively, but the angles differ by 122° (or 238°), which indicates that the fast magnetization direction changes in space.

We use the local, isolated, stronger magnetic sources (a magnetization intensity stronger than 0.15 A/m, such as the magnetic sources shown in Figure 9 of the magnetic vector structure) to conduct numerical simulations, and the magnetization directions and intensities inside the magnetic sources are uniform. With magnetic sources, we calculate the B_θ, B_ϕ, and B_r components at different locations in the subsurface space to obtain the simulation results of the angles at different locations. As shown in Figure 9d, we obtain the inclination correlation between the simulation results and the magnetic vector structure in an 89.1°E slice. The results show that the global correlation is relatively high, and the correlation near the magnetic source (the black solid line) can reach more than 75% due to the small influence of other magnetic sources. And, the location of the black dashed line, which is influenced by the magnetic source in the 84.6°E slice, also shows a high correlation. Based on the above results, we believe that the calculated magnetic vector structure conforms to the physical distribution and is relatively accurate.

4.2. Mare Crisium

A topographic map of Mare Crisium from the Lunar Orbiter Laser Altimeter topographic data (the data were downloaded from https://planetarymaps.usgs.gov/mosaic/Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.tif (accessed on 21 October 2023)) is shown in Figure 10a, which shows the presence of a large impact basin, Crisium, with a number of smaller craters in the region. The geologic map in Figure 10b shows that most of the region is covered by basalt. Figure 10c–f is based on the magnetic anomaly data of the lunar satellites solved by Hood et al. [58] at the altitude of 30 km.

The magnetic vector structures in Figure 11 and Figure 12 show that the maximum magnetization intensity in the Mare Crisium region is 0.01 A/m, and the lunar magnetization model of Carley et al. [57] shows a minimum magnetization intensity of 0.0005 A/m, so the magnetization intensity of 0.01 A/m is reasonable.

From Figure 11, it can be seen that the depths of sources I and III are close to each other, ranging from 15–150 km, and the depth of source II ranges from 30–85 km. The inclination direction of magnetic source I is in the range of 45°–60°, that of magnetic source III is in the range of 0°–15°, and that of magnetic source II is around 0°, and the magnetization direction of all three magnetic sources points to the west.

As can be seen in Figure 12, the declinations of the magnetic sources I, II, and III are close to each other, and all of them are oriented in the northwest direction. By analyzing the magnetization direction, it is possible to see that a stable dipole dynamo existed during the formation of the Crisium magnetic anomaly and that other small impact craters may not have produced the magnetic material that caused serious damage to the Crisium impact magnetic field.

5. Discussion

5.1. Causes of Magnetic Anomalies in the Mare Australe Region

Magnetic field simulations by Crawford [39] show that a projectile impact with a velocity of 20 km/s and an orientation of 45° from a horizontal standpoint can produce an impact crater with a diameter of 340 km and a maximum near-surface residual field of only 15 nT. They additionally show that, even in the region of high magnetization, it must be 10 km thick to produce an observable field at orbital altitude. The craters in our study area are all much smaller than this diameter but all produce magnetic fields much larger than 15 nT, and our source is not 10 km thick, so the impact field may not be responsible for the magnetic field in this region.

As can be seen in Figure 8 and Figure 9a–c, the magnetization direction is relatively consistent within the magnetic source in the Mare Australe region, and the magnetization directions between different magnetic sources are inconsistent. With regard to the differences of magnetization directions between different magnetic sources on the Moon, the lunar dynamo theory believes that it may be caused by the polar wander of the ancient geomagnetic pole. Arkani-Hamed [59] explained the polar wander of the Moon. He suggested that the Moon experienced two different polar wanders: one occurred during the Aitkenian period (4.3–3.92 Ga), which moved ~68°, and another event occurred after the Aitken period, when the formation of a large Mascon basin forced the Moon to experience a further ~52° polar wander. However, in our results, the declinations between magnetic sources II and III with a latitude and longitude distance of only 0.3° differ by 122° (or 238°) and are both in the Aitkenian period, which is a large difference from the polar wander of 68°.

Takahashi and Tsunakawa [60] found through numerical simulations that thermal convection is dominated by westward zonal flow in the farside part of the Moon with no radial flow. The westward flow in the farside part and the eastward flow in the nearside part converge to form frequent and strong radial flows. This peculiar flow structure generates a magnetic field, and because of the complex convection phenomena that is generated, it might result in a very high frequency of lunar polarity reversal. In Figure 8, we can also observe that the magnetic source is located near the boundary between the far and near sides of the moon, which may be affected by frequent polarity reversal and could result in differences in magnetization direction. If the cooling timescale of the whole object is long enough to capture one or more geomagnetic reversals, then the high reversal rate could easily lead to an inconsistent direction of magnetization in the two near-neighboring sources. The simulation also found that the Moon has relatively weak magnetization, and we also found through inversion that the magnetization in this area is much lower than that of the Earth. In addition to this, intermittent or episodic dynamo mechanisms are also possible causes of magnetization direction inconsistencies, such as impact-driven precession dynamos [61] or founding dynamos of cold accretion to the lunar core–mantle boundary after magma–ocean solidification [62]. Due to the short lifespan of this type of dynamo, their behavior may change significantly, even within the cooling timeframe of a single region of the lunar crust, which may be an alternative explanation for the Australe magnetic field. In addition, due to the existence of many basaltic rocks in the region, which suggests the existence of intense magmatic activity, we suggest that the shallow magnetic sources formed in the region are related to igneous intrusions.

Figure 9d shows that the correlation between the simulated magnetization directions and the actual magnetization directions is high in and near the magnetic sources, with a correlation of more than 75%, indicating that the inversion yields a more accurate distribution of magnetization directions. Combined with previous studies and discoveries as well as with the results of this paper, we suggest that weak magnetic anomalies in the lunar crust in the Mare Australe region are likely caused by the core dynamo and are not related to meteorite impacts.

5.2. Causes of Magnetic Anomalies in the Mare Crisium Region

Baek et al. [63] calculated Crisium’s magnetic anomaly and showed that the deepest sources can be up to 55–70 km, whereas our results can be up to 150 km deep, suggesting that Crisium’s magnetic sources exist at deeper locations. Experimental results from Crawford [39] showed that impact craters with diameters of 540 km produced by impacts can generate magnetic fields at a subsurface of 200 km in depth to generate a magnetic field. However, the depths of magnetic sources Ⅰ and Ⅲ in Figure 10 are close to each other, ranging from about 15–150 km, and the depths of magnetic source Ⅱ range from 30–85 km, and there are depth differences between the sources. Baek et al. [63] suggested that different depths may cause differences in cooling times and record different magnetic fields and that the impact magnetic field is difficult to record at different depths simultaneously because it is generated transiently. From this point of view, the magnetic anomaly of the Crisium impact crater may originate from the dynamo magnetic field.

As can be seen in Figure 11 and Figure 12, the declinations of the three sources are very close to each other, suggesting that a dipole dynamo may have existed, and the change in inclination may be due to the different positions and depths.

6. Conclusions

We developed a fast magnetization vector inversion method for the inversion of magnetic data with an undulating observation surface in a spherical coordinate system. The theoretical model test proves that the proposed method can effectively recover the magnetic vector model, and the computational efficiency of this paper’s method is 23 times higher than the conventional method with a 60 × 60 × 10 grid division. Furthermore, the improvement of inversion efficiency is more obvious with the increase in the number of grid divisions.

The three-dimensional magnetic vector structures of the Mare Australe region (27°S–39°S, 75°E–105°E) and the Mare Crisium region (2°N–30°N, 45°E–70°E) are obtained by the developed fast magnetic vector inversion method for undulating terrain in a spherical coordinate system, which can provide an analytical basis for the study of the origins of lunar magnetic anomalies.

In the Mare Australe region, the magnetization vector results indicate that the depths of the magnetic sources are less than 10 km, probably due to igneous intrusion. In addition, the magnetization direction of a single magnetic source is relatively concentrated. However, the magnetization directions of adjacent magnetic sources are not consistent, and there are significant differences between certain magnetic sources (differences in declination of 122° or 238°), which may be related to the rapid reversal of the geomagnetism as well as the short period generators. Based on the above, we suggest that the magnetic anomaly in the Mare Australe region is formed by dipole dynamo.

In the Mare Crisium region, the magnetic sources are buried deeper, and there are differences in depths between the different sources. This phenomenon suggests that the magnetic sources at different depths are formed at different times, which is difficult to explain by the transient magnetic field generated by the impact. Moreover, the declinations of different sources are close to each other, which indicates that a section dipole dynamo existed during the formation period of the sources in this region, and no polar drift occurred. We believe that the magnetic anomalies in the region are origins of generator magnetic fields.

In addition to the above characteristics, regardless of the depth of the magnetic sources in the Mare Australe or the Mare Crisium regions, the magnetization directions of the magnetic sources are uniformly distributed, which does not reflect the dispersion of the magnetization direction of the meteorite impact magnetic field.

In conclusion, our method is effective for the construction of large region magnetic vector structures, and the study of magnetic anomalies in two regions of the Moon reveals that the generator magnetic field may be the main origin of the current lunar magnetic anomaly.