CN114357788A - Low-orbit giant constellation deviation evolution analysis method and device - Google Patents
Low-orbit giant constellation deviation evolution analysis method and device Download PDFInfo
- Publication number
- CN114357788A CN114357788A CN202210022728.9A CN202210022728A CN114357788A CN 114357788 A CN114357788 A CN 114357788A CN 202210022728 A CN202210022728 A CN 202210022728A CN 114357788 A CN114357788 A CN 114357788A
- Authority
- CN
- China
- Prior art keywords
- constellation
- analyzed
- flat
- perturbation
- target
- 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.)
- Granted
Links
Images
Landscapes
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
The invention provides a low-orbit giant constellation deviation evolution analysis method and a device, which relate to the technical field of satellites and are used for acquiring target perturbation; based on a semi-analytical algorithm, creating a constellation dynamics model according to the target perturbation, wherein parameters of the constellation dynamics model comprise a flat number, and the constellation dynamics model is used for representing the motion state of a constellation; acquiring N initial flat roots of a constellation to be analyzed, and inputting the N initial flat roots into a constellation dynamics model; obtaining the number of to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model to obtain N number of to-be-analyzed flat roots; fitting N flat roots to be analyzed to obtain a target curved surface based on a Monte Carlo algorithm; and carrying out evolution analysis on the motion state of the constellation to be analyzed according to the target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit, and carrying out evolution analysis on the satellite by combining a semi-analysis algorithm and a Monte Carlo algorithm, so that the efficiency of the evolution analysis is improved.
Description
Technical Field
The invention relates to the technical field of satellites, in particular to a low-orbit giant constellation deviation evolution analysis method and device.
Background
With the development of satellite miniaturization, light weight and one-rocket-multi-satellite technology, it is a great trend to build giant constellations in low-orbit space. The problem of constellation operation safety brought by the problem is paid much attention to, and any small collision in the celestial body can possibly cause catastrophic cascade collision. Due to the influence of complex perturbation, divergence may occur in the measured constellation phase due to errors in the measured data, the dynamical model, the solving algorithm and the actuator. Therefore, the constellation needs to be analyzed and researched during orbital operation, for example, the deviation evolution of the constellation configuration needs to be analyzed and researched.
In the related art, for a general random power orbit system, an evolution equation of probability density along with time satisfies a Fokker-plane equation, and for an orbit dynamics system (generally 6 dimensions) with higher dimension and nonlinear perturbation, solving of the Fokker-plane equation becomes very difficult, so that probability density evolution often needs to perform local linearization assumption and evolution analysis only on low-order moments (such as a mean value and a covariance matrix). In the evolution of the deviations of orbital dynamics, the same row assumes that the initial error is gaussian distributed, the gaussian error ellipsoid is rotationally distorted over time and may become non-gaussian due to non-linear factors in the dynamics. The main existing bias evolution methods include a polynomial chaotic expansion method, a state transition tensor method, a differential algebra method, a mixed Gaussian model method and the like.
However, in the related art, the process of evolution of the constellation deviation is complex in calculation and large in calculation amount, so that the analysis efficiency is low.
In order to achieve the above object, the present invention provides a method and a device for analyzing evolution of low-orbit giant constellation deviation, aiming to simplify the evolution analysis process and improve the evolution analysis efficiency, the specific technical scheme is as follows:
in a first aspect of the embodiments of the present invention, a method for analyzing evolution of low-orbit giant constellation deviation is provided, in which a target perturbation is obtained; based on a semi-analysis algorithm, creating a constellation dynamics model according to the target perturbation, wherein parameters of the constellation dynamics model comprise a flat root number, and the constellation dynamics model is used for representing the motion state of a constellation; acquiring N initial flat roots of a constellation to be analyzed, and inputting the N initial flat roots into the constellation kinetic model; obtaining the number of to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model to obtain N number of to-be-analyzed flat roots; fitting the N planar numbers to be analyzed to obtain a target curved surface based on a Monte Carlo algorithm; and carrying out evolution analysis on the motion state of the constellation to be analyzed according to the target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit.
Optionally, the acquiring the target perturbation comprises: acquiring the non-spherical gravitational perturbation of the earth and acquiring the atmospheric resistance perturbation.
Optionally, the acquiring the non-spherical gravitational perturbation of the earth comprises: and acquiring a first-order long-term of the global non-spherical gravity perturbation and a second-order long-term of the global non-spherical gravity perturbation.
Optionally, the acquiring the atmospheric resistance perturbation further comprises: acquiring a long term of the atmospheric resistance perturbation.
Optionally, the fitting the N planar numbers to be analyzed to the target curved surface based on a Monte Carlo algorithm includes: fitting the N planar numbers to be analyzed to form a middle curved surface based on the Monte Carlo algorithm; and fitting the intermediate curved surface by adopting a least square method to obtain the target curved surface.
Optionally, the obtaining of the N initial flat roots of the constellation to be analyzed includes: acquiring the actual flat root of the first constellation to be analyzed and the actual flat root of the second constellation to be analyzed; converting the actual flat root into an actual transient root through a Brouwer flat transient conversion algorithm, wherein the actual transient root comprises the running time of the constellation to be analyzed and a target standard deviation between the first constellation to be analyzed and the second constellation to be analyzed; obtaining a deviation transient root according to the target standard deviation and the actual transient root, wherein the target standard deviation comprises a position standard deviation between the first constellation to be analyzed and the second constellation to be analyzed and/or a velocity standard deviation between the first constellation to be analyzed and the second constellation to be analyzed; generating a transient root normal distribution relation according to the deviation transient root, wherein the transient root normal distribution relation represents the corresponding relation between the constellation operation time to be analyzed and the target standard deviation; acquiring N deviation transient numbers from the normal distribution relation; and converting the N deviation transient numbers into the N initial flat numbers according to the Brouwer flat transient conversion algorithm.
Optionally, converting the N deviation transient numbers into the N initial flat numbers according to the Brouwer flat transient conversion algorithm, including: acquiring a first-order periodic item of the earth non-spherical gravity perturbation and acquiring an additional perturbation item of a coordinate system; and obtaining the N initial flat root numbers according to the first-order long-term of the earth aspheric gravity perturbation, the second-order long-term of the earth aspheric gravity perturbation, the first-order periodic term of the earth aspheric gravity perturbation, the long-term of the atmospheric resistance perturbation, the coordinate system additional perturbation term and the N deviation transient root numbers on the basis of the Brouwer flat transient conversion algorithm.
Optionally, the target standard deviation comprises a position standard deviation and/or a velocity standard deviation.
Optionally, the target standard deviation further comprises a phase standard deviation.
In a first aspect of the embodiments of the present invention, there is provided a low-orbit giant constellation deviation evolution analysis apparatus, including: the acquisition module is used for acquiring the target perturbation; the creating module is used for creating a constellation dynamics model according to the target perturbation based on a semi-analysis algorithm, wherein parameters of the constellation dynamics model comprise a flat number, and the constellation dynamics model is used for representing the motion state of a constellation; the input module is used for acquiring N initial flat roots of a constellation to be analyzed and inputting the N initial flat roots into the constellation kinetic model; the output module is used for obtaining the to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model to obtain N to-be-analyzed flat roots; the fitting module is used for fitting the N flat roots to be analyzed to form a target curved surface based on a Monte Carlo algorithm; and the analysis module is used for carrying out evolution analysis on the motion state of the constellation to be analyzed according to the target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit.
Compared with the prior art, the invention provides a low-orbit giant constellation deviation evolution analysis method and a device, which have the following beneficial effects:
firstly, acquiring target perturbation; based on a semi-analysis algorithm, a constellation dynamics model is created according to the target perturbation, wherein parameters of the constellation dynamics model comprise a flat number, and the constellation dynamics model is used for representing the motion state of a constellation; acquiring N initial flat roots of a constellation to be analyzed, and inputting the N initial flat roots into a constellation dynamics model; obtaining the number of to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model to obtain N number of to-be-analyzed flat roots; fitting N flat roots to be analyzed to obtain a target curved surface based on a Monte Carlo algorithm; the method comprises the steps of carrying out evolution analysis on the motion state of a constellation to be analyzed according to a target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit, and a semi-analysis algorithm has higher calculation efficiency, so that a constellation dynamic model can be quickly obtained, and a Monte Carlo algorithm does not need to operate and process the obtained constellation dynamic model.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below. It is obvious that the drawings in the following description are only some embodiments of the invention, and that for a person skilled in the art, other drawings can be derived from them without inventive effort.
Fig. 1 is a flowchart illustrating a method for analyzing evolution of low-orbit giant constellation deviation according to an embodiment of the present disclosure;
FIG. 2 is a flow chart illustrating the sub-steps of step S130 of FIG. 1 of the present application;
fig. 3 is a flowchart illustrating a method for analyzing evolution of low-orbit giant constellation deviation according to another embodiment of the present application.
Detailed Description
The technical solution in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. It is to be understood that the embodiments described are only a few embodiments of the present invention, and not all embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
With the development of satellite miniaturization, light weight and one-rocket-multi-satellite technology, it is a great trend to build giant constellations in low-orbit space. The problem of constellation operation safety brought by the problem is paid much attention to, and any small collision in the celestial body can possibly cause catastrophic cascade collision. Due to the influence of complex perturbation, the constellation phase may diverge due to the existence of errors in the measured data, the dynamical model, the solving algorithm and the actuator. Therefore, the constellation needs to be analyzed and researched during orbital operation, for example, the deviation evolution of the constellation configuration needs to be analyzed and researched.
In the related art, for a general random dynamic system, an evolution equation of probability density with time satisfies a Fokker-plane equation, and for an orbital dynamic system (generally 6 dimensions) with higher dimension and nonlinear perturbation, solving of the Fokker-plane equation becomes very difficult, so that probability density evolution often needs to perform local linearization assumption and evolution analysis only on low-order moments (such as a mean value matrix and a covariance matrix). In the evolution of the deviations of orbital dynamics, the same row assumes that the initial error is gaussian distributed, the gaussian error ellipsoid is rotationally distorted over time and may become non-gaussian due to non-linear factors in the dynamics. The main existing bias evolution methods include a polynomial chaotic expansion method, a state transition tensor method, a differential algebra method, a mixed Gaussian model method and the like.
However, in the related art, the process of evolution of the constellation deviation is complex in calculation and large in calculation amount, so that the analysis efficiency is low.
In view of this, the present invention provides a method and an apparatus for analyzing evolution of low-orbit giant constellation deviation, which are provided by the following embodiments and are intended to simplify the analysis process and improve the analysis efficiency.
Fig. 1 is a flowchart illustrating a method for analyzing evolution of low-orbit giant constellation deviation according to an embodiment of the present application, please refer to fig. 1, where the method for analyzing evolution of low-orbit giant constellation deviation includes the following steps:
and step S110, acquiring the target perturbation.
The perturbation refers to a deviation generated on a track where a celestial body runs due to the attraction of other celestial bodies or the influence of other factors when the celestial body moves around another celestial body according to the rule of a two-body problem, and the perturbation is understandably a perturbation force acting on the celestial body.
In this embodiment, the object perturbation represents the influence of other celestial bodies on the forces caused by the constellation to be analyzed, and the constellation to be analyzed may refer to an artificial giant satellite launched in space and used for military object detection, mobile phone communication and the like. Moreover, the constellation to be analyzed in this embodiment may be a low-orbit constellation, wherein the low-orbit constellation is closer to the ground than a conventional communication satellite, for example, the low-orbit constellation may run on an orbit 200 to 2000 km from the ground.
In the present embodiment, the object perturbation includes an earth aspherical gravity perturbation J2 and an atmospheric resistance perturbation, the earth aspherical gravity perturbation is acquired, and the atmospheric resistance perturbation is acquired. The non-spherical gravitational perturbation of the earth refers to perturbation generated by non-uniform density distribution in the earth due to the fact that the earth is not a regular sphere but an ellipsoid. The atmospheric resistance perturbation refers to the perturbation generated by the atmospheric resistance on the earth surface when a constellation to be analyzed runs on a spatial orbit close to the earth.
Optionally, after the target perturbation is obtained, a target perturbation item may be obtained according to the target perturbation, where the target perturbation item represents an index of the target perturbation, and under the effect of the target perturbation, elements such as coordinates, a motion speed, a motion orbit, and the like of a constellation to be analyzed are changed, so that a component of the change becomes the target perturbation item.
In one embodiment, when the target perturbation comprises an earth aspheric gravitational perturbation, acquiring the earth aspheric gravitational perturbation term comprises acquiring a first order long term of the earth aspheric gravitational perturbation and a second order long term of the earth aspheric gravitational perturbation.
In another embodiment, when the target perturbation includes an atmospheric resistance perturbation, acquiring an atmospheric resistance perturbation term includes acquiring a long term of the atmospheric resistance perturbation.
The target perturbation may be a gravitational perturbation, in addition to the above-described non-gravitational perturbation such as an earth non-spherical gravitational perturbation and an atmospheric resistance perturbation, and is not particularly limited herein. The gravitational perturbation may include, but is not limited to, an earth gravitational perturbation.
And S120, based on a semi-analysis algorithm, creating a constellation dynamics model according to the target perturbation, wherein parameters of the constellation dynamics model comprise a flat number, the constellation dynamics model is used for representing the motion state of a constellation, and the constellation to be analyzed is subjected to dynamic analysis through the constellation dynamics model.
Among the tools for the semi-analytical method are STELA (semi-analytical Satellite Life Analysis kit for End of Life Analysis) and DSST (Draper semi-analytical Tool kit for health) and the like.
In one embodiment, the target perturbation comprises an earth non-spherical gravity perturbation J2Perturbation of atmospheric resistance, perturbation of J by global non-spherical attraction2And atmospheric resistance perturbation creating a constellation dynamics model. As one mode, acquiring the aspheric gravity perturbation J of the earth2First-order long-term of (J), global non-spherical gravity perturbation2And a second-order long-term of the atmospheric resistance perturbation, and creating a constellation dynamics model according to the three.
Step S130, obtaining N initial flat roots of the constellation to be analyzed, and inputting the N initial flat roots into the constellation dynamics model.
In one embodiment, the constellation to be analyzed includes a first constellation to be analyzed and a second constellation to be analyzed that is in the same orbital plane as the first constellation to be analyzed, where a phase difference exists between the first constellation to be analyzed and the second constellation to be analyzed, and the phase difference is Δ Φ, that is, the phase of the first constellation to be analyzed is different from the phase of the second constellation to be analyzed. The actual average numbers of the first constellation to be analyzed and the second constellation to be analyzed are the same, the first constellation to be analyzed is used as a reference constellation, and when the constellation to be analyzed includes the first constellation to be analyzed and the second constellation to be analyzed, the low-orbit giant constellation deviation evolution analysis method in this embodiment analyzes the change of the motion state between the first constellation to be analyzed and the second constellation to be analyzed along with the motion of the constellation, fig. 2 shows a flowchart of the substep of step S130 in fig. 1 of the present application, and step S130 includes the substeps of:
and a substep S131, obtaining the actual flat root of the first constellation to be analyzed and the actual flat root of the second constellation to be analyzed.
The flat root represents an average value of motion parameters of the constellation to be analyzed, for example, when the flat root is a motion speed of the constellation to be analyzed in the space, the flat root represents an average speed of the constellation to be analyzed in the space, which runs around an orbit for one circle; for another example, when the flat root is a motion period of the constellation to be analyzed in the space, the flat root represents an average period of the constellation to be analyzed in the space; for another example, when the flat root is the angular velocity of the constellation to be analyzed in the space, the flat root represents the average angular velocity of the constellation to be analyzed that orbits around the space for one circle. The actual flat number of the first constellation to be analyzed is the same as the actual flat number of the second constellation to be analyzed, and it can be understood that the initial motion states of both the first constellation to be analyzed and the second constellation to be analyzed are the same.
And a substep S132 of converting the actual flat root into an actual transient root through a Brouwer flat transient conversion algorithm, wherein the actual transient root includes the running time of the constellation to be analyzed and the target standard deviation between the first constellation to be analyzed and the second constellation to be analyzed.
It is understood that the transient root represents an instantaneous value of a motion parameter of the constellation to be analyzed, for example, when the transient root represents a motion speed of the constellation to be analyzed in space, the transient root represents an instantaneous speed of the constellation to be analyzed orbiting to a certain position in space; for another example, when the transient number is a motion period of the constellation to be analyzed in space, the transient number represents a time when the constellation to be analyzed moves to a certain position on a trajectory; for another example, when the instantaneous number is the angular velocity of the constellation to be analyzed in the space, the instantaneous number represents the instantaneous angular velocity of the constellation to be analyzed moving to a certain position around the orbit in the space.
The Brouwer flat transient conversion algorithm can realize flat transient conversion, that is, conversion between a flat root and an instantaneous root, that is, the flat root can be converted into the instantaneous root through the Brouwer flat transient conversion algorithm, and the instantaneous root can also be converted into the flat root through the Brouwer flat transient conversion algorithm.
As an implementation mode, the method obtains the earth non-spherical gravity perturbation J2The first order short period term of (1), using the Brouwer flat transient conversion algorithm J2The first order short period term of (2) is smoothly transformed.
In the present embodiment, the actual flat root is converted into the actual transient root, and the transient root is obtained such that the actual transient root is an unbiased transient root.
It should be noted that, in addition to the Brouwer flat transient transformation algorithm, a software Tool can be used to perform the flat transient transformation, for example, the Tool can be STK9 (Satellite toolkit, Satellite toolkit Tool 9).
And a substep S133, obtaining a deviation transient number according to the target standard deviation and the actual transient number, wherein the target standard deviation includes a position standard deviation between the first constellation to be analyzed and the second constellation to be analyzed and/or a velocity standard deviation between the first constellation to be analyzed and the second constellation to be analyzed.
Wherein the deviation transient is deviated from the actual transient.
In one embodiment, when the target standard deviation includes a position standard deviation between the first constellation to be analyzed and the second constellation to be analyzed, the deviation transient number includes a position component, and the position standard deviation is added to the position components of the two satellites of the first constellation to be analyzed and the second constellation to be analyzed, so as to obtain the deviation transient numbers corresponding to the two satellites.
In another embodiment, when the target standard deviation includes a velocity standard deviation between the first constellation to be analyzed and the second constellation to be analyzed, the deviation transient number includes a velocity component, and the velocity standard deviation is added to the respective velocity components of the two satellites of the first constellation to be analyzed and the second constellation to be analyzed, so as to obtain the deviation transient numbers corresponding to the two satellites.
And a substep S134, generating a transient root normal distribution relation according to the deviation transient root, wherein the transient root normal distribution relation represents the corresponding relation between the constellation operation time to be analyzed and the target standard deviation.
Optionally, the target standard deviation further includes a satellite phase standard deviation, the satellite phase is expressed by a sum of an angle distance of an apogee and an angle of the apogee, and the transient normal distribution relationship represents a corresponding relationship among the running time, the position component, the velocity component and the phase standard deviation of the constellation to be analyzed.
And a substep S135, acquiring N deviation transient numbers from the normal distribution relation.
Acquiring at least one of the running time, the position component, the speed component and the phase standard deviation of the constellation to be analyzed corresponding to N positions of the orbit of the constellation to be analyzed, and acquiring other deviation transient values corresponding to each position of the N positions according to the normal distribution relation. For example, when the running time of the constellation to be analyzed is obtained, based on the normal distribution relationship, the position component, the velocity component and the phase standard deviation corresponding to the running time of the constellation to be analyzed are obtained.
And a substep S136, converting the N deviation transient numbers into the N initial flat numbers according to the Brouwer flat transient conversion algorithm. As an embodiment, acquiring a first-order periodic item of the earth non-spherical gravity perturbation, and acquiring a coordinate system additional perturbation item; and obtaining the N initial flat root numbers according to the first-order long-term of the earth aspheric gravity perturbation, the second-order long-term of the earth aspheric gravity perturbation, the first-order periodic term of the earth aspheric gravity perturbation, the long-term of the atmospheric resistance perturbation, the coordinate system additional perturbation term and the N deviation transient root numbers on the basis of the Brouwer flat transient conversion algorithm.
It should be noted that, the Brouwer flat transient conversion algorithm can implement flat transient conversion, which has been described in detail in the sub-step S132, and please refer to the sub-step S132, which is not described herein again.
Step S140, obtaining the numbers of the to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model, and obtaining the N numbers of the to-be-analyzed flat roots.
The N initial flat roots are input into the constellation dynamics model and used for analyzing the motion parameters through the constellation dynamics model, the N initial flat roots output by the constellation dynamics model correspond to the flat roots to be analyzed, each initial flat root corresponds to one flat root to be analyzed, and correspondingly, the N initial flat roots correspond to the N flat roots to be analyzed.
And S150, fitting the N flat roots to be analyzed to form a target curved surface based on a Monte Carlo algorithm.
In one embodiment, when the transient deviation number includes a position component, based on Monte Carlo algorithm, a position component-to-be-analyzed constellation runtime-phase standard deviation surface Γ (t, σ (r), σ (Δ Φ)) is fitted according to the position component, the to-be-analyzed constellation runtime, and the phase standard deviation, and it can be understood that the surface Γ (t, σ (r), σ (Δ Φ)) represents a mapping relationship among the position component, the to-be-analyzed constellation runtime, and the phase standard deviation.
In another embodiment, when the deviation transient includes a velocity component, a position component-constellation runtime to be analyzed-phase standard deviation surface Γ (t, σ (r), σ (Δ Φ)) is fitted based on the velocity component, the constellation runtime to be analyzed, and the phase standard deviation based on a Monte Carlo algorithm. It can be understood that Γ (t, σ (v), σ (Δ Φ)) characterizes the mapping of the velocity component, the constellation runtimes to be analyzed, and the phase standard deviation.
Assuming that only position errors exist in constellation station position keeping control, 4000 normally distributed samples are generated, instantaneous transformation is used for adding the samples into a flat root number, t-day track extrapolation is carried out in a dynamic model, the standard deviation of an initial sample is changed, and Monte Carlo simulation is repeated for multiple times.
In yet another embodiment, the N number of flat surfaces to be analyzed are first fitted to an intermediate surface based on the Monte Carlo algorithm (i.e., Monte Carlo algorithm). And fitting the intermediate curved surface by adopting a least square method to obtain the target curved surface. The Monte Carlo algorithm is also called a statistical simulation method, and solves the calculation problem by using random numbers (or pseudo-random numbers), and is an important numerical calculation method.
Firstly, acquiring the flat root of a constellation to be analyzed, establishing two satellites (namely a first constellation to be analyzed and a second constellation to be analyzed) with fixed phase difference by using the flat root, and perturbing J by using the non-spherical gravity of the earth2First order short period term of term κspAnd a first order long period term κipConstructing a nonlinear equation system F:
wherein
The mean number of satellites and κ (t) the instantaneous number.
And then, Newton iteration gives a transient number corresponding to the flat number. KappaspThe expression of each Kepler root number in the formula is as follows:
wherein (a, e, i, omega, M) are respectively a track semimajor axis, a track inclination angle, a rising intersection declination, a perigee angular distance and a perigee angle,
is a semi-major axis J2The first order short period term is used,
is eccentricity J2The first order short period term is used,
at an angle of inclination of the track J2The first order short period term is used,
j as the right ascension of the ascending crossing point2The first order short period term is used,
is a distance of a point2The first order short period term is used,
at a mean angle of approach J2A first order short period term, f is the true periapical angle, J2 is 0.00108263, p is a (1-e)2) Is a half-diameter. J. the design is a square2First order long period term
The expression of the number of the elements is as follows:
wherein (a, e, i, omega, M) are respectively a track semimajor axis, a track inclination angle, a rising intersection declination, a perigee angular distance and a perigee angle,
is a semi-major axis J2A first order long period term is included,
Secondly, given the standard deviation σ of position or velocity under the assumption that the station position holding control has only position or velocity errorrAnd σvGenerating a sample N (0, σ) in which the position error or the velocity error is normally distributedr) And N (0, σ)v). The biased position velocity is then converted into a biased transient number, and further into a biased flat number. With J2The first-order second-order long-term and the atmospheric resistance second-order long-term are dynamic models, and samples of two stars are extrapolated for t days. Wherein J2First order long term
Comprises the following steps:
wherein, (a, e, i, omega, M) is respectively a track semimajor axis, a track inclination angle, a rising intersection declination, an apogee angular distance and an apogee angle, J2＝0.00108263，p＝a(1-e2) Is a half-path, n is the average motion angular velocity of the target satellite,
is a semi-major axis J2The first-order long-term,
is eccentricity J2The first-order long-term,
at an angle of inclination of the track J2The first-order long-term,
j as the right ascension of the ascending crossing point2The first-order long-term,
is a distance of a point2The first-order long-term,
at a mean angle of approach J2The first order long term, i, is the track dip.
J2Second order long term
Comprises the following steps:
wherein, (a, e, i, omega, M) are respectively a track semimajor axis, a track inclination angle, a rising intersection declination, an angle distance of a near point and an angle of a near point,
is a semi-major axis J2The second-order long-term,
is eccentricity J2The second-order long-term,
at an angle of inclination of the track J2The second-order long-term,
Atmospheric resistance long term κdrag，secComprises the following steps:
wherein, (a, e, I, omega, M) is respectively a track semimajor axis, a track inclination angle, a rising intersection declination, an angle distance between a near point and an angle between a near point, I0、I1、I2、I3、I4First, second, third and fourth order first class imaginary variable Bessel function, asec，dragLong term atmospheric resistance for semi-major axis, esec，dragLong term of atmospheric resistance, i, of eccentricitysec，dragLong term of atmospheric resistance, omega, for track inclinationsec，dragLong term of atmospheric resistance, omega, for the ascension crossing right ascensionsec，dragLong term of atmospheric resistance for angular distance of near-earth, Msec，dragThe long term for atmospheric resistance at the mean anomaly angle. I ism(z) is a Bessel function of a first type of imaginary variable, z being its argument
B1、B2F, C are intermediate variables:
wherein C isDIs a drag coefficient, S is the cross-sectional area of the satellite, m is the mass of the satellite, v is 0.1 the variability of the atmospheric density elevation, ne＝7.2921158553×10-5rad/s is the Earth rotation speed, (a)0，e0，i0，ω0) Is the initial flat number, pP0And HP0Is the atmospheric density and density elevation at the initial close location, rP0And vP0Is the position and velocity at the initial near point.
The two-star phase difference distribution obtained by Monte Carlo algorithm is used for counting the standard deviation sigma (delta phi), and the extrapolation time t and the initial position speed standard deviation sigma are drawnrOr σvAnd the standard deviation sigma (delta phi) of the relative phase. The time sequence is x, the initial position or speed error is y, the standard deviation of the relative phase difference is z, and the (x, y, z) are row vectors, so as to
For the basis function, a least squares surface fit is performed based on the following function:
construction with cpqSurface cluster Γ (x, y) which is a coefficient array:
cpqexpressed in the form of a matrix C:
C＝(BTB)-1BTzG(GTG)-1
wherein the expressions of the matrices B and G are
In the above, an approximate relationship of the constellation initial state error to the terminal state error is established.
And S160, carrying out evolution analysis on the motion state of the constellation to be analyzed according to the target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit.
The flat roots include a plurality of types, and when the low-orbit giant constellation deviation evolution analysis is performed, according to a flat root (for example, time t1, t2), other flat roots (for example, phase difference and the like) corresponding to the flat root are obtained from the target curved surface, so that the deviation between the first generation analysis constellation and the second constellation to be analyzed after the time t2-t1 is obtained.
Optionally, the evolution analysis result may be displayed on a display screen of the electronic device, so as to be convenient for a researcher to browse.
In the low-orbit giant constellation deviation evolution analysis method provided by the embodiment, a target perturbation is obtained first; based on a semi-analysis algorithm, a constellation dynamics model is created according to the target perturbation, wherein parameters of the constellation dynamics model comprise a flat number, and the constellation dynamics model is used for representing the motion state of a constellation; acquiring N initial flat roots of a constellation to be analyzed, and inputting the N initial flat roots into a constellation dynamics model; obtaining the number of to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model to obtain N number of to-be-analyzed flat roots; fitting N flat roots to be analyzed to obtain a target curved surface based on a Monte Carlo algorithm; the method comprises the steps of carrying out evolution analysis on the motion state of a constellation to be analyzed according to a target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit, and a semi-analysis algorithm has higher calculation efficiency, so that a constellation dynamic model can be quickly obtained, and a Monte Carlo algorithm does not need to operate and process the obtained constellation dynamic model.
Fig. 3 shows a method for analyzing evolution of low-orbit giant constellation deviations according to an embodiment of the present application.
To implement the above method embodiments, this embodiment provides an apparatus for analyzing evolution of low-orbit giant constellation deviation, where the apparatus includes: the acquisition module is used for acquiring the target perturbation;
the creating module is used for creating a constellation dynamics model according to the target perturbation based on a semi-analysis algorithm, wherein parameters of the constellation dynamics model comprise a flat number, and the constellation dynamics model is used for representing the motion state of a constellation;
the input module is used for acquiring N initial flat roots of a constellation to be analyzed and inputting the N initial flat roots into the constellation kinetic model;
the output module is used for obtaining the to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model to obtain N to-be-analyzed flat roots;
the fitting module is used for fitting the N flat roots to be analyzed to form a target curved surface based on a Monte Carlo algorithm;
and the analysis module is used for carrying out evolution analysis on the motion state of the constellation to be analyzed according to the target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit.
It can be clearly understood by those skilled in the art that, for convenience and brevity of description, the specific working processes of the modules/units/sub-units/components in the above-described apparatus may refer to the corresponding processes in the foregoing method embodiments, and are not described herein again.
In the several embodiments provided in the present application, the coupling or direct coupling or communication connection between the modules shown or discussed may be through some interfaces, and the indirect coupling or communication connection between the devices or modules may be in an electrical, mechanical or other form.
In addition, functional modules in the embodiments of the present application may be integrated into one processing module, or each of the modules may exist alone physically, or two or more modules are integrated into one module. The integrated module can be realized in a hardware mode, and can also be realized in a software functional module mode.
In summary, the method and the device for analyzing evolution of low-orbit giant constellation deviation provided by the application first obtain a target perturbation; based on a semi-analysis algorithm, a constellation dynamics model is created according to the target perturbation, wherein parameters of the constellation dynamics model comprise a flat number, and the constellation dynamics model is used for representing the motion state of a constellation; acquiring N initial flat roots of a constellation to be analyzed, and inputting the N initial flat roots into a constellation dynamics model; obtaining the number of to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model to obtain N number of to-be-analyzed flat roots; fitting N flat roots to be analyzed to obtain a target curved surface based on a Monte Carlo algorithm; the method comprises the steps of carrying out evolution analysis on the motion state of a constellation to be analyzed according to a target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit, and a semi-analysis algorithm has higher calculation efficiency, so that a constellation dynamic model can be quickly obtained, and a Monte Carlo algorithm does not need to operate and process the obtained constellation dynamic model.
It is noted that, herein, relational terms such as first and second, and the like may be used solely to distinguish one entity or action from another entity or action without necessarily requiring or implying any actual such relationship or order between such entities or actions. Also, the terms "comprises," "comprising," or any other variation thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising an … …" does not exclude the presence of other identical elements in a process, method, article, or apparatus that comprises the element.
All the embodiments in the present specification are described in a related manner, and the same and similar parts among the embodiments may be referred to each other, and each embodiment focuses on the differences from the other embodiments. In particular, for the system embodiment, since it is substantially similar to the method embodiment, the description is simple, and for the relevant points, reference may be made to the partial description of the method embodiment.
The above description is only for the preferred embodiment of the present invention, and is not intended to limit the scope of the present invention. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention shall fall within the protection scope of the present invention.
Claims (10)
1. A method for analyzing evolution of low-orbit giant constellation deviation, the method comprising:
acquiring a target perturbation;
based on a semi-analysis algorithm, creating a constellation dynamics model according to the target perturbation, wherein parameters of the constellation dynamics model comprise a flat root number, and the constellation dynamics model is used for representing the motion state of a constellation;
acquiring N initial flat roots of a constellation to be analyzed, and inputting the N initial flat roots into the constellation kinetic model;
obtaining the number of to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model to obtain N number of to-be-analyzed flat roots;
fitting the N planar numbers to be analyzed to obtain a target curved surface based on a Monte Carlo algorithm;
and carrying out evolution analysis on the motion state of the constellation to be analyzed according to the target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit.
2. The method of claim 1, wherein the obtaining the perturbation of the target comprises:
acquiring the non-spherical gravitational perturbation of the earth and acquiring the atmospheric resistance perturbation.
3. The method of analyzing evolution of low-orbit giant constellation deviation as claimed in claim 2, wherein the obtaining perturbation of the aspheric gravitational field of the earth comprises:
and acquiring a first-order long-term of the global non-spherical gravity perturbation and a second-order long-term of the global non-spherical gravity perturbation.
4. The method of evolution analysis of low-orbit giant constellation deviation of claim 3, wherein the obtaining of atmospheric resistance perturbation further comprises:
acquiring a long term of the atmospheric resistance perturbation.
5. The method of claim 4, wherein the constellation to be analyzed includes a first constellation to be analyzed and a second constellation to be analyzed located on the same orbital plane as the first constellation to be analyzed, the first constellation to be analyzed and the second constellation to be analyzed have a phase difference therebetween and have the same actual flat number, and the method using the first constellation to be analyzed as a reference constellation and obtaining N initial flat numbers of the constellation to be analyzed comprises:
acquiring the actual flat root of the first constellation to be analyzed and the actual flat root of the second constellation to be analyzed;
converting the actual flat root into an actual transient root through a Brouwer flat transient conversion algorithm, wherein the actual transient root comprises the running time of the constellation to be analyzed and a target standard deviation between the first constellation to be analyzed and the second constellation to be analyzed;
obtaining a deviation transient root according to the target standard deviation and the actual transient root, wherein the target standard deviation comprises a position standard deviation between the first constellation to be analyzed and the second constellation to be analyzed and/or a velocity standard deviation between the first constellation to be analyzed and the second constellation to be analyzed;
generating a transient root normal distribution relation according to the deviation transient root, wherein the transient root normal distribution relation represents the corresponding relation between the constellation operation time to be analyzed and the target standard deviation;
acquiring N deviation transient numbers from the normal distribution relation;
and converting the N deviation transient numbers into the N initial flat numbers according to the Brouwer flat transient conversion algorithm.
6. The method of claim 5, wherein transforming the N bias transient numbers into the N initial flat numbers according to the Brouwer flat transient transformation algorithm comprises:
acquiring a first-order periodic item of the earth non-spherical gravity perturbation and acquiring an additional perturbation item of a coordinate system;
and obtaining the N initial flat root numbers according to the first-order long-term of the earth aspheric gravity perturbation, the second-order long-term of the earth aspheric gravity perturbation, the first-order periodic term of the earth aspheric gravity perturbation, the long-term of the atmospheric resistance perturbation, the coordinate system additional perturbation term and the N deviation transient root numbers on the basis of the Brouwer flat transient conversion algorithm.
7. The method of claim 6, wherein the target standard deviation comprises a position standard deviation and/or a velocity standard deviation.
8. The method of claim 7, wherein the target standard deviation further comprises a phase standard deviation.
9. The method for evolution analysis of low-orbit giant constellation deviation according to any one of claims 1 to 8, wherein the fitting of the N number of flat elements to be analyzed to a target surface based on Monte Carlo algorithm comprises:
fitting the N planar numbers to be analyzed to form a middle curved surface based on the Monte Carlo algorithm;
and fitting the intermediate curved surface by adopting a least square method to obtain the target curved surface.
10. An apparatus for analyzing evolution of constellation deviations in giant low-orbit, the apparatus comprising:
the acquisition module is used for acquiring the target perturbation;
the creating module is used for creating a constellation dynamics model according to the target perturbation based on a semi-analysis algorithm, wherein parameters of the constellation dynamics model comprise a flat number, and the constellation dynamics model is used for representing the motion state of a constellation;
the input module is used for acquiring N initial flat roots of a constellation to be analyzed and inputting the N initial flat roots into the constellation kinetic model;
the output module is used for obtaining the to-be-analyzed flat roots corresponding to the N initial flat roots output by the constellation dynamics model to obtain N to-be-analyzed flat roots;
the fitting module is used for fitting the N flat roots to be analyzed to form a target curved surface based on a Monte Carlo algorithm;
and the analysis module is used for carrying out evolution analysis on the motion state of the constellation to be analyzed according to the target curved surface, wherein the constellation to be analyzed is a satellite in a low orbit.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210022728.9A CN114357788B (en) | 2022-01-10 | 2022-01-10 | Low-orbit giant constellation deviation evolution analysis method and device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210022728.9A CN114357788B (en) | 2022-01-10 | 2022-01-10 | Low-orbit giant constellation deviation evolution analysis method and device |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114357788A true CN114357788A (en) | 2022-04-15 |
| CN114357788B CN114357788B (en) | 2023-08-01 |
Family
ID=81109809
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210022728.9A Active CN114357788B (en) | 2022-01-10 | 2022-01-10 | Low-orbit giant constellation deviation evolution analysis method and device |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114357788B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115098828A (en) * | 2022-08-26 | 2022-09-23 | 北京控制工程研究所 | Method and device for calculating low-orbit satellite orbit in near circle |
| CN115806060A (en) * | 2022-11-10 | 2023-03-17 | 北京航天驭星科技有限公司 | Modeling method, modeling method and obtaining method of satellite relative phase holding strategy model |
| CN115892516A (en) * | 2022-11-10 | 2023-04-04 | 北京航天驭星科技有限公司 | Modeling method, modeling method and obtaining method of satellite relative phase holding strategy model |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160009425A1 (en) * | 2014-07-10 | 2016-01-14 | The Aerospace Corporation | Systems and Methods for Optimizing Satellite Constellation Deployment |
| CN107451319A (en) * | 2017-05-05 | 2017-12-08 | 中国科学院国家天文台 | A kind of modeling method of space debris environment long-term evolution model |
| CN109613574A (en) * | 2018-11-13 | 2019-04-12 | 中国人民解放军战略支援部队航天工程大学 | Calculate the method that rail satellite grave track passes through other Global Satellite Navigation System orbit times earliest in Beidou |
| CN111854765A (en) * | 2020-06-08 | 2020-10-30 | 中国人民解放军战略支援部队航天工程大学 | Medium-orbit navigation satellite orbit long-term forecasting method |
-
2022
- 2022-01-10 CN CN202210022728.9A patent/CN114357788B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160009425A1 (en) * | 2014-07-10 | 2016-01-14 | The Aerospace Corporation | Systems and Methods for Optimizing Satellite Constellation Deployment |
| CN107451319A (en) * | 2017-05-05 | 2017-12-08 | 中国科学院国家天文台 | A kind of modeling method of space debris environment long-term evolution model |
| CN109613574A (en) * | 2018-11-13 | 2019-04-12 | 中国人民解放军战略支援部队航天工程大学 | Calculate the method that rail satellite grave track passes through other Global Satellite Navigation System orbit times earliest in Beidou |
| CN111854765A (en) * | 2020-06-08 | 2020-10-30 | 中国人民解放军战略支援部队航天工程大学 | Medium-orbit navigation satellite orbit long-term forecasting method |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115098828A (en) * | 2022-08-26 | 2022-09-23 | 北京控制工程研究所 | Method and device for calculating low-orbit satellite orbit in near circle |
| CN115098828B (en) * | 2022-08-26 | 2022-11-04 | 北京控制工程研究所 | Method and device for calculating low-orbit satellite orbit in near circle |
| CN115806060A (en) * | 2022-11-10 | 2023-03-17 | 北京航天驭星科技有限公司 | Modeling method, modeling method and obtaining method of satellite relative phase holding strategy model |
| CN115892516A (en) * | 2022-11-10 | 2023-04-04 | 北京航天驭星科技有限公司 | Modeling method, modeling method and obtaining method of satellite relative phase holding strategy model |
| CN115806060B (en) * | 2022-11-10 | 2023-05-19 | 北京航天驭星科技有限公司 | Modeling method, model and acquisition method of satellite relative phase maintaining strategy model |
| CN115892516B (en) * | 2022-11-10 | 2023-06-20 | 北京航天驭星科技有限公司 | Modeling method, model and acquisition method of satellite relative phase maintaining strategy model |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114357788B (en) | 2023-08-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114357788A (en) | Low-orbit giant constellation deviation evolution analysis method and device | |
| CN109255096B (en) | Geosynchronous satellite orbit uncertain evolution method based on differential algebra | |
| CN110595485A (en) | Low-orbit satellite long-term orbit forecasting method based on two-line root number | |
| Shi et al. | Global search for periodic orbits in the irregular gravity field of a binary asteroid system | |
| Olson et al. | Spin state estimation of tumbling small bodies | |
| Alessi et al. | Desaturation manoeuvres and precise orbit determination for the BepiColombo mission | |
| Teja Nallapu et al. | Attitude control of spacecraft swarms for visual mapping of planetary bodies | |
| Merisio et al. | Characterization of ballistic capture corridors aiming at autonomous ballistic capture at Mars | |
| FREY | Evolution and hazard analysis of orbital fragmentation continua | |
| Dirkx et al. | Propagation and estimation of the dynamical behaviour of gravitationally interacting rigid bodies | |
| CN115242333A (en) | Method for quantitatively characterizing plasma environment of near-earth space satellite group ionized layer | |
| Vittaldev et al. | Unified State Model theory and application in Astrodynamics | |
| CN108959665B (en) | Orbit prediction error empirical model generation method and system suitable for low-orbit satellite | |
| Moore et al. | Solar sails for perturbation relief: Application to asteroids | |
| Feng et al. | Orbit propagation in irregular and uncertain gravity field using differential algebra | |
| Kelly et al. | Geostationary debris mitigation using minimum time solar sail trajectories with eclipse constraints | |
| CN114386282A (en) | Low-orbit giant constellation orbit dynamics analysis method and device of semi-analysis method | |
| Broschart et al. | The small-body dynamics toolkit and associated close-proximity navigation analysis tools at JPL | |
| Zhang et al. | A fast method of estimating the gravitational acceleration of asteroids with irregular shapes | |
| Wright | Orbit Determination Using Vinti's Solution | |
| Cope | Semi-Analytic Approach to Rapid Orbit Determination for the Perturbed Two-Body Problem | |
| Cangahuala | Augmentations to the polyhedral gravity model to facilitate small body navigations | |
| Charls | HERMES CubeSat constellation enhancement by novel pointing strategies and cislunar space extension investigation | |
| Neves | Multi-fidelity modelling of low-energy trajectories for space mission design | |
| Jain et al. | DOCKS Propagator: An Open-source Adaptive Time-step Trajectory Propagator for CubeSat Missions |
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 |