Elevation datum plane construction method and system based on low-order cosine correction Stokes kernel function, medium and product

Abstract

The invention discloses an elevation datum plane construction method and system based on a low-order cosine correction Stokes kernel function, a medium and a product, and relates to the technical field of physical geodetic survey. The method comprises the following steps: calculating by utilizing an earth gravitational field model, topographic data and a high-frequency topographic influence combined algorithm to obtain model data, direct influence of topographic compression on ground gravity anomaly and indirect influence of topographic compression on an elevation datum plane; determining an empirical value and a test value of a correction order of the low-order cosine correction Stokes kernel function based on the accumulative error of the model gravity anomaly and the test frequency bandwidth; calculating an empirical elevation datum plane and a test elevation datum plane based on the empirical value and the test value of the correction order by adopting a low-order cosine correction Stokes kernel function; and determining a final elevation datum based on the empirical elevation datum and the test elevation datum. The construction precision of the elevation datum plane can be improved.

Description

Technical Field

[0001] The present invention relates to the field of physical geodetic technology, and in particular to a method, system, medium and product for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function. Background Art

[0002] The two elevation systems, orthometric height and normal height, are widely used in different countries and regions. Their elevation datums are the geoid and quasi-geoid, respectively. Determining a geoid with 1 cm accuracy is an important goal of geodetic work in the 21st century. Over the past few decades, many geodetic researchers have improved the accuracy of (quasi-)geoids by continuously optimizing the geoid construction method. The third boundary value problem in the field of physical geodesy is the theoretical basis for the construction of regional (quasi-)geoids. The integral used in the third boundary value theory is the Stokes integral, and the physical input of the Stokes integral is the gravity anomaly data. Gravity anomaly is the difference between the actual gravity measurement value and the normal gravity value caused by the irregular distribution of the earth's mass. It has important value in geodesy and geophysical research and application.

[0003] In theory, the Stokes integral is a global integral. However, since global gravity data is difficult to obtain and the global integral calculation is large, the actual calculation usually uses the removal and recovery technique that only requires local integration. In the removal and recovery technique mode, the gravity field model is used as the reference gravity field to participate in the boundary value solution, so that the Stokes integral actually includes two types of gravity data: measured gravity anomalies and reference model gravity anomalies. The actual gravity measurement error and the gridding error of discrete measurement points make the long-wave accuracy of the measured gravity data relatively low, while the accuracy of the gravity field model in the low-frequency band is very high. With the increase of the spherical harmonic order, the accuracy of the reference model gravity data gradually decreases, while the relative accuracy of the measured gravity data gradually increases. The ideal Stokes kernel function should adjust the weights (or contribution rates) of the two types of gravity data in the generated elevation anomaly / geoid fluctuation according to the relative accuracy of the two types of gravity data. When using the standard Stokes kernel function, the conversion rate of measured data in the low-frequency band is high and the conversion rate of reference data is low, which makes the calculated low-order elevation anomaly / geoid fluctuation error large; the far-field truncation error of gravity data measured in the high-frequency band is large, which makes the calculated high-order elevation anomaly / geoid fluctuation error also large. Therefore, some geodetic scholars have carried out research on the optimization method of the Stokes kernel function. Correcting the Stokes kernel function can effectively improve the solution accuracy of the (quasi) geoid in the removal recovery mode, but there are relatively few research results on the correction method of the Stokes kernel function and the analysis of spectral characteristics. Summary of the invention

[0004] The purpose of the present invention is to provide a method, system, medium and product for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function, which can improve the construction accuracy of a geoid and a quasi-geoid.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A method for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function includes:

[0007] Obtaining the earth's gravity field model, ground gravity anomaly data, terrain data, GNSS level point data, the type of elevation reference surface and the test frequency bandwidth; the types of elevation reference surfaces include geoid and quasi-geoid; the geoid is the reference surface of the orthometric elevation system, and the quasi-geoid is the reference surface of the normal height elevation system;

[0008] The model data, the direct effect of terrain compression on ground gravity anomaly and the indirect effect of terrain compression on elevation datum are calculated by using the earth gravity field model, terrain data and high-frequency terrain influence combined algorithm; the model data includes one of the model geoid height and the model elevation anomaly and the model gravity anomaly; the indirect effect of terrain compression on elevation datum is the indirect effect of terrain compression on geoid or the indirect effect of terrain compression on quasi-geoid;

[0009] Determine the empirical value and test value of the correction order of the low-order cosine-corrected Stokes kernel function based on the accumulated error of the model gravity anomaly and the test frequency bandwidth; the correction order includes an initial correction order and a cutoff correction order;

[0010] The empirical elevation datum and the test elevation datum are calculated by using a low-order cosine-corrected Stokes kernel function based on ground gravity anomaly data, terrain data, model data, the direct effect of terrain compression on ground gravity anomaly, the indirect effect of terrain compression on elevation datum, the empirical value of the correction order and the test value of the correction order; the empirical elevation datum is an empirical geoid or an empirical quasi-geoid; the test elevation datum is a test geoid or a test quasi-geoid;

[0011] The final elevation datum is determined based on GNSS level point data, the empirical elevation datum and the test elevation datum; the final elevation datum is the final geoid or the final quasi-geoid.

[0012] A computer system comprises: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of any one of the above-mentioned methods for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function.

[0013] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of any of the above-mentioned methods for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function.

[0014] A computer program product comprises a computer program, which, when executed by a processor, implements the steps of any of the above-mentioned methods for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function.

[0015] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0016] The invention discloses a method, system, medium and product for constructing an elevation datum based on a low-order cosine corrected Stokes kernel function. The method comprises the following steps: obtaining an earth gravity field model, ground gravity anomaly data, terrain data, GNSS level point data, the type of elevation datum and a test frequency bandwidth; calculating model data, direct influence of terrain compression on ground gravity anomaly and indirect influence of terrain compression on elevation datum by using an earth gravity field model, terrain data and a high-frequency terrain influence combined algorithm; determining an empirical value of a correction order of a low-order cosine corrected Stokes kernel function and a test value of the correction order based on the accumulated error of the model gravity anomaly and the test frequency bandwidth; calculating an empirical elevation datum and a test elevation datum based on ground gravity anomaly data, terrain data, model data, direct influence of terrain compression on ground gravity anomaly, indirect influence of terrain compression on elevation datum, empirical value of correction order and test value of correction order by using a low-order cosine corrected Stokes kernel function; and determining a final elevation datum based on the GNSS level point data, the empirical elevation datum and the test elevation datum. The present invention can improve the construction accuracy of the elevation reference surface. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0018] Figure 1 Schematic diagram of the normalized far-zone cutoff coefficient of the standard Stokes kernel function;

[0019] Figure 2 Schematic diagram of the conversion rate of gravity data when using the standard Stokes kernel function;

[0020] Figure 3 Schematic diagram of the normalized far-zone cutoff coefficient of the low-order cosine modified Stokes kernel function;

[0021] Figure 4 Schematic diagram of the conversion rate of gravity data when using low-order cosine-corrected Stokes kernel function;

[0022] Figure 5 Schematic diagram of the process of constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function. DETAILED DESCRIPTION

[0023] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0024] The purpose of the present invention is to provide a method, system, medium and product for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function, aiming to improve the construction accuracy of a geoid and a quasi-geoid.

[0025] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0026] Example 1

[0027] like Figure 5 As shown, a method for constructing an elevation reference surface based on a low-order cosine modified Stokes kernel function in this embodiment includes:

[0028] Step 101: Obtain the earth gravity field model, ground gravity anomaly data, terrain data, GNSS level point data, the type of elevation reference surface and the test frequency bandwidth. The types of elevation reference surfaces include geoid and quasi-geoid. The geoid is the reference surface of the orthometric elevation system. The quasi-geoid is the reference surface of the normal elevation system.

[0029] Step 102: Calculate the model data, the direct effect of terrain compression on ground gravity anomaly and the indirect effect of terrain compression on elevation datum using the earth gravity field model, terrain data and high-frequency terrain influence combined algorithm. The model data includes one of the model geoid height and the model elevation anomaly and the model gravity anomaly. The indirect effect of terrain compression on elevation datum is the indirect effect of terrain compression on geoid or the indirect effect of terrain compression on quasi-geoid.

[0030] If the type of the elevation datum is the geoid, the model geoid height and the ground model gravity anomaly are calculated using the earth gravity field model and terrain data.

[0031] A high-frequency terrain influence combined algorithm is used to calculate the direct effect of terrain compression on ground gravity anomaly and the indirect effect of terrain compression on the geoid based on terrain data.

[0032] If the type of the elevation datum is a quasi-geoid, the model elevation anomaly and the ground model gravity anomaly are calculated using the earth gravity field model and terrain data.

[0033] A high-frequency terrain influence combined algorithm is used to calculate the direct effect of terrain compression on ground gravity anomaly and the indirect effect of terrain compression on the quasi-geoid based on terrain data.

[0034] Step 103: Determine the empirical value and test value of the correction order of the low-order cosine-corrected Stokes kernel function based on the cumulative error of the model gravity anomaly and the test frequency bandwidth. The correction order includes an initial correction order and a cutoff correction order.

[0035] The cumulative errors of each order of model gravity anomaly are calculated based on the calculation formula of the cumulative error of model gravity anomaly.

[0036] The calculation formula of the cumulative error of the model gravity anomaly is:

[0037]

[0038] Where L is the cutoff order of the reference gravity field model, δΔg G (L) represents the cumulative error of the model gravity anomaly truncated to the Lth order, represents the gravity anomaly error of the n-order model, γ represents normal gravity, δC nm Represents the spherical harmonic coefficient term C of the gravity field model nm The error, δS nm Represents the spherical harmonic coefficient term S of the gravity field model nm The error, σ represents the standard deviation, n represents the order, and m represents the degree.

[0039] The cumulative errors of each order of model gravity anomaly are calculated based on the calculation formula of the cumulative error of model gravity anomaly.

[0040] The lowest order corresponding to the cumulative error of the model gravity anomaly exceeding 0.01 mGal is set as the empirical value of the initial correction order.

[0041] The lowest order corresponding to the cumulative error of the model gravity anomaly exceeding 1 mGal is set as the empirical value of the cutoff correction order.

[0042] The empirical value of the initial correction order is added to or subtracted from a multiple of 10 or 100 within the frequency band of the test frequency bandwidth and set as the test value of the initial correction order.

[0043] The empirical value of the cutoff correction order is added to or subtracted from a multiple of 10 or 100 within the frequency band of the test frequency bandwidth and set as the test value of the cutoff correction order.

[0044] Step 104: Calculate the empirical elevation datum and the test elevation datum using a low-order cosine-corrected Stokes kernel function based on ground gravity anomaly data, terrain data, model data, the direct effect of terrain compression on ground gravity anomaly, the indirect effect of terrain compression on elevation datum, the empirical value of the correction order, and the test value of the correction order. The empirical elevation datum is an empirical geoid or an empirical quasi-geoid. The test elevation datum is a test geoid or a test quasi-geoid.

[0045] The low-order cosine modified Stokes kernel function includes two forms: non-extended and extended, specifically:

[0046]

[0047] Where S(ψ) represents the non-extended low-order cosine modified Stokes kernel function; S(r P ,ψ) represents the extended low-order cosine modified Stokes kernel function; rp represents the geocentric radius of the calculation point; ψ represents the angular distance between the calculation point and the flow integration point; represents the n-order spectral weight of the low-order cosine modified Stokes kernel function; n represents the order; R represents the radius of the earth; P n (cosψ) represents the nth-order Legendre polynomial; L 0 is the initial correction order, L 1 is the cutoff correction order.

[0048] The ground gravity anomaly is subtracted from the ground model gravity anomaly to obtain the gravity anomaly difference.

[0049] The gravity anomaly difference is added to the direct effect of the terrain compression on the ground gravity anomaly to obtain the residual ground Helmert gravity anomaly.

[0050] The remaining ground Helmert gravity anomaly is extended to obtain the processed remaining boundary surface Helmert gravity anomaly.

[0051] If the type of the elevation datum is the geoid, a non-extended low-order cosine-corrected Stokes kernel function is used to perform Stokes integral processing on the processed residual boundary surface Helmert gravity anomaly based on the empirical value of the correction order and the test value of the correction order to obtain the empirical residual geoid height and the test residual geoid height.

[0052] The empirical residual geoid height, the indirect effect of the terrain compression on the geoid and the model geoid height are added together to obtain the empirical geoid.

[0053] The test residual geoid height, the indirect effect of the terrain compression on the geoid and the model geoid height are added together to obtain the test geoid.

[0054] If the type of the elevation datum is a quasi-geoid, an extended form of a low-order cosine-corrected Stokes kernel function is used to perform Stokes integral processing on the processed residual boundary surface Helmert gravity anomaly based on the empirical value of the correction order and the test value of the correction order to obtain the empirical residual elevation anomaly and the test residual elevation anomaly.

[0055] The empirical residual elevation anomaly, the indirect effect of the terrain compression on the quasi-geoid and the model elevation anomaly are added together to obtain the empirical quasi-geoid.

[0056] The test residual elevation anomaly, the indirect effect of the terrain compression on the quasi-geoid and the model elevation anomaly are added together to obtain the test quasi-geoid.

[0057] Step 105: Determine a final elevation reference surface based on the empirical elevation reference surface and the test elevation reference surface; the final elevation reference surface is a final geoid or a final quasi-geoid.

[0058] If the type of the elevation datum is the geoid, the final elevation datum is determined based on the empirical geoid and the test geoid using GNSS level point data.

[0059] If the type of the elevation datum is a quasi-geoid, the final elevation datum is determined based on the empirical quasi-geoid and the test quasi-geoid using GNSS level point data.

[0060] Standard Stokes kernel function

[0061] Stokes integral is a theoretical tool for calculating geoid height or elevation anomaly using gravity anomaly data in the third boundary value problem. Stokes integral includes two forms: non-extended and extended. The non-extended Stokes integral is a theoretical tool for calculating geoid height, while the extended Stokes integral is a theoretical tool for calculating elevation anomaly. The standard Stokes kernel function in the extended form is

[0062]

[0063] Where S ST (r P ,ψ) is the extended form of the Stokes standard kernel function, R is the average radius of the earth, r P is the radius from the center of the earth at the calculation point in spherical coordinates, P n (cosψ) is the nth-order Legendre polynomial, ψ is the angular distance between the flow point and the calculation point, l P0 To calculate the distance from the point to the center of the sphere flow point, the formula is:

[0064] The non-expanded form of the standard Stokes kernel function is

[0065]

[0066] Where S ST (ψ) is the standard Stokes kernel function in its non-expanded form, and its symbol t is expressed as sin(ψ / 2). ST (ψ) is S ST (r P ,ψ) in the parameter r P This is a special case when R is taken, so the present invention only describes the correction method of the extended kernel function.

[0067] The Stokes kernel function can be expressed in the following general spectral form.

[0068]

[0069] Where S(r P ,ψ) represents the Stokes kernel function, ω n represents the n-order spectral weight (or correction factor). The difference between different types of Stokes kernel functions lies in the different values ​​of spectral weight. P ,ψ) specifically represents the standard Stokes kernel function, the spectral weights of all orders in equation (10) are all 1. P ,ψ) indicates that when the Stokes kernel function is modified, the value range of the spectral weight is usually [0,1]. The calculation formula of the spectral weight of the cosine low-order modified Stokes kernel function is:

[0070]

[0071] in, represents the n-order spectral weight of the low-order cosine modified Stokes kernel function; n represents the order; L 0 is the low-order initial correction order, L 1 is the low-order cutoff correction order. The cosine correction kernel function includes two forms: low-order correction only and high- and low-order correction. The spectral weight of the cosine high- and low-order corrected Stokes kernel function is:

[0072]

[0073] In the formula represents the n-order spectral weight of the cosine high and low order modified Stokes kernel function, L 2 is the high-order initial correction order, L 3 is the high-order cutoff correction order.

[0074] Introducing a piecewise function It is similar to the Stokes kernel function S(r P ,ψ) is

[0075]

[0076] Piecewise Function The Legendre series form of

[0077]

[0078] Where Q n is the n-th order far-zone truncation error coefficient of the Stokes kernel function, and its calculation formula is

[0079]

[0080] According to the expansion theory of spherical functions, we have

[0081]

[0082] Where Δg n is the nth-order gravity anomaly on the spherical surface, and Δg is the gravity anomaly on the spherical surface. The Stokes integral in the near zone can be expressed as the following general spectral form

[0083]

[0084] Where ζ is the elevation anomaly, which is the vertical distance between the quasi-geoid and the reference ellipsoid. The elevation anomaly data grid is usually used to represent the regional quasi-geoid. 1 It means that the calculation point is the center, ψ 1is the radius of the near zone, γ is the normal gravity of the ground calculation point, is the normalized n-th order far zone truncation error coefficient, which is calculated as follows:

[0085]

[0086] The spectral relationship between elevation anomaly and gravity anomaly is:

[0087]

[0088] Where ζ n It is the n-th order elevation anomaly.

[0089] Using equations (15) and (17), the general spectral decomposition formula of the Stokes integral under the recovery mode can be deduced as

[0090]

[0091] Where Δg T Measure gravity anomalies on a spherical surface, is the nth-order gravity anomaly, Δg G is the spherical reference gravity anomaly value, is its nth-order gravity anomaly, ζ G is the reference elevation anomaly, M is the highest order of measured gravity anomaly (usually determined by the resolution of the ground gravity grid according to the Nyquist sampling theorem), and L is the truncation order of the reference gravity field model. 2 ≥M,L 1 ≤L, and the spectral weights above the Lth order are all 1. Since both the reference and measured gravity data actually contain errors, the formula (18) It represents the conversion rate of n-order measured gravity data and its error into elevation anomaly. It represents the conversion rate of n-order reference gravity data and its error into elevation anomaly. For different types of Stokes kernel functions, it is only necessary to replace the normalized far-zone cutoff coefficient in equation (18) with and spectral weight factor ω n Just adjust it to the coefficients and spectral weights of the corresponding Stokes kernel function.

[0092] In the removal and recovery mode, the kernel function affects the accuracy of low- and high-frequency elevation anomalies differently. In the low-frequency band (2≤n<L) of removal and recovery, both the measured and reference gravity data contribute to the elevation anomaly. The sum of the conversion rates of each order of the low-frequency band measurement and reference gravity data is always 1, indicating that there is no far-field truncation error in the low-frequency band. Since the errors of the two gravity data are also converted into the elevation anomaly with the Stokes integral, the error of the low-frequency elevation anomaly comes from the low-frequency error of the measured and reference gravity data. Different Stokes kernel functions have different far-field truncation error coefficients, so when using different kernel functions, the proportion of the errors of the two gravity data in the conversion into low-order elevation anomaly errors is also different. In the high-frequency band (L+1~M) of removal and recovery, the conversion rate of the measured data is Reference data does not participate in the conversion. Since the high-order elevation anomaly is only contributed by the measured data and the high-order conversion rate is not 1, the error sources of the high-order elevation anomaly are the high-frequency error of the measured gravity data and the far-field truncation error.

[0093] In formula (18), r P Replace with R, γ replaced with the normal gravity γ of the boundary surface 0 , the elevation anomaly ζ is replaced by the geoid height N, and the corresponding spectral decomposition formula for calculating the Stokes integral of the geoid height in the removal recovery mode can be obtained.

[0094] The following is an analysis of the spectral characteristics of the Stokes integral using the standard and cosine modified Stokes kernel functions based on the spectral decomposition of the Stokes integral in the removal recovery mode. In the spectral characteristic analysis of the present invention, the calculation point elevation is 1000m, the integral radius is 1°, and the highest order of the reference model is 2160.

[0095] Depend on Figure 1 It can be seen that the normalized far-zone cutoff coefficient of the standard kernel function has a relatively large change amplitude. Among them, the maximum amplitude of the normalized far-zone cutoff coefficient of the standard kernel function in the medium and high frequency bands is about 0.2. Figure 2 In the lower frequency band, the conversion rate of the measured data is larger. Since the measured data contains more long-wave errors, the long-wave errors of the elevation anomaly calculated using the standard Stokes kernel function are larger. In the higher frequency band, the amplitude of the conversion rate of the measured data when using the standard Stokes kernel function is about 1±0.2, which deviates relatively greatly from the ideal value of 1.

[0096] The modified Stokes kernel function includes two types: low-order correction and high-low order correction. Since the difference in the spectral characteristics of the two is only in the terminal frequency band of the measurement data, and the amount of elevation anomaly information in this frequency band is very small, the spectrum analysis of the present invention is only carried out for the low-order modified Stokes kernel function. In order to illustrate the spectral characteristics of the cosine-corrected Stokes kernel function under different correction frequency band conditions, the normalized far-zone truncation coefficient and conversion rate of the cosine-corrected Stokes kernel function in the frequency bands of 100-300, 200-400, 2-720, and 2-2160 are shown in Figure 1. Figure 3 to Figure 4 .

[0097] Figure 3 It shows that under the four correction conditions, the normalized far-zone truncation coefficient of the Stokes kernel function in the 100-300 frequency band of the cosine low-order correction is the largest, followed by the coefficient in the 200-400 frequency band, and the far-zone truncation coefficient in the 2-720 and 2-2160 frequency bands is very small. Figure 4 In the low-frequency band before correction, the conversion rate of the reference gravity data is approximately 1, and the conversion rate of the measured gravity data is approximately 0, indicating that the low-frequency elevation anomaly information mainly comes from the contribution of the reference model gravity data. In the correction frequency band, with the increase of the spherical harmonic order, the conversion rate of the reference data decreases and the conversion rate of the measured data increases. At the intermediate order of correction, the conversion rate of the reference data and the measured data is close to 0.5. At the cutoff frequency of the correction, the conversion rate of the reference data is approximately 1, and the conversion rate of the measured data is approximately 0.

[0098] In order to improve the accuracy of the elevation anomaly or geoid height calculated by Stokes integral, the conversion rate of the reference model gravity data in the ultra-low frequency band with very high accuracy should be close to 1, and the conversion rate in the high frequency band with significantly lower accuracy than the measured gravity data should be close to 0. The present invention sets the thresholds of the corresponding reference gravity data errors to 0.01mGal and 1mGal, respectively, that is, when the reference model gravity data error is less than 0.01mGal, the conversion rate of the reference data should be 1, and when the reference model gravity data error is greater than 1mGal, the conversion rate of the reference data should be 0. To this end, the lowest order corresponding to the cumulative error of the model gravity anomaly exceeding 0.01mGal can be set as the initial correction order L 0 The minimum order corresponding to the cumulative error of the model gravity anomaly exceeding 1mGal is set as the cutoff correction order L 1 After calculation, when using the EIGEN-6C4 reference model, L 0 The experience value is 104, L 1 The experience value is 305.

[0099] In order to test the practical application effect of different Stokes kernel functions, a quasi-geoid construction (i.e., elevation anomaly grid data calculation) experiment was carried out in the test area in central China. The test area is a mountainous area with a latitude and longitude range of 106°~116°E, 26°~34°N. The terrain data used in this experiment is SRTM data. The terrain of the test area with a 2′ resolution is 3423.56m at the highest, 11.08m at the lowest, 537.40m at the average, and 714.51m at the root mean square. The regional quasi-geoid is constructed according to the Stokes-Helmert boundary value theory method, and the downward continuation method is used for analytical continuation. The continuation integral radius and the Stokes integral radius are both 1°. The range of the constructed regional gravity quasi-geoid is 108°~114°E, 28°~32°N, and the reference model used is the 2190-order EIGEN-6C4 gravity field model. The accuracy check of the gravity quasi-geoid is carried out by measuring GNSS level point data. The 2′ resolution gravity quasi-geoid solution results based on different Stokes kernel functions are shown in Table 1. The parameters L of the high and low order cosine corrected Stokes kernel functions are 0 Take 104, L 1 Take 305, L 2 Take 5400, L 3 Take 5600.

[0100] Table 1 Quasi-geoid accuracy (unit: m) and calculation time (unit: s) calculated by different Stokes kernel functions

[0101]

[0102] According to Table 1, among several statistics, the standard deviation index (SD) is usually used to characterize the accuracy of the quasi-geoid that does not contain systematic deviations. In addition to the kernel function correction method and the model gravity field error, the accuracy of the quasi-geoid construction is also affected by factors such as the local gravity measurement data error and the GNSS leveling point error. As can be seen from Table 1, the quasi-geoid calculated by the standard Stokes kernel function has the largest error, so it is not suitable for solving the (quasi) geoid. The accuracy of the quasi-geoid calculated using the low-order cosine correction kernel function is improved by 47.61% compared with the traditional standard kernel function. Compared with the corresponding kernel function with low-order correction, the accuracy of the quasi-geoid in the test area solved by the Stokes kernel function with high and low order cosine correction changes very little (0.01cm), but the calculation time is greatly increased. This experiment shows that it is unnecessary to make high-order corrections to the Stokes kernel function.

[0103] Finally, in order to illustrate the influence of the correction frequency band on the accuracy of the quasi-geoid construction, the gravity quasi-geoid is calculated based on multiple sets of low-order cosine-corrected Stokes kernel function correction order test values. Table 2 shows the accuracy statistics of the gravity quasi-geoid.

[0104] Table 2 Quasi-geoid accuracy based on test values ​​of different correction orders (unit: m)

[0105] Correction frequency band Maximum Minimum average value RMS Standard Deviation 2～300 0.1844 -0.1339 0.0878 0.1055 0.0590 2～400 0.1907 -0.0911 0.0881 0.1042 0.0561 100～300 0.1797 -0.1009 0.0872 0.1034 0.0561 200～400 0.2142 -0.0603 0.0890 0.1091 0.0636

[0106] According to Table 2, the correction frequency band has a great influence on the accuracy of the quasi-geoid. The conversion rates of measured and reference gravity data are different under different correction frequency bands. Under the correction frequency bands of 2-300, 2-400, 100-300, and 200-400, the conversion rates of reference and measured gravity data are all 0.5, and the corresponding orders are about 151, 200, 200, and 300, respectively. As shown in Table 3, the solution accuracy of the correction frequency bands of 2-400 and 100-300 is consistent (5.61cm), which is better than the solution accuracy of the correction frequency bands of 2-300 and 200-400 (5.90cm and 6.36cm). This shows that when the spherical harmonic order increases to around 200, it is more reasonable to reduce the conversion rate of reference gravity data to about 0.5 and increase the conversion rate of measured gravity data to about 0.5. The intermediate order of the correction frequency band (i.e., the average of the initial and cutoff correction orders) should theoretically be close to the order corresponding to the equal order error of the measured gravity anomaly and the model gravity anomaly. When the correction kernel function is actually applied, multiple groups of correction order test values ​​can be selected based on the empirical values ​​of the correction parameters, combined with the order of the reference gravity field model used and the accuracy of the test area measurement and reference gravity data, to conduct multiple tests to determine the Stokes kernel function correction frequency band that is more suitable for the local (quasi) geoid.

[0107] The spectrum analysis and boundary value solution test of this embodiment are both performed on the extended Stokes kernel function. The non-extended Stokes kernel function is a special case of the extended Stokes kernel function. Therefore, the analysis and conclusions of the present invention are also applicable to the non-extended Stokes kernel function for calculating the geoid fluctuation.

[0108] Example 2

[0109] A computer system includes: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the method for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function in embodiment 1.

[0110] Example 3

[0111] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of the method for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function in embodiment 1.

[0112] Example 4

[0113] A computer program product includes a computer program, which, when executed by a processor, implements the steps of the method for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function in embodiment 1.

[0114] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0115] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A method for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function, characterized in that: The method comprises: Obtaining the earth's gravity field model, ground gravity anomaly data, terrain data, GNSS level point data, the type of elevation reference surface and the test frequency bandwidth; the types of elevation reference surfaces include geoid and quasi-geoid; the geoid is the reference surface of the orthometric elevation system, and the quasi-geoid is the reference surface of the normal height elevation system; The model data, the direct effect of terrain compression on ground gravity anomaly and the indirect effect of terrain compression on elevation datum are calculated by using the earth gravity field model, terrain data and high-frequency terrain influence combined algorithm; the model data includes one of the model geoid height and the model elevation anomaly and the model gravity anomaly; the indirect effect of terrain compression on elevation datum is the indirect effect of terrain compression on geoid or the indirect effect of terrain compression on quasi-geoid; Determine the empirical value and test value of the correction order of the low-order cosine-corrected Stokes kernel function based on the accumulated error of the model gravity anomaly and the test frequency bandwidth; the correction order includes an initial correction order and a cutoff correction order; The empirical elevation datum and the test elevation datum are calculated by using a low-order cosine-corrected Stokes kernel function based on ground gravity anomaly data, terrain data, model data, the direct effect of terrain compression on ground gravity anomaly, the indirect effect of terrain compression on elevation datum, the empirical value of the correction order and the test value of the correction order; the empirical elevation datum is an empirical geoid or an empirical quasi-geoid; the test elevation datum is a test geoid or a test quasi-geoid; The final elevation datum is determined based on GNSS level point data, the empirical elevation datum and the test elevation datum; the final elevation datum is the final geoid or the final quasi-geoid.

2. The method for constructing an elevation reference surface based on a low-order cosine-corrected Stokes kernel function according to claim 1, characterized in that: The earth gravity field model, terrain data and high-frequency terrain influence combined algorithm are used to calculate the model data, the direct impact of terrain compression on ground gravity anomaly and the indirect impact of terrain compression on elevation datum, including: If the type of the elevation datum is the geoid, the model geoid height and the ground model gravity anomaly are calculated using the earth gravity field model and terrain data; A high-frequency terrain impact combined algorithm is used to calculate the direct impact of terrain compression on ground gravity anomaly and the indirect impact of terrain compression on the geoid based on terrain data; If the type of the elevation datum is a quasi-geoid, the model elevation anomaly and the ground model gravity anomaly are calculated using the earth gravity field model and terrain data; A high-frequency terrain influence combined algorithm is used to calculate the direct effect of terrain compression on ground gravity anomaly and the indirect effect of terrain compression on the quasi-geoid based on terrain data.

3. The method for constructing an elevation reference surface based on a low-order cosine-corrected Stokes kernel function according to claim 1, characterized in that: The low-order cosine modified Stokes kernel function includes two forms: non-extended and extended, specifically: Where S(ψ) represents the non-extended low-order cosine modified Stokes kernel function; S(r P ,ψ) represents the extended low-order cosine modified Stokes kernel function; rp represents the geocentric radius of the calculation point; ψ represents the angular distance between the calculation point and the flow integration point; represents the n-order spectral weight of the low-order cosine modified Stokes kernel function; n represents the order; R represents the radius of the earth; P n (cosψ) represents the nth-order Legendre polynomial; L0 is the initial correction order, and L1 is the cutoff correction order.

4. The method for constructing an elevation reference surface based on a low-order cosine-corrected Stokes kernel function according to claim 1, characterized in that: Based on the cumulative error of the model gravity anomaly and the test frequency bandwidth, the empirical value and test value of the correction order of the low-order cosine-corrected Stokes kernel function are determined, including: Calculate the cumulative errors of each order of the model gravity anomaly based on the calculation formula of the cumulative error of the model gravity anomaly; The lowest order corresponding to the cumulative error of the model gravity anomaly exceeding 0.01mGal is set as the empirical value of the initial correction order; The lowest order corresponding to the cumulative error of the model gravity anomaly exceeding 1mGal is set as the empirical value of the cutoff correction order; The empirical value of the initial correction order plus or minus a multiple of 10 or 100 within the frequency band of the test frequency bandwidth is set as the test value of the initial correction order; The empirical value of the cutoff correction order is added to or subtracted from a multiple of 10 or 100 within the frequency band of the test frequency bandwidth and set as the test value of the cutoff correction order.

5. The method for constructing an elevation reference surface based on a low-order cosine-corrected Stokes kernel function according to claim 4, characterized in that: The calculation formula of the cumulative error of the model gravity anomaly is: Where L is the cutoff order of the reference gravity field model, δΔg G (L) represents the cumulative error of the model gravity anomaly truncated to the Lth order, represents the gravity anomaly error of the n-order model, γ represents normal gravity, δC nm Represents the spherical harmonic coefficient term C of the gravity field model nm The error, δS nm Represents the spherical harmonic coefficient term S of the gravity field model nm The error, σ represents the standard deviation, n represents the order, and m represents the degree.

6. The method for constructing an elevation reference surface based on a low-order cosine-corrected Stokes kernel function according to claim 1, characterized in that: The low-order cosine-corrected Stokes kernel function is used to calculate the empirical elevation datum and the test elevation datum based on ground gravity anomaly data, terrain data, model data, the direct impact of terrain compression on ground gravity anomaly, the indirect impact of terrain compression on elevation datum, the empirical value of the correction order and the test value of the correction order, including: Subtracting the ground gravity anomaly from the ground model gravity anomaly to obtain a gravity anomaly difference; Adding the gravity anomaly difference and the direct effect of the terrain compression on the ground gravity anomaly to obtain the residual ground Helmert gravity anomaly; Performing continuation processing on the remaining ground Helmert gravity anomaly to obtain a processed remaining boundary surface Helmert gravity anomaly; If the type of the elevation datum is the geoid, a non-extended low-order cosine-corrected Stokes kernel function is used to perform Stokes integral processing on the processed residual boundary surface Helmert gravity anomaly based on the empirical value of the correction order and the test value of the correction order to obtain the empirical residual geoid height and the test residual geoid height; Adding the empirical residual geoid height, the indirect effect of terrain compression on the geoid and the model geoid height to obtain an empirical geoid; Adding the test residual geoid height, the indirect effect of the terrain compression on the geoid and the model geoid height to obtain a test geoid; If the type of the elevation reference surface is a quasi-geoid, an extended low-order cosine-corrected Stokes kernel function is used to perform Stokes integral processing on the processed residual boundary surface Helmert gravity anomaly based on the empirical value of the correction order and the test value of the correction order to obtain the empirical residual elevation anomaly and the test residual elevation anomaly; Adding the empirical residual elevation anomaly, the indirect effect of the terrain compression on the quasi-geoid and the model elevation anomaly to obtain an empirical quasi-geoid; The test residual elevation anomaly, the indirect effect of the terrain compression on the quasi-geoid and the model elevation anomaly are added together to obtain the test quasi-geoid.

7. The method for constructing an elevation reference surface based on a low-order cosine-corrected Stokes kernel function according to claim 1, characterized in that: The final elevation datum is determined based on the empirical elevation datum and the test elevation datum, including: If the type of the elevation datum is a geoid, the final elevation datum is determined based on the empirical geoid and the test geoid using GNSS level point data; If the type of the elevation datum is a quasi-geoid, the final elevation datum is determined based on the empirical quasi-geoid and the test quasi-geoid using GNSS level point data.

8. A computer system comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that the processor executes the computer program to implement the steps of the elevation datum construction method based on the low-order cosine corrected Stokes kernel function as described in any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the method for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function as described in any one of claims 1 to 7 are implemented.

10. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the steps of the method for constructing an elevation datum based on a low-order cosine-corrected Stokes kernel function as described in any one of claims 1 to 7 are implemented.

Images