Gravity rapid inversion method and system based on multivariate geometric equivalence

Abstract

The invention is suitable for the technical field of gravity rapid inversion, and provides a gravity rapid inversion method and system based on multivariate geometric equivalence, and the method comprises the following steps: dividing an underground region by using cuboids with the same length in x and y directions; calculating a kernel function matrix between the single hexahedron subdivision unit with the node equivalent response and the earth surface observation point by using a point equivalent method; equivalently expanding the response of the subdivision unit to an initial layer through an application body equivalent relationship to obtain a kernel function matrix of all the subdivision units of the layer to earth surface observation points, equivalently transferring the response of nodes in the layer to a next layer, and performing layer-by-layer calculation to obtain a kernel function matrix of the whole region; and performing inversion by adopting a conjugate gradient method, and obtaining gravity physical property distribution of an inversion region by utilizing fast Fourier transform. According to the method, a large amount of repeated calculation is reduced, rapid calculation of the gravity kernel function is achieved, and compared with a traditional inversion scheme, inversion precision and resolution cannot be lost.

Description

Technical Field

[0001] The present invention belongs to the technical field of gravity rapid inversion, and particularly relates to a gravity rapid inversion method and system based on multi - element geometric equivalence. Background Art

[0002] The density inversion of gravity data usually discretizes the underground space into units, establishes the mapping relationship between each unit and the observation points through the forward formula, calculates the physical property distribution of each unit in the underground space, and delimits the distribution of mineral resources according to its physical properties by establishing a three - dimensional underground density structure. However, in density inversion, the calculation of the kernel function is a very time - consuming process, which will affect the exploration efficiency; therefore, developing a fast and accurate gravity density kernel function calculation method has become an urgent problem to be solved.

[0003] There are many methods to achieve fast forward calculation in gravity inversion, such as forward modeling methods in the Fourier domain, parallel computing, etc. However, the forward calculation in the Fourier domain will distort the sensitivity, thus reducing the inversion accuracy; parallel computing has certain hardware limitations, for example, it requires a graphics processing unit (GPU) to support the Compute Unified Device Architecture platform and specific hardware performance requirements.

[0004] The equivalent geometric lattice technology reduces memory usage and hardware requirements without reducing accuracy, can reduce the repeated calculation of the sensitivity matrix elements, and improves the calculation efficiency. Therefore, it has been widely used. In addition, when using an equal - amount uniform grid, the sensitivity matrix adopts a special block - Toeplitz matrix; however, this method only considers the equivalence of block - based equivalent units and ignores the equivalent properties existing between other structures, resulting in limited efficiency improvement. Summary of the Invention

[0005] The purpose of the embodiments of the present invention is to provide a gravity rapid inversion method based on multi - element geometric equivalence, aiming to solve the problems proposed in the above - mentioned background art.

[0006] The embodiments of the present invention are implemented as follows. The gravity rapid inversion method based on multi - element geometric equivalence includes the following steps:

[0007] Grid the gravity anomaly data, and divide the underground area with a cuboid with equal lengths in the x and y directions.

[0008] Apply the point - equivalence method to calculate the relationship function between the nodes of the underground divided grid and the surface observation points, and calculate the kernel function matrix between a single hexahedron divided unit with node - equivalent response and the surface observation points.

[0009] By applying the volume equivalent relationship, the response of the subdivision unit is equivalently extended to the initial layer to obtain the kernel function matrix of all subdivision units in this layer for the surface observation points, and the response of the nodes in this layer is equivalently transferred to the next layer, calculating layer by layer to obtain the kernel function matrix of the whole area;

[0010] The conjugate gradient method is used for inversion, and the fast Fourier transform is utilized to obtain the gravity physical property distribution of the inversion area.

[0011] Another object of the embodiments of the present invention is to provide a gravity rapid inversion system based on multi - element geometric equivalence for implementing the above - mentioned gravity rapid inversion method based on multi - element geometric equivalence, including:

[0012] A subdivision module for meshing gravity anomaly data and subdividing the underground area with cuboids having equal lengths in the x and y directions;

[0013] A point equivalence calculation module for calculating the relationship function between the nodes of the underground subdivision grid and the surface observation points by applying the point equivalence method, and calculating the kernel function matrix between a single hexahedron subdivision unit with node equivalent response and the surface observation points;

[0014] A volume equivalent technology module for equivalently extending the response of the subdivision unit to the initial layer by applying the volume equivalent relationship to obtain the kernel function matrix of all subdivision units in this layer for the surface observation points, and equivalently transferring the response of the nodes in this layer to the next layer, calculating layer by layer to obtain the kernel function matrix of the whole area;

[0015] An inversion module for performing inversion by using the conjugate gradient method and obtaining the gravity physical property distribution of the inversion area by utilizing the fast Fourier transform.

[0016] The gravity rapid inversion method based on multi - element geometric equivalence provided by the embodiments of the present invention reduces a large amount of repeated calculations, can achieve the rapid calculation of the gravity kernel function, and does not lose the inversion accuracy and resolution compared with the traditional inversion scheme. Description of the Drawings

[0017] Figure 1 It is a flow chart of the gravity rapid inversion method based on multi - element geometric equivalence provided by the embodiments of the present invention;

[0018] Figure 2 It is a schematic diagram of the overall subdivision of the grid provided by the embodiments of the present invention;

[0019] Figure 3 It is a schematic diagram of the point equivalence calculation provided by the embodiments of the present invention;

[0020] Figure 4 It is a schematic diagram of the volume equivalence calculation provided by the embodiments of the present invention;

[0021] Figure 5 Schematic diagram for transmitting node responses when calculating the next layer provided by an embodiment of the present invention. Detailed implementation manners

[0022] In order to make the objectives, technical solutions and advantages of the present invention clearer and more understandable, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0023] The existing gravity equivalent fast inversion method only considers the dissected blocks as the basic equivalent units, but ignores the equivalence of nodes when calculating the responses of the dissected blocks. Therefore, a large number of repeated calculations will occur in the calculation, thus increasing the amount of calculation. On the same horizontal plane, there is a certain equivalence among nodes with the same distance to the observation point. At the same time, when calculating the upper and lower layers, there are shared nodes, and there is also a certain equivalence. At the same time, on the same horizontal plane, for dissected blocks with the same distance to the observation point, the same equivalence is also satisfied. The embodiments of the present invention combine the equivalent relationship of nodes with the equivalent relationship of bodies, so as to replace repeated calculations with equivalent relationships, reduce the amount of calculation, and thus improve the calculation efficiency.

[0024] The following describes the specific implementation of the present invention in detail with specific embodiments.

[0025] As Figure 1 shown, it is a flowchart of a gravity fast inversion method based on multi - element geometric equivalence provided by an embodiment of the present invention, including the following steps:

[0026] S1. Grid the gravity anomaly data, and dissect the underground area with cuboids having equal lengths in the x and y directions:

[0027] Dissect the underground area with cuboids having equal lengths in the x and y directions. The distribution of the dissection grid is consistent with the terrain undulation, and the observation points correspond one - to - one with the dissection units. The dissection schematic diagram is as Figure 2 shown.

[0028] S2. Apply the point - equivalence method to calculate the relationship function between the nodes of the underground dissection grid and the surface observation points, and calculate the kernel function matrix between a single hexahedron dissection unit with node - equivalent response and the surface observation points:

[0029] Define the functional relationship between the nodes of a single dissection unit and the observation points as:

[0030]

[0031] Among them, O ij represents the observation point at position (i, j) in the observation - point grid, and its coordinates are (x i , y j, z0); Q r，p，q represents the node in the subdivision element located at the coordinate (ξ r , η p , ζ q ); R represents the subdivision node Q r，p，q to the observation point O i，j 's distance, and the value ranges of r, p, and q are all 0 or 1;

[0032] For four nodes on the same horizontal plane, as Figure 3 shown, there is an equivalent relationship, and the formula is expressed as:

[0033]

[0034] At the same time, there is also an identity relationship, and the formula is expressed as:

[0035] Cop(Q r，p，q (ξ r , η p , ζ q )); O i，j (x i , y j , z0)) = Cop(Q p，r，q (ξ p , η r , ζ q )); O j，i (x j , y i , z0))(3);

[0036] The observation point directly above the subdivision element and the four nodes on the same horizontal plane of the element have the following relationship:

[0037] Cop(Q r，p，q (ξ r , η p , ζ q )); O i，j (x i , y j , z0)) + Cop(Q r+1，p，q (ξ r+1 , η p , ζ q )); O i，j (x i , y j , z0)) + Cop(Q r，p+1，q (ξ r , η p+1 , ζ q )); O i，j (x i , y j , z0)) + Cop(Q r+1，p+1，q(ξ r+1 , η p+1 ,ζ q );O i，j (x i ,y j , z0))=0(4);

[0038] By applying equations (2)-(4) to calculate the functional relationship between nodes and observation points, more than 62.5% of the computational effort can be saved.

[0039] The calculated functional relationship is combined into the kernel function of the subdivision block and the observation point:

[0040]

[0041] S3. By applying the volume equivalence relation, the response of the subdivision unit is equivalently extended to the initial layer, and the kernel function matrix of all subdivision units in the layer to the surface observation point is obtained. The response of the nodes in the layer is equivalently transferred to the next layer, and the kernel function matrix of the entire area is obtained by calculating layer by layer.

[0042] After obtaining the kernel function of the subdivision unit to the observation point, the property of volume equivalence is applied, where the coordinates of the subdivision unit center point P are expressed as (α u , β v ,γ w ), the subdivision unit is at the observation point O i，j (x i ,y j ,z0) is M(P u，v (α u , β v , γ w );O i，j (x i ,y j , z0));

[0043] The translation equivalence of volume equivalence is as follows Figure 4 As shown, the formula is expressed as:

[0044] M(P 1，1 (α1, β1, γ w );Q i，j (x i ,y j , Z0))=M(P u，v (α u , β v , γ w );Q i+u-1，j+v-1 (x i+u-1 ,y j+v-1 , z0)) (6);

[0045] At the same time, there is cross equivalence, which is expressed as:

[0046] M(P u，v (α u ,β v ,γ w );Q i，j (x i ,y j ,z0)) = M(P i，j (x i ,y j ,γ w );Q u，v (α u ,β v ,z0)) (7);

[0047] Combining the above two properties, the calculated Mg response value can be written as a row vector and expanded in the form of a Toeplitz matrix, so as to obtain the response of all the subdivision units in this layer to the observation point;

[0048] When calculating the next layer, due to the existence of common nodes, the relationship functions of some nodes can be inherited to the next layer to reduce the amount of calculation, as Figure 5 shown:

[0049] Cop block-1 (Q r，p，q+1 (ξ r ,η p ,ζ q+1 );O i，j (x i ,y j ,z0)) = Cop block-2 (Q r，p,q+1 (ξ r ,η p ,ζ q+1 );O i，j (x i ,y j ,z0))

[0050] The amount of calculation can be further reduced by about half on this basis;

[0051] Using equations (1) - (6) to perform forward calculation on the kernel function can significantly improve the forward calculation speed of the kernel function, thus improving the inversion efficiency.

[0052] S4. Use the conjugate gradient method for inversion and utilize the fast Fourier transform to obtain the gravity physical property distribution in the inversion area.

[0053] In actual calculations, under the same meshing conditions, the embodiments of the present invention not only maintain the same calculation accuracy as the traditional method, but also greatly reduce the calculation time. In the model experiment, when performing the forward calculation of the kernel function with a grid of 50×50×50 on a computer with a CPU main frequency of 2.20 GHz and a memory size of 64 GB in MATLAB, the time used is only 0.05 s, and the calculation speed is increased by about 2×10 4 times, and the inversion result is completely consistent with the traditional inversion method, ensuring the accuracy;

[0054] Using the gravity rapid inversion method provided by the embodiments of the present invention to invert the gravity data in the Qihe area of Shandong effectively improves the inversion speed, compresses the forward calculation time to about 14 s, and the inversion speed is increased by 3.7×10 5 or so.

[0055] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A rapid gravity inversion method based on multi - element geometric equivalence, characterized in that, It includes the following steps: Grid the gravity anomaly data, and divide the underground area with cuboids of equal length in the x and y directions; Apply the point equivalent method to calculate the relationship function between the nodes of the underground dissection grid and the surface observation points, and calculate the kernel function matrix between a single hexahedron dissection unit with node equivalent response and the surface observation points; By applying the volume equivalent relationship, extend the response of the dissection unit equivalently to the initial layer, obtain the kernel function matrix of all dissection units in this layer for the surface observation points, and equivalently transfer the response of the nodes in this layer to the next layer, and calculate layer by layer to obtain the kernel function matrix of the whole area; Adopt the conjugate gradient method for inversion, and use the fast Fourier transform to obtain the gravity physical property distribution of the inversion area.

2. The gravity rapid inversion method based on multivariate geometric equivalence according to claim 1, characterized in that In the step of gridding the gravity anomaly data and dividing the underground area with cuboids of equal length in the x and y directions, the distribution of the dissection grid is consistent with the terrain undulation, and the observation points and the dissection units are in one-to-one correspondence.

3. The method for rapid gravity inversion based on multivariate geometric equivalence according to claim 1, characterized in that: The step of applying the point equivalent method to calculate the relationship function between the nodes of the underground dissection grid and the surface observation points, and calculating the kernel function matrix between a single hexahedron dissection unit with node equivalent response and the surface observation points specifically includes: Define the functional relationship between the nodes of a single dissection unit and the observation points as: Among them, O i，j represents the observation point at position (i, j) in the observation point grid, and its coordinates are (x i , y j , z0); Q r，p，q represents the node at coordinates (ξ r , η p , ζ q ) in the subdivision element; R represents the distance from the subdivision node Q r，p，q to the observation point O i，j . The value ranges of r, p, and q are all 0 or 1; For four nodes on the same horizontal plane, there is an equivalent relationship, and the formula is expressed as: There is an identity relationship, and the formula is expressed as: Cop(Q r，p，q (ξ r ,η p ,ζ q )); O i，j (x i ,y j ,z0)) = Cop(Q p，r，q (ξ p ,η r ,ζ q )); O j，i (x j ,y i ,z0)); The observation point directly above the dissection unit has the following relationship with the four nodes on the same horizontal plane of this unit: Cop(Q r，p，q (ξ r ,η p ,ζ q );O i，j (x i ,y j ,z0)) + Cop(Q r+1，p，q (ξ r+1 ,η p ,ζ q );O i，j (x i ,y j ,z0)) + Cop(Q r，p+1，q (ξ r ,η p+1 ,ζ q );O i，j (x i ,y j ,z0)) + Cop(Q r+1，p+1，q (ξ r+1 ,η p+1 ,ζ q );O i，j (x i ,y j ,z0)) = 0; Combine the calculated functional relationship into the kernel function of this dissection block and the observation points as:

4. The gravity rapid inversion method based on multi - element geometric equivalence according to claim 3, characterized in that, The step of applying the volume equivalent relationship to equivalently extend the response of the dissection unit to the initial layer, obtain the kernel function matrix of all dissection units in this layer for the surface observation points, and equivalently transfer the response of the nodes in this layer to the next layer, and calculate layer by layer to obtain the kernel function matrix of the whole area specifically includes: The coordinates of the center point P of the subdivision unit are expressed as (α u , β v , γ w ). The response of the subdivision unit at the observation point I i，j (x i , y j , z0) is M(P u，v (α u , β v , γ w );O i，j (x i , y j , z0)); There is the translational equivalence of volume equivalence, and the formula is expressed as: M(P 1，1 (α1, β1, γ w )); Q i，j (x i , y j , z0)) = M(P u，v (α u , β v , γ w )); Q i+u-1，j+v-1 (x i+u-1 , y j+v-1 , z0)); There is cross equivalence, and the formula is expressed as: M(P u，v (α u ,β v ,γ w )); Q i，j (x j ,y j ,z0)) = M(P i，j (x i ,y j ,γ w )); Q u，v (α u ,β v ,z0)); Output the calculated M g The response value is a row vector and is expanded in the form of a Toeplitz matrix to obtain the response of all subdivision units in this layer to the observation points; When calculating the next layer, inherit the relationship function of some nodes to the next layer: Cop block-1 (Q r，p，q+1 (ξ r ,η p ,ζ q+1 );O i，j (x i ,y j ,z0)) = Cop block-2 (Q r，p，q+1 ( ξ r,η p ,ζ q+1 );O i，j (x i ,y j ,z0)).

5. A gravity rapid inversion system based on multi - element geometric equivalence, which is used to implement the gravity rapid inversion method based on multi - element geometric equivalence as described in any one of claims 1 - 4, and is characterized in that It includes: A dissection module for gridding the gravity anomaly data and dividing the underground area with cuboids of equal length in the x and y directions; A point equivalent calculation module for applying the point equivalent method to calculate the relationship function between the nodes of the underground dissection grid and the surface observation points, and calculating the kernel function matrix between a single hexahedron dissection unit with node equivalent response and the surface observation points; A volume equivalent technology module for equivalently extending the response of the dissection unit to the initial layer by applying the volume equivalent relationship, obtaining the kernel function matrix of all dissection units in this layer for the surface observation points, and equivalently transferring the response of the nodes in this layer to the next layer, and calculating layer by layer to obtain the kernel function matrix of the whole area; An inversion module for adopting the conjugate gradient method for inversion and using the fast Fourier transform to obtain the gravity physical property distribution of the inversion area.

Images