An Improved 3D Magnetization Inversion Based on Smoothness Constraints in Spherical Coordinates

Abstract

In the inverse problem, the traditional way to obtain a stable solution is based on the maximum smoothness criteria. However, this approach cannot generate clearer and more focused images. In this study, we propose an improved inversion method based on the smoothness constraints. In the algorithm, the model weighting functions are updated by adding a model’s total gradient module matrix, which can effectively constrain the boundary of the recovery model in the iterative operation. We invert the 3D magnetization intensity for the three-component magnetic data in the spatial domain by spherical coordinates. The preconditional conjugate gradient algorithm is introduced to improve the efficiency of the solutions. We design two sets of synthetic examples to evaluate the inversion effects, which show that the improved method is more reliable than the smoothness constraint method. The boundary of the magnetic bodies is more precise, and the magnetization ranges are more focused. The method does not rely on the initial model and is suitable for magnetic vector data inversion. We also apply the algorithm to a set of Dabie orogen three-component magnetic data derived from a geomagnetic field model and verify the effectiveness of the inversion method.

1. Introduction

The successive launches of CHAMP, GRACE, GOCE, Swarm, and other potential field-observing satellites have brought unprecedented opportunities and challenges to the processing and applications of gravity and magnetic data. Several high-order geomagnetic field models were developed, such as NGDC-720 and EMM2015, compiled from different data sources (satellite, oceanic, aeromagnetic, and geomagnetic survey data). The magnetic data derived from the models are one of the powerful instruments for the regional tectono-geophysical investigation of the Earth’s crust and upper mantle [1], and are being increasingly mapped for large regions geological analysis [2,3,4,5,6].

Another application for magnetic data is 3D physical property (magnetization or susceptibility) inversion, which visually outlines the spatial shapes and distributed features of magnetic field sources [7,8,9,10]. As one of the main physical property inversion methods, a considerable amount of in-depth research has been conducted regarding 3D magnetization intensity inversion techniques. Wang et al. [11] proposed a total magnetization vector inversion method to invert 2D layered models, and obtained the magnitude and direction of the magnetization. In another related study, Liu et al. [12] carried out magnetization vector inversion research on 2D borehole magnetic survey data. Furthermore, various researchers have determined the three-dimensional distribution features of magnetization based on magnetic anomaly amplitude data [13,14]. The magnetization intensity inversion methods of those spatial domains are generally carried out in Cartesian coordinates. They are not directly applicable to the inversion of global magnetic data, which are inherently in the spherical coordinates. The inversion methods extended to the spherical coordinates provide technical support for regional physical property inversion [8,15,16,17,18]. It is necessary to develop an applicable magnetization intensity inversion method to adapt to the rapidly increasing satellite-derived magnetic data.

The inverse problem in spherical coordinates is similar to that in Cartesian coordinates. One of the difficulties that the inverse problem faced was the non-uniqueness of the model. Many methods have been proposed to solve this problem in the Cartesian coordinates. These methods included smooth and small model inversion [19,20], focusing inversion [21,22], cokriging [23,24], binary and multinary inversion [25,26], and c-means clustering [27]. The available methods showed that they could decrease the non-uniqueness of models. Additionally, one possible solution to solve the problem is to adopt the Occam inversion method [28], which is developed based on Tikhonov regularization [29]. It is less dependent on the initial models and is now widely used in gravity and electromagnetic research [30,31]. This method can determine the smoothest solution for multi-layered geophysical models. The roughness of a model is minimized under a certain standard of a misfit. Due to the inherent properties of the method, the border of recovered models is not clear, and the images are not focused, which may be inapplicable under some geological issues. We propose an improved inversion method to improve the situation and describe it in detail later.

The magnetic data derived from geomagnetic field models are the magnetization responses generated at the lithospheric scale. Little is known about the remnant magnetization structures in the Earth’s continents. The induced magnetization and remnant magnetization in the magnetic anomalies of the continental crust cannot be decoupled/distinguished effectively [32,33]. Compared with the oceanic crust, the induced magnetic anomalies in the continental crust usually play major roles, and the remnant magnetization is usually ignored [34,35,36]. Moreover, due to many observation grids, field source space cells, and magnetic parameters, the regional magnetization intensity inversion in the spherical coordinates will be very time-consuming. Consequently, it can be considered that the lithospheric magnetic field is mainly dominated by induced magnetization to simplify the parameter space of the models and actual data.

The remainder of this paper is as follows: we first review the forward modeling and inversion of the magnetic problem in spherical coordinates and propose an improved inversion method based on the smoothness constraints. The third section evaluates the method with two synthetic examples of regional models. We apply the method to the magnetic vector data of the Dabie orogen to test the actual application effect in the fourth section. Finally, in the fifth section, we present the conclusions of this study.

2. Methodology

2.1. Forward Modeling in Spherical Coordinates

Different calculation methods can obtain magnetic potential field responses. In space domain magnetic vector forwarding methods, the magnetic potential and vector fields were usually directly calculated based on the magnetic dipole effects [37,38,39]. Blakely [40] proposes an alternative to link the magnetic potential and vector fields with the gravitational fields, referred to as Poisson’s relation. The two means were determined to be mathematically equivalent. We adopted the second method to obtain the magnetic vector field. We have to consider the gravitational potential and related derivatives calculations for the first step.

The gravitational potential $V\left(P\right)$ of a solid body with volume $\mathrm{\Omega }$ and density $\rho$ at the observation point $P$ outside of $\mathrm{\Omega }$ is given by Newton’s integral [41]

$$V\left(P\right)=G\rho {\displaystyle \underset{\mathrm{\Omega }}{\iiint }\frac{1}{\ell }d\mathrm{\Omega }}$$

(1)

The numerical evaluations of Equation (1) rely on mass discretization. It is needed to discretize the model’s underground space into a specific geometric shape and density. Using prisms is especially beneficial for local applications in Cartesian coordinates. For larger application areas, the Earth’s curvature has to be taken into account. It has been determined that a tesseroid model (Figure 1) with a constant height delineated by geographic grid lines, which are directly linked to the curvature of the Earth, is an ideal volume unit for forwarding.

The gravity potential based on tesseroid model has no analytical solution; it must be calculated approximately by integral method. There are mainly two kinds of algorithms: Gauss Legendre Quadrature (GLQ) and Taylor series expansion method based on integral kernel function. Ku [42] determines the distribution of anomalies using the least squares method based on equivalent point- mass. The interference point- mass in the anomaly body can be estimated according to the GLQ in the gravity field expression. He proposes an empirical standard in the method, that is, the distance between the points is not greater than the distance to the calculated point, and the accuracy of the solution basically remains unchanged. M. F. Asgharzadeh et al. [43] give an efficient and high-precision calculation method for gravity potential, gravity field, and gravity tensor of spherical unit body under the standard of Ku [42]. The minimum number of nodes can be given during calculation, and the distance to the observation points is greater than the nodes spacing, which reduces the number of calculation units and improves the efficiency. However, when the calculation point is close to the tesseroid element, the solutions of GLQ are unstable. Heck and Seitz [44] compare the gravity field calculation approaches of tesseroid, prism and point-mass. The third-order Taylor series expansion of the integral kernel function is carried out to approximate the gravity potential and its radial derivative. The results show that the zero-order form is equivalent to the formula of the point- mass, and the calculation efficiency and accuracy of the tesseroid are higher than those of the other two elements. To some extent, tesseroid has advantages in calculating the gravitational field.

A tesseroid is composed of radial, latitude, and longitude $\left(r,\varphi ,\lambda \right)$, which was used to calculate the potential field and gradient tensors in the spherical coordinates. The gravitational potential, gravitational acceleration, and Marussi tensor, respectively, generated by the homogeneous body at the observation point $P$ were described by

$$\left\{\begin{array}{c}V\left(P\right)\\ g_{\alpha }\left(P\right)\\ g_{\alpha \beta }\left(P\right)\end{array}\right\}=G\rho {\displaystyle {\int }_{{\lambda }_1}^{{\lambda }_2}{\displaystyle {\int }_{{\varphi }_1}^{{\varphi }_2}{\displaystyle {\int }_{r_1}^{r_2}\frac{1}{{\ell }^3}\left\{\begin{array}{c}{\ell }^2\\ {\Delta }_{\alpha }\\ 3{\Delta }_{\alpha }{\Delta }_{\beta }{\ell }^{-2}-{\delta }_{\alpha \beta }\end{array}\right\}\kappa d{r}^{\prime }d{\varphi }^{\prime }d{\lambda }^{\prime }}}}$$

(2)

where

$\begin{array}{l}{\Delta }_x=r^{\prime }\left(\cos \varphi \sin {\varphi }^{\prime }-\sin \varphi \cos {\varphi }^{\prime }\cos \left({\lambda }^{\prime }-\lambda \right)\right)\\ {\Delta }_y=r^{\prime }\cos {\varphi }^{\prime }\sin \left({\lambda }^{\prime }-\lambda \right)\\ {\Delta }_z=r^{\prime }\cos \psi -r\end{array}$

$\alpha ,\beta \in \left\{x,y,z\right\}$, in which $x,y,z$-, respectively, represent the north- east- up direction of the local north-oriented frame (LNOF), $\ell =\sqrt{{r^{\prime }}^2+r^2-2r^{\prime }r\cos \psi }$, $\kappa ={r^{\prime }}^2\cos {\varphi }^{\prime }$, and ${\delta }_{\alpha \beta }$ are Kronecker’s Delta, and $\cos \psi =\sin \varphi \sin {\varphi }^{\prime }+\cos \varphi \cos {\varphi }^{\prime }\cos \left({\lambda }^{\prime }-\lambda \right)$.

We adopted the Gauss- Legendre quadrature (GLQ) integral to calculate the gravity field and gradient tensors [44,45] and truncated the order of each component of the GLQ to 2 to improve the calculating efficiency.

According to Poisson’s relation [40], an object with uniform magnetization and uniform density has the following relation between the magnetic vector $B$ and Marussi tensor $T$

$$B=-C_{m}T\cdot M$$

(3)

where

$T=\left[\begin{array}{lll}{g}_{xx}\left(r_P,{\varphi }_P,{\lambda }_P\right)\hfill & g_{xy}\left(r_P,{\varphi }_P,{\lambda }_P\right)\hfill & g_{xz}\left(r_P,{\varphi }_P,{\lambda }_P\right)\hfill \\ g_{yx}\left(r_P,{\varphi }_P,{\lambda }_P\right)\hfill & g_{yy}\left(r_P,{\varphi }_P,{\lambda }_P\right)\hfill & g_{yz}\left(r_P,{\varphi }_P,{\lambda }_P\right)\hfill \\ g_{zx}\left(r_P,{\varphi }_P,{\lambda }_P\right)\hfill & g_{zy}\left(r_P,{\varphi }_P,{\lambda }_P\right)\hfill & g_{zz}\left(r_P,{\varphi }_P,{\lambda }_P\right)\hfill \end{array}\right]$

the constant $C_m={\mu }_0/\left(4\pi G\rho \right)$, in which $G=6.674\times {10}^{-11}{\mathrm{m}}^3\cdot {\mathrm{kg}}^{-1}{\mathrm{s}}^{-2}$ is the gravitational constant and ${\mu }_0=4\pi \times {10}^{-7}\mathrm{H}\cdot {\mathrm{m}}^{-1}$ is the magnetic permeability of the free space, $M={\left[\begin{array}{ccc}{M}_x& M_y& M_z\end{array}\right]}^T$ refers to a matrix or vector composed of three orthogonal components of total magnetization intensity $M$.

The relations among each component of the $M$, magnetic inclination $I$, and magnetic declination $A$ were as follows

$$\left\{\begin{array}{l}{M}_x=M\cos I\cos A\\ M_y=M\cos I\sin A\\ M_z=M\sin I\end{array}\right.$$

(4)

Here, when only considering the induced magnetization of the magnetic body, $I$ indicates the local geomagnetic inclination, and $A$ is the angle between the direction of the magnetization direction and the magnetic north.

Based on the principle of superposition, the effect of the whole magnetization distribution can be approximated by the sum of the impact over all individual magnetic bodies. The parameterized forward modeling problem is expressed as

$$d=Gm$$

(5)

where, $d$ is the $n$ dimensional data vector of observed magnetic anomaly, $m$ is the $m$ dimensional magnetization vector, and $G$ refers to the nonlinear function that maps parameters $m$ to magnetic field data $d$, which is expressed as the forwarding kernel matrix of the $n\times m$ size.

2.2. Occam Inversion Method

The general form of constructing the objective function and minimizing it can be described as

$$\min :\mathrm{\Gamma }\left(m\right)={\varphi }_d\left(m\right)+\lambda {\varphi }_m\left(m\right)$$

(6)

where the Lagrange multiplier $\lambda$ is the regularization factor of the preconditioned conjugate gradient (PCG) algorithm to balance the data misfit function ${\varphi }_d\left(m\right)$ and the model objective function ${\varphi }_m\left(m\right)$.

We introduced a weighting function, and the re-weighted data misfit function was expressed as the ${\ell }_2$ norm of the observed and predicted data

$${\varphi }_d\left(m\right)={\left(\Delta d-G\Delta m\right)}^T{W}_d^T{W}_d\left(\Delta d-G\Delta m\right)$$

(7)

where $\Delta d$ is the corresponding modification of the observed data column vector, $\Delta m=m-m_{ref}$ represents the modification of the model parameter vector $m$ to the initial reference model $m_{ref}$, and $W_d=diag\left\{\left.1/{\sigma }_1,1/{\sigma }_2,\cdots ,1/{\sigma }_m\right\}\right.$ refers to the data space weighted matrix, which is a normalized diagonal matrix representing the standard deviation of the $i-th$ observation data.

In order to ensure that the model is smooth in three spatial directions, it is necessary to select a suitable model objective function to impose constraints on the inversion model. Li and Oldenburg [19] and Li [46] designed the model objective functions of the maximum smoothing method. However, those functions were only applicable to the Cartesian coordinates. Du et al. [8] and Liang et al. [15] introduced the smoothest model objective functions into the spherical coordinates. Considering the offset effects of the recovered models relative to the reference models on the sphere, we introduced a weighting function to constrain the models smoothly. Constructing the model objective functions as

$$\begin{array}{l}{\varphi }_m\left(m\right)={\alpha }_s{\displaystyle {\int }_v{\left[w\left(r\right)\left(\Delta m\right)\right]}^2}dv+{\alpha }_r{\displaystyle {\int }_v{\left[\frac{\partial w\left(r\right)\left(\Delta m\right)}{\partial r}\right]}^{2}dv}\\ \begin{array}{ccc}& & +\end{array}{\alpha }_{\varphi }{\displaystyle {\int }_v{\left[\frac{\partial w\left(r\right)\left(\Delta m\right)}{r\partial \varphi }\right]}^{2}dv}+{\alpha }_{\lambda }{\displaystyle {\int }_v{\left[\frac{\partial w\left(r\right)\left(\Delta m\right)}{r\cos \varphi \partial \lambda }\right]}^{2}dv}\end{array}$$

(8)

where ${\alpha }_s$ represents the weighting function of the minimum model constraint, $\left({\alpha }_r,{\alpha }_{\varphi },{\alpha }_{\lambda }\right)$ is the relative weighting of the smoothest model objective function along three different directions on the sphere. In particular, $w\left(r\right)$ indicates the radial weighting function.

In the unconstrained inversion of the potential field data, the kernel functions are attenuated sharply with increased depth. It will cause the recovered model to be concentrated near the surface. Li and Oldenburg [19,20] first introduced the depth weighting function $w\left(z\right)={\left(z+z_0\right)}^{-\beta /2}$ into the model constraint terms of 3D gravity and magnetic inversions to overcome the attenuation effect. They applied it in 3D magnetic susceptibility imaging. Du et al. [8], Liang et al. [15], and Wang et al. [16], respectively, extended the depth weighting function to spherical coordinates and proposed different radial weighting functions. We used the radial weighting function of Wang et al. [16] as follows

$$w\left(r\right)=\frac{1}{{\left(H+R-r\right)}^{-\frac{\beta }{2}}}\cdot \sqrt{\frac{r}{R}}$$

(9)

where $R=6371.2$ km is the average radius of the Earth, $H$ is the average height of the observation point, and $r$ refers to the radial distance from the center of the model cell to the center of the Earth.

After that, the finite difference was substituted for the partial differential to discretize the model objective function, which could be described as a matrix

$$\begin{array}{ll}{\varphi }_m\left(m\right)& =\Delta m^T\left(W_s^T{W}_s+W_r^T{W}_r+W_{\varphi }^T{W}_{\varphi }+W_{\lambda }^T{W}_{\lambda }\right)\Delta m\\ & =\Delta m^T{W}_m^T{W}_m\Delta m\end{array}$$

(10)

where $W_m={\alpha }_m{R}_{m}D$ represents the weighted matrix after discretization: ${\alpha }_m$ is the weighting coefficient of each item in the model objective function, $R_m$ indicates the differential operator along each direction of $(r,\varphi ,\lambda )$, and $D$ refers to the discrete matrix of the radial weighting function.

The data misfit and the model objective function were substituted into the objective function, which was expressed as

$$\mathrm{\Gamma }\left(m\right)={\left(\Delta d-G\Delta m\right)}^T{W}_d^T{W}_d\left(\Delta d-G\Delta m\right)+\lambda \Delta m^T{W}_m^T{W}_m\Delta m$$

(11)

The expression could be written as

$$\left[\begin{array}{c}G\\ \sqrt{\lambda }{W}_m\end{array}\right]\Delta m=\left[\begin{array}{c}d\\ 0\end{array}\right]$$

(12)

Let $A=\left[\begin{array}{c}G\\ \sqrt{\lambda }{W}_m\end{array}\right]$ represent the Jacobian Matrix $b=\left[\begin{array}{c}d\\ 0\end{array}\right]$, and Equation (12) could be simplified to

$$A\Delta m=b$$

(13)

The inversions of three-component magnetic anomalies come down to solving Equation (13), whose coefficient matrices are symmetric and positive-definite. The convergence speed of Equation (13) is determined by the condition numbers of the coefficient matrices. However, the condition numbers of the coefficient matrices are generally huge, which may seriously affect the computational efficiency of iterative convergence. This problem can be solved by indirect (iterative) methods such as the conjugate gradient (CG) method [47]. It has been found that the CG method is very effective in solving large-scale linear equations and is widely used in geophysical inversion processes [48,49,50]. The CG method combines the gradient direction and the conjugate direction. It searches for the optimal solution along a set of conjugate directions constructed by the gradient of the initial point.

We adopted the PCG algorithm in the matrix operation, which was developed based on the CG method. The preconditioned matrix was utilized to improve the condition numbers of the equations, and both ends of the equation were multiplied by the preconditioned matrix [48]. By implementing the approaches mentioned above, the eigenvalues of Equation (13) were concentrated along the diagonal, thereby improving the iteration efficiency.

2.3. Improved Inversion Method

In the regularization method, the general scheme to obtain the solution under certain conditions is to constrain the regularization terms. Many model constraint methods are proposed to obtain a suitable solution. These methods include minimum norm constraint [29], compact constraint [51], minimum moment of inertia constraint [52], minimum gradient constraint [21], and convexity constraint [53]. In this study, we need to obtain a stable solution and expect to effectively describe the boundary of the geological body.

In the minimum moment of the inertia constraint method, the weighting matrix $W_m$ is made of a reference model that model parameters are headed toward it during the iterative inversion process. Inspired by this method, we consider constraining and processing the weighting matrix of Occam’s model objective function. However, this process is not carried out directly on the weighting matrix. The first step is to calculate the total gradient module of the reference model

$$\widehat{m}=\sqrt{{\left(\frac{\partial m}{\partial r}\right)}^2+{\left(\frac{\partial m}{r\partial \varphi }\right)}^2+{\left(\frac{\partial m}{r\cos \varphi \partial \lambda }\right)}^2}$$

(14)

The updated weighting matrix is expressed as the inner product of two matrices

$${\widehat{W}}_m=W_m{K}_m$$

(15)

where $K_m=diag\left\{{\left.{\widehat{m}}_1,{\widehat{m}}_2,\cdots ,{\widehat{m}}_m\right\}}^{-\gamma }\right.$ is a diagonal matrix.

Finally, Equation (15) was substituted into Equation (11) for iterative calculation.

3. Synthesis Model Analysis

3.1. Models’ Geometry and Magnetic Parameters

Two sets of synthesis models were utilized to test the effectiveness of the inversion method. The first set included two hexahedron models with the same magnetic parameters (Models A and B). The second set included four hexahedron models with different magnetic and geometric parameters (Models C, D, E, and F). The model’s parameters and 3D visualization are shown in

Table 1. Synthesis models’ geometry and magnetic parameters.

Model Id	Longitude	Latitude	Depth	Inclination	Declination
A	114.0–115.0	30.1–30.9	30–100	45	45
B	116.0–117.0	30.1–30.9	30–100	45	45
C	114.2–114.8	31.2–31.8	30–70	45	45
D	116.0–117.0	31.3–31.7	40–100	60	45
E	114.1–114.9	29.5–30.3	40–90	45	60
F	116.3–116.7	29.3–30.3	30–70	90	30

Note: In the table, the units of longitude, latitude, inclination, and declination are shown in degrees; the units of depth are in km.

and Figure 2. The height of forwarding calculations was set as 4 km, and the observational grid size measured $0.1°\times 0.1°$. During the inversion process, the field source spaces were subdivided into tesseroid models, and the sizes of the subdivided cells were $0.1°\times 0.1°\times 5$ km. The total magnetization intensity is 10$A\cdot m^{-1}$.

3.2. Analysis of Inversion Results

3.2.1. Synthesis Model of the Same Magnetic Parameters

Figure 3 describes the three-component magnetic anomalies of the first set of models. Figure 4 and Figure 5 detail the magnetization imaging results of the Occam method and the improved method, respectively. The inversion results of the Occam method indicate that the magnetization boundary is fuzzy, and the field source range is difficult to define. The imaging results of the improved method can fit the actual model better. The magnetization ranges are relatively focused, and the field source’s boundaries are clear.

The two methods are both robust in the presence of noise. Some tiny false magnetic sources appeared near the surface of $\Delta B\mathrm{x}-$ component and $\Delta B\mathrm{y}-$ component imaging results.

It should be noted that the values obtained by both inversion methods deviate several times from the design models. The main reason is affected by the volume effect. Other reasons are that the recovered magnetic sources are still not focused enough (especially in the depth direction), and the model space is not finely divided. In comparison, the imaging results of the improved method are closer to the actual values.

Figure 3. Three-component magnetic anomalies, $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$, of the first set of models: (a–c) cases with noise-free; (d–f) cases with 5% noise signal. The black borders represent the position of the actual models.

Figure 3. Three-component magnetic anomalies, $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$, of the first set of models: (a–c) cases with noise-free; (d–f) cases with 5% noise signal. The black borders represent the position of the actual models.

Figure 4. The inversion results (the profile at latitude = 30.5°) for three-component magnetic anomalies using the Occam method: (a–c) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$ of noise-free; (d–f) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$ with 5% noise signal.

Figure 4. The inversion results (the profile at latitude = 30.5°) for three-component magnetic anomalies using the Occam method: (a–c) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$ of noise-free; (d–f) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$ with 5% noise signal.

Figure 5. The inversion results (the profile at latitude = 30.5°) for three-component magnetic anomalies using the improved method: (a–c) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$ of noise-free; (d–f) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$ with 5% noise signal.

Figure 5. The inversion results (the profile at latitude = 30.5°) for three-component magnetic anomalies using the improved method: (a–c) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$ of noise-free; (d–f) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$ with 5% noise signal.

3.2.2. Synthesis Model of Different Magnetic Parameters

We still need to test the anti-noise ability of both inversion algorithms in the second set of models. The 5% noise signal was added to the magnetic three-component anomalies, and the noise levels are as shown in Figure 6, Figure 7 and Figure 8 are the magnetization imaging results of the Occam method and the improved method, respectively. In both figures, (a) to (f) are the vertical slices of the latitude directions, and (g) to (l) are the vertical slices of the longitude directions.

Comparing the imaging effects of the two inversion methods, the inversion results based on the improved method are in better agreement with the designed actual models. The ranges of magnetization distribution are small, and the properties boundaries are clear. The imaging result of the Occam method has a lower resolution, especially in the depth direction.

Similar to the conclusions of the first set of models, the magnetization intensity values of the actual model are several times higher than those obtained by the two inversion methods. The inversion results of the improved method are closer to the actual model. Noise has little effect on the $\Delta Bz$-component.

Overall, the improved inversion method clarifies the recovered magnetic model boundary and improves the depth resolution. For the three-component magnetic data, the inversion results of the vertical component are relatively stable.

In the model test, we obtained a group of weighted parameters ($\beta$ takes 3, $\gamma$ is 0.5) with good results, and used the same parameter values for the following magnetic data inversion.

4. Magnetic Data Tests

The Dabie orogen is located in the middle section of the Qinling- Dabie- Sulu orogen in China, which is formed by the subduction of the northern margin of the South China Craton during the Triassic period and the continental collision with the North China Craton [54,55]. It is one of the largest high-pressure (HP) and ultrahigh-pressure (UHP) metamorphic rock belts worldwide. The orogen is composed of a series of metamorphic units bounded by faults (Figure 9). Tectonically, the NE- SW striking fault zone F6 divides the orogen into eastern Dabie and western Dabie. Four metamorphic units are further divided in eastern Dabie: northern Huaiyang (NHY) tectonic belt, northern Dabie complex unit, UHP metamorphic belt, and HP metamorphic belt. The NHY and HP belts are two relatively low-grade metamorphic rock units. Within the orogen, the rock types are mainly lower-degree to higher-degree metamorphic rock.

Many deep seismic investigations, electromagnetic explorations, and experimental gravity studies have been carried out to reveal the deeper physical properties information of the lithosphere. Seismic refraction profile results [57] show that the crust can be divided into upper, middle, and lower layers according to the velocity difference, and the thickness of the crust gradually increases from 35km in the south to 41km in the north. From the results of gravity spectrum analysis [58], there is obvious lateral heterogeneity in the middle and upper crust of the orogen. The inversion result of the magnetotelluric profile [59] shows that there are several isolated high-conductivity anomaly bodies in the middle crust. To sum up, there are low-density, high-resistivity, and low-velocity anomalies in Dabie orogen. However, researches regarding the lithospheric magnetic structures is still lacking. Yang and Li [58] previously used continuous wavelet transform to delineate two magnetic sources based on the total field magnetic anomalies along the profile line (the white line of dashes in Figure 10), which provided a basis for the comparative study in the area.

The lithospheric three-component magnetic anomalies of the Dabie orogen have been calculated at an aeromagnetic altitude (Figure 10). The magnetic survey datasets came from EMM2017 (Enhanced Magnetic Model, data can be found online: https://www.ngdc.noaa.gov/geomag/EMM/index.html, accessed on 1 September 2022). Spherical harmonic functions of the model are up to the orders of 790; that is, the spatial resolution is about 0.25 degrees. Moreover, the resolvable spatial wavelength is as fine as 51 km. The method used to calculate the lithosphere’s magnetic anomalies included stripping the core magnetic field (orders $n\le 15$) from the total magnetic field [60]. The observational grid size measured $0.1°\times 0.1°$ and the size of each model space cell was $0.1°\times 0.1°\times 5$ km in the applications of the forward modelling and inversion.

Figure 11 and Figure 12 show the vertical slices of the magnetization inversion results for the three-component magnetic anomalies of the Dabie orogen using Occam and the improved method, respectively.

There are three magnetic sources in the profile. However, the Occam method cannot determine the depth and extension of the magnetic field sources. In particular, the inversion results for $\Delta B\mathrm{x}-$ and $\Delta B\mathrm{y}-$ components data do not focus, even diverge on the boundary.

The improved method obtains a better inversion effect. We can clearly outline the contours of the three magnetic sources in Figure 12. The two magnetic sources below the surface near 30° N and 31° N, respectively, are consistent with the depth and horizontal extensions of the results of Yang and Li [58]. Comparing the inversion effects of the three-component magnetic data, it can be seen that the $\Delta B\mathrm{z}-$ component is more stable.

The formation of UHP rocks indicates that the crustal rocks were subducted to depths of more than 100 km and then rapidly exhumed [61,62]. The upper part of the obvious magnetic source in Figure 12 is the F5 fault. Deep seismic wide-angle reflection/refraction studies have revealed that the Moho is cut off at 41 km below the F5 fault [63]. It may indicate the location of the magnetic source intrusions. The mafic-ultramafic intrusions exposed along the F5 included both alpine-type peridotite and pyroxene-gabbro intrusions produced by crust-mantle interactions. These are products of syn-collision intrusion processes [64], which may be related to the generation period of the magnetic source. At present, we know little about the information of magnetic structure, and we need to do further research in the future.

Figure 11. The inversion results for three-component magnetic anomalies of the Dabie orogen using the Occam method: (a–c) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$, respectively.

Figure 11. The inversion results for three-component magnetic anomalies of the Dabie orogen using the Occam method: (a–c) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$, respectively.

Figure 12. The inversion results for three-component magnetic anomalies of the Dabie orogen using the improved method: (a–c) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$, respectively.

Figure 12. The inversion results for three-component magnetic anomalies of the Dabie orogen using the improved method: (a–c) the inversion results for $\Delta B\mathrm{x}$, $\Delta By$, and $\Delta Bz$, respectively.

5. Conclusions

We introduced the Occam inversion method into spherical coordinates and inverted the 3D magnetization intensity for the three-component magnetic anomalies. We designed two sets of synthetic models to test the inversion effects and noise impact. The results indicated that the Occam method could reflect the approximate position of the magnetic sources. At the same time, the boundary could not be delineated, and the resolution in the depth direction was low. Under the influence of noise, the inversion result of the $\Delta B\mathrm{z}-$ component is relatively stable.

We proposed an improved inversion method and compared the imaging effects of the improved method and the Occam method on the two sets of synthetic models. The results showed that the developed method improved the focus and depth resolution.

The magnetization amplitudes of both algorithms deviated from the actual model, mainly caused by the volume effect of the spherical coordinates. We applied the improved method to the magnetic data of the Dabie orogen. The consistency of these results with those obtained from the continuous wavelet transform method verifies the newly developed algorithm.

Furthermore, to improve the computational efficiency, the observed grids or field source space cells of the forward modeling and inversion subdivision were still not sufficiently fine. It will be necessary to develop a fast algorithm suitable for larger regions or large-scale inversion in the future.