A satellite celestial navigation calculation method using attitude sensor

Abstract

The present invention proposes a satellite astronomical navigation calculation method using attitude sensors. The method of the present invention uses attitude sensors on satellites: star sensors and sun sensors to perform satellite autonomous navigation, without the need to add other hardware devices, thereby reducing the cost and weight of the satellite navigation system. In particular, in the past decade or so, with the development of electronic technology, the measurement accuracy of attitude sensors has been greatly improved, and using star sensors and sun sensors to autonomously determine the orbital information of satellites has great engineering significance.

Description

Technical Field

[0001] The invention belongs to the field of aerospace technology, is applied to autonomous navigation of satellites, uses star sensors and sun sensors to calculate the position of satellites, and particularly relates to a satellite astronomical navigation calculation method using an attitude sensor. Background Art

[0002] Navigation and positioning technology plays an increasingly important role in various civil and military fields. Low-orbit micro-satellites require orbit determination systems that are light in weight, low in power consumption, highly reliable, and low in cost. However, the orbit determination and prediction of satellites currently generally rely on the support of ground measurement and control, and are not very autonomous. Some satellites achieve autonomous orbit determination by carrying GPS receivers, but the cost is high and there is a problem of human interference. Among the commonly used navigation technologies, there are many shortcomings in the use of a single technology: the inertial navigation system has high accuracy in a short period of time, but due to the cumulative error of inertial devices, it is not suitable for long-term navigation; radio navigation equipment is simple, reliable, and highly accurate, but the signal is easily shielded and interfered, and the navigation base station itself is also vulnerable to attack. Therefore, it is very necessary to study efficient and low-cost autonomous orbit determination technology.

[0003] Astronomical navigation is an autonomous navigation method that can obtain real-time and accurate navigation information by observing celestial information. In recent years, astronomical navigation technology that uses solar and lunar information for autonomous satellite orbit determination has received great attention. The sun is the natural celestial body most familiar to humans, with the clearest movement laws and the most prominent and easily identifiable optical features, which can be used for satellite position navigation.

[0004] The present invention utilizes the attitude sensors on the satellite: star sensors and sun sensors for autonomous satellite navigation, without the need to add other hardware devices, thereby reducing the cost and weight of the satellite navigation system. In particular, in the past decade or so, with the development of electronic technology, the measurement accuracy of attitude sensors has been greatly improved, and using star sensors and sun sensors to autonomously determine the orbit information of satellites has great engineering significance. Summary of the invention

[0005] The purpose of the present invention is to solve the problems in the prior art and propose a satellite celestial navigation calculation method using an attitude sensor.

[0006] The present invention is implemented by the following technical scheme. The present invention proposes a satellite celestial navigation calculation method using an attitude sensor, and the method comprises the following steps:

[0007] Step 1: Describe the spatial position of the satellite;

[0008] Step 2: Determine the satellite orbit period T, semi-major axis a and inclination i;

[0009] Step 3: Determine the sun's position;

[0010] Step 4: Determine the date;

[0011] Step 5: Determine the number of orbits to complete the calculation.

[0012] Furthermore, in step 1, the equatorial inertial coordinate system O is defined i X i Y i Z i :Coordinate origin O i Located at the center of the Earth, X i The axis is along the intersection of the Earth's equatorial plane and the ecliptic plane, pointing to the vernal equinox Y; Z i The Y axis points to the North Pole; i The axis is perpendicular to the X axis in the equatorial plane. i axis.

[0013] Furthermore, in step 2, one 01 solar sensor is installed on each of the six faces of the satellite to ensure that at least one 01 solar sensor can see the sun in the sunlit area; the orbital period T is calculated by the time interval between two times when the satellite enters the sunlit area from the earth's shadow area; according to Kepler's third law, the relationship between the satellite orbital period and the semi-major axis a is:

[0014]

[0015] Thus, the semi-major axis a of the orbit can be obtained.

[0016] Furthermore, in step 2, the satellite orbit plane is subjected to the non-spherical gravitational perturbation of the Earth, and only the harmonic term J is considered. 2 The long-term perturbation of the ascending node is The average value over one orbit is:

[0017]

[0018] Where: n is the average orbital speed, gravitational constant μ = 3.986006×10 5 km 3 / s 2 , with harmonic perturbation coefficient J 2 =0.001082, R e =6378.145km, i is the orbital inclination, a is the orbital semi-major axis, and e is the orbital eccentricity.

[0019] Further, in step 2, for the sun-synchronous orbit, the precession angular velocity of the satellite orbital plane is equal to the angular velocity of the mean sun moving on the equator; the mean sun makes an annual apparent motion along the equator, and the length of the tropical year is 365.2422 mean solar days. Therefore, the angular velocity of the mean sun moving on the equator is:

[0020]

[0021] Combining equations (2) and (3), assuming that the orbital eccentricity is e = 0, that is, the orbit is considered to be circular, the change in the right ascension Ω of the ascending node of the sun-synchronous orbit in one day is:

[0022]

[0023] Substituting the semi-major axis a into equation (4) yields the orbital inclination i.

[0024] Furthermore, in step 3, one digital sun sensor is installed on each of the six faces of the satellite to ensure that at least one digital sun sensor can see the sun in the sunlit area;

[0025] The digital sun sensitivity measurement value includes vector X and vector Y; vector X is defined as the sun direction vector The angle between the projection of the XOZ plane in the digital Taimin coordinate system and the +Z axis; the vector Y is defined as In the digital Taimin coordinate system, the angle between the projection of the YOZ plane and the +Z axis is The relationship is:

[0026]

[0027] According to the installation position of the digital Taimin on the satellite The components of the sun direction vector in the satellite body coordinate system are:

[0028]

[0029] The satellite is equipped with a star sensor to measure the conversion quaternion from the equatorial inertial coordinate system to the star sensor coordinate system. According to the installation position of the star sensor on the satellite Get the equatorial inertial coordinate system to the satellite body coordinate system O b X b Y b Z b The attitude conversion quaternion is:

[0030]

[0031] Combining equations (6) and (7), the components of the sun direction vector in the equatorial inertial coordinate system are:

[0032]

[0033] To keep the satellite in a stable inertial space, Is a constant value.

[0034] Furthermore, in step 4, the earth rotates around the sun from west to east, and the intersection angle between the ecliptic plane and the equatorial plane is α=23.43°; the geocentric ecliptic coordinate system O is defined s X s Y s Z s :Coordinate origin O s Located at the center of the Earth, X s The axis is along the intersection of the Earth's equatorial plane and the ecliptic plane, pointing to the vernal equinox Y; Z s The Y axis points to the north pole; s The axis is perpendicular to the X axis in the ecliptic plane. s axis;

[0035] The attitude conversion quaternion from the equatorial inertial coordinate system to the geocentric ecliptic coordinate system is:

[0036]

[0037] Combining equations (8) and (9), the components of the sun direction vector in the geocentric ecliptic coordinate system are:

[0038]

[0039] In X s O s Y s Rotation in the plane, with O s X s The angle between the axes is:

[0040]

[0041] θ∈(-180°,180°], θ=0° means it is the vernal equinox, and the current date can be calculated based on θ.

[0042] Further, in step 5, according to the orbital period T, the total number of circles the satellite orbits around the earth every day is N = 24 / T; from formula (4), it can be seen that the change in the right ascension of the ascending node of the sun-synchronous orbit in one day is ΔΩ = 0.9856°, so the change in the right ascension of the ascending node of the satellite per circle is:

[0043]

[0044] Combining equations (11) and (12), according to the sun direction vector With O s X s The increment dθ of the axis angle θ can be used to calculate the circle number of the satellite, which is the dθ / dΩ circle.

[0045] Beneficial effects of the present invention:

[0046] The present invention proposes a satellite astronomical navigation calculation method using an attitude sensor. The method of the present invention is designed to use a star sensor and a sun sensor for satellite autonomous navigation, and can effectively calculate the satellite orbit period T, semi-major axis a, inclination i and the current date, and the number of orbits. The method is simple in principle, easy to implement, and can be applied to engineering practice. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.

[0048] Figure 1 Schematic diagram of the satellite's spatial position.

[0049] Figure 2 Schematic diagram of the satellite's inertial space stable attitude.

[0050] Figure 3 Schematic diagram of the position relationship between the earth and the sun.

[0051] Figure 4 Schematic diagram of the on-orbit measurement accuracy of the digital Taimin.

[0052] Figure 5 This is a schematic diagram of digital Taimin vector data.

[0053] Figure 6 t 0 Schematic diagram of digital time-sensitive vector data.

[0054] Figure 7 Schematic diagram of the change in the angle of the sun direction vector.

[0055] Figure 8 Schematic diagram of satellite orbits. DETAILED DESCRIPTION

[0056] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0057] Combination Figure 1-Figure 8 The present invention proposes a satellite celestial navigation calculation method using an attitude sensor, the method comprising the following steps:

[0058] Step 1: Describe the spatial position of the satellite;

[0059] Step 2: Determine the satellite orbit period T, semi-major axis a and inclination i;

[0060] Step 3: Determine the sun's position;

[0061] Step 4: Determine the date;

[0062] Step 5: Determine the number of orbits to complete the calculation.

[0063] In step 1, define the equatorial inertial coordinate system O i X i Y i Z i :Coordinate origin O i Located at the center of the Earth, X i The axis is along the intersection of the Earth's equatorial plane and the ecliptic plane, pointing to the vernal equinox Y; Z i The Y axis points to the North Pole; i The axis is perpendicular to the X axis in the equatorial plane. i Axis, the position of the satellite in the equatorial inertial coordinate system is described as Figure 1 As shown in the figure, ON is the node line of the satellite orbit, N is the ascending node; Ω is the right ascension of the ascending node, and the node line ON and X i The angle between the two axes; S is the position of the satellite; P is the perigee of the satellite orbit; and i is the orbital inclination.

[0064] In step 2, one 01 solar sensor is installed on each of the six faces of the satellite to ensure that at least one 01 solar sensor can see the sun in the sunlit area. The orbital period T is calculated by the time interval between two times when the satellite enters the sunlit area from the earth's shadow area. According to Kepler's third law, the relationship between the satellite orbital period and the semi-major axis a is:

[0065]

[0066] Thus, the semi-major axis a of the orbit can be obtained.

[0067] In step 2, the satellite orbit plane is subjected to the non-spherical gravitational perturbation of the Earth, and only the harmonic term J is considered. 2 The long-term perturbation of the ascending node is The average value over one orbit is:

[0068]

[0069] Where: n is the average orbital speed, gravitational constant μ = 3.986006×10 5 km 3 / s 2 , with harmonic perturbation coefficient J 2 =0.001082, Re =6378.145km, i is the orbital inclination, a is the orbital semi-major axis, and e is the orbital eccentricity.

[0070] In step 2, for a sun-synchronous orbit, the precession angular velocity of the satellite orbital plane is equal to the angular velocity of the mean sun moving on the equator; the mean sun makes an annual apparent motion along the equator, and the length of the tropical year is 365.2422 mean solar days. Therefore, the angular velocity of the mean sun moving on the equator is:

[0071]

[0072] Combining equations (2) and (3), assuming that the orbital eccentricity is e = 0, that is, the orbit is considered to be circular, the change in the right ascension Ω of the ascending node of the sun-synchronous orbit in one day is:

[0073]

[0074] Substituting the semi-major axis a into equation (4) yields the orbital inclination i.

[0075] In step 3, one digital sun sensor is installed on each of the six faces of the satellite to ensure that at least one digital sun sensor can see the sun in the sunlit area;

[0076] The digital sun sensitivity measurement value includes vector X and vector Y; vector X is defined as the sun direction vector The angle between the projection of the XOZ plane in the digital Taimin coordinate system and the +Z axis; the vector Y is defined as In the digital Taimin coordinate system, the angle between the projection of the YOZ plane and the +Z axis is The relationship is:

[0077]

[0078] According to the installation position of the digital Taimin on the satellite The components of the sun direction vector in the satellite body coordinate system are:

[0079]

[0080] The satellite is equipped with a star sensor to measure the conversion quaternion from the equatorial inertial coordinate system to the star sensor coordinate system. According to the installation position of the star sensor on the satellite Get the equatorial inertial coordinate system to the satellite body coordinate system O b X b Y b Z b The attitude conversion quaternion is:

[0081]

[0082] Combining equations (6) and (7), the components of the sun direction vector in the equatorial inertial coordinate system are:

[0083]

[0084] like Figure 2 As shown, if the satellite maintains a stable attitude in inertial space, then Is a constant value.

[0085] In step 4, if Figure 3 As shown, the earth rotates around the sun from west to east, and the intersection angle between the ecliptic plane and the equatorial plane is α = 23.43°; the geocentric ecliptic coordinate system O is defined s X s Y s Z s :Coordinate origin O s Located at the center of the Earth, X s The axis is along the intersection of the Earth's equatorial plane and the ecliptic plane, pointing to the vernal equinox Y; Z s The Y axis points to the north pole; s The axis is perpendicular to the X axis in the ecliptic plane. s axis;

[0086] The attitude conversion quaternion from the equatorial inertial coordinate system to the geocentric ecliptic coordinate system is:

[0087]

[0088] Combining equations (8) and (9), the components of the sun direction vector in the geocentric ecliptic coordinate system are:

[0089]

[0090] In X s O s Y s Rotation in the plane, with O s X s The angle between the axes is:

[0091]

[0092] θ∈(-180°,180°], θ=0° means it is the vernal equinox, and the current date can be calculated based on θ.

[0093] In step 5, according to the orbital period T (h), the total number of circles the satellite orbits around the earth every day is N = 24 / T; from formula (4), it can be seen that the change in the right ascension of the ascending node of the sun-synchronous orbit in one day is ΔΩ = 0.9856°, so the change in the right ascension of the ascending node of the satellite per circle is:

[0094]

[0095] Combining equations (11) and (12), according to the sun direction vector With O s X s The increment dθ of the axis angle θ can be used to calculate the circle number of the satellite, which is the dθ / dΩ circle.

[0096] For example, a low-orbit sun-synchronous orbit satellite orbits the earth 15 times a day, and the right ascension of the ascending node changes by about 0.0657° per circle. Figure 4 As shown, from Figure 4 It can be seen that the high-frequency noise is about 0.05°. After filtering, the accuracy is further improved, and the change in the right ascension of the ascending node of the satellite per orbit can be identified.

[0097] Example

[0098] The first time the simulated satellite enters the sunlit area from the earth's shadow is 4:16:48 on September 16, 2023, and the second time it enters the sunlit area from the earth's shadow is 5:52:08 on September 16, 2023. The time difference between the two is the orbital period T = 5720s. Substituting it into formula (1), the semi-major axis of the orbit is:

[0099] a=(μ(T / 2π) 2 ) 1 / 3 =6912.84km (13)

[0100] Substituting the semi-major axis a into equation (4) yields the orbital inclination i

[0101]

[0102] Installation position of the star body + X-plane digital sensor on the star The measurement range of vector X and vector Y is [-60°, 60°], and the data is valid from September 16 to October 31, 2023. Figure 5 shown.

[0103] Starting time t 0 It is 4:17 on September 16, 2023 (UTC time: 748109820), at which time the digital sensitive vector X = 19.01°, vector Y = -13.9°, such as Figure 6 shown.

[0104] Installation location of star sensor on the satellite The measured value is Quaternion for attitude conversion from equatorial inertial coordinate system to satellite body coordinate system Quaternion for attitude conversion from equatorial inertial coordinate system to geocentric ecliptic coordinate system Then the attitude conversion quaternion from the digital Taimin coordinate system to the geocentric ecliptic coordinate system is:

[0105]

[0106] The sun direction vector is in the geocentric ecliptic coordinate system and O s X s The change of the axis angle θ is as follows Figure 7 As shown, when θ = 180°, the sun is at -X s The axis is the autumnal equinox (September 23).

[0107] like Figure 8 As shown in the figure, the angle θ is selected based on the time when the satellite enters the sunlit area from the shadow area. 0 =163.3376° is the 0th orbit of the satellite. At this time, the angle change from the autumnal equinox is 17°. The date is 17 / 0.9856≈17 days after the autumnal equinox, that is, October 10, 2023, which is consistent with the simulation (UTC time: 750217800). 1 =163.2724°, increment Δθ=θ 1 -θ 0 =-0.0652°; the second circle angle θ 2 =163.2066°, increment Δθ=θ 2 -θ 1 =-0.0658°, which is consistent with the theoretical analysis.

[0108] The above is a detailed introduction to a satellite celestial navigation calculation method using an attitude sensor proposed in the present invention. Specific examples are used in this article to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core idea. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation method and application scope. In summary, the content of this specification should not be understood as a limitation on the present invention.

Claims

1. A satellite celestial navigation calculation method using an attitude sensor, characterized in that: The method comprises the following steps: Step 1: Describe the spatial position of the satellite; Step 2: Determine the satellite orbit period T, semi-major axis a and inclination i; Step 3: Determine the sun's position; Step 4: Determine the date; Step 5: Determine the number of orbital circles to complete the calculation; In step 2, the satellite orbital plane is subjected to the non-spherical gravitational perturbation of the Earth, and only the long-term perturbation with the harmonic term J2 is considered. Then the rate of change of the right ascension of the ascending node is The average value over one orbit is: Where: n is the average orbital speed, gravitational constant μ = 3.986006×10 5 km 3 / s 2 , with harmonic perturbation coefficient J2=0.001082, R e =6378.145km, i is the orbital inclination, a is the orbital semi-major axis, and e is the orbital eccentricity; In step 2, for a sun-synchronous orbit, the precession angular velocity of the satellite orbital plane is equal to the angular velocity of the mean sun moving on the equator; the mean sun makes an annual apparent motion along the equator, and the length of the tropical year is 365.2422 mean solar days. Therefore, the angular velocity of the mean sun moving on the equator is: Combining equations (2) and (3), assuming that the orbital eccentricity is e = 0, that is, the orbit is considered to be circular, the change in the right ascension Ω of the ascending node of the sun-synchronous orbit in one day is: Substituting the semi-major axis a into equation (4) yields the orbital inclination i; In step 3, one digital sun sensor is installed on each of the six faces of the satellite to ensure that at least one digital sun sensor can see the sun in the sunlit area; The digital sun sensitivity measurement value includes vector X and vector Y; vector X is defined as the sun direction vector The angle between the projection of the XOZ plane in the digital Taimin coordinate system and the +Z axis; the vector Y is defined as In the digital Taimin coordinate system, the angle between the projection of the YOZ plane and the +Z axis is The relationship is: According to the installation position of the digital Taimin on the satellite The components of the sun direction vector in the satellite body coordinate system are: The satellite is equipped with a star sensor to measure the conversion quaternion from the equatorial inertial coordinate system to the star sensor coordinate system. According to the installation position of the star sensor on the satellite Get the equatorial inertial coordinate system to the satellite body coordinate system O b X b Y b Z b The attitude conversion quaternion is: Combining equations (6) and (7), the components of the sun direction vector in the equatorial inertial coordinate system are: To keep the satellite in a stable inertial space, Is a constant value.

2. The method according to claim 1, characterized in that In step 1, define the equatorial inertial coordinate system O i X i Y i Z i :Coordinate origin O i Located at the center of the Earth, X i The axis is along the intersection of the Earth's equatorial plane and the ecliptic plane, pointing to the vernal equinox Y; Z i The Y axis points to the North Pole; i The axis is perpendicular to the X axis in the equatorial plane. i axis.

3. The method according to claim 2, characterized in that In step 2, one 01 solar sensor is installed on each of the six faces of the satellite to ensure that at least one 01 solar sensor can see the sun in the sunlit area. The orbital period T is calculated by the time interval between two times when the satellite enters the sunlit area from the earth's shadow area. According to Kepler's third law, the relationship between the satellite orbital period and the semi-major axis a is: Thus, the semi-major axis a of the orbit can be obtained.

4. The method according to claim 1, characterized in that: In step 4, the earth rotates around the sun from west to east, and the intersection angle between the ecliptic plane and the equatorial plane is α = 23.43°; define the geocentric ecliptic coordinate system O s X s Y s Z s :Coordinate origin O s Located at the center of the Earth, X s The axis is along the intersection of the Earth's equatorial plane and the ecliptic plane, pointing to the vernal equinox Y; Z s The Y axis points to the north pole; s The axis is perpendicular to the X axis in the ecliptic plane. s axis; The attitude conversion quaternion from the equatorial inertial coordinate system to the geocentric ecliptic coordinate system is: Combining equations (8) and (9), the components of the sun direction vector in the geocentric ecliptic coordinate system are: In X s O s Y s Rotation in the plane, with O s X s The angle between the axes is: θ∈(-180°,180°], θ=0° means it is the vernal equinox, and the current date can be calculated based on θ.

5. The method according to claim 4, characterized in that In step 5, according to the orbital period T, the total number of circles the satellite orbits around the earth every day is N = 24 / T; from formula (4), it can be seen that the change in the right ascension of the ascending node of the sun-synchronous orbit in one day is ΔΩ = 0.9856°, so the change in the right ascension of the ascending node of the satellite per circle is: Combining equations (11) and (12), according to the sun direction vector With O s X s The increment dθ of the axis angle θ can be used to calculate the circle number of the satellite, which is the dθ / dΩ circle.

Images