Calculation method of attitude perturbation of space debris

Abstract

The space debris attitude perturbation calculation method disclosed by the present invention includes calculating the aerodynamic force and torque on the space debris, and then calculating the debris attitude precession angular velocity under the action of the aerodynamic torque; calculating the gravity gradient torque on the space debris, and then calculating the debris attitude precession angular velocity under the action of the gravity gradient torque; calculating the geomagnetic torque on the space debris, and then calculating the debris attitude precession angular velocity under the action of the geomagnetic torque; calculating the solar pressure torque on the space debris, and then calculating the debris attitude precession angular velocity and the debris attitude nutation angular velocity under the action of the solar pressure torque; adding the obtained debris attitude precession angular velocity and nutation angular velocity, and integrating to obtain the space debris attitude orientation composed of the precession angle, nutation angle, and spin angle. The space debris attitude perturbation calculation method of the present invention solves the problem that the existing space debris attitude perturbation calculation method cannot directly obtain the attitude orientation, precession angle, and nutation angle of the space debris.

Description

Technical Field

[0001] The invention belongs to the technical field of aerospace engineering, and in particular relates to a space debris attitude perturbation calculation method. Background Art

[0002] Space debris is generally in a spinning state. Under the action of attitude perturbation torque, the inertial space direction of the spin axis, spin angular rate, precession and nutation of the space debris will change. For space debris in low-Earth orbit, aerodynamic torque, gravity gradient torque, geomagnetic torque and solar pressure torque all play an important role. For space debris in geosynchronous orbit, aerodynamic torque can be ignored, and solar pressure torque, gravity gradient torque and geomagnetic torque play the main role. The existing space debris attitude perturbation calculation, by substituting the space debris attitude torque into the momentum moment equation, solves the variable as momentum moment, and cannot directly obtain the attitude direction, precession angle and nutation angle of the space debris. In addition, the integral calculation is large and the intuitiveness is low. Summary of the invention

[0003] The purpose of the present invention is to provide a space debris attitude perturbation calculation method, which solves the problem that the existing space debris attitude perturbation calculation method cannot directly obtain the attitude pointing, precession angle and nutation angle of the space debris.

[0004] The technical solution adopted by the present invention is: a space debris attitude perturbation calculation method, comprising the following steps:

[0005] Step 1: Calculate the aerodynamic force and torque on the space debris, and bring them into the attitude perturbation equation to calculate the attitude precession angular velocity of the debris under the action of the aerodynamic torque;

[0006] Step 2: Calculate the gravity gradient torque on the space debris, and bring it into the attitude perturbation equation to calculate the attitude precession angular velocity of the debris under the action of the gravity gradient torque;

[0007] Step 3: Calculate the geomagnetic torque on the space debris, and bring it into the attitude perturbation equation to calculate the attitude precession angular velocity of the debris under the action of the geomagnetic torque;

[0008] Step 4, calculate the solar pressure moment on the space debris, and bring it into the attitude perturbation equation to calculate the debris attitude precession angular velocity and debris attitude nutation angular velocity under the action of the solar pressure moment;

[0009] Step 5: Add the debris attitude precession angular velocity obtained in steps 1 to 4 and the debris attitude nutation angular velocity obtained in step 4, and integrate to obtain the space debris attitude direction composed of the precession angle, nutation angle, and spin angle.

[0010] The present invention is also characterized in that:

[0011] In step 1, the aerodynamic force on the space debris is calculated by formula (1):

[0012]

[0013] In formula (1), c d is the aerodynamic coefficient, which is between 2.1 and 2.2, ρ a is the atmospheric density, S is the headwind area, v a is the average velocity vector of the incident airflow relative to the headwind area element, v a The direction is the flying direction of the debris, v a is the velocity vector v a size;

[0014] The aerodynamic torque on space debris is calculated by formula (2):

[0015]

[0016] In formula (2), r a The radius vector from the debris center of mass to the aerodynamic pressure center, r a =[r ax r ay r az ] T , v a =[v ax v ay v az ] T .

[0017] In step 1, the angular velocity of the debris attitude precession under the action of aerodynamic torque is calculated by formula (3):

[0018]

[0019] In formula (3), N is the angular velocity of the debris attitude precession under the action of aerodynamic torque. The coordinate system Oξηζ of the debris's center of mass is defined as rotating around the η axis by an angle of σ, and then rotating around the new position of the ξ axis by an angle of ρ, so that the η axis points to the direction of the momentum vector H. The atmospheric density at the perigee of the debris orbit is ρ π , the atmospheric density at the orbit with true anomaly f is ρ π h(f), l is the force arm of the aerodynamic moment, e is the eccentricity, and p is the semi-diameter. Under the action of the aerodynamic moment, the spin axis points to the η axis of the perigee coordinate system and performs forced regular precession, satisfying ρ=ρ0;μ=398600.44km 3 / s 2 is the Earth's gravitational constant.

[0020] In step 2, the gravity gradient torque on the space debris is calculated by formula (4):

[0021]

[0022] In formula (4), μ is the gravitational constant of the earth, r is the distance from the debris to the center of the earth, r is the vector from the center of the earth to the debris, E = -r / r, and E is expressed in the reference coordinate system as E = [E x E y E z ] T , I is the inertia matrix of the fragment relative to the center of mass, is the adjoint matrix of E,

[0023] Assume T g In the reference coordinate system, it is represented as T g =[T gx T gy T gz ] T ,but:

[0024]

[0025] If the reference coordinate system is chosen as the fragment inertial principal axis system, then I xy =I yz =I zx =0, then the gravity gradient torque is calculated as follows:

[0026]

[0027] In step 2, the angular velocity of the debris attitude precession under the action of the gravity gradient torque is calculated by formula (7):

[0028]

[0029] In formula (7), the perigee coordinate system Oξηζ of the center of mass rotates by an angle σ around the ζ axis, and then rotates by an angle ρ around the new position of the η axis. The ζ axis points to the direction of the momentum H. Under the action of the gravity gradient torque, the spin axis points to the normal of the orbital plane for forced regular precession. The angular velocity is N, ω c is the angular velocity of the Earth's rotation, σ0 and ρ0 are the values ​​of σ and ρ when f = 0, i.e. when the debris moves to the perigee, satisfying ρ=ρ0; the inertia of the star along the three principal axes of inertia are A, B and C respectively, and ψ is the precession angle.

[0030] In step 3, the geomagnetic torque on the space debris is calculated by formula (8):

[0031] T m =M m ×B (8)

[0032] In formula (8), the unit of geomagnetic moment is N·m, M m is the fragment magnetic moment, in A·m 2 ; B is the magnetic field strength, unit is Wb / m 2 , the magnetic field strength B is calculated according to formulas (9) and (10):

[0033]

[0034]

[0035] The average geomagnetic reference radius is a = 6371.2 km, θ and is the latitude and longitude, t is the time, n and m are degrees and orders, and Gaussian coefficient of the Earth's magnetic field, is the normalized Legendre function, N gm is the maximum geomagnetic field model order used, and r is the distance from the fragment to the center of the Earth.

[0036] In step 3, the angular velocity of the debris attitude precession under the action of the geomagnetic torque is calculated by formula (11):

[0037]

[0038] Calculate the location of the singularity:

[0039]

[0040]

[0041] Determine the Z through the singular points S1 and S2 s axis:

[0042] Z s = sinρ s Y I +cosρ s Z I (14)

[0043] In formula (14), Under the influence of the geomagnetic moment, the fragment spin axis points to the average position Z around the geomagnetic field strength B. s The axis precesses, and the precession angular velocity is N; in formula (11), μ m is the Earth's magnetic moment constant, μ is the Earth's gravitational constant, I m is the modulus of the star's magnetic moment, i is the orbital inclination, H is the magnitude of the star's angular momentum, and p is the semi-diameter.

[0044] In step 4, the solar pressure moment on the space debris is calculated by formula (15):

[0045] Ts =r s ×F s (15)

[0046] In formula (15), r s Fs is the radius vector from the debris mass center to the light pressure center, and is the sunlight pressure on the space debris, which is calculated by formula (16):

[0047]

[0048] In formula (16), n and t are the unit vectors of the normal and tangent directions of the illuminated surface, respectively; I is the solar constant; c is the speed of light; S is the effective illuminated surface area; ν is the surface reflection coefficient; σ is the angle of incidence, i.e., the angle between the incident light and n.

[0049] In step 4, the debris attitude precession angular velocity and debris attitude nutation angular velocity under the action of the solar pressure moment are calculated by formulas (17) and (18) respectively:

[0050]

[0051]

[0052] In formula (17), is the angular velocity of the fragment spin axis pointing to the forced regular precession around the ζ axis of the perigee coordinate system, that is, the angular velocity of the fragment spin axis pointing to the forced regular precession around the orbital normal vector; is the unit vector pointing from the center of the earth to the sun, is the unit vector of the momentum vector direction, c L is the distance from the center of mass to the center of pressure, H is the magnitude of the momentum vector, e is the orbital eccentricity, f is the true anomaly, p is the semi-diameter, and μ is the earth's gravitational constant; suppose the perigee coordinate system Oξηζ of the center of mass rotates by an angle σ around the ζ axis and then rotates by an angle ρ around the new position of the η axis, then the ζ axis points in the direction of the momentum H; the upper and lower areas of the fragment are S1, and the surrounding area is S2.

[0053] Step 5 specifically includes: integrating the kinematic equation of the space debris principal inertia coordinate system relative to the mass center inertia coordinate system through formula (19):

[0054]

[0055] In formula (19), ψ is the precession angle, θ is the nutation angle, is the spin angle, and the attitude dynamics equation is integrated using formula (20):

[0056]

[0057] The component of the unit vector of the spin axis of space debris in the geocentric inertial coordinate system is expressed by formula (21):

[0058]

[0059] The right ascension α and declination δ of the spin axis of the space debris attitude rotation and their change rate are calculated by formula (22):

[0060]

[0061] The beneficial effect of the present invention is that the space debris attitude perturbation calculation method of the present invention can directly calculate the attitude pointing, precession angle and nutation angle of the space debris, which is helpful to estimate and calculate the debris attitude in advance during the process of capturing and removing space debris. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 It is a schematic flow chart of the space debris attitude perturbation calculation method of the present invention;

[0063] Figure 2 is a diagram showing the change of gravity gradient moment of each axis of the principal inertia coordinate system of an embodiment of the present invention;

[0064] Figure 3 is a diagram showing changes in the angular velocity of each axis of the principal inertia coordinate system of an embodiment of the present invention;

[0065] Figure 4 is a diagram of the change of the precession angle of an embodiment of the present invention;

[0066] Figure 5 is a diagram showing a gradual increase in the nutation angle of an embodiment of the present invention;

[0067] Figure 6 is a diagram showing a gradual decrease in the nutation angle according to an embodiment of the present invention;

[0068] Figure 7 2 is a diagram of right ascension change according to an embodiment of the present invention. DETAILED DESCRIPTION

[0069] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0070] The present invention provides a space debris attitude perturbation calculation method, such as Figure 1 As shown, the following steps are included:

[0071] Step 1: Calculate the aerodynamic force and torque on the space debris, and then calculate the attitude precession angular velocity of the debris under the action of the aerodynamic torque;

[0072] Specifically, the aerodynamic force on space debris is calculated by formula (1):

[0073]

[0074] In formula (1), c d is the aerodynamic coefficient, which is between 2.1 and 2.2, ρ a is the atmospheric density, S is the headwind area, v a is the average velocity vector of the incident airflow relative to the headwind area element, v a The direction is the flying direction of the debris, v a is the velocity vector v a size.

[0075] The aerodynamic torque on space debris is calculated by formula (2):

[0076]

[0077] In formula (2), r a The radius vector from the debris center of mass to the aerodynamic pressure center, r a =[r ax r ay r az ] T , v a =[v ax v ay v az ] T .

[0078] The angular velocity of the debris attitude precession under the action of aerodynamic torque is calculated by formula (3):

[0079]

[0080] In formula (3), N is the angular velocity of the debris attitude precession under the action of aerodynamic torque. The coordinate system Oξηζ of the debris's center of mass is defined as rotating around the η axis by an angle of σ, and then rotating around the new position of the ξ axis by an angle of ρ, so that the η axis points to the direction of the momentum vector H. The atmospheric density at the perigee of the debris orbit is ρ π , the atmospheric density at the orbit with true anomaly f is ρ π h(f), l is the force arm of the aerodynamic moment, e is the eccentricity, and p is the semi-diameter. Under the action of the aerodynamic moment, the spin axis points to the η axis of the perigee coordinate system and performs forced regular precession, satisfying ρ=ρ0;μ=398600.44km 3 / s 2 is the Earth's gravitational constant.

[0081] Step 2: Calculate the gravity gradient torque on the space debris, and then calculate the attitude precession angular velocity of the debris under the gravity gradient torque;

[0082] Specifically, the gravity gradient torque on the space debris is calculated by formula (4):

[0083]

[0084] In formula (4), μ is the gravitational constant of the earth, r is the distance from the debris to the center of the earth, r is the vector from the center of the earth to the debris, E = -r / r, and E is expressed in the reference coordinate system as E = [E x E y E z ] T , I is the inertia matrix of the fragment relative to the center of mass, is the adjoint matrix of E,

[0085] Assume T g In the reference coordinate system, it is represented as T g =[T gx T gy T gz ] T ,but:

[0086]

[0087] If the reference coordinate system is chosen as the fragment inertial principal axis system, then I xy =I yz =I zx =0, then the gravity gradient torque is calculated as follows:

[0088]

[0089] The angular velocity of the debris attitude precession under the action of the gravity gradient torque is calculated by formula (7):

[0090]

[0091] In formula (7), the perigee coordinate system Oξηζ of the center of mass rotates by an angle σ around the ζ axis, and then rotates by an angle ρ around the new position of the η axis. The ζ axis points to the direction of the momentum H. Under the action of the gravity gradient torque, the spin axis points to the normal of the orbital plane for forced regular precession. The angular velocity is N, ω c is the angular velocity of the Earth's rotation, σ0 and ρ0 are the values ​​of σ and ρ when f = 0, i.e. when the debris moves to the perigee, satisfying ρ=ρ0; the inertia of the star along the three principal axes of inertia are A, B and C respectively, and ψ is the precession angle.

[0092] Step 3, calculating the geomagnetic torque exerted on the space debris, and then calculating the attitude precession angular velocity of the debris under the action of the geomagnetic torque;

[0093] Specifically, the geomagnetic torque exerted on space debris is calculated by formula (8):

[0094] T m =M m ×B (8)

[0095] In formula (8), the unit of geomagnetic moment is N·m, M m is the fragment magnetic moment, in A·m 2 ; B is the magnetic field strength, unit is Wb / m 2 , the magnetic field strength B is calculated according to formulas (9) and (10):

[0096]

[0097]

[0098] The average geomagnetic reference radius is a = 6371.2 km, θ and is the latitude and longitude, t is the time, n and m are degrees and orders, and Gaussian coefficient of the Earth's magnetic field, is the normalized Legendre function, N gm is the maximum geomagnetic field model order used, and r is the distance from the fragment to the center of the Earth.

[0099] The angular velocity of the debris attitude precession under the action of the geomagnetic torque is calculated by formula (11):

[0100]

[0101] Calculate the location of the singularity:

[0102]

[0103]

[0104] Determine the Z through the singular points S1 and S2 s axis:

[0105] Z s = sinρ s Y I +cosρ s Z I (14)

[0106] In formula (14), Under the influence of the geomagnetic moment, the fragment spin axis points to the average position Z around the geomagnetic field strength B. s The axis precesses, and the precession angular velocity is N; in formula (11), μ m is the Earth's magnetic moment constant, μ is the Earth's gravitational constant, I m is the modulus of the star's magnetic moment, i is the orbital inclination, H is the magnitude of the star's angular momentum, and p is the semi-diameter.

[0107] Step 4, calculating the solar pressure moment exerted on the space debris, and then calculating the debris attitude precession angular velocity and debris attitude nutation angular velocity under the action of the solar pressure moment;

[0108] Specifically, the solar pressure moment exerted on space debris is calculated by formula (15):

[0109] T s =r s ×F s (15)

[0110] In formula (15), r s Fs is the radius vector from the debris mass center to the light pressure center, and is the sunlight pressure on the space debris, which is calculated by formula (16):

[0111]

[0112] In formula (16), n and t are the unit vectors of the normal and tangent directions of the illuminated surface, respectively; I is the solar constant; c is the speed of light; S is the effective illuminated surface area; ν is the surface reflection coefficient; σ is the angle of incidence, i.e., the angle between the incident light and n.

[0113] The debris attitude precession angular velocity and debris attitude nutation angular velocity under the action of the solar pressure moment are calculated by formulas (17) and (18) respectively:

[0114]

[0115]

[0116] In formula (17), is the angular velocity of the fragment spin axis pointing to the forced regular precession around the ζ axis of the perigee coordinate system, that is, the angular velocity of the fragment spin axis pointing to the forced regular precession around the orbital normal vector; is the unit vector pointing from the center of the earth to the sun, is the unit vector of the momentum vector direction, c L is the distance from the center of mass to the center of pressure, H is the magnitude of the momentum vector, e is the orbital eccentricity, f is the true anomaly, p is the semi-diameter, and μ is the earth's gravitational constant; suppose the perigee coordinate system Oξηζ of the center of mass rotates by an angle σ around the ζ axis and then rotates by an angle ρ around the new position of the η axis, then the ζ axis points in the direction of the momentum H; the upper and lower areas of the fragment are S1, and the surrounding area is S2.

[0117] Step 5, adding the debris attitude precession angular velocity obtained in steps 1 to 4 and the debris attitude nutation angular velocity obtained in step 4, integrating to obtain the space debris attitude orientation composed of the precession angle, nutation angle, and spin angle, specifically including:

[0118] The kinematic equation of the principal inertia coordinate system of the space debris relative to the mass center inertia coordinate system is integrated by formula (19):

[0119]

[0120] In formula (19), ψ is the precession angle, θ is the nutation angle, is the spin angle, and the attitude dynamics equation is integrated using formula (20):

[0121]

[0122] The component of the unit vector of the spin axis of space debris in the geocentric inertial coordinate system is expressed by formula (21):

[0123]

[0124] The right ascension α and declination δ of the spin axis of the space debris attitude rotation and their rate of change are calculated by formula (22) (radians):

[0125]

[0126] Example

[0127] The following initial values ​​are used for calculation

[0128]

[0129] Extrapolate the gravity gradient moment, attitude angular velocity, precession angle, nutation angle, right ascension and declination of the spin axis for 2 days, with a step length of 0.5 seconds, output a set of data every minute, and draw a graph of its changing law. The gravity gradient moment of each axis of the main inertia coordinate system is shown in Figure 2 , the angular velocity of each axis of the principal inertia coordinate system is Figure 3 , the precession angle changes basically linearly. Figure 4 The calculation results show that the nutation angle basically increases gradually. Figure 5 , and the change rate of the nutation angle is much smaller than that of the precession angle. In two days, the change of the precession angle is about 2.3×10 4 The increase in nutation angle is due to the fact that the satellite rotates around the axis of minimum principal inertia at the initial moment. z =1300kg·m 2 , the other initial values ​​remain unchanged, the satellite rotates around the maximum principal inertia axis, and the nutation angle basically decreases gradually, such as Figure 6The linear change of right ascension is shown in Figure 7 , the declination gradually decreases, and the speed of declination change is much smaller than that of right ascension. In 2 days, the change of right ascension is about 2.3×10 4 Compared with the traditional method, this method has high calculation accuracy and can directly calculate and output the attitude pointing, precession angle and nutation angle of space debris.

Claims

1. A space debris attitude perturbation calculation method, characterized in that: The following steps are involved: Step 1: Calculate the aerodynamic force and torque on the space debris, and bring them into the attitude perturbation equation to calculate the attitude precession angular velocity of the debris under the action of the aerodynamic torque; The aerodynamic force on space debris is calculated by formula (1): In formula (1), c d is the aerodynamic coefficient, which is between 2.1 and 2.2, ρ a is the atmospheric density, S is the headwind area, v a is the average velocity vector of the incident airflow relative to the headwind area element, v a The direction is the flying direction of the debris, v a is the velocity vector v a size; The aerodynamic torque on space debris is calculated by formula (2): In formula (2), r a The radius vector from the debris center of mass to the aerodynamic pressure center, r a =[r ax r ay r az ] T , v a =[v ax v ay v az ] T ; The angular velocity of the debris attitude precession under the action of aerodynamic torque is calculated by formula (3): In formula (3), N is the angular velocity of the debris attitude precession under the action of aerodynamic torque. The coordinate system Oξηζ of the debris's center of mass is defined as rotating around the η axis by an angle of σ, and then rotating around the new position of the ξ axis by an angle of ρ, so that the η axis points to the direction of the momentum vector H. The atmospheric density at the perigee of the debris orbit is ρ π , the atmospheric density at the orbit with true anomaly f is ρ π h(f), l is the force arm of the aerodynamic moment, e is the eccentricity, and p is the semi-diameter. Under the action of the aerodynamic moment, the spin axis points to the η axis of the perigee coordinate system and performs forced regular precession, satisfying ρ=ρ0;μ=398600.44km 3 / s 2 is the Earth's gravitational constant; Step 2: Calculate the gravity gradient torque on the space debris, and bring it into the attitude perturbation equation to calculate the attitude precession angular velocity of the debris under the action of the gravity gradient torque; The gravity gradient torque on space debris is calculated by formula (4): In formula (4), μ is the gravitational constant of the earth, r is the distance from the debris to the center of the earth, r is the vector from the center of the earth to the debris, E = -r / r, and E is expressed in the reference coordinate system as E = [E x E y E z ] T , I is the inertia matrix of the fragment relative to the center of mass, is the adjoint matrix of E, Assume T g In the reference coordinate system, it is represented as T g =[T gx T gy T gz ] T ,but: If the reference coordinate system is chosen as the fragment inertial principal axis system, then I xy =I yz =I zx =0, then the gravity gradient torque is calculated as follows: The angular velocity of the debris attitude precession under the action of the gravity gradient torque is calculated by formula (7): In formula (7), the perigee coordinate system Oξηζ of the center of mass rotates by an angle σ around the ζ axis, and then rotates by an angle ρ around the new position of the η axis. The ζ axis points to the direction of the momentum H. Under the action of the gravity gradient torque, the spin axis points to the normal of the orbital plane for forced regular precession. The angular velocity is N, ω c is the angular velocity of the Earth's rotation, σ0 and ρ0 are the values ​​of σ and ρ when f = 0, i.e. when the debris moves to the perigee, satisfying ρ=ρ0; the inertia of the star along the three main inertia axes are A, B and C respectively, and ψ is the precession angle; Step 3: Calculate the geomagnetic torque on the space debris, and bring it into the attitude perturbation equation to calculate the attitude precession angular velocity of the debris under the action of the geomagnetic torque; The geomagnetic torque on space debris is calculated by formula (8): T m =M m ×B(8) In formula (8), the unit of geomagnetic moment is N·m, M m is the fragment magnetic moment, in A·m 2 ; B is the magnetic field strength, unit is Wb / m 2 , the magnetic field strength B is calculated according to formulas (9) and (10): The average geomagnetic reference radius is a = 6371.2 km, θ and is the latitude and longitude, t is the time, n and m are degrees and orders, and Gaussian coefficient of the Earth's magnetic field, is the normalized Legendre function, N gm is the maximum geomagnetic field model order used, r is the distance from the fragment to the center of the Earth; The angular velocity of the debris attitude precession under the action of the geomagnetic torque is calculated by formula (11): Calculate the location of the singularity: Determine the Z through the singular points S1 and S2 s axis: Z s =sinρ s Y I +cosρ s Z I (14) In formula (14), Under the influence of the geomagnetic moment, the fragment spin axis points to the average position Z around the geomagnetic field strength B. s The axis precesses, and the precession angular velocity is N; in formula (11), μ m is the Earth's magnetic moment constant, μ is the Earth's gravitational constant, I m is the magnitude of the star's magnetic moment, i is the orbital inclination, H is the magnitude of the star's momentum moment, and p is the semi-diameter; Step 4, calculate the solar pressure moment on the space debris, and bring it into the attitude perturbation equation to calculate the debris attitude precession angular velocity and debris attitude nutation angular velocity under the action of the solar pressure moment; The solar pressure moment on space debris is calculated by formula (15): T s =r s ×F s (15) In formula (15), r s Fs is the radius vector from the debris mass center to the light pressure center, and is the sunlight pressure on the space debris, which is calculated by formula (16): In formula (16), n and t are the unit vectors of the normal and tangent directions of the illuminated surface, respectively; I is the solar constant; c is the speed of light; S is the effective illuminated surface area; ν is the surface reflection coefficient; σ is the angle of incidence, i.e., the angle between the incident light and n. The debris attitude precession angular velocity and debris attitude nutation angular velocity under the action of sunlight pressure moment are calculated by formula (17) and (18): In formula (17), is the angular velocity of the fragment spin axis pointing to the forced regular precession around the ζ axis of the perigee coordinate system, that is, the angular velocity of the fragment spin axis pointing to the forced regular precession around the orbital normal vector; is the unit vector pointing from the center of the earth to the sun, is the unit vector of the momentum vector direction, c L is the distance from the center of mass to the center of pressure, H is the magnitude of the momentum vector, e is the orbital eccentricity, f is the true anomaly angle, p is the semi-diameter, and μ is the Earth's gravitational constant; suppose the perigee coordinate system Oξηζ of the center of mass rotates by an angle σ around the ζ axis and then rotates by an angle ρ around the new position of the η axis, so that the ζ axis points in the direction of the momentum H; the upper and lower areas of the fragment are S1, and the surrounding area is S2; Step 5, adding the debris attitude precession angular velocity obtained in steps 1 to 4 and the debris attitude nutation angular velocity obtained in step 4, and integrating to obtain the space debris attitude orientation composed of the precession angle, the nutation angle, and the spin angle; The kinematic equation of the principal inertia coordinate system of the space debris relative to the mass center inertia coordinate system is integrated by formula (19): In formula (19), ψ is the precession angle, θ is the nutation angle, is the spin angle, and the attitude dynamics equation is integrated using formula (20): The component of the unit vector of the spin axis of space debris in the geocentric inertial coordinate system is expressed by formula (21): The right ascension α and declination δ of the spin axis of the space debris attitude rotation and their change rate are calculated by formula (22):

Images