A method and device for calculating an unloading attitude of a momentum wheel of a halo orbit satellite

Abstract

The application discloses a kind of Halo orbit satellite momentum wheel unloading posture calculation method and device, comprising: the velocity increment generated in historical unloading process is estimated by precise orbit determination, and the relationship between its and historical unloading moment of angular momentum is fitted to obtain fitting coefficient.Nominal control parameter required for satellite to maintain Halo orbit configuration and nominal position and velocity after first crossing XZ plane of geolunar rotating coordinate system are determined.Velocity increment generated by unloading is determined based on unloading moment of angular momentum and fitting coefficient, and unloading position and velocity after first crossing XZ plane after unloading are determined in combination with unloading posture.Bias between nominal orbit after control and unloading orbit after unloading is obtained according to nominal position and velocity and unloading position and velocity after unloading, and based on coordinate wheel method, unloading posture is optimized to make bias minimum, and unloading optimal posture quaternion is obtained, to realize that velocity increment generated by unloading is used for satellite Halo orbit configuration maintenance, offset Halo orbit divergence tendency, reduce maintenance frequency.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of space dynamics and control, and particularly relates to a Halo orbit satellite momentum wheel unloading attitude calculation method and device. BACKGROUND

[0002] In the field of space dynamics, Halo orbits in the Earth-Moon space have wide applications. Satellites located in Halo orbits need to frequently implement orbit control to maintain Halo orbit configurations in long-term operation due to the instability of Halo orbits and long-term interference of solar radiation pressure, gravity of large celestial bodies and the like, which will greatly consume the remaining propellant of the satellites.

[0003] Modern satellites usually use momentum wheels as actuators of attitude control systems, and the momentum wheels need to be unloaded to restore to the initial rotation speed when the rotation speed reaches a critical value. Satellites located in Halo orbits are long-term subjected to interference torques generated by solar radiation pressure, and the satellites absorb the interference torques through momentum wheels to maintain target attitudes, and the momentum wheels need to consume propellant for unloading after saturation, and the speed increment generated by unloading will further disturb the orbit.

[0004] Therefore, how to utilize the inevitable momentum wheel unloading of satellites to reduce the orbit maintenance frequency and reduce the propellant consumption is a problem to be solved at present. SUMMARY

[0005] Embodiments of the present application provide a Halo orbit satellite momentum wheel unloading attitude calculation method, which is used for actively unloading the momentum wheels at appropriate attitudes, uses the speed increment generated by unloading for satellite orbit configuration maintenance, realizes the use of momentum wheel unloading to offset the divergence trend of Halo orbits, reduces the orbit maintenance frequency, and reduces the propellant consumption of the satellites.

[0006] In a first aspect, embodiments of the present application provide a Halo orbit satellite momentum wheel unloading attitude calculation method, comprising:

[0007] planning a satellite unloading time, determining a nominal control direction required for the satellite to maintain a halo orbit configuration at the unloading time and a nominal controlled orbit, and determining a nominal position and velocity in an inertial system when the satellite first crosses an XZ plane of an Earth-Moon rotating coordinate system on the nominal controlled orbit based on the nominal controlled orbit;

[0008] establishing an angular momentum coordinate system of the satellite based on an angular momentum vector in the inertial system at the unloading time and the nominal control direction, and determining an angular optimization variable of the satellite, the angular optimization variable being Euler angles of rotation of the satellite around each coordinate axis of the angular momentum coordinate system according to a rotation sequence;

[0009] determine a velocity increment generated by the satellite during the unloading process based on the angular momentum vector and the fitting coefficient, and determine a post-unloading orbit of the satellite, and a post-unloading position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the Earth-Moon rotating coordinate system on the post-unloading orbit; the fitting coefficient represents a fitting relationship between the angular momentum of the satellite at the unloading moment and the velocity increment generated during the unloading process;

[0010] obtain a deviation between the post-control nominal orbit and the post-unloading orbit according to the post-unloading position and velocity and the post-control nominal position and velocity, and optimize the deviation according to the angle optimization variable based on the coordinate permutation method to obtain an optimal attitude quaternion of the satellite momentum wheel unloading, so that the satellite performs momentum wheel unloading according to the optimal attitude quaternion, and the goal of maintaining the halo orbit configuration through active unloading is achieved.

[0011] In the above technical solution, first, the unloading moment of the satellite is planned according to the rotation speed of the satellite momentum wheel. The halo orbit can also be referred to as a Halo orbit. Then, the nominal control direction and the post-control nominal orbit of the satellite are determined according to the halo orbit maintenance control method of the satellite. Then, the post-control nominal position and velocity of the satellite when the satellite first crosses the XZ plane of the Earth-Moon rotating coordinate system are determined according to the position of the satellite on the post-control nominal orbit. The post-control nominal position and velocity include the post-control nominal position and post-control nominal velocity of the satellite when the satellite crosses the XZ plane. The satellite momentum wheel generates an angular momentum vector when unloading, so the velocity increment generated by the satellite during the unloading process can be determined based on the angular momentum vector and the fitting coefficient, and the post-unloading orbit of the satellite can be determined based on the velocity increment of the satellite. The fitting coefficient represents a fitting relationship between the angular momentum of the satellite at the unloading moment and the velocity increment generated during the unloading process, and the product of the angular momentum and the fitting coefficient is the velocity increment. Then, the post-unloading position and velocity of the satellite when the satellite first crosses the XZ plane of the Earth-Moon rotating coordinate system are determined according to the position of the satellite on the post-unloading orbit. The post-unloading position and velocity include the post-unloading position and post-unloading velocity of the satellite when the satellite crosses the XZ plane. Finally, the deviation between the post-control nominal orbit and the post-unloading orbit is obtained according to the deviation between the post-unloading position and velocity and the post-control nominal position and velocity. Finally, the deviation between the post-control nominal orbit and the post-unloading orbit is optimized by adjusting the attitude of the satellite at the unloading moment to obtain the optimal attitude quaternion of the satellite momentum wheel unloading, so that the velocity increment generated by the unloading is used for satellite orbit configuration maintenance, the orbit maintenance frequency is reduced, and the consumption of propellant of the satellite is reduced.

[0012] Optionally, before planning the unloading moment of the satellite, the method further comprises:

[0013] acquire a historical unloading moment attitude quaternion, a historical unloading moment angular momentum and a corresponding historical unloading speed increment; the historical unloading speed increment is a speed increment generated by the satellite in a historical unloading process in an inertial system, the historical unloading speed increment is obtained according to precise orbit determination results before and after unloading of the satellite, and the historical unloading moment angular momentum is three-axis angular momentum of the satellite in a body system at the unloading moment;

[0014] determine a first conversion matrix converted from the inertial system to the body system at the unloading moment based on the historical unloading moment attitude quaternion;

[0015] convert the historical unloading speed increment from the inertial system to the body system at the unloading moment based on the first conversion matrix, and obtain a fitting coefficient according to a fitting relationship between the historical unloading speed increment and the historical unloading moment angular momentum in the body system at the unloading moment.

[0016] In the technical solution, it can be understood that the satellite is located at different positions at different historical unloading moments, and therefore the body system at the unloading moment of the satellite is also different.

[0017] Optionally, the nominal position and speed of the satellite in the inertial system when the satellite first crosses the XZ plane of the geolunar rotating coordinate system on the nominal controlled orbit are determined based on the nominal controlled orbit, and the nominal position and speed of the satellite in the inertial system when the satellite first crosses the XZ plane of the geolunar rotating coordinate system on the nominal controlled orbit include:

[0018] The position and speed of the satellite in the inertial system are dynamically integrated based on the nominal controlled orbit according to a specified step length, to obtain position and speed of the satellite at the kth moment, and a rotation matrix converted from the inertial system to the geolunar rotating coordinate system at the kth moment is calculated; k is an integer greater than or equal to 1, and the first moment is located after the unloading moment and is separated from the unloading moment by a preset time period;

[0019] The position and speed of the satellite at the kth moment are converted from the inertial system to the geolunar rotating coordinate system at the kth moment based on the rotation matrix, to obtain position coordinates of the satellite at the kth moment in the geolunar rotating coordinate system at the kth moment;

[0020] The product of the Y-axis component of the position coordinates of the satellite at the kth moment and the Y-axis component of the position coordinates of the satellite at the k+1th moment is calculated;

[0021] When the product is less than 0 for the first time, the satellite first crosses the XZ plane of the geolunar rotating coordinate system on the nominal controlled orbit, and the kth moment and the k+1th moment are accurately iterated to obtain the nominal position and speed of the satellite in the inertial system when the satellite first crosses the XZ plane of the geolunar rotating coordinate system on the nominal controlled orbit.

[0022] Optionally, an angular momentum coordinate system of the satellite is established based on the angular momentum vector in the inertial system at the unloading moment and the nominal control direction, and an angle optimization variable of the satellite is determined, including:

[0023] A direction of the angular momentum vector is taken as a Y axis of the angular momentum coordinate system, a direction of a cross product result of the Y axis and the nominal control direction is taken as an X axis, a Z axis is determined based on a right-hand rule, the angular momentum coordinate system is established, and a second conversion matrix converted from the inertial system to the angular momentum coordinate system is determined;

[0024] The angle optimization variable is determined based on angles of rotation of the satellite around each coordinate axis of the angular momentum coordinate system, and a third conversion matrix converted from the angular momentum coordinate system to the body system at the unloading moment is determined;

[0025] A fourth conversion matrix converted from the inertial system to the body system at the unloading moment is obtained based on the second conversion matrix and the third conversion matrix.

[0026] Optionally, a velocity increment generated by the satellite in an unloading process is determined based on the angular momentum vector and a fitting coefficient, and a post-unloading orbit and a post-unloading position and velocity of the satellite in the inertial system at a first time when the satellite crosses an XZ plane of a geolunar rotation coordinate system are determined, including:

[0027] The angular momentum vector is converted from the inertial system to the body system at the unloading moment based on the fourth conversion matrix, and a velocity increment generated by the satellite in the unloading process in the body system at the unloading moment is obtained based on the angular momentum vector in the body system at the unloading moment and the fitting coefficient;

[0028] The velocity increment generated in the unloading process is converted from the body system at the unloading moment to the inertial system based on an inverse matrix of the fourth conversion matrix, and the post-unloading orbit is determined based on the velocity increment generated in the unloading process;

[0029] The position and velocity of the satellite in the inertial system are dynamically integrated based on the post-unloading orbit according to a specified step length, a position and velocity of the satellite at an mth moment are obtained, and a rotation matrix converted from the inertial system to a geolunar rotation coordinate system at the mth moment is calculated; m is an integer greater than or equal to 1, and the first moment is located after the unloading moment and is separated from the unloading moment by a preset time period;

[0030] The position and velocity at the mth moment are converted from the inertial system to the geolunar rotation coordinate system at the mth moment based on the rotation matrix, and a position coordinate at the mth moment in the geolunar rotation coordinate system at the mth moment is obtained;

[0031] A product result of a Y axis component of the position coordinate at the mth moment and a Y axis component of a position coordinate at an (m+1) th moment is calculated;

[0032] When the product result is first less than 0, the satellite first crosses the XZ plane of the geolunar rotation coordinate system on the post-unloading orbit, and the m-th moment and the m+1-th moment when the product result is first less than 0 are accurately iterated to obtain the post-unloading position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the geolunar rotation coordinate system on the post-unloading orbit.

[0033] Optionally, the deviation between the nominal post-control orbit and the post-unloading orbit is obtained according to the nominal position and velocity and the post-unloading position and velocity, and includes:

[0034] A product result of a modulus value of a vector difference between a satellite position in the post-unloading position and velocity and a satellite position in the nominal position and velocity and a position weight is calculated as a position deviation;

[0035] A product result of a modulus value of a vector difference between a satellite velocity in the post-unloading position and velocity and a satellite velocity in the nominal position and velocity and a velocity weight is calculated as a velocity deviation; the position weight is less than the velocity weight, and the position weight and the velocity weight are 1;

[0036] A summation result of the position deviation and the velocity deviation is calculated as the deviation between the nominal post-control orbit and the post-unloading orbit.

[0037] Optionally, the deviation is optimized according to the angle optimization variable based on a coordinate permutation method to obtain an optimal attitude quaternion of the satellite momentum wheel unloading, and includes:

[0038] The X axis of the angular momentum coordinate system is divided into a positive direction and a negative direction to obtain two half spaces, and a value range of the angle optimization variable in each half space is determined;

[0039] For any half space in the two half spaces, the deviation is optimized according to the angle optimization variable based on the coordinate permutation method according to the value range of the angle optimization variable in the half space until the deviation is minimum to obtain an optimal value of the angle optimization variable corresponding to the half space and a minimum deviation;

[0040] The optimal value of the angle optimization variable corresponding to the smaller minimum deviation of the two half spaces is taken as an optimal value of the angle optimization variable of the satellite.

[0041] The optimal attitude quaternion of the satellite momentum wheel unloading is obtained based on the fourth conversion matrix corresponding to the optimal value of the angle optimization variable of the satellite.

[0042] In a second aspect, an embodiment of the present application provides a Halo orbit satellite momentum wheel unloading attitude calculation device, including:

[0043] an acquisition module configured to plan a satellite unloading time, determine a nominal control direction and a nominal post-control orbit required for the satellite to maintain a halo orbit configuration at the unloading time, and determine a nominal position and velocity of the satellite in an inertial system when the satellite first crosses an XZ plane of a geocentric inertial coordinate system on the nominal post-control orbit based on the nominal post-control orbit;

[0044] a processing module configured to establish an angular momentum coordinate system of the satellite based on an angular momentum vector in the inertial system at the unloading time and the nominal control direction, and determine an angular optimization variable of the satellite, the angular optimization variable being Euler angles of the satellite rotating around each coordinate axis of the angular momentum coordinate system in a rotation sequence;

[0045] determine a velocity increment generated by the satellite during the unloading process based on the angular momentum vector and a fitting coefficient, and determine a post-unloading orbit and a post-unloading position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the geocentric inertial coordinate system on the post-unloading orbit, the fitting coefficient representing a fitting relationship between the angular momentum of the satellite at the unloading time and the velocity increment generated by the satellite during the unloading process;

[0046] obtain a deviation between the nominal post-control orbit and the post-unloading orbit according to the nominal position and velocity and the post-unloading position and velocity, and optimize the deviation according to the angular optimization variable based on a coordinate permutation method to obtain an optimal attitude quaternion of the satellite momentum wheel unloading, so that the satellite performs the momentum wheel unloading according to the optimal attitude quaternion to achieve the goal of maintaining the halo orbit configuration through active unloading.

[0047] Optionally, the processing module is further configured to:

[0048] acquire a historical unloading time attitude quaternion, a historical unloading time angular momentum, and a corresponding historical unloading velocity increment, the historical unloading velocity increment being a velocity increment generated by the satellite in the inertial system during a historical unloading process, the historical unloading velocity increment being obtained according to precise orbit determination results before and after the satellite unloading, and the historical unloading time angular momentum being three-axis angular momentum of the satellite in the body system at the unloading time;

[0049] determine a first conversion matrix from the inertial system to the body system at the unloading time based on the historical unloading time attitude quaternion;

[0050] convert the historical unloading velocity increment from the inertial system to the body system at the unloading time based on the first conversion matrix, and obtain a fitting coefficient according to a fitting relationship between the historical unloading velocity increment and the historical unloading time angular momentum in the body system at the unloading time.

[0051] Optionally, the processing module is specifically configured to:

[0052] performing dynamic integration on the position and velocity of the satellite in the inertial system based on the nominal post-control orbit according to a specified step size, obtaining position and velocity of the satellite at the kth moment, and calculating a rotation matrix converted from the inertial system to a geocentric rotating coordinate system at the kth moment; k is an integer greater than or equal to 1, and the first moment is located after the unloading moment and is separated from the unloading moment by a preset time period;

[0053] converting the position and velocity at the kth moment from the inertial system to the geocentric rotating coordinate system at the kth moment based on the rotation matrix, and obtaining position coordinates of the satellite at the kth moment in the geocentric rotating coordinate system at the kth moment;

[0054] calculating a product result of a Y-axis component of the position coordinates at the kth moment and a Y-axis component of position coordinates at the k+1th moment;

[0055] when the product result is less than 0 for the first time, the satellite crosses the XZ plane of the geocentric rotating coordinate system for the first time on the nominal post-control orbit, performing accurate iteration on the kth moment and the k+1th moment, and obtaining a nominal position and velocity of the satellite in the inertial system when the satellite crosses the XZ plane of the geocentric rotating coordinate system for the first time on the nominal post-control orbit.

[0056] Optionally, the processing module is specifically configured to:

[0057] taking the direction of the angular momentum vector as a Y-axis of the angular momentum coordinate system, taking the direction of the cross product result of the Y-axis and the nominal control direction as an X-axis, determining a Z-axis based on the right-hand rule, establishing the angular momentum coordinate system, and determining a second conversion matrix converted from the inertial system to the angular momentum coordinate system;

[0058] determining the angle optimization variable based on the angle of rotation of the satellite around each coordinate axis of the angular momentum coordinate system, and determining a third conversion matrix from the angular momentum coordinate system to the body system at the unloading moment;

[0059] based on the second conversion matrix and the third conversion matrix, obtaining a fourth conversion matrix from the inertial system to the body system at the unloading moment.

[0060] Optionally, the processing module is specifically configured to:

[0061] based on the fourth conversion matrix, converting the angular momentum vector from the inertial system to the body system at the unloading moment, and based on the angular momentum vector in the body system at the unloading moment and the fitting coefficient, obtaining a velocity increment generated by the satellite in the unloading process in the body system at the unloading moment;

[0062] based on an inverse matrix of the fourth conversion matrix, converting the velocity increment generated in the unloading process from the body system at the unloading moment to the inertial system, and based on the velocity increment generated in the unloading process, determining a post-unloading orbit.

[0063] performing dynamic integration on the position and velocity of the satellite in the inertial system based on the unloaded orbit according to a specified step size, to obtain the position and velocity of the satellite at the mth moment, and calculating a rotation matrix converted from the inertial system to the geocentric rotating coordinate system at the mth moment; m is an integer greater than or equal to 1, and the 1st moment is located after the unloading moment and is separated from the unloading moment by a preset time period;

[0064] converting the position and velocity at the mth moment from the inertial system to the geocentric rotating coordinate system at the mth moment based on the rotation matrix, to obtain the position coordinates at the mth moment in the geocentric rotating coordinate system at the mth moment;

[0065] calculating the product of the Y-axis component of the position coordinates at the mth moment and the Y-axis component of the position coordinates at the m+1th moment;

[0066] when the product result is less than 0 for the first time, the satellite crosses the XZ plane of the geocentric rotating coordinate system for the first time on the unloaded orbit, and the mth moment and the m+1th moment when the product result is less than 0 for the first time are accurately iterated to obtain the unloaded position and velocity in the inertial system when the satellite crosses the XZ plane of the geocentric rotating coordinate system for the first time on the unloaded orbit.

[0067] Optionally, the processing module is specifically configured to:

[0068] calculate the product of the modulus of the vector difference between the position of the satellite in the unloaded position and velocity and the position of the satellite in the nominal position and velocity and the position weight as a position deviation;

[0069] calculate the product of the modulus of the vector difference between the velocity of the satellite in the unloaded position and velocity and the velocity of the satellite in the nominal position and velocity and the velocity weight as a velocity deviation; the position weight is less than the velocity weight, and the sum of the position weight and the velocity weight is 1;

[0070] calculate the sum of the position deviation and the velocity deviation as the deviation between the nominal controlled orbit and the unloaded orbit.

[0071] Optionally, the processing module is specifically configured to:

[0072] divide the X-axis positive direction and the X-axis negative direction of the angular momentum coordinate system to obtain two half spaces, and determine the value range of the angular optimization variable in each half space;

[0073] For any one of the two half spaces, based on coordinate transformation method, according to the value range of the angle optimization variable in the half space, the deviation is optimized according to the angle optimization variable until the deviation is minimum, and the optimal value of the angle optimization variable corresponding to the half space and the minimum deviation are obtained;

[0074] Among the minimum deviations corresponding to the two half spaces, the optimal value of the angle optimization variable corresponding to the smaller minimum deviation is taken as the optimal value of the angle optimization variable of the satellite;

[0075] Based on the fourth conversion matrix corresponding to the optimal value of the angle optimization variable of the satellite, the optimal attitude quaternion of the momentum wheel unloading of the satellite is obtained.

[0076] In a third aspect, an embodiment of the present application further provides a computer device, comprising:

[0077] a memory for storing program instructions;

[0078] a processor for calling the program instructions stored in the memory, and performing the above-mentioned calculation method of the Halo orbit satellite momentum wheel unloading attitude according to the obtained program.

[0079] In a fourth aspect, an embodiment of the present application further provides a computer readable storage medium, which stores computer executable instructions, and the computer executable instructions are used for making a computer execute the above-mentioned calculation method of the Halo orbit satellite momentum wheel unloading attitude.

[0080] In a fifth aspect, an embodiment of the present application further provides a computer program product, which comprises an executable program, and the executable program is used for executing the above-mentioned calculation method of the Halo orbit satellite momentum wheel unloading attitude by a processor. BRIEF DESCRIPTION OF DRAWINGS

[0081] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed in the embodiment description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0082] Figure 1 A system architecture schematic diagram is provided for the embodiments of the present application;

[0083] Figure 2 A flowchart of a Halo orbit satellite momentum wheel unloading attitude calculation method is provided for the embodiments of the present application;

[0084] Figure 3A structural schematic diagram of a Halo orbit satellite momentum wheel unloading attitude calculation method provided by an embodiment of the present application. DETAILED DESCRIPTION

[0085] In order to make the objects, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with reference to the drawings. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments of the present application, all other embodiments obtained by those skilled in the art without creative work fall within the scope of protection of the present application.

[0086] The application scenarios described in the embodiments of the present application are used to more clearly illustrate the technical solutions protected by the embodiments of the present application, and do not constitute a limitation on the technical solutions provided by the embodiments of the present application. Those skilled in the art can know that, with the emergence of new application scenarios, the technical solutions provided by the embodiments of the present application are also applicable to similar technical problems. The terms "first" and "second" in the specification and claims of the present application and the above-mentioned drawings are used to distinguish different objects, rather than to describe a specific order. In the description of the present application, unless otherwise specified, the meaning of "multiple" is two or more.

[0087] Before introducing the Halo orbit satellite momentum wheel unloading attitude calculation method provided by the embodiments of the present application, in order to facilitate understanding, first, the following terms related to the embodiments of the present application are introduced.

[0088] Halo orbit (Halo Orbit): also known as halo orbit, halo orbit, is an important concept in the three-body problem of orbit mechanics. It is a periodic three-dimensional orbit around one of the Lagrange points L1, L2 and L3. The Lagrange point is an imaginary point in space, but it has the unique feature of being surrounded by a spacecraft. Halo orbits exist in any three-body system, such as the sun-earth satellite system and the earth-moon satellite system. In the vicinity of each Lagrange point, there are continuous southward or northward halo orbit curve clusters.

[0089] Orbital elements: are the basic parameters required to describe the motion of a spacecraft in orbit, which together determine the motion trajectory and position of the spacecraft. The orbital elements include the following specific parameters:

[0090] 1. Semi-major axis (a): half of the major axis of an elliptical orbit, used to describe the size of the orbit. The larger the semi-major axis, the higher the energy of the orbit.

[0091] 2. Eccentricity (e): represents the shape of the orbit, which is a quantitative indicator of the ellipticity of the orbit. When e = 0, it is a circular orbit, when 0 < e < 1, it is an elliptical orbit, and when e >= 1, it is a parabolic or hyperbolic orbit.

[0092] 3. Inclination (i): The angle between the orbital plane and the reference plane (usually the Earth's equatorial plane), which determines the degree of tilt of the orbit.

[0093] 4. Right ascension of the ascending node (Ω): The angle of the ascending node (the point where the orbit crosses the equator from south to north) relative to the reference direction (usually the vernal equinox).

[0094] 5. Argument of perigee (ω): The angle from the ascending node to the orbit's perigee, describing the orientation of the orbit's ellipse.

[0095] 6. True anomaly The instantaneous angle measured counterclockwise from the perigee to the satellite's position, used to determine the specific position of the spacecraft on the orbit.

[0096] Normal: In the literal sense, the direction of the normal is perpendicular to the tangent, i.e., the direction of the tangent. Both the tangent and the normal are relative to interfaces, trajectories, etc.

[0097] Lunar orbit: An astronomical term referring to the projection of the moon's orbit onto the celestial sphere. The intersection angle between the lunar orbit and the ecliptic varies between 4°57' and 5°19', with an average value of about 5°09', and a cycle of about 173 days. Due to the sun's gravitational influence on the moon, the line connecting the two nodes moves westward along the ecliptic in the opposite direction of the moon's motion, a phenomenon known as node regression.

[0098] Satellite attitude quaternion: A parameter describing the orientation of the satellite's body frame relative to the reference frame, used to determine the satellite's attitude.

[0099] Figure 1 An exemplary system architecture suitable for embodiments of the present application is shown, which includes a server 100, which can include a processor 110, a communication interface 120, and a memory 130.

[0100] The communication interface 120 is used for transmitting data.

[0101] The processor 110 is the control center of the server 100, connecting all parts of the server 100 through various interfaces and routes, executing various functions of the server 100 and processing data by running or executing software programs / modules stored in the memory 130 and calling data stored in the memory 130. Optionally, the processor 110 can include one or more processing units.

[0102] The memory 130 can be used to store software programs and modules, and the processor 110 executes various function applications and data processing by running the software programs and modules stored in the memory 130. The memory 130 can mainly include a program storage area and a data storage area, wherein the program storage area can store an operating system, application programs required by at least one function, etc., and the data storage area can store data created according to business processing, etc. In addition, the memory 130 can include a high-speed random access memory, and can also include a non-volatile memory, such as at least one magnetic disk storage device, a flash memory device, or other volatile solid-state memory device.

[0103] It should be noted that the structure shown in the above Figure 1 is only an example, and the embodiments of the present application are not limited thereto.

[0104] Based on the above description, Figure 2 An exemplary flow diagram of a Halo orbit satellite momentum wheel unloading attitude calculation method provided by the embodiments of the present application is shown, and the flow can be executed by a Halo orbit satellite momentum wheel unloading attitude calculation device.

[0105] As Figure 2 shown, the flow specifically includes:

[0106] Step 210: planning a satellite unloading time, determining a nominal control direction and a nominal controlled orbit required for the satellite to maintain a halo orbit configuration at the unloading time, and determining a nominal position and velocity in an inertial system when the satellite first crosses an XZ plane of a geocentric-inertial coordinate system on the nominal controlled orbit.

[0107] In the embodiments of the present application, the momentum wheel of the satellite will continuously accumulate energy, and before the rotation speed of the momentum wheel reaches the limit, momentum wheel unloading is needed to restore the initial rotation speed of the momentum wheel, so the time of momentum wheel unloading needs to be planned to avoid the rotation speed of the momentum wheel exceeding the limit due to too long an interval of satellite unloading time, or the unloading frequency being too high due to too short an interval of satellite unloading time, thereby wasting propellant. After the satellite unloading time is planned, the nominal control direction and the nominal controlled orbit required for the satellite to maintain a halo orbit configuration at the unloading time need to be determined. The halo orbit can also be referred to as a Halo orbit. Exemplarily, at the satellite unloading time, the nominal control information of the satellite maintaining the Halo orbit configuration is calculated according to the precise orbit determination result of the satellite by a Halo orbit maintenance control method. The Halo orbit maintenance control method in the present embodiment can be any method, and the Halo orbit maintenance control method is not specifically limited herein. The nominal control information can include a nominal control direction, a nominal velocity increment, and a nominal controlled orbit, wherein the nominal controlled orbit includes an epoch of the orbit and six elements of the orbit.

[0108] Then, based on the nominal post-control orbit, the nominal position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the geolunar rotating coordinate system on the nominal post-control orbit are determined. Specifically, first, the position and velocity of the satellite in the inertial system are dynamically integrated based on the nominal post-control orbit according to a specified step size, to obtain the position and velocity of the satellite at the kth time, where k is an integer greater than or equal to 1, and the first time is located after the unloading time and is separated from the unloading time by a preset time period. The specified step size is a value set according to experience, for example, 2 hours, 1 hour, etc., which is not specifically limited here. It can be understood that, in order to prevent the distance between the nominal control point of the satellite on the orbit and the crossing point of the XZ plane of the geolunar rotating coordinate system from being too close, resulting in a decrease in calculation accuracy, the first time needs to be located after the unloading time and separated from the unloading time by a preset time period, where the preset time period is a value set according to experience, for example, 2 days, and the value of the preset time period is not specifically limited here. Exemplarily, the orbit prediction is continuously performed according to the specified step size dT based on the nominal post-control orbit, the satellite is subjected to forces including the geolunar gravity, the gravity of other large celestial bodies, and the solar pressure, and the dynamic model is as follows:

[0109]

[0110] The position and velocity of the satellite at the unloading time are obtained based on the orbital elements of the nominal post-control orbit at the unloading time, and then the position and velocity of the satellite at the kth time are obtained by dynamically integrating the position and velocity of the satellite at the unloading time. It can be understood that, because the dynamic integration is performed according to the specified step size, the kth time and the k+1th time are separated by the specified step size.

[0111] Then, the rotation matrix of the inertial system converted to the geolunar rotating coordinate system at the kth time is calculated. Exemplarily, the halo orbit around the Lagrange point L2 is taken as an example in this paper, so the geolunar rotating coordinate system is the geolunar L2 point rotating coordinate system, and the geolunar L2 point rotating coordinate system at the kth time is taken as an example in the following, which is defined as follows: the origin O is the geolunar L2 point, the X axis is the direction of the earth pointing to the moon at the kth time, the Y axis is parallel to the direction of the moon's orbit tangent at the kth time, and the Z axis is parallel to the normal direction of the white line at the kth time. The rotation matrix L of the J2000 inertial system to the geolunar rotating coordinate system is calculated ri and the rotation angular velocity ω r , and the formula is as follows:

[0112]

[0113] wherein, is the relative position between the moon and the earth, is the relative velocity between the moon and the earth, m represents the moon, and e represents the earth, and are the position and velocity of the moon in the inertial system at the kth time, and are the position and velocity of the Earth at the kth time in the inertial system, and X, Y, Z are the unit vectors in the X, Y, Z directions of the Earth-Moon L2 point rotating coordinate system.

[0114] Then, based on the rotation matrix, the position and velocity at the kth time are converted from the inertial system to the Earth-Moon rotating coordinate system at the kth time, to obtain the position coordinates at the kth time in the Earth-Moon rotating coordinate system at the kth time. It can be understood that in addition to obtaining the position coordinates, the velocity at the kth time also needs to be converted from the inertial system to the Earth-Moon rotating coordinate system at the kth time, to obtain the velocity at the kth time in the Earth-Moon rotating coordinate system at the kth time. For example, the Earth-Moon distance rate of change V em :

[0115]

[0116] Then the position and velocity of the center of mass of the Earth-Moon system are determined:

[0117]

[0118] where M m is the weight of the Moon, M e is the weight of the Earth, is the position of the center of mass of the Earth-Moon system, is the velocity of the center of mass of the Earth-Moon system.

[0119] The position and velocity of the satellite are rotated from the inertial system to the rotating coordinate system:

[0120]

[0121] where, is the position of the satellite in the rotating coordinate system, is the velocity of the satellite in the rotating coordinate system, is the position of the satellite in the inertial system, is the velocity of the satellite in the inertial system.

[0122] Finally, the position and velocity of the satellite are rotated from the inertial system to the Earth-Moon L2 point rotating coordinate system:

[0123]

[0124] where, is the position of the satellite in the Earth-Moon L2 point rotating coordinate system, is the velocity of the satellite in the Earth-Moon L2 point rotating coordinate system, and L2 is the distance of the normalized L2 point relative to the center of mass of the Earth-Moon system.

[0125] Then, a product result of the Y-axis component of the position coordinate at the kth moment and the Y-axis component of the position coordinate at the k+1th moment is calculated. Exemplarily, the calculation formula is shown as follows:

[0126] flag=r eml，k (y)·r eml，k+1 (y)

[0127] wherein eml is a geocentric inertial coordinate system, r eml , k (y) is the Y-axis component of the position coordinate at the kth moment in the geocentric inertial coordinate system at the kth moment, r eml , k+1 (y) is the Y-axis component of the position coordinate at the k+1th moment in the geocentric inertial coordinate system at the k+1th moment, and flag is the product result.

[0128] When the product result is less than 0 for the first time, the satellite crosses the XZ plane of the geocentric inertial coordinate system for the first time on the nominal controlled post-orbit, and the kth moment and the k+1th moment are accurately iterated to obtain the nominal position and velocity of the satellite in the inertial system when the satellite crosses the XZ plane of the geocentric inertial coordinate system for the first time on the nominal controlled post-orbit. The method for accurately iterating the kth moment and the k+1th moment can be selected according to experience, for example, the 618 method, and the method for accurately iterating is not specifically limited herein. Exemplarily, whether the product result of the Y-axis component of the position coordinate at the 1th moment and the Y-axis component of the position coordinate at the 2th moment is less than 0 for the first time is judged, if not, whether the product result of the Y-axis component of the position coordinate at the 2th moment and the Y-axis component of the position coordinate at the 3th moment is less than 0 for the first time is judged, if not, the judgment is continued in turn until the product result is less than 0 for the first time. If the product result of the Y-axis component of the position coordinate at the 6th moment and the Y-axis component of the position coordinate at the 7th moment is less than 0 for the first time, it indicates that the satellite crosses the XZ plane of the geocentric inertial coordinate system for the first time on the nominal controlled post-orbit between the 6th moment and the 7th moment. Then, the 618 method is used to accurately iterate the 6th moment and the 7th moment, the time convergence threshold is set to 1.0e-2 seconds, and finally the nominal position and velocity of the satellite in the inertial system when the satellite crosses the XZ plane of the geocentric inertial coordinate system for the first time on the nominal controlled post-orbit are obtained

[0129] In a possible implementation, before planning the satellite unloading moment, the fitting coefficient of the satellite also needs to be determined. It can be understood that generally the Halo orbit satellite has two force-couple thrusters, and unloading almost does not produce velocity increment and has no effect on the orbit. The other direction is a non-force-couple thruster, which offsets the angular momentum by producing a velocity increment. It is assumed that the Z-axis direction is the non-force-couple thruster, and the X and Y axes are the force-couple thrusters. The angular momentum of the satellite in the body system at the unloading moment is H b =[H bx Hby H bz ],wherein b represents the unloading moment, the fitting relationship between the velocity increment and the angular momentum of the satellite in the body system at the unloading moment is:

[0130]

[0131] The determination process of the fitting coefficient is as follows: first, the historical unloading moment attitude quaternion, the historical unloading moment angular momentum and the corresponding historical unloading velocity increment are obtained. The historical unloading moment angular momentum is the three-axis angular momentum of the satellite in the body system at the unloading moment. The historical unloading velocity increment is the velocity increment generated by the satellite in the inertial system during the historical unloading process, and the historical unloading velocity increment is obtained according to the precise orbit determination results before and after the unloading of the satellite.

[0132] Then, based on the historical unloading moment attitude quaternion, the first conversion matrix from the inertial system to the body system at the unloading moment is determined. For example, according to the satellite attitude quaternion Q = [Q c , Q1, Q2, Q3] at the unloading moment, the conversion matrix M bi from the inertial system to the body system is calculated, and the formula is as follows:

[0133]

[0134] Finally, the historical unloading velocity increment is converted from the inertial system to the body system at the unloading moment based on the first conversion matrix, and the fitting coefficient is obtained according to the fitting relationship between the historical unloading velocity increment and the historical unloading moment angular momentum in the body system at the unloading moment. For example, the historical unloading velocity increment in the inertial system is converted to the body system at the unloading moment as Then let V b,orb = V b , that is, Finally, the fitting coefficient K can be obtained.

[0135] It can be understood that in order to ensure the accuracy of the calculated fitting coefficient, a large amount of historical unloading moment data needs to be used for calculation, that is, for any historical unloading moment, the fitting coefficient corresponding to the historical unloading moment is calculated according to the historical unloading moment attitude quaternion, the historical unloading moment angular momentum and the corresponding historical unloading velocity increment at the historical unloading moment. Finally, the fitting coefficients corresponding to all historical unloading moments are fitted by the least square method to determine the final fitting coefficient K. And through the above calculation process, it can be seen that the fitting coefficient represents the fitting relationship between the angular momentum of the satellite in the body system at the unloading moment and the velocity increment generated during the unloading process.

[0136] Step 220, establishing an angular momentum coordinate system of the satellite based on the angular momentum vector in the inertial system at the unloading moment and the nominal control direction, and determining an angle optimization variable of the satellite, the angle optimization variable being Euler angles of rotation of the satellite around each coordinate axis of the angular momentum coordinate system according to a rotation sequence.

[0137] In the embodiment of the present application, the angular momentum vector in the inertial system at the unloading moment can be obtained according to the rotation speed of the momentum wheel at the unloading moment. The angle optimization variable of the satellite is Euler angles of rotation of the satellite around each coordinate axis of the angular momentum coordinate system according to a rotation sequence. The rotation sequence of the Euler angles is set according to experience, for example, the rotation sequence is 123 or 213, which is not limited herein. Specifically, the direction of the angular momentum vector in the inertial system at the unloading moment is taken as the Y axis of the angular momentum coordinate system, which is expressed as The cross product of the Y axis and the nominal control direction is taken as the X axis, which is expressed as The Z axis is determined based on the right-hand rule, which is expressed as z=x×y; the angular momentum coordinate system of the satellite is established, and a second conversion matrix from the inertial system to the angular momentum coordinate system is determined, which is expressed as M jdl,i =[x y z] T , jdl represents the angular momentum coordinate system.

[0138] The angle optimization variable of the satellite is determined based on the rotation angle of the satellite around each coordinate axis of the angular momentum coordinate system, and a third conversion matrix from the angular momentum coordinate system to the body system at the unloading moment is determined. For example, first, the satellite is rotated around the X axis of the angular momentum coordinate system by The conversion matrix is Then, the satellite is rotated around the Y axis by The conversion matrix is Finally, the satellite is rotated around the Z axis by The unloading attitude is obtained, and the conversion matrix is Wherein, M x , M y , and M z are elementary conversion matrices. The angle optimization variable is determined as According to the Euler angle 123 rotation sequence, the satellite is rotated in sequence Then, the third conversion matrix M b,jdl from the angular momentum coordinate system to the body system at the unloading moment is:

[0139]

[0140] Finally, based on the second conversion matrix and the third conversion matrix, a fourth conversion matrix from the inertial system to the body system at the unloading moment is obtained. For example, the cross product between the second conversion matrix M jdl,i and the third conversion matrix M b,jdl is calculated, which is taken as the fourth conversion matrix M b,i =Mb,jdl ×M jdl,i .

[0141] In step 230, the velocity increment generated by the satellite during the unloading process is determined based on the angular momentum vector and the fitting coefficient, and the post-unloading orbit and the post-unloading position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the Earth-Moon rotating coordinate system on the post-unloading orbit are determined; the fitting coefficient represents the fitting relationship between the angular momentum of the satellite at the unloading moment and the velocity increment generated during the unloading process.

[0142] In the embodiment of the application, it can be known through the foregoing calculation process of the fitting coefficient that the fitting coefficient represents the fitting relationship between the angular momentum of the satellite at the unloading moment and the velocity increment generated during the unloading process in the local system at the unloading moment, and therefore, when the velocity increment generated during the unloading process is determined based on the angular momentum vector at the unloading moment and the fitting coefficient, coordinate system conversion is required. Specifically, the angular momentum vector is converted from the inertial system to the local system at the unloading moment based on the fourth conversion matrix, and the velocity increment generated during the unloading process of the satellite in the local system at the unloading moment is obtained based on the angular momentum vector in the local system at the unloading moment and the fitting coefficient. For example, the angular momentum vector H i in the inertial system is calculated as b =M b,i ×H i =(H b,x H b,y H b,z ). Then, the velocity increment in the local system at the unloading moment is calculated based on the fitting coefficient K and the angular momentum vector H b in the local system at the unloading moment.

[0143] Then, the velocity increment generated during the unloading process is converted from the local system at the unloading moment to the inertial system based on the inverse matrix of the fourth conversion matrix, as represented by the formula , and the post-unloading orbit is determined based on the velocity increment generated during the unloading process. The post-unloading orbit includes the epoch of the orbit and the six elements of the orbit.

[0144] Finally, the post-unloading position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the lunar rotating coordinate system on the post-unloading orbit are determined based on the post-unloading orbit. Specifically, the position and velocity of the satellite in the inertial system are dynamically integrated based on the post-unloading orbit according to a specified step length to obtain the position and velocity of the satellite at the mth moment, and a rotation matrix for converting from the inertial system to the lunar rotating coordinate system at the mth moment is calculated. Here, m is an integer greater than or equal to 1, and the 1st moment is located after the unloading moment and is separated from the unloading moment by a preset time period. Based on the rotation matrix, the position and velocity at the mth moment are converted from the inertial system to the lunar rotating coordinate system at the mth moment to obtain the position coordinates of the satellite in the lunar rotating coordinate system at the mth moment. The product of the Y-axis component of the position coordinates at the mth moment and the Y-axis component of the position coordinates at the (m+1)th moment is calculated. When the product is less than 0 for the first time, the satellite first crosses the XZ plane of the lunar rotating coordinate system on the post-unloading orbit, and the mth moment and the (m+1)th moment when the product is less than 0 for the first time are accurately iterated to obtain the post-unloading position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the lunar rotating coordinate system on the post-unloading orbit. It can be understood that the calculation process of the post-unloading position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the lunar rotating coordinate system on the post-unloading orbit is the same as the calculation process of the nominal position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the lunar rotating coordinate system on the nominal post-control orbit in step 210 described above, and thus will not be described here. Finally, the post-unloading position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the lunar rotating coordinate system on the post-unloading orbit is obtained

[0145] In step 240, the deviation between the nominal post-control orbit and the post-unloading orbit is obtained according to the nominal position and velocity and the post-unloading position and velocity, and the deviation is optimized according to the angle optimization variable based on the coordinate permutation method to obtain the optimal attitude quaternion of the satellite momentum wheel unloading, so that the satellite performs momentum wheel unloading according to the optimal attitude quaternion, and the goal of maintaining the halo orbit configuration by active unloading is achieved.

[0146] In the embodiment of the present application, the deviation between the nominal control orbit and the unloading orbit can be obtained through the deviation between the nominal position velocity and the position velocity after unloading, that is, the deviation between the nominal control and the unloading control of the satellite. Specifically, the product of the modulus of the vector difference between the satellite position in the position velocity after unloading and the satellite position in the nominal position velocity and the position weight is calculated as the position deviation. The product of the modulus of the vector difference between the satellite velocity in the position velocity after unloading and the satellite velocity in the nominal position velocity and the velocity weight is calculated as the velocity deviation. The position weight is less than the velocity weight, and the sum of the position weight and the velocity weight is 1. It can be understood that, since the nominal position velocity and the position velocity after unloading are both the position velocities when the satellite first crosses the XZ plane of the geolunar rotation coordinate system, the XZ plane of the geolunar rotation coordinate system is the same, and therefore the weight of the position deviation is smaller. Finally, the sum of the position deviation and the velocity deviation is calculated as the deviation between the nominal control orbit and the unloading orbit. For example, the deviation between the nominal control orbit and the unloading orbit is calculated by the following formula:

[0147]

[0148] wherein S err is the deviation between the nominal control orbit and the unloading orbit when the satellite is at the rotation angle n is the velocity weight; (1-n) is the position weight; is the satellite position in the position velocity after unloading; is the satellite position in the nominal position velocity; is the satellite velocity in the position velocity after unloading; is the satellite velocity in the nominal position velocity.

[0149] Then, based on the coordinate permutation method, the deviation between the nominal control orbit and the unloading orbit is optimized according to the angle optimization variable to obtain the optimal attitude quaternion of the satellite momentum wheel unloading. Specifically, the X axis of the angular momentum coordinate system is divided into positive and negative directions to obtain two half spaces, and the value range of the angle optimization variable in each half space is determined. For example, the value range of the angle optimization variable in half space 1 is the value range of the angle optimization variable in half space 2 is For any one of the two half spaces, based on the coordinate permutation method, the deviation between the nominal control orbit and the unloading orbit is optimized according to the value range of the angle optimization variable in the half space until the deviation is minimized to obtain the optimal value of the angle optimization variable corresponding to the half space and the minimum deviation. Finally, the optimal value of the angle optimization variable corresponding to the smaller minimum deviation of the two half spaces is taken as the optimal value of the angle optimization variable of the satellite, such as

[0150]

[0151] Based on the fourth transformation matrix corresponding to the optimal value of the satellite's angle optimization variable, the optimal attitude quaternion for unloading the satellite's momentum wheel is obtained. For example, the optimal value of the satellite's angle optimization variable... The corresponding fourth transformation matrix is Transform the fourth transformation matrix corresponding to the optimal value of the angle optimization variable to obtain the optimal attitude quaternion for momentum wheel unloading. By unloading the momentum wheel according to the optimal attitude quaternion, the satellite can achieve the goal of maintaining the Halo orbit configuration through active unloading.

[0152] The following is about Figure 2 The following examples illustrate the implementation:

[0153] Taking a Halo-orbited satellite at the Earth-Moon L2 point as an example, we first fit the relationship between the satellite's three-axis angular momentum at the unloading moment and the velocity increment generated during the unloading process. Using the satellite's unloading data on January 3, 2024, at 11:30 AM, the satellite's intrinsic angular momentum at that time is:

[0154] X axis Y axis Z axis System angular momentum 0.848391 -2.809106 -0.28722

[0155] The velocity increment caused by unloading under the satellite's own system is V. b,z =K(H b,x +H b,y = K × 3.657497.

[0156] The attitude quaternion of the satellite at this unloading moment is:

[0157] Attitude quaternion Value QC 0.28699928 Q1 -0.66184163 Q2 -0.61272231 Q3 -0.32275136

[0158] The transformation matrix of the satellite from the inertial frame to its own frame, calculated using attitude quaternions, is:

[0159]

[0160] Based on the precise orbit determination results before and after satellite unloading, estimate the velocity increment V generated in the inertial frame by satellite unloading. i,orb =(0.0012 + 0.0125 - 0.01) T Based on the above transformation matrix, it is transformed to the current system as follows:

[0161]

[0162] Let V b,orb =V b The fitting coefficient K = 0.00438893 can then be obtained.

[0163] The satellite long-term on-orbit unloading angular momentum data and precise orbit determination results before and after unloading are calculated, the relationship between the satellite three-axis angular momentum at the unloading time and the velocity increment generated by unloading is fitted, and the fitting coefficient K is obtained.

[0164] Taking the satellite unloading on December 5, 2024 as an example for calculation. The satellite orbit elements and initial attitude quaternion are as follows, epoch 2024-12-05T10:00:00.

[0165] Orbital elements Value Semi-major axis (m) 841364788.8626 Eccentricity 0.4781437406329836 Orbital inclination (°) 30.4744 Longitude of ascending node (°) 356.1490 Argument of perigee (°) 314.4349 Longitude of ascending node (°) 354.5827

[0166]

[0167]

[0168] The planned unloading time is 2024-12-05T12:30:00, and the three-axis telemetry angular momentum at this time is (1.224229-2.362641-1.686074). If the satellite is orbitally maintained at this time, the nominal control attitude quaternion is as follows, and the nominal velocity increment is 0.0154 m / s.

[0169] Argument of perigee (°) Nominal control attitude quaternion Value 0.30159259 Q1 0.72125952 Q2 0.22798981 Q3 0.58038544

[0170] The nominal control orbit is as follows:

[0171] QC Orbital elements Value 859972472.1059 Semi-major axis (m) 0.4903719701563828 Eccentricity 30.4881 Orbital inclination (°) 356.0979 Longitude of ascending node (°) 315.4140 Argument of perigee (°) 354.9126

[0172] Taking the nominal control orbit as the reference, the step is taken as 1 hour, and the time when the satellite first crosses the XZ plane of the earth-moon rotating coordinate system and the position and velocity in the L2 center rotating coordinate system are calculated as follows:

[0173]

[0174] The time when the satellite first crosses the XZ plane of the earth-moon rotating coordinate system and the position and velocity in the L2 center rotating coordinate system are accurately iterated by the 618 method, and the results are as follows:

[0175] Argument of perigee (°) Item (unit) Value 2024-12-09T10:21:02.56 Epoch -15560976.552 X (m) -0.258 Y (m) 10710243.526 Z (m) 0.311948 VX (m / s) 203.354995 VY (m / s) -0.660270

[0176] The satellite angular momentum coordinate system is established according to the contents in the above step 220, and the conversion matrix from the inertial system to the angular momentum coordinate system is calculated.

[0177] The deviation between the nominal control orbit and the unloading orbit is obtained by the method in the above steps 230 and 240. Based on the coordinate permutation method, the deviation is optimized, and the different angle optimization variables are rotated in turn according to the Euler angle 123 rotation sequence in the angular momentum coordinate system The optimal value of the satellite angular optimization variable is obtained For:

[0178]

[0179] At this time, the optimal attitude quaternion of the satellite is as follows:

[0180] VZ (m / s) Attitude quaternion Value 0.72379567 Q1 -0.64775000 Q2 -0.16537506 Q3 0.17085330

[0181] In the unloading attitude coordinate system, the satellite angular momentum vector is (-1.236071 -1.235012 -1.753562), and the velocity increment generated after unloading is 0.024073 m / s. If the satellite maintains the orbit according to the optimal attitude quaternion, the required velocity increment is 0.024082 m / s, so the residual velocity increment of the satellite after unloading is only 9e-6 m / s.

[0182] In the embodiment of the application, based on the long-term on-orbit unloading angular momentum data of the satellite and the precise orbit determination results before and after unloading, the relationship between the angular momentum of the satellite at the unloading time and the velocity increment generated during the unloading process is fitted. When the satellite is unloaded, the nominal control parameters for maintaining the Halo orbit configuration are calculated, and the nominal position and velocity when passing through the XZ plane of the geocentric lunar rotating coordinate system for the first time after control are calculated. According to the angular momentum vector at the unloading time and the nominal control direction, an angular momentum coordinate system is defined, and the unloading attitude is obtained by rotating the angle optimization variable in the order of Euler angles 123. The three-axis distribution of the angular momentum under the unloading attitude is calculated, the velocity increment generated during the unloading is calculated according to the fitting relationship, the post-unloading orbit under the unloading attitude is obtained, the post-unloading position and velocity when passing through the XZ plane of the geocentric lunar rotating coordinate system for the first time after control are calculated, and the deviation between the post-unloading orbit and the nominal post-control orbit is evaluated based on the post-unloading position and velocity and the nominal position and velocity. The optimal attitude quaternion is calculated by using the coordinate wheel method, so that the deviation between the post-unloading orbit and the nominal post-control orbit is minimized. The satellite can maintain the Halo orbit configuration to the greatest extent by using the optimal attitude quaternion for long-term momentum wheel unloading, so as to realize the use of the velocity increment generated during the unloading for the satellite orbit configuration maintenance, offset the divergence trend of the Halo orbit by using the momentum wheel unloading, reduce the orbit maintenance frequency, and reduce the consumption of propellant of the satellite.

[0183] Based on the same technical concept, QC An exemplary structure schematic diagram of a Halo orbit satellite momentum wheel unloading attitude calculation device provided by the embodiment of the application is shown, and the device can execute the flow of the Halo orbit satellite momentum wheel unloading attitude calculation method.

[0184] As Figure 3 shown, the device specifically includes:

[0185] The acquisition module 310 is configured to plan a satellite unloading time, determine a nominal control direction and a nominal post-control orbit required for the satellite to maintain a halo orbit configuration at the unloading time, and determine a nominal position and velocity of the satellite in an inertial system when the satellite first crosses an XZ plane of a geocentric lunar rotating coordinate system on the nominal post-control orbit based on the nominal post-control orbit.

[0186] The processing module 320 is configured to establish an angular momentum coordinate system of the satellite based on an angular momentum vector in the inertial system at the unloading time and the nominal control direction, and determine an angular optimization variable of the satellite, the angular optimization variable being Euler angles of the satellite rotating around coordinate axes of the angular momentum coordinate system in a rotation order.

[0187] The processing module 320 is configured to determine a velocity increment generated by the satellite during the unloading process based on the angular momentum vector and the fitting coefficient, and determine a post-unloading orbit and a post-unloading position and velocity of the satellite in the inertial system when the satellite first crosses the XZ plane of the geocentric lunar rotating coordinate system on the post-unloading orbit, the fitting coefficient representing a fitting relationship between the angular momentum of the satellite at the unloading time and the velocity increment generated by the satellite during the unloading process.

[0188] The processing module 320 is configured to obtain a deviation between the nominal post-control orbit and the post-unloading orbit according to the nominal position and velocity and the post-unloading position and velocity, and optimize the deviation based on a coordinate permutation method according to the angular optimization variable to obtain an optimal attitude quaternion of the satellite momentum wheel unloading, so that the satellite performs the momentum wheel unloading according to the optimal attitude quaternion, and achieves the goal of maintaining the halo orbit configuration through active unloading.

[0189] Optionally, the processing module 320 is further configured to:

[0190] The processing module 320 is configured to obtain a historical unloading time attitude quaternion, a historical unloading time angular momentum, and a corresponding historical unloading velocity increment, the historical unloading velocity increment being a velocity increment generated by the satellite in the inertial system during a historical unloading process, the historical unloading velocity increment being obtained according to precise orbit determination results before and after unloading of the satellite, and the historical unloading time angular momentum being three-axis angular momentum of the satellite in a local system at the unloading time.

[0191] The processing module 320 is configured to determine a first conversion matrix from the inertial system to the local system at the unloading time based on the historical unloading time attitude quaternion.

[0192] The processing module 320 is configured to convert the historical unloading velocity increment from the inertial system to the local system at the unloading time based on the first conversion matrix, and obtain a fitting coefficient according to a fitting relationship between the historical unloading velocity increment and the historical unloading time angular momentum in the local system at the unloading time.

[0193] Optionally, the processing module 320 is specifically configured to:

[0194] performing dynamic integration on the position and velocity of the satellite in the inertial system based on the nominal post-control orbit according to a specified step size, obtaining position and velocity of the satellite at the kth moment, and calculating a rotation matrix converted from the inertial system to the geolunar rotation coordinate system at the kth moment; k is an integer greater than or equal to 1, and the first moment is located after the unloading moment and is separated from the unloading moment by a preset time period;

[0195] converting the position and velocity at the kth moment from the inertial system to the geolunar rotation coordinate system at the kth moment based on the rotation matrix, and obtaining position coordinates of the satellite at the kth moment in the geolunar rotation coordinate system at the kth moment;

[0196] calculating a product result of a Y-axis component of the position coordinates at the kth moment and a Y-axis component of the position coordinates at the k+1th moment;

[0197] when the product result is less than 0 for the first time, the satellite crosses the XZ plane of the geolunar rotation coordinate system for the first time on the nominal post-control orbit, and the kth moment and the k+1th moment are accurately iterated to obtain a nominal position and velocity of the satellite in the inertial system when the satellite crosses the XZ plane of the geolunar rotation coordinate system for the first time on the nominal post-control orbit.

[0198] Optionally, the processing module 320 is specifically configured to:

[0199] taking the direction of the angular momentum vector as a Y-axis of the angular momentum coordinate system, taking the direction of the cross product result of the Y-axis and the nominal control direction as an X-axis, determining a Z-axis based on the right-hand rule, establishing the angular momentum coordinate system, and determining a second conversion matrix converted from the inertial system to the angular momentum coordinate system;

[0200] determining the angle optimization variable based on the angle of rotation of the satellite around each coordinate axis of the angular momentum coordinate system, and determining a third conversion matrix from the angular momentum coordinate system to the body coordinate system at the unloading moment;

[0201] based on the second conversion matrix and the third conversion matrix, obtaining a fourth conversion matrix from the inertial system to the body coordinate system at the unloading moment.

[0202] Optionally, the processing module 320 is specifically configured to:

[0203] based on the fourth conversion matrix, converting the angular momentum vector from the inertial system to the body coordinate system at the unloading moment, and based on the angular momentum vector in the body coordinate system at the unloading moment and the fitting coefficient, obtaining a velocity increment generated by the satellite in the unloading process in the body coordinate system at the unloading moment;

[0204] Converting the velocity increment generated in the unloading process from the body system at the unloading time to the inertial system based on an inverse matrix of the fourth conversion matrix, and determining an orbit after unloading based on the velocity increment generated in the unloading process;

[0205] Performing dynamic integration on the position and velocity of the satellite in the inertial system based on the orbit after unloading according to a specified step length to obtain position and velocity of the satellite at the mth time, and calculating a rotation matrix for converting from the inertial system to a geocentric inertial coordinate system at the mth time; m is an integer greater than or equal to 1, and the 1st time is located after the unloading time and is separated from the unloading time by a preset time period;

[0206] Converting the position and velocity at the mth time from the inertial system to the geocentric inertial coordinate system at the mth time based on the rotation matrix to obtain position coordinates of the satellite at the mth time in the geocentric inertial coordinate system at the mth time;

[0207] Calculating a product result of a Y-axis component of the position coordinates at the mth time and a Y-axis component of position coordinates at the m+1th time;

[0208] When the product result is less than 0 for the first time, the satellite crosses the XZ plane of the geocentric inertial coordinate system for the first time on the orbit after unloading, and the mth time and the m+1th time when the product result is less than 0 for the first time are iterated accurately to obtain the position and velocity of the satellite in the inertial system after unloading when the satellite crosses the XZ plane of the geocentric inertial coordinate system for the first time on the orbit after unloading.

[0209] Optionally, the processing module 320 is specifically configured to:

[0210] Calculate a product result between a modulus value of a vector difference between a position of the satellite in the position and velocity after unloading and a position of the satellite in the nominal position and velocity and a position weight as a position deviation;

[0211] Calculate a product result between a modulus value of a vector difference between a velocity of the satellite in the position and velocity after unloading and a velocity of the satellite in the nominal position and velocity and a velocity weight as a velocity deviation; the position weight is less than the velocity weight, and the position weight and the velocity weight are equal to 1;

[0212] Calculate a summation result of the position deviation and the velocity deviation as a deviation between the nominal control orbit and the orbit after unloading.

[0213] Optionally, the processing module 320 is specifically configured to:

[0214] Divide the X-axis positive direction and the X-axis negative direction of the angular momentum coordinate system to obtain two half spaces, and determine a value range of the angular optimization variable in each half space;

[0215] For any one of the two half spaces, based on coordinate transformation method, according to the value range of the angle optimization variable in the half space, the bias is optimized according to the angle optimization variable until the bias is minimum, and the optimal value of the angle optimization variable corresponding to the half space and the minimum bias are obtained;

[0216] Among the minimum biases corresponding to the two half spaces, the optimal value of the angle optimization variable corresponding to the smaller minimum bias is taken as the optimal value of the angle optimization variable of the satellite;

[0217] Based on the fourth conversion matrix corresponding to the optimal value of the angle optimization variable of the satellite, the optimal attitude quaternion of the momentum wheel unloading of the satellite is obtained.

[0218] Based on the same technical concept, the embodiment of the present application further provides a computer device, comprising:

[0219] a memory for storing program instructions;

[0220] a processor for calling the program instructions stored in the memory, and performing the method for calculating the momentum wheel unloading attitude of the Halo orbit satellite according to the obtained program.

[0221] Based on the same technical concept, the embodiment of the present application further provides a computer readable storage medium, which stores computer executable instructions, and the computer executable instructions are used for making the computer execute the method for calculating the momentum wheel unloading attitude of the Halo orbit satellite.

[0222] Based on the same technical concept, the embodiment of the present application further provides a computer program product, characterized in that the computer program product comprises an executable program, and the executable program is used for executing the method for calculating the momentum wheel unloading attitude of the Halo orbit satellite by the processor.

[0223] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, a system or a computer program product. Therefore, the present application can adopt a complete hardware embodiment, a complete software embodiment or an embodiment combining software and hardware aspects. Moreover, the present application can adopt a computer program product implemented on one or more computer usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer usable program codes.

[0224] The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart block or blocks. Figure 3 one or more flow or blocks Figure 1 means for functionally implementing the steps in one or more flow or blocks

[0225] These computer program instructions can also be stored in a computer- readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instructions which implement the function specified in the flowchart block or blocks. Figure 1 one or more flow or blocks Figure 1 means for functionally implementing the steps in one or more flow or blocks

[0226] The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart block or blocks. Figure 1 one or more flow or blocks Figure 1 Figure 1 means for functionally implementing the steps in one or more flow or blocks

[0227] Obviously, numerous modifications and variations of the present application are possible in light of the above teachings. It is therefore to be understood that within the scope of the appended claims and their equivalents, the application can be practiced otherwise than as specifically described.

Claims

1. A method for calculating the momentum wheel unloading attitude of a Halo-orbiting satellite, characterized in that, The method comprises the following steps: planning a satellite unloading time, determining a nominal control direction and a nominal post-control orbit required for the satellite to maintain a halo orbit configuration at the unloading time, and determining a nominal position and velocity in an inertial system when the satellite first crosses an XZ plane of a geolunar rotating coordinate system on the nominal post-control orbit based on the nominal post-control orbit; establishing an angular momentum coordinate system of the satellite based on an angular momentum vector in the inertial system at the unloading time and the nominal control direction, and determining an angular optimization variable of the satellite, the angular optimization variable being Euler angles of the satellite rotating around each coordinate axis of the angular momentum coordinate system according to a rotation sequence; determining a velocity increment generated by the satellite during the unloading process based on the angular momentum vector and a fitting coefficient, and determining a post-unloading orbit and a post-unloading position and velocity in the inertial system when the satellite first crosses the XZ plane of the geolunar rotating coordinate system on the post-unloading orbit; the fitting coefficient represents a fitting relationship between the angular momentum of the satellite at the unloading time and the velocity increment generated during the unloading process; obtaining a deviation between the nominal post-control orbit and the post-unloading orbit according to the nominal position and velocity and the post-unloading position and velocity, and optimizing the deviation according to the angular optimization variable based on a coordinate permutation method to obtain an optimal attitude quaternion of a momentum wheel unloading of the satellite, so that the satellite performs the momentum wheel unloading according to the optimal attitude quaternion, and the goal of maintaining the halo orbit configuration through active unloading is achieved.

2. The method of claim 1, wherein, Before planning the satellite unloading time, the method further comprises the following steps: obtaining a historical unloading time attitude quaternion, a historical unloading time angular momentum and a corresponding historical unloading velocity increment; the historical unloading velocity increment is a velocity increment generated by the satellite in the inertial system during a historical unloading process, the historical unloading velocity increment is obtained according to precise orbit determination results before and after the unloading of the satellite, and the historical unloading time angular momentum is three-axis angular momentum in a local system of the satellite at the unloading time; determining a first conversion matrix from the inertial system to the local system at the unloading time based on the historical unloading time attitude quaternion; converting the historical unloading velocity increment from the inertial system to the local system at the unloading time based on the first conversion matrix, and obtaining a fitting coefficient according to a fitting relationship between the historical unloading velocity increment and the historical unloading time angular momentum in the local system at the unloading time.

3. The method of claim 1, wherein, determining a nominal position and velocity in an inertial system when the satellite first crosses an XZ plane of a geolunar rotating coordinate system on the nominal post-control orbit based on the nominal post-control orbit, comprises the following steps: performing dynamic integration on the position and velocity of the satellite in the inertial system based on the nominal post-control orbit according to a specified step length to obtain a position and velocity of the satellite at a kth time, and calculating a rotation matrix from the inertial system to a geolunar rotating coordinate system at the kth time; k is an integer greater than or equal to 1, and the 1st time is located after the unloading time and is separated from the unloading time by a preset time period; converting the position and velocity at the kth time from the inertial system to the geolunar rotating coordinate system at the kth time based on the rotation matrix to obtain position coordinates of the kth time in the geolunar rotating coordinate system at the kth time; a product of a Y-axis component of the position coordinate at the kth moment and a Y-axis component of the position coordinate at the k+1th moment; when the product result is less than 0 for the first time, the satellite crosses the XZ plane of the geolunar rotating coordinate system for the first time on the nominal controlled post-orbit, and the kth moment and the k+1th moment are accurately iterated to obtain the nominal position and velocity of the satellite in the inertial system when the satellite crosses the XZ plane of the geolunar rotating coordinate system for the first time on the nominal controlled post-orbit.

4. The method of claim 1, wherein, an angular momentum coordinate system of the satellite is established based on the angular momentum vector in the inertial system at the unloading moment and the nominal control direction, and an angular optimization variable of the satellite is determined, including: the direction of the angular momentum vector is taken as a Y-axis of the angular momentum coordinate system, the direction of a cross product result of the Y-axis and the nominal control direction is taken as an X-axis, a Z-axis is determined based on a right-hand rule, the angular momentum coordinate system is established, and a second conversion matrix converted from the inertial system to the angular momentum coordinate system is determined; the angular optimization variable is determined based on angles of rotation of the satellite around each coordinate axis of the angular momentum coordinate system, and a third conversion matrix converted from the angular momentum coordinate system to the body system at the unloading moment is determined; a fourth conversion matrix converted from the inertial system to the body system at the unloading moment is obtained based on the second conversion matrix and the third conversion matrix.

5. The method of claim 4, wherein, a velocity increment generated in the unloading process of the satellite is determined based on the angular momentum vector and the fitting coefficient, and a post-unloading orbit and a post-unloading position and velocity of the satellite in the inertial system when the satellite crosses the XZ plane of the geolunar rotating coordinate system for the first time on the post-unloading orbit are determined, including: the angular momentum vector is converted from the inertial system to the body system at the unloading moment based on the fourth conversion matrix, and a velocity increment generated in the unloading process of the satellite in the body system at the unloading moment is obtained based on the angular momentum vector in the body system at the unloading moment and the fitting coefficient; the velocity increment generated in the unloading process is converted from the body system at the unloading moment to the inertial system based on an inverse matrix of the fourth conversion matrix, and the post-unloading orbit is determined based on the velocity increment generated in the unloading process; a position and velocity of the satellite in the inertial system are dynamically integrated based on the post-unloading orbit according to a specified step length to obtain a position and velocity of the satellite at an mth moment, and a rotation matrix converted from the inertial system to a geolunar rotating coordinate system at the mth moment is calculated; m is an integer greater than or equal to 1, and the 1st moment is located after the unloading moment and is separated from the unloading moment by a preset time period; the position and velocity at the mth moment are converted from the inertial system to the geolunar rotating coordinate system at the mth moment based on the rotation matrix to obtain a position coordinate at the mth moment in the geolunar rotating coordinate system at the mth moment; a product of a Y-axis component of the position coordinate at the mth moment and a Y-axis component of the position coordinate at the m+1th moment is calculated. When the product result is less than 0 for the first time, the satellite crosses the XZ plane of the geolunar rotation coordinate system for the first time on the post-unloading orbit, and the mth moment and the m+1th moment when the product result is less than 0 for the first time are accurately iterated to obtain the post-unloading position and velocity of the satellite in the inertial system when the satellite crosses the XZ plane of the geolunar rotation coordinate system for the first time on the post-unloading orbit.

6. The method according to any one of claims 1 to 4, wherein, According to the nominal position and velocity, the post-unloading position and velocity, the deviation between the nominal post-control orbit and the post-unloading orbit is obtained, including: The product result of the modulus value of the vector difference between the satellite position in the post-unloading position and velocity and the satellite position in the nominal position and velocity and the position weight is calculated as the position deviation; The product result of the modulus value of the vector difference between the satellite velocity in the post-unloading position and velocity and the satellite velocity in the nominal position and velocity and the velocity weight is calculated as the velocity deviation; the position weight is less than the velocity weight, and the sum of the position weight and the velocity weight is 1; The sum result of the position deviation and the velocity deviation is calculated as the deviation between the nominal post-control orbit and the post-unloading orbit.

7. The method of claim 4, wherein, Based on the coordinate permutation method, the deviation is optimized according to the angle optimization variable to obtain the optimal attitude quaternion of the satellite momentum wheel unloading, including: According to the positive direction and the negative direction of the X axis of the angular momentum coordinate system, two half spaces are obtained, and the value range of the angle optimization variable in each half space is determined; For any half space in the two half spaces, based on the coordinate permutation method, according to the value range of the angle optimization variable in the half space, the deviation is optimized according to the angle optimization variable until the deviation is minimized to obtain the optimal value of the angle optimization variable corresponding to the half space and the minimum deviation; The optimal value of the angle optimization variable corresponding to the smaller minimum deviation of the two half spaces is taken as the optimal value of the angle optimization variable of the satellite; Based on the fourth conversion matrix corresponding to the optimal value of the angle optimization variable of the satellite, the optimal attitude quaternion of the satellite momentum wheel unloading is obtained.

8. A calculation device for the momentum wheel unloading attitude of a Halo orbit satellite, characterized in that, Including: The acquisition module is configured to plan a satellite unloading time, determine a nominal control direction and a nominal post-control orbit required for the satellite to maintain a halo orbit configuration at the unloading time, and determine a nominal position and velocity in an inertial system when the satellite crosses an XZ plane of a geolunar rotation coordinate system for the first time on the nominal post-control orbit based on the nominal post-control orbit. The processing module is configured to establish an angular momentum coordinate system of the satellite based on an angular momentum vector in the inertial system at the unloading time and the nominal control direction, and determine an angle optimization variable of the satellite, the angle optimization variable being Euler angles of rotation of the satellite around each coordinate axis of the angular momentum coordinate system in a rotation sequence. The processing module is configured to determine a velocity increment generated by the satellite during the unloading process based on the angular momentum vector and a fitting coefficient, and determine a post-unloading orbit and a post-unloading position and velocity in the inertial system when the satellite crosses the XZ plane of the geolunar rotation coordinate system for the first time on the post-unloading orbit; the fitting coefficient represents a fitting relationship between the angular momentum of the satellite at the unloading time and the velocity increment generated during the unloading process. According to the nominal position velocity, the post-unloading position velocity, a deviation between the nominal post-control orbit and the post-unloading orbit is obtained, and based on a coordinate wheel method, the deviation is optimized according to the angle optimization variable, so as to obtain an optimal attitude quaternion of the satellite momentum wheel unloading, so that the satellite performs momentum wheel unloading according to the optimal attitude quaternion, and a target of maintaining the halo orbit configuration through active unloading is achieved.

9. A computer device, comprising: The method comprises the following steps: a memory for storing program instructions; a processor for calling the program instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer readable storage medium stores computer executable instructions for causing a computer to execute the method according to any one of claims 1 to 7.

11. A computer program product, characterised in that, The computer program product comprises an executable program which is executed by a processor to implement the method according to any one of claims 1 to 7.

Images