An Accelerated Algorithm for 3D Inversion of Gravity Data Based on Improved Conjugate Gradient Method

: The 3D inversion algorithm for gravity data based on a smooth model constraint has been proven to yield a reasonable density distribution. However, as the amount of observed data and model parameters increases, the algorithm experiences issues with high memory consumption and prolonged computation time. Therefore, the corresponding problem in interpreting gravity inversion lies in developing a fast inversion algorithm. The conventional smooth model constraint inversion algorithm, based on regularization theory, requires the introduction of a model weighting function with a large matrix, and involves storage and operation of a large matrix with intermediate variables during inversion iteration, contributing signiﬁcantly to the prolonged computation time. In this paper, a diagonal weight matrix is represented by vectorization, and the intermediate variable of the large matrix type in the iteration is replaced with the combination of a small matrix and a vector. Additionally, the algorithm ﬂow of the conjugate gradient method is further optimized to minimize the number of vectors that need to be stored during iteration. As a result of these optimizations, the memory consumption of the algorithm during the operation process is successfully reduced. Finally, the experiments demonstrate the successful development of a fast 3D inversion algorithm for gravity data. Speciﬁcally, for a 80 × 80 × 20 mesh number inversion, our accelerated algorithm achieves an average speed of ~0.5 s per iteration, and the iterative process speeds up by a factor of 1000, providing an effective strategy for the fast inversion of large-scale data.


Introduction
The 3D inversion of gravity data can determine the geometric location of subterranean anomalous bodies, and can quantitatively calculate their shape, volume, size and density distribution values, providing a strong basis for subsequent geological interpretation [1][2][3][4].The inversion of gravity data is often ill-conditioned, and Tikhonov regularization theory is typically used to reduce the multiple solutions of inversion and ensure the uniqueness and stability of inversion results.Under the regularization theory, a lot of work has been completed to address the issue of non-uniqueness in inversion.Last and Kubik [5] introduced the minimum volume condition on the anomalous body as a means to reduce non-uniqueness in the inversion and named it the compact inversion method.Li and Oldenburg [1] introduced the roughness matrix and a depth weighting function into the objective function of gravity inversion to effectively address the skin effect problem and named this method smooth inversion.Portniaguine and Zhdanov [6] introduced the minimum gradient support function, based on the compact inversion process, to obtain the distribution of anomalous bodies with sharp boundaries.This method was named focused inversion and was further improved in subsequent years [7][8][9].Boulanger and Chouteau [10] examined and discussed the setting of mesh size in the inversion process by analyzing the kernel function and applied constraints such as minimum distance, Appl.Sci.2023, 13, 10265 2 of 15 flatness and smoothness to gravity data inversion through the use of the Lagrange formula.Pilkington [11] organized the methods and techniques of 3D physical property inversion for gravity and magnetic data and proposed a magnetic data inversion method based on sparsity constraints.Smooth inversion is generally used to address the inversion problem of geological bodies with smooth boundaries, while focused inversion has been demonstrated to have improved results in restoring the 3D physical properties of geological bodies with sharp boundaries.Both of these two methods have been widely used in gravity and magnetic data inversion.
The 3D inversion algorithm for gravity data can produce a reasonable density distribution and provide reference for further geological interpretation.Based on the development of fast mobile platform exploration technology (airborne and shipborne), the volume of gravity data is increasing and the accuracy of measurements are becoming more reliable.However, as the amount of observed data and model parameters increases, the inversion process faces problems of high memory consumption and long calculation time.Therefore, it is necessary to develop fast algorithms that can handle large-scale data.Using more effective optimization techniques [12][13][14] or compression methods [10,15,16] to extract more relevant information from large data sets is essential for improving inversion efficiency.Given the relative maturity of preconditioned decomposition technology, the preconditioned conjugate gradient method (PCG) has been effectively used and has produced favorable results in geophysical applications [17][18][19].These methods can decrease the number of inversion iterations, thus reducing calculation time and enhancing inversion efficiency.Moreover, CUDA is a general parallel computing architecture introduced by NVIDIA, which has become a powerful tool for performing large-scale data parallel computing.Thus, CUDA has become a widely adopted tool for parallel computing particularly in the field of fast inversion of gravity or magnetic data [20][21][22][23][24].
Based on the traditional 3D inversion algorithm for gravity data that uses smooth constraints, a specific optimization scheme approach is proposed to address the low efficiency of the traditional conjugate gradient algorithm in processing large-scale data, and the application effect of the proposed method is comparable to that of the existing matrix compression methods.The proposed method not only retains the inversion capabilities of the original algorithm, but also eliminates the need for storing and processing a large matrix throughout the program execution, thus reducing the memory usage, increasing computational efficiency and achieving the goal of fast 3D inversion.We use Python tools to implement the acceleration algorithm and develop 3D gravity data inversion software, in which the user graphical interface (GUI), data management and visualization are constructed using PyQt5.Finally, this paper demonstrates the inversion performance of the accelerated algorithm using both synthetic and real measurement data.It also presents a comparison of the calculation time of the accelerated algorithm and the traditional algorithm for different mesh numbers, thus clearly explaining the improved calculation speed of our algorithm.

Smooth Inversion of Gravity Data Based on Conjugate Gradient Algorithm
In 3D inversion, the subsurface is typically divided into a finite number of prism units with equal sizes, uniform density distribution and close arrangements.The formula for gravity anomaly response generated by a single prism is as follows: Among them, G = 6.67 × 10 −11 m 3 •kg −1 •s −2 is the gravitational constant, ρ is the density contrast, (ξ 1 , η 1 , ζ 1 ), (ξ 2 , η 2 , ζ 2 ) are the diagonal coordinate of the prism and (x, y, z) is the coordinate of the observation point.Equation (1) has various analytical solutions.The solution used in this paper is [25]: Extracting the density ρ in Equation (2), that is, define the sensitivity matrix: Then, the relation between vector → ρ composed of all the density value of all prisms in the whole underground space and all observed data ∆ → g is: where, S ij represents the sensitivity matrix.Equation ( 4) can be expressed in matrix form as [6]: where, d represents the column vector of observation data, m is the model parameter column vector, S is the sensitivity matrix of type M × N and M << N.
Since M << N, meaning that Equation ( 5) is an underdetermined system of equations, and given that there is always inevitable noise interference in the observed data, the inversion problem is often ill-conditioned, and the solution obtained through inversion has multiple solutions.The Tikhonov regularization method, as a conventional processing method, can effectively solve this problem.In other words, the aforementioned problems can be transformed into finding the optimal solution of the objective function.The objective function is: In the formula, φ d represents data fitting term, φ m represents the model constraint term, α represents regularization parameter.Data fitting function and model constraint function are defined as: The data weighting matrix, W d and the model weighting matrix, W m , used in this paper are diagonal square matrices, with the diagonal element defined as [26]: (10) After incorporating the data weighting matrix and the model weighting matrix, Equation ( 6) can be re-expressed as: We can solve the problem in Equation ( 11) using the reweighting algorithm proposed by Zhdanov [7], as described in the following formulas: In this paper, the conjugate gradient algorithm is used to find the optimal solution m w , which minimizes the objective function in (12).The iterative process of the conjugate gradient algorithm requires the calculation of the first derivative of the objective function φ with respect to the weighted density parameter m w , and the formula for calculating the first derivative is given as follows: An adaptive method [7] is selected to select regularization parameters, where the convergence of the inversion results gradually decreases the contribution value of the model constraint function and the value of α is selected through iteration is: In the next iteration, the decay value of α is determined by Equation (18): Additionally, constrains on the upper and lower limits of the density parameter have been added to control the inversion results within a reasonable range, as expressed in the following formula: Based on the above considerations, the gravity data inversion algorithm based on the traditional conjugate gradient algorithm is presented as follows.

An Accelerated 3D Fast Inversion Algorithm for Gravity Data
The efficiency of a computer is related to the memory occupied by data.Under the condition of achieving the same calculation results, finding a way to make the data occupy less memory is the key for improving the calculation efficiency.For the traditional conjugate gradient algorithm of gravity data inversion, it can be observed that three large matrices of size N × N, Q, W m and W m −1 are involved in the inversion iteration process.Their existence consumes a considerable amount of computer memory.It is precisely due to their frequent involvement in the iterative calculation process that the algorithm's operating efficiency is very low.It also results in considerable memory usage.To address the issues present in the existing algorithm, we propose three improvements from the standpoint of reducing matrix dimension and optimizing the calculation flow of the algorithm.(a) the large matrix Q involved in the traditional algorithm appears to simplify the calculation process by calculating and storing S T S in advance, but it actually increases the burden on the computer, and this approach is unnecessary.In fact, Q is the coefficient matrix of m w in Equation ( 16).Although it has mathematical significance, in computer operation, we only need to expand Equation ( 16) directly and replace the key item S T w S w m w in the calculation order.In other words, we can change the order of calculation by using the associative law and transform it to S T w (S w m w ), thus successfully avoiding the generation and participation of the large matrix in the computer calculation.Therefore, in the process of algorithm improvement, we no longer store Q in advance for step (4) of the traditional algorithm, but instead only redistribute the storage of q.And q is defined as other items excluding S T w (S w m w ) in the expansion of Equation ( 16), namely, q is redefined as: In this way, the calculation expression involving the negative gradient in step ( 5) and ( 8) is modified as follows: Essentially, the change made to the denominator of the search step length in step ( 10) is simply a rearrangement of the calculation order of parameters, resulting in the final improvement shown as follows: For the above improvements, although the formal expression becomes complicated, the actual computer operation process avoids N × N oversized array, thus improving the computation efficiency.
(b) We consider the improvement in compression for the large matrix W m .In fact, for the matrices W m and W m −1 , even though they are large matrices of type N × N, W m is essentially a diagonal square matrix with the useful information being distributed along the diagonal, and the values in other positions are 0. In other words, the valuable information in W m is 1  N of its total matrix size.As a result, when W m is stored as a matrix in a computer, it consumes a significant amount of memory space, even though the crucial information only takes up a fraction of the storage space, this leads to an inefficient utilization of memory.In this study, we suggest expressing the original large matrix W m in the form of vector using linear representation: The vector representation still retains all of the useful information from the original matrix and requires only 1  N of the storage space.We can do the same processing with M × M diagonal square W d : To convert the two weighted matrices W d M×M and W m N×N into their vector forms, simply let them serve the same purpose as the original two matrices.In fact, we can analyze this from the perspective of matrix transformation.In Equation ( 15), S w essentially performs a two-step transformation: initially applying a row transformation to the matrix S. To be more specific, it multiplies each row of the S M×N matrix by W d M×M (i, i) (i = 1, 2, . . .M). Subsequently, it conducts a column transformation by multiplying each column with W d N×N (j, j) (j = 1, 2, . . .N).This process yields the resulting matrix S w .Therefore, rewriting the weighted matrix in the algorithm into vector form can also achieve the purpose of Equation (15).We only need to multiply the ith row of matrix S M×N by the ith component accordingly, and then multiply the jth column of the obtained matrix by the jth component W m N×1 (j) (j = 1, 2, . . .N) of vector W m N×1 to get S w .Other processing methods involving the diagonal matrix to participate in the calculation are similar to Equation (15).In this way, we can rewrite the large matrix involved in the algorithm, especially W m , into vector form, effectively releasing computer memory and improving running speed.(c) In each iteration of the traditional algorithm, the long vector r k , p k with dimension N are calculated and stored several times, especially in the step (9), where two large vectors r k and r k−1 must be stored simultaneously; however, this is completely unnecessary.In fact, in further algorithm improvement, we only need to store the inner product value of the negative gradient vector calculated at each iteration and use it in the next iteration stage.In this way, the simultaneous storage of the two long vectors r k and r k−1 can be avoided, releasing more computer memory space, and making the algorithm run more efficiently.
In this paper, we make several improvements to give the optimized conjugate gradient method, which improves operation efficiency and achieves the goal of fast inversion.The improved algorithm is as follows.
Under the condition of the same equipment, the acceleration algorithm proposed in this paper is faster than the traditional algorithm and can handle a larger data scale.It is worth noting that the new algorithm is consistent with the traditional method in calculating the sensitivity matrix, so it is limited by the size of the sensitivity matrix like the traditional method.In other words, due to the limited performance memory of the computer device used, the size of the storage sensitivity matrix exceeds the limit, so the method in this paper will also be ineffective.

Software Overview
As the core of the detailed inversion algorithm described above, we have developed user-oriented gravity data inversion software using Python tools, in which the user graphical interface (GUI), data management and visualization are constructed using PyQt5.The software is divided into two modules: forward modeling and inversion.In the forward modeling module, users can build models to their specifications and obtain synthetic gravity data for relevant synthetic model test.In the inversion module, there are both traditional and acceleration algorithms, for comparison of the acceleration algorithm in efficiency.More importantly, users can use the acceleration algorithm module to quickly obtain inversion results for large scale gravity data.The performance test is conducted on a DELL R740 workstation (CPU: 2 Intel Xeon 6254; GPU: 1 NVIDIA T4; Memory: 16 × 64 GB) to verify the validity and practicality of the accelerated inversion method.The inversion software is available now upon request and will be accessible on the authors Github site (https://github.com/Ashendudu/soft,accessed on 4 September 2023).

Synthetic Data Testing
To demonstrate the computational efficiency of the proposed method, we designed the model based on the published article by Hou et al. [24].Regarding the accelerated algorithm we presented, we aim to compare the acceleration advantages of the proposed method with traditional algorithms and the GPU parallel algorithm proposed by Hou et al. [24].To better explain the impact of our method, we built two models to test the effect of the improved algorithm.The parameters of the models are shown in Table 1.In this paper, the steps (1) to (6) in Algorithms 1 and 2 are referred as the preprocessing stage, and the steps (7) to (11) are referred to as the iteration stage.
2 ) (3) Calculate the weighted parameter: 5) Calculate the initial search direction: p 0 = r 0 = q − Qm W0 , and calculate the initial search step size: Update the initial solution and calculate the error: Set the algorithm iteration condition: while loss > ε or k < N max : (8) Upset parameter Q, and then update negative gradient r k : Calculate the intermediate parameter β, and then update the search direction p k : (10) Update the search step size: 11) Update solution m k and return step (7) to determine whether it iterates: (3) Calculate the weighted parameter: 5) Calculate the initial search direction: p 0 = r 0 = q 0 − S T W (S W m W0 ), set β up = r T 0 r 0 , and calculate the initial search step size: t 0 = βup (S w p 0 ) T (S w p 0 +α 0 p T 0 p 0 (6) Update the initial solution and calculate the error: Set the algorithm iteration condition: while loss > ε or k < N max : (8) Upset parameter q k , and then update negative gradient r k : Calculate the intermediate parameter β, and then update the search direction p k : 11) Update solution m k and return step (7) to determine whether it iterates: In Model I, the underground space is divided into 40 × 40 × 20 = 32,000 regular prisms with a size of 100 m × 100 m × 100 m.The density contrast of the two prisms is 0.5 g/cm 3 , and they have dimensions of 1000 m × 1500 m × 600 m; the top buried depth of the prisms is 600 m, and the horizontal distance of two prisms is 1000 m as shown in Figure 1.The observation plane is located on z = 0, and the measurement area covers a range of 4000 m × 4000 m with a sampling interval of 100 m.The maximum inversion depth is set to 2000 m.To ensure the stability of the algorithm, 5% Gaussian random noise is added to the synthetic gravity data.Figure 2 shows the gravity anomaly data with and without noise.The proposed algorithm in this paper is applied to invert the gravity data anomaly with noise, and the inversion parameters were set as follows: the density constraint range was controlled at [−0.1, 0.6] g/cm 3 , the maximum number of iterations k = 100, and the cutoff error was set at ε = 0.002.After 17 inversion iterations, the data fitting difference reached 0.001, thus ending the iteration.In the inversion process, the pre-processing stage took 6.27 s, and the iteration process took 1.36 s, resulting in a single iteration time of 0.08 s.As shown in Figure 3, we select a cross section along y = 2000 m and z = 900 m to display the inversion results, in which the black boxes in the section represents the true boundary of the established model.It can be seen that the inversion results are in good agreement with the real model position, indicating that the improved algorithm has anti-noise performance.Since our algorithm is based on smooth inversion, the density values are lower than the real values of the model, but we can still determine the approximate range of the density distribution of the model relatively.
Appl.Sci.2023, 13, x FOR PEER REVIEW 9 of 16 0.5 g cm ⁄ , and they have dimensions of 1000 m 1500 m 600 m; the top buried depth of the prisms is 600 m, and the horizontal distance of two prisms is 1000 m as shown in Figure 1.The observation plane is located on  0, and the measurement area covers a range of 4000 m 4000 m with a sampling interval of 100 m.The maximum inversion depth is set to 2000 m.To ensure the stability of the algorithm, 5% Gaussian random noise is added to the synthetic gravity data.Figure 2 shows the gravity anomaly data with and without noise.The proposed algorithm in this paper is applied to invert the gravity data anomaly with noise, and the inversion parameters were set as follows: the density constraint range was controlled at 0.1, 0.6 g cm ⁄ , the maximum number of iterations k = 100, and the cutoff error was set at  0.002.After 17 inversion iterations, the data fitting difference reached 0.001, thus ending the iteration.In the inversion process, the pre-processing stage took 6.27 s, and the iteration process took 1.36 s, resulting in a single iteration time of 0.08 s.As shown in Figure 3    Appl.Sci.2023, 13, x FOR PEER REVIEW 9 of 16 0.5 g cm ⁄ , and they have dimensions of 1000 m 1500 m 600 m; the top buried depth of the prisms is 600 m, and the horizontal distance of two prisms is 1000 m as shown in Figure 1.The observation plane is located on  0, and the measurement area covers a range of 4000 m 4000 m with a sampling interval of 100 m.The maximum inversion depth is set to 2000 m.To ensure the stability of the algorithm, 5% Gaussian random noise is added to the synthetic gravity data.Figure 2 shows the gravity anomaly data with and without noise.The proposed algorithm in this paper is applied to invert the gravity data anomaly with noise, and the inversion parameters were set as follows: the density constraint range was controlled at 0.1, 0.6 g cm ⁄ , the maximum number of iterations k = 100, and the cutoff error was set at  0.002.After 17 inversion iterations, the data fitting difference reached 0.001, thus ending the iteration.In the inversion process, the pre-processing stage took 6.27 s, and the iteration process took 1.36 s, resulting in a single iteration time of 0.08 s.As shown in Figure 3     In order to further test the algorithm's inversion effect on the complex model, we established Model II as shown in Figure 4, and its specific parameters are shown in Table 1.The fourth and fifth prisms are combined together to form an L-shaped target.The observation range on both the x and y axes is 0 to 2000 m, the observation interval is 25 m and the maximum inversion depth is set as 1000 m.If the underground space is divided into 80 80 40, then the size of each prism unit is 25 m 25 m 25 m.Similarly, 5% Gaussian random noise is added to the synthetic gravity data (Figure 5), and the inversion parameters are set as follows: The density constraint range is controlled within 0.1,1.1 g cm ⁄ , the maximum number of iterations is set to k = 100 and the cutoff error is set to  0.002.After 22 inversion iterations, the inversion result is obtained, and the total inversion time was 359.01 s.In the inversion process, the pre-processing stage took 338.16s to run, and the iteration process took 20.85 s, resulting in an iteration time of 0.95 s.Four cross sections are shown in Figure 6, and the true boundaries of the prisms are represented by the black boxes.The results indicate that the density distribution is in good agreement with the boundaries represented by the black boxes.In order to further test the algorithm's inversion effect on the complex model, we established Model II as shown in Figure 4, and its specific parameters are shown in Table 1.The fourth and fifth prisms are combined together to form an L-shaped target.The observation range on both the x and y axes is 0 to 2000 m, the observation interval is 25 m and the maximum inversion depth is set as 1000 m.If the underground space is divided into 80 × 80 × 40, then the size of each prism unit is 25 m × 25 m × 25 m.Similarly, 5% Gaussian random noise is added to the synthetic gravity data (Figure 5), and the inversion parameters are set as follows: The density constraint range is controlled within [−0.1, 1.1]g/cm 3 , the maximum number of iterations is set to k = 100 and the cutoff error is set to ε = 0.002.After 22 inversion iterations, the inversion result is obtained, and the total inversion time was 359.01 s.In the inversion process, the pre-processing stage took 338.16s to run, and the iteration process took 20.85 s, resulting in an iteration time of 0.95 s.Four cross sections are shown in Figure 6, and the true boundaries of the prisms are represented by the black boxes.The results indicate that the density distribution is in good agreement with the boundaries represented by the black boxes.              2 and 3, respectively, record the calculation time comparisons for data of different scales in the pre-processing stage and the single calculation time comparisons in the iterative process.As can be seen from the results, due to limitations in computer performance (such as memory), the traditional inversion algorithm was not calculated for the 100 × 100 × 20 mesh number, but the new algorithm still achieved a substantial inversion speed under the same computer configuration.It is evident that the improved algorithm has a significant acceleration in the running time of both pre-processing stage and the iteration stage.In the pre-processing stage, the improved acceleration algorithm is typically more than 10 times faster than the traditional algorithm (Table 2), while in the iterative process, the acceleration algorithm has more obvious advantages.As can be seen from Table 3, the fastest acceleration ratio can reach over 1000 times, and the acceleration ratio in the iterative process of the acceleration algorithm typically exceeds 400 times, fully demonstrating the advantages of the improved algorithm in terms of inversion speed.We recorded the maximum memory footprint of 40 × 40 × 20 grid data in both traditional and accelerated algorithms.The maximum memory usage of the accelerated algorithm is 909.6 MiB.For the same number of grids, the maximum memory usage of the traditional algorithm is 24,295.9MiB.In addition, we also performed a simple comparison between our proposed acceleration algorithm and the GPU parallel acceleration inversion algorithm developed by Hou et al. [24].The GPU parallel algorithm developed by Hou et al. [24] takes approximately 6 s for a single iteration with the 100 × 100 × 20 mesh number.The time required for a single iteration of our algorithm at this partition scale is usually no more than 2 s.This explains our method can reduce computer memory usage and computation time for large-scale data inversion.

Real Data Testing
The proposed method is applied to the gravity anomaly data of the Gonghe Basin in the northeastern Qinghai-Tibet Plateau (Figure 7).The inner and peripheral areas of the Gonghe Basin show abundant geothermal resources [27].On the east side of the basin, the temperature of the Zhacang and Qunlianhai hot springs reached 93-96 • C, exceeding the local boiling point.In 2017, a high-temperature dry hot rock mass of 236 • C was discovered at 3705 m depth in Qiapuqia in the Gonghe Basin, northeast of the Tibetan Plateau.This is the first time that the shallowest and hottest hot dry rock has been discovered in China, further confirming that the geothermal resources in Gonghe Basin are rich and are attracting renewed interest in evaluating geothermal resources in this region [28].Several hot springs and geothermal wells are located in and around the basin [29,30], and a granite basement with hot dry rock is approximately 2 km underground [31,32].
area was 33 25 km as shown in Figure 7.In this paper, the inversion maximum depth is set at 9 km, and the underground space of the inversion target area is divided into 66 50 60 closely spaced upright prism units, each with a size of 500 m 500 m 150 m.The density constraint range of inversion is controlled at 1, 1 g cm ⁄ , the maximum number of iterations is set at k = 200, and the cutoff error is set at  0.01.After 110 iterations the inversion was terminated, with a total time of 221.8 s, including 114.7 s for the iteration time.Figure 9 illustrates the 3D imaging display of gravity survey data inversion results.The cap layer has weak density characteristics, and its thickness is generally less than 3 km.The results of our processing of real measured data are consistent with previous studies, indicating the practical application capabilities of our algorithm.The Institute of Geological Survey of Jilin University collected the gravity anomaly data of the Gonghe region used in this study (Figure 8).After grid processing, the measured data revealed that the grid spacing was 500 m, and the range of the entire measuring area was 33 × 25 km as shown in Figure 7.In this paper, the inversion maximum depth is set at 9 km, and the underground space of the inversion target area is divided into 66 × 50 × 60 closely spaced upright prism units, each with a size of 500 m × 500 m × 150 m.The density constraint range of inversion is controlled at [−1, 1]g/cm 3 , the maximum number of iterations is set at k = 200, and the cutoff error is set at ε = 0.01.After 110 iterations the inversion was terminated, with a total time of 221.8 s, including 114.7 s for the iteration time.Figure 9 illustrates the 3D imaging display of gravity survey data inversion results.The cap layer has weak density characteristics, and its thickness is generally less than 3 km.The results of our processing of real measured data are consistent with previous studies, indicating the practical application capabilities of our algorithm.

Conclusions
The gravity exploration method plays an important role in ore exploration and oil and gas exploration.For example, in the case of metal minerals, gravity exploration can be used to find ore directly (such as chromite, magnetite, etc.) according to the characteristics of high density.In the geological exploration of oil, gas and coal fields, gravity exploration can be used to delineate sedimentary basins, study the thickness of sedimentary layers, and find the structures related to oil, gas and coal fields according to the characteristics of low density.The 3D inversion algorithm for gravity data can produce a reasonable density distribution and provide reference for further geological interpretation.However, as the number of observed data and model parameters increases, the inversion process faces problems of high memory consumption and long calculation time.Therefore, it is necessary to develop fast algorithms that can handle large-scale data.
We study a 3D fast smoothing constraint algorithm for gravity data based on an optimized conjugate gradient flow.The algorithm is implemented in Python language, and PyQt gravity data inversion software is developed.In the data tests, the inversion obtains satisfactory 3D density distribution imaging results.It can be seen from the performance test records that the acceleration algorithm is scalable and can accelerate the inversion quickly.The new algorithm avoids generating an N × N large matrix in the whole operation process, reducing memory consumption, which leads to a faster speed in the iterative process.Compared to other GPU-based parallel algorithms, our proposed method is

Conclusions
The gravity exploration method plays an important role in ore exploration and oil and gas exploration.For example, in the case of metal minerals, gravity exploration can be used to find ore directly (such as chromite, magnetite, etc.) according to the characteristics of high density.In the geological exploration of oil, gas and coal fields, gravity exploration can be used to delineate sedimentary basins, study the thickness of sedimentary layers, and find the structures related to oil, gas and coal fields according to the characteristics of low density.The 3D inversion algorithm for gravity data can produce a reasonable density distribution and provide reference for further geological interpretation.However, as the number of observed data and model parameters increases, the inversion process faces problems of high memory consumption and long calculation time.Therefore, it is necessary to develop fast algorithms that can handle large-scale data.
We study a 3D fast smoothing constraint algorithm for gravity data based on an optimized conjugate gradient flow.The algorithm is implemented in Python language, and PyQt gravity data inversion software is developed.In the data tests, the inversion obtains satisfactory 3D density distribution imaging results.It can be seen from the performance test records that the acceleration algorithm is scalable and can accelerate the inversion quickly.The new algorithm avoids generating an N × N large matrix in the whole operation process, reducing memory consumption, which leads to a faster speed in the iterative process.Compared to other GPU-based parallel algorithms, our proposed method is equally efficient computationally.Moreover, the calculation time required for a single iteration is at least 350 times faster than that of the traditional algorithm, making our proposed acceleration algorithm efficient and practical for large-scale gravity data inversion.Therefore, this paper provides a new strategy for the rapid inversion of large-scale gravity data, saves the computational cost (time and memory) of data processing and provides important guidance for further ore and oil and gas exploration.

Algorithm 2
The workflow of accelerated 3D fast inversion algorithm (1) Input: S, d Set: k = 0, m 0 = 0 (2) Calculate the weighted parameter matrix: , we select a cross section along y = 2000 m and z = 900 m to display the inversion results, in which the black boxes in the section represents the true boundary of the established model.It can be seen that the inversion results are in good agreement with the real model position, indicating that the improved algorithm has anti-noise performance.Since our algorithm is based on smooth inversion, the density values are lower than the real values of the model, but we can still determine the approximate range of the density distribution of the model relatively.

Figure 1 .
Figure 1.Model I with two prisms.

Figure 2 .
Figure 2. Model I related gravity data.(a) Gravity data; (b) gravity data added noise.

Figure 1 .
Figure 1.Model I with two prisms.
, we select a cross section along y = 2000 m and z = 900 m to display the inversion results, in which the black boxes in the section represents the true boundary of the established model.It can be seen that the inversion results are in good agreement with the real model position, indicating that the improved algorithm has anti-noise performance.Since our algorithm is based on smooth inversion, the density values are lower than the real values of the model, but we can still determine the approximate range of the density distribution of the model relatively.

Figure 1 .
Figure 1.Model I with two prisms.

Figure 2 .
Figure 2. Model I related gravity data.(a) Gravity data; (b) gravity data added noise.

Figure 2 .
Figure 2. Model I related gravity data.(a) Gravity data; (b) gravity data added noise.

Figure 3 .
Figure 3. Model I inversion results.(a) Cross section along y = 2000 m; (b) cross section along z = 900 m.

Figure 3 .
Figure 3. Model I inversion results.(a) Cross section along y = 2000 m; (b) cross section along z = 900 m.

Figure 5 .
Figure 5. Model II related gravity data.(a) Gravity data; (b) gravity data added noise.

Figure 5 .
Figure 5. Model II related gravity data.(a) Gravity data; (b) gravity data added noise.

Figure 5 .
Figure 5. Model II related gravity data.(a) Gravity data; (b) gravity data added noise.

Figure 5 .
Figure 5. Model II related gravity data.(a) Gravity data; (b) gravity data added noise.

Figure 6 .
Figure 6.Model II inversion results.(a) Cross section along y = 300 m; (b) cross section along y = 1700 m; (c) cross section along x = 300 m; (d) cross section along z = 500 m.Finally, to demonstrate the acceleration advantage of the improved algorithm compared to the traditional algorithm, we use the gravity data obtained from Model I and Model II to perform inversions with mesh numbers of 40 × 40 × 20, 50 × 50 × 20, 80 × 80 × 20,

Figure 7 .
Figure 7. (a) Distribution of faults, wells and hot springs in the Gonghe Basin and surrounding areas; the red boxes represent the study and inversion region (revised from Wang et al. [33]); (b) location of the research area.

Figure 7 .
Figure 7. (a) Distribution of faults, wells and hot springs in the Gonghe Basin and surrounding areas; the red boxes represent the study and inversion region (revised from Wang et al. [33]); (b) location of the research area.

16 Figure 8 .
Figure 8. Real measured gravity data corresponding to Gonghe Basin.Figure 8. Real measured gravity data corresponding to Gonghe Basin.

Figure 8 .
Figure 8. Real measured gravity data corresponding to Gonghe Basin.Figure 8. Real measured gravity data corresponding to Gonghe Basin.

Figure 8 .
Figure 8. Real measured gravity data corresponding to Gonghe Basin.

Figure 9 .
Figure 9. Real measured gravity data inversion results.

Figure 9 .
Figure 9. Real measured gravity data inversion results.

Table 1 .
Related parameters of Model I and Model II.

Table 2 .
Comparison of calculation time of different mesh number in the pre-processing stage.

Table 3 .
Comparison of calculation time of single iteration for different mesh number.