The Efficient 3D Gravity Focusing Density Inversion Based on Preconditioned JFNK Method under Undulating Terrain: A Case Study from Huayangchuan, Shaanxi Province, China

Since polymetallic ores show higher anomalies in gravity exploration methods, we usually obtain the position and range of ore bodies by density inversion of gravity data. The three-dimensional (3D) gravity focusing density inversion is a common interpretation method in mineral exploration, which can directly and quantitatively obtain the density distribution of subsurface targets. However, in actual cases, it is computation inefficient. We proposed the preconditioned Jacobian-free Newton-Krylov (JFNK) method to accomplish the focusing inversion. The JFNK method is an efficient algorithm in solving large sparse systems of nonlinear equations, and we further accelerate the inversion process by the preconditioned technique. In the actual area, the gravity anomalies are distributed on the naturally undulating surface. Nowadays, the gravity inversion under undulating terrain was mainly achieved by discretizing the ground into unstructured meshes, but it is complicated and time-consuming. To improve the practicality, we presented an equivalent-dimensional method that incorporates unstructured meshes with structured meshes in gravity inversion, and the horizontal size is determined by the gradient of observed gravity and terrain data. The small size meshes are adopted at the position where the terrain or gravity gradient is large. We used synthetic data with undulating-terrain to test our new method. The results indicated that the recovered model obtained by this method was similar to the inversion method of unstructured meshes, and the new method computes faster. We also applied the method to field data in Huayangchuan, Shaanxi Province. The survey area has complicated terrain conditions and contains multiple polymetallic ores. Based on the high-density characteristics of polymetallic ore bodies in the area, we calculate the field data into 3D density models of the subsurface by the preconditioned JFNK method and infer six polymetallic ores.


Introduction
The gravity exploration method is used to detect polymetallic ores because of its higher density feature, and they can use the density inversion method of gravity anomaly to obtain the approximate horizontal position of the ores. The most commonly gravity inversion method is three-dimensional gravity inversion, which can directly obtain the subsurface density distribution to interpret the

Methodology of Preconditioned JFNK Gravity Focusing Inversion
The gravity inversion involved discretizing the subsurface into a finite number of rectangular units, and the discrete map is shown in Figure 1. There are n observation points and m rectangular units, the forward expression of the gravity anomaly could be expressed as: where g n×1 is the gravity anomaly, G n×m is the contribution to the nth datum of a unit density in the mth cell called kernel function matrix, and ρ m×1 is the density matrix of rectangular units.

Methodology of Preconditioned JFNK Gravity Focusing Inversion
The gravity inversion involved discretizing the subsurface into a finite number of rectangular units, and the discrete map is shown in Figure 1. There are n observation points and m rectangular units, the forward expression of the gravity anomaly could be expressed as: where 1 n g × is the gravity anomaly, n m G × is the contribution to the nth datum of a unit density in the mth cell called kernel function matrix, and 1 m ρ × is the density matrix of rectangular units.  Performing gravity inversion requires solving Equation (1) where g φ is the square norm of the difference between the observed anomaly( obs g ) and the calculated anomaly( Gρ ), and it represents the fitting functional of the data. In the expression, ρ φ is the stabilizing functional that constrained the inversion results to the real conditions. u is a regularization parameter. The value of the regularization parameter is usually 10 n , and the n is determined by a "trial and error" method. Wz is the depth-weighting function which was used to counteract the inherent decay of the kernel function [27]. z W ρ is the weighted solution variable, The kernel function matrix is [25]: µ ijk x i ln y j + r ijk + y j ln x i + r ijk − z k arctan In the expression, γ is the gravitational constant, (x, y, z) are the coordinates of the observation points, and (ξ 1 , η 1 , ζ 1 ) and (ξ 2 , η 2 , ζ 2 ) are the coordinates of the minimum and maximum corner points of the cells.
Performing gravity inversion requires solving Equation (1) for ρ m×1 . Since the number of underground cells was greater than the data, the Equation (1) was underdetermined.
We solved this underdetermined equation in the following manner [26]: where φ g is the square norm of the difference between the observed anomaly(g obs ) and the calculated anomaly(Gρ), and it represents the fitting functional of the data. In the expression, φ ρ is the stabilizing functional that constrained the inversion results to the real conditions. u is a regularization parameter. The value of the regularization parameter is usually 10 n , and the n is determined by a "trial and error" method. W z is the depth-weighting function which was used to counteract the inherent decay of the kernel function [27]. W z ρ is the weighted solution variable, we replaced it with ρ w . The Equation (4) is the minimum support function [4]. In general, the e is 0.1, which is a focusing parameter to determine the sharpness of the inversion results. The derivative formula of Equation (3) is: We can substitute the solution of Equation (3) with the minimization of Equation (5). Finally, we used the JFNK method to obtain the ρ w of Equation (5), and the results ρ can be obtained by removing the depth-weighting function.
The JFNK method was introduced as follows. The f (x) = 0 can be processed as: The Krylov subspace method was created to solve the Equation (7), and we used the CG method which is a common solution technique with the Krylov subspace [26]. The matrix-vector product required by CG method was approximated by forward finite difference schemes: where v is a vector of the Krylov subspace. The process of the preconditioned CG method used to solve Equation (7) is: where k is the number of iterations, and as the number of iteration increased, the residual r was minimized by moving a distance α in search direction p. The S is a preconditioner of the search direction, and we could obtain the new search direction. The optimum preconditioning is (A T A) -1 . It produced the least-squares solution to Equation (7) in a single iteration. Since we were unable to obtain the optimum preconditioning directly, we substituted the optimum preconditioning with a depth-weighting function [10]. All the other parameters were intermediate variables with no significance.
Under the condition of the natural surface, the structured mesh method was not applicable, and the unstructured mesh method w triangulation (shown in Figure 2a). The structured method is computationally intensive. We proposed the adaptive equivalent-dimensional mesh method which combines the unstructured and structure mesh to achieve gravity inversion with the undulating terrain. Only the surface layer (the yellow area in Figure 2b) was discretized by the triangular unstructured mesh method, and the other subsurface was discretized into structure rectangular mesh. The three blue cells in Figure 2b have the same height, and it is shown that the structure rectangular mesh was the n dimension in each horizontal position. Therefore, many judgment processes were avoided in the calculation, and the calculation efficiency was improved; the memory was not increased in the case of fitting the surface as far as possible. This method could efficiently and accurately obtain the 3D density models, so it was more applicable to the gravity inversion with a large amount of data.
preconditioning with a depth-weighting function [10]. All the other parameters were intermediate variables with no significance.
Under the condition of the natural surface, the structured mesh method was not applicable, and the unstructured mesh method w triangulation (shown in Figure 2a). The structured method is computationally intensive. We proposed the adaptive equivalent-dimensional mesh method which combines the unstructured and structure mesh to achieve gravity inversion with the undulating terrain. Only the surface layer (the yellow area in Figure 2b) was discretized by the triangular unstructured mesh method, and the other subsurface was discretized into structure rectangular mesh. The three blue cells in Figure 2b have the same height, and it is shown that the structure rectangular mesh was the n dimension in each horizontal position. Therefore, many judgment processes were avoided in the calculation, and the calculation efficiency was improved; the memory was not increased in the case of fitting the surface as far as possible. This method could efficiently and accurately obtain the 3D density models, so it was more applicable to the gravity inversion with a large amount of data. As the adaptive equivalent-dimensional mesh method, we derived the equations to determine intervals (dx, dy) of the mesh according to the terrain and anomaly:  As the adaptive equivalent-dimensional mesh method, we derived the equations to determine intervals (dx, dy) of the mesh according to the terrain and anomaly: In the expression, X and Y represent the total length of the measuring area along x and y directions. THD g and THD h represent the total horizontal derivative (THD) of gravity anomaly g and elevation of measure points h, respectively. Its equation is shown in Equation (12). In this way, the dx and dy of the cells would be 0.5~1.5 times the average length. Therefore, the intervals of mesh were close in the area where the gravity anomaly or terrain varied greatly.

Gravity Model Tests
We tested our method in the inversion consisting of typical models of 150 × 150 × 100 m with undulating terrain, and the results are shown in Figure 3. The buried depth of the center of models Minerals 2020, 10, 741 6 of 13 was 100 m from the surface, and the central coordinates of two models were (225, 250) m and (725, 250) m. As we can see in Figure 3b, the extreme value position of the gravity effect was asymmetric because of the asymmetry of the undulating terrain. Therefore, the x-coordinates of the two extreme value positions of the gravity anomaly were 200 and 800 m, respectively. There was deviation from the x-coordinates of 225 and 725 m of the two actual models. When the terrain fitting accuracy was consistent, the terrain grid size was positively correlated with the fineness of the subsurface mesh. Figure 3c shows the run time of programs in different terrain grid sizes with twenty-five obverse points. It is shown in Figure 3c that the run time of programs with equivalent-dimensional mesh was obvious shorter than with triangulation mesh. elevation of measure points h, respectively. Its equation is shown in Equation (12). In this way, the dx and dy of the cells would be 0.5~1.5 times the average length. Therefore, the intervals of mesh were close in the area where the gravity anomaly or terrain varied greatly.

Gravity Model Tests
We tested our method in the inversion consisting of typical models of 150 × 150 × 100 m with undulating terrain, and the results are shown in Figure 3. The buried depth of the center of models was 100 m from the surface, and the central coordinates of two models were (225, 250) m and (725, 250) m. As we can see in Figure 3b, the extreme value position of the gravity effect was asymmetric because of the asymmetry of the undulating terrain. Therefore, the x-coordinates of the two extreme value positions of the gravity anomaly were 200 and 800 m, respectively. There was deviation from the x-coordinates of 225 and 725 m of the two actual models. When the terrain fitting accuracy was consistent, the terrain grid size was positively correlated with the fineness of the subsurface mesh. Figure 3c shows the run time of programs in different terrain grid sizes with twenty-five obverse points. It is shown in Figure 3c that the run time of programs with equivalent-dimensional mesh was obvious shorter than with triangulation mesh.  The 3D gravity inversion in Figure 3d was performed based on equivalent-dimensional mesh. Figure 3d shows the slice of the inversion result at y = 250 m, and the red curve is consistent with the undulating terrain and represents the top of the subsurface, and the black rectangles represent real positions of typical models. It can be seen in Figure 3d that the high-density units in the inversion result were similar to the actual range of models. Therefore, this method could obtain the information of models directly by the inversion with equivalent-dimensional mesh. Figure 3e shows the 3D gravity inversion with triangulation mesh, and the result of the inversion with equivalent-dimensional mesh was similar to the inversion with triangulation mesh, comparing Figure 3d,e. Therefore, the new meshing method proposed by us can obtain gravity inversion results consistent with the traditional method in the regions with undulating terrain efficiently. Figure 4a shows the gravity anomaly of Figure 3a with noise, and the anomaly showed some distortion after the addition of 30 signal noise ratio noise. Figure 4b is the slice of the inversion result, and Figure 4c is the 3D view of slices. In Figure 4d, blue lines represent the edge of the models, and the yellow and blue blocks represent densities greater than 0.28 and 0.16 g/cm 3 in the inversion results,  Figure 3c. Therefore, the method was stabilized and could obtain an accurate range of models after adding noise. dimensional mesh was similar to the inversion with triangulation mesh, comparing Figure 3d,e. Therefore, the new meshing method proposed by us can obtain gravity inversion results consistent with the traditional method in the regions with undulating terrain efficiently. Figure 4a shows the gravity anomaly of Figure 3a with noise, and the anomaly showed some distortion after the addition of 30 signal noise ratio noise. Figure 4b is the slice of the inversion result, and Figure 4c is the 3D view of slices. In Figure 4d, blue lines represent the edge of the models, and the yellow and blue blocks represent densities greater than 0.28 and 0.16 g/cm 3 in the inversion results, respectively. The inversion results shown in Figure 4b-d are close to the result in Figure 3c. Therefore, the method was stabilized and could obtain an accurate range of models after adding noise. We processed a gravity anomaly with complex terrain to verify the method of equivalentdimensional mesh. There were two peaks and one valley-shown in Figure 5b. In the 100 m depth from the surface, there were three prisms of 100 × 100 × 100 m; the central coordinates of which were (250, 250), (250, 750), and (750, 500) m, respectively. Figure 5a is the gravity anomaly of models, and the gravity containing 30 SNR noise is shown in Figure 5c. Because of the complex terrain's effect, the coordinates of the extreme value points shown in Figure 5a were (250, 200), (250, 800), and (750, 550) m, and the coordinates were not consistent with the models' coordinates. We achieved the density inversion based on the mesh with the equivalent dimension and obtained the slice diagram of the inversion results shown in Figure 5e-h by cutting along the white dashed lines numbered 1 and 2 shown in Figure 5a-c. Figure 5d shows the high-density blocks of a density greater than 0.3 g/cm 3 , and the blue lines represent the edge of the models. We processed a gravity anomaly with complex terrain to verify the method of equivalentdimensional mesh. There were two peaks and one valley-shown in Figure 5b. In the 100 m depth from the surface, there were three prisms of 100 × 100 × 100 m; the central coordinates of which were (250, 250), (250, 750), and (750, 500) m, respectively. Figure 5a is the gravity anomaly of models, and the gravity containing 30 SNR noise is shown in Figure 5c. Because of the complex terrain's effect, the coordinates of the extreme value points shown in Figure 5a were (250, 200), (250, 800), and (750, 550) m, and the coordinates were not consistent with the models' coordinates. We achieved the density inversion based on the mesh with the equivalent dimension and obtained the slice diagram of the inversion results shown in Figure 5e-h by cutting along the white dashed lines numbered 1 and 2 shown in Figure 5a-c. Figure 5d shows the high-density blocks of a density greater than 0.3 g/cm 3 , and the blue lines represent the edge of the models.
Although the extreme values' coordinates of the gravity anomaly deviated from the model positions at lines 1 and 2, the central position of the inversion results obtained by the gravity inversion based on equivalent-dimensional mesh were still consistent with the model position. The results shown in Figure 5d-h indicate that this method is applicable and stabilized in complex terrain.  Figure 5d-h indicate that this method is applicable and stabilized in complex terrain.

Actual Data Processing
In order to validate the applicability of the adaptive equivalent-dimensional mesh method in the field data, we processed the gravity data of polymetallic mining areas in Shaanxi province of China. The area is located in Huayangchuan in Shaanxi province. Its tectonic location belongs to Xiaoqinling intracontinental orogenic belt in the southern margin of north China landmass. The main outcrops are the Neoarchean erathem Taihua group (deep metamorphic crystalline complex), the Mesoproterozoic erathem Xionger group (volcanic-sedimentary), Gaoshanhe formation (clasolite), and the Luonan group (carbonatite); some areas have the Sinian series and the Cambrian system clasolite. Paleogene, Neogene, and Quaternary are deposited in Cenozoic basins, and the Ordovician, Silurian, and Paleozoic to Mesozoic stratums are missing [28,29]. The geological condition of the survey area is shown in Figure 6a.
clasolite. Paleogene, Neogene, and Quaternary are deposited in Cenozoic basins, and the Ordovician, Silurian, and Paleozoic to Mesozoic stratums are missing [28,29]. The geological condition of the survey area is shown in Figure 6a. Figure 6b shows the Bouguer gravity anomaly of the survey area. The center area obviously had a high-value gravity anomaly in Figure 6b, so there was a high-density body in the subsurface. There were polymetallic mines in the vicinity of this area, and the geophysical characteristics of the metal minerals were high values of density. Therefore, we inferred that the high-value gravity in this area was a polymetallic deposit. Metallic minerals were produced by magma intruding into the shallow strata of the subsurface through the faults, so the location of faults had great value for us to infer the range of minerals.   Figure 6b shows the Bouguer gravity anomaly of the survey area. The center area obviously had a high-value gravity anomaly in Figure 6b, so there was a high-density body in the subsurface. There were polymetallic mines in the vicinity of this area, and the geophysical characteristics of the metal minerals were high values of density. Therefore, we inferred that the high-value gravity in this area was a polymetallic deposit. Metallic minerals were produced by magma intruding into the shallow strata of the subsurface through the faults, so the location of faults had great value for us to infer the range of minerals.
The total horizontal derivative, of which the maximum value indicates the location of the faults, is commonly used to detect faults in gravity interpretation [30]. Therefore, by comprehensively analyzing the geological map and the gravity data of this region, we finally divided the faults of this region into the total horizontal derivative and local gravity field map shown in Figure 7a. The total horizontal derivative, of which the maximum value indicates the location of the faults, is commonly used to detect faults in gravity interpretation [30]. Therefore, by comprehensively analyzing the geological map and the gravity data of this region, we finally divided the faults of this region into the total horizontal derivative and local gravity field map shown in Figure 7a. The Bouguer gravity anomaly is the comprehensive response of the density at all depths. The 3D density inversion depth range was an altitude of −1500 m from the surface, so we needed to perform separation of the potential field by using the matched filter method. The expression of the matched filter is [31,32]: where ( ) LnE ω is the average logarithmic power of the potential field, H is the depth of the layer corresponding to the anomaly, ω is the wave number, and B is a constant related to the deep or shallow source anomaly.
In this way, we removed the regional field anomaly corresponding to the density below −1500m and extracted the shallow response of the gravity anomaly, namely the local field. The result is shown in Figure 7b, indicating that the local anomalies in this area were separated by faults; this is consistent The Bouguer gravity anomaly is the comprehensive response of the density at all depths. The 3D density inversion depth range was an altitude of −1500 m from the surface, so we needed to perform separation of the potential field by using the matched filter method. The expression of the matched filter is [31,32]: where LnE(ω) is the average logarithmic power of the potential field, H is the depth of the layer corresponding to the anomaly, ω is the wave number, and B is a constant related to the deep or shallow source anomaly.
In this way, we removed the regional field anomaly corresponding to the density below −1500 m and extracted the shallow response of the gravity anomaly, namely the local field. The result is shown in Figure 7b, indicating that the local anomalies in this area were separated by faults; this is consistent with the geological characteristics of the ores distributed around the faults, which shows the accuracy of the matched filter method. The terrain in this region was complex with an average altitude of 1480 m and a maximum height difference of 509 m, as shown in Figure 8a. Considering the regional terrain relief, our adaptive equivalent-dimensional mesh method was suitable for this real case.
with the geological characteristics of the ores distributed around the faults, which shows the accuracy of the matched filter method. The terrain in this region was complex with an average altitude of 1480 m and a maximum height difference of 509 m, as shown in Figure 8a. Considering the regional terrain relief, our adaptive equivalent-dimensional mesh method was suitable for this real case. The depth of −1500m extending to surface of this region was meshed by the proposed method in this paper, and 3D gravity inversion was performed to obtain the density in the subsurface, as shown in Figure 8c, d. Figure 8b is the local field of gravity, in which the white dotted lines represent the horizontal position of the slices, and Figure 8d shows the slices of the 3D density model obtained from the gravity inversion results shown in Figure 8c. The density model in each section show the high-density units, which represent the position of the explained minerals.
Through the recognition of high-density units in the 3D density inversion model, the range of metallic ore bodies in this region was interpreted in Figure 9 by a red range. The ore bodies distributed near the faults that conform to the principle of the metal minerals were generated in the magma which intruded along the faults. There were six ore bodies explained in this area, and the number shown in Figure 9 around the red range is the height above sea level of these ore bodies. The depth of −1500 m extending to surface of this region was meshed by the proposed method in this paper, and 3D gravity inversion was performed to obtain the density in the subsurface, as shown in Figure 8c,d. Figure 8b is the local field of gravity, in which the white dotted lines represent the horizontal position of the slices, and Figure 8d shows the slices of the 3D density model obtained from the gravity inversion results shown in Figure 8c. The density model in each section show the high-density units, which represent the position of the explained minerals.
Through the recognition of high-density units in the 3D density inversion model, the range of metallic ore bodies in this region was interpreted in Figure 9 by a red range. The ore bodies distributed near the faults that conform to the principle of the metal minerals were generated in the magma which intruded along the faults. There were six ore bodies explained in this area, and the number shown in Figure 9 around the red range is the height above sea level of these ore bodies.
Minerals 2020, 10, x FOR PEER REVIEW 12 of 14 Figure 9. The range of ore bodies interpreted from 3D inversion results.

Conclusions
We proposed the preconditioned JFNK method to perform the intensive inversion computation; therefore, the 3D focusing inversion algorithm became more efficient. In addition, we combined unstructured and structured meshes to achieve gravity inversion with undulating terrain. Therefore, Figure 9. The range of ore bodies interpreted from 3D inversion results.

Conclusions
We proposed the preconditioned JFNK method to perform the intensive inversion computation; therefore, the 3D focusing inversion algorithm became more efficient. In addition, we combined unstructured and structured meshes to achieve gravity inversion with undulating terrain. Therefore, the computational complexity is less than that of the unstructured mesh method in the forward operation of the kernel function. In order to maintain the efficiency, we improved the terrain fitting and fineness of meshes by the adaptive equivalent-dimensional method. The fineness of meshes is positively correlated with the gradient of gravity and terrain. Synthetic tests showed that our method was suitable for gravity inversion under undulating terrain, and the algorithm was fast with accurate inversion results. Finally, we applied it to the field data processing in Huayangchuan, Shaanxi Province, and the distribution range and depth of the ore bodies were inferred from recovered density models, which further verified the stability and practicability of the new method.