Inversion of high-amplitude magnetic total field anomaly: an application to the Mengku iron-ore deposit, northwest China

In magnetic prospecting, the total field anomaly formula that represents the projection of the magnetic anomaly vector on the geomagnetic field is widely used because it simplifies the calculation of forward modelling and inversion of magnetic data. However, the projection anomaly yields errors relative to the true observed magnetic anomaly, especially for high-amplitude magnetic anomalies such as in iron orebody and unexploded ordnance prospecting. In this study, we analyse the difference between the projection anomaly and observed modulus difference anomaly with physical parameters, and propose to directly invert for the modulus difference anomaly by constructing a nonlinear matrix equation between the model corrections and data corrections. The inversion is then implemented using a preconditioned conjugate gradient algorithm. Synthetic and field magnetic data were used to test the inversion method. Comparison of the two types of total field anomalies shows that the error of the projection anomaly increased with increasing total-field magnetic anomaly. When the total-field magnetic anomaly was < 5,000 nT, the difference between the projection anomaly and modulus difference anomaly results can be ignored. For high-amplitude magnetic anomalies, the modulus difference anomaly inversion produced more accurate representations of both the shape and location of the magnetic sources.

www.nature.com/scientificreports/ relative-error (the difference/the magnitude magnetic anomaly) reach 10% with spherical model. They suggest projection anomaly is not credible especially as high-precision of quantitative processing and interpretation are required and propose a method using the error yielded by the projection anomaly. Coleman and Li 29 studied the difference between the errors of the total-field anomaly and magnetic amplitude data and found that the errors in the three orthogonal components converted from the total-field anomaly had similar standard deviations. Based on the total magnetic anomaly is no longer approximately equal to projection anomaly in the highly magnetic environments, Sun et al. 30 used the information from the borehole data and structural orientation as constraints and get an acceptable result. In our study, we define the error between projection anomaly and observed modulus difference anomaly firstly and carefully discuss the error caused by amplitude, inclination and declination individually and together between the projection anomaly results (data processing, inversion and interpretation) with those of the modulus difference anomaly for strong magnetic anomalies. Second, we propose directly inverting for the modulus difference anomaly in strong magnetic bodies under Cartesian coordinate. Third, the operator, given by modulus difference anomaly, is nonlinear which is different with projection anomaly. Finally, the synthetic and field data are used to test the inversion method with two types of total field data.

Total field anomalies
Two types of total field anomalies: projection anomaly and modulus difference anomaly. Generally, the obtained anomaly T at any point is the vector sum produced by all the magnetizations, comprising two contributions: the normal geomagnetic field T 0 and the magnetic anomaly vector T a (Fig. 1). The modulus difference anomaly (ΔT true ) is the modulus difference between T and T 0 : and Experimental data [31][32][33] confirmed that when the magnetic amplitude data T a is far less than T 0 , the magnetic anomaly can be simplified for forward modelling. In this way, the modulus difference anomaly can be regarded as the projection anomaly (ΔT pro ). To compare the difference between the modulus difference anomaly and projection anomaly, we define the error as: where θ is the angle between the normal geomagnetic field and the magnetic anomaly vector. The Eq. (3) demonstrates that the error arises mainly from the value of the magnitude magnetic anomaly and the angle θ.

Comparison of the projection and modulus difference anomalies
We first test the method through the cube model, which has a length of 100 m and a central burial depth of 100 m. The horizontal coordinate of the centre of the model is (500, 500). The geomagnetic field T 0 is 50,000 nT, the total inclination I is 45° and the total declination D is 0°. The magnetic anomalies can be changed by changing the magnetization intensity (m) of the model.
The projection anomalies, modulus difference anomalies and errors for m = 1, 50 and 500 A/m are shown in Fig. 2. When m = 1 A/m (Fig. 2a) or m = 50 A/m (Fig. 2b), the magnetic anomalies are less than 5,000 nT and the differences between ΔT pro and ΔT true are small. Additionally, the magnetic anomalies calculated by these two methods are very close in shape and value. When ΔT true = 1,794 nT and ΔT pro = 1,783 nT, the maximum www.nature.com/scientificreports/ error is 45 nT, which is less than 2% of the total-field magnetic anomaly (Fig. 2b). In this case, the error that arises from using ΔT pro to replace ΔT true is small and may be neglected. However, when m = 500 A/m (Fig. 2c), ΔT true = 18,920 nT and ΔT pro = 17,830 nT, which results in a big difference between them. The maximum error at this time is 3,945 nT, accounting for 21% of the total-field magnetic anomaly's peak value. In Hinze's research 24 , the error can reach 1,667 nT as T 0 = 30,000 nT and T a = 10,000 nT, which is close to our numerical simulation result. Such a big error (3,945 nT) has a deep negative impact on the interpretation and even absolutely mistake. The error between ΔT pro and ΔT true should be considered carefully when the data indicate high magnetic amplitudes 30 . To understand the difference between ΔT pro and ΔT true , we converted the magnetic anomalies to reduction to the pole (RTP) values for m = 50 A/m (Fig. 3a, b) and m = 500 A/m (Fig. 3c, d). The RTP result of real data (Fig. 3d) cannot eliminate the asymmetry of magnetic anomalies caused by magnetization field, which is unacceptable for interpretation. Both the amplitude and shape of the RTP results in the strong magnetic environments were significantly different (Fig. 3c, d). The error is so obvious that projection anomaly is no longer suitable for interpretation of high-precision measurements. This indicates that ΔT pro cannot be regarded as equal to ΔT true in strong magnetic bodies because the error between them will be further enlarged during data processing and interpretation. Figure 4 shows the error (ΔT true − ΔT pro ) for observed modulus difference anomaly with different amplitudes. When the total-field magnetic anomaly is less than 1,000 nT, the maximum error is less than 10 nT and can be ignored. When the total-field magnetic anomaly reaches 2,000 nT, the maximum error is 54 nT (see Error in www.nature.com/scientificreports/ Table 1). As the magnetic anomaly increases to 5,000 nT and 10,000 nT, the maximum error is close to 335 nT and 1,301 nT, respectively (see Error in Table 1). It is far beyond the precise magnetic allowance. What's more, the error increases exponentially as the total-field magnetic anomaly rises and the order begins to approach the order of the total-field magnetic anomaly. For magnetic bodies of high susceptibility, the real data (observed modulus difference anomaly) have a high amplitude in ground or borehole magnetic surveys 34,35 , especially for shallow  www.nature.com/scientificreports/ ores bodies, which is more than 10,000nT. When we use projection anomaly to fit the real data approximately, it will do harm to exploration resource. For instance, it may lead to spending more time and money on drilling to verification, which causes the waste of resource. According to Eq. (3), the angle θ, between the magnetic anomaly vector T a and the geomagnetic field T 0 , is also a factor affecting the error. The angle θ depends mainly on the magnetic inclination (I) and declination (D). Therefore, we looked at the error (ΔT true − ΔT pro ) caused by the magnetic direction for magnetic inclination values from − 90° to 90° and a magnetic declination range of 0°-360°. Figure 5a shows the error distribution calculated for ΔT true and ΔT pro without remanence. The amplitude of the total field magnetic anomaly is 28,500 nT and the maximum error is 3,992 nT, which is 14% of the total field magnetic anomaly. The influence of the magnetic declination is negligible and the error distribution has good symmetry. When the error is calculated for data with the magnetic remanence ( Fig. 5b; I 0 = 45°, D 0 = 0°), the amplitude of the magnetic anomaly is 24,721 nT and the maximum error is 6,935 nT, which is 28% of the total field magnetic anomaly. The error under the influence of remanence is almost one time bigger than those without remanence. There is no doubt that the difference between ΔT true and ΔT pro is enlarged under the influence of remanence. In practice, most of magnetic ores are affected by remanence. So it is high time that we should take the impact of approximation by projection anomaly in data process into account.
Moreover, the error decreases when I changes from − 90° to 0° and then increases when I changes from 0° to 90°. The influence of D on the error is weak. It depends mainly on I, so we discuss the error caused only by the effect of I. In the cube model, D = 0° and I is set as 0°, 30°, 45° and 90°. We model the effect of I without and with remanence (Figs. 6, 7), respectively. Figure 6a-d shows ΔT true − ΔT pro without remanence for I = 0°, 30°, 45° and 90°, respectively. The error distribution shows good symmetry with respect to the centre of the model, especially when the magnetization is horizontal or vertical. When the remanence is taken into account (I 0 = 45°, D 0 = 0°), the error pattern changes (Fig. 7) in both amplitude and shape.
As discussed above, the magnetic anomalies and directions affect the errors independently. However, the effect on ΔT true − ΔT pro , which combined the amplitude and the direction of the magnetic anomaly, is unclear. Therefore, we use the trigonometric relationship between T, T 0 , T a and θ to clarify this point. Equation (3) can be transformed into:  www.nature.com/scientificreports/ The E values for θ between 0° and 180° are plotted for various T a values and T 0 = 50,000 nT in Fig. 8. The E curves peak at 90° < θ < 120° and E increases as the magnetic anomaly increases. As Eq. (4) shows, the maximum error is reached when ΔT true = 0 nT (see E in Table 1); in this case, the corresponding direction (θ*) meets the following relationship: The geometry of the three parameters is shown in Fig. 9a. The maximum error is plotted in Fig. 9b for θ* ranging from 90° to 120°.
Through the analysis above, we can find there is a huge error between ΔT true and ΔT pro in strong magnetic anomaly. Furthermore, it is more obvious that data-processing and interpretation, such as RTP, are affected by error. It takes disadvantage of our practical exploration. In addition, it is indicated that the error is enlarged and even doubled in the case of remanence. Indeed, most of the ore deposits exist the natural remanence. We should www.nature.com/scientificreports/ take the error seriously and use modulus difference anomaly rather than projection approximation for forward modelling to meet the needs of high-precision and quantitative process and interpretation in mineral exploration.

Inversion of the projection and modulus difference anomalies
The direction of the normal geomagnetic filed is almost constant within a certain range, then it means the geomagnetic inclination I 0 and declination D 0 are certain. The magnitude magnetic anomaly T a can be calculated from its three components (H ax , H ay and Z a ) in x-, y-and z-directions. At the same time, we assume that the components of the normal geomagnetic filed T 0 in x-, y-and z-directions are H ox , H oy and Z o . The components of T 0 are fixed due to the geomagnetic direction is certain. The relationship between the components of magnitude magnetic anomaly T a and magnetization intensity is linear.   www.nature.com/scientificreports/ H i represents one of the three components (H ax , H ay and Z a ) of T a at the ith observation point. G Hi is the sensitivity matrix of H i , whose element is G H(i,j) . m is a matrix including the magnetization intensity of related model. N is the number of models. We assume that each mesh cell has a homogeneous magnetization intensity. ΔT pro can be treated as a harmonic function, while ΔT true is not an additive and harmonic function. Projection anomaly and modulus difference anomaly caused by the grid element of the jth unit magnetic intensity at the ith observation point can be shown as: To solve the inverse problem, we minimize the following objective function: and Subject to 0 < m < m Max .where d represents the geophysical observation data-vector, G, the m × n-dimensional sensitivity matrix, is a function operator and m represents the vector of magnetization intensity of the corresponding model.
The conjugate gradient method is one of the most effective methods for solving inversion problem and has been widely applied, such as geoelectrical [36][37][38] , gravity [39][40][41][42] and magnetic 43,44 data. The conjugate gradient method searches for the optimal solution after finite iteration convergence along the direction of the conjugate gradient. All the operations of the conjugate gradient method are vector operations. Variables can be used repeatedly during the cyclic iteration process, which reduces the operation time and storage space.
Inversion of potential data indeed faces serious non-uniqueness problem, which is a big challenge. There is no doubt that the non-uniqueness affects the accuracy and reliability of inversion results 45,46 . With m observation points data and the subsurface divided into n mesh cells, the proposed method has the similar degree of non-uniqueness as the traditional inversion methods [47][48][49][50][51] . The modulus difference anomaly inversion has the same unknow number of parameters (n) with the magnetization inversion. In other words, the non-uniqueness degree of the proposed method does not increase drastically compared with previous work. We add a simple magnetization intensity boundary constraint. Furthermore, we use the preconditioner instead of model constraint. The preconditioner 34 plays the same role with depth weighting 52 , acting as a depth-weighted function in recovering the magnetization intensity distribution. It can also lower the condition number of coefficient matrix and increase the convergence rate.
where the preconditioner P = z β I . z is the depth of divided cell and β is a constant. The proper value β ranges from 4 to 6 in 3D inversion pointed by Liu 34 . It is concern with the distance between observation and mesh cell. I is a unit matrix. The preconditioned conjugate gradient algorithm follows that of Liu and Hu 53 . First, the initial model (m 0 ) is given. An appropriate model is of advantage to reduce the uncertainty. Then, the correction (Δm) of the model's parameters is calculated. The model's parameters are updated by solving Eq. (13) with the preconditioned conjugate gradient method. The optimal solution is obtained after repeated iterations.

Synthetic examples
Cuboid model. In this experiment we conducted two tests-one with moderate and one with high magnetization intensity. The source was a cuboid model, 200 m × 150 m × 80 m, with its top buried at a depth of 50 m. The total inclination and declination were 45° and 0°, respectively. Figures 10a and 11a show the magnetic anomalies computed for m = 80 A/m and 110 A/m, respectively; the amplitudes of the total-field magnetic anomalies were 4,691 nT and 12,633 nT, respectively. Here, the magnetic anomaly is not displayed by the projection of the magnetic anomaly vector on the geomagnetic field but by the modulus difference.   Fig. 10b-d. Figure 10b, c shows the horizontal cross-section at a depth of 90 m calculated using ΔT pro and ΔT true , respectively. Figure 10d, e shows the vertical cross-section at easting = 500 m calculated using ΔT pro and ΔT true , respectively. There is no significant difference between the inversion results calculated by the two methods.
The inversion results for m = 110 A/m are shown in Fig. 11. The maximum misfit data of ΔT true is 36 nT and the maximum misfit data of ΔT pro is 159 nT. Figure 11d, e shows the vertical cross-section at easting = 500 m calculated using ΔT pro and ΔT true , respectively. Figure 11b, c shows the horizontal cross-section at a depth of 90 m calculated using ΔT pro and ΔT true , respectively. The red lines show the position of the true model. Here, the  www.nature.com/scientificreports/ differences between the two datasets are clear. The inverted magnetization intensity distributions based on ΔT true (Fig. 11c, e) are in good agreement with the true model and well describe the boundary of the target body. The inversion result shows it has favourable symmetry capability. However, the results based on ΔT pro (Fig. 11b, d) are distorted. The results of the northing end show some fake anomalies at cross-section. The inverted magnetization intensity boundary is distorted; hence the recovered depths do not agree with the true depths. The inversion result based on ΔT true are much closer to the source than those based on ΔT pro . the dipping dike model. The dipping model is shown in Fig. 12a. Without considering the magnetic remanence, the total inclination and declination of the magnetization vector are 45° and 0°, respectively. The total magnetization intensity was set as 150 A/m and the geomagnetic field is 50,000 nT. The observed data are arranged in a 10 m × 20 m grid comprising 5,151 cells (Fig. 12b). The subsurface is divided into regular cells of size 25 × 25 × 25 m in the X-, Y-, Z-directions. The inversion results are shown in Fig. 12c-f. Figure 12c, d repre- www.nature.com/scientificreports/ sents the 3-D inversion magnetization intensity distribution based on ΔT pro and ΔT true , respectively. Figure 12e, f shows the vertical cross-section map at northing = 500 m derived from the 3-D inversion calculated using ΔT pro and ΔT true , respectively. Figure 12c, e shows the magnetization intensity distributions by calculated using ΔT pro . The inverted results show the target has a buried inclined geological structure. However, it differs from the true model, which we can see in the cross-section and 3D result. The map of the vertical cross-section at northing = 500 m (Fig. 12e) shows there is no magnetization at easting of 600-800 m with the depth of 350-450 m. it results in a significant discrepancy between the shapes of the inversion results and the model. Conversely, the boundary of the dipping dike is well displayed by magnetization intensity distributions of ΔT true compared with those of ΔT pro . The inversion results based on ΔT true accurately reproduced the shape of the source (Fig. 12d, f). and were generally consistent with Fig. 12a except for a little deviation at the depth 350-450 m caused by the deeply buried magnetic source. It is a common scene for inverted property to have a trend concentrating near the surface in potential field.
In terms of shape recovery and characterization of targets, modulus difference anomaly has advantages, which can better describe boundary of magnetic bodies at the same situation. It can show some characteristics of the model and realistically restore the shape of the model. These results show a marked improvement over the ΔT pro inversion results in the intensity, shape and position of the source.

Field example
The Mengku iron deposit is one of the largest iron deposits in Xinjiang Province, China. It is located in the Aletia active margin of the Siberian plate and the central region of the north-eastern margin of the Late Palaeozoic Maizi inland rift basin. The geomagnetic field intensity is T 0 = 58,110 nT, with a geomagnetic inclination I 0 = 67° and declination D 0 = 3°. The main rock is amphibole oligo gneiss formed by regional metamorphism, which is locally known as fault migmatite. It is a polygenetic composite iron deposit formed by layered deposition, metamorphism and magmatic hydrothermal superimposition. Magnetite is the main mineral resource in the Mengku deposit. The rich magnetite has a strong magnetism with the susceptibility of 0.8-1.8 SI and the magnetization intensity is ranged from 40 to 90 A/m 35,54 . The magnetite ores show an amplitude change of − 6,000 nT from 22,749 nT, which provides data that are very suitable for testing inversion methods. The total-field magnetic anomaly is shown in Fig. 13. There are also some drillholes in line 135 and 143. Additionally, studies have shown that Mengku iron ore has strong remanence; the ratio between the induced magnetization and remanence is about 2.3 and the average residual ratio of its surrounding rocks is about 2.8 35 . Remanent magnetization changes the direction of the total field magnetization. Therefore, we used the correlation between the total magnitude anomaly and RTP to obtain the average magnetization direction 55 , where I = 68° and D = 85°. After obtaining the total field magnetization direction, it is used in the inversion of both projection anomaly and modulus difference anomaly to reduce the effect by remanence. We recovered the 3-D total magnetization intensity distribution. The subsurface was divided into 50 × 50 × 50 regular cells, with all the cells having the same magnetization intensity. A boundary constraint for the magnetization intensity of 0 A/m ≤ M ≤ 100 A/m was applied in the inversion. The inversion results based on ΔT true and ΔT pro are shown in Fig. 14a, b, respectively. www.nature.com/scientificreports/ The magnetization intensity calculated using the ΔT true inversion mainly ranges from 30 to 80 A/m (Fig. 14a), which is in good agreement with the magnetic properties of field samples (from poor magnetite orebodies to rich magnetite orebodies) 54 . However, most of the inverted magnetization intensities based on ΔT pro lie within the range of 25-50 A/m (Fig. 14b). This again demonstrates the advantage of using ΔT true to calculate the inverted magnetization intensity distribution. Figure 14 shows that the ΔT true -based anomalies have almost continuous highly-amplitude values while the ΔT pro -based inversion results have intermittent low-amplitude magnetization. At the southeast end of the model, there is an anomaly of almost 6,000 nT (Fig. 13), which is not well reproduced by the ΔT pro inversion (Fig. 14b).
According to the information from bore logs of drill holes, what's more, we compare the difference inverting for between ΔT true and ΔT pro on the vertical cross-section between line 135 (Fig. 15) and 143 (Fig. 16). The results show inverting for ΔT true (Fig. 15a), the presence of more than one magnetite belts in NE direction. However, the results from using ΔT pro show significantly poor magnetization intensity (Fig. 15b), which is out of accord with the records of drillholes ZK135-2 and ZK135-3 (Fig. 15c). Thus, using ΔT pro does not produce a good representation of the deposit shape and reserves, which leads to an inaccurate interpretation of the data. The drillhole ZK143-1 probes the magnetite belt at depth of 67-178 m and ZK143-4 do not intercept magnetite ores (Fig. 16c). The profile of line 143 inverted by ΔT true (Fig. 16a) shows there are magnetite ores from 50 to 250 m with high susceptibility. Figure 16b shows the result using ΔT pro , the depth of the magnetite ores does not www.nature.com/scientificreports/ exceed 40 m, which is less than the results intercepted by drillholes. As discuss above, it indicated that directly inverting for ΔT true produce more close results with the depth, shape and occurrences of mineral bodies. We therefore conclude that using ΔT true in the inversion scheme to recover the magnetization distribution and location of orebodies will produce more reliable and acceptable results.

conclusions
For a strong magnetic anomaly, using projection anomaly to replace the modulus difference anomaly will yield errors. The error from using the projection anomaly as an approximation comes from two elements: the value of the magnetic anomaly and the magnetization direction. Our numerical modelling shows that the error increases exponentially as the magnetic anomaly increases, which is far beyond the precise magnetic allowance and no longer acceptable for high-precise interpretation; the influence of the magnetization direction on the error is mainly in the magnetic inclination, especially without remanence magnetisation. Furthermore, we found that when the error is the largest, the direction between the geomagnetic field and magnetic anomaly vector satisfies the mathematical relationship cos θ * = −T a /2T 0 . Through synthetic and field examples, we found that magnetic data inversion based on the projection anomaly do not produce a good fit with geological bodies in strong magnetic environments, while directly inverting for modulus difference anomaly data yield a more accurate representation of both the shape and location of the magnetic source, which is helpful to promote the horizonal and vertical resolution and reduce the ambiguity in the highly magnetic environment inversion. www.nature.com/scientificreports/ As the requirements of higher accuracy in data processing and interpretation of magnetic sources are rising, the error of the magnetic anomaly should be reduced. Hence, using projection anomaly data is not appropriate. In cases of high-amplitude magnetic anomaly, we recommend using modulus difference anomaly data for reliable results.