Earth Gravity Calculation Method and Related Device for On-Orbit Real-Time Navigation of Low-Earth Orbit Satellites

Abstract

The present invention discloses a method for calculating the earth's gravity and related devices for real-time on-orbit navigation of low-earth orbit satellites, belonging to the technical field of satellite real-time navigation; the earth's gravity is split into non-spherical gravity with an order ≤ #imgabs0#, central gravity, and non-spherical gravity with an order > #imgabs1#; grid points are evenly distributed globally, and the fitting coefficients of all grid points are calculated and stored; according to the current position of the low-earth orbit satellite, the grid points within the adjacent range are searched, and combined with the fitting coefficients of all grid points, the central gravity and non-spherical gravity with an order > #imgabs2# are calculated; according to the current position of the low-earth orbit satellite, the non-spherical gravity with an order ≤ #imgabs3# is calculated through the recurrence formula of the spherical harmonic model; the non-spherical gravity with an order ≤ #imgabs4#, the central gravity, and the non-spherical gravity with an order > #imgabs5# are added together to obtain the total earth's gravity. The present invention mixes the calculation method using the spherical harmonic model and the calculation method of equivalent center fitting, and balances the calculation efficiency and calculation accuracy at the same time, meeting the requirements of real-time on-orbit navigation processing of low-earth orbit satellites.

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 GNSS L band on the satellite and combined 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 a high computing efficiency on the premise of ensuring orbit accuracy, and of course, it should also take into account a reasonable storage capacity 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 perform integral calculations of the orbit and the state transition matrix. Integration 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 up a large proportion of the computing 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 a relatively low orbit altitude (≤500 km), the traditional spherical harmonic model faces a great contradiction: 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 a reasonable coefficient storage space, which is the key to the engineering implementation of on-orbit real-time navigation of low-earth orbit satellite on-board GNSS.

[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 objective of the present invention is to provide a method and related device for calculating the earth's gravity for on-orbit real-time navigation of low-earth orbit satellites, so as to solve the technical problem in the prior art that it is difficult to balance the calculation accuracy and calculation efficiency in the method for calculating the earth's gravity during satellite navigation.

[0006] To achieve the above objective, the present invention adopts the following technical solutions:

[0007] In a first aspect, the present invention provides a method for calculating the earth's gravity for on-orbit real-time navigation of low-earth orbit satellites, including the following steps:

[0008] 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;

[0009] Uniformly distribute grid points globally, calculate and store the fitting coefficients of all grid points;

[0010] 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 > ;

[0011] 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;

[0012] 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 the real-time navigation of low-earth orbit satellites.

[0013] Furthermore, the step of uniformly distributing grid points globally, calculating and storing the fitting coefficients of all grid points specifically includes:

[0014] Uniformly distribute grid points globally in 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 corresponding to the low-earth orbit satellite at different orbital heights ;

[0015] 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:

[0016]

[0017] In the formula, 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; 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;

[0018] 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];

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

[0020] 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:

[0021] Arrange the 9 fitting coefficients of all grid points globally to form 9 different types of coefficient arrays, obtaining:

[0022]

[0023] 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; 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; Indicates The first fitting coefficient in the corresponding coefficient array; Indicates The second fitting coefficient in the corresponding coefficient array; Indicates The N th fitting coefficient in the corresponding coefficient array;

[0024] 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 occurrences to obtain a sequence, and let there be values in the sequence, and the number of times each value appears is , and ;

[0025] The number of occurrences of the first values with the most occurrences in the sequence is . If , then values are encoded as in the -ary number system, i.e., 0, 1, …, Indicates the maximum integer value that can take when it does not exceed 255;

[0026] If , then there is no need to compress and the compression process stops;

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

[0028] 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:

[0029] Search for grid points within the adjacent range according to the current position of the low-earth orbit satellite, 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;

[0030] 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 , calculate the equivalent center position vector corresponding to each grid point;

[0031] Based on the equivalent center position vector 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;

[0032] According to the equivalent center 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 > of the satellite.

[0033] Furthermore, the step of searching for grid points within the vicinity range according to the current position of the low-earth orbit 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:

[0034] According to the current position of the low-earth orbit satellite, calculate its geocentric longitude , geocentric latitude and geocentric altitude ;

[0035] Search for the nearest grid points around the low-earth orbit satellite, and calculate the spherical search radius centered on the low-earth orbit satellite. The specific calculation formula is:

[0036]

[0037] In the formula, is the spherical search radius; represents the number of grid points in the global range;

[0038] Determine the latitude change range and longitude change range of the grid points:

[0039]

[0040] 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;

[0041] According to the minimum latitude value and the maximum latitude value of the grid points, 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 , the specific calculation formula is:

[0042]

[0043] In the formula, is the ceiling function; is the floor function;

[0044] 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 formula is:

[0045]

[0046] In the formula, is the golden ratio; is the n th geocentric latitude of the alternative grid point; is the n th geocentric longitude of the alternative grid point; is the n th coordinate value of the alternative grid point in the x direction; is the n th coordinate value of the alternative grid point in the y direction; is the n th coordinate value of the alternative grid point in the z direction;

[0047] According to the geocentric longitude and geocentric latitude of the alternative grid points, 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.

[0048] 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 double-independent-variable quadratic surface fitting method based on the equivalent center position vector corresponding to each grid point specifically includes:

[0049] 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 nth grid point when it is at the same geocentric altitude as the satellite;

[0050] Taking the spherical distance from the low Earth orbit satellite and the spherical azimuth angle as the two independent variables, and the equivalent center position vector

[0051]

[0052] In the formula, represents the x component of the equivalent center position vector in the direction; y represents the component of the equivalent center position vector in the z direction; represents 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 ;

[0053] 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:

[0054]

[0055] 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; Indicates the P The polar diameter of the polar coordinates of the grid points required for interpolation; P The polar angle of the polar coordinates of the grid points required for interpolation; Indicates that the first grid point required for interpolation is x The components of the equivalent center position vector of the direction; Indicates that the second interpolation grid point is x The components of the equivalent center position vector of the direction; Indicates the P The grid points required for interpolation are x The components of the equivalent center position vector of the direction; Indicates that the first grid point required for interpolation is y The components of the equivalent center position vector of the direction; Indicates that the second interpolation grid point is y The components of the equivalent center position vector of the direction; Indicates the P The grid points required for interpolation are y The components of the equivalent center position vector of the direction; Indicates that the first grid point required for interpolation is z The components of the equivalent center position vector of the direction; Indicates that the second interpolation grid point is z The components of the equivalent center position vector of the direction; Indicates the P The grid points required for interpolation are z The components of the equivalent center position vector of the direction; represents a matrix; Indicates that all grid points required for interpolation are in x A column vector consisting of the components of the equivalent center position vector of the direction; Indicates that all grid points required for interpolation are in y A column vector consisting of the components of the equivalent center position vector of the direction; Indicates that all grid points required for interpolation are in z A column vector consisting of the components of the equivalent center position vector of the direction;

[0056] Convert the calculation equation of the fitting coefficient into a determinant:

[0057]

[0058] The determinant Convert to upper triangle form , calculate the coefficients of the surface equation , get the equivalent center position vector corresponding to the low-orbit satellite position:

[0059]

[0060] 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 ;

[0061] Furthermore, the step of calculating the non-spherical gravity with order ≤ according to the current position of the low-earth orbit satellite through the recurrence formula of the spherical harmonic model specifically includes:

[0062] Obtain the spherical harmonic coefficients with order ≤ , combine with the current position of the low-earth orbit satellite, and use the recurrence formula of the spherical harmonic model to calculate the spherical harmonic function with order ≤ , calculate the non-spherical gravity of each order, and the specific calculation formula is:

[0063]

[0064] 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

[0065] the component of the gravitational acceleration generated by the gravitational potential of order in the direction; .

[0066] In a second aspect, the present invention provides an earth gravity calculation system for on-orbit real-time navigation of a low-earth orbit satellite, including:

[0067] The gravitational force splitting module is used to split the Earth's gravitational force into non-spherical gravitational forces with orders ≤ , central gravitational force, and non-spherical gravitational forces with orders > ; the is a preset order threshold;

[0068] The grid point layout module is used to uniformly layout grid points globally, calculate and store the fitting coefficients of all grid points;

[0069] The first gravitational force calculation module is used to search for grid points within the adjacent range according to the current position of the low-Earth orbit satellite, and calculate the central gravitational force and non-spherical gravitational forces with orders > by combining the fitting coefficients of all grid points;

[0070] The second gravitational force calculation module is used to calculate the non-spherical gravitational forces with orders ≤ according to the current position of the low-Earth orbit satellite through the spherical harmonic model recurrence formula;

[0071] The total gravitational force calculation module is used to add the non-spherical gravitational forces with orders ≤ , the central gravitational force, and the non-spherical gravitational forces with orders > to obtain the total Earth's gravitational force, and the total Earth's gravitational force is used for real-time navigation of the low-Earth orbit satellite.

[0072] 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 gravitational force for real-time in-orbit navigation of a low-Earth orbit satellite are implemented.

[0073] 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 gravitational force for real-time in-orbit navigation of a low-Earth orbit satellite are implemented.

[0074] Compared with the prior art, the present invention has the following beneficial effects:

[0075] The present invention discloses a method and related device for calculating the earth's gravity for on-orbit real-time 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 convert 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 on-orbit real-time 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 integral of the orbit and 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.

[0076] Further, in the process of obtaining the fitting coefficients of the present invention, the Fibonacci grid arrangement method 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 deployed by the method of equally spacing the grid points by longitude and latitude are relatively sparse in the low-latitude region and relatively dense in the high-latitude region or polar region. The Fibonacci grid arrangement method 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 under the same number of grid points. Similarly, under the same accuracy requirements, fewer grid points are required by using the Fibonacci grid arrangement method. 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.

[0077] Further, the present invention constructs a quadratic surface equation for fitting with the equivalent center position vector as the dependent variable and the spherical polar coordinates of the grid point relative to the position of the satellite as the two independent variables. The constant term of the equation is the equivalent center position vector of the position where the satellite is located. Compared with the conventional inverse distance weighted interpolation method, this fitting and interpolation method more fully considers the difference in gravity in different planar regions of the earth, which is beneficial to improving the calculation accuracy of higher-order gravity. Using the elementary row transformation of the determinant instead of matrix inversion to calculate the constant term of the quadratic surface equation also effectively improves the calculation efficiency. Description of the Drawings

[0078] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings required for 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, without creative efforts, other related drawings can also be obtained based on these drawings.

[0079] Figure 1 It is a flowchart of the method of the present invention;

[0080] Figure 2 It is a schematic diagram of the principle of the system of the present invention;

[0081] Figure 3 It is a schematic diagram of the fitting coefficient generation process of the embodiments of the present invention;

[0082] Figure 4 It is a schematic diagram of the calculation process of the total gravitational force of the earth in the embodiments of the present invention;

[0083] Figure 5 It is a schematic diagram of the on-orbit real-time navigation processing process of low-earth orbit satellites in the embodiments of the present invention;

[0084] Figure 6a It is a comparison chart of the calculation accuracy between the method of the present invention and the traditional spherical harmonic model in the orbit altitude range of 300 km to 330 km in the embodiments of the present invention;

[0085] Figure 6b It is a comparison chart of the calculation efficiency between the method of the present invention and the traditional spherical harmonic model in the orbit altitude range of 300 km to 330 km in the embodiments of the present invention;

[0086] Figure 6c It is a comparison chart of the coefficient storage space between the method of the present invention and the traditional spherical harmonic model in the orbit altitude range of 300 km to 330 km in the embodiments of the present invention;

[0087] Figure 7a It is a schematic diagram of the orbit altitude change of the GRACE-A satellite (Gravity Recovery and Climate Experiment–A) in the embodiments of the present invention;

[0088] Figure 7b It is a comparison chart of the gravitational calculation error when using the method of the present invention and the spherical harmonic models of 55×55 order and 70×70 order in the embodiments of the present invention;

[0089] Figure 7c It is a comparison chart of the gravitational calculation error when using the method of the present invention and the spherical harmonic model of 110×110 order in the embodiments of the present invention;

[0090] Figure 8a This is a comparison graph of the GPS (Global Positioning System) real-time navigation orbit results error when the GRACE-A satellite in the embodiment of the present invention uses the method of the present invention and the spherical harmonic models of order 55×55 and 70×70.

[0091] Figure 8b This is a comparison graph of the GPS real-time navigation orbit results error when the GRACE-A satellite in the embodiment of the present invention uses the method of the present invention and the spherical harmonic model of order 110×110. Detailed implementation manners

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

[0093] 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 art 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 implementation manners according to the present invention.

[0094] See Figure 1 , the embodiment of the present invention discloses a method for calculating the earth's gravity for on-orbit real-time navigation of low-earth orbit satellites, including the following steps:

[0095] 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;

[0096] In the spherical harmonic model, the gravitational potential of the earth on an external space point and the gravitational acceleration can be expressed as:

[0097] (1)

[0098] where 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; is the normalized spherical harmonic coefficient, which is 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 space point. is the position of the spatial point Related functions are:

[0099] (2)

[0100] in, represents the Dirac function; is the gravitational potential corresponding to spherical harmonics of different orders; gravitational potential of different orders By taking the partial derivative of the position of a point in space, we can calculate the corresponding gravitational acceleration at that order. ; Gravitational acceleration of different orders The accumulation of is the total gravitational acceleration.

[0101] According to formula (1), the overall gravitational acceleration can be divided into the central gravity , order ≤ Non-spherical gravity and order Non-spherical gravity:

[0102] (3)

[0103] Order ≤ Non-spherical gravity The calculation requires input order ≤ Normalized spherical harmonic coefficients of , and at the same time, the order ≤ of Carry out recursive calculation. Order ≤ Normalized spherical harmonic coefficients of This is one of the coefficients that needs to be stored on board.

[0104] S2, evenly distributes grid points around the globe, calculates and stores the fitting coefficients of all grid points, such as Figure 3 As shown;

[0105] The central gravity of the space point in formula (3) and order> The non-spherical gravitational equivalent of is expressed as the two-body gravitational force of a particle with the same mass as the Earth at a certain position inside the Earth on a point in space:

[0106] (4)

[0107] Where, Indicates central gravity; Represents a position vector inside the earth;

[0108] This point is called the equivalence center, where The position vector representing the equivalent center. For spatial points at different positions, different position vectors of the equivalent center can be generated.

[0109] Arrange grid points evenly worldwide according to the Fibonacci grid arrangement, and sequentially set serial numbers for each grid point from small to large; let a total of grid points be arranged, then the position and geocentric longitude and latitude of the th grid point are:

[0110] (5)

[0111] 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 <X 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 point positions, with a low coefficient of point positions in low-latitude regions and dense point positions in high-latitude regions, wasting a large number of grid points under the same precision; while the Fibonacci grid arrangement method can ensure that the grid points cover the Earth's surface as evenly as possible.

[0112] Above each grid point, within the preset orbital altitude range, calculate the corresponding position vectors of the equivalent center of the low-Earth orbit satellite at different orbital altitudes ; 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 altitude, within a certain range, the corresponding position vectors of the equivalent center at different geocentric heights can be calculated .

[0113] At each grid point, according to the three components xyz of the position vector of the equivalent center , using the orbital altitude as the independent variable, perform quadratic polynomial fitting to obtain:

[0114] (6)

[0115] In the formula, represents the minimum geocentric height; 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;

[0116] Each grid point corresponds to 9 fitting coefficients , , , , , , , and , and the 9 fitting coefficients are all constrained to be integers in the range of [-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.

[0117] 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, it includes:

[0118] 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:

[0119] (7)

[0120] 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; represents The first fitting coefficient in the corresponding coefficient array; represents The second fitting coefficient in the corresponding coefficient array; Indicates The N th fitting coefficient in the corresponding coefficient array; Indicates The first fitting coefficient in the corresponding coefficient array; Indicates The second fitting coefficient in the corresponding coefficient array; Indicates The N th fitting coefficient in the corresponding coefficient array;

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

[0122] 2) 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 in the sequence, and the number of occurrences of each value is ;

[0123] 3) The first values with the most occurrences in the sequence, that is, the first group values in the sequence, if (where represents the maximum integer value that can take when not exceeding 255), then these values are worth compression storage, and these values can be encoded as the numbers 0, 1,..., in the -ary number system; 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 of the original;

[0124] 4) Sequentially obtain the number of occurrences of the next group of values. If , then repeat the compression process in step 3); if , then it means there is no need for compression, and the compression process stops.

[0125] 5) Store the compressed elements. For the elements that are not compressed, record their positions and values and directly store them using an integer array.

[0126] S3. According to the current position of the low Earth orbit satellite, search for the grid points within the adjacent range, and calculate the central gravity and the non-spherical gravity of the order combining the fitting coefficients of all grid points. of the non-spherical gravity;

[0127] S301. According to the current position of the low Earth orbit satellite, search for the 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.

[0128] 1) Calculate its geocentric longitude according to the current position of the low Earth orbit satellite. and geocentric latitude and geocentric altitude ;

[0129] 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:

[0130] (8)

[0131] In the formula, is the spherical search radius; represents the number of grid points in the global range;

[0132] 3) Determine the latitude change range and longitude change range of the grid points:

[0133] (9)

[0134] In the formula, is the minimum latitude value of the grid points; is the maximum latitude value of the grid points; is the minimum longitude value of the grid points; is the maximum longitude value of the grid points;

[0135] 4) According to the minimum latitude value and the maximum latitude value of the grid points, 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. The specific calculation formula is:

[0136] (10)

[0137] Wherein, is the ceiling function; is the floor function;

[0138] 5) According to the serial number range of the alternative grid points , calculate the geocentric longitude and geocentric latitude of these alternative grid points through formula (5); thereby 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 as the low-earth orbit satellite . 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.

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

[0140] S303. Based on the equivalent center position vectors corresponding to each grid point, adopt 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;

[0141] 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;

[0142] 2) Using the spherical distance and the spherical azimuth angle relative to the low-earth orbit satellite as two independent variables, and the equivalent center position vector as the dependent variable, construct a quadratic surface equation with two independent variables:

[0143] (11)

[0144] Wherein, represents x the equivalent center position vector in the Components of Denote y Equivalent center position vector in the direction; Denote z Equivalent center position vector in the direction; , , , , , , , , , , , , , , , , and are all coefficients of the surface equation; is the rectangular coordinate corresponding to the polar coordinate;

[0145] 3) 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 a calculation equation for the fitting coefficients:

[0146] (12)

[0147] In the formula, denotes the polar radius of the polar coordinate of the first interpolation required grid point; denotes the polar angle of the polar coordinate of the first interpolation required grid point; denotes the polar radius of the polar coordinate of the second interpolation required grid point; denotes the polar angle of the polar coordinate of the second interpolation required grid point; denotes the P polar radius of the polar coordinate of the P th interpolation required grid point; denotes the polar angle of the polar coordinate of the x th interpolation required grid point; [[ID=8 of the equivalent center position vector of the first interpolation required grid point in the direction; x denotes the component of the equivalent center position vector of the second interpolation required grid point in the direction; P denotes the x component of the equivalent center position vector of the th interpolation required grid point in the yComponents of the equivalent central position vector in the direction; Indicates the grid point required for the second interpolation in y Components of the equivalent central position vector in the direction; Indicates the P th grid point required for interpolation in y Components of the equivalent central position vector in the direction; Indicates the grid point required for the first interpolation in z Components of the equivalent central position vector in the direction; Indicates the grid point required for the second interpolation in z Components of the equivalent central position vector in the direction; Indicates the P th grid point required for interpolation in z Components of the equivalent central position vector in the direction; Indicates the matrix; Indicates the column vector composed of the components of the equivalent central position vector of all the grid points required for interpolation in x the direction; Indicates the column vector composed of the components of the equivalent central position vector of all the grid points required for interpolation in y the direction; Indicates the column vector composed of the components of the equivalent central position vector of all the grid points required for interpolation in z the direction;

[0148] 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 ; Substituting into Equation (11), 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, only needs to be calculated. Further transforming Equation (12) to form a determinant:

[0149] (13)

[0150] 5) After multiple elementary row transformations, the determinant is transformed into an 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:

[0151] (14)

[0152] In the formula, indicates the element in the 6th row and 7th column of the upper triangular form ; indicates the upper triangular form The element in the 6th row and 6th column of ; Represents the upper triangle form The element at row 6 and column 8 of ; Represents the upper triangle form The element at row 6 and column 9 of .

[0153] S304, based on the equivalent central position vector corresponding to the position of the low-orbit satellite, use the two-body gravity formula to calculate the central gravity and order of the satellite. Non-spherical gravity.

[0154] S4, according to the current position of the low-orbit satellite, calculate the order ≤ The non-spherical gravity of

[0155] Get order ≤ Spherical harmonic coefficients , combined with the current position of the low-orbit satellite, the spherical harmonic model recursive formula is used to calculate the order ≤ Spherical harmonics of , calculate the non-spherical gravity of each order , the specific calculation formula is:

[0156] (15)

[0157] Where, is the Earth's gravitational constant; is the average radius of the Earth; and Both represent the order of the spherical harmonic model; represents the Dirac function; express The gravitational acceleration generated by the order gravitational potential is x Directional component; express The gravitational acceleration generated by the order gravitational potential is y Directional component; express The gravitational acceleration generated by the order gravitational potential is z Directional component;

[0158] The non-spherical gravity of each order The order obtained by cumulative calculation is ≤ Non-spherical gravity .

[0159] S5, set the order ≤ Non-spherical gravity, central gravity and order> The non-spherical gravitational forces are added to obtain the total gravitational force of the Earth. The total gravitational force of the Earth is used for real-time navigation of low-Earth orbit satellites, enabling time update processing of the real-time navigation filter for low-Earth orbit satellites, and further enabling the calculation of the real-time orbit and clock bias of low-Earth orbit satellites, as Figure 5 shown.

[0160] Embodiment:

[0161] As Figure 6a 、 Figure 6b 、 Figure 6c shown, the method of the present invention is used to calculate the Earth's gravitational force in the orbital altitude range of 300 km to 330 km by using the spherical harmonic models of traditional 55×55 order, 70×70 order, and 110×110 order, and the performance of these methods in terms of calculation accuracy, calculation efficiency, and coefficient storage space is compared. After each step shown in Figure 1 of the present invention, taking the Earth's gravitational force calculated by the spherical harmonic model of EGM (Earth Gravitational Model) 2008 500×500 order as the reference true value, the central gravitational force and non-spherical gravitational forces with orders > 45 are fitted in the interval of 300 km to 330 km, the fitting coefficients are generated and compressed for storage, and at the same time, the spherical harmonic coefficients of non-spherical gravitational forces with orders ≤ 45 are stored to obtain all the coefficients for calculating the Earth's gravitational force of the present invention. After each step shown in Figure 4 of the present invention, the Earth's gravitational force in the orbital altitude range of 300 km to 330 km is calculated. At the same time, taking the Earth's gravitational force calculated by the spherical harmonic model of EGM 2008 500×500 order as the reference true value, the three-dimensional error of the gravitational acceleration is output and RMS (Root Mean Square) statistics are performed as the gravitational calculation accuracy; in addition, the calculation time is also recorded 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 of the method of the present invention is less than that of the spherical harmonic model of 55×55 order, 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.

[0162] As Figure 7a 、 Figure 7b 、 Figure 7cAs shown in the figure, 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 errors 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 smaller than the errors of the gravitational acceleration calculated by the traditional spherical harmonic models of order 55×55, order 70×70, and order 110×110.

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

[0164] Refer to Figure 2 , the embodiment of the present invention discloses a system for calculating the earth's gravity for on-orbit real-time navigation of low-earth orbit satellites, 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:

[0165] Among them, the gravity splitting module is used to 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; the grid point layout module is used to evenly layout grid points globally, calculate and store the fitting coefficients of all grid points; the first gravity calculation module is used to search for grid points within the adjacent range according to the current position of the low-earth orbit satellite, and calculate the central gravity and non-spherical gravity with order > in combination with the fitting coefficients of all grid points; the second gravity calculation module is used to calculate the non-spherical gravity with order ≤ according to the current position of the low-earth orbit satellite through the spherical harmonic model recurrence formula; the total gravity calculation module is used to add the non-spherical gravity with order ≤ , the central gravity, and the non-spherical gravity with order > to obtain the total earth's gravity, and the total earth's gravity is used for the real-time navigation of low-earth orbit satellites.

[0166] 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, and 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 gravity for on-orbit real-time navigation of low-earth orbit satellites.

[0167] 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 and is 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 the operating system of the terminal is stored in this storage space. And, one or more instructions suitable for being loaded and executed by the processor are also stored in this storage space. 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 on-orbit real-time navigation of low-earth orbit satellites in the above embodiment.

[0168] 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 memory, CD-ROM (Compact Disc Read-Only Memory), optical memory, etc.) that contain computer-usable program code.

[0169] The present invention is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to embodiments of the present invention. It should be understood that each flow and / or block in the flowchart and / or block diagram, as well as the combination of flows and / or blocks in the flowchart and / or block diagram, 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, such that the instructions executed by the processor of the computer or other programmable data processing devices generate means for implementing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0170] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing devices to work in a specific manner, such that the instructions stored in the computer-readable memory generate a manufactured article including instruction means that implement the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0171] These computer program instructions can also be loaded onto a computer or other programmable data processing devices, such that a series of operation steps are executed on the computer or other programmable devices to generate a computer-implemented process, and thus the instructions executed on the computer or other programmable devices provide steps for implementing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0172] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not 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: still can modify the specific implementation manners of the present invention or make equivalent replacements, and any modification or equivalent replacement 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 on-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 an order ≤ , central gravity, and non-spherical gravity with an order > ; the is a preset order threshold; Uniformly distribute grid points globally, calculate and store the fitting coefficients of all grid points; Search for grid points within the vicinity based on the current position of the low-earth orbit satellite, and calculate the central gravity and non-spherical gravity of order > using the fitting coefficients of all grid points; According to the current position of the low-earth orbit satellite, calculate the non-spherical gravity with the order ≤ using the recurrence formula of the spherical harmonic model; Add the non-spherical gravitational force, central gravitational force with order ≤ and the non-spherical gravitational force 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 on-orbit navigation of low-earth orbit satellites according to claim 1, wherein The step of uniformly distributing grid points globally, calculating and storing the fitting coefficients of all grid points specifically includes: Uniformly deploy grid points globally in the arrangement according to the Fibonacci grid, 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 corresponding equivalent central position vector of the low-earth orbit satellite at different orbital heights ; 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, wherein, 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, and sort them in descending order of the number of occurrences to obtain a sequence. 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 the values are encoded as the decimal numbers 0, 1, …, ; other values in the sequence are defaulted to 0 for compression; denotes 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 on-orbit real-time navigation of low-earth orbit satellites according to claim 1, characterized in that, Search for grid points within the vicinity based on the current position of the LEO satellite, and calculate the central gravitational force and the non-spherical gravitational force of the order by combining the fitting coefficients of all grid points. The steps specifically include: 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 height of the low-earth orbit satellite , and calculate the equivalent center position vector corresponding to each grid point; Based on the equivalent center position vector corresponding to each grid point, use the double-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 LEO satellite, the central gravitational force and the non-spherical gravitational force with an order greater than acting on the satellite are calculated using the two-body gravitational formula.

5. The method for calculating the earth's gravity for real-time on-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 a low Earth orbit (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: 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. The specific calculation formulas are as follows: Wherein, 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 real-time on-orbit navigation of low-earth orbit satellites according to claim 4, characterized in that, The step of based on the equivalent center position vector corresponding to each grid point, using the double-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 grid points required for interpolation around the position of the LEO satellite be: , where represents the spherical distance between the th grid point required for interpolation 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 central 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-orbit satellite and the spherical azimuth angle as 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: ​ In the formula, 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 the 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 the 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 the 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-earth 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 real-time in-orbit navigation of low-earth orbit satellites according to claim 1, characterized in that, The step of 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 includes: Obtain spherical harmonic coefficients with orders ≤ 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 ≤ Calculate the non-spherical gravity of each order The specific calculation formula is: The specific calculation formula is: In the formula, is the gravitational constant of the Earth; is the average radius of the Earth; and both represent the degree of the spherical harmonic model; represents the Dirac function; represents the component of the gravitational acceleration generated by the gravitational potential of degree x in the direction; the component of the gravitational acceleration generated by the gravitational potential of degree y in the direction; the component of the gravitational acceleration generated by the gravitational potential of degree z in the Accumulate the non-spherical gravity of each order to obtain the non-spherical gravity with order ≤ . .

8. A system for calculating the Earth's gravity for real-time navigation of low-Earth orbit satellites, characterized in that, It includes: A gravitational 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 > ; 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 according to the current position of the low-orbit satellite, and calculate the central gravitational force and the non-spherical gravitational force of the 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 less than or equal to 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 real-time in-orbit 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 real-time in-orbit navigation of a low-earth orbit satellite as described in any one of claims 1-7.

Images