The estimated magnetization inclinations for rectangular model using correlation method.
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.
- magnetic anomaly
- inversion and interpretation
- remanent magnetization
- magnetic amplitudes
- mineral resources
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 ; 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 , 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  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 ; Bilim and Ates estimated the magnetization direction by searching for the maximum correlation between pseudo-gravity and gravity anomalies . 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 . Nicolosi et al. computed the magnetization direction of crustal structures using an equivalent source algorithm . Dannemiller and Li estimated the total magnetization direction based on the correlation between the vertical gradient and the total gradient of the RTP field . Gerovska et al. inverted the magnetization direction by correlating RTP and the magnitude magnetic anomalies . Li et al. estimated the magnetization direction of magnetic anomalies through the correlation between normalized source strength and RTP . 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 , 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 . 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 . 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 . 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 .
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 . 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 . The self-demagnetization widely exists in magnetic exploration  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 . 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 . 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 . Liu et al. inverted the 2D magnetization magnitude and direction distributions of a high-susceptibility dike model using borehole magnetization vector inversion . Krahenbuhl and Li implemented the inversion of multiple source bodies and complex structures exhibiting strong self-demagnetization based on the magnetic amplitude data .
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.1. 2D magnitude magnetic anomaly
The magnitude magnetic anomaly (i.e., ) belongs to the magnitude transforms and has some differences compared with the total field anomaly (i.e., ). The 2D magnitude magnetic anomaly is defined as
where and are the horizontal and vertical components of magnetic anomalies, respectively. Given that and anomalies satisfy the linear superposition principle, can be added by anomalies from each magnetic cell numbered
For single magnetic cell, the magnitude anomaly can be written as
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:
where is a constant, is the residual density,is the effective magnetization intensity, and , , , and 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.,
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 = 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.
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., ) belongs to the transformed quantity which is computed from observed total field anomaly (i.e., ) in frequency domain . Initially, we calculate the frequency spectrum of by implementing fast Fourier transform (FFT):
where represents the frequency spectrum of . Then, the frequency spectrums of horizontal and vertical components (i.e., and ) are achieved by, respectively, multiplying a transformed factor on the frequency spectrum:
where and are the frequency spectrums of and components; and are the frequency factors transforming to and components and expressed as:
After obtaining the and 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 and 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 is nonlinearly related to the magnetization intensity vector
The minimum error solution of Eq. (15) is equivalent to solving the symmetric positive definite equation:
We multiply a matrix
where is the buried depth of mesh cells, is a constant related to the magnetic anomalies’ attenuation rate with the increase of distances between cells and observation point, and
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 to with a certain step:
where is the magnetization inclination, is the correlation coefficient between observed and predicted total field anomalies, is the observed total field anomalies, is the predicted total field anomalies when the magnetization inclination is set to be , and , are the covariance of observed and predicted data, 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,
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 to . 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
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).
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.,
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
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.
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.,
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.
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
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 (
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
Which, however, has low sensitivity to the direction of magnetization. Stavrev and Gerovska  and Pilkington and Beiki  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
where is the magnitude magnetic transform;
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.
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
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.
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).