Three-dimensional gravity rapid inversion optimization method, system, storage medium and electronic device

Abstract

The present invention relates to a three-dimensional gravity rapid inversion optimization method, system, storage medium and electronic device. The method includes: calculating the forward kernel matrix of each depth layer in the three-dimensional model, combining every two forward kernel matrices and performing a fast Fourier transform on each combined forward kernel matrix to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix; based on the size of the forward kernel matrix, padding zeros to extend a two-dimensional vector and performing a fast Fourier transform to obtain a two-dimensional matrix spectrum; successively substituting each combined forward kernel matrix spectrum and the two-dimensional matrix spectrum into a preset composite spectrum calculation formula to obtain G<supgt;T< / supgt;u. The present invention improves the calculation efficiency through the symmetry of the fast Fourier transform; at the same time, it further saves time. According to the characteristics of the method flow, parallel computing can also be used to further improve the calculation efficiency, reduce the calculation time in actual work, and increase the practicability of three-dimensional gravity inversion.

Description

Background Art

[0002] The inversion optimization algorithm is the conjugate gradient algorithm (CG). The CG algorithm is one of the fast iterative algorithms for solving linear equations. Compared with the steepest descent method and the Newton method, the CG algorithm has more advantages in terms of convergence speed. The storage cost and computational cost required by the CG algorithm are relatively low. Since only matrix-vector products and vector inner product operations are needed and there is no need to explicitly generate large matrices, it is widely used in geophysical and other field inversions.

[0003] It can be found in the CG algorithm that there is no need to explicitly generate and store the inverse matrix of the coefficient matrix during the calculation process, which can greatly save calculation time and storage space for large-scale gravity data calculation. However, when performing actual three-dimensional inversion calculations, as the amount of observed data or the number of discrete grids increases, even the matrix in the data space will become quite large and the calculation will still be very difficult. Therefore, there is an urgent need to propose an inversion optimization method to solve the above problems. Summary of the Invention

[0004] To solve the above technical problems, the present invention provides a three-dimensional gravity fast inversion optimization method, system, storage medium, and electronic device.

[0005] The technical solution of the three-dimensional gravity fast inversion optimization method of the present invention is as follows:

[0006] Calculate the forward kernel matrix of each depth layer in the three-dimensional model And combine every two forward kernel matrices to obtain multiple combined forward kernel matrices

[0007] For each combined forward kernel matrix Perform a fast Fourier transform to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix

[0008] Based on the size of the forward kernel matrix, perform zero-padding extension and fast Fourier transform on the two-dimensional matrix in sequence to obtain the two-dimensional matrix spectrum

[0009] Substitute each combined forward kernel matrix spectrum and the two-dimensional matrix spectrum into a preset composite spectrum calculation formula, obtain and based on the and of every two depth layers, obtain G T u, where G T u is an intermediate parameter in the three-dimensional gravity inversion process, k is the k-th depth layer of the three-dimensional model, l is the l-th depth layer of the three-dimensional model, and j 2= -1。

[0010] The beneficial effects of the three-dimensional gravity rapid inversion optimization method of the present invention are as follows:

[0011] Based on the three-dimensional gravity rapid inversion method, the three-dimensional gravity rapid inversion optimization method further proposed by the present invention reduces the amount of calculation according to the symmetry of the fast Fourier transform, further improves the calculation efficiency, and achieves an acceleration ratio of about 1.6 times. At the same time, when calculating the inversion of a large amount of gravity data, using this method will further save valuable time. According to the characteristics of the method process, parallel computing can also be used to further improve the calculation efficiency, reduce the calculation time in actual work, and increase the practicality of three-dimensional gravity inversion.

[0012] On the basis of the above solution, the three-dimensional gravity rapid inversion optimization method of the present invention can also be improved as follows.

[0013] Further, the preset composite spectrum formula is:

[0014]

[0015] Wherein, is the forward kernel matrix of the k-th depth layer in the three-dimensional model, is the forward kernel matrix of the l-th depth layer in the three-dimensional model, (1≤k,l≤p), p is the total number of depth layers in the three-dimensional model, F represents the forward transform of the fast Fourier transform, F -1 represents the inverse transform of the fast Fourier transform, u is the single-layer model, j 2 = -1, m is the number of grid divisions in the x (north) direction of the three-dimensional model, n is the number of grid divisions in the y (east) direction of the three-dimensional model, and H is the composite spectrum signal.

[0016] Further, the calculation process of and for every two depth layers is:

[0017]

[0018]

[0019] Wherein, Imag[H] is the imaginary part of the composite spectrum signal, and Real[H] is the real part of the composite spectrum signal.

[0020] Further, the specific form of G T u is:

[0021]

[0022] Further, perform a fast Fourier transform on each combined forward kernel matrix to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix After that, it further includes:

[0023] Perform a conjugate symmetry transform on each combined forward kernel matrix spectrum to obtain the transformed combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix spectrum

[0024] Based on the size of the forward kernel matrix of each depth layer, zero-padding extension is sequentially performed on the density matrix of each depth layer to obtain the extended density matrix corresponding to each density matrix, and based on the combination method of every two forward kernel matrices, every two extended model matrices are combined and subjected to a fast Fourier transform to obtain a plurality of combined density matrix spectra

[0025] Sequentially accumulate the product of the transformed combined forward kernel matrix spectra of every two combinations of the same depth layer and the corresponding combined density matrix spectra to obtain a first calculation result

[0026] Perform a conjugate symmetry transform on the first calculation result to obtain a second calculation result and substitute the first calculation result and the second calculation result into a preset forward spectrum formula to obtain and perform an inverse Fourier transform on the sum of the forward spectra of each depth layer of the three-dimensional model to obtain the forward result of the three-dimensional model.

[0027] Further, the preset forward spectrum formula is:

[0028]

[0029] where is the sum of the forward spectra of each depth layer of the three-dimensional model, k is the k-th depth layer of the three-dimensional model, l is the l-th depth layer of the three-dimensional model, j 2 = -1, (1 ≤ k, l ≤ p), p is the total number of depth layers in the three-dimensional model, F represents the forward transform of the fast Fourier transform, m is the number of grid partitions in the x (north) direction of the three-dimensional model, and n is the number of grid partitions in the y (east) direction of the three-dimensional model.

[0030] The technical solution of a three-dimensional gravity fast inversion optimization system of the present invention is as follows:

[0031] It includes: a first processing module, a second processing module, a third processing module, and an operation module;

[0032] The first processing module is used for: calculating the forward kernel matrix of each depth layer in the three-dimensional model and combining every two forward kernel matrices to obtain a plurality of combined forward kernel matrices

[0033] The second processing module is used for: for each combined forward kernel matrix performing a fast Fourier transform to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix

[0034] The third processing module is used for: based on the size of the forward kernel matrix, performing zero-padding extension and fast Fourier transform on the two-dimensional matrix in sequence to obtain the two-dimensional matrix spectrum performing zero-padding extension and fast Fourier transform on the two-dimensional matrix in sequence to obtain the two-dimensional matrix spectrum

[0035] The operation module is used for: sequentially substituting each combined forward kernel matrix spectrum and the two-dimensional matrix spectrum into a preset composite spectrum calculation formula, obtaining and according to every two depth layers and obtaining G T u, where G T u is an intermediate parameter in the three-dimensional gravity inversion process, k is the k-th depth layer of the three-dimensional model, l is the l-th depth layer of the three-dimensional model, and j 2 = -1.

[0036] The beneficial effects of the three-dimensional gravity fast inversion optimization system of the present invention are as follows:

[0037] Based on the three-dimensional gravity fast inversion method, the system of the present invention further proposes a three-dimensional gravity fast inversion optimization method. According to the symmetry of the fast Fourier transform, by reducing the amount of calculation, the calculation efficiency is further improved, and an acceleration ratio of about 1.6 times is achieved. At the same time, when performing the calculation of mass gravity data inversion, using this method will further save valuable time. According to the characteristics of the method flow, parallel computing can also be used to further improve the calculation efficiency, reduce the calculation time in actual work, and increase the practicability of three-dimensional gravity inversion.

[0038] On the basis of the above solution, the three-dimensional gravity fast inversion optimization system of the present invention can also be improved as follows.

[0039] Further, the preset composite spectrum formula is:

[0040]

[0041] Among them, is the forward kernel matrix of the k-th depth layer in the three-dimensional model, is the forward kernel matrix of the l-th depth layer in the three-dimensional model, (1 ≤ k, l ≤ p), where p is the total number of depth layers in the three-dimensional model, F represents the forward transform of the fast Fourier transform, and F -1 represents the inverse transform of the fast Fourier transform, u is the single-layer model, and j 2 = -1, m is the number of grid divisions in the x (north) direction of the three-dimensional model, n is the number of grid divisions in the y (east) direction of the three-dimensional model, and H is the composite frequency spectrum signal.

[0042] The technical solution of a storage medium of the present invention is as follows:

[0043] Instructions are stored in the storage medium. When a computer reads the instructions, the computer is caused to execute the steps of the three-dimensional gravity fast inversion optimization method of the present invention.

[0044] The technical solution of an electronic device of the present invention is as follows:

[0045] It includes a memory, a processor, and a computer program stored on the memory and executable on the processor. It is characterized in that when the processor executes the computer program, the computer is caused to execute the steps of the three-dimensional gravity fast inversion optimization method of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 is a schematic flow chart of the three-dimensional gravity fast inversion optimization method according to an embodiment of the present invention;

[0047] Figure 2 is a three-dimensional model diagram in the three-dimensional gravity fast inversion optimization method according to an embodiment of the present invention;

[0048] Figure 3 is a schematic diagram of the inversion result of the three-dimensional model in the three-dimensional gravity fast inversion optimization method according to an embodiment of the present invention;

[0049] Figure 4 is a schematic structural diagram of the three-dimensional gravity fast inversion optimization system according to an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0050] As Figure 1 shown, the three-dimensional gravity fast inversion optimization method according to an embodiment of the present invention includes the following steps:

[0051] S1. Calculate the forward kernel matrix of each depth layer in the three-dimensional model and combine every two forward kernel matrices to obtain a plurality of combined forward kernel matrices

[0052] Among them, the three-dimensional model is a three-dimensional model of the underground half-space. Figure 2 It is a structural diagram of the three-dimensional model of the underground half-space. The three-dimensional model is divided into M (M = m × n × p) prismatic grids, where m and n respectively represent the number of grid divisions in the x (north) direction and the y (east) direction, and p represents the number of grid divisions in the z direction (the total number of depth layers). The surface observed gravity has N (N = m × n) measuring points.

[0053] Among them, the forward kernel matrix is a matrix of (2m - 1) × (2n - 1).

[0054] Among them, the combination of every two forward kernel matrices is specifically as follows: By default, the forward kernel matrices of two depth layers are sequentially taken from top to bottom in the depth layer for combination; when the total number of layers of the three-dimensional model is odd, the calculation of the last remaining bottom depth layer adopts the original fast algorithm.

[0055] It should be noted that it can also be combined from bottom to top or any two depth layers, and there is no restriction here.

[0056] Specifically, take the forward kernel matrices of two layers and respectively combine the forward kernel matrices of the two layers into a combined forward kernel matrix

[0057] S2. Perform a fast Fourier transform on each combined forward kernel matrix to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix

[0058] Specifically, convert each combined forward kernel matrix by fast Fourier transform into a combined forward kernel matrix spectrum and then store it to obtain a plurality of combined forward kernel matrix spectra.

[0059] S3. Based on the size of the forward kernel matrix, perform zero-padding extension and fast Fourier transform on the two-dimensional matrix in sequence to obtain a two-dimensional matrix spectrum

[0060] Among them, the two-dimensional matrix has the same size as the forward kernel matrix, both are two-dimensional matrices, and the size is (2m - 1) × (2n - 1).

[0061] S4. Sequentially multiply each combined forward kernel matrix spectrum with the two-dimensional matrix spectrum Substitute into the preset composite spectrum calculation formula, obtain and according to every two depth layers and to get G T u, where G T u is an intermediate parameter in the three-dimensional gravity inversion process, k is the k-th depth layer of the three-dimensional model, l is the l-th depth layer of the three-dimensional model, and j 2 = -1.

[0062] Among them, substitute each combined forward kernel matrix spectrum and the two-dimensional matrix spectrum into the preset composite spectrum calculation formula for calculation, and obtain the and corresponding to two depth layers of the combined forward kernel matrix spectrum, and calculate layer by layer until the end to obtain G T u.

[0063] Among them, when using the original algorithm and the calculation amount is that of 2 inverse Fourier transforms. When using this method, because and can be regarded as known, so when using the method of the present invention, only 1 inverse Fourier transform calculation amount and auxiliary addition and subtraction operations are required. Compared with the original algorithm, the theoretical calculation amount is approximately reduced by half.

[0064] It should be noted that in actual calculation, for the convenience of calculation, generally take l = k + 1.

[0065] Preferably, the preset composite spectrum formula is:

[0066]

[0067] Among them, is the forward kernel matrix of the k-th depth layer in the three-dimensional model, is the forward kernel matrix of the l-th depth layer in the three-dimensional model, (1 ≤ k, l ≤ p), p is the total number of depth layers in the three-dimensional model, F represents the forward transform of the fast Fourier transform, and F -1 represents the inverse transform of the fast Fourier transform, u is the single-layer model, j 2 = -1, m is the number of grid divisions in the x (north) direction of the three-dimensional model, n is the number of grid divisions in the y (east) direction of the three-dimensional model, and H is the composite spectrum signal.

[0068] Preferably, the calculation process of and for every two depth layers is:

[0069]

[0070]

[0071] Wherein, Imag[H] is the imaginary part of the composite spectrum signal, and Real[H] is the real part of the composite spectrum signal.

[0072] Preferably, the specific form of the G T u is as follows:

[0073]

[0074] Preferably, after performing a fast Fourier transform on each combined forward kernel matrix to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix the following steps are further included:

[0075] Performing a conjugate symmetry transformation on each combined forward kernel matrix spectrum to obtain the transformed combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix spectrum

[0076] Based on the size of the forward kernel matrix of each depth layer, zero-padding expansion is sequentially performed on the density matrix of each depth layer to obtain the expanded density matrix corresponding to each density matrix, and based on the combination method of every two forward kernel matrices, every two expanded model matrices are combined and subjected to a fast Fourier transform to obtain a plurality of combined density matrix spectra

[0077] Sequentially accumulating the product of the transformed combined forward kernel matrix spectrum and the corresponding combined density matrix spectrum of every two combinations at the same depth layer to obtain a first calculation result

[0078] Performing a conjugate symmetry transformation on the first calculation result to obtain a second calculation result and substituting the first calculation result and the second calculation result into a preset forward spectrum formula to obtain and perform an inverse Fourier transform on the sum of the forward spectra of each depth layer of the three-dimensional model to obtain the forward result of the three-dimensional model.

[0079] Specifically, first calculate the forward kernel matrix of each depth layer then obtain the combined forward kernel matrix and convert the combined forward kernel matrix into a combined forward kernel matrix spectrum through FFT Finally, the transformed combined forward kernel matrix is obtained through the conjugate symmetric transformation of the array. and then stored. For the model matrices of each depth layer They are respectively zero-padded and extended to the extended model matrices of size (2m - 1)×(2n - 1), and then every two extended model matrices are combined to obtain multiple combined model matrices. And they are transformed into the spectrum of the combined model matrix through FFT.

[0080] According to the specific calculation steps, it can be seen that the optimization algorithm (preset forward spectrum formula) first sums and then performs array transformation, saving the computational amount of p - 1 array transformations. Therefore, the computational amount of the optimization algorithm is only the computational amount of p / 2 Fourier forward transforms and the computational amount of 1 array transformation. When the vertical number of layers p is much larger than 2, the computational time of 1 array transformation can be ignored. The gravity forward calculation according to Equation 2 - 8 can be understood as calculating the sum of the forward spectra of the density models of 2 density layers, and only 1 Fourier forward transform needs to be calculated. Therefore, the method of the present invention can further reduce the computational amount of gravity forward calculation and improve the computational efficiency. The actual computational acceleration effect will be verified and analyzed through model experiments.

[0081] Note that during actual calculation, generally l = k + 1 for convenience of calculation. When the total number of layers p of the model is odd, the last layer is calculated using the fast algorithm, and the remaining layers use the optimized fast algorithm.

[0082] Preferably, the preset forward spectrum formula is:

[0083]

[0084] where is the sum of the forward spectra of each depth layer of the three-dimensional model, k is the k-th depth layer of the three-dimensional model, l is the l-th depth layer of the three-dimensional model, j 2 = -1, (1 ≤ k, l ≤ p), p is the total number of depth layers in the three-dimensional model, F represents the forward transform of the fast Fourier transform, m is the number of grid divisions in the x (north) direction of the three-dimensional model, and n is the number of grid divisions in the y (east) direction of the three-dimensional model.

[0085] Specifically, the derivation process of the preset forward spectrum formula is as follows:

[0086] Take and of any two layers and of the corresponding layers respectively, and combine the two groups to obtain and Take the real and imaginary parts of the Fourier transform of these two sets of signals and get G Re , G Im , Re and ρ Im ,in,

[0087] Where: G Re and G Im Denote the combined forward kernel matrix The real and imaginary parts of the Fourier transform; ρ Re and ρ Im Represent the combined model matrix Real and imaginary parts of the Fourier transform; j 2 =-1. According to conjugate symmetry, correspond The conjugate even symmetric component of correspond The conjugate odd-symmetric component of .

[0088] For the convenience of expression, the conjugate symmetric transformation G′(u,v) of the two-dimensional matrix is ​​defined as:

[0089] Where: u=0,1,2,……,2m-2,v=0,1,2,……,2n-2; represents the conjugation of G.

[0090] According to the above definition, it is not difficult to find that the conjugate symmetric transformation of a two-dimensional matrix has the following properties:

[0091]

[0092]

[0093]

[0094] The same can be said and

[0095] Therefore, the sum of the spectra of the forward fields of the k-layer and l-layer density models is obtained as

[0096]

[0097] From the above formula, we can find that for the calculation of the convolution sum of two sets of signals, if the original gravity three-dimensional fast algorithm is used, then and When known, the main calculation amount is: calculation and The computational load of the forward Fourier transform is performed a total of 2 times. If the optimized 3D fast forward gravity algorithm of the present invention is adopted, then and are known, and the fast Fourier transform is directly used, then the main computational load is: calculating The computational load of the forward fast Fourier transform is performed a total of 1 time and from to The computational load of the array transformation is performed a total of 1 time. During actual calculation, the calculation time of 1 forward fast Fourier transform and the calculation time of 1 from to to have a small difference. Therefore, for the case where the number of layers p of the model is very small, the improvement effect of the optimized algorithm is limited. However, for the case where the number of layers p of the model is often much greater than 2, the sum of the forward spectra of the density models of each layer can be further simplified to obtain as: The second calculation result is obtained by conjugate symmetric transformation of the array for the first calculation result and the forward spectrum formula is preset according to the formula to obtain Finally, the forward result of the 3D model is obtained through inverse Fourier transform.

[0098] It should be noted that here is equivalent to G v , G v is explained in the following text.

[0099] In the present invention, mainly the calculation of G T u is optimized to reduce the calculation time. It should be noted that the inversion calculation is to find the solution of the minimum value of the objective function. The gradient at the minimum value of the objective function must be 0, and the fixed-point iteration equations of the model ρ can be obtained in the model space and the data space respectively. The calculation expressions of the two methods are as follows.

[0100] The model ρ in the model space is: ρ = ρ 0 +(G T D -1 G + αW -1 ) -1 G T D -1 (d - Gρ 0 ), the model ρ in the data space is: ρ = ρ 0 + αWG T (D + GWG T ) -1 (d - Gρ 0 ). The most difficult part of the inversion calculation in the model space and the data space is to find (G T D -1 G + αW-1 ) and (D + GWG T )'s inverse matrix.

[0101] The inversion optimization method in the present invention is the conjugate gradient algorithm (CG). The CG algorithm is one of the fast iterative algorithms for solving linear equations. Compared with the convergence speeds of the steepest descent method and Newton's method, the CG algorithm has more advantages. The storage cost and computational cost required by the CG algorithm are relatively low because only matrix-vector products and vector inner product operations are needed, and there is no need to explicitly generate large matrices. Therefore, it is widely applied to inversions in geophysics and other fields. The main process of the CG algorithm for gravity inversion is shown in Table 1.

[0102] Table 1 CG Algorithm

[0103]

[0104] In the table: m k is the result of the k-th iteration, r k is the gradient of the objective function, k is the number of iterations, p k represents the search direction of the iteration, and α k represents the step size of the iteration search direction.

[0105] It can be seen from the conjugate gradient algorithm that the main computational amount of a single conjugate gradient inversion is the product operation of matrix A and vector p k . When the model space inversion method is adopted, the coefficient matrix A is an M×M square matrix (G T D -1 G + αW -1 ), and when the data space inversion method is adopted, the coefficient matrix A is an N×N square matrix (D + GWG T ). It can be found in the CG algorithm that there is no need to explicitly generate and store the inverse matrix of the coefficient matrix A during the calculation process, which can greatly save computational time and storage space for large-scale gravity data calculations. However, when performing actual three-dimensional inversion calculations, as the amount of observed data or the number of discrete grids increases, even the matrix A in the data space will become quite large, and the calculation will still be very difficult.

[0106] Further analysis of the above conjugate gradient algorithm shows that the main computational amount is the product of matrix A and vector p. The computational forms involved in the calculation are Gv and G T u, where v is a column vector of M×1 and u is a column vector of N×1. The column vectors u and v are intermediate variables and do not have fixed physical meanings. For the calculation of the Gv term, if the vector v is regarded as the vector m, then this term calculation can be regarded as a forward calculation; for the calculation of the G T u term, it is necessary to calculate according to G TThe characteristics of the matrix are further analyzed. Since the observed data is regularly distributed in a grid on the horizontal plane, and the model dissection units and the data grid are in a one-to-one correspondence in the plane projection, the matrix G can be divided into p sub-matrices according to the depth layer of the model: G = [G 1 G 2 , …, G k G p 1 ≤ k ≤ p; where any coefficient G k of G i，j，k represents the gravitational response value of the j-th gravity observation point when the i-th model density in the model unit of the k-th layer is 1 unit, and the coefficient G j，i，k represents the gravitational response value of the i-th gravity observation point when the j-th model density in the model unit of the k-th layer is 1 unit. The model unit body and the gravity observation point are in a one-to-one correspondence in the horizontal position. According to the translational equivalence and interchange symmetry, G i，j，k = G j，i，k Therefore, G k is a symmetric matrix, that is For the operation of G T u, it can also be regarded as calculating the forward modeling values at different heights for the single-layer model u. Since only 1 Fourier forward transform of the vector u and p Fourier inverse transforms of G T u are required for the calculation of G k u, the time complexity of the calculation of G T u is the same as that of the forward modeling calculation of Gv, and it can also be regarded as 1 forward modeling calculation of the model.

[0107] Here, the calculation decomposition processes of two calculation methods in the model space and the data space are given. When performing inversion in the model space, the product of the matrix and the vector is Ap:

[0108]

[0109] In the formula: the vector p is a one-dimensional column vector, and its size is the same as that of the model vector ρ. The g and u in this formula are temporary intermediate variables. The Gp term and the G T u (only representing the computational amount of the current term) term can be regarded as 1 forward modeling calculation respectively. In the formula, u = D - 1 The g term can be regarded as weighting the forward modeling result. If the matrix D is a main diagonal matrix, the weighting calculation can be simplified to the corresponding multiplication of two two-dimensional arrays. The W -1 p term is equivalent to performing calculations such as weighting and derivation on the "model p". The D -1 g term and the W - 1 p term have a computational amount compared to G T D -1The computational load of the Gp term can be ignored. Therefore, the calculation of Ap, which is the most difficult part to store and calculate in the inversion calculation, can be mainly decomposed into two forward calculations.

[0110] When performing inversion in the data space, the product of the matrix and the vector Ap:

[0111]

[0112] In the formula: the vector p is a one-dimensional column vector, whose size is the same as that of the observed data vector g. The m and u in this formula are temporary intermediate variables. G T The Gp term and the Gu term can be regarded as one forward calculation respectively. The u = Wm term in the formula can be regarded as calculations such as weighting and differentiating the model. The Dp term in the formula is equivalent to weighting the forward result. The computational loads of the Dp term and the Wm term are compared with those of GWG T The computational load of the Gp term can be ignored. Therefore, the calculation of Ap in the data space inversion can also be mainly decomposed into two forward calculations. It can be seen from the single inversion decomposition of the conjugate gradient that whether using model space inversion or data space inversion, their computational loads can be mainly decomposed into two forward calculations, which is the key basis for realizing fast inversion.

[0113] Therefore, the present invention realizes fast inversion calculation and improves the efficiency of fast inversion by quickly solving G T u and Gv.

[0114] In order to verify the optimization of three-dimensional gravity inversion by the method of the present invention, the following method is used for verification. For example, regarding the use of memory space: When performing three-dimensional regularization inversion of large-scale gravity data, whether in the model space or in the data space, if the matrix A is explicitly generated, the memory requirement for conjugate gradient inversion calculation is extremely large. The geometric lattice method, the fast algorithm, and the optimized fast algorithm (the method of the present invention) do not need to explicitly generate the matrix A, so the memory usage is greatly reduced. The horizontal and vertical directions of the forward kernel matrix of the fast algorithm and the optimized fast algorithm are about twice that of the model matrix. The vertical layer number of the forward kernel matrix of the fast algorithm is the same as that of the model, both being p layers; while the vertical layer number of the forward kernel matrix of the optimized fast algorithm is about half of that of the model vertical layer number, but in order to ensure the calculation speed, the optimized algorithm generally requires two matrices, so the total layer number is still p layers. Because the forward kernel matrices of the fast algorithm and the optimized fast algorithm are both complex matrices, the memory requirement of the forward kernel matrices of the fast algorithm and the optimized fast algorithm is about 8 times that of the equivalent geometric lattice method.

[0115] For computational efficiency: Whether directly generating matrix A, or indirectly generating matrix A using the geometric lattice method or the optimization method of this paper, when using the calculation method of multiplying matrix A step by step with a vector, the computational cost of one CG iteration for the four methods is the same (Table 2), which is approximately twice the computational cost of one forward modeling of the model. Since the fast algorithm uses the FFT algorithm, the calculation speed increases exponentially, so the computational efficiency can be greatly improved. For the further improved optimization algorithm, when not considering the computational time consumed by matrix addition and subtraction operations, theoretically, compared with the time complexity of the fast algorithm, the time complexity of the optimized fast algorithm can be reduced by half, that is, the computational efficiency is doubled.

[0116] Table 2 Comparison of computational costs of various methods

[0117]

[0118] To verify whether the optimization algorithm of this paper has the ability to quickly process large-scale gravity data, a computational performance analysis was carried out based on theoretical model data. The factors affecting the computational cost of 3D gravity inversion mainly include: the accuracy requirements of the inversion results, the constraint conditions, and different inversion strategies. The common point of the above factors is that they affect the number of iterations of conjugate gradient inversion. When the grid division is fixed, the computational cost of each CG iteration is almost unchanged, and the main computational cost in CG iteration is 1 Gv calculation and 1 G T u calculation. For the computational efficiency of G T u, the computational efficiency of the 3D gravity fast forward modeling algorithm and the optimized 3D gravity fast forward modeling algorithm (the method of the present invention) is an exponential acceleration compared to the computational efficiency of the conventional forward modeling algorithm.

[0119] The experiments still used the 12 models in Table 3. The average time taken for 200 forward modelings of the 3D gravity fast forward modeling algorithm and the optimized 3D gravity fast forward modeling algorithm for all models was used as the evaluation basis. It can be seen from the numerical simulation results (Table 4) that the average value of the acceleration ratios of the two methods is 1.592, and the acceleration ratio of Model 12 is 1.474, and the average of the two is about 1.5. Therefore, based on the experimental results, it can be considered that the optimized 3D gravity fast forward modeling algorithm can achieve a computational efficiency of about 1.5 times.

[0120] It can be found that the optimized optimization algorithm has a higher acceleration for Gv than for G T u. Analyzing the reasons: The main computational cost of Gv is the forward Fourier transform, and the main computational cost of G T u is the inverse Fourier transform. However, comparing the computational time of the fast algorithm, it can be seen that when using the fast algorithm, the time taken to calculate G T u is generally less than the time taken to calculate Gv; when using the optimized fast algorithm, the time taken to calculate G TThe time taken for u is generally greater than that for calculating Gv. However, when the matrix reaches a certain scale, the built-in Fourier transform subroutine has no pattern in the relative lengths of time taken for the forward and inverse Fourier transforms. Therefore, it can be inferred that the factors affecting the speedup ratio should be caused by other steps in the calculation program, and this will not be delved into further here. What can be determined through experiments is that the optimized fast algorithm has a certain acceleration effect on the calculations of both Gv and G T u.

[0121] Table 3 Model Discretization Table

[0122]

[0123]

[0124] Table 4 Comparison Table of Computational Efficiencies of 3D Gravity Forward Modeling Algorithms

[0125]

[0126] The test computer has a main frequency of 2.20 GHz and a memory of 64 GB, which are the current mainstream conventional computing conditions. The inversion method is the spatial domain conjugate gradient focusing inversion method. The inclination angle of the model is 45°, the residual density is 1 g / cm3, and the number of observation data points is 1024×1024. The underground space discretization method is 1024×1024×312, with a total of 327155712 cubic grid cells, and the side length of the cube is 100 m. Only the depth weighting factor and the focusing factor (reweighting matrix) are added in the inversion calculation.

[0127] After 171 iterations of inversion calculation, the data fitting error decreased from the initial 100% to the final fitting error of 2.3%. Without considering the calculation time of the forward kernel matrix (1.126 h), the total time taken was approximately 4.37 hours, and the average time for 1 conjugate iteration was approximately 92 seconds. Figure 3 (a) and Figure 3 (b) are respectively the horizontal slice plan view and the vertical slice profile view of the inversion result. The red frame in the figure is the position of the inclined plate-like body model. Figure 3 (c) and Figure 3 (d) are respectively Figure 3 (a) and Figure 3 (b) partial detail display diagrams. It can be seen that the inversion result of the focusing inversion algorithm in the data space has a good degree of coincidence with the model body. From the numerical experiments, it can be seen that the optimized algorithm in this paper can quickly complete the 3D inversion of massive data without sacrificing the calculation accuracy, proving the effectiveness and feasibility of this method.

[0128] Based on the three-dimensional gravity rapid inversion method, the technical solution of this embodiment further proposes an optimized method for three-dimensional gravity rapid inversion. According to the symmetry of the fast Fourier transform, this method further improves the calculation efficiency by reducing the amount of calculation, achieving an acceleration ratio of about 1.6 times. At the same time, when calculating the inversion of a large amount of gravity data, using this method will further save valuable time. According to the characteristics of the method flow, parallel computing can also be used to further improve the calculation efficiency, reduce the calculation time in actual work, and increase the practicability of three-dimensional gravity inversion.

[0129] As Figure 4 shown, an optimized system 200 for three-dimensional gravity rapid inversion according to an embodiment of the present invention includes: a first processing module 210, a second processing module 220, a third processing module 230, and an operation module 240;

[0130] The first processing module 210 is configured to: calculate the forward kernel matrix of each depth layer in the three-dimensional model and combine every two forward kernel matrices to obtain a plurality of combined forward kernel matrices

[0131] The second processing module 220 is configured to: perform a fast Fourier transform on each combined forward kernel matrix to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix

[0132] The third processing module 230 is configured to: based on the size of the forward kernel matrix, perform zero-padding extension and fast Fourier transform on the two-dimensional matrix in sequence to obtain the two-dimensional matrix spectrum

[0133] The operation module 240 is configured to: sequentially substitute each combined forward kernel matrix spectrum and the two-dimensional matrix spectrum into a preset composite spectrum calculation formula, obtain and according to the and of every two depth layers to obtain G T u, where G T u is an intermediate parameter in the three-dimensional gravity inversion process, k is the k-th depth layer of the three-dimensional model, l is the l-th depth layer of the three-dimensional model, and j 2 = -1.

[0134] Preferably, the preset composite spectrum formula is:

[0135]

[0136] Where is the forward kernel matrix of the k-th depth layer in the three-dimensional model, is the forward kernel matrix of the l-th depth layer in the three-dimensional model, (1 ≤ k, l ≤ p), where p is the total number of depth layers in the three-dimensional model, F represents the forward transform of the fast Fourier transform, F -1 represents the inverse transform of the fast Fourier transform, u is the single-layer model, j 2 = -1, m is the number of grid divisions in the x (north) direction of the three-dimensional model, n is the number of grid divisions in the y (east) direction of the three-dimensional model, and H is the composite frequency spectrum signal.

[0137] Based on the three-dimensional gravity fast inversion method, the technical solution of this embodiment further proposes a three-dimensional gravity fast inversion optimization method. According to the symmetry of the fast Fourier transform, by reducing the amount of calculation, the calculation efficiency is further improved, achieving an acceleration ratio of about 1.6 times. At the same time, when performing the calculation of inverting a large amount of gravity data, using this method will further save valuable time. According to the characteristics of the method flow, parallel computing can also be used to further improve the calculation efficiency, reduce the calculation time in actual work, and increase the practicability of three-dimensional gravity inversion.

[0138] For the parameters and the steps of each module in the three-dimensional gravity fast inversion optimization system 200 of this embodiment to achieve the corresponding functions, reference can be made to the parameters and steps in the embodiment of a three-dimensional gravity fast inversion optimization method in the above text, which will not be elaborated here.

[0139] A storage medium provided by an embodiment of the present invention includes: instructions are stored in the storage medium. When a computer reads the instructions, the computer is made to execute the steps of the three-dimensional gravity fast inversion optimization method. Specifically, reference can be made to the parameters and steps in the embodiment of the three-dimensional gravity fast inversion optimization method in the above text, which will not be elaborated here.

[0140] Computer storage media such as: USB flash drives, external hard drives, etc.

[0141] An electronic device provided by an embodiment of the present invention includes a memory, a processor, and a computer program stored on the memory and executable on the processor. It is characterized in that when the processor executes the computer program, the computer is made to execute the steps of the three-dimensional gravity fast inversion optimization method. Specifically, reference can be made to the parameters and steps in the embodiment of the three-dimensional gravity fast inversion optimization method in the above text, which will not be elaborated here.

[0142] Those skilled in the art of the relevant technical field know that the present invention can be implemented as a method, a system, a storage medium, and an electronic device.

[0143] Accordingly, the present invention may be embodied in the following forms, namely: it may be entirely hardware, entirely software (including firmware, resident software, microcode, etc.), or a combination of hardware and software, which is generally referred to herein as a "circuit", "module" or "system". In addition, in some embodiments, the present invention may also be embodied in the form of a computer program product in one or more computer-readable media, which contain computer-readable program code. Any combination of one or more computer-readable media may be employed. The computer-readable media may be a computer-readable signal medium or a computer-readable storage medium. The computer-readable storage medium may be, for example, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination of the foregoing. More specific examples (non-exhaustive list) of the computer-readable storage medium include: an electrical connection having one or more wires, a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In this document, the computer-readable storage medium may be any tangible medium that contains or stores a program that can be used by or in connection with an instruction execution system, apparatus, or device. Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those of ordinary skill in the art can make changes, modifications, substitutions, and variations to the above embodiments within the scope of the present invention.

Claims

1. A three-dimensional gravity rapid inversion optimization method, characterized in that, comprising: Calculate the forward kernel matrix for each depth layer in the 3D model And combine every two forward kernel matrices to obtain multiple combined forward kernel matrices Perform a fast Fourier transform on each combined forward kernel matrix to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix Based on the size of the forward kernel matrix, perform zero-padding expansion and fast Fourier transform on the two-dimensional matrix sequentially to obtain the two-dimensional matrix spectrum Successively, the spectrum of each forward kernel matrix of the combination and the spectrum of the two-dimensional matrix are substituted into a preset composite spectrum calculation formula, and G and is obtained for each pair of depth layers and based on this, Gu T is obtained, where Gu T is an intermediate parameter in the three-dimensional gravity inversion process, k is the k-th depth layer of the three-dimensional model, l is the l-th depth layer of the three-dimensional model, and j 2 = -1.

2. The three-dimensional gravity rapid inversion optimization method according to claim 1, characterized in that, the preset composite spectrum formula is: Among them, is the forward kernel matrix of the k-th depth layer in the three-dimensional model, is the forward kernel matrix of the l-th depth layer in the three-dimensional model, (1 ≤ k, l ≤ p), where p is the total number of depth layers in the three-dimensional model, F represents the forward transform of the fast Fourier transform, and F -1 represents the inverse transform of the fast Fourier transform, u is the single-layer model, and j 2 = -1, m is the number of grid divisions in the x (north) direction in the three-dimensional model, n is the number of grid divisions in the y (east) direction in the three-dimensional model, and H is the composite frequency spectrum signal.

3. The three-dimensional gravity rapid inversion optimization method according to claim 2, characterized in that, For every two depth layers, and the calculation process is as follows: wherein, Imag[H] is the imaginary part of the composite spectrum signal, and Real[H] is the real part of the composite spectrum signal.

4. The three-dimensional gravity rapid inversion optimization method according to claim 3, characterized in that, The said G T The specific form of u is as follows:

5. The three-dimensional gravity rapid inversion optimization method according to any one of claims 1-4, characterized in that, For each combined forward kernel matrix perform a fast Fourier transform to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix After that, it further includes: For the spectrum of each combined forward kernel matrix perform a conjugate symmetry transformation to obtain the transformed combined forward kernel matrix spectrum corresponding to the spectrum of each combined forward kernel matrix Based on the size of the forward kernel matrix of each depth layer, successively zero-pad the density matrix of each depth layer to obtain an extended density matrix corresponding to each density matrix, and based on the combination method of every two forward kernel matrices, combine every two extended model matrices and perform fast Fourier transform to obtain multiple combined density matrix spectra Accumulate the spectra of the forward kernel matrices of the transformation combinations of every two combinations at the same depth layer in sequence with the spectra of the corresponding combined density matrices to obtain the first calculation result Perform a conjugate symmetry transformation on the first calculation result to obtain a second calculation result and substitute the first calculation result and the second calculation result into a preset forward Fourier spectrum formula to obtain and sum the forward Fourier spectra of each depth layer of the three-dimensional model Perform an inverse Fourier transform on it to obtain the forward result of the three-dimensional model.

6. The three-dimensional gravity rapid inversion optimization method according to claim 5, characterized in that, the preset forward spectrum formula is: Among them, is the sum of the forward spectra of each depth layer of the three-dimensional model, k is the k-th depth layer of the three-dimensional model, l is the l-th depth layer of the three-dimensional model, j 2 = -1, (1 ≤ k, l ≤ p), p is the total number of depth layers in the three-dimensional model, F represents the forward transform of the fast Fourier transform, m is the number of grid divisions in the x (north) direction of the three-dimensional model, and n is the number of grid divisions in the y (east) direction of the three-dimensional model.

7. A three-dimensional gravity rapid inversion optimization system, characterized in that, comprising: a first processing module, a second processing module, a third processing module and an operation module; The first processing module is used to: calculate the forward kernel matrix of each depth layer in the three-dimensional model and combine every two forward kernel matrices to obtain a plurality of combined forward kernel matrices The second processing module is used to: perform a fast Fourier transform on each combined forward kernel matrix to obtain the combined forward kernel matrix spectrum corresponding to each combined forward kernel matrix The third processing module is configured to: based on the size of the forward kernel matrix, perform zero-padding expansion and fast Fourier transform on the two-dimensional matrix in sequence to obtain the two-dimensional matrix spectrum The operation module is used to: sequentially substitute each of the combined forward kernel matrix spectra and the two-dimensional matrix spectrum into a preset composite spectrum calculation formula, obtain and based on the and of every two depth layers, obtain G T u, where G T u is an intermediate parameter in the three-dimensional gravity inversion process, k is the k-th depth layer of the three-dimensional model, l is the l-th depth layer of the three-dimensional model, and j 2 = -1.

8. The three-dimensional gravity rapid inversion optimization system according to claim 7, characterized in that, the preset composite spectrum formula is: Among them, is the forward kernel matrix of the k-th depth layer in the three-dimensional model, is the forward kernel matrix of the l-th depth layer in the three-dimensional model, (1 ≤ k, l ≤ p), where p is the total number of depth layers in the three-dimensional model, F represents the forward transform of the fast Fourier transform, F -1 represents the inverse transform of the fast Fourier transform, u is the single-layer model, j 2 = -1, m is the number of grid partitions in the x (north) direction of the three-dimensional model, n is the number of grid partitions in the y (east) direction of the three-dimensional model, and H is the composite frequency spectrum signal.

9. A storage medium, characterized in that, instructions are stored in the storage medium, and when the computer reads the instructions, the computer executes the three-dimensional gravity rapid inversion optimization method according to any one of claims 1 to 6.

10. An electronic device, comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that, when the processor executes the computer program, the computer executes the three-dimensional gravity rapid inversion optimization method according to any one of claims 1 to 6.

Images