CN115659579A - Three-dimensional gravity field numerical simulation method and system based on 3D AS-FT - Google Patents
Three-dimensional gravity field numerical simulation method and system based on 3D AS-FT Download PDFInfo
- Publication number
- CN115659579A CN115659579A CN202211076505.7A CN202211076505A CN115659579A CN 115659579 A CN115659579 A CN 115659579A CN 202211076505 A CN202211076505 A CN 202211076505A CN 115659579 A CN115659579 A CN 115659579A
- Authority
- CN
- China
- Prior art keywords
- dimensional
- domain
- wave number
- fourier transform
- gravity field
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 230000005484 gravity Effects 0.000 title claims abstract description 81
- 238000000034 method Methods 0.000 title claims abstract description 59
- 238000004088 simulation Methods 0.000 title claims abstract description 35
- 238000005070 sampling Methods 0.000 claims abstract description 57
- 238000001228 spectrum Methods 0.000 claims abstract description 22
- 230000006870 function Effects 0.000 claims description 19
- 230000009466 transformation Effects 0.000 claims description 15
- 230000002159 abnormal effect Effects 0.000 claims description 14
- 238000012887 quadratic function Methods 0.000 claims description 6
- 238000004590 computer program Methods 0.000 claims description 4
- 230000001131 transforming effect Effects 0.000 claims description 4
- 230000014509 gene expression Effects 0.000 claims description 3
- 230000010354 integration Effects 0.000 claims description 3
- 230000003595 spectral effect Effects 0.000 abstract description 11
- 238000004422 calculation algorithm Methods 0.000 abstract description 5
- 238000004364 calculation method Methods 0.000 description 13
- 230000008859 change Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 238000003384 imaging method Methods 0.000 description 3
- 230000000737 periodic effect Effects 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000004611 spectroscopical analysis Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 238000010561 standard procedure Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 239000013598 vector Substances 0.000 description 1
- 239000002699 waste material Substances 0.000 description 1
Images
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Complex Calculations (AREA)
Abstract
The invention discloses a three-dimensional gravity field numerical simulation method and a three-dimensional gravity field numerical simulation system based on 3D random sampling Fourier transform (AS-FT), which are characterized in that Fourier transform of three-dimensional random sampling is used for carrying out Fourier transform on a Poisson equation meeting the gravity potential, the gravity potential is directly obtained in a spectrum domain, the solution is simple, the spectrum domain gravity field is obtained by using the relation between the spectrum domain gravity potential and the gravity field, and the space domain gravity field can be obtained by carrying out three-dimensional Fourier inverse transform on the gravity field. The algorithm has good parallelism and occupies less memory. Therefore, the spectral domain gravity numerical simulation that both the spatial domain and the wavenumber domain can be sampled at will is realized, and the method has good guiding significance for the condition that non-uniform sampling is required. By means of Fourier transform capable of sampling at will, the algorithm has high precision and efficiency.
Description
Technical Field
The invention relates to the field of gravity field numerical simulation, in particular to a three-dimensional gravity field numerical simulation method and system based on 3D AS-FT.
Background
The method is used as a means for solving differential equations, can directly convert partial derivatives into products, and has the advantages of high precision, high convergence speed and the like. The spectroscopic method is essentially a generalization of the standard technique of separating variables. Generally, chebyshev polynomials and Legendre polynomials are taken as basis functions for approximating the expansion. For an equation satisfying periodic boundary conditions, fourier series and harmonic series are simple and convenient. The spectroscopic methods are classified as Fourier methods, chebyshev or Legendre methods. The former applies to periodic problems, while the latter two apply to non-periodic problems. The basis of the methods is to establish a spatial basis function, and the precision of the method directly depends on the number of expansion terms, so that a plurality of terms are needed to achieve better precision, thereby causing large calculation amount and low efficiency.
The Fourier spectrum method is adopted, and the calculation efficiency of the Fourier spectrum method depends on the selected Fourier transform method. However, the accuracy of fourier forward modeling is relatively low due to inherent drawbacks when using Discrete Fourier Transforms (DFT) instead of continuous fourier transforms, such as aliasing, edge effects, added periodicity, and truncation effects. If the standard FFT is adopted, the edge expanding treatment is needed to weaken the influence of the boundary effect, but the calculation scale is increased; if Gauss-FFT is used for numerical simulation, the efficiency is greatly reduced. And the standard FFT or the Gauss-FFT can only be used for uniform sampling, thereby causing the waste of computing power.
The method for performing gravity field numerical simulation in a spectral domain by applying Fourier transform has related application. The method comprises the following steps that documents (Three-dimensional numerical modeling of gravity and magnetic analysis in a mixed space-wave number domain, geophilcs, 2018, shikun Dai, dongdong Zhao, et al) perform numerical simulation on a gravity field by respectively using a standard FFT method and a Gauss-FFT method from a gravity differential equation, and documents (Three-dimensional numerical simulation of a gravity abnormal field space-wave number mixed domain, geophysical report, 2020, daiky, chen California and the like) perform numerical simulation on the gravity field by respectively using the standard FFT method and the Gauss-FFT method from an integral equation; however, it is pointed out that the standard FFT must be edge-extended to achieve the same accuracy as the Gauss-FFT, but the Gauss-FFT has relatively low efficiency although it has high accuracy and small boundary effect, and the standard FFT and Gauss-FFT methods are both non-uniform sampling and inflexible. In the literature (High-accuracy 3D Fourier transformer of gradient field b filed on the Gauss-FFT technique, journal of Applied geomatics, 2018, guangdong ZHao, bo Chen, et al), 3D Gauss-FFT is adopted for numerical simulation, and the method realizes the High-precision numerical simulation of the gravity three-dimensional spectral domain, but the Gauss-FFT occupies a higher memory and can only carry out non-uniform sampling.
Disclosure of Invention
The invention provides a three-dimensional gravity field numerical simulation method and system based on 3D AS-FT, which are used for solving the technical problem that the efficiency and the accuracy of the conventional gravity field numerical simulation method are incompatible.
In order to solve the technical problems, the technical scheme provided by the invention is as follows:
a three-dimensional gravity field numerical simulation method based on 3D AS-FT comprises the following steps:
constructing a three-dimensional space numerical model of the target area based on the actual topography of the target area and the actual distribution condition of the abnormal body; assigning residual density to each node in the three-dimensional space numerical model, and constructing a three-dimensional Poisson equation of the target area based on the three-dimensional space numerical model after assignment;
transforming the three-dimensional Poisson equation to a wave number domain by adopting three-dimensional Fourier transform of arbitrary sampling, solving the wave number domain gravity potential of the target area, and solving the wave number domain gravity field of the target area according to the relationship between the wave number domain gravity potential and the wave number domain gravity field;
and performing randomly sampled three-dimensional Fourier inverse transformation on the wave number domain gravity field of the target area to obtain a space domain gravity field of the target area.
Preferably, the arbitrarily sampled three-dimensional fourier positive transform is:
carrying out x-direction one-dimensional Fourier forward transformation on the three-dimensional Poisson equation;
wherein x, y, z represent three mutually perpendicular directions; k is a radical of x Representing the wave number in the x direction, F is a spatial domain function, and F is a wave number spectrum;
to F (k) x Y, z) performing a one-dimensional fourier transform in the y direction:
to F xy (k x ,k y Z) performing a z-direction one-dimensional fourier transform:
preferably, the arbitrarily sampled three-dimensional fourier transform is a three-dimensional arbitrary fourier transform based on shape function interpolation.
Preferably, the one-dimensional fourier transform specifically includes:
let the continuous one-dimensional fourier transform be respectively expressed as:
the continuous one-dimensional Fourier transform integral is dispersed to obtain:
wherein M represents the number of cells, e j Represents the jth cell, where i is an imaginary number;
interpolating f (x) with a quadratic function:
when the quadratic interpolation shape function fitting is adopted in the unit, the coordinates of three nodes in any unit are setIs other than x 1 、x 2 、x 3 ，x 2 Is a midpoint, satisfies x 1 +x 3 ＝2x 2 The value at each node is f (x) 1 )、f(x 2 )、f(x 3 ) F (x) is expressed by a quadratic function, which can be expressed as:
f(x)＝N 1 f(x 1 )+N 2 f(x 2 )+N 3 f(x 3 )
wherein,
the above formula can be written as:
order to
For the in-cell fourier transform node coefficients, the above equation is abbreviated as:
when wave number k x When not 0, substituting the above formula into W 1 、W 2 、W 3 In the above, the in-cell fourier transform node coefficients are obtained:
W 1 、W 2 、W 3 the integral kernel functions all comprise
Which is in [ x ] 1 ,x 3 ]The upper unit integral is resolved into:
thus, k can be obtained x W when not 0 1 、W 2 、W 3 The semi-analytic solution is:
when wave number k x When the value is 0, the number of the first electrode is,
the Fourier transform node coefficient under zero wave number can be obtained by simple integration as follows:
and accumulating the analytical expressions of different units to obtain a final one-dimensional Fourier forward transformation result.
Preferably, the wave number domain gravitational potential is obtained through three-dimensional Fourier transform of arbitrary sampling
Satisfies the equation:
wherein,
for the wave number domain residual density distribution, the number domain gravitational potential can be obtained by solving according to the formula
Preferably, the wave number domain gravitational potential
Sum wave number domain gravity field
The relationship of (1) is:
according to the formula, the wave number domain gravity field can be obtained by solving
In milli-gal.
Preferably, the three-dimensional arbitrary sampling inverse fourier transform formula is:
in the formula k x 、k y 、k z Representing the wavenumber, F (x, y, z) is a function of the spatial domain, F (k) x ,k y ,k z ) Representing a wavenumber spectrum.
Preferably, the three-dimensional poisson equation is:
wherein, U g Is a space domain gravitational potential, pi is a circumferential rate, G is a universal gravitation constant, and rho is a space domain residual density with the unit of kg/m 3 。
Preferably, the space domain of the target model is subdivided in any one of the following ways to construct a three-dimensional numerical model of the target region,
the first method is as follows: non-uniform subdivision is carried out on a preset first area, and encryption is carried out; wherein the first region satisfies the following formula:
ρ i is the residual density, rho, of the corresponding node of the first region j The residual density of the jth node around the corresponding node of the first area is shown, and n is the number of the nodes around the first area; omega is weight and the value range is (0, 1).
The second method comprises the following steps: sparse sampling is carried out on a preset second area, wherein the second area meets the following formula:
ρ i′ the residual density, rho, of the corresponding node of the second region j′ The residual density of the jth node around the corresponding node of the second area is obtained, and n is the number of the nodes around the first area; omega is weight and the value range is (0, 1);
the third method comprises the following steps: the spatial domain samples uniformly in all three directions.
A computer system comprising a memory, a processor and a computer program stored on the memory and executable on the processor, the processor implementing the steps of the method when executing the computer program.
The invention has the following beneficial effects:
the method applies a three-dimensional random sampling Fourier transform method (3 DAS-FT), transforms a three-dimensional gravity differential equation into a spectrum domain to be directly solved, and obtains a solution of a space domain through three-dimensional random sampling Fourier inverse transform. By means of high precision and high efficiency of the 3D AS-FT, the Fourier spectrum method of the three-dimensional gravity field can achieve high precision, has high calculation efficiency, can sample at will, and has important significance for accurate imaging of gravity three-dimensional inversion.
In addition, the spatial domain is transformed to the wave number domain by using three-dimensional random sampling Fourier forward transform, so that the x, y and z directions of the spatial domain can be non-uniformly subdivided. The abnormal body can be better fitted, and the calculation precision is improved. And sparse sampling is performed in the area which does not need to be encrypted and subdivided, so that the calculation efficiency is improved.
In addition, the inverse transform of the wave number domain back to the space domain applies a three-dimensional arbitrary sampling Fourier inverse transform method, and arbitrary interval sampling of wave numbers is also realized. The non-uniform sampling of wave number can enable the region with intense spectral change and strong energy to be encrypted and divided, and sparse sampling is conducted in the region with slow spectral change and weak energy, so that the accuracy and the efficiency of the algorithm are improved, the speed of spectral change can be integrally evaluated from a spectrogram, the region with intense spectral change and strong energy is generally close to the energy peak of the spectrogram, the spectral energy of a gravity field is concentrated in a wavelet number interval, the frequency spectrum is symmetrical about the wave number of 0, and the energy is weakened and the spectral change is slowed down along with the increase of the wave number. Therefore, encryption is generally performed in a small wave number range, and as the wave number increases, sampling can be sparse.
In addition, the gravity potential is solved by performing three-dimensional random sampling Fourier transform on the Poisson equation satisfied by the gravity potential to a spectrum domain and then converting the Fourier transform of random sampling to a solution of a space domain. The method has the advantages that the equation is simple, and the problems of overlarge calculation amount and low efficiency of the traditional Fourier spectrum method for solving the gravity field are greatly improved by means of high efficiency and high precision of any Fourier transform.
In addition to the objects, features and advantages described above, other objects, features and advantages of the present invention are also provided. The present invention will be described in further detail below with reference to the accompanying drawings.
Drawings
The accompanying drawings, which are incorporated in and constitute a part of this application, are included to provide a further understanding of the invention, and are incorporated in and constitute a part of this specification. In the drawings:
FIG. 1 is a schematic view of a calculation region and an abnormal region in a preferred embodiment of the present invention;
FIG. 2 is a schematic view of a model subdivision in a preferred embodiment of the present invention, in which (a) represents a model in which x, y, and z are uniformly sampled; (b) A model representing x, y non-uniform sampling and z uniform sampling; (c) a model in which x, y, z are all non-uniformly sampled;
FIG. 3 is a schematic diagram of the sectional uniform sampling of the wavenumber domain in the preferred embodiment of the present invention;
FIG. 4 is a schematic diagram of a cell and a node within the cell in accordance with a preferred embodiment of the present invention;
FIG. 5 is a comparison graph of a numerical solution and an analytical solution of a sphere in a preferred embodiment of the present invention, wherein (a) represents g x Numerical solution, (b) graph represents g x Analytic solution, (c) graph shows g x Absolute error of numerical solution compared to analytic solution; (d) Denotes g y Numerical solution, (e) diagram representation g y Analytic solution, (f) graph representationg y Absolute error of the numerical solution compared to the analytic solution; (g) Denotes g z Numerical solution, (h) diagram represents g z Analytically solving, (i) graph represents g z Absolute error of numerical solution compared to analytic solution;
fig. 6 is a flow chart of a three-dimensional gravity field numerical simulation method based on 3D AS-FT according to a preferred embodiment of the present invention.
Detailed Description
The embodiments of the invention will be described in detail below with reference to the drawings, but the invention can be implemented in many different ways as defined and covered by the claims.
The first embodiment is as follows:
in the actual geophysical gravity inversion, the actual geological model is reversely deduced from the measured data, and numerical simulation is needed in the inversion process, so that the precision and the efficiency of the selected numerical simulation method are very important. The numerical simulation and inversion are the reverse processes, namely the actual geological anomaly model is known, and field value distribution is obtained. Therefore, the gravity numerical simulation needs to establish a geological model which is consistent with the actual situation, such as the residual density distribution situation of underground three-dimensional in a certain region, and then obtains the gravity field value distribution of any point in the three-dimensional space of the model region according to the established residual density distribution, thereby providing guarantee for the high-precision inversion of gravity and being beneficial to the high-precision imaging of residual density in gravity exploration.
The 3D AS-FT method is applied to gravity field numerical simulation, the space domain and the wave number domain can be sampled at will, the number of points can be encrypted correspondingly at the position where the abnormal body or the spectrum changes violently, and the number of sampling points can be reduced at the position where the abnormal body or the spectrum changes slowly, so that higher precision and efficiency are achieved.
The research scheme of the invention is as follows:
the invention provides a three-dimensional gravity field Fourier spectrum method of 3D AS-FT based on shape function interpolation, which comprises the following steps:
the method comprises the following steps: model building
And completing geological modeling work on the numerical simulation calculation area. The size of the whole calculation area is determined, and then the distribution of the abnormal body is determined, wherein the abnormal body can be any complex condition and any complex shape, and the abnormal body is required to be in the calculation area. A schematic diagram of a simple model is shown in FIG. 1, in which the anomaly is a sphere.
Step two: model subdivision
Modeling in a spatial domain:
after the model is built, the model is subdivided, and the sampling points in the x direction, the y direction and the z direction are Nx, ny direction and Nz direction respectively. One of the advantages of the invention is that the model subdivision is arbitrary in the x, y and z directions, non-uniform subdivision can be adopted in the first area where the abnormal body of the model changes rapidly, encryption is carried out, and sparse sampling is carried out in the second area where the abnormal body changes slowly or does not change. It is also possible that all three directions are uniformly sampled. The sampling three-dimensional schematic diagrams are respectively shown in FIGS. 2a, b and c, wherein a is uniform sampling in three directions; b is x, y direction non-uniform sampling, z direction uniform sampling; c is non-uniform sampling in x, y and z directions.
Determining wave number k according to the space domain subdivision x ，k y ，k z Cut-off frequency (k) of x ，k y ，k z Maximum positive value and minimum negative value) and k x ，k y ，k z The sampling manner of (1).
The cut-off frequency is related to the minimum subdivision interval in the corresponding direction of the space domain, and the minimum subdivision interval in the x direction is set to be Deltax min The minimum split in the y-direction is Δ y min The minimum split interval in the z direction is Δ z min Then the corresponding cutoff frequency is
Sampling within the cut-off frequency ensures that all spectral information is sampled. After the cut-off frequency is determined, the number of samples is determined, assuming k x ，k y ，k z The number of samples of (1) is Nkx, nky, nkz, respectively.
The uniform sampling, i.e. k, can be chosen x ，k y ，k z The arrangement intervals are the same; alternatively, uniform sampling in the log domain may be used.
When sampling in logarithmic interval, setting wave number selection range as [ -k ] max ,k max ]The sampling point number in the wave number domain is 2M +1, and sampling is carried out in the log domain at equal intervals, wherein the sampling interval is
Wherein k is min Is a decimal fraction, generally 10 -6 ～10 -3 。
Wave number is arranged in [ -k [ ] max ,0]Is arranged at
The wave number is arranged in [0, k ] max ]Is provided with
k x ，k y ，k z The logarithmic domain sampling of (2) can be performed by a sum sampling mode, thereby giving the arrangement of the space domain x, y, z and the wave number domain kx, ky, kz.
Wherein the first region satisfies the following formula:
ρ i is the residual density, rho, of the corresponding node of the first region j Correspond to the first regionThe residual density of the jth node around the node, and n is the number of nodes around the first area; omega is weight and the value range is (0, 1).
The second region satisfies the following formula:
ρ i′ the residual density, rho, of the corresponding node of the second region j′ The residual density of the jth node around the corresponding node of the second area is shown, and n is the number of the nodes around the first area; omega is weight and the value range is (0, 1); omega takes a value of 0.2;
and a segmented uniform subdivision mode can be also adopted in the wave number domain. As shown in fig. 3.
Step three: model parameter residual density assignment
Assigning a residual density to the nodes in fig. 2 or fig. 3. The residual density is expressed as rho, and is a scalar quantity with the unit of kg/m 3 。
Step four: by means of a counterweight U g The three-dimensional Poisson equation is subjected to three-dimensional Fourier transform to obtain the wave number domain gravity potential
Spatial domain gravitational potential U g And the residual density p satisfy the poisson equation
Wherein, U g Is a space domain gravitational potential,. Pi.is a circumferential rate, G is a universal gravitational constant, and has a value of 6.672 x 10 -11 m 3 /(kg·s 2 ) Rho is the residual density of the space domain and has the unit of kg/m 3 。
And performing three-dimensional Fourier transform on the formula, wherein the Fourier transform adopts three-dimensional arbitrary Fourier transform based on shape function interpolation, and the principle of the three-dimensional Fourier transform is as follows:
the three-dimensional Fourier transform formula is
In the formula k x 、k y 、k z Representing wave number, F (x, y, z) as a function of the spatial domain, F (k) x ,k y ,k z ) Representing a wavenumber spectrum.
The three-dimensional transformation is completed by three times of one-dimensional Fourier forward transformation, and firstly, the principle of the one-dimensional Fourier forward transformation is introduced.
The one-dimensional fourier forward transform can be represented as:
wherein k is x Representing wave number, F (x) being a function of the spatial domain, F (k) x ) Is a wavenumber spectrum.
The positive transformation integral in the formula is dispersed to obtain
Wherein M represents the number of cells, e j Denotes the jth cell, where i is an imaginary number.
F (x) is interpolated by a quadratic function. When the quadratic interpolation shape function fitting is adopted in the unit, the coordinates of three nodes in any unit are respectively set as x 1 、x 2 、x 3 ，x 2 Is a midpoint, satisfies x 1 +x 3 ＝2x 2 The nodes in the unit are shown in fig. 4:
the value at each node is f (x) 1 )、f(x 2 )、f(x 3 ) F (x) is expressed by a quadratic function
f(x)＝N 1 f(x 1 )+N 2 f(x 2 )+N 3 f(x 3 ) (11)
Wherein,
can be written as
Order to
For the in-cell Fourier transform node coefficients, the formula is abbreviated as
When wave number k x When not 0, substituting W with formula 1 、W 2 、W 3 In, the in-cell Fourier transform node coefficients are obtained
W 1 、W 2 、W 3 The integral kernel functions all comprise
Which is in [ x ] 1 ,x 3 ]Upper unit integral resolution into
Thus, k can be obtained x W when not 0 1 、W 2 、W 3 Semi-analytic solution to
When wave number k x When the average molecular weight is 0, the average molecular weight,
the Fourier transform node coefficient under zero wave number can be obtained by simple integration
And accumulating the analytical expressions of different units to obtain a final one-dimensional Fourier forward transformation result. It is easy to know that when the space domain and the frequency domain are divided invariably, the Fourier transform node coefficient W 1 、W 2 、W 3 And W 1 0 、
The Fourier transform coefficients are all unchanged, the Fourier transform coefficients are calculated and stored in advance, repeated calculation can be reduced, and the algorithm efficiency is improved, which is one of the advantages of the algorithm.
The three-dimensional Fourier transform is to perform one-dimensional Fourier transform on x
Then pair F (k) x Y, z) performs a y-direction one-dimensional fourier transform:
last pair F xy (k x ,k y Z) performing a z-direction one-dimensional fourier transform:
the principle of the three-time one-dimensional Fourier transform is completely the same as the process, and therefore, the description is omitted.
Obtaining a wave number domain through three-dimensional Fourier transform of arbitrary sampling
The following equation is satisfied:
wherein,
is the wavenumber domain residual density distribution.
The wave number domain gravity potential can be obtained by directly solving the formula
Step five: according to the relation between the gravity potential and the gravity field, solving the gravity field of the wave number domain
Space domain gravitational potential U g In relation to the spatial domain gravitational field g
Wherein
i, j, k are unit vectors in x, y, z directions, respectively. Thus obtaining the wave number domain gravity potential
Sum wave number domain gravitational field
In a relationship of
According to the formula, the wave number domain gravity field can be obtained by solving
The unit is milliGal (mGal).
Step six: and obtaining a spatial domain gravity field g through three-dimensional Fourier inverse transformation of arbitrary sampling.
The application of the randomly sampled three-dimensional inverse Fourier transform is also a great innovation of the invention, and the randomly sampled three-dimensional inverse Fourier transform can be ensured when the three-dimensional inverse Fourier transform returns to the space domain during the gravity numerical simulation, so that the precision and the efficiency are improved.
The three-dimensional arbitrary sampling Fourier inverse transformation formula is
In the formula k x 、k y 、k z Representing the wavenumber, F (x, y, z) is a function of the spatial domain, F (k) x ,k y ,k z ) Representing a wavenumber spectrum. The reverse transformation formula has the same form and principle as the forward transformation formula, and is not described again.
The precision and efficiency of the three-dimensional gravity Fourier spectrum method based on the three-dimensional Fourier transform of arbitrary sampling provided by the invention are tested.
The test computer is configured to i7-11800, the main frequency is 2.30GHz, and the memory is 32GB.
Design sphere model, 1000m × 1000m × 1000m, range: x-500-500m, y-500 m and z-0-1000 m. The center of the abnormal sphere model is (0m, 0m and 500m), the sphere radius is 250m, and the residual density of the sphere is set to be 2000kg/m 3 The density of the surrounding rock is 0kg/m 3 . And uniformly dividing the space domain in three directions, and taking 101 nodes. The wave number domain sampling range is based on the sampling theorem, and the maximum wave number is the cut-off frequency
The deltax is a fixed value of 10m, so the maximum wave number is pi/10, a mode of logarithmic domain uniform sampling is adopted, and the minimum logarithmic number is 10 -4 And the number of sampling points in three directions is 101. Calculated, the ground field value g x ,g y ,g z The relative root mean square errors are respectively 0.39%, 0.39% and 0.41% by comparing the numerical solution of (2) with the analytic solution of the sphere, as shown in fig. 5, the absolute errors of the relative root mean square errors and the analytic solution of (5) have more than two orders of magnitude difference, and the general precision requirement is met. Occupies 0.8GB of the memory and takes 0.9s.
In summary, AS shown in fig. 6, the present application provides a three-dimensional gravity field numerical simulation method based on 3D AS-FT, including the following steps: constructing a three-dimensional space numerical model of the target area based on the actual terrain of the target area and the actual distribution condition of the abnormal body; assigning residual density to each node in the three-dimensional space numerical model, and constructing a three-dimensional Poisson equation of the target area based on the three-dimensional space numerical model after assignment; transforming the three-dimensional Poisson equation to a wave number domain by adopting three-dimensional Fourier transform of arbitrary sampling, solving the wave number domain gravity potential of the target area, and solving the wave number domain gravity field of the target area according to the relationship between the wave number domain gravity potential and the wave number domain gravity field; and performing randomly sampled three-dimensional Fourier inverse transformation on the wave number domain gravity field of the target area to obtain a space domain gravity field of the target area. And transforming the three-dimensional gravity differential equation into a spectral domain for direct solution, and then obtaining a solution of a spatial domain through three-dimensional random sampling Fourier inversion. By means of the high precision and high efficiency of the 3D AS-FT, the Fourier spectrum method of the three-dimensional gravity field can achieve high precision, has high calculation efficiency, can sample randomly, and has important significance for accurate imaging of gravity three-dimensional inversion.
The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims (10)
1. A three-dimensional gravity field numerical simulation method based on 3D AS-FT is characterized by comprising the following steps:
constructing a three-dimensional space numerical model of the target area based on the actual terrain of the target area and the actual distribution condition of the abnormal body;
assigning residual density to each node in the three-dimensional space numerical model, and constructing a three-dimensional Poisson equation of the target area based on the assigned three-dimensional space numerical model;
transforming the three-dimensional Poisson equation to a wave number domain by adopting three-dimensional Fourier forward transform of arbitrary sampling, solving a wave number domain gravity potential of the target area, and solving a wave number domain gravity field of the target area according to the relationship between the wave number domain gravity potential and the wave number domain gravity field;
and performing randomly sampled three-dimensional Fourier inverse transformation on the wave number domain gravity field of the target area to obtain a space domain gravity field of the target area.
2. The 3D AS-FT based three-dimensional gravitational field numerical simulation method of claim 1, wherein the arbitrarily sampled three-dimensional Fourier transform is:
carrying out x-direction one-dimensional Fourier forward transform on the three-dimensional Poisson equation;
wherein x, y, z represent three mutually perpendicular directions; k is a radical of formula x Representing the wave number in the x direction, F is a spatial domain function, and F is a wave number spectrum;
to F (k) x Y, z) is subjected to a one-dimensional Fourier transform in the y direction, k y Wavenumber representing y direction:
to F is aligned with xy (k x ,k y Z) one-dimensional in the z-directionFourier transform, k z Wavenumber representing z direction:
3. the 3D AS-FT based three-dimensional gravitational field numerical simulation method of claim 2, wherein said arbitrary sampled three-dimensional Fourier transform is a three-dimensional arbitrary Fourier transform using shape function interpolation.
4. The 3D AS-FT based three-dimensional gravitational field numerical simulation method of claim 3, wherein the one-dimensional Fourier transform is specifically:
let the continuous one-dimensional fourier transform be expressed as:
dispersing continuous one-dimensional Fourier transform integral to obtain:
wherein M represents the number of cells, e j Denotes the jth cell, where i is an imaginary number, k x Represents the x-direction wave number;
interpolating f (x) with a quadratic function:
when the quadratic interpolation shape function fitting is adopted in the unit, the coordinates of three nodes in any unit are respectively set as x 1 、x 2 、x 3 ，x 2 Is a midpoint, satisfies x 1 +x 3 ＝2x 2 The value at each node is f (x) 1 )、f(x 2 )、f(x 3 ) And f (x) is expressed by a quadratic function:
f(x)＝N 1 f(x 1 )+N 2 f(x 2 )+N 3 f(x 3 )
wherein,
the above formula is written as:
order to
For the in-cell fourier transform node coefficients, the above equation is abbreviated as:
when wave number k x When not 0, substituting the above formula into W 1 、W 2 、W 3 Obtaining an in-cell Fourier transform node coefficient:
W 1 、W 2 、W 3 the integral kernel functions all comprise
Which is in [ x ] 1 ,x 3 ]The upper unit integral is resolved into:
thus, k can be obtained x W when not 0 1 、W 2 、W 3 The semi-analytic solution is:
when wave number k x When the average molecular weight is 0, the average molecular weight,
the simple integration is carried out to obtain the Fourier transform node coefficient under the zero wave number as follows:
and accumulating the analytical expressions of different units to obtain a final one-dimensional Fourier forward transformation result.
5. The 3D AS-FT-based three-dimensional gravity field numerical simulation method of claim 3, wherein a wave number domain gravity potential is obtained through three-dimensional Fourier forward transform of arbitrary sampling
Satisfies the equation:
wherein,
solving to obtain the number domain gravitational potential for the wave number domain residual density distribution according to the formula
6. The 3D AS-FT based three-dimensional gravity field numerical simulation method of claim 4, characterized in that wave number domain gravity potential
Sum wave number domain gravitational field
The relationship of (c) is:
according to the formula, the wave number domain gravity field is obtained through solving
In milli-gal.
7. The 3D AS-FT based three-dimensional gravity field numerical simulation method of claim 6, wherein a three-dimensional arbitrary sampling inverse Fourier transform formula is AS follows:
in the formula k x 、k y 、k z Representing wave number, F (x, y, z) as a function of the spatial domain, F (k) x ,k y ,k z ) Representing a wavenumber spectrum.
8. The 3D AS-FT based three-dimensional gravitational field numerical simulation method according to any one of claims 1-7, wherein the three-dimensional Poisson equation is:
wherein, U g Is space domain gravitational potential, pi is circumferential rate, G is universal gravitation constant, and rho is space domain residual density with unit of kg/m 3 。
9. The 3D AS-FT based three-dimensional gravity field numerical simulation method according to any one of claims 1-7, characterized in that the spatial domain of the target model is subdivided in any one of the following ways to construct a three-dimensional spatial numerical model of the target region,
the first method is as follows: non-uniform subdivision is carried out on a preset first area, and encryption is carried out; wherein the first region satisfies the following formula:
ρ i residual density, rho, of nodes corresponding to the first region j The residual density of the jth node around the corresponding node of the first area is obtained, and n is the number of the nodes around the first area; omega is weight and the value range is (0, 1);
the second method comprises the following steps: sparse sampling is carried out on a preset second region, wherein the second region meets the following formula:
ρ i′ the residual density, rho, of the corresponding node of the second region j′ The residual density of the jth node around the corresponding node of the second area is shown, and n is the number of the nodes around the first area; omega is weight and the value range is (0, 1);
the third method comprises the following steps: the spatial domain samples uniformly in all three directions.
10. A computer system comprising a memory, a processor and a computer program stored on the memory and executable on the processor, wherein the processor implements the steps of the method of any one of claims 1 to 9 when executing the computer program.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211076505.7A CN115659579B (en) | 2022-09-05 | 2022-09-05 | Three-dimensional gravity field numerical simulation method and system based on 3D AS-FT |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211076505.7A CN115659579B (en) | 2022-09-05 | 2022-09-05 | Three-dimensional gravity field numerical simulation method and system based on 3D AS-FT |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115659579A true CN115659579A (en) | 2023-01-31 |
| CN115659579B CN115659579B (en) | 2023-07-04 |
Family
ID=84983377
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211076505.7A Active CN115659579B (en) | 2022-09-05 | 2022-09-05 | Three-dimensional gravity field numerical simulation method and system based on 3D AS-FT |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115659579B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116187107A (en) * | 2023-04-27 | 2023-05-30 | 中南大学 | Three-dimensional ground temperature field dynamic numerical simulation method, equipment and medium |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130166212A1 (en) * | 2011-12-21 | 2013-06-27 | Technoimaging, Llc | Method of terrain correction for potential field geophysical survey data |
| CN108984939A (en) * | 2018-07-30 | 2018-12-11 | 中南大学 | Three-dimensional Gravity field of force forward modeling method based on 3D Gauss-FFT |
| CN109100804A (en) * | 2018-06-11 | 2018-12-28 | 中国石油天然气集团有限公司 | Improve data re-establishing method, the apparatus and system of seismic data spatial sampling attribute |
| CN111505596A (en) * | 2020-04-16 | 2020-08-07 | 北京理工大学重庆创新中心 | Three-dimensional wind field inversion method based on non-uniform sampling correction VAD technology |
| CN111967169A (en) * | 2020-10-21 | 2020-11-20 | 中南大学 | Two-degree body weight abnormal product decomposition numerical simulation method and device |
| CN112287534A (en) * | 2020-10-21 | 2021-01-29 | 中南大学 | NUFFT-based two-dimensional magnetic anomaly fast forward modeling method and device |
| CN112800657A (en) * | 2021-04-15 | 2021-05-14 | 中南大学 | Gravity field numerical simulation method and device based on complex terrain and computer equipment |
| CN113626757A (en) * | 2021-08-12 | 2021-11-09 | 南京工程学院 | Two-dimensional discrete data smoothing method based on fast Fourier transform |
| CN113779498A (en) * | 2021-08-03 | 2021-12-10 | 中国电子产品可靠性与环境试验研究所（（工业和信息化部电子第五研究所）（中国赛宝实验室）） | Discrete Fourier matrix reconstruction method, device, equipment and storage medium |
-
2022
- 2022-09-05 CN CN202211076505.7A patent/CN115659579B/en active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130166212A1 (en) * | 2011-12-21 | 2013-06-27 | Technoimaging, Llc | Method of terrain correction for potential field geophysical survey data |
| CN109100804A (en) * | 2018-06-11 | 2018-12-28 | 中国石油天然气集团有限公司 | Improve data re-establishing method, the apparatus and system of seismic data spatial sampling attribute |
| CN108984939A (en) * | 2018-07-30 | 2018-12-11 | 中南大学 | Three-dimensional Gravity field of force forward modeling method based on 3D Gauss-FFT |
| CN111505596A (en) * | 2020-04-16 | 2020-08-07 | 北京理工大学重庆创新中心 | Three-dimensional wind field inversion method based on non-uniform sampling correction VAD technology |
| CN111967169A (en) * | 2020-10-21 | 2020-11-20 | 中南大学 | Two-degree body weight abnormal product decomposition numerical simulation method and device |
| CN112287534A (en) * | 2020-10-21 | 2021-01-29 | 中南大学 | NUFFT-based two-dimensional magnetic anomaly fast forward modeling method and device |
| CN112800657A (en) * | 2021-04-15 | 2021-05-14 | 中南大学 | Gravity field numerical simulation method and device based on complex terrain and computer equipment |
| CN113779498A (en) * | 2021-08-03 | 2021-12-10 | 中国电子产品可靠性与环境试验研究所（（工业和信息化部电子第五研究所）（中国赛宝实验室）） | Discrete Fourier matrix reconstruction method, device, equipment and storage medium |
| CN113626757A (en) * | 2021-08-12 | 2021-11-09 | 南京工程学院 | Two-dimensional discrete data smoothing method based on fast Fourier transform |
Non-Patent Citations (5)
| Title |
|---|
| ALDEN S.JURLING ET AL: "Techniques for arbitrary sampling in two-dimensional Fourier transforms", JOURNAL OF THE OPTICAL SOCIETY OF AMERICA A, vol. 35, no. 11, pages 1784 - 1796 * |
| 周印明;戴世坤;李昆;何展翔;胡晓颖;王金海;: "基于样条插值的FFT及其在重磁场正演中的应用", 石油地球物理勘探, vol. 55, no. 4, pages 915 - 922 * |
| 周印明;戴世坤;李昆;凌嘉宣;胡晓颖;熊彬;: "复杂形体重磁位场三维高效高精度数值模拟", 石油地球物理勘探, vol. 55, no. 05, pages 1149 - 1159 * |
| 戴世坤;王旭龙;赵东东;刘志伟;张钱江;孙金飞;: "空间―波数混合域二度体重力异常正演", 石油地球物理勘探, no. 06, pages 18 - 19 * |
| 戴世坤;陈轻蕊;李昆;凌嘉宣;: "重力异常场空间-波数混合域三维数值模拟", 地球物理学报, no. 05, pages 2107 - 2119 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116187107A (en) * | 2023-04-27 | 2023-05-30 | 中南大学 | Three-dimensional ground temperature field dynamic numerical simulation method, equipment and medium |
| CN116187107B (en) * | 2023-04-27 | 2023-08-11 | 中南大学 | Three-dimensional ground temperature field dynamic numerical simulation method, equipment and medium |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115659579B (en) | 2023-07-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109856689B (en) | Noise suppression processing method and system for superconducting aeromagnetic gradient tensor data | |
| CN106646645B (en) | A kind of gravity forward modeling accelerated method | |
| CN112287534A (en) | NUFFT-based two-dimensional magnetic anomaly fast forward modeling method and device | |
| CN115292973B (en) | Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system | |
| CN110852025B (en) | Three-dimensional electromagnetic slow diffusion numerical simulation method based on super-convergence interpolation approximation | |
| Mirzaei et al. | Direct meshless local Petrov–Galerkin method for elastodynamic analysis | |
| CN115659579A (en) | Three-dimensional gravity field numerical simulation method and system based on 3D AS-FT | |
| Wu | Comparison of 3-D Fourier forward algorithms for gravity modelling of prismatic bodies with polynomial density distribution | |
| CN113608262B (en) | Seismic data processing method and device for calculating rotation component by using translation component | |
| Roux et al. | Analysis of numerically induced oscillations in two-dimensional finite-element shallow-water models part II: Free planetary waves | |
| Qiu et al. | Gravity field of a tesseroid by variable-order Gauss–Legendre quadrature | |
| Chan et al. | Volume statistics as a probe of large-scale structure | |
| CN113792445B (en) | Three-dimensional magnetotelluric numerical simulation method based on integral equation method | |
| Peter et al. | Surface wave tomography: global membrane waves and adjoint methods | |
| Galewski | Spectrum-based modal parameters identification with Particle Swarm Optimization | |
| Wang et al. | Fast Nonlinear generalized inversion of gravity data with application to the three-dimensional crustal density structure of Sichuan Basin, Southwest China | |
| CN106353801B (en) | The three-dimensional domain Laplace ACOUSTIC WAVE EQUATION method for numerical simulation and device | |
| Loeffler et al. | Comparison between the formulation of the boundary element method that uses fundamental solution dependent of frequency and the direct radial basis boundary element formulation for solution of Helmholtz problems | |
| CN116244877B (en) | Three-dimensional magnetic field numerical simulation method and system based on 3D Fourier transform | |
| Geng et al. | Reconstruction and analysis of a non-stationary sound field in a uniformly flow medium | |
| CN114325870A (en) | Method and system for calculating potential field gradient tensor based on cubic spline function | |
| Helland et al. | Bispectra and energy transfer in grid-generated turbulence | |
| Wu et al. | Experimental study of the mapping relationship based near-field acoustic holography with spherical fundamental solutions | |
| Guo et al. | Microwave inversion for sparse data using descent learning technique | |
| Baxansky et al. | Calculating geometric properties of three-dimensional objects from the spherical harmonic representation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |