A method for analytically predicting parameter periodic variation of a planar weakly stable boundary orbit

Abstract

The application discloses a kind of plane weak stable boundary orbit parameter periodic variation analytical prediction method, belong to aerospace technology field.The application implementation method is: for the orbit plane and two celestial bodies coplanar case, the orbit parameter of spacecraft is obviously changed at apocenter, apocenter orbit parameter is the main factor influencing analytical solution, using apocenter orbit parameter can improve the accuracy of analytical solution.Based on Gauss perturbation equation, the function of eccentricity about mean anomaly variation is derived, the traditional analytical solution solving formula is obtained by integrating eccentricity function in orbit period, apocenter semi-major axis is solved using potential energy method and apocenter eccentricity is solved using iterative method, and the improved analytical solution of eccentricity periodic variation can be obtained by substituting into traditional analytical solution solving formula.According to the formula of pericenter distance, the analytical solution of pericenter distance periodic variation is obtained by substituting into eccentricity analytical solution, and the rapid and accurate prediction of eccentricity and pericenter distance variation of weak stable boundary orbit in deep space exploration mission is realized.

Description

Technical Field

[0001] This invention relates to a method for predicting spacecraft orbital parameters, and more particularly to an analytical method for predicting large eccentricity orbital parameters such as plane weakly stable boundary orbits, belonging to the field of aerospace technology. Background Technology

[0002] Weakly stable boundary orbits are a crucial approach for achieving low-energy transfer and capture in deep space exploration. Due to the sensitivity of weakly stable boundary orbits to parameter variations, studying analytical solutions for periodic changes in orbital parameters can improve the efficiency of orbit determination and reduce the computational demands of deep space probes. Weakly stable boundary orbits are based on the Circular Restricted Three-Body Problem (CRTBP), which lacks a direct analytical solution and often requires orbital integration to derive the variation patterns of orbital parameters. For deep space exploration missions, the vast distance from Earth and the resulting time delay in signal transmission necessitate the spacecraft's ability to autonomously plan and correct its orbit. Therefore, improving the efficiency of solving for parameter variations in weakly stable boundary orbits is a key challenge in deep space exploration missions.

[0003] In the existing analytical solution research, the first technology [1] (Cook G E. Luni-solar perturbations of the orbit of an Earth satellite[J]. Geophysical Journal International, 1962, 6(3): 271-291) proposed the traditional analytical solution method for CRTBP. This method is based on the two-body problem, considers the third-body perturbation of celestial bodies such as the Sun, obtains the Gaussian perturbation equation, and then obtains the differential equations of each orbital parameter as a function of the true anomaly angle, using 1 / Δ 3 The equation is simplified using a first-order Taylor expansion, where Δ refers to the distance between the spacecraft and the sun. Then, by integrating the differential equation over the initial orbital period, an analytical solution for the periodic changes in orbital parameters can be obtained. This method has the advantage of high computational efficiency, but its accuracy is low for orbits with large eccentricities, such as those with weakly stable boundary conditions.

[0004] Unlike the previous technique, Cook is for 1 / Δ 3 Performing a first-order Taylor expansion, in prior art [2] (Paul SN. Orbital perturbations for space situational awareness [D]. Purdue University Graduate School, 2020.), for 1 / Δ 3A second-order Taylor expansion is performed to obtain a second-order analytical solution for the periodic variation of orbital parameters, which contains nine integral correlation coefficients. The advantage of this method is that it can improve the accuracy of long-term prediction analytical solutions. However, the equations contain many integral correlation coefficients, making the solution process complex. Furthermore, this method only studies the analytical solution for eccentricity and is not very meaningful for solving the analytical solution for the periodic variation of orbital parameters in weakly stable boundary conditions. Summary of the Invention

[0005] To address the problems of large errors and low efficiency in existing analytical solutions for the periodic variations of orbital parameters in weakly stable boundary orbits, this invention aims to provide an analytical prediction method for the periodic variations of orbital parameters in a plane with weakly stable boundary orbits. Based on the Gaussian perturbation equation, differential equations for the variation of each orbital parameter with respect to the true anomaly angle are derived. The differential equation for eccentricity is integrated over the initial orbital period to obtain the traditional analytical solution formula. The potential energy method is used to solve for the semi-major axis at the apocentric point, and the iterative method is used to solve for the apocentric point eccentricity. Substituting these into the traditional analytical solution formula proposed by Cook yields an improved analytical solution for the periodic variation of eccentricity. Finally, an improved analytical solution for the periodic variation of pericentricity is obtained using the pericentricity formula. This enables rapid and accurate prediction of the eccentricity and pericentricity variations of weakly stable boundary orbits in deep space exploration missions. This invention has the advantages of high accuracy and high efficiency in predicting the periodic variations of orbital parameters in a plane with weakly stable boundary orbits. The plane with weakly stable boundary orbital parameters include eccentricity and pericentricity.

[0006] The objective of this invention is achieved through the following technical solution:

[0007] This invention discloses an analytical prediction method for the periodic variation of orbital parameters in a weakly stable boundary orbit, considering only the case where the orbital plane is coplanar with two celestial bodies. Since the orbital parameters of a spacecraft change significantly at the apocentric point, the apocentric orbital parameters are the main factor affecting the analytical solution; using the apocentric orbital parameters can improve the accuracy of the analytical solution. Based on the Gaussian perturbation equation, a function of eccentricity with respect to the true anomaly angle is derived. Integrating the eccentricity function over the orbital period yields the traditional analytical solution formula. The semi-major axis at the apocentric point is solved using the potential energy method, and the apocentric eccentricity is solved using an iterative method. Substituting these solutions into the traditional analytical solution formula yields an improved analytical solution for the periodic variation of eccentricity. Based on the pericentric distance calculation formula, substituting these solutions into the analytical solution for eccentricity, an analytical solution for the periodic variation of pericentric distance is obtained, enabling rapid and accurate prediction of the eccentricity and pericentric distance variations of weakly stable boundary orbits in deep space exploration missions.

[0008] This invention discloses an analytical prediction method for the periodic variation of orbital parameters in a weakly stable planar boundary, comprising the following steps:

[0009] Step 1: Construct Gaussian perturbation equations in the second satellite orbital coordinate system. Based on these equations, derive the differential equations for the orbital parameters' variations with respect to the true anomaly angle. These equations include the semi-major axis, eccentricity, right ascension of the ascending node, orbital inclination, argument of perigee, and mean anomaly angle. Considering only the case where the orbital plane is coplanar with both celestial bodies, integrate the semi-major axis function over the orbital period to obtain the analytical solution for the semi-major axis period variation. Similarly, integrate the eccentricity function over the orbital period to obtain the analytical solution for the eccentricity period variation.

[0010] The origin of the second satellite orbital coordinate system is at the satellite's center of mass, and the unit vectors of the three coordinate axes are u... r ,u t ,u n , where u r Along the direction of satellite geocentric distance; u t Perpendicular to u in the instantaneous orbital plane of the satellite r , pointing in the direction of satellite velocity; u n Parallel to the normal to the instantaneous orbital plane.

[0011] The Gaussian perturbation equation is as follows:

[0012]

[0013] Where a, e, Ω, i, ω, M, and f represent the semi-major axis, eccentricity, right ascension of the ascending node, orbital inclination, argument of perigee, mean perigee, and true perigee, respectively; r and n represent the spacecraft's distance from the planet and its orbital angular velocity around the planet, respectively; and F... r F t F n These represent the components of the perturbation force along the three coordinate axes of the second satellite orbit coordinate system.

[0014] The perturbation equation is expressed as a function of the change in the true anomaly angle, and the changes in the orbital parameters before and after the period are solved by integration.

[0015] Will Substituting into equation (1), we get:

[0016]

[0017] The coordinates of the spacecraft are represented as (x, y, z), and the coordinates of the perturbed celestial body are represented as (x, y, z). d ,y d ,z d If ), then the disturbance function is expressed as:

[0018]

[0019] Among them, M d The mass r of the perturbed celestial body dThis represents the distance between the center of mass of the perturbed celestial body and the center of mass of the planet. Δ represents the distance between the spacecraft and the center of mass of the perturbed celestial body, calculated as follows:

[0020] Δ 2 =(xx) d ) 2 +(yy d ) 2 +(zz d ) 2 (4)

[0021] The coordinates of the perturbed celestial body are represented as follows:

[0022]

[0023] In the formula, Ω d i represents the right ascension of the ascending node of the perturbed celestial body. d u represents the orbital inclination of a perturbed celestial body. d This represents the angular difference between the perturbed celestial body and the right ascension of the ascending node.

[0024] The differential of the perturbation function represents the components of the perturbation force along the OX, OY, and OZ directions:

[0025]

[0026] Let φ be the angle between the position vectors of the perturbed celestial body and the spacecraft. The distance between the centers of mass of the spacecraft and the perturbed celestial body can be expressed using the law of cosines:

[0027]

[0028] therefore

[0029]

[0030] The cosine of angle φ is expressed as:

[0031]

[0032] Where u = ω + f,

[0033]

[0034] Perturbation force radial F r tangential to F t Normal direction F n They are represented as follows:

[0035]

[0036] in,

[0037] The semi-major axis varies with the true anterior angle as follows:

[0038]

[0039] Substituting u=ω+f into equation (12) and integrating it over one period, we obtain the change in the semi-major axis Δa over that period. During integration, the motion of the disturbance is considered a small quantity, meaning the position of the disturbance is fixed within the integration period. All integral values ​​shown in equation (13) are zero.

[0040]

[0041] Therefore, we get:

[0042]

[0043] The result obtained after calculation is as follows:

[0044] Δa=0 (15)

[0045] The variation of eccentricity with the true anterior angle is as follows:

[0046]

[0047] Integrating equation (16) yields the periodic change in eccentricity Δe, as shown in equation (17):

[0048]

[0049] Step two: In the two-body system, the spacecraft's mechanical energy is conserved during its orbit around the planet. Therefore, the perturbation force of the sun on the spacecraft is equivalent to an external force, and the work done by this force on the spacecraft is equal to the change in the spacecraft's mechanical energy in the two-body system. By treating the effect of the sun's perturbation force on the spacecraft as work done by the gravitational field, the change in the spacecraft's mechanical energy from its initial position to its apocentric position is calculated. The semi-major axis parameter of the spacecraft at the apocentric point is then obtained using the two-body energy formula.

[0050] For a two-body system, the mechanical energy of the spacecraft is conserved during its orbit around the planet. Therefore, the perturbation force of the solar three-body system can be considered as an external force. The work done by this force on the spacecraft is equal to the change in the spacecraft's mechanical energy within the two-body system. The potential energy method is used to determine the work done by the perturbation force on the spacecraft. The specific process is as follows:

[0051] Analysis of the relative perturbation force F between the solar three-body system and the spacecraft in the two-body system per The relative perturbation force of the three solar bodies in the centroidal inertial coordinate system is obtained from the relative acceleration in theoretical mechanics as shown in equation (18):

[0052] F per =F sat -F earth (18)

[0053] In equation (18), F sat This indicates that the spacecraft is subject to the Sun's gravitational pull in the heliocentric inertial coordinate system, F. earth This indicates that the Earth is subject to the gravitational pull of the Sun in a heliocentric inertial coordinate system.

[0054] F sat and F earth The calculation formula is shown in equation (19):

[0055]

[0056] In equation (19), G represents the gravitational constant, and M... sun Represents the mass of the sun, r d The vector r represents the spacecraft pointing towards the center of the heliosphere. e This represents the vector whose center of gravity points towards the heliocenter.

[0057] The formula for calculating the work done by the solar perturbation on the spacecraft is:

[0058]

[0059] In equation (20), r0 represents the initial position of the spacecraft in the geocentric inertial coordinate system, r ap This indicates the position of the spacecraft at its centroid in the geocentric inertial coordinate system.

[0060] In the following analysis, the position of the Sun is assumed to be constant, corresponding to the same position when the evolution time is a quarter of the orbital period.

[0061] F sat The size is related to the distance between the spacecraft and the center of the heliosphere, therefore F sat The effect is considered as work done by the point source field, which is determined by the distance between the spacecraft and the heliocenter before and after the work is done. The expression (21) is as follows:

[0062]

[0063] In equation (21), s0 and s ap Let be the distances between the spacecraft and the heliocenter at the initial moment and the apocentric moment, respectively. The calculation formula is:

[0064]

[0065] Where, r e The distances from the Sun to the planets, a0 and a0, represent the distances from the Sun to the planets. ap Let e0 and θ0 represent the semi-major axes of the orbit at the initial moment and at the apocentric moment, respectively. Let e0 represent the initial eccentricity of the orbit, and θ0 represent the angle between the spacecraft and the constant position of the sun relative to the center of gravity at the initial moment.

[0066] Since the motion of the sun relative to the earth is equivalent to uniform circular motion, F earth Since the magnitude of the force is constant and the position of the sun is constant, the direction of the force is also definite. Its effect is equivalent to work done by a uniform field, yielding F. earth The expression for the work done on the spacecraft is as follows:

[0067] M earth =F earth (x ap -x0) (23)

[0068] In equation (23), x0 and x ap Let these represent the projected coordinates of the spacecraft's initial position and apocentric position in the direction pointing towards the heliocenter, respectively. The calculation formula is:

[0069]

[0070] The energy of a spacecraft in a two-body system is represented by E:

[0071]

[0072] Therefore, the spacecraft's energy transformation during the entire process is as follows:

[0073] ΔE=M=M sat -M earth (26)

[0074] Based on the energy change and the relationship between the semi-major axis and energy (25), the change in the semi-major axis can be obtained, and then the semi-major axis a at the centroid can be calculated. ap .

[0075] Step 3: When the spacecraft reaches the apocentric point and the orbital parameters change significantly, the analytical solution for the periodic variation of eccentricity is improved using the apocentric point parameters. The semi-major axis of the apocentric point obtained in Step 2 is substituted into the analytical solution formula for the periodic variation of eccentricity to calculate the periodic variation of eccentricity. The eccentricity at the apocentric point is approximately the average of the eccentricities at the initial and final moments. The apocentric eccentricity is substituted into the analytical solution formula for the periodic variation of eccentricity for iterative calculation until the preset convergence condition is met, thus obtaining the analytical solution for the periodic variation of eccentricity and achieving analytical prediction of the periodic variation of eccentricity.

[0076] The main changes in orbital parameters occur when the true anomaly angle is π, i.e., when the spacecraft reaches the apocenter. At this moment, the semi-major axis and eccentricity of the large eccentricity orbit have changed significantly compared to the initial point. Therefore, the semi-major axis 'a' of the orbital parameters at the apocenter is adopted. ap and eccentricity e ap Improve the analytical solution.

[0077] The orbital eccentricity changes monotonically within the period, and the eccentricity value corresponding to the centroid is close to the average eccentricity corresponding to the initial and final states. Therefore, the average eccentricity corresponding to the preset initial state and the final state of the period is equivalent to the eccentricity at the centroid.

[0078] According to equation (17), the periodic change in eccentricity can be expressed as a function that depends only on the eccentricity at the distal point:

[0079]

[0080] First time solving for the eccentricity e at the centroid ap Substituting the initial eccentricity as the centroidal eccentricity into equation (27) yields the periodic change in eccentricity (Δe)1 and the eccentricity e1 at the end of the period:

[0081] (Δe)1=f(e0) (28)

[0082] e1=e0+(Δe)1 (29)

[0083] Using eccentricity e1 as the eccentricity at the end of the period, and combining it with the initial eccentricity e0, the eccentricity at the centroid is calculated:

[0084]

[0085] Substituting into equation (27) again, we get:

[0086] (Δe)2=f(e ap (31)

[0087] Repeat the above process until the preset convergence condition is met. The result is the eccentricity of the state at the end of the cycle, which means the prediction of the eccentricity of the state at the end of the cycle is achieved.

[0088] Step four: The periodic change of the semi-major axis is zero, meaning the semi-major axis at the end of the period is equal to that at the beginning. Substitute the periodic change of the eccentricity obtained in step three into the formula for calculating the proximal distance to obtain an analytical solution for the periodic change of the proximal distance.

[0089] Proximal distance r p It can be directly obtained from the eccentricity and semi-major axis:

[0090] r p =a(1-e) (32)

[0091] Considering that the periodic variation of the semi-major axis is 0, therefore

[0092] Δr p =-a0Δe (33)

[0093] The improved analytical solution of the perturbation equation of the eccentricity is obtained from step three. At the same time, the semi-major axis state quantity at the end of the period is equal to the semi-major axis state quantity at the initial point. Then, substituting a0 and Δe into equation (33) yields the analytical solution of the periodic change of the pericentric distance, thus realizing the prediction of the eccentricity at the end of the period.

[0094] Step 5: Based on the semi-major axis of the centroid obtained in Step 2, the analytical solution of the periodic variation of the eccentricity is obtained by iterative method in Step 3, and the analytical solution of the periodic variation of the pericentric distance is obtained by the pericentric distance formula in Step 4, thereby realizing the analytical prediction of the periodic variation of the orbital parameters of the plane weakly stable boundary.

[0095] Beneficial effects:

[0096] 1. This invention discloses an analytical prediction method for the periodic variation of orbital parameters in a plane with weak stability. It improves the traditional analytical solution by using the orbital parameters at the centroid point. The semi-major axis and eccentricity of the centroid point are substituted for the semi-major axis and eccentricity of the initial point and used in the traditional analytical solution formula to obtain the analytical solution for the periodic variation of eccentricity. The analytical solution for the periodic variation of pericentricity is obtained by using the pericentricity formula, thereby achieving efficient and accurate prediction of the periodic variation of orbital parameters in a plane with weak stability.

[0097] 2. The present invention discloses an improved analytical solution for the periodic variation of planar weakly stable boundary orbit parameters. Based on the conservation of mechanical energy during the spacecraft's orbit around the planet, the perturbation effect of the sun on the spacecraft is regarded as the work done by the gravitational field. The change in mechanical energy of the spacecraft from the initial position to the apocentric position is obtained, and the semi-major axis parameters of the spacecraft at the apocentric point are obtained according to the two-body energy formula.

[0098] 3. The present invention discloses an improved analytical solution for the periodic variation of the orbital parameters of a plane with weak stability. Based on the fact that the eccentricity period is monotonically changing and the parameter value at the centroid is approximately equal to the average value of the parameters at the beginning and end, the eccentricity parameter value at the centroid is calculated using the obtained periodic variation of the eccentricity, and then substituted into the formula for calculating the periodic variation of the eccentricity for iterative solution until the preset convergence condition is met. Attached Figure Description

[0099] Figure 1 This is a flowchart of an analytical prediction method for periodic changes in planar weakly stable boundary orbit parameters according to the present invention.

[0100] Figure 2 This is a schematic diagram of the second satellite orbit coordinate system in which the Gaussian perturbation equation of step one of the present invention is located;

[0101] Figure 3 This is a schematic diagram illustrating the variation law of orbital parameters within a period in step two of this invention, wherein... Figure 3 (a) is a graph showing the rate of change of the semi-major axis as a function of the true anterior angle. Figure 3(b) is a graph showing the rate of change of eccentricity as a function of the true anterior angle. Figure 3 (c) is a graph showing the rate of change of the right ascension of the ascending node as a function of the true anomaly angle. Figure 3 (d) is a graph showing the rate of change of orbital inclination as a function of the true anomaly angle;

[0102] Figure 4 This is a schematic diagram illustrating the solution of the semi-major axis error at the telecentric point using the potential energy method in step three of this invention. Figure 4 (a) is a schematic diagram when the perigee argument is 0°. Figure 4 (b) is a schematic diagram when the perigee argument is 30°. Figure 4 (c) is a schematic diagram when the perigee argument is 60°. Figure 4 (d) is a schematic diagram when the perigee argument is 90°. Figure 4 (e) is a schematic diagram when the perigee argument is 120°. Figure 4 (f) is a schematic diagram when the perigee argument is 150°;

[0103] Figure 5 This is a schematic diagram illustrating the variation law of eccentricity within the period in step four of this invention, wherein... Figure 5 (a) is a schematic diagram when the perigee argument is 0°. Figure 5 (b) is a schematic diagram when the perigee argument is 30°. Figure 5 (c) is a schematic diagram when the perigee argument is 60°. Figure 5 (d) is a schematic diagram when the perigee argument is 90°. Figure 5 (e) is a schematic diagram when the perigee argument is 120°. Figure 5 (f) is a schematic diagram when the perigee argument is 150°;

[0104] Figure 6 This is a schematic diagram showing the results of calculating the periodic change in eccentricity using different methods in step four of this invention. Detailed Implementation

[0105] To better illustrate the purpose and advantages of the present invention, the invention will be further described below in conjunction with the accompanying drawings and examples.

[0106] like Figure 1 As shown in the figure, the analytical prediction method for periodic variations of plane weakly stable boundary orbit parameters disclosed in this embodiment includes the following steps:

[0107] Taking an orbit whose initial point is located in the ecliptic plane and whose orbital plane coincides with the ecliptic plane as an example, assuming that the initial point pericentric distance is 8000km, the initial perigee argument is 45°, the initial eccentricity is 0.98, and the initial semi-major axis is 400000km, the changes in orbital parameters such as eccentricity and pericentric distance are calculated using the improved solution method proposed in this patent.

[0108] Step 1: Construct Gaussian perturbation equations in the second satellite orbital coordinate system. Based on these equations, derive the differential equations for the orbital parameters' variations with respect to the true anomaly angle. These equations include the semi-major axis, eccentricity, right ascension of the ascending node, orbital inclination, argument of perigee, and mean anomaly angle. Considering only the case where the orbital plane is coplanar with both celestial bodies, integrate the semi-major axis function over the orbital period to obtain the analytical solution for the semi-major axis period variation. Similarly, integrate the eccentricity function over the orbital period to obtain the analytical solution for the eccentricity period variation.

[0109] The origin of the second satellite orbital coordinate system is at the satellite's center of mass, and the unit vectors of the three coordinate axes are u... r ,u t ,u n , where u r Along the direction of satellite geocentric distance; u t Perpendicular to u in the instantaneous orbital plane of the satellite r , pointing in the direction of satellite velocity; u n Parallel to the normal to the instantaneous orbital plane.

[0110] The Gaussian perturbation equation is as follows:

[0111]

[0112] Where a, e, Ω, i, ω, M, and f represent the semi-major axis, eccentricity, right ascension of the ascending node, orbital inclination, argument of perigee, mean perigee, and true perigee, respectively; r and n represent the spacecraft's distance from the planet and its orbital angular velocity around the planet, respectively; and F... r F t F n These represent the components of the perturbation force along the three coordinate axes of the second satellite orbit coordinate system.

[0113] The perturbation equation is expressed as a function of the change in the true anomaly angle, and the changes in the orbital parameters before and after the period are solved by integration.

[0114] Will Substituting into equation (34), we get:

[0115]

[0116] The coordinates of the spacecraft are represented as (x, y, z), and the coordinates of the perturbed celestial body are represented as (x, y, z). d ,y d ,z d If ), then the disturbance function is expressed as:

[0117]

[0118] Among them, M d The mass r of the perturbed celestial bodyd This represents the distance between the center of mass of the perturbed celestial body and the center of mass of the planet. Δ represents the distance between the spacecraft and the center of mass of the perturbed celestial body, calculated as follows:

[0119] Δ 2 =(xx) d ) 2 +(yy d ) 2 +(zz d ) 2 (37)

[0120] The coordinates of the perturbed celestial body are represented as follows:

[0121]

[0122] In the formula, Ω d i represents the right ascension of the ascending node of the perturbed celestial body. d u represents the orbital inclination of a perturbed celestial body. d This represents the angular difference between the perturbed celestial body and the right ascension of the ascending node.

[0123] The differential of the perturbation function represents the components of the perturbation force along the OX, OY, and OZ directions:

[0124]

[0125] Let φ be the angle between the position vectors of the perturbed celestial body and the spacecraft. The distance between the centers of mass of the spacecraft and the perturbed celestial body can be expressed using the law of cosines:

[0126]

[0127] therefore

[0128]

[0129] The cosine of angle φ is expressed as:

[0130]

[0131] Where u = ω + f,

[0132]

[0133] Perturbation force radial F r tangential to F t Normal direction F n They are represented as follows:

[0134]

[0135] in,

[0136] The semi-major axis varies with the true anterior angle as follows:

[0137]

[0138] Substituting u=ω+f into equation (45) and integrating it over one period, we obtain the change in the semi-major axis Δa over that period. During integration, the motion of the disturbance is considered a small quantity, i.e., the position of the disturbance is fixed within the integration period. All integral values ​​shown in equation (46) are zero.

[0139]

[0140] Therefore, we get:

[0141]

[0142] The result obtained after calculation is as follows:

[0143] Δa=0 (48)

[0144] The variation of eccentricity with the true anterior angle is as follows:

[0145]

[0146] Integrating equation (49) yields the periodic change in eccentricity Δe, as shown in equation (50):

[0147]

[0148] Step two: In the two-body system, the spacecraft's mechanical energy is conserved during its orbit around the planet. Therefore, the perturbation force of the sun on the spacecraft is equivalent to an external force, and the work done by this force on the spacecraft is equal to the change in the spacecraft's mechanical energy in the two-body system. By treating the effect of the sun's perturbation force on the spacecraft as work done by the gravitational field, the change in the spacecraft's mechanical energy from its initial position to its apocentric position is calculated. The semi-major axis parameter of the spacecraft at the apocentric point is then obtained using the two-body energy formula.

[0149] For a two-body system, the mechanical energy of the spacecraft is conserved during its orbit around the planet. Therefore, the perturbation force of the solar three-body system can be considered as an external force. The work done by this force on the spacecraft is equal to the change in the spacecraft's mechanical energy within the two-body system. The potential energy method is used to determine the work done by the perturbation force on the spacecraft. The specific process is as follows:

[0150] Analysis of the relative perturbation force F between the solar three-body system and the spacecraft in the two-body system per The relative perturbation force of the three solar bodies in the centroidal inertial coordinate system is obtained from the relative acceleration in theoretical mechanics as shown in equation (51):

[0151] F per =F sat -F earth(51)

[0152] In equation (51), F sat This indicates that the spacecraft is subject to the Sun's gravitational pull in the heliocentric inertial coordinate system, F. earth This indicates that the Earth is subject to the gravitational pull of the Sun in a heliocentric inertial coordinate system.

[0153] F sat and F earth The calculation formula is shown in equation (52):

[0154]

[0155] In equation (52), G represents the gravitational constant, and M... sun Represents the mass of the sun, r d The vector r represents the spacecraft pointing towards the center of the heliosphere. e This represents the vector whose center of gravity points towards the heliocenter.

[0156] The formula for calculating the work done by the solar perturbation on the spacecraft is:

[0157]

[0158] In equation (53), r0 represents the initial position of the spacecraft in the geocentric inertial coordinate system, r ap This indicates the position of the spacecraft at its centroid in the geocentric inertial coordinate system.

[0159] In the following analysis, the position of the Sun is assumed to be constant, corresponding to the same position when the evolution time is a quarter of the orbital period.

[0160] F sat The size is related to the distance between the spacecraft and the center of the heliosphere, therefore F sat The effect is considered as work done by the point source field, which is determined by the distance between the spacecraft and the heliocenter before and after the work is done. The expression (54) is as follows:

[0161]

[0162] In equation (54), s0 and s ap Let be the distances between the spacecraft and the heliocenter at the initial moment and the apocentric moment, respectively. The calculation formula is:

[0163]

[0164] Where, r e The distances from the Sun to the planets, a0 and a0, represent the distances from the Sun to the planets. ap Let e0 and θ0 represent the semi-major axes of the orbit at the initial moment and at the apocentric moment, respectively. Let e0 represent the initial eccentricity of the orbit, and θ0 represent the angle between the spacecraft and the constant position of the sun relative to the center of gravity at the initial moment.

[0165] Since the motion of the sun relative to the earth is equivalent to uniform circular motion, F earth Since the magnitude of the force is constant and the position of the sun is constant, the direction of the force is also definite. Its effect is equivalent to work done by a uniform field, yielding F. earth The expression for the work done on the spacecraft is as follows:

[0166] M earth =F earth (x ap -x0) (56)

[0167] In equation (56), x0 and x ap Let these represent the projected coordinates of the spacecraft's initial position and apocentric position in the direction pointing towards the heliocenter, respectively. The calculation formula is:

[0168]

[0169] The energy of a spacecraft in a two-body system is represented by E:

[0170]

[0171] Therefore, the spacecraft's energy transformation during the entire process is as follows:

[0172] ΔE=M=M sat -M earth (59)

[0173] Based on the energy change and the relationship between the semi-major axis and energy (59), the change in the semi-major axis can be obtained, and then the semi-major axis a at the centroid can be calculated. ap =401230km.

[0174] Step 3: When the spacecraft reaches the apocentric point and the orbital parameters change significantly, the analytical solution for the periodic variation of eccentricity is improved using the apocentric point parameters. The semi-major axis of the apocentric point obtained in Step 2 is substituted into the analytical solution formula for the periodic variation of eccentricity to calculate the periodic variation of eccentricity. The eccentricity at the apocentric point is approximately the average of the eccentricities at the initial and final moments. The apocentric eccentricity is substituted into the analytical solution formula for the periodic variation of eccentricity for iterative calculation until the preset convergence condition is met, thus obtaining the analytical solution for the periodic variation of eccentricity and achieving analytical prediction of the periodic variation of eccentricity.

[0175] The main changes in orbital parameters occur when the true anomaly angle is π, i.e., when the spacecraft reaches the apocenter. At this moment, the semi-major axis and eccentricity of the large eccentricity orbit have changed significantly compared to the initial point. Therefore, the semi-major axis 'a' of the orbital parameters at the apocenter is adopted. ap and eccentricity e ap Improve the analytical solution.

[0176] The orbital eccentricity changes monotonically within the period, and the eccentricity value corresponding to the centroid is close to the average eccentricity corresponding to the initial and final states. Therefore, the average eccentricity corresponding to the preset initial state and the final state of the period is equivalent to the eccentricity at the centroid.

[0177] According to equation (17), the periodic change in eccentricity can be expressed as a function that depends only on the eccentricity at the distal point:

[0178]

[0179] First time solving for the eccentricity e at the centroid ap Substituting the initial eccentricity as the centroidal eccentricity into equation (60) yields the periodic change in eccentricity (Δe)1 and the eccentricity e1 at the end of the period:

[0180] (Δe)1=f(e0)=-0.0254 (61)

[0181] e1=e0+(Δe)1=0.9546 (62)

[0182] Using eccentricity e1 as the eccentricity at the end of the period, and combining it with the initial eccentricity e0, the eccentricity at the centroid is calculated:

[0183]

[0184] Substituting back into equation (60), we get:

[0185] (Δe)2=f(e ap (64)

[0186] Repeating the above process seven times satisfies the convergence condition |(Δe) n -(Δe) n-1 If | < 0.00001, the result is the eccentricity at the end of the period, e = 0.9631, and the change in eccentricity Δe = -0.0169.

[0187] Step four: The periodic change of the semi-major axis is zero, meaning the semi-major axis at the end of the period is equal to that at the beginning. Substitute the periodic change of the eccentricity obtained in step three into the formula for calculating the proximal distance to obtain an analytical solution for the periodic change of the proximal distance.

[0188] Proximal distance r p It can be directly obtained from the eccentricity and semi-major axis:

[0189] r p =a(1-e) (65)

[0190] Considering that the periodic variation of the semi-major axis is 0, therefore

[0191] Δr p=-a0Δe=13520km (66)

[0192] The improved analytical solution of the perturbation equation of the eccentricity is obtained from step three. At the same time, the semi-major axis state quantity at the end of the period is equal to the semi-major axis state quantity at the initial point. Then, substituting a0 and Δe into equation (66) yields the analytical solution of the periodic change of the pericentric distance, thus realizing the prediction of the eccentricity at the end of the period.

[0193] Step 5: Based on the semi-major axis of the centroid obtained in Step 2, the analytical solution of the periodic variation of the eccentricity is obtained by iterative method in Step 3, and the analytical solution of the periodic variation of the pericentric distance is obtained by the pericentric distance formula in Step 4, thereby realizing the analytical prediction of the periodic variation of the orbital parameters of the plane weakly stable boundary.

[0194] Step Six: Based on the analytical prediction results of the periodic changes of the plane weakly stable boundary orbit parameters obtained in Step Five, spacecraft conducting deep space exploration missions can provide autonomous trajectory planning and trajectory correction in deep space, thereby improving the efficiency and accuracy of deep space exploration missions.

[0195] The above detailed description further illustrates the purpose, technical solution, and beneficial effects of the invention. It should be understood that the above description is only a specific embodiment of the present invention and is 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. An analytical prediction method for periodic variations of orbital parameters in a weakly stable plane boundary, characterized in that: Includes the following steps, Step 1: Construct Gaussian perturbation equations in the second satellite orbital coordinate system. Based on these equations, derive differential equations for the orbital parameters with respect to the true anomaly angle, including the semi-major axis, eccentricity, right ascension of the ascending node, orbital inclination, argument of perigee, and mean anomaly angle. Considering only the case where the orbital plane is coplanar with the two celestial bodies, integrate the semi-major axis function over the orbital period to obtain the analytical solution for the semi-major axis period variation. Integrate the eccentricity function over the orbital period to obtain the analytical solution for the eccentricity period variation. Step 2: In the two-body system, the mechanical energy of the spacecraft is conserved during its orbit around the planet. Therefore, the perturbation force of the sun on the spacecraft is equivalent to an external force. The work done by this force on the spacecraft is equal to the change in the mechanical energy of the spacecraft in the two-body system. By treating the effect of the sun's perturbation force on the spacecraft as work done by the gravitational field, the change in the mechanical energy of the spacecraft from the initial position to the apocentric position is obtained. The semi-major axis parameter of the spacecraft at the apocentric point is obtained according to the two-body energy formula. In step two, For a two-body system, the mechanical energy of the spacecraft is conserved during its orbit around the planet. Therefore, the perturbation force of the solar three-body system can be considered as an external force. The work done by this force on the spacecraft is equal to the change in the spacecraft's mechanical energy in the two-body system. The potential energy method is used to solve for the work done by the perturbation force on the spacecraft. The specific implementation method is as follows: Analysis of the relative perturbation forces exerted on the spacecraft by the solar three-body system in the two-body system The relative perturbation force of the three solar bodies in the centroidal inertial coordinate system is obtained from the relative acceleration in theoretical mechanics as shown in equation (18): Mode middle, This indicates that the spacecraft is subject to the Sun's gravity in a heliocentric inertial coordinate system. This indicates that the Earth experiences the gravitational pull of the Sun in a heliocentric inertial coordinate system. and The calculation formula is shown in equation (19): Mode middle, Represents the gravitational constant. Indicates the mass of the sun. The vector representing the spacecraft pointing towards the center of the heliosphere. The vector representing the direction from the center of the orbit to the center of the heliosphere; Indicates the distance the spacecraft points towards the center of the sun. Indicates the distance between the center of the orbit and the center of the sun; The formula for calculating the work done by the solar perturbation on the spacecraft is: In equation (20), This indicates the initial position of the spacecraft in the geocentric inertial coordinate system. This indicates the position of the spacecraft at its centroid in the geocentric inertial coordinate system. In the following analysis, the position of the Sun is assumed to be constant, corresponding to the same position when the evolution time is a quarter of the orbital period; The size is related to the distance between the spacecraft and the center of the heliosphere, therefore... The effect is considered as work done by the point source field, which is determined by the distance between the spacecraft and the heliocenter before and after the work is done. The expression (21) is as follows: Mode middle, and Let be the distances between the spacecraft and the heliocenter at the initial moment and the apocentric moment, respectively. The calculation formula is: in, It represents the distance from the sun to a planet. and Let these represent the semi-major axes of the orbit at the initial time and the time at the centroid, respectively. This represents the initial eccentricity of the orbit. This represents the angle between the spacecraft and the sun's constant position relative to the center of gravity at the initial moment; Since the Sun's motion relative to the Earth is equivalent to uniform circular motion... The magnitude of the force is constant, and the position of the sun is constant, therefore the direction of the force is also definite. Its effect is equivalent to work done by a uniform field, resulting in... The expression for the work done on the spacecraft is as follows: Mode middle, and Let these represent the projected coordinates of the spacecraft's initial position and apocentric position in the direction pointing towards the heliocenter, respectively. The calculation formula is: The spacecraft's energy in the two-body system is used To indicate: in, Let represent the semi-major axis of the orbit. Therefore, the spacecraft's energy transformation throughout the process is as follows: The change in the semi-major axis can be obtained by considering the energy change and the relationship between the semi-major axis and energy (25), and then the semi-major axis at the centroid can be determined. ; Step 3: When the spacecraft reaches the apocentric point and the orbital parameters change significantly, the analytical solution for the periodic change of eccentricity is improved using the apocentric point parameters. The semi-major axis of the apocentric point obtained in Step 2 is substituted into the analytical solution formula for the periodic change of eccentricity to calculate the periodic change of eccentricity. The eccentricity of the apocentric point is approximately the average value of the eccentricity at the beginning and end. The eccentricity of the apocentric point is substituted into the analytical solution formula for the periodic change of eccentricity for iterative calculation until the preset convergence condition is met to obtain the analytical solution for the periodic change of eccentricity, thereby realizing the analytical prediction of the periodic change of eccentricity. Step 4: The semi-major axis periodic change is zero, that is, the semi-major axis at the end of the period is equal to that at the beginning of the period; substitute the eccentricity periodic change obtained in Step 3 into the proximal distance calculation formula to solve the analytical solution of the proximal distance periodic change. Step 5: Based on the semi-major axis of the centroid obtained in Step 2, the analytical solution of the periodic variation of the eccentricity is obtained by iterative method in Step 3, and the analytical solution of the periodic variation of the pericentric distance is obtained by the pericentric distance formula in Step 4, thereby realizing the analytical prediction of the periodic variation of the orbital parameters of the plane weakly stable boundary.

2. The analytical prediction method for periodic variations of plane weakly stable boundary orbit parameters as described in claim 1, characterized in that: The implementation method for step one is as follows: The origin of the second satellite orbital coordinate system is at the satellite's center of mass, and the unit vectors of the three coordinate axes are respectively... ,in Along the direction of satellite geocentric distance; Perpendicular to the satellite's instantaneous orbital plane , pointing in the direction of the satellite's speed; Parallel to the normal to the instantaneous orbital plane; The Gaussian perturbation equation is as follows: in, Let these represent the semi-major axis, eccentricity, right ascension of the ascending node, inclination, argument of perigee, mean anomaly, and true anomaly, respectively. These represent the spacecraft's distance from the planet and its angular velocity around the planet, respectively. These represent the components of the perturbation force along the three coordinate axes of the second satellite orbit coordinate system, respectively; The perturbation equation is expressed as a function of the true anomaly angle, and the changes in orbital parameters before and after the period are solved by integration. Will Substitution get: The coordinates of the spacecraft are represented as The coordinates of the perturbed celestial body are represented as The disturbance function is then expressed as: in, This indicates the mass of the disturbed celestial body. This represents the distance between the center of mass of the perturbed celestial body and the center of mass of the planet; The distance between the spacecraft and the center of mass of the perturbed celestial body is expressed by the following formula: The coordinates of the perturbed celestial body are represented as follows: In the formula, The right ascension of the ascending node of the perturbed celestial body. Indicates the orbital inclination of the perturbed celestial body. This represents the angular difference between the perturbed celestial body and the right ascension of the ascending node; The differential representation of the perturbation function along the perturbation force , and Components of direction: The angle between the position vectors of the perturbed celestial body and the spacecraft is The distance between the center of mass of the spacecraft and the perturbed celestial body is expressed using the law of cosines: therefore angle The cosine of is expressed as: in , Perturbation radial tangential , legal direction They are represented as follows: in, ; The semi-major axis varies with the true anterior angle as follows: Will Substituting into equation (12) and integrating it over one period, we obtain the change in the semi-major axis over that period. During the integration process, the motion of the disturbance is considered as a small quantity, that is, the position of the disturbance is fixed during the integration period; all integral values ​​as shown in equation (13) are zero: Therefore, we get: The result obtained after calculation is as follows: The eccentricity varies with the true anterior angle as follows: Pair Integrating the eccentricity periodic change yields the result. As shown in equation (17): 。 3. The analytical prediction method for periodic variations of plane weakly stable boundary orbit parameters as described in claim 2, characterized in that: The method for implementing step three is as follows: The main changes in orbital parameters occur at the true anomaly angle of 0. This refers to the point at which the spacecraft reaches its apocenter. At this moment, the semi-major axis and eccentricity of the orbit with large eccentricity have changed significantly compared to the initial point. Therefore, the semi-major axis of the orbital parameters at the apocenter is used. and eccentricity Improve the analytical solution; The changes in orbital eccentricity within the period are all monotonic, and the eccentricity value corresponding to the centroid point is close to the average value of the eccentricity corresponding to the initial and final states. Therefore, the average value of the eccentricity corresponding to the preset initial state and the final state of the period is equivalent to the eccentricity at the centroid point. According to the formula The periodic variation of eccentricity is expressed as a function that depends only on the eccentricity at the distal point: First time solving for the eccentricity of the centroid The initial eccentricity is then substituted into equation (27) as the centroid eccentricity to calculate the periodic change in eccentricity. and the eccentricity at the end of the period : Eccentricity As the eccentricity at the end of the period, combined with the eccentricity at the beginning of the period. Find the eccentricity of the centroid: Substituting back into equation (27), we get: Repeat the above process until the preset convergence condition is met. The result is the eccentricity of the state at the end of the cycle, which means the prediction of the eccentricity of the state at the end of the cycle is achieved.

4. The analytical prediction method for periodic variations of plane weakly stable boundary orbit parameters as described in claim 3, characterized in that: Step four is implemented as follows: Proximal distance It can be directly obtained from the eccentricity and semi-major axis: Considering that the periodic variation of the semi-major axis is 0, therefore The improved analytical solution of the eccentricity perturbation equation is obtained from step three. Simultaneously, the semi-major axis state variable at the end of the period is equal to the semi-major axis state variable at the initial point. Therefore, and Substituting into equation (33) yields the analytical solution for the periodic change of the proximal distance, thus enabling the prediction of the eccentricity rate at the end of the period.

Images