Characterization parameter analysis method of earth-moon collinear translation point orbit

Abstract

The invention relates to a characterization parameter analysis method for an earth-moon collinear translation point orbit. The method comprises the following steps: describing a circular restrictive three-body problem as a chaotic Hamiltonian power system, converting a motion equation into a new coordinate system for series expansion, and taking a quadratic term of a Hamiltonian function in a polynomial form as a linear model of a translation point in the circular restrictive three-body problem; carrying out complex transformation and regular transformation on the linearized model transformation to realize decoupling of a hyperbolic unstable direction and a central direction of the central manifold, and constructing a mapping relation between CRTBP coordinates and characterization parameters to carry out parameter characterization; analyzing the two motion characteristics in the hyperbolic unstable direction and the central direction and the corresponding characterization parameters, determining the orbit of the target spacecraft through the phase difference of the angular variables in the characterization parameters, and obtaining the orbit transfer information of the target spacecraft through the change condition of the hyperbolic motion components in the characterization parameters. By adopting the method, characterization parameter analysis can be realized.

Description

Technical Field

[0001] This application relates to the technical field of spatial characterization parameter analysis, and particularly to a method for analyzing characterization parameters of the Earth-Moon collinear libration point orbit. Background Art

[0002] In recent years, with the exploration activities of the moon carried out by various countries, the number of Earth-Moon space targets has been increasing year by year. Due to its special position and dynamic characteristics, the Earth-Moon libration points are important resources for future Earth-Moon space development. The ARTEMIS mission was launched in 2007, deploying lunar exploration spacecraft in the Lissajous orbits at the Earth-Moon L1 and L2 points; CE5-T1 was launched in 2014 and entered the Lissajous orbit at the Earth-Moon L2 point to perform exploration tasks; Queqiao entered the Halo orbit at the Earth-Moon L2 point in 2018 to perform communication relay tasks; the CAPSTONE mission verified the autonomous navigation technology of the Lunar Reconnaissance Orbiter in the NRHO orbit, reducing the dependence on ground control; in addition, NASA's GATEWAY mission was launched in 2024, and various robotic and human lunar landing missions will be gradually realized. The development of the Earth-Moon space will lead to the future libration point area becoming crowded. In order to improve the Earth-Moon space situation awareness ability, it is necessary to explore an efficient and intuitive target cataloging method for the Earth-Moon space libration point area.

[0003] Space target cataloging refers to the process of systematically identifying, tracking, recording, and classifying space objects in orbit. These space objects include satellites, space debris, abandoned spacecraft, and other natural and artificial celestial bodies, etc. The purpose of space target cataloging is to maintain a detailed database recording information such as the orbital parameters, size, shape, and use of these objects, in order to facilitate space situation awareness, collision warning, space traffic management, etc. Existing space target cataloging methods mainly form a recognized standard format based on Keplerian orbital elements: Two-Line Orbital Element (TLE). In addition, there are also methods using regular Delaunay variables to describe orbital motion. Considering the perturbed two-body problem, the Keplerian orbital elements are no longer integrals of motion, and at this time, the method of circular orbital elements can also be considered.

[0004] Due to the strong gravitational forces from the Earth and the Moon in the Earth-Moon space, mission analysts cannot easily rely on typical two-body Keplerian orbital elements to determine the orbit of a spacecraft. To better describe the motion of an object at the Earth-Moon libration points, the Circular Restricted Three-Body Problem (CRTBP) model is often adopted, but this causes great difficulties in cataloging the objects at the Earth-Moon libration points. On the one hand, the two-body problem is integrable, while the restricted three-body problem is non-integrable, making it impossible to extract more first integrals in the restricted three-body problem in Cartesian coordinates except for energy conservation. Therefore, theoretically, it does not have a cataloging basis similar to the two-body problem. On the other hand, there is no generally recognized parametric description method for the special dynamic structures near the libration points, such as periodic orbits, quasi-periodic orbits, hyperbolic invariant manifolds, etc. It is difficult to directly correlate the physical characteristics, such as amplitude and period, with the CRTBP dynamics.

[0005] Therefore, the core of the problem of cataloging the Earth-Moon libration points lies in how to transform the CRTBP into an integrable system and extract the corresponding characteristic parameters, and it is necessary to analyze what parameters can be used as characteristic parameters. Summary of the Invention

[0006] Based on this, in view of the above technical problems, it is necessary to provide a method for analyzing the characteristic parameters of the Earth-Moon collinear libration point orbit that can achieve the analysis of characteristic parameters.

[0007] A method for analyzing the characteristic parameters of the Earth-Moon collinear libration point orbit, the method includes:

[0008] Obtain the motion equation of the spacecraft in the rendezvous coordinate system under the assumption of the Circular Restricted Three-Body Problem, describe the Circular Restricted Three-Body Problem as a chaotic Hamiltonian dynamical system, transform the motion equation to a new coordinate system, and obtain the Hamiltonian function of the corresponding system after the coordinate transformation;

[0009] Use the Legendre expansion method to perform a series expansion on the non-linear terms in the Hamiltonian function to obtain a Hamiltonian function in polynomial form; use the quadratic term of the Hamiltonian function in polynomial form as the linearized model of the libration point in the Circular Restricted Three-Body Problem;

[0010] According to the real linear symplectic transformation matrix, transform the linearized model into a real canonical form and perform a complex transformation on the central part of the real canonical form to obtain a linear complex canonical form; perform a canonical transformation on the non-linear terms of the linear complex canonical form to decouple the hyperbolic unstable direction and the central direction of the central manifold, and obtain the final Hamiltonian function; refer to the linearly integrable part of the final Hamiltonian function, define local action-angle variables according to the motion mode of the libration point to describe the motion of the spacecraft on the central manifold, and construct a mapping relationship between the CRTBP coordinates and the characteristic parameters for parameter characterization;

[0011] Analyze the two motion characteristics of the hyperbolic unstable direction and the central direction and the corresponding characterization parameters. The orbit of the target spacecraft can be determined by the phase difference of the angular variables in the characterization parameters, and the orbit transfer information of the target spacecraft can be obtained by the change of the hyperbolic motion component in the characterization parameters; the orbit transfer information includes whether the target spacecraft enters / leaves the libration point periodic orbit and which invariant manifold it enters.

[0012] The above method for analyzing the characterization parameters of the Earth-Moon collinear libration point orbit first obtains the motion equation of the spacecraft in the rendezvous coordinate system under the assumption of the circular restricted three-body problem, and describes it as a chaotic Hamiltonian dynamical system. The motion equation is transformed to a new coordinate system to obtain the corresponding Hamiltonian function, which is convenient for subsequent transformations and analyses. In this way, the motion of the spacecraft is related to the Hamiltonian function, enabling the study of the motion characteristics of the spacecraft from the perspective of the Hamiltonian function and laying a foundation for extracting the characterization parameters. The Legendre expansion method is used to expand the nonlinear terms in the Hamiltonian function into a series, transforming the complex nonlinear relationship into a polynomial form for easy analysis and processing. The quadratic term of the Hamiltonian function in polynomial form serves as the linearized model of the libration point in the circular restricted three-body problem. Near the libration point, the quadratic term can well approximate the local behavior of the system. The linearized model is relatively simple and easy to analyze. According to the real linear symplectic transformation matrix, the linearized model is transformed into a real canonical form, and a complex transformation is performed on the central part of the real canonical form to obtain a linear complex canonical form, transforming the equation of the system into a form more convenient for analysis. The real canonical form and the linear complex canonical form have specific structures and properties, which can clearly separate the different motion characteristics of the system, enabling targeted analysis of the motion characteristics of different parts. A canonical transformation is performed on the nonlinear terms of the linear complex canonical form to decouple the hyperbolic unstable direction and the central direction of the center manifold, obtaining the final Hamiltonian function. Decoupling can enable the motion characteristics of the hyperbolic unstable direction and the central direction to be clearly shown separately, and corresponding characterization parameters can be defined and extracted respectively for the motion characteristics of different directions, making the description of the spacecraft motion more accurate and detailed. Referring to the linearly integrable part of the final Hamiltonian function, local action-angle variables are defined according to the motion mode of the libration point to describe the motion of the spacecraft on the center manifold, and a mapping relationship between the CRTBP coordinates and the characterization parameters is constructed for parameter characterization. By constructing the mapping relationship, the CRTBP coordinates are related to the characterization parameters, enabling the extraction of more physically meaningful and analytically valuable characterization parameters from the original coordinate information, thereby realizing the effective description and analysis of the spacecraft motion. Finally, the two motion characteristics of the hyperbolic unstable direction and the central direction and the corresponding characterization parameters are analyzed. The orbit of the target spacecraft can be determined by the phase difference of the angular variables in the characterization parameters, and the orbit change information of the target spacecraft can be obtained by the change of the hyperbolic motion component in the characterization parameters. The phase difference of the angular variables reflects the relative position and motion state of the spacecraft on the orbit, and different phase differences correspond to different orbit characteristics, and the orbit where the spacecraft is located can be determined by the phase difference. The change of the hyperbolic motion component is closely related to the orbit change behavior of the spacecraft. By analyzing its change, it can be judged whether the spacecraft enters or leaves the libration point periodic orbit and which branch of the invariant manifold it enters, thereby realizing the effective acquisition and analysis of the spacecraft orbit change information. Description of the Drawings

[0013] Figure 1 It is a schematic flow chart of a method for analyzing the characterization parameters of a collinear Earth-Moon libration point orbit in an embodiment;

[0014] Figure 2 It is a schematic diagram of a south Halo orbit at the L1 point obtained through numerical calculation in an embodiment;

[0015] Figure 3 It is a schematic diagram of the characterization parameters of the south Halo orbit at the L1 point in an embodiment;

[0016] Figure 4 It is a schematic diagram of the Poincaré section at the L1 point and the corresponding periodic / quasi-periodic orbit when C = 1 in another embodiment;

[0017] Figure 5 It is a schematic diagram of the Poincaré section at the L1 and L2 points at different energy levels and the corresponding orbit families in an embodiment; (a) is a schematic diagram of the Poincaré section at the L1 point at different energy levels and the corresponding orbit families, and (b) is a schematic diagram of the Poincaré section at the L2 point and the corresponding orbit families;

[0018] Figure 6 It is a schematic diagram of the stable / unstable manifolds of the Lyapunov orbit at the L1 point and the corresponding changes in characterization parameters in an embodiment. Detailed implementation manners

[0019] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0020] In one embodiment, as Figure 1 shown, a method for analyzing the characterization parameters of a collinear Earth-Moon libration point orbit is provided, including the following steps:

[0021] Step 102, obtain the motion equation of the spacecraft in the rendezvous coordinate system under the assumption of the circular restricted three-body problem, describe the circular restricted three-body problem as a chaotic Hamiltonian dynamical system, and transform the motion equation to a new coordinate system to obtain the Hamiltonian function of the corresponding system after the coordinate system transformation.

[0022] Under the assumption of the circular restricted three-body problem (CRTBP), the motion equation of the spacecraft in the rendezvous coordinate system can be expressed as:

[0023]

[0024] where ρ = (X, Y, Z) T, in the formula, length, time, and mass are normalized with the Earth-Moon distance, the reciprocal of the angular velocity of the Moon's revolution, and the sum of the Earth-Moon masses as the benchmarks respectively. Ω is the equivalent potential energy, expressed as

[0025]

[0026] where r1 and r2 are the distances from the spacecraft to the Earth and the Moon respectively, and μ = 0.012150568 is the normalized mass of the Moon. For the convenience of studying the dynamics near the libration points and to make the series expansion have better numerical characteristics, the origin of the coordinate system is translated from the center of mass of the Earth-Moon system to the corresponding collinear libration point L i as shown in Figure 2 and the length dimension is further normalized with the distance γ i from the libration point L i to the nearest large celestial body. Let X i be the position of the collinear libration point L i in the rendezvous coordinate system. Then the coordinates (x, y, z) in the new coordinate system and the coordinates (X, Y, Z) in the original rendezvous coordinate system have the following transformation relationship:

[0027]

[0028] The Hamiltonian function of the corresponding system after the coordinate transformation is

[0029]

[0030] where the momentum of the system is defined as

[0031]

[0032] Step 104: Use the Legendre expansion method to perform a series expansion on the non-linear terms in the Hamiltonian function to obtain a Hamiltonian function in polynomial form; take the quadratic term of the Hamiltonian function in polynomial form as the linearized model of the libration point in the circular restricted three-body problem.

[0033] In order to be able to perform a Lie transformation on this system to transform it into the Birkhoff–Gustavson normal form, it is necessary to perform a series expansion on the non-linear terms (1 - μ) / r1 and μ / r2 in the Hamiltonian function (4) to make it into the form of a polynomial. Here, the Legendre expansion method is applied:

[0034]

[0035] where P n is the Legendre polynomial of order n. By expanding the non-linear terms through the above transformation method, the Hamiltonian function can be transformed into the form of a polynomial of order n:

[0036]

[0037] Among them, H n is an nth-order homogeneous polynomial, and c n (μ) is a real coefficient, and its value depends on μ and the selected libration point L i .

[0038] The linearized model near the libration point can be given by the quadratic term of the Hamiltonian function (7), that is

[0039]

[0040] Step 106: Transform the linearized model into a real normal form according to the real linear symplectic transformation matrix and perform a complex transformation on the central part of the real normal form to obtain a linear complex normal form; perform a canonical transformation on the non-linear terms of the linear complex normal form to decouple the hyperbolic unstable direction and the central direction of the central manifold, and obtain the final Hamiltonian function; refer to the linearly integrable part of the final Hamiltonian function, define local action-angle variables according to the motion mode of the libration point to describe the motion of the spacecraft on the central manifold, and construct the mapping relationship between the CRTBP coordinates and the characterization parameters for parameter characterization.

[0041] The linearized dynamics of the libration point in CRTBP can be transformed into a real normal form by a real linear symplectic transformation matrix for Equation (8):

[0042]

[0043] The specific transformation process is as follows: The canonical equations of the system corresponding to Equation (8) are:

[0044]

[0045] Among them, J and M are defined by the following equations:

[0046]

[0047] I3 is the 3rd-order identity matrix, then the characteristic polynomial of the linear system (10) is

[0048]

[0049] Let η = λ 2 , then the roots of p(λ) = 0 are

[0050]

[0051] Since c2 > 1, then η1 < 0, η2 > 0, η3 > 0, which also indicates that the collinear libration points have a saddle×center×center dynamical structure. Define

[0052]

[0053] Suppose the matrix formed by scaling with the eigenvectors u i , v i of M is

[0054]

[0055] To make the real linear symplectic transformation matrix C satisfy the discriminant condition of the symplectic matrix: C T JC = J, the numerical values of the corresponding scaling coefficients s1, s2, s3 can be calculated. To sum up, the expression of the real linear symplectic transformation matrix is as follows:

[0056]

[0057] Therefore, through the transformation

[0058] (xyzp x p y p z ) T = C(q1q2q3p1p2p3) T (18)

[0059] The Hamiltonian function can be transformed into the real normal form (9). From the canonical equations, the linearized dynamical equation corresponding to H2 is:

[0060]

[0061] The solution of this system of equations is:

[0062]

[0063] Further transform H2 into the complex normal form. Since the saddle part already has a suitable form, only the complex transformation of the center part is needed here:

[0064]

[0065] Through the transformation, the Hamiltonian function (9) is further reduced to the linear complex normal form:

[0066] H2 = λq1p1 + iω p q2p2 + iω v q3p3 (22)

[0067] Now the action-angle variables can be defined: I j = q jp j ,θ j =arctan(q j / p j ). Therefore, the linear standard form is integrable with the integral of motion directly corresponding to the central manifold. In these action-angle variables, the equation of motion becomes a constant action and a linearly varying angle. In order to analyze the characteristics of the family of collinear libration point orbits, it is necessary to decouple the hyperbolic unstable direction of the central manifold from the central direction while maintaining the CRTBP dynamic characteristics as much as possible. Here, it is necessary to perform a canonical transformation on the nonlinear terms of the Hamiltonian function. The following describes how to use the Lie transformation method to deal with the high-order terms of the Hamiltonian function.

[0068] Simplification to the central manifold is a semi-standard form processing method that can simplify and analyze the nonlinear terms of the Hamiltonian function. This process requires a canonical transformation of the variables. A canonical transformation is a transformation that can maintain the form of the Hamiltonian canonical equation unchanged. Compared with the 6-dimensional differential equations used to study CRTBP, only one-dimensional Hamiltonian functions need to be studied after the canonical transformation. Usually, it is difficult to obtain a canonical transformation, and the generating function corresponding to the transformation needs to be calculated. The main methods for determining the generating function are implicit transformation (von Zeipel transformation) and explicit transformation (Lie transformation). In order to facilitate coordinate transformation, the Lie transformation method is used here.

[0069] Lie transform is a commonly used technique that can simplify the dynamics while maintaining accuracy. The theory of Lie transform relies on the following two facts: first, if (q, p) transforms to (Q, P) with the phase flow of the Hamiltonian system H over time t, then (q, p) → (Q, P) is a canonical transformation; second, assuming that a function f is a function of the phase flow (q, p) of the Hamiltonian system H, denoted by f(q, p), then the derivative of the function f with respect to t can be expressed as:

[0070]

[0071] The calculation of the Poisson bracket {F,G} is defined as

[0072]

[0073] Similarly, the second-order derivative of function f can be expressed as

[0074]

[0075] The n-fold Poisson bracket of function f in Hamiltonian system H is Then the nth-order derivative of f can be expressed as

[0076]

[0077] Suppose there is a Hamiltonian system G to be constructed, which can transform the original Hamiltonian function H(q, p) of the system into and further endow it with certain special properties (such as eliminating certain special terms). This process can be regarded as the canonical coordinates at time t = 0 moving to (q, p) at time t = 1 under the action of the Hamiltonian system G. By performing a Taylor series expansion at time t = 0, we can obtain the explicit transformation of:

[0078]

[0079] where i = 1, 2, 3, and the Hamiltonian function of the corresponding system is transformed into

[0080]

[0081] In this Lie transformation, the Hamiltonian system G is also called the generating function. It can be proved that if we want to achieve the inverse transformation we only need to select the generating function -G. Therefore, the core problem of the Lie transformation is how to select a suitable generating function to transform the Hamiltonian function of the system into a more concise form. In this problem, it is to eliminate the high-order terms in the Hamiltonian function. The following introduces how to implement this process.

[0082] Since we seek to decouple the hyperbolic unstable direction and the central direction in the high-order terms of the Hamiltonian function, we will choose generating functions G of different orders to keep the integrity of the saddle point part q1p1 in the high-order terms of the Hamiltonian function of the system. First, consider how to handle the third-order terms in the Hamiltonian function after performing the complex transformation (21). According to Equation (28), it can be seen that the transformed third-order terms are

[0083]

[0084] where G3 represents that the generating function G is a third-order homogeneous polynomial. Since H2 and H3 can be obtained from the original Hamiltonian function (7) after the complex transformation, and H2 is in the linear canonical form, so let we can get

[0085]

[0086] where h i.j is the coefficient of H3: The purpose of this application is to decouple the hyperbolic direction and the central direction, so we only need to eliminate Specific terms in it. The method of eliminating terms is not unique. Here we choose to keep the saddle point part as a whole q1p1, that is, eliminate the terms where the powers of q1 and p1 are not equal. Therefore, when calculating the generating function G3 according to (30), adding the restriction condition i1≠j1 can achieve the above purpose. After determining the generating function G3, the third-order terms of the transformed Hamiltonian function already meet the requirement that the powers of q1 and p1 in each term are equal. However, the terms higher than the third order will become new forms due to the transformation of G3. Specifically, it can be calculated from equation (28):

[0087]

[0088] Next, only need to determine G n (n≥4) in turn to process the higher-order terms to separate the hyperbolic unstable direction from the central direction. Similarly, the calculation method of the higher-order generating function is

[0089]

[0090] where i1≠j1. Implement the above process to N - order accuracy, and the Hamiltonian function will be transformed into

[0091]

[0092] where R N (q,p) is the remainder higher than N - order, is the term between the second order and the N - order, and all are in the form of q1p1. For the convenience of analysis, perform the inverse transformation of the complex transformation (21) on equation (33), and the complex variables can be converted into real variables. So far, through the canonical transformation of the Hamiltonian function, the decoupling of the hyperbolic unstable direction and the central direction within N - order accuracy has been achieved.

[0093] On the basis of completing the decoupling in the previous section, referring to the linearly integrable part of the Hamiltonian function, define the local action - angle variables according to the motion mode of the equilibrium point. Different ways are adopted to define the central integral and the saddle - point integral:

[0094]

[0095] where the subscript c represents the central (center) motion mode, and s represents the saddle (saddle) motion mode. For the collinear equilibrium points, their motion mode is saddle×center×center. Therefore, for the definition of the motion integral and the angle variables (I1,θ1,I2,θ2,I3,θ3), the definition method can adopt (I s ,θ s ) to define (I1,θ1), (I c ,θ cDefine \((I_2, \theta_2)\) and \((I_3, \theta_3)\). However, when selecting the parameters of the orbital representation, it is more inclined to retain \(q_1\) and \(p_1\) rather than representing them with their integrals of motion and angular variables. The main reasons are as follows: First, the definition method of the angular variables at the saddle point involves the generation of complex variables. Although the quadrants of its complex variables contain certain physical meanings, they are relatively abstract in practical applications. For specific references, please refer to the description in Appendix

[26] ; Second, the above physical meanings can be fully represented by \(q_1\) and \(p_1\) alone, that is, the degree to which the spacecraft cuts into the unstable and stable manifolds. It can be seen from the solution (20) of the collinear libration point linearized dynamics that the change of \(q_1\) rises exponentially, and the change of \(p_1\) decreases exponentially. Their changes correspond to the trend of the target moving along the invariant manifold: when the target leaves the libration point along the unstable manifold, the value of \(q_1\) will rise rapidly; when the target enters the vicinity of the libration point along the stable manifold, the value of \(p_1\) will approach 0 infinitely. And \(q_2, p_2, q_3, p_3\) correspond to the motion law of the target on the central manifold, and their integrals of motion and angular variables can be defined by \((I c , \theta c ) and have certain physical meanings: \(I_2\) represents the motion amplitude of the target in the XY plane of the rendezvous coordinate system, \(I_3\) represents the motion amplitude of the target in the Z direction, and \(\theta_2, \theta_3\) have physical meanings similar to phases. When \(\theta_2, \theta_3\) have specific relationships, periodic / quasi-periodic orbits will be generated. Therefore, make the following transformation:

[0096]

[0097] where \(j = 2, 3\). The Hamiltonian function at this time can be expressed as

[0098] \(H = H_2(q_1p_1, I_2, I_3)+H N (q_1p_1, I_2, I_3, \theta_2, \theta_3)+R N (36)

[0099] where \(R N is the remainder term of higher than order N, and \(H N is the term between order 2 and order N. To sum up, the catalog parameters of the target are selected as \((q_1, p_1, I_2, \theta_2, I_3, \theta_3)\). According to the previous derivation process, the mapping relationship between the state of the target in the rendezvous coordinate system and the catalog parameters can be constructed:

[0100]

[0101] The process of mutual conversion is further summarized as follows:

[0102]

[0103] Among them, A represents the process of coordinate system translation and scaling shown in formula (3), B represents the process of constructing a Hamiltonian system and redefining the generalized coordinates and generalized momenta shown in formula (5), C represents the process of transforming the linear term of the Hamiltonian function into a real canonical form shown in formula (18), D represents the complex transformation process shown in formula (21), and G n represents the process of eliminating the nth-order specific term through a canonical transformation shown in (32), and E represents the process of defining the characterization parameters according to (35). Since each of the above transformation processes is reversible, the transformation process from the characterization parameters to the coordinates in the rendezvous coordinate system can be similarly deduced to realize parameter characterization.

[0104] Step 108: Analyze the two motion characteristics of the hyperbolic unstable direction and the central direction and their corresponding characterization parameters. The orbit of the target spacecraft can be determined by the phase difference of the angular variables in the characterization parameters, and the orbit change information of the target spacecraft can be obtained from the change of the hyperbolic motion component in the characterization parameters; the orbit change information includes whether the target spacecraft enters / leaves the libration point periodic orbit and which invariant manifold it enters.

[0105] After decoupling the hyperbolic unstable direction and the central direction, analyze the two motion characteristics and their corresponding characterization parameters respectively. The last step transformation in (38) does not change the structure of the hyperbolic components q1 and p1, so their variation law still satisfies the canonical equations:

[0106]

[0107] If only the motion characteristics of the central direction are analyzed, the terms containing hyperbolic motion components can be simply excluded, that is, let q1p1 = 0. Ignoring the truncation terms higher than the Nth order, the Hamiltonian function can be expressed as

[0108] H CM = H2(I2, I3)+H N (I2, I3, θ2, θ3) (40)

[0109] At this time, the Hamiltonian function is only a function of the action-angle variables. According to the canonical transformation theory, the derivatives of the action-angle variables are

[0110]

[0111] where j = 2, 3. It can be seen that H2(I2, I3) is the linear system part of the original system, and this part is integrable. Since it does not explicitly contain the angular variables θ2, θ3, there is {I j , H2} = 0 in the linear part, that is, there is a cyclic integral I j is a constant. However, only the hyperbolic unstable direction and the central direction are separated in the high-order terms, so the high-order terms do not satisfy the integrable property. Therefore, Ij The change can be regarded as adding a perturbation term with long-term oscillation on the basis of the integration constant of the linear system. Similarly, for the angular variable, the linear part gives the law of linear change of the angular quantity: θ = ωt + θ0, and the high-order term will cause small-amplitude vibrations of the actual angular quantity on the basis of the linear change.

[0112] To verify the accuracy of the above results, a numerical method was used to calculate and test a Halo orbit south of the L1 point, as Figure 2 shown. The following method was used for the test: The states at different positions on the Halo orbit were selected and transformed into characterization parameters, and these were used as a benchmark. Since the change of the parameters is fundamentally governed by the CRTBP dynamic equation, as Figure 3 shown by the blue curve in. Next, the initial state of the Halo orbit was selected and transformed into a characterization parameter. Using this parameter as the initial value, numerical integration was performed using (39), (41), and (42), and the result was the red curve in Figure 3 . CM represents the central manifold. Since the change of the parameters is essentially governed by the canonical equation corresponding to the central manifold. It can be seen that among the characterization parameters corresponding to the central direction, the operation result of CRTBP is consistent with the operation result of the central manifold canonical equation, proving the consistency of the parameter transformation. However, there are differences in the parameters (q1 and p1) in the hyperbolic direction. This is because the numerical accuracy is limited when calculating the Halo orbit, and there are truncation errors higher than the Nth order when expanding the Hamiltonian function. Therefore, there is a deviation of about 10 -4 orders of magnitude in q1 and p1 during the parameter conversion process (in an ideal periodic orbit, since there is no motion in the hyperbolic direction, q1 and p1 should remain 0). Therefore, under the action of (39), q1 and p1 are amplified / diminished exponentially. This also reflects the necessity of decoupling the hyperbolic direction and the central direction in the previous section: The rapid change of the hyperbolic direction parameters does not affect the change of the central direction parameters. Therefore, it has a certain robustness in the presence of a series of errors. The mission analysts can focus on the evolution law of the characterization parameters in the central direction and study the characteristics such as the period and amplitude of the orbit family, which is also the aspect that target cataloging pays more attention to.

[0113] The motion of the spacecraft in the central direction can be regarded as the motion on a four-dimensional phase space manifold. To facilitate the selection of the integration initial value, the method of Poincaré section is used to reduce the dimension of the phase space. Here, first, θ2 = 0 is selected as the section because it is found through simulation that most of the phase flow of the central manifold intersects with this section. Since the Hamiltonian function (40) does not explicitly contain time t, the Hamiltonian (that is, the energy in the broad sense) is a conserved quantity during the motion. Fixing the energy level C can obtain the corresponding Poincaré section:

[0114] H CM (I2, 0, I3, θ3) = C (43)

[0115] By uniformly selecting lattice points on the cross-section to calculate the initial values corresponding to the center manifold and integrating backward, it can be found that the phase flow will cross the Poincaré cross-section multiple times, and the crossing positions are Figure 4 marked. These positions together will form a cross-sectional diagram. Orbits of different orbit families will form cross-sectional diagrams with different characteristics: the cross-sectional diagram of the Lissajous orbit will pass through the entire cross-section as θ3 varies; the cross-section corresponding to the Lyapunov orbit intersects with the I2 - θ3 plane, and the cross-section corresponding to the vertical Lyapunov orbit intersects with the I3 - θ3 cross-section; the quasi-Halo orbit will form a ring-like structure locally; in particular, when the ring shrinks continuously to the center point, a Halo orbit will be formed; during the process of θ3 varying in [0, 2π), two ring-like structures will appear, and their centers are located at θ3 = π / 2 and θ3 = 3π / 2 respectively. By transforming to the rendezvous coordinate system, it can be found that the Halo orbit families corresponding to these two ring-like structures are the north family and the south family orbits respectively. The north family and the south family orbits are symmetric about the x - y plane in the rendezvous coordinate system, and in the center manifold coordinates, they correspond to phases with θ3 differing by π. The above is the Poincaré cross-section formed at a fixed energy level. If the energy changes, the cross-section will change, and the cross-sectional diagrams formed by each orbit family on the cross-section will also change. For further convenience of representation, a plane that can contain all orbit families (especially the centers of the ring-like structures, that is, the Halo orbit family) is further selected to intercept the cross-section of the above three-dimensional space. Here we choose the cross-section of θ3 = π / 2 and draw the situation including all energy levels, forming a Poincaré cross-sectional diagram that includes all orbit families near the equilibrium point, as Figure 5 shown. It can be seen that the Lyapunov orbit family and the vertical Lyapunov orbit family correspond to the cases of I3 = 0 and I2 = 0 respectively. The asymptotic Lyapunov orbit (referred to as Large quasi-Halo in some literature) is the boundary between the Lissajous orbit and the quasi-Halo orbit. When the quasi-Halo orbit meets specific conditions, a Halo orbit will be formed, and its position is Figure 5It has been marked accordingly. Different energy levels C are marked by contour lines in the figure. It can be seen that when the energy level is low, there is no Halo orbit, which also reflects the origin of the Halo orbit: the Halo orbit bifurcates during the extension of the Lyapunov orbit from a low energy level to a high energy level. Since the cross-sectional diagrams of θ3 = π / 2 and θ3 = 3π / 2 are exactly the same, the Halo orbits of the north family and the south family cannot be distinguished in the figure, which is also a manifestation of the symmetry between the north family orbit and the south family orbit. In practical applications, if it is necessary to determine whether the target is specifically on the Halo orbit of the south family or the north family, it is only necessary to observe the phase difference between θ3 and θ2 to make a judgment.

[0116] Since the process of reducing to the center manifold has decoupled the hyperbolic unstable direction and the center direction through a canonical transformation, it is possible to introduce a small unstable component ε in the q1 and p1 directions and obtain the initial integration state in the rendezvous coordinate system through a coordinate transformation. According to the definitions of the stable manifold and the unstable manifold, the stable and unstable manifolds in the center manifold coordinates can be defined as

[0117] W s ={(q1,p1,I2,θ2,I3,θ3)|H CM =C, q1 = 0} (44)

[0118] W u ={(q1,p1,I2,θ2,I3,θ3)|H CM =C, p1 = 0} (45)

[0119] Among them, W s represents the stable manifold, and W u represents the unstable manifold. When q1 = 0 and p1 ≠ 0, p1 will exponentially approach 0 over time, and the corresponding stable manifold will infinitely approach the libration point orbit. Similarly, when p1 = 0 and q1 ≠ 0, q1 will exponentially increase over time, and the corresponding unstable manifold will gradually move infinitely away from the libration point orbit. The initial state of the libration point orbit in the center manifold coordinates can be expressed as

[0120] σ center =(00I2θ2I3θ3) T (46)

[0121] If it is necessary to calculate the unstable manifold or the stable manifold, it is only necessary to add an unstable component ε to the corresponding characterization parameter q1 or p1. Taking the unstable manifold as an example, that is

[0122] σ u =σ center +(±ε00000) T (47)

[0123] Among them, the sign of ε represents the unstable manifolds in different directions of the left and right branches. Then σ u is transformed into the rendezvous coordinate system and integrated forward with it as the initial value. Similarly, after adding an unstable component to p1 and transforming it into the rendezvous coordinate system, the stable manifold can be obtained by integrating it backward with it as the initial value.

[0124] This method is reminiscent of the process of calculating the stable and unstable manifolds of Lyapunov orbits or Halo orbits: First, calculate the Monodromy matrix of the periodic orbit, and select the corresponding eigenvectors according to the eigenvalues of the Monodromy matrix. The direction of the eigenvector corresponding to the eigenvalue greater than 1 is the direction of the unstable manifold, and the direction of the eigenvector corresponding to the eigenvalue less than 1 is the direction of the stable manifold. However, this method has certain limitations. When it comes to quasi-periodic orbits, the Monodromy matrix cannot be clearly defined. The method described in this application can make up for the above deficiencies. In addition to calculating the stable and unstable manifolds, the characterization parameters in the hyperbolic direction can also be used as key parameters for monitoring the entry / exit of a spacecraft from / to a libration point. As Figure 6 shown, it is the change of the stable / unstable manifolds and the corresponding characterization parameters of the Lyapunov orbit at the L1 point. Taking the Lyapunov orbit at the L1 point as an example, mission analysts only need to observe the changes of the characterization parameters q1 and p1 in the hyperbolic direction to obtain information such as whether the target spacecraft enters / leaves the libration point periodic orbit and which invariant manifold it enters. Considering fuel savings in low-energy transfers, most transfers of spacecraft in the Earth-Moon space are achieved using invariant manifolds. Therefore, this parameter has important reference value for judging events such as spacecraft orbit changes and transfers. This is also the reason for directly selecting q1 and p1 instead of their action-angle variables as the characterization parameters.

[0125] The above method for analyzing the characterization parameters of the Earth-Moon collinear libration point orbit first obtains the motion equation of the spacecraft in the rendezvous coordinate system under the assumption of the circular restricted three-body problem, and describes it as a chaotic Hamiltonian dynamical system. The motion equation is transformed to a new coordinate system to obtain the corresponding Hamiltonian function, which is convenient for subsequent transformations and analyses. In this way, the motion of the spacecraft is related to the Hamiltonian function, enabling the study of the motion characteristics of the spacecraft from the perspective of the Hamiltonian function and laying a foundation for extracting the characterization parameters. The Legendre expansion method is used to expand the nonlinear terms in the Hamiltonian function into a series, transforming the complex nonlinear relationship into a polynomial form for easy analysis and processing. The quadratic term of the Hamiltonian function in polynomial form serves as the linearized model of the libration point in the circular restricted three-body problem. Near the libration point, the quadratic term can well approximate the local behavior of the system. The linearized model is relatively simple and easy to analyze. According to the real linear symplectic transformation matrix, the linearized model is transformed into a real canonical form, and a complex transformation is performed on the central part of the real canonical form to obtain a linear complex canonical form, transforming the system equation into a form more convenient for analysis. The real canonical form and the linear complex canonical form have specific structures and properties, which can clearly separate the different motion characteristics of the system, enabling targeted analysis of the motion characteristics of different parts. A canonical transformation is performed on the nonlinear terms of the linear complex canonical form to decouple the hyperbolic unstable direction and the central direction of the center manifold, obtaining the final Hamiltonian function. Decoupling can enable the motion characteristics of the hyperbolic unstable direction and the central direction to be clearly displayed separately, and corresponding characterization parameters can be defined and extracted respectively for the motion characteristics of different directions, making the description of the spacecraft motion more accurate and detailed. Referring to the linearly integrable part of the final Hamiltonian function, local action-angle variables are defined according to the motion mode of the libration point to describe the motion of the spacecraft on the center manifold, and a mapping relationship between the CRTBP coordinates and the characterization parameters is constructed for parameter characterization. By constructing the mapping relationship, the CRTBP coordinates are related to the characterization parameters, enabling the extraction of more physically meaningful and analytically valuable characterization parameters from the original coordinate information, thereby realizing the effective description and analysis of the spacecraft motion. Finally, the motion characteristics of the hyperbolic unstable direction and the central direction and the corresponding characterization parameters are analyzed. The orbit of the target spacecraft can be determined by the phase difference of the angular variables in the characterization parameters, and the orbit change information of the target spacecraft can be obtained by the change of the hyperbolic motion component in the characterization parameters. The phase difference of the angular variables reflects the relative position and motion state of the spacecraft on the orbit. Different phase differences correspond to different orbit characteristics, and the orbit where the spacecraft is located can be determined by the phase difference. The change of the hyperbolic motion component is closely related to the orbit change behavior of the spacecraft. By analyzing its change, it can be judged whether the spacecraft enters or leaves the libration point periodic orbit and which invariant manifold it enters, thereby realizing the effective acquisition and analysis of the spacecraft orbit change information.

[0126] In one embodiment, the equations of motion are transformed into a new coordinate system to obtain the Hamiltonian function of the corresponding system after the coordinate transformation, including:

[0127] The equations of motion are transformed into a new coordinate system, and the Hamiltonian function of the corresponding system after the coordinate transformation is:

[0128]

[0129] where r1 and r2 are the distances from the spacecraft to the Earth and the Moon respectively, μ = 0.012150568 is the normalized mass of the Moon, and the momentum of the system is defined as:

[0130]

[0131] (x, y, z) represents the coordinates in the new coordinate system.

[0132] In one embodiment, the linearized model is transformed into a real canonical form according to the real linear symplectic transformation matrix, including:

[0133] The linearized model is transformed into a real canonical form according to the real linear symplectic transformation matrix as:

[0134]

[0135] where, (xyzp x p y p z ) T = C(q1q2q3p1p2p3) T , C represents the real linear symplectic transformation matrix, λ, ω p , ω v represent the coefficients of each term in the canonical form after the symplectic transformation.

[0136] In one embodiment, a complex transformation is performed on the central part of the real canonical form to obtain a linear complex canonical form, including:

[0137] A complex transformation is performed on the central part of the real canonical form, and the linear complex canonical form obtained is:

[0138] H2 = λq1p1 + iω p q2p2 + iω v q3p3.

[0139] In one embodiment, a canonical transformation is performed on the non - linear terms of the linear complex canonical form to decouple the hyperbolic unstable direction and the central direction of the center manifold, obtaining the final Hamiltonian function, including:

[0140] Perform a canonical transformation on the nonlinear terms of the linear complex normal form to decouple the hyperbolic unstable direction and the center direction of the center manifold, and obtain the final Hamiltonian function as

[0141]

[0142] where R N (q, p) is the remainder term greater than the Nth order, is the term between the 2nd order and the Nth order, and all are in the form of q1p1. R N .

[0143] In one embodiment, the Hamiltonian function describing the motion of the spacecraft on the center manifold is:

[0144] H = H2(q1p1, I2, I3) + H N (q1p1, I2, I3, θ2, θ3) + R N

[0145] where H2 represents the linear complex normal form, and R N (q, p) is the remainder term greater than the Nth order, is the term between the 2nd order and the Nth order, and all are in the form of q1p1. R N .

[0146] In one embodiment, construct the mapping relationship between the CRTBP coordinates and the characterization parameters, including:

[0147] Select the catalog parameters of the target as (q1, p1, I2, θ2, I3, θ3), and construct the mapping relationship between the state of the target in the rendezvous coordinate system and the catalog parameters:

[0148]

[0149] where (X, Y, Z) represents the position of the spacecraft in the rendezvous coordinate system, represents the velocity of the spacecraft in the rendezvous coordinate system, q1 and p1 respectively characterize the degree of the spacecraft cutting into the unstable manifold and the stable manifold, I2 represents the motion amplitude of the target in the XY plane of the rendezvous coordinate system, I3 represents the motion amplitude of the target in the Z direction, and θ2, θ3 have the physical meaning of phase. When θ2, θ3 have a specific relationship, periodic / quasi-periodic orbits will be generated; the physical meaning of the phase includes the variation law of the generalized coordinates and the generalized momenta.

[0150] In one embodiment, analyze the two motion characteristics of the hyperbolic unstable direction and the center direction and the corresponding characterization parameters. The orbit of the target spacecraft can be determined by the phase difference of the angular variables in the characterization parameters, and the orbit change information of the target spacecraft can be obtained by the change of the hyperbolic motion components in the characterization parameters, including:

[0151] Analyze the motion characteristics in the central direction and the corresponding characterization parameters, exclude the terms containing hyperbolic motion components, that is, let q1p1 = 0, and ignore the truncation terms higher than the Nth order. The Hamiltonian function is expressed as:

[0152] H CM = H2(I2, I3) + H N (I2, I3, θ2, θ3)

[0153] The motion of the spacecraft in the central direction is regarded as the motion on a four-dimensional phase space manifold. The Poincaré section method is used to reduce the dimension of the phase space. Select θ2 = 0 as the section. Since the Hamiltonian function does not explicitly contain time t, the Hamiltonian is a conserved quantity during the motion. Fixing the energy level C, the corresponding Poincaré section can be obtained as:

[0154] H CM (I2, 0, I3, θ3) = C

[0155] Uniformly select lattice points on the section to calculate the initial values corresponding to the central manifold and integrate backward. The phase flow will cross the Poincaré section multiple times, causing the Halo orbit to form a ring-like structure locally. When the ring shrinks continuously to the center point, a Halo orbit is formed; during the process of θ3 changing in [0, 2π), two ring-like structures appear, and their centers are located at θ3 = π / 2 and θ3 = 3π / 2 respectively. By transforming to the rendezvous coordinate system, the Halo orbit families corresponding to these two ring-like structures are determined to be the north family and the south family orbits. The north family and the south family orbits are symmetric about the x - y plane in the rendezvous coordinate system, and in the central manifold coordinates, they correspond to a phase difference of π in θ3. Therefore, the orbit of the target spacecraft can be determined by the phase difference of the angular variables in the characterization parameters.

[0156] In one embodiment, according to the definitions of the stable manifold and the unstable manifold, the stable and unstable manifolds in the central manifold coordinates are defined as

[0157] W s ={(q1, p1, I2, θ2, I3, θ3)|H CM = C, q1 = 0}

[0158] W u ={(q1, p1, I2, θ2, I3, θ3)|H CM = C, p1 = 0}

[0159] Among them, W s represents the stable manifold, and W u represents the unstable manifold;

[0160] When q1 = 0 and p1 ≠ 0, p1 will exponentially approach 0 over time, and the corresponding stable manifold will infinitely approach the libration point orbit. Similarly, when p1 = 0 and q1 ≠ 0, q1 will exponentially increase over time, and the corresponding unstable manifold will gradually move infinitely away from the libration point orbit. The initial state of the libration point orbit in the central manifold coordinates can be expressed as

[0161] σ center =(0 0 I2 θ2 I3 θ3) T

[0162] After adding the unstable component ε to the corresponding characterization parameter q1 or p1, σ u is transformed to the rendezvous coordinate system, and integrating it forward with this as the initial value can obtain the unstable manifold. Similarly, after adding the unstable component to p1 and transforming it to the rendezvous coordinate system, integrating it backward with this as the initial value can obtain the stable manifold. Therefore, the orbit change information of the target spacecraft can be obtained by the change of the hyperbolic motion component in the characterization parameter.

[0163] It should be understood that although Figure 1 the steps in the flowchart of Figure 1 are shown in sequence according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise clearly stated in this article, there is no strict order restriction for the execution of these steps, and these steps can be executed in other orders. Moreover,[[]] Figure 1 at least a part of the steps in

[0164] can include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily executed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential either, but can be executed alternately or in turn with at least a part of other steps or sub-steps or stages of other steps.

[0164] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as the scope recorded in this specification.

[0165] The above-described embodiments only represent several implementation manners of the present application. The description is relatively specific and detailed, but it cannot be understood as a limitation to the scope of the invention. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present application, several modifications and improvements can still be made, and these all belong to the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the appended claims.

Claims

1. A method for characterizing parameters of Earth-Moon collinear libration point orbits, characterized in that: The method comprises: Obtain the motion equation of the spacecraft in the rendezvous coordinate system under the assumption of a circular restricted three-body problem, describe the circular restricted three-body problem as a chaotic Hamiltonian dynamic system, transform the motion equation into a new coordinate system, and obtain the Hamiltonian function of the corresponding system after the coordinate system transformation; The nonlinear terms in the Hamiltonian function are expanded in series by using the Legendre expansion method to obtain a Hamiltonian function in polynomial form; the quadratic terms of the Hamiltonian function in polynomial form are used as a linearized model of the translation point in the circular restricted three-body problem; The linearized model is transformed into a real standard form according to a real linear symplectic transformation matrix and a central part of the real standard form is complex transformed to obtain a linear complex standard form; a canonical transformation is performed on the nonlinear terms of the linear complex standard form to achieve decoupling of the hyperbolic unstable direction of the central manifold from the central direction to obtain a final Hamiltonian function; with reference to the linear integrable part of the final Hamiltonian function, a local action-angle variable is defined according to the motion mode of the libration point to describe the motion of the spacecraft on the central manifold, and a mapping relationship between CRTBP coordinates and characterization parameters is constructed for parameter characterization; The two motion characteristics of the hyperbolic unstable direction and the central direction and the corresponding characterization parameters are analyzed. The orbit of the target spacecraft can be determined by the phase difference of the angular variable in the characterization parameters, and the orbit change information of the target spacecraft can be obtained by the change of the hyperbolic motion component in the characterization parameters; the orbit change information includes whether the target spacecraft enters / leaves the periodic orbit of the libration point and which invariant manifold it enters.

2. The method according to claim 1, characterized in that The motion equation is transformed into a new coordinate system to obtain the Hamiltonian function of the corresponding system after the coordinate system transformation, including: The motion equation is transformed into the new coordinate system, and the Hamiltonian function of the corresponding system after the coordinate system transformation is obtained: Among them, r1 and r2 are the distances from the spacecraft to the earth and the moon, respectively, μ = 0.012150568 is the normalized mass of the moon, and the momentum of the system is defined as: (x, y, z) represents the coordinates in the new coordinate system.

3. The method according to claim 1, characterized in that Transforming the linearized model into a real standard form according to a real linear symplectic transformation matrix includes: The linearized model is transformed into a real standard form according to the real linear symplectic transformation matrix: Among them, (xyzp x p y p z ) T =C(q1 q2 q3 p1 p2 p3) T , C represents the real linear symplectic transformation matrix, λ, ω p ,ω v Represents the coefficients of the standard form terms after the symplectic transformation.

4. The method according to claim 3, characterized in that A central part of the real standard form is subjected to a complex transformation to obtain a linear complex standard form, including: The central part of the real standard form is complex transformed to obtain the linear complex standard form: H2=λq1p1+iω p q2p2+iω v q3p3.

5. The method according to claim 4, characterized in that The nonlinear terms of the linear complex standard form are subjected to canonical transformation to achieve decoupling of the hyperbolic unstable direction of the central manifold from the central direction, and the final Hamiltonian function is obtained, including: The nonlinear terms of the linear complex standard form are canonically transformed to decouple the hyperbolic unstable direction of the central manifold from the central direction, and the final Hamiltonian function is obtained as Among them, R N (q,p) is a remainder greater than order N, is a term between 2nd and Nth order, and all are of the form q1p1R N .

6. The method according to claim 1, characterized in that The Hamiltonian function describing the motion of the spacecraft on the central manifold is: H=H2(q1p1,I2,I3)+H N (q1p1,I2,I3,θ2,θ3)+R N Among them, H2 represents the linear complex standard form, R N (q,p) is a remainder greater than order N, is a term between 2nd and Nth order, and all are of the form q1p1R N .

7. The method according to claim 1, characterized in that Construct the mapping relationship between CRTBP coordinates and characterization parameters, including: The target catalog parameters are selected as (q1, p1, I2, θ2, I3, θ3), and the mapping relationship between the target state in the rendezvous coordinate system and the catalog parameters is constructed: Among them, (X, Y, Z) represents the position of the spacecraft in the rendezvous coordinate system, It represents the speed of the spacecraft in the rendezvous coordinate system, q1 and p1 respectively characterize the degree of the spacecraft's entry into the unstable manifold and the stable manifold, I2 represents the motion amplitude of the target in the XY plane of the rendezvous coordinate system, I3 represents the motion amplitude of the target in the Z direction, θ2 and θ3 have the physical meaning of phase. When θ2 and θ3 have a specific relationship, a periodic / quasi-periodic orbit will be generated; the physical meaning of the phase includes the changing laws of generalized coordinates and generalized momentum.

8. The method according to claim 1, characterized in that The two motion characteristics of the hyperbolic unstable direction and the central direction and the corresponding characterization parameters are analyzed. The orbit of the target spacecraft can be determined by the phase difference of the angular variable in the characterization parameters, and the orbit change information of the target spacecraft can be obtained by the change of the hyperbolic motion component in the characterization parameters, including: The motion characteristics of the center direction and the corresponding characterization parameters are analyzed, and the terms containing the hyperbolic motion components are excluded, that is, q1p1=0, and the truncation terms greater than order N are ignored. The Hamiltonian function is expressed as: H CM =H2(I2,I3)+H N (I2,I3,θ2,θ3) The motion of the spacecraft in the center direction is regarded as the motion on the four-dimensional phase space manifold. The Poincare section method is used to reduce the dimension of the phase space, and θ2 = 0 is selected as the section. Since the Hamiltonian function does not explicitly contain time t, the Hamiltonian is a conserved quantity during the motion process. The corresponding Poincare section can be obtained for a fixed energy level C: H CM (I2,0,I3,θ3)=C By uniformly selecting point arrays on the cross section to calculate the initial values ​​corresponding to the central manifold and integrating backward, the phase flow will cross the Poincare cross section multiple times, so that the Halo orbit forms a ring-shaped structure locally. When the ring is continuously reduced to the center point, a Halo orbit is formed. In the process of θ3 changing in [0, 2π), two ring structures appear, and their centers are located at θ3=π / 2 and θ3=3π / 2 respectively. By transforming to the conjunction coordinate system, it is determined that the Halo orbit families corresponding to these two ring structures are the northern and southern orbits respectively. The northern and southern orbits are symmetric about the xy plane in the conjunction coordinate system, and the corresponding phase of θ3 in the central manifold coordinate is π. Therefore, the orbit of the target spacecraft can be determined by the phase difference of the angular variable in the characterization parameter.

9. The method according to claim 8, characterized in that The method further comprises: According to the definition of stable manifold and unstable manifold, the stable and unstable manifolds in the central manifold coordinates are defined as W s ={(q1,p1,I2,θ2,I3,θ3)|H CM =C,q1=0} W u ={(q1,p1,I2,θ2,I3,θ3)|H CM =C,p1=0} Among them, W s represents a stable manifold, W u represents an unstable manifold; When q1=0,p1≠0,p1 will approach 0 exponentially with time, and the corresponding stable manifold will be infinitely close to the orbit of the libration point. Similarly, when p1=0,q1≠0,q1 will increase exponentially with time, and the corresponding unstable manifold will gradually move away from the orbit of the libration point. The initial state of the orbit of the libration point in the coordinates of the central manifold can be expressed as s center =(00I2θ2I3θ3) T After adding the unstable component ε to the corresponding characterization parameter q1 or p1, σ u Transform it to the rendezvous coordinate system and use it as the initial value for forward integration to obtain the unstable manifold. Similarly, add the unstable component to p1 and transform it to the rendezvous coordinate system. Use it as the initial value for reverse integration to obtain the stable manifold. Therefore, the trajectory change information of the target spacecraft is obtained by characterizing the changes in the hyperbolic motion components in the parameters.

Images