A satellite atmospheric perturbation measurement method and device

Abstract

The present invention discloses a satellite atmospheric perturbation measurement method and device. The method includes: obtaining the satellite coordinates, satellite speed and current time of the satellite based on Beidou positioning and timing functions; obtaining the satellite measurement attitude information and solar sailboard rotation angle information of the satellite; determining the total perturbation of the satellite in the inertial coordinate system according to the satellite coordinates, satellite speed and current time; determining the non-spherical gravitational perturbation corresponding to the satellite according to the satellite position and the third-party gravitational field model; determining the daily gravitational perturbation and lunar gravitational perturbation corresponding to the satellite according to the solar coordinates and lunar coordinates at the current time; determining the solar light pressure perturbation corresponding to the satellite according to the current time, satellite coordinates, solar coordinates, satellite measurement attitude information and solar sailboard rotation angle information; determining the satellite atmospheric perturbation of the satellite according to the total perturbation, non-spherical gravitational perturbation, daily gravitational perturbation, lunar gravitational perturbation and solar light pressure perturbation. The present invention can improve the calculation accuracy of the atmospheric perturbation.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite atmospheric perturbation measurement and calculation, and in particular to a satellite atmospheric perturbation measurement and calculation method and device. Background Art

[0002] The orbital perturbations of low-orbit satellites (orbital altitude 100-1500km) in space are mainly non-spherical gravitational perturbations and atmospheric perturbations, followed by solar and lunar gravitational perturbations and solar pressure perturbations. Orbital prediction requires accurate estimation of various perturbations, and non-spherical gravity, solar and lunar gravity, and solar pressure can all be solved well: for non-spherical gravity, it can be accurately solved using a high-order Earth gravity field model; for solar and lunar gravity, it is also easy to determine using the ephemeris or solar and lunar orbital parameters; for solar pressure, it is only necessary to determine the lighting conditions, the satellite area-mass ratio, and the satellite surface reflection coefficient to solve it. Atmospheric perturbations are the most complex, and there is great uncertainty in the atmospheric density, atmospheric drag coefficient, and atmospheric movement speed that are directly related to the perturbation acceleration:

[0003] 1) The atmospheric density ρ is greatly affected by various space environments (such as geomagnetic activity and solar activity). There are also factors such as diurnal terms and semi-diurnal terms. It can only be estimated through empirical formulas, and it is difficult to get an accurate value.

[0004] 2) Atmospheric drag coefficient C d Its size is difficult to estimate accurately, as it is affected by the shape of the spacecraft, its surface material, and the composition of the atmosphere and changes in real time. However, in engineering, a constant value of 2.2 is usually taken, which causes errors in the calculation of atmospheric drag perturbations.

[0005] 3) The motion of the atmosphere at high altitudes (above 100 km) is complex and its mechanism is unclear. Orbital perturbation analysis usually simplifies it and assumes that the atmosphere rotates synchronously with the earth, which also causes errors in the calculation of the relative motion speed between the satellite and the atmosphere. Summary of the invention

[0006] The technical problem solved by the present invention is to overcome the deficiencies of the prior art and provide a satellite atmospheric perturbation measurement method and device.

[0007] The technical solution of the present invention is:

[0008] In a first aspect, an embodiment of the present invention provides a satellite atmospheric perturbation measurement method, comprising:

[0009] Based on Beidou positioning and timing functions, obtain the satellite coordinates, satellite speed and current time of the satellite;

[0010] Acquiring satellite measurement attitude information and solar panel rotation angle information of the satellite;

[0011] Determine the sum of perturbations of the satellite in an inertial coordinate system according to the satellite coordinates, the satellite speed and the current time;

[0012] Determining a non-spherical gravitational perturbation corresponding to the satellite according to the satellite position and a third-party gravitational field model;

[0013] Determine the daily gravitational perturbation and the lunar gravitational perturbation corresponding to the satellite according to the solar coordinates and the lunar coordinates at the current moment;

[0014] Determine the solar light pressure perturbation corresponding to the satellite according to the current time, the satellite coordinates, the sun coordinates, the satellite measurement attitude information and the solar sail panel rotation angle information;

[0015] The satellite atmospheric perturbation of the satellite is determined according to the perturbation sum, the non-spherical gravitational perturbation, the daily gravitational perturbation, the lunar gravitational perturbation, and the solar light pressure perturbation.

[0016] Optionally, determining the sum of perturbations of the satellite in an inertial coordinate system according to the satellite coordinates, the satellite speed and the current time includes:

[0017] Obtain the J2000 epoch of the current moment;

[0018] According to the J2000 epoch at the current moment, a coordinate system transfer matrix between the earth-fixed system and the J2000 geocentric equator is calculated;

[0019] Calculate the target position vector and target velocity vector of the satellite in the J2000 geocentric equatorial coordinate system according to the coordinate system transfer matrix, the position vector corresponding to the satellite coordinates, and the velocity vector corresponding to the satellite velocity;

[0020] Calculating six satellite orbit numbers of the satellite according to the target position vector and the target velocity vector;

[0021] According to the six satellite orbit numbers and the current time, the theoretical position of the satellite under two-body motion is calculated;

[0022] The perturbation sum of the satellite in the inertial coordinate system is calculated according to the satellite theoretical position, the satellite position of the satellite in the earth fixed system and the current time.

[0023] Optionally, determining the non-spherical gravitational perturbation corresponding to the satellite according to the satellite position and a third-party gravitational field model includes:

[0024] According to the satellite position, calculate the earth's gravity of the satellite under two-body motion, and the satellite longitude, satellite latitude and distance from the center of the earth corresponding to the satellite;

[0025] Calculate the reference earth gravity corresponding to the satellite according to the satellite longitude, the satellite latitude, the distance from the center of the earth and the third-party gravity field model;

[0026] The non-spherical gravitational perturbation corresponding to the satellite is calculated according to the earth's gravity and the reference earth's gravity.

[0027] Optionally, determining the daily gravitational perturbation and the lunar gravitational perturbation corresponding to the satellite according to the solar coordinates and the lunar coordinates at the current moment includes:

[0028] According to the current time, the solar coordinates of the sun at the current time and the lunar coordinates of the moon at the current time are calculated;

[0029] According to the satellite theoretical position under two-body motion, the solar coordinates and the lunar coordinates, the daily gravitational perturbation and the lunar gravitational perturbation corresponding to the satellite are calculated.

[0030] Optionally, determining the satellite atmospheric perturbation of the satellite according to the perturbation sum, the non-spherical gravitational perturbation, the daily gravitational perturbation, the lunar gravitational perturbation, and the solar light pressure perturbation includes:

[0031] Calculating and obtaining the perturbation sum value between the daily gravitational perturbation, the lunar gravitational perturbation and the solar light pressure perturbation;

[0032] The difference between the perturbation sum and the perturbation sum is calculated, and the difference is used as the satellite atmosphere perturbation.

[0033] In a second aspect, an embodiment of the present invention provides a satellite atmospheric perturbation measurement device, comprising:

[0034] Satellite information acquisition module, used to obtain the satellite coordinates, satellite speed and current time of the satellite based on Beidou positioning and timing functions;

[0035] A satellite attitude acquisition module, used to acquire satellite measurement attitude information and solar panel rotation angle information of the satellite;

[0036] A perturbation sum determination module, used to determine the perturbation sum of the satellite in an inertial coordinate system according to the satellite coordinates, the satellite speed and the current time;

[0037] A non-spherical gravitational perturbation determination module, used to determine the non-spherical gravitational perturbation corresponding to the satellite according to the satellite position and a third-party gravitational field model;

[0038] A solar-lunar gravitational perturbation determination module, used to determine the solar gravitational perturbation and lunar gravitational perturbation corresponding to the satellite according to the solar coordinates and the lunar coordinates at the current moment;

[0039] A solar light pressure perturbation determination module, configured to determine the solar light pressure perturbation corresponding to the satellite according to the current time, the satellite coordinates, the sun coordinates, the satellite measurement attitude information and the solar sail panel rotation angle information;

[0040] The satellite atmospheric perturbation determination module is used to determine the satellite atmospheric perturbation of the satellite according to the perturbation sum, the non-spherical gravitational perturbation, the daily gravitational perturbation, the lunar gravitational perturbation, and the solar light pressure perturbation.

[0041] Optionally, the perturbation sum determination module includes:

[0042] A J2000 epoch acquisition unit, used to acquire the J2000 epoch of the current moment;

[0043] A coordinate system transfer matrix calculation unit, used for calculating the coordinate system transfer matrix between the earth fixed system and the J2000 geocentric equator according to the J2000 epoch at the current moment;

[0044] a target position velocity vector calculation unit, configured to calculate a target position vector and a target velocity vector of the satellite in a J2000 geocentric equatorial coordinate system according to the coordinate system transfer matrix, the position vector corresponding to the satellite coordinates, and the velocity vector corresponding to the satellite velocity;

[0045] A satellite orbit six-element number calculation unit, used to calculate the satellite orbit six-element number of the satellite according to the target position vector and the target velocity vector;

[0046] A satellite theoretical position calculation unit, used for calculating the satellite theoretical position of the satellite under two-body motion according to the six satellite orbit numbers and the current time;

[0047] The perturbation sum calculation unit is used to calculate the perturbation sum of the satellite in the inertial coordinate system according to the satellite theoretical position, the satellite position of the satellite in the earth fixed system and the current time.

[0048] Optionally, the non-spherical gravitational perturbation determination module includes:

[0049] An earth gravity calculation unit, used to calculate the earth gravity of the satellite under two-body motion, and the satellite longitude, satellite latitude and distance from the center of the earth corresponding to the satellite according to the satellite position;

[0050] A reference earth gravity calculation unit, used to calculate the reference earth gravity corresponding to the satellite according to the satellite longitude, the satellite latitude, the distance from the center of the earth and the third-party gravity field model;

[0051] The non-spherical gravitational perturbation calculation unit is used to calculate the non-spherical gravitational perturbation corresponding to the satellite according to the earth's gravity and the reference earth's gravity.

[0052] Optionally, the solar and lunar gravitational perturbation determination module includes:

[0053] A sun-moon coordinate calculation unit, used to calculate the sun coordinates of the sun at the current moment and the moon coordinates of the moon at the current moment according to the current moment;

[0054] The solar-lunar gravitational perturbation calculation unit is used to calculate the solar gravitational perturbation and lunar gravitational perturbation corresponding to the satellite according to the satellite theoretical position under two-body motion, the solar coordinates and the lunar coordinates.

[0055] Optionally, the satellite atmospheric perturbation determination module includes:

[0056] A perturbation sum value calculation unit, used for calculating the perturbation sum value between the daily gravitational perturbation, the lunar gravitational perturbation and the solar light pressure perturbation;

[0057] The satellite atmospheric perturbation acquisition unit is used to calculate the difference between the perturbation sum and the perturbation sum value, and use the difference as the satellite atmospheric perturbation.

[0058] The advantages of the present invention compared with the prior art are:

[0059] 1. By utilizing the timing and navigation functions of the BeiDou system and jointly solving multiple perturbation factors, the calculation errors caused by atmospheric uncertainty are avoided, and accurate calculation of atmospheric perturbations is achieved;

[0060] 2. By utilizing the short message communication function of the Beidou system, the data acquisition link is freed from the dependence on satellite telemetry channels, achieving full-arc and all-weather data acquisition and real-time calculation of perturbations;

[0061] 3. By designing the computing platform software, the full process of automatic solution of orbital atmospheric perturbations is realized. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 A flowchart of a method for measuring satellite atmospheric perturbations provided by an embodiment of the present invention;

[0063] Figure 2 A schematic diagram of a satellite atmospheric perturbation measurement system provided by an embodiment of the present invention;

[0064] Figure 3 A schematic diagram of an atmospheric perturbation calculation process in a computing platform provided by an embodiment of the present invention;

[0065] Figure 4 A schematic diagram of a computing platform provided by an embodiment of the present invention;

[0066] Figure 5 A schematic diagram of satellite orbit elements in an equatorial inertial coordinate system provided by an embodiment of the present invention;

[0067] Figure 6 A schematic diagram of the structure of a satellite atmospheric perturbation measurement device provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0068] Embodiment 1

[0069] Reference Figure 1 , shows a flowchart of the steps of a satellite atmosphere measurement method provided by an embodiment of the present invention, such as Figure 1 As shown, the satellite atmosphere measurement method may include the following steps:

[0070] Step 101: Based on the Beidou positioning and timing functions, obtain the satellite coordinates, satellite speed and current time of the satellite.

[0071] The embodiment of the present invention is driven by information such as coordinates, speed, time and the satellite's own measured attitude obtained by the satellite through Beidou navigation, and at the same time utilizes the short message communication function of the Beidou system to ensure real-time feedback of information, and builds a software platform to perform full-process calculations on the atmospheric perturbations experienced by the satellite.

[0072] The testing method principle of the present invention is as follows Figure 2 As shown, there are 6 data streams in total. The test process is carried out in the order of the data stream numbers. The specific process is as follows:

[0073] Data stream 1 is four initial information obtained by the Beidou terminal and the satellite attitude control system: satellite coordinates, satellite speed, current time, satellite measured attitude and solar panel rotation angle. These four information are the inputs of data streams 2 to 6. The satellite-borne Beidou navigation receiver obtains three pieces of information: satellite coordinates, satellite speed and current time. The satellite attitude control system obtains satellite measured attitude and solar panel rotation angle information. The satellite-borne Beidou short message terminal transmits these four pieces of information to the ground Beidou short message terminal.

[0074] The input of data stream 2 is satellite coordinates, satellite speed, and current time, and the output is the total perturbation of the satellite, which serves as the input of data stream 6.

[0075] The input of data stream 3 is the satellite coordinates, and the output is the non-spherical gravitational perturbation exerted on the satellite, which serves as the input of data stream 6.

[0076] The input of data stream 4 is the satellite coordinates and the current time, and the output is the gravitational perturbations of the sun and the moon acting on the satellite, which serves as the input of data stream 6.

[0077] The input of data stream 5 is the satellite coordinates, current time, satellite measurement attitude and solar panel rotation angle, and the output is the sunlight pressure perturbation on the satellite, which serves as the input of data stream 6.

[0078] The input of data stream 6 is the sum of satellite perturbations, non-spherical gravitational perturbations, solar and lunar gravitational perturbations, and solar pressure perturbations (i.e., the output of data streams 2, 3, 4, and 5), and the output is the satellite's atmospheric perturbations.

[0079] The principle of the test method provided in this embodiment can be as follows Figure 2 As shown, it is divided into two parts: data acquisition and computing platform. The data acquisition part can include: a satellite-borne Beidou receiver and a satellite attitude control system. The satellite-borne Beidou receiver can use the Beidou positioning and timing functions to obtain the satellite coordinates, satellite speed and current time of the satellite.

[0080] Step 102: Obtain satellite measurement attitude information and solar panel rotation angle information of the satellite.

[0081] The satellite attitude control system can obtain satellite measurement attitude information and solar panel rotation angle information.

[0082] like Figure 2 As shown, after the satellite coordinates, satellite speed and current time of the satellite are obtained by the onboard Beidou receiver, and the satellite measurement attitude information and solar panel angle information are obtained by the satellite attitude control system, the onboard Beidou receiver can send the satellite coordinates, satellite speed and current time to the onboard Beidou short message terminal, and the satellite attitude control system can send the satellite measurement attitude information and solar panel angle information to the onboard Beidou short message terminal, and then, the onboard Beidou short message terminal can send it to the ground Beidou short message terminal via a short message, and the ground Beidou short message terminal can send this information to the computing platform.

[0083] In this example, the description of satellite coordinates and speed is 3 variables, attitude quaternion is 4 variables, current time and sailboard angle are 1 variable, to ensure sufficient accuracy, satellite coordinates, speed, attitude, sailboard angle are 4 bytes, current time is 8 bytes, a total of 52 bytes. The regional communication capacity of Beidou short message is 1750 bytes per time, and the global communication capacity is 70 bytes per time. The capacity of a single short message fully meets the data transmission requirements.

[0084] The embodiment of the present invention utilizes the short message communication function of the Beidou system to ensure the real-time return of information, and can realize the real-time measurement of satellite atmospheric perturbations.

[0085] After obtaining the satellite coordinates, satellite speed and current time of the satellite, as well as the satellite measurement attitude information and solar panel rotation angle information of the satellite, step 103 is executed.

[0086] Step 103: Determine the sum of perturbations of the satellite in an inertial coordinate system according to the satellite coordinates, the satellite speed and the current time.

[0087] After obtaining the satellite coordinates, satellite speed and current time of the satellite, the sum of the perturbations of the satellite in the inertial coordinate system can be determined based on the satellite coordinates, satellite speed and current time of the satellite. Specifically, it can be described in detail in combination with the following specific implementation method.

[0088] First, some formulas involved in the embodiments of the present invention are described.

[0089] like Figure 5 As shown, Υ-vernal equinox, N-ascending node of satellite orbit, S-satellite position, P-perigee of satellite orbit, Ω-ascending node right ascension, ω-argument of perigee, f-true anomaly, i-orbital inclination, e-eccentricity vector, pointing from the center of the earth to the perigee, with a length of e, W-unit vector of the normal to the orbital plane, and r-satellite radius.

[0090] In a specific implementation of the present invention, the above step 103 may include:

[0091] Sub-step A1: Obtain the J2000 epoch of the current moment.

[0092] In the embodiment of the present invention, after the computing platform receives the above information of the satellite, it can perform corresponding calculations. The calculation process can be as follows: Figure 2 The right part, Figure 3 and Figure 4 shown.

[0093] First, we can remember the initial data: the satellite coordinates in the ground-fixed system are r e =(x e ,y e ,z e ) T , the satellite speed in the ground-fixed system is The current time is t0.

[0094] You can get the current epoch of J2000. Specifically, the Julian century number of the current moment, starting from the standard epoch of J2000, is calculated using the following formula:

[0095]

[0096] JD(t0) is the Julian day corresponding to time t0.

[0097] Sub-step A2: According to the J2000 epoch at the current moment, a coordinate system transfer matrix between the Earth-fixed system and the J2000 geocentric equator is calculated.

[0098] After obtaining the J2000 epoch at the current moment, the coordinate system transfer matrix between the earth-fixed system and the J2000 geocentric equator can be calculated according to the J2000 epoch at the current moment. Specifically, the calculation process can refer to the following formula.

[0099] The expression of the vector or coordinate p in the earth-fixed system to the inertial system p' is:

[0100] p'=(PR) T (NR) T (ER) T (EP) T p=(IR)p (2-2)

[0101] In the above formula, EP is the polar motion, ER is the direction of the Earth's rotation, NR is the nutation, and PR is the precession. The expressions of each matrix are

[0102] as follows:

[0103] EP=R y (-x p )R x (-y p ) (2-3)

[0104] ER=R z (S G ) (2-4)

[0105] NR=R x (-Δε)R y (Δθ)R z (-Δμ) (2-5)

[0106] PR=R z (-z A )R y (θ A )R z (-ξ A ) (2-6)

[0107] Where R x (θ), R y (θ), R z (θ) is the coordinate transfer matrix, which is expressed as:

[0108]

[0109] Since the polar shift generally does not exceed 0".8, its influence is ignored in the calculation. In formula (5), S G is the true sidereal hour angle of Greenwich, and its expression is:

[0110]

[0111] In formula (2-8), Δμ is the precession and nutation of right ascension, and the calculation formula is shown in formula (13): For Greenwich mean sidereal time, the expression is as follows:

[0112]

[0113] The precession model adopts IAU1976. In formula (2-6), ξ A 、z A ,θ A is the precession angle, which is expressed as follows:

[0114]

[0115] Nutation includes ecliptic longitude nutation ΔΨ, right ascension nutation Δμ, declination nutation Δθ, and angle nutation Δε. The nutation sequence adopts IAU1980, and there are 106 items related to the gravity of the sun and the moon, and its expression is as follows:

[0116]

[0117] In the sequence of formula (2-11), α i The five basic variables represented by (t) are the lunar mean anomaly l, the solar mean anomaly l', the lunar mean ascending node angle F, the lunar mean angular distance D and the ecliptic longitude Ω of the lunar orbit ascending node. The expressions are as follows:

[0118]

[0119] The corresponding right ascension and declination nutation Δμ and Δθ are:

[0120] Δμ=Δψcosε (2-13)

[0121] Δθ=Δψsinε (2-14)

[0122] The calculation formula of the obliquity of the ecliptic is as follows:

[0123] ε=84381.448”-46.8150”t-0.00059”t 2 +0.001813”t 3 (2-15)

[0124] In formulas (2-9) to (2-15), t is obtained from (2-1).

[0125] Sub-step A3: Calculate the target position vector and target velocity vector of the satellite in the J2000 geocentric equatorial coordinate system based on the coordinate system transfer matrix, the position vector corresponding to the satellite coordinates, and the velocity vector corresponding to the satellite velocity.

[0126] The target position vector refers to the position vector of the satellite in the J2000 geocentric equatorial coordinate system.

[0127] The target velocity vector refers to the velocity vector of the satellite in the J2000 geocentric equatorial coordinate system.

[0128] After the coordinate system transfer matrix is ​​calculated, the target position vector and target velocity vector of the satellite in the J2000 geocentric equatorial coordinate system can be calculated based on the coordinate system transfer matrix, the position vector corresponding to the satellite coordinates, and the velocity vector corresponding to the satellite velocity. The specific calculation process can be referred to the description of the following formula.

[0129] The satellite coordinates in the J2000 geocentric equatorial coordinate system are r i =(x i ,y i ,z i ) T , the satellite speed is From formula (2-2), we can get

[0130] v i =(IR)g(v e +v e0 ) (2-16)

[0131] r i =(IR)·r e (2-17)

[0132] v e0 is the velocity component of the earth's solid system formed by the earth's rotation,

[0133]

[0134] Sub-step A4: Calculate six satellite orbit numbers of the satellite based on the target position vector and the target velocity vector.

[0135] After the target position vector and the target velocity vector are calculated, the six satellite orbit numbers of the satellite can be calculated according to the target position vector and the target velocity vector. The specific calculation process can refer to the description of the following formula.

[0136] The satellite orbit is described by six numbers a, e, i, Ω, ω, and f. The six numbers are defined as follows: Figure 3 As shown, the position velocity vector calculated in 2-c is needed to obtain the orbital elements.

[0137] In the J2000 geocentric equatorial coordinate system, the satellite momentum is h = r i ×v i .

[0138] The mechanical energy of the satellite (kinetic energy + gravitational potential energy) is Where v i is the satellite motion speed, r i is the distance from the satellite to the center of the Earth, μ is the earth's gravitational constant, which can be taken as μ=398600.4415km 3 / s 2 The following is the formula for solving the six orbital numbers:

[0139] The semi-major axis of the orbit is determined by the mechanical energy and is

[0140] The orbital inclination is determined by the momentum, which is Where u z =(0,0,1) is the unit vector of the +Z coordinate axis.

[0141] The right ascension of the orbital node is also determined by the momentum, which is Ω = arccos (Ngu x ), where N = (u z ×h) / |u z ×h| is the orbital pitch line unit vector, u x =(0,0,1) is the unit vector of the +X coordinate axis.

[0142] The orbital eccentricity is determined by the semi-major axis and the semi-diameter, which is p is the semi-diameter of the track,

[0143] The true anomaly of the satellite is determined by the semi-diameter, the distance from the center of the earth, and the eccentricity:

[0144] The orbital perigee argument is determined by the satellite position vector and the true anomaly angle, which is

[0145] The anomaly angle can also be expressed as the partial anomaly angle or the mean anomaly angle. The relationship between the true anomaly angle and the partial anomaly angle E is:

[0146]

[0147] The mean anomaly angle can be obtained from the eccentric anomaly angle

[0148] M=E-esinE (2-20)

[0149] Sub-step A5: Calculate the theoretical position of the satellite under two-body motion according to the six satellite orbit numbers and the current time.

[0150] After calculating the six numbers of the satellite orbit, the theoretical position of the satellite under two-body motion (i.e. the actual position of the satellite) can be calculated based on the six numbers of the satellite orbit and the current time. The specific calculation process can be referred to the description of the following formula.

[0151] Two-body motion treats both the earth and the satellite as point masses. Let the time corresponding to one frame of input data be t0, and the time corresponding to the next frame of data be t1. Taking the satellite positioning data at time t0 as the starting point (input), ignoring all perturbations, the theoretical satellite position at time t1 is calculated according to the two-body motion law.

[0152] The orbital parameters of the satellite at time t0 calculated according to the above steps are recorded as a, e, i, Ω, ω, f0, E0, M0. Then the orbital parameters at time t1 under two-body motion are a, e, i, Ω, ω, f b , E b ,M b , it can be seen that the five orbital parameters remain unchanged, and the method for calculating the anomaly angle is as follows.

[0153] Mean anomaly angle M b For uniform change:

[0154]

[0155] Eccentric angle E b Obtained by multiple iterations of the mean anomaly angle:

[0156] E1=M b +esinE0

[0157] E2=M b +esinE1

[0158] L

[0159] E n =M b +esinE n-1

[0160] E b =M b +esinE n (2-22)

[0161] True anomaly f b The true anomaly angle can be derived from the eccentric anomaly angle according to formula (2-19) to derive the theoretical coordinates of the satellite in the J2000 geocentric equatorial coordinate system: ib for:

[0162]

[0163] Sub-step A6: Calculate the sum of the perturbations of the satellite in the inertial coordinate system according to the theoretical position of the satellite, the position of the satellite in the earth-fixed system and the current time.

[0164] After the satellite theoretical position is calculated, the sum of the perturbations of the satellite in the inertial coordinate system can be calculated based on the satellite theoretical position, the satellite position in the earth fixed system and the current time. The specific calculation process can be described in the following formula.

[0165] The actual position of the satellite obtained by the Beidou system at time t1 is r i1 =[x i1 ,y i1 ,z i1 ] T

[0166] The sum of the perturbation acceleration in the inertial coordinate system is W all =[x all ,y all ,z all ] T , then

[0167]

[0168] Step 104: Determine the non-spherical gravitational perturbation corresponding to the satellite according to the satellite position and the third-party gravitational field model.

[0169] After obtaining the satellite position, the non-spherical gravitational perturbation corresponding to the satellite can be determined according to the satellite position and the third-party gravitational field model. Specifically, it can be combined with the description of the following specific implementation method.

[0170] In another specific implementation of the present invention, the above step 104 may include:

[0171] Sub-step B1: According to the satellite position, calculate the earth's gravity of the satellite under two-body motion, and the satellite longitude, satellite latitude and distance from the center of the earth corresponding to the satellite.

[0172] In an embodiment of the present invention, after the satellite position of the satellite is obtained, the earth's gravity of the satellite under two-body motion, and the satellite longitude, satellite latitude and distance from the center of the earth corresponding to the satellite can be calculated based on the satellite position. Specifically, the description of the following formula can be referred to.

[0173] 1. Calculate the Earth's gravity under two-body motion

[0174] Under two-body motion, the simplified gravitational force G on the satellite is a The components in each direction are: radial East (longitude) G aλ = 0, north (latitude) μ is still the gravitational constant of the Earth.

[0175] 2. Calculate satellite longitude, latitude, and distance from the center of the earth

[0176] From the satellite coordinates fixed on the ground, the satellite longitude is

[0177] The latitude of the Earth's center is

[0178] The distance from the center of the earth is

[0179] Sub-step B2: Calculate the reference earth gravity corresponding to the satellite based on the satellite longitude, the satellite latitude, the distance from the center of the earth and the third-party gravity field model.

[0180] After calculating the satellite longitude, satellite latitude, distance from the geocenter and the earth's gravity, the reference earth's gravity corresponding to the satellite, that is, the precise earth's gravity of the satellite, can be calculated based on the satellite longitude, satellite latitude, distance from the geocenter and the earth's gravity in combination with a third-party gravity field model. The specific calculation process can be referred to the description of the following formula.

[0181] According to the normalized celestial body gravitational field potential function and gravitational field model, the precise earth gravity acting on the satellite can be written as follows:

[0182]

[0183] Where μ is the Earth's gravitational constant a e is the equatorial radius of the reference ellipsoid, is the normalized Legendre polynomial, and its recursive formula is:

[0184]

[0185] The order lmax and coefficients of the gravitational field potential function and Provided by the adopted gravitational field model, the currently publicly available gravitational field models range from dozens to thousands of orders and can be selected according to needs.

[0186] Sub-step B3: Calculate the non-spherical gravitational perturbation corresponding to the satellite based on the earth's gravity and the reference earth's gravity.

[0187] After obtaining the earth's gravity and the reference earth's gravity, the non-spherical gravitational perturbation corresponding to the satellite can be calculated based on the earth's gravity and the reference earth's gravity. The specific calculation process can be described in the following formula.

[0188] According to the partial derivative of the gravitational field function in the above steps, we can get the precise earth gravity W on the satellite. a Components in the radial, easting, and northing directions:

[0189]

[0190]

[0191]

[0192] The derivatives involving Legendre polynomials have the following recursive formula:

[0193]

[0194] The gravitational perturbation of the satellite is W ga =W a -G a , its radial, eastward, and northward components are:

[0195]

[0196]

[0197]

[0198] Step 105: Determine the daily gravitational perturbation and the lunar gravitational perturbation corresponding to the satellite according to the solar coordinates and the lunar coordinates at the current moment.

[0199] After obtaining the current time, the solar coordinates and lunar coordinates of the current time can be obtained, and the daily gravitational perturbation and lunar gravitational perturbation corresponding to the satellite can be determined based on the solar coordinates and lunar coordinates. Specifically, it can be described in detail in combination with the following specific implementation method.

[0200] In another specific implementation of the present invention, the above step 105 may include:

[0201] Sub-step C1: According to the current moment, the solar coordinates of the sun at the current moment and the lunar coordinates of the moon at the current moment are calculated.

[0202] Sub-step C2: Calculate the daily gravitational perturbation and the lunar gravitational perturbation corresponding to the satellite according to the theoretical position of the satellite under two-body motion, the solar coordinates and the lunar coordinates.

[0203] In the embodiment of the present invention, there are two methods for calculating the sun and moon coordinates: using the average orbital parameters of the moon and the earth to obtain approximate coordinates or reading the astronomical ephemeris. The former has the advantage of occupying very little space when implemented in software, and the latter has the advantage of being more accurate. The average orbital parameter method is selected here.

[0204] Among the lunar average orbital elements used, the calculation formulas for the semi-major axis, eccentricity, inclination, and ascending node ecliptic longitude are:

[0205] a l =383397.7916+0.0038t

[0206] e l =0.055545526-0.000000016t

[0207] i l =5°.15668983-0″.00008t+0″.02966t 2 -0″.000042t 3 -0″.00000013t 4

[0208] Ω l =125°.04455504-6967919″.3622t+6″.3622t 2 +0″.007625t 3 -0″.00003586t 4

[0209] The calculation formulas for the lunar perigee argument and mean anomaly are:

[0210]

[0211] in,

[0212]

[0213] λ l =218°.31665436+1732559343″.73604t-5″.8883t 2 +0″.006604t 3 -0″.00003169t 4

[0214] The time unit is still the Julian century

[0215] The reference system is the J2000 geocentric ecliptic coordinate system. The coordinates of the moon in this coordinate system are solved by the orbital elements. Refer to formulas (2-19), (2-22), (2-23) to obtain the coordinates of the moon as r l0 =(x l0 ,y l0 ,z l0 ) T .

[0216] The coordinates of the moon in the J2000 geocentric equatorial coordinate system are r l =(x l ,y l ,z l ) T =r l0 R x (-ε0), where ε0 is the obliquity of the mean ecliptic at J2000.0, ε0=23°26′21″.448, R x The function is shown in formula 2-7.

[0217] In the average orbital elements of the Earth around the Sun, the coordinate system is the J2000 heliocentric ecliptic coordinate system, and the calculation formulas for the semi-major axis, eccentricity, inclination, and ascending node ecliptic longitude are:

[0218] a e =1.0000010178au

[0219] e e =0.0167086342-0.0004203654t-0.0000126734t 2 +144410 -10 t 3 -210 -10 t 4 +310 - 10 t 5

[0220] i e =469″.97289t-3″.35053t 2 -0″.12374t 3 +0″.00027t 4 -0″.00001t 5 +0″.00001t 6

[0221] Ω e =174°.87317577-8679″.27034t+15″.34191t 2 +0″.00532t 3 -0″.03734t 4 -0″.00073t 5 +0″.00004t 6

[0222] Among them, au is the astronomical unit, 1au = 1.49597870710 8 km.

[0223] The calculation formulas for the argument of perigee and mean anomaly are:

[0224]

[0225] in:

[0226]

[0227] λ e =100°.46645683+1295977422″.83429t-2″.04411t 2 -0″.00523t 3

[0228] The time unit here is the Julian millennium

[0229] The coordinates of the Earth in the J2000 heliocentric ecliptic coordinate system are solved by the orbital elements. Refer to formulas (2-19), (2-22), (2-23) to obtain the coordinates of the Earth as r e0 =(x e0 ,y e0 ,z e0 ) T The sun's coordinates in the J2000 geocentric ecliptic coordinate system are r s0 =(x s0 ,y s0 ,z s0 ) T =-r e0 The solar coordinates in the J2000 geocentric equatorial coordinate system are

[0230] r s =(x s ,y s ,z s ) T =r s0 R x (-ε0).

[0231] According to the above steps, the vector is:

[0232] Vr s =r i -r s , Vr l =r i -r l (2-30)

[0233] The gravitational perturbation of the sun is:

[0234]

[0235] The gravitational perturbation of the moon is:

[0236]

[0237] G is the gravitational constant, m s is the mass of the sun, m l is the mass of the moon. The gravitational perturbation of the sun and the moon is:

[0238] W sl =W s +W l (2-33)

[0239] Step 106: Determine the solar pressure perturbation corresponding to the satellite according to the current time, the satellite coordinates, the sun coordinates, the satellite measurement attitude information and the solar sail panel rotation angle information.

[0240] After obtaining the relevant parameters of the satellite, the solar pressure perturbation corresponding to the satellite can be determined according to the current time, satellite coordinates, sun coordinates, satellite measurement attitude information and solar sail angle information, and the additional push:

[0241] The satellite attitude quaternion relative to the J2000 geocentric equatorial coordinate system is q = [q1, q2, q3, q4] T , the solar sail panel rotation angle is a p is the satellite information obtained through Beidou short message. The area of ​​the solar sail is S p , the reflection coefficient is η p , is a known quantity. Assume that the sailboard rotation axis is parallel to the +Y axis of the star, and when the sailboard rotation angle is 0, its normal is parallel to the -Z axis of the star.

[0242] The satellite's attitude matrix is:

[0243]

[0244] In the J2000 geocentric equatorial coordinate system, the unit vector of the satellite solar sail normal is:

[0245]

[0246] The unit vector of the direction of the sun's source is:

[0247]

[0248] The angle between the normal direction of the sailboard and the direction of the sun light source (sunlight incident angle) is:

[0249]

[0250] The solar radiation pressure is taken as ρ f =4.5×10 -6 N / m.

[0251] The solar pressure on the sailboard is:

[0252]

[0253] The approximate equivalent cross-sectional area of ​​a star is:

[0254] S x S y S z are the cross-sectional areas of the XYZ plane of the star, and the reflection coefficient is η s , the solar pressure on the star is estimated to be:

[0255]

[0256] Satellite solar panels have a large illuminated area and are always pointed toward the sun, absorbing and reflecting a high amount of solar energy, so they are the main part of the solar light pressure perturbation. The illumination of a star is greatly affected by its shape, shielding by the panels, and its posture, and the illuminated area is usually significantly smaller than the panels, so a simplified estimate of its light pressure is made. The magnitude of the solar light pressure perturbation of a low-orbit satellite is usually 10 -8 , which is much lower than the non-spherical gravity (10 -3 )、Sun and Moon Gravitation(10 -7 ), atmospheric resistance (10 -6 ), and the light pressure on the star does not account for a major part of the total light pressure, so this estimate is acceptable.

[0257] If the satellite mass is known to be m s , then the solar pressure perturbation received by the satellite is

[0258]

[0259] In the specific implementation, this calculation is performed only when the satellite is illuminated. If the satellite is in the shadow of the Earth or the Moon, the solar pressure perturbation W f =0.

[0260] Step 107: Determine the satellite atmospheric perturbation of the satellite according to the perturbation sum, the non-spherical gravitational perturbation, the daily gravitational perturbation, the lunar gravitational perturbation, and the solar pressure perturbation.

[0261] After the perturbation sum, non-spherical gravitational perturbation, daily gravitational perturbation, lunar gravitational perturbation, and solar light pressure perturbation are obtained through the above steps, the satellite atmospheric perturbation of the satellite can be determined according to the perturbation sum, non-spherical gravitational perturbation, daily gravitational perturbation, lunar gravitational perturbation, and solar light pressure perturbation. Specifically, the perturbation sum value between the daily gravitational perturbation, lunar gravitational perturbation, and solar light pressure perturbation can be calculated, and then the difference between the perturbation sum and the perturbation sum value is calculated, and the difference is used as the satellite atmospheric perturbation. The specific calculation process can be shown in the following formula:

[0262] W A =W all -W ga -W sl -W f (2-37)

[0263] The embodiment of the present invention utilizes the timing and navigation functions of the Beidou system and jointly solves multiple perturbation factors to avoid calculation errors caused by atmospheric uncertainty and achieve accurate calculation of atmospheric perturbations.

[0264] Embodiment 2

[0265] Reference Figure 6 , shows a schematic diagram of the structure of a satellite atmosphere measurement device provided by an embodiment of the present invention, such as Figure 6 As shown, the satellite atmosphere measurement device may include the following modules:

[0266] The satellite information acquisition module 610 is used to obtain the satellite coordinates, satellite speed and current time of the satellite based on the Beidou positioning and timing functions;

[0267] A satellite attitude acquisition module 620 is used to acquire satellite measurement attitude information and solar panel rotation angle information of the satellite;

[0268] A perturbation sum determination module 630, configured to determine a perturbation sum of the satellite in an inertial coordinate system according to the satellite coordinates, the satellite speed and the current time;

[0269] A non-spherical gravitational perturbation determination module 640, configured to determine the non-spherical gravitational perturbation corresponding to the satellite according to the satellite position and a third-party gravitational field model;

[0270] A solar-lunar gravitational perturbation determination module 650 is used to determine the solar gravitational perturbation and the lunar gravitational perturbation corresponding to the satellite according to the solar coordinates and the lunar coordinates at the current moment;

[0271] A solar light pressure perturbation determination module 660, configured to determine the solar light pressure perturbation corresponding to the satellite according to the current time, the satellite coordinates, the sun coordinates, the satellite measurement attitude information and the solar sail panel rotation angle information;

[0272] The satellite atmospheric perturbation determination module 670 is used to determine the satellite atmospheric perturbation of the satellite according to the perturbation sum, the non-spherical gravitational perturbation, the daily gravitational perturbation, the lunar gravitational perturbation, and the solar pressure perturbation.

[0273] Optionally, the perturbation sum determination module includes:

[0274] A J2000 epoch acquisition unit, used to acquire the J2000 epoch of the current moment;

[0275] A coordinate system transfer matrix calculation unit, used for calculating the coordinate system transfer matrix between the earth fixed system and the J2000 geocentric equator according to the J2000 epoch at the current moment;

[0276] a target position velocity vector calculation unit, configured to calculate a target position vector and a target velocity vector of the satellite in a J2000 geocentric equatorial coordinate system according to the coordinate system transfer matrix, the position vector corresponding to the satellite coordinates, and the velocity vector corresponding to the satellite velocity;

[0277] A satellite orbit six-element number calculation unit, used for calculating the satellite orbit six-element number of the satellite according to the target position vector and the target velocity vector;

[0278] A satellite theoretical position calculation unit, used for calculating the satellite theoretical position of the satellite under two-body motion according to the six satellite orbit numbers and the current time;

[0279] The perturbation sum calculation unit is used to calculate the perturbation sum of the satellite in the inertial coordinate system according to the satellite theoretical position, the satellite position of the satellite in the earth fixed system and the current time.

[0280] Optionally, the non-spherical gravitational perturbation determination module includes:

[0281] An earth gravity calculation unit, used to calculate the earth gravity of the satellite under two-body motion, and the satellite longitude, satellite latitude and distance from the center of the earth corresponding to the satellite according to the satellite position;

[0282] A reference earth gravity calculation unit, used to calculate the reference earth gravity corresponding to the satellite according to the satellite longitude, the satellite latitude, the distance from the center of the earth and the third-party gravity field model;

[0283] The non-spherical gravitational perturbation calculation unit is used to calculate the non-spherical gravitational perturbation corresponding to the satellite according to the earth's gravity and the reference earth's gravity.

[0284] Optionally, the solar and lunar gravitational perturbation determination module includes:

[0285] A sun-moon coordinate calculation unit, used to calculate the sun coordinates of the sun at the current moment and the moon coordinates of the moon at the current moment according to the current moment;

[0286] The solar-lunar gravitational perturbation calculation unit is used to calculate the solar gravitational perturbation and lunar gravitational perturbation corresponding to the satellite according to the satellite theoretical position under two-body motion, the solar coordinates and the lunar coordinates.

[0287] Optionally, the satellite atmospheric perturbation determination module includes:

[0288] A perturbation sum value calculation unit, used for calculating the perturbation sum value between the daily gravitational perturbation, the lunar gravitational perturbation and the solar light pressure perturbation;

[0289] The satellite atmospheric perturbation acquisition unit is used to calculate the difference between the perturbation sum and the perturbation sum value, and use the difference as the satellite atmospheric perturbation.

[0290] The specific implementation methods described in this application can enable those skilled in the art to understand this application more comprehensively, but do not limit this application in any way. Therefore, those skilled in the art should understand that this application can still be modified or replaced by equivalents; and all technical solutions and improvements that do not deviate from the spirit and technical essence of this application should be included in the protection scope of the patent of this application.

[0291] The contents not described in detail in the specification of the present invention belong to the common knowledge of those skilled in the art.

Claims

1. A satellite atmospheric perturbation measurement method, characterized in that: include: Based on Beidou positioning and timing functions, obtain the satellite coordinates, satellite speed and current time of the satellite; Acquiring satellite measurement attitude information and solar panel rotation angle information of the satellite; Determine the sum of perturbations of the satellite in an inertial coordinate system according to the satellite coordinates, the satellite speed and the current time; Determining a non-spherical gravitational perturbation corresponding to the satellite according to the satellite coordinates and a third-party gravitational field model; Determine the daily gravitational perturbation and the lunar gravitational perturbation corresponding to the satellite according to the solar coordinates and the lunar coordinates at the current moment; Determine the solar light pressure perturbation corresponding to the satellite according to the current time, the satellite coordinates, the sun coordinates, the satellite measurement attitude information and the solar sail panel rotation angle information; Determine the satellite atmosphere perturbation of the satellite according to the perturbation sum, the non-spherical gravitational perturbation, the daily gravitational perturbation, the lunar gravitational perturbation, and the solar light pressure perturbation; The step of determining the sum of perturbations of the satellite in an inertial coordinate system according to the satellite coordinates, the satellite speed and the current time includes: Obtain the J2000 epoch of the current moment; According to the J2000 epoch at the current moment, a coordinate system transfer matrix between the earth-fixed system and the J2000 geocentric equator is calculated; Calculate the target position vector and target velocity vector of the satellite in the J2000 geocentric equatorial coordinate system according to the coordinate system transfer matrix, the position vector corresponding to the satellite coordinates, and the velocity vector corresponding to the satellite velocity; Calculating six satellite orbit numbers of the satellite according to the target position vector and the target velocity vector; According to the six satellite orbit numbers and the current time, the theoretical position of the satellite under two-body motion is calculated; Calculate the sum of the perturbations of the satellite in the inertial coordinate system according to the theoretical position of the satellite, the satellite position in the earth-fixed system and the current time; The step of determining the non-spherical gravitational perturbation corresponding to the satellite according to the satellite coordinates and the third-party gravitational field model includes: According to the satellite coordinates, calculate the earth's gravity of the satellite under two-body motion, and the satellite longitude, satellite latitude and distance from the center of the earth corresponding to the satellite; Calculate the reference earth gravity corresponding to the satellite according to the satellite longitude, the satellite latitude, the distance from the center of the earth and the third-party gravity field model; Calculating a non-spherical gravitational perturbation corresponding to the satellite according to the earth's gravity and the reference earth's gravity; Determining the daily gravitational perturbation and the lunar gravitational perturbation corresponding to the satellite according to the solar coordinates and the lunar coordinates at the current moment includes: According to the current time, the solar coordinates of the sun at the current time and the lunar coordinates of the moon at the current time are calculated; According to the satellite theoretical position under two-body motion, the solar coordinates and the lunar coordinates, the daily gravitational perturbation and the lunar gravitational perturbation corresponding to the satellite are calculated.

2. The method according to claim 1, characterized in that Determining the satellite atmospheric perturbation of the satellite according to the perturbation sum, the non-spherical gravitational perturbation, the daily gravitational perturbation, the lunar gravitational perturbation, and the solar light pressure perturbation comprises: Calculating and obtaining the perturbation sum value between the daily gravitational perturbation, the lunar gravitational perturbation and the solar light pressure perturbation; The difference between the perturbation sum and the perturbation sum is calculated, and the difference is used as the satellite atmosphere perturbation.

3. A satellite atmospheric perturbation measurement device, characterized in that: include: Satellite information acquisition module, used to obtain the satellite coordinates, satellite speed and current time of the satellite based on Beidou positioning and timing functions; A satellite attitude acquisition module, used to acquire satellite measurement attitude information and solar panel rotation angle information of the satellite; A perturbation sum determination module, used to determine the perturbation sum of the satellite in an inertial coordinate system according to the satellite coordinates, the satellite speed and the current time; A non-spherical gravitational perturbation determination module, used to determine the non-spherical gravitational perturbation corresponding to the satellite according to the satellite coordinates and a third-party gravitational field model; A solar-lunar gravitational perturbation determination module, used to determine the solar gravitational perturbation and lunar gravitational perturbation corresponding to the satellite according to the solar coordinates and the lunar coordinates at the current moment; A solar light pressure perturbation determination module, used to determine the solar light pressure perturbation corresponding to the satellite according to the current time, the satellite coordinates, the sun coordinates, the satellite measurement attitude information and the solar sail panel rotation angle information; A satellite atmospheric perturbation determination module, configured to determine the satellite atmospheric perturbation of the satellite according to the perturbation sum, the non-spherical gravitational perturbation, the daily gravitational perturbation, the monthly gravitational perturbation, and the solar light pressure perturbation; The perturbation sum determination module comprises: A J2000 epoch acquisition unit, used to acquire the J2000 epoch of the current moment; A coordinate system transfer matrix calculation unit, used for calculating the coordinate system transfer matrix between the earth fixed system and the J2000 geocentric equator according to the J2000 epoch at the current moment; a target position velocity vector calculation unit, configured to calculate a target position vector and a target velocity vector of the satellite in a J2000 geocentric equatorial coordinate system according to the coordinate system transfer matrix, the position vector corresponding to the satellite coordinates, and the velocity vector corresponding to the satellite velocity; A satellite orbit six-element number calculation unit, used for calculating the satellite orbit six-element number of the satellite according to the target position vector and the target velocity vector; A satellite theoretical position calculation unit, used for calculating the satellite theoretical position of the satellite under two-body motion according to the six satellite orbit numbers and the current time; a perturbation sum calculation unit, configured to calculate the perturbation sum of the satellite in the inertial coordinate system according to the satellite theoretical position, the satellite position in the earth fixed system and the current time; The non-spherical gravitational perturbation determination module comprises: An earth gravity calculation unit, used to calculate the earth gravity of the satellite under two-body motion, and the satellite longitude, satellite latitude and distance from the center of the earth corresponding to the satellite according to the satellite coordinates; A reference earth gravity calculation unit, configured to calculate a reference earth gravity corresponding to the satellite according to the satellite longitude, the satellite latitude, the distance from the center of the earth and the third-party gravity field model; A non-spherical gravitational perturbation calculation unit, used for calculating the non-spherical gravitational perturbation corresponding to the satellite according to the earth gravity and the reference earth gravity; The solar and lunar gravitational perturbation determination module comprises: A sun-moon coordinate calculation unit, used to calculate the sun coordinates of the sun at the current moment and the moon coordinates of the moon at the current moment according to the current moment; The solar-lunar gravitational perturbation calculation unit is used to calculate the solar gravitational perturbation and lunar gravitational perturbation corresponding to the satellite according to the satellite theoretical position under two-body motion, the solar coordinates and the lunar coordinates.

4. The device according to claim 3, characterized in that The satellite atmospheric perturbation determination module comprises: A perturbation sum value calculation unit, used for calculating the perturbation sum value between the daily gravitational perturbation, the lunar gravitational perturbation and the solar light pressure perturbation; The satellite atmospheric perturbation acquisition unit is used to calculate the difference between the perturbation sum and the perturbation sum value, and use the difference as the satellite atmospheric perturbation.

Images