A fast prediction method for lunar orbit rendezvous guidance

Abstract

This invention provides a rapid prediction method for lunar orbit rendezvous guidance, achieving high-precision and rapid prediction to meet mission requirements. The prediction method of this invention applies the gravitational field acceleration calculation method used for GPS autonomous orbit determination of Earth orbit satellites to the field of lunar orbit rendezvous. Addressing the problems that the algorithm cannot be directly applied to the design of lunar orbit rendezvous guidance strategies, improvements and upgrades are made one by one, resulting in a high-precision and rapid algorithm suitable for lunar orbit rendezvous flight control missions to meet mission requirements. Specifically, this invention is the first to transplant this algorithm to the field of lunar orbit applications, which is unprecedented. Considering the special characteristics of the lunar gravitational field, where the acceleration measures of various gravitational perturbations are of similar order and cannot be approximated by the J2 term alone, the J22 term perturbation acceleration is added to the algorithm, thereby significantly improving the computational accuracy and numerical stability of the algorithm.

Description

Technical Field

[0001] This invention relates to the field of deep space exploration orbit design technology, and specifically to a rapid prediction method for lunar orbit rendezvous guidance. Background Technology

[0002] Near-circular lunar orbits at altitudes of 100km to 200km are the most commonly used orbit type in lunar exploration missions. The primary perturbation force affecting these orbits is the non-spherical gravitational field of the Moon. Gravitational perturbation acceleration calculation models typically use spherical harmonic functions, and each acceleration calculation requires a set of Legendre polynomial recursive formulas. The computational complexity of these recursions increases geometrically with the order of the gravitational field. Due to the inhomogeneity of the Moon's gravitational field, it lacks a J2 principal term similar to those of Earth or Mars. To ensure computational accuracy, higher-order gravitational field models are usually introduced for calculations.

[0003] Engineering experience shows that for orbit design, a 32x32 order is typically required to ensure the accuracy of the design results; while for orbit determination and flight control maneuvering strategies, a gravitational field order of at least 100x100 is needed to meet mission requirements. Increasing the gravitational field order leads to a geometric increase in computational load, resulting in a significant increase in computation time. During the flight control implementation phase of the lunar orbit rendezvous and docking remote guidance and maneuvering strategy, the online design, generation, and injection of the orbit control strategy need to be completed within a short time based on real-time orbit determination results. Due to the large number of control variables, long integration time, and high accuracy requirements of the gravitational field model, traditional numerical integration methods, due to their slow computation speed, are difficult to adapt to the requirements of high-precision, rapid strategy reconstruction design under emergency conditions. Summary of the Invention

[0004] In view of this, the present invention provides a rapid prediction method for lunar orbit rendezvous guidance, which can achieve high-precision and rapid prediction and meet the requirements of the mission.

[0005] To achieve the above objectives, the technical solution of the present invention is as follows:

[0006] The present invention provides a rapid prediction method for lunar orbit rendezvous guidance. First, the gravitational field acceleration is pruned to exclude the central gravitational field and the J2 term. 22The residual acceleration after gravitational acceleration from the residual gravitational field is converted into an equivalent "pseudo-center" parameter, thus nonlinearly mapping the residual gravitational field acceleration into a "pseudo-center". Then, within a given longitude and latitude point and a specified altitude range, polynomial interpolation is used to calculate the interpolation polynomial coefficients of the "pseudo-center" function with respect to altitude within that altitude range. The lunar surface is then divided into several grid points according to longitude and latitude, and the fitting polynomial coefficients of the "pseudo-center" function for each longitude and latitude point within a given altitude range are calculated and stored. Finally, for any spatial location point, the pseudo-center corresponding to that point is obtained by calculating the "pseudo-center" corresponding to the instantaneous altitude of its six nearest longitude and latitude grid points through two-dimensional interpolation, and the true gravitational field acceleration corresponding to that point is calculated through nonlinear mapping. The partial derivative matrix of the acceleration vector with respect to the position vector is calculated. Based on the true gravitational field acceleration and the partial derivative matrix of the acceleration vector with respect to the position vector, the position and velocity of the probe are calculated.

[0007] This includes the following steps:

[0008] Step 1: Given the initial position and velocity of the probe; calculate the acceleration and its partial derivatives for the two-body problem in the lunar central gravitational field;

[0009] Step 2, calculate the corresponding J2 terms and J 22 The term perturbation acceleration and partial derivative;

[0010] Step 3, calculate the sum of terms J2 and J2 after removing them. 22 The pseudo-center of the item;

[0011] Step 4: Decompose the three-dimensional space into three directions: longitude, latitude, and altitude. Discretize the pseudo-center to obtain a polynomial fit for the pseudo-center in the altitude direction.

[0012] Step 5: Based on the given detector position, obtain discrete grid points of height h, longitude and latitude; perform two-dimensional six-point interpolation on the two-dimensional discrete grid points of height h, longitude and latitude to obtain the pseudo-center interpolation results of the corresponding grid points;

[0013] Step 6: Calculate the corresponding gravitational acceleration based on the interpolation results, and add the J2 term and J... 22 The effect of the term ultimately yields the true acceleration;

[0014] Step 8: Based on the partial derivatives of the calculated state transition matrix, obtain the position, velocity, and state transition matrix at any given time. This completes the solution to the orbit prediction problem.

[0015] In step 1, the acceleration 'a' in the two-body problem of the central gravitational field 2B Represented as:

[0016]

[0017] Where μ is the Earth's gravitational field constant, and x, y, and z are the three directional components of the detector's position vector, r = (xy z), where r is the magnitude of the position vector.

[0018] a 2B The partial derivative with respect to the position vector is:

[0019]

[0020] In step 2, the perturbation acceleration of term J2 is:

[0021]

[0022] Where J2 is the coefficient of the second-order band harmonic term, R e The radius of the Earth's equator;

[0023] The partial derivative of acceleration in the J2 term is:

[0024]

[0025]

[0026]

[0027]

[0028] J 22 The perturbation acceleration of the term is

[0029]

[0030] Among them, J 22 λ is the second-order quadratic harmonic coefficient. 22 The geographical longitude of the principal axis of the second-order harmonic term, i.e., the azimuth angle of Changzhou on the equatorial ellipse, is a constant; Λ 22 For intermediate computational quantities, Λ 22 =xcosλ 22 +ysinλ 22 ;

[0031] J 22 The partial derivative of the perturbation acceleration with respect to the position vector is

[0032]

[0033]

[0034]

[0035]

[0036] In step 3, the mapping between any point position r and gravitational acceleration a0 in the lunar-centered fixed coordinate system and the pseudo-center c0 is as follows: Where r is the actual distance from the lunar center of the probe, representing the pseudo-lunar center distance ρ plus the pseudo-center c; μ is the Earth's gravitational field constant; a0 is the distance from the lunar center of the probe after removing the pseudo-center c. and The subsequent perturbation gravitational acceleration satisfies the following equation:

[0037]

[0038] Among them, a CBF It is the gravitational acceleration of the probe relative to the lunar-centered fixed coordinate system, calculated using a non-spherical gravitational field acceleration function. It includes the central gravitational field acceleration and perturbation acceleration, where CBF represents the central celestial fixed coordinate system; the pseudo-center c and the perturbation acceleration a... CBF The relationship is

[0039] Beneficial effects:

[0040] 1. The prediction method of this invention applies the gravitational field acceleration calculation method used for GPS autonomous orbit determination of Earth orbit satellites to the field of lunar orbit rendezvous. Addressing the problem that the algorithm cannot be directly applied to the design of lunar orbit rendezvous guidance strategies, improvements and upgrades are made one by one, resulting in a high-precision, fast algorithm suitable for lunar orbit rendezvous flight control missions to meet mission requirements. Specifically, this invention is the first to transplant the algorithm to the lunar orbit application field, which is unprecedented. Considering the special characteristics of the lunar gravitational field, where the acceleration measures of various gravitational perturbations are of similar order and cannot be approximated by the J2 term alone, the J22 term perturbation acceleration is added to the algorithm, thereby significantly improving the calculation accuracy and numerical stability of the algorithm.

[0041] 2. To address the problem that the original algorithm could not calculate the sensitivity matrix and therefore could not be used for iterative correction of guidance strategies, this invention adds an interpolation calculation method for the partial derivative of acceleration to the original algorithm, thereby meeting the needs of lunar orbit rendezvous guidance law design.

[0042] 3. The method of this invention addresses the needs of remote guidance missions for lunar orbit rendezvous and docking. For the first time, it introduces lunar gravitational field calculation into the GAAF algorithm, which is traditionally used to calculate gravitational acceleration in the Earth's gravitational field. Furthermore, it improves and upgrades the algorithm to address the differences between the lunar and Earth's gravitational fields and the requirements for calculating lunar orbit guidance strategies. Without sacrificing accuracy, it can effectively increase the predicted velocity of the lunar orbit by up to two orders of magnitude, meeting the requirements of real-time on-orbit mission planning. It is particularly suitable for the requirement of rapid reconfiguration of orbit control strategies under emergency conditions. Attached Figure Description

[0043] Figure 1 This is a flowchart of a method according to an embodiment of the present invention.

[0044] Figure 2 This is a schematic diagram of the two-dimensional six-point interpolation method of the present invention. Detailed Implementation

[0045] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0046] This invention provides a rapid prediction method for lunar orbit rendezvous guidance. Addressing the needs of long-range lunar orbit rendezvous and docking guidance missions, it introduces lunar gravitational field calculations for the first time into the traditional GAAF algorithm for calculating gravitational acceleration in the Earth's gravitational field. Furthermore, the algorithm is improved and upgraded to account for the differences between the lunar and Earth's gravitational fields and the requirements for calculating lunar orbit guidance strategies. This is achieved by using J... 22 The gravitational field coefficients are introduced into the GAAF algorithm, and the formula for the difference of partial derivatives of GAAF gravitational acceleration is derived for rapid calculation of the state transition matrix, enabling high-precision and rapid prediction of detector position and velocity. The specific process of this invention is as follows:

[0047] First, remove the central gravitational field and the J2 term from the gravitational field acceleration. 22 The residual acceleration after the residual gravitational field acceleration is converted into an equivalent "pseudo-center" parameter, thus nonlinearly mapping the residual gravitational field acceleration into a "pseudo-center".

[0048] Then, within a given longitude and latitude point and a specified altitude range, the interpolation polynomial coefficients of the "pseudo-center" as a function of altitude within that altitude range are calculated using polynomial interpolation.

[0049] The lunar surface is then divided into several grid points according to longitude and latitude, and the fitting polynomial coefficients of the "pseudo-center" function of each longitude and latitude point within a given altitude range are calculated and saved.

[0050] Finally, for any spatial location, the pseudo-center corresponding to the instantaneous height of the six nearby latitude and longitude grid points is calculated by two-dimensional interpolation. The actual gravitational field acceleration corresponding to the point is then calculated by nonlinear mapping. The partial derivative matrix of the acceleration vector with respect to the position vector is calculated. Based on the actual gravitational field acceleration and the partial derivative matrix of the acceleration vector with respect to the position vector, the position and velocity of the detector are calculated.

[0051] Specifically, the principle of the method of the present invention is as follows:

[0052] Given the spacecraft's position and velocity (r0, v0) at time t0 and the spacecraft's acceleration function a 3×1After (r,v,t), the problem of predicting the lunar orbit of a spacecraft is defined as the initial value problem of the following 42-dimensional first-order ordinary differential equation system:

[0053]

[0054] Where (r, v) represents the spacecraft's position and velocity at time t, Φ 6×6 (t,t0) is the state transition matrix from time t0 to time t.

[0055] The initial values ​​of the system of equations are:

[0056]

[0057] The acceleration experienced by a lunar orbiter mainly includes the lunar gravitational field acceleration α. CBF Earth's three-body perturbation acceleration a Earth3B The solar three-body perturbation acceleration a Sun3B and solar radiation pressure acceleration a SRP ,Right now:

[0058] a 3×1 (r,v,t)=a CBF +a Earth3B +a Sun3B +a SRP (3)

[0059] For a low-Earth orbit spacecraft in a near-circular lunar orbit, the lunar gravitational acceleration a CBF and its partial derivatives (Gravitational acceleration depends only on the position vector r, therefore) The calculation of higher-order Legendre polynomials is a major factor affecting computational speed. This invention provides the gravitational field acceleration a. CBF and the corresponding partial derivative matrix A fast interpolation calculation method is proposed. Considering the unique characteristics of the lunar gravitational field, this paper proposes two main terms for its perturbation acceleration: the J2 term (second-order zonal harmonic term) for gravitational perturbation acceleration. and J 22 Gravitational acceleration (second-order quadratic harmonic term) and the corresponding partial derivatives of acceleration and The remaining accelerations a0 are given using an analytical expression method:

[0060]

[0061] After mapping it to a pseudo-center, it is discretized in the altitude, latitude, and longitude directions. The altitude discretization uses polynomial fitting coefficients for storage, while the latitude and longitude discretization uses equally spaced latitude and longitude discrete points. When calculating the true acceleration, the two main acceleration terms are first calculated analytically. and Then, the residual acceleration a0 is obtained through the interpolation discrete data inversion method. By adding the two, the total gravitational perturbation acceleration a can be calculated. CBF The corresponding partial derivatives of acceleration The same approach is used, which will not be elaborated here. The calculated a... CBF and Substituting into equation (1), under the initial value (2), orbit prediction can be performed. The calculation speed of orbit prediction is mainly determined by the calculation speed of acceleration in the right function. Since the a0 part is processed by numerical interpolation, the calculation speed is greatly improved, thus greatly improving the speed of orbit prediction. Moreover, since the interpolation calculation retains the characteristics of all higher-order perturbation accelerations, the algorithm will not reduce the calculation accuracy due to the speed improvement brought by interpolation. Actual tests show that compared with the traditional orbit prediction model of 100x100 order gravitational field, the orbit prediction position error of this algorithm is less than 100m, the velocity error is less than 0.01m / s, the state transition matrix error is less than 0.1%, and the calculation speed is improved by nearly two orders of magnitude.

[0062] The flowchart of the orbit prediction method in this embodiment is as follows: Figure 1 As shown, it includes the following steps:

[0063] Step 1, give the detector the initial position and velocity Calculate the acceleration and its partial derivatives in the two-body problem of the lunar central gravitational field. Wherein, the acceleration *a* in the two-body problem of the central gravitational field is... 2B It can be represented as:

[0064]

[0065] Where μ is the Earth's gravitational field constant, and x, y, and z are the three directional components of the detector's position vector, r = (xy z), where r is the magnitude of the position vector.

[0066] a 2B The partial derivative with respect to the position vector is:

[0067]

[0068] Step 2, calculate the corresponding J2 terms and J 22 The perturbation acceleration and partial derivatives of the term; where the perturbation acceleration of the J2 term is:

[0069]

[0070] Where J2 is the coefficient of the second-order band harmonic term, R e This is the radius of the Earth's equator.

[0071] The partial derivative of acceleration in the J2 term is:

[0072]

[0073]

[0074]

[0075]

[0076] J 22 The perturbation acceleration of the term is

[0077]

[0078] Among them, J 22 λ is the second-order quadratic harmonic coefficient. 22 The geographical longitude of the principal axis of the second-order harmonic term, i.e., the azimuth angle of Changzhou on the equatorial ellipse, is a constant; Λ 22 For intermediate computational quantities, Λ 22 =xcosλ 22 +ysinλ 22 .

[0079] J 22 The partial derivative of the perturbation acceleration with respect to the position vector is

[0080]

[0081]

[0082]

[0083]

[0084] Step 3, calculate the values ​​after removing J2 and J 22 The pseudo-center of the item;

[0085] Among them, gravitational acceleration can be expressed as:

[0086]

[0087] Among them, a CBFThis is the gravitational acceleration of the probe relative to the lunar-centered fixed coordinate system, calculated using a non-spherical gravitational field acceleration function. It includes the central gravitational field acceleration and perturbation acceleration. CBF represents the central fixed coordinate system. ρ is defined as the probe's relative gravitational acceleration a. CBF The pseudo-lunar distance (Pseudo-Radius) is ρ = |ρ|.

[0088] The actual lunar center distance r of the probe is defined as the pseudo-lunar center distance plus the pseudo-center (c):

[0089] c = r - ρ (18)

[0090] pseudo-center c and perturbation acceleration a CBF The relationship is:

[0091]

[0092] Calculation removed and The subsequent perturbation gravitational acceleration a0 is

[0093]

[0094] Furthermore, the mapping between any point position r and gravitational acceleration a0 in the lunar-centered fixed coordinate system and the pseudo-center c0 (c0 is the pseudo-center corresponding to a0) in three-dimensional space can be obtained as follows:

[0095]

[0096] Step 4: Decompose the three-dimensional space into three directions: longitude, latitude, and altitude. Discretize the pseudo-center to obtain a polynomial fitting of the pseudo-center in the altitude direction. The specific process is as follows:

[0097] In the longitude and latitude directions, the sphere is discretized into several grid points. For a given grid point coordinate... Its pseudo-center in the height direction is c0 = [c x c y c z Fitting was performed using rational polynomials:

[0098]

[0099] i = x, y, z represent the x, y, and z directions of the pseudo-center vector, D k Let N be the denominator polynomial. k The numerator polynomial is a1…a n-1 b1…b n-1 All coefficients are polynomials. The fitting method uses the common Pade approximation function.

[0100]

[0101] H max and H min These represent the fitted height ranges. For a low-orbit, near-circular lunar orbit, H is typically taken as... min =10km, H max =350km can meet the requirements of most tasks. Different values ​​of m and n will affect the fitting accuracy and calculation speed. Through extensive calculations and testing, n=7 and m=6 can achieve a good trade-off between storage capacity, calculation accuracy and calculation speed.

[0102] Step 5: Based on the given detector location, obtain discrete grid points with altitude h, longitude, and latitude. Perform two-dimensional six-point interpolation on these two-dimensional discrete grid points to obtain the pseudo-center interpolation results for the corresponding grid points. The bivariate six-point interpolation method is used for any point... The interpolation algorithm is as follows:

[0103]

[0104] Where p and q are the normalized parameters for longitude and latitude:

[0105]

[0106] λ0 is any point The longitude of the corresponding bottom-left discrete grid point. For any point The corresponding latitude, Δλ, and lower left discrete grid point The width of the discrete grid points for longitude and latitude is typically taken as Δλ = 0.5°. c(h) is given in formula (22).

[0107] c 0,0 c 0,1 c 1,0 c 1,1 c 0,-1 and c -1,0 For six discrete points in the latitude and longitude directions at the same altitude as c0, see details. Figure 2 The specific expression is:

[0108]

[0109] Step 6: Calculate the corresponding gravitational acceleration based on the interpolation results, and add J2 and J... 22 The effect of the term ultimately yields the true acceleration, and the specific process is as follows:

[0110] Based on the interpolation result c0, calculate the corresponding pseudo-center distance ρ0:

[0111] ρ0=r-c0 (27)

[0112] Then calculate the corresponding residual acceleration a0:

[0113]

[0114] Add J2 and J 22 The effect of the term ultimately yields gravitational acceleration.

[0115]

[0116] Step 7: Calculate the partial derivative matrix of the acceleration vector with respect to the position vector. The specific process is as follows:

[0117] First, the partial derivative of the pseudo-center c with respect to h is calculated as follows:

[0118]

[0119] Since the fitting function differs for each height interval at different interpolation points, the partial derivative of the corresponding pseudo-center with respect to x can be calculated for each of the six two-dimensional interpolation points.

[0120]

[0121] Where k = x, y, z represents the x, y, z directions of the pseudo-center vector.

[0122] It can be calculated

[0123]

[0124] For the six interpolation points in the latitude and longitude directions, the interpolation values ​​can be obtained separately.

[0125] The formula for calculating the partial derivative of total acceleration is as follows:

[0126]

[0127] in,

[0128]

[0129] In the above formula,

[0130]

[0131]

[0132]

[0133] Calculate using formula (6); Calculate according to formula (34); Calculate according to the formula; Calculate according to formula (13).

[0134] Step 8, based on the calculated partial derivatives of the state transition matrix Substituting into formula (1) and using formula (2) as the initial value, we can solve the initial value problem of the corresponding ordinary differential equation system and obtain the position r(t), velocity v(t), and state transition matrix Φ(t,t0) at any time t. This completes the solution of the orbit prediction problem.

[0135] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A fast prediction method for lunar orbit rendezvous guidance, characterized in that, First, exclude the central gravitational field from the gravitational field acceleration. item, The residual acceleration after gravitational acceleration from the residual gravitational field is converted into an equivalent "pseudo-center" parameter, thus nonlinearly mapping the residual gravitational field acceleration into a "pseudo-center". Then, within a given longitude and latitude point and a specified altitude range, polynomial interpolation is used to calculate the interpolation polynomial coefficients of the "pseudo-center" function with respect to altitude within that altitude range. Next, the lunar surface is divided into several grid points according to longitude and latitude, and the fitting polynomial coefficients of the "pseudo-center" function for each longitude and latitude point within a given altitude range are calculated and saved. Finally, for any spatial location point, the "pseudo-center" corresponding to the instantaneous altitude of the six nearby longitude and latitude grid points is calculated. Then, the "pseudo-center" corresponding to the location point is obtained by two-dimensional interpolation of the calculated "pseudo-centers" of the six longitude and latitude grid points, and the real gravitational field acceleration corresponding to the location point is calculated through nonlinear mapping. Calculate the partial derivative matrix of the acceleration vector with respect to the position vector; based on the actual gravitational field acceleration and the partial derivative matrix of the acceleration vector with respect to the position vector, calculate the position and velocity of the detector; Includes the following steps: Step 1: Given the initial position and velocity of the probe; calculate the acceleration and its partial derivatives for the two-body problem in the lunar central gravitational field; Step 2, compute the corresponding J2 term and J 22 terms perturbation acceleration and partial derivatives; Step 3, calculate the sum of terms J2 and J2 after removing them. 22 The "pseudo-center" of the item; Step 4: Decompose the three-dimensional space into three directions: longitude, latitude, and altitude. Discretize the "pseudo-center" to obtain a polynomial fit of the "pseudo-center" in the altitude direction. Step 5: Obtain the altitude based on the given detector position. Discrete grid points for longitude and latitude; for altitude Two-dimensional discrete grid points of longitude and latitude are interpolated using two-dimensional six-point interpolation to obtain the "pseudo-center" interpolation results of the corresponding grid points; Step 6, according to the interpolation result, the corresponding gravitational acceleration is calculated, and the influence of J2 term and J 22 term is added, and the final real acceleration is obtained; Step 7: Based on the partial derivatives of the calculated state transition matrix, obtain the position, velocity, and state transition matrix at any given time. This completes the solution to the orbit prediction problem. In step 2, The perturbation acceleration of the term is: wherein is the second order harmonic term coefficient, is the Earth equatorial radius; The partial derivative of the acceleration term is: The perturbation acceleration of an item is in, These are the second-order quadratic harmonic coefficients. The geographical longitude of the principal axis of the second-order harmonic term is the azimuth angle of Changzhou on the equatorial ellipse, which is a constant. This is for intermediate calculations. ; The partial derivative of the perturbation acceleration with respect to the position vector is ； In step 3, the position of any point in three-dimensional space under the lunar-centered fixed coordinate system. and gravitational acceleration With pseudo-center The mapping is ,in The actual distance from the lunar center to the probe represents the pseudo distance from the lunar center to the probe. Add pseudo-center ; The Earth's gravitational field constant; To remove and The subsequent perturbation gravitational acceleration satisfies the following equation: in, It is the gravitational acceleration of the probe relative to the lunar-centered fixed coordinate system, calculated using a non-spherical gravitational field acceleration function, including the central gravitational field acceleration and perturbation acceleration, where CBF represents the central celestial body fixed coordinate system; the pseudo-center With perturbation acceleration The relationship is .

2. The prediction method of claim 1, wherein, The acceleration of the central gravitational field two-body problem in step 1 is represented as: in, is the Earth's gravitational field constant, where x , y , z These are the three directional components of the detector position vector. , , ; The partial derivative with respect to the position vector is: 。

Images