Satellite attitude guidance angle calculation method and system based on GNSS absolute positioning data

Abstract

The present invention provides a method and system for calculating satellite attitude and guidance angles based on GNSS absolute positioning data, including: step S1: converting the time, position, and velocity in the GNSS real-time absolute positioning data into an inertial coordinate system to obtain the position and velocity in the inertial coordinate system; step S2: recursively integrating the position and velocity at the time of the GNSS real-time absolute positioning to the current satellite time to obtain the position and velocity at the current satellite time; step S3: converting the position and velocity information at the current time to obtain the orbital parameters at the corresponding time; and step S4: calculating the yaw and pitch guidance angles at the current time based on the orbital parameters at the current satellite time. This method can calculate satellite attitude and guidance angles in real time with high precision, effectively supporting high-precision attitude control, effectively compensating for image deformation caused by the Earth's rotation, and improving imaging quality.

Description

Technical Field

[0001] The present invention relates to the field of satellite attitude guidance planning and calculation, and in particular to a method and system for calculating satellite attitude guidance angles based on GNSS absolute positioning data. Background Art

[0002] Satellites in orbit are affected by the Earth's rotation and elliptical orbits, resulting in image displacement and distortion. For SAR satellites in particular, the Doppler center can shift due to the Earth's curvature and rotation. Typical offsets exceed the system's pulse repetition frequency, affecting image quality. Therefore, two-dimensional attitude guidance angle calculation is required to minimize these effects, improve image quality, and simplify post-processing.

[0003] Huai Chao, Wang Wenyan. Research on zero-Doppler attitude guidance of InSAR formation satellites[J]. Shanghai Aerospace, 2014, 31(6). A zero-Doppler attitude guidance method for InSAR formation satellites was derived. The analytical expressions of yaw angle and pitch angle were given based on the Doppler center frequency of the auxiliary satellite.

[0004] Guo Deming, Xu Huaping, Li Jingwen. Attitude Guidance of High-Resolution Spaceborne Squint SAR[J]. Acta Astronautics, 2011(5). Based on the spatial geometric relationship between the satellite and the ground and the attitude guidance principle, the analytical expressions of the yaw angle and pitch angle for achieving full zero Doppler guidance are derived, and the attitude guidance law that minimizes the variation of the Doppler center frequency with distance is further given and simulated.

[0005] However, none of the above literatures involve the processing and introduction of GNSS absolute positioning data in engineering practice.

[0006] Patent CN106843249A, "A Two-Dimensional Guidance Attitude Control Method," describes a high-precision, high-stability two-dimensional guidance attitude control method. This method enables rapid two-dimensional guidance access control at any satellite position and improves system control accuracy. However, it does not involve guidance angle calculation.

[0007] Patent document CN106915477B (application number: CN201710128823.6) discloses an attitude control method that includes the following steps: using pulse-per-second signals to align onboard time, which affects satellite attitude accuracy; correcting star sensor errors used for satellite attitude measurement in real time; using dynamic coupling to compensate for the gravity gradient interference torque of obliquely flying satellites; using input shaping control to suppress satellite flexibility during attitude control; implementing high-precision and high-stability attitude guidance control using an attitude control law that adds angular velocity feedforward commands to a position-velocity dual loop and a position correction loop; and implementing a saturated sliding mode variable structure control algorithm to achieve rapid attitude maneuvers for high-inertia satellites. However, this invention does not address the processing and integration of GNSS absolute positioning data in practical engineering applications. Summary of the Invention

[0008] In view of the defects in the prior art, the purpose of the present invention is to provide a method and system for calculating satellite attitude guidance angles based on GNSS absolute positioning data.

[0009] According to the present invention, a method for calculating satellite attitude guidance angle based on GNSS absolute positioning data includes:

[0010] Step S1: converting the time, position, and velocity in the GNSS real-time absolute positioning data into an inertial coordinate system to obtain the position and velocity in the inertial coordinate system;

[0011] Step S2: recursively integrate the position and velocity at the time of GNSS real-time absolute positioning to the current satellite time to obtain the position and velocity at the current satellite time;

[0012] Step S3: Convert the current position and velocity information to obtain the orbit parameters at the corresponding moment;

[0013] Step S4: Calculate the yaw steering angle and pitch steering angle at the current moment according to the orbit parameters at the current satellite time.

[0014] Preferably, in step S1:

[0015] Step S1.1:

[0016] Calculate the coordinate transformation matrix from the J2000 Earth-centered inertial system to the Earth-fixed coordinate system based on time T0:

[0017] M ECI2ECF =EP·ER·NR·PR

[0018] Among them, EP is the polar motion matrix, which is simplified to a 3×3 unit matrix; ER is the Earth rotation matrix, NR is the nutation matrix, and PR is the precession matrix;

[0019] Step S1.2:

[0020] The position and velocity information of the GNSS absolute positioning data is converted from the WGS84 coordinate system to the J2000.0 inertial coordinate system to obtain the position r ECI0 , speed v ECI0 for:

[0021] r ECI0 =(M ECI2ECF ) -1 ·r ECF0

[0022] v ECI0 =(M ECI2ECF ) -1 v ECF0 +ω e ×r ECF0

[0023] Among them, r ECFO is the position in the WGS84 coordinate system at time T0, v ECFO is the velocity in the WGS84 coordinate system at time T0; The T in the above table represents transposition, and the · in the above table represents the derivative with respect to time.

[0024] Preferably, in step S2:

[0025] According to the position velocity at time T0, the fourth-order Runge-Kutta formula is used to numerically integrate and recursively carry out the current star time T to obtain the position r at time T. ECI , speed v ECI .

[0026] Preferably, in step S3:

[0027] Step S3.1: Calculate process parameters:

[0028] Position scalar r = |r ECI |, velocity scalar v = |v ECI |

[0029] Conversion matrix from orbital coordinate system to J2000.0 inertial coordinate system:

[0030]

[0031] Where Y = v ECI ×r ECI , X=r ECI ×Y, are all 3×1 column vectors; |X| represents the modulus of X;

[0032] Orbital energy constant:

[0033]

[0034] Among them, the Earth's gravitational constant μ = 398600km 3 / s 2 ;

[0035] Moment of momentum:

[0036] H=r ECI ×v ECI

[0037] The modulus of the moment of momentum:

[0038] H=|H|

[0039] Semi-diameter:

[0040] p=H 2 / μ

[0041] Laplace vector:

[0042]

[0043] Among them, L x For L y For L z The components of the Laplace vector L on the three coordinate axes in the inertial coordinate system respectively;

[0044] Step S3.2: Calculate the number of orbital elements at time T:

[0045] Semi-major axis:

[0046]

[0047] Eccentricity:

[0048]

[0049] Orbital inclination:

[0050] i=arccos(H zi / H)

[0051] Among them, H zi is the component of the angular momentum H in the Z direction;

[0052] Latitude Argument:

[0053] u=arctan2(-M orb2ECI33 ,M orb2ECI31 )

[0054] Among them, M orb2ECI33 is the matrix M orb2ECI The value in row 3 and column 3, M orb2ECI31 is the matrix M orb2ECI The value in row 3 and column 1;

[0055] Ascending node right ascension:

[0056] Ω=arctan2(-M orb2ECI12 ,M orb2ECI22 )

[0057] Among them, M orb2ECI12 is the matrix M orb2ECI The value in row 1 and column 2, M orb2ECI22 is the matrix M orb2ECI The value in row 2 and column 2;

[0058] Argument of perigee:

[0059]

[0060] True anomaly:

[0061] f=u-ω。

[0062] Preferably, in step S4:

[0063] The yaw steering angle ψ is:

[0064]

[0065] Where: e is the Earth's rotation angular velocity, a is the semi-major axis, r is the position scalar, u is the latitude argument, i is the orbital inclination, e is the eccentricity, and f is the true anomaly;

[0066]

[0067] μ=398600km 3 / s 2 ;

[0068] The pitch steering angle θ is:

[0069]

[0070] According to the present invention, a satellite attitude and guidance angle calculation system based on GNSS absolute positioning data is provided, comprising:

[0071] Module M1: Converts the time, position, and speed in the GNSS real-time absolute positioning data into the inertial coordinate system to obtain the position and speed in the inertial coordinate system;

[0072] Module M2: recursively integrate the position and velocity at the time of GNSS real-time absolute positioning to the current satellite time to obtain the position and velocity at the current satellite time;

[0073] Module M3: Convert the current position and velocity information to obtain the orbit parameters at the corresponding moment;

[0074] Module M4: Calculates the yaw and pitch steering angles at the current moment based on the orbital parameters of the current satellite time.

[0075] Preferably, in the module M1:

[0076] Module M1.1:

[0077] Calculate the coordinate transformation matrix from the J2000 Earth-centered inertial system to the Earth-fixed coordinate system based on time T0:

[0078] M ECI2ECF =EP·ER·NR·PR

[0079] Among them, EP is the polar motion matrix, which is simplified to a 3×3 unit matrix; ER is the Earth rotation matrix, NR is the nutation matrix, and PR is the precession matrix;

[0080] Module M1.2:

[0081] The position and velocity information of the GNSS absolute positioning data is converted from the WGS84 coordinate system to the J2000.0 inertial coordinate system to obtain the position r ECI0 , speed v ECI0 for:

[0082] r ECI0 =(M ECI2ECF ) -1 ·r ECF0

[0083] v ECI0 =(M ECI2ECF ) -1 v ECF0 +ω e ×r ECF0

[0084] Among them, r ECFO is the position in the WGS84 coordinate system at time T0, v ECFO is the velocity in the WGS84 coordinate system at time T0; The T in the above table represents transposition, and the · in the above table represents the derivative with respect to time.

[0085] Preferably, in the module M2:

[0086] According to the position velocity at time T0, the fourth-order Runge-Kutta formula is used to numerically integrate and recursively carry out the current star time T to obtain the position r at time T. ECI , speed v ECI .

[0087] Preferably, in the module M3:

[0088] Module M3.1: Calculation process parameters:

[0089] Position scalar r = |r ECI |, velocity scalar v = |v ECI |

[0090] Conversion matrix from orbital coordinate system to J2000.0 inertial coordinate system:

[0091]

[0092] Where Y = v ECI ×r ECI , X=r ECI ×Y, are all 3×1 column vectors; |X| represents the modulus of X;

[0093] Orbital energy constant:

[0094]

[0095] Among them, the Earth's gravitational constant μ = 398600km 3 / s 2 ;

[0096] Moment of momentum:

[0097] H=r ECI ×v ECI

[0098] The modulus of the moment of momentum:

[0099] H=|H|

[0100] Semi-diameter:

[0101] p=H 2 / μ

[0102] Laplace vector:

[0103]

[0104] Among them, L x For L y For L z The components of the Laplace vector L on the three coordinate axes in the inertial coordinate system respectively;

[0105] Module M3.2: Calculate the number of orbital elements at time T:

[0106] Semi-major axis:

[0107]

[0108] Eccentricity:

[0109]

[0110] Orbital inclination:

[0111] i=arccos(H zi / H)

[0112] Among them, H zi is the component of the angular momentum H in the Z direction;

[0113] Latitude Argument:

[0114] u=arctan2(-M orb2ECI33 ,M orb2ECI31 )

[0115] Among them, M orb2ECI33 is the matrix M orb2ECI The value in row 3 and column 3, M orb2ECI31 is the matrix M orb2ECI The value in row 3 and column 1;

[0116] Ascending node right ascension:

[0117] Ω=arctan2(-M orb2ECI12 ,M orb2ECI22 )

[0118] Among them, M orb2ECI12 is the matrix M orb2ECI The value in row 1 and column 2, M orb2ECI22 is the matrix M orb2ECI The value in row 2 and column 2;

[0119] Argument of perigee:

[0120]

[0121] True anomaly:

[0122] f=u-ω。

[0123] Preferably, in the module M4:

[0124] The yaw steering angle ψ is:

[0125]

[0126] Where: e is the Earth's rotation angular velocity, a is the semi-major axis, r is the position scalar, u is the latitude argument, i is the orbital inclination, e is the eccentricity, and f is the true anomaly;

[0127]

[0128] μ=398600km 3 / s 2 ;

[0129] The pitch steering angle θ is:

[0130]

[0131] Compared with the prior art, the present invention has the following beneficial effects:

[0132] 1. This calculation method is highly effective;

[0133] 2. The present invention combines engineering practice to calculate satellite attitude guidance angles in real time with high precision;

[0134] 3. This method can calculate the satellite attitude guidance angle in real time with high precision, effectively support high-precision attitude control, effectively compensate for image deformation caused by the rotation of the earth, and improve imaging quality. BRIEF DESCRIPTION OF THE DRAWINGS

[0135] Other features, objects and advantages of the present invention will become more apparent upon reading the detailed description of non-limiting embodiments with reference to the following drawings:

[0136] Figure 1 This is a flow chart of the satellite attitude guidance angle calculation based on GNSS absolute positioning data in the present invention;

[0137] Figure 2 The attitude steering angle of a certain type of satellite is calculated according to the method of the present invention, wherein the rolling direction steering angle is 0. DETAILED DESCRIPTION

[0138] The present invention will be described in detail below with reference to specific embodiments. The following examples will help those skilled in the art to further understand the present invention, but are not intended to limit the present invention in any form. It should be noted that, for those skilled in the art, several changes and improvements can be made without departing from the scope of the present invention. These all fall within the scope of protection of the present invention.

[0139] Example 1:

[0140] Satellites in orbit rely on the GNSS system to provide high-precision absolute positioning data that is updated at a certain frequency in real time. The present invention uses this data to calculate the current attitude guidance angle in real time with high precision, effectively supporting high-precision attitude control and thus improving imaging quality.

[0141] According to the present invention, a satellite attitude guidance angle calculation method based on GNSS absolute positioning data is provided. Figure 1-Figure 2 Shown, including:

[0142] Step S1: converting the time, position, and velocity in the GNSS real-time absolute positioning data into an inertial coordinate system to obtain the position and velocity in the inertial coordinate system;

[0143] Specifically, in step S1:

[0144] Step S1.1:

[0145] Calculate the coordinate transformation matrix from the J2000 Earth-centered inertial system to the Earth-fixed coordinate system based on time T0:

[0146] M ECI2ECF =EP·ER·NR·PR

[0147] Among them, EP is the polar motion matrix, which is simplified to a 3×3 unit matrix; ER is the Earth rotation matrix, NR is the nutation matrix, and PR is the precession matrix;

[0148] Step S1.2:

[0149] The position and velocity information of the GNSS absolute positioning data is converted from the WGS84 coordinate system to the J2000.0 inertial coordinate system to obtain the position r ECI0 , speed v ECI0 for:

[0150] r ECI0 =(M ECI2ECF ) -1 ·r ECF0

[0151] v ECI0 =(M ECI2ECF ) -1 v ECF0 +ω e ×r ECF0

[0152] Among them, r ECFO is the position in the WGS84 coordinate system at time T0, v ECFO is the velocity in the WGS84 coordinate system at time T0; The T in the above table represents transposition, and the · in the above table represents the derivative with respect to time.

[0153] Step S2: recursively integrate the position and velocity at the time of GNSS real-time absolute positioning to the current satellite time to obtain the position and velocity at the current satellite time;

[0154] Specifically, in step S2:

[0155] According to the position velocity at time T0, the fourth-order Runge-Kutta formula is used to numerically integrate and recursively carry out the current star time T to obtain the position r at time T. ECI , speed v ECI .

[0156] Step S3: Convert the current position and velocity information to obtain the orbit parameters at the corresponding moment;

[0157] Specifically, in step S3:

[0158] Step S3.1: Calculate process parameters:

[0159] Position scalar r = |r ECI |, velocity scalar v = |v ECI |

[0160] Conversion matrix from orbital coordinate system to J2000.0 inertial coordinate system:

[0161]

[0162] Where Y = v ECI ×r ECI , X=r ECI ×Y, are all 3×1 column vectors; |X| represents the modulus of X;

[0163] Orbital energy constant:

[0164]

[0165] Among them, the Earth's gravitational constant μ = 398600km 3 / s 2 ;

[0166] Moment of momentum:

[0167] H=r ECI ×v ECI

[0168] The modulus of the moment of momentum:

[0169] H=|H|

[0170] Semi-diameter:

[0171] p=H 2 / μ

[0172] Laplace vector:

[0173]

[0174] Among them, L x For L y For L z The components of the Laplace vector L on the three coordinate axes in the inertial coordinate system respectively;

[0175] Step S3.2: Calculate the number of orbital elements at time T:

[0176] Semi-major axis:

[0177]

[0178] Eccentricity:

[0179]

[0180] Orbital inclination:

[0181] i=arccos(H zi / H)

[0182] Among them, H zi is the component of the angular momentum H in the Z direction;

[0183] Latitude Argument:

[0184] u=arctan2(-M orb2ECI33 ,M orb2ECI31 )

[0185] Among them, M orb2ECI33 is the matrix M orb2ECI The value in row 3 and column 3, M orb2ECI31 is the matrix M orb2ECI The value in row 3 and column 1;

[0186] Ascending node right ascension:

[0187] Ω=arctan2(-M orb2ECI12 ,M orb2ECI22 )

[0188] Among them, M orb2ECI12 is the matrix M orb2ECI The value in row 1 and column 2, M orb2ECI22 is the matrix M orb2ECI The value in row 2 and column 2;

[0189] Argument of perigee:

[0190]

[0191] True anomaly:

[0192] f=u-ω。

[0193] Step S4: Calculate the yaw steering angle and pitch steering angle at the current moment according to the orbit parameters at the current satellite time.

[0194] Specifically, in step S4:

[0195] The yaw steering angle ψ is:

[0196]

[0197] Where: e is the Earth's rotation angular velocity, a is the semi-major axis, r is the position scalar, u is the latitude argument, i is the orbital inclination, e is the eccentricity, and f is the true anomaly;

[0198]

[0199] μ=398600km 3 / s 2 ;

[0200] The pitch steering angle θ is:

[0201]

[0202] Example 2:

[0203] Example 2 is a preferred example of Example 1 and is used to illustrate the present invention in more detail.

[0204] Those skilled in the art may understand the satellite attitude steering angle calculation method based on GNSS absolute positioning data provided by the present invention as a specific implementation of the satellite attitude steering angle calculation system based on GNSS absolute positioning data, that is, the satellite attitude steering angle calculation system based on GNSS absolute positioning data can be implemented by executing the step flow of the satellite attitude steering angle calculation method based on GNSS absolute positioning data.

[0205] According to the present invention, a satellite attitude and guidance angle calculation system based on GNSS absolute positioning data is provided, comprising:

[0206] Module M1: Converts the time, position, and speed in the GNSS real-time absolute positioning data into the inertial coordinate system to obtain the position and speed in the inertial coordinate system;

[0207] Specifically, in the module M1:

[0208] Module M1.1:

[0209] Calculate the coordinate transformation matrix from the J2000 Earth-centered inertial system to the Earth-fixed coordinate system based on time T0:

[0210] M ECI2ECF =EP·ER·AR·PR

[0211] Among them, EP is the polar motion matrix, which is simplified to a 3×3 unit matrix; ER is the Earth rotation matrix, NR is the nutation matrix, and PR is the precession matrix;

[0212] Module M1.2:

[0213] The position and velocity information of the GNSS absolute positioning data is converted from the WGS84 coordinate system to the J2000.0 inertial coordinate system to obtain the position r ECI0 , speed v ECI0 for:

[0214] r ECI0 =(MECI2ECF ) -1 ·r ECF0

[0215] v ECI0 =(M ECI2ECF ) -1 v ECF0 +ω e ×r ECF0

[0216] Among them, r ECFO is the position in the WGS84 coordinate system at time T0, v ECFO is the velocity in the WGS84 coordinate system at time T0; The T in the above table represents transposition, and the · in the above table represents the derivative with respect to time.

[0217] Module M2: recursively integrate the position and velocity at the time of GNSS real-time absolute positioning to the current satellite time to obtain the position and velocity at the current satellite time;

[0218] Specifically, in the module M2:

[0219] According to the position velocity at time T0, the fourth-order Runge-Kutta formula is used to numerically integrate and recursively carry out the current star time T to obtain the position r at time T. ECI , speed v ECI .

[0220] Module M3: Convert the current position and velocity information to obtain the orbit parameters at the corresponding moment;

[0221] Specifically, in the module M3:

[0222] Module M3.1: Calculation process parameters:

[0223] Position scalar r = |r ECI |, velocity scalar v = |v ECI |

[0224] Conversion matrix from orbital coordinate system to J2000.0 inertial coordinate system:

[0225]

[0226] Where Y = v ECI ×r ECI , X=r ECI ×Y, are all 3×1 column vectors; |X| represents the modulus of X;

[0227] Orbital energy constant:

[0228]

[0229] Among them, the Earth's gravitational constant μ = 398600km 3 / s 2 ;

[0230] Moment of momentum:

[0231] H=r ECI ×v ECI

[0232] The modulus of the moment of momentum:

[0233] H=|H|

[0234] Semi-diameter:

[0235] p=H 2 / μ

[0236] Laplace vector:

[0237]

[0238] Among them, L x For L y For L z The components of the Laplace vector L on the three coordinate axes in the inertial coordinate system respectively;

[0239] Module M3.2: Calculate the number of orbital elements at time T:

[0240] Semi-major axis:

[0241]

[0242] Eccentricity:

[0243]

[0244] Orbital inclination:

[0245] i=arccos(H zi / H)

[0246] Among them, H zi is the component of the angular momentum H in the Z direction;

[0247] Latitude Argument:

[0248] u=arctan2(-M orb2ECI33 ,M orb2ECI31 )

[0249] Among them, M orb2ECI33 is the matrix M orb2ECI The value in row 3 and column 3, M orb2ECI31 is the matrix M orb2ECI The value in row 3 and column 1;

[0250] Ascending node right ascension:

[0251] Ω=arctan2(-M orb2ECI12 ,M orb2ECI22 )

[0252] Among them, M orb2ECI12 is the matrix M orb2ECI The value in row 1 and column 2, M orb2ECI22 is the matrix M orb2ECI The value in row 2 and column 2;

[0253] Argument of perigee:

[0254]

[0255] True anomaly:

[0256] f=u-ω。

[0257] Module M4: Calculates the yaw and pitch steering angles at the current moment based on the orbital parameters of the current satellite time.

[0258] Specifically, in the module M4:

[0259] The yaw steering angle ψ is:

[0260]

[0261] Where: e is the Earth's rotation angular velocity, a is the semi-major axis, r is the position scalar, u is the latitude argument, i is the orbital inclination, e is the eccentricity, and f is the true anomaly;

[0262]

[0263] μ=398600km 3 / s 2 ;

[0264] The pitch steering angle θ is:

[0265]

[0266] Example 3:

[0267] Example 3 is a preferred example of Example 1 and is used to illustrate the present invention in more detail.

[0268] To address the need for real-time calculation of two-dimensional guidance angles for satellites in orbit to improve the quality of satellite payload imaging, the present invention discloses a method for calculating satellite attitude guidance angles based on GNSS absolute positioning data. By converting and processing the absolute positioning data output by the GNSS in real time, the satellite's current orbital elements are obtained, and the satellite's yaw and pitch guidance angles are calculated. This aspect, combined with practical engineering practices, enables high-precision, real-time calculation of satellite attitude guidance angles. To achieve this objective, the present invention is implemented through the following technical solutions, specifically comprising the following steps.

[0269] A method for calculating satellite attitude guidance angle based on GNSS absolute positioning data comprises the following steps:

[0270] Step 1: According to the time T0 and position r in the GNSS real-time absolute positioning data ECF0 , speed v ECF0 Convert to the inertial coordinate system to describe the position r ECI0 , speed v ECI0 .

[0271] Step 2: Numerically integrate the position and velocity at time T0 to the current star time T to obtain the position r at time T ECI , speed v ECI .

[0272] Step 3: Convert the current position and velocity information to obtain the orbital parameters at the corresponding moment.

[0273] Step 4: Based on the orbit parameters at time T, calculate the yaw steering angle and pitch steering angle at the current moment respectively.

[0274] The step 1 includes the following:

[0275] Step 1-1: Calculate the coordinate transformation matrix M from the J2000 Earth-centered inertial system to the Earth-fixed coordinate system based on time T0 ECI2ECF =EP·ER·NR·PR, EP is the polar motion matrix, simplified to a 3×3 unit matrix; ER is the Earth rotation matrix, NR is the nutation matrix, and PR is the precession matrix.

[0276] Step 1-2: Convert the GNSS absolute positioning data position and velocity information from the WGS84 coordinate system to the J2000.0 inertial coordinate system to obtain the position r ECI0 , speed v ECI0 for

[0277] r ECI0 =(M ECI2ECF ) -1 ·r ECF0

[0278] v ECI0 =(M ECI2ECF )-1 v ECF0 +ω e ×r ECF0

[0279] in, The T in the above table represents transposition, and the · in the above table represents the derivative with respect to time.

[0280] According to the position velocity at time T0, the fourth-order Runge-Kutta formula is used to numerically integrate and recursively carry out the numerical integration to the current star time T, and the position r at time T is obtained. ECI , speed v ECI .

[0281] Step 3 includes the following:

[0282] Step 3-1: Calculate process parameters

[0283] Position scalar r = |r ECI |, velocity scalar v = |v ECI |

[0284] Conversion matrix from orbital coordinate system to J2000.0 inertial coordinate system Where Y = v ECI ×r ECI , X=r ECI ×Y, are both 3×1 column vectors. , |X| represents the modulus of X.

[0285] Orbital energy constant Moment of momentum H = r ECI ×v ECI , the modulus of the moment of momentum H = |H|, the semi-diameter p = H 2 / μ, Laplace vector Earth's gravitational constant μ = 398,600 km 3 / s 2 .

[0286] Step 3-2: Calculate the number of orbital elements at time T

[0287] semi-major axis

[0288] Eccentricity

[0289] Orbital inclination i = arccos(H zi / H), H zi is the component of the angular momentum H in the Z direction.

[0290] Latitude argument u=arctan2(-M orb2ECI33 ,M orb2ECI31 ).

[0291] Ascending node right ascension Ω=arctan2(-Morb2ECI12 ,M orb2ECI22 ).

[0292] Argument of perigee

[0293] True anomaly f=u-ω

[0294] The yaw steering angle ψ is:

[0295]

[0296] in: ω e is the Earth's rotational angular velocity

[0297] The pitch steering angle θ is:

[0298]

[0299] Those skilled in the art will appreciate that, in addition to implementing the system, device, and various modules provided by the present invention in purely computer-readable program code, it is entirely possible to implement the same program in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, embedded microcontrollers, and the like by logically programming the method steps. Therefore, the system, device, and various modules provided by the present invention can be considered a hardware component, and the modules included therein for implementing various programs can also be considered structures within the hardware component; the modules for implementing various functions can also be considered both software programs for implementing the method and structures within the hardware component.

[0300] The above describes specific embodiments of the present invention. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art may make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. The embodiments of this application and the features in the embodiments may be combined with each other in any manner unless there is a conflict.

Claims

1. A method for calculating satellite attitude and guidance angle based on GNSS absolute positioning data, characterized in that: include: Step S1: converting the time, position, and velocity in the GNSS real-time absolute positioning data into an inertial coordinate system to obtain the position and velocity in the inertial coordinate system; Step S2: recursively integrate the position and velocity at the time of GNSS real-time absolute positioning to the current satellite time to obtain the position and velocity at the current satellite time; Step S3: Convert the position and velocity information of the current star time to obtain the orbital parameters of the corresponding time; Step S4: Calculate the yaw steering angle and pitch steering angle at the current satellite time according to the orbit parameters at the current satellite time; In step S1: Step S1.1: Calculate the coordinate transformation matrix from the J2000 Earth-centered inertial system to the Earth-fixed coordinate system based on time T0: M ECI2ECF =EP ER NR PR Among them, EP is the polar motion matrix, which is simplified to a 3×3 unit matrix; ER is the Earth rotation matrix, NR is the nutation matrix, and PR is the precession matrix; Step S1.2: Convert the GNSS absolute positioning data position and velocity information from the WGS84 coordinate system to the J2000 inertial coordinate system to obtain the position r ECI0 , speed v ECI0 for: r ECI0 ＝(M ECI2ECF ) -1 ·r ECF0 v ECI0 =(M ECI2ECF ) -1 v ECF0 +ω e × r ECF0 Among them, r ECFO is the position in the WGS84 coordinate system at time T0, v ECFO is the velocity in the WGS84 coordinate system at time T0; The superscript T indicates transposition, and the superscript · indicates the time derivative; In step S2: According to the position velocity at time T0, the fourth-order Runge-Kutta formula is used to numerically integrate and recursively carry out the current star time T to obtain the position r at time T. ECI , speed v ECI ; In step S3: Step S3.1: Calculate process parameters: Position scalar r = |r ECI |, velocity scalar v = |v ECI | Conversion matrix from orbital coordinate system to J2000 inertial coordinate system: Where Y = v ECI ×r ECI , X=r ECI ×Y, are all 3×1 column vectors; |X| represents the modulus of X; orbital energy constant: Among them, the Earth's gravitational constant μ = 398600km 3 / s 2 ; Moment of momentum: H=r ECI ×v ECI The modulus of the moment of momentum: H=|H| Semi-diameter: p=H 2 / μ Laplace vector: Among them, L x 、L y 、L z are the components of the three coordinate axes of the Laplace vector L in the inertial coordinate system; Step S3.2: Calculate the number of orbital elements at time T: Semi-major axis: Eccentricity: Orbital inclination: i=arccos(H zi / H) Among them, H zi is the component of the angular momentum H in the Z direction; Latitude Argument: u=arctan2(-M orb2ECI33 ,M orb2ECI31 ) Among them, M orb2ECI33 is the matrix M orb2ECI The value in row 3 and column 3, M orb2ECI31 is the matrix M orb2ECI The value in row 3 and column 1; Ascending node right ascension: Ω=arctan2(-M orb2ECI12 ,M orb2ECI22 ) Among them, M orb2ECI12 is the matrix M orb2ECI The value in row 1 and column 2, M orb2ECI22 is the matrix M orb2ECI The value in row 2 and column 2; Argument of perigee: True anomaly: f=u-ω。 2. The satellite attitude guidance angle calculation method based on GNSS absolute positioning data according to claim 1, characterized in that: In step S4: The yaw steering angle ψ is: Where: e is the Earth's rotation angular velocity, a is the semi-major axis, r is the position scalar, u is the latitude argument, i is the orbital inclination, e is the eccentricity, and f is the true anomaly; μ=398600km 3 / s 2 ; The pitch steering angle θ is:

3. A satellite attitude and guidance angle calculation system based on GNSS absolute positioning data, characterized in that: include: Module M1: Converts the time, position, and speed in the GNSS real-time absolute positioning data into the inertial coordinate system to obtain the position and speed in the inertial coordinate system; Module M2: recursively integrate the position and velocity at the time of GNSS real-time absolute positioning to the current satellite time to obtain the position and velocity at the current satellite time; Module M3: Convert the position and velocity information of the current star time into the orbital parameters of the corresponding time; Module M4: Calculate the yaw and pitch angles of the current satellite time based on the orbital parameters of the current satellite time; In the module M1: Module M1.1: Calculate the coordinate transformation matrix from the J2000 Earth-centered inertial system to the Earth-fixed coordinate system based on time T0: M ECI2ECF =EP ER NR PR Among them, EP is the polar motion matrix, which is simplified to a 3×3 unit matrix; ER is the Earth rotation matrix, NR is the nutation matrix, and PR is the precession matrix; Module M1.2: Convert the GNSS absolute positioning data position and velocity information from the WGS84 coordinate system to the J2000 inertial coordinate system to obtain the position r ECI0 , speed v ECI0 for: r ECI0 ＝(M ECI2ECF ) -1 ·r ECF0 v ECI0 =(M ECI2ECF ) -1 v ECF0 +ω e ×r ECF0 Among them, r ECFO is the position in the WGS84 coordinate system at time T0, v ECFO is the velocity in the WGS84 coordinate system at time T0; The superscript T indicates transposition, and the superscript · indicates the time derivative; In the module M2: According to the position velocity at time T0, the fourth-order Runge-Kutta formula is used to numerically integrate and recursively carry out the current star time T to obtain the position r at time T. ECI , speed v ECI ; In the module M3: Module M3.1: Calculation process parameters: Position scalar r = |r ECI |, velocity scalar v = |v ECI | Conversion matrix from orbital coordinate system to J2000 inertial coordinate system: Where Y = v ECI ×r ECI , X=r ECI ×Y, are all 3×1 column vectors; |X| represents the modulus of X; Orbital energy constant: Among them, the Earth's gravitational constant μ = 398600km 3 / s 2 ; Moment of momentum: H=r ECI ×v ECI The modulus of the moment of momentum: H=|H| Semi-diameter: p=H 2 / μ Laplace vector: Among them, L x 、L y 、L z are the components of the three coordinate axes of the Laplace vector L in the inertial coordinate system; Module M3.2: Calculate the number of orbital elements at time T: Semi-major axis: Eccentricity: Orbital inclination: i=arccos(H zi / H) Among them, H zi is the component of the angular momentum H in the Z direction; Latitude Argument: u=arctan2(-M orb2ECI33 ,M orb2ECI31 ) Among them, M orb2ECI33 is the matrix M orb2ECI The value in row 3 and column 3, M orb2ECI31 is the matrix M orb2ECI The value in row 3 and column 1; Ascending node right ascension: Ω=arctan2(-M orb2ECI12 ,M orb2ECI22 ) Among them, M orb2ECI12 is the matrix M orb2ECI The value in row 1 and column 2, M orb2ECI22 is the matrix M orb2ECI The value in row 2 and column 2; Argument of perigee: True anomaly: f=u-ω。 4. The satellite attitude and guidance angle calculation system based on GNSS absolute positioning data according to claim 3, characterized in that: In the module M4: The yaw steering angle ψ is: Where: e is the Earth's rotation angular velocity, a is the semi-major axis, r is the position scalar, u is the latitude argument, i is the orbital inclination, e is the eccentricity, and f is the true anomaly; μ=398600km 3 / s 2 ; The pitch steering angle θ is:

Images