Earth gravity calculation method for on-orbit real-time navigation of low-orbit satellite and related device

Abstract

The invention discloses a terrestrial gravity calculation method for on-orbit real-time navigation of a low-orbit satellite and a related device, and belongs to the technical field of real-time navigation of satellites. The gravitational force is split into non-spherical gravitational force with the order smaller than or equal to # imgabs0 #, central gravitational force and the order gt; the non-spherical gravitation of the # imgabs 1 #; grid points are uniformly distributed in a global range, and fitting coefficients of all the grid points are calculated and stored; searching grid points in an adjacent range according to the current position of the low earth orbit satellite, and calculating a central gravitation and order gt by combining fitting coefficients of all the grid points; a non-spherical gravitational force of # imgabs2 #; according to the current position of the low earth orbit satellite, through a spherical harmonic model recursion formula, calculating non-spherical gravity of which the order is less than or equal to # imgabs3 #; the non-spherical gravitational force, the central gravitational force and the order gt of which the order is smaller than or equal to # imgabs4 #; and adding the non-spherical gravitational forces of # imgabs5 # to obtain the total gravitational force of the earth. According to the method, the calculation mode of the spherical harmonic model and the calculation mode of equivalent center fitting are mixed for use, the calculation efficiency and the calculation precision are balanced at the same time, and the requirement for on-orbit real-time navigation processing of the low-orbit satellite is met.

Description

Technical Field

[0001] The present invention belongs to the technical field of satellite real-time navigation, and relates to a method for calculating the earth's gravity and related devices for on-orbit real-time navigation of low-earth orbit satellites. Background Art

[0002] A low-earth orbit satellite is equipped with a GNSS (Global Navigation Satellite System) receiver. By means of precise measurement in the L band of the on-board GNSS and combining with satellite dynamics information, high-precision orbit parameters of the satellite itself can be obtained in real time. This is the mainstream technical means for low-earth orbit satellites to achieve on-orbit real-time navigation. The on-orbit real-time navigation processing of low-earth orbit satellites is completed on an on-board platform with very limited computing power. This requires that the entire processing algorithm has high computing efficiency on the premise of ensuring orbit accuracy, and of course, reasonable storage capacity should also be taken into account as much as possible.

[0003] When a low-earth orbit satellite conducts on-orbit real-time navigation processing based on on-board GNSS measurement, it is necessary to continuously calculate the integral of the orbit and the state transition matrix. The integral requires multiple calculations of the right function and the forces acting on the satellite. Among the forces acting on the low-earth orbit satellite, the earth's gravity plays an absolute dominant role. In the commonly used gravity field model, the gravitational potential is expressed in the form of spherical harmonic function expansion, which is called the spherical harmonic model. The earth's gravity is calculated according to the partial derivative of the gravitational potential. When using the spherical harmonic model, multiple recursive calculations of the spherical harmonic function are required. Existing research shows that the calculation of the earth's gravity based on the spherical harmonic model takes a large proportion of the calculation time in the on-orbit real-time navigation processing of low-earth orbit satellites. At the same time, the accuracy of gravity calculation directly affects the overall accuracy of the dynamic model and the accuracy of the final real-time navigation result. For low-earth orbit satellites with relatively low orbit altitudes (≤500 km), the traditional spherical harmonic model faces great contradictions: on the one hand, to ensure the accuracy of gravity calculation and not reduce the real-time navigation accuracy, a higher-order spherical harmonic model must be required; on the other hand, the higher the order of the spherical harmonic model, the computing time consumed increases exponentially with the increase of the order of the spherical harmonic model, which is restricted by the limited computing power of the on-board platform. How to make a reasonable balance between calculation accuracy and calculation efficiency, not only meet the accuracy requirements of real-time navigation, but also ensure real-time and efficient operation on the on-board platform, and in addition, take into account reasonable coefficient storage space, which is the key to the engineering implementation of on-orbit real-time navigation of on-board GNSS for low-earth orbit satellites.

[0004] In summary, the traditional spherical harmonic model of the earth's gravity field cannot meet the requirements of on-orbit real-time navigation of satellites with relatively low orbits. There is an urgent need for an optimized method for calculating the earth's gravity that is applicable to on-orbit real-time navigation processing of satellites and takes into account both calculation accuracy and calculation efficiency. Summary of the Invention

[0005] The purpose of the present invention is to provide a method and related device for calculating the earth's gravity for real-time in-orbit navigation of low-earth orbit satellites, so as to solve the technical problem that the existing method for calculating the earth's gravity in satellite navigation is difficult to balance calculation accuracy and calculation efficiency.

[0006] To achieve the above object, the present invention adopts the following technical solutions: In the first aspect, the present invention provides a method for calculating the earth's gravity for real-time in-orbit navigation of low-earth orbit satellites, including the following steps: Split the earth's gravity into non-spherical gravity with order ≤ , central gravity, and non-spherical gravity with order > ; the is a preset order threshold; Uniformly distribute grid points globally, calculate and store the fitting coefficients of all grid points; According to the current position of the low-earth orbit satellite, search for grid points within the adjacent range, and combine the fitting coefficients of all grid points to calculate the central gravity and non-spherical gravity with order > ; According to the current position of the low-earth orbit satellite, calculate the non-spherical gravity with order ≤ through the spherical harmonic model recurrence formula; Add the non-spherical gravity with order ≤ , the central gravity, and the non-spherical gravity with order > to obtain the total earth's gravity, which is used for real-time navigation of low-earth orbit satellites.

[0007] Further, the step of uniformly distributing grid points globally, calculating and storing the fitting coefficients of all grid points specifically includes: Uniformly distribute grid points globally according to the arrangement of Fibonacci grids, and sequentially set serial numbers for each grid point from small to large; at each grid point, within the preset orbital height range, calculate the equivalent central position vector of the low-earth orbit satellite at different orbital heights ; At each grid point, according to the three components xyz of the equivalent central position vector , fit a quadratic polynomial with the orbital height as the independent variable, and the specific calculation formula is:

[0008] In the formula, represents the minimum geocentric height; represents the maximum geocentric height; represents the geocentric height; represents the orbital height; representsx The equivalent center position vector in the direction; denote y the equivalent center position vector in the direction; denote z the equivalent center position vector in the direction; , , , , , , , and are all fitting coefficients; Each grid point corresponds to 9 fitting coefficients , , , , , , , and , and all 9 fitting coefficients are constrained to be integers in the range of [-128, 127]; Arrange the 9 fitting coefficients of all grid points globally to form 9 different types of coefficient arrays, and use a compression algorithm to compress and store the coefficient arrays.

[0009] Furthermore, the step of arranging the 9 fitting coefficients of all grid points globally to form 9 different types of coefficient arrays and using a compression algorithm to compress and store the coefficient arrays specifically includes: Arrange the 9 fitting coefficients of all grid points globally to form 9 different types of coefficient arrays, obtaining:

[0010] wherein, denotes the number of grid points globally; denote the first fitting coefficient in the corresponding coefficient array; denote the second fitting coefficient in the corresponding coefficient array; denote the N th fitting coefficient in the corresponding coefficient array; denote the first fitting coefficient in the corresponding coefficient array; denote the second fitting coefficient in the corresponding coefficient array; denotes the N th fitting coefficient in the corresponding coefficient array; denotes the first fitting coefficient in the corresponding coefficient array; denotes the second fitting coefficient in the corresponding coefficient array; denotes the N th fitting coefficient in the corresponding coefficient array; denotes the first fitting coefficient in the corresponding coefficient array; denotes the second fitting coefficient in the corresponding coefficient array; denotes the N th fitting coefficient in the corresponding coefficient array; denotes the first fitting coefficient in the corresponding coefficient array; denotes the second fitting coefficient in the corresponding coefficient array; denotes the N th fitting coefficient in the corresponding coefficient array; denotes the first fitting coefficient in the corresponding coefficient array; denotes the second fitting coefficient in the corresponding coefficient array; denotes the N th fitting coefficient in the corresponding coefficient array; denotes the first fitting coefficient in the corresponding coefficient array; denotes the second fitting coefficient in the corresponding coefficient array; denotes the N th fitting coefficient in the corresponding coefficient array; denotes the first fitting coefficient in the corresponding coefficient array; denotes the second fitting coefficient in the corresponding coefficient array; denotes the N th fitting coefficient in the corresponding coefficient array; denotes the first fitting coefficient in the corresponding coefficient array; denotes the second fitting coefficient in the corresponding coefficient array; Indicates the N th fitting coefficient in the corresponding coefficient array; Statistically analyze all value types of the elements in each coefficient array and the number of times each value appears, sort them in descending order of the number of occurrences to obtain a sequence, and let there be values in the sequence, the number of times each value appears is , and ; For the top values with the most occurrences in the sequence, the number of occurrences is . If , then encode the values into the binary numbers 0, 1, …, in the -ary number system; default the other values in the sequence to 0 and perform compression; Indicates the maximum integer value that can take when it does not exceed 255; If , then there is no need to compress and the compression process stops; Store the compressed elements. For the elements that did not participate in the compression, record their positions and values and directly store them using an integer array.

[0011] Furthermore, the step of searching for grid points within the adjacent range according to the current position of the low-earth orbit satellite, combining the fitting coefficients of all grid points, and calculating the central gravitational force and the non-spherical gravitational force of order specifically includes: According to the current position of the low-earth orbit satellite, search for grid points within the adjacent range, obtain the serial numbers and geocentric longitude and latitude of the grid points within the adjacent range, and determine the grid points required for interpolation; According to the serial numbers of the grid points required for interpolation, use the decompression algorithm to obtain the fitting coefficients of the grid points required for interpolation, and combine the geocentric altitude where the low-earth orbit satellite is located to calculate the equivalent central position vector corresponding to each grid point; Based on the equivalent central position vectors corresponding to each grid point, use the double-independent-variable quadratic surface fitting method to interpolate and calculate the equivalent central position vector corresponding to the position where the low-earth orbit satellite is located; According to the equivalent central position vector corresponding to the position where the low-earth orbit satellite is located, use the two-body gravitational formula to calculate the central gravitational force and the non-spherical gravitational force of order acting on the satellite.

[0012] Further, the step of searching for grid points within the vicinity range based on the current position of the LEO satellite, obtaining the serial numbers and geocentric longitude and latitude of the grid points within the vicinity range, and determining the grid points required for interpolation specifically includes: Calculate its geocentric longitude based on the current position of the LEO satellite , geocentric latitude and geocentric height ; Search for the nearest grid points around the LEO satellite and calculate the spherical search radius centered on the LEO satellite , and the specific calculation formula is:

[0013] In the formula, is the spherical search radius; represents the number of grid points within the global range; Determine the latitude change range and longitude change range :

[0014] In the formula, is the minimum latitude value of the grid point; is the maximum latitude value of the grid point; is the minimum longitude value of the grid point; is the maximum longitude value of the grid point; According to the minimum latitude value and the maximum latitude value of the grid point, obtain the minimum serial number and the maximum serial number of the alternative grid points within the latitude range, and obtain the serial number range of the alternative grid points, and the specific calculation formula is:

[0015] In the formula, is the ceiling function; is the floor function; According to the serial number range of the alternative grid points, calculate the geocentric longitude and geocentric latitude of the alternative grid points, and the specific calculation formula is:

[0016] In the formula, is the golden ratio; is the geocentric latitude of the n th alternative grid point; is the nThe geocentric longitude of the For the n th alternative grid point in x direction coordinate value; For the n th alternative grid point in y direction coordinate value; For the n th alternative grid point in z direction coordinate value; According to the geocentric longitude and geocentric latitude of the alternative grid point, calculate the spherical distance and azimuth angle between the alternative grid point and the low-earth orbit satellite when they are at the same geocentric height ; if the spherical distance is less than the spherical search radius , it means that this alternative grid point is the grid point required for interpolation.

[0017] Furthermore, the step of interpolating and calculating the equivalent center position vector corresponding to the position of the low-earth orbit satellite by using the quadratic surface fitting method with two independent variables based on the equivalent center position vector corresponding to each grid point specifically includes: Let the basic information of the interpolation required grid points around the position of the low-earth orbit satellite be: , where represents the spherical distance between the th interpolation required grid point and the low-earth orbit satellite when they are at the same geocentric height, represents the azimuth angle between the th grid point and the low-earth orbit satellite when they are at the same geocentric height, represents the equivalent center position vector corresponding to the th grid point when it is at the same geocentric height as the satellite; Taking the spherical distance relative to the low-earth orbit satellite and the spherical azimuth angle as two independent variables, and the equivalent center position vector as the dependent variable, construct a quadratic surface equation with two independent variables:

[0018] In the formula, represents the component of the equivalent center position vector x in the direction; represents the component of the equivalent center position vector y in the direction; represents the component of the equivalent center position vector z in the direction; , , , , , , , , , , , , , , , , and are all coefficients of the surface equation; is the polar coordinate corresponding rectangular coordinate; Substitute the basic information of the interpolation required grid points around the position of the low-earth orbit satellite into the quadratic surface equation with two independent variables to establish the calculation equation of the fitting coefficients:

[0019] In the formula, represents the polar radius of the polar coordinate of the first interpolation required grid point; represents the polar angle of the polar coordinate of the first interpolation required grid point; represents the polar radius of the polar coordinate of the second interpolation required grid point; represents the polar angle of the polar coordinate of the second interpolation required grid point; represents the P polar radius of the polar coordinate of the P th interpolation required grid point; represents the polar angle of the polar coordinate of the first interpolation required grid point in the x direction of the equivalent center position vector component; represents the component of the equivalent center position vector of the second interpolation required grid point in the x direction; represents the P th interpolation required grid point in the x direction of the equivalent center position vector component; represents the component of the equivalent center position vector of the first interpolation required grid point in the y direction; represents the component of the equivalent center position vector of the second interpolation required grid point in the y direction; represents the P th interpolation required grid point in the y direction of the equivalent center position vector component; Denote the component of the equivalent center position vector of the first interpolation required grid point in the z direction; Denote the component of the equivalent center position vector of the second interpolation required grid point in the z direction; Denote the P th component of the equivalent center position vector of the interpolation required grid point in the z direction; Denote the matrix; Denote the column vector composed of the components of the equivalent center position vectors of all interpolation required grid points in the x direction; Denote the column vector composed of the components of the equivalent center position vectors of all interpolation required grid points in the y direction; Denote the column vector composed of the components of the equivalent center position vectors of all interpolation required grid points in the z direction; Convert the calculation equation of the fitting coefficient into a determinant:

[0020] For the determinant Convert it to upper triangular form , calculate the coefficients of the surface equation , and obtain the equivalent center position vector corresponding to the position of the low-earth orbit satellite:

[0021] In the formula, Denote the element in the 6th row and 7th column of the upper triangular form ; Denote the element in the 6th row and 6th column of the upper triangular form ; Denote the element in the 6th row and 8th column of the upper triangular form ; Denote the element in the 6th row and 9th column of the upper triangular form ;

[0022] Furthermore, the step of calculating the non-spherical gravity with order ≤ according to the current position of the low-earth orbit satellite through the spherical harmonic model recurrence formula specifically includes: Obtain the spherical harmonic coefficients with order ≤ , combine with the current position of the low-earth orbit satellite, and use the spherical harmonic model recurrence formula to calculate the spherical harmonic function with order ≤ , calculate the non-spherical gravity of each order, and the specific calculation formula is:

[0023] In the formula, is the gravitational constant of the earth; is the average radius of the earth; and both represent the order of the spherical harmonic model; represents the Dirac function; represents the component of the gravitational acceleration generated by the gravitational potential of order x in the direction; represents y the component of the gravitational acceleration generated by the gravitational potential of order in the direction; z represents the component of the gravitational acceleration generated by the gravitational potential of order in the direction; .

[0024] In a second aspect, the present invention provides a system for calculating the earth's gravity for real-time in-orbit navigation of a low-earth orbit satellite, including: A gravity splitting module, configured to split the earth's gravity into non-spherical gravity with an order ≤ , central gravity, and non-spherical gravity with an order > ; the is a preset order threshold; A grid point layout module, configured to uniformly layout grid points globally, calculate and store the fitting coefficients of all grid points; A first gravity calculation module, configured to search for grid points within a neighboring range according to the current position of the low-earth orbit satellite, and calculate the central gravity and non-spherical gravity with an order > in combination with the fitting coefficients of all grid points; A second gravity calculation module, configured to calculate the non-spherical gravity with an order ≤ according to the current position of the low-earth orbit satellite through the spherical harmonic model recurrence formula; A total gravity calculation module, configured to add the non-spherical gravity with an order ≤ , the central gravity, and the non-spherical gravity with an order > to obtain the total earth's gravity, and the total earth's gravity is used for real-time navigation of the low-earth orbit satellite.

[0025] In a third aspect, the present invention provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the steps of the method for calculating the earth's gravity for real-time on-orbit navigation of low-earth orbit satellites are implemented.

[0026] In a fourth aspect, the present invention provides a computer-readable storage medium storing a computer program, and when the computer program is executed by a processor, the steps of the method for calculating the earth's gravity for real-time on-orbit navigation of low-earth orbit satellites are implemented.

[0027] Compared with the prior art, the present invention has the following beneficial effects: The present invention discloses a method and related device for calculating the earth's gravity for real-time on-orbit navigation of low-earth orbit satellites. First, the earth's gravity is split into lower-order non-spherical gravity, central gravity, and higher-order non-spherical gravity. For the lower-order non-spherical gravity, since the number of recursion times is small, the spherical harmonic model recursion formula is still used for calculation to ensure its calculation accuracy. For the higher-order non-spherical gravity and central gravity, since the required number of recursion times is large, the equivalent center fitting method is used to transform it into the calculation method of two-body gravity, thereby effectively improving the calculation efficiency. The two methods are used in combination to balance the requirements of calculation efficiency and calculation accuracy at the same time. The present invention ensures the calculation accuracy of the earth's gravity while effectively improving the calculation efficiency of the earth's gravity. In the process of real-time on-orbit navigation processing of low-earth orbit satellites, using the method of the present invention to calculate the earth's gravity can ensure the high-efficiency and high-precision calculation of the integration of the orbit and the state transition matrix, thereby ensuring the real-time and efficient processing of the entire real-time navigation algorithm on the spaceborne platform, and at the same time being able to improve or at least not reduce the result accuracy of real-time navigation.

[0028] Furthermore, in the process of obtaining the fitting coefficients of the present invention, the Fibonacci grid arrangement is used to deploy grid points globally, rather than the traditional method of equally spacing the grid points by longitude and latitude. The grid points arranged by the equally spacing method of longitude and latitude are sparse in the low-latitude region and dense in the high-latitude region or polar region. The Fibonacci grid arrangement can deploy grid points as evenly as possible globally, thereby realizing the balanced expression of gravity in different regions of the world and achieving higher interpolation accuracy with the same number of grid points. Similarly, with the same accuracy requirement, fewer grid points are required using the Fibonacci grid arrangement. At the same time, the fitting coefficients at each grid point are represented by an integer array and compressed for storage, thereby minimizing the storage space required for the coefficients on the satellite.

[0029] Furthermore, in the present invention, the equivalent center position vector is used as the dependent variable, and the spherical polar coordinates of the grid points relative to the position of the satellite are used as the two independent variables to construct a quadratic surface equation for fitting. The constant term of the equation is the equivalent center position vector of the satellite's location. Compared with the conventional inverse distance weighted interpolation method, this fitting and interpolation method takes into account the differences in gravity in different planar regions of the earth more fully, which is beneficial to improving the calculation accuracy of high-order gravity. Calculating the constant term of the quadratic surface equation by using elementary row transformation of the determinant instead of matrix inversion also effectively improves the calculation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as limiting the scope. For those of ordinary skill in the art, other related drawings can be obtained based on these drawings without creative efforts.

[0031] Figure 1 is a flowchart of the method of the present invention; Figure 2 is a schematic diagram of the principle of the system of the present invention; Figure 3 is a schematic diagram of the fitting coefficient generation process of the embodiment of the present invention; Figure 4 is a schematic diagram of the calculation process of the total gravity of the earth in the embodiment of the present invention; Figure 5 is a schematic diagram of the on-orbit real-time navigation processing process of a low-orbit satellite in the embodiment of the present invention; Figure 6a is a comparison diagram of the calculation accuracy between the method of the present invention and the traditional spherical harmonic model in the orbital altitude range of 300 km to 330 km in the embodiment of the present invention; Figure 6b is a comparison diagram of the calculation efficiency between the method of the present invention and the traditional spherical harmonic model in the orbital altitude range of 300 km to 330 km in the embodiment of the present invention; Figure 6c is a comparison diagram of the coefficient storage space between the method of the present invention and the traditional spherical harmonic model in the orbital altitude range of 300 km to 330 km in the embodiment of the present invention; Figure 7a is a schematic diagram of the orbital altitude change of the GRACE-A satellite (Gravity Recovery and Climate Experiment–A) in the embodiment of the present invention; Figure 7bGravitational calculation error comparison diagram when the method of the present invention is used with spherical harmonic models of order 55×55 and 70×70 in the embodiments of the present invention; Figure 7c Gravitational calculation error comparison diagram when the method of the present invention is used with a spherical harmonic model of order 110×110 in the embodiments of the present invention; Figure 8a GPS (Global Positioning System) real-time navigation orbit result error comparison diagram when the GRACE-A satellite in the embodiments of the present invention uses the method of the present invention with spherical harmonic models of order 55×55 and 70×70; Figure 8b GPS real-time navigation orbit result error comparison diagram when the GRACE-A satellite in the embodiments of the present invention uses the method of the present invention with a spherical harmonic model of order 110×110. Detailed implementation manners

[0032] The present invention will be described in detail below with reference to the drawings and in conjunction with embodiments. It should be noted that, without conflict, the embodiments in the present application and the features in the embodiments can be combined with each other.

[0033] The following detailed descriptions are all exemplary descriptions, aiming to provide further detailed descriptions of the present invention. Unless otherwise specified, all technical terms adopted by the present invention have the same meaning as commonly understood by those of ordinary skill in the technical field to which the present application belongs. The terms used in the present invention are only for describing specific implementation manners and are not intended to limit the exemplary embodiments according to the present invention.

[0034] See Figure 1 , the embodiments of the present invention disclose a method for calculating the earth's gravity for on-orbit real-time navigation of low-earth orbit satellites, including the following steps: S1, splitting the earth's gravity into non-spherical gravity with order ≤ , central gravity, and non-spherical gravity with order > ; the is a preset order threshold; In the spherical harmonic model, the gravitational potential of the earth on an external space point and the gravitational acceleration can be expressed as: (1) Among them, is the earth's gravitational constant, and both represent the order of the spherical harmonic model, is the vector length of the space point, is the geocentric latitude and geocentric longitude of the space point, is the average radius of the earth; are the normalized spherical harmonic coefficients given by the corresponding gravity field model. These coefficients depend on the mass distribution inside the Earth and are independent of the position of the spatial point; is a function related to the position of the spatial point and has: (2) where represents the Dirac function; is the gravitational potential corresponding to spherical harmonic terms of different orders; the gravitational potentials of different orders are differentiated with respect to the position of the spatial point to calculate the corresponding gravitational acceleration at this order ; the summation of the gravitational accelerations of different orders is the total gravitational acceleration.

[0035] According to Equation (1), the total gravitational acceleration can be decomposed into the central gravity , the non-spherical gravity with order ≤ and the non-spherical gravity with order > : (3) The non-spherical gravity with order ≤ requires the input of the normalized spherical harmonic coefficients with order ≤ . At the same time, a recursive calculation needs to be carried out for with order ≤ . The normalized spherical harmonic coefficients with order ≤ are one of the coefficients that need to be stored on the satellite.

[0036] S2, uniformly distribute grid points globally, calculate and store the fitting coefficients of all grid points, as shown in Figure 3 ; The central gravity of the spatial point in Equation (3) and the non-spherical gravity with order > are equivalently represented as the two-body gravity of a mass point with the same mass as the Earth at a certain position inside the Earth for the spatial point: (4) In the formula, represents the central gravity; represents a position vector inside the Earth; This point is called the equivalent center, where represents the position vector of the equivalent center. For spatial points at different positions, different equivalent center position vectors can be generated.

[0037] The grid points are evenly arranged globally according to the Fibonacci grid, and each grid point is sequentially numbered from small to large; let there be a total of grid points, then the position and geocentric longitude and latitude of the th grid point are: (5) where is the golden ratio; is the geocentric latitude of the n th alternative grid point; is the geocentric longitude of the n th alternative grid point; is the coordinate value of the n th alternative grid point in the x direction; is the coordinate value of the n th alternative grid point in the y direction; is the coordinate value of the n th alternative grid point in the z direction. The conventional equidistant arrangement of longitude and latitude will lead to uneven distribution of points, with a low point coefficient in low-latitude regions and dense points in high-latitude regions, wasting a large number of grid points under the same precision; while the Fibonacci grid layout can ensure that the grid points cover the earth's surface as evenly as possible.

[0038] Above each grid point, within the preset orbital height range, calculate the equivalent center position vector corresponding to different orbital heights of the low-earth orbit satellite; for a certain spatial point on the earth, fixing its planar position (geocentric longitude and geocentric latitude), when the low-earth orbit satellite is at this spatial point and continuously changes its geocentric height, that is, the orbital height, within a certain range, the equivalent center position vector corresponding to different geocentric heights can be calculated.

[0039] At each grid point, according to the three components of the equivalent center position vector xyz , using a quadratic polynomial for fitting with the orbital height as the independent variable, we get: (6) In the formula, where represents the minimum geocentric height; represents the maximum geocentric height; represents the geocentric height; represents the orbital height; represents x the component of the equivalent center position vector in the direction;y The equivalent center position vector in the direction component; denote z The equivalent center position vector in the direction component; , , , , , , , and are all fitting coefficients; Each grid point corresponds to 9 fitting coefficients , , , , , , , and , and the 9 fitting coefficients are all constrained to be integers in the interval [-128, 127]; thus ensuring that each coefficient can be stored in an integer variable that only occupies 1 byte, so as to minimize the on-board storage space of these coefficients.

[0040] Arrange the 9 fitting coefficients of all grid points globally to form 9 different types of coefficient arrays, and use a compression algorithm to compress and store the coefficient arrays. Specifically include: 1) Layout grid points, each grid point has 9 fitting coefficients, and each fitting coefficient is represented by an integer variable of 1 byte, then a total of bytes of storage space are required. For further compressed storage, first, the 9 fitting coefficients of all grid points globally can form 9 different integer arrays: (7) In the formula, represents the number of grid points globally; denote the first fitting coefficient in the corresponding coefficient array; denote the second fitting coefficient in the corresponding coefficient array; denote the N th fitting coefficient in the corresponding coefficient array; denote the first fitting coefficient in the corresponding coefficient array; denote the second fitting coefficient in the corresponding coefficient array; denote The N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents The N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents The N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents The N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents The N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents The N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents The N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents the N th fitting coefficient in the corresponding coefficient array; For each coefficient array, compression storage can be performed separately. The value range of the elements in each coefficient array is [-128, 127], and the number of elements is , which is generally much larger than 256, meaning that many element values in the array will appear frequently. Based on this element repeatability, a compression algorithm can be designed.

[0041] 2) Count all the value types of the elements in each coefficient array and the number of times each value appears, sort them in descending order of the number of appearances to obtain a sequence, and let there be a total of values in the sequence, and the number of times each value appears is ; 3) The top values with the most appearances in the sequence, that is, the first group values appearance times , if (where represents the maximum integer value that can take when not exceeding 255), then these values have the value of compression storage, that is, the values can be encoded as in the number system of the number system 0, 1,..., ; and the other values in the sequence can be defaulted to 0. Thus, each value in the sequence can be compressed into an integer not exceeding 255, meaning that the storage space becomes 4) Obtain the appearance times of the next group of values in sequence. If , then repeat the compression process in step 3); if , it means that no compression is required and the compression process stops.

[0042] 5) Store the compressed elements. For the elements that do not participate in compression, record their positions and values, and directly store them using an integer array.

[0043] S3. According to the current position of the low-earth orbit satellite, search for the grid points within the adjacent range, and combine the fitting coefficients of all grid points to calculate the central gravitational force and the non-spherical gravitational force of the order ; S301. Search for grid points within the vicinity based on the current position of the low Earth orbit satellite, obtain the sequence numbers and geocentric longitude and latitude of the grid points within the vicinity, and determine the grid points required for interpolation. 1) Calculate its geocentric longitude based on the current position of the low Earth orbit satellite , geocentric latitude and geocentric height . 2) Search for the nearest grid points around, and calculate the spherical search radius (in radians) centered on the low Earth orbit satellite. The specific calculation formula is: (8) In the formula, is the spherical search radius; represents the number of grid points globally; 3) Determine the latitude change range and longitude change range of the grid points: (9) In the formula, is the minimum latitude value of the grid point; is the maximum latitude value of the grid point; is the minimum longitude value of the grid point; is the maximum longitude value of the grid point; 4) Based on the minimum latitude value and the maximum latitude value of the grid points, obtain the minimum sequence number and the maximum sequence number of the alternative grid points within the latitude range, and get the sequence number range of the alternative grid points. The specific calculation formula is: (10) In the formula, is the ceiling function; is the floor function; 5) Based on the sequence number range of the alternative grid points, calculate the geocentric longitude and geocentric latitude of these alternative grid points through formula (5); and further calculate the spherical distance and azimuth angle between the alternative grid points and the low Earth orbit satellite when they are at the same geocentric height . If the spherical distance is less than the spherical search radius , it means that this alternative grid point is the grid point required for interpolation.

[0044] S302. According to the serial numbers of the grid points required for interpolation, use the decompression algorithm to obtain the fitting coefficients of the grid points required for interpolation, and combine the geocentric height at which the low-earth orbit satellite is located , and calculate the equivalent center position vector corresponding to each grid point.

[0045] S303. Based on the equivalent center position vectors corresponding to each grid point, use the quadratic surface fitting method with two independent variables to interpolate and calculate the equivalent center position vector corresponding to the position of the low-earth orbit satellite; 1) Let the basic information of the grid points required for interpolation around the position of the low-earth orbit satellite be: , where represents the spherical distance between the -th grid point required for interpolation and the low-earth orbit satellite when they are at the same geocentric height, represents the azimuth angle between the -th grid point and the low-earth orbit satellite when they are at the same geocentric height, represents the equivalent center position vector corresponding to the -th grid point when it is at the same geocentric height as the satellite; 2) Using the spherical distance relative to the low-earth orbit satellite and the spherical azimuth angle as the two independent variables, and the equivalent center position vector as the dependent variable, construct a quadratic surface equation with two independent variables: (11) In the formula, represents the component of the equivalent center position vector x in the direction; represents the component of the equivalent center position vector y in the direction; represents the component of the equivalent center position vector z in the direction; , , , , , , , , , , , , , , , , and All are the coefficients of the surface equation; is the polar coordinate The corresponding rectangular coordinate; 3) Substitute the basic information of the grid points required for interpolation around the position of the low-earth orbit satellite into the quadratic surface equation with two independent variables to establish an equation for calculating the fitting coefficients: (12) In the formula, represents the radial distance of the polar coordinate of the first grid point required for interpolation; represents the polar angle of the polar coordinate of the first grid point required for interpolation; represents the radial distance of the polar coordinate of the second grid point required for interpolation; represents the polar angle of the polar coordinate of the second grid point required for interpolation; represents the P radial distance of the polar coordinate of the P th grid point required for interpolation; represents the polar angle of the polar coordinate of the x th grid point required for interpolation; represents the component of the equivalent central position vector of the first grid point required for interpolation in the x direction; represents the component of the equivalent central position vector of the second grid point required for interpolation in the P th direction; x represents the component of the equivalent central position vector of the th grid point required for interpolation in the y direction; represents the component of the equivalent central position vector of the first grid point required for interpolation in the y direction; represents the component of the equivalent central position vector of the second grid point required for interpolation in the P th direction; y represents the component of the equivalent central position vector of the th grid point required for interpolation in the z direction; represents the component of the equivalent central position vector of the first grid point required for interpolation in the z direction; represents the component of the equivalent central position vector of the second grid point required for interpolation in the P th direction; z represents the component of the equivalent central position vector of the th grid point required for interpolation in the direction; represents the matrix; x represents the column vector composed of the components of the equivalent central position vectors of all grid points required for interpolation in the Denote the column vector composed of the equivalent central position vector components of all the grid points required for interpolation in the y direction; Denote the column vector composed of the equivalent central position vector components of all the grid points required for interpolation in the z direction; 4) In the spherical polar coordinate system with the low-earth orbit satellite as the origin, the spherical distance and spherical azimuth angle of the low-earth orbit satellite are ; Substitute it into Equation (11), and the equivalent central position vector can be calculated as . Therefore, when calculating the coefficients of the quadratic surface equation according to Equation (12), it is not necessary to calculate all the coefficients, and only needs to be calculated. Further transform Equation (12) to form a determinant: (13) 5) After multiple elementary row transformations, transform the determinant into the upper triangular form , and the coefficients of the surface equation can be directly calculated, that is, the corresponding pseudo-central position vector of the low-earth orbit satellite is obtained: (14) In the formula, represents the element in the 6th row and 7th column of the upper triangular form ; represents the element in the 6th row and 6th column of the upper triangular form ; represents the element in the 6th row and 8th column of the upper triangular form ; represents the element in the 6th row and 9th column of the upper triangular form .

[0046] S304. According to the equivalent central position vector corresponding to the position of the low-earth orbit satellite, use the two-body gravitational formula to calculate the central gravitational force and the non-spherical gravitational force with order > acting on the satellite.

[0047] S4. According to the current position of the low-earth orbit satellite, use the spherical harmonic model recurrence formula to calculate the non-spherical gravitational force with order ≤ ; Obtain the spherical harmonic coefficients with order ≤ , combine the current position of the low-earth orbit satellite, and use the spherical harmonic model recurrence formula to calculate the spherical harmonic function with order ≤ , calculate the non-spherical gravitational force for each order, and the specific calculation formula is: (15) In the formula, is the gravitational constant of the Earth; is the average radius of the Earth; and both represent the order of the spherical harmonic model; represents the Dirac function; represents the component of the gravitational acceleration generated by the gravitational potential of order x in the direction; represents y the component of the gravitational acceleration generated by the gravitational potential of order in the direction; z represents the component of the gravitational acceleration generated by the gravitational potential of order in the direction; .

[0048] S5, add the non-spherical gravity of order ≤ , the central gravity, and the non-spherical gravity of order > to obtain the total Earth gravity, which is used for real-time navigation of low-orbit satellites, can realize the time update process of the real-time navigation filter of low-orbit satellites, and further realize the solution of the real-time orbit and clock error of low-orbit satellites, as Figure 5 shown.

[0049] Example: As Figure 6a , Figure 6b , Figure 6c shown, using the method of the present invention and the traditional spherical harmonic models of 55×55 order, 70×70 order, and 110×110 order to calculate the Earth gravity in the orbit height range of 300 km to 330 km, and compare the performance of these methods in terms of calculation accuracy, calculation efficiency, and coefficient storage space. After each step shown in Figure 1 of the present invention, taking the Earth gravity calculated by the spherical harmonic model of EGM (Earth Gravitational Model) 2008 500×500 order as the reference true value, fitting the central gravity and the non-spherical gravity of order > 45 in the range of 300 km to 330 km to generate fitting coefficients and compress and store them, and at the same time store the spherical harmonic coefficients of the non-spherical gravity of order ≤ 45 to obtain all the coefficients for calculating the Earth gravity of the present invention. After Figure 4For each of the steps shown, calculate the Earth's gravity within the orbital altitude range of 300 km to 330 km. At the same time, taking the Earth's gravity calculated by the spherical harmonic model of EGM 2008 at the 500×500 order as the reference true value, output the three-dimensional error of the gravitational acceleration, and perform RMS (Root Mean Square) statistics, which is used as the gravitational calculation accuracy. In addition, record the calculation time consumed to evaluate the calculation efficiency. From the perspective of calculation accuracy, the method of the present invention is not only superior to the spherical harmonic models of 55×55 order and 70×70 order, but also superior to the spherical harmonic model of 110×110 order, with the highest calculation accuracy. From the perspective of calculation efficiency, the calculation time consumed by the method of the present invention is less than that of the 55×55 order spherical harmonic model, with the least calculation time. In addition, considering the required storage space, although the coefficient storage space required by the method of the present invention is more than that of the spherical harmonic model, it is still only 381.3 Kb, within the range that can be borne by the spaceborne platform.

[0050] As Figure 7a , Figure 7b , Figure 7c shown, it is a comparison chart of the orbital altitude change of the GRACE-A satellite in the embodiment of the present invention on October 23, 2017 and the comparison of the truncation error in gravitational calculation using the method of the present invention and the traditional spherical harmonic model respectively. It can be clearly seen that the error of the gravitational acceleration calculated by the method of the present invention is significantly less than the error of the gravitational acceleration calculated by the traditional spherical harmonic models of 55×55 order, 70×70 order and 110×110 order.

[0051] As Figure 8a , Figure 8b shown, it is a comparison chart of the three-dimensional error of the orbital results of the GRACE-A satellite in the embodiment of the present invention on October 23, 2017 when performing spaceborne GPS real-time navigation processing using the method of the present invention and the traditional spherical harmonic model respectively. It can be clearly seen that the real-time orbital error obtained by using the method of the present invention is significantly less than the real-time orbital error at 55×55 order and 70×70 order, and is also slightly less than the real-time orbital error at 110×110 order.

[0052] See Figure 2 , the embodiment of the present invention discloses a system for calculating the Earth's gravity for on-orbit real-time navigation of a low-Earth orbit satellite, including a gravity splitting module, a grid point layout module, a first gravity calculation module, a second gravity calculation module and a total gravity calculation module: Among them, the gravity splitting module is used to split the Earth's gravity into non-spherical gravity with an order ≤ , central gravity and non-spherical gravity with an order > ; the is a preset order threshold; a grid point layout module, configured to evenly layout grid points globally, calculate and store the fitting coefficients of all grid points; a first gravitational force calculation module, configured to search for grid points within a neighboring range according to the current position of the low-earth orbit satellite, and calculate the central gravitational force and the non-spherical gravitational force with an order > of the non-spherical gravitational force; a second gravitational force calculation module, configured to calculate the non-spherical gravitational force with an order ≤ according to the current position of the low-earth orbit satellite through the recurrence formula of the spherical harmonic model; a total gravitational force calculation module, configured to add the non-spherical gravitational force with an order ≤ of the non-spherical gravitational force, the central gravitational force and the non-spherical gravitational force with an order > of the non-spherical gravitational force to obtain the total earth gravitational force, and the total earth gravitational force is used for real-time navigation of the low-earth orbit satellite.

[0053] In one embodiment of the present invention, a computer device is provided. The computer device includes a processor and a memory. The memory is used to store a computer program. The computer program includes program instructions. The processor is used to execute the program instructions stored in the computer storage medium. The processor may be a central processing unit (CPU), or may also be other general-purpose processors, digital signal processors (DSPs), application specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. It is the computing core and control core of the terminal, and is suitable for implementing one or more instructions. Specifically, it is suitable for loading and executing one or more instructions in the computer storage medium to implement the corresponding method flow or corresponding function; the processor described in the embodiment of the present invention can be used for the operation of the method for calculating the earth's gravitational force for real-time navigation of a low-earth orbit satellite.

[0054] The present invention also provides a storage medium, specifically a computer-readable storage medium (Memory). The computer-readable storage medium is a memory device in a computer device, used to store programs and data. It can be understood that the computer-readable storage medium here can include both the built-in storage medium in the computer device and, of course, the extended storage medium supported by the computer device. The computer-readable storage medium provides a storage space, and this storage space stores the operating system of the terminal. Moreover, one or more instructions suitable for being loaded and executed by the processor are stored in this storage space, and these instructions can be one or more computer programs (including program codes). It should be noted that the computer-readable storage medium here can be high-speed RAM (Random Access Memory), or non-volatile memory, such as at least one disk memory. One or more instructions stored in the computer-readable storage medium can be loaded and executed by the processor to implement the corresponding steps of the method for calculating the earth's gravity for real-time in-orbit navigation of low-earth orbit satellites in the above embodiments.

[0055] Those skilled in the art should understand that the embodiments of the present invention can be provided as a method, a system, or a computer program product. Therefore, the present invention can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk memories, CD-ROMs (Compact Disc Read-Only Memories), optical memories, etc.) containing computer-usable program codes.

[0056] The present invention is described with reference to the flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to the embodiments of the present invention. It should be understood that each flow and / or block in the flowcharts and / or block diagrams, as well as the combination of flows and / or blocks in the flowcharts and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate means for implementing the functions specified in one Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0057] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to work in a particular manner, such that the instructions stored in the computer-readable memory produce a manufacture including an instruction device that implements the functions specified in one or more of the procedures Figure 1 one or more procedures and / or blocks Figure 1 specified in one or more of the blocks.

[0058] These computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process, so that the instructions executed on the computer or other programmable apparatus provide steps for implementing the functions specified in one or more of the procedures Figure 1 one or more procedures and / or blocks Figure 1 specified in one or more of the blocks.

[0059] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the above embodiments, those of ordinary skill in the art should understand that: modifications or equivalent substitutions can still be made to the specific embodiments of the present invention. Any modification or equivalent substitution that does not depart from the spirit and scope of the present invention shall be covered by the protection scope of the claims of the present invention.

Claims

1. A method for calculating the earth's gravity for real-time in-orbit navigation of low-earth orbit satellites, characterized in that, It includes the following steps: Split the earth's gravity into non-spherical gravity with order ≤ , central gravity, and non-spherical gravity with order > ; the is a preset order threshold; Uniformly distribute grid points globally, calculate and store the fitting coefficients of all grid points; According to the current position of the low-earth orbit satellite, search for grid points within the adjacent range, and combine the fitting coefficients of all grid points to calculate the central gravity and the non-spherical gravity of order > of the spherical harmonic coefficients; According to the current position of the low-earth orbit satellite, calculate the non-spherical gravity with order ≤ by using the recurrence formula of the spherical harmonic model; Add the non-spherical gravity with order ≤ , the central gravity, and the non-spherical gravity with order > to obtain the total gravitational force of the Earth, which is used for real-time navigation of low-orbit satellites.

2. The method for calculating the earth's gravity for real-time in-orbit navigation of low-earth orbit satellites according to claim 1, characterized in that, The step of uniformly distributing grid points globally, calculating and storing the fitting coefficients of all grid points specifically includes: The grid points are evenly arranged globally according to the Fibonacci grid, and serial numbers are sequentially set for each grid point from small to large; at each grid point, within the preset orbit height range, the equivalent central position vectors corresponding to the low-earth orbit satellites at different orbit heights are calculated ; At each grid point, according to the equivalent center position vector of xyz the three components , taking the orbital altitude as the independent variable, a quadratic polynomial is fitted, and the specific calculation formula is: In the formula, represents the minimum geocentric altitude; represents the maximum geocentric altitude; represents the geocentric altitude; represents the orbital altitude; represents x the component of the equivalent center position vector in the direction; represents y the component of the equivalent center position vector in the direction; represents z the component of the equivalent center position vector in the direction; , , , , , , , and are all fitting coefficients; Each grid point corresponds to nine fitting coefficients 、 、 、 、 、 、 、 and , and all nine fitting coefficients are constrained to be integers within the range of [-128, 127]; Arrange the 9 fitting coefficients of all grid points globally to form 9 different types of coefficient arrays, and use a compression algorithm to compress and store the coefficient arrays.

3. The method for calculating the earth's gravity for real-time in-orbit navigation of low-earth orbit satellites according to claim 2, characterized in that, The step of arranging the 9 fitting coefficients of all grid points globally to form 9 different types of coefficient arrays and using a compression algorithm to compress and store the coefficient arrays specifically includes: Arrange the 9 fitting coefficients of all grid points globally to form 9 different types of coefficient arrays, and obtain: In the formula, represents the number of grid points globally; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents the N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents the N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents the N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents the N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents the N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents the N th fitting coefficient in the corresponding coefficient array; represents the first fitting coefficient in the corresponding coefficient array; represents the second fitting coefficient in the corresponding coefficient array; represents the N th fitting coefficient in the corresponding coefficient array; denotes the first fitting coefficient in the corresponding coefficient array; denotes the second fitting coefficient in the corresponding coefficient array; denotes the N th fitting coefficient in the corresponding coefficient array; denotes the first fitting coefficient in the corresponding coefficient array; denotes the second fitting coefficient in the corresponding coefficient array; denotes the N th fitting coefficient in the corresponding coefficient array; Count all the value types of the elements in each coefficient array and the number of occurrences of each value, sort them in descending order of the number of occurrences to obtain a sequence, and let there be a total of values , the number of occurrences of each value is , and ; The occurrence count of the top values that appear most frequently in the sequence is . If , then values are encoded as the decimal numbers 0, 1, …, in the decimal number system; the other values in the sequence are defaulted to 0 for compression; represents the maximum integer value that can take when it does not exceed 255; If , there is no need for compression and the compression process stops; Store the compressed elements. For the elements that do not participate in compression, record their positions and values and directly store them using an integer array.

4. The method for calculating the earth's gravity for real-time in-orbit navigation of low-earth orbit satellites according to claim 1, characterized in that, The step of searching for grid points within a neighboring range according to the current position of the LEO satellite, and calculating the central gravity and the non-spherical gravity of an order greater than specifically includes: According to the current position of the low-earth orbit satellite, search for grid points within the adjacent range, obtain the serial numbers and geocentric longitude and latitude of the grid points within the adjacent range, and determine the grid points required for interpolation; According to the serial numbers of the grid points required for interpolation, the decompression algorithm is used to obtain the fitting coefficients of the grid points required for interpolation, and combined with the geocentric height of the low-earth orbit satellite , calculate the equivalent center position vector corresponding to each grid point; Based on the equivalent center position vectors corresponding to each grid point, use the double-independent-variable quadratic surface fitting method to interpolate and calculate the equivalent center position vector corresponding to the position of the low-earth orbit satellite; According to the equivalent central position vector corresponding to the position of the low-earth orbit satellite, use the two-body gravitational formula to calculate the central gravitational force and the non-spherical gravitational force with an order greater than that the satellite is subjected to.

5. The method for calculating the earth's gravity for real-time in-orbit navigation of low-earth orbit satellites according to claim 4, characterized in that, The step of according to the current position of the low-earth orbit satellite, searching for grid points within the adjacent range, obtaining the serial numbers and geocentric longitude and latitude of the grid points within the adjacent range, and determining the grid points required for interpolation specifically includes: Calculate the geocentric longitude, geocentric latitude, and geocentric altitude of the LEO satellite based on its current position; Search for the nearest grid points around the LEO satellite and calculate the spherical search radius centered on the LEO satellite , and the specific calculation formula is: In the formula, is the spherical search radius; represents the number of grid points globally; Determine the latitude variation range of grid points and the longitude variation range : In the formula, is the minimum latitude value of the grid point; is the maximum latitude value of the grid point; is the minimum longitude value of the grid point; is the maximum longitude value of the grid point; According to the minimum latitude value of the grid points and the maximum latitude value of the grid points , obtain the minimum serial number of the alternative grid points within the latitude range and the maximum serial number , to obtain the serial number range of the alternative grid points , the specific calculation formula is: wherein is the ceiling function; is the floor function; According to the serial number range of the alternative grid points , calculate the geocentric longitude and geocentric latitude of the alternative grid points. The specific calculation formulas are as follows: In the formula, is the golden ratio; is the geocentric latitude of the n th alternative grid point; is the geocentric longitude of the n th alternative grid point; is the coordinate value of the n th alternative grid point in the x direction; is the coordinate value of the n th alternative grid point in the y direction; is the coordinate value of the n th alternative grid point in the z direction; According to the geocentric longitude and geocentric latitude of the alternative grid points, calculate the spherical distance and azimuth angle between the alternative grid points and the low-earth orbit satellite when they are at the same geocentric altitude ; if the spherical distance is less than the spherical search radius , it means that the alternative grid point is the grid point required for interpolation.

6. The method for calculating the earth's gravity for on-orbit real-time navigation of low-earth orbit satellites according to claim 4, wherein, The step of based on the equivalent center position vectors corresponding to each grid point, using the double-independent-variable quadratic surface fitting method to interpolate and calculate the equivalent center position vector corresponding to the position of the low-earth orbit satellite specifically includes: Let the basic information of the interpolation required grid points around the position of the LEO satellite be: , where represents the spherical distance between the -th interpolation required grid point and the LEO satellite when they are at the same geocentric altitude, represents the azimuth angle between the -th grid point and the LEO satellite when they are at the same geocentric altitude, represents the equivalent center position vector corresponding to the -th grid point when it is at the same geocentric altitude as the satellite; In terms of the spherical distance relative to the low-earth orbit satellite and the spherical azimuth angle as the two independent variables, and the equivalent center position vector as the dependent variable, a quadratic surface equation with two independent variables is constructed as follows: In the formula, represents x the component of the equivalent center position vector in the direction; represents y the component of the equivalent center position vector in the direction; represents z the component of the equivalent center position vector in the direction; , , , , , , , , , , , , , , , , and are all coefficients of the surface equation; is the rectangular coordinate corresponding to the polar coordinate ; Substitute the basic information of the grid points required for interpolation around the position of the low-earth orbit satellite into the quadratic surface equation with two independent variables to establish a calculation equation for the fitting coefficients: ​ Wherein, represents the radial distance of the polar coordinates of the first grid point required for interpolation; represents the polar angle of the polar coordinates of the first grid point required for interpolation; represents the radial distance of the polar coordinates of the second grid point required for interpolation; represents the polar angle of the polar coordinates of the second grid point required for interpolation; represents the P radial distance of the polar coordinates of the P th grid point required for interpolation; represents the polar angle of the polar coordinates of the x th grid point required for interpolation; represents the component of the equivalent center position vector of the first grid point required for interpolation in the x direction; represents the component of the equivalent center position vector of the second grid point required for interpolation in the P direction; x represents the component of the equivalent center position vector of the th grid point required for interpolation in the y direction; represents the component of the equivalent center position vector of the first grid point required for interpolation in the y direction; represents the component of the equivalent center position vector of the second grid point required for interpolation in the P direction; y represents the component of the equivalent center position vector of the th grid point required for interpolation in the z direction; represents the component of the equivalent center position vector of the first grid point required for interpolation in the z direction; represents the component of the equivalent center position vector of the second grid point required for interpolation in the P direction; z represents the component of the equivalent center position vector of the th grid point required for interpolation in the direction; represents the matrix; x represents the column vector composed of the components of the equivalent center position vectors of all grid points required for interpolation in the direction; y represents the column vector composed of the components of the equivalent center position vectors of all grid points required for interpolation in the direction; z represents the column vector composed of the components of the equivalent center position vectors of all grid points required for interpolation in the Convert the calculation equation of the fitting coefficient into a determinant: Convert the determinant to upper triangular form , calculate the coefficients of the surface equation , and obtain the equivalent central position vector corresponding to the position of the low-orbit satellite: In the formula, represents the element in the 6th row and 7th column of the upper triangular form ; represents the element in the 6th row and 6th column of the upper triangular form ; represents the element in the 6th row and 8th column of the upper triangular form ; represents the element in the 6th row and 9th column of the upper triangular form .

7. The method for calculating the earth's gravity for on-orbit real-time navigation of low-earth orbit satellites according to claim 1, wherein, Calculating the non-spherical gravity with an order less than or equal to according to the current position of the LEO satellite through the recurrence formula of the spherical harmonic model, specifically including: Obtain spherical harmonic coefficients with orders ≤ , and combined with the current position of the low-earth orbit satellite, use the recurrence formula of the spherical harmonic model to calculate spherical harmonic functions with orders ≤ , and calculate the non-spherical gravity of each order . The specific calculation formula is: , and calculate the non-spherical gravity of each order . The specific calculation formula is as follows: In the formula, is the Earth's gravitational constant; is the average radius of the Earth; and All represent the order of the spherical harmonic model; represents the Dirac function; express The gravitational acceleration produced by the gravitational potential of order x Directional component; express The gravitational acceleration produced by the gravitational potential of order y Directional component; express The gravitational acceleration produced by the gravitational potential of order z Directional component; Accumulatively calculate the non-spherical gravity of each order to obtain the non-spherical gravity with the order ≤ . .

8. The earth's gravity calculation system for on-orbit real-time navigation of low-earth orbit satellites, wherein, It includes: A gravitational splitting module for splitting the Earth's gravity into non-spherical gravity with an order ≤ , central gravity, and non-spherical gravity with an order > ; where the is a preset order threshold; A grid point distribution module for uniformly distributing grid points globally, calculating and storing the fitting coefficients of all grid points; The first gravitational force calculation module is used to search for grid points within the adjacent range based on the current position of the low-orbit satellite, and calculate the central gravitational force and the non-spherical gravitational force of order > by combining the fitting coefficients of all grid points; The second gravitational force calculation module is used to calculate the non-spherical gravitational force with an order ≤ based on the current position of the low-earth orbit satellite through the recurrence formula of the spherical harmonic model; The total gravitational force calculation module is used to add the non-spherical gravitational force with an order ≤ , the central gravitational force, and the non-spherical gravitational force with an order > to obtain the total gravitational force of the Earth, and the total gravitational force of the Earth is used for real-time navigation of low-orbit satellites.

9. A computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the method for calculating the earth's gravity for on-orbit real-time navigation of a low-earth orbit satellite as described in any one of claims 1-7.

10. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by the processor, it implements the steps of the method for calculating the earth's gravity for on-orbit real-time navigation of a low-earth orbit satellite as described in any one of claims 1-7.

Images