Geoid construction method, system, medium and product based on improved Molodensky modified Stokes kernel function

Abstract

The invention discloses a geoid construction method, system, medium and product based on an improved Molodensky modified Stokes kernel function, and relates to the field of physical geodetic measurement. The method comprises obtaining ground grid gravity anomaly based on ground discrete gravity point data; calculating model gravity anomaly, model geoid height, direct influence and indirect influence; determining initial and cutoff correction orders according to the accumulated error of model gravity anomaly; calculating an improved Molodensky modified Stokes kernel function according to the initial and cutoff correction orders; using the improved Molodensky modified Stokes kernel function to determine the regional geoid based on ground grid gravity anomaly, direct influence, indirect influence, model gravity anomaly, model geoid height, initial correction order, cutoff correction order and GNSS level point data. The invention can improve the construction accuracy of the regional geoid.

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 a geoid based on an improved Molodensky modified Stokes kernel function. Background Art

[0002] Gravity measurement information has important application value in many disciplines such as geodesy, space science, oceanography, geophysics, geodynamics, etc. In the field of physical geodesy research, geoid refinement and vertical deviation modeling are two important aspects of gravity measurement data application.

[0003] As the starting surface of the orthometric elevation system, the geoid plays an irreplaceable role in the construction of the modern elevation datum system. The precise geoid model will fundamentally change the operation mode of elevation measurement. In addition to the direct measurement method based on geometric relationships, the refined geoid can also be calculated by indirect calculation based on physical relationships. The indirect method refers to establishing a geodetic boundary value problem about the gravity disturbance potential based on one or more gravity data as boundary conditions, determining the disturbance potential function by solving the boundary value problem, and then converting it into the geoid height by the Bruns formula.

[0004] The research on the solution theory and method of geodetic boundary value problem has gone through three major development stages in the world, and three famous geodetic boundary value theories have been born, namely Stokes geoid boundary value theory, Molodensky earth natural surface boundary value theory and Bjerhammar virtual sphere boundary value theory. In order to ensure that the disturbance is outside the geoid and meets the harmonic hypothesis, the Stokes boundary value theory requires that there is no terrain mass outside the geoid, so the terrain mass outside must be removed. Because no matter what method is used to remove the terrain mass, the geoid will change and produce indirect effects, so the geoid calculated by Stokes integral is the adjusted geoid. At present, the terrain mass adjustment scheme that is widely used in practice is the Helmert second compression method, and its corresponding boundary value problem solution is called the Stokes-Helmert method.

[0005] In theory, Stokes integral is a global integral, but in practice, the removal and recovery technique that only requires near-area integration is usually used. The removal and recovery technique can effectively reduce the far-area truncation error of the near-area Stokes integral, thereby improving the accuracy of (quasi) geoid construction. When applying the removal and recovery technique in Stokes integral, the residual high-frequency gravity anomaly is first subtracted from the measured gravity anomaly data to obtain the residual high-frequency gravity anomaly, and then the Stokes integral is performed on it. Finally, the residual high-frequency geoid height obtained by the Stokes integral is summed with the model geoid height to obtain the geoid. The kernel function of the Stokes integral has an important influence on the accuracy of geoid construction. In order to improve the accuracy of geoid construction, some geodetic scholars have carried out research on the optimization method of the Stokes kernel function. The Molodensky correction method is an important kernel function correction method, but the traditional Molodensky correction method is a correction method for non-removal and recovery mode and is not applicable to the removal and recovery mode. Therefore, it is necessary to study the Molodensky correction method suitable for the removal and recovery mode to improve the accuracy of geoid construction. Summary of the invention

[0006] The purpose of the present invention is to provide a geoid construction method, system, medium and product based on an improved Molodensky modified Stokes kernel function, which can improve the construction accuracy of the geoid.

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

[0008] A method for constructing a geoid based on an improved Molodensky modified Stokes kernel function includes:

[0009] Obtain regional ground discrete gravity point data, terrain data, different earth gravity field reference models and GNSS level point data;

[0010] Grid processing is performed on the discrete ground gravity point data to obtain the ground grid gravity anomaly;

[0011] The model gravity anomaly and model geoid height are calculated using different reference models of the Earth's gravity field;

[0012] Calculate direct and indirect impacts using terrain data;

[0013] Determine the initial correction order and cut-off correction order of different earth gravity field reference models according to the accumulated error of model gravity anomaly;

[0014] The improved Molodensky modified Stokes kernel function is calculated based on the initial correction order and the cutoff correction order;

[0015] The regional geoid is determined by using an improved Molodensky modified Stokes kernel function based on ground grid gravity anomaly, direct influence, indirect influence, model gravity anomaly, model geoid height, initial correction order, cutoff correction order and GNSS level point data.

[0016] 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 a geoid based on an improved Molodensky modified Stokes kernel function.

[0017] 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 a geoid based on an improved Molodensky modified Stokes kernel function.

[0018] 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 a geoid based on an improved Molodensky modified Stokes kernel function.

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

[0020] The invention discloses a geoid construction method, system, medium and product based on an improved Molodensky modified Stokes kernel function, and relates to the technical field of physical geodetic measurement. The method comprises the following steps: obtaining regional ground discrete gravity point data, terrain data, different earth gravity field reference models and GNSS level point data; performing gridding processing on the ground discrete gravity point data to obtain ground grid gravity anomaly; calculating model gravity anomaly and model geoid height by using different earth gravity field reference models; calculating direct influence and indirect influence by using terrain data; determining initial correction order and cutoff correction order of different earth gravity field reference models according to the accumulated error of the model gravity anomaly; calculating an improved Molodensky modified Stokes kernel function according to the initial correction order and the cutoff correction order; and determining the regional geoid by using the improved Molodensky modified Stokes kernel function based on ground grid gravity anomaly, direct influence, indirect influence, model gravity anomaly, model geoid height, initial correction order, cutoff correction order and GNSS level point data. The present invention can improve the construction accuracy of the regional geoid. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] 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.

[0022] Figure 1 Schematic diagram of the spectral weights of the Molodensky modified Stokes kernel function;

[0023] Figure 2 Schematic diagram of the normalized far-zone cutoff coefficient of the Molodensky modified Stokes kernel function;

[0024] Figure 3 This is a schematic diagram of the conversion rate of gravity data when using the Molodensky modified Stokes kernel function;

[0025] Figure 4 This is a schematic diagram of the cumulative error of the model gravity anomaly;

[0026] Figure 5 This is the topographic map of the test area;

[0027] Figure 6 Measure gravity anomaly distribution map for the test area;

[0028] Figure 7 The quasi-geoid and GNSS level point distribution map of the test area;

[0029] Figure 8 Schematic diagram of the process of constructing the geoid based on the improved Molodensky modified Stokes kernel function. DETAILED DESCRIPTION

[0030] 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.

[0031] The purpose of the present invention is to provide a geoid construction method, system, medium and product based on an improved Molodensky modified Stokes kernel function, aiming to improve the accuracy of geoid construction.

[0032] 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.

[0033] Example 1

[0034] like Figure 8 As shown, a method for constructing a geoid based on an improved Molodensky modified Stokes kernel function in this embodiment includes:

[0035] Step 101: Acquire regional ground discrete gravity point data, terrain data, different earth gravity field reference models and GNSS level point data.

[0036] Step 102: Gridding the ground discrete gravity point data to obtain ground grid gravity anomalies.

[0037] The observed gravity values ​​of the ground discrete gravity points are corrected spatially, interlayer-wise and locally, to obtain corrected ground discrete gravity points.

[0038] The corrected ground discrete gravity point is subtracted from the preset gravity value to obtain the complete Bouguer gravity anomaly of the ground discrete gravity point; the preset gravity value is the gravity of a normal ellipsoid that is close to the earth and has a regular mass distribution, which describes the regular part of the earth's gravity field, and the value is calculated using the normal gravity formula.

[0039] The complete Bouguer gravity anomaly of the discrete gravity points on the ground is interpolated to obtain the complete Bouguer gravity anomaly of the grid points.

[0040] The complete Bouguer gravity anomaly of the grid points is restored with spatial correction, inter-layer correction and local terrain correction to obtain the ground grid gravity anomaly.

[0041] Step 103: Calculate the model gravity anomaly and the model geoid height using different earth gravity field reference models.

[0042] The disturbance potential model is obtained by subtracting the 2nd, 4th, 6th, 8th, and 10th order band harmonic coefficients in the spherical harmonic coefficients of different earth gravity field models from the 2nd, 4th, 6th, 8th, and 10th order band harmonic coefficients in the preset gravity potential. The preset gravity potential is the gravity potential generated by the normal ellipsoid.

[0043] Calculate the model gravity anomaly and model geoid height based on the disturbance potential model:

[0044]

[0045]

[0046] Where, Δg M is the model gravity anomaly, n is the order, m is the number, N M is the model geoid height, r P is the radius from the center of the earth at the calculation point, r Eis the geocentric radius of the projection point of the calculation point on the reference ellipsoid, GM is the gravitational constant, a represents the major semi-axis of the reference ellipsoid, L is the truncation order of the gravity field model, is the geocentric latitude of the calculation point, λ is the geocentric longitude of the calculation point, is a fully normalized adjoint Legendre function of order n and degree m, The disturbance potential model The spherical harmonic coefficients of the term, The disturbance potential model N0 is the spherical harmonic coefficient of the geoid height, C is the correction of elevation anomaly-geoid height, and C is calculated by the spherical harmonic coefficient of elevation anomaly-geoid height correction truncated to the Lth order.

[0047] Step 104: Calculate direct impact and indirect impact using terrain data.

[0048] Calculates the second, third, and fourth power spherical harmonic coefficients of elevation based on terrain data.

[0049] Based on the quadratic spherical harmonic coefficient of elevation, the cubic spherical harmonic coefficient of elevation, the quartic spherical harmonic coefficient of elevation and terrain data, the direct impact and the indirect impact are calculated according to the high-frequency terrain impact combination algorithm; the direct impact is the direct impact of terrain compression on ground gravity anomaly, and the indirect impact is the indirect impact of terrain compression on the geoid.

[0050] Step 105: Determine the initial correction order and the cutoff correction order of different earth gravity field reference models according to the accumulated error of the model gravity anomaly.

[0051] Calculate the cumulative error of the model gravity anomaly:

[0052]

[0053] Among them, L is the cutoff order of the gravity field model, σ is the standard deviation, n is the order, m is the number, δΔg G (L) represents the cumulative error of the model gravity anomaly truncated to the Lth order, γ 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 of error.

[0054] When the cumulative error of the model gravity anomaly exceeds 0.01 mGal, the corresponding lowest order is determined as the initial correction order.

[0055] When the cumulative error of the model gravity anomaly exceeds 1 mGal, the corresponding lowest order is determined as the cutoff correction order.

[0056] Step 106: Calculate an improved Molodensky modified Stokes kernel function according to the initial correction order and the cutoff correction order.

[0057] The improved Molodensky modified Stokes kernel function is:

[0058]

[0059] Among them, S MOL (ψ) represents the improved Molodensky modified Stokes kernel function, t n represents the n-order Molodensky correction coefficient, P n (cosψ) represents the nth-order Legendre polynomial, ψ represents the angular distance between the flow point and the calculation point, S BL (ψ) represents the low-order truncated Stokes kernel function, L0 represents the initial correction order, L1 represents the cutoff correction order, and ψ1 is the near-zone integration radius.

[0060] S BL The calculation formula of (ψ) is

[0061]

[0062] Here, s represents sin(ψ / 2).

[0063] The calculation formula of Molodensky correction coefficient is:

[0064]

[0065] in, is the m-th order far-field truncation coefficient of the non-expanded low-order truncated Stokes kernel function, e mn (ψ1) is the coefficient matrix element. mn (ψ1), The calculation formulas are

[0066]

[0067] Step 107: Determine the regional geoid using the improved Molodensky modified Stokes kernel function based on ground grid gravity anomaly, direct impact, indirect impact, model gravity anomaly, model geoid height, initial correction order, cutoff correction order and GNSS level point data.

[0068] The residual ground gravity anomaly is obtained by subtracting the ground grid gravity anomaly from the model gravity anomaly.

[0069] The residual ground gravity anomaly is added to the direct impact to obtain the residual ground Helmert gravity anomaly.

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

[0071] The improved Molodensky modified Stokes kernel function is used to perform Stokes integral processing on the processed residual boundary surface Helmert gravity anomaly to obtain the residual geoid height.

[0072] The residual geoid height, the indirect impact and the model geoid height are added together to obtain a geoid based on different earth gravity field reference models; the earth gravity field model used for the model geoid height is consistent with the earth gravity field model used for the model gravity anomaly.

[0073] The GNSS level point data is used to verify the accuracy of the geoid based on different earth gravity field reference models and determine the regional geoid.

[0074] Standard Stokes kernel function

[0075] Stokes integral includes two forms: non-extended and extended. The non-extended Stokes integral is used to calculate the geoid height, while the extended Stokes integral is used to calculate the elevation anomaly. The standard Stokes kernel function in the extended form is:

[0076]

[0077] Where S ST (r P ,ψ) is the extended form of the 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:

[0078] The standard Stokes kernel function in non-expanded form is:

[0079]

[0080] Where S ST (ψ) is the non-extended standard kernel function, and its symbol s is expressed as sin(ψ / 2). ST (ψ) is S ST (r P ,ψ) in the parameter r PThis is a special case when R is taken, so the present invention takes the extended kernel function as an example to study the Stokes kernel function.

[0081] Truncated Stokes kernel function

[0082] At present, the truncated Stokes kernel function is widely used in constructing (quasi) geoid based on the removal and restoration technology. The truncated kernel function includes two forms: high-order and low-order truncation and low-order truncation only. The geoid accuracy constructed by the Stokes kernel function of the two truncated kernel functions is basically the same, but the construction efficiency of the low-order truncated Stokes kernel function is much higher than that of the high-order and low-order truncated Stokes kernel function. Therefore, the present invention only gives the expression of the extended form of the low-order truncated Stokes kernel function, that is,

[0083]

[0084] Where S BL (r P ,ψ) is the extended low-order truncated Stokes kernel function, and L1 is its low-order truncation order. P Taking R as R, we can get the following non-expanded low-order truncated Stokes kernel function:

[0085]

[0086] Where S BL (ψ) is the non-extended low-order truncated Stokes kernel function, and s represents sin(ψ / 2).

[0087] According to the relationship between the low-order truncation order and the truncation order (L) of the gravity field model, the truncated Stokes kernel function can be divided into two types: the truncated kernel function with L1=L and the truncated kernel function with L1<L.

[0088] Molodensky modified Stokes kernel function

[0089] The expression of the traditional extended form of the Molodensky modified Stokes kernel function is:

[0090]

[0091] In the formula is the traditional Molodensky modified Stokes kernel function, L1 is the highest Molodensky modified order, t n is the n-order Molodensky correction coefficient, which is determined according to the principle of minimizing the far-zone truncation error. Its calculation formula is:

[0092]

[0093] Where ψ1 is the near zone integral radius, For S BL (r P ,ψ) The m-order far-zone cutoff coefficient of the kernel function, e mn (ψ1) is the coefficient matrix element calculated by the m-order and n-order Legendre polynomials. mn (ψ1), The calculation formulas are:

[0094]

[0095] Based on the traditional Molodensky correction method, the present invention proposes the following extended form of improved Molodensky modified Stokes kernel function:

[0096]

[0097] In the formula It is an extended form of the improved Molodensky modified Stokes kernel function, with an initial correction order of L0, a cutoff correction order of L1, and a correction coefficient t n Calculated by the following system of equations:

[0098]

[0099] The parameter r in the extended Molodensky modified Stokes kernel function P Taking R as R, we can get the following non-expanded Molodensky modified Stokes kernel function:

[0100]

[0101] Where s represents sin(ψ / 2), t n is the n-order Molodensky correction coefficient, and its calculation formula is:

[0102]

[0103] in, is the m-order far-zone truncation coefficient of the non-extended low-order truncated Stokes kernel function, and its calculation formula is:

[0104]

[0105] Comparing the traditional and improved Molodensky modified Stokes kernel functions, it can be seen that the traditional Molodensky modified Stokes kernel function is a special case of the improved Molodensky modified Stokes kernel function when the initial correction order is 2.

[0106] The Stokes kernel function can be expressed in the following general spectral form:

[0107]

[0108] Where S(r P ,ψ) represents the Stokes kernel function, ω n Represents the n-order spectral weight. The difference between different types of Stokes kernel functions lies in the different values ​​of the spectral weight. The spectral weights of each order of the standard Stokes kernel function are 1, while the spectral weight of the removed order in the truncated Stokes kernel function is 0 and the spectral weights of the remaining orders are 1. The spectral weight of the traditional Molodensky modified Stokes kernel function is:

[0109]

[0110] The spectral weight of the improved Molodensky modified Stokes kernel function is:

[0111]

[0112] in, represents the n-order spectral weight of the improved Molodensky modified Stokes kernel function, where n represents the order and r p represents the diameter of the calculation point from the center of the earth, R represents the radius of the earth, ψ represents the angular distance between the calculation point and the flow integration point, L0 represents the low-order initial correction order, and L1 represents the low-order cutoff correction order.

[0113] The near-zone Stokes integral can be expressed in the following general spectral form:

[0114]

[0115] Where ζ is the elevation anomaly, C1 represents the near area with the calculation point as the center and ψ1 as the radius, γ is the normal gravity of the ground calculation point, and Δg n is the nth-order spherical gravity anomaly, Δg is the spherical gravity anomaly, is the normalized n-order far-zone truncation error coefficient, which is calculated as follows:

[0116]

[0117] It can be further deduced that the general spectral decomposition formula of Stokes integral under the removal recovery mode is:

[0118]

[0119] 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. 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. Substituting the spectral weight factors of different kernel functions into the formula can obtain the specific conversion rate of measured and reference gravity data when using different kernel functions. P By replacing R, γ with the normal gravity γ0 of the boundary surface, and the elevation anomaly ζ with the geoid height N, we can obtain the corresponding spectral decomposition formula for the Stokes integral for calculating the geoid height in the removal recovery mode.

[0120] In order to compare and illustrate the spectral characteristics of the traditional and improved Molodensky modified Stokes kernel functions, the spectral weights, normalized far-zone cutoff coefficients and conversion rates of the Molodensky modified Stokes kernel functions in the frequency bands of 2-300, 2-400, 2-720, 2-2160, 100-300 and 200-400 are shown in Fig. Figures 1 to 3 In the spectrum 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.

[0121] like Figure 1 As shown in the figure, the spectral weight of the Stokes kernel function of the Molodensky correction 2~2160 frequency band seriously deviates from the normal value range of the spectral weight ([0,1]). The minimum spectral weight of the Molodensky correction 2~2160 frequency band is -1.34, the conversion rate of the corresponding order of reference gravity data is 2.34, and the conversion rate of the measured gravity data is -1.34. Since the errors of the two gravity data are different, the converted elevation anomaly contains more errors. Therefore, the Molodensky correction order should not be too high. Among the several correction frequency bands, the kernel function starting from 2 is the traditional Molodensky corrected Stokes kernel function, and the kernel function of the Molodensky correction 100~300 and 200~400 frequency bands is the improved Molodensky corrected Stokes kernel function. From Figure 1It can also be clearly seen that the L0 (starting correction order) spectral weight of the kernel function of the Molodensky correction 100-300 and 200-400 frequency bands is greater than 0, and the L1 (cutoff correction order) spectral weight is less than 1. Figure 2 It can be seen that in the correction frequency band, the far-zone cutoff coefficient of the Molodensky correction kernel function is approximately 0. In the non-correction frequency band, the normalized far-zone cutoff coefficient of the traditional Molodensky correction Stokes kernel function is very small, while the normalized far-zone cutoff coefficient of the improved Molodensky correction Stokes kernel function is relatively large. Figure 3 It can be seen that in the low-frequency band, the conversion rate of reference data decreases rapidly when the traditional Molodensky correction is in the 2-300 and 2-400 frequency bands. Therefore, compared with the traditional Molodensky modified Stokes kernel function, the improved Molodensky modified Stokes kernel function can increase the contribution of reference gravity data in low-order elevation anomalies / geoid fluctuations.

[0122] In the boundary value problem based on the removal recovery technique, the gravity field model provides long-wave reference gravity field information. At present, different research institutions around the world have released several gravity field models with different accuracy. The calculation formula for the error of the model gravity anomaly calculated by the gravity field model is:

[0123]

[0124] Where L is the cutoff order of the reference gravity field model, σ(δΔg G ) represents the cumulative error of the model gravity anomaly truncated to the Lth order, represents the gravity anomaly error of the n-order model, The formula is:

[0125]

[0126] Where GM is the gravitational constant, δC nm ,δS nm They are the spherical harmonic coefficients of the gravity field model C nm , S nm of error. Figure 4 The accumulated errors of various orders of the EGM2008 model released by the National Geospatial-Intelligence Agency of the United States, the EIGEN-6C4 model released by the German Geoscience Center GFZ, and the XGM2019e model released by the Technical University of Munich are plotted.

[0127] from Figure 4It can be seen that the ultra-long wave error of the model gravity data is very small, among which the ultra-long wave accuracy of the XGM2019e and EIGEN-6C4 models is significantly better than that of EGM2008. As the order of the gravity field model increases, the error of the model gravity data tends to increase. It should be noted that the reliability of the spherical harmonic coefficients of the gravity field model usually decreases with the increase of the spherical harmonic order.

[0128] 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 exceeds 1mGal, the conversion rate of the reference data should be 0. For this reason, the present invention sets the lowest order corresponding to when the cumulative error of the reference model gravity data exceeds 0.01mGal to L0, and the lowest order corresponding to when the cumulative error of the reference model gravity data exceeds 1mGal to L1. After calculation, the correction orders of the EGM2008 model are 45 (L0) and 140 (L1), the correction orders of the EIGEN-6C4 model are 104 (L0) and 305 (L1), and the correction orders of the XGM2019e model are 121 (L0) and 438 (L1).

[0129] In order to test the practical application effects 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 test area includes 34138 measured discrete gravity points, which are converted into 2′ resolution grid gravity anomaly data through the gridding process. The maximum value of the 2′ resolution measured gravity anomaly data is 184.17mGal, the minimum value is -121.32mGal, the average value is -21.20mGal, and the root mean square value is 36.41mGal. The terrain data used in this experiment is 3″×3″ resolution SRTM data, and the converted 2′ resolution terrain data has a maximum value of 3423.56m, a minimum value of 11.08m, an average value of 537.40m, and a root mean square value of 714.51m. Figure 5 , Figure 6 The fluctuations of the topographic data of the test area and the measured gravity anomaly data were plotted respectively.

[0130] In this experiment, the regional quasi-geoid is constructed according to the Stokes-Helmert boundary value theory method using the removal recovery technique, and the radius of the near zone of the Stokes integral is 1°. Before the Stokes integral, the surface gravity anomaly needs to be converted into the boundary surface gravity anomaly data through the downward continuation process. In the experiment, the downward continuation is realized using an analytical continuation algorithm without the need for inversion or iteration. The calculation of the radial derivative of the gravity anomaly in the analytical continuation algorithm uses a 1° near zone range. The parameters of the kernel function in the experiment are the empirical values ​​given above. When using the EGM2008 reference model, L0 is 45 and L1 is 140; when using EIGEN-6C4, L0 is 104 and L1 is 305; when using XGM2019e, L0 is 121 and L1 is 438. The truncation order of the three reference gravity field models is 2190. The range of the constructed regional gravity quasi-geoid is 108°~114°E, 28°~32°N, and its accuracy is verified by 63 measured GNSS level point data. The GNSS level points contain unknown centimeter-level precision system errors. The 2′ resolution gravity quasi-geoid solution results obtained based on different Stokes kernel functions are shown in Table 1.

[0131] Table 1 Accuracy of gravity quasi-geoid calculated by different Stokes kernel functions (unit: m)

[0132]

[0133]

[0134] 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 local gravity measurement data errors and GNSS leveling point errors. As can be seen from Table 1, the quasi-geoid calculated by the standard Stokes kernel function has the lowest accuracy, followed by the truncated kernel function, and the quasi-geoid calculated by the improved Molodensky corrected Stokes kernel function has the highest accuracy. Under the condition that the EIGEN-6C4 model is used as the reference model, the calculation accuracy of the improved Molodensky corrected Stokes kernel function is improved by about 44.81% compared with the standard kernel function, and by about 15.05% compared with the truncated kernel function. Among the three reference gravity field models, the boundary value solution accuracy of the EGM2008 model as a reference model is significantly lower than that of the EIGEN-6C4 model and the XGM2019e model, which is consistent with the above. Figure 4 The phenomenon reflected in the long-wave accuracy of the EGM2008 model is relatively low. When the EGM2008 model reference model is used for boundary value solution, the long-wave accuracy of the EGM2008 model can be improved by combining the EGM2008 model and the high-precision pure satellite model. Figure 7 Specifically, it shows the distribution of GNSS level point data and the gravity quasi-geoid of the test area constructed by the improved Molodensky modified Stokes kernel function under the condition of EIGEN-6C4 as the reference model.

[0135] In order to compare and illustrate the application effects of the traditional and improved Molodensky modified Stokes kernel functions, Table 2 shows the accuracy of the gravity quasi-geoid calculated using the traditional Molodensky modified Stokes kernel function, where the reference gravity field model is the 2190-order EIGEN-6C4 model, and the parameter L1 of the traditional Molodensky modified kernel function is 104 and 305.

[0136] Table 2 Quasi-geoid accuracy calculated by traditional Molodensky modified Stokes kernel function (unit: m)

[0137] Correction frequency band Maximum Minimum average value RMS Standard Deviation 2～104 0.2430 -0.1675 0.0916 0.1201 0.0782 2～305 0.2070 -0.1442 0.0897 0.1108 0.0655

[0138] According to Table 2, the accuracy of the quasi-geoid calculated by the Stokes kernel function in the 2-104 and 2-305 frequency bands corrected by the traditional Molodensky method is 7.82cm and 6.55cm respectively. The accuracy of the quasi-geoid calculated by the Stokes kernel function in the 104-305 frequency band corrected by the improved Molodensky method in Table 1 is 5.59cm, which is 28.52% and 14.66% higher than the traditional algorithm of the traditional Molodensky correction of the 2-104 and 2-305 frequency bands respectively. This is because compared with the traditional Molodensky correction method, the improved Molodensky correction Stokes kernel function improves the conversion rate of the reference gravity data in the low-frequency band with higher accuracy. The results of Tables 1 and 2 show that it is effective to improve the traditional Molodensky correction method when constructing the (quasi) geoid using the removal recovery technique.

[0139] Although the spectrum analysis and boundary value solution test of the present invention are both conducted on the extended Stokes kernel function, the non-extended Stokes kernel function is a special case of the extended Stokes kernel function, and thus the analysis and conclusion of the present invention are applicable to the non-extended Stokes kernel function for calculating the geoid. It should be noted that the time taken to construct a quasi-geoid using the Molodensky modified kernel function is much longer than that of the standard kernel function and the truncated kernel function. This is because when calculating elevation anomalies, the correction coefficient t of the Molodensky modified kernel function is nIt is related to the elevation of the calculation point, so the Molodensky correction coefficient for each calculation point needs to be solved by matrix inversion. However, when the Molodensky correction kernel function is used for geoid calculation, the Molodensky correction coefficient is independent of the elevation of the calculation point. At this time, only one correction coefficient needs to be calculated. Therefore, the efficiency of constructing a geoid based on the Molodensky correction kernel function is significantly higher than that of constructing a quasi-geoid based on the Molodensky correction kernel function.

[0140] In the process of constructing a regional (quasi) geoid according to the geodetic boundary value theory based on the removal and recovery technique, the Stokes integral kernel function has an important influence on the accuracy of the boundary value solution. The present invention first gives the expression of the Molodensky modified Stokes kernel function on the basis of the standard and truncated Stokes kernel functions. The starting order of the correction coefficient of the traditional Molodensky correction kernel function is the lowest order of the disturbed potential field. The present invention further proposes an improved Molodensky correction method that starts the correction from the non-lowest order (L0 order). Compared with the traditional Molodensky correction method, the improved Molodensky correction method increases the contribution of the reference model gravity data in the constructed geoid in the low-frequency band.

[0141] Example 2

[0142] 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 the method for constructing a geoid based on an improved Molodensky modified Stokes kernel function in embodiment 1.

[0143] Example 3

[0144] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of the method for constructing a geoid based on an improved Molodensky modified Stokes kernel function in embodiment 1.

[0145] Example 4

[0146] A computer program product includes a computer program, which, when executed by a processor, implements the steps of the method for constructing a geoid based on an improved Molodensky modified Stokes kernel function in embodiment 1.

[0147] It should be noted that the object information and data involved in the present invention (including but not limited to data used for analysis, stored data, displayed data, etc.) are all information and data authorized by the object or fully authorized by all parties, and the collection, use and processing of relevant data must comply with relevant laws, regulations and standards of relevant countries and regions.

[0148] Those skilled in the art can understand that all or part of the processes in the above-mentioned embodiments can be completed by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to the memory, database or other medium used in the embodiments provided by the present invention can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetoresistive random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. As an illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM). The database involved in each embodiment provided by the present invention may include at least one of a relational database and a non-relational database. Non-relational databases may include distributed databases based on blockchains, etc., but are not limited thereto. The processor involved in each embodiment provided by the present invention may be a general-purpose processor, a central processing unit, a graphics processor, a digital signal processor, a programmable logic device, a data processing logic device based on quantum computing, etc., but are not limited thereto.

[0149] 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.

[0150] 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 a geoid based on an improved Molodensky modified Stokes kernel function, characterized in that: The method comprises: Obtain regional ground discrete gravity point data, terrain data, different earth gravity field reference models and GNSS level point data; Grid processing is performed on the discrete ground gravity point data to obtain the ground grid gravity anomaly; The model gravity anomaly and model geoid height are calculated using different reference models of the Earth's gravity field; Calculate direct and indirect impacts using terrain data; Determine the initial correction order and cut-off correction order of different earth gravity field reference models according to the accumulated error of model gravity anomaly; The improved Molodensky modified Stokes kernel function is calculated based on the initial correction order and the cutoff correction order; The regional geoid is determined by using an improved Molodensky modified Stokes kernel function based on ground grid gravity anomaly, direct influence, indirect influence, model gravity anomaly, model geoid height, initial correction order, cutoff correction order and GNSS level point data.

2. The method for constructing a geoid based on an improved Molodensky modified Stokes kernel function according to claim 1, characterized in that: The discrete ground gravity point data are gridded to obtain the ground grid gravity anomaly, including: Performing spatial correction, inter-layer correction and local terrain correction on the observed gravity values ​​of the ground discrete gravity points to obtain corrected ground discrete gravity points; Subtract the corrected discrete gravity point from the preset gravity value to obtain the complete Bouguer gravity anomaly of the discrete gravity point on the ground; The complete Bouguer gravity anomaly of the discrete gravity points on the ground is interpolated to obtain the complete Bouguer gravity anomaly of the grid points; The complete Bouguer gravity anomaly of the grid points is restored with spatial correction, inter-layer correction and local terrain correction to obtain the ground grid gravity anomaly.

3. The method for constructing a geoid based on an improved Molodensky modified Stokes kernel function according to claim 1, characterized in that: The model gravity anomaly and model geoid height are calculated using different reference models of the Earth's gravity field, including: Subtract the 2nd, 4th, 6th, 8th, and 10th order band harmonic coefficients in the spherical harmonic coefficients of different earth gravity field models from the 2nd, 4th, 6th, 8th, and 10th order band harmonic coefficients in the preset gravity potential to obtain the disturbance potential model; Calculate the model gravity anomaly and model geoid height based on the disturbance potential model: Where, Δg M is the model gravity anomaly, n is the order, m is the number, N M is the model geoid height, r P is the radius from the center of the earth at the calculation point, r E is the geocentric radius of the projection point of the calculation point on the reference ellipsoid, GM is the gravitational constant, a represents the major semi-axis of the reference ellipsoid, L is the truncation order of the gravity field model, is the geocentric latitude of the calculation point, λ is the geocentric longitude of the calculation point, is a fully normalized adjoint Legendre function of order n and degree m, The disturbance potential model The spherical harmonic coefficients of the term, The disturbance potential model N0 is the spherical harmonic coefficient of the geoid height, C is the correction of elevation anomaly-geoid height, and C is calculated by the spherical harmonic coefficient of elevation anomaly-geoid height correction truncated to the Lth order.

4. The method for constructing a geoid based on an improved Molodensky modified Stokes kernel function according to claim 1, characterized in that: Direct and indirect impacts are calculated using terrain data, including: Calculate the second, third and fourth power spherical harmonic coefficients of elevation based on terrain data; Based on the quadratic spherical harmonic coefficient of elevation, the cubic spherical harmonic coefficient of elevation, the quartic spherical harmonic coefficient of elevation and terrain data, the direct and indirect impacts are calculated according to the high-frequency terrain influence combination algorithm; the direct impact is the direct impact of terrain compression on ground gravity anomaly, and the indirect impact is the indirect impact of terrain compression on the geoid.

5. The method for constructing a geoid based on an improved Molodensky modified Stokes kernel function according to claim 1, characterized in that: The initial correction order and the cutoff correction order of different reference models of the Earth's gravity field are determined according to the cumulative error of the model gravity anomaly, including: Calculate the cumulative error of the model gravity anomaly: Among them, L is the cutoff order of the gravity field model, σ is the standard deviation, n is the order, m is the number, δΔg G (L) represents the cumulative error of the model gravity anomaly truncated to the Lth order, γ 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 of When the cumulative error of the model gravity anomaly exceeds 0.01 mGal, the corresponding lowest order is determined as the initial correction order; When the cumulative error of the model gravity anomaly exceeds 1 mGal, the corresponding lowest order is determined as the cutoff correction order.

6. The method for constructing a geoid based on an improved Molodensky modified Stokes kernel function according to claim 1, characterized in that: The improved Molodensky modified Stokes kernel function is calculated based on the initial correction order and the cutoff correction order, specifically including: The improved Molodensky modified Stokes kernel function is: Among them, S MOL (ψ) represents the improved Molodensky modified Stokes kernel function, t n represents the n-order Molodensky correction coefficient, P n (cosψ) represents the nth-order Legendre polynomial, ψ represents the angular distance between the flow point and the calculation point, S BL (ψ) represents the low-order truncated Stokes kernel function, L0 represents the initial correction order, L1 represents the cutoff correction order, and ψ1 is the near-zone integration radius.

7. The method for constructing a geoid based on an improved Molodensky modified Stokes kernel function according to claim 1, characterized in that: The improved Molodensky modified Stokes kernel function is used to determine the regional geoid based on ground grid gravity anomaly, direct impact, indirect impact, model gravity anomaly, model geoid height, initial correction order, cutoff correction order and GNSS level point data, including: Subtracting the ground grid gravity anomaly from the model gravity anomaly to obtain a residual ground gravity anomaly; Adding the residual ground gravity anomaly to the direct impact to obtain the residual ground Helmert gravity anomaly; The remaining ground Helmert gravity anomaly is processed downward to obtain a processed remaining boundary surface Helmert gravity anomaly; The improved Molodensky modified Stokes kernel function is used to perform Stokes integral processing on the Helmert gravity anomaly of the residual boundary surface after processing to obtain the residual geoid height. Adding the residual geoid height, the indirect impact, and the model geoid height to obtain a geoid based on different earth gravity field reference models; the earth gravity field model used in the model geoid height is consistent with the earth gravity field model used in the model gravity anomaly; The GNSS level point data is used to verify the accuracy of the geoid based on different earth gravity field reference models and determine the regional geoid.

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 method for constructing a geoid based on the improved Molodensky modified 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 a geoid based on an improved Molodensky modified 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 a geoid based on an improved Molodensky modified Stokes kernel function as described in any one of claims 1 to 7 are implemented.

Images