Calculation Method for Collision Avoidance Region of Space Targets in Sun-Synchronous Orbits

Abstract

The method for calculating the collision avoidance area of space targets in a sun-synchronous orbit disclosed by the present invention calculates the areas of various space targets in a sun-synchronous orbit, takes into account the semi-major axis orbit determination deviation and the semi-major axis decay rate caused by atmospheric drag, uses the latest space target catalog data, and iteratively calculates the area of the collision warning of space targets in a sun-synchronous orbit, and finally gives the area for orbit control. The present invention takes into account the orbit perturbation caused by the particularity of the sun-synchronous orbit, and derives the collision avoidance area of space targets in a sun-synchronous orbit. Compared with the traditional calculation of the collision avoidance area of space targets in a sun-synchronous orbit, more accurate results can be obtained.

Description

Technical Field

[0001] The present invention belongs to the technical field of aerospace control methods, and particularly relates to a method for calculating the collision avoidance region of space targets in a sun-synchronous orbit. Background Technique

[0002] On-orbit space targets include spacecraft, space debris, etc. In order to prevent collisions between space targets and spacecraft, which may cause catastrophic consequences, it is necessary to pre-calculate the region during the movement of space targets. The region is a long and narrow shape along the flight direction of the space target, and other space targets outside the region will not collide with this space target within a certain period of time. Traditional calculation of the collision avoidance region of space targets in a sun-synchronous orbit only considers the three-dimensional distance of the space target, ignoring the difference in the collision probability of the normal distance and the in-plane distance of the orbital plane, resulting in inaccurate calculation results. Summary of the Invention

[0003] The purpose of the present invention is to provide a method for calculating the collision avoidance region of space targets in a sun-synchronous orbit, which solves the problem of inaccurate calculation results of existing methods.

[0004] The technical solution adopted by the present invention is: a method for calculating the collision avoidance region of space targets in a sun-synchronous orbit, including the following steps:

[0005] Step 1: Calculate the deviation in the y direction caused by two quantities, namely the semi-major axis orbit determination deviation and the semi-major axis decay rate caused by atmospheric drag, to obtain the region length;

[0006] Step 2: Extrapolate the space target orbit;

[0007] Step 3: Calculate the long-term change rate of the orbital inclination of the space target caused by the gravitational perturbation of the sun and the moon;

[0008] Step 4: Calculate the normal drift of the space target in the sun-synchronous orbit caused by complex perturbing forces;

[0009] Step 5: Calculate the normal radius of the region and the amount by which the normal direction of the region center deviates from the reference orbit to obtain the shape of the collision avoidance region.

[0010] The characteristics of the present invention also lie in that

[0011] Step 1 specifically includes: establishing an orbital coordinate system Cxyz with the centroid of the space target C as the origin, where the direction of the Cx axis is the same as the direction from the Earth's center to the centroid of the space target, the Cz axis points to the normal direction of the space target's orbital plane, and the Cy axis forms a right-handed orthogonal coordinate system with the Cz axis and the Cx axis; let r represent the distance from the space target to the Earth's center, a represent the semi-major axis, e represent the eccentricity, i represent the inclination, Ω represent the right ascension of the ascending node, ω represent the argument of perigee, f represent the true anomaly, M represent the mean anomaly, u = ω + f, λ = ω + M; let the subscripts "c" and "d" represent the space target and the regional boundary point respectively, and the space target is in a nearly circular orbit; let δa0 = a d -a c , δi0 = i d -i c , δΩ0 = Ω d -Ω c , δu0 = u d -u c , δM0 = M d -M c is the maximum value of the orbit determination deviation; δr0 = r d -r c , r c = [r c 0 0] T is the component representation of the position vector of the space target relative to the Earth's center in the space target's orbital coordinate system, r d = [x d y d z d T is the component representation of the position vector of the regional boundary point relative to the Earth's center in the space target's orbital coordinate system; let the orbit determination result be represented by a, e, i, Ω, ω, M, and the orbit determination error be represented by ±δa, ±δe, ±δi, ±δΩ, ±δω, ±δM; for a space target in a low Earth orbit, read the daily decay rate of the semi-major axis in the orbit determination result, da / (dDay), with a negative sign and the unit of meter per day;

[0012] The decay rate of the semi-major axis per second is:

[0013]

[0014] The deviation in the y direction caused by the two quantities of the semi-major axis orbit determination deviation and the semi-major axis decay rate caused by atmospheric drag is calculated according to the following formula:

[0015]

[0016] In formula (2), t is time, T is the orbital period, δa0 is the semi-major axis orbit determination deviation, is the semi-major axis decay rate, π is the pi, y + is the deviation along the +y direction, y​- is the deviation along the -y direction; initial orbit extrapolation, with the extrapolated orbit as the reference benchmark; y + and y - are the two sides of the region; is the region length, is the amount by which the region center deviates from the reference orbit.

[0017] Step 2 specifically includes: using the predictor-corrector formula (3) for orbit extrapolation:

[0018]

[0019] In formula (3), x j is the independent variable at the j-th step, y j is the function value at the j-th step, f is the integral function, h is the step size, is the predicted value, y n+1 is the corrected value.

[0020] In Step 3, the long-term change rate of the space target orbit inclination caused by the sun and moon gravitational perturbations is obtained by summing the long-term change rate of the space target orbit inclination of the sun-synchronous orbit under the action of the sun gravitational perturbation and the long-term change rate of the space target orbit inclination caused by the moon gravitational perturbation;

[0021] The long-term change rate of the space target orbit inclination of the sun-synchronous orbit under the action of the sun gravitational perturbation is calculated by formula (4):

[0022]

[0023] In formula (4), t is the time, i s is the ecliptic inclination, n s is the angular rate of the earth's revolution around the sun, β s is the ecliptic longitude of the sun's apparent motion;

[0024] The long-term change rate of the space target orbit inclination caused by the moon gravitational perturbation is calculated by formula (5):

[0025]

[0026] In formula (5) m m is the moon mass, m e is the earth mass, i m is the inclination of the moon's orbit, Ω m is the right ascension of the ascending node of the moon's orbit, n m is the angular rate of the moon's motion around the earth, r em is the mean earth-moon distance.

[0027] In step 4, the normal drift of the space target in the sun-synchronous orbit is calculated by formula (6):

[0028]

[0029] In formula (6), the unit of calculation is radian.

[0030] Step 5 specifically includes: First, from the normal view, the deviation of the area from the initial reference orbit is:

[0031]

[0032] Then the normal radius of the area is calculated by formula (8):

[0033]

[0034] The amount of deviation of the normal direction of the area center from the reference orbit is calculated by formula (9):

[0035]

[0036] Let K = 2.0 be the area isolation coefficient, d = Kd, the area is appropriately enlarged, and the radius of curvature of the area is

[0037] The physical quantities required for the collision avoidance area include: the amount of deviation z of the area from the initial reference orbit, the normal radius R of the area p , the amount of deviation C of the normal direction of the area center from the reference orbit drift , the length of the area The amount of deviation of the area center from the reference orbit

[0038] The beneficial effects of the present invention are as follows: The method for calculating the collision avoidance area of space targets in the sun-synchronous orbit of the present invention calculates the areas of various space targets in the sun-synchronous orbit, takes into account the semi-major axis orbit determination deviation and the semi-major axis decay rate caused by atmospheric drag, uses the latest space target catalog data, and iteratively calculates the areas for collision warning of space targets in the sun-synchronous orbit, and finally gives the areas for orbit control. The present invention takes into account the orbit perturbation caused by the particularity of the sun-synchronous orbit, and deduces the collision avoidance area of space targets in the sun-synchronous orbit. Compared with the traditional method for calculating the collision avoidance area of space targets in the sun-synchronous orbit, more accurate results can be obtained. Description of the Drawings

[0039] Figure 1 is the flow chart of the method for calculating the collision avoidance area of space targets in the sun-synchronous orbit of the present invention;

[0040] Figure 2Schematic diagram of the collision avoidance region for space targets in a sun-synchronous orbit obtained in an embodiment of the present invention. Detailed implementation manners

[0041] The present invention will be described in detail below with reference to the accompanying drawings and specific implementation manners.

[0042] The present invention provides a method for calculating the collision avoidance region of space targets in a sun-synchronous orbit, as Figure 1 shown, including the following steps:

[0043] Step 1. Calculate the deviation in the y direction caused by two quantities, namely the semi-major axis orbit determination deviation and the semi-major axis decay rate caused by atmospheric drag, to obtain the region length. An orbit coordinate system Cxyz is established with the centroid of space target C as the origin. The direction of the Cx axis is the same as the direction from the earth center to the centroid of the space target, the Cz axis points to the normal direction of the space target orbit plane, and the Cy axis forms a right-handed orthogonal coordinate system with the Cz axis and the Cx axis. Let r represent the distance from the space target to the earth center, a represent the semi-major axis, e represent the eccentricity, i represent the inclination, Ω represent the right ascension of the ascending node, ω represent the argument of perigee, f represent the true anomaly, M represent the mean anomaly, u = ω + f, and λ = ω + M. Let the subscripts "c" and "d" represent the space target and the regional boundary point respectively, and the space target is in a nearly circular orbit. Let δa0 = a d -a c ,δi0 = i d -i c ,δΩ0 = Ω d -Ω c ,δu0 = u d -u c ,δM0 = M d -M c be the maximum value of the orbit determination deviation. δr0 = r d -r c ,r c =[r c 0 0] T be the component representation of the position vector of the space target relative to the earth center in the space target orbit coordinate system, and r d =[x d y d z d T be the component representation of the position vector of the regional boundary point relative to the earth center in the space target orbit coordinate system. Let the orbit determination result be represented by a, e, i, Ω, ω, M, and the orbit determination error be represented by ±δa, ±δe, ±δi, ±δΩ, ±δω, ±δM. For space targets in low-earth orbits, read the daily decay rate of the semi-major axis, da / (dDay), with a negative sign and a unit of meters per day.

[0044] The decay rate of the semi-major axis per second is ​

[0045]

[0046] The deviation in the y direction caused by the two quantities of the semi-major axis orbit determination deviation and the semi-major axis decay rate caused by atmospheric drag is calculated according to the following formula

[0047]

[0048] where t is the time, T is the orbital period, δa0 is the semi-major axis orbit determination deviation, is the semi-major axis decay rate, π is the pi, y + is the deviation along the +y direction, y - is the deviation along the -y direction. For the initial orbit extrapolation, the extrapolated orbit is used as the reference benchmark. y + and y - are the two sides of the region. is the region length, is the amount by which the region center deviates from the reference orbit. Thus, the region length is obtained.

[0049] Step 2. Space target orbit extrapolation. The orbit extrapolation is carried out using the predictor-corrector formula (3):

[0050]

[0051] where, x j is the independent variable at the j-th step, y j is the function value at the j-th step, f is the integral function, h is the step size, is the predicted value, y n+1 is the corrected value.

[0052] The initial orbit is substituted into the instantaneous orbital elements for extrapolation. Considering the 10*10 order (or other orders such as 20*20) gravitational field model, the sun-moon gravity, and the solar radiation pressure, the atmospheric drag is not considered. The extrapolation is carried out with the instantaneous orbital elements as the reference values. Since the deviation of the calculation region center has considered the atmospheric drag, the atmospheric drag is not considered during extrapolation.

[0053] Step 3. Calculate the long-term change rate of the space target orbit inclination caused by the sun-moon gravity perturbation. The long-term change rate of the space target orbit inclination in the sun-synchronous orbit under the action of the sun gravity perturbation caused by the resonance effect is calculated according to the following formula

[0054]

[0055] where, t is the time, i s is the ecliptic inclination, n s is the angular rate of the earth's revolution around the sun, β sis the ecliptic longitude of the sun's apparent motion. At this time, the long-term change rate of the orbital inclination of the constellation space target is independent of the right ascension of the ascending node of the space target and is related to the local time of the descending node of the space target. (β s -Ω) is a quantity determined by the local time of the descending node and has no direct relationship with the right ascension of the ascending node Ω. If the local mean solar time of the descending node of a certain space target in the quasi-sun-synchronous orbit constellation is 12:00, then β s -Ω is equal to 180° after being converted to the interval [0, 360).

[0056] The long-term change rate of the orbital inclination of the space target caused by the lunar gravitational perturbation is calculated according to the following formula

[0057]

[0058] where m m is the mass of the moon, m e is the mass of the earth, i m is the inclination of the lunar orbit, Ω m is the right ascension of the ascending node of the lunar orbit, n m is the angular rate of the moon's motion around the earth, r em is the mean earth-moon distance.

[0059] Step 4. Calculate the normal drift of the sun-synchronous orbit space target caused by complex perturbation forces. Due to the orbit determination deviation, the semi-major axis decay rate caused by atmospheric drag, and the inclination perturbation caused by the sun and moon gravitation, the normal drift of the space target caused by these quantities is calculated according to the following formula:

[0060]

[0061] The unit of calculation is radians.

[0062] Step 5. Calculate the regional normal radius and the amount by which the regional center normal direction deviates from the reference orbit. Looking from the normal direction, the region deviates from the initial reference orbit (extrapolated by the two-body plus J2 term) by

[0063]

[0064] The regional normal radius is calculated according to the following formula

[0065]

[0066] The amount by which the regional center normal direction deviates from the reference orbit is calculated according to the following formula

[0067]

[0068] Output the two sides y + and y - , the regional center Region length |y + -y - |, region radius R p , normal deviation C of the region center from the reference orbit drift ; Let K = 2.0 be the region isolation coefficient, d = Kd, the region is appropriately enlarged, and the radius of curvature of the region is

[0069] Obtain the shape of the collision avoidance region. The physical quantities required for the collision avoidance region include: the amount z by which the region deviates from the initial reference orbit, the normal radius R of the region p , the amount C by which the normal direction of the region center deviates from the reference orbit drift , the region length The amount of tangential deviation of the region center from the reference orbit Among them, the region length and the amount of tangential deviation of the region center from the reference orbit are obtained from step 1; the amount by which the region deviates from the initial reference orbit, the normal radius of the region, and the amount by which the normal direction of the region center deviates from the reference orbit are obtained from steps 4 and 5.

[0070] Embodiment

[0071] The local time of the descending node of a certain sun-synchronous orbit satellite is 6:30 local time. The semi-major axis of the satellite orbit is 7180000 m, the orbit inclination is 98.6 deg, the right ascension of the ascending node at the initial moment is 120 deg, the semi-major axis decay rate is 10 m per day. Assume that the maximum deviation of the semi-major axis for satellite orbit determination is 6 m, the maximum deviation of the inclination for orbit determination is 0.001 deg, the obliquity of the ecliptic is 18.3 deg, the right ascension of the ascending node of the lunar orbit is -12 deg, and the region isolation coefficient is 2.0. Calculate the size, shape, and region center position of the safe region after 1 day.

[0072] According to the above calculation steps, the length of the collision avoidance region is 1.6139 km, and the length direction of the collision avoidance region is along the satellite running direction. The normal radius of the collision avoidance region is 27.8045 m, and the normal drift of the region is -0.1570 m (facing the satellite running direction, positive to the right and negative to the left). Figure 2 The schematic diagram of the safe region after 1 day is given. For the sake of visualization, Figure 2 the scale sizes of different coordinate axes are inconsistent. From Figure 2 it can be seen that the collision avoidance region is a very slender cylinder.

Claims

1. A method for calculating the collision avoidance area of space targets in a sun-synchronous orbit, characterized in that, It includes the following steps: Step 1: Calculate the deviation in the y direction caused by the two quantities of the semi-major axis orbit determination deviation and the semi-major axis decay rate caused by atmospheric drag, and obtain the regional length; Step 2: Space target orbit extrapolation; Step 3: Calculate the long-term change rate of the space target orbit inclination angle caused by the sun and moon gravitational perturbations; Step 4: Calculate the normal drift of the sun-synchronous orbit space target caused by complex perturbation forces; Step 5: Calculate the regional normal radius and the amount by which the regional center normal direction deviates from the reference orbit, and obtain the shape of the collision avoidance area.

2. The method for calculating the collision avoidance area of space targets in a sun-synchronous orbit according to claim 1, characterized in that, The specific steps of Step 1 include: establishing an orbital coordinate system Cxyz with the centroid of the space target C as the origin, where the direction of the Cx axis is the same as the direction from the earth's center to the centroid of the space target, the Cz axis points to the normal direction of the space target's orbital plane, and the Cy axis forms a right-handed orthogonal coordinate system with the Cz axis and the Cx axis; let r represent the distance from the space target to the earth's center, a represent the semi-major axis, e represent the eccentricity, i represent the inclination, Ω represent the right ascension of the ascending node, ω represent the argument of perigee, f represent the true anomaly, M represent the mean anomaly, u = ω + f, and λ = ω + M; let the subscripts "c" and "d" represent the space target and the regional boundary point respectively, and the space target is in a nearly circular orbit; let δa0 = a d -a c , δi0 = i d -i c , δΩ0 = Ω d -Ω c , δu0 = u d -u c , δM0 = M d -M c as the maximum value of the orbit determination deviation; δr0 = r d -r c , r c = [r c 0 0] T is the component representation of the position vector of the space target relative to the earth's center in the space target's orbital coordinate system, and r d = [x d y d z d T is the component representation of the position vector of the regional boundary point relative to the earth's center in the space target's orbital coordinate system; let the orbit determination result be represented by a, e, i, Ω, ω, M, and the orbit determination error be represented by ±δa, ±δe, ±δi, ±δΩ, ±δω, ±δM; for a space target in a low-earth orbit, read the daily decay rate of the semi-major axis in the orbit determination result, da / (dDay), with a negative sign and the unit of meter per day;​ The semi-major axis decay rate per second is: The deviation in the y direction caused by the two quantities of the semi-major axis orbit determination deviation and the semi-major axis decay rate caused by atmospheric drag is calculated according to the following formula: In Equation (2), t is time, T is the orbital period, δa0 is the deviation of the semi-major axis for orbit determination, is the semi-major axis decay rate, π is pi, y + is the deviation along the +y direction, y - is the deviation along the -y direction; for the initial orbit extrapolation, the extrapolated orbit is used as the reference benchmark; y + and y - are the two sides of the region; is the region length, is the amount by which the region center deviates from the reference orbit.

3. The method for calculating the collision avoidance area of space targets in a sun-synchronous orbit according to claim 2, characterized in that, The specific content of Step 2 includes: Using the predictor-corrector formula (3) for orbit extrapolation: In formula (3), x j is the independent variable at the j-th step, y j is the function value at the j-th step, f is the integral function, h is the step size, is the predicted value, and y n+1 is the corrected value.

4. The method for calculating the collision avoidance area of space targets in a sun-synchronous orbit according to claim 3, characterized in that, In Step 3, the long-term change rate of the space target orbit inclination angle caused by the sun and moon gravitational perturbations is obtained by summing the long-term change rate of the sun-synchronous orbit space target orbit inclination angle under the action of the sun gravitational perturbation and the long-term change rate of the space target orbit inclination angle caused by the moon gravitational perturbation; The long-term change rate of the sun-synchronous orbit space target orbit inclination angle under the action of the sun gravitational perturbation is calculated by formula (4): In Equation (4), t is time, i s is the obliquity of the ecliptic, n s is the angular rate of the Earth's revolution around the Sun, β s is the ecliptic longitude of the apparent solar motion; The long-term change rate of the space target orbit inclination angle caused by the moon gravitational perturbation is calculated by formula (5): In Equation (5) m m is the mass of the moon, m e is the mass of the earth, i m is the inclination of the lunar orbit, Ω m is the right ascension of the ascending node of the lunar orbit, n m is the angular rate of the moon's motion around the earth, r em is the average distance between the earth and the moon.

5. The method for calculating the collision avoidance area of space targets in a sun-synchronous orbit according to claim 4, characterized in that, In Step 4, the normal drift of the sun-synchronous orbit space target is calculated by formula (6): In formula (6), the calculation unit is radians.

6. The method for calculating the collision avoidance area of space targets in a sun-synchronous orbit according to claim 5, characterized in that, The specific content of Step 5 includes: First, looking from the normal direction, the region deviates from the initial reference orbit by: Then the regional normal radius is calculated by formula (8): The amount by which the regional center normal direction deviates from the reference orbit is calculated by formula (9): Let K = 2.0 be the regional isolation coefficient, d = Kd, the region is appropriately enlarged, and the radius of curvature of the region is The physical quantities required for the collision avoidance region include: the amount z by which the region deviates from the initial reference orbit, and the normal radius R of the region p , the amount C by which the normal direction of the region center deviates from the reference orbit drift , the length of the region The amount by which the region center deviates from the reference orbit

Images