The non-uniqueness and instability of the solutions have typically focused on previous magnetic anomaly inversion studies, as they are important factors restricting magnetic inversion development. Tikhonov and Arsenin (1977) proposed the regularization theory to solve such problems. The forward sum function, A, decreases rapidly with increasing depth, yielding an inversion result near the ground (Li and Oldenburg, 1998; Zhdanov, 2002). The depth weighting of this regularization method depends on the selection of the regularization factor, α. The adaptive regularization factor selection method proposed by Hu et al. (2019) was applied in this study. We followed Li and Oldenburg (1996), using the reciprocal of the vertical surface component of the total magnetic anomaly data as the data-weighting factor. If the model sensitivity matrix is directly incorporated into data-fitting (Zhdanov, 2002; Hu et al., 2019; Zhao et al., 2020), this problem can be solved. After model weighting, the total sensitivities of the weighted model parameter of all grid cells are equal, and the contribution of observed data is essentially identical. The objective function is defined as follows: P α ( m w ) = Φ ( m w ) + α S ( m w )   = ‖ W d ( A w m w − d ) ‖ 2 + α ‖ ( m w − m a p r w ) ‖ 2 (1) where ϕ(m ^w) is the data-fitting function and S (m ^w) is the model-stabilizing function.

The model weighted matrix (W _m) is: W m = d i a g ∑ A 2 = d i a g ( A T A ) 1 / 2 (2) where m a p r w is the weighted prior reference model, A w = AW m − 1 is the weighted forward operator, and m w = W m m is the weighted model parameter.

The uniqueness of the 3-D inversion of the total-field magnetic anomalies is the forward operator, A, with different magnetic inclination and declination in different grids (Bhattacharyya, 1964). A ( x , y , z ) = μ 0 4 π { k 1 ⋅ l n [ r + ( ξ − x ) ] + k 2 ⋅ l n [ r + ( η − y ) ] + k 3 ⋅ l n [ r + ( ζ − z ) ] + k 4 ⋅ a r c t a n ( ξ − x ) ( η − y ) ( ξ − x ) 2 + r ⋅ ( ζ − z ) + ( ζ − z ) 2 + k 5 ⋅ a r c t a n ( ξ − x ) ( η − y ) ( η − y ) 2 + r ⋅ ( ζ − z ) + ( ζ − z ) 2 + k 6 ⋅ a r c t a n ( ξ − x ) ( η − y ) r ⋅ ( ζ − z ) } | x 0 + d x 2 x 0 − d x 2 | y 0 + d y 2 y 0 − d y 2 | z 0 + d z 2 z 0 − d z 2 (3) where r = [ ( ξ − x ) 2 + ( η − y ) 2 + ( ζ − z ) 2 ] 1 / 2 , μ ₀ is the vacuum permeability; dx, dy, and dz are the dimensions of the grid in the x, y, and z directions, respectively; (x ₀ , y ₀ , z ₀ ) are the grid center coordinates   k 1 = cos ( I 0 ) ⋅ sin ( π 2 − D 0 ) ⋅ sin ( I ) + cos ( I ) ⋅ sin ( π 2 − D ) ⋅ sin ( I 0 ) , k2 = cos ( I 0 ) ⋅ cos ( π 2 − D 0 ) ⋅ sin ( I ) + cos ( I ) ⋅ o s ( π 2 − D ) ⋅ sin ( I 0 ) , k3 = cos ( I 0 ) ⋅ cos ( π 2 − D 0 ) ⋅ cos ( I ) ⋅ sin ( π 2 − D ) + cos ( I ) ⋅ sin ( π 2 − D ) ⋅ cos ( I 0 ) ⋅ cos ( π 2 − D 0 ) , k4 = cos ( I 0 ) ⋅ cos ( π 2 − D 0 ) ⋅ cos ( I ) ⋅ cos ( π 2 − D ) , k5 = cos ( I 0 ) ⋅ sin ( π 2 − D 0 ) ⋅ cos ( I ) ⋅ sins ( π 2 − D ) , and k6 = - sin ( I 0 ) ⋅ sin ( I ) .

If there is no remanence, and assuming that the inclination (I ₀ ) and declination (D ₀ ) of each data grid are known, then the oblique magnetization forward matrix A can be directly calculated. In this study, the inclination and declination were directly incorporated into the total magnetic intensity inversion process, ΔT. Zhdanov (2002) reported the iteration process of the conjugate gradient.

High-precision magnetic anomaly data were acquired in the study area using five G-856 and four G-858 magnetometers. We set up nine magnetic observation stations for diurnal variations in Jiangyou, Shuijingbao, Sedi, Lianghekou, Hongyuan, Jiuzhi, Maqu, Diebu, and Dangchang to adequately and effectively control the measurement area and improve the accuracy of each measurement point. The locations were free from human interference and active magneticity, with the convenience of use, clear markings, and accessible storage. Joint measurements were conducted at these stations from 5–8 am and 5–10 pm (when the magnetic field was relatively stable) to achieve magnetic field normalization. The data recorded during these two periods were used to obtain the basic magnetic field value of the observation stations. The radius for the joint measurements was controlled to within 50 km. The data collected included the total magnetic anomaly △T, with a total root-mean-square error of the magnetic measurements <5nT. The field survey was set up with 13 survey lines spaced 30 km apart. The measurement points were arranged on each line at 2-km intervals, resulting in 1,985 points (see Table 1).

Data from other areas were collated and supplemented using Satellite Magnetic Anomaly Data (https://www.ngdc.noaa.gov/geomag/emag2.html) (Figure 2A) from the Earth Magnetic Anomaly Grid (EMAG2) database published by the Commission for the Geological Map of the World (Maus et al., 2009). This database of crustal magnetic anomalies was accumulated from airborne, offshore, and satellite magnetic surveys over 50 years at a 2’ × 2’ global grid resolution. As shown by the magnetic anomaly map of the Songpan–Aba area, the object under investigation in this study was characterized by stable weak magnetic anomalies, varying from −30 to −50 nT in a steady, smooth magnetic field. The magnetic field in the eastern part of the survey area varied drastically with local anomalies. The measurement area was divided into the following zones according to the strength of the magnetic anomalies:

1) North China plate: This region exhibited varying positive and negative magnetic anomalies. Numerous rock outcrops primarily formed in the Indosinian period, and a some formed during the Hercynian and Yanshanian periods, exist in this region. The distribution of exposed rocks was well correlated with the magnetic anomalies. The magnetic anomaly was in the north-northeast direction and positive, with the highest value >200 nT. The rocks formed during the Indosinian period exhibited various magnetic properties, ranging from weak to strongly magnetic, with their magneticity related to the magnetic minerals content of the rocks. The ferromagnetic mineral content is low in acidic rocks, which exhibit weak to moderate magneticity, whereas intermediate mafic rocks are more magnetic and show high magnetic anomalies.

First, a theoretical model was designed to verify the practicality and stability of the proposed algorithm, then a model consisting of three blocks of different depths and sizes was designed, as shown in (Figure 3A). The model space was 40 × 40 × 21 grid cells, that were each 1 × 1 × 1 km in size, and the magnetization of the model was 1 A/m. Figure 3B shows the synthetic data with the I ₀ = 52°and D ₀ = −1°, where the data space was 40 × 40 grid cells, and each cell was 1 × 1 km. A profile selected in this study is indicated by the pink line in Figure 3B. Regularized inversion was then performed, and the results are shown in Figure 3C. The profile was selected for a detailed comparison of the model and inversion results, as indicated in Figure 3D.

Further, the number of grids in the east–west, north–south, and vertical directions were 97, 81, and 40, respectively, for the 3D inversion of the measurement data obtained from the 5 × 5 km horizontal-dimension grids in the study area. The vertical depth of the forward modeling and inversion was on average, divided into 40 grids covering 40 km (the vertical interval for the grid is 1 km). One hundred iterations were performed in the inversion, in which the data-fitting error was reduced from an initial value of six root-mean-square (RMS) to approximately two RMS. Figure 4 shows the results of the regularized inversion based on the weighted model. We derived the trend and distribution of the magnetic anomaly in the study area from the iso-surface images.

A method for solving the non-uniqueness of geophysical inverse problems is described in the Tikhonov and Arsenin (1977) regularization theory that uses both data misfit and a model-stabilizing functional to construct the parametric objective functional. The objective functional of the regularized inversion can be written as P ( m,d ) = ‖ A ( m ) − d ‖ 2 + α φ ( m ) (4) where p (m,d) is the parametric objective functional, m is a vector consisting of model parameters, A(m) is the forward operator, d is a vector containing the observation data, ||·|| is the L2-norm, φ(m) is the model-stabilizing functional, and α is the regularization parameter.

A stabilizing functional is used to select the appropriate class of models from the model space of all possible models that are a correction set for solving the inverse problem (Zhdanov, 2002). The minimum support (MS) stabilizing functional was proposed to obtain a model that minimizes the volume with anomalous model parameters to improve the resolution of the sharp electrical interface (Last and Kubik, 1983). It focuses on where/how the model parameters are different with respect to the given a priori model. Here, the MS stabilizing functional can be expressed as φ M S ( m , β ) = ‖ m − m a p r ( m − m a p r ) 2 + β 2 ‖ 2 . (5)

Xiang et al. (2017) and Zhang et al. (2020) analyzed the advantages and disadvantages of several commonly used families of stabilizing functionals, such as the smooth, MS or minimum gradient support (MSG). We then applied a new stabilizing functional, the MSG that combines smooth and sharp-boundary constraints with a unified form for the objective functional in the 2-D regularized inversion. The MSG stabilizing functional is defined as: φ M S G ( m , β ) = ‖ ∇ ( m − m a p r ) ( m − m a p r ) 2 + β 2 ‖ 2 . (6)

The spatial gradient is calculated after the MS stabilizer is obtained using Eq. 3, thus, we refer to this new stabilizing functional as the MSG, which uses an MS value instead of m—m _apr, leading to a stable constraint that focuses on sharp boundaries for inversion. This functional acts on the class of models such as the domain with an anomalous parameter distribution that occupies the minimum volume, while avoiding any instability caused by simultaneously focusing parameter tuning. The MSG stabilizing functional is minimally affected by the a priori model. Zhdanov (2002) proved that the MS and MGS satisfy the Tikhonov criterion for a regularization stabilizer. Since the MSG is a spatial gradient of the MS, it also satisfies the regularization criterion.

Four survey lines (Lines 1, 2, 3, and 4) were identified for MT data acquisition in the Songpan–Aba area. Lines one and two were parallel and trending in a SE−NW direction with an ∼45° azimuth, whereas Lines three and four were trending in an SW−NE direction with the same azimuth. The survey lines spanned the entire Songpan–Aba area. Field testing was conducted using Model V-2000 devices in a five-component tensor-impedance format, in 2002. In total, 809 points were collected from the 4 MT profiles in the region. The detailed line layout and workload requirements are shown in Figure 2B and Table 2, respectively. The survey was fully implemented using fixed-station and remote-reference methods throughout the survey area. The effective recording bands were 320–0.001 Hz, the electrode distance was 100 m, and the ground resistance was <1 KΩ. The horizontal magnetic rods were buried below 50 cm and the vertical magnetic rods were inserted approximately three-quarters into the ground. Differential GPS positioning was applied to determine the locations of the survey points, maintaining the precision of the plane coordinates within ±10 m. The measurement network was maintained in a regular form, with point distance and lateral deviations within 400 m. Maintenance was conducted to ensure that at least 85% of the physical points were working during the survey.

First, a theoretical model was designed (Figure 5A) with the resistivity of wall rock at 100 Ω m and the resistivity of a high-conductivity strip at 5 Ω m. The above-mentioned two-dimensional MT-regularized inversion method, with inversion parameters based on the constraints of the MSG model, was used to perform the inversion for this theoretical model. The horizontal grids were set up on a 10-km scale and the vertical grids were divided into 86 layers, which were equally spaced in a log scale up to 300 km. The model response incorporated 5% Gaussian random noise, with 40 frequencies from 320–0.00055 Hz. The initial model was an infinite half-space of 100 Ω m. The inversion results are shown in Figure 5B, which illustrates that the resistivity and distribution of the high-conductivity strip can be accurately delineated in the model, with the fitting error of the inversion data reduced to 1.06.

Further, this two-dimensional MT-regularized inversion method was applied to the inversion of the four profiles in the Songpan area. The inversion parameters were as follows: the horizontal grid-scale was the spacing of the measurement points, and the vertical grids were divided into 86 layers, which were equally spaced in a log scale up to 300 km. The maximum number of iterations in the inversion was 50, and the error level was 5% (i.e., errors with an initial value < 5% were set to 5%, while >5% remained unchanged). The initial model was an infinite half-space of 100 Ω m. The four profiles converged to 1.5, 1.8, 1.3, and 1.2, respectively, using the chi-square test error as the convergence criterion. As shown in Figures 6−9 B, the resistivity structure of the four profiles revealed a vertical distribution of high-resistance-low-resistance-high-resistance. A high-resistance layer ∼20 km thick existed in the shallow ground, composed of sedimentary materials and substances from the upper crust. A high-conductivity layer was located 20–60 km below the study area, while heterogeneous high-resistance and high-conductivity anomalies appeared below 60 km. A connection was found between the high-conductivity layer (20–60 km) and the high-conductivity anomalies below 60 km; therefore, it was concluded that the high-conductivity layer between 20 and 60 km might be formed by a substance from the deep mantle which was brought to the crust through upwelling channels. A comprehensive elaboration on this conclusion, integrating the magnetic anomaly inversion results will be presented in the following section.