CN103941749A - East-west position keeping and angular momentum unloading combined control method - Google Patents

Abstract

The joint control method of east-west position keeping and angular momentum unloading firstly analyzes the influence of geospherical perturbation and solar light pressure perturbation on satellite orbits for geostationary orbit satellites, and at the same time estimates the impact of angular momentum unloading on satellite orbits; then, Compare the effects of the earth's non-spherical perturbation and angular momentum unloading on the longitude drift rate, and determine the unloading time; and consider the impact on eccentricity, eliminate the possibility of multiple solutions at the unloading time, and determine the optimal unloading time, which is beneficial to the east-west position Keep, and make the eccentricity as small as possible. The method of the invention is suitable for geostationary orbit satellites that need to use thrusters to unload angular momentum intensively. Using the method described in the invention helps to reduce the times of maintaining the east-west positions, save propellant, increase the mission time of satellites, and prolong the service life of satellites.

Description

Joint Control Method of East-West Position Keeping and Angular Momentum Unloading

technical field

The invention relates to a combined control method for east-west position keeping and angular momentum unloading, and belongs to the technical field of geostationary orbit satellite orbit maintenance.

Background technique

Most satellites in geostationary orbit, especially communication satellites, detect that the angular momentum exceeds a certain threshold according to the measurement of the rotational speed of the momentum wheel, and then perform jet unloading. Geostationary orbit navigation satellites and earth observation satellites have strict restrictions on jet. This type of satellite uses momentum wheels to absorb the angular momentum generated by interference such as solar light pressure, and performs centralized unloading when unloading is allowed or when the maximum envelope of angular momentum is reached.

Other than geostationary orbit satellite orbit maintenance design, little consideration is given to the impact of angular momentum unloading on the satellite orbit, and there is no design to use the impulse generated by angular momentum unloading to maintain the east-west position.

For the satellites whose angular momentum is concentrated and unloaded, especially those with asymmetric sail panels or with large accessories, the angular momentum accumulates more, and the impact of jet unloading on the satellite orbit cannot be ignored. If the effects of the earth's non-spherical perturbation, solar light pressure perturbation, and angular momentum unloading are not considered comprehensively, unreasonable selection of unloading time, and reasonable selection of thrusters, in order to overcome the impact of angular momentum unloading on satellite orbits, more propulsion will be consumed Agents are used to maintain satellite orbits and take up more mission time.

Contents of the invention

The technical problem of the present invention is: aiming at the impact of angular momentum unloading on the satellite orbit, rationally select the thruster, design the unloading time, and make full use of the impulse generated by angular momentum unloading to maintain the east-west position. This can effectively save propellant, and at the same time reduce the number of times to maintain the east-west position and increase the mission time.

Technical solution of the present invention is:

The joint control method of east-west position keeping and angular momentum unloading, the steps are as follows:

(1) According to the fixed-point position of the satellite, determine the non-spherical perturbation of the earth that the satellite is subjected to, and further determine the daily variation ΔD _E of the flat longitude drift rate caused by the perturbation of the earth;

(2) According to the satellite structure and quality, determine the solar light pressure perturbation on the satellite, and obtain the components of the eccentricity change rate vector in the X and Y directions of the geocentric inertial coordinate system

(3) According to the satellite structure and quality, determine The relationship between satellite right ascension l and solar argument u _s , where, is the angle between the angular momentum to be unloaded in the orbital plane and the Z-axis of the satellite body;

(4) Calculate the flat longitude drift rate ΔD caused by angular momentum unloading, and the components Δe _x and Δe _y of the eccentricity variation caused by angular momentum unloading in the X and Y directions of the earth-centered inertial coordinate system;

(5) Comparing the daily variation ΔD _E of the flat longitude drift rate caused by the earth perturbation with the magnitude of the flat longitude drift rate ΔD, and then determine value;

(6) The influence of combined solar light pressure and angular momentum unloading on satellite eccentricity, excluding (5) Multiple solutions of , determine the unique value;

(7) According to The value of , combined with the relationship between the angular momentum vector and the sun vector, determines the angle θ between the sun in the orbit plane and the orbit Z axis when unloading, and according to θ, the sub-satellite point local time is obtained;

(8) Perform angular momentum unloading at the sub-satellite point obtained in the step (7), thereby realizing east-west position maintenance.

In the step (3), The relationship between l and u _s is:

When the angular momentum vector lags behind the sun vector by 90°, The angular momentum vector leads the sun vector by 90°,

In the step (4), the calculation of ΔD, Δe _x , Δe _y is specifically:

When the angular momentum vector lags behind the sun vector by 90°,

When the angular momentum vector leads the sun vector by 90°,

Among them, ΔH is the angular momentum that needs to be unloaded every day, T _t is the unloading torque of the X-direction thruster, F _t is the force of the X-direction thruster, T _r is the unloading torque of the Z-direction thruster, F _r is The working force of the Z-direction thruster, m is the mass of the satellite, and V _s is the velocity of the geosynchronous orbit.

In the step (5), Determination of value:

If |ΔD| ≤ ΔD _E |, or

If |ΔD| > |ΔD _E |, Find the inverse cosine the value of is the flat longitude perturbation acceleration caused by the earth perturbation.

In described step (6), get rid of The multiple solutions of are specifically:

Calculate the following formula, considering the joint effect of solar light pressure and angular momentum unloading, and find the corresponding eccentricity minimum solution, sure Unique values for :

When the angular momentum vector lags behind the sun vector by 90°,

When the angular momentum vector leads the sun vector by 90°

Among them, e _x (t), e _y (t) are the components of the eccentricity in the X and Y directions of the geocentric inertial coordinate system at time t, and e _x (t ₀ ), e _y (t ₀ ) are the eccentricity at time t ₀ The components of the rate in the X and Y directions of the geocentric inertial coordinate system, $A=\frac{T_{\mathrm{year}}\&Center\; Dot;\mathrm{\ΔH}}{86400\·2\mathrm{\π}\·V_{\mathrm{the\; s}}\·m},$ T _Year is the number of seconds in a year, k is the solar light pressure constant, A _eff is the effective area of the satellite's solar radiation, and _{is is} the angle between the satellite's yellow and red.

In the step (7), the time when the subsatellite point is the angular momentum unloading moment:

When the angular momentum vector lags behind the sun vector by 90°, Taking the main value of θ between 0° and 360°, then The value of is the sub-satellite point local time;

When the angular momentum vector leads the sun vector by 90°, Taking the main value of θ between 0° and 360°, then The value of is the sub-satellite point local time.

The beneficial effect of the present invention compared with prior art is:

By adopting the method of the invention, the momentum generated by angular momentum unloading can be fully utilized to maintain satellite orbits, realize east-west position maintenance, save propellant, and increase satellite mission time.

Description of drawings

Fig. 1 is a flowchart of the present invention;

Figure 2 shows the annual variation of eccentricity under the joint action of solar light pressure and angular momentum unloading jet, where the abscissa is the component of eccentricity in the X direction of the earth-centered inertial coordinate system, and the ordinate is the eccentricity in the earth-centered inertial coordinate system component in the Y direction.

Detailed ways

The specific implementation manners of the present invention will be further described in detail below in conjunction with the accompanying drawings.

As shown in Fig. 1, the present invention provides a combined control method for east-west position maintenance and angular momentum unloading, the steps are:

(1) According to the fixed-point position of the satellite, determine the non-spherical perturbation of the earth that the satellite is subjected to, and further determine the daily variation ΔD _E of the flat longitude drift rate caused by the perturbation of the earth;

The value of ΔD _E is specifically determined as follows:

${\mathrm{<span\; class="notranslate">\; \&Delta;D</span>}}_{\mathrm{<span\; class="notranslate">\; E.</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 3</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&Center\; Dot;</span>\frac{{\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}}_{\mathrm{<span\; class="notranslate">\; E.</span>}}}{\mathrm{<span\; class="notranslate">\; 540</span>}}$

in, is the flat longitude perturbation acceleration caused by the earth perturbation ( The relationship between the value of and the longitude of a fixed point, see page 248 of "Geostationary Orbit Handbook", National Defense Industry Press, 1999).

(2) According to the satellite structure and quality, determine the solar light pressure perturbation on the satellite, and obtain the components of the eccentricity change rate vector in the X and Y directions of the geocentric inertial coordinate system

The value of is specifically determined as follows:

${\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; e</span>}}}_{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; 3</span>}{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; eff</span>}}\mathrm{<span\; class="notranslate">\; k</span>}}{{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}\mathrm{<span\; class="notranslate">\; m</span>}}\mathrm{<span\; class="notranslate">\; sin</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}\mathrm{<span\; class="notranslate">\; cos</span>}{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}$

${\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; e</span>}}}_{\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; 3</span>}{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; eff</span>}}\mathrm{<span\; class="notranslate">\; k</span>}}{{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}\mathrm{<span\; class="notranslate">\; m</span>}}\mathrm{<span\; class="notranslate">\; cos</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}$

Among them, m is the mass of the satellite, k is the solar light pressure constant, A _ef f is the effective area of solar radiation of the satellite, V _s is the synchronous orbit velocity, u _s is the solar argument, and i _s is the satellite eclipse angle.

(3) According to the satellite structure and quality, determine The relationship between satellite right ascension l and solar argument u _s , where, is the angle between the angular momentum to be unloaded in the orbital plane and the Z-axis of the satellite body;

The relationship between l and u _s is determined in the following way:

When the sun irradiates the surface of the satellite, because of the asymmetry of the solar panels, or there is only one solar panel, the disturbance torque generated acts on the angular momentum vector accumulated on the satellite in the orbital plane, and the included angle with the sun vector is 90°. If the interference force caused by solar light pressure on the +Y plane of the satellite is greater than that on the -Y plane, the angular momentum vector lags behind the sun vector by 90°; otherwise, the angular momentum vector leads the sun vector by 90°. That is, when the angular momentum vector lags behind the sun vector by 90°, The angular momentum vector leads the sun vector by 90°,

(4) Calculate the flat longitude drift rate ΔD caused by angular momentum unloading, and the components Δe _x and Δe _y of the eccentricity variation caused by angular momentum unloading in the X and Y directions of the earth-centered inertial coordinate system;

ΔD, Δe _x , Δe _y are specifically determined as follows:

The increment equation of the orbital elements controlled by the pulse control in the geostationary orbit per unit mass is (see "Satellite Orbital Attitude Dynamics and Control" written by Zhang Ren, Beijing University of Aeronautics and Astronautics Press, 1998)

$\mathrm{<span\; class="notranslate">\; \&Delta;D</span>}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; 3</span>}}{{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}}{\mathrm{<span\; class="notranslate">\; \&Delta;V</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}$

${\mathrm{<span\; class="notranslate">\; \&Delta;\; e</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}}<span\; class="notranslate">\; (</span>{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; \&Delta;V</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}}\mathrm{<span\; class="notranslate">\; siml</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&Delta;V</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}\mathrm{<span\; class="notranslate">\; cos</span>}\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; )</span>$

${\mathrm{<span\; class="notranslate">\; \&Delta;e</span>}}_{\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&Delta;V</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}}\mathrm{<span\; class="notranslate">\; cos</span>}\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&Delta;V</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}\mathrm{<span\; class="notranslate">\; sin</span>}\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; )</span>$

In the formula, ΔV _t and ΔV _r are the orbital tangential and radial pulse velocity increments respectively, and l is the right ascension of the satellite when the orbital pulse acts. ΔD is the flat longitude drift rate caused by angular momentum unloading; Δe _x and _Δey are the components of the eccentricity change caused by angular momentum unloading in the X and Y directions of the geocentric inertial coordinate system, and V _s is the synchronous orbit velocity.

Assuming that the angular momentum that needs to be unloaded every day is ΔH, the unloading moment is T, the force generated during unloading is F, and the positional mass is m, then the velocity increment generated by unloading is

When unloading, the angle between the angular momentum to be unloaded and the Z-axis of the satellite body is , the X-direction thruster works, the unloading moment is T _t , the force generated during unloading is F _t , the Z-direction thruster works, the unloading moment is T _r , the unloading force is F _r , and the mass of the satellite is m, then, The velocity increment produced by one unloading is

Substituting the above formula into the incremental equation and considering the relationship between unloading time and satellite right ascension, then

When the angular momentum vector lags behind the sun vector by 90°,

When the angular momentum vector leads the sun vector by 90°,

(5) Comparing the daily variation ΔD _E of the flat longitude drift rate caused by the earth perturbation with the magnitude of the flat longitude drift rate ΔD, and then determine value;

The value of is specifically determined as follows:

Comparing the effect of the earth's non-spherical perturbation ΔD _E and angular momentum unloading on the orbital longitude ΔD,

If |ΔD| ≤ |ΔD _E |, or

If |ΔD|>|ΔD _E |, let but

Find the inverse cosine value. The value of , there are four solutions.

(6) The influence of combined solar light pressure and angular momentum unloading on satellite eccentricity, excluding (5) Multiple solutions of , determine the unique value;

exclude The multiple solutions of are determined in the following way:

Calculate the following formula to find the minimum corresponding eccentricity solution, sure unique value of

When the angular momentum vector lags behind the sun vector by 90°,

When the angular momentum vector leads the sun vector by 90°

Among them, e _x (t), e _y (t) are the components of the eccentricity in the X and Y directions of the geocentric inertial coordinate system at time t, and e _x (t ₀ ), e _y (t ₀ ) are the eccentricity at time t ₀ The components of the rate in the X and Y directions of the geocentric inertial coordinate system, $A=\frac{T_{\mathrm{year}}\&Center\; Dot;\mathrm{\&Delta;H}}{86400\&Center\; Dot;2\mathrm{\&pi;}\&Center\; Dot;V_{\mathrm{the\; s}}\&Center\; Dot;m},$ T _Year is the number of seconds in a year.

(7) According to The value of , combined with the relationship between the angular momentum vector and the sun vector, determines the angle θ between the sun in the orbit plane and the orbit Z axis when unloading, and according to θ, the sub-satellite point local time is obtained;

The sub-satellite location time of the satellite is specifically determined in the following manner:

When the angular momentum vector lags behind the sun vector by 90°, Taking the main value of θ between 0° and 360°, then The value of is the sub-satellite point local time;

When the angular momentum vector leads the sun vector by 90°, Taking the main value of θ between 0° and 360°, then The value of is the sub-satellite point local time.

(8) Perform angular momentum unloading at the sub-satellite point obtained in the step (7), thereby realizing east-west position maintenance.

Example

For a certain satellite, the mass of the satellite is m=2366.1kg, the effective area of solar radiation A _eff of the satellite is 26.335m ² , and the fixed point is at 112° east longitude. The accumulated angular momentum ΔH=30.4Nms.

(1) According to the fixed-point position of the satellite, determine the non-spherical perturbation of the earth on the satellite, and further determine the daily variation of the flat longitude drift rate ΔD _E = -0.0020°/d caused by the perturbation of the earth;

(2) According to the satellite structure and quality, determine the solar light pressure perturbation on the satellite, and obtain the components of the eccentricity change rate vector in the X and Y directions of the geocentric inertial coordinate system

${\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; e</span>}}}_{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2.2707</span>}<span\; class="notranslate">\; \&times;</span>{\mathrm{<span\; class="notranslate">\; 10</span>}}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 11</span>}}<span\; class="notranslate">\; \&Center\; Dot;</span>\mathrm{<span\; class="notranslate">\; sin</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}$

${\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; e</span>}}}_{\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2.4760</span>}<span\; class="notranslate">\; \&times;</span>{\mathrm{<span\; class="notranslate">\; 10</span>}}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 11</span>}}<span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; cos</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}$

Wherein, the sunlight pressure constant k=4.56×10 ^-6 N/m ² .

(3) According to the satellite structure and quality, determine The relationship between satellite right ascension l and solar argument u _s ;

The interference force caused by the solar light pressure on the +Y plane of the satellite is smaller than that on the -Y plane, and the angular momentum vector is ahead of the sun vector by 90°. but The relationship between satellite right ascension l and solar argument u _s is

According to the relationship between the solar argument u _s and the right ascension l of the satellite, the sub-satellite local time of the unloaded satellite can be determined.

(4) Calculate the flat longitude drift rate ΔD caused by angular momentum unloading, and the components Δe _x and Δe _y of the eccentricity variation caused by angular momentum unloading in the X and Y directions of the earth-centered inertial coordinate system;

Since the fixed-point position is 112° east longitude, and the satellite is perturbed eastward by the earth, when unloading, the X-direction thruster is required to provide -X-direction thrust. So, the moment of unloading is T _t =15Nm, the force generated during unloading is F _t =-9.40N, the Z-direction thruster works, the moment of unloading is T _r =19Nm, and the force generated during unloading is F _r =8.66N, then , the velocity increment produced by one unloading is

Substituting the above formula into the incremental equation and considering the relationship between unloading time and satellite right ascension, then

ΔD＝0.0028°/d

(5) Comparing the daily variation ΔD _E of the flat longitude drift rate caused by the earth perturbation with the magnitude of the flat longitude drift rate ΔD, and then determine value;

Since |ΔD|>|ΔD _E |, then

between (-180°180°], The values of can be -133.86°, -46.14°, 46.14°, 133.86°.

(6) The influence of combined solar light pressure and angular momentum unloading on satellite eccentricity, excluding (5) Multiple solutions of , determine the unique value;

The joint effect of unloading and solar pressure on eccentricity is

At the beginning, let _ex (t ₀ ) and e _y (t ₀ ) be zero, and calculate the eccentricity according to the above formula, as shown in Figure 2. As can be seen from Figure 2, When , under the joint action of solar light pressure and angular momentum unloading jet, the variation range of eccentricity is the smallest. then choose uninstall.

(7) According to The value of , combined with the relationship between the angular momentum vector and the sun vector, determines the angle θ between the sun in the orbit plane and the orbit Z axis when unloading, and according to θ, the sub-satellite point local time is obtained;

, considering that the angular momentum vector is ahead of the sun vector by 90°, the angle between the sun and the orbital Z-axis in the orbital plane is 316.14° when unloading, which corresponds to 21.076 hours local time at the subsatellite point of the satellite. That is, around 21:04:33.6 local time at the subsatellite point of the satellite.

(8) Perform angular momentum unloading at the sub-satellite point obtained in the step (7), thereby realizing east-west position maintenance.

The content that is not described in detail in the specification of the present invention belongs to the well-known technology of those skilled in the art.

Claims (6)

1. The east-west position keeping and angular momentum unloading combined control method is characterized by comprising the following steps of:

(1) according to the fixed point position of the satellite, determining the global non-spherical perturbation suffered by the satellite, and further determining the daily change quantity delta D of the flatness drift rate caused by the global perturbation_E；

(2) According to the structure and the quality of the satellite, determining the solar pressure perturbation on the satellite to obtain the components of the eccentricity ratio change rate vector in the X and Y directions of the geocentric inertial coordinate system

(3) Determining according to satellite structure and qualityRight ascension l and solar amplitude u of satellite_sIn the above-mentioned formula, wherein,the included angle between the angular momentum to be unloaded in the orbital plane and the Z axis of the satellite body is shown;

(4) calculating the flatness drift rate Delta D generated by angular momentum unloading, and the component Delta e of the eccentricity variation caused by the angular momentum unloading in the X and Y directions of the geocentric inertial coordinate system_xAnd Δ e_y；

(5) Comparing daily change quantity delta D of flatness drift rate caused by perturbation of the earth_EAnd the magnitude of the flatness drift rate Delta D, and then determiningA value of (d);

(6) the effect of the combined solar light pressure and angular momentum unloading on the satellite eccentricity is eliminated in (5)Determining a unique value;

(7) according toDetermining the included angle theta between the sun and the Z axis of the orbit in the orbit plane during unloading by combining the relationship between the angular momentum vector and the sun vector, and obtaining the satellite sub-satellite point local time according to the theta;

(8) and (4) carrying out angular momentum unloading at the subsatellite point position of the satellite obtained in the step (7), thereby realizing east-west position maintenance.

2. East-west position-keeping and angular momentum unloading combined control according to claim 1The method is characterized in that: in the step (3), the step (c),l、u_sthe relation of (A) is as follows:

when the angular momentum vector lags the sun vector by 90,the angular momentum vector leads the sun vector by 90 degrees,

3. the east-west position preserving and angular momentum unloading combined control method according to claim 1, wherein: in the step (4), Δ D and Δ e are calculated_x、Δe_yThe method specifically comprises the following steps:

when the angular momentum vector lags the sun vector by 90,

when the angular momentum vector leads the sun vector by 90,

where Δ H is the angular momentum required to be unloaded per day, T_tMoment of unloading for operation of the X-thruster, F_tActing force for the operation of the X-thruster, T_rMoment of unloading for Z-direction thruster operation, F_rIs the acting force of the Z-direction thruster, m is the satellite mass, V_sIs the geosynchronous orbital velocity.

4. The east-west position preserving and angular momentum unloading combined control method according to claim 1, wherein: in the step (5), the step (c),determination of the value:

if | Delta D | ≦ | Delta D_E|，Or

If | Δ D |>|ΔD_E|，Obtaining by inverse cosineThe value of (a) is,acceleration of the planar longitude perturbation caused by the perturbation of the earth.

5. The east-west position preserving and angular momentum unloading combined control method according to claim 1, wherein: in the step (6), excludingThe multiple solutions are as follows:

calculating the following formula, and searching the one with the minimum corresponding eccentricity ratio under the combined action of sunlight pressure and angular momentum unloadingTo determineUnique value of (c):

when the angular momentum vector lags the sun vector by 90,

when the angular momentum vector leads the sun vector by 90 DEG

Wherein e is_x(t)、e_y(t) is the component of eccentricity in X and Y directions of the geocentric inertial coordinate system at time t, e_x(t₀)、e_y(t₀) Is t₀The components of the moment eccentricity in the X and Y directions of the geocentric inertial frame, <math> <mrow> <mi>A</mi> <mo>=</mo> <mfrac> <mrow> <msub> <mi>T</mi> <mi>Year</mi> </msub> <mo>&CenterDot;</mo> <mi>&Delta;H</mi> </mrow> <mrow> <mn>86400</mn> <mo>&CenterDot;</mo> <mn>2</mn> <mi>&pi;</mi> <mo>&CenterDot;</mo> <msub> <mi>V</mi> <mi>s</mi> </msub> <mo>&CenterDot;</mo> <mi>m</mi> </mrow> </mfrac> <mo>,</mo> </mrow> </math> T_Yearis the number of seconds of a year, k is the solar light pressure constant, A_effIs the effective area of solar radiation of the satellite, i_sThe included angle between yellow and red of the satellite.

6. The east-west position preserving and angular momentum unloading combined control method according to claim 1, wherein: the satellite lower point in the step (7) is the angular momentum unloading moment:

when the angular momentum vector lags the sun vector by 90,for theta between 0 DEG and 360 DEG, thenThe value of (1) is the satellite sub-satellite point local time;

when the angular momentum vector leads the sun vector by 90,for theta between 0 DEG and 360 DEG, thenThe value of (c) is the satellite sub-satellite point-in-place time.