Inversion and Interpretation of Magnetic Anomaly in the Presence of Significant Remanence and Self-Demagnetization Based on Magnetic Amplitude

Abstract

Remanent magnetization and self-demagnetization effects of high-susceptibility body distort the intensity and direction of internal magnetization and hence complicate the inversion and interpretation of magnetic anomaly. The magnitude magnetic anomaly, which is weakly sensitive to the magnetization direction, provides an indirect way to investigate these complex anomalies. We study the sensitivity characteristics of 2D magnitude magnetic anomaly to magnetization direction and source shapes, implement the magnetization intensity inversion, and further estimate the magnetization direction by inverting for the total field data. The magnetic amplitude inversion is tested by the use of synthetic data, which are caused by prism models with strong remanent magnetization and high susceptibility. It is also applied to the field data of an iron-ore deposit in South Australia. The primary advantage of magnitude anomaly inversion is that the magnetization directions are not assumed to parallel the geomagnetic field. The magnetization intensity inversion and magnetization direction estimation make full use of the amplitude and phase information of magnetic anomalies. Magnetic amplitude inversion including other amplitude quantities such as normalized source strength and analytic signal offers an effective approach to investigate and interpret the magnetic anomalies affected by complicated remanence and self-demagnetization.

1. Introduction

Susceptibility is the primary parameter used to represent the magnetic property of rocks and ores. Thus, numerous studies of magnetic data inversion were devoted to recover the susceptibility distributions and to infer the positions and shapes of magnetic sources. For example, Li and Oldenburg presented the techniques of depth-weighted, wavelet-transform compression and joint inversion of surface and borehole magnetic data [1, 2, 3]; Pilkington utilized the preconditioned conjugate gradient algorithm to solve the matrix equation [4]; Portniaguine and Zhdanov presented the image focusing techniques based on the minimum gradient support functions [5, 6]; Fedi presented the depth from extreme point method based on upward continuation theory [7], and so on. Nonetheless, the remanent magnetization also is an important part of rocks’ and ores’ magnetic property. It originates from conditions at their time of formation and widely exists in many real examples. The remanence alters the strength and direction of internal magnetization and exhibits large extents of uncertainty and regionality, which complicates the interpretation of magnetic data. Apart from the remanence, the self-demagnetization effect of high-susceptibility field sources also changes the magnitude and direction of internal magnetization [8, 9, 10, 11]. Inversion of magnetic anomaly in the presence of remanence and self-demagnetization has become a hot topic in recent years.

Some strategies have been proposed to deal with the remanence problem as Clark [12] summarized. The first kind of approach is to estimate the magnetization direction before recovering the physical property distributions using a standard magnetic inversion. The methods of estimating the magnetization direction are many. For example, Fedi et al. proposed the max-min method of reduced-to-the-pole (RTP) to obtain the magnetization direction [13]; Bilim and Ates estimated the magnetization direction by searching for the maximum correlation between pseudo-gravity and gravity anomalies [14]. Phillips used Helbig’s integrals for estimating the vector components of the magnetic dipole moment from the first-order moments of the vector magnetic field components [15]. Nicolosi et al. computed the magnetization direction of crustal structures using an equivalent source algorithm [16]. Dannemiller and Li estimated the total magnetization direction based on the correlation between the vertical gradient and the total gradient of the RTP field [17]. Gerovska et al. inverted the magnetization direction by correlating RTP and the magnitude magnetic anomalies [18]. Li et al. estimated the magnetization direction of magnetic anomalies through the correlation between normalized source strength and RTP [19]. The above magnetization direction estimation approaches are more amenable for simple and isolated anomalies because usually a unique magnetization direction is achieved. Due to the fact that only an averaged magnetization direction can be achieved, the above methods are more applicable for some simple and isolated anomalies.

In addition, an alternative method is presented to directly invert for some kinds of amplitude anomalies which are low sensitive to magnetization directions such as the total gradient data [20], the magnitude magnetic anomaly [21, 22, 23, 24, 25], the normalized source strength [26, 27, 28], and the analytic signal [29, 30]. This approach is more effective when the magnetization direction is highly variable due to the factors such as structural changes. The third method is the magnetization vector inversion. Wang et al. recovered a three-component Cartesian magnetization model and inverted the three components of total magnetization. However, their approach was more applicable in determining the total magnetization of separated, homogeneous bodies [31]. Lelièvre and Oldenburg improved their methods and calculated the three components of magnetization in a Cartesian and spherical framework, which served more complicated scenarios and had widespread applications in magnetic data inversion under the influences of significant remanent magnetization [32]. Similarly, Ellis et al. established the matrix equations between the magnetization components and magnetic anomalies and then optimized the objective function to obtain three components of magnetization vector [33]. These proposed methods were more applied to the inversion in the presence of remanence, but seldom used to invert the high-susceptibility distributions when the self-demagnetization effect is considered. Liu et al. inverted the 2D magnetization vector distributions of a high-susceptibility prism model based on the borehole magnitude magnetic anomalies [8].

The physical principles of self-demagnetization are different with remanence, but they have similar response in that both change the strength and direction of internal magnetization vectors. In cases where susceptibility of magnet is <0.1 SI, the effects of demagnetization are insignificant and can be neglected in forward modeling. However, such effects are important when modeling bodies with high susceptibility [11]. It is demonstrated that the self-demagnetization effect tends to reduce the magnitude and biases the direction of internal magnetization, thereby distorting the amplitudes and shapes of magnetic anomalies [34]. The self-demagnetization widely exists in magnetic exploration [10] and engineering prospecting [9, 35]. A large number of forward modeling studies in relation to self-demagnetization have been carried out [9, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50], but it is still difficult to invert the property distributions considering the implications of self-demagnetization. The earlier approaches dealing with this problem involve correcting the magnetic anomaly by the use of the demagnetization factors for some simple models such as the 3D sphere and 2D elliptic cylinder [10]. Obviously, this method is only suitable for some simple geological conditions. Being similar to the electrical methods, in addition, Lelièvre and Oldenburg directly solved the Maxwell’s equations using a finite volume discretization to recover the 3D high susceptibility distributions [9]. This is an effective way to solve the self-demagnetization problem, but the algorithms for solving the partial differential equations are difficult to implement under the complicated boundary conditions and rugged topography. Krahenbuhl and Li proposed an amplitude inversion method, and the study gave good results under the influence of self-demagnetization at high magnetic susceptibility [35]. Liu et al. inverted the 2D magnetization magnitude and direction distributions of a high-susceptibility dike model using borehole magnetization vector inversion [8]. Krahenbuhl and Li implemented the inversion of multiple source bodies and complex structures exhibiting strong self-demagnetization based on the magnetic amplitude data [50].

Taking the 2D magnetic anomaly as an example, we discuss the characteristics of magnitude magnetic anomaly to magnetization direction and shapes of magnetic sources, and introduce the computations of magnitude magnetic anomaly in frequency domain. We recover the magnetization intensity distributions from the magnitude magnetic anomaly. Then, the total field anomalies are computed to estimate the magnetization direction. We simulate the magnetic field responses of high-susceptibility source under the self-demagnetization effect using the finite element method (FEM) and use the synthetic prism models with significant remanent magnetization and high susceptibility to test the amplitude inversion, respectively. Finally, amplitude inversion is applied to the field data of an iron-ore deposit in Southern Australia.

2. Methodology

2.1. 2D magnitude magnetic anomaly

The magnitude magnetic anomaly (i.e., $T_a$) belongs to the magnitude transforms and has some differences compared with the total field anomaly (i.e., $∆T$). The 2D magnitude magnetic anomaly is defined as

$$T_a={\left(H_{\mathit{ax}}^2+Z_a^2\right)}^{\frac{1}{2}},$$E1

where $H_{\mathit{ax}}$ and $Z_a$ are the horizontal and vertical components of magnetic anomalies, respectively. Given that $H_{\mathit{ax}}$and $Z_a$ anomalies satisfy the linear superposition principle, $T_a$can be added by anomalies from each magnetic cell numbered i:

$$T_a={\left[{\left({\sum }_i{H}_{\mathit{axi}}\right)}^2+{\left({\sum }_i{Z}_{\mathit{ai}}\right)}^2\right]}^{\frac{1}{2}}.$$E2

For single magnetic cell, the magnitude anomaly can be written as

$$T_{\mathit{ai}}={\left(H_{\mathit{axi}}^2+Z_{\mathit{ai}}^2\right)}^{\frac{1}{2}}.$$E3

Based on Eqs. (2) and (3), the following equation can be deduced:

$$T_a\le {\sum }_i{T}_{\mathit{ai}},$$E4

which indicates that the magnitude magnetic anomaly, the first difference with total field, is nonlinear relative to the magnetization intensity. It complicates the forward modeling and inversion. For example, the magnitude magnetic anomalies cannot be computed by adding single mesh cell’s anomalies. Also, their sensitivity matrix is more complex to calculate than that of total field anomalies.

According to the 2D Poisson formula of magnetic field, the horizontal and vertical component anomalies are given by:

$$\left\{\begin{array}{c}{H}_{\mathit{ax}}=\frac{{\mu }_0}{4\mathrm{\pi G\sigma }}\left[M_x{V}_x+M_z{V}_z\right]\\ Z_a=\frac{{\mu }_0}{4\mathrm{\pi G\sigma }}\left[M_x{V}_x+M_z{V}_z\right]\end{array}\right.,$$E5

where V is the gravitational potential. Thus, $T_a$can be written as

$$T_a={\left(H_{\mathit{ax}}^2+Z_a^2\right)}^{\frac{1}{2}}=\frac{{\mu }_0{M}_s}{4\mathrm{\pi G\sigma }}{\left[V_{\mathit{xx}}^2+V_{\mathit{xz}}^2\right]}^{\frac{1}{2}}=\frac{{\mu }_0{M}_s}{4\mathrm{\pi G\sigma }}{\left[V_{\mathit{zz}}^2+V_{\mathit{zx}}^2\right]}^{\frac{1}{2}},$$E6

where $\frac{{\mu }_0{M}_s}{4\mathit{\pi G\sigma }}$ is a constant, $\sigma$ is the residual density,$M_s$ is the effective magnetization intensity, and $V_{\mathit{xx}}$, $V_{\mathit{zz}}$, $V_{\mathit{xz}}$, and $V_{\mathit{zx}}$ are the second-order partial derivatives of gravitational potential, none of which are dependent on the direction of magnetization. Eq. (6) demonstrates that the magnitude magnetic anomaly is not dependent on the magnetization orientation.

Figure 1 shows the examples of total field anomalies and magnitude magnetic anomalies of a rectangular prism model magnetized by different magnetization inclinations (i.e., I = 0°, 30°, 60°, 90°, 135°), of which the magnetization intensity is $M$ = 100 A/m and geomagnetic inclination and declination are I₀ = 45° and D₀ = 0°, respectively. For the different magnetization inclinations, the total field anomalies (colored solid lines) are different, but their corresponding magnitude magnetic anomalies (black solid line) are completely similar. In fact, the observed total field anomalies contain amplitude and phase information. The magnitude anomalies reflect their amplitude information, which is related to the intensity of magnetic bodies’ magnetization or susceptibility. The phase information is related to the magnetization direction. Therefore, the primary advantage of using amplitude anomalies of total field anomalies to invert the magnetization intensity or susceptibility distributions is that it is not essential to input the magnetization inclination or not to assume the magnetization inclination paralleling the Earth’s magnetic field.

Figure 1.

(a) The rectangular prism model (b) corresponding total field anomalies and magnitude magnetic anomalies.

Besides, the magnitude magnetic anomalies show better discrimination to the occurrences of the magnetic bodies. As for tabular models with different inclinations, the shapes of magnitude magnetic anomalies are different. In Figure 2, the magnitude anomalies of vertical tabular bodies (i.e., α = 90°) are symmetrical with similar decrease rate at the two sides of the anomalies. As for dipping tabular bodies with inclinations $\alpha$ = 15°, 30°, 45°, 60° and 75°, the magnitude anomalies at inclined direction decrease slightly, while they are steep at another direction. Therefore, the magnitude magnetic anomalies are conveniently utilized to preliminarily determine the inclined direction of magnetic bodies.

Figure 2.

The magnitude magnetic anomalies of the dipping prisms with different inclinations.

Overall, magnitude magnetic anomaly has two advantages. First, it does not depend on the magnetization direction. The magnitude anomaly has the approximate resolution to recover the magnetic sources as total field anomaly. Second, magnitude anomaly performs more accurately in estimating the occurrence of magnetic bodies. Owing to the influences of inclined magnetization, it is difficult to determine the dipping direction using total field anomaly.

2.2. Computation of magnitude magnetic anomaly

Magnitude magnetic anomaly (i.e., $T_a$) belongs to the transformed quantity which is computed from observed total field anomaly (i.e., $∆T$) in frequency domain [51]. Initially, we calculate the frequency spectrum of $∆T$ by implementing fast Fourier transform (FFT):

$$∆T{\displaystyle \stackrel{\mathit{FFT}}{\to }}{\mathcal{F}}_{∆T},$$E7

where ${\mathcal{F}}_{∆T}$represents the frequency spectrum of $∆T$. Then, the frequency spectrums of horizontal and vertical components (i.e., $H_{\mathit{ax}}$ and $Z_a$) are achieved by, respectively, multiplying a transformed factor on the $∆T$ frequency spectrum:

$${\mathcal{F}}_{H_{\mathit{ax}}}={\omega }_{∆T\to H_{\mathit{ax}}}{\mathcal{F}}_{∆T},$$E8

and

$${\mathcal{F}}_{Z_a}={\omega }_{∆T\to Z_a}{\mathcal{F}}_{∆T},$$E9

where ${\mathcal{F}}_{H_{\mathit{ax}}}$ and ${\mathcal{F}}_{Z_a}$ are the frequency spectrums of $H_{\mathit{ax}}$and $Z_a$ components; ${\omega }_{∆T\to H_{\mathit{ax}}}$ and ${\omega }_{∆T\to Z_a}$ are the frequency factors transforming $∆T$ to $H_{\mathit{ax}}$ and $Z_a$ components and expressed as:

$${\omega }_{∆T\to H_{\mathit{ax}}}\left(\kappa \right)=\frac{i}{\mathit{\text{sin}}{I}_0+\mathit{\text{icos}}{I}_0\mathit{\text{cos}}{A}_0},$$E10

and

$${\omega }_{∆T\to Z_a}\left(\kappa \right)=\frac{1}{\mathit{\text{sin}}{I}_0+\mathit{\text{icos}}{I}_0\mathit{\text{cos}}{A}_0},$$E11

where I₀ and A₀ denote the inclination of the Earth’s magnetic field and the azimuth of profile, and i is the imaginary number. Hence, the horizontal and vertical component anomalies are obtained by carrying out the inverse frequent Fourier transform (IFFT):

$${\mathcal{F}}_{H_{\mathit{ax}}}{\displaystyle \stackrel{\mathit{\text{IFFT}}}{\to }}{H}_{\mathit{ax}},$$E12

and

$${\mathcal{F}}_{Z_a}{\displaystyle \stackrel{\mathit{\text{IFFT}}}{\to }}{Z}_a.$$E13

After obtaining the $H_{\mathit{ax}}$and $Z_a$ components, the magnitude anomalies are computed by the use of Eq. (1).

Eqs. (7)–(13) summarize the calculation processes of magnitude magnetic anomaly. The critical processes are based on the computations of $H_{\mathit{ax}}$ and $Z_a$ components in frequency domain. Regardless of geomagnetic inclination and profile’s azimuth, it does not need inputting another parameter during the whole calculation processes.

2.3. Magnetization intensity inversion

Eq. (4) indicates that the magnitude magnetic anomaly vector ${\mathbf{T}}_a$ is nonlinearly related to the magnetization intensity vector m. Their relation is given by

$${\mathbf{T}}_a=T_a\text{(}\mathbf{m}\text{)}$$E14

where $T_a$(m) is a nonlinear function. Using first-order Taylor expansions on Eq. (14), we obtain the following matrix equation:

$${\mathrm{\Delta }\mathbf{T}}_a={\mathbf{J}}_{T_a}\mathrm{\Delta }\mathbf{m}$$E15

where ${\mathrm{\Delta }\mathbf{T}}_a$ and$\mathrm{\Delta }$m are corrections of ${\mathbf{T}}_a$and m, respectively, and ${\mathbf{J}}_{T_a}$ is named sensitivity matrix whose element $J_{T_a\left(i,j\right)}$ is the partial derivative of $T_a$at the ith observation point to the jth model parameter m_j [23]. Thus,

$$J_{T_a\left(i,j\right)}=\frac{\partial T_{\mathrm{ai}}}{\partial m_j}=\frac{\partial }{\partial m_j}{\left(H_{\mathit{axi}}^2+Z_{\mathit{ai}}^2\right)}^{\frac{1}{2}}=\frac{1}{T_{\mathit{ai}}}\left(H_{\mathit{axi}}{G}_{H_{\mathit{ax}\left(i,j\right)}}+Z_{\mathit{ai}}{G}_{Z_{\mathit{a}\left(i,j\right)}}\right),$$E16

where i (i = 1,2,…,m) denotes the ith observation point, j (j = 1,2,…,n) denotes the jth mesh cell, m and n are the total numbers of observation points and mesh cells, respectively, $$and $G_{H_{\mathit{ax}\left(i,j\right)}}$and $G_{Z_{a\left(i,j\right)}}$ are the elements of ${\mathit{G}}_{H_{\mathit{ax}}}$ and ${\mathit{G}}_{Z_a}$. Here, ${\mathit{G}}_{H_{\mathit{ax}}}$ and ${\mathit{G}}_{Z_a}$ are the constant sensitivity matrices of $H_{\mathit{ax}}$ and $Z_a$ anomalies. This equation demonstrates that because of the nonlinearity relation between magnitude anomalies and magnetization intensity, it is more complicated to compute the sensitive matrix of $T_a$than that of $H_{\mathit{ax}}$, $Z_a,$and $∆T$ anomalies. The sensitivity matrix of $T_a$is related to the sensitive matrices of $H_{\mathit{ax}}$, $Z_a$ and the calculated $H_{\mathit{ax}}$, $Z_a$, and $T_a$anomalies.

The minimum error solution of Eq. (15) is equivalent to solving the symmetric positive definite equation:

$${\mathbf{J}}_{T_a}^T{\mathbf{J}}_{T_a}\mathrm{\Delta }\mathbf{m}={\mathbf{J}}_{T_a}^T{\mathrm{\Delta }\mathbf{T}}_a.$$E17

We multiply a matrix P at both sides of Eq. (17); hence,

$$\mathbf{P}\left({\mathbf{J}}_{T_a}^T{\mathbf{J}}_{T_a}\mathrm{\Delta }\mathbf{m}\right)=\mathbf{P}\left({\mathbf{J}}_{T_a}^T{\mathrm{\Delta }\mathbf{T}}_a\right),$$E18

where P usually is a diagonal matrix named preconditioner, which is used to reduce the condition number of Eq. (17) and promote the convergence rate. When Pilkington [4] and Liu et al. [8] used preconditioned conjugate gradient method to invert the magnetic anomaly, the preconditioner is given by.

$$\mathbf{P}=z^{\beta }\mathbf{I},$$E19

where $z$ is the buried depth of mesh cells, $\beta$ is a constant related to the magnetic anomalies’ attenuation rate with the increase of distances between cells and observation point, and I is the unit matrix.

2.4. Estimation of effective magnetization direction

After obtaining the magnetization intensity distributions, we can regard them as known information and then calculate the total field anomalies using different magnetization directions. Thus, if the magnetization direction is given appropriately, the computed total field anomalies should fit the observed total field anomalies and their correlation coefficients get to maximum. Therefore, we compute the correlation coefficients between the observed and predicted total field anomalies of which the magnetization inclinations rotate a cycle from $0$ to ${360}^{°}$ with a certain step:

$$R\left(I\right)=\frac{C\left[∆T^{\mathit{obs}},∆T^{\mathit{pre}}\right]}{\sqrt{C\left[∆T^{\mathit{obs}},∆T^{\mathit{obs}}\right]C\left[∆T^{\mathit{pre}},∆T^{\mathit{pre}}\right]}}\left(0^{°}\le I\le {360}^{°}\right),$$E20

where $I$ is the magnetization inclination, $R\left(I\right)$ is the correlation coefficient between observed and predicted total field anomalies, $∆T^{\mathit{obs}}$ is the observed total field anomalies, $∆T^{\mathit{pre}}$is the predicted total field anomalies when the magnetization inclination is set to be $I$, and $C\left[∆T^{\mathit{obs}},∆T^{\mathit{obs}}\right]$, $C\left[∆T^{\mathit{pre}},∆T^{\mathit{pre}}\right]$ are the covariance of observed and predicted data, $C\left[∆T^{\mathit{obs}},∆T^{\mathit{pre}}\right]$ is the cross-covariance between observed and predicted data. Therefore, the most appropriate magnetization inclination is that when the correlation coefficients of Eq. (20) get to the maximal values, it can be expressed as,

$$I_{\text{best}}=\text{arg}\mathit{max}\left[R\left(I\right)\right].$$E21

Eqs. (20) and (21) are the principle formulas used to estimate the magnetization direction. Based on the recovered magnetization distributions, we calculate the predicted magnetic anomalies magnetized by different magnetization inclinations varied from $0$ to ${360}^{°}$. The magnetization inclination with the largest correlation coefficients between the observed and predicted anomalies is defined as the most appropriate magnetization direction. In essence, the method makes use of the phase information of total field anomalies to determine the magnetization direction.

3. Synthetic examples: magnetic amplitude inversion with significant remanence

3.1. Rectangular prism with different magnetization inclinations

We firstly test the method by the use of the 2D rectangular prism in Figure 1a, of which the top buried depth is 150 m and the length and width are 150 and 200 m, respectively. The Earth’s magnetic field intensity is T₀ = 50,000 nT, with inclination I₀ = 45° and declination D₀ = 0°. The rectangular prism is magnetized by a constant magnetization with intensity M = 100 A/m and declination D = 0°, but the magnetization inclinations are assumed to be I = 0°, 30°, 45°, 60°, 90°, and 135° under influence of remanent magnetization. The magnetic observation point spacing is 20 m, and there are 51 points in total. As shown in Figure 1b, the total field anomalies are changed for different magnetization inclinations. However, for different magnetization inclinations varied from 0 to 135°, their amplitude anomalies always are the black solid line of the Figure 1b. The 2D amplitude data are strictly invariant with the magnetization direction. Therefore, using the amplitude anomalies to recover the magnetization, intensity distributions reduces the errors resulting from the incorrect magnetization direction. When inverting for the total field data, it is necessary to input the correct magnetization direction.

We invert for the amplitude data of the rectangular prism model in Figure 1b. The subsurface is divided into 800 (20 rows × 40 columns) mesh cells with size of 25 × 25 m. The preconditioned conjugate gradient algorithm converges stably after hundreds of iterations, and the predicted amplitude data accurately fit the observed data. The recovered magnetization distributions including the position and shape of magnetic sources yield a good approximation with the true model (Figure 3a). The magnitude anomalies show similar resolution to the recovery of physical property distributions compared with the total field anomalies (Figure 3b).

Figure 3.

The magnetization intensity inversion results of the rectangular model using (a) magnitude magnetic anomalies and (b) total field anomalies.

With the known magnetization intensity distribution of Figure 3a, subsequently, the total field anomalies are computed of which the magnetization inclinations are varied a cycle from 0 to 360° by a step of 0.5°. Then, we calculate the correlations between the observed and predicted total field anomalies for each magnetization inclination. As shown in Figure 4, the six colored solid lines, respectively, represent the correlation curves of the six synthetic magnetization inclinations (i.e., I = 0°, 30°, 45°, 60°, 90°, and 135°), of which the horizontal and vertical axes denote the magnetization inclination and the value of correlation. The correlation curve fluctuates between −1 and 1 shaped as the sine or cosine function. Each curve has one maximum point, and the position of the maximum point represents the optimal magnetization direction. For example, for the red solid curve when the magnetization direction is horizontal (i.e., I = 0°), the maximal peak locates at I = 0.5° with maximal correlation R = 0.9986, and therefore, the estimated optimal magnetization is I_best = 0.5°. The results of magnetization direction determination for the six synthetic magnetization inclinations are shown in Table 1. The synthetic magnetization inclinations are 0°, 30°, 45°, 60°, 90°, and 135°, and their corresponding estimated magnetization inclinations are 0.5°, 30.0°, 45.0°, 60.0°, 90.0°, and 135.0°. The maximal correlations resulted to 0.9986. The estimated magnetization inclinations are in agreement with the true values (error < 0.5°).

Figure 4.

The correlation curves of the rectangular model with different magnetization inclinations. The positions of maximal values reflect the magnetization directions.

Synthetic (degree)	0	30	45	60	90	135
Estimated (degree)	0.5	30.0	45.0	60.0	90.0	135.0
Correlation	0.9987	0.9987	0.9988	0.9987	0.9987	0.9986

Table 1.

The estimated magnetization inclinations for rectangular model using correlation method.

3.2. Complicated prisms with the same magnetization inclination

We design four 2D prism models, of which the cross sections are the dipping tabular, syncline tabular, cut tabular, and reproduction tabular (Figure 5). The Earth’s magnetic field intensity is T₀ = 50,000 nT, with inclination I₀ = 45° and declination D₀ = 0°. The models are magnetized by a constant magnetization magnitude M = 100 A/m with inclination I = 60° and declination D = 0°. The magnetic observation point spacing is 20 m, and there are 51 points in total. The observed total field anomalies and the transformed magnitude anomalies are shown in Figure 6. Being similar to amplitude inversion of the Figure 3, the subsurface is divided into 800 (20 rows × 40 columns) rectangular cells. The preconditioned conjugate gradient converges after 100 iterations on average, and simultaneously, the predicted magnitude anomalies accurately fit the observed magnitude data (Figure 7). The recovered magnetization intensity distributions are shown in Figure 5. The inversion results are basically in accordance with the synthetic models. Owing to the lower resolution due to deep-buried magnetic sources, the deeper sources of combinational models including the cut prisms (Figure 5c) and reproduction prisms (Figure 5d) are not clearly distinguished.

Figure 5.

The magnetization intensity inversion results using magnitude magnetic anomalies for the synthetic prism model: (a) dipping prism, (b) syncline prism, (c) cut prism, and (d) reproduction prism.

Figure 6.

The observed and predicted total field anomalies and magnitude magnetic anomalies of the synthetic prism models: (a) dipping prism, (b) syncline prism, (c) cut prism, and (d) reproduction prism.

Figure 7.

The correlation coefficient curves and the computed magnetization directions of the synthetic prism model: (a) dipping prism, (b) syncline prism, (c) cut prism, and (d) reproduction prism.

After recovering the magnetization intensity distributions (Figure 5), we estimated the magnetization directions based on the known magnetization intensity distributions. The correlation curves and the estimated magnetization inclination are shown in Figure 7 and Table 2. Except the syncline prism model with error of 12°, other determined magnetization inclinations are close to 60°. The amplitude data do not clearly distinguish the closed magnetic bodies leading to the predicted total field data not fitting the observed total field data accurately for combinational syncline and cut prisms. Besides, compared with the traditional correlation methods for estimating the magnetization direction [14, 15, 16, 17, 18], this method considers magnetic sources’ shapes and positions, which improves the precision of magnetization direction determination.

Model	a	B	c	D
Synthetic (degree)	60	60	60	60
Estimated (degree)	63.0	59.5	72.0	62.0
Correlation	0.9978	0.8717	0.9446	0.9994

Table 2.

The magnetization inclination determination results for four prism models: (a) dipping prism, (b) syncline prism, (c) cut prism, and (d) reproduction prism.

Taking the dipping prism as an example, we invert for the magnetization distributions from total field anomalies (true magnetization inclination is 45°) by giving different magnetization inclinations (i.e., I = 0°, 45°, 90°, 135°), assuming that we do not know the real magnetization direction of prism model. Only when the magnetization direction is given to be correct (true magnetization inclination is I = 45°), the inversion results coincide with the true model (Figure 8b). Otherwise, the inversion results have large errors compared with the true model (Figure 8a, c, and d). And even some of the iterations cannot be converged. Therefore, it is essential to give the correct magnetization inclination when using total field anomalies to recover the magnetization distributions. Otherwise, the inversion results may be distorted by the inaccurate magnetization direction.

Figure 8.

The magnetization intensity inversion results for the total field anomalies of the dipping prism model using different magnetization inclinations: (a) I = 0°, (b) I = 45°, (c) I = 90°, and (d) I = 135°.

4. Field examples: Brennand iron-ore deposit, Eyre peninsula, South Australia

The Brennand iron-ore deposit lies in the Eyre Peninsula, South Australia, located at longitude: 135° 52′ 00″ E and latitude: 34° 24′ 00″ S. The banded iron formation (BIF) ore bodies have high magnetic susceptibility and produce strong magnetic anomalies. The aeromagnetic anomalies strike northeast-southwest direction with amplitudes from −1000 to 5000 nT (Figure 9a). We use the data of Line 6 that traverses the center of mining area to test the method. The point spacing of this profile is 20 m, and there are 107 observation points in total.

Figure 9.

The magnetization intensity and direction inversion results for the Line 6 of the Brennand iron-ore deposit, Eyre Peninsula, South Australia: (a) the contour map of total field anomalies, (b) the recovered magnetization intensity distributions, and (c) the estimated magnetization direction.

First, we divided the subsurface into 1060 (20 rows × 53 columns) square cells with size of 40 m. The recovered magnetization intensity distributions indicate that the magnetic bodies are inclined to northwest about 60° and extend downward around 300 m (Figure 9b). Besides, the correlation curve demonstrates that the magnetization inclination is 264.0° (Figure 9c). The declination and inclination of the Earth’s magnetic field in the mining area are NE6.8° and −67.1°. In the profile of Line 6, the effective magnetization inclination increases to −72.4°. Therefore, the magnetization direction deflects 23.6° (i.e., 360–264 − 72.4°) from Earth’s magnetic field.

5. Synthetic example: magnetic amplitude inversion of high-susceptibility body

We design a 2D dike model with high susceptibility = 10.0 SI, of which the four vertices’ coordinates are A (350–300 m), B (450–300 m), C (550–100 m), and D (450–100 m) (Figure 10). The geomagnetic field intensity is T₀ = 50,000 nT with an inclination of 45°. First, the finite element method (FEM) is used to simulate the surface total field anomalies (i.e., $∆T$), the magnitude magnetic anomaly (i.e., $T_a$), the true internal magnetization intensity, and the inclination distributions. Theoretical simulation reveals that the internal magnetization distribution is inhomogeneous with average strength (50–60 A/m) (Figure 10a), which is far less than 397 A/m (i.e., $M=\kappa H_0=\kappa T_0/{\mu }_0$) when demagnetization is not considered. The self-demagnetization effect decreases the internal magnetization intensity. Moreover, the internal magnetization inclination ranges from 50 to 65° with average of 60° (Figure 10b). The self-demagnetization effect deflects the magnetization direction about 15° (i.e., 60–45°) to the long-axis direction of the dike model. The computed total field anomaly and magnitude anomaly are shown with red and blue solid curves in Figure 11. The point spacing is 3 m.

Figure 10.

The magnetization inversion results of the high-susceptibility (κ = 10.0 SI) dike model: (a) the real internal magnetization intensity, (b) direction distributions when the self-demagnetization effect is not negligible, (c) the inverted magnetization intensity distributions using magnetic amplitude inversion, and (d) the estimated magnetization inclination using correlation coefficient method.

Figure 11.

The observed and predicted total field anomalies and magnitude magnetic anomalies.

We inverted for the magnitude magnetic anomaly (i.e., blue curve in Figure 11), and the cross section was divided into 40 × 20 square cells with size of 25 m. The preconditioned conjugate gradient method was converged after 200 times of iterations, and the predicted magnitude anomaly fit the observed anomaly accurately (Figure 11). The inverted magnetization intensity distributions are shown in Figure 10c. The range of the inverted magnetization intensity is (40–65 A/m), of which the amplitudes and shapes are in accordance with the real distribution of magnetization in the presence of demagnetization (Figure 10a). The correlation coefficients (with step of 0.2° varied from 0 to 90°) between the observed and predicted data reach the maximum at A (I = 52.0°, R = 0. 997,696). Hence, for this high-susceptibility prism, the estimated internal magnetization inclination is 52° (Figure 10d). This result is in good agreement with the theoretical distributions of magnetization inclination (Figure 10b). Therefore, from the magnetization vector inversion results, it can be concluded that the demagnetization effect biases the magnetization direction to the dike’s long-axis direction 7° (i.e., 52–45°). They reduce the magnetization intensity and magnetic anomalies by about seven times (i.e., 397/60). The relative error of magnetic anomalies caused by demagnetization effects reaches 85% (i.e., (397–60)/397). Theoretically, for high-susceptibility magnetic bodies of which the susceptibility is 10 SI, the magnetic anomalies’ error caused by demagnetization effects reaches 83% [8].

6. Weak sensitivity of 3D magnitude magnetic anomaly to magnetization direction

The 2D magnitude magnetic anomaly is totally independent of the magnetization direction and has high centricity with the magnetic source’s position, which provides an idea of magnitude magnetic transform to investigate the inversion and interpretation of magnetic anomaly [21, 52]. For 3D cases, magnitude magnetic anomaly is written as

$$T_a={\left(H_{\mathit{ax}}^2+H_{\mathit{ay}}^2+Z_a^2\right)}^{\frac{1}{2}},$$E22

Which, however, has low sensitivity to the direction of magnetization. Stavrev and Gerovska [52] and Pilkington and Beiki [26] used a variable to evaluate the sensitivity of magnitude magnetic transform to magnetization direction by comparing with the field of vertical magnetization direction, expressed as

$$S\left(I,D\right)=\iint \left|F\left(I,D\right)-F\left(I_b,D_b\right)\right|\mathit{\text{dxdy}}/\iint \left|F\left(I_b,D_b\right)\right|\mathit{\text{dxdy}},$$E23

where $F$ is the magnitude magnetic transform; I and D are inclination and declination of magnetization direction, respectively; I_b = 90° and D_b = 90°.

Figure 12 shows the sensitivities of magnitude magnetic anomaly and total field anomaly to the magnetization inclination and declination. It is revealed that both the magnitude magnetic anomaly and total field anomaly are mainly sensitive to magnetization inclination. And the magnitude magnetic anomaly shows far weaker sensitivity than total field. The sensitivity value in Eq. (23) of total field is up to 1.3 when the magnetization direction is horizontal, while it is only 0.3 for magnitude magnetic anomaly. When the magnetization direction is horizontal, for magnitude magnetic data and total field data, the sensitivity reaches a maximal value.

Figure 12.

Sensitivities of (a, c) magnitude magnetic anomaly and (b, d) total field anomaly to the total magnetization direction. Plot (a) shows the sensitivity of magnitude magnetic anomaly to magnetization inclination and declination. Plot (b) shows the sensitivity of total field anomaly to magnetization inclination and declination. Plot (c) shows the sensitivity of magnitude magnetic anomaly to magnetization inclination (i.e., cross section at declination = 50° in plot (a)). Plot (d) shows the sensitivity of total field anomaly to magnetization inclination (i.e., cross section at declination = 50° in plot (b)).

The weak sensitivity feature of 3D magnitude magnetic anomaly to magnetization direction impacts the magnetization intensity inversion results. Figure 13a and b shows the total field anomaly and magnitude magnetic anomaly of the 3D synthetic model when the total magnetization direction is horizontal (magnetization intensity M = 1 A/m, inclination I = 0°, declination = 30°; geomagnetic inclination I₀ = 45°, declination D₀ = 0°). The total fields include positive and negative anomalies showing a relatively complicated feature. The magnitude magnetic anomaly only has positive values, and its extreme points show some offsets from the horizontal centers of magnetic source (Figure 13b). When inverting the magnitude magnetic anomaly, as shown in Figure 13(c–f), the recovered magnetization distributions also have some offsets to north compared with the true models. The deep shapes of recovered magnetization intensity distributions in particularly for sources A and C have big differences with the true models, which also would increase the error of magnetization direction estimation.

Figure 13.

Magnetization inversion results of magnitude magnetic anomaly of a 3D synthetic model with horizontal total magnetization direction: (a) total field anomaly, (b) magnitude magnetic anomaly, (c) horizontal cross section at depth = −175 m, and vertical cross sections at northing = (d) 725 m, (e) 525 m, and (f) 275 m.

7. Conclusions

Remanent magnetization is prevalent in many mining areas, but their directions usually are unknown because of the difficulty in collecting the oriented samples. For high-susceptibility source, the influence of self-demagnetization effect also cannot be ignored. The magnitude anomalies are frequently transformed from the observed total field anomalies. The primary advantage of magnitude anomalies inversion is that the magnetization directions are not assumed to parallel the geomagnetic field, or it is not necessary to input the magnetization directions as the total field data inversion. The magnitude anomalies have similar resolution to the recovery of physical property distributions and perform higher sensitivity to the occurrences of magnetic bodies compared with the total field anomaly. Based on the known magnetization intensity distributions, the total field anomalies are computed by the use of different magnetization directions. Thus, the magnetization direction of maximal correlation with the observed total field anomaly is deemed as the most appropriate magnetization direction. This strategy considers the influences of magnetic sources’ shapes and obtains an average magnetization direction compared with another correlation approaches. The inverted magnetization intensity and direction help to study the influences of remanence and self-demagnetization. The method makes full use of the amplitude and phase information of total field anomaly to determine the magnetization intensity and direction, respectively. The amplitude anomaly is more related to the intensity of magnetization vector, while the phase is mainly dependent on the direction of magnetization vector. The 2D magnitude magnetic anomaly provides an idea of magnitude magnetic transform to implement the amplitude data inversion, but in 3D case, magnitude magnetic anomaly is weakly sensitive to magnetization direction, which brings some errors for the inversion of magnetization intensity. Other magnetic anomaly quantities such as analytic signal and normalized source strength have weak sensitivity to magnetization direction, so they also can be inverted to study the remanent magnetization and self-demagnetization. The amplitude data inversion provides an effective approach to investigate the complex remanence and self-demagnetization.

Acknowledgments

This work was financially supported by the National Natural Sciences Foundation of China (41604087, 41630317), the Natural Sciences Foundation of Hubei Province (2016CFB122), the China Postdoctoral Science Foundation (2016M590132), and the Fundamental Research Funds for the Central Universities, China University of Geosciences (Wuhan) (CUGL170407, CUG160609).