High-Resolution Cooperate Density-Integrated Inversion Method of Airborne Gravity and Its Gradient Data

Airborne (or satellite) gravity measurement is a commonly used remote sensing method to obtain the underground density distribution. Airborne gravity gradiometry data have a higher horizontal resolution to shallower causative sources than airborne gravity anomaly, so joint exploration of airborne gravity and its gradient data can simultaneously obtain the anomaly feature of sources with different depths. The most commonly used joint inversion method of gravity and its gradient data is the data combined method, which is to combine all the components into a data matrix as mutual constraints to reduce ambiguity and non-uniqueness. In order to obtain higher resolution results, we proposed a cooperate density-integrated inversion method of airborne gravity and its gradient data, which firstly carried out the joint inversion using cross-gradient constraints to obtain two density structures, and then fused two recovered models into a result through Fourier transform; finally, data combined joint inversion of airborne gravity, and gradient data were reperformed to achieve high-resolution density result using fused density results as a reference model. Compared to the data combined joint inversion method, the proposed cooperate density-integrated inversion method can obtain higher resolution and more accurate density distribution of shallow and deep bodies meanwhile. We also applied it to real data in the mining area of western Liaoning Province, China. The results showed that the depth of the skarn-type iron mine in the region is about 900–1300 m and gives a more specific distribution compared to the geological results, which provided reliable data for the next exploration plan.


Introduction
Airborne gravity measurement is a usual used remote sensing method to obtain the underground density structure. However, it is often ill-posed and requires certain prior information and constraints to guarantee the results that are unique and stable. Many researchers made numerous contributions based on the Tikhonov regularization method for solving this problem [1][2][3][4][5][6][7][8][9]. Compared to gravity exploration, the principal advantage of gravity gradiometry has a higher resolution in the near surface, which helps to obtain a better-characterized target body and more accurate shape and orientation delineation of geological bodies [10]. Chen obtained the gravity gradient tensors by building a matrix equation of the gravity vector and its neighbors by Taylor series expansion. Compared to the real gravity gradient tensors data from the Otway Basin in Australia, the results obtained by the proposed approach had a relative error [11]. Dai proposed a new approach of using the component of gravity change and GRACE data to reveal the high spatial resolution and confirms that the GRACE gravity and gravity gradient changes agree well with seismic model spectra [12]. Gravity gradiometer can measure the small change of gravity at two points, which contains more abundant navigation and positioning information than gravity. Gao proposed an aided navigation method based on strapdown gravity gradiometer. The performance of aided navigation is analyzed and evaluated from six aspects, which contained the high-resolution results [13]. Cattin combined the ground gravity data with satellite gradients in a joint inversion to assess the location and the geometry of transition and obtained a ca. 10 km-wide transition zone located at the western border of Bhutan [14]. Plasman inverted the six tensor components of GOCE gravity gradient data and concluded that the simultaneous inversion of several components displayed a significant improvement for the global tensor recovery. The proposed method was successfully verified in complex subduction cases with gradient and gravity data [15]. Residual errors in the terrain correction could lead to errors in data interpretation. If a desired terrain correction error could be given, we could select an optimum survey flying height over a known terrain [16]. Kass provided the basis for optimizing the terrain correction to improve efficiency is not only gravity gradiometry terrain corrections and forward modeling [17]. So, the joint exploration of gravity and its gradient data has further expanded the application range [18][19][20][21][22][23].
The joint inversion method of gravity and its gradient data has two methods that are based on data combined and structure constrained. The structure constrained joint inversion method is mostly based on cross-gradient function [24][25][26][27]. Moorkamp presented a 3-D joint inversion framework for seismic, magnetotelluric (MT), and scalar and tensorial gravity data. Using large-scale optimization methods, parallel forward solvers to investigate different coupling approaches for the various physical parameters [28]. Gallardo summarized the role played by the structural gradients-based approach for coupling fundamentally different physical fields in geophysical inversion [29]. Zhang extended the data-space joint inversion algorithm of magnetotelluric, gravity, and magnetic data to include first-arrival seismic travel-time and normalized cross-gradient constraints. This method could effectively improve the computational speed and greatly reduce memory requirements [30]. Jiang developed a 3D joint seismic waveform and gravity gradiometry inversion method. A case study in the South China Sea showed that joint inversion results are consistent with the pre-stack depth migration section, and the shape of the salt body is well resolved [31].
However, this method presents two inversion results that cannot reflect the data correlation. So, the most joint inversion method of gravity and its gradient data are based on data combined. Zhdanov combined gravity curvature, gravity data, and gravity gradient data for the joint inversion [5]. Wu proposed an adaptive weighting method to achieve joint inversion of gravity and gravity gradient data [32]. Wan and Zhdanov developed the iterative migration of gravity and gravity gradient data, which can obtain a more accurate density distribution [33]. Capriotti and Li used an adaptive weighting function and provided a general equation for the joint inversion of gravity and gravity gradient data [34]. Zhdanov and Lin proposed an adaptive polynomial inversion based on a regularized conjugate gradient inversion method and achieved suitable results [35]. Liu proposed a joint density inversion for gravity gradient data based on the Tikhonov regularization method [36]. The structure constrained joint inversion method is mostly based on cross-gradient function, which is an effective method for the joint inversion of different geophysical data.
Since the cross-gradient function is mostly used in the joint inversion method between different physical parameters, we introduce structural constraints for different components of the gravity data and obtain two density inversion results. Data fusion can combine the respective advantages of different geophysical data to obtain high-resolution results. Therefore, we apply data fusion to the above two density inversion results and combine their respective advantages to obtain the final result. Data fusion is realized by decomposing data into different frequency information and combine them to realize data recovery [37][38][39][40][41]. People take fused data into the inversion process to ensure that the information is complete [42].
We presented a cooperate density-integrated inversion method, which used the crossgradient and fusion methods to improve the resolution of the joint inversion of airborne gravity and its gradient data. The proposed method is also applicable to joint inversion of ground gravity and airborne gravity gradients. The high-resolution method is verified on synthetic and real data.

High-Resolution Cooperate Density-Integrated Inversion Method
The specific inversion objective function of the data combined joint inversion of airborne gravity and its gradient data computed by Tikhonov regularized density inversion method is given by [43] where G z represents the kernel function matrix of airborne gravity data, G xx , G xy , G xz , G yz, and G zz represent the kernel function matrix of various airborne gravity gradient data. dz stands for gravity anomaly, d xx , d xy , d xz , d yz, and d zz represent gravity gradient components. m is a vector of density parameter. The magnitude of the regularization coefficient is crucial in defining the model smoothness and the data misfit. Previous studies have demonstrated that it is efficient if the inversion begins with a large regularization parameter [44,45]. We select regularization parameter α as 10 n (n is an integer), which is according to experience. At each iteration k, the parameter α is reduced slowly dependent on parameter γ, using α (k) = α (k−1) γ, 0 ≤ γ < 1. m 0 is a reference model, if available. W m = 1 (z+z 0 ) β is depth weighting designed to counteract the decay of sensitivities [3,46], and the parameter of z 0 can be obtained by matching the depth weighting function with the kernel function beneath the observation point, and β is the empirical constant. The airborne gradient data has a higher resolution to shallow sources because the density inversion is sensitive to the data with a higher rate. The joint inversion result has a higher resolution to shallow sources. The specific calculating process can be referenced as follow [47,48], Set G W = GW −1 m and m W = W m m, and the object function can be summarized as The derivative of the above objective function is After the optimal solution m W is obtained, m can be calculated by In order to improve the resolution to deep sources, we present a high-resolution cooperate density-integrated inversion method of integrated airborne gravity and its gradient data, and the specific computation processes are shown in Figure 1. Firstly, we use the cross-gradient function to accomplish the joint inversion of airborne gravity and its gradient data, and this method can effectively keep the characteristics of different types of data. The cross-gradient function is mostly used for joint inversion between different physical parameters, and the cross-product of the gradient vectors of different physical parameters is used to reflect the structural similarity of the target body. The cross-gradient term is introduced into the inversion objective function to achieve structural constraints. So, the inversion result computed by the cross-gradient method is more convergent. The cross-gradient function of the two physical parameters m (1) and m (2) are defined as follows [49,50] Φ cross (m (1) , m (2) t y (m (1) , m (2) ) = ∂m (2) where m (1) and m (2) represent the gradient of two physical parameters. Although the airborne gravity and its gradient data both reflect the density change, the characteristics of the two data are different. The cross-gradient constraint function is considered for structural constraints to achieve joint inversion of different forms of data under the same physical property. The objective functions of airborne gravity and its vertical gradient cross-gradient inversion are as follows where G (1) and G (2) represent the kernel function matrices corresponding to airborne gravity data and its gradient data, respectively. m (1) and m (2) represent the density parameters corresponding to airborne gravity data and its gradient data, respectively. d (1) and d (2) represent gravity anomaly and its gradient anomaly. W d is the data-weighting matrix. α is the regularization coefficient that is optimized at each iteration to minimize the error weighted root mean square error. λ is the coefficient of the cross-gradient terms. The amount of structural similarity obtained through the joint inversion algorithm can be adjusted using different choices of λ [51]. At each iteration k, the parameter λ is reduced slowly dependent on parameter γ, using λ (k) = λ (k−1) γ, where 0 ≤ γ < 1. The cross-gradient method can effectively reflect the features that the original airborne gravity data are sensitive to deep sources, and the airborne gravity gradient data can highlight the shallow sources. Secondly, we merge two different density results computed by the cross-gradient method into a result under the condition of keeping the characteristic of different results, and the Fourier and wavelet transforms are two commonly used fusion methods [52,53]. The different wavelet decomposition orders will affect the fusion results, and it is hard to ascertain the reasonable order for real data, so we adopt the Fourier transform fusion method, and the expression is where F represents Fourier transform, F −1 represents inverse Fourier transform. F g (ω) and F g zz (ω) are the airborne gravity data, and its gradient data in the frequency domain after Fourier transform, and H lowpass represents the low pass filter, which remains the low-frequency information of gravity data. Finally, we take the fusion results m fusion as the reference model and use the format of the Tikhonov regularized density inversion method to obtain the final density result by the combination of airborne gravity and its gradient data, and the computation process solved by the conjugate gradient algorithm is k is the number of iteration, and k = 0 means the initial condition, and the iteration termination condition is that data fitting error is less than the given threshold.

Theoretical Model Tests
Airborne gravity gradiometry data has a higher horizontal resolution to shallower causative sources than airborne gravity anomaly, so joint exploration of airborne gravity and its gradient data can simultaneously obtain the anomaly feature of sources with different depths. In order to verify the effectiveness of the high-resolution cooperate density-integrated inversion method, we build models of two prisms in different depths with the density contrast of 1000 kg/m 3, which have depths from 2 to 4 km and 4 to 7 km, and sizes of models are 2 × 2 × 2 km and 4 × 2 × 3 km, as shown in Figure 2a. Compared with Figure 2b,c, airborne gravity gradiometry data better highlight shallow information than airborne gravity anomaly, and airborne gravity data better displays deep information.
To verify the high-resolution results of the proposed cooperate density-integrated inversion method, we come up with the density results by common Tikhonov regularized method, data combined joint inversion method, structure constrained joint inversion method, and cooperate density-integrated inversion method. The subsurface is divided into 30 × 20 × 10 cubic prisms with an edge length of 1 km. Figure 3a shows the density inversion result of gravity anomaly by common Tikhonov regularized method with y-slice of 10 km, and the results have a lower resolution to the shallower body. The root mean square misfit of data is 0.15. Figure 3b represents the density result of airborne gravity gradient data by Tikhonov regularized inversion method with y-slice of 10 km, which can clearly show the location of the shallow source. The root mean square misfit of data is 0.01. Figure 3c,d shows the density slice and 3D density distribution by data combined joint inversion method of airborne gravity and its gradient data. The root mean square misfit of data is 0.008. We find that the joint inversion improves the resolution and accuracy of the shallow source, but it has a lower resolution to the deep source. Figure 3e,f show the density inversion results of airborne gravity and gradient data computed by structure constrained joint inversion method, respectively, and it can recover a more compact model than separate inversion. The root mean square misfits of data are 0.004 and 0.002. The cooperate density-integrated inversion method is introduced and combined with the advantages of the results by different data, the vertical slice and 3D perspective view with a value larger than 350 kg/m 3 are shown in Figure 3g,h, and we obviously find that the shape of recovered density models is refined. The root mean square misfit of data is 0.001. It is proved that the cooperate density-integrated inversion method of airborne gravity and its gradient data have a higher horizontal and vertical resolution that is more convergent compared to the data combined joint inversion method.
During the density inversion, the existence of noise will have a great impact on the inversion results, and noise is inevitable in actual situations. To test the anti-noisy of the high-resolution cooperate density-integrated inversion method and simulate the actual noise condition, we build the same model contaminated by 5% Gaussian noise in Figure  2a. Figure 4b displays noise-corrupted airborne gravity anomaly and its gradient anomaly, and the noise has an obvious impact on the anomaly shape.   The airborne gravity and its gradient data are affected by noise. The density result with a slice of y = 10 km computed by the data combined joint inversion method is distorted, as shown in Figure 5a. So, the 3D density distribution with larger than 350 kg/m 3 cannot obtain the results with high resolution. Figure 5c,d show the density results with a slice of y = 10 km and the 3D density distribution with larger than 350 kg/m 3 computed by the presented cooperative inversion method. Although the airborne gravity and its vertical gradient data are affected by noise, the density results obtained by the presented cooperative inversion method can obtain better recovery of physical parameters and more accurate position descriptions. For testing the universality of the proposed method, we build models with two prisms buried at the same depth, and the density contrast of the prism is 1000 kg/m 3, which have depths from 2 to 4 km, and the size of models are 4 × 2 × 2 km and 4 × 2 × 2 km, as shown in Figure 6a. The airborne gravity anomaly and gravity gradient anomaly are as shown in Figure 6b,c. Previous studies have shown that a gradiometer is better than a gravimeter in detecting short-wavelength anomalies at the same depth [54]. For the models buried at the same depth, airborne gravity gradiometry data achieve a more detailed description and a better description of the shape of target bodies compared to the airborne gravity data. The subsurface is divided into 20 × 20 × 10 cubic prisms with an edge length of 1 km.
We obtain the density results with a slice of y = 10 km and the 3D density distribution with larger than 370 kg/m 3 computed by the data combined joint inversion method as shown in Figure 7a,b. The density results are not focused enough and cannot describe the distribution of models clearly. The density results with a slice of y = 10 km and the 3D density distribution larger than 370 kg/m 3 computed by the presented cooperate density-integrated inversion method are shown in Figure 7c,d. Compared to the density results computed by the data combined joint inversion method, the results computed by the presented cooperative inversion method are more convergent, and both horizontal and vertical resolution is higher, and the density distributions are closer to the true value. So, the presented cooperative inversion method can obtain high-resolution results for the models buried at the same depth. In actual situations, the distributions of iron mines are complex. To verify the acceptability of the proposed method under the actual situation, we build complex models with three prisms, and the density contrasts of the prisms are 800, 1000, and 1200 kg/m 3, which have depths from 2 to 4 km, 4 to 7 km, and 2 to 4 km, as shown in Figure 8a, and the sizes of models are 4 × 4 × 2 km, 6 × 2 × 3 km, and 4 × 4 × 2 km. The airborne gravity anomaly and gravity gradient anomaly are as shown in Figure 8b,c. Also, the subsurface is divided into 40 × 20 × 10 cubic prisms with an edge length of 1 km.
The density results with a slice of x = 10 km and y = 10 km computed by the data combined joint inversion method as shown in Figure 9a,c, and the black lines represent the true distribution of sources. The results can obtain better recovery of physical parameters to shallower causative sources. However, the results cannot obtain a suitable description of deep causative sources. The density results with a slice of x = 10 km and y = 10 km computed by the presented cooperative inversion method are shown in Figure 9b,d. The 3D density distributions with larger than 370 kg/m 3 computed by the data combined joint inversion method and the presented cooperative inversion method are shown in Figure 9e,f. The results computed by the presented cooperative inversion method have a higher horizontal and vertical resolution, which is more convergent compared to the current inversion method, and the density distributions are closer to the true value. So, the presented cooperative inversion method can obtain high-resolution results for the complex model tests. (c) Density slice (y = 10 km) by data combined joint inversion method of airborne gravity and its gradient data. (d) Density slice (y = 10 km) by the cooperate density-integrated inversion method. (e) 3D density distribution (larger than 370 kg/m 3 ) by data combined joint inversion method of airborne gravity and its gradient data. (f) 3D density distribution (larger than 370 kg/m 3 ) by the cooperate density-integrated inversion method of airborne gravity and its gradient data.
To test the accuracy of the proposed cooperate density-integrated inversion method depends upon the airborne's height, we build the same complex model in Figure 10a. With the airborne's height increase, the maximum of gravity and its vertical gradient data reduce significantly, and the characteristics are very different from the form in Figure 10b,c. For density inversion results calculated at 150 m altitude, the recovery of the shallow bodies is missing, as shown in Figure 10c,d, and the deep body is out of convergence. Figure 10g,h are density slice (y = 10 km) by cooperate density-integrated inversion method at 200 m altitude and 3D density distribution (larger than 370 kg/m 3 ). For density inversion results calculated at 200 m altitude, it is far from the real model. So, with the increase in altitude, it has a great influence on the final inversion results, and the results are more deviated from the real model.
When the airborne gravity gradient data are not measured, as an alternative, the airborne gravity vertical gradient data can be calculated from airborne gravity data, and we calculate airborne gravity vertical gradient data by Fourier transform method. In order to verify the feasibility of this method, we carried out the following work. Firstly, we compare the calculated airborne gravity gradient data and theoretical airborne gravity gradient data ( Figure 8c). Then, we obtain the density slice by cooperate density-integrated inversion method of airborne gravity and its calculated vertical gradient data. Finally, we achieve the 3D density distribution (larger than 370 kg/m 3 ) by cooperate density-integrated inversion method.
We find that the theoretical airborne gravity gradient data ( Figure 8c) and calculated airborne gravity gradient data (Figure 11a) are almost the same. The density results ( Figure 11b) with a slice of y = 10 km computed by the proposed cooperate densityintegrated inversion method are close to the same results in Figure 9d. Three-dimensional density distributions (Figure 11c) with larger than 370 kg/m 3 computed by the presented cooperative inversion method are consistent with the same distribution in Figure 9f. Compared with the results, we confirm that the calculated vertical gradient data is also available when missing the real gravity gradient data, and the inversion results by cooperate densityintegrated inversion method are almost the same. So, we can use the calculated gravity gradient data to replace the real gravity gradient data under the condition of missing data.

Real Data Application
The western Liaoning Province is rich in iron ore bodies, which shows the high value of airborne gravity anomaly. In order to obtain the clear distribution of the mineral-induced anomalous source, we use the joint exploration of airborne gravity and its vertical gradient data with a sampling interval of 250 m, and the vertical gradient data are obtained by calculating the vertical derivative of gravity data. The geological map is shown in Figure 12a, and the distributions of the iron mines are predicted by the regional structural feature, and both fault and Quaternary coverage areas are identified by geological information [55,56]. The airborne gravity anomaly and its vertical gradient anomaly are shown in Figure 12b,c, and there is a suitable corresponding relationship between the high values of anomalies and the predicted locations of iron mine, and the gravity gradient data has a higher resolution. Based on the drilling data W1 as shown in Figure 12d, we can verify the accuracy of the inversion results, and the drilling data also provide the accurate density contrast for obtaining the delineation of iron ore bodies.  Through the demonstration in Figure 11, calculated gravity gradient data can replace the actual gravity gradient data, and the inversion result can obtain high-resolution and accurate results by the proposed cooperate density-integrated inversion method. Density inversion results are carried out for airborne gravity data, and its calculated vertical gradient data in the western Liaoning area by data combined joint inversion method and the proposed cooperate density-integrated inversion method. There has a drilled borehole in the survey area, and the borehole is located inside the high-value area. Based on the drilling data W1, the high-density magnetite ore bodies are almost 3120 kg/m 3 , and surrounding skarn are almost 2900 kg/m 3 . It is believed that high-density bodies that are greater than or equal to 200 kg/m 3 correspond well to the range of iron mines. So, we retain density bodies with a value greater than 200 kg/m 3 , as shown in Figure 13a,c.
In order to confirm the buried depth of the density results, we obtain the horizontal slices of 3D density results by data combined joint inversion method from 0 to 1500 m as shown in Figure 13b and obtain the horizontal slices of 3D density results by proposed cooperative density inversion method from 0 to 1500 m as shown in Figure 13d. We estimate that the high-density bodies were buried from 900 to 1300 m. According to the cooperative density inversion result, we delineate the horizontal range of high-density bodies, as shown in Figure 14. Compared to the inferred range of iron mines by geological information, the density inversion result by the proposed cooperative method has a better agreement than the density by data combined joint inversion method. Therefore, cooperate density-integrated inversion method has the ability to achieve higher resolution density inversion result by gravity data and its gradient data, and the proposed method presents more detailed features compared to the geological results. Figure 14b displays that the results computed by the proposed method are more accurate and have higher resolution and more specific detailed descriptions of the iron bodies. Finally, we identify five iron mines depending on the inversion results of the proposed method, named I-V. The buried depth ranges of the five iron bodies are determined according to the two profiles. The buried depth of the iron mine area I is 900-1300 m. The buried depth of the iron mine area II is 1000-1100 m, the iron mine area III is 900-1200 m, the iron mine area IV is 950-1200 m, and the iron mine area V is 1000-1200 m.

Conclusions
We propose a cooperate density-integrated inversion method of airborne gravity and its gradient data, which used the cross-gradient and fusion tools to obtain higher resolution results, especially that the proposed method can better obtain the distribution of the deep source. We tested the proposed method by synthetic model tests and real data. Both model tests verified that the cooperate density-integrated inversion method of airborne gravity and its gradient data can obtain high-resolution results, which are buried at the same depth and different depths. Moreover, inversion results of complex models showed the feasibility of the proposed method. This method could recover higher resolution density models compared to the data combined joint inversion. In practice, calculated airborne gravity vertical gradient data can replace the theoretical airborne gravity vertical gradient data when the airborne gravity gradient data are not measured.
The inversion results of real data provide the location of the mining area accurately and delineate five iron ore bodies buried from 900 to 1300 m. The cooperate density-integrated inversion method provides a better tool to obtain high-resolution density distribution of airborne gravity and its gradient data.