Gravitational field correction method, system and device and medium

Abstract

The invention discloses a gravity field correction method. The method comprises the following steps: 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 a building around the gravity actual measurement point to the gravity actual measurement point; and obtaining an actual gravity field value at the gravity actual measurement point based on the gravity correction value, and the method solves the defects of low automation degree, complex field measurement work and low terrain correction precision of a conventional method.

Description

Technical Field

[0001] This specification relates to the technical field of urban microgravity exploration, and particularly to a gravity field correction method, system, device, and medium. Background Art

[0002] With the continuous development of urban construction and the increasing environmental noise, methods such as shallow seismic method, high-density electrical method, and electromagnetic sounding method, which had good application effects in the past, are increasingly restricted by urban construction sites and cannot be carried out. The microgravity exploration method has the advantages of low exploration cost, high construction efficiency, little environmental interference, and few construction condition restrictions. It should be the preferred method in three-dimensional geophysical exploration. However, due to the complex "quasi-topography" gravity field interference caused by urban-specific building complexes such as high-rise buildings, subways, and pipe networks, the high-precision gravity method often fails to meet the basic requirements of the total gravity observation accuracy (±40×10-8m·s-2) in the Urban Engineering Geophysical Exploration Code (CJJ7-2016) when working in cities.

[0003] Therefore, how to quantitatively calculate the influence of various urban building complexes on the gravity field value and correct it is a key scientific problem that must be solved in the current urban microgravity exploration method. Summary of the Invention

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

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

[0006] In some embodiments, obtaining the spatial contour and elevation data of the building complex based on the satellite remote sensing data includes: preprocessing the satellite remote sensing data to obtain processed data, where the preprocessing includes at least one of image fusion, custom coordinate system, orthorectification, 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; obtaining the 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; and obtaining the spatial contour and elevation data of the building complex based on the sample statistical results of the data samples.

[0007] In some embodiments, constructing the three-dimensional data model based on the spatial contour and elevation data includes: constructing an initial data model based on the spatial contour and elevation data; clustering the buildings included in the initial data model based on the graphic attributes corresponding to each building to obtain a number of clustering clusters; where 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, obtaining the density value of the buildings it contains, and taking the density value as the density attribute of the buildings in the clustering cluster; and assigning the density value to the buildings under each clustering cluster in the initial model to obtain the three-dimensional data model with density attributes.

[0008] In some embodiments, obtaining the density value of the buildings it contains for each clustering cluster includes: for each clustering cluster, obtaining one building as a building sample; for the building sample, determining multiple candidate densities through density hypothesis; performing three-dimensional forward fitting based on the candidate densities to obtain calibration value simulation curves corresponding to the multiple candidate densities; performing field measurement for each building sample to obtain a calibration value measured curve corresponding to each building sample; taking the candidate density corresponding to the calibration value simulation curve with the highest similarity to the calibration value measured curve as the density value of the building sample; and for the clustering cluster where the building sample is located, the density values of all its buildings are taken as the density value.

[0009] In some embodiments, processing 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 measurement point at the gravity measurement point includes: obtaining the gravity correction value of each tetrahedral unit in the three-dimensional data model at the gravity measurement point based on a preset algorithm; and accumulating the gravity correction values corresponding to all the tetrahedral units 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, obtaining the gravity correction value of each tetrahedral element in the three-dimensional data model at the actual gravity measurement point based on a preset algorithm includes:

[0011] For any tetrahedral element ABCD, use formula (1): Determine the gravity correction value of the tetrahedral element at the actual gravity measurement point; where the coordinates of the four vertices of the tetrahedral element ABCD are expressed as: A(x 1 , y 1 , z 1 ), B(x 2 , y 2 , z 2 ), C(x 3 , y 3 , z 3 ), and D(x 4 , y 4 , z 4 ); Δ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, and the centroid of the tetrahedral element ABCD can be determined based on the coordinates of the four vertices; use Gauss's formula Convert the three-dimensional volume integral of formula (1) 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 a solid constant, combined with Gauss's formula, formula (1) is changed to formula (2): Where j represents the jth triangular face of the tetrahedral element ABCD, and γ j Represents the normal vectors of the respective triangular faces of the tetrahedral element ABCD; the normal vectors of the respective triangular faces are obtained from the coordinates of the vertices of the triangle.

[0012] The present invention also discloses a gravity field correction system, including 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 complex in a target area; the three-dimensional data model is a model based on a tetrahedral 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 the gravity correction value of the buildings around the actual gravity measurement point at the actual gravity measurement point; the data calculation module is configured to obtain the actual gravity field value at the actual gravity measurement point based on the gravity correction value.

[0013] The present invention also discloses a gravity field correction device, including a processor, and the processor is used to execute the above-mentioned gravity field correction method.

[0014] The present invention also discloses a computer-readable storage medium storing computer instructions. When a computer reads the computer instructions in the storage medium, the computer executes the gravity field correction method as described above.

[0015] Beneficial effects: Based on the gravity field correction method of the present invention, it overcomes the inefficient working mode of traditional manual measurement, and can obtain more accurate boundary coordinates and elevation data of the building complex. In addition, the use of the irregular tetrahedron subdivision technology replaces the traditional artificial approximate calculation, and for buildings with complex shapes, it can be achieved by continuously densifying the tetrahedron subdivision. Finally, through the integration of relevant core algorithms, high-precision quantitative calculation of the gravity field of urban building complexes can be realized under the computer platform, solving the deficiencies of the previous methods such as low automation level, complex field measurement work, and low terrain correction accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] This specification will further illustrate in the form of exemplary embodiments, and these exemplary embodiments will be described in detail through the drawings. These embodiments are not restrictive. In these embodiments, the same numbers represent the same structures, where:

[0017] Figure 1 is a schematic diagram of the module structure of the gravity field correction system shown in some embodiments of this specification;

[0018] Figure 2 is an exemplary flowchart of the gravity field correction method shown in some embodiments of this specification;

[0019] Figure 3 is a schematic diagram of the three-dimensional digital elevation model of urban buildings shown in some embodiments of this specification;

[0020] Figure 4 is a schematic diagram of the dissection of a building into a prismoid-like shape shown in some embodiments of this specification;

[0021] Figure 5 is a schematic diagram of a standard triangular prism and a prismoid-like shape shown in some embodiments of this specification;

[0022] Figure 6 is a schematic diagram of dissecting a prismoid-like element into tetrahedron units shown in some embodiments of this specification;

[0023] Figure 7 is a schematic diagram of a tetrahedron unit ABCD and a gravity measurement point P shown in some embodiments of this specification;

[0024] Figure 8 is a schematic diagram of rotating the triangular belt ABC in the tetrahedron unit ABCD in a plane rectangular coordinate system shown in some embodiments of this specification;

[0025] Figure 9 is a location plan of the microgravity exploration profile in the test area shown in some embodiments of this specification;

[0026] Figure 10 is a WorldView-II remote sensing image map of the test area shown in some embodiments of this specification;

[0027] Figure 11 is a DSM plan view of the test area shown in some embodiments of this specification;

[0028] Figure 12 is a surface building model diagram of the test area shown in some embodiments of this specification;

[0029] Figure 13 is a schematic diagram showing the resistivity and shear wave velocity characteristics of the main strata lithology in the test area shown in some embodiments of this specification;

[0030] Figure 14 is a schematic diagram showing the high-precision gravity exploration results of the Sumadou fault in the test area shown in some embodiments of this specification;

[0031] Figure 15 is a schematic diagram showing the shallow seismic exploration results of the Sumadou fault in the test area shown in some embodiments of this specification;

[0032] Figure 16 is a comparison schematic diagram of the results of the equivalent reverse magnetic flux transient electromagnetic and shallow seismic exploration of the Matou fault in the test area shown in some embodiments of this specification Figure 1 ;

[0033] Figure 17 is a comparison schematic diagram of the results of the equivalent reverse magnetic flux transient electromagnetic and shallow seismic exploration of the Matou fault in the test area shown in some embodiments of this specification Figure 2 ;

[0034] Figure 18 is a comparison schematic diagram of the results of the microtremor and equivalent reverse magnetic flux transient electromagnetic exploration of the Sumadou fault in the test area shown in some embodiments of this specification Figure 1 ;

[0035] Figure 19 is a comparison schematic diagram of the results of the microtremor and shallow seismic exploration of the Sumadou fault in the test area shown in some embodiments of this 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 available for microgravity exploration work shown in some embodiments of this specification;

[0037] Figure 21 It is a schematic diagram of the main interface of the intelligent extraction module for surface features of computer application software that can be used for microgravity exploration work as shown in some embodiments of this specification;

[0038] Figure 22 It is a schematic diagram of a pseudo triangular prism and its corresponding special cases as shown in some embodiments of this specification;

[0039] Figure 23 It is a schematic diagram of the cross-section of a funnel model as shown in some embodiments of this specification. Detailed implementation manners

[0040] To more clearly illustrate the technical solutions of the embodiments of this specification, the accompanying drawings required for the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings in the following description are only some examples or embodiments of this specification. For those of ordinary skill in the art, without creative efforts, this specification can also be applied to other similar scenarios based on these drawings. Unless obvious from the language context or otherwise stated, the same reference numerals in the figures 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, portions or assemblies at different levels. However, if other words can achieve the same purpose, the said words can be replaced by other expressions.

[0042] As shown in this specification and the claims, unless the context clearly indicates an exception, words such as "a", "an", "one" and / or "the" are not specifically singular and may also include plural. Generally speaking, the terms "comprising" and "including" only indicate the inclusion of the clearly identified steps and elements, and these steps and elements do not constitute an exclusive list. The method or device may also include other steps or elements.

[0043] Currently, there is no effective topographic correction method for urban building groups. In the past, urban gravity exploration work was carried out with reference to the "Large Scale Gravity Exploration Specification [DZ / T 0171-1997]" and the "Technical Specification for Gravity Survey (1:50000) [DZ / T 0004-2015]" of the Geological and Mineral Industry Standard of the People's Republic of China. In the range of 0-20m in the near area, a laser rangefinder was used to measure the elevation of each "mass body" within the range of 0-20m at intervals of 8 to 16 azimuths, and approximate regular shapes (cones, fans, inclined planes, steps, etc.) were used to replace and calculate the topographic correction value in the near area.

[0044] In the range of 20m-2km in the middle area, the "mass body" within this range is simply gridded, and the elevation value of each grid node is the elevation of the building complex on the node. When calculating the terrain correction in the middle area, the elevation of the measuring point is translated to the four adjacent elevation grid nodes around it, and the middle area terrain correction values ​​on the four elevation nodes are calculated respectively using the complex trapezoidal integral formula. Finally, according to the plane position of the measuring point relative to the four adjacent elevation nodes, the middle area terrain correction value of the actual measuring point is approximated by bilinear interpolation.

[0045] In the context of various complex buildings in the city, the above traditional methods cannot achieve accurate measurement of the elevation data of the "mass body" of all buildings. At the same time, the simple approximate shape replacement calculation method itself has a lot of artificial approximations, and ultimately cannot guarantee the effective accuracy of the terrain correction calculation value.

[0046] Therefore, this specification provides a gravity field correction method and system, which realizes high-precision quantitative calculation of the gravity field of complex urban buildings by combining cutting-edge high-resolution satellite remote sensing technology, high-precision three-dimensional numerical forward modeling technology of gravity and supercomputing cluster parallel computing technology.

[0047] Figure 1 It is a schematic diagram of the module structure of the gravity field correction system shown in some embodiments of this specification.

[0048] like Figure 1 As 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 the building complex in the 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 the gravity measurement point for buildings surrounding the gravity measurement point;

[0051] The data calculation module is configured to obtain an actual gravity field value at a 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 contours and elevation data of the building complex based on the satellite remote sensing data; and construct the three-dimensional data model based on the spatial contours and elevation data.

[0053] In some embodiments, the data acquisition module is further configured to: preprocess the satellite remote sensing data to obtain processed data, where the preprocessing includes at least one of image fusion, custom coordinate system, orthorectification, and atmospheric correction; adjust the segmentation scale and 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 the data features corresponding to the segmented data, and based on the merging scale, perform band merging on the segmented data according to the band features corresponding to the data features to obtain merged data; select data samples from the merged data; based on the sample statistical results of the data samples, obtain the spatial contour and elevation data of the building complex.

[0054] In some embodiments, the data acquisition module is further configured to: based on the spatial contour and elevation data, construct an initial data model; cluster the buildings included in the initial data model based on the graphic attributes corresponding to each building to obtain several clustering clusters; where 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, obtain the density value of the buildings it contains, and use the density value as the density attribute of the buildings in the clustering cluster; assign density values to the buildings under each clustering cluster in the initial model to obtain the three-dimensional data model with density attributes.

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

[0056] In some embodiments, the data processing module is further configured to: based on a preset algorithm, obtain the gravity correction value of each tetrahedral unit in the three-dimensional data model at the gravity measurement point; accumulate the gravity correction values corresponding to all the tetrahedral units within the 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.

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

[0058] For any tetrahedron element ABCD, use formula (1): Determine the gravity correction value of the tetrahedron element at the gravity measurement point; among them, the coordinates of the four vertices of the tetrahedron element ABCD are expressed as: A(x 1 , y 1 , z 1 ), B(x 2 , y 2 , z 2 ), C(x 3 , y 3 , z 3 ), and D(x 4 , y 4 , z 4 ); Δg represents the gravity correction value, G represents the universal gravitational constant; σ represents the density value of the building corresponding to the tetrahedron element; (x, y, z) represents the coordinates of the centroid of the tetrahedron element ABCD, and the centroid of the tetrahedron element ABCD can be determined based on the coordinates of the four vertices; use Gauss's formula Convert the three-dimensional volume integral of formula (1) 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 element ABCD are all a solid constant, combined with Gauss's formula, formula (1) becomes formula (2): Among them, j represents the j-th triangular face of the tetrahedron element ABCD, and γ j represents the normal vectors of the respective triangular faces of the tetrahedron element ABCD; the normal vectors of the respective triangular faces are obtained from the coordinates of the vertices of the triangle

[0059] For the specific content of the above embodiments, reference can be made to Figures 2 - 3 and related descriptions.

[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 through software integration of relevant core algorithms based on the functions of its respective modules. Through this software, according to the coordinate position (x, y, z) of the gravity measurement point, the spatial information of all buildings within a preset range (such as a range of 0 - 2 km) around the measurement point can be automatically retrieved from the three-dimensional spatial data model of the urban building complex in the TEN structure, and the tetrahedron element can be automatically dissected and the three-dimensional forward numerical calculation under parallel computing conditions can be performed, and finally the high-precision gravity terrain correction value at the gravity measurement point position can be obtained.

[0061] Only as an example, the human-computer interaction interface of this software can be designed based on the C# winform library, the software development platform can use Microsoft Visual Studio 2019, and the software modules and data processing flow can be divided into 8 sub-modules, such asFigure 20 As shown, it may include a data warehouse management sub-module, a surface terrain extraction sub-module, various data correction sub-modules, a survey network profile projection sub-module, a potential field transformation processing sub-module, a two-dimensional profile inversion sub-module, a three-dimensional space inversion sub-module, an inversion result mapping sub-module, etc.

[0062] Among them, as Figure 21 Shown is a schematic diagram of the main interface of the surface feature intelligent extraction module corresponding to the surface terrain extraction sub-module. The surface terrain extraction sub-module can load the surface building extraction result data (*.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 through algorithms such as cluster analysis for the calculation and processing of high-precision terrain correction of the later measured gravity data, eliminating the influence of surrounding buildings on the gravity field of the measured points, and obtaining high-precision Bouguer gravity anomalies. Each data correction sub-module is responsible for the processing of zero drift correction, terrain correction, normal field correction, Bouguer correction, etc. of the gravity measurement data. The terrain correction sub-interface in it will call the aforementioned core algorithm of the present invention to implement the 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 elements, in this system 100 and its corresponding software, a collaborative heterogeneous parallel computing technology of MPI+GPU is further introduced.

[0064] The MPI (Message Passing Interface) parallel computing technology realizes multi-machine collaborative computing 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 hosts and multiple CPU processes. The result is to realize a multi-machine CPU parallel computing cluster.

[0065] GPU (Graphics Processing Unit) parallel computing technology improves the computing efficiency of a single machine by incorporating the graphics processing unit (graphics card) of a computer into the overall computing process. GPU uses CUDA (Compute Unit Device Architecture) as the software and hardware architecture system for parallel computing devices. Currently, mainstream manufacturers (NVIDIA, Intel, AMD) have all implemented support for the OpenCL open standard, and GPU devices (graphics cards) of mainstream manufacturers can all develop parallel programs relying on the API toolkit of OpenCL. The CPU+GPU programming model for a single machine uses the CPU as the host, responsible for transactions with strong logic and serial computing; and uses the GPU as the device, responsible for executing highly parallel threaded tasks.

[0066] The MPI+GPU programming model combines the CPU+GPU programming models of multiple single machines into a multi-machine CPU+GPU parallel computing cluster through communication between CPU processes. The optimization algorithm based on MPI+GPU heterogeneous parallel technology in the present invention reasonably classifies subroutines or functions in the high-precision terrain correction algorithm for gravity. It assigns computations with strong logic (such as retrieval and calculation of building model units) to the CPU and highly parallel computations (3D forward modeling of tetrahedral element sets) to the GPU, and realizes scheduling, operation, storage, communication, etc. between multiple machines and between the CPU and the GPU by calling the APIs in the MPI and GPU software packages, thereby effectively solving the problems of large computational volume, large memory consumption, and long calculation time.

[0067] It should be noted that the above description of the gravity field correction system and its platform is only for convenience of description and does not limit this specification to the scope of the examples given. It can be understood that for those skilled in the art, after understanding the principle of the system, they may, without departing from this principle, make any combination of each platform, or form a subsystem and connect it to other platforms. In some embodiments, Figure 1 the data acquisition module 110, data processing module 120, and data calculation module 130 disclosed in may be different platforms in a system, or a single platform may implement the functions of two or more of the above platforms. For example, each platform may share a storage database, or each platform may have its own storage database separately. Such variations are all within the protection scope of this specification.

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

[0069] Step 210, obtain a three-dimensional data model of the building complex in the target area.

[0070] Among them, the target area may be the area where the gravity measurement point for determining the actual gravity field value is located, such as the area covered by the preset range (such as 10 KM) around the gravity measurement point, or the administrative area where the gravity measurement point is located, etc.

[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 includes: obtaining satellite remote sensing data of the target area; based on the satellite remote sensing data, obtaining the spatial contour and elevation data of the building complex; based on the spatial contour and elevation data, constructing the three-dimensional data model.

[0073] In some embodiments, the satellite remote sensing data of the target area may be high-resolution remote sensing data. High-resolution remote sensing data can use high-resolution panchromatic images and stereo pair data. For example, the WorldView-III data sold globally by DigitalGlobe of the United States, whose panchromatic image resolution is 0.3 meters and whose multispectral image resolution is 1.24 meters. The extraction of the urban building contour is realized by using an object-oriented method on the eCognition remote sensing image processing software.

[0074] In some embodiments, the obtaining the spatial contour and elevation data of the building complex based on the satellite remote sensing data includes: preprocessing the satellite remote sensing data to obtain processed data, and the preprocessing includes at least one of image fusion, custom coordinate system, orthorectification, 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; obtaining the 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 contour and elevation data of the building complex.

[0075] By way of example only, system 100 can extract building elevation information through the DSM (Digital Surface Model) extraction module in ENVI remote sensing image processing software. The main processing procedures include: (a) Loading the high-resolution stereo image pair data of the target area into the ENVI software, selecting 10 - 15 evenly distributed field measurement points within the target area as ground control points in the DSM extraction module, and selecting 50 pairs of high-precision matching points on the images to extract the DSM; (b) Setting the resolution of the target image to 0.5m × 0.5m and outputting the DSM image of the target area; (c) Using the layer operation tool in the ArcGIS software to subtract the DSM of the target area from the DEM to obtain the elevation distribution layer of the surface objects in the target area; (d) Spatially overlaying the previously extracted building contour data with the surface object height distribution layer to obtain the height of each building object. Finally, in the ArcGIS software, a three-dimensional digital elevation model of urban buildings as shown in Figure 3 is obtained, and the generated result is exported in the shapefile format file including the spatial geometry and attribute feature information of the buildings using the software export function.

[0076] In some embodiments, in order to convert the processing results of commercial remote sensing software into a spatial data model for subsequent terrain improvement algorithms, system 100 can convert objectively real buildings into object models that are convenient for computer organization, storage, retrieval, and processing, so as to obtain the data structures of various geometric elements (volumes, surfaces, arcs, points) in the object model and the topological relationships between adjacent object models.

[0077] By way of example only, in order to minimize redundancy in spatial entity subdivision as much as possible, when constructing the building object model, system 100 first dissects the building into a combination of several prismoid-like bodies, as shown in Figure 4 is a schematic diagram of dissecting a building into prismoid-like bodies. Among them, the prismoid-like body is different from the standard triangular prism. The prismoid-like body is a triangular prism with non-parallel upper and lower triangles, unequal and non-parallel side edge lengths, and the four vertices of the side quadrilateral not in the same plane. As shown in Figure 5 is a comparison schematic diagram of the standard triangular prism and the prismoid-like body. The advantage of the prismoid-like body is that it can achieve a high-precision approximation of any complex-shaped building body with as few prismoid-like body elements as possible, and at the same time, it will greatly reduce the storage data volume of the entire model.

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

[0079] In some embodiments, the system 100 can describe any complex urban buildings by using triangular prisms and their special combinations, such as Figure 22 The figure shows a schematic diagram of a triangular prism and its corresponding special cases. For example, when the six vertices of a standard triangular prism overlap, four special cases of triangular prisms can be generated, namely: Figure 22 Special case (a), special case (b), special case (c), special case (d) in .

[0080] In order to ensure the accurate "fitting" of the building's outer contour features while saving storage space as much as possible, a triangular prism is used as the voxel object when modeling urban buildings. Figure 6 The figure shows a schematic diagram of dividing a triangular prism-like element into tetrahedral elements. Figure 6 As shown in the figure, when performing gravity terrain calculations, one triangular prism-like element can be converted into three tetrahedral units, also called tetrahedral elements (i.e., TEN structures), by connecting the diagonals of the three side quadrilaterals of the triangular prism-like element. This can easily construct a three-dimensional spatial data model of urban buildings with a TEN structure, thus achieving a seamless link with the three-dimensional high-precision gravity forward calculation of "irregular tetrahedral partitioning".

[0081] In some embodiments, constructing the three-dimensional data model based on the spatial contour and elevation data includes: constructing an initial data model based on the spatial contour and elevation data; clustering the buildings included in the initial data model based on the graphic attributes corresponding to each building to obtain a plurality of clusters; wherein a building under a cluster corresponds to a building type, and the graphic attributes include at least one of surface area and image features; for each cluster, obtaining the density value of the buildings contained therein, and using the density value as the density attribute of the buildings in the cluster; assigning density values ​​to the buildings under each cluster in the initial model to obtain the three-dimensional data model with density attributes.

[0082] In some embodiments, the system 100 can perform a clustering analysis on the attribute information (such as surface area, image features, etc.) related to graphic elements in a shapefile format vector data file, classify the same type of modeled objects into the same category, incorporate the average density attribute object of the building complex into the urban building complex object based on the data structure of a single building object, and finally attach the density attribute to the data structure of the urban building.

[0083] Merely by way of example, the system 100 can use an open-source algorithm library (such as shapelib) to read the attribute information related to building graphic elements from *.shx and *.dbf files, group the buildings with the same or similar attributes into the same category of buildings, and then assign the average density values of different types of buildings to the density object data structure, thereby realizing the construction of the 3D spatial attribute data model of urban buildings.

[0084] In some embodiments, obtaining the density value of the buildings included in each clustering cluster includes: for each clustering cluster, obtaining one building as a building sample; for the building sample, determining multiple candidate densities through density assumptions; performing 3D forward modeling fitting based on the candidate densities to obtain the corrected value simulation curves corresponding to the multiple candidate densities; conducting field measurements for each building sample to obtain the corrected value measured curves corresponding to each building sample; taking the candidate density corresponding to the corrected value simulation curve with the highest similarity to the corrected value measured curve as the density value of the building sample; for all the buildings in the clustering cluster where the building sample is located, the density values of all the buildings are taken as the density value. Among them, the abscissa of the corrected value measured curve and the corrected value simulation curve is the horizontal distance between the test point and the gravity measurement point, and the ordinate is the corrected value.

[0085] Step 220: Process the 3D data model based on the 3D forward numerical simulation method to obtain the gravity correction value of the buildings around the gravity measurement point at the gravity measurement point.

[0086] In some embodiments, processing the 3D data model based on the 3D forward numerical simulation method to obtain the gravity correction value of the buildings around the gravity measurement point at the gravity measurement point includes: based on a preset algorithm, obtaining the gravity correction value of each tetrahedral unit in the 3D data model at the gravity measurement point; accumulating the gravity correction values corresponding to all the tetrahedral units 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.

[0087] In some embodiments, the system 100 may also adopt the method of on-site measurement plus three-dimensional forward modeling fitting to obtain the average density of a certain type of building complex. For example, on-site measurement of the gravity profile is carried out near an isolated single building, and the method of trial-and-error fitting of the measured data is used to determine the average density of this type of building, and density assignment is performed on other buildings of the same type in the three-dimensional data model. By using the method of on-site measurement plus three-dimensional forward modeling fitting for each typical type of building in the exploration area to obtain the average density and perform density assignment, the construction of a three-dimensional spatial attribute data model of urban buildings with density attributes can be achieved.

[0088] Only as an example, the spatial geometry and average density data of all buildings within the 0-2 km central area can be retrieved from the three-dimensional spatial attribute data model of urban buildings according to the measured point coordinates and the size of the central area radius. Since the building data model based on the TEN structure itself is a combination of multiple irregular tetrahedral elements, if the terrain correction 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 this measuring point can be obtained through the superposition calculation of different tetrahedral elements. The analytical solution calculation method of the tetrahedral element can refer to the gravity calculation formula of the uniform density tetrahedral element derived by M. Okabe.

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

[0090] For any tetrahedral unit ABCD, use formula (1): Determine the gravity correction value of the tetrahedral unit at the gravity measurement point; where the coordinates of the four vertices of the tetrahedral unit ABCD are expressed as: A(x 1 ,y 1 ,z 1 ), B(x 2 ,y 2 ,z 2 ), C(x 3 ,y 3 ,z 3 ), and D(x 4 ,y 4 ,z 4 ); Δg represents the gravity correction value, G represents the gravitational constant; σ represents the density value of the building corresponding to the tetrahedral unit; (x, y, z) represents the coordinates of the centroid of the tetrahedral unit ABCD, and the centroid of the tetrahedral unit ABCD can be determined based on the coordinates of the four vertices; using Gauss's formula Convert the three-dimensional volume integral in formula (1) 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 a solid constant, combined with the Gauss formula, formula (1) is changed to formula (2): where j represents the j-th triangular face of the tetrahedral element ABCD, γ j represents the normal vectors of the respective triangular faces of the tetrahedral element ABCD; the normal vectors of the respective triangular faces are obtained from the coordinates of the vertices of the triangle.

[0091] As Figure 7 shown is a schematic diagram of the tetrahedral element ABCD and the gravity measurement point P. As Figure 7 shown, for any tetrahedron ABCD in the rectangular coordinate system, each of its faces is a triangle, and the coordinates of each vertex can be expressed as A(x 1 , y 1 , z 1 ), B(x 2 , y 2 , z 2 ), C(x 3 , y 3 , z 3 ), and D(x 4 , y 4 , z 4 ). At the gravity terrain correction calculation formula at the measurement point P(x p , y p , z p ) can be expressed as Δg, that is, formula 1: where G is the universal gravitational constant; σ is the average density of the building.

[0092] Using the Gauss formula the three-dimensional volume integral can be converted into a two-dimensional surface integral. Among them, γ in is the normal vector of Σ at the point (x, y, z). Since the normal vectors of the 4 triangles of the tetrahedron are all a solid constant, formula 1 can be changed to formula 2:

[0093]

[0094] where j represents each face of the tetrahedron, γ j represents the normal vectors of each face of the tetrahedron. This formula shows that calculating the volume integral of the tetrahedron can be obtained by multiplying the surface integral of each triangular face enclosing the tetrahedron by the direction cosine of the normal vector of each face.

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

[0096] In the formula

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

[0098]

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

[0100] Formula 6:

[0101] Formula 7:

[0102] Formula 8:

[0103] Then rotate the coordinates of the XOY plane. Suppose it rotates counterclockwise by an angle so that the new Y-axis is consistent with the outer normal direction of a certain side, obtaining a new coordinate system εoη. The coordinate values in the new coordinate system can be calculated by formula 9. Formula 9: Among them,

[0104] In the new coordinate system εoη, the surface integral can be reduced to a line integral, and thus formula 10 can be obtained: In formula 10, Among them, v can be obtained from formula 11. Formula 11:

[0105] Through integration by parts operation on formula 11, formula 12 can be obtained:

[0106]

[0107] Substituting formula 5 into formula 12, formula 13 can be obtained:

[0108]

[0109]

[0110] In summary, substituting Equation 9 and Equation 13 into Equation 10, the area integral Δg of the triangle can be obtained. ΔABC 。

[0111] Step 230: Based on the gravity correction value, obtain the actual gravity field value at the gravity measurement point.

[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 illustration and example, and does not limit the scope of application of this specification. For those skilled in the art, various modifications and changes can be made to the process 200 under the guidance of this specification. However, these modifications and changes are still within the scope of this specification.

[0114] The following takes a specific practical application scenario as an example to further illustrate the method of the present invention:

[0115] In some embodiments, to further verify the correctness and accuracy of the algorithm, a funnel model with an exact theoretical analytical solution is also selected for example comparison. As Figure 23 is a schematic cross-sectional view of the funnel model. As Figure 23 shown, in the example, the near-field terrain correction area is set as a circular domain, the terrain correction radius is 20 m, and the slope is 45°. The differences between the algorithm calculation results and their analytical solutions are analyzed under different triangular mesh subdivisions.

[0116] The analytical formula of the funnel model is: Δg 漏斗 = 2πGσR 0 (1 - cosθ). When θ = 45° and R 0 = 20 m, Δg 漏斗 = 0.655852×10 -5 m / s 2 . It should be further noted that to characterize the fitting degree between the triangular mesh and the actual terrain, the difference between the normal vector of the triangular mesh and the normal vector of the true terrain is used for quantification. The calculation formula is as Equation 14: where r i is the normal vector of the triangular mesh, and γ i is the normal vector of the true terrain.

[0117] In this funnel model, the true terrain normal vector corresponding to any of its sub-triangulation meshes is 45°. The trial calculation results are shown in Table 1. From the comparison results of the trial calculations, it can be seen that when the sub-triangulation mesh is divided finer, the fitting error between the triangulation mesh and the true terrain is smaller, and the algorithm calculation result is closer to its theoretical value. From the 6th column of Table 2, it can be seen that the calculation accuracy of this algorithm has reached the "microgal 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 its calculation accuracy is less than 1.0×10 -8 m / s 2 .

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

[0119] In the process of 3D digital modeling of buildings, errors inevitably exist. Therefore, on the basis of the above-mentioned monoclinic model, Gaussian random noise is simulated and added for trial calculations. The magnitude of the Gaussian noise is set to have a mean of 0 m and a standard deviation of 1 m. The comparison calculation results are shown in Table 2

[0120] Table 2 shows the comparison results of the trial calculations under the monoclinic model in the Gaussian noise mode

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

[0122] Taking the application of this solution to a specific test area as an example, the technical effects that can be achieved by this solution are described. As Figure 9 shown is the location plan of the microgravity exploration profile in the test area. As Figure 10 is the WorldView-II remote sensing image of the test area. As Figure 11 is the DSM plan of the test area. As Figure 12 is the surface building model diagram of the test area

[0123] During the test, first, three national C-level GPS control points near the gravity test area were used as a reference to establish gravity and GPS measurement control base points in the working area. Through combined measurement and adjustment, the 2000 coordinates (X, Y, Z) of the GPS measurement spatial points in the working area were calculated. Based on this, characteristic buildings in the method test area were selected, and the plane coordinates and heights of the characteristic points were accurately measured using RTK GPS at the intersections of ground roads and the rooftops of buildings respectively. Then, data such as the plane coordinates and heights were compared with the spatial geographical data of the ground building complex extracted by the intelligent algorithm to compare and correct relevant errors.

[0124] During the test, the typical Sumadou Fault and Longjiageng Fault in the method test area were taken as the test research objects. Five parallel experimental gravity profiles were deployed perpendicular to the fault strike. The average length of the survey lines was 1.2 km, and there were many surface buildings around the survey lines. An attempt was made to solve the geological problems of the plane spatial distribution and the deep underground dip of the Sumadou Fault through gravity anomalies.

[0125] The test area covers various working conditions and building objects such as roads, lakes, rivers, under-construction buildings, and completed residential areas, with all typical elements of urban complex building groups. The test area is a typical urban under-construction area, with numerous construction sites scattered in the area. Medium and large engineering vehicles travel frequently on the main roads for a long time, resulting in strong background vibration interference. The surface buildings are dense, including completed residential areas and under-construction office buildings, with serious site restrictions. There are many high-voltage transmission lines on the surface, with a large electromagnetic interference background. The concrete pavement is exposed over a large area on the surface, making grounding difficult.

[0126] During the specific work, the coordinates (x, y, z) of the measured points were accurately measured using RTK GPS, with the requirement that the absolute measurement errors in horizontal and vertical heights ≤ 0.05 m. Two different types of gravimeters, namely CG-5 quartz spring and Burris metal spring, were used to collect the gravity data of the measured points. The acquisition parameters of the instrument data (duration, reading method, etc.) were set according to the optimal results of the data acquisition quality tests under different working conditions. During this period, 3% of the checkpoints were arranged, and repeated measurements were carried out in the way of "the same measured point, different times, different instruments, different operators" to statistically analyze the errors. Finally, the accuracy requirements of each error were required to meet the accuracy requirements of the "Code for Urban Engineering Geophysical Exploration CJJ7-2016".

[0127] As shown in Table 3, the total accuracy and the accuracy distribution of the high-precision gravity measurement data correction calculation during the test are as follows:

[0128] During the test, the urban gravity high-precision measurement data processing software developed based on the core algorithm of the present invention was used to process 5 measured gravity data profiles collected. The processing contents included: (I) preprocessing of the data collected by the instrument; (II) conventional correction processing such as solid tide correction, normal field correction, and Bouguer correction; (III) terrain correction processing for near, middle, and far areas according to natural geomorphic features; (IV) ground object correction processing for near, middle, and far areas due to the influence of surface buildings; (V) gravity field correction processing for known special underground buildings; (VI) abnormal separation processing of regional field and residual field; (VII) two-dimensional and three-dimensional data inversion processing.

[0129] Among them, when performing ground object correction processing for urban building groups, a three-dimensional spatial attribute data model of the building groups in the 25 km² method test area was constructed for calculation. Among them, when performing terrain correction processing for natural geomorphic features, for the middle area terrain correction calculation (20 m - 2 km), 1:10,000 scale digital elevation data (DEM) purchased from the local surveying and mapping management department was used, and for the far area terrain correction calculation (2 km - 166.7 km), ASTER GDEM V2 data with a resolution of 1″×1″ was used for calculation.

[0130] According to the relevant geological background information of the test area, the Quaternary strata exposed on the surface are gravel, pebble layers of the floodplain facies, in unconformable contact with the underlying strata, and the stratum thickness is 20 - 40 m; the underlying strata are the Guankou Formation K2g of the upper Cretaceous series, with lithology of reddish-brown mudstone and siltstone mudstone, containing mirabilite and gypsum, and the stratum thickness > 149 m. The statistical table of the lithology parameters of the main strata is shown in Table 4 and Figure 13 as follows.

[0131] At the same time, shallow seismic, equivalent reverse flux transient electromagnetic sounding, and microtremor exploration work were deployed for the concealed fault of Sumadou in the test area. Through the collected data such as engineering geological boreholes, core physical properties, seismic, and electrical methods, the physical property parameters such as resistivity and density of the main strata in the test area were obtained. By comparing the actual detection effects of gravity, electrical method, and seismic exploration methods on the same exploration line profile, the actual effects of each method were evaluated from indicators such as "exploration cost, work efficiency, effective depth, detection accuracy, and minimum resolution". As Figure 14 shown is the high-precision gravity exploration result of the Sumadou fault in the test area. As Figure 14 shown in the density inversion result of the high-precision gravity profile of line L2, there is an obvious low-density abnormal body at the position of the Sumadou fault, and high-density bodies are on both sides of the fault position.

[0132] The high-resolution reflection wave seismic exploration method was used to carry out shallow seismic detection work on the concealed Sumadou fault spreading in Tianfu New Area of Chengdu. A receiving point spacing of 2 m, a shot point spacing of 4 m, 512 receiving channels, 128 coverage times, and a middle-shot symmetric receiving observation system were adopted. AsFigure 15 As shown in the figure, the Sumadou Fault is mainly characterized by a pop-up structure. A pop-up structure or a 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 Figure 3 shows the spatial distribution position of the Sumadou Fault. Two or more faults with opposite dips are arranged in a "V" shape or an "8" shape. The "V" shape combination shows that one or more faults are thrusting against a large thrust fault and are buried beneath it, forming a pop-up structure; the "8" 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 Sumadou 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 are manifested as 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.

[0133] Two equivalent anti-flux transient electromagnetic profiles were completed near the above-mentioned shallow seismic profile location, controlling the Sumadou 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 displacement of this fault is between 10 and 490 m, and the displacement gradually decreases from southwest to northeast, and the attitude is steeper 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 overturned. 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 basement-cover interface and the fault fracture zone within 200 m, and at the same time, it also has a certain ability to identify multiple mined-out areas within 100 m, with good resolution.

[0134] As Figure 18 and Figure 19 shown, two microtremor profile explorations were completed near the above-mentioned equivalent anti-flux transient electromagnetic sounding profile location. In the case of strong interference, microtremor measurement can use the passive source surface wave background noise for detection and 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.

[0135] In summary, through the comparison of various exploration methods such as microgravity, reflection seismology, transient electromagnetic method, and microtremor method at positions adjacent to the Sumadou fault, the correctness, reliability, and effectiveness of the application of the present invention in urban microgravity exploration work are demonstrated. At the same time, it is also demonstrated that the microgravity anomaly information obtained by using the present invention can be used to identify the spatial distribution of hidden faults underground in the city.

[0136] The basic concepts have been described above. Obviously, for those skilled in the art, the above detailed disclosure is only an example and does not constitute a limitation to this specification. Although not explicitly stated here, those skilled in the art may make various modifications, improvements, and corrections to this specification. Such modifications, improvements, and corrections are suggested in this specification, so such modifications, improvements, and corrections still fall within the spirit and scope of the exemplary embodiments of this specification.

[0137] At the same time, this specification uses specific terms to describe the embodiments of this specification. For example, "some embodiments" means a certain feature, structure, or characteristic related to at least one embodiment of this specification. Therefore, it should be emphasized and noted that "some embodiments" mentioned twice or more at different positions in this specification are not necessarily referring to the same embodiment. In addition, certain features, structures, or characteristics in one or more embodiments of this specification can be appropriately combined.

[0138] In addition, the order of the processing elements and sequences, the use of numbers and letters, or the use of other names described in this specification are not used to limit the order of the processes and methods of this specification. Although some currently considered useful invention embodiments are discussed through various examples in the above disclosure, it should be understood that such details only serve the purpose of illustration and are not limited to the disclosed embodiments. For example, although the system components described above can be implemented by hardware devices, they can also be implemented only through software solutions, such as installing the described system on existing servers or mobile devices.

[0139] Similarly, it should be noted that in order to simplify the expression of the disclosure of this specification and thus help the understanding of one or more invention embodiments, in the description of the embodiments of this specification above, sometimes multiple features are merged into one embodiment, drawing, or description thereof. In fact, the features of the embodiments are less than all the features of the individual embodiments disclosed above.

[0140] In some embodiments, numbers are used to describe components and the quantity of attributes. It should be understood that such numbers used in the description of embodiments are, in some examples, modified by the modifiers "about", "approximately", or "substantially". Unless otherwise specified, "about", "approximately", or "substantially" indicate that the stated number allows for a variation of ±20%. Accordingly, in some embodiments, the numerical parameters used in the specification are approximate values, and such approximate values may vary according to the characteristics required by individual embodiments. In some embodiments, the numerical parameters should consider the specified significant digits and adopt the method of retaining the general number of digits. Although the numerical ranges and parameters used in some embodiments of this specification to confirm the breadth of their scope are approximate values, in specific embodiments, such numerical settings are made as precise as possible within the feasible range.

[0141] For each patent, patent application, patent application publication, and other materials cited in this specification, such as articles, books, specifications, publications, documents, etc., their entire contents are hereby incorporated into this specification by reference. This excludes the application history files that are inconsistent with or conflict with the content of this specification, as well as the files that limit the broadest scope of the claims of this specification (currently or subsequently appended to this specification). It should be noted that if there are inconsistencies or conflicts between the descriptions, definitions, and / or uses of terms in the supplementary materials of this specification and the content described in this specification, the descriptions, definitions, and / or uses of terms in this specification shall prevail.

[0142] Finally, it should be understood that the embodiments described in this specification are only used to illustrate the principles of the embodiments of this specification. Other variations may also fall within the scope of this specification. Therefore, by way of example and not limitation, alternative configurations of the embodiments of this specification may be considered to be consistent with the teachings of this specification. Accordingly, the embodiments of this specification are not limited to the embodiments explicitly introduced and described in this specification.

Claims

1. A gravity field correction method, characterized in that: The method comprises: Acquire a three-dimensional data model of the building complex in the target area; the three-dimensional data model is a model based on a tetrahedral grid structure; The three-dimensional data model is processed based on a three-dimensional forward numerical simulation method to obtain a gravity correction value of the gravity measurement point for buildings surrounding the gravity measurement point; Based on the gravity correction value, the actual gravity field value at the gravity measurement point is obtained.

2. The method according to claim 1, characterized in that: The three-dimensional data model of the building complex in the target area is obtained by: Acquiring satellite remote sensing data of the target area; Based on the satellite remote sensing data, obtaining the spatial outline and elevation data of the building complex; The three-dimensional data model is constructed based on the spatial contour and elevation data.

3. The method according to claim 2, characterized in that: The acquiring of the spatial outline and elevation data of the building complex based on the satellite remote sensing data comprises: Preprocessing the satellite remote sensing data to obtain processed data, wherein the preprocessing includes at least one of image fusion, custom coordinate system, orthorectification, 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; Acquire data features corresponding to the segmented data, and based on the merging scale and according to the band features corresponding to the data features, perform band merging on the segmented data to obtain merged data; Selecting a data sample from the merged data; Based on the sample statistical results of the data samples, the spatial outline and elevation data of the building complex are obtained.

4. The method according to claim 3, characterized in that: The constructing of the three-dimensional data model based on the spatial profile and elevation data comprises: Based on the spatial contour and elevation data, construct an initial data model; Based on the graphic attributes corresponding to each building, the buildings included in the initial data model are clustered to obtain a plurality of clusters; wherein the buildings under a cluster correspond to a building type, and the graphic attributes include at least one of surface area and image features; For each cluster, obtain the density value of the buildings contained in the cluster, and use the density value as the density attribute of the buildings in the cluster; In the initial model, density values ​​are assigned to the buildings under each cluster to obtain the three-dimensional data model with density attributes.

5. The method according to claim 4, characterized in that: For each cluster, obtaining the density value of the buildings contained therein includes: For each cluster, one of the buildings is taken as a building sample; For the building sample, a plurality of candidate densities are determined through density hypothesis; Perform three-dimensional forward fitting based on the candidate densities to obtain correction value simulation curves corresponding to multiple candidate densities; Conduct field measurements for each building sample to obtain a calibration value measurement curve corresponding to each building sample; The candidate density corresponding to the correction value simulation curve having the greatest similarity to the correction value measured curve is used as the density value of the building sample; For the cluster where the building sample is located, the density values ​​of all its buildings are taken as the density value.

6. The method according to claim 5, characterized in that The three-dimensional data model is processed based on the three-dimensional forward numerical simulation method to obtain the gravity correction value of the gravity measurement point for the buildings around the gravity measurement point, including: Based on a preset algorithm, obtaining a gravity correction value of each tetrahedral unit in the three-dimensional data model at the gravity measurement point; The gravity correction values ​​corresponding to all the tetrahedral units within a target range from the gravity measurement point are accumulated to obtain the gravity correction value of the gravity measurement point for the buildings surrounding the gravity measurement point.

7. The method according to claim 6, characterized in that The step of obtaining the gravity correction value of each tetrahedral unit in the three-dimensional data model at the gravity measurement point based on a preset algorithm includes: For any tetrahedral unit ABCD, using formula (1): Determining the gravity correction value of the tetrahedral element at the gravity measurement point; Wherein, the coordinates of the four vertices of the tetrahedral unit ABCD are expressed 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 unit; (x, y, z) represents the coordinates of the center of gravity of the tetrahedral unit ABCD, and the center of gravity of the tetrahedral unit ABCD can be determined based on the coordinates of the four vertices; Using Gauss's formula Convert the three-dimensional volume integral of formula (1) 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 unit ABCD are all solid constants, combined with the Gaussian formula, formula (1) is transformed into formula (2): Where j represents the jth triangular face of the tetrahedral unit ABCD, γ j Represents the normal vectors of each triangular face of the tetrahedral unit ABCD; the normal vectors of each triangular face are obtained from the coordinates of each vertex of the triangle.

8. A gravity field correction system, characterized in that: It includes data acquisition module, data processing module and data calculation module; The data acquisition module is configured to obtain a three-dimensional data model of the building complex in the target area; the three-dimensional data model is a model based on a tetrahedral 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 the gravity measurement point for buildings surrounding the gravity measurement point; The data calculation module is configured to obtain an actual gravity field value at a gravity measurement point based on the gravity correction value.

9. A gravity field correction device, comprising a processor, characterized in that: The processor is used to execute the gravity field correction method as described in any one of claims 1-7.

10. A computer-readable storage medium storing computer instructions, characterized in that: After the computer reads the computer instructions in the storage medium, the computer executes the gravity field correction method as described in any one of claims 1-7.

Images