A gravity field interpolation method and device, and storage medium

Abstract

The present invention discloses a gravity field interpolation method and device, and a storage medium, including: step S1, a server generates sampling points in a grid corresponding to a satellite orbit; step S2, the server obtains the coefficients of a polynomial according to the measured values ​​of the sampling points; step S3, a user terminal converts the ECI coordinates at a reference time into the RTN coordinates corresponding to the grid, and normalizes the RTN coordinates; the normalized RTN coordinates are brought into the polynomial, the gravitational acceleration before the spherical harmonic gravity field term is calculated using the spherical harmonic method, and the gravitational acceleration after the spherical harmonic gravity field term is calculated using the polynomial interpolation method. The technical solution of the present invention can improve the computational efficiency of the real-time orbit determination algorithm and reduce the computational complexity of on-orbit data processing.

Description

Technical Field

[0001] The present invention belongs to the field of satellite dynamics technology, and in particular relates to a gravity field interpolation method and device, and a storage medium suitable for autonomous orbit determination of a satellite-borne GNSS. Background Art

[0002] With the development of aerospace science and technology, aerospace missions such as earth observation and real-time geocoding by high-resolution remote sensing satellites and formation flying of scientific exploration satellites have put forward higher requirements on the accuracy, real-time and autonomy of satellite orbit determination.

[0003] In the real-time orbit determination algorithm of satellite-borne GPS, satellite-borne GPS pseudo-range observations are usually combined with a simplified dynamic model, and the extended Kalman filter is used for on-orbit data processing to obtain high-precision satellite orbit parameters. In the extended filter processing, the dynamic orbit integral prediction generally uses the single-step integration method of the Runge-Kutta series, and each integration requires multiple calculations of the earth's gravitational acceleration. In order to improve the real-time orbit determination accuracy of low-orbit satellites, the earth's gravity field model needs to be accurate to 50×50 or even 70×70 orders. The higher the model order, the greater the amount of calculation required to calculate the earth's gravitational acceleration using the traditional spherical harmonic function recursive algorithm.

[0004] At present, the computing power of the internal processor of the satellite GPS receiver and the satellite processor is much lower than that of the ground PC processor. In order to make the high-precision satellite GPS real-time orbit determination algorithm enter the engineering application, it is urgent to improve the computational efficiency of the real-time orbit determination algorithm, especially to optimize the dynamic orbit integration algorithm and reduce the computational complexity of on-orbit data processing. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a gravity field interpolation method and device, and a storage medium suitable for satellite-borne GNSS autonomous orbit determination, so as to improve the calculation efficiency of the real-time orbit determination algorithm and reduce the calculation amount of on-orbit data processing.

[0006] To achieve the above object, the present invention adopts the following technical solution:

[0007] A gravity field interpolation method, comprising:

[0008] Step S1, the server generates sampling points in the grid corresponding to the satellite orbit;

[0009] Step S2: The server obtains the coefficients of the polynomial according to the measured values ​​of the sampling points;

[0010] Step S3: The user terminal converts the ECI coordinates at the reference time into the RTN coordinates corresponding to the grid, and normalizes the RTN coordinates; the normalized RTN coordinates are introduced into the polynomial, and the spherical harmonic gravity field is calculated using the spherical harmonic method. The gravitational acceleration before the term is used to calculate the spherical harmonic gravity field using the polynomial interpolation method The gravitational acceleration after the term.

[0011] Preferably, step S1 comprises:

[0012] Step S11, determining the grid corresponding to the satellite orbit;

[0013] Step S12: Design sampling points in the grid according to Chebyshev node distribution.

[0014] Preferably, step S2 comprises:

[0015] Step S21, generating candidate polynomials;

[0016] Step S22, obtaining a selected polynomial according to the candidate polynomials;

[0017] Step S23, obtaining coefficients of the polynomial according to the measured values ​​corresponding to the sampling points and the selected polynomial;

[0018] Step S24, verify whether the internal accuracy of the sampling points meets the standard; if the internal accuracy meets the standard, store the coefficients of the polynomial; if the required accuracy is not achieved, increase the polynomial order and repeat the above steps S21-S23.

[0019] The present invention also provides a gravity field interpolation device, comprising:

[0020] The first processing module is used for the server to generate sampling points in the grid corresponding to the satellite orbit;

[0021] The second processing module is used for the server to obtain the coefficients of the polynomial according to the measurement values ​​of the sampling points;

[0022] The third processing module is used for the user end to convert the ECI coordinates at the reference time into the RTN coordinates corresponding to the grid, and normalize the RTN coordinates; the normalized RTN coordinates are brought into the polynomial, and the spherical harmonic gravity field is calculated using the spherical harmonic method. The gravitational acceleration before the term is used to calculate the spherical harmonic gravity field using the polynomial interpolation method The gravitational acceleration after the term.

[0023] Preferably, the first processing module comprises:

[0024] A determination unit, used to determine the grid corresponding to the satellite orbit;

[0025] The design unit is used to design the sampling points in the grid according to the Chebyshev node distribution.

[0026] Preferably, the second processing module comprises:

[0027] A generating unit, used for generating a candidate polynomial;

[0028] A selection unit, used for obtaining a selected polynomial according to the candidate polynomials;

[0029] A calculation unit, used for obtaining coefficients of the polynomial according to the measured values ​​corresponding to the sampling points and the selected polynomial;

[0030] The inspection unit is used to check whether the internal accuracy of the sampling points meets the standard; if the internal accuracy meets the standard, the coefficients of the polynomial are stored.

[0031] The present invention also provides a storage medium, on which a computer program is stored, and the computer program executes the gravity field interpolation method when running.

[0032] The present invention introduces a gravity field interpolation algorithm and applies it to the gravity field calculation of satellite-borne GPS autonomous orbit determination. In order to save the time of spherical harmonic function recursive calculation, it is proposed to divide the space into multiple regions. For each small "cage", a polynomial is selected in advance to fit the high-order terms of gravity acceleration in the region, and the coefficients of the polynomial are stored in advance. The user end only needs to convert the coordinates and call the polynomial to calculate the satellite gravity acceleration here. When calculating the local gravity field, the present invention exchanges storage space for calculation time, improves the efficiency of gravity field calculation, and the adaptive polynomial algorithm of the present invention also ensures the calculation accuracy of the gravity field. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] 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 or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.

[0034] Figure 1 This is a flow chart of a gravity field interpolation method according to an embodiment of the present invention;

[0035] Figure 2 Schematic diagram of grid spatial distribution. DETAILED DESCRIPTION

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

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

[0038] Embodiment 1:

[0039] like Figure 1 As shown, an embodiment of the present invention provides a gravity field interpolation method suitable for autonomous orbit determination of space-borne GNSS, including:

[0040] Step S1, the server generates sampling points in the grid corresponding to the satellite orbit;

[0041] Step S2: The server obtains the coefficients of the polynomial according to the measured values ​​of the sampling points;

[0042] Step S3: The user terminal converts the ECI coordinates at the reference time into the RTN coordinates corresponding to the grid, and normalizes the RTN coordinates; the normalized RTN coordinates are introduced into the polynomial, and the spherical harmonic gravity field is calculated using the spherical harmonic method. The gravitational acceleration before the term is used to calculate the spherical harmonic gravity field using the polynomial interpolation method The gravitational acceleration after the term.

[0043] As an implementation of the embodiment of the present invention, step S1 includes:

[0044] Step S11: Determine the grid corresponding to the satellite orbit

[0045] Starting from the satellite perigee, the true anomaly angle is taken at every 1 degree, and 360 points are obtained on one orbit. These points are rotated around the z-axis, and the coordinates after each rotation of 1 degree are recorded, and finally a 360×360 grid is obtained. The geometric shape of the grid is a rectangular parallelepiped, and the directions of its three sides are parallel to the three axes of the RTN coordinate system of the current position of the satellite, such as Figure 2 The grid size should ensure that there is no singular area. The lengths in the T and N directions should be as small as possible while ensuring that there is no singular area. The length in the R direction can be taken to a smaller value. The main reason is that when the satellite is in different grids, its R direction is always relatively stable.

[0046] Step S12: Design sampling points within the grid

[0047] According to the Chebyshev node distribution, the sampling points in the grid are designed, and the sampling equation is as follows:

[0048] ;

[0049] in, is the total number of sampling points, and in each dimension points; each dimension is independent of each other, and this location selection helps to minimize the Runge phenomenon and increase the robustness of the fit.

[0050] As an implementation of the embodiment of the present invention, step S2 includes:

[0051] Step S21: Generate candidate polynomials

[0052] The candidate polynomial is a polynomial that may be used for interpolation fitting of the gravity field, and its specific form is:

[0053] ;

[0054] in, , , are the three components of the normalized coordinates of the point to be calculated, are the polynomial coefficients, express An element in is the total accumulated cycle counter, It represents the total degree of the polynomial, and a, b, and c are the variables of the accumulation process.

[0055] Step S22: Get the selected polynomial based on the candidate polynomials

[0056] Since the gravitational acceleration changes slowly in the radial direction, a lower order can well fit the change of acceleration in the radial direction, so the selected polynomial can remove some radial high-order terms to reduce the storage burden.

[0057] by For example: When selecting candidate polynomials, you can consider removing The polynomial of the term That is, the radial direction can choose a lower term fit.

[0058] Step S23: Obtain the coefficients of the polynomial according to the measured values ​​corresponding to the sampling points and the selected polynomial.

[0059] Each sampling point corresponds to a measurement value , solve for the polynomial coefficients as follows:

[0060] ;

[0061] in, is the matrix of Chebyshev sampling points; is a column vector, each element of which corresponds to the measurement value of a sampling point.

[0062] Since there is sampling positions, so there is one sampling points are used to solve the polynomial coefficients. Since the order is , so the matrix Y is a The matrix of N rows and N columns. The coefficients of the polynomial to be solved are is a column vector of N rows, and the least squares formula is used to find The solution is:

[0063] ;

[0064] in, For the matrix Transpose

[0065] Step S24, calculate the gravitational acceleration of the fitting polynomial at the sampling point, compare the calculation results of the spherical harmonic method, obtain the modeling residual and RMS of the sampling point, and verify whether the accuracy within the sampling point meets the standard. If the accuracy of the verification result meets the task requirements, the candidate polynomial can be retained as the interpolation fitting polynomial. If the accuracy is not enough, replace the candidate polynomial and repeat the above steps S21-S23.

[0066] As an implementation of the embodiment of the present invention, step S3 includes:

[0067] Step S31: Update current time With reference time , find the current time Rotation matrix from inertial system to Earth-centered Earth-fixed system and The rotation matrix from the inertial frame to the Earth-centered Earth-fixed coordinate system at the reference time ;

[0068] Step S32: ECI coordinates are transformed using the rotation matrix Convert the coordinates to the Earth-centered Earth-fixed system, and then use the rotation matrix The inverse matrix of is used to transform the coordinates into the ECI coordinate system at the reference time; at this time, the state quantity after coordinate transformation at the reference time is used to calculate the true anomaly angle and the right ascension of the ascending node, and the corresponding grid position in the reference coordinate system is found;

[0069] Step S33: convert the ECI coordinates at the reference time into the RTN coordinates corresponding to the grid, and normalize the RTN coordinates.

[0070] Step S34: Substitute the normalized RTN coordinates into the polynomial and use the spherical harmonic method to calculate the spherical harmonic gravity field. The gravitational acceleration before the term is used to calculate the spherical harmonic gravity field using the polynomial interpolation method The gravitational acceleration after the term.

[0071] Furthermore, the spherical harmonic gravity field for: ;

[0072] in, and Spherical harmonic gravity field The gravitational acceleration before the term, Spherical Harmonic Gravity Field The gravitational acceleration after the term.

[0073] and The analytical solution is:

[0074] ;

[0075] in, is the reference gravitational constant, is the radius of the Earth, and are the magnitude of the position vector and the geocentric latitude, respectively. is a second-order Legendre polynomial.

[0076] Obtained using polynomial fitting.

[0077] The embodiment of the present invention adopts a grid interpolation method to calculate the satellite gravity acceleration, which greatly reduces the time complexity and reduces the satellite computing power requirements for real-time orbit determination. The embodiment of the present invention adopts an adaptive polynomial to ensure that the polynomial meets the accuracy requirements while reducing the coefficient requirements, and reduces the storage occupation compared with other grid interpolation algorithms.

[0078] Embodiment 2:

[0079] The embodiment of the present invention further provides a gravity field interpolation device, comprising:

[0080] The first processing module is used for the server to generate sampling points in the grid corresponding to the satellite orbit;

[0081] The second processing module is used for the server to obtain the coefficients of the polynomial according to the measurement values ​​of the sampling points;

[0082] The third processing module is used for the user end to convert the ECI coordinates at the reference time into the RTN coordinates corresponding to the grid, and normalize the RTN coordinates; the normalized RTN coordinates are brought into the polynomial, and the spherical harmonic gravity field is calculated using the spherical harmonic method. The gravitational acceleration before the term is used to calculate the spherical harmonic gravity field using the polynomial interpolation method The gravitational acceleration after the term.

[0083] As an implementation of the embodiment of the present invention, the first processing module includes:

[0084] A determination unit, used to determine the grid corresponding to the satellite orbit;

[0085] The design unit is used to design the sampling points in the grid according to the Chebyshev node distribution.

[0086] As an implementation of the embodiment of the present invention, the second processing module includes:

[0087] A generating unit, used for generating a candidate polynomial;

[0088] A selection unit, used for obtaining a selected polynomial according to the candidate polynomials;

[0089] A calculation unit, used for obtaining coefficients of the polynomial according to the measured values ​​corresponding to the sampling points and the selected polynomial;

[0090] The inspection unit is used to check whether the internal accuracy of the sampling points meets the standard; if the internal accuracy meets the standard, the coefficients of the polynomial are stored.

[0091] Embodiment 3:

[0092] The present invention also provides a storage medium, on which a computer program is stored, and the computer program executes the gravity field interpolation method when running.

[0093] The embodiments described above are only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A gravity field interpolation method, characterized in that: include: Step S1, the server generates sampling points in the grid corresponding to the satellite orbit; Step S2: The server obtains the coefficients of the polynomial according to the measured values ​​of the sampling points; Step S3: The user terminal converts the ECI coordinates at the reference time into the RTN coordinates corresponding to the grid, and normalizes the RTN coordinates; Substitute the normalized RTN coordinates into the polynomial and use the spherical harmonic method to calculate the spherical harmonic gravity field The gravitational acceleration before the term is used to calculate the spherical harmonic gravity field using the polynomial interpolation method The gravitational acceleration after the term; Step S1 includes: Step S11, determining the grid corresponding to the satellite orbit; Step S12: Design sampling points in the grid according to Chebyshev node distribution ; The sampling equation is as follows: ； in, is the total number of sampling points, taking m points in each dimension; Step S2 includes: Step S21: Generate a candidate polynomial, which is: ； in, , , are the three components of the normalized coordinates of the point to be calculated, are the polynomial coefficients, express An element in is the total cumulative loop counter, d represents the total number of polynomials, and a, b, c are variables of the accumulation process; Step S22, obtaining a selected polynomial according to the candidate polynomials; Step S23, obtaining coefficients of the polynomial according to the measured values ​​corresponding to the sampling points and the selected polynomial; Step S24, verify whether the accuracy of the sampling points meets the standard; if the accuracy of the sampling points meets the standard, store the coefficients of the polynomial; if the required accuracy is not achieved, increase the polynomial order and repeat the above steps S21 to S23; In step S3, the spherical harmonic gravity field for: ; in, and Spherical harmonic gravity field The gravitational acceleration before the term, Spherical Harmonic Gravity Field The gravitational acceleration after the term; and The analytical solution is: ； in, is the reference gravitational constant, is the radius of the Earth, and are the magnitude of the position vector and the geocentric latitude, respectively. is the second-order Legendre polynomial; Obtained using polynomial fitting.

2. A gravity field interpolation device for implementing the gravity field interpolation method according to claim 1, characterized in that: include: The first processing module is used for the server to generate sampling points in the grid corresponding to the satellite orbit; The second processing module is used for the server to obtain the coefficients of the polynomial according to the measurement values ​​of the sampling points; The third processing module is used for the user end to convert the ECI coordinates at the reference time into the RTN coordinates corresponding to the grid and normalize the RTN coordinates; Substitute the normalized RTN coordinates into the polynomial and use the spherical harmonic method to calculate the spherical harmonic gravity field The gravitational acceleration before the term is used to calculate the spherical harmonic gravity field using the polynomial interpolation method The gravitational acceleration after the term.

3. A storage medium, characterized in that: The storage medium stores a computer program, which executes the gravity field interpolation method as claimed in claim 1 when running.

Images