A method for maintaining the position of a geosynchronous orbit satellite using omnidirectional thrust coupling

Abstract

This invention belongs to the field of spacecraft control and relates to an omnidirectional thrust-coupled position-holding method for geostationary orbit satellites. The invention employs a multi-pulse control scheme considering thrust coupling, applicable to situations where north-south and east-west position-holding have radial thrust coupling effects. The invention includes processes such as calculating north-south position-holding control quantities, eccentricity control quantities, radial thrust compensation, iterative calculation of longitude drift rate control quantities, calculation of two-pulse control parameters, and calculation of three-pulse control parameters. Through optimization of position-holding ignition timing and iterative solution of control quantities, joint thrust-coupled position-holding control is achieved.

Description

Technical Field

[0001] This invention relates to a method for maintaining the position of a geosynchronous orbit satellite using omnidirectional thrust coupling, which belongs to the field of spacecraft control. Background Technology

[0002] Because the new generation geostationary orbit satellite platform (hereinafter referred to as the "new platform") adopts a modular design concept, separating the payload module and the propulsion service module, the thrusters can only be arranged in the propulsion service module located on one side of the satellite-rocket docking ring. This means that the installation height of all thrusters is lower than the center of mass height of the satellite during its on-orbit phase. Due to this characteristic, compared with the 10N thruster layout scheme of traditional geostationary orbit satellite platforms, the 10N thrusters of the new platform always have a radial thrust component during position maintenance. This causes the north-south tilt control to be coupled with the east-west eccentricity and longitude control, bringing new challenges to the position maintenance strategy design of the 10N thrusters.

[0003] In traditional geostationary orbit satellite 10N thruster position-keeping strategies, the coupling between north-south and east-west position-keeping is mainly caused by control errors and is generally treated as a disturbance. No specific strategy optimization is needed; a post-control, pre-control orbit measurement and real-time independent control strategy suffices to meet mission requirements. However, the new platform itself features a high degree of coupling between north-south and east-west position-keeping thrust control, requiring the control strategy to be optimized to address this coupling factor and ensure fuel efficiency. Summary of the Invention

[0004] The technical problem solved by this invention is to provide an omnidirectional thrust-coupled position-keeping method for geostationary orbit satellites. Through control quantity calculation, radial thrust compensation, position-keeping timing optimization, and iterative solution of control quantities, joint control for position-keeping under thrust coupling is achieved.

[0005] The technical solution of this invention is: a method for maintaining the omnidirectional thrust coupling position of a geosynchronous orbit satellite, comprising the following steps:

[0006] (1) Calculate the orbital angular velocity n and mean longitude based on the orbital parameters at the current moment. Longitude deviation Longitude drift rate D, eccentricity vector Inclination vector And perform control calculation and judgment. When it is judged to be a waiting track, suspend control; when it is judged to perform north-south position-protection control calculation, go to step (2); when it is judged to perform east-west position-protection control calculation, go to step (3).

[0007] (2) According to the Julian century number T in J2000 C Calculate the real-time Julian century number T NCJ2000 tilt angle variation coefficient T Ci Real-time tilt angle change coefficient T NCi Then calculate the long-term term of the tilt perturbation [di] x ,di y ];

[0008] Based on the above calculation of the tilt perturbation long term [di] x ,di y ] Calculate the ignition reference right ascension L F And tilt angle control quantity Δi;

[0009] The north-south positioning system employs a two-pulse control, setting the first ignition's right ascension to L. F Thus, the ignition times t of the two pulses with north and south position-preserving characteristics are obtained. ns1 t ns2 ;

[0010] (3) Calculate the eccentricity control quantity Δe x Δe y and drift rate control quantity And calculate the drift rate control speed increment ΔV based on the calculation results. d Eccentricity control speed increment ΔV e ;

[0011] Based on the allowable width Δλ of the longitude drift ring and the drift rate control amount Drift rate control speed increment ΔV d Determine the pulse control type and the control quantity for each pulse control type.

[0012] Preferably, the control calculation and judgment are performed in the following manner:

[0013] Every preset N C Calculations for north-south positional protection control were performed.

[0014] if or or Then perform east-west position-protection control calculations; otherwise, it is a waiting track and will not be controlled for the time being.

[0015] The aforementioned Δλ represents the width of the allowable longitude drift ring. For longitude drift acceleration, e f To control the target radius using eccentricity, the function norm(.) is used to calculate the magnitude of a two-dimensional vector in a plane.

[0016] Preferably, the two ignition times t of the north-south position-preserving dual-pulse system. ns1 t ns2 The calculation formula is as follows:

[0017]

[0018]

[0019] Where k is an integer, and its value should be such that the actual first ignition time is after the current time and is the closest; ω is the perigee argument; Ω is the right ascension of the ascending node; t p This refers to the moment when the satellite passes through its closest point to Earth.

[0020] Preferably, the pulse control type is determined in the following manner:

[0021] According to the deviation of longitude Dual-pulse drift rate control quantity Drift rate control speed increment ΔV d The nominal radius r of the synchronous orbit s Calculate the pulse control threshold value for longitude.

[0022] like If two pulses are used, then two pulses control is used; otherwise, three pulses control is used. This can be categorized into two types:

[0023] like That is, for GEO satellites that drift westward due to natural perturbations:

[0024] like If the left boundary is used, then three-pulse control is used; otherwise, three-pulse control is used.

[0025] like That is, for GEO satellites that drift eastward due to natural perturbations.

[0026] like If the right boundary three-pulse control is used, then the left boundary three-pulse control is used; otherwise, the left boundary three-pulse control is used.

[0027] Preferably, the calculation steps for the control quantity of the left boundary three-pulse control are as follows:

[0028] Calculate the target longitude difference:

[0029] Calculate the target drift rate:

[0030] Calculate the drift rate control value:

[0031] The calculation takes into account the radial component causing longitude drift:

[0032] Calculate the change in longitude:

[0033] Then, the three-pulse control ignition parameters are determined: that is, the three ignition times are determined, and the east-west position-maintaining three-pulse control speed increment is calculated based on the calculated change in longitude difference.

[0034] Preferably, the calculation steps for the control quantity of the right boundary three-pulse control are as follows:

[0035] Calculate the target longitude difference:

[0036] Calculate the target drift rate: D f =0

[0037] Calculate the drift rate control value:

[0038] The calculation takes into account the radial component causing longitude drift:

[0039]

[0040] Calculate the change in longitude:

[0041] Then, the three-pulse control ignition parameters are determined: that is, the three ignition times are determined, and the east-west position-maintaining three-pulse control speed increment is calculated based on the calculated change in longitude difference.

[0042] Preferably, the speed increment ΔV is controlled during the east-west position-preserving control calculation. k Limitation, i.e., when When, take ΔV k =0;

[0043] Where abs(.) represents taking the absolute value, This is the minimum speed increment limit.

[0044] The beneficial effects of this invention compared to the prior art are:

[0045] This invention proposes an omnidirectional thrust-coupled position-keeping method for geostationary orbit satellites. Considering the characteristics of the new platform, it proposes a coordinated position-keeping control strategy involving north-south and east-west control. The position-keeping timing adopts a "step-by-step execution, unified calculation" approach, meaning that the ignition for north-south and east-west position-keeping is executed separately. The impact of north-south position-keeping ignition on east-west position-keeping is considered, and coordinated control is achieved through control quantity compensation and control timing optimization.

[0046] The north-south position maintainer operates in two separate 180-degree intervals at the ascending and descending nodes to eliminate eccentricity interference caused by the radial thrust component. The east-west position maintainer employs combined longitude and eccentricity control, incorporating different strategies such as single-pulse, double-pulse, or triple-pulse control. This corrects for the radial component effects of the east-west position maintainer while simultaneously controlling longitude drift and eccentricity vector.

[0047] The omnidirectional thrust-coupled position-keeping method for geostationary orbit satellites proposed in this invention is applicable to various north-south and east-west position-keeping scenarios with radial thrust coupling. Through position-keeping timing optimization and iterative solution of the longitude control variables, thrust coupling compensation is achieved, enabling joint control under dual-pulse and triple-pulse operating conditions. Compared with existing technologies, it has the following advantages:

[0048] (1) A uniform interval control strategy for north-south position maintenance was proposed, and based on this, an optimization strategy for the timing selection of north-south position maintenance was designed, thereby extending the east-west position maintenance interval and avoiding increased position maintenance fuel consumption due to control coupling.

[0049] (2) An iterative calculation process for the east-west longitude control quantity considering radial thrust coupling is proposed. Through method design, the iterative calculation is made convergent and can adapt to dual-pulse and triple-pulse operating conditions. This method achieves thrust coupling compensation through a finite number of iterations and has the advantages of small computational load and good convergence. Attached Figure Description

[0050] Figure 1 This is a flowchart of the calculation process;

[0051] Figure 2 Control effect after adding eccentricity control at fixed position 120°E

[0052] Figure 3 Control effect after adding eccentricity control at fixed position 160°E

[0053] Figure 4 To improve the track tilt control effect after relative position closed-loop control Detailed Implementation

[0054] This invention provides a method for maintaining the position of a geosynchronous orbit satellite using omnidirectional thrust coupling, achieving onboard autonomous control performance indicators while meeting engineering constraints such as fuel consumption, software computation, and energy limitations.

[0055] A method for maintaining the position of a geostationary orbit satellite using omnidirectional thrust coupling, comprising the following steps under certain computational conditions:

[0056] (1) Input parameter calculation

[0057] Input the calculation parameters for the current time t. The meanings and some values ​​of each parameter are as follows:

[0058] a semi-major axis

[0059] e Eccentricity

[0060] i. Track Inclination

[0061] ω Perigee Argument

[0062] Ω Right Ascension of Ascending Node

[0063] M (approximate point angle)

[0064] tp last perigee time

[0065] SG Greenwich Mean Time

[0066] us is the solar direction angle measured from the vernal equinox.

[0067] is the obliquity of the ecliptic, and its value is...

[0068] r s The nominal radius of the geostationary orbit is 42165.7 km.

[0069] V s The geosynchronous orbit velocity is 3.0746 km / s.

[0070] μ is the Earth's gravitational constant, with values ​​ranging from...

[0071] Td East-West Position-Maintaining Drift Period

[0072] λ f Longitude of fixed target

[0073] e f Eccentricity control target radius

[0074] θ e Eccentricity vector control target angle

[0075] Longitude drift acceleration

[0076] Δλ is the width of the allowable longitude drift ring.

[0077] N C The north-south positioning interval is 7 days by default.

[0078] Calculate other parameters based on the input parameters, using the following formula:

[0079] Calculate the orbital angular velocity n:

[0080] Calculate the mean longitude:

[0081] Calculate the deviation of longitude:

[0082] Calculate the longitude drift rate: D = -3π(ar) s ) / r s

[0083] Calculate the eccentricity vector Two of the components are:

[0084] e x =e*cos(ω+Ω),e y = e*sin(ω+Ω)

[0085] Calculate the tilt vector Two of the components are:

[0086] i x =i*cos(Ω),i y =i*sin(Ω)

[0087] (2) Control calculation judgment:

[0088] Every N C The system calls step (3) once to perform north-south positional protection control calculations.

[0089] if or or Then step (4) is called once to perform the east-west position protection control calculation;

[0090] otherwise

[0091] It is a waiting track and will not be controlled for the time being.

[0092] (3) Calculation of North-South Position Preservation Control

[0093] Enter the Julian century number relative to 12:00:00 on January 1, 2000, i.e., J2000 Julian century number T. C ;

[0094] Calculate the real-time Julian century number T NC J2000 tilt angle variation coefficient T Ci Real-time tilt angle change coefficient T NCi :

[0095] T NC =T C +N C / 36525.0

[0096] T Ci =-2.182439197+33.75704461*T C

[0097] T NCi =-2.182439197+33.75704461*T NC

[0098] Calculate the long-term term of tilt perturbation [di] x ,di y ]:

[0099] di x =36525.0*[0.106938839*10 -4 *(cos(T Ci )-cos(T NCi ))-0.2081091193*10- 6 *(cos 2 (T Ci )-cos 2 (T NCi ))]

[0100] di y =36525.0*[-0.2329662483*10 -2 *(TT NC -0.7967206128*10 -5 *(sin(T Ci )-sin(T NCi ))+0.1909363838*10 -6 *(cos(T Ci )*sin(T Ci )-cos(T NCi )*sin(T NCi ))]

[0101] Calculate the ignition reference right ascension L F And tilt control quantity Δi:

[0102]

[0103]

[0104] Δi=-norm([Δi x ,Δi y ])

[0105] L F =arg([-Δi x ,-Δi y ])

[0106] Where Δi x , Δi y These are the X and Y components of the tilt control quantity in the inertial frame, respectively. x di y These are the X and Y components of the long-term term of the tilt perturbation in the inertial frame, respectively, K. i The tilt control coefficient is set to 0.15; i x i yLet X and Y be the tilt vector in the inertial frame at the current moment. The function norm(.) calculates the magnitude of a two-dimensional planar vector, and arg(.) calculates the argument of a two-dimensional planar vector.

[0107] The north and south position maintainers use two pulse controls, with a half-day interval between the two controls.

[0108] Let the first ignition be of right ascension L. F Thus, the ignition times t of the two pulses with north and south position-preserving characteristics are obtained. ns1 t ns2 They are respectively:

[0109]

[0110]

[0111] Where k is an integer, and its value should be such that the actual first ignition time is after the current time and is the closest; ω is the perigee argument; Ω is the right ascension of the ascending node; t p This refers to the moment when the satellite passes through its closest point to Earth.

[0112] (4) Calculation of East-West Position Preservation Control

[0113] (4-1) Calculate the eccentricity control quantity

[0114]

[0115]

[0116] Where η is the eccentricity control bias coefficient.

[0117] (4-2) Calculate the double-pulse control quantity:

[0118] Let Δλ be the longitudinal drift caused by radial disturbance. R Take Δλ R The initial value is 0. The following three steps (a) to (c) are iteratively calculated 10 times, and the drift rate control value is output.

[0119]

[0120]

[0121]

[0122] Where σ is the longitude control offset coefficient and β is the thruster radial installation angle.

[0123] Then calculate the drift rate control speed increment ΔV. d Eccentricity control speed increment ΔV e :

[0124]

[0125]

[0126] (4-3) Determine the pulse control type:

[0127] Calculate the longitude pulse control threshold

[0128] like Then dual-pulse control is used, rotating (4-4);

[0129] Otherwise, three-pulse control is used, which falls into two categories:

[0130] like That is, for GEO satellites that drift westward due to natural perturbations:

[0131] like Then, the left boundary three-pulse control is adopted, and the rotation is (4-5);

[0132] Otherwise, use right boundary three-pulse control and turn (4-6).

[0133] like That is, for GEO satellites that drift eastward due to natural perturbations.

[0134] like Then, the right boundary three-pulse control is adopted, and the rotation is (4-6);

[0135] Otherwise, use left boundary three-pulse control and turn (4-5).

[0136] (4-4) Dual-pulse control ignition parameters

[0137] Dual-pulse control is used, with a half-day interval between the two control cycles.

[0138] The first ignition at right ascension l1 is:

[0139] l1=arg(±[Δe x Δe y ])

[0140] This yields the two ignition times t of the east-west position-preserving double pulse. ew1 t ew2 They are, in order:

[0141]

[0142]

[0143] Where k takes an appropriate integer, and the choice of the ± sign in the l1 calculation formula, as well as the value of k, should ensure that the actual first ignition time is after the current time and is the closest. Then, determine the ± in the following formula according to the correspondence. Symbols, calculating the two velocity increments ΔV in the east-west position-preserving double-pulse control. 21 ΔV 22 They are, in order:

[0144]

[0145]

[0146] Proceed to step (5).

[0147] (4-5) Calculation of the three-pulse control quantity at the left boundary

[0148] Calculate the target longitude difference:

[0149] Calculate the target drift rate:

[0150] Calculate the drift rate control value:

[0151] The calculation takes into account the radial component causing longitude drift:

[0152]

[0153] Calculate the change in longitude difference of the three pulses: Turn (4-7);

[0154] (4-6) Calculation of the three-pulse control quantity at the right boundary

[0155] Calculate the target longitude difference:

[0156] Calculate the target drift rate: D f =0

[0157] Calculate the drift rate control value:

[0158] The calculation takes into account the radial component causing longitude drift:

[0159]

[0160] Calculate the change in longitude difference of the three pulses: Turn (4-7);

[0161] (4-7) Three-pulse control ignition parameters

[0162] The first ignition is equal to the right ascension:

[0163] l1=arg(±[Δe x Δe y ])

[0164] This yields the three ignition times t. ew1 t ew2 t ew3 They are, in order:

[0165]

[0166]

[0167] Where k takes an appropriate integer, and the choice of ± sign in the l1 calculation formula and the value of k should ensure that the actual first ignition time is after the current time and is the closest.

[0168] Then, determine ± in the following formula according to the correspondence. Symbol, calculate the velocity increment ΔV of the three pulse control for east-west position maintenance. 31 ΔV 32 ΔV 33 They are, in order:

[0169]

[0170]

[0171]

[0172] In the formula, nt = M + n·t p average angular velocity

[0173] Proceed to step (5);

[0174] (5) Control speed increment limit

[0175] When ΔV k When (k = 21, 22 or k = 31, 32, 33) is very small, this control is not performed. That is: when When ΔVk = 0, take ΔVk = 0.

[0176] Where abs(.) represents taking the absolute value, The minimum speed increment limit is set to 0.01 m / s.

[0177] (6) Calculation ends

[0178] Example

[0179] The parameters set in the calculation example are shown in Table 1:

[0180] Table 1 Example Setting Parameters

[0181]

[0182] Table 1 Initial orbital elements

[0183]

[0184] Simulation results show that after adding eccentricity control, the satellite eccentricity remains within 0.0002 throughout the year, resulting in minimal longitude diurnal fluctuations and maintaining the satellite's east-west positional accuracy within ±0.05°. Table 3 shows the total velocity increments required for both longitude and eccentricity control throughout the year, and the control process curves are shown in [Figure 1]. Figure 2 and Figure 3 .

[0185] Table 3. Statistics on the Increment of Control Speed ​​Based on Longitude and Eccentricity

[0186] Add eccentricity control at fixed position 120°E 2.562 Add eccentricity control at fixed position 160°E 2.400

[0187] Simulation results show that after increasing the tilt angle closed-loop control variable, the satellite orbital tilt angle remains within 0.01° throughout the year, and the tilt angle vector drift loop stabilizes within a very small range. The statistics of the tilt angle control velocity increment over the year are shown in Table 4, and the control process curves are shown in... Figure 4 .

[0188] Table 4. Statistics on Speed ​​Increment for Tilt Angle Control

[0189] Perturbation control quantity + Closed-loop control quantity 49.464

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

Claims

1. A method for maintaining the position of a geosynchronous orbit satellite using omnidirectional thrust coupling, characterized in that... Includes the following steps: (1) Calculate the orbital angular velocity n and mean longitude based on the orbital parameters at the current moment. Longitude deviation Longitude drift rate D, eccentricity vector Inclination vector The system performs control calculations and makes judgments. When the system is judged to be a waiting track, the system pauses control. When the system is judged to be performing north-south position-keeping control calculations, the system proceeds to step (2). When the system is judged to be performing east-west position-keeping control calculations, the system proceeds to step (3). (2) According to the Julian century number J2000 Calculate real-time Julian centuries J2000 tilt angle variation coefficient Real-time tilt angle change coefficient Then calculate the long-term term of the tilt perturbation. , ]; Based on the above calculation of the tilt perturbation long term [ , ] Calculate the reference right ascension for ignition and tilt control amount ; The north-south positioning system employs a two-pulse control mechanism, ensuring that the first ignition is at the same right ascension. This yields the ignition times of the two dual-pulse firing cycles in the north and south positions. , ; (3) Calculate the eccentricity control value , and drift rate control quantity And calculate the drift rate control speed increment based on the calculation results. Eccentricity control speed increment ; Based on the width of the allowed longitude drift ring Drift rate control quantity Drift rate control speed increment Determine the pulse control type and the control quantity for each pulse control type; The control calculation and judgment are performed in the following manner: Every preset Calculations for north-south positional protection control were performed. if or or If the trajectory is positive, then east-west position-protection control calculations will be performed; otherwise, it is a waiting track and will not be controlled for the time being. The above For the width of the allowable longitude drift ring, For longitude drift acceleration, To control the target radius using eccentricity, the function This represents finding the magnitude of a two-dimensional vector in a plane. North and South Position-Preserving Dual-Pulse Ignition Timing , The calculation formula is as follows: Where k is an integer, and its value should be such that the actual first ignition time is after the current time and is the closest; ω is the perigee argument; Ω is the right ascension of the ascending node; The time of the satellite's closest approach to Earth; The pulse control type is determined as follows: According to the deviation of longitude Dual-pulse drift rate control quantity Drift rate control speed increment nominal radius of synchronous orbit Calculate the pulse control threshold value for longitude. ; like If the condition is met, then dual-pulse control is used; otherwise, three-pulse control is used. This is divided into two categories: like That is, for GEO satellites that drift westward due to natural perturbations: like If the left boundary three-pulse control is used, then the right boundary three-pulse control is used; otherwise, the left boundary three-pulse control is used. like That is, for GEO satellites that drift eastward due to natural perturbations, like If the condition is met, then the right boundary three-pulse control is used; otherwise, the left boundary three-pulse control is used.

2. The method according to claim 1, characterized in that: The calculation steps for the control quantity of the left boundary three-pulse control are as follows: Calculate the target longitude difference: Calculate the target drift rate: Calculate the drift rate control value: The calculation takes into account the radial component causing longitude drift: Calculate the change in longitude: Then, the three-pulse control ignition parameters are determined: that is, the three ignition times are determined, and the east-west position-maintaining three-pulse control speed increment is calculated based on the calculated change in longitude difference.

3. The method according to claim 1, characterized in that: The calculation steps for the control quantity of the right boundary three-pulse control are as follows: Calculate the target longitude difference: Calculate the target drift rate: Calculate the drift rate control value: The calculation takes into account the radial component causing longitude drift: Calculate the change in longitude: Then, the three-pulse control ignition parameters are determined: that is, the three ignition times are determined, and the east-west position-maintaining three-pulse control speed increment is calculated based on the calculated change in longitude difference.

4. The method according to claim 1, characterized in that: During the east-west position-preserving control calculation, the control velocity increment ΔV is calculated. k Limiting, i.e., when abs(ΔV) k )< When, take ΔV k =0; Where abs(.) represents taking the absolute value, This is the minimum speed increment limit.

5. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by a processor, it implements the steps of the method as described in any one of claims 1 to 4.

6. A processing apparatus, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, it implements the steps of the method as described in any one of claims 1 to 4.

Images