Geostationary satellite multi-pulse high-precision autonomous position keeping control method

Abstract

The application discloses a kind of geostationary satellite multi-pulse high-precision autonomous position keeping control methods, comprising: according to the instantaneous orbital elements of satellite, the average orbital elements of satellite are determined;According to the average orbital elements of satellite, the in-plane and out-of-plane pulse control sequence of satellite orbit is determined;According to the in-plane and out-of-plane pulse control sequence of satellite orbit, the switch logic of thruster is determined;According to the switch logic of thruster, the switch of thruster is controlled, realizes geostationary satellite multi-pulse high-precision autonomous position keeping control.The method of the application realizes geostationary satellite multi-pulse high-precision autonomous position keeping control.

Description

Technical Field

[0001] This invention belongs to the field of geostationary orbit satellite control technology, and particularly relates to a multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites. Background Technology

[0002] Geostationary meteorological satellites require high-precision attitude control and high-precision orbital position maintenance control to facilitate payload calibration. Limited by the number of thruster control cycles and control time, and affected by Earth's non-spherical perturbations, solar radiation pressure perturbations, and lunar gravitational perturbations, the satellite needs to perform high-precision position maintenance control daily.

[0003] Currently, methods for autonomous position-keeping control of satellites mainly focus on control methods using conical electric thrusters. These methods typically do not consider the effects of Earth's non-spherical perturbations, solar radiation pressure perturbations, and gravitational perturbations from the Sun and Moon. Furthermore, the control process requires solving optimization functions, resulting in a significant computational burden. Summary of the Invention

[0004] The technical problem solved by this invention is to overcome the shortcomings of the prior art and provide a multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites, aiming to achieve multi-pulse high-precision autonomous position-keeping control for geostationary orbit satellites.

[0005] To address the aforementioned technical problems, this invention discloses a multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites, comprising:

[0006] Determine the satellite's average orbital elements based on its instantaneous orbital elements;

[0007] Based on the average orbital features of the satellite, determine the pulse control sequences in and out of the satellite's orbital plane;

[0008] The on / off logic of the thruster is determined based on the pulse control sequences inside and outside the satellite orbital plane;

[0009] Based on the thruster's on / off logic, the thruster's on / off state is controlled to achieve multi-pulse high-precision autonomous position-keeping control for geostationary orbit satellites.

[0010] The present invention has the following advantages:

[0011] This invention discloses a multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites. It employs averaged orbital elements for position-keeping control and eliminates the effects of Earth's non-spherical gravity, lunar and solar gravitational forces, and solar radiation pressure perturbations. It controls only the long-term and long-period terms of orbital elements within one orbital period, improving orbital control accuracy and effectively avoiding the additional fuel consumption caused by short-period terms of orbital perturbations, thus significantly reducing fuel consumption and thruster exhaust frequency. The related algorithm is easily implemented in onboard code. This method can directly serve the high-precision autonomous position-keeping control of my country's high-resolution Earth imaging satellites and also provides valuable insights for other similar geostationary orbit satellites. Attached Figure Description

[0012] Figure 1 This is a flowchart illustrating the steps of a multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites according to an embodiment of the present invention.

[0013] Figure 2 This is a schematic diagram of the mathematical simulation verification results of a multi-pulse high-precision autonomous position-keeping control method for geostationary satellites in an embodiment of the present invention. Detailed Implementation

[0014] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments disclosed in the present invention will be described in further detail below with reference to the accompanying drawings.

[0015] One of the core ideas of this invention is to disclose a method for joint control of a thruster with a non-conical layout (the thruster nozzle is installed along the +X, -X, +Y, and -Y axes of the spacecraft; during operation, the thruster only generates thrust along the corresponding axes of the spacecraft and not along other axes). This method uses the spacecraft's average orbital elements as control variables, simultaneously performing drift rate, eccentricity vector, and phase control within the orbital plane, as well as inclination vector control outside the orbital plane, within each orbital period. By using average orbital elements for position-keeping control and eliminating the effects of Earth's non-spherical gravity, solar radiation pressure, and lunar gravitational perturbations, the method controls only the long-term and long-period terms of the orbital elements within one orbital period, improving orbital control accuracy and effectively avoiding the additional fuel consumption caused by short-period terms of orbital perturbations, thus effectively reducing fuel consumption and thruster exhaust frequency. The control algorithm in this method is directly solved algebraically, facilitating on-board code implementation.

[0016] Reference Figure 1 In this embodiment, the multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellite includes:

[0017] Step 1: Determine the average orbital elements of the satellite based on its instantaneous orbital elements.

[0018] In this embodiment, the specific implementation process of step 1 is as follows: Obtain the instantaneous orbital elements of the satellite: instantaneous semi-major axis a Osc Instantaneous eccentricity e Osc Instantaneous ascending node right ascension Ω Osc Instantaneous orbital inclination angle i Osc instantaneous perigee argument ω Osc and instantaneous approximate angle M Osc Based on the satellite's instantaneous orbital elements, the short-period terms of orbital element perturbations caused by J2, solar gravity, lunar gravity, and solar radiation pressure are calculated. Based on these factors, the satellite's average orbital elements are determined: the average semi-major axis a. Mean Mean eccentricity e Mean Mean ascending node right ascension Ω Mean Average orbital inclination i Mean Mean perigee argument ω Mean and mean aperimeter angle M Mean .

[0019] Preferably, the formula for calculating the average orbital elements of a satellite is as follows:

[0020] a Mean =a Osc -a SJ2 -a SS -a SL -a SSP

[0021] e Mean =e Osc -e SJ2 -e SS -e SL -e SSP

[0022] Ω Mean =Ω Osc -Ω SJ2 -Ω SS -Ω SL -Ω SSP

[0023] i Mean =i Osc -i SJ2 -i SS -i SL -i SSP

[0024] ωMean =ω Osc -ω SJ2 -ω SS -ω SL -ω SSP

[0025] M Mean =M Osc -M SJ2 -M SS -M SL -M SSP

[0026] Among them, a SJ2 e SJ2 i SJ2 Ω SJ2 ω SJ2 and M SJ2 Let a represent the short-period terms that cause perturbations in the semi-major axis, eccentricity, inclination, right ascension of the ascending node, argument of perigee, and mean perigee angle caused by J2; SS e SS i SS Ω SS ω SS and M SS Let a represent the short-period terms representing the perturbations of the orbit's semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and mean perigee caused by solar gravitational force; SL e SL i SL Ω SL ω SL and M SL Let a represent the short-period terms of the perturbations caused by lunar gravity on the semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and mean perigee angle; SSP e SSP i SSP Ω SSP ω SSP and M SSP These represent the short-period terms of the perturbations of the orbital semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and mean perigee angle caused by solar radiation pressure.

[0027] Preferably, the formula for calculating the short-period term of the orbital element perturbation caused by J2 is as follows:

[0028]

[0029]

[0030] Where J2 represents the second-order band harmonic coefficient of the Earth's gravitational potential; r Osc This represents the geocentric distance vector corresponding to the instantaneous orbital element. fOsc p represents the true anomaly angle corresponding to the instantaneous orbital element. Osc This represents the orbital semi-circle corresponding to the instantaneous orbital element.

[0031] Preferably, the formula for calculating the short-period term of the orbital element perturbation caused by solar gravity is as follows:

[0032]

[0033] Where, β S G represents the perturbation coefficient of solar gravity. S1 G represents the function related to the semi-major axis in the short-period terms of the orbital element perturbations caused by solar gravity. S2 G represents the function related to eccentricity in the short-period term of the orbital perturbation caused by solar gravity. S3 The first term G, representing the orbital inclination-related function in the short-period terms of the orbital element perturbation caused by solar gravity, is... S4 G represents the second term in the short-period term of the orbital element perturbation caused by solar gravity, which is related to the orbital inclination. S5 G represents the function related to the right ascension of the ascending node in the short-period term of the orbital perturbation caused by solar gravitational force. S6 G represents the function relating the perigee angle to the short-period term of the orbital perturbation caused by solar gravitational forces. S7 This represents a function relating the mean anomaly angle to the short-period terms of orbital element perturbations caused by solar gravity.

[0034]

[0035] G S3 =S2G S8 +S3G S9

[0036]

[0037] Where, m S The mass r of the sun S S1 represents the distance from the Sun to the central celestial body; S2 represents the first coefficient related to the angular relationship between the Sun and the central celestial body; S3 represents the third coefficient related to the angular relationship between the Sun and the central celestial body; E Osc G represents the anomalous angle corresponding to the instantaneous orbital elements of the satellite. S8 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S3 G S4 G S5 G S7 The first function used in the function procedure, G S9 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity.S3 G S4 G S5 G S7 The second function used in the function procedure, G S10 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S3 G S4 G S5 G S7 The third function used in the function procedure; A S B represents the first coefficient used in calculating the correlation coefficients S1, S2, and S3 between the sun and the central celestial body's angular term. S This represents the second coefficient used in calculating the correlation coefficients S1, S2, and S3 between the sun and the central celestial body's angular term; A S1 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S4 G S5 The first function of class A used in the function procedure, A S2 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S4 G S5 The second function of type A used in the function procedure; B S1 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S4 G S5 The first function of class B used in the function procedure, B S2 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S4 G S5 The second function of class B used in the function procedure.

[0038]

[0039] S3 = A S B S

[0040]

[0041]

[0042] Among them, i S θ represents the orbital inclination of the Sun relative to the central celestial body. SOsc θ represents the difference in right ascension between the satellite and the Sun relative to the central celestial body at the ascending node. SOsc =Ω Osc -Ω S Ω S Indicates the right ascension of the Sun relative to the ascending node of the central celestial body; u S u represents the latitudinal argument of the Sun relative to the central celestial body. S =fS +ω S ω S f represents the angle of perigee of the Sun relative to the central celestial body. S It represents the true angle of abduction of the Sun relative to the central celestial body.

[0043] Preferably, the formula for calculating the short-period term of the orbital element perturbation caused by lunar gravity is as follows:

[0044]

[0045] Where, β L G represents the perturbation coefficient of the moon's gravity. L1 G represents the function related to the semi-major axis in the short-period term of the orbital perturbation caused by lunar gravity. L2 G represents the function related to the eccentricity in the short-period term of the orbital perturbation caused by lunar gravity. L3 G represents the first term of the short-period term related to the orbital inclination in the perturbation of orbital elements caused by lunar gravity. L4 G represents the second term in the short-period term of the orbital inclination function related to the perturbation of orbital elements caused by lunar gravity. L5 G represents the function related to the right ascension of the ascending node in the short-period term of the orbital perturbation caused by lunar gravity. L6 G represents the function relating the perigee angle to the short-period term of the orbital perturbation caused by lunar gravity. L7 This represents a function relating the mean anomaly angle to the short-period term of the orbital element perturbation caused by lunar gravity.

[0046]

[0047]

[0048] G L3 =S2G L8 +S3G L9

[0049]

[0050] Where, m L The mass r of the moon L G represents the distance from the Moon to the central body; L1 represents the first coefficient related to the angular relationship between the Moon and the central body, L2 represents the second coefficient related to the angular relationship between the Moon and the central body, and L3 represents the third coefficient related to the angular relationship between the Moon and the central body; L8 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L3 G L4 G L5 G L7 The first function used in the function procedure, G L9G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L3 G L4 G L5 G L7 The second function used in the function procedure, G L10 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L3 G L4 G L5 G L7 The third function used in the function procedure; A L B represents the first coefficient used in calculating the correlation coefficients L1, L2, and L3 between the Moon and the central celestial body's angular term. L This refers to the second coefficient used in calculating the correlation coefficients L1, L2, and L3 between the Moon and the central celestial body; A L1 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L4 G L5 The first function of class A used in the function procedure, A L2 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L4 G L5 The second function of type A used in the function procedure; B L1 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L4 G L5 The first function of class B used in the function procedure, B L2 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L4 G L5 The second function of class B used in the function procedure.

[0051]

[0052] L3 = A L B L

[0053]

[0054]

[0055] Among them, i L θ represents the orbital inclination of the Moon relative to the central celestial body. LOsc θ represents the difference in right ascension between the satellite and the moon relative to the central celestial body at their ascending nodes. LOsc =Ω Osc -Ω L Ω L Indicates the right ascension of the Moon's ascending node relative to the central celestial body; u Lu represents the latitudinal argument of the Moon relative to the central celestial body. L =f L +ω L ω L f represents the angle of perigee of the Moon relative to the central celestial body. L It represents the true angle of anomaly of the Moon relative to the central celestial body.

[0056] Preferably, the formula for calculating the short-period term of the orbital element perturbation caused by solar radiation pressure is as follows:

[0057]

[0058]

[0059] Where, β SP The perturbation coefficient, A, represents the solar radiation pressure. SP B represents the first term of the function used to calculate the short-period terms of orbital element perturbations caused by solar radiation pressure. SP H1 represents the second term of the short-period function used to calculate the perturbation of orbital elements caused by solar radiation pressure, and H2 represents the first auxiliary variable and the second auxiliary variable.

[0060]

[0061]

[0062] H1 = sini Osc [(1+cosi S cos(ω) Osc -Ω Osc +u S )

[0063] +(1-cosi S cos(ω) Osc -Ω Osc -u S )]

[0064] -sini Osc [(1+cosi S cos(ω) Osc +Ω Osc -u S )

[0065] +(1-cosi S cos(ω) Osc +Ω Osc +u S )]

[0066] +2cosi Osc sini S [cos(ωOsc -u S )-cos(ω Osc +u S )]

[0067] H2 = sinini Osc [(1+cosi S sin(ω) Osc -Ω Osc +u S )

[0068] +(1-cosi S sin(ω) Osc -Ω Osc -u S )]

[0069] -sini Osc [(1+cosi S sin(ω) Osc +Ω Osc -u S )

[0070] +(1-cosi S sin(ω) Osc +Ω Osc +u S )]

[0071] +2cosi Osc sini S [sin(ω Osc -u S )-sin(ω Osc +u S )]

[0072] Where, k SP ρ represents the light pressure coefficient, S represents the illuminated area of ​​the satellite, and ρ represents the radiation pressure coefficient. S Indicates light pressure intensity, Δ S The distance from the sun to the satellite is represented by m, and the mass of the satellite is represented by m.

[0073] Step 2: Determine the pulse control sequences in and out of the satellite's orbital plane based on the satellite's average orbital features.

[0074] In this embodiment, the specific implementation process of step 2 is as follows: Based on the average orbital elements of the satellite, the orbital element control quantities are calculated: orbital drift rate control quantity ΔD, x-direction eccentricity vector control quantity Δe. x y-direction eccentricity vector control quantity Δe y x-direction orbital tilt vector control quantity Δi x y-direction orbital tilt vector control quantity Δi yAnd longitude control quantity Δλ; based on orbital element control quantities and the satellite's mean longitude l at the current moment. Mean The time t of the first in-plane pulse control was calculated. AEM1 The corresponding jet pulse quantity Δv AEM1 The second in-plane pulse control time t AEM2 The corresponding jet pulse quantity Δv AEM2 and the third in-plane pulse control time t AEM3 The corresponding jet pulse quantity Δv AEM3 ; and, the out-of-plane pulse control time t is calculated. I1 The corresponding jet pulse quantity Δv I1 According to Δv AEM1 Δv AEM2 Δv AEM3 and Δv I1 The pulse control sequences inside and outside the satellite orbital plane were constructed.

[0075] Preferably, ΔD, Δe x Δe y , Δi x , Δi y , Δλ and l Mean The calculation formula is as follows:

[0076]

[0077] Δe x =-e Mean cos(Ω Mean +ω Mean )

[0078] Δe y =-e Mean sin(Ω Mean +ω Mean )

[0079] Δi x =-i Mean cosΩ Mean

[0080] Δi y =-i Mean sinΩ Mean

[0081] Δλ=λ0-(Ω Mean +ω Mean +M Mean -α E )

[0082] l Mean =Ω Mean +ω Mean +MMean

[0083] Among them, a S ω represents the actual orbital radius of a geosynchronous satellite. E λ represents the Earth's rotational angular velocity, λ0 represents the longitude of the satellite's fixed point, and α... E This represents the sidereal hour angle corresponding to the current moment.

[0084] Preferably, based on the orbital control parameters and the satellite's current longitude... Mean The time t of the first in-plane pulse control was calculated. AEM1 The corresponding jet pulse quantity Δv AEM1 The second in-plane pulse control time t AEM2 The corresponding jet pulse quantity Δv AEM2 and the third in-plane pulse control time t AEM3 The corresponding jet pulse quantity Δv AEM3 ; and, the out-of-plane pulse control time t is calculated. I1 The corresponding jet pulse quantity Δv I1 ,include:

[0085] According to Δe x and Δe y The first mean longitude l for eccentricity control is calculated. ECtrl1 Second longitude l ECtrl2 :

[0086] l ECtrl1 =arctan2(Δe y ,Δe x )

[0087] l ECtrl2 =arctan2(-Δe y ,-Δe x )

[0088] Among them, l ECtrl1 and l ECtrl2 The domain of the definition is [0, 2π).

[0089] When the satellite's current longitude is l Mean Less than l ECtrlMin At that time, there were:

[0090]

[0091] Δv AEM1 =Δv xD +Δv xe +Δv xλ

[0092] Δv AEM2 =ΔvxD -Δv xe

[0093] Δv AEM3 =-Δv xλ

[0094] l AEM1 =l ECtrlMin

[0095] Where t represents the satellite's current time; ω E Indicates the Earth's angular velocity of rotation; l ECtrlMin For l ECtrl1 l ECtrl2 The smaller value; Δv xD This represents the pulse control quantity related to the orbital drift rate. v S Δv represents the actual velocity of a satellite in geosynchronous orbit. xe This represents the pulse control quantity related to the eccentricity vector; Δv xλ This represents the pulse control quantity related to longitude. l AEM1 This indicates the longitude at which the first in-plane pulse control is performed. When in l ECtrl1 When performing the first in-plane pulse control, When in l ECtrl2 When performing the first in-plane pulse control,

[0096] When the satellite's current longitude is l Mean Greater than l ECtrlMin And less than l ECtrlMax At that time, there were:

[0097]

[0098] Δv AEM1 =Δv xD +Δv xe +Δv xλ

[0099] Δv AEM2 =Δv xD -Δv xe

[0100] Δv AEM3 =-Δv xλ

[0101] l AEM1 =l ECtrlMax

[0102] Among them, l ECtrlMax For l ECtrl1 l ECtrl2 The larger value.

[0103] When the satellite's current longitude is l Mean Greater than l ECtrlMax At that time, there were:

[0104]

[0105]

[0106] Δv AEM1 =Δv xD +Δv xe +Δv xλ

[0107] Δv AEM2 =Δv xD -Δv xe

[0108] Δv AEM3 =-Δv xλ

[0109] l AEM1 =l ECtrlMin +2π

[0110] According to Δi x and Δi y The longitude l for tilt control is calculated. ICtrl :

[0111] l ICtrl =arctan2(Δi y ,Δi x )

[0112] When the satellite's current longitude is l Mean Less than l ICtrl At that time, there were:

[0113]

[0114] When the satellite's current longitude is l Mean Greater than l ICtrl At that time, there were:

[0115]

[0116] Step 3: Determine the thruster's on / off logic based on the pulse control sequences inside and outside the satellite's orbital plane.

[0117] In this embodiment, the specific implementation process of step 2 is as follows: Based on the pulse control sequences inside and outside the satellite orbital plane, the start time t of the first start-up operation of the satellite's X-axis thruster is calculated sequentially. Fire1Begin and deadline t Fire1EndThe start time t of the second activation of the satellite X-axis thruster Fire2Begin and deadline t Fire2End The start time t of the third activation of the satellite X-axis thruster Fire3Begin and deadline t Fire3End The start time t of the satellite-Y thruster's operation Fire4Begin and deadline t Fire4End According to t Fire1Begin t Fire1End t Fire2Begin t Fire2End t Fire3Begin t Fire3End t Fire4Begin and t Fire4End Determine the on / off logic of the thruster.

[0118]

[0119] Among them, F X F represents the magnitude of the satellite's thrust in the +X direction. InvX F represents the magnitude of the satellite's thrust in the X direction. InvY This indicates the magnitude of the satellite's thrust in the Y direction.

[0120] Step 4: Based on the thruster's on / off logic, control the thruster switch to achieve multi-pulse high-precision autonomous position-keeping control for geostationary orbit satellites.

[0121] like Figure 2 As shown in the simulation results, the east-west and north-south position control can achieve high accuracy using the method described in this invention.

[0122] In summary, this invention discloses a multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites. It employs numerical calculation methods, considering perturbations from Earth's non-spherical gravitational force, lunar and solar gravitational forces, and solar radiation pressure. Based on the satellite's instantaneous orbital elements, it calculates the average orbital elements, eliminating the influence of short-period terms of orbital perturbations while retaining the influence of long-term and long-period terms. A joint control method using drift rate, eccentricity, and phase is employed to formulate an east-west position-keeping control strategy within one orbital period, using average orbital elements to control the drift rate, eccentricity vector, orbital inclination vector, and phase error to near zero. An orbital inclination vector control method is used to formulate a north-south position-keeping control strategy within one orbital period. The method outputs the jetting time and duration for east-west and north-south position-keeping control.

[0123] Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make possible changes and modifications to the technical solutions of the present invention by utilizing the methods and techniques disclosed above without departing from the spirit and scope of the present invention. Therefore, any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solutions of the present invention shall fall within the protection scope of the technical solutions of the present invention.

[0124] The contents not described in detail in this specification are common knowledge to those skilled in the art.

Claims

1. A multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites, characterized in that, include: Determine the satellite's average orbital elements based on its instantaneous orbital elements; Based on the average orbital elements of the satellite, determine the pulse control sequences in and out of the satellite's orbital plane; based on the pulse control sequences in and out of the satellite's orbital plane, determine the thruster's on / off logic; based on the thruster's on / off logic, control the thruster's on / off state to achieve multi-pulse high-precision autonomous position-keeping control for geostationary orbit satellites; Based on the average orbital features of the satellite, determine the in-plane and out-of-plane pulse control sequences, including: calculating the orbital feature control quantities based on the average orbital features of the satellite: orbital drift rate control quantity ΔD, x-direction eccentricity vector control quantity Δe. x y-direction eccentricity vector control quantity Δe y x-direction orbital tilt vector control quantity Δi x y-direction orbital tilt vector control quantity Δi y And longitude control quantity Δλ; based on orbital element control quantities and the satellite's mean longitude l at the current moment. Mean The time t of the first in-plane pulse control was calculated. AEM1 The corresponding jet pulse quantity Δv AEM1 The second in-plane pulse control time t AEM2 The corresponding jet pulse quantity Δv AEM2 and the third in-plane pulse control time t AEM3 The corresponding jet pulse quantity Δv AEM3 ; and, the out-of-plane pulse control time t is calculated. I1 The corresponding jet pulse quantity Δv I1 According to Δv AEM1 Δv AEM2 Δv AEM3 and Δv I1 The pulse control sequences inside and outside the satellite orbital plane were constructed. Based on the orbital control parameters and the satellite's current longitude l Mean The time t of the first in-plane pulse control was calculated. AEM1 The corresponding jet pulse quantity Δv AEM1 The second in-plane pulse control time t AEM2 The corresponding jet pulse quantity Δv AEM2 and the third in-plane pulse control time t AEM3 The corresponding jet pulse quantity Δv AEM3 ; and, the out-of-plane pulse control time t is calculated. I1 The corresponding jet pulse quantity Δv I1 ,include: According to Δe x and Δe y The first mean longitude l for eccentricity control is calculated. ECtrl1 Second longitude l ECtrl2 :l ECtrl1 =arctan2(Δe y ,Δe x ), l ECtrl2 =arctan2(-Δe y ,-Δe x ); where l ECtrl1 and l ECtrl2 The domain is [0, 2π); When the satellite's current longitude is l Mean Less than l ECtrlMin At that time, there were: Δv AEM1 =Δv xD +Δv xe +Δv xλ Δv AEM2 =Δv xD -Δv xe ,Δv AEM3 =-Δv xλ ,l AEM1 =l ECtrlMin Where t represents the satellite's current time; ω E Indicates the Earth's angular velocity of rotation; l ECtrlMin For l ECtrl1 l ECtrl2 The smaller value; Δv xD The pulse control quantity, v, represents the orbital drift rate. S Δv represents the actual velocity of a satellite in geosynchronous orbit. xe This represents the pulse control quantity related to the eccentricity vector; Δv xλ The pulse control quantity related to longitude, l AEM1 Indicates the longitude at which the first in-plane pulse control is performed; When the satellite's current longitude is l Mean Greater than l ECtrlMin And less than l ECtrlMax At that time, there were: Δv AEM1 =Δv xD +Δv xe +Δv xλ Δv AEM2 =Δv xD -Δv xe ,Δv AEM3 =-Δv xλ ,l AEM1 =l ECtrlMax Among them, l ECtrlMax For l ECtrl1 l ECtrl2 The larger value; When the satellite's current longitude is l Mean Greater than l ECtrlMax At that time, there were: Δv AEM1 =Δv xD +Δv xe +Δv xλ Δv AEM2 =Δv xD -Δv xe ,Δv AEM3 =-Δv xλ ,l AEM1 =l ECtrlMin +2p According to Δi x and Δi y The longitude l for tilt control is calculated. ICtrl : l ICtrl =arctan2(Δi y ,Δi x ) When the satellite's current longitude is l Mean Less than l ICtrl At that time, there were: When the satellite's current longitude is l Mean Greater than l ICtrl At that time, there were:

2. The multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites according to claim 1, characterized in that, Based on the satellite's instantaneous orbital features, determine the satellite's average orbital features, including: Obtain the satellite's instantaneous orbital elements: instantaneous semi-major axis a Osc Instantaneous eccentricity e Osc Instantaneous ascending node right ascension Ω Osc Instantaneous orbital inclination angle i Osc instantaneous perigee argument ω Osc and instantaneous approximate angle M Osc ; Based on the satellite's instantaneous orbital elements, the short-period terms of orbital element perturbations caused by J2, the short-period terms of orbital element perturbations caused by solar gravity, the short-period terms of orbital element perturbations caused by lunar gravity, and the short-period terms of orbital element perturbations caused by solar radiation pressure were calculated. Based on the satellite's instantaneous orbital elements, the short-period terms of orbital element perturbations caused by J2, the short-period terms of orbital element perturbations caused by solar gravity, the short-period terms of orbital element perturbations caused by lunar gravity, and the short-period terms of orbital element perturbations caused by solar radiation pressure, the satellite's average orbital elements are determined: the average semi-major axis a Mean Mean eccentricity e Mean Mean ascending node right ascension Ω Mean Average orbital inclination i Mean Mean perigee argument ω Mean and mean aperimeter angle M Mean .

3. The multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites according to claim 2, characterized in that, The formula for calculating the average orbital elements of a satellite is as follows: a Mean =a Osc -a SJ2 -a SS -a SL -a SSP And Mean =and Osc -And SJ2 -And SS -And SL -And SSP Oh Mean =Oh Osc -Oh SJ2 -Oh SS -Oh SL -Oh SSP i Mean =i Osc -i SJ2 -i SS -i SL -i SSP oh Mean =ω Osc -oh SJ2 -oh SS -oh SL -oh SSP M Mean =M Osc -M SJ2 -M SS -M SL -M SSP Among them, a SJ2 e SJ2 i SJ2 Ω SJ2 ω SJ2 and M SJ2 Let a represent the short-period terms that cause perturbations in the semi-major axis, eccentricity, inclination, right ascension of the ascending node, argument of perigee, and mean perigee angle caused by J2; SS e SS i SS Ω SS ω SS and M SS Let a represent the short-period terms representing the perturbations of the orbit's semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and mean perigee caused by solar gravitational force; SL e SL i SL Ω SL ω SL and M SL Let a represent the short-period terms of the perturbations caused by lunar gravity on the semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and mean perigee angle; SSP e SSP i SSP Ω SSP ω SSP and M SSP These represent the short-period terms of the perturbations of the orbital semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and mean perigee angle caused by solar radiation pressure.

4. The multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites according to claim 3, characterized in that, The formula for calculating the short-period term of the orbital element perturbation caused by J2 is as follows: Where J2 represents the second-order band harmonic coefficient of the Earth's gravitational potential; r Osc This represents the geocentric distance vector corresponding to the instantaneous orbital element. f Osc p represents the true anomaly angle corresponding to the instantaneous orbital element. Osc This represents the orbital semi-circle corresponding to the instantaneous orbital element. The formula for calculating the short-period term of the orbital element perturbation caused by solar gravity is as follows: Where, β S G represents the perturbation coefficient of solar gravity. S1 G represents the function related to the semi-major axis in the short-period terms of the orbital element perturbations caused by solar gravity. S2 G represents the function related to eccentricity in the short-period term of the orbital perturbation caused by solar gravity. S3 The first term G, representing the orbital inclination-related function in the short-period terms of the orbital element perturbation caused by solar gravity, is... S4 G represents the second term in the short-period term of the orbital element perturbation caused by solar gravity, which is related to the orbital inclination. S5 G represents the function related to the right ascension of the ascending node in the short-period term of the orbital perturbation caused by solar gravitational force. S6 G represents the function relating the perigee angle to the short-period term of the orbital perturbation caused by solar gravitational forces. S7 A function representing the short-period term related to the mean anomaly angle in the orbital element perturbation caused by solar gravity; G S3 =S2G S8 +S3G S9 Where, m S The mass r of the sun S S1 represents the distance from the Sun to the central celestial body; S2 represents the first coefficient related to the angular relationship between the Sun and the central celestial body; S3 represents the third coefficient related to the angular relationship between the Sun and the central celestial body; E Osc G represents the anomalous angle corresponding to the instantaneous orbital elements of a satellite. S8 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S3 G S4 G S5 G S7 The first function used in the function procedure, G S9 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S3 G S4 G S5 G S7 The second function used in the function procedure, G S10 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S3 G S4 G S5 G S7 The third function used in the function procedure; A S B represents the first coefficient used in calculating the correlation coefficients S1, S2, and S3 between the sun and the central celestial body's angular term. S This represents the second coefficient used in calculating the correlation coefficients S1, S2, and S3 between the sun and the central celestial body's angular term; A S1 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S4 G S5 The first function of class A used in the function procedure, A S2 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S4 G S5 The second term function of type A used in the function procedure; B S1 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S4 G S5 The first function of class B used in the function procedure, B S2 G represents the short-period term in the calculation of orbital element perturbations caused by solar gravity. S4 G S5 The second function of class B used in the function procedure; S3=A S B S Among them, i S θ represents the orbital inclination of the Sun relative to the central celestial body. SOsc θ represents the difference in right ascension between the satellite and the Sun relative to the central celestial body at the ascending node. SOsc =Ω Osc -Ω S Ω S Indicates the right ascension of the Sun relative to the ascending node of the central celestial body; u S u represents the latitudinal argument of the Sun relative to the central celestial body. S =f S +ω S ω S f represents the angle of perigee of the Sun relative to the central celestial body. S This represents the true angle of abduction of the Sun relative to the central celestial body. The formula for calculating the short-period term of the orbital element perturbation caused by lunar gravity is as follows: Where, β L G represents the perturbation coefficient of the moon's gravity. L1 G represents the function related to the semi-major axis in the short-period term of the orbital element perturbation caused by lunar gravity. L2 G represents the function related to eccentricity in the short-period term of the orbital perturbation caused by lunar gravity. L3 G represents the first term of the short-period term related to the orbital inclination in the perturbation of orbital elements caused by lunar gravity. L4 G represents the second term in the short-period term of the orbital inclination function related to the perturbation of orbital elements caused by lunar gravity. L5 G represents the function related to the right ascension of the ascending node in the short-period term of the orbital perturbation caused by lunar gravity. L6 G represents the function relating the perigee angle to the short-period term of the orbital perturbation caused by lunar gravity. L7 A function relating the mean anomaly angle in the short-period terms of orbital element perturbations caused by lunar gravity; G L3 =S2G L8 +S3G L9 Where, m L The mass r of the moon L G represents the distance from the Moon to the central body; L1 represents the first coefficient related to the angular relationship between the Moon and the central body, L2 represents the second coefficient related to the angular relationship between the Moon and the central body, and L3 represents the third coefficient related to the angular relationship between the Moon and the central body; L8 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L3 G L4 G L5 G L7 The first function used in the function procedure, G L9 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L3 G L4 G L5 G L7 The second function used in the function procedure, G L10 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L3 G L4 G L5 G L7 The third function used in the function procedure; A L B represents the first coefficient used in calculating the correlation coefficients L1, L2, and L3 between the Moon and the central celestial body's angular term. L This refers to the second coefficient used in calculating the correlation coefficients L1, L2, and L3 between the Moon and the central celestial body; A L1 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L4 G L5 The first function of class A used in the function procedure, A L2 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L4 G L5 The second term function of type A used in the function procedure; B L1 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L4 G L5 The first function of class B used in the function procedure, B L2 G represents the short-period term in the calculation of orbital element perturbations caused by lunar gravity. L4 G L5 The second function of class B used in the function procedure; L3=A L B L Among them, i L θ represents the orbital inclination of the Moon relative to the central celestial body. LOsc θ represents the difference in right ascension between the satellite and the moon relative to the central celestial body at their ascending nodes. LOsc =Ω Osc -Ω L Ω L Indicates the right ascension of the Moon's ascending node relative to the central celestial body; u L u represents the latitudinal argument of the Moon relative to the central celestial body. L =f L +ω L ω L f represents the angle of perigee of the Moon relative to the central celestial body. L This represents the true anomaly of the Moon relative to the central celestial body; The formula for calculating the short-period term of the orbital element perturbation caused by solar radiation pressure is as follows: Where, β SP The perturbation coefficient, A, represents the solar radiation pressure. SP B represents the first term of the function used to calculate the short-period terms of orbital element perturbations caused by solar radiation pressure. SP H1 represents the second term of the short-period function used to calculate the perturbation of orbital elements caused by solar radiation pressure, and H2 represents the first auxiliary variable and the second auxiliary variable. H1=sini Osc [(1+so S )cos(ω Osc -Ω Osc +u S ) +(1-so S )cos(ω Osc -Ω Osc -u S )] -sini Osc [(1+so S )cos(ω Osc +Ω Osc -u S ) +(1-so S )cos(ω Osc +Ω Osc +u S )] +2 things Osc son S [cos(ω Osc -in S )-cos(ω Osc +in S )] H2=sini Osc [(1+cos S )sin(ω Osc -Oh Osc +u S ) +(1-cos S )sin(ω Osc -Oh Osc -u S )] -sini Osc [(1+cos S )sin(ω Osc +Oh Osc -u S ) +(1-cos S )sin(ω Osc +Oh Osc +u S )] +2so Osc this S [sin(ω Osc -u S )-sin(ω Osc +u S )] Where, k SP ρ represents the light pressure coefficient, S represents the illuminated area of ​​the satellite, and ρ represents the radiation pressure coefficient. S Indicates light pressure intensity, Δ S The distance from the sun to the satellite is represented by m, and the mass of the satellite is represented by m.

5. The multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites according to claim 4, characterized in that, ΔD, Δe x Δe y , Δi x , Δi y , Δλ and l Mean The calculation formula is as follows: No x =-e Mean cos(Ω Mean +oh Mean ) No y =-e Mean sin(Ω Mean +oh Mean ) Δi x =-i Mean cosΩ Mean Δi y =-i Mean sinΩ Mean Δλ=λ0-(Ω Mean +oh Mean +M Mean -a E ) l Mean =Oh Mean +oh Mean +M Mean Among them, a S ω represents the actual orbital radius of a geostationary satellite. E λ represents the Earth's rotational angular velocity, λ0 represents the longitude of the satellite's fixed point, and α... E This represents the sidereal hour angle corresponding to the current moment.

6. The multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites according to claim 5, characterized in that, When in l ECtrl1 When performing the first in-plane pulse control, When in l ECtrl2 When performing the first in-plane pulse control, 7. The multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites according to claim 6, characterized in that, Based on the pulse control sequences inside and outside the satellite orbital plane, the thruster's on / off logic is determined, including: Based on the pulse control sequences inside and outside the satellite's orbital plane, the start time t of the satellite's X-axis thruster's first activation was calculated sequentially. Fire1Begin and deadline t Fire1End The start time t of the second activation of the satellite X-axis thruster Fire2Begin and deadline t Fire2End The start time t of the third activation of the satellite X-axis thruster Fire3Begin and deadline t Fire3End The start time t of the satellite-Y thruster's operation Fire4Begin and deadline t Fire4End ; According to t Fire1Begin t Fire1End t Fire2Begin t Fire2End t Fire3Begin t Fire3End t Fire4Begin and t Fire4End Determine the on / off logic of the thruster.

8. The multi-pulse high-precision autonomous position-keeping control method for geostationary orbit satellites according to claim 7, characterized in that, Among them, F X F represents the magnitude of the satellite's thrust in the +X direction. InvX F represents the magnitude of the satellite's thrust in the X direction. InvY This indicates the magnitude of the satellite's thrust in the Y direction.

Images