A method for determining land area height anomaly

Abstract

The present application belongs to the field of geodesy, and particularly relates to a method for determining land area height anomaly. The present application calculates the regional vertical deflection by using gravity and terrain data, and then calculates and constructs the land area height anomaly model based on the vertical deflection, so that the calculated height anomaly can effectively absorb the high-frequency information of the terrain, and the grid resolution of the height anomaly model is not constrained by the grid resolution of the gravity data, thereby effectively improving the construction precision of the height anomaly model, i.e. improving the accuracy of the determined land area height anomaly. Furthermore, the present application performs plane approximation on the grid, establishes a local coordinate system with the current calculation point as the center, obtains the relationship between the residual vertical deflection and the height anomaly of the corresponding grid in the central area, and then combines the double quadratic polynomial interpolation function of the residual vertical deflection constructed for the central area to obtain the contribution of the residual vertical deflection of 3x3 grids to the height anomaly, so as to more accurately calculate the height anomaly value corresponding to the central area.

Description

Technical Field

[0001] This invention belongs to the field of geodesy, and specifically relates to a method for determining elevation anomalies in land areas. Background Technology

[0002] my country has one of the most complex gravity fields in the world, but currently there are only about 1 million gravity measurement points on land, mainly from geological departments. Fewer than 200,000 gravity points are used for geodetic surveying. The most significant problem is the highly uneven distribution of these points; they are denser in the eastern regions and sparser in the western and mountainous areas. Approximately 40% of the 5′×5′ grid lacks actual measured gravity points. Therefore, the construction of basic gravity anomaly grid data in my country, especially in the western and mountainous regions, primarily relies on high-resolution topographic data.

[0003] The common method for constructing land elevation anomaly models is as follows: first, based on Molodensky theory, gravity and topographic data are directly calculated, and then the systematic difference is fitted using GNSS / leveling. Since elevation anomalies cannot effectively absorb high-frequency information such as topography, it is difficult to significantly improve the accuracy. The method for determining elevation anomalies mainly relies on dense GNSS / leveling, which is not only costly, but also results in a low grid resolution of the elevation anomaly model due to the low resolution of the gravity data grid. Therefore, the accuracy of the determined elevation anomalies is relatively low. Summary of the Invention

[0004] The purpose of this invention is to provide a method for determining elevation anomalies in land areas, which solves the problem that existing methods for determining elevation anomalies cannot effectively absorb high-frequency information such as topography and the grid resolution of elevation anomaly models is low, resulting in relatively low accuracy of the determined elevation anomalies.

[0005] To achieve the above objectives, the present invention provides a method for determining elevation anomalies in land areas, which calculates the model elevation anomaly values ​​of points in each grid in the target area by referring to the Earth's gravity field model.

[0006] The remaining vertical deviation of the ground at each point in the grid within the target area is calculated using grid gravity anomaly and terrain data.

[0007] Based on the relationship between the residual vertical deviation and the geoid, the residual geoid value of each point in each grid is obtained by using the residual vertical deviation value of the ground at each point in the target area.

[0008] The beneficial effects of the above technical solution are as follows: by using gravity and terrain data to calculate the vertical deviation of the region, and using the vertical deviation to participate in the calculation and construction of the land area elevation anomaly model, the final calculated elevation anomaly can effectively absorb high-frequency information such as terrain, and the grid resolution of the elevation anomaly model is not constrained by the grid resolution of gravity data, thereby effectively improving the construction accuracy of the elevation anomaly model, that is, improving the accuracy of the determined land area elevation anomaly.

[0009] Furthermore, the method for obtaining the remaining elevation anomalies corresponding to points in each grid is as follows:

[0010] For each current calculation point, for the term corresponding to the central area in the formula relating the remaining vertical deviation and the elevation anomaly, the total value of the term corresponding to the central area is calculated using the relationship between the remaining vertical deviation and the elevation anomaly of the grid corresponding to the central area; the current calculation point is a point in one of the grids in the target area;

[0011] For the term corresponding to the non-central area in the formula relating the residual vertical deviation and the geoid eccentricity, substitute the ground residual vertical deviation value of each point in the target area into the term corresponding to the non-central area to obtain the value of each term corresponding to the non-central area.

[0012] Add the total value of the items corresponding to the central area to the values ​​of all items corresponding to the non-central areas to obtain the remaining elevation anomaly value corresponding to the current calculation point;

[0013] The central region refers to the set number of grid regions centered on the current calculation point and closest to it within the integration region corresponding to the current calculation point.

[0014] The beneficial effect of the above technical solution is that it can avoid the problem that when the integral flow point approaches the current calculation point, the relationship formula between the final residual vertical deviation and the elevation anomaly will become singular, causing the elevation anomaly to be uncalculated.

[0015] Furthermore, by analyzing the relationship between the remaining vertical deviation of the grid corresponding to the central area and the geoid undulation, the total value of the corresponding term for the central area is calculated as follows:

[0016] A biquadratic polynomial interpolation function for the remaining vertical deviation is constructed for the central region; based on the biquadratic polynomial interpolation function for the remaining vertical deviation, the total value of the terms corresponding to the central region is calculated.

[0017] Furthermore, the relationship between the remaining vertical deviation of the grid corresponding to the central area and the geoid undulation is as follows:

[0018]

[0019] In the formula, x and y are the coordinates of the approximate grid obtained by approximating the grid in the central area using a local coordinate system with the current calculation point as the center, where the x-axis points north and the y-axis points east. This represents the remaining vertical deviation of the grid. This represents the elevation anomaly corresponding to the central region.

[0020] Furthermore, the central region is the 3×3 grid closest to the current calculation point; the biquadratic polynomial interpolation function for the remaining perpendicular deviation constructed for the 3×3 grid is:

[0021]

[0022] The total value of the terms corresponding to the central region, obtained from the biquadratic polynomial interpolation function of the remaining vertical deviations of the 3×3 grids, is:

[0023]

[0024] in, , , , All are interpolation coefficients corresponding to the biquadratic polynomial interpolation function for the remaining perpendicular deviations of the 3×3 grid; and

[0025]

[0026] In the formula, , The step size of the planar approximate mesh is approximately taken as... , ; For the latitudinal resolution of the grid gravity anomaly, This represents the meridional resolution of the grid gravity anomaly.

[0027] Furthermore, the method for calculating the remaining vertical deviation of the ground at each point in the target area using grid gravity anomaly and terrain data is as follows:

[0028]

[0029] in, The contribution term to the residual gravity anomaly, its single-band or multi-band Fourier formula is as follows:

[0030]

[0031] In the formula, Here, k represents the geocentric latitude and longitude of the current calculation point, which is a point in the grid at row i and column j within the target region; k is the possible value of row number i in the target region. For the latitudinal resolution of the grid gravity anomaly, The meridional resolution of the grid gravity anomaly. To address the Faye gravity anomaly, , They represent the Discrete Fourier Transform and the Inverse Fourier Transform, respectively. This represents the average normal gravity value of the target area; The model gravity anomaly value for the gravity field model used is calculated using the following formula:

[0032]

[0033] in, Let R be the geocentric radius of the current calculation point, and R be the reference radius of the spherical harmonic expansion. Let N be the Earth's gravitational constant, and N be the highest order of the Earth's gravity field model used. and These are the spherical harmonic coefficients of the perturbation potential corresponding to the gravity field model. For the normalized associative Legendre function, n is the order in the spherical harmonic expansion, and m is the degree; The kernel function is calculated using the following formula:

[0034]

[0035] , respectively, represent the geocentric latitudes of the integration flow point and the current calculation point; S is the Stokes kernel function, which has a correlation coefficient with respect to variables. The expression for the derivative is:

[0036]

[0037] It is the azimuth angle of the integration flow point relative to the current calculation point, and its magnitude is calculated by the following formula:

[0038]

[0039] For elevation-related correction terms, the single-zone or multi-zone Fourier formula for this term is:

[0040]

[0041] in The average normal gravity value for the target area. Indicates the elevation of the integration flow point. This represents the radial gradient of the complete Bouguer anomaly at the integral flow point.

[0042] Furthermore, by referencing the Earth's gravity field model, the method for calculating the model elevation anomalies of points in each grid within the target region is as follows:

[0043] The model's geoid undulation is calculated using a first-order Taylor series expansion along the elevation direction, as shown in the following formula:

[0044]

[0045] In the formula, i and j are the row and column numbers of the grid where the current calculation point is located, respectively. Let be the first-order Taylor series term of the model geoid undulation at the i-th row and j-th column of the grid, expanded along the elevation direction. Let R be the geodetic height of a point in the grid, and R be the reference radius of the spherical harmonic expansion. The geocentric radius of the current calculation point. Let be the current calculation point, which is a point in the grid at row i and column j in the target region. These are the geocentric latitude and longitude of the current calculation point, respectively. Let be the geocentric radius vector of the i-th latitude circle on the ellipsoid;

[0046] To obtain the model elevation undulation up to the 0th order term, the calculation formula is:

[0047]

[0048] in, The average normal gravity value for the target area. R is the geocentric radius at the current calculation point, and R is the reference radius of the spherical harmonic expansion, the size of which is related to the gravity field model used. Let N be the Earth's gravitational constant, and N be the highest order of the Earth's gravity field model used. and These are the spherical harmonic coefficients of the perturbation potential corresponding to the gravity field model. For the normalized associative Legendre function, n is the order in the spherical harmonic expansion, and m is the degree;

[0049] .

[0050] The beneficial effects of the above technical solution are as follows: by performing a first-order Taylor series expansion of the model's geoid undulation calculation in the elevation direction, the grid geoid undulation values ​​of the Earth's surface in the determined area can be calculated quickly, thereby improving computational efficiency.

[0051] Furthermore, the formula relating the remaining vertical deviation to the geoid undulation is:

[0052]

[0053] In the formula, Let be the meridional and lateral components of the remaining perpendicular deviation corresponding to the grid in the i-th row and j-th column, respectively. Let i be the area element size of the grid in the i-th row and j-th column. Remaining elevation anomaly; Here, represents the geocentric latitude and longitude of the current calculation point, respectively; R is the reference radius of the spherical harmonic expansion. This is the azimuth angle corresponding to the vertical deviation; The integral kernel function is expressed as follows:

[0054]

[0055] in The angular distance between the center of the sphere and the current calculation point is defined as the angular distance between the integral flow point and the current calculation point. The integration region corresponding to the integral flow point is determined based on the geographical grid resolution of the vertical deviation data, combined with the minimum integration radius, with the current calculation point as the center. Attached Figure Description

[0056] Figure 1 This is a flowchart illustrating the method and verification process for determining elevation anomalies in a land area, as described in an embodiment of the method for determining elevation anomalies in a land area according to the present invention.

[0057] Figure 2 This is a schematic diagram of gravity anomaly in a 4°×6° experimental area in an embodiment of the method for determining elevation anomalies in land areas according to the present invention.

[0058] Figure 3 This is a schematic diagram of the elevation of a 4°×6° experimental area in an embodiment of the method for determining elevation anomalies in land areas according to the present invention.

[0059] Figure 4a This is a schematic diagram of the meridian component of the vertical deviation of a 4°×6° experimental area in an embodiment of the method for determining elevation anomalies in land areas according to the present invention.

[0060] Figure 4b This is a schematic diagram of the vertical deviation of the zonal component in a 4°×6° experimental area in an embodiment of the method for determining elevation anomalies in land areas according to the present invention.

[0061] Figure 5 In an embodiment of the method for determining elevation anomalies in land areas according to the present invention, the elevation anomaly is determined according to the method for determining elevation anomalies based on vertical deviation calculation of the present invention. A schematic diagram of elevation anomalies;

[0062] Figure 6 In an embodiment of the method for determining elevation anomalies in land areas according to the present invention, the elevation anomaly is determined according to the prior art method for determining elevation anomalies. A schematic diagram of elevation anomalies. Detailed Implementation

[0063] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0064] Implementation Examples of Methods for Determining Elevation Anomalies in Land Areas

[0065] This embodiment provides a technical solution for determining elevation anomalies in land areas, referring to... Figure 1 Mainly includes:

[0066] By referencing the Earth's gravity field model, the model elevation anomaly values ​​of points in each grid within the target area are calculated.

[0067] The remaining vertical deviation of the ground at each point in the grid within the target area is calculated using grid gravity anomaly and terrain data.

[0068] Based on the relationship between residual vertical deviation and geostationary anomaly, the residual geostationary anomaly value corresponding to each point in each grid is obtained by using the residual vertical deviation value of the ground at each point in the target area. To obtain more optimized calculation results, in this embodiment, "points in each grid" refers to the midpoint of the grid.

[0069] Specifically, the method for calculating the model elevation anomalies of points in each grid within the target area by referencing the Earth's gravity field model is as follows:

[0070] The model's geoid undulation is calculated using a first-order Taylor series expansion along the elevation direction, as shown in the following formula:

[0071]

[0072] In the formula, i and j are the row and column numbers of the grid where the current calculation point is located, respectively. Let be the first-order Taylor series term of the model geoid undulation at the i-th row and j-th column of the grid, expanded along the elevation direction. Let R be the geodetic height of a point in the grid, and R be the reference radius of the spherical harmonic expansion. The geocentric radius of the current calculation point. Let be the current calculation point, which is a point in the grid at row i and column j in the target region. These are the geocentric latitude and longitude of the current calculation point, respectively. Let be the geocentric radius of the i-th latitude circle on the ellipsoid; thus, the Legendre function only needs to be calculated once for the same latitude circle, thereby speeding up the calculation. The sine and cosine function values ​​of longitude are also obtained by recursion to speed up the calculation.

[0073] To obtain the model elevation undulation up to the 0th order term, the calculation formula is:

[0074]

[0075] in, The average normal gravity value for the target area. Let R be the geocentric radius of the current calculation point, and R be the reference radius of the spherical harmonic expansion. Let N be the Earth's gravitational constant, and N be the highest order of the Earth's gravity field model used. and These are the spherical harmonic coefficients of the perturbation potential corresponding to the gravity field model. For the normalized associative Legendre function, n is the order in the spherical harmonic expansion, and m is the degree;

[0076]

[0077] in, The reference radius for spherical harmonic expansion is usually taken as the major radius of the Earth's reference ellipsoid. Its size is related to the gravity field model used. For example, the EGM2008 gravity field model uses this parameter with a size of 6378136.3m.

[0078] In this embodiment, using grid gravity anomalies and terrain data, the two components of the vertical deviation are calculated according to the Molodensky theory that takes into account the first-order term of the terrain, thereby obtaining the ground residual vertical deviation value of each grid point in the target area. That is, the method of calculating the ground residual vertical deviation value of each grid point in the target area using grid gravity anomalies and terrain data is as follows:

[0079]

[0080] in, The contribution term to the residual gravity anomaly, its single-band or multi-band Fourier formula is as follows:

[0081]

[0082] In the formula, The coordinates are the geocentric latitude and longitude of the current calculation point, which is a point in the grid at row i and column j within the target area; k is the possible value of row i in the target area, i.e., an integer not exceeding the maximum number of rows in the target area; For the latitudinal resolution of the grid gravity anomaly, The meridional resolution of the grid gravity anomaly. To address the Faye gravity anomaly, , They represent the Discrete Fourier Transform and the Inverse Fourier Transform, respectively. This represents the average normal gravity value of the target area; The model gravity anomaly value for the gravity field model used is calculated using the following formula:

[0083]

[0084] in, Let R be the geocentric radius of the current calculation point, and R be the reference radius of the spherical harmonic expansion. Let N be the Earth's gravitational constant, and N be the highest order of the Earth's gravity field model used. and These are the spherical harmonic coefficients of the perturbation potential corresponding to the gravity field model. For the normalized associative Legendre function, n is the order in the spherical harmonic expansion, and m is the degree; The kernel function is calculated using the following formula:

[0085]

[0086] , respectively, represent the geocentric latitude of the integration flow point and the current calculation point; where the integration flow point represents the location of each vertical deviation data point that needs to be accumulated after discretization of the integration operation, and the range of the integration flow point is usually the range for calculating the remaining vertical deviation data; S is the Stokes kernel function, which affects the variables The expression for the derivative is:

[0087]

[0088] It is the azimuth angle of the integration flow point relative to the current calculation point, and its magnitude is calculated by the following formula:

[0089]

[0090] For elevation-related correction terms, the single-zone or multi-zone Fourier formula for this term is:

[0091]

[0092] in The average normal gravity value for the target area. Indicates the elevation of the integration flow point. This represents the radial gradient of the complete Bouguer anomaly at the integral flow point.

[0093] After obtaining the residual vertical deviation values ​​of the ground at each grid point in the target area, the influence of the residual vertical deviation on the geoid can be used to achieve a fast and high-precision geoid solution. Specifically, this embodiment uses the inverse Vening-Meinesz formula and the Bruns formula to obtain the relationship between the geoid and the vertical deviation of the grid points.

[0094]

[0095] In the formula, On a sphere of radius R The perpendicular deviation at point cistern and the components of the east-west axis. It is the azimuth angle corresponding to the vertical deviation, and the integral kernel function. The expression is as follows

[0096]

[0097] in This is the angular distance between the center of the sphere and the current calculation point.

[0098] It is worth noting that the above integral assumes the total mass of the normal gravitational field is equal to the actual total mass of the Earth, thus avoiding systematic errors caused by the inconsistency between the two. To avoid [further issues related to the azimuth angle]... Quadrant determination, rewriting the above formula, yields the final formula relating the remaining vertical deviation to the geoid undulation:

[0099]

[0100] In the formula, Let be the meridional and lateral components of the remaining perpendicular deviation corresponding to the grid in the i-th row and j-th column, respectively. Let i be the area element size of the grid in the i-th row and j-th column. Remaining elevation anomaly; Here, represents the geocentric latitude and longitude of the current calculation point, respectively; R is the reference radius of the spherical harmonic expansion. This is the azimuth angle corresponding to the vertical deviation; The integral kernel function is expressed as follows:

[0101]

[0102] in For the integral flow point With the current calculation point The angular distance between the centers of the spheres. The integration region corresponding to the integration flow point is centered on the calculation point and determined based on the geographic grid resolution of the vertical deviation data, combined with the minimum integration radius; for example, for residual vertical deviation grid data with a resolution of 2′, assuming an integration radius of 2°, the calculation point is determined. Then, expand by 2° in each of the four directions: east, west, south, and north. That is, in addition to the 2′ geographic grid where the calculation point is located, expand by 60 grids in each of the four directions. Therefore, i and j in the formula are taken as 1 to 121. If the calculation point is located in the 61st row and 61st column grid, it is necessary to use the remaining vertical deviation of the 121st row × 121st column grid to integrate and sum according to the above formula to calculate the remaining elevation anomaly at the calculation point located in the 61st row and 61st column grid.

[0103] However, when the integration point approaches the calculation point, the formula relating the final residual vertical deviation to the geoid becomes singular. Therefore, the calculation for the case where the integration point approaches the calculation point requires a different formula; that is, because...

[0104]

[0105] Therefore, when When the kernel function is very small, it can be approximated as:

[0106]

[0107] in Let be the Euclidean distance between the integration flow point and the calculation point. Therefore, for each calculation point, if the location of the integration flow point is too close to the calculation point, the denominator in the integral summation term will be approximately 0, leading to singular or large errors in the calculation. Therefore, in actual calculations, a region within a certain latitude and longitude range from the calculation point is usually selected as the central region. When using the formula relating the residual vertical deviation and the geoid, the formula is expanded. Terms not corresponding to the central region are substituted into the residual vertical deviation value already calculated in the final formula relating the residual vertical deviation and the geoid. Terms corresponding to the central region are calculated separately according to the modified formula for the central region. Finally, the terms corresponding to the non-central region and the terms corresponding to the central region are added together to obtain the geoid of a given calculation point.

[0108] In this embodiment, the central region refers to the grid region within the integration region corresponding to the current calculation point that is the closest grid region to the current calculation point. For example, when the resolution is 2′, the radius of the central region in this embodiment is 2′, which means taking one row of latitude grids in the north and south and one column of longitude grids in the east and west, that is, the central 3×3 grids are the central region.

[0109] Based on the relationship between residual vertical deviation and geoid eccentricity, the method for obtaining the residual geoid eccentricity value of each point in each grid within the target area by using the residual vertical deviation of the ground at each point in the grid is as follows:

[0110] For each current calculation point, for the term corresponding to the central area in the formula relating the remaining vertical deviation and the geoid undulation, the total value of the term corresponding to the central area is calculated by using the relationship between the remaining vertical deviation and the geoid undulation of the grid corresponding to the central area; here, the current calculation point refers to a point in a certain grid in the target area.

[0111] For the term corresponding to the non-central area in the formula relating the residual vertical deviation and the geoid eclipse, substitute the ground residual vertical deviation value of each point in the grid of the target area into the term corresponding to the non-central area to obtain the value of each term corresponding to the non-central area.

[0112] The total value of the items corresponding to the central area is added together with the values ​​of all items corresponding to the non-central areas to obtain the remaining elevation anomaly value corresponding to the current calculation point.

[0113] To more accurately calculate the elevation anomaly corresponding to the central area, this embodiment performs the following planar approximation on the grid; specifically, a local coordinate system is established with the current calculation point as the center, the x-axis pointing north and the y-axis pointing east, with...

[0114]

[0115] Under the plane approximation condition, the relationship between the residual vertical line deviation of the grid corresponding to the central area and the geoid undulation is as follows:

[0116]

[0117] In the formula, x and y are the coordinates of the planar approximation grid obtained after approximating the grid in the central area by a local coordinate system with the calculation point as the center, where the x-axis points north and the y-axis points east. This represents the remaining vertical deviation of the grid. This represents the elevation anomaly corresponding to the central region.

[0118] The total value of the corresponding term in the central area is calculated by relating the residual vertical deviation of the grid to the geoid undulation as follows:

[0119] Construct a biquadratic polynomial interpolation function for the residual vertical deviation in the central region; calculate the total value of the terms for the corresponding central region based on the biquadratic polynomial interpolation function for the residual vertical deviation.

[0120] In this embodiment, the central region is the 3×3 grid closest to the calculation point; the biquadratic polynomial interpolation function for the remaining vertical deviation constructed for the 3×3 grid in the central region is:

[0121]

[0122] The interpolation coefficients for the remaining vertical deviations of the 3×3 grid can be obtained. Thus, the contribution of these remaining vertical deviations of the 3×3 grid to the geoid undulation can be expressed as:

[0123]

[0124] Taking advantage of the symmetry of the integration region, the integral result of the terms with odd powers of x and y is zero, and finally only... , , , The relevant terms are retained and obtained through integration:

[0125]

[0126] in, , , , All are interpolation coefficients corresponding to the biquadratic polynomial interpolation function for the remaining perpendicular deviations of the 3×3 grid; and

[0127]

[0128] In the formula, , The step size of the planar approximate mesh is approximately taken as... , ; For the latitudinal resolution of the grid gravity anomaly, This represents the meridional resolution of the grid gravity anomaly.

[0129] The contribution of the remaining vertical deviation of the 3×3 grid to the geoid undulation. This is the total value of the terms in the corresponding central region obtained from the biquadratic polynomial interpolation function based on the remaining vertical deviations of the 3×3 grids.

[0130] Below, this embodiment takes the construction of a quasi-geoid as an example. Using gravity anomaly data and topographic data of the same grid in a certain area, the quasi-geoid is calculated using the traditional Molodensky theory and the new method proposed in this patent, respectively. Then, the accuracy of the quasi-geoid constructed by the two methods is evaluated using high-precision GNSS / leveling data, thereby illustrating the advantages and innovation of the method proposed in this patent.

[0131] This embodiment selects a 4°×6° experimental area. The gravity anomaly and elevation of this area are as follows: Figure 2 , Figure 3 As shown, the two components of the vertical deviation are calculated according to the Molodensky theory, which takes into account the first-order terms of the terrain. The magnitudes of the numerical changes are as follows: Figure 4a , 4b As shown, Figure 4a shows the meridional component, Figure 4b The components are Mao and You; the data resolution is... .

[0132] Using known Figure 4a , Figure 4b Based on the vertical deviation data within a range of 4° × 6°, and following the elevation anomaly determination method in this embodiment, using the inverse Vening-Meinesz formula and the removal recovery technique (which divides the calculation of elevation anomalies into two parts: elevation anomaly model values ​​and residual elevation anomalies), a set of elevation anomalies with a shrinkage of 1.5° (i.e., 1° × 3°) is obtained; then, based on the vertical deviation calculation... Elevation anomalies such as Figure 5 As shown.

[0133] Similarly, directly using known information Figure 2 , Figure 3 Gravity anomaly and topographic data within a 4°×6° area were used to calculate the experimental area using the Molodensky theory and the removal recovery technique, with an integration radius of 1.5°. The results of the elevation anomaly are shown below. Figure 6 .

[0134] Depend on Figure 5 and Figure 6 A comparison reveals that the overall trend of the elevation anomaly calculated based on vertical deviation in this embodiment is similar to that of the elevation anomaly calculated based on gravity anomaly in the prior art. The elevation anomalies calculated using both methods were interpolated to the elevation anomalies of several known GNSS / leveling points, and the difference was calculated between these interpolations and the known high-precision GNSS / leveling elevation anomalies (considered the true values ​​here). The statistical analysis results are shown in Table 1.

[0135] Table 1. Data analysis of elevation anomalies and true values ​​using different algorithms / m

[0136]

[0137] As can be clearly seen from Table 1, the elevation anomaly obtained based on vertical deviation is about 5 cm better than the elevation anomaly obtained based on gravity anomaly in the traditional method. Moreover, the average and maximum values ​​of the elevation anomaly obtained by the former are smaller than the corresponding average and maximum values ​​of the latter. This proves that the accuracy of the method for determining the elevation anomaly of land areas in this embodiment is significantly better than that of the prior art.

[0138] This invention has the following characteristics:

[0139] 1) Calculate the vertical deviation of the region using gravity and topographic data, and then construct a land area elevation anomaly model based on the vertical deviation. This allows the calculated elevation anomaly to effectively absorb high-frequency information such as topography, and makes the grid resolution of the elevation anomaly model not constrained by the grid resolution of gravity data, thereby effectively improving the construction accuracy of the elevation anomaly model, that is, improving the accuracy of the determined land area elevation anomaly.

[0140] 2) The grid was approximated in a plane. A local coordinate system was established with the current calculation point as the center, with the x-axis pointing north and the y-axis pointing east. The relationship between the residual vertical deviation of the grid corresponding to the central area and the geoid was obtained. Then, combined with the biquadratic polynomial interpolation function of the residual vertical deviation constructed for the central area, the contribution of the residual vertical deviation of the 3×3 grids to the geoid was obtained, so as to calculate the geoid value corresponding to the central area more accurately.

[0141] 3) The calculation of the model's geoid undulation is performed by first-order Taylor series expansion in the elevation direction to quickly calculate the geoid undulation value of the grid model on the Earth's surface in the determined area, thereby improving computational efficiency.

[0142] It should be understood that the above-described specific embodiments of the present invention are merely illustrative or explanatory of the principles of the present invention, and do not constitute a limitation thereof.

Claims

1. A method for determining elevation anomalies in a land area, characterized in that, By referencing the Earth's gravity field model, the model elevation anomaly values ​​of points in each grid within the target area are calculated. The remaining vertical deviation of the ground at each point in the grid within the target area is calculated using grid gravity anomaly and terrain data. For each current calculation point, for the term corresponding to the central region in the formula relating the remaining vertical deviation and the geoid, a biquadratic polynomial interpolation function for the remaining vertical deviation is constructed for the central region; based on the biquadratic polynomial interpolation function for the remaining vertical deviation, the total value of the term corresponding to the central region is calculated. The current calculation point is a point within a grid cell of the target region; For the term corresponding to the non-central area in the formula relating the residual vertical deviation and the geoid, substitute the ground residual vertical deviation value of each point in the target area into the term corresponding to the non-central area to obtain the value of each term corresponding to the non-central area; add the total value of the term corresponding to the central area and the values ​​of each term corresponding to the non-central area to obtain the residual geoid value corresponding to the current calculation point. The central region refers to the set number of grid areas centered on the current calculation point and closest to it within the integration region corresponding to that calculation point.

2. The method for determining elevation anomalies in land areas according to claim 1, characterized in that, The relationship between the remaining vertical deviation of the grid in the central area and the geoid is as follows: In the formula, x and y are the coordinates of the approximate grid obtained by approximating the grid in the central area using a local coordinate system with the current calculation point as the center, where the x-axis points north and the y-axis points east. This represents the remaining vertical deviation of the grid. This represents the elevation anomaly corresponding to the central region.

3. The method for determining elevation anomalies in land areas according to claim 2, characterized in that, The central region is the 3×3 grid closest to the current calculation point; the biquadratic polynomial interpolation function for the remaining vertical deviation constructed for the 3×3 grid is: The total value of the terms corresponding to the central region, obtained from the biquadratic polynomial interpolation function of the remaining vertical deviations of the 3×3 grids, is: in, , , , All are interpolation coefficients corresponding to the biquadratic polynomial interpolation function for the remaining perpendicular deviations of the 3×3 grid; and In the formula, , The step size of the planar approximate mesh is approximately taken as... , ; For the latitudinal resolution of the grid gravity anomaly, This represents the meridional resolution of the grid gravity anomaly.

4. The method for determining elevation anomalies in land areas according to any one of claims 1-3, characterized in that, The method for calculating the remaining vertical deviation of the ground at each point in the target area using grid gravity anomaly and terrain data is as follows: in, The contribution term to the residual gravity anomaly, its single-band or multi-band Fourier formula is as follows: In the formula, Here, k represents the geocentric latitude and longitude of the current calculation point, which is a point in the grid at row i and column j within the target region; k is a possible value for row number i in the target region. For the latitudinal resolution of the grid gravity anomaly, The meridional resolution of the grid gravity anomaly. To address the Faye gravity anomaly, , They represent the Discrete Fourier Transform and the Inverse Fourier Transform, respectively. This represents the average normal gravity value of the target area; The model gravity anomaly value for the gravity field model used is calculated using the following formula: in, Let R be the geocentric radius of the current calculation point, and R be the reference radius of the spherical harmonic expansion. Let N be the Earth's gravitational constant, and N be the highest order of the Earth's gravity field model used. and These are the spherical harmonic coefficients of the perturbation potential corresponding to the gravity field model. For the normalized associative Legendre function, n is the order in the spherical harmonic expansion, and m is the degree; The kernel function is calculated using the following formula: , respectively, represent the geocentric latitudes of the integration flow point and the current calculation point; S is the Stokes kernel function, which has a correlation coefficient with respect to variables. The expression for the derivative is: It is the azimuth angle of the integration flow point relative to the current calculation point, and its magnitude is calculated by the following formula: For elevation-related correction terms, the single-zone or multi-zone Fourier formula for this term is: in The average normal gravity value for the target area. Indicates the elevation of the integration flow point. This represents the radial gradient of the complete Bouguer anomaly at the integral flow point.

5. The method for determining elevation anomalies in land areas according to any one of claims 1-3, characterized in that, The method for calculating the model elevation anomalies of points in each grid within the target area by referencing the Earth's gravity field model is as follows: The model's geoid undulation is calculated using a first-order Taylor series expansion along the elevation direction, as shown in the following formula: In the formula, i and j are the row and column numbers of the grid where the current calculation point is located, respectively. Let be the first-order Taylor series term of the model geoid undulation at the i-th row and j-th column of the grid, expanded along the elevation direction. Let R be the geodetic height of a point in the grid, and R be the reference radius of the spherical harmonic expansion. The geocentric radius of the current calculation point. Let be the current calculation point, which is a point in the grid at row i and column j in the target region. These are the geocentric latitude and longitude of the current calculation point, respectively. Let be the geocentric radius vector of the i-th latitude circle on the ellipsoid; To obtain the model elevation undulation up to the 0th order term, the calculation formula is: in, The average normal gravity value for the target area. R is the geocentric radius at the current calculation point, and R is the reference radius of the spherical harmonic expansion, the size of which is related to the gravity field model used. Let N be the Earth's gravitational constant, and N be the highest order of the Earth's gravity field model used. and These are the spherical harmonic coefficients of the perturbation potential corresponding to the gravity field model. For the normalized associative Legendre function, n is the order in the spherical harmonic expansion, and m is the degree; 。 6. The method for determining elevation anomalies in a land area according to any one of claims 1-3, characterized in that, The formula relating the residual vertical deviation to the geoid is: In the formula, Let be the meridional and lateral components of the remaining perpendicular deviation corresponding to the grid in the i-th row and j-th column, respectively. Let be the area element size of the grid in the i-th row and j-th column, and , Remaining elevation anomaly; Here, represents the geocentric latitude and longitude of the current calculation point, respectively; R is the reference radius of the spherical harmonic expansion. This is the azimuth angle corresponding to the vertical deviation; The integral kernel function is expressed as follows: in For the integral flow point With the current calculation point The angular distance between the centers of the spheres; the integration region corresponding to the integration flow point is determined with the current calculation point as the center, based on the geographical grid resolution of the vertical deviation data, combined with the minimum integration radius.

Images