A gravity terrain correction accuracy constraint method, system, storage medium and product

Abstract

The application provides a gravity terrain correction precision constraint method, system, storage medium and product, relates to the field of electric digital data processing, and the method comprises the following steps: generating an initial non-uniform calculation grid according to the local terrain complexity of a digital elevation model; performing geometric fitting on the real terrain data in each grid cell in the initial calculation grid to obtain a first physical model and a second physical model; calculating the gravity effect of each grid cell on an observation point based on the two physical models respectively to obtain two gravity terrain correction influence values and an approximate degree residual corresponding to the difference between the two gravity terrain correction influence values; when the approximate degree residual exceeds a preset precision threshold, refining and dividing the grid cell, and refining and dividing each new calculation grid generated again until the new approximate degree residual does not exceed the preset precision threshold. The method can improve the overall calculation efficiency on the premise of ensuring the global calculation precision, thereby solving the calculation redundancy problem of the traditional uniform grid method.

Description

Technical Field

[0001] This application relates to the field of electronic digital data processing, and in particular to a gravity terrain correction accuracy constraint method, system, storage medium and product. Background Technology

[0002] In geophysical exploration, geodesy, and national basic mapping, high-precision gravity data is fundamental for geological structure interpretation, resource exploration, and geoid refinement. However, the gravity values ​​directly observed on the ground are the result of a combination of factors, including the distribution of materials within the Earth, the Earth's rotation, and the topographic mass between the observation point and the reference ellipsoid. Among these, the influence of topographic mass on gravity observations, known as gravity-topographic correction, is a crucial correction that must be accurately subtracted. Therefore, performing high-precision gravity-topographic correction is essential for extracting effective gravity signals that reflect anomalies in subsurface density structure.

[0003] In related technologies, a grid-based calculation method based on a Digital Elevation Model (DEM) is commonly used. Specifically, this method divides the DEM data covering the entire survey area into regular and uniform rectangular grids. For each grid cell, the complex terrain within it is simplified into a regular geometric shape; for example, it is typically approximated as a flat-topped rectangular prism with a defined height. The height of this prism is generally determined by taking the average of all DEM elevation points within the grid cell or the elevation of its center point. Then, the gravity effect generated by this prism at a specified gravity observation point is calculated, and the gravity effect of all grid cells across the entire survey area is integrated and summed to obtain the total gravity terrain correction value. To ensure calculation accuracy, this method tends to use the finest possible, uniform resolution grid for subdivision throughout the entire survey area, aiming to more accurately approximate the actual terrain.

[0004] However, to ensure the overall accuracy of the final calculation results, related technologies generally select grid resolution based on accurately describing the most complex and undulating areas of the survey area (such as steep ridges or deep canyons). This is because if the grid is too coarse in these critical areas, it will lead to huge approximation errors in the terrain model, thus seriously affecting the reliability of the entire correction result. However, since this technical solution uses a uniform grid, the high-resolution grid that meets the accuracy requirements of locally complex areas is applied indiscriminately to the entire survey area. This results in the same high-intensity calculations with the same density as in complex areas in vast areas with relatively flat and gently undulating terrain (such as plains and basins). In these flat areas, the same accuracy requirements could be achieved using a much sparser grid. Therefore, this implementation method using a uniform and fine grid results in a large amount of computational resources being wasted on unnecessary redundant calculations, thus making the overall efficiency of the calculation process low. Summary of the Invention

[0005] This application provides a gravity terrain correction accuracy constraint method, system, storage medium, and product to address the problem of low overall computational efficiency caused by the implementation method using uniform fine mesh in related technologies, which needs to meet the accuracy requirements of local complex regions.

[0006] In a first aspect, this application provides a gravity terrain correction accuracy constraint method, applied to a gravity terrain correction accuracy constraint system, the method comprising:

[0007] A non-uniform initial computational grid is generated based on the local terrain complexity of the digital elevation model, wherein the digital elevation model is obtained by fusing terrain elevation data of at least two different resolutions, and the size of the initial computational grid is inversely correlated with the corresponding local terrain complexity.

[0008] Geometric fitting is performed on the real terrain data in each grid cell of the initial computational grid to obtain a parameterized first physical model and a second physical model. The first physical model and the second physical model have different complexities in the same grid cell.

[0009] The gravity effect of each grid cell on the observation point is calculated based on the first physical model and the second physical model respectively, to obtain the first gravity terrain correction influence value and the second gravity terrain correction influence value;

[0010] The approximation residual of the physical model is calculated based on the difference between the first gravity terrain correction influence value and the second gravity terrain correction influence value;

[0011] When the approximation residual exceeds a preset accuracy threshold, the grid cells are refined and divided to generate multiple new computational grids;

[0012] The steps of geometric fitting, gravity effect calculation, approximation residual and refinement segmentation are repeated for each new computational grid until the approximation residual corresponding to all grid cells does not exceed the preset accuracy threshold.

[0013] By adopting the above technical solution, the system first generates a non-uniform initial computational grid based on the local terrain complexity, making the grid size inversely correlated with the terrain complexity, thus achieving preliminary optimized allocation of computational resources. Then, it performs dual modeling by constructing first and second physical models with different complexities and calculates the difference in gravity effects between the two to obtain the approximation residual. This difference calculation mechanism can more accurately quantify the model approximation error of each grid cell. When the residual exceeds a preset accuracy threshold, the system adaptively refines the grid and iterates until the accuracy requirements are met. This adaptive refinement mechanism ensures that a denser grid is used in complex terrain areas to guarantee accuracy, while a sparser grid is used in simple terrain areas to save computational resources. Ultimately, it can improve overall computational efficiency while ensuring global computational accuracy, thereby addressing the computational redundancy problem of traditional uniform grid methods.

[0014] In some embodiments, the step of refining the grid cells to generate multiple new computational grids specifically includes:

[0015] If the current subdivision level corresponding to the grid cell is less than the preset maximum subdivision level, the grid cell is cut into four sub-grid cells with the same geometric shape and equal area. The maximum subdivision level is set according to the accuracy of the digital elevation model.

[0016] The subgrid cell is identified as the new computing grid, and the current subdivision level corresponding to the new computing grid is updated;

[0017] If the current subdivision level corresponding to the grid cell is greater than or equal to the preset maximum subdivision level, the grid cell is marked as a non-convergent cell, and further subdivision of the non-convergent cell is terminated.

[0018] By adopting the above technical solution, the system introduces a hierarchical control mechanism during mesh refinement. When the subdivision level has not reached its maximum value, the mesh is uniformly divided into four sub-meshes. This quadtree partitioning method ensures the regularity and controllability of mesh refinement. Simultaneously, setting the maximum subdivision level as a constraint avoids infinite subdivision in extremely complex regions, preventing excessive consumption of computational resources. For cells that still do not meet the accuracy requirements even after reaching the maximum level, the system marks them as non-convergent cells and terminates the partitioning. This marking mechanism avoids invalid computational loops and provides clear location information for subsequent accuracy analysis and optimization.

[0019] In some embodiments, the step of geometrically fitting the real terrain data within each grid cell of the initial computational grid to obtain a parameterized first physical model and a second physical model specifically includes:

[0020] Obtain sampling point data of the digital elevation model within the area covered by the grid cell;

[0021] Based on the sampling point data, the terrain within the grid cell is simplified into a flat-topped prism with uniform height to obtain the first physical model. The height value of the flat-topped prism is the arithmetic mean of the elevation values ​​of the sampling points.

[0022] The sampling point data is fitted with a high-order surface to generate a second physical model constructed from surface prisms.

[0023] By adopting the above technical solution, the system first acquires sampling point data of the digital elevation model within the grid cells, and then constructs two physical models with different complexities: the first physical model simplifies the terrain into flat-topped prisms and uses the arithmetic mean as the height, achieving a fast approximation; while the second physical model, which generates curved prisms through high-order surface fitting, more accurately describes the terrain undulation characteristics. This dual-model mechanism allows the system to accurately assess the magnitude of the error caused by model simplification by comparing the differences between the two. The first model has high computational efficiency but lower accuracy, while the second model has high accuracy but is computationally complex. The combination of the two provides a basis for judgment in subsequent adaptive accuracy control, achieving an effective balance between computational efficiency and accuracy assessment.

[0024] In some embodiments, the step of calculating the gravity effect of each grid cell on the observation point based on the first physical model and the second physical model respectively, and obtaining the first gravity terrain correction influence value and the second gravity terrain correction influence value, specifically includes:

[0025] Based on the distance between the grid cells and the observation points, the computational domain is divided into the near zone, the middle zone, the far zone 1, and the far zone 2.

[0026] For the grid cells in the near zone, middle zone, and far zone, a high-precision integral algorithm is used to calculate the gravity terrain correction influence value based on the first physical model and the second physical model, respectively. The gravity terrain correction influence value includes the first gravity terrain correction influence value and the second gravity terrain correction influence value.

[0027] For the far second zone grid cell, a fast approximation algorithm is used to calculate the gravity terrain correction value based on the first physical model and the second physical model respectively.

[0028] By adopting the above technical solution, the system intelligently divides the computational region into near zone, middle zone, far zone 1, and far zone 2 based on the distance between the grid cells and the observation points, realizing a differentiated processing strategy for each zone. A high-precision integration algorithm is used for the near, middle, and far zones 1, ensuring computational accuracy in the near-field region where the gravitational field has a significant impact. For the far zone 2, a fast approximation algorithm is used, leveraging the rapid decay of the far-field effect to significantly reduce computational complexity while maintaining accuracy. This zoned computational strategy fully considers the spatial decay characteristics of the gravitational field and optimizes the efficiency of global computation through the appropriate selection of algorithms, reducing computation time without sacrificing overall accuracy. It is particularly suitable for large-scale terrain correction calculations.

[0029] In some embodiments, after the step of calculating the approximation residual of the physical model based on the difference between the first gravity terrain correction influence value and the second gravity terrain correction influence value, the method further includes:

[0030] When the approximation residual does not exceed the preset accuracy threshold and meets the preset refinement conditions, the current grid cell is temporarily divided into four sub-grid cells to form a set of sub-grid cells.

[0031] The gravity effect value of the observation point corresponding to each sub-grid cell is calculated based on the first physical model corresponding to each sub-grid cell. The calculation method of the gravity effect value is the same as the calculation method of the first gravity terrain correction influence value.

[0032] When the difference between the sum of the gravity effect values ​​and the first gravity terrain correction effect value is greater than a preset verification threshold, the sub-grid cell corresponding to the target gravity effect value whose difference from the mean of the gravity effect values ​​exceeds a preset difference threshold is marked as a cell to be refined.

[0033] When the difference between the sum of the gravity effect values ​​and the first gravity terrain correction effect value is less than or equal to a preset verification threshold, further segmentation of the current grid cell is terminated.

[0034] By adopting the above technical solution, the system introduces a secondary verification mechanism when the approximation residual meets the requirements but still meets the refinement conditions: by comparing the differences in terrain correction values ​​between the current grid and its four sub-grids, the system accurately identifies areas that need further refinement. This verification mechanism avoids misjudgments that may be caused by relying solely on the approximation residual, especially when the local terrain changes are drastic but the overall residual is not large, and can accurately capture the key locations that need refinement. The segmentation terminates when the differences in correction values ​​of all sub-grids are within a threshold, preventing unnecessary over-refinement. This dual judgment mechanism improves the accuracy and efficiency of adaptive grid refinement, ensuring that computational resources are accurately allocated to areas that truly require high-precision computation.

[0035] In some embodiments, the step of meeting the preset refinement conditions specifically includes:

[0036] Calculate the roughness index of the terrain elevation data within the grid cell, wherein the roughness index is determined based on the standard deviation of the terrain elevation data;

[0037] When the roughness index is less than the preset roughness threshold, it is determined that the mesh element does not meet the preset refinement condition;

[0038] When the roughness index is not less than the roughness threshold, the mesh cell is determined to meet the preset refinement conditions.

[0039] By adopting the above technical solution, the system introduces a roughness index as the criterion for determining refinement conditions. This index, calculated based on the standard deviation of terrain elevation data, can quantitatively reflect the degree of terrain undulation within a grid cell. When the roughness is below a threshold, refinement is deemed unnecessary. This mechanism can quickly identify flat terrain areas and avoid unnecessary calculations. The introduction of the roughness index allows the system to pre-screen grids that do not require refinement through simple statistical calculations before performing complex dual-model calculations, significantly reducing the number of grids requiring full accuracy evaluation. This pre-screening mechanism significantly improves overall computational efficiency, especially in survey areas containing large flat regions.

[0040] In some embodiments, after the step of repeatedly performing the steps of geometric fitting, gravity effect calculation, approximation residual, and refinement segmentation on each of the new computational grids until the approximation residuals corresponding to all grid cells do not exceed the preset accuracy threshold, the method further includes:

[0041] The gravity terrain correction influence values ​​corresponding to all grid cells that meet the preset accuracy threshold are summarized to generate a global terrain correction influence value matrix. The correction influence value matrix is ​​the final gravity terrain correction result covering the entire calculation area.

[0042] A standardized terrain correction data file is generated based on the correction influence value matrix. The standardized terrain correction data file is an output file that conforms to the data format specifications of GIS or geophysical software.

[0043] An accuracy certificate report is generated based on the correction influence value matrix and the final cell grid distribution. The accuracy certificate report includes a non-uniform grid distribution map showing the differences in grid density in different regions, an accuracy heatmap representing the magnitude of calculation errors in different regions, and key performance indicators.

[0044] By adopting the above technical solution, the system aggregates the terrain correction values ​​of all grid cells that meet the accuracy requirements to generate a global correction matrix, effectively integrating the distributed calculation results. The generated standardized terrain correction data file conforms to the format specifications of GIS and geophysical software, ensuring that the calculation results can be seamlessly integrated into subsequent data processing and analysis workflows. The non-uniform grid distribution map included in the accuracy certificate report visually demonstrates the differences in calculation density in different regions, the accuracy heatmap quantitatively shows the distribution of calculation errors in each region, and the key performance indicators provide a quantitative assessment of the overall calculation quality. This complete output system not only provides high-quality calculation results but also provides detailed accuracy verification and quality assurance information, greatly enhancing the credibility and practicality of gravity terrain correction results.

[0045] Secondly, this application provides a gravity terrain correction accuracy constraint system, the system comprising: one or more processors and a memory;

[0046] The memory is coupled to the one or more processors. The memory is used to store computer program code, which includes computer instructions. The one or more processors call the computer instructions so that the system can implement the gravity terrain correction accuracy constraint method provided in the above embodiments, which will not be described in detail here.

[0047] Thirdly, this application provides a computer-readable storage medium including instructions that, when executed on a gravity terrain correction accuracy constraint system, enable the system to implement a gravity terrain correction accuracy constraint method provided in the above embodiments, which will not be elaborated here.

[0048] Fourthly, this application provides a computer program product that, when run on a gravity terrain correction accuracy constraint system, enables the system to implement a gravity terrain correction accuracy constraint method provided in the above embodiments, which will not be elaborated here.

[0049] One or more technical solutions provided in the embodiments of this application have at least the following technical effects or advantages:

[0050] 1. By establishing an inverse correlation between grid size and local terrain complexity, and using the approximation residual calculated by the dual-physical model as the criterion for grid refinement, on-demand allocation of computing resources is achieved. This can significantly improve overall computing efficiency while ensuring global accuracy, effectively solving the computational redundancy problem of traditional uniform grid methods.

[0051] 2. The computational region is divided into different partitions based on distance and a differentiated algorithm is adopted, which makes full use of the physical property that the gravitational field decays with distance. At the same time, the introduction of mechanisms such as secondary verification and hierarchical control enables the algorithm to achieve the optimal balance between accuracy requirements and computational efficiency under different spatial scales and complex scenarios.

[0052] 3. From dual-model error assessment to secondary verification, and finally to the accuracy certificate report containing visualization charts and key indicators, this application constructs a complete quality control closed loop from the calculation process to the presentation of results, which greatly improves the credibility, reliability and traceability of gravity terrain correction results. Attached Figure Description

[0053] Figure 1 This is a flowchart illustrating a gravity terrain correction accuracy constraint method in an embodiment of this application.

[0054] Figure 2 This is another flowchart illustrating a gravity terrain correction accuracy constraint method in the embodiments of this application;

[0055] Figure 3 This is a schematic diagram of a physical device structure of a gravity terrain correction accuracy constraint system in the embodiments of this application. Detailed Implementation

[0056] The terminology used in the following embodiments of this application is for the purpose of describing particular embodiments only and is not intended to be limiting of this application. As used in the specification and appended claims of this application, the singular expressions “a,” “an,” “the,” “the,” “the,” and “this” are intended to include the plural expressions as well, unless the context clearly indicates otherwise. It should also be understood that the term “and / or” as used in this application refers to any or all possible combinations including one or more of the listed items.

[0057] Hereinafter, the terms "first" and "second" are used for descriptive purposes only and should not be construed as implying or suggesting relative importance or implicitly indicating the number of indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature, and in the description of the embodiments of this application, unless otherwise stated, "multiple" means two or more.

[0058] For ease of understanding, the method provided in this implementation is described in process below. Please refer to [link / reference]. Figure 1 This is a flowchart illustrating a gravity terrain correction accuracy constraint method in an embodiment of this application.

[0059] S101. Generate a non-uniform initial computational grid based on the local terrain complexity of the digital elevation model.

[0060] Among them, a digital elevation model (DEM) refers to a model that can digitally represent the elevation information of the Earth's surface, obtained by fusing topographic elevation data from at least two different resolutions. The data sources for fusion include, but are not limited to, the RGIS2006 elevation database (using a 1km*1km planar grid for the range of 2-20km and a 5′5′ planar grid for the range of 20-166.7km) and the STRM high-resolution global elevation database (using a 30m30m planar grid, with an elevation accuracy of up to 10 meters); local topographic complexity refers to the characteristic index used to represent the degree of undulation in a local area of ​​the terrain, usually quantified by slope, curvature, topographic roughness (determined based on the standard deviation of the topographic elevation data) or topographic maturity index (such as RMS slope); non-uniform initial computation grid refers to the initial grid set in which the size of the grid cell is inversely correlated with the corresponding local topographic complexity, that is, the grid density is high (small cells) in complex terrain areas and low (large cells) in flat terrain areas; topographic elevation data refers to the basic data used to construct the digital elevation model, which contains the elevation information of various points on the Earth's surface at different resolutions.

[0061] Specifically, the system first calls the preprocessed multi-source fused digital elevation model—the system can first convert the format of the raw data from the RGIS2006 elevation database and the STRM high-resolution elevation database (supporting GeoTIFF, LAS, etc.), unify the coordinate system, and filter noise, and then fuse data of different resolutions through data interpolation technology to generate a high-precision DEM covering the entire calculation area. The system then performs local terrain complexity analysis on this DEM:

[0062] Optionally, the DEM is traversed using a sliding window method to calculate the slope, curvature, or RMS slope (terrain maturity index) within each window. Simultaneously, a terrain roughness index is calculated based on the standard deviation of elevation data. These parameters are then combined to quantify the terrain complexity of each local area. Finally, based on the inverse correlation rule of "higher complexity, smaller grid," the system assigns corresponding grid parameters to areas of different complexities—for example, in rugged areas with steep slopes and high curvature, a small-size grid with a resolution close to 30m from the STRM database is used; in plain areas with gentle slopes and low roughness, a large-size grid with a resolution close to 1km from the RGIS2006 database is used, ultimately generating a non-uniform initial computational grid covering the entire region.

[0063] In one specific embodiment, the system imports preprocessed RGIS2006 and STRM fused DEM data and delineates the boundary of the calculation area; it uses a 500m500m sliding window to traverse the DEM, calculates the average slope and elevation standard deviation within each window, and sums the two with weights of 0.6 and 0.4 to obtain the quantized value of complexity; it presets a complexity threshold range (e.g., [0, 0.3] corresponds to a 1km1km grid, [0.3, 0.7] corresponds to a 100m100m grid, and [0.7, 1.0] corresponds to a 30m30m grid), matches the grid size according to the quantized value of each window, and generates a non-uniform initial calculation grid.

[0064] It should be noted that adopting a multi-source fusion digital elevation model is the preferred way to improve data quality in this field, but the core of the accuracy constraint method proposed in this application lies in the adaptive grid refinement and error evaluation mechanism, so it is also applicable to digital elevation models generated from a single data source.

[0065] S102. Perform geometric fitting on the real terrain data in each grid cell of the initial computational grid to obtain the parameterized first physical model and second physical model.

[0066] Among them, grid cell refers to the smallest indivisible independent computational unit in a non-uniform initial computational grid, and its size dynamically changes with the local terrain complexity; real terrain data refers to the elevation data of sampling points extracted from the fused DEM within the area covered by the grid cell (the sampling point density matches the DEM resolution, such as 1 / 30m sampling point for a 30m grid); geometric fitting refers to the process of transforming discrete real terrain data into continuous geometric shapes through mathematical methods, which is used to simplify terrain features to adapt to gravity effect calculations; the parameterized first physical model refers to the model that simplifies the grid cell terrain into a uniformly high flat-topped prism, whose core parameters are the length, width (consistent with the grid cell), and height (arithmetic mean of sampling point elevations) of the prism; the parameterized second physical model refers to the curved prism model obtained by fitting the grid cell terrain with a higher-order curved surface, whose core parameters are the coefficients of the surface equation (such as the a, b, and c coefficients of a quadratic surface) and the base area of ​​the prism (consistent with the grid cell).

[0067] Specifically, the system processes each grid cell in the initial calculation grid one by one: first, it extracts the real terrain data within the area covered by the grid cell from the fused DEM—and determines the number of sampling points according to the grid cell size (e.g., extract 9 uniformly distributed sampling points for a 30m30m grid, and extract 25 uniformly distributed sampling points for a 1km1km grid) to ensure that the sampling points can reflect the overall terrain within the cell. The system then performs geometric fitting in two steps: First, a first physical model is constructed. The arithmetic mean of the elevations of all sampling points is calculated, and this mean is used as the height. A regular flat-topped prism is constructed using the length and width of the grid cells as the base dimensions. The model is parameterized using the length, width, and height of the prism. This model prioritizes computational efficiency and ignores subtle terrain undulations. Second, a second physical model is constructed. A higher-order surface fitting algorithm (such as quadratic polynomial surface fitting) is used to fit the sampling point data. The coefficients of the surface equation (such as z=ax²+by²+cxy+dx+ey+f) are solved using the least squares method. This surface is used as the top face of the prism, and the length and width of the grid cells are used as the base dimensions to construct a surface prism. The model is parameterized using the surface equation coefficients and the base dimensions. This model prioritizes terrain realism and preserves subtle undulations. Ultimately, each grid cell corresponds to two parameterized physical models.

[0068] S103. Calculate the gravity effect of each grid cell on the observation point based on the first physical model and the second physical model respectively, and obtain the first gravity terrain correction influence value and the second gravity terrain correction influence value.

[0069] Among them, gravity effect refers to the additional influence of the terrain quality represented by the grid cell (calculated based on the physical model) on the gravity observation value at the observation point, which needs to be calculated using the gravity field theory formula; observation point refers to the specific location point for gravity data acquisition, whose coordinates (latitude, longitude, elevation) must be consistent with the DEM coordinate system, commonly found at mineral exploration points, geological disaster monitoring points, etc.; the first gravity terrain correction influence value refers to the correction value of the grid cell on the gravity observation value of the observation point calculated based on the first physical model (flat-topped prism), reflecting the gravity influence of simplified terrain; the second gravity terrain correction influence value refers to the correction value of the grid cell on the gravity observation value of the observation point calculated based on the second physical model (curved prism); high-precision integration algorithm refers to the accurate method (such as prism volume integration) used to calculate the gravity effect of near-to-mid-distance grid cells, suitable for areas with significant gravity influence; fast approximation algorithm refers to the simplified method (such as FFT fast Fourier transform algorithm) used to calculate the gravity effect of far-distance grid cells, improving computational efficiency while ensuring accuracy.

[0070] Specifically, the system determines the coordinate information (latitude, longitude, and elevation) of the observation point and unifies it with the coordinate system of the non-uniform initial calculation grid (e.g., converting it to the WGS84 coordinate system). Then, based on the straight-line distance between the grid cell and the observation point, the system divides all grid cells into four calculation regions: near zone (0-50m), middle zone (50m ≤ distance < 2km), far zone one (2km ≤ distance < 20km), and far zone two (20km ≤ distance < 166.7km). This partitioning is based on the characteristic of gravity field attenuation with distance and in accordance with the "Regional Gravity Survey Specification." The near and middle zones have a significant impact on the gravity of the observation point and require high-precision calculations, while the far zone two has a weak impact and can be simplified in calculations.

[0071] Next, the system calculates the gravity effect by region and model: For the grid cells in the near, middle, and far regions, the high-precision integration algorithm in the core computing engine (such as the prism volume integral method with integrated GPU acceleration) is called, and the length, width, and height parameters of the first physical model (flat-topped prism) and the surface equation coefficients of the second physical model (curved prism) are substituted to calculate the gravity influence value of each model on the observation point; For the grid cells in the far region, the fast approximation algorithm (such as the FFT algorithm) is called, and the calculation is simplified based on the equivalent mass and distance parameters of the model to obtain the gravity influence values ​​of the two models respectively.

[0072] In one specific embodiment, the system acquires the coordinates (X0, Y0, Z0) of the observation point and calculates the straight-line distance d between the center point of each grid cell and the observation point. For grid cells with d < 2km, a prism volume integral algorithm is used: the length L, width W, and height H of the first physical model (flat-topped prism) are substituted into the formula to calculate the gravity influence value; the surface equation z = f(x, y) of the second physical model (curved prism) is substituted into the formula to calculate the gravity influence value through double integral. For grid cells with d ≥ 2km, a fast FFT approximation algorithm is used: the first and second physical models are equivalent to point masses (mass = base area × average height × ρ), and the gravity influence value is calculated through the formula.

[0073] In another embodiment, the system divides the computational region into a near zone (0m≤d<20m), a middle zone (20m≤d<500m), and a far zone (500m≤d<2km) based on distance, omitting the far zone division. For the near zone grid cells, a high-precision GPU-accelerated integration algorithm is used, while simultaneously calculating the gravity influence values ​​of the first and second physical models (parallel computing improves efficiency). For the middle zone grid cells, a simplified integration algorithm is used (ignoring minor undulations such as the tilt angle of the model top surface being <5°). For the far zone grid cells, a spherical approximation algorithm is used (incorporating the curvature of the Earth into the calculation). Error corrections are applied to the calculation results for each region (e.g., correcting for terrain density inhomogeneity in the near zone and correcting for spherical approximation errors in the far zone).

[0074] It is understandable that other methods can be used to calculate the influence value of gravity terrain correction, such as combining Monte Carlo simulation to correct calculation errors, which is not limited here.

[0075] S104. The approximation residual of the physical model is calculated based on the difference between the influence value of the first gravity terrain correction and the influence value of the second gravity terrain correction.

[0076] Among them, the approximation residual of the physical model refers to the index used to quantify the degree of approximation error of the first physical model to the real terrain. The smaller the residual, the higher the fit between the first physical model and the real terrain, and the closer the calculation accuracy is to the second physical model. Conversely, the larger the residual, the larger the approximation error of the first physical model, which may require refinement of the grid cells.

[0077] Specifically, the system retrieves the first and second gravity terrain correction influence values ​​calculated for each grid cell. Since the first physical model is a simplified fit to the terrain, while the second physical model is closer to the real terrain, the difference between the two calculation results directly reflects the approximation error of the simplified model. The system performs numerical calculations on these two influence values, subtracting the first gravity terrain correction influence value from the second gravity terrain correction influence value to obtain the difference. To more accurately quantify the impact of this difference on the model's approximation, the system can further process the difference by combining the terrain features corresponding to the grid cell (such as terrain roughness, slope, etc.) and preset error weighting coefficients, ultimately obtaining the physical model's approximation residual. During the calculation process, the system automatically filters out extreme differences caused by data anomalies (such as excessively large elevation deviations at individual sampling points), ensuring the reliability of the residual results through smoothing processing and avoiding abnormal data from misleading subsequent grid refinement decisions.

[0078] S105. When the approximation residual exceeds the preset accuracy threshold, the grid cells are refined and divided to generate multiple new computational grids.

[0079] Among them, the preset accuracy threshold refers to the value set in advance by the system to determine whether the approximation residual is within an acceptable range. This threshold is determined according to the specific gravity exploration task requirements (such as mineral exploration requiring higher accuracy, so the threshold is set smaller; conventional topographic mapping requires relatively lower accuracy, so the threshold is set larger) and industry gravity standards. For example, in geological hazard assessment, the preset accuracy threshold can be set to 0.01mGal.

[0080] Specifically, the system compares the approximation residual calculated for each grid cell with a preset accuracy threshold. If the approximation residual of a grid cell is less than or equal to the preset accuracy threshold, it means that the approximation error of the first physical model of that grid cell is within an acceptable range and no further refinement is required. The system marks this grid cell as an "accuracy-compliant cell." If the approximation residual exceeds the preset accuracy threshold, the system determines that the current size of the grid cell cannot meet the accuracy requirements and further refinement is necessary.

[0081] Before refining the grid cell, the system first obtains the initial size, location coordinates, and distribution of surrounding grid cells to ensure seamless integration between the new grid cells and their neighbors, avoiding overlaps or gaps. Then, the system uses a regular partitioning algorithm (such as quadtree partitioning) to uniformly divide the original grid cell into multiple new computational grids of the same size. The number of partitions is determined by the size of the original grid cell and the preset refinement level (e.g., dividing a 100m×100m grid cell into four 50m×50m grids). After partitioning, the system assigns a unique identifier to each new computational grid, recording its location coordinates, the information of its original grid cell, and the current subdivision level for subsequent data tracking and further refinement decisions.

[0082] S106. Repeat the steps of geometric fitting, gravity effect calculation, approximation residual and refinement segmentation for each new computational grid until the approximation residual corresponding to all grid cells does not exceed the preset accuracy threshold.

[0083] Each new computational grid refers to a smaller grid cell generated by refining and dividing the original grid cells whose approximation residuals exceeded the limit.

[0084] Specifically, the system treats each newly generated computational grid as an independent processing unit. Following the same process as the original grid unit, it sequentially performs real terrain data extraction, dual physical model (first physical model and second physical model) construction, gravity terrain correction impact value calculation, and approximation residual calculation to obtain the new approximation residual for each new computational grid. Subsequently, the system compares the new approximation residual of each new computational grid with a preset accuracy threshold. If the new approximation residual of a new computational grid is less than or equal to the preset accuracy threshold, it indicates that the accuracy of the new computational grid has met the standard, and it is marked as a "final qualified grid unit" without further refinement. If the new approximation residual still exceeds the preset accuracy threshold, the system determines that the size of the new computational grid still cannot meet the accuracy requirements, and it needs to be further refined and segmented to generate a new generation of smaller computational grids. During the iterative refinement process, the system records the number of refinements (i.e., subdivision levels) for each new computational grid and sets a maximum subdivision level limit (this limit is determined based on the accuracy of the digital elevation model; for example, when the accuracy of the digital elevation model is 30m, the grid size corresponding to the maximum subdivision level is not less than 30m). If a new computational grid reaches the maximum subdivision level but the new approximation residual still exceeds the limit, the system will mark it as a "precision non-convergence unit" and terminate further subdivision of that unit. The system will also note the location and error of this unit in subsequent reports for manual review. The system continues this iterative process until the new approximation residuals of all new computational grids do not exceed the preset accuracy threshold, or all non-compliant new computational grids have reached the maximum subdivision level. At this point, the refinement operation stops, and the system proceeds to the subsequent results aggregation stage.

[0085] In the above embodiment, the system first generates a non-uniform initial computational grid based on the local terrain complexity, making the grid size inversely correlated with the terrain complexity, thus achieving preliminary optimized allocation of computational resources. Then, it performs dual modeling by constructing first and second physical models with different complexities and calculates the difference in gravity effects between the two to obtain the approximation residual. This difference calculation mechanism can more accurately quantify the model approximation error of each grid cell. When the residual exceeds a preset accuracy threshold, the system adaptively refines the grid and iterates until the accuracy requirements are met. This adaptive refinement mechanism ensures that a denser grid is used in complex terrain areas to guarantee accuracy, while a sparser grid is used in simple terrain areas to save computational resources. Ultimately, it can improve overall computational efficiency while ensuring global computational accuracy, thereby addressing the computational redundancy problem of traditional uniform grid methods.

[0086] The following provides a more detailed description of the process of the method provided in this implementation. Please refer to [link / reference]. Figure 2 This is another flowchart illustrating a gravity terrain correction accuracy constraint method in an embodiment of this application.

[0087] S201. Obtain sampling point data of the digital elevation model within the grid cell coverage area.

[0088] Among them, sampling point data refers to the discretized terrain feature data extracted from the digital elevation model and distributed within the grid cell coverage area. The core data includes the three-dimensional coordinates (latitude, longitude, and elevation) of the sampling points, and may also include derived information such as slope and aspect, which are used to accurately characterize the terrain details within the grid cell.

[0089] Specifically, the system determines the boundary range of the target grid cell (including latitude, longitude, and elevation range) and performs consistency verification with the coordinate system of the digital elevation model to ensure the accuracy of data extraction. Then, based on the size of the grid cell and the terrain complexity, the system adaptively determines the density and distribution of sampling points—grid cells with high terrain complexity (such as steep slopes and canyon areas) will have a higher sampling point density, while grid cells with gentle terrain will have a lower density, in order to control the data volume while ensuring data representativeness. Next, the system uses spatial indexing technology to quickly locate the corresponding region of the grid cell in the digital elevation model, extracts all sampling point data that meet the distribution requirements within that region, and performs preliminary data cleaning to remove invalid sampling points with abnormal coordinates or elevation values ​​outside the reasonable range, ultimately obtaining a high-quality sampling point dataset for geometric fitting.

[0090] S202. Based on the sampling point data, the terrain within the grid cell is simplified into a flat-topped prism with uniform height to obtain the first physical model. Then, the sampling point data is fitted with a high-order surface to generate a second physical model constructed from the surface prism.

[0091] S203. Based on the distance between the grid cells and the observation points, the computational region is divided into the near zone, the middle zone, the far zone 1, and the far zone 2.

[0092] Specifically, the system first acquires the coordinates of the center point of each grid cell and the coordinates of the observation point. Using spatial coordinate calculation methods, it accurately calculates the straight-line distance between the center point of each grid cell and the observation point. Then, based on pre-defined distance division criteria, the system divides the spatial range of all grid cells into near zone, middle zone, far zone one, and far zone two. The division criteria comprehensively consider the physical characteristics of gravitational field decay with distance, the degree of influence of different regions on the gravity value of the observation point, and the accuracy and efficiency requirements of subsequent calculation algorithms. For example, based on extensive experimental data and theoretical analysis, the system defines the near zone as 0-50m from the observation point, the middle zone as 50m ≤ distance < 2km, the far zone one as 2km ≤ distance < 20km, and the far zone two as 20km ≤ distance < 166.7km and beyond. This regional division lays the foundation for using calculation algorithms with different accuracies for different regions, maximizing computational efficiency while ensuring computational accuracy.

[0093] S204. For the grid cells in the near zone, middle zone, and far zone, a high-precision integral algorithm is used to calculate the gravity terrain correction value based on the first physical model and the second physical model, respectively.

[0094] Specifically, for each grid cell in the near, middle, and far zones, the system retrieves parameter information from the previously constructed first physical model (flat-topped prism model) and second physical model (curved prism model), including the model's dimensions, height, and surface equation coefficients. Then, the high-precision integration algorithm module is activated, substituting the physical model parameters into the gravity field integration calculation formula. During the calculation, the system finely divides the integration interval and increases the number of integration nodes to improve the accuracy of the integration calculation and reduce errors caused by integration approximation. For the first physical model, which is a regular flat-topped prism, the integration calculation is relatively simple, and the system can quickly complete the calculation of the gravity terrain correction influence value. For the second physical model, which is a curved prism with a relatively complex shape, the system adopts a more refined integration strategy, such as an adaptive integration method, dynamically adjusting the integration node density according to the changes in the surface to ensure the accuracy of the integration calculation. By calculating the first and second gravity terrain correction influence values ​​for each grid cell separately, the system obtains the first and second gravity terrain correction influence values.

[0095] In some embodiments, this step can be implemented in multiple ways: Optionally, the system retrieves the first and second physical model parameters of the near-field, middle-field, and far-field grid cells, including the length, width, and height of the flat-topped prism and the surface equation coefficients of the curved prism; it starts the Gaussian integration module in the high-precision integration algorithm, sets the integration order to 20 to ensure calculation accuracy, substitutes the physical model parameters into the gravity field integration formula, and performs integration calculations on the two physical models of each grid cell separately; during the calculation process, the system monitors the integration convergence in real time, and if the integration result does not reach the preset convergence accuracy, it automatically increases the integration order until the accuracy requirement is met, and finally outputs the first and second gravity terrain correction influence values ​​of each grid cell.

[0096] Optionally, the system employs GPU acceleration technology to process grid cell data in the near, middle, and far regions in parallel, while simultaneously retrieving physical model parameters from multiple grid cells. It utilizes a high-precision Monte Carlo integration algorithm to generate a large number of random sampling points and perform random sampling calculations on the gravity field integration region of each physical model. A threshold for the number of sampling points is set; when the standard deviation of the calculation results is less than a preset accuracy threshold, sampling stops and the final gravity terrain correction impact value is output, thus obtaining the calculation results for each grid cell based on two physical models.

[0097] It is understandable that this step can also be achieved in other ways, such as combining the finite element analysis method to discretize the physical model into a finite number of elements for integral calculation, which is not limited here.

[0098] It should be noted that, in some cases, even if the overall approximation residual of a grid cell meets the accuracy requirements, there may still be local terrain abrupt changes within it. Relying solely on the overall residual may miss these critical areas that need refinement. Therefore, this application also provides a secondary verification mechanism (as shown in steps S205-S208) as a supplement and enhancement to the main accuracy judgment.

[0099] S205. When the approximation residual of the physical model does not exceed the preset accuracy threshold and meets the preset refinement conditions, the current grid cell is temporarily divided into four sub-grid cells to form a set of sub-grid cells.

[0100] Specifically, the system determines whether the approximation residual of each grid cell does not exceed a preset accuracy threshold. If this condition is met, the system further determines whether the grid cell meets preset refinement conditions. When determining the preset refinement conditions, the system calculates the roughness index of the terrain elevation data within the grid cell. This index is determined based on the standard deviation of the terrain elevation data. When the roughness index is not less than a preset roughness threshold, the system determines that the preset refinement conditions are met. After meeting the conditions, the system temporarily and uniformly divides the current grid cell into four sub-grid cells with the same geometric shape and area, forming a set of sub-grid cells, and assigns a unique identifier to each sub-grid cell.

[0101] Next, the system calculates the gravity effect value of each sub-grid cell on the observation point based on the first physical model corresponding to each sub-grid cell (constructed in the same way as the first physical model of the original grid cell, i.e., simplified to a flat-topped prism). It uses the same calculation method as for calculating the first gravity terrain correction influence value of the original grid cell. After calculation, the system adds the gravity effect values ​​of all sub-grid cells to obtain the sum of gravity effect values ​​and compares it with the first gravity terrain correction influence value of the original grid cell. If the difference is greater than a preset verification threshold, it indicates that the simplified model of the original grid cell may not accurately reflect the influence of local terrain on gravity. In this case, the system calculates the mean of the gravity effect values ​​of all sub-grid cells, compares the gravity effect value of each sub-grid cell with the mean, and marks sub-grid cells with differences exceeding a preset difference threshold as cells to be refined. If the difference is less than or equal to the preset verification threshold, it indicates that the simplified model of the original grid cell can reflect the influence of terrain on gravity well, and no further refinement is needed. The system then terminates further segmentation of the current grid cell.

[0102] S206. Calculate the gravity effect value of the observation point corresponding to each sub-grid cell based on the first physical model corresponding to each sub-grid cell.

[0103] The gravity effect value refers to the additional influence of the terrain quality represented by the sub-grid cell on the gravity observation value at the observation point. The calculation logic of this value is completely consistent with the calculation logic of the first gravity terrain correction influence value of the initial grid cell.

[0104] Specifically, for each sub-grid cell, the system extracts sampling point data from the digital elevation model within its coverage area. The density of the sampling points matches the size of the sub-grid cell to ensure that the terrain features within the sub-grid cell can be represented. For example, a 30m × 30m sub-grid cell corresponds to 9 evenly distributed sampling points. Next, following the rules for constructing the first physical model, the system calculates the arithmetic mean of the elevation values ​​of the sampling points within each sub-grid cell. This average value is used as the height of the flat-topped prism, and the length and width of the sub-grid cell are used as the base dimensions of the flat-topped prism to complete the construction of the first physical model for each sub-grid cell. Subsequently, the system determines the region (near zone, middle zone, far zone 1, or far zone 2) to which the sub-grid cell belongs in terms of distance from the observation point. Using the same algorithm as the initial grid cell for calculating the first gravity terrain correction influence value (high-precision integral algorithm for near zone, middle zone, and far zone 1, and fast approximation algorithm for far zone 2), the system substitutes the parameters (length, width, and height) of the first physical model of the sub-grid cell to calculate the gravity effect value of each sub-grid cell on the observation point, ensuring the consistency of the calculation logic and the comparability of the results.

[0105] S207. When the difference between the sum of gravity effect values ​​and the first gravity terrain correction influence value is greater than the preset verification threshold, the sub-grid unit corresponding to the target gravity effect value whose difference from the mean of gravity effect values ​​exceeds the preset difference threshold is marked as a unit to be refined.

[0106] Specifically, the system sums the gravity effect values ​​of the four sub-grid cells to obtain the total gravity effect value. Then, it retrieves the first gravity terrain correction influence value of the original grid cell and calculates the difference between the two (taking the absolute value to eliminate the influence of positive and negative signs). Next, it compares this difference with a preset verification threshold: if the difference is greater than the preset verification threshold, it indicates that the simplification of the local terrain in the first physical model (flat-topped prism) of the original grid cell may be biased and cannot accurately reflect the influence of the local terrain on gravity. At this time, the system calculates the arithmetic mean of the gravity effect values ​​of the four sub-grid cells, and then calculates the difference between the gravity effect value of each sub-grid cell and this mean (again, taking the absolute value). Next, each difference is compared with a preset difference threshold. Gravity effect values ​​with differences exceeding the preset difference threshold are identified as target gravity effect values, and their corresponding sub-grid cells are the cells to be refined. The system adds a special label (such as "to be refined-001") to these cells for further geometric fitting, gravity effect calculation, and refinement segmentation operations. If the difference is less than or equal to the preset verification threshold, it means that the simplified model of the original grid cells can better reflect the influence of the overall and local terrain on gravity. There is no need to label the cells to be refined, and the process can directly proceed to the subsequent termination segmentation process.

[0107] S208. When the difference between the sum of gravity effect values ​​and the first gravity terrain correction effect value is less than or equal to the preset verification threshold, the further segmentation of the current grid cell is terminated.

[0108] Specifically, after calculating and summing the gravity effect values ​​of the sub-grid units, the system calculates the difference (taking the absolute value) between the sum and the first gravity terrain correction influence value of the original grid unit. Then, the system retrieves a preset verification threshold (this threshold is set before calculation based on the accuracy level of the gravity exploration task, such as 0.01mGal for mineral exploration tasks and 0.03mGal for conventional topographic mapping), and compares the calculated difference with the threshold. If the difference is less than or equal to the preset verification threshold, the system determines that the first physical model of the original grid unit can fully represent the terrain features within its coverage area, the calculation accuracy of the gravity effect at the observation point meets the requirements, and there is no accuracy deviation due to terrain simplification in the local area. At this point, the system triggers a "Terminate Segmentation" command, marks the current grid unit's status as "Verified and Meets Standard," and stores its corresponding first gravity terrain correction influence value in the "Meets Standard Grid Data Pool." Subsequent processing only requires summarizing this data; no further refinement of the current grid unit is needed, thus saving computational resources and accelerating the overall process.

[0109] S209. Summarize the gravity terrain correction influence values ​​corresponding to all grid cells that meet the preset accuracy threshold, and generate a global terrain correction influence value matrix.

[0110] Specifically, the system first retrieves all grid cell data that meet the preset accuracy threshold from the "compliant grid data pool," including the geographic coordinate range (longitude range, latitude range), unique identifier (e.g., "G001 - R01C01," representing the first grid cell in the first row and first column), and the corresponding gravity terrain correction value for each grid cell. Then, based on the geographic coordinate range of all grid cells, the system determines the boundary of the global calculation area (taking the minimum and maximum longitude and latitude values ​​of all grid cells as the global boundary), and assigns each grid cell a row and column position in a two-dimensional matrix according to the spatial arrangement of the grid cells (e.g., longitude from west to east, latitude from south to north). For example, the grid cell with the smallest longitude and latitude corresponds to the first row and first column of the matrix; for each increase in longitude by one grid width, the column number increments by 1, and for each increase in latitude by one grid width, the row number increments by 1. Next, the system fills the correction impact value of each grid cell into its corresponding matrix row and column positions. For "positional gaps" (where the grid cell does not actually exist) that may exist due to the non-uniform distribution of the grid, the system will mark them as "no data" or fill them with a default value (such as 0, which will be noted in subsequent explanations). Finally, after completing the filling of data for all qualified grid cells, a global terrain correction impact value matrix covering the entire calculation area and with a one-to-one correspondence between location and correction value is generated.

[0111] S210. Generate standardized terrain correction data files based on the correction influence value matrix.

[0112] Standardized terrain correction data files refer to files generated after format conversion and data encapsulation of the global terrain correction impact value matrix according to the common data format specifications of Geographic Information Systems (GIS) or geophysical software. These files can be directly read, parsed, and used by mainstream GIS software (such as ArcGIS and QGIS) or geophysical data processing software (such as GMT and Surfer) without additional format conversion. GIS or geophysical software data format specifications refer to industry-standard data storage standards, including but not limited to Shapefile format (including .shp, .shx, .dbf, etc.), GeoTIFF format (raster image format with geographic coordinate information), ASCII Grid format (raster data in plain text format), etc. The specifications define the data storage structure, coordinate information, attribute fields (such as grid cell number and correction impact value), etc.

[0113] Specifically, the system first reads the data content of the global terrain correction influence value matrix, including the number of rows and columns of the matrix, the geographic coordinates (longitude and latitude) of each element, and the correction influence value. Simultaneously, it extracts key spatial information such as the global boundary and grid cell size. Then, based on the user-preset target software type (or the default mainstream format, such as GeoTIFF), the system calls the corresponding format conversion module. If the target format is GeoTIFF, the system treats the correction influence value matrix as the "pixel values" of the raster data, and writes geographic coordinate information (such as the coordinate system WGS84) and grid cell size into the GeoTIFF file's header, ensuring the file contains complete spatial reference information. If the target format is Shapefile, the system converts each grid cell into a "polygon feature," writes the grid cell number and correction influence value as attribute fields into a .dbf file, and generates a .shp file to store the geometric coordinates of the polygon features and a .shx file to store index information. During the format conversion process, the system verifies the data to ensure the integrity of attribute fields (e.g., each grid cell has a corresponding correction impact value) and the accuracy of coordinate information (consistent with the geographic coordinates of the original matrix), thus avoiding data loss or deviation due to format conversion. Finally, a standardized terrain correction data file conforming to the target format specification is generated and stored in the specified path.

[0114] S211. Generate an accuracy certificate report based on the correction influence value matrix and the final cell grid distribution.

[0115] Specifically, the system first integrates two types of core data: one is the error data associated with the global terrain correction influence value matrix (the approximation residual of each grid cell), and the other is the final cell grid distribution data (grid size, subdivision level, coordinate range). Then, the system generates report content in modules:

[0116] The non-uniform grid distribution map generation module can draw the boundary of each grid unit on the geographic base map (such as a simplified administrative division map) according to the coordinate range and size of the final grid, distinguish the subdivision level with different colors (such as black for the initial grid, blue for the first refinement, and green for the second refinement), mark the grid size of key areas (such as mountainous areas with complex terrain and flat plains), and intuitively display the non-uniform distribution characteristics.

[0117] The precision heatmap generation module can associate the approximation residual of each grid cell with the corresponding geographic coordinates, establish a "coordinate-error" mapping relationship, and use a color gradient algorithm (such as Jet color mapping) to convert the error value into the corresponding color (blue for error ≤ 0.005mGal, yellow for 0.005-0.01mGal, and red for > 0.01mGal). The heatmap is then rendered on the geographic base map, clearly marking the error range corresponding to the color.

[0118] The key performance indicator calculation module can automatically calculate quantitative indicators such as "average calculation error across the entire domain" (summing up the approximation residuals of all grid cells and dividing by the number of grid cells), "percentage of qualified grid cells" (number of qualified cells / total number of cells × 100%), "total calculation time" (total time from initial grid generation to matrix output), and "maximum subdivision level" (the highest subdivision level among all grid cells). At the same time, it can calculate the average error for different terrain types (mountainous areas, plains, hills) and form an indicator statistics table.

[0119] Finally, the system integrates the above-mentioned graphics, statistical tables, and calculation instructions (such as data sources and algorithm selection criteria) according to a preset template (such as Word or PDF format) to generate a complete accuracy certificate report. The report includes a cover (project name, report generation time), table of contents, main text (grid distribution analysis, accuracy analysis, performance indicators), appendix (original data list), etc., which can be directly exported for user use.

[0120] It should be noted that the specific values ​​of the various preset thresholds involved in this embodiment, such as the preset accuracy threshold, the maximum subdivision level, and the preset verification threshold, are not fixed. Those skilled in the art can flexibly set them or obtain them through experimental calibration based on factors such as the specific application scenario, the accuracy requirements of the exploration task (e.g., detailed mineral exploration is higher than regional general exploration), the data quality of the digital elevation model used, and industry operating standards.

[0121] The gravity terrain correction accuracy constraint system of this invention is applied to electronic devices. Figure 3 A schematic diagram of the architecture of an electronic device suitable for implementing embodiments of the present invention is shown.

[0122] It should be noted that, Figure 3 The electronic device shown is merely an example and should not be construed as limiting the functionality and scope of use of the embodiments of the present invention.

[0123] Those skilled in the art will understand that all or part of the steps in the various methods of the above embodiments can be implemented by instructions (computer programs), or by instructions (computer programs) controlling related hardware. These instructions can be stored in a computer-readable storage medium and loaded and executed by a processor. The electronic device of this embodiment includes a storage medium and a processor, wherein the storage medium stores multiple instructions that can be loaded by the processor to execute any step of the method provided in the embodiments of the present invention.

[0124] Specifically, the storage medium and the processor are electrically connected directly or indirectly to enable data transmission or interaction. For example, these components can be electrically connected to each other via one or more signal lines. The storage medium stores computer-executable instructions that implement data access control methods, including at least one software functional module that can be stored in the storage medium in the form of software or firmware. The processor executes various functional applications and data processing by running the software program and module stored in the storage medium. The storage medium can be, but is not limited to, Random Access Memory (RAM), Read-Only Memory (ROM), Programmable Read-Only Memory (PROM), Erasable Programmable Read-Only Memory (EPROM), Electrically Erasable Programmable Read-Only Memory (EEPROM), etc. The storage medium stores the program, and the processor executes the program after receiving the execution instructions.

[0125] Furthermore, the software programs and modules within the aforementioned storage medium may also include an operating system, which may include various software components and / or drivers for managing system tasks (e.g., memory management, storage device control, power management, etc.) and can communicate with various hardware or software components to provide an operating environment for other software components. The processor may be an integrated circuit chip with signal processing capabilities. The aforementioned processor may be a general-purpose processor, including a Central Processing Unit (CPU), a Network Processor (NP), etc., which can implement or execute the methods, steps, and logic flowcharts disclosed in this embodiment. The general-purpose processor may be a microprocessor or any conventional processor.

[0126] Since the instructions stored in the storage medium can execute the steps in any of the methods provided in the embodiments of the present invention, the beneficial effects of any of the methods provided in the embodiments of the present invention can be achieved, as detailed in the preceding embodiments, and will not be repeated here.

[0127] The above description is merely a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A gravity terrain correction accuracy constraint method, applied to a gravity terrain correction accuracy constraint system, characterized in that, The method includes: A non-uniform initial computational grid is generated based on the local terrain complexity of the digital elevation model, wherein the digital elevation model is obtained by fusing terrain elevation data of at least two different resolutions, and the size of the initial computational grid is inversely correlated with the corresponding local terrain complexity. Geometric fitting is performed on the real terrain data in each grid cell of the initial computational grid to obtain a parameterized first physical model and a second physical model. The first physical model and the second physical model have different complexities in the same grid cell. The step of geometrically fitting the real terrain data within each grid cell of the initial computational grid to obtain a parameterized first physical model and a second physical model includes: acquiring sampling point data of the digital elevation model within the area covered by the grid cell; simplifying the terrain within the grid cell into a flat-topped prism with uniform height based on the sampling point data to obtain the first physical model, wherein the height of the flat-topped prism is the arithmetic mean of the elevation values ​​of the sampling points; and performing high-order surface fitting on the sampling point data to generate a second physical model constructed from the surface prism. The gravity effect of each grid cell on the observation point is calculated based on the first physical model and the second physical model respectively, to obtain the first gravity terrain correction influence value and the second gravity terrain correction influence value; The steps of calculating the gravity effect of each grid cell on the observation point based on the first physical model and the second physical model to obtain the first gravity terrain correction influence value and the second gravity terrain correction influence value include: dividing the calculation area into near zone, middle zone, far zone 1, and far zone 2 according to the distance between the grid cell and the observation point; for the grid cells in the near zone, middle zone, and far zone 1, using a high-precision integration algorithm to calculate the gravity terrain correction influence value based on the first physical model and the second physical model respectively, wherein the gravity terrain correction influence value includes the first gravity terrain correction influence value and the second gravity terrain correction influence value; for the grid cells in far zone 2, using a fast approximation algorithm to calculate the gravity terrain correction influence value based on the first physical model and the second physical model respectively. The approximation residual of the physical model is calculated based on the difference between the first gravity terrain correction influence value and the second gravity terrain correction influence value; When the approximation residual exceeds a preset accuracy threshold, the grid cells are refined and divided to generate multiple new computational grids; The steps of geometric fitting, gravity effect calculation, approximation residual and refinement segmentation are repeated for each new computational grid until the approximation residual corresponding to all grid cells does not exceed the preset accuracy threshold.

2. The method according to claim 1, characterized in that, The step of refining the grid cells to generate multiple new computational grids specifically includes: If the current subdivision level corresponding to the grid cell is less than the preset maximum subdivision level, the grid cell is cut into four sub-grid cells with the same geometric shape and equal area. The maximum subdivision level is set according to the accuracy of the digital elevation model. The subgrid cell is identified as the new computing grid, and the current subdivision level corresponding to the new computing grid is updated; If the current subdivision level corresponding to the grid cell is greater than or equal to the preset maximum subdivision level, the grid cell is marked as a non-convergent cell, and further subdivision of the non-convergent cell is terminated.

3. The method according to claim 1, characterized in that, After the step of calculating the approximation residual of the physical model based on the difference between the first gravity terrain correction influence value and the second gravity terrain correction influence value, the method further includes: When the approximation residual does not exceed the preset accuracy threshold and meets the preset refinement conditions, the current grid cell is temporarily divided into four sub-grid cells to form a set of sub-grid cells. The gravity effect value of the observation point corresponding to each sub-grid cell is calculated based on the first physical model corresponding to each sub-grid cell. The calculation method of the gravity effect value is the same as the calculation method of the first gravity terrain correction influence value. When the difference between the sum of the gravity effect values ​​and the first gravity terrain correction effect value is greater than a preset verification threshold, the sub-grid cell corresponding to the target gravity effect value whose difference from the mean of the gravity effect values ​​exceeds a preset difference threshold is marked as a cell to be refined. When the difference between the sum of the gravity effect values ​​and the first gravity terrain correction effect value is less than or equal to a preset verification threshold, further segmentation of the current grid cell is terminated.

4. The method according to claim 3, characterized in that, The steps that meet the preset refinement conditions specifically include: Calculate the roughness index of the terrain elevation data within the grid cell, wherein the roughness index is determined based on the standard deviation of the terrain elevation data; When the roughness index is less than the preset roughness threshold, it is determined that the mesh element does not meet the preset refinement condition; When the roughness index is not less than the roughness threshold, the mesh cell is determined to meet the preset refinement conditions.

5. The method according to claim 1, characterized in that, After repeatedly performing the steps of geometric fitting, gravity effect calculation, approximation residual, and refinement segmentation on each new computational grid until the approximation residual corresponding to all grid cells does not exceed the preset accuracy threshold, the method further includes: The gravity terrain correction influence values ​​corresponding to all grid cells that meet the preset accuracy threshold are summarized to generate a global terrain correction influence value matrix. The correction influence value matrix is ​​the final gravity terrain correction result covering the entire calculation area. A standardized terrain correction data file is generated based on the correction influence value matrix. The standardized terrain correction data file is an output file that conforms to the data format specifications of GIS or geophysical software. An accuracy certificate report is generated based on the correction influence value matrix and the final cell grid distribution. The accuracy certificate report includes a non-uniform grid distribution map showing the differences in grid density in different regions, an accuracy heatmap representing the magnitude of calculation errors in different regions, and key performance indicators.

6. A gravity terrain correction accuracy constraint system, characterized in that, The system includes: one or more processors and memory; The memory is coupled to the one or more processors, the memory being used to store computer program code, the computer program code including computer instructions, the one or more processors invoking the computer instructions to cause the system to perform the method as described in any one of claims 1-5.

7. A computer-readable storage medium comprising instructions, characterized in that, When the instruction is executed on a gravity terrain correction accuracy constraint system, the system performs the method as described in any one of claims 1-5.

8. A computer program product, characterized in that, When the computer program product is run on a gravity terrain correction accuracy constraint system, the system performs the method as described in any one of claims 1-5.

Images