Intersection Constraint Weighting (ICW) Method: High-Resolution Joint Magnetic Susceptibility Inversion of Aeromagnetic and Gradient Data

: Aeromagnetic gradient data have higher horizontal resolution on shallow sources, a magnetic anomaly that can better reﬂect the features of deeper sources. Therefore, we used a joint survey of aeromagnetic and gradient data to obtain the distribution of sources with different depths. In this paper, we propose an intersection constraint weighting (ICW) method based on aeromagnetic and gradient data for original and gradient anomalies with inconsistent high-low variation characteristics. The ICW method can effectively improve the resolution of inversion results and can more accurately obtain the distribution of magnetic bodies via cross-gradient by gradually adding a gradient component and applying a normalized property weighting function. Our theoretical model tests indicate that the distribution of the recovered magnetic susceptibility model of the ICW method was similar to that of the true model. In addition, the anomaly containing noise with different signal-to-noise ratios veriﬁed that the ICW method had a stronger anti-noise ability than other methods. We also inverted real data in the Zhurihe area of the Inner Mongolia Autonomous Region in northern China. The inversion result showed that the main trend of high magnetic bodies was in the northeast direction, that the shallowest depth of high magnetic bodies was 100 m, and that the greatest depth was 960 m.


Introduction
Metal mineral resources generally have high magnetic susceptibility, so they appear as high-value anomalies in magnetic exploration.Aeromagnetic surveys are unrestricted by terrain, fast and economical, and widely employed in practice [1,2].The magnetic susceptibility inversion results of magnetic anomalies can be used to delineate the distribution of magnetic ore bodies [2][3][4][5][6][7][8].Bott [9] applied the matrix linear solution to determine the magnetic susceptibility distribution.Arkani-Hamed and Stranqway [10] developed a technique to invert a scalar magnetic anomaly map into a magnetic susceptibility anomaly map.Well-defined susceptibility highs are associated with the Kursk iron formation and small-scale shields and plates, and well-defined susceptibility lows delineate collision zones in the Alpine-Himalayan belt.Li and Oldenburg [11] elaborated the regularization solution method of 3D physical property inversion of magnetic data and obtained a 3D magnetic susceptibility distribution by introducing model constraints and a depth-weighting function.Pilkington [12] presented a preconditioned conjugate gradient approach to solve the 3D magnetic inversion problem with the least amount of prior information to determine the magnetic susceptibility distribution of a given magnetic anomaly.Li and Oldenburg [13] developed a fast inversion method based on wavelet compression and a logarithmic barrier of minimization to recover the 3D distribution of magnetic susceptibility.Pilkington [14] introduced the 3D sparse inversion of aeromagnetic data combined with the conjugate gradient algorithm to improve the resolution of the inversion results.Li et al. [15] developed a level-set approach to determine the magnetic susceptibility value through joint inversion of surface and borehole magnetic data, which can be utilized as a probing parameter to assess the uncertainty in the spatial extent of causative bodies.
Different gravity gradient tensor components can enhance different characteristics of the field source [16].Different components have different resolutions for the source, and due to the lack of a calculation method, different components will produce different inversion results [17][18][19].Certain scholars have improved the effect of model reconstruction by adding anomaly components during inversion [20,21].Qin et al. [22] directly integrated different gravity gradient components into a matrix to improve the resolution of inversion results.Magnetic gradient data have a higher horizontal resolution than magnetic anomalies, effectively highlighting shallow geological bodies, while original magnetic data can more effectively retain deep field source information.Therefore, the joint exploration of aeromagnetic fields and their gradients is commonly performed to delineate metal mineral resources [23][24][25][26][27][28].
The conventional joint inversion method is data joint inversion.which superimposes the original magnetic anomaly onto the magnetic gradient anomaly in a matrix to improve the precision of inversion results using weight constraints of different anomaly components [29].Joint inversion based on structural constraints is another way to improve inversion resolution, such as the cross-gradient method, which is mostly utilized for the inversion of different types of data [30][31][32][33][34][35][36].The structure factor of the cross-gradient has a better constraining effect on the magnetic and magnetic gradient data of the same source and different resolutions, which can improve the inversion precision and satisfy the consistency of the magnetic and magnetic gradient inversion results of the final results of the ICW method.
To improve the resolution and precision of the joint inversion of aeromagnetic and gradient data and to better reveal the distribution of the magnetic sources, this paper proposes an intersection constraint weighting (ICW) method.The ICW method involves a data constraint strategy, model weighting function design, and a step-by-step, cross-gradient constraint calculation process to obtain higher resolution inversion results.The precision and practicability of the method are proven by theoretical model tests and real data application.

Methodology 2.1. Data Joint Inversion
Magnetic susceptibility inversion divides the underground into grids according to the measurement point distance, calculates the magnetic susceptibility of each grid, and then obtains the distribution of the underground field source.The number of divided grids generally exceeds the number of observation points.There are two main ways to realize the joint magnetic susceptibility inversion of aeromagnetic and gradient data.A more commonly employed inversion method is putting different parameter data in a matrix to achieve magnetic susceptibility constrained inversion, which is named data joint inversion [29].The objective function is where d ∆T is the magnetic anomaly and d ∆T x , d ∆T y , and d ∆T z are the magnetic gradients in the x-direction, y-direction, and z-direction, respectively.The matrices G ∆T , G ∆T x , G ∆T y , and G ∆T z are the forward operators that connect the observation data with the physical parameter [12].The column vector m 2 represents the physical parameter (magnetic susceptibility) inverted by the combined matrix d ∆T ; d ∆T x ; d ∆T y ; d ∆T z .δ is the regularization parameter, and we use the L-curve method to determine the regularization parameters [37].W m DJ is the model weighting matrix [38,39], and W d DJ is the data weighting matrix [11].

Cross-Gradient Inversion
cross-gradient inversion is another method that can be utilized for joint magnetic susceptibility inversion of aeromagnetic and gradient data by using the gradient change in physical properties as structural constraints to improve the resolution of inversion results.The objective functions are where γ is the coefficient of the cross-gradient term, Φ cg = t(x, y, z) 2  2 is the cross-gradient term [31], and t(x, y, z) is calculated by t(x, y, z) = ∇m ∆T × ∇m ∆T g = t x , t y , t z , where

, y) ∂y
To verify the resolution and depth information of different components of magnetic and gradient data to different depths, we introduce depth resolution plots (DRPs) [40] and eigenvalue spectra for kernel function matrices [16].The source volume has dimensions of 3 × 3 × 1 km, which is discretized with n x = n y = 30 and n z = 10; that is, the solution domain consists of n = 9000 cells.The data grid covers an area of 3 km × 3 km, and they are measured at ground level.For the first case, which uses a combination of magnetic gradients in the x-direction, y-direction, and z-direction, the data points p of each gradient anomaly are arranged in a 30 × 30 grid (p = 900).For the magnetic anomaly component, the data points are arranged in a 54 × 50 grid (p = 2700).
Consider a model that contains two prisms with the same burial depths, where the top burial depth and bottom burial depth of prism 1 and prism 2 are 300 m and 500 m and 300 m and 400 m, respectively, and the real magnetic susceptibility is 0.25.The inclination and declination are 45 • and 60 • , respectively.
Figure 1a shows a 3D view of the model and the corresponding magnetic anomaly with a flight height of 100 m, and Figure 1b-d shows gradient anomalies in the x-direction, y-direction and z-direction.The dotted line indicates the real position of the geological body.
Figure 1e,f are the DRP of the magnetic and combination matrix of gradient components for the first case.As the depth increases, the depth resolution decreases, and when the depth is shallower, the resolution of the combination matrix of gradient components is relatively high.Figure 1e-g shows the corresponding eigenvalue spectra for the kernel function matrices.In the first 900, the depth information carried by the combination matrix of gradient components is slightly lower than that of the combination matrix of magnetic and gradient components.In the latter 1800, the depth information carried by the combination matrix of gradient components decays more slowly.It can be seen from the inversion results in Figure 1h-j that the resolution of the combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ) inversion is worse and that the magnetic anomaly (d ∆T ) inversion results are similar to the combined matrix ( d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ) inversion, indicating that the addition of gradient components cannot improve the resolution of the inversion results, which is consistent with the conclusion of Paoletti et al. [16].
For the second case, which uses a combination of magnetic gradients in the x-direction, y-direction, and z-direction, the data points of each gradient anomaly are arranged in a 62 × 50 grid (p = 3100).For the magnetic anomaly component, the data points are arranged in a 100 × 93 grid (p = 9300).In Figure 2e,f, we discovered that the depth resolution of the combination matrix of gradients is higher than that of the magnetic matrix when the depth is shallow, which is consistent with the conclusion of Figure 1.In Figure 1g, the depth information of the combination matrix of gradients decays the slowest, and the decay trend of the depth information of the magnetic component is consistent with the decay trend of the combination matrix of the magnetic and gradient components.The figure also shows that the depth information carried in each singular vector is almost independent of the number p.In Figure 2h-j, the resolution of the combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ) inversion is slightly improved with an increase in the number of data points, and the results of the magnetic anomaly (d ∆T ) inversion and combination matrix ( d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ) inversion are similar to those shown in Figure 1.In general, the addition of gradient components does not contribute to the improvement in inversion resolution with the same burial depth and similar characteristics of magnetic and gradient anomalies.Next, we will discuss the situation of different burial depths of subsurface geological bodies and different characteristics of magnetic and gradient anomalies.;; ) inversion by the data constraint joint inversion method.(j) Combination matrix (

;
; ; d d d d ) inversion by the data constraint joint inversion method.A model that contains two prisms with different burial depths, different sizes, and identical physical properties was designed to test the effect of the existing method.The top burial depth and bottom burial depth of prism 1 and prism 2 were 300 m and 500 m and 400 m and 800 m, respectively, and the real magnetic susceptibility was 0.25.The inclination and declination were 45 • and 60 • , respectively.
Figure 3a is a 3D view of the model and the corresponding magnetic anomaly with a flight height of 100 m, and Figure 3b-d  ).

Intersection Constraint Weighting (ICW) Method
To obtain higher resolution and precision inversion results, we propose an ICW Figure 3e shows a magnetic susceptibility slice (y = 1900 m) by regularization inversion of the magnetic anomaly.The result can better reflect the distribution of deeper geological bodies, but the resolution is lower.Figure 3f is the magnetic susceptibility slice (y = 1900 m) obtained by data constraint joint inversion of the combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ).The result better highlights the shallower geological body.It has a higher horizontal resolution, but the deviation between the inversion result of the deep field source and the real model is larger.Compared with Figure 3e,f, the inversion results of the deep geological body with the original anomaly are better than those of the gradient anomaly, while the inversion results of the shallow geological body with the gradient anomaly are better than those of the original anomaly, and the resolution is higher than that of the original anomaly.Figure 3g is the magnetic susceptibility slice (y = 1900 m) obtained by data constraint joint inversion of the combination matrix ( d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ), which further improves the effect of shallow geological bodies but remains poor for deeper geological bodies.Integrating different components into a matrix for inversion is an effective way to improve the resolution of the inversion results [22].
Figure 3h,i show magnetic susceptibility slices (y = 1900 m) obtained by cross-gradient inversion of the magnetic anomaly and combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ), respectively.Compared with Figure 3e,f, the horizontal resolution of deep and shallow geological bodies is higher.However, the distinct distribution of geological bodies of different depths still cannot be obtained simultaneously, and the value of magnetic susceptibility is still quite different from the true magnetic susceptibility.

Intersection Constraint Weighting (ICW) Method
To obtain higher resolution and precision inversion results, we propose an ICW method that gradually adds the gradient data, normalizes the inversion result as a weight function, and uses the cross-gradient to realize the structure constraint.First, the regularization inversion result and data constraint joint inversion result are normalized as the weight function of magnetic anomaly inversion and the combination matrix ([d ∆T ; d ∆T x ]) inversion, respectively.The objective functions are where weight function of magnetic anomaly inversion and the combination matrix ( ) inversion, respectively.The objective functions are ( ) ( ) is the number of data, M is the number of cells in the discretized model, and . 0.5 1.5   , the value for  should usually be near 1 [39]. 1 m and 2 m are the magnetic susceptibility results of regularization inversion and data constraint joint inversion, respectively, in (1).Here, the result with a better shallow inversion effect is applied as the weight to highlight the deep anomaly, and the result with a better deep inversion effect is selected as the weight to highlight the shallow anomaly.The cross use of this result weight can improve the resolution of field sources at different depths.
Second, the gradient anomaly in the y -direction is added in the inversion, ; ; The objective functions are ( ) ; ; The objective functions are , j = 1, . . ., M, weight function of magnetic anomaly inversion and the combination matrix ( ) inversion, respectively.The objective functions are ( ) ( ) is the number of data, M is the number of cells in the discretized model, and . 0.5 1.5   , the value for  should usually be near 1 [39]. 1 m and 2 m are the magnetic susceptibility results of regularization inversion and data constraint joint inversion, respectively, in (1).Here, the result with a better shallow inversion effect is applied as the weight to highlight the deep anomaly, and the result with a better deep inversion effect is selected as the weight to highlight the shallow anomaly.The cross use of this result weight can improve the resolution of field sources at different depths.
Second, the gradient anomaly in the y -direction is added in the inversion, ; ; The objective functions are ( ) ; ; The objective functions are , j = 1, . . ., M, N is the number of data, M is the number of cells in the discretized model, and 0.5 < β < 1.5, the value for β should usually be near 1 [39].m 1 and m 2 are the magnetic susceptibility results of regularization inversion and data constraint joint inversion, respectively, in (1).Here, the result with a better shallow inversion effect is applied as the weight to highlight the deep anomaly, and the result with a better deep inversion effect is selected as the weight to highlight the shallow anomaly.
The cross use of this result weight can improve the resolution of field sources at different depths.Second, the gradient anomaly in the y-direction is added in the inversion, m 1∆T 1 calculated in the previous step is taken as the weight of the magnetic anomaly inversion, and m 1∆T is taken as the weight of the combination matrix ( d ∆T ; d ∆T x ; d ∆T y ) inversion.The objective functions are (7) where m 2∆T and m 2∆T 2 are the magnetic susceptibility inversion results.
The gradient anomaly in the z-direction is added to the inversion, m 2∆T 2 calculated in the previous step, is used as the weight of the magnetic anomaly inversion, and m 2∆T is used as the weight of the combination matrix ( d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ) inversion.The objective functions are (8) where m 3∆T 3 is the final result by the ICW method, which can simultaneously achieve data and structure constraints, effectively balance the field source responses at different depths, and improve the resolution and precision of the inversion result.The specific implementation process is shown in Figure 4.
Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 23 ''' ''' ''' ''' min where 3 T m  is the final result by the ICW method, which can simultaneously achieve data and structure constraints, effectively balance the field source responses at different depths, and improve the resolution and precision of the inversion result.The specific implementation process is shown in Figure 4.

Theoretical Model Tests
The model shown in Figure 3 was used to verify the effect of the ICW method.To highlight the gradual improvement in the inversion performance at each stage of the method, we set the same truncated magnetic susceptibility.Figure 5a ) inversion of the second stage by the ICW method with a magnetic susceptibility larger than 0.075.Compared with the inversion results of the first stage, the resolution is further improved, and the value of magnetic susceptibility is somewhat closer to the true magnetic susceptibility; however, there is still a difference with the real location of the model.Figure 5g shows the magnetic susceptibility slice (y = 1900 m) of the final result by the ICW method.Figure 5h shows the volume-rendered image of the final result by the ICW method with a magnetic susceptibility larger than 0.075.Compared with the inversion results of the second stage, the value of magnetic susceptibility is closer to the true magnetic susceptibility, and the position and shape of the recovered magnetic susceptibility distribution are more consistent with the real model.Variable k represents the truncated magnetic susceptibility.The dotted prism identifies the real model.The burial depth of the two prisms included in the model in Figure 3 is simultaneously increased by 100 m, and the 3D distribution of the model and the corresponding anomalies are obtained, as shown in Figure 6a-d.The magnitude of the magnetic and gradient anomalies is further reduced, which will greatly affect the resolution and precision of the inversion.Figure 6e shows the volume-rendered images of the data constraint joint inversion of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ]) with a magnetic susceptibility larger than 0.028.When the truncated magnetic susceptibility is similar to the real magnetic susceptibility to the greatest extent, it can recover the deep geological body better but cannot restore the shallow geological body.Figure 6f,g show the volume-rendered images of the magnetic anomaly inversion and combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]) inversion by the cross-gradient method.Compared with Figure 6e, Figure 6f recovers the deep geological body better but still cannot recover the shallow geological body.In Figure 6g, when the truncated magnetic susceptibility is larger, more deep geological bodies can be recovered, and shallow geological bodies can be recovered to a certain extent, but the recovered model is quite different from the real model position (dotted box).Compared with Figure 6e-g, the ICW method recovers the model best, and the value of magnetic susceptibility is closer to the real magnetic susceptibility.However, in general, when the distance between the measurement and the field source is increased, the resolution of the inversion will inevitably be reduced, which will adversely affect the recovery of the model.Compared with other methods, the ICW method still has absolute advantages in this case.
Remote Sens. 2022, 14, x FOR PEER REVIEW 14 of 23 gradient anomalies is further reduced, which will greatly affect the resolution and precision of the inversion.Figure 6e shows the volume-rendered images of the data constraint joint inversion of the combination matrix ([ d d d d ) with a magnetic suscep- tibility larger than 0.028.When the truncated magnetic susceptibility is similar to the real magnetic susceptibility to the greatest extent, it can recover the deep geological body better but cannot restore the shallow geological body.Figure 6f,g show the volume-rendered images of the magnetic anomaly inversion and combination matrix ([ version by the cross-gradient method.Compared with Figure 6e, Figure 6f recovers the deep geological body better but still cannot recover the shallow geological body.In Figure 6g, when the truncated magnetic susceptibility is larger, more deep geological bodies can be recovered, and shallow geological bodies can be recovered to a certain extent, but the recovered model is quite different from the real model position (dotted box).Compared with Figure 6e-g, the ICW method recovers the model best, and the value of magnetic susceptibility is closer to the real magnetic susceptibility.However, in general, when the distance between the measurement and the field source is increased, the resolution of the inversion will inevitably be reduced, which will adversely affect the recovery of the model.Compared with other methods, the ICW method still has absolute advantages in this case.In actual situations, anomalous data will inevitably be disturbed by noise.We add Gaussian noise with a signal-to-noise ratio of 20 to the magnetic and gradient anomalies in Figure 3, as shown in Figure 7a-d  To further verify the sensitivity of the ICW method to noise intensity, we added Gaussian noise with a signal-to-noise ratio of 10 to the magnetic and gradient anomalies in Figure 3, as shown in Figure 8a-d.Using the same truncated magnetic susceptibility as the methods corresponding to Figure 7, the volume-rendered images of Figure 8 were obtained.We discovered that the magnetic susceptibility models recovered by the data constraint joint inversion method and cross-gradient method in Figure 8e-g were worse and that the precision was further reduced.Only the magnetic susceptibility model recovered by the ICW method in Figure 8h is almost the same as that without noise, which further verifies that the ICW method is extremely insensitive to noise.inversion by the data constraint joint inversion method with a magnetic susceptibility larger than 0.032, which is quite different from the real model.Figure 7f,g show the volume-rendered images obtained by cross-gradient inversion of the magnetic anomaly and combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]), respectively.Figure 7f compares Figure 7e,g with the truncated magnetic susceptibility closer to the real magnetic susceptibility, which restores a better magnetic susceptibility model for deep and shallow geological bodies.Although Figure 7g also restores the magnetic susceptibility model, it completely deviates from the real model position.Figure 7h shows the volume-rendered image of the final result by the ICW method with a magnetic susceptibility larger than 0.075.Compared with the inversion results of the previous methods, the value of magnetic susceptibility obtained by the ICW method is closer to the true magnetic susceptibility (κ = 0.25), and the magnetic susceptibility model recovered when the truncated magnetic susceptibility was most similar to the real model.Compared with the inversion results in Figure 5, the value of magnetic susceptibility is higher and the recovered magnetic susceptibility model is better, indicating that the ICW method has strong noise immunity.
To further verify the sensitivity of the ICW method to noise intensity, we added Gaussian noise with a signal-to-noise ratio of 10 to the magnetic and gradient anomalies in Figure 3, as shown in Figure 8a-d.Using the same truncated magnetic susceptibility as the methods corresponding to Figure 7, the volume-rendered images of Figure 8 were obtained.We discovered that the magnetic susceptibility models recovered by the data constraint joint inversion method and cross-gradient method in Figure 8e-g were worse and that the precision was further reduced.Only the magnetic susceptibility model recovered by the ICW method in Figure 8h is almost the same as that without noise, which further verifies that the ICW method is extremely insensitive to noise.Considering the complexity of the actual underground situation, we designed a model that contains four anomalous bodies with different burial depths to test the applicability of the ICW method.The top and bottom burial depths of prisms 1 and 3 and Considering the complexity of the actual underground situation, we designed a model that contains four anomalous bodies with different burial depths to test the applicability of the ICW method.The top and bottom burial depths of prisms 1 and 3 and prisms 2 and 4 were 200 m and 400 m and 300 m and 700 m, respectively.The magnetic susceptibility was 0.25, and the inclination and deflection were 45  Figure 10a shows the volume-rendered image of the combination matrix inversion by the data constraint joint inversion method with magnetic susceptibility value larger than 0.078.The data constraint joint inversion has better recovery for shallow geological bodies.Nevertheless, it is quite different from the actual situation.Figure 10b shows the volume-rendered image of the magnetic anomaly inversion by the cross-gradient method.Compared with Figure 10a, the magnetic susceptibility distributions of the deep and shallow parts are more accurate but are more different from the true magnetic susceptibility value.Figure 10c shows the volume-rendered images of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]) inversion by the cross-gradient method with a magnetic susceptibility larger than 0.1 and a low vertical resolution of the magnetic susceptibility distribution.Figure 10d shows the volume-rendered images of the final result by the ICW method with a value of magnetic susceptibility larger than 0.1.The resolution and precision of the inversion results are higher, the recovered magnetic susceptibility distribution is concentrated, there is a distinct boundary of physical property, and the value of magnetic susceptibility is more similar to the true magnetic susceptibility.

Real Data Application
Zhurihe is in the central part of the Inner Mongolia Autonomous Region, China, which is located on the northern margin of the North China Plate and spans two major tectonic units: the North China Craton and Central Asian Orogenic Belt.Bordered by the Bayan Obo-Chifeng fault, with the North China Craton in the south and the Xing-Meng orogenic belt in the north, the Zhurihe area has undergone complicated accretionary orogenic processes and multistage tectonic deformation events, and its tectonic evolution is very complicated [41][42][43].There are different rock types and genetic types, including Triassic monzogranite, Permian granite quartz diorite, Carboniferous granite, and Permian granite [44].Figure 11a shows the geological map of the survey area and its surrounding areas.
To clarify the position and range of regional metal mineral resources and aeromagnetic and gradient surveys were carried out outside the current mining area from April to July 2019.The distances between the survey lines and the survey points were 100 m and 40 m, respectively.The average flight height was 100 m, and the measurement error of the aeromagnetic field and its gradient anomaly was 1.91 nT and 0.02 nT/m, respectively.Aeromagnetic anomalies and gradient anomalies in the x, y, and z-directions are shown in Figure 11b-e   According to World Magnetic Model (WMM) 2015, the inclination and declination are 59 • and −3 • , respectively.The ICW method is used to carry out the joint inversion of aeromagnetic and gradient data.Slices were cut to analyze the characteristics of ore body trends.In Figure 12a, L1-L4 are the 4 slices extracted from the 3D inversion results by the data constraint joint inversion of the combination matrix ( d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ).The value of magnetic susceptibility recovered from the inversion result is low, making it difficult to delineate high magnetic bodies.In Figure 12b, L1-L4 are the 4 slices extracted from the 3D inversion results by the ICW method, and the dashed circle represents the magnetic susceptibility distribution with a magnetic susceptibility larger than 0.15.In Figure 12b, high magnetic body I inclines toward the south, especially in L2, with a larger inclination angle and wider extension.Therefore, we believe that I is the mining area with the most mining potential.Figure 12c shows the 3D inversion result by the ICW method with a magnetic susceptibility larger than 0.15.The inversion results indicate that the survey area mainly contains one large high magnetic body and four small high magnetic bodies.The volume of high magnetic body I is the largest, distributed in the northeast direction, and the top burial depth and bottom burial depth are 100 m and 960 m, respectively.Figure 12c shows that the largest high-value area is affected by the fault, which corresponds to high magnetic body I, so we infer that high magnetic body I is a large intrusive deposit.The depth of high magnetic bodies II and III is 180 m, the depth of high magnetic body IV is 200 m, and the depth of high magnetic body V is 100 m.Compared with magnetic body I, these high-magnetic bodies are smaller in volume, shallower in depth, and have minimal mining potential.
are 59° and −3°, respectively.The ICW method is used to carry out the joint inversion of aeromagnetic and gradient data.Slices were cut to analyze the characteristics of ore body trends.In Figure 12a, L1-L4 are the 4 slices extracted from the 3D inversion results by the data constraint joint inversion of the combination matrix (

;
; ; d d d d ).The value of magnetic susceptibility recovered from the inversion result is low, making it difficult to delineate high magnetic bodies.In Figure 12b, L1-L4 are the 4 slices extracted from the 3D inversion results by the ICW method, and the dashed circle represents the magnetic susceptibility distribution with a magnetic susceptibility larger than 0.15.In Figure 12b, high magnetic body I inclines toward the south, especially in L2, with a larger inclination angle and wider extension.Therefore, we believe that I is the mining area with the most mining potential.Figure 12c shows the 3D inversion result by the ICW method with a magnetic susceptibility larger than 0.15.The inversion results indicate that the survey area mainly contains one large high magnetic body and four small high magnetic bodies.The volume of high magnetic body I is the largest, distributed in the northeast direction, and the top burial depth and bottom burial depth are 100 m and 960 m, respectively.Figure 12c shows that the largest high-value area is affected by the fault, which corresponds to high magnetic body I, so we infer that high magnetic body I is a large intrusive deposit.The depth of high magnetic bodies II and III is 180 m, the depth of high magnetic body IV is 200 m, and the depth of high magnetic body V is 100 m.Compared with magnetic body I, these high-magnetic bodies are smaller in volume, shallower in depth, and have minimal mining potential.
Figure 12d shows the horizontal distribution of high magnetic bodies.The depth range of certain areas of high magnetic body I has been given, which has important guiding information for the future exploration and development of deep concealed ore bodies.Figure 12d shows the horizontal distribution of high magnetic bodies.The depth range of certain areas of high magnetic body I has been given, which has important guiding information for the future exploration and development of deep concealed ore bodies.

Figure 1 .
Figure 1.DRP, eigenvalue spectra for kernel function matrices and magnetic susceptibility slices ( 1900m y = ) obtained by magnetic anomaly ( T  d ) inversion and combination matrix (

Figure 1 .
Figure 1.DRP, eigenvalue spectra for kernel function matrices and magnetic susceptibility slices (y = 1900 m) obtained by magnetic anomaly (d ∆T ) inversion and combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ) inversion.(a) 3D view of 2 prisms and magnetic anomaly.(b) Magnetic gradient anomaly in the x-direction.(c) Magnetic gradient anomaly in the y-direction.(d) Magnetic gradient anomaly in the z-direction.(e,f) DRP of the magnetic anomaly and combination matrix of magnetic gradient anomalies.(g) Eigenvalue spectra for kernel function matrices.(h) Magnetic anomaly (d ∆T ) inversion by the regularization inversion method.(i) Combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ) inversion by the data constraint joint inversion method.(j) Combination matrix ( d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ) inversion by the data constraint joint inversion method.

Figure 2 .
Figure 2. DRP, eigenvalue spectra for kernel function matrices and magnetic susceptibility slices ( 1900m y = ) obtained by magnetic anomaly ( T  d ) inversion and combination matrix (

Figure 2 .
Figure 2. DRP, eigenvalue spectra for kernel function matrices and magnetic susceptibility slices (y = 1900 m) obtained by magnetic anomaly (d ∆T ) inversion and combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ) inversion.(a) 3D view of 2 prisms and magnetic anomaly.(b) Magnetic gradient anomaly in the x-direction.(c) Magnetic gradient anomaly in the y-direction.(d) Magnetic gradient anomaly in the z-direction.(e,f) DRP of the magnetic anomaly and combination matrix of magnetic gradient anomalies.(g) Eigenvalue spectra for kernel function matrices.(h) Magnetic anomaly (d ∆T ) inversion by the regularization inversion method.(i) Combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ) inversion by the data constraint joint inversion method.(j) Combination matrix ( d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ) inversion by the data constraint joint inversion method.

23 Figure 3 .
Figure 3. Magnetic and gradient anomalies and magnetic susceptibility slices ( 1900m y = ) obtained by different inversion methods.(a) 3D view of 2 prisms and magnetic anomaly.(b) Magnetic gradient anomaly in the x -direction. (c) Magnetic gradient anomaly in the y -direction. (d) Mag- netic gradient anomaly in the z -direction. (e) Magnetic susceptibility slice ( 1900m y = ) ob- 1900m y = ) by data constraint joint inversion of the combination matrix ( h) Magnetic susceptibility slice ( 1900m y = ) by cross-gradient inversion of the magnetic anomaly ( T  d ).(i) Magnetic susceptibility slice ( inversion of the combination matrix (

Figure 3 .
Figure 3. Magnetic and gradient anomalies and magnetic susceptibility slices (y = 1900 m) obtained by different inversion methods.(a) 3D view of 2 prisms and magnetic anomaly.(b) Magnetic gradient anomaly in the x-direction.(c) Magnetic gradient anomaly in the y-direction.(d) Magnetic gradient anomaly in the z-direction.(e) Magnetic susceptibility slice (y = 1900 m) obtained by regularization inversion of the magnetic anomaly (d ∆T ).(f) Magnetic susceptibility slice (y = 1900 m) by data constraint joint inversion of the combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ).(g) Magnetic susceptibility slice (y = 1900 m) by data constraint joint inversion of the combination matrix ( d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ).(h) Magnetic susceptibility slice (y = 1900 m) by cross-gradient inversion of the magnetic anomaly (d ∆T ). (i) Magnetic susceptibility slice (y = 1900 m) by cross-gradient inversion of the combination matrix ( d ∆T x ; d ∆T y ; d ∆T z ).
previous step is taken as the weight of the magnetic anomaly inversion, and 1 T  m is taken as the weight of the combination matrix (

m
are the magnetic susceptibility inversion results.The gradient anomaly in the z -direction is added to the inversion, previous step, is used as the weight of the magnetic anomaly inversion, and 2 T  m is used as the weight of the combination matrix ( ; previous step is taken as the weight of the magnetic anomaly inversion, the weight of the combination matrix (

m
are the magnetic susceptibility inversion results.The gradient anomaly in the z -direction is added to the inversion, previous step, is used as the weight of the magnetic anomaly inversion, and 2 T  m is used as the weight of the combination matrix ( ;

Figure 4 .
Figure 4. Flowchart of the ICW method.

Figure 4 .
Figure 4. Flowchart of the ICW method.
,b show the magnetic susceptibility slices (y = 1900 m) of the first stage of magnetic anomaly inversion and combination matrix ([d ∆T ; d ∆T x ]) inversion by the ICW method.Figure 5c shows the volume-rendered image of the combination matrix ([d ∆T ; d ∆T x ]) inversion of the first stage by the ICW method with a magnetic susceptibility larger than 0.075.The resolution and precision of the inversion results in Figure 5a,b are significantly better than those of the data constraint joint Inversion and the cross-gradient inversion, but the inversion results in Figure 5c deviate greatly from the real model.Figure 5d,e show the magnetic susceptibility slices (y = 1900 m) of the second stage of magnetic anomaly inversion and combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]) inversion by the ICW method.Figure 5f shows the volume-rendered image of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]

Figure 5 .
Figure 5. Magnetic susceptibility slice (y = 1900 m) and volume-rendered images by the ICW method.(a) Magnetic susceptibility slice (y = 1900 m) of the first stage of magnetic anomaly (d ∆T ) inversion by the ICW method.(b) Magnetic susceptibility slice (y = 1900 m) of the first stage of combination matrix ([d ∆T ; d ∆T x ]) inversion by the ICW method.(c) volume-rendered image of the combination matrix ([d ∆T ; d ∆T x ]) inversion of the first stage by the ICW method.(d) Magnetic susceptibility slice (y = 1900 m) of the second stage of magnetic anomaly (d ∆T ) inversion by the ICW method.(e) Magnetic susceptibility slice (y = 1900 m) of the second stage of combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]) inversion by the ICW method.(f) volume-rendered image of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]) inversion of the second stage by the ICW method.(g) Magnetic susceptibility slice (y = 1900 m) of the final result by the ICW method.(h) volume-rendered image of the final result by the ICW method.

Figure 6 .
Figure 6.Magnetic and gradient anomalies (compared with the model in Figure 2, the buried depths of the prisms are increased by 100 m) and the volume-rendered images by different methods.(a) 3D view of 2 prisms and magnetic anomaly.(b) Magnetic gradient anomaly in the x -direction. (c)

Figure 6 .
Figure 6.Magnetic and gradient anomalies (compared with the model in Figure 2, the buried depths of the prisms are increased by 100 m) and the volume-rendered images by different methods.(a) 3D view of 2 prisms and magnetic anomaly.(b) Magnetic gradient anomaly in the x-direction.(c) Magnetic gradient anomaly in the y-direction.(d) Magnetic gradient anomaly in the z-direction.(e) volume-rendered

Figure 7 .
Figure 7. Magnetic and gradient anomalies with Gaussian noise and volume-rendered images by different inversion methods.(a-d) Noise-corrupted magnetic and gradient anomalies with Gaussian noise of a signal-to-noise ratio of 20.(e) Volume-rendered image of combination matrix ( [ ; ; ; ] x y z T T T T     d d d d ) inversion by the data constraint joint inversion method.(f) Volume-ren- dered image of magnetic anomaly ( T  d ) inversion by the cross-gradient method.(g) Volume-ren- dered image of the combination matrix ( [ ; ; ] y T T x T    d d d

Figure 7 .
Figure 7. Magnetic and gradient anomalies with Gaussian noise and volume-rendered images by different inversion methods.(a-d) Noise-corrupted magnetic and gradient anomalies with Gaussian noise of a signal-to-noise ratio of 20.(e) volume-rendered image of combination matrix ([d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ]) inversion by the data constraint joint inversion method.(f) volume-rendered image of magnetic anomaly (d ∆T ) inversion by the cross-gradient method.(g) volume-rendered image of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]) inversion by the cross-gradient method.(h) volume-rendered image of the final result by the ICW method.

Figure
Figure 7e shows the volume-rendered image of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ])inversion by the data constraint joint inversion method with a magnetic susceptibility larger than 0.032, which is quite different from the real model.Figure7f,g show the volume-rendered images obtained by cross-gradient inversion of the magnetic anomaly and combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]), respectively.Figure7fcompares Figure7e,g with the truncated magnetic susceptibility closer to the real magnetic susceptibility, which restores a better magnetic susceptibility model for deep and shallow geological bodies.Although Figure7galso restores the magnetic susceptibility model, it Remote Sens. 2022, 14, x FOR PEER REVIEW 17 of 23

Figure 8 .)
Figure 8. Magnetic and gradient anomalies with Gaussian noise and volume-rendered images by different inversion methods.(a-d) Noise-corrupted magnetic and gradient anomalies with Gaussian noise of a signal-to-noise ratio of 10.(e) Volume-rendered image of combination matrix ( [ ; ; ; ] x y z T T T T     d d d d ) inversion by the data constraint joint inversion method.(f) Volume-ren- dered image of magnetic anomaly ( T  d ) inversion by the cross-gradient method.(g) Volume-ren- dered image of the combination matrix ( [ ; ; ] y T T x T    d d d

Figure 8 .
Figure 8. Magnetic and gradient anomalies with Gaussian noise and volume-rendered images by different inversion methods.(a-d) Noise-corrupted magnetic and gradient anomalies with Gaussian noise of a signal-to-noise ratio of 10.(e) volume-rendered image of combination matrix ([d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ]) inversion by the data constraint joint inversion method.(f) volume-rendered image of magnetic anomaly (d ∆T ) inversion by the cross-gradient method.(g) volume-rendered image of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]) inversion by the cross-gradient method.(h) volume-rendered image of the final result by the ICW method.

1 Figure 9 .
Figure 9. 3D view of the 4 prisms and corresponding magnetic and gradient anomalies.(a) 3D view of the 4 models and the magnetic anomaly.(b) Magnetic gradient anomaly in the x-direction.(c) Magnetic gradient anomaly in the y-direction.(d) Magnetic gradient anomaly in the z-direction.

Figure
Figure 9b-d show the magnetic gradient anomalies in the x-direction, y-direction, and z-direction with a flight height of 100 m.The gradient anomalies can highlight shallow geological bodies.Figure10ashows the volume-rendered image of the combination matrix inversion by the data constraint joint inversion method with magnetic susceptibility value larger than 0.078.The data constraint joint inversion has better recovery for shallow geological bodies.Nevertheless, it is quite different from the actual situation.Figure10bshows the volume-rendered image of the magnetic anomaly inversion by the cross-gradient method.Compared with Figure10a, the magnetic susceptibility distributions of the deep and shallow parts are more accurate but are more different from the true magnetic susceptibility value.Figure10cshows the volume-rendered images of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]) inversion by the cross-gradient method with a magnetic susceptibility larger than 0.1 and a low vertical resolution of the magnetic susceptibility distribution.Figure10dshows the volume-rendered images of the final result by the ICW method with a value of magnetic susceptibility larger than 0.1.The resolution and precision of the inversion results are higher, the recovered magnetic susceptibility distribution is concentrated, there is a distinct boundary of physical property, and the value of magnetic susceptibility is more similar to the true magnetic susceptibility.

Figure 10 .
Figure 10.volume-rendered images by different inversion methods.(a) volume-rendered image of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ]) inversion by the data constraint joint inversion method.(b) volume-rendered image of magnetic anomaly (d ∆T ) inversion by the cross-gradient method.(c) volume-rendered image of combination matrix ([d ∆T ; d ∆T x ; d ∆T y ]) inversion by the cross-gradient method.(d) volume-rendered image of the final result by the ICW method.
. The aeromagnetic and gradient anomalies are mainly distributed in the northeast direction, and the direction of the survey line is distributed in the north-south direction.The scope of the survey area is 2250 × 1470 m. 40 m, respectively.The average flight height was 100 m, and the measurement error of the aeromagnetic field and its gradient anomaly was 1.91 nT and 0.02 nT/m, respectively.Aeromagnetic anomalies and gradient anomalies in the x , y , and z -directions are shown in Figure 11b-e.The aeromagnetic and gradient anomalies are mainly distributed in the northeast direction, and the direction of the survey line is distributed in the northsouth direction.The scope of the survey area is 2250 × 1470 m.

Figure 11 .
Figure 11.Geographical location and geological conditions and the measured aeromagnetic and gradient anomalies.(a) Geographical location and geological conditions of the survey area.(b) Measured aeromagnetic anomaly.(c) Measured aeromagnetic gradient anomaly in the x -direc-

Figure 11 .
Figure 11.Geographical location and geological conditions and the measured aeromagnetic and gradient anomalies.(a) Geographical location and geological conditions of the survey area.(b) Measured aeromagnetic anomaly.(c) Measured aeromagnetic gradient anomaly in the x-direction.(d) Measured aeromagnetic gradient anomaly in the y-direction.(e) Measured aeromagnetic gradient anomaly in the z-direction.

Figure 12 .
Figure 12.Inversion result of the measured data.(a) Magnetic susceptibility slice by the data constraint joint inversion of the combination matrix ( [ ; ; ; ] x y z T T T T     d d d d ).(b) Magnetic suscepti- bility slice by the ICW method.(c) 3D inversion result by the ICW method.(d) Horizontal range distribution of the high magnetic body.

Figure 12 .
Figure 12.Inversion result of the measured data.(a) Magnetic susceptibility slice by the data constraint joint inversion of the combination matrix ([d ∆T ; d ∆T x ; d ∆T y ; d ∆T z ]).(b) Magnetic susceptibility slice by the ICW method.(c) 3D inversion result by the ICW method.(d) Horizontal range distribution of the high magnetic body. .