CN109444973A - Gravity forward modeling accelerated method under a kind of spherical coordinate system - Google Patents

Abstract

The invention discloses a gravity forward modeling acceleration method under a spherical coordinate system, comprising: dividing a global three-dimensional model into a series of non-overlapping spherical shell elements; establishing corresponding octree structures for observation points and spherical shell elements respectively; For each observation point octree cell, divide the spherical shell unit octree (source octree) cells into the near area and the far area; calculate the gravity parameters of the near area and the far area separately; based on the spherical shell elements contained in the near area Using the Gaussian numerical integral formula of , calculates the near-area integral, and obtains the near-area gravity parameters; and based on the observation point octree and the source octree, adopts the adaptive multi-layer fast multipole algorithm to calculate the far-area integral, and obtains the far-area gravity parameter; The calculated gravity parameters of the near area and the far area are summed to obtain the gravity parameters of the observation point. The invention can carry out the gravity forward calculation for the satellite gravity observation points distributed arbitrarily based on the division of the spherical shell unit, and greatly improves the efficiency of the satellite gravity forward calculation.

Description

A Gravity Forward Acceleration Method in Spherical Coordinate System

technical field

The invention belongs to the technical field of geophysics, and in particular relates to a gravity forward acceleration method in a spherical coordinate system.

Background technique

Since Germany launched CHAMP, a dedicated Earth-gravity detection satellite, in 2000, the United States and Germany have successively launched GRACE (2002) and GOCE (2009) Earth-gravity satellites. my country is also expected to launch its own gravity satellite this year. The arrival of a large number of satellite gravity data provides the possibility to study the global three-dimensional density distribution, which will revolutionize the study of the Earth's lithosphere, hydrosphere and atmosphere and their interactions in the near future. However, how to use large-scale satellite gravity data inversion to obtain a reliable global three-dimensional density distribution is facing challenges such as large computing scale and large amount of observation data. First, in order to invert the Earth's arbitrarily complex density distribution, we need to divide it into discrete units. In the spatial domain, we can use different types of geometric discretization elements, including: polyhedrons, prisms, etc. However, due to the influence of the earth's curvature, the gravitational field error caused by the discrete method of polyhedral and prismatic elements can reach several mGal, and the Zhang Li error of the gravitational gradient can reach 0.01E, which has greatly exceeded the observation error of satellite gravity. In addition, satellite gravity data is globally covered and observations will be in the millions. More than 233 million discrete spherical crust units (based on the ETOP1 Topography-Bathymetry model released by NOAA) are required to accurately describe the upper crust and topography. If the mass information of the lower crust, mantle, and core materials and their The uneven changes in the horizontal longitude and latitude direction, the number of units should be in the tens of billions. For N satellite gravity observation data (millions) and M discrete mass units of the earth (tens of billions), the time complexity of one forward modeling is O(NM), and one forward modeling using the current mainstream computing workstation requires at least 600 For many days, the inversion usually takes more than a few decades. Forward modeling is the engine of inversion, and the accuracy and efficiency of forward modeling largely determine the accuracy and efficiency of inversion.

The current mainstream gravity forward acceleration methods, such as wavelet compression (Wavelet), fast Fourier transform (FFT) and fast multipole expansion (FMM), are all aimed at local gravity exploration problems based on polyhedral or prismatic elements. For example, as disclosed in CN201611252162.X, the grid is discretized based on tetrahedron. For spheres such as the earth, the calculation error is very large. Therefore, for the global gravity problem based on satellite gravity data, we must use spherical shell units to better simulate the earth. surface. But Wavelet and FFT algorithms cannot deal with the spherical shell element problem due to the limitations of the algorithms themselves, while the FMM algorithm has better expansibility and applicability than other accelerated algorithms. But at present, there is no relevant report on the application of the FMM acceleration algorithm in the satellite gravity forward modeling problem based on spherical shell elements.

Therefore, the development of a fast and high-precision satellite gravity forward modeling engine based on spherical shell elements and FMM algorithm will lay the foundation for the acquisition of a global reliable 3D density model.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide a gravity forward modeling acceleration method in a spherical coordinate system, which can greatly improve the satellite gravity forward modeling by dividing spherical shell elements and performing key forward modeling calculations on any distributed satellite gravity observation point. computational efficiency. And the method is also applicable to other spherical planets, not only limited to the earth.

A gravity forward acceleration method in a spherical coordinate system, comprising the following steps:

Step 1: Region discretization and octree establishment:

In the spherical coordinate system, if the density distribution of the earth Ω is ρ(s), the gravity potential φ(o), the gravity gradient g(o) and the gravity gradient tensor T(o) generated at the observation point o can be expressed as :

Among them, ▽ is the gradient operator, R=|os| is the distance from the observation point o to the source point s, ρ(s) is the density of the earth Ω at the source point s, G is the gravitational constant, G=6.674×10 ^{− 11} m ³ kg ^-1 s ^-2 ;

Model establishment: Use the latest continental-ocean terrain relief model, crust, mantle and core layer to establish a global three-dimensional geometric model, as shown in Figure 1, and design model parameters (model coordinates and density distribution) according to the global density distribution characteristics; Point design observation point parameters (coordinates and number of observation points). For example, in practical applications, the observation point is a satellite.

Regional discretization: use spherical shell elements to discretize the global 3D model, and discretize the earth into multiple spherical shell elements Among them, T _i represents the i-th spherical shell unit, and N represents the number of spherical shell units contained in the earth. The left figure in Figure 4 is a schematic diagram of a spherical shell element, which consists of three surface pairs in a spherical coordinate system: one surface pair has constant radial distances r ₁ and r ₂ , and one surface pair has constant polar Angles (or constant dimensions in geographic coordinate systems) θ ₁ and θ ₂ , a surface pair with constant azimuth (or longitude in geographic coordinate systems) and A spherical shell element is represented by its center point.

Octree establishment: Octree is a tree-like data structure used to describe three-dimensional space. The present invention adopts the data construction method of octree to establish the octree structure of observation points and global three-dimensional discrete units, namely observation points Octree and spherical shell unit (source) octree, Figure 2 is a quadtree structure division similar to octree division;

Establish a corresponding octree structure for the observation points: build a cube unit that just contains all the observation points, called the highest layer cell, that is, the 0th layer cell, and then take the 0th layer cell as the mother cell, and divide it into 8 Each daughter cell is recorded as the first layer of cells; continue to divide each daughter cell in the above manner, and the cells obtained from each subsequent layer are recorded as the i-th (i=0,1,2,...) layer of cells in turn, Until the number of observation points contained in each cell is less than the set number or the number of observation point octree layers i reaches the set number of layers, and delete all cells that do not contain observation points; cells without child cells are called leaf cells;

Build an octree on the source:

Establish a corresponding octree structure for the spherical shell element of the global three-dimensional model Ω: construct a cubic element that just contains the spherical shell element of the global three-dimensional model Ω, which is called the highest layer cell, that is, the 0th layer cell, and then start with the 0th layer. The layer of cells is the mother cell, and it is divided into 8 daughter cells, which are recorded as the first layer of cells; continue to divide each daughter cell in the above manner, and the daughter cells obtained by each subsequent layer of division are recorded as the jth (j= 0, 1, 2, ...) layer cells, until the number of spherical shell elements contained in each cell is less than the set number or the source octree layer j reaches the set number of layers, and delete all cells that do not contain Cells of spherical shell cells; cells without daughter cells are called leaf cells;

Step 2. The division of the near area and the far area:

The conventional satellite gravity acceleration method is to integrate each spherical shell element once at each observation point, and then superimpose and sum to obtain the satellite gravity response of the observation point. The purpose of the adaptive fast multipole expansion algorithm is to divide the integral area into the far area and the near area for each observation cell. Referring to the division of the near area and the far area disclosed in CN201611252162.X, for the near area integration of each observation point, use The conventional satellite gravitational acceleration methods are the same, and since the realization process of the near-area integration is a means of the prior art, the present invention will not describe it too much. For the integration of the far area of the observation point, by selecting the far area expansion point, the integral of each observation point and all the far area spherical shell elements is converted into the connection between the observation point and the far area expansion center point, thereby reducing the number of integral calculations and achieving purpose of acceleration. In order to achieve the purpose of this acceleration, the connection between the observation point and the spherical shell unit is divided into three connection relationships: the first is the connection between the center point of the child cell of the spherical shell unit and the center point of the mother cell, which is connected by M2M conversion; The second is the connection between the center point of the parent cell of the spherical shell unit and the center point of the parent cell of the observation point, which is connected by M2L conversion; the third is the connection between the center point of the mother cell of the observation point and the center point of the daughter cell, which is connected by L2L conversion. As shown in Figure 3.

For any pair of observation point cells and source cells By calculating the relative position metric factor λ, and comparing the relative position metric factor λ with the preset distance value λ ₀ , the result of dividing the near area and the far area is obtained. Among them, the side lengths of the cell cube at the observation point and the source cell cube are recorded as L _o and L _s respectively, the distance between the center point of the observation point cell and the center point of the source cell is recorded as L _os , and the relative position measurement factor λ measures :

If λ>λ ₀ , the source cell is considered cell at the observation point The far region of , denoted as T _far ; otherwise, the source cell Cells at the observation point The near area of λ is denoted as T _near , wherein, λ ₀ is a preset value, and preferably its value range is: λ ₀ =[5,10]. For each observed cell there are:

Ω＝T _near +T _far (3)

Among them, the near area and the far area of the leaf cell of the observation point are the near area and the far area of the observation point in the observation point cell.

The present invention adopts the above ratio method to divide the near area and the far area. Compared with the distance judgment of the near and far area division of CN201611252162.X, the near area and the far area of the large observation cell can be more accurately divided, so that the final area can be improved. The forward calculation accuracy of .

Step 3. Calculate the gravity parameters of the near area and the far area respectively:

The integral of the satellite gravity forward modeling is divided into the integral of the near area and the integral of the far area; the integral of the near area is directly calculated based on the Gaussian numerical integration of the spherical shell elements contained in the near area, and the gravity parameters of the near area are obtained; and based on the observation For the point octree and the source octree, an adaptive fast multi-pole expansion algorithm is used to calculate the integral of the far area, and the gravity parameter of the far area is obtained; and then the satellite gravity forward modeling result is obtained by summing.

The near-area gravity parameters include the near-area gravity potential, near-area gravity gradient and near-area gravity gradient tensor generated at the observation point o, respectively denoted as and The far-area gravity parameters include the far-area gravity potential, the far-area gravity gradient, and the far-area gravity gradient tensor generated at the observation point o, respectively denoted as and

Step 4: Summing up the near-area gravity parameters and the far-area gravity parameters of the same type of gravity parameters of the same observation point calculated in step 3 respectively to obtain various types of gravity parameters of each observation point.

Among them, the near-area gravity parameters calculated in step 3 and the far-area gravity parameters are summed to obtain the gravity parameters of the observation point, including the gravity potential φ(o), the gravity gradient g(o) and the gravity generated at the observation point o The gradient tensor T(o), its calculation formula is:

The following first describes in detail the principles and steps of using the adaptive fast multi-pole expansion algorithm to calculate the integral of the far area to obtain the gravity parameters of the far area:

In order to use the fast multipole expansion algorithm, the integration point needs to be transformed, that is, the integration kernel function 1/|o-s| is expanded. When it is located in the far region, the integration kernel function can be expanded as:

in, and are the spherical coordinates of the observation point o and the integration point s, respectively, s _c is the center point of the far region Ω _far , and P is the integral kernel function The multipole expansion order of , the angle Θ satisfies R _nm and S _nm are conventional and singular spherical harmonics, respectively, which depend on n-order and m-order associated Legendre functions ^[1] , is the complex conjugate of S _nm .

Substituting Equation (5) into Equation (1) can obtain the expansion of the gravitational potential generated by the far-area spherical shell element:

where ρ(s) varies with the position in the far region Ω _far , is the multipole distance of the gravitational potential at the extended center point s _c of Ω _far in the far region. From equation (6), it can be seen that the multipole distance The integral kernel function has nothing to do with the observation point, which means that when multiple observation points correspond to the same far area, I only need to calculate the far area integral and some simple algebraic operations, which can greatly improve the gravity forward calculation. efficiency.

The principle of the acceleration process of the adaptive fast multipole expansion algorithm is described in detail below by taking the simple observation point octree structure and the spherical shell element octree structure in FIG. 3 as an example.

Source octree multipole distance calculation:

Consider first the bottom cell (leaf cell) of the spherical shell cell tree Assume for leaf cells the center of the source point Indicates the integration point, the points in the entire spherical shell element are the source points, therefore, the far area leaf cells The multipole distance can be expressed as:

Among them, we assume that leaf cells There are t spherical shell elements, namely

In addition, r, Represents the radial distance, azimuth, or latitude, inclination, or longitude of a geographic coordinate system in a spherical coordinate system, respectively. The spherical shell element can be transformed into a unit cube element as shown in the right figure of Figure 4 by the following coordinate transformation formula, in which the local variables η and γ satisfy -1≤ η,γ≤1,

Substitute (8) into (7) to get:

Since the above integral is defined on the unit cube, and the integral kernel function term R _nm is is also non-singular, so it can be solved using the standard Gaussian numerical integral formula ^[5] .

Source Octree Multipole Distance Upward Shift (M2M): Calculate the other leaf cells in Figure 3 in advance using formula (9) multipole distance. The multipole distance of its mother cell is obtained by transferring the multipole distance of the daughter cells, which can be calculated by the following recursive formula ^[2] :

Among them, j=7 and l=9, 10, 11, 12 represent the source octree layer 4 multipole transition to the third layer multipole distance, j=4 and l=6, 7, 8 represent the third layer The multipole pitch shifts to layer 2. Therefore, Equation (10) is also known as the transfer of multipole moments to multipole moments in the upper layers of the source-observation tree (Multipole moments to multipole moments, M2M, see Fig. 3). The multipoles of all cells on the source octree can be obtained by transferring the multipoles of leaf cells to the multipoles of root cells.

Observation point octree local coefficient calculation:

Now our task is to transfer the multipole distance calculated in the source octree to each cell of the observation octree. In order to achieve this conversion with observation points, we first need to introduce a The local expansion point _oc of , _oc is the observation cell the center point. Therefore, similar to formula (5), the integral kernel function 1/R can be expanded at the local expansion point _oc :

in, _s∈Ωfar .

Similarly, by substituting Eq. (11) into Eq. (1), another form of the gravitational potential produced by the remote geological body can be obtained:

Through the above formula, we can calculate the gravitational potential generated by the far-area geological body at the observation point o, where is defined as the local coefficient (ie, the integral contribution of the far-zone geological body at the extended center point _oc ).

Transfer of source octree multipole distance to observation point octree local coefficients (M2L): by means of the calculated multipole distance on the source octree The local coefficient in equation (12) can be quickly calculated by the following conversion formula ^[2]

in, For the multipole of the far-area geological body Ω _far at its extended center point sc, equation (13) is called the transfer of multipole moments to local coefficients ( _M2L , see Fig. 3).

Observation point octree local coefficient shift down (L2L): if the observation cell contains mother cells Then some of its local coefficients need to be calculated from the local coefficients of its parent cell ^[2] :

in, is the observation point mother cell The above formula is called the transfer of local coefficients from the upper layer of the observation tree to the local coefficients of the lower layer (local coefficients to local coefficients, L2L, see Figure 3).

Based on the above-mentioned principle analysis, it can be known that the near and far regions of the observation point in the present invention are the near and far regions of the leaf cells of the observation point, and the gravitational position at the observation point corresponds to the gravitational position of the leaf cell of the observation point. Therefore, if in the above formula 11 and formula 12, _oc is the observed leaf cell The local expansion point of , that is, the center point, and o is the observation point, then formula 12 represents the far-area gravity position of the leaf cell at the observation point, that is, the far-area gravity position of the observation point. It can be seen from formula 12 that to obtain the far-area gravity position of the observation point, the local coefficient of the leaf cell of the observation point needs to be calculated. In the same way, according to formula 1 and formula 4, the far area gravity parameter of the observation point to be obtained by the present invention can be expressed as follows ^[4] , it can be seen from the following formula, whether it is to calculate the far area gravity position Far Gravity Gradient or the far-region gravity gradient tensor All need to first calculate the local coefficients of the leaf cells at the observation point.

Based on the above theoretical analysis, the present invention uses the above theory to obtain the local coefficients at each observation point o The process is as follows:

A: Calculate the multipole moment of each source cell in the source octree in order from the bottommost cell of the source octree to the topmost cell;

First, the multipole moments of leaf cells in the source octree are calculated, and then the multipole moments of other source cells are calculated according to the M2M transfer formula. Sum of transfer values of M2M of polar moments.

As shown in the source octree in Figure 3, first use Equation 9 to calculate The multipole moments of the leaf cells in the four source octrees are then calculated according to Equation 10 The transfer value of M2M of leaf cells in each source octree in , and then summed to get the multipole moment. By analogy, the multipole moment of each source cell in the source octree is calculated. It should be noted that when the far regions are the same, even for different observation points, their multipole moments are the same. The present invention first calculates the multipole moment of each source cell in the source octree, which can improve the efficiency of the whole algorithm, because this scheme only needs to obtain the multipole moment of the corresponding far area when calculating the local coefficient of any observation point. It is enough to perform the conversion, and there is no need to judge whether the multipole moment in the far region already exists; therefore, compared with the conventional method of obtaining the local coefficient of an observation point, a path identification is performed, and then the matching source cells are judged in turn according to the path identification result. Whether the multipole moment exists or not is calculated, the method of the present invention can improve the operation efficiency and reduce the logical judgment.

B: Calculate the local coefficients of each observation point in the order of the observation point octree from the top cell to the bottom cell;

First, based on the multipole moment of the source cell on the source octree and according to the M2L transfer formula, the local coefficient of the top cell on the observation point octree is calculated, and then based on the M2L and L2L transfer formulas, the branch where the leaf cell of the observation point is located is calculated in turn. The local coefficients of each cell on the road are obtained until the local coefficients of the leaf cells at the observation point are obtained.

Among them, if the observation point octree has n layers, and the source octree has N layers, let the topmost cell be the first layer of the octree, and the observation point cell in the observation point octree The process of obtaining the local coefficients is briefly described as follows:

a: If the point cell is observed At level 1 in the observation point octree:

a1: Determine whether the first layer of cells in the source octree is an observation point cell the distant origin cells;

If so, calculate the direction from the first layer cell in the source octree to the observation point cell The transfer value of M2L is obtained from observation point cells The local coefficient of ;

If not, obtain the distant source cells in the second layer of cells, and then calculate the direction of each distant source cell in the second layer of cells in the source octree to the observation point cell. The sum of the transfer values of M2L obtains the local coefficient of the observation point cell;

That is, for the topmost cell on the observation point octree, obtain the observation point cell from the top two layers (1st and 2nd layer cells) in the source cell octree The distant region of origin cells, re-identified. It can be seen from this that in this implementation, it is not necessary to identify whether each source cell is a near area or a far area for an observation point cell, but only the number that matches the algorithm needs to be identified. For example, for the first layer of observation point cells , you need to identify the source cells in layers 1 and 2 in the source octree.

b: If the point cell is observed not leaf cells, observation point cells and observation point cells mother cell They are located at the i-th layer and the i-1-th layer in the observation point octree respectively, where 1 in 2 is at the observation point and i and 1: such source cells in the source octree have not all been detected and the transfer has been completed.

b1: Determine whether the source cell is the observation point cell in the i-th layer of cells in the source cell octree cells of distant origin and not mother cells The far region source cell of , if it exists, obtain the far region source cell and execute b2; if it does not exist, execute b2;

b2: Determine whether there are observation point cells in the i-th layer cells in the source octree of proximal cells;

If it exists and the near-region source cell is not a leaf cell, determine whether there is an observation point cell in the daughter cell of the near-region source cell The far region source cell, if there is, obtain the far region source cell, and perform b3; if there is a near region source cell and the near region source cell is a leaf cell; In the presence of distant source cells, execute b3;

b3: Calculate the direction of all distant source cells to the observation point cells The sum of the M2L transfer values and the calculated blast to the observation point cell The L2L transfer value is calculated, and then the sum of the M2L transfer value and the L2L transfer value is calculated to obtain the observation point cell The local coefficients of ; 2 local coefficients; to;

In the source octree, if the parent cell is in the far region, the source cells under any branch below the parent cell are in the far region; if the parent cell is in the near region, there must be a source cell in the daughter cell that is in the near region. Therefore, the present invention first identifies whether an active cell in the near region of the parent cell of the observation point becomes the far region of the daughter cell of the current observation point. The multipole moment is not transferred to the parent cell of the observation point; if the near-region source cell of the parent cell of the observation point is still the near-region source cell of the daughter cell of the observation point, but the daughter cell of the near-region source cell may be the parent cell of the observation point. The far-region source cells of the point cell and the multipole moment of this layer of cells are not included in the local coefficient of the parent cell of the observation point, that is, part of the local coefficient of the observed point daughter cell is transferred from the distant source cell M2L, so it is necessary to The multipole moments of the distant-derived cells are transferred to the observation point daughter cells through M2L.

c: If the point cell is observed and observation point cells mother cell They are located at the i-th layer, i-1-th layer, and 2 in the observation point octree, respectively, and the observation point cell is at the time of view and i and at view time The local coefficients of are derived only from the mother cell L2L transfer of the local coefficient, the mother cell was calculated using the L2L transfer formula to the observation point cell The L2L transfer value is obtained from the observation point cell local coefficients. This class is for all source cells in the source octree whose multipole moments have been transferred, and the local coefficient of the observation point cell is only derived from the L2L transfer of its parent cell.

d: If the point cell is observed The leaf cell located in the observation point octree, the initial value of i is the observation point cell Layer number, observation point cell The mother cell is expressed as And n and fine time:

d1: Determine whether the source cell is the observation point cell in the i-th layer of cells in the source cell octree cells of distant origin and not mother cells If there is a distant source cell, obtain the distant source cell and execute d2; if it does not exist, execute d2;

d2: Determine whether there is an observation point cell in the i-th layer of cells in the source octree of proximal cells;

If it exists and the source cell in the near region is not a leaf cell, then set i=i+1, and return to step d1; wherein, the value of i after each cycle cannot exceed N;

If it exists and the source cell in the near region is a leaf cell, perform step d3;

If there is no near-region source cell, perform step d3;

d3: Calculate the sum of the M2L transfer values of all distant source cells obtained and calculate the mother cell to the observation point cell The L2L transfer value is calculated, and then the sum of the M2L transfer value and the L2L transfer value is calculated to obtain the observation point cell local coefficients.

The present invention identifies from which source cells the cells at each observation point have also undergone M2L transfer through the above logic, which is more logical and can prevent the repeated transfer of multipole moments of the source cells.

Through the above local coefficient calculation process, the local coefficients of all cells on the octree at the observation point can be calculated from top to bottom, and finally the local coefficient of the leaf cell at the observation point can be obtained. Far Region Gravity Gradient and Far Region Gravity Gradient Tensor. Finally, the gravity response of each observation point can be obtained by substituting the near-region gravity response and the calculated far-region gravity response calculated by Gaussian numerical integration into Eq. (4).

beneficial effect

1. The traditional satellite gravity forward acceleration method needs to integrate all spherical shell elements for each observation point, and then superimpose and sum to obtain the gravity value of the observation point; if the number of observation points is M, the spherical shell element is If the number is N, the time complexity of the satellite gravity forward modeling of the traditional method is O(MN). The invention adopts an adaptive multi-layer fast multi-pole expansion algorithm, divides the integral area into a near area and a far area, and adopts the traditional gravity acceleration method for the integral of the spherical shell element in the near area; The integration of the observation point and all the spherical shell elements in the far region is converted into the connection between the observation point and the extended center point in the far region, and the integration between the extended center point in the far region and the spherical shell element in the far region, so as to realize the integration between the observation point and the spherical shell element in the far region. Decoupling of units. If the integral of the far-region extension center point and the far-region spherical shell element is calculated, for the far-region integration of the new observation point, only some simple algebraic operations are required to reduce the number of integration calculations, and the time complexity will be O(MN ) is reduced to O(MlogN), which can greatly improve its computational efficiency.

2. The present invention improves on the basis of the prior art CN201611252162.X, uses spherical shell elements for discrete, and constructs a corresponding coordinate system for calculation, and further improves the calculation accuracy for spherical applications, greatly reducing the The error of applying CN201611252162.X to a sphere is reduced, and the FMM acceleration algorithm is applied to the gravity forward modeling problem based on spherical shell elements for the first time.

3. The present invention also provides a process for deriving local coefficients of observation points, which first calculates the multipole moment of each source cell in the entire source octree, and then selects the multipole moment of the corresponding source cell for any observation point cell. At the same time, since the far and near regions of the observed point cell are different from those of the parent cell, the changes are irregular for different source cells. In the calculation of multipole moment transfer, it is very easy to repeat the transfer of multipole moments, and the transfer analysis of the local coefficient of the observation point from the source cell provided by the present invention is more logical and more rigorous, which can prevent the occurrence of repeated multipole moments. transfer problem.

Description of drawings

Figure 1 is a schematic diagram of a global three-dimensional density model.

FIG. 2 is a schematic diagram of the division of the quadtree structure, in which the black dots represent sector-shaped ring units.

Figure 3 is a schematic diagram of the Multipole moments to Multipole moments (M2M) transformation, Multipole moments to local coefficients (M2L) transformation, and Localcoefficients to local coefficients (L2L) transformation in the observation point octree and the source octree.

FIG. 4 is a schematic diagram of the mapping of the spherical shell element in the spherical coordinate system to the unit cube in the local Cartesian coordinate system.

Figure 5 is a schematic diagram of the 1/8 Preliminary Reference Earth Model (PREM) spherical shell element division. The number of elements in the entire PREM model is 11,145,600.

Detailed ways

The present invention will be further described below with reference to the embodiments.

S1: Design of measuring points and source files: determine the number and coordinates of measuring points, determine the density of the global 3D model and the corresponding position coordinate file; divide the 3D geometric model into spherical shell elements; The shell element establishes a corresponding octree structure, and for each observation cell, the integration area is divided into its own near area and far area.

S2: Calculation of gravity parameters in the far area: Calculate the far area integral through an adaptive fast multi-pole expansion acceleration algorithm.

S3: Calculation of gravity parameters in the near area: Calculate the near area integral by Gaussian numerical integration based on spherical shell elements.

S4: Calculate the gravity parameter of the observation point: The gravity parameter of the observation point is obtained by summing the near-area gravity parameter and the far-area gravity parameter of the same type of gravity parameter at the same observation point.

The following is an example of calculating the Preliminary Reference Earth Model (PREM) of the present invention.

As shown in Figure 5, the Preliminary Reference Earth Model (PREM) consists of 11 concentric sphere layers, namely the inner core, outer core, D″ layer, upper mantle, two inner transition zones, one outer transition zone, seismic low velocity zone, Inner crust, outer crust and ocean layer. The PREM model is discretized into 11,145,600 spherical shell elements along the three directions of radius, longitude and latitude. The specific grid and density distribution parameters are shown in Table 1, which is the grid of the discrete PREM model. and the density distribution parameter table, R ₁ , R ₂ represent the outer radius and inner radius, respectively, represent the number of divisions in longitude and latitude, respectively. The observation points are evenly distributed on the sphere above the earth (longitude from -180° to 180°, latitude from -90° to 90°), the height of the sphere and the surface of the PREM model is 250km, and the total number of observation points is 259,200.

Table 1

Due to the spherical symmetry of the PREM model, the gravitational responses (gravity potential, gravitational field and gravitational gradient tensor) generated by the PREM model are independent of the longitude and latitude of the observation point, only related to its radius. Therefore, the gravitational response observed on the observation plane of h=250km is constant, and is calculated by the analytical formula: gravitational potential Gravity Field Gravity Gradient Tensor The rest of the gravity gradient tensor components are zero.

Next, we use the adaptive multi-layer fast multipole algorithm (AMFMM) proposed by the present invention to calculate the gravity anomaly generated by the PREM model on the observation plane of h=250km, and compare the results with the analytical solution and the Gaussian numerical integral solution, Verify the accuracy and speed of the algorithm. Set the calculation parameters λ = 8, P = 8 and q = 8, and use the AMFMM algorithm to calculate the statistical results of the gravitational response of all the observation points on the observation surface as shown in Table 2, which is the gravitational response generated by the PREM model on the observation plane. AMFMM calculation result statistics table. Comparing the statistical results in Table 2 with the analytical solutions, it can be seen that the absolute maximum relative errors of the gravitational potential, gravitational field and gravitational gradient tensors are 1.6×10 ^-6 %, 2.2×10 ^-6 % and 3.2×10 ^{- 5} %. The maximum absolute error of the gravitational field gr is _0.02mGal , which is far less than the 1mGal accuracy of the instrument designed by the GOCE satellite [3]. The absolute value of the maximum relative error is 0.001E, which is less than 0.01E of the GOCE satellite instrument noise (Cesare et al., 2010)0. The results calculated by the AMFMM algorithm above are highly consistent with the analytical solution, which verifies the accuracy of the AMFMM algorithm. At the same time, for the satellite gravity forward modeling model with 259,200 observation points and 11,145,600 spherical shell elements, the calculation time of the AMFMM algorithm is 2056.4s. If the traditional gravity forward modeling method is used, the calculation time will be calculated in months.

Table 2 The statistical table of the AMFMM calculation results of the gravitational response generated by the PREM model on the observation plane.

It can be seen from the above calculation examples that the fast acceleration method for satellite gravity forward inversion proposed by the present invention will make it possible to obtain a global reliable three-dimensional density inversion model. It should be emphasized that the examples described in the present invention are illustrative rather than restrictive, so the present invention is not limited to the examples described in the specific implementation manner, and all the examples obtained by those skilled in the art according to the technical solutions of the present invention Other embodiments that do not depart from the spirit and scope of the present invention, whether modified or replaced, also belong to the protection scope of the present invention.

references:

[1] Greengard, L., and V. Rokhlin (1987), A fast algorithm for particleimulations, Journal of Computational Physics, 73(2), 325-348.

[2] Ken-ichi, Y. (2001), Applications of fast multipole method to boundary integral equationmethod, Ph.D.thesis, Dept.of Global EnvironmentEngrg.Kyoto Univ.,Japan March.

[3] Cesare, S., M. Aguirre, A. Allasio, B. Leone, L. Massotti, D. Muzi, and P. Silvestrin (2010), The measurement of earth gravity field after the gocemission, ActaAstronautica, 67(7 ),702-712,

[4] Ren, Z., T. Kalscheuer, S. Greenhalgh, and H. Maurer (2013), Boundaryelement solutions for broad-band 3-D geo-electromagnetic problems accelerated by an adaptive multilevelfast multipole method, Geophysical Journal International, 192(2) , 473-499.

[5] Jin, J.-M. (2014), The finite element method in electromagnetics, John Wiley & Sons.

Claims (8)

1. A gravity forward modeling acceleration method under a spherical coordinate system is characterized in that: the method comprises the following steps:

step 1: establishing a region discrete and octree:

firstly, determining a geometric model and density distribution of a global three-dimensional model omega, and designing the number and coordinates of observation points; wherein, the point in the global three-dimensional model Ω is called a source point;

dividing a global three-dimensional model omega into a plurality of spherical shell unitsWherein N is the number of spherical shell units, and the spherical shell units are represented by the central points thereof;

wherein the spherical shell unit consists of 3 surface pairs in a spherical coordinate system, one surface having a constant radial distance r₁And r₂A surface pair having a constant polar angle or a constant dimension theta in a geographic coordinate system₁And theta₂A surface pair having a constant azimuth or longitude in a geographic coordinate systemAnd

then, establishing corresponding octree structures for the observation points and the spherical shell units of the global three-dimensional model omega to obtain an observation point octree and a source octree;

establishing an octree structure for the observation points: constructing a cubic unit which just contains all observation points, namely a highest layer cell, namely a layer 0 cell, and dividing the layer 0 cell into 8 daughter cells which are marked as a layer 1 cell by taking the layer 0 cell as a mother cell; continuously dividing each daughter cell in the above manner, sequentially recording the cells obtained in each divided layer as the cells in the ith (i is 0,1,2, …) layer until the number of observation points contained in each cell is less than the set number or the number i of observation point octree layers reaches the set number, and deleting all the cells not containing the observation points; cells without daughter cells are called leaf cells;

establishing an octree structure for a spherical shell unit of a global three-dimensional model omega: constructing a cubic unit just containing a global three-dimensional model omega spherical shell unit, namely a highest layer cell, namely a layer 0 cell, dividing the layer 0 cell into 8 daughter cells by taking the layer 0 cell as a mother cell, and marking as the layer 1 cell; continuously subdividing each sub-cell in the manner, sequentially recording the sub-cells obtained in each subdivision layer as the j (j is 0,1,2, …) th layer cells until the number of spherical shell units contained in each cell is less than the set number or the number j of source octree layers reaches the set number, and deleting all cells not containing the spherical shell units; cells without daughter cells are called leaf cells;

step 2: dividing a global three-dimensional model omega into a near zone and a far zone:

confirming a near zone and a far zone of each observation point cell, wherein the near zone and the far zone of the leaf cell of the observation point are the near zone and the far zone of the observation point;

any pair of observation point cells ○_oAnd source cell ○_sCalculating a relative position metric factor lambda according to the following formula, and comparing the relative position metric factor lambda with a preset distance value lambda₀Comparing to obtain a near zone and a far zone division result;

wherein λ > λ₀Indicates the source cell ○_sCell ○ at observation point_oA distal region of (a); otherwise, the position is in the near zone;

wherein L is_oAnd L_sRespectively observation point cells ○_oThe cube and the source cell ○_sLength of the cube, L_osCell ○ as observation point_oCentral point and source cell ○_sThe distance between the center points;

and step 3: respectively calculating near-zone and far-zone gravity parameters of an observation point, wherein the gravity parameters comprise a gravity potential, a gravity gradient and a gravity gradient tensor;

calculating the integral of the near zone by adopting a Gaussian numerical integration formula of a spherical shell unit contained in the near zone to obtain a near zone gravity parameter; calculating the integral of the remote area by adopting a self-adaptive rapid multipole expansion algorithm based on the observation point octree and the source octree to obtain a remote area gravity parameter;

the gravitational potential, gravitational gradient and gravitational gradient tensors are respectively expressed as follows:

wherein phi (o), g (o), and gamma (o) respectively represent the gravitational potential, gravitational gradient, and gravitational gradient tensor at observation point o,the gradient operator is R ═ o-s |, which is the distance from the observation point o to the source point s, ρ(s) is the density of the global three-dimensional model Ω at the source point s, and G is the gravity constant;

and 4, step 4: and (3) summing the near-region gravity parameters and the far-region gravity parameters of the same observation point calculated in the step (3) respectively to obtain various gravity parameters of each observation point.

2. The method of claim 1, wherein: let o_cTaking the central point of leaf cell ○ o at the observation point, wherein o is the observation point, and the calculation formula for calculating the integral of the far zone by adopting the adaptive fast multipole expansion algorithm based on the octree at the observation point and the source octree in step 3 to obtain the gravity parameter of the far zone is as follows:

wherein,respectively representing the remote gravity potential, the remote gravity gradient and the remote gravity gradient tensor at the observation point o, wherein P is an integral kernel functionMultipole expansion order of R_nmIs a conventional spherical harmonic function, which depends on an n-th and m-th order banded legendre function,is a local coefficient and represents the geologic body central point o of a far zone_cThe integral contribution of (a).

3. The method of claim 2, wherein: obtaining local coefficients at each observation point oThe process is as follows:

a: sequentially calculating the multipole moment of each source cell in the source octree from the bottom layer cell to the top layer cell of the source octree;

firstly, calculating the multipole moment of leaf cells in a source octree, and then respectively calculating the multipole moments of other source cells according to an M2M transfer formula; wherein M2M represents the shift of the daughter cell multipole distance to the mother cell multipole distance in a source octree, the multipole moment of the source cell of all non-leaf cells being equal to the sum of the shift values of M2M of its daughter cell multipole moment;

b: calculating local coefficients of all observation points according to the sequence of the observation point octree from the topmost cell to the bottommost cell;

firstly, calculating local coefficients of topmost cells on an observation point octree based on a multipole moment of source cells on the source octree and according to an M2L transfer formula, and then sequentially calculating local coefficients of all cells on branches where leaf cells of the observation point are located based on M2L and L2L transfer formulas until local coefficients of the leaf cells of the observation point are obtained, wherein the local coefficients of the leaf cells of the observation point are the local coefficients of the observation point;

where M2L represents the transition of the source octree cell multipole distance to the octree cell local coefficient at the observation point, and L2L represents the transition of the mother cell local coefficient to the daughter cell local coefficient in the octree at the observation point.

4. The method of claim 3, wherein if the observation point octree has N layers, the source octree has N layers, the topmost cell is layer 1 of the octree, and the observation point cell ○ in the observation point octree_oThe local coefficient obtaining process is as follows:

a cells ○ if observation point_oLayer 1 in the observation point octree:

a1 determination of whether layer 1 cells in source octree are observation point cells ○_oThe remote source cell of (3);

if yes, calculate layer 1 cell in source octree to observation point cell ○_oM2L metastasis results in observation point cell ○_oThe local coefficient of (a);

if not, obtaining far-region source cells in the layer 2 cells, and calculating the direction of each far-region source cell in the layer 2 cells in the source octree to the observation point cell ○_oThe sum of the transfer values of M2L gives the local coefficients of the cells at the observation point;

b-cell ○ if observed Point_oObservation site cells ○, not leaf cells_oAnd observation site cell ○_oMother cell of (2)Respectively located at the ith layer and the (i-1) th layer in the observation point octree, wherein in 2, the observation point octree, i and view are as follows:

b1 judging whether the source cell exists in the ith layer of cells in the source cell octree, wherein the source cell is an observation point cell ○_oIs not a mother cellIf present, obtaining said remote source cell and performing b 2; if not, execute b 2;

b2 judging whether observation point cells exist in the ith layer of cells in the source octree ○_oA cell of near zone origin;

if the cells exist and the near region source cells are not leaf cells, judging whether observation point cells ○ exist in the sub-cells of the near region source cells or not_oIf present, obtaining said remote source cell, and performing b 3; b3 if there is a near region source cell and the near region source cell is a leaf cell and if there is no near region source cell and there is no far region source cell in the daughter cell;

b3 calculating the direction of all far-region source cells to observation point cells ○_oThe sum of the transfer values of M2L and calculation of the mother cellCells ○ to observation point_oThen calculating the sum of the transfer value of M2L and the sum of the transfer value of L2L to obtain observation point cell ○_oThe local coefficient of (a); 2, local coefficient; to;

c-if observation point cell ○_oAnd observation site cell ○_oMother cell of (2)At levels i, i-1 and 2 of the observation point octree, and at view and i and at view, the observation point cell ○_oThe local coefficients of (A) are derived only from the mother cellsThe local coefficient L2L transfer is used to calculate the mother cell by L2L transfer formulaCells ○ to observation point_oTransfer of L2L to give observation point cell ○_oThe local coefficient of (a);

d if observation point cell ○_oLeaf cells in the observation point OctreeImage and the number of layers is represented by i, observation point cell ○_oRepresentation of mother cellIs composed ofAnd when n is small:

d1 judging whether the source cell is present in the i-th layer of cells in the source cell octree, and if so, determining that the source cell is the observation point cell ○_oIs not a mother cellIf the far-region source cell exists, obtaining the far-region source cell and executing d2, if the far-region source cell does not exist, executing d2, d2, judging whether observation point cells ○ exist in the ith layer of cells in the source octree_oA cell of near zone origin;

if the cells exist and the source cells in the near region are not leaf cells, the step d1 is returned to if i is equal to i + 1; wherein the value of i after each cycle cannot exceed N;

if present and the near zone source cell is a leaf cell, performing step d 3;

if no near-zone source cell exists, performing step d 3;

d3 calculating the direction of all far-region source cells to observation point cells ○_oThe sum of the transfer values of M2L and calculation of the mother cellCells ○ to observation point_oThen calculating the sum of the transfer value of M2L and the sum of the transfer value of L2L to obtain observation point cell ○_oThe local coefficient of (2).

5. The method of claim 3, wherein: in the source octree, letAre leaf cellsCenter, source point ofLeaf cells in the source octree were calculated according to the following formulaMulti polar moment of

Wherein t is derived from leaf cells in octreeNumber of middle spherical shell elements, i denotes the ith spherical shell element, ρ_iThe density at the ith spherical shell unit is shown as r,Respectively, a radial distance, an azimuth angle, or a latitude, an inclination angle, or a longitude of a geographic coordinate system in a spherical coordinate system.

6. The method of claim 3, wherein: in the source octree, letIs a daughter cellMother cell of (2)Center of (2), transfer value of one daughter cell's multipole moment to M2M of mother cell's multipole momentThe following were used:

in the formula,representing daughter cellsThe multipole moment of (a).

7. The method of claim 4, wherein: in the observation point octree and the source octree, let o_cCell ○ as observation point_oCentral point of(s)_cCell ○ as observation point_oDistal cell ○_sFrom the source cell ○_sMulti-polar moment to observation point cell ○_oTransition value of M2L for local coefficientThe following were used:

in the formula,indicating source cell ○_sThe multipole moment of (a).

8. The method of claim 4, wherein: in the observation point octree, a point is setCells as observation pointsCentral point of (1), point o_cCell ○ as observation point_oIs measured at a central point of the beam,is ○_oThe mother cell of (2), from the observation pointLocal coefficients to observation point cell ○_oL2L branch value of local coefficientThe following were used:

in the formula,cells as observation pointsMidpointLocal coefficient of (2)