Three-dimensional gravity and gravity tensor fast forward optimization method and system

Abstract

This invention discloses a fast forward modeling optimization method for three-dimensional gravity and gravity tensor, and a further optimization method for the system. The fast forward modeling method divides the underground three-dimensional space into multiple horizontal layered media and constructs a discrete model of the underground space and observation points. Based on the model division, the forward modeling coefficient matrix and coefficient spectrum matrix of gravity and gravity tensor are calculated. Based on the coefficient spectrum matrix, a fast algorithm is used to calculate the gravity value and gravity tensor value of each depth layer corresponding to the density model. The forward modeling values ​​of each layer are superimposed to obtain the gravity forward modeling anomaly of the density model. Furthermore, based on the parity characteristics of the fast Fourier transform of real number sequences, the spatial complexity of the forward modeling coefficient matrix is ​​reduced while ensuring fast calculation efficiency, providing support for further increasing the practicality of three-dimensional gravity inversion.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration, and in particular relates to a fast forward modeling optimization method and system for three-dimensional gravity and gravity tensor. Background Technology

[0002] Seeking a fast and efficient constraint inversion algorithm is crucial for the practical application of gravity three-dimensional inversion methods. Conventional gravity three-dimensional inversion algorithms require a large amount of computation time and storage space. To address this issue, a newly proposed fast gravity three-dimensional inversion algorithm has achieved exponential acceleration without sacrificing computational accuracy. This method is a combinatorial optimization algorithm: (1) it utilizes the precision of the spatial domain and the speed of the frequency domain to achieve high-precision and fast three-dimensional forward modeling; (2) it utilizes the symmetry in the inversion matrix to decompose the main calculations in gravity conjugate gradient inversion (CG) into two forward modeling calculations, thereby optimizing the computational efficiency of three-dimensional gravity inversion.

[0003] The algorithm described above can significantly reduce the memory footprint of the forward coefficient matrix. However, when the model is finely divided, the spectral matrix of the forward coefficient matrix will still occupy a large portion of the memory. For example, when the underground space is divided into 1200×1200×300 grids, the model matrix will occupy 3.45GB of memory, while the spectral matrix of the forward coefficient matrix, being a complex matrix four times the size of the model matrix, will occupy 27.65GB of memory. If the model is further expanded by 10% in all directions, the grid of the underground space becomes (1200+120*2)×(1200+120*2)×300. At this point, the model matrix will occupy 4.98GB of memory, and the spectral matrix of the forward coefficients will occupy 39.81GB of memory, making computation impossible on a typical mainstream computer with only 64GB of memory. Therefore, further optimization of the memory footprint of the forward coefficient matrix is ​​of great significance for the practical application of gravity forward and inverse modeling. Summary of the Invention

[0004] In view of this, this invention proposes an optimization algorithm that can further reduce the space complexity of the forward coefficient matrix. Based on the parity characteristics of the Fast Fourier Transform of real even sequences, this algorithm compresses the space complexity of the forward coefficient spectrum matrix to one-eighth of its original size without affecting computational efficiency. By reducing the space complexity of the forward coefficient matrix, the practicality of three-dimensional gravity inversion is further increased.

[0005] In a first aspect, the present invention provides a fast forward modeling optimization method for three-dimensional gravity and gravity tensor, comprising the following steps:

[0006] The underground three-dimensional space is divided into multiple horizontal layered media and a discrete model of the underground space and observation points is constructed. The discrete model is then divided and the gravity and gravity tensor forward modeling coefficient matrices are calculated.

[0007] Convert the forward coefficient matrix into a forward coefficient frequency spectrum matrix through the fast Fourier transform;

[0008] Classify and compressively store the forward coefficient frequency spectrum matrix;

[0009] Calculate the model forward field based on the compressed forward coefficient frequency spectrum matrix to complete the three-dimensional gravity fast forward calculation.

[0010] As an optimization of the above solution, divide the underground three-dimensional space into multiple horizontal layered media and construct a discrete model of the underground space and observation points, including:

[0011] Construct a Cartesian coordinate system, with the x-axis and y-axis as the horizontal directions and the z-axis as the vertically downward direction, corresponding to the eastward, northward, and underground depth directions respectively. Let (ξ, η, ζ) be the coordinates of any volume element dv = dζdηdζ in the anomaly body. Then the expression of the gravitational potential formula dV of the mass unit dm = ρ(ξ, η, ζ)dv at any point (x, y, z) in space is:

[0012]

[0013] where γ is the gravitational constant, and its value is 6.67×10 -11 m 3 / (kg·s 2 ); r = [(x - ζ) 2 +(y - η) 2 +(z - ζ) 2 1 / 2 , and r is the distance from the mass unit to any point (x, y, z) in space;

[0014] Integrate the formula (2-1) in the underground half-space according to the prism to obtain the gravitational potential expression V(x, y, z):

[0015]

[0016] where ρ(ξ, η, ζ) represents the spatial density distribution of the underground half-space; a and b are respectively half of the lengths of the single-layer medium in the x-direction and y-direction; L is the thickness of the underground half-space; and H represents the top surface depth.

[0017] As an optimization of the above solution, in the spatial domain, divide the underground three-dimensional space into multiple horizontal layered media. After vertically differentiating the expression (2-2), the gravity forward expression g(x, y, z) of a single-layer medium with a thickness of l (l < L) at the horizontal height z can be obtained:

[0018]

[0019] ​Where ρ(ξ,η) represents the lateral density distribution of the monolayer medium; a and b are half the length of the monolayer medium along the x and y directions, respectively. h represents the top surface depth, and l represents the thickness of the monolayer medium. For a monolayer medium, when the density function ρ is kept constant along the longitudinal direction, and h and l are fixed, G(x-ξ,y-η) is a function of x, ξ,y, and η, and the density function ρ(ξ,η) is a function of ξ and η.

[0020] Formula (2-1) can be regarded as a two-dimensional convolution of two signals. The expression for signal convolution is g = G * ρ', where ρ' is the density distribution matrix of the plate-like body in the horizontal space, i.e., the model signal; G can be regarded as the gravity forward modeling kernel matrix of a vertical line body with depth h, thickness l, and density 1. "*" indicates convolution, and g indicates forward modeling matrix.

[0021] As a preferred embodiment of the above scheme, the forward coefficient matrix is ​​transformed into a forward coefficient spectrum matrix through a fast Fourier transform, including:

[0022] Taking the forward modeling of a one-dimensional prism model as an example, the forward modeling calculation is performed by calculating the forward modeling values ​​of the corresponding observation points from a set of prisms, i.e., m prisms. The forward modeling calculation process is as follows:

[0023] Calculate the first forward kernel matrix G, and G1 to G2 in the first forward kernel matrix. m When the density value ρ of the m-th prism m For unit density, the forward gravity values ​​at each observation point, and G in the matrix. m+1 ~G 2m-1 and G1~G m It's about G m They are mutually symmetrical, i.e., G1 to G2. 2m-1 (1≤i≤m);

[0024] Perform circular convolution calculations on the model matrix. Zero-padding is applied to expand the matrix to the size of the forward kernel matrix. The convolution calculation process involves first inverting the model matrix, then multiplying it by the corresponding values ​​of the forward coefficient matrix and summing the results to obtain g. 2m-1 Then, the model matrix is ​​shifted one position to the left and multiplied by the corresponding forward coefficient matrix, and the results are summed to obtain g. 2m-2 Finally, repeat the above steps to complete one iteration of the calculation and obtain the forward anomaly matrix g. The g in matrix g... m ~g 2m-1 For the forward modeling result, g1~g in the matrix m-1 This is a useless calculation.

[0025] As a preferred embodiment of the above scheme, the forward coefficient spectrum matrix is ​​classified and compressed for storage, including:

[0026] The underground space is pre-divided into M prism meshes, where m and n represent the number of meshes in the x and y directions, respectively, and p represents the number of meshes in the z direction. Surface gravity observations are conducted from N observation points, where M = m × n × p and N = m × n. The forward modeling process of the three-dimensional prism model can be decomposed as follows:

[0027] Calculate the second forward kernel matrix for each depth layer. And convert it to a spectrum through Fast Fourier Transform. Post-storage;

[0028] Model matrix for each depth layer After extending the matrix to size (2m-1)×2n-1 by adding zeros, the spectrum is calculated by fast Fourier transform. And the spectrum of the second forward kernel matrix corresponding to the depth layer. Perform product operations;

[0029] Summing the product results of each depth layer and then performing an inverse Fourier transform to the spatial domain yields the fast forward anomaly formula:

[0030]

[0031] In the formula: d represents the model forward anomaly, F represents the Fast Fourier Transform, F -1 The symbol represents the inverse fast Fourier transform, and "·" represents the Hadamard product, which is the operation of multiplying corresponding elements in matrices of the same order.

[0032] As a preferred embodiment of the above scheme, the forward modeling field is calculated based on the compressed forward modeling coefficient spectrum matrix to complete the rapid forward modeling of three-dimensional gravity, including:

[0033]

[0034]

[0035] According to the Fourier transform, when x(t) is a real number sequence, its spectrum formula is:

[0036]

[0037] In the above formula, the real part and imaginary part of the spectrum of the real number sequence are respectively:

[0038]

[0039] When the real sequence x(t) is an even function, the real part of the spectrum F Re (ω) remains an even function, while the imaginary part of the spectrum F Im (ω) is then zero;

[0040] When the real sequence x(t) is an odd function, the real part of the spectrum F Re (ω) is zero, while the imaginary part of the spectrum F Im (ω) is an even function.

[0041] As a preferred embodiment of the above scheme, the two-dimensional matrix of the k-th depth layer is expressed as follows, based on the forward kernel matrix G corresponding to gravity and the gravity tensor:

[0042]

[0043] The elements of this two-dimensional matrix are related to G. 0,0 The rows or columns in which they are located are symmetrically distributed as even functions or symmetrically distributed as odd functions;

[0044] In actual spectrum calculation, matrix G needs to be... k To fill in zeros, we get...

[0045] matrix After Fast Fourier Transform, the spectrum matrix is ​​obtained.

[0046] As a preferred embodiment of the above scheme, based on the parity of the matrix, the forward coefficient matrix and its spectrum matrix corresponding to gravity and gravity tensor can be obtained. Divided into three categories:

[0047] The first type consists of coefficient matrices that are even functions along both the x and y directions, and the coefficient matrix includes V. xx V yy V zz The corresponding coefficient spectrum matrix The imaginary part of the coefficient spectrum matrix is ​​zero, and the real part of the coefficient spectrum matrix is ​​still an even function. The real part of the coefficient spectrum matrix is ​​related to... When the rows or columns containing the elements are symmetrically distributed, the coefficient spectrum matrix is... This can be expressed as:

[0048]

[0049] Compressed storage as At this time the matrix For real numbers:

[0050]

[0051] Where real() means taking the real part of the matrix;

[0052] The second type involves coefficient matrices where the function is even along one of the x or y directions and odd along the other. The coefficient matrix includes V.xz V yz The corresponding coefficient spectrum matrix The real part of the coefficient spectrum matrix is ​​zero, and the imaginary part of the coefficient spectrum matrix is ​​an even function along the even function direction and an odd function along the odd function direction.

[0053] Taking the x-direction as an odd function and the y-direction as an even function as an example, the coefficient spectrum matrix is ​​as follows: This can be expressed as:

[0054]

[0055] Compressed storage as

[0056]

[0057] In the formula: imag() represents taking the imaginary part of the matrix; when restoring this type of coefficient matrix, it needs to be set to an imaginary number;

[0058] The third type consists of coefficient matrices that are odd functions along both the x and y directions. These coefficient matrices include V. xy The corresponding coefficient spectrum matrix When the imaginary part is zero, and the real part of the coefficient spectrum matrix is ​​also an odd function in both the x and y directions, the coefficient spectrum matrix... This can be expressed as:

[0059]

[0060] Compressed storage as

[0061]

[0062] In the matrix The elements in the current row and column are all zero, and the elements in the first row and first column are also all zero.

[0063] As a preferred embodiment of the above scheme, the calculation process for gravity and the forward modeling of the gravity tensor includes:

[0064] Step 1: Calculate the second forward kernel matrix of the depth layer k to determine gravity and the gravity tensor G. k ((2m-1)×(2n-1)),(1≤k≤p), G k Add zero to extend Where p represents the number of vertical subdivisions of the model;

[0065] Step 2: Calculation Spectrum matrix Then compress and store the data to obtain a matrix.

[0066] Step 3: Expand each density matrix mk to 2m×2n by adding zeros:

[0067]

[0068] Then calculate the matrix. The spectrum of the real Fourier transform (FFT):

[0069]

[0070] In the formula: fft2 represents the two-dimensional forward Fourier transform;

[0071] Step 4: Through Recovering the spectrum of the forward coefficient matrix and the model spectrum matrix After multiplying and summing, g is obtained by inverse Fourier transform. ext :

[0072]

[0073] In the formula: ifft2 represents the two-dimensional inverse Fourier transform;

[0074] Step 5: Extract g ext The last m×n elements are used to obtain the gravity forward modeling result g of the density model.

[0075] Secondly, the present invention also provides a fast forward modeling optimization system for three-dimensional gravity and gravity tensor, comprising:

[0076] The model building module is used to divide the underground three-dimensional space into multiple horizontal layered media and construct a discrete model of the underground space and observation points. It also divides the discrete model and calculates the gravity and gravity tensor forward modeling coefficient matrix.

[0077] The forward modeling calculation module is used to transform the forward modeling coefficient matrix into a forward modeling coefficient spectrum matrix through a fast Fourier transform.

[0078] The compressed storage module is used to classify and compress the forward coefficient spectrum matrix.

[0079] The classification analysis module is used to calculate the forward modeling field based on the compressed forward modeling coefficient spectrum matrix to complete the rapid forward modeling of three-dimensional gravity.

[0080] This invention provides a fast forward modeling optimization method and system for three-dimensional gravity and gravity tensor. It involves dividing the underground three-dimensional space into multiple horizontal layered media and constructing a discrete model of the underground space and observation points. The discrete model is then divided and the forward modeling coefficient matrices of gravity and gravity tensor are calculated. These coefficient matrices are transformed into forward modeling coefficient spectrum matrices using a fast Fourier transform. The forward modeling coefficient spectrum matrices are then categorized and compressed for storage. Based on the compressed forward modeling coefficient spectrum matrix, the forward modeling field is calculated to complete the fast three-dimensional gravity forward modeling. Furthermore, an optimization algorithm for the space complexity of the forward modeling coefficient matrix can be developed. This algorithm, based on the parity characteristics of the fast Fourier transform of real number sequences, compresses the space complexity of the forward modeling coefficient spectrum matrix to one-eighth of its original size while ensuring fast computational efficiency. This reduces the space complexity of the forward modeling coefficient matrix and further enhances the practicality of three-dimensional gravity inversion. Attached Figure Description

[0081] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0082] Figure 1 This is a flowchart of the fast forward modeling optimization method for three-dimensional gravity and gravity tensor of the present invention;

[0083] Figure 2 This is a diagram showing the correspondence between the model unit mesh and the observation point mesh of the present invention;

[0084] Figure 3 This is a schematic diagram of the forward convolution calculation of the model in this invention;

[0085] Figure 4 This is a block diagram of the fast forward modeling optimization system for three-dimensional gravity and gravity tensor of the present invention. Detailed Implementation

[0086] Embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0087] See Figure 1 This invention provides a fast forward modeling optimization method for three-dimensional gravity and gravity tensor, comprising the following steps:

[0088] S1: Divide the underground three-dimensional space into multiple horizontal layered media and construct a discrete model of the underground space and observation points. Dissect the discrete model and calculate the forward coefficient matrices of gravity and gravity tensor.

[0089] S2: Convert the forward coefficient matrices into forward coefficient spectral matrices through fast Fourier transform.

[0090] S3: Classify and compressively store the forward coefficient spectral matrices.

[0091] S4: Calculate the forward field of the model based on the compressed forward coefficient spectral matrices to complete the three-dimensional gravity fast forward calculation.

[0092] In this embodiment, dividing the underground three-dimensional space into multiple horizontal layered media and constructing a discrete model of the underground space and observation points includes: constructing a Cartesian coordinate system, with the x-axis and y-axis as the horizontal directions and the z-axis as the vertically downward direction, corresponding to the eastward, northward, and underground depth directions respectively. Let (ξ, η, ζ) be the coordinates of any volume element dv = dξdηdζ in the abnormal body, then the expression of the gravitational potential formula dV of this mass element dm = ρ(ξ, η, ζ)dv at any point (x, y, z) in space is:

[0093]

[0094] where γ is the gravitational constant, and its value is 6.67×10 -11 m 3 / (kg·s 2 ); r = [(x - ξ) 2 +(y - η) 2 +(z - ζ) 2 1 / 2 , and r is the distance from the mass element (ξ, η, ζ) to any point (x, y, z) in space;

[0095] Integrate the formula (2-1) over the underground half-space according to the prism to obtain the gravitational potential expression V(x, y, z):

[0096]

[0097] where ρ(ξ, η, ζ) represents the spatial density distribution of the underground half-space; a and b are respectively half of the lengths of the single-layer medium along the x-direction and y-direction; L is the thickness of the underground half-space; H represents the depth of the top surface.

[0098] It should be noted that in the spatial domain, when the underground three-dimensional space is divided into multiple horizontal layered media and the expression (2-2) is vertically differentiated, the forward gravity expression g(x, y, z) of a single-layer medium with a thickness of l (l < L) at the horizontal height z can be obtained:

[0099] ​

[0100] Where ρ(ξ,η) represents the lateral density distribution of the monolayer medium; a and b are half the length of the monolayer medium along the x and y directions, respectively. h represents the top surface depth, and l represents the thickness of the monolayer medium. For a monolayer medium, when the density function ρ is kept constant along the longitudinal direction, and h and l are fixed, G(x-ξ,y-η) is a function of x, ξ,y, and η, and the density function ρ(ξ,η) is a function of ξ and η.

[0101] Formula (2-1) can be regarded as a two-dimensional convolution of two signals. The expression for signal convolution is g = G * ρ', where ρ' is the density distribution matrix of the plate-like body in the horizontal space, i.e., the model signal; G can be regarded as the gravity forward modeling kernel matrix of a vertical line body with depth h, thickness l, and density 1. "*" indicates convolution, and g indicates forward modeling matrix.

[0102] It should be understood that in the spatial domain, to illustrate the relationship between forward modeling and convolution, the underground three-dimensional space is divided into multiple horizontal layered media. In actual calculations, both the underground space and the observation points are discrete, with a one-to-one correspondence between gravity observation points and the horizontal positions of the prism. Figure 2 As shown in the figure, this ensures that the model resolution is high enough and facilitates the design of targeted computational strategies. Convolution can be divided into linear convolution and circular convolution. Only circular convolution can use the Fast Fourier Transform (FFT) algorithm to improve computational speed. Further optimization algorithms can be developed to reduce the space complexity of the forward coefficient matrix. Based on the parity characteristics of the FFT for real even sequences, this algorithm compresses the space complexity of the forward coefficient spectrum matrix to one-eighth of its original size without affecting computational efficiency. By reducing the space complexity of the forward coefficient matrix, the practicality of three-dimensional gravity inversion is further increased.

[0103] Optionally, the forward coefficient matrix is ​​transformed into a forward coefficient spectrum matrix via a fast Fourier transform, including:

[0104] Taking the forward modeling of a one-dimensional prism model as an example, the forward modeling calculation is performed by calculating the forward modeling values ​​(m forward values) of the corresponding observation points from a set of prisms (m prisms). The forward modeling calculation process is as follows:

[0105] Calculate the first forward kernel matrix G, and G1 to G2 in the first forward kernel matrix. m When the density value ρ of the m-th prism m For unit density, the forward gravity values ​​at each observation point, and G in the matrix. m+1 ~G 2m-1 and G1~G m It's about G m They are mutually symmetrical, i.e., G1 to G2.2m-1 (1≤i≤m);

[0106] Perform circular convolution calculations on the model matrix. Zero-padding is applied to expand the matrix to the size of the forward kernel matrix. The convolution calculation process involves first inverting the model matrix, then multiplying it by the corresponding values ​​of the forward coefficient matrix and summing the results to obtain g. 2m-1 Then, the model matrix is ​​shifted one position to the left and multiplied by the corresponding forward coefficient matrix, and the results are summed to obtain g. 2m-2 Finally, repeat the above steps to complete one iteration of the calculation and obtain the forward anomaly matrix g. In matrix g, g... m ~g 2m-1 For the forward modeling result, g1~g in the matrix m-1 This is a useless calculation.

[0107] In this embodiment, see Figure 3 To illustrate the relationship between forward modeling and circular convolution in a single-layer model, we'll take the forward modeling of a one-dimensional prism model as an example. This forward modeling calculation involves using a set of prisms (m prisms) to obtain the forward modeling values ​​(m forward values) for the corresponding observation points. The forward modeling process can be divided into two steps:

[0108] Step 1: Calculate the first forward kernel matrix G( Figure 3 (a) in the matrix, G1~G m When the density ρ of the m-th prism m For unit density, the forward gravity values ​​at each observation point, and G in the matrix. m+1 ~G 2m-1 and G1~G m It's about G m They are mutually symmetrical, i.e., G1 to G2. 2m-1 (1≤i≤m).

[0109] Step 2: Perform circular convolution calculation. Circular convolution requires the two matrices to have the same length; therefore, the model matrix needs to be... Zero-padding expands the matrix to the size of the forward kernel. The convolution calculation process begins by inverting the model matrix. Figure 3 After part (b) in the formula, multiply the corresponding parts with the forward coefficient matrix and sum them to obtain g. 2m-1 Then, the model matrix is ​​shifted one position to the left in a circular motion. Figure 3 The (c) part of the matrix is ​​multiplied by the corresponding forward coefficient matrix and summed to obtain g. 2m-2 Finally, repeat the above steps. Figure 3 Part (b) of the middle Figure 3 The (f) part of the algorithm completes one iteration of the calculation to obtain the forward anomaly matrix g. The matrix g contains g... m ~g 2m-1 For the forward modeling result, g1~g in the matrix m-1In this invention, useless calculations are performed ( Figure 3 (parts (d) to (e) in the text).

[0110] Optionally, the forward coefficient spectrum matrix is ​​classified and compressed for storage, including:

[0111] The underground space is presumably divided into M (M = m × n × p) prism meshes, where m and n represent the number of meshes in the x (north) and y (east) directions, respectively, and p represents the number of meshes in the z (vertical) direction. The surface gravity observations are conducted from N (N = m × n) observation points. The forward modeling process of the three-dimensional prism model can be decomposed as follows:

[0112] Calculate the second forward kernel matrix for each depth layer. And convert it to a spectrum through Fast Fourier Transform. Post-storage;

[0113] Model matrix for each depth layer After extending the matrix to size (2m-1)×2n-1 by adding zeros, the spectrum is calculated by fast Fourier transform. And the spectrum of the second forward kernel matrix corresponding to the depth layer. Perform product operations;

[0114] Summing the product results of each depth layer and then performing an inverse Fourier transform to the spatial domain yields the fast forward anomaly formula:

[0115]

[0116] In the formula: d represents the model forward anomaly, F represents the Fast Fourier Transform, F -1 The symbol "☐" represents the inverse fast Fourier transform; "·" represents the Hadamard product, which is the operation of multiplying corresponding elements of matrices of the same order.

[0117] In this embodiment, when the structure of the model is determined, the spectrum of the forward kernel matrix is... This is also certain; therefore, for cases requiring multiple forward modeling calculations, only one calculation is needed. Then you can save it. Because... The matrix is ​​obtained from the spatial domain kernel matrix through FFT, therefore the theoretical forward modeling accuracy of this algorithm is the same as that of the spatial domain forward modeling. In actual calculations, there is no spectral truncation error, only numerical truncation error, which will be proven in numerical simulations.

[0118] It should be noted that when performing forward modeling only in the spatial domain, with the first forward kernel matrix G already calculated, the time complexity of the forward modeling is O(M×N). However, when using the fast forward modeling formula (Equation 2-4), the time complexity is O(M×N) after the first forward kernel matrix G has been calculated and transformed to the frequency domain. In this case, it is necessary to Performing p forward Fourier transforms with a time complexity of O(4N×log2(4N)) and one inverse Fourier transform on the product with a time complexity of O(4N×log2(4N)), when p>>1, the time complexity of the frequency domain forward modeling can be approximated as O(4M×log2(4N)). Therefore, the theoretical speedup of the fast algorithm is approximately N / (4log2(4N)). The speedup effect of this algorithm is similar to the relationship between the fast Fourier transform and the Fourier transform; therefore, the speedup increases with the increase of the observed data N.

[0119] Optionally, the forward modeling field is calculated based on the compressed forward modeling coefficient spectrum matrix to complete the three-dimensional gravity fast forward modeling, including:

[0120]

[0121] According to the Fourier transform, when x(t) is a real number sequence, its spectrum formula is:

[0122]

[0123] In the above formula, the real part and imaginary part of the spectrum of the real number sequence are respectively:

[0124]

[0125] When the real sequence x(t) is an even function, the real part of the spectrum F Re (ω) remains an even function, while the imaginary part of the spectrum F Im (ω) is then zero;

[0126] When the real sequence x(t) is an odd function, the real part of the spectrum F Re (ω) is zero, while the imaginary part of the spectrum F Im (ω) is an even function.

[0127] In this embodiment, the two-dimensional matrix of the k-th depth layer is represented by the first forward kernel matrix G corresponding to gravity and the gravity tensor as follows:

[0128]

[0129] The elements of this two-dimensional matrix are related to G. 0,0 The rows or columns in which they are located are symmetrically distributed as even functions or symmetrically distributed as odd functions;

[0130] In actual spectrum calculation, matrix G needs to be... k To fill in zeros, we get...

[0131] matrix After Fast Fourier Transform, the spectrum matrix is ​​obtained.

[0132] It should be noted that, based on the parity of the matrix, the forward coefficient matrix and its spectrum matrix corresponding to gravity and the gravity tensor can be obtained. Divided into three categories:

[0133] The first type consists of coefficient matrices that are even functions along both the x and y directions, and the coefficient matrix includes V. xx V yy V zz The corresponding coefficient spectrum matrix The imaginary part of the coefficient spectrum matrix is ​​zero, and the real part of the coefficient spectrum matrix is ​​still an even function. The real part of the coefficient spectrum matrix is ​​related to... When the rows or columns containing the elements are symmetrically distributed, the coefficient spectrum matrix is... This can be expressed as:

[0134]

[0135] Compressed storage as At this time the matrix For real numbers:

[0136]

[0137] Where real() means taking the real part of the matrix;

[0138] The second type involves coefficient matrices where the function is even along one of the x or y directions and odd along the other. The coefficient matrix includes V. xz V yz The corresponding coefficient spectrum matrix The real part of the coefficient spectrum matrix is ​​zero, and the imaginary part of the coefficient spectrum matrix is ​​an even function along the even function direction and an odd function along the odd function direction.

[0139] Taking the x-direction as an odd function and the y-direction as an even function as an example, the coefficient spectrum matrix is ​​as follows: This can be expressed as:

[0140]

[0141] Compressed storage as

[0142]

[0143] In the formula: imag() represents taking the imaginary part of the matrix; when restoring this type of coefficient matrix, it needs to be set to an imaginary number;

[0144] The third type consists of coefficient matrices that are odd functions along both the x and y directions. These coefficient matrices include V. xy The corresponding coefficient spectrum matrix The imaginary part of the coefficient spectrum matrix is ​​zero, and the real part of the coefficient spectrum matrix is ​​also an odd function in both the x and y directions. This can be expressed as:

[0145]

[0146] Compressed storage as

[0147]

[0148] In the matrix The elements in the current row and column are all zero, and the elements in the first row and first column are also all zero.

[0149] The above analysis shows that the spectrum of gravity and the forward coefficient matrix of the gravity tensor contains a large number of zero values ​​and repeating values. In actual calculations, only a small number of these zero values ​​and repeating values ​​need to be extracted. The real or imaginary parts of the first (m+1)×(n+1) elements are stored as a real matrix. The original spectrum can then be completely recovered. And the real matrix... The memory occupied is only an imaginary matrix. It occupies one-eighth of the memory size, thus greatly saving the space complexity of the forward coefficient matrix.

[0150] Optionally, the calculation process for gravity and the forward modeling of the gravity tensor includes:

[0151] Step 1: Calculate the second forward kernel matrix of the depth layer k to determine gravity and the gravity tensor G. k ((2m-1)×(2n-1)),(1≤k≤p), G k Add zero to extend Where p represents the number of longitudinal subdivisions of the model, and those with "~" (i.e., tilde) represent the spectral matrix, while those without "~" (i.e., tilde) represent the real matrix.

[0152] Step 2: Calculation Spectrum matrix Then compress and store the data to obtain a matrix.

[0153] Step 3: Expand each density matrix mk to 2m×2n by adding zeros:

[0154]

[0155] Then calculate the matrix. The spectrum of the real Fourier transform (FFT):

[0156]

[0157] In the formula: fft2 represents the two-dimensional forward Fourier transform;

[0158] Step 4: Through Recovering the spectrum of the forward coefficient matrix and the model spectrum matrix After multiplying and summing, g is obtained by inverse Fourier transform. ext :

[0159]

[0160] In the formula: ifft2 represents the two-dimensional inverse Fourier transform;

[0161] Step 5: Extract g ext The last m×n elements are used to obtain the gravity forward modeling result g of the density model.

[0162] It should be noted that when a large number of forward calculations are required during the inversion calculation of gravity and gravity tensor, steps one and two only need to be calculated once, and steps three to five only need to be repeated.

[0163] See Figure 4 The present invention also provides a fast forward modeling optimization system for three-dimensional gravity and gravity tensor, comprising:

[0164] The model building module is used to divide the underground three-dimensional space into multiple horizontal layered media and construct a discrete model of the underground space and observation points. It also divides the discrete model and calculates the gravity and gravity tensor forward modeling coefficient matrix.

[0165] The forward modeling calculation module is used to transform the forward modeling coefficient matrix into a forward modeling coefficient spectrum matrix through a fast Fourier transform.

[0166] The compressed storage module is used to classify and compress the forward coefficient spectrum matrix.

[0167] The classification analysis module is used to calculate the forward modeling field based on the compressed forward modeling coefficient spectrum matrix to complete the rapid forward modeling of three-dimensional gravity.

[0168] In this embodiment, for the case where the subdivision grid of the underground space is (1200+120*2)×(1200+120*2)×300, the model matrix will occupy 4.98G of memory space. At this time, the spectrum matrix of the forward coefficients also only occupies about 4.98G of memory, so it can be calculated on a computer with 64G of memory. By dividing the underground three-dimensional space into multiple horizontal layered media and constructing a discrete model of the underground space and observation points, the gravity and gravity tensor forward modeling coefficient matrices are calculated. The forward modeling coefficient matrices are then transformed into forward modeling coefficient spectrum matrices through fast Fourier transform. The forward modeling coefficient spectrum matrices are classified and compressed for storage. Based on the compressed forward modeling coefficient spectrum matrix, the forward modeling field of the model is calculated to complete the fast forward modeling of three-dimensional gravity. Further optimization algorithms for the space complexity of the forward modeling coefficient matrix can be developed. Based on the parity characteristics of the fast Fourier transform of real number sequences, this algorithm compresses the space complexity of the forward modeling coefficient spectrum matrix to one-eighth of its original size while ensuring fast computational efficiency. This reduces the space complexity of the forward modeling coefficient matrix and further enhances the practicality of three-dimensional gravity inversion.

[0169] In all examples shown and described herein, any specific values ​​should be interpreted as merely exemplary and not as limitations; therefore, other examples of exemplary embodiments may have different values.

[0170] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0171] The above-described embodiments are merely illustrative of several implementations of the present invention, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of the invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention.

Claims

1. A fast forward modeling optimization method for three-dimensional gravity and gravity tensor, characterized in that, Includes the following steps: The underground three-dimensional space is divided into multiple horizontal layered media and a discrete model of the underground space and observation points is constructed. The discrete model is then divided and the gravity and gravity tensor forward modeling coefficient matrices are calculated. The forward coefficient matrix is ​​transformed into a forward coefficient spectrum matrix using a fast Fourier transform; The forward coefficient spectrum matrix is ​​classified and compressed for storage; including: pre-defined subsurface space subdivision. A prismatic mesh, where m and n represent the number of mesh subdivisions in the x and y directions, respectively, and p represents the number of mesh subdivisions in the z direction. The observed surface gravity is... There are 1 observation point, among which, , ; (2-5) (2-6) (2-7) (2-8) (2-9) (2-10) As can be seen from the Fourier transform, when When the sequence is a real number sequence, its spectrum formula is: （2-11） In the above formula, the real part and imaginary part of the spectrum of the real number sequence are respectively: (2-12) (2-13) When the real number sequence When it is an even function, the real part of the spectrum It remains an even function, but the imaginary part of the spectrum... Then it is zero; When the real number sequence When it is an odd function, the real part of the spectrum The value is zero, while the imaginary part of the spectrum is zero. It is then an even function; Based on gravity and the first forward kernel matrix corresponding to the gravity tensor. The two-dimensional matrix of the k-th depth layer is represented as: (2-14) The elements of this two-dimensional matrix are related to The rows or columns in which they are located are symmetrically distributed as even functions or symmetrically distributed as odd functions; In actual spectrum calculation, it is necessary to use the matrix To fill in zeros, we get... : (2-15) matrix After Fast Fourier Transform, the spectrum matrix is ​​obtained. (2-16) in, Represents a two-dimensional forward Fourier transform; Based on the parity of the matrix, the forward coefficient matrix and its spectrum matrix corresponding to gravity and gravity tensor can be obtained. : The first type consists of coefficient matrices that are even functions along both the x and y directions. The coefficient matrix includes... , , The corresponding coefficient spectrum matrix The imaginary part of the coefficient spectrum matrix is ​​zero, the real part of the coefficient spectrum matrix is ​​an even function, and the real part of the coefficient spectrum matrix is ​​related to the fact that the imaginary part of the coefficient spectrum matrix is ​​zero. When the rows or columns containing the elements are symmetrically distributed, the coefficient spectrum matrix is... (2-17) Compressed storage as At this time, the matrix For real numbers: in, This indicates taking the real part of the matrix; The second type involves coefficient matrices where the function is even along one of the x or y directions and odd along the other. The coefficient matrix includes... The corresponding coefficient spectrum matrix The real part of the coefficient spectrum matrix is ​​zero, and the imaginary part of the coefficient spectrum matrix is ​​an even function along the even function direction and an odd function along the odd function direction. When the function is odd along the x-direction and even along the y-direction, the coefficient spectrum matrix is... The expression is as follows: （2-18） Compressed storage as : In the formula: This indicates taking the imaginary part of the matrix; when restoring the coefficient matrix, it is set to an imaginary number. The third type consists of coefficient matrices that are odd functions along both the x and y directions. The coefficient matrix includes... The corresponding coefficient spectrum matrix The imaginary part of the coefficient spectrum matrix is ​​zero, and the real part of the coefficient spectrum matrix is ​​also an odd function in both the x and y directions. The expression is as follows: （2-19） Compressed storage as : The elements in the current row and column are all zero, and the elements in the first row and first column are also all zero; The underground three-dimensional space is divided into multiple horizontal layered media, and a discrete model of the underground space and observation points is constructed, including: Construct a Cartesian coordinate system, with the x and y axes representing the horizontal direction and the z axis representing the vertical downward direction, corresponding to the east, north, and underground depth directions, respectively. Let... Let be the coordinates of any volume element within the anomalous body. Then, the mass element at any point in space... Gravitational potential formula The expression is: (2-1) The expression for the arbitrary volume element within the abnormal body is: The expression for the mass element is: , It is the gravitational constant; r is the distance from a mass element to any point in space. The distance; Integrating equation (2-1) into the subsurface half-space of a prism yields the expression for the gravitational potential. : (2-2) in, This represents the spatial density distribution of the underground half-space; a and b are half the length of the single-layer medium along the x and y directions, respectively; L is the thickness of the underground half-space; H represents the depth of the top surface. In the spatial domain, the underground three-dimensional space is divided into multiple horizontally layered media. Taking the vertical derivative of expression (2-2), we can obtain a certain thickness of... Forward gravity modeling expression for a single-layer medium at horizontal height z : (2-3) in, , This represents the lateral density distribution of a monolayer medium; a and b are half the length of the monolayer medium along the x and y directions, respectively, resulting in... h represents the depth of the top surface. This represents the thickness of a single-layer medium, where the density function is constrained for a single-layer medium. Invariable along the longitudinal direction, and h and When fixed, yes The function, density function yes The function; Treating equation (2-3) as a two-dimensional convolution of two signals, the expression for signal convolution is: ,in Let G be the density distribution matrix of the plate-like body in the transverse space, i.e., the model signal; G is considered as having a depth of h and a thickness of... A gravity forward modeling kernel matrix for a vertical line body with density of 1, where "*" indicates convolution. Represents the orthogonal matrix; The forward model is calculated based on the compressed forward coefficient spectrum matrix to complete the rapid forward modeling of three-dimensional gravity.

2. The fast forward modeling optimization method for three-dimensional gravity and gravity tensor according to claim 1, characterized in that, The forward coefficient matrix is ​​transformed into a forward coefficient spectrum matrix using a fast Fourier transform, including: For forward modeling of a one-dimensional prism model, the forward modeling calculation is performed using a set of prisms, i.e. The forward modeling value for the corresponding observation point is obtained by calculating the value of each prism. There are several forward values, and the forward calculation process is as follows: Calculate the first forward kernel matrix In the first forward kernel matrix Is the first Density value of a prism Forward gravity values ​​at each observation point when the density is unit; in the matrix and It is about The corresponding forward values ​​of gravity are symmetrical. ; Perform circular convolution calculations on the model matrix. Zero-padding is applied to expand the matrix to the size of the forward kernel matrix. The convolution calculation process involves first inverting the model matrix, then multiplying it by the corresponding values ​​of the forward coefficient matrix and summing the results. Then, the model matrix is ​​shifted one position to the left and multiplied by the corresponding forward coefficient matrix, and the results are summed to obtain the model matrix. Finally, repeat the above steps to complete one iteration of the calculation and obtain the forward anomaly matrix. ,matrix In For the forward modeling result, the matrix contains This is a useless calculation.

3. The fast forward modeling optimization method for three-dimensional gravity and gravity tensor according to claim 1, characterized in that, The forward coefficient spectrum matrix is ​​classified and compressed for storage, including: Pre-defined underground space division A prismatic mesh, where m and n represent the number of mesh subdivisions in the x and y directions, respectively, and p represents the number of mesh subdivisions in the z direction. The observed surface gravity is... There are 1 observation point, among which, , The forward modeling process of the three-dimensional prism model can be decomposed into: Calculate the second forward kernel matrix for each depth layer. And converted into a spectrum through Fast Fourier Transform. Post-storage; Model matrix for each depth layer Expanded by adding zeros respectively After obtaining the matrix of size, the spectrum is then calculated using a Fast Fourier Transform. And the spectrum of the second forward kernel matrix corresponding to the depth layer. Perform product operations; Summing the product results of each depth layer and then performing an inverse Fourier transform to the spatial domain yields the fast forward anomaly formula: (2-4) In the formula: This indicates a forward modeling anomaly, where F represents the forward Fast Fourier Transform. Indicates the inverse fast Fourier transform; "" represents the Hadamard product, which is the operation of multiplying corresponding elements in matrices of the same order.

4. The fast forward modeling optimization method for three-dimensional gravity and gravity tensor according to claim 1, characterized in that, The calculation process for gravity and the forward modeling of the gravity tensor includes: Step 1: Calculation The second forward kernel matrix of the depth layer is used to determine the gravity and gravity tensor, where the expressions for gravity and gravity tensor are: ,Will Add zero to extend , This indicates the number of vertical subdivisions in the model. This represents a two-dimensional matrix representing the k-th depth layer; Step 2: Calculation Spectrum matrix Then compress and store the data to obtain a matrix. ; Step 3: Assemble the density layer matrices Add zeros to extend to : （2-20） Then calculate the matrix. The spectrum of the real Fourier transform (FFT): （2-21） In the formula: fft2 represents the two-dimensional forward Fourier transform; Step 4: Through Recovering the spectrum of the forward coefficient matrix and with After multiplying and summing, the result is obtained by inverse Fourier transform. : （2-22） In the formula: ifft2 represents the two-dimensional inverse Fourier transform; Step 5: Extraction After The elements are used to obtain the gravity forward modeling results of the density model. .

5. A fast forward modeling optimization system for three-dimensional gravity and gravity tensor, applied to the fast forward modeling optimization method for three-dimensional gravity and gravity tensor as described in any one of claims 1-4, characterized in that, include: The model building module is used to divide the underground three-dimensional space into multiple horizontal layered media and construct a discrete model of the underground space and observation points. It also divides the discrete model and calculates the gravity and gravity tensor forward modeling coefficient matrix. The forward modeling calculation module is used to transform the forward modeling coefficient matrix into a forward modeling coefficient spectrum matrix through a fast Fourier transform. The compressed storage module is used to classify and compress the forward coefficient spectrum matrix. The classification analysis module is used to calculate the forward modeling field based on the compressed forward modeling coefficient spectrum matrix to complete the rapid forward modeling of three-dimensional gravity.

Images