CN115659579B - 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
- CN115659579B CN115659579B CN202211076505.7A CN202211076505A CN115659579B CN 115659579 B CN115659579 B CN 115659579B CN 202211076505 A CN202211076505 A CN 202211076505A CN 115659579 B CN115659579 B CN 115659579B
- Authority
- CN
- China
- Prior art keywords
- dimensional
- domain
- wave number
- gravity field
- gravity
- 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.)
- Active
Links
- 230000005484 gravity Effects 0.000 title claims abstract description 81
- 238000000034 method Methods 0.000 title claims abstract description 54
- 238000004088 simulation Methods 0.000 title claims abstract description 36
- 238000005070 sampling Methods 0.000 claims abstract description 58
- 230000009466 transformation Effects 0.000 claims abstract description 33
- 238000001228 spectrum Methods 0.000 claims description 22
- 230000006870 function Effects 0.000 claims description 20
- 230000002159 abnormal effect Effects 0.000 claims description 13
- 238000012887 quadratic function Methods 0.000 claims description 6
- 238000004458 analytical method Methods 0.000 claims description 5
- 238000004590 computer program Methods 0.000 claims description 4
- 230000010354 integration Effects 0.000 claims description 4
- 230000001131 transforming effect Effects 0.000 claims description 4
- 230000014509 gene expression Effects 0.000 claims description 3
- 238000012876 topography Methods 0.000 claims description 3
- 230000003595 spectral effect Effects 0.000 abstract description 13
- 238000004422 calculation algorithm Methods 0.000 abstract description 5
- 238000004364 calculation method Methods 0.000 description 13
- 230000008859 change Effects 0.000 description 10
- 238000010586 diagram Methods 0.000 description 4
- 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
- 210000000746 body region Anatomy 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000004611 spectroscopical analysis Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000011426 transformation method 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 system based on 3D random sampling Fourier transform (Fourier transform of arbitrary sampling, abbreviated AS AS-FT), which are characterized in that a Poisson equation satisfied by gravity is Fourier transformed through the three-dimensional random sampling Fourier transform, the gravity is directly obtained in a spectral domain, the solution is simple, the relation between the gravity of the spectral domain and the gravity field is used for obtaining the gravity field of the spectral domain, and the gravity field of the spatial domain can be solved by performing three-dimensional Fourier inverse transformation on the gravity field. The algorithm has good parallelism and occupies less memory. Therefore, the spectral domain gravity numerical simulation which can be arbitrarily sampled in the space domain and the wave number domain is realized, and the method has good guiding significance for the condition of non-uniform sampling. The algorithm has higher precision and efficiency by means of Fourier transformation which can be sampled arbitrarily.
Description
Technical Field
The invention relates to the field of gravity field numerical simulation, in particular to a 3D AS-FT-based three-dimensional gravity field numerical simulation method and system.
Background
The potential field high-precision numerical simulation has important significance for geophysical interpretation and inversion, and the spectral method is a method commonly used in gravitational field numerical simulation, and is characterized in that a solution is approximately expanded into a finite series expansion of a smooth function (generally an orthogonal polynomial), so that an approximate spectral expansion of the solution is obtained, and then an equation set of expansion coefficients is solved according to the expansion and an original equation. Spectral methods are essentially a generalization of standard separation variant techniques. Typically Chebyshev polynomials and Legendre polynomials are taken as basis functions for the approximate expansion. For equations meeting periodic boundary conditions, the Fourier series and the harmonic series are simpler. The spectroscopic methods are classified as Fourier methods, chebyshev or Legendre methods. The former is applicable to periodic problems, and the latter two are applicable to non-periodic problems. The basis of the methods is to establish a space basis function, and the precision of the method directly depends on the number of the expansion terms, so that a plurality of terms are needed to achieve better precision, thereby causing large calculation amount and low efficiency.
The patent adopts a Fourier spectrum method, 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 such as aliasing, edge effects, added periodicity, and truncation effects when using a Discrete Fourier Transform (DFT) instead of a continuous fourier transform. If standard FFT is adopted, edge expansion processing is needed to weaken the influence of boundary effect, but the cost is increased in calculation scale; if Gauss-FFT is used for numerical simulation, the efficiency is greatly reduced. And whether standard FFT or Gauss-FFT is used, the sampling can only be uniformly performed, and the waste of calculation capacity is further caused.
Methods for applying fourier transform to gravitational field numerical simulation in the spectral domain have been used in related applications. The literature (Three-dimensional numerical modeling ofgravity and magnetic anomaly in a mixed space-wavenumber domain, geophysics,2018,Shikun Dai,Dongdong Zhao,etal) carries out numerical simulation on gravitational fields by respectively applying a standard FFT and Gauss-FFT method from a gravity differential equation, and the literature (g-force abnormal field space-wave number mixed domain Three-dimensional numerical simulation, geophysics theory report, 2020, dai Shikun, chen Qingrui and the like) carries out numerical simulation on gravitational fields by respectively applying a standard FFT and Gauss-FFT method from an integral equation; however, it is pointed out that the standard FFT must spread the edges to achieve the same precision as the Gauss-FFT, whereas Gauss-FFT, although having higher precision, has smaller boundary effects, is relatively inefficient, and both standard FFT and Gauss-FFT methods can only sample unevenly, which is not flexible enough. The literature (High-accuracy 3D Fourier forwardmodeling of gravity field bFIed on the Gauss-FFT technique, journal of Applied Geophysics,2018,Guangdong Zhao,Bo Chen,etal) adopts 3D Gauss-FFT to perform numerical simulation, and the method realizes High-precision numerical simulation of a gravity three-dimensional spectral domain, but Gauss-FFT occupies higher memory and can only perform non-uniform sampling.
Disclosure of Invention
The invention provides a 3D AS-FT-based three-dimensional gravity field numerical simulation method and a system, which are used for solving the technical problem that the efficiency and the accuracy of the existing gravity field numerical simulation method are incompatible.
In order to solve the technical problems, the technical scheme provided by the invention is as follows:
A3D AS-FT-based three-dimensional gravity field numerical simulation method 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 densities 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 transformation of arbitrary sampling, solving a wave number domain gravity position of a target area, and solving a wave number domain gravity field of the target area according to the relation between the wave number domain gravity position and the wave number domain gravity field;
and performing three-dimensional Fourier inverse transformation on the arbitrary sampling wave number domain gravitational field of the target area to obtain the spatial domain gravitational field of the target area.
Preferably, the three-dimensional fourier transform of the arbitrary samples 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 (k) x The wave number in the x direction is represented, F is a spatial domain function, and F is a wave number spectrum;
for F (k) x Y, z) performs a one-dimensional fourier positive transform in the y-direction:
for F xy (k x ,k y Z) performing a one-dimensional fourier positive transform in the z direction:
preferably, the three-dimensional fourier positive transformation of the arbitrary sample is a three-dimensional arbitrary fourier positive transformation using shape function-based interpolation.
Preferably, the one-dimensional fourier positive transform is specifically:
let the continuous one-dimensional fourier positive transforms be expressed as:
discrete one-dimensional fourier positive transform integration is obtained:
wherein M represents the number of units, e j Represents the j-th cell, where i is an imaginary number;
interpolation of f (x) with a quadratic function:
when the quadratic interpolation shape function fitting is adopted in the units, three node coordinates in any unit are respectively set as x 1 、x 2 、x 3 ，x 2 Is the midpoint, satisfy x 1 +x 3 ＝2x 2 Each node has a value f (x 1 )、f(x 2 )、f(x 3 ) F (x) is expressed by using a quadratic function:
f(x)＝N 1 f(x 1 )+N 2 f(x 2 )+N 3 f(x 3 )
wherein,,
the above can be written as:
order the
For the intra-cell fourier transform node coefficients, the above abbreviation is:
when the wave number is k x When the value is not 0, substituting the above formula into W 1 、W 2 、W 3 In (2), the intra-cell fourier transform node coefficients can be obtained:
W 1 、W 2 、W 3 the integral kernel functions all include
It is at [ x ] 1 ,x 3 ]The upper unit integral analysis solution is:
thus, k can be obtained x W when it is not 0 1 、W 2 、W 3 The semi-analytical solution is:
when the wave number is k x When the number of the organic light emitting diode is 0,
and (3) carrying out simple integration to obtain the Fourier transform node coefficient under zero wave number as follows:
and accumulating the analysis expressions of different units to obtain a final one-dimensional Fourier positive transformation result.
Preferably, the wave number domain gravity position is obtained by three-dimensional Fourier positive transformation of arbitrary sampling
The equation is satisfied:
wherein,,
for the wavenumber domain residual density distributionAccording to the above, the gravity position of the number domain can be obtained by solving
Preferably, the wave number domain gravity position
And wave number domain gravity field->
The relation of (2) is:
according to the above, the wave number domain gravity field can be obtained by solving
In milligamma.
Preferably, the fourier transform formula of the three-dimensional arbitrary sample is:
k in x 、k y 、k z Representing wavenumber, F (x, y, z) is a spatial domain function, F (k) x ,k y ,k z ) Representing the wavenumber spectrum.
Preferably, the three-dimensional poisson equation is:
wherein U is g Is space domain gravitational potential, pi is circumference rate, G is gravitational constant, ρ is space domain residual density, and unit is kg/m 3 。
Preferably, the spatial domain of the target model is split in any one of the following ways to construct a three-dimensional spatial numerical model of the target region,
mode one: non-uniform subdivision is adopted for a preset first area, and encryption is carried out; wherein the first region satisfies the following formula:
ρ i for the residual density ρ 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 given, and n is the number of the nodes around the first area; omega is weight, and the value range is (0, 1).
Mode two: sparse sampling is carried out on a preset second area, wherein the second area meets the following formula:
ρ i′ for the remaining density of the corresponding node of the second region, ρ j′ The remaining density of the jth node around the corresponding node of the second area is given, and n is the number of the nodes around the first area; omega is weight, and the value range is (0, 1);
mode three: 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 described above when the computer program is executed.
The invention has the following beneficial effects:
the invention applies a three-dimensional arbitrary sampling Fourier transform method (3 DAS-FT), transforms a three-dimensional gravity differential equation into a spectral domain to directly solve, and obtains a solution of a spatial domain through three-dimensional arbitrary 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 higher precision, has higher calculation efficiency, can be randomly sampled, and has important significance for accurate imaging of gravity three-dimensional inversion.
In addition, the spatial domain is transformed into the wave number domain by using three-dimensional arbitrary sampling Fourier forward transformation, so that the spatial domain x, y and z directions can be split unevenly. The abnormal body can be better fitted, and the calculation accuracy is improved. Sparse sampling is performed in the area where encryption subdivision is not needed, so that the calculation efficiency is improved.
In addition, the wave number domain inverse transformation back to the space domain in the invention applies a three-dimensional arbitrary sampling Fourier inverse transformation method, and the arbitrary interval sampling of wave numbers is realized. The non-uniform sampling of the wave numbers can enable the region with strong energy and strong spectrum change to be encrypted and split, and the region with weak energy and weak spectrum change is sparsely sampled, so that the accuracy and efficiency of an algorithm are improved, the speed of spectrum change can be integrally estimated from a spectrogram, the vicinity of an energy peak of the spectrogram is generally a region with strong spectrum change and strong energy, the spectrum energy of a gravitational field is concentrated in a small wave number region, the spectrum is symmetrical about the wave number 0, and the energy weakens and slows down along with the increase of the wave number. Therefore, encryption is generally performed between the wavenumbers of the cells, and as the wavenumber increases, sampling may be sparse.
In addition, the gravity bit solving method is to perform three-dimensional arbitrary sampling Fourier transform on a Poisson equation satisfied by the gravity bit to a spectral domain, and then to change the three-dimensional arbitrary sampling Fourier transform to a solution of a spatial domain. The equation is simple, and the problems of overlarge calculation amount and lower efficiency of solving the gravitational field by the traditional Fourier spectrum method 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, the present invention has other objects, features and advantages. The invention will be described in further detail with reference to the accompanying drawings.
Drawings
The accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description serve to explain the invention. In the drawings:
FIG. 1 is a schematic view of a calculation region and an abnormal body region in a preferred embodiment of the present invention;
FIG. 2 is a schematic diagram of a model subdivision in a preferred embodiment of the present invention, wherein (a) represents a model in which x, y, z are uniformly sampled; (b) Representing a model of non-uniform sampling of x and y and uniform sampling of z; (c) a model in which x, y, z are non-uniformly sampled;
FIG. 3 is a schematic view of the wave number domain segment uniform sampling in the preferred embodiment of the present invention;
FIG. 4 is a schematic diagram of a unit and nodes within the unit in a preferred embodiment of the invention;
FIG. 5 is a graph showing a numerical solution and an analytical solution of the sphere in the preferred embodiment of the present invention, wherein (a) represents g x Numerical solution, (b) graph represents g x Resolving solutions, (c) graph represents g x Absolute error of a numerical solution compared to an analytical solution; (d) G represents g y Numerical solution, (e) graph represents g y Resolving solutions, the (f) graph represents g y Absolute error of a numerical solution compared to an analytical solution; (g) G represents g z Numerical solution, (h) graph represents g z Resolving solutions, (i) graph represents g z Absolute error of a numerical solution compared to an analytical solution;
FIG. 6 is a flow chart of a 3D AS-FT based three-dimensional gravity field numerical simulation method in a preferred embodiment of the present invention.
Detailed Description
Embodiments of the invention are described in detail below with reference to the attached drawings, but the invention can be implemented in a number of different ways, which are defined and covered by the claims.
Embodiment one:
in the actual geophysical gravity inversion, the actual geological model is reversely deduced from measured data, numerical simulation is needed in the inversion process, and therefore the accuracy and the efficiency of the selected numerical simulation method are important. Numerical simulation and inversion are the opposite processes, i.e., the actual geologic anomaly model is known, resulting in a field value distribution. Therefore, the gravity numerical simulation needs to build a geological model conforming to the actual situation, such as the situation of the three-dimensional residual density distribution of the underground of a certain area, and then obtains the gravity field value distribution of any point in the three-dimensional space of the model area according to the built residual density distribution, thereby providing guarantee for high-precision inversion of gravity and being beneficial to high-precision imaging of the residual density in gravity exploration.
According to the invention, the 3D AS-FT method is applied to gravity field numerical simulation, the spatial domain and the wave number domain can be sampled at will, the number of encryption points can be correspondingly increased at the position of the abnormal body or the position of the spectrum with intense change, and the number of sampling points can be reduced at the position of the abnormal body or the position of the spectrum with slow change, so that higher precision and efficiency are achieved.
The research scheme of the invention is as follows:
the three-dimensional gravity field Fourier spectrum method of the 3D AS-FT based on shape function interpolation comprises the following steps:
step one: model building
And (5) performing geological modeling on the calculated area of the numerical simulation. The size of the whole calculation region is firstly determined, and then the distribution of abnormal bodies, which can be any complex condition and any complex shape, are determined, wherein the abnormal bodies are required to be in the calculation region. 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 split, and the sampling points in the x, y and z directions are Nx, ny and Nz respectively. One of the advantages of the invention is that the model subdivision is arbitrary in the x, y and z directions, the non-uniform subdivision can be adopted for encryption in a first area with a faster abnormal body change of the model, and sparse sampling can be carried out in a second area with a slower abnormal body change or without change. It is also possible that all three directions are uniformly sampled. The three-dimensional schematic diagrams of the sampling are respectively shown in fig. 2a, b and c, wherein a is uniformly sampled in three directions; b is non-uniform sampling in the x-y direction and uniform sampling in the z-direction; c is the non-uniform sampling in the x, y and z directions.
Determining wavenumber k from spatial domain subdivision x ，k y ，k z The cut-off frequency (k) x ，k y ，k z Maximum positive and minimum negative of (a) and k x ，k y ，k z Is a sampling mode of the system.
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 as Deltax min The minimum subdivision interval in the y direction is delta y min The z-direction minimum subdivision interval is deltaz min The corresponding cut-off frequency is
Sampling within the cut-off frequency can ensure that all spectral information is sampled. After determining the cut-off frequency, the number of samples is determined again, assuming k x ，k y ，k z The number of samples is Nkx, nky, nkz, respectively.
Can select uniform sampling, namely k x ，k y ，k z The arrangement intervals are the same; alternatively, the sampling may be uniform in the logarithmic domain.
When sampling logarithmic interval, the wave number is set to be in the range of [ -k max ,k max ]The number of the wave number domain sampling points is 2M+1, the sampling is carried out on the digital domain sampling at equal intervals, and the sampling interval is
Wherein k is min Is a decimal, typically 10 -6 ～10 -3 。
The wave number is arranged at [ -k max ,0]Upper is
Wavenumbers are arranged in [0, k ] max ]Upper is
k x ，k y ，k z Is sum of logarithmic domain samplesThe sampling pattern thus gives an arrangement of the spatial domain x, y, z and the wavenumber domain kx, ky, kz.
Wherein the first region satisfies the following formula:
ρ i for the residual density ρ 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 given, and n is the number of the nodes around the first area; omega is weight, and the value range is (0, 1).
The second region satisfies the following formula:
ρ i′ for the remaining density of the corresponding node of the second region, ρ j′ The remaining density of the jth node around the corresponding node of the second area is given, 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;
the method can also sample the segmentation uniform subdivision mode in the wave number domain. As shown in fig. 3.
Step three: model parameter residual density assignment
The remaining density assignment is made to the nodes in fig. 2 or fig. 3. The residual density is expressed as ρ, a scalar in kg/m 3 。
Step four: by aligning gravity position U g Performing three-dimensional Fourier transform on the satisfied three-dimensional poisson equation to obtain a wave number domain gravity position
Space domain gravity position U g And the residual density ρ satisfies poisson's equation
Wherein U is g Is space domain gravitational potential, pi is circumference ratio, G is gravitational constant, its value is 6.672×10 -11 m 3 /(kg·s 2 ) ρ is the residual density of the spatial domain in kg/m 3 。
And performing three-dimensional Fourier forward transformation, wherein the three-dimensional Fourier forward transformation adopts three-dimensional arbitrary Fourier forward transformation based on shape function interpolation, and the three-dimensional Fourier forward transformation principle is as follows:
the three-dimensional Fourier positive transformation formula is
K in x 、k y 、k z Representing wavenumber, F (x, y, z) is a spatial domain function, F (k) x ,k y ,k z ) Representing the wavenumber spectrum.
The three-dimensional transformation is completed through three one-dimensional Fourier positive transformation, and the principle of the one-dimensional Fourier positive transformation is described first.
The one-dimensional fourier positive transforms can be expressed as:
wherein k is x Representing wavenumber, F (x) is a spatial domain function, F (k) x ) Is wave number spectrum.
Discrete-obtaining the positive transform integral in
Wherein M represents the number of units, e j Represents the j-th 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 units, three node coordinates in any unit are respectively set as x 1 、x 2 、x 3 ，x 2 Is the midpoint, satisfy x 1 +x 3 ＝2x 2 The intra-cell nodes are as shown in fig. 4:
each node has a value f (x 1 )、f(x 2 )、f(x 3 ) F (x) is obtained by using quadratic function to represent
f(x)＝N 1 f(x 1 )+N 2 f(x 2 )+N 3 f(x 3 ) (11)
Wherein,,
can be written as
Order the
Is the intra-unit Fourier transform node coefficient, abbreviated as
When the wave number is k x When the value is not 0, substituting W 1 、W 2 、W 3 In the method, the intra-cell Fourier transform node coefficients can be obtained
W 1 、W 2 、W 3 The integral kernel functions all include
It is at [ x ] 1 ,x 3 ]Upper unit integral resolution
Thus, k can be obtained x W when it is not 0 1 、W 2 、W 3 Semi-analytical solution into
When the wave number is k x When the number of the organic light emitting diode is 0,
simple integration is carried out to obtain the Fourier transform node coefficient under zero wave number as
And accumulating the analysis expressions of different units to obtain a final one-dimensional Fourier positive transformation result. It is easy to know that when the space domain and frequency domain subdivision is unchanged, the Fourier transform node coefficient W 1 、W 2 、W 3 And W is 1 0 、
The Fourier transform coefficients are calculated and stored in advance, so that repeated calculation can be reduced, the algorithm efficiency is improved, and the method is one of the advantages of the algorithm.
The three-dimensional Fourier transform is that the one-dimensional Fourier transform is completed on x
After that for F (k) x Y, z) performs a one-dimensional fourier transform in the y-direction:
finally to F xy (k x ,k y Z) performing a z-direction one-dimensional fourier transform:
the principle of the three-dimensional fourier transform is exactly the same as that of the process, and thus will not be described in detail.
Obtaining wave number domain through three-dimensional Fourier transform of arbitrary sampling
The following equation is satisfied:
wherein,,
the density profile remains for the wavenumber domain.
Directly solving to obtain the wave number domain gravity position
Step five: solving the wave number domain gravity field according to the relation between the gravity position and the gravity field
Space domain gravity position U g The relation with the gravitational field g of the spatial domain is +.>
Wherein->
i, j, k are unit vectors in the x, y, z directions, respectively. Thus the wave number domain gravity bit is available>
And wave number domain gravity field->
The relation of (2) is that
According to the method, the wave number domain gravity field can be obtained by solving
In milligamma (mGal).
Step six: and obtaining the space domain gravity field g through the three-dimensional Fourier inverse transformation of arbitrary sampling.
The application of arbitrary sampling three-dimensional Fourier inverse transformation is also a great innovation of the invention, and the arbitrary sampling can be ensured when the gravity numerical simulation of the invention is inversely transformed back to the space domain, thereby improving the precision and the efficiency.
The three-dimensional arbitrary sampling Fourier inverse transformation formula is
K in x 、k y 、k z Representing wavenumber, F (x, y, z) is a spatial domain function, F (k) x ,k y ,k z ) Representing the wavenumber spectrum. The form of the reverse transformation formula is identical to that of the forward transformation formula, the principle is also identical, and the description is omitted.
The accuracy 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 checked.
The test computer is configured as i7-11800, the main frequency is 2.30GHz, and the memory is 32GB.
Designing a sphere model, 1000m×1000m, and ranging: x direction-500 m, y direction-500 m, and z direction 0-1000 m. The abnormal sphere model center was located at (0 m,500 m), the sphere radius was 250m, and the sphere residual density was 2000kg/m 3 The density of surrounding rock is 0kg/m 3 . The three directions of the space domain are evenly split, and 101 nodes are taken. The sampling range of the wave number domain is based on the sampling theorem, and the maximum wave number is the cut-off frequency
Δx is a constant value of 10m, and thus the maximum wave number isPi/10, adopting a mode of uniformly sampling logarithmic domain, wherein the logarithmic minimum is 10 -4 And the number of sampling points in three directions is 101. Calculated ground field value g x ,g y ,g z Compared with the analytic solution of the sphere, the relative root mean square error is 0.39%, 0.39% and 0.41%, respectively, and as shown in fig. 5, the absolute error is more than two orders of magnitude different from the analytic solution, and the general precision requirement is met. Occupies 0.8GB of memory and takes 0.9s.
In summary, AS shown in fig. 6, the present application provides a 3D AS-FT-based three-dimensional gravity field numerical simulation method, which includes 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 densities 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 transformation of arbitrary sampling, solving a wave number domain gravity position of a target area, and solving a wave number domain gravity field of the target area according to the relation between the wave number domain gravity position and the wave number domain gravity field; and performing three-dimensional Fourier inverse transformation on the arbitrary sampling wave number domain gravitational field of the target area to obtain the spatial domain gravitational field of the target area. And transforming the three-dimensional gravity differential equation into a spectral domain to directly solve, and obtaining a solution of a spatial domain through three-dimensional arbitrary sampling Fourier inverse transformation. 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 higher precision, has higher calculation efficiency, can be randomly sampled, and has important significance for accurate imaging of gravity three-dimensional inversion.
The above description is only of the preferred embodiments of the present invention and is not intended to limit the present invention, but various modifications and variations can be made to the present invention by those skilled in the art. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims (10)
1. The 3D AS-FT-based three-dimensional gravity field numerical simulation method is characterized by comprising the following steps of:
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 densities 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 transformation of arbitrary sampling, solving a wave number domain gravity position of a target area, and solving a wave number domain gravity field of the target area according to the relation between the wave number domain gravity position and the wave number domain gravity field;
and performing three-dimensional Fourier inverse transformation on the arbitrary sampling wave number domain gravitational field of the target area to obtain the spatial domain gravitational field of the target area.
2. The 3D AS-FT based three-dimensional gravity field numerical simulation method according to claim 1, wherein the arbitrarily sampled three-dimensional fourier 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 (k) x The wave number in the x direction is represented, F is a spatial domain function, and F is a wave number spectrum;
for F (k) x Y, z) performs a one-dimensional fourier positive transform in the y-direction, k y Wavenumber representing y-direction:
for F xy (k x ,k y Z) performing a one-dimensional Fourier forward transform in the z direction, k z Wavenumber representing z-direction:
3. the 3D AS-FT based three-dimensional gravity field numerical simulation method according to claim 2, wherein the arbitrary sampled three-dimensional fourier transform is a three-dimensional arbitrary fourier transform using shape function interpolation.
4. A 3D AS-FT based three-dimensional gravity field numerical simulation method according to claim 3, characterized in that the one-dimensional fourier positive transformation is specifically:
let the continuous one-dimensional fourier positive transforms be expressed as:
discrete the continuous one-dimensional fourier transform integral to obtain:
wherein M represents the number of units, e j Represents the jth cell, where i is an imaginary number, k x Representing the wave number in the x-direction;
interpolation of f (x) with a quadratic function:
when the quadratic interpolation shape function fitting is adopted in the units, three node coordinates in any unit are respectively set as x1, x2, x3 and x2 as midpoints, so that x is satisfied 1 +x 3 ＝2x 2 Each node has a value f (x 1 )、f(x 2 )、f(x 3 ) F (x) is expressed by using a quadratic function:
f(x)＝N 1 f(x 1 )+N 2 f(x 2 )+N 3 f(x 3 ) (6)
wherein,,
the above formula (5) is written as:
order the
For the intra-cell fourier transform node coefficients, then equation (8) above is abbreviated as:
when the wave number is k x When the value is not 0, substituting the above formula (7) into W 1 、W 2 、W 3 Obtaining the intra-unit Fourier transform node coefficient:
W 1 、W 2 、W 3 the integral kernel functions all include
It is at [ x ] 1 ,x 3 ]The upper unit integral analysis solution is:
thus, k can be obtained x W when it is not 0 1 、W 2 、W 3 The semi-analytical solution is:
when the wave number is k x When the number of the organic light emitting diode is 0,
simple integration is performed to obtain zeroThe fourier transform node coefficients at wavenumber are:
and accumulating the analysis expressions of the different units to obtain a final one-dimensional Fourier positive transformation result.
5. A 3D AS-FT based three-dimensional gravity field numerical simulation method according to claim 3, characterized in that the wave number domain gravity bits are obtained by arbitrary sampling three-dimensional fourier positive transformation
The equation is satisfied:
wherein,,
for the residual density distribution of the wave number domain, according to the formula (14), solving to obtain the gravity position of the wave number domain
6. The 3D AS-FT based three-dimensional gravity field numerical simulation method according to claim 4, wherein the wave number domain gravity bit
And wave number domain gravity field->
The relation of (2) is:
according to the above (15), solving to obtain the wave number domain gravity field
In milligamma.
7. The 3D AS-FT based three-dimensional gravity field numerical simulation method according to claim 6, wherein the fourier transform formula of three-dimensional arbitrary sampling is:
k in x 、k y 、k z Representing wavenumber, F (x, y, z) is a spatial domain function, F (k) x ,k y ,k z ) Representing the wavenumber spectrum.
8. The 3D AS-FT based three-dimensional gravity field numerical simulation method according to any one of claims 1-7, wherein the three-dimensional poisson equation is:
wherein U is g Is space domain gravitational potential, pi is circumference rate, G is gravitational constant, ρ is space domain residual density, and unit is kg/m 3 。
9. The 3D AS-FT based three-dimensional gravity field numerical simulation method according to any one of claims 1-7, wherein the spatial domain of the target model is split in any one of the following ways to construct a three-dimensional spatial numerical model of the target region,
mode one: non-uniform subdivision is adopted for a preset first area, and encryption is carried out; wherein the first region satisfies the following formula:
ρ i for the residual density ρ of the corresponding node of the first region j For the remaining density of the jth node around the corresponding node of the first region, n 1 The number of nodes around the first area; omega is weight, and the value range is (0, 1);
mode two: sparse sampling is carried out on a preset second area, wherein the second area meets the following formula:
ρ i′ for the remaining density of the corresponding node of the second region, ρ j′ For the remaining density of the jth node around the corresponding node of the second region, n 2 The number of nodes around the second area; omega is weight, and the value range is (0, 1);
mode three: 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, characterized in that the processor implements the steps of the method of any of the preceding claims 1 to 9 when the computer program is executed.
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 CN115659579A (en) | 2023-01-31 |
| CN115659579B true 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) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116187107B (en) * | 2023-04-27 | 2023-08-11 | 中南大学 | Three-dimensional ground temperature field dynamic numerical simulation method, equipment and medium |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111505596A (en) * | 2020-04-16 | 2020-08-07 | 北京理工大学重庆创新中心 | Three-dimensional wind field inversion method based on non-uniform sampling correction VAD technology |
| CN113779498A (en) * | 2021-08-03 | 2021-12-10 | 中国电子产品可靠性与环境试验研究所（（工业和信息化部电子第五研究所）（中国赛宝实验室）） | Discrete Fourier matrix reconstruction method, device, equipment and storage medium |
Family Cites Families (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9964653B2 (en) * | 2011-12-21 | 2018-05-08 | Technoimaging, Llc | Method of terrain correction for potential field geophysical survey data |
| CN109100804B (en) * | 2018-06-11 | 2020-10-09 | 中国石油天然气集团有限公司 | Data reconstruction method, device and system for improving seismic data space sampling attribute |
| CN108984939A (en) * | 2018-07-30 | 2018-12-11 | 中南大学 | Three-dimensional Gravity field of force forward modeling method based on 3D Gauss-FFT |
| CN111967169B (en) * | 2020-10-21 | 2021-01-01 | 中南大学 | Two-degree body weight abnormal product decomposition numerical simulation method and device |
| CN112287534B (en) * | 2020-10-21 | 2022-05-13 | 中南大学 | NUFFT-based two-dimensional magnetic anomaly fast forward modeling method and device |
| CN112800657B (en) * | 2021-04-15 | 2021-06-18 | 中南大学 | 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 |
-
2022
- 2022-09-05 CN CN202211076505.7A patent/CN115659579B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111505596A (en) * | 2020-04-16 | 2020-08-07 | 北京理工大学重庆创新中心 | Three-dimensional wind field inversion method based on non-uniform sampling correction VAD technology |
| CN113779498A (en) * | 2021-08-03 | 2021-12-10 | 中国电子产品可靠性与环境试验研究所（（工业和信息化部电子第五研究所）（中国赛宝实验室）） | Discrete Fourier matrix reconstruction method, device, equipment and storage medium |
Non-Patent Citations (1)
| Title |
|---|
| 空间―波数混合域二度体重力异常正演;戴世坤;王旭龙;赵东东;刘志伟;张钱江;孙金飞;;石油地球物理勘探(06);18-19+223-229 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115659579A (en) | 2023-01-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Sharan et al. | Application of the multiquadric method for numerical solution of elliptic partial differential equations | |
| CN112287534B (en) | NUFFT-based two-dimensional magnetic anomaly fast forward modeling method and device | |
| CN109856689B (en) | Noise suppression processing method and system for superconducting aeromagnetic gradient tensor data | |
| CN115292973B (en) | Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system | |
| CN115659579B (en) | Three-dimensional gravity field numerical simulation method and system based on 3D AS-FT | |
| van Driel et al. | Accelerating numerical wave propagation using wavefield adapted meshes. Part I: forward and adjoint modelling | |
| CN114065586B (en) | Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method | |
| Roux et al. | Analysis of numerically induced oscillations in two-dimensional finite-element shallow-water models part II: Free planetary waves | |
| CN111967169B (en) | Two-degree body weight abnormal product decomposition numerical simulation method and device | |
| CN113608262B (en) | Seismic data processing method and device for calculating rotation component by using translation component | |
| CN114167511B (en) | Bit field data rapid inversion method based on continuous expansion downward continuation | |
| Papadakis et al. | Spatial filtering in a 6D hybrid-Vlasov scheme to alleviate adaptive mesh refinement artifacts: a case study with Vlasiator (versions 5.0, 5.1, and 5.2. 1) | |
| JPWO2019171660A1 (en) | Magnetic field analyzer, analysis method, and program | |
| Hou et al. | 3D inversion of vertical gravity gradient with multiple graphics processing units based on matrix compression | |
| CN106353801B (en) | The three-dimensional domain Laplace ACOUSTIC WAVE EQUATION method for numerical simulation and device | |
| CN116244877B (en) | Three-dimensional magnetic field numerical simulation method and system based on 3D Fourier transform | |
| WO2023280122A1 (en) | Density determination method and apparatus, and electronic device | |
| CN114325870A (en) | Method and system for calculating potential field gradient tensor based on cubic spline function | |
| CN113267830A (en) | Two-dimensional gravity gradient and seismic data joint inversion method based on non-structural grid | |
| Guo et al. | Microwave inversion for sparse data using descent learning technique | |
| Franklin et al. | A high-order fast marching scheme for the linearized eikonal equation | |
| CN115795231B (en) | Space wave number mixed domain three-dimensional strong magnetic field iteration method and system | |
| CN111665550A (en) | Underground medium density information inversion method | |
| Kupiainen et al. | GPU-acceleration of a high order finite difference code using curvilinear coordinates | |
| CN107562690A (en) | A kind of AIM PO mixed methods based on GPU |
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 |