CN113343369B - Perturbation analysis method for spacecraft aerodynamic fusion orbit - Google Patents

Perturbation analysis method for spacecraft aerodynamic fusion orbit Download PDF

Info

Publication number
CN113343369B
CN113343369B CN202110899370.3A CN202110899370A CN113343369B CN 113343369 B CN113343369 B CN 113343369B CN 202110899370 A CN202110899370 A CN 202110899370A CN 113343369 B CN113343369 B CN 113343369B
Authority
CN
China
Prior art keywords
spacecraft
orbit
aerodynamic
initial
attitude
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110899370.3A
Other languages
Chinese (zh)
Other versions
CN113343369A (en
Inventor
高兴龙
李志辉
陈钦
彭傲平
丁娣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Equipment Design and Testing Technology Research Institute of China Aerodynamics Research and Development Center
Original Assignee
Equipment Design and Testing Technology Research Institute of China Aerodynamics Research and Development Center
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Equipment Design and Testing Technology Research Institute of China Aerodynamics Research and Development Center filed Critical Equipment Design and Testing Technology Research Institute of China Aerodynamics Research and Development Center
Priority to CN202110899370.3A priority Critical patent/CN113343369B/en
Publication of CN113343369A publication Critical patent/CN113343369A/en
Application granted granted Critical
Publication of CN113343369B publication Critical patent/CN113343369B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Abstract

The invention discloses a perturbation analysis method for a spacecraft aerodynamic fusion orbit, which relates to the field of aerospace forecasting, and the perturbation analysis method comprehensively considers the influences of a refined cross-basin aerodynamic modeling and a posture rotation effect and forecasts the information of a degradation orbit of the spacecraft in advance, overcomes the defect that the traditional technical method does not take pneumatic influence factors into consideration, can remarkably improve the forecasting precision of the large-scale uncontrolled spacecraft in an orbit flight and reentry process, can prejudge in advance, thus providing a basis for forecasting of information such as reentry time position and the like, and being beneficial to developing problem solving forecasting, hazard analysis and risk assessment of the reentry process of the large-scale uncontrolled spacecraft; meanwhile, the technology can support multidisciplinary coupling research and optimization design of the problems of uncontrolled reentry of a large spacecraft, maneuvering orbital transfer and the like.

Description

Perturbation analysis method for spacecraft aerodynamic fusion orbit
Technical Field
The invention relates to the field of space forecasting, in particular to a perturbation analysis method for a spacecraft aerodynamic fusion orbit.
Background
The process of returning and re-entering the earth atmosphere by the large spacecraft at the end of the service life can be divided into uncontrolled re-entry and controlled off-orbit re-entry. Large-scale spacecrafts without controlled reentry are usually out of control due to communication faults and finally naturally derail and reenter the atmosphere. Historically, large spacecraft have failed and then entered the atmosphere in many cases, such as the "Skylab" space station in the united states. After the service life of the large-scale spacecraft is ended, the orbit height is gradually reduced to finally enter the atmosphere under the influence of perturbation factors such as atmospheric resistance, gravitational attraction and the like in a low-earth orbit environment, the flight attitude, the trajectory and the falling area of the spacecraft are not controllable under the condition, and the risk of accidents caused by disintegration fragments to a population dense area must be avoided by means of forecasting in advance. The method can accurately and reliably track and forecast the orbit evolution of the spacecraft in the orbit degradation process, can obviously reduce the harmfulness of a reentry disintegration zone, improves the reentry forecast precision, and has very important significance for developing the orbit forecast, collision early warning and on-orbit service technology of space non-cooperative targets.
Generally, the whole natural meteority process of the spacecraft is divided into an orbit attenuation process, a reentry disintegration process and a disintegration debris fragment distribution process, and demarcation points of all stages can be respectively assumed to be an orbit descent reentry point and a spacecraft disintegration process. In order to reliably forecast the time, position and speed information of the spacecraft in the derailing reentry, a countermeasure must be made in advance and feasible derailing input parameters are provided for the disintegration process and the calculation of the debris scattering, which needs a reliable orbit forecasting model. The internationally common method for medium-and-long-term forecasting of low-orbit spacecraft is to forecast by combining two-line root (TLE) and a simplified root perturbation model. The TLE (Two-Line Element) data format is a Two-Line orbit number which is internationally and generally used at present and used for describing satellite orbit parameters, the conventional TLE data comprises parameters such as ephemeris of a satellite, a flat orbit number and the like, and does not comprise orbit error information, and errors of the TLE data are influenced by factors such as observation equipment, an orbit model, observation data, space environment conditions and the like. Because the number of TLE tracks is the average number of tracks, if the instantaneous number of tracks is to be obtained, extrapolation is carried out according to a set track forecasting model, namely, the track forecasting models such as SGP4 or SGP8, and corresponding disturbance terms are added, and the SGP4 track forecasting model is selected as a reference model to carry out track number extrapolation calculation. The SGP4 model was developed in 1966 based on Kozai's analytic theory. The error source of the prediction accuracy of the TLE + SGP4 model mainly comprises two aspects: the method is characterized in that TLE root errors exist, and an index atmospheric density model with atmospheric density changing along with altitude is adopted by an SGP4 model. The main unknown item in the TLE number is the aerodynamic resistance item, and the conventional method is to use historical data to solve the TLE number and other orbit numbers together to obtain an averaged calculation initial value, and the prediction process is an invariant item, which may cause large error accumulation. The resistance coefficient items obtained by simply depending on historical data fitting still cannot meet the forecasting precision.
The atmospheric resistance has obvious influence on the orbit prediction of the low-orbit spacecraft, 5% of resistance error is considered, and the maximum error of the 24-hour position calculation precision of the low-orbit space target can reach 64 km. For the process that a large-scale spacecraft returns from a space orbit and reenters, the large-scale spacecraft can experience the stages of in-orbit free molecular flow, thin transition flow, high-altitude near-continuous slip flow and the like in sequence, and the process is a cross-watershed multi-scale non-equilibrium change process. How to accurately simulate the problem of cross-basin complex multi-physical-field unbalanced streaming in the process of the large spacecraft returning to the reentry in an uncontrolled manner at extremely high speed plays a crucial role in the on-orbit flight and reentry forecasting precision of the large spacecraft.
In the existing internationally adopted spacecraft orbit forecasting methods, a method of fitting a ballistic coefficient is adopted for aerodynamic consideration, even a constant resistance coefficient mode is directly adopted aerodynamically, the influence of the near-earth orbit cross-basin aerodynamic characteristics on the aerodynamic force and moment coefficient of a large spacecraft under different orbit heights and attitudes is not considered, and the influence on the aerodynamic characteristic change and the perturbation of the medium and long-term orbit caused by the attitude change of the large uncontrolled spacecraft is not fully considered. However, the flight and reentry process of the actual large uncontrolled spacecraft is subjected to complex cross-basin aerodynamic perturbation under the conditions of different orbit heights and flight attitudes. As a non-conservative force, the energy dissipation of the large-scale spacecraft caused by aerodynamic force in the process of in-orbit flight is difficult to ignore, and even is a main factor influencing the process of converting the elliptical orbit into the near-circular orbit. The existing technical method cannot fully consider the influence of the aerodynamic characteristics and attitude changes of the near-earth orbit cross-basin unbalanced streaming on the perturbation of the orbit of the large spacecraft, so that the orbit prediction precision is insufficient, and misleading influence can be even generated on the decision of flight control.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides a perturbation analysis method for a spacecraft aerodynamic fusion orbit, and the method can improve the calculation precision of on-orbit flight and reentry forecast.
In order to achieve the above object, the present invention provides a perturbation analysis method for a spacecraft aerodynamic fusion orbit, which comprises:
step 1: generating first two-row root orbit root information of a spacecraft based on ephemeris data and instantaneous orbit root information in observation data corresponding to the initial moment of the spacecraft;
step 2: calculating to obtain the initial orbit height and the initial velocity vector of the spacecraft based on the state variable of the first two-row root orbit root information of the spacecraft; establishing a coordinate rotation matrix for describing the attitude change of the spacecraft based on the initial attitude observation data of the spacecraft, and converting the initial attitude observation data according to the coordinate rotation matrix to obtain an initial attack angle and an initial sideslip angle;
and step 3: calculating and obtaining an initial aerodynamic drag coefficient and an initial aerodynamic disturbance moment coefficient of the spacecraft based on the initial orbit altitude, the initial velocity vector, the initial angle of attack and the initial sideslip angle;
and 4, step 4: calculating to obtain a ballistic coefficient based on the initial aerodynamic resistance coefficient and the surface-to-mass ratio of the spacecraft, and replacing a resistance item in the first two-row radical orbit radical information by using the ballistic coefficient to obtain second two-row radical orbit radical information;
and 5: performing track root information extrapolation calculation based on the second two-line root track root information and the track forecasting model to obtain the track height and the velocity vector of the spacecraft corresponding to the next moment; if the orbit altitude obtained in the step 5 reaches the reentry orbit altitude, outputting ephemeris data and orbit information of the spacecraft corresponding to the next moment; if the height of the track obtained in the step 5 does not reach the height of the reentry track, executing a step 6;
step 6: calculating attitude angle data corresponding to the spacecraft at the next moment by using a state equation of the spacecraft based on the initial attitude observation data and the initial aerodynamic disturbance moment coefficient, and converting the attitude angle data according to the coordinate rotation matrix to obtain an attack angle and a sideslip angle corresponding to the next moment;
and 7: and substituting the orbit altitude, the velocity vector and the attack angle and the sideslip angle obtained in the step 5 and obtained in the step 6 into the step 3 to calculate and obtain the aerodynamic resistance coefficient and the aerodynamic disturbance moment coefficient of the spacecraft corresponding to the next moment, and continuing to iteratively execute the steps 4 to 7 based on the aerodynamic resistance coefficient and the aerodynamic disturbance moment coefficient of the spacecraft corresponding to the next moment until the orbit altitude obtained in the step 5 reaches the re-entering orbit altitude.
The orbit attenuation simulation of the large uncontrolled spacecraft in the low-earth orbit environment needs to fully consider the influence of cross-basin aerodynamic characteristics on the orbit attenuation process, meanwhile, most uncontrolled reentry spacecrafts are in an attitude out-of-control state, and the influence of the rotating attitude of the spacecraft needs to be simultaneously considered when the orbit dynamics characteristics of the large uncontrolled spacecraft in the low-orbit or even ultra-low orbit condition at the end of the service life are researched. Therefore, a higher-precision pneumatic coefficient modeling method needs to be integrated into an orbit dynamics model by combining the attitude and the cross-basin pneumatic environment, so that the on-orbit flight and reentry forecast precision is improved. The method considers the influence of the near-earth orbit cross-basin aerodynamic characteristics on the aerodynamic force and moment coefficient of the large spacecraft under different orbit heights and attitudes, and fully considers the influence on the aerodynamic characteristic change caused by the attitude change of the large uncontrolled spacecraft and the perturbation of the medium and long-term orbit, so that the method considers the attitude rotation effect of the large uncontrolled spacecraft, establishes a pneumatic localization fast engineering algorithm based on the cross-basin aerodynamic characteristics, establishes a pneumatic fusion orbit perturbation analysis method by combining with an average orbit root number model, and can improve the calculation precision of on-orbit flight and reentry prediction.
Preferably, in the method, the step 1 specifically includes: acquiring observation data corresponding to the initial time of the spacecraft in a J2000 geocentric inertial coordinate system, and generating the first two-row root orbit root information in a TEME instant true equator spring deciduous point coordinate system based on ephemeris data and instant orbit root information in the observation data.
Preferably, in the method, the orbit prediction model is an SGP4 model.
Preferably, in the method, the state variable of the two-row root track root information is represented as:
Figure 100002_DEST_PATH_IMAGE001
wherein the content of the first and second substances,
Figure 100002_DEST_PATH_IMAGE003
as a matter of time, the time is,
Figure 100002_DEST_PATH_IMAGE005
the average number of tracks in the form of a two-row number,
Figure 974137DEST_PATH_IMAGE006
is time of day
Figure 100002_DEST_PATH_IMAGE007
Including position and velocity,
Figure 51814DEST_PATH_IMAGE008
for SGP4 model function, initial time
Figure 100002_DEST_PATH_IMAGE009
Corresponding variable
Figure 192683DEST_PATH_IMAGE010
The average number of tracks in the form of a two-row number,
Figure 100002_DEST_PATH_IMAGE011
including eccentricity
Figure 141048DEST_PATH_IMAGE012
Inclination of the track
Figure 100002_DEST_PATH_IMAGE013
The right ascension channel
Figure 577845DEST_PATH_IMAGE014
Argument of near place
Figure 100002_DEST_PATH_IMAGE015
Flat near point angle
Figure 142819DEST_PATH_IMAGE016
And average angular velocity
Figure 100002_DEST_PATH_IMAGE017
；
Resistance term in two rows of radicals
Figure 995368DEST_PATH_IMAGE018
For normalized atmospheric drag coefficients, the form is calculated as:
Figure 100002_DEST_PATH_IMAGE019
wherein the coefficient of ballistic trajectory
Figure 562354DEST_PATH_IMAGE020
The calculation form is:
Figure 100002_DEST_PATH_IMAGE021
wherein the content of the first and second substances,
Figure 170052DEST_PATH_IMAGE022
is a dimensionless aerodynamic drag coefficient,
Figure 100002_DEST_PATH_IMAGE023
for space flightThe surface-to-quality ratio of the device;
Figure 956743DEST_PATH_IMAGE024
in order to be the reference value of the atmospheric density,
Figure 100002_DEST_PATH_IMAGE025
，
Figure 675300DEST_PATH_IMAGE026
the radius of the earth.
Preferably, in the method, the aerodynamic drag coefficient
Figure 394994DEST_PATH_IMAGE022
The calculation method is as follows:
Figure 100002_DEST_PATH_IMAGE027
wherein the content of the first and second substances,
Figure 173595DEST_PATH_IMAGE028
is the aerodynamic force applied to the spacecraft,
Figure 100002_DEST_PATH_IMAGE029
is the average atmospheric density at the altitude where the spacecraft is located,
Figure 713160DEST_PATH_IMAGE030
is the characteristic area of the spacecraft,
Figure 100002_DEST_PATH_IMAGE031
is the velocity of the spacecraft relative to the atmosphere.
Preferably, a coordinate rotation matrix describing the change of the spacecraft attitude is established based on the initial attitude observation data of the spacecraft, and the initial attitude observation data is converted according to the coordinate rotation matrix to obtain an initial attack angle and an initial sideslip angle, which specifically includes:
step a: calculating to obtain initial attitude angle data of the spacecraft based on the initial attitude observation data of the spacecraft;
step b: describing the initial attitude angle data in a quaternion form, and establishing a coordinate rotation matrix for describing attitude change;
step c: and converting the initial attitude angle data of the spacecraft into an airflow coordinate system to convert the initial attitude angle data and the coordinate rotation matrix to obtain an initial attack angle and an initial sideslip angle.
Preferably, in order to avoid the singularity of the solution, the method describes the attitude change of the spacecraft in a quaternion mode, and the state equation of the spacecraft in the step 6 is obtained by combining an attitude kinematics equation of the spacecraft and an attitude dynamics equation of the spacecraft;
describing the attitude change of the spacecraft in a quaternion mode, wherein the attitude kinematic equation of the spacecraft described by the quaternion is obtained by the following steps:
Figure 474224DEST_PATH_IMAGE032
wherein the content of the first and second substances,
Figure 100002_DEST_PATH_IMAGE033
is the differential quantity of the attitude quaternion,
Figure 782845DEST_PATH_IMAGE034
is the quaternion of the attitude,
Figure 100002_DEST_PATH_IMAGE035
for the cross-product sign of the vector matrix,
Figure 201188DEST_PATH_IMAGE036
representing the angular velocity of the extended projectile coordinate system relative to the orbital coordinate system and the inertial angular velocity of the projectile
Figure 100002_DEST_PATH_IMAGE037
The relationship of (1) is:
Figure 962471DEST_PATH_IMAGE038
wherein the content of the first and second substances,
Figure 100002_DEST_PATH_IMAGE039
the rotating angular velocity of the spacecraft projectile system,
Figure 22831DEST_PATH_IMAGE040
the symbols are transposed for the matrix,
Figure 100002_DEST_PATH_IMAGE041
a pose matrix of the projectile coordinate system to the orbit coordinate system described by a quaternion,
Figure 185959DEST_PATH_IMAGE042
for the track angular velocity, a coordinate rotation matrix for describing the attitude change in a quaternion form can be obtained by combining the specified coordinate axis rotation sequence;
the attitude dynamics equation of the spacecraft is as follows:
Figure 100002_DEST_PATH_IMAGE043
wherein the content of the first and second substances,
Figure 804897DEST_PATH_IMAGE044
is a matrix of the moment of inertia of the spacecraft,
Figure 100002_DEST_PATH_IMAGE045
is the angular velocity vector of the rotation of the spacecraft,
Figure 53476DEST_PATH_IMAGE046
is the differential quantity of the rotation angular velocity vector of the spacecraft,
Figure 100002_DEST_PATH_IMAGE047
in order to be a gravity gradient moment,
Figure 917526DEST_PATH_IMAGE048
is a pneumatic disturbance moment;
Figure 100002_DEST_PATH_IMAGE049
is marked as
Figure 669582DEST_PATH_IMAGE050
，
Figure 492044DEST_PATH_IMAGE050
Is gyro moment;
the state equation of the spacecraft is expressed as:
Figure 100002_DEST_PATH_IMAGE051
wherein
Figure 227919DEST_PATH_IMAGE052
In order to be a state variable, the state variable,
Figure 100002_DEST_PATH_IMAGE053
is a generalized matrix of quaternion differential quantities,
Figure 630082DEST_PATH_IMAGE054
，
Figure 100002_DEST_PATH_IMAGE055
is the angular velocity of rotation of the spacecraft about the x-axis,
Figure 758DEST_PATH_IMAGE056
is the differential quantity of the rotation angular velocity of the spacecraft around the x-axis,
Figure 100002_DEST_PATH_IMAGE057
is the angular velocity of the spacecraft about the y-axis,
Figure 462963DEST_PATH_IMAGE058
is the differential quantity of the rotation angular velocity of the spacecraft around the y-axis,
Figure 100002_DEST_PATH_IMAGE059
is the angular velocity of the spacecraft about the z-axis,
Figure 686134DEST_PATH_IMAGE060
is a differential quantity of the rotation angular velocity of the spacecraft around the z-axis,
Figure 100002_DEST_PATH_IMAGE061
、
Figure 626408DEST_PATH_IMAGE062
、
Figure 100002_DEST_PATH_IMAGE063
and
Figure 618635DEST_PATH_IMAGE064
the differential quantities are respectively corresponding to four elements of a quaternion;
and integrating the state variable by using the state equation of the spacecraft to obtain the attitude angle data of the spacecraft at the required moment.
Preferably, the aerodynamic disturbance torque is used in the method
Figure 782900DEST_PATH_IMAGE066
The method comprises the following steps:
Figure 227788DEST_PATH_IMAGE068
、
Figure 502912DEST_PATH_IMAGE070
and
Figure 880803DEST_PATH_IMAGE072
，
Figure 448926DEST_PATH_IMAGE068
、
Figure 912268DEST_PATH_IMAGE070
and
Figure 459924DEST_PATH_IMAGE072
the calculation method is as follows:
Figure 100002_DEST_PATH_IMAGE073
wherein the roll moment coefficient is
Figure 100002_DEST_PATH_IMAGE075
Coefficient of pitching moment of
Figure 100002_DEST_PATH_IMAGE077
Yaw moment coefficient of
Figure 100002_DEST_PATH_IMAGE079
The characteristic length of the spacecraft is
Figure 100002_DEST_PATH_IMAGE081
，
Figure 98847DEST_PATH_IMAGE082
In order to be a gravity gradient moment,
Figure 100002_DEST_PATH_IMAGE083
for aerodynamic disturbance moments on the x-axis,
Figure 73756DEST_PATH_IMAGE084
for the aerodynamic disturbance moment on the y-axis,
Figure 100002_DEST_PATH_IMAGE085
for aerodynamic disturbance moments in the z-axis,
Figure 493236DEST_PATH_IMAGE086
is the average atmospheric density at the altitude where the spacecraft is located,
Figure 602838DEST_PATH_IMAGE030
is the characteristic area of the spacecraft,
Figure 100002_DEST_PATH_IMAGE087
is the velocity of the spacecraft relative to the atmosphere.
Preferably, in the method, the aerodynamic disturbance torque
Figure 158584DEST_PATH_IMAGE088
The calculation method is as follows:
Figure 100002_DEST_PATH_IMAGE089
wherein the content of the first and second substances,
Figure 569974DEST_PATH_IMAGE090
is the radial diameter of the aerodynamic pressure center relative to the mass center of the spacecraft,
Figure 100002_DEST_PATH_IMAGE091
is the aerodynamic force that the spacecraft is subjected to.
One or more technical schemes provided by the invention at least have the following technical effects or advantages:
the invention provides an orbit perturbation analysis method of a large-scale uncontrolled spacecraft, which comprehensively considers the influences of fine cross-basin aerodynamic modeling and attitude rotation effect, and forecasts the degradation orbit information of the spacecraft in advance, overcomes the defect that the traditional technical method does not take pneumatic influence factors into consideration, can remarkably improve the forecasting precision of the in-orbit flight and reentry process of the large-scale uncontrolled spacecraft, can forecast in advance, provides a basis for forecasting of information such as reentry time and position and the like, and is favorable for developing problem solving forecasting, hazard analysis and risk assessment of the reentry process of the large-scale uncontrolled spacecraft. Meanwhile, the technology can support multidisciplinary coupling research and optimization design of the problems of uncontrolled reentry of a large spacecraft, maneuvering orbital transfer and the like.
Drawings
The accompanying drawings, which are included to provide a further understanding of the embodiments of the invention and are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description serve to explain the principles of the invention;
FIG. 1 is a schematic flow chart of a perturbation analysis method for a spacecraft aerodynamic fusion orbit in the invention.
Detailed Description
In order that the above objects, features and advantages of the present invention can be more clearly understood, a more particular description of the invention will be rendered by reference to the appended drawings. It should be noted that the embodiments of the present invention and features of the embodiments may be combined with each other without conflicting with each other.
In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present invention, however, the present invention may be practiced in other ways than those specifically described and thus the scope of the present invention is not limited by the specific embodiments disclosed below.
Example one
Referring to fig. 1, fig. 1 is a schematic flow chart of a method for analyzing perturbation of a spacecraft aerodynamic fusion orbit, and an embodiment of the present invention provides a method for analyzing perturbation of a spacecraft aerodynamic fusion orbit, where the method includes:
step 1: generating first two-row root orbit root information of a spacecraft based on ephemeris data and instantaneous orbit root information in observation data corresponding to the initial moment of the spacecraft;
step 2: calculating to obtain the initial orbit height and the initial velocity vector of the spacecraft based on the state variable of the first two-row root orbit root information of the spacecraft; establishing a coordinate rotation matrix for describing the attitude change of the spacecraft based on the initial attitude observation data of the spacecraft, and converting the initial attitude observation data according to the coordinate rotation matrix to obtain an initial attack angle and an initial sideslip angle;
and step 3: calculating and obtaining an initial aerodynamic drag coefficient and an initial aerodynamic disturbance moment coefficient of the spacecraft based on the initial orbit altitude, the initial velocity vector, the initial angle of attack and the initial sideslip angle;
and 4, step 4: calculating to obtain a ballistic coefficient based on the initial aerodynamic resistance coefficient and the surface-to-mass ratio of the spacecraft, and replacing a resistance item in the first two-row radical orbit radical information by using the ballistic coefficient to obtain second two-row radical orbit radical information;
and 5: performing track root information extrapolation calculation based on the second two-line root track root information and the track forecasting model to obtain the track height and the velocity vector of the spacecraft corresponding to the next moment; if the orbit altitude obtained in the step 5 reaches the reentry orbit altitude, outputting ephemeris data and orbit information of the spacecraft corresponding to the next moment; if the height of the track obtained in the step 5 does not reach the height of the reentry track, executing a step 6;
step 6: calculating attitude angle data corresponding to the spacecraft at the next moment by using a state equation of the spacecraft based on the initial attitude observation data and the initial aerodynamic disturbance moment coefficient, and converting the attitude angle data according to the coordinate rotation matrix to obtain an attack angle and a sideslip angle corresponding to the next moment;
and 7: and substituting the orbit altitude, the velocity vector and the attack angle and the sideslip angle obtained in the step 5 and obtained in the step 6 into the step 3 to calculate and obtain the aerodynamic resistance coefficient and the aerodynamic disturbance moment coefficient of the spacecraft corresponding to the next moment, and continuing to iteratively execute the steps 4 to 7 based on the aerodynamic resistance coefficient and the aerodynamic disturbance moment coefficient of the spacecraft corresponding to the next moment until the orbit altitude obtained in the step 5 reaches the re-entering orbit altitude.
In the first embodiment of the present invention, the method specifically includes the following steps:
and generating TLE two-row root orbit information under a TEME instantaneous true equator vernal equinox coordinate system by using ephemeris and instantaneous orbit root information contained in the initial-moment observation data of the large-scale uncontrolled spacecraft under the existing J2000 geocentric inertial coordinate system.
The method comprises the steps of calculating to obtain an initial orbit height and a velocity vector by utilizing initial instantaneous orbit root information of the large-scale uncontrolled spacecraft, obtaining an initial attack angle and a sideslip angle by utilizing initial attitude observation data of the large-scale uncontrolled spacecraft, calling a pneumatic localization fast engineering algorithm, and calculating to obtain an initial pneumatic resistance coefficient of the large-scale uncontrolled spacecraft and a pneumatic moment coefficient for solving an attitude dynamic equation.
Calculating a ballistic coefficient by utilizing the initial aerodynamic resistance coefficient and the surface-to-mass ratio of the large uncontrolled spacecraft, and replacing the resistance item in two lines of TLE (total line of space) with the calculated ballistic coefficient
Figure 476750DEST_PATH_IMAGE018
。
And performing orbit extrapolation calculation by combining two lines of TLE elements containing updated resistance terms with an SGP4 model to obtain the position and speed information of the large-scale uncontrolled spacecraft at the next moment, calling a pneumatic localization fast engineering algorithm and an attitude dynamics calculation program, and integrating to obtain the pneumatic resistance coefficient, the pneumatic moment coefficient and the attitude angle information at the next moment.
And converting the orbit height according to the position and speed information at the next moment, judging whether the given reentry orbit height is reached, if so, terminating the calculation, and outputting the ephemeris and the orbit information at the moment.
If the specified re-entry height is not reached, the pneumatic resistance coefficient at the next moment is also used to update the TLE two-line number of resistance terms
Figure 631788DEST_PATH_IMAGE018
And returning to the step of calculating the orbit extrapolation, and sequentially and iteratively calculating until the given height of the reentry orbit is reached, terminating the calculation and outputting the ephemeris and the orbit information at the termination moment.
The method comprises the following specific technical scheme:
firstly, the instantaneous orbit root and ephemeris information obtained by observation data under a J2000 geocentric inertial coordinate system need to be converted into a TLE two-row root model. The SGP4 model adopts a true equator vernal equinox coordinate system, and the state variable of two-row root of TLE can be expressed as
Figure 573199DEST_PATH_IMAGE001
(1)
Wherein the content of the first and second substances,
Figure 889911DEST_PATH_IMAGE006
is time of day
Figure 549563DEST_PATH_IMAGE007
The state vector (position and velocity) of (c),
Figure 508291DEST_PATH_IMAGE008
for SGP4 model function, corresponding to initial time
Figure 569788DEST_PATH_IMAGE009
Variable of (2)
Figure 87095DEST_PATH_IMAGE010
Then the average number of tracks in TLE form
Figure 968463DEST_PATH_IMAGE011
Including eccentricity
Figure 465304DEST_PATH_IMAGE012
Inclination of the track
Figure 381307DEST_PATH_IMAGE013
The right ascension channel
Figure 305401DEST_PATH_IMAGE014
Argument of near place
Figure 939644DEST_PATH_IMAGE015
Flat near point angle
Figure 240176DEST_PATH_IMAGE016
And average angular velocity
Figure 10686DEST_PATH_IMAGE017
In particular the resistance term in the two-line root of TLE
Figure 371260DEST_PATH_IMAGE018
For normalized atmospheric drag coefficients, the form is calculated as follows:
Figure 492800DEST_PATH_IMAGE019
(2)
wherein the coefficient of ballistic trajectory
Figure 65863DEST_PATH_IMAGE020
The calculation form is as follows:
Figure 956459DEST_PATH_IMAGE021
(3)
wherein the content of the first and second substances,
Figure 487934DEST_PATH_IMAGE022
is a dimensionless aerodynamic drag coefficient,
Figure 532989DEST_PATH_IMAGE023
the surface-to-mass ratio of the large uncontrolled spacecraft is obtained;
Figure 175323DEST_PATH_IMAGE024
in order to be the reference value of the atmospheric density,
Figure 920425DEST_PATH_IMAGE092
，
Figure 622801DEST_PATH_IMAGE026
the radius of the earth. Solving equation (1) can obtain the position and speed information at any moment.
And then, a localized fast engineering algorithm of aerodynamic force and moment coefficients of the large uncontrolled spacecraft needs to be established based on the cross-basin aerodynamic characteristics. Because the aerodynamic coefficient of the large uncontrolled spacecraft obviously changes along with the height and the attitude of the orbit, the high efficiency of the calculation of the aerodynamic characteristics is also considered according to the quick response of the return and reentry processes of the spacecraft in the orbit forecasting process.
The large uncontrolled spacecraft orbit flight process is a very high supersonic flow problem of a high rarefied free molecular flow state of hundreds to tens of orders of magnitude Knudsen and a complex configuration of a multi-physics field, and the cross-basin gas flow problem experienced in the low orbit flight process can be solved by adopting a unified algorithm. The free molecular flow theory is adopted to control the highly thin flow region, and Nocilla wall reflection models corrected based on different model materials can be adopted to calculate the pressure coefficient and the friction coefficient. Under the assumption of a Maxwell equilibrium state gas molecule velocity distribution function, the pressure and friction coefficient of each surface element can be calculated, and verification analysis correctness is confirmed. In the continuous flow region where Knudsen numbers tend to be very small, it can be based on modified Newton's theory of non-viscous flow. However, aiming at the long orbit decay of the large uncontrolled spacecraft, the large uncontrolled spacecraft undergoes the unbalanced reentry process of multiple physical fields and multiple scales in different flight attitudes, in order to reliably capture the influence of aerodynamic characteristics on different arc-segment orbits, the method specifically adopts a literature that the aerodynamic characteristics of the cross-basin are summarized and described in the general literature of integrated modeling and calculation research (manned space) of the low-orbit control aerodynamic characteristics of the spacecraft, and the method establishes the pneumatic localization fast engineering algorithm based on the cross-basin aerodynamic characteristics without repeated description, and calculating typical streaming states of the spacecraft in the height areas (340 km,280km, 200km and 120 km) in the process of descending the spacecraft from the orbit to the reentry respectively at intervals of 10km, 5km and 2.5km by using a cross-streaming area aerodynamic numerical method, and taking the result as the boundary value of each sub-area. And for the middle transition zone connecting the boundaries of each sub-zone, a modified Boettcher/Legge asymmetric bridge function theory is adopted to develop a bridge function related to asymmetric pressure and friction coefficient which can be described in a segmented manner, so that an integrated rapid calculation technology for aerodynamic characteristics of the low-orbit cross-flow zone is formed, and a specific construction method can refer to the document.
The pressure coefficient calculated by using the modified Newtonian non-viscous flow theory for the continuous flow region is as follows:
Figure DEST_PATH_IMAGE093
（4）
wherein
Figure 187775DEST_PATH_IMAGE094
、
Figure DEST_PATH_IMAGE095
And
Figure 837062DEST_PATH_IMAGE096
the pressure, density and speed of the free incoming flow,
Figure DEST_PATH_IMAGE097
is an angle of an object plane, and is a plane,
Figure 436671DEST_PATH_IMAGE098
at maximum surface pressure
Figure DEST_PATH_IMAGE099
The pressure coefficient of the post-shock stationary point is specifically calculated as follows:
Figure 542905DEST_PATH_IMAGE100
(5)
wherein
Figure DEST_PATH_IMAGE101
In order to obtain a mach number of the incoming flow,
Figure 329595DEST_PATH_IMAGE102
is the specific heat ratio of the gas. The free molecular flow theory is adopted to control the highly thin flow region, a Nocilla wall reflection model corrected based on different model materials can be adopted to calculate the pressure coefficient and the friction coefficient, and the pressure of each surface element is calculated under the assumption of Maxwell equilibrium state gas distribution
Figure DEST_PATH_IMAGE103
And coefficient of friction
Figure 48152DEST_PATH_IMAGE104
Comprises the following steps:
Figure DEST_PATH_IMAGE105
(6)
Figure 502268DEST_PATH_IMAGE106
(7)
wherein the content of the first and second substances,
Figure DEST_PATH_IMAGE107
in order to be an exponential function of the,
Figure 15288DEST_PATH_IMAGE108
、
Figure DEST_PATH_IMAGE109
respectively the object plane and the incoming flow temperature,
Figure 554854DEST_PATH_IMAGE110
and
Figure DEST_PATH_IMAGE111
normal and tangential object plane momentum adaptation coefficients respectively,
Figure 77102DEST_PATH_IMAGE112
in order to be a function of the error,
Figure DEST_PATH_IMAGE113
is the speed ratio, defined as:
Figure 624539DEST_PATH_IMAGE114
(8)
wherein
Figure DEST_PATH_IMAGE115
Is the universal gas constant. And smoothly overlapping the continuous flow coefficient and the free molecular flow coefficient in a transition flow region by adopting a semi-empirical local bridge function method, so that the continuous flow coefficient and the free molecular flow coefficient respectively approach to a flow function. At a given object plane angle
Figure 574041DEST_PATH_IMAGE116
When is coming into contact withThe pressure and friction coefficients on the surface elements are:
Figure DEST_PATH_IMAGE117
(9)
Figure 335323DEST_PATH_IMAGE118
(10)
wherein
Figure DEST_PATH_IMAGE119
And
Figure 661263DEST_PATH_IMAGE120
pressure and friction bridge functions, respectively, dependent on independent parameters
Figure DEST_PATH_IMAGE121
、
Figure 558811DEST_PATH_IMAGE122
And
Figure DEST_PATH_IMAGE123
，
Figure 679214DEST_PATH_IMAGE121
knudsen number. The irregular object of the large uncontrolled spacecraft can be represented by a triangular unstructured grid surface element, the large uncontrolled spacecraft is subjected to object forming processing by adopting a general approximate surface element method, when the surface elements are divided into small enough, the normal force and tangential force integrals on all the surface elements are integrated along the whole object surface of the large uncontrolled spacecraft geometry, and therefore the error caused by the overall aerodynamic force coefficient is obtained to be smaller. And finally, a high-performance computer is fully utilized to develop DSMC and solve a Boltzmann model equation unified algorithm to carry out fine calculation on aerodynamic characteristics, so that a localized fast engineering algorithm suitable for the aerodynamic characteristics of the large uncontrolled spacecraft in a specific geometric configuration is determined by correction, and the corresponding aerodynamic coefficients of the large uncontrolled spacecraft under the conditions of different orbit heights and attack angles are obtained.
In addition, the attitude information of the large-scale uncontrolled spacecraft, including the angle of attack and the sideslip angle, needs to be synchronously calculated in real time by using the attitude dynamics model. The attack angle and the sideslip angle are included angles of a large-scale uncontrolled spacecraft projectile coordinate system relative to a speed coordinate system, and the attack angle and the sideslip angle can be obtained only by solving attitude angle changes of the large-scale uncontrolled spacecraft projectile coordinate system relative to an orbit coordinate system. The attitude angle is generally defined as the euler angle of the large uncontrolled spacecraft in rotation around the x, y and z axes, respectively:
Figure 458951DEST_PATH_IMAGE124
-a roll angle of the roll-off,
Figure DEST_PATH_IMAGE125
-a pitch angle,
Figure 555958DEST_PATH_IMAGE126
-yaw angle.
In order to avoid the singularity of the solution, the attitude change of the large-scale uncontrolled spacecraft is described in a quaternion mode. The quaternion representation is directly related to the euler angle representation. The following four elements are defined:
Figure DEST_PATH_IMAGE127
（11）
wherein q is0、q1、q2、q3Four elements of a quaternion, q is a vector of three imaginary elements in the quaternion,
Figure 573593DEST_PATH_IMAGE128
the symbols are transposed for the matrix, n is the rotation axis, the first of the four elements in the quaternion represents the euler rotation angle, the last three represent the direction of the euler rotation axis. The four elements are combined in the following vector form:
Figure DEST_PATH_IMAGE129
（12）
it is clear that these four elements satisfy the normalization condition. According to the Euler rotation formula and the definition of quaternion, the attitude matrix represented by quaternion is:
Figure 130476DEST_PATH_IMAGE130
（13）
wherein the content of the first and second substances,
Figure DEST_PATH_IMAGE131
is a rotational inertia matrix of a large uncontrolled spacecraft,
Figure 600771DEST_PATH_IMAGE132
an augmentation matrix being a quaternion matrix q;
in order to establish the relationship between the instantaneous change rate of the attitude parameters and the rotation angular velocity, it is first required to solve the attitude kinematics equation. If the rotating speed of the attitude of the large uncontrolled spacecraft relative to the reference coordinate is
Figure DEST_PATH_IMAGE133
The rotating shaft is
Figure 2934DEST_PATH_IMAGE134
I.e. by
Figure DEST_PATH_IMAGE135
If at the moment
Figure 875075DEST_PATH_IMAGE136
Attitude matrix is
Figure DEST_PATH_IMAGE137
In a
Figure 602860DEST_PATH_IMAGE138
The time attitude matrix is
Figure DEST_PATH_IMAGE139
And then:
Figure 793407DEST_PATH_IMAGE140
（14）
wherein
Figure 530419DEST_PATH_IMAGE141
To around a rotating shaft
Figure 788225DEST_PATH_IMAGE134
Is rotated over
Figure DEST_PATH_IMAGE142
The angle rotation matrix applies the relation between the attitude quaternion and the direction cosine to obtain the change equation of the attitude quaternion, which is as follows:
Figure 421332DEST_PATH_IMAGE143
（15）
wherein the content of the first and second substances,
Figure 397378DEST_PATH_IMAGE056
is the angular velocity of rotation of the spacecraft about the x-axis,
Figure 406922DEST_PATH_IMAGE058
is the angular velocity of the spacecraft about the y-axis,
Figure 519235DEST_PATH_IMAGE060
is the angular velocity of the spacecraft about the z-axis,
Figure DEST_PATH_IMAGE144
、
Figure 588822DEST_PATH_IMAGE145
、
Figure 786585DEST_PATH_IMAGE063
and
Figure 865400DEST_PATH_IMAGE064
the differential quantities are corresponding to four elements of a quaternion.
The attitude kinematics equation resulting in the quaternion description is as follows:
Figure 832219DEST_PATH_IMAGE032
（16）
wherein the content of the first and second substances,
Figure 40084DEST_PATH_IMAGE033
is the differential quantity of the attitude quaternion,
Figure 725143DEST_PATH_IMAGE034
is the quaternion of the attitude,
Figure DEST_PATH_IMAGE146
for the cross-product sign of the vector matrix,
Figure 76490DEST_PATH_IMAGE036
representing the angular velocity of the extended projectile coordinate system relative to the orbital coordinate system and the inertial angular velocity of the projectile
Figure 897815DEST_PATH_IMAGE037
The relationship of (a) to (b) is as follows:
Figure 309205DEST_PATH_IMAGE038
（17）
in the formula
Figure 404198DEST_PATH_IMAGE041
A pose matrix of the projectile coordinate system to the orbit coordinate system described by a quaternion,
Figure 824816DEST_PATH_IMAGE042
is the track angular velocity.
On the basis of the attitude kinematics equation, an attitude kinematics equation can be further established for solving the attitude angle in an integral mode. Considering the attitude angle change of the large-scale uncontrolled spacecraft under the action of the perturbation moments such as the aerodynamic disturbance moment geocentric non-spherical gravitation and the like, an attitude dynamics equation is established as follows:
Figure 766227DEST_PATH_IMAGE147
（18）
wherein
Figure 348518DEST_PATH_IMAGE044
Is a rotational inertia matrix of a large uncontrolled spacecraft,
Figure 8169DEST_PATH_IMAGE049
an item is generally denoted as
Figure 966898DEST_PATH_IMAGE050
By angular velocity
Figure 28395DEST_PATH_IMAGE045
Of the representation
Figure 47166DEST_PATH_IMAGE050
Comprises the following steps:
Figure DEST_PATH_IMAGE148
（19）
the state equation of the attitude motion of the large-scale uncontrolled spacecraft can be expressed by combining an attitude kinematics equation as follows:
Figure 397376DEST_PATH_IMAGE149
（20）
wherein
Figure 159796DEST_PATH_IMAGE052
Is a state variable.
To the right of equation (20)
Figure 75799DEST_PATH_IMAGE048
For the aerodynamic disturbance torque, the aerodynamic disturbance torque model can be expressed as:
Figure 265472DEST_PATH_IMAGE089
（21）
wherein
Figure 899716DEST_PATH_IMAGE090
Is the radius of the mass center of the aerodynamic pressure center relative to the large-scale uncontrolled spacecraft,
Figure 433203DEST_PATH_IMAGE091
the aerodynamic force of a large uncontrolled spacecraft is expressed as follows:
Figure 203713DEST_PATH_IMAGE027
（22）
wherein
Figure 564287DEST_PATH_IMAGE022
In order to obtain the coefficient of aerodynamic drag,
Figure 420248DEST_PATH_IMAGE029
is the average atmospheric density of the height of the large uncontrolled spacecraft,
Figure 524470DEST_PATH_IMAGE030
is a characteristic area of a large uncontrolled spacecraft,
Figure 149486DEST_PATH_IMAGE087
the speed of the large uncontrolled spacecraft relative to the atmosphere.
And (3) calculating the result according to the localization rapid engineering algorithm in the step 2: roll moment coefficient of
Figure 680962DEST_PATH_IMAGE150
Coefficient of pitching moment of
Figure DEST_PATH_IMAGE151
Yaw moment coefficient of
Figure 758639DEST_PATH_IMAGE152
The characteristic length of the large-scale uncontrolled spacecraft is
Figure DEST_PATH_IMAGE153
Then the aerodynamic moment at the three-axis component of the projectile coordinate is:
Figure 135394DEST_PATH_IMAGE073
（23）
to the right of equation (20)
Figure 880496DEST_PATH_IMAGE047
The gravity gradient moment is generated by different gravitations of the earth to parts of the large-scale uncontrolled spacecraft.
Converting the attitude angle obtained by solving the attitude dynamics equation into an airflow coordinate system to obtain an attack angle and a sideslip angle through conversion, converting the height and the velocity vector of the orbit by using the position and the velocity obtained by solving the equation (1), calling the pneumatic localization fast engineering algorithm of the step 2, and obtaining the corresponding pneumatic resistance coefficient under the conditions of corresponding orbit height and attack angle through integral equations (9) and (10)
Figure 317294DEST_PATH_IMAGE022
。
The aerodynamic drag coefficient is obtained by calculation based on a refined aerodynamic localization quick engineering algorithm and by considering the attitude rotation effect
Figure 646382DEST_PATH_IMAGE022
Calculating and updating the ballistic coefficient in equation (3)
Figure 92407DEST_PATH_IMAGE020
And calculating a resistance term in the two-line root of the TLE according to equation (2)
Figure 957594DEST_PATH_IMAGE018
Meanwhile, the position and speed calculation of the next moment is carried out by combining the SGP4 model, and the solution is iterated in sequence until the given reentry is reachedAnd (4) outputting the time, the position and the speed information corresponding to the reentry altitude and the previously experienced orbit ephemeris data to complete the on-orbit prediction and the reentry prediction of the large-scale uncontrolled spacecraft.
The method established by the invention is used for carrying out orbit perturbation analysis and on-orbit flight and reentry forecast on the similar Tiangong I, the given reentry height is 120km, the initial time is 3 months in 2018, 18 months in 2018, 12 days in 18 months, 16 minutes and 1 second, the obtained specific position, speed and attitude information is shown in table 2, and the reentry forecast information is as follows:
reentry time (year, month, day, hour, minute, second):
2018 4 3 19 54 1。
then enter longitude, latitude, altitude (°, km):
-144.76 16.18 119.95。
while preferred embodiments of the present invention have been described, additional variations and modifications in those embodiments may occur to those skilled in the art once they learn of the basic inventive concepts. Therefore, it is intended that the appended claims be interpreted as including preferred embodiments and all such alterations and modifications as fall within the scope of the invention.
It will be apparent to those skilled in the art that various changes and modifications may be made in the present invention without departing from the spirit and scope of the invention. Thus, if such modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include such modifications and variations.

Claims (9)

1. A perturbation analysis method for a spacecraft aerodynamic fusion orbit is characterized by comprising the following steps:
step 1: generating first two-row root orbit root information of a spacecraft based on ephemeris data and instantaneous orbit root information in observation data corresponding to the initial moment of the spacecraft;
step 2: calculating to obtain the initial orbit height and the initial velocity vector of the spacecraft based on the state variable of the first two-row root orbit root information of the spacecraft; establishing a coordinate rotation matrix for describing the attitude change of the spacecraft based on the initial attitude observation data of the spacecraft, and converting the initial attitude observation data according to the coordinate rotation matrix to obtain an initial attack angle and an initial sideslip angle;
and step 3: calculating and obtaining an initial aerodynamic drag coefficient and an initial aerodynamic disturbance moment coefficient of the spacecraft based on the initial orbit altitude, the initial velocity vector, the initial angle of attack and the initial sideslip angle;
and 4, step 4: calculating to obtain a ballistic coefficient based on the initial aerodynamic resistance coefficient and the surface-to-mass ratio of the spacecraft, and replacing a resistance item in the first two-row radical orbit radical information by using the ballistic coefficient to obtain second two-row radical orbit radical information;
and 5: performing track root information extrapolation calculation based on the second two-line root track root information and the track forecasting model to obtain the track height and the velocity vector of the spacecraft corresponding to the next moment; if the orbit altitude obtained in the step 5 reaches the reentry orbit altitude, outputting ephemeris data and orbit information of the spacecraft corresponding to the next moment; if the height of the track obtained in the step 5 does not reach the height of the reentry track, executing a step 6;
step 6: calculating attitude angle data corresponding to the spacecraft at the next moment by using a state equation of the spacecraft based on the initial attitude observation data and the initial aerodynamic disturbance moment coefficient, and converting the attitude angle data according to the coordinate rotation matrix to obtain an attack angle and a sideslip angle corresponding to the next moment;
and 7: and substituting the orbit altitude, the velocity vector and the attack angle and the sideslip angle obtained in the step 5 and obtained in the step 6 into the step 3 to calculate and obtain the aerodynamic resistance coefficient and the aerodynamic disturbance moment coefficient of the spacecraft corresponding to the next moment, and continuing to iteratively execute the steps 4 to 7 based on the aerodynamic resistance coefficient and the aerodynamic disturbance moment coefficient of the spacecraft corresponding to the next moment until the orbit altitude obtained in the step 5 reaches the re-entering orbit altitude.
2. The perturbation analysis method for the spacecraft aerodynamic fusion orbit according to claim 1, wherein the step 1 specifically comprises: acquiring observation data corresponding to the initial time of the spacecraft in a J2000 geocentric inertial coordinate system, and generating the first two-row root orbit root information in a TEME instant true equator spring deciduous point coordinate system based on ephemeris data and instant orbit root information in the observation data.
3. A spacecraft aerodynamic fusion orbit perturbation analysis method according to claim 1, wherein the orbit prediction model is an SGP4 model.
4. A spacecraft aerodynamic fusion orbit perturbation analysis method according to claim 1, wherein the state variables of two rows of root orbit root information are expressed as:
Figure DEST_PATH_IMAGE001
wherein the content of the first and second substances,
Figure DEST_PATH_IMAGE003
as a matter of time, the time is,
Figure DEST_PATH_IMAGE005
the average number of tracks in the form of a two-row number,
Figure 908230DEST_PATH_IMAGE006
is time of day
Figure DEST_PATH_IMAGE007
Including position and velocity,
Figure 231895DEST_PATH_IMAGE008
for SGP4 model function, initial time
Figure DEST_PATH_IMAGE009
Corresponding variable
Figure 683736DEST_PATH_IMAGE010
The average number of tracks in the form of a two-row number,
Figure DEST_PATH_IMAGE011
including eccentricity
Figure 610104DEST_PATH_IMAGE012
Inclination of the track
Figure DEST_PATH_IMAGE013
The right ascension channel
Figure 831001DEST_PATH_IMAGE014
Argument of near place
Figure DEST_PATH_IMAGE015
Flat near point angle
Figure 450201DEST_PATH_IMAGE016
And average angular velocity
Figure DEST_PATH_IMAGE017
；
Resistance term in two rows of radicals
Figure 622294DEST_PATH_IMAGE018
For normalized atmospheric drag coefficients, the form is calculated as:
Figure DEST_PATH_IMAGE019
wherein the coefficient of ballistic trajectory
Figure 352353DEST_PATH_IMAGE020
The calculation form is:
Figure DEST_PATH_IMAGE021
wherein the content of the first and second substances,
Figure 693335DEST_PATH_IMAGE022
is a dimensionless aerodynamic drag coefficient,
Figure DEST_PATH_IMAGE023
is the surface-to-mass ratio of the spacecraft;
Figure 93224DEST_PATH_IMAGE024
in order to be the reference value of the atmospheric density,
Figure DEST_PATH_IMAGE025
，
Figure 378712DEST_PATH_IMAGE026
the radius of the earth.
5. A spacecraft aerodynamic fusion orbit perturbation analysis method according to claim 1, wherein aerodynamic drag coefficient
Figure 522248DEST_PATH_IMAGE022
The calculation method is as follows:
Figure DEST_PATH_IMAGE027
wherein the content of the first and second substances,
Figure 576792DEST_PATH_IMAGE028
is the aerodynamic force applied to the spacecraft,
Figure DEST_PATH_IMAGE029
is the average atmospheric density at the altitude where the spacecraft is located,
Figure 177275DEST_PATH_IMAGE030
is the characteristic area of the spacecraft,
Figure DEST_PATH_IMAGE031
is the velocity of the spacecraft relative to the atmosphere.
6. The method for analyzing perturbation of a spacecraft aerodynamic fusion orbit according to claim 1, wherein a coordinate rotation matrix describing changes of the spacecraft attitude is established based on initial attitude observation data of the spacecraft, and the initial attitude observation data is converted according to the coordinate rotation matrix to obtain an initial attack angle and an initial sideslip angle, and specifically comprises:
step a: calculating to obtain initial attitude angle data of the spacecraft based on the initial attitude observation data of the spacecraft;
step b: describing the initial attitude angle data in a quaternion form, and establishing a coordinate rotation matrix for describing attitude change;
step c: and converting the initial attitude angle data of the spacecraft into an airflow coordinate system to convert the initial attitude angle data and the coordinate rotation matrix to obtain an initial attack angle and an initial sideslip angle.
7. A perturbation analysis method for a spacecraft aerodynamic fusion orbit according to claim 6, wherein the state equation of the spacecraft in the step 6 is obtained by combining an attitude kinematics equation of the spacecraft and an attitude dynamics equation of the spacecraft;
describing the attitude change of the spacecraft in a quaternion mode, wherein the attitude kinematic equation of the spacecraft described by the quaternion is obtained by the following steps:
Figure 950059DEST_PATH_IMAGE032
wherein the content of the first and second substances,
Figure DEST_PATH_IMAGE033
is the differential quantity of the attitude quaternion,
Figure 631707DEST_PATH_IMAGE034
is the quaternion of the attitude,
Figure DEST_PATH_IMAGE035
for the cross-product sign of the vector matrix,
Figure 947282DEST_PATH_IMAGE036
representing the angular velocity of the extended projectile coordinate system relative to the orbital coordinate system and the inertial angular velocity of the projectile
Figure DEST_PATH_IMAGE037
The relationship of (1) is:
Figure 220132DEST_PATH_IMAGE038
wherein the content of the first and second substances,
Figure DEST_PATH_IMAGE039
the rotating angular velocity of the spacecraft projectile system,
Figure 480212DEST_PATH_IMAGE040
the symbols are transposed for the matrix,
Figure DEST_PATH_IMAGE041
a pose matrix of the projectile coordinate system to the orbit coordinate system described by a quaternion,
Figure 469945DEST_PATH_IMAGE042
is the track angular velocity;
the attitude dynamics equation of the spacecraft is as follows:
Figure DEST_PATH_IMAGE043
wherein the content of the first and second substances,
Figure 499081DEST_PATH_IMAGE044
is a matrix of the moment of inertia of the spacecraft,
Figure DEST_PATH_IMAGE045
is the angular velocity vector of the rotation of the spacecraft,
Figure 942832DEST_PATH_IMAGE046
is the differential quantity of the rotation angular velocity vector of the spacecraft,
Figure DEST_PATH_IMAGE047
in order to be a gravity gradient moment,
Figure 424629DEST_PATH_IMAGE048
is a pneumatic disturbance moment;
Figure DEST_PATH_IMAGE049
is marked as
Figure 713659DEST_PATH_IMAGE050
，
Figure 597301DEST_PATH_IMAGE050
Is gyro moment;
the state equation of the spacecraft is expressed as:
Figure DEST_PATH_IMAGE051
wherein
Figure 946374DEST_PATH_IMAGE052
In order to be a state variable, the state variable,
Figure DEST_PATH_IMAGE053
is a generalized matrix of quaternion differential quantities,
Figure 915467DEST_PATH_IMAGE054
，
Figure DEST_PATH_IMAGE055
is the angular velocity of rotation of the spacecraft about the x-axis,
Figure 506723DEST_PATH_IMAGE056
is the differential quantity of the rotation angular velocity of the spacecraft around the x-axis,
Figure DEST_PATH_IMAGE057
is the angular velocity of the spacecraft about the y-axis,
Figure 385817DEST_PATH_IMAGE058
is the differential quantity of the rotation angular velocity of the spacecraft around the y-axis,
Figure DEST_PATH_IMAGE059
is the angular velocity of the spacecraft about the z-axis,
Figure 296005DEST_PATH_IMAGE060
is a differential quantity of the rotation angular velocity of the spacecraft around the z-axis,
Figure DEST_PATH_IMAGE061
、
Figure 627760DEST_PATH_IMAGE062
、
Figure DEST_PATH_IMAGE063
and
Figure 383226DEST_PATH_IMAGE064
the differential quantities are respectively corresponding to four elements of a quaternion;
and integrating the state variable by using the state equation of the spacecraft to obtain the attitude angle data of the spacecraft at the required moment.
8. A spacecraft aerodynamic fusion orbit perturbation analysis method according to claim 7, wherein the aerodynamic disturbance torque
Figure 585669DEST_PATH_IMAGE066
The method comprises the following steps:
Figure 135599DEST_PATH_IMAGE068
、
Figure 79284DEST_PATH_IMAGE070
and
Figure 746763DEST_PATH_IMAGE072
，
Figure 193925DEST_PATH_IMAGE068
、
Figure 649177DEST_PATH_IMAGE070
and
Figure 814580DEST_PATH_IMAGE072
the calculation method is as follows:
Figure DEST_PATH_IMAGE073
wherein the roll moment coefficient is
Figure DEST_PATH_IMAGE075
Coefficient of pitching moment of
Figure DEST_PATH_IMAGE077
Yaw moment coefficient of
Figure DEST_PATH_IMAGE079
The characteristic length of the spacecraft is
Figure DEST_PATH_IMAGE081
，
Figure 990477DEST_PATH_IMAGE082
In order to be a gravity gradient moment,
Figure DEST_PATH_IMAGE083
for aerodynamic disturbance moments on the x-axis,
Figure 698670DEST_PATH_IMAGE084
for the aerodynamic disturbance moment on the y-axis,
Figure DEST_PATH_IMAGE085
for aerodynamic disturbance moments in the z-axis,
Figure 855982DEST_PATH_IMAGE086
is the average atmospheric density at the altitude where the spacecraft is located,
Figure DEST_PATH_IMAGE087
is the characteristic area of the spacecraft,
Figure 413740DEST_PATH_IMAGE088
is the velocity of the spacecraft relative to the atmosphere.
9. A spacecraft aerodynamic fusion orbit perturbation analysis method according to claim 7, wherein aerodynamic disturbance torque
Figure DEST_PATH_IMAGE089
The calculation method is as follows:
Figure 190066DEST_PATH_IMAGE090
wherein the content of the first and second substances,
Figure DEST_PATH_IMAGE091
is the radial diameter of the aerodynamic pressure center relative to the mass center of the spacecraft,
Figure 877399DEST_PATH_IMAGE092
is the aerodynamic force that the spacecraft is subjected to.
CN202110899370.3A 2021-08-06 2021-08-06 Perturbation analysis method for spacecraft aerodynamic fusion orbit Active CN113343369B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110899370.3A CN113343369B (en) 2021-08-06 2021-08-06 Perturbation analysis method for spacecraft aerodynamic fusion orbit

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110899370.3A CN113343369B (en) 2021-08-06 2021-08-06 Perturbation analysis method for spacecraft aerodynamic fusion orbit

Publications (2)

Publication Number Publication Date
CN113343369A CN113343369A (en) 2021-09-03
CN113343369B true CN113343369B (en) 2021-10-01

Family

ID=77480890

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110899370.3A Active CN113343369B (en) 2021-08-06 2021-08-06 Perturbation analysis method for spacecraft aerodynamic fusion orbit

Country Status (1)

Country Link
CN (1) CN113343369B (en)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114219187B (en) * 2022-02-22 2022-05-24 中国人民解放军32035部队 Reentry target forecasting method and device based on combination of double stars and foundation and electronic equipment
CN114330080B (en) * 2022-03-04 2022-05-13 中国空气动力研究与发展中心计算空气动力研究所 Prediction method for aircraft surface-symmetric cross-basin flow field
CN114580224B (en) * 2022-05-09 2022-07-05 中国空气动力研究与发展中心设备设计与测试技术研究所 Distributed pneumatic fusion track coupling attitude perturbation analysis method
CN115238836B (en) * 2022-09-23 2023-04-28 中国空气动力研究与发展中心计算空气动力研究所 Fusion method based on correlation degree of pneumatic data and physical model
CN116384599B (en) * 2023-06-06 2023-08-22 中国空气动力研究与发展中心超高速空气动力研究所 Spacecraft LEO circular orbit attenuation process parameter forecasting method based on energy analysis

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109460055A (en) * 2018-10-30 2019-03-12 中国运载火箭技术研究院 A kind of flying vehicles control capacity judging method, device and electronic equipment
CN111310277A (en) * 2020-03-13 2020-06-19 中国运载火箭技术研究院 Modeling method for pipeline transfer characteristics of atmospheric data sensing system, aircraft and storage medium
CN111680462A (en) * 2020-08-11 2020-09-18 北京控制与电子技术研究所 Guidance method and system based on position change of space target in optical phase plane
CN111680354A (en) * 2020-04-20 2020-09-18 北京航空航天大学 Method for calculating self-intersection point of orbit of near-earth regression orbit satellite subsatellite point and photographing point
CN111874267A (en) * 2020-04-30 2020-11-03 中国人民解放军战略支援部队航天工程大学 Low-orbit satellite off-orbit control method and system based on particle swarm optimization
CN112363410A (en) * 2020-11-13 2021-02-12 浙江大学 Intelligent autonomous control research and verification system for spacecraft
CN112389672A (en) * 2020-11-09 2021-02-23 中国运载火箭技术研究院 Aerospace vehicle transverse course stability control characteristic design method based on optimal performance
CN113091731A (en) * 2021-03-03 2021-07-09 北京控制工程研究所 Spacecraft autonomous navigation method based on star sight relativistic effect
CN113128032A (en) * 2021-04-01 2021-07-16 北京航空航天大学 Intersection point time and position solving algorithm based on orbit analysis perturbation solution

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040046691A1 (en) * 2002-08-22 2004-03-11 Barker Lee A. Method for optimizing spacecraft yaw pointing to minimize geometric pointing error in antenna systems

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109460055A (en) * 2018-10-30 2019-03-12 中国运载火箭技术研究院 A kind of flying vehicles control capacity judging method, device and electronic equipment
CN111310277A (en) * 2020-03-13 2020-06-19 中国运载火箭技术研究院 Modeling method for pipeline transfer characteristics of atmospheric data sensing system, aircraft and storage medium
CN111680354A (en) * 2020-04-20 2020-09-18 北京航空航天大学 Method for calculating self-intersection point of orbit of near-earth regression orbit satellite subsatellite point and photographing point
CN111874267A (en) * 2020-04-30 2020-11-03 中国人民解放军战略支援部队航天工程大学 Low-orbit satellite off-orbit control method and system based on particle swarm optimization
CN111680462A (en) * 2020-08-11 2020-09-18 北京控制与电子技术研究所 Guidance method and system based on position change of space target in optical phase plane
CN112389672A (en) * 2020-11-09 2021-02-23 中国运载火箭技术研究院 Aerospace vehicle transverse course stability control characteristic design method based on optimal performance
CN112363410A (en) * 2020-11-13 2021-02-12 浙江大学 Intelligent autonomous control research and verification system for spacecraft
CN113091731A (en) * 2021-03-03 2021-07-09 北京控制工程研究所 Spacecraft autonomous navigation method based on star sight relativistic effect
CN113128032A (en) * 2021-04-01 2021-07-16 北京航空航天大学 Intersection point time and position solving algorithm based on orbit analysis perturbation solution

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
大型航天器无控飞行再入时间短期预报的轨道摄动方法研究;高兴龙 等;《载人航天》;20201031;第26卷(第5期);第566-573页 *
大型航天器无控飞行轨道衰降预报初步研究;高兴龙 等;《飞行力学》;20201231;第38卷(第6期);第70-76页 *

Also Published As

Publication number Publication date
CN113343369A (en) 2021-09-03

Similar Documents

Publication Publication Date Title
CN113343369B (en) Perturbation analysis method for spacecraft aerodynamic fusion orbit
CN110595485A (en) Low-orbit satellite long-term orbit forecasting method based on two-line root number
Thomson Introduction to space dynamics
CN109255096B (en) Geosynchronous satellite orbit uncertain evolution method based on differential algebra
CN109592079A (en) A kind of spacecraft coplanar encounter of limiting time becomes rail strategy and determines method
CN107861517A (en) The online trajectory planning method of guidance of great-jump-forward reentry vehicle based on linear pseudo- spectrum
Carrara An open source satellite attitude and orbit simulator toolbox for Matlab
CN110068845B (en) Method for determining theoretical orbit of satellite based on flat root theory
CN114580224B (en) Distributed pneumatic fusion track coupling attitude perturbation analysis method
Piñeros et al. Analysis of the orbit lifetime of CubeSats in low Earth orbits including periodic variation in drag due to attitude motion
CN102116630A (en) Mars probe on-board quick and high-precision determination method
Davis et al. Lunar impact probability for spacecraft in near rectilinear halo orbits
Sun et al. Satellite attitude identification and prediction based on neural network compensation
CN105799949B (en) A kind of pressure heart design method, attitude control method and the system of Asia orbiter
Desai et al. Entry, descent, and landing operations analysis for the genesis entry capsule
CN112389679B (en) Moonlet constant thrust orbit recursion method considering multiple perturbation forces
Batcheldor et al. Lunar regolith trajectories as a result of plume surface interactions
Falcone et al. Closed-Form Trajectory Solution for Shallow, High-Altitude Atmospheric Flight
CN111382514A (en) Mars detection flight orbit accurate calculation method and system based on supervised learning
Hematulin et al. Cooperative Motion Planning for multiple uavs via the Bézier curve guided line of sight techniques
CN114386282B (en) Low-orbit giant constellation orbit dynamics analysis method and device by semi-analysis method
CN113282097B (en) Control method for relative position non-spherical perturbation error of GEO game spacecraft
Moreschi et al. Aerodynamic resistance in upper atmosphere: case of the last stage Delta rocket fall in Argentina
CN112034704B (en) Dynamic gain online estimation method
Wu et al. Nonlinear dynamic modeling and simulation of the re-entry crew return vehicle

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant