An angular momentum offloading method and device based on trajectory prediction

Abstract

The application discloses a kind of angular momentum unloading method and device based on trajectory prediction, belong to the field of spacecraft attitude control.Method includes: according to the inertia characteristic of satellite, geometry configuration and solar angle, the combined interference torque of satellite is calculated by the solar pressure resultant moment and gravity gradient moment in body coordinate system;According to the numerical integral average of X axis and Z axis interference torque of satellite under different solar angles and the numerical integral function maximum of Y axis interference torque predicted by combined interference torque, and according to the predicted result, the angular momentum pre-bias value of three axes of satellite under different solar angles is obtained;According to angular momentum pre-bias value, the target angular momentum corresponding to the solar angle of satellite at preset unloading time is found, and according to the search result, the jet control moment of thruster of satellite is adjusted, to unload satellite angular momentum to target angular momentum.The application can effectively reduce the fluctuation range of the accumulated angular momentum of star body and the rotation speed of momentum wheel.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of spacecraft attitude control, and particularly relates to a method and device for angular momentum unloading based on trajectory prediction. BACKGROUND

[0002] The angular momentum variation law of asymmetrically configured spacecraft in orbit is complex under the action of gravity gradient and solar radiation pressure. Generally, the daily cumulative angular momentum of a satellite with good symmetry is about 10 Nms, and the daily cumulative angular momentum of a complexly configured spacecraft, especially a satellite with a large antenna, is about 90-100 Nms. Therefore, the angular momentum pre-bias method is needed to reduce the fluctuation range of the satellite wheel system angular momentum of a large asymmetrically configured spacecraft in orbit, and the data set is used to set the target angular momentum of the satellite to adapt to the phase angle and elevation angle changes of the solar vector.

[0003] In the related art, a fixed unloading value is usually used to unload the angular momentum of a simple configuration spacecraft, but this method cannot describe the variation law of the pre-biased angular momentum when facing an asymmetrically configured spacecraft, which greatly increases the fluctuation range of the satellite cumulative angular momentum and the fluctuation range of the momentum wheel rotation speed.

[0004] Therefore, there is an urgent need for a method and device for angular momentum unloading based on trajectory prediction to solve the above technical problems. SUMMARY

[0005] The present application provides a method and device for angular momentum unloading based on trajectory prediction, which can effectively reduce the fluctuation range of the satellite cumulative angular momentum and the momentum wheel rotation speed. The technical solution is as follows:

[0006] On the one hand, a method for angular momentum unloading based on trajectory prediction is provided, which comprises:

[0007] According to the inertia characteristics, geometric configuration and solar angle of the satellite, the combined interference torque of the satellite in the body coordinate system composed of the solar radiation pressure resultant moment and the gravity gradient moment is calculated;

[0008] The numerical integral average values of the X-axis and Z-axis interference torques and the numerical integral function maximum value of the Y-axis interference torque of the satellite under different solar angles are predicted according to the combined interference torque, and the angular momentum pre-bias values of the three axes of the satellite under different solar angles are calculated according to the prediction results;

[0009] The target angular momentum corresponding to the solar angle of the satellite at the preset unloading time is found according to the angular momentum pre-bias values, and the control torque of the satellite is adjusted according to the finding result to unload the satellite angular momentum to the target angular momentum.

[0010] On the other hand, an angular momentum unloading device based on trajectory prediction is provided, the device comprising:

[0011] The calculation module is used to calculate the combined disturbance torque of the satellite in its body coordinate system, which is composed of the solar radiation pressure torque and the gravity gradient torque, based on the satellite's inertial characteristics, geometric configuration and solar angle.

[0012] The prediction module is used to predict the numerical integral average value of the interference torque of the satellite on the X-axis and Z-axis under different solar angles, as well as the maximum and minimum values ​​of the numerical integral function of the interference torque on the Y-axis, and to calculate the angular momentum pre-bias values ​​of the satellite on the three axes under different solar angles based on the prediction results.

[0013] The unloading module is used to find the target angular momentum of the satellite at the solar angle at the preset unloading time based on the angular momentum pre-bias value, and adjust the thruster jet control torque of the satellite according to the search result to unload the satellite's angular momentum to the target angular momentum.

[0014] On the other hand, a computer device is provided, the computer device including a memory and a processor, the memory for storing a computer program, and the processor for executing the computer program stored in the memory to implement the steps of the trajectory prediction-based angular momentum unloading method described above.

[0015] On the other hand, a computer-readable storage medium is provided, wherein a computer program is stored therein, and when the computer program is executed by a processor, it implements the steps of the above-described trajectory prediction-based angular momentum unloading method.

[0016] On the other hand, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps of the trajectory prediction-based angular momentum unloading method described above.

[0017] The technical solution provided by this invention can bring at least the following beneficial effects: First, the combined disturbance torque, consisting of the solar radiation pressure torque and the gravitational gradient torque, is calculated on the satellite based on the solar vector azimuth. Then, a pre-bias method is used to predict and set the satellite's angular momentum pre-bias value at different solar angles based on the combined disturbance torque. Finally, the target angular momentum is obtained by finding the corresponding angular momentum pre-bias value based on the solar angle at the target unloading time, and the satellite's angular momentum is unloaded to the target angular momentum using a centralized, timed unloading method with thrusters. This method uses centralized, timed unloading logic to unload the satellite's angular momentum to the target angular momentum, effectively reducing the fluctuation range of the momentum wheel speed. This method can be extended to the angular momentum management design of other large asymmetric configuration spacecraft in geosynchronous orbit. Attached Figure Description

[0018] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the following will briefly introduce the drawings needed to be used in the embodiments or prior art description. Obviously, the drawings described below are some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor on the basis of these drawings.

[0019] Figure 1 is a flow chart of a method for angular momentum offloading based on trajectory prediction provided by an embodiment of the present application;

[0020] Figure 2 is a schematic diagram of a complex satellite structure provided by an embodiment of the present application;

[0021] Figure 3 is a calculation result of an X-axis angular momentum pre-bias value provided by an embodiment of the present application;

[0022] Figure 4 is a calculation result of a Y-axis angular momentum pre-bias value provided by an embodiment of the present application;

[0023] Figure 5 is a calculation result of a Z-axis angular momentum pre-bias value provided by an embodiment of the present application;

[0024] Figure 6 is a simulation result schematic diagram of offloading according to the angular momentum pre-bias value provided by an embodiment of the present application;

[0025] Figure 7 is a simulation result schematic diagram of not offloading by using the angular momentum pre-bias value provided by an embodiment of the present application;

[0026] Figure 8 is a structural diagram of an angular momentum offloading device based on trajectory prediction provided by an embodiment of the present application;

[0027] Figure 9 is a hardware architecture diagram of a computer device provided by an embodiment of the present application. DETAILED DESCRIPTION

[0028] In order to make the objects, technical solutions and advantages of the embodiments of the present application clearer, the following will combine the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are some embodiments of the present application, but not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the present application.

[0029] As described above, the traditional angular momentum unloading method often leads to the increase of the fluctuation range of the accumulated angular momentum of the satellite and the fluctuation range of the momentum wheel rotation speed in the unloading process when facing large asymmetric configuration spacecrafts due to the inconsistent angular momentum requirements of the spacecrafts.

[0030] Based on this, the concept of the present application is to determine the angular momentum pre-offset value by using the orbit prediction method, and to unload the angular momentum of the satellite to the target angular momentum by using the timing unloading logic, thereby effectively reducing the fluctuation range of the momentum wheel rotation speed.

[0031] The specific implementation of the above concept is described below.

[0032] Please refer to Figure 1 The angular momentum unloading method based on orbit prediction provided by the embodiment of the present application comprises the following steps:

[0033] In step 100, the combined interference torque of the satellite under the body coordinate system composed of the combined torque of the solar pressure and the gravity gradient torque is calculated according to the inertia characteristics, geometric configuration and solar angle of the satellite.

[0034] In step 102, the numerical integral average values of the X-axis and Z-axis interference torques and the maximum values of the numerical integral function of the Y-axis interference torque of the satellite under different solar angles are predicted according to the combined interference torque, and the angular momentum pre-offset values of the three axes of the satellite under different solar angles are calculated according to the prediction results.

[0035] In step 104, the target angular momentum corresponding to the solar angle of the satellite at the preset unloading time is found according to the angular momentum pre-offset value, and the jet control torque of the thruster of the satellite is adjusted according to the finding result, so as to unload the angular momentum of the satellite to the target angular momentum.

[0036] In the embodiment of the present application, first, the combined interference torque of the satellite composed of the combined torque of the solar pressure and the gravity gradient torque is calculated according to the solar vector direction, then the angular momentum pre-offset values of the satellite under different solar angles are predicted and set according to the combined interference torque by using the pre-offset method, finally, the unloading target angular momentum is obtained by finding the corresponding angular momentum pre-offset value according to the solar angle at the target unloading time, and the angular momentum of the satellite is unloaded to the target angular momentum by using the concentrated timing unloading mode of the thruster. This method uses the concentrated timing unloading logic to unload the angular momentum of the satellite to the target angular momentum, effectively reduces the fluctuation range of the momentum wheel rotation speed, and the related method can be applied to the angular momentum management design of other large asymmetric configuration spacecrafts in geosynchronous orbit.

[0037] The execution modes of each step shown in the following description. Figure 1

[0038] ​First, for step 100, according to the inertia characteristics, geometric configuration and solar angle of the satellite, a combined interference torque of the satellite in the body coordinate system is calculated, which is composed of a solar radiation pressure resultant torque and a gravity gradient torque.

[0039] In the embodiment of the present application, the combined interference torque is calculated by the following steps: according to the phase angle and elevation angle of the sun in the satellite orbit system, a unit direction vector of the sun in the satellite orbit system is calculated; a single solar radiation pressure torque suffered by each configuration panel of the satellite is calculated according to the unit direction vector, and a solar radiation pressure resultant torque suffered by the satellite is calculated according to the single solar radiation pressure torque; a gravity gradient torque suffered by the satellite is calculated according to the orbital angular velocity and rotational inertia of the satellite; and the solar radiation pressure resultant torque and the gravity gradient torque are fused to obtain the combined interference torque of the satellite in the body coordinate system.

[0040] Specifically, as shown in Figure 2 , Figure 2 A complex non-symmetrical configuration satellite system structure schematic diagram is provided for the embodiment, and the satellite is composed of a satellite central body, solar wings and a reflector antenna. The surfaces involved in the light pressure calculation in the present application include east panels, west panels, south panels, north panels, earth panels, back earth panels, south solar wings, north solar wings, upper surfaces of the reflector antenna, lower surfaces of the reflector antenna, side surfaces 1-18 of the reflector antenna. The satellite body coordinate system O B X B Y B Z B of the present application has its coordinate origin O B located at the satellite center of mass; the O B Z B axis passes through the coordinate origin O B and is perpendicular to the satellite-rocket separation surface, along the longitudinal axis direction of the satellite, and the positive direction points to the earth-facing surface direction; the O B X B axis passes through the coordinate origin O B and is located in the satellite-rocket separation surface, parallel to the east panel theoretical normal direction, the positive direction is consistent with the east panel outer normal direction, and points to the normal flight direction of the satellite; the O B Y B axis passes through the coordinate origin O B and is located in the satellite-rocket separation surface, forming a right-hand system with the O B X B axis and the O O Z O axis. The satellite orbit system O O X O Y O Z O of the present application has its coordinate origin O O located at the satellite center of mass, and the O OY O axis points to the negative normal direction of the orbital plane, O O X O axis points to the positive normal direction of the orbital plane, O O Y O , O O Z O constitutes a right-hand system (points to the forward direction of the satellite).

[0041] First, the unit position vector i of the sun in the satellite orbital system is calculated by the following formula Sun :

[0042]

[0043] In the formula, i Sun is the unit position vector of the sun in the satellite orbital coordinate system, which points to the sun; α Sun is the sun phase angle, which is the projection of the unit position vector i Sun in the satellite orbital coordinate system X O O O Z O plane and the angle between the O O Z O axis, unit rad, positive to the +O O X O direction; β Sun is the sun elevation angle, which is the angle between the unit position vector i Sun and the satellite orbital coordinate system X O O O Z O plane, unit rad, positive to the +O O Y O direction.

[0044] Then, the single solar pressure moment of each configuration panel of the satellite is calculated according to the unit position vector.

[0045] The solar pressure moment T E (α Sun ,β Sun ) of the satellite on the east panel is calculated by the following formula:

[0046] T E (α Sun ,β Sun ) = (r EB ) × F EB (α Sun ,β Sun ) S E

[0047] wherein () × is the vector cross product matrix, and T E(α Sun ,β Sun ) is the solar radiation pressure torque on the east panel in the satellite body coordinate system, unit Nm, which is a function of α Sun and β Sun ; r EB is the position vector from the satellite center of mass to the geometric center of the east panel in the satellite body coordinate system, unit m; S E is the area of the east panel, unit m 2 ; F EB (α Sun ,β Sun ) is the solar radiation pressure on the east panel in the satellite body coordinate system, unit N, and the specific calculation formula is as follows:

[0048] F EB (α Sun ,β Sun )=EH(cosθ E )(1-η E )cosθ E ((1-C RSE )s B

[0049] -2(C RSE cosθ E +C RDE / 3)n BE ) / C

[0050] In the formula, E is the solar radiation intensity; C is the speed of light; H(x) is the Heaviside function, which is 1 when x≥0 and 0 when x<0; η E is the light transmittance of the east panel; s B is the unit vector of the satellite body coordinate system, equal to -i Sun ; n BE is the outer normal vector of the east panel in the satellite body coordinate system; θ E is the angle between the unit vector s B of the sun's incidence and the inner normal vector -n BE of the east panel; C RSE is the mirror reflection coefficient of the east panel, and C RDE is the diffuse reflection coefficient of the east panel.

[0051] The solar radiation pressure torque T W (α Sun ,β Sun ) on the west panel is calculated by the following formula:

[0052] T W (α Sun ,β Sun )=(r WB )× F WB (α Sun ,β Sun )S W

[0053] wherein r WB is the position vector from the satellite center of mass to the geometric center of the west panel in the satellite body coordinate system, with unit m; S W is the area of the west panel, with unit m 2 ; F WB (α Sun ,β Sun ) is the solar radiation pressure on the west panel of the satellite in the satellite body coordinate system, with unit N, and the specific calculation formula is as follows:

[0054] F WB (α Sun ,β Sun ) = EH (cosθ W ) (1-η W ) cosθ W ((1-C RSW )s B -2(C RSW cosθ W +C RDW / 3)n BW ) / C

[0055] wherein η W is the light transmittance of the west panel; n BW is the outer normal vector of the west panel in the satellite body coordinate system; θ W is the angle between the unit vector s B of the solar incidence and the inner normal vector -n BW of the west panel, C RSW is the mirror reflection coefficient of the west panel, and C RDW is the diffuse reflection coefficient of the west panel.

[0056] The solar radiation pressure moment T S (α Sun ,β Sun ) on the south panel of the satellite is calculated by the following formula:

[0057] T S (α Sun ,β Sun ) = (r SB ) × F SB (α Sun ,β Sun ) S S

[0058] wherein r SBis the position vector from the satellite center of mass to the geometric center of the south panel in the satellite body coordinate system, with the unit of m; S S is the area of the south panel, with the unit of m 2 ; F SB (α Sun ,β Sun ) is the solar radiation pressure on the south panel of the satellite in the satellite body coordinate system, with the unit of N, and the specific calculation formula is as follows:

[0059] F SB (α Sun ,β Sun ) = EH (cosθ S ) (1-η S ) cosθ S ((1-C RSS ) s B

[0060] -2(C RSS cosθ S +C RDS / 3) n BS ) / C

[0061] wherein, η S is the light transmittance of the south panel; n BS is the outer normal vector of the south panel in the satellite body coordinate system; θ S is the angle between the unit vector s B of the solar incidence and the inner normal vector -n BS of the south panel, C RSS is the mirror reflection coefficient of the south panel, and C RDS is the diffuse reflection coefficient of the south panel.

[0062] The solar radiation pressure moment T N (α Sun ,β Sun ) on the north panel of the satellite is calculated by the following formula:

[0063] T N (α Sun ,β Sun ) = (r NB ) × F NB (α Sun ,β Sun ) S N

[0064] wherein, r NB is the position vector from the satellite center of mass to the geometric center of the north panel in the satellite body coordinate system, with the unit of m; S N is the area of the north panel, with the unit of m 2 ; F NB (αSun ,β Sun ) is the solar radiation pressure on the north panel in the satellite body coordinate system, with the unit of N, and the specific calculation formula is as follows:

[0065] F NB (α Sun ,β Sun )=EH(cosθ N )(1-η N )cosθ N ((1-C RSN )s B

[0066] -2(C RSN cosθ N +C RDN / 3)n BN ) / C

[0067] In the formula, η N is the light transmittance of the north panel; n BN is the outer normal vector of the north panel in the satellite body coordinate system; θ N is the angle between the unit vector s B of the solar incidence and the inner normal vector -n BN of the north panel, C RSN is the mirror reflection coefficient of the north panel, and C RDN is the diffuse reflection coefficient of the north panel.

[0068] The solar radiation pressure moment T D (α Sun ,β Sun ) on the earth panel is calculated by the following formula:

[0069] T D (α Sun ,β Sun )=(r DB ) × F DB (α Sun ,β Sun )S D

[0070] Wherein, r DB is the position vector from the satellite center of mass to the geometric center of the earth panel in the satellite body coordinate system, with the unit of m; S D is the area of the earth panel, with the unit of m 2 ; and F DB (α Sun ,β Sun ) is the solar radiation pressure on the earth panel in the satellite body coordinate system, with the unit of N, and the specific calculation formula is as follows:

[0071] F DB (α Sun ,β Sun )=EH(cosθ D )(1-η D cosθ D ((1-C RSD )s B

[0072] -2(C RSD cosθ D +C RDD / 3)n BD ) / C

[0073] In the formula, η D The light transmittance of the floor panel; n BD θ is the outward normal vector of the Earth panel in the satellite body coordinate system; N s is the unit vector of solar incidence. B With respect to the inner normal vector of the ground panel -n BD The included angle, C RSD C is the specular reflection coefficient of the floor panel. RDD denoted as the diffuse reflection coefficient of the floor panel.

[0074] The solar radiation pressure torque T exerted on the satellite by the back panel is calculated using the following formula. U (α Sun ,β Sun ):

[0075] T U (α Sun ,β Sun )=(r UB ) × F UB (α Sun ,β Sun )S U

[0076] Where, r UB S is the position vector from the satellite's center of mass to the geometric center of the back panel in the satellite's body coordinate system, in meters (m). U The area of ​​the floor panel is expressed in meters (m²). 2 ;F UB (α Sun ,β Sun The solar radiation pressure exerted on the satellite by the back panel in the satellite's body coordinate system is expressed in nanometers (N). The specific calculation formula is as follows:

[0077] F UB (α Sun ,β Sun )=EH(cosθ U )(1-ηU )cosθ U ((1-C RSU )s B

[0078] -2(C RSU cosθ U +C RDU / 3)n BU ) / C

[0079] η U =transmittance of the backside panel;n BU =the outer normal vector of the backside panel in the satellite body coordinate system;θ U =the angle between the unit vector s B of the sun's incidence and the inner normal vector -n BU of the backside panel;C RSU =the specular reflection coefficient of the backside panel;C RDU =the diffuse reflection coefficient of the backside panel.

[0080] The solar pressure moment T SA (α Sun ,β Sun ) on the satellite is calculated by the following formula:

[0081] T SA (α Sun ,β Sun ) = (r SAB ) × F SAB (α Sun ,β Sun ) S SA

[0082] wherein r SAB is the position vector from the satellite center of mass to the geometric center of the south solar wing in the satellite body coordinate system, with a unit of m; S SA is the area of the south solar wing, with a unit of m 2 ; F SAB (α Sun ,β Sun ) is the solar pressure on the satellite at the south solar wing in the satellite body coordinate system, with a unit of N, and the specific calculation formula is as follows:

[0083] F SAB (α Sun ,β Sun ) = EH(cosθ SA )(1-η SA )cosθ SA ((1-C RSSA )s B

[0084] -2(C RSSA cosθ SA +C RDSA / 3)n BSA ) / C

[0085] wherein η SA is the transmittance of the south solar wing; n BSA is the outer normal vector of the south solar wing in the satellite body coordinate system; θ SA is the angle between the unit vector s B of the sun incidence and the inner normal vector -n BSA of the south solar wing; C RSSA is the mirror reflection coefficient of the south solar wing; and C RDSA is the diffuse reflection coefficient of the south solar wing.

[0086] The solar pressure moment T NA (α Sun ,β Sun ) of the satellite on the north solar wing is calculated by the following formula:

[0087] T NA (α Sun ,β Sun ) = (r NAB ) × F NAB (α Sun ,β Sun ) S NA

[0088] wherein r NAB is the position vector from the satellite center of mass to the geometric center of the north solar wing in the satellite body coordinate system, with a unit of m; S NA is the area of the south solar wing, with a unit of m 2 ; and F NAB (α Sun ,β Sun ) is the solar pressure of the satellite on the south solar wing in the satellite body coordinate system, with a unit of N, and the specific calculation formula is as follows:

[0089] F NAB (α Sun ,β Sun ) = EH(cosθ NA )(1-η NA )cosθ NA ((1-C RSNA )s B

[0090] -2(C RSNA cosθ NA +C RDNA / 3)n BNA ) / C

[0091] In the formula, η NA n represents the light transmittance of the north solar array. BNA θ is the outward normal vector of the north solar array in the satellite body coordinate system. NA s is the unit vector of solar incidence. B With the inner normal vector of the north solar wing -n BNA The included angle, C RSNA C is the specular reflectance coefficient of the north solar array. RDNA This represents the diffuse reflectance coefficient of the north solar array.

[0092] The solar pressure torque T exerted on the upper surface of the reflector antenna by the satellite is calculated using the following formula. UAT (α Sun ,β Sun ):

[0093] T UAT (α Sun ,β Sun )=(r UATB ) × F UATB (α Sun ,β Sun )S UAT

[0094] Where, r UATB S is the position vector in the satellite's body coordinate system from the satellite's center of mass to the geometric center of the upper surface of the reflector antenna, in meters (m). UAT The area of ​​the upper surface of the reflector antenna, in meters. 2 ;F UATB (α Sun ,β Sun The solar radiation pressure exerted on the upper surface of the reflector antenna by the satellite in the satellite's body coordinate system is expressed in nanometers (N). The specific calculation formula is as follows:

[0095] F UATB (α Sun ,β Sun )=EH(cosθ UAT )(1-η UAT cosθ UAT ((1-C RSUAT )s B

[0096] -2(C RSUAT cosθ UAT +C RDUAT / 3)n BUAT ) / C

[0097] In the formula, η UAT n is the transmittance of the upper surface of the reflector antenna;BUAT is the outer normal vector of the upper surface of the reflector antenna in the satellite body coordinate system; θ UAT is the unit vector of the sun incidence s B is the inner normal vector of the upper surface of the reflector antenna -n BUAT is the included angle of the inner normal vector of the upper surface of the reflector antenna and the outer normal vector of the upper surface of the reflector antenna, C RSUAT is the mirror reflection coefficient of the upper surface of the reflector antenna, C RDUAT is the diffuse reflection coefficient of the upper surface of the reflector antenna.

[0098] The solar pressure moment T experienced by the satellite on the lower surface of the reflector antenna is calculated by the following formula DAT (α Sun ,β Sun ):

[0099] T DAT (α Sun ,β Sun )=(r DATB ) × F DATB (α Sun ,β Sun )S DAT

[0100] wherein r DATB is the position vector from the satellite center of mass to the geometric center of the lower surface of the reflector antenna in the satellite body coordinate system, with a unit of m; S DAT is the area of the lower surface of the reflector antenna, with a unit of m 2 ; F DATB (α Sun ,β Sun ) is the solar pressure experienced by the satellite on the lower surface of the reflector antenna in the satellite body coordinate system, with a unit of N, and the specific calculation formula is as follows:

[0101] F DATB (α Sun ,β Sun )=EH(cosθ DAT )(1-η DAT )cosθ DAT ((1-C RSDAT )s B

[0102] -2(C RSDAT cosθ DAT +C RDDAT / 3)n BDAT ) / C

[0103] In the formula, η DAT is the light transmittance of the lower surface of the reflector antenna; n BDAT is the outer normal vector of the lower surface of the reflector antenna in the satellite body coordinate system; θDAT is the unit vector of the sun's incidence s B is the angle between the inner normal vector -n of the lower surface of the reflector antenna BDAT C RSDAT is the specular reflection coefficient of the lower surface of the reflector antenna RDDAT is the diffuse reflection coefficient of the lower surface of the reflector antenna

[0104] The solar pressure moment T of the satellite on the reflector antenna side surface 1, the reflector antenna side surface 2, …, the reflector antenna side surface 18 is calculated by the following formula SATi (α Sun ,β Sun )(i = 1, 2, …, 18):

[0105] T SATi (α Sun ,β Sun ) = (r SATiB ) × F SATiB (α Sun ,β Sun ) S SATi

[0106] wherein r SATiB is the position vector from the satellite mass center to the geometric center of the reflector antenna side surface i in the satellite body coordinate system, with a unit of m; S SATi is the area of the reflector antenna side surface i, with a unit of m 2 ; F SATiB (α Sun ,β Sun ) is the solar pressure of the satellite on the reflector antenna side surface i in the satellite body coordinate system, with a unit of N, and the specific calculation formula is as follows:

[0107] F SATiB (α Sun ,β Sun ) = EH (cosθ SATi )(1-η SATi ) cosθ SATi ((1-C RSSATi )s B

[0108] -2(C RSSATi cosθ SATi +C RDSATi / 3)n BSATi ) / C

[0109] In the formula, η SATi is the light transmittance of the reflector antenna side surface i; n BSATi is the outer normal vector of the reflector antenna side surface i in the satellite body coordinate system; θSATi is the unit vector of the sun's incidence s B is the angle between the inner normal vector -n of the reflector antenna side surface i BSATi C is the angle between the inner normal vector -n of the reflector antenna side surface i RSSATi C is the mirror reflection coefficient of the reflector antenna side surface i RDSATi C is the diffuse reflection coefficient of the reflector antenna side surface i

[0110] The resultant torque T of the solar radiation pressure on the satellite is calculated according to the above formula SP (α Sun ,β Sun ) :

[0111]

[0112] Further, the resultant torque T of the gravity gradient force on the satellite is calculated according to the orbital angular velocity and the moment of inertia of the satellite G :

[0113]

[0114] In the formula, ω Orb is the orbital angular velocity, with the unit of rad / s; k B is the unit vector pointing to the center of the earth in the satellite body coordinate system; I B is the moment of inertia of the satellite, with the unit of kgm 2 .

[0115] Finally, the resultant interference torque T of the satellite in the body coordinate system is calculated according to the following formula B (α Sun ,β Sun ) :

[0116]

[0117] In the formula, T B (α Sun ,β Sun ) is the resultant interference torque of the satellite in the body coordinate system under the action of the solar radiation pressure and the gravity gradient, with the unit of Nm; T Bx (α Sun ,β Sun ), T By (α Sun ,β Sun ), and T Bz (α Sun ,β Sun ) are the components of the resultant interference torque of the satellite in the body coordinate system under the action of the solar radiation pressure and the gravity gradient along the X, Y, and Z axes of the satellite body, with the unit of Nm.

[0118] Then, for step 102, the numerical integral average values of the X-axis and Z-axis interference torques of the satellite under different sun angles are predicted according to the synthetic interference torques, and the extreme values of the numerical integral function of the Y-axis interference torque are predicted, and the angular momentum pre-offset values of the three axes of the satellite under different sun angles are calculated according to the prediction results.

[0119] Due to the weak intensity of the earth magnetic field near the geosynchronous orbit, the satellite in the orbit usually uses a thruster for centralized timed unloading to eliminate the accumulated angular momentum of the satellite. Therefore, the trajectory prediction is introduced into the process of angular momentum unloading in the embodiment of the application, that is, the angular momentum pre-offset value required by the satellite to be unloaded at the corresponding moment is predicted according to the interference torques suffered by the satellite at different moments when the satellite is in different running trajectories.

[0120] Specifically, the sun phase angle α Sun is divided into discrete values corresponding to different orbit positions by the following formula:

[0121]

[0122] Wherein, j is the serial number of the sun phase angle, and m is the number of divisions of the sun phase angle.

[0123] The sun elevation angle β Sun is divided into discrete values corresponding to different orbit positions by the following formula:

[0124]

[0125] Wherein, β Sunk is the threshold value of the sun angle β Sun , greater than 0, unit rad, k is the serial number of the sun elevation angle, and n is the number of divisions of the sun elevation angle.

[0126] Further, according to the component models T Bx (α Sun ,β Sun ), T By (α Sun ,β Sun ), T Bz (α Sun ,β Sun ) of the synthetic interference torques along the X, Y and Z axes of the satellite, the numerical integral average values of the X-axis and Z-axis interference torques suffered by the satellite in the body coordinate system within one orbit period when the sun angle β Sun takes β Sun1 , β Sun2 , …, β Sunn are calculated in sequence.

[0127]

[0128] wherein, are the different discrete values of the sun elevation angle Sunk are the numerical integral average values of the X and Z axis disturbance torques that the satellite receives in the body coordinate system in one orbit period, units Nm, when the sun phase angle

[0129] Further, according to the component model, the sun phase angle and the sun elevation angle take different discrete values, i.e. α Sun and β Sun take values α Sun1 and β Sun1 , α Sun1 and β Sun2 , …, α Sun1 and β Sunn , α Sun2 and β Sun1 , α Sun2 and β Sun2 , …, α Sun2 and β Sunn , …, α Sunm and β Sun1 , α Sunm and β Sun2 , …, α Sunm and β Sunn , the maximum and minimum values of the numerical integral function of the Y axis disturbance torque that the satellite receives in the body coordinate system in the next orbit period are:

[0130]

[0131] wherein, h ymaxjk , h yminjk are the maximum and minimum values of the numerical integral function of the Y axis disturbance torque that the satellite receives in the body coordinate system in the next orbit period when the sun angle takes each sun angle α Sunj and β Sunk , units Nms; max() is the maximum value function, and min() is the minimum value function.

[0132] Finally, the three-axis angular momentum pre-bias values of the satellite under different sun angles are calculated according to the following formula, as shown in FIG. 4: Figures 3-5

[0133]

[0134] In the formula, H xtk is the pre-bias value of the angular momentum along the X axis of the satellite in the body coordinate system when the sun angle takes value β Sunk ; H ytjk is the pre-bias value of the angular momentum along the Y axis of the satellite in the body coordinate system when the sun angle takes values α Sunj and β Sunk ; H​ztk The solar angle is set to β. Sunk The pre-offset value of the angular momentum along the Z-axis of the satellite in the satellite body coordinate system.

[0135] For step 104, the target angular momentum of the satellite at the solar angle at the preset unloading time is found based on the angular momentum pre-bias value, and the thruster jet control torque of the satellite is adjusted according to the search result to unload the satellite's angular momentum to the target angular momentum.

[0136] In this embodiment of the invention, the target angular momentum is obtained as follows: the solar phase angle and solar elevation angle at the unloading time are calculated based on the unit azimuth vector of the sun at the unloading time; the solar phase angle and solar elevation angle at the unloading time are rounded to obtain the discrete index of the solar angle corresponding to the unloading time; the angular momentum pre-bias value of the three axes is found based on the discrete index of the solar angle at the unloading time to obtain the corresponding target angular momentum.

[0137] Specifically, when the satellite is at the set unloading time, based on the unit vector i of the sun on the satellite in the satellite orbital coordinate system... SunSat =[i SunSatx i SunSaty i SunSatz ] T Calculate the solar angle α SunSat and β SunSat :

[0138] α SunSat =mod(arctan2(i SunSatx i SunSatz ),2π)

[0139] β SunSat =arcsin(i SunSaty )

[0140] Among them, i SunSatx i SunSaty and i SunSatz The following are the unit vectors i of the Sun on the star in the satellite's orbital coordinate system. SunSat Components along the X, Y, and Z axes in the satellite orbit coordinate system; α SunSat The solar azimuth vector i for the unloading time calculated on the satellite SunSat In the satellite orbit coordinate system X O O O Z O Projection in the plane and O O Z O Angle between axes, in rad, deflection +O O X O The direction is positive; β SunSatThe sun position vector i at the unloading time for the on-board computation Sun The satellite orbit coordinate system X O O O Z O The angle between the plane and the direction of +O, unit rad, the bias +Y O Y O The direction is positive; mod() is the integer function.

[0141] Further, according to the sun phase angle α SunSat and the sun elevation angle β SunSat of the satellite at the set unloading time, the discrete index values j SunSat and k SunSat of the pre-bias values H txk , H tyjk , H tzk of the angular momentum along the X, Y and Z axes of the satellite body in the satellite body coordinate system in the pre-bias angular momentum calculation module 2 are calculated:

[0142]

[0143] Where int() is the integer function.

[0144] Finally, according to the index values j SunSat and k SunSat , the corresponding unloading target angular momentum is selected from the pre-bias values H txk , H tyjk , H tzk of the angular momentum along the X, Y and Z axes of the satellite body in the satellite body coordinate system in the pre-bias angular momentum calculation module 2:

[0145]

[0146] Where H tx , H ty and H tz are the unloading target angular momentum along the X, Y and Z axes of the satellite body in the satellite body coordinate system.

[0147] After obtaining the unloading target angular momentum, the thruster is used to exert a control torque on the satellite body to unload the on-board computed angular momentum to the target angular momentum H tx , H ty and H tz , thereby completing the unloading of the angular momentum of the large asymmetrically configured satellite.

[0148] As shown in Figure 6 and Figure 7 , Figure 6Simulation results of the pre-biased angular momentum provided by the embodiment for angular momentum unloading. The calculation results show that in the satellite body coordinate system, the fluctuation range of the X-axis angular momentum is-19.93 Nms~ -5.7 Nms, the fluctuation range of the Y-axis angular momentum is-43.26 Nms~ 44.59 Nms, and the fluctuation range of the Z-axis angular momentum is-16.69 Nms~ -1.431 Nms; Figure 7 Simulation results of the pre-biased angular momentum (X: 0 Nms, Y: 0 Nms, Z: 0 Nms) for angular momentum unloading. The calculation results show that in the satellite body coordinate system, the fluctuation range of the X-axis angular momentum is-32.31 Nms~ 6.748 Nms, the fluctuation range of the Y-axis angular momentum is-62.11 Nms~ 25.73 Nms, and the fluctuation range of the Z-axis angular momentum is-31.09 Nms~ 11 Nms.

[0149] According to the calculation results of Figure 6 and Figure 7 , the method described in the embodiment can effectively reduce the fluctuation range of the accumulated angular momentum of the satellite body, and further reduce the fluctuation range of the momentum wheel rotation speed.

[0150] Please refer to Figure 8 , the embodiment of the application provides an angular momentum unloading device based on trajectory prediction, which comprises:

[0151] A calculation module 800 is configured to calculate a combined interference torque of a satellite in a body coordinate system, which is composed of a solar pressure resultant torque and a gravity gradient torque, according to the inertia characteristics, geometric configuration and sun angle of the satellite;

[0152] A prediction module 802 is configured to predict the numerical integral average of the X-axis and Z-axis interference torques and the numerical integral function maximum of the Y-axis interference torque of the satellite at different sun angles according to the combined interference torque, and to calculate the pre-biased angular momentum values of the three axes of the satellite at different sun angles according to the prediction results;

[0153] An unloading module 804 is configured to find a target angular momentum corresponding to the sun angle of the satellite at a preset unloading time according to the pre-biased angular momentum values, and to adjust the jet control torque of the thruster of the satellite according to the finding results, so as to unload the angular momentum of the satellite to the target angular momentum.

[0154] In the embodiment of the present application, the calculation module 800 is specifically configured to perform the following operations when calculating the combined interference torque of the satellite in the body coordinate system formed by the combined solar pressure torque and the gravity gradient torque according to the satellite inertia characteristics, the geometric configuration and the solar angle: calculating a unit azimuth vector of the sun in the satellite orbit system according to the phase angle and the elevation angle of the sun in the satellite orbit system; calculating a single solar pressure torque of each configuration panel of the satellite according to the unit azimuth vector, and calculating the combined solar pressure torque of the satellite according to the single solar pressure torque; wherein the configuration panel includes an east panel, a west panel, a south panel, a north panel, a ground panel, a back ground panel, a south solar wing, a north solar wing, an upper surface of a reflector antenna, a lower surface of the reflector antenna and a side surface of the reflector antenna; calculating the gravity gradient torque of the satellite according to the orbit angular velocity and the moment of inertia of the satellite; and fusing the combined solar pressure torque and the gravity gradient torque to obtain the combined interference torque of the satellite in the body coordinate system.

[0155] In the embodiment of the present application, the prediction module 802 is specifically configured to perform the following operations when predicting the numerical integral average values of the X-axis and Z-axis interference torques of the satellite and the maximum value of the numerical integral function of the Y-axis interference torque of the satellite under different solar angles according to the combined interference torque: equally dividing the solar phase angle and the solar elevation angle to obtain a plurality of solar phase angle discrete values and a plurality of solar elevation angle discrete values in sequence; calculating the numerical integral average values of the X-axis interference torque and the Z-axis interference torque of the satellite in the next orbit period of the satellite in the body coordinate system under different solar elevation angles according to the component model of the combined interference torque along the X, Y and Z axes of the satellite and the numerical integral average values of the Z-axis interference torque

[0156]

[0157] wherein k is the serial number of the solar elevation angle, k = 1, 2, …, n, n is the number of equal divisions of the solar elevation angle; α Sun is the solar phase angle; β Sun is the solar elevation angle; T Bx (α Sun , β Sun ) is the X-axis component model of the combined interference torque; T Bz (α Sun , β Sun ) is the Z-axis component model of the combined interference torque; ω Orb is the orbit angular velocity.

[0158] According to the component model, the maximum value h ymaxjkand the minimum value h yminjk :

[0159]

[0160] wherein T By (α Sun ,β Sun ) is a Y-axis component model of the synthesized interference torque; j is a serial number of the sun phase angle, j = 1, 2, …, m, and m is the number of equal divisions of the sun phase angle.

[0161] In the embodiment of the application, the angular momentum pre-bias values of the three axes of the satellite under different sun angles are calculated by the following formula:

[0162]

[0163] In the formula, H xtk is the pre-bias value of the angular momentum along the X-axis of the satellite body coordinate system when the sun angle is β Sunk ; H ytjk is the pre-bias value of the angular momentum along the Y-axis of the satellite body coordinate system when the sun angle is α Sunj and β Sunk ; and H ztk is the pre-bias value of the angular momentum along the Z-axis of the satellite body coordinate system when the sun angle is β Sunk .

[0164] In the embodiment of the application, when the unloading module finds the target angular momentum corresponding to the sun angle of the satellite at a preset unloading moment according to the angular momentum pre-bias value, it is specifically configured to perform the following operations: the sun phase angle and the sun elevation angle at the unloading moment are calculated according to the unit direction vector of the sun at the unloading moment; the sun phase angle and the sun elevation angle at the unloading moment are calculated by rounding to obtain the discrete serial number of the sun angle corresponding to the unloading moment; and the angular momentum pre-bias value of the three axes is found according to the discrete serial number of the sun angle at the unloading moment to obtain the corresponding target angular momentum.

[0165] In the embodiment of the application, the discrete serial number of the sun angle corresponding to the unloading moment is calculated by the following formula:

[0166]

[0167] In the formula, α SunSat is the sun phase angle at the unloading moment; β SunSat is the sun elevation angle at the unloading moment; int() is the rounding function; j SunSat is the discrete serial number corresponding to the sun phase angle at the unloading moment; and k SunSat is the discrete serial number corresponding to the sun elevation angle at the unloading moment.

[0168] It should be noted that the trajectory prediction-based angular momentum offloading apparatus provided in the above embodiments is only exemplified by the division of the above functional modules, and in actual application, the above functions can be completed by different functional modules according to needs, that is, the internal structure of the apparatus is divided into different functional modules to complete all or part of the functions described above. In addition, the trajectory prediction-based angular momentum offloading apparatus provided in the above embodiments and the trajectory prediction-based angular momentum offloading method embodiments belong to the same concept, and the specific implementation process is detailed in the method embodiments, which will not be described here.

[0169] Embodiments of the present application also provide a computer device, which refers to Figure 9 The computer device includes a processor and a memory, and the memory stores at least one instruction, at least one program, a code set or an instruction set, which is loaded and executed by the processor to implement the trajectory prediction-based angular momentum offloading method provided by each method embodiment.

[0170] Embodiments of the present application also provide a computer readable storage medium, which stores at least one instruction, at least one program, a code set or an instruction set, which is loaded and executed by the processor to implement the trajectory prediction-based angular momentum offloading method provided by each method embodiment.

[0171] Embodiments of the present application also provide a computer program product, which includes a computer program, and the processor of the computer device reads the computer program from the computer readable storage medium, and the processor executes the computer program to enable the computer device to execute the trajectory prediction-based angular momentum offloading method described in any of the above embodiments.

[0172] Finally, it should also be noted that in this document, relational terms such as first, second, third and fourth are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply that there is any such actual relationship or order between these entities or operations. Moreover, the terms "include", "contain" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements not only includes those elements, but also includes other elements not explicitly listed or inherent to such a process, method, article or device. Without more limitations, the element defined by the statement "including a" does not exclude the presence of additional identical elements in the process, method, article or device including the element.

[0173] The above merely describes the preferred embodiments of the present application, and it should be pointed out that, for those skilled in the art, some improvements and refinements can be made without departing from the principles of the present application, and these improvements and refinements should also be considered as the protection scope of the present application.

Claims

1. A method for unloading angular momentum based on trajectory prediction, characterized in that, The method includes: Based on the satellite's inertia characteristics, geometric configuration, and solar angle, calculate the combined disturbance torque of the satellite in its body coordinate system, which consists of the solar radiation pressure torque and the gravity gradient torque. Based on the synthesized interference moment, the numerical integral average of the X-axis and Z-axis interference moments of the satellite at different solar angles is predicted, as well as the maximum and minimum values ​​of the numerical integral function of the Y-axis interference moment, including: By dividing the solar phase angle and solar elevation angle into equal parts, multiple discrete values ​​of the solar phase angle and multiple discrete values ​​of the solar elevation angle are obtained sequentially. Based on the component model of the synthetic interference moment along the X, Y, and Z axes of the satellite, the numerical integral average value of the X-axis interference moment experienced by the satellite in the next orbital period in the body coordinate system under different solar elevation angles is calculated. The numerical integral average of the Z-axis disturbance torque : Where k is the index of the solar elevation angle, k = 1, 2, ..., n, and n is the number of equal divisions of the solar elevation angle; The solar phase angle; The angle of elevation of the sun; The X-axis component model of the synthesized disturbance torque; The Z-axis component model of the synthesized disturbance torque; It is the orbital angular velocity; Based on the component model, the maximum value of the numerical integral function of the Y-axis disturbance torque experienced by the satellite in the next orbital period in the body coordinate system under different solar phase angles and solar elevation angles is calculated. and minimum value : in, The Y-axis component model of the synthesized disturbance torque; j is the index of the solar phase angle, j=1, 2, ..., m, m is the number of equal divisions of the solar phase angle; And based on the prediction results, the angular momentum pre-bias values ​​of the satellite's three axes under different solar angles were calculated; The target angular momentum of the satellite at the solar angle at the preset unloading time is obtained by finding the angular momentum pre-bias value, and the thruster jet control torque of the satellite is adjusted according to the search result to unload the satellite's angular momentum to the target angular momentum.

2. The method as described in claim 1, characterized in that, The calculation of the combined disturbance torque of the satellite in its body coordinate system, consisting of the solar radiation pressure torque and the gravity gradient torque, based on the satellite's inertia characteristics, geometric configuration, and solar angle, includes: Calculate the unit azimuth vector of the sun in the satellite orbit system based on the phase angle and elevation angle of the sun in the satellite orbit system; The single solar pressure torque on each configuration panel of the satellite is calculated based on the unit azimuth vector, and the solar pressure torque on the satellite is calculated based on the single solar pressure torque; wherein, the configuration panels include an east panel, a west panel, a south panel, a north panel, a ground-facing panel, a ground-reflecting panel, a south solar wing, a north solar wing, an upper surface of the reflector antenna, a lower surface of the reflector antenna, and a side surface of the reflector antenna; The gravitational gradient torque acting on the satellite is calculated based on the satellite's orbital angular velocity and moment of inertia. By fusing the solar radiation pressure torque and the gravity gradient torque, the combined interference torque of the satellite in its body coordinate system is obtained.

3. The method as described in claim 1, characterized in that, The angular momentum pre-bias values ​​of the satellite's three axes at different solar angles are calculated using the following formula: In the formula, The value of the solar angle is... Pre-offset value of the angular momentum along the X-axis of the satellite in the satellite body coordinate system; The value of the solar angle is... and Pre-offset value of the angular momentum along the Y-axis of the satellite in the satellite body coordinate system; The value of the solar angle is... The pre-offset value of the angular momentum along the Z-axis of the satellite in the satellite body coordinate system.

4. The method as described in claim 1, characterized in that, The step of finding the target angular momentum of the satellite at the solar angle corresponding to the preset unloading time based on the angular momentum pre-bias value includes: Based on the unit azimuth vector of the sun at the unloading time, the solar phase angle and solar elevation angle at the unloading time are calculated. The solar phase angle and solar elevation angle at the unloading time are rounded to obtain the discrete sequence number of the solar angle corresponding to the unloading time. The target angular momentum is obtained by finding the angular momentum pre-bias value of the three axes based on the discrete index of the solar angle at the unloading time.

5. The method as described in claim 4, characterized in that, The discrete index of the solar angle corresponding to the unloading time is calculated using the following formula: In the formula, The solar phase angle at the moment of unloading; The solar elevation angle at the moment of unloading; int() is the integer function; This is the discrete index corresponding to the solar phase angle at the unloading moment; This is the discrete index corresponding to the solar elevation angle at the unloading time.

6. An angular momentum unloading device based on trajectory prediction, characterized in that, The apparatus, used in the method of any one of claims 1-5, comprises: The calculation module is used to calculate the combined disturbance torque of the satellite in its body coordinate system, which is composed of the solar radiation pressure torque and the gravity gradient torque, based on the satellite's inertial characteristics, geometric configuration and solar angle. The prediction module is used to predict the numerical integral average value of the interference torque of the satellite on the X-axis and Z-axis under different solar angles, as well as the maximum and minimum values ​​of the numerical integral function of the interference torque on the Y-axis, and to calculate the angular momentum pre-bias values ​​of the satellite on the three axes under different solar angles based on the prediction results. The unloading module is used to find the target angular momentum of the satellite at the solar angle at the preset unloading time based on the angular momentum pre-bias value, and adjust the control torque of the satellite according to the search result to unload the satellite angular momentum to the target angular momentum.

7. A computer device, characterized in that, The computer device includes a memory and a processor. The memory is used to store computer programs, and the processor is used to execute the computer programs stored in the memory to implement the steps of the method according to any one of claims 1-5.

8. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, which, when executed by a processor, implements the steps of the method described in any one of claims 1-5.

9. A computer program product, characterized in that, Includes a computer program, which, when executed by a processor, implements the steps of the method according to any one of claims 1-5.

Images