Two-Dimensional Compact Inversion of Magnetic Data in the Presence of Remanent Magnetization

Remnant magnetization causes a change in the direction and intensity of the magnetization vector. If inversion is performed regardless of remnance, in some cases it may have unreliable and misleading results. For inversion with respect to remnant magnetization, several solutions have been proposed so far, one of which is to convert the data of total magnetic field into data that is independent of the direction of magnetization. In this study, the transformation of Total Field Anomaly (TFA) into Total Magnitude Anomaly (TMA) is used. The inversion algorithm is based on improving compact inversion method and is just two-dimensional. In compact inversion, anomalies may concentrate on the surface of the earth, and thus the response is unreliable. To solve this problem, a combination of matrices and weighting functions have been used, including elements such as magnetic susceptibility and depth function. The resulting model can be smooth or compact (with sharp edges) based on changing compactness factor. The method has been tested using several synthetic and real data. The synthetic data are a 2D tabular prism, of which the top buried-depth is 20 m and the length and width are 40 to 20 m, a dipping prism with a vertical tabular nearby. The real example is magnetic data over Galinge iron-ore deposit in Qinghai province of China, and the data of four profiles have been considered for 2D inversion. Inversion even smooth or sharp, have been conducted with all models, and especially sharper models are consistent with the known geologic attributes of the magnetic sources. Remnant Magnetization


INTRODUCTION
In some magnetic observations, it can be assumed that there is no remnant magnetization or can be neglected. In these cases, the direction of remnant magnetization is assumed to be parallel with the direction of the Earth's magnetization, and some modeling can be done with this hypothesis. But there is remnant magnetization in most cases, which, if severe, causes uncertainty in inversion of magnetic data. Although this problem was present from the beginning of the magnetic measurement, it was not important for two reasons: first, in many magnetic studies, the amount of remnant magnetization was considered to be negligible; secondly, many magnetic data were interpreted qualitatively; It was performed only for the purpose of estimating depth, in which there was no need to know the direction of magnetization [1]. However, remnant magnetization is severe and inevitable in many applications of magnetic method, including mineral exploration, regional explorations of the earth's crust and archeology. Given the above mentioned reasons, inversion result is unreliable without attention to the remnant magnetization; so many researchers have offered numerous ways to solve this problem. In general, there are three methods for inversion of magnetic data with respect to remnant magnetization: 1) inversion method with estimation for the direction of remnant magnetization, 2) inversion by converting magnetic data into the data that is independent of magnetization and 3) Magnetic Vector Inversion (MVI).
If the direction of remnant magnetization is somehow estimated or calculated, it can be easily inverted by replacing it into the forward equations. This method is the simplest and most precise method for inversion of magnetic data [2,3,4,5]. Li et al [6] have used this method to model magnetic data in the North America. There are several methods for estimating the direction of magnetization, which can be done by Helbig's method [7,8,9,10], wavelet multiscale edge method [11], And Cross-Correlation Method [12,13].
If the magnetic data are converted to another kind of data so that it is independent of the direction of magnetization, modeling becomes easier and more reliable. As we know, the combination of vertical and horizontal derivatives of the magnetic field in two dimensions is called Analytical Signal (AS), which is introduced by [14]. The analytical signal value, which is same as vector magnitude of magnetic gradients, is almost independent of the magnetization and can be used throughout the profile to estimate the depth of the anomaly without attention to remnant magnetization. Shearer and Li [15] and Shearer [1] have developed a method based on analytical signal data (total gradient vector magnitude) and taking into account the inversion with positivity constraints. If it is possible to estimate the components of the magnetic field in three directions of Cartesian or to measure three components of the magnetic field (in some aeromagnetic measurements, each of the three components of the magnetic field is taken), it is possible to measure Total Magnitude Anomaly (TMA). The general method for obtaining three components of the magnetic field is the use of the wavenumber method [16]. Inversion of TMA data extensively is used by Li et al [6] and twodimensionally by Liu et al [17,18] for ground magnetic data and magnetometric drilling logs.
The inversion method with estimation of the remnant magnetization direction, despite its simplicity, cannot be used in many cases or is not reliable. Because the direction of remnant magnetization may not be same throughout the magnetic body, or due to tectonic activity (e.g. faulting), a part of the body is displaced or rotated, so there is more than one magnetization direction.
According to the limitations mentioned above, another method called magnetic vector inversion is proposed by Ellis et al [19] and MacLoad and Ellis [20]. Both the magnitude and the magnetic vector are estimated as two unknown parameters. This approach has been proposed with some differences by Kubota and Uchyama [21], Lelivre and Oldenburg [22] and Pratt et al [23]. The method dramatically reduces the risk and error as compared with estimation of the direction of remnant magnetization methods, and is a suitable method for inversion of local and regional magnetic data, but due to increasing model parameters triply, the non-uniqueness of this method will increase.
Last and Kubik [23], in particular, developed a method for inversion of gravity data called 'compact inversion', explaining the observed anomaly by structures of minimum volume. Their strategy is to minimize the area of the model, so this is equivalent to maximizing its compactness. This method strongly depends on the number of model parameters, and by increasing the model parameters, the anomalies tend to concentrate near the surface. Guillen and Menichetti [25] applied an approach which includes the search for solutions, minimizing the moment of inertia with respect either to the center of gravity or to an axis of a given dip line passing through it. The method works properly for a single gravity source, but the problem is dealing with multisource and complicated anomalies which do not lie in one point or one axe. Barbosa and Silva [26] generalized the compact inversion method to compact along several axes using Tikhonov's regularization. They improved the method offered by Guillen and Menichetti [25] for multisource and complicated anomalies. The most recent compact method is called 'interactive inversion' that was developed by Barbosa and Silva [27] for magnetic data, which estimates the location and geometry of several anomalies. They simplified their old method [26] for computational performing. The method is suitable for multi-source and even complicated anomalies depending on the quantity and quality of 'a priori' information. Silva et al [28,29] also developed interactive inversion successfully for 3-D gravity data closest to pre-specified geometric elements such as axes and points. Ghalehnoee et al [30] have provided a more comprehensive and general method compare with the above methods for two-dimensional gravity inversion. They introduced three weighting functions such as depth function, compactness function, and weighting function based on the kernel matrix combined together to solve the problem.
We modified Ghalehnoee et al [30] method for inversion 2-D magnetic data. We changed weighting functions of the gravity method to a depth weighting function and compactness weighting function for magnetic method, therefore, the kernel weighting function has been eliminated compared with the gravity inversion algorithm due to lower lateral dependence than the gravity method. The elements of compactness weighting matrix are the susceptibility or magnetization of the model. The Total Field Anomaly data (TFA) also converted to Total Magnitude Anomaly data (TMA) to reduce remnant magnetization effect, so the algorithm runs on TMA data.

METHODOLOGY
Total magnitude anomaly (TMA) data is calculated from the total field magnetic anomaly (TFA) through a Fourier domain calculation or equivalent source by converting to the three orthogonal components of the total field anomaly, Bx, By, and Bz, and then calculating the magnitude of the quantities, B: For two-dimensional data, By will be eliminated. To find TMA data, the components of the total field anomaly, Bx, By, and Bz, are calculated in the wavenumber domain [16] or by the equivalent source technique [31]. If the ground beneath data measurements is subdivided into two-dimensional prisms ( Fig. 1), the total magnitude anomaly (TMA) in the ith data is [24]: Where is the susceptibility (SI) or magnetization intensity (A/m) of the jth prism, is the noise associated with ith data point, and is kernel matrix, the elements of which represent the influence of the jth prism on the ith magnetic value.
The data equation can be rewritten in matrix notation (3) The inversion method here is linear, like other linear geophysical inversions: given, the observed magnetic data (B), finding a density distribution 'm' which explains 'B', with a certain noise level. The solution of the system in kth iteration can be in the least-square problem in matrix notation [32]: Where is model weighting matrix, is a noise weighting matrix in kth iteration, both of which are diagonal, and . Mu ( ) is damping factor or regularization parameter to get rid of matrix singularity, having the positive value and depending on the noise level of the data measurements. The less value of damping factor refers to the less noise level of the data. The most common and simplest method for estimating damping factor is the use of L-curve. The L-curve is a log-log plot of the norm of a regularized solution ‖ ‖ versus the norm of the corresponding residual norm ‖ ‖ . It is a convenient graphical tool for displaying the trade-off between the size of a regularized solution and its fit to the given data, as the regularization parameter varies. The L-curve thus gives insight into the regularizing properties of the underlying regularization method, and it is an aid in choosing an appropriate regularization parameter or damping factor for the given data [33]. After fitting the curve on the data, the point with the greatest curvature is considered as the desired value of the damping factor. Figure 2 schematically shows L-curve for different values of damping factor; the optimum value for damping factor is the greatest curvature in the curve (red circle in Fig. 2).
We introduce a general weighting function including the compactness weighting Wc and depth weighting Wz matrices in kth iteration: (6) Noise weighting matrix , can be simply written as [24]: Magnetic data, like other potential field data, have no inherent resolution. When minimizing, ‖ ‖ ∫ , final model tends to concentrate near the surface regardless of the true depth of the causative bodies. This arises because the constructed model is a linear combination of kernels that decays sharply rapidly with depth [34].
A depth weighting has been presented here to be appropriate for 2D magnetic data: Where is the depth of the jth prism, is the measurement height for each data from surface, and β is a constant value.
The value of is usually chosen to reproduce the exponential decay of the magnetic response of a sphere with distance [35]. Li and Oldenburg [3] estimated the value of β = 3 and Liu et al [18] used . Accordingly, the value same β = 3.5 has been chosen in this study.
The compactness functional (minimum area or minimum support) was first introduced by Last and Kubik [24], who suggest seeking the source distribution with the minimum area to model the anomaly. This concept is illustrated as: if d and h are the prism dimensions, a definition of area for 2D model is [24]: This leads to choose the compactness weighting function (11) Where is compactness factor and the parameter is a small number of 10 -10 which is introduced to provide stability as . This weight is not appropriate when reference model exists due to geological information, therefore, further developed by Portniaguine and Zhdanov [5] who used the term "minimum structure". They added a reference model (m o ) to Eq. (11); therefore, the final compactness function is yielded | | The reference model m o may be a general background model that is estimated from previous investigations, or it could be zero (null).
The value of has been chosen as a constant value by many authors. Last and Kubic [24] and Vatankhah et al [36] presented ; Guillen and Menichetti [25], Barbosa et al [26], Barbosa and Silva [27], Silva et al [28,29] and Grandis and Dahrin [37] have chosen . In this study the value of is variable and depends on the situations coming from the depth of causative bodies and prior information. When the source lies at large depths, the large value of can be opted and vice-versa. If there is a subsurface geological information such as drilling or known minimum and maximum bounds, the value of compactness factor can also be large to have a sharp model. Finally, it should be pointed out that the appropriate range of value is for 2D TMA data inversion.

Inversion Procedure
For the beginning, if there is no priori information, the reference model is chosen , then Wm (Eq. 6) and We (Eq. 7) all of which are calculated, based on value of compactness factor. Then, we will have iteration procedure with least-square solution as Eqs. (4 and 5). The inversion process is shown as a flowchart in Figure 3. The procedure should be performed for many damping factors to plot Lcurve, thus yielding the model with the best value of damping factor.
Besides, the approach explained here, does not depend strongly on the number of parameters (M), unlike compact inversion; it is advised to have fewer number of parameters by increasing the prism dimensions as far as possible thus having the most appropriate model and fast computational running.

Bound and Positivity Constraints
To produce a physically meaningful model, refer to prior information obtained from geological maps, well-loggings, and laboratory tests, the susceptibility of each prism must satisfy (14) Where is lower and is upper bound. If at an iteration mj exceeds one of its bounds, then it will be fixed at the violated bound. Instead of being calculated by Eqs. (6 and 7), the corresponding weight wj will be set at a very large value. The large value of wj will force the prism density contrast to be 'frozen' at the violated boundary, at least temporarily during the first few iterations. This way, the response of the modified parameter estimates will not fit the observations, which will in turn trigger the necessity for further parameter perturbation as a function of the misfit [25,26,37]. For positivity constraints, the lower bound should be equal or close to zero.

INVERSION OF SYNTHETIC DATA
To demonstrate the capability of the method, the algorithm is tested with two synthetic data. For this purpose, a simple model including a tabular prism with the top depth of 20 meters, width and length, respectively, of 20 and 40 meters is considered first (Fig. 4). Effective susceptibility k = 0.23 SI and Earth's magnetic field intensity are 50000 nT with a 45 o angle of inclination (Ii). Therefore, the induction magnetization Mi = 4 A/m is parallel to the direction of the Earth's magnetization. It is assumed that the source has the remnant magnetization with the magnetization of Mr = 6 A/m. The remnance inclinations (Ir) are considered 0 and 90 degrees. As shown in this figure, the TFA value is different at these two inclination angles, but the TMA value is completely constant, so using TMA data for inversion provides more reliable results than TFA data. If TFA data is used, it is necessary to use both Ir and Ii angels correctly. If this is neglected, assuming that there is no remnance, the inversion may not be acceptable. Figure 5 shows the inversion results of TFA data regardless of the remnant magnetization, that is, the only geomagnetic inclination has been included. As seen, all two models differ from true anomalies in both data set of TFA data for remanence inclination of 0 and 90 degree.

Figure 5. The inversion results obtained from the TFA data, regardless of the direction of remanent magnetization when the direction of the remanent magnetization is a) 0 and b) 90 degrees; the rectangle shows the true magnetic source.
The effective susceptibility (SI) can be directly computed by [22]: (15) Where keff is the effective susceptibility, M is total magnetization, To is the intensity of the Earth's induced magnetic field and μo is permeability of free space ( ).
Inversion of TMA data is performed using the present algorithm with square prism of 5 meter dimensions. Fig. 6 shows the inversion results using a damping factor of 0.1 and compactness factor of 0, 0.5, and 1 after 30 iterations.
As seen in the figure, with the highest compactness factor (α = 1), the anomaly edges are sharper and the actual depth is also more accurately modeled, while reducing the value of compactness factor, the depth of the anomaly becomes deeper and the amount of calculated susceptibility is reduced. For the case where the sharp edges or blocky model are considered, the compactness factor should opt a large value, and also it is necessary to perform inversion using bound constraining. In this example, for all models, the positivity constraining is used; also, for α = 1, the upper bound constraining is considered 0 ≤ m ≤ 0.23 SI. The second example is related to two anomalies, including a dipping prism with a dip of 45 ° and a nearby vertical tabular. The intensity of the earth's magnetic field intensity is 50,000 nT, the effective magnetic susceptibility is considered 0.112 SI (4.5 A/m) and the geomagnetic inclination, and remnant inclination are 45 ° and -45 °, respectively. Fig. 7 shows the noise-free data and noise-added data with Gaussian noise level of 10%. The inversion of noise-free data was performed with square-cell of 5m dimensions using compactness factor of 0, 0.5 and 1 by damping factor of 2, after 10 iterations (Fig. 8). As shown in Fig. 8, by increasing the value of compactness factor, the model is more similar to its actual anomaly and also the modeled magnetic susceptibility is closer to the real value.  Therefore, for inversion of TMA data when there is geology information, the compactness factor is suggested to be α = 1. Thus, for the next inversions, only compact factor of 1 is used.
Data is also affected by 10% Gaussian noise. In real data measurements we often have some noise and therefore, modeling such data will encounter problems. If inversion is done regardless of noise, unwanted anomalies may be created without geological origin, or the models may not be in actual depth or differ with real geometry; in other words, uncertainty may occur. The reliability of the inversion with noisy data can be done by selecting the best value of damping factor or regularization parameter. If a small value of damping factor is taken, the effect of noise is strongly demonstrated, and models with no relation to geology may be yielded, and if damping factor is large, then smaller anomalies may not be modeled, and the bodies may be modeled greater than their actual depths. Therefore, using L-curve is the simplest way to select the best value of damping factor. According to the above-mentioned statements, the inversion of noisy data has been done with the same cell dimensions and 10 iterations using different damping factor separately with only the compactness factor of 1. Fig. 9a shows L-curve for different values of µ. As seen, the best value of damping factor is 35, which compared with noise free data (µ = 2) is a much larger. The inversion for a small damping factor of 8 and a large damping factor of 90 is also done (Fig.  10). As shown in Fig. 10, by using a small amount for damping factor, the effect of noise strongly affects the model and creates unwanted anomalies with unrealistic geometry, but for a large value of damping factor, the bodies are modeled at more than their actual depths.

INVERSION OF REAL DATA
In synthetic examples, the generality and comprehensiveness of the method for anomalies with different geometric shapes and different depths was discussed. In this study, the data of the Galinge Iron Ore deposit are in focus.
The Galinge Iron Ore deposit is located in Qinghai Province, northwest China (longitude: 92° 09' 00″-92° 19' 00″ east and latitude: 37° 04' 00″-37° 08' 00″ north). The deposit contains 117 to 210 meters overburden of quaternary gravel sediments, and therefore the mineral veins are relatively deep. A total of 16 exploratory boreholes have been drilled up to a maximum depth of 600 meters. Borehole drillings show that there are a total of 8 iron belts in the Ordovician rocks of the Tanjianshan group (Fig. 11). This deposit contains dispersed and dense gray-black magnetite with some pyrite and limonite. The surrounding rocks as Tanjianshan group, include the mud siliceous rocks, diopsidite intercalated with clastic rocks, intermediate volcanic rocks, and marble [18]. The magnetic properties of the collected samples show that the susceptibility of iron ores are much higher than the surroundings and significant remnant magnetization occurs. The susceptibility of iron ore in drilling samples is between 0.5 and 6.25 in SI units, and the remnant magnetization is from 20 to 60 A/m. Susceptibility of the surrounding rocks is less than 0.011 SI. The ratio of remnant to induced magnetization varies from 0.24 to 1.26, with the largest variation related to dispersed mineralization and the least value related to massive mineralization. The direction of remnant magnetization is undetermined due to difficulty to collect the oriented samples in the boreholes [18]. Figure 12a shows the residual magnetic data (IGRF is removed) of Galinge deposit which varies from 400 to 1600 nT. As it can be seen, the general mineralization trend is almost east-west.
Liu et al [18] proved that the difference between inclination of remnant and Earth's magnetization is 50 degrees, and therefore, to reduce the remnance effect, TFA data is converted to TMA data (Fig. 12b). The magnetic anomaly of the TMA data is in agreement with the image of the ore deposit on the surface. For inversion, TMA data of four sections: p196, p204, p212 and p220, are being prepared along with drilling results. Fig. 11 shows the geological sections of these profiles. As it is shown, the Galinge iron deposit consists of several iron belts with a high degree dip to the southwest.
Inversion is done using 25 meter square cells with compactness factor 1 and different damping factors for all four sections in 20 iterations. Fig.  13 shows the L-curve of the sections. In this figure, the most appropriate value of damping factor corresponding to the highest curvature is different for each of the four sections, so for p196, p206, p212 and p220 the best damping factor values are 3, 6.75, 0.58 and 2 respectively. Fig. 14 shows the results of inversion. Since α = 1, the models are fairly sharp, this value is used for all sections. It is worth noting that due to the close proximity of the iron belts, inversion cannot separate the belts from each other. In the inversion results, the maximum susceptibility is 0.233 SI, while the average susceptibility of the ore is approximately 1.15 SI from the drilling data. The reason is that the mineralization has occurred in the case of alternating belts, and the modeled value of susceptibility is the result of the average mineralization and surrounding rocks, and so it is natural that the susceptibility of the final model decreases. Fig. 15 shows the TMA data of all four profiles with the corresponding response that are fitted. In this figure, it is seen that the amount of data noise is low and there is a good fit between observational data and calculated data.

CONCLUSIONS
According to the above statements, the algorithm used here is to improve the compact method using a depth weight matrix to obtain the distribution of the magnetization intensity or susceptibility of the magnetic anomalies. In this method, the remnant magnetization is taken into consideration, and in order to reduce the effect of remnance, the TFA data is converted to TMA data, which is almost independent of the magnetization direction. In this method, the compactness factor is introduced, which is determined based on the depth of the magnetic field or the availability of geological information. If the depth of the body is large or the exploration information is such that at least the bound constraint can enter into the model, then the compactness factor is chosen great value and vice versa, if the anomaly is at a low depth, or there is less geological information, it should be used a lower or zero compactness factor to make the model more reliable. This method can also be used for three-dimensional inversion with some changes, with the fact that the inversion of 3D data is more time consuming than the twodimensional inversion due to the increasing number of parameters and data.  L-curves (Fig. 13); the polygons show the mineralization belts. Figure 15. The fit between observed TMA data (obs) and calculated (cal) of the four sections p196, p204, p212 and p220; the observed and calculated data show with the solid circle and solid line respectively.