CN111400654B - Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix - Google Patents

Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix Download PDF

Info

Publication number
CN111400654B
CN111400654B CN202010174859.XA CN202010174859A CN111400654B CN 111400654 B CN111400654 B CN 111400654B CN 202010174859 A CN202010174859 A CN 202010174859A CN 111400654 B CN111400654 B CN 111400654B
Authority
CN
China
Prior art keywords
matrix
vector
expression
rho
sub
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
Application number
CN202010174859.XA
Other languages
Chinese (zh)
Other versions
CN111400654A (en
Inventor
赵广东
柳建新
陈波
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202010174859.XA priority Critical patent/CN111400654B/en
Publication of CN111400654A publication Critical patent/CN111400654A/en
Application granted granted Critical
Publication of CN111400654B publication Critical patent/CN111400654B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computing Systems (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Discrete Mathematics (AREA)
  • Complex Calculations (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention provides a three-dimensional gravity field fast forward modeling method under a spherical coordinate system, which is characterized in that a forward modeling kernel matrix is divided into a plurality of Toplitz sub-matrixes, then the blocks are quickly calculated by adopting frequency domain FFT, and finally, summation is carried out, time-consuming matrix-vector product operation is converted into summation operation, so that the calculation time of the kernel matrix in forward modeling is reduced, the repetitive calculation efficiency of multiplication of the kernel matrix and a density vector is improved, and the overall efficiency of forward modeling calculation is improved; the method adopts a special subdivision mode, does not need to calculate all kernel matrix elements, only needs to calculate the kernel matrix elements generated by the first column of tesseroid unit bodies of the underground field source, and maps the values of other elements through the equivalent relation, thereby greatly reducing the calculation cost of the kernel matrix and reducing the occupation of the memory. The invention also provides an inversion method, which is combined with the forward modeling method, does not need to store a large Jacobian matrix, only needs to store a final matrix, and reduces the calculation and storage of intermediate variables.

Description

Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix
Technical Field
The invention relates to the technical field of geophysical and exploration, in particular to a gravity field fast forward modeling method and an inversion method based on a Toplitz nuclear matrix.
Background
With the popularization of global gravity satellite data with high precision and high coverage rate, the method provides possibility for exploration of gravity exploration, deep geological structure, shell mantle deformation, moon mars and other stars.
The gravity field forward inversion is an important technical means and is widely applied to underground three-dimensional density imaging research. The gravity field forward and backward inversion can be carried out at a conventional exploration scale, is used for research of resource exploration (mineral products, petroleum and natural gas) and hydrological environments and the like, is mostly carried out on the basis of a rectangular coordinate system, and common rapid algorithms comprise an FFT (fast Fourier transform algorithm), a Gauss-FFT (Gaussian-fast Fourier transform) method, a multi-stage expansion method and the like, so that the calculation efficiency and the precision of the inversion are ensured. However, when we study the gravity anomaly in a certain region or even in a global scale, the influence of the curvature of the earth must be considered, and the conventional method for forward and backward modeling of the gravity field in the rectangular coordinate system is no longer applicable. The corresponding data processing, the forward algorithm and the inversion are carried out under the spherical coordinate system.
Since the spherical three-dimensional integral has no analytic solution, discrete solution must be performed by a numerical method, and common numerical methods include a two-dimensional gaussian legendre integral method, a three-dimensional gaussian legendre integral method, a taylor series expansion method and the like. The traditional methods have low computational efficiency, are difficult to adapt to the requirements of large-scale three-dimensional inversion, and forward computational efficiency directly restricts large-scale high-resolution gravity field inversion. Under the background that the computing performance of the current computer is hard to be qualitatively leap, only a new forward modeling method can be searched for in order to greatly improve the computing efficiency. Secondly, in the traditional three-dimensional gravity field inversion, the memory occupation of a large Jacobian matrix is a problem which seriously restricts the inversion resolution.
On the other hand, since the 21 st century, the development of satellite gravity measurement technology has revolutionized not only the research of the earth gravitational field, but also the research of detecting the internal structure and material migration of the earth by using the information of the earth gravitational field, and the research of detecting the mass tumor basin and the lunar valance elevation inside the moon by using the GRAIL satellite. The increasing abundance of satellite gravity (including satellite altimetry), marine gravity, land gravity and GPS measurement data greatly improves the resolution and resolution of the earth gravitational field. The global-coverage high-resolution and high-precision gravity data creates conditions for obtaining a regional and global mantle high-resolution density disturbance model, is also important basic data for solid earth science research, and provides a new opportunity for global large-scale three-dimensional satellite gravity inversion. The conventional gravitational field forward and backward evolution method has the basic requirement of difficult use of mass high-resolution data, and a set of efficient high-precision forward and backward evolution method is urgently needed to be developed.
Disclosure of Invention
The invention provides a gravity field fast forward modeling method based on a Toplitz nuclear matrix (Toplitz nuclear matrix), which enables complex calculation of multiplication of a gravity field nuclear matrix and a density vector to be converted into a frequency domain, improves the calculation efficiency by three orders of magnitude, reduces the storage of a large sparse matrix, is beneficial to large-scale high-resolution gravity field inversion, and has the following specific technical scheme:
a gravity field fast forward modeling method based on a Toplitz nuclear matrix comprises the following steps:
step one, dividing an underground three-dimensional field source into N equally spaced in the longitude direction, the latitude direction and the depth directionλ、
Figure BDA0002410461060000021
And NrSegment, in total
Figure BDA0002410461060000022
Each tesseoid unit body has a division pitch of delta lambda,
Figure BDA0002410461060000023
And Δ r; the observation surface is positioned right above the three-dimensional field source, the position distribution of the observation points corresponds to the central point positions of the tesseroiid unit bodies divided by the three-dimensional field source one by one, and the number of the observation points is counted
Figure BDA0002410461060000024
Step two, expressing the matrix-vector product of the gravity anomaly generated by the underground three-dimensional field source tesseloid unit body on the observation point by adopting an expression 1):
K·ρ＝g 1)；
wherein: k is a gravity anomaly kernel matrix, rho is a density vector of an underground tesseoid unit body, and g is a gravity anomaly vector on an observation point;
step three, obtaining a sub-matrix of the gravity anomaly kernel matrix K, a sub-vector of the density vector rho and a sub-vector of the observation gravity anomaly vector g;
the specific sub-matrix for obtaining the gravity anomaly kernel matrix K is as follows:
partitioning a kernel matrix K into
Figure BDA0002410461060000025
Sub-matrices, each sub-matrix having a size of Nλ×NλAs in expression 4):
Figure BDA0002410461060000026
wherein: each sub-matrix K of the kernel matrix Ki,jAre both toeplitz matrices, as shown in expression 5),
Figure BDA0002410461060000027
Figure BDA0002410461060000028
Figure BDA0002410461060000031
definition ki,j-toeplitz (t), denotes ki,jIs a Toplitz matrix, t ═ t1-n,…,t-1,t0,t1,…,tn-1) Is the eigenvector of the Topritz matrix and has ta＝t-a，a＝0,1,…,n-1，n＝NλIs the number of elements of the feature vector t;
dividing the density vector rho and the observed gravity anomaly vector g into subvectors, such as expression 6):
Figure BDA0002410461060000032
wherein: each subvector ρ of the density vector ρjAnd observing each subvector g of the gravity anomaly vector giAre all NλAnd is and
Figure BDA0002410461060000033
and
Figure BDA0002410461060000034
step four, finding out the submatrix K of the kernel matrix Ki,jSubvectors rho of corresponding density vectors rhojAnd a sub-vector g for observing the gravity anomaly vector giAs shown in expression 7):
Figure BDA0002410461060000035
the sub-matrix K of the kernel matrix Ki,jSubvectors rho corresponding to density vectors rhojThe product of (d) is written as expression 8):
ki,jρj＝gi 8)；
step five, calling a fast Fourier transform algorithm to calculate and obtain a sub-vector g of the observation gravity anomaly vector giThe method specifically comprises the following steps:
implementing the sub-matrix k by invoking a fast Fourier transform algorithmi,jSubvectors rho of the density vector rhojThe convolution operation of (1) is as in expression 11):
Figure BDA0002410461060000041
wherein:
Figure BDA0002410461060000042
submatrix K being a kernel matrix Ki,jThe size of the spreading matrix of (2 n × 2 n); t is textIs an extension vector of the vector t, with a length of 2 n;
Figure BDA0002410461060000043
and
Figure BDA0002410461060000044
is a sub-vector g of the observed gravity anomaly vector giAnd a subvector of the density vector rhojThe extension vectors of (1) are all 2n in length, the first n elements of the extension vectors and giAnd ρjThe same, the last n elements are 0, 0 represents zero vector and the length is n; s is an intermediate variable matrix used for storing a temporary calculation result, and the size of the intermediate variable matrix is n multiplied by n;
Figure BDA0002410461060000045
and
Figure BDA0002410461060000046
respectively representing one-dimensional inverse Fourier transform and one-dimensional forward Fourier transform;
will be provided with
Figure BDA0002410461060000047
Taking the first n items to obtain a sub-vector g of the observed gravity anomaly vector gi；
Taking i as i +1 or j as j +1, and if i is less than or equal to i
Figure BDA0002410461060000048
Or j is less than or equal to
Figure BDA0002410461060000049
Returning to the fourth step; obtaining each submatrix K of the kernel matrix Ki,jCorresponding gi；
Step seven, g obtained in step sixiArranging in the order of expression 6) to obtain an observed gravity anomaly vector g.
Preferably, in the above technical solution, the obtaining of the submatrix of the gravity anomaly kernel matrix K in the third step is to calculate to obtain the gravity anomaly kernel matrix K, specifically:
adopting an expression 2) to calculate a gravity anomaly kernel matrix K:
Figure BDA00024104610600000410
wherein: r is0Is the radial component of the observation points numbered p and q,
Figure BDA00024104610600000411
is the latitude component of the observation point, λpIs the longitude component of the observation point, p is 1,2, …, Nλ，
Figure BDA00024104610600000412
rl′、
Figure BDA00024104610600000413
And λ'mThree components of the coordinates of the center point of tesseoid unit cells numbered l, N and m, l being 1,2, …, Nr，m＝1，
Figure BDA00024104610600000414
G is the gravitational constant; r
Figure BDA00024104610600000415
And λ' are integral variables in the radial, latitudinal and longitudinal directions, respectively, and the integration intervals are rl' -0.5 Δ r to rl′+0.5Δr、
Figure BDA00024104610600000416
To
Figure BDA00024104610600000417
And λ'm-0.5 Δ λ to λ'm+0.5Δλ；
Figure BDA00024104610600000418
The cos ψ is the cosine of the intermediate angle from the observation point to any point in the tesseroid unit body, and is calculated by adopting an expression 3) as the distance from the observation point to any point in the tesseroid unit body:
Figure BDA0002410461060000051
preferably, in the above technical solution, in the fourth step:
due to the Topritz matrix ki,jThe product of the matrix and the vector is a one-dimensional discrete convolution operation with the element gkWith expression 9):
Figure BDA0002410461060000052
wherein: k is 0,1, …, n-1; t is the Toplitz matrix ki,jCharacteristic vector of (d), tk-bIs the k-b element of the vector t; rhojIs a density vectorSubvectors of rho, rhobIs the density subvector ρjJ is 0,1, …, n-1, b is 0,1, …, n-1; n is the number of elements of the feature vector t;
will Toplitz matrix ki,jConstructed as a cyclic matrix
Figure BDA0002410461060000053
Where the n × n circulant matrix C has expression 10);
Figure BDA0002410461060000054
wherein: c ═ c0,c1,…,cn-1) A feature vector representing the circulant matrix C;
will Toplitz matrix ki,jConversion into a 2n x 2n circulant matrix
Figure BDA0002410461060000055
text＝(t0,t1,…,tn-1,0,t1-n,…,t-1) Representing a circulant matrix
Figure BDA0002410461060000056
The feature vector of (2).
The invention also provides a gravity field rapid inversion method based on the Toplitz nuclear matrix, which specifically comprises the following steps:
obtaining an inversion target function phi (rho) in a spherical coordinate system as an expression 13):
Figure BDA0002410461060000057
wherein: wdThe method comprises the following steps of (1) obtaining a diagonal matrix of observation data weight, wherein each diagonal element represents the credibility of the observation data; wρIs a large sparse model weight matrix; beta is a regularization coefficient used to weigh the specific gravity between the model objective function and the data objective function; rhorefTo representA reference model; k is a gravity anomaly kernel matrix, and rho is a density vector of an underground tesserioid unit body;
the following system of linear equations, as expressed in expression 14, is obtained by minimizing the inverse objective function Φ (ρ):
Figure BDA0002410461060000061
wherein, KTIs the transpose of K, Wd TIs WdTranspose of (W)ρ TIs WρTransposing;
substituting the observed gravity anomaly vector g obtained by the forward modeling method into expression 14) to obtain the density vector rho of the underground tesseoid unit body.
The technical scheme of the invention has the following beneficial effects:
1. the invention provides a novel three-dimensional gravity field fast forward modeling method under a spherical coordinate system, which divides a forward modeling kernel matrix into a plurality of Toplitz submatrices (Toplitz submatrices), then blocks are quickly computed by adopting frequency domain FFT, and finally, summation is accumulated, time-consuming matrix-vector product operation is converted into summation operation, so that the computation time of the kernel matrix in forward modeling is reduced by 2 orders of magnitude, the repetitive computation efficiency of multiplying the kernel matrix and a density vector is improved by more than 20 times, and the forward modeling computation efficiency is integrally improved by 3 orders of magnitude.
2. When the gravity field nuclear matrix is calculated, due to the adoption of a special subdivision mode, all nuclear matrix elements do not need to be calculated, only the nuclear matrix elements generated by the first column of tesseroid unit bodies of the underground field source are calculated, and the values of other elements can be mapped out through an equivalent relationship. The method greatly reduces the calculation cost of the kernel matrix and simultaneously reduces the memory occupation. In the inversion, the large Jacobian matrix does not need to be stored, and only the final matrix (beta W) needs to be storedρ T Wρρ and β Wρ T Wρρref) That is, the calculation and storage of intermediate variables are reduced.
3. In order to prove the correctness of the method provided by the invention, the model comprises 360 × 720 × 10 unit bodies and the number of observation points is 360 × 720, the forward modeling method provided by the invention is used in 1965 seconds, while the calculation time of the traditional method is 890964 seconds, and the calculation efficiency is improved by 432 times. In addition, the maximum relative error of the forward modeling method and the inversion method is 0.024%, and the effectiveness of the method provided by the invention is proved.
Then, the forward modeling method and the inversion method are applied to inversion of lunar mass nodules, and the inversion result is basically consistent with the depth information of the Mohuo surface obtained by the predecessor, so that the accuracy of the method is further proved.
The forward modeling and inversion method provided by the invention is suitable for resource exploration of exploration scale and hydrological environment monitoring, is also suitable for large-scale three-dimensional density imaging, researching earth crust mantle structure, structure movement, rock ring transition and the like, and has wide application range.
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 drawings.
Drawings
The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate an embodiment of the invention and, together with the description, serve to explain the invention and not to limit the invention. In the drawings:
FIG. 1 is a schematic diagram of a tesseroiid unit body subdivision in a spherical coordinate system;
FIG. 2 is a schematic diagram of the relationship of a layer of 3 × 4 tesseoid units to observation points and the equivalence relationship of a kernel matrix, where: (a) a kernel matrix and a submatrix thereof; (b) and (c) is a schematic diagram of the equivalence of the elements of the kernel matrix, wherein the values of the elements of the kernel matrix represented by the same line type are equal;
FIG. 3 is a graph of the maximum relative error with respect to the subdivision interval for the proposed method and the conventional method;
FIG. 4 is a graph of relative computation time with respect to subdivision intervals for the proposed method and the conventional method;
FIG. 5 is a process picture for processing lunar global gravity data according to a preferred embodiment of the present invention, wherein: (a) a lunar global topographic map; (b) free air gravity anomaly; (c) terrain correction; (d) the Booth gravity is abnormal, (d) Chinese and English words are the names of large mass tumor basins on the moon, and 6 sections are represented by black solid lines;
FIG. 6 is a density profile at various depths, where: (a) 5km, (b) 15km, (c) 25km, (d) 35km, (e) 45km, (f) 55km, (g) 65km, (h) 75km, (i) 85km, (j) 95 km;
FIG. 7 is the inverted density distribution of the 6 sections on (d) in FIG. 5, numbered 1-6; the solid bold black lines in the figure indicate the depth of the mojo face previously found.
Detailed Description
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, which are defined and covered by the claims.
The invention provides a gravity field fast forward modeling method (namely a nuclear matrix and density vector product calculation method) under a spherical coordinate system, which is applied to gravity field three-dimensional density inversion and has the following details:
the method for fast forward modeling of the gravity field in the spherical coordinate system based on the frequency domain FFT algorithm specifically comprises the following steps:
firstly, dividing an underground three-dimensional field source into a plurality of tesserioid unit bodies, wherein the details are as follows:
for exploration scale, the forward and backward rotation of the gravity field is usually performed in a rectangular coordinate system, however, when the research area is large or the whole earth rock circle is researched, the influence of the earth curvature has to be considered, and therefore, the forward and backward rotation of the gravity field needs to be performed in a spherical coordinate system. In a spherical coordinate system, the field sources are usually subdivided into tesseoid unit cells as shown in fig. 1. Dividing an underground three-dimensional field source into N equally spaced in the longitude direction, the latitude direction and the depth directionλ、
Figure BDA0002410461060000081
And NrSegment, in total
Figure BDA0002410461060000082
Each tesseoid unit body has a division pitch of delta lambda,
Figure BDA0002410461060000083
And Δ r. The observation surface is positioned right above the three-dimensional field source, the central point position of the observation point corresponds to the central point position of the dissected tesseoid unit body of the three-dimensional field source one by one, and the number of the observation points is counted
Figure BDA0002410461060000084
As shown in fig. 2.
Secondly, obtaining a matrix-vector product of gravity anomaly generated by the underground three-dimensional field source tesseloid unit body on an observation point, and specifically adopting an expression 1) to represent:
K·ρ＝g 1)；
where K is the gravity anomaly kernel matrix with dimension Nobs×Ntess(ii) a ρ is the density vector of the subsurface tesseoid unit cell, whose dimension is Ntess(ii) a g is the gravity anomaly vector at the observation point, with dimension Nobs。
Adopting an expression 2) to calculate a gravity anomaly kernel matrix K:
Figure BDA0002410461060000085
wherein: r is0Is the radial component of the observation points numbered p and q,
Figure BDA0002410461060000086
is the latitude component of the observation point, λpIs the longitude component of the observation point, p is 1,2, …, Nλ，
Figure BDA0002410461060000087
rl′、
Figure BDA0002410461060000088
And λ'mThree components of the coordinates of the center point of tesseoid unit cells numbered l, N and m, l being 1,2, …, Nr，m＝1，
Figure BDA0002410461060000089
G is the gravitational constant; r
Figure BDA00024104610600000810
And λ' are integral variables in the radial, latitudinal and longitudinal directions, respectively, and the integration intervals are rl' -0.5 Δ r to rl′+0.5Δr、
Figure BDA00024104610600000811
To
Figure BDA00024104610600000812
And λ'm-0.5 Δ λ to λ'm+0.5Δλ；
Figure BDA00024104610600000813
The cos ψ is the cosine of the intermediate angle from the observation point to any point in the tesseroid unit body, and is calculated by adopting an expression 3) as the distance from the observation point to any point in the tesseroid unit body:
Figure BDA00024104610600000814
the theoretical formula expression 2) has no analytic solution, and must be solved by a numerical method, a two-dimensional Gaussian legendre method, a three-dimensional Gaussian legendre method or a Taylor series expansion method can be selected, because the Taylor series expansion method has analytic singularities at poles, the two-dimensional Gaussian legendre method has complex calculation formula and consumes time, and the three-dimensional Gaussian legendre numerical method with 2 points in three directions is selected to solve the nuclear matrix. The more the number of Gaussian points is, the higher the calculation precision is, but the corresponding calculation time is multiplied. And 2-point Gaussian coefficients are selected, so that the high-precision forward modeling requirement can be met, and the maximum relative error is lower than 0.1%.
Thirdly, obtaining a sub-matrix of the gravity anomaly kernel matrix K, a sub-vector of the density vector rho and a sub-vector of the observation gravity anomaly vector g, wherein the details are as follows:
as can be seen from fig. 2, the kernel matrix K is a large partitioned Toplitze matrix, i.e. a large number of equal elements exist in the kernel matrix. Thus, when the kernel matrix is computed only, it is only necessary for all p to be 1,2, …, Nλ，
Figure BDA0002410461060000091
l＝1,2,…,Nr，m＝1，
Figure BDA0002410461060000092
And (4) calculating. Where m is taken to be only 1, this means that only the kernel matrix elements resulting from the first column of tesseoid unit cells need to be computed to reduce duplicate computations. By this strategy, N is reducedλThe calculation amount and the memory occupation amount are multiplied, and the strategy has more remarkable advantages in large-scale three-dimensional forward inversion.
Dividing the calculated kernel matrix K into
Figure BDA0002410461060000093
Sub-matrices, each sub-matrix having a size of Nλ×NλExpressed by expression 4):
Figure BDA0002410461060000094
wherein: each sub-matrix K of the kernel matrix Ki,jAre both toplitz matrices, as in expression 5),
Figure BDA0002410461060000095
Figure BDA0002410461060000096
Figure BDA0002410461060000097
definition ki,j-toeplitz (t), denotes ki,jIs a Toplitz matrix, t ═ t1-n,…,t-1,t0,t1,…,tn-1) Is the eigenvector of the Topritz matrix and has ta＝t-aFor a ═ 0,1, …, N-1, and N ═ NλIs the number of elements of the feature vector t;
dividing the density vector rho and the observation gravity anomaly vector g into subvectors as shown in expression 6):
Figure BDA0002410461060000101
wherein: each subvector ρ of the density vector ρjAnd observing each subvector g of the gravity anomaly vector giAre all NλAnd is and
Figure BDA0002410461060000102
and
Figure BDA0002410461060000103
the fourth step is to find out the sub-matrix K of the kernel matrix Ki,jSubvectors rho of corresponding density vectors rhojAnd a sub-vector g for observing the gravity anomaly vector giAs in expression 7):
Figure BDA0002410461060000104
then the kernel matrix submatrix ki,jAnd density vector subvector ρjThe product of (d) is written as expression 8):
ki,jρj＝gi 8)。
due to the Topritz matrix ki,jThe product of the matrix and the vector is a one-dimensional discrete convolution operation with the element gkWith expression 9):
Figure BDA0002410461060000105
wherein: k is 0,1, …, n-1; t is the Toplitz matrix ki,jCharacteristic vector of (d), tk-bIs the k-b element of the vector t; rhojIs a sub-vector of the density vector p, pbIs the density subvector ρjJ is 0,1, …, n-1, b is 0,1, …, n-1; n is the number of elements of the feature vector t;
will Toplitz matrix ki,jConstructed as a cyclic matrix
Figure BDA0002410461060000111
Where n × n circulant (C) is expression 10):
Figure BDA0002410461060000112
wherein: c ═ c0,c1,…,cn-1) A feature vector representing the circulant matrix C;
will Toplitz matrix ki,jConversion into a 2n x 2n circulant matrix
Figure BDA0002410461060000113
text＝(t0,t1,…,tn-1,0,t1-n,…,t-1) Representing a circulant matrix
Figure BDA0002410461060000114
The feature vector of (2).
Fifthly, realizing the sub-matrix k by calling a fast Fourier transform algorithmi,jSubvectors rho of the density vector rhojThe convolution operation of (1) is as in expression 11):
Figure BDA0002410461060000115
wherein:
Figure BDA0002410461060000116
submatrix K being a kernel matrix Ki,jThe size of the spreading matrix of (2 n × 2 n); vector textIs an extension vector of the vector t, with a length of 2 n;
Figure BDA0002410461060000117
and
Figure BDA0002410461060000118
respectively, a subvector g for the observed gravity anomaly vector giAnd a subvector of the density vector rhojThe extension vectors of (1) are all 2n in length, the first n elements of the extension vectors and giAnd ρjThe same, the last n elements are 0; 0 represents a zero vector with a length of n; s is an intermediate variable matrix used for storing a temporary calculation result, and the size of the intermediate variable matrix is n multiplied by n;
Figure BDA0002410461060000119
and
Figure BDA00024104610600001110
respectively representing one-dimensional inverse Fourier transform and one-dimensional forward Fourier transform;
will be provided with
Figure BDA00024104610600001111
Taking the first n items to obtain a subvector gi。
So far, the product of the sub-matrix and the density vector is quickly calculated by the frequency domain FFT algorithm, and then the product of each sub-matrix and the corresponding density sub-vector only needs to be calculated (specifically, after the fourth step to the fifth step are completed, i ═ i +1 or j ═ j +1 is taken, if i is less than or equal to i ≦ j ≦ 1
Figure BDA00024104610600001112
Or j is less than or equal to
Figure BDA00024104610600001113
Returning to the fourth step), and then arranging the vectors in the order of expression 6) to obtain an observed gravity anomaly vector g. Calculating the product of each nuclear matrix sub-matrix and the corresponding density sub-vector, and accumulating and summing to obtain the gravity anomaly sub-vector giUsing the following expression 12) for transportationCalculating:
Figure BDA0002410461060000121
a method for quickly inverting a gravity field in a spherical coordinate system based on a frequency domain FFT algorithm specifically comprises the following steps:
constructing an inversion target function phi (rho) as an expression 13 under a spherical coordinate system by using a Gihonov regularization inversion method:
Figure BDA0002410461060000122
wherein: wdThe method comprises the following steps of (1) obtaining a diagonal matrix of observation data weight, wherein each diagonal element represents the credibility of the observation data; wρIs a large sparse model weight matrix; beta is a regularization coefficient used to weigh the specific gravity between the model objective function and the data objective function; rhorefRepresenting a reference model, typically taking p ref0; k is the gravity anomaly kernel matrix and ρ is the density vector of the subsurface tesserioid unit volume.
By minimizing the inverse objective function, a system of linear equations like expression 14) can be obtained:
Figure BDA0002410461060000123
wherein: kTIs the transpose of K, Wd TIs WdTranspose of (W)ρ TIs WρThe transposing of (1). The present invention solves this system of equations using a well-established conjugate gradient method.
It is noted that the above fast forward method for the gravity field in the spherical coordinate system based on the topelitz (Toplitze) kernel matrix frequency domain FFT algorithm can be directly used for calculating K · ρ in expression 14). Due to KTAlso a Toplitz matrix, so that W can be calculated firstd TWdK.rho and Wd TWdg, obtaining the bidirectional quantityThen, K is calculated by utilizing the gravity field fast forward modeling method under the spherical coordinate system based on the Toplitz (Toplitz) kernel matrix frequency domain FFT algorithmTThe product of these two vectors. And without storing the large matrix Wρ TAnd WρOnly the final matrix beta W needs to be storedρ T Wρρ and β Wρ T WρρrefAnd large-scale high-resolution inversion is possible.
When the rapid inversion method based on the Toplitz matrix in the spherical coordinate system is utilized, the specific implementation steps are as follows: firstly, inputting observation data from a gridding grid file, automatically judging the number and the initial range of the observation data by a program, and calculating subdivision intervals, coordinate information of each measuring point and the like; then setting an observation height and setting an error mean square error of observation data; then, inversion parameters are set, and a program needs to input an initial range of inversion depth, wherein the range can be obtained according to geology, regional structure and previous research results of the region; and finally, running a program, and operating according to the formula to obtain a final inversion result.
Example 1:
the present embodiment uses a lunar spherical shell forward model to verify the correctness and efficiency of the method proposed by the present invention, since only the spherical shell model has an analytic solution. The thickness of the spherical shell layer is 100km, the initial range is 1638 km to 1738km (the actual radius of the moon), and the density of the spherical shell is set to be 1000kg/m3The observation surface degree is a spherical surface which is 10km away from the outer spherical shell, and the spherical shell layer is divided into 10 sections in the radial direction at intervals of 10 km. The interval between the warp direction and the weft direction is 0.5-10 degrees, and the relationship between the maximum relative error and the calculation time along with the subdivision interval is counted.
The relationship between the maximum relative error and the subdivision interval is shown in fig. 3, and the relationship between the relative calculation time and the subdivision interval is shown in fig. 4. In order to more intuitively feel the high efficiency of the method provided by the invention, table 1 counts the maximum relative error when the subdivision interval is 0.5 ° (i.e. the number of subdivision units is 360 × 720 × 10), the calculation time of the kernel matrix, the calculation time of the product of the kernel matrix and the density vector, and the comparison of the total time.
In the forward modeling test, when the number of the subdivision units is 360 × 720 × 10, the calculation time based on the method provided by the invention is 1965 seconds and about 0.5 hour, while the calculation time based on the traditional forward modeling method is 890964 seconds and about 248 hours, and the calculation efficiency is improved by about 500 times. In addition, the method provided by the invention can also ensure high precision, and the maximum relative error of the method is the same as that of the traditional method and is lower than 0.1 percent.
TABLE 1 comparison of maximum relative error, time of kernel matrix calculation, time of kernel matrix multiplied by density vector, and total time
Figure BDA0002410461060000131
Example 2:
in order to further prove the effectiveness and high efficiency of the rapid forward modeling method for the gravity field of the spherical coordinate system based on the Toplitz matrix frequency domain FFT algorithm, which is provided by the application of the invention, next, an actual data inversion case is provided in the embodiment, and the proposed method is applied to lunar global scale three-dimensional density imaging research, and details are as follows:
the lunar global terrain is solved by a lunar terrain harmonic model LRO _ LTM05_2050, which is the latest lunar terrain harmonic model, as shown in fig. 5 (a).
The present invention uses the lunar latest gravity field model GL1500E to calculate the free air gravity anomaly, as shown in fig. 5 (b), and the terrain correction results are shown in fig. 5 (c). The bogey anomaly is obtained by subtracting the terrain correction (fig. 5 (c)) from the free air gravity anomaly (fig. 5 (b)), as shown in fig. 5 (d). Here, the density used for terrain correction was 2560kg/m3The observation height was 10 km. With the first 180 steps of the model, corresponding to a spatial resolution of 1.0.
The invention divides the underground field source into 180 x 360 tesseoid unit bodies with equal intervals, the model range is set to be 0km to 100km in the depth direction according to the research foundation of the predecessor, the tesseoid unit bodies are divided into 10 sections with the interval of 10km, the total number of the tesseoid unit bodies is 180 x 360 x 10, and the number of observation points is 180 x 360. The three-dimensional gravity field inversion is carried out by using the Booth gravity anomaly, and the density distribution of anomaly at each depth is shown in figure 6. The density distribution along the 6 sections shown in (d) in fig. 5 is shown in fig. 7. It can be seen that the method provided by the invention can correctly invert the positions and depths of the abnormal bodies, and the high-density abnormal bodies are mainly distributed at the bottoms of the mass basins, and the distribution range of the high-density abnormal bodies in the depths is about 15-60 km. In addition, the top interface of these high density anomalies is consistent with the mohuo face relief depth previously studied.
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 (4)

1. A gravity field fast forward modeling method based on a Toplitz nuclear matrix is characterized by comprising the following steps:
step one, dividing an underground three-dimensional field source into N equally spaced in the longitude direction, the latitude direction and the depth directionλ、
Figure FDA00034760993300000111
And NrSegment, in total
Figure FDA0003476099330000012
Each tesseoid unit body has a division pitch of delta lambda,
Figure FDA0003476099330000019
And Δ r; the observation surface is positioned right above the three-dimensional field source, the position distribution of the observation points corresponds to the central point positions of the tesseroiid unit bodies divided by the three-dimensional field source one by one, and the number of the observation points is counted
Figure FDA00034760993300000110
A plurality of;
step two, expressing the matrix-vector product of the gravity anomaly generated by the underground three-dimensional field source tesseloid unit body on the observation point by adopting an expression 1):
K·ρ＝g 1)；
wherein: k is a gravity anomaly kernel matrix, rho is a density vector of an underground tesseroid unit body, and g is an observation gravity anomaly vector;
step three, obtaining a sub-matrix of the gravity anomaly kernel matrix K, a sub-vector of the density vector rho and a sub-vector of the observation gravity anomaly vector g;
the specific sub-matrix for obtaining the gravity anomaly kernel matrix K is as follows:
partitioning a kernel matrix K into
Figure FDA00034760993300000112
Sub-matrices, each sub-matrix having a size of Nλ×NλAs in expression 4):
Figure FDA0003476099330000016
wherein: each sub-matrix K of the kernel matrix Ki,jAre both toeplitz matrices, as shown in expression 5),
Figure FDA0003476099330000017
Figure FDA0003476099330000018
Figure FDA0003476099330000021
definition ki,j-toeplitz (t), denotes ki,jIs a Toplitz matrix, t ═ t1-n,…,t-1,t0,t1,…,tn-1) Is the eigenvector of the Topritz matrix and has ta＝t-a，a＝0,1,…,n-1，n＝NλIs the number of elements of the feature vector t;
dividing the density vector rho and the observed gravity anomaly vector g into subvectors, such as expression 6):
Figure FDA0003476099330000022
wherein: each subvector ρ of the density vector ρjAnd observing each subvector g of the gravity anomaly vector giAre all NλAnd is and
Figure FDA0003476099330000023
and
Figure FDA0003476099330000024
step four, finding out the submatrix K of the kernel matrix Ki,jSubvectors rho of corresponding density vectors rhojAnd a sub-vector g for observing the gravity anomaly vector giAs shown in expression 7):
Figure FDA0003476099330000025
the sub-matrix K of the kernel matrix Ki,jSubvectors rho corresponding to density vectors rhojThe product of (d) is written as expression 8):
ki,jρj＝gi 8)；
step five, calling a fast Fourier transform algorithm to calculate and obtain a sub-vector g of the observation gravity anomaly vector giThe method specifically comprises the following steps:
implementing the sub-matrix k by invoking a fast Fourier transform algorithmi,jSubvectors rho of the density vector rhojThe convolution operation of (1) is as in expression 11):
Figure FDA0003476099330000031
wherein:
Figure FDA0003476099330000032
submatrix K being a kernel matrix Ki,jThe size of the spreading matrix of (2 n × 2 n); t is textIs an extension vector of the vector t, with a length of 2 n; gi extAnd ρj extRespectively, a subvector g for the observed gravity anomaly vector giAnd a subvector of the density vector rhojThe extension vectors of (1) are all 2n in length, the first n elements of the extension vectors and giAnd ρjThe same, the last n elements are 0, 0 represents zero vector and the length is n; s is an intermediate variable matrix used for storing a temporary calculation result, and the size of the intermediate variable matrix is n multiplied by n;
Figure FDA00034760993300000312
and
Figure FDA00034760993300000311
respectively representing one-dimensional inverse Fourier transform and one-dimensional forward Fourier transform;
g is prepared fromi extTaking the first n items to obtain a sub-vector g of the observed gravity anomaly vector gi；
Taking i as i +1 or j as j +1, and if i is less than or equal to i
Figure FDA00034760993300000315
Or j is less than or equal to
Figure FDA00034760993300000314
Returning to the fourth step; obtaining each submatrix K of the kernel matrix Ki,jCorresponding gi；
Step seven, g obtained in step sixiArranging in the order of expression 6) to obtain an observed gravity anomaly vector g.
2. The gravity field fast forward modeling method based on Toplize kernel matrix according to claim 1, wherein the obtaining of the submatrix of the gravity anomaly kernel matrix K in the third step is to obtain the gravity anomaly kernel matrix K by calculation, specifically:
adopting an expression 2) to calculate a gravity anomaly kernel matrix K:
Figure FDA0003476099330000033
wherein: r is0Is the radial component of the observation points numbered p and q,
Figure FDA00034760993300000316
is the latitude component of the observation point, λpIs the longitude component of the observation point, p is 1,2, …, Nλ，
Figure FDA0003476099330000035
r′l、
Figure FDA0003476099330000036
And λ'mThree components of the coordinates of the center point of tesseoid unit cells numbered l, N and m, l being 1,2, …, Nr，m＝1，
Figure FDA0003476099330000037
G is the gravitational constant; r
Figure FDA0003476099330000038
And λ 'are integral variables in the radial direction, the latitude direction and the longitude direction, and the integral intervals are r'l-0.5 Δ r to r'l+0.5Δr、
Figure FDA0003476099330000039
To
Figure FDA00034760993300000310
And λ'm-0.5 Δ λ to λ'm+0.5 Δ λ; l is the distance from the observation point to any point in the tesseroid unit body, and cos ψ is the distance from the observation point to any point in the tesseroid unit bodyThe cosine of the intermediate angle is calculated by adopting an expression 3):
Figure FDA0003476099330000041
3. the Toplitz-kernel-matrix-based gravity field fast forward method according to claim 2, wherein in the fourth step:
due to the Topritz matrix ki,jThe product of the matrix and the vector is a one-dimensional discrete convolution operation with the element gkWith expression 9):
Figure FDA0003476099330000042
wherein: k is 0,1, …, n-1; t is the Toplitz matrix ki,jCharacteristic vector of (d), tk-bIs the k-b element of the vector t; rhojIs a sub-vector of the density vector p, pbIs the density subvector ρjJ is 0,1, …, n-1, b is 0,1, …, n-1; n is the number of elements of the feature vector t;
will Toplitz matrix ki,jConstructed as a cyclic matrix
Figure FDA0003476099330000043
Where n × n cyclic matrix C is circular (C) expression 10);
Figure FDA0003476099330000044
wherein: c ═ c0,c1,…,cn-1) A feature vector representing the circulant matrix C;
will Toplitz matrix ki,jConversion into a 2n x 2n circulant matrix
Figure FDA0003476099330000045
text＝(t0,t1,…,tn-1,0,t1-n,…,t-1) Representing a circulant matrix
Figure FDA0003476099330000046
The feature vector of (2).
4. A gravity field fast inversion method based on a Toplitz nuclear matrix is characterized by comprising the following steps:
obtaining an inversion target function phi (rho) in a spherical coordinate system as an expression 13):
Figure FDA0003476099330000047
wherein: wdThe method comprises the following steps of (1) obtaining a diagonal matrix of observation data weight, wherein each diagonal element represents the credibility of observation data; wρIs a large sparse model weight matrix; beta is a regularization coefficient used to weigh the specific gravity between the model objective function and the data objective function; rhorefRepresenting a reference model; g is an observation gravity anomaly vector; k is a gravity anomaly kernel matrix, and rho is a density vector of an underground tesserioid unit body;
the following system of linear equations, as expressed in expression 14, is obtained by minimizing the inverse objective function Φ (ρ):
Figure FDA0003476099330000051
wherein, KTIs the transpose of K, Wd TIs WdTranspose of (W)ρ TIs WρTransposing;
substituting the observed gravity anomaly vector g obtained in any one of claims 1-3 into expression 14) to obtain a density vector p for the subsurface tesserioid unit cell.
CN202010174859.XA 2020-03-13 2020-03-13 Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix Active CN111400654B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010174859.XA CN111400654B (en) 2020-03-13 2020-03-13 Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010174859.XA CN111400654B (en) 2020-03-13 2020-03-13 Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix

Publications (2)

Publication Number Publication Date
CN111400654A CN111400654A (en) 2020-07-10
CN111400654B true CN111400654B (en) 2022-03-18

Family

ID=71428778

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010174859.XA Active CN111400654B (en) 2020-03-13 2020-03-13 Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix

Country Status (1)

Country Link
CN (1) CN111400654B (en)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112596106B (en) * 2020-11-09 2024-04-26 中国人民武装警察部队黄金第五支队 Method for carrying out gravity and vibration joint inversion on density interface distribution under spherical coordinate system
CN113204057B (en) * 2021-05-07 2022-07-12 湖南科技大学 Multilayer-method-based gravity-magnetic fast inversion method
CN113238284B (en) * 2021-05-07 2022-09-27 湖南科技大学 Gravity and magnetic fast forward modeling method
CN113109883B (en) * 2021-05-26 2023-01-03 湖南科技大学 Method for forward modeling satellite gravity field based on isoparametric transformation global discrete grid spherical coordinates
CN113311494B (en) * 2021-05-26 2022-04-26 中南大学 Satellite gravity field inversion method
CN113740915B (en) * 2021-07-28 2024-04-19 中国人民武装警察部队黄金第五支队 Method for jointly inverting crust structure parameters by spherical coordinate system gravity and receiving function
CN113673163B (en) * 2021-08-25 2023-09-12 中南大学 Three-dimensional magnetic abnormal constant rapid forward modeling method, device and computer equipment
CN113420487B (en) * 2021-08-25 2021-10-29 中南大学 Three-dimensional gravity gradient tensor numerical simulation method, device, equipment and medium

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130258810A1 (en) * 2012-03-30 2013-10-03 Wenyi Hu Method and System for Tomographic Inversion
CN110045432B (en) * 2018-12-05 2020-06-30 中南大学 Gravity field forward modeling method and three-dimensional inversion method under spherical coordinate system based on 3D-GLQ
CN109375280B (en) * 2018-12-10 2020-03-31 中南大学 Gravity field rapid high-precision forward modeling method under spherical coordinate system

Also Published As

Publication number Publication date
CN111400654A (en) 2020-07-10

Similar Documents

Publication Publication Date Title
CN111400654B (en) Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix
CN110045432B (en) Gravity field forward modeling method and three-dimensional inversion method under spherical coordinate system based on 3D-GLQ
Weiss et al. High‐resolution surface velocities and strain for Anatolia from Sentinel‐1 InSAR and GNSS data
Golub et al. Large-scale geodetic least-squares adjustment by dissection and orthogonal decomposition
CN109375280B (en) Gravity field rapid high-precision forward modeling method under spherical coordinate system
Hill et al. Combination of geodetic observations and models for glacial isostatic adjustment fields in Fennoscandia
Robbins et al. A new global database of Mars impact craters≥ 1 km: 1. Database creation, properties, and parameters
Zhao et al. Efficient 3‐D large‐scale forward modeling and inversion of gravitational fields in spherical coordinates with application to lunar mascons
Boy et al. A comparison of tidal ocean loading models using superconducting gravimeter data
Doganalp et al. Local geoid determination in strip area projects by using polynomials, least-squares collocation and radial basis functions
Ophaug et al. A comparative assessment of coastal mean dynamic topography in N orway by geodetic and ocean approaches
Roussel et al. Complete gravity field of an ellipsoidal prism by Gauss–Legendre quadrature
Li et al. The contribution of the GRAV‐D airborne gravity to geoid determination in the Great Lakes region
Işık et al. High-resolution geoid modeling using least squares modification of Stokes and Hotine formulas in Colorado
Saraswati et al. New analytical solution and associated software for computing full-tensor gravitational field due to irregularly shaped bodies
Shen et al. Improved geoid determination based on the shallow-layer method: a case study using EGM08 and CRUST2. 0 in the Xinjiang and Tibetan regions
Grafarend et al. Geodesy-the Challenge of the 3rd Millennium
Tsoulis et al. A spectral assessment review of current satellite‐only and combined Earth gravity models
Haji-Aghajany et al. Assessment of InSAR tropospheric signal correction methods
Shi et al. Inferring decelerated land subsidence and groundwater storage dynamics in Tianjin–Langfang using Sentinel-1 InSAR
Rapp et al. Determination of the geoid: present and future
Fašková et al. Finite element method for solving geodetic boundary value problems
Fadil et al. New Zealand 20th century sea level rise: Resolving the vertical land motion using space geodetic and geological data
Agoshkov et al. Variational data assimilation technique in mathematical modeling of ocean dynamics
Coss et al. Channel water storage anomaly: A new remotely sensed quantity for global river analysis

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