Orbit identification method based on earth-moon collinear translation point orbit parameter characterization

Abstract

The invention relates to an orbit identification method based on earth-moon collinear translation point orbit parameter characterization. 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; and carrying out complex transformation and regular transformation on the linearized model transformation to realize the decoupling of the hyperbolic unstable direction and the central direction of the central manifold, and constructing a mapping relation between CRTBP coordinates and characterization parameters to carry out parameter characterization. Selecting spacecraft coordinates from parameter characterization to calculate characterization parameters of a reference orbit and integral initial values of the coordinates, and then calculating mean square errors of a real orbit and the reference orbit; and by taking the mean square error as a target function and a preset constraint condition, constructing an orbit identification solving model and solving to obtain an orbit identification result. By adopting the method, the track identification accuracy can be improved.

Description

Technical Field

[0001] The present application relates to the technical field of space orbit recognition, and particularly to an orbit recognition method based on the characterization of the orbital parameters of the Earth-Moon collinear libration point orbit. Background Art

[0002] In recent years, with the development of lunar exploration activities 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 point is an important resource for future Earth-Moon space development.

[0003] Space object 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 object cataloging is to maintain a detailed database to record information such as the orbital parameters, size, shape, and use of these objects, so as to conduct space situation awareness, collision warning, space traffic management, etc. The existing space object cataloging methods mainly form a recognized standard format based on the Keplerian orbital elements: Two-Line Orbital Element (TLE). In addition, there are also methods using regular Delaunay variables to describe orbital motion. Considering that in the perturbed two-body problem, the Keplerian orbital elements are no longer integrals of motion, 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 judge the orbit of a spacecraft. To better describe the motion of a target at the Earth-Moon libration point, the Circular Restricted Three-Body Problem (CRTBP) model is usually adopted, but this brings great difficulties to the cataloging of Earth-Moon libration point targets. On the one hand, the two-body problem is integrable, while the restricted three-body problem is non-integrable. This makes it impossible to extract more first integrals for the restricted three-body problem in the Cartesian coordinate system except for energy conservation, and thus it does not have a cataloging basis similar to the two-body problem in theory. On the other hand, for the special dynamic structures near the libration point, such as periodic orbits, quasi-periodic orbits, hyperbolic invariant manifolds, etc., there is currently no recognized parameterization description method, and its physical characteristics, such as amplitude, period, etc., are difficult to be directly related to 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. Just as there is a mapping relationship between Cartesian coordinates and Keplerian orbital elements, there is also a mapping relationship between Cartesian coordinates and local action-angle variables, which meets the parameter representation requirements of the Earth-Moon libration point targets. However, in the process of transforming the high-order expansion of dynamics into a normal form, the resonance terms of the Hamiltonian function are removed, and these resonance terms explain the bifurcation reasons of the halo orbit family. Therefore, some orbits cannot be explained by this high-order expansion method theoretically. At the same time, due to the elimination of many terms in the Hamiltonian function during the high-order expansion process, the effective approximation region of this method is small. For orbits in the region far from the libration points, it is difficult to maintain the integrable dynamic characteristics due to the influence of truncation errors. Therefore, when this method is applied to the cataloging work of collinear libration point targets, it only has a good cataloging effect on the periodic / quasi-periodic orbits in a small area near the libration points. When the orbit amplitude is large, the oscillation of the parameters will become very strong, resulting in low orbit recognition efficiency and accuracy. The complex dynamic model and the difficulty of parameterization not only make the orbit data processing time-consuming and laborious, but also lead to inaccurate recognition results, seriously affecting the application and development of the cataloging of space targets in the Earth-Moon libration point region. Summary of the Invention

[0006] Based on this, it is necessary to provide an orbit recognition method based on the parameter representation of the Earth-Moon collinear libration point orbit to improve the accuracy of orbit recognition for the above technical problems.

[0007] An orbit recognition method based on the parameter representation 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 dynamic system, transform the motion equation to a new coordinate system, and obtain the Hamiltonian function of the corresponding system after the coordinate system 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] The linearized model is transformed into a real canonical form according to the real linear symplectic transformation matrix, and a complex transformation is performed on the central part of the real canonical form to obtain a linear complex canonical form; a normal 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 central manifold, and the final Hamiltonian function is obtained; referring to the linearly integrable part of the final Hamiltonian function, local action-angle variables are defined according to the motion mode of the equilibrium point to describe the motion of the spacecraft on the central manifold, and the mapping relationship between the CRTBP coordinates and the characterization parameters is constructed for parameter characterization;

[0011] According to the Poincaré section, the characterization parameters of the reference orbit and the integral initial values of the coordinates are selected from the parameter characterization for the spacecraft coordinates, and the mean square error between the true orbit and the reference orbit is calculated using the characterization parameters of the reference orbit and the integral initial values of the coordinates; an orbit identification solution model is constructed with the mean square error as the objective function and the preset constraint conditions; the orbit identification solution model is solved according to the Bayesian optimization algorithm to obtain the orbit identification result.

[0012] The above-mentioned orbit recognition method based on the parameter characterization of the Earth-Moon collinear libration point orbit obtains the motion equation of the spacecraft in the rendezvous coordinate system, describes the circular restricted three-body problem as a chaotic Hamiltonian dynamical system, and then transforms it into the Hamiltonian function of the corresponding system in the new coordinate system. This transforms the originally complex and difficult-to-directly-analyze three-body problem into a form that can be deeply studied using Hamiltonian theory. Through this transformation, the motion law of the spacecraft can be systematically understood from the perspectives of physical quantities such as energy and momentum, laying a foundation for accurately identifying the orbit subsequently. The non-linear terms in the Hamiltonian function are expanded using the Legendre series to obtain the Hamiltonian function in polynomial form, and its quadratic term is used as the linearized model of the libration point. Through the linearized model, the motion of the spacecraft near the libration point can be initially approximately analyzed. The linearized model can quickly and relatively accurately describe its motion, providing an initial and relatively accurate reference for the overall orbit recognition. Then, according to the real linear symplectic transformation matrix, the linearized model is transformed into the real normal form, and then a complex transformation is performed on the central part of the real normal form to obtain the linear complex normal form. The normal form transformation can transform a complex linear system into a more standardized and easy-to-analyze form. Through this transformation, different motion modes and characteristics in the system can be more clearly separated, which is conducive to accurately judging the type of the spacecraft orbit and improving the recognition accuracy. The non-linear terms of the linear complex normal form are subjected to a canonical transformation to decouple the hyperbolic unstable direction and the central direction of the center manifold, obtaining the final Hamiltonian function. The decoupling operation enables the separate treatment of motions with different properties in the system. In the three-body problem, the motion on the center manifold contains the main motion characteristics of the spacecraft near the libration point, while the motion in the hyperbolic unstable direction may perturb the orbit. Through decoupling, these two motions with different properties can be accurately analyzed and described respectively. During orbit recognition, the stable motion part of the spacecraft on the center manifold and the part affected by the hyperbolic unstable direction can be more accurately grasped, thereby improving the recognition accuracy of the overall orbit. 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 the mapping relationship between the CRTBP coordinates and the characterization parameters is constructed. The local action-angle variables can more intuitively reflect the motion characteristics of the spacecraft near the libration point. By constructing the mapping relationship, the complex CRTBP coordinates are transformed into more physically meaningful characterization parameters, which can more accurately depict the characteristics of the orbit. According to the Poincaré section, the characterization parameters and the integral initial values of the coordinates of the reference orbit are selected from the parameter characterization for the spacecraft coordinates, and the mean square error between the real orbit and the reference orbit is calculated using these values. An orbit recognition solution model is constructed with the mean square error as the objective function and the preset constraint conditions, and then the model is solved according to the Bayesian optimization algorithm. The Poincaré section can transform the continuous orbit motion into a discrete point set, facilitating the selection of representative spacecraft coordinates from it.By calculating the mean square error, the difference between the true orbit and the reference orbit can be quantified. The smaller the mean square error, the closer the true orbit is to the reference orbit. Taking the mean square error as the objective function and combining the constraint conditions to construct a solution model, the Bayesian optimization algorithm can quickly and accurately search for the orbit parameters that best match the true orbit among many possible orbits, thereby obtaining an accurate orbit recognition result. This method based on the optimization algorithm and the quantified error can efficiently find the optimal solution in the complex three-body problem orbit space, greatly improving the accuracy of orbit recognition. Brief Description of the Drawings

[0013] Figure 1 FIG. is a schematic flow chart of an orbit recognition method based on the characterization of the collinear libration point orbit parameters in an embodiment;

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

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

[0016] Figure 4 FIG. 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 FIG. is a schematic diagram of the Poincaré section at the L1 point and the L2 point at different energy levels and the corresponding orbit families; (a) is a schematic diagram of the Poincaré section at the L1 point and the corresponding orbit families at different energy levels, and (b) is a schematic diagram of the Poincaré section at the L2 point and the corresponding orbit families;

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

[0019] Figure 7 FIG. is a schematic diagram of the true orbit and the reference orbit of the spacecraft at the L2 point in the rendezvous coordinate system in an embodiment;

[0020] Figure 8 FIG. is a schematic diagram of the characterization parameters corresponding to the true orbit and the reference orbit in an embodiment;

[0021] Figure 9 FIG. is a schematic diagram showing the Bayesian optimization process in an embodiment;

[0022] Figure 10 FIG. is a schematic diagram of the position of the reference orbit of the spacecraft in the Poincaré section at the L2 point in an embodiment. Detailed Description of the Invention

[0023] In order to make the objectives, technical solutions and advantages of the present application clearer and more understandable, the present application will be further described in detail below in conjunction with 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.

[0024] In one embodiment, as Figure 1 shown, an orbit recognition method based on the characterization of the collinear Earth-Moon libration point orbit parameters is provided, including the following steps:

[0025] 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, transform the motion equation to a new coordinate system, and obtain the Hamiltonian function of the corresponding system after the coordinate system transformation.

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

[0027]

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

[0029]

[0030] 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. In order to facilitate the study of the dynamics near the libration point and make the series expansion have good 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 Figure 2 shown, 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 coordinate (x, y, z) in the new coordinate system has the following transformation relationship with the coordinate (X, Y, Z) in the original rendezvous coordinate system:

[0031]

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

[0033]

[0034] Among them, the momentum of the system is defined as

[0035]

[0036] In order to perform a Lie transformation on this system to transform it into the Birkhoff–Gustavson normal form, it is necessary to expand the nonlinear terms \((1 - \mu) / r_1\) and \(\mu / r_2\) in the Hamiltonian function (4) into a series to make it into the form of a polynomial. Here, the Legendre expansion method is applied:

[0037]

[0038] Among them, P n is the Legendre polynomial of order \(n\). By expanding the nonlinear terms through the above transformation method, the Hamiltonian function can be transformed into the form of a polynomial of order \(n\):

[0039]

[0040] Among them, \(H\) n is a homogeneous polynomial of order \(n\), \(c\) n (\(\mu\)) is a real coefficient, and its value depends on \(\mu\) and the selected libration point \(L\) i .

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

[0042]

[0043] 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 nonlinear terms of the linear complex normal form to decouple the hyperbolic unstable direction and the central direction of the central manifold to 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.

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

[0045]

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

[0047]

[0048] Among them, J and M are defined by the following formula:

[0049]

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

[0051]

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

[0053]

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

[0055]

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

[0057]

[0058] In order 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:

[0059]

[0060] Therefore, through the transformation

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

[0062] The Hamiltonian function can be transformed into the real canonical form (9). Through the canonical equations, it can be known that the linearized dynamical equation corresponding to H2 is:

[0063]

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

[0065]

[0066] Further transform H2 into complex standard form. Since the saddle point part already has a suitable form, only the central part needs to be complex transformed:

[0067]

[0068] Through transformation, the Hamiltonian function (9) is further transformed into a linear complex standard form:

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

[0070] Now we can define the action-angle variable: I j =q j p 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.

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

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

[0073]

[0074] Among them, the calculation definition of the Poisson bracket {F, G} is

[0075]

[0076] Similarly, the second derivative of the function f can be expressed as

[0077]

[0078] Denote the n-fold Poisson bracket of the function f in the Hamiltonian system H as Then the nth derivative of f can be expressed as

[0079]

[0080] Suppose there exists a Hamiltonian system G to be constructed, which can transform the original Hamiltonian function H(q, p) of the system into Furthermore, make it have some special properties (such as eliminating some 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, and performing a Taylor series expansion at time t = 0, we can obtain The explicit transformation of:

[0081]

[0082] Among them, i = 1, 2, 3, and the Hamiltonian function of the corresponding system is transformed into

[0083]

[0084] In this Lie transformation, the Hamiltonian system G is also called the generating function. It can be proved that if we hope 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 achieve this process.

[0085] Since we seek to decouple the hyperbolic unstable direction and the central direction in the high-order terms of the Hamiltonian function, we will select 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 deal with 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

[0086]

[0087] Among them, 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) through a complex transformation, and H2 is in the linear canonical form, thus let It can be obtained that

[0088]

[0089] where h i.j is the coefficient of H3: The purpose of this application is to decouple the hyperbolic direction and the central direction. Therefore, only specific terms in need to be eliminated. 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, which can be specifically calculated by formula (28) as follows:

[0090]

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

[0092]

[0093] where, i1≠j1. Implementing the above process to the Nth order accuracy, the Hamiltonian function will be transformed into

[0094]

[0095] where R N (q,p) is the remainder term higher than the Nth order, is the term between the second order and the Nth order, and all are in the form of q1p1. For the convenience of analysis, performing the inverse transformation of the complex transformation (21) on formula (33) can convert the complex variables into real variables. Thus, through the canonical transformation of the Hamiltonian function, the decoupling of the hyperbolic unstable direction and the central direction within the Nth order accuracy is achieved.

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

[0097]

[0098] Among them, the subscript c represents the center motion mode, and s represents the saddle motion mode. For collinear libration points, their motion modes are saddle × center × center. Therefore, for the definition of the motion integrals and angular variables (I1, θ1, I2, θ2, I3, θ3), it can be defined as (I s , θ s ) to define (I1, θ1), (I c , θ c ) to define (I2, θ2) and (I3, θ3). However, when selecting the parameters for orbit characterization, it is more inclined to retain q1 and p1 rather than characterizing them with their motion integrals and angular variables, mainly for the following reasons: First, the definition of the angular variable of the saddle point involves the generation of complex variables. Although the quadrants of its complex variables contain certain physical meanings, they are rather abstract in practical applications. For details, refer to the description in Appendix

[26] ; Second, the above physical meanings can be fully characterized by q1 and p1 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 linearized dynamics of collinear libration points that the change of q1 shows an exponential increase, and the change of p1 shows an exponential decrease. 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 q1 will increase rapidly; when the target enters the vicinity of the libration point along the stable manifold, the value of p1 will approach 0 infinitely. While q2, p2, q3, p3 correspond to the motion law of the target in the central manifold, and their motion integrals and angular variables can be defined as (I c , θ c ) and have certain physical meanings: 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 physical meanings similar to phases. When θ2, θ3 have specific relationships, periodic / quasi-periodic orbits will be generated. Therefore, make the following transformation:

[0099]

[0100] where j = 2, 3. The Hamiltonian function at this time can be expressed as

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

[0102] where R N is the remainder of higher than N order, and H Nis a term between the second order and the Nth order. To sum up, the cataloging parameters of the target are selected as (q1, p1, I2, θ2, I3, θ3). According to the previously derived process, the mapping relationship between the state of the target in the rendezvous coordinate system and the cataloging parameters can be constructed:

[0103]

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

[0105]

[0106] Among them, A represents the process of coordinate translation and scaling shown in formula (3), B represents the process of constructing the Hamiltonian system and redefining the generalized coordinates and generalized momenta shown in formula (5), C represents the process of converting the linear term of the Hamiltonian function into the real canonical form shown in formula (18), D represents the complex transformation process shown in formula (21), G n represents the process of eliminating the specific terms of the nth order through the 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.

[0107] Step 108: Select the spacecraft coordinates from the parameter characterization according to the Poincaré section to calculate the characterization parameters of the reference orbit and the integral initial values of the coordinates, and calculate the mean square error between the true orbit and the reference orbit by using the characterization parameters of the reference orbit and the integral initial values of the coordinates; construct an orbit identification and solution model with the mean square error as the objective function and the preset constraint conditions; solve the orbit identification and solution model according to the Bayesian optimization algorithm to obtain the orbit identification result.

[0108] After decoupling the hyperbolic unstable direction and the central direction, analyze the two kinds of 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 laws still satisfy the canonical equations:

[0109]

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

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

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

[0113]

[0114] 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, so in the linear part {I j , H2} = 0, that is, there exists a cyclic integral I j is a constant. However, only the hyperbolic unstable direction and the center direction are separated in the high-order terms, so the high-order terms do not satisfy the integrable property. Therefore, the change of I j can be regarded as an additional secular oscillation perturbation term on the basis of the integral constant of the linear system. Similarly, for the angular variables, the linear part gives the law that the angular quantity changes linearly: θ = ωt + θ0, and the high-order terms will cause small-amplitude vibrations of the actual angular quantity on the basis of the linear change.

[0115] To verify the accuracy of the above results, a numerical method was used to calculate a south Halo orbit at the L1 point for testing, as Figure 2 shown. The following method was used for testing: The states at different positions on the Halo orbit were selected and transformed into characteristic parameters, and these were used as the benchmarks. Because the change of the parameters is fundamentally governed by the CRTBP dynamic equations, as Figure 3 shown by the blue curve in. Next, the initial state of the Halo orbit was selected and transformed into characteristic parameters. Taking this parameter as the initial value, numerical integration was carried out using (39), (41), (42), and the result is Figure 3 the red curve in. CM represents the central manifold. Because the change of the parameters is essentially governed by the canonical equations corresponding to the central manifold. It can be seen that among the characteristic parameters corresponding to the center direction, the operation results of CRTBP are consistent with the operation results of the canonical equations of the central manifold, 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 an approximately 10 -4Magnitude deviation (in an ideal periodic orbit, since there is no motion in the hyperbolic direction, q1 and p1 should remain 0), so 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 hyperbolic direction parameters will not affect the change of central direction parameters, so it has a certain robustness in the presence of a series of errors. Task analysts can focus on the evolution law of central direction characterization parameters 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.

[0116] 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 integration initial values, 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 flows of the central manifold intersect with this section. Since the Hamiltonian function (40) does not explicitly contain time t, the Hamiltonian (that is, the energy in a broad sense) is a conserved quantity during the motion. Fixing the energy level C can obtain the corresponding Poincaré section:

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

[0118] By uniformly selecting lattice points on the section to calculate the initial values corresponding to the central manifold and integrating backward, it can be found that the phase flow will cross the Poincaré section multiple times, and the crossing positions are Figure 4It has been marked that these positions will form a cross-sectional view together. The orbits of different orbit families will form cross-sectional views with different characteristics: the cross-sectional view of Lissajous orbits will pass through the entire cross-section as θ3 varies; the intersection line of the cross-section corresponding to Lyapunov orbits with the I2-θ3 plane is perpendicular to the intersection line of the cross-section corresponding to Lyapunov orbits with the I3-θ3 cross-section; quasi-Halo orbits will form a ring-shaped 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-shaped structures will appear, with their centers located at the positions of θ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-shaped 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 central manifold coordinates, they correspond to a phase difference of π in θ3. 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 views 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 center of the ring-shaped structure, that is, the Halo orbit family) is further selected to intercept the cross-section of the above three-dimensional space. Here we select the cross-section of θ3 = π / 2 and draw the situation including all energy levels, forming a Poincaré cross-sectional view that contains all orbit families near the libration point, as Figure 5 shown. It can be seen that the Lyapunov orbit family and the perpendicular Lyapunov orbit family correspond to the cases of I3 = 0 and I2 = 0 respectively. The asymptotic Lyapunov orbit (which is called Large quasi-Halo in some literature) is the boundary between Lissajous orbits and quasi-Halo orbits. When the quasi-Halo orbit satisfies specific conditions, a Halo orbit will be formed, and its position has been marked correspondingly in Figure 5 . Different energy levels C are marked by contour lines in the figure. It can be seen that at lower energy levels, there is no Halo orbit, which also reflects the origin of the Halo orbit: the Halo orbit bifurcates from the extension of the Lyapunov orbit from a lower energy level to a higher energy level. Since the cross-sectional views 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, only the phase difference between θ3 and θ2 needs to be observed for judgment.

[0119] Based on the above, orbit identification is carried out. Assume that the state of the spacecraft at the collinear libration point of the Earth-Moon system in the rendezvous coordinate system for a period of time is obtained:

[0120] (t1, t2,..., tn ) → (X1, X2, ..., X n ) (44)

[0121] Among them First, calculate the characterization parameter σ corresponding to each moment t according to the characterization parameter conversion method described in Equation (38). i The corresponding characterization parameter σ i = (q1 (i) , p1 (i) , I2 (i) , θ2 (i) , I3 (i) , θ3 (i) ) T , and then select coordinates from the Poincaré section diagram of Figure 5 to calculate the characterization parameter of the reference orbit. According to the previous description, there is a correspondence between the coordinates (I2 (0) , I3 (0) ) in the Poincaré section and the reference orbit. Specifically, the initial integration values can be obtained according to the coordinates Then, perform orbit integration according to Equations (39), (41), and (42), and the result can be expressed by the phase flow function: represents the result of numerical integration of the initial value σ0 from t0 to t i . Since the time of the spacecraft orbit is not aligned with the time of the reference orbit, t0 needs to be optimized. For a periodic orbit, t0 only needs to be found within one period; for a quasi-periodic orbit, t0 may need to span a longer time period to roughly cover the motion of the spacecraft on the entire quasi-periodic orbit. The method of quantifying the real orbit and the reference orbit can be represented by the mean square error:

[0122]

[0123] From Equation (45), we can see that only the errors of the two action variables I2 and I3 are used as indicators when defining the mean square error. This is because through the transformation of action-angle variables, the characteristics of the central manifold are concentrated on the integration of the linear part, and the angle variables, which are periodic variables, can be ignored. This also reflects the purpose of only selecting some action-angle variables when defining the characterization parameter. MSE describes the deviation between the characterization parameters of the real orbit and the reference orbit. When the deviation reaches the minimum value, it is considered that the orbit corresponding to (I2 (0) , I3 (0) ) in the section diagram can be regarded as the reference orbit corresponding to the real orbit. Therefore, the problem of orbit identification can be expressed as the solution of the following optimization model:

[0124]

[0125] Among them, the first constraint relationship represents the search range of coordinates on the Poincaré section diagram. This is because the change rule of I j is to make small fluctuations above and below a reference value. Therefore, the characteristic parameters of the true orbit can be arbitrarily selected as the center point and searched within the range of I max . The second constraint relationship represents the search range at the initial moment, which is mainly used to align the time of the true orbit and the reference orbit. According to tests, taking I max = 0.1 and T max = 50 can achieve a good orbit recognition effect while ensuring the calculation efficiency.

[0126] Since a long-time orbit integration is required during the optimization process, in order to further improve the operation efficiency, Bayesian optimization is selected as the optimization method to solve the optimization problem described in Equation (46). Bayesian optimization is a global optimization strategy based on Bayes' theorem, mainly used to find the optimal solution of complex functions, especially suitable for the cases where the calculation cost of the objective function is high, it is difficult to directly take the derivative or there is no explicit expression. Its core idea is to construct a prior probability model of the objective function, and then continuously update the posterior probability distribution according to the existing observed data, so as to predict the value of the objective function in the unexplored area, and then guide the search process to efficiently approach the global optimal solution and obtain the orbit recognition result.

[0127] The following shows the process of orbit recognition through an example. Assume that the spacecraft is moving near the L2 point, and the state information of the spacecraft in the rendezvous coordinate system within a period of time is returned according to the inter-satellite positioning, as Figure 7 shown. First, the state sequence and time are made dimensionless; then the dimensionless state variables are converted into characteristic parameters, and the results are as Figure 8 shown; the action (I2, I3) at time 0 is selected as the center point, and the Bayesian optimization method is used to solve the optimization model (46). The optimization process is as Figure 9 shown. It can be seen that the Bayesian method can find the global optimal point within 30 function operations, and the efficiency is very high. After optimization, the position of the reference orbit corresponding to the spacecraft trajectory in the Poincaré section diagram can be obtained, as Figure 10As shown, its orbit in the rendezvous coordinate system and the corresponding characterization parameters have been marked in the previous figures. It can be seen that affected by various perturbations, the true orbit of the spacecraft will not completely coincide with the reference orbit, but move near the reference orbit, which also makes the characterization parameters corresponding to the true orbit not strictly coincide with those of the reference orbit. However, we find that when the change in the action function of the two orbits is sufficiently close, that is, when the mean square error described by Equation (45) reaches the minimum value, the reference orbit is already very close to the true orbit. At this time, the orbital parameters of the true orbit, such as the period, amplitude, etc., can be characterized by the parameters of its reference orbit. For the target cataloging and storage management, six characterization parameters can be used, and there is a one-to-one correspondence between the characterization parameters and the target state variables. By the position of the reference orbit on the Poincaré section diagram and in combination with the phase of θ3, we can also easily determine that the spacecraft is on the north family Halo orbit at the L2 point. For the characteristics of the target orbit, such as the amplitude, period, energy, etc., they can also be judged by its reference orbit. Since the motion of the target near L2 can be represented by a coordinate in the Poincaré section diagram shown in Figure 10 , therefore Figure 5 's Poincaré section diagram can be used as a diagram of the distribution situation of the libration point targets, for mission analysts to quickly judge the target distribution near the L2 point.

[0128] The above-mentioned orbit recognition method based on the parameter characterization of the Earth-Moon collinear libration point orbit obtains the motion equation of the spacecraft in the rendezvous coordinate system, describes the circular restricted three-body problem as a chaotic Hamiltonian dynamical system, and then transforms it to the new coordinate system to obtain the Hamiltonian function of the corresponding system, so that the originally complex and difficult-to-directly-analyze three-body problem is transformed into a form that can be deeply studied using Hamiltonian theory. Through this transformation, the motion law of the spacecraft can be systematically understood from the perspectives of physical quantities such as energy and momentum, laying a foundation for accurately identifying the orbit subsequently. Expand the non-linear term in the Hamiltonian function by Legendre series to obtain the Hamiltonian function in polynomial form, and use its quadratic term as the linearized model of the libration point. Through the linearized model, the motion of the spacecraft near the libration point can be initially approximately analyzed. The linearized model can quickly and relatively accurately describe its motion, providing an initial and relatively accurate reference for the overall orbit recognition. Then, according to the real linear symplectic transformation matrix, transform the linearized model into a real canonical form, and then perform a complex transformation on the central part of the real canonical form to obtain a linear complex canonical form. The canonical form transformation can transform a complex linear system into a more canonical and easy-to-analyze form. Through this transformation, different motion modes and characteristics in the system can be more clearly separated, which is conducive to accurately judging the type of the spacecraft orbit and improving the recognition accuracy. Perform a canonical transformation on the non-linear term 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. The decoupling operation enables the separate treatment of motions with different properties in the system. In the three-body problem, the motion on the center manifold contains the main motion characteristics of the spacecraft near the libration point, while the motion in the hyperbolic unstable direction may perturb the orbit. Through decoupling, these two motions with different properties can be accurately analyzed and described respectively. During orbit recognition, the stable motion part of the spacecraft on the center manifold and the part affected by the hyperbolic unstable direction can be more accurately grasped, thereby improving the recognition accuracy of the overall orbit. 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 center manifold, and construct the mapping relationship between the CRTBP coordinates and the characterization parameters. The local action-angle variables can more intuitively reflect the motion characteristics of the spacecraft near the libration point. By constructing the mapping relationship, the complex CRTBP coordinates are transformed into more physically meaningful characterization parameters, which can more accurately depict the characteristics of the orbit. Select the spacecraft coordinates from the parameter characterization according to the Poincaré section to calculate the characterization parameters and the integral initial values of the coordinates of the reference orbit, and use these values to calculate the mean square error between the real orbit and the reference orbit. Construct an orbit recognition solution model with the mean square error as the objective function and the preset constraint conditions, and then solve the model according to the Bayesian optimization algorithm. The Poincaré section can transform the continuous orbit motion into a discrete point set, facilitating the selection of representative spacecraft coordinates from it.By calculating the mean square error, the difference between the true orbit and the reference orbit can be quantified. The smaller the mean square error, the closer the true orbit is to the reference orbit. Taking the mean square error as the objective function and combining the constraint conditions to construct a solution model, the Bayesian optimization algorithm can quickly and accurately search for the orbital parameters that best match the true orbit among numerous possible orbits, thus obtaining an accurate orbit identification result. This method based on the optimization algorithm and the quantified error can efficiently find the optimal solution in the complex orbital space of the three-body problem, greatly improving the accuracy of orbit identification.

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

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

[0131]

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

[0133]

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

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

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

[0137]

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

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

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

[0141] H2 = λq1p1 + iω p q2p2 + iω v q3p3。

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

[0143] Performing 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 center manifold, the final Hamiltonian function is

[0144]

[0145] where, R N (q, p) is a remainder term of order higher than N, is a term between the second order and the Nth order, and both are in the form of q1p1 R N .

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

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

[0148] where, H2 represents the linear complex normal form, R N (q, p) is a remainder term of order higher than N, is a term between the second order and the Nth order, and both are in the form of q1p1 R N .

[0149] In one embodiment, constructing the mapping relationship between the CRTBP coordinates and the characterization parameters includes:

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

[0151]

[0152] where, (X, Y, Z) represents the position of the spacecraft in the rendezvous coordinate system, denotes the velocity of the spacecraft in the rendezvous coordinate system. q1 and p1 respectively characterize the degree to which the spacecraft cuts 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, and 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, periodic / quasi-periodic orbits will be generated. The physical meaning of phase includes the variation laws of generalized coordinates and generalized momenta.

[0153] In one embodiment, the characterization parameters of the reference orbit and the integral initial values of the coordinates are selected from the parameter characterization according to the Poincaré section, including:

[0154] According to the Poincaré section, the spacecraft coordinates are selected from the parameter characterization to calculate the characterization parameters corresponding to each moment t i of the reference orbit as σ i =(q1 (i) , p1 (i) , I2 (i) , θ2 (i) , I3 (i) , θ3 (i) ) T and the integral initial values of the coordinates are where i represents the i-th moment and T represents the transpose operation.

[0155] In one embodiment, the mean square error between the true orbit and the reference orbit is calculated using the characterization parameters of the reference orbit and the integral initial values of the coordinates, including:

[0156] The mean square error between the true orbit and the reference orbit calculated using the characterization parameters of the reference orbit and the integral initial values of the coordinates is:

[0157]

[0158] where denotes the phase flow function after orbit integration of the integral initial values at moment t i , representing the result of numerical integration of the initial value σ0 from t0 to t i . t0 represents the initial moment and n represents the number of data at n moments.

[0159] In one embodiment, an orbit identification and solution model is constructed with the mean square error as the objective function and pre-set constraint conditions, including:

[0160] The orbit identification and solution model constructed with the mean square error as the objective function and pre-set constraint conditions is:

[0161]

[0162] Among them, the first constraint relationship represents the search range of coordinates on the Poincaré section diagram, and any characteristic parameter of the true orbit is arbitrarily selected as the center point Search within the range of I max The range of I max represents the range of coordinate search in the situation diagram, and T max represents the time range for adjustment to align the data. The second constraint relationship represents the search range at the initial moment, which is mainly used to align the time of the true orbit and the reference orbit.

[0163] It should be understood that although Figure 1 The steps in the flowchart of are shown sequentially according to the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless there is a clear description in this article, there is no strict order limit 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 may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily executed at the same moment, but can be executed at different moments. The execution order of these sub-steps or stages is not necessarily sequential, 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. Their descriptions are relatively specific and detailed, but they should not be construed as limiting 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. An orbit identification method based on the characterization of orbital parameters of the Earth-Moon collinear libration point, 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 spacecraft coordinates are selected from the parameter characterization according to the Poincare section to calculate the characterization parameters of the reference orbit and the initial value of the integral of the coordinates; the mean square error between the real orbit and the reference orbit is calculated using the characterization parameters of the reference orbit and the initial value of the integral of the coordinates; an orbit identification solution model is constructed with the mean square error as the objective function and pre-set constraints; the orbit identification solution model is solved according to the Bayesian optimization algorithm to obtain an orbit identification result.

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 7, characterized in that According to the Poincare section, the spacecraft coordinates are selected from the parameter characterization to calculate the characterization parameters of the reference orbit and the initial value of the integral of the coordinates, including: According to the Poincare section, the spacecraft coordinates are selected from the parameter representation to calculate the reference orbit at each time t i The corresponding characterization parameter is σ i =(q1 (i) ,p1 (i) ,I2 (i) ,θ2 (i) ,I3 (i) ,θ3 (i) ) T The initial value of the integral of the coordinates is Among them, i represents the i-th moment, and T represents the transposition operation.

9. The method according to claim 8, characterized in that Calculating the mean square error between the real orbit and the reference orbit using the characterization parameters of the reference orbit and the initial value of the integration of the coordinates, including: The mean square error between the real orbit and the reference orbit is calculated using the characterization parameters of the reference orbit and the initial value of the integral of the coordinates: in, Indicates t i The phase flow function after orbital integration of the initial value of the integral at time t represents the initial value σ0 from t0 to t i The result of numerical integration, t0 represents the initial time, and n represents the data of a total of n time points.

10. The method according to claim 9, characterized in that A track identification solution model is constructed with the mean square error as the objective function and pre-set constraints, including: The track identification solution model is constructed with the mean square error as the objective function and the preset constraints as follows: The first constraint relationship represents the search range of the coordinates on the Poincare section diagram, and the characterization parameters of the real orbit are arbitrarily selected as the center point in I max Search within the range, I max Indicates the range of coordinate search in the situation map, T max It represents the time range for adjusting the data to align. The second constraint relationship represents the search range of the initial time, which is mainly used to align the time of the real orbit with the reference orbit.

Images