A gravity field correction method and system, device and medium

Abstract

The present specification discloses a gravity field correction method, the method comprises: obtaining a three-dimensional data model of a building group in a target area; the three-dimensional data model is a model based on a tetrahedral grid structure; processing the three-dimensional data model based on a three-dimensional forward numerical simulation method to obtain a gravity correction value of surrounding buildings at a gravity measurement point to the gravity measurement point; based on the gravity correction value, the actual gravity field value at the gravity measurement point is obtained, which solves the problems of low automation, complex field measurement and low terrain correction accuracy in the prior art.

Description

TECHNICAL FIELD

[0001] The present specification relates to the technical field of urban micro-gravity exploration, and in particular to a gravity field correction method and system, device and medium. BACKGROUND

[0002] With the continuous development of urban construction and the continuous increase of environmental noise, the methods such as shallow seismic method, high-density electrical method and electromagnetic sounding method, which have been applied well in the past, are more and more restricted by urban construction sites and cannot be carried out. Micro-gravity exploration method has the advantages of low exploration cost, high construction efficiency, small environmental interference and less construction condition restriction, and should be the preferred method in three-dimensional geophysical exploration. However, the complex "terrain-like" gravity field interference caused by the urban unique building groups such as high-rise buildings, subways and pipe networks makes it difficult for high-precision gravity method to meet the basic requirement of total gravity observation accuracy (±40×10-8m·s-2) in urban engineering geophysical exploration specification (CJJ7-2016) when it is carried out in the city.

[0003] Therefore, how to quantitatively calculate the influence of various building groups in the city on the gravity field value and correct it is a key scientific problem that must be solved for current urban micro-gravity exploration method. SUMMARY

[0004] The present specification discloses a gravity field correction method, which comprises: obtaining a three-dimensional data model of a building group in a target area; the three-dimensional data model is a model based on a tetrahedral grid structure; processing the three-dimensional data model based on a three-dimensional forward numerical simulation method to obtain a gravity correction value of buildings around a gravity measurement point on the gravity measurement point; and obtaining an actual gravity field value at the gravity measurement point based on the gravity correction value.

[0005] In some embodiments, the obtaining of the three-dimensional data model of the building group in the target area comprises: obtaining satellite remote sensing data of the target area; obtaining spatial profile and elevation data of the building group based on the satellite remote sensing data; and constructing the three-dimensional data model based on the spatial profile and elevation data.

[0006] In some embodiments, the obtaining the spatial profile and the elevation data of the building group based on the satellite remote sensing data comprises: preprocessing the satellite remote sensing data to obtain processed data, the preprocessing comprising at least one of image fusion, self-defined coordinate system, orthorectification, and atmospheric correction; adjusting a segmentation scale and a merging scale corresponding to the processed data, and performing data segmentation on the processed data based on the segmentation scale to obtain segmented data; obtaining data features corresponding to the segmented data, performing band merging on the segmented data according to band features corresponding to the data features based on the merging scale, to obtain merged data; selecting data samples from the merged data; and obtaining the spatial profile and the elevation data of the building group based on sample statistical results of the data samples.

[0007] In some embodiments, the constructing the three-dimensional data model based on the spatial profile and the elevation data comprises: constructing an initial data model based on the spatial profile and the elevation data; clustering buildings included in the initial data model based on graphic attributes corresponding to each building to obtain a plurality of clustering clusters; wherein the buildings under one clustering cluster correspond to one building type, and the graphic attributes comprise at least one of surface area and image features; for each clustering cluster, obtaining a density value of the buildings included in the clustering cluster, and taking the density value as a density attribute of the buildings in the clustering cluster; and assigning density values to the buildings under each clustering cluster in the initial model to obtain the three-dimensional data model with the density attribute.

[0008] In some embodiments, the obtaining, for each clustering cluster, the density value of the buildings included in the clustering cluster comprises: for each clustering cluster, taking one of the buildings as a building sample; for the building sample, determining a plurality of candidate densities through density assumption; performing three-dimensional forward fitting based on the candidate densities to obtain correction value simulation curves corresponding to the plurality of candidate densities; performing field measurement on each building sample to obtain a correction value measured curve corresponding to each building sample; taking the candidate density corresponding to the correction value simulation curve with the greatest similarity to the correction value measured curve as the density value of the building sample; and taking the density value of the building sample as the density values of all the buildings in the clustering cluster where the building sample is located.

[0009] In some embodiments, the processing the three-dimensional data model based on the three-dimensional forward numerical simulation method to obtain a gravity correction value of buildings around a gravity measurement point at the gravity measurement point comprises: obtaining, based on a preset algorithm, gravity correction values of each tetrahedral element in the three-dimensional data model at the gravity measurement point; and accumulating the gravity correction values corresponding to all the tetrahedral elements within a target range from the gravity measurement point to obtain the gravity correction value of the buildings around the gravity measurement point at the gravity measurement point.

[0010] In some embodiments, the obtaining, based on a preset algorithm, of the gravity correction value of each tetrahedron unit in the three-dimensional data model at the gravity measurement point comprises:

[0011] For any one tetrahedron unit ABCD, the formula (1) is used: to determine the gravity correction value of the tetrahedron unit at the gravity measurement point; wherein the coordinates of the four vertices of the tetrahedron unit ABCD are represented as: A(x1, y1, z1), B(x2, y2, z2), C(x3, y3, z3) and D(x4, y4, z4); Δg represents the gravity correction value, G represents the universal gravitational constant; σ represents the density value of the building corresponding to the tetrahedron unit; (x, y, z) represents the coordinates of the center of gravity of the tetrahedron unit ABCD, and the center of gravity of the tetrahedron unit ABCD can be determined based on the coordinates of the four vertices; the three-dimensional volume integral of formula (1) is converted into a two-dimensional curved surface integral by using the Gauss formula wherein γ is the normal vector of ∑ at the point (x, y, z); since the normal vectors of the four triangles of the tetrahedron unit ABCD are all a solid constant, formula (1) is changed into formula (2) by combining the Gauss formula: wherein j represents the jth triangle of the tetrahedron unit ABCD, γ j the normal vector of each triangle of the tetrahedron unit ABCD is obtained from the coordinates of the vertices of the triangle.

[0012] The application further discloses a gravity field correction system, comprising a data acquisition module, a data processing module and a data calculation module; the data acquisition module is configured to acquire a three-dimensional data model of a building group in a target area; the three-dimensional data model is a model based on a tetrahedron grid structure; the data processing module is configured to process the three-dimensional data model based on a three-dimensional forward numerical simulation method to obtain a gravity correction value of buildings around a gravity measurement point at the gravity measurement point; and the data calculation module is configured to obtain an actual gravity field value at the gravity measurement point based on the gravity correction value.

[0013] The application further discloses a gravity field correction device comprising a processor, wherein the processor is used to execute the gravity field correction method.

[0014] The application further discloses a computer readable storage medium, wherein the storage medium stores computer instructions, and when a computer reads the computer instructions in the storage medium, the computer executes the gravity field correction method.

[0015] Beneficial effects: Based on the gravity field correction method of the present application, the low-efficiency working mode of traditional manual measurement is overcome, and more accurate building group boundary coordinate and elevation data can be obtained. In addition, the irregular tetrahedral subdivision technology is used instead of the traditional artificial approximate calculation, and for complex building, the subdivision precision of tetrahedron can be realized by continuous encryption. Finally, the high-precision quantitative calculation of urban building group gravity field on the computer platform can be realized through the integration of related core algorithms, and the problems of low automation, complex field measurement and low terrain correction precision in the previous methods are solved. BRIEF DESCRIPTION OF DRAWINGS

[0016] The present specification will be further illustrated in the form of exemplary embodiments, which will be described in detail with reference to the accompanying drawings. These embodiments are not limiting, and in these embodiments, the same numbers represent the same structures, wherein:

[0017] Figure 1 is a module structure schematic diagram of the gravity field correction system according to some embodiments of the present specification;

[0018] Figure 2 is an exemplary flowchart of the gravity field correction method according to some embodiments of the present specification;

[0019] Figure 3 is a schematic diagram of the three-dimensional digital elevation model of urban buildings according to some embodiments of the present specification;

[0020] Figure 4 is a schematic diagram of building quasi-prism subdivision according to some embodiments of the present specification;

[0021] Figure 5 is a schematic diagram of standard triangular prism and quasi-prism according to some embodiments of the present specification;

[0022] Figure 6 is a schematic diagram of subdividing quasi-prism units into tetrahedral units according to some embodiments of the present specification;

[0023] Figure 7 is a schematic diagram of tetrahedral unit ABCD and gravity measurement point P according to some embodiments of the present specification;

[0024] Figure 8 is a schematic diagram of rotating the triangular band ABC in the tetrahedral unit ABCD in the plane rectangular coordinate system according to some embodiments of the present specification;

[0025] Figure 9 is a location plan of the micro-gravity exploration profile of the test area according to some embodiments of the present specification;

[0026] Figure 10 is a WorldView-II remote sensing image of the test area according to some embodiments of the present specification;

[0027] Figure 11 is a DSM plan of the test area according to some embodiments of the present specification;

[0028] Figure 12 is a model map of surface buildings of the test area according to some embodiments of the present specification;

[0029] Figure 13 is a schematic diagram of the main lithologic resistivity and shear wave velocity characteristics of the test area according to some embodiments of the present specification;

[0030] Figure 14 is a schematic diagram of the high-precision gravity exploration results of the Sumei fault of the test area according to some embodiments of the present specification;

[0031] Figure 15 is a schematic diagram of the shallow seismic detection results of the Sumei fault of the test area according to some embodiments of the present specification;

[0032] Figure 16 is a schematic diagram of the comparison of the equivalent reverse magnetic flux transient electromagnetic and shallow seismic detection results of the Sumei fault of the test area according to some embodiments of the present specification; Figure 1 ;

[0033] Figure 17 is a schematic diagram of the comparison of the equivalent reverse magnetic flux transient electromagnetic and shallow seismic detection results of the Sumei fault of the test area according to some embodiments of the present specification; Figure 2 ;

[0034] Figure 18 is a schematic diagram of the comparison of the microseismic and equivalent reverse magnetic flux transient electromagnetic detection results of the Sumei fault of the test area according to some embodiments of the present specification; Figure 1 ;

[0035] Figure 19 is a schematic diagram of the comparison of the microseismic and shallow seismic detection results of the Sumei fault of the test area according to some embodiments of the present specification; Figure 2 ;

[0036] Figure 20 is a schematic diagram of the main interface of the work area data management module of the computer application software for microgravity exploration work according to some embodiments of the present specification;

[0037] Figure 21 is a schematic diagram of the main interface of the surface feature intelligent extraction module of the computer application software for microgravity exploration work according to some embodiments of the present specification;

[0038] Figure 22 According to some embodiments of the present specification, a schematic diagram of a pseudo-triangular prism and its corresponding special case is shown.

[0039] Figure 23 According to some embodiments of the present specification, a schematic diagram of a funnel model cross section is shown. DETAILED DESCRIPTION

[0040] In order to more clearly illustrate the technical solutions of the embodiments of the present specification, the drawings needed to be used in the embodiment description will be briefly introduced as follows. Obviously, the drawings in the following description are only some examples or embodiments of the present specification, and for those skilled in the art, the present specification can also be applied to other similar scenarios without creative labor on the basis of these drawings. Unless it is clear from the language environment or otherwise stated, the same reference numbers in the drawings represent the same structures or operations.

[0041] It should be understood that the "system", "device", "unit" and / or "module" used herein is a method for distinguishing different components, elements, parts, sections or assemblies at different levels. However, if other words can achieve the same purpose, the words can be replaced by other expressions.

[0042] As shown in the specification and claims, unless the context clearly indicates otherwise, "a", "one", "an", and / or "the" do not refer to the singular, but can also include the plural. Generally speaking, the terms "comprise" and "include" only indicate the inclusion of the steps and elements explicitly identified, and these steps and elements do not constitute an exclusive list, and the method or device can also include other steps or elements.

[0043] There is no effective terrain correction method for urban building groups at present. In the past, urban gravity exploration work was carried out in reference to the People's Republic of China Industry Standard for Geological and Mineral Resources "Large-scale Gravity Exploration Specification

DZ / T0171-1997

DZ / T 0004-2015

[0044] In the middle zone 20m-2km range, the "mass body" in the range is simply meshed, and the elevation value of each mesh node is the elevation of the building group on the node. When calculating the middle zone terrain correction, the elevation of the measurement point is shifted to the four adjacent elevation grid nodes, and the middle zone terrain correction value of the four elevation nodes is calculated using the complex trapezoidal integral formula. Finally, the measured point's middle zone terrain correction value is calculated by bilinear interpolation according to the plane position of the measurement point relative to the four adjacent elevation nodes.

[0045] Under the complex working conditions of various buildings in the city, the above traditional method cannot achieve accurate measurement of the elevation data of all building "mass bodies". At the same time, the simple approximate shape substitution calculation method itself has a lot of human approximation, and finally cannot guarantee the effective accuracy of the terrain correction calculation value.

[0046] Therefore, the present specification provides a gravity field correction method and system, which combines the leading high-resolution satellite remote sensing technology, high-precision gravity three-dimensional numerical forward technology and supercomputer cluster parallel computing technology to realize high-precision quantitative calculation of the gravity field of complex buildings in the city.

[0047] Figure 1 It is a module structure schematic diagram of the gravity field correction system according to some embodiments of the present specification.

[0048] As Figure 1 shown, the gravity field correction system 100 (hereinafter referred to as system 100) includes a data acquisition module 110, a data processing module 120, and a data calculation module 130.

[0049] The data acquisition module is configured to obtain a three-dimensional data model of a building group in a target area; the three-dimensional data model is a model based on a tetrahedral grid structure.

[0050] The data processing module is configured to process the three-dimensional data model based on a three-dimensional forward numerical simulation method to obtain a gravity correction value of a gravity measurement point at the gravity measurement point by surrounding buildings;

[0051] The data calculation module is configured to obtain an actual gravity field value at the gravity measurement point based on the gravity correction value.

[0052] In some embodiments, the data acquisition module is further configured to: obtain satellite remote sensing data of the target area; obtain spatial contour and elevation data of the building group based on the satellite remote sensing data; and construct the three-dimensional data model based on the spatial contour and elevation data.

[0053] In some embodiments, the data acquisition module is further configured to: pre-process the satellite remote sensing data to obtain processed data, the pre-processing including at least one of image fusion, self-defined coordinate system, orthorectification, and atmospheric correction; adjust a segmentation scale and a merging scale corresponding to the processed data, and perform data segmentation on the processed data based on the segmentation scale to obtain segmented data; obtain data features corresponding to the segmented data, perform band merging on the segmented data based on the merging scale and band features corresponding to the data features to obtain merged data; select data samples from the merged data; and obtain spatial contours and elevation data of the building group based on sample statistical results of the data samples.

[0054] In some embodiments, the data acquisition module is further configured to: construct an initial data model based on the spatial contours and the elevation data; cluster buildings included in the initial data model based on graphical attributes corresponding to each building to obtain a plurality of clustering clusters; wherein the buildings under one clustering cluster correspond to one building type, and the graphical attributes include at least one of surface area and image features; for each clustering cluster, obtain a density value of the buildings included therein, and take the density value as a density attribute of the buildings in the clustering cluster; and assign density values to the buildings under each clustering cluster in the initial model to obtain the three-dimensional data model with the density attribute.

[0055] In some embodiments, the data acquisition module is further configured to: for each clustering cluster, obtain one of the buildings therein as a building sample; for the building sample, determine a plurality of candidate densities through density assumption; perform three-dimensional forward fitting based on the candidate densities to obtain correction value simulation curves corresponding to the candidate densities; perform field measurement for each building sample to obtain a correction value measured curve corresponding to each building sample; take the candidate density corresponding to the correction value simulation curve with the highest similarity to the correction value measured curve as a density value of the building sample; and take the density value as the density values of all the buildings in the clustering cluster where the building sample is located.

[0056] In some embodiments, the data processing module is further configured to: obtain gravity correction values of each tetrahedral unit in the three-dimensional data model at the gravity measurement point based on a preset algorithm; and accumulate the gravity correction values corresponding to all the tetrahedral units within a target range from the gravity measurement point to obtain a gravity correction value of buildings around the gravity measurement point at the gravity measurement point.

[0057] In some embodiments, the data processing module is further configured to:

[0058] For any tetrahedral element ABCD, using formula (1): Determine the gravity correction value of the tetrahedral element at the gravity measurement point; wherein, the coordinates of the four vertices of the tetrahedral element ABCD are represented as: A(x1,y1,z1), B(x2,y2,z2), C(x3,y3,z3), and D(x4,y4,z4); Δg represents the gravity correction value, G represents the gravitational constant; σ represents the density value of the building corresponding to the tetrahedral element; (x,y,z) represents the coordinates of the centroid of the tetrahedral element ABCD, which can be determined based on the coordinates of the four vertices; using Gauss's formula... The three-dimensional volume integral of formula (1) is converted into a two-dimensional surface integral, where γ is the normal vector of ∑ at point (x,y,z); since the normal vectors of the four triangles of the tetrahedral element ABCD are all solid constants, combined with Gauss's formula, formula (1) is transformed into formula (2): Where j represents the j-th triangular face of the tetrahedral element ABCD, γ j This represents the normal vector of each triangular face of the tetrahedral element ABCD; the normal vector of each triangular face is obtained from the coordinates of the vertices of the triangle.

[0059] For details regarding the above embodiments, please refer to [link / reference]. Figures 2-3 And related explanations.

[0060] Some embodiments of this specification, based on the system 100 of the present invention, can form a set of computer application software for microgravity exploration work by integrating the functions of its various modules and related core algorithms. This software can automatically retrieve the spatial information of all buildings within a preset range (e.g., 0-2km) around the gravity measuring point from the three-dimensional spatial data model of the urban building complex with a TEN structure, based on the coordinate position (x, y, z) of the gravity measuring point. It then automatically performs tetrahedral element subdivision and three-dimensional forward numerical calculation under parallel computing conditions, ultimately obtaining a high-precision gravity terrain correction value for the gravity measuring point location.

[0061] As an example only, the human-computer interaction interface of this software can be designed based on the C# WinForms library, and the software development platform can use Microsoft Visual Studio 2019. The software modules and data processing flow can be divided into 8 sub-modules, such as... Figure 20 As shown, it may include a data warehouse management submodule, a surface topography extraction submodule, various data correction submodules, a survey network profile projection submodule, a potential field transformation processing submodule, a profile 2D inversion submodule, a spatial 3D inversion submodule, and an inversion result mapping submodule, etc.

[0062] Among them, such as Figure 21The surface feature extraction sub-module can load the surface building extraction result data (in the *.shapefile format) of professional remote sensing image processing software such as ENVI and eCognition into the software system, and construct a vector three-dimensional data model of the building group in the test area range through clustering analysis and other algorithms, so as to provide calculation and processing of high-precision terrain correction of the measured gravity data, eliminate the influence of surrounding buildings on the measured point gravity field, and obtain high-precision Bouguer gravity anomaly. The data correction sub-module is responsible for zero drift correction, terrain correction, normal field correction, Bouguer correction and other processing of gravity measurement data. The terrain correction sub-interface will call the aforementioned core algorithm of the application to realize quantitative correction calculation of the gravity field value.

[0063] In some embodiments, in order to improve the algorithm calculation efficiency problem under the fine subdivision of tetrahedral units, the system 100 and its corresponding software further introduce the cooperative heterogeneous parallel computing technology of MPI+GPU.

[0064] MPI (Message Passing Interface) parallel computing technology is a multi-machine cooperative computing method by decomposing a serial task into multiple computing nodes of a parallel computer cluster to execute one or several processes. Its essence is a message passing function library, and its core is the mutual communication between multiple CPUs of multiple hosts. The result is to realize a multi-machine CPU parallel computing cluster.

[0065] GPU (Graphics Processing Unit) parallel computing technology is a method of improving the computing efficiency of a single machine by calling the graphics processor (graphics card) on the computer to join the entire calculation. GPU uses CUDA (Compute Unit Device Architecture) as the software and hardware architecture system of parallel computing devices. Currently, mainstream manufacturers (NVIDIA, Intel, AMD) have implemented support for the OpenCL open standard. GPU devices (graphics cards) of mainstream manufacturers can rely on the OpenCL API tool development package to develop parallel programs. The CPU+GPU programming model of a single machine uses CPU as the host (host) responsible for logical transactions and serial calculations, and uses GPU as the device (device) responsible for executing highly parallel thread tasks.

[0066] The MPI+GPU programming model is a CPU+GPU programming model of multiple single computers combined into a multi-computer CPU+GPU parallel computing cluster through inter-process communication of the CPU.

[0067] It should be noted that the above description of the gravity field correction system and its platform is for convenience of description only and cannot limit the present specification to the scope of the embodiments. It can be understood that, for those skilled in the art, after understanding the principle of the system, any combination of the platforms or connection of the subsystems with other platforms can be made without departing from the principle. In some embodiments, Figure 1 The data acquisition module 110, the data processing module 120, and the data calculation module 130 disclosed in the present specification can be different platforms in a system or can be a platform realizing the functions of two or more than two platforms. For example, the platforms can share a storage database, and the platforms can also have their own storage databases. Variations such as this are within the protection scope of the present specification.

[0068] Figure 2 is an exemplary flowchart of the gravity field correction method according to some embodiments of the present specification. As shown in Figure 2 , the flow 200 includes the following steps 210-230. In some embodiments, the flow 200 can be executed by the system 100.

[0069] Step 210, obtaining a three-dimensional data model of a building group in a target area.

[0070] The target area can be an area where a gravity measurement point is located, such as a predetermined range (e.g., 10 KM) around the gravity measurement point, or an administrative area where the gravity measurement point is located.

[0071] In some embodiments, the three-dimensional data model is a model based on a tetrahedral grid structure.

[0072] In some embodiments, step 210 specifically comprises: obtaining satellite remote sensing data of the target region; obtaining spatial profile and elevation data of the building group based on the satellite remote sensing data; and constructing the three-dimensional data model based on the spatial profile and elevation data.

[0073] In some embodiments, the satellite remote sensing data of the target region can be high-resolution remote sensing data, which can use high-resolution panchromatic images and stereo image pair data, such as WorldView-III data sold by the U.S. DigitalGlobe company for the world, which has a panchromatic image resolution of 0.3 meters and a multispectral image resolution of 1.24 meters. The extraction of the spatial profile of the urban building is realized by using the eCognition remote sensing image processing software based on the object-oriented method.

[0074] In some embodiments, the obtaining of the spatial profile and the elevation data of the building group based on the satellite remote sensing data comprises: pre-processing the satellite remote sensing data to obtain processed data, the pre-processing including at least one of image fusion, self-defined coordinate system, orthorectification, and atmospheric correction; adjusting a segmentation scale and a merging scale corresponding to the processed data, and performing data segmentation on the processed data based on the segmentation scale to obtain segmented data; obtaining data features corresponding to the segmented data, performing band merging on the segmented data according to band features corresponding to the data features based on the merging scale, to obtain merged data; selecting data samples from the merged data; and obtaining the spatial profile and the elevation data of the building group based on sample statistical results of the data samples.

[0075] For example only, the system 100 can realize the extraction of the building elevation information through the DSM (Digital Surface Model) extraction module on the ENVI remote sensing image processing software. The main processing procedure includes: (a) loading the high-resolution stereo image pair data of the target region into the ENVI software, selecting 10-15 evenly distributed field measurement points in the target region as ground control points in the DSM extraction module, and selecting 50 high-precision matching points on the image to extract the DSM; (b) setting the resolution of the target image to 0.5 m x 0.5 m, and outputting the target region DSM image; (c) subtracting the DEM from the DSM of the target region in the ArcGIS software using the layer operation tool to obtain the ground object elevation distribution layer of the target region; (d) spatially superimposing the building profile data extracted before and the ground object height distribution layer to obtain the height of each building object. Finally, the height distribution of the buildings in the target region is obtained in the ArcGIS software as shown in FIG. 2B. Figure 3The shown city building three-dimensional digital elevation model, and use software export function to export the generated results into shapefile format file including building space geometry and attribute information.

[0076] In some embodiments, in order to realize the conversion of the processing results of commercial remote sensing software into a spatial data model for subsequent terrain improvement algorithm calculation, the system 100 can realize the conversion of objective and real buildings into object models that are convenient for computer organization, storage, retrieval and processing, so as to obtain the data structure of each geometric element (volume, surface, arc segment, point) in the object model and the topological relationship between adjacent object models.

[0077] For example, in order to minimize redundancy in spatial entity subdivision, the system 100 first subdivides the building into a combination of several quasi-prismatic bodies when constructing the building object model, such as Figure 4 The figure shows a schematic diagram of a building subdivided into quasi-prismatic bodies. The quasi-prismatic body is different from the standard prismatic body. The quasi-prismatic body is a prismatic body whose upper and lower triangles are not parallel to each other, the side edges are not equal in length and not parallel to each other, and the four vertices of the side quadrilateral are not in the same plane. For example Figure 5 The figure shows a comparison between a standard prismatic body and a quasi-prismatic body. The advantage of the quasi-prismatic body is that it can achieve high-precision approximation of any complex modeling building body with as few quasi-prismatic body elements as possible, while greatly reducing the storage data volume of the entire model.

[0078] In some embodiments, the system 100 can use the TIN formed by the upper and lower triangular sets of the quasi-prismatic body to express the upper and lower top and bottom surfaces of the building, use the quasi-prismatic body set to construct the internal entity between two adjacent top and bottom surfaces, and use the side quadrilateral to describe the spatial relationship between the layers. The quasi-prismatic body-based building three-dimensional modeling is based on the quasi-prismatic body element as the bottom basic unit. The geometric elements of a standard quasi-prismatic body element can be decomposed into: 1 prismatic body, 2 triangles, 3 side quadrilaterals, 3 quadrilateral edges, 6 triangular edges, and 6 vertices.

[0079] In some embodiments, the system 100 can use quasi-prismatic bodies and their special combination bodies to describe any complex modeling city building, such as Figure 22 The figure shows a quasi-prismatic body and its corresponding special case diagram. For example, when the 6 vertices of the standard quasi-prismatic body partially coincide, 4 special cases of quasi-prismatic bodies can be generated, namely Figure 22 Special case (a), special case (b), special case (c), and special case (d) in

[0080] To ensure the precise "fit" building peripheral profile features while trying to save storage space, therefore, in the modeling of urban buildings are like three-prism as a body object. As shown in Figure 6 Figure 1 shows a schematic diagram of the like three-prism body element is divided into tetrahedron element. As shown in Figure 6 Figure 2, when the gravity terrain calculation, can be converted by connecting the like three-prism three side quadrilateral diagonal, 1 like three-prism element into three tetrahedral units, also known as tetrahedral body element (ie TEN structure), thus easily build TEN structure of urban building group three-dimensional spatial data model, realized with "irregular tetrahedral subdivision" three-dimensional high-precision gravity forward calculation seamless link.

[0081] In some embodiments, the based on the spatial profile and elevation data, the construction of the three-dimensional data model comprises: based on the spatial profile and elevation data, the construction of the initial data model; based on the respective building corresponding graphical attributes, the buildings included in the initial data model are clustered to obtain a number of clustering clusters; wherein the buildings under a clustering cluster correspond to a building type, the graphical attributes include at least one of the surface area, the image features; for each clustering cluster, the density value of the building contained therein is obtained, and the density value is taken as the density attribute of the building in the clustering cluster; the buildings under each clustering cluster in the initial model are assigned a density value, and the three-dimensional data model with density attributes is obtained.

[0082] In some embodiments, the system 100 can cluster the attribute information (surface area, image features, etc.) related to the graphical elements in the shapefile format vector data file, and classify the same type of modeling object into the same category. On the basis of a single building object data structure, the average density attribute object of the building group is included in the urban building group object, and finally the urban building data structure with the density attribute is attached.

[0083] For example only, the system 100 can use an open source algorithm library (such as shapelib) to read the attribute information related to the building graphical elements from the *.shx and *.dbf files, merge the buildings with the same or similar attributes into the same category of buildings, and then assign the average density value of different types of buildings to the density object data structure, thereby realizing the construction of the urban building three-dimensional spatial attribute data model.

[0084] In some embodiments, the obtaining, for each cluster, the density value of the building contained therein comprises: for each cluster, obtaining one of the buildings as a building sample; for the building sample, determining a plurality of candidate densities by a density hypothesis; performing a three-dimensional forward fitting based on the candidate densities to obtain a correction value simulation curve corresponding to each candidate density; performing field measurement for each building sample to obtain a correction value measured curve corresponding to each building sample; taking the candidate density corresponding to the correction value simulation curve with the highest similarity to the correction value measured curve as the density value of the building sample; and taking the density value of the building sample as the density value of all the buildings in the cluster where the building sample is located. The horizontal coordinate of the correction value measured curve and the correction value simulation curve is the horizontal distance between the test point and the gravity measured point, and the vertical coordinate is the correction value.

[0085] In step 220, the three-dimensional data model is processed based on a three-dimensional forward numerical simulation method to obtain the gravity correction value of the buildings around the gravity measured point to the gravity measured point.

[0086] In some embodiments, the processing of the three-dimensional data model based on the three-dimensional forward numerical simulation method to obtain the gravity correction value of the buildings around the gravity measured point to the gravity measured point comprises: obtaining the gravity correction value of each tetrahedron unit in the three-dimensional data model at the gravity measured point based on a preset algorithm; and accumulating the gravity correction values of all the tetrahedron units within a target range from the gravity measured point to obtain the gravity correction value of the buildings around the gravity measured point to the gravity measured point.

[0087] In some embodiments, the system 100 can also obtain the average density of a certain type of building group by using the method of field measurement plus three-dimensional forward fitting. For example, the gravity profile measurement is carried out near the isolated single building, and the method of trial and error fitting of measured data is used to determine the average density of the type of building, and the density of other buildings of the same type in the three-dimensional data model is assigned. The method of field measurement plus three-dimensional forward fitting is carried out for each typical type of building in the survey area to obtain the average density and assign the density, so as to realize the construction of the three-dimensional spatial attribute data model of the city building with density attribute.

[0088] As an example only, the spatial geometry and average density data of all buildings within a 0-2km central area can be retrieved from the three-dimensional spatial attribute data model of urban buildings based on the coordinates of the measured points and the radius of the central area. Since the building data model based on the TEN structure is itself a combination of multiple irregular tetrahedral elements, if the ground change value of each tetrahedral element relative to the measuring point is accurately calculated, then the quantitative correction value of the gravity field of the surrounding urban building complex at the measuring point can be obtained by superimposing different tetrahedral elements. The analytical solution calculation method for tetrahedral elements can be found in the gravity calculation formula for uniform density tetrahedral elements derived by M. Okabe.

[0089] In some embodiments, obtaining the gravity correction value of each tetrahedral element in the three-dimensional data model at the gravity measurement point based on a preset algorithm includes:

[0090] For any tetrahedral element ABCD, using formula (1): Determine the gravity correction value of the tetrahedral element at the gravity measurement point; wherein, the coordinates of the four vertices of the tetrahedral element ABCD are represented as: A(x1,y1,z1), B(x2,y2,z2), C(x3,y3,z3), and D(x4,y4,z4); Δg represents the gravity correction value, G represents the gravitational constant; σ represents the density value of the building corresponding to the tetrahedral element; (x,y,z) represents the coordinates of the centroid of the tetrahedral element ABCD, which can be determined based on the coordinates of the four vertices; using Gauss's formula... The three-dimensional volume integral of formula (1) is converted into a two-dimensional surface integral, where γ is the normal vector of Σ at the point (x,y,z); since the normal vectors of the four triangles of the tetrahedral element ABCD are all solid constants, combined with Gauss's formula, formula (1) is transformed into formula (2): Where j represents the j-th triangular face of the tetrahedral element ABCD, γ j This represents the normal vector of each triangular face of the tetrahedral element ABCD; the normal vector of each triangular face is obtained from the coordinates of each vertex of the triangle.

[0091] like Figure 7 The diagram shown is a schematic representation of the tetrahedral element ABCD and the gravity measurement point P. Figure 7 As shown, any tetrahedron ABCD in a rectangular coordinate system has each face as a triangle, and the coordinates of each vertex can be represented as A(x1,y1,z1), B(x2,y2,z2), C(x3,y3,z3), and D(x4,y4,z4). At the measured point P(x... p ,y p ,z p The formula for calculating the gravity change at point () can be expressed as Δg, i.e., Formula 1: Where G is the gravitational constant; σ is the average density of the building.

[0092] Using Gauss's formula The three-dimensional volume integral can be converted into a two-dimensional surface integral, where γ is the normal vector of Σ at the point (x,y,z). Since the normal vectors of the four triangles of the tetrahedron are all solid constants, Equation 1 can be transformed into Equation 2:

[0093]

[0094] Where j represents each face of the tetrahedron, γ j Let represent the normal vectors of each face of the tetrahedron. This formula shows that the volume integral of the tetrahedron can be obtained by multiplying the surface integrals of the triangular faces that enclose the tetrahedron by the direction cosines of the normal vectors of each face.

[0095] The normal vector of each triangular face can be obtained from the coordinates of each vertex of the triangle. Taking the face of △ABC as an example, the corresponding calculation formula is Formula 4:

[0096] In the formula

[0097] For solving the surface integral of a triangle in a rectangular coordinate system, the analytical solution formula can be derived through coordinate system rotation transformation. For example... Figure 8 As shown, taking the plane of △ABC as an example, first, rotate the x-axis and y-axis around the z-axis, and let the rotation be counterclockwise by an angle θ, so that the direction of the x-axis is consistent with the projection direction of the external normal of △ABC onto the xoy plane. Then, rotate the z-axis and the new x-axis around the new y-axis, and let the rotation be clockwise by an angle φ, so that the new z-axis is consistent with the direction of the external normal of △ABC. Thus, the coordinate values ​​under the new coordinates can be calculated by formula 5, where x, y, and z are the coordinates under the original coordinates, and X, Y, and Z are the coordinates under the new coordinates. Formula 5 is as follows:

[0098]

[0099] Among them, cosφ, sinφ, cosθ, and sinθ can all be calculated from the coordinates of each point in △ABC, and the specific calculation formulas are as follows: Formulas 6-8:

[0100] Formula 6:

[0101] Formula 7:

[0102] Formula 8:

[0103] Then rotate the XOY plane coordinate system, assuming it rotates counterclockwise. The new Y axis is consistent with the outer normal direction of an edge, a new coordinate system εoη is obtained, and the coordinate value in the new coordinate system can be calculated by formula 9. wherein,

[0104] In the new coordinate system εoη, the area integral can be converted into a line integral, and formula 10 can be obtained. In formula 10, wherein, v can be obtained by formula 11.

[0105] Formula 11 can be obtained by formula 12 through the operation of partial integral.

[0106]

[0107] Formula 5 is substituted into formula 12 to obtain formula 13.

[0108]

[0109]

[0110] As described above, formula 9 and formula 13 are substituted into formula 10 to obtain the area integral Δg of the triangle. ΔABC .

[0111] In step 230, the actual gravity field value at the gravity observation point is obtained based on the gravity correction value.

[0112] In some embodiments, after obtaining the gravity correction value, the system 100 can adjust and correct the gravity observation value based on the gravity correction value.

[0113] It should be noted that the above description of the process 200 is only for example and illustration, and does not limit the scope of the present specification. Those skilled in the art can make various modifications and changes to the process 200 under the guidance of the present specification. However, these modifications and changes are still within the scope of the present specification.

[0114] The method of the present application is further described below with reference to a specific practical application scenario:

[0115] In some embodiments, to further verify the correctness and accuracy of the algorithm, a funnel model with an accurate theoretical analytical solution is also selected for example comparison. For example, Figure 23 is a schematic view of the cross section of the funnel model, as shown in Figure 23 In the example, the near-zone geodetic correction area is set as a circular area, the geodetic correction radius is 20 m, and the slope is 45°. The difference between the calculation result of the algorithm and the analytical solution under different triangular meshing conditions is analyzed.

[0116] The analytical expression of the funnel model is: Δg 漏斗 = 2πGσR0(1-cosθ), when θ = 45°, R0= 20m, Δg 漏斗 = 0.655852 x 10 -5 m / s 2 . It needs to be further explained that, in order to characterize the fitting degree of the triangular net and the actual terrain, the difference between the normal vector of the triangular net and the normal vector of the real terrain is used for quantification, and the calculation formula is as formula 14: Wherein, r i is the normal vector of the triangular net, and γ i is the normal vector of the real terrain.

[0117] In this funnel model, the real terrain normal vector corresponding to any one of the sub-triangular nets is 45°, and the trial results are shown in Table 1. From the comparison results of the trial, when the sub-triangular net is divided more finely, the fitting error of the triangular net and the real terrain is smaller, and the calculation result of the algorithm is closer to the theoretical value. From the 6th column of Table 2, it can be seen that the calculation accuracy of the algorithm has reached the "micro-g level", and the relative error is less than 1%. As shown in the 5th row of Table 2, when the number of azimuths = 32 and the number of rings = 3, the order of magnitude of the calculation accuracy is less than 1.0 x 10 -8 m / s 2 .

[0118] Table 1 is the comparison results of the trial under the funnel model

[0119]

[0120] In the process of three-dimensional digital modeling of buildings, errors are inevitably generated, therefore, on the basis of the above single-inclined model, Gaussian random noise is simulated and added for trial, the size of the Gaussian noise is set to be 0m in mean value and 1m in standard deviation, and the comparison calculation results are shown in Table 2.

[0121] Table 2 is the comparison results of the trial under the single-inclined model under the Gaussian noise mode.

[0122]

[0123] By comparing the trial results in Table 3, it can be seen that: under the Gaussian measurement error with the mean value of 0m and the standard deviation of 1m, the order of magnitude of the calculation error of the terrain correction of the building can reach 20 x 10 -8 m / s 2 , which can meet the accuracy required in the relevant specifications.

[0124] The following takes the application of the scheme in a specific test area as an example to illustrate the technical effects that can be achieved by the scheme. For example Figure 9The position plan of the micro-gravity exploration profile in the test area is shown in Figure 1, and the gravity data of the profile is shown in Figure 2. Figure 10 The WorldView-II remote sensing image of the test area is shown in Figure 3. Figure 11 The DSM plan of the test area is shown in Figure 4. Figure 12 The model of the surface buildings in the test area is shown in Figure 5.

[0125] During the test, first, three national C-level GPS control points were used to establish gravity and GPS measurement control points in the working area, and the 2000 coordinates (X, Y, Z) of the GPS measurement space points in the working area were calculated through joint measurement and adjustment. Then, the characteristic buildings in the method test area were selected, and the plane coordinates and height of the characteristic points were accurately measured at the intersection of the ground road and the roof of the building by using RTK GPS. Finally, the plane coordinates and height were compared with the spatial geographic data of the ground building group extracted by the intelligent algorithm, and the related errors were corrected.

[0126] During the test, the typical Sumei fault and Longjiangeng fault in the method test area were taken as the research objects, and five parallel test gravity profiles were deployed perpendicular to the fault strike. The average length of the measuring lines was 1.2 kilometers, and the surrounding positions of the measuring lines were densely covered with surface buildings. It was attempted to solve the plane spatial strike distribution and deep underground spatial inclination of the Sumei fault by using gravity anomalies.

[0127] The test area covers various working conditions and building objects such as roads, lakes, rivers, under-construction buildings, and completed residential areas, and has various typical elements of complex urban building groups. The test area belongs to a typical urban under-construction area, and the construction sites are densely distributed in the area. Large and medium-sized engineering vehicles frequently come and go on the main roads for a long time, and the background vibration interference is strong. The surface is densely covered with buildings, and there are completed residential areas and under-construction office buildings, etc. The site is severely restricted. The surface is densely covered with high-voltage power transmission networks, and there is a strong electromagnetic interference background. There is a large area of exposed concrete pavement on the surface, and the grounding is difficult.

[0128] During the specific work, the coordinates (x, y, z) of the measured points were accurately measured by using RTK GPS, and the absolute measurement error of the horizontal and vertical height was required to be less than or equal to 0.05 m. Two different types of gravity meters, CG-5 quartz spring and Burris metal spring, were used to collect gravity data of the measured points, and the instrument data collection parameters (time length, reading mode, etc.) were set according to the optimal results of the data collection quality test under different working conditions. During the test, 3% of the number of inspection points were arranged, and the repeated measurement was carried out by using the method of "same point, different time, different instrument, and different operator", and the error was counted. Finally, the error precision of each item met the precision requirements of the "CJJ7-2016 Urban Engineering Geophysical Exploration Specification".

[0129] The total accuracy and accuracy distribution of the high-precision gravity measurement data correction calculation during the test are shown in Table 3:

[0130] During the test, the urban high-precision gravity measurement data processing software developed based on the core algorithm of the application is used to process the collected 5 measured gravity data profiles, and the processing content includes: (I) pre-processing of instrument collected data; (II) solid tide correction, normal field correction, Bouguer correction and other conventional correction processing; (III) near, medium and far zone topographic correction processing for natural landform features; (IV) near, medium and far zone building correction processing for the influence of surface buildings; (V) gravity field correction processing for known special underground buildings; (VI) anomaly separation processing of regional field and residual field; (VII) two-dimensional and three-dimensional data inversion processing.

[0131] Among them, when the building correction processing for urban building groups is carried out, the three-dimensional spatial attribute data model of the building group of 25 km2 in the method test area is used for calculation. Among them, when the topographic correction processing for natural landform features is carried out, the medium zone georectification calculation (20m-2km) uses the 1:10000 scale digital elevation data (DEM) purchased from the local surveying and mapping management department, and the far zone georectification calculation (2km-166.7km) uses ASTER GDEM V2 data with a resolution of 1"×1" for calculation.

[0132] According to the relevant geological background information of the test area, the surface exposed Quaternary stratum is a sand and gravel pebble layer of river floodplain facies, which is in unconformable contact with the underlying stratum, and the stratum thickness is 20-40m; the underlying stratum is the upper Cretaceous Guankou Group K2g, the lithology is brownish red mudstone and silty mudstone, containing mirabilite and gypsum, and the stratum thickness is >149m. The lithology parameter statistics table of the main stratum is shown in Table 4 and Figure 13 .

[0133] At the same time, shallow earthquakes, equivalent inverse magnetic flux transient electromagnetic sounding and micro-oscillation exploration are deployed for the Sutou buried fault in the test area. Through the collection of engineering geological drilling, core physical properties, seismic and electrical data, the resistivity, density and other physical parameters of the main strata in the test area are obtained. By comparing the actual exploration effects of gravity, electrical and seismic exploration methods on the same exploration line profile, the actual effects of each method are evaluated from the indexes of "exploration cost, work efficiency, effective depth, exploration accuracy and minimum resolution". As Figure 14 shown in the high-precision gravity exploration results of the Sutou fault in the test area. As Figure 14 shown in the density inversion results of the high-precision gravity profile of the L2 line, a low-density anomaly body appears obviously at the position of the Sutou fault, and high-density bodies appear on both sides of the fault position.

[0134] High-resolution reflection wave seismic exploration method was used to conduct shallow seismic detection on the Subaotou buried fault distributed in Chengdu Tianfu New Area. A middle-excitation symmetric receiving observation system was adopted with a receiver interval of 2m, an excitation point interval of 4m, 512 receiving channels, and 128 coverage times. As Figure 15 shown, from the seismic profile, it can be seen that the Subaotou fault is mainly characterized by a pop-up structure. A pop-up structure or thrust triangle structure is formed between the backthrust fault and the thrust fault, which is manifested as an obvious wedge-shaped feature on the seismic profile. F 1-2 is the spatial distribution position of the Subaotou fault. Two or more faults with opposite dips are arranged in a "V" shape or "Λ" shape combination. The "V" shape combination shows that one or more faults are thrusting against a large thrust fault and are hidden beneath it, forming a pop-up structure; the "Λ" shape fault combination is that two faults are thrusting against each other, and a compressive tectonic background of a structural triangle zone is formed in the area bounded by the faults. Affected by the NW-trending main compressive stress, the fault strike is mainly NE, the dip is mainly SE, and the dip of some faults is NW. The Subaotou fault and the buried fault F 1-9 play a decisive role in the overall structural form, forming a "flat fault - slope fault - flat fault" structural style. The secondary faults and the main fault show a "Y" shape combination fault on the section, forming a structural triangle zone and a pop-up structure, resulting in the repetition or absence of some strata.

[0135] Two equivalent anti-flux transient electromagnetic profiles were completed at positions adjacent to the above shallow seismic profile, controlling the Subaotou fault in the area. As Figure 16 and Figure 17 shown, through comprehensive comparative analysis of the exploration results of these survey lines, it is considered that the horizontal throw of this fault is between 10 and 490m, and the throw gradually decreases from southwest to northeast, and the dip is relatively steep at the fault. The fault is mainly distributed in the shallow hill area, and most of the fracture zones are covered by the Quaternary system. However, due to the influence of the fault, traction phenomena often appear in the strata beside the fault, and the strata are upright or inverted. Under the working conditions of the under-construction area with strong interference, equivalent anti-flux transient electromagnetic can collect good data. Combining borehole information and the inverted resistivity contour map can trace and depict geological structures such as the bedrock-cover interface and fault fracture zones within 200m, and at the same time, it also has a certain ability to identify multi-layer mined-out areas within 100m, with good resolution.

[0136] As Figure 18 and Figure 19 shown, two microtremor profile explorations were completed at positions adjacent to the above equivalent anti-flux transient electromagnetic sounding profile. In the case of strong interference, microtremor measurement can use the passive source surface wave background noise for detection imaging, and the purpose of increasing the detection depth and accuracy can be achieved by increasing the radius of the observation array and the number of arrays. Comparing the exploration effects of the equivalent anti-flux transient electromagnetic method and shallow seismic, the microtremor method has a certain ability to identify the overall structure, position, dip, etc. of the fault.

[0137] In summary, through comparison of microgravity, reflection seismic, transient electromagnetic method, microtremor method and other different exploration methods in the vicinity of Sujiaochang fault, it is proved that the application of the present application to urban microgravity exploration is correct, reliable and effective. It is also proved that the microgravity anomaly information obtained by the present application can be used to find out the spatial distribution of the hidden fault in the city.

[0138] The foregoing description has been directed to certain embodiments. This description is not intended to limit the application. Out of an abundance of caution, no disclosure of any kind is intended or should be inferred from any suggestions of alternative aspects or modifications to aspects described herein to realize alternate aspects of the application. It is also to be understood that a variety of modifications and changes can be made to the aspects described and equivalents employed without departing from the spirit of the application. It is expressly intended that all such modifications, changes and equivalents fall within the spirit and scope of the application as described herein.

[0139] Also, the present application has been described with particular terminology to convey the essence of the application. For example, the term "some embodiments" means some, but not necessarily all implementations according to the present application. Thus, use of the phrase some embodiments throughout the specification is not necessarily intended to refer to the same embodiment. Furthermore, the terms "some embodiments" and "one or more embodiments" are used synonymously, unless context dictates otherwise.

[0140] In addition, the order of presentation of the processing elements and sequences, the use of numerical notations, or the use of other designations, are not intended to limit the order of the processes and methods of the present application. Although the above disclosure discusses some presently preferred embodiments of the application by way of example, it is understood that the details as disclosed are merely for the purposes of exemplification and should not be considered limiting. For example, although the system components described above can be implemented by hardware devices, they can also be implemented by software solutions, such as installing the described system on existing servers or mobile devices.

[0141] Similarly, it should be noted that, in order to simplify the presentation of the present disclosure and to aid in the understanding of one or more embodiments of the application, the foregoing description of the embodiments of the present application sometimes refers to multiple features being combined into a single embodiment, drawing or description of the embodiments. In fact, the features of the embodiments are less than all of the features of the single embodiment disclosed.

[0142] In some embodiments, numbers that describe amounts, dimensions, and so forth, are used in the description of the embodiments. It should be understood that such numbers are used only to illustrate examples and that the embodiments can be practiced with other numbers, depending on the application(s) at hand. In some examples, the numbers are modified by the modifier "about" or "approximately" to indicate that the described value allows some leeway. Unless otherwise stated, "about" or "approximately" indicates that the described value allows a ±20% variation. Accordingly, in some embodiments, numerical parameters in the specification are approximations that can vary depending on the desired properties sought to be obtained by the individual embodiments. In some embodiments, numerical parameters are determined by the use of common rounding techniques. Although the numerical ranges and parameters setting forth the broad scope of the embodiments of the application are approximations, the numerical values set forth in the specific examples are reported as precisely as practicable. The numerical values set forth in the specific examples are provided to be as precise as practicable.

[0143] Each patent, patent application, publication, and other material cited in this specification is hereby incorporated by reference in its entirety for the purpose of describing and disclosing, by reference, the chemicals, instruments, articles, or processes described in such patent, patent application, publication, and other material as of the filing date of this patent application. Nothing herein is to be construed as an admission that the application is not entitled to antedate such patent, patent application, publication, and other material by virtue of prior application. In the event of a conflict between the descriptions, definitions, and / or terminology in this specification and that of the above-mentioned materials incorporated by reference, the description, definitions, and / or terminology in this specification controls.

[0144] Finally, it should be understood that the embodiments described herein are merely examples of embodiments of the application. Other variations of the embodiments described herein can also be possible. As such, the alternative configurations of the embodiments of the application, as described and depicted herein, are not to be considered limiting of the scope of the application. Accordingly, the embodiments of the application are not to be considered limited by the embodiments that are described and depicted, but include any alternatives and modifications that fall within the scope of the present application.

Claims

1. A gravity field correction method characterized by, The method comprises: acquiring a three-dimensional data model of a building group in a target area; the three-dimensional data model is a model based on a tetrahedral grid structure; processing the three-dimensional data model based on a three-dimensional forward numerical simulation method to obtain a gravity correction value of buildings around a gravity measurement point at the gravity measurement point; based on the gravity correction value, obtaining an actual gravity field value at the gravity measurement point; the acquisition of the three-dimensional data model of the building group in the target area comprises: acquiring satellite remote sensing data of the target area; the satellite remote sensing data of the target area is high-resolution remote sensing data; based on the satellite remote sensing data, acquiring spatial profile and elevation data of the building group; based on the spatial profile and elevation data, constructing the three-dimensional data model; the acquisition of the spatial profile and elevation data of the building group based on the satellite remote sensing data comprises: preprocessing the satellite remote sensing data to obtain processed data, the preprocessing comprising at least one of image fusion, self-defined coordinate system, orthographic correction, and atmospheric correction; adjusting the segmentation scale and merging scale corresponding to the processed data, and performing data segmentation on the processed data based on the segmentation scale to obtain segmented data; acquiring data features corresponding to the segmented data, and based on the merging scale, performing band merging on the segmented data according to the band features corresponding to the data features to obtain merged data; selecting data samples from the merged data; based on the sample statistical results of the data samples, obtaining the spatial profile and elevation data of the building group; the construction of the three-dimensional data model based on the spatial profile and elevation data comprises: based on the spatial profile and elevation data, constructing an initial data model; based on the graphic attributes corresponding to each building, clustering the buildings included in the initial data model to obtain a plurality of clustering clusters; wherein the buildings under one clustering cluster correspond to one building type, and the graphic attributes include at least one of surface area and image features; for each clustering cluster, acquiring the density value of the buildings included therein, and taking the density value as the density attribute of the buildings in the clustering cluster; assigning the density values of the buildings under each clustering cluster in the initial model to obtain the three-dimensional data model with density attributes; the acquisition of the density value of the buildings included in each clustering cluster comprises: for each clustering cluster, taking one of the buildings therein as a building sample; for the building sample, determining a plurality of candidate densities through density assumption; based on the candidate densities, performing three-dimensional forward fitting to obtain correction value simulation curves corresponding to the candidate densities; for each building sample, performing field measurement to obtain a correction value measured curve corresponding to each building sample; taking the candidate density corresponding to the correction value simulation curve with the greatest similarity to the correction value measured curve as the density value of the building sample; for the clustering cluster in which the building sample is located, the density values of all the buildings therein are all taken as the density value.

2. The method of claim 1, wherein, The three-dimensional forward numerical simulation method is used to process the three-dimensional data model to obtain the gravity correction value of the surrounding buildings of the gravity measurement point at the gravity measurement point. The gravity correction value of each tetrahedron unit in the three-dimensional data model at the gravity measurement point is obtained based on a preset algorithm. The gravity correction values corresponding to all tetrahedron units within the target range from the gravity measurement point are accumulated to obtain the gravity correction value of the surrounding buildings of the gravity measurement point at the gravity measurement point.

3. The method of claim 2, wherein, The gravity correction value of each tetrahedron unit in the three-dimensional data model at the gravity measurement point is obtained based on a preset algorithm. For any tetrahedron unit ABCD, the formula (1) is used: The gravity correction value of the tetrahedron unit at the gravity observation point is determined. Wherein, the coordinates of the four vertices of the tetrahedron unit ABCD are represented as: A(x1, y1, z1), B(x2, y2, z2), C(x3, y3, z3) and D(x4, y4, z4); Δg represents the gravity correction value, G represents the gravitational constant; σ represents the density value of the building corresponding to the tetrahedron unit; (x, y, z) represents the coordinates of the center of gravity of the tetrahedron unit ABCD, which can be determined based on the coordinates of the four vertices; Using the Gauss formula The three-dimensional volume integral of equation (1) is converted to a two-dimensional surface integral, where γ is the normal vector to ∑ at point (x, y, z). Since the normal vectors of the four triangles of the tetrahedron unit ABCD are all a solid constant, formula (1) is changed to formula (2) by combining the Gauss formula: where j denotes the jth triangular face of the tetrahedral unit ABCD, γ j denotes the normal vector of each triangular face of the tetrahedral unit ABCD; the normal vector of each triangular face is calculated from the coordinates of the vertices of the triangular face.

4. A gravity field correction system characterized by, It comprises a data acquisition module, a data processing module and a data calculation module. The data acquisition module is configured to obtain a three-dimensional data model of a building group in a target area; the three-dimensional data model is a model based on a tetrahedron grid structure; The data processing module is configured to process the three-dimensional data model based on a three-dimensional forward numerical simulation method to obtain a gravity correction value of surrounding buildings of a gravity measurement point at the gravity measurement point; The data calculation module is configured to obtain an actual gravity field value at the gravity measurement point based on the gravity correction value; The three-dimensional data model of the building group in the target area comprises: Satellite remote sensing data of the target area is obtained; the satellite remote sensing data of the target area is high-resolution remote sensing data; Based on the satellite remote sensing data, the spatial profile and elevation data of the building group are obtained; Based on the spatial profile and elevation data, the three-dimensional data model is constructed; The satellite remote sensing data is preprocessed to obtain processed data, and the preprocessing includes at least one of image fusion, user-defined coordinate system, orthographic correction, and atmospheric correction; The segmentation scale and the merging scale corresponding to the processed data are adjusted, and the processed data is segmented based on the segmentation scale to obtain segmented data; The data features corresponding to the segmented data are obtained, and based on the merging scale, the segmented data is merged according to the band features corresponding to the data features to obtain merged data; Data samples are selected from the merged data; Based on the sample statistical results of the data samples, the spatial profile and elevation data of the building group are obtained; ​ constructing the three-dimensional data model based on the spatial profile and elevation data comprises: constructing an initial data model based on the spatial profile and elevation data; clustering buildings included in the initial data model based on graphic attributes corresponding to respective buildings to obtain a plurality of clustering clusters; wherein buildings under one clustering cluster correspond to one building type, and the graphic attributes include at least one of surface area and image features; for each clustering cluster, obtaining a density value of buildings included in the clustering cluster, and taking the density value as a density attribute of buildings in the clustering cluster; assigning density values to buildings under respective clustering clusters in the initial model to obtain the three-dimensional data model with density attributes; the obtaining of the density value of the buildings included in each clustering cluster comprises: for each clustering cluster, taking one of the buildings as a building sample; for the building sample, determining a plurality of candidate densities through density assumption; based on the candidate densities, performing three-dimensional forward fitting to obtain correction value simulation curves corresponding to the candidate densities; for each building sample, performing field measurement to obtain a correction value measured curve corresponding to each building sample; taking the candidate density corresponding to the correction value simulation curve with the greatest similarity to the correction value measured curve as the density value of the building sample; for the clustering cluster in which the building sample is located, the density values of all the buildings are all taken as the density value.

5. A gravity field correction device comprising a processor, characterized in that, The processor is configured to execute the gravity field correction method according to any one of claims 1-3.

6. A computer-readable storage medium storing computer instructions, wherein, When a computer reads computer instructions in a storage medium, the computer executes the gravity field correction method according to any one of claims 1-3.

Images