Land gravity far area terrain correction method and system based on digital elevation model

Abstract

The invention discloses a land gravity far-area terrain correction method and system based on a digital elevation model. The method comprises the following steps of: performing weighted average calculation on a far first region terrain correction value and a far second region generalized terrain correction value of regular distribution measuring points near irregular distribution measuring points by adopting a bilinear interpolation algorithm on the irregular distribution measuring points to obtain a far region terrain correction value; on the basis of a calculation algorithm of a terrain unit gravitational potential vertical partial derivative of a local rectangular coordinate integral kernel of a measuring point under a spherical coordinate system, high-precision calculation and accumulation are carried out on terrain units one by one, the influence of earth surface bending is considered, and a high-precision and high-resolution digital elevation model can be utilized to the greatest extent. A high-precision digital elevation model represented by geographic latitude and longitude is uniformly subdivided for terrain interfaces of a far first area and a far second area, so that digital elevation units with different resolutions are ensured not to be overlapped or have gaps to the greatest extent, and a high-precision and high-resolution digital elevation model is adapted; and subsequent updating and maintenance of the digital elevation model and the remote terrain correction value library are facilitated.

Description

Technical Field

[0001] The present invention belongs to the technical field of regional gravity exploration, and in particular relates to a land gravity far-zone terrain correction method and system based on a digital elevation model. Background Art

[0002] Gravity exploration reflects the density changes of underground geological bodies through gravity anomalies. The original field observation values ​​contain anomalies caused by differences in air mass, observation planes, and near-surface undulating terrain. Therefore, a series of corrections are required for the field observation values. Among them, the purpose of terrain correction is to eliminate the influence of undulating terrain. The "Regional Gravity Survey Specification" (DZ / T0082-2021) stipulates that the terrain correction range is 166.7km, which is specifically divided into four parts: near zone (0m~50m or 100m), middle zone (50m or 100m~2km), far zone 1 (2km~20km) and far zone 2 (20km~166.7km); the accuracy of far zone 1 and far zone 2 of 1:250,000 scale gravity exploration is ±0.114×10 -5 m / s 2 and ±0.213×10 -5 m / s 2 This accuracy is a rough value obtained by comparing the results of manual calculations and computer calculations at several measuring points distributed in different geographical locations.

[0003] Initially, surveyors treated the geoid as an infinite plane and developed algorithms and programs based on plane rectangular coordinates and the principles of fast Fourier transforms. Subsequently, to account for the spherical effects of the actual Earth, a generalized terrain correction algorithm for the far second zone was designed based on spherical coordinates. A numerical library for this correction was also developed. This library has been integrated into the RGIS software, and the algorithm has been incorporated into the regional gravity exploration industry standard. RGIS's far first zone terrain correction algorithm is based on plane coordinates and utilizes a terrain database with a spatial resolution of 1 km × 1 km. RGIS's far second zone generalized terrain calculation is based on spherical approximation and rotational symmetry, deriving an analytical formula for the gravitational force exerted by the mass of a spherical ring or spherical shell at a measuring point on the Earth's rotational axis. For measuring points on non-rotating axes, coordinate transformations are used to obtain new coordinate values ​​in a spherical polar coordinate system with the measuring point as the pole. Therefore, the RGIS far second zone generalized terrain correction algorithm requires, for any measuring point, to reconstruct the coordinates of the terrain to be calculated. By presetting the appropriate number of rings and orientations, the elevation values ​​of the spherical shells after terrain coordinate reconstruction are estimated, and the gravitational field of the mass of these spherical shells is calculated using analytical expressions, thereby estimating the far second zone terrain correction value (terrain correction value) of the measuring point.

[0004] The numerical database and algorithm for far-field terrain correction have the following shortcomings: First, the spatial resolution of the far-field second-area terrain database integrated by RGIS is 5 arc-min × 5 arc-min. Due to the limitations of the spatial resolution of the terrain database at that time and the limitations of computing performance, digital elevation model data with higher spatial resolution was not used; second, according to the pre-set number of rings and orientations, after the terrain coordinates are reconstructed, some spherical shells are large in size, and their elevation values ​​may be the average of the elevation values ​​of multiple original terrain units (such as Figure 1 (As shown), although the RGIS Far Zone 2 terrain correction algorithm reduces the amount of calculation, it also reduces the accuracy of the calculation. Third, the spatial resolution of the digital terrain used for terrain correction in different ranges is different, and the calculation methods used are different, so there is an interface problem. The calculation process of the interface between the Far Zone 1 and the Central Zone only uses the average value of the elevation nodes within the range, which to a certain extent reduces the accuracy of the calculation results. Similarly, the data organization form of the terrain units used in the Far Zone 1 and the Far Zone 2 is different. The former is represented by plane coordinates, and the latter is represented by spherical coordinates. The connection between the two is a jagged circle composed of square terrain units, and there are overlaps or gaps in the terrain units involved in the calculation. Summary of the Invention

[0005] In view of this, the present invention provides a land gravity far-zone terrain correction method and system based on a digital elevation model, which not only takes into account the influence of the curvature of the earth's surface, but also makes maximum use of the high-precision, high-resolution digital elevation model, solves the terrain correction calculation problem at the interface between the middle zone and the far zone one, and the far zone one and the far zone two, and ensures to the greatest extent that there is no overlap or gap between digital elevation model units of different resolutions, adapts to high-precision, high-resolution digital elevation models, and facilitates subsequent updates and maintenance of digital elevation models and far-zone terrain correction numerical libraries. The following technical solutions are specifically adopted to achieve this.

[0006] In a first aspect, the present invention provides a method for correcting land gravity far-zone terrain based on a digital elevation model, comprising the following steps:

[0007] Obtain digital elevation models and regularly distributed measurement points, and construct digital terrain models within the far-field terrain correction range;

[0008] Extracting, from the digital terrain model, a first digital elevation unit (DEU) that satisfies a first-far-zone terrain correction range and a second DEU that satisfies a second-far-zone generalized terrain correction range of any measuring point in the regularly distributed measuring points based on the position information of the measuring point;

[0009] Obtaining a preset formula for calculating the vertical first-order partial derivative of the gravitational potential of a single digital terrain unit that takes into account the curvature of the earth's surface at the arbitrary measuring point, and calculating the vertical first-order partial derivatives of the gravitational potential of all digital terrain units in the far first zone and the far second zone of the arbitrary measuring point based on the first digital elevation unit and the second digital elevation unit, and accumulating and summing the calculated values ​​to obtain a far first zone terrain correction value and a far second zone generalized terrain correction value for the arbitrary measuring point, wherein the digital terrain model includes multiple digital terrain units;

[0010] Using a bilinear interpolation algorithm on irregularly distributed measuring points in the digital terrain model, a weighted average of the far first area terrain correction value and the far second area generalized terrain correction value of the regularly distributed measuring points near the irregularly distributed measuring points is calculated to obtain the far area terrain correction value of the irregularly distributed measuring points;

[0011] The land gravity far zone terrain correction value is determined based on the far zone one terrain correction value of the regularly distributed measuring points and the far zone two generalized terrain correction value.

[0012] As a preferred embodiment of the above technical solution, a digital elevation model and regularly distributed measurement points are obtained to construct a digital terrain model within the far-area terrain correction range, including:

[0013] Obtaining a high-precision, high-spatial-resolution digital elevation model and the geographic coordinates of regularly distributed measurement points, as well as a far-zone terrain correction range radius, wherein the far-zone terrain correction range radius includes a far-zone 1 terrain correction radius and a far-zone 2 generalized terrain correction radius;

[0014] Calculating a range of digital terrain data with high spatial resolution covering the regularly distributed measuring points, and calculating a range of digital terrain data with low spatial resolution covering the regularly distributed measuring points;

[0015] The digital terrain units corresponding to the high spatial resolution and the digital terrain units corresponding to the low spatial resolution are extracted respectively.

[0016] As a preferred embodiment of the above technical solution, the high spatial resolution is 15 arc-sec, the low spatial resolution is 1 arc-min, the first digital elevation unit is a high spatial resolution terrain unit, and the second digital elevation unit is a low spatial resolution terrain unit.

[0017] As a preferred embodiment of the above technical solution, based on the position information of any measuring point in the regularly distributed measuring points, extracting from the digital terrain model a first digital elevation unit that satisfies a far first zone terrain correction range of the arbitrary measuring point and a second digital elevation unit that satisfies a far second zone generalized terrain correction range, comprising:

[0018] Obtaining the position information of the regularly distributed measuring points, the spatial resolution of the far first area terrain correction range, the inner radius and outer radius of the far first area terrain correction range, and the outer radius of the far second area generalized terrain correction range;

[0019] Determine digital terrain units within the terrain correction range of the regularly distributed measuring points based on the information, wherein the digital terrain units include high spatial resolution terrain units covered by the annular terrain in the first far zone and low spatial resolution terrain units covered by the circular terrain in the second far zone;

[0020] The terrain of the far zone 1 and the middle zone is a circular interface, the terrain of the far zone 1 and the far zone 2 is a circular interface, and the minimum digital terrain unit of the digital elevation model is divided by longitude and latitude lines. The high spatial resolution terrain unit covered by the circular interface between the far zone 1 and the middle zone is uniformly subdivided to obtain a first terrain subunit, and the far zone terrain correction calculation is performed on the first terrain subunit distributed in the far zone 1 range of the measuring point;

[0021] The low spatial resolution terrain unit covered by the circular interface of the far zone 1 and the far zone 2 is subdivided to obtain the second terrain subunit, and the far zone 2 generalized terrain correction calculation is performed on the second terrain subunit distributed in the far zone 2 of the measuring point.

[0022] As a preferred embodiment of the above technical solution, the terrain of the far zone 1 and the middle zone is a circular interface, the terrain of the far zone 1 and the far zone 2 is a circular interface, the minimum digital terrain unit of the digital elevation model is divided by longitude and latitude lines, and the high spatial resolution terrain unit covered by the circular interface of the far zone 1 and the middle zone is uniformly subdivided to obtain a first terrain subunit, including:

[0023] The high spatial resolution terrain unit covered by the circular interface between the far first area and the middle area is uniformly divided to obtain the first terrain subunit;

[0024] The first terrain sub-unit that falls outside the inner radius of the circular terrain in the far zone is involved in the terrain correction calculation in the far zone.

[0025] As a preferred embodiment of the above technical solution, the low spatial resolution terrain unit covered by the circular interface of the far zone 1 and the far zone 2 is subdivided to obtain the second terrain sub-unit, and the far zone 2 generalized terrain correction calculation is performed on the second terrain sub-unit distributed in the far zone 2 range of the measuring point, including:

[0026] The low spatial resolution terrain units covered by the circular interface of the far zone 1 and the far zone 2 are divided to obtain low spatial resolution terrain sub-units;

[0027] The low spatial resolution terrain subunits falling outside the inner radius of the circular terrain in the far second zone and the low spatial resolution terrain subunits within the outer radius of the circular terrain in the far second zone are involved in the generalized terrain correction calculation of the far second zone to obtain the second terrain subunit.

[0028] As a preferred embodiment of the above technical solution, the far-zone terrain correction value is: the vertical first-order partial derivative value of the first gravitational potential generated when the elevation value of the digital terrain unit is greater than the surplus mass corresponding to the elevation value of the measuring point, or when the elevation value of the digital terrain unit is less than the deficit mass corresponding to the elevation value of the measuring point, wherein the measuring points include regularly distributed measuring points and irregularly distributed measuring points;

[0029] The generalized terrain correction value of the far second zone is: the vertical first-order partial derivative value of the second gravitational potential generated by the excess mass or deficit mass between the reference ellipsoid corresponding to the digital terrain unit and the top surface of the digital terrain unit.

[0030] As a preferred embodiment of the above technical solution, a calculation formula for the vertical first-order partial derivative of the gravitational potential of a preset single digital terrain unit at the arbitrary measuring point taking into account the curvature of the earth's surface is obtained, including:

[0031] A pair of longitude planes, a pair of latitude planes and a pair of concentric spheres are preset, and the corresponding integral space is [r1, r2] × [θ1, θ2] × [λ1, λ2]. The terrain correction value calculation expression of a single terrain unit is: [r1, r2] represents the distance from the bottom and top surfaces of the digital terrain unit to the center of the earth, [θ1, θ2] represents a pair of latitude values ​​that limit the digital terrain unit, and [λ1, λ2] represents a pair of longitude values ​​that limit the digital terrain unit;

[0032]

[0033] Among them, λ represents the longitude of the measuring point, θ represents the latitude of the measuring point, and r represents the distance from the measuring point to the center of the earth. represents the first-order partial derivative of the gravitational potential, i = 1, 2, 3 represent the north direction, east direction and radial direction of the local rectangular coordinate system of the measuring point, t represents the integration node in the longitude direction, j represents the integration node in the latitude direction, k represents the integration node in the radial direction, Δλ represents the longitude dimension of the digital terrain unit, Δθ represents the latitude dimension of the digital terrain unit, Δr represents the radial dimension of the digital terrain unit, G represents the gravitational constant, ρ represents the density; P(λ,θ,r) is the spherical coordinate of the measuring point, ρ s is the constant density of the digital terrain unit, (r′ k ,θ′ j ,λ′ t ) is the spherical coordinate of the integration node, N λ 、N θ and N rThey represent the maximum order used when applying the triple Gauss-Legendre integral algorithm, ω λt 、ω θj 、ω rk is the weight, according to the vertical partial derivative V3 of the calculated gravitational potential, l is used to represent the Euclidean distance between the integration node and the measuring point, then I 3D (r′ k ,θ′ j ,λ′ t ) is:

[0034]

[0035] Among them, I 3D (r′ k ,θ′ j ,λ′ t ) represents the integral kernel function of the local rectangular coordinate system of the measuring point, r′ k represents the radial distance from the integration node to the center of the earth, θ′ j represents the latitude of the integration node, λ′ t Indicates the longitude of the integration node.

[0036] As a preferred embodiment of the above technical solution, the bilinear interpolation algorithm includes:

[0037]

[0038] Among them, x and y represent the easting and northing coordinates in the plane rectangular coordinate system, respectively; the geographic coordinates of the irregular measuring point are (λ, θ), and δg(x, y) is the far-zone correction value; g(x1, y1), g(x2, y1), g(x1, y2), and g(x2, y2) are the far-zone terrain correction values ​​of the four elevation nodes near the irregular measuring point in the plane rectangular coordinate system; in the spherical coordinate system, the correspondence between the longitude and latitude of the four elevation nodes and the plane coordinates is x1→λ1, x2→λ2, y1→θ1, and y2→θ2, respectively.

[0039] In a second aspect, the present invention further provides a land gravity far-zone terrain correction system based on a digital elevation model, which is applied to the above-mentioned land gravity far-zone terrain correction method based on a digital elevation model, comprising:

[0040] Terrain data acquisition module, used to obtain digital elevation models and regularly distributed measurement points, and construct digital terrain models within the far-area terrain correction range;

[0041] a digital elevation extraction module for extracting, from the digital terrain model, a first digital elevation unit satisfying a far first zone terrain correction range and a second digital elevation unit satisfying a far second zone generalized terrain correction range of any measuring point in the regularly distributed measuring points based on the position information of the measuring point;

[0042] a terrain improvement value calculation module, configured to obtain a preset calculation formula for the vertical first-order partial derivative of the gravitational potential of a single digital terrain unit at the arbitrary measuring point, taking into account the curvature of the earth's surface; and, based on the first digital elevation unit and the second digital elevation unit, respectively calculate the vertical first-order partial derivatives of the gravitational potential of all digital terrain units in the far first zone and the far second zone of the arbitrary measuring point and cumulatively sum them to obtain a far first zone terrain correction value and a far second zone generalized terrain correction value for the arbitrary measuring point, wherein the digital terrain model includes multiple digital terrain units;

[0043] a bilinear interpolation calculation module, configured to use a bilinear interpolation algorithm to calculate the weighted average of the far first zone terrain correction value and the far second zone generalized terrain correction value near the irregularly distributed measuring points in the digital terrain model to obtain the far zone terrain correction value of the irregularly distributed measuring points;

[0044] The far-zone terrain correction value determination module is used to determine the far-zone terrain correction value of land gravity based on the far-zone one terrain correction value of the regularly distributed measuring points and the far-zone two generalized terrain correction value.

[0045] The present invention provides a method and system for land gravity far-zone terrain correction based on a digital elevation model. A formula for calculating the vertical first-order partial derivative of the gravitational potential of a single digital terrain unit at any measuring point is used to calculate the vertical first-order partial derivatives of the gravitational potential of all digital terrain units in the far first and far second zones of any measuring point based on a first digital elevation unit and a second digital elevation unit, and the first-order partial derivatives are accumulated and summed to obtain a far-zone terrain correction value and a far-zone generalized terrain correction value for the arbitrary measuring point. A bilinear interpolation algorithm is used for irregularly distributed measuring points in the digital terrain model to perform a weighted average calculation of the far-zone terrain correction values ​​and the far-zone generalized terrain correction values ​​of regularly distributed measuring points near the irregularly distributed measuring points to obtain the far-zone terrain correction value for the irregularly distributed measuring points. An algorithm for calculating the vertical partial derivative of the gravitational potential of terrain units based on the local rectangular coordinate integral kernel of the measuring points in a spherical coordinate system performs high-precision calculations and accumulations on each terrain unit within the terrain correction range. This method not only takes into account the influence of the curvature of the earth's surface but also maximizes the use of the high-precision, high-resolution digital elevation model. The terrain interfaces of Far Zone 1 and Far Zone 2 are uniformly divided using high-precision digital elevation models represented by geographic longitude and latitude to solve the terrain correction calculation problems at the interfaces between the central zone and Far Zone 1, and between Far Zone 1 and Far Zone 2. This ensures to the greatest extent that there is no overlap or gap between digital elevation units of different resolutions, adapts to high-precision and high-resolution digital elevation models, and facilitates subsequent updates and maintenance of digital elevation models and the far zone terrain correction numerical library. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments. It should be understood that the following drawings only illustrate certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other relevant drawings can be obtained based on these drawings without paying any creative work.

[0047] Figure 1 This is a schematic diagram of the generalized terrain correction algorithm for the far second area of ​​RGIS;

[0048] Figure 2 A flowchart of the land gravity far-zone terrain correction method based on the digital elevation model provided by the present invention;

[0049] Figure 3 Schematic diagram of the generalized terrain of the Far 1st Area and Far 2nd Area of ​​gravity exploration provided by the present invention;

[0050] Figure 4 A geometric diagram of the spherical shell piece and the local rectangular coordinate system of the calculation point provided by the present invention;

[0051] Figure 5 A schematic diagram of the horizontal adaptive subdivision of a spherical shell provided by the present invention;

[0052] Figure 6 A schematic diagram of the bilinear interpolation method provided by the present invention;

[0053] Figure 7 This is a structural block diagram of the land gravity far-zone terrain correction system based on the digital elevation model provided by the present invention. DETAILED DESCRIPTION

[0054] The following describes embodiments of the present invention in detail. Examples of the embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended only to explain the present invention and are not to be construed as limiting the present invention.

[0055] It should be noted that when an element is referred to as being "fixed to" another element, it may be directly on the other element or there may be an intermediate element. When an element is considered to be "connected to" another element, it may be directly connected to the other element or there may be an intermediate element. Conversely, when an element is referred to as being "directly on" another element, there is no intermediate element. The terms "vertical," "horizontal," "left," "right," and similar expressions used herein are for illustrative purposes only.

[0056] See Figure 2The present invention provides a method for correcting land gravity far-zone terrain based on a digital elevation model, comprising the following steps:

[0057] S1: Obtain digital elevation model and regularly distributed measurement points, and construct digital terrain model within the far-area terrain correction range;

[0058] S2: extracting, from the digital terrain model, a first digital elevation unit (DEU) that satisfies a first-far-zone terrain correction range and a second DEU that satisfies a second-far-zone generalized terrain correction range of any measuring point in the regularly distributed measuring points according to the position information of the measuring point;

[0059] S3: obtaining a preset formula for calculating the vertical first-order partial derivative of the gravitational potential of a single digital terrain unit taking into account the curvature of the earth's surface at the arbitrary measuring point; and calculating the vertical first-order partial derivatives of the gravitational potential of all digital terrain units in the far first zone and the far second zone of the arbitrary measuring point based on the first digital elevation unit and the second digital elevation unit, and accumulating and summing the calculated values ​​to obtain a far first zone terrain correction value and a far second zone generalized terrain correction value of the arbitrary measuring point, wherein the digital terrain model includes a plurality of digital terrain units;

[0060] S4: using a bilinear interpolation algorithm for the irregularly distributed measuring points in the digital terrain model, performing weighted average calculation on the far first area terrain correction value and the far second area generalized terrain correction value of the regularly distributed measuring points near the irregularly distributed measuring points to obtain the far area terrain correction value of the irregularly distributed measuring points;

[0061] S5: Determine the land gravity far zone terrain correction value based on the far zone one terrain correction value of the regularly distributed measuring points and the far zone two generalized terrain correction value.

[0062] In this embodiment, obtaining a digital elevation model and regularly distributed measurement points to construct a digital terrain model within a far-field terrain correction range includes: obtaining a high-precision, high-spatial-resolution digital elevation model and the geographic coordinates of the regularly distributed measurement points, as well as a far-field terrain correction range radius, wherein the far-field terrain correction range radius includes a far-field first-area terrain correction radius and a far-field second-area generalized terrain correction radius; calculating a range of high-spatial-resolution digital terrain data covering the regularly distributed measurement points, and calculating a range of low-spatial-resolution digital terrain data covering the regularly distributed measurement points; and extracting digital terrain units corresponding to the high and low spatial resolutions, respectively. The high spatial resolution is 15 arc-second, the low spatial resolution is 1 arc-minute, the first digital elevation unit is a high-spatial-resolution terrain unit, and the second digital elevation unit is a low-spatial-resolution terrain unit. A digital elevation model in geography or geology refers to a representation of the actual undulating surface using a set of regularly spaced numbers. A digital terrain model (or digital terrain data) contains a lot of discrete information. The smallest unit of information is called a digital elevation unit. The data it contains represents the elevation of a certain spatial location, also known as a digital terrain unit. A regularly distributed measuring point refers to a measuring point located at the top center of a digital terrain unit. A digital terrain unit is the smallest unit of a digital elevation model: a "single" digital terrain unit emphasizes the number of this smallest unit.

[0063] It should be noted that the above-mentioned land gravity far-zone terrain correction process mainly includes: obtaining a high-resolution digital elevation model and obtaining regularly distributed measuring points, constructing digital terrain units within the far-zone terrain correction range that meet the regularly distributed measuring points, including digital terrain data with a spatial resolution of 15 arc-sec and digital terrain data with a spatial resolution of 1 arc-min; according to the position information of any measuring point P in the regularly distributed measuring points, extracting a digital elevation unit with a spatial resolution of 15 arc-sec that meets the far-zone 1 terrain correction of the measuring point P from the pre-constructed terrain data, and extracting a digital elevation unit with a spatial resolution of 15 arc-sec that meets the far-zone 2 generalized terrain correction of the measuring point P. The digital elevation unit with an interval resolution of 1 arc-min is used; the vertical first-order partial derivative calculation formula of the gravitational potential of a single digital elevation unit at the measuring point P, which takes into account the curvature of the earth's surface, is used to calculate the vertical first-order partial derivatives of the gravitational potential of all digital elevation units in the far first and far second zones of the measuring point respectively and accumulate the sums to obtain the far first zone terrain correction value and the far second zone generalized terrain correction value of the measuring point P; for irregularly distributed measuring points, the bilinear interpolation method is used to calculate the far first zone terrain correction value and the far second zone generalized terrain correction value of the center position of the four neighboring digital elevation units in a 2×2 area near the measuring point through weighted average.

[0064] It should be understood that the calculation formula for the vertical first-order partial derivative of the gravitational potential of a single digital terrain unit at any measuring point is used to calculate the vertical first-order partial derivatives of the gravitational potential of all digital terrain units in the far first and far second zones of any measuring point based on the first digital elevation unit and the second digital elevation unit, respectively, and the cumulative sums are taken to obtain the far first zone terrain correction value and the far second zone generalized terrain correction value of the arbitrary measuring point. The far zone terrain correction value of the irregularly distributed measuring point in the digital terrain model is obtained by performing a weighted average calculation using a bilinear interpolation algorithm on the far first zone terrain correction value and the far second zone generalized terrain correction value of the regularly distributed measuring point near the irregularly distributed measuring point. The calculation algorithm for the vertical partial derivative of the gravitational potential of the terrain unit based on the local rectangular coordinate integral kernel of the measuring point in the spherical coordinate system performs high-precision calculations on each terrain unit within the terrain correction range, which not only takes into account the influence of the curvature of the earth's surface but also maximizes the use of the high-precision, high-resolution digital elevation model. A high-precision digital elevation model expressed in geographic longitude and latitude is uniformly adopted for the Far Zone 1 and Far Zone 2. The uniform division of the digital elevation model units solves the problem of terrain correction calculation at the interface between the central zone and Far Zone 1, and between Far Zone 1 and Far Zone 2. This ensures to the greatest extent that there is no overlap or gap between digital elevation units of different resolutions, adapts to high-precision and high-resolution digital elevation models, and facilitates subsequent updates and maintenance of digital elevation models and the far zone terrain correction numerical library.

[0065] Optionally, based on the position information of any measuring point in the regularly distributed measuring points, extracting from the digital terrain model a first digital elevation unit that satisfies a far first zone terrain correction range of the arbitrary measuring point and a second digital elevation unit that satisfies a far second zone generalized terrain correction range, including:

[0066] S10: Acquire the position information of the regularly distributed measuring points, the spatial resolution of the far first area terrain correction range, the inner radius and outer radius of the far first area terrain correction range, and the outer radius of the far second area generalized terrain correction range;

[0067] S11: Determine digital terrain units within a terrain correction range of regularly distributed measurement points based on the information, wherein the digital terrain units include high spatial resolution terrain units covered by the annular terrain in the first far zone and low spatial resolution terrain units covered by the circular terrain in the second far zone;

[0068] S12: The terrain of the far zone 1 and the middle zone is a circular interface, the terrain of the far zone 1 and the far zone 2 is a circular interface, and the minimum digital terrain unit of the digital elevation model is divided by longitude and latitude lines. The high spatial resolution terrain unit covered by the circular interface between the far zone 1 and the middle zone is uniformly subdivided to obtain a first terrain sub-unit, and the far zone terrain correction calculation is performed on the first terrain sub-unit distributed within the far zone 1 range of the measuring point;

[0069] S13: Subdivide the low spatial resolution terrain unit covered by the circular interface of the far zone 1 and the far zone 2 to obtain the second terrain sub-unit, and perform the far zone 2 generalized terrain correction calculation on the second terrain sub-unit distributed in the far zone 2 of the measuring point.

[0070] In this embodiment, the terrain of the far first zone and the middle zone is a circular interface, the terrain of the far first zone and the far second zone is a circular interface, and the smallest digital terrain unit of the digital elevation model is divided by longitude and latitude lines. The high spatial resolution terrain units covered by the circular interface between the far first zone and the middle zone are uniformly subdivided to obtain first terrain subunits, including: uniformly subdividing the high spatial resolution terrain units covered by the circular interface between the far first zone and the middle zone to obtain first terrain subunits; involving the first terrain subunits falling outside the inner radius of the circular ring terrain of the far first zone in the terrain correction calculation of the far first zone; and involving the high spatial resolution terrain units falling within the outer radius of the circular ring terrain of the far first zone in the terrain correction calculation of the far first zone with the terrain of the far first zone and the far second zone as a circular interface. The low spatial resolution terrain cells covered by the circular interface of the far first zone and the far second zone are subdivided to obtain second terrain sub-cells, and the far second zone generalized terrain correction calculation is performed on the second terrain sub-cells distributed within the far second zone of the measuring point, including: subdividing the low spatial resolution terrain cells covered by the circular interface of the far first zone and the far second zone to obtain low spatial resolution terrain sub-cells; and involving the low spatial resolution terrain sub-cells falling outside the inner radius of the circular ring terrain of the far second zone and the low spatial resolution terrain cells within the outer radius of the circular ring terrain of the far second zone in the far second zone generalized terrain correction calculation.

[0071] It should be noted that the far first zone terrain correction value is: the first vertical first-order partial derivative value of the gravitational potential generated when the elevation value of the digital terrain unit is greater than the excess mass corresponding to the elevation value of the measuring point, or when the elevation value of the digital terrain unit is less than the deficit mass corresponding to the elevation value of the measuring point, wherein the measuring points include regularly distributed measuring points and irregularly distributed measuring points; the far second zone generalized terrain correction value is: the second vertical first-order partial derivative value of the gravitational potential generated by the excess mass or deficit mass between the reference ellipsoid corresponding to the digital terrain unit and the top surface of the digital terrain unit. In other words, the above S1 mainly includes: obtaining a high-precision digital elevation model and the outer boundary geographic coordinate information of regularly distributed measuring points to obtain the far-zone terrain correction range; calculating the 15arc-sec digital terrain data range covering the regular measuring network, which is used to calculate the far-zone one terrain correction; calculating the 1arc-min digital terrain data range covering the regular measuring network, which is used to calculate the far-zone two generalized terrain correction, where the generalized terrain correction is the sum of the intermediate layer correction and the terrain correction; extracting the corresponding 15arc-sec and 1arc-min digital terrain data according to the calculated coverage range, thereby minimizing the number of files for reading data terrain and computer memory consumption.

[0072] Specifically, the above steps after S2 and before S3 mainly include: obtaining the geographic coordinates of any regularly distributed measuring points and the spatial resolution of the terrain data of the far zone, obtaining the inner radius and outer radius of the terrain correction ring range of the far zone, and the outer radius of the general terrain correction range of the far zone; obtaining the digital terrain units covered by the inner radius and outer radius of the circular terrain of the far zone, and the spatial position of the digital terrain units relative to the measuring point; the terrain of the far zone and the middle zone is a circular interface, and the digital elevation model data is represented as a square domain, so the digital terrain units covered by the inner radius of the far zone are evenly divided to obtain a terrain with a smaller volume. Subunits (first terrain subunits) that have been segmented and fall outside the inner radius of the Far Zone 1 participate in the terrain correction calculation of the Far Zone 1. The terrain of the Far Zone 1 and the Far Zone 2 forms a circular interface. The terrain units that fall within the outer radius of the Far Zone 1 participate in the terrain correction calculation of the Far Zone 1. The digital terrain units covered by the inner radius of the Far Zone 2 are segmented to obtain smaller terrain subunits (second terrain subunits). The segmented terrain subunits that fall outside the inner radius of the Far Zone 2 participate in the generalized terrain correction calculation of the Far Zone 2. The terrain units that fall within the outer radius of the Far Zone 2 participate in the generalized terrain correction calculation of the Far Zone 2.

[0073] Optionally, obtaining a preset calculation formula for the vertical first-order partial derivative of the gravitational potential of a single digital terrain unit at the arbitrary measuring point taking into account the curvature of the earth's surface comprises:

[0074] A pair of longitude planes, a pair of latitude planes and a pair of concentric spheres are preset, and the corresponding integral space is [r1, r2] × [θ1, θ2] × [λ1, λ2]. The terrain correction value calculation expression of a single terrain unit is: [r1, r2] represents the distance from the bottom and top surfaces of the digital terrain unit to the center of the earth, [θ1, θ2] represents a pair of latitude values ​​that limit the digital terrain unit, and [λ1, λ2] represents a pair of longitude values ​​that limit the digital terrain unit;

[0075]

[0076] Among them, λ represents the longitude of the measuring point, θ represents the latitude of the measuring point, and r represents the distance from the measuring point to the center of the earth. represents the first-order partial derivative of the gravitational potential, i = 1, 2, 3 represent the north direction, east direction and radial direction of the local rectangular coordinate system of the measuring point, t represents the integration node in the longitude direction, j represents the integration node in the latitude direction, k represents the integration node in the radial direction, Δλ represents the longitude dimension of the digital terrain unit, Δθ represents the latitude dimension of the digital terrain unit, Δr represents the radial dimension of the digital terrain unit, G represents the gravitational constant, ρ represents the density; P(λ,θ,r) is the spherical coordinate of the measuring point, ρ sis the constant density of the digital terrain unit, (r′ k ,θ′ j ,λ′ t ) is the spherical coordinate of the integration node, N λ 、N θ and N r are the orders used in the triple Gauss-Legendre integral algorithm, ω λt 、ω θj 、ω rk is the weight, according to the vertical partial derivative V3 of the calculated gravitational potential, l is used to represent the Euclidean distance between the integration node and the measuring point, then I 3D (r′ k ,θ′ j ,λ′ t ) is:

[0077]

[0078] Among them, I 3D (r′ k ,θ′ j ,λ′ t ) represents the integral kernel function of the local rectangular coordinate system of the measuring point, r′ k represents the radial distance from the integration node to the center of the earth, θ′ j represents the latitude of the integration node, λ′ t Indicates the longitude of the integration node.

[0079] In this embodiment, the above S3 mainly includes: the calculated value of the terrain correction in the far first zone is the vertical first-order partial derivative value of the gravitational potential generated by the surplus mass or the deficit mass in the digital terrain unit that is greater than the elevation value of the measuring point (this value is a positive number); the calculated value of the generalized terrain in the far second zone is the vertical first-order partial derivative value of the gravitational potential generated by the surplus mass or the deficit mass between the reference ellipsoid corresponding to the terrain unit and the top surface of the terrain unit (this value can be a positive number); the preset terrain correction value calculation formula (1) for a single terrain unit with an integral space of [r1, r2]×[θ1, θ2]×[λ1, λ2] defined by a pair of longitude planes, a pair of latitude planes, and a pair of concentric spheres.

[0080] It should be noted that the bilinear interpolation algorithm includes:

[0081]

[0082] Among them, x and y represent the easting and northing coordinates in the plane rectangular coordinate system, respectively; the geographic coordinates of the irregular measuring point are (λ, θ), and δg(x, y) is the far-zone correction value; g(x1, y1), g(x2, y1), g(x1, y2), and g(x2, y2) are the far-zone terrain correction values ​​of the four elevation nodes near the irregular measuring point in the plane rectangular coordinate system; in the spherical coordinate system, the correspondence between the longitude and latitude of the four elevation nodes and the plane coordinates is x1→λ1, x2→λ2, y1→θ1, and y2→θ2, respectively.

[0083] Specifically, see Figure 3 The Earth's reference ellipsoid is approximated as a sphere, with the center of the sphere as the Earth's center of mass and the spherical surface as the reference ellipsoid. Circles are constructed with an observation point P on the actual undulating surface as the center and radii of 2 km, 20 km, and 166.7 km on the spherical surface. In regional gravity specifications, far-zone terrain corrections include two parts: Far Zone 1 and Far Zone 2. Within the range of 2 km to 20 km, the gray area defined by the actual undulating surface relative to the plane passing through point P is considered Far Zone 1 terrain. Within the range of 20 km to 166.7 km, the striped area defined by the actual undulating surface and the reference ellipsoid passing through points M, M1, and M2 is considered Far Zone 2 generalized terrain.

[0084] Among them, the data of the current high-precision digital elevation model is generally expressed in LON (geodetic longitude), LAT (geodetic latitude) and H (altitude above sea level), and stored in GRID or GeoTiff format. In order to utilize these high-precision digital elevation models and take into account the influence of the curvature of the earth's surface, the terrain quality between 2km and 166.7km in the far area is divided according to the plane grid position of the corresponding digital elevation model. After the division, each (digital) terrain unit has the elevation value of the corresponding position of the digital elevation model. By adding the geoid disturbance value, the ellipsoid height of the (digital) terrain unit is obtained. Under the spherical approximation condition, each terrain unit is approximately regarded as a spherical shell, which is represented by a pair of geocentric longitude line planes (λ1, λ2), a pair of geocentric latitude line planes (θ1, θ2) and a pair of concentric spheres (r1, r2), such as Figure 4 As shown in the figure, r1 and r2 represent the distances from the bottom and top of the spherical shell to the center of the sphere, respectively. Since the terrain of each region can be decomposed into several such spherical shells, the gravitational force exerted by each spherical shell at the measuring point can be calculated and then added together to obtain the gravitational value of the terrain mass in the region at the measuring point.

[0085] Where, under the spherical approximation, the density is a constant ρ S The integral space is [r1, r2]×[θ1, θ2]×[λ1, λ2], point S is the geometric center of the spherical shell, the gravitational constant G, and the external gravitational potential expression at P(λ, θ, r) is:

[0086]

[0087] in, Represents the Euclidean distance between the measurement point and the source point, and cosψ=sinθsinθ S +cosθcosθ S cos(λ S -λ), the local rectangular coordinates of P(λ,θ,r) are defined as (f1,f2,f3), with point P as the origin, and f1, f2, and f3 are distributed as unit vectors pointing due north, due east, and radially outward from the center of the earth. Then any vector q in space is expressed in this local rectangular coordinate system as: q=q1f1+q2f2+q3f3. Among them, dλ S , dθ S , dr S is the integration variable, λ S represents the longitude coordinate of the integration point, θ S represents the latitude coordinate of the integration point, r S Represents the radial coordinate of the integration point, q1 represents the north coordinate of the local rectangular coordinate system of the measuring point P, q2 represents the east coordinate of the local rectangular coordinate system, and q3 represents the radial coordinate of the local rectangular coordinate system.

[0088] The partial derivative of the gravitational potential is obtained at the local rectangular coordinates (f1, f2, f3) of point P. It is necessary to convert the global coordinate system l PS Convert to the local rectangular coordinate system of point P. The unit vector of the global rectangular coordinate system and the unit vector of the local rectangular coordinate system of point P satisfy the following conversion:

[0089] (f1,f2,f3) T =K y R y (90°-θ)R z (180°-λ)(e x ,e y ,e z ) T (5)

[0090] Where R z (α) is the transformation matrix around the z-axis. If the coordinates before the transformation are regarded as a whole, then the matrix represents the clockwise rotation of the entire coordinate around the z-axis with an angle of α. The rotation matrix R y (β) realizes the clockwise rotation of the entire coordinate around the y-axis, with a rotation angle of β; K y The effect is to reverse the positive direction of the y-axis, e x 、e y 、e zRespectively represent the coordinates in the x, y, and z directions of the global rectangular coordinate system. The order of action of these three matrices is R z (α), R y (β), and finally K y Their specific expressions are listed below, R z The specific expression of (α):

[0091]

[0092] R y The specific expression of (β):

[0093]

[0094] and K y The specific expression of:

[0095]

[0096] The product K of three matrices y R y (90°-θ)R z The result of (180°-λ) is recorded as In the global coordinate system l PS Converted to the local rectangular coordinate system of point P, the result is marked as (△q1,△q2,△q3), which satisfies the following expression:

[0097]

[0098] From formula (4), we can know that the gravitational potential V generated by the spherical shell at P(λ,θ,r) is sph The first-order partial derivative of (λ,θ,r) in the local rectangular coordinate system is expressed as:

[0099]

[0100] Where:

[0101]

[0102] Among them, △q1, △q2, △q3 represent the distance / PS In the local rectangular coordinate system of point P, the northing, easting and radial components are represented by grad, V1 represents the first-order derivative of the potential field scalar toward the north, V2 represents the first-order derivative of the potential field scalar toward the east, and V3 represents the first-order derivative of the potential field scalar toward the radial outward direction.

[0103] Then, the partial derivative of the gravitational potential of any constant density spherical shell based on the local rectangular coordinate system is:

[0104]

[0105] In the formula, i∈1,2,3, where the subscripts 1 and 2 refer to the north and east directions respectively, and 3 refers to the radial outward direction. The triple Gauss-Legendre integral method is used to approximate the first-order partial derivative of the gravitational potential of any constant-density spherical shell at the measuring point:

[0106]

[0107] Where N λ 、N θ and N r are the orders used in the triple Gauss-Legendre integration method; ω λi 、ω θj 、ω rk is the weight coefficient; for calculating the vertical partial derivative of the gravitational potential V3, I 3D (r′ k ,θ′ j ,λ′ t )The specific expression is:

[0108]

[0109] The parameter (r′ k ,θ′ j ,λ′ t ) are the spherical coordinates of the Gaussian nodes or equivalent point sources, which are:

[0110]

[0111] Based on formula (10) and the superposition principle, the mass of all spherical shells within a certain range of the measuring point can be obtained by summing. When obtaining the far-field terrain of regularly distributed (grid) measuring points, the elevation value of the measuring point is generally set to the ellipsoid height of the digital elevation model at the same plane position, that is, the undulating terrain represented by the digital elevation model. In this case, the calculation accuracy of the gravity simulation results obtained directly using the Gauss-Legendre integral method is low. In order to solve this problem, it is necessary to subdivide the model into smaller models, which is equivalent to increasing the ratio parameter D of the measuring point-source point distance to the source body size, thereby obtaining higher calculation accuracy. The present invention uses horizontal adaptive subdivision technology, and the size L of the spherical shell (topographic unit) is λ 、L θ The definition is as follows:

[0112]

[0113] Whether a (digital) terrain unit needs horizontal adaptive subdivision can be determined according to the following inequality, where distance is used to represent the distance between the measuring point and the source body, L λ 、L θThey represent the horizontal dimensions of the spherical shell. When the inequality (15) holds, no subdivision is performed. If the inequality does not hold, the subdivision process is performed. The subdivision process is as follows: Figure 5 shown.

[0114]

[0115] The distance from the measuring point to the center of the top surface of the terrain unit is expressed as:

[0116]

[0117] Where, L i The specific expression is shown in formula (16), where i∈{λ,θ}, cosψ′ represents the cosine of the angle between the vector from the center of the earth to the calculation point and the vector from the center of the earth to the center of the top surface of the spherical shell. The horizontal adaptive subdivision technology combined with the Gauss-Legendre integral can ensure the calculation accuracy of the gravity simulation under the condition that the calculation point is close to the spherical shell. In the process of designing the calculation program, a "stack" space is opened up and the geometric information of the subdivided spherical shell is stored in sequence by "pushing"; these subdivided spherical shells are called in sequence, and the need for further subdivision is determined based on their size information. If not, the Gauss-Legendre integral method is directly applied for calculation; according to the above principle, the gravity values ​​of all subdivided spherical shells of the original spherical shell are superimposed to form the gravity value of the original spherical shell. The above calculation process requires presetting the capacity of this "stack" and the minimum horizontal size to terminate the subdivision process. Based on many practices, the "stack" capacity in the program design is generally set to 50, and the minimum horizontal size is generally set to 1 / 8 or 1 / 16 of the spatial resolution of the digital elevation model.

[0118] In the plane rectangular coordinate system, for irregular measuring points, first obtain the far-field terrain correction values ​​g(x1,y1), g(x2,y1), g(x1,y2), g(x2,y2) of the four regular elevation nodes near the measuring point, and their distribution relationship with the measuring point (such as Figure 6 As shown). In the spherical coordinate system, the corresponding relationship between the longitude and latitude of the above four elevation nodes and the horizontal coordinates is x1→λ1, x2→λ2, y1→θ1, y2→θ2. The geographic coordinates of the irregularly distributed measurement points are (λ,θ). Its far-zone correction value δg(x,y) is calculated using the following bilinear interpolation formula:

[0119]

[0120] See Figure 7 The present invention further provides a land gravity far-zone terrain correction system based on a digital elevation model, which is applied to the above-mentioned land gravity far-zone terrain correction method based on a digital elevation model, comprising:

[0121] Terrain data acquisition module, used to obtain digital elevation models and regularly distributed measurement points, and construct digital terrain models within the far-area terrain correction range;

[0122] a digital elevation extraction module for extracting, from the digital terrain model, a first digital elevation unit satisfying a far first zone terrain correction range and a second digital elevation unit satisfying a far second zone generalized terrain correction range of any measuring point in the regularly distributed measuring points based on the position information of the measuring point;

[0123] a terrain improvement value calculation module, configured to obtain a preset calculation formula for the vertical first-order partial derivative of the gravitational potential of a single digital terrain unit at the arbitrary measuring point, taking into account the curvature of the earth's surface; and, based on the first digital elevation unit and the second digital elevation unit, respectively calculate the vertical first-order partial derivatives of the gravitational potential of all digital terrain units in the far first zone and the far second zone of the arbitrary measuring point and cumulatively sum them to obtain a far first zone terrain correction value and a far second zone generalized terrain correction value for the arbitrary measuring point, wherein the digital terrain model includes multiple digital terrain units;

[0124] a bilinear interpolation calculation module, configured to use a bilinear interpolation algorithm to calculate the weighted average of the far first zone terrain correction value and the far second zone generalized terrain correction value near the irregularly distributed measuring points in the digital terrain model to obtain the far zone terrain correction value of the irregularly distributed measuring points;

[0125] The far-zone terrain correction value determination module is used to determine the far-zone terrain correction value of land gravity based on the far-zone one terrain correction value of the regularly distributed measuring points and the far-zone two generalized terrain correction value.

[0126] In all examples shown and described herein, any specific values ​​should be interpreted as merely exemplary and not limiting, and thus other examples of the exemplary embodiments may have different values.

[0127] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, it does not need to be further defined or explained in subsequent drawings.

[0128] The above-described embodiments merely illustrate several implementations of the present invention, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the present invention. It should be noted that variations and modifications are possible without departing from the scope of the present invention, and such variations and modifications are fully within the scope of protection of the present invention.

Claims

1. A method for land gravity far-zone terrain correction based on digital elevation model, characterized in that: The following steps are involved: Obtain digital elevation models and regularly distributed measurement points, and construct digital terrain models within the far-field terrain correction range; Extracting, from the digital terrain model, a first digital elevation unit (DEU) that satisfies a first-far-zone terrain correction range and a second DEU that satisfies a second-far-zone generalized terrain correction range of any measuring point in the regularly distributed measuring points based on the position information of the measuring point; Obtaining a preset formula for calculating the vertical first-order partial derivative of the gravitational potential of a single digital terrain unit that takes into account the curvature of the earth's surface at the arbitrary measuring point, and calculating the vertical first-order partial derivatives of the gravitational potential of all digital terrain units in the far first zone and the far second zone of the arbitrary measuring point based on the first digital elevation unit and the second digital elevation unit, and accumulating and summing the calculated values ​​to obtain a far first zone terrain correction value and a far second zone generalized terrain correction value for the arbitrary measuring point, wherein the digital terrain model includes multiple digital terrain units; Using a bilinear interpolation algorithm on irregularly distributed measuring points in the digital terrain model, a weighted average of the far first area terrain correction value and the far second area generalized terrain correction value of the regularly distributed measuring points near the irregularly distributed measuring points is calculated to obtain the far area terrain correction value of the irregularly distributed measuring points; The land gravity far zone terrain correction value is determined based on the far zone one terrain correction value of the regularly distributed measuring points and the far zone two generalized terrain correction value.

2. The method for correcting land gravity far-zone terrain based on a digital elevation model according to claim 1, characterized in that: Obtain digital elevation models and regularly distributed measurement points to construct a digital terrain model within the far-field terrain correction range, including: Obtaining a high-precision, high-spatial-resolution digital elevation model and the geographic coordinates of regularly distributed measurement points, as well as a far-zone terrain correction range radius, wherein the far-zone terrain correction range radius includes a far-zone 1 terrain correction radius and a far-zone 2 generalized terrain correction radius; Calculating a range of digital terrain data with high spatial resolution covering the regularly distributed measuring points, and calculating a range of digital terrain data with low spatial resolution covering the regularly distributed measuring points; The digital terrain units corresponding to the high spatial resolution and the digital terrain units corresponding to the low spatial resolution are extracted respectively.

3. The method for correcting land gravity far-zone terrain based on a digital elevation model according to claim 2, characterized in that: The high spatial resolution is 15 arc-sec, the low spatial resolution is 1 arc-min, the first digital elevation unit is a high spatial resolution terrain unit, and the second digital elevation unit is a low spatial resolution terrain unit.

4. The method for correcting land gravity far-zone terrain based on a digital elevation model according to claim 1, wherein: Extracting, from the digital terrain model, according to the position information of any measuring point among the regularly distributed measuring points, a first digital elevation unit satisfying a far first zone terrain correction range of the arbitrary measuring point and a second digital elevation unit satisfying a far second zone generalized terrain correction range, comprising: Obtaining the position information of the regularly distributed measuring points, the spatial resolution of the far first area terrain correction range, the inner radius and outer radius of the far first area terrain correction range, and the outer radius of the far second area generalized terrain correction range; Determine digital terrain units within the terrain correction range of the regularly distributed measuring points based on the information, wherein the digital terrain units include high spatial resolution terrain units covered by the annular terrain in the first far zone and low spatial resolution terrain units covered by the circular terrain in the second far zone; The terrain of the far zone 1 and the middle zone is a circular interface, the terrain of the far zone 1 and the far zone 2 is a circular interface, and the minimum digital terrain unit of the digital elevation model is divided by longitude and latitude lines. The high spatial resolution terrain unit covered by the circular interface between the far zone 1 and the middle zone is uniformly subdivided to obtain a first terrain subunit, and the far zone terrain correction calculation is performed on the first terrain subunit distributed in the far zone 1 range of the measuring point; The low spatial resolution terrain unit covered by the circular interface of the far zone 1 and the far zone 2 is subdivided to obtain the second terrain subunit, and the far zone 2 generalized terrain correction calculation is performed on the second terrain subunit distributed in the far zone 2 of the measuring point.

5. The method for correcting land gravity far-zone terrain based on a digital elevation model according to claim 4, characterized in that: The terrain of the far first zone and the middle zone is a circular interface, the terrain of the far first zone and the far second zone is a circular interface, the minimum digital terrain unit of the digital elevation model is divided by longitude and latitude lines, and the high spatial resolution terrain unit covered by the circular interface of the far first zone and the middle zone is uniformly subdivided to obtain a first terrain subunit, including: The high spatial resolution terrain unit covered by the circular interface between the far first area and the middle area is uniformly divided to obtain the first terrain subunit; The first terrain sub-unit that falls outside the inner radius of the circular terrain in the far zone is involved in the terrain correction calculation in the far zone.

6. The method for correcting land gravity far-zone terrain based on a digital elevation model according to claim 4, characterized in that: The low spatial resolution terrain unit covered by the circular interface of the far zone 1 and far zone 2 is subdivided to obtain the second terrain sub-unit, and the far zone 2 generalized terrain correction calculation is performed on the second terrain sub-unit distributed in the far zone 2 range of the measuring point, including: The low spatial resolution terrain units covered by the circular interface of the far zone 1 and the far zone 2 are divided to obtain low spatial resolution terrain sub-units; The low spatial resolution terrain sub-units falling outside the inner radius of the circular terrain in the Far Second Zone and the low spatial resolution terrain units within the outer radius of the circular terrain in the Far Second Zone are involved in the calculation of the generalized terrain correction in the Far Second Zone.

7. The method for correcting land gravity far-zone terrain based on a digital elevation model according to claim 6, characterized in that: The far-field terrain correction value is: the vertical first-order partial derivative value of the first gravitational potential generated when the elevation value of the digital terrain unit is greater than the surplus mass corresponding to the elevation value of the measuring point, or when the elevation value of the digital terrain unit is less than the deficit mass corresponding to the elevation value of the measuring point, wherein the measuring points include regularly distributed measuring points and irregularly distributed measuring points; The generalized terrain correction value of the far second zone is: the vertical first-order partial derivative value of the second gravitational potential generated by the excess mass or deficit mass between the reference ellipsoid corresponding to the digital terrain unit and the top surface of the digital terrain unit.

8. The method for correcting land gravity far-zone terrain based on a digital elevation model according to claim 7, characterized in that: Obtaining a calculation formula for the vertical first-order partial derivative of the gravitational potential of a preset single digital terrain unit taking into account the curvature of the earth's surface at the arbitrary measuring point, including: A pair of longitude planes, a pair of latitude planes and a pair of concentric spheres are preset, and the corresponding integral space is [r1, r2] × [θ1, θ2] × [λ1, λ2]. The terrain correction value calculation expression of a single terrain unit is: [r1, r2] represents the distance from the bottom and top surfaces of the digital terrain unit to the center of the earth, [θ1, θ2] represents a pair of latitude values ​​that limit the digital terrain unit, and [λ1, λ2] represents a pair of longitude values ​​that limit the digital terrain unit; Among them, λ represents the longitude of the measuring point, θ represents the latitude of the measuring point, and r represents the distance from the measuring point to the center of the earth. represents the first-order partial derivative of the gravitational potential, i = 1, 2, 3 represent the north direction, east direction and radial direction of the local rectangular coordinate system of the measuring point, t represents the integration node in the longitude direction, j represents the integration node in the latitude direction, k represents the integration node in the radial direction, Δλ represents the longitude dimension of the digital terrain unit, Δθ represents the latitude dimension of the digital terrain unit, Δr represents the radial dimension of the digital terrain unit, G represents the gravitational constant, ρ represents the density; (λ, θ, r) are the spherical coordinates of the measuring point, (r′ k ,θ′ j ,λ′ t ) is the spherical coordinate of the integration node, N λ 、N θ and N r They represent the maximum order used when applying the triple Gauss-Legendre integral algorithm, ω λt 、ω θj 、ω rk As weight, according to the vertical partial derivative V3 of the calculated gravitational potential, l is used to represent the Euclidean distance between the integration node and the measuring point, then I3D(r′ k ,θ′ j ,λ′ t ) is: Among them, I3D(r′ k ,θ′ j ,λ′ t ) represents the integral kernel function of the local rectangular coordinate system of the measuring point, r′ k represents the radial distance from the integration node to the center of the earth, θ′ j represents the latitude of the integration node, λ′ t Indicates the longitude of the integration node.

9. The method for correcting land gravity far-zone terrain based on a digital elevation model according to claim 1, wherein: The bilinear interpolation algorithm includes: Among them, x and y represent the easting and northing coordinates in the plane rectangular coordinate system, respectively; the geographic coordinates of the irregular measuring point are (λ, θ), and δg(x, y) is the far-zone correction value; g(x1, y1), g(x2, y1), g(x1, y2), and g(x2, y2) are the far-zone terrain correction values ​​of the four elevation nodes near the irregular measuring point in the plane rectangular coordinate system; in the spherical coordinate system, the correspondence between the longitude and latitude of the four elevation nodes and the plane coordinates is x1→λ1, x2→λ2, y1→θ1, and y2→θ2, respectively.

10. A land gravity far-zone terrain correction system based on digital elevation model, characterized in that: The method for correcting land gravity far-zone terrain based on a digital elevation model as claimed in any one of claims 1 to 9 comprises: Terrain data acquisition module, used to obtain digital elevation models and regularly distributed measurement points, and construct digital terrain models within the far-area terrain correction range; a digital elevation extraction module for extracting, from the digital terrain model, a first digital elevation unit satisfying a far first zone terrain correction range and a second digital elevation unit satisfying a far second zone generalized terrain correction range of any measuring point in the regularly distributed measuring points based on the position information of the measuring point; a terrain improvement value calculation module, configured to obtain a preset calculation formula for the vertical first-order partial derivative of the gravitational potential of a single digital terrain unit at the arbitrary measuring point, taking into account the curvature of the earth's surface; and, based on the first digital elevation unit and the second digital elevation unit, respectively calculate the vertical first-order partial derivatives of the gravitational potential of all digital terrain units in the far first zone and the far second zone of the arbitrary measuring point and cumulatively sum them to obtain a far first zone terrain correction value and a far second zone generalized terrain correction value for the arbitrary measuring point, wherein the digital terrain model includes multiple digital terrain units; a bilinear interpolation calculation module, configured to use a bilinear interpolation algorithm to calculate the weighted average of the far first zone terrain correction value and the far second zone generalized terrain correction value near the irregularly distributed measuring points in the digital terrain model to obtain the far zone terrain correction value of the irregularly distributed measuring points; The far-zone terrain correction value determination module is used to determine the far-zone terrain correction value of land gravity based on the far-zone one terrain correction value of the regularly distributed measuring points and the far-zone two generalized terrain correction value.

Images